[apps/regression] Levenberg Marquardt regression algorithm

apps/regression/Makefile

@@ -12,10 +12,23 @@ app_objs += $(addprefix apps/regression/,\
   graph_controller.o\
   graph_view.o\
   initialisation_parameter_controller.o\
+  levenberg_marquardt.o\
   prediction_parameter_controller.o\
+  regression_context.o\
   store.o\
   store_controller.o\
-  regression_context.o\
+)
+
+app_objs += $(addprefix apps/regression/model/,\
+  cubic_model.o\
+  exponential_model.o\
+  linear_model.o\
+  logarithmic_model.o\
+  logistic_model.o\
+  power_model.o\
+  quadratic_model.o\
+  quartic_model.o\
+  trigonometric_model.o\
 )
 
 i18n_files += $(addprefix apps/regression/,\

apps/regression/calculation_controller.cpp

@@ -6,6 +6,7 @@
 #include "../../poincare/src/layout/vertical_offset_layout.h"
 #include <poincare.h>
 #include <assert.h>
+#include "levenberg_marquardt.h"
 
 using namespace Poincare;
 using namespace Shared;
@@ -47,6 +48,9 @@ bool CalculationController::handleEvent(Ion::Events::Event event) {
 }
 
 void CalculationController::didBecomeFirstResponder() {
+  LevenbergMarquartd lm(m_store);
+  double modelCoefficients[] = {1.0, 1.0, 1.0, 1.0, 1.0};
+  lm.fit(modelCoefficients, const_cast<AppsContainer *>(static_cast<const AppsContainer *>(app()->container()))->globalContext());
   if (selectedRow() == -1) {
     selectCellAtLocation(1, 0);
   } else {

apps/regression/levenberg_marquardt.cpp

@@ -0,0 +1,133 @@
+#include "levenberg_marquardt.h"
+#include <poincare/complex.h>
+#include <math.h>
+
+using namespace Poincare;
+using namespace Shared;
+
+namespace Regression {
+
+void LevenbergMarquartd::fit(double * modelCoefficients, Context * context) {
+  double currentChi2 = chi2(modelCoefficients);
+  double lambda = k_initialLambda;
+  int n = m_model.numberOfCoefficients(); // n unknown coefficients
+  int smallChi2ChangeCounts = 0;
+  int iterationCount = 0;
+  while (smallChi2ChangeCounts < k_consecutiveSmallChi2ChangesLimit && iterationCount < k_maxIterations) {
+    // Compute modelCoefficients step
+    // Create the alpha matrix
+    Expression * coefficientsAPrime[RegressionModel::k_maxNumberOfCoefficients][RegressionModel::k_maxNumberOfCoefficients]; // TODO find a way not to use so much space, we only need Expression * coefficientsAPrime[n][n]
+    // TODO compute only half the matrix, it is symmetric.
+    for (int i = 0; i < n; i++) {
+      for (int j = 0; j < n; j++) {
+        coefficientsAPrime[i][j] = new Complex<double>(Complex<double>::Float((alphaPrimeCoefficient(modelCoefficients, i, j, lambda))));
+      }
+    }
+    // Create the beta matrix
+    Expression ** operandsB = new Expression * [n];
+    for (int j = 0; j < n; j++) {
+      operandsB[j] = new Complex<double>(Complex<double>::Float(betaCoefficient(modelCoefficients, j)));
+    }
+    double modelCoefficientSteps[n];
+    solveLinearSystem(modelCoefficientSteps, coefficientsAPrime, operandsB, n, context);
+
+    // Compare new chi2 with the previous value
+    double newModelCoefficients[n];
+    for (int i = 0; i < n; i++) {
+      newModelCoefficients[i] = modelCoefficients[i] + modelCoefficientSteps[i];
+    }
+    double newChi2 = chi2(newModelCoefficients);
+    smallChi2ChangeCounts = (fabs(currentChi2 - newChi2) > k_chi2ChangeCondition) ? 0 : smallChi2ChangeCounts + 1;
+    if (newChi2 >= currentChi2) {
+      lambda*= k_lambdaFactor;
+    } else {
+      lambda/= k_lambdaFactor;
+      for (int i = 0; i < n; i++) {
+        modelCoefficients[i] = newModelCoefficients[i];
+      }
+      currentChi2 = newChi2;
+    }
+    iterationCount++;
+  }
+}
+
+double LevenbergMarquartd::chi2(double * modelCoefficients) const {
+  double result = 0.0;
+  for (int i = 0; i < m_store->numberOfPairsOfSeries(m_series); i++) {
+    double xi = m_store->get(m_series, 0, i);
+    double yi = m_store->get(m_series, 1, i);
+    double difference = yi - m_model.evaluate(modelCoefficients, xi);
+    result += difference * difference;
+  }
+  return result;
+}
+
+// a'(k,k) = a(k,k) * (1 + lambda)
+// a'(k,l) = a(l,k) when (k != l)
+double LevenbergMarquartd::alphaPrimeCoefficient(double * modelCoefficients, int k, int l, double lambda) const {
+  assert(k >= 0 && k < m_model.numberOfCoefficients());
+  assert(l >= 0 && l < m_model.numberOfCoefficients());
+  double result = k == l ? alphaCoefficient(modelCoefficients, k, l)*(1+lambda) : alphaCoefficient(modelCoefficients, l, k);
+  return result;
+}
+
+// a(k,l) = sum(0, N-1, derivate(y(xi|a), ak) * derivate(y(xi|a), a))
+double LevenbergMarquartd::alphaCoefficient(double * modelCoefficients, int k, int l) const {
+  assert(k >= 0 && k < m_model.numberOfCoefficients());
+  assert(l >= 0 && l < m_model.numberOfCoefficients());
+  double result = 0.0;
+  int m = m_store->numberOfPairsOfSeries(m_series);
+  for (int i = 0; i < m; i++) {
+    double xi = m_store->get(m_series, 0, i);
+    result += m_model.partialDerivate(modelCoefficients, k, xi) * m_model.partialDerivate(modelCoefficients, l, xi);
+  }
+  return result;
+}
+
+// b(k) = sum(0, N-1, (yi - y(xi|a)) * derivate(y(xi|a), ak))
+double LevenbergMarquartd::betaCoefficient(double * modelCoefficients, int k) const {
+  assert(k >= 0 && k < m_model.numberOfCoefficients());
+  double result = 0.0;
+  int m = m_store->numberOfPairsOfSeries(m_series); // m equations
+  for (int i = 0; i < m; i++) {
+    double xi = m_store->get(m_series, 0, i);
+    double yi = m_store->get(m_series, 1, i);
+    result += (yi - m_model.evaluate(modelCoefficients, xi)) * m_model.partialDerivate(modelCoefficients, k, xi);
+  }
+  return result;
+}
+
+
+// TODO should return an error if no solution ?
+void LevenbergMarquartd::solveLinearSystem(double * solutions, Expression * coefficients[RegressionModel::k_maxNumberOfCoefficients][RegressionModel::k_maxNumberOfCoefficients], Expression * * constants, int solutionDimension, Context * context) {
+  int n = solutionDimension;
+  const Expression ** operandsA = new const Expression * [n*n];
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      operandsA[i*n+j] = coefficients[i][j];
+    }
+  }
+  Matrix * AMatrix = new Matrix(operandsA, n, n, false);
+  delete[] operandsA;
+  Matrix * BMatrix = new Matrix(constants, n, 1, false);
+  Expression::AngleUnit angleUnit = Preferences::sharedPreferences()->angleUnit();
+  Matrix * matrixInverse = AMatrix->createApproximateInverse<double>();
+  assert(matrixInverse != nullptr);
+  Matrix * result = Multiplication::computeOnMatrices<double>(const_cast<const Matrix *>(matrixInverse), const_cast<const Matrix *>(BMatrix));
+  Expression * exactSolutions[n];
+  for (int i = 0; i < n; i++) {
+    Expression * sol = result->matrixOperand(i,0);
+    exactSolutions[i] = sol;
+    result->detachOperand(sol);
+    Expression::Simplify(&exactSolutions[i], *context, angleUnit);
+    assert(exactSolutions[i] != nullptr);
+    solutions[i] = exactSolutions[i]->approximateToScalar<double>(*context, angleUnit);
+    delete exactSolutions[i];
+  }
+  delete result;
+  delete matrixInverse;
+  delete BMatrix;
+  delete AMatrix;
+}
+
+}

apps/regression/levenberg_marquardt.h

@@ -0,0 +1,49 @@
+#ifndef REGRESSION_LEVENBERG_MARQUARDT_H
+#define REGRESSION_LEVENBERG_MARQUARDT_H
+
+#include <poincare.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/regression_model.h"
+#include "model/trigonometric_model.h"
+#include "store.h"
+
+namespace Regression {
+
+class LevenbergMarquartd {
+public:
+  LevenbergMarquartd(Store * store) :
+    m_store(store),
+    m_series(0),
+    m_model()
+  {
+  }
+  void fit(double * modelCoefficients, Poincare::Context * context);
+private:
+  static constexpr double k_maxIterations = 1000;
+  static constexpr double k_initialLambda = 0.001;
+  static constexpr double k_lambdaFactor = 10;
+  static constexpr double k_chi2ChangeCondition = 0.0005;
+  static constexpr int k_consecutiveSmallChi2ChangesLimit = 10;
+
+  double chi2(double * modelCoefficients) const;
+  double alphaPrimeCoefficient(double * modelCoefficients, int k, int l, double lambda) const;
+  double alphaCoefficient(double * modelCoefficients, int k, int l) const;
+  double betaCoefficient(double * modelCoefficients, int k) const;
+  void solveLinearSystem(double * solutions, Poincare::Expression * coefficients[RegressionModel::k_maxNumberOfCoefficients][RegressionModel::k_maxNumberOfCoefficients], Poincare::Expression * * constants, int solutionDimension, Poincare::Context * context);
+
+  Store * m_store;
+  int m_series; //TODO put as argument
+  LogisticModel m_model; //TODO put as argument
+};
+
+}
+
+
+#endif

apps/regression/model/cubic_model.cpp

@@ -0,0 +1,36 @@
+#include "cubic_model.h"
+#include <math.h>
+#include <assert.h>
+
+namespace Regression {
+
+double CubicModel::evaluate(double * modelCoefficients, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  double c = modelCoefficients[2];
+  double d = modelCoefficients[3];
+  return a*x*x*x+b*x*x+c*x+d;
+}
+
+double CubicModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
+  if (derivateCoefficientIndex == 0) {
+    // Derivate: x^3
+    return x*x*x;
+  }
+  if (derivateCoefficientIndex == 1) {
+    // Derivate: x^2
+    return x*x;
+  }
+  if (derivateCoefficientIndex == 2) {
+    // Derivate: x
+    return x;
+  }
+  if (derivateCoefficientIndex == 3) {
+    // Derivate: 1
+    return 1;
+  }
+  assert(false);
+  return 0.0;
+}
+
+}

