Commit d382b88f authored by Laurent Belcour's avatar Laurent Belcour
Browse files

Adding a working plugin format

parent dcdaffac
TEMPLATE = subdirs
SUBDIRS = rational_function \
rational_data \
rational_fitter_cgal
rational_fitter_cgal.depends = rational_function rational_data
#include "rational_data.h"
#include <boost/regex.hpp>
#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
void rational_data::load(const std::string& filename)
{
load(filename, -std::numeric_limits<double>::max(), std::numeric_limits<double>::max()) ;
}
void rational_data::load(const std::string& filename, double min, double max)
{
std::ifstream file(filename) ;
_min = std::numeric_limits<double>::max() ;
_max = -std::numeric_limits<double>::max() ;
if(!file.is_open())
{
std::cerr << "<<ERROR>> unable to open file \"" << filename << "\"" << std::endl ;
}
// N-Floats regexp
boost::regex e ("^([0-9]*\.?[0-9]+[\\t ]?)+");
double x, y, dy ;
while(file.good())
{
std::string line ;
std::getline(file, line) ;
std::stringstream linestream(line) ;
// Discard incorrect lines
if(!boost::regex_match(line,e))
{
continue ;
}
linestream >> x >> y ;
if(linestream.good()) {
linestream >> dy ;
} else {
// TODO Specify the delta in case
dy = 0.01f ;
}
if(x <= max && x >= min)
{
std::vector<double> v ;
v.push_back(x) ;
v.push_back(y-dy) ;
v.push_back(y+dy) ;
_data.push_back(v) ;
// Update min and max
_min = std::min(_min, x) ;
_max = std::max(_max, x) ;
}
}
//TODO Test for small data
/* std::vector<std::vector<double> > temp ;
std::ofstream temp_out("temp.gnuplot", std::ios_base::trunc) ;
for(int i=0; i<20; ++i)
{
int k = (i * _data.size()) / 20 ;
temp_out << _data[k][0] << "\t" << _data[k][1] << "\t" << 0.5*(_data[k][2] - _data[k][1]) << std::endl ;
temp.push_back(_data[k]) ;
}
_data = temp ;
*/
// Sort the vector
std::sort(_data.begin(), _data.end(), [](const std::vector<double>& a, const std::vector<double>& b){return (a[0]<b[0]);});
std::cout << "<<INFO>> loaded file \"" << filename << "\"" << std::endl ;
std::cout << "<<INFO>> data inside [" << _min << ", " << _max << "]" << std::endl ;
std::cout << "<<INFO>> " << _data.size() << " elements to fit" << std::endl ;
}
bool rational_data::get(int i, double& x, double& yl, double& yu) const
{
if(i >= (int)_data.size())
{
return false ;
}
x = _data[i][0] ;
yl = _data[i][1] ;
yu = _data[i][2] ;
return true ;
}
int rational_data::input_dimension() const
{
return 1 ;
}
int rational_data::output_dimension() const
{
return 1 ;
}
const std::vector<double>& rational_data::operator[](int i) const
{
return _data[i] ;
}
const std::vector<double>& rational_data::get(int i) const
{
return _data[i] ;
}
int rational_data::size() const
{
return _data.size() ;
}
double rational_data::min() const
{
return _min ;
}
double rational_data::max() const
{
return _max ;
}
Q_EXPORT_PLUGIN2(rational_data, rational_data)
#pragma once
// Include STL
#include <functional>
#include <vector>
#include <string>
#include <tuple>
// Interface
#include <QObject>
#include <core/function.h>
#include <core/data.h>
#include <core/fitter.h>
#include <core/args.h>
class rational_data : public QObject, public data
{
Q_OBJECT
Q_INTERFACES(data)
public: // methods
// Load data from a file
virtual void load(const std::string& filename) ;
virtual void load(const std::string& filename, double min, double max) ;
// Acces to data
virtual bool get(int i, double& x, double& yl, double& yu) const ;
virtual const std::vector<double>& get(int i) const ;
virtual const std::vector<double>& operator[](int i) const ;
// Get data size
virtual int size() const ;
// Get min and max input parameters
virtual double min() const ;
virtual double max() const ;
// Get the dimension of the input and output space
// WARNING: this dimension is defined after loading
// the data!
virtual int input_dimension() const ;
virtual int output_dimension() const ;
private: // data
// Store for each point of data, the upper
// and lower value
std::vector<std::vector<double> > _data ;
// Store the min and max value on the input
// domain
double _min, _max ;
} ;
TEMPLATE = lib
CONFIG += plugin
INCLUDEPATH += ../..
HEADERS = rational_data.h
SOURCES = rational_data.cpp
LIBS += -lboost_regex
QMAKE_CXXFLAGS += -std=c++11 -frounding-math
#include "rational_1d_fitter_cgal.h"
#include <CGAL/basic.h>
#include <CGAL/QP_models.h>
#include <CGAL/QP_functions.h>
#include <CGAL/MP_Float.h>
#include <Eigen/SVD>
#include <boost/regex.hpp>
#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
#include <rational_function.h>
#include <rational_data.h>
typedef CGAL::MP_Float ET ;
typedef CGAL::Quadratic_program<ET> Program ;
typedef CGAL::Quadratic_program_solution<ET> Solution ;
typedef CGAL::Quadratic_program_options Options ;
bool rational_1d_fitter_cgal::fit_data(const data* d, function*& fit)
{
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)
{
if(fit_data(d, temp_np, temp_nq, fit))
{
return true ;
}
std::cout << "<<INFO>> fitt using np = " << temp_np << " & nq = " << temp_nq << " failed\r" ;
if(temp_np <= _max_np)
{
++temp_np ;
}
if(temp_nq <= _max_nq)
{
++temp_nq ;
}
}
return false ;
}
void rational_1d_fitter_cgal::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_1d_fitter_cgal::fit_data(const data* d, int np, int nq, function*& rf)
{
// by default, we have a nonnegative QP with Ax - b >= 0
Program qp (CGAL::LARGER, false, 0, false, 0) ;
rational_function* r = dynamic_cast<rational_function*>(rf) ;
if(r == nullptr)
{
std::cerr << "<<ERROR>> not passing the correct class to the fitter" << std::endl ;
return false ;
}
// Select the size of the result vector to
// be equal to the dimension of p + q
for(int i=0; i<np+nq; ++i)
{
qp.set_d(i, i, 1.0) ;
}
// Each constraint (fitting interval or point
// add another dimension to the constraint
// matrix
Eigen::MatrixXd CI(np+nq, 2*d->size()) ;
Eigen::VectorXd ci(2*d->size()) ;
for(int i=0; i<d->size(); ++i)
{
// Norm of the row vector
double a0_norm = 0.0 ;
double a1_norm = 0.0 ;
// 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((*d)[i][0], j) ;
a0_norm += pi*pi ;
a1_norm += pi*pi ;
qp.set_a(j, i, ET(pi)) ;
qp.set_a(j, i+d->size(), ET(-pi)) ;
CI(j, i) = pi ;
CI(j, i+d->size()) = -pi ;
}
// Filling the q part
else
{
const double qi = r->q((*d)[i][0], j-np) ;
a0_norm += qi*qi * ((*d)[i][1]*(*d)[i][1]) ;
qp.set_a(j, i, ET(-(*d)[i][1] * qi)) ;
CI(j, i) = -(*d)[i][1] * qi ;
a1_norm += qi*qi * ((*d)[i][2]*(*d)[i][2]) ;
qp.set_a(j, i+d->size(), ET((*d)[i][2] * qi)) ;
CI(j, i+d->size()) = (*d)[i][2] * 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) ;
}
// 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
const double delta = sigma_m / sigma_M ;
if(isnan(delta) || (abs(delta) == std::numeric_limits<double>::infinity()))
{
#ifdef DEBUG
std::cerr << "<<ERROR>> delta factor is NaN of Inf" << std::endl ;
#endif
return false ;
}
#ifdef DEBUG
std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
#endif
for(int i=0; i<2*d->size(); ++i)
{
qp.set_b(i, ET(delta * ci(i))) ;
}
#ifdef DEBUG
// Export some informations on the problem to solve
std::cout << "<<DEBUG>> " << qp.get_n() << " variables" << std::endl ;
std::cout << "<<DEBUG>> " << qp.get_m() << " constraints" << std::endl ;
#endif
if(qp.is_linear() && !qp.is_nonnegative())
{
std::cerr << "<<ERROR>> the problem should not be linear or negative!" << std::endl ;
return false ;
}
#ifdef EXTREM_DEBUG
// Iterate over the rows
std::cout << std::endl ;
std::cout << "A = [" ;
for(int j=0; j<qp.get_m(); ++j)
{
if(j == 0) std::cout << " " ;
// Iterate over the columns
for(int i=0; i<qp.get_n(); ++i)
{
if(i > 0) std::cout << ",\t" ;
std::cout << *(*(qp.get_a()+i) +j) ;
}
std::cout << ";" ;
if(j < qp.get_n()-1) std::cout << std::endl ;
}
std::cout << "]" << std::endl << std::endl ;
std::cout << std::endl ;
std::cout << "D = [" ;
for(int j=0; j<np+nq; ++j)
{
if(j == 0) std::cout << " " ;
// Iterate over the columns
for(int i=0; i<np+nq; ++i)
{
if(i > 0) std::cout << ",\t" ;
std::cout << *(*(qp.get_d()+i) +j) ;
}
std::cout << ";" ;
}
std::cout << "]" << std::endl << std::endl ;
#endif
// solve the program, using ET as the exact type
Options o ;
o.set_auto_validation(true) ;
Solution s = CGAL::solve_quadratic_program(qp, ET(), o) ;
#ifdef DEBUG
if(s.is_infeasible())
{
std::cout << "<<DEBUG>> the current program is infeasible" << std::endl ;
}
#endif
bool solves_qp = !s.is_infeasible() && s.solves_quadratic_program(qp) ;
for(int i=0; i<np+nq; ++i)
{
const double v = (double)CGAL::to_double(*(s.variable_numerators_begin()+i)) ;
solves_qp = solves_qp && !isnan(v) && (v != std::numeric_limits<double>::infinity()) ;
}
if(solves_qp)
{
// Recopy the vector d
std::vector<double> p, q;
for(int i=0; i<np+nq; ++i)
{
const double v = CGAL::to_double(*(s.variable_numerators_begin()+i)) ;
if(i < np)
{
p.push_back(v) ;
}
else
{
q.push_back(v) ;
}
}
if(r != nullptr)
{
delete r ;
}
r = new rational_function(p, q);
return true;
}
else
{
return false;
}
}
Q_EXPORT_PLUGIN2(rational_plugin_cgal, rational_1d_fitter_cgal)
#pragma once
// Include STL
#include <functional>
#include <vector>
#include <string>
#include <tuple>
// Interface
#include <QObject>
#include <core/function.h>
#include <core/data.h>
#include <core/fitter.h>
#include <core/args.h>
class rational_1d_fitter_cgal : public QObject, public fitter
{
Q_OBJECT
Q_INTERFACES(fitter)
public: // methods
// Fitting a data object
virtual bool fit_data(const data* d, function*& fit) ;
// Fitting a data object using np elements in the numerator and nq
// elements in the denominator
virtual bool fit_data(const data* d, int np, int nq, function*& fit) ;
virtual void set_parameters(const arguments& args) ;
protected: // data
// min and Max usable np and nq values for the fitting
int _max_np, _max_nq ;
int _min_np, _min_nq ;
} ;
TEMPLATE = lib
CONFIG += plugin
INCLUDEPATH += ../rational_function ../rational_data ../.. /home/belcour/Sources/Eigen/include/eigen3
HEADERS = rational_1d_fitter_cgal.h
SOURCES = rational_1d_fitter_cgal.cpp
LIBS += -lCGAL -lboost_regex \
-L../rational_function -L../rational_data \
-lrational_function -lrational_data
QMAKE_CXXFLAGS += -std=c++11 -frounding-math
#include "rational_function.h"
#include <boost/regex.hpp>
#include <string>
#include <iostream>
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
rational_function::rational_function()
{
}
rational_function::rational_function(const std::vector<double>& a,
const std::vector<double>& b) :
a(a), b(b)
{
}
rational_function::~rational_function()
{
}
// Overload the function operator
double rational_function::operator()(double x) const
{
double p = 0.0f ;
double q = 0.0f ;
for(int i=a.size()-1; i>=0; --i)
{
p = x*p + a[i] ;
}
for(int i=b.size()-1; i>=0; --i)
{
q = x*q + b[i] ;
}
return p/q ;
}
// Get the p_i and q_j function
double rational_function::p(double x, int i) const
{
return pow(x, i) ;
}
double rational_function::q(double x, int j) const
{
return pow(x, j) ;
}
// IO function to text files
void rational_function::load(const std::string& filename)
{
}
void rational_function::save() const
{
}
std::ostream& operator<< (std::ostream& out, const rational_function& r)
{
std::cout << "p = [" ;
for(int i=0; i<r.a.size(); ++i)
{