grind.cc

Go to the documentation of this file.

/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2015 - 2023 by Kay F. Jahnke                 */
/*                                                                      */
/*    The git repository for this software is at                        */
/*                                                                      */
/*    https://bitbucket.org/kfj/vspline                                 */
/*                                                                      */
/*    Please direct questions, bug reports, and contributions to        */
/*                                                                      */
/*    kfjahnke+vspline@gmail.com                                        */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/
 
/*! \file grind.cc
 
    \brief performance test. The test is twofold: one aspect is the
    'fidelity' of the filtering operations: if we repeatedly prefilter
    and restore the data, how badly do they suffer, given different data
    types? The second question is, how long does a prefilter-restore
    cycle take (roughly)?
*/
    
#include <vigra/multi_array.hxx>
#include <vigra/accumulator.hxx>
#include <vigra/multi_math.hxx>
#include <iostream>
#include <typeinfo>
#include <random>
#include <ctime>
#include <chrono>
 
#include <vspline/vspline.h>
 
bool verbose = true ; // false ;
 
// 'condense' aggregate types (TinyVectors etc.) into a single value
 
template < typename T >
double condense ( const T & t , std::true_type )
{
  return std::abs ( t ) ;
}
 
template < typename T >
double condense ( const T & t , std::false_type )
{
  return sqrt ( sum ( t * t ) ) / t.size() ;
}
 
template < typename T >
double condense ( const std::complex<T> & t , std::false_type )
{
  return std::abs ( t ) ;
}
 
template < class T >
using is_singular = typename
  std::conditional
  < std::is_fundamental < T > :: value ,
    std::true_type ,
    std::false_type
  > :: type ;
 
template < typename T >
double condense ( const T & t )
{
  return condense ( t , is_singular<T>() ) ;
}
 
template < int dim , typename T >
double check_diff ( vigra::MultiArrayView < dim , T > & reference ,
                    vigra::MultiArrayView < dim , T > & candidate )
{
  using namespace vigra::multi_math ;
  using namespace vigra::acc;
  
  assert ( reference.shape() == candidate.shape() ) ;
  
  vigra::MultiArray < 1 , double >
    error_array ( vigra::Shape1(reference.size() ) ) ;
    
  for ( int i = 0 ; i < reference.size() ; i++ )
  {
    auto help = candidate[i] - reference[i] ;
//     std::cerr << reference[i] << " <> " << candidate[i]
//               << " CFD " << help << std::endl ;
    error_array [ i ] = condense ( help ) ;
  }
 
  AccumulatorChain < double , Select < Mean, Maximum > > ac ;
  extractFeatures ( error_array.begin() , error_array.end() , ac ) ;
  double mean = get<Mean>(ac) ;
  double max = get<Maximum>(ac) ;
  if ( verbose )
  {
    std::cout << "delta Mean:    "
              << mean << std::endl;
    std::cout << "delta Maximum: "
              << max << std::endl;
  }
  return max ;
}
 
/// create an array of random data ditributed between -1 and 1.
/// Then repeatedly prefilter and restore the data, and check the
/// difference between the restored and original data, plus the
/// time needed to do the processing.
/// This test shows that for higher-degree splines, arithmetic
/// precision becomes an issue, and using floats becomes futile.
/// But it also shows that using double data and operating in
/// double precision does not cost too much more processing time.
/// keeping the data in float and using double maths does not
/// help.
 
template < int dim , typename T , typename math_ele_type > void
grind_test ( vigra::TinyVector < int , dim > shape ,
              vspline::bc_code bc ,
              int spline_degree )
{
  typedef vigra::MultiArray < dim , T > array_type ;
  typedef vspline::bspline < T , dim > spline_type ;
  
  vigra::TinyVector < vspline::bc_code , dim > bcv { bc } ;
  spline_type bsp  ( shape , spline_degree , bcv ) ; // , 0.0 ) ;
 
//   if ( verbose )
//     std::cout << "grind: created b-spline:" << std::endl << bsp << std::endl ;
  
  std::random_device rd ;
  std::mt19937 gen ( rd() ) ;
//   gen.seed ( 765 ) ;              // level playing field
  std::uniform_real_distribution<> dis ( -1.0 , 1.0 ) ;
  for ( auto & e : bsp.core )
    e = dis ( gen ) ;         // fill core with random data
  
  array_type reference = bsp.core ; // hold a copy of these data
 
  int times ;
  
  std::chrono::system_clock::time_point start
    = std::chrono::system_clock::now() ;
  
  for ( times = 0 ; times < 1000 ; )
  {
    bsp.prefilter() ;
    vspline::restore < dim , T , math_ele_type > ( bsp ) ;
    times++ ;
    if ( times == 1 || times == 10 || times == 100 || times == 1000 )
    {
      std::cout << "after " << times
                << " prefilter-restore cycles" << std::endl ;
      double emax = check_diff < dim , T > ( reference , bsp.core ) ;  
      if ( verbose )
      {
        // print a summary, can use '| grep CF' on cout
        std::cout
        << typeid(T()).name()
        << " CF "
        << " D " << dim
        << " " << bsp.core.shape()
        << " BC " << vspline::bc_name[bc]
        << " DG " << spline_degree
        << " TOL " << bsp.tolerance
        << " EMAX "
        << emax << std::endl ;
      }
    }
  }
  
  std::chrono::system_clock::time_point end
    = std::chrono::system_clock::now() ;
    
  std::cout << "average cycle duration: "
      << std::chrono::duration_cast<std::chrono::milliseconds>
        ( end - start ) . count() / float ( times )
      << " ms" << std::endl ;
       
  
  std::cout << "---------------------------------------------------" << std::endl ;
  std::cout << std::endl ;
  
}
 
int main ( int argc , char * argv[] )
{
  if ( argc < 2 )
  {
    std::cerr << "pass the spline degree as parameter" << std::endl ;
    exit ( -1 ) ;
  }
  
  int degree = std::atoi ( argv[1] ) ;
  
  std::cout << "performing arithmetics in single precision" << std::endl ;
  grind_test < 2 , float , float > ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
  grind_test < 2 , double , float > ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
  grind_test < 2 , long double , float > ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
  std::cout << "performing arithmetics in double precision" << std::endl ;
  grind_test < 2 , float , double > ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
  grind_test < 2 , double , double > ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
  grind_test < 2 , long double , double > ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
  std::cout << "performing arithmetics in long double precision" << std::endl ;
  grind_test < 2 , float , long double > ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
  grind_test < 2 , double , long double > ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
  grind_test < 2 , long double , long double > ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
}