/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2015 - 2018 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 bootstrap.cc
///
/// \brief Code to calculate b-spline basis function values at half
/// unit steps, and prefilter poles.
///
/// These calculations are done using GNU GSL, BLAS and GNU GMP
/// which are not used in other parts of vspline. this file also has
/// code to verify that the precomputed data in <vspline/poles.h>
/// are valid. TODO make that a separate program
///
/// In contrast to my previous implementation of generating the
/// precomputed values (prefilter_poles.cc, now gone), I don't
/// compute the values one-by-one using recursive functions, but
/// instead build them layer by layer, each layer using it's
/// predecessor. This makes the process very fast, even for
/// large spline degrees.
///
/// In my tests, I found that with the very precise basis function
/// values obtained with GMP, I could get GSL to provide raw values
/// for the prefilter poles up to degree 45, above which it would
/// miss roots. So this is where I set the limit, and I think it's
/// unlikely anyone would ever want more than degree 45. If they
/// do, they'll have to find other polynomial root code.
///
/// compile: g++ -O3 -std=gnu++11 bootstrap.cc -obootstrap \
///              -lgmpxx -lgmp -lgsl -lblas -pthread
///
/// run by simply issuing './bootstrap' on the command line. you may
/// want to redirect output to a file. The output should be equal
/// to the values in vspline/poles.h

#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>
#include <numeric>
#include <assert.h>
#include <gsl/gsl_poly.h>
#include <gmpxx.h>

// #include <eigen3/unsupported/Eigen/Polynomials>

#include <vspline/vspline.h>

// In my trials, I could take the spline degree up to 45 when calculating
// the prefilter poles. The basis function value can be found for arbitrary
// degrees.

#define degree_limit 46

// when using mpf_class from GMP, we use 512 bits of precision.

#define mpf_bits 512

// we want to build a table of basis function values at multiples
// of 1/2 mimicking the Cox-de Boor recursion. We only need one half
// of the table since the basis function is symmetric (b(x) == b(-x))
// - except for degree 0, which needs special treatment.
// Note that we 'invert' the Cox-de Boor recursion and generate
// the result iteratively, starting with degree 0. This results in
// every single difference relation being executed excactly once,
// which makes the iterative solution very fast.
// Note the use of x2 throughout, this is the doubled argument,
// so that we can operate in integral maths throughout the basis
// function calculations (we use mpq_class for arbitrary-sized
// fractions, so there is no error at all)
// We make the table one column wider than strictly necessary to
// allow access to x2 outside the range, in which case we want to
// access a location containing 0.

mpq_class basis_table [ degree_limit ] [ degree_limit + 1 ] ;

// we want to access the table both for reading and writing, mapping
// negative x2 to positive ones and treating access to the first row
// (for degree 0) specially. Once the table is filled, we can access
// the table at arbitrary x2.

mpq_class & access_table ( int degree , int x2 )
{
  if ( degree >= degree_limit || degree < 0 )
    throw std::invalid_argument ( "access with degree out of bounds" ) ;

  int column ;

  // at degree 0, we have no symmetry: the basis function is
  // 1 for x2 == -1 and x2 == 0.
  if ( degree == 0 )
  {
    if ( x2 == -1 || x2 == 0 )
      column = 0 ;
    else
      column = 1 ;
  }
  // for all other degrees the basis function is symmetric
  else
    column = x2 >= 0 ? x2 : -x2 ;

  // for out-of-range x2, land in the last column, which is 0.
  column = std::min ( column , degree_limit ) ;
  
  return basis_table [ degree ] [ column ] ; 
}

// subroutine filling the table at a specific x2 and degree.
// we rely on the table being filled correctly up to degree - 1

