portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2024-11-10 00:45:36 +00:00
Code Releases Activity

support for variable naming for function, model and data

Browse Source
This commit is contained in:
Antonin Portelli 2016-03-31 12:12:30 +01:00
parent 8ef69a4843
commit b3c2b17eef
13 changed files with 206 additions and 60 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
examples/exFit.cpp
View File

@ -24,9 +24,9 @@ int main(void)
DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
data.addXDim("x", nPoint1);
data.addXDim("off", nPoint2);
data.addYDim("y");
data.addXDim(nPoint1);
data.addXDim(nPoint2);
data.addYDim();
for (Index i1 = 0; i1 < nPoint1; ++i1)
{
xBuf[0] = i1*dx1;
@ -52,11 +52,11 @@ int main(void)
FitResult p;
MinuitMinimizer minimizer;
f.parName().setName(0, "m");
f.parName().setName(1, "A");
minimizer.setVerbosity(Minimizer::Verbosity::Normal);
p = data.fit(minimizer, init, f);
cout << "a= " << p(0) << " b= " << p(1) << endl;
cout << "chi^2/ndof= " << p.getChi2PerDof() << endl;
cout << "p-value= " << p.getPValue() <<endl;
p.print();
return EXIT_SUCCESS;
}

17
examples/exFitSample.cpp
View File

@ -24,9 +24,9 @@ int main(void)
DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
data.addXDim("x", nPoint1);
data.addXDim("off", nPoint2);
data.addYDim("y");
data.addXDim(nPoint1);
data.addXDim(nPoint2);
data.addYDim();
for (Index s = central; s < nSample; ++s)
{
for (Index i1 = 0; i1 < nPoint1; ++i1)
@ -47,16 +47,13 @@ int main(void)
// fit
DVec init = DVec::Constant(2, 0.1);
DMat err;
SampleFitResult p;
MinuitMinimizer minimizer;
p = data.fit(minimizer, init, f);
err = p.variance().cwiseSqrt();
cout << "a= " << p[central](0) << " +/- " << err(0) << endl;
cout << "b= " << p[central](1) << " +/- " << err(1) << endl;
cout << "chi^2/ndof= " << p.getChi2PerDof() << endl;
cout << "p-value= " << p.getPValue() << endl;
f.parName().setName(0, "m");
f.parName().setName(1, "A");
p = data.fit(minimizer, init, f);
p.print();
return EXIT_SUCCESS;
}

52
lib/FitInterface.cpp
View File

