restore_test.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                    */
/*                                                                      */
/*    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.                                   */
/*                                                                      */
/************************************************************************/
 
/// restore_test - create b-splines from random data and restore the original
/// data from the spline. This has grown to function as a unit test
/// instantiating all sorts of splines and evaluators and using the evaluators
/// with the functions in transform.h.
///
/// Ideally, this test should result in restored data which are identical to
/// the initial random data, but several factors come into play:
///
/// the precision of the coefficient calculation is parametrizable in vspline,
/// via the 'tolerance' parameter, or, in some places, the 'horizon' parameter.
///
/// the data type of the spline is important - single precision data are unfit
/// to produce coefficients capable of reproducing the data very precisely
///
/// the spline degree is important. High degrees need wide horizons, but wide
/// horizons also need many calculation steps, which may introduce errors.
///
/// The dimension of the spline is also important. Since higher-dimension
/// splines need more calculations for prefiltering and evaluation, results are
/// less precise. This test program can test up to 4D.
///
/// Please note that, since this program uses a good many different data types
/// and tries to test as much of vspline's code as possible, compile times may
/// be very long, especially if the higher-D tests are not commented out.
///
/// With dimensions > 2 running this program in full takes a long time, since
/// it tries to be comprehensive to catch any corner cases which might have
/// escaped scrutiny. Run time of this program is not a very good indicator
/// of vspline's speed, since many operations (like statistics, array copy
/// operations) are neither multithreaded nor vectorized. For speed tests,
/// 'roundtrip' or 'grind' are better.
///
/// A quick way to test is to run only the 1D tests and only look at the
/// maximum error for each data type, like:
///
/// ./hwy_restore_test_clang++ 1 | grep 'max err'
///
/// This should create output like:
///
/// test for type FfvE: max error = +0.000001370906829834
/// test for type FdvE: max error = +0.000000000000002887
/// test for type FevE: max error = +0.000000000000000001
/// test for type FSt7complexIfEvE: max error = +0.000001788139343262
/// test for type FSt7complexIdEvE: max error = +0.000000000000003442
/// test for type FSt7complexIeEvE: max error = +0.000000000000000002
/// test for type FN5vigra8RGBValueIfLj0ELj1ELj2EEEvE: max error = +0.000000963557215769
/// test for type FN5vigra8RGBValueIdLj0ELj1ELj2EEEvE: max error = +0.000000000000001859
/// test for type FN5vigra8RGBValueIeLj0ELj1ELj2EEEvE: max error = +0.000000000000000001
/// test for type FN5vigra10TinyVectorIfLi3EEEvE: max error = +0.000001170033762005
/// test for type FN5vigra10TinyVectorIdLi3EEEvE: max error = +0.000000000000001602
/// test for type FN5vigra10TinyVectorIeLi3EEEvE: max error = +0.000000000000000001
/// test for type FN5vigra10TinyVectorIfLi2EEEvE: max error = +0.000000927230667240
/// test for type FN5vigra10TinyVectorIdLi2EEEvE: max error = +0.000000000000002355
/// test for type FN5vigra10TinyVectorIfLi1EEEvE: max error = +0.000001609325408936
/// test for type FN5vigra10TinyVectorIdLi1EEEvE: max error = +0.000000000000002776
 
#include <vigra/multi_array.hxx>
#include <vigra/accumulator.hxx>
#include <vigra/multi_math.hxx>
#include <iostream>
#include <typeinfo>
#include <random>
 
#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>() ) ;
}
 
// compare two arrays and calculate the mean and maximum difference
 
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 ;
}
 
/// do a restore test. This test fills the array that's
/// passed in with small random data, constructs a b-spline with the requested
/// parameters over the data, then calls vspline::restore(), which evaluates the
/// spline at discrete locations.
/// While this test fails to address several aspects (derivatives, behaviour
/// at locations which aren't discrete), it does make sure that prefiltering
/// has produced a correct result and reconstruction succeeds: If the spline
/// coefficients were wrong, reconstruction would fail just as it would if
/// the coefficients were right and the evaluation code was wrong. That both
/// should be wrong and accidentally produce a correct result is highly unlikely.
///
/// This routine has grown to be more of a unit test for all of vspline,
/// additional vspline functions are executed and the results are inspected.
/// This way we can assure that the transform-type routines are usable with
/// all supported data types and vectorized and unvectorized results are
/// consistent.
 
