From f826f30e82a7eef4bd3be44b22ec28d7db626d75 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 23 Jan 2020 19:11:01 +0000 Subject: [PATCH] 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