mirror of
https://github.com/aportelli/LatAnalyze.git
synced 2025-05-01 01:25:56 +01:00
code cleaning
This commit is contained in:
parent
9b9c86cf72
commit
d699e9e564
@ -1,325 +0,0 @@
|
|||||||
/*
|
|
||||||
* Chi2Function.cpp, part of LatAnalyze 3
|
|
||||||
*
|
|
||||||
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#include <LatAnalyze/Chi2Function.hpp>
|
|
||||||
#include <LatAnalyze/includes.hpp>
|
|
||||||
#include <LatAnalyze/XYStatData.hpp>
|
|
||||||
|
|
||||||
using namespace std;
|
|
||||||
using namespace Latan;
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* Chi2Function implementation *
|
|
||||||
******************************************************************************/
|
|
||||||
// constructors ////////////////////////////////////////////////////////////////
|
|
||||||
Chi2Function::Chi2Function(const XYStatData &data)
|
|
||||||
: data_(&data)
|
|
||||||
, buffer_(new Chi2FunctionBuffer)
|
|
||||||
{
|
|
||||||
resizeBuffer();
|
|
||||||
}
|
|
||||||
|
|
||||||
Chi2Function::Chi2Function(const XYStatData &data,
|
|
||||||
const vector<const DoubleModel *> &modelVector)
|
|
||||||
: Chi2Function(data)
|
|
||||||
{
|
|
||||||
setModel(modelVector);
|
|
||||||
}
|
|
||||||
|
|
||||||
// access //////////////////////////////////////////////////////////////////////
|
|
||||||
Index Chi2Function::getNArg(void) const
|
|
||||||
{
|
|
||||||
if (nPar_ < 0)
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Memory, "no model set");
|
|
||||||
}
|
|
||||||
|
|
||||||
return nPar_ + data_->getStatXDim()*data_->getNFitPoint();
|
|
||||||
}
|
|
||||||
|
|
||||||
Index Chi2Function::getNDof(void) const
|
|
||||||
{
|
|
||||||
if (nPar_ < 0)
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Memory, "no model set");
|
|
||||||
}
|
|
||||||
|
|
||||||
return data_->getYDim()*data_->getNFitPoint() - nPar_;
|
|
||||||
}
|
|
||||||
|
|
||||||
Index Chi2Function::getNPar(void) const
|
|
||||||
{
|
|
||||||
if (nPar_ < 0)
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Memory, "no model set");
|
|
||||||
}
|
|
||||||
|
|
||||||
return nPar_;
|
|
||||||
}
|
|
||||||
|
|
||||||
void Chi2Function::setModel(const DoubleModel &model, const Index j)
|
|
||||||
{
|
|
||||||
typedef decltype(model_.size()) size_type;
|
|
||||||
|
|
||||||
if (static_cast<Index>(model_.size()) != data_->getYDim())
|
|
||||||
{
|
|
||||||
model_.resize(static_cast<size_type>(data_->getYDim()));
|
|
||||||
}
|
|
||||||
if (model.getNArg() != data_->getXDim())
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
|
|
||||||
}
|
|
||||||
for (unsigned int l = 0; l < data_->getYDim(); ++l)
|
|
||||||
{
|
|
||||||
if (model_[l]&&(l != j))
|
|
||||||
{
|
|
||||||
if (model_[l]->getNPar() != model.getNPar())
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Size, "model number of parameter mismatch");
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
model_[static_cast<size_type>(j)] = &model;
|
|
||||||
nPar_ = model.getNPar();
|
|
||||||
}
|
|
||||||
|
|
||||||
void Chi2Function::setModel(const vector<const DoubleModel *> &modelVector)
|
|
||||||
{
|
|
||||||
typedef decltype(model_.size()) size_type;
|
|
||||||
|
|
||||||
if (static_cast<Index>(model_.size()) != data_->getYDim())
|
|
||||||
{
|
|
||||||
model_.resize(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");
|
|
||||||
}
|
|
||||||
for (unsigned int j = 0; j < model_.size(); ++j)
|
|
||||||
{
|
|
||||||
if (!modelVector[j])
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Memory, "trying to set a null model");
|
|
||||||
}
|
|
||||||
if (modelVector[j]->getNArg() != data_->getXDim())
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
|
|
||||||
}
|
|
||||||
model_[static_cast<size_type>(j)] = modelVector[j];
|
|
||||||
if (modelVector[j]->getNPar() != modelVector[0]->getNPar())
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Size, "model number of parameter mismatch");
|
|
||||||
}
|
|
||||||
}
|
|
||||||
nPar_ = modelVector[0]->getNPar();
|
|
||||||
}
|
|
||||||
|
|
||||||
void Chi2Function::requestInit(void) const
|
|
||||||
{
|
|
||||||
buffer_->isInit = false;
|
|
||||||
}
|
|
||||||
|
|
||||||
void Chi2Function::resizeBuffer(void) const
|
|
||||||
{
|
|
||||||
Index size;
|
|
||||||
|
|
||||||
size = (data_->getYDim() + data_->getStatXDim())*data_->getNFitPoint();
|
|
||||||
buffer_->v.setConstant(size, 0.0);
|
|
||||||
buffer_->x.setConstant(data_->getXDim(), 0.0);
|
|
||||||
buffer_->invVar.setConstant(size, size, 0.0);
|
|
||||||
buffer_->xInd.setConstant(data_->getStatXDim(), 0);
|
|
||||||
buffer_->dInd.setConstant(data_->getNFitPoint(), 0);
|
|
||||||
}
|
|
||||||
|
|
||||||
// compute variance matrix inverse /////////////////////////////////////////////
|
|
||||||
void Chi2Function::setVarianceBlock(const Index l1, const Index l2,
|
|
||||||
ConstBlock<MatBase<double>> m) const
|
|
||||||
{
|
|
||||||
const Index nPoint = data_->getNFitPoint();
|
|
||||||
|
|
||||||
FOR_VEC(buffer_->dInd, k2)
|
|
||||||
FOR_VEC(buffer_->dInd, k1)
|
|
||||||
{
|
|
||||||
if (data_->isDataCorrelated(buffer_->dInd(k1), buffer_->dInd(k2)))
|
|
||||||
{
|
|
||||||
buffer_->invVar(l1*nPoint + k1, l2*nPoint + k2) =
|
|
||||||
m(buffer_->dInd(k1), buffer_->dInd(k2));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void Chi2Function::initBuffer(void) const
|
|
||||||
{
|
|
||||||
const Index xDim = data_->getXDim();
|
|
||||||
const Index statXDim = data_->getStatXDim();
|
|
||||||
const Index yDim = data_->getYDim();
|
|
||||||
const Index nData = data_->getNData();
|
|
||||||
const Index nPoint = data_->getNFitPoint();
|
|
||||||
Index is, kf;
|
|
||||||
|
|
||||||
// resize buffer
|
|
||||||
resizeBuffer();
|
|
||||||
|
|
||||||
// build index tables
|
|
||||||
is = 0;
|
|
||||||
for (Index i = 0; i < xDim; ++i)
|
|
||||||
{
|
|
||||||
if (!data_->isXExact(i))
|
|
||||||
{
|
|
||||||
buffer_->xInd(is) = i;
|
|
||||||
is++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
kf = 0;
|
|
||||||
for (Index k = 0; k < nData; ++k)
|
|
||||||
{
|
|
||||||
if (data_->isFitPoint(k))
|
|
||||||
{
|
|
||||||
buffer_->dInd(kf) = k;
|
|
||||||
kf++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// set y/y variance matrix
|
|
||||||
for (Index j2 = 0; j2 < yDim; ++j2)
|
|
||||||
for (Index j1 = 0; j1 < yDim; ++j1)
|
|
||||||
{
|
|
||||||
if (data_->isYYCorrelated(j1, j2))
|
|
||||||
{
|
|
||||||
setVarianceBlock(j1, j2, data_->yyVar(j1, j2));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// set x/x variance matrix
|
|
||||||
FOR_VEC(buffer_->xInd, i2)
|
|
||||||
FOR_VEC(buffer_->xInd, i1)
|
|
||||||
{
|
|
||||||
if (data_->isXXCorrelated(buffer_->xInd(i1), buffer_->xInd(i2)))
|
|
||||||
{
|
|
||||||
setVarianceBlock(i1 + yDim, i2 + yDim,
|
|
||||||
data_->xxVar(buffer_->xInd(i1), buffer_->xInd(i2)));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// set y/x variance matrix
|
|
||||||
FOR_VEC(buffer_->xInd, i)
|
|
||||||
for (Index j = 0; j < yDim; ++j)
|
|
||||||
{
|
|
||||||
if (data_->isYXCorrelated(j, buffer_->xInd(i)))
|
|
||||||
{
|
|
||||||
setVarianceBlock(j, i + yDim, data_->yxVar(j, buffer_->xInd(i)));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
auto lowerYX = buffer_->invVar.block(yDim*nPoint, 0, yDim*statXDim,
|
|
||||||
yDim*nPoint);
|
|
||||||
auto upperYX = buffer_->invVar.block(0, yDim*nPoint, yDim*nPoint,
|
|
||||||
yDim*statXDim);
|
|
||||||
lowerYX = upperYX.transpose();
|
|
||||||
|
|
||||||
// inversion
|
|
||||||
buffer_->invVar = buffer_->invVar.pInverse(data_->getSvdTolerance());
|
|
||||||
buffer_->isInit = true;
|
|
||||||
}
|
|
||||||
|
|
||||||
// function call ///////////////////////////////////////////////////////////////
|
|
||||||
double Chi2Function::operator()(const double *arg) const
|
|
||||||
{
|
|
||||||
typedef decltype(model_.size()) size_type;
|
|
||||||
|
|
||||||
if (!model_[0])
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Memory, "null model");
|
|
||||||
}
|
|
||||||
|
|
||||||
const Index xDim = data_->getXDim();
|
|
||||||
const Index yDim = data_->getYDim();
|
|
||||||
const Index nPoint = data_->getNFitPoint();
|
|
||||||
Index is;
|
|
||||||
ConstMap<DVec> xi(arg + nPar_, getNArg() - nPar_, 1);
|
|
||||||
double res;
|
|
||||||
|
|
||||||
// initialize buffers if necessary
|
|
||||||
if (!buffer_->isInit)
|
|
||||||
{
|
|
||||||
initBuffer();
|
|
||||||
}
|
|
||||||
|
|
||||||
// set the upper part of v: f_j(xi, p) - y_j
|
|
||||||
for (Index j = 0; j < yDim; ++j)
|
|
||||||
FOR_VEC(buffer_->dInd, k)
|
|
||||||
{
|
|
||||||
const DoubleModel *f = model_[static_cast<size_type>(j)];
|
|
||||||
double f_jk, y_jk = data_->y(j, buffer_->dInd(k));
|
|
||||||
|
|
||||||
if (!f)
|
|
||||||
{
|
|
||||||
LATAN_ERROR(Memory, "null model");
|
|
||||||
}
|
|
||||||
is = 0;
|
|
||||||
for (Index i = 0; i < xDim; ++i)
|
|
||||||
{
|
|
||||||
if (data_->isXExact(i))
|
|
||||||
{
|
|
||||||
buffer_->x(i) = data_->x(i, buffer_->dInd(k));
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
buffer_->x(i) = xi(is*nPoint + k);
|
|
||||||
is++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// call model on double pointers to avoid any copy
|
|
||||||
f_jk = (*f)(buffer_->x.data(), arg);
|
|
||||||
buffer_->v(j*nPoint + k) = f_jk - y_jk;
|
|
||||||
}
|
|
||||||
|
|
||||||
// set the down part of v: xi_ik - x_ik
|
|
||||||
FOR_VEC(buffer_->xInd, i)
|
|
||||||
FOR_VEC(buffer_->dInd, k)
|
|
||||||
{
|
|
||||||
double x_ik = data_->x(buffer_->xInd(i), buffer_->dInd(k));
|
|
||||||
double xi_ik = xi(i*nPoint + k);
|
|
||||||
|
|
||||||
buffer_->v(yDim*nPoint + i*nPoint + k) = xi_ik - x_ik;
|
|
||||||
}
|
|
||||||
|
|
||||||
// compute result
|
|
||||||
res = buffer_->v.dot(buffer_->invVar*buffer_->v);
|
|
||||||
|
|
||||||
return res;
|
|
||||||
}
|
|
||||||
|
|
||||||
// DoubleFunction factory //////////////////////////////////////////////////////
|
|
||||||
DoubleFunction Chi2Function::makeFunction(const bool makeHardCopy) const
|
|
||||||
{
|
|
||||||
DoubleFunction res;
|
|
||||||
|
|
||||||
if (makeHardCopy)
|
|
||||||
{
|
|
||||||
Chi2Function copy(*this);
|
|
||||||
|
|
||||||
res.setFunction([copy](const double *p){return copy(p);}, getNArg());
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
res.setFunction([this](const double *p){return (*this)(p);}, getNArg());
|
|
||||||
}
|
|
||||||
|
|
||||||
return res;
|
|
||||||
}
|
|
||||||
@ -1,80 +0,0 @@
|
|||||||
/*
|
|
||||||
* Chi2Function.hpp, part of LatAnalyze 3
|
|
||||||
*
|
|
||||||
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#ifndef Latan_Chi2Function_hpp_
|
|
||||||
#define Latan_Chi2Function_hpp_
|
|
||||||
|
|
||||||
#include <LatAnalyze/Global.hpp>
|
|
||||||
#include <LatAnalyze/Function.hpp>
|
|
||||||
#include <LatAnalyze/Model.hpp>
|
|
||||||
|
|
||||||
BEGIN_LATAN_NAMESPACE
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* <Class> *
|
|
||||||
******************************************************************************/
|
|
||||||
// forward declaration of XYStatData
|
|
||||||
class XYStatData;
|
|
||||||
|
|
||||||
class Chi2Function: public DoubleFunctionFactory
|
|
||||||
{
|
|
||||||
private:
|
|
||||||
struct Chi2FunctionBuffer
|
|
||||||
{
|
|
||||||
DVec v, x;
|
|
||||||
DMat invVar;
|
|
||||||
bool isInit{false};
|
|
||||||
Vec<Index> xInd, dInd;
|
|
||||||
};
|
|
||||||
public:
|
|
||||||
// constructor
|
|
||||||
explicit Chi2Function(const XYStatData &data);
|
|
||||||
Chi2Function(const XYStatData &data,
|
|
||||||
const std::vector<const DoubleModel *> &modelVector);
|
|
||||||
// destructor
|
|
||||||
virtual ~Chi2Function(void) = default;
|
|
||||||
// access
|
|
||||||
virtual Index getNArg(void) const;
|
|
||||||
Index getNDof(void) const;
|
|
||||||
Index getNPar(void) const;
|
|
||||||
void setModel(const DoubleModel &model, const Index j = 0);
|
|
||||||
void setModel(const std::vector<const DoubleModel *> &modelVector);
|
|
||||||
void requestInit(void) const;
|
|
||||||
// function call
|
|
||||||
double operator()(const double *arg) const;
|
|
||||||
// factory
|
|
||||||
virtual DoubleFunction makeFunction(const bool makeHardCopy = true) const;
|
|
||||||
private:
|
|
||||||
// access
|
|
||||||
void resizeBuffer(void) const;
|
|
||||||
|
|
||||||
// compute variance matrix inverse
|
|
||||||
void setVarianceBlock(const Index l1, const Index l2,
|
|
||||||
ConstBlock<MatBase<double>> m) const;
|
|
||||||
void initBuffer(void) const;
|
|
||||||
private:
|
|
||||||
const XYStatData *data_;
|
|
||||||
std::shared_ptr<Chi2FunctionBuffer> buffer_;
|
|
||||||
std::vector<const DoubleModel *> model_;
|
|
||||||
Index nPar_{-1};
|
|
||||||
};
|
|
||||||
|
|
||||||
END_LATAN_NAMESPACE
|
|
||||||
|
|
||||||
#endif // Latan_Chi2Function_hpp_
|
|
||||||
@ -1,209 +0,0 @@
|
|||||||
/*
|
|
||||||
* FitInterface.cpp, part of LatAnalyze 3
|
|
||||||
*
|
|
||||||
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#include <LatAnalyze/FitInterface.hpp>
|
|
||||||
#include <LatAnalyze/includes.hpp>
|
|
||||||
|
|
||||||
using namespace std;
|
|
||||||
using namespace Latan;
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* FitInterface implementation *
|
|
||||||
******************************************************************************/
|
|
||||||
// constructors ////////////////////////////////////////////////////////////////
|
|
||||||
FitInterface::FitInterface(const Index nData, const Index xDim,
|
|
||||||
const Index yDim)
|
|
||||||
{
|
|
||||||
resize(nData, xDim, yDim);
|
|
||||||
}
|
|
||||||
|
|
||||||
// access //////////////////////////////////////////////////////////////////////
|
|
||||||
void FitInterface::assumeXExact(const Index i, const bool isExact)
|
|
||||||
{
|
|
||||||
isXExact_(i) = (isExact) ? 1 : 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::assumeXXCorrelated(const Index i1, const Index i2,
|
|
||||||
const bool isCorrelated)
|
|
||||||
{
|
|
||||||
int val = (isCorrelated) ? 1 : 0;
|
|
||||||
|
|
||||||
isXXCorr_(i1, i2) = val;
|
|
||||||
isXXCorr_(i2, i1) = val;
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::assumeYYCorrelated(const Index j1, const Index j2,
|
|
||||||
const bool isCorrelated)
|
|
||||||
{
|
|
||||||
int val = (isCorrelated) ? 1 : 0;
|
|
||||||
|
|
||||||
isYYCorr_(j1, j2) = val;
|
|
||||||
isYYCorr_(j2, j1) = val;
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::assumeYXCorrelated(const Index j, const Index i,
|
|
||||||
const bool isCorrelated)
|
|
||||||
{
|
|
||||||
int val = (isCorrelated) ? 1 : 0;
|
|
||||||
|
|
||||||
isYXCorr_(j, i) = val;
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::assumeDataCorrelated(const Index k1, const Index k2,
|
|
||||||
const bool isCorrelated)
|
|
||||||
{
|
|
||||||
int val = (isCorrelated) ? 1 : 0;
|
|
||||||
|
|
||||||
isDataCorr_(k1, k2) = val;
|
|
||||||
isDataCorr_(k2, k1) = val;
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::assumeDataCorrelated(const bool isCorrelated)
|
|
||||||
{
|
|
||||||
FOR_MAT(isDataCorr_, k1, k2)
|
|
||||||
{
|
|
||||||
if (k1 != k2)
|
|
||||||
{
|
|
||||||
assumeDataCorrelated(k1, k2, isCorrelated);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::fitPoint(const Index i, const bool isFitPoint)
|
|
||||||
{
|
|
||||||
isFitPoint_(i) = (isFitPoint) ? 1 : 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::fitPointRange(const Index k1, const Index k2,
|
|
||||||
const bool isFitPoint)
|
|
||||||
{
|
|
||||||
int size = static_cast<int>(k2-k1+1);
|
|
||||||
|
|
||||||
isFitPoint_.segment(k1, size) = IVec::Constant(size, (isFitPoint) ? 1 : 0);
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::fitAllPoints(const bool isFitPoint)
|
|
||||||
{
|
|
||||||
fitPointRange(0, getNData()-1, isFitPoint);
|
|
||||||
}
|
|
||||||
|
|
||||||
Index FitInterface::getNData(void) const
|
|
||||||
{
|
|
||||||
return isFitPoint_.size();
|
|
||||||
}
|
|
||||||
|
|
||||||
Index FitInterface::getNFitPoint(void) const
|
|
||||||
{
|
|
||||||
return isFitPoint_.sum();
|
|
||||||
}
|
|
||||||
|
|
||||||
Index FitInterface::getXDim(void) const
|
|
||||||
{
|
|
||||||
return isXXCorr_.rows();
|
|
||||||
}
|
|
||||||
|
|
||||||
Index FitInterface::getYDim(void) const
|
|
||||||
{
|
|
||||||
return isYYCorr_.rows();
|
|
||||||
}
|
|
||||||
|
|
||||||
Index FitInterface::getStatXDim(void) const
|
|
||||||
{
|
|
||||||
return isXExact_.size() - isXExact_.sum();
|
|
||||||
}
|
|
||||||
|
|
||||||
double FitInterface::getSvdTolerance(void) const
|
|
||||||
{
|
|
||||||
return svdTolerance_;
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::setFitInterface(const FitInterface &fitInterface)
|
|
||||||
{
|
|
||||||
if (&fitInterface != this)
|
|
||||||
{
|
|
||||||
isXExact_ = fitInterface.isXExact_;
|
|
||||||
isFitPoint_ = fitInterface.isFitPoint_;
|
|
||||||
isXXCorr_ = fitInterface.isXXCorr_;
|
|
||||||
isYYCorr_ = fitInterface.isYYCorr_;
|
|
||||||
isYXCorr_ = fitInterface.isYXCorr_;
|
|
||||||
isDataCorr_ = fitInterface.isDataCorr_;
|
|
||||||
svdTolerance_ = fitInterface.svdTolerance_;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::setNData(const Index nData)
|
|
||||||
{
|
|
||||||
resize(nData, getXDim(), getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::setXDim(const Index xDim)
|
|
||||||
{
|
|
||||||
resize(getNData(), xDim, getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::setYDim(const Index yDim)
|
|
||||||
{
|
|
||||||
resize(getNData(), getXDim(), yDim);
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::setSvdTolerance(const double tolerance)
|
|
||||||
{
|
|
||||||
svdTolerance_ = tolerance;
|
|
||||||
}
|
|
||||||
|
|
||||||
void FitInterface::resize(const Index nData, const Index xDim, const Index yDim)
|
|
||||||
{
|
|
||||||
isXExact_.setConstant(xDim, 0);
|
|
||||||
isFitPoint_.setConstant(nData, 0);
|
|
||||||
isXXCorr_.setIdentity(xDim, xDim);
|
|
||||||
isYYCorr_.setIdentity(yDim, yDim);
|
|
||||||
isYXCorr_.setConstant(yDim, xDim, 0);
|
|
||||||
isDataCorr_.setIdentity(nData, nData);
|
|
||||||
}
|
|
||||||
|
|
||||||
// test ////////////////////////////////////////////////////////////////////////
|
|
||||||
bool FitInterface::isFitPoint(const Index k) const
|
|
||||||
{
|
|
||||||
return (isFitPoint_[k] == 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
bool FitInterface::isXExact(const Index i) const
|
|
||||||
{
|
|
||||||
return (isXExact_[i] == 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
bool FitInterface::isXXCorrelated(const Index i1, const Index i2) const
|
|
||||||
{
|
|
||||||
return (isXXCorr_(i1, i2) == 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
bool FitInterface::isYYCorrelated(const Index j1, const Index j2) const
|
|
||||||
{
|
|
||||||
return (isYYCorr_(j1, j2) == 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
bool FitInterface::isYXCorrelated(const Index j, const Index i) const
|
|
||||||
{
|
|
||||||
return (isYXCorr_(j, i) == 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
bool FitInterface::isDataCorrelated(const Index k1, const Index k2) const
|
|
||||||
{
|
|
||||||
return (isDataCorr_(k1, k2) == 1);
|
|
||||||
}
|
|
||||||
@ -1,83 +0,0 @@
|
|||||||
/*
|
|
||||||
* FitInterface.hpp, part of LatAnalyze 3
|
|
||||||
*
|
|
||||||
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#ifndef Latan_FitInterface_hpp_
|
|
||||||
#define Latan_FitInterface_hpp_
|
|
||||||
|
|
||||||
#include <LatAnalyze/Global.hpp>
|
|
||||||
#include <LatAnalyze/Minimizer.hpp>
|
|
||||||
|
|
||||||
BEGIN_LATAN_NAMESPACE
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* base class for data fit *
|
|
||||||
******************************************************************************/
|
|
||||||
class FitInterface
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
typedef Minimizer::Verbosity FitVerbosity;
|
|
||||||
public:
|
|
||||||
// constructors
|
|
||||||
FitInterface(void) = default;
|
|
||||||
FitInterface(const Index nData, const Index xDim, const Index yDim);
|
|
||||||
// destructor
|
|
||||||
virtual ~FitInterface(void) = default;
|
|
||||||
// access
|
|
||||||
void assumeXExact(const Index i, const bool isExact = true);
|
|
||||||
void assumeXXCorrelated(const Index i1, const Index i2,
|
|
||||||
const bool isCorrelated = true);
|
|
||||||
void assumeYYCorrelated(const Index j1, const Index j2,
|
|
||||||
const bool isCorrelated = true);
|
|
||||||
void assumeYXCorrelated(const Index j, const Index i,
|
|
||||||
const bool isCorrelated = true);
|
|
||||||
void assumeDataCorrelated(const Index k1, const Index k2,
|
|
||||||
const bool isCorrelated = true);
|
|
||||||
void assumeDataCorrelated(const bool isCorrelated = true);
|
|
||||||
void fitPoint(const Index k, const bool isFitPoint = true);
|
|
||||||
void fitPointRange(const Index k1, const Index k2,
|
|
||||||
const bool isFitPoint = true);
|
|
||||||
void fitAllPoints(const bool isFitPoint = true);
|
|
||||||
Index getNData(void) const;
|
|
||||||
Index getNFitPoint(void) const;
|
|
||||||
Index getXDim(void) const;
|
|
||||||
Index getYDim(void) const;
|
|
||||||
Index getStatXDim(void) const;
|
|
||||||
double getSvdTolerance(void) const;
|
|
||||||
void setFitInterface(const FitInterface &fitInterface);
|
|
||||||
void setNData(const Index nData);
|
|
||||||
void setXDim(const Index xDim);
|
|
||||||
void setYDim(const Index yDim);
|
|
||||||
void setSvdTolerance(const double tolerance);
|
|
||||||
void resize(const Index nData, const Index xDim, const Index yDim);
|
|
||||||
// test
|
|
||||||
bool isFitPoint(const Index k) const;
|
|
||||||
bool isXExact(const Index i) const;
|
|
||||||
bool isXXCorrelated(const Index i1, const Index i2) const;
|
|
||||||
bool isYYCorrelated(const Index j1, const Index j2) const;
|
|
||||||
bool isYXCorrelated(const Index j, const Index i) const;
|
|
||||||
bool isDataCorrelated(const Index k1, const Index k2) const;
|
|
||||||
private:
|
|
||||||
IVec isXExact_, isFitPoint_;
|
|
||||||
IMat isXXCorr_, isYYCorr_, isYXCorr_, isDataCorr_;
|
|
||||||
double svdTolerance_{1.0e-10};
|
|
||||||
};
|
|
||||||
|
|
||||||
END_LATAN_NAMESPACE
|
|
||||||
|
|
||||||
#endif // Latan_FitInterface_hpp_
|
|
||||||
@ -29,7 +29,7 @@ static constexpr double initErr = 0.1;
|
|||||||
static constexpr unsigned int maxTry = 10u;
|
static constexpr unsigned int maxTry = 10u;
|
||||||
|
|
||||||
/******************************************************************************
|
/******************************************************************************
|
||||||
* MinuitMinimizer implementation *
|
* MinuitMinimizer implementation *
|
||||||
******************************************************************************/
|
******************************************************************************/
|
||||||
// constructors ////////////////////////////////////////////////////////////////
|
// constructors ////////////////////////////////////////////////////////////////
|
||||||
MinuitMinimizer::MinuitMinimizer(const Algorithm algorithm)
|
MinuitMinimizer::MinuitMinimizer(const Algorithm algorithm)
|
||||||
|
|||||||
@ -27,7 +27,7 @@
|
|||||||
BEGIN_LATAN_NAMESPACE
|
BEGIN_LATAN_NAMESPACE
|
||||||
|
|
||||||
/******************************************************************************
|
/******************************************************************************
|
||||||
* MinuitMinimizer *
|
* interface to CERN Minuit minimizer (http://www.cern.ch/minuit) *
|
||||||
******************************************************************************/
|
******************************************************************************/
|
||||||
class MinuitMinimizer: public Minimizer
|
class MinuitMinimizer: public Minimizer
|
||||||
{
|
{
|
||||||
|
|||||||
@ -1,374 +0,0 @@
|
|||||||
/*
|
|
||||||
* XYSampleData.cpp, part of LatAnalyze 3
|
|
||||||
*
|
|
||||||
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#include <LatAnalyze/XYSampleData.hpp>
|
|
||||||
#include <LatAnalyze/Math.hpp>
|
|
||||||
#include <LatAnalyze/includes.hpp>
|
|
||||||
|
|
||||||
using namespace std;
|
|
||||||
using namespace Latan;
|
|
||||||
using namespace Math;
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* 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<double>(nDof_);
|
|
||||||
}
|
|
||||||
|
|
||||||
double SampleFitResult::getPValue(const Index s) const
|
|
||||||
{
|
|
||||||
return chi2PValue(getChi2(s), getNDof());
|
|
||||||
}
|
|
||||||
|
|
||||||
const DoubleFunction & SampleFitResult::getModel(const Index s,
|
|
||||||
const Index j) const
|
|
||||||
{
|
|
||||||
return model_[static_cast<unsigned int>(j)][s];
|
|
||||||
}
|
|
||||||
|
|
||||||
const DoubleFunctionSample & SampleFitResult::getModel(
|
|
||||||
const PlaceHolder ph __dumb,
|
|
||||||
const Index j) const
|
|
||||||
{
|
|
||||||
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 *
|
|
||||||
******************************************************************************/
|
|
||||||
// constructors ////////////////////////////////////////////////////////////////
|
|
||||||
XYSampleData::XYSampleData(const Index nData, const Index xDim,
|
|
||||||
const Index yDim, const Index nSample)
|
|
||||||
{
|
|
||||||
resize(nData, xDim, yDim, nSample);
|
|
||||||
}
|
|
||||||
|
|
||||||
// access //////////////////////////////////////////////////////////////////////
|
|
||||||
const XYStatData & XYSampleData::getData(const Index s)
|
|
||||||
{
|
|
||||||
setDataToSample(s);
|
|
||||||
|
|
||||||
return data_;
|
|
||||||
}
|
|
||||||
|
|
||||||
void XYSampleData::resize(const Index nData, const Index xDim,
|
|
||||||
const Index yDim, const Index nSample)
|
|
||||||
{
|
|
||||||
FitInterface::resize(nData, xDim, yDim);
|
|
||||||
x_.resize(nSample);
|
|
||||||
x_.resizeMat(nData, xDim);
|
|
||||||
y_.resize(nSample);
|
|
||||||
y_.resizeMat(nData, yDim);
|
|
||||||
data_.resize(nData, xDim, yDim);
|
|
||||||
isCovarianceInit_ = false;
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::SampleBlock XYSampleData::x(const PlaceHolder ph1 __dumb,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
{
|
|
||||||
isCovarianceInit_ = false;
|
|
||||||
|
|
||||||
return x_.block(0, 0, getNData(), getXDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::ConstSampleBlock XYSampleData::x(const PlaceHolder ph1 __dumb,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return x_.block(0, 0, getNData(), getXDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::SampleBlock XYSampleData::x(const Index i,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
{
|
|
||||||
isCovarianceInit_ = false;
|
|
||||||
|
|
||||||
return x_.block(0, i, getNData(), 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::ConstSampleBlock XYSampleData::x(const Index i,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return x_.block(0, i, getNData(), 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::SampleBlock XYSampleData::x(const PlaceHolder ph1 __dumb,
|
|
||||||
const Index k)
|
|
||||||
{
|
|
||||||
isCovarianceInit_ = false;
|
|
||||||
|
|
||||||
return x_.block(k, 0, 1, getXDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::ConstSampleBlock XYSampleData::x(const PlaceHolder ph1 __dumb,
|
|
||||||
const Index k) const
|
|
||||||
{
|
|
||||||
return x_.block(k, 0, 1, getXDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::SampleBlock XYSampleData::x(const Index i, const Index k)
|
|
||||||
{
|
|
||||||
isCovarianceInit_ = false;
|
|
||||||
|
|
||||||
return x_.block(k, i, 1, 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::ConstSampleBlock XYSampleData::x(const Index i, const Index k)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return x_.block(k, i, 1, 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::SampleBlock XYSampleData::y(const PlaceHolder ph1 __dumb,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
{
|
|
||||||
isCovarianceInit_ = false;
|
|
||||||
|
|
||||||
return y_.block(0, 0, getNData(), getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::ConstSampleBlock XYSampleData::y(const PlaceHolder ph1 __dumb,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return y_.block(0, 0, getNData(), getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::SampleBlock XYSampleData::y(const Index j,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
{
|
|
||||||
isCovarianceInit_ = false;
|
|
||||||
|
|
||||||
return y_.block(0, j, getNData(), 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::ConstSampleBlock XYSampleData::y(const Index j,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return y_.block(0, j, getNData(), 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::SampleBlock XYSampleData::y(const PlaceHolder ph1 __dumb,
|
|
||||||
const Index k)
|
|
||||||
{
|
|
||||||
isCovarianceInit_ = false;
|
|
||||||
|
|
||||||
return y_.block(k, 0, 1, getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::ConstSampleBlock XYSampleData::y(const PlaceHolder ph1 __dumb,
|
|
||||||
const Index k) const
|
|
||||||
{
|
|
||||||
return y_.block(k, 0, 1, getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::SampleBlock XYSampleData::y(const Index j, const Index k)
|
|
||||||
{
|
|
||||||
isCovarianceInit_ = false;
|
|
||||||
|
|
||||||
return y_.block(k, j, 1, 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
XYSampleData::ConstSampleBlock XYSampleData::y(const Index j, const Index k)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return y_.block(k, j, 1, 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
// fit /////////////////////////////////////////////////////////////////////////
|
|
||||||
SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init,
|
|
||||||
const std::vector<const DoubleModel *> &modelVector)
|
|
||||||
{
|
|
||||||
const Index nSample = x_.size();
|
|
||||||
FitResult sampleResult;
|
|
||||||
SampleFitResult result;
|
|
||||||
DVec initBuf = init;
|
|
||||||
|
|
||||||
result.resize(nSample);
|
|
||||||
result.chi2_.resize(nSample);
|
|
||||||
FOR_STAT_ARRAY(x_, s)
|
|
||||||
{
|
|
||||||
// reinit chi^2 for central value only
|
|
||||||
if (s == central)
|
|
||||||
{
|
|
||||||
data_.reinitChi2(true);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
data_.reinitChi2(false);
|
|
||||||
}
|
|
||||||
|
|
||||||
// set data
|
|
||||||
setDataToSample(s);
|
|
||||||
|
|
||||||
// fit
|
|
||||||
sampleResult = data_.fit(minimizer, initBuf, modelVector);
|
|
||||||
if (s == central)
|
|
||||||
{
|
|
||||||
initBuf = sampleResult;
|
|
||||||
}
|
|
||||||
|
|
||||||
// store result
|
|
||||||
result[s] = sampleResult;
|
|
||||||
result.chi2_[s] = sampleResult.getChi2();
|
|
||||||
result.nDof_ = sampleResult.getNDof();
|
|
||||||
result.model_.resize(modelVector.size());
|
|
||||||
for (unsigned int j = 0; j < modelVector.size(); ++j)
|
|
||||||
{
|
|
||||||
result.model_[j].resize(nSample);
|
|
||||||
result.model_[j][s] = sampleResult.getModel(j);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
return result;
|
|
||||||
}
|
|
||||||
|
|
||||||
void XYSampleData::setDataToSample(const Index s)
|
|
||||||
{
|
|
||||||
// compute covariance matrices if necessary
|
|
||||||
if (!isCovarianceInit_)
|
|
||||||
{
|
|
||||||
DMatSample buf1, buf2;
|
|
||||||
|
|
||||||
for (Index i2 = 0; i2 < getXDim(); ++i2)
|
|
||||||
for (Index i1 = 0; i1 < getXDim(); ++i1)
|
|
||||||
{
|
|
||||||
buf1 = x(i1);
|
|
||||||
buf2 = x(i2);
|
|
||||||
data_.xxVar(i1, i2) = buf1.covarianceMatrix(buf2);
|
|
||||||
}
|
|
||||||
for (Index j2 = 0; j2 < getYDim(); ++j2)
|
|
||||||
for (Index j1 = 0; j1 < getYDim(); ++j1)
|
|
||||||
{
|
|
||||||
buf1 = y(j1);
|
|
||||||
buf2 = y(j2);
|
|
||||||
data_.yyVar(j1, j2) = buf1.covarianceMatrix(buf2);
|
|
||||||
}
|
|
||||||
for (Index i = 0; i < getXDim(); ++i)
|
|
||||||
for (Index j = 0; j < getYDim(); ++j)
|
|
||||||
{
|
|
||||||
buf1 = y(j);
|
|
||||||
buf2 = x(i);
|
|
||||||
data_.yxVar(j, i) = buf1.covarianceMatrix(buf2);
|
|
||||||
}
|
|
||||||
isCovarianceInit_ = true;
|
|
||||||
}
|
|
||||||
|
|
||||||
// copy interface to sample data
|
|
||||||
data_.setFitInterface(*this);
|
|
||||||
|
|
||||||
// set data
|
|
||||||
data_.x() = x_[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;
|
|
||||||
}
|
|
||||||
@ -1,136 +0,0 @@
|
|||||||
/*
|
|
||||||
* XYSampleData.hpp, part of LatAnalyze 3
|
|
||||||
*
|
|
||||||
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#ifndef Latan_XYSampleData_hpp_
|
|
||||||
#define Latan_XYSampleData_hpp_
|
|
||||||
|
|
||||||
#include <LatAnalyze/Global.hpp>
|
|
||||||
#include <LatAnalyze/FitInterface.hpp>
|
|
||||||
#include <LatAnalyze/Function.hpp>
|
|
||||||
#include <LatAnalyze/MatSample.hpp>
|
|
||||||
#include <LatAnalyze/Minimizer.hpp>
|
|
||||||
#include <LatAnalyze/XYStatData.hpp>
|
|
||||||
|
|
||||||
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<DoubleFunctionSample> model_;
|
|
||||||
};
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* XYSampleData *
|
|
||||||
******************************************************************************
|
|
||||||
* index convention: i: X, j: Y, k: data
|
|
||||||
*/
|
|
||||||
class XYSampleData: public FitInterface
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
typedef DMatSample::Block SampleBlock;
|
|
||||||
typedef DMatSample::ConstBlock ConstSampleBlock;
|
|
||||||
public:
|
|
||||||
// constructors
|
|
||||||
XYSampleData(void) = default;
|
|
||||||
XYSampleData(const Index nData, const Index xDim, const Index yDim,
|
|
||||||
const Index nSample);
|
|
||||||
// destructor
|
|
||||||
virtual ~XYSampleData(void) = default;
|
|
||||||
// access
|
|
||||||
const XYStatData & getData(const Index s = central);
|
|
||||||
void resize(const Index nData, const Index xDim,
|
|
||||||
const Index yDim, const Index nSample);
|
|
||||||
SampleBlock x(const PlaceHolder ph1 = _, const PlaceHolder ph2 = _);
|
|
||||||
ConstSampleBlock x(const PlaceHolder ph1 = _,
|
|
||||||
const PlaceHolder ph2 = _) const;
|
|
||||||
SampleBlock x(const Index i, const PlaceHolder ph2 = _);
|
|
||||||
ConstSampleBlock x(const Index i, const PlaceHolder ph2 = _) const;
|
|
||||||
SampleBlock x(const PlaceHolder ph1, const Index k);
|
|
||||||
ConstSampleBlock x(const PlaceHolder ph1, const Index k) const;
|
|
||||||
SampleBlock x(const Index i, const Index k);
|
|
||||||
ConstSampleBlock x(const Index i, const Index k) const;
|
|
||||||
SampleBlock y(const PlaceHolder ph1 = _, const PlaceHolder ph2 = _);
|
|
||||||
ConstSampleBlock y(const PlaceHolder ph1 = _,
|
|
||||||
const PlaceHolder ph2 = _) const;
|
|
||||||
SampleBlock y(const Index j, const PlaceHolder ph2 = _);
|
|
||||||
ConstSampleBlock y(const Index j, const PlaceHolder ph2 = _) const;
|
|
||||||
SampleBlock y(const PlaceHolder ph1, const Index k);
|
|
||||||
ConstSampleBlock y(const PlaceHolder ph1, const Index k) const;
|
|
||||||
SampleBlock y(const Index j, const Index k);
|
|
||||||
ConstSampleBlock y(const Index j, const Index k) const;
|
|
||||||
// fit
|
|
||||||
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
|
|
||||||
const std::vector<const DoubleModel *> &modelVector);
|
|
||||||
template <typename... Mods>
|
|
||||||
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
|
|
||||||
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:
|
|
||||||
void setDataToSample(const Index s);
|
|
||||||
private:
|
|
||||||
bool isCovarianceInit_;
|
|
||||||
DMatSample x_, y_;
|
|
||||||
XYStatData data_;
|
|
||||||
};
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* XYSampleData template implementation *
|
|
||||||
******************************************************************************/
|
|
||||||
template <typename... Ts>
|
|
||||||
SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init,
|
|
||||||
const DoubleModel &model, const Ts... models)
|
|
||||||
{
|
|
||||||
static_assert(static_or<std::is_assignable<DoubleModel &, Ts>::value...>::value,
|
|
||||||
"model arguments are not compatible with DoubleModel &");
|
|
||||||
|
|
||||||
std::vector<const DoubleModel *> modelVector{&model, &models...};
|
|
||||||
|
|
||||||
return fit(minimizer, init, modelVector);
|
|
||||||
}
|
|
||||||
|
|
||||||
END_LATAN_NAMESPACE
|
|
||||||
|
|
||||||
#endif // Latan_XYSampleData_hpp_
|
|
||||||
@ -1,315 +0,0 @@
|
|||||||
/*
|
|
||||||
* XYStatData.cpp, part of LatAnalyze 3
|
|
||||||
*
|
|
||||||
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#include <LatAnalyze/XYStatData.hpp>
|
|
||||||
#include <LatAnalyze/includes.hpp>
|
|
||||||
#include <LatAnalyze/Math.hpp>
|
|
||||||
|
|
||||||
using namespace std;
|
|
||||||
using namespace Latan;
|
|
||||||
using namespace Math;
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* FitResult implementation *
|
|
||||||
******************************************************************************/
|
|
||||||
// access //////////////////////////////////////////////////////////////////////
|
|
||||||
double FitResult::getChi2(void) const
|
|
||||||
{
|
|
||||||
return chi2_;
|
|
||||||
}
|
|
||||||
|
|
||||||
double FitResult::getChi2PerDof(void) const
|
|
||||||
{
|
|
||||||
return chi2_/getNDof();
|
|
||||||
}
|
|
||||||
|
|
||||||
double FitResult::getNDof(void) const
|
|
||||||
{
|
|
||||||
return static_cast<double>(nDof_);
|
|
||||||
}
|
|
||||||
|
|
||||||
double FitResult::getPValue(void) const
|
|
||||||
{
|
|
||||||
return chi2PValue(getChi2(), getNDof());;
|
|
||||||
}
|
|
||||||
|
|
||||||
const DoubleFunction & FitResult::getModel(const Index j) const
|
|
||||||
{
|
|
||||||
return model_[static_cast<unsigned int>(j)];
|
|
||||||
}
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* XYStatData implementation *
|
|
||||||
******************************************************************************/
|
|
||||||
// constructor /////////////////////////////////////////////////////////////////
|
|
||||||
XYStatData::XYStatData(void)
|
|
||||||
: chi2_(*this)
|
|
||||||
{}
|
|
||||||
|
|
||||||
XYStatData::XYStatData(const Index nData, const Index xDim, const Index yDim)
|
|
||||||
: XYStatData()
|
|
||||||
{
|
|
||||||
resize(nData, xDim, yDim);
|
|
||||||
}
|
|
||||||
|
|
||||||
// access //////////////////////////////////////////////////////////////////////
|
|
||||||
void XYStatData::resize(const Index nData, const Index xDim, const Index yDim)
|
|
||||||
{
|
|
||||||
FitInterface::resize(nData, xDim, yDim);
|
|
||||||
x_.resize(nData, xDim);
|
|
||||||
y_.resize(nData, yDim);
|
|
||||||
var_[xx].resize(xDim, xDim);
|
|
||||||
var_[yy].resize(yDim, yDim);
|
|
||||||
var_[yx].resize(yDim, xDim);
|
|
||||||
FOR_MAT(var_[xx], i1, i2)
|
|
||||||
{
|
|
||||||
var_[xx](i1, i2).resize(nData, nData);
|
|
||||||
}
|
|
||||||
FOR_MAT(var_[yy], j1, j2)
|
|
||||||
{
|
|
||||||
var_[yy](j1, j2).resize(nData, nData);
|
|
||||||
}
|
|
||||||
FOR_MAT(var_[yx], j, i)
|
|
||||||
{
|
|
||||||
var_[yx](j, i).resize(nData, nData);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void XYStatData::reinitChi2(const bool doReinit)
|
|
||||||
{
|
|
||||||
reinitChi2_ = doReinit;
|
|
||||||
}
|
|
||||||
|
|
||||||
Block<MatBase<double>> XYStatData::x(const PlaceHolder ph1 __dumb,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
{
|
|
||||||
return x_.block(0, 0, getNData(), getXDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
ConstBlock<MatBase<double>> XYStatData::x(const PlaceHolder ph1 __dumb,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return x_.block(0, 0, getNData(), getXDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
Block<MatBase<double>> XYStatData::x(const Index i,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
{
|
|
||||||
return x_.block(0, i, getNData(), 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
ConstBlock<MatBase<double>> XYStatData::x(const Index i,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return x_.block(0, i, getNData(), 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
Block<MatBase<double>> XYStatData::x(const PlaceHolder ph1 __dumb,
|
|
||||||
const Index k)
|
|
||||||
{
|
|
||||||
return x_.block(k, 0, 1, getXDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
ConstBlock<MatBase<double>> XYStatData::x(const PlaceHolder ph1 __dumb,
|
|
||||||
const Index k) const
|
|
||||||
{
|
|
||||||
return x_.block(k, 0, 1, getXDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
double & XYStatData::x(const Index i, const Index k)
|
|
||||||
{
|
|
||||||
return x_(k, i);
|
|
||||||
}
|
|
||||||
|
|
||||||
const double & XYStatData::x(const Index i, const Index k) const
|
|
||||||
{
|
|
||||||
return x_(k, i);
|
|
||||||
}
|
|
||||||
|
|
||||||
Block<MatBase<double>> XYStatData::y(const PlaceHolder ph1 __dumb,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
{
|
|
||||||
return y_.block(0, 0, getNData(), getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
ConstBlock<MatBase<double>> XYStatData::y(const PlaceHolder ph1 __dumb,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return y_.block(0, 0, getNData(), getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
Block<MatBase<double>> XYStatData::y(const Index j,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
{
|
|
||||||
return y_.block(0, j, getNData(), 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
ConstBlock<MatBase<double>> XYStatData::y(const Index j,
|
|
||||||
const PlaceHolder ph2 __dumb)
|
|
||||||
const
|
|
||||||
{
|
|
||||||
return y_.block(0, j, getNData(), 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
Block<MatBase<double>> XYStatData::y(const PlaceHolder ph1 __dumb, const Index k)
|
|
||||||
{
|
|
||||||
return y_.block(k, 0, 1, getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
ConstBlock<MatBase<double>> XYStatData::y(const PlaceHolder ph1 __dumb,
|
|
||||||
const Index k) const
|
|
||||||
{
|
|
||||||
return y_.block(k, 0, 1, getYDim());
|
|
||||||
}
|
|
||||||
|
|
||||||
double & XYStatData::y(const Index j, const Index k)
|
|
||||||
{
|
|
||||||
return y_(k, j);
|
|
||||||
}
|
|
||||||
|
|
||||||
const double & XYStatData::y(const Index j, const Index k) const
|
|
||||||
{
|
|
||||||
return y_(k, j);
|
|
||||||
}
|
|
||||||
|
|
||||||
#define FULL_BLOCK(m) (m).block(0, 0, (m).rows(), (m).cols())
|
|
||||||
|
|
||||||
Block<MatBase<double>> XYStatData::xxVar(const Index i1, const Index i2)
|
|
||||||
{
|
|
||||||
return FULL_BLOCK(var_[xx](i1, i2));
|
|
||||||
}
|
|
||||||
|
|
||||||
ConstBlock<MatBase<double>> XYStatData::xxVar(const Index i1,
|
|
||||||
const Index i2) const
|
|
||||||
{
|
|
||||||
return FULL_BLOCK(var_[xx](i1, i2));
|
|
||||||
}
|
|
||||||
|
|
||||||
Block<MatBase<double>> XYStatData::yyVar(const Index j1, const Index j2)
|
|
||||||
{
|
|
||||||
return FULL_BLOCK(var_[yy](j1, j2));
|
|
||||||
}
|
|
||||||
|
|
||||||
ConstBlock<MatBase<double>> XYStatData::yyVar(const Index j1,
|
|
||||||
const Index j2) const
|
|
||||||
{
|
|
||||||
return FULL_BLOCK(var_[yy](j1, j2));
|
|
||||||
}
|
|
||||||
|
|
||||||
Block<MatBase<double>> XYStatData::yxVar(const Index j, const Index i)
|
|
||||||
{
|
|
||||||
return FULL_BLOCK(var_[yx](j, i));
|
|
||||||
}
|
|
||||||
|
|
||||||
ConstBlock<MatBase<double>> XYStatData::yxVar(const Index j,
|
|
||||||
const Index i) const
|
|
||||||
{
|
|
||||||
return FULL_BLOCK(var_[yx](j, i));
|
|
||||||
}
|
|
||||||
|
|
||||||
// fit /////////////////////////////////////////////////////////////////////////
|
|
||||||
FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
|
|
||||||
const vector<const DoubleModel *> &modelVector)
|
|
||||||
{
|
|
||||||
// initialization
|
|
||||||
chi2_.setModel(modelVector);
|
|
||||||
if (reinitChi2_)
|
|
||||||
{
|
|
||||||
chi2_.requestInit();
|
|
||||||
}
|
|
||||||
// initial parameters
|
|
||||||
const Index nPoint = getNFitPoint();
|
|
||||||
DVec fullInit = init;
|
|
||||||
Index is = 0, kf = 0;
|
|
||||||
|
|
||||||
fullInit.conservativeResize(chi2_.getNArg());
|
|
||||||
for (Index i = 0; i < getXDim(); ++i)
|
|
||||||
{
|
|
||||||
if (!isXExact(i))
|
|
||||||
{
|
|
||||||
kf = 0;
|
|
||||||
for (Index k = 0; k < getNData(); ++k)
|
|
||||||
{
|
|
||||||
if (isFitPoint(k))
|
|
||||||
{
|
|
||||||
fullInit(chi2_.getNPar() + nPoint*is + kf) = x(i, k);
|
|
||||||
kf++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
is++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
minimizer.setInit(fullInit);
|
|
||||||
|
|
||||||
// fit
|
|
||||||
DoubleFunction chi2 = chi2_.makeFunction(false);
|
|
||||||
FitResult result;
|
|
||||||
|
|
||||||
result = minimizer(chi2);
|
|
||||||
result.chi2_ = chi2(result);
|
|
||||||
result.nDof_ = chi2_.getNDof();
|
|
||||||
result.model_.resize(modelVector.size());
|
|
||||||
for (unsigned int j = 0; j < modelVector.size(); ++j)
|
|
||||||
{
|
|
||||||
result.model_[j] = modelVector[j]->fixPar(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;
|
|
||||||
}
|
|
||||||
@ -1,150 +0,0 @@
|
|||||||
/*
|
|
||||||
* XYData.hpp, part of LatAnalyze 3
|
|
||||||
*
|
|
||||||
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
|
|
||||||
*/
|
|
||||||
|
|
||||||
#ifndef Latan_XYData_hpp_
|
|
||||||
#define Latan_XYData_hpp_
|
|
||||||
|
|
||||||
#include <LatAnalyze/Global.hpp>
|
|
||||||
#include <LatAnalyze/Chi2Function.hpp>
|
|
||||||
#include <LatAnalyze/FitInterface.hpp>
|
|
||||||
#include <LatAnalyze/Function.hpp>
|
|
||||||
#include <LatAnalyze/Mat.hpp>
|
|
||||||
#include <LatAnalyze/Minimizer.hpp>
|
|
||||||
#include <LatAnalyze/Model.hpp>
|
|
||||||
#include <vector>
|
|
||||||
|
|
||||||
BEGIN_LATAN_NAMESPACE
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* object for fit result *
|
|
||||||
******************************************************************************/
|
|
||||||
class FitResult: public DVec
|
|
||||||
{
|
|
||||||
friend class XYStatData;
|
|
||||||
friend class SampleFitResult;
|
|
||||||
public:
|
|
||||||
// constructors
|
|
||||||
FitResult(void) = default;
|
|
||||||
EIGEN_EXPR_CTOR(FitResult, FitResult, Base, MatExpr)
|
|
||||||
// destructor
|
|
||||||
virtual ~FitResult(void) = default;
|
|
||||||
// access
|
|
||||||
double getChi2(void) const;
|
|
||||||
double getChi2PerDof(void) const;
|
|
||||||
double getNDof(void) const;
|
|
||||||
double getPValue(void) const;
|
|
||||||
const DoubleFunction & getModel(const Index j = 0) const;
|
|
||||||
private:
|
|
||||||
double chi2_{0.0};
|
|
||||||
Index nDof_{0};
|
|
||||||
std::vector<DoubleFunction> model_;
|
|
||||||
};
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* object for X vs. Y statistical data *
|
|
||||||
******************************************************************************
|
|
||||||
* index convention: i: X, j: Y, k: data
|
|
||||||
*/
|
|
||||||
class XYStatData: public FitInterface
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
enum
|
|
||||||
{
|
|
||||||
xx = 0,
|
|
||||||
yy = 1,
|
|
||||||
yx = 2
|
|
||||||
};
|
|
||||||
public:
|
|
||||||
// constructors
|
|
||||||
XYStatData(void);
|
|
||||||
XYStatData(const Index nData, const Index nXDim, const Index nYDim);
|
|
||||||
// destructor
|
|
||||||
virtual ~XYStatData(void) = default;
|
|
||||||
// access
|
|
||||||
void resize(const Index nData, const Index xDim,
|
|
||||||
const Index yDim);
|
|
||||||
void reinitChi2(const bool doReinit = true);
|
|
||||||
Block<MatBase<double>> x(const PlaceHolder ph1 = _,
|
|
||||||
const PlaceHolder ph2 = _);
|
|
||||||
ConstBlock<MatBase<double>> x(const PlaceHolder ph1 = _,
|
|
||||||
const PlaceHolder ph2 = _) const;
|
|
||||||
Block<MatBase<double>> x(const Index i, const PlaceHolder ph2 = _);
|
|
||||||
ConstBlock<MatBase<double>> x(const Index i,
|
|
||||||
const PlaceHolder ph2 = _) const;
|
|
||||||
Block<MatBase<double>> x(const PlaceHolder ph1, const Index k);
|
|
||||||
ConstBlock<MatBase<double>> x(const PlaceHolder ph1,
|
|
||||||
const Index k) const;
|
|
||||||
double & x(const Index i, const Index k);
|
|
||||||
const double & x(const Index i, const Index k) const;
|
|
||||||
Block<MatBase<double>> y(const PlaceHolder ph1 = _,
|
|
||||||
const PlaceHolder ph2 = _);
|
|
||||||
ConstBlock<MatBase<double>> y(const PlaceHolder ph1 = _,
|
|
||||||
const PlaceHolder ph2 = _) const;
|
|
||||||
Block<MatBase<double>> y(const Index i, const PlaceHolder ph2 = _);
|
|
||||||
ConstBlock<MatBase<double>> y(const Index i,
|
|
||||||
const PlaceHolder ph2 = _) const;
|
|
||||||
Block<MatBase<double>> y(const PlaceHolder ph1, const Index k);
|
|
||||||
ConstBlock<MatBase<double>> y(const PlaceHolder ph1,
|
|
||||||
const Index k) const;
|
|
||||||
double & y(const Index i, const Index k);
|
|
||||||
const double & y(const Index i, const Index k) const;
|
|
||||||
Block<MatBase<double>> xxVar(const Index i1, const Index i2);
|
|
||||||
ConstBlock<MatBase<double>> xxVar(const Index i1, const Index i2) const;
|
|
||||||
Block<MatBase<double>> yyVar(const Index j1, const Index j2);
|
|
||||||
ConstBlock<MatBase<double>> yyVar(const Index j1, const Index j2) const;
|
|
||||||
Block<MatBase<double>> yxVar(const Index j, const Index i);
|
|
||||||
ConstBlock<MatBase<double>> yxVar(const Index j, const Index i) const;
|
|
||||||
// fit
|
|
||||||
FitResult fit(Minimizer &minimizer, const DVec &init,
|
|
||||||
const std::vector<const DoubleModel *> &modelVector);
|
|
||||||
template <typename... Ts>
|
|
||||||
FitResult fit(Minimizer &minimizer, const DVec &init,
|
|
||||||
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:
|
|
||||||
DMat x_, y_;
|
|
||||||
Mat<DMat> var_[3];
|
|
||||||
IVec isXExact_, isFitPoint_;
|
|
||||||
Chi2Function chi2_;
|
|
||||||
bool reinitChi2_{true};
|
|
||||||
};
|
|
||||||
|
|
||||||
/******************************************************************************
|
|
||||||
* XYStatData template implementation *
|
|
||||||
******************************************************************************/
|
|
||||||
template <typename... Ts>
|
|
||||||
FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
|
|
||||||
const DoubleModel &model, const Ts... models)
|
|
||||||
{
|
|
||||||
static_assert(static_or<std::is_assignable<DoubleModel &, Ts>::value...>::value,
|
|
||||||
"model arguments are not compatible with DoubleModel &");
|
|
||||||
|
|
||||||
std::vector<const DoubleModel *> modelVector{&model, &models...};
|
|
||||||
|
|
||||||
return fit(minimizer, init, modelVector);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
END_LATAN_NAMESPACE
|
|
||||||
|
|
||||||
#endif // Latan_XYData_hpp_
|
|
||||||
Loading…
x
Reference in New Issue
Block a user