Commit 0c877e68 authored by Laurent Belcour's avatar Laurent Belcour

Cleaning of the parallel rational fitter.

Added an interface tool to select if it has to compute the delta or not.
parent b63dbe6d
......@@ -11,8 +11,8 @@ 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)
quadratic_program(int np, int nq, bool compute_delta = false) :
_np(np), _nq(nq), _compute_delta(compute_delta), CI(0.0, _np+_nq, 0)
{ }
//! \brief Remove the already defined constraints
......@@ -122,16 +122,17 @@ class quadratic_program
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) ;
delta = sigma_M / sigma_m ;
/*/
delta = 1.0;
//*/
if(_compute_delta)
{
// 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) ;
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)
......@@ -255,6 +256,7 @@ class quadratic_program
protected:
int _np, _nq;
bool _compute_delta;
QuadProgPP::Matrix<double> CI;
//! Contains the indices of the vertical segment unused during the
......
......@@ -56,7 +56,8 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
const int nb_starting_points = args.get_int("nb-starting-points", 100);
std::cout << "<<INFO>> number of data point used in start: " << nb_starting_points << std::endl;
const int step = args.get_int("np-step", 1);
const int step = args.get_int("np-step", 1);
const bool use_delta = args.is_defined("use_delta");
for(int i=_min_np; i<=_max_np; i+=step)
{
......@@ -74,7 +75,8 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
omp_set_num_threads(nb_cores) ;
#endif
double min_delta = std::numeric_limits<double>::max();
double min_delta = std::numeric_limits<double>::max();
double min_l2_dist = std::numeric_limits<double>::max();
double mean_delta = 0.0;
int nb_sol_found = 0;
int nb_sol_tested = 0;
......@@ -112,8 +114,9 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
// Set the rational function size
rk->setSize(temp_np, temp_nq);
double delta, linf_dist, l2_dist;
bool is_fitted = fit_data(d, temp_np, temp_nq, rk, args, p, q, delta, linf_dist, l2_dist);
double delta = 1.0;
double linf_dist, l2_dist;
bool is_fitted = fit_data(d, temp_np, temp_nq, rk, args, delta, linf_dist, l2_dist);
if(is_fitted)
{
#pragma omp critical (nb_sol_found)
......@@ -127,11 +130,12 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
std::cout << "<<INFO>> delta = " << delta << std::endl;
std::cout << std::endl;
// Get the solution with the minimum delta, and update the main
// rational function r.
if(delta < min_delta)
// Get the solution with the minimum delta or the minimum L2 distance,
// and update the main rational function r.
if((use_delta && delta < min_delta) || (!use_delta && l2_dist < min_l2_dist))
{
min_delta = delta ;
min_delta = delta ;
min_l2_dist = l2_dist ;
r->setSize(temp_np, temp_nq);
for(int y=0; y<r->dimY(); ++y)
{
......@@ -171,14 +175,14 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
return false ;
}
void rational_fitter_parallel::set_parameters(const arguments& args)
void rational_fitter_parallel::set_parameters(const arguments&)
{
}
bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int nq,
rational_function* r, const arguments &args,
vec& P, vec& Q, double& delta, double& linf_dist, double& l2_dist)
rational_function* r, const arguments &args,
double& delta, double& linf_dist, double& l2_dist)
{
// Fit the different output dimension independantly
for(int j=0; j<d->dimY(); ++j)
......@@ -186,7 +190,8 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n
vec p(np), q(nq);
rational_function_1d* rf = r->get(j);
rf->resize(np, nq);
if(!fit_data(d, np, nq, j, rf, p, q, delta))
if(!fit_data(d, np, nq, j, rf, args, p, q, delta))
{
return false ;
}
......@@ -205,12 +210,13 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n
// y is the dimension to fit on the y-data (e.g. R, G or B for RGB signals)
// the function returns a rational BRDF function and a boolean
bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int nq, int ny,
rational_function_1d* r, vec& p, vec& q, double& delta)
rational_function_1d* r, const arguments& args,
vec& p, vec& q, double& delta)
{
const int m = d->size(); // 2*m = number of constraints
const int n = np+nq; // n = np+nq
quadratic_program qp(np, nq);
quadratic_program qp(np, nq, args.is_defined("use_delta"));
// Starting with only a nb_starting_points vertical segments
std::list<unsigned int> training_set;
......
......@@ -12,12 +12,13 @@
#include <core/fitter.h>
#include <core/args.h>
/*! \brief A vertical segment fitter for rational functions that can run in parallel
* using QuadProg++ quadratic solver.
/*! \brief A vertical segment fitter for rational functions that search for a solution
* for a fixed number of coefficient. This plugin can run in parallel with OpenMP and
* is using QuadProg++ quadratic solver.
* \ingroup plugins
*
* \details
* You can find QuadProg++ here: http://quadprog.sourceforge.net/
* You can find QuadProg++ <a href="http://quadprog.sourceforge.net/">here</a>.
* <br />
*
* <h3>Plugin parameters</h3>
......@@ -33,6 +34,9 @@
* of the rational function. By default, this number is 1.</li>
* <li><b>--nb-cores</b> <em>[int]</em> number of core allocated to perform
* the seach. By default, this is equal to the number of processors.</li>
* <li><b>--use_delta</b> use the strategy of Pacanowski et al. [2012] to
* modify the constraint vector by the condition number of the constraint
* matrix. We did not experience any benefit from using it.</li>
* </ul>
*/
class rational_fitter_parallel : public fitter
......@@ -56,11 +60,10 @@ class rational_fitter_parallel : public fitter
// elements in the denominator
virtual bool fit_data(const vertical_segment* d, int np, int nq,
rational_function* fit, const arguments &args,
vec& p, vec& q, double& delta,
double& linf_dist,double& l2_dist) ;
double& delta, double& linf_dist,double& l2_dist) ;
virtual bool fit_data(const vertical_segment* dat, int np, int nq,
int ny, rational_function_1d* fit, vec& p, vec& q,
double& delta) ;
int ny, rational_function_1d* fit, const arguments& args,
vec& p, vec& q, double& delta) ;
//! \brief Create a constraint vector given its index i in the data
//! object and the rational function object to fit. This function
......
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