NewFitterGSL Class Reference

A kinematic fitter using the Newton-Raphson method to solve the equations. More...

#include <NewFitterGSL.h>

Inheritance diagram for NewFitterGSL:

Public Types
enum	{ NPARMAX =50, NCONMAX =10, NUNMMAX =10 }
enum	{ NITMAX = 100 }

Public Member Functions
	NewFitterGSL ()
	Constructor.
virtual	~NewFitterGSL ()
	Virtual destructor.
virtual double	fit ()
	The fit method, returns the fit probability.
virtual int	getError () const
	Get the error code of the last fit: 0=OK, 1=failed.
virtual double	getProbability () const
	Get the fit probability of the last fit.
virtual double	getChi2 () const
	Get the chi**2 of the last fit.
virtual int	getDoF () const
	Get the number of degrees of freedom of the last fit.
virtual int	getIterations () const
	Get the number of iterations of the last fit.
virtual int	getNcon () const
	Get the number of hard constraints of the last fit.
virtual int	getNsoft () const
	Get the number of soft constraints of the last fit.
virtual int	getNpar () const
	Get the number of all parameters of the last fit.
virtual int	getNunm () const
	Get the number of unmeasured parameters of the last fit.
virtual bool	initialize ()
	Initialize the fitter.
virtual void	setDebug (int debuglevel)
	Set the Debug Level.
virtual void	determineLambdas (gsl_vector *vecxnew, const gsl_matrix *MatM, const gsl_vector *vecx, gsl_matrix *MatW, gsl_vector *vecw, double eps=0)
	Determine best lambda values. More...
virtual void	calc2ndOrderCorr (gsl_vector *vecdxhat, const gsl_vector *vecxnew, const gsl_matrix *MatM, gsl_matrix *MatW, gsl_vector *vecw, double eps=0)
	Calculate 2nd order correction step. More...
virtual gsl_matrix_view	calcZ (int &rankA, gsl_matrix *MatW1, gsl_matrix *MatW2, gsl_vector *vecw1, gsl_vector *vecw2, gsl_permutation *permW, double eps=0)
	Calculate null space of constraints; return value is a matrix view of MatW1, with column vectors spanning null(A) More...
virtual gsl_matrix_view	calcReducedHessian (int &rankH, gsl_matrix *MatW1, const gsl_vector *vecx, gsl_matrix *MatW2, gsl_matrix *MatW3, gsl_vector *vecw1, gsl_vector *vecw2, gsl_permutation *permW, double eps=0)
	Calculate reduced Hessian; return value is a matrix view of MatW1 giving the reduced Hessian. More...
virtual gsl_vector_view	calcReducedHessianEigenvalues (int &rankH, gsl_matrix *MatW1, const gsl_vector *vecx, gsl_matrix *MatW2, gsl_matrix *MatW3, gsl_vector *vecw1, gsl_vector *vecw2, gsl_permutation *permW, gsl_eigen_symm_workspace *eigenws, double eps=0)
virtual double	calcChi2 ()
	Calculate the chi2.
