From f826f30e82a7eef4bd3be44b22ec28d7db626d75 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 23 Jan 2020 19:11:01 +0000 Subject: [PATCH 1/5] first stab at correlator fit utility classes --- lib/Makefile.am | 4 + lib/Physics/CorrelatorFitter.cpp | 331 +++++++++++++++++++++++++++++++ lib/Physics/CorrelatorFitter.hpp | 79 ++++++++ lib/Physics/EffectiveMass.cpp | 103 ++++++++++ lib/Physics/EffectiveMass.hpp | 49 +++++ 5 files changed, 566 insertions(+) create mode 100644 lib/Physics/CorrelatorFitter.cpp create mode 100644 lib/Physics/CorrelatorFitter.hpp create mode 100644 lib/Physics/EffectiveMass.cpp create mode 100644 lib/Physics/EffectiveMass.hpp diff --git a/lib/Makefile.am b/lib/Makefile.am index 48c9401..ddb23b6 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -54,6 +54,8 @@ libLatAnalyze_la_SOURCES = \ Numerical/Minimizer.cpp \ Numerical/RootFinder.cpp \ Numerical/Solver.cpp \ + Physics/CorrelatorFitter.cpp \ + Physics/EffectiveMass.cpp \ Statistics/FitInterface.cpp \ Statistics/Histogram.cpp \ Statistics/Random.cpp \ @@ -97,6 +99,8 @@ HPPFILES = \ Numerical/Minimizer.hpp \ Numerical/RootFinder.hpp \ Numerical/Solver.hpp \ + Physics/CorrelatorFitter.hpp \ + Physics/EffectiveMass.hpp \ Statistics/Dataset.hpp \ Statistics/FitInterface.hpp \ Statistics/Histogram.hpp \ diff --git a/lib/Physics/CorrelatorFitter.cpp b/lib/Physics/CorrelatorFitter.cpp new file mode 100644 index 0000000..6677978 --- /dev/null +++ b/lib/Physics/CorrelatorFitter.cpp @@ -0,0 +1,331 @@ +/* + * CorrelatorFitter.cpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2020 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 . + */ + +#include +#include + +using namespace std; +using namespace Latan; + +DoubleModel CorrelatorModels::makeExpModel(const Index nState) +{ + DoubleModel mod; + + mod.setFunction([nState](const double *x, const double *p) + { + double res = 0.; + + for (unsigned int i = 0; i < nState; ++i) + { + res += p[2*i + 1]*exp(-p[2*i]*x[0]); + } + + return res; + }, 1, 2*nState); + for (unsigned int i = 0; i < nState; ++i) + { + mod.parName().setName(2*i, "E_" + strFrom(i)); + mod.parName().setName(2*i + 1, "Z_" + strFrom(i)); + } + + return mod; +} + +DoubleModel CorrelatorModels::makeCoshModel(const Index nState, const Index nt) +{ + DoubleModel mod; + + mod.setFunction([nState, nt](const double *x, const double *p) + { + double res = 0.; + + for (unsigned int i = 0; i < nState; ++i) + { + res += p[2*i + 1]*(exp(-p[2*i]*x[0]) + exp(-p[2*i]*(nt - x[0]))); + } + + return res; + }, 1, 2*nState); + for (unsigned int i = 0; i < nState; ++i) + { + mod.parName().setName(2*i, "E_" + strFrom(i)); + mod.parName().setName(2*i + 1, "Z_" + strFrom(i)); + } + + return mod; +} + +DoubleModel CorrelatorModels::makeSinhModel(const Index nState, const Index nt) +{ + DoubleModel mod; + + mod.setFunction([nState, nt](const double *x, const double *p) + { + double res = 0.; + + for (unsigned int i = 0; i < nState; ++i) + { + res += p[2*i + 1]*(exp(-p[2*i]*x[0]) - exp(-p[2*i]*(nt - x[0]))); + } + + return res; + }, 1, 2*nState); + for (unsigned int i = 0; i < nState; ++i) + { + mod.parName().setName(2*i, "E_" + strFrom(i)); + mod.parName().setName(2*i + 1, "Z_" + strFrom(i)); + } + + return mod; +} + +DoubleModel CorrelatorModels::makeConstModel(void) +{ + DoubleModel mod; + + mod.setFunction([](const double *x, const double *p __dumb) + { + return x[0]; + }, 1, 1); + mod.parName().setName(0, "cst"); + + return mod; +} + +DoubleModel CorrelatorModels::makeLinearModel(void) +{ + DoubleModel mod; + + mod.setFunction([](const double *x, const double *p) + { + return p[0] + p[1]*x[0]; + }, 1, 2); + + return mod; +} + +CorrelatorModels::ModelPar CorrelatorModels::parseModel(const string s) +{ + smatch sm; + ModelPar par; + + if (regex_match(s, sm, regex("exp([0-9]+)"))) + { + par.type = CorrelatorType::exp; + par.nState = strTo(sm[0].str()); + } + else if (regex_match(s, sm, regex("cosh([0-9]+)"))) + { + par.type = CorrelatorType::cosh; + par.nState = strTo(sm[0].str()); + } + else if (regex_match(s, sm, regex("sinh([0-9]+)"))) + { + par.type = CorrelatorType::sinh; + par.nState = strTo(sm[0].str()); + } + else if (s == "linear") + { + par.type = CorrelatorType::linear; + par.nState = 1; + } + else if (s == "cst") + { + par.type = CorrelatorType::cst; + par.nState = 1; + } + else + { + par.type = CorrelatorType::undefined; + par.nState = 0; + } + + return par; +} + +DoubleModel CorrelatorModels::makeModel(const CorrelatorModels::ModelPar par, + const Index nt) +{ + switch (par.type) + { + case CorrelatorType::undefined: + LATAN_ERROR(Argument, "correlator type is undefined"); + break; + case CorrelatorType::exp: + return makeExpModel(par.nState); + break; + case CorrelatorType::cosh: + return makeCoshModel(par.nState, nt); + break; + case CorrelatorType::sinh: + return makeSinhModel(par.nState, nt); + break; + case CorrelatorType::linear: + return makeLinearModel(); + break; + case CorrelatorType::cst: + return makeConstModel(); + break; + } +} + +DVec CorrelatorModels::parameterGuess(const DMatSample &corr, + const ModelPar par) +{ + DVec init; + Index nt = corr.size(); + + switch (par.type) + { + case CorrelatorType::undefined: + LATAN_ERROR(Argument, "correlator type is undefined"); + break; + case CorrelatorType::exp: + case CorrelatorType::cosh: + case CorrelatorType::sinh: + init.resize(2*par.nState); + init(0) = log(corr[central](nt/4)/corr[central](nt/4 + 1)); + init(1) = corr[central](nt/4)/(exp(-init(0)*nt/4)); + for (Index p = 2; p < init.size(); p += 2) + { + init(p) = 2*init(p - 2); + init(p + 1) = init(p - 1)/2.; + } + break; + case CorrelatorType::linear: + init.resize(2); + init(0) = corr[central](nt/4) - corr[central](nt/4 + 1, 0); + init(1) = corr[central](nt/4, 0) + nt/4*init(0); + break; + case CorrelatorType::cst: + init.resize(1); + init(0) = corr[central](nt/4); + break; + default: + break; + } + + return init; +} + +CorrelatorFitter::CorrelatorFitter(const DMatSample &corr) +{ + setCorrelator(corr); +} + +CorrelatorFitter::CorrelatorFitter(const std::vector &corr) +{ + setCorrelators(corr); +} + +void CorrelatorFitter::setCorrelator(const DMatSample &corr) +{ + std::vector vec; + + vec.push_back(corr); + setCorrelators(vec); +} + +void CorrelatorFitter::setCorrelators(const std::vector &corr) +{ + Index nSample = corr[0].size(); + DMatSample tVec(nSample); + std::vector ptVec; + + nt_ = corr[0][central].rows(); + tVec.fill(DVec::LinSpaced(nt_, 0, nt_ - 1)); + for (auto &c: corr) + { + ptVec.push_back(&c); + } + data_.reset(new XYSampleData(corr[0].size())); + data_->addXDim(nt_, "t/a", true); + for (unsigned int i = 0; i < corr.size(); ++i) + { + data_->addYDim("C_" + strFrom(i) + "(t)"); + } + data_->setUnidimData(tVec, ptVec); + model_.resize(corr.size()); + range_.resize(corr.size(), make_pair(0, nt_ - 1)); + thinning_.resize(corr.size(), 1); +} + +void CorrelatorFitter::setModel(const DoubleModel &model, const Index i) +{ + model_[i] = model; +} + +const DoubleModel & CorrelatorFitter::getModel(const Index i) const +{ + return model_.at(i); +} + +void CorrelatorFitter::setFitRange(const Index tMin, const Index tMax, + const Index i) +{ + range_[i] = make_pair(tMin, tMax); + refreshRanges(); +} + +void CorrelatorFitter::setCorrelation(const bool isCorrelated, const Index i, + const Index j) +{ + data_->assumeYYCorrelated(isCorrelated, i, j); +} + +DMat CorrelatorFitter::getVarianceMatrix(void) +{ + return data_->getFitVarMat(); +} + +void CorrelatorFitter::setThinning(const Index thinning, const Index i) +{ + thinning_[i] = thinning; + refreshRanges(); +} + +SampleFitResult CorrelatorFitter::fit(Minimizer &minimizer, const DVec &init) +{ + vector vecPt = {&minimizer}; + + return fit(vecPt, init); +} + +SampleFitResult CorrelatorFitter::fit(vector &minimizer, + const DVec &init) +{ + vector vecPt(model_.size()); + + for (unsigned int i = 0; i < model_.size(); ++i) + { + vecPt[i] = &(model_[i]); + } + + return data_->fit(minimizer, init, vecPt); +} + +void CorrelatorFitter::refreshRanges(void) +{ + for (unsigned int i = 0; i < range_.size(); ++i) + for (Index t = 0; t < nt_; ++t) + { + data_->fitPoint((t >= range_[i].first) and (t <= range_[i].second) + and ((t - range_[i].first) % thinning_[i] == 0), t); + } +} diff --git a/lib/Physics/CorrelatorFitter.hpp b/lib/Physics/CorrelatorFitter.hpp new file mode 100644 index 0000000..fe4071e --- /dev/null +++ b/lib/Physics/CorrelatorFitter.hpp @@ -0,0 +1,79 @@ +/* + * CorrelatorFitter.hpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2020 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 . + */ + +#ifndef Latan_CorrelatorFitter_hpp_ +#define Latan_CorrelatorFitter_hpp_ + +#include +#include +#include + +BEGIN_LATAN_NAMESPACE + +enum class CorrelatorType {undefined, exp, cosh, sinh, linear, cst}; + +namespace CorrelatorModels +{ + struct ModelPar + { + CorrelatorType type; + Index nState; + }; + + DoubleModel makeExpModel(const Index nState); + DoubleModel makeCoshModel(const Index nState, const Index nt); + DoubleModel makeSinhModel(const Index nState, const Index nt); + DoubleModel makeConstModel(void); + DoubleModel makeLinearModel(void); + ModelPar parseModel(const std::string s); + DoubleModel makeModel(const ModelPar par, const Index nt); + DVec parameterGuess(const DMatSample &corr, const ModelPar par); +}; + +class CorrelatorFitter +{ +public: + CorrelatorFitter(const DMatSample &corr); + CorrelatorFitter(const std::vector &corr); + void setCorrelator(const DMatSample &corr); + void setCorrelators(const std::vector &corr); + const DMatSample & getCorrelator(const Index i = 0) const; + const std::vector & getCorrelators(void) const; + void setModel(const DoubleModel &model, const Index i = 0); + const DoubleModel & getModel(const Index i = 0) const; + void setFitRange(const Index tMin, const Index tMax, const Index i = 0); + void setCorrelation(const bool isCorrelated, const Index i = 0, + const Index j = 0); + DMat getVarianceMatrix(void); + void setThinning(const Index thinning, const Index i = 0); + SampleFitResult fit(Minimizer &minimizer, const DVec &init); + SampleFitResult fit(std::vector &minimizer, const DVec &init); +private: + void refreshRanges(void); +private: + Index nt_{0}; + std::unique_ptr data_; + std::vector model_; + std::vector> range_; + std::vector thinning_; +}; + +END_LATAN_NAMESPACE + +#endif // Latan_CorrelatorFitter_hpp_ \ No newline at end of file diff --git a/lib/Physics/EffectiveMass.cpp b/lib/Physics/EffectiveMass.cpp new file mode 100644 index 0000000..672643f --- /dev/null +++ b/lib/Physics/EffectiveMass.cpp @@ -0,0 +1,103 @@ +/* + * EffectiveMass.cpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2020 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 . + */ + +#include +#include + +using namespace std; +using namespace Latan; + +// constructors //////////////////////////////////////////////////////////////// +EffectiveMass::EffectiveMass(const CorrelatorType type) +: type_(type) +{} + +// access ////////////////////////////////////////////////////////////////////// +CorrelatorType EffectiveMass::getType(void) const +{ + return type_; +} + +void EffectiveMass::setType(const CorrelatorType type) +{ + type_ = type; +} + +// compute effective mass ////////////////////////////////////////////////////// +DVec EffectiveMass::operator()(const DVec &corr) const +{ + Index nt = corr.size(); + DVec em; + + if (nt < 2) + { + LATAN_ERROR(Size, "input vector has less than 2 elements"); + } + switch (type_) + { + case CorrelatorType::undefined: + LATAN_ERROR(Argument, "correlator type is undefined"); + break; + case CorrelatorType::exp: + em.resize(nt - 1); + for (Index t = 1; t < nt; ++t) + { + em(t - 1) = log(corr(t - 1)/corr(t)); + } + break; + case CorrelatorType::cosh: + em.resize(nt - 2); + for (Index t = 1; t < nt - 1; ++t) + { + em(t - 1) = acosh((corr(t - 1) + corr(t + 1))/(2.*corr(t))); + } + break; + case CorrelatorType::sinh: + em.resize(nt - 2); + for (Index t = 1; t < nt - 1; ++t) + { + em(t - 1) = acosh((corr(t - 1) + corr(t + 1))/(2.*corr(t))); + } + break; + case CorrelatorType::linear: + em.resize(nt - 1); + for (Index t = 0; t < nt - 1; ++t) + { + em(t) = corr(t) - corr(t - 1); + } + break; + case CorrelatorType::cst: + em = corr; + break; + } + + return em; +} + +DMatSample EffectiveMass::operator()(const DMatSample &corr) const +{ + DMatSample em(corr.size()); + + FOR_STAT_ARRAY(em, s) + { + em[s] = (*this)(corr[s]); + } + + return em; +} diff --git a/lib/Physics/EffectiveMass.hpp b/lib/Physics/EffectiveMass.hpp new file mode 100644 index 0000000..500dab7 --- /dev/null +++ b/lib/Physics/EffectiveMass.hpp @@ -0,0 +1,49 @@ +/* + * EffectiveMass.hpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2020 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 . + */ + +#ifndef Latan_EffectiveMass_hpp_ +#define Latan_EffectiveMass_hpp_ + +#include +#include +#include + +BEGIN_LATAN_NAMESPACE + +/****************************************************************************** + * Effective mass class * + ******************************************************************************/ +class EffectiveMass +{ +public: + // constructors + EffectiveMass(const CorrelatorType type = CorrelatorType::exp); + // access + CorrelatorType getType(void) const; + void setType(const CorrelatorType type); + // compute effective mass + DVec operator()(const DVec &corr) const; + DMatSample operator()(const DMatSample &corr) const; +private: + CorrelatorType type_; +}; + +END_LATAN_NAMESPACE + +#endif // Latan_EffectiveMass_hpp_ \ No newline at end of file From 1bde8822b1b9b10447229f67aa644a9fee5434eb Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 28 Jan 2020 17:34:00 +0000 Subject: [PATCH 2/5] band plot does not rely anymore on tac & tail -r commands --- lib/Core/Plot.cpp | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/lib/Core/Plot.cpp b/lib/Core/Plot.cpp index 25d9286..fd695f7 100644 --- a/lib/Core/Plot.cpp +++ b/lib/Core/Plot.cpp @@ -236,14 +236,21 @@ PlotFunction::PlotFunction(const DoubleFunction &function, const double xMin, // PlotPredBand constructor //////////////////////////////////////////////////// void PlotPredBand::makePredBand(const DMat &low, const DMat &high, const double opacity) { - string lowFileName, highFileName; + string lowFileName, highFileName, contFileName; + DMat contour(low.rows() + high.rows() + 1, 2); - lowFileName = dumpToTmpFile(low); - highFileName = dumpToTmpFile(high); - pushTmpFile(lowFileName); - pushTmpFile(highFileName); - setCommand("'< (cat " + lowFileName + "; tac " + highFileName + - "; head -n1 " + lowFileName + ")' u 1:2 w filledcurves closed" + + FOR_MAT(low, i, j) + { + contour(i, j) = low(i, j); + } + FOR_MAT(high, i, j) + { + contour(low.rows() + i, j) = high(high.rows() - i - 1, j); + } + contour.row(low.rows() + high.rows()) = low.row(0); + contFileName = dumpToTmpFile(contour); + pushTmpFile(contFileName); + setCommand("'" + contFileName + "' u 1:2 w filledcurves closed" + " fs solid " + strFrom(opacity) + " noborder"); } From 685d4330329c9ef3d07cefad1335bcd6acb5d33f Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 28 Jan 2020 17:34:37 +0000 Subject: [PATCH 3/5] minor fixes in new physics classes --- lib/Physics/CorrelatorFitter.cpp | 13 +++++++++---- lib/Physics/CorrelatorFitter.hpp | 2 ++ lib/Physics/EffectiveMass.cpp | 30 ++++++++++++++++++++++++++++-- lib/Physics/EffectiveMass.hpp | 1 + 4 files changed, 40 insertions(+), 6 deletions(-) diff --git a/lib/Physics/CorrelatorFitter.cpp b/lib/Physics/CorrelatorFitter.cpp index 6677978..d48699f 100644 --- a/lib/Physics/CorrelatorFitter.cpp +++ b/lib/Physics/CorrelatorFitter.cpp @@ -128,17 +128,17 @@ CorrelatorModels::ModelPar CorrelatorModels::parseModel(const string s) if (regex_match(s, sm, regex("exp([0-9]+)"))) { par.type = CorrelatorType::exp; - par.nState = strTo(sm[0].str()); + par.nState = strTo(sm[1].str()); } else if (regex_match(s, sm, regex("cosh([0-9]+)"))) { par.type = CorrelatorType::cosh; - par.nState = strTo(sm[0].str()); + par.nState = strTo(sm[1].str()); } else if (regex_match(s, sm, regex("sinh([0-9]+)"))) { par.type = CorrelatorType::sinh; - par.nState = strTo(sm[0].str()); + par.nState = strTo(sm[1].str()); } else if (s == "linear") { @@ -189,7 +189,7 @@ DVec CorrelatorModels::parameterGuess(const DMatSample &corr, const ModelPar par) { DVec init; - Index nt = corr.size(); + Index nt = corr[central].size(); switch (par.type) { @@ -234,6 +234,11 @@ CorrelatorFitter::CorrelatorFitter(const std::vector &corr) setCorrelators(corr); } +XYSampleData & CorrelatorFitter::data(void) +{ + return *data_; +} + void CorrelatorFitter::setCorrelator(const DMatSample &corr) { std::vector vec; diff --git a/lib/Physics/CorrelatorFitter.hpp b/lib/Physics/CorrelatorFitter.hpp index fe4071e..018cb87 100644 --- a/lib/Physics/CorrelatorFitter.hpp +++ b/lib/Physics/CorrelatorFitter.hpp @@ -51,6 +51,8 @@ class CorrelatorFitter public: CorrelatorFitter(const DMatSample &corr); CorrelatorFitter(const std::vector &corr); + virtual ~CorrelatorFitter(void) = default; + XYSampleData & data(void); void setCorrelator(const DMatSample &corr); void setCorrelators(const std::vector &corr); const DMatSample & getCorrelator(const Index i = 0) const; diff --git a/lib/Physics/EffectiveMass.cpp b/lib/Physics/EffectiveMass.cpp index 672643f..aa21126 100644 --- a/lib/Physics/EffectiveMass.cpp +++ b/lib/Physics/EffectiveMass.cpp @@ -25,8 +25,9 @@ using namespace Latan; // constructors //////////////////////////////////////////////////////////////// EffectiveMass::EffectiveMass(const CorrelatorType type) -: type_(type) -{} +{ + setType(type); +} // access ////////////////////////////////////////////////////////////////////// CorrelatorType EffectiveMass::getType(void) const @@ -39,6 +40,31 @@ void EffectiveMass::setType(const CorrelatorType type) type_ = type; } +DVec EffectiveMass::getTime(const Index nt) const +{ + DVec tvec; + + switch (type_) + { + case CorrelatorType::undefined: + LATAN_ERROR(Argument, "correlator type is undefined"); + break; + case CorrelatorType::exp: + case CorrelatorType::linear: + tvec = DVec::LinSpaced(nt - 1, 0, nt - 2); + break; + case CorrelatorType::cosh: + case CorrelatorType::sinh: + tvec = DVec::LinSpaced(nt - 2, 1, nt - 2); + break; + case CorrelatorType::cst: + tvec = DVec::LinSpaced(nt, 0, nt - 1); + break; + } + + return tvec; +} + // compute effective mass ////////////////////////////////////////////////////// DVec EffectiveMass::operator()(const DVec &corr) const { diff --git a/lib/Physics/EffectiveMass.hpp b/lib/Physics/EffectiveMass.hpp index 500dab7..ed42c1e 100644 --- a/lib/Physics/EffectiveMass.hpp +++ b/lib/Physics/EffectiveMass.hpp @@ -37,6 +37,7 @@ public: // access CorrelatorType getType(void) const; void setType(const CorrelatorType type); + DVec getTime(const Index nt) const; // compute effective mass DVec operator()(const DVec &corr) const; DMatSample operator()(const DMatSample &corr) const; From 1775f4992bcdc9382f86d529de1f9c3529812d0f Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 28 Jan 2020 17:35:07 +0000 Subject: [PATCH 4/5] rewrite of the 2pt fitter using the new physics classes --- physics/2pt-fit.cpp | 321 ++++++++++++-------------------------------- 1 file changed, 84 insertions(+), 237 deletions(-) diff --git a/physics/2pt-fit.cpp b/physics/2pt-fit.cpp index 98ebdb8..4e86ee1 100644 --- a/physics/2pt-fit.cpp +++ b/physics/2pt-fit.cpp @@ -1,11 +1,13 @@ +#include #include +#include #include #include -#include -#include #include #include -#include +#include +#include +#include #include using namespace std; @@ -17,17 +19,6 @@ struct TwoPtFit Index tMin, tMax; }; -void setFitRange(XYSampleData &data, const Index ti, const Index tf, - const Index thinning, const Index nt) -{ - for (Index t = 0; t < nt; ++t) - { - data.fitPoint((t >= ti) and (t <= tf) - and ((t - ti) % thinning == 0), t); - } -} - - int main(int argc, char *argv[]) { // parse arguments ///////////////////////////////////////////////////////// @@ -47,7 +38,7 @@ int main(int argc, char *argv[]) opt.addOption("s", "shift" , OptParser::OptType::value , true, "time variable shift", "0"); opt.addOption("m", "model" , OptParser::OptType::value , true, - "fit model (exp|exp2|exp3|sinh|cosh|cosh2|cosh3|explin|const|)", "cosh"); + "fit model (exp|sinh|cosh|linear|cst|)", "exp1"); opt.addOption("" , "nPar" , OptParser::OptType::value , true, "number of model parameters for custom models " "(-1 if irrelevant)", "-1"); @@ -138,91 +129,15 @@ int main(int argc, char *argv[]) } } - // make models ///////////////////////////////////////////////////////////// - DoubleModel mod; - bool sinhModel = false, coshModel = false, linearModel = false, constModel = false; + // make model ////////////////////////////////////////////////////////////// + CorrelatorFitter fitter(corr); + DoubleModel mod; + auto modelPar = CorrelatorModels::parseModel(model); - if ((model == "exp") or (model == "exp1")) + if (modelPar.type != CorrelatorType::undefined) { - nPar = 2; - mod.setFunction([](const double *x, const double *p) - { - return p[1]*exp(-p[0]*x[0]); - }, 1, nPar); - } - else if (model == "exp2") - { - nPar = 4; - mod.setFunction([](const double *x, const double *p) - { - return p[1]*exp(-p[0]*x[0]) + p[3]*exp(-p[2]*x[0]); - }, 1, nPar); - } - else if (model == "exp3") - { - nPar = 6; - mod.setFunction([](const double *x, const double *p) - { - return p[1]*exp(-p[0]*x[0]) + p[3]*exp(-p[2]*x[0]) - + p[5]*exp(-p[4]*x[0]); - }, 1, nPar); - } - else if (model == "sinh") - { - sinhModel = true; - nPar = 2; - mod.setFunction([nt](const double *x, const double *p) - { - return p[1]*(exp(-p[0]*x[0])-exp(-p[0]*(nt-x[0]))); - }, 1, nPar); - } - else if ((model == "cosh") or (model == "cosh1")) - { - coshModel = true; - nPar = 2; - mod.setFunction([nt](const double *x, const double *p) - { - return p[1]*(exp(-p[0]*x[0])+exp(-p[0]*(nt-x[0]))); - }, 1, nPar); - } - else if (model == "cosh2") - { - coshModel = true; - nPar = 4; - mod.setFunction([nt](const double *x, const double *p) - { - return p[1]*(exp(-p[0]*x[0])+exp(-p[0]*(nt-x[0]))) - + p[3]*(exp(-p[2]*x[0])+exp(-p[2]*(nt-x[0]))); - }, 1, nPar); - } - else if (model == "cosh3") - { - coshModel = true; - nPar = 6; - mod.setFunction([nt](const double *x, const double *p) - { - return p[1]*(exp(-p[0]*x[0])+exp(-p[0]*(nt-x[0]))) - + p[3]*(exp(-p[2]*x[0])+exp(-p[2]*(nt-x[0]))) - + p[5]*(exp(-p[2]*x[0])+exp(-p[4]*(nt-x[0]))); - }, 1, nPar); - } - else if (model == "explin") - { - linearModel = true; - nPar = 2; - mod.setFunction([](const double *x, const double *p) - { - return p[1] - p[0]*x[0]; - }, 1, nPar); - } - else if (model == "const") - { - constModel = true; - nPar = 1; - mod.setFunction([](const double *x __dumb, const double *p) - { - return p[0]; - }, 1, nPar); + mod = CorrelatorModels::makeModel(modelPar, nt); + nPar = mod.getNPar(); } else { @@ -240,81 +155,44 @@ int main(int argc, char *argv[]) } // fit ///////////////////////////////////////////////////////////////////// - DMatSample tvec(nSample); - XYSampleData data(nSample); DVec init(nPar); NloptMinimizer globMin(NloptMinimizer::Algorithm::GN_CRS2_LM); MinuitMinimizer locMin; vector unCorrMin{&globMin, &locMin}; - FOR_STAT_ARRAY(tvec, s) - { - tvec[s] = DVec::LinSpaced(nt, 0, nt - 1); - } - data.addXDim(nt, "t/a", true); - data.addYDim("C(t)"); - data.setUnidimData(tvec, corr); - // set parameter name ****************************************************** - if(constModel) - { - mod.parName().setName(0, "const"); - } - else - { - for (Index p = 0; p < nPar; p += 2) - { - mod.parName().setName(p, "E_" + strFrom(p/2)); - mod.parName().setName(p + 1, "Z_" + strFrom(p/2)); - } - } - // set initial values ****************************************************** - if (linearModel) - { - init(0) = data.y(nt/4, 0)[central] - data.y(nt/4 + 1, 0)[central]; - init(1) = data.y(nt/4, 0)[central] + nt/4*init(0); - } - else if(constModel) - { - init(0) = data.y(nt/4, 0)[central]; + // set fitter ************************************************************** + fitter.setModel(mod); + fitter.data().setSvdTolerance(svdTol); + fitter.setThinning(thinning); + // set initial values ****************************************************** + if (modelPar.type != CorrelatorType::undefined) + { + init = CorrelatorModels::parameterGuess(corr, modelPar); } else { - init(0) = log(data.y(nt/4, 0)[central]/data.y(nt/4 + 1, 0)[central]); - init(1) = data.y(nt/4, 0)[central]/(exp(-init(0)*nt/4)); - } - for (Index p = 2; p < nPar; p += 2) - { - init(p) = 2*init(p - 2); - init(p + 1) = init(p - 1)/2.; - + init.fill(0.1); } + // set limits for minimisers *********************************************** for (Index p = 0; p < nPar; p += 2) { - if (linearModel) - { - globMin.setLowLimit(p, -10.*fabs(init(p))); - globMin.setHighLimit(p, 10.*fabs(init(p))); - } - else if(constModel) - { - globMin.setLowLimit(p, -10*fabs(init(0))); - locMin.setLowLimit(p, -10*fabs(init(0))); - globMin.setHighLimit(p, 10*fabs(init(0))); - locMin.setHighLimit(p, 10*fabs(init(0))); - } - else + if ((modelPar.type == CorrelatorType::exp) or + (modelPar.type == CorrelatorType::cosh) or + (modelPar.type == CorrelatorType::sinh)) { globMin.setLowLimit(p, 0.); + locMin.setLowLimit(p, 0.); globMin.setHighLimit(p, 10.*init(p)); - } - if(!constModel) - { globMin.setLowLimit(p + 1, -10.*fabs(init(p + 1))); globMin.setHighLimit(p + 1, 10.*fabs(init(p + 1))); } - + else + { + globMin.setLowLimit(p, -10*fabs(init(p))); + globMin.setHighLimit(p, 10*fabs(init(p))); + } } globMin.setPrecision(0.001); globMin.setMaxIteration(100000); @@ -322,28 +200,28 @@ int main(int argc, char *argv[]) locMin.setMaxIteration(1000000); locMin.setVerbosity(verbosity); - // fit ///////////////////////////////////////////////////////////////////// + // standard fit //////////////////////////////////////////////////////////// if (!doScan) { + // fit ***************************************************************** SampleFitResult fit; - setFitRange(data, ti, tf, thinning, nt); + fitter.setFitRange(ti, tf); if (doCorr) { cout << "-- uncorrelated fit..." << endl; } cout << "using model '" << model << "'" << endl; - data.setSvdTolerance(svdTol); - data.assumeYYCorrelated(false, 0, 0); - fit = data.fit(unCorrMin, init, mod); + fitter.setCorrelation(false); + fit = fitter.fit(unCorrMin, init); fit.print(); if (doCorr) { cout << "-- correlated fit..." << endl; cout << "using model '" << model << "'" << endl; init = fit[central]; - data.assumeYYCorrelated(true, 0, 0); - fit = data.fit(locMin, init, mod); + fitter.setCorrelation(true); + fit = fitter.fit(locMin, init); fit.print(); } if (!outFileName.empty()) @@ -353,84 +231,50 @@ int main(int argc, char *argv[]) // plots *************************************************************** if (doPlot) { - if (!constModel) + DMatSample tvec(nSample); + + tvec.fill(DVec::LinSpaced(nt, 0, nt - 1)); + if (modelPar.type != CorrelatorType::cst) { Plot p; p << PlotRange(Axis::x, 0, nt - 1); - if (!linearModel and !constModel) + if ((modelPar.type == CorrelatorType::exp) or + (modelPar.type == CorrelatorType::cosh) or + (modelPar.type == CorrelatorType::sinh)) { p << LogScale(Axis::y); } p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1); p << Color("rgb 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1); - p << Color("rgb 'red'") << PlotData(data.getData()); + p << Color("rgb 'red'") << PlotData(fitter.data().getData()); + p << Label("t/a", Axis::x) << Caption("Correlator"); p.display(); if(savePlot != "") { p.save(savePlot + "_corr"); } } + if (modelPar.type != CorrelatorType::undefined) { - Plot p; - DMatSample effMass(nSample); - DVec effMassT, fitErr; - Index maxT = (coshModel) ? (nt - 2) : (nt - 1); - double e0, e0Err; + Plot p; + EffectiveMass effMass(modelPar.type); + DMatSample em; + DVec fitErr, emtvec; + double e0, e0Err; - effMass.resizeMat(maxT, 1); - effMassT.setLinSpaced(maxT, 0, maxT-1); + emtvec = effMass.getTime(nt); + em = effMass(corr); fitErr = fit.variance().cwiseSqrt(); e0 = fit[central](0); e0Err = fitErr(0); - if (coshModel or sinhModel) - { - FOR_STAT_ARRAY(effMass, s) - { - for (Index t = 1; t < nt - 1; ++t) - { - effMass[s](t - 1) = acosh((corr[s](t-1) + corr[s](t+1)) - /(2.*corr[s](t))); - } - } - } - else if (linearModel) - { - FOR_STAT_ARRAY(effMass, s) - { - for (Index t = 0; t < nt - 1; ++t) - { - effMass[s](t) = corr[s](t) - corr[s](t+1); - } - } - } - else if (constModel) - { - FOR_STAT_ARRAY(effMass, s) - { - for (Index t = 0; t < nt - 1; ++t) - { - effMass[s](t) = corr[s](t); - } - } - } - else - { - FOR_STAT_ARRAY(effMass, s) - { - for (Index t = 1; t < nt; ++t) - { - effMass[s](t - 1) = log(corr[s](t-1)/corr[s](t)); - } - } - } p.reset(); - p << PlotRange(Axis::x, 0, maxT); - p << PlotRange(Axis::y, e0 - 20.*e0Err, e0 + 20.*e0Err); - p << Color("rgb 'blue'") << PlotBand(0, maxT, e0 - e0Err, e0 + e0Err); + p << PlotRange(Axis::x, 0, nt - 1); + p << PlotRange(Axis::y, e0 - 30.*e0Err, e0 + 30.*e0Err); + p << Color("rgb 'blue'") << PlotBand(0, nt - 1, e0 - e0Err, e0 + e0Err); p << Color("rgb 'blue'") << PlotHLine(e0); - p << Color("rgb 'red'") << PlotData(effMassT, effMass); - p << Caption("Effective Mass"); + p << Color("rgb 'red'") << PlotData(emtvec, em); + p << Label("t/a", Axis::x) << Caption("Effective Mass"); p.display(); if(savePlot != "") { @@ -440,16 +284,19 @@ int main(int argc, char *argv[]) if (doHeatmap) { Plot p; - Index n = data.getFitVarMat().rows(); - DMat id = DMat::Identity(n, n); + Index n = fitter.data().getFitVarMat().rows(); + DMat id = DMat::Identity(n, n), + var = fitter.data().getFitVarMat(); - p << PlotMatrix(Math::varToCorr(data.getFitVarMat())); + p << PlotMatrix(Math::varToCorr(var)); p << Caption("correlation matrix"); p.display(); if (svdTol > 0.) { + DMat proj = id - var*fitter.data().getFitVarMatPInv(); + p.reset(); - p << PlotMatrix(id - data.getFitVarMat()*data.getFitVarMatPInv()); + p << PlotMatrix(proj); p << Caption("singular space projector"); p.display(); } @@ -460,8 +307,9 @@ int main(int argc, char *argv[]) // scan fits /////////////////////////////////////////////////////////////// else { + // fits **************************************************************** Index nFit = 0, f = 0, ti0 = ti + (tf - ti)/4, tf0 = tf - (tf - ti)/4, - matSize = tf - ti - nPar + 1; + matSize = tf - ti + 1; DMat err, pVal(matSize, matSize), relErr(matSize, matSize), ccdf(matSize, matSize), val(matSize, matSize); map fit; @@ -474,14 +322,13 @@ int main(int argc, char *argv[]) << endl; thinning = 1; } - setFitRange(data, ti0, tf0, thinning, nt); - data.setSvdTolerance(svdTol); - data.assumeYYCorrelated(false, 0, 0); - tmpFit = data.fit(unCorrMin, init, mod); + fitter.setFitRange(ti0, tf0); + fitter.setCorrelation(false); + tmpFit = fitter.fit(unCorrMin, init); tmpFit.print(); cout << "-- scanning all possible fit ranges..." << endl; init = tmpFit[central]; - data.assumeYYCorrelated(doCorr, 0, 0); + fitter.setCorrelation(doCorr); pVal.fill(Math::nan); relErr.fill(Math::nan); val.fill(Math::nan); @@ -496,8 +343,8 @@ int main(int argc, char *argv[]) { Index i = ta - ti, j = tb - ti; - setFitRange(data, ta, tb, thinning, nt); - tmpFit = data.fit(locMin, init, mod); + fitter.setFitRange(ta, tb); + tmpFit = fitter.fit(locMin, init); err = tmpFit.variance().cwiseSqrt(); pVal(i, j) = tmpFit.getPValue(); ccdf(i, j) = tmpFit.getCcdf(); @@ -531,8 +378,8 @@ int main(int argc, char *argv[]) p << PlotMatrix(pVal); p << Caption("p-value matrix"); - p << Label("tMin - " + strFrom(ti), Axis::x); - p << Label("tMax - " + strFrom(ti), Axis::y); + p << Label("tMax - " + strFrom(ti), Axis::x); + p << Label("tMin - " + strFrom(ti), Axis::y); p.display(); if(savePlot != "") { @@ -541,8 +388,8 @@ int main(int argc, char *argv[]) p.reset(); p << PlotMatrix(relErr); p << Caption("Relative error matrix"); - p << Label("tMin - " + strFrom(ti), Axis::x); - p << Label("tMax - " + strFrom(ti), Axis::y); + p << Label("tMax - " + strFrom(ti), Axis::x); + p << Label("tMin - " + strFrom(ti), Axis::y); p.display(); if(savePlot != "") { @@ -551,8 +398,8 @@ int main(int argc, char *argv[]) p.reset(); p << PlotMatrix(val); p << Caption("Fit result matrix"); - p << Label("tMin - " + strFrom(ti), Axis::x); - p << Label("tMax - " + strFrom(ti), Axis::y); + p << Label("tMax - " + strFrom(ti), Axis::x); + p << Label("tMin - " + strFrom(ti), Axis::y); p.display(); if(savePlot != "") { @@ -561,8 +408,8 @@ int main(int argc, char *argv[]) p.reset(); p << PlotMatrix(ccdf); p << Caption("chi^2 CCDF matrix"); - p << Label("tMin - " + strFrom(ti), Axis::x); - p << Label("tMax - " + strFrom(ti), Axis::y); + p << Label("tMax - " + strFrom(ti), Axis::x); + p << Label("tMin - " + strFrom(ti), Axis::y); p.display(); if(savePlot != "") { From 1d6a66263d45c84f18bd5d9fe108c525cd31561c Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 28 Jan 2020 17:55:47 +0000 Subject: [PATCH 5/5] code cleaning --- lib/Physics/CorrelatorFitter.cpp | 12 +++++++++++- lib/Physics/CorrelatorFitter.hpp | 13 ++++++++++++- lib/Physics/EffectiveMass.cpp | 3 +++ 3 files changed, 26 insertions(+), 2 deletions(-) diff --git a/lib/Physics/CorrelatorFitter.cpp b/lib/Physics/CorrelatorFitter.cpp index d48699f..c997272 100644 --- a/lib/Physics/CorrelatorFitter.cpp +++ b/lib/Physics/CorrelatorFitter.cpp @@ -23,6 +23,9 @@ using namespace std; using namespace Latan; +/****************************************************************************** + * Correlator models * + ******************************************************************************/ DoubleModel CorrelatorModels::makeExpModel(const Index nState) { DoubleModel mod; @@ -224,6 +227,10 @@ DVec CorrelatorModels::parameterGuess(const DMatSample &corr, return init; } +/****************************************************************************** + * CorrelatorFitter implementation * + ******************************************************************************/ +// constructors //////////////////////////////////////////////////////////////// CorrelatorFitter::CorrelatorFitter(const DMatSample &corr) { setCorrelator(corr); @@ -234,6 +241,7 @@ CorrelatorFitter::CorrelatorFitter(const std::vector &corr) setCorrelators(corr); } +// access ////////////////////////////////////////////////////////////////////// XYSampleData & CorrelatorFitter::data(void) { return *data_; @@ -294,7 +302,7 @@ void CorrelatorFitter::setCorrelation(const bool isCorrelated, const Index i, data_->assumeYYCorrelated(isCorrelated, i, j); } -DMat CorrelatorFitter::getVarianceMatrix(void) +DMat CorrelatorFitter::getVarianceMatrix(void) const { return data_->getFitVarMat(); } @@ -305,6 +313,7 @@ void CorrelatorFitter::setThinning(const Index thinning, const Index i) refreshRanges(); } +// fit functions /////////////////////////////////////////////////////////////// SampleFitResult CorrelatorFitter::fit(Minimizer &minimizer, const DVec &init) { vector vecPt = {&minimizer}; @@ -325,6 +334,7 @@ SampleFitResult CorrelatorFitter::fit(vector &minimizer, return data_->fit(minimizer, init, vecPt); } +// internal function to refresh fit ranges ///////////////////////////////////// void CorrelatorFitter::refreshRanges(void) { for (unsigned int i = 0; i < range_.size(); ++i) diff --git a/lib/Physics/CorrelatorFitter.hpp b/lib/Physics/CorrelatorFitter.hpp index 018cb87..9acadb7 100644 --- a/lib/Physics/CorrelatorFitter.hpp +++ b/lib/Physics/CorrelatorFitter.hpp @@ -26,6 +26,9 @@ BEGIN_LATAN_NAMESPACE +/****************************************************************************** + * Correlator types & models * + ******************************************************************************/ enum class CorrelatorType {undefined, exp, cosh, sinh, linear, cst}; namespace CorrelatorModels @@ -46,12 +49,18 @@ namespace CorrelatorModels DVec parameterGuess(const DMatSample &corr, const ModelPar par); }; +/****************************************************************************** + * Correlator fit utility class * + ******************************************************************************/ class CorrelatorFitter { public: + // constructors CorrelatorFitter(const DMatSample &corr); CorrelatorFitter(const std::vector &corr); + // destructor virtual ~CorrelatorFitter(void) = default; + // access XYSampleData & data(void); void setCorrelator(const DMatSample &corr); void setCorrelators(const std::vector &corr); @@ -62,11 +71,13 @@ public: void setFitRange(const Index tMin, const Index tMax, const Index i = 0); void setCorrelation(const bool isCorrelated, const Index i = 0, const Index j = 0); - DMat getVarianceMatrix(void); + DMat getVarianceMatrix(void) const; void setThinning(const Index thinning, const Index i = 0); + // fit functions SampleFitResult fit(Minimizer &minimizer, const DVec &init); SampleFitResult fit(std::vector &minimizer, const DVec &init); private: + // internal function to refresh fit ranges void refreshRanges(void); private: Index nt_{0}; diff --git a/lib/Physics/EffectiveMass.cpp b/lib/Physics/EffectiveMass.cpp index aa21126..7899843 100644 --- a/lib/Physics/EffectiveMass.cpp +++ b/lib/Physics/EffectiveMass.cpp @@ -23,6 +23,9 @@ using namespace std; using namespace Latan; +/****************************************************************************** + * EffectiveMass implementation * + ******************************************************************************/ // constructors //////////////////////////////////////////////////////////////// EffectiveMass::EffectiveMass(const CorrelatorType type) {