source: trunk/MagicSoft/Mars/manalysis/MSimulatedAnnealing.cc@ 2811

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