Merge 269d0c338e into 05138baa08

fix of potential conflict with std::size
Update build-macos.yml
2025-10-27 13:39:32 +00:00 · 2023-06-20 15:48:04 +01:00 · 2023-06-20 15:47:58 +01:00 · 2023-03-20 14:04:17 +00:00 · 2022-08-03 14:24:36 +01:00 · 2022-08-03 14:24:20 +01:00
122 changed files with 3104 additions and 911 deletions
--- a/.github/workflows/build-macos.yml
+++ b/.github/workflows/build-macos.yml
@@ -0,0 +1,26 @@
+name: Build macOS
+
+on: [push, workflow_dispatch]
+
+jobs:
+  build:
+
+    runs-on: macos-11
+
+    steps:
+    - name: Checkout
+      uses: actions/checkout@v2
+    - name: Install basic dependencies
+      run: brew install automake autoconf libtool bison flex
+    - name: Build dependencies
+      shell: bash
+      run: |
+        export PATH=/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:${PATH}
+        cd ci-scripts
+        ./install-deps.sh prefix 6
+    - name: Build LatAnalyze
+      shell: bash
+      run: |
+        export PATH=/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:${PATH}
+        cd ci-scripts
+        ./install-latan.sh prefix 6
--- a/.github/workflows/build-ubuntu.yml
+++ b/.github/workflows/build-ubuntu.yml
@@ -0,0 +1,26 @@
+name: Build Ubuntu
+
+on: [push]
+
+jobs:
+  build:
+
+    runs-on: ubuntu-20.04
+
+    steps:
+    - name: Checkout
+      uses: actions/checkout@v2
+    - name: Install basic dependencies
+      run: |
+        sudo bash -c "$(wget -O - https://apt.llvm.org/llvm.sh)"
+        sudo apt install cmake bison flex
+    - name: Build dependencies
+      shell: bash
+      run: |
+        cd ci-scripts
+        CC=clang CXX=clang++ ./install-deps.sh prefix 6
+    - name: Build LatAnalyze
+      shell: bash
+      run: |
+        cd ci-scripts
+        CC=clang CXX=clang++ ./install-latan.sh prefix 6
--- a/.gitignore
+++ b/.gitignore
@@ -16,6 +16,7 @@ autom4te.cache/*
 *.in~
 config.h*
 configure
+configure~
 .buildutils/*
 aclocal.m4

--- a/Readme.md
+++ b/Readme.md
@@ -1,21 +1,6 @@
 # LatAnalyze 

-License: GNU General Public License v3
-
-<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>
+[![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)

 ## 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.**
--- a/bootstrap.sh
+++ b/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
--- a/ci-scripts/install-deps.sh
+++ b/ci-scripts/install-deps.sh
@@ -1,15 +1,16 @@
 #!/usr/bin/env bash

-if (( $# != 1 )); then
-  echo "usage: `basename $0` <prefix>" 1>&2
+if (( $# != 2 )); then
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
  exit 1
 fi
 PREFIX=$1
+NTASKS=$2

 set -ex
 mkdir -p local/build
 for d in gsl nlopt minuit hdf5; do
  if [ ! -e local/.built.${d} ]; then
-    ./install-${d}.sh ${PREFIX}
+    ./install-${d}.sh ${PREFIX} ${NTASKS}
  fi
 done
--- a/ci-scripts/install-gsl.sh
+++ b/ci-scripts/install-gsl.sh
@@ -2,21 +2,25 @@

 NAME='gsl-2.6'

-if (( $# != 1 )); then
-  echo "usage: `basename $0` <prefix>" 1>&2
+if (( $# != 2 )); then
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
  exit 1
 fi
 PREFIX=$1
+NTASKS=$2

 set -ex
-INITDIR=`pwd`
-cd local/build
+INITDIR=$(pwd -P)
+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
 cd ${NAME}/build
 ../configure --prefix=${PREFIX}
-make -j4 
+make -j${NTASKS}
 make install
 cd ${INITDIR}/local
 touch .built.gsl
--- a/ci-scripts/install-hdf5.sh
+++ b/ci-scripts/install-hdf5.sh
@@ -1,22 +1,26 @@
 #!/usr/bin/env bash

-NAME='hdf5-1.10.5'
+NAME='hdf5-1.10.8'

-if (( $# != 1 )); then
-  echo "usage: `basename $0` <prefix>" 1>&2
+if (( $# != 2 )); then
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
  exit 1
 fi
 PREFIX=$1
+NTASKS=$2

 set -ex
-INITDIR=`pwd`
-cd local/build
+INITDIR=$(pwd -P)
+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
 cd ${NAME}/build
 ../configure --prefix=${PREFIX} --enable-cxx
-make -j4 
+make -j${NTASKS}
 make install
 cd ${INITDIR}/local
 touch .built.hdf5
--- a/ci-scripts/install-latan.sh
+++ b/ci-scripts/install-latan.sh
@@ -1,18 +1,23 @@
 #!/usr/bin/env bash

-if (( $# != 1 )); then
-  echo "usage: `basename $0` <prefix>" 1>&2
+if (( $# != 2 )); then
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
  exit 1
 fi
 PREFIX=$1
-OS=$2
+NTASKS=$2

 set -ex
-./install-deps.sh ${PREFIX}
+INITDIR=$(pwd -P)
+mkdir -p ${PREFIX}
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}
+./install-deps.sh ${PREFIX} ${NTASKS}
 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"
-make -j4
+../configure --prefix=${PREFIX} --with-minuit=${PREFIX} --with-nlopt=${PREFIX} --with-hdf5=${PREFIX} --with-gsl=${PREFIX} CXXFLAGS="${CXXFLAGS} -O3 -march=native -mtune=native"
+make -j${NTASKS}
 make install
--- a/ci-scripts/install-minuit.sh
+++ b/ci-scripts/install-minuit.sh
@@ -1,22 +1,25 @@
 #!/usr/bin/env bash

-NAME='Minuit2-5.34.14'

-if (( $# != 1 )); then
-  echo "usage: `basename $0` <prefix>" 1>&2
+if (( $# != 2 )); then
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
  exit 1
 fi
 PREFIX=$1
+NTASKS=$2

 set -ex
-INITDIR=`pwd`
-cd 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
-cd ${NAME}/build
-../configure --prefix=${PREFIX} --disable-openmp
-make -j4 
+INITDIR=$(pwd -P)
+mkdir -p ${PREFIX}
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}/local/build
+rm -rf root
+git clone https://github.com/root-project/root.git
+cd root/math/minuit2/
+mkdir build; cd build
+cmake .. -Dminuit2_standalone=ON -DCMAKE_INSTALL_PREFIX=${PREFIX}
+make -j${NTASKS}
 make install
 cd ${INITDIR}/local
 touch .built.minuit
--- a/ci-scripts/install-nlopt.sh
+++ b/ci-scripts/install-nlopt.sh
@@ -2,22 +2,26 @@

 NAME='2.6.1'

-if (( $# != 1 )); then
-  echo "usage: `basename $0` <prefix>" 1>&2
+if (( $# != 2 )); then
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
  exit 1
 fi
 PREFIX=$1
+NTASKS=$2

 set -ex
-INITDIR=`pwd`
-cd local/build
+INITDIR=$(pwd -P)
+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} ..
-make -j4 
+cmake -DCMAKE_INSTALL_PREFIX=${PREFIX} -DCMAKE_BUILD_WITH_INSTALL_NAME_DIR=TRUE -DCMAKE_INSTALL_NAME_DIR="${PREFIX}/lib" ..
+make -j${NTASKS}
 make install
 cd ${INITDIR}/local
 touch .built.nlopt
--- a/configure.ac
+++ b/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-dev],[antonin.portelli@me.com],[LatAnalyze])
 AC_CONFIG_AUX_DIR([.buildutils])
 AC_CONFIG_SRCDIR([lib/Global.cpp])
 AC_CONFIG_SRCDIR([utils/sample_read.cpp])
@@ -36,7 +36,7 @@ AC_ARG_WITH([gsl],
 AC_ARG_WITH([minuit],
    [AS_HELP_STRING([--with-minuit=prefix],
 		[try this for a non-standard install prefix of the Minuit2 library])],
-    [AM_CXXFLAGS="$AM_CXXFLAGS -I$with_minuit/include"]
+    [AM_CXXFLAGS="$AM_CXXFLAGS -I$with_minuit/include -I$with_minuit/include/Minuit2 -I$with_minuit/include/Fit"]
 	[AM_LDFLAGS="$AM_LDFLAGS -L$with_minuit/lib"])
 AC_ARG_WITH([nlopt],
    [AS_HELP_STRING([--with-nlopt=prefix],
@@ -74,6 +74,7 @@ CXXFLAGS_CPY=$CXXFLAGS
 LDFLAGS_CPY=$LDFLAGS
 CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS"
 LDFLAGS="$AM_LDFLAGS $LDFLAGS"
+AC_CHECK_LIB([pthread],[pthread_create],[],[AC_MSG_ERROR([pthread library not found])])
 AC_CHECK_LIB([m],[cos],[],[AC_MSG_ERROR([libm library not found])])
 AC_CHECK_LIB([gslcblas],[cblas_dgemm],[],
             [AC_MSG_ERROR([GSL CBLAS library not found])])
@@ -90,10 +91,10 @@ AC_CHECK_LIB([hdf5_cpp],[H5Fopen],
    [AC_MSG_ERROR([HDF5 library not found])], [-lhdf5])
 SAVED_LDFLAGS=$LDFLAGS
 LDFLAGS="$LDFLAGS -lMinuit2"
-AC_MSG_CHECKING([for ROOT::Minuit2::BasicMinimumError in -lMinuit2]);
+AC_MSG_CHECKING([for ROOT::Minuit2::VariableMetricMinimizer in -lMinuit2]);
 AC_LINK_IFELSE(
-	[AC_LANG_PROGRAM([#include <Minuit2/BasicMinimumError.h>],
-	[ROOT::Minuit2::BasicMinimumError dummy(0)])],
+	[AC_LANG_PROGRAM([#include <Minuit2/VariableMetricMinimizer.h>],
+	[ROOT::Minuit2::VariableMetricMinimizer dummy()])],
 	[LIBS="$LIBS -lMinuit2"]
 	[AC_DEFINE([HAVE_MINUIT2],
 				[1],
@@ -103,6 +104,20 @@ AC_LINK_IFELSE(
 	[have_minuit=false]
    [AC_MSG_RESULT([no])])
 AM_CONDITIONAL([HAVE_MINUIT], [test x$have_minuit = xtrue])
+LDFLAGS="$LDFLAGS -lMinuit2Math"
+AC_MSG_CHECKING([for ROOT::Math::MinimizerOptions in -lMinuit2Math]);
+AC_LINK_IFELSE(
+	[AC_LANG_PROGRAM([#include <Minuit2/Math/MinimizerOptions.h>],
+	[ROOT::Math::MinimizerOptions dummy()])],
+	[LIBS="$LIBS -lMinuit2Math"]
+	[AC_DEFINE([HAVE_MINUIT2MATH],
+				[1],
+				[Define to 1 if you have the `Minuit2Math' library (-lMinuit2Math).])]
+    [have_minuitmath=true]
+	[AC_MSG_RESULT([yes])],
+	[have_minuitmath=false]
+    [AC_MSG_RESULT([no])])
+AM_CONDITIONAL([HAVE_MINUITMATH], [test x$have_minuit = xtrue])
 LDFLAGS=$SAVED_LDFLAGS
 CXXFLAGS=$CXXFLAGS_CPY
 LDFLAGS=$LDFLAGS_CPY
--- a/eigen-3.3.7.tar.bz2
+++ b/eigen-3.3.7.tar.bz2
--- a/eigen-3.3.8.tar.bz2
+++ b/eigen-3.3.8.tar.bz2
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -9,6 +9,7 @@ endif
 noinst_PROGRAMS =           \
    exCompiledDoubleFunction\
    exDerivative            \
+    exDWT                   \
    exFit                   \
    exFitSample             \
    exIntegrator            \
@@ -19,7 +20,8 @@ noinst_PROGRAMS =           \
    exPlot                  \
    exPValue                \
    exRand                  \
-    exRootFinder
+    exRootFinder            \
+    exThreadPool

 exCompiledDoubleFunction_SOURCES  = exCompiledDoubleFunction.cpp
 exCompiledDoubleFunction_CXXFLAGS = $(COM_CXXFLAGS)
@@ -29,6 +31,10 @@ exDerivative_SOURCES              = exDerivative.cpp
 exDerivative_CXXFLAGS             = $(COM_CXXFLAGS)
 exDerivative_LDFLAGS              = -L../lib/.libs -lLatAnalyze

+exDWT_SOURCES                     = exDWT.cpp
+exDWT_CXXFLAGS                    = $(COM_CXXFLAGS)
+exDWT_LDFLAGS                     = -L../lib/.libs -lLatAnalyze
+
 exFit_SOURCES                     = exFit.cpp
 exFit_CXXFLAGS                    = $(COM_CXXFLAGS)
 exFit_LDFLAGS                     = -L../lib/.libs -lLatAnalyze
@@ -73,4 +79,8 @@ exRootFinder_SOURCES              = exRootFinder.cpp
 exRootFinder_CXXFLAGS             = $(COM_CXXFLAGS)
 exRootFinder_LDFLAGS              = -L../lib/.libs -lLatAnalyze

+exThreadPool_SOURCES              = exThreadPool.cpp
+exThreadPool_CXXFLAGS             = $(COM_CXXFLAGS)
+exThreadPool_LDFLAGS              = -L../lib/.libs -lLatAnalyze
+
 ACLOCAL_AMFLAGS = -I .buildutils/m4
--- a/examples/exDWT.cpp
+++ b/examples/exDWT.cpp
@@ -0,0 +1,28 @@
+#include <LatAnalyze/Numerical/DWT.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(void)
+{
+    DVec                  data, dataRec;
+    vector<DWT::DWTLevel> dataDWT;
+    DWT                   dwt(DWTFilters::db3);
+
+    cout << "-- random data" << endl;
+    data.setRandom(16);
+    cout << data.transpose() << endl;
+    cout << "-- compute Daubechies 3 DWT" << endl;
+    dataDWT = dwt.forward(data, 4);
+    for (unsigned int l = 0; l < dataDWT.size(); ++l)
+    {
+        cout << "* level " << l << endl;
+        cout << "L= " << dataDWT[l].first.transpose() << endl;
+        cout << "H= " << dataDWT[l].second.transpose() << endl;
+    }
+    cout << "-- check inverse DWT" << endl;
+    dataRec = dwt.backward(dataDWT);
+    cout << "rel diff = " << 2.*(data - dataRec).norm()/(data + dataRec).norm() << endl;
+
+    return EXIT_SUCCESS;
+}
--- a/examples/exRand.cpp
+++ b/examples/exRand.cpp
@@ -7,7 +7,7 @@
 using namespace std;
 using namespace Latan;

-constexpr Index  size          = 8;
+constexpr Index  n             = 8;
 constexpr Index  nDraw         = 20000;
 constexpr Index  nSample       = 2000;
 const     string stateFileName = "exRand.seed";
@@ -40,14 +40,14 @@ int main(void)
    p << PlotFunction(compile("return exp(-x_0^2/2)/sqrt(2*pi);", 1), -5., 5.);
    p.display();

-    DMat       var(size, size);
-    DVec       mean(size);
-    DMatSample sample(nSample, size, 1);
+    DMat       var(n, n);
+    DVec       mean(n);
+    DMatSample sample(nSample, n, 1);

    cout << "-- generating " << nSample << " Gaussian random vectors..." << endl;
-    var   = DMat::Random(size, size);
+    var   = DMat::Random(n, n);
    var  *= var.adjoint();
-    mean  = DVec::Random(size); 
+    mean  = DVec::Random(n); 
    RandomNormal mgauss(mean, var, rd());
    sample[central] = mgauss();
    FOR_STAT_ARRAY(sample, s)
--- a/examples/exThreadPool.cpp
+++ b/examples/exThreadPool.cpp
@@ -0,0 +1,29 @@
+#include <LatAnalyze/Core/ThreadPool.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(void)
+{
+    ThreadPool pool;
+
+    cout << "Using " << pool.getThreadNum() << " threads" << endl;
+    for (unsigned int i = 1; i <= 20; ++i)
+    {
+        pool.addJob([i, &pool](void)
+        {
+            pool.critical([i](void)
+            {
+                cout << "job " << i << " wait for " << i*100 << " ms" << endl;
+            });
+            this_thread::sleep_for(chrono::milliseconds(i*100));
+            pool.critical([i](void)
+            {
+                cout << "job " << i << " done" << endl;
+            });
+        });
+    }
+    pool.terminate();
+
+    return EXIT_SUCCESS;
+}
--- a/lib/Core/Eigen.hpp
+++ b/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
--- a/lib/Core/EigenPlugin.hpp
+++ b/lib/Core/EigenPlugin.hpp
@@ -17,7 +17,7 @@
 * along with LatAnalyze.  If not, see <http://www.gnu.org/licenses/>.
 */

