diff --git a/examples/exFit.cpp b/examples/exFit.cpp index e8db6e7..be5b04f 100644 --- a/examples/exFit.cpp +++ b/examples/exFit.cpp @@ -1,46 +1,51 @@ #include +#include +#include +#include +#include #include 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(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() <> yDataIndex_; std::set> xxCorr_, yyCorr_, xyCorr_; Index maxDataIndex_{1}; - bool initLayout_{true}; + bool initLayout_{true}, initVarMat_{true}; }; std::ostream & operator<<(std::ostream &out, FitInterface &f); diff --git a/lib/XYStatData.cpp b/lib/XYStatData.cpp index f870bdc..7d5dfd9 100644 --- a/lib/XYStatData.cpp +++ b/lib/XYStatData.cpp @@ -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 &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; } diff --git a/lib/XYStatData.hpp b/lib/XYStatData.hpp index b8680a6..03b0c8d 100644 --- a/lib/XYStatData.hpp +++ b/lib/XYStatData.hpp @@ -83,6 +83,9 @@ public: // fit FitResult fit(Minimizer &minimizer, const DVec &init, const std::vector &v); + template + 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 xxVar_, yyVar_, xyVar_; DMat fitVar_, fitVarInv_; DVec chi2DataVec_, chi2ModVec_, chi2Vec_; - bool initVarMat_{true}; bool initXMap_{true}; bool initChi2DataVec_{true}; }; +/****************************************************************************** + * XYStatData template implementation * + ******************************************************************************/ +template +FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init, + const DoubleModel &model, const Ts... models) +{ + static_assert(static_or::value...>::value, + "model arguments are not compatible with DoubleModel"); + + std::vector modelVector{&model, &models...}; + + return fit(minimizer, init, modelVector); +} + /****************************************************************************** * error check macros * ******************************************************************************/