source: trunk/MagicSoft/Mars/mtools/MSimulatedAnnealing.cc@ 8070

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