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

fit: lookup table optimisation after profiling

Browse Source
This commit is contained in:
Antonin Portelli 2016-04-06 18:34:33 +01:00
parent 69ac2d5c82
commit ed8398145a
4 changed files with 41 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
lib/FitInterface.cpp
View File

@ -632,6 +632,7 @@ void FitInterface::updateLayout(void) const
l.totalSize = layout.totalXSize + layout.totalYSize; l.totalSize = layout.totalXSize + layout.totalYSize;
l.nXFitDim = static_cast<Index>(layout.xSize.size()); l.nXFitDim = static_cast<Index>(layout.xSize.size());
l.nYFitDim = static_cast<Index>(layout.ySize.size()); l.nYFitDim = static_cast<Index>(layout.ySize.size());
l.xIndFromData.resize(getMaxDataIndex());
for (Index k: layout.dataIndexSet) for (Index k: layout.dataIndexSet)
{ {
v = dataCoord(k); v = dataCoord(k);

3
lib/FitInterface.hpp
View File

@ -51,7 +51,8 @@ private:
std::vector<Index> xDim, yDim, xFitDim, yFitDim; std::vector<Index> xDim, yDim, xFitDim, yFitDim;
std::vector<std::vector<Index>> x, y, data, xFit, yFit; std::vector<std::vector<Index>> x, y, data, xFit, yFit;
std::vector<std::map<Index, Index>> yFitFromData; std::vector<std::map<Index, Index>> yFitFromData;
std::map<Index, std::vector<Index>> xIndFromData; // no map here for fit performance
std::vector<std::vector<Index>> xIndFromData;
} Layout; } Layout;
public: public:
// constructor // constructor

49
lib/XYStatData.cpp
View File

@ -267,19 +267,31 @@ FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
// get number of parameters // get number of parameters
Index nPar = v[0]->getNPar(); Index nPar = v[0]->getNPar();
Index nXDim = getNXDim();
Index totalNPar = nPar + layout.totalXSize; Index totalNPar = nPar + layout.totalXSize;
// chi^2 function // chi^2 functions
auto chi2Func = [this, totalNPar, &v](const double *x)->double auto corrChi2Func = [this, nPar, nXDim, totalNPar, &v](const double *x)->double
{ {
ConstMap<DVec> p(x, totalNPar); ConstMap<DVec> p(x, totalNPar);
updateChi2ModVec(p, v); updateChi2ModVec(p, v, nPar, nXDim);
chi2Vec_ = (chi2ModVec_ - chi2DataVec_); chi2Vec_ = (chi2ModVec_ - chi2DataVec_);
return (chi2Vec_.transpose()*fitVarInv_).dot(chi2Vec_); return chi2Vec_.dot(fitVarInv_*chi2Vec_);
}; };
DoubleFunction chi2(chi2Func, totalNPar); DoubleFunction corrChi2(corrChi2Func, totalNPar);
auto uncorrChi2Func = [this, nPar, nXDim, totalNPar, &v](const double *x)->double
{
ConstMap<DVec> p(x, totalNPar);
updateChi2ModVec(p, v, nPar, nXDim);
chi2Vec_ = (chi2ModVec_ - chi2DataVec_);
return chi2Vec_.dot(chi2Vec_.cwiseQuotient(fitVar_.diagonal()));
};
DoubleFunction uncorrChi2(uncorrChi2Func, totalNPar);
DoubleFunction &chi2 = hasCorrelations() ? corrChi2 : uncorrChi2;
for (Index p = 0; p < nPar; ++p) for (Index p = 0; p < nPar; ++p)
{ {
@ -523,6 +535,7 @@ void XYStatData::updateXMap(void) const
XYStatData * modThis = const_cast<XYStatData *>(this); XYStatData * modThis = const_cast<XYStatData *>(this);
modThis->xMap_.clear(); modThis->xMap_.clear();
modThis->xMap_.resize(getMaxDataIndex());
for (auto k: getDataIndexSet()) for (auto k: getDataIndexSet())
{ {
modThis->xMap_[k] = DVec(getNXDim()); modThis->xMap_[k] = DVec(getNXDim());
@ -563,26 +576,32 @@ void XYStatData::updateChi2DataVec(void)
} }
} }
// WARNING: updateChi2ModVec is heavily called by fit
void XYStatData::updateChi2ModVec(const DVec p, void XYStatData::updateChi2ModVec(const DVec p,
const vector<const DoubleModel *> &v) const vector<const DoubleModel *> &v,
const Index nPar, const Index nXDim)
{ {
updateLayout(); updateLayout();
updateXMap();
Index nPar = v[0]->getNPar(), a = 0, j, k, ind; Index a = 0, j, k, ind;
auto &par = p.segment(0, nPar), &xsi = p.segment(nPar, layout.totalXSize); auto &par = p.segment(0, nPar), &xsi = p.segment(nPar, layout.totalXSize);
for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit) for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
{ {
j = layout.yDim[jfit]; j = layout.yDim[jfit];
k = layout.data[jfit][sfit]; for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
for (Index i = 0; i < getNXDim(); ++i)
{ {
ind = layout.xIndFromData[k][i] - layout.totalYSize;
xBuf_(i) = (ind >= 0) ? xsi(ind) : x(k)(i); k = layout.data[jfit][sfit];
for (Index i = 0; i < nXDim; ++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_(a) = (*v[j])(xBuf_.data(), par.data());
a++;
} }
chi2ModVec_.segment(a, layout.totalXSize) = xsi; chi2ModVec_.segment(a, layout.totalXSize) = xsi;
} }

6
lib/XYStatData.hpp
View File

@ -121,11 +121,13 @@ private:
// buffer chi^2 vectors // buffer chi^2 vectors
void updateChi2DataVec(void); void updateChi2DataVec(void);
void updateChi2ModVec(const DVec p, void updateChi2ModVec(const DVec p,
const std::vector<const DoubleModel *> &v); const std::vector<const DoubleModel *> &v,
const Index nPar, const Index nXDim);
private: private:
std::vector<std::map<Index, double>> yData_; std::vector<std::map<Index, double>> yData_;
// no map here for fit performance
std::vector<DVec> xData_; std::vector<DVec> xData_;
std::map<Index, DVec> xMap_; std::vector<DVec> xMap_;
Mat<DMat> xxVar_, yyVar_, xyVar_; Mat<DMat> xxVar_, yyVar_, xyVar_;
DMat fitVar_, fitVarInv_; DMat fitVar_, fitVarInv_;
DVec chi2DataVec_, chi2ModVec_, chi2Vec_; DVec chi2DataVec_, chi2ModVec_, chi2Vec_;