/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2018 - 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 bls.cpp
///
/// \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 periodic 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, and you'll see a cumulated offset of zero for this
/// last iteration.
///
/// The test is done with a periodic signal to avoid margin effects.
/// The initial sequence is created by performing an IFFT on random
/// data up to a given cutoff frequency, producing a band-limited signal.
/// Alternatively, by modifying the initial values which are put into
/// the frequency-domain representation of the test signal, more specific
/// scenarios can be investigated (try setting the values to 1 uniformly etc)
/// With the cutoff factor approaching 1.0, the signal contains more and
/// higher high frequencies, making it ever more difficult for the spline
/// to represent the data faithfully. With a cutoff of <= .5, fidelity is
/// very high for high-degree b-splines; even after many iterations
/// the spline is nearly unchanged and faithfully represents the
/// original signal.
///
/// There are two factors to observe: high-degree splines cope better with
/// high frequency components, but they suffer from the large magnitude of
/// coefficient values, making the evaluation error larger - especially since
/// high-degree splines have wide support.
///
/// An important conclusion is that even though all input signals produced
/// with this test are band-limited by design, this quality is not sufficient
/// to produce a b-spline from them which is 'stable' in the sense of being
/// immune to the multiple shifting operation, unless the signal is free from
/// frequencies above *half the Nyquist frequency*.
///
/// If the input contains high frequencies, and yet a faithful b-spline
/// representation has to be provided which is immune to shifting, one
/// solution is to create a high-degree spline over the signal, sample it
/// at half unit steps and use the resulting signal to create the 'working'
/// spline. By creating the high-order 'super-spline', a signal is created
/// which is low in high-frequency content and can be faithfully represented
/// by the working spline. The upsampling can also be done in any other
/// way, like FFT/IFFT. See n_shift.cc for an implementation of upsampling
/// using a high-degree b-spline as 'super-spline'.
///
/// An aside: this program shows that we can easily shift and resample a
/// high-degree periodic spline for thousands of times without much signal
/// degradation, which is quite remarkable, but requires the signal to be
/// band-limited to *half the Nyquist frequency*. But if we look at the
/// fourier transformation of each of the shifted variants of the signal,
/// we don't - as we might expect - get precisely the same magnitudes for
/// the partial frequencies. I'm not quite sure what to make of this.
///
/// note that this example requires fftw3.
///
/// compile with: clang++ -O3 -std=c++11 -obls bls.cpp \
///                       -pthread -lvigraimpex -lfftw3
///
/// to use highway, add -DUSE_HWY and -lhwy
/// to use Vc, add -DUSE_VC and -lVc
/// to use std::simd, use -DUSE_STDSIMD and --std=c++17
// TODO: comes out wrong with clang++ and std::simd on my debian11 system
///
/// invoke with: bls <spline degree> <number of iterations> \
///                 [ <frequency cutoff> ]
///
/// Try different spline degrees and cutoff frequencies. With cutoff
/// frequencies of 0.5 or less you'll notice that even with relatively
/// low spline degree, the 'immunity' to resampling is quite good: even
/// with thousands of iterations, the original signal can be regained
/// with minimal error. Raising the spline degree, the cutoff frequency
/// can even be raised above 0.5 and the spline remains quite 'immune'.

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

// a functor to add an offset to an incoming value. This is
// used further down to create a combined functor producing
// resamplings of the initial signal.

struct offset_f : public vspline::unary_functor < double >
{
  const double offset ;

  offset_f ( const double & _offset )
  : offset ( _offset )
  { }

  template < typename IN , typename OUT >
  void eval ( const IN & in , OUT & out ) const
  {
    out = in + offset ;
  }
} ;

