source: branches/Corsika7500Compatibility/mtools/MSimulatedAnnealing.cc

Last change on this file was 8907, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 21.7 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): Markus Gaug 10/2002 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2002
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26// //
27// MSimulatedAnnealing //
28// //
29// class to perform a Simulated Annealing minimization on an n-dimensional //
30// simplex of a function 'FunctionToMinimize(TArrayF &)' in multi- //
31// dimensional parameter space. //
32// (The code is adapted from Numerical Recipies in C++, 2nd ed., //
33// pp. 457-459) //
34// //
35// Classes can inherit from MSimulatedAnnealing //
36// and use the function: //
37// //
38// RunSimulatedAnnealing(); //
39// //
40// They HAVE TO initialize the following input arguments //
41// (with ndim being the parameter dimension (max. 20)): //
42// //
43// 1) a TMatrix p(ndim+1,ndim) //
44// holding the start simplex //
45// 2) a TArrayF y(ndim+1) //
46// whose components must be pre-initialized to the values of //
47// FunctionToMinimize evaluated at the fNdim+1 vertices (rows) of p //
48// 3) a TArrayF p0(ndim) //
49// whose components contain the lower simplex borders //
50// 4) a TArrayF p1(ndim) //
51// whose components contain the upper simplex borders //
52// (The simplex will not get reflected out of these borders !!!) //
53// //
54// These arrays have to be initialized with a call to: //
55// Initialize(TMatrix \&, TArrayF \&, TArrayF \&, TArrayF \&) //
56// //
57// 5) a virtual function FunctionToMinimize(TArrayF &) //
58// acting on a TArrayF(ndim) array of parameter values //
59// //
60// Additionally, a global start temperature can be chosen with: //
61// //
62// SetStartTemperature(Float_t temp) //
63// (default is: 10) //
64// //
65// A total number of total moves (watch out for the CPU time!!!) with: //
66// //
67// SetNumberOfMoves(Float_t totalMoves) //
68// (default is: 200) //
69// //
70// The temperature is reduced after evaluation step like: //
71// CurrentTemperature = StartTemperature*(1-currentMove/totalMoves) //
72// where currentMove is the cumulative number of moves so far //
73// //
74// WARNING: The start temperature and number of moves has to be optimized //
75// for each individual problem. //
76// It is not straightforward using the defaults! //
77// In case, you omit this important step, //
78// you will get local minima without even noticing it!! //
79// //
80// You may define the following variables: //
81// //
82// 1) A global convergence criterium fTol //
83// which determines an early return for: //
84// //
85// max(FunctionToMinimize(p))-min(FunctionToMinimize(p)) //
86// ----------------------------------------------------- \< fTol //
87// max(FunctionToMinimize(p))+min(FunctionToMinimize(p)) //
88// //
89// ModifyTolerance(Float_t) //
90// //
91// 2) A verbose level for prints to *fLog //
92// //
93// SetVerbosityLevel(Verbosity_t) //
94// //
95// 3) A bit if you want to have stored //
96// the full simplex after every call to Amebsa: //
97// //
98// SetFullStorage() //
99// //
100// 4) The random number generator //
101// e.g. if you want to test the stability of the output //
102// //
103// SetRandom(TRandom *rand) //
104// //
105// //
106// Output containers: //
107// //
108// MHSimulatedAnnealing //
109// //
110// Use: //
111// GetResult()->Draw(Option_t *o) //
112// or //
113// GetResult()->DrawClone(Option_t *o) //
114// //
115// to retrieve the output histograms //
116// //
117//////////////////////////////////////////////////////////////////////////////
118#include "MSimulatedAnnealing.h"
119
120#include <fstream>
121#include <iostream>
122
123#include <TMath.h>
124#include <TRandom.h>
125
126#include "MLog.h"
127#include "MLogManip.h"
128
129#include "MHSimulatedAnnealing.h"
130
131const Float_t MSimulatedAnnealing::gsYtryStr = 10000000;
132const Float_t MSimulatedAnnealing::gsYtryCon = 20000000;
133const Int_t MSimulatedAnnealing::gsMaxDim = 20;
134const Int_t MSimulatedAnnealing::gsMaxStep = 50;
135
136ClassImp(MSimulatedAnnealing);
137
138using namespace std;
139
140// ---------------------------------------------------------------------------
141//
142// Default Constructor
143// Initializes random number generator and default variables
144//
145MSimulatedAnnealing::MSimulatedAnnealing()
146 : fResult(NULL), fTolerance(0.0001),
147 fNdim(0), fNumberOfMoves(200),
148 fStartTemperature(10), fFullStorage(kFALSE),
149 fInit(kFALSE),
150 fP(gsMaxDim, gsMaxDim), fP0(gsMaxDim),
151 fP1(gsMaxDim), fY(gsMaxDim), fYb(gsMaxDim), fYconv(gsMaxDim),
152 fPb(gsMaxDim), fPconv(gsMaxDim),
153 fBorder(kEStrictBorder), fVerbose(kEDefault)
154{
155
156 // random number generator
157 fRandom = gRandom;
158}
159
160// --------------------------------------------------------------------------
161//
162// Destructor.
163//
164MSimulatedAnnealing::~MSimulatedAnnealing()
165{
166 if (fResult)
167 delete fResult;
168}
169
170// ---------------------------------------------------------------------------
171//
172// Initialization needs the following four members:
173//
174// 1) a TMatrix p(ndim+1,ndim)
175// holding the start simplex
176// 2) a TVector y(ndim+1)
177// whose components must be pre-initialized to the values of
178// FunctionToMinimize evaluated at the fNdim+1 vertices (rows) of fP
179// 3) a TVector p0(ndim)
180// whose components contain the lower simplex borders
181// 4) a TVector p1(ndim)
182// whose components contain the upper simplex borders
183// (The simplex will not get reflected out of these borders !!!)
184//
185// It is possible to perform an initialization and
186// a subsequent RunMinimization several times.
187// Each time, a new class MHSimulatedAnnealing will get created
188// (and destroyed).
189//
190Bool_t MSimulatedAnnealing::Initialize(const TMatrix &p, const TVector &y,
191 const TVector &p0, const TVector &p1)
192{
193
194 fNdim = p.GetNcols();
195 fMpts = p.GetNrows();
196
197 //
198 // many necessary checks ...
199 //
200 if (fMpts > gsMaxDim)
201 {
202 gLog << err << "Dimension of Matrix fP is too big ... aborting." << endl;
203 return kFALSE;
204 }
205 if (fNdim+1 != fMpts)
206 {
207 gLog << err << "Matrix fP does not have the right dimensions ... aborting." << endl;
208 return kFALSE;
209 }
210 if (y.GetNrows() != fMpts)
211 {
212 gLog << err << "Array fY has not the right dimension ... aborting." << endl;
213 return kFALSE;
214 }
215 if (p0.GetNrows() != fNdim)
216 {
217 gLog << err << "Array fP0 has not the right dimension ... aborting." << endl;
218 return kFALSE;
219 }
220 if (p1.GetNrows() != fNdim)
221 {
222 gLog << err << "Array fP1 has not the right dimension ... aborting." << endl;
223 return kFALSE;
224 }
225
226 //
227 // In order to allow multiple use of the class
228 // without need to construct the class every time new
229 // delete the old fResult and create a new one in RunMinimization
230 //
231 if (fResult)
232 delete fResult;
233
234 fY.ResizeTo(fMpts);
235
236 fPsum.ResizeTo(fNdim);
237 fPconv.ResizeTo(fNdim);
238
239 fP0.ResizeTo(fNdim);
240 fP1.ResizeTo(fNdim);
241 fPb.ResizeTo(fNdim);
242
243 fP.ResizeTo(fMpts,fNdim);
244
245 fY = y;
246 fP = p;
247 fP0 = p0;
248 fP1 = p1;
249 fPconv.Zero();
250
251 fInit = kTRUE;
252 fYconv = 0.;
253
254 return kTRUE;
255}
256
257
258// ---------------------------------------------------------------------------
259//
260// RunMinimization:
261//
262// Runs only eafter a call to Initialize(const TMatrix \&, const TVector \&,
263// const TVector \&, const TVector \&)
264//
265// Temperature and number of moves should have been set
266// (default: StartTemperature = 10, NumberOfMoves = 200
267//
268//
269// It is possible to perform an initialization and
270// a subsequent RunMinimization several times.
271// Each time, a new class MHSimulatedAnnealing will get created
272// (and destroyed).
273Bool_t MSimulatedAnnealing::RunMinimization()
274{
275 if (!fInit)
276 {
277 gLog << err << "No succesful initialization performed yet... aborting." << endl;
278 return kFALSE;
279 }
280
281 Int_t iter = 0;
282 UShort_t iret = 0;
283 UShort_t currentMove = 0;
284 Real_t currentTemp = fStartTemperature;
285
286 fResult = new MHSimulatedAnnealing(fNumberOfMoves,fNdim);
287 if (fFullStorage)
288 fResult->InitFullSimplex();
289
290 while(1)
291 {
292 if (iter > 0)
293 {
294 gLog << "Convergence at move: " << currentMove ;
295 gLog << " and temperature: " << currentTemp << endl;
296 break;
297 }
298
299 if (currentTemp > 0.)
300 {
301 //
302 // Reduce the temperature
303 //
304 // FIXME: Maybe it is necessary to also incorporate other
305 // ways to reduce the temperature (T0*(1-k/K)**alpha)
306 //
307 currentTemp = fStartTemperature*(1.-(float)currentMove++/fNumberOfMoves);
308 iter = 1;
309 }
310 else
311 {
312 // Make sure that now, the program will return only on convergence !
313 // The program returns to here only after gsMaxStep moves
314 // If we have not reached convergence until then, we assume that an infinite
315 // loop has occurred and quit.
316 if (iret != 0)
317 {
318 gLog << warn << "No Convergence at the end ! " << endl;
319 fY.Zero();
320
321 break;
322 }
323 iter = 150;
324 iret++;
325 currentMove++;
326 }
327
328 if (fVerbose==2) {
329 gLog << dbginf << " current..." << endl;
330 gLog << " - move: " << currentMove << endl;
331 gLog << " - temperature: " << currentTemp << endl;
332 gLog << " - best function evaluation: " << fYb << endl;
333 }
334
335 iter = Amebsa(iter, currentTemp);
336
337 // Store the current best values in the histograms
338 fResult->StoreBestValueEver(fPb,fYb,currentMove);
339
340 // Store the complete simplex if we have full storage
341 if (fFullStorage)
342 fResult->StoreFullSimplex(fP,currentMove);
343 }
344
345 //
346 // Now, the matrizes and vectors have all the same value,
347 // Need to initialize again to allow a new Minimization
348 //
349 fInit = kFALSE;
350
351 return kTRUE;
352}
353
354// ---------------------------------------------------------------------------
355//
356// Amebsa
357//
358// This is the (adjusted) amebsa function from
359// Numerical Recipies (pp. 457-458)
360//
361// The routine makes iter function evaluations at an annealing
362// temperature fCurrentTemp, then returns. If iter is returned
363// with a poisitive value, then early convergence has occurred.
364//
365Int_t MSimulatedAnnealing::Amebsa(Int_t iter, const Real_t temp)
366{
367 GetPsum();
368
369 while (1)
370 {
371 UShort_t ihi = 0; // Simplex point with highest function evaluation
372 UShort_t ilo = 1; // Simplex point with lowest function evaluation
373
374 // Function eval. at ilo (with random fluctuations)
375 Real_t ylo = fY(0) + gRandom->Exp(temp);
376
377 // Function eval. at ihi (with random fluctuations)
378 Real_t yhi = fY(1) + gRandom->Exp(temp);
379
380 // The function evaluation at next highest point
381 Real_t ynhi = ylo;
382
383 if (ylo > yhi)
384 {
385 // Determine which point is the highest (worst),
386 // next-highest and lowest (best)
387 ynhi = yhi;
388 yhi = ylo;
389 ylo = ynhi;
390 }
391
392 // By looping over the points in the simplex
393 for (UShort_t i=2;i<fMpts;i++)
394 {
395 const Real_t yt = fY(i) + gRandom->Exp(temp);
396
397 if (yt <= ylo)
398 {
399 ilo = i;
400 ylo = yt;
401 }
402
403 if (yt > yhi)
404 {
405 ynhi = yhi;
406 ihi = i;
407 yhi = yt;
408 }
409 else
410 if (yt > ynhi)
411 ynhi = yt;
412 }
413
414 // Now, fY(ilo) is smallest and fY(ihi) is at biggest value
415 if (iter < 0)
416 {
417 // Enough looping with this temperature, go to decrease it
418 // First put best point and value in slot 0
419
420 Real_t dum = fY(0);
421 fY(0) = fY(ilo);
422 fY(ilo) = dum;
423
424 for (UShort_t n=0;n<fNdim;n++)
425 {
426 dum = fP(0,n);
427 fP(0,n) = fP(ilo,n);
428 fP(ilo,n) = dum;
429 }
430
431 break;
432 }
433
434 // Compute the fractional range from highest to lowest and
435 // return if satisfactory
436 Real_t tol = fabs(yhi) + fabs(ylo);
437 if (tol != 0)
438 tol = 2.0*fabs(yhi-ylo)/tol;
439
440 if (tol<fTolerance)
441 {
442 // Put best point and value in fPconv
443 fYconv = fY(ilo);
444 for (UShort_t n=0; n<fNdim; n++)
445 fPconv(n) = fP(ilo, n);
446
447 break;
448 }
449 iter -= 2;
450
451 // Begin new Iteration. First extrapolate by a factor of -1 through
452 // the face of the simplex across from the high point, i.e. reflect
453 // the simplex from the high point
454 Real_t ytry = Amotsa(-1.0, ihi, yhi,temp);
455
456 if (ytry <= ylo)
457 {
458 // cout << " !!!!!!!!!!!!!! E X P A N D !!!!!!!!!!!!!!" << endl;
459 // Gives a result better than the best point, so try an additional
460 // extrapolation by a factor of 2
461 ytry = Amotsa(2.0, ihi, yhi,temp);
462 continue;
463 }
464
465 if (ytry < ynhi)
466 {
467 iter++;
468 continue;
469 }
470
471 // cout << " !!!!!!!!!!!! R E F L E C T !!!!!!!!!!!!!!!!!!!!" << endl;
472 // The reflected point is worse than the second-highest, so look for an
473 // intermediate lower point, for (a one-dimensional contraction */
474 const Real_t fYsave = yhi;
475 ytry = Amotsa(0.5, ihi, yhi,temp);
476
477 if (ytry < fYsave)
478 continue;
479
480 // cout << " !!!!!!!!!!!! R E F L E C T !!!!!!!!!!!!!!!!!!!!" << endl;
481 // The reflected point is worse than the second-highest, so look for an
482 // intermediate lower point, for (a one-dimensional contraction */
483 const Real_t ysave = yhi;
484 ytry = Amotsa(0.5, ihi, yhi,temp);
485
486 if (ytry < ysave)
487 continue;
488
489 // cout << " !!!!!!!!!!!! C O N T R A C T !!!!!!!!!!!!!!!!!!" << endl;
490 // Cannot seem to get rid of that point, better contract around the
491 // lowest (best) point
492 for (UShort_t i=0; i<fMpts; i++)
493 {
494 if (i != ilo)
495 {
496 for (UShort_t j=0;j<fNdim;j++)
497 {
498 fPsum(j) = 0.5*(fP(i, j) + fP(ilo, j));
499
500 // Create new cutvalues
501 fP(i, j) = fPsum(j);
502 }
503 fY(i) = FunctionToMinimize(fPsum);
504 }
505 }
506
507 iter -= fNdim;
508 GetPsum();
509 }
510 return iter;
511}
512
513void MSimulatedAnnealing::GetPsum()
514{
515 for (Int_t n=0; n<fNdim; n++)
516 {
517 Real_t sum=0.0;
518 for (Int_t m=0;m<fMpts;m++)
519 sum += fP(m,n);
520
521 fPsum(n) = sum;
522 }
523}
524
525
526Real_t MSimulatedAnnealing::Amotsa(const Float_t fac, const UShort_t ihi,
527 Real_t &yhi, const Real_t temp)
528{
529
530 const Real_t fac1 = (1.-fac)/fNdim;
531 const Real_t fac2 = fac1 - fac;
532
533 Int_t borderflag = 0;
534 TVector ptry(fNdim);
535 TVector cols(fMpts);
536
537 for (Int_t j=0; j<fNdim; j++)
538 {
539 ptry(j) = fPsum(j)*fac1 - fP(ihi, j)*fac2;
540
541 // Check that the simplex does not go to infinite values,
542 // in case of: reflect it
543 const Real_t newcut = ptry(j);
544
545 if (fP1(j) > fP0(j))
546 {
547 if (newcut > fP1(j))
548 {
549 ptry(j) = fP1(j);
550 borderflag = 1;
551 }
552 else
553 if (newcut < fP0(j))
554 {
555 ptry(j) = fP0(j);
556 borderflag = 1;
557 }
558 }
559
560 else
561 {
562 if (newcut < fP1(j))
563 {
564 ptry(j) = fP1(j);
565 borderflag = 1;
566 }
567 else
568 if (newcut > fP0(j))
569 {
570 ptry(j) = fP0(j);
571 borderflag = 1;
572 }
573 }
574 }
575
576 Real_t faccompare = 0.5;
577 Real_t ytry = 0;
578
579 switch (borderflag)
580 {
581 case kENoBorder:
582 ytry = FunctionToMinimize(fPsum);
583 break;
584
585 case kEStrictBorder:
586 ytry = FunctionToMinimize(fPsum) + gsYtryStr;
587 break;
588
589 case kEContractBorder:
590 ytry = fac == faccompare ? gsYtryCon : gsYtryStr;
591 break;
592 }
593
594 if (ytry < fYb)
595 {
596 fPb = ptry;
597 fYb = ytry;
598 }
599
600 const Real_t yflu = ytry + gRandom->Exp(temp);
601
602 if (yflu >= yhi)
603 return yflu;
604
605 fY(ihi) = ytry;
606 yhi = yflu;
607
608 for(Int_t j=0; j<fNdim; j++)
609 {
610 fPsum(j) += ptry(j)-fP(ihi, j);
611 fP(ihi, j) = ptry(j);
612 }
613
614 return yflu;
615}
616
617// ---------------------------------------------------------------------------
618//
619// Dummy FunctionToMinimize
620//
621// A class inheriting from MSimulatedAnnealing NEEDS to contain a similiar
622//
623// virtual Float_t FunctionToMinimize(const TVector \&)
624//
625// The TVector contains the n parameters (dimensions) of the function
626//
627Float_t MSimulatedAnnealing::FunctionToMinimize(const TVector &arr)
628{
629 return 0.0;
630}
Note: See TracBrowser for help on using the repository browser.