template < int dim , typename T >
double restore_test ( vigra::MultiArrayView < dim , T > & arr ,
                      vspline::bc_code bc ,
                      int spline_degree )
{
  if ( verbose )
    std::cout << "**************************************" << std::endl
              << "testing type " << typeid(T()).name() << std::endl ;
              
  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  ( arr.shape() , spline_degree , bcv , 0.0 ) ;
 
  if ( verbose )
    std::cout << "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 : arr )
    e = dis ( gen ) ;
  array_type reference = arr ;
  
  bsp.prefilter ( arr ) ;
  
  vspline::restore < dim , T > ( bsp , arr ) ;
 
  if ( verbose )
    std::cout << "after restoration of original data:" << std::endl ;
  
  double emax = check_diff < dim , T > ( reference , arr ) ;
  
  if ( verbose )
  {
    // print a summary, can use '| grep CF' on cout
    std::cout
    << typeid(T()).name()
    << " CF "
    << " D " << dim
    << " " << arr.shape()
    << " BC " << vspline::bc_name[bc]
    << " DG " << spline_degree
    << " TOL " << bsp.tolerance
    << " EMAX "
    << emax << std::endl ;
  }
  
//   // have: coeffs in bsp, original data in reference, restored data in arr
//
// the next bit of code implements a trial of 'spline polishing':
// the coefficients are filtered with the reconstruction kernel, and
// the original data are subtracted from resulting 'restored' data,
// producing an array of errors. A spline is erected over the
// errors and it's coefficents are subtracted from the coefficients
// in the first spline, 'augmenting' it to better represent the
// original data. Seems to squash the error pretty much to arithmetic
// precision.
  
  // difference original / restored
  arr -= reference ;
  // create another bspline
  spline_type polish ( arr.shape() , spline_degree , bcv , 0.0 ) ;
  // prefilter it, sucking in the difference
  polish.prefilter ( arr ) ;
  // and layer the resulting coeffs on top of the previous set
  bsp.core -= polish.core ;
  // brace the 'augmented' spline
  bsp.brace() ;
  // and check it's quality
  vspline::restore < dim , T > ( bsp , arr ) ;
 
  if ( verbose )
    std::cout << "after polishing of spline:" << std::endl ;
 
  double emax2 = check_diff < dim , T > ( reference , arr ) ;
  
  if ( verbose )
  {
    // print a summary, can use '| grep CF' on cout
    std::cout
    << typeid(T()).name()
    << " CF "
    << " D " << dim
    << " " << arr.shape()
    << " BC " << vspline::bc_name[bc]
    << " DG " << spline_degree
    << " TOL " << bsp.tolerance
    << " EMAX "
    << emax2 << std::endl ;
  }
  
  // test the factory functions make_evaluator and make_safe_evaluator
 
  auto raw_ev =  vspline::make_evaluator
                      < spline_type , float > ( bsp ) ;
                      
  typedef typename decltype ( raw_ev ) :: in_type fc_type ;
  typedef typename decltype ( raw_ev ) :: in_nd_ele_type nd_fc_type ;
  
  // fc_type is what the evaluator expects as input type. This may be a plain
  // fundamental, so we produce an nD view to it for uniform handling.
  
  fc_type cs ( 0 ) ;
  nd_fc_type & nd_cs ( reinterpret_cast < nd_fc_type & > ( cs ) ) ;
  
  auto rs = raw_ev ( cs ) ;
  raw_ev.eval ( cs , rs ) ;
 
  // try evaluating the spline at it's lower and upper limit. We assign the
  // spline's limit values for each dimension via the nD view to cs:
 
  for ( int d = 0 ; d < dim ; d++ )
    nd_cs[d] = bsp.lower_limit ( d ) ;
  
  raw_ev.eval ( cs , rs ) ;
  
  if ( verbose )
    std::cout << cs << " -> " << rs << std::endl ;
  
  for ( int d = 0 ; d < dim ; d++ )
    nd_cs[d] = bsp.upper_limit ( d ) ;
 
  raw_ev.eval ( cs , rs ) ;
  
  if ( verbose )
    std::cout << cs << " -> " << rs << std::endl ;
  
  // additionally, we perform a test with a 'safe evaluator' and random
  // coordinates. For a change we use double precision coordinates
  
  auto _ev = vspline::make_safe_evaluator
                      < spline_type , double > ( bsp ) ;
 
  enum { vsize = decltype ( _ev ) :: vsize } ;
  
  typedef typename decltype ( _ev ) :: in_type coordinate_type ;
  typedef typename decltype ( _ev ) :: in_ele_type rc_type ;
  typedef typename decltype ( _ev ) :: out_ele_type ele_type ;
 
  // throw a domain in to test that as well:
 
  auto dom = vspline::domain
             < coordinate_type , spline_type , vsize >
               ( coordinate_type(-.3377) ,
                 coordinate_type(3.11) ,
                 bsp ) ;
 
  auto ev = vspline::grok ( dom + _ev ) ;
  
  coordinate_type c ;
  
  vigra::MultiArray < 1 , coordinate_type > ca ( vigra::Shape1 ( 10003 ) ) ;
  vigra::MultiArray < 1 , T > ra ( vigra::Shape1 ( 10003 ) ) ;
  vigra::MultiArray < 1 , T > ra2 ( vigra::Shape1 ( 10003 ) ) ;
  
  auto pc = (rc_type*) &c ;
  // make sure we can evaluate at the lower and upper limit
  int k = 0 ;
  {
    for ( int e = 0 ; e < dim ; e++ )
      pc[e] = 0.0 ;
    ra[k] = ev ( c ) ;
    ca[k] = c ;
  }
  k = 1 ;
  {
    for ( int e = 0 ; e < dim ; e++ )
      pc[e] = 1.0 ;
    ra[k] = ev ( c ) ;
    ca[k] = c ;
  }
  k = 2 ;
  {
    for ( int e = 0 ; e < dim ; e++ )
      pc[e] = -2.0 ;
    ra[k] = ev ( c ) ;
    ca[k] = c ;
  }
  k = 3 ;
  {
    for ( int e = 2.0 ; e < dim ; e++ )
      pc[e] = 1.0 ;
    ra[k] = ev ( c ) ;
    ca[k] = c ;
  }
  
  // the remaining coordinates are picked randomly
  
  std::uniform_real_distribution<> dis2 ( -2.371 , 2.1113 ) ;
  for ( k = 4 ; k < 10003 ; k++ )
  {
    for ( int e = 0 ; e < dim ; e++ )
      pc[e] = dis2 ( gen ) ;
    ra[k] = ev ( c ) ;
    ca[k] = c ;
  }
  
  // run an index-based transform. With a domain in action, this
  // does *not* recreate the original data
 
  vspline::transform ( ev , arr ) ;
 
  // run an array-based transform. result should be identical
  // to single-value-eval above, within arithmetic precision limits
  // given the specific optimization level.
  
  // we can produce different views on the data to make sure feeding
  // the coordinates still works correctly:
  
  vigra::MultiArrayView < 2 , coordinate_type > ca2d
    ( vigra::Shape2 ( 99 , 97 ) , ca.data() ) ;
  vigra::MultiArrayView < 2 , T > ra2d
    ( vigra::Shape2 ( 99 , 97 ) , ra2.data() ) ;
    
  vigra::MultiArrayView < 2 , coordinate_type > ca2ds
    ( vigra::Shape2 ( 49 , 49 ) , vigra::Shape2 ( 2 , 200 ) , ca.data() ) ;
  vigra::MultiArrayView < 2 , T > ra2ds
    ( vigra::Shape2 ( 49 , 49 ) , vigra::Shape2 ( 2 , 200 ) , ra2.data() ) ;
    
  vigra::MultiArrayView < 3 , coordinate_type > ca3d
    ( vigra::Shape3 ( 20 , 21 , 19 ) , ca.data() ) ;
  vigra::MultiArrayView < 3 , T > ra3d
    ( vigra::Shape3 ( 20 , 21 , 19 ) , ra2.data() ) ;
  
  vspline::transform ( ev , ca , ra2 ) ;
  vspline::transform ( ev , ca2d , ra2d ) ;
  vspline::transform ( ev , ca2ds , ra2ds ) ;
  vspline::transform ( ev , ca3d , ra3d ) ;
 
  // vectorized and unvectorized operation may produce slightly
  // different results in optimized code.
  // usually assert ( ra == ra2 ) holds, but we don't consider
  // it an error if it doesn't and allow for a small difference
  
  auto dsv = check_diff < 1 , T > ( ra , ra2 ) ;
  auto tolerance = 10 * std::numeric_limits<ele_type>::epsilon() ;
  
  if ( dsv > tolerance )
    std::cout << vspline::bc_name[bc]
              << " - max difference single/vectorized eval "
              << dsv << std::endl ;
              
  if ( dsv > .0001 )
  {
    std::cout << vspline::bc_name[bc]
              << " - excessive difference single/vectorized eval "
              << dsv << std::endl ;
    for ( int k = 0 ; k < 10003 ; k++ )
    {
      if ( ra[k] != ra2[k] )
      {
        std::cout << "excessive at k = " << k
                  << ": " << ca[k] << " -> "
                  << ra[k] << ", " << ra2[k] << std::endl ;
      }
    }
  }
  
  return emax ;
}
 
