/************************************************************************/
/*                                                                      */
/*    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.                                   */
/*                                                                      */
/************************************************************************/

/*! \file scope_test.cc

    \brief tries to fathom vspline's scope
  
    vspline is very flexible when it comes to types and other
    compile-time issues. It's not really feasible to instantiate
    every single possible combination of compile-time parameters,
    but this program is an attempt to test a few relevant ones,
    specifically for vspline's 'higher' API, namely classes
    bspline and evaluator, and the functions in transform.h.
    
    We test data in 1, 2 and 3 dimensions. For individual values,
    we use plain fundamentals and TinyVectors of two and three
    fundamentals, with the fundamental types float, double and
    long double. Wherever possible, we perform the tests
    with several vectorization widths, to make sure that the
    wielding code performs as expected.
    
    For every such combination, we perform this test sequence:
    - produce random data and random coordinates
    - create a b-spline with given degree and boundary conditions
      from the data
    - create a 'safe evaluator' for the spline
    - use the evaluator to bulk-evaluate the random coordinates
    - and to evaluate them one-by-one; compare the result
    - restore the original data from the coefficients using
      an index-based transform and also using vspline::restore;
      compare the results.
      
    On my system, this is just about what the compiler can handle,
    and the compilation takes a long time. The test can be limited
    to smaller parameter sets, a good place to narrow the scope is
    by picking less rc_type/ele_type combinations in doit().
    
    This test passes if
    - the code compiles at all
    - the program runs until completion
    - the errors are near the resolution you'd expect from
      the types involved, so you'd expect output like
    
    ...
    vsz +5 chn +2 trg_dim +3 cf_dim +3
    float crd, float value 
    d Mean: +0.000000763713326440d Max: +0.000007158763497749
    d Mean: +0.000000757349428671d Max: +0.000007168434323489
    double crd, double value 
    d Mean: +0.000000000000009688d Max: +0.000000000000899253
    d Mean: +0.000000000000009681d Max: +0.000000000000899253
    long double crd, long double value 
    d Mean: +0.000000000000000001d Max: +0.000000000000000006
    d Mean: +0.000000000000000001d Max: +0.000000000000000006
    ...
    
    This test does not focus on different spline degrees or
    boundary conditions, it's purpose is to test compile-time
    parametrization. But it can easily be modified to use
    different spline degrees and boundary conditions by passing
    the relevant parameters to test(); I have limited the
    choice to a few 'typical' combinations. Another aspect which
    is not tested is variation of the extents of the arrays of
    data involved. Other tests, like 'restore_test.cc' and
    'roundtrip.cc' focus more on these run-time parameters.
    
    Testing varying vectorization widths is for completeness'
    sake, normally using the defaults should give good performance.
*/

#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
// here we also take the 'value_extent' - the largest coefficient
// value we have used - to put the error in relation, so that the high
// coefficient values we use for integral coefficients don't produce
// results which look wronger than they are ;)

template < unsigned int dim , typename T >
double check_diff ( vigra::MultiArrayView < dim , T > & reference ,
                    vigra::MultiArrayView < dim , T > & candidate ,
                    double value_extent
                  )
{
  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] ;
    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 << "rel. error Mean: " << mean / value_extent
              << " rel. error Max: "  << max / value_extent << std::endl;
  }
  return max ;
}

using namespace vspline ;

#define ast(x) std::integral_constant<int,x>

// with this test routine, the idea is to test all of vspline's
// higher functions with all possible compile time parameters.

template < unsigned int cf_dim ,  // dimension of spline/coefficient array
           unsigned int trg_dim , // dimension of target (result) array
           int chn ,              // number of channels in spline's data type
           int vsz ,              // vectorization width 
           typename rc_type ,     // elementary type of a real coordinate
           typename ele_type ,    // elementary type of coefficients
           typename math_ele_type > // elementary type used for arithmetics
void test ( int spline_degree = 3 ,
            vspline::bc_code bc = vspline::MIRROR )
{
  // TODO: when operating on integral values, vectorized and unvectorized
  // operation sometimes produces results differing by at most 1, hence
  // this expression for 'tolerance'. I'm not sure why this happens, it
  // looks like different results of rounding/truncation, should track
  // this down make sure the logic isn't flawed somewhere
  
  double tolerance = std::is_integral < ele_type > :: value
                     ? 1.0 : 0.000001 ;
                     
  // for integral coefficients, we use high knot point values to
  // provide sufficient dynamic range. the exact values are chosen
  // in a slightly ad-hoc manner, but they should demonstrate that
  // more bits can provide better results.

  double value_extent = std::is_integral < ele_type > :: value
                        ? sizeof ( ele_type ) == 2
                          ? 1000.0
                          : sizeof ( ele_type ) == 4
                            ? 1000000.0
                            : sizeof ( ele_type ) == 8
                              ? 1000000000000.0
                              : 100.0
                        : 1.0 ;
                         
  typedef typename vigra::MultiArrayShape<cf_dim>::type cf_shape_t ;
  typedef typename vigra::MultiArrayShape<trg_dim>::type trg_shape_t ;
  
  typedef typename std::conditional
                   < chn == 1 ,
                     ele_type ,
                     vigra::TinyVector < ele_type , chn >
                   > :: type dtype ;
  
  typedef typename std::conditional
                   < cf_dim == 1 ,
                     rc_type ,
                     vigra::TinyVector < rc_type , cf_dim >
                   > :: type crd_type ;
  
  // allocate storage for original data
                   
  cf_shape_t cf_shape { 99 } ;
  
  vigra::MultiArray < cf_dim , dtype >
    original ( cf_shape ) ;
    
  vigra::MultiArray < cf_dim , dtype >
    restored ( cf_shape ) ;

  // and storage for coordinates and results
    
  trg_shape_t trg_shape { 101 } ;
  
  vigra::MultiArray < trg_dim , crd_type >
    coordinates ( trg_shape ) ;
    
  vigra::MultiArray < 1 , rc_type >
    crd_1d { 101 } ;
  
  vigra::MultiArray < trg_dim , dtype >
    target ( trg_shape ) ;
  
  vigra::MultiArray < trg_dim , double >
    d_target ( trg_shape ) ;
  
  cf_shape_t gcf_shape { 101 } ;
  vigra::MultiArray < cf_dim , dtype >
    ge_target ( gcf_shape ) ;
  
  vigra::MultiArray < cf_dim , dtype >
    ge_target_2 ( gcf_shape ) ;
  
  // produce random original data. We produce some very high values
  // here, since we want to use them also for testing processing of
  // integral data, where we want to exhaust the dynamic range of the
  // integral type, since we can't have postcomma digits. The results
  // for processing floats will therefore have errors in the order of
  // magnitude of 1, rather than the usual errors around 1e-5 for
  // float data, when the test data are small numbers in [0,1]
    
  std::random_device rd ;
  std::mt19937 gen ( rd() ) ;
  std::uniform_real_distribution<> dis ( - value_extent , value_extent ) ;
  auto data_ele_view = original.expandElements ( 0 ) ;
  for ( auto & e : data_ele_view )
    e = dis ( gen ) ;

  // produce random coordinates. Note that many of these coordinates
  // will be out-of-range. Thy are folded into the range by the
  // 'safe evaluator', so this feature is also tested.
  // we deliberately overshoot the defined range by more than the first
  // reflection to be sure the extrapolation works as intended.
  
  std::uniform_real_distribution<> crd_dis ( -300.0 , 300.0 ) ;
  auto crd_ele_view = coordinates.expandElements ( 0 ) ;
  for ( auto & e : crd_ele_view )
    e = crd_dis ( gen ) ;
  
  // produce a set of 'safe' 1D coordinates to test grid eval
  // grid eval can't handle out-of-range coordinates, we have to
  // stay within the spline's safe range.
  
  std::uniform_real_distribution<> crd1_dis ( 98.0 ) ;
  for ( auto & e : crd_1d )
    e = crd1_dis ( gen ) ;

  grid_spec < cf_dim , rc_type > grid ;
  for ( int d = 0 ; d < cf_dim ; d++ )
    grid [ d ] = crd_1d ;
  
  // create a b-spline over the data, prefilter it and create
  // an evaluator. note how we pass rc_type to make_safe_evaluator.
  // we test with the same boundary conditions along all axes.
  
  vigra::TinyVector < vspline::bc_code , cf_dim > bcv ( bc ) ;
  
  typedef vspline::bspline < dtype , cf_dim > spline_type ;
  spline_type bspl ( cf_shape ,
                     spline_degree ,
                     bcv ) ;
  
  // we instantiate prefilter with the given vsz to make sure that
  // we can indeed pick arbitrary vectorization widths.
  
  bspl.template prefilter < dtype , math_ele_type , vsz >
    ( original ) ;
  
  // we instantiate ev with the given vsz to make sure that
  // we can indeed pick arbitrary vectorization widths
    
  enum { ev_vsz = vsz } ;
  
  auto ev = vspline::make_safe_evaluator
            < spline_type , rc_type , ev_vsz , math_ele_type >
              ( bspl ) ;

  // run an array-based transform with the coordinates  
  
  vspline::transform ( ev , coordinates , target ) ;
  
  // check the result against single-value evaluation.
  // here we expect to see precisely equal results.
  // note that in this test, we can't know what the
  // correct result should be, we only make sure that
  // vectorized evaluation and unvectorized evaluation
  // produce near identical results.
  
  auto it = target.begin() ;
  
  for ( auto const & e : coordinates )
  {
    // TODO when working on int coefficients, I don't get total equality
    assert ( condense ( *it - ev ( e ) ) <= tolerance ) ;
    ++it ;
  }
  
  // create a safe evaluator with specified target type 'double'
  // and process the coordinates in 'coordinates' with it. Check
  // single-eval results against the result of 'transform'.
  
  auto evd = vspline::make_safe_evaluator
            < spline_type , rc_type , ev_vsz ,
              math_ele_type , double >
              ( bspl ) ;

  vspline::transform ( evd , coordinates , d_target ) ;
  
  auto itd = d_target.begin() ;  
  
  for ( auto const & e : coordinates )
  {
    auto cnd = condense ( *itd - evd ( e ) ) ;
    if ( cnd > tolerance )
      std::cout << "**** consensed: " << cnd
                << " gt. tolerance " << tolerance << std::endl ;
    assert ( cnd <= tolerance ) ;
    ++itd ;
  }
  
  // create an evaluator over the spline to pass it to grid eval.
  // here we use a 'raw' evaluator, not the type make_safe_evaluator
  // or make_evaluator would produce.
  
  auto gev = vspline::evaluator
             < crd_type , dtype , ev_vsz , -1 , math_ele_type >
               ( bspl ) ;
            
  // call grid-based transform, first with 'ev', which is
  // of vspline::grok_type, next with 'gev', which is of type
  // vspline::evaluator and uses different (and faster) code.
  
  vspline::transform ( grid , ev , ge_target ) ; 
  vspline::transform ( grid , gev , ge_target_2 ) ; 

  // Now we 'manually' iterate over the coordinates in the grid
  // and compare the result to the content of ge_target to make sure
  // grid eval has worked correctly.
  
  crd_type c ;
  auto & cc = reinterpret_cast
       < vigra::TinyVector < rc_type , cf_dim > & > ( c ) ;
       
  vigra::MultiCoordinateIterator < cf_dim > mci ( gcf_shape ) ;
  auto itg2 = ge_target_2.begin() ;
  
  for ( auto & ref : ge_target )
  {
    // fill in the coordinate

    for ( int d = 0 ; d < cf_dim ; d++ )
      cc [ d ] = grid [ d ] [ (*mci) [ d ] ] ;
   
    // evaluate at the coordinate and compare to grid eval's output.
    // we expect equality of results here, since the arithmetic
    // is near identical for both ways of generating the result
    // TODO fails with -Ofast or -ffast-math, hence less strict

//     assert ( ref == gev ( c ) ) ;
//     assert ( ref == *itg2 ) ;

    assert ( condense ( ref - gev(c) ) <= tolerance ) ;
    assert ( condense ( ref - *itg2 ) <= tolerance ) ;

    ++ mci ; 
    ++itg2 ;
  }
  
  // restore original data, first by using an index-based
  // transform, then by calling 'restore' (using convolution)
  // directly on the spline, producing the restored data in
  // the spline's core.
  // Note how when working on integral data, the index-based
  // transform produced wildly wrong results due to the fact
  // that the evaluation was done entirely in the integral type.
  // The new evaluation code can, if math_ele_type is real,
  // produce reasonable results which only suffer from
  // quantization errors, which was only possible with
  // restoration by convolution in the last version.
  
  vspline::transform ( ev , restored ) ;
  
  // we instantiate restore with the given vsz to make sure that
  // we can indeed pick arbitrary vectorization widths.
  
  vspline::restore < cf_dim , dtype , math_ele_type , vsz > ( bspl ) ;
  
  // compare the data obtained from the two restoration methods.
  // due to the slightly different arithmetic involved, we expect
  // similar, but nor necessarily identical results.
  
  check_diff ( original , restored , value_extent ) ;
  check_diff ( original , bspl.core , value_extent ) ;
}

