MarlinKinfit  0.4.0
Public Types | Public Member Functions | Static Public Member Functions | Public Attributes | List of all members
NewFitterGSL Class Reference

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

#include <NewFitterGSL.h>

Inheritance diagram for NewFitterGSL:
BaseFitter

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 BaseTracergetTracer ()
 
virtual const BaseTracergetTracer () 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.
 
BaseFitteroperator= (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.
 
BaseTracertracer
 

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:

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
vecdxhatthe correction step
vecxnewthe current state vector (parameters must be set to vecxnew!)
MatMThe matrix of the last Newton step
MatWwork matrix
vecwwork vector
epsSingular 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
alphaOutput value alpha
muValue of mu for merit function
vecxnewNew vector x
imodemode: 0=Armijo, 1=Wolfe, 2=Goldstein
vecxCurrent vector x
vecdxhatthe 2nd order correction step, if any
vecdxUpdate vector dx
vecdxscalResult: Update vector dx, scaled
veceError vector e
MatMMatrix M
MatMscalMatrix M, scaled
MatWWork matrix
vecwWork 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
vecxCurrent vector x
veceCurrent errors x
vecdxCurrent step dx
vecdxscalCurrent step dx, scaled
xnewNew vector x
MatMCurrent matrix M
MatMscalCurrent scaled matrix M
vecwWork 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
vecdxResult: Update vector dx
vecdxscalResult: Update vector dx, scaled
vecxCurrent vector x
veceCurrent ``error'' set of x
MatMMatrix M
MatMscalMatrix M, scaled
vecyVector y
vecyscalVector y, scaled,
MatWWork matrix
MatW2Work matrix
permWWork permutation vector
vecwWork 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
vecdxCurrent step dx
MatMCurrent matrix M
vecwWork 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
rankHdimension of H
MatW1work matrix, contains H at the end
vecxvector with current x values
MatW2work matrix, contains Z at the end
MatW3work matrix
vecw1work vector
vecw2work vector
permWWork permutation vector
epsSingular 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
rankHdimension of H
MatW1work matrix, contains H at the end
vecxvector with current x values
MatW2work matrix, contains Z at the end
MatW3work matrix
vecw1work vector
vecw2work vector
permWWork permutation vector
eigenwsWork space for eigenvalue calculation
epsSingular 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
rankArank of A (= number of lin. indep. constraints)
MatW1work matrix, contains Z at the end
MatW2work matrix
vecw1work vector
vecw2work vector
permWWork permutation vector
epsSingular 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
vecxnewvector with new lambda values
MatMmatrix with constraint derivatives
vecxvector with current x values
MatWwork matrix
vecwwork vector
epsSingular 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
alphaOutput value alpha
vecxnewNew vector x
imodemode: 0=Armijo, 1=Wolfe, 2=Goldstein
phi0Merit function for alpha=0
dphi0Directional derivative of merit function for alpha=0
phiRMerit function for given alpha
etaConstant for Armijo's rule
zetaConstant for Wolfe's or Goldstein's rule
muValue of mu for merit function
vecxCurrent vector x
vecdxUpdate vector dx
veceError vector e
vecwWork 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
muValue of mu
vecxCurrent vector x
veceCurrent 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
muValue of mu
vecxCurrent vector x
veceCurrent errors x
vecdxCurrent update vector dx
vecwwork 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
AinvResult: m x n matrix A+
AInput: n x m matrix A, n >= m (is destroyed!)
WWork matrix, at least m x m
wWork vector w, at least m
epsSingular values < eps*(max(abs(s_i))) are set to 0

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