Commit 61221e31 authored by Laurent Belcour's avatar Laurent Belcour

Working on the non-linear module.

parent 1aec0bbb
......@@ -3,6 +3,7 @@
#include <iostream>
#include <limits>
#include <string>
#include <fstream>
clustering::clustering(const data* d, const arguments& args)
{
......@@ -23,7 +24,6 @@ clustering::clustering(const data* d, const arguments& args)
assert(indices.size() > 0 && indices.size() < d->dimX());
// Fit the slice into the domain of definition of the data
int compressed_size = d->dimX()-indices.size();
vec compressed_vec = args.get_vec("cluster-slice", d->dimX(), std::numeric_limits<float>::max());
for(int i=0; i<d->dimX(); ++i)
{
......@@ -70,6 +70,23 @@ clustering::clustering(const data* d, const arguments& args)
_data.push_back(e);
}
std::cout << "<<INFO>> clustering left " << _data.size() << " elements" << std::endl;
save("cluster.txt");
}
void clustering::save(const std::string& filename)
{
std::ofstream file(filename.c_str(), std::ios_base::trunc);
file << "#DIM " << _nX << " " << _nY << std::endl;
for(int i=0; i<size(); ++i)
{
for(int j=0; j<_nX+_nY; ++j)
{
file << _data[i][j] << "\t";
}
file << std::endl;
}
file.close();
}
vec clustering::get(int i) const
......
......@@ -13,6 +13,9 @@ class clustering : public data
//! it into a low dimension one.
clustering(const data* d, const arguments& args);
//! \brief the clustering class can save a clustering to a file.
virtual void save(const std::string& filename);
//! \brief the clustering class cannot load from a file. It requires
//! a complete object to be created.
virtual void load(const std::string& filename)
......
......@@ -2,14 +2,14 @@
#include <functional>
#include <string>
#include <fstream>
#include <QtPlugin>
#include "common.h"
#include "args.h"
#include "params.h"
class data ;
#include "data.h"
/*! \brief A representation of an analytical function.
* \ingroup core
......@@ -123,9 +123,29 @@ class function
//! the data object to output the function at the input location only.
virtual void save_gnuplot(const std::string& filename, const data* d,
const arguments& args) const
{
NOT_IMPLEMENTED();
}
{
#ifdef OLD
std::ofstream file(filename.c_str(), std::ios_base::trunc);
for(int i=0; i<d->size(); ++i)
{
vec v = d->get(i) ;
// vec y1 ; y1.assign(d->dimY(), 0.0) ;
// for(int k=0; k<d->dimY(); ++k) { y1[k] = v[d->dimX() + k] ; }
vec y2 = value(v) ;
for(int u=0; u<d->dimX(); ++u)
file << v[u] << "\t" ;
for(int u=0; u<d->dimY(); ++u)
file << y2[u] << "\t" ;
file << std::endl ;
}
file.close();
#else
NOT_IMPLEMENTED();
#endif
}
//! \brief Output the rational function using a C++ function formating.
virtual void save_cpp(const std::string& filename, const arguments& args) const
......
......@@ -570,7 +570,7 @@ void rational_function::save(const std::string& filename) const
file << "#DIM " << _nX << " " << _nY << std::endl ;
file << "#NP " << a.size() / _nY << std::endl ;
file << "#NQ " << b.size() / _nY << std::endl ;
file << "#BASIS LEGENDRE" << std::endl ;
file << "#BASIS LEGENDRE" << std::endl ;
unsigned int np = a.size() / _nY ;
unsigned int nq = b.size() / _nY ;
......
......@@ -35,13 +35,6 @@ void phong_function::load(const std::string& filename)
throw;
}
//! Save the current function to a specific file type
void phong_function::save(const std::string& filename, const arguments& args) const
{
std::cerr << "Not implemented " << __FILE__ << ":" << __LINE__ << std::endl;
throw;
}
//! Number of parameters to this non-linear function
int phong_function::nbParameters() const
{
......@@ -77,8 +70,32 @@ void phong_function::setParameters(const vec& p)
//! parameters.
vec phong_function::parametersJacobian(const vec& x) const
{
std::cerr << "Not implemented " << __FILE__ << ":" << __LINE__ << std::endl;
throw;
vec jac(dimY()*nbParameters());
for(int i=0; i<dimY(); ++i)
{
for(int j=0; j<dimY(); ++j)
{
if(i == j)
{
// df / dk_d
jac[i+0*dimY() + j*nbParameters()] = 1.0;
// df / dk_s
jac[i+1*dimY() + j*nbParameters()] = std::pow(x[0], _N[j]);
// df / dN
jac[i+2*dimY() + j*nbParameters()] = _N[j] * _ks[i] * std::pow(x[0], _N[j]);
}
else
{
jac[i+0*dimY() + j*nbParameters()] = 0.0;
jac[i+1*dimY() + j*nbParameters()] = 0.0;
jac[i+2*dimY() + j*nbParameters()] = 0.0;
}
}
}
return jac;
}
Q_EXPORT_PLUGIN2(phong_function, phong_function)
......@@ -14,7 +14,7 @@
#include <core/common.h>
/*! \brief A phong lobe class. It is provided for testing with the nonlinear
/*! \brief A phong lobe class. It is provided for testing with the nonlinear
* fitting algorithms.
*
* \details
......@@ -27,17 +27,30 @@ class phong_function : public nonlinear_function, public QObject
// Q_OBJECT
Q_INTERFACES(function)
public: // methods
public: // methods
// Overload the function operator
virtual vec operator()(const vec& x) const ;
virtual vec value(const vec& x) const ;
// Overload the function operator
virtual vec operator()(const vec& x) const ;
virtual vec value(const vec& x) const ;
virtual vec value(const vec& x, const vec& p) const
{
// Test input parameters for correct size
assert(p.size() == nbParameters());
vec res(dimY());
for(int i=0; i<dimY(); ++i)
{
const double kd = p[i*3 + 0];
const double ks = p[i*3 + 1];
const double N = p[i*3 + 2];
res[i] = kd + ks * std::pow(x[0], N);
}
return res;
}
//! \brief Load function specific files
virtual void load(const std::string& filename) ;
//! \brief Save the current function to a specific file type
virtual void save(const std::string& filename, const arguments& args) const ;
//! \brief Number of parameters to this non-linear function
virtual int nbParameters() const ;
......@@ -80,6 +93,35 @@ class phong_function : public nonlinear_function, public QObject
_N.resize(_nY) ;
}
protected: // methods
virtual void save(const std::string& filename) const
{
std::ofstream file(filename.c_str(), std::ios_base::trunc);
file << "#DIM " << _nX << " " << _nY << std::endl ;
file << "#kd" << std::endl;
for(int i=0; i<_nY; ++i)
{
file << _kd[i] << std::endl;
}
file << std::endl;
file << "#ks" << std::endl;
for(int i=0; i<_nY; ++i)
{
file << _ks[i] << std::endl;
}
file << std::endl;
file << "#N" << std::endl;
for(int i=0; i<_nY; ++i)
{
file << _N[i] << std::endl;
}
file << std::endl;
}
private: // data
//! \brief The phong lobe data
......
......@@ -20,21 +20,31 @@ fitter* provide_fitter()
struct EigenFunctor
{
EigenFunctor(nonlinear_function* f) : _f(f)
EigenFunctor(nonlinear_function* f, const data* d) : _f(f)
{
}
inline int inputs() const { return _f->dimX(); }
inline int values() const { return _f->dimY(); }
inline int inputs() const { return _f->nbParameters(); }
inline int values() const { return _f->dimY(); }
int operator()(const Eigen::VectorXd& x, Eigen::VectorXd& y) const
{
vec _x(inputs());
for(int i=0; i<inputs(); ++i) _x[i] = x(i);
vec _y = (*_f)(_x);
// Update the parameters vector
vec _p(inputs());
for(int i=0; i<inputs(); ++i)
{
_p[i] = x(i);
}
vec _x = _d->get(0);
vec _di = vec(_f->dimY());
for(int i=0; i<_f->dimY(); ++i)
_di[i] = _x[_f->dimX() + i];
vec _y = (*_f)(_x) - _di;
for(int i=0; i<values(); ++i) y(i) = _y[i];
std::cout << "f(x) = " << y << std::endl;
return 0;
}
......@@ -50,11 +60,14 @@ struct EigenFunctor
fjac(i,j) = _jac[values()*j + i];
}
std::cout << "J = " << fjac << std::endl;
return 0;
}
nonlinear_function* _f;
const data* _d;
};
nonlinear_fitter_eigen::nonlinear_fitter_eigen() : QObject()
......@@ -66,37 +79,47 @@ nonlinear_fitter_eigen::~nonlinear_fitter_eigen()
bool nonlinear_fitter_eigen::fit_data(const data* d, function* fit, const arguments &args)
{
// I need to set the dimension of the resulting function to be equal
// to the dimension of my fitting problem
fit->setDimX(d->dimX()) ;
fit->setDimY(d->dimY()) ;
fit->setMin(d->min()) ;
fit->setMax(d->max()) ;
if(dynamic_cast<nonlinear_function*>(fit) == NULL)
{
// I need to set the dimension of the resulting function to be equal
// to the dimension of my fitting problem
fit->setDimX(d->dimX()) ;
fit->setDimY(d->dimY()) ;
fit->setMin(d->min()) ;
fit->setMax(d->max()) ;
if(dynamic_cast<nonlinear_function*>(fit) == NULL)
{
std::cerr << "<<ERROR>> the function is not a non-linear function" << std::endl;
return false;
}
nonlinear_function* nf = dynamic_cast<nonlinear_function*>(fit);
return false;
}
nonlinear_function* nf = dynamic_cast<nonlinear_function*>(fit);
/* the following starting values provide a rough fit. */
int info;
Eigen::VectorXd x;
x.setConstant(nf->parameters().size(), 1.);
/* the following starting values provide a rough fit. */
int info;
Eigen::VectorXd x;
x.setConstant(nf->nbParameters(), 1.);
EigenFunctor functor(nf);
Eigen::LevenbergMarquardt<EigenFunctor> lm(functor);
info = lm.minimize(x);
vec p(nf->parameters().size());
for(int i=0; i<p.size(); ++i)
{
p[i] = x(i);
}
nf->setParameters(p);
EigenFunctor functor(nf, d);
Eigen::LevenbergMarquardt<EigenFunctor> lm(functor);
info = lm.minimize(x);
if(info == Eigen::LevenbergMarquardtSpace::ImproperInputParameters)
{
std::cerr << "<<ERROR>> incorrect parameters" << std::endl;
return false;
}
vec p(nf->nbParameters());
for(int i=0; i<p.size(); ++i)
{
p[i] = x(i);
}
std::cout << "<<INFO>> found parameters: " << p << std::endl;
nf->setParameters(p);
return true;
return info == 1 ;
}
void nonlinear_fitter_eigen::set_parameters(const arguments& args)
......
TEMPLATE = subdirs
SUBDIRS = \
rational_fitter_cgal \
# rational_fitter_cgal \
rational_fitter_quadprog \
rational_fitter_parallel \
# rational_fitter_quadproge \
......
......@@ -11,7 +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) : _np(np), _nq(nq), CI(0.0, _np+_nq, 0), _boundary_factor(1.0) { }
quadratic_program(int np, int nq, double boundary_factor) : _np(np), _nq(nq), CI(0.0, _np+_nq, 0), _boundary_factor(boundary_factor) { }
//! \brief Remove the already defined constraints
void clear_constraints()
......@@ -163,6 +164,25 @@ class quadratic_program
vec x, yl, yu;
data->get(n, x, yl, yu);
// <boundary conditions>
if(_boundary_factor != 1.0)
{
bool is_boundary = false;
for(int i=0; i<data->dimX(); ++i)
{
is_boundary = is_boundary || (x[i] <= data->min()[i]) || (x[i] >= data->max()[i]);
}
if(is_boundary)
{
vec mean = 0.5*(yl+yu);
yl = mean + _boundary_factor * (yl - mean);
yu = mean + _boundary_factor * (yu - mean);
}
}
// </boundary condition>
vec y = r->value(x);
bool fail_upper = y > yu;
bool fail_lower = y < yl;
......@@ -238,4 +258,7 @@ class quadratic_program
int _np, _nq ;
QuadProgPP::Matrix<double> CI;
// Boundary conditions factor
double _boundary_factor;
};
......@@ -241,6 +241,7 @@ bool rational_fitter_parallel::fit_data(const data* dat, function* fit, const ar
void rational_fitter_parallel::set_parameters(const arguments& args)
{
_boundary_factor = args.get_float("boundary-constraint", 1.0f);
}
......@@ -295,7 +296,7 @@ bool rational_fitter_parallel::fit_data(const vertical_segment* d, int np, int n
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, _boundary_factor);
#ifndef TODO_PUT_IN_METHOD
for(int i=0; i<d->size()/100; ++i)
......
......@@ -43,7 +43,7 @@ class rational_fitter_parallel : public QObject, public fitter
virtual data* provide_data() const ;
virtual function* provide_function() const ;
protected: // data
protected: // methods
// Fitting a data object using np elements in the numerator and nq
// elements in the denominator
......@@ -53,5 +53,10 @@ class rational_fitter_parallel : public QObject, public fitter
//! \brief Create a constraint vector given its index i in the data
// object and the rational function object to fit.
virtual void get_constraint(int i, int np, int nq, int ny, const vertical_segment* data, const rational_function* func, vec& cu, vec& cl);
protected: // data
// Adding a factor to fit more aggresively the boundaries
double _boundary_factor;
} ;
......@@ -82,11 +82,11 @@ bool rational_fitter_quadprog::fit_data(const data* dat, function* fit, const ar
std::cout << "<<INFO>> fit using np = " << temp_np << " & nq = " << temp_nq << " failed\r" ;
std::cout.flush() ;
if(temp_np <= _max_np)
if(temp_np < _max_np)
{
++temp_np ;
}
if(temp_nq <= _max_nq)
if(temp_nq < _max_nq)
{
++temp_nq ;
}
......@@ -101,6 +101,10 @@ void rational_fitter_quadprog::set_parameters(const arguments& args)
_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) ;
_max_np = std::max<int>(_max_np, _min_np);
_max_nq = std::max<int>(_max_nq, _min_nq);
_boundary = args.get_float("boundary-constraint", 1.0f);
}
......
......@@ -58,7 +58,7 @@ int main(int argc, char** argv)
int sec = (msec / 1000) % 60 ;
int min = (msec / 60000) % 60 ;
int hour = (msec / 3600000) ;
// Display the result
if(is_fitted)
......@@ -81,8 +81,8 @@ int main(int argc, char** argv)
}
/*/
f->save(args["output"], args) ;
#ifndef OLD // use brdf2gnuplot
f->save(args["output"], args) ;
#ifdef OLD // use brdf2gnuplot
size_t n = args["output"].find('.') ;
std::string gnuplot_filename = args["output"].substr(0,n);
gnuplot_filename.append(".gnuplot") ;
......
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