From 7b292980e5665d191dccfe3460504be8b91e373c Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Wed, 9 Mar 2016 16:28:32 +0100 Subject: [PATCH] added the class PFitterFcnDKS which eventually will handle the DKS/GPU binding. Currently it is without real functionality --- src/classes/Makefile.am | 4 +- src/classes/PFitter.cpp | 144 ++++++++++++++++++++++++++++------ src/classes/PFitterFcn.cpp | 8 +- src/classes/PFitterFcnDKS.cpp | 107 +++++++++++++++++++++++++ src/include/PFitter.h | 3 + src/include/PFitterFcn.h | 2 +- src/include/PFitterFcnDKS.h | 61 ++++++++++++++ 7 files changed, 299 insertions(+), 30 deletions(-) create mode 100644 src/classes/PFitterFcnDKS.cpp create mode 100644 src/include/PFitterFcnDKS.h diff --git a/src/classes/Makefile.am b/src/classes/Makefile.am index 3cd9792b..2afc9dfa 100644 --- a/src/classes/Makefile.am +++ b/src/classes/Makefile.am @@ -1,8 +1,9 @@ ## Process this file with automake to create Makefile.in h_sources = \ - ../include/PFitterFcn.h \ ../include/PFitter.h \ + ../include/PFitterFcn.h \ + ../include/PFitterFcnDKS.h \ ../include/PFourier.h \ ../include/PFourierCanvas.h \ ../include/PFunctionGrammar.h \ @@ -50,6 +51,7 @@ dict_h_sources_userFcn = \ cpp_sources = \ PFitter.cpp \ PFitterFcn.cpp \ + PFitterFcnDKS.cpp \ PFourier.cpp \ PFourierCanvas.cpp \ PFunction.cpp \ diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index e366bda1..c8e34adc 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -38,6 +38,7 @@ #include #include #include +#include #include #include @@ -82,6 +83,7 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bo fChisqOnly(chisq_only), fRunInfo(runInfo), fRunListCollection(runListCollection) { // initialize variables + fDKSReady = false; fIsScanOnly = true; fConverged = false; fUseChi2 = true; // chi^2 is the default @@ -93,6 +95,7 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bo // init class variables fFitterFcn = 0; + fFitterFcnDKS = 0; fFcnMin = 0; fScanAll = true; @@ -117,11 +120,26 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bo return; } - // create fit function object - fFitterFcn = new PFitterFcn(runListCollection, fUseChi2); - if (!fFitterFcn) { - fIsValid = false; - return; + // check if the theory function can already run on the GPU + string theo = fRunInfo->GetDKSTheoryString(); + if (!theo.compare("??")) { // theory not yet DKS ready + cout << endl << ">> PFitter::PFitter(): **INFO** theory not yet DKS/GPU ready. Will run on the CPU." << endl; + } else { + fDKSReady = true; + cout << endl << ">> PFitter::PFitter(): **INFO** theory DKS/GPU ready. Will run on the GPU." << endl; + } + + // create fit function object depending whether DKS/GPU can be used or not + if (fDKSReady) { // run on the GPU + fFitterFcnDKS = new PFitterFcnDKS(runListCollection, fUseChi2); + if (!fFitterFcnDKS) { + fIsValid = false; + } + } else { // run on the CPU + fFitterFcn = new PFitterFcn(runListCollection, fUseChi2); + if (!fFitterFcn) { + fIsValid = false; + } } } @@ -144,6 +162,11 @@ PFitter::~PFitter() fFcnMin = 0; } + if (fFitterFcnDKS) { + delete fFitterFcnDKS; + fFitterFcnDKS = 0; + } + if (fFitterFcn) { delete fFitterFcn; fFitterFcn = 0; @@ -172,13 +195,24 @@ Bool_t PFitter::DoFit() if (error[i] != 0.0) usedParams++; } - Int_t ndf = fFitterFcn->GetTotalNoOfFittedBins() - usedParams; - Double_t val = (*fFitterFcn)(param); + Int_t ndf = 0; + Double_t val = 0.0; + if (fDKSReady) { + ndf = fFitterFcnDKS->GetTotalNoOfFittedBins() - usedParams; + val = (*fFitterFcnDKS)(param); + } else { + ndf = fFitterFcn->GetTotalNoOfFittedBins() - usedParams; + val = (*fFitterFcn)(param); + } + if (fUseChi2) { // calculate expected chisq Double_t totalExpectedChisq = 0.0; std::vector expectedChisqPerHisto; - fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + if (fDKSReady) + fFitterFcnDKS->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + else + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); // calculate chisq per run std::vector chisqPerHisto; for (UInt_t i=0; iGetMsrRunList()->size(); i++) { @@ -191,14 +225,21 @@ Bool_t PFitter::DoFit() cout << endl << ">> expected chisq = " << totalExpectedChisq << ", NDF = " << ndf << ", expected chisq/NDF = " << totalExpectedChisq/ndf; UInt_t ndf_histo = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (fDKSReady) + ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + else + ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); if (ndf_histo > 0) cout << endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerHisto[i]/ndf_histo << "/" << expectedChisqPerHisto[i]/ndf_histo << ")"; } } else if (chisqPerHisto.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits UInt_t ndf_histo = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (fDKSReady) + ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + else + ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_histo > 0) cout << endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerHisto[i]/ndf_histo << ")"; } @@ -993,7 +1034,13 @@ Bool_t PFitter::ExecuteContours() return false; } - ROOT::Minuit2::MnContours contours((*fFitterFcn), *fFcnMin); + ROOT::Minuit2::FCNBase *fcn = 0; + if (fDKSReady) + fcn = fFitterFcnDKS; + else + fcn = fFitterFcn; + + ROOT::Minuit2::MnContours contours((*fcn), *fFcnMin); fScanData = contours(fScanParameter[0], fScanParameter[1], fScanNoPoints); @@ -1131,8 +1178,13 @@ Bool_t PFitter::ExecuteHesse() // call hesse Double_t start=0.0, end=0.0; + ROOT::Minuit2::FCNBase *fcn = 0; + if (fDKSReady) + fcn = fFitterFcnDKS; + else + fcn = fFitterFcn; start=MilliTime(); - ROOT::Minuit2::MnUserParameterState mnState = hesse((*fFitterFcn), fMnUserParams, maxfcn); + ROOT::Minuit2::MnUserParameterState mnState = hesse((*fcn), fMnUserParams, maxfcn); end=MilliTime(); cout << ">> PFitter::ExecuteMinimize(): execution time for Hesse = " << setprecision(3) << (end-start)/1.0e3 << " sec." << endl; TString str = TString::Format("Hesse: %.3f sec", (end-start)/1.0e3); @@ -1169,7 +1221,12 @@ Bool_t PFitter::ExecuteMigrad() // create migrad object // strategy is by default = 'default' - ROOT::Minuit2::MnMigrad migrad((*fFitterFcn), fMnUserParams, fStrategy); + ROOT::Minuit2::FCNBase *fcn = 0; + if (fDKSReady) + fcn = fFitterFcnDKS; + else + fcn = fFitterFcn; + ROOT::Minuit2::MnMigrad migrad((*fcn), fMnUserParams, fStrategy); // minimize // maxfcn is MINUIT2 Default maxfcn @@ -1211,7 +1268,11 @@ Bool_t PFitter::ExecuteMigrad() // handle statistics Double_t minVal = min.Fval(); - UInt_t ndf = fFitterFcn->GetTotalNoOfFittedBins(); + UInt_t ndf = 0.0; + if (fDKSReady) + ndf = fFitterFcnDKS->GetTotalNoOfFittedBins(); + else + ndf = fFitterFcn->GetTotalNoOfFittedBins(); // subtract number of varied parameters from total no of fitted bins -> ndf for (UInt_t i=0; iGetTotalNoOfFittedBins(); + UInt_t ndf = 0.0; + if (fDKSReady) + ndf = fFitterFcnDKS->GetTotalNoOfFittedBins(); + else + ndf = fFitterFcn->GetTotalNoOfFittedBins(); // subtract number of varied parameters from total no of fitted bins -> ndf for (UInt_t i=0; i> PFitter::ExecuteScan(): will call scan ..." << endl; - ROOT::Minuit2::MnScan scan((*fFitterFcn), fMnUserParams); + ROOT::Minuit2::FCNBase *fcn = 0; + if (fDKSReady) + fcn = fFitterFcnDKS; + else + fcn = fFitterFcn; + ROOT::Minuit2::MnScan scan((*fcn), fMnUserParams); if (fScanAll) { // not clear at the moment what to be done here // TO BE IMPLEMENTED @@ -1576,7 +1656,10 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) for (UInt_t i=0; iCalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + if (fDKSReady) + fFitterFcnDKS->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + else + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); // calculate chisq per run std::vector chisqPerHisto; @@ -1588,7 +1671,10 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) // get the ndf's of the histos UInt_t ndf_histo; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (fDKSReady) + ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + else + ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); ndfPerHisto.push_back(ndf_histo); } @@ -1603,7 +1689,10 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) } else if (chisqPerHisto.size() > 1) { // in case expected chisq is not applicable like for asymmetry fits UInt_t ndf_histo = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (fDKSReady) + ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + else + ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); ndfPerHisto.push_back(ndf_histo); } @@ -1901,7 +1990,12 @@ Bool_t PFitter::ExecuteSimplex() // create minimizer object // strategy is by default = 'default' - ROOT::Minuit2::MnSimplex simplex((*fFitterFcn), fMnUserParams, fStrategy); + ROOT::Minuit2::FCNBase *fcn = 0; + if (fDKSReady) + fcn = fFitterFcnDKS; + else + fcn = fFitterFcn; + ROOT::Minuit2::MnSimplex simplex((*fcn), fMnUserParams, fStrategy); // minimize // maxfcn is 10*MINUIT2 Default maxfcn @@ -1943,7 +2037,11 @@ Bool_t PFitter::ExecuteSimplex() // handle statistics Double_t minVal = min.Fval(); - UInt_t ndf = fFitterFcn->GetTotalNoOfFittedBins(); + UInt_t ndf = 0.0; + if (fDKSReady) + fFitterFcnDKS->GetTotalNoOfFittedBins(); + else + fFitterFcn->GetTotalNoOfFittedBins(); // subtract number of varied parameters from total no of fitted bins -> ndf for (UInt_t i=0; i +using namespace std; + +#include "PFitterFcnDKS.h" + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor. + * + * \param runList run list collection + * \param useChi2 if true, a chisq fit will be performed, otherwise a log max-likelihood fit will be carried out. + */ +PFitterFcnDKS::PFitterFcnDKS(PRunListCollection *runList, Bool_t useChi2) : + fRunListCollection(runList), + fUseChi2(useChi2) +{ + if (fUseChi2) + fUp = 1.0; + else + fUp = 0.5; +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Destructor + */ +PFitterFcnDKS::~PFitterFcnDKS() +{ + +} + +//-------------------------------------------------------------------------- +// operator() +//-------------------------------------------------------------------------- +/** + *

