23 #include <gsl/gsl_multifit_nlinear.h>
24 #include <gsl/gsl_blas.h>
35 return std::make_shared<GSLEngine>(max_iterations);
41 m_itmax{itmax}, m_xtol{xtol}, m_gtol{gtol}, m_ftol{ftol}, m_delta{delta} {
43 gsl_set_error_handler_off();
70 return gsl_vector_get(
m_v,
m_i);
74 return *gsl_vector_ptr(
m_v,
m_i);
81 const gsl_vector *
m_v;
102 return gsl_vector_get(
m_v,
m_i);
111 auto adata =
std::tie(parameter_manager, residual_estimator);
114 const gsl_multifit_nlinear_type *type = gsl_multifit_nlinear_trust;
119 gsl_multifit_nlinear_parameters params = gsl_multifit_nlinear_default_parameters();
123 params.trs = gsl_multifit_nlinear_trs_lm;
126 params.scale = gsl_multifit_nlinear_scale_levenberg;
128 params.solver = gsl_multifit_nlinear_solver_cholesky;
130 params.fdtype = GSL_MULTIFIT_NLINEAR_FWDIFF;
135 gsl_multifit_nlinear_workspace *workspace = gsl_multifit_nlinear_alloc(
139 if (workspace ==
nullptr) {
146 gsl_vector_view gsl_param_view = gsl_vector_view_array(param_values.
data(), param_values.
size());
149 auto function = [](
const gsl_vector*
x,
void *extra, gsl_vector *f) ->
int {
150 auto *extra_ptr = (decltype(adata) *) extra;
157 gsl_multifit_nlinear_fdf fdf;
166 gsl_multifit_nlinear_init(&gsl_param_view.vector, &fdf, workspace);
170 gsl_vector *residual = gsl_multifit_nlinear_residual(workspace);
171 gsl_blas_ddot(residual, residual, &chisq0);
175 int ret = gsl_multifit_nlinear_driver(
188 gsl_blas_ddot(residual, residual, &chisq);
195 summary.
iteration_no = gsl_multifit_nlinear_niter(workspace);
199 gsl_matrix *J = gsl_multifit_nlinear_jac(workspace);
200 gsl_matrix_view covar = gsl_matrix_view_array(covariance_matrix.
data(), parameter_manager.
numberOfParameters(),
202 gsl_multifit_nlinear_covar(J, 0.0, &covar.matrix);
210 int levmar_reason = 0;
211 if (ret == GSL_SUCCESS) {
212 levmar_reason = (info == 1) ? 2 : 1;
214 else if (ret == GSL_EMAXITER) {
221 gsl_blas_dnrm2(workspace->g),
222 gsl_blas_dnrm2(workspace->dx),
225 static_cast<double>(levmar_reason),
226 static_cast<double>(fdf.nevalf),
227 static_cast<double>(fdf.nevaldf),
232 gsl_multifit_nlinear_free(workspace);