apps/regression/model/cubic_model.h

@@ -0,0 +1,19 @@
+#ifndef REGRESSION_CUBIC_MODEL_H
+#define REGRESSION_CUBIC_MODEL_H
+
+#include "regression_model.h"
+
+namespace Regression {
+
+class CubicModel : public RegressionModel {
+public:
+  using RegressionModel::RegressionModel;
+  double evaluate(double * modelCoefficients, double x) const override;
+  double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
+  double numberOfCoefficients() const override { return 4; }
+};
+
+}
+
+
+#endif

apps/regression/model/exponential_model.cpp

@@ -0,0 +1,28 @@
+#include "exponential_model.h"
+#include <math.h>
+#include <assert.h>
+
+namespace Regression {
+
+double ExponentialModel::evaluate(double * modelCoefficients, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  return a*exp(b*x);
+}
+
+double ExponentialModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  if (derivateCoefficientIndex == 0) {
+    // Derivate: exp(b*x)
+    return exp(b*x);
+  }
+  if (derivateCoefficientIndex == 1) {
+    // Derivate: a*x*exp(b*x)
+    return a*x*exp(b*x);
+  }
+  assert(false);
+  return 0.0;
+}
+
+}

apps/regression/model/exponential_model.h

@@ -0,0 +1,19 @@
+#ifndef REGRESSION_EXPONENTIAL_MODEL_H
+#define REGRESSION_EXPONENTIAL_MODEL_H
+
+#include "regression_model.h"
+
+namespace Regression {
+
+class ExponentialModel : public RegressionModel {
+public:
+  using RegressionModel::RegressionModel;
+  double evaluate(double * modelCoefficients, double x) const override;
+  double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
+  double numberOfCoefficients() const override { return 2; }
+};
+
+}
+
+
+#endif

