foxleg.cpp File Reference

#include "geoslib_old_f.h"
#include "Matrix/MatrixSquareSymmetric.hpp"
#include "Basic/VectorHelper.hpp"
#include "Basic/Utilities.hpp"
#include "Basic/String.hpp"
#include "Basic/OptDbg.hpp"
#include <math.h>

Functions
static void	st_gradient (VectorDouble &param, VectorDouble &lower, VectorDouble &upper, VectorDouble &scale, VectorDouble &tabwgt, MatrixRectangular &Jr, VectorDouble &param1, VectorDouble &param2, VectorDouble &tabmod1, VectorDouble &tabmod2)
static double	st_residuals (VectorDouble &param, VectorDouble &tabexp, VectorDouble &tabwgt, VectorDouble &tabmod, VectorDouble &residuals)
static void	st_determine_gauss (MatrixRectangular &Jr, MatrixSquareGeneral &gauss)
static double	st_norm_hgn (VectorDouble &hgn, VectorDouble &scale)
static double	st_essai (VectorDouble &hgnadm, VectorDouble &grad_red, MatrixSquareGeneral &gauss_red)
static int	st_solve_hgnc (int npar, const VectorDouble &grad, const MatrixSquareGeneral &gauss, VectorDouble &hgnc, int flaginvsign)
static void	st_fill_constraints (const MatrixRectangular &acont, VectorDouble &grad, MatrixSquareGeneral &gauss)
static int	st_calcul0 (VectorDouble &param, VectorDouble &lower, VectorDouble &upper, VectorDouble &scale, const MatrixRectangular &acont, VectorDouble &tabwgt, VectorDouble &residuals, MatrixRectangular &Jr, VectorDouble &grad, MatrixSquareGeneral &gauss, VectorDouble &hgnc, VectorDouble &param1, VectorDouble &param2, VectorDouble &tabmod1, VectorDouble &tabmod2)
static int	st_possibilities (int npar, MatrixRectangular &bords, VectorDouble &ai, VectorDouble &hgnc, VectorInt &flag, VectorDouble &temp)
static int	st_define_constraints (int mode, MatrixRectangular &bords_red, VectorDouble &ai_red, VectorDouble &hgnc, MatrixRectangular &consts, VectorInt &flag, VectorDouble &temp)
static void	st_minimum (VectorInt &, VectorInt &flag, MatrixRectangular &bords_red, const VectorDouble &top, const VectorDouble &bot, VectorDouble &hgnc, VectorDouble &hgnadm)
static void	st_update_bords (MatrixRectangular &bords, VectorInt &ind_util, MatrixRectangular &bords_red)
static int	st_suppress_unused_constraints (MatrixRectangular &bords, VectorDouble &ai, VectorDouble &grad, MatrixSquareGeneral &gauss, VectorDouble &hgnc, VectorInt &ind_util, MatrixRectangular &bords_red, VectorDouble &ai_red, VectorDouble &grad_red, MatrixSquareGeneral &gauss_red, VectorInt &flag1, VectorInt &flag2, VectorDouble &temp)
static int	st_establish_minimization (int nactive, VectorInt &ind_util, VectorInt &flag_active, MatrixRectangular &bords_red, VectorDouble &ai_red, VectorDouble &grad_red, MatrixSquareGeneral &gauss_red, int *lambda_neg, VectorDouble &hgnc, MatrixSquareGeneral &a, VectorDouble &b, VectorDouble &temp)
static void	st_check (VectorInt &ind_util, VectorDouble &hgnc, const MatrixRectangular &acont)
static int	st_minimization_under_constraints (VectorInt &ind_util, MatrixRectangular &bords_red, VectorDouble &ai_red, VectorDouble &grad_red, MatrixSquareGeneral &gauss_red, MatrixRectangular &consts, VectorDouble &hgnc, VectorDouble &hgnadm, VectorInt &flag_active, VectorInt &flag_actaux, MatrixSquareGeneral &a, VectorDouble &b1, VectorDouble &b2, VectorDouble &b3, VectorDouble &temp, const MatrixRectangular &acont)
static void	st_constraints_init (VectorInt &ind_util, VectorDouble &ai)
static void	st_define_bounds (VectorDouble &param, VectorDouble &lower, VectorDouble &upper, VectorDouble &scale, double delta, MatrixRectangular &bords)
static void	st_foxleg_debug_title (void)
static void	st_foxleg_debug_current (double mscur, double delta, VectorDouble &param)
static void	st_foxleg_score (const Option_AutoFit &mauto, double mscur, double delta, double arret)
static void	st_linear_interpolate (double mscur, VectorDouble &param, const MatrixRectangular &acont, VectorDouble &tabexp, VectorDouble &tabwgt, MatrixRectangular &bords, VectorDouble &grad, double *msaux, VectorDouble &paramaux, VectorDouble &residuals, VectorDouble &tabmod1)
static int	st_check_param (VectorDouble &param, VectorDouble &lower, VectorDouble &upper)
int	foxleg_f (int ndat, int npar, int ncont, const MatrixRectangular &acont, VectorDouble &param, VectorDouble &lower, VectorDouble &upper, VectorDouble &scale, const Option_AutoFit &mauto, int flag_title, void(*func_evaluate)(int ndat, int npar, VectorDouble &param, VectorDouble &work), VectorDouble &tabexp, VectorDouble &tabwgt)
int	add_sill_constraints (Constraints &constraints, double constantSill)
int	add_unit_sill_constraints (Constraints &constraints)

