UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-03-18 21:30:38 +01:00
Code Issues Packages Projects Releases Wiki Activity

[apps/reg] Put levenberg-marquardt methods in Model

Browse Source
This commit is contained in:
Léa Saviot
2018-06-06 16:57:53 +02:00
committed by Émilie Feral
parent 86c482aff3
commit 1c8ff9bd83
5 changed files with 59 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
apps/regression/Makefile
View File

@@ -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\

4
apps/regression/calculation_controller.cpp
View File

@@ -6,7 +6,6 @@
#include "../../poincare/src/layout/vertical_offset_layout.h"
#include <poincare.h>
#include <assert.h>
#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<AppsContainer *>(static_cast<const AppsContainer *>(app()->container()))->globalContext());
if (selectedRow() == -1) {
selectCellAtLocation(1, 0);
} else {

49
apps/regression/levenberg_marquardt.h
View File

@@ -1,49 +0,0 @@
#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/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

66
apps/regression/levenberg_marquardt.cpp → apps/regression/model/model.cpp
View File

@@ -1,5 +1,8 @@
#include "levenberg_marquardt.h"
#include "model.h"
#include "../store.h"
#include <poincare/complex.h>
#include <poincare/matrix.h>
#include <poincare/multiplication.h>
#include <math.h>
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<double> alphaPrime = Complex<double>::Float(alphaPrimeCoefficient(modelCoefficients, i, j, lambda));
Complex<double> alphaPrime = Complex<double>::Float(alphaPrimeCoefficient(store, series, modelCoefficients, i, j, lambda));
coefficientsAPrime[i][j] = new Complex<double>(alphaPrime);
if (i != j) {
coefficientsAPrime[j][i] = new Complex<double>(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<double>(Complex<double>::Float(betaCoefficient(modelCoefficients, j)));
operandsB[j] = new Complex<double>(Complex<double>::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
}
}

21
apps/regression/model/model.h
View File

@@ -2,9 +2,13 @@
#define REGRESSION_MODEL_H
#include <stdint.h>
#include <poincare/context.h>
#include <poincare/expression.h>
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);
};
}