void	fillx (gsl_vector *vecx)
void	fillperr (gsl_vector *vece)
void	assembleM (gsl_matrix *MatM, const gsl_vector *vecx, bool errorpropagation=false)
void	assembleG (gsl_matrix *MatM, const gsl_vector *vecx)
void	scaleM (gsl_matrix *MatMscal, const gsl_matrix *MatM, const gsl_vector *vece)
void	assembley (gsl_vector *vecy, const gsl_vector *vecx)
void	scaley (gsl_vector *vecyscal, const gsl_vector *vecy, const gsl_vector *vece)
void	assembleChi2Der (gsl_vector *vecy)
void	addConstraints (gsl_vector *vecy)
void	assembleConstDer (gsl_matrix *MatM)
int	calcNewtonDx (gsl_vector *vecdx, gsl_vector *vecdxscal, gsl_vector *vecx, const gsl_vector *vece, gsl_matrix *MatM, gsl_matrix *MatMscal, gsl_vector *vecy, gsl_vector *vecyscal, gsl_matrix *MatW, gsl_matrix *MatW2, gsl_permutation *permW, gsl_vector *vecw)
int	calcLimitedDx (double &alpha, double &mu, gsl_vector *vecxnew, int imode, gsl_vector *vecx, gsl_vector *vecdxhat, const gsl_vector *vecdx, const gsl_vector *vecdxscal, const gsl_vector *vece, const gsl_matrix *MatM, const gsl_matrix *MatMscal, gsl_matrix *MatW, gsl_vector *vecw)
int	doLineSearch (double &alpha, gsl_vector *vecxnew, int imode, double phi0, double dphi0, double phiR, double eta, double zeta, double mu, const gsl_vector *vecx, const gsl_vector *vecdx, const gsl_vector *vece, gsl_vector *vecw)
double	calcMu (const gsl_vector *vecx, const gsl_vector *vece, const gsl_vector *vecdx, const gsl_vector *vecdxscal, const gsl_vector *xnew, const gsl_matrix *MatM, const gsl_matrix *MatMscal, gsl_vector *vecw)
double	meritFunction (double mu, const gsl_vector *vecx, const gsl_vector *vece)
double	meritFunctionDeriv (double mu, const gsl_vector *vecx, const gsl_vector *vece, const gsl_vector *vecdx, gsl_vector *vecw)
bool	updateParams (gsl_vector *vecx)
int	invertM ()
void	calcCovMatrix (gsl_matrix *MatW, gsl_permutation *permW, gsl_vector *vecx)
double	calcpTLp (const gsl_vector *vecdx, const gsl_matrix *MatM, gsl_vector *vecw)
int	solveSystem (gsl_vector *vecdxscal, double &detW, const gsl_vector *vecyscal, const gsl_matrix *MatMscal, gsl_matrix *MatW, gsl_matrix *MatW2, gsl_vector *vecw, double epsLU, double epsSV)
	solve system of equations Mscal*dxscal = yscal
int	solveSystemLU (gsl_vector *vecdxscal, double &detW, const gsl_vector *vecyscal, const gsl_matrix *MatMscal, gsl_matrix *MatW, gsl_vector *vecw, double eps)
	solve system of equations Mscal*dxscal = yscal using LU decomposition
int	solveSystemSVD (gsl_vector *vecdxscal, double &detW, const gsl_vector *vecyscal, const gsl_matrix *MatMscal, gsl_matrix *MatW, gsl_matrix *MatW2, gsl_vector *vecw, double eps)
	solve system of equations Mscal*dxscal = yscal using SVD decomposition
Public Member Functions inherited from BaseFitter
virtual void	addFitObject (BaseFitObject *fitobject_)
virtual void	addFitObject (BaseFitObject &fitobject_)
virtual void	addConstraint (BaseConstraint *constraint_)
virtual void	addConstraint (BaseConstraint &constraint_)
virtual void	addHardConstraint (BaseHardConstraint *constraint_)
virtual void	addHardConstraint (BaseHardConstraint &constraint_)
virtual void	addSoftConstraint (BaseSoftConstraint *constraint_)
virtual void	addSoftConstraint (BaseSoftConstraint &constraint_)
virtual std::vector < BaseFitObject * > *	getFitObjects ()
virtual std::vector < BaseHardConstraint * > *	getConstraints ()
virtual std::vector < BaseSoftConstraint * > *	getSoftConstraints ()
virtual void	reset ()
virtual BaseTracer *	getTracer ()
virtual const BaseTracer *	getTracer () const
virtual void	setTracer (BaseTracer *newTracer)
virtual void	setTracer (BaseTracer &newTracer)
virtual const double *	getGlobalCovarianceMatrix (int &idim) const
virtual double *	getGlobalCovarianceMatrix (int &idim)

Static Public Member Functions
static void	ini_gsl_permutation (gsl_permutation *&p, unsigned int size)
static void	ini_gsl_vector (gsl_vector *&v, int unsigned size)
static void	ini_gsl_matrix (gsl_matrix *&m, int unsigned size1, unsigned int size2)
static void	debug_print (const gsl_matrix *m, const char *name)
static void	debug_print (const gsl_vector *v, const char *name)
static void	add (gsl_vector *vecz, const gsl_vector *vecx, double a, const gsl_vector *vecy)
static bool	isfinite (const gsl_vector *vec)
static bool	isfinite (const gsl_matrix *mat)
static void	MoorePenroseInverse (gsl_matrix *Ainv, gsl_matrix *A, gsl_matrix *W, gsl_vector *w, double eps=0)
	Compute the Moore-Penrose pseudo-inverse A+ of A, using SVD. More...