@ -26,6 +26,12 @@ using namespace Latan;
/******************************************************************************
* FitInterface implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
FitInterface::FitInterface(void)
: xName_("x")
, yName_("y")
{}
// copy object (not as a constructor to be accessed from derived class) ////////
void FitInterface::copyInterface(const FitInterface &d)
{
@ -33,7 +39,7 @@ void FitInterface::copyInterface(const FitInterface &d)
}
// add dimensions //////////////////////////////////////////////////////////////
void FitInterface::addXDim(const string name, const Index nData,
void FitInterface::addXDim(const Index nData, const string name,
const bool isExact)
{
if (getYSize() != 0)
@ -43,34 +49,38 @@ void FitInterface::addXDim(const string name, const Index nData,
}
else
{
xDimName_.push_back(name);
xSize_.push_back(nData);
xIsExact_.push_back(isExact);
xDimIndex_[name] = xDimName_.size();
maxDataIndex_ *= nData;
maxDataIndex_ *= nData;
createXData(name, nData);
scheduleLayoutInit();
if (!name.empty())
{
xName().setName(getNXDim(), name);
}
}
}
void FitInterface::addYDim(const string name)
{
yDimName_.push_back(name);
yDataIndex_.push_back(map<Index, bool>());
yDimIndex_[name] = yDimName_.size();
createYData(name);
scheduleLayoutInit();
if (!name.empty())
{
yName().setName(getNYDim(), name);
}
}
// size access /////////////////////////////////////////////////////////////////
// access //////////////////////////////////////////////////////////////////////
Index FitInterface::getNXDim(void) const
{
return xDimName_.size();
return xSize_.size();
}
Index FitInterface::getNYDim(void) const
{
return yDimName_.size();
return yDataIndex_.size();
}
Index FitInterface::getXSize(void) const
@ -175,6 +185,26 @@ Index FitInterface::getMaxDataIndex(void) const
return maxDataIndex_;
}
VarName & FitInterface::xName(void)
{
return xName_;
}
const VarName & FitInterface::xName(void) const
{
return xName_;
}
VarName & FitInterface::yName(void)
{
return yName_;
}
const VarName & FitInterface::yName(void) const
{
return yName_;
}
// Y dimension index helper ////////////////////////////////////////////////////
Index FitInterface::dataIndex(const vector<Index> &v) const
{
@ -554,7 +584,7 @@ ostream & Latan::operator<<(ostream &out, FitInterface &f)
out << "X dimensions: " << f.getNXDim() << endl;
for (Index i = 0; i < f.getNXDim(); ++i)
{
out << " * " << i << " \"" << f.xDimName_[i] << "\": ";
out << " * " << i << " \"" << f.xName().getName(i) << "\": ";
out << f.getXSize(i) << " value(s)";
if (f.xIsExact_[i])
{
@ -565,7 +595,7 @@ ostream & Latan::operator<<(ostream &out, FitInterface &f)
out << "Y dimensions: " << f.getNYDim() << endl;
for (Index j = 0; j < f.getNYDim(); ++j)
{
out << " * " << j << " \"" << f.yDimName_[j] << "\": ";
out << " * " << j << " \"" << f.yName().getName(j) << "\": ";
out << f.getYSize(j) << " value(s)" << endl;
for (auto &p: f.yDataIndex_[j])
{

37
lib/FitInterface.hpp
View File

@ -55,27 +55,31 @@ private:
} Layout;
public:
// constructor
FitInterface(void) = default;
FitInterface(void);
// destructor
virtual ~FitInterface(void) = default;
// copy object (not as a constructor to be accessed from derived class)
void copyInterface(const FitInterface &d);
// add dimensions
void addXDim(const std::string name, const Index nData,
void addXDim(const Index nData, const std::string name = "",
const bool isExact = false);
void addYDim(const std::string name);
// size access
Index getNXDim(void) const;
Index getNYDim(void) const;
Index getXSize(void) const;
Index getXSize(const Index i) const;
Index getYSize(void) const;
Index getYSize(const Index j) const;
Index getXFitSize(void) const;
Index getXFitSize(const Index i) const;
Index getYFitSize(void) const;
Index getYFitSize(const Index j) const;
Index getMaxDataIndex(void) const;
void addYDim(const std::string name = "");
// access
Index getNXDim(void) const;
Index getNYDim(void) const;
Index getXSize(void) const;
Index getXSize(const Index i) const;
Index getYSize(void) const;
Index getYSize(const Index j) const;
Index getXFitSize(void) const;
Index getXFitSize(const Index i) const;
Index getYFitSize(void) const;
Index getYFitSize(const Index j) const;
Index getMaxDataIndex(void) const;
VarName & xName(void);
const VarName & xName(void) const;
VarName & yName(void);
const VarName & yName(void) const;
// Y dimension index helper
template <typename... Ts>
Index dataIndex(const Ts... is) const;
@ -120,8 +124,7 @@ protected:
protected:
Layout layout;
private:
std::vector<std::string> xDimName_, yDimName_;
std::map<std::string, Index> xDimIndex_, yDimIndex_;
VarName xName_, yName_;
std::vector<Index> xSize_;
std::vector<bool> xIsExact_;
std::vector<std::map<Index, bool>> yDataIndex_;

11
lib/Function.cpp
View File

@ -29,6 +29,7 @@ using namespace Latan;
// constructor /////////////////////////////////////////////////////////////////
DoubleFunction::DoubleFunction(const vecFunc &f, const Index nArg)
: buffer_(new DVec)
, varName_("x")
{
setFunction(f, nArg);
}
@ -45,6 +46,16 @@ void DoubleFunction::setFunction(const vecFunc &f, const Index nArg)
f_ = f;
}
VarName & DoubleFunction::varName(void)
{
return varName_;
}
const VarName & DoubleFunction::varName(void) const
{
return varName_;
}
// error checking //////////////////////////////////////////////////////////////
void DoubleFunction::checkSize(const Index nPar) const
{

7
lib/Function.hpp
View File

@ -43,8 +43,10 @@ public:
// destructor
virtual ~DoubleFunction(void) = default;
// access
virtual Index getNArg(void) const;
void setFunction(const vecFunc &f, const Index nArg);
virtual Index getNArg(void) const;
void setFunction(const vecFunc &f, const Index nArg);
VarName & varName(void);
const VarName & varName(void) const;
// function call
double operator()(const double *arg) const;
double operator()(const DVec &arg) const;
@ -75,6 +77,7 @@ private:
void checkSize(const Index nPar) const;
private:
std::shared_ptr<DVec> buffer_{nullptr};
VarName varName_;
vecFunc f_;
};

10
lib/MinuitMinimizer.cpp
View File

@ -111,7 +111,7 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
min->SetFunction(minuitF);
for (Index i = 0; i < x.size(); ++i)
{
name = "x_" + strFrom(i);
name = f.varName().getName(i);
val = x(i);
step = (fabs(x(i)) != 0.) ? initErr*fabs(x(i)) : 1.;
if (hasHighLimit(i) and !hasLowLimit(i))
@ -139,7 +139,8 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
if (getVerbosity() >= Verbosity::Normal)
{
cout << "========== Minuit pre-minimization " << endl;
cout << "========== Minuit minimization, pass #1";
cout << " ==========" << endl;
}
min->SetStrategy(0);
min->Minimize();
@ -148,7 +149,8 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
n++;
if (getVerbosity() >= Verbosity::Normal)
{
cout << "========== Minuit minimization, try #" << n << endl;
cout << "========== Minuit minimization, pass #" << n + 1;
cout << " ==========" << endl;
}
min->SetStrategy(2);
min->Minimize();
@ -156,7 +158,7 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
} while (status and (n < maxTry));
if (getVerbosity() >= Verbosity::Normal)
{
cout << "==============================" << endl;
cout << "======================================" << endl;
}
switch (status)
{

22
lib/Model.cpp
View File

@ -31,6 +31,8 @@ using namespace Latan;
// constructor /////////////////////////////////////////////////////////////////
DoubleModel::DoubleModel(const vecFunc &f, const Index nArg, const Index nPar)
: size_(new ModelSize)
, varName_("x")
, parName_("p")
{
setFunction(f, nArg, nPar);
}
@ -54,6 +56,26 @@ void DoubleModel::setFunction(const vecFunc &f, const Index nArg,
f_ = f;
}
VarName & DoubleModel::varName(void)
{
return varName_;
}
const VarName & DoubleModel::varName(void) const
{
return varName_;
}
VarName & DoubleModel::parName(void)
{
return parName_;
}
const VarName & DoubleModel::parName(void) const
{
return parName_;
}
// error checking //////////////////////////////////////////////////////////////
void DoubleModel::checkSize(const Index nArg, const Index nPar) const
{

13
lib/Model.hpp
View File

@ -43,10 +43,14 @@ public:
// destructor
virtual ~DoubleModel(void) = default;
// access
virtual Index getNArg(void) const;
virtual Index getNPar(void) const;
void setFunction(const vecFunc &f, const Index nArg,
const Index nPar);
virtual Index getNArg(void) const;
virtual Index getNPar(void) const;
void setFunction(const vecFunc &f, const Index nArg,
const Index nPar);
VarName & varName(void);
const VarName & varName(void) const;
VarName & parName(void);
const VarName & parName(void) const;
// function call
double operator()(const DVec &data, const DVec &par) const;
double operator()(const std::vector<double> &data,
@ -61,6 +65,7 @@ private:
void checkSize(const Index nArg, const Index nPar) const;
private:
std::shared_ptr<ModelSize> size_;
VarName varName_, parName_;
vecFunc f_;
};

33
lib/XYSampleData.cpp
View File

@ -52,6 +52,11 @@ double SampleFitResult::getNDof(void) const
return static_cast<double>(nDof_);
}
Index SampleFitResult::getNPar(void) const
{
return nPar_;
}
double SampleFitResult::getPValue(const Index s) const
{
return Math::chi2PValue(getChi2(s), getNDof());
@ -86,6 +91,24 @@ FitResult SampleFitResult::getFitResult(const Index s) const
return fit;
}
// IO //////////////////////////////////////////////////////////////////////////
void SampleFitResult::print(const bool printXsi, ostream &out) const
{
char buf[256];
Index pMax = printXsi ? size() : nPar_;
DMat err = this->variance().cwiseSqrt();
sprintf(buf, "chi^2/dof= %.1f/%d= %.2f -- p-value= %.2e", getChi2(),
static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
out << buf << endl;
for (Index p = 0; p < pMax; ++p)
{
sprintf(buf, "%8s= % e +/- %e", parName_[p].c_str(),
(*this)[central](p), err(p));
out << buf << endl;
}
}
/******************************************************************************
* XYSampleData implementation *
******************************************************************************/
@ -238,6 +261,7 @@ SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init,
result.resize(nSample_);
result.chi2_.resize(nSample_);
result.model_.resize(v.size());
FOR_STAT_ARRAY(result, s)
{
setDataToSample(s);
@ -245,14 +269,17 @@ SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init,
initCopy = sampleResult.segment(0, initCopy.size());
result[s] = sampleResult;
result.chi2_[s] = sampleResult.getChi2();
result.nDof_ = sampleResult.getNDof();
result.model_.resize(v.size());
for (unsigned int j = 0; j < v.size(); ++j)
{
result.model_[j].resize(nSample_);
result.model_[j][s] = sampleResult.getModel(j);
}
}
result.nPar_ = sampleResult.getNPar();
result.nDof_ = sampleResult.getNDof();
result.parName_ = sampleResult.parName_;
return result;
}
@ -363,7 +390,7 @@ void XYSampleData::computeVarMat(void)
// create data /////////////////////////////////////////////////////////////////
void XYSampleData::createXData(const string name, const Index nData)
{
data_.addXDim(name, nData);
data_.addXDim(nData, name);
xData_.push_back(vector<DSample>(nData));
}

