Put the Status outside of the class, it really does not depend on the

FunctorType or Scalar template parameters.
This commit is contained in:
Thomas Capricelli 2010-01-28 10:29:26 +01:00
parent 27cf1b3a97
commit 98a584ceb8
2 changed files with 134 additions and 127 deletions

View File

@ -28,6 +28,19 @@
#ifndef EIGEN_HYBRIDNONLINEARSOLVER_H #ifndef EIGEN_HYBRIDNONLINEARSOLVER_H
#define EIGEN_HYBRIDNONLINEARSOLVER_H #define EIGEN_HYBRIDNONLINEARSOLVER_H
namespace HybridNonLinearSolverSpace {
enum Status {
Running = -1,
ImproperInputParameters = 0,
RelativeErrorTooSmall = 1,
TooManyFunctionEvaluation = 2,
TolTooSmall = 3,
NotMakingProgressJacobian = 4,
NotMakingProgressIterations = 5,
UserAksed = 6
};
}
/** /**
* \ingroup NonLinearOptimization_Module * \ingroup NonLinearOptimization_Module
* \brief Finds a zero of a system of n * \brief Finds a zero of a system of n
@ -46,17 +59,6 @@ public:
HybridNonLinearSolver(FunctorType &_functor) HybridNonLinearSolver(FunctorType &_functor)
: functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; } : functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; }
enum Status {
Running = -1,
ImproperInputParameters = 0,
RelativeErrorTooSmall = 1,
TooManyFunctionEvaluation = 2,
TolTooSmall = 3,
NotMakingProgressJacobian = 4,
NotMakingProgressIterations = 5,
UserAksed = 6
};
struct Parameters { struct Parameters {
Parameters() Parameters()
: factor(Scalar(100.)) : factor(Scalar(100.))
@ -77,38 +79,38 @@ public:
/* TODO: if eigen provides a triangular storage, use it here */ /* TODO: if eigen provides a triangular storage, use it here */
typedef Matrix< Scalar, Dynamic, Dynamic > UpperTriangularType; typedef Matrix< Scalar, Dynamic, Dynamic > UpperTriangularType;
Status hybrj1( HybridNonLinearSolverSpace::Status hybrj1(
FVectorType &x, FVectorType &x,
const Scalar tol = ei_sqrt(epsilon<Scalar>()) const Scalar tol = ei_sqrt(epsilon<Scalar>())
); );
Status solveInit( HybridNonLinearSolverSpace::Status solveInit(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
Status solveOneStep( HybridNonLinearSolverSpace::Status solveOneStep(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
Status solve( HybridNonLinearSolverSpace::Status solve(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
Status hybrd1( HybridNonLinearSolverSpace::Status hybrd1(
FVectorType &x, FVectorType &x,
const Scalar tol = ei_sqrt(epsilon<Scalar>()) const Scalar tol = ei_sqrt(epsilon<Scalar>())
); );
Status solveNumericalDiffInit( HybridNonLinearSolverSpace::Status solveNumericalDiffInit(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
Status solveNumericalDiffOneStep( HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
Status solveNumericalDiff( HybridNonLinearSolverSpace::Status solveNumericalDiff(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
@ -142,7 +144,7 @@ private:
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename HybridNonLinearSolver<FunctorType,Scalar>::Status HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::hybrj1( HybridNonLinearSolver<FunctorType,Scalar>::hybrj1(
FVectorType &x, FVectorType &x,
const Scalar tol const Scalar tol
@ -152,7 +154,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrj1(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || tol < 0.) if (n <= 0 || tol < 0.)
return ImproperInputParameters; return HybridNonLinearSolverSpace::ImproperInputParameters;
resetParameters(); resetParameters();
parameters.maxfev = 100*(n+1); parameters.maxfev = 100*(n+1);
@ -165,7 +167,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrj1(
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename HybridNonLinearSolver<FunctorType,Scalar>::Status HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveInit( HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
FVectorType &x, FVectorType &x,
const int mode const int mode
@ -187,17 +189,17 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0. ) if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0. )
return ImproperInputParameters; return HybridNonLinearSolverSpace::ImproperInputParameters;
if (mode == 2) if (mode == 2)
for (int j = 0; j < n; ++j) for (int j = 0; j < n; ++j)
if (diag[j] <= 0.) if (diag[j] <= 0.)
return ImproperInputParameters; return HybridNonLinearSolverSpace::ImproperInputParameters;
/* evaluate the function at the starting point */ /* evaluate the function at the starting point */
/* and calculate its norm. */ /* and calculate its norm. */
nfev = 1; nfev = 1;
if ( functor(x, fvec) < 0) if ( functor(x, fvec) < 0)
return UserAksed; return HybridNonLinearSolverSpace::UserAksed;
fnorm = fvec.stableNorm(); fnorm = fvec.stableNorm();
/* initialize iteration counter and monitors. */ /* initialize iteration counter and monitors. */
@ -207,11 +209,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
nslow1 = 0; nslow1 = 0;
nslow2 = 0; nslow2 = 0;
return Running; return HybridNonLinearSolverSpace::Running;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename HybridNonLinearSolver<FunctorType,Scalar>::Status HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep( HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
FVectorType &x, FVectorType &x,
const int mode const int mode
@ -224,7 +226,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
/* calculate the jacobian matrix. */ /* calculate the jacobian matrix. */
if ( functor.df(x, fjac) < 0) if ( functor.df(x, fjac) < 0)
return UserAksed; return HybridNonLinearSolverSpace::UserAksed;
++njev; ++njev;
wa2 = fjac.colwise().blueNorm(); wa2 = fjac.colwise().blueNorm();
@ -276,7 +278,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
if ( functor(wa2, wa4) < 0) if ( functor(wa2, wa4) < 0)
return UserAksed; return HybridNonLinearSolverSpace::UserAksed;
++nfev; ++nfev;
fnorm1 = wa4.stableNorm(); fnorm1 = wa4.stableNorm();
@ -334,17 +336,17 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
/* test for convergence. */ /* test for convergence. */
if (delta <= parameters.xtol * xnorm || fnorm == 0.) if (delta <= parameters.xtol * xnorm || fnorm == 0.)
return RelativeErrorTooSmall; return HybridNonLinearSolverSpace::RelativeErrorTooSmall;
/* tests for termination and stringent tolerances. */ /* tests for termination and stringent tolerances. */
if (nfev >= parameters.maxfev) if (nfev >= parameters.maxfev)
return TooManyFunctionEvaluation; return HybridNonLinearSolverSpace::TooManyFunctionEvaluation;
if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm) if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
return TolTooSmall; return HybridNonLinearSolverSpace::TolTooSmall;
if (nslow2 == 5) if (nslow2 == 5)
return NotMakingProgressJacobian; return HybridNonLinearSolverSpace::NotMakingProgressJacobian;
if (nslow1 == 10) if (nslow1 == 10)
return NotMakingProgressIterations; return HybridNonLinearSolverSpace::NotMakingProgressIterations;
/* criterion for recalculating jacobian. */ /* criterion for recalculating jacobian. */
if (ncfail == 2) if (ncfail == 2)
@ -365,18 +367,18 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
jeval = false; jeval = false;
} }
return Running; return HybridNonLinearSolverSpace::Running;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename HybridNonLinearSolver<FunctorType,Scalar>::Status HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solve( HybridNonLinearSolver<FunctorType,Scalar>::solve(
FVectorType &x, FVectorType &x,
const int mode const int mode
) )
{ {
Status status = solveInit(x, mode); HybridNonLinearSolverSpace::Status status = solveInit(x, mode);
while (status==Running) while (status==HybridNonLinearSolverSpace::Running)
status = solveOneStep(x, mode); status = solveOneStep(x, mode);
return status; return status;
} }
@ -384,7 +386,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solve(
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename HybridNonLinearSolver<FunctorType,Scalar>::Status HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::hybrd1( HybridNonLinearSolver<FunctorType,Scalar>::hybrd1(
FVectorType &x, FVectorType &x,
const Scalar tol const Scalar tol
@ -394,7 +396,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrd1(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || tol < 0.) if (n <= 0 || tol < 0.)
return ImproperInputParameters; return HybridNonLinearSolverSpace::ImproperInputParameters;
resetParameters(); resetParameters();
parameters.maxfev = 200*(n+1); parameters.maxfev = 200*(n+1);
@ -408,7 +410,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrd1(
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename HybridNonLinearSolver<FunctorType,Scalar>::Status HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit( HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
FVectorType &x, FVectorType &x,
const int mode const int mode
@ -434,17 +436,17 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.nb_of_subdiagonals< 0 || parameters.nb_of_superdiagonals< 0 || parameters.factor <= 0. ) if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.nb_of_subdiagonals< 0 || parameters.nb_of_superdiagonals< 0 || parameters.factor <= 0. )
return ImproperInputParameters; return HybridNonLinearSolverSpace::ImproperInputParameters;
if (mode == 2) if (mode == 2)
for (int j = 0; j < n; ++j) for (int j = 0; j < n; ++j)
if (diag[j] <= 0.) if (diag[j] <= 0.)
return ImproperInputParameters; return HybridNonLinearSolverSpace::ImproperInputParameters;
/* evaluate the function at the starting point */ /* evaluate the function at the starting point */
/* and calculate its norm. */ /* and calculate its norm. */
nfev = 1; nfev = 1;
if ( functor(x, fvec) < 0) if ( functor(x, fvec) < 0)
return UserAksed; return HybridNonLinearSolverSpace::UserAksed;
fnorm = fvec.stableNorm(); fnorm = fvec.stableNorm();
/* initialize iteration counter and monitors. */ /* initialize iteration counter and monitors. */
@ -454,11 +456,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
nslow1 = 0; nslow1 = 0;
nslow2 = 0; nslow2 = 0;
return Running; return HybridNonLinearSolverSpace::Running;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename HybridNonLinearSolver<FunctorType,Scalar>::Status HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep( HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
FVectorType &x, FVectorType &x,
const int mode const int mode
@ -473,7 +475,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
/* calculate the jacobian matrix. */ /* calculate the jacobian matrix. */
if (ei_fdjac1(functor, x, fvec, fjac, parameters.nb_of_subdiagonals, parameters.nb_of_superdiagonals, parameters.epsfcn) <0) if (ei_fdjac1(functor, x, fvec, fjac, parameters.nb_of_subdiagonals, parameters.nb_of_superdiagonals, parameters.epsfcn) <0)
return UserAksed; return HybridNonLinearSolverSpace::UserAksed;
nfev += std::min(parameters.nb_of_subdiagonals+parameters.nb_of_superdiagonals+ 1, n); nfev += std::min(parameters.nb_of_subdiagonals+parameters.nb_of_superdiagonals+ 1, n);
wa2 = fjac.colwise().blueNorm(); wa2 = fjac.colwise().blueNorm();
@ -525,7 +527,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
if ( functor(wa2, wa4) < 0) if ( functor(wa2, wa4) < 0)
return UserAksed; return HybridNonLinearSolverSpace::UserAksed;
++nfev; ++nfev;
fnorm1 = wa4.stableNorm(); fnorm1 = wa4.stableNorm();
@ -583,17 +585,17 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
/* test for convergence. */ /* test for convergence. */
if (delta <= parameters.xtol * xnorm || fnorm == 0.) if (delta <= parameters.xtol * xnorm || fnorm == 0.)
return RelativeErrorTooSmall; return HybridNonLinearSolverSpace::RelativeErrorTooSmall;
/* tests for termination and stringent tolerances. */ /* tests for termination and stringent tolerances. */
if (nfev >= parameters.maxfev) if (nfev >= parameters.maxfev)
return TooManyFunctionEvaluation; return HybridNonLinearSolverSpace::TooManyFunctionEvaluation;
if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm) if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
return TolTooSmall; return HybridNonLinearSolverSpace::TolTooSmall;
if (nslow2 == 5) if (nslow2 == 5)
return NotMakingProgressJacobian; return HybridNonLinearSolverSpace::NotMakingProgressJacobian;
if (nslow1 == 10) if (nslow1 == 10)
return NotMakingProgressIterations; return HybridNonLinearSolverSpace::NotMakingProgressIterations;
/* criterion for recalculating jacobian. */ /* criterion for recalculating jacobian. */
if (ncfail == 2) if (ncfail == 2)
@ -614,18 +616,18 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
jeval = false; jeval = false;
} }
return Running; return HybridNonLinearSolverSpace::Running;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename HybridNonLinearSolver<FunctorType,Scalar>::Status HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff( HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(
FVectorType &x, FVectorType &x,
const int mode const int mode
) )
{ {
Status status = solveNumericalDiffInit(x, mode); HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x, mode);
while (status==Running) while (status==HybridNonLinearSolverSpace::Running)
status = solveNumericalDiffOneStep(x, mode); status = solveNumericalDiffOneStep(x, mode);
return status; return status;
} }

View File

@ -28,21 +28,8 @@
#ifndef EIGEN_LEVENBERGMARQUARDT__H #ifndef EIGEN_LEVENBERGMARQUARDT__H
#define EIGEN_LEVENBERGMARQUARDT__H #define EIGEN_LEVENBERGMARQUARDT__H
/**
* \ingroup NonLinearOptimization_Module
* \brief Performs non linear optimization over a non-linear function,
* using a variant of the Levenberg Marquardt algorithm.
*
* Check wikipedia for more information.
* http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
*/
template<typename FunctorType, typename Scalar=double>
class LevenbergMarquardt
{
public:
LevenbergMarquardt(FunctorType &_functor)
: functor(_functor) { nfev = njev = iter = 0; fnorm=gnorm = 0.; }
namespace LevenbergMarquardtSpace {
enum Status { enum Status {
NotStarted = -2, NotStarted = -2,
Running = -1, Running = -1,
@ -57,6 +44,24 @@ public:
GtolTooSmall = 8, GtolTooSmall = 8,
UserAsked = 9 UserAsked = 9
}; };
}
/**
* \ingroup NonLinearOptimization_Module
* \brief Performs non linear optimization over a non-linear function,
* using a variant of the Levenberg Marquardt algorithm.
*
* Check wikipedia for more information.
* http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
*/
template<typename FunctorType, typename Scalar=double>
class LevenbergMarquardt
{
public:
LevenbergMarquardt(FunctorType &_functor)
: functor(_functor) { nfev = njev = iter = 0; fnorm=gnorm = 0.; }
struct Parameters { struct Parameters {
Parameters() Parameters()
@ -77,45 +82,45 @@ public:
typedef Matrix< Scalar, Dynamic, 1 > FVectorType; typedef Matrix< Scalar, Dynamic, 1 > FVectorType;
typedef Matrix< Scalar, Dynamic, Dynamic > JacobianType; typedef Matrix< Scalar, Dynamic, Dynamic > JacobianType;
Status lmder1( LevenbergMarquardtSpace::Status lmder1(
FVectorType &x, FVectorType &x,
const Scalar tol = ei_sqrt(epsilon<Scalar>()) const Scalar tol = ei_sqrt(epsilon<Scalar>())
); );
Status minimize( LevenbergMarquardtSpace::Status minimize(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
Status minimizeInit( LevenbergMarquardtSpace::Status minimizeInit(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
Status minimizeOneStep( LevenbergMarquardtSpace::Status minimizeOneStep(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
static Status lmdif1( static LevenbergMarquardtSpace::Status lmdif1(
FunctorType &functor, FunctorType &functor,
FVectorType &x, FVectorType &x,
int *nfev, int *nfev,
const Scalar tol = ei_sqrt(epsilon<Scalar>()) const Scalar tol = ei_sqrt(epsilon<Scalar>())
); );
Status lmstr1( LevenbergMarquardtSpace::Status lmstr1(
FVectorType &x, FVectorType &x,
const Scalar tol = ei_sqrt(epsilon<Scalar>()) const Scalar tol = ei_sqrt(epsilon<Scalar>())
); );
Status minimizeOptimumStorage( LevenbergMarquardtSpace::Status minimizeOptimumStorage(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
Status minimizeOptimumStorageInit( LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
Status minimizeOptimumStorageOneStep( LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(
FVectorType &x, FVectorType &x,
const int mode=1 const int mode=1
); );
@ -146,7 +151,7 @@ private:
}; };
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename LevenbergMarquardt<FunctorType,Scalar>::Status LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::lmder1( LevenbergMarquardt<FunctorType,Scalar>::lmder1(
FVectorType &x, FVectorType &x,
const Scalar tol const Scalar tol
@ -157,7 +162,7 @@ LevenbergMarquardt<FunctorType,Scalar>::lmder1(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || m < n || tol < 0.) if (n <= 0 || m < n || tol < 0.)
return ImproperInputParameters; return LevenbergMarquardtSpace::ImproperInputParameters;
resetParameters(); resetParameters();
parameters.ftol = tol; parameters.ftol = tol;
@ -169,21 +174,21 @@ LevenbergMarquardt<FunctorType,Scalar>::lmder1(
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename LevenbergMarquardt<FunctorType,Scalar>::Status LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimize( LevenbergMarquardt<FunctorType,Scalar>::minimize(
FVectorType &x, FVectorType &x,
const int mode const int mode
) )
{ {
Status status = minimizeInit(x, mode); LevenbergMarquardtSpace::Status status = minimizeInit(x, mode);
do { do {
status = minimizeOneStep(x, mode); status = minimizeOneStep(x, mode);
} while (status==Running); } while (status==LevenbergMarquardtSpace::Running);
return status; return status;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename LevenbergMarquardt<FunctorType,Scalar>::Status LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeInit( LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
FVectorType &x, FVectorType &x,
const int mode const int mode
@ -207,29 +212,29 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.) if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
return ImproperInputParameters; return LevenbergMarquardtSpace::ImproperInputParameters;
if (mode == 2) if (mode == 2)
for (int j = 0; j < n; ++j) for (int j = 0; j < n; ++j)
if (diag[j] <= 0.) if (diag[j] <= 0.)
return ImproperInputParameters; return LevenbergMarquardtSpace::ImproperInputParameters;
/* evaluate the function at the starting point */ /* evaluate the function at the starting point */
/* and calculate its norm. */ /* and calculate its norm. */
nfev = 1; nfev = 1;
if ( functor(x, fvec) < 0) if ( functor(x, fvec) < 0)
return UserAsked; return LevenbergMarquardtSpace::UserAsked;
fnorm = fvec.stableNorm(); fnorm = fvec.stableNorm();
/* initialize levenberg-marquardt parameter and iteration counter. */ /* initialize levenberg-marquardt parameter and iteration counter. */
par = 0.; par = 0.;
iter = 1; iter = 1;
return NotStarted; return LevenbergMarquardtSpace::NotStarted;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename LevenbergMarquardt<FunctorType,Scalar>::Status LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep( LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
FVectorType &x, FVectorType &x,
const int mode const int mode
@ -240,7 +245,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
/* calculate the jacobian matrix. */ /* calculate the jacobian matrix. */
int df_ret = functor.df(x, fjac); int df_ret = functor.df(x, fjac);
if (df_ret<0) if (df_ret<0)
return UserAsked; return LevenbergMarquardtSpace::UserAsked;
if (df_ret>0) if (df_ret>0)
// numerical diff, we evaluated the function df_ret times // numerical diff, we evaluated the function df_ret times
nfev += df_ret; nfev += df_ret;
@ -282,7 +287,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
/* test for convergence of the gradient norm. */ /* test for convergence of the gradient norm. */
if (gnorm <= parameters.gtol) if (gnorm <= parameters.gtol)
return CosinusTooSmall; return LevenbergMarquardtSpace::CosinusTooSmall;
/* rescale if necessary. */ /* rescale if necessary. */
if (mode != 2) if (mode != 2)
@ -304,7 +309,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
if ( functor(wa2, wa4) < 0) if ( functor(wa2, wa4) < 0)
return UserAsked; return LevenbergMarquardtSpace::UserAsked;
++nfev; ++nfev;
fnorm1 = wa4.stableNorm(); fnorm1 = wa4.stableNorm();
@ -356,29 +361,29 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
/* tests for convergence. */ /* tests for convergence. */
if (ei_abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm) if (ei_abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm)
return RelativeErrorAndReductionTooSmall; return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall;
if (ei_abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.) if (ei_abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.)
return RelativeReductionTooSmall; return LevenbergMarquardtSpace::RelativeReductionTooSmall;
if (delta <= parameters.xtol * xnorm) if (delta <= parameters.xtol * xnorm)
return RelativeErrorTooSmall; return LevenbergMarquardtSpace::RelativeErrorTooSmall;
/* tests for termination and stringent tolerances. */ /* tests for termination and stringent tolerances. */
if (nfev >= parameters.maxfev) if (nfev >= parameters.maxfev)
return TooManyFunctionEvaluation; return LevenbergMarquardtSpace::TooManyFunctionEvaluation;
if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.) if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.)
return FtolTooSmall; return LevenbergMarquardtSpace::FtolTooSmall;
if (delta <= epsilon<Scalar>() * xnorm) if (delta <= epsilon<Scalar>() * xnorm)
return XtolTooSmall; return LevenbergMarquardtSpace::XtolTooSmall;
if (gnorm <= epsilon<Scalar>()) if (gnorm <= epsilon<Scalar>())
return GtolTooSmall; return LevenbergMarquardtSpace::GtolTooSmall;
} while (ratio < Scalar(1e-4)); } while (ratio < Scalar(1e-4));
return Running; return LevenbergMarquardtSpace::Running;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename LevenbergMarquardt<FunctorType,Scalar>::Status LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::lmstr1( LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
FVectorType &x, FVectorType &x,
const Scalar tol const Scalar tol
@ -389,7 +394,7 @@ LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || m < n || tol < 0.) if (n <= 0 || m < n || tol < 0.)
return ImproperInputParameters; return LevenbergMarquardtSpace::ImproperInputParameters;
resetParameters(); resetParameters();
parameters.ftol = tol; parameters.ftol = tol;
@ -400,7 +405,7 @@ LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename LevenbergMarquardt<FunctorType,Scalar>::Status LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit( LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
FVectorType &x, FVectorType &x,
const int mode const int mode
@ -424,30 +429,30 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.) if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
return ImproperInputParameters; return LevenbergMarquardtSpace::ImproperInputParameters;
if (mode == 2) if (mode == 2)
for (int j = 0; j < n; ++j) for (int j = 0; j < n; ++j)
if (diag[j] <= 0.) if (diag[j] <= 0.)
return ImproperInputParameters; return LevenbergMarquardtSpace::ImproperInputParameters;
/* evaluate the function at the starting point */ /* evaluate the function at the starting point */
/* and calculate its norm. */ /* and calculate its norm. */
nfev = 1; nfev = 1;
if ( functor(x, fvec) < 0) if ( functor(x, fvec) < 0)
return UserAsked; return LevenbergMarquardtSpace::UserAsked;
fnorm = fvec.stableNorm(); fnorm = fvec.stableNorm();
/* initialize levenberg-marquardt parameter and iteration counter. */ /* initialize levenberg-marquardt parameter and iteration counter. */
par = 0.; par = 0.;
iter = 1; iter = 1;
return NotStarted; return LevenbergMarquardtSpace::NotStarted;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename LevenbergMarquardt<FunctorType,Scalar>::Status LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep( LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
FVectorType &x, FVectorType &x,
const int mode const int mode
@ -464,7 +469,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
fjac.fill(0.); fjac.fill(0.);
int rownb = 2; int rownb = 2;
for (i = 0; i < m; ++i) { for (i = 0; i < m; ++i) {
if (functor.df(x, wa3, rownb) < 0) return UserAsked; if (functor.df(x, wa3, rownb) < 0) return LevenbergMarquardtSpace::UserAsked;
ei_rwupdt<Scalar>(fjac, wa3, qtf, fvec[i]); ei_rwupdt<Scalar>(fjac, wa3, qtf, fvec[i]);
++rownb; ++rownb;
} }
@ -528,7 +533,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
/* test for convergence of the gradient norm. */ /* test for convergence of the gradient norm. */
if (gnorm <= parameters.gtol) if (gnorm <= parameters.gtol)
return CosinusTooSmall; return LevenbergMarquardtSpace::CosinusTooSmall;
/* rescale if necessary. */ /* rescale if necessary. */
if (mode != 2) if (mode != 2)
@ -550,7 +555,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
if ( functor(wa2, wa4) < 0) if ( functor(wa2, wa4) < 0)
return UserAsked; return LevenbergMarquardtSpace::UserAsked;
++nfev; ++nfev;
fnorm1 = wa4.stableNorm(); fnorm1 = wa4.stableNorm();
@ -602,43 +607,43 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
/* tests for convergence. */ /* tests for convergence. */
if (ei_abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm) if (ei_abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1. && delta <= parameters.xtol * xnorm)
return RelativeErrorAndReductionTooSmall; return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall;
if (ei_abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.) if (ei_abs(actred) <= parameters.ftol && prered <= parameters.ftol && Scalar(.5) * ratio <= 1.)
return RelativeReductionTooSmall; return LevenbergMarquardtSpace::RelativeReductionTooSmall;
if (delta <= parameters.xtol * xnorm) if (delta <= parameters.xtol * xnorm)
return RelativeErrorTooSmall; return LevenbergMarquardtSpace::RelativeErrorTooSmall;
/* tests for termination and stringent tolerances. */ /* tests for termination and stringent tolerances. */
if (nfev >= parameters.maxfev) if (nfev >= parameters.maxfev)
return TooManyFunctionEvaluation; return LevenbergMarquardtSpace::TooManyFunctionEvaluation;
if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.) if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.)
return FtolTooSmall; return LevenbergMarquardtSpace::FtolTooSmall;
if (delta <= epsilon<Scalar>() * xnorm) if (delta <= epsilon<Scalar>() * xnorm)
return XtolTooSmall; return LevenbergMarquardtSpace::XtolTooSmall;
if (gnorm <= epsilon<Scalar>()) if (gnorm <= epsilon<Scalar>())
return GtolTooSmall; return LevenbergMarquardtSpace::GtolTooSmall;
} while (ratio < Scalar(1e-4)); } while (ratio < Scalar(1e-4));
return Running; return LevenbergMarquardtSpace::Running;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename LevenbergMarquardt<FunctorType,Scalar>::Status LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage( LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(
FVectorType &x, FVectorType &x,
const int mode const int mode
) )
{ {
Status status = minimizeOptimumStorageInit(x, mode); LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x, mode);
do { do {
status = minimizeOptimumStorageOneStep(x, mode); status = minimizeOptimumStorageOneStep(x, mode);
} while (status==Running); } while (status==LevenbergMarquardtSpace::Running);
return status; return status;
} }
template<typename FunctorType, typename Scalar> template<typename FunctorType, typename Scalar>
typename LevenbergMarquardt<FunctorType,Scalar>::Status LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::lmdif1( LevenbergMarquardt<FunctorType,Scalar>::lmdif1(
FunctorType &functor, FunctorType &functor,
FVectorType &x, FVectorType &x,
@ -651,7 +656,7 @@ LevenbergMarquardt<FunctorType,Scalar>::lmdif1(
/* check the input parameters for errors. */ /* check the input parameters for errors. */
if (n <= 0 || m < n || tol < 0.) if (n <= 0 || m < n || tol < 0.)
return ImproperInputParameters; return LevenbergMarquardtSpace::ImproperInputParameters;
NumericalDiff<FunctorType> numDiff(functor); NumericalDiff<FunctorType> numDiff(functor);
// embedded LevenbergMarquardt // embedded LevenbergMarquardt
@ -660,7 +665,7 @@ LevenbergMarquardt<FunctorType,Scalar>::lmdif1(
lm.parameters.xtol = tol; lm.parameters.xtol = tol;
lm.parameters.maxfev = 200*(n+1); lm.parameters.maxfev = 200*(n+1);
Status info = Status(lm.minimize(x)); LevenbergMarquardtSpace::Status info = LevenbergMarquardtSpace::Status(lm.minimize(x));
if (nfev) if (nfev)
* nfev = lm.nfev; * nfev = lm.nfev;
return info; return info;