Public Attributes
int	npar
	total number of parameters
int	ncon
	total number of hard constraints
int	nsoft
	total number of soft constraints
int	nunm
	total number of unmeasured parameters
int	ierr
	Error status.
int	nit
	Number of iterations.
double	fitprob
	fit probability
double	chi2
	final chi2
unsigned int	idim
gsl_vector *	x
gsl_vector *	xold
gsl_vector *	xnew
gsl_vector *	dx
gsl_vector *	dxscal
gsl_vector *	y
gsl_vector *	yscal
gsl_vector *	perr
gsl_vector *	v1
gsl_vector *	v2
gsl_matrix *	M
gsl_matrix *	Mscal
gsl_matrix *	W
gsl_matrix *	W2
gsl_matrix *	W3
gsl_matrix *	M1
gsl_matrix *	M2
gsl_matrix *	M3
gsl_matrix *	M4
gsl_matrix *	M5
gsl_matrix *	CC
gsl_matrix *	CC1
gsl_matrix *	CCinv
gsl_permutation *	permW
gsl_eigen_symm_workspace *	eigenws
unsigned int	eigenwsdim
double	chi2best
double	chi2new
double	chi2old
double	fvalbest
double	scale
double	scalebest
double	stepsize
double	stepbest
double	scalevals [NITMAX]
double	fvals [NITMAX]
int	imerit
bool	try2ndOrderCorr
int	debug
Public Attributes inherited from BaseFitter
std::map< std::string, double >	traceValues

Additional Inherited Members
Protected Types inherited from BaseFitter
typedef std::vector < BaseFitObject * >	FitObjectContainer
typedef std::vector < BaseHardConstraint * >	ConstraintContainer
typedef std::vector < BaseSoftConstraint * >	SoftConstraintContainer
typedef FitObjectContainer::iterator	FitObjectIterator
typedef ConstraintContainer::iterator	ConstraintIterator
typedef SoftConstraintContainer::iterator	SoftConstraintIterator
Protected Member Functions inherited from BaseFitter
	BaseFitter (const BaseFitter &rhs)
	Copy constructor disabled.
BaseFitter &	operator= (const BaseFitter &rhs)
	Assignment disabled.
Protected Attributes inherited from BaseFitter
FitObjectContainer	fitobjects
ConstraintContainer	constraints
SoftConstraintContainer	softconstraints
int	covDim
	dimension of global covariance matrix
double *	cov
	global covariance matrix of last fit problem
bool	covValid
	Flag whether global covariance is valid.
BaseTracer *	tracer

Detailed Description

A kinematic fitter using the Newton-Raphson method to solve the equations.

This class implements a kinematic fitter using the Newton-Raphson method to solve the system of equations arising from the Lagrange multiplier method

Author: Benno List Last update:

Date:

2011/05/03 13:16:41

by:

Author:

blist

Changelog:

• 15.11.2010 First version

Member Function Documentation

void NewFitterGSL::calc2ndOrderCorr ( gsl_vector *  vecdxhat, const gsl_vector *  vecxnew, const gsl_matrix *  MatM, gsl_matrix *  MatW, gsl_vector *  vecw, double  eps = 0  )

virtual

Calculate 2nd order correction step.

Parameters

vecdxhat	the correction step
vecxnew	the current state vector (parameters must be set to vecxnew!)
MatM	The matrix of the last Newton step
MatW	work matrix
vecw	work vector
eps	Singular values < eps*(max(abs(s_i))) are set to 0

References ncon, and npar.

Referenced by calcLimitedDx().

int NewFitterGSL::calcLimitedDx	(	double &	alpha,
		double &	mu,
		gsl_vector *	vecxnew,
		int	imode,
		gsl_vector *	vecx,
		gsl_vector *	vecdxhat,
		const gsl_vector *	vecdx,
		const gsl_vector *	vecdxscal,
		const gsl_vector *	vece,
		const gsl_matrix *	MatM,
		const gsl_matrix *	MatMscal,
		gsl_matrix *	MatW,
		gsl_vector *	vecw
	)

