portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2025-11-04 08:04:32 +00:00
Code Releases Activity

XYStatData: various fixes and improvement, fit is now working

Browse Source
This commit is contained in:
Antonin Portelli
2016-03-16 13:05:56 +00:00
parent 9f67ead605
commit 9bc9288dfb
5 changed files with 87 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files

73
examples/exFit.cpp
View File

@@ -1,46 +1,51 @@
#include <iostream>
#include <cmath>
#include <LatAnalyze/CompiledModel.hpp>
#include <LatAnalyze/MinuitMinimizer.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <LatAnalyze/XYStatData.hpp>
using namespace std;
using namespace Latan;
const Index nPoint = 30;
const double xErr = .01, yErr = .1;
const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast<double>(nPoint);
int main(void)
{
XYStatData f;
// generate fake data
XYStatData data;
RandGen rg;
double x_k, y_k;
DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0]);}, 1, 2);
f.addYDim("q1");
f.addYDim("q2");
f.addXDim("x1", 6);
f.addXDim("x2", 5);
f.addXDim("x3", 5);
f.y(f.dataIndex(0,0,0), 0) = 2;
f.y(f.dataIndex(1,1,1), 0) = 4;
f.y(f.dataIndex(2,2,2), 0) = 5;
f.y(f.dataIndex(2,3,3), 0) = 10;
f.y(f.dataIndex(0,0,0), 1) = 1;
f.y(f.dataIndex(1,1,1), 1) = 2;
f.y(f.dataIndex(2,2,3), 1) = 4;
f.fitPoint(false, f.dataIndex(2,2,2), 0);
f.fitPoint(false, f.dataIndex(1,1,1), 1);
f.assumeXXCorrelated(true, 0, 0, 1, 0);
f.assumeXXCorrelated(true, 0, 1, 1, 1);
f.assumeXXCorrelated(true, 0, 2, 1, 2);
f.assumeXXCorrelated(true, 0, 0, 1, 2);
f.assumeXXCorrelated(true, 3, 2, 4, 2);
f.assumeYYCorrelated(true, f.dataIndex(0,0,0), 0, f.dataIndex(2,3,3), 0);
f.assumeYYCorrelated(true, f.dataIndex(0,0,0), 1, f.dataIndex(2,2,3), 1);
f.assumeXYCorrelated(true, 0, 0, f.dataIndex(1,1,1), 0);
f.assumeXExact(true, 0);
f.assumeXExact(true, 1);
f.assumeXExact(true, 2);
cout << f << endl;
f.setXXVar(0, 0, DMat::Identity(6, 6));
f.setXXVar(0, 2, DMat::Identity(6, 5));
f.setXXVar(1, 1, DMat::Identity(5, 5));
f.setXXVar(2, 2, DMat::Identity(5, 5));
DEBUG_MAT(f.makeCorrFilter());
DEBUG_MAT(f.getFitVar());
DEBUG_MAT(f.getFitVar().cwiseProduct(f.makeCorrFilter()));
data.addXDim("x", nPoint);
data.addYDim("y");
for (Index k = 0; k < nPoint; ++k)
{
x_k = k*dx + rg.gaussian(0.0, xErr);
y_k = f(&x_k, exactPar) + rg.gaussian(0.0, yErr);
printf("% 8e % 8e % 8e % 8e\n", x_k, xErr, y_k, yErr);
data.x(k) = x_k;
data.y(k) = y_k;
}
cout << endl;
data.setXError(0, DVec::Constant(nPoint, xErr));
data.setYError(0, DVec::Constant(nPoint, yErr));
cout << data << endl;
// fit
DVec init = DVec::Constant(2, 0.5);
FitResult p;
MinuitMinimizer minimizer;
minimizer.setVerbosity(MinuitMinimizer::Verbosity::Debug);
p = data.fit(minimizer, init, f);
cout << "a= " << p(0) << " b= " << p(1);
cout << " chi^2/ndof= " << p.getChi2PerDof();
cout << " p-value= " << p.getPValue() <<endl;
return EXIT_SUCCESS;
}

17
lib/FitInterface.cpp
View File