Variables
static int	VERBOSE_GQO = 0
static int	NPAR
static int	NPAR2
static int	NPARAC
static int	NPARAC2
static int	NDAT
static int	NCONT
static int	NPCT
static int	NPCT2
static int	ITERATION
static int	SOUSITER
static void(*	FUNC_EVALUATE )(int ndat, int npar, VectorDouble &param, VectorDouble &work)

Function Documentation

add_sill_constraints()

int add_sill_constraints	(	Constraints &	constraints,
		double	constantSill
	)

Add constraints to the Option_AutoFit structure

Returns

Error return code

Parameters

[in]	constraints	Constraints structure
[in]	constantSill	Constant value for the Sill as a constraint

add_unit_sill_constraints()

int add_unit_sill_constraints

(

Constraints &

constraints

)

Add constraints (all equal to 1) to the Option_AutoFit structure

Returns

Error return code

Parameters

[in]

constraints

Constraints structure

foxleg_f()

int foxleg_f	(	int	ndat,
		int	npar,
		int	ncont,
		const MatrixRectangular &	acont,
		VectorDouble &	param,
		VectorDouble &	lower,
		VectorDouble &	upper,
		VectorDouble &	scale,
		const Option_AutoFit &	mauto,
		int	flag_title,
		void(*)(int ndat, int npar, VectorDouble &param, VectorDouble &work)	func_evaluate,
		VectorDouble &	tabexp,
		VectorDouble &	tabwgt
	)

Foxleg algorithm

Returns

Error returned code

0 : Correct termination

-1 : The convergence has not been reached

1 : Core problem

Parameters

[in]	ndat	Number of control points
[in]	npar	Number of parameters to estimate
[in]	ncont	Number of additional constraints
[in]	acont	Matrix of additional constraints (Dimension = ncont * npar)
[in]	param	Current values of the parameters
[in]	lower	Array of lower values
[in]	upper	Array of upper values
[in]	scale	Array of scale
[in]	mauto	Option_AutoFit structure
[in]	flag_title	Print the title after func_evaluate()
[in]	func_evaluate	Function for evaluating the model
[in]	tabexp	Array of values at control points
[in]	tabwgt	Array of weights at control points

Remarks

The evaluation function func_evaluate() is called with the

following arguments:

ndat : Number of control points

npar : Number of parameters

param : Vector of current values for the parameters

work : Output vector for the values at control points

evaluated using the current parameter values

Some additional constraints may be applied on the parameters

When not used, we must set: ncont=0, acont=empty()

st_calcul0()

static int st_calcul0 ( VectorDouble &  param, VectorDouble &  lower, VectorDouble &  upper, VectorDouble &  scale, const MatrixRectangular &  acont, VectorDouble &  tabwgt, VectorDouble &  residuals, MatrixRectangular &  Jr, VectorDouble &  grad, MatrixSquareGeneral &  gauss, VectorDouble &  hgnc, VectorDouble &  param1, VectorDouble &  param2, VectorDouble &  tabmod1, VectorDouble &  tabmod2  )

static

Evaluate Gradient and Hermitian matrices

Parameters