Parameters

alpha	Output value alpha
mu	Value of mu for merit function
vecxnew	New vector x
imode	mode: 0=Armijo, 1=Wolfe, 2=Goldstein
vecx	Current vector x
vecdxhat	the 2nd order correction step, if any
vecdx	Update vector dx
vecdxscal	Result: Update vector dx, scaled
vece	Error vector e
MatM	Matrix M
MatMscal	Matrix M, scaled
MatW	Work matrix
vecw	Work vector w1

References calc2ndOrderCorr(), calcChi2(), calcMu(), doLineSearch(), meritFunction(), meritFunctionDeriv(), and BaseTracer::substep().

Referenced by fit().

double NewFitterGSL::calcMu	(	const gsl_vector *	vecx,
		const gsl_vector *	vece,
		const gsl_vector *	vecdx,
		const gsl_vector *	vecdxscal,
		const gsl_vector *	xnew,
		const gsl_matrix *	MatM,
		const gsl_matrix *	MatMscal,
		gsl_vector *	vecw
	)

Parameters

vecx	Current vector x
vece	Current errors x
vecdx	Current step dx
vecdxscal	Current step dx, scaled
xnew	New vector x
MatM	Current matrix M
MatMscal	Current scaled matrix M
vecw	Work vector w

References calcpTLp(), ncon, and npar.

Referenced by calcLimitedDx().

int NewFitterGSL::calcNewtonDx	(	gsl_vector *	vecdx,
		gsl_vector *	vecdxscal,
		gsl_vector *	vecx,
		const gsl_vector *	vece,
		gsl_matrix *	MatM,
		gsl_matrix *	MatMscal,
		gsl_vector *	vecy,
		gsl_vector *	vecyscal,
		gsl_matrix *	MatW,
		gsl_matrix *	MatW2,
		gsl_permutation *	permW,
		gsl_vector *	vecw
	)

Parameters

vecdx	Result: Update vector dx
vecdxscal	Result: Update vector dx, scaled
vecx	Current vector x
vece	Current ``error'' set of x
MatM	Matrix M
MatMscal	Matrix M, scaled
vecy	Vector y
vecyscal	Vector y, scaled,
MatW	Work matrix
MatW2	Work matrix
permW	Work permutation vector
vecw	Work vector

References calcpTLp(), determineLambdas(), ncon, npar, and solveSystem().

Referenced by fit().

double NewFitterGSL::calcpTLp	(	const gsl_vector *	vecdx,
		const gsl_matrix *	MatM,
		gsl_vector *	vecw
	)

Parameters

vecdx	Current step dx
MatM	Current matrix M
vecw	Work vector w

References npar.

Referenced by calcMu(), and calcNewtonDx().

gsl_matrix_view NewFitterGSL::calcReducedHessian ( int &  rankH, gsl_matrix *  MatW1, const gsl_vector *  vecx, gsl_matrix *  MatW2, gsl_matrix *  MatW3, gsl_vector *  vecw1, gsl_vector *  vecw2, gsl_permutation *  permW, double  eps = 0  )

virtual

Calculate reduced Hessian; return value is a matrix view of MatW1 giving the reduced Hessian.

Parameters

rankH	dimension of H
MatW1	work matrix, contains H at the end
vecx	vector with current x values
MatW2	work matrix, contains Z at the end
MatW3	work matrix
vecw1	work vector
vecw2	work vector
permW	Work permutation vector
eps	Singular values < eps*(max(abs(s_i))) are set to 0

References calcZ(), and npar.

Referenced by calcReducedHessianEigenvalues().

gsl_vector_view NewFitterGSL::calcReducedHessianEigenvalues ( int &  rankH, gsl_matrix *  MatW1, const gsl_vector *  vecx, gsl_matrix *  MatW2, gsl_matrix *  MatW3, gsl_vector *  vecw1, gsl_vector *  vecw2, gsl_permutation *  permW, gsl_eigen_symm_workspace *  eigenws, double  eps = 0  )

virtual

Parameters

