/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2017 - 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 self_test.cc
///
/// \brief test consistency of precomputed data, prefiltering and evaluation
///
/// the self-test uses three entities: a unit pulse, the reconstruction
/// kernel (which is a unit-spaced sampling of the b-spline basis function's
/// values at integer arguments) and the prefilter. These data have a few
/// fundamental relations we can test:
/// - prefiltering the reconstruction kernel results in a unit pulse
/// - producing a unit-spaced sampling of a spline with only one single
///   unit-valued coefficient produces the reconstruction kernel
///
/// Performing the tests also assures us that the evaluation machinery with
/// it's 'weight matrix' does what's intended, and that access to the basis
/// function and it's derivatives (see basis_functor) functions correctly.
///   
/// With rising spline degree, the test is ever more demanding. this is
/// reflected in the maximum error returned for every degree: it rises
/// with the spline degree. With the complex operations involved, seeing
/// a maximal error in the order of magnitude of 1e-12 for working with
/// long doubles seems reasonable enough (on my system).
///
/// I assume that the loss of precision with high degrees is mainly due to
/// the filter's overall gain becoming very large. Since the gain is 
/// applied as a factor before or during prefiltering and prefiltering
/// has the reverse effect, in the sum we end up having the effect of first
/// multiplying with, then dividing by a very large number, 'crushing'
/// the values to less precision. In bootstrap.cc, I perfom the test
/// with GMP high precision floats and there I can avoid the problem, since
/// the magnitude of the numbers I use there is well beyond the magnitude
/// of the gain occuring with high spline degrees. So the conclusion is that
/// high spline degrees can be used, but may not produce very precise results,
/// since the normal C++ types are hard pressed to cope with the dynamic
/// range covered by the filter.
///
/// The most time is spent calculating the basis function values recursively
/// using cdb_bspline_basis, for cross-reference.
///
/// compile: clang++ -O3 -std=c++11 self_test.cc -o self_test -pthread

#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>
#include <numeric>
#include <assert.h>
#include <random>
#include <vspline/vspline.h>

long double circular_test_previous ;