[in]	param	Current values of the parameters
[in]	lower	Array of lower values
[in]	upper	Array of upper values
[in]	scale	Array of scaling values
[in]	acont	Array of constraints
[in]	tabwgt	Array of weights at control points
[out]	residuals	Array of residuals
[out]	Jr	Array of gradients
[out]	grad	Gradient matrix
[out]	gauss	Gauss matrix
[out]	hgnc	Resulting hgnc array
[out]	param1	Working array (Dimension: NPAR)
[out]	param2	Working array (Dimension: NPAR)
[out]	tabmod1	Working array (Dimension: NDAT)
[out]	tabmod2	Working array (Dimension: NDAT)

st_check()

static void st_check ( VectorInt &  ind_util, VectorDouble &  hgnc, const MatrixRectangular &  acont  )

static

Check that the algorithm is valid

Parameters

[in]	ind_util	List of retained constraint indices
[in]	hgnc	Working vector
[in]	acont	Matrix of additional constraints

st_check_param()

static int st_check_param ( VectorDouble &  param, VectorDouble &  lower, VectorDouble &  upper  )

static

Check if the default, lower and upper bounds of the parameters are consistent or not

Returns

Error returned code

Parameters

[in]	param	Current values of the parameters
[in]	lower	Array of lower values
[in]	upper	Array of upper values

st_constraints_init()

static void st_constraints_init ( VectorInt &  ind_util, VectorDouble &  ai  )

static

Initialize the constraints

Parameters

[in]	ind_util	List of retained constraint indices
[in]	ai	AI matrix

st_define_bounds()

static void st_define_bounds ( VectorDouble &  param, VectorDouble &  lower, VectorDouble &  upper, VectorDouble &  scale, double  delta, MatrixRectangular &  bords  )

static

Evaluate the bounds and check if one constraint is on its boudn

Parameters

[in]	param	Current values of the parameters
[in]	lower	Array of lower values
[in]	upper	Array of upper values
[in]	scale	Array of scaling values
[in]	delta	Radius of the trusting area
[out]	bords	Value for the bounds

st_define_constraints()

static int st_define_constraints ( int  mode, MatrixRectangular &  bords_red, VectorDouble &  ai_red, VectorDouble &  hgnc, MatrixRectangular &  consts, VectorInt &  flag, VectorDouble &  temp  )

static

Calculate the number of constraints

Parameters

[in]	mode	Type of constraints -1 : the constraints which are non positive 0 : the constraints must be zero 1 : the constraints must be zero (cumul)
[in]	bords_red	Reduced array containing the bounds
[in]	ai_red	Reduced AI matrix
[in]	hgnc	Resulting hgnc array
[out]	consts	Array of constraints
[out]	flag	Array of indices with zero valid constraint
[out]	temp	Working array

st_determine_gauss()

static void st_determine_gauss ( MatrixRectangular &  Jr, MatrixSquareGeneral &  gauss  )

static

Calculate the Gauss matrix

Parameters

[in]	Jr	Array of gradients
[out]	gauss	Gaussian matrix

st_essai()

static double st_essai ( VectorDouble &  hgnadm, VectorDouble &  grad_red, MatrixSquareGeneral &  gauss_red  )

static

Score of the Minimization under constraints

Parameters

[in]	hgnadm	Admissible vector
[in]	grad_red	Reduced Gradient matrix
[in]	gauss_red	Reduced Gauss matrix

st_establish_minimization()

static int st_establish_minimization ( int  nactive, VectorInt &  ind_util, VectorInt &  flag_active, MatrixRectangular &  bords_red, VectorDouble &  ai_red, VectorDouble &  grad_red, MatrixSquareGeneral &  gauss_red, int *  lambda_neg, VectorDouble &  hgnc, MatrixSquareGeneral &  a, VectorDouble &  b, VectorDouble &  temp  )

static

Minimization under constraints

Returns

Error returned code

Parameters

[in]	nactive	Number of active constraints
[in]	ind_util	List of retained constraint indices
[in]	flag_active	Array of indices with zero valid constraint
[in]	bords_red	Reduced array containing the bounds
[in]	ai_red	Reduced AI matrix
[in]	grad_red	Reduced Gradient matrix
[in]	gauss_red	Reduced Gauss matrix
[out]	lambda_neg	Index of the first negative lambda value
[out]	hgnc	Resulting Hgnc array
[out]	a	Minimization L.H.S. matrix
[out]	b	Minimization R.H.S. matrix
[out]	temp	Working array

