Merge branch 'feature/correlator-fitter' into develop

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);
+    FOR_MAT(low, i, j)
-    highFileName = dumpToTmpFile(high);
+    {
-    pushTmpFile(lowFileName);
+        contour(i, j) = low(i, j);
-    pushTmpFile(highFileName);
+    }
-    setCommand("'< (cat " + lowFileName + "; tac " + highFileName +
+    FOR_MAT(high, i, j)
-               "; head -n1 " + lowFileName + ")' u 1:2 w filledcurves closed" +
+    {
+        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");
 }
 

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         \

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 <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Physics/CorrelatorFitter.hpp>
+#include <LatAnalyze/includes.hpp>
+
+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<Index>(sm[1].str());
+    }
+    else if (regex_match(s, sm, regex("cosh([0-9]+)")))
+    {
+        par.type   = CorrelatorType::cosh;
+        par.nState = strTo<Index>(sm[1].str());
+    }
+    else if (regex_match(s, sm, regex("sinh([0-9]+)")))
+    {
+        par.type   = CorrelatorType::sinh;
+        par.nState = strTo<Index>(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<DMatSample> &corr)
+{
+    setCorrelators(corr);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+XYSampleData & CorrelatorFitter::data(void)
+{
+    return *data_;
+}
+
+void CorrelatorFitter::setCorrelator(const DMatSample &corr)
+{
+    std::vector<DMatSample> vec;
+
+    vec.push_back(corr);
+    setCorrelators(vec);
+}   
+
+void CorrelatorFitter::setCorrelators(const std::vector<DMatSample> &corr)
+{
+    Index                           nSample = corr[0].size();
+    DMatSample                      tVec(nSample);
+    std::vector<const DMatSample *> 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<Minimizer *> vecPt = {&minimizer};
+
+    return fit(vecPt, init);
+}
+
+SampleFitResult CorrelatorFitter::fit(vector<Minimizer *> &minimizer,
+                                      const DVec &init)
+{
+    vector<const DoubleModel *> 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);
+    }
+}

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 <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_CorrelatorFitter_hpp_
+#define Latan_CorrelatorFitter_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Functional/Model.hpp>
+#include <LatAnalyze/Statistics/XYSampleData.hpp>
+
+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<DMatSample> &corr);
+    // destructor
+    virtual ~CorrelatorFitter(void) = default;
+    // access
+    XYSampleData & data(void);
+    void setCorrelator(const DMatSample &corr);
+    void setCorrelators(const std::vector<DMatSample> &corr);
+    const DMatSample & getCorrelator(const Index i = 0) const;
+    const std::vector<DMatSample> & 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 *> &minimizer, const DVec &init);
+private:
+    // internal function to refresh fit ranges
+    void refreshRanges(void);
+private:
+    Index                                nt_{0};
+    std::unique_ptr<XYSampleData>        data_;
+    std::vector<DoubleModel>             model_;
+    std::vector<std::pair<Index, Index>> range_;
+    std::vector<Index>                   thinning_;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_CorrelatorFitter_hpp_

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 <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Physics/EffectiveMass.hpp>
+#include <LatAnalyze/includes.hpp>
+
+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;
+}

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 <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_EffectiveMass_hpp_
+#define Latan_EffectiveMass_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
+#include <LatAnalyze/Physics/CorrelatorFitter.hpp>
+
+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_

physics/2pt-fit.cpp

@@ -1,11 +1,13 @@
+#include <LatAnalyze/Core/Math.hpp>
 #include <LatAnalyze/Core/OptParser.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
 #include <LatAnalyze/Functional/CompiledModel.hpp>
 #include <LatAnalyze/Io/Io.hpp>
-#include <LatAnalyze/Statistics/MatSample.hpp>
-#include <LatAnalyze/Core/Math.hpp>
 #include <LatAnalyze/Numerical/MinuitMinimizer.hpp>
 #include <LatAnalyze/Numerical/NloptMinimizer.hpp>
-#include <LatAnalyze/Core/Plot.hpp>
+#include <LatAnalyze/Physics/CorrelatorFitter.hpp>
+#include <LatAnalyze/Physics/EffectiveMass.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
 #include <LatAnalyze/Statistics/XYSampleData.hpp>
 
 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|<interpreter code>)", "cosh");