rankH	dimension of H
MatW1	work matrix, contains H at the end
vecx	vector with current x values
MatW2	work matrix, contains Z at the end
MatW3	work matrix
vecw1	work vector
vecw2	work vector
permW	Work permutation vector
eigenws	Work space for eigenvalue calculation
eps	Singular values < eps*(max(abs(s_i))) are set to 0

References calcReducedHessian().

gsl_matrix_view NewFitterGSL::calcZ ( int &  rankA, gsl_matrix *  MatW1, gsl_matrix *  MatW2, gsl_vector *  vecw1, gsl_vector *  vecw2, gsl_permutation *  permW, double  eps = 0  )

virtual

Calculate null space of constraints; return value is a matrix view of MatW1, with column vectors spanning null(A)

Parameters

rankA	rank of A (= number of lin. indep. constraints)
MatW1	work matrix, contains Z at the end
MatW2	work matrix
vecw1	work vector
vecw2	work vector
permW	Work permutation vector
eps	Singular values < eps*(max(abs(s_i))) are set to 0

References ncon, and npar.

Referenced by calcReducedHessian().

void NewFitterGSL::determineLambdas ( gsl_vector *  vecxnew, const gsl_matrix *  MatM, const gsl_vector *  vecx, gsl_matrix *  MatW, gsl_vector *  vecw, double  eps = 0  )

virtual

Determine best lambda values.

Parameters

vecxnew	vector with new lambda values
MatM	matrix with constraint derivatives
vecx	vector with current x values
MatW	work matrix
vecw	work vector
eps	Singular values < eps*(max(abs(s_i))) are set to 0

References ncon, and npar.

Referenced by calcNewtonDx(), and fit().

int NewFitterGSL::doLineSearch	(	double &	alpha,
		gsl_vector *	vecxnew,
		int	imode,
		double	phi0,
		double	dphi0,
		double	phiR,
		double	eta,
		double	zeta,
		double	mu,
		const gsl_vector *	vecx,
		const gsl_vector *	vecdx,
		const gsl_vector *	vece,
		gsl_vector *	vecw
	)

Parameters

alpha	Output value alpha
vecxnew	New vector x
imode	mode: 0=Armijo, 1=Wolfe, 2=Goldstein
phi0	Merit function for alpha=0
dphi0	Directional derivative of merit function for alpha=0
phiR	Merit function for given alpha
eta	Constant for Armijo's rule
zeta	Constant for Wolfe's or Goldstein's rule
mu	Value of mu for merit function
vecx	Current vector x
vecdx	Update vector dx
vece	Error vector e
vecw	Work vector w

References meritFunction(), meritFunctionDeriv(), nit, and BaseTracer::substep().

Referenced by calcLimitedDx().

double NewFitterGSL::meritFunction	(	double	mu,
		const gsl_vector *	vecx,
		const gsl_vector *	vece
	)

Parameters

mu	Value of mu
vecx	Current vector x
vece	Current errors x

References calcChi2(), BaseHardConstraint::getGlobalNum(), and BaseHardConstraint::getValue().

Referenced by calcLimitedDx(), and doLineSearch().

double NewFitterGSL::meritFunctionDeriv	(	double	mu,
		const gsl_vector *	vecx,
		const gsl_vector *	vece,
		const gsl_vector *	vecdx,
		gsl_vector *	vecw
	)

Parameters

mu	Value of mu
vecx	Current vector x
vece	Current errors x
vecdx	Current update vector dx
vecw	work vector

References BaseHardConstraint::getGlobalNum(), BaseHardConstraint::getValue(), and npar.

Referenced by calcLimitedDx(), and doLineSearch().

void NewFitterGSL::MoorePenroseInverse ( gsl_matrix *  Ainv, gsl_matrix *  A, gsl_matrix *  W, gsl_vector *  w, double  eps = 0  )

static

Compute the Moore-Penrose pseudo-inverse A+ of A, using SVD.

Parameters

Ainv	Result: m x n matrix A+
A	Input: n x m matrix A, n >= m (is destroyed!)
W	Work matrix, at least m x m
w	Work vector w, at least m
eps	Singular values < eps*(max(abs(s_i))) are set to 0

---

The documentation for this class was generated from the following files:

• NewFitterGSL.h
• NewFitterGSL.cc