void _fill_basis_table ( int x2 , int degree )
{
  // look up the neighbours one degree down and one to the
  // left/right

  mpq_class & f1 ( access_table ( degree - 1 , x2 + 1 ) ) ;
  mpq_class & f2 ( access_table ( degree - 1 , x2 - 1 ) ) ;
  
  // use the recursion to calculate the current value and store
  // that value at the appropriate location in the table

  access_table ( degree , x2 )
    =   (   f1 * ( degree + 1 + x2 )
          + f2 * ( degree + 1 - x2 ) )
      / ( 2 * degree ) ;
}

// to fill the whole table, we initialize the single value
// in the first row, then fill the remainder by calling the
// single-value subroutine for all non-zero locations

void fill_basis_table()
{
  access_table ( 0 , 0 ) = 1 ;
  for ( int degree = 1 ; degree < degree_limit ; degree++ )
  {
    for ( int x2 = 0 ; x2 <= degree ; x2++ )
      _fill_basis_table ( x2 , degree ) ;
  }
}

// routine for one iteration of Newton's method applied to
// a polynomial with ncoeff coefficients stored at pcoeff.
// the polynomial and it's first derivative are computed
// at x, then the quotient of function value and derivative
// is subtracted from x, yielding the result, which is
// closer to the zero we're looking for.

template < typename dtype >
void newton ( dtype & result , dtype x , int ncoeff , dtype * pcoeff )
{
  dtype * pk = pcoeff + ncoeff - 1 ; // point to last pcoefficient
  dtype power = 1 ;                  // powers of x
  int xd = 1 ;                       // ex derivative, n in (x^n)' = n(x^(n-1))
  dtype fx = 0 ;                     // f(x)
  dtype dfx = 0 ;                    // f'(x)
  
  // we work our way back to front, starting with x^0 and raising
  // the power with each iteration
  
  fx += power * *pk ;
  for ( ; ; )
  {
    --pk ;
    if ( pk == pcoeff )
      break ;
    dfx += power * xd * *pk ;
    xd++ ;
    power *= x ;
    fx += power * *pk ;
  }
  power *= x ;
  fx += power * *pk ;
  result = x - fx / dfx ;
}

// calculate_prefilter_poles relies on the table with basis function
// values having been filled with fill_basis_table(). It extracts
// the basis function values at even x2 (corresponding to whole x)
// as double precision values, takes these to be coefficients of
// a polynomial, and uses gsl and blas to obtain the roots of this
// polynomial. The real parts of those roots which are less than one
// are the raw values of the filter poles we're after.
// Currently I don't have an arbitrary-precision root finder, so I
// use the double precision one from GSL and submit the 'raw' result
// to polishing code in high precision arithmetic. This only works
// as far as the root finder can go, I found from degree 46 onwards
// it fails to find some roots, so that's how far I take it for now.
// The poles are stored in a std::vector of mpf_class and returned.
// These values are precise to mpf_bits, since after using gsl's
// root finder in double precision, they are polished as mpf_class
// values with the newton method, so the very many postcomma digits
// we print out below are all proper significant digits.