+                  "fit model (exp<n>|sinh<n>|cosh<n>|linear|cst|<interpreter code>)", "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 /////////////////////////////////////////////////////////////
+    // make model //////////////////////////////////////////////////////////////
-    DoubleModel mod;
+    CorrelatorFitter fitter(corr);
-    bool        sinhModel = false, coshModel = false, linearModel = false, constModel = false;
+    DoubleModel      mod;
+    auto             modelPar = CorrelatorModels::parseModel(model);
     
-    if ((model == "exp") or (model == "exp1"))
+    if (modelPar.type != CorrelatorType::undefined)
     {
-        nPar = 2;
+        mod  = CorrelatorModels::makeModel(modelPar, nt);
-        mod.setFunction([](const double *x, const double *p)
+        nPar = mod.getNPar();
-                        {
-                            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);
     }
     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<Minimizer *> unCorrMin{&globMin, &locMin};
 
-    FOR_STAT_ARRAY(tvec, s)
+    // set fitter **************************************************************
-    {
+    fitter.setModel(mod);
-        tvec[s] = DVec::LinSpaced(nt, 0, nt - 1);
+    fitter.data().setSvdTolerance(svdTol);
-    }
+    fitter.setThinning(thinning);
-    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 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.fill(0.1);
-        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.;
-        
     }
+
     // set limits for minimisers ***********************************************
     for (Index p = 0; p < nPar; p += 2)
     {
-        if (linearModel)
+        if ((modelPar.type == CorrelatorType::exp) or
-        {
+            (modelPar.type == CorrelatorType::cosh) or
-            globMin.setLowLimit(p, -10.*fabs(init(p)));
+            (modelPar.type == CorrelatorType::sinh))
-            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
         {
             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);
+        fitter.setCorrelation(false);
-        data.assumeYYCorrelated(false, 0, 0);
+        fit = fitter.fit(unCorrMin, init);
-        fit = data.fit(unCorrMin, init, mod);
         fit.print();
         if (doCorr)
         {
             cout << "-- correlated fit..." << endl;
             cout << "using model '" << model << "'" << endl;
             init = fit[central];
-            data.assumeYYCorrelated(true, 0, 0);
+            fitter.setCorrelation(true);
-            fit = data.fit(locMin, init, mod);
+            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;
+                Plot          p;
-                DMatSample effMass(nSample);
+                EffectiveMass effMass(modelPar.type);
-                DVec       effMassT, fitErr;
+                DMatSample    em;
-                Index      maxT = (coshModel) ? (nt - 2) : (nt - 1);
+                DVec          fitErr, emtvec;
-                double     e0, e0Err;
+                double        e0, e0Err;
 
-                effMass.resizeMat(maxT, 1);
+                emtvec = effMass.getTime(nt);
-                effMassT.setLinSpaced(maxT, 0, maxT-1);
+                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::x, 0, nt - 1);
-                p << PlotRange(Axis::y, e0 - 20.*e0Err, e0 + 20.*e0Err);
+                p << PlotRange(Axis::y, e0 - 30.*e0Err, e0 + 30.*e0Err);
-                p << Color("rgb 'blue'") << PlotBand(0, maxT, e0 - e0Err, e0 + 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 << Color("rgb 'red'") << PlotData(emtvec, em);
-                p << Caption("Effective Mass");
+                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();
+                Index n  = fitter.data().getFitVarMat().rows();
-                DMat  id = DMat::Identity(n, n);
+                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<double, TwoPtFit> fit;
@@ -474,14 +322,13 @@ int main(int argc, char *argv[])
                  << endl;
             thinning = 1;
         }
-        setFitRange(data, ti0, tf0, thinning, nt);
+        fitter.setFitRange(ti0, tf0);
-        data.setSvdTolerance(svdTol);
+        fitter.setCorrelation(false);
-        data.assumeYYCorrelated(false, 0, 0);
+        tmpFit = fitter.fit(unCorrMin, init);
-        tmpFit = data.fit(unCorrMin, init, mod);
         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);
+            fitter.setFitRange(ta, tb);
-            tmpFit              = data.fit(locMin, init, mod);
+            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::x);
-            p << Label("tMax - " + strFrom(ti), Axis::y);
+            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::x);
-            p << Label("tMax - " + strFrom(ti), Axis::y);
+            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::x);
-            p << Label("tMax - " + strFrom(ti), Axis::y);
+            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::x);
-            p << Label("tMax - " + strFrom(ti), Axis::y);
+            p << Label("tMin - " + strFrom(ti), Axis::y);
             p.display();
             if(savePlot != "")
             {