@@ -253,7 +253,7 @@ void FitInterface::assumeXXCorrelated(const bool isCorr, const Index r1,
{
addCorr(xxCorr_, isCorr, c);
}
scheduleLayoutInit();
scheduleFitVarMatInit();
}
void FitInterface::assumeYYCorrelated(const bool isCorr, const Index k1,
@@ -268,7 +268,7 @@ void FitInterface::assumeYYCorrelated(const bool isCorr, const Index k1,
{
addCorr(yyCorr_, isCorr, c);
}
scheduleLayoutInit();
scheduleFitVarMatInit();
}
void FitInterface::assumeXYCorrelated(const bool isCorr, const Index r,
@@ -280,7 +280,7 @@ void FitInterface::assumeXYCorrelated(const bool isCorr, const Index r,
checkXIndex(r, i);
checkPoint(k, j);
addCorr(xyCorr_, isCorr, c);
scheduleLayoutInit();
scheduleFitVarMatInit();
}
// tests ///////////////////////////////////////////////////////////////////////
@@ -387,6 +387,17 @@ void FitInterface::registerDataPoint(const Index k, const Index j)
void FitInterface::scheduleLayoutInit(void)
{
initLayout_ = true;
scheduleFitVarMatInit();
}
void FitInterface::scheduleFitVarMatInit(const bool init)
{
initVarMat_ = init;
}
bool FitInterface::initVarMat(void)
{
return initVarMat_;
}
void FitInterface::updateLayout(void)

4
lib/FitInterface.hpp
View File

@@ -105,6 +105,8 @@ protected:
virtual void createYData(void) = 0;
// global layout management
void scheduleLayoutInit(void);
void scheduleFitVarMatInit(const bool init = true);
bool initVarMat(void);
void updateLayout(void);
Index indX(const Index r, const Index i) const;
Index indY(const Index k, const Index j) const;
@@ -118,7 +120,7 @@ private:
std::vector<std::map<Index, bool>> yDataIndex_;
std::set<std::array<Index, 4>> xxCorr_, yyCorr_, xyCorr_;
Index maxDataIndex_{1};
bool initLayout_{true};
bool initLayout_{true}, initVarMat_{true};
};
std::ostream & operator<<(std::ostream &out, FitInterface &f);

24
lib/XYStatData.cpp
View File

@@ -205,6 +205,11 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
// check model consistency
checkModelVec(v);
// buffering
updateLayout();
updateFitVarMat();
updateChi2DataVec();
// get number of parameters
Index nPar = v[0]->getNPar();
Index totalNPar = nPar + layout.totalXSize;
@@ -225,8 +230,6 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
FitResult result;
DVec totalInit(totalNPar);
updateFitVarMat();
updateChi2DataVec();
totalInit.segment(0, nPar) = init;
totalInit.segment(nPar, layout.totalXSize) =
chi2DataVec_.segment(layout.totalYSize, layout.totalXSize);
@@ -280,11 +283,6 @@ void XYStatData::resizeVarMat(void)
}
// schedule buffer computation /////////////////////////////////////////////////
void XYStatData::scheduleFitVarMatInit(void)
{
initVarMat_ = true;
}
void XYStatData::scheduleXMapInit(void)
{
initXMap_ = true;
@@ -298,7 +296,7 @@ void XYStatData::scheduleChi2DataVecInit(void)
// buffer total fit variance matrix ////////////////////////////////////////////
void XYStatData::updateFitVarMat(void)
{
if (initVarMat_)
if (initVarMat())
{
updateLayout();
@@ -378,8 +376,9 @@ void XYStatData::updateFitVarMat(void)
chi2DataVec_.resize(layout.totalSize);
chi2ModVec_.resize(layout.totalSize);
chi2Vec_.resize(layout.totalSize);
fitVarInv_ = fitVar_.pInverse();
initVarMat_ = false;
fitVar_ = fitVar_.cwiseProduct(makeCorrFilter());
fitVarInv_ = fitVar_.pInverse();
scheduleFitVarMatInit(false);
}
}
@@ -419,6 +418,7 @@ void XYStatData::updateChi2DataVec(void)
{
Index a = 0, j, k, i, r;
updateLayout();
for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
{
@@ -442,6 +442,8 @@ void XYStatData::updateChi2DataVec(void)
void XYStatData::updateChi2ModVec(const DVec p,
const vector<const DoubleModel *> &v)
{
updateLayout();
Index nPar = v[0]->getNPar(), a = 0, j, k;
auto &par = p.segment(0, nPar), &xsi = p.segment(nPar, layout.totalXSize);
@@ -454,7 +456,7 @@ void XYStatData::updateChi2ModVec(const DVec p,
chi2ModVec_(a) = (*v[j])(xMap_[k].data(), par.data());
a++;
}
chi2ModVec_.segment(a, layout.totalSize) = xsi;
chi2ModVec_.segment(a, layout.totalXSize) = xsi;
}

20
lib/XYStatData.hpp
View File

@@ -83,6 +83,9 @@ public:
// fit
FitResult fit(Minimizer &minimizer, const DVec &init,
const std::vector<const DoubleModel *> &v);
template <typename... Ts>
FitResult fit(Minimizer &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models);
protected:
// create data
virtual void createXData(const Index nData);
@@ -90,7 +93,6 @@ protected:
void resizeVarMat(void);
private:
// schedule buffer computation
void scheduleFitVarMatInit(void);
void scheduleXMapInit(void);
void scheduleChi2DataVecInit(void);
// buffer total fit variance matrix
@@ -108,11 +110,25 @@ private:
Mat<DMat> xxVar_, yyVar_, xyVar_;
DMat fitVar_, fitVarInv_;
DVec chi2DataVec_, chi2ModVec_, chi2Vec_;
bool initVarMat_{true};
bool initXMap_{true};
bool initChi2DataVec_{true};
};
/******************************************************************************
* XYStatData template implementation *
******************************************************************************/
template <typename... Ts>
FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models)
{
static_assert(static_or<std::is_assignable<DoubleModel &, Ts>::value...>::value,
"model arguments are not compatible with DoubleModel");
std::vector<const DoubleModel *> modelVector{&model, &models...};
return fit(minimizer, init, modelVector);
}
/******************************************************************************
* error check macros *
******************************************************************************/