void calculate_prefilter_poles ( std::vector<mpf_class> &res ,
                                 int degree )
{
  const int r = degree / 2 ;
  
  // we need storage for the coefficients of the polynomial
  // first in double precision to feed gsl's root finder,
  // then in mpf_class, to do the polishing
  
  double coeffs [ 2 * r + 1 ] ;
  mpf_class mpf_coeffs [ 2 * r + 1 ] ;
  
  // here is space for the roots we want to compute
  
  double roots [ 4 * r + 2 ] ;
  
  // we fetch the coefficients from the table of b-spline basis
  // function values at even x2, corresponding to whole x
  
  double * pcoeffs = coeffs ;
  mpf_class * pmpf_coeffs = mpf_coeffs ;
  
  for ( int x2 = -r * 2 ;
        x2 <= r * 2 ;
        x2 += 2 , pcoeffs++ , pmpf_coeffs++ )
  {
    *pmpf_coeffs = access_table ( degree , x2 ) ;
    *pcoeffs = access_table ( degree , x2 ) . get_d() ;
  }
  
  // we set up the environment gsl needs to find the roots
  
  gsl_poly_complex_workspace * w 
          = gsl_poly_complex_workspace_alloc ( 2 * r + 1 ) ;
          
  // now we call gsl's root finder
          
  gsl_poly_complex_solve ( coeffs , 2 * r + 1 , w , roots ) ;
  
  // and release it's workspace
  
  gsl_poly_complex_workspace_free ( w ) ;

  // we only look at the real parts of the roots, which are stored
  // interleaved real/imag. And we take them back to front, even though
  // it doesn't matter which end we start with - but conventionally
  // Pole[0] is the root with the largest absolute, so I stick with that.
  
  // I tried using eigen alternatively, but it found less roots
  // than gsl/blas, so I stick with the GNU code, but for the reference,
  // here is the code to use with eigen3 - for degrees up to 23 it was okay
  // when I last tested. TODO: maybe the problem is that I did not get
  // long double values to start with, which is tricky from GMP
//   {
//     using namespace Eigen ;
//     
//     Eigen::PolynomialSolver<long double, Eigen::Dynamic> solver;
//     Eigen::Matrix<long double,Dynamic, 1 >  coeff(2*r+1);
//     
//     for ( int i = 0 ; i < 2*r+1 ; i++ )
//     {
//       // this is stupid: can't get a long double for an mpq_type
//       // TODO take the long route via a string (sigh...)
//       coeff[i] = mpf_coeffs[i].get_d() ;
//     }
//     
//     solver.compute(coeff);
//     
//     const Eigen::PolynomialSolver<long double, Eigen::Dynamic>::RootsType & r
//       = solver.roots();
//     
//     for ( int i = r.rows() - 1 ; i >= 0 ; i-- )
//     {
//       if ( std::fabs ( r[i].real() ) < 1 )
//       {
//         std::cout << "eigen gets " << r[i].real() << std::endl ;
//       }
//     }
//   }

  for ( int i = 2 * r - 2 ; i >= 0 ; i -= 2 )
  {
    if ( std::abs ( roots[i] ) < 1.0 )
    {
      // fetch a double precision root from gsl's output
      // converting to high-precision
      
      mpf_class root ( roots[i] ) ;
      
      // for the polishing process, we need a few more
      // high-precision values
      
      mpf_class pa ,     // current value of the iteration
                pb ,     // previous value
                pdelta ; // delta we want to go below
      
      pa = root , pb = 0 ;
      pdelta = "1e-300" ;
      
      // while we're not yet below delta
      
      while ( ( pa - pb ) * ( pa - pb ) >= pdelta )
      {
        // assign current to previous
        
        pb = pa ;
        
        // polish current value
        
        newton ( pa , pa , 2*r+1 , mpf_coeffs ) ;
      }
      
      // polishing iteration terminated, we're content and store the value
      
      root = pa ;
      res.push_back ( root ) ;
    }
  }
}

// forward recursive filter, used for testing

void forward ( mpf_class * x , int size , mpf_class pole )
{
  mpf_class X = 0 ;

  for ( int n = 1 ; n < size ; n++ )
  {
    x[n] = X = x[n] + pole * X ;
  }
}

// backward recursive filter, used for testing

void backward ( mpf_class * x , int size , mpf_class pole )
{
  mpf_class X = 0 ;
  
  for ( int n = size - 2 ; n >= 0 ; n-- )
  {
    x[n] = X = pole * ( X - x[n] );
  }
}

// to test the poles, we place the reconstruction kernel - a unit-spaced
// sampling of the basis function at integral x - into the center of a large
// buffer which is otherwise filled with zeroes. Then we apply all poles
// in sequence with a forward and a backward run of the prefilter.
// Afterwards, we expect to find a unit pulse in the center of the buffer.
// Since we use very high precision (see mpf_bits) we can set a conservative
// limit for the maximum error, here I use 1e-150. If this test throws up
// warnings, there might be something wrong with the code.

