first stab at correlator fit utility classes

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,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 <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Physics/CorrelatorFitter.hpp>
+#include <LatAnalyze/includes.hpp>
+
+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<Index>(sm[0].str());
+    }
+    else if (regex_match(s, sm, regex("cosh([0-9]+)")))
+    {
+        par.type   = CorrelatorType::cosh;
+        par.nState = strTo<Index>(sm[0].str());
+    }
+    else if (regex_match(s, sm, regex("sinh([0-9]+)")))
+    {
+        par.type   = CorrelatorType::sinh;
+        par.nState = strTo<Index>(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<DMatSample> &corr)
+{
+    setCorrelators(corr);
+}
+
+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)
+{
+    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<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);
+}
+
+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,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 <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
+
+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<DMatSample> &corr);
+    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);
+    void setThinning(const Index thinning, const Index i = 0);
+    SampleFitResult fit(Minimizer &minimizer, const DVec &init);
+    SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init);
+private:
+    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,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 <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Physics/EffectiveMass.hpp>
+#include <LatAnalyze/includes.hpp>
+
+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;
+}

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 <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);
+    // 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_