-Derived pInverse(const double tolerance = 1.0e-10)
+Derived pInverse(const double tolerance = 1.0e-10) const
 {
    auto         svd   = jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV);
    const auto   u     = svd.matrixU();
@@ -52,7 +52,7 @@ Derived pInverse(const double tolerance = 1.0e-10)
    return v*s.asDiagonal()*u.transpose();
 }

-Derived singularValues(void)
+Derived singularValues(void) const
 {
    auto svd = jacobiSvd();
    
--- a/lib/Core/Exceptions.cpp
+++ b/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
--- a/lib/Core/Exceptions.hpp
+++ b/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
--- a/lib/Core/Mat.cpp
+++ b/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
--- a/lib/Core/Mat.hpp
+++ b/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
--- a/lib/Core/Math.cpp
+++ b/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
@@ -29,7 +29,8 @@ using namespace Latan;
 ******************************************************************************/
 DMat MATH_NAMESPACE::varToCorr(const DMat &var)
 {
-    DMat res = var, invDiag = res.diagonal();
+    DMat res = var;
+    DVec invDiag = res.diagonal();
    
    invDiag = invDiag.cwiseInverse().cwiseSqrt();
    res     = (invDiag*invDiag.transpose()).cwiseProduct(res);
@@ -37,6 +38,28 @@ DMat MATH_NAMESPACE::varToCorr(const DMat &var)
    return res;
 }

+DMat MATH_NAMESPACE::corrToVar(const DMat &corr, const DVec &varDiag)
+{
+    DMat res = corr;
+    DVec varSqrtDiag = varDiag.cwiseSqrt();
+    
+    res = (varSqrtDiag*varSqrtDiag.transpose()).cwiseProduct(res);
+
+    return res;
+}
+
+double MATH_NAMESPACE::svdDynamicRange(const DMat &mat)
+{
+    DVec s = mat.singularValues();
+
+    return s.maxCoeff()/s.minCoeff();
+}
+
+double MATH_NAMESPACE::svdDynamicRangeDb(const DMat &mat)
+{
+    return 10.*log10(svdDynamicRange(mat));
+}
+
 /******************************************************************************
 *                             Standard C functions                           *
 ******************************************************************************/
--- a/lib/Core/Math.hpp
+++ b/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
@@ -70,6 +70,11 @@ namespace MATH_NAMESPACE
    
    // convert variance matrix to correlation matrix
    DMat varToCorr(const DMat &var);
+    DMat corrToVar(const DMat &corr, const DVec &varDiag);
+
+    // matrix SVD dynamic range
+    double svdDynamicRange(const DMat &mat);
+    double svdDynamicRangeDb(const DMat &mat);
    
    // Constants
    constexpr double pi  = 3.1415926535897932384626433832795028841970;
--- a/lib/Core/MathInterpreter.cpp
+++ b/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
@@ -217,146 +217,25 @@ void RunContext::reset(void)
 #define CODE_WIDTH 6
 #define CODE_MOD   setw(CODE_WIDTH) << left

-// Instruction operator ////////////////////////////////////////////////////////
-ostream &Latan::operator<<(ostream& out, const Instruction& ins)
+auto readConstant(Program::const_iterator ip)
+    -> std::tuple<double, Program::const_iterator>
 {
-    ins.print(out);
-    
-    return out;
+    double value = 0.0;
+    std::copy(ip, ip + sizeof(double), reinterpret_cast<std::uint8_t*>(&value));
+
+    return std::make_tuple(value, ip + sizeof(double));
 }

-// Push constructors ///////////////////////////////////////////////////////////
-Push::Push(const double val)
-: type_(ArgType::Constant)
-, val_(val)
-, address_(0)
-, name_("")
-{}
-
-Push::Push(const unsigned int address, const string &name)
-: type_(ArgType::Variable)
-, val_(0.0)
-, address_(address)
-, name_(name)
-{}
-
-// Push execution //////////////////////////////////////////////////////////////
-void Push::operator()(RunContext &context) const
+auto readAddress(Program::const_iterator ip)
+    -> std::tuple<unsigned int, Program::const_iterator>
 {
-    if (type_ == ArgType::Constant)
-    {
-        context.stack().push(val_);
-    }
-    else
-    {
-        context.stack().push(context.getVariable(address_));
-    }
-    context.incrementInsIndex();
+    unsigned int address = 0.0;
+    const auto end = ip + sizeof(unsigned int);
+    std::copy(ip, end, reinterpret_cast<std::uint8_t*>(&address));
+
+    return std::make_tuple(address, end);
 }

-// Push print //////////////////////////////////////////////////////////////////
-void Push::print(ostream &out) const
-{
-    out << CODE_MOD << "push";
-    if (type_ == ArgType::Constant)
-    {
-        out << CODE_MOD << val_;
-    }
-    else
-    {
-        out << CODE_MOD << name_ << " @v" << address_;
-    }
-}
-
-// Pop constructor /////////////////////////////////////////////////////////////
-Pop::Pop(const unsigned int address, const string &name)
-: address_(address)
-, name_(name)
-{}
-
-// Pop execution ///////////////////////////////////////////////////////////////
-void Pop::operator()(RunContext &context) const
-{
-    if (!name_.empty())
-    {
-        context.setVariable(address_, context.stack().top());
-    }
-    context.stack().pop();
-    context.incrementInsIndex();
-}
-
-// Pop print ///////////////////////////////////////////////////////////////////
-void Pop::print(ostream &out) const
-{
-    out << CODE_MOD << "pop" << CODE_MOD << name_ << " @v" << address_;
-}
-
-// Store constructor ///////////////////////////////////////////////////////////
-Store::Store(const unsigned int address, const string &name)
-: address_(address)
-, name_(name)
-{}
-
-// Store execution /////////////////////////////////////////////////////////////
-void Store::operator()(RunContext &context) const
-{
-    if (!name_.empty())
-    {
-        context.setVariable(address_, context.stack().top());
-    }
-    context.incrementInsIndex();
-}
-
-// Store print /////////////////////////////////////////////////////////////////
-void Store::print(ostream &out) const
-{
-    out << CODE_MOD << "store" << CODE_MOD << name_ << " @v" << address_;
-}
-
-// Call constructor ////////////////////////////////////////////////////////////
-Call::Call(const unsigned int address, const string &name)
-: address_(address)
-, name_(name)
-{}
-
-// Call execution //////////////////////////////////////////////////////////////
-void Call::operator()(RunContext &context) const
-{
-    context.stack().push((*context.getFunction(address_))(context.stack()));
-    context.incrementInsIndex();
-}
-
-// Call print //////////////////////////////////////////////////////////////////
-void Call::print(ostream &out) const
-{
-    out << CODE_MOD << "call" << CODE_MOD << name_ << " @f" << address_;
-}
-
-// Math operations /////////////////////////////////////////////////////////////
-#define DEF_OP(name, nArg, exp, insName)\
-void name::operator()(RunContext &context) const\
-{\
-    double x[nArg];\
-    for (int i = 0; i < nArg; ++i)\
-    {\
-        x[nArg-1-i] = context.stack().top();\
-        context.stack().pop();\
-    }\
-    context.stack().push(exp);\
-    context.incrementInsIndex();\
-}\
-void name::print(ostream &out) const\
-{\
-    out << CODE_MOD << insName;\
-}
-
-DEF_OP(Neg, 1, -x[0],          "neg")
-DEF_OP(Add, 2, x[0] + x[1],    "add")
-DEF_OP(Sub, 2, x[0] - x[1],    "sub")
-DEF_OP(Mul, 2, x[0]*x[1],      "mul")
-DEF_OP(Div, 2, x[0]/x[1],      "div")
-DEF_OP(Pow, 2, pow(x[0],x[1]), "pow")
-
 /******************************************************************************
 *                        ExprNode implementation                             *
 ******************************************************************************/
@@ -442,29 +321,35 @@ ostream &Latan::operator<<(ostream &out, const ExprNode &n)
    return out;
 }

