Commit bf1704a9 authored by Laurent Belcour's avatar Laurent Belcour

I can fit a 1D phong

parent 61221e31
......@@ -124,7 +124,7 @@ class function
virtual void save_gnuplot(const std::string& filename, const data* d,
const arguments& args) const
{
#ifdef OLD
#ifndef OLD
std::ofstream file(filename.c_str(), std::ios_base::trunc);
for(int i=0; i<d->size(); ++i)
{
......
......@@ -72,28 +72,19 @@ vec phong_function::parametersJacobian(const vec& x) const
{
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_d
jac[i*nbParameters() + 0] = 1.0;
// df / dk_s
jac[i+1*dimY() + j*nbParameters()] = std::pow(x[0], _N[j]);
// df / dk_s
jac[i*nbParameters() + 1] = std::pow(x[0], _N[0]);
// 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;
}
}
}
// df / dN
if(x[0] == 0.0)
jac[i*nbParameters() + 2] = 0.0;
else
jac[i*nbParameters() + 2] = _ks[0] * std::log(x[0]) * std::pow(x[0], _N[0]);
}
return jac;
}
......
......@@ -2,7 +2,7 @@
#include <Eigen/Dense>
#include <Eigen/SVD>
#include <unsupported/Eigen/NonLinearOptimization>
#include <unsupported/Eigen/LevenbergMarquardt>
#include <string>
#include <iostream>
......@@ -18,49 +18,61 @@ fitter* provide_fitter()
return new nonlinear_fitter_eigen();
}
struct EigenFunctor
struct EigenFunctor: Eigen::DenseFunctor<double>
{
EigenFunctor(nonlinear_function* f, const data* d) : _f(f)
EigenFunctor(nonlinear_function* f, const data* d) :
DenseFunctor<double>(f->nbParameters(), d->size()), _f(f), _d(d)
{
#ifndef DEBUG
std::cout << "<<DEBUG>> constructing an EigenFunctor for n=" << inputs() << " parameters and m=" << values() << " points" << std::endl ;
#endif
}
inline int inputs() const { return _f->nbParameters(); }
inline int values() const { return _f->dimY(); }
int operator()(const Eigen::VectorXd& x, Eigen::VectorXd& y) const
{
// 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;
std::cout << "input vector = " << x << std::endl;
// Update the parameters vector
vec _p(inputs());
for(int i=0; i<inputs(); ++i)
{
_p[i] = x(i);
}
_f->setParameters(_p);
for(int s=0; s<values(); ++s)
{
vec _x = _d->get(s);
vec _di = vec(_f->dimY());
for(int i=0; i<_f->dimY(); ++i)
_di[i] = _x[_f->dimX() + i];
// Should add the resulting vector completely
vec _y = _di - (*_f)(_x);
y(s) = _y[0];
}
return 0;
}
int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const
{
vec _x(inputs());
for(int i=0; i<inputs(); ++i) _x[i] = x(i);
vec _p(inputs());
for(int i=0; i<inputs(); ++i) _p[i] = x(i);
_f->setParameters(_p);
for(int i=0; i<values(); ++i)
{
vec xi = _d->get(i);
vec _jac = _f->parametersJacobian(_x);
for(int j=0; j<_f->nbParameters(); ++j)
for(int i=0; i<values(); ++i)
vec _jac = _f->parametersJacobian(xi);
for(int j=0; j<_f->nbParameters(); ++j)
{
fjac(i,j) = _jac[values()*j + i];
fjac(i,j) = -_jac[j];
}
std::cout << "J = " << fjac << std::endl;
}
return 0;
}
......@@ -93,15 +105,19 @@ bool nonlinear_fitter_eigen::fit_data(const data* d, function* fit, const argume
}
nonlinear_function* nf = dynamic_cast<nonlinear_function*>(fit);
#ifndef DEBUG
std::cout << "<<DEBUG>> number of parameters: " << nf->nbParameters() << std::endl;
#endif
/* the following starting values provide a rough fit. */
int info;
Eigen::VectorXd x;
x.setConstant(nf->nbParameters(), 1.);
EigenFunctor functor(nf, d);
Eigen::LevenbergMarquardt<EigenFunctor> lm(functor);
info = lm.minimize(x);
info = lm.minimize(x);
if(info == Eigen::LevenbergMarquardtSpace::ImproperInputParameters)
{
......@@ -117,6 +133,9 @@ bool nonlinear_fitter_eigen::fit_data(const data* d, function* fit, const argume
p[i] = x(i);
}
std::cout << "<<INFO>> found parameters: " << p << std::endl;
#ifndef DEBUG
std::cout << "<<DEBUG>> using " << lm.iterations() << " iterations" << std::endl;
#endif
nf->setParameters(p);
return true;
......
......@@ -73,7 +73,19 @@ int main(int argc, char** argv)
for(int i=0; i<nbx; ++i)
{
const float x = i / (float)nbx ;
const float z = exp(-10.0 * x*x) + 0.1 ;
const float z = 0.1 + 0.5 * std::pow(x, 1.5) ;
f << x << "\t" << z << std::endl ;
}
}
else if(k == 6)
{
f << "#DIM 1 1" << std::endl ;
f << "#PARAM_IN COS_TH" << std::endl;
for(int i=0; i<nbx; ++i)
{
const float x = i / (float)nbx ;
const float z = 1.0 - x ;
f << x << "\t" << z << std::endl ;
}
......
......@@ -82,7 +82,7 @@ int main(int argc, char** argv)
/*/
f->save(args["output"], args) ;
#ifdef OLD // use brdf2gnuplot
#ifndef 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