/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2017 - 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 n_shift.cc
///
/// \brief fidelity test
/// 
/// This is a test to see how much a signal degrades when it is submitted
/// to a sequence of operations:
///
/// - create a b-spline over the signal
/// - evaluate the spline at unit-spaced locations with an arbitrary offset
/// - yielding a shifted signal, for which the process is repeated
///
/// Finally, a last shift is performed which samples the penultimate version
/// of the signal at points coinciding with coordinates 0...N-1 of the
/// original signal. This last iteration should ideally recreate the
/// original sequence.
///
/// The test is done with a periodic signal to avoid margin effects.
/// The initial sequence is created by evaluating a periodic high-degree
/// b-spline of half the size at steps n/2. This way we start out
/// with a signal with low high-frequency content - a signal which can
/// be approximated well with a b-spline. Optionally, the supersampling
/// factor can be passed on the command line to experiment with different
/// values than the default of 2. Supersampling factors should be whole
/// numbers (halves also work, but less precise) - so that the knot
/// points of the original signal coincide with knot points of the
/// supersampled signal. This way, every interval contains a partial
/// polynomial with no discontinuities, and the spline can faithfully
/// represent the signal.
/// Both this program and bls.cpp show that there is a significant
/// threshold involved: if the initial signal contains *no frequencies
/// above half the Nyquist frequency*, a b-spline of sufficiently
/// large degree becomes *immune to resampling*. See also bls.cpp,
/// which produces test signals using IFFT, which is a clearer way
/// of excluding part of the frequency spectrum, because the signal
/// generated by IFT can - by definition - only contain frequencies
/// which have non-zero corresponding coefficients in the FT.
///
/// compile with: clang++ -pthread -O3 -std=c++11 n_shift.cc -o n_shift
///
/// invoke like: n_shift 17 500

#include <iostream>
#include <random>
#include <vigra/accumulator.hxx>
#include <vigra/multi_math.hxx>
#include <vspline/vspline.h>