-#define PUSH_INS(program, type, ...)\
-program.push_back(unique_ptr<type>(new type(__VA_ARGS__)))
-#define GET_ADDRESS(address, table, name)\
-try\
-{\
-    address = (table).at(name);\
-}\
-catch (out_of_range)\
-{\
-    address         = (table).size();\
-    (table)[(name)] = address;\
-}\
+// Bytecode helper functions ///////////////////////////////////////////////////
+void pushInstruction(Program &program, const Instruction instruction) {
+    program.push_back(static_cast<std::uint8_t>(instruction));
+}
+
+void pushAddress(Program &program, const unsigned int address) {
+    const auto address_ptr = reinterpret_cast<const std::uint8_t*>(&address);
+    const auto size = sizeof(unsigned int);
+    program.insert(program.end(), address_ptr, address_ptr + size);
+}
+
+void pushConstant(Program &program, const double value) {
+    const auto value_ptr = reinterpret_cast<const std::uint8_t*>(&value);
+    const auto size = sizeof(double);
+    program.insert(program.end(), value_ptr, value_ptr + size);
+}

 // VarNode compile /////////////////////////////////////////////////////////////
 void VarNode::compile(Program &program, RunContext &context) const
 {
-    PUSH_INS(program, Push, context.getVariableAddress(getName()), getName());
+    pushInstruction(program, Instruction::LOAD);
+    pushAddress(program, context.getVariableAddress(getName()));
 }

 // CstNode compile /////////////////////////////////////////////////////////////
 void CstNode::compile(Program &program, RunContext &context __dumb) const
 {
-    PUSH_INS(program, Push, strTo<double>(getName()));
+    pushInstruction(program, Instruction::CONST);
+    pushConstant(program, strTo<double>(getName()));
 }

 // SemicolonNode compile ///////////////////////////////////////////////////////
@@ -482,6 +367,10 @@ void SemicolonNode::compile(Program &program, RunContext &context) const
        {
            n[i].compile(program, context);
        }
+        // Where a variable has just been assigned, pop it off the stack.
+        if (isAssign) {
+            pushInstruction(program, Instruction::POP);
+        }
    }
 }

@@ -492,19 +381,10 @@ void AssignNode::compile(Program &program, RunContext &context) const
    
    if (isDerivedFrom<VarNode>(&n[0]))
    {
-        bool hasSemicolonParent = isDerivedFrom<SemicolonNode>(getParent());
-        unsigned int address;
-        
        n[1].compile(program, context);
-        address = context.addVariable(n[0].getName());
-        if (hasSemicolonParent)
-        {
-            PUSH_INS(program, Pop, address, n[0].getName());
-        }
-        else
-        {
-            PUSH_INS(program, Store, address, n[0].getName());
-        }
+        const unsigned int address = context.addVariable(n[0].getName());
+        pushInstruction(program, Instruction::STORE);
+        pushAddress(program, address);
    }
    else
    {
@@ -519,19 +399,37 @@ void AssignNode::compile(Program &program, RunContext &context) const

 void MathOpNode::compile(Program &program, RunContext &context) const
 {
+#define PUSH_BINARY_OP(op, instruction) \
+    case op: \
+        pushInstruction(program, Instruction::instruction); \
+        break;
+
    auto &n = *this;
    
    for (Index i = 0; i < n.getNArg(); ++i)
    {
        n[i].compile(program, context);
    }
-    IFNODE("-", 1)   PUSH_INS(program, Neg,);
-    ELIFNODE("+", 2) PUSH_INS(program, Add,);
-    ELIFNODE("-", 2) PUSH_INS(program, Sub,);
-    ELIFNODE("*", 2) PUSH_INS(program, Mul,);
-    ELIFNODE("/", 2) PUSH_INS(program, Div,);
-    ELIFNODE("^", 2) PUSH_INS(program, Pow,);
-    ELSE LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+    if (n.getName() == "-" and n.getNArg() == 1) {
+        pushInstruction(program, Instruction::NEG);
+        return;
+    }
+
+    if (getNArg() != 2) {
+        LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+    }
+
+    switch (getName()[0]) {
+    PUSH_BINARY_OP('+', ADD)
+    PUSH_BINARY_OP('-', SUB)
+    PUSH_BINARY_OP('*', MUL)
+    PUSH_BINARY_OP('/', DIV)
+    PUSH_BINARY_OP('^', POW)
+    default:
+        LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+    }
+
+#undef PUSH_BINARY_OP
 }

 // FuncNode compile ////////////////////////////////////////////////////////////
@@ -543,7 +441,8 @@ void FuncNode::compile(Program &program, RunContext &context) const
    {
        n[i].compile(program, context);
    }
-    PUSH_INS(program, Call, context.getFunctionAddress(getName()), getName());
+    pushInstruction(program, Instruction::CALL);
+    pushAddress(program, context.getFunctionAddress(getName()));
 }

 // ReturnNode compile ////////////////////////////////////////////////////////////
@@ -552,7 +451,7 @@ void ReturnNode::compile(Program &program, RunContext &context) const
    auto &n = *this;
    
    n[0].compile(program, context);
-    program.push_back(nullptr);
+    pushInstruction(program, Instruction::RET);
 }

 /******************************************************************************
@@ -582,7 +481,7 @@ MathInterpreter::MathInterpreter(const std::string &code)
 // access //////////////////////////////////////////////////////////////////////
 const Instruction * MathInterpreter::operator[](const Index i) const
 {
-    return program_[i].get();
+    return reinterpret_cast<const Instruction*>(&program_[i]);
 }

 const ExprNode * MathInterpreter::getAST(void) const
@@ -592,7 +491,7 @@ const ExprNode * MathInterpreter::getAST(void) const

 void MathInterpreter::push(const Instruction *i)
 {
-    program_.push_back(unique_ptr<const Instruction>(i));
+    pushInstruction(program_, *i);
 }

 // initialization //////////////////////////////////////////////////////////////
@@ -695,7 +594,7 @@ void MathInterpreter::compile(RunContext &context)
            root_->compile(program_, context);
            for (unsigned int i = 0; i < program_.size(); ++i)
            {
-                if (!program_[i])
+                if (static_cast<Instruction>(program_[i]) == Instruction::RET)
                {
                    gotReturn = true;
                    program_.resize(i);
@@ -726,20 +625,145 @@ void MathInterpreter::operator()(RunContext &context)

 void MathInterpreter::execute(RunContext &context) const
 {
-    context.setInsIndex(0);
-    while (context.getInsIndex() != program_.size())
-    {
-        (*(program_[context.getInsIndex()]))(context);
+#define BINARY_OP_CASE(instruction, expr) \
+        case Instruction::instruction: { \
+            const auto second = context.stack().top(); \
+            context.stack().pop(); \
+            const auto first = context.stack().top(); \
+            context.stack().pop(); \
+            context.stack().push(expr); \
+            break; \
+        }
+    auto ip = program_.begin();
+
+    while (ip != program_.end()) {
+        const auto instruction = static_cast<Instruction>(*ip);
+        ip++;
+
+        switch (instruction) {
+        BINARY_OP_CASE(ADD, first + second)
+        BINARY_OP_CASE(SUB, first - second)
+        BINARY_OP_CASE(MUL, first * second)
+        BINARY_OP_CASE(DIV, first / second)
+        BINARY_OP_CASE(POW, std::pow(first, second))
+        case Instruction::NEG: {
+            const auto operand = context.stack().top();
+            context.stack().pop();
+            context.stack().push(-operand);
+            break;
+        }
+        case Instruction::CONST: {
+            double value = 0.0;
+            std::tie(value, ip) = readConstant(ip);
+            context.stack().push(value);
+            break;
+        }
+        case Instruction::POP:
+            context.stack().pop();
+            break;
+        case Instruction::LOAD: {
+            unsigned int address = 0;
+            std::tie(address, ip) = readAddress(ip);
+            context.stack().push(context.getVariable(address));
+            break;
+        }
+        case Instruction::STORE: {
+            unsigned int address = 0;
+            std::tie(address, ip) = readAddress(ip);
+            context.setVariable(address, context.stack().top());
+            break;
+        }
+        case Instruction::CALL: {
+            unsigned int address = 0;
+            std::tie(address, ip) = readAddress(ip);
+            auto& stack = context.stack();
+            stack.push((*context.getFunction(address))(stack));
+            break;
+        }
+        case Instruction::RET:
+            break;
+        }
    }
+
+#undef BINARY_OP_CASE
+}
+
+Program::const_iterator instructionToStream(
+    ostream &out, Program::const_iterator ip)
+{
+    const auto instruction = static_cast<Instruction>(*ip);
+    ip++;
+
+    switch (instruction) {
+    case Instruction::ADD:
+        out << "ADD";
+        break;
+    case Instruction::SUB:
+        out << "SUB";
+        break;
+    case Instruction::MUL:
+        out << "MUL";
+        break;
+    case Instruction::DIV:
+        out << "DIV";
+        break;
+    case Instruction::POW:
+        out << "POW";
+        break;
+    case Instruction::NEG:
+        out << "NEG";
+        break;
+    case Instruction::CONST: {
+        double value = 0.0;
+        std::tie(value, ip) = readConstant(ip);
+        out << CODE_MOD << setfill(' ') << "CONST" << value;
+        break;
+    }
+    case Instruction::POP:
+        out << "POP";
+        break;
+    case Instruction::LOAD: {
+        unsigned int address = 0;
+        std::tie(address, ip) = readAddress(ip);
+        out << CODE_MOD << setfill(' ') << "LOAD" << address;
+        break;
+    }
+    case Instruction::STORE: {
+        unsigned int address = 0;
+        std::tie(address, ip) = readAddress(ip);
+        out << CODE_MOD << setfill(' ') << "STORE" << address;
+        break;
+    }
+    case Instruction::CALL: {
+        unsigned int address = 0;
+        std::tie(address, ip) = readAddress(ip);
+        out << CODE_MOD << setfill(' ') << "CALL" << address;
+        break;
+    }
+    case Instruction::RET:
+        out << "RET";
+        break;
+    }
+
+    return ip;
+}
+
+ostream &programToStream(ostream &out, const Program &program)
+{
+    auto ip = program.begin();
+
+    while (ip != program.end()) {
+        const auto i = std::distance(program.begin(), ip);
+        cout << setw(4) << setfill('0') << right << i << "  ";
+        ip = instructionToStream(out, ip);
+        out << '\n';
+    }
+
+    return out;
 }

 // IO //////////////////////////////////////////////////////////////////////////
 ostream &Latan::operator<<(ostream &out, const MathInterpreter &program)
 {
-    for (unsigned int i = 0; i < program.program_.size(); ++i)
-    {
-        out << *(program.program_[i]) << endl;
-    }
-    
-    return out;
+    return programToStream(out, program.program_);
 }
--- a/lib/Core/MathInterpreter.hpp
+++ b/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
@@ -20,6 +20,8 @@
 #ifndef Latan_MathInterpreter_hpp_
 #define	Latan_MathInterpreter_hpp_

+#include <cstdint>
+
 #include <LatAnalyze/Functional/Function.hpp>
 #include <LatAnalyze/Global.hpp>
 #include <LatAnalyze/Core/ParserState.hpp>
@@ -79,108 +81,26 @@ private:
 *                         Instruction classes                                *
 ******************************************************************************/
 // Abstract base
