From ee8ed05b81514b07a6f41fad52fd2c6a3496bdee Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 23 Mar 2016 17:07:31 +0000 Subject: [PATCH] various critical fixes and improvements to the new fit interface --- lib/FitInterface.cpp | 15 +++++++++++++-- lib/FitInterface.hpp | 6 +++++- lib/XYSampleData.cpp | 6 ++++-- lib/XYStatData.cpp | 10 ++++++++-- lib/XYStatData.hpp | 1 + 5 files changed, 31 insertions(+), 7 deletions(-) diff --git a/lib/FitInterface.cpp b/lib/FitInterface.cpp index ab39dfa..0bef6a1 100644 --- a/lib/FitInterface.cpp +++ b/lib/FitInterface.cpp @@ -407,7 +407,8 @@ void FitInterface::updateLayout(void) { if (initLayout_) { - Index size, ifit; + Index size, ifit; + vector v; layout.nXFitDim = 0; layout.nYFitDim = 0; @@ -416,6 +417,7 @@ void FitInterface::updateLayout(void) layout.totalYSize = 0; layout.xSize.clear(); layout.ySize.clear(); + layout.dataIndexSet.clear(); layout.xDim.clear(); layout.yDim.clear(); layout.xFitDim.clear(); @@ -479,6 +481,7 @@ void FitInterface::updateLayout(void) { if (p.second) { + layout.dataIndexSet.insert(p.first); layout.y[j].push_back(s); layout.yFit[j].push_back(layout.y[j].size() - 1); layout.data[j].push_back(p.first); @@ -495,7 +498,15 @@ void FitInterface::updateLayout(void) layout.totalSize = layout.totalXSize + layout.totalYSize; layout.nXFitDim = static_cast(layout.xSize.size()); layout.nYFitDim = static_cast(layout.ySize.size()); - initLayout_ = false; + for (Index k: layout.dataIndexSet) + { + v = dataCoord(k); + for (Index i = 0; i < getNXDim(); ++i) + { + layout.xIndFromData[k].push_back(indX(v[i], i)); + } + } + initLayout_ = false; } } diff --git a/lib/FitInterface.hpp b/lib/FitInterface.hpp index 53802c0..242262c 100644 --- a/lib/FitInterface.hpp +++ b/lib/FitInterface.hpp @@ -38,16 +38,20 @@ private: Index totalSize, totalXSize, totalYSize; // size of each X/Y dimension std::vector xSize, ySize; + // set of active data indices + std::set dataIndexSet; // lookup tables // xDim : x fit dim ifit -> x dim i // x : x fit point ifit,rfit -> x point r // xFitDim : x dim i -> x fit dim ifit (-1 if empty) // xFit : x point i,r -> x fit point rfit (-1 if empty) // data : y fit point jfit,sfit -> y point index k - // yFitFromData: y point indec k,j -> y fit point sfit (-1 if empty) + // yFitFromData: y point index k,j -> y fit point sfit (-1 if empty) + // xIndFromData: data index k -> index of coordinates of associated x std::vector xDim, yDim, xFitDim, yFitDim; std::vector> x, y, data, xFit, yFit; std::vector> yFitFromData; + std::map> xIndFromData; } Layout; public: // constructor diff --git a/lib/XYSampleData.cpp b/lib/XYSampleData.cpp index 1f09bfe..6ffd73f 100644 --- a/lib/XYSampleData.cpp +++ b/lib/XYSampleData.cpp @@ -207,7 +207,7 @@ void XYSampleData::setDataToSample(const Index s) { data_.x(r, i) = xData_[i][r][s]; } - for (Index j = 0; j < getNXDim(); ++j) + for (Index j = 0; j < getNYDim(); ++j) for (auto &p: yData_[j]) { data_.y(p.first, j) = p.second[s]; @@ -234,13 +234,15 @@ SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init, SampleFitResult result; FitResult sampleResult; + DVec initCopy = init; result.resize(nSample_); result.chi2_.resize(nSample_); FOR_STAT_ARRAY(result, s) { setDataToSample(s); - sampleResult = data_.fit(minimizer, init, v); + sampleResult = data_.fit(minimizer, initCopy, v); + initCopy = sampleResult.segment(0, initCopy.size()); result[s] = sampleResult; result.chi2_[s] = sampleResult.getChi2(); result.nDof_ = sampleResult.getNDof(); diff --git a/lib/XYStatData.cpp b/lib/XYStatData.cpp index af11f55..41ac711 100644 --- a/lib/XYStatData.cpp +++ b/lib/XYStatData.cpp @@ -250,6 +250,7 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init, void XYStatData::createXData(const std::string name __dumb, const Index nData) { xData_.push_back(DVec::Zero(nData)); + xBuf_.resize(xData_.size()); resizeVarMat(); } @@ -444,7 +445,7 @@ void XYStatData::updateChi2ModVec(const DVec p, { updateLayout(); - Index nPar = v[0]->getNPar(), a = 0, j, k; + Index nPar = v[0]->getNPar(), a = 0, j, k, ind; auto &par = p.segment(0, nPar), &xsi = p.segment(nPar, layout.totalXSize); updateXMap(); @@ -453,7 +454,12 @@ void XYStatData::updateChi2ModVec(const DVec p, { j = layout.yDim[jfit]; k = layout.data[jfit][sfit]; - chi2ModVec_(a) = (*v[j])(xMap_[k].data(), par.data()); + for (Index i = 0; i < getNXDim(); ++i) + { + ind = layout.xIndFromData[k][i] - layout.totalYSize; + xBuf_(i) = (ind >= 0) ? xsi(ind) : xMap_[k](i); + } + chi2ModVec_(a) = (*v[j])(xBuf_.data(), par.data()); a++; } chi2ModVec_.segment(a, layout.totalXSize) = xsi; diff --git a/lib/XYStatData.hpp b/lib/XYStatData.hpp index 68af292..22dd1ec 100644 --- a/lib/XYStatData.hpp +++ b/lib/XYStatData.hpp @@ -110,6 +110,7 @@ private: Mat xxVar_, yyVar_, xyVar_; DMat fitVar_, fitVarInv_; DVec chi2DataVec_, chi2ModVec_, chi2Vec_; + DVec xBuf_; bool initXMap_{true}; bool initChi2DataVec_{true}; };