Mise à jour terminée. Pour connaître les apports de la version 13.8.4 par rapport à notre ancienne version vous pouvez lire les "Release Notes" suivantes :
https://about.gitlab.com/releases/2021/02/11/security-release-gitlab-13-8-4-released/
https://about.gitlab.com/releases/2021/02/05/gitlab-13-8-3-released/

Commit 04fed075 authored by Laurent Belcour's avatar Laurent Belcour

Working Levenberg-Marquardt algorithm for multidimensional data

parent bf1704a9
......@@ -47,9 +47,9 @@ vec phong_function::parameters() const
vec res(3*dimY());
for(int i=0; i<dimY(); ++i)
{
res[i] = _kd[i];
res[i+ dimY()] = _ks[i];
res[i+2*dimY()] = _N[i];
res[i*3 + 0] = _kd[i];
res[i*3 + 1] = _ks[i];
res[i*3 + 2] = _N[i];
}
return res;
......@@ -60,9 +60,9 @@ void phong_function::setParameters(const vec& p)
{
for(int i=0; i<dimY(); ++i)
{
_kd[i] = p[i];
_ks[i] = p[i+ dimY()];
_N[i] = p[i+2*dimY()];
_kd[i] = p[i*3 + 0];
_ks[i] = p[i*3 + 1];
_N[i] = p[i*3 + 2];
}
}
......@@ -71,20 +71,30 @@ void phong_function::setParameters(const vec& p)
vec phong_function::parametersJacobian(const vec& x) const
{
vec jac(dimY()*nbParameters());
for(int i=0; i<dimY(); ++i)
{
// df / dk_d
jac[i*nbParameters() + 0] = 1.0;
for(int i=0; i<dimY(); ++i)
for(int j=0; j<dimY(); ++j)
{
if(i == j)
{
// df / dk_d
jac[i*nbParameters() + j*3+0] = 1.0;
// df / dk_s
jac[i*nbParameters() + 1] = std::pow(x[0], _N[0]);
// df / dk_s
jac[i*nbParameters() + j*3+1] = std::pow(x[0], _N[j]);
// 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]);
}
// df / dN
if(x[0] == 0.0)
jac[i*nbParameters() + j*3+2] = 0.0;
else
jac[i*nbParameters() + j*3+2] = _ks[j] * std::log(x[0]) * std::pow(x[0], _N[j]);
}
else
{
jac[i*nbParameters() + j*3+0] = 0.0;
jac[i*nbParameters() + j*3+1] = 0.0;
jac[i*nbParameters() + j*3+2] = 0.0;
}
}
return jac;
}
......
......@@ -21,7 +21,7 @@ fitter* provide_fitter()
struct EigenFunctor: Eigen::DenseFunctor<double>
{
EigenFunctor(nonlinear_function* f, const data* d) :
DenseFunctor<double>(f->nbParameters(), d->size()), _f(f), _d(d)
DenseFunctor<double>(f->nbParameters(), d->dimY()*d->size()), _f(f), _d(d)
{
#ifndef DEBUG
std::cout << "<<DEBUG>> constructing an EigenFunctor for n=" << inputs() << " parameters and m=" << values() << " points" << std::endl ;
......@@ -30,18 +30,16 @@ struct EigenFunctor: Eigen::DenseFunctor<double>
int operator()(const Eigen::VectorXd& x, Eigen::VectorXd& y) const
{
std::cout << "input vector = " << x << std::endl;
#ifndef DEBUG
std::cout << "parameters:" << std::endl << x << std::endl << std::endl ;
#endif
// Update the parameters vector
vec _p(inputs());
for(int i=0; i<inputs(); ++i)
{
_p[i] = x(i);
}
for(int i=0; i<inputs(); ++i) { _p[i] = x(i); }
_f->setParameters(_p);
for(int s=0; s<values(); ++s)
for(int s=0; s<_d->size(); ++s)
{
vec _x = _d->get(s);
......@@ -51,29 +49,56 @@ struct EigenFunctor: Eigen::DenseFunctor<double>
// Should add the resulting vector completely
vec _y = _di - (*_f)(_x);
y(s) = _y[0];
for(int i=0; i<_f->dimY(); ++i)
y(i*_d->size() + s) = _y[i];
}
#ifndef DEBUG
std::cout << "diff vector:" << std::endl << y << std::endl << std::endl ;
#endif
return 0;
}
int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const
{
// Update the paramters
vec _p(inputs());
for(int i=0; i<inputs(); ++i) _p[i] = x(i);
for(int i=0; i<inputs(); ++i) { _p[i] = x(i); }
_f->setParameters(_p);
for(int i=0; i<values(); ++i)
// For each element to fit, fill the rows of the matrix
for(int s=0; s<_d->size(); ++s)
{
vec xi = _d->get(i);
// Get the position
vec xi = _d->get(s);
// Get the associated jacobian
vec _jac = _f->parametersJacobian(xi);
// Fill the columns of the matrix
#ifdef DEBUG
Eigen::MatrixXd temp (_f->dimY(), _f->nbParameters());
#endif
for(int j=0; j<_f->nbParameters(); ++j)
{
fjac(i,j) = -_jac[j];
// For each output channel, update the subpart of the
// vector row
for(int i=0; i<_f->dimY(); ++i)
{
fjac(i*_d->size() + s, j) = -_jac[i*_f->nbParameters() + j];
#ifdef DEBUG
temp(i, j) = _jac[i*_f->nbParameters() + j];
#endif
}
}
}
#ifdef DEBUG
std::cout << temp << std::endl << std::endl ;
#endif
}
#ifndef DEBUG
std::cout << "jacobian :" << std::endl << fjac << std::endl << std::endl;
#endif
return 0;
}
......@@ -111,7 +136,7 @@ bool nonlinear_fitter_eigen::fit_data(const data* d, function* fit, const argume
/* the following starting values provide a rough fit. */
int info;
Eigen::VectorXd x;
Eigen::VectorXd x(6);
x.setConstant(nf->nbParameters(), 1.);
EigenFunctor functor(nf, d);
......
......@@ -68,14 +68,16 @@ int main(int argc, char** argv)
}
else if(k == 5)
{
f << "#DIM 1 1" << std::endl ;
f << "#DIM 1 3" << 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 = 0.1 + 0.5 * std::pow(x, 1.5) ;
const float z1 = 0.1 + 0.5 * std::pow(x, 1.5) ;
const float z2 = 0.0 - 0.1 * std::pow(x, 4.0) ;
const float z3 = 0.5 + 0.7 * std::pow(x, 1.0) ;
f << x << "\t" << z << std::endl ;
f << x << "\t" << z1 << "\t" << z2 << "\t" << z3 << std::endl ;
}
}
else if(k == 6)
......@@ -85,7 +87,7 @@ int main(int argc, char** argv)
for(int i=0; i<nbx; ++i)
{
const float x = i / (float)nbx ;
const float z = 1.0 - x ;
const float z = 0.1 + 0.5 * std::pow(x, 1.5) ;
f << x << "\t" << z << std::endl ;
}
......
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