LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

added the class PFitterFcnDKS which eventually will handle the DKS/GPU binding. Currently it is without real functionality

Browse Source
This commit is contained in:
suter_a 2016-03-09 16:28:32 +01:00
parent d9782c55c8
commit 7b292980e5
7 changed files with 299 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/classes/Makefile.am
View File

@ -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 \

144
src/classes/PFitter.cpp
View File

@ -38,6 +38,7 @@
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <limits>
#include <cmath>
@ -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<Double_t> 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<Double_t> chisqPerHisto;
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->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; i<expectedChisqPerHisto.size(); i++) {
ndf_histo = fFitterFcn->GetNoOfFittedBins(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; i<chisqPerHisto.size(); i++) {
ndf_histo = fFitterFcn->GetNoOfFittedBins(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; i<fParams.size(); i++) {
if ((min.UserState().Error(i) != 0.0) && !fMnUserParams.Parameters().at(i).IsFixed())
@ -1244,7 +1305,12 @@ Bool_t PFitter::ExecuteMinimize()
// create minimizer object
// strategy is by default = 'default'
ROOT::Minuit2::MnMinimize minimize((*fFitterFcn), fMnUserParams, fStrategy);
ROOT::Minuit2::FCNBase *fcn = 0;
if (fDKSReady)
fcn = fFitterFcnDKS;
else
fcn = fFitterFcn;
ROOT::Minuit2::MnMinimize minimize((*fcn), fMnUserParams, fStrategy);
// minimize
// maxfcn is MINUIT2 Default maxfcn
@ -1287,7 +1353,11 @@ Bool_t PFitter::ExecuteMinimize()
// 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; i<fParams.size(); i++) {
if ((min.UserState().Error(i) != 0.0) && !fMnUserParams.Parameters().at(i).IsFixed())
@ -1334,8 +1404,13 @@ Bool_t PFitter::ExecuteMinos()
// make minos analysis
Double_t start=0.0, end=0.0;
ROOT::Minuit2::FCNBase *fcn = 0;
if (fDKSReady)
fcn = fFitterFcnDKS;
else
fcn = fFitterFcn;
start=MilliTime();
ROOT::Minuit2::MnMinos minos((*fFitterFcn), (*fFcnMin));
ROOT::Minuit2::MnMinos minos((*fcn), (*fFcnMin));
for (UInt_t i=0; i<fParams.size(); i++) {
// only try to call minos if the parameter is not fixed!!
@ -1523,7 +1598,12 @@ Bool_t PFitter::ExecuteScan()
{
cout << ">> 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; i<fParams.size(); i++)
param.push_back(fParams[i].fValue);
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto);
if (fDKSReady)
fFitterFcnDKS->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto);
else
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto);
// calculate chisq per run
std::vector<Double_t> 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; i<expectedChisqPerHisto.size(); i++) {
ndf_histo = fFitterFcn->GetNoOfFittedBins(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; i<chisqPerHisto.size(); i++) {
ndf_histo = fFitterFcn->GetNoOfFittedBins(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<fParams.size(); i++) {
if ((min.UserState().Error(i) != 0.0) && !fMnUserParams.Parameters().at(i).IsFixed())

8
src/classes/PFitterFcn.cpp
View File

@ -41,16 +41,14 @@ using namespace std;
* \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.
*/
PFitterFcn::PFitterFcn(PRunListCollection *runList, Bool_t useChi2)
PFitterFcn::PFitterFcn(PRunListCollection *runList, Bool_t useChi2) :
fRunListCollection(runList),
fUseChi2(useChi2)
{
fUseChi2 = useChi2;
if (fUseChi2)
fUp = 1.0;
else
fUp = 0.5;
fRunListCollection = runList;
}
//--------------------------------------------------------------------------

107
src/classes/PFitterFcnDKS.cpp Normal file
View File

@ -0,0 +1,107 @@
/***************************************************************************
PFitterFcnDKS.cpp
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. *
***************************************************************************/
#include <iostream>
using namespace std;
#include "PFitterFcnDKS.h"
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p>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
//--------------------------------------------------------------------------
/**
* <p>Destructor
*/
PFitterFcnDKS::~PFitterFcnDKS()
{
}
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
* <p>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<Double_t>& par) const
{
Double_t value = 0.0;
return value;
}
//--------------------------------------------------------------------------
// CalcExpectedChiSquare()
//--------------------------------------------------------------------------
/**
* <p>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<Double_t> &par, Double_t &totalExpectedChisq, std::vector<Double_t> &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; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
value = fRunListCollection->GetSingleHistoChisqExpected(par, i); // calculate the expected chisq for single histo run block 'i'
expectedChisqPerRun.push_back(value);
totalExpectedChisq += value;
}
}
}

3
src/include/PFitter.h
View File

@ -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

2
src/include/PFitterFcn.h
View File

@ -37,7 +37,7 @@
#include "PRunListCollection.h"
/**
* <p>This is the minuit2 interface function class porviding the function to be optimized (chisq or log max-likelihood).
* <p>This is the minuit2 interface function class providing the function to be optimized (chisq or log max-likelihood).
*/
class PFitterFcn : public ROOT::Minuit2::FCNBase
{

61
src/include/PFitterFcnDKS.h Normal file
View File

@ -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 <vector>
#include "Minuit2/FCNBase.h"
#include "PRunListCollection.h"
/**
* <p>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<Double_t> &par) const;
UInt_t GetTotalNoOfFittedBins() { return fRunListCollection->GetTotalNoOfBinsFitted(); }
UInt_t GetNoOfFittedBins(const UInt_t idx) { return fRunListCollection->GetNoOfBinsFitted(idx); }
void CalcExpectedChiSquare(const std::vector<Double_t> &par, Double_t &totalExpectedChisq, std::vector<Double_t> &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_