st_fill_constraints()

static void st_fill_constraints ( const MatrixRectangular &  acont, VectorDouble &  grad, MatrixSquareGeneral &  gauss  )

static

Add ncont linear constraints (AX=B) to the linear system (of size npar) used for quadratic optimization

Parameters

[in]	acont	Matrix A
[out]	grad	left hand-side of the system
[out]	gauss	matrix of the system

st_foxleg_debug_current()

static void st_foxleg_debug_current ( double  mscur, double  delta, VectorDouble &  param  )

static

Display the current status for FOXLEG trace

st_foxleg_debug_title()

static void st_foxleg_debug_title ( void  )

static

Display the title for FOXLEG trace

st_foxleg_score()

static void st_foxleg_score ( const Option_AutoFit &  mauto, double  mscur, double  delta, double  arret  )

static

Display the FOXLEG score

st_gradient()

static void st_gradient ( VectorDouble &  param, VectorDouble &  lower, VectorDouble &  upper, VectorDouble &  scale, VectorDouble &  tabwgt, MatrixRectangular &  Jr, VectorDouble &  param1, VectorDouble &  param2, VectorDouble &  tabmod1, VectorDouble &  tabmod2  )

static

Calculate the gradient

Parameters

[in]	param	Current values of the parameters
[in]	lower	Array of lower values
[in]	upper	Array of upper values
[in]	scale	Array of scaling values
[in]	tabwgt	Array of weights
[out]	Jr	Array of gradients
[out]	param1	Working array (Dimension: NPAR)
[out]	param2	Working array (Dimension: NPAR)
[out]	tabmod1	Working array (Dimension: NDAT)
[out]	tabmod2	Working array (Dimension: NDAT)

st_linear_interpolate()

static void st_linear_interpolate ( double  mscur, VectorDouble &  param, const MatrixRectangular &  acont, VectorDouble &  tabexp, VectorDouble &  tabwgt, MatrixRectangular &  bords, VectorDouble &  grad, double *  msaux, VectorDouble &  paramaux, VectorDouble &  residuals, VectorDouble &  tabmod1  )

static

Interpolate linearly the vector of parameters

Parameters

[in]	mscur	Current minimization value
[in]	param	Current values of the parameters
[in]	acont	Array of constraints
[in]	tabexp	Array of values at control points
[in]	tabwgt	Array of weights at control points
[in]	bords	Array containing the bounds
[in]	grad	Gradient matrix
[out]	msaux	New minimization value
[out]	paramaux	New vector of parameters
[out]	residuals	Array of residuals
[out]	tabmod1	Working array (Dimension: NDAT)

st_minimization_under_constraints()

static int st_minimization_under_constraints ( VectorInt &  ind_util, MatrixRectangular &  bords_red, VectorDouble &  ai_red, VectorDouble &  grad_red, MatrixSquareGeneral &  gauss_red, MatrixRectangular &  consts, VectorDouble &  hgnc, VectorDouble &  hgnadm, VectorInt &  flag_active, VectorInt &  flag_actaux, MatrixSquareGeneral &  a, VectorDouble &  b1, VectorDouble &  b2, VectorDouble &  b3, VectorDouble &  temp, const MatrixRectangular &  acont  )

static

Minimization under constraints

Returns

1 if the optimum has not been reached (various reasons)

Parameters

[in]	ind_util	List of retained constraint indices
[in]	bords_red	Reduced array containing the bounds
[in]	ai_red	Reduced AI matrix
[in]	grad_red	Reduced Gradient matrix
[in]	gauss_red	Reduced Gauss matrix
[out]	consts	Array of constraints
[out]	hgnc	Working vector
[out]	hgnadm	Working array
[out]	flag_active	Array of indices with zero valid constraint
[out]	flag_actaux	Array of indices with negative valid constraint
[out]	a	Minimization L.H.S. matrix
[out]	b1	Minimization R.H.S. matrix
[out]	b2	Minimization R.H.S. matrix
[out]	b3	Minimization R.H.S. matrix
[out]	temp	Working array
[out]	acont	Constraint array

st_minimum()

static void st_minimum ( VectorInt &  , VectorInt &  flag, MatrixRectangular &  bords_red, const VectorDouble &  top, const VectorDouble &  bot, VectorDouble &  hgnc, VectorDouble &  hgnadm  )

static

