source: releases/Mars.2014.05.26/manalysis/MMarquardt.cc@ 18481

Last change on this file since 18481 was 2475, checked in by wittek, 21 years ago
*** empty log message ***
File size: 9.3 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Wolfgang Wittek 10/2003 <mailto:wittek@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26// //
27// MMarquardt //
28// //
29// Marquardt method of solving nonlinear least-squares problems //
30// //
31// (see Numerical recipes (2nd ed.), W.H.Press et al., p.688 ff) //
32// //
33/////////////////////////////////////////////////////////////////////////////
34#include "MMarquardt.h"
35
36#include <math.h> // fabs
37
38#include <TVectorD.h>
39#include <TMatrixD.h>
40
41#include <TStopwatch.h>
42
43#include "MLog.h"
44#include "MLogManip.h"
45#include "MParContainer.h"
46
47ClassImp(MMarquardt);
48
49using namespace std;
50
51
52// --------------------------------------------------------------------------
53//
54// Default constructor.
55//
56MMarquardt::MMarquardt(const char *name, const char *title)
57{
58 fName = name ? name : "MMarquardt";
59 fTitle = title ? title : "Marquardt minimization";
60}
61
62// -----------------------------------------------------------------------
63//
64// Set - the number of parameters
65// - the maximum number of steps allowed in the minimization and
66// - the change in chi2 signaling convergence
67
68void MMarquardt::SetNpar(Int_t numpar, Int_t numstepmax, Double_t loopchi2)
69{
70 fNpar = numpar;
71 fNumStepMax = numstepmax;
72 fLoopChi2 = loopchi2;
73
74 fdParam.ResizeTo(fNpar);
75
76 fParam.ResizeTo(fNpar);
77 fGrad.ResizeTo(fNpar);
78 fCovar.ResizeTo(fNpar, fNpar);
79
80 fmyParam.ResizeTo(fNpar);
81 fmyGrad.ResizeTo(fNpar);
82 fmyCovar.ResizeTo(fNpar, fNpar);
83
84 fIxc.ResizeTo(fNpar);
85 fIxr.ResizeTo(fNpar);
86 fIp.ResizeTo(fNpar);
87}
88
89
90// -----------------------------------------------------------------------
91//
92// do the minimization
93//
94// fcn is the function which calculates for a given set of parameters
95// - the function L to be minimized
96// - beta_k = -1/2 * dL/da_k (a kind of gradient of L)
97// - alfa_kl = 1/2 * dL/(da_k da_l) (a kind of 2nd derivative of L)
98//
99// Vinit contains the starting values of the parameters
100//
101
102Int_t MMarquardt::Loop(
103 Bool_t (*fcn)(TVectorD &, TMatrixD &, TVectorD &, Double_t &),
104 TVectorD &Vinit)
105{
106 fFunc = fcn;
107
108 // set the initial parameter values
109 for (Int_t i=0; i<fNpar; i++)
110 fParam(i) = Vinit(i);
111
112 //-------------------------------------------
113 // first call of the function func
114 Bool_t rcfirst = FirstCall();
115 if (!rcfirst)
116 {
117 *fLog << "MMarquardt::Loop; first call of function failed " << endl;
118 return -1;
119 }
120
121 Double_t oldChi2 = fChi2;
122 Double_t fdChi2 = 1.e10;
123 Int_t fNumStep = 0;
124
125 //-------------------------------------------
126 // do the next step in the minimization
127 Bool_t rcnext;
128 do
129 {
130 fNumStep++;
131
132 rcnext = NextStep();
133 if (!rcnext) break;
134
135 fdChi2 = fabs(oldChi2-fChi2);
136 oldChi2 = fChi2;
137 } while (fdChi2 > fLoopChi2 && fNumStep < fNumStepMax);
138
139 //-------------------------------------------
140 // do the final calculations
141 if (!rcnext)
142 {
143 *fLog << "MMarquardt::Loop; function call failed in step " << fNumStep
144 << endl;
145 return -2;
146 }
147
148 if (fdChi2 > fLoopChi2)
149 {
150 *fLog << "MMarquardt::Loop; minimization has not converged, fChi2, fdChi2 = "
151 << fChi2 << ", " << fdChi2 << endl;
152 return -3;
153 }
154
155 *fLog << "MMarquardt::Loop; minimization has converged, fChi2, fdChi2, fNumStep = "
156 << fChi2 << ", " << fdChi2 << ", " << fNumStep << endl;
157
158
159 Bool_t rccov = CovMatrix();
160 if (!rccov)
161 {
162 *fLog << "MMarquardt::Loop; calculation of covariance matrix failed "
163 << endl;
164 return 1;
165 }
166
167 //-------------------------------------------
168 // results
169
170 *fLog << "MMarquardt::Loop; Results of fit : fChi2, fNumStep, fdChi2 ="
171 << fChi2 << ", " << fNumStep << ", " << fdChi2 << endl;
172
173 for (Int_t i=0; i<fNpar; i++)
174 fdParam(i) = sqrt(fCovar(i,i));
175
176 *fLog << "MMarquardt::Loop; i, Param(i), dParam(i)" << endl;
177 for (Int_t i=0; i<fNpar; i++)
178 {
179 *fLog << i << " " << fParam(i) << " " << fdParam(i) << endl;
180 }
181
182 *fLog << "MMarquardt::Loop; Covariance matrix" << endl;
183 for (Int_t i=0; i<fNpar; i++)
184 {
185 *fLog << i;
186 for (Int_t j=0; j<fNpar; j++)
187 {
188 *fLog << fCovar(i,j) << " ";
189 }
190 *fLog << endl;
191 }
192
193 return 0;
194}
195
196
197// -----------------------------------------------------------------------
198//
199// do 1st step of the minimization
200//
201
202Bool_t MMarquardt::FirstCall()
203{
204 fLambda = 0.001;
205 Bool_t rc = (*fFunc)(fParam, fCovar, fGrad, fChi2);
206 if (!rc) return kFALSE;
207
208 fCHIq = fChi2;
209 for (Int_t j=0; j<fNpar; j++)
210 fmyParam(j) = fParam(j);
211
212 return kTRUE;
213}
214
215
216// -----------------------------------------------------------------------
217//
218// do one step of the minimization
219//
220
221Bool_t MMarquardt::NextStep()
222{
223 for (Int_t j=0; j<fNpar; j++)
224 {
225 for (Int_t k=0; k<fNpar; k++)
226 fmyCovar(j,k) = fCovar(j,k);
227
228 fmyCovar(j,j) *= (1.0 + fLambda);
229 fmyGrad(j) = fGrad(j);
230 }
231
232 Bool_t rgj = GaussJordan(fNpar, fmyCovar, fmyGrad);
233 if(!rgj) return kFALSE;
234
235 for (Int_t j=0; j<fNpar; j++)
236 fmyParam(j) = fParam(j) + fmyGrad(j);
237
238 Bool_t rc = (*fFunc)(fmyParam, fmyCovar, fmyGrad, fChi2);
239 if(!rc) return kFALSE;
240
241 if (fChi2 < fCHIq)
242 {
243 fLambda *= 0.1;
244 fCHIq = fChi2;
245
246 for (Int_t j=0; j<fNpar; j++)
247 {
248 for (Int_t k=0; k<fNpar; k++)
249 fCovar(j,k) = fmyCovar(j,k);
250
251 fGrad(j) = fmyGrad(j);
252 fParam(j) = fmyParam(j);
253 }
254 }
255 else
256 fLambda *= 10.0;
257
258
259 return kTRUE;
260}
261
262// -----------------------------------------------------------------------
263//
264// calculate error matrix of fitted parameters
265//
266
267Bool_t MMarquardt::CovMatrix()
268{
269 Bool_t rc = (*fFunc)(fParam, fCovar, fGrad, fChi2);
270 if(!rc) return kFALSE;
271
272 for (Int_t j=0; j<fNpar; j++)
273 {
274 for (Int_t k=0; k<fNpar; k++)
275 fmyCovar(j,k) = fCovar(j,k);
276
277 fmyCovar(j,j) *= (1.0 + fLambda);
278 fmyGrad(j) = fGrad(j);
279 }
280
281 Bool_t rgj = GaussJordan(fNpar, fmyCovar, fmyGrad);
282 if(!rgj) return kFALSE;
283
284 for (Int_t j=0; j<fNpar; j++)
285 {
286 for (Int_t k=0; k<fNpar; k++)
287 fCovar(j,k) = fmyCovar(j,k);
288 }
289
290 return kTRUE;
291}
292
293// -----------------------------------------------------------------------
294//
295// solve normal equations
296//
297// sum(covar_kl * x_l) = beta_k (k=0,... (n-1))
298//
299// by the Gauss-Jordan method
300// (see Numerical recipes (2nd ed.), W.H.Press et al., p.39)
301//
302// on return : covar contains the inverse of the input matrix covar
303// beta contains the result for x
304//
305// return value = kTRUE means OK
306// kFALSE means singular matrix
307//
308
309Bool_t MMarquardt::GaussJordan(Int_t &n, TMatrixD &covar, TVectorD &beta)
310{
311 Int_t i, j, k, l, ll;
312 Int_t ic = 0;
313 Int_t ir = 0;
314 Double_t h, d, p;
315
316 for (j=0; j<n; j++)
317 fIp(j) = 0;
318
319 for (i=0; i<n; i++)
320 {
321 h = 0.0;
322 for (j=0; j<n; j++)
323 {
324 if (fIp(j) != 1)
325 {
326 for (k=0; k<n; k++)
327 {
328 if (fIp(k) == 0)
329 {
330 if (fabs(covar(j,k)) >= h)
331 {
332 h = fabs(covar(j,k));
333 ir = j;
334 ic = k;
335 }
336 }
337 else
338 {
339 if (fIp(k) > 1) return kFALSE;
340 }
341 }
342 }
343 }
344
345 fIp(ic)++;
346 if (ir != ic)
347 {
348 for (l=0; l<n; l++)
349 {
350 d = covar(ir,l);
351 covar(ir,l) = covar(ic,l);
352 covar(ic,l) = d;
353 }
354 d = beta(ir);
355 beta(ir) = beta(ic);
356 beta(ic) = d;
357 }
358
359 fIxr(i) = ir;
360 fIxc(i) = ic;
361 if (covar(ic,ic) == 0.0) return kFALSE;
362 p = 1.0 / covar(ic,ic);
363 covar(ic,ic) = 1.0;
364
365 for (l=0; l<n; l++)
366 covar(ic,l) *= p;
367
368 beta(ic) *= p;
369
370 for (ll=0; ll<n; ll++)
371 {
372 if (ll!= ic)
373 {
374 d = covar(ll,ic);
375 covar(ll,ic) = 0.0;
376
377 for (l=0; l<n; l++)
378 covar(ll,l) -= covar(ic,l) * d;
379
380 beta(ll) -= beta(ic) * d;
381 }
382 }
383 }
384
385 for (l=n-1; l>=0; l--)
386 {
387 if (fIxr(l) != fIxc(l))
388 {
389 for (k=0; k<n; k++)
390 {
391 d = covar(k,fIxr(l));
392 covar(k,fIxr(l)) = covar(k,fIxc(l));
393 covar(k,fIxc(l)) = d;
394 }
395 }
396 }
397
398 return kTRUE;
399}
400//=========================================================================
401
402
403
404
405
406
407
408
409
410
411
412
413
414
Note: See TracBrowser for help on using the repository browser.