Commit 90e1f850 authored by Laurent Belcour's avatar Laurent Belcour

Got a compiling version of QuadProg/Eigen

parent 72d606cc
TEMPLATE = subdirs
SUBDIRS = rational_function \
rational_data \
rational_fitter_cgal \
rational_fitter_quadprog \
rational_fitter_eigen \
SUBDIRS = rational_function \
rational_data \
rational_fitter_cgal \
rational_fitter_quadprog \
rational_fitter_quadproge \
rational_fitter_eigen \
rational_fitter_matlab
rational_fitter_cgal.depends = rational_function rational_data
rational_fitter_quadprog.depends = rational_function rational_data
rational_fitter_eigen.depends = rational_function rational_data
rational_fitter_matlab.depends = rational_function rational_data
rational_fitter_cgal.depends = rational_function rational_data
rational_fitter_quadprog.depends = rational_function rational_data
rational_fitter_quadproge.depends = rational_function rational_data
rational_fitter_eigen.depends = rational_function rational_data
rational_fitter_matlab.depends = rational_function rational_data
This diff is collapsed.
#include "rational_fitter.h"
#include <Eigen/Dense>
#include <Eigen/SVD>
#include "eiquadprog.hpp"
#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
#include <QTime>
#include <rational_function.h>
#include <rational_data.h>
data* rational_fitter_quadprog::provide_data() const
{
return new rational_data() ;
}
function* rational_fitter_quadprog::provide_function() const
{
return new rational_function() ;
}
rational_fitter_quadprog::rational_fitter_quadprog() : QObject()
{
}
rational_fitter_quadprog::~rational_fitter_quadprog()
{
}
bool rational_fitter_quadprog::fit_data(const data* dat, function* fit)
{
rational_function* r = dynamic_cast<rational_function*>(fit) ;
const rational_data* d = dynamic_cast<const rational_data*>(dat) ;
if(r == NULL || d == NULL)
{
std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
return false ;
}
// I need to set the dimension of the resulting function to be equal
// to the dimension of my fitting problem
r->setDimX(d->dimX()) ;
r->setDimY(d->dimY()) ;
std::cout << "<<INFO>> np in [" << _min_np << ", " << _max_np
<< "] & nq in [" << _min_nq << ", " << _max_nq << "]" << std::endl ;
int temp_np = _min_np, temp_nq = _min_nq ;
while(temp_np <= _max_np || temp_nq <= _max_nq)
{
QTime time ;
time.start() ;
if(fit_data(d, temp_np, temp_nq, r))
{
int msec = time.elapsed() ;
int sec = (msec / 1000) % 60 ;
int min = (msec / 60000) % 60 ;
int hour = (msec / 3600000) ;
std::cout << "<<INFO>> got a fit using np = " << temp_np << " & nq = " << temp_nq << " " << std::endl ;
std::cout << "<<INFO>> it took " << hour << "h " << min << "m " << sec << "s" << std::endl ;
return true ;
}
std::cout << "<<INFO>> fit using np = " << temp_np << " & nq = " << temp_nq << " failed\r" ;
std::cout.flush() ;
if(temp_np <= _max_np)
{
++temp_np ;
}
if(temp_nq <= _max_nq)
{
++temp_nq ;
}
}
return false ;
}
void rational_fitter_quadprog::set_parameters(const arguments& args)
{
_max_np = args.get_float("np", 10) ;
_max_nq = args.get_float("nq", 10) ;
_min_np = args.get_float("min-np", _max_np) ;
_min_nq = args.get_float("min-nq", _max_nq) ;
}
bool rational_fitter_quadprog::fit_data(const rational_data* d, int np, int nq, rational_function* r)
{
// Multidimensional coefficients
std::vector<double> Pn ; Pn.reserve(d->dimY()*np) ;
std::vector<double> Qn ; Qn.reserve(d->dimY()*nq) ;
for(int j=0; j<d->dimY(); ++j)
{
if(!fit_data(d, np, nq, j, r))
return false ;
for(int i=0; i<np; ++i) { Pn.push_back(r->getP(i)) ; }
for(int i=0; i<nq; ++i) { Qn.push_back(r->getQ(i)) ; }
}
r->update(Pn, Qn) ;
return true ;
}
// dat is the data object, it contains all the points to fit
// np and nq are the degree of the RP to fit to the data
// y is the dimension to fit on the y-data (e.g. R, G or B for RGB signals)
// the function return a ration BRDF function and a boolean
bool rational_fitter_quadprog::fit_data(const rational_data* dat, int np, int nq, int ny, rational_function* rf)
{
rational_function* r = dynamic_cast<rational_function*>(rf) ;
const rational_data* d = dynamic_cast<const rational_data*>(dat) ;
if(r == NULL || d == NULL)
{
std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
return false ;
}
// I need to set the dimension of the resulting function to be equal
// to the dimension of my fitting problem
r->setDimX(d->dimX()) ;
r->setDimY(d->dimY()) ;
// 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() ;
// Matrices of the problem
Eigen::MatrixXd G (np+nq, np+nq) ;
Eigen::VectorXd g (np+nq) ;
Eigen::MatrixXd CI(np+nq, 2*d->size()) ;
Eigen::VectorXd ci(2*d->size()) ;
Eigen::MatrixXd CE(np+nq, 2*d->size()) ;
Eigen::VectorXd ce(2*d->size()) ;
Eigen::MatrixXd eCI(np+nq, 2*d->size()) ;
// Select the size of the result vector to
// be equal to the dimension of p + q
for(int i=0; i<np+nq; ++i)
{
G(i,i) = 1.0 ;
}
// Each constraint (fitting interval or point
// add another dimension to the constraint
// matrix
for(int i=0; i<d->size(); ++i)
{
// Norm of the row vector
double a0_norm = 0.0 ;
double a1_norm = 0.0 ;
vec xi = d->get(i) ;
for(int k=0; k<d->dimX(); ++k)
{
xi[k] /= dmax[k] ;
}
// A row of the constraint matrix has this
// form: [p_{0}(x_i), .., p_{np}(x_i), -f(x_i) q_{0}(x_i), .., -f(x_i) q_{nq}(x_i)]
// For the lower constraint and negated for
// the upper constraint
for(int j=0; j<np+nq; ++j)
{
// Filling the p part
if(j<np)
{
const double pi = r->p(xi, j) ;
a0_norm += pi*pi ;
a1_norm += pi*pi ;
CI(j,i) = pi ;
CI(j,i+d->size()) = -pi ;
}
// Filling the q part
else
{
vec yl, yu ;
d->get(i, yl, yu) ;
const double qi = r->q(xi, j-np) ;
a0_norm += qi*qi * (yl[ny]*yl[ny]) ;
CI(j,i) = -yl[ny] * qi ;
a1_norm += qi*qi * (yu[ny]*yu[ny]) ;
CI(j,i+d->size()) = yu[ny] * qi ;
}
}
// Set the c vector, will later be updated using the
// delta parameter.
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)
{
for(int i=0; i<np+nq; ++i)
{
std::cout << CI(i,j) ;
if(i != np+nq-1) std::cout << ", ";
}
if(j != d->size()*2-1)
std::cout << ";" << std::endl;
else
std::cout << "]" << std::endl ;
}
#endif
// Update the ci column with the delta parameter
// (See Celis et al. 2007 p.12)
Eigen::JacobiSVD<Eigen::MatrixXd> svd(CI);
const double sigma_m = svd.singularValues()(std::min(2*d->size(), np+nq)-1) ;
const double sigma_M = svd.singularValues()(0) ;
#ifdef DEBUG
std::cout << "<<DEBUG>> SVD = [ " ;
for(int i=0; i<std::min(2*d->size(), np+nq); ++i)
{
std::cout << svd.singularValues()(i) << ", " ;
}
std::cout << " ]" << std::endl ;
#endif
double delta = sigma_m / sigma_M ;
if(std::isnan(delta) || (std::abs(delta) == std::numeric_limits<double>::infinity()))
{
#ifdef DEBUG
std::cerr << "<<ERROR>> delta factor is NaN of Inf" << std::endl ;
#endif
return false ;
}
else if(delta == 0.0)
{
delta = 1.0 ;
}
#ifndef DEBUG
std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
#endif
for(int i=0; i<2*d->size(); ++i)
{
ci(i) = ci(i) * delta ;
}
// Compute the solution
Eigen::VectorXd x;
double cost = solve_quadprog(G, g, CE, ce, CI, ci, x);
bool solves_qp = !(cost == std::numeric_limits<double>::infinity());
for(int i=0; i<np+nq; ++i)
{
const double v = x[i];
solves_qp = solves_qp && !std::isnan(v) && (v != std::numeric_limits<double>::infinity()) ;
}
if(solves_qp)
{
// Recopy the vector d
std::vector<double> p, q;
double norm = 0.0 ;
for(int i=0; i<np+nq; ++i)
{
const double v = x[i];
norm += v*v ;
if(i < np)
{
p.push_back(v / r->p(dmax, i)) ;
}
else
{
q.push_back(v / r->q(dmax, i-np)) ;
}
}
r->update(p, q);
std::cout << "<<INFO>> got solution " << *r << std::endl ;
return norm > 0.0;
}
else
{
return false;
}
}
Q_EXPORT_PLUGIN2(rational_fitter_quadprog, rational_fitter_quadprog)
#pragma once
// Include STL
#include <vector>
#include <string>
// Interface
#include <QObject>
#include <core/function.h>
#include <core/data.h>
#include <core/fitter.h>
#include <core/args.h>
#include <rational_function.h>
#include <rational_data.h>
class rational_fitter_quadprog : public QObject, public fitter
{
Q_OBJECT
Q_INTERFACES(fitter)
public: // methods
rational_fitter_quadprog() ;
virtual ~rational_fitter_quadprog() ;
// Fitting a data object
//
virtual bool fit_data(const data* d, function* fit) ;
// Provide user parameters to the fitter
//
virtual void set_parameters(const arguments& args) ;
// Obtain associated data and functions
//
virtual data* provide_data() const ;
virtual function* provide_function() const ;
protected: // data
// Fitting a data object using np elements in the numerator and nq
// elements in the denominator
virtual bool fit_data(const rational_data* d, int np, int nq, rational_function* fit) ;
virtual bool fit_data(const rational_data* dat, int np, int nq, int ny, rational_function* fit) ;
// min and Max usable np and nq values for the fitting
int _max_np, _max_nq ;
int _min_np, _min_nq ;
} ;
TARGET = rational_fitter_quadproge
TEMPLATE = lib
CONFIG *= qt \
plugin \
eigen
DESTDIR = ../../build
INCLUDEPATH += ../rational_function \
../rational_data \
../..
HEADERS = rational_fitter.h
SOURCES = rational_fitter.cpp
LIBS += -L../../build \
-lrational_function \
-lrational_data
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