/************************************************************************/
/*                                                                      */
/*    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 eval.cc
///
/// \brief simple demonstration of creation and evaluation of a b-spline
///
/// takes a set of knot point values from std::cin, calculates a 1D b-spline
/// over them, and evaluates it at coordinates taken from std::cin.
/// The output shows how the coordinate is split into integral and real
/// part and the result of evaluating the spline at this point.
/// Note how the coordinate is automatically folded into the defined range.
///
/// Two evaluations are demonstrated, using unvectorized and vectorized
/// input/output.
///
/// Since this is a convenient place to add code for testing evaluation
/// speed, you may pass a number on the command line. eval.cc will then
/// perform as many evaluations and print the time it took.
/// The evaluations will be distributed to several threads, so there is
/// quite a bit of overhead; pass numbers from 1000000 up to get realistic
/// timings.
///
/// compile: ./examples.sh eval.cc
///
/// compare speed test results from different compilers/back-ends, using
/// 'echo' to provide the arguments normally taken from cin. This invocation
/// runs the test with a periodic degree-19 spline over knot points 0 and 1
/// and only shows the name of the binary and the result of the speed test:
///
/// for f in *eval*++; do echo $f;echo 19 2 0 1 . .5 | ./$f 100000000 | grep took; done
///
/// this will produce output like:
/// hwy_eval_clang++
/// 100000000 evaluations took 2229 ms
/// hwy_eval_g++
/// 100000000 evaluations took 2527 ms
/// stds_eval_clang++
/// 100000000 evaluations took 3401 ms
/// stds_eval_g++
/// 100000000 evaluations took 2351 ms
/// vc_eval_clang++
/// 100000000 evaluations took 3281 ms
/// vc_eval_g++
/// 100000000 evaluations took 3385 ms
/// vs_eval_clang++
/// 100000000 evaluations took 3883 ms
/// vs_eval_g++
/// 100000000 evaluations took 5479 ms
/// 
/// note how both the back-end and the compiler make a difference, highway
/// and clang++ currently coming out on top on my system. This is a moving
/// target, as both the back-ends and the compilers evolve.

#include <iostream>
#include <iomanip>
#include <ctime>
#include <chrono>

#include <vspline/vspline.h>

// you can use float, but then can't use very high spline degrees.
// If you use long double, the code won't be vectorized in hardware.

typedef double dtype ;
typedef double rc_type ;

enum { vsize = vspline::vector_traits<dtype>::vsize } ;

// speed_test will perform as many vectorized evaluations as the
// range it receives spans. The evaluation is always at the same place,
// trying to lower all memory access effects to the minimum.

template < class ev_type >
void speed_test ( ev_type & ev , std::size_t nr_ev )
{
  typename ev_type::in_v in = dtype(.111) ;
  typename ev_type::out_v out ;
  
  while ( nr_ev-- )
    ev.eval ( in , out ) ;
}

int main ( int argc , char * argv[] )
{
  long TIMES = 0 ;
  
  if ( argc > 1 )
    TIMES = std::atoi ( argv[1] ) ;
    
  // get the spline degree and boundary conditions from the console

  std::cout << "enter spline degree: " ;
  int spline_degree ;
  std::cin >> spline_degree ;
  
  int bci = -1 ;
  vspline::bc_code bc ;
  
  while ( bci < 1 || bci > 4 )
  {
    std::cout << "choose boundary condition" << std::endl ;
    std::cout << "1) MIRROR" << std::endl ;
    std::cout << "2) PERIODIC" << std::endl ;
    std::cout << "3) REFLECT" << std::endl ;
    std::cout << "4) NATURAL" << std::endl ;
    std::cin >> bci ;
  }
  
  switch ( bci )
  {
    case 1 :
      bc = vspline::MIRROR ;
      break ;
    case 2 :
      bc = vspline::PERIODIC ;
      break ;
    case 3 :
      bc = vspline::REFLECT ;
      break ;
    case 4 :
      bc = vspline::NATURAL ;
      break ;
  }
  
  // obtain knot point values

  dtype v ;
  std::vector<dtype> dv ;
  std::cout << "enter knot point values (end with single '.')" << std::endl ;
  while ( std::cin >> v )
    dv.push_back ( v ) ;

  std::cin.clear() ;
  std::cin.ignore() ;
  
  // fix the type for the bspline object
  
  typedef vspline::bspline < dtype , 1 > spline_type ;
  spline_type bsp ( dv.size() , spline_degree , bc ) ;
  
  std::cout << "created bspline object:" << std::endl << bsp << std::endl ;

  // fill the data into the spline's 'core' area
  
  for ( size_t i = 0 ; i < dv.size() ; i++ )
    bsp.core[i] = dv[i] ;

  // prefilter the data
  
  bsp.prefilter() ;
  
  std::cout << std::fixed << std::showpoint << std::setprecision(12) ;
  std::cout << "spline coefficients (with frame)" << std::endl ;
  for ( auto& coeff : bsp.container )
    std::cout << " " << coeff << std::endl ;

  // obtain a 'safe' evaluator which folds incoming coordinates
  // into the defined range
  
  auto ev = vspline::make_safe_evaluator ( bsp ) ;
  typedef  vspline::evaluator < int , double > iev_t ;
  auto iev = iev_t ( bsp ) ;

  int ic ;
  rc_type rc ;
  dtype res ;

  std::cout << "enter coordinates to evaluate (end with EOF)"
            << std::endl ;
  while ( std::cin >> rc )
  {
    if ( rc < bsp.lower_limit(0) || rc > bsp.upper_limit(0) )
    {
      std::cout << "warning: " << rc
                << " is outside the spline's defined range."
                << std::endl ;
      std::cout << "using automatic folding to process this value."
                << std::endl ;
    }

    // evaluate the spline at this location

    ev.eval ( rc , res ) ;

    std::cout << rc << " -> " << res << std::endl ;
    
    if ( TIMES )
    {
      std::chrono::system_clock::time_point start
        = std::chrono::system_clock::now() ;

      // for the speed test we build a plain evaluator; we'll
      // only be evaluating the spline near 0.0 repeatedly, so
      // we don't need folding into the safe range. We're not
      // fixing 'specialize', so evaluation will be general
      // b-spline evaluation, even for degree 0 and 1 splines.

      auto ev = vspline::evaluator < dtype , dtype , vsize > ( bsp ) ;

      std::size_t share = ( TIMES / vsize ) / vspline::default_njobs ;

      std::function < void() > payload =
      [&]() { speed_test ( ev , share ) ; } ;

      vspline::multithread ( payload , vspline::default_njobs ) ;
  
      std::chrono::system_clock::time_point end
        = std::chrono::system_clock::now() ;

      std::cout << TIMES << " evaluations took "
                << std::chrono::duration_cast<std::chrono::milliseconds>
                    ( end - start ) . count()
                << " ms" << std::endl ;
    }
  }
}