apps/regression/model/linear_model.cpp

@@ -0,0 +1,26 @@
+#include "linear_model.h"
+#include <math.h>
+#include <assert.h>
+
+namespace Regression {
+
+double LinearModel::evaluate(double * modelCoefficients, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  return a*x+b;
+}
+
+double LinearModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
+  if (derivateCoefficientIndex == 0) {
+    // Derivate: x
+    return x;
+  }
+  if (derivateCoefficientIndex == 1) {
+    // Derivate: 1;
+    return 1;
+  }
+  assert(false);
+  return 0.0;
+}
+
+}

apps/regression/model/linear_model.h

@@ -0,0 +1,19 @@
+#ifndef REGRESSION_LINEAR_MODEL_H
+#define REGRESSION_LINEAR_MODEL_H
+
+#include "regression_model.h"
+
+namespace Regression {
+
+class LinearModel : public RegressionModel {
+public:
+  using RegressionModel::RegressionModel;
+  double evaluate(double * modelCoefficients, double x) const override;
+  double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
+  double numberOfCoefficients() const override { return 2; }
+};
+
+}
+
+
+#endif

apps/regression/model/logarithmic_model.cpp

@@ -0,0 +1,26 @@
+#include "logarithmic_model.h"
+#include <math.h>
+#include <assert.h>
+
+namespace Regression {
+
+double LogarithmicModel::evaluate(double * modelCoefficients, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  return a*log(x)+b;
+}
+
+double LogarithmicModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
+  if (derivateCoefficientIndex == 0) {
+    // Derivate: ln(x)
+    return log(x);
+  }
+  if (derivateCoefficientIndex == 1) {
+    // Derivate: 1
+    return 1;
+  }
+  assert(false);
+  return 0.0;
+}
+
+}

