Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 2680) +++ /trunk/MagicSoft/Mars/Changelog (revision 2681) @@ -5,4 +5,19 @@ -*-*- END OF LINE -*-*- + 2003/12/12: Markus Gaug + + * manalysis/MSimulatedAnnealing.[h,cc] + * mhist/MHSimulatedAnnealing.[h,cc] + - new classes to do a minimization after the Simulated Annealing + procedure + Please do make dox and look into the class documentation to know + how to use it + + * manalysis/Makefile + * mhist/Makefile + * manalysis/AnalysisLinkDef.h + * mhist/HistLinkDef.h + - added the SimulatedAnnealing Classes + 2003/12/12: Markus Gaug / Michele Doro Index: /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2680) +++ /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2681) @@ -83,4 +83,6 @@ #pragma link C++ class MMcCalibrationCalc+; +#pragma link C++ class MSimulatedAnnealing+; + #endif Index: /trunk/MagicSoft/Mars/manalysis/MSimulatedAnnealing.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSimulatedAnnealing.cc (revision 2681) +++ /trunk/MagicSoft/Mars/manalysis/MSimulatedAnnealing.cc (revision 2681) @@ -0,0 +1,625 @@ +/* ======================================================================== *\ +! +! * +! * 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" + +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 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/manalysis/Makefile =================================================================== --- /trunk/MagicSoft/Mars/manalysis/Makefile (revision 2680) +++ /trunk/MagicSoft/Mars/manalysis/Makefile (revision 2681) @@ -86,5 +86,6 @@ MArrivalTime.cc \ MArrivalTimeCalc.cc \ - MMcCalibrationCalc.cc + MMcCalibrationCalc.cc \ + MSimulatedAnnealing.cc SRCS = $(SRCFILES) Index: /trunk/MagicSoft/Mars/mhist/HistLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 2680) +++ /trunk/MagicSoft/Mars/mhist/HistLinkDef.h (revision 2681) @@ -56,4 +56,6 @@ #pragma link C++ class MHCalibrationPixel+; +#pragma link C++ class MHSimulatedAnnealing+; + #endif Index: /trunk/MagicSoft/Mars/mhist/MHSimulatedAnnealing.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHSimulatedAnnealing.cc (revision 2681) +++ /trunk/MagicSoft/Mars/mhist/MHSimulatedAnnealing.cc (revision 2681) @@ -0,0 +1,246 @@ +/* ======================================================================== *\ +! +! * +! * 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 "MBinning.h" + +#include +#include +#include + +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 2681) +++ /trunk/MagicSoft/Mars/mhist/MHSimulatedAnnealing.h (revision 2681) @@ -0,0 +1,67 @@ +#ifndef MARS_MHSimulatedannealing +#define MARS_MHSimulatedannealing +/////////////////////////////////////////////////////////////////////////////// +// +// MHSimulatedAnnealing +// +// Output container of MSimulatedAnnealing +//////////////////////////////////////////////////////////////////////////////////////////////////// +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TMatrix +#include "TMatrix.h" +#endif + +#ifndef ROOT_TVector +#include "TVector.h" +#endif + +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + +#ifndef ROOT_TH1 +#include "TH1.h" +#endif + +#ifndef ROOT_TObjarray +#include "TObjArray.h" +#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 2680) +++ /trunk/MagicSoft/Mars/mhist/Makefile (revision 2681) @@ -69,5 +69,7 @@ MHCalibrationBlindPixel.cc \ MHCalibrationPINDiode.cc \ - MHSupercuts.cc + MHCalibrationPINDiode.cc \ + MHSupercuts.cc \ + MHSimulatedAnnealing.cc # MHCurrents.cc \