#include "Function.hpp"

#include <cassert>
#include <QFile>
#include <QTextStream>
#include <stdio.h>
#include <QString>
#include <QStringList>

#include <Information.hpp>

#if defined(__APPLE__) || defined(linux)
#  include <sys/file.h>
#endif

#define REPORT_ERROR(filename, line_nb, error) \
  error_occured = true; \
  err << "Error in file " << filename << " on line " << line_nb << error << endl

static QTextStream err(stderr);

namespace mgx
{

  /** @brief Initialization function that should be called from every constructor
  */
  void Function::init()
  {
    samples = 100000;
    cache_valid = false;
    error_occured = false;
  }

  /** @brief Default constructor. */
  Function::Function() : pts(), max(), min(), scaling_x(1), scaling_y(1), shift_x(0), shift_y(0)
  {
    init();
  }

  /** @brief Constructor with initialisation from file.
    @bug   If the specified file contains an incomplete or improper
    specification, the result of construction will be
    unknown.  Most probably, the object will unuseable;
    however, this is not currently reported.
    @param filename The function file to load.
    */
  Function::Function(const QString &file) : filename(file), pts(), max(), min(), 
    scaling_x(1), scaling_y(1), shift_x(0), shift_y(0)
  {
    init();
    reRead();
  }

  /** @brief Copy constructor. */
  Function::Function(const Function& copy) : filename(copy.filename), pts(copy.pts), max(copy.max), 
    min(copy.min), samples(copy.samples), scaling_x(copy.scaling_x), scaling_y(copy.scaling_y), 
    shift_x(copy.shift_x), shift_y(copy.shift_y), cache_valid(copy.cache_valid), cache(copy.cache),
    error_occured(copy.error_occured) {}

  Function& Function::operator=(const Function& f) 
  {
    filename = f.filename;
    pts = f.pts;
    max = f.max;
    min = f.min;
    samples = f.samples;
    cache_valid = f.cache_valid;
    cache = f.cache;
    error_occured = f.error_occured;
    scaling_x = f.scaling_x;
    scaling_y = f.scaling_y;
    shift_x = f.shift_x;
    shift_y = f.shift_y;
    return *this;
  }

  bool Function::error() { return error_occured; }

  bool Function::setFile(const QString &file)
  {
    init();
    filename = file;
    if(filename.size() > 0)
      reRead();

    return true;
  }