int main ( int argc , char * argv[] )
{
  std::cout << std::fixed << std::showpoint
            << std::setprecision
                (std::numeric_limits<double>::max_digits10) ;

  if ( argc < 3 )
  {
    std::cerr << "pass the spline's degree and the number of iterations" 
              << std::endl
              << "and, optionally, the cutoff frequency"
              << std::endl ;
    exit ( -1 ) ;
  }

  int degree = std::atoi ( argv[1] ) ;
  
  assert ( degree >= 0 && degree <= vspline_constants::max_degree ) ;
  
  int iterations = std::max ( 1 , std::atoi ( argv[2] ) ) ;
  
  double f_cutoff = .5 ;
  if ( argc == 4 )
  {
    f_cutoff = std::atof ( argv[3] ) ;
    std::cout << "using frequency cutoff " << f_cutoff << std::endl ;
  }
  
  const std::size_t sz = 1024 ;
  
  vigra::MultiArray < 1 , double > signal ( sz ) ;
  vigra::MultiArray < 1 , double > target ( sz ) ;

  vigra::MultiArray < 1 , vigra::FFTWComplex<double> >
    fourier ( signal.shape() / 2 + 1 ) ;

  std::random_device rd ;
  std::mt19937 gen ( rd() ) ;
  gen.seed ( 42 ) ;              // level playing field
  std::uniform_real_distribution<> dis ( -1 , 1 ) ;
  
  // we fill the coefficients in, filling in random values, until
  // we reach the cutoff point, then set the higher ones to zero.

  int fill = 0 ;
  for ( auto & e : fourier )
  {
    if ( fill < f_cutoff * sz / 2.0 )
      e = vigra::FFTWComplex<double> ( dis(gen) , dis(gen) ) ;
    else
      e = vigra::FFTWComplex<double> ( 0.0 , 0.0 ) ;
    ++fill ;
  }
  
  // an inverse FFT produces the signal we'll run the test on

  vigra::fourierTransformInverse ( fourier , signal ) ;

  // now we set up the working spline

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

  // we 'pull in' the signal data we've just generated via the spline's
  // 'prefilter' function which accepts an array of knot points

  bspw.prefilter ( signal ) ;
  
  // the next few lines are just statistics to better judge what's going on

  using namespace vigra::multi_math ;
  using namespace vigra::acc;
  
  vigra::MultiArray < 1 , double > error_array
    ( vigra::multi_math::squaredNorm ( bspw.core ) ) ;

  AccumulatorChain < double , Select < Mean, Maximum > > ac ;
  extractFeatures ( error_array.begin() , error_array.end() , ac ) ;
  {
    std::cout << "coefficients Mean:    "
              << sqrt(get<Mean>(ac)) << std::endl;
    std::cout << "coefficients Maximum: "
              << sqrt(get<Maximum>(ac)) << std::endl;
  }
  
  // now we create an evaluator to obtain interpolated values,
  // using 'make_safe_evaluator' which will create a periodized
  // spline because we've specified PERIODIC boundary conditions.
  // So the resulting evaluator will accept arbitrary coordinates
  // and map them to the defined range.
  
  auto ev = vspline::make_safe_evaluator ( bspw ) ;

  // we cumulate the offsets so we can 'undo' the cumulated offset
  // in the last iteration
    
  double cumulated_offset = 0.0 ;
  
  for ( int n = 0 ; n < iterations ; n++ )
  {
    // use a random, largish offset (+/- 1000). any offset
    // will do, since we have a periodic spline, mapping the
    // coordinates for evaluation into the spline's range
    
    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. Note the functional
    // composition: offset_f creates the offsetting functor, adding the
    // given offset to incoming discrete coordinates, ev is the periodic
    // spline functor we made before, and 'chaining' the two with '+'
    // creates a combined functor which feeds the first functor's output
    // to the second functor's input. So the resulting functor 'f' combines
    // three steps: add an offset to an incoming coordinate, map the
    // result into the spline's defined range, and evaluate it there.

    auto f = offset_f ( offset ) + ev ;

    // use vspline::transform to 'roll out' the combined functor.
    // vspline::transform with a single array argument 'feeds' discrete
    // coordinates to the functor. This is such a common operation that
    // vspline::transform has a specific overload for the purpose. The
    // discrete coordinates are simply those coordinates at which the
    // values calculated by the functor should be stored.

    vspline::transform ( f , target ) ;

    // 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 (ev) 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 - and to see how close we get to the
    // original signal after n iterations, when we use a last offset to get
    // the sampling locations back to 0, 1, ..., before the 'final' result,
    // we echo the statistics for the last ten iterations. The final
    // iteration resamples the result of the penultimate one so that it
    // should recreate the original signal, if the spline is indeed
    // 'immune' to resampling. This will be the case if the spline is
    // of sufficiently large degree, and the signal is band-limited to half
    // the Nyquist frequency.
      
    vigra::MultiArray < 1 , double > error_array
      ( vigra::multi_math::squaredNorm ( target - signal ) ) ;

    AccumulatorChain < 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 signal 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;
    }
  }
}