Minuit2 interface function call routine. This is the function which should be minimized. + * + * \param par a vector with all the parameters of the function + */ +Double_t PFitterFcnDKS::operator()(const std::vector& par) const +{ + Double_t value = 0.0; + + return value; +} + +//-------------------------------------------------------------------------- +// CalcExpectedChiSquare() +//-------------------------------------------------------------------------- +/** + *

Calculates the expected chisq, expected chisq per run, and chisq per run, if applicable. + * + * \param par + * \param totalExpectedChisq expected chisq for all run blocks + * \param expectedChisqPerRun expected chisq vector for all the run blocks + */ +void PFitterFcnDKS::CalcExpectedChiSquare(const std::vector &par, Double_t &totalExpectedChisq, std::vector &expectedChisqPerRun) +{ + // init expected chisq related variables + totalExpectedChisq = 0.0; + expectedChisqPerRun.clear(); + + // only do something for chisq + if (fUseChi2) { + Double_t value = 0.0; + + // single histo + for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { + value = fRunListCollection->GetSingleHistoChisqExpected(par, i); // calculate the expected chisq for single histo run block 'i' + expectedChisqPerRun.push_back(value); + totalExpectedChisq += value; + } + } +} diff --git a/src/include/PFitter.h b/src/include/PFitter.h index c9a09014..dad9675e 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -37,6 +37,7 @@ #include "PMsrHandler.h" #include "PRunListCollection.h" #include "PFitterFcn.h" +#include "PFitterFcnDKS.h" #define PMN_INTERACTIVE 0 #define PMN_CONTOURS 1 @@ -74,6 +75,7 @@ class PFitter Bool_t DoFit(); private: + Bool_t fDKSReady; ///< flag. true: fit via DKS/GPU. false: fit on CPU Bool_t fIsValid; ///< flag. true: the fit is valid. Bool_t fIsScanOnly; ///< flag. true: scan along some parameters (no fitting). Bool_t fConverged; ///< flag. true: the fit has converged. @@ -92,6 +94,7 @@ class PFitter PIntPairVector fCmdList; ///< command list, first=cmd, second=cmd line index PFitterFcn *fFitterFcn; ///< pointer to the fitter function object + PFitterFcnDKS *fFitterFcnDKS; ///< pointer to the DKS fitter function object ROOT::Minuit2::MnUserParameters fMnUserParams; ///< minuit2 input parameter list ROOT::Minuit2::FunctionMinimum *fFcnMin; ///< function minimum object diff --git a/src/include/PFitterFcn.h b/src/include/PFitterFcn.h index 94686e00..e21f44e0 100644 --- a/src/include/PFitterFcn.h +++ b/src/include/PFitterFcn.h @@ -37,7 +37,7 @@ #include "PRunListCollection.h" /** - *

This is the minuit2 interface function class porviding the function to be optimized (chisq or log max-likelihood). + *

This is the minuit2 interface function class providing the function to be optimized (chisq or log max-likelihood). */ class PFitterFcn : public ROOT::Minuit2::FCNBase { diff --git a/src/include/PFitterFcnDKS.h b/src/include/PFitterFcnDKS.h new file mode 100644 index 00000000..37bb4348 --- /dev/null +++ b/src/include/PFitterFcnDKS.h @@ -0,0 +1,61 @@ +/*************************************************************************** + + PFitterFcnDKS.h + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2016 by Andreas Suter * + * andreas.suter@psi.ch * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + +#ifndef _PFITTERFCNDKS_H_ +#define _PFITTERFCNDKS_H_ + +#include + +#include "Minuit2/FCNBase.h" + +#include "PRunListCollection.h" + +/** + *

This is the minuit2 interface function class providing the function to be optimized (chisq or log max-likelihood). + */ +class PFitterFcnDKS : public ROOT::Minuit2::FCNBase +{ + public: + PFitterFcnDKS(PRunListCollection *runList, Bool_t useChi2); + ~PFitterFcnDKS(); + + Double_t Up() const { return fUp; } + Double_t operator()(const std::vector &par) const; + + UInt_t GetTotalNoOfFittedBins() { return fRunListCollection->GetTotalNoOfBinsFitted(); } + UInt_t GetNoOfFittedBins(const UInt_t idx) { return fRunListCollection->GetNoOfBinsFitted(idx); } + void CalcExpectedChiSquare(const std::vector &par, Double_t &totalExpectedChisq, std::vector &expectedChisqPerRun); + + private: + Double_t fUp; ///< for chisq == 1.0, i.e. errors are 1 std. deviation errors. for log max-likelihood == 0.5, i.e. errors are 1 std. deviation errors (for details see the minuit2 user manual). + Bool_t fUseChi2; ///< true = chisq fit, false = log max-likelihood fit + PRunListCollection *fRunListCollection; ///< pre-processed data to be fitted +}; + +#endif // _PFITTERFCNDKS_H_