apps/regression/model/logarithmic_model.h

@@ -0,0 +1,19 @@
+#ifndef REGRESSION_LOGARITHMIC_MODEL_H
+#define REGRESSION_LOGARITHMIC_MODEL_H
+
+#include "regression_model.h"
+
+namespace Regression {
+
+class LogarithmicModel : public RegressionModel {
+public:
+  using RegressionModel::RegressionModel;
+  double evaluate(double * modelCoefficients, double x) const override;
+  double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
+  double numberOfCoefficients() const override { return 2; }
+};
+
+}
+
+
+#endif

apps/regression/model/logistic_model.cpp

@@ -0,0 +1,35 @@
+#include "logistic_model.h"
+#include <math.h>
+#include <assert.h>
+
+namespace Regression {
+
+double LogisticModel::evaluate(double * modelCoefficients, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  double c = modelCoefficients[2];
+  return c/(1.0+a*exp(-b*x));
+}
+
+double LogisticModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  double c = modelCoefficients[2];
+  double denominator = 1.0+a*exp(-b*x);
+  if (derivateCoefficientIndex == 0) {
+    // Derivate: exp(-b*x)*(-1 * c/(1.0+a*exp(-b*x))^2)
+    return -exp(-b*x) * c/(denominator * denominator);
+  }
+  if (derivateCoefficientIndex == 1) {
+    // Derivate: (-x)*a*exp(-b*x)*(-1/(1.0+a*exp(-b*x))^2)
+    return x*a*exp(-b*x)*c/(denominator * denominator);
+  }
+  if (derivateCoefficientIndex == 2) {
+    // Derivate: (-x)*a*exp(-b*x)*(-1/(1.0+a*exp(-b*x))^2)
+    return 1.0/denominator;
+  }
+  assert(false);
+  return 0.0;
+}
+
+}