void test ( int size ,
            const std::vector<mpf_class> & poles ,
            int degree )
{
  mpf_class buffer [ size ] ;
  for ( int k = 0 ; k < size ; k++ )
  {
    buffer[k] = 0 ;
  }
  
  // calculate overall gain of the filter
  
  mpf_class lambda = 1 ;

  for ( int k = 0 ; k < poles.size() ; k++ )

    lambda = lambda * ( 1 - poles[k] ) * ( 1 - 1 / poles[k] ) ;

  int center = size / 2 ;
  
  // put basis function samples into the buffer, applying gain
  
  for ( int x2 = 0 ; x2 <= degree ; x2 += 2 )
      buffer [ center - x2/2 ]
    = buffer [ center + x2/2 ]
    = lambda * access_table ( degree , x2 ) ;
  
  // filter

  for ( int k = 0 ; k < poles.size() ; k++ )
  {
    forward ( buffer , size , poles[k] ) ;
    backward ( buffer , size , poles[k] ) ;
  }

  // test

  mpf_class error = 0 ;
  mpf_class max_error = 0 ;

  for ( int x2 = 0 ; x2 <= degree ; x2 += 2 )
  {
    error = abs ( buffer [ center - x2/2 ] - ( x2 == 0 ? 1 : 0 ) ) ;
    if ( error > max_error )
      max_error = error ;
    error = abs ( buffer [ center + x2/2 ] - ( x2 == 0 ? 1 : 0 ) ) ;
    if ( error > max_error )
      max_error = error ;
  }
  mpf_class limit = 1e-150 ;
  if ( max_error > limit )
    std::cerr << "WARNING: degree " << degree
              << " error exceeds limit 1e-150: "
              << max_error << std::endl ;
}

// we want to do the test in long double as well, to see how much
// downcasting the data to long double affects precision:

vspline::xlf_type get_xlf ( const mpf_class & _x )
{
  long exp ;
  std::string sign ( "+" ) ;
  mpf_class x ( _x ) ;
  
  if ( x < 0 )
  {
    x = -x ;
    sign = std::string ( "-" ) ;
  }
  
  auto str = x.get_str ( exp , 10 , 64 ) ;

  std::string res =   sign + std::string ( "." ) + str
                    + std::string ( "e" ) + std::to_string ( exp )
                    + std::string ( "l" ) ;
  
  vspline::xlf_type ld = std::stold ( res ) ;
  
  return ld ;
} ;

// forward recursive filter, used for testing

void ldforward ( vspline::xlf_type * x , int size , vspline::xlf_type pole )
{
  vspline::xlf_type X = 0 ;

  for ( int n = 1 ; n < size ; n++ )
  {
    x[n] = X = x[n] + pole * X ;
  }
}

// backward recursive filter, used for testing

void ldbackward ( vspline::xlf_type * x , int size , vspline::xlf_type pole )
{
  vspline::xlf_type X = 0 ;
  
  for ( int n = size - 2 ; n >= 0 ; n-- )
  {
    x[n] = X = pole * ( X - x[n] );
  }
}

// test with values cast down to long double