Calculate the minimum of the criterion and update hgnadm

Parameters

[in]	flag	Array of indices with negative valid constraint
[in]	bords_red	Reduced array containing the bounds
[in]	top	Calculated top of the fraction
[in]	bot	Calculated bottom of the fraction
[in]	hgnc	Hgnc array
[out]	hgnadm	Admissible Hgn array

st_norm_hgn()

static double st_norm_hgn ( VectorDouble &  hgn, VectorDouble &  scale  )

static

Calculate the Norm of the HGN vector

Returns

The norm value

Parameters

[in]	hgn	Working vector
[in]	scale	Scaling values

st_possibilities()

static int st_possibilities ( int  npar, MatrixRectangular &  bords, VectorDouble &  ai, VectorDouble &  hgnc, VectorInt &  flag, VectorDouble &  temp  )

static

Calculate the number of new inactive constraints

Returns

Number of highlited constraints

Parameters

[in]	npar	Current number of parameters
[in]	bords	Array containing the bounds
[in]	ai	AI matrix
[in]	hgnc	Resulting hgnc array
[out]	flag	Working array
[out]	temp	Working array

st_residuals()

static double st_residuals ( VectorDouble &  param, VectorDouble &  tabexp, VectorDouble &  tabwgt, VectorDouble &  tabmod, VectorDouble &  residuals  )

static

Calculate the residuals between the model and the experimental

Returns

The weighted mean squared error

Parameters

[in]	param	Current values of the parameters
[in]	tabexp	Array of values at control points
[in]	tabwgt	Array of weights at control points
[out]	tabmod	Working array (Dimension: NDAT)
[out]	residuals	Array of residuals

st_solve_hgnc()

static int st_solve_hgnc ( int  npar, const VectorDouble &  grad, const MatrixSquareGeneral &  gauss, VectorDouble &  hgnc, int  flaginvsign  )

static

Solve the direct minimization problem

Returns

Error return code

Parameters

[in]	npar	Current number of parameters
[in]	grad	Gradient matrix
[in]	gauss	Gauss matrix
[in]	hgnc	Resulting hgnc array
[in]	flaginvsign	if 1, the result is multiplied by -1

st_suppress_unused_constraints()

static int st_suppress_unused_constraints ( MatrixRectangular &  bords, VectorDouble &  ai, VectorDouble &  grad, MatrixSquareGeneral &  gauss, VectorDouble &  hgnc, VectorInt &  ind_util, MatrixRectangular &  bords_red, VectorDouble &  ai_red, VectorDouble &  grad_red, MatrixSquareGeneral &  gauss_red, VectorInt &  flag1, VectorInt &  flag2, VectorDouble &  temp  )

static

Eliminate the useless constraints

Returns

1 if the number of possibilities is reduced down to zero

or when the Generalized Inverse calculation failed

Parameters

[in]	bords	Array containing the bounds
[in]	ai	AI matrix
[in]	grad	Gradient matrix
[in]	gauss	Gaussian matrix
[in]	hgnc	hgnc array
[out]	ind_util	List of retained constraint indices
[out]	bords_red	Reduced Bounds array
[out]	ai_red	Reduced AI matrix
[out]	grad_red	Reduced Gradient matrix
[out]	gauss_red	Reduced Gauss matrix
[out]	flag1	Working array
[out]	flag2	Working array
[out]	temp	Working array

st_update_bords()

static void st_update_bords ( MatrixRectangular &  bords, VectorInt &  ind_util, MatrixRectangular &  bords_red  )

static

Update the constraints when no move has been performed

Parameters

[in]	bords	Array containing the bounds
[in]	ind_util	List of retained constraint indices
[out]	bords_red	Reduced Bounds array

Variable Documentation

FUNC_EVALUATE

void(* FUNC_EVALUATE) (int ndat, int npar, VectorDouble &param, VectorDouble &work) ( int  ndat, int  npar, VectorDouble &  param, VectorDouble &  work  )

static

ITERATION

int ITERATION

static

NCONT

int NCONT

static

NDAT

int NDAT

static

NPAR

int NPAR

static

NPAR2

int NPAR2

static

NPARAC

int NPARAC

static

NPARAC2

int NPARAC2

static

NPCT

int NPCT

static

NPCT2

int NPCT2

static

SOUSITER

int SOUSITER

static

VERBOSE_GQO

int VERBOSE_GQO = 0

static