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

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