Update Readme.md

Changed linear model's parameters

Readme.md

@@ -1,21 +1,6 @@
 # LatAnalyze 
 
-License: GNU General Public License v3
+[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![DOI](https://zenodo.org/badge/10201777.svg)](https://zenodo.org/badge/latestdoi/10201777) [![Build Ubuntu](https://github.com/aportelli/LatAnalyze/actions/workflows/build-ubuntu.yml/badge.svg)](https://github.com/aportelli/LatAnalyze/actions/workflows/build-ubuntu.yml) [![Build macOS](https://github.com/aportelli/LatAnalyze/actions/workflows/build-macos.yml/badge.svg)](https://github.com/aportelli/LatAnalyze/actions/workflows/build-macos.yml)
-
-<table>
-<tr>
-    <td>Last stable release</td>
-    <td><a href="https://travis-ci.org/aportelli/LatAnalyze">
-    <img src="https://travis-ci.org/aportelli/LatAnalyze.svg?branch=master"></a>
-    </td>
-</tr>
-<tr>
-    <td>Development branch</td>
-    <td><a href="https://travis-ci.org/aportelli/LatAnalyze">
-    <img src="https://travis-ci.org/aportelli/LatAnalyze.svg?branch=develop"></a>
-    </td>
-</tr>
-</table>
 
 ## Description
 LatAnalyze is a C++11 library for statistical data analysis based on bootstrap
@@ -71,6 +56,9 @@ make install
 where `<prefix>` should be replaced by the desired prefix for LatAnalyze installation, and `<n>` is the number of parallel build processes (typically twice your number of cores).
 
 ## History
+#### v3.5.1
+Various fixes and cleaning of outdated code.
+
 #### v3.5
 Additions:
 * 'Impulse' & line type plots
@@ -161,4 +149,4 @@ Fixes:
 #### v3.0
 Commit `7b4f2884a5e99bbfab4d4bd7623f609a55403c39`.  
 First 'stable' version of LatAnalyze in C++. The v2.0 refers to the [C version](https://github.com/aportelli/LatAnalyze-legacy) and v1.0 to an old undistributed version.  
-**This version compiles fine on OS X with clang but does have many portability issues to other platforms/compilers, v3.1 is the first real release.**
+**This version compiles fine on OS X with clang but does have many portability issues to other platforms/compilers, v3.1 is the first real release.**

bootstrap.sh

@@ -2,5 +2,5 @@
 
 rm -rf .buildutils
 mkdir -p .buildutils/m4
-./update_eigen.sh eigen-3.3.7.tar.bz2
+./update_eigen.sh eigen-3.3.8.tar.bz2
 autoreconf -fvi

ci-scripts/install-gsl.sh

@@ -9,8 +9,11 @@ fi
 PREFIX=$1
 
 set -ex
-INITDIR=`pwd`
+INITDIR=$(pwd -P)
-cd local/build
+mkdir -p ${PREFIX}
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}/local/build
 wget http://ftpmirror.gnu.org/gsl/${NAME}.tar.gz
 tar -xzvf ${NAME}.tar.gz
 mkdir -p ${NAME}/build

ci-scripts/install-hdf5.sh

@@ -9,8 +9,11 @@ fi
 PREFIX=$1
 
 set -ex
-INITDIR=`pwd`
+INITDIR=$(pwd -P)
-cd local/build
+mkdir -p ${PREFIX}
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}/local/build
 wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/${NAME}/src/${NAME}.tar.gz
 tar -xzvf ${NAME}.tar.gz
 mkdir ${NAME}/build

ci-scripts/install-latan.sh

@@ -5,14 +5,18 @@ if (( $# != 1 )); then
   exit 1
 fi
 PREFIX=$1
-OS=$2
 
 set -ex
+INITDIR=$(pwd -P)
+mkdir -p ${PREFIX}
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}
 ./install-deps.sh ${PREFIX}
 cd ..
 ./bootstrap.sh
 mkdir -p build
 cd build
-../configure --prefix=$PREFIX --with-minuit=$PREFIX --with-nlopt=$PREFIX --with-latcore=$PREFIX --with-hdf5=$PREFIX --with-gsl=$PREFIX CXXFLAGS="${CXXFLAGS} -O3 -march=native -mtune=native"
+../configure --prefix=${PREFIX} --with-minuit=${PREFIX} --with-nlopt=${PREFIX} --with-hdf5=${PREFIX} --with-gsl=${PREFIX} CXXFLAGS="${CXXFLAGS} -O3 -march=haswell -mtune=haswell"
 make -j4
 make install

ci-scripts/install-minuit.sh

@@ -9,8 +9,11 @@ fi
 PREFIX=$1
 
 set -ex
-INITDIR=`pwd`
+INITDIR=$(pwd -P)
-cd local/build
+mkdir -p ${PREFIX}
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}/local/build
 wget http://www.cern.ch/mathlibs/sw/5_34_14/Minuit2/${NAME}.tar.gz
 tar -xzvf ${NAME}.tar.gz
 mkdir -p ${NAME}/build

ci-scripts/install-nlopt.sh

@@ -9,14 +9,17 @@ fi
 PREFIX=$1
 
 set -ex
-INITDIR=`pwd`
+INITDIR=$(pwd -P)
-cd local/build
+mkdir -p ${PREFIX}
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}/local/build
 wget https://github.com/stevengj/nlopt/archive/v${NAME}.tar.gz
 tar -xzvf v${NAME}.tar.gz
 NAME=nlopt-${NAME}
 mkdir -p ${NAME}/build
 cd ${NAME}/build
