Merge pull request #19 from AndrewYongZhenNing/develop

Corrected Sinh model's definition, init fit parameter & fold definition.

.github/workflows/build-macos.yml

@@ -0,0 +1,26 @@
+name: Build macOS
+
+on: [push]
+
+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

.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

.gitignore

@@ -16,6 +16,7 @@ autom4te.cache/*
 *.in~
 config.h*
 configure
+configure~
 .buildutils/*
 aclocal.m4
 

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

bootstrap.sh

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

ci-scripts/install-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

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

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

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

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

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

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

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

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;
+}

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;
+}

lib/Core/Eigen.hpp

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

lib/Core/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();
     

lib/Core/Exceptions.cpp

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

lib/Core/Exceptions.hpp

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

lib/Core/Mat.cpp

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

lib/Core/Mat.hpp

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

lib/Core/Math.cpp

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

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;

lib/Core/MathInterpreter.cpp

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

lib/Core/MathInterpreter.hpp

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

lib/Core/MathLexer.lpp

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

lib/Core/MathParser.ypp

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

lib/Core/ParserState.hpp

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

lib/Core/Plot.cpp

@@ -1,7 +1,7 @@
 /*
  * Plot.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -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;
@@ -692,23 +807,26 @@ void Plot::save(string dirName)
     
     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))
     {
         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;

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:

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();
+}

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

lib/Core/Utilities.cpp

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

lib/Core/Utilities.hpp

@@ -1,7 +1,7 @@
 /*
  * Utilities.hpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2015 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -108,6 +108,23 @@ inline std::string strFrom(const T x)
 }
 
 // specialization for vectors
+template<>
+inline std::vector<Index> strTo<std::vector<Index>>(const std::string &str)
+{
+    std::vector<Index>  res;
+    std::vector<double> vbuf;
+    double              buf;
+    std::istringstream  stream(str);
+    
+    while (!stream.eof())
+    {
+        stream >> buf;
+        res.push_back(buf);
+    }
+    
+    return res;
+}
+
 template<>
 inline DVec strTo<DVec>(const std::string &str)
 {

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>

lib/Functional/CompiledFunction.cpp

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

lib/Functional/CompiledFunction.hpp

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

lib/Functional/CompiledModel.cpp

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

lib/Functional/CompiledModel.hpp

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

lib/Functional/Function.cpp

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

lib/Functional/Function.hpp

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

lib/Functional/Model.cpp

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

lib/Functional/Model.hpp

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

lib/Functional/TabFunction.cpp

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

lib/Functional/TabFunction.hpp

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

lib/Global.cpp

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

lib/Global.hpp

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

lib/Io/AsciiFile.cpp

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

lib/Io/AsciiFile.hpp

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

lib/Io/AsciiLexer.lpp

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

lib/Io/AsciiParser.ypp

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

lib/Io/File.cpp

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

lib/Io/File.hpp

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

lib/Io/Hdf5File.cpp

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

lib/Io/Hdf5File.hpp

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

lib/Io/Io.cpp

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

lib/Io/Io.hpp

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

lib/Io/IoObject.hpp

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

lib/Io/XmlReader.cpp

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

lib/Io/XmlReader.hpp

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

lib/Makefile.am

@@ -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         \

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;
+}

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_

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}
+};

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_

lib/Numerical/Derivative.cpp

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

lib/Numerical/Derivative.hpp

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

lib/Numerical/FFT.hpp

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

lib/Numerical/GslFFT.cpp

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

lib/Numerical/GslFFT.hpp

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

lib/Numerical/GslHybridRootFinder.cpp

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

lib/Numerical/GslHybridRootFinder.hpp

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

lib/Numerical/GslMinimizer.cpp

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

lib/Numerical/GslMinimizer.hpp

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

lib/Numerical/GslQagsIntegrator.cpp

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

lib/Numerical/GslQagsIntegrator.hpp

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

lib/Numerical/Integrator.hpp

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

lib/Numerical/Minimizer.cpp

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

lib/Numerical/Minimizer.hpp

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

lib/Numerical/MinuitMinimizer.cpp

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

lib/Numerical/MinuitMinimizer.hpp

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

lib/Numerical/NloptMinimizer.cpp

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

lib/Numerical/NloptMinimizer.hpp

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

lib/Numerical/RootFinder.cpp

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

lib/Numerical/RootFinder.hpp

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

lib/Numerical/Solver.cpp

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

lib/Numerical/Solver.hpp

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

lib/Physics/CorrelatorFitter.cpp

@@ -0,0 +1,440 @@
+/*
+ * 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 CorrelatorModels::ModelPar &par)
+{
+    const Index nt  = c[central].rows();
+    DMatSample  buf = c;
+    int         sign;
+    bool        fold = false;
+    
+    switch (par.type)
+    {
+        case CorrelatorType::cosh:
+        case CorrelatorType::cst:
+            sign = 1;
+            fold = true;
+            break;
+        case CorrelatorType::sinh:
+            sign = -1;
+            fold = true;
+            break;
+        case CorrelatorType::linear:
+            cout << "Linear model is asymmetric: will not fold." << endl;
+            break;
+        default:
+            break;
+    }
+
+    if (fold)
+    {
+        FOR_STAT_ARRAY(buf, s)
+        {
+            for (Index t = 0; t < nt; ++t)
+            {
+                buf[s](t) = 0.5*(c[s](t) + sign*c[s]((nt - t) % nt));
+            }
+        }
+    }
+
+    return buf;
+}
+
+DMatSample CorrelatorUtils::fourierTransform(const DMatSample &c, FFT &fft, 
+                                             const unsigned int dir)
+{
+    const Index nSample   = c.size();
+    const Index nt        = c[central].rows();
+    bool        isComplex = (c[central].cols() > 1);
+    CMatSample  buf(nSample, nt, 1);
+    DMatSample  out(nSample, nt, 2);
+
+    fft.resize(nt);
+    FOR_STAT_ARRAY(buf, s)
+    {
+        buf[s].real() = c[s].col(0);
+        if (isComplex)
+        {
+            buf[s].imag() = c[s].col(1);
+        }
+        else
+        {
+            buf[s].imag() = DVec::Constant(nt, 0.);
+        }
+        fft(buf[s], dir);
+        out[s].col(0) = buf[s].real();
+        out[s].col(1) = buf[s].imag();
+    }
+
+    return out;    
+}
+
+/******************************************************************************
+ *                      CorrelatorFitter implementation                       *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+CorrelatorFitter::CorrelatorFitter(const DMatSample &corr)
+{
+    setCorrelator(corr);
+}
+
+CorrelatorFitter::CorrelatorFitter(const std::vector<DMatSample> &corr)
+{
+    setCorrelators(corr);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+XYSampleData & CorrelatorFitter::data(void)
+{
+    return *data_;
+}
+
+void CorrelatorFitter::setCorrelator(const DMatSample &corr)
+{
+    std::vector<DMatSample> vec;
+
+    vec.push_back(corr);
+    setCorrelators(vec);
+}   
+
+void CorrelatorFitter::setCorrelators(const std::vector<DMatSample> &corr)
+{
+    Index                           nSample = corr[0].size();
+    DMatSample                      tVec(nSample);
+    std::vector<const DMatSample *> ptVec;
+
+    nt_ = corr[0][central].rows();
+    tVec.fill(DVec::LinSpaced(nt_, 0, nt_ - 1));
+    for (auto &c: corr)
+    {
+        ptVec.push_back(&c);
+    }
+    data_.reset(new XYSampleData(corr[0].size()));
+    data_->addXDim(nt_, "t/a", true);
+    for (unsigned int i = 0; i < corr.size(); ++i)
+    {
+        data_->addYDim("C_" + strFrom(i) + "(t)");
+    }
+    data_->setUnidimData(tVec, ptVec);
+    model_.resize(corr.size());
+    range_.resize(corr.size(), make_pair(0, nt_ - 1));
+    thinning_.resize(corr.size(), 1);
+}
+
+void CorrelatorFitter::setModel(const DoubleModel &model, const Index i)
+{
+    model_[i] = model;
+}
+
+const DoubleModel & CorrelatorFitter::getModel(const Index i) const
+{
+    return model_.at(i);
+}
+
+void CorrelatorFitter::setFitRange(const Index tMin, const Index tMax, 
+                                   const Index i)
+{
+    range_[i] = make_pair(tMin, tMax);
+    refreshRanges();
+}
+
+void CorrelatorFitter::setCorrelation(const bool isCorrelated, const Index i, 
+                                      const Index j)
+{
+    data_->assumeYYCorrelated(isCorrelated, i, j);
+}
+
+DMat CorrelatorFitter::getVarianceMatrix(void) const
+{
+    return data_->getFitVarMat();
+}
+
+void CorrelatorFitter::setThinning(const Index thinning, const Index i)
+{
+    thinning_[i] = thinning;
+    refreshRanges();
+}
+
+// fit functions ///////////////////////////////////////////////////////////////
+SampleFitResult CorrelatorFitter::fit(Minimizer &minimizer, const DVec &init)
+{
+    vector<Minimizer *> vecPt = {&minimizer};
+
+    return fit(vecPt, init);
+}
+
+SampleFitResult CorrelatorFitter::fit(vector<Minimizer *> &minimizer,
+                                      const DVec &init)
+{
+    vector<const DoubleModel *> vecPt(model_.size());
+    
+    for (unsigned int i = 0; i < model_.size(); ++i)
+    {
+        vecPt[i] = &(model_[i]);
+    }
+
+    return data_->fit(minimizer, init, vecPt);
+}
+
+// internal function to refresh fit ranges /////////////////////////////////////
+void CorrelatorFitter::refreshRanges(void)
+{
+    for (unsigned int i = 0; i < range_.size(); ++i)
+    for (Index t = 0; t < nt_; ++t)
+    {
+        data_->fitPoint((t >= range_[i].first) and (t <= range_[i].second)
+                         and ((t - range_[i].first) % thinning_[i] == 0), t);
+    }
+}

lib/Physics/CorrelatorFitter.hpp

@@ -0,0 +1,104 @@
+/*
+ * CorrelatorFitter.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_CorrelatorFitter_hpp_
+#define Latan_CorrelatorFitter_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Functional/Model.hpp>
+#include <LatAnalyze/Numerical/FFT.hpp>
+#include <LatAnalyze/Statistics/XYSampleData.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                         Correlator types & models                          *
+ ******************************************************************************/
+enum class CorrelatorType {undefined, exp, cosh, sinh, linear, cst};
+
+namespace CorrelatorModels
+{
+    struct ModelPar
+    {
+        CorrelatorType type;
+        Index          nState;
+    };
+
+    DoubleModel makeExpModel(const Index nState);
+    DoubleModel makeCoshModel(const Index nState, const Index nt);
+    DoubleModel makeSinhModel(const Index nState, const Index nt);
+    DoubleModel makeConstModel(void);
+    DoubleModel makeLinearModel(void);
+    ModelPar    parseModel(const std::string s);
+    DoubleModel makeModel(const ModelPar par, const Index nt);
+    DVec        parameterGuess(const DMatSample &corr, const ModelPar par);
+};
+
+/******************************************************************************
+ *                         Correlator utilities                               *
+ ******************************************************************************/
+namespace CorrelatorUtils
+{
+    DMatSample shift(const DMatSample &c, const Index ts);
+    DMatSample fold(const DMatSample &c, const CorrelatorModels::ModelPar &par);
+    DMatSample fourierTransform(const DMatSample &c, FFT &fft, 
+                                const unsigned int dir = FFT::Forward);
+};
+
+/******************************************************************************
+ *                       Correlator fit utility class                         *
+ ******************************************************************************/
+class CorrelatorFitter
+{
+public:
+    // constructors
+    CorrelatorFitter(const DMatSample &corr);
+    CorrelatorFitter(const std::vector<DMatSample> &corr);
+    // destructor
+    virtual ~CorrelatorFitter(void) = default;
+    // access
+    XYSampleData & data(void);
+    void setCorrelator(const DMatSample &corr);
+    void setCorrelators(const std::vector<DMatSample> &corr);
+    const DMatSample & getCorrelator(const Index i = 0) const;
+    const std::vector<DMatSample> & getCorrelators(void) const;
+    void setModel(const DoubleModel &model, const Index i = 0);
+    const DoubleModel & getModel(const Index i = 0) const;
+    void setFitRange(const Index tMin, const Index tMax, const Index i = 0);
+    void setCorrelation(const bool isCorrelated, const Index i = 0, 
+                        const Index j = 0);
+    DMat getVarianceMatrix(void) const;
+    void setThinning(const Index thinning, const Index i = 0);
+    // fit functions
+    SampleFitResult fit(Minimizer &minimizer, const DVec &init);
+    SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init);
+private:
+    // internal function to refresh fit ranges
+    void refreshRanges(void);
+private:
+    Index                                nt_{0};
+    std::unique_ptr<XYSampleData>        data_;
+    std::vector<DoubleModel>             model_;
+    std::vector<std::pair<Index, Index>> range_;
+    std::vector<Index>                   thinning_;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_CorrelatorFitter_hpp_