using namespace vspline ;
 
template < class arr_t >
double view_test ( arr_t & arr )
{
  double emax = 0.0 ;
  
  enum { dimension = arr_t::actual_dimension } ;
  typedef typename arr_t::value_type value_type ;
  vspline::bc_code bc_seq[] { PERIODIC , MIRROR , REFLECT , NATURAL } ;
 
  // TODO: for degree-1 splines, I sometimes get different results
  // for unvectorized and vectorized operation. why?
  
  for ( int spline_degree = 0 ; spline_degree < 8 ; spline_degree++ )
  {
    for ( auto bc : bc_seq )
    {
      auto e = restore_test < dimension , value_type >
                ( arr , bc , spline_degree ) ;
      if ( e > emax )
        emax = e ;
    }
  }
  return emax ;
}
 
int d0[] { 1 , 2 , 3 , 5 , 8 , 13 , 16 , 21 ,
           34 , 55 , 123 , 128 , 289 , 500 , 1031 ,
           2001 , 4999 } ;
int d1[] { 1 , 2 , 3 , 5 , 8 , 13 , 16 , 21 ,
           34 , 89 , 160 , 713 } ;
int d2[] { 2 , 3 , 5 , 8 , 13 , 21 } ;
int d3[] { 2 , 3 , 5 , 8 , 13 } ;
 
