From 1c8ff9bd83fbd29b57bf6fece82199d3165d7a17 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?L=C3=A9a=20Saviot?= Date: Wed, 6 Jun 2018 16:57:53 +0200 Subject: [PATCH] [apps/reg] Put levenberg-marquardt methods in Model --- apps/regression/Makefile | 2 +- apps/regression/calculation_controller.cpp | 4 -- apps/regression/levenberg_marquardt.h | 49 -------------- .../model.cpp} | 66 +++++++++++-------- apps/regression/model/model.h | 21 ++++++ 5 files changed, 59 insertions(+), 83 deletions(-) delete mode 100644 apps/regression/levenberg_marquardt.h rename apps/regression/{levenberg_marquardt.cpp => model/model.cpp} (59%) diff --git a/apps/regression/Makefile b/apps/regression/Makefile index 4c1c202a1..7acd5c7d7 100644 --- a/apps/regression/Makefile +++ b/apps/regression/Makefile @@ -12,7 +12,6 @@ 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\ regression_controller.o\ @@ -27,6 +26,7 @@ app_objs += $(addprefix apps/regression/model/,\ linear_model.o\ logarithmic_model.o\ logistic_model.o\ + model.o\ power_model.o\ quadratic_model.o\ quartic_model.o\ diff --git a/apps/regression/calculation_controller.cpp b/apps/regression/calculation_controller.cpp index bd200d79a..522eb77cd 100644 --- a/apps/regression/calculation_controller.cpp +++ b/apps/regression/calculation_controller.cpp @@ -6,7 +6,6 @@ #include "../../poincare/src/layout/vertical_offset_layout.h" #include #include -#include "levenberg_marquardt.h" using namespace Poincare; using namespace Shared; @@ -48,9 +47,6 @@ 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(static_cast(app()->container()))->globalContext()); if (selectedRow() == -1) { selectCellAtLocation(1, 0); } else { diff --git a/apps/regression/levenberg_marquardt.h b/apps/regression/levenberg_marquardt.h deleted file mode 100644 index 99985ef0f..000000000 --- a/apps/regression/levenberg_marquardt.h +++ /dev/null @@ -1,49 +0,0 @@ -#ifndef REGRESSION_LEVENBERG_MARQUARDT_H -#define REGRESSION_LEVENBERG_MARQUARDT_H - -#include -#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/model.h" -#include "model/power_model.h" -#include "model/quadratic_model.h" -#include "model/quartic_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[Model::k_maxNumberOfCoefficients][Model::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 diff --git a/apps/regression/levenberg_marquardt.cpp b/apps/regression/model/model.cpp similarity index 59% rename from apps/regression/levenberg_marquardt.cpp rename to apps/regression/model/model.cpp index 391d27b94..6c3107760 100644 --- a/apps/regression/levenberg_marquardt.cpp +++ b/apps/regression/model/model.cpp @@ -1,5 +1,8 @@ -#include "levenberg_marquardt.h" +#include "model.h" +#include "../store.h" #include +#include +#include #include using namespace Poincare; @@ -7,10 +10,14 @@ using namespace Shared; namespace Regression { -void LevenbergMarquartd::fit(double * modelCoefficients, Context * context) { - double currentChi2 = chi2(modelCoefficients); +void Model::fit(Store * store, int series, double * modelCoefficients, Poincare::Context * context) { + fitLevenbergMarquardt(store, series, modelCoefficients, context); +} + +void Model::fitLevenbergMarquardt(Store * store, int series, double * modelCoefficients, Context * context) { + double currentChi2 = chi2(store, series, modelCoefficients); double lambda = k_initialLambda; - int n = m_model.numberOfCoefficients(); // n unknown coefficients + int n = numberOfCoefficients(); // n unknown coefficients int smallChi2ChangeCounts = 0; int iterationCount = 0; while (smallChi2ChangeCounts < k_consecutiveSmallChi2ChangesLimit && iterationCount < k_maxIterations) { @@ -19,7 +26,7 @@ void LevenbergMarquartd::fit(double * modelCoefficients, Context * context) { Expression * coefficientsAPrime[Model::k_maxNumberOfCoefficients][Model::k_maxNumberOfCoefficients]; // TODO find a way not to use so much space, we only need Expression * coefficientsAPrime[n][n] for (int i = 0; i < n; i++) { for (int j = i; j < n; j++) { - Complex alphaPrime = Complex::Float(alphaPrimeCoefficient(modelCoefficients, i, j, lambda)); + Complex alphaPrime = Complex::Float(alphaPrimeCoefficient(store, series, modelCoefficients, i, j, lambda)); coefficientsAPrime[i][j] = new Complex(alphaPrime); if (i != j) { coefficientsAPrime[j][i] = new Complex(alphaPrime); @@ -29,7 +36,7 @@ void LevenbergMarquartd::fit(double * modelCoefficients, Context * context) { // Create the beta matrix Expression ** operandsB = new Expression * [n]; for (int j = 0; j < n; j++) { - operandsB[j] = new Complex(Complex::Float(betaCoefficient(modelCoefficients, j))); + operandsB[j] = new Complex(Complex::Float(betaCoefficient(store, series, modelCoefficients, j))); } double modelCoefficientSteps[n]; solveLinearSystem(modelCoefficientSteps, coefficientsAPrime, operandsB, n, context); @@ -39,7 +46,7 @@ void LevenbergMarquartd::fit(double * modelCoefficients, Context * context) { for (int i = 0; i < n; i++) { newModelCoefficients[i] = modelCoefficients[i] + modelCoefficientSteps[i]; } - double newChi2 = chi2(newModelCoefficients); + double newChi2 = chi2(store, series, newModelCoefficients); smallChi2ChangeCounts = (fabs(currentChi2 - newChi2) > k_chi2ChangeCondition) ? 0 : smallChi2ChangeCounts + 1; if (newChi2 >= currentChi2) { lambda*= k_lambdaFactor; @@ -54,12 +61,12 @@ void LevenbergMarquartd::fit(double * modelCoefficients, Context * context) { } } -double LevenbergMarquartd::chi2(double * modelCoefficients) const { +double Model::chi2(Store * store, int series, 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); + for (int i = 0; i < store->numberOfPairsOfSeries(series); i++) { + double xi = store->get(series, 0, i); + double yi = store->get(series, 1, i); + double difference = yi - evaluate(modelCoefficients, xi); result += difference * difference; } return result; @@ -67,42 +74,42 @@ double LevenbergMarquartd::chi2(double * modelCoefficients) const { // 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); +double Model::alphaPrimeCoefficient(Store * store, int series, double * modelCoefficients, int k, int l, double lambda) const { + assert(k >= 0 && k < numberOfCoefficients()); + assert(l >= 0 && l < numberOfCoefficients()); + double result = k == l ? alphaCoefficient(store, series, modelCoefficients, k, l)*(1+lambda) : alphaCoefficient(store, series, 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 Model::alphaCoefficient(Store * store, int series, double * modelCoefficients, int k, int l) const { + assert(k >= 0 && k < numberOfCoefficients()); + assert(l >= 0 && l < numberOfCoefficients()); double result = 0.0; - int m = m_store->numberOfPairsOfSeries(m_series); + int m = store->numberOfPairsOfSeries(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); + double xi = store->get(series, 0, i); + result += partialDerivate(modelCoefficients, k, xi) * 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 Model::betaCoefficient(Store * store, int series, double * modelCoefficients, int k) const { + assert(k >= 0 && k < numberOfCoefficients()); double result = 0.0; - int m = m_store->numberOfPairsOfSeries(m_series); // m equations + int m = store->numberOfPairsOfSeries(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); + double xi = store->get(series, 0, i); + double yi = store->get(series, 1, i); + result += (yi - evaluate(modelCoefficients, xi)) * partialDerivate(modelCoefficients, k, xi); } return result; } // TODO should return an error if no solution ? -void LevenbergMarquartd::solveLinearSystem(double * solutions, Expression * coefficients[Model::k_maxNumberOfCoefficients][Model::k_maxNumberOfCoefficients], Expression * * constants, int solutionDimension, Context * context) { +void Model::solveLinearSystem(double * solutions, Expression * coefficients[Model::k_maxNumberOfCoefficients][Model::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++) { @@ -134,3 +141,4 @@ void LevenbergMarquartd::solveLinearSystem(double * solutions, Expression * coef } } + diff --git a/apps/regression/model/model.h b/apps/regression/model/model.h index a807516bc..3b2577e95 100644 --- a/apps/regression/model/model.h +++ b/apps/regression/model/model.h @@ -2,9 +2,13 @@ #define REGRESSION_MODEL_H #include +#include +#include namespace Regression { +class Store; + class Model { public: enum class Type : uint8_t { @@ -18,10 +22,27 @@ public: Trigonometric = 7, Logistic = 8 }; + static constexpr int k_numberOfModels = 9; static constexpr int k_maxNumberOfCoefficients = 5; virtual double evaluate(double * modelCoefficients, double x) const = 0; + virtual void fit(Store * store, int series, double * modelCoefficients, Poincare::Context * context); +private: + // Model attributes virtual double partialDerivate(double * modelCoefficients, int derivateCoefficientIndex, double x) const = 0; virtual double numberOfCoefficients() const = 0; + + // Levenberg-Marquardt + 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; + void fitLevenbergMarquardt(Store * store, int series, double * modelCoefficients, Poincare::Context * context); + double chi2(Store * store, int series, double * modelCoefficients) const; + double alphaPrimeCoefficient(Store * store, int series, double * modelCoefficients, int k, int l, double lambda) const; + double alphaCoefficient(Store * store, int series, double * modelCoefficients, int k, int l) const; + double betaCoefficient(Store * store, int series, double * modelCoefficients, int k) const; + void solveLinearSystem(double * solutions, Poincare::Expression * coefficients[Model::k_maxNumberOfCoefficients][Model::k_maxNumberOfCoefficients], Poincare::Expression * * constants, int solutionDimension, Poincare::Context * context); }; }