/**
file: interpol.cpp
date: 2018-03-28
author: Michal Krelina
purpose:
  method for interpolation of data from grid

  Please cite the following paper if you use the code with vector mesons wave function of total cross section grids:
  J.Cepila, M.Krelina, J.Nemchik, R.Pasechnik, Phys.Rev.DXX:XXXX,2018 [arXiv:XXXX.XXXX].
**/


//...basic system include....
#include "interpol.h"

namespace interpol {

  void iniFile(data & s, std::string file, int print)
  {

    s.file = file;

    //...open file
    std::ifstream f(s.file.c_str());
    if (!f.is_open()) { // if cannot open
      std::cerr << "ERROR: cannot open file: " << s.file << std::endl;
      exit(1);
    }

    

    // counting of nonempty lines
    std::string line;
    std::string firstLine;
    while (getline(f, line)) {
      if (s.N == 0)
        firstLine = line;
      if (!line.empty())
        s.N++;
    }

    // counting of cols from firstLine
    std::stringstream ss;
    ss.clear();
    ss << firstLine;
    double temp;
    int ncols = 0;
    while (ss >> temp)
    {
      ncols++;
    }
    if (ncols == 2) {
      s.xDim = 1;
      s.yDim = 1;
    }
    else if (ncols == 3) {
      s.xDim = 2;
      s.yDim = 1;
    }
    else {
      std::cerr << "ERROR: wrong format of the grid file. " << std::endl;
    }

    //...initialise some arrays arrays
    s.y = new double*[s.yDim]; // 2D
    for (int i = 0; i < s.yDim; ++i) {
      s.y[i] = new double[s.N];
    }

    s.xNum = new int[s.xDim]; // 1D
    s.xLimits = new double*[s.xDim]; // 2D
    for (int i = 0; i < s.xDim; ++i) {
      s.xLimits[i] = new double[2];
    }

    //...temporary arrays
    // will keep all x values
    double **t = new double*[s.xDim]; //2D
    for (int i = 0; i < s.xDim; ++i) {
      t[i] = new double[s.N];
    }
    // will keep y values
    double **ty = new double*[s.yDim]; //2D
    for (int i = 0; i < s.yDim; ++i) {
      ty[i] = new double[s.N];
    }
    // will keep position in array ty
    int **ti = new int*[s.N]; //2D
    for (int i = 0; i < s.N; ++i) {
      ti[i] = new int[s.xDim];
    }



    //...read data
    f.clear(); // retrun to the begining of the file
    f.seekg(0, std::ios::beg); // retrun to the begining of the file
    double* xx = new double[s.xDim];
    for (int i = 0; i < s.N; i++)
    {
      for (int j = 0; j < s.xDim; j++) {
        f >> t[j][i];
      }
      for (int j = 0; j < s.yDim; j++) {
        f >> ty[j][i];
      }
    }

    f.close();

    //...counting number of each x value
    int tmpx;
    int xTot = 0;
    for (int i = 0; i < s.xDim; i++)
    {
      tmpx = 0;
      double xx;
      for (int j = 0; j < s.N; j++) {
        if (j == 0) {
          xx = t[i][j];
          tmpx++;
        }
        else {
          if (t[i][j] > xx) {
            xx = t[i][j];
            tmpx++;
          }
        }
      }
      s.xNum[i] = tmpx;
      xTot += tmpx;
    }

    // this array keep number of xvalues from more cols
    // necessary for array position 
    int* ni = new int[s.xDim];
    for (int i = s.xDim - 1; i >= 0; i--) {
      ni[i] = 1;
      for (int j = s.xDim - 1; j > i; j--) {
        ni[i] *= s.xNum[j];
      }
    }

    // calculation of xND array of ty
    for (int i = 0; i < s.xDim; i++) {
      int n = 0;
      int m = 0;
      for (int j = 0; j < s.N; j++)
      {
        ti[j][i] = m;
        n++;
        if (n > ni[i] - 1) {
          m++;
          n = 0;
        }

        if (m >= s.xNum[i]) {
          m = 0;
          n = 0;
        }
      }
    }


    //...define x values array
    s.x = new double[xTot]; // 1D
    // fill the x values array
    // and find the limits
    int tx = 0;
    for (int i = 0; i < s.xDim; i++)
    {
      double xx;
      for (int j = 0; j < s.N; j++) {
        if (j == 0) {
          xx = t[i][j];
          s.x[tx] = xx;
          tx++;
          s.xLimits[i][0] = xx;
        }
        else {
          if (t[i][j] > xx) {
            xx = t[i][j];
            s.x[tx] = xx;
            tx++;
          }
        }
      }
      s.xLimits[i][1] = xx;
    }

    // next temporary array for changing position
    for (int i = 0; i < s.xDim; i++) {
      ni[i] = 1;
      for (int j = 0; j < i; j++) {
        ni[i] *= s.xNum[j];
      }
    }

    // reposition of values of y for right using with function get
    for (int i = 0; i < s.N; i++) {
      int pos = 0;
      for (int j = 0; j < s.xDim; j++) {
        pos += ti[i][j] * ni[j];
      }
      for (int j = 0; j < s.yDim; j++) {
        s.y[j][pos] = ty[j][i];
      }

    }


    // print info
    if (print == 1) {
      std::cout << "Data loaded:" << std::endl;
      std::cout << "..file: " << s.file << std::endl;
      std::cout << "..number of X values: " << s.xDim << std::endl;
      std::cout << "..number of Y values: " << s.yDim << std::endl;
      std::cout << "..number of lines: " << s.N << std::endl;
      std::cout << "..x values: (i, number, lower limit, upper limit)" << std::endl;
      for (int i = 0; i < s.xDim; i++)
      {
        std::cout << "  .." << i << "\t" << s.xNum[i] << "\t" << s.xLimits[i][0] << "\t" << s.xLimits[i][1] << std::endl;
      }
    }

    // data loaded flag
    s.isData = 1;

    //...clean temporary arrs
    delete[] t;
    delete[] ty;
    delete[] ti;
    delete[] ni;

  }