7
lib/XYSampleData.hpp
View File

@ -47,16 +47,21 @@ public:
double getChi2PerDof(const Index s = central) const;
DSample getChi2PerDof(const PlaceHolder ph) const;
double getNDof(void) const;
Index getNPar(void) const;
double getPValue(const Index s = central) const;
const DoubleFunction & getModel(const Index s = central,
const Index j = 0) const;
const DoubleFunctionSample & getModel(const PlaceHolder ph,
const Index j = 0) const;
FitResult getFitResult(const Index s = central) const;
// IO
void print(const bool printXsi = false,
std::ostream &out = std::cout) const;
private:
DSample chi2_;
double nDof_{0.};
Index nDof_{0}, nPar_{0};
std::vector<DoubleFunctionSample> model_;
std::vector<std::string> parName_;
};
/******************************************************************************

35
lib/XYStatData.cpp
View File

@ -43,6 +43,11 @@ double FitResult::getNDof(void) const
return static_cast<double>(nDof_);
}
Index FitResult::getNPar(void) const
{
return nPar_;
}
double FitResult::getPValue(void) const
{
return Math::chi2PValue(getChi2(), getNDof());;
@ -53,6 +58,22 @@ const DoubleFunction & FitResult::getModel(const Index j) const
return model_[static_cast<unsigned int>(j)];
}
// IO //////////////////////////////////////////////////////////////////////////
void FitResult::print(const bool printXsi, ostream &out) const
{
char buf[256];
Index pMax = printXsi ? size() : nPar_;
sprintf(buf, "chi^2/dof= %.1f/%d= %.2f -- p-value= %.2e", getChi2(),
static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
out << buf << endl;
for (Index p = 0; p < pMax; ++p)
{
sprintf(buf, "%8s= %e", parName_[p].c_str(), (*this)(p));
out << buf << endl;
}
}
/******************************************************************************
* XYStatData implementation *
******************************************************************************/
@ -226,6 +247,15 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
};
DoubleFunction chi2(chi2Func, totalNPar);
for (Index p = 0; p < nPar; ++p)
{
chi2.varName().setName(p, v[0]->parName().getName(p));
}
for (Index p = 0; p < totalNPar - nPar; ++p)
{
chi2.varName().setName(p + nPar, "xsi_" + strFrom(p));
}
// minimization
FitResult result;
DVec totalInit(totalNPar);
@ -236,12 +266,17 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
minimizer.setInit(totalInit);
result = minimizer(chi2);
result.chi2_ = chi2(result);
result.nPar_ = nPar;
result.nDof_ = layout.totalYSize - nPar;
result.model_.resize(v.size());
for (unsigned int j = 0; j < v.size(); ++j)
{
result.model_[j] = v[j]->fixPar(result);
}
for (Index p = 0; p < totalNPar; ++p)
{
result.parName_.push_back(chi2.varName().getName(p));
}
return result;
}

10
lib/XYStatData.hpp
View File

@ -33,6 +33,7 @@ BEGIN_LATAN_NAMESPACE
class FitResult: public DVec
{
friend class XYStatData;
friend class XYSampleData;
friend class SampleFitResult;
public:
// constructors
@ -44,12 +45,17 @@ public:
double getChi2(void) const;
double getChi2PerDof(void) const;
double getNDof(void) const;
Index getNPar(void) const;
double getPValue(void) const;
const DoubleFunction & getModel(const Index j = 0) const;
// IO
void print(const bool printXsi = false,
std::ostream &out = std::cout) const;
private:
double chi2_{0.0};
Index nDof_{0};
double chi2_{0.};
Index nDof_{0}, nPar_{0};
std::vector<DoubleFunction> model_;
std::vector<std::string> parName_;
};
/******************************************************************************