From 1f2150a42a90ff9849a6451421688eeaf7654276 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 16 Mar 2016 18:33:51 +0000 Subject: [PATCH] first implementation of the new XYSampleData (to be tested) --- lib/Makefile.am | 2 + lib/XYSampleData.cpp | 361 +++++++++++++++++++++++++++++++++++++++++++ lib/XYSampleData.hpp | 129 ++++++++++++++++ 3 files changed, 492 insertions(+) create mode 100644 lib/XYSampleData.cpp create mode 100644 lib/XYSampleData.hpp diff --git a/lib/Makefile.am b/lib/Makefile.am index ab6148c..c25aa8c 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -45,6 +45,7 @@ libLatAnalyze_la_SOURCES = \ RootFinder.cpp \ Solver.cpp \ TabFunction.cpp \ + XYSampleData.cpp \ XYStatData.cpp \ ../config.h libLatAnalyze_ladir = $(pkgincludedir) @@ -79,6 +80,7 @@ libLatAnalyze_la_HEADERS = \ TabFunction.hpp \ Solver.hpp \ StatArray.hpp \ + XYSampleData.hpp \ XYStatData.hpp if HAVE_MINUIT libLatAnalyze_la_SOURCES += MinuitMinimizer.cpp diff --git a/lib/XYSampleData.cpp b/lib/XYSampleData.cpp new file mode 100644 index 0000000..9a22d28 --- /dev/null +++ b/lib/XYSampleData.cpp @@ -0,0 +1,361 @@ +/* + * XYSampleData.cpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2016 Antonin Portelli + * + * LatAnalyze 3 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 3 of the License, or + * (at your option) any later version. + * + * LatAnalyze 3 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 LatAnalyze 3. If not, see . + */ + +#include +#include +#include + +using namespace std; +using namespace Latan; + +/****************************************************************************** + * SampleFitResult implementation * + ******************************************************************************/ +double SampleFitResult::getChi2(const Index s) const +{ + return chi2_[s]; +} + +const DSample & SampleFitResult::getChi2(const PlaceHolder ph __dumb) const +{ + return chi2_; +} + +double SampleFitResult::getChi2PerDof(const Index s) const +{ + return chi2_[s]/getNDof(); +} + +DSample SampleFitResult::getChi2PerDof(const PlaceHolder ph __dumb) const +{ + return chi2_/getNDof(); +} + +double SampleFitResult::getNDof(void) const +{ + return static_cast(nDof_); +} + +double SampleFitResult::getPValue(const Index s) const +{ + return Math::chi2PValue(getChi2(s), getNDof()); +} + +const DoubleFunction & SampleFitResult::getModel(const Index s, + const Index j) const +{ + return model_[static_cast(j)][s]; +} + +const DoubleFunctionSample & SampleFitResult::getModel( + const PlaceHolder ph __dumb, + const Index j) const +{ + return model_[static_cast(j)]; +} + +FitResult SampleFitResult::getFitResult(const Index s) const +{ + FitResult fit; + + fit = (*this)[s]; + fit.chi2_ = getChi2(); + fit.nDof_ = static_cast(getNDof()); + fit.model_.resize(model_.size()); + for (unsigned int k = 0; k < model_.size(); ++k) + { + fit.model_[k] = model_[k][s]; + } + + return fit; +} + +/****************************************************************************** + * XYSampleData implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +XYSampleData::XYSampleData(const Index nSample) +: nSample_(nSample) +{} + +// data access ///////////////////////////////////////////////////////////////// +DSample & XYSampleData::x(const Index r, const Index i) +{ + checkXIndex(r, i); + scheduleDataInit(); + scheduleComputeVarMat(); + + return xData_[i][r]; +} + +const DSample & XYSampleData::x(const Index r, const Index i) const +{ + checkXIndex(r, i); + + return xData_[i][r]; +} + +DSample & XYSampleData::y(const Index k, const Index j) +{ + checkYDim(j); + if (!pointExists(k, j)) + { + registerDataPoint(k, j); + } + scheduleDataInit(); + scheduleComputeVarMat(); + + return yData_[j][k]; +} + +const DSample & XYSampleData::y(const Index k, const Index j) const +{ + checkPoint(k, j); + + return yData_[j].at(k); +} + +const DMat & XYSampleData::getXXVar(const Index i1, const Index i2) +{ + checkXDim(i1); + checkXDim(i2); + computeVarMat(); + + return data_.getXXVar(i1, i2); +} + +const DMat & XYSampleData::getYYVar(const Index j1, const Index j2) +{ + checkYDim(j1); + checkYDim(j2); + computeVarMat(); + + return data_.getYYVar(j1, j2); +} + +const DMat & XYSampleData::getXYVar(const Index i, const Index j) +{ + checkXDim(i); + checkYDim(j); + computeVarMat(); + + return data_.getXYVar(i, j); +} + +DVec XYSampleData::getXError(const Index i) +{ + checkXDim(i); + computeVarMat(); + + return data_.getXError(i); +} + +DVec XYSampleData::getYError(const Index j) +{ + checkYDim(j); + computeVarMat(); + + return data_.getYError(j); +} + +// get total fit variance matrix and its pseudo-inverse //////////////////////// +const DMat & XYSampleData::getFitVarMat(void) +{ + computeVarMat(); + + return data_.getFitVarMat(); +} + +const DMat & XYSampleData::getFitVarMatPInv(void) +{ + computeVarMat(); + + return data_.getFitVarMatPInv(); +} + +// set data to a particular sample ///////////////////////////////////////////// +void XYSampleData::setDataToSample(const Index s) +{ + if (initData_ or (s != dataSample_)) + { + data_.copyInterface(*this); + for (Index i = 0; i < getNXDim(); ++i) + for (Index r = 0; r < getXSize(i); ++r) + { + data_.x(r, i) = xData_[i][r][s]; + } + for (Index j = 0; j < getNXDim(); ++j) + for (auto &p: yData_[j]) + { + data_.y(p.first, j) = p.second[s]; + } + dataSample_ = s; + initData_ = false; + } +} + +// get internal XYStatData ///////////////////////////////////////////////////// +const XYStatData & XYSampleData::getData(void) +{ + setDataToSample(central); + + return data_; +} + +// fit ///////////////////////////////////////////////////////////////////////// +SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init, + const std::vector &v) +{ + computeVarMat(); + + SampleFitResult result; + FitResult sampleResult; + + result.resize(nSample_); + result.chi2_.resize(nSample_); + FOR_STAT_ARRAY(result, s) + { + setDataToSample(s); + sampleResult = data_.fit(minimizer, init, v); + 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); + } + } + + return result; +} + +// schedule data initilization from samples //////////////////////////////////// +void XYSampleData::scheduleDataInit(void) +{ + initData_ = true; +} + +// variance matrix computation ///////////////////////////////////////////////// +void XYSampleData::scheduleComputeVarMat(void) +{ + computeVarMat_ = true; +} + +void XYSampleData::computeVarMat(void) +{ + if (computeVarMat_) + { + // initialize data if necessary + setDataToSample(central); + + // compute relevant sizes + Index size = 0, ySize = 0; + + for (Index j = 0; j < getNYDim(); ++j) + { + size += getYSize(j); + } + ySize = size; + for (Index i = 0; i < getNXDim(); ++i) + { + size += getXSize(i); + } + + // compute total matrix + DMatSample z(nSample_, size, 1); + DMat var; + Index a = 0; + + FOR_STAT_ARRAY(z, s) + { + for (Index j = 0; j < getNYDim(); ++j) + for (auto &p: yData_[j]) + { + z[s](a, 0) = p.second[s]; + a++; + } + for (Index i = 0; i < getNXDim(); ++i) + for (Index r = 0; r < getXSize(i); ++r) + { + z[s](a, 0) = xData_[i][r][s]; + a++; + } + } + var = z.varianceMatrix(); + + // assign blocks to data + Index a1, a2; + + a1 = ySize; + for (Index i1 = 0; i1 < getNXDim(); ++i1) + { + a2 = ySize; + for (Index i2 = 0; i2 < getNXDim(); ++i2) + { + data_.setXXVar(i1, i2, + var.block(a1, getXSize(i1), a2, getXSize(i2))); + a2 += getXSize(i2); + } + a1 += getXSize(i1); + } + a1 = 0; + for (Index j1 = 0; j1 < getNYDim(); ++j1) + { + a2 = 0; + for (Index j2 = 0; j2 < getNYDim(); ++j2) + { + data_.setYYVar(j1, j2, + var.block(a1, getYSize(j1), a2, getYSize(j2))); + a2 += getYSize(j2); + } + a1 += getYSize(j1); + } + a1 = ySize; + for (Index i = 0; i < getNXDim(); ++i) + { + a2 = 0; + for (Index j = 0; j < getNXDim(); ++j) + { + data_.setXYVar(i, j, + var.block(a1, getXSize(i), a2, getYSize(j))); + a2 += getYSize(j); + } + a1 += getXSize(i); + } + computeVarMat_ = false; + } + if (initVarMat()) + { + data_.scheduleFitVarMatInit(); + scheduleFitVarMatInit(false); + } +} + +// create data ///////////////////////////////////////////////////////////////// +void XYSampleData::createXData(const string name __dumb, const Index nData) +{ + xData_.push_back(vector(nData)); +} + +void XYSampleData::createYData(const string name __dumb) +{ + yData_.push_back(map()); +} diff --git a/lib/XYSampleData.hpp b/lib/XYSampleData.hpp new file mode 100644 index 0000000..b9b0025 --- /dev/null +++ b/lib/XYSampleData.hpp @@ -0,0 +1,129 @@ +/* + * XYSampleData.hpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2016 Antonin Portelli + * + * LatAnalyze 3 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 3 of the License, or + * (at your option) any later version. + * + * LatAnalyze 3 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 LatAnalyze 3. If not, see . + */ + +#ifndef Latan_XYSampleData_hpp_ +#define Latan_XYSampleData_hpp_ + +#include +#include +#include +#include +#include +#include + +BEGIN_LATAN_NAMESPACE + +/****************************************************************************** + * object for fit result * + ******************************************************************************/ +class SampleFitResult: public DMatSample +{ + friend class XYSampleData; +public: + // constructors + SampleFitResult(void) = default; + EIGEN_EXPR_CTOR(SampleFitResult, SampleFitResult, DMatSample, ArrayExpr) + // destructor + virtual ~SampleFitResult(void) = default; + // access + double getChi2(const Index s = central) const; + const DSample & getChi2(const PlaceHolder ph) const; + double getChi2PerDof(const Index s = central) const; + DSample getChi2PerDof(const PlaceHolder ph) const; + double getNDof(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; +private: + DSample chi2_; + double nDof_{0.}; + std::vector model_; +}; + +/****************************************************************************** + * XYSampleData * + ******************************************************************************/ +class XYSampleData: public FitInterface +{ +public: + // constructor + explicit XYSampleData(const Index nSample); + // destructor + virtual ~XYSampleData(void) = default; + // data access + DSample & x(const Index r, const Index i = 0); + const DSample & x(const Index r, const Index i = 0) const; + DSample & y(const Index k, const Index j = 0); + const DSample & y(const Index k, const Index j = 0) const; + const DMat & getXXVar(const Index i1, const Index i2); + const DMat & getYYVar(const Index j1, const Index j2); + const DMat & getXYVar(const Index i, const Index j); + DVec getXError(const Index i); + DVec getYError(const Index j); + // get total fit variance matrix and its pseudo-inverse + const DMat & getFitVarMat(void); + const DMat & getFitVarMatPInv(void); + // set data to a particular sample + void setDataToSample(const Index s); + // get internal XYStatData + const XYStatData & getData(void); + // fit + SampleFitResult fit(Minimizer &minimizer, const DVec &init, + const std::vector &v); + template + SampleFitResult fit(Minimizer &minimizer, const DVec &init, + const DoubleModel &model, const Mods... models); +private: + // schedule data initilization from samples + void scheduleDataInit(void); + // variance matrix computation + void scheduleComputeVarMat(void); + void computeVarMat(void); + // create data + virtual void createXData(const std::string name, const Index nData); + virtual void createYData(const std::string name); +private: + std::vector> yData_; + std::vector> xData_; + XYStatData data_; + Index nSample_, dataSample_{central}; + bool initData_{true}, computeVarMat_{true}; +}; + +/****************************************************************************** + * XYSampleData template implementation * + ******************************************************************************/ +template +SampleFitResult XYSampleData::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); +} + +END_LATAN_NAMESPACE + +#endif // Latan_XYSampleData_hpp_