-cmake -DCMAKE_INSTALL_PREFIX=${PREFIX} ..
+cmake -DCMAKE_INSTALL_PREFIX=${PREFIX} -DCMAKE_BUILD_WITH_INSTALL_NAME_DIR=TRUE -DCMAKE_INSTALL_NAME_DIR="${PREFIX}/lib" ..
 make -j4 
 make install
 cd ${INITDIR}/local

configure.ac

@@ -2,7 +2,7 @@
 
 # Initialization
 AC_PREREQ([2.63])
-AC_INIT([LatAnalyze],[3.5],[antonin.portelli@me.com],[LatAnalyze])
+AC_INIT([LatAnalyze],[3.5.1],[antonin.portelli@me.com],[LatAnalyze])
 AC_CONFIG_AUX_DIR([.buildutils])
 AC_CONFIG_SRCDIR([lib/Global.cpp])
 AC_CONFIG_SRCDIR([utils/sample_read.cpp])

lib/Core/Eigen.hpp

@@ -1,7 +1,7 @@
 /*
  * Eigen.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/Exceptions.cpp

@@ -1,7 +1,7 @@
 /*
  * Exceptions.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/Exceptions.hpp

@@ -1,7 +1,7 @@
 /*
  * Exceptions.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/Mat.cpp

@@ -1,7 +1,7 @@
 /*
  * Mat.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/Mat.hpp

@@ -1,7 +1,7 @@
 /*
  * Mat.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/Math.cpp

@@ -1,7 +1,7 @@
 /*
  * Math.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/Math.hpp

@@ -1,7 +1,7 @@
 /*
  * Math.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/MathInterpreter.cpp

@@ -1,7 +1,7 @@
 /*
  * MathInterpreter.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/MathInterpreter.hpp

@@ -1,7 +1,7 @@
 /*
  * MathInterpreter.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/MathLexer.lpp

@@ -1,7 +1,7 @@
 /*
  * MathLexer.lpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/MathParser.ypp

@@ -1,7 +1,7 @@
 /*
  * MathParser.ypp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/ParserState.hpp

@@ -1,7 +1,7 @@
 /*
  * ParserState.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Core/Plot.cpp

@@ -1,7 +1,7 @@
 /*
  * Plot.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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
@@ -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");
 }
 
@@ -500,7 +507,7 @@ Plot::Plot(void)
 // default options /////////////////////////////////////////////////////////////
 void Plot::initOptions(void)
 {
-    options_.terminal     = "qt";
+    options_.terminal     = "qt noenhanced font 'Arial,12'";
     options_.output       = "";
     options_.caption      = "";
     options_.title        = "";
@@ -580,57 +587,48 @@ Plot & Plot::operator<<(PlotModifier &&modifier)
 }
 
 // find gnuplot ////////////////////////////////////////////////////////////////
-void Plot::getProgramPath(void)
+#define SEARCH_DIR(dir) \
+sprintf(buf, "%s/%s", dir, gnuplotBin_.c_str());\
+if (access(buf, X_OK) == 0)\
+{\
+    return dir;\
+}
+
+std::string Plot::getProgramPath(void)
 {
     int         i, j, lg;
     char        *path;
     static char buf[MAX_PATH_LENGTH];
     
-    // trivial case: try in CWD
-    sprintf(buf,"./%s", gnuplotBin_.c_str()) ;
-    if (access(buf, X_OK) == 0)
-    {
-        sprintf(buf,".");
-        gnuplotPath_ = buf;
-    }
     // try out in all paths given in the PATH variable
-    else
+    buf[0] = 0;
+    path = getenv("PATH") ;
+    if (path)
     {
-        buf[0] = 0;
+        for (i=0;path[i];)
-        path = getenv("PATH") ;
-        if (path)
         {
-            for (i=0;path[i];)
+            for (j=i;(path[j]) and (path[j]!=':');j++);
+            lg = j - i;
+            strncpy(buf,path + i,(size_t)(lg));
+            if (lg == 0)
             {
-                for (j=i;(path[j]) and (path[j]!=':');j++);
+                buf[lg++] = '.';
-                lg = j - i;
-                strncpy(buf,path + i,(size_t)(lg));
-                if (lg == 0)
-                {
-                    buf[lg++] = '.';
-                }
-                buf[lg++] = '/';
-                strcpy(buf + lg, gnuplotBin_.c_str());
-                if (access(buf, X_OK) == 0)
-                {
-                    // found it!
-                    break ;
-                }
-                buf[0] = 0;
-                i = j;
-                if (path[i] == ':') i++ ;
             }
+            buf[lg++] = '/';
+            strcpy(buf + lg, gnuplotBin_.c_str());
+            if (access(buf, X_OK) == 0)
+            {
+                // found it!
+                break ;
+            }
+            buf[0] = 0;
+            i = j;
+            if (path[i] == ':') i++ ;
         }
-        else
+    }
-        {
+    // if the buffer is still empty, the command was not found
-            LATAN_ERROR(System, "PATH variable not set");
+    if (buf[0] != 0)
-        }
+    {
-        // if the buffer is still empty, the command was not found
-        if (buf[0] == 0)
-        {
-            LATAN_ERROR(System, "cannot find gnuplot in $PATH");
-        }
-        // otherwise truncate the command name to yield path only
         lg = (int)(strlen(buf) - 1);
         while (buf[lg]!='/')
         {
@@ -639,7 +637,19 @@ void Plot::getProgramPath(void)
         }
         buf[lg] = 0;
         gnuplotPath_ = buf;
+
+        return gnuplotPath_;
     }
+
+    // try in CWD, /usr/bin & /usr/local/bin
+    SEARCH_DIR(".");
+    SEARCH_DIR("/usr/bin");
+    SEARCH_DIR("/usr/local/bin");
+
+    // if this code is reached nothing was found
+    LATAN_ERROR(System, "cannot find gnuplot");
+
+    return "";
 }
 
 // plot parsing and output /////////////////////////////////////////////////////
@@ -653,10 +663,6 @@ void Plot::display(void)
         string        command;
         ostringstream scriptBuf;
 
-        if (!getenv("DISPLAY"))
-        {
-            LATAN_ERROR(System, "cannot find DISPLAY variable: is it set ?");
-        }
         getProgramPath();
         command     = gnuplotPath_ + "/" + gnuplotBin_ + " 2>/dev/null";
         gnuplotPipe = popen(command.c_str(), "w");
@@ -683,7 +689,7 @@ void Plot::display(void)
     }
 }
 
-void Plot::save(string dirName)
+void Plot::save(string dirName, bool savePdf)
 {
     vector<string> commandBack;
     string         path, terminalBack, outputBack, gpCommand, scriptName;
@@ -691,11 +697,6 @@ void Plot::save(string dirName)
     ofstream       script;
     
     mode755 = S_IRWXU|S_IRGRP|S_IXGRP|S_IROTH|S_IXOTH;
-    
-    // backup I/O parameters
-    terminalBack = options_.terminal;
-    outputBack   = options_.output;
-    commandBack  = plotCommand_;
 
     // generate directory
     if (mkdir(dirName))
@@ -703,12 +704,20 @@ void Plot::save(string dirName)
         LATAN_ERROR(Io, "impossible to create directory '" + dirName + "'");
     }
     
+    // backup I/O parameters
+    terminalBack = options_.terminal;
+    outputBack   = options_.output;
+    commandBack  = plotCommand_;
+
     // save PDF
-    options_.terminal = "pdf";
+    if (savePdf)
-    options_.output   = dirName + "/plot.pdf";
+    {
-    display();
+        options_.terminal = "pdf";
-    options_.terminal = terminalBack;
+        options_.output   = dirName + "/plot.pdf";
-    options_.output   = outputBack;
+        display();
+        options_.terminal = terminalBack;
+        options_.output   = outputBack;
+    }
     
     // save script and datafiles
     for (unsigned int i = 0; i < tmpFileName_.size(); ++i)
@@ -776,18 +785,25 @@ ostream & Latan::operator<<(ostream &out, const Plot &plot)
         out << "yMin = " << plot.options_.scale[y].min << endl;
         out << "yMax = " << plot.options_.scale[y].max << endl;
     }
-    if (!plot.options_.title.empty())
+    out << "unset xrange" << endl;
-    {
-        out << "set title '" << plot.options_.title << "'" << endl;
-    }
     if (plot.options_.scaleMode[x] & Plot::Scale::manual)
     {
         out << "set xrange [xMin:xMax]" << endl;
     }
+    else
+    {
+        out << "set xrange [:]" << endl;
+    }
+    out << "unset yrange" << endl;
     if (plot.options_.scaleMode[y] & Plot::Scale::manual)
     {
         out << "set yrange [yMin:yMax]" << endl;
     }
+    else
+    {
+        out << "set yrange [:]" << endl;
+    }
+    out << "unset log" << endl;
     if (plot.options_.scaleMode[x] & Plot::Scale::log)
     {
         out << "set log x" << endl;

lib/Core/Plot.hpp

@@ -1,7 +1,7 @@
 /*
  * Plot.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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
@@ -377,13 +377,13 @@ public:
     Plot & operator<<(PlotModifier &&modifier);
     // plot parsing and output
     void display(void);
-    void save(std::string dirName);
+    void save(std::string dirName, bool savePdf = true);
     friend std::ostream & operator<<(std::ostream &out, const Plot &plot);
     // plot reset
     void reset(void);
-private:
     // find gnuplot
-    void getProgramPath(void);
+    std::string getProgramPath(void);
+private:
     // default options
     void initOptions(void);
 private:

lib/Core/Utilities.cpp

@@ -1,7 +1,7 @@
 /*
  * Utilities.cpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2015 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

lib/Core/Utilities.hpp

@@ -1,7 +1,7 @@
 /*
  * Utilities.hpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2015 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

lib/Core/stdincludes.hpp

@@ -1,7 +1,7 @@
 /*
  * stdincludes.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/CompiledFunction.cpp

@@ -1,7 +1,7 @@
 /*
  * CompiledFunction.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/CompiledFunction.hpp

@@ -1,7 +1,7 @@
 /*
  * CompiledFunction.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/CompiledModel.cpp

@@ -1,7 +1,7 @@
 /*
  * CompiledModel.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/CompiledModel.hpp

@@ -1,7 +1,7 @@
 /*
  * CompiledModel.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/Function.cpp

@@ -1,7 +1,7 @@
 /*
  * Function.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/Function.hpp

@@ -1,7 +1,7 @@
 /*
  * Function.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/Model.cpp

@@ -1,7 +1,7 @@
 /*
  * Model.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/Model.hpp

@@ -1,7 +1,7 @@
 /*
  * Model.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/TabFunction.cpp

@@ -1,7 +1,7 @@
 /*
  * TabFunction.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Functional/TabFunction.hpp

@@ -1,7 +1,7 @@
 /*
  * TabFunction.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Global.cpp

@@ -1,7 +1,7 @@
 /*
  * Global.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Global.hpp

@@ -1,7 +1,7 @@
 /*
  * Global.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/AsciiFile.cpp

@@ -1,7 +1,7 @@
 /*
  * AsciiFile.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/AsciiFile.hpp

@@ -1,7 +1,7 @@
 /*
  * AsciiFile.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/AsciiLexer.lpp

@@ -1,7 +1,7 @@
 /*
  * AsciiLexer.lpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/AsciiParser.ypp

@@ -1,7 +1,7 @@
 /*
  * AsciiParser.ypp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/File.cpp

@@ -1,7 +1,7 @@
 /*
  * File.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/File.hpp

@@ -1,7 +1,7 @@
 /*
  * File.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/Hdf5File.cpp

@@ -1,7 +1,7 @@
 /*
  * Hdf5File.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
+ * Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs
  *
  * 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

lib/Io/Hdf5File.hpp

@@ -1,7 +1,7 @@
 /*
  * Hdf5File.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
+ * Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs
  *
  * 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

lib/Io/Io.cpp

@@ -1,7 +1,7 @@
 /*
  * Io.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/Io.hpp

@@ -1,7 +1,7 @@
 /*
  * Io.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/IoObject.hpp

@@ -1,7 +1,7 @@
 /*
  * IoObject.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Io/XmlReader.cpp

@@ -1,7 +1,7 @@
 /*
  * XmlReader.cpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2014 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

lib/Io/XmlReader.hpp

@@ -1,7 +1,7 @@
 /*
  * XmlReader.hpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2015 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

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/Numerical/Derivative.cpp

@@ -1,7 +1,7 @@
 /*
  * Derivative.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/Derivative.hpp

@@ -1,7 +1,7 @@
 /*
  * Derivative.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/FFT.hpp

@@ -1,7 +1,7 @@
 /*
  * FFT.hpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2017 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

lib/Numerical/GslFFT.cpp

@@ -1,7 +1,7 @@
 /*
  * GslFFT.cpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2017 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -43,6 +43,7 @@ void GslFFT::resize(const Index size)
 {
     if (size_ != size)
     {
+        clear();
         size_      = size;
         wavetable_ = gsl_fft_complex_wavetable_alloc(size_);
         workspace_ = gsl_fft_complex_workspace_alloc(size_);

lib/Numerical/GslFFT.hpp

@@ -1,7 +1,7 @@
 /*
  * GslFFT.hpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2017 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

lib/Numerical/GslHybridRootFinder.cpp

@@ -1,7 +1,7 @@
 /*
  * GslHybridRootFinder.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/GslHybridRootFinder.hpp

@@ -1,7 +1,7 @@
 /*
  * GslHybridRootFinder.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/GslMinimizer.cpp

@@ -1,7 +1,7 @@
 /*
  * GslMinimizer.cpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

lib/Numerical/GslMinimizer.hpp

@@ -1,7 +1,7 @@
 /*
  * GslMinimizer.hpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by

lib/Numerical/GslQagsIntegrator.cpp

@@ -1,7 +1,7 @@
 /*
  * GslQagsIntegrator.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/GslQagsIntegrator.hpp

@@ -1,7 +1,7 @@
 /*
  * GslQagsIntegrator.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/Integrator.hpp

@@ -1,7 +1,7 @@
 /*
  * Integrator.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/Minimizer.cpp

@@ -1,7 +1,7 @@
 /*
  * Minimizer.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/Minimizer.hpp

@@ -1,7 +1,7 @@
 /*
  * Minimizer.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/MinuitMinimizer.cpp

@@ -1,7 +1,7 @@
 /*
  * MinuitMinimizer.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/MinuitMinimizer.hpp

@@ -1,7 +1,7 @@
 /*
  * MinuitMinimizer.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/NloptMinimizer.cpp

@@ -1,7 +1,7 @@
 /*
  * NloptMinimizer.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/NloptMinimizer.hpp

@@ -1,7 +1,7 @@
 /*
  * NloptMinimizer.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/RootFinder.cpp

@@ -1,7 +1,7 @@
 /*
  * RootFinder.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/RootFinder.hpp

@@ -1,7 +1,7 @@
 /*
  * RootFinder.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/Solver.cpp

@@ -1,7 +1,7 @@
 /*
  * Solver.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Numerical/Solver.hpp

@@ -1,7 +1,7 @@
 /*
  * Solver.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Physics/CorrelatorFitter.cpp

@@ -0,0 +1,417 @@
+/*
+ * 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 __dumb, const double *p __dumb)
+    {
+        return p[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[1] + p[0]*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;
+}
+
+/******************************************************************************
+ *                         Correlator utilities                               *
+ ******************************************************************************/
+DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts)
+{
+    if (ts != 0)
+    {
+        const Index nt  = c[central].rows();
+        DMatSample  buf = c;
+
+        FOR_STAT_ARRAY(buf, s)
+        {
+            for (Index t = 0; t < nt; ++t)
+            {
+                buf[s]((t - ts + nt)%nt) = c[s](t);
+            }
+        }
+
+        return buf;
+    }
+    else
+    {
+        return c;
+    }
+}
+
+DMatSample CorrelatorUtils::fold(const DMatSample &c)
+{
+    const Index nt  = c[central].rows();
+    DMatSample  buf = c;
+
+    FOR_STAT_ARRAY(buf, s)
+    {
+        for (Index t = 0; t < nt; ++t)
+        {
+            buf[s](t) = 0.5*(c[s](t) + c[s]((nt - t) % nt));
+        }
+    }
+
+    return buf;
+}
+
+DMatSample CorrelatorUtils::fourierTransform(const DMatSample &c, FFT &fft, 
+                                             const unsigned int dir)
+{
+    const Index nSample   = c.size();
+    const Index nt        = c[central].rows();
+    bool        isComplex = (c[central].cols() > 1);
+    CMatSample  buf(nSample, nt, 1);
+    DMatSample  out(nSample, nt, 2);
+
+    fft.resize(nt);
+    FOR_STAT_ARRAY(buf, s)
+    {
+        buf[s].real() = c[s].col(0);
+        if (isComplex)
+        {
+            buf[s].imag() = c[s].col(1);
+        }
+        else
+        {
+            buf[s].imag() = DVec::Constant(nt, 0.);
+        }
+        fft(buf[s], dir);
+        out[s].col(0) = buf[s].real();
+        out[s].col(1) = buf[s].imag();
+    }
+
+    return out;    
+}
+
+/******************************************************************************
+ *                      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,104 @@
+/*
+ * 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/Numerical/FFT.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 utilities                               *
+ ******************************************************************************/
+namespace CorrelatorUtils
+{
+    DMatSample shift(const DMatSample &c, const Index ts);
+    DMatSample fold(const DMatSample &c);
+    DMatSample fourierTransform(const DMatSample &c, FFT &fft, 
+                                const unsigned int dir = FFT::Forward);
+};
+
+/******************************************************************************
+ *                       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 = 1; t < nt; ++t)
+        {
+            em(t - 1) = 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_

lib/Statistics/Dataset.hpp

@@ -1,7 +1,7 @@
 /*
  * Dataset.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/FitInterface.cpp

@@ -1,7 +1,7 @@
 /*
  * FitInterface.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/FitInterface.hpp

@@ -1,7 +1,7 @@
 /*
  * FitInterface.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/Histogram.cpp

@@ -1,7 +1,7 @@
 /*
  * Histogram.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/Histogram.hpp

@@ -1,7 +1,7 @@
 /*
  * Histogram.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/MatSample.hpp

@@ -1,7 +1,7 @@
 /*
  * MatSample.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/Random.cpp

@@ -1,7 +1,7 @@
 /*
  * Random.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/Random.hpp

@@ -1,7 +1,7 @@
 /*
  * Random.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2019 Antonin Portelli
+ * 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

lib/Statistics/StatArray.cpp

@@ -1,7 +1,7 @@
 /*
  * StatArray.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/StatArray.hpp

@@ -1,7 +1,7 @@
 /*
  * StatArray.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/XYSampleData.cpp

@@ -1,7 +1,7 @@
 /*
  * XYSampleData.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/XYSampleData.hpp

@@ -1,7 +1,7 @@
 /*
  * XYSampleData.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/XYStatData.cpp

@@ -1,7 +1,7 @@
 /*
  * XYStatData.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/Statistics/XYStatData.hpp

@@ -1,7 +1,7 @@
 /*
  * XYStatData.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

lib/includes.hpp

@@ -1,7 +1,7 @@
 /*
  * includes.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

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");
@@ -114,115 +105,25 @@ int main(int argc, char *argv[])
     DMatSample tmp, corr;
     Index      nSample, nt;
     
-    tmp     = Io::load<DMatSample>(corrFileName);
+    corr    = Io::load<DMatSample>(corrFileName);
-    nSample = tmp.size();
+    nSample = corr.size();
-    nt      = tmp[central].rows();
+    nt      = corr[central].rows();
-    tmp     = tmp.block(0, 0, nt, 1);
+    corr    = corr.block(0, 0, nt, 1);
-    corr    = tmp;
+    corr    = CorrelatorUtils::shift(corr, shift);
-    FOR_STAT_ARRAY(corr, s)
-    {
-        for (Index t = 0; t < nt; ++t)
-        {
-            corr[s]((t - shift + nt)%nt) = tmp[s](t);
-        }
-    }
     if (fold)
     {
-        tmp = corr;
+        corr = CorrelatorUtils::fold(corr);
-        FOR_STAT_ARRAY(corr, s)
-        {
-            for (Index t = 0; t < nt; ++t)
-            {
-                corr[s](t) = 0.5*(tmp[s](t) + tmp[s]((nt - t) % nt));
-            }
-        }
     }
     
-    // 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 +141,52 @@ 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 if(modelPar.type == CorrelatorType::linear)
+        {
+            globMin.setLowLimit(p, -10.*fabs(init(p)));
+            locMin.setLowLimit(p, -10.*fabs(init(p)));
+            globMin.setHighLimit(p, 10.*init(p));
+            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 +194,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 +225,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 +278,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 +301,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 +316,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 +337,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 +372,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 +382,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 +392,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 +402,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 != "")
             {

update_eigen.sh

@@ -6,7 +6,7 @@ if (( $# != 1 )); then
 fi
 ARC=$1
 
-INITDIR=`pwd`
+INITDIR=$(pwd -P)
 rm -rf lib/Eigen
 ARCDIR=`tar -tf ${ARC} | head -n1 | sed -e 's@/.*@@'`
 tar -xf ${ARC}

utils/resample.cpp

@@ -1,7 +1,7 @@
 /*
  * resample.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

utils/sample-combine.cpp

@@ -1,7 +1,7 @@
 /*
  * sample_combine.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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

utils/sample-element.cpp

@@ -1,7 +1,7 @@
 /*
  * sample-element.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
+ * Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs
  *
  * 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

utils/sample-fake.cpp

@@ -1,7 +1,7 @@
 /*
  * make_fake_sample.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * 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