Commit c0951e87 authored by Laurent Belcour's avatar Laurent Belcour

Adding a cleaner interface to the rational_1d fitter example

parent f9b35cb9
......@@ -60,17 +60,17 @@ class arguments
return std::string() ;
}
} ;
float get_float(const std::string& key, float default_value = 0.0f)
float get_float(const std::string& key, float default_value = 0.0f) const
{
std::string value = _map[key] ;
std::string value = _map.at(key) ;
if(value.empty())
return default_value ;
else
return atof(value.c_str()) ;
} ;
int get_int(const std::string& key, int default_value = 0)
int get_int(const std::string& key, int default_value = 0) const
{
std::string value = _map[key] ;
std::string value = _map.at(key) ;
if(value.empty())
return default_value ;
else
......
......@@ -191,10 +191,36 @@ double rational_1d_data::max() const
return _max ;
}
void rational_1d_fitter::set_parameters(int min_np, int max_np, int min_nq, int max_nq)
bool rational_1d_fitter::fit_data(const rational_1d_data& data, rational_1d& fit)
{
_min_np = min_np ;
_min_nq = min_nq ;
_max_np = max_np ;
_max_nq = max_nq ;
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(data, 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::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) ;
}
......@@ -10,6 +10,7 @@
#include <core/function.h>
#include <core/data.h>
#include <core/fitter.h>
#include <core/args.h>
class rational_1d : public function
{
......@@ -76,13 +77,13 @@ class rational_1d_fitter //: public fitter
public: // methods
// Fitting a data object
virtual bool fit_data(const rational_1d_data& data, rational_1d& fit) = 0;
virtual bool fit_data(const rational_1d_data& data, rational_1d& fit) ;
// Fitting a data object using np elements in the numerator and nq
// elements in the denominator
virtual bool fit_data(const rational_1d_data& data, int np, int nq, rational_1d& fit) = 0;
virtual bool fit_data(const rational_1d_data& data, int np, int nq, rational_1d& fit) = 0 ;
virtual void set_parameters(int min_np, int max_np, int min_nq, int max_nq) ;
virtual void set_parameters(const arguments& args) ;
protected: // data
......
......@@ -20,32 +20,6 @@ typedef CGAL::Quadratic_program<ET> Program ;
typedef CGAL::Quadratic_program_solution<ET> Solution ;
typedef CGAL::Quadratic_program_options Options ;
// Fitting a data
bool rational_1d_fitter_cgal::fit_data(const rational_1d_data& data, rational_1d& 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(data, 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 ;
}
bool rational_1d_fitter_cgal::fit_data(const rational_1d_data& data, int np, int nq, rational_1d& r)
{
......@@ -93,7 +67,7 @@ bool rational_1d_fitter_cgal::fit_data(const rational_1d_data& data, int np, int
const double qi = r.q(data[i][0], j-np) ;
a0_norm += qi*qi * (data[i][1]*data[i][1]) ;
qp.set_a(j, i, ET(-data[i][1] * qi)) ;
CI(j, i+data.size()) = -data[i][1] * qi ;
CI(j, i) = -data[i][1] * qi ;
a1_norm += qi*qi * (data[i][2]*data[i][2]) ;
qp.set_a(j, i+data.size(), ET(data[i][2] * qi)) ;
......
......@@ -14,9 +14,6 @@ class rational_1d_fitter_cgal : public rational_1d_fitter
{
public: // methods
// Fitting a data object
virtual bool fit_data(const rational_1d_data& data, rational_1d& fit) ;
// Fitting a data object using np elements in the numerator and nq
// elements in the denominator.
virtual bool fit_data(const rational_1d_data& data, int np, int nq, rational_1d& fit) ;
......
......@@ -7,12 +7,6 @@
#include <iostream>
#include <cfenv>
// Fitting a data
bool rational_1d_fitter_eigen::fit_data(const rational_1d_data& data, rational_1d& fit)
{
return fit_data(data, 10, 10, fit) ;
}
bool rational_1d_fitter_eigen::fit_data(const rational_1d_data& data, int np, int nq, rational_1d& r)
{
// Select the size of the result vector to be equal to the dimension
......@@ -35,8 +29,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
double a0_norm = 0.0f ;
double a1_norm = 0.0f ;
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)]
......@@ -50,25 +44,25 @@ bool rational_1d_fitter_eigen::fit_data(const rational_1d_data& data, int np, in
const double pi = r.p(data[i][0], j) ;
a0_norm += pi*pi ;
a1_norm += pi*pi ;
CI(j, 2*i+0) = pi ;
CI(j, 2*i+1) = - pi ;
CI(j, i) = pi ;
CI(j, i+data.size()) = - pi ;
}
// Filling the q part
else
{
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 ;
CI(j, i) = - data[i][1] * qi ;
a1_norm += qi*qi * (data[i][2]*data[i][2]) ;
CI(j, 2*i+1) = data[i][2] * qi ;
CI(j, i+data.size()) = data[i][2] * qi ;
}
}
// Set the c vector, will later be updated using the
// delta parameter.
ci(2*i+0) = -sqrt(a0_norm) ;
ci(2*i+1) = -sqrt(a1_norm) ;
ci(i) = -sqrt(a0_norm) ;
ci(i+data.size()) = -sqrt(a1_norm) ;
}
// Update the ci column with the delta parameter
......@@ -76,7 +70,25 @@ bool rational_1d_fitter_eigen::fit_data(const rational_1d_data& data, int np, in
Eigen::JacobiSVD<Eigen::MatrixXd> svd(CI);
const double sigma_m = svd.singularValues()(std::min(2*data.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*data.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
......@@ -92,7 +104,7 @@ 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 solves_qp = cost != std::numeric_limits<double>::infinity() ;
bool solves_qp = abs(cost) != std::numeric_limits<double>::infinity() /*&& cost > 0*/ ;
// Check the data
for(int i=0; i<np+nq; ++i)
......
......@@ -6,9 +6,6 @@ class rational_1d_fitter_eigen : public rational_1d_fitter
{
public: // methods
// Fitting a data object
virtual bool fit_data(const rational_1d_data& data, rational_1d& fit) ;
// Fitting a data object using np elements
// in the numerator and nq elements in the
// denominator
......
......@@ -12,7 +12,13 @@ int main(int argc, char** argv)
{
arguments args(argc, argv) ;
if(args.is_defined("help")) {
std::cout << argv[0] << " --np <int> --nq <int> --input <filename> --output <filename>" << std::endl ;
std::cout << argv[0] << " --algorithm {cgal|eigen} --min <float> --max <float> --min-np <int> --min-nq <int> --np <int> --nq <int> --input <filename> --output <filename>" << std::endl << std::endl ;
std::cout << " ....\"min\" defines the left boundary of the domain" << std::endl ;
std::cout << " ....\"max\" defines the right boundary of the domain" << std::endl ;
std::cout << " .\"min-np\" defines the minimum number of elements at the numerator" << std::endl ;
std::cout << " .....\"np\" defines the maxmimum number of elements at the numerator" << std::endl ;
std::cout << " .\"min-nq\" defines the minimum number of elements at the denominator" << std::endl ;
std::cout << " .....\"nq\" defines the maxmimum number of elements at the denominator" << std::endl ;
return 0 ;
}
......@@ -56,7 +62,7 @@ int main(int argc, char** argv)
std::cout << "<<INFO>> using CGAL method" << std::endl ;
fitter = new rational_1d_fitter_cgal() ;
}
fitter->set_parameters(args.get_int("min-np", 10), args.get_int("np", 10), args.get_int("min-nq", 10), args.get_int("nq", 10)) ;
fitter->set_parameters(args) ;
rational_1d r ;
bool is_fitted = fitter->fit_data(data, r) ;
......
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