#include "store.h"
  #include "model/cubic_model.h"
  #include "model/exponential_model.h"
  #include "model/linear_model.h"
  #include "model/logarithmic_model.h"
  #include "model/logistic_model.h"
  #include "model/power_model.h"
  #include "model/quadratic_model.h"
  #include "model/quartic_model.h"
  #include "model/trigonometric_model.h"
  #include "apps/apps_container.h"
  #include <poincare/preferences.h>
  #include <assert.h>
  #include <float.h>
  #include <cmath>
  #include <string.h>
  
  using namespace Shared;
  
  namespace Regression {
  
  static inline float max(float x, float y) { return (x>y ? x : y); }
  static inline float min(float x, float y) { return (x<y ? x : y); }
  
  static_assert(Model::k_numberOfModels == 9, "Number of models changed, Regression::Store() needs to adapt");
  static_assert(Store::k_numberOfSeries == 3, "Number of series changed, Regression::Store() needs to adapt (m_seriesChecksum)");
  
  Store::Store() :
    InteractiveCurveViewRange(nullptr),
    DoublePairStore(),
    m_seriesChecksum{0, 0, 0},
    m_angleUnit(Poincare::Expression::AngleUnit::Degree)
  {
    for (int i = 0; i < k_numberOfSeries; i++) {
      m_regressionTypes[i] = Model::Type::Linear;
      m_regressionChanged[i] = false;
    }
    m_regressionModels[0] = new LinearModel();
    m_regressionModels[1] = new QuadraticModel();
    m_regressionModels[2] = new CubicModel();
    m_regressionModels[3] = new QuarticModel();
    m_regressionModels[4] = new LogarithmicModel();
    m_regressionModels[5] = new ExponentialModel();
    m_regressionModels[6] = new PowerModel();
    m_regressionModels[7] = new TrigonometricModel();
    m_regressionModels[8] = new LogisticModel();
  }
  
  Store::~Store() {
    for (int i = 0; i < Model::k_numberOfModels; i++) {
      delete m_regressionModels[i];
    }
  }
  
  /* Regressions */
  void Store::setSeriesRegressionType(int series, Model::Type type) {
    assert(series >= 0 && series < k_numberOfSeries);
    if (m_regressionTypes[series] != type) {
      m_regressionTypes[series] = type;
      m_regressionChanged[series] = true;
    }
  }
  
  int Store::closestVerticalRegression(int direction, double x, double y, int currentRegressionSeries, Poincare::Context * globalContext) {
    int regressionSeries = -1;
    float closestDistance = INFINITY;
    /* The conditions to test on all the regressions are in this order:
     * - the current regression is not the current regression
     * - the next regression point should be within the window abscissa bounds
     * - it is the closest one in abscissa to x
     * - it is above y if direction > 0 and below otherwise */
    for (int series = 0; series < k_numberOfSeries; series ++) {
      if (!seriesIsEmpty(series) && series != currentRegressionSeries) {
        double regressionY = yValueForXValue(series, x, globalContext);
        if ((m_yMin <= regressionY && regressionY <= m_yMax)
            && (std::fabs(regressionY - y) < closestDistance)
            && (regressionY - y > 0) == (direction > 0)) {
          closestDistance = std::fabs(regressionY - y);
          regressionSeries = series;
        }
      }
    }
    return regressionSeries;
  }
  
  /* Dots */
  
  int Store::closestVerticalDot(int direction, double x, double y, int currentSeries, int currentDot, int * nextSeries, Poincare::Context * globalContext) {
    double nextX = INFINITY;
    double nextY = INFINITY;
    int selectedDot = -1;
    /* The conditions to test on all dots are in this order:
     * - if the currentDot is valid, the next series should not be the current series
     * - the next dot should not be the current dot
     * - the next dot should be within the window abscissa bounds
     * - the next dot is the closest one in abscissa to x
     * - the next dot is above the regression curve if direction == 1 and below
     * otherwise
     * - the next dot is above/under y
     * - if the current dot is valid, do not select a dot of the same series */
    for (int series = 0; series < k_numberOfSeries; series ++) {
      if (!seriesIsEmpty(series) && (currentDot < 0 || currentSeries != series)) {
        for (int index = 0; index < numberOfPairsOfSeries(series); index++) {
          if ((currentSeries != series) || (index != currentDot)) {
            double currentDataX = m_data[series][0][index];
            double currentDataY = m_data[series][1][index];
            if ((m_xMin <= currentDataX && currentDataX <= m_xMax) &&
                (std::fabs(currentDataX - x) <= std::fabs(nextX - x)) &&
                ((currentDataY - yValueForXValue(currentSeries, currentDataX, globalContext) >= 0) == (direction > 0)) &&
                ((currentDataY > y) == (direction > 0))) {
              // Handle edge case: if 2 dots have the same abscissa but different ordinates
              if (nextX != currentDataX || ((nextY - currentDataY >= 0) == (direction > 0))) {
                nextX = currentDataX;
                nextY = currentDataY;
                selectedDot = index;
                *nextSeries = series;
              }
            }
          }
        }
        // Compare with the mean dot
        if ((currentSeries != series) || (numberOfPairsOfSeries(series) != currentDot)) {
          double meanX = meanOfColumn(series, 0);
          double meanY = meanOfColumn(series, 1);
          if (m_xMin <= meanX && meanX <= m_xMax &&
              (std::fabs(meanX - x) <= std::fabs(nextX - x)) &&
              ((meanY - yValueForXValue(currentSeries, meanX, globalContext) >= 0) == (direction > 0)) &&
              ((meanY > y) == (direction > 0))) {
            if (nextX != meanX || ((nextY - meanY >= 0) == (direction > 0))) {
              selectedDot = numberOfPairsOfSeries(series);
              *nextSeries = series;
            }
          }
        }
      }
    }
    return selectedDot;
  }
  
  int Store::nextDot(int series, int direction, int dot) {
    float nextX = INFINITY;
    int selectedDot = -1;
    double meanX = meanOfColumn(series, 0);
    float x = meanX;
    if (dot >= 0 && dot < numberOfPairsOfSeries(series)) {
      x = get(series, 0, dot);
    }
    /* We have to scan the Store in opposite ways for the 2 directions to ensure to
     * select all dots (even with equal abscissa) */
    if (direction > 0) {
      for (int index = 0; index < numberOfPairsOfSeries(series); index++) {
        /* The conditions to test are in this order:
         * - the next dot is the closest one in abscissa to x
         * - the next dot is not the same as the selected one
         * - the next dot is at the right of the selected one */
        if (std::fabs(m_data[series][0][index] - x) < std::fabs(nextX - x) &&
            (index != dot) &&
            (m_data[series][0][index] >= x)) {
          // Handle edge case: 2 dots have same abscissa
          if (m_data[series][0][index] != x || (index > dot)) {
            nextX = m_data[series][0][index];
            selectedDot = index;
          }
        }
      }
      // Compare with the mean dot
      if (std::fabs(meanX - x) < std::fabs(nextX - x) &&
          (numberOfPairsOfSeries(series) != dot) &&
          (meanX >= x)) {
        if (meanX != x || (numberOfPairsOfSeries(series) > dot)) {
          selectedDot = numberOfPairsOfSeries(series);
        }
      }
    } else {
      // Compare with the mean dot
      if (std::fabs(meanX - x) < std::fabs(nextX - x) &&
          (numberOfPairsOfSeries(series) != dot) &&
          (meanX <= x)) {
        if ((meanX != x) || (numberOfPairsOfSeries(series) < dot)) {
          nextX = meanX;
          selectedDot = numberOfPairsOfSeries(series);
        }
      }
      for (int index = numberOfPairsOfSeries(series)-1; index >= 0; index--) {
        if (std::fabs(m_data[series][0][index] - x) < std::fabs(nextX - x) &&
            (index != dot) &&
            (m_data[series][0][index] <= x)) {
          // Handle edge case: 2 dots have same abscissa
          if (m_data[series][0][index] != x || (index < dot)) {
            nextX = m_data[series][0][index];
            selectedDot = index;
          }
        }
      }
    }
    return selectedDot;
  }
  
  /* Window */
  
  void Store::setDefault() {
    float minX = FLT_MAX;
    float maxX = -FLT_MAX;
    for (int series = 0; series < k_numberOfSeries; series++) {
      if (!seriesIsEmpty(series)) {
        minX = min(minX, minValueOfColumn(series, 0));
        maxX = max(maxX, maxValueOfColumn(series, 0));
      }
    }
    float range = maxX - minX;
    setXMin(minX - k_displayHorizontalMarginRatio*range);
    setXMax(maxX + k_displayHorizontalMarginRatio*range);
    setYAuto(true);
  }
  
  /* Series */
  
  bool Store::seriesIsEmpty(int series) const {
    return numberOfPairsOfSeries(series) < 2;
  }
  
  /* Calculations */
  
  double * Store::coefficientsForSeries(int series, Poincare::Context * globalContext) {
    assert(series >= 0 && series <= k_numberOfSeries);
    assert(!seriesIsEmpty(series));
    uint32_t storeChecksumSeries = storeChecksumForSeries(series);
    Poincare::Expression::AngleUnit currentAngleUnit = Poincare::Preferences::sharedPreferences()->angleUnit();
    if (m_angleUnit != currentAngleUnit) {
      m_angleUnit = currentAngleUnit;
      for (int i = 0; i < k_numberOfSeries; i++) {
        if (m_regressionTypes[i] == Model::Type::Trigonometric) {
          m_regressionChanged[i] = true;
        }
      }
    }
    if (m_regressionChanged[series] || (m_seriesChecksum[series] != storeChecksumSeries)) {
      Model * seriesModel = modelForSeries(series);
      seriesModel->fit(this, series, m_regressionCoefficients[series], globalContext);
      m_regressionChanged[series] = false;
      m_seriesChecksum[series] = storeChecksumSeries;
    }
    return m_regressionCoefficients[series];
  }
  
  double Store::doubleCastedNumberOfPairsOfSeries(int series) const {
    return DoublePairStore::numberOfPairsOfSeries(series);
  }
  
  float Store::maxValueOfColumn(int series, int i) const {
    float maxColumn = -FLT_MAX;
    for (int k = 0; k < numberOfPairsOfSeries(series); k++) {
      maxColumn = max(maxColumn, m_data[series][i][k]);
    }
    return maxColumn;
  }
  
  float Store::minValueOfColumn(int series, int i) const {
    float minColumn = FLT_MAX;
    for (int k = 0; k < numberOfPairsOfSeries(series); k++) {
      minColumn = min(minColumn, m_data[series][i][k]);
    }
    return minColumn;
  }
  
  double Store::squaredValueSumOfColumn(int series, int i) const {
    double result = 0;
    for (int k = 0; k < numberOfPairsOfSeries(series); k++) {
      result += m_data[series][i][k]*m_data[series][i][k];
    }
    return result;
  }
  
  double Store::columnProductSum(int series) const {
    double result = 0;
    for (int k = 0; k < numberOfPairsOfSeries(series); k++) {
      result += m_data[series][0][k]*m_data[series][1][k];
    }
    return result;
  }
  
  double Store::meanOfColumn(int series, int i) const {
    return numberOfPairsOfSeries(series) == 0 ? 0 : sumOfColumn(series, i)/numberOfPairsOfSeries(series);
  }
  
  double Store::varianceOfColumn(int series, int i) const {
    double mean = meanOfColumn(series, i);
    return squaredValueSumOfColumn(series, i)/numberOfPairsOfSeries(series) - mean*mean;
  }
  
  double Store::standardDeviationOfColumn(int series, int i) const {
    return std::sqrt(varianceOfColumn(series, i));
  }
  
  double Store::covariance(int series) const {
    return columnProductSum(series)/numberOfPairsOfSeries(series) - meanOfColumn(series, 0)*meanOfColumn(series, 1);
  }
  
  double Store::slope(int series) const {
    return covariance(series)/varianceOfColumn(series, 0);
  }
  
  double Store::yIntercept(int series) const {
    return meanOfColumn(series, 1) - slope(series)*meanOfColumn(series, 0);
  }
  
  double Store::yValueForXValue(int series, double x, Poincare::Context * globalContext) {
    Model * model = m_regressionModels[(int)m_regressionTypes[series]];
    double * coefficients = coefficientsForSeries(series, globalContext);
    return model->evaluate(coefficients, x);
  }
  
  double Store::xValueForYValue(int series, double y, Poincare::Context * globalContext) {
    Model * model = m_regressionModels[(int)m_regressionTypes[series]];
    double * coefficients = coefficientsForSeries(series, globalContext);
    return model->levelSet(coefficients, xMin(), xGridUnit()/10.0, xMax(), y, globalContext);
  }
  
  double Store::correlationCoefficient(int series) const {
    double sd0 = standardDeviationOfColumn(series, 0);
    double sd1 = standardDeviationOfColumn(series, 1);
    return (sd0 == 0.0 || sd1 == 0.0) ? 1.0 : covariance(series)/(sd0*sd1);
  }
  
  double Store::squaredCorrelationCoefficient(int series) const {
    double cov = covariance(series);
    double v0 = varianceOfColumn(series, 0);
    double v1 = varianceOfColumn(series, 1);
    return (v0 == 0.0 || v1 == 0.0) ? 1.0 : cov*cov/(v0*v1);
  }
  
  }