-class Instruction
-{
-public:
-    // destructor
-    virtual ~Instruction(void) = default;
-    // instruction execution
-    virtual void operator()(RunContext &context) const = 0;
-    friend std::ostream & operator<<(std::ostream &out, const Instruction &ins);
-private:
-    virtual void print(std::ostream &out) const = 0;
+
+enum class Instruction : std::uint8_t {
+    ADD,
+    SUB,
+    MUL,
+    DIV,
+    POW,
+    NEG,
+    CONST,
+    POP,
+    LOAD,
+    STORE,
+    CALL,
+    RET,
 };

 std::ostream & operator<<(std::ostream &out, const Instruction &ins);

 // Instruction container
-typedef std::vector<std::unique_ptr<const Instruction>> Program;
-
-// Push
-class Push: public Instruction
-{
-private:
-    enum class ArgType
-    {
-        Constant = 0,
-        Variable = 1
-    };
-public:
-    //constructors
-    explicit Push(const double val);
-    explicit Push(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    ArgType      type_;
-    double       val_;
-    unsigned int address_;
-    std::string  name_;
-};
-
-// Pop
-class Pop: public Instruction
-{
-public:
-    //constructor
-    explicit Pop(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    unsigned int address_;
-    std::string name_;
-};
-
-// Store
-class Store: public Instruction
-{
-public:
-    //constructor
-    explicit Store(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    unsigned int address_;
-    std::string name_;
-};
-
-// Call function
-class Call: public Instruction
-{
-public:
-    //constructor
-    explicit Call(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    unsigned int address_;
-    std::string name_;
-};
-
-// Floating point operations
-#define DECL_OP(name)\
-class name: public Instruction\
-{\
-public:\
-    virtual void operator()(RunContext &context) const;\
-private:\
-    virtual void print(std::ostream &out) const;\
-}
-
-DECL_OP(Neg);
-DECL_OP(Add);
-DECL_OP(Sub);
-DECL_OP(Mul);
-DECL_OP(Div);
-DECL_OP(Pow);
+typedef std::vector<std::uint8_t> Program;

 /******************************************************************************
 *                            Expression node classes                         *
--- a/lib/Core/MathLexer.lpp
+++ b/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
--- a/lib/Core/MathParser.ypp
+++ b/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
--- a/lib/Core/ParserState.hpp
+++ b/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
--- a/lib/Core/Plot.cpp
+++ b/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
@@ -112,7 +112,7 @@ PlotHeadCommand::PlotHeadCommand(const string &command)
 }

 // PlotData constructor ////////////////////////////////////////////////////////
-PlotData::PlotData(const DMatSample &x, const DMatSample &y)
+PlotData::PlotData(const DMatSample &x, const DMatSample &y, const bool abs)
 {
    if (x[central].rows() != y[central].rows())
    {
@@ -122,16 +122,23 @@ PlotData::PlotData(const DMatSample &x, const DMatSample &y)
    DMat d(x[central].rows(), 4);
    string usingCmd, tmpFileName;

-    d.col(0)    = x[central];
-    d.col(2)    = y[central];
-    d.col(1)    = x.variance().cwiseSqrt();
-    d.col(3)    = y.variance().cwiseSqrt();
+    d.col(0)    = x[central].col(0);
+    d.col(2)    = y[central].col(0);
+    d.col(1)    = x.variance().cwiseSqrt().col(0);
+    d.col(3)    = y.variance().cwiseSqrt().col(0);
    tmpFileName = dumpToTmpFile(d);
    pushTmpFile(tmpFileName);
-    setCommand("'" + tmpFileName + "' u 1:3:2:4 w xyerr");
+    if (!abs)
+    {
+        setCommand("'" + tmpFileName + "' u 1:3:2:4 w xyerr");
+    }
+    else
+    {
+        setCommand("'" + tmpFileName + "' u 1:(abs($3)):2:4 w xyerr");
+    }
 }

-PlotData::PlotData(const DVec &x, const DMatSample &y)
+PlotData::PlotData(const DVec &x, const DMatSample &y, const bool abs)
 {
    if (x.rows() != y[central].rows())
    {
@@ -142,14 +149,21 @@ PlotData::PlotData(const DVec &x, const DMatSample &y)
    string usingCmd, tmpFileName;

    d.col(0)    = x;
-    d.col(1)    = y[central];
-    d.col(2)    = y.variance().cwiseSqrt();
+    d.col(1)    = y[central].col(0);
+    d.col(2)    = y.variance().cwiseSqrt().col(0);
    tmpFileName = dumpToTmpFile(d);
    pushTmpFile(tmpFileName);
-    setCommand("'" + tmpFileName + "' u 1:2:3 w yerr");
+    if (!abs)
+    {
+        setCommand("'" + tmpFileName + "' u 1:2:3 w yerr");
+    }
+    else
+    {
+        setCommand("'" + tmpFileName + "' u 1:(abs($2)):3 w yerr");
+    }
 }

-PlotData::PlotData(const DMatSample &x, const DVec &y)
+PlotData::PlotData(const DMatSample &x, const DVec &y, const bool abs)
 {
    if (x[central].rows() != y.rows())
    {
@@ -159,24 +173,93 @@ PlotData::PlotData(const DMatSample &x, const DVec &y)
    DMat d(x[central].rows(), 3), xerr, yerr;
    string usingCmd, tmpFileName;

-    d.col(0)    = x[central];
+    d.col(0)    = x[central].col(0);
    d.col(2)    = y;
-    d.col(1)    = x.variance().cwiseSqrt();
+    d.col(1)    = x.variance().cwiseSqrt().col(0);
+    tmpFileName = dumpToTmpFile(d);
+    pushTmpFile(tmpFileName);
+    if (!abs)
+    {
+        setCommand("'" + tmpFileName + "' u 1:3:2 w xerr");
+    }
+    else
+    {
+        setCommand("'" + tmpFileName + "' u 1:(abs($3)):2 w xerr");
+    }
+}
+
+PlotData::PlotData(const XYStatData &data, const Index i, const Index j, const bool abs)
+{
+    string usingCmd, tmpFileName;
+    
+    if (!abs)
+    {
+        usingCmd = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:3:2:4 w xyerr";
+    }
+    else
+    {
+        usingCmd = (data.isXExact(i)) ? "u 1:(abs($3)):4 w yerr" : "u 1:(abs($3)):2:4 w xyerr";
+    }
+   
+    tmpFileName = dumpToTmpFile(data.getTable(i, j));
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' " + usingCmd);
+}
+
+// PlotPoint constructor ///////////////////////////////////////////////////////
+PlotPoint::PlotPoint(const double x, const double y)
+{
+    DMat d(1, 2);
+    string usingCmd, tmpFileName;
+
+    d(0, 0)     = x;
+    d(0, 1)     = y;
+    tmpFileName = dumpToTmpFile(d);
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' u 1:2");
+}
+
+PlotPoint::PlotPoint(const DSample &x, const double y)
+{
+    DMat d(1, 3);
+    string usingCmd, tmpFileName;
+
+    d(0, 0)     = x[central];
+    d(0, 2)     = y;
+    d(0, 1)     = sqrt(x.variance());
    tmpFileName = dumpToTmpFile(d);
    pushTmpFile(tmpFileName);
    setCommand("'" + tmpFileName + "' u 1:3:2 w xerr");
 }

-PlotData::PlotData(const XYStatData &data, const Index i, const Index j)
+PlotPoint::PlotPoint(const double x, const DSample &y)
 {
+    DMat d(1, 3);
    string usingCmd, tmpFileName;
-    
-    usingCmd    = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:3:2:4 w xyerr";
-    tmpFileName = dumpToTmpFile(data.getTable(i, j));
+
+    d(0, 0)     = x;
+    d(0, 1)     = y[central];
+    d(0, 2)     = sqrt(y.variance());
+    tmpFileName = dumpToTmpFile(d);
    pushTmpFile(tmpFileName);
-    setCommand("'" + tmpFileName + "' " + usingCmd);
+    setCommand("'" + tmpFileName + "' u 1:2:3 w yerr");
 }

+PlotPoint::PlotPoint(const DSample &x, const DSample &y)
+{
+    DMat d(1, 4);
+    string usingCmd, tmpFileName;
+
+    d(0, 0)     = x[central];
+    d(0, 2)     = y[central];
+    d(0, 1)     = sqrt(x.variance());
+    d(0, 3)     = sqrt(y.variance());
+    tmpFileName = dumpToTmpFile(d);
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' u 1:3:2:4 w xyerr");
+}
+
+
 // PlotLine constructor ////////////////////////////////////////////////////////
 PlotLine::PlotLine(const DVec &x, const DVec &y)
 {
@@ -195,6 +278,24 @@ PlotLine::PlotLine(const DVec &x, const DVec &y)
    setCommand("'" + tmpFileName + "' u 1:2 w lines");
 }

+// PlotPoints constructor ////////////////////////////////////////////////////////
+PlotPoints::PlotPoints(const DVec &x, const DVec &y)
+{
+    if (x.size() != y.size())
+    {
+        LATAN_ERROR(Size, "x and y vectors do not have the same size");
+    }
+
+    DMat d(x.size(), 2);
+    string usingCmd, tmpFileName;
+
+    d.col(0)    = x;
+    d.col(1)    = y;
+    tmpFileName = dumpToTmpFile(d);
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' u 1:2");
+}
+
 // PlotHLine constructor ///////////////////////////////////////////////////////
 PlotHLine::PlotHLine(const double y)
 {
@@ -217,7 +318,8 @@ PlotBand::PlotBand(const double xMin, const double xMax, const double yMin,

 // PlotFunction constructor ////////////////////////////////////////////////////
 PlotFunction::PlotFunction(const DoubleFunction &function, const double xMin,
-                           const double xMax, const unsigned int nPoint)
+                           const double xMax, const unsigned int nPoint,
+                           const bool abs)
 {
    DMat   d(nPoint, 2);
    string tmpFileName;
@@ -230,20 +332,34 @@ PlotFunction::PlotFunction(const DoubleFunction &function, const double xMin,
    }
    tmpFileName = dumpToTmpFile(d);
    pushTmpFile(tmpFileName);
-    setCommand("'" + tmpFileName + "' u 1:2 w lines");
+    if (!abs)
+    {
+        setCommand("'" + tmpFileName + "' u 1:2 w lines");
+    }
+    else
+    {
+        setCommand("'" + tmpFileName + "' u 1:(abs($2)) w lines");
+    }
 }

 // PlotPredBand constructor ////////////////////////////////////////////////////
 void PlotPredBand::makePredBand(const DMat &low, const DMat &high, const double opacity)
 {
-    string lowFileName, highFileName;
+    string lowFileName, highFileName, contFileName;
+    DMat   contour(low.rows() + high.rows() + 1, 2);

-    lowFileName  = dumpToTmpFile(low);
-    highFileName = dumpToTmpFile(high);
-    pushTmpFile(lowFileName);
-    pushTmpFile(highFileName);
-    setCommand("'< (cat " + lowFileName + "; tac " + highFileName +
-               "; head -n1 " + lowFileName + ")' u 1:2 w filledcurves closed" +
+    FOR_MAT(low, i, j)
+    {
+        contour(i, j) = low(i, j);
+    }
+    FOR_MAT(high, i, j)
+    {
+        contour(low.rows() + i, j) = high(high.rows() - i - 1, j);
+    }
+    contour.row(low.rows() + high.rows()) = low.row(0);
+    contFileName = dumpToTmpFile(contour);
+    pushTmpFile(contFileName);
+    setCommand("'" + contFileName + "' u 1:2 w filledcurves closed" +
               " fs solid " + strFrom(opacity) + " noborder");
 }

@@ -500,7 +616,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 +696,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;
-        path = getenv("PATH") ;
-        if (path)
+        for (i=0;path[i];)
        {
-            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++);
-                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++] = '.';
            }
+            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
-        {
-            LATAN_ERROR(System, "PATH variable not set");
-        }
-        // 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
+    }
+    // if the buffer is still empty, the command was not found
+    if (buf[0] != 0)
+    {
        lg = (int)(strlen(buf) - 1);
        while (buf[lg]!='/')
        {
@@ -639,7 +746,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 +772,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 +798,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 +806,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 +813,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";
-    options_.output   = dirName + "/plot.pdf";
-    display();
-    options_.terminal = terminalBack;
-    options_.output   = outputBack;
+    if (savePdf)
+    {
+        options_.terminal = "pdf";
+        options_.output   = dirName + "/plot.pdf";
+        display();
+        options_.terminal = terminalBack;
+        options_.output   = outputBack;
+    }
    
    // save script and datafiles
    for (unsigned int i = 0; i < tmpFileName_.size(); ++i)
@@ -776,18 +894,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 << "set title '" << plot.options_.title << "'" << endl;
-    }
+    out << "unset xrange" << 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;
--- a/lib/Core/Plot.hpp
+++ b/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
@@ -89,14 +89,27 @@ class PlotData: public PlotObject
 {
 public:
    // constructor
-    PlotData(const DMatSample &x, const DMatSample &y);
-    PlotData(const DVec       &x, const DMatSample &y);
-    PlotData(const DMatSample &x, const DVec       &y);
-    PlotData(const XYStatData &data, const Index i = 0, const Index j = 0);
+    PlotData(const DMatSample &x, const DMatSample &y, const bool abs = false);
+    PlotData(const DVec       &x, const DMatSample &y, const bool abs = false);
+    PlotData(const DMatSample &x, const DVec       &y, const bool abs = false);
+    PlotData(const XYStatData &data, const Index i = 0, const Index j = 0, 
+             const bool abs = false);
    // destructor
    virtual ~PlotData(void) = default;
 };

+class PlotPoint: public PlotObject
+{
+public:
+    // constructor
+    PlotPoint(const double  x,  const double  y);
+    PlotPoint(const DSample &x, const double  y);
+    PlotPoint(const double  x,  const DSample &y);
+    PlotPoint(const DSample &x, const DSample &y);
+    // destructor
+    virtual ~PlotPoint(void) = default;
+};
+
 class PlotHLine: public PlotObject
 {
 public:
@@ -115,6 +128,15 @@ public:
    virtual ~PlotLine(void) = default;
 };

+class PlotPoints: public PlotObject
+{
+public:
+    // constructor
+    PlotPoints(const DVec &x, const DVec &y);
+    // destructor
+    virtual ~PlotPoints(void) = default;
+};
+
 class PlotBand: public PlotObject
 {
 public:
@@ -130,7 +152,8 @@ class PlotFunction: public PlotObject
 public:
    // constructor
    PlotFunction(const DoubleFunction &function, const double xMin,
-                 const double xMax, const unsigned int nPoint = 1000);
+                 const double xMax, const unsigned int nPoint = 1000, 
+                 const bool abs = false);
    // destructor
    virtual ~PlotFunction(void) = default;
 };
@@ -182,6 +205,11 @@ PlotRange(Axis::x, -.5, (m).cols() - .5) <<\
 PlotRange(Axis::y, (m).rows() - .5, -.5) <<\
 PlotMatrixNoRange(m)

+#define PlotCorrMatrix(m)\
+PlotHeadCommand("set cbrange [-1:1]") <<\
+PlotHeadCommand("set palette defined (0 'blue', 1 'white', 2 'red')") <<\
+PlotMatrix(m)
+
 /******************************************************************************
 *                             Plot modifiers                                 *
 ******************************************************************************/
@@ -377,13 +405,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:
--- a/lib/Core/ThreadPool.cpp
+++ b/lib/Core/ThreadPool.cpp
@@ -0,0 +1,117 @@
+/*
+ * ThreadPool.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2021 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/Core/ThreadPool.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                         ThreadPool implementation                          *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+ThreadPool::ThreadPool(void)
+: ThreadPool(std::thread::hardware_concurrency())
+{}
+
+ThreadPool::ThreadPool(const unsigned int nThreads)
+: nThreads_(nThreads)
+{
+    for (unsigned int t = 0; t < nThreads_; ++t)
+    {
+        threads_.push_back(thread(&ThreadPool::workerLoop, this));
+    }
+}
+
+// destructor //////////////////////////////////////////////////////////////////
+ThreadPool::~ThreadPool(void)
+{
+    terminate();
+}
+
+// get the number of threads ///////////////////////////////////////////////////
+unsigned int ThreadPool::getThreadNum(void) const
+{
+    return nThreads_;
+}
+
+// get the pool mutex for synchronisation //////////////////////////////////////
+std::mutex & ThreadPool::getMutex(void)
+{
+    return mutex_;
+}
+
+// worker loop /////////////////////////////////////////////////////////////////
+void ThreadPool::workerLoop(void)
+{
+    while (true)
+    {
+        Job job;
+        {
+            unique_lock<mutex> lock(mutex_);
+
+            condition_.wait(lock, [this](){
+                return !queue_.empty() || terminatePool_;
+            });
+            if (terminatePool_ and queue_.empty())
+            {
+                return;
+            }
+            job = queue_.front();
+            queue_.pop();
+        }
+        job();
+    }
+}
+
+// add jobs ////////////////////////////////////////////////////////////////////
+void ThreadPool::addJob(Job newJob)
+{
+    {
+        unique_lock<mutex> lock(mutex_);
+
+        queue_.push(newJob);
+    }
+    condition_.notify_one();
+}
+
+// critical section ////////////////////////////////////////////////////////////
+void ThreadPool::critical(Job fn)
+{
+    unique_lock<mutex> lock(mutex_);
+
+    fn();
+}
+
+// wait for completion /////////////////////////////////////////////////////////
+void ThreadPool::terminate(void)
+{
+    {
+        unique_lock<mutex> lock(mutex_);
+
+        terminatePool_ = true;
+    }
+    condition_.notify_all();
+    for (auto &thread: threads_)
+    {
+        thread.join();
+    }
+    threads_.clear();
+}
--- a/lib/Core/ThreadPool.hpp
+++ b/lib/Core/ThreadPool.hpp
@@ -0,0 +1,56 @@
+/*
+ * ThreadPool.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2021 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_ThreadPool_hpp_
+#define Latan_ThreadPool_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+class ThreadPool
+{
+public:
+    typedef std::function<void(void)> Job;
+public:
+    // constructors/destructor
+    ThreadPool(void);
+    ThreadPool(const unsigned int nThreads);
+    virtual ~ThreadPool(void);
+    // get the number of threads
+    unsigned int getThreadNum(void) const;
+    // get the pool mutex for synchronisation
+    std::mutex & getMutex(void);
+    // add jobs
+    void addJob(Job newJob);
+    // critical section
+    void critical(Job fn);
+    // wait for completion and terminate
+    void terminate(void);
+private:
+    // worker loop
+    void workerLoop(void);
+private:
+    unsigned int             nThreads_;
+    std::condition_variable  condition_;
+    std::vector<std::thread> threads_;
+    bool                     terminatePool_{false};
+    std::queue<Job>          queue_;
+    std::mutex               mutex_;
+};
+
+#endif
--- a/lib/Core/Utilities.cpp
+++ b/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
--- a/lib/Core/Utilities.hpp
+++ b/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
--- a/lib/Core/stdincludes.hpp
+++ b/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
@@ -24,6 +24,7 @@
 #include <array>
 #include <chrono>
 #include <complex>
+#include <condition_variable>
 #include <fstream>
 #include <functional>
 #include <iostream>
@@ -40,6 +41,7 @@
 #include <stack>
 #include <string>
 #include <sstream>
+#include <thread>
 #include <type_traits>
 #include <unordered_map>
 #include <utility>
--- a/lib/Functional/CompiledFunction.cpp
+++ b/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
--- a/lib/Functional/CompiledFunction.hpp
+++ b/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
--- a/lib/Functional/CompiledModel.cpp
+++ b/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
--- a/lib/Functional/CompiledModel.hpp
+++ b/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
--- a/lib/Functional/Function.cpp
+++ b/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
--- a/lib/Functional/Function.hpp
+++ b/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
--- a/lib/Functional/Model.cpp
+++ b/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
--- a/lib/Functional/Model.hpp
+++ b/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
--- a/lib/Functional/TabFunction.cpp
+++ b/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
--- a/lib/Functional/TabFunction.hpp
+++ b/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
--- a/lib/Global.cpp
+++ b/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
--- a/lib/Global.hpp
+++ b/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
--- a/lib/Io/AsciiFile.cpp
+++ b/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
--- a/lib/Io/AsciiFile.hpp
+++ b/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
--- a/lib/Io/AsciiLexer.lpp
+++ b/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
--- a/lib/Io/AsciiParser.ypp
+++ b/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
--- a/lib/Io/File.cpp
+++ b/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
--- a/lib/Io/File.hpp
+++ b/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
--- a/lib/Io/Hdf5File.cpp
+++ b/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
--- a/lib/Io/Hdf5File.hpp
+++ b/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
--- a/lib/Io/Io.cpp
+++ b/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
--- a/lib/Io/Io.hpp
+++ b/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
--- a/lib/Io/IoObject.hpp
+++ b/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
--- a/lib/Io/XmlReader.cpp
+++ b/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
--- a/lib/Io/XmlReader.hpp
+++ b/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
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -32,6 +32,7 @@ libLatAnalyze_la_SOURCES =           \
    Core/MathParser.ypp              \
    Core/OptParser.cpp               \
    Core/Plot.cpp                    \
+    Core/ThreadPool.cpp              \
    Core/Utilities.cpp               \
    Functional/CompiledFunction.cpp  \
    Functional/CompiledModel.cpp     \
@@ -47,6 +48,8 @@ libLatAnalyze_la_SOURCES =           \
    Io/XmlReader.cpp                 \
    Io/Xml/tinyxml2.cpp              \
    Numerical/Derivative.cpp         \
+    Numerical/DWT.cpp                \
+    Numerical/DWTFilters.cpp         \
    Numerical/GslFFT.cpp             \
    Numerical/GslHybridRootFinder.cpp\
    Numerical/GslMinimizer.cpp       \
@@ -54,6 +57,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            \
@@ -73,6 +78,7 @@ HPPFILES =                           \
    Core/OptParser.hpp               \
    Core/ParserState.hpp             \
    Core/Plot.hpp                    \
+    Core/ThreadPool.hpp              \
    Core/stdincludes.hpp             \
    Core/Utilities.hpp               \
    Functional/CompiledFunction.hpp  \
@@ -88,6 +94,8 @@ HPPFILES =                           \
    Io/IoObject.hpp                  \
    Io/XmlReader.hpp                 \
    Numerical/Derivative.hpp         \
+    Numerical/DWT.hpp                \
+    Numerical/DWTFilters.hpp         \
    Numerical/FFT.hpp                \
    Numerical/GslFFT.hpp             \
    Numerical/GslHybridRootFinder.hpp\
@@ -97,6 +105,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         \
--- a/lib/Numerical/DWT.cpp
+++ b/lib/Numerical/DWT.cpp
@@ -0,0 +1,137 @@
+/*
+ * DWT.cpp, part of LatAnalyze
+ *
+ * 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
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 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.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Numerical/DWT.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                            DWT implementation                              *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+DWT::DWT(const DWTFilter &filter)
+: filter_(filter)
+{}
+
+// convolution primitive ///////////////////////////////////////////////////////
+void DWT::filterConvolution(DVec &out, const DVec &data, 
+                            const std::vector<double> &filter, const Index offset)
+{
+    Index n = data.size(), nf = n*filter.size();
+
+    out.resize(n);
+    out.fill(0.);
+    for (unsigned int i = 0; i < filter.size(); ++i)
+    {
+        FOR_VEC(out, j)
+        {
+            out(j) += filter[i]*data((j + i + nf - offset) % n);
+        }
+    }
+}
+
+// downsampling/upsampling primitives //////////////////////////////////////////
+void DWT::downsample(DVec &out, const DVec &in)
+{
+    if (out.size() < in.size()/2)
+    {
+        LATAN_ERROR(Size, "output vector smaller than half the input vector size");
+    }
+    for (Index i = 0; i < in.size(); i += 2)
+    {
+        out(i/2) = in(i);
+    }
+}
+
+void DWT::upsample(DVec &out, const DVec &in)
+{
+    if (out.size() < 2*in.size())
+    {
+        LATAN_ERROR(Size, "output vector smaller than twice the input vector size");
+    }
+    out.segment(0, 2*in.size()).fill(0.);
+    for (Index i = 0; i < in.size(); i ++)
+    {
+        out(2*i) = in(i);
+    }
+}
+
+// DWT /////////////////////////////////////////////////////////////////////////
+std::vector<DWT::DWTLevel> 
+DWT::forward(const DVec &data, const unsigned int level) const
+{
+    std::vector<DWTLevel> dwt(level);
+    DVec                  *finePt = const_cast<DVec *>(&data);
+    DVec                  tmp; 
+    Index                 n = data.size(), o = filter_.fwdL.size()/2, minSize;
+
+    minSize = 1;
+    for (unsigned int l = 0; l < level; ++l) minSize *= 2;
+    if (n < minSize)
+    {
+        LATAN_ERROR(Size, "data vector too small for a " + strFrom(level) 
+                          + "-level DWT (data size is " + strFrom(n) + ")");
+    }
+    for (unsigned int l = 0; l < level; ++l)
+    {
+        n /= 2;
+        dwt[l].first.resize(n);
+        dwt[l].second.resize(n);
+        filterConvolution(tmp, *finePt, filter_.fwdL, o);
+        downsample(dwt[l].first, tmp);
+        filterConvolution(tmp, *finePt, filter_.fwdH, o);
+        downsample(dwt[l].second, tmp);
+        finePt = &dwt[l].first;
+    }
+
+    return dwt;
+}
+
+DVec DWT::backward(const std::vector<DWTLevel>& dwt) const
+{
+    unsigned int level = dwt.size();
+    Index        n = dwt.back().second.size(), o = filter_.bwdL.size()/2 - 1;
+    DVec         res, tmp, conv;
+
+    res = dwt.back().first;
+    for (int l = level - 2; l >= 0; --l)
+    {
+        n *= 2;
+        if (dwt[l].second.size() != n)
+        {
+            LATAN_ERROR(Size, "DWT result size mismatch");
+        }
+    }
+    n = dwt.back().second.size();
+    for (int l = level - 1; l >= 0; --l)
+    {
+        n *= 2;
+        tmp.resize(n);
+        upsample(tmp, res);
+        filterConvolution(conv, tmp, filter_.bwdL, o);
+        res = conv;
+        upsample(tmp, dwt[l].second);
+        filterConvolution(conv, tmp, filter_.bwdH, o);
+        res += conv;
+    }
+
+    return res;
+}
--- a/lib/Numerical/DWT.hpp
+++ b/lib/Numerical/DWT.hpp
@@ -0,0 +1,55 @@
+/*
+ * DWT.hpp, part of LatAnalyze
+ *
+ * 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
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 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.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_DWT_hpp_
+#define Latan_DWT_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Numerical/DWTFilters.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                   Discrete wavelet transform class                         *
+ ******************************************************************************/
+class DWT
+{
+public:
+  typedef std::pair<DVec, DVec> DWTLevel;
+public:
+    // constructor
+    DWT(const DWTFilter &filter);
+    // destructor
+    virtual ~DWT(void) = default;
+    // convolution primitive
+    static void filterConvolution(DVec &out, const DVec &data, 
+                                  const std::vector<double> &filter, const Index offset);
+    // downsampling/upsampling primitives
+    static void downsample(DVec &out, const DVec &in);
+    static void upsample(DVec &out, const DVec &in);
+    // DWT
+    std::vector<DWTLevel> forward(const DVec &data, const unsigned int level) const;
+    DVec                  backward(const std::vector<DWTLevel>& dwt) const;
+private:
+    DWTFilter filter_;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_DWT_hpp_
--- a/lib/Numerical/DWTFilters.cpp
+++ b/lib/Numerical/DWTFilters.cpp
@@ -0,0 +1,528 @@
+/*
+ * DWTFilters.cpp, part of LatAnalyze
+ *
+ * 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
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 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.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Numerical/DWTFilters.hpp>
+#include <LatAnalyze/includes.hpp>
+
+// cf. http://wavelets.pybytes.com
+// *here we implement the reverse filters more convenient for convolutions*
+
+using namespace std;
+using namespace Latan;
+
+#define FILTDICT(x) {#x, &DWTFilters::x}
+
+std::map<std::string, const DWTFilter *> DWTFilters::fromName = {
+    FILTDICT(haar),
+    FILTDICT(db2),
+    FILTDICT(db3),
+    FILTDICT(db4),
+    FILTDICT(db5),
+    FILTDICT(db6),
+    FILTDICT(bior13),
+    FILTDICT(bior15),
+    FILTDICT(bior22),
+    FILTDICT(bior24),
+    FILTDICT(bior31),
+    FILTDICT(bior33),
+    FILTDICT(bior35)
+};
+
+DWTFilter DWTFilters::haar = {
+    // fwdL
+    {0.7071067811865476,
+    0.7071067811865476},
+    // fwdH
+    {0.7071067811865476,
+    -0.7071067811865476},
+    // bwdL
+    {0.7071067811865476,
+    0.7071067811865476},
+    // bwdH
+    {-0.7071067811865476,
+    0.7071067811865476}
+};
+
+DWTFilter DWTFilters::db2 = {
+    // fwdL
+    {0.48296291314469025,
+    0.836516303737469,
+    0.22414386804185735,
+    -0.12940952255092145},
+    // fwdH
+    {-0.12940952255092145,
+    -0.22414386804185735,
+    0.836516303737469,
+    -0.48296291314469025},
+    // bwdL
+    {-0.12940952255092145,
+    0.22414386804185735,
+    0.836516303737469,
+    0.48296291314469025},
+    // bwdH
+    {-0.48296291314469025,
+    0.836516303737469,
+    -0.22414386804185735,
+    -0.12940952255092145}
+};
+
+DWTFilter DWTFilters::db3 = {
+    // fwdL
+    {0.3326705529509569,
+    0.8068915093133388,
+    0.4598775021193313,
+    -0.13501102001039084,
+    -0.08544127388224149,
+    0.035226291882100656},
+    // fwdH
+    {0.035226291882100656,
+    0.08544127388224149,
+    -0.13501102001039084,
+    -0.4598775021193313,
+    0.8068915093133388,
+    -0.3326705529509569},
+    // bwdL
+    {0.035226291882100656,
+    -0.08544127388224149,
+    -0.13501102001039084,
+    0.4598775021193313,
+    0.8068915093133388,
+    0.3326705529509569},
+    // bwdH
+    {-0.3326705529509569,
+    0.8068915093133388,
+    -0.4598775021193313,
+    -0.13501102001039084,
+    0.08544127388224149,
+    0.035226291882100656}
+};
+
+DWTFilter DWTFilters::db4 = {
+    // fwdL
+    {0.23037781330885523,
+    0.7148465705525415,
+    0.6308807679295904,
+    -0.02798376941698385,
+    -0.18703481171888114,
+    0.030841381835986965,
+    0.032883011666982945,
+    -0.010597401784997278},
+    // fwdH
+    {-0.010597401784997278,
+    -0.032883011666982945,
+    0.030841381835986965,
+    0.18703481171888114,
+    -0.02798376941698385,
+    -0.6308807679295904,
+    0.7148465705525415,
+    -0.23037781330885523},
+    // bwdL
+    {-0.010597401784997278,
+    0.032883011666982945,
+    0.030841381835986965,
+    -0.18703481171888114,
+    -0.02798376941698385,
+    0.6308807679295904,
+    0.7148465705525415,
+    0.23037781330885523},
+    // bwdH
+    {-0.23037781330885523,
+    0.7148465705525415,
+    -0.6308807679295904,
+    -0.02798376941698385,
+    0.18703481171888114,
+    0.030841381835986965,
+    -0.032883011666982945,
+    -0.010597401784997278}
+};
+
+DWTFilter DWTFilters::db5 = {
+    // fwdL
+    {0.160102397974125,
+    0.6038292697974729,
+    0.7243085284385744,
+    0.13842814590110342,
+    -0.24229488706619015,
+    -0.03224486958502952,
+    0.07757149384006515,
+    -0.006241490213011705,
+    -0.012580751999015526,
+    0.003335725285001549},
+    // fwdH
+    {0.003335725285001549,
+    0.012580751999015526,
+    -0.006241490213011705,
+    -0.07757149384006515,
+    -0.03224486958502952,
+    0.24229488706619015,
+    0.13842814590110342,
+    -0.7243085284385744,
+    0.6038292697974729,
+    -0.160102397974125},
+    // bwdL
+    {0.003335725285001549,
+    -0.012580751999015526,
+    -0.006241490213011705,
+    0.07757149384006515,
+    -0.03224486958502952,
+    -0.24229488706619015,
+    0.13842814590110342,
+    0.7243085284385744,
+    0.6038292697974729,
+    0.160102397974125},
+    // bwdH
+    {-0.160102397974125,
+    0.6038292697974729,
+    -0.7243085284385744,
+    0.13842814590110342,
+    0.24229488706619015,
+    -0.03224486958502952,
+    -0.07757149384006515,
+    -0.006241490213011705,
+    0.012580751999015526,
+    0.003335725285001549}
+};
+
+DWTFilter DWTFilters::db6 = {
+    // fwdL
+    {0.11154074335008017,
+    0.4946238903983854,
+    0.7511339080215775,
+    0.3152503517092432,
+    -0.22626469396516913,
+    -0.12976686756709563,
+    0.09750160558707936,
+    0.02752286553001629,
+    -0.031582039318031156,
+    0.0005538422009938016,
+    0.004777257511010651,
+    -0.00107730108499558},
+    // fwdH
+    {-0.00107730108499558,
+    -0.004777257511010651,
+    0.0005538422009938016,
+    0.031582039318031156,
+    0.02752286553001629,
+    -0.09750160558707936,
+    -0.12976686756709563,
+    0.22626469396516913,
+    0.3152503517092432,
+    -0.7511339080215775,
+    0.4946238903983854,
+    -0.11154074335008017},
+    // bwdL
+    {-0.00107730108499558,
+    0.004777257511010651,
+    0.0005538422009938016,
+    -0.031582039318031156,
+    0.02752286553001629,
+    0.09750160558707936,
+    -0.12976686756709563,
+    -0.22626469396516913,
+    0.3152503517092432,
+    0.7511339080215775,
+    0.4946238903983854,
+    0.11154074335008017},
+    // bwdH
+    {-0.11154074335008017,
+    0.4946238903983854,
+    -0.7511339080215775,
+    0.3152503517092432,
+    0.22626469396516913,
+    -0.12976686756709563,
+    -0.09750160558707936,
+    0.02752286553001629,
+    0.031582039318031156,
+    0.0005538422009938016,
+    -0.004777257511010651,
+    -0.00107730108499558}
+};
+
+DWTFilter DWTFilters::bior13 = {
+    // fwdL
+    {-0.08838834764831845,
+    0.08838834764831845,
+    0.7071067811865476,
+    0.7071067811865476,
+    0.08838834764831845,
+    -0.08838834764831845},
+    // fwdH
+    {0.0,
+    0.0,
+    0.7071067811865476,
+    -0.7071067811865476,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.7071067811865476,
+    0.7071067811865476,
+    0.0,
+    0.0},
+    // bwdH
+    {0.08838834764831845,
+    0.08838834764831845,
+    -0.7071067811865476,
+    0.7071067811865476,
+    -0.08838834764831845,
+    -0.08838834764831845}
+};
+
+DWTFilter DWTFilters::bior15 = {
+    // fwdL
+    {0.01657281518405971,
+    -0.01657281518405971,
+    -0.12153397801643787,
+    0.12153397801643787,
+    0.7071067811865476,
+    0.7071067811865476,
+    0.12153397801643787,
+    -0.12153397801643787,
+    -0.01657281518405971,
+    0.01657281518405971},
+    // fwdH
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.7071067811865476,
+    -0.7071067811865476,
+    0.0,
+    0.0,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.7071067811865476,
+    0.7071067811865476,
+    0.0,
+    0.0,
+    0.0,
+    0.0},
+    // bwdH
+    {-0.01657281518405971,
+    -0.01657281518405971,
+    0.12153397801643787,
+    0.12153397801643787,
+    -0.7071067811865476,
+    0.7071067811865476,
+    -0.12153397801643787,
+    -0.12153397801643787,
+    0.01657281518405971,
+    0.01657281518405971}
+};
+
+DWTFilter DWTFilters::bior22 = {
+    // fwdL
+    {-0.1767766952966369,
+    0.3535533905932738,
+    1.0606601717798214,
+    0.3535533905932738,
+    -0.1767766952966369,
+    0.0},
+    // fwdH
+    {0.0,
+    0.0,
+    0.3535533905932738,
+    -0.7071067811865476,
+    0.3535533905932738,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.3535533905932738,
+    0.7071067811865476,
+    0.3535533905932738,
+    0.0},
+    // bwdH
+    {0.1767766952966369,
+    0.3535533905932738,
+    -1.0606601717798214,
+    0.3535533905932738,
+    0.1767766952966369,
+    0.0}
+};
+
+DWTFilter DWTFilters::bior24 = {
+    // fwdL
+    {0.03314563036811942,
+    -0.06629126073623884,
+    -0.1767766952966369,
+    0.4198446513295126,
+    0.9943689110435825,
+    0.4198446513295126,
+    -0.1767766952966369,
+    -0.06629126073623884,
+    0.03314563036811942,
+    0.0},
+    // fwdH
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.3535533905932738,
+    -0.7071067811865476,
+    0.3535533905932738,
+    0.0,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.3535533905932738,
+    0.7071067811865476,
+    0.3535533905932738,
+    0.0,
+    0.0,
+    0.0},
+    // bwdH
+    {-0.03314563036811942,
+    -0.06629126073623884,
+    0.1767766952966369,
+    0.4198446513295126,
+    -0.9943689110435825,
+    0.4198446513295126,
+    0.1767766952966369,
+    -0.06629126073623884,
+    -0.03314563036811942,
+    0.0}
+};
+
+DWTFilter DWTFilters::bior31 = {
+    // fwdL
+    {-0.3535533905932738,
+    1.0606601717798214,
+    1.0606601717798214,
+    -0.3535533905932738},
+    // fwdH
+    {0.1767766952966369,
+    -0.5303300858899107,
+    0.5303300858899107,
+    -0.1767766952966369},
+    // bwdL
+    {0.1767766952966369,
+    0.5303300858899107,
+    0.5303300858899107,
+    0.1767766952966369},
+    // bwdH
+    {0.3535533905932738,
+    1.0606601717798214,
+    -1.0606601717798214,
+    -0.3535533905932738}
+};
+
+DWTFilter DWTFilters::bior33 = {
+    // fwdL
+    {0.06629126073623884,
+    -0.19887378220871652,
+    -0.15467960838455727,
+    0.9943689110435825,
+    0.9943689110435825,
+    -0.15467960838455727,
+    -0.19887378220871652,
+    0.06629126073623884},
+    // fwdH
+    {0.0,
+    0.0,
+    0.1767766952966369,
+    -0.5303300858899107,
+    0.5303300858899107,
+    -0.1767766952966369,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.1767766952966369,
+    0.5303300858899107,
+    0.5303300858899107,
+    0.1767766952966369,
+    0.0,
+    0.0},
+    // bwdH
+    {-0.06629126073623884,
+    -0.19887378220871652,
+    0.15467960838455727,
+    0.9943689110435825,
+    -0.9943689110435825,
+    -0.15467960838455727,
+    0.19887378220871652,
+    0.06629126073623884}
+};
+
+DWTFilter DWTFilters::bior35 = {
+    // fwdL
+    {-0.013810679320049757,
+    0.04143203796014927,
+    0.052480581416189075,
+    -0.26792717880896527,
+    -0.07181553246425874,
+    0.966747552403483,
+    0.966747552403483,
+    -0.07181553246425874,
+    -0.26792717880896527,
+    0.052480581416189075,
+    0.04143203796014927,
+    -0.013810679320049757},
+    // fwdH
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.1767766952966369,
+    -0.5303300858899107,
+    0.5303300858899107,
+    -0.1767766952966369,
+    0.0,
+    0.0,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.1767766952966369,
+    0.5303300858899107,
+    0.5303300858899107,
+    0.1767766952966369,
+    0.0,
+    0.0,
+    0.0,
+    0.0},
+    // bwdH
+    {0.013810679320049757,
+    0.04143203796014927,
+    -0.052480581416189075,
+    -0.26792717880896527,
+    0.07181553246425874,
+    0.966747552403483,
+    -0.966747552403483,
+    -0.07181553246425874,
+    0.26792717880896527,
+    0.052480581416189075,
+    -0.04143203796014927,
+    -0.013810679320049757}
+};
--- a/lib/Numerical/DWTFilters.hpp
+++ b/lib/Numerical/DWTFilters.hpp
@@ -0,0 +1,53 @@
+/*
+ * DWTFilters.hpp, part of LatAnalyze
+ *
+ * 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
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 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.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_DWTFilters_hpp_
+#define Latan_DWTFilters_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+struct DWTFilter
+{
+    const std::vector<double> fwdL, fwdH, bwdL, bwdH;
+};
+
+namespace DWTFilters
+{
+    extern DWTFilter haar;
+    extern DWTFilter db2;
+    extern DWTFilter db3;
+    extern DWTFilter db4;
+    extern DWTFilter db5;
+    extern DWTFilter db6;
+    extern DWTFilter bior13;
+    extern DWTFilter bior15;
+    extern DWTFilter bior22;
+    extern DWTFilter bior24;
+    extern DWTFilter bior31;
+    extern DWTFilter bior33;
+    extern DWTFilter bior35;
+
+    extern std::map<std::string, const DWTFilter *> fromName;
+}
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_DWTFilters_hpp_
--- a/lib/Numerical/Derivative.cpp
+++ b/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
--- a/lib/Numerical/Derivative.hpp
+++ b/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
--- a/lib/Numerical/FFT.hpp
+++ b/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
--- a/lib/Numerical/GslFFT.cpp
+++ b/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_);
--- a/lib/Numerical/GslFFT.hpp
+++ b/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
--- a/lib/Numerical/GslHybridRootFinder.cpp
+++ b/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
--- a/lib/Numerical/GslHybridRootFinder.hpp
+++ b/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
--- a/lib/Numerical/GslMinimizer.cpp
+++ b/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
--- a/lib/Numerical/GslMinimizer.hpp
+++ b/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
--- a/lib/Numerical/GslQagsIntegrator.cpp
+++ b/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
--- a/lib/Numerical/GslQagsIntegrator.hpp
+++ b/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
--- a/lib/Numerical/Integrator.hpp
+++ b/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
--- a/lib/Numerical/Minimizer.cpp
+++ b/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
--- a/lib/Numerical/Minimizer.hpp
+++ b/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
--- a/lib/Numerical/MinuitMinimizer.cpp
+++ b/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
@@ -19,6 +19,20 @@

 #include <LatAnalyze/Numerical/MinuitMinimizer.hpp>
 #include <LatAnalyze/includes.hpp>
+
+// forward declaration necessary in the ROOT-based version of Minuit2
+namespace ROOT
+{
+    namespace Fit
+    {
+        class ParameterSettings;
+    };
+};
+
+// macros necessary in the ROOT-based version of Minuit2
+#define ROOT_Math_VecTypes
+#define MATHCORE_STANDALONE
+
 #include <Minuit2/Minuit2Minimizer.h>
 #include <Math/Functor.h>

--- a/lib/Numerical/MinuitMinimizer.hpp
+++ b/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
--- a/lib/Numerical/NloptMinimizer.cpp
+++ b/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
--- a/lib/Numerical/NloptMinimizer.hpp
+++ b/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
--- a/lib/Numerical/RootFinder.cpp
+++ b/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
--- a/lib/Numerical/RootFinder.hpp
+++ b/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
--- a/lib/Numerical/Solver.cpp
+++ b/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
--- a/lib/Numerical/Solver.hpp
+++ b/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
--- a/lib/Physics/CorrelatorFitter.cpp
+++ b/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);
+    }
+}
--- a/lib/Physics/CorrelatorFitter.hpp
+++ b/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_
--- a/lib/Physics/EffectiveMass.cpp
+++ b/lib/Physics/EffectiveMass.cpp
@@ -0,0 +1,132 @@
+/*
+ * EffectiveMass.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <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;
+}
--- a/lib/Physics/EffectiveMass.hpp
+++ b/lib/Physics/EffectiveMass.hpp
@@ -0,0 +1,50 @@
+/*
+ * EffectiveMass.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <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_
--- a/lib/Statistics/Dataset.hpp
+++ b/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
--- a/lib/Statistics/FitInterface.cpp
+++ b/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
--- a/lib/Statistics/FitInterface.hpp
+++ b/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
--- a/lib/Statistics/Histogram.cpp
+++ b/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
--- a/lib/Statistics/Histogram.hpp
+++ b/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
--- a/lib/Statistics/MatSample.hpp
+++ b/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
@@ -103,6 +103,10 @@ public:
                     const Index nCol);
    // resize all matrices
    void resizeMat(const Index nRow, const Index nCol);
+    // covariance matrix
+    Mat<T> covarianceMatrix(const MatSample<T> &sample) const;
+    Mat<T> varianceMatrix(void) const;
+    Mat<T> correlationMatrix(void) const;
 };

 // non-member operators
@@ -379,6 +383,78 @@ void MatSample<T>::resizeMat(const Index nRow, const Index nCol)
    }
 }

+// covariance matrix ///////////////////////////////////////////////////////////
+template <typename T>
+Mat<T> MatSample<T>::covarianceMatrix(const MatSample<T> &sample) const
+{
+    if (((*this)[central].cols() != 1) or (sample[central].cols() != 1))
+    {
+        LATAN_ERROR(Size, "samples have more than one column");
+    }
+
+    Index  n1 = (*this)[central].rows(), n2 = sample[central].rows();
+    Index  nSample = this->size();
+    Mat<T> tmp1(n1, nSample), tmp2(n2, nSample), res(n1, n2);
+    Mat<T> s1(n1, 1), s2(n2, 1), one(nSample, 1);
+
+    one.fill(1.);
+    s1.fill(0.);
+    s2.fill(0.);
+    for (unsigned int s = 0; s < nSample; ++s)
+    {
+        s1          += (*this)[s];
+        tmp1.col(s)  = (*this)[s];
+    }
+    tmp1 -= s1*one.transpose()/static_cast<double>(nSample);
+    for (unsigned int s = 0; s < nSample; ++s)
+    {
+        s2          += sample[s];
+        tmp2.col(s)  = sample[s];
+    }
+    tmp2 -= s2*one.transpose()/static_cast<double>(nSample);
+    res   = tmp1*tmp2.transpose()/static_cast<double>(nSample - 1);
+
+    return res;
+}
+
+template <typename T>
+Mat<T> MatSample<T>::varianceMatrix(void) const
+{
+    if ((*this)[central].cols() != 1)
+    {
+        LATAN_ERROR(Size, "samples have more than one column");
+    }
+
+    Index  n1 = (*this)[central].rows();
+    Index  nSample = this->size();
+    Mat<T> tmp1(n1, nSample), res(n1, n1);
+    Mat<T> s1(n1, 1), one(nSample, 1);
+
+    one.fill(1.);
+    s1.fill(0.);
+    for (unsigned int s = 0; s < nSample; ++s)
+    {
+        s1          += (*this)[s];
+        tmp1.col(s)  = (*this)[s];
+    }
+    tmp1 -= s1*one.transpose()/static_cast<double>(nSample);
+    res   = tmp1*tmp1.transpose()/static_cast<double>(nSample - 1);
+
+    return res;
+}
+
+template <typename T>
+Mat<T> MatSample<T>::correlationMatrix(void) const
+{
+    Mat<T> res = varianceMatrix();
+    Mat<T> invDiag(res.rows(), 1);
+
+    invDiag = res.diagonal();
+    invDiag = invDiag.cwiseInverse().cwiseSqrt();
+    res     = (invDiag*invDiag.transpose()).cwiseProduct(res);
+
+    return res;
+}

 END_LATAN_NAMESPACE

--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Matt Spraggs	b796bfbd68	Merge `269d0c338e` into `05138baa08`	2023-06-20 15:48:04 +01:00
Antonin Portelli	05138baa08	fix of potential conflict with std::size	2023-06-20 15:47:58 +01:00
Antonin Portelli	a0bdbfd9dd	Update build-macos.yml	2023-03-20 14:04:17 +00:00
Antonin Portelli	7fd31d1fcc	object to plot single point	2022-08-03 14:24:36 +01:00
Antonin Portelli	f0c3fd4d7d	PlotData fix	2022-08-03 14:24:20 +01:00
Antonin Portelli	470aff3b4a	Laplace filter in 2-pt fit	2022-04-29 16:45:02 +01:00
Antonin Portelli	c9ea23dc92	Merge branch 'master' into develop	2022-03-10 08:34:53 +00:00
Antonin Portelli	58a355478a	Merge branch 'feature/dwt' into develop	2022-03-10 08:21:40 +00:00
Antonin Portelli	4f919bc007	tool to merge samples	2022-03-10 08:21:26 +00:00
Antonin Portelli	9455e2c66e	correlation matrix plot compute dynamic range	2022-03-10 08:21:10 +00:00
Antonin Portelli	43dd295f94	data plots compatible with multi-column arrays	2022-03-10 08:20:38 +00:00
Antonin Portelli	9afd40a1ad	utility to compute sample DWT	2022-03-10 08:16:40 +00:00
Antonin Portelli	9e78b96260	DWT working and tested	2022-02-18 14:06:52 +00:00
Antonin Portelli	65a656f257	first skeleton for DWT	2022-02-17 19:24:34 +00:00
Antonin Portelli	47d0b3f040	correlation dynamic range renaming	2022-02-16 19:03:19 +00:00
Antonin Portelli	35f6733292	sample-plot-corr better palette	2022-02-16 18:55:24 +00:00
Antonin Portelli	ebc1bd4c2e	fit: stable variance inversion and SVD dynamic range	2022-02-16 18:55:08 +00:00
Antonin Portelli	857a8e59c9	corr to var and SVD dynamic range	2022-02-16 18:54:16 +00:00
Antonin Portelli	0de8091f3c	Eigen plugin const	2022-02-16 18:53:45 +00:00
Antonin Portelli	e4cefae515	special macro to plot correlation matrices	2021-12-28 12:17:02 +01:00
Antonin Portelli	8cd29c2bee	XYStatData fixes	2021-12-26 22:21:03 +01:00
Antonin Portelli	bac8356de5	Ignore verbose minimiser for samples	2021-12-21 18:18:14 +01:00
Antonin Portelli	60d91cbff5	plot data with points	2021-12-20 01:30:26 +01:00
Antonin Portelli	adf2c9cc69	normalised residuals routines	2021-12-20 01:30:03 +01:00
Antonin Portelli	24a7b9c203	remove new covariance routine regression code	2021-12-20 01:29:29 +01:00
Antonin Portelli	57c6004797	significant optimisation of covariance routines + checks	2021-12-20 01:25:13 +01:00
Antonin Portelli	c796187d1e	Absolute value plot (useful for log scale)	2021-12-14 13:13:36 +00:00
Antonin Portelli	b92fb84e9d	adding pthread link	2021-11-29 00:03:19 +00:00
Antonin Portelli	5e04a0321e	Merge branch 'master' into develop	2021-11-28 23:57:15 +00:00
Antonin Portelli	78351a9b76	better interface for critical threaded sections	2021-11-28 23:51:56 +00:00
Antonin Portelli	fe8c6c6630	Merge branch 'develop' of github.com:aportelli/LatAnalyze into develop	2021-11-28 23:26:20 +00:00
Antonin Portelli	5f192ad30f	Thread pool implementation	2021-11-28 23:26:17 +00:00
Matt Spraggs	269d0c338e	MathInterpreter: Replace Instruction polymorphism with bytecode	2021-11-19 23:59:31 +00:00
Antonin Portelli	ccb837a244	Update build-macos.yml	2021-11-17 16:37:05 +00:00
Antonin Portelli	499e173bac	Update Readme.md	2021-11-17 16:20:56 +00:00
Antonin Portelli	75485219d8	CI uses HDF5 1.10.8 hoping to solve compilation issues on Big Sur	2021-11-17 16:17:45 +00:00
Antonin Portelli	a3054a0f44	first GitHub CI	2021-11-17 16:06:26 +00:00
Antonin Portelli	d6e5ba724d	CI scripts expose number of build tasks	2021-11-16 21:35:35 +00:00
Antonin Portelli	9341a31cf4	gitignore update	2021-11-16 21:34:39 +00:00
Antonin Portelli	b4b6bd22fa	compatible with Minuit2 in ROOT	2021-11-16 21:34:25 +00:00
Antonin Portelli	68d22eca11	Badges update	2020-11-27 17:09:31 +00:00
Antonin Portelli	d4704267d6	Dev version number	2020-11-27 13:34:58 +00:00
Antonin Portelli	2de70a2775	Merge branch 'release/v3.5.1'	2020-11-27 13:34:26 +00:00
Antonin Portelli	d67a25245e	Merge tag 'v3.5.1' into develop	2020-11-27 13:34:26 +00:00
Antonin Portelli	cfb2c0c5e8	version number bump	2020-11-27 13:33:39 +00:00
Antonin Portelli	0160c88c29	bootstrap script fix	2020-11-17 11:56:10 +00:00
Antonin Portelli	376bdfc38b	Remove archaic error is DISPLAY is not set	2020-11-17 11:50:52 +00:00
Antonin Portelli	0cec36dded	Eigen update	2020-11-17 11:20:16 +00:00
Antonin Portelli	0cda9e20cd	Warning suppression	2020-11-17 11:19:51 +00:00
Antonin Portelli	2b9508c20f	Merge pull request #18 from AndrewYongZhenNing/develop Changed linear model's parameters	2020-05-08 18:53:05 +01:00
Andrew Zhen Ning Yong	417e068485	Swapped p[0] -> p[1] such that p[0] is always the effective mass parameter.	2020-05-07 17:08:12 +01:00
Antonin Portelli	524c11d2ba	Merge pull request #17 from AndrewYongZhenNing/develop Develop	2020-04-30 17:22:10 +01:00
Andrew Zhen Ning Yong	3602bbd368	Merging with upstream before pull request.	2020-04-30 10:50:54 +01:00
Andrew Zhen Ning Yong	bb4e9e1d42	Added else if block for CorrelatorType::linear to set its limits.	2020-04-30 10:50:08 +01:00
Antonin Portelli	4347511e39	Merge pull request #16 from AndrewYongZhenNing/develop For 'cst' model in makeConstModel(), changed bug from x[0] -> p[0] .	2020-02-21 15:22:22 +01:00
Andrew Zhen Ning Yong	e944f9c4aa	For 'cst' model, changed bug from x[0] -> p[0].	2020-02-21 11:18:02 +00:00
Antonin Portelli	3e70792a06	correlator transform utilities	2020-02-02 19:43:14 +01:00
Antonin Portelli	f014003593	minor FFT leak	2020-02-02 19:42:54 +01:00
Antonin Portelli	c37e6e1bfd	effective mass fix	2020-01-28 20:15:28 +00:00
Antonin Portelli	0e8b9d2a8f	Merge branch 'feature/correlator-fitter' into develop	2020-01-28 17:56:08 +00:00
Antonin Portelli	1d6a66263d	code cleaning	2020-01-28 17:55:47 +00:00
Antonin Portelli	1775f4992b	rewrite of the 2pt fitter using the new physics classes	2020-01-28 17:35:07 +00:00
Antonin Portelli	685d433032	minor fixes in new physics classes	2020-01-28 17:34:37 +00:00
Antonin Portelli	1bde8822b1	band plot does not rely anymore on tac & tail -r commands	2020-01-28 17:34:00 +00:00
Antonin Portelli	f826f30e82	first stab at correlator fit utility classes	2020-01-23 19:11:01 +00:00
Antonin Portelli	51efb2a81f	nlopt CI linking fix	2020-01-13 16:40:04 +00:00
Antonin Portelli	5f5ecf241f	Copyright year update	2020-01-13 09:57:06 +00:00
Antonin Portelli	d85860a2ef	CI scripts update	2020-01-11 20:16:56 +00:00
Antonin Portelli	d894aa185c	bigger font in plots	2020-01-11 17:23:02 +00:00
Antonin Portelli	0fbe00da0d	gnuplot path search fix	2020-01-11 16:48:36 +00:00
Antonin Portelli	caeb78b143	more default options when trying to localise gnuplot	2020-01-11 16:34:05 +00:00
Antonin Portelli	9f98ed42c3	Haswell is the default architecture for CI builds	2020-01-10 17:10:33 +00:00
Antonin Portelli	d43197ccc7	public access to gnuplot path	2020-01-10 15:15:46 +00:00
Antonin Portelli	4b5ad9014c	more script tuning	2020-01-10 14:15:40 +00:00
Antonin Portelli	2e2e676196	more CI script tuning	2020-01-10 12:09:39 +00:00
Antonin Portelli	c0dac8063e	scripts fine-tuning	2020-01-10 11:59:07 +00:00
Antonin Portelli	29863a348b	CI build script minor fixes	2020-01-09 17:35:53 +00:00
Antonin Portelli	a1fae4356e	CI scripts build a static version of NLOpt	2020-01-09 13:43:52 +00:00
Antonin Portelli	5e891063e1	plot fixes when reploting the same Plot object	2020-01-09 13:43:27 +00:00
Antonin Portelli	97267c196f	saving PDF is optional when saving plot	2020-01-09 13:42:59 +00:00
Antonin Portelli	c81316ef32	development version number	2019-12-12 18:24:29 +00:00
Antonin Portelli	267bd33a97	Merge tag '3.5' into develop no message	2019-12-12 18:23:50 +00:00