template < typename vsz_t ,
           typename chn_t ,
           typename trg_dim_t ,
           typename cf_dim_t ,
           typename ... other_types >
void doit()
{
  enum { vsz = vsz_t::value } ;
  enum { chn = chn_t::value } ;
  enum { trg_dim = trg_dim_t::value } ;
  enum { cf_dim = cf_dim_t::value } ;
  
  std::cout << "vsz " << vsz
            << " chn " << chn
            << " trg_dim " << trg_dim
            << " cf_dim " << cf_dim << std::endl ;
  
  // we do a few exemplary runs - if we tried to exhaust all
  // possible combinations of data types, spline degrees and
  // boundary conditions, this would take forever, and we already
  // have code in restore_test exploring that way. Here we're
  // more interested in making sure that 1D and nD data and
  // real and integral types are processed as expected and that
  // all 'high-level' capabilities of vspline are invoked.
  
  // a very 'standard' scenario: cubic b-spline, all in float:
  
  std::cout << "float crd, float value " << std::endl ;
  
  test < cf_dim , trg_dim , chn , vsz , float , float , float >
    ( 3 , vspline::MIRROR ) ;

  // testing an integral-valued spline
    
  std::cout << "float crd, short value " << std::endl ;
  
  test < cf_dim , trg_dim , chn , vsz , float , short , float >
    ( 2 , vspline::MIRROR ) ;
  
  std::cout << "float crd, int value " << std::endl ;
  
  test < cf_dim , trg_dim , chn , vsz , float , int , float >
    ( 4 , vspline::PERIODIC ) ;
  
  // we can use long-valued coefficients, but evaluation for these
  // types won't be vectorized with Vc; vspline will use it's own
  // fallback type simd_tv
    
  std::cout << "double crd, long value " << std::endl ;
  
  test < cf_dim , trg_dim , chn , vsz , float , long , double >
    ( 3 , vspline::PERIODIC ) ;
    
  // all data in double. This will be vectorized, so it's quite
  // fast, and yet very precise.
    
  std::cout << "double crd, double value " << std::endl ;
  
  test < cf_dim , trg_dim , chn , vsz , double , double , double >
    ( 5 , vspline::REFLECT ) ;
    
  // finally we test with all data in long double. This won't
  // be vectorized at all, but it's extremely precise.
    
  std::cout << "long double crd, long double value " << std::endl ;
  
  test < cf_dim , trg_dim , chn , vsz ,
         long double , long double , long double >
    ( 3 , vspline::NATURAL ) ;
}