void ldtest ( int size ,
              const std::vector<mpf_class> & poles ,
              int degree )
{
  int nbpoles = degree / 2 ;
  
  vspline::xlf_type buffer [ size ] ;
  for ( int k = 0 ; k < size ; k++ )
  {
    buffer[k] = 0 ;
  }
  
  // calculate overall gain of the filter

  vspline::xlf_type lambda = 1 ;

  for ( int k = 0 ; k < nbpoles ; k++ )
  {
    auto ldp = get_xlf ( poles[k] ) ;
    lambda = lambda * ( 1 - ldp ) * ( 1 - 1 / ldp ) ;
  }
  
  int center = size / 2 ;
  
  // put basis function samples into the buffer, applying gain
  
  for ( int x2 = 0 ; x2 <= degree ; x2 += 2 )
  {
    vspline::xlf_type bf = get_xlf ( access_table ( degree , x2 ) ) ;
    
    buffer [ center - x2/2 ]
    = buffer [ center + x2/2 ]
    = lambda * bf ;
  }
  
  // filter

  for ( int k = 0 ; k < nbpoles ; k++ )
  {
    auto ldp = get_xlf ( poles[k] ) ;
    ldforward ( buffer , size , ldp ) ;
    ldbackward ( buffer , size , ldp ) ;
  }

  // test

  vspline::xlf_type error = 0 ;
  vspline::xlf_type max_error = 0 ;

  for ( int x2 = 0 ; x2 <= degree ; x2 += 2 )
  {
    error = std::abs ( buffer [ center - x2/2 ] - ( x2 == 0 ? 1 : 0 ) ) ;
    if ( error > max_error )
      max_error = error ;
    error = std::abs ( buffer [ center + x2/2 ] - ( x2 == 0 ? 1 : 0 ) ) ;
    if ( error > max_error )
      max_error = error ;
  }
  vspline::xlf_type limit = 1e-12 ;
  
  if ( max_error > limit )
    std::cerr << "WARNING: degree " << degree
              << " ld error exceeds limit 1e-12: "
              << max_error << std::endl ;
}

#include <random>

int main ( int argc , char * argv[] )
{
  mpf_set_default_prec ( mpf_bits ) ;
  
  bool print_basis_function = true ;
  bool print_prefilter_poles = true ;
  bool print_umbrella_structures = true ;

  std::cout << std::setprecision ( 64 ) ;
  
  fill_basis_table() ;
  mpf_class value ;
                 
  if ( print_basis_function )
  {
    for ( int degree = 0 ; degree < degree_limit ; degree++ )
    {
      std::cout << "const vspline::xlf_type K"
                << degree
                << "[] = {"
                << std::endl ;
      for ( int x2 = 0 ; x2 <= degree ; x2++ )
      {
        value = access_table ( degree , x2 ) ;
        std::cout << "  XLF("
                  << value
                  << ") ,"
                  << std::endl ;
        get_xlf ( value ) ;
      }
      std::cout << "} ;"
                << std::endl ;
    }
  }
  
  for ( int degree = 2 ; degree < degree_limit ; degree++ )
  {
    std::vector<mpf_class> res ;
    calculate_prefilter_poles ( res , degree ) ;
    
    if ( res.size() != degree / 2 )
    {
      std::cerr << "not enough poles for degree " << degree << std::endl ;
      continue ;
    }
    
    if ( print_prefilter_poles )
    {
      std::cout << "const vspline::xlf_type Poles_"
                << degree
                << "[] = {"
                << std::endl ;

      for ( int p = 0 ; p < degree / 2 ; p++ )
      {
        value = res[p] ;
        std::cout << "  XLF("
                  << value
                  << ") ,"
                  << std::endl ;
      }
      std::cout << "} ;"
                << std::endl ;
    }
    // test if the roots are good
    test ( 10000 , res , degree ) ;
    // optional, to see what happens when data are cast down to long double
    ldtest ( 10000 , res , degree ) ;
  }
  
  if ( print_umbrella_structures )
  {
    std::cout << std::noshowpos ;
    std::cout << "const vspline::xlf_type* const precomputed_poles[] = {"
              << std::endl ;
    std::cout << "  0, " << std::endl ;
    std::cout << "  0, " << std::endl ;
    for ( int i = 2 ; i < degree_limit ; i++ )
      std::cout << "  Poles_" << i << ", " << std::endl ;
    std::cout << "} ;" << std::endl ;
    std::cout << "const vspline::xlf_type* const precomputed_basis_function_values[] = {"
              << std::endl ;
    for ( int i = 0 ; i < degree_limit ; i++ )
      std::cout << "  K" << i << ", " << std::endl ;
    std::cout << "} ;" << std::endl ;
  }
}