template < typename dtype >
dtype self_test ( int degree , dtype threshold , dtype strict_threshold )
{
  if ( degree == 0 )
    circular_test_previous = 1 ;
  
  dtype max_error = 0 ;
  
  // self-test for plausibility. we know that using the b-spline
  // prefilter of given degree on a set of unit-spaced samples of
  // the basis function (at 0, +/-k) should yield a unit pulse.
  
  // we create a bspline object for 100 core coefficients
  
  typedef vspline::bspline < dtype , 1 > spline_type ;
  spline_type bspl ( 100 , degree ) ;
            
  // next we create an evaluator for this spline.
  // Note how, to be as precise as possible for this test,
  // we specify 'double' as elementary type for coordinates.
  // This is overkill, but so what...
  
  auto ev = vspline::make_evaluator < spline_type , double >
                    ( bspl ) ;
  
  // and two arrays with the same size as the spline's 'core'.
            
  vigra::MultiArray<1,dtype> result ( 100 ) ;
  vigra::MultiArray<1,dtype> reference ( 100 ) ;
  
  // we obtain a pointer to the reference array's center
  
  dtype * p_center = &(reference[50]) ;

  // we can obtain the reconstruction kernel by accessing precomputed
  // basis function values via bspline_basis_2()
  
  for ( int x = - degree / 2 ; x <= degree / 2 ; x++ )
  {
    p_center[x] = vspline::bspline_basis_2<dtype> ( x+x , degree ) ;
  }
  
  // alternatively we can put a unit pulse into the center of the
  // coefficient array, transform and assign back. Transforming
  // the unit pulse produces the reconstruction kernel. Doing so
  // additionally assures us that the evaluation machinery with
  // it's 'weight matrix' is functioning correctly.
  
  // we obtain a pointer to the coefficient array's center
  
  p_center = &(bspl.core[50]) ;
  
  *p_center = 1 ;
  
  vspline::transform ( ev , result ) ;

  bspl.core = result ;
  
  // we compare the two versions of the reconstruction kernel we
  // have to make sure they agree. the data should be identical.
  // we also compare with the result of a vspline::basis_functor,
  // which uses the same method of evaluating a b-spline with a
  // single unit-valued coefficient.
  // Here we expect complete agreement.
  
  vspline::basis_functor<dtype> bf ( degree ) ;
  
  for ( int x = 50 - degree / 2 ; x <= 50 + degree / 2 ; x++ )
  {
    assert ( result[x] == reference[x] ) ;
    assert ( bf ( x - 50 ) == reference[x] ) ;
  }
  
  // now we apply the prefilter, expecting that afterwards we have
  // a single unit pulse coinciding with the location where we
  // put the center of the kernel. This test will exceed the strict
  // threshold, but the ordinary one will hold.
  
  bspl.prefilter() ;
  
  // we test our predicate
  
  for ( int x = - degree / 2 ; x <= degree / 2 ; x++ )
  {
    dtype error ;
    if ( x == 0 )
    {
      // at the origin we expect a value of 1.0
      error = std::fabs ( p_center [ x ] - 1.0 ) ;
    }
    else
    {
      // off the origin we expect a value of 0.0
      error = std::fabs ( p_center [ x ] ) ;
    }
    if ( error > threshold )
      std::cout << "unit pulse test, x = " << x << ", error = "
                << error << std::endl ;
    max_error = std::max ( max_error , error ) ;
  }
  
  // test bspline_basis() at k/2, k E N against precomputed values
  // while bspline_basis at whole arguments has delta == 0 and hence
  // makes no use of rows beyond row 0 of the weight matrix, arguments
  // at X.5 use all these rows. We can test against bspline_basis_2,
  // which provides precomputed values for half unit steps.
  // we run this test with strict_threshold.
  
  int xmin = - degree - 1 ;
  int xmax = degree + 1 ;
  
  for ( int x2 = xmin ; x2 <= xmax ; x2++ )
  {
    auto a = bf ( x2 / 2.0L ) ;
    auto b = vspline::bspline_basis_2<dtype> ( x2 , degree ) ;
    auto error = std::abs ( a - b ) ;
    if ( error > strict_threshold )
      std::cout << "bfx2: " << x2 / 2.0 << " : "
                << a << " <--> " << b
                << " error " << error << std::endl << std::endl ;
    max_error = std::max ( max_error , error ) ;
  }
  
  // set all coefficients to 1, evaluate. result should be 1,
  // because every set of weights is a partition of unity
  // this is a nice test, because it requires participation of
  // all rows in the weight matrix, since the random arguments
  // produce arbitrary delta. we run this test with strict_threshold.
  {
    std::random_device rd ;
    std::mt19937 gen ( rd() ) ;
    std::uniform_real_distribution<>
       dis ( 50 - degree -1 , 50 + degree + 1 ) ;
       
    bspl.container = 1 ;
    for ( int k = 0 ; k < 1000 ; k++ )
    {
      double x = dis ( gen ) ;
      dtype y = ev ( x ) ;
      dtype error = std::abs ( y - 1 ) ;
      if ( error > strict_threshold )
        std::cout << "partition of unity test, d0: " << x << " : "
                  << y << " <--> " << 1
                  << " error " << error << std::endl << std::endl ;
      max_error = std::max ( max_error , error ) ;
    }
    
    vigra::TinyVector < int , 1 > deriv_spec ;
    
    // we also test evaluation of derivatives up to 2.
    // Here, with the coefficients all equal, we expect 0 as a result.
    
    if ( degree > 1 )
    {
      deriv_spec[0] = 1 ;
      auto dev = vspline::make_evaluator < spline_type , double >
                          ( bspl , deriv_spec ) ;
      for ( int k = 0 ; k < 1000 ; k++ )
      {
        double x = dis ( gen ) ;
        dtype y = dev ( x ) ;
        dtype error = std::abs ( y ) ;
        if ( error > strict_threshold )
          std::cout << "partition of unity test, d1: " << x << " : "
                    << y << " <--> " << 0
                    << " error " << error << std::endl << std::endl ;
        max_error = std::max ( max_error , error ) ;
      }
    }
    
    if ( degree > 2 )
    {
      deriv_spec[0] = 2 ;
      auto ddev = vspline::make_evaluator< spline_type , double >
                           ( bspl , deriv_spec ) ;
      for ( int k = 0 ; k < 1000 ; k++ )
      {
        double x = dis ( gen ) ;
        dtype y = ddev ( x ) ;
        dtype error = std::abs ( y ) ;
        if ( error > strict_threshold )
          std::cout << "partition of unity test, d2: " << x << " : "
                    << y << " <--> " << 0
                    << " error " << error << std::endl << std::endl ;
        max_error = std::max ( max_error , error ) ;
      }
    }
  }
  
  // for smaller degrees, the cdb recursion is usable, so we can
  // doublecheck basis_functor for a few sample x. The results here
  // are also very precise, so we use strict_threshold. Initially
  // I took this up to degree 19, but now it's only up to 13, which
  // should be convincing enough and is much faster.
  
  if ( degree < 13 ) // was: 19
  {
    std::random_device rd ;
    std::mt19937 gen ( rd() ) ;
    std::uniform_real_distribution<> dis ( - degree -1 , degree + 1 ) ;
    for ( int k = 0 ; k < 1000 ; k++ )
    {
      dtype x = dis ( gen ) ;
      dtype a = bf ( x ) ;
      dtype b = vspline::cdb_bspline_basis ( x , degree ) ;
      dtype error = std::abs ( a - b ) ;
      if ( error > strict_threshold )
        std::cout << "bf: " << x << " : "
                  << a << " <--> " << b
                  << " error " << error << std::endl << std::endl ;
      max_error = std::max ( max_error , error ) ;
    }
  }
  
  if ( degree > 1 && degree < 13 )
  {
    vspline::basis_functor<dtype> dbf ( degree , 1 ) ;
    std::random_device rd ;
    std::mt19937 gen ( rd() ) ;
    std::uniform_real_distribution<> dis ( - degree -1 , degree + 1 ) ;
    for ( int k = 0 ; k < 1000 ; k++ )
    {
      dtype x = dis ( gen ) ;
      dtype a = dbf ( x ) ;
      dtype b = vspline::cdb_bspline_basis ( x , degree , 1 ) ;
      dtype error = std::abs ( a - b ) ;
      if ( error > strict_threshold )
        std::cout << "dbf: " << x << " : "
                  << a << " <--> " << b
                  << " error " << error << std::endl << std::endl ;
      max_error = std::max ( max_error , error ) ;
    }
  }
  
  if ( degree > 2 && degree < 13 )
  {
    vspline::basis_functor<dtype> ddbf ( degree , 2 ) ;
    std::random_device rd ;
    std::mt19937 gen ( rd() ) ;
    std::uniform_real_distribution<> dis ( - degree -1 , degree + 1 ) ;
    for ( int k = 0 ; k < 1000 ; k++ )
    {
      dtype x = dis ( gen ) ;
      dtype a = ddbf ( x ) ;
      dtype b = vspline::cdb_bspline_basis ( x , degree , 2 ) ;
      dtype error = std::abs ( a - b ) ;
      if ( error > strict_threshold )
        std::cout << "ddbf: " << x << " : "
                  << a << " <--> " << b
                  << " error " << error << std::endl << std::endl ;
      max_error = std::max ( max_error , error ) ;
   }
  }
  
  if ( degree > 0 )
  {
    std::random_device rd ;
    std::mt19937 gen ( rd() ) ;
    std::uniform_real_distribution<>
       dis ( -1 , 1 ) ;
    
    dtype circle_error = 0 ;
    dtype avg_circle_error = 0 ;
    
    // consider a spline with a single 1.0 coefficient at the origin

    // reference is the spline's value at ( 1 , 0 ), which is
    // certainly on the unit circle
    
    dtype v2 = bf ( 1 ) * bf ( 0 ) ;
    
    // let's assume 10000 evaluations is a good enough
    // statistical base

    for ( int k = 0 ; k < 10000 ; k++ )
    {
      // take a random x and y coordinate
      
      double x = dis ( gen ) ;
      double y = dis ( gen ) ;
      
      // normalize to unit circle
      
      double s = sqrt ( x * x + y * y ) ;
      
      x /= s ;
      y /= s ;
      
      // and take the value of the spline there, which is
      // the product of the basis function values
      
      dtype v1 = bf ( x ) * bf ( y ) ;
    
      // we assume that with rising spline degree, the difference
      // between these two values will become ever smaller, as the
      // equipotential lines of the splines become more and
      // more circular
      
      dtype error = std::abs ( v1 - v2 ) ;
      circle_error = std::max ( circle_error , error ) ;
      avg_circle_error += error ;
    }
    
    assert ( circle_error < circular_test_previous ) ;
    circular_test_previous = circle_error ;
    
    // in my tests, circle_error goes down to ca 7.4e-7,
    // so with degree 45 evaluations on the unit circle
    // differ very little from each other.
  
//     std::cout << "unit circle test, degree " << degree
//               << " emax = " << circle_error
//               << " avg(e) = " << avg_circle_error / 10000
//               << " value: " << v2 << std::endl ;
  }
//   std::cout << "max error for degree " << degree
//             << ": " << max_error << std::endl ;
//             
  return max_error ;
}