// when choosing vsize, the vectorization width, we don't go through
// all values from 1 to 32, but only pick a few representative ones.
// most of the time, vsize won't be picked 'manually' - and vectorization
// widths which aren't a multiple of the hardware vector size are
// quite futile, but here the point is to see if the code still
// performs correctly even with 'odd' vsize values.
// TODO: it doesn't - 3 fails below

template < typename ... other_types >
void set_vsz ( int vsz )
{
  switch ( vsz )
  {
    case 1:
      doit < ast(1) , other_types ... >() ;
      break ;
    case 3:
      doit < ast(3) , other_types ... >() ;
      break ;
    case 4:
      doit < ast(4) , other_types ... >() ;
      break ;
    case 8:
      doit < ast(8) , other_types ... >() ;
      break ;
    case 16:
      doit < ast(16) , other_types ... >() ;
      break ;
    default:
      break ;
  }
}

template < typename ... other_types >
void set_chn ( int chn , int vsz )
{
  switch ( chn )
  {
    case 1:
      set_vsz < ast(1) , other_types ... > ( vsz ) ;
      break ;
    case 2:
      set_vsz < ast(2) , other_types ... > ( vsz ) ;
      break ;
    case 3:
      set_vsz < ast(3) , other_types ... > ( vsz ) ;
      break ;
    default:
      break ;
  }
}

template < typename ... other_types >
void set_trg_dim ( int trg_dim , int chn , int vsz )
{
  switch ( trg_dim )
  {
    case 1:
      set_chn < ast(1) , other_types ... > ( chn , vsz ) ;
      break ;
    case 2:
      set_chn < ast(2) , other_types ... > ( chn , vsz ) ;
      break ;
    default:
      break ;
  }
}

template < typename ... other_types >
void set_cf_dim ( int cf_dim , int trg_dim , int chn , int vsz )
{
  switch ( cf_dim )
  {
    case 1:
      set_trg_dim < ast(1) , other_types ... > ( trg_dim , chn , vsz ) ;
      break ;
    case 2:
      set_trg_dim < ast(2) , other_types ... > ( trg_dim , chn , vsz ) ;
      break ;
    default:
      break ;
  }
}

#define CF_DIM_MAX 2
#define TRG_DIM_MAX 2
#define CHN_MAX 2

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);
  
  for ( int cf_dim = 1 ; cf_dim <= CF_DIM_MAX ; cf_dim++ )
  {
    for ( int trg_dim = 1 ; trg_dim <= TRG_DIM_MAX ; trg_dim++ )
    {
      for ( int chn = 1 ; chn <= CHN_MAX ; chn++ )
      {
        set_cf_dim ( cf_dim , trg_dim , chn , 1 ) ;
        set_cf_dim ( cf_dim , trg_dim , chn , 3 ) ;
        set_cf_dim ( cf_dim , trg_dim , chn , 8 ) ;
        set_cf_dim ( cf_dim , trg_dim , chn , 16 ) ;
      }
    }
  }

  std::cout << "terminating" << std::endl ;
}