  void Function::reRead() 
  {
    QFile f(filename);
    error_occured = false;
    if(!f.open(QIODevice::ReadOnly))
    {
      error_occured = true;
      throw(QString("Error opening file: %1, %2").arg(filename).arg(f.errorString()));
    }
#if defined(__APPLE__) || defined(linux)
    flock(f.handle(), LOCK_SH);
#endif
    QTextStream ts(&f);
    int ptid = 0, line_nb = 0;

    max = Vector<2,double>(-HUGE_VAL,-HUGE_VAL);
    min = Vector<2,double>(HUGE_VAL,HUGE_VAL);

    while(!ts.atEnd())
    {
      QString line = ts.readLine().trimmed();
      line_nb++;
      if(line.startsWith("fver"))
      {
        QStringList fields = line.split(" ");
        if(fields.size() != 3)
        {
          REPORT_ERROR(filename, line_nb, 
                "cannot find major and minor in function version (fver)");
          return;
        }
        bool ok;
        int major = fields[1].toInt(&ok);
        if(!ok)
        {
          REPORT_ERROR(filename, line_nb, 
                "major version number (" << fields[1] << ") is not a valid integer");
          return;
        }
        int minor = fields[2].toInt(&ok);
        if(!ok)
        {
          REPORT_ERROR(filename, line_nb, "minor version number (" << fields[2] 
                << ") is not a valid integer");
          return;
        }
        if(major != 1 or minor != 1)
        {
          err << "Warning, this version usually deals with fver 1 1. "
                 "Processing will continue but may fail" << endl;
        }
      }
      else if(line.startsWith("points:"))
      {
        QStringList fields = line.split(":");
        if(fields.size() != 2)
        {
          REPORT_ERROR(filename, line_nb, 
            "Invalid point number specification, there is more than one ':'");
          return;
        }
        bool ok;
        unsigned int nbpts = fields[1].toUInt(&ok);
        if(!ok)
        {
          REPORT_ERROR(filename, line_nb, "Number of points (" << fields[1] 
                                                       << ") is not a valid integer");
          return;
        }
        pts.resize(nbpts);
      }
      else if(!line.contains(":")) // specification of a points
      {
        if(ptid >= (int)pts.size())
        {
          REPORT_ERROR(filename, line_nb, "There are more points than specified");
          return;
        }
        QStringList coords = line.split(" ");
        if(coords.size() != 2)
        {
          REPORT_ERROR(filename, line_nb, "There are " << coords.size() 
                                                    << " coordinates instead of 2");
          return;
        }
        double x, y;
        bool ok;
        x = coords[0].toDouble(&ok);
        if(!ok)
        {
          REPORT_ERROR(filename, line_nb, coords[0] << "' is not a valid double");
          return;
        }
        y = coords[1].toDouble(&ok);
        if(!ok)
        {
          REPORT_ERROR(filename, line_nb, coords[1] << "' is not a valid double");
          return;
        }
        if(x < min.x()) min.x() = x;
        if(y < min.y()) min.y() = y;
        if(x > max.x()) max.x() = x;
        if(y > max.y()) max.y() = y;
        Vector<2,double> p(x, y);
        pts[ptid++] = p;
      }
      // Otherwise, just ignore the line
    }
    if(ptid != (int)pts.size())
    {
      error_occured = true;
      err << "Error in file " << filename
        << ", there was " << pts.size() << " points announced, but only " << ptid 
                                                    << " points declared" << endl;
      return;
    }

    cache_valid = false;
#if defined(__APPLE__) || defined(linux)
    flock(f.handle(), LOCK_UN);
#endif
    f.close();
    //  ifstream in(filename.c_str());
    //  string buffer;

    //  if (!in || !in.good() || in.eof()) {
    //    error_occured = true;
    //    std::cerr << "Function - Cannot open file '" << filename << "'\n";
    //    return;
    //  }
    //  error_occured = false;

    //  in >> ws >> buffer;

    //  if (buffer == string("range:")) {
    //    // old version

    //    double rmin, rmax;
    //    in >> rmin >> rmax;
    //  }
    //  else if (buffer == string("fver")) {
    //    int major, minor;
    //    in >> major >> minor >> ws;

    //    if (major == 1 && minor == 1) {
    //      getline(in, buffer); // name
    //      getline(in, buffer); // samples
    //      getline(in, buffer); // flip
    //    }
    //  }

    //  unsigned int num;
    //  in >> buffer >> num;

    //  pts.reserve(num);

    //  bool first = true;
    //  for (unsigned int i = 0; i < num; i++) {
    //    double x, y;
    //    in >> x >> y;

    //    if (first) {
    //      max.set(x, y);
    //      min.set(x, y);
    //      first = false;
    //    }
    //    else {
    //      if (x > max.x()) max.x(x);
    //      if (y > max.y()) max.y(y);
    //      if (x < min.x()) min.x(x);
    //      if (y < min.y()) min.y(y);
    //    }

    //    Vector<2,double> p(x, y);
    //    pts.push_back(p);
    //  }

    //  cache_valid = false;
  }

  /** @brief sets the number of samples to be used
    @param n The number of samples to use

    This sets the new number of samples and invalidates the cache.
    */
  void Function::setSamples (size_t n)
  {
    samples = (unsigned int)n;
    cache_valid = false;
  }

  /** @brief Function operator.
    @param x The position in the function.

    This returns the y-value of the function for a given x-value.  If x is outside the domain of
    the function, the value of the function start or end point is returned. Dynamic caching is
    used to speed things up. Basically, the function is only evaluated at a number of places
    (specified in 'samples'), and otherwise the values are linearly interpolated inbetween. The
    real values are obtained by calling getVal ().
    */
  double Function::operator()(double x) {
    x += shift_x;
    x *= scaling_x;
    if (x <= min.x()) return scaling_y*pts[0].y() + shift_y;
    if (x >= max.x()) return scaling_y*pts[pts.size() - 1].y() + shift_y;

    // check for cache
    if (! cache_valid) {
      cache_valid = true;
      cache.resize (samples+1);
      for (size_t i = 0 ; i < samples ; i ++)
        cache[i].valid = false;
    }
    // lower index
    double dx = (max.x() - min.x()) / samples;
    size_t lo = (size_t) ((x - min.x()) / dx);
    size_t hi = lo + 1;
    // make sure we have both lo and hi
    double lox = min.x() + dx * lo;
    double hix = lox + dx;
    if (! cache [lo].valid) {
      cache [lo].valid = true;
      cache [lo].val = getVal (lox);
    }
    if (! cache [hi].valid) {
      cache [hi].valid = true;
      cache [hi].val = getVal (hix);
    }
    // linearly interpolate
    double r = (x - lox) / dx;
    double ret = (1-r) * cache[lo].val + r * cache[hi].val;

    return scaling_y*ret + shift_y;
  }

