Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2816) +++ trunk/MagicSoft/Mars/Changelog (revision 2817) @@ -23,5 +23,7 @@ * manalysis/MSimulatedAnnealing.[h,cc] + * manalysis/Makefile * mhist/MHSimulatedAnnealing.[h,cc] + * mhist/Makefile - moved to directory mtools - MSimulatedAnnealing now inherits from TObject Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 2816) +++ trunk/MagicSoft/Mars/NEWS (revision 2817) @@ -1,4 +1,11 @@ -*-*- END -*-*- *** Version + + - added all class to perform the calibration, see macro calibration.C + + - added class MFFT to perform Fast Fourier Transforms + + - added class MSimulatedAnnealing to perform simulated annealing + minimizations - added new macro bootcampstandardanalysis.C which holds the skeleton Index: trunk/MagicSoft/Mars/manalysis/MSimulatedAnnealing.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MSimulatedAnnealing.cc (revision 2816) +++ (revision ) @@ -1,627 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Markus Gaug 10/2002 -! -! Copyright: MAGIC Software Development, 2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MSimulatedAnnealing // -// // -// class to perform a Simulated Annealing minimization on an n-dimensional // -// simplex of a function 'FunctionToMinimize(TArrayF &)' in multi- // -// dimensional parameter space. // -// (The code is adapted from Numerical Recipies in C++, 2nd ed., // -// pp. 457-459) // -// // -// Classes can inherit from MSimulatedAnnealing // -// and use the function: // -// // -// RunSimulatedAnnealing(); // -// // -// They HAVE TO initialize the following input arguments // -// (with ndim being the parameter dimension (max. 20)): // -// // -// 1) a TMatrix p(ndim+1,ndim) // -// holding the start simplex // -// 2) a TArrayF y(ndim+1) // -// whose components must be pre-initialized to the values of // -// FunctionToMinimize evaluated at the fNdim+1 vertices (rows) of fP // -// 3) a TArrayF p0(ndim) // -// whose components contain the lower simplex borders // -// 4) a TArrayF p1(ndim) // -// whose components contain the upper simplex borders // -// (The simplex will not get reflected out of these borders !!!) // -// // -// These arrays have to be initialized with a call to: // -// Initialize(TMatrix \&, TArrayF \&, TArrayF \&, TArrayF \&) // -// // -// 5) a virtual function FunctionToMinimize(TArrayF &) // -// acting on a TArrayF(ndim) array of parameter values // -// // -// Additionally, a global start temperature can be chosen with: // -// // -// SetStartTemperature(Float_t temp) // -// (default is: 10) // -// // -// A total number of total moves (watch out for the CPU time!!!) with: // -// // -// SetNumberOfMoves(Float_t totalMoves) // -// (default is: 200) // -// // -// The temperature is reduced after evaluation step like: // -// CurrentTemperature = StartTemperature*(1-currentMove/totalMoves) // -// where currentMove is the cumulative number of moves so far // -// // -// WARNING: The start temperature and number of moves has to be optimized // -// for each individual problem. // -// It is not straightforward using the defaults! // -// In case, you omit this important step, // -// you will get local minima without even noticing it!! // -// // -// You may define the following variables: // -// // -// 1) A global convergence criterium fTol // -// which determines an early return for: // -// // -// max(FunctionToMinimize(p))-min(FunctionToMinimize(p)) // -// ----------------------------------------------------- \< fTol // -// max(FunctionToMinimize(p))+min(FunctionToMinimize(p)) // -// // -// ModifyTolerance(Float_t) // -// // -// 2) A verbose level for prints to *fLog // -// // -// SetVerbosityLevel(Verbosity_t) // -// // -// 3) A bit if you want to have stored // -// the full simplex after every call to Amebsa: // -// // -// SetFullStorage() // -// // -// 4) The random number generator // -// e.g. if you want to test the stability of the output // -// // -// SetRandom(TRandom *rand) // -// // -// // -// Output containers: // -// // -// MHSimulatedAnnealing // -// // -// Use: // -// GetResult()->Draw(Option_t *o) // -// or // -// GetResult()->DrawClone(Option_t *o) // -// // -// to retrieve the output histograms // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MSimulatedAnnealing.h" - -#include - -#include "MLog.h" -#include "MLogManip.h" - -#include "MHSimulatedAnnealing.h" - -const Float_t MSimulatedAnnealing::gsYtryStr = 10000000; -const Float_t MSimulatedAnnealing::gsYtryCon = 20000000; -const Int_t MSimulatedAnnealing::gsMaxDim = 20; -const Int_t MSimulatedAnnealing::gsMaxStep = 50; - -ClassImp(MSimulatedAnnealing); - -using namespace std; - -// --------------------------------------------------------------------------- -// -// Default Constructor -// Initializes random number generator and default variables -// -MSimulatedAnnealing::MSimulatedAnnealing() - : fResult(NULL), fTolerance(0.0001), - fNdim(0), fNumberOfMoves(200), - fStartTemperature(10), fFullStorage(kFALSE), - fInit(kFALSE), - fP(gsMaxDim, gsMaxDim), fP0(gsMaxDim), - fP1(gsMaxDim), fY(gsMaxDim), fYb(gsMaxDim), fYconv(gsMaxDim), - fPb(gsMaxDim), fPconv(gsMaxDim), - fBorder(kEStrictBorder), fVerbose(kEDefault) -{ - - // random number generator - fRandom = gRandom; -} - -// -------------------------------------------------------------------------- -// -// Destructor. -// -MSimulatedAnnealing::~MSimulatedAnnealing() -{ - if (fResult) - delete fResult; -} - -// --------------------------------------------------------------------------- -// -// Initialization needs the following four members: -// -// 1) a TMatrix p(ndim+1,ndim) -// holding the start simplex -// 2) a TVector y(ndim+1) -// whose components must be pre-initialized to the values of -// FunctionToMinimize evaluated at the fNdim+1 vertices (rows) of fP -// 3) a TVector p0(ndim) -// whose components contain the lower simplex borders -// 4) a TVector p1(ndim) -// whose components contain the upper simplex borders -// (The simplex will not get reflected out of these borders !!!) -// -// It is possible to perform an initialization and -// a subsequent RunMinimization several times. -// Each time, a new class MHSimulatedAnnealing will get created -// (and destroyed). -// -Bool_t MSimulatedAnnealing::Initialize(const TMatrix &p, const TVector &y, - const TVector &p0, const TVector &p1) -{ - - fNdim = p.GetNcols(); - fMpts = p.GetNrows(); - - // - // many necessary checks ... - // - if (fMpts > gsMaxDim) - { - *fLog << err << "Dimension of Matrix fP is too big ... aborting." << endl; - return kFALSE; - } - if (fNdim+1 != fMpts) - { - *fLog << err << "Matrix fP does not have the right dimensions ... aborting." << endl; - return kFALSE; - } - if (y.GetNrows() != fMpts) - { - *fLog << err << "Array fY has not the right dimension ... aborting." << endl; - return kFALSE; - } - if (p0.GetNrows() != fNdim) - { - *fLog << err << "Array fP0 has not the right dimension ... aborting." << endl; - return kFALSE; - } - if (p1.GetNrows() != fNdim) - { - *fLog << err << "Array fP1 has not the right dimension ... aborting." << endl; - return kFALSE; - } - - // - // In order to allow multiple use of the class - // without need to construct the class every time new - // delete the old fResult and create a new one in RunMinimization - // - if (fResult) - delete fResult; - - fY.ResizeTo(fMpts); - - fPsum.ResizeTo(fNdim); - fPconv.ResizeTo(fNdim); - - fP0.ResizeTo(fNdim); - fP1.ResizeTo(fNdim); - fPb.ResizeTo(fNdim); - - fP.ResizeTo(fMpts,fNdim); - - fY = y; - fP = p; - fP0 = p0; - fP1 = p1; - fPconv.Zero(); - - fInit = kTRUE; - fYconv = 0.; - - return kTRUE; -} - - -// --------------------------------------------------------------------------- -// -// RunMinimization: -// -// Runs only eafter a call to Initialize(const TMatrix \&, const TVector \&, -// const TVector \&, const TVector \&) -// -// Temperature and number of moves should have been set -// (default: StartTemperature = 10, NumberOfMoves = 200 -// -// -// It is possible to perform an initialization and -// a subsequent RunMinimization several times. -// Each time, a new class MHSimulatedAnnealing will get created -// (and destroyed). -Bool_t MSimulatedAnnealing::RunMinimization() -{ - if (!fInit) - { - *fLog << err << "No succesful initialization performed yet... aborting." << endl; - return kFALSE; - } - - Int_t iter = 0; - UShort_t iret = 0; - UShort_t currentMove = 0; - Real_t currentTemp = fStartTemperature; - - fResult = new MHSimulatedAnnealing(fNumberOfMoves,fNdim); - if (fFullStorage) - fResult->InitFullSimplex(); - - while(1) - { - if (iter > 0) - { - *fLog << "Convergence at move: " << currentMove ; - *fLog << " and temperature: " << currentTemp << endl; - break; - } - - if (currentTemp > 0.) - { - // - // Reduce the temperature - // - // FIXME: Maybe it is necessary to also incorporate other - // ways to reduce the temperature (T0*(1-k/K)**alpha) - // - currentTemp = fStartTemperature*(1.-(float)currentMove++/fNumberOfMoves); - iter = 1; - } - else - { - // Make sure that now, the program will return only on convergence ! - // The program returns to here only after gsMaxStep moves - // If we have not reached convergence until then, we assume that an infinite - // loop has occurred and quit. - if (iret != 0) - { - *fLog << warn << "No Convergence at the end ! " << endl; - fY.Zero(); - - break; - } - iter = 150; - iret++; - currentMove++; - } - - if (fVerbose==2) { - *fLog << dbginf << " current..." << endl; - *fLog << " - move: " << currentMove << endl; - *fLog << " - temperature: " << currentTemp << endl; - *fLog << " - best function evaluation: " << fYb << endl; - } - - iter = Amebsa(iter, currentTemp); - - // Store the current best values in the histograms - fResult->StoreBestValueEver(fPb,fYb,currentMove); - - // Store the complete simplex if we have full storage - if (fFullStorage) - fResult->StoreFullSimplex(fP,currentMove); - } - - // - // Now, the matrizes and vectors have all the same value, - // Need to initialize again to allow a new Minimization - // - fInit = kFALSE; - - return kTRUE; -} - -// --------------------------------------------------------------------------- -// -// Amebsa -// -// This is the (adjusted) amebsa function from -// Numerical Recipies (pp. 457-458) -// -// The routine makes iter function evaluations at an annealing -// temperature fCurrentTemp, then returns. If iter is returned -// with a poisitive value, then early convergence has occurred. -// -Int_t MSimulatedAnnealing::Amebsa(Int_t iter, const Real_t temp) -{ - GetPsum(); - - while (1) - { - UShort_t ihi = 0; // Simplex point with highest function evaluation - UShort_t ilo = 1; // Simplex point with lowest function evaluation - - // Function eval. at ilo (with random fluctuations) - Real_t ylo = fY(0) + gRandom->Exp(temp); - - // Function eval. at ihi (with random fluctuations) - Real_t yhi = fY(1) + gRandom->Exp(temp); - - // The function evaluation at next highest point - Real_t ynhi = ylo; - - if (ylo > yhi) - { - // Determine which point is the highest (worst), - // next-highest and lowest (best) - ynhi = yhi; - yhi = ylo; - ylo = ynhi; - } - - // By looping over the points in the simplex - for (UShort_t i=2;iExp(temp); - - if (yt <= ylo) - { - ilo = i; - ylo = yt; - } - - if (yt > yhi) - { - ynhi = yhi; - ihi = i; - yhi = yt; - } - else - if (yt > ynhi) - ynhi = yt; - } - - // Now, fY(ilo) is smallest and fY(ihi) is at biggest value - if (iter < 0) - { - // Enough looping with this temperature, go to decrease it - // First put best point and value in slot 0 - - Real_t dum = fY(0); - fY(0) = fY(ilo); - fY(ilo) = dum; - - for (UShort_t n=0;n fP0(j)) - { - if (newcut > fP1(j)) - { - ptry(j) = fP1(j); - borderflag = 1; - } - else - if (newcut < fP0(j)) - { - ptry(j) = fP0(j); - borderflag = 1; - } - } - - else - { - if (newcut < fP1(j)) - { - ptry(j) = fP1(j); - borderflag = 1; - } - else - if (newcut > fP0(j)) - { - ptry(j) = fP0(j); - borderflag = 1; - } - } - } - - Real_t faccompare = 0.5; - Real_t ytry = 0; - - switch (borderflag) - { - case kENoBorder: - ytry = FunctionToMinimize(fPsum); - break; - - case kEStrictBorder: - ytry = FunctionToMinimize(fPsum) + gsYtryStr; - break; - - case kEContractBorder: - ytry = fac == faccompare ? gsYtryCon : gsYtryStr; - break; - } - - if (ytry < fYb) - { - fPb = ptry; - fYb = ytry; - } - - const Real_t yflu = ytry + gRandom->Exp(temp); - - if (yflu >= yhi) - return yflu; - - fY(ihi) = ytry; - yhi = yflu; - - for(Int_t j=0; j -#endif - -class MHSimulatedAnnealing; -class TRandom; - -class MSimulatedAnnealing : public MParContainer -{ -private: - - static const Float_t gsYtryStr; // Fixed high value to keep the simplex inside the borders - static const Float_t gsYtryCon; // Fixed high value to keep the simplex inside the borders - static const Int_t gsMaxDim; // Fixed maximum number of dimensions - static const Int_t gsMaxStep; // Fixed maximum number of loops with temperature=0 - - MHSimulatedAnnealing *fResult; //! The histogram output container - - TRandom *fRandom; // The random number generator -> random numbers between 0 and 1 - - Real_t fTolerance; // The convergence break condition - - UShort_t fNdim; // The number of parameters - UShort_t fMpts; // The number of simplex points (=fNdim+1) - - UShort_t fNumberOfMoves; // The total number of moves (== CPU time) - - Real_t fStartTemperature; // The start temperature -> will slowly get decreased to 0 - - Bool_t fFullStorage; // kTRUE -> the whole simplex gets stored in MHSimlutedAnnealing - Bool_t fInit; // kTRUE -> initialization was succesful - - TMatrix fP; // The (ndim+1,ndim) matrix containing the simplex - - TVector fPsum; // The sum of each point of the simplex - - TVector fP0; // The boundary conditions on the weak side - TVector fP1; // The boundary conditions on the strong side - TVector fY; // The array containing the function evaluation results - - Real_t fYb; // The best function evaluation value ever found - Real_t fYconv; // The function evaluation value at the convergence point - - TVector fPb; // The parameters belonging to fYb - TVector fPconv; // The parameters belonging to fYconv - - Int_t Amebsa(Int_t iter, - const Real_t temp); // The function deciding if the simplex has to get reflected, expanded or contracted - - Real_t Amotsa(const Float_t fac, - const UShort_t ihi, - Real_t &yhi, - const Real_t temp); // The function reflecting, expanding and contracting the simplex: fac=-1 -> reflection, fac=0.5 -> contraction, fac=2.0 -> expansion - - void GetPsum(); - -protected: - - virtual Float_t FunctionToMinimize(const TVector &arr); // The optimization function - -public: - enum BorderFlag_t { kENoBorder, kEStrictBorder, kEContractBorder }; - enum Verbosity_t { kEDefault, kEVerbose, kEDebug }; - -private: - BorderFlag_t fBorder; - Verbosity_t fVerbose; - -public: - MSimulatedAnnealing(); - virtual ~MSimulatedAnnealing(); - - void ModifyTolerance(Float_t tol) { fTolerance = tol; } - void ModifyBorderFlag(BorderFlag_t border) { fBorder = border; } - - Bool_t Initialize(const TMatrix &p, const TVector &y, - const TVector &p0, const TVector &p1); - - void SetNumberOfMoves(UShort_t moves) { fNumberOfMoves = moves; } - void SetStartTemperature(Float_t temp) { fStartTemperature = temp; } - void SetFullStorage() { fFullStorage = kTRUE; } - void SetVerbosityLevel(Verbosity_t level) { fVerbose = level; } - void SetRandom(TRandom *rand) { fRandom = rand; } - - const TVector &GetPb() const { return fPb; } - Float_t GetYb() const { return fYb; } - - const TVector &GetPconv() const { return fPconv; } - Float_t GetYconv() const { return fYconv; } - - - MHSimulatedAnnealing *GetResult() { return fResult; } - - Bool_t RunMinimization(); - - ClassDef(MSimulatedAnnealing,1) // Class to perform a Simulated Annealing Minimization -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/MHSimulatedAnnealing.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSimulatedAnnealing.cc (revision 2816) +++ (revision ) @@ -1,249 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Markus Gaug -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -/////////////////////////////////////////////////////////////////////// -// -// MHSimulatedAnnealing -// -// This class contains different histograms of the Simulated Annealing -// Snapshort during optimization and the final results -// -/////////////////////////////////////////////////////////////////////// -#include "MHSimulatedAnnealing.h" - -#include -#include - -#include -#include - -#include "MBinning.h" - -ClassImp(MHSimulatedAnnealing); - -// -------------------------------------------------------------------------- -// -// Setup histograms -// -MHSimulatedAnnealing::MHSimulatedAnnealing(UShort_t moves, UShort_t ndim, - const char *name, - const char *title) - : fDim(ndim), fMoves(moves), fTimeEvolution(NULL), fBestEver(), fBestFuncEval() -{ - - // - // set the name and title of this object - // - fName = name ? name : "MHSimulatedAnnealing"; - fTitle = title ? title : "Output Histograms of a Simulated Annealing Run"; - - fBestEver.SetName("Hist_BestEver"); - fBestEver.SetTitle("Best Param. Combinations"); - fBestEver.SetXTitle("Run Duration"); - fBestEver.SetYTitle("Parameter Nr."); - fBestEver.SetZTitle("Parameter Value"); - fBestEver.SetDirectory(NULL); - - fBestFuncEval.SetName("Hist_BestFuncEval"); - fBestFuncEval.SetTitle("Best Function Evaluation"); - fBestFuncEval.SetXTitle("Run Duration"); - fBestFuncEval.SetYTitle("Function Value"); - fBestFuncEval.SetDirectory(NULL); - - MBinning binsx, binsy; - binsx.SetEdges(fMoves+1, 0, 1); - binsy.SetEdges(fDim, 0.5, fDim+0.5); - MH::SetBinning(&fBestEver, &binsx, &binsy); - - // For better visibility, omit the first entry in fBestFuncEval - // It has nothing significant, anyway - binsx.SetEdges(fMoves,1./(fMoves+1), 1); - binsx.Apply(fBestFuncEval); - -} - -void MHSimulatedAnnealing::InitFullSimplex() -{ - if (fTimeEvolution) - delete fTimeEvolution; - - fTimeEvolution = new TObjArray; - fTimeEvolution->SetOwner(); - - for (Int_t i=0;iSetXTitle("Run Duration"); - hist->SetYTitle("Point Nr. Simplex"); - hist->SetZTitle(Form("Value of Parameter %d",i)); - fTimeEvolution->Add(hist); - } -} - -Bool_t MHSimulatedAnnealing::StoreFullSimplex(const TMatrix &p, const UShort_t move) -{ - - if (!fTimeEvolution) - return kFALSE; - - Int_t idx=0; - const Axis_t bin = (move-0.5)/(fMoves+1); - - TIter Next(fTimeEvolution); - TH2F *hist=NULL; - while ((hist=(TH2F*)Next())) - { - for (Int_t i=0;iFill(bin,i,p(i,idx)); - idx++; - } - return kTRUE; -} - -Bool_t MHSimulatedAnnealing::StoreBestValueEver(const TVector &y, const Float_t yb, const UShort_t move) -{ - if (y.GetNrows() != fDim) - return kFALSE; - - const Axis_t bin = (move-0.5)/(fMoves+1); - - for (Int_t i=0;iAt(index))) - return kFALSE; - - hist->SetNameTitle(Form("Hist_%s",title),title); - hist->SetYTitle(Form("Value of Parameter %s",title)); - - return kTRUE; -} - -void MHSimulatedAnnealing::ChangeFuncTitle(const char* title) -{ - fBestFuncEval.SetTitle(title); -} - -TObject *MHSimulatedAnnealing::DrawClone(Option_t *opt) const -{ - UShort_t i=2; - - TCanvas *c = MakeDefCanvas(this, 720, 810); - if (fTimeEvolution) - c->Divide(2,(int)(fDim/2.)+1); - else - gPad->Divide(1,2); - - gROOT->SetSelectedPad(NULL); - - c->cd(1); - gStyle->SetOptStat(0); - ((TH1&)fBestFuncEval).DrawCopy(); - - c->cd(2); - gStyle->SetOptStat(10); - ((TH2&)fBestEver).DrawCopy(opt); - - if (fTimeEvolution) - { - TH2F *hist = NULL; - TIter Next(fTimeEvolution); - while ((hist=(TH2F*)Next())) - { - c->cd(++i); - hist->DrawCopy(opt); - } - } - c->Modified(); - c->Update(); - - return c; -} - - - -// -------------------------------------------------------------------------- -// -// Draw all histograms. -// -void MHSimulatedAnnealing::Draw(Option_t *opt) -{ - UShort_t i=2; - - if (!gPad) - MakeDefCanvas(this,780,940); - - if (fTimeEvolution) - gPad->Divide(2,(int)(fDim/2.)+1); - else - gPad->Divide(1,2); - - gPad->cd(1); - gStyle->SetOptStat(0); - fBestFuncEval.Draw(); - gPad->Modified(); - gPad->Update(); - - gPad->cd(2); - gStyle->SetOptStat(10); - fBestEver.Draw(opt); - gPad->Modified(); - gPad->Update(); - - if (!fTimeEvolution) - return; - - TH2F *hist = NULL; - TIter Next(fTimeEvolution); - while ((hist=(TH2F*)Next())) - { - gPad->cd(++i); - hist->Draw(opt); - } - gPad->Modified(); - gPad->Update(); -} - -// -------------------------------------------------------------------------- -// -// Delete the histograms. -// -MHSimulatedAnnealing::~MHSimulatedAnnealing() -{ - if (fTimeEvolution) - delete fTimeEvolution; -} - Index: trunk/MagicSoft/Mars/mhist/MHSimulatedAnnealing.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHSimulatedAnnealing.h (revision 2816) +++ (revision ) @@ -1,51 +1,0 @@ -#ifndef MARS_MHSimulatedAnnealing -#define MARS_MHSimulatedAnnealing -/////////////////////////////////////////////////////////////////////////////// -// -// MHSimulatedAnnealing -// -// Output container of MSimulatedAnnealing -/////////////////////////////////////////////////////////////////////////////// -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TH2 -#include -#endif - -class MHSimulatedAnnealing : public MH -{ -private: - UShort_t fDim; // The dimension of the whole thing - UShort_t fMoves; // The total number of moves - - TObjArray *fTimeEvolution; //-> Display the time evolution of the simplex in TH1D's - - TH2F fBestEver; // The best values ever found during search - TH1F fBestFuncEval; // The best function values ever found during search - -public: - - MHSimulatedAnnealing(UShort_t moves = 0,UShort_t ndim = 0, - const char *name=NULL, const char *title=NULL); - ~MHSimulatedAnnealing(); - - void InitFullSimplex(); - Bool_t StoreFullSimplex(const TMatrix &p, const UShort_t move); - Bool_t StoreBestValueEver(const TVector &y, const Float_t yb, const UShort_t move); - - Bool_t ChangeTitle(const UShort_t index, const char* title); - void ChangeFuncTitle(const char* title); - - TObjArray *GetTimeEvolution() const { return fTimeEvolution; } - const TH2F &GetBestEver() const { return fBestEver; } - const TH1F &GetBestFuncEval() const { return fBestFuncEval; } - - void Draw(Option_t *opt=NULL); - TObject *DrawClone(Option_t *opt=NULL) const; - - ClassDef(MHSimulatedAnnealing,1) // Storage Histogram Container for Cuteval Results -}; - -#endif Index: trunk/MagicSoft/Mars/mhist/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhist/Makefile (revision 2816) +++ trunk/MagicSoft/Mars/mhist/Makefile (revision 2817) @@ -59,6 +59,5 @@ MHCT1Supercuts.cc \ MHCamera.cc \ - MHSupercuts.cc \ - MHSimulatedAnnealing.cc + MHSupercuts.cc # MHCurrents.cc \ Index: trunk/MagicSoft/Mars/mtools/MHSimulatedAnnealing.cc =================================================================== --- trunk/MagicSoft/Mars/mtools/MHSimulatedAnnealing.cc (revision 2817) +++ trunk/MagicSoft/Mars/mtools/MHSimulatedAnnealing.cc (revision 2817) @@ -0,0 +1,249 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Markus Gaug +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////// +// +// MHSimulatedAnnealing +// +// This class contains different histograms of the Simulated Annealing +// Snapshort during optimization and the final results +// +/////////////////////////////////////////////////////////////////////// +#include "MHSimulatedAnnealing.h" + +#include +#include + +#include +#include + +#include "MBinning.h" + +ClassImp(MHSimulatedAnnealing); + +// -------------------------------------------------------------------------- +// +// Setup histograms +// +MHSimulatedAnnealing::MHSimulatedAnnealing(UShort_t moves, UShort_t ndim, + const char *name, + const char *title) + : fDim(ndim), fMoves(moves), fTimeEvolution(NULL), fBestEver(), fBestFuncEval() +{ + + // + // set the name and title of this object + // + fName = name ? name : "MHSimulatedAnnealing"; + fTitle = title ? title : "Output Histograms of a Simulated Annealing Run"; + + fBestEver.SetName("Hist_BestEver"); + fBestEver.SetTitle("Best Param. Combinations"); + fBestEver.SetXTitle("Run Duration"); + fBestEver.SetYTitle("Parameter Nr."); + fBestEver.SetZTitle("Parameter Value"); + fBestEver.SetDirectory(NULL); + + fBestFuncEval.SetName("Hist_BestFuncEval"); + fBestFuncEval.SetTitle("Best Function Evaluation"); + fBestFuncEval.SetXTitle("Run Duration"); + fBestFuncEval.SetYTitle("Function Value"); + fBestFuncEval.SetDirectory(NULL); + + MBinning binsx, binsy; + binsx.SetEdges(fMoves+1, 0, 1); + binsy.SetEdges(fDim, 0.5, fDim+0.5); + MH::SetBinning(&fBestEver, &binsx, &binsy); + + // For better visibility, omit the first entry in fBestFuncEval + // It has nothing significant, anyway + binsx.SetEdges(fMoves,1./(fMoves+1), 1); + binsx.Apply(fBestFuncEval); + +} + +void MHSimulatedAnnealing::InitFullSimplex() +{ + if (fTimeEvolution) + delete fTimeEvolution; + + fTimeEvolution = new TObjArray; + fTimeEvolution->SetOwner(); + + for (Int_t i=0;iSetXTitle("Run Duration"); + hist->SetYTitle("Point Nr. Simplex"); + hist->SetZTitle(Form("Value of Parameter %d",i)); + fTimeEvolution->Add(hist); + } +} + +Bool_t MHSimulatedAnnealing::StoreFullSimplex(const TMatrix &p, const UShort_t move) +{ + + if (!fTimeEvolution) + return kFALSE; + + Int_t idx=0; + const Axis_t bin = (move-0.5)/(fMoves+1); + + TIter Next(fTimeEvolution); + TH2F *hist=NULL; + while ((hist=(TH2F*)Next())) + { + for (Int_t i=0;iFill(bin,i,p(i,idx)); + idx++; + } + return kTRUE; +} + +Bool_t MHSimulatedAnnealing::StoreBestValueEver(const TVector &y, const Float_t yb, const UShort_t move) +{ + if (y.GetNrows() != fDim) + return kFALSE; + + const Axis_t bin = (move-0.5)/(fMoves+1); + + for (Int_t i=0;iAt(index))) + return kFALSE; + + hist->SetNameTitle(Form("Hist_%s",title),title); + hist->SetYTitle(Form("Value of Parameter %s",title)); + + return kTRUE; +} + +void MHSimulatedAnnealing::ChangeFuncTitle(const char* title) +{ + fBestFuncEval.SetTitle(title); +} + +TObject *MHSimulatedAnnealing::DrawClone(Option_t *opt) const +{ + UShort_t i=2; + + TCanvas *c = MakeDefCanvas(this, 720, 810); + if (fTimeEvolution) + c->Divide(2,(int)(fDim/2.)+1); + else + gPad->Divide(1,2); + + gROOT->SetSelectedPad(NULL); + + c->cd(1); + gStyle->SetOptStat(0); + ((TH1&)fBestFuncEval).DrawCopy(); + + c->cd(2); + gStyle->SetOptStat(10); + ((TH2&)fBestEver).DrawCopy(opt); + + if (fTimeEvolution) + { + TH2F *hist = NULL; + TIter Next(fTimeEvolution); + while ((hist=(TH2F*)Next())) + { + c->cd(++i); + hist->DrawCopy(opt); + } + } + c->Modified(); + c->Update(); + + return c; +} + + + +// -------------------------------------------------------------------------- +// +// Draw all histograms. +// +void MHSimulatedAnnealing::Draw(Option_t *opt) +{ + UShort_t i=2; + + if (!gPad) + MakeDefCanvas(this,780,940); + + if (fTimeEvolution) + gPad->Divide(2,(int)(fDim/2.)+1); + else + gPad->Divide(1,2); + + gPad->cd(1); + gStyle->SetOptStat(0); + fBestFuncEval.Draw(); + gPad->Modified(); + gPad->Update(); + + gPad->cd(2); + gStyle->SetOptStat(10); + fBestEver.Draw(opt); + gPad->Modified(); + gPad->Update(); + + if (!fTimeEvolution) + return; + + TH2F *hist = NULL; + TIter Next(fTimeEvolution); + while ((hist=(TH2F*)Next())) + { + gPad->cd(++i); + hist->Draw(opt); + } + gPad->Modified(); + gPad->Update(); +} + +// -------------------------------------------------------------------------- +// +// Delete the histograms. +// +MHSimulatedAnnealing::~MHSimulatedAnnealing() +{ + if (fTimeEvolution) + delete fTimeEvolution; +} + Index: trunk/MagicSoft/Mars/mtools/MHSimulatedAnnealing.h =================================================================== --- trunk/MagicSoft/Mars/mtools/MHSimulatedAnnealing.h (revision 2817) +++ trunk/MagicSoft/Mars/mtools/MHSimulatedAnnealing.h (revision 2817) @@ -0,0 +1,51 @@ +#ifndef MARS_MHSimulatedAnnealing +#define MARS_MHSimulatedAnnealing +/////////////////////////////////////////////////////////////////////////////// +// +// MHSimulatedAnnealing +// +// Output container of MSimulatedAnnealing +/////////////////////////////////////////////////////////////////////////////// +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH2 +#include +#endif + +class MHSimulatedAnnealing : public MH +{ +private: + UShort_t fDim; // The dimension of the whole thing + UShort_t fMoves; // The total number of moves + + TObjArray *fTimeEvolution; //-> Display the time evolution of the simplex in TH1D's + + TH2F fBestEver; // The best values ever found during search + TH1F fBestFuncEval; // The best function values ever found during search + +public: + + MHSimulatedAnnealing(UShort_t moves = 0,UShort_t ndim = 0, + const char *name=NULL, const char *title=NULL); + ~MHSimulatedAnnealing(); + + void InitFullSimplex(); + Bool_t StoreFullSimplex(const TMatrix &p, const UShort_t move); + Bool_t StoreBestValueEver(const TVector &y, const Float_t yb, const UShort_t move); + + Bool_t ChangeTitle(const UShort_t index, const char* title); + void ChangeFuncTitle(const char* title); + + TObjArray *GetTimeEvolution() const { return fTimeEvolution; } + const TH2F &GetBestEver() const { return fBestEver; } + const TH1F &GetBestFuncEval() const { return fBestFuncEval; } + + void Draw(Option_t *opt=NULL); + TObject *DrawClone(Option_t *opt=NULL) const; + + ClassDef(MHSimulatedAnnealing,1) // Storage Histogram Container for Cuteval Results +}; + +#endif Index: trunk/MagicSoft/Mars/mtools/MSimulatedAnnealing.cc =================================================================== --- trunk/MagicSoft/Mars/mtools/MSimulatedAnnealing.cc (revision 2817) +++ trunk/MagicSoft/Mars/mtools/MSimulatedAnnealing.cc (revision 2817) @@ -0,0 +1,628 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Markus Gaug 10/2002 +! +! Copyright: MAGIC Software Development, 2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MSimulatedAnnealing // +// // +// class to perform a Simulated Annealing minimization on an n-dimensional // +// simplex of a function 'FunctionToMinimize(TArrayF &)' in multi- // +// dimensional parameter space. // +// (The code is adapted from Numerical Recipies in C++, 2nd ed., // +// pp. 457-459) // +// // +// Classes can inherit from MSimulatedAnnealing // +// and use the function: // +// // +// RunSimulatedAnnealing(); // +// // +// They HAVE TO initialize the following input arguments // +// (with ndim being the parameter dimension (max. 20)): // +// // +// 1) a TMatrix p(ndim+1,ndim) // +// holding the start simplex // +// 2) a TArrayF y(ndim+1) // +// whose components must be pre-initialized to the values of // +// FunctionToMinimize evaluated at the fNdim+1 vertices (rows) of fP // +// 3) a TArrayF p0(ndim) // +// whose components contain the lower simplex borders // +// 4) a TArrayF p1(ndim) // +// whose components contain the upper simplex borders // +// (The simplex will not get reflected out of these borders !!!) // +// // +// These arrays have to be initialized with a call to: // +// Initialize(TMatrix \&, TArrayF \&, TArrayF \&, TArrayF \&) // +// // +// 5) a virtual function FunctionToMinimize(TArrayF &) // +// acting on a TArrayF(ndim) array of parameter values // +// // +// Additionally, a global start temperature can be chosen with: // +// // +// SetStartTemperature(Float_t temp) // +// (default is: 10) // +// // +// A total number of total moves (watch out for the CPU time!!!) with: // +// // +// SetNumberOfMoves(Float_t totalMoves) // +// (default is: 200) // +// // +// The temperature is reduced after evaluation step like: // +// CurrentTemperature = StartTemperature*(1-currentMove/totalMoves) // +// where currentMove is the cumulative number of moves so far // +// // +// WARNING: The start temperature and number of moves has to be optimized // +// for each individual problem. // +// It is not straightforward using the defaults! // +// In case, you omit this important step, // +// you will get local minima without even noticing it!! // +// // +// You may define the following variables: // +// // +// 1) A global convergence criterium fTol // +// which determines an early return for: // +// // +// max(FunctionToMinimize(p))-min(FunctionToMinimize(p)) // +// ----------------------------------------------------- \< fTol // +// max(FunctionToMinimize(p))+min(FunctionToMinimize(p)) // +// // +// ModifyTolerance(Float_t) // +// // +// 2) A verbose level for prints to *fLog // +// // +// SetVerbosityLevel(Verbosity_t) // +// // +// 3) A bit if you want to have stored // +// the full simplex after every call to Amebsa: // +// // +// SetFullStorage() // +// // +// 4) The random number generator // +// e.g. if you want to test the stability of the output // +// // +// SetRandom(TRandom *rand) // +// // +// // +// Output containers: // +// // +// MHSimulatedAnnealing // +// // +// Use: // +// GetResult()->Draw(Option_t *o) // +// or // +// GetResult()->DrawClone(Option_t *o) // +// // +// to retrieve the output histograms // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MSimulatedAnnealing.h" +#include "MHSimulatedAnnealing.h" + +#include +#include + +#include + +#include "MLog.h" +#include "MLogManip.h" + +const Float_t MSimulatedAnnealing::gsYtryStr = 10000000; +const Float_t MSimulatedAnnealing::gsYtryCon = 20000000; +const Int_t MSimulatedAnnealing::gsMaxDim = 20; +const Int_t MSimulatedAnnealing::gsMaxStep = 50; + +ClassImp(MSimulatedAnnealing); + +using namespace std; + +// --------------------------------------------------------------------------- +// +// Default Constructor +// Initializes random number generator and default variables +// +MSimulatedAnnealing::MSimulatedAnnealing() + : fResult(NULL), fTolerance(0.0001), + fNdim(0), fNumberOfMoves(200), + fStartTemperature(10), fFullStorage(kFALSE), + fInit(kFALSE), + fP(gsMaxDim, gsMaxDim), fP0(gsMaxDim), + fP1(gsMaxDim), fY(gsMaxDim), fYb(gsMaxDim), fYconv(gsMaxDim), + fPb(gsMaxDim), fPconv(gsMaxDim), + fBorder(kEStrictBorder), fVerbose(kEDefault) +{ + + // random number generator + fRandom = gRandom; +} + +// -------------------------------------------------------------------------- +// +// Destructor. +// +MSimulatedAnnealing::~MSimulatedAnnealing() +{ + if (fResult) + delete fResult; +} + +// --------------------------------------------------------------------------- +// +// Initialization needs the following four members: +// +// 1) a TMatrix p(ndim+1,ndim) +// holding the start simplex +// 2) a TVector y(ndim+1) +// whose components must be pre-initialized to the values of +// FunctionToMinimize evaluated at the fNdim+1 vertices (rows) of fP +// 3) a TVector p0(ndim) +// whose components contain the lower simplex borders +// 4) a TVector p1(ndim) +// whose components contain the upper simplex borders +// (The simplex will not get reflected out of these borders !!!) +// +// It is possible to perform an initialization and +// a subsequent RunMinimization several times. +// Each time, a new class MHSimulatedAnnealing will get created +// (and destroyed). +// +Bool_t MSimulatedAnnealing::Initialize(const TMatrix &p, const TVector &y, + const TVector &p0, const TVector &p1) +{ + + fNdim = p.GetNcols(); + fMpts = p.GetNrows(); + + // + // many necessary checks ... + // + if (fMpts > gsMaxDim) + { + gLog << err << "Dimension of Matrix fP is too big ... aborting." << endl; + return kFALSE; + } + if (fNdim+1 != fMpts) + { + gLog << err << "Matrix fP does not have the right dimensions ... aborting." << endl; + return kFALSE; + } + if (y.GetNrows() != fMpts) + { + gLog << err << "Array fY has not the right dimension ... aborting." << endl; + return kFALSE; + } + if (p0.GetNrows() != fNdim) + { + gLog << err << "Array fP0 has not the right dimension ... aborting." << endl; + return kFALSE; + } + if (p1.GetNrows() != fNdim) + { + gLog << err << "Array fP1 has not the right dimension ... aborting." << endl; + return kFALSE; + } + + // + // In order to allow multiple use of the class + // without need to construct the class every time new + // delete the old fResult and create a new one in RunMinimization + // + if (fResult) + delete fResult; + + fY.ResizeTo(fMpts); + + fPsum.ResizeTo(fNdim); + fPconv.ResizeTo(fNdim); + + fP0.ResizeTo(fNdim); + fP1.ResizeTo(fNdim); + fPb.ResizeTo(fNdim); + + fP.ResizeTo(fMpts,fNdim); + + fY = y; + fP = p; + fP0 = p0; + fP1 = p1; + fPconv.Zero(); + + fInit = kTRUE; + fYconv = 0.; + + return kTRUE; +} + + +// --------------------------------------------------------------------------- +// +// RunMinimization: +// +// Runs only eafter a call to Initialize(const TMatrix \&, const TVector \&, +// const TVector \&, const TVector \&) +// +// Temperature and number of moves should have been set +// (default: StartTemperature = 10, NumberOfMoves = 200 +// +// +// It is possible to perform an initialization and +// a subsequent RunMinimization several times. +// Each time, a new class MHSimulatedAnnealing will get created +// (and destroyed). +Bool_t MSimulatedAnnealing::RunMinimization() +{ + if (!fInit) + { + gLog << err << "No succesful initialization performed yet... aborting." << endl; + return kFALSE; + } + + Int_t iter = 0; + UShort_t iret = 0; + UShort_t currentMove = 0; + Real_t currentTemp = fStartTemperature; + + fResult = new MHSimulatedAnnealing(fNumberOfMoves,fNdim); + if (fFullStorage) + fResult->InitFullSimplex(); + + while(1) + { + if (iter > 0) + { + gLog << "Convergence at move: " << currentMove ; + gLog << " and temperature: " << currentTemp << endl; + break; + } + + if (currentTemp > 0.) + { + // + // Reduce the temperature + // + // FIXME: Maybe it is necessary to also incorporate other + // ways to reduce the temperature (T0*(1-k/K)**alpha) + // + currentTemp = fStartTemperature*(1.-(float)currentMove++/fNumberOfMoves); + iter = 1; + } + else + { + // Make sure that now, the program will return only on convergence ! + // The program returns to here only after gsMaxStep moves + // If we have not reached convergence until then, we assume that an infinite + // loop has occurred and quit. + if (iret != 0) + { + gLog << warn << "No Convergence at the end ! " << endl; + fY.Zero(); + + break; + } + iter = 150; + iret++; + currentMove++; + } + + if (fVerbose==2) { + gLog << dbginf << " current..." << endl; + gLog << " - move: " << currentMove << endl; + gLog << " - temperature: " << currentTemp << endl; + gLog << " - best function evaluation: " << fYb << endl; + } + + iter = Amebsa(iter, currentTemp); + + // Store the current best values in the histograms + fResult->StoreBestValueEver(fPb,fYb,currentMove); + + // Store the complete simplex if we have full storage + if (fFullStorage) + fResult->StoreFullSimplex(fP,currentMove); + } + + // + // Now, the matrizes and vectors have all the same value, + // Need to initialize again to allow a new Minimization + // + fInit = kFALSE; + + return kTRUE; +} + +// --------------------------------------------------------------------------- +// +// Amebsa +// +// This is the (adjusted) amebsa function from +// Numerical Recipies (pp. 457-458) +// +// The routine makes iter function evaluations at an annealing +// temperature fCurrentTemp, then returns. If iter is returned +// with a poisitive value, then early convergence has occurred. +// +Int_t MSimulatedAnnealing::Amebsa(Int_t iter, const Real_t temp) +{ + GetPsum(); + + while (1) + { + UShort_t ihi = 0; // Simplex point with highest function evaluation + UShort_t ilo = 1; // Simplex point with lowest function evaluation + + // Function eval. at ilo (with random fluctuations) + Real_t ylo = fY(0) + gRandom->Exp(temp); + + // Function eval. at ihi (with random fluctuations) + Real_t yhi = fY(1) + gRandom->Exp(temp); + + // The function evaluation at next highest point + Real_t ynhi = ylo; + + if (ylo > yhi) + { + // Determine which point is the highest (worst), + // next-highest and lowest (best) + ynhi = yhi; + yhi = ylo; + ylo = ynhi; + } + + // By looping over the points in the simplex + for (UShort_t i=2;iExp(temp); + + if (yt <= ylo) + { + ilo = i; + ylo = yt; + } + + if (yt > yhi) + { + ynhi = yhi; + ihi = i; + yhi = yt; + } + else + if (yt > ynhi) + ynhi = yt; + } + + // Now, fY(ilo) is smallest and fY(ihi) is at biggest value + if (iter < 0) + { + // Enough looping with this temperature, go to decrease it + // First put best point and value in slot 0 + + Real_t dum = fY(0); + fY(0) = fY(ilo); + fY(ilo) = dum; + + for (UShort_t n=0;n fP0(j)) + { + if (newcut > fP1(j)) + { + ptry(j) = fP1(j); + borderflag = 1; + } + else + if (newcut < fP0(j)) + { + ptry(j) = fP0(j); + borderflag = 1; + } + } + + else + { + if (newcut < fP1(j)) + { + ptry(j) = fP1(j); + borderflag = 1; + } + else + if (newcut > fP0(j)) + { + ptry(j) = fP0(j); + borderflag = 1; + } + } + } + + Real_t faccompare = 0.5; + Real_t ytry = 0; + + switch (borderflag) + { + case kENoBorder: + ytry = FunctionToMinimize(fPsum); + break; + + case kEStrictBorder: + ytry = FunctionToMinimize(fPsum) + gsYtryStr; + break; + + case kEContractBorder: + ytry = fac == faccompare ? gsYtryCon : gsYtryStr; + break; + } + + if (ytry < fYb) + { + fPb = ptry; + fYb = ytry; + } + + const Real_t yflu = ytry + gRandom->Exp(temp); + + if (yflu >= yhi) + return yflu; + + fY(ihi) = ytry; + yhi = yflu; + + for(Int_t j=0; j +#endif + +class MHSimulatedAnnealing; +class TRandom; + +class MSimulatedAnnealing : public TObject +{ +private: + + static const Float_t gsYtryStr; // Fixed high value to keep the simplex inside the borders + static const Float_t gsYtryCon; // Fixed high value to keep the simplex inside the borders + static const Int_t gsMaxDim; // Fixed maximum number of dimensions + static const Int_t gsMaxStep; // Fixed maximum number of loops with temperature=0 + + MHSimulatedAnnealing *fResult; //! The histogram output container + + TRandom *fRandom; // The random number generator -> random numbers between 0 and 1 + + Real_t fTolerance; // The convergence break condition + + UShort_t fNdim; // The number of parameters + UShort_t fMpts; // The number of simplex points (=fNdim+1) + + UShort_t fNumberOfMoves; // The total number of moves (== CPU time) + + Real_t fStartTemperature; // The start temperature -> will slowly get decreased to 0 + + Bool_t fFullStorage; // kTRUE -> the whole simplex gets stored in MHSimlutedAnnealing + Bool_t fInit; // kTRUE -> initialization was succesful + + TMatrix fP; // The (ndim+1,ndim) matrix containing the simplex + + TVector fPsum; // The sum of each point of the simplex + + TVector fP0; // The boundary conditions on the weak side + TVector fP1; // The boundary conditions on the strong side + TVector fY; // The array containing the function evaluation results + + Real_t fYb; // The best function evaluation value ever found + Real_t fYconv; // The function evaluation value at the convergence point + + TVector fPb; // The parameters belonging to fYb + TVector fPconv; // The parameters belonging to fYconv + + Int_t Amebsa(Int_t iter, + const Real_t temp); // The function deciding if the simplex has to get reflected, expanded or contracted + + Real_t Amotsa(const Float_t fac, + const UShort_t ihi, + Real_t &yhi, + const Real_t temp); // The function reflecting, expanding and contracting the simplex: fac=-1 -> reflection, fac=0.5 -> contraction, fac=2.0 -> expansion + + void GetPsum(); + +protected: + + virtual Float_t FunctionToMinimize(const TVector &arr); // The optimization function + +public: + enum BorderFlag_t { kENoBorder, kEStrictBorder, kEContractBorder }; + enum Verbosity_t { kEDefault, kEVerbose, kEDebug }; + +private: + BorderFlag_t fBorder; + Verbosity_t fVerbose; + +public: + MSimulatedAnnealing(); + virtual ~MSimulatedAnnealing(); + + void ModifyTolerance(Float_t tol) { fTolerance = tol; } + void ModifyBorderFlag(BorderFlag_t border) { fBorder = border; } + + Bool_t Initialize(const TMatrix &p, const TVector &y, + const TVector &p0, const TVector &p1); + + void SetNumberOfMoves(UShort_t moves) { fNumberOfMoves = moves; } + void SetStartTemperature(Float_t temp) { fStartTemperature = temp; } + void SetFullStorage() { fFullStorage = kTRUE; } + void SetVerbosityLevel(Verbosity_t level) { fVerbose = level; } + void SetRandom(TRandom *rand) { fRandom = rand; } + + const TVector &GetPb() const { return fPb; } + Float_t GetYb() const { return fYb; } + + const TVector &GetPconv() const { return fPconv; } + Float_t GetYconv() const { return fYconv; } + + + MHSimulatedAnnealing *GetResult() { return fResult; } + + Bool_t RunMinimization(); + + ClassDef(MSimulatedAnnealing,1) // Class to perform a Simulated Annealing Minimization +}; + +#endif