apps/regression/model/logistic_model.h

@@ -0,0 +1,19 @@
+#ifndef REGRESSION_LOGISTIC_MODEL_H
+#define REGRESSION_LOGISTIC_MODEL_H
+
+#include "regression_model.h"
+
+namespace Regression {
+
+class LogisticModel : public RegressionModel {
+public:
+  using RegressionModel::RegressionModel;
+  double evaluate(double * modelCoefficients, double x) const override;
+  double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
+  double numberOfCoefficients() const override { return 3; }
+};
+
+}
+
+
+#endif

apps/regression/model/power_model.cpp

@@ -0,0 +1,32 @@
+#include "power_model.h"
+#include <math.h>
+#include <assert.h>
+
+namespace Regression {
+
+double PowerModel::evaluate(double * modelCoefficients, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  return a*pow(x,b);
+}
+
+double PowerModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  if (derivateCoefficientIndex == 0) {
+    // Derivate: pow(x,b)
+    return pow(x,b);
+  }
+  if (derivateCoefficientIndex == 1) {
+    assert(x >= 0);
+    /* We assume all xi are positive.
+     * For x = 0, a*pow(x,b) = 0, the partial derivate along b is 0
+     * For x > 0, a*pow(x,b) = a*exp(b*ln(x)), the partial derivate along b is
+     *   ln(x)*a*pow(x,b) */
+    return x == 0 ? 0 : a*log(x)*exp(b*x);
+  }
+  assert(false);
+  return 0.0;
+}
+
+}

apps/regression/model/power_model.h

@@ -0,0 +1,19 @@
+#ifndef REGRESSION_POWER_MODEL_H
+#define REGRESSION_POWER_MODEL_H
+
+#include "regression_model.h"
+
+namespace Regression {
+
+class PowerModel : public RegressionModel {
+public:
+  using RegressionModel::RegressionModel;
+  double evaluate(double * modelCoefficients, double x) const override;
+  double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
+  double numberOfCoefficients() const override { return 2; }
+};
+
+}
+
+
+#endif

apps/regression/model/quadratic_model.cpp

@@ -0,0 +1,31 @@
+#include "quadratic_model.h"
+#include <math.h>
+#include <assert.h>
+
+namespace Regression {
+
+double QuadraticModel::evaluate(double * modelCoefficients, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  double c = modelCoefficients[2];
+  return a*x*x+b*x+c;
+}
+
+double QuadraticModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
+  if (derivateCoefficientIndex == 0) {
+    // Derivate: x^2
+    return x*x;
+  }
+  if (derivateCoefficientIndex == 1) {
+    // Derivate: x
+    return x;
+  }
+  if (derivateCoefficientIndex == 2) {
+    // Derivate: 1
+    return 1;
+  }
+  assert(false);
+  return 0.0;
+}
+
+}

apps/regression/model/quadratic_model.h

@@ -0,0 +1,19 @@
+#ifndef REGRESSION_QUADRATIC_MODEL_H
+#define REGRESSION_QUADRATIC_MODEL_H
+
+#include "regression_model.h"
+
+namespace Regression {
+
+class QuadraticModel : public RegressionModel {
+public:
+  using RegressionModel::RegressionModel;
+  double evaluate(double * modelCoefficients, double x) const override;
+  double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
+  double numberOfCoefficients() const override { return 3; }
+};
+
+}
+
+
+#endif