int main ( int argc , char * argv[] )
{
  if ( argc < 3 )
  {
    std::cerr << "pass the spline's degree, the number of iterations" 
              << std::endl
              << "and optionally the supersampling factor"
              << std::endl ;
    exit ( -1 ) ;
  }

  int degree = std::atoi ( argv[1] ) ;
  
  assert ( degree >= 0 && degree <= vspline_constants::max_degree ) ;
  
  int iterations = 1 + std::atoi ( argv[2] ) ;
  
  const int sz = 1024 ;
  
  long double widen = 2.0 ;
  
  if ( argc > 3 )
    widen = atof ( argv[3] ) ;
  
  assert ( widen >= 1.0 ) ;
  
  int wsz = sz * widen ;
  
  vigra::MultiArray < 1 , long double > original ( wsz ) ;
  vigra::MultiArray < 1 , long double > target ( wsz ) ;
  
  // we start out by filling the first bit of 'original' with random data
  // between -1 and 1
  
  std::random_device rd ;
  std::mt19937 gen ( rd() ) ;
//   gen.seed ( 1 ) ;              // if desired, level playing field
  std::uniform_real_distribution<> dis ( -1 , 1 ) ;
  for ( int x = 0 ; x < sz ; x++ )
    original [ x ] = dis ( gen ) ;

  // create the bspline object to produce the data we'll work with
  
  vspline::bspline < long double ,   // the spline's data type
                     1 >             // one dimension only
    bsp ( sz ,                  // sz values
          20 ,                  // high degree for smoothness
          vspline::PERIODIC ,   // periodic boundary conditions
          0.0 ) ;               // no tolerance
          
  vigra::MultiArrayView < 1 , long double > initial
    ( vigra::Shape1(sz) , original.data() ) ;
  
  // pull in the data while prefiltering
    
  bsp.prefilter ( initial ) ;
  
  // create an evaluator to obtain interpolated values
  
  typedef vspline::evaluator < long double , long double > ev_type ;
  
  ev_type ev ( bsp ) ; // from the bspline object we just made
  
  // now we evaluate at 1/widen steps, into target
  
  for ( int x = 0 ; x < wsz ; x++ )
    target [ x ] = ev ( (long double) ( x ) / (long double) ( widen ) ) ;
  
  // we take this as our original signal. Since this is a sampling
  // of a periodic signal (the spline in bsp) representing a full
  // period, we assume that a b-spline over this signal will, within
  // the spline's capacity, approximate the signal in 'original'.
  // what we want to see is how sampling at offsetted positions
  // and recreating a spline over the offsetted signal will degrade
  // the signal with different-degree b-splines and different numbers
  // of iterations.

  original = target ;
  
  // now we set up the working spline

  vspline::bspline < long double ,    // spline's data type
                     1 >         // one dimension
    bspw ( wsz ,                 // wsz values
           degree ,              // degree as per command line
           vspline::PERIODIC ,   // periodic boundary conditions
           0.0 ) ;               // no tolerance

  // we pull in the working data we've just generated

  bspw.prefilter ( original ) ;
  
  // and set up the evaluator for the test
  
  ev_type evw ( bspw ) ;
  
  // we want to map the incoming coordinates into the defined range.
  // Since we're using a periodic spline, the range is from 0...N,
  // rather than 0...N-1 for non-periodic splines

  auto gate = vspline::periodic ( 0.0L , (long double)(wsz) ) ;

  // now we do a bit of functional programming.
  // we chain gate and evaluator:

  auto periodic_ev = gate + evw ;
  
  // we cumulate the offsets so we can 'undo' the cumulated offset
  // in the last iteration
    
  long double cumulated_offset = 0.0 ;
  
  for ( int n = 0 ; n < iterations ; n++ )
  {
    using namespace vigra::multi_math ;
    using namespace vigra::acc;
  
    // use a random, largish offset (+/- 1000). any offset
    // will do, since we have a periodic gate, mapping the
    // coordinates for evaluation into the spline's range
    
    long double offset = 1000.0 * dis ( gen ) ;
    
    // with the last iteration, we shift back to the original
    // 0-based locations. This last shift should recreate the
    // original signal as best as a spline of this degree can
    // do after so many iterations.

    if ( n == iterations - 1 )
      offset = - cumulated_offset ;

    cumulated_offset += offset ;

    if ( n > ( iterations - 10 ) )
      std::cout << "iteration " << n << " offset " << offset
                << " cumulated offset " << cumulated_offset << std::endl ;
    
    // we evaluate the spline at unit-stepped offsetted locations,
    // so, 0 + offset , 1 + offset ...
    // in the last iteration, this should ideally reproduce the original
    // signal.
    
    for ( int x = 0 ; x < wsz ; x++ )
    {
      auto arg = x + offset ;
      target [ x ] = periodic_ev ( arg ) ;
    }
    
    // now we create a new spline over target, reusing bspw
    // note how this merely changes the coefficients of the spline,
    // the container for the coefficients is reused, and therefore
    // the evaluator (evw) will look at the new set of coefficients.
    // So we don't need to create a new evaluator.

    bspw.prefilter ( target ) ;
    
    // to convince ourselves that we really are working on a different
    // sampling of the signal signal - and to see how close we get to the
    // original signal after n iterations, when we use an offset to get
    // the sampling locations back to 0, 1, ...
      
    vigra::MultiArray < 1 , long double > error_array
      ( vigra::multi_math::squaredNorm ( target - original ) ) ;

    AccumulatorChain < long double , Select < Mean, Maximum > > ac ;
    extractFeatures ( error_array.begin() , error_array.end() , ac ) ;
    
    if ( n > ( iterations - 10 ) )
    {
      if ( n == iterations - 1 )
        std::cout << "final result, evaluating at original unit steps"
                  << std::endl ;
      std::cout << "signal difference Mean:    "
                << sqrt(get<Mean>(ac)) << std::endl;
      std::cout << "signal difference Maximum: "
                << sqrt(get<Maximum>(ac)) << std::endl;
    }
  }
}