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

Making the 3D Yellow data working.

parent c3eb5492
......@@ -7,8 +7,7 @@ CONFIG -= qt
DESTDIR = ../build
unix{
# QMAKE_CXXFLAGS += -std=c++0x
QMAKE_CXXFLAGS += -Wall -pedantic
QMAKE_CXXFLAGS += -Wall -pedantic -fPIC
LIBS += -ldl
}
......
......@@ -30,6 +30,69 @@ bool rational_function_1d::load(std::istream&)
return true;
}
void rational_function_1d::save_body(std::ostream& out, const arguments& args) const
{
bool is_matlab = args["export"] == "matlab";
bool is_alta = !args.is_defined("export") || args["export"] == "alta";
if(is_alta)
{
for(int i=0; i<a.size(); ++i)
{
std::vector<int> index = this->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
out << index[j] << "\t" ;
}
out << a[i] << std::endl ;
}
for(int i=0; i<b.size(); ++i)
{
std::vector<int> index = this->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
out << index[j] << "\t" ;
}
out << b[i] << std::endl ;
}
}
else if(is_matlab)
{
out << "(";
for(int i=0; i<a.size(); ++i)
{
out << a[i];
std::vector<int> indices = this->index2degree(i);
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << indices[k];
}
if(i != a.size()-1) { out << " + "; }
}
out << ") / (";
for(int i=0; i<b.size(); ++i)
{
out << b[i] << "x.^" << i;
std::vector<int> indices = this->index2degree(i);
for(int k=0; k<dimX(); ++k)
{
if(k != dimX()-1) { out << ".*"; }
out << "x(k).^" << indices[k];
}
if(i != b.size()-1) { out << " + "; }
}
out << ")";
}
else
{
NOT_IMPLEMENTED();
}
}
void rational_function_1d::update(const vec& in_a,
const vec& in_b)
{
......@@ -472,7 +535,7 @@ bool rational_function::load(std::istream& in)
}
void rational_function::save_call(std::ostream& out, const arguments&) const
void rational_function::save_call(std::ostream& out, const arguments& args) const
{
out << "#FUNC rational_function" << std::endl;
out << "#NP " << np << std::endl ;
......@@ -483,28 +546,7 @@ void rational_function::save_call(std::ostream& out, const arguments&) const
for(int k=0; k<_nY; ++k)
{
const rational_function_1d* rf = get(k);
vec a = rf->getP();
vec b = rf->getQ();
for(int i=0; i<np; ++i)
{
std::vector<int> index = rf->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
out << index[j] << "\t" ;
}
out << a[i] << std::endl ;
}
for(int i=0; i<nq; ++i)
{
std::vector<int> index = rf->index2degree(i) ;
for(unsigned int j=0; j<index.size(); ++j)
{
out << index[j] << "\t" ;
}
out << b[i] << std::endl ;
}
rf->save_body(out, args);
}
}
......
......@@ -25,10 +25,25 @@ class rational_function_1d : public function
rational_function_1d(const vec& a, const vec& b) ;
virtual ~rational_function_1d() {}
/* FUNCTION INHERITANCE */
// Overload the function operator
virtual vec value(const vec& x) const ;
virtual vec operator()(const vec& x) const { return value(x) ; }
// IO function to text files
virtual bool load(std::istream& in);
//! \brief Save the rational function expansion. It should
//! not be store in any variable (e.g. "y = rf(x);") as the
//! nD rational function can append factor to the 1D rational
//! function.
virtual void save_body(std::ostream&, const arguments&) const;
/* RATIONAL FUNCTION SPECIFIC */
// Get the numerator (p) and denominator (q) functions
virtual vec p(const vec& x) const ;
virtual vec q(const vec& x) const ;
......@@ -37,9 +52,6 @@ class rational_function_1d : public function
virtual double p(const vec& x, int i) const ;
virtual double q(const vec& x, int j) const ;
// IO function to text files
virtual bool load(std::istream& in);
// Update the function
virtual void update(const vec& in_a,
const vec& in_b) ;
......@@ -58,7 +70,6 @@ class rational_function_1d : public function
friend std::ostream& operator<< (std::ostream& out,
const rational_function_1d& r) ;
//! Convert a 1D index into a vector of degree for a
//! multinomial coeffcient. The resulting vector v should
//! be used as prod_k x[k]^v[k] for the monomial basis
......@@ -67,6 +78,7 @@ class rational_function_1d : public function
protected: // functions
static int estimate_dk(int k, int d);
static void populate(std::vector<int>& vec, int N, int M, int j);
......@@ -76,10 +88,10 @@ class rational_function_1d : public function
// the functions to be polynomials.
vec a, b ;
//! Is the function separable with respect to its input dimensions?
//! \todo Make possible to have only part of the dimensions
//! separable.
bool _separable;
//! Is the function separable with respect to its input dimensions?
//! \todo Make possible to have only part of the dimensions
//! separable.
bool _separable;
} ;
class rational_function : public function
......
......@@ -154,11 +154,14 @@ class quadratic_program
return solves_qp;
}
#define PACANOWSKI2012
//! \brief Test all the constraints of the data.
//! Add the sample that is farest away from the function.
bool test_constraints(int ny, const rational_function_1d* r, const vertical_segment* data)
{
int nb_failed = 0;
#ifdef PACANOWSKI2012
int nb_failed = 0;
double max_dev = 0.0; // Maximum absolute distance of the current
// solution to the data.
vec cu, cl;
......@@ -200,6 +203,26 @@ class quadratic_program
{
return true;
}
#else
int n = next_unmatching_constraint(0, ny, r, data);
if(n < data->size())
{
vec x, yl, yu;
data->get(n, x, yl, yu);
vec cu, cl;
get_constraint(x, yl, yu, ny, r, cu, cl);
add_constraints(cu);
add_constraints(cl);
return false;
}
else
{
return true;
}
#endif
}
//! \brief Generate two constraint vectors from a vertical segment and a
......@@ -257,7 +280,7 @@ int quadratic_program::next_unmatching_constraint(int i, int ny, const rational_
data->get(n, x, yl, yu);
vec y = r->value(x);
if(y[ny] < yl[ny] || y[ny] > yu[ny])
if(y[ny] < yl[ny] || y[ny] > yu[ny])
{
return n;
}
......
......@@ -21,7 +21,7 @@ ALTA_DLL_EXPORT fitter* provide_fitter()
return new rational_fitter_parallel();
}
rational_fitter_parallel::rational_fitter_parallel()
rational_fitter_parallel::rational_fitter_parallel() : nb_starting_points(100)
{
}
rational_fitter_parallel::~rational_fitter_parallel()
......@@ -49,6 +49,9 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
const int _max_np = args.get_int("np", _min_np);
std::cout << "<<INFO>> N in [" << _min_np << ", " << _max_np << "]" << std::endl ;
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;
for(int i=_min_np; i<=_max_np; ++i)
{
std::cout << "<<INFO>> fit using np+nq = " << i << std::endl ;
......@@ -56,153 +59,90 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
timer time ;
time.start() ;
#ifdef NOT_WORKING
// Allocate enough processor to run fully in parallel, but account for
// the fact that each thread will require its own memory.
size_t need_memory = ((3*i+1)*d->size()*2 + 2*i + 2*d->dimY()*d->size()) * sizeof(double);
size_t avai_memory = plugins_manager::get_system_memory();
int nb_cores = (60 * avai_memory) / (100 * need_memory) ;
nb_cores = std::min<int>(nb_cores, omp_get_num_procs());
if(nb_cores == 0)
{
std::cerr << "<<ERROR>> not enough memory to perform the fit" << std::endl ;
#ifdef DEBUG
std::cout << "<<DEBUG>> " << need_memory / (1024*1024) << "MB required / "
<< avai_memory / (1024*1024) << "MB available" << std::endl;
#endif
//return false;
nb_cores = 1;
}
#endif
int nb_cores = args.get_int("nbcores", omp_get_num_procs());
int nb_cores = args.get_int("nb-cores", omp_get_num_procs());
#ifdef DEBUG
std::cout << "<<DEBUG>> will use " << nb_cores << " threads to compute the quadratic programs" << std::endl ;
#endif
omp_set_num_threads(nb_cores) ;
std::vector<rational_function*> rs;
for(int j=0; j<nb_cores; ++j)
{
rational_function* rj = dynamic_cast<rational_function*>(plugins_manager::get_function(args));
//! \todo I need to do a check of compatibility here, but the
//! signature of the check_compatibility function is not good.
rj->setParametrization(fit->input_parametrization());
rj->setParametrization(fit->output_parametrization());
rj->setDimX(d->dimX()) ;
rj->setDimY(d->dimY()) ;
rj->setMin(d->min()) ;
rj->setMax(d->max()) ;
if(rj == NULL)
{
std::cerr << "<<ERROR>> unable to obtain a rational function from the plugins manager" << std::endl;
return false;
}
rs.push_back(rj);
}
// Solution for the case of optimizing the L2 norm
double min_l2 = std::numeric_limits<double>::max();
rational_function* min_l2_fun = NULL;
// Solution for the case of optimizing the LINF norm
double min_linf = std::numeric_limits<double>::max();
rational_function* min_linf_fun = NULL;
double min_delta = std::numeric_limits<double>::max();
double mean_delta = 0.0;
int nb_sol_found = 0;
int nb_sol_tested = 0;
#pragma omp parallel for
for(int j=1; j<i; ++j)
{
int temp_np = i - j;
int temp_nq = j;
#pragma omp parallel for shared(nb_sol_found, nb_sol_tested, min_delta, mean_delta), schedule(dynamic,1)
for(int j=1; j<i; ++j)
{
// Compute the number of coefficients in the numerator and in the denominator
// from the current number of coefficients i and the current index in the
// loop j.
int temp_np = i - j;
int temp_nq = j;
vec p(temp_np*r->dimY()), q(temp_nq*r->dimY());
// Allocate a rational function and set it to the correct size, dimensions
// and parametrizations.
rational_function* rk = dynamic_cast<rational_function*>(plugins_manager::get_function(args));
rk->setParametrization(r->input_parametrization());
rk->setParametrization(r->output_parametrization());
rk->setDimX(r->dimX()) ;
rk->setDimY(r->dimY()) ;
rk->setMin(r->min()) ;
rk->setMax(r->max()) ;
if(rk == NULL)
{
std::cerr << "<<ERROR>> unable to obtain a rational function from the plugins manager" << std::endl;
throw;
}
vec p(temp_np*r->dimY()), q(temp_nq*r->dimY());
rs[omp_get_thread_num()]->setSize(temp_np, temp_nq);
// 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);
if(is_fitted)
{
#pragma omp critical (nb_sol_found)
{
++nb_sol_found ;
mean_delta += delta ;
std::cout << "<<INFO>> found a solution with np=" << temp_np << ", nq = " << temp_nq << std::endl;
std::cout << "<<INFO>> Linf error = " << linf_dist << std::endl;
std::cout << "<<INFO>> L2 error = " << l2_dist << std::endl;
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)
{
min_delta = delta ;
r->setSize(temp_np, temp_nq);
for(int y=0; y<r->dimY(); ++y)
{
r->update(y, rk->get(y));
}
}
}
}
double delta, linf_dist, l2_dist;
bool is_fitted = fit_data(d, temp_np, temp_nq, rs[omp_get_thread_num()], args, p, q, delta, linf_dist, l2_dist);
if(is_fitted)
{
#pragma omp critical
{
++nb_sol_found ;
mean_delta += delta ;
std::cout << "<<INFO>> found a solution with np=" << temp_np << ", nq = " << temp_nq << std::endl;
std::cout << "<<INFO>> Linf error = " << linf_dist << std::endl;
std::cout << "<<INFO>> L2 error = " << l2_dist << std::endl;
std::cout << "<<INFO>> delta = " << delta << std::endl;
std::cout << std::endl;
// Get the solution with the minimum delta
if(delta < min_delta)
{
min_delta = delta ;
r->setSize(temp_np, temp_nq);
for(int y=0; y<r->dimY(); ++y)
{
r->update(y, rs[omp_get_thread_num()]->get(y));
}
}
// Get the solution with the minumum L2 norm
if(l2_dist < min_l2)
{
min_l2 = l2_dist;
min_l2_fun = rs[omp_get_thread_num()];
}
// Get the solution with the minumum LINF norm
if(linf_dist < min_linf)
{
min_linf = linf_dist;
min_linf_fun = rs[omp_get_thread_num()];
}
}
}
#pragma omp critical
{
// Update the solution
nb_sol_tested++;
std::cout << "<<DEBUG>> nb solutions tested: " << nb_sol_tested << " / " << i << "\r";
std::cout.flush();
}
}
delete rk; // memory clean
/*
// Clean memory
for(int j=0; j<nb_cores; ++j)
{
delete rs[j];
}
*/
std::cout << " \r";
if(min_l2 < std::numeric_limits<double>::max())
{
std::cout << "<<INFO>> min L2 = " << min_l2 << std::endl;
std::cout << *min_l2_fun << std::endl;
std::cout << std::endl;
}
// delete min_l2_fun;
#pragma omp critical (nb_sol_tested)
{
// Update the solution
nb_sol_tested++;
if(min_linf < std::numeric_limits<double>::max())
{
std::cout << "<<INFO>> min Linf = " << min_linf << std::endl;
std::cout << *min_linf_fun << std::endl;
std::cout << std::endl;
}
// delete min_linf_fun;
std::cout << "<<DEBUG>> nb solutions tested: " << nb_sol_tested << " / " << i << "\r";
std::cout.flush();
}
}
if(min_delta < std::numeric_limits<double>::max())
if(min_delta < std::numeric_limits<double>::max())
{
std::cout << "<<INFO>> mean delta = " << mean_delta/nb_sol_found << std::endl;
std::cout << "<<INFO>> min delta = " << min_delta << std::endl;
......@@ -261,8 +201,8 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n
quadratic_program qp(np, nq);
// Starting with only a 100 vertical segments (meaning a 200xn matrix)
const int di = (m-1) / 99;
// Starting with only a nb_starting_points vertical segments
const int di = (m-1) / (nb_starting_points-1);
for(int i=0; i<m; i+=di)
{
......@@ -274,7 +214,7 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n
qp.add_constraints(c2);
}
while(qp.nb_constraints() < 2*m)
while(qp.nb_constraints() < 2*m)
{
#ifdef DEBUG
std::cout << "<<DEBUG>> thread " << omp_get_thread_num() << ", number of intervals tested = " << qp.nb_constraints()/2 << std::endl ;
......@@ -290,24 +230,6 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n
#ifdef DEBUG
std::cout << "<<INFO>> got solution " << *r << std::endl ;
#endif
/*
int current = 0, i=0;
while(i < 100 && current < m)
{
int next = quadratic_program::next_unmatching_constraint(current, ny, );
// Create two vector of constraints
vec c1(n), c2(n);
get_constraint(next, np, nq, ny, d, r, c1, c2);
qp.add_constraints(c1);
qp.add_constraints(c2);
++i;
current = next;
}
*/
return true;
}
}
......@@ -320,7 +242,6 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n
}
}
return false;
}
......
......@@ -58,5 +58,8 @@ class rational_fitter_parallel : public fitter
protected: // data
//! Number of points used when starting the adaptive interpolation
//! procedure. By default, this value is 100.
int nb_starting_points;
} ;
......@@ -19,4 +19,5 @@ LIBS += -L../../build \
-lcore
#QMAKE_CXXFLAGS += -fPIC
QMAKE_CXXFLAGS += -g
......@@ -18,4 +18,4 @@ plot "../papers/retro/mesures/original/3M_jaune/3d/633nm/Fichiers definitifs/den
# output rational fits
set output "yellow_retro_rat.tex"
plot "../papers/retro/mesures/original/3M_jaune/3d/633nm/Fichiers definitifs/densify_helmholtz/3M_jaune_3D+3DS+3DR__BRDF_min_retro_lobe_dense.alta" using 2:($3 > 0.0 && $3 < 0.005 ? $4 : 1/0) title "Yellow cloth data", "./results/3d/retro/half/3M_jaune_rat.dat" using 1:($2 > 0.0 && $2 < 0.005 ? $3 : 1/0) title "rational interpolation"
plot "../papers/retro/mesures/original/3M_jaune/3d/633nm/Fichiers definitifs/densify_helmholtz/3M_jaune_3D+3DS+3DR__BRDF_min_retro_lobe_dense.alta" using 2:($3 > 0.0 && $3 < 0.005 ? $4 : 1/0) title "Yellow cloth data", "./results/3d/retro/half/3M_jaune_rat.dat" using 2:($3 > -0.01 && $3 < 0.05 ? $3 : 1/0) title "rational interpolation"
......@@ -41,13 +41,14 @@
<!-- Define the ftting procedure to use -->
<!--<plugin type="fitter" name="rational_fitter_quadprog" />-->
<plugin type="fitter" name="rational_fitter_parallel" />
<parameter name="min-np" value="100" />
<parameter name="min-np" value="60" />
<!--<parameter name="min-np" value="53" />-->
<parameter name="min-nq" value="10" />
<parameter name="np" value="100" />
<parameter name="nq" value="100" />
<parameter name="dt" value="0.5" />
<parameter name="data-positive" value="" />
<parameter name="nb-starting-points" value="2000" />
<!--<parameter name="dt-relative" value="" />-->
</action>
......@@ -58,7 +59,7 @@
<input name="./results/3d/retro/half/3M_jaune_rat.brdf" />
<output name="./results/3d/retro/half/3M_jaune_rat.dat" />
<parameter name="data" value="../papers/retro/mesures/original/3M_jaune/3d/633nm/Fichiers\ definitifs/densify_helmholtz/3M_jaune_3D+3DS+3DR_dense__nbsgrid_81.alta" />
<parameter name="data" value="../papers/retro/mesures/original/3M_jaune/3d/633nm/Fichiers\ definitifs/densify_helmholtz/3M_jaune_3D+3DS+3DR_dense__nbsgrid_162.alta" />
<!--<parameter name="data" value="/tmp/yellow_slice_inc30.dat" />-->
</action>
</alta>
<?xml version="1.0"?>
<alta>
<!-- This script file compute the fitting of the retro-reflecting materials
present in the ALTA library. This script file should be executed in
the sources directory of repository as all directories are relative.
It is also necessary to create a results/3d/retro directory to store
the resulting fits and exports to BRDF-explorer and matlab.