apps/regression/model/quartic_model.cpp

@@ -0,0 +1,41 @@
+#include "quartic_model.h"
+#include <math.h>
+#include <assert.h>
+
+namespace Regression {
+
+double QuarticModel::evaluate(double * modelCoefficients, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  double c = modelCoefficients[2];
+  double d = modelCoefficients[3];
+  double e = modelCoefficients[4];
+  return a*x*x*x*x+b*x*x*x+c*x*x+d*x+e;
+}
+
+double QuarticModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
+  if (derivateCoefficientIndex == 0) {
+    // Derivate: x^4
+    return x*x*x*x;
+  }
+  if (derivateCoefficientIndex == 1) {
+    // Derivate: x^3
+    return x*x*x;
+  }
+  if (derivateCoefficientIndex == 2) {
+    // Derivate: x^2
+    return x*x;
+  }
+  if (derivateCoefficientIndex == 3) {
+    // Derivate: x
+    return x;
+  }
+  if (derivateCoefficientIndex == 4) {
+    // Derivate: 1
+    return 1;
+  }
+  assert(false);
+  return 0.0;
+}
+
+}

apps/regression/model/quartic_model.h

@@ -0,0 +1,19 @@
+#ifndef REGRESSION_QUARTIC_MODEL_H
+#define REGRESSION_QUARTIC_MODEL_H
+
+#include "regression_model.h"
+
+namespace Regression {
+
+class QuarticModel : public RegressionModel {
+public:
+  using RegressionModel::RegressionModel;
+  double evaluate(double * modelCoefficients, double x) const override;
+  double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
+  double numberOfCoefficients() const override { return 5; }
+};
+
+}
+
+
+#endif

apps/regression/model/regression_model.h

@@ -0,0 +1,16 @@
+#ifndef REGRESSION_REGRESSION_MODEL_H
+#define REGRESSION_REGRESSION_MODEL_H
+
+namespace Regression {
+
+class RegressionModel {
+public:
+  static constexpr int k_maxNumberOfCoefficients = 5;
+  virtual double evaluate(double * modelCoefficients, double x) const = 0;
+  virtual double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const = 0;
+  virtual double numberOfCoefficients() const = 0;
+};
+
+}
+
+#endif

apps/regression/model/trigonometric_model.cpp

@@ -0,0 +1,34 @@
+#include "trigonometric_model.h"
+#include <math.h>
+#include <assert.h>
+
+namespace Regression {
+
+double TrigonometricModel::evaluate(double * modelCoefficients, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  double c = modelCoefficients[2];
+  return a*sin(b*x+c);
+}
+
+double TrigonometricModel::partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const {
+  double a = modelCoefficients[0];
+  double b = modelCoefficients[1];
+  double c = modelCoefficients[2];
+  if (derivateCoefficientIndex == 0) {
+    // Derivate: sin(b*x+c)
+    return sin(b*x+c);
+  }
+  if (derivateCoefficientIndex == 1) {
+    // Derivate: x*a*cos(b*x+c);
+    return x*a*cos(b*x+c);
+  }
+  if (derivateCoefficientIndex == 2) {
+    // Derivate: a*cos(b*x+c)
+    return a*cos(b*x+c);
+  }
+  assert(false);
+  return 0.0;
+}
+
+}

apps/regression/model/trigonometric_model.h

@@ -0,0 +1,19 @@
+#ifndef REGRESSION_TRIGONOMETRIC_MODEL_H
+#define REGRESSION_TRIGONOMETRIC_MODEL_H
+
+#include "regression_model.h"
+
+namespace Regression {
+
+class TrigonometricModel : public RegressionModel {
+public:
+  using RegressionModel::RegressionModel;
+  double evaluate(double * modelCoefficients, double x) const override;
+  double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const override;
+  double numberOfCoefficients() const override { return 3; }
+};
+
+}
+
+
+#endif