// run a self test of vspline's constants, prefiltering and evaluation.
// This tests if a set of common operations produces larger errors than
// anticipated, to alert us if something has gone amiss.
// The thresholds are fixed heuristically to be quite close to the actually
// occuring maximum error.

int main ( int argc , char * argv[] )
{
  long double max_error_l = 0 ;
  
  for ( int degree = 0 ; degree <= vspline_constants::max_degree ; degree++ )
  {
    max_error_l = std::max ( max_error_l ,
                             self_test<long double>
                             ( degree , 4e-13l , 1e-18 ) ) ;
  }
  
  std::cout << "maximal error of tests with long double precision: "
  
  << max_error_l << std::endl ;
  
  double max_error_d = 0 ;
  
  for ( int degree = 0 ; degree <= vspline_constants::max_degree ; degree++ )
  {
    max_error_d = std::max ( max_error_d ,
                             self_test<double>
                             ( degree , 1e-9 , 7e-16 ) ) ;
  }
  
  std::cout << "maximal error of tests with double precision: "
  << max_error_d << std::endl ;
  
  float max_error_f = 0 ;
  
  // test float only up to degree 15.

  for ( int degree = 0 ; degree < 16 ; degree++ )
  {
    max_error_f = std::max ( max_error_f ,
                             self_test<float>
                             ( degree , 3e-6 , 4e-7 ) ) ;
  }
  
  std::cout << "maximal error of tests with float precision: "
  << max_error_f << std::endl ;
  
  if (    max_error_l < 4e-13
       && max_error_d < 1e9
       && max_error_f < 3e-6 )
    
    std::cout << "test passed" << std::endl ;
  else
    std::cout << "test failed" << std::endl ;
}