Commit 7af7ad29 authored by Laurent Belcour's avatar Laurent Belcour

Adding a quadratic program class to handle progressive algorithm

parent 46b4a7bd
......@@ -23,8 +23,8 @@ class plugins_manager
//! \brief Create the object, parse the argument and load all the plugins
plugins_manager(const arguments& args) ;
//! Get instances of the function, the data and the fitter. Select the first
//! in the map,
//! \brief Get instances of the function, the data and the fitter. Select
//! the first in the map,
function* get_function() const ;
data* get_data() const ;
fitter* get_fitter() const ;
......
#pragma once
#include <Eigen/SVD>
#include <Array.hh>
#include <QuadProg++.hh>
class quadratic_program
{
public:
//! \brief Constructor need to specify the number of coefficients
quadratic_program(int np, int nq) : _np(np), _nq(nq), CI(0.0, _np+_nq, 0) { }
//! \brief Remove the already defined constraints
void clear_constraints()
{
CI.resize(_np+_nq, 0);
}
//! \brief Add a constraint by specifying the vector
void add_constraints(const QuadProgPP::Vector<double> v)
{
const int m = CI.nrows();
const int n = CI.ncols();
CI.resize(m, n+1);
for(int i=0; i<m; ++i)
{
CI[i][n] = v[i];
}
}
//! \brief Solves the quadratic program
bool solve_program(QuadProgPP::Vector<double>& v)
{
const int m = CI.nrows();
const int n = CI.ncols();
QuadProgPP::Matrix<double> G (0.0, m, m) ;
QuadProgPP::Vector<double> g (0.0, m) ;
QuadProgPP::Vector<double> ci(0.0, n) ;
QuadProgPP::Matrix<double> CE(0.0, m, 0) ;
QuadProgPP::Vector<double> ce(0.0, 0) ;
// Update the ci column with the delta parameter
// (See Celis et al. 2007 p.12)
Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::HouseholderQRPreconditioner> svd(Eigen::MatrixXd::Map(&CI[0][0], m, n));
const double sigma_m = svd.singularValues()(std::min(m, n)-1) ;
const double sigma_M = svd.singularValues()(0) ;
const double delta = sigma_M / sigma_m ;
// Select the size of the result vector to
// be equal to the dimension of p + q
for(int i=0; i<m; ++i)
{
G[i][i] = 1.0 ;
}
// Each constraint (fitting interval or point
// add another dimension to the constraint
// matrix
for(int i=0; i<n; ++i)
{
// Norm of the row vector
double norm = 0.0 ;
for(int j=0; j<m; ++j)
{
norm += CI[j][i]*CI[j][i] ; ;
}
// Set the c vector, will later be updated using the
// delta parameter.
ci[i] = -delta * sqrt(norm) ;
}
// Compute the solution
const double cost = QuadProgPP::solve_quadprog(G, g, CE, ce, CI, ci, v);
bool solves_qp = !(cost == std::numeric_limits<double>::infinity());
return solves_qp;
}
protected:
int _np, _nq ;
QuadProgPP::Matrix<double> CI;
};
......@@ -12,9 +12,11 @@
#include <limits>
#include <algorithm>
#include <cmath>
#include <string>
#include <QTime>
#include "quadratic_program.h"
data* rational_fitter_parallel::provide_data() const
{
......@@ -83,7 +85,7 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit)
double min_delta = std::numeric_limits<double>::max();
int nb_sol_found = 0;
int np, nq ;
#pragma omp parallel for
// #pragma omp parallel for
for(int j=1; j<i; ++j)
{
int temp_np = i - j;
......@@ -168,6 +170,7 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
return false ;
}
#ifdef OLD
// Get the maximum value in data to scale the input parameter space
// so that it reduces the values of the polynomial
vec dmax = d->max() ;
......@@ -243,6 +246,7 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
ci[i] = -sqrt(a0_norm) ;
ci[i+d->size()] = -sqrt(a1_norm) ;
}
#ifdef DEBUG
std::cout << "CI = [" ;
for(int j=0; j<d->size()*2; ++j)
......@@ -313,7 +317,43 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* dat, int np, int
const double v = x[i];
solves_qp = solves_qp && !std::isnan(v) && (v != std::numeric_limits<double>::infinity()) ;
}
#else
quadratic_program qp(np, nq);
for(int i=0; i<d->size(); ++i)
{
vec xi = d->get(i) ;
// Create two vector of constraints
QuadProgPP::Vector<double> c1(np+nq), c2(np+nq);
for(int j=0; j<np+nq; ++j)
{
// Filling the p part
if(j<np)
{
const double pi = r->p(xi, j) ;
c1[j] = pi ;
c2[j] = -pi ;
}
// Filling the q part
else
{
vec yl, yu ;
d->get(i, yl, yu) ;
const double qi = r->q(xi, j-np) ;
c1[j] = -yl[ny] * qi ;
c2[j] = yu[ny] * qi ;
}
}
qp.add_constraints(c1);
qp.add_constraints(c2);
}
QuadProgPP::Vector<double> x(np+nq);
bool solves_qp = qp.solve_program(x);
#endif
if(solves_qp)
{
// Recopy the vector d
......
......@@ -12,7 +12,8 @@ INCLUDEPATH += ../rational_function \
../rational_data \
../..
HEADERS = rational_fitter.h
HEADERS = rational_fitter.h \
quadratic_program.h
SOURCES = rational_fitter.cpp
LIBS += -L../../build \
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment