MarlinKinfit  0.4.0
NewFitterGSL.h
Go to the documentation of this file.
1 
27 #ifndef __NEWFITTERGSL_H
28 #define __NEWFITTERGSL_H
29 
30 #include "BaseFitter.h"
31 
32 #include <gsl/gsl_vector.h>
33 #include <gsl/gsl_matrix.h>
34 #include <gsl/gsl_permutation.h>
35 #include <gsl/gsl_eigen.h>
36 
37 // Class NewFitterGSL
39 
53 class NewFitterGSL : public BaseFitter {
54  public:
56  NewFitterGSL();
58  virtual ~NewFitterGSL();
59 
61  virtual double fit();
62 
64  virtual int getError() const;
65 
67  virtual double getProbability() const;
69  virtual double getChi2() const;
71  virtual int getDoF() const;
73  virtual int getIterations() const;
74 
76  virtual int getNcon() const;
77 
79  virtual int getNsoft() const;
80 
82  virtual int getNpar() const;
83 
85  virtual int getNunm() const;
86 
88  virtual bool initialize();
89 
91  virtual void setDebug (int debuglevel);
92 
94  virtual void determineLambdas (gsl_vector *vecxnew,
95  const gsl_matrix *MatM,
96  const gsl_vector *vecx,
97  gsl_matrix *MatW,
98  gsl_vector *vecw,
99  double eps = 0
100  );
101 
103  virtual void calc2ndOrderCorr ( gsl_vector *vecdxhat,
104  const gsl_vector *vecxnew,
105  const gsl_matrix *MatM,
106  gsl_matrix *MatW,
107  gsl_vector *vecw,
108  double eps = 0
109  );
110 
112  virtual gsl_matrix_view calcZ (int& rankA,
113  gsl_matrix *MatW1,
114  gsl_matrix *MatW2,
115  gsl_vector *vecw1,
116  gsl_vector *vecw2,
117  gsl_permutation *permW,
118  double eps = 0
119  );
120 
122  virtual gsl_matrix_view calcReducedHessian ( int& rankH,
123  gsl_matrix *MatW1,
124  const gsl_vector *vecx,
125  gsl_matrix *MatW2,
126  gsl_matrix *MatW3,
127  gsl_vector *vecw1,
128  gsl_vector *vecw2,
129  gsl_permutation *permW,
130  double eps = 0
131  );
132 
133  virtual gsl_vector_view calcReducedHessianEigenvalues ( int& rankH,
134  gsl_matrix *MatW1,
135  const gsl_vector *vecx,
136  gsl_matrix *MatW2,
137  gsl_matrix *MatW3,
138  gsl_vector *vecw1,
139  gsl_vector *vecw2,
140  gsl_permutation *permW,
141  gsl_eigen_symm_workspace *eigenws,
142  double eps = 0
143  );
144  public:
146  virtual double calcChi2();
147 
148  void fillx (gsl_vector *vecx);
149  void fillperr(gsl_vector *vece);
150 
151  // Fill matrix MatM, using lambdas from vecx
152  void assembleM (gsl_matrix *MatM, const gsl_vector *vecx, bool errorpropagation = false);
153 
154  // Fill matrix MatM with 2nd derivative of Lagrangian, using lambdas from vecx
155  void assembleG (gsl_matrix *MatM, const gsl_vector *vecx);
156 
157  // Scale matrix MatM to MatMscal using errors from vece
158  void scaleM (gsl_matrix *MatMscal, const gsl_matrix *MatM, const gsl_vector *vece);
159 
160  // Fill vector y, using lambdas from vecx
161  void assembley(gsl_vector *vecy, const gsl_vector *vecx);
162 
163  // Scale vector vecy to vecyscal using errors from vece
164  void scaley (gsl_vector *vecyscal, const gsl_vector *vecy, const gsl_vector *vece);
165 
166  // Fill chi2 derivatives into vector y
167  void assembleChi2Der(gsl_vector *vecy);
168 
169  // Fill vector y with values of constraints
170  void addConstraints(gsl_vector *vecy);
171 
172  // Fill constraint derivatives into Matrix M
173  void assembleConstDer(gsl_matrix *MatM);
174 
175  // Calculate Newton step update vector vecdx from current point vecx and errors vece
176  int calcNewtonDx ( gsl_vector *vecdx,
177  gsl_vector *vecdxscal,
178  gsl_vector *vecx,
179  const gsl_vector *vece,
180  gsl_matrix *MatM,
181  gsl_matrix *MatMscal,
182  gsl_vector *vecy,
183  gsl_vector *vecyscal,
184  gsl_matrix *MatW,
185  gsl_matrix *MatW2,
186  gsl_permutation *permW,
187  gsl_vector *vecw
188  );
189 
190  // Calculate limited step after linesearch
191  int calcLimitedDx ( double& alpha,
192  double& mu,
193  gsl_vector *vecxnew,
194  int imode,
195  gsl_vector *vecx,
196  gsl_vector *vecdxhat,
197  const gsl_vector *vecdx,
198  const gsl_vector *vecdxscal,
199  const gsl_vector *vece,
200  const gsl_matrix *MatM,
201  const gsl_matrix *MatMscal,
202  gsl_matrix *MatW,
203  gsl_vector *vecw
204  );
205 
206  // Perform a line search
207  int doLineSearch ( double& alpha,
208  gsl_vector *vecxnew,
209  int imode,
210  double phi0,
211  double dphi0,
212  double phiR,
213  double eta,
214  double zeta,
215  double mu,
216  const gsl_vector *vecx,
217  const gsl_vector *vecdx,
218  const gsl_vector *vece,
219  gsl_vector *vecw
220  );
221 
222  // Calculate mu for the merit function
223  double calcMu (const gsl_vector *vecx,
224  const gsl_vector *vece,
225  const gsl_vector *vecdx,
226  const gsl_vector *vecdxscal,
227  const gsl_vector *xnew,
228  const gsl_matrix *MatM,
229  const gsl_matrix *MatMscal,
230  gsl_vector *vecw
231  );
232 
233  // Calculate the merit function at x
234  double meritFunction ( double mu,
235  const gsl_vector *vecx,
236  const gsl_vector *vece
237  );
238  // Calculate the directional derivative of the merit function at x
239  double meritFunctionDeriv ( double mu,
240  const gsl_vector *vecx,
241  const gsl_vector *vece,
242  const gsl_vector *vecdx,
243  gsl_vector *vecw
244  );
245 
246  // Transfer values from vecx to FitObjects
247  bool updateParams (gsl_vector *vecx);
248 
249 
250  int invertM();
251 
252  void calcCovMatrix(gsl_matrix *MatW, gsl_permutation *permW, gsl_vector *vecx);
253 
254  enum {NPARMAX=50, NCONMAX=10, NUNMMAX=10};
255 
256  int npar;
257  int ncon;
258  int nsoft;
259  int nunm;
260  int ierr;
261  int nit;
262 
263  double fitprob;
264  double chi2;
265 
266  static void ini_gsl_permutation (gsl_permutation *&p, unsigned int size);
267  static void ini_gsl_vector (gsl_vector *&v, int unsigned size);
268  static void ini_gsl_matrix (gsl_matrix *&m, int unsigned size1, unsigned int size2);
269 
270  static void debug_print (const gsl_matrix *m, const char *name);
271  static void debug_print (const gsl_vector *v, const char *name);
272 
273  // z = x + a y
274  static void add (gsl_vector *vecz, const gsl_vector *vecx, double a, const gsl_vector *vecy);
275 
276  // Check whether all elements are finite
277  static bool isfinite (const gsl_vector *vec);
278 
279  // Check whether all elements are finite
280  static bool isfinite (const gsl_matrix *mat);
281 
283  static void MoorePenroseInverse (gsl_matrix *Ainv,
284  gsl_matrix *A,
285  gsl_matrix *W,
286  gsl_vector *w,
287  double eps = 0
288  );
289  double calcpTLp (const gsl_vector *vecdx,
290  const gsl_matrix *MatM,
291  gsl_vector *vecw
292  );
293 
295  int solveSystem ( gsl_vector *vecdxscal,
296  double& detW,
297  const gsl_vector *vecyscal,
298  const gsl_matrix *MatMscal,
299  gsl_matrix *MatW,
300  gsl_matrix *MatW2,
301  gsl_vector *vecw,
302  double epsLU,
303  double epsSV
304  );
305 
307  int solveSystemLU ( gsl_vector *vecdxscal,
308  double& detW,
309  const gsl_vector *vecyscal,
310  const gsl_matrix *MatMscal,
311  gsl_matrix *MatW,
312  gsl_vector *vecw,
313  double eps
314  );
315 
317  int solveSystemSVD ( gsl_vector *vecdxscal,
318  double& detW,
319  const gsl_vector *vecyscal,
320  const gsl_matrix *MatMscal,
321  gsl_matrix *MatW,
322  gsl_matrix *MatW2,
323  gsl_vector *vecw,
324  double eps
325  );
326 
327  public:
328  unsigned int idim;
329  gsl_vector *x;
330  gsl_vector *xold;
331  gsl_vector *xnew;
332 // gsl_vector *xbest;
333  gsl_vector *dx;
334  gsl_vector *dxscal;
335 // gsl_vector *grad;
336  gsl_vector *y;
337  gsl_vector *yscal;
338  gsl_vector *perr;
339  gsl_vector *v1;
340  gsl_vector *v2;
341 // gsl_vector *Meval;
342 
343  gsl_matrix *M;
344  gsl_matrix *Mscal;
345  gsl_matrix *W;
346  gsl_matrix *W2;
347  gsl_matrix *W3;
348  // these are only used locally in calcCovMatrix
349  gsl_matrix *M1;
350  gsl_matrix *M2;
351  gsl_matrix *M3;
352  gsl_matrix *M4;
353  gsl_matrix *M5;
354 // gsl_matrix *Mevec;
355  gsl_matrix *CC;
356  gsl_matrix *CC1;
357  gsl_matrix *CCinv;
358 
359  gsl_permutation *permW;
360  gsl_eigen_symm_workspace *eigenws;
361 // gsl_eigen_symmv_workspace *eigenwsv;
362  unsigned int eigenwsdim;
363 
364  double chi2best;
365  double chi2new;
366  double chi2old;
367  double fvalbest;
368  double scale;
369  double scalebest;
370  double stepsize;
371  double stepbest;
372  enum {NITMAX = 100};
373  double scalevals[NITMAX];
374  double fvals[NITMAX];
375 
376  int imerit;
377  bool try2ndOrderCorr;
378 
379  int debug;
380 };
381 
382 #endif // __NEWFITTERGSL_H
virtual int getNcon() const
Get the number of hard constraints of the last fit.
Definition: NewFitterGSL.cc:463
int nunm
total number of unmeasured parameters
Definition: NewFitterGSL.h:259
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)
Definition: NewFitterGSL.cc:1885
double chi2
final chi2
Definition: NewFitterGSL.h:264
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.
Definition: NewFitterGSL.cc:1456
Abstract base class for fitting engines of kinematic fits.
Definition: BaseFitter.h:63
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 span...
Definition: NewFitterGSL.cc:1807
int ncon
total number of hard constraints
Definition: NewFitterGSL.h:257
virtual double fit()
The fit method, returns the fit probability.
Definition: NewFitterGSL.cc:114
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
Definition: NewFitterGSL.cc:1680
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
Definition: NewFitterGSL.cc:1714
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)
Definition: NewFitterGSL.cc:767
double meritFunction(double mu, const gsl_vector *vecx, const gsl_vector *vece)
Definition: NewFitterGSL.cc:1228
double calcpTLp(const gsl_vector *vecdx, const gsl_matrix *MatM, gsl_vector *vecw)
Definition: NewFitterGSL.cc:1590
Declares class BaseFitter.
virtual double getChi2() const
Get the chi**2 of the last fit.
Definition: NewFitterGSL.cc:398
A kinematic fitter using the Newton-Raphson method to solve the equations.
Definition: NewFitterGSL.h:53
virtual int getNpar() const
Get the number of all parameters of the last fit.
Definition: NewFitterGSL.cc:466
int nsoft
total number of soft constraints
Definition: NewFitterGSL.h:258
virtual bool initialize()
Initialize the fitter.
Definition: NewFitterGSL.cc:300
int ierr
Error status.
Definition: NewFitterGSL.h:260
double fitprob
fit probability
Definition: NewFitterGSL.h:263
virtual int getNunm() const
Get the number of unmeasured parameters of the last fit.
Definition: NewFitterGSL.cc:465
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)
Definition: NewFitterGSL.cc:880
int nit
Number of iterations.
Definition: NewFitterGSL.h:261
NewFitterGSL()
Constructor.
Definition: NewFitterGSL.cc:58
virtual int getDoF() const
Get the number of degrees of freedom of the last fit.
Definition: NewFitterGSL.cc:399
virtual ~NewFitterGSL()
Virtual destructor.
Definition: NewFitterGSL.cc:79
virtual double calcChi2()
Calculate the chi2.
Definition: NewFitterGSL.cc:381
double meritFunctionDeriv(double mu, const gsl_vector *vecx, const gsl_vector *vece, const gsl_vector *vecdx, gsl_vector *vecw)
Definition: NewFitterGSL.cc:1264
int npar
total number of parameters
Definition: NewFitterGSL.h:256
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.
Definition: NewFitterGSL.cc:1609
virtual void setDebug(int debuglevel)
Set the Debug Level.
Definition: NewFitterGSL.cc:1342
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)
Definition: NewFitterGSL.cc:1018
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.
Definition: NewFitterGSL.cc:1547
virtual int getError() const
Get the error code of the last fit: 0=OK, 1=failed.
Definition: NewFitterGSL.cc:396
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...
Definition: NewFitterGSL.cc:1842
virtual double getProbability() const
Get the fit probability of the last fit.
Definition: NewFitterGSL.cc:397
virtual int getIterations() const
Get the number of iterations of the last fit.
Definition: NewFitterGSL.cc:400
virtual int getNsoft() const
Get the number of soft constraints of the last fit.
Definition: NewFitterGSL.cc:464
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
Definition: NewFitterGSL.cc:1767
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)
Definition: NewFitterGSL.cc:1136