35 #ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_LEVMARQFITTER1D_H
36 #define OPENMS_TRANSFORMATIONS_FEATUREFINDER_LEVMARQFITTER1D_H
40 #include <gsl/gsl_rng.h>
41 #include <gsl/gsl_randist.h>
42 #include <gsl/gsl_vector.h>
43 #include <gsl/gsl_multifit_nlin.h>
44 #include <gsl/gsl_blas.h>
64 this->defaults_.setValue(
"max_iteration", 500,
"Maximum number of iterations using by Levenberg-Marquardt algorithm.",
StringList::create(
"advanced"));
65 this->defaults_.setValue(
"deltaAbsError", 0.0001,
"Absolute error used by the Levenberg-Marquardt algorithm.",
StringList::create(
"advanced"));
66 this->defaults_.setValue(
"deltaRelError", 0.0001,
"Relative error used by the Levenberg-Marquardt algorithm.",
StringList::create(
"advanced"));
72 max_iteration_(source.max_iteration_),
73 abs_error_(source.abs_error_),
74 rel_error_(source.rel_error_)
86 if (&source ==
this)
return *
this;
113 virtual void printState_(
Int iter, gsl_multifit_fdfsolver * s) = 0;
118 return gsl_strerror(gsl_status_);
127 Int (*
residual)(
const gsl_vector * x,
void * params, gsl_vector * f),
128 Int (*
jacobian)(
const gsl_vector * x,
void * params, gsl_matrix * J),
129 Int (*
evaluate)(
const gsl_vector * x,
void * params, gsl_vector * f, gsl_matrix * J),
130 void * advanced_params
134 const gsl_multifit_fdfsolver_type * T;
135 gsl_multifit_fdfsolver * s;
146 if (n < p)
throw Exception::UnableToFit(__FILE__, __LINE__, __PRETTY_FUNCTION__,
"UnableToFit-FinalSet",
"Skipping feature, gsl always expects N>=p");
149 gsl_matrix * covar = gsl_matrix_alloc(p, p);
150 gsl_multifit_function_fdf f;
152 gsl_vector_view x = gsl_vector_view_array(x_init, p);
162 f.params = advanced_params;
164 T = gsl_multifit_fdfsolver_lmsder;
165 s = gsl_multifit_fdfsolver_alloc(T, n, p);
166 gsl_multifit_fdfsolver_set(s, &f, &x.vector);
168 #ifdef DEBUG_FEATUREFINDER
169 printState_(iter, s);
178 status = gsl_multifit_fdfsolver_iterate(s);
180 #ifdef DEBUG_FEATUREFINDER
182 printState_(iter, s);
189 status = gsl_multifit_test_delta(s->dx, s->x, abs_error_, rel_error_);
191 while (status == GSL_CONTINUE && iter < max_iteration_);
195 gsl_multifit_covar(s->J, 0.0, covar);
197 #ifdef DEBUG_FEATUREFINDER
198 gsl_matrix_fprintf(stdout, covar,
"covar %g");
201 #define FIT(i) gsl_vector_get(s->x, i)
202 #define ERR(i) sqrt(gsl_matrix_get(covar, i, i))
205 gsl_status_ = status;
207 #ifdef DEBUG_FEATUREFINDER
214 printf(
"chisq/dof = %g\n", pow(chi, 2.0) / dof);
216 for (
Size i = 0; i < p; ++i)
219 printf(
".Parameter = %.5f +/- %.5f\n",
FIT(i), c *
ERR(i));
225 for (
Size i = 0; i < p; ++i)
230 gsl_multifit_fdfsolver_free(s);
231 gsl_matrix_free(covar);
238 max_iteration_ = this->param_.getValue(
"max_iteration");
239 abs_error_ = this->param_.getValue(
"deltaAbsError");
240 rel_error_ = this->param_.getValue(
"deltaRelError");
246 #endif // OPENMS_TRANSFORMATIONS_FEATUREFINDER_LEVMARQFITTER1D_H
Abstract class for 1D-model fitter using Levenberg-Marquardt algorithm for parameter optimization...
Definition: LevMarqFitter1D.h:52
A more convenient string class.
Definition: String.h:56
virtual ~LevMarqFitter1D()
destructor
Definition: LevMarqFitter1D.h:79
Int max_iteration_
Maximum number of iterations.
Definition: LevMarqFitter1D.h:103
virtual LevMarqFitter1D & operator=(const LevMarqFitter1D &source)
assignment operator
Definition: LevMarqFitter1D.h:84
LevMarqFitter1D(const LevMarqFitter1D &source)
copy constructor
Definition: LevMarqFitter1D.h:70
unsigned int UInt
Unsigned integer type.
Definition: Types.h:92
virtual Fitter1D & operator=(const Fitter1D &source)
assignment operator
CoordinateType rel_error_
Relative error.
Definition: LevMarqFitter1D.h:108
int evaluate(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J)
Driver function for the evaluation of function and jacobian.
void optimize_(const RawDataArrayType &set, Int num_params, CoordinateType x_init[], Int(*residual)(const gsl_vector *x, void *params, gsl_vector *f), Int(*jacobian)(const gsl_vector *x, void *params, gsl_matrix *J), Int(*evaluate)(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J), void *advanced_params)
Optimize start parameter.
Definition: LevMarqFitter1D.h:126
LevMarqFitter1D()
Default constructor.
Definition: LevMarqFitter1D.h:61
bool symmetric_
Parameter indicates symmetric peaks.
Definition: LevMarqFitter1D.h:101
static StringList create(const String &list, const char splitter= ',')
Returns a list that is created by splitting the given (comma-separated) string (String are not trimme...
const String getGslStatus_()
Return GSL status as string.
Definition: LevMarqFitter1D.h:116
int residual(const gsl_vector *x, void *params, gsl_vector *f)
Evaluation of the target function for nonlinear optimization.
Exception used if an error occurred while fitting a model to a given dataset.
Definition: Exception.h:662
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
Abstract base class for all 1D-dimensional model fitter.
Definition: Fitter1D.h:59
int jacobian(const gsl_vector *x, void *params, gsl_matrix *J)
Compute the Jacobian of the residual, where each row of the matrix corresponds to a point in the data...
std::vector< double > ContainerType
Definition: LevMarqFitter1D.h:58
virtual void updateMembers_()
This method is used to update extra member variables at the end of the setParameters() method...
int Int
Signed integer type.
Definition: Types.h:100
CoordinateType abs_error_
Absolute error.
Definition: LevMarqFitter1D.h:106
Int gsl_status_
GSL status.
Definition: LevMarqFitter1D.h:99
std::vector< PeakType > RawDataArrayType
Raw data container type using for the temporary storage of the input data.
Definition: Fitter1D.h:76
void updateMembers_()
This method is used to update extra member variables at the end of the setParameters() method...
Definition: LevMarqFitter1D.h:235