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

Functions
static void	st_gradient (VectorDouble &param, VectorDouble &lower, VectorDouble &upper, VectorDouble &scale, VectorDouble &tabwgt, MatrixDense &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 (MatrixDense &Jr, MatrixSquare &gauss)

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

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

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

static void	st_fill_constraints (const MatrixDense &acont, VectorDouble &grad, MatrixSquare &gauss)

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

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

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

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

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

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

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

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

static int	st_minimization_under_constraints (VectorInt &ind_util, MatrixDense &bords_red, VectorDouble &ai_red, VectorDouble &grad_red, MatrixSquare &gauss_red, MatrixDense &consts, VectorDouble &hgnc, VectorDouble &hgnadm, VectorInt &flag_active, VectorInt &flag_actaux, MatrixSquare &a, VectorDouble &b1, VectorDouble &b2, VectorDouble &b3, VectorDouble &temp, const MatrixDense &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, MatrixDense &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 MatrixDense &acont, VectorDouble &tabexp, VectorDouble &tabwgt, MatrixDense &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 MatrixDense &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)

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

◆ foxleg_f()

int foxleg_f	(	int	ndat,
		int	npar,
		int	ncont,
		const MatrixDense &	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 MatrixDense &	acont,
		VectorDouble &	tabwgt,
		VectorDouble &	residuals,
		MatrixDense &	Jr,
		VectorDouble &	grad,
		MatrixSquare &	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 MatrixDense &	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,
		MatrixDense &	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,
		MatrixDense &	bords_red,
		VectorDouble &	ai_red,
		VectorDouble &	hgnc,
		MatrixDense &	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	(	MatrixDense &	Jr,
		MatrixSquare &	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,
		MatrixSquare &	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,
		MatrixDense &	bords_red,
		VectorDouble &	ai_red,
		VectorDouble &	grad_red,
		MatrixSquare &	gauss_red,
		int *	lambda_neg,
		VectorDouble &	hgnc,
		MatrixSquare &	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 MatrixDense &	acont,
		VectorDouble &	grad,
		MatrixSquare &	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,
		MatrixDense &	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 MatrixDense &	acont,
		VectorDouble &	tabexp,
		VectorDouble &	tabwgt,
		MatrixDense &	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,
		MatrixDense &	bords_red,
		VectorDouble &	ai_red,
		VectorDouble &	grad_red,
		MatrixSquare &	gauss_red,
		MatrixDense &	consts,
		VectorDouble &	hgnc,
		VectorDouble &	hgnadm,
		VectorInt &	flag_active,
		VectorInt &	flag_actaux,
		MatrixSquare &	a,
		VectorDouble &	b1,
		VectorDouble &	b2,
		VectorDouble &	b3,
		VectorDouble &	temp,
		const MatrixDense &	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,
		MatrixDense &	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,
		MatrixDense &	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 MatrixSquare &	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	(	MatrixDense &	bords,
		VectorDouble &	ai,
		VectorDouble &	grad,
		MatrixSquare &	gauss,
		VectorDouble &	hgnc,
		VectorInt &	ind_util,
		MatrixDense &	bords_red,
		VectorDouble &	ai_red,
		VectorDouble &	grad_red,
		MatrixSquare &	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	(	MatrixDense &	bords,
		VectorInt &	ind_util,
		MatrixDense &	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

Functions

Variables

Function Documentation

◆ foxleg_f()

◆ st_calcul0()

◆ st_check()

◆ st_check_param()

◆ st_constraints_init()

◆ st_define_bounds()

◆ st_define_constraints()

◆ st_determine_gauss()

◆ st_essai()

◆ st_establish_minimization()

◆ st_fill_constraints()

◆ st_foxleg_debug_current()

◆ st_foxleg_debug_title()

◆ st_foxleg_score()

◆ st_gradient()

◆ st_linear_interpolate()

◆ st_minimization_under_constraints()

◆ st_minimum()

◆ st_norm_hgn()

◆ st_possibilities()

◆ st_residuals()

◆ st_solve_hgnc()

◆ st_suppress_unused_constraints()

◆ st_update_bords()

Variable Documentation

◆ FUNC_EVALUATE

◆ ITERATION

◆ NCONT

◆ NDAT

◆ NPAR

◆ NPAR2

◆ NPARAC

◆ NPARAC2

◆ NPCT

◆ NPCT2

◆ SOUSITER

◆ VERBOSE_GQO