lib/Physics/EffectiveMass.cpp

@@ -0,0 +1,132 @@
+/*
+ * EffectiveMass.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Physics/EffectiveMass.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                       EffectiveMass implementation                         *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+EffectiveMass::EffectiveMass(const CorrelatorType type)
+{
+    setType(type);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+CorrelatorType EffectiveMass::getType(void) const
+{
+    return type_;
+}
+
+void EffectiveMass::setType(const CorrelatorType type)
+{
+    type_ = type;
+}
+
+DVec EffectiveMass::getTime(const Index nt) const
+{
+    DVec tvec;
+
+    switch (type_)
+    {
+    case CorrelatorType::undefined:
+        LATAN_ERROR(Argument, "correlator type is undefined");
+        break;
+    case CorrelatorType::exp:
+    case CorrelatorType::linear:
+        tvec = DVec::LinSpaced(nt - 1, 0, nt - 2);
+        break;
+    case CorrelatorType::cosh:
+    case CorrelatorType::sinh:
+        tvec = DVec::LinSpaced(nt - 2, 1, nt - 2);
+        break;
+    case CorrelatorType::cst:
+        tvec = DVec::LinSpaced(nt, 0, nt - 1);
+        break;
+    }
+
+    return tvec;
+}
+
+// compute effective mass //////////////////////////////////////////////////////
+DVec EffectiveMass::operator()(const DVec &corr) const
+{
+    Index nt = corr.size();
+    DVec  em;
+
+    if (nt < 2)
+    {
+        LATAN_ERROR(Size, "input vector has less than 2 elements");
+    }
+    switch (type_)
+    {
+    case CorrelatorType::undefined:
+        LATAN_ERROR(Argument, "correlator type is undefined");
+        break;
+    case CorrelatorType::exp:
+        em.resize(nt - 1);
+        for (Index t = 1; t < nt; ++t)
+        {
+            em(t - 1) = log(corr(t - 1)/corr(t));
+        }
+        break;
+    case CorrelatorType::cosh:
+        em.resize(nt - 2);
+        for (Index t = 1; t < nt - 1; ++t)
+        {
+            em(t - 1) = acosh((corr(t - 1) + corr(t + 1))/(2.*corr(t)));
+        }
+        break;
+    case CorrelatorType::sinh:
+        em.resize(nt - 2);
+        for (Index t = 1; t < nt - 1; ++t)
+        {
+            em(t - 1) = acosh((corr(t - 1) + corr(t + 1))/(2.*corr(t)));
+        }
+        break;
+    case CorrelatorType::linear:
+        em.resize(nt - 1);
+        for (Index t = 1; t < nt; ++t)
+        {
+            em(t - 1) = corr(t) - corr(t - 1);
+        }
+        break;
+    case CorrelatorType::cst:
+        em = corr;
+        break;
+    }
+
+    return em;
+}
+
+DMatSample EffectiveMass::operator()(const DMatSample &corr) const
+{
+    DMatSample em(corr.size());
+
+    FOR_STAT_ARRAY(em, s)
+    {
+        em[s] = (*this)(corr[s]);
+    }
+
+    return em;
+}

lib/Physics/EffectiveMass.hpp

@@ -0,0 +1,50 @@
+/*
+ * EffectiveMass.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_EffectiveMass_hpp_
+#define Latan_EffectiveMass_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
+#include <LatAnalyze/Physics/CorrelatorFitter.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                          Effective mass class                              *
+ ******************************************************************************/
+class EffectiveMass
+{
+public:
+    // constructors
+    EffectiveMass(const CorrelatorType type = CorrelatorType::exp);
+    // access
+    CorrelatorType getType(void) const;
+    void           setType(const CorrelatorType type);
+    DVec           getTime(const Index nt) const;
+    // compute effective mass
+    DVec       operator()(const DVec &corr) const;
+    DMatSample operator()(const DMatSample &corr) const;
+private:
+    CorrelatorType type_;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_EffectiveMass_hpp_

lib/Statistics/Dataset.hpp

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

lib/Statistics/FitInterface.cpp

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

lib/Statistics/FitInterface.hpp

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

lib/Statistics/Histogram.cpp

@@ -1,7 +1,7 @@
 /*
  * Histogram.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -146,6 +146,16 @@ double Histogram::getX(const Index i) const
     return x_(i);
 }
 
+double Histogram::getXMin(void) const
+{
+    return xMin_;
+}
+
+double Histogram::getXMax(void) const
+{
+    return xMax_;
+}
+
 double Histogram::operator[](const Index i) const
 {
     return bin_(i)*(isNormalized() ? norm_ : 1.);

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
@@ -52,6 +52,8 @@ public:
     const StatArray<double> & getData(void) const;
     const StatArray<double> & getWeight(void) const;
     double                    getX(const Index i) const;
+    double                    getXMin(void) const;
+    double                    getXMax(void) const;
     double                    operator[](const Index i) const;
     double                    operator()(const double x) const;
     // percentiles & confidence interval

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
 

lib/Statistics/Random.cpp

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