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

XY datatypes: tools to get residuals or partial residuals from a fit

Browse Source
This commit is contained in:
Antonin Portelli 2015-07-07 18:50:03 +01:00
parent 59bb3fb78c
commit cebf2334fa
6 changed files with 167 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
lib/Chi2Function.cpp
View File

@ -29,7 +29,7 @@ using namespace Latan;
******************************************************************************/ ******************************************************************************/
// constructors //////////////////////////////////////////////////////////////// // constructors ////////////////////////////////////////////////////////////////
Chi2Function::Chi2Function(const XYStatData &data) Chi2Function::Chi2Function(const XYStatData &data)
: data_(data) : data_(&data)
, buffer_(new Chi2FunctionBuffer) , buffer_(new Chi2FunctionBuffer)
{ {
resizeBuffer(); resizeBuffer();
@ -50,7 +50,7 @@ Index Chi2Function::getNArg(void) const
LATAN_ERROR(Memory, "no model set"); LATAN_ERROR(Memory, "no model set");
} }
return nPar_ + data_.getStatXDim()*data_.getNFitPoint(); return nPar_ + data_->getStatXDim()*data_->getNFitPoint();
} }
Index Chi2Function::getNDof(void) const Index Chi2Function::getNDof(void) const
@ -60,7 +60,7 @@ Index Chi2Function::getNDof(void) const
LATAN_ERROR(Memory, "no model set"); LATAN_ERROR(Memory, "no model set");
} }
return data_.getYDim()*data_.getNFitPoint() - nPar_; return data_->getYDim()*data_->getNFitPoint() - nPar_;
} }
Index Chi2Function::getNPar(void) const Index Chi2Function::getNPar(void) const
@ -77,15 +77,15 @@ void Chi2Function::setModel(const DoubleModel &model, const Index j)
{ {
typedef decltype(model_.size()) size_type; typedef decltype(model_.size()) size_type;
if (static_cast<Index>(model_.size()) != data_.getYDim()) if (static_cast<Index>(model_.size()) != data_->getYDim())
{ {
model_.resize(static_cast<size_type>(data_.getYDim())); model_.resize(static_cast<size_type>(data_->getYDim()));
} }
if (model.getNArg() != data_.getXDim()) if (model.getNArg() != data_->getXDim())
{ {
LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch"); LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
} }
for (unsigned int l = 0; l < data_.getYDim(); ++l) for (unsigned int l = 0; l < data_->getYDim(); ++l)
{ {
if (model_[l]&&(l != j)) if (model_[l]&&(l != j))
{ {
@ -103,11 +103,11 @@ void Chi2Function::setModel(const vector<const DoubleModel *> &modelVector)
{ {
typedef decltype(model_.size()) size_type; typedef decltype(model_.size()) size_type;
if (static_cast<Index>(model_.size()) != data_.getYDim()) if (static_cast<Index>(model_.size()) != data_->getYDim())
{ {
model_.resize(static_cast<size_type>(data_.getYDim())); model_.resize(static_cast<size_type>(data_->getYDim()));
} }
if (modelVector.size() != static_cast<size_type>(data_.getYDim())) if (modelVector.size() != static_cast<size_type>(data_->getYDim()))
{ {
LATAN_ERROR(Size, "number of models and y-dimension mismatch"); LATAN_ERROR(Size, "number of models and y-dimension mismatch");
} }
@ -117,7 +117,7 @@ void Chi2Function::setModel(const vector<const DoubleModel *> &modelVector)
{ {
LATAN_ERROR(Memory, "trying to set a null model"); LATAN_ERROR(Memory, "trying to set a null model");
} }
if (modelVector[j]->getNArg() != data_.getXDim()) if (modelVector[j]->getNArg() != data_->getXDim())
{ {
LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch"); LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
} }
@ -139,24 +139,24 @@ void Chi2Function::resizeBuffer(void) const
{ {
Index size; Index size;
size = (data_.getYDim() + data_.getStatXDim())*data_.getNFitPoint(); size = (data_->getYDim() + data_->getStatXDim())*data_->getNFitPoint();
buffer_->v.setConstant(size, 0.0); buffer_->v.setConstant(size, 0.0);
buffer_->x.setConstant(data_.getXDim(), 0.0); buffer_->x.setConstant(data_->getXDim(), 0.0);
buffer_->invVar.setConstant(size, size, 0.0); buffer_->invVar.setConstant(size, size, 0.0);
buffer_->xInd.setConstant(data_.getStatXDim(), 0); buffer_->xInd.setConstant(data_->getStatXDim(), 0);
buffer_->dInd.setConstant(data_.getNFitPoint(), 0); buffer_->dInd.setConstant(data_->getNFitPoint(), 0);
} }
// compute variance matrix inverse ///////////////////////////////////////////// // compute variance matrix inverse /////////////////////////////////////////////
void Chi2Function::setVarianceBlock(const Index l1, const Index l2, void Chi2Function::setVarianceBlock(const Index l1, const Index l2,
ConstBlock<MatBase<double>> m) const ConstBlock<MatBase<double>> m) const
{ {
const Index nPoint = data_.getNFitPoint(); const Index nPoint = data_->getNFitPoint();
FOR_VEC(buffer_->dInd, k2) FOR_VEC(buffer_->dInd, k2)
FOR_VEC(buffer_->dInd, k1) FOR_VEC(buffer_->dInd, k1)
{ {
if (data_.isDataCorrelated(buffer_->dInd(k1), buffer_->dInd(k2))) if (data_->isDataCorrelated(buffer_->dInd(k1), buffer_->dInd(k2)))
{ {
buffer_->invVar(l1*nPoint + k1, l2*nPoint + k2) = buffer_->invVar(l1*nPoint + k1, l2*nPoint + k2) =
m(buffer_->dInd(k1), buffer_->dInd(k2)); m(buffer_->dInd(k1), buffer_->dInd(k2));
@ -166,11 +166,11 @@ void Chi2Function::setVarianceBlock(const Index l1, const Index l2,
void Chi2Function::initBuffer(void) const void Chi2Function::initBuffer(void) const
{ {
const Index xDim = data_.getXDim(); const Index xDim = data_->getXDim();
const Index statXDim = data_.getStatXDim(); const Index statXDim = data_->getStatXDim();
const Index yDim = data_.getYDim(); const Index yDim = data_->getYDim();
const Index nData = data_.getNData(); const Index nData = data_->getNData();
const Index nPoint = data_.getNFitPoint(); const Index nPoint = data_->getNFitPoint();
Index is, kf; Index is, kf;
// resize buffer // resize buffer
@ -180,7 +180,7 @@ void Chi2Function::initBuffer(void) const
is = 0; is = 0;
for (Index i = 0; i < xDim; ++i) for (Index i = 0; i < xDim; ++i)
{ {
if (!data_.isXExact(i)) if (!data_->isXExact(i))
{ {
buffer_->xInd(is) = i; buffer_->xInd(is) = i;
is++; is++;
@ -189,7 +189,7 @@ void Chi2Function::initBuffer(void) const
kf = 0; kf = 0;
for (Index k = 0; k < nData; ++k) for (Index k = 0; k < nData; ++k)
{ {
if (data_.isFitPoint(k)) if (data_->isFitPoint(k))
{ {
buffer_->dInd(kf) = k; buffer_->dInd(kf) = k;
kf++; kf++;
@ -200,9 +200,9 @@ void Chi2Function::initBuffer(void) const
for (Index j2 = 0; j2 < yDim; ++j2) for (Index j2 = 0; j2 < yDim; ++j2)
for (Index j1 = 0; j1 < yDim; ++j1) for (Index j1 = 0; j1 < yDim; ++j1)
{ {
if (data_.isYYCorrelated(j1, j2)) if (data_->isYYCorrelated(j1, j2))
{ {
setVarianceBlock(j1, j2, data_.yyVar(j1, j2)); setVarianceBlock(j1, j2, data_->yyVar(j1, j2));
} }
} }
@ -210,10 +210,10 @@ void Chi2Function::initBuffer(void) const
FOR_VEC(buffer_->xInd, i2) FOR_VEC(buffer_->xInd, i2)
FOR_VEC(buffer_->xInd, i1) FOR_VEC(buffer_->xInd, i1)
{ {
if (data_.isXXCorrelated(buffer_->xInd(i1), buffer_->xInd(i2))) if (data_->isXXCorrelated(buffer_->xInd(i1), buffer_->xInd(i2)))
{ {
setVarianceBlock(i1 + yDim, i2 + yDim, setVarianceBlock(i1 + yDim, i2 + yDim,
data_.xxVar(buffer_->xInd(i1), buffer_->xInd(i2))); data_->xxVar(buffer_->xInd(i1), buffer_->xInd(i2)));
} }
} }
@ -221,9 +221,9 @@ void Chi2Function::initBuffer(void) const
FOR_VEC(buffer_->xInd, i) FOR_VEC(buffer_->xInd, i)
for (Index j = 0; j < yDim; ++j) for (Index j = 0; j < yDim; ++j)
{ {
if (data_.isYXCorrelated(j, buffer_->xInd(i))) if (data_->isYXCorrelated(j, buffer_->xInd(i)))
{ {
setVarianceBlock(j, i + yDim, data_.yxVar(j, buffer_->xInd(i))); setVarianceBlock(j, i + yDim, data_->yxVar(j, buffer_->xInd(i)));
} }
} }
auto lowerYX = buffer_->invVar.block(yDim*nPoint, 0, yDim*statXDim, auto lowerYX = buffer_->invVar.block(yDim*nPoint, 0, yDim*statXDim,
@ -247,9 +247,9 @@ double Chi2Function::operator()(const double *arg) const
LATAN_ERROR(Memory, "null model"); LATAN_ERROR(Memory, "null model");
} }
const Index xDim = data_.getXDim(); const Index xDim = data_->getXDim();
const Index yDim = data_.getYDim(); const Index yDim = data_->getYDim();
const Index nPoint = data_.getNFitPoint(); const Index nPoint = data_->getNFitPoint();
Index is; Index is;
ConstMap<DVec> xi(arg + nPar_, getNArg() - nPar_, 1); ConstMap<DVec> xi(arg + nPar_, getNArg() - nPar_, 1);
double res; double res;
@ -265,7 +265,7 @@ double Chi2Function::operator()(const double *arg) const
FOR_VEC(buffer_->dInd, k) FOR_VEC(buffer_->dInd, k)
{ {
const DoubleModel *f = model_[static_cast<size_type>(j)]; const DoubleModel *f = model_[static_cast<size_type>(j)];
double f_jk, y_jk = data_.y(j, buffer_->dInd(k)); double f_jk, y_jk = data_->y(j, buffer_->dInd(k));
if (!f) if (!f)
{ {
@ -274,9 +274,9 @@ double Chi2Function::operator()(const double *arg) const
is = 0; is = 0;
for (Index i = 0; i < xDim; ++i) for (Index i = 0; i < xDim; ++i)
{ {
if (data_.isXExact(i)) if (data_->isXExact(i))
{ {
buffer_->x(i) = data_.x(i, buffer_->dInd(k)); buffer_->x(i) = data_->x(i, buffer_->dInd(k));
} }
else else
{ {
@ -293,7 +293,7 @@ double Chi2Function::operator()(const double *arg) const
FOR_VEC(buffer_->xInd, i) FOR_VEC(buffer_->xInd, i)
FOR_VEC(buffer_->dInd, k) FOR_VEC(buffer_->dInd, k)
{ {
double x_ik = data_.x(buffer_->xInd(i), buffer_->dInd(k)); double x_ik = data_->x(buffer_->xInd(i), buffer_->dInd(k));
double xi_ik = xi(i*nPoint + k); double xi_ik = xi(i*nPoint + k);
buffer_->v(yDim*nPoint + i*nPoint + k) = xi_ik - x_ik; buffer_->v(yDim*nPoint + i*nPoint + k) = xi_ik - x_ik;

2
lib/Chi2Function.hpp
View File

@ -69,7 +69,7 @@ private:
ConstBlock<MatBase<double>> m) const; ConstBlock<MatBase<double>> m) const;
void initBuffer(void) const; void initBuffer(void) const;
private: private:
const XYStatData &data_; const XYStatData *data_;
std::shared_ptr<Chi2FunctionBuffer> buffer_; std::shared_ptr<Chi2FunctionBuffer> buffer_;
std::vector<const DoubleModel *> model_; std::vector<const DoubleModel *> model_;
Index nPar_{-1}; Index nPar_{-1};

71
lib/XYSampleData.cpp
View File

@ -64,6 +64,22 @@ const DoubleFunctionSample & SampleFitResult::getModel(
return model_[static_cast<unsigned int>(j)]; return model_[static_cast<unsigned int>(j)];
} }
FitResult SampleFitResult::getFitResult(const Index s) const
{
FitResult fit;
fit = (*this)[s];
fit.chi2_ = getChi2();
fit.nDof_ = static_cast<Index>(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 * * XYSampleData implementation *
******************************************************************************/ ******************************************************************************/
@ -294,3 +310,58 @@ void XYSampleData::setDataToSample(const Index s)
data_.x() = x_[s]; data_.x() = x_[s];
data_.y() = y_[s]; data_.y() = y_[s];
} }
// residuals ///////////////////////////////////////////////////////////////////
XYSampleData XYSampleData::getResiduals(const SampleFitResult &fit) const
{
const Index nSample = x_.size();
XYSampleData res(*this);
DMatSample xBuf(nSample, getXDim(), 1), tmp(nSample, 1, 1);
for (Index j = 0; j < res.getYDim(); ++j)
{
const DoubleFunctionSample &f = fit.getModel(_, j);
for (Index k = 0; k < res.getNData(); ++k)
{
xBuf = this->x(_, k);
tmp = this->y(j, k);
FOR_STAT_ARRAY(xBuf, s)
{
tmp[s](0) -= f[s](xBuf[s].transpose());
}
res.y(j, k) = tmp;
}
}
return res;
}
XYSampleData XYSampleData::getPartialResiduals(const SampleFitResult &fit,
const DVec &x,
const Index i) const
{
const Index nSample = x_.size();
XYSampleData res(*this);
DMatSample xBuf(nSample, getXDim(), 1), tmp(nSample, 1, 1);
DVec buf(x);
for (Index j = 0; j < res.getYDim(); ++j)
{
const DoubleFunctionSample &f = fit.getModel(_, j);
for (Index k = 0; k < res.getNData(); ++k)
{
xBuf = this->x(_, k);
tmp = this->y(j, k);
FOR_STAT_ARRAY(xBuf, s)
{
buf(i) = xBuf[s](i);
tmp[s](0) -= f[s](xBuf[s].transpose()) - f[s](buf);
}
res.y(j, k) = tmp;
}
}
return res;
}

17
lib/XYSampleData.hpp
View File

@ -51,6 +51,7 @@ public:
const Index j = 0) const; const Index j = 0) const;
const DoubleFunctionSample & getModel(const PlaceHolder ph, const DoubleFunctionSample & getModel(const PlaceHolder ph,
const Index j = 0) const; const Index j = 0) const;
FitResult getFitResult(const Index s = central) const;
private: private:
DSample chi2_; DSample chi2_;
double nDof_{0.}; double nDof_{0.};
@ -59,7 +60,9 @@ private:
/****************************************************************************** /******************************************************************************
* XYSampleData * * XYSampleData *
******************************************************************************/ ******************************************************************************
* index convention: i: X, j: Y, k: data
*/
class XYSampleData: public FitInterface class XYSampleData: public FitInterface
{ {
public: public:
@ -88,18 +91,22 @@ public:
SampleBlock y(const PlaceHolder ph1 = _, const PlaceHolder ph2 = _); SampleBlock y(const PlaceHolder ph1 = _, const PlaceHolder ph2 = _);
ConstSampleBlock y(const PlaceHolder ph1 = _, ConstSampleBlock y(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _) const; const PlaceHolder ph2 = _) const;
SampleBlock y(const Index i, const PlaceHolder ph2 = _); SampleBlock y(const Index j, const PlaceHolder ph2 = _);
ConstSampleBlock y(const Index i, const PlaceHolder ph2 = _) const; ConstSampleBlock y(const Index j, const PlaceHolder ph2 = _) const;
SampleBlock y(const PlaceHolder ph1, const Index k); SampleBlock y(const PlaceHolder ph1, const Index k);
ConstSampleBlock y(const PlaceHolder ph1, const Index k) const; ConstSampleBlock y(const PlaceHolder ph1, const Index k) const;
SampleBlock y(const Index i, const Index k); SampleBlock y(const Index j, const Index k);
ConstSampleBlock y(const Index i, const Index k) const; ConstSampleBlock y(const Index j, const Index k) const;
// fit // fit
SampleFitResult fit(Minimizer &minimizer, const DVec &init, SampleFitResult fit(Minimizer &minimizer, const DVec &init,
const std::vector<const DoubleModel *> &modelVector); const std::vector<const DoubleModel *> &modelVector);
template <typename... Mods> template <typename... Mods>
SampleFitResult fit(Minimizer &minimizer, const DVec &init, SampleFitResult fit(Minimizer &minimizer, const DVec &init,
const DoubleModel &model, const Mods... models); const DoubleModel &model, const Mods... models);
// residuals
XYSampleData getResiduals(const SampleFitResult &fit) const;
XYSampleData getPartialResiduals(const SampleFitResult &fit, const DVec &x,
const Index j) const;
private: private:
void setDataToSample(const Index s); void setDataToSample(const Index s);
private: private:

38
lib/XYStatData.cpp
View File

@ -268,3 +268,41 @@ FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
return result; return result;
} }
// residuals ///////////////////////////////////////////////////////////////////
XYStatData XYStatData::getResiduals(const FitResult &fit) const
{
XYStatData res(*this);
for (Index j = 0; j < res.getYDim(); ++j)
{
const DoubleFunction &f = fit.getModel(j);
for (Index k = 0; k < res.getNData(); ++k)
{
res.y(j, k) -= f(res.x(_, k).transpose());
}
}
return res;
}
XYStatData XYStatData::getPartialResiduals(const FitResult &fit, const DVec &x,
const Index i) const
{
XYStatData res(*this);
DVec buf(x), xk;
for (Index j = 0; j < res.getYDim(); ++j)
{
const DoubleFunction &f = fit.getModel(j);
for (Index k = 0; k < res.getNData(); ++k)
{
buf(i) = res.x(i, k);
res.y(j, k) -= f(res.x(_, k).transpose()) - f(buf);
}
}
return res;
}

9
lib/XYStatData.hpp
View File

@ -37,6 +37,7 @@ BEGIN_LATAN_NAMESPACE
class FitResult: public DVec class FitResult: public DVec
{ {
friend class XYStatData; friend class XYStatData;
friend class SampleFitResult;
public: public:
// constructors // constructors
FitResult(void) = default; FitResult(void) = default;
@ -56,7 +57,9 @@ private:
/****************************************************************************** /******************************************************************************
* object for X vs. Y statistical data * * object for X vs. Y statistical data *
******************************************************************************/ ******************************************************************************
* index convention: i: X, j: Y, k: data
*/
class XYStatData: public FitInterface class XYStatData: public FitInterface
{ {
public: public:
@ -112,6 +115,10 @@ public:
template <typename... Ts> template <typename... Ts>
FitResult fit(Minimizer &minimizer, const DVec &init, FitResult fit(Minimizer &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models); const DoubleModel &model, const Ts... models);
// residuals
XYStatData getResiduals(const FitResult &fit) const;
XYStatData getPartialResiduals(const FitResult &fit, const DVec &x,
const Index j) const;
private: private:
DMat x_, y_; DMat x_, y_;