mirror of
https://github.com/aportelli/LatAnalyze.git
synced 2024-11-10 00:45:36 +00:00
new fit: partial residuals and plots
This commit is contained in:
parent
f82b20dc73
commit
07bacf0e4b
@ -4,6 +4,7 @@
|
||||
#include <LatAnalyze/Io.hpp>
|
||||
#include <LatAnalyze/MinuitMinimizer.hpp>
|
||||
#include <LatAnalyze/NloptMinimizer.hpp>
|
||||
#include <LatAnalyze/Plot.hpp>
|
||||
#include <LatAnalyze/RandGen.hpp>
|
||||
#include <LatAnalyze/XYStatData.hpp>
|
||||
|
||||
@ -11,7 +12,7 @@ using namespace std;
|
||||
using namespace Latan;
|
||||
|
||||
const Index nPoint1 = 10, nPoint2 = 10;
|
||||
const double xErr = .1, yErr = .1;
|
||||
const double xErr = .1, yErr = .3;
|
||||
const double exactPar[2] = {0.5,5.};
|
||||
const double dx1 = 10.0/static_cast<double>(nPoint1);
|
||||
const double dx2 = 5.0/static_cast<double>(nPoint2);
|
||||
@ -25,6 +26,7 @@ int main(void)
|
||||
DoubleModel f([](const double *x, const double *p)
|
||||
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
|
||||
|
||||
cout << "-- generating fake data..." << endl;
|
||||
data.addXDim(nPoint1);
|
||||
data.addXDim(nPoint2);
|
||||
data.addYDim();
|
||||
@ -34,10 +36,10 @@ int main(void)
|
||||
data.x(i1, 0) = rg.gaussian(xBuf[0], xErr);
|
||||
for (Index i2 = 0; i2 < nPoint2; ++i2)
|
||||
{
|
||||
xBuf[1] = i2*dx2;
|
||||
data.x(i2, 1) = xBuf[1];
|
||||
data.y(data.dataIndex(i1, i2)) = rg.gaussian(f(xBuf, exactPar),
|
||||
yErr);
|
||||
xBuf[1] = i2*dx2;
|
||||
data.x(i2, 1) = xBuf[1];
|
||||
data.y(data.dataIndex(i1, i2), 0) = rg.gaussian(f(xBuf, exactPar),
|
||||
yErr);
|
||||
}
|
||||
}
|
||||
data.setXError(0, DVec::Constant(data.getXSize(0), xErr));
|
||||
@ -52,7 +54,7 @@ int main(void)
|
||||
vector<Minimizer *> min{&globalMin, &localMin};
|
||||
|
||||
globalMin.setVerbosity(Minimizer::Verbosity::Normal);
|
||||
globalMin.setPrecision(0.01);
|
||||
globalMin.setPrecision(0.1);
|
||||
globalMin.setMaxIteration(10000);
|
||||
globalMin.useLowLimit(0);
|
||||
globalMin.setLowLimit(0, 0.);
|
||||
@ -65,10 +67,26 @@ int main(void)
|
||||
localMin.setVerbosity(Minimizer::Verbosity::Normal);
|
||||
|
||||
// fit
|
||||
cout << "-- fit..." << endl;
|
||||
f.parName().setName(0, "m");
|
||||
f.parName().setName(1, "A");
|
||||
p = data.fit(min, init, f);
|
||||
p.print();
|
||||
|
||||
// plot
|
||||
Plot plot;
|
||||
DVec ref(2);
|
||||
XYStatData res;
|
||||
|
||||
cout << "-- generating plots..." << endl;
|
||||
ref(1) = 0.;
|
||||
res = data.getPartialResiduals(p, ref, 0);
|
||||
plot << PlotRange(Axis::x, 0., 10.);
|
||||
plot << Color("rgb 'blue'");
|
||||
plot << PlotFunction(p.getModel().bind(0, ref), 0., 10.);
|
||||
plot << Color("rgb 'red'");
|
||||
plot << PlotData(res);
|
||||
plot.display();
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
@ -3,6 +3,7 @@
|
||||
#include <LatAnalyze/CompiledModel.hpp>
|
||||
#include <LatAnalyze/MinuitMinimizer.hpp>
|
||||
#include <LatAnalyze/NloptMinimizer.hpp>
|
||||
#include <LatAnalyze/Plot.hpp>
|
||||
#include <LatAnalyze/RandGen.hpp>
|
||||
#include <LatAnalyze/XYSampleData.hpp>
|
||||
|
||||
@ -11,7 +12,7 @@ using namespace Latan;
|
||||
|
||||
const Index nPoint1 = 10, nPoint2 = 10;
|
||||
const Index nSample = 1000;
|
||||
const double xErr = .1, yErr = .1;
|
||||
const double xErr = .1, yErr = .3;
|
||||
const double exactPar[2] = {0.5,5.};
|
||||
const double dx1 = 10.0/static_cast<double>(nPoint1);
|
||||
const double dx2 = 5.0/static_cast<double>(nPoint2);
|
||||
@ -25,6 +26,7 @@ int main(void)
|
||||
DoubleModel f([](const double *x, const double *p)
|
||||
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
|
||||
|
||||
cout << "-- generating fake data..." << endl;
|
||||
data.addXDim(nPoint1);
|
||||
data.addXDim(nPoint2);
|
||||
data.addYDim();
|
||||
@ -36,9 +38,9 @@ int main(void)
|
||||
data.x(i1, 0)[s] = rg.gaussian(xBuf[0], xErr);
|
||||
for (Index i2 = 0; i2 < nPoint2; ++i2)
|
||||
{
|
||||
xBuf[1] = i2*dx2;
|
||||
data.x(i2, 1)[s] = xBuf[1];
|
||||
data.y(data.dataIndex(i1, i2))[s] =
|
||||
xBuf[1] = i2*dx2;
|
||||
data.x(i2, 1)[s] = xBuf[1];
|
||||
data.y(data.dataIndex(i1, i2), 0)[s] =
|
||||
rg.gaussian(f(xBuf, exactPar), yErr);
|
||||
}
|
||||
}
|
||||
@ -52,7 +54,7 @@ int main(void)
|
||||
MinuitMinimizer localMin;
|
||||
vector<Minimizer *> min{&globalMin, &localMin};
|
||||
|
||||
globalMin.setPrecision(0.01);
|
||||
globalMin.setPrecision(0.1);
|
||||
globalMin.setMaxIteration(10000);
|
||||
globalMin.useLowLimit(0);
|
||||
globalMin.setLowLimit(0, 0.);
|
||||
@ -64,10 +66,28 @@ int main(void)
|
||||
globalMin.setHighLimit(1, 20.);
|
||||
|
||||
// fit
|
||||
cout << "-- fit..." << endl;
|
||||
f.parName().setName(0, "m");
|
||||
f.parName().setName(1, "A");
|
||||
p = data.fit(min, init, f);
|
||||
p.print();
|
||||
|
||||
// plot
|
||||
Plot plot;
|
||||
DVec ref(2);
|
||||
XYStatData res;
|
||||
|
||||
cout << "-- generating plots..." << endl;
|
||||
ref(1) = 0.;
|
||||
res = data.getPartialResiduals(p, ref, 0).getData();
|
||||
plot << PlotRange(Axis::x, 0., 10.);
|
||||
plot << Color("rgb 'blue'");
|
||||
plot << PlotPredBand(p.getModel(_).bind(0, ref), 0., 10.);
|
||||
plot << Color("rgb 'blue'");
|
||||
plot << PlotFunction(p.getModel().bind(0, ref), 0., 10.);
|
||||
plot << Color("rgb 'red'");
|
||||
plot << PlotData(res);
|
||||
plot.display();
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
||||
@ -54,6 +54,7 @@ void FitInterface::addXDim(const Index nData, const string name,
|
||||
maxDataIndex_ *= nData;
|
||||
createXData(name, nData);
|
||||
scheduleLayoutInit();
|
||||
scheduleDataCoordInit();
|
||||
if (!name.empty())
|
||||
{
|
||||
xName().setName(getNXDim(), name);
|
||||
@ -121,7 +122,7 @@ Index FitInterface::getYSize(const Index j) const
|
||||
return static_cast<Index>(yDataIndex_[j].size());
|
||||
}
|
||||
|
||||
Index FitInterface::getXFitSize(void) const
|
||||
Index FitInterface::getXFitSize(void)
|
||||
{
|
||||
Index size = 0;
|
||||
|
||||
@ -133,7 +134,7 @@ Index FitInterface::getXFitSize(void) const
|
||||
return size;
|
||||
}
|
||||
|
||||
Index FitInterface::getXFitSize(const Index i) const
|
||||
Index FitInterface::getXFitSize(const Index i)
|
||||
{
|
||||
set<Index> fitCoord;
|
||||
vector<Index> v;
|
||||
@ -185,6 +186,11 @@ Index FitInterface::getMaxDataIndex(void) const
|
||||
return maxDataIndex_;
|
||||
}
|
||||
|
||||
const set<Index> & FitInterface::getDataIndexSet(void) const
|
||||
{
|
||||
return dataIndexSet_;
|
||||
}
|
||||
|
||||
VarName & FitInterface::xName(void)
|
||||
{
|
||||
return xName_;
|
||||
@ -221,22 +227,13 @@ Index FitInterface::dataIndex(const vector<Index> &v) const
|
||||
return k;
|
||||
}
|
||||
|
||||
vector<Index> FitInterface::dataCoord(const Index k) const
|
||||
const vector<Index> & FitInterface::dataCoord(const Index k)
|
||||
{
|
||||
vector<Index> v(getNXDim());
|
||||
Index buf, dimProd;
|
||||
|
||||
checkDataIndex(k);
|
||||
buf = k;
|
||||
dimProd = 1;
|
||||
for (Index d = getNXDim() - 1; d >= 0; --d)
|
||||
{
|
||||
v[d] = (buf/dimProd)%xSize_[d];
|
||||
buf -= dimProd*v[d];
|
||||
dimProd *= xSize_[d];
|
||||
}
|
||||
|
||||
return v;
|
||||
updateDataCoord();
|
||||
|
||||
return dataCoord_.at(k);
|
||||
}
|
||||
|
||||
// enable fit points ///////////////////////////////////////////////////////////
|
||||
@ -335,7 +332,12 @@ bool FitInterface::pointExists(const Index k, const Index j) const
|
||||
return !(yDataIndex_[j].find(k) == yDataIndex_[j].end());
|
||||
}
|
||||
|
||||
bool FitInterface::isXUsed(const Index r, const Index i, const bool inFit) const
|
||||
bool FitInterface::isXExact(const Index i) const
|
||||
{
|
||||
return xIsExact_[i];
|
||||
}
|
||||
|
||||
bool FitInterface::isXUsed(const Index r, const Index i, const bool inFit)
|
||||
{
|
||||
vector<Index> v;
|
||||
|
||||
@ -346,7 +348,7 @@ bool FitInterface::isXUsed(const Index r, const Index i, const bool inFit) const
|
||||
{
|
||||
if (p.second or !inFit)
|
||||
{
|
||||
v = dataCoord(p.first);
|
||||
v = dataCoord(p.first);
|
||||
if (v[i] == r)
|
||||
{
|
||||
return true;
|
||||
@ -418,9 +420,29 @@ void FitInterface::registerDataPoint(const Index k, const Index j)
|
||||
{
|
||||
checkYDim(j);
|
||||
yDataIndex_[j][k] = true;
|
||||
dataIndexSet_.insert(k);
|
||||
scheduleLayoutInit();
|
||||
}
|
||||
|
||||
// coordinate buffering ////////////////////////////////////////////////////////
|
||||
void FitInterface::scheduleDataCoordInit(void)
|
||||
{
|
||||
initDataCoord_ = true;
|
||||
}
|
||||
|
||||
void FitInterface::updateDataCoord(void)
|
||||
{
|
||||
if (initDataCoord_)
|
||||
{
|
||||
dataCoord_.clear();
|
||||
for (auto k: getDataIndexSet())
|
||||
{
|
||||
dataCoord_[k] = rowMajToCoord(k);
|
||||
}
|
||||
initDataCoord_ = false;
|
||||
}
|
||||
}
|
||||
|
||||
// global layout management ////////////////////////////////////////////////////
|
||||
void FitInterface::scheduleLayoutInit(void)
|
||||
{
|
||||
@ -578,6 +600,25 @@ Index FitInterface::indY(const Index k, const Index j) const
|
||||
return ind;
|
||||
}
|
||||
|
||||
// function to convert an row-major index into coordinates /////////////////////
|
||||
vector<Index> FitInterface::rowMajToCoord(const Index k) const
|
||||
{
|
||||
vector<Index> v(getNXDim());
|
||||
Index buf, dimProd;
|
||||
|
||||
checkDataIndex(k);
|
||||
buf = k;
|
||||
dimProd = 1;
|
||||
for (Index d = getNXDim() - 1; d >= 0; --d)
|
||||
{
|
||||
v[d] = (buf/dimProd)%xSize_[d];
|
||||
buf -= dimProd*v[d];
|
||||
dimProd *= xSize_[d];
|
||||
}
|
||||
|
||||
return v;
|
||||
}
|
||||
|
||||
// IO //////////////////////////////////////////////////////////////////////////
|
||||
ostream & Latan::operator<<(ostream &out, FitInterface &f)
|
||||
{
|
||||
|
||||
@ -65,26 +65,27 @@ public:
|
||||
const bool isExact = false);
|
||||
void addYDim(const std::string name = "");
|
||||
// access
|
||||
Index getNXDim(void) const;
|
||||
Index getNYDim(void) const;
|
||||
Index getXSize(void) const;
|
||||
Index getXSize(const Index i) const;
|
||||
Index getYSize(void) const;
|
||||
Index getYSize(const Index j) const;
|
||||
Index getXFitSize(void) const;
|
||||
Index getXFitSize(const Index i) const;
|
||||
Index getYFitSize(void) const;
|
||||
Index getYFitSize(const Index j) const;
|
||||
Index getMaxDataIndex(void) const;
|
||||
VarName & xName(void);
|
||||
const VarName & xName(void) const;
|
||||
VarName & yName(void);
|
||||
const VarName & yName(void) const;
|
||||
Index getNXDim(void) const;
|
||||
Index getNYDim(void) const;
|
||||
Index getXSize(void) const;
|
||||
Index getXSize(const Index i) const;
|
||||
Index getYSize(void) const;
|
||||
Index getYSize(const Index j) const;
|
||||
Index getXFitSize(void);
|
||||
Index getXFitSize(const Index i);
|
||||
Index getYFitSize(void) const;
|
||||
Index getYFitSize(const Index j) const;
|
||||
Index getMaxDataIndex(void) const;
|
||||
const std::set<Index> & getDataIndexSet(void) const;
|
||||
VarName & xName(void);
|
||||
const VarName & xName(void) const;
|
||||
VarName & yName(void);
|
||||
const VarName & yName(void) const;
|
||||
// Y dimension index helper
|
||||
template <typename... Ts>
|
||||
Index dataIndex(const Ts... is) const;
|
||||
Index dataIndex(const std::vector<Index> &v) const;
|
||||
std::vector<Index> dataCoord(const Index k) const;
|
||||
Index dataIndex(const Ts... is) const;
|
||||
Index dataIndex(const std::vector<Index> &v) const;
|
||||
const std::vector<Index> & dataCoord(const Index k);
|
||||
// enable fit points
|
||||
void fitPoint(const bool isFitPoint, const Index k, const Index j = 0);
|
||||
// variance interface
|
||||
@ -98,7 +99,8 @@ public:
|
||||
// tests
|
||||
bool pointExists(const Index k) const;
|
||||
bool pointExists(const Index k, const Index j) const;
|
||||
bool isXUsed(const Index r, const Index i, const bool inFit = true) const;
|
||||
bool isXExact(const Index i) const;
|
||||
bool isXUsed(const Index r, const Index i, const bool inFit = true);
|
||||
bool isFitPoint(const Index k, const Index j) const;
|
||||
// make correlation filter for fit variance matrix
|
||||
DMat makeCorrFilter(void);
|
||||
@ -115,22 +117,31 @@ protected:
|
||||
// abstract methods to create data containers
|
||||
virtual void createXData(const std::string name, const Index nData) = 0;
|
||||
virtual void createYData(const std::string name) = 0;
|
||||
// coordinate buffering
|
||||
void scheduleDataCoordInit(void);
|
||||
void updateDataCoord(void);
|
||||
// global layout management
|
||||
void scheduleLayoutInit(void);
|
||||
bool initVarMat(void);
|
||||
void updateLayout(void);
|
||||
Index indX(const Index r, const Index i) const;
|
||||
Index indY(const Index k, const Index j) const;
|
||||
private:
|
||||
// function to convert an row-major index into coordinates
|
||||
std::vector<Index> rowMajToCoord(const Index k) const;
|
||||
protected:
|
||||
Layout layout;
|
||||
private:
|
||||
VarName xName_, yName_;
|
||||
std::vector<Index> xSize_;
|
||||
std::vector<bool> xIsExact_;
|
||||
std::vector<std::map<Index, bool>> yDataIndex_;
|
||||
std::set<std::array<Index, 4>> xxCorr_, yyCorr_, xyCorr_;
|
||||
Index maxDataIndex_{1};
|
||||
bool initLayout_{true}, initVarMat_{true};
|
||||
VarName xName_, yName_;
|
||||
std::vector<Index> xSize_;
|
||||
std::vector<bool> xIsExact_;
|
||||
std::map<Index, std::vector<Index>> dataCoord_;
|
||||
std::set<Index> dataIndexSet_;
|
||||
std::vector<std::map<Index, bool>> yDataIndex_;
|
||||
std::set<std::array<Index, 4>> xxCorr_, yyCorr_, xyCorr_;
|
||||
Index maxDataIndex_{1};
|
||||
bool initLayout_{true}, initVarMat_{true};
|
||||
bool initDataCoord_{true};
|
||||
};
|
||||
|
||||
std::ostream & operator<<(std::ostream &out, FitInterface &f);
|
||||
|
||||
23
lib/Plot.cpp
23
lib/Plot.cpp
@ -167,20 +167,15 @@ PlotData::PlotData(const DMatSample &x, const DVec &y)
|
||||
setCommand("'" + tmpFileName + "' u 1:3:2 w xerr");
|
||||
}
|
||||
|
||||
//PlotData::PlotData(const XYStatData &data, const Index i, const Index j)
|
||||
//{
|
||||
// DMat d(data.getNData(), 4);
|
||||
// string usingCmd, tmpFileName;
|
||||
//
|
||||
// d.col(0) = data.x(i);
|
||||
// d.col(2) = data.y(j);
|
||||
// d.col(1) = data.xxVar(i, i).diagonal().array().sqrt();
|
||||
// d.col(3) = data.yyVar(j, j).diagonal().array().sqrt();
|
||||
// usingCmd = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:3:2:4 w xyerr";
|
||||
// tmpFileName = dumpToTmpFile(d);
|
||||
// pushTmpFile(tmpFileName);
|
||||
// setCommand("'" + tmpFileName + "' " + usingCmd);
|
||||
//}
|
||||
PlotData::PlotData(XYStatData &data, const Index i, const Index j)
|
||||
{
|
||||
string usingCmd, tmpFileName;
|
||||
|
||||
usingCmd = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:3:2:4 w xyerr";
|
||||
tmpFileName = dumpToTmpFile(data.getTable(i, j));
|
||||
pushTmpFile(tmpFileName);
|
||||
setCommand("'" + tmpFileName + "' " + usingCmd);
|
||||
}
|
||||
|
||||
// PlotHLine constructor ///////////////////////////////////////////////////////
|
||||
PlotHLine::PlotHLine(const double y)
|
||||
|
||||
@ -25,7 +25,7 @@
|
||||
#include <LatAnalyze/Mat.hpp>
|
||||
#include <LatAnalyze/MatSample.hpp>
|
||||
#include <LatAnalyze/Histogram.hpp>
|
||||
//#include <LatAnalyze/XYStatData.hpp>
|
||||
#include <LatAnalyze/XYStatData.hpp>
|
||||
#include <sstream>
|
||||
#include <vector>
|
||||
|
||||
@ -94,7 +94,7 @@ public:
|
||||
PlotData(const DMatSample &x, const DMatSample &y);
|
||||
PlotData(const DVec &x, const DMatSample &y);
|
||||
PlotData(const DMatSample &x, const DVec &y);
|
||||
//PlotData(const XYStatData &data, const Index i = 0, const Index j = 0);
|
||||
PlotData(XYStatData &data, const Index i = 0, const Index j = 0);
|
||||
// destructor
|
||||
virtual ~PlotData(void) = default;
|
||||
};
|
||||
|
||||
@ -138,6 +138,15 @@ const DSample & XYSampleData::x(const Index r, const Index i) const
|
||||
return xData_[i][r];
|
||||
}
|
||||
|
||||
const DMatSample & XYSampleData::x(const Index k)
|
||||
{
|
||||
checkDataIndex(k);
|
||||
|
||||
updateXMap();
|
||||
|
||||
return xMap_.at(k);
|
||||
}
|
||||
|
||||
DSample & XYSampleData::y(const Index k, const Index j)
|
||||
{
|
||||
checkYDim(j);
|
||||
@ -299,6 +308,74 @@ SampleFitResult XYSampleData::fit(Minimizer &minimizer,
|
||||
return fit(mv, init, v);
|
||||
}
|
||||
|
||||
// residuals ///////////////////////////////////////////////////////////////////
|
||||
XYSampleData XYSampleData::getResiduals(const SampleFitResult &fit)
|
||||
{
|
||||
XYSampleData res(*this);
|
||||
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
const DoubleFunctionSample &f = fit.getModel(_, j);
|
||||
|
||||
for (auto &p: yData_[j])
|
||||
{
|
||||
res.y(p.first, j) -= f(x(p.first));
|
||||
}
|
||||
}
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
XYSampleData XYSampleData::getPartialResiduals(const SampleFitResult &fit,
|
||||
const DVec &ref, const Index i)
|
||||
{
|
||||
XYSampleData res(*this);
|
||||
DMatSample buf(nSample_);
|
||||
|
||||
buf.fill(ref);
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
const DoubleFunctionSample &f = fit.getModel(_, j);
|
||||
|
||||
for (auto &p: yData_[j])
|
||||
{
|
||||
FOR_STAT_ARRAY(buf, s)
|
||||
{
|
||||
buf[s](i) = x(p.first)[s](i);
|
||||
}
|
||||
res.y(p.first, j) -= f(x(p.first)) - f(buf);
|
||||
}
|
||||
}
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
// buffer list of x vectors ////////////////////////////////////////////////////
|
||||
void XYSampleData::scheduleXMapInit(void)
|
||||
{
|
||||
initXMap_ = true;
|
||||
}
|
||||
|
||||
void XYSampleData::updateXMap(void)
|
||||
{
|
||||
if (initXMap_)
|
||||
{
|
||||
for (Index s = central; s < nSample_; ++s)
|
||||
{
|
||||
setDataToSample(s);
|
||||
for (auto k: getDataIndexSet())
|
||||
{
|
||||
if (s == central)
|
||||
{
|
||||
xMap_[k].resize(nSample_);
|
||||
}
|
||||
xMap_[k][s] = data_.x(k);
|
||||
}
|
||||
}
|
||||
initXMap_ = false;
|
||||
}
|
||||
}
|
||||
|
||||
// schedule data initilization from samples ////////////////////////////////////
|
||||
void XYSampleData::scheduleDataInit(void)
|
||||
{
|
||||
|
||||
@ -75,15 +75,16 @@ public:
|
||||
// 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);
|
||||
DSample & x(const Index r, const Index i);
|
||||
const DSample & x(const Index r, const Index i) const;
|
||||
const DMatSample & x(const Index k);
|
||||
DSample & y(const Index k, const Index j);
|
||||
const DSample & y(const Index k, const Index j) 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);
|
||||
@ -102,7 +103,14 @@ public:
|
||||
template <typename... Ts>
|
||||
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
|
||||
const DoubleModel &model, const Ts... models);
|
||||
// residuals
|
||||
XYSampleData getResiduals(const SampleFitResult &fit);
|
||||
XYSampleData getPartialResiduals(const SampleFitResult &fit, const DVec &x,
|
||||
const Index i);
|
||||
private:
|
||||
// buffer list of x vectors
|
||||
void scheduleXMapInit(void);
|
||||
void updateXMap(void);
|
||||
// schedule data initilization from samples
|
||||
void scheduleDataInit(void);
|
||||
// variance matrix computation
|
||||
@ -114,9 +122,11 @@ private:
|
||||
private:
|
||||
std::vector<std::map<Index, DSample>> yData_;
|
||||
std::vector<std::vector<DSample>> xData_;
|
||||
std::map<Index, DMatSample> xMap_;
|
||||
XYStatData data_;
|
||||
Index nSample_, dataSample_{central};
|
||||
bool initData_{true}, computeVarMat_{true};
|
||||
bool initXMap_{true};
|
||||
};
|
||||
|
||||
/******************************************************************************
|
||||
|
||||
@ -96,6 +96,15 @@ const double & XYStatData::x(const Index r, const Index i) const
|
||||
return xData_[i](r);
|
||||
}
|
||||
|
||||
const DVec & XYStatData::x(const Index k)
|
||||
{
|
||||
checkDataIndex(k);
|
||||
|
||||
updateXMap();
|
||||
|
||||
return xMap_.at(k);
|
||||
}
|
||||
|
||||
double & XYStatData::y(const Index k, const Index j)
|
||||
{
|
||||
checkYDim(j);
|
||||
@ -206,6 +215,29 @@ DVec XYStatData::getYError(const Index j) const
|
||||
return yyVar_(j, j).diagonal().cwiseSqrt();
|
||||
}
|
||||
|
||||
DMat XYStatData::getTable(const Index i, const Index j)
|
||||
{
|
||||
checkXDim(i);
|
||||
checkYDim(j);
|
||||
|
||||
DMat table(getYSize(j), 4);
|
||||
Index row = 0;
|
||||
|
||||
for (auto &p: yData_[j])
|
||||
{
|
||||
Index k = p.first;
|
||||
Index r = dataCoord(k)[i];
|
||||
|
||||
table(row, 0) = x(k)(i);
|
||||
table(row, 2) = p.second;
|
||||
table(row, 1) = xxVar_(i, i).diagonal().cwiseSqrt()(r);
|
||||
table(row, 3) = yyVar_(j, j).diagonal().cwiseSqrt()(row);
|
||||
row++;
|
||||
}
|
||||
|
||||
return table;
|
||||
}
|
||||
|
||||
// get total fit variance matrix ///////////////////////////////////////////////
|
||||
const DMat & XYStatData::getFitVarMat(void)
|
||||
{
|
||||
@ -308,6 +340,44 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
|
||||
return fit(mv, init, v);
|
||||
}
|
||||
|
||||
// residuals ///////////////////////////////////////////////////////////////////
|
||||
XYStatData XYStatData::getResiduals(const FitResult &fit)
|
||||
{
|
||||
XYStatData res(*this);
|
||||
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
const DoubleFunction &f = fit.getModel(j);
|
||||
|
||||
for (auto &p: yData_[j])
|
||||
{
|
||||
res.y(p.first, j) -= f(x(p.first));
|
||||
}
|
||||
}
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
XYStatData XYStatData::getPartialResiduals(const FitResult &fit,
|
||||
const DVec &ref, const Index i)
|
||||
{
|
||||
XYStatData res(*this);
|
||||
DVec buf(ref);
|
||||
|
||||
for (Index j = 0; j < res.getNYDim(); ++j)
|
||||
{
|
||||
const DoubleFunction &f = fit.getModel(j);
|
||||
|
||||
for (auto &p: yData_[j])
|
||||
{
|
||||
buf(i) = x(p.first)(i);
|
||||
res.y(p.first, j) -= f(x(p.first)) - f(buf);
|
||||
}
|
||||
}
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
// create data /////////////////////////////////////////////////////////////////
|
||||
void XYStatData::createXData(const std::string name __dumb, const Index nData)
|
||||
{
|
||||
@ -450,24 +520,13 @@ void XYStatData::updateXMap(void)
|
||||
{
|
||||
if (initXMap_)
|
||||
{
|
||||
vector<Index> v;
|
||||
set<Index> indSet;
|
||||
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
for (auto &p: yData_[j])
|
||||
{
|
||||
indSet.insert(p.first);
|
||||
}
|
||||
}
|
||||
xMap_.clear();
|
||||
for (auto k: indSet)
|
||||
for (auto k: getDataIndexSet())
|
||||
{
|
||||
xMap_[k] = DVec(getNXDim());
|
||||
v = dataCoord(k);
|
||||
xMap_[k] = DVec(getNXDim());
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
{
|
||||
xMap_[k](i) = xData_[i](v[i]);
|
||||
xMap_[k](i) = xData_[i](dataCoord(k)[i]);
|
||||
}
|
||||
}
|
||||
initXMap_ = false;
|
||||
@ -510,7 +569,6 @@ void XYStatData::updateChi2ModVec(const DVec p,
|
||||
Index nPar = v[0]->getNPar(), a = 0, j, k, ind;
|
||||
auto &par = p.segment(0, nPar), &xsi = p.segment(nPar, layout.totalXSize);
|
||||
|
||||
updateXMap();
|
||||
for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
|
||||
for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
|
||||
{
|
||||
@ -519,7 +577,7 @@ void XYStatData::updateChi2ModVec(const DVec p,
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
{
|
||||
ind = layout.xIndFromData[k][i] - layout.totalYSize;
|
||||
xBuf_(i) = (ind >= 0) ? xsi(ind) : xMap_[k](i);
|
||||
xBuf_(i) = (ind >= 0) ? xsi(ind) : x(k)(i);
|
||||
}
|
||||
chi2ModVec_(a) = (*v[j])(xBuf_.data(), par.data());
|
||||
a++;
|
||||
|
||||
@ -69,10 +69,11 @@ public:
|
||||
// destructor
|
||||
virtual ~XYStatData(void) = default;
|
||||
// data access
|
||||
double & x(const Index r, const Index i = 0);
|
||||
const double & x(const Index r, const Index i = 0) const;
|
||||
double & y(const Index k, const Index j = 0);
|
||||
const double & y(const Index k, const Index j = 0) const;
|
||||
double & x(const Index r, const Index);
|
||||
const double & x(const Index r, const Index) const;
|
||||
const DVec & x(const Index k);
|
||||
double & y(const Index k, const Index);
|
||||
const double & y(const Index k, const Index) const;
|
||||
void setXXVar(const Index i1, const Index i2, const DMat &m);
|
||||
void setYYVar(const Index j1, const Index j2, const DMat &m);
|
||||
void setXYVar(const Index i, const Index j, const DMat &m);
|
||||
@ -83,6 +84,7 @@ public:
|
||||
const DMat & getXYVar(const Index i, const Index j) const;
|
||||
DVec getXError(const Index i) const;
|
||||
DVec getYError(const Index j) const;
|
||||
DMat getTable(const Index i, const Index j);
|
||||
// get total fit variance matrix and its pseudo-inverse
|
||||
const DMat & getFitVarMat(void);
|
||||
const DMat & getFitVarMatPInv(void);
|
||||
@ -97,6 +99,10 @@ public:
|
||||
template <typename... Ts>
|
||||
FitResult fit(Minimizer &minimizer, const DVec &init,
|
||||
const DoubleModel &model, const Ts... models);
|
||||
// residuals
|
||||
XYStatData getResiduals(const FitResult &fit);
|
||||
XYStatData getPartialResiduals(const FitResult &fit, const DVec &ref,
|
||||
const Index i);
|
||||
protected:
|
||||
// create data
|
||||
virtual void createXData(const std::string name, const Index nData);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user