s2_ftrans.hpp

/* s2_ftrans.hpp */
/* $Id: s2_ftrans.hpp 21267 2011-11-23 15:19:39Z guser1 $ */
#ifndef __s2_ftrans_hpp
#define __s2_ftrans_hpp

#include <stdio.h>
#include <vector>
#include <math.h>

#ifdef GE_BUILD
#include <least_squares_2d.hpp>
#else
#include <mth/least_squares_2d.hpp>
#endif

namespace FTrans{
  
  typedef double arg_type;
  
  class BaseFunc {
    
  public:
    
    virtual ~BaseFunc() {};
    
    arg_type operator ()(arg_type x) { return calc(x); };
    
  protected:
    
    virtual arg_type calc(arg_type x) const = 0;
  };
  
  class InverseFunc : public BaseFunc {
  protected:
    
    virtual arg_type calc(arg_type x) const;
    
  };
  
  template<typename TR> 
  void ftrans(std::vector<arg_type> &func, arg_type &xmn, arg_type &xmx)
  {
#define FTRANS_REAL_EPS 1.e-6
    int debug = 0;
    TR trans;
    
    /* Recalculate argument limits */
    arg_type zmn = trans(xmn);
    arg_type zmx = trans(xmx);
    if (debug>0)
      printf("INFO: ftrans: new arg limits %lg->%lg\n", zmn, zmx);
    
    std::vector<arg_type> args(func.size());
    arg_type dx = (xmx-xmn) / (func.size()-1.);
    
    /* Recalculate arguments */
    arg_type x_cur = xmn;
    for(int i=0; i<func.size(); ++i) {
      args[i] = trans(x_cur);
      if (debug>0)
    printf("INFO: ftrans: arg[%ld] = %lg -> %lg\n", i, x_cur, args[i]);
      x_cur += dx;
    }
    if (debug>0)
      printf("INFO: ftrans: func recalculated...\n");
    
    std::vector<arg_type> func_new;
    
#if 0
    /* Save function on new irregular grid */
    FILE* f = fopen("ftrans_res1_f.dat", "w");
    for(int i=0; i<func.size(); ++i)
      fprintf(f,"%g\t%g\n", args[i], func[i]);
    fclose(f);
#endif
    
    std::size_t i_max = func.size()-1;
    
    /* Calculate regular argument grid step */
    arg_type dz = (zmx-zmn) / (func.size()-1.);
#if 1
    /* Find the minimal step in new irregular argument grid, and use it as the step of new regular grid. */
    for(int i=1; i<func.size(); ++i) {
      if (fabs(dz) > fabs(args[i]-args[i-1]))
    dz = args[i]-args[i-1];
    }
    i_max = int(fabs((zmx-zmn)/dz));
#endif
    if (debug>0)
      printf("INFO: ftrans: dz=%lg...\n", dz);
    arg_type z_cur = zmn+dz;
    
    /* Add first function value (immovable) */
    func_new.push_back(func[0]);
    
    for(std::size_t i=1; i<i_max; ++i) {
      std::size_t j=0;
      bool bingo = false;
      
      /* Search two values of argument (irregular grid) on both sides from the current argument value (regular grid) */
      for(; j<func.size()-1; ++j) {
    //double tmp = double(args[j]-z_cur);
    if (fabs(args[j]-z_cur) < FTRANS_REAL_EPS) { bingo = true; break; }
    if (args[j]>z_cur && args[j+1]<z_cur) break;
    if (args[j]<z_cur && args[j+1]>z_cur) break;
      }
      if (j == func.size()-1) {
    printf("ERROR: arg=%lg was not found...\n", z_cur);
    z_cur+=dz;
    continue;
      }
      else
    if (debug>0)
      printf("INFO: ftrans: arg=%lg limits found, j0=%zu...\n", z_cur, j);
      
      if (bingo) { /* The current argument value is equal to some value of argument of irregular grid */
    func_new.push_back(func[j]);
    z_cur+=dz;
    continue;
      }
      
      /* Form arrays of points (irregular grid argument values) in the neighbourhood of current argument value.
     Mean goal - to form the polynom and calculate function value with respect to current argument value.
     Number of selected points = degree of selected polynom + 1. Default polynom - parabolic. */
      unsigned int p_num = 3;
      std::vector<double> xp(p_num); //xp[0] = args[j]; xp[1] = args[j+1]; xp[2] = (j<(func.size()-2)) ? args[j+2] : args[j-1];
      std::vector<double> yp(p_num); //yp[0] = func[j]; yp[1] = func[j+1]; yp[2] = (j<(func.size()-2)) ? func[j+2] : func[j-1];
      
      if ((j+p_num-1) >= func.size())
    j = func.size()-p_num;
      
      for (std::size_t i=0; i<p_num; ++i) {
    xp[i] = args[j];
    yp[i] = func[j];
    if (debug>0)
      printf("INFO: ftrans: arg=%lg       func[%zu] = %lg\n", z_cur, j, yp[i]);
    ++j;
      }
      
      unsigned int deg = p_num-1;
      //deg = 2;
      LeastSquares2D::PolyCoeffs coeffs(deg+1);
      bool ok = LeastSquares2D::getPoly(xp, yp, coeffs);
      if (!ok) {
    printf("ERROR: ftrans: arg=%lg : failed to form poly...", z_cur);
      }
      else {
    double res = LeastSquares2D::calcPoly(coeffs, (double)z_cur);
    func_new.push_back(res);
      
    if (debug>0)
      printf("INFO: ftrans: arg=%lg process complete... val=%lg\n", z_cur, res);
      }
      
      z_cur+=dz;
    } // end of the cycle along argument range
    
    
    /* Add last function value (immovable) */
    func_new.push_back(func[func.size()-1]);
    if (debug>0)
      printf("INFO: ftrans: new func formed...\n");
    
    /* Clear initial function data */
    func.clear();
    
    /* Form result */
    for(std::size_t i=0; i<func_new.size(); ++i) {
      if (zmn>zmx)
    func.push_back(func_new[func_new.size()-i-1]);
      else
    func.push_back(func_new[i]);
      if (debug>0)
    printf("INFO: ftrans: new func [%zu] = %lg\n", i, func[i]);
    }
    
    /* Define correct argument limits sequence */
    if (zmn>zmx) {
      arg_type tmp = zmn;
      zmn = zmx;
      zmx = tmp;
    }
    
    xmn = zmn;
    xmx = zmx;
    //printf("INFO: ftrans: xmn=%lg xmx=%lg\n", xmn, xmx);
    if (debug>0)
      printf("INFO: ftrans: complete...\n");
#undef FTRANS_REAL_EPS
  };
      
};


#endif /* s2_ftrans.hpp */