int* dext[] { d0 , d1 , d2 , d3 } ;
int dsz[] { sizeof(d0) / sizeof(int) ,
            sizeof(d1) / sizeof(int) ,
            sizeof(d2) / sizeof(int) ,
            sizeof(d3) / sizeof(int) } ;
 
template < int dim , typename T >
struct test
{
  typedef vigra::TinyVector < int , dim > shape_type ;
  double emax = 0.0 ;
  
  double operator() ()
  {
    shape_type dshape ; 
    for ( int d = 0 ; d < dim ; d++ )
      dshape[d] = dsz[d] ;
    
    vigra::MultiCoordinateIterator<dim> i ( dshape ) ,
    end = i.getEndIterator() ;  
    
    while ( i != end )
    {
      shape_type shape ;
      for ( int d = 0 ; d < dim ; d++ )
        shape[d] = * ( dext[d] + (*i)[d] ) ;
      
      vigra::MultiArray < dim , T > _arr ( 2 * shape + 1 ) ;
      auto stride = _arr.stride() * 2 ;
      vigra::MultiArrayView < dim , T >
        arr ( shape , stride , _arr.data() + long ( sum ( stride ) ) ) ;
      auto e = view_test ( arr ) ;
      if ( e > emax )
        emax = e ;
      ++i ;
      
      // make sure that we have only written back to 'arr', leaving
      // _arr untouched
      for ( auto & e : arr )
        e = T(0.0) ;
      for ( auto e : _arr )
        assert ( e == T(0.0) ) ;
    }
    return emax ;
  }
} ;
 
template < typename T >
struct test < 0 , T >
{ 
  double operator() ()
  { 
    return 0.0 ;
  } ;
} ;
 
template < int dim ,
           typename tuple_type ,
           int ntypes = std::tuple_size<tuple_type>::value >
struct multitest
{
  void operator() ()
  {
    typedef typename std::tuple_element<ntypes-1,tuple_type>::type T ;
    auto e = test < dim , T >() () ;
    std::cout << "test for type " << typeid(T()).name()
              << ": max error = " << e << std::endl ;
    multitest < dim , tuple_type , ntypes - 1 >() () ;
  }
} ;
 
template < int dim ,
           typename tuple_type >
struct multitest < dim , tuple_type , 0 >
{
  void operator() ()
  {
  }
} ;
 
int main ( int argc , char * argv[] )
{
  std::cout << std::fixed << std::showpos << std::showpoint
            << std::setprecision(18);
  std::cerr << std::fixed << std::showpos << std::showpoint
            << std::setprecision(18);
 
  int test_dim = 2 ;
  if ( argc > 1 )
    test_dim = std::atoi ( argv[1] ) ;
  if ( test_dim > 4 )
    test_dim = 4 ;
  
  std::cout << "testing with " << test_dim << " dimensions" << std::endl ;
  
  typedef std::tuple <
                       vigra::TinyVector < double , 1 > ,
                       vigra::TinyVector < float , 1 > ,
                       vigra::TinyVector < double , 2 > ,
                       vigra::TinyVector < float , 2 > ,
                       vigra::TinyVector < long double , 3 > ,
                       vigra::TinyVector < double , 3 > ,
                       vigra::TinyVector < float , 3 >  ,
                       vigra::RGBValue<long double> ,
                       vigra::RGBValue<double> ,
                       vigra::RGBValue<float,0,1,2> ,
                       std::complex < long double > ,
                       std::complex < double > ,
                       std::complex < float > ,
                       long double ,
                       double ,
                       float
 > tuple_type ;
 
 switch ( test_dim )
 {
   case 1:
    multitest < 1 , tuple_type >() () ;
    break ;
   case 2:
    multitest < 2 , tuple_type >() () ;
    break ;
//    case 3:
//     multitest < 3 , tuple_type >() () ;
//     break ;
//    case 4:
//     multitest < 4 , tuple_type >() () ;
//     break ;
   default:
     break ;
 }
 std::cout << "terminating" << std::endl ;
}