  /**
   * Make sure the x axis spans from 0 to 1
   * \param shift If true, it will also change the shift of the axis
   */
  void Function::normalizeX(bool shift)
  {
    if(shift)
    {
      shift_x = min.x();
      scaling_x = max.x() - min.x();
    }
    else
    {
      scaling_x = max.x();
    }
  }

  /**
   * Make sure the y axis spans from 0 to 1
   * \param shift If true, it will also change the shift of the axis
   */
  void Function::normalizeY(bool shift)
  {
    double ymin = HUGE_VAL, ymax = -HUGE_VAL;
    double dx = (max.x() - min.x()) / samples;
    double x0 = min.x();
    if (! cache_valid) {
      cache_valid = true;
      cache.resize (samples+1);
      for (size_t i = 0 ; i < samples ; i ++)
        cache[i].valid = false;
    }
    for(size_t i = 0 ; i < samples+1 ; ++i)
    {
      double y;
      if(cache[i].valid)
      {
        y = cache[i].val;
      }
      else
      {
        y = getVal(x0+dx*i);
        cache[i].valid = true;
        cache[i].val = y;
      }
      if(y < ymin) ymin = y;
      if(y > ymax) ymax = y;
    }
    if(shift)
    {
      shift_y = -ymin;
      if(ymax > ymin)
        scaling_y = 1/(ymax-ymin);
      else
        scaling_y = 1;
    }
    else
    {
      if(ymax != 0)
      {
        scaling_y = 1/ymax;
      }
      else
      {
        scaling_y = 1;
      }
    }
  }

  /** @brief Function operator.
    @param x The position in the function.

    This returns the y-value of the function for a given x-value. The parameter x is assumed to be
    between minx and maxx.
    */
  double Function::getVal (double x) const {
    const double tofind = x;
    const double MaxError = 0.00001; // TBD: this should be user controlled
    double low = 0.0;
    double high = 1.0;
    double check = 0.5;
    Vector<2,double> tst;
    int counter = 0;
    do {
      check = (low + high) / 2.0;
      tst = P(check);
      if (tst.x() > tofind)
        high = check;
      else 
        low = check;
      counter++;
    } while (fabs(tst.x() - tofind) > MaxError);
    return tst.y();
  }

  /** @brief P basis function */
  Vector<2,double> Function::P(double x) const {
    const int n = (unsigned int)pts.size() - 1;
    const int t = 4;
    double u = x * (double(n - t) + 2.0);

    Vector<2,double> sum;

    for (int k = 0; k <= n; k++) {
      double coeff = N(k, t, u);
      sum += pts[k] * coeff;
    }

    return sum;
  }

  /** @brief A basis function. */
  double Function::N(int k, int t, double u) const {
    double res = 0.0;
    if (1==t)
      res = Nk1(k, u);
    else
      res = Nkt(k, t, u);
    return res;
  }


  /** @brief A basis function. */
  double Function::Nk1(int k, double u) const {
    if (Uk(k) <= u)
    {
      if (u < Uk(k + 1))
        return 1.0;
    }
    return 0.0;
  }

  /** @brief A basis function. */
  double Function::Nkt(int k, int t, double u) const {
    double sum = 0.0;
    int div = Uk(k + t - 1) - Uk(k);
    if (0 != div)
      sum = (u - Uk(k)) / div * N(k, t - 1, u);

    div = Uk(k + t) - Uk(k + 1);
    if (0 != div)
      sum += (Uk(k + t) - u) / div * N(k + 1, t - 1, u);

    return sum;
  }

  /** @brief A basis function. */
  int Function::Uk(int j) const {
    const int n = (unsigned int)pts.size() - 1;
    const int t = 4;
    if (j < t)
      return 0;
    if (j > n)
      return n - t + 2;
    return j - t + 1;
  }

  /** @brief Return the maximum x and y values. */
  const Vector<2,double>& Function::getMax() const {
    return max;
  }

  /** @brief Return the minimum x and y values. */
  const Vector<2,double>& Function::getMin() const {
    return min;
  }

}