  double get(data & s, double* x)
  {
    double y[1];
    if (s.isData != 1) {
      std::cerr << "ERROR: no data loaded - need to call ini()" << std::endl;
      return 0.;
    }

    //... check x values and set min/max
    for (int i = 0; i < s.xDim; i++) {
      if (x[i] <= s.xLimits[i][0])
        x[i] = s.xLimits[i][0];
      else if (x[i] >= s.xLimits[i][1])
        x[i] = s.xLimits[i][1];
    }

    for (int i = 0; i < s.yDim; i++) {
      y[i] = 0.;
    }


    //...interpoluji
    //...THE INTERPOLATION ROUTINE (CERN LIBRARY ROUTINE E104)

    int* NCOMB = new int[s.xDim];
    double* D = new double[s.xDim];
    int* IENT = new int[s.xDim];
    int KD, M, J, JR, JA, JB;

    KD = 1;
    M = 1;
    JA = 1;

    for (int I = 1; I <= s.xDim; I++)
    {
      NCOMB[I - 1] = 1;
      JB = JA - 1 + s.xNum[I - 1];

      for (J = JA; J <= JB; J++) {
        if (x[I - 1] <= s.x[J - 1])
          goto g3;
        continue;
      }
      J = JB;

    g3:
      if (J != JA)
        goto g4;

      J = J + 1;

    g4:
      JR = J - 1;
      //std::cout << "D > " << s.x[J - 1] << " - " << x[I - 1] << "\t" << s.x[J - 1] << " - " << s.x[JR - 1] << std::endl;
      D[I - 1] = (s.x[J - 1] - x[I - 1]) / (s.x[J - 1] - s.x[JR - 1]);
      //std::cout << "DD > " << D[I - 1] << std::endl;
      IENT[I - 1] = J - JA;
      KD = KD + IENT[I - 1] * M;
      M = M * s.xNum[I - 1];

    g5:
      JA = JB + 1;

    }

    double FAC;
    int IADR, IFADR, IL;

  g10: 
    FAC = 1.;
    IADR = KD;
    IFADR = 1;
    for (int I = 1; I <= s.xDim; I++)
    {
      if (NCOMB[I - 1] == 0)
        goto g12;
      //std::cout << "1 > " << I << "\t" << FAC << std::endl;
      FAC = FAC*(1. - D[I - 1]);
      //std::cout << "2 > " << I << "\t" << FAC << "\t" << D[I - 1] << std::endl;
      goto g15;

    g12:
      FAC = FAC*D[I - 1];
      //std::cout << "3 > " << I << "\t" << FAC << std::endl;
      IADR = IADR - IFADR;

    g15:
      IFADR = IFADR*s.xNum[I - 1];
    }

    for (int k = 0; k < s.yDim; k++) {
      //std::cout << "int: "<< k << "\t" << y[k] << "\t" << FAC << "\t" << s.y[k][IADR - 1] << "\t" << IADR - 1 << std::endl;
      y[k] = y[k] + FAC*s.y[k][IADR - 1];
    }
    IL = s.xDim;

  g40: // GO TO 40
    if (NCOMB[IL - 1] == 0)
      goto g80;
    NCOMB[IL - 1] = 0;

    if (IL == s.xDim)
      goto g10;
    IL = IL + 1;

    for (int K = IL; K <= s.xDim; K++) {
      NCOMB[K - 1] = 1;
    }
    goto g10;

  g80: // GO TO 80
    IL = IL - 1;

    if (IL != 0)
      goto g40;
  
    // cleaning...
    delete[] NCOMB;
    delete[] D;
    delete[] IENT;
    return y[0];
  }

  void del(data & s)
  {
    if (s.isData != 1) {
      std::cerr << "ERROR: no data loaded - need to call ini()" << std::endl;
      return;
    }
    // clean all data
    delete[] s.x;
    delete[] s.y;
    delete[] s.xLimits;
    delete[] s.xNum;
    s.isData = 0;
  }

}