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

Working CGAL version with check of the solution.

parent 19e4e31c
......@@ -11,8 +11,8 @@ rational_1d::rational_1d()
{
}
rational_1d::rational_1d(const std::vector<float>& a,
const std::vector<float>& b) :
rational_1d::rational_1d(const std::vector<double>& a,
const std::vector<double>& b) :
a(a), b(b)
{
}
......@@ -21,10 +21,10 @@ rational_1d::~rational_1d()
}
// Overload the function operator
float rational_1d::operator()(float x) const
double rational_1d::operator()(double x) const
{
float p = 0.0f ;
float q = 0.0f ;
double p = 0.0f ;
double q = 0.0f ;
for(int i=a.size()-1; i>=0; --i)
{
......@@ -40,11 +40,11 @@ float rational_1d::operator()(float x) const
}
// Get the p_i and q_j function
float rational_1d::p(float x, int i) const
double rational_1d::p(double x, int i) const
{
return pow(x, i) ;
}
float rational_1d::q(float x, int j) const
double rational_1d::q(double x, int j) const
{
return pow(x, j) ;
}
......@@ -85,15 +85,15 @@ std::ostream& operator<< (std::ostream& out, const rational_1d& r)
void rational_1d_data::load(const std::string& filename)
{
load(filename, -std::numeric_limits<float>::max(), std::numeric_limits<float>::max()) ;
load(filename, -std::numeric_limits<double>::max(), std::numeric_limits<double>::max()) ;
}
void rational_1d_data::load(const std::string& filename, float min, float max)
void rational_1d_data::load(const std::string& filename, double min, double max)
{
std::ifstream file(filename) ;
_min = std::numeric_limits<float>::max() ;
_max = -std::numeric_limits<float>::max() ;
_min = std::numeric_limits<double>::max() ;
_max = -std::numeric_limits<double>::max() ;
if(!file.is_open())
{
......@@ -103,7 +103,7 @@ void rational_1d_data::load(const std::string& filename, float min, float max)
// N-Floats regexp
boost::regex e ("^([0-9]*\.?[0-9]+[\\t ]?)+");
float x, y, dy ;
double x, y, dy ;
while(file.good())
{
std::string line ;
......@@ -121,12 +121,12 @@ void rational_1d_data::load(const std::string& filename, float min, float max)
linestream >> dy ;
} else {
// TODO Specify the delta in case
dy = 0.001f ;
dy = 0.1f ;
}
if(x <= max && x >= min)
{
std::vector<float> v ;
std::vector<double> v ;
v.push_back(x) ;
v.push_back(y-dy) ;
v.push_back(y+dy) ;
......@@ -139,13 +139,13 @@ void rational_1d_data::load(const std::string& filename, float min, float max)
}
// Sort the vector
std::sort(_data.begin(), _data.end(), [](const std::vector<float>& a, const std::vector<float>& b){return (a[0]<b[0]);});
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 ;
}
bool rational_1d_data::get(int i, float& x, float& yl, float& yu) const
bool rational_1d_data::get(int i, double& x, double& yl, double& yu) const
{
if(i >= (int)_data.size())
{
......@@ -159,7 +159,7 @@ bool rational_1d_data::get(int i, float& x, float& yl, float& yu) const
return true ;
}
const std::vector<float>& rational_1d_data::operator[](int i) const
const std::vector<double>& rational_1d_data::operator[](int i) const
{
return _data[i] ;
}
......@@ -169,12 +169,12 @@ int rational_1d_data::size() const
return _data.size() ;
}
float rational_1d_data::min() const
double rational_1d_data::min() const
{
return _min ;
}
float rational_1d_data::max() const
double rational_1d_data::max() const
{
return _max ;
}
......@@ -6,20 +6,20 @@
#include <string>
#include <tuple>
class rational_1d : public std::function<float(float)>
class rational_1d : public std::function<double(double)>
{
public: // methods
rational_1d() ;
rational_1d(const std::vector<float>& a, const std::vector<float>& b) ;
rational_1d(const std::vector<double>& a, const std::vector<double>& b) ;
virtual ~rational_1d() ;
// Overload the function operator
virtual float operator()(float x) const ;
virtual double operator()(double x) const ;
// Get the p_i and q_j function
virtual float p(float x, int i) const ;
virtual float q(float x, int j) const ;
virtual double p(double x, int i) const ;
virtual double q(double x, int j) const ;
// IO function to text files
void load(const std::string& filename) ;
......@@ -32,8 +32,8 @@ class rational_1d : public std::function<float(float)>
// Store the coefficients for the moment, I assume
// the functions to be polynomials.
std::vector<float> a ;
std::vector<float> b ;
std::vector<double> a ;
std::vector<double> b ;
} ;
class rational_1d_data // : public fitting_data
......@@ -42,28 +42,28 @@ class rational_1d_data // : public fitting_data
// Load data from a file
void load(const std::string& filename) ;
void load(const std::string& filename, float min, float max) ;
void load(const std::string& filename, double min, double max) ;
// Acces to data
bool get(int i, float& x, float& yl, float& yu) const ;
const std::vector<float>& operator[](int i) const ;
bool get(int i, double& x, double& yl, double& yu) const ;
const std::vector<double>& operator[](int i) const ;
// Get data size
int size() const ;
// Get min and max input parameters
float min() const ;
float max() const ;
double min() const ;
double max() const ;
private: // data
// Store for each point of data, the upper
// and lower value
std::vector<std::vector<float> > _data ;
std::vector<std::vector<double> > _data ;
// Store the min and max value on the input
// domain
float _min, _max ;
double _min, _max ;
} ;
class rational_1d_fitter // : public fitting_algorithm
......
......@@ -12,6 +12,7 @@
#include <fstream>
#include <limits>
#include <algorithm>
#include <cmath>
typedef CGAL::MP_Float ET ;
......@@ -93,11 +94,12 @@ bool rational_1d_fitter_cgal::fit_data(const rational_1d_data& data, int np, int
#ifdef DEBUG
std::cout << "<<DEBUG>> delta factor: " << sigma_m << " / " << sigma_M << " = " << delta << std::endl ;
#endif
/* for(int i=0; i<2*data.size(); ++i)
//*
for(int i=0; i<2*data.size(); ++i)
{
qp.set_b(i, delta * ci(i)) ;
}
*/
//*/
#ifdef DEBUG
// Export some informations on the problem to solve
std::cout << "<<DEBUG>> " << qp.get_n() << " variables" << std::endl ;
......@@ -153,13 +155,20 @@ bool rational_1d_fitter_cgal::fit_data(const rational_1d_data& data, int np, int
Solution s = CGAL::solve_nonnegative_quadratic_program(qp, ET());
//*/
if(s.solves_quadratic_program(qp))
bool solves_qp = 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 data
std::vector<float> p, q;
std::vector<double> p, q;
for(int i=0; i<np+nq; ++i)
{
const float v = (float)CGAL::to_double(*(s.variable_numerators_begin()+i)) ;
const double v = CGAL::to_double(*(s.variable_numerators_begin()+i)) ;
if(i < np)
{
......
......@@ -35,8 +35,8 @@ bool rational_1d_fitter_eigen::fit_data(const rational_1d_data& data, int np, in
for(int i=0; i<data.size(); ++i)
{
// Norm of the row vector
float a0_norm = 0.0f ;
float a1_norm = 0.0f ;
double a0_norm = 0.0f ;
double a1_norm = 0.0f ;
// 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)]
......@@ -47,7 +47,7 @@ bool rational_1d_fitter_eigen::fit_data(const rational_1d_data& data, int np, in
// Filling the p part
if(j<np)
{
const float pi = r.p(data[i][0], j) ;
const double pi = r.p(data[i][0], j) ;
a0_norm += pi*pi ;
a1_norm += pi*pi ;
CI(j, 2*i+0) = pi ;
......@@ -56,7 +56,7 @@ bool rational_1d_fitter_eigen::fit_data(const rational_1d_data& data, int np, in
// Filling the q part
else
{
const float qi = r.q(data[i][0], j-np) ;
const double qi = r.q(data[i][0], j-np) ;
a0_norm += qi*qi * (data[i][1]*data[i][1]) ;
CI(j, 2*i+0) = -data[i][1] * qi ;
......@@ -92,19 +92,25 @@ bool rational_1d_fitter_eigen::fit_data(const rational_1d_data& data, int np, in
Eigen::VectorXd x(np+nq) ;
double cost = solve_quadprog(G, g0, CE, ce, CI, ci, x) ;
bool prog_solved = cost != std::numeric_limits<double>::infinity() ;
bool solves_qp = cost != std::numeric_limits<double>::infinity() ;
// Check the data
for(int i=0; i<np+nq; ++i)
{
solves_qp = solves_qp && !isnan(x(i)) && (x(i) != std::numeric_limits<double>::infinity()) ;
}
#ifdef DEBUG
std::cout << "<<DEBUG>> cost: " << cost << std::endl ;
#endif
if(prog_solved)
if(solves_qp)
{
// Recopy the vector data
std::vector<float> p, q;
std::vector<double> p, q;
for(int i=0; i<np+nq; ++i)
{
const float v = x(i) ;
const double v = x(i) ;
if(i < np)
{
......
......@@ -11,7 +11,7 @@ int main(int argc, int argv)
const float x = i / (float)10.0f ;
const float y = exp(-10.0 * x*x) * x*x - 0.1 *x*x*x ;
f << x << "\t" << y << "\t" << 0.01f << std::endl ;
f << x << "\t" << y << "\t" << 1.0f << std::endl ;
}
return 0 ;
......
......@@ -39,10 +39,12 @@ int main(int argc, char** argv)
rational_1d_fitter* fitter ;
if(args.is_defined("algorithm") && args["algorithm"] == std::string("eigen"))
{
std::cout << "<<INFO>> using Eigen method" << std::endl ;
fitter = new rational_1d_fitter_eigen() ;
}
else
{
std::cout << "<<INFO>> using CGAL method" << std::endl ;
fitter = new rational_1d_fitter_cgal() ;
}
......@@ -56,8 +58,8 @@ int main(int argc, char** argv)
//*
std::ofstream file(args["output"], std::ios_base::trunc);
const float dt = (data.max() - data.min()) / 100.0f ;
for(float x=data.min(); x<=data.max(); x+=dt)
const double dt = (data.max() - data.min()) / 100.0f ;
for(double x=data.min(); x<=data.max(); x+=dt)
{
file << x << "\t" << r(x) << std::endl ;
}
......
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