7 #include <gsl/gsl_multifit_nlinear.h>
8 #include <gsl/gsl_blas.h>
15 namespace ModelFitting {
19 return std::make_shared<GSLEngine>(max_iterations);
25 m_itmax{itmax}, m_xtol{xtol}, m_gtol{gtol}, m_ftol{ftol}, m_delta{delta} {
27 gsl_set_error_handler_off();
54 return gsl_vector_get(
m_v,
m_i);
58 return *gsl_vector_ptr(
m_v,
m_i);
65 const gsl_vector *
m_v;
86 return gsl_vector_get(
m_v,
m_i);
95 auto adata =
std::tie(parameter_manager, residual_estimator);
98 const gsl_multifit_nlinear_type *type = gsl_multifit_nlinear_trust;
103 gsl_multifit_nlinear_parameters params = gsl_multifit_nlinear_default_parameters();
107 params.trs = gsl_multifit_nlinear_trs_lm;
110 params.scale = gsl_multifit_nlinear_scale_levenberg;
112 params.solver = gsl_multifit_nlinear_solver_cholesky;
114 params.fdtype = GSL_MULTIFIT_NLINEAR_FWDIFF;
119 gsl_multifit_nlinear_workspace *workspace = gsl_multifit_nlinear_alloc(
123 if (workspace ==
nullptr) {
130 gsl_vector_view gsl_param_view = gsl_vector_view_array(param_values.data(), param_values.size());
133 auto function = [](
const gsl_vector*
x,
void *extra, gsl_vector *f) ->
int {
134 auto *extra_ptr = (decltype(adata) *) extra;
141 gsl_multifit_nlinear_fdf fdf;
150 gsl_multifit_nlinear_init(&gsl_param_view.vector, &fdf, workspace);
154 gsl_vector *residual = gsl_multifit_nlinear_residual(workspace);
155 gsl_blas_ddot(residual, residual, &chisq0);
159 int ret = gsl_multifit_nlinear_driver(
172 gsl_blas_ddot(residual, residual, &chisq);
178 summary.success_flag = (ret == GSL_SUCCESS);
179 summary.iteration_no = gsl_multifit_nlinear_niter(workspace);
180 summary.parameter_sigmas = {};
183 gsl_matrix *J = gsl_multifit_nlinear_jac(workspace);
184 gsl_matrix_view covar = gsl_matrix_view_array(covariance_matrix.data(), parameter_manager.
numberOfParameters(),
186 gsl_multifit_nlinear_covar(J, 0.0, &covar.matrix);
194 int levmar_reason = 0;
195 if (ret == GSL_SUCCESS) {
196 levmar_reason = (info == 1) ? 2 : 1;
198 else if (ret == GSL_EMAXITER) {
205 gsl_blas_dnrm2(workspace->g),
206 gsl_blas_dnrm2(workspace->dx),
208 static_cast<double>(summary.iteration_no),
209 static_cast<double>(levmar_reason),
210 static_cast<double>(fdf.nevalf),
211 static_cast<double>(fdf.nevaldf),
216 gsl_multifit_nlinear_free(workspace);
GslVectorConstIterator(const gsl_vector *v)
GslVectorIterator & operator++()
Class containing the summary information of solving a least square minimization problem.
void populateResiduals(DoubleIter output_iter) const
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
void getEngineValues(DoubleIter output_iter) const
Returns the engine values of the managed parameters.
GSLEngine(int itmax=1000, double xtol=1e-8, double gtol=1e-8, double ftol=1e-8, double delta=1e-4)
Constructs a new instance of the engine.
LeastSquareSummary solveProblem(EngineParameterManager ¶meter_manager, ResidualEstimator &residual_estimator) override
GslVectorConstIterator & operator++()
void updateEngineValues(DoubleIter new_values_iter)
Updates the managed parameters with the given engine values.
static std::shared_ptr< LeastSquareEngine > createLevmarEngine(unsigned max_iterations)
std::size_t numberOfResiduals() const
GslVectorIterator(gsl_vector *v)
static LeastSquareEngineManager::StaticEngine levmar_engine
GslVectorIterator operator++(int)
Class responsible for managing the parameters the least square engine minimizes.
std::vector< double > convertCovarianceMatrixToWorldSpace(std::vector< double > covariance_matrix) const
GslVectorConstIterator operator++(int)
Provides to the LeastSquareEngine the residual values.
std::size_t numberOfParameters()
Returns the number of parameters managed by the manager.