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

#include <NewFitterGSL.h>

Inheritance diagram for NewFitterGSL:

Collaboration diagram for NewFitterGSL:

[legend]

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

enum	{ NITMAX = 100 }

Public Member Functions
	NewFitterGSL ()
	Constructor.

virtual	~NewFitterGSL ()
	Virtual destructor.

virtual double	fit () override
	The fit method, returns the fit probability.

virtual int	getError () const override
	Get the error code of the last fit: 0=OK, 1=failed.

virtual double	getProbability () const override
	Get the fit probability of the last fit.

virtual double	getChi2 () const override
	Get the chi**2 of the last fit.

virtual int	getDoF () const override
	Get the number of degrees of freedom of the last fit.

virtual int	getIterations () const override
	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 () override
	Initialize the fitter.

virtual void	setDebug (int debuglevel)
	Set the Debug Level.

virtual int	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)

int	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 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 ()

int	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, 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

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

std::map< std::string, double >	traceValues

Protected Types
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 Attributes
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

Definition at line 51 of file NewFitterGSL.h.

Member Function Documentation

◆ calc2ndOrderCorr()

void 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

Definition at line 1633 of file NewFitterGSL.cc.

     {
       assert(vecdxhat);
       assert(vecdxhat->size == idim);
       assert(vecxnew);
       assert(vecxnew->size == idim);
       assert(MatM);
       assert(MatM->size1 == idim && MatM->size2 == idim);
       assert(MatW);
       assert(MatW->size1 == idim && MatW->size2 == idim);
       assert(vecw);
       assert(vecw->size == idim);
       assert(idim == static_cast<unsigned int>(npar + ncon));
  
  
       // Calculate 2nd order correction
       // see Nocedal&Wright (15.36)
       gsl_matrix_const_view AT(gsl_matrix_const_submatrix(MatM, 0,    npar, npar, ncon));
       gsl_matrix_view AAT(gsl_matrix_submatrix(MatW, npar, npar, ncon, ncon));
       gsl_matrix_view ATcopy(gsl_matrix_submatrix(MatW, 0,    npar, npar, ncon));
       gsl_matrix_view& V = AAT;
  
       gsl_vector_set_zero(vecdxhat);
       addConstraints(vecdxhat);
       gsl_vector_view c(gsl_vector_subvector(vecdxhat, npar, ncon));
       gsl_vector_view phat = (gsl_vector_subvector(vecdxhat, 0, npar));
  
       gsl_vector_set_zero(vecw);
       gsl_vector_view AATinvc(gsl_vector_subvector(vecw, npar, ncon));
       gsl_vector_view& s = AATinvc;
  
  
       // AAT = 1*A*A^T + 0*AAT
       gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1, &AT.matrix, &AT.matrix, 0, &AAT.matrix);
  
       // solve AAT * AATinvc = c using the Cholsky factorization method
       gsl_error_handler_t* old_handler =  gsl_set_error_handler_off();
       int cholesky_result = gsl_linalg_cholesky_decomp(&AAT.matrix);
       gsl_set_error_handler(old_handler);
       if (cholesky_result) {
         B2INFO("NewFitterGSL::calc2ndOrderCorr: resorting to SVD");
         // AAT is not positive definite, i.e. A does not have full column rank
         // => use the SVD of AT to solve A lambdanew = gradf
  
         // ATcopy = AT
         gsl_matrix_memcpy(&ATcopy.matrix, &AT.matrix);
  
         // SVD decomposition of Acopy
         gsl_linalg_SV_decomp_jacobi(&ATcopy.matrix, &V.matrix, &s.vector);
         // set small values to zero
         double mins = eps * std::fabs(gsl_vector_get(&s.vector, 0));
         for (int i = 0; i < ncon; ++i) {
           if (std::fabs(gsl_vector_get(&s.vector, i)) <= mins)
             gsl_vector_set(&s.vector, i, 0);
         }
         gsl_linalg_SV_solve(&ATcopy.matrix, &V.matrix, &s.vector, &c.vector, &AATinvc.vector);
       } else {
         gsl_linalg_cholesky_solve(&AAT.matrix, &c.vector, &AATinvc.vector);
       }
  
       // phat = -1*A^T*AATinvc+ 0*phat
       gsl_blas_dgemv(CblasNoTrans, -1, &AT.matrix, &AATinvc.vector, 0, &phat.vector);
       gsl_vector_set_zero(&c.vector);
  
     }

◆ calcLimitedDx()

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
	)

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

Definition at line 913 of file NewFitterGSL.cc.

◆ calcMu()

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
	)

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

Definition at line 1161 of file NewFitterGSL.cc.

◆ calcNewtonDx()

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
	)

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

Definition at line 802 of file NewFitterGSL.cc.

◆ calcpTLp()

double 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

Definition at line 1613 of file NewFitterGSL.cc.

◆ calcReducedHessian()

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`
	)

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

Definition at line 1870 of file NewFitterGSL.cc.

◆ calcReducedHessianEigenvalues()

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

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

Definition at line 1914 of file NewFitterGSL.cc.

◆ calcZ()

gsl_matrix_view 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

Definition at line 1834 of file NewFitterGSL.cc.

◆ determineLambdas()

int 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

Definition at line 1482 of file NewFitterGSL.cc.

◆ doLineSearch()

int doLineSearch	(	double &	alpha,
		gsl_vector *	vecxnew,
		int	imode,
		double	phi0,
		double	dphi0,
		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
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

Definition at line 1047 of file NewFitterGSL.cc.

◆ getGlobalCovarianceMatrix() [1/2]

double * getGlobalCovarianceMatrix ( int & idim )

virtualinherited

Parameters

idim	1st dimension of global covariance matrix

Definition at line 158 of file BaseFitter.cc.

◆ getGlobalCovarianceMatrix() [2/2]

const double * getGlobalCovarianceMatrix ( int & idim ) const

virtualinherited

Parameters

idim	1st dimension of global covariance matrix

Definition at line 148 of file BaseFitter.cc.

◆ meritFunction()

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

Parameters

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

Definition at line 1249 of file NewFitterGSL.cc.

◆ meritFunctionDeriv()

double 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

Definition at line 1284 of file NewFitterGSL.cc.

◆ MoorePenroseInverse()

void 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

Definition at line 1569 of file NewFitterGSL.cc.

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

analysis/OrcaKinFit/include/NewFitterGSL.h
analysis/OrcaKinFit/src/NewFitterGSL.cc

Public Types

Public Member Functions

Static Public Member Functions

Public Attributes

Protected Types

Protected Attributes

Detailed Description

Member Function Documentation

◆ calc2ndOrderCorr()

◆ calcLimitedDx()

◆ calcMu()

◆ calcNewtonDx()

◆ calcpTLp()

◆ calcReducedHessian()

◆ calcReducedHessianEigenvalues()

◆ calcZ()

◆ determineLambdas()

◆ doLineSearch()

◆ getGlobalCovarianceMatrix() [1/2]

◆ getGlobalCovarianceMatrix() [2/2]

◆ meritFunction()

◆ meritFunctionDeriv()

◆ MoorePenroseInverse()