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");
}
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..c997272
--- /dev/null
+++ b/lib/Physics/CorrelatorFitter.cpp
@@ -0,0 +1,346 @@
+/*
+ * 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;
+
+/******************************************************************************
+ * Correlator models *
+ ******************************************************************************/
+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[1].str());
+ }
+ else if (regex_match(s, sm, regex("cosh([0-9]+)")))
+ {
+ par.type = CorrelatorType::cosh;
+ par.nState = strTo(sm[1].str());
+ }
+ else if (regex_match(s, sm, regex("sinh([0-9]+)")))
+ {
+ par.type = CorrelatorType::sinh;
+ par.nState = strTo(sm[1].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[central].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 implementation *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+CorrelatorFitter::CorrelatorFitter(const DMatSample &corr)
+{
+ setCorrelator(corr);
+}
+
+CorrelatorFitter::CorrelatorFitter(const std::vector &corr)
+{
+ setCorrelators(corr);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+XYSampleData & CorrelatorFitter::data(void)
+{
+ return *data_;
+}
+
+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) const
+{
+ return data_->getFitVarMat();
+}
+
+void CorrelatorFitter::setThinning(const Index thinning, const Index i)
+{
+ thinning_[i] = thinning;
+ refreshRanges();
+}
+
+// fit functions ///////////////////////////////////////////////////////////////
+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);
+}
+
+// internal function to refresh fit ranges /////////////////////////////////////
+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..9acadb7
--- /dev/null
+++ b/lib/Physics/CorrelatorFitter.hpp
@@ -0,0 +1,92 @@
+/*
+ * 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
+
+/******************************************************************************
+ * Correlator types & models *
+ ******************************************************************************/
+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);
+};
+
+/******************************************************************************
+ * 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);
+ 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) 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};
+ 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..7899843
--- /dev/null
+++ b/lib/Physics/EffectiveMass.cpp
@@ -0,0 +1,132 @@
+/*
+ * 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;
+
+/******************************************************************************
+ * EffectiveMass implementation *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+EffectiveMass::EffectiveMass(const CorrelatorType type)
+{
+ setType(type);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+CorrelatorType EffectiveMass::getType(void) const
+{
+ return type_;
+}
+
+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
+{
+ 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..ed42c1e
--- /dev/null
+++ b/lib/Physics/EffectiveMass.hpp
@@ -0,0 +1,50 @@
+/*
+ * 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);
+ DVec getTime(const Index nt) const;
+ // 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
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 != "")
{