mirror of
https://github.com/aportelli/LatAnalyze.git
synced 2025-06-16 22:47:04 +01:00
reintegration of LatCore & folder restructuration
This commit is contained in:
152
lib/Statistics/Dataset.hpp
Normal file
152
lib/Statistics/Dataset.hpp
Normal file
@ -0,0 +1,152 @@
|
||||
/*
|
||||
* Dataset.hpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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_Dataset_hpp_
|
||||
#define Latan_Dataset_hpp_
|
||||
|
||||
#include <LatAnalyze/Global.hpp>
|
||||
#include <LatAnalyze/Io/File.hpp>
|
||||
#include <LatAnalyze/Statistics/StatArray.hpp>
|
||||
|
||||
BEGIN_LATAN_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* Dataset class *
|
||||
******************************************************************************/
|
||||
template <typename T>
|
||||
class Dataset: public StatArray<T>
|
||||
{
|
||||
public:
|
||||
typedef std::random_device::result_type SeedType;
|
||||
public:
|
||||
// constructors
|
||||
Dataset(void) = default;
|
||||
Dataset(const Index size);
|
||||
EIGEN_EXPR_CTOR(Dataset, Dataset<T>, StatArray<T>, ArrayExpr)
|
||||
// destructor
|
||||
virtual ~Dataset(void) = default;
|
||||
// IO
|
||||
template <typename FileType>
|
||||
void load(const std::string &listFileName, const std::string &dataName);
|
||||
// resampling
|
||||
Sample<T> bootstrapMean(const Index nSample, const SeedType seed);
|
||||
Sample<T> bootstrapMean(const Index nSample);
|
||||
void dumpBootstrapSeq(std::ostream &out, const Index nSample,
|
||||
const SeedType seed);
|
||||
private:
|
||||
// mean from pointer vector for resampling
|
||||
void ptVectorMean(T &m, const std::vector<const T *> &v);
|
||||
};
|
||||
|
||||
/******************************************************************************
|
||||
* Dataset template implementation *
|
||||
******************************************************************************/
|
||||
// constructors ////////////////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
Dataset<T>::Dataset(const Index size)
|
||||
: StatArray<T>(size)
|
||||
{}
|
||||
|
||||
// IO //////////////////////////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
template <typename FileType>
|
||||
void Dataset<T>::load(const std::string &listFileName,
|
||||
const std::string &dataName)
|
||||
{
|
||||
FileType file;
|
||||
std::vector<std::string> dataFileName;
|
||||
|
||||
dataFileName = readManifest(listFileName);
|
||||
this->resize(dataFileName.size());
|
||||
for (Index i = 0; i < static_cast<Index>(dataFileName.size()); ++i)
|
||||
{
|
||||
file.open(dataFileName[i], File::Mode::read);
|
||||
(*this)[i] = file.template read<T>(dataName);
|
||||
file.close();
|
||||
}
|
||||
}
|
||||
|
||||
// resampling //////////////////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
Sample<T> Dataset<T>::bootstrapMean(const Index nSample, const SeedType seed)
|
||||
{
|
||||
std::vector<const T *> data(this->size());
|
||||
Sample<T> s(nSample);
|
||||
std::mt19937 gen(seed);
|
||||
std::uniform_int_distribution<Index> dis(0, this->size() - 1);
|
||||
|
||||
for (unsigned int j = 0; j < this->size(); ++j)
|
||||
{
|
||||
data[j] = &((*this)[static_cast<Index>(j)]);
|
||||
}
|
||||
ptVectorMean(s[central], data);
|
||||
for (Index i = 0; i < nSample; ++i)
|
||||
{
|
||||
for (unsigned int j = 0; j < this->size(); ++j)
|
||||
{
|
||||
data[j] = &((*this)[dis(gen)]);
|
||||
}
|
||||
ptVectorMean(s[i], data);
|
||||
}
|
||||
|
||||
return s;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
Sample<T> Dataset<T>::bootstrapMean(const Index nSample)
|
||||
{
|
||||
std::random_device rd;
|
||||
|
||||
return bootstrapMean(nSample, rd());
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void Dataset<T>::dumpBootstrapSeq(std::ostream &out, const Index nSample,
|
||||
const SeedType seed)
|
||||
{
|
||||
std::mt19937 gen(seed);
|
||||
std::uniform_int_distribution<Index> dis(0, this->size() - 1);
|
||||
|
||||
for (Index i = 0; i < nSample; ++i)
|
||||
{
|
||||
for (unsigned int j = 0; j < this->size(); ++j)
|
||||
{
|
||||
out << dis(gen) << " " << std::endl;
|
||||
}
|
||||
out << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void Dataset<T>::ptVectorMean(T &m, const std::vector<const T *> &v)
|
||||
{
|
||||
if (v.size())
|
||||
{
|
||||
m = *(v[0]);
|
||||
for (unsigned int i = 1; i < v.size(); ++i)
|
||||
{
|
||||
m += *(v[i]);
|
||||
}
|
||||
m /= static_cast<double>(v.size());
|
||||
}
|
||||
}
|
||||
|
||||
END_LATAN_NAMESPACE
|
||||
|
||||
#endif // Latan_Dataset_hpp_
|
||||
794
lib/Statistics/FitInterface.cpp
Normal file
794
lib/Statistics/FitInterface.cpp
Normal file
@ -0,0 +1,794 @@
|
||||
/*
|
||||
* FitInterface.cpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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/Statistics/FitInterface.hpp>
|
||||
#include <LatAnalyze/includes.hpp>
|
||||
|
||||
using namespace std;
|
||||
using namespace Latan;
|
||||
|
||||
/******************************************************************************
|
||||
* FitInterface implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
FitInterface::FitInterface(void)
|
||||
: xName_("x")
|
||||
, yName_("y")
|
||||
{}
|
||||
|
||||
// copy object (not as a constructor to be accessed from derived class) ////////
|
||||
void FitInterface::copyInterface(const FitInterface &d)
|
||||
{
|
||||
*this = d;
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
// add dimensions //////////////////////////////////////////////////////////////
|
||||
void FitInterface::addXDim(const Index nData, const string name,
|
||||
const bool isExact)
|
||||
{
|
||||
if (getYSize() != 0)
|
||||
{
|
||||
LATAN_ERROR(Logic, "cannot add an X dimension if fit data is "
|
||||
"not empty");
|
||||
}
|
||||
else
|
||||
{
|
||||
xSize_.push_back(nData);
|
||||
xIsExact_.push_back(isExact);
|
||||
maxDataIndex_ *= nData;
|
||||
createXData(name, nData);
|
||||
scheduleLayoutInit();
|
||||
scheduleDataCoordInit();
|
||||
if (!name.empty())
|
||||
{
|
||||
xName().setName(getNXDim() - 1, name);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void FitInterface::addYDim(const string name)
|
||||
{
|
||||
yDataIndex_.push_back(map<Index, bool>());
|
||||
createYData(name);
|
||||
scheduleLayoutInit();
|
||||
if (!name.empty())
|
||||
{
|
||||
yName().setName(getNYDim() - 1, name);
|
||||
}
|
||||
}
|
||||
|
||||
// access //////////////////////////////////////////////////////////////////////
|
||||
Index FitInterface::getNXDim(void) const
|
||||
{
|
||||
return xSize_.size();
|
||||
}
|
||||
|
||||
Index FitInterface::getNYDim(void) const
|
||||
{
|
||||
return yDataIndex_.size();
|
||||
}
|
||||
|
||||
Index FitInterface::getXSize(void) const
|
||||
{
|
||||
Index size = 0;
|
||||
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
{
|
||||
size += getXSize(i);
|
||||
}
|
||||
|
||||
return size;
|
||||
}
|
||||
|
||||
Index FitInterface::getXSize(const Index i) const
|
||||
{
|
||||
checkXDim(i);
|
||||
|
||||
return xSize_[i];
|
||||
}
|
||||
|
||||
Index FitInterface::getYSize(void) const
|
||||
{
|
||||
Index size = 0;
|
||||
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
size += getYSize(j);
|
||||
}
|
||||
|
||||
return size;
|
||||
}
|
||||
|
||||
Index FitInterface::getYSize(const Index j) const
|
||||
{
|
||||
checkYDim(j);
|
||||
|
||||
return static_cast<Index>(yDataIndex_[j].size());
|
||||
}
|
||||
|
||||
Index FitInterface::getXFitSize(void) const
|
||||
{
|
||||
Index size = 0;
|
||||
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
{
|
||||
size += getXFitSize(i);
|
||||
}
|
||||
|
||||
return size;
|
||||
}
|
||||
|
||||
Index FitInterface::getXFitSize(const Index i) const
|
||||
{
|
||||
set<Index> fitCoord;
|
||||
vector<Index> v;
|
||||
|
||||
checkXDim(i);
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
for (auto &p: yDataIndex_[j])
|
||||
{
|
||||
if (p.second)
|
||||
{
|
||||
v = dataCoord(p.first);
|
||||
fitCoord.insert(v[i]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return fitCoord.size();
|
||||
}
|
||||
|
||||
Index FitInterface::getYFitSize(void) const
|
||||
{
|
||||
Index size = 0;
|
||||
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
size += getYFitSize(j);
|
||||
}
|
||||
|
||||
return size;
|
||||
}
|
||||
|
||||
Index FitInterface::getYFitSize(const Index j) const
|
||||
{
|
||||
Index size;
|
||||
auto pred = [](const pair<Index, bool> &p)
|
||||
{
|
||||
return p.second;
|
||||
};
|
||||
|
||||
checkYDim(j);
|
||||
size = count_if(yDataIndex_[j].begin(), yDataIndex_[j].end(), pred);
|
||||
|
||||
return size;
|
||||
}
|
||||
|
||||
Index FitInterface::getMaxDataIndex(void) const
|
||||
{
|
||||
return maxDataIndex_;
|
||||
}
|
||||
|
||||
const set<Index> & FitInterface::getDataIndexSet(void) const
|
||||
{
|
||||
return dataIndexSet_;
|
||||
}
|
||||
|
||||
double FitInterface::getSvdTolerance(void) const
|
||||
{
|
||||
return svdTol_;
|
||||
}
|
||||
|
||||
void FitInterface::setSvdTolerance(const double &tol)
|
||||
{
|
||||
svdTol_ = tol;
|
||||
scheduleLayoutInit();
|
||||
}
|
||||
|
||||
VarName & FitInterface::xName(void)
|
||||
{
|
||||
return xName_;
|
||||
}
|
||||
|
||||
const VarName & FitInterface::xName(void) const
|
||||
{
|
||||
return xName_;
|
||||
}
|
||||
|
||||
VarName & FitInterface::yName(void)
|
||||
{
|
||||
return yName_;
|
||||
}
|
||||
|
||||
const VarName & FitInterface::yName(void) const
|
||||
{
|
||||
return yName_;
|
||||
}
|
||||
|
||||
// Y dimension index helper ////////////////////////////////////////////////////
|
||||
Index FitInterface::dataIndex(const vector<Index> &v) const
|
||||
{
|
||||
Index k, n = v.size();
|
||||
|
||||
checkDataCoord(v);
|
||||
k = xSize_[1]*v[0];
|
||||
for (unsigned int d = 1; d < n-1; ++d)
|
||||
{
|
||||
k = xSize_[d+1]*(v[d] + k);
|
||||
}
|
||||
k += v[n-1];
|
||||
|
||||
return k;
|
||||
}
|
||||
|
||||
const vector<Index> & FitInterface::dataCoord(const Index k) const
|
||||
{
|
||||
checkDataIndex(k);
|
||||
|
||||
updateDataCoord();
|
||||
|
||||
return dataCoord_.at(k);
|
||||
}
|
||||
|
||||
// enable fit points ///////////////////////////////////////////////////////////
|
||||
void FitInterface::fitPoint(const bool isFitPoint, const Index k, const Index j)
|
||||
{
|
||||
checkPoint(k, j);
|
||||
yDataIndex_[j][k] = isFitPoint;
|
||||
scheduleLayoutInit();
|
||||
}
|
||||
|
||||
// variance interface //////////////////////////////////////////////////////////
|
||||
void FitInterface::assumeXExact(const bool isExact, const Index i)
|
||||
{
|
||||
checkXDim(i);
|
||||
xIsExact_[i] = isExact;
|
||||
scheduleLayoutInit();
|
||||
}
|
||||
|
||||
void FitInterface::addCorr(set<array<Index, 4>> &s, const bool isCorr,
|
||||
const array<Index, 4> &c)
|
||||
{
|
||||
if (isCorr)
|
||||
{
|
||||
s.insert(c);
|
||||
}
|
||||
else
|
||||
{
|
||||
auto it = s.find(c);
|
||||
|
||||
if (it != s.end())
|
||||
{
|
||||
s.erase(it);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void FitInterface::assumeXXCorrelated(const bool isCorr, const Index r1,
|
||||
const Index i1, const Index r2,
|
||||
const Index i2)
|
||||
{
|
||||
array<Index, 4> c{{r1, i1, r2, i2}};
|
||||
|
||||
checkXIndex(r1, i1);
|
||||
checkXIndex(r2, i2);
|
||||
if ((i1 != i2) or (r1 != r2))
|
||||
{
|
||||
addCorr(xxCorr_, isCorr, c);
|
||||
}
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
void FitInterface::assumeXXCorrelated(const bool isCorr, const Index i1,
|
||||
const Index i2)
|
||||
{
|
||||
for (Index r1 = 0; r1 < getXSize(i1); ++r1)
|
||||
for (Index r2 = 0; r2 < getXSize(i2); ++r2)
|
||||
{
|
||||
assumeXXCorrelated(isCorr, r1, i1, r2, i2);
|
||||
}
|
||||
}
|
||||
|
||||
void FitInterface::assumeYYCorrelated(const bool isCorr, const Index k1,
|
||||
const Index j1, const Index k2,
|
||||
const Index j2)
|
||||
{
|
||||
array<Index, 4> c{{k1, j1, k2, j2}};
|
||||
|
||||
checkPoint(k1, j1);
|
||||
checkPoint(k2, j2);
|
||||
if ((j1 != j2) or (k1 != k2))
|
||||
{
|
||||
addCorr(yyCorr_, isCorr, c);
|
||||
}
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
void FitInterface::assumeYYCorrelated(const bool isCorr, const Index j1,
|
||||
const Index j2)
|
||||
{
|
||||
checkYDim(j1);
|
||||
checkYDim(j2);
|
||||
for (auto &p1: yDataIndex_[j1])
|
||||
for (auto &p2: yDataIndex_[j2])
|
||||
{
|
||||
assumeYYCorrelated(isCorr, p1.first, j1, p2.first, j2);
|
||||
}
|
||||
}
|
||||
|
||||
void FitInterface::assumeXYCorrelated(const bool isCorr, const Index r,
|
||||
const Index i, const Index k,
|
||||
const Index j)
|
||||
{
|
||||
array<Index, 4> c{{r, i, k, j}};
|
||||
|
||||
checkXIndex(r, i);
|
||||
checkPoint(k, j);
|
||||
addCorr(xyCorr_, isCorr, c);
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
void FitInterface::assumeXYCorrelated(const bool isCorr, const Index i,
|
||||
const Index j)
|
||||
{
|
||||
checkYDim(j);
|
||||
for (Index r = 0; r < getXSize(i); ++r)
|
||||
for (auto &p: yDataIndex_[j])
|
||||
{
|
||||
assumeXYCorrelated(isCorr, r, i, p.first, j);
|
||||
}
|
||||
}
|
||||
|
||||
// tests ///////////////////////////////////////////////////////////////////////
|
||||
bool FitInterface::pointExists(const Index k) const
|
||||
{
|
||||
bool isUsed = false;
|
||||
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
isUsed = isUsed or pointExists(k, j);
|
||||
}
|
||||
|
||||
return isUsed;
|
||||
}
|
||||
|
||||
bool FitInterface::pointExists(const Index k, const Index j) const
|
||||
{
|
||||
checkDataIndex(k);
|
||||
checkYDim(j);
|
||||
|
||||
return !(yDataIndex_[j].find(k) == yDataIndex_[j].end());
|
||||
}
|
||||
|
||||
bool FitInterface::isXExact(const Index i) const
|
||||
{
|
||||
checkXDim(i);
|
||||
|
||||
return xIsExact_[i];
|
||||
}
|
||||
|
||||
bool FitInterface::isXUsed(const Index r, const Index i, const bool inFit) const
|
||||
{
|
||||
vector<Index> v;
|
||||
|
||||
checkXDim(i);
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
for (auto &p: yDataIndex_[j])
|
||||
{
|
||||
if (p.second or !inFit)
|
||||
{
|
||||
v = dataCoord(p.first);
|
||||
if (v[i] == r)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
bool FitInterface::isFitPoint(const Index k, const Index j) const
|
||||
{
|
||||
checkPoint(k, j);
|
||||
|
||||
return yDataIndex_[j].at(k);
|
||||
}
|
||||
|
||||
bool FitInterface::isXXCorrelated(const Index r1, const Index i1,
|
||||
const Index r2, const Index i2) const
|
||||
{
|
||||
array<Index, 4> c{{r1, i1, r2, i2}};
|
||||
auto it = xxCorr_.find(c);
|
||||
|
||||
return (it != xxCorr_.end());
|
||||
}
|
||||
|
||||
bool FitInterface::isYYCorrelated(const Index k1, const Index j1,
|
||||
const Index k2, const Index j2) const
|
||||
{
|
||||
array<Index, 4> c{{k1, j1, k2, j2}};
|
||||
auto it = yyCorr_.find(c);
|
||||
|
||||
return (it != yyCorr_.end());
|
||||
}
|
||||
|
||||
bool FitInterface::isXYCorrelated(const Index r, const Index i,
|
||||
const Index k, const Index j) const
|
||||
{
|
||||
array<Index, 4> c{{r, i, k, j}};
|
||||
auto it = xyCorr_.find(c);
|
||||
|
||||
return (it != xyCorr_.end());
|
||||
}
|
||||
|
||||
bool FitInterface::hasCorrelations(void) const
|
||||
{
|
||||
return ((xxCorr_.size() != 0) or (yyCorr_.size() != 0)
|
||||
or (xyCorr_.size() != 0));
|
||||
}
|
||||
|
||||
// make correlation filter for fit variance matrix /////////////////////////////
|
||||
DMat FitInterface::makeCorrFilter(void)
|
||||
{
|
||||
updateLayout();
|
||||
|
||||
DMat f = DMat::Identity(layout.totalSize, layout.totalSize);
|
||||
Index row, col;
|
||||
|
||||
for (auto &c: xxCorr_)
|
||||
{
|
||||
row = indX(c[0], c[1]);
|
||||
col = indX(c[2], c[3]);
|
||||
if ((row != -1) and (col != -1))
|
||||
{
|
||||
f(row, col) = 1.;
|
||||
f(col, row) = 1.;
|
||||
}
|
||||
}
|
||||
for (auto &c: yyCorr_)
|
||||
{
|
||||
row = indY(c[0], c[1]);
|
||||
col = indY(c[2], c[3]);
|
||||
if ((row != -1) and (col != -1))
|
||||
{
|
||||
f(row, col) = 1.;
|
||||
f(col, row) = 1.;
|
||||
}
|
||||
}
|
||||
for (auto &c: xyCorr_)
|
||||
{
|
||||
row = indX(c[0], c[1]);
|
||||
col = indY(c[2], c[3]);
|
||||
if ((row != -1) and (col != -1))
|
||||
{
|
||||
f(row, col) = 1.;
|
||||
f(col, row) = 1.;
|
||||
}
|
||||
}
|
||||
|
||||
return f;
|
||||
}
|
||||
|
||||
// schedule variance matrix initialization /////////////////////////////////////
|
||||
void FitInterface::scheduleFitVarMatInit(const bool init)
|
||||
{
|
||||
initVarMat_ = init;
|
||||
}
|
||||
|
||||
// register a data point ///////////////////////////////////////////////////////
|
||||
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;
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
void FitInterface::updateDataCoord(void) const
|
||||
{
|
||||
FitInterface * modThis = const_cast<FitInterface *>(this);
|
||||
|
||||
if (initDataCoord_)
|
||||
{
|
||||
modThis->dataCoord_.clear();
|
||||
for (auto k: getDataIndexSet())
|
||||
{
|
||||
modThis->dataCoord_[k] = rowMajToCoord(k);
|
||||
}
|
||||
modThis->initDataCoord_ = false;
|
||||
}
|
||||
}
|
||||
|
||||
// global layout management ////////////////////////////////////////////////////
|
||||
void FitInterface::scheduleLayoutInit(void)
|
||||
{
|
||||
initLayout_ = true;
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
bool FitInterface::initVarMat(void) const
|
||||
{
|
||||
return initVarMat_;
|
||||
}
|
||||
|
||||
void FitInterface::updateLayout(void) const
|
||||
{
|
||||
if (initLayout_)
|
||||
{
|
||||
FitInterface * modThis = const_cast<FitInterface *>(this);
|
||||
Layout & l = modThis->layout;
|
||||
Index size, ifit;
|
||||
vector<Index> v;
|
||||
|
||||
l.nXFitDim = 0;
|
||||
l.nYFitDim = 0;
|
||||
l.totalSize = 0;
|
||||
l.totalXSize = 0;
|
||||
l.totalYSize = 0;
|
||||
l.xSize.clear();
|
||||
l.ySize.clear();
|
||||
l.dataIndexSet.clear();
|
||||
l.xDim.clear();
|
||||
l.yDim.clear();
|
||||
l.xFitDim.clear();
|
||||
l.yFitDim.clear();
|
||||
l.x.clear();
|
||||
l.y.clear();
|
||||
l.xFit.clear();
|
||||
l.yFit.clear();
|
||||
ifit = 0;
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
{
|
||||
if (!xIsExact_[i])
|
||||
{
|
||||
l.nXFitDim++;
|
||||
size = getXFitSize(i);
|
||||
l.xSize.push_back(size);
|
||||
l.totalXSize += size;
|
||||
l.xDim.push_back(i);
|
||||
l.xFitDim.push_back(layout.xDim.size() - 1);
|
||||
l.x.push_back(vector<Index>());
|
||||
l.xFit.push_back(vector<Index>());
|
||||
for (Index r = 0; r < getXSize(i); ++r)
|
||||
{
|
||||
if (isXUsed(r, i))
|
||||
{
|
||||
l.x[ifit].push_back(r);
|
||||
l.xFit[i].push_back(layout.x[ifit].size() - 1);
|
||||
}
|
||||
else
|
||||
{
|
||||
l.xFit[i].push_back(-1);
|
||||
}
|
||||
}
|
||||
ifit++;
|
||||
}
|
||||
else
|
||||
{
|
||||
l.xFitDim.push_back(-1);
|
||||
l.xFit.push_back(vector<Index>());
|
||||
for (Index r = 0; r < getXSize(i); ++r)
|
||||
{
|
||||
l.xFit[i].push_back(-1);
|
||||
}
|
||||
}
|
||||
}
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
Index s = 0;
|
||||
|
||||
l.nYFitDim++;
|
||||
size = getYFitSize(j);
|
||||
l.ySize.push_back(size);
|
||||
l.totalYSize += size;
|
||||
l.yDim.push_back(j);
|
||||
l.yFitDim.push_back(layout.yDim.size() - 1);
|
||||
l.y.push_back(vector<Index>());
|
||||
l.yFit.push_back(vector<Index>());
|
||||
l.data.push_back(vector<Index>());
|
||||
l.yFitFromData.push_back(map<Index, Index>());
|
||||
for (auto &p: yDataIndex_[j])
|
||||
{
|
||||
if (p.second)
|
||||
{
|
||||
l.dataIndexSet.insert(p.first);
|
||||
l.y[j].push_back(s);
|
||||
l.yFit[j].push_back(layout.y[j].size() - 1);
|
||||
l.data[j].push_back(p.first);
|
||||
l.yFitFromData[j][p.first] = layout.y[j].size() - 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
l.yFit[j].push_back(-1);
|
||||
l.yFitFromData[j][p.first] = -1;
|
||||
}
|
||||
s++;
|
||||
}
|
||||
}
|
||||
l.totalSize = layout.totalXSize + layout.totalYSize;
|
||||
l.nXFitDim = static_cast<Index>(layout.xSize.size());
|
||||
l.nYFitDim = static_cast<Index>(layout.ySize.size());
|
||||
l.xIndFromData.resize(getMaxDataIndex());
|
||||
for (Index k: layout.dataIndexSet)
|
||||
{
|
||||
v = dataCoord(k);
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
{
|
||||
l.xIndFromData[k].push_back(indX(v[i], i));
|
||||
}
|
||||
}
|
||||
modThis->initLayout_ = false;
|
||||
}
|
||||
}
|
||||
|
||||
Index FitInterface::indX(const Index r, const Index i) const
|
||||
{
|
||||
Index ind = -1;
|
||||
|
||||
if (layout.xFit[i][r] != -1)
|
||||
{
|
||||
Index ifit = layout.xFitDim[i], rfit = layout.xFit[i][r];
|
||||
|
||||
ind = layout.totalYSize;
|
||||
for (Index a = 0; a < ifit; ++a)
|
||||
{
|
||||
ind += layout.xSize[a];
|
||||
}
|
||||
ind += rfit;
|
||||
}
|
||||
|
||||
return ind;
|
||||
}
|
||||
|
||||
Index FitInterface::indY(const Index k, const Index j) const
|
||||
{
|
||||
Index ind = -1;
|
||||
|
||||
if (layout.yFitFromData[j].at(k) != -1)
|
||||
{
|
||||
Index jfit = layout.yFitDim[j], sfit = layout.yFitFromData[j].at(k);
|
||||
|
||||
ind = 0;
|
||||
for (Index b = 0; b < jfit; ++b)
|
||||
{
|
||||
ind += layout.ySize[b];
|
||||
}
|
||||
ind += sfit;
|
||||
}
|
||||
|
||||
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)
|
||||
{
|
||||
out << "X dimensions: " << f.getNXDim() << endl;
|
||||
for (Index i = 0; i < f.getNXDim(); ++i)
|
||||
{
|
||||
out << " * " << i << " \"" << f.xName().getName(i) << "\": ";
|
||||
out << f.getXSize(i) << " value(s)";
|
||||
if (f.isXExact(i))
|
||||
{
|
||||
out << " (assumed exact)";
|
||||
}
|
||||
out << endl;
|
||||
}
|
||||
out << "Y dimensions: " << f.getNYDim() << endl;
|
||||
for (Index j = 0; j < f.getNYDim(); ++j)
|
||||
{
|
||||
out << " * " << j << " \"" << f.yName().getName(j) << "\": ";
|
||||
out << f.getYSize(j) << " value(s)" << endl;
|
||||
for (auto &p: f.yDataIndex_[j])
|
||||
{
|
||||
out << " " << setw(3) << p.first << " (";
|
||||
for (auto vi: f.dataCoord(p.first))
|
||||
{
|
||||
out << vi << ",";
|
||||
}
|
||||
out << "\b) fit: " << (p.second ? "true" : "false") << endl;
|
||||
}
|
||||
}
|
||||
out << "X/X correlations (r1 i1 r2 i2): ";
|
||||
if (f.xxCorr_.empty())
|
||||
{
|
||||
out << "no" << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
out << endl;
|
||||
for (auto &c: f.xxCorr_)
|
||||
{
|
||||
out << " * ";
|
||||
for (auto i: c)
|
||||
{
|
||||
out << i << " ";
|
||||
}
|
||||
out << endl;
|
||||
}
|
||||
}
|
||||
out << "Y/Y correlations (k1 j1 k2 j2): ";
|
||||
if (f.yyCorr_.empty())
|
||||
{
|
||||
out << "no" << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
out << endl;
|
||||
for (auto &c: f.yyCorr_)
|
||||
{
|
||||
out << " * ";
|
||||
for (auto i: c)
|
||||
{
|
||||
out << i << " ";
|
||||
}
|
||||
out << endl;
|
||||
}
|
||||
}
|
||||
out << "X/Y correlations (r i k j): ";
|
||||
if (f.xyCorr_.empty())
|
||||
{
|
||||
out << "no";
|
||||
}
|
||||
else
|
||||
{
|
||||
out << endl;
|
||||
for (auto &c: f.xyCorr_)
|
||||
{
|
||||
out << " * ";
|
||||
for (auto i: c)
|
||||
{
|
||||
out << i << " ";
|
||||
}
|
||||
out << endl;
|
||||
}
|
||||
}
|
||||
|
||||
return out;
|
||||
}
|
||||
228
lib/Statistics/FitInterface.hpp
Normal file
228
lib/Statistics/FitInterface.hpp
Normal file
@ -0,0 +1,228 @@
|
||||
/*
|
||||
* FitInterface.hpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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/Core/Mat.hpp>
|
||||
|
||||
BEGIN_LATAN_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* FitInterface *
|
||||
******************************************************************************/
|
||||
class FitInterface
|
||||
{
|
||||
private:
|
||||
typedef struct
|
||||
{
|
||||
Index nXFitDim, nYFitDim;
|
||||
// X/Y block sizes
|
||||
Index totalSize, totalXSize, totalYSize;
|
||||
// size of each X/Y dimension
|
||||
std::vector<Index> xSize, ySize;
|
||||
// set of active data indices
|
||||
std::set<Index> dataIndexSet;
|
||||
// lookup tables
|
||||
// xDim : x fit dim ifit -> x dim i
|
||||
// x : x fit point ifit,rfit -> x point r
|
||||
// xFitDim : x dim i -> x fit dim ifit (-1 if empty)
|
||||
// xFit : x point i,r -> x fit point rfit (-1 if empty)
|
||||
// data : y fit point jfit,sfit -> y point index k
|
||||
// yFitFromData: y point index k,j -> y fit point sfit (-1 if empty)
|
||||
// xIndFromData: data index k -> index of coordinates of associated x
|
||||
std::vector<Index> xDim, yDim, xFitDim, yFitDim;
|
||||
std::vector<std::vector<Index>> x, y, data, xFit, yFit;
|
||||
std::vector<std::map<Index, Index>> yFitFromData;
|
||||
// no map here for fit performance
|
||||
std::vector<std::vector<Index>> xIndFromData;
|
||||
} Layout;
|
||||
public:
|
||||
// constructor
|
||||
FitInterface(void);
|
||||
// destructor
|
||||
virtual ~FitInterface(void) = default;
|
||||
// copy object (not as a constructor to be accessed from derived class)
|
||||
void copyInterface(const FitInterface &d);
|
||||
// add dimensions
|
||||
void addXDim(const Index nData, const std::string name = "",
|
||||
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;
|
||||
const std::set<Index> & getDataIndexSet(void) const;
|
||||
double getSvdTolerance(void) const;
|
||||
void setSvdTolerance(const double &tol);
|
||||
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;
|
||||
const std::vector<Index> & dataCoord(const Index k) const;
|
||||
// enable fit points
|
||||
void fitPoint(const bool isFitPoint, const Index k, const Index j = 0);
|
||||
// variance interface
|
||||
void assumeXExact(const bool isExact, const Index i);
|
||||
void assumeXXCorrelated(const bool isCorr, const Index r1, const Index i1,
|
||||
const Index r2, const Index i2);
|
||||
void assumeXXCorrelated(const bool isCorr, const Index i1, const Index i2);
|
||||
void assumeYYCorrelated(const bool isCorr, const Index k1, const Index j1,
|
||||
const Index k2, const Index j2);
|
||||
void assumeYYCorrelated(const bool isCorr, const Index j1, const Index j2);
|
||||
void assumeXYCorrelated(const bool isCorr, const Index r, const Index i,
|
||||
const Index k, const Index j);
|
||||
void assumeXYCorrelated(const bool isCorr, const Index i, const Index j);
|
||||
// tests
|
||||
bool pointExists(const Index k) const;
|
||||
bool pointExists(const Index k, const Index j) const;
|
||||
bool isXExact(const Index i) const;
|
||||
bool isXUsed(const Index r, const Index i, const bool inFit = true) const;
|
||||
bool isFitPoint(const Index k, const Index j) const;
|
||||
bool isXXCorrelated(const Index r1, const Index i1, const Index r2,
|
||||
const Index i2) const;
|
||||
bool isYYCorrelated(const Index k1, const Index j1, const Index k2,
|
||||
const Index j2) const;
|
||||
bool isXYCorrelated(const Index r, const Index i, const Index k,
|
||||
const Index j) const;
|
||||
bool hasCorrelations(void) const;
|
||||
// make correlation filter for fit variance matrix
|
||||
DMat makeCorrFilter(void);
|
||||
// schedule variance matrix initialization
|
||||
void scheduleFitVarMatInit(const bool init = true);
|
||||
// IO
|
||||
friend std::ostream & operator<<(std::ostream &out, FitInterface &f);
|
||||
protected:
|
||||
// register a data point
|
||||
void registerDataPoint(const Index k, const Index j = 0);
|
||||
// add correlation to a set
|
||||
static void addCorr(std::set<std::array<Index, 4>> &s, const bool isCorr,
|
||||
const std::array<Index, 4> &c);
|
||||
// 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) const;
|
||||
// global layout management
|
||||
void scheduleLayoutInit(void);
|
||||
bool initVarMat(void) const;
|
||||
void updateLayout(void) const;
|
||||
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::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};
|
||||
bool initVarMat_{true};
|
||||
bool initDataCoord_{true};
|
||||
double svdTol_{1.e-10};
|
||||
};
|
||||
|
||||
std::ostream & operator<<(std::ostream &out, FitInterface &f);
|
||||
|
||||
/******************************************************************************
|
||||
* FitInterface template implementation *
|
||||
******************************************************************************/
|
||||
// Y dimension index helper ////////////////////////////////////////////////////
|
||||
template <typename... Ts>
|
||||
Index FitInterface::dataIndex(const Ts... coords) const
|
||||
{
|
||||
static_assert(static_or<std::is_convertible<Index, Ts>::value...>::value,
|
||||
"fitPoint arguments are not compatible with Index");
|
||||
|
||||
const std::vector<Index> coord = {coords...};
|
||||
|
||||
return dataIndex(coord);
|
||||
}
|
||||
|
||||
/******************************************************************************
|
||||
* error check macros *
|
||||
******************************************************************************/
|
||||
#define checkXDim(i)\
|
||||
if ((i) >= getNXDim())\
|
||||
{\
|
||||
LATAN_ERROR(Range, "X dimension " + strFrom(i) + " out of range");\
|
||||
}
|
||||
|
||||
#define checkXIndex(vi, i)\
|
||||
if ((vi) >= getXSize(i))\
|
||||
{\
|
||||
LATAN_ERROR(Range, "index " + strFrom(vi) + " in X dimension "\
|
||||
+ strFrom(i) + " out of range");\
|
||||
}
|
||||
|
||||
#define checkYDim(j)\
|
||||
if ((j) >= getNYDim())\
|
||||
{\
|
||||
LATAN_ERROR(Range, "Y dimension " + strFrom(j) + " out of range");\
|
||||
}
|
||||
|
||||
#define checkDataIndex(k)\
|
||||
if ((k) >= getMaxDataIndex())\
|
||||
{\
|
||||
LATAN_ERROR(Range, "data point index " + strFrom(k) + " invalid");\
|
||||
}
|
||||
|
||||
#define checkDataCoord(v)\
|
||||
if (static_cast<Index>((v).size()) != getNXDim())\
|
||||
{\
|
||||
LATAN_ERROR(Size, "number of coordinates and number of X dimensions "\
|
||||
"mismatch");\
|
||||
}\
|
||||
for (unsigned int i_ = 0; i_ < (v).size(); ++i_)\
|
||||
{\
|
||||
checkXIndex((v)[i_], i_);\
|
||||
}
|
||||
|
||||
#define checkPoint(k, j)\
|
||||
if (!pointExists(k, j))\
|
||||
{\
|
||||
LATAN_ERROR(Range, "no data point in Y dimension " + strFrom(j)\
|
||||
+ " with index " + strFrom(k));\
|
||||
}
|
||||
|
||||
END_LATAN_NAMESPACE
|
||||
|
||||
#endif // Latan_FitInterface_hpp_
|
||||
228
lib/Statistics/Histogram.cpp
Normal file
228
lib/Statistics/Histogram.cpp
Normal file
@ -0,0 +1,228 @@
|
||||
/*
|
||||
* Histogram.cpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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/Statistics/Histogram.hpp>
|
||||
#include <LatAnalyze/includes.hpp>
|
||||
#include <gsl/gsl_histogram.h>
|
||||
#include <gsl/gsl_sf.h>
|
||||
#include <gsl/gsl_sort.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Latan;
|
||||
|
||||
#define DECL_GSL_HIST(h) \
|
||||
gsl_histogram h{static_cast<size_t>(bin_.size()), x_.data(), bin_.data()}
|
||||
#define DECL_CONST_GSL_HIST(h) \
|
||||
const gsl_histogram h{static_cast<size_t>(bin_.size()),\
|
||||
const_cast<double *>(x_.data()),\
|
||||
const_cast<double *>(bin_.data())}
|
||||
|
||||
/******************************************************************************
|
||||
* Histogram implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
Histogram::Histogram(const DVec &data, const double xMin, const double xMax,
|
||||
const Index nBin)
|
||||
: Histogram()
|
||||
{
|
||||
setFromData(data, xMin, xMax, nBin);
|
||||
}
|
||||
|
||||
Histogram::Histogram(const DVec &data, const DVec &w, const double xMin,
|
||||
const double xMax, const Index nBin)
|
||||
: Histogram()
|
||||
{
|
||||
setFromData(data, w, xMin, xMax, nBin);
|
||||
}
|
||||
|
||||
// resize //////////////////////////////////////////////////////////////////////
|
||||
void Histogram::resize(const Index nBin)
|
||||
{
|
||||
x_.resize(nBin + 1);
|
||||
bin_.resize(nBin);
|
||||
}
|
||||
|
||||
// generate from data //////////////////////////////////////////////////////////
|
||||
void Histogram::setFromData(const DVec &data, const DVec &w, const double xMin,
|
||||
const double xMax, const Index nBin)
|
||||
{
|
||||
if (data.size() != w.size())
|
||||
{
|
||||
LATAN_ERROR(Size, "data vector and weight vector size mismatch");
|
||||
}
|
||||
resize(nBin);
|
||||
data_ = data.array();
|
||||
w_ = w.array();
|
||||
xMax_ = xMax;
|
||||
xMin_ = xMin;
|
||||
makeHistogram();
|
||||
}
|
||||
|
||||
void Histogram::setFromData(const DVec &data, const double xMin,
|
||||
const double xMax, const Index nBin)
|
||||
{
|
||||
resize(nBin);
|
||||
data_ = data.array();
|
||||
xMax_ = xMax;
|
||||
xMin_ = xMin;
|
||||
w_.resize(data.size());
|
||||
w_.fill(1.);
|
||||
makeHistogram();
|
||||
}
|
||||
|
||||
// histogram calculation ///////////////////////////////////////////////////////
|
||||
void Histogram::makeHistogram(void)
|
||||
{
|
||||
DECL_GSL_HIST(h);
|
||||
|
||||
gsl_histogram_set_ranges_uniform(&h, xMin_, xMax_);
|
||||
FOR_STAT_ARRAY(data_, i)
|
||||
{
|
||||
gsl_histogram_accumulate(&h, data_[i], w_[i]);
|
||||
}
|
||||
total_ = w_.sum();
|
||||
sortIndices();
|
||||
computeNorm();
|
||||
}
|
||||
|
||||
// generate sorted indices /////////////////////////////////////////////////////
|
||||
void Histogram::sortIndices(void)
|
||||
{
|
||||
sInd_.resize(data_.size());
|
||||
gsl_sort_index(sInd_.data(), data_.data(), 1, data_.size());
|
||||
}
|
||||
|
||||
// compute normalization factor ////////////////////////////////////////////////
|
||||
void Histogram::computeNorm(void)
|
||||
{
|
||||
norm_ = static_cast<double>(bin_.size())/(total_*(xMax_ - xMin_));
|
||||
}
|
||||
|
||||
// normalize as a probablility /////////////////////////////////////////////////
|
||||
void Histogram::normalize(const bool n)
|
||||
{
|
||||
normalize_ = n;
|
||||
}
|
||||
|
||||
bool Histogram::isNormalized(void) const
|
||||
{
|
||||
return normalize_;
|
||||
}
|
||||
|
||||
// access //////////////////////////////////////////////////////////////////////
|
||||
Index Histogram::size(void) const
|
||||
{
|
||||
return bin_.size();
|
||||
}
|
||||
|
||||
const StatArray<double> & Histogram::getData(void) const
|
||||
{
|
||||
return data_;
|
||||
}
|
||||
|
||||
const StatArray<double> & Histogram::getWeight(void) const
|
||||
{
|
||||
return w_;
|
||||
}
|
||||
|
||||
double Histogram::getX(const Index i) const
|
||||
{
|
||||
return x_(i);
|
||||
}
|
||||
|
||||
double Histogram::operator[](const Index i) const
|
||||
{
|
||||
return bin_(i)*(isNormalized() ? norm_ : 1.);
|
||||
}
|
||||
|
||||
double Histogram::operator()(const double x) const
|
||||
{
|
||||
size_t i;
|
||||
|
||||
DECL_CONST_GSL_HIST(h);
|
||||
|
||||
gsl_histogram_find(&h, x, &i);
|
||||
|
||||
return (*this)[static_cast<Index>(i)];
|
||||
}
|
||||
|
||||
// percentiles & confidence interval ///////////////////////////////////////////
|
||||
double Histogram::percentile(const double p) const
|
||||
{
|
||||
if ((p < 0.0) or (p > 100.0))
|
||||
{
|
||||
LATAN_ERROR(Range, "percentile (" + strFrom(p) + ")"
|
||||
" is outside the [0, 100] range");
|
||||
}
|
||||
|
||||
// cf. http://en.wikipedia.org/wiki/Percentile
|
||||
double wPSum, p_i, p_im1, w_i, res = 0.;
|
||||
bool haveResult;
|
||||
|
||||
wPSum = w_[sInd_[0]];
|
||||
p_i = (100./total_)*wPSum*0.5;
|
||||
if (p < p_i)
|
||||
{
|
||||
res = data_[sInd_[0]];
|
||||
}
|
||||
else
|
||||
{
|
||||
haveResult = false;
|
||||
p_im1 = p_i;
|
||||
for (Index i = 1; i < data_.size(); ++i)
|
||||
{
|
||||
w_i = w_[sInd_[i]];
|
||||
wPSum += w_i;
|
||||
p_i = (100./total_)*(wPSum-0.5*w_i);
|
||||
if ((p >= p_im1) and (p < p_i))
|
||||
{
|
||||
double d_i = data_[sInd_[i]], d_im1 = data_[sInd_[i-1]];
|
||||
|
||||
res = d_im1 + (p-p_im1)/(p_i-p_im1)*(d_i-d_im1);
|
||||
haveResult = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (!haveResult)
|
||||
{
|
||||
res = data_[sInd_[data_.size()-1]];
|
||||
}
|
||||
}
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
double Histogram::median(void) const
|
||||
{
|
||||
return percentile(50.);
|
||||
}
|
||||
|
||||
pair<double, double> Histogram::confidenceInterval(const double nSigma) const
|
||||
{
|
||||
pair<double, double> interval, p;
|
||||
double cl;
|
||||
|
||||
cl = gsl_sf_erf(nSigma/sqrt(2.));
|
||||
p.first = 50.*(1. - cl);
|
||||
p.second = 50.*(1. + cl);
|
||||
interval.first = percentile(p.first);
|
||||
interval.second = percentile(p.second);
|
||||
|
||||
return interval;
|
||||
}
|
||||
80
lib/Statistics/Histogram.hpp
Normal file
80
lib/Statistics/Histogram.hpp
Normal file
@ -0,0 +1,80 @@
|
||||
/*
|
||||
* Histogram.hpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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_Histogram_hpp_
|
||||
#define Latan_Histogram_hpp_
|
||||
|
||||
#include <LatAnalyze/Global.hpp>
|
||||
#include <LatAnalyze/Statistics/StatArray.hpp>
|
||||
|
||||
BEGIN_LATAN_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* Histogram class *
|
||||
******************************************************************************/
|
||||
class Histogram
|
||||
{
|
||||
public:
|
||||
// constructor
|
||||
Histogram(void) = default;
|
||||
Histogram(const DVec &data, const double xMin, const double xMax,
|
||||
const Index nBin);
|
||||
Histogram(const DVec &data, const DVec &w, const double xMin,
|
||||
const double xMax, const Index nBin);
|
||||
// destructor
|
||||
virtual ~Histogram(void) = default;
|
||||
// generate from data
|
||||
void setFromData(const DVec &data, const double xMin, const double xMax,
|
||||
const Index nBin);
|
||||
void setFromData(const DVec &data, const DVec &w, const double xMin,
|
||||
const double xMax, const Index nBin);
|
||||
// normalize as a probablility
|
||||
void normalize(const bool n = true);
|
||||
bool isNormalized(void) const;
|
||||
// access
|
||||
Index size(void) const;
|
||||
const StatArray<double> & getData(void) const;
|
||||
const StatArray<double> & getWeight(void) const;
|
||||
double getX(const Index i) const;
|
||||
double operator[](const Index i) const;
|
||||
double operator()(const double x) const;
|
||||
// percentiles & confidence interval
|
||||
double percentile(const double p) const;
|
||||
double median(void) const;
|
||||
std::pair<double, double> confidenceInterval(const double nSigma) const;
|
||||
private:
|
||||
// resize
|
||||
void resize(const Index nBin);
|
||||
// histogram calculation
|
||||
void makeHistogram(void);
|
||||
// generate sorted indices
|
||||
void sortIndices(void);
|
||||
// compute normalization factor
|
||||
void computeNorm(void);
|
||||
private:
|
||||
StatArray<double> data_, w_;
|
||||
DVec x_, bin_;
|
||||
Vec<size_t> sInd_;
|
||||
double total_, norm_, xMax_, xMin_;
|
||||
bool normalize_{false};
|
||||
};
|
||||
|
||||
END_LATAN_NAMESPACE
|
||||
|
||||
#endif // Latan_Histogram_hpp_
|
||||
385
lib/Statistics/MatSample.hpp
Normal file
385
lib/Statistics/MatSample.hpp
Normal file
@ -0,0 +1,385 @@
|
||||
/*
|
||||
* MatSample.hpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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_MatSample_hpp_
|
||||
#define Latan_MatSample_hpp_
|
||||
|
||||
#include <LatAnalyze/Global.hpp>
|
||||
#include <LatAnalyze/Core/Mat.hpp>
|
||||
#include <LatAnalyze/Statistics/StatArray.hpp>
|
||||
|
||||
BEGIN_LATAN_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* matrix sample class *
|
||||
******************************************************************************/
|
||||
#define SCAL_OP_RETURN(op, s, x) s.unaryExpr(\
|
||||
std::bind(MatSample<T>::scalar##op,\
|
||||
std::placeholders::_1, x))
|
||||
|
||||
template <typename T>
|
||||
class MatSample: public Sample<Mat<T>>
|
||||
{
|
||||
public:
|
||||
// block type template
|
||||
template <class S>
|
||||
class BlockTemplate
|
||||
{
|
||||
private:
|
||||
typedef typename std::remove_const<S>::type NonConstType;
|
||||
public:
|
||||
// constructors
|
||||
BlockTemplate(S &sample, const Index i, const Index j, const Index nRow,
|
||||
const Index nCol);
|
||||
BlockTemplate(BlockTemplate<NonConstType> &b);
|
||||
BlockTemplate(BlockTemplate<NonConstType> &&b);
|
||||
// destructor
|
||||
~BlockTemplate(void) = default;
|
||||
// access
|
||||
S & getSample(void);
|
||||
const S & getSample(void) const;
|
||||
Index getStartRow(void) const;
|
||||
Index getStartCol(void) const;
|
||||
Index getNRow(void) const;
|
||||
Index getNCol(void) const;
|
||||
// assignement operators
|
||||
BlockTemplate<S> & operator=(const S &sample);
|
||||
BlockTemplate<S> & operator=(const S &&sample);
|
||||
private:
|
||||
S &sample_;
|
||||
const Index i_, j_, nRow_, nCol_;
|
||||
};
|
||||
// block types
|
||||
typedef BlockTemplate<Sample<Mat<T>>> Block;
|
||||
typedef const BlockTemplate<const Sample<Mat<T>>> ConstBlock;
|
||||
public:
|
||||
// constructors
|
||||
MatSample(void) = default;
|
||||
MatSample(const Index nSample);
|
||||
MatSample(const Index nSample, const Index nRow, const Index nCol);
|
||||
MatSample(ConstBlock &sampleBlock);
|
||||
MatSample(ConstBlock &&sampleBlock);
|
||||
EIGEN_EXPR_CTOR(MatSample, MatSample<T>, Sample<Mat<T>>, ArrayExpr)
|
||||
// destructor
|
||||
virtual ~MatSample(void) = default;
|
||||
// assignement operator
|
||||
MatSample<T> & operator=(Block &sampleBlock);
|
||||
MatSample<T> & operator=(Block &&sampleBlock);
|
||||
MatSample<T> & operator=(ConstBlock &sampleBlock);
|
||||
MatSample<T> & operator=(ConstBlock &&sampleBlock);
|
||||
// product/division by scalar operators (not provided by Eigen)
|
||||
static inline Mat<T> scalarMul(const Mat<T> &m, const T &x)
|
||||
{
|
||||
return m*x;
|
||||
}
|
||||
static inline Mat<T> scalarDiv(const Mat<T> &m, const T &x)
|
||||
{
|
||||
return m/x;
|
||||
}
|
||||
MatSample<T> & operator*=(const T &x);
|
||||
MatSample<T> & operator*=(const T &&x);
|
||||
MatSample<T> & operator/=(const T &x);
|
||||
MatSample<T> & operator/=(const T &&x);
|
||||
// block access
|
||||
ConstBlock block(const Index i, const Index j, const Index nRow,
|
||||
const Index nCol) const;
|
||||
Block block(const Index i, const Index j, const Index nRow,
|
||||
const Index nCol);
|
||||
// resize all matrices
|
||||
void resizeMat(const Index nRow, const Index nCol);
|
||||
};
|
||||
|
||||
// non-member operators
|
||||
template <typename T>
|
||||
inline auto operator*(MatSample<T> s, const T &x)
|
||||
->decltype(SCAL_OP_RETURN(Mul, s, x))
|
||||
{
|
||||
return SCAL_OP_RETURN(Mul, s, x);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline auto operator*(MatSample<T> s, const T &&x)
|
||||
->decltype(SCAL_OP_RETURN(Mul, s, x))
|
||||
{
|
||||
return SCAL_OP_RETURN(Mul, s, x);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline auto operator*(const T &x, MatSample<T> s)->decltype(s*x)
|
||||
{
|
||||
return s*x;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline auto operator*(const T &&x, MatSample<T> s)->decltype(s*x)
|
||||
{
|
||||
return s*x;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline auto operator/(MatSample<T> s, const T &x)
|
||||
->decltype(SCAL_OP_RETURN(Div, s, x))
|
||||
{
|
||||
return SCAL_OP_RETURN(Div, s, x);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline auto operator/(MatSample<T> s, const T &&x)
|
||||
->decltype(SCAL_OP_RETURN(Div, s, x))
|
||||
{
|
||||
return SCAL_OP_RETURN(Div, s, x);
|
||||
}
|
||||
|
||||
// type aliases
|
||||
typedef MatSample<double> DMatSample;
|
||||
typedef MatSample<std::complex<double>> CMatSample;
|
||||
|
||||
/******************************************************************************
|
||||
* Block template implementation *
|
||||
******************************************************************************/
|
||||
// constructors ////////////////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
template <class S>
|
||||
MatSample<T>::BlockTemplate<S>::BlockTemplate(S &sample, const Index i,
|
||||
const Index j, const Index nRow,
|
||||
const Index nCol)
|
||||
: sample_(sample)
|
||||
, i_(i)
|
||||
, j_(j)
|
||||
, nRow_(nRow)
|
||||
, nCol_(nCol)
|
||||
{}
|
||||
|
||||
template <typename T>
|
||||
template <class S>
|
||||
MatSample<T>::BlockTemplate<S>::BlockTemplate(BlockTemplate<NonConstType> &b)
|
||||
: sample_(b.getSample())
|
||||
, i_(b.getStartRow())
|
||||
, j_(b.getStartCol())
|
||||
, nRow_(b.getNRow())
|
||||
, nCol_(b.getNCol())
|
||||
{}
|
||||
|
||||
template <typename T>
|
||||
template <class S>
|
||||
MatSample<T>::BlockTemplate<S>::BlockTemplate(BlockTemplate<NonConstType> &&b)
|
||||
: BlockTemplate(b)
|
||||
{}
|
||||
|
||||
// access //////////////////////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
template <class S>
|
||||
S & MatSample<T>::BlockTemplate<S>::getSample(void)
|
||||
{
|
||||
return sample_;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
template <class S>
|
||||
const S & MatSample<T>::BlockTemplate<S>::getSample(void) const
|
||||
{
|
||||
return sample_;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
template <class S>
|
||||
Index MatSample<T>::BlockTemplate<S>::getStartRow(void) const
|
||||
{
|
||||
return i_;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
template <class S>
|
||||
Index MatSample<T>::BlockTemplate<S>::getStartCol(void) const
|
||||
{
|
||||
return j_;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
template <class S>
|
||||
Index MatSample<T>::BlockTemplate<S>::getNRow(void) const
|
||||
{
|
||||
return nRow_;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
template <class S>
|
||||
Index MatSample<T>::BlockTemplate<S>::getNCol(void) const
|
||||
{
|
||||
return nCol_;
|
||||
}
|
||||
|
||||
// assignement operators ///////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
template <class S>
|
||||
typename MatSample<T>::template BlockTemplate<S> &
|
||||
MatSample<T>::BlockTemplate<S>::operator=(const S &sample)
|
||||
{
|
||||
FOR_STAT_ARRAY(sample_, s)
|
||||
{
|
||||
sample_[s].block(i_, j_, nRow_, nCol_) = sample[s];
|
||||
}
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
template <class S>
|
||||
typename MatSample<T>::template BlockTemplate<S> &
|
||||
MatSample<T>::BlockTemplate<S>::operator=(const S &&sample)
|
||||
{
|
||||
*this = sample;
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
/******************************************************************************
|
||||
* DMatSample implementation *
|
||||
******************************************************************************/
|
||||
// constructors ////////////////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
MatSample<T>::MatSample(const Index nSample)
|
||||
: Sample<Mat<T>>(nSample)
|
||||
{}
|
||||
|
||||
template <typename T>
|
||||
MatSample<T>::MatSample(const Index nSample, const Index nRow,
|
||||
const Index nCol)
|
||||
: MatSample(nSample)
|
||||
{
|
||||
resizeMat(nRow, nCol);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
MatSample<T>::MatSample(ConstBlock &sampleBlock)
|
||||
: MatSample(sampleBlock.getSample().size(), sampleBlock.getNRow(),
|
||||
sampleBlock.getNCol())
|
||||
{
|
||||
const MatSample<T> &sample = sampleBlock.getSample();
|
||||
|
||||
this->resize(sample.size());
|
||||
FOR_STAT_ARRAY(*this, s)
|
||||
{
|
||||
(*this)[s] = sample[s].block(sampleBlock.getStartRow(),
|
||||
sampleBlock.getStartCol(),
|
||||
sampleBlock.getNRow(),
|
||||
sampleBlock.getNCol());
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
MatSample<T>::MatSample(ConstBlock &&sampleBlock)
|
||||
: MatSample(sampleBlock)
|
||||
{}
|
||||
|
||||
// assignement operator ////////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
MatSample<T> & MatSample<T>::operator=(Block &sampleBlock)
|
||||
{
|
||||
MatSample<T> tmp(sampleBlock);
|
||||
|
||||
this->swap(tmp);
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
MatSample<T> & MatSample<T>::operator=(Block &&sampleBlock)
|
||||
{
|
||||
*this = sampleBlock;
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
MatSample<T> & MatSample<T>::operator=(ConstBlock &sampleBlock)
|
||||
{
|
||||
MatSample<T> tmp(sampleBlock);
|
||||
|
||||
this->swap(tmp);
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
MatSample<T> & MatSample<T>::operator=(ConstBlock &&sampleBlock)
|
||||
{
|
||||
*this = sampleBlock;
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
// product/division by scalar operators (not provided by Eigen) ////////////////
|
||||
template <typename T>
|
||||
MatSample<T> & MatSample<T>::operator*=(const T &x)
|
||||
{
|
||||
return *this = (*this)*x;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
MatSample<T> & MatSample<T>::operator*=(const T &&x)
|
||||
{
|
||||
return *this = (*this)*x;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
MatSample<T> & MatSample<T>::operator/=(const T &x)
|
||||
{
|
||||
return *this = (*this)/x;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
MatSample<T> & MatSample<T>::operator/=(const T &&x)
|
||||
{
|
||||
return *this = (*this)/x;
|
||||
}
|
||||
|
||||
// block access ////////////////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
typename MatSample<T>::ConstBlock MatSample<T>::block(const Index i,
|
||||
const Index j,
|
||||
const Index nRow,
|
||||
const Index nCol) const
|
||||
{
|
||||
return ConstBlock(*this, i, j, nRow, nCol);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
typename MatSample<T>::Block MatSample<T>::block(const Index i,
|
||||
const Index j,
|
||||
const Index nRow,
|
||||
const Index nCol)
|
||||
{
|
||||
return Block(*this, i, j, nRow, nCol);
|
||||
}
|
||||
|
||||
// resize all matrices /////////////////////////////////////////////////////////
|
||||
template <typename T>
|
||||
void MatSample<T>::resizeMat(const Index nRow, const Index nCol)
|
||||
{
|
||||
FOR_STAT_ARRAY(*this, s)
|
||||
{
|
||||
(*this)[s].resize(nRow, nCol);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
END_LATAN_NAMESPACE
|
||||
|
||||
#endif // Latan_MatSample_hpp_
|
||||
32
lib/Statistics/StatArray.cpp
Normal file
32
lib/Statistics/StatArray.cpp
Normal file
@ -0,0 +1,32 @@
|
||||
/*
|
||||
* StatArray.cpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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/Statistics/StatArray.hpp>
|
||||
#include <LatAnalyze/includes.hpp>
|
||||
|
||||
using namespace std;
|
||||
|
||||
namespace Latan
|
||||
{
|
||||
template <>
|
||||
IoObject::IoType StatArray<Mat<double>, -1>::getType(void) const
|
||||
{
|
||||
return IoType::dMatSample;
|
||||
}
|
||||
}
|
||||
284
lib/Statistics/StatArray.hpp
Normal file
284
lib/Statistics/StatArray.hpp
Normal file
@ -0,0 +1,284 @@
|
||||
/*
|
||||
* StatArray.hpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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_StatArray_hpp_
|
||||
#define Latan_StatArray_hpp_
|
||||
|
||||
#include <LatAnalyze/Global.hpp>
|
||||
#include <LatAnalyze/Core/Mat.hpp>
|
||||
|
||||
#define FOR_STAT_ARRAY(ar, i) \
|
||||
for (Latan::Index i = -(ar).offset; i < (ar).size(); ++i)
|
||||
|
||||
BEGIN_LATAN_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* Array class with statistics *
|
||||
******************************************************************************/
|
||||
template <typename T, Index os = 0>
|
||||
class StatArray: public Array<T, dynamic, 1>, public IoObject
|
||||
{
|
||||
protected:
|
||||
typedef Array<T, dynamic, 1> Base;
|
||||
public:
|
||||
// constructors
|
||||
StatArray(void);
|
||||
explicit StatArray(const Index size);
|
||||
EIGEN_EXPR_CTOR(StatArray, unique_arg(StatArray<T, os>), Base, ArrayExpr)
|
||||
// destructor
|
||||
virtual ~StatArray(void) = default;
|
||||
// access
|
||||
Index size(void) const;
|
||||
void resize(const Index size);
|
||||
// operators
|
||||
T & operator[](const Index s);
|
||||
const T & operator[](const Index s) const;
|
||||
// statistics
|
||||
void bin(Index binSize);
|
||||
T mean(const Index pos = 0, const Index n = -1) const;
|
||||
T covariance(const StatArray<T, os> &array, const Index pos = 0,
|
||||
const Index n = -1) const;
|
||||
T covarianceMatrix(const StatArray<T, os> &array, const Index pos = 0,
|
||||
const Index n = -1) const;
|
||||
T variance(const Index pos = 0, const Index n = -1) const;
|
||||
T varianceMatrix(const Index pos = 0, const Index n = -1) const;
|
||||
T correlationMatrix(const Index pos = 0, const Index n = -1) const;
|
||||
// IO type
|
||||
virtual IoType getType(void) const;
|
||||
public:
|
||||
static constexpr Index offset = os;
|
||||
};
|
||||
|
||||
// reduction operations
|
||||
namespace ReducOp
|
||||
{
|
||||
// general templates
|
||||
template <typename T>
|
||||
inline T prod(const T &a, const T &b);
|
||||
template <typename T>
|
||||
inline T tensProd(const T &v1, const T &v2);
|
||||
template <typename T>
|
||||
inline T sum(const T &a, const T &b);
|
||||
}
|
||||
|
||||
// Sample types
|
||||
const int central = -1;
|
||||
|
||||
template <typename T>
|
||||
using Sample = StatArray<T, 1>;
|
||||
|
||||
typedef Sample<double> DSample;
|
||||
typedef Sample<std::complex<double>> CSample;
|
||||
|
||||
/******************************************************************************
|
||||
* StatArray class template implementation *
|
||||
******************************************************************************/
|
||||
// constructors ////////////////////////////////////////////////////////////////
|
||||
template <typename T, Index os>
|
||||
StatArray<T, os>::StatArray(void)
|
||||
: Base(static_cast<typename Base::Index>(os))
|
||||
{}
|
||||
|
||||
template <typename T, Index os>
|
||||
StatArray<T, os>::StatArray(const Index size)
|
||||
: Base(static_cast<typename Base::Index>(size + os))
|
||||
{}
|
||||
|
||||
// access //////////////////////////////////////////////////////////////////////
|
||||
template <typename T, Index os>
|
||||
Index StatArray<T, os>::size(void) const
|
||||
{
|
||||
return Base::size() - os;
|
||||
}
|
||||
|
||||
template <typename T, Index os>
|
||||
void StatArray<T, os>::resize(const Index size)
|
||||
{
|
||||
Base::resize(size + os);
|
||||
}
|
||||
|
||||
// operators ///////////////////////////////////////////////////////////////////
|
||||
template <typename T, Index os>
|
||||
T & StatArray<T, os>::operator[](const Index s)
|
||||
{
|
||||
return Base::operator[](s + os);
|
||||
}
|
||||
|
||||
template <typename T, Index os>
|
||||
const T & StatArray<T, os>::operator[](const Index s) const
|
||||
{
|
||||
return Base::operator[](s + os);
|
||||
}
|
||||
|
||||
|
||||
// statistics //////////////////////////////////////////////////////////////////
|
||||
template <typename T, Index os>
|
||||
void StatArray<T, os>::bin(Index binSize)
|
||||
{
|
||||
Index q = size()/binSize, r = size()%binSize;
|
||||
|
||||
for (Index i = 0; i < q; ++i)
|
||||
{
|
||||
(*this)[i] = mean(i*binSize, binSize);
|
||||
}
|
||||
if (r != 0)
|
||||
{
|
||||
(*this)[q] = mean(q*binSize, r);
|
||||
this->conservativeResize(os + q + 1);
|
||||
}
|
||||
else
|
||||
{
|
||||
this->conservativeResize(os + q);
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T, Index os>
|
||||
T StatArray<T, os>::mean(const Index pos, const Index n) const
|
||||
{
|
||||
T result = T();
|
||||
const Index m = (n >= 0) ? n : size();
|
||||
|
||||
if (m)
|
||||
{
|
||||
result = this->segment(pos+os, m).redux(&ReducOp::sum<T>);
|
||||
}
|
||||
return result/static_cast<double>(m);
|
||||
}
|
||||
|
||||
template <typename T, Index os>
|
||||
T StatArray<T, os>::covariance(const StatArray<T, os> &array, const Index pos,
|
||||
const Index n) const
|
||||
{
|
||||
T s1, s2, prs, res = T();
|
||||
const Index m = (n >= 0) ? n : size();
|
||||
|
||||
if (m)
|
||||
{
|
||||
auto arraySeg = array.segment(pos+os, m);
|
||||
auto thisSeg = this->segment(pos+os, m);
|
||||
|
||||
s1 = thisSeg.redux(&ReducOp::sum<T>);
|
||||
s2 = arraySeg.redux(&ReducOp::sum<T>);
|
||||
prs = thisSeg.binaryExpr(arraySeg, &ReducOp::prod<T>)
|
||||
.redux(&ReducOp::sum<T>);
|
||||
res = prs - ReducOp::prod(s1, s2)/static_cast<double>(m);
|
||||
}
|
||||
|
||||
return res/static_cast<double>(m - 1);
|
||||
}
|
||||
|
||||
template <typename T, Index os>
|
||||
T StatArray<T, os>::covarianceMatrix(const StatArray<T, os> &array,
|
||||
const Index pos, const Index n) const
|
||||
{
|
||||
T s1, s2, prs, res = T();
|
||||
const Index m = (n >= 0) ? n : size();
|
||||
|
||||
if (m)
|
||||
{
|
||||
auto arraySeg = array.segment(pos+os, m);
|
||||
auto thisSeg = this->segment(pos+os, m);
|
||||
|
||||
s1 = thisSeg.redux(&ReducOp::sum<T>);
|
||||
s2 = arraySeg.redux(&ReducOp::sum<T>);
|
||||
prs = thisSeg.binaryExpr(arraySeg, &ReducOp::tensProd<T>)
|
||||
.redux(&ReducOp::sum<T>);
|
||||
res = prs - ReducOp::tensProd(s1, s2)/static_cast<double>(m);
|
||||
}
|
||||
|
||||
return res/static_cast<double>(m - 1);
|
||||
}
|
||||
|
||||
template <typename T, Index os>
|
||||
T StatArray<T, os>::variance(const Index pos, const Index n) const
|
||||
{
|
||||
return covariance(*this, pos, n);
|
||||
}
|
||||
|
||||
template <typename T, Index os>
|
||||
T StatArray<T, os>::varianceMatrix(const Index pos, const Index n) const
|
||||
{
|
||||
return covarianceMatrix(*this, pos, n);
|
||||
}
|
||||
|
||||
template <typename T, Index os>
|
||||
T StatArray<T, os>::correlationMatrix(const Index pos, const Index n) const
|
||||
{
|
||||
T res = varianceMatrix(pos, n);
|
||||
T invDiag(res.rows(), 1);
|
||||
|
||||
invDiag = res.diagonal();
|
||||
invDiag = invDiag.cwiseInverse().cwiseSqrt();
|
||||
res = (invDiag*invDiag.transpose()).cwiseProduct(res);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
// reduction operations ////////////////////////////////////////////////////////
|
||||
namespace ReducOp
|
||||
{
|
||||
template <typename T>
|
||||
inline T sum(const T &a, const T &b)
|
||||
{
|
||||
return a + b;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline T prod(const T &a, const T &b)
|
||||
{
|
||||
return a*b;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
inline T tensProd(const T &v1 __dumb, const T &v2 __dumb)
|
||||
{
|
||||
LATAN_ERROR(Implementation,
|
||||
"tensorial product not implemented for this type");
|
||||
}
|
||||
|
||||
template <>
|
||||
inline Mat<double> prod(const Mat<double> &a, const Mat<double> &b)
|
||||
{
|
||||
return a.cwiseProduct(b);
|
||||
}
|
||||
|
||||
template <>
|
||||
inline Mat<double> tensProd(const Mat<double> &v1,
|
||||
const Mat<double> &v2)
|
||||
{
|
||||
if ((v1.cols() != 1) or (v2.cols() != 1))
|
||||
{
|
||||
LATAN_ERROR(Size,
|
||||
"tensorial product is only valid with column vectors");
|
||||
}
|
||||
|
||||
return v1*v2.transpose();
|
||||
}
|
||||
}
|
||||
|
||||
// IO type /////////////////////////////////////////////////////////////////////
|
||||
template <typename T, Index os>
|
||||
IoObject::IoType StatArray<T, os>::getType(void) const
|
||||
{
|
||||
return IoType::noType;
|
||||
}
|
||||
|
||||
END_LATAN_NAMESPACE
|
||||
|
||||
#endif // Latan_StatArray_hpp_
|
||||
507
lib/Statistics/XYSampleData.cpp
Normal file
507
lib/Statistics/XYSampleData.cpp
Normal file
@ -0,0 +1,507 @@
|
||||
/*
|
||||
* XYSampleData.cpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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/Statistics/XYSampleData.hpp>
|
||||
#include <LatAnalyze/includes.hpp>
|
||||
#include <LatAnalyze/Core/Math.hpp>
|
||||
|
||||
using namespace std;
|
||||
using namespace Latan;
|
||||
|
||||
/******************************************************************************
|
||||
* 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_);
|
||||
}
|
||||
|
||||
Index SampleFitResult::getNPar(void) const
|
||||
{
|
||||
return nPar_;
|
||||
}
|
||||
|
||||
double SampleFitResult::getPValue(const Index s) const
|
||||
{
|
||||
return Math::chi2PValue(getChi2(s), getNDof());
|
||||
}
|
||||
|
||||
const DoubleFunction & SampleFitResult::getModel(const Index s,
|
||||
const Index j) const
|
||||
{
|
||||
return model_[j][s];
|
||||
}
|
||||
|
||||
const DoubleFunctionSample & SampleFitResult::getModel(
|
||||
const PlaceHolder ph __dumb,
|
||||
const Index j) const
|
||||
{
|
||||
return model_[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;
|
||||
}
|
||||
|
||||
// IO //////////////////////////////////////////////////////////////////////////
|
||||
void SampleFitResult::print(const bool printXsi, ostream &out) const
|
||||
{
|
||||
char buf[256];
|
||||
Index pMax = printXsi ? size() : nPar_;
|
||||
DMat err = this->variance().cwiseSqrt();
|
||||
|
||||
sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- p-value= %.2e", getChi2(),
|
||||
static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
|
||||
out << buf << endl;
|
||||
for (Index p = 0; p < pMax; ++p)
|
||||
{
|
||||
sprintf(buf, "%8s= % e +/- %e", parName_[p].c_str(),
|
||||
(*this)[central](p), err(p));
|
||||
out << buf << endl;
|
||||
}
|
||||
}
|
||||
|
||||
/******************************************************************************
|
||||
* XYSampleData implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
XYSampleData::XYSampleData(const Index nSample)
|
||||
: nSample_(nSample)
|
||||
{}
|
||||
|
||||
// data access /////////////////////////////////////////////////////////////////
|
||||
DSample & XYSampleData::x(const Index r, const Index i)
|
||||
{
|
||||
checkXIndex(r, i);
|
||||
scheduleDataInit();
|
||||
scheduleComputeVarMat();
|
||||
if (xData_[i][r].size() == 0)
|
||||
{
|
||||
xData_[i][r].resize(nSample_);
|
||||
}
|
||||
|
||||
return xData_[i][r];
|
||||
}
|
||||
|
||||
const DSample & XYSampleData::x(const Index r, const Index i) const
|
||||
{
|
||||
checkXIndex(r, i);
|
||||
|
||||
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);
|
||||
if (!pointExists(k, j))
|
||||
{
|
||||
registerDataPoint(k, j);
|
||||
}
|
||||
scheduleDataInit();
|
||||
scheduleComputeVarMat();
|
||||
if (yData_[j][k].size() == 0)
|
||||
{
|
||||
yData_[j][k].resize(nSample_);
|
||||
}
|
||||
|
||||
return yData_[j][k];
|
||||
}
|
||||
|
||||
const DSample & XYSampleData::y(const Index k, const Index j) const
|
||||
{
|
||||
checkPoint(k, j);
|
||||
|
||||
return yData_[j].at(k);
|
||||
}
|
||||
|
||||
void XYSampleData::setUnidimData(const DMatSample &xData,
|
||||
const vector<const DMatSample *> &v)
|
||||
{
|
||||
FOR_STAT_ARRAY(xData, s)
|
||||
FOR_VEC(xData[central], r)
|
||||
{
|
||||
x(r, 0)[s] = xData[s](r);
|
||||
for (unsigned int j = 0; j < v.size(); ++j)
|
||||
{
|
||||
y(r, j)[s] = (*(v[j]))[s](r);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const DMat & XYSampleData::getXXVar(const Index i1, const Index i2)
|
||||
{
|
||||
checkXDim(i1);
|
||||
checkXDim(i2);
|
||||
computeVarMat();
|
||||
|
||||
return data_.getXXVar(i1, i2);
|
||||
}
|
||||
|
||||
const DMat & XYSampleData::getYYVar(const Index j1, const Index j2)
|
||||
{
|
||||
checkYDim(j1);
|
||||
checkYDim(j2);
|
||||
computeVarMat();
|
||||
|
||||
return data_.getYYVar(j1, j2);
|
||||
}
|
||||
|
||||
const DMat & XYSampleData::getXYVar(const Index i, const Index j)
|
||||
{
|
||||
checkXDim(i);
|
||||
checkYDim(j);
|
||||
computeVarMat();
|
||||
|
||||
return data_.getXYVar(i, j);
|
||||
}
|
||||
|
||||
DVec XYSampleData::getXError(const Index i)
|
||||
{
|
||||
checkXDim(i);
|
||||
computeVarMat();
|
||||
|
||||
return data_.getXError(i);
|
||||
}
|
||||
|
||||
DVec XYSampleData::getYError(const Index j)
|
||||
{
|
||||
checkYDim(j);
|
||||
computeVarMat();
|
||||
|
||||
return data_.getYError(j);
|
||||
}
|
||||
|
||||
// get total fit variance matrix and its pseudo-inverse ////////////////////////
|
||||
const DMat & XYSampleData::getFitVarMat(void)
|
||||
{
|
||||
computeVarMat();
|
||||
|
||||
return data_.getFitVarMat();
|
||||
}
|
||||
|
||||
const DMat & XYSampleData::getFitVarMatPInv(void)
|
||||
{
|
||||
computeVarMat();
|
||||
|
||||
return data_.getFitVarMatPInv();
|
||||
}
|
||||
|
||||
// set data to a particular sample /////////////////////////////////////////////
|
||||
void XYSampleData::setDataToSample(const Index s)
|
||||
{
|
||||
if (initData_ or (s != dataSample_))
|
||||
{
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
for (Index r = 0; r < getXSize(i); ++r)
|
||||
{
|
||||
data_.x(r, i) = xData_[i][r][s];
|
||||
}
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
for (auto &p: yData_[j])
|
||||
{
|
||||
data_.y(p.first, j) = p.second[s];
|
||||
}
|
||||
dataSample_ = s;
|
||||
initData_ = false;
|
||||
}
|
||||
}
|
||||
|
||||
// get internal XYStatData /////////////////////////////////////////////////////
|
||||
const XYStatData & XYSampleData::getData(void)
|
||||
{
|
||||
setDataToSample(central);
|
||||
computeVarMat();
|
||||
|
||||
return data_;
|
||||
}
|
||||
|
||||
// fit /////////////////////////////////////////////////////////////////////////
|
||||
SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
|
||||
const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v)
|
||||
{
|
||||
computeVarMat();
|
||||
|
||||
SampleFitResult result;
|
||||
FitResult sampleResult;
|
||||
DVec initCopy = init;
|
||||
|
||||
result.resize(nSample_);
|
||||
result.chi2_.resize(nSample_);
|
||||
result.model_.resize(v.size());
|
||||
FOR_STAT_ARRAY(result, s)
|
||||
{
|
||||
setDataToSample(s);
|
||||
if (s == central)
|
||||
{
|
||||
sampleResult = data_.fit(minimizer, initCopy, v);
|
||||
initCopy = sampleResult.segment(0, initCopy.size());
|
||||
}
|
||||
else
|
||||
{
|
||||
sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
|
||||
}
|
||||
result[s] = sampleResult;
|
||||
result.chi2_[s] = sampleResult.getChi2();
|
||||
for (unsigned int j = 0; j < v.size(); ++j)
|
||||
{
|
||||
result.model_[j].resize(nSample_);
|
||||
result.model_[j][s] = sampleResult.getModel(j);
|
||||
}
|
||||
}
|
||||
result.nPar_ = sampleResult.getNPar();
|
||||
result.nDof_ = sampleResult.nDof_;
|
||||
result.parName_ = sampleResult.parName_;
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
SampleFitResult XYSampleData::fit(Minimizer &minimizer,
|
||||
const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v)
|
||||
{
|
||||
vector<Minimizer *> mv{&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)
|
||||
{
|
||||
initData_ = true;
|
||||
}
|
||||
|
||||
// variance matrix computation /////////////////////////////////////////////////
|
||||
void XYSampleData::scheduleComputeVarMat(void)
|
||||
{
|
||||
computeVarMat_ = true;
|
||||
}
|
||||
|
||||
void XYSampleData::computeVarMat(void)
|
||||
{
|
||||
if (computeVarMat_)
|
||||
{
|
||||
// initialize data if necessary
|
||||
setDataToSample(central);
|
||||
|
||||
// compute relevant sizes
|
||||
Index size = 0, ySize = 0;
|
||||
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
size += getYSize(j);
|
||||
}
|
||||
ySize = size;
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
{
|
||||
size += getXSize(i);
|
||||
}
|
||||
|
||||
// compute total matrix
|
||||
DMatSample z(nSample_, size, 1);
|
||||
DMat var;
|
||||
Index a;
|
||||
|
||||
FOR_STAT_ARRAY(z, s)
|
||||
{
|
||||
a = 0;
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
for (auto &p: yData_[j])
|
||||
{
|
||||
z[s](a, 0) = p.second[s];
|
||||
a++;
|
||||
}
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
for (Index r = 0; r < getXSize(i); ++r)
|
||||
{
|
||||
z[s](a, 0) = xData_[i][r][s];
|
||||
a++;
|
||||
}
|
||||
}
|
||||
var = z.varianceMatrix();
|
||||
|
||||
// assign blocks to data
|
||||
Index a1, a2;
|
||||
|
||||
a1 = ySize;
|
||||
for (Index i1 = 0; i1 < getNXDim(); ++i1)
|
||||
{
|
||||
a2 = ySize;
|
||||
for (Index i2 = 0; i2 < getNXDim(); ++i2)
|
||||
{
|
||||
data_.setXXVar(i1, i2,
|
||||
var.block(a1, a2, getXSize(i1), getXSize(i2)));
|
||||
a2 += getXSize(i2);
|
||||
}
|
||||
a1 += getXSize(i1);
|
||||
}
|
||||
a1 = 0;
|
||||
for (Index j1 = 0; j1 < getNYDim(); ++j1)
|
||||
{
|
||||
a2 = 0;
|
||||
for (Index j2 = 0; j2 < getNYDim(); ++j2)
|
||||
{
|
||||
data_.setYYVar(j1, j2,
|
||||
var.block(a1, a2, getYSize(j1), getYSize(j2)));
|
||||
a2 += getYSize(j2);
|
||||
}
|
||||
a1 += getYSize(j1);
|
||||
}
|
||||
a1 = ySize;
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
{
|
||||
a2 = 0;
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
data_.setXYVar(i, j,
|
||||
var.block(a1, a2, getXSize(i), getYSize(j)));
|
||||
a2 += getYSize(j);
|
||||
}
|
||||
a1 += getXSize(i);
|
||||
}
|
||||
computeVarMat_ = false;
|
||||
}
|
||||
if (initVarMat())
|
||||
{
|
||||
data_.copyInterface(*this);
|
||||
scheduleFitVarMatInit(false);
|
||||
}
|
||||
}
|
||||
|
||||
// create data /////////////////////////////////////////////////////////////////
|
||||
void XYSampleData::createXData(const string name, const Index nData)
|
||||
{
|
||||
data_.addXDim(nData, name);
|
||||
xData_.push_back(vector<DSample>(nData));
|
||||
}
|
||||
|
||||
void XYSampleData::createYData(const string name)
|
||||
{
|
||||
data_.addYDim(name);
|
||||
yData_.push_back(map<Index, DSample>());
|
||||
}
|
||||
178
lib/Statistics/XYSampleData.hpp
Normal file
178
lib/Statistics/XYSampleData.hpp
Normal file
@ -0,0 +1,178 @@
|
||||
/*
|
||||
* XYSampleData.hpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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/Statistics/FitInterface.hpp>
|
||||
#include <LatAnalyze/Numerical/Minimizer.hpp>
|
||||
#include <LatAnalyze/Statistics/MatSample.hpp>
|
||||
#include <LatAnalyze/Functional/Model.hpp>
|
||||
#include <LatAnalyze/Statistics/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;
|
||||
Index getNPar(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;
|
||||
// IO
|
||||
void print(const bool printXsi = false,
|
||||
std::ostream &out = std::cout) const;
|
||||
private:
|
||||
DSample chi2_;
|
||||
Index nDof_{0}, nPar_{0};
|
||||
std::vector<DoubleFunctionSample> model_;
|
||||
std::vector<std::string> parName_;
|
||||
};
|
||||
|
||||
/******************************************************************************
|
||||
* XYSampleData *
|
||||
******************************************************************************/
|
||||
class XYSampleData: public FitInterface
|
||||
{
|
||||
public:
|
||||
// constructor
|
||||
explicit XYSampleData(const Index nSample);
|
||||
// destructor
|
||||
virtual ~XYSampleData(void) = default;
|
||||
// data access
|
||||
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;
|
||||
void setUnidimData(const DMatSample &xData,
|
||||
const std::vector<const DMatSample *> &v);
|
||||
template <typename... Ts>
|
||||
void setUnidimData(const DMatSample &xData,
|
||||
const Ts & ...yDatas);
|
||||
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);
|
||||
// set data to a particular sample
|
||||
void setDataToSample(const Index s);
|
||||
// get internal XYStatData
|
||||
const XYStatData & getData(void);
|
||||
// fit
|
||||
SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v);
|
||||
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v);
|
||||
template <typename... Ts>
|
||||
SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
|
||||
const DoubleModel &model, const Ts... models);
|
||||
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
|
||||
void scheduleComputeVarMat(void);
|
||||
void computeVarMat(void);
|
||||
// create data
|
||||
virtual void createXData(const std::string name, const Index nData);
|
||||
virtual void createYData(const std::string name);
|
||||
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};
|
||||
};
|
||||
|
||||
/******************************************************************************
|
||||
* XYSampleData template implementation *
|
||||
******************************************************************************/
|
||||
template <typename... Ts>
|
||||
void XYSampleData::setUnidimData(const DMatSample &xData, const Ts & ...yDatas)
|
||||
{
|
||||
static_assert(static_or<std::is_assignable<DMatSample, Ts>::value...>::value,
|
||||
"y data arguments are not compatible with DMatSample");
|
||||
|
||||
std::vector<const DMatSample *> v{&yDatas...};
|
||||
|
||||
setUnidimData(xData, v);
|
||||
}
|
||||
|
||||
template <typename... Ts>
|
||||
SampleFitResult XYSampleData::fit(std::vector<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);
|
||||
}
|
||||
|
||||
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<Minimizer *> mv{&minimizer};
|
||||
|
||||
return fit(mv, init, model, models...);
|
||||
}
|
||||
|
||||
END_LATAN_NAMESPACE
|
||||
|
||||
#endif // Latan_XYSampleData_hpp_
|
||||
610
lib/Statistics/XYStatData.cpp
Normal file
610
lib/Statistics/XYStatData.cpp
Normal file
@ -0,0 +1,610 @@
|
||||
/*
|
||||
* XYStatData.cpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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/Statistics/XYStatData.hpp>
|
||||
#include <LatAnalyze/includes.hpp>
|
||||
#include <LatAnalyze/Core/Math.hpp>
|
||||
|
||||
using namespace std;
|
||||
using namespace Latan;
|
||||
|
||||
static constexpr double maxXsiDev = 10.;
|
||||
|
||||
/******************************************************************************
|
||||
* 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_);
|
||||
}
|
||||
|
||||
Index FitResult::getNPar(void) const
|
||||
{
|
||||
return nPar_;
|
||||
}
|
||||
|
||||
double FitResult::getPValue(void) const
|
||||
{
|
||||
return Math::chi2PValue(getChi2(), getNDof());;
|
||||
}
|
||||
|
||||
const DoubleFunction & FitResult::getModel(const Index j) const
|
||||
{
|
||||
return model_[j];
|
||||
}
|
||||
|
||||
// IO //////////////////////////////////////////////////////////////////////////
|
||||
void FitResult::print(const bool printXsi, ostream &out) const
|
||||
{
|
||||
char buf[256];
|
||||
Index pMax = printXsi ? size() : nPar_;
|
||||
|
||||
sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- p-value= %.2e", getChi2(),
|
||||
static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
|
||||
out << buf << endl;
|
||||
for (Index p = 0; p < pMax; ++p)
|
||||
{
|
||||
sprintf(buf, "%8s= %e", parName_[p].c_str(), (*this)(p));
|
||||
out << buf << endl;
|
||||
}
|
||||
}
|
||||
|
||||
/******************************************************************************
|
||||
* XYStatData implementation *
|
||||
******************************************************************************/
|
||||
// data access /////////////////////////////////////////////////////////////////
|
||||
double & XYStatData::x(const Index r, const Index i)
|
||||
{
|
||||
checkXIndex(r, i);
|
||||
scheduleXMapInit();
|
||||
scheduleChi2DataVecInit();
|
||||
|
||||
return xData_[i](r);
|
||||
}
|
||||
|
||||
const double & XYStatData::x(const Index r, const Index i) const
|
||||
{
|
||||
checkXIndex(r, i);
|
||||
|
||||
return xData_[i](r);
|
||||
}
|
||||
|
||||
const DVec & XYStatData::x(const Index k) const
|
||||
{
|
||||
checkDataIndex(k);
|
||||
|
||||
updateXMap();
|
||||
|
||||
return xMap_.at(k);
|
||||
}
|
||||
|
||||
double & XYStatData::y(const Index k, const Index j)
|
||||
{
|
||||
checkYDim(j);
|
||||
if (!pointExists(k, j))
|
||||
{
|
||||
registerDataPoint(k, j);
|
||||
resizeVarMat();
|
||||
}
|
||||
scheduleXMapInit();
|
||||
scheduleChi2DataVecInit();
|
||||
|
||||
return yData_[j][k];
|
||||
}
|
||||
|
||||
const double & XYStatData::y(const Index k, const Index j) const
|
||||
{
|
||||
checkPoint(k, j);
|
||||
|
||||
return yData_[j].at(k);
|
||||
}
|
||||
|
||||
void XYStatData::setXXVar(const Index i1, const Index i2, const DMat &m)
|
||||
{
|
||||
checkXDim(i1);
|
||||
checkXDim(i2);
|
||||
checkVarMat(m, xxVar_(i1, i2));
|
||||
xxVar_(i1, i2) = m;
|
||||
if (i1 != i2)
|
||||
{
|
||||
xxVar_(i2, i1) = m.transpose();
|
||||
}
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
void XYStatData::setYYVar(const Index j1, const Index j2, const DMat &m)
|
||||
{
|
||||
checkYDim(j1);
|
||||
checkYDim(j2);
|
||||
checkVarMat(m, yyVar_(j1, j2));
|
||||
yyVar_(j1, j2) = m;
|
||||
if (j1 != j2)
|
||||
{
|
||||
yyVar_(j2, j1) = m.transpose();
|
||||
}
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
void XYStatData::setXYVar(const Index i, const Index j, const DMat &m)
|
||||
{
|
||||
checkXDim(i);
|
||||
checkYDim(j);
|
||||
checkVarMat(m, xyVar_(i, j));
|
||||
xyVar_(i, j) = m;
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
void XYStatData::setXError(const Index i, const DVec &err)
|
||||
{
|
||||
checkXDim(i);
|
||||
checkErrVec(err, xxVar_(i, i));
|
||||
xxVar_(i, i).diagonal() = err.cwiseProduct(err);
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
void XYStatData::setYError(const Index j, const DVec &err)
|
||||
{
|
||||
checkXDim(j);
|
||||
checkErrVec(err, yyVar_(j, j));
|
||||
yyVar_(j, j).diagonal() = err.cwiseProduct(err);
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
const DMat & XYStatData::getXXVar(const Index i1, const Index i2) const
|
||||
{
|
||||
checkXDim(i1);
|
||||
checkXDim(i2);
|
||||
|
||||
return xxVar_(i1, i2);
|
||||
}
|
||||
|
||||
const DMat & XYStatData::getYYVar(const Index j1, const Index j2) const
|
||||
{
|
||||
checkYDim(j1);
|
||||
checkYDim(j2);
|
||||
|
||||
return yyVar_(j1, j2);
|
||||
}
|
||||
|
||||
const DMat & XYStatData::getXYVar(const Index i, const Index j) const
|
||||
{
|
||||
checkXDim(i);
|
||||
checkYDim(j);
|
||||
|
||||
return xyVar_(i, j);
|
||||
}
|
||||
|
||||
DVec XYStatData::getXError(const Index i) const
|
||||
{
|
||||
checkXDim(i);
|
||||
|
||||
return xxVar_(i, i).diagonal().cwiseSqrt();
|
||||
}
|
||||
|
||||
DVec XYStatData::getYError(const Index j) const
|
||||
{
|
||||
checkXDim(j);
|
||||
|
||||
return yyVar_(j, j).diagonal().cwiseSqrt();
|
||||
}
|
||||
|
||||
DMat XYStatData::getTable(const Index i, const Index j) const
|
||||
{
|
||||
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)
|
||||
{
|
||||
updateFitVarMat();
|
||||
|
||||
return fitVar_;
|
||||
}
|
||||
|
||||
const DMat & XYStatData::getFitVarMatPInv(void)
|
||||
{
|
||||
updateFitVarMat();
|
||||
|
||||
return fitVarInv_;
|
||||
}
|
||||
|
||||
// fit /////////////////////////////////////////////////////////////////////////
|
||||
FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
|
||||
const vector<const DoubleModel *> &v)
|
||||
{
|
||||
// check model consistency
|
||||
checkModelVec(v);
|
||||
|
||||
// buffering
|
||||
updateLayout();
|
||||
updateFitVarMat();
|
||||
updateChi2DataVec();
|
||||
|
||||
// get number of parameters
|
||||
Index nPar = v[0]->getNPar();
|
||||
Index nXDim = getNXDim();
|
||||
Index totalNPar = nPar + layout.totalXSize;
|
||||
|
||||
// chi^2 functions
|
||||
auto corrChi2Func = [this, nPar, nXDim, totalNPar, &v](const double *x)->double
|
||||
{
|
||||
ConstMap<DVec> p(x, totalNPar);
|
||||
|
||||
updateChi2ModVec(p, v, nPar, nXDim);
|
||||
chi2Vec_ = (chi2ModVec_ - chi2DataVec_);
|
||||
|
||||
return chi2Vec_.dot(fitVarInv_*chi2Vec_);
|
||||
};
|
||||
DoubleFunction corrChi2(corrChi2Func, totalNPar);
|
||||
auto uncorrChi2Func = [this, nPar, nXDim, totalNPar, &v](const double *x)->double
|
||||
{
|
||||
ConstMap<DVec> p(x, totalNPar);
|
||||
|
||||
updateChi2ModVec(p, v, nPar, nXDim);
|
||||
chi2Vec_ = (chi2ModVec_ - chi2DataVec_);
|
||||
|
||||
return chi2Vec_.dot(chi2Vec_.cwiseQuotient(fitVar_.diagonal()));
|
||||
};
|
||||
DoubleFunction uncorrChi2(uncorrChi2Func, totalNPar);
|
||||
DoubleFunction &chi2 = hasCorrelations() ? corrChi2 : uncorrChi2;
|
||||
|
||||
for (Index p = 0; p < nPar; ++p)
|
||||
{
|
||||
chi2.varName().setName(p, v[0]->parName().getName(p));
|
||||
}
|
||||
for (Index p = 0; p < totalNPar - nPar; ++p)
|
||||
{
|
||||
chi2.varName().setName(p + nPar, "xsi_" + strFrom(p));
|
||||
}
|
||||
|
||||
// minimization
|
||||
FitResult result;
|
||||
DVec totalInit(totalNPar);
|
||||
|
||||
//// set total init vector
|
||||
totalInit.segment(0, nPar) = init;
|
||||
totalInit.segment(nPar, layout.totalXSize) =
|
||||
chi2DataVec_.segment(layout.totalYSize, layout.totalXSize);
|
||||
for (auto &m: minimizer)
|
||||
{
|
||||
m->setInit(totalInit);
|
||||
if (m->supportLimits())
|
||||
{
|
||||
//// do not allow more than maxXsiDev std. deviations on the x-axis
|
||||
for (Index p = nPar; p < totalNPar; ++p)
|
||||
{
|
||||
double err;
|
||||
|
||||
err = sqrt(fitVar_.diagonal()(layout.totalYSize + p - nPar));
|
||||
m->useLowLimit(p);
|
||||
m->useHighLimit(p);
|
||||
m->setLowLimit(p, totalInit(p) - maxXsiDev*err);
|
||||
m->setHighLimit(p, totalInit(p) + maxXsiDev*err);
|
||||
}
|
||||
}
|
||||
//// minimize and store results
|
||||
result = (*m)(chi2);
|
||||
totalInit = result;
|
||||
}
|
||||
result.chi2_ = chi2(result);
|
||||
result.nPar_ = nPar;
|
||||
result.nDof_ = layout.totalYSize - nPar;
|
||||
result.model_.resize(v.size());
|
||||
for (unsigned int j = 0; j < v.size(); ++j)
|
||||
{
|
||||
result.model_[j] = v[j]->fixPar(result);
|
||||
}
|
||||
for (Index p = 0; p < totalNPar; ++p)
|
||||
{
|
||||
result.parName_.push_back(chi2.varName().getName(p));
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
|
||||
const vector<const DoubleModel *> &v)
|
||||
{
|
||||
vector<Minimizer *> mv{&minimizer};
|
||||
|
||||
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)
|
||||
{
|
||||
xData_.push_back(DVec::Zero(nData));
|
||||
xBuf_.resize(xData_.size());
|
||||
resizeVarMat();
|
||||
}
|
||||
|
||||
void XYStatData::createYData(const std::string name __dumb)
|
||||
{
|
||||
yData_.push_back(map<Index, double>());
|
||||
resizeVarMat();
|
||||
}
|
||||
|
||||
void XYStatData::resizeVarMat(void)
|
||||
{
|
||||
xxVar_.conservativeResize(getNXDim(), getNXDim());
|
||||
for (Index i1 = 0; i1 < getNXDim(); ++i1)
|
||||
for (Index i2 = 0; i2 < getNXDim(); ++i2)
|
||||
{
|
||||
xxVar_(i1, i2).conservativeResize(getXSize(i1), getXSize(i2));
|
||||
}
|
||||
yyVar_.conservativeResize(getNYDim(), getNYDim());
|
||||
for (Index j1 = 0; j1 < getNYDim(); ++j1)
|
||||
for (Index j2 = 0; j2 < getNYDim(); ++j2)
|
||||
{
|
||||
yyVar_(j1, j2).conservativeResize(getYSize(j1), getYSize(j2));
|
||||
}
|
||||
xyVar_.conservativeResize(getNXDim(), getNYDim());
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
for (Index j = 0; j < getNYDim(); ++j)
|
||||
{
|
||||
xyVar_(i, j).conservativeResize(getXSize(i), getYSize(j));
|
||||
}
|
||||
scheduleFitVarMatInit();
|
||||
}
|
||||
|
||||
// schedule buffer computation /////////////////////////////////////////////////
|
||||
void XYStatData::scheduleXMapInit(void)
|
||||
{
|
||||
initXMap_ = true;
|
||||
}
|
||||
|
||||
void XYStatData::scheduleChi2DataVecInit(void)
|
||||
{
|
||||
initChi2DataVec_ = true;
|
||||
}
|
||||
|
||||
// buffer total fit variance matrix ////////////////////////////////////////////
|
||||
void XYStatData::updateFitVarMat(void)
|
||||
{
|
||||
if (initVarMat())
|
||||
{
|
||||
updateLayout();
|
||||
|
||||
DMat &v = fitVar_;
|
||||
Index roffs, coffs;
|
||||
|
||||
v.resize(layout.totalSize, layout.totalSize);
|
||||
roffs = layout.totalYSize;
|
||||
for (Index ifit1 = 0; ifit1 < layout.nXFitDim; ++ifit1)
|
||||
{
|
||||
coffs = layout.totalYSize;
|
||||
for (Index ifit2 = 0; ifit2 < layout.nXFitDim; ++ifit2)
|
||||
{
|
||||
for (Index rfit1 = 0; rfit1 < layout.xSize[ifit1]; ++rfit1)
|
||||
for (Index rfit2 = 0; rfit2 < layout.xSize[ifit2]; ++rfit2)
|
||||
{
|
||||
Index i1, i2, r1, r2;
|
||||
|
||||
i1 = layout.xDim[ifit1];
|
||||
i2 = layout.xDim[ifit2];
|
||||
r1 = layout.x[ifit1][rfit1];
|
||||
r2 = layout.x[ifit2][rfit2];
|
||||
|
||||
v(roffs+rfit1, coffs+rfit2) = xxVar_(i1, i2)(r1, r2);
|
||||
v(coffs+rfit2, roffs+rfit1) = v(roffs+rfit1, coffs+rfit2);
|
||||
}
|
||||
coffs += layout.xSize[ifit2];
|
||||
}
|
||||
roffs += layout.xSize[ifit1];
|
||||
}
|
||||
roffs = 0;
|
||||
for (Index jfit1 = 0; jfit1 < layout.nYFitDim; ++jfit1)
|
||||
{
|
||||
coffs = 0;
|
||||
for (Index jfit2 = 0; jfit2 < layout.nYFitDim; ++jfit2)
|
||||
{
|
||||
for (Index sfit1 = 0; sfit1 < layout.ySize[jfit1]; ++sfit1)
|
||||
for (Index sfit2 = 0; sfit2 < layout.ySize[jfit2]; ++sfit2)
|
||||
{
|
||||
Index j1, j2, s1, s2;
|
||||
|
||||
j1 = layout.yDim[jfit1];
|
||||
j2 = layout.yDim[jfit2];
|
||||
s1 = layout.y[jfit1][sfit1];
|
||||
s2 = layout.y[jfit2][sfit2];
|
||||
|
||||
v(roffs+sfit1, coffs+sfit2) = yyVar_(j1, j2)(s1, s2);
|
||||
v(coffs+sfit2, roffs+sfit1) = v(roffs+sfit1, coffs+sfit2);
|
||||
}
|
||||
coffs += layout.ySize[jfit2];
|
||||
}
|
||||
roffs += layout.ySize[jfit1];
|
||||
}
|
||||
roffs = layout.totalYSize;
|
||||
for (Index ifit = 0; ifit < layout.nXFitDim; ++ifit)
|
||||
{
|
||||
coffs = 0;
|
||||
for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
|
||||
{
|
||||
for (Index rfit = 0; rfit < layout.xSize[ifit]; ++rfit)
|
||||
for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
|
||||
{
|
||||
Index i, j, r, s;
|
||||
|
||||
i = layout.xDim[ifit];
|
||||
j = layout.yDim[jfit];
|
||||
r = layout.x[ifit][rfit];
|
||||
s = layout.y[jfit][sfit];
|
||||
|
||||
v(roffs+rfit, coffs+sfit) = xyVar_(i, j)(r, s);
|
||||
v(coffs+sfit, roffs+rfit) = v(roffs+rfit, coffs+sfit);
|
||||
}
|
||||
coffs += layout.ySize[jfit];
|
||||
}
|
||||
roffs += layout.xSize[ifit];
|
||||
}
|
||||
chi2DataVec_.resize(layout.totalSize);
|
||||
chi2ModVec_.resize(layout.totalSize);
|
||||
chi2Vec_.resize(layout.totalSize);
|
||||
fitVar_ = fitVar_.cwiseProduct(makeCorrFilter());
|
||||
fitVarInv_ = fitVar_.pInverse(getSvdTolerance());
|
||||
scheduleFitVarMatInit(false);
|
||||
}
|
||||
}
|
||||
|
||||
// buffer list of x vectors ////////////////////////////////////////////////////
|
||||
void XYStatData::updateXMap(void) const
|
||||
{
|
||||
if (initXMap_)
|
||||
{
|
||||
XYStatData * modThis = const_cast<XYStatData *>(this);
|
||||
|
||||
modThis->xMap_.clear();
|
||||
modThis->xMap_.resize(getMaxDataIndex());
|
||||
for (auto k: getDataIndexSet())
|
||||
{
|
||||
modThis->xMap_[k] = DVec(getNXDim());
|
||||
for (Index i = 0; i < getNXDim(); ++i)
|
||||
{
|
||||
modThis->xMap_[k](i) = xData_[i](dataCoord(k)[i]);
|
||||
}
|
||||
}
|
||||
modThis->initXMap_ = false;
|
||||
}
|
||||
}
|
||||
|
||||
// buffer chi^2 vectors ////////////////////////////////////////////////////////
|
||||
void XYStatData::updateChi2DataVec(void)
|
||||
{
|
||||
if (initChi2DataVec_)
|
||||
{
|
||||
Index a = 0, j, k, i, r;
|
||||
|
||||
updateLayout();
|
||||
for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
|
||||
for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
|
||||
{
|
||||
j = layout.yDim[jfit];
|
||||
k = layout.data[jfit][sfit];
|
||||
chi2DataVec_(a) = yData_[j][k];
|
||||
a++;
|
||||
}
|
||||
for (Index ifit = 0; ifit < layout.nXFitDim; ++ifit)
|
||||
for (Index rfit = 0; rfit < layout.xSize[ifit]; ++rfit)
|
||||
{
|
||||
i = layout.xDim[ifit];
|
||||
r = layout.x[ifit][rfit];
|
||||
chi2DataVec_(a) = xData_[i](r);
|
||||
a++;
|
||||
}
|
||||
initChi2DataVec_ = false;
|
||||
}
|
||||
}
|
||||
|
||||
// WARNING: updateChi2ModVec is heavily called by fit
|
||||
void XYStatData::updateChi2ModVec(const DVec p,
|
||||
const vector<const DoubleModel *> &v,
|
||||
const Index nPar, const Index nXDim)
|
||||
{
|
||||
updateLayout();
|
||||
updateXMap();
|
||||
|
||||
Index a = 0, j, k, ind;
|
||||
auto &par = p.segment(0, nPar), &xsi = p.segment(nPar, layout.totalXSize);
|
||||
|
||||
for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
|
||||
{
|
||||
j = layout.yDim[jfit];
|
||||
for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
|
||||
{
|
||||
|
||||
k = layout.data[jfit][sfit];
|
||||
for (Index i = 0; i < nXDim; ++i)
|
||||
{
|
||||
ind = layout.xIndFromData[k][i] - layout.totalYSize;
|
||||
xBuf_(i) = (ind >= 0) ? xsi(ind) : xMap_[k](i);
|
||||
}
|
||||
chi2ModVec_(a) = (*v[j])(xBuf_.data(), par.data());
|
||||
a++;
|
||||
}
|
||||
}
|
||||
chi2ModVec_.segment(a, layout.totalXSize) = xsi;
|
||||
}
|
||||
236
lib/Statistics/XYStatData.hpp
Normal file
236
lib/Statistics/XYStatData.hpp
Normal file
@ -0,0 +1,236 @@
|
||||
/*
|
||||
* XYStatData.hpp, part of LatAnalyze 3
|
||||
*
|
||||
* Copyright (C) 2013 - 2016 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_XYStatData_hpp_
|
||||
#define Latan_XYStatData_hpp_
|
||||
|
||||
#include <LatAnalyze/Global.hpp>
|
||||
#include <LatAnalyze/Statistics/FitInterface.hpp>
|
||||
#include <LatAnalyze/Numerical/Minimizer.hpp>
|
||||
#include <LatAnalyze/Functional/Model.hpp>
|
||||
|
||||
BEGIN_LATAN_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* object for fit result *
|
||||
******************************************************************************/
|
||||
class FitResult: public DVec
|
||||
{
|
||||
friend class XYStatData;
|
||||
friend class XYSampleData;
|
||||
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;
|
||||
Index getNPar(void) const;
|
||||
double getPValue(void) const;
|
||||
const DoubleFunction & getModel(const Index j = 0) const;
|
||||
// IO
|
||||
void print(const bool printXsi = false,
|
||||
std::ostream &out = std::cout) const;
|
||||
private:
|
||||
double chi2_{0.};
|
||||
Index nDof_{0}, nPar_{0};
|
||||
std::vector<DoubleFunction> model_;
|
||||
std::vector<std::string> parName_;
|
||||
};
|
||||
|
||||
/******************************************************************************
|
||||
* class for X vs. Y statistical data *
|
||||
******************************************************************************/
|
||||
class XYStatData: public FitInterface
|
||||
{
|
||||
public:
|
||||
// constructor
|
||||
XYStatData(void) = default;
|
||||
// destructor
|
||||
virtual ~XYStatData(void) = default;
|
||||
// data access
|
||||
double & x(const Index r, const Index i);
|
||||
const double & x(const Index r, const Index i) const;
|
||||
const DVec & x(const Index k) const;
|
||||
double & y(const Index k, const Index j);
|
||||
const double & y(const Index k, const Index j) 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);
|
||||
void setXError(const Index i, const DVec &err);
|
||||
void setYError(const Index j, const DVec &err);
|
||||
template <typename... Ts>
|
||||
void setUnidimData(const DMat &xData, const Ts & ...yDatas);
|
||||
const DMat & getXXVar(const Index i1, const Index i2) const;
|
||||
const DMat & getYYVar(const Index j1, const Index j2) const;
|
||||
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) const;
|
||||
// get total fit variance matrix and its pseudo-inverse
|
||||
const DMat & getFitVarMat(void);
|
||||
const DMat & getFitVarMatPInv(void);
|
||||
// fit
|
||||
FitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v);
|
||||
FitResult fit(Minimizer &minimizer, const DVec &init,
|
||||
const std::vector<const DoubleModel *> &v);
|
||||
template <typename... Ts>
|
||||
FitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
|
||||
const DoubleModel &model, const Ts... models);
|
||||
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);
|
||||
virtual void createYData(const std::string name);
|
||||
void resizeVarMat(void);
|
||||
private:
|
||||
// schedule buffer computation
|
||||
void scheduleXMapInit(void);
|
||||
void scheduleChi2DataVecInit(void);
|
||||
// buffer total fit variance matrix
|
||||
void updateFitVarMat(void);
|
||||
// buffer list of x vectors
|
||||
void updateXMap(void) const;
|
||||
// buffer chi^2 vectors
|
||||
void updateChi2DataVec(void);
|
||||
void updateChi2ModVec(const DVec p,
|
||||
const std::vector<const DoubleModel *> &v,
|
||||
const Index nPar, const Index nXDim);
|
||||
private:
|
||||
std::vector<std::map<Index, double>> yData_;
|
||||
// no map here for fit performance
|
||||
std::vector<DVec> xData_;
|
||||
std::vector<DVec> xMap_;
|
||||
Mat<DMat> xxVar_, yyVar_, xyVar_;
|
||||
DMat fitVar_, fitVarInv_;
|
||||
DVec chi2DataVec_, chi2ModVec_, chi2Vec_;
|
||||
DVec xBuf_;
|
||||
bool initXMap_{true};
|
||||
bool initChi2DataVec_{true};
|
||||
};
|
||||
|
||||
/******************************************************************************
|
||||
* XYStatData template implementation *
|
||||
******************************************************************************/
|
||||
template <typename... Ts>
|
||||
void XYStatData::setUnidimData(const DMat &xData, const Ts & ...yDatas)
|
||||
{
|
||||
static_assert(static_or<std::is_assignable<DMat, Ts>::value...>::value,
|
||||
"y data arguments are not compatible with DMat");
|
||||
|
||||
std::vector<const DMat *> yData{&yDatas...};
|
||||
|
||||
FOR_VEC(xData, r)
|
||||
{
|
||||
x(r, 0) = xData(r);
|
||||
for (unsigned int j = 0; j < yData.size(); ++j)
|
||||
{
|
||||
y(r, j) = (*(yData[j]))(r);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename... Ts>
|
||||
FitResult XYStatData::fit(std::vector<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);
|
||||
}
|
||||
|
||||
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<Minimizer *> mv{&minimizer};
|
||||
|
||||
return fit(mv, init, model, models...);
|
||||
}
|
||||
|
||||
/******************************************************************************
|
||||
* error check macros *
|
||||
******************************************************************************/
|
||||
#define checkVarMat(m, var)\
|
||||
if (((m).rows() != (var).rows()) or ((m).cols() != (var).cols()))\
|
||||
{\
|
||||
LATAN_ERROR(Size, "provided variance matrix has a wrong size"\
|
||||
" (expected " + strFrom((var).rows()) + "x"\
|
||||
+ strFrom((var).cols()) + ", got " + strFrom((m).rows())\
|
||||
+ "x" + strFrom((m).cols()) + ")");\
|
||||
}
|
||||
|
||||
|
||||
#define checkErrVec(err, var)\
|
||||
if ((err).size() != (var).rows())\
|
||||
{\
|
||||
LATAN_ERROR(Size, "provided error vector has a wrong size"\
|
||||
" (expected " + strFrom((var).rows()) + ", got "\
|
||||
+ strFrom((err).size()) + ")");\
|
||||
}
|
||||
|
||||
#define checkModelVec(v)\
|
||||
if (static_cast<Index>((v).size()) != getNYDim())\
|
||||
{\
|
||||
LATAN_ERROR(Size, "provided model vector has a wrong size"\
|
||||
" (expected " + strFrom(getNYDim()) + ", got "\
|
||||
+ strFrom((v).size()) + ")");\
|
||||
}\
|
||||
for (unsigned int _i = 1; _i < (v).size(); ++_i)\
|
||||
{\
|
||||
if ((v)[_i]->getNArg() != getNXDim())\
|
||||
{\
|
||||
LATAN_ERROR(Size, "model " + strFrom(_i) + " has a wrong"\
|
||||
+ " number of argument (expected " + strFrom(getNXDim())\
|
||||
+ ", got " + strFrom((v)[_i]->getNArg()));\
|
||||
}\
|
||||
}\
|
||||
{\
|
||||
Index _nPar = (v)[0]->getNPar();\
|
||||
for (unsigned int _i = 1; _i < (v).size(); ++_i)\
|
||||
{\
|
||||
if ((v)[_i]->getNPar() != _nPar)\
|
||||
{\
|
||||
LATAN_ERROR(Size, "model " + strFrom(_i) + " has a wrong"\
|
||||
+ " number of parameter (expected " + strFrom(_nPar)\
|
||||
+ ", got " + strFrom((v)[_i]->getNPar()));\
|
||||
}\
|
||||
}\
|
||||
}
|
||||
|
||||
END_LATAN_NAMESPACE
|
||||
|
||||
#endif // Latan_XYStatData_hpp_
|
||||
Reference in New Issue
Block a user