#include "clipm_priv_optimize.h"
#include "clipm_math.h"
#include "clipm_compatibility_replacements.h"
#include "clipm_priv_error.h"
#include "clipm_priv_matrix.h"
#include "clipm_priv_vector.h"

Functions

cpl_error_code	clipm_priv_optimize_gaussian (const cpl_vector xvalues, const cpl_vector xsigmas, const cpl_vector yvalues, const cpl_vector ysigmas, cpl_fit_mode fit_opts, int max_iterations, double x0, double x0_uncertainty, double sigma, double sigma_uncertainty, double fwhm, double fwhm_uncertainty, double area, double offset, double mse, double red_chisq, cpl_matrix *covariance, cpl_size lower_ndx, cpl_size *upper_ndx)
	Apply a Gaussian fitting until it is successful. More...

cpl_error_code	clipm_priv_optimize_downhill_simplex (cpl_vector fit_params, const cpl_vector start_delta, const cpl_vector precision,(cpl_vector params, void other), void other_func_par, int max_iterations, double optimum, int done_iterations)
	Optimize a set of parameters using the Downhill Simplex method by Nelder and Mead (also known as the "Amoeba" method). More...

Function Documentation

cpl_error_code clipm_priv_optimize_downhill_simplex	(	cpl_vector *	fit_params,
		const cpl_vector *	start_delta,
		const cpl_vector *	precision,
		(cpl_vector params,void other)	,
		void *	other_func_par,
		int	max_iterations,
		double *	optimum,
		int *	done_iterations
	)

Optimize a set of parameters using the Downhill Simplex method by Nelder and Mead (also known as the "Amoeba" method).

Parameters

fit_params	Vector containing the parameters to be optimized, which must be set to an initial guess, and will contain the result at the end
start_delta	Vector containing a set of parameter step sizes for the start of the optimization
precision	Optional vector containing respective precision limits (> 0) which must be met for convergence, can be NULL
evalfunc	Pointer to function which evaluates the parameters and returns a value which must be minimal in the optimum (for example an error value, or the inverse or negative of a matching), taking as first argument the parameter vector, and as second argument an arbitrary pointer (to provide any other data)
other_func_par	Arbitrary pointer used as the input for the second argument of evalfunc, for example a pointer to a dataset or a reference to a struct of pointers and further parameters
max_iterations	The maximum number of iterations
optimum	Optional output return value of evalfunc with optimized parameters, can be NULL
done_iterations	Optional output of number of performed iterations, can be NULL

Returns: CPL error code

Indication:

This function optimizes a set of paramaters, which can for example be used in a fitting process. The feature of the Downhill Simplex method is that it works without an analytical derivative of the evaluation function. If a derivative is available, the usage of cpl_fit_lvmq() (Levenberg-Marquard fitting) should be preferred.

Usage:

Two modes of operation are provided:

If precision is given, then fit_params is optimized until the variations of the parameters stay below their respective precision entry. This is defined as convergence. If max_iterations is reached before convergence, fit_params will contain the latest parameter values, and the error code CPL_ERROR_CONTINUE is set and returned.
If precision is NULL, then max_iterations iterations are performed, and fit_params contains the result at the end.

Evaluation Function:

If evalfunc sets a CPL error, this is interpreted as a real error, and the optimization is stopped. In this case, this error code is set and returned. fit_params contains the latest parameters.
If there are hard limits for the parameters, these must be implemented in the way that evalfunc returns a very bad (high) value, but does not set an error.

Note

Although this function works also with only one dimension (i.e. a parameter vector of size 1), this is not recommended since there are better optimization methods for the univariate case.
The number of optimization iterations does not equal the number of calls of evalfunc. Some iteration steps require several evaluations, i.e. the number of evalfunc calls is higher than the number of iterations.
done_iterations, if not NULL, returns the number of performed iterations. This is intended for instance for debugging purposes, or for the optimization of methods using this function.

Limitations:

It is possible that a local minimum is found. There are different approaches to solve this issue, one possibility is for example to restart the optimization with a reasonable high initial step size (start_delta), or to apply the optimization to a cluster of parameter guesses. But it should be tried out whether this is necessary at all, before taking the effort. This depends strongly on the data, and a general answer to this topic is hardly possible.

Error Handling:

If an error occurs, but some iteration could be made, then the latest parameter values are stored in fit_params.

The following error codes can be set and returned:

CPL_ERROR_NULL_INPUT: fit_params or start_delta or evalfunc is/are NULL
CPL_ERROR_INCOMPATIBLE_INPUT: fit_params and start_delta or precision (if not NULL) have different sizes
CPL_ERROR_ILLEGAL_INPUT: an entry of precision (if not NULL) is <= 0
CPL_ERROR_CONTINUE: the optimization could not converge to precision (see above)
Any error that is set by evalfunc

See Also: Nelder, J. A. & Mead, R. A: A Simplex Method for Function Minimization. Computer Journal, 1965, 7, 308-313

cpl_error_code clipm_priv_optimize_gaussian	(	const cpl_vector *	xvalues,
		const cpl_vector *	xsigmas,
		const cpl_vector *	yvalues,
		const cpl_vector *	ysigmas,
		cpl_fit_mode	fit_opts,
		int	max_iterations,
		double *	x0,
		double *	x0_uncertainty,
		double *	sigma,
		double *	sigma_uncertainty,
		double *	fwhm,
		double *	fwhm_uncertainty,
		double *	area,
		double *	offset,
		double *	mse,
		double *	red_chisq,
		cpl_matrix **	covariance,
		cpl_size *	lower_ndx,
		cpl_size *	upper_ndx
	)

Apply a Gaussian fitting until it is successful.

Parameters

xvalues	Positions to fit
xsigmas	Uncertainties of positions to fit (NOT supported)
yvalues	Values to fit
ysigmas	Uncertainties of values to fit
fit_opts	Fitting options: bitwise combination of CPL_FIT_CENTROID, CPL_FIT_STDEV, CPL_FIT_AREA, CPL_FIT_OFFSET or CPL_FIT_ALL (for all above)
max_iterations	Maximum number of iterations
x0	(output) Center of best fit gaussian. If CPL_FIT_CENTROID is not set, this is also an input parameter.
x0_uncertainty	(output) Optional uncertainty of x0.
sigma	(output) Width of best fit gaussian. A positive number on success. If CPL_FIT_STDEV is not set, this is also an input parameter.
sigma_uncertainty	(output) Optional uncertainty of sigma.
fwhm	(output) FWHM of best fit gaussian. A positive number on success.
fwhm_uncertainty	(output) Optional uncertainty of fwhm.
area	(output) Area of gaussian. A positive number on success. If CPL_FIT_AREA is not set, this is also an input parameter.
offset	(output) Fitted background level. If CPL_FIT_OFFSET is not set, this is also an input parameter.
mse	(output) If non-NULL, the mean squared error of the best fit is returned.
red_chisq	(output) If non-NULL, the reduced chi square of the best fit is returned. This requires the ysigmas to be specified.
covariance	(output) If non-NULL, the formal covariance matrix of the best fit is returned. This requires ysigmas to be specified. The order of fit parameters in the covariance matrix is defined as (x0, sigma, area, offset), for example the (3,3) element of the matrix (counting from zero) is the variance of the fitted offset. The matrix must be deallocated by calling cpl_matrix_delete() . On error, NULL is returned.
lower_ndx	Output lower index of iteration window, can be NULL
upper_ndx	Output upper index of iteration window, can be NULL

Returns: CPL error code

Principle:: A gaussian fitting is performed. If max_iterations > 0, a narrow window around the found peak is defined, and the fitting repeated in this data window. This is intended to exclude noise influences far away from the peak.

Constraints:: The vector xvalues must be sorted. If it is not, the iteration might fail.

Iteration Details:

Iteration window definition:
- If a gaussian fitting was successful, the next window is defined to cover the +/- 5 sigma region around the peak (using the determined sigma).
- If a gaussian fitting was not successful, the +/- 1 sigma region is used, using the centre and sigma guess values from cpl_vector_fit_gaussian(). This guessed sigma is usually much wider than the real one, that's why.
- The iteration window is always at least 5 data points wide.
The iteration is continued until either the window is stable or until max_iteration steps are performed.

Output Parameter Details:

The output parameters and returned error codes are defined to be the same as from function cpl_vector_fit_gaussian(), with the following extensions:

cpl_fit_gaussian() is called, and on failure a fitting window is defined to constrain the data and try another fitting, as long as the maximum number of trials max_iterations is not reached.
If max_iterations allows, on success a fitting is tried once more in a window defined as the +/- 5 sigma region, to improve the result.
xvalues does not need to be provided. If it is NULL, the indizes are taken as positions, starting with 0.

Error Handling:

The following error cases might happen, and the corresponding error code set:

CPL_ERROR_NULL_INPUT: yvalues or x0 or sigma or area or offset is NULL
CPL_ERROR_UNSUPPORTED_MODE:
- xsigmas != NULL
- the specified fit_opts is not a bitwise combination of the allowed values (e.g. 0 or 1)
CPL_ERROR_ILLEGAL_INPUT:
- max_iterations < 0
- the length of yvalues < 5
- any input noise values, sigma or area (if input) is non-positive
CPL_ERROR_INCOMPATIBLE_INPUT: any input vectors have different sizes
CPL_ERROR_CONTINUE: the optimization could not converge
CPL_ERROR_SINGULAR_MATRIX: the covariance matrix could not be calculated

Todo:

Enable cpl_vector_fit_gaussian() error history again when fixed [DFS05410].
Return area and offset uncertainties.
Fix iteration window for non-equidistant xvalues.

Functions

Function Documentation