update build helper scripts

.github/workflows/build-macos.yml

@@ -0,0 +1,26 @@
+name: Build macOS
+
+on: [push, workflow_dispatch]
+
+jobs:
+  build:
+
+    runs-on: macos-11
+
+    steps:
+    - name: Checkout
+      uses: actions/checkout@v2
+    - name: Install basic dependencies
+      run: brew install automake autoconf libtool bison flex
+    - name: Build dependencies
+      shell: bash
+      run: |
+        export PATH=/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:${PATH}
+        cd ci-scripts
+        ./install-deps.sh prefix 6
+    - name: Build LatAnalyze
+      shell: bash
+      run: |
+        export PATH=/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:${PATH}
+        cd ci-scripts
+        ./install-latan.sh prefix 6

.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
 
@@ -24,9 +25,10 @@ lib/*Lexer.cpp
 lib/*Parser.cpp
 lib/*Parser.hpp
 
-# Eigen headers
-lib/Eigen/*
+# Eigen headers and archives
+lib/Eigen
 lib/eigen_files.mk
+eigen-*.tar.bz2
 
 # CI builds
 ci-scripts/local/*

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
@@ -164,4 +149,4 @@ Fixes:
 #### v3.0
 Commit `7b4f2884a5e99bbfab4d4bd7623f609a55403c39`.  
 First 'stable' version of LatAnalyze in C++. The v2.0 refers to the [C version](https://github.com/aportelli/LatAnalyze-legacy) and v1.0 to an old undistributed version.  
-**This version compiles fine on OS X with clang but does have many portability issues to other platforms/compilers, v3.1 is the first real release.**
+**This version compiles fine on OS X with clang but does have many portability issues to other platforms/compilers, v3.1 is the first real release.**

bootstrap.sh

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

build.sh

@@ -1,24 +0,0 @@
-#!/bin/bash
-
-set -e
-
-PREFIX=`cat Makefile | grep '^prefix =' | awk '{print $3}'`
-case $1 in
-    '')
-        echo '-- building...'
-        make -j8
-        echo '-- installing...'
-        make uninstall 1>/dev/null
-        make install 1>/dev/null;;
-#        if [[ `basename \`pwd\`` == "lib" ]]
-#        then
-#            echo '-- creating debug symbols...'
-#            dsymutil .libs/libLatAnalyze.0.dylib -o ${PREFIX}/lib/libLatAnalyze.0.dylib.dSYM
-#        fi;;
-    'clean')
-        echo '-- cleaning...'
-        make -j8 clean;;
-    *)
-        echo 'error: unknown action' 1>&2
-        exit 1;;
-esac

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,11 +2,12 @@
 
 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 -P)
@@ -19,7 +20,7 @@ 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,12 +1,13 @@
 #!/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 -P)
@@ -19,7 +20,7 @@ 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,10 +1,11 @@
 #!/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
 INITDIR=$(pwd -P)
@@ -12,11 +13,11 @@ mkdir -p ${PREFIX}
 cd ${PREFIX}
 PREFIX=$(pwd -P)
 cd ${INITDIR}
-./install-deps.sh ${PREFIX}
+./install-deps.sh ${PREFIX} ${NTASKS}
 cd ..
 ./bootstrap.sh
 mkdir -p build
 cd build
-../configure --prefix=${PREFIX} --with-minuit=${PREFIX} --with-nlopt=${PREFIX} --with-hdf5=${PREFIX} --with-gsl=${PREFIX} CXXFLAGS="${CXXFLAGS} -O3 -march=haswell -mtune=haswell"
-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,12 +1,12 @@
 #!/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 -P)
@@ -14,12 +14,12 @@ mkdir -p ${PREFIX}
 cd ${PREFIX}
 PREFIX=$(pwd -P)
 cd ${INITDIR}/local/build
-wget http://www.cern.ch/mathlibs/sw/5_34_14/Minuit2/${NAME}.tar.gz
-tar -xzvf ${NAME}.tar.gz
-mkdir -p ${NAME}/build
-cd ${NAME}/build
-../configure --prefix=${PREFIX} --disable-openmp
-make -j4 
+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,11 +2,12 @@
 
 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 -P)
@@ -20,7 +21,7 @@ NAME=nlopt-${NAME}
 mkdir -p ${NAME}/build
 cd ${NAME}/build
 cmake -DCMAKE_INSTALL_PREFIX=${PREFIX} -DCMAKE_BUILD_WITH_INSTALL_NAME_DIR=TRUE -DCMAKE_INSTALL_NAME_DIR="${PREFIX}/lib" ..
-make -j4 
+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.1],[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/exRand.cpp

@@ -7,7 +7,7 @@
 using namespace std;
 using namespace Latan;
 
-constexpr Index  size          = 8;
+constexpr Index  n             = 8;
 constexpr Index  nDraw         = 20000;
 constexpr Index  nSample       = 2000;
 const     string stateFileName = "exRand.seed";
@@ -40,14 +40,14 @@ int main(void)
     p << PlotFunction(compile("return exp(-x_0^2/2)/sqrt(2*pi);", 1), -5., 5.);
     p.display();
 
-    DMat       var(size, size);
-    DVec       mean(size);
-    DMatSample sample(nSample, size, 1);
+    DMat       var(n, n);
+    DVec       mean(n);
+    DMatSample sample(nSample, n, 1);
 
     cout << "-- generating " << nSample << " Gaussian random vectors..." << endl;
-    var   = DMat::Random(size, size);
+    var   = DMat::Random(n, n);
     var  *= var.adjoint();
-    mean  = DVec::Random(size); 
+    mean  = DVec::Random(n); 
     RandomNormal mgauss(mean, var, rd());
     sample[central] = mgauss();
     FOR_STAT_ARRAY(sample, s)

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/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/Math.cpp

@@ -18,6 +18,7 @@
  */
 
 #include <LatAnalyze/Core/Math.hpp>
+#include <LatAnalyze/Numerical/GslFFT.hpp>
 #include <LatAnalyze/includes.hpp>
 #include <gsl/gsl_cdf.h>
 
@@ -29,7 +30,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 +39,54 @@ 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::conditionNumber(const DMat &mat)
+{
+    DVec s = mat.singularValues();
+
+    return s.maxCoeff()/s.minCoeff();
+}
+
+double MATH_NAMESPACE::cdr(const DMat &mat)
+{
+    return 10.*log10(conditionNumber(mat));
+}
+
+template <typename FFT>
+double nsdr(const DMat &m)
+{
+    Index  n = m.rows();
+    FFT    fft(n);
+    CMat   buf(n, 1);
+
+    FOR_VEC(buf, i)
+    {
+        buf(i) = 0.;
+        for (Index j = 0; j < n; ++j)
+        {
+            buf(i) += m(j, (i+j) % n);
+        }
+        buf(i) /= n;
+    }
+    fft(buf, FFT::Forward);
+
+    return 10.*log10(buf.real().maxCoeff()/buf.real().minCoeff());
+}
+
+double MATH_NAMESPACE::nsdr(const DMat &mat)
+{
+    return ::nsdr<GslFFT>(mat);
+}
+
 /******************************************************************************
  *                             Standard C functions                           *
  ******************************************************************************/

lib/Core/Math.hpp

@@ -70,6 +70,12 @@ 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 conditionNumber(const DMat &mat);
+    double cdr(const DMat &mat);
+    double nsdr(const DMat &mat);
     
     // Constants
     constexpr double pi  = 3.1415926535897932384626433832795028841970;

lib/Core/Plot.cpp

@@ -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,7 +332,14 @@ 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 ////////////////////////////////////////////////////
@@ -406,14 +515,16 @@ void Dash::operator()(PlotOptions &option) const
 }
 
 // LogScale constructor ////////////////////////////////////////////////////////
-LogScale::LogScale(const Axis axis)
+LogScale::LogScale(const Axis axis, const double basis)
 : axis_(axis)
+, basis_(basis)
 {}
 
 // Logscale modifier ///////////////////////////////////////////////////////////
 void LogScale::operator()(PlotOptions &option) const
 {
-    option.scaleMode[static_cast<int>(axis_)] |= Plot::Scale::log;
+    option.scaleMode[static_cast<int>(axis_)]     |= Plot::Scale::log;
+    option.logScaleBasis[static_cast<int>(axis_)]  = basis_;
 }
 
 // PlotRange constructors //////////////////////////////////////////////////////
@@ -806,11 +917,11 @@ ostream & Latan::operator<<(ostream &out, const Plot &plot)
     out << "unset log" << endl;
     if (plot.options_.scaleMode[x] & Plot::Scale::log)
     {
-        out << "set log x" << endl;
+        out << "set log x " << plot.options_.logScaleBasis[x] << endl;;
     }
     if (plot.options_.scaleMode[y] & Plot::Scale::log)
     {
-        out << "set log y" << endl;
+        out << "set log y " << plot.options_.logScaleBasis[y] << endl;
     }
     if (!plot.options_.label[x].empty())
     {

lib/Core/Plot.hpp

@@ -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                                 *
  ******************************************************************************/
@@ -199,6 +227,7 @@ struct PlotOptions
     std::string              caption;
     std::string              title;
     unsigned int             scaleMode[2];
+    double                   logScaleBasis[2];
     Range                    scale[2];
     std::string              label[2];
     std::string              lineColor;
@@ -286,13 +315,14 @@ class LogScale: public PlotModifier
 {
 public:
     // constructor
-    explicit LogScale(const Axis axis);
+    explicit LogScale(const Axis axis, const double basis = 10);
     // destructor
     virtual ~LogScale(void) = default;
     // modifier
     virtual void operator()(PlotOptions &option) const;
 private:
     const Axis axis_;
+    const double basis_;
 };
 
 class PlotRange: public PlotModifier

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/stdincludes.hpp

@@ -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/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       \
@@ -55,6 +58,7 @@ libLatAnalyze_la_SOURCES =           \
     Numerical/RootFinder.cpp         \
     Numerical/Solver.cpp             \
     Physics/CorrelatorFitter.cpp     \
+    Physics/DataFilter.cpp           \
     Physics/EffectiveMass.cpp        \
     Statistics/FitInterface.cpp      \
     Statistics/Histogram.cpp         \
@@ -75,6 +79,7 @@ HPPFILES =                           \
     Core/OptParser.hpp               \
     Core/ParserState.hpp             \
     Core/Plot.hpp                    \
+    Core/ThreadPool.hpp              \
     Core/stdincludes.hpp             \
     Core/Utilities.hpp               \
     Functional/CompiledFunction.hpp  \
@@ -90,6 +95,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\
@@ -100,6 +107,7 @@ HPPFILES =                           \
     Numerical/RootFinder.hpp         \
     Numerical/Solver.hpp             \
     Physics/CorrelatorFitter.hpp     \
+    Physics/DataFilter.hpp           \
     Physics/EffectiveMass.hpp        \
     Statistics/Dataset.hpp           \
     Statistics/FitInterface.hpp      \

lib/Numerical/DWT.cpp

@@ -0,0 +1,205 @@
+/*
+ * 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 ///////////////////////////////////////////////////////
+template <typename MatType>
+void filterConvolution(MatType &out, const MatType &data, 
+                       const std::vector<double> &filter, const Index offset)
+{
+    Index n = data.rows(), nf = n*filter.size();
+
+    out.resizeLike(data);
+    out.fill(0.);
+    for (unsigned int i = 0; i < filter.size(); ++i)
+    {
+        FOR_MAT(out, j, k)
+        {
+            out(j, k) += filter[i]*data((j + i + nf - offset) % n, k);
+        }
+    }
+}
+
+void DWT::filterConvolution(DVec &out, const DVec &data, 
+                            const std::vector<double> &filter, const Index offset)
+{
+    ::filterConvolution(out, data, filter, offset);
+}
+
+void DWT::filterConvolution(DMat &out, const DMat &data, 
+                            const std::vector<double> &filter, const Index offset)
+{
+    ::filterConvolution(out, data, filter, offset);
+}
+
+// downsampling/upsampling primitives //////////////////////////////////////////
+template <typename MatType>
+void downsample(MatType &out, const MatType &in)
+{
+    if (out.rows() < in.rows()/2)
+    {
+        LATAN_ERROR(Size, "output rows smaller than half the input vector rows");
+    }
+    if (out.cols() != in.cols())
+    {
+        LATAN_ERROR(Size, "output and input number of columns mismatch");
+    }
+    for (Index j = 0; j < in.cols(); j++)
+    for (Index i = 0; i < in.rows(); i += 2)
+    {
+        out(i/2, j) = in(i, j);
+    }
+}
+
+void DWT::downsample(DVec &out, const DVec &in)
+{
+    ::downsample(out, in);
+}
+
+void DWT::downsample(DMat &out, const DMat &in)
+{
+    ::downsample(out, in);
+}
+
+template <typename MatType>
+void upsample(MatType &out, const MatType &in)
+{
+    if (out.size() < 2*in.size())
+    {
+        LATAN_ERROR(Size, "output rows smaller than twice the input rows");
+    }
+    if (out.cols() != in.cols())
+    {
+        LATAN_ERROR(Size, "output and input number of columns mismatch");
+    }
+    out.block(0, 0, 2*in.size(), out.cols()).fill(0.);
+    for (Index j = 0; j < in.cols(); j++)
+    for (Index i = 0; i < in.size(); i ++)
+    {
+        out(2*i, j) = in(i, j);
+    }
+}
+
+void DWT::upsample(DVec &out, const DVec &in)
+{
+    ::upsample(out, in);
+}
+
+void DWT::upsample(DMat &out, const DMat &in)
+{
+    ::upsample(out, in);
+}
+
+// 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;
+}
+
+// concatenate levels //////////////////////////////////////////////////////////
+DVec DWT::concat(const std::vector<DWTLevel> &dwt, const int maxLevel, const bool dropLow)
+{
+    unsigned int level = ((maxLevel >= 0) ? (maxLevel + 1) : dwt.size());
+    Index nlast = dwt[level - 1].first.size();
+    Index n = 2*dwt.front().first.size() - ((dropLow) ? nlast : 0);
+    Index pt = n, nl;
+    DVec  res(n);
+
+    for (unsigned int l = 0; l < level; ++l)
+    {
+        nl  = dwt[l].second.size();
+        pt -= nl;
+        res.segment(pt, nl) = dwt[l].second;
+    }
+    if (!dropLow)
+    {
+        res.segment(0, nl) = dwt[level-1].first;
+    }
+
+    return res;
+}

lib/Numerical/DWT.hpp

@@ -0,0 +1,62 @@
+/*
+ * 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>
+#include <LatAnalyze/Core/Mat.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);
+    static void filterConvolution(DMat &out, const DMat &data, 
+                                  const std::vector<double> &filter, const Index offset);
+    // downsampling/upsampling primitives
+    static void downsample(DVec &out, const DVec &in);
+    static void downsample(DMat &out, const DMat &in);
+    static void upsample(DVec &out, const DVec &in);
+    static void upsample(DMat &out, const DMat &in);
+    // DWT
+    std::vector<DWTLevel> forward(const DVec &data, const unsigned int level) const;
+    DVec                  backward(const std::vector<DWTLevel>& dwt) const;
+    // concatenate levels
+    static DVec concat(const std::vector<DWTLevel>& dwt, const int maxLevel = -1, const bool dropLow = false);
+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/MinuitMinimizer.cpp

@@ -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/Physics/DataFilter.cpp

@@ -0,0 +1,83 @@
+/*
+ * DataFilter.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/DataFilter.hpp>
+#include <LatAnalyze/includes.hpp>
+#include <LatAnalyze/Numerical/DWT.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                        DataFilter implementation                           *
+ ******************************************************************************/
+// constructor ////////////////////////////////////////////////////////////////
+DataFilter::DataFilter(const vector<double> &filter, const bool downsample)
+: filter_(filter), downsample_(downsample)
+{}
+
+// filtering //////////////////////////////////////////////////////////////////
+template <typename MatType>
+void filter(MatType &out, const MatType &in, const vector<double> &filter, 
+            const bool downsample, MatType &buf)
+{
+    if (!downsample)
+    {
+        out.resizeLike(in);
+        DWT::filterConvolution(out, in, filter, filter.size()/2);
+    }
+    else
+    {
+        out.resize(in.rows()/2, in.cols());
+        buf.resizeLike(in);
+        DWT::filterConvolution(buf, in, filter, filter.size()/2);
+        DWT::downsample(out, buf);
+    }
+}
+
+void DataFilter::operator()(DVec &out, const DVec &in)
+{
+    filter(out, in, filter_, downsample_, vBuf_);
+}
+
+void DataFilter::operator()(DMat &out, const DMat &in)
+{
+    filter(out, in, filter_, downsample_, mBuf_);
+}
+
+/******************************************************************************
+ *                     LaplaceDataFilter implementation                       *
+ ******************************************************************************/
+// constructor ////////////////////////////////////////////////////////////////
+LaplaceDataFilter::LaplaceDataFilter(const bool downsample)
+: DataFilter({1., -2. , 1.}, downsample)
+{}
+
+// filtering //////////////////////////////////////////////////////////////////
+void LaplaceDataFilter::operator()(DVec &out, const DVec &in, const double lambda)
+{
+    filter_[1] = -2. - lambda;
+    DataFilter::operator()(out, in);
+}
+
+void LaplaceDataFilter::operator()(DMat &out, const DMat &in, const double lambda)
+{
+    filter_[1] = -2. - lambda;
+    DataFilter::operator()(out, in);
+}

lib/Physics/DataFilter.hpp

@@ -0,0 +1,139 @@
+/*
+ * DataFilter.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_DataFilter_hpp_
+#define Latan_DataFilter_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Core/Math.hpp>
+#include <LatAnalyze/Statistics/StatArray.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
+#include <LatAnalyze/Numerical/Minimizer.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                    Generic convolution filter class                        *
+ ******************************************************************************/
+class DataFilter
+{
+public:
+    // constructor
+    DataFilter(const std::vector<double> &filter, const bool downsample = false);
+    // filtering
+    void operator()(DVec &out, const DVec &in);
+    void operator()(DMat &out, const DMat &in);
+    template <typename MatType, Index o>
+    void operator()(StatArray<MatType, o> &out, const StatArray<MatType, o> &in);
+protected:
+    std::vector<double> filter_;
+private:
+    bool                downsample_;
+    DVec                vBuf_;
+    DMat                mBuf_;
+};
+
+/******************************************************************************
+ *                         Laplacian filter class                             *
+ ******************************************************************************/
+class LaplaceDataFilter: public DataFilter
+{
+public:
+    // constructor
+    LaplaceDataFilter(const bool downsample = false);
+    // filtering
+    void operator()(DVec &out, const DVec &in, const double lambda = 0.);
+    void operator()(DMat &out, const DMat &in, const double lambda = 0.);
+    template <typename MatType, Index o>
+    void operator()(StatArray<MatType, o> &out, const StatArray<MatType, o> &in, 
+                    const double lambda = 0.);
+    // correlation optimisation
+    template <typename MatType, Index o>
+    double optimiseCdr(const StatArray<MatType, o> &data, Minimizer &min,
+                       const unsigned int nPass = 3);
+};
+
+/******************************************************************************
+ *                 DataFilter class template implementation                   *
+ ******************************************************************************/
+// filtering //////////////////////////////////////////////////////////////////
+template <typename MatType, Index o>
+void DataFilter::operator()(StatArray<MatType, o> &out, const StatArray<MatType, o> &in)
+{
+    FOR_STAT_ARRAY(in, s)
+    {
+        (*this)(out[s], in[s]);
+    }
+}
+
+/******************************************************************************
+ *            LaplaceDataFilter class template implementation                 *
+ ******************************************************************************/
+// filtering //////////////////////////////////////////////////////////////////
+template <typename MatType, Index o>
+void LaplaceDataFilter::operator()(StatArray<MatType, o> &out, 
+                                   const StatArray<MatType, o> &in, const double lambda)
+{
+    FOR_STAT_ARRAY(in, s)
+    {
+        (*this)(out[s], in[s], lambda);
+    }
+}
+
+// correlation optimisation ///////////////////////////////////////////////////
+template <typename MatType, Index o>
+double LaplaceDataFilter::optimiseCdr(const StatArray<MatType, o> &data, 
+                                      Minimizer &min, const unsigned int nPass)
+{
+    StatArray<MatType, o> fdata(data.size());
+    DVec init(1);
+    double reg, prec;
+    DoubleFunction cdr([&data, &fdata, this](const double *x)
+    {
+        double res;
+        (*this)(fdata, data, x[0]);
+        res = Math::cdr(fdata.correlationMatrix());
+        
+        return res;
+    }, 1);
+
+    min.setLowLimit(0., -0.1);
+    min.setHighLimit(0., 100000.);
+    init(0) = 0.1;
+    min.setInit(init);
+    prec = 0.1;
+    min.setPrecision(prec);
+    reg = min(cdr)(0);
+    for (unsigned int pass = 0; pass < nPass; pass++)
+    {
+        min.setLowLimit(0., (1.-10.*prec)*reg);
+        min.setHighLimit(0., (1.+10.*prec)*reg);
+        init(0) = reg;
+        min.setInit(init);
+        prec *= 0.1;
+        min.setPrecision(prec);
+        reg = min(cdr)(0);
+    }
+
+    return reg;
+}
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_DataFilter_hpp_

lib/Statistics/MatSample.hpp

@@ -379,7 +379,6 @@ void MatSample<T>::resizeMat(const Index nRow, const Index nCol)
     }
 }
 
-
 END_LATAN_NAMESPACE
 
 #endif // Latan_MatSample_hpp_

lib/Statistics/StatArray.hpp

@@ -51,14 +51,14 @@ public:
     const T & operator[](const Index s) const;
     // statistics
     void bin(Index binSize);
+    T    sum(const Index pos = 0, const Index n = -1) const;
     T    mean(const Index pos = 0, const Index n = -1) const;
-    T    covariance(const StatArray<T, os> &array, const Index pos = 0,
-                    const Index n = -1) const;
-    T    covarianceMatrix(const StatArray<T, os> &array, const Index pos = 0,
-                          const Index n = -1) const;
-    T    variance(const Index pos = 0, const Index n = -1) const;
-    T    varianceMatrix(const Index pos = 0, const Index n = -1) const;
-    T    correlationMatrix(const Index pos = 0, const Index n = -1) const;
+    T    covariance(const StatArray<T, os> &array) const;
+    T    variance(void) const;
+    T    covarianceMatrix(const StatArray<T, os> &data) const;
+    T    varianceMatrix(void) const;
+    T    correlationMatrix(void) const;
+
     // IO type
     virtual IoType getType(void) const;
 public:
@@ -66,7 +66,7 @@ public:
 };
 
 // reduction operations
-namespace ReducOp
+namespace StatOp
 {
     // general templates
     template <typename T>
@@ -148,80 +148,118 @@ void StatArray<T, os>::bin(Index binSize)
     }
 }
 
+template <typename T, Index os>
+T StatArray<T, os>::sum(const Index pos, const Index n) const
+{
+    T           result;
+    const Index m = (n >= 0) ? n : size();
+
+    result = (*this)[pos];
+    for (Index i = pos + 1; i < pos + m; ++i)
+    {
+        result += (*this)[i];
+    }
+
+    return result;
+}
+
 template <typename T, Index os>
 T StatArray<T, os>::mean(const Index pos, const Index n) const
 {
-    T           result = T();
     const Index m = (n >= 0) ? n : size();
 
-    if (m)
+    return sum(pos, n)/static_cast<double>(m);
+}
+
+template <typename T, Index os>
+T StatArray<T, os>::covariance(const StatArray<T, os> &array) const
+{
+    T s1, s2, res;
+
+    s1  = array.sum();
+    s2  = this->sum();
+    res = StatOp::prod<T>(array[0], (*this)[0]);
+    for (Index i = 1; i < size(); ++i)
     {
-        result = this->segment(pos+os, m).redux(&ReducOp::sum<T>);
+        res += StatOp::prod<T>(array[i], (*this)[i]);
     }
-    return result/static_cast<double>(m);
+    res -= StatOp::prod<T>(s1, s2)/static_cast<double>(size());
+    res /= static_cast<double>(size() - 1);
+
+    return res;
 }
 
 template <typename T, Index os>
-T StatArray<T, os>::covariance(const StatArray<T, os> &array, const Index pos,
-                               const Index n) const
+T StatArray<T, os>::variance(void) const
 {
-    T           s1, s2, prs, res = T();
-    const Index m = (n >= 0) ? n : size();
-    
-    if (m)
+    return covariance(*this);
+}
+
+template <typename MatType, Index os>
+MatType StatArray<MatType, os>::covarianceMatrix(
+    const StatArray<MatType, os> &data) const
+{
+    if (((*this)[central].cols() != 1) or (data[central].cols() != 1))
     {
-        auto arraySeg = array.segment(pos+os, m);
-        auto thisSeg  = this->segment(pos+os, m);
-        
-        s1  = thisSeg.redux(&ReducOp::sum<T>);
-        s2  = arraySeg.redux(&ReducOp::sum<T>);
-        prs = thisSeg.binaryExpr(arraySeg, &ReducOp::prod<T>)
-                     .redux(&ReducOp::sum<T>);
-        res = prs - ReducOp::prod(s1, s2)/static_cast<double>(m);
+        LATAN_ERROR(Size, "samples have more than one column");
     }
-    
-    return res/static_cast<double>(m - 1);
-}
 
-template <typename T, Index os>
-T StatArray<T, os>::covarianceMatrix(const StatArray<T, os> &array,
-                                     const Index pos, const Index n) const
-{
-    T           s1, s2, prs, res = T();
-    const Index m = (n >= 0) ? n : size();
-    
-    if (m)
+    Index   n1 = (*this)[central].rows(), n2 = data[central].rows();
+    Index   nSample = this->size();
+    MatType tmp1(n1, nSample), tmp2(n2, nSample), res(n1, n2);
+    MatType 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)
     {
-        auto arraySeg = array.segment(pos+os, m);
-        auto thisSeg  = this->segment(pos+os, m);
-
-        s1  = thisSeg.redux(&ReducOp::sum<T>);
-        s2  = arraySeg.redux(&ReducOp::sum<T>);
-        prs = thisSeg.binaryExpr(arraySeg, &ReducOp::tensProd<T>)
-                     .redux(&ReducOp::sum<T>);
-        res = prs - ReducOp::tensProd(s1, s2)/static_cast<double>(m);
+        s1          += (*this)[s];
+        tmp1.col(s)  = (*this)[s];
     }
-    
-    return res/static_cast<double>(m - 1);
+    tmp1 -= s1*one.transpose()/static_cast<double>(nSample);
+    for (unsigned int s = 0; s < nSample; ++s)
+    {
+        s2          += data[s];
+        tmp2.col(s)  = data[s];
+    }
+    tmp2 -= s2*one.transpose()/static_cast<double>(nSample);
+    res   = tmp1*tmp2.transpose()/static_cast<double>(nSample - 1);
+
+    return res;
 }
 
-template <typename T, Index os>
-T StatArray<T, os>::variance(const Index pos, const Index n) const
+template <typename MatType, Index os>
+MatType StatArray<MatType, os>::varianceMatrix(void) const
 {
-    return covariance(*this, pos, n);
+    if ((*this)[0].cols() != 1)
+    {
+        LATAN_ERROR(Size, "samples have more than one column");
+    }
+
+    Index   n1 = (*this)[0].rows();
+    Index   nSample = this->size();
+    MatType tmp1(n1, nSample), res(n1, n1);
+    MatType 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, Index os>
-T StatArray<T, os>::varianceMatrix(const Index pos, const Index n) const
+template <typename MatType, Index os>
+MatType StatArray<MatType, os>::correlationMatrix(void) const
 {
-    return covarianceMatrix(*this, pos, n);
-}
-
-template <typename T, Index os>
-T StatArray<T, os>::correlationMatrix(const Index pos, const Index n) const
-{
-    T res = varianceMatrix(pos, n);
-    T invDiag(res.rows(), 1);
+    MatType res = varianceMatrix();
+    MatType invDiag(res.rows(), 1);
 
     invDiag = res.diagonal();
     invDiag = invDiag.cwiseInverse().cwiseSqrt();
@@ -231,45 +269,19 @@ T StatArray<T, os>::correlationMatrix(const Index pos, const Index n) const
 }
 
 // reduction operations ////////////////////////////////////////////////////////
-namespace ReducOp
+namespace StatOp
 {
-    template <typename T>
-    inline T sum(const T &a, const T &b)
-    {
-        return a + b;
-    }
-
     template <typename T>
     inline T prod(const T &a, const T &b)
     {
         return a*b;
     }
 
-    template <typename T>
-    inline T tensProd(const T &v1 __dumb, const T &v2 __dumb)
-    {
-        LATAN_ERROR(Implementation,
-                    "tensorial product not implemented for this type");
-    }
-
     template <>
     inline Mat<double> prod(const Mat<double>  &a, const Mat<double>  &b)
     {
         return a.cwiseProduct(b);
     }
-
-    template <>
-    inline Mat<double> tensProd(const Mat<double>  &v1,
-                                const Mat<double>  &v2)
-    {
-        if ((v1.cols() != 1) or (v2.cols() != 1))
-        {
-            LATAN_ERROR(Size,
-                        "tensorial product is only valid with column vectors");
-        }
-        
-        return v1*v2.transpose();
-    }
 }
 
 // IO type /////////////////////////////////////////////////////////////////////

lib/Statistics/XYSampleData.cpp

@@ -62,6 +62,11 @@ double SampleFitResult::getPValue(const Index s) const
     return Math::chi2PValue(getChi2(s), getNDof());
 }
 
+double SampleFitResult::getCorrRangeDb(void) const
+{
+    return corrRangeDb_;
+}
+
 double SampleFitResult::getCcdf(const Index s) const
 {
     return Math::chi2Ccdf(getChi2(s), getNDof());
@@ -107,9 +112,11 @@ void SampleFitResult::print(const bool printXsi, ostream &out) const
         getChi2(), static_cast<int>(getNDof()), getChi2PerDof(), getCcdf(), 
         getPValue());
     out << buf << endl;
+    sprintf(buf, "correlation dynamic range= %.1f dB", getCorrRangeDb());
+    out << buf << endl;
     for (Index p = 0; p < pMax; ++p)
     {
-        sprintf(buf, "%8s= % e +/- %e", parName_[p].c_str(),
+        sprintf(buf, "%12s= % e +/- %e", parName_[p].c_str(),
                 (*this)[central](p), err(p));
         out << buf << endl;
     }
@@ -249,6 +256,20 @@ const DMat & XYSampleData::getFitVarMatPInv(void)
     return data_.getFitVarMatPInv();
 }
 
+const DMat & XYSampleData::getFitCorrMat(void)
+{
+    computeVarMat();
+    
+    return data_.getFitCorrMat();
+}
+
+const DMat & XYSampleData::getFitCorrMatPInv(void)
+{
+    computeVarMat();
+    
+    return data_.getFitCorrMatPInv();
+}
+
 // set data to a particular sample /////////////////////////////////////////////
 void XYSampleData::setDataToSample(const Index s)
 {
@@ -285,9 +306,10 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
 {
     computeVarMat();
     
-    SampleFitResult result;
-    FitResult       sampleResult;
-    DVec            initCopy = init;
+    SampleFitResult      result;
+    FitResult            sampleResult;
+    DVec                 initCopy = init;
+    Minimizer::Verbosity verbCopy = minimizer.back()->getVerbosity();
     
     result.resize(nSample_);
     result.chi2_.resize(nSample_);
@@ -299,9 +321,14 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
         {
             sampleResult = data_.fit(minimizer, initCopy, v);
             initCopy     = sampleResult.segment(0, initCopy.size());
+            if (verbCopy != Minimizer::Verbosity::Debug)
+            {
+                minimizer.back()->setVerbosity(Minimizer::Verbosity::Silent);
+            }
         }
         else
         {
+            
             sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
         }
         result[s]       = sampleResult;
@@ -312,9 +339,11 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
             result.model_[j][s] = sampleResult.getModel(j);
         }
     }
-    result.nPar_    = sampleResult.getNPar();
-    result.nDof_    = sampleResult.nDof_;
-    result.parName_ = sampleResult.parName_;
+    minimizer.back()->setVerbosity(verbCopy);
+    result.nPar_       = sampleResult.getNPar();
+    result.nDof_       = sampleResult.nDof_;
+    result.parName_    = sampleResult.parName_;
+    result.corrRangeDb_ = Math::cdr(getFitCorrMat());
     
     return result;
 }
@@ -346,6 +375,29 @@ XYSampleData XYSampleData::getResiduals(const SampleFitResult &fit)
     return res;
 }
 
+XYSampleData XYSampleData::getNormalisedResiduals(const SampleFitResult &fit)
+{
+    XYSampleData res(*this);
+    
+    for (Index j = 0; j < getNYDim(); ++j)
+    {
+        const DoubleFunctionSample &f   = fit.getModel(_, j);
+        
+        for (auto &p: yData_[j])
+        {
+            res.y(p.first, j) -= f(x(p.first));
+        }
+
+        const DMat &var = res.getYYVar(j, j);
+        for (auto &p: yData_[j])
+        {
+            res.y(p.first, j) /= sqrt(var(p.first, p.first));
+        }
+    }
+    
+    return res;
+}
+
 XYSampleData XYSampleData::getPartialResiduals(const SampleFitResult &fit,
                                                const DVec &ref, const Index i)
 {

lib/Statistics/XYSampleData.hpp

@@ -49,6 +49,7 @@ public:
     double                       getNDof(void) const;
     Index                        getNPar(void) const;
     double                       getPValue(const Index s = central) const;
+    double                       getCorrRangeDb(void) const;
     double                       getCcdf(const Index s = central) const;
     const DoubleFunction &       getModel(const Index s = central,
                                           const Index j = 0) const;
@@ -60,6 +61,7 @@ public:
                std::ostream &out = std::cout) const;
 private:
     DSample                           chi2_;
+    double                            corrRangeDb_{0.};
     Index                             nDof_{0}, nPar_{0};
     std::vector<DoubleFunctionSample> model_;
     std::vector<std::string>          parName_;
@@ -91,9 +93,11 @@ public:
     const DMat &       getXYVar(const Index i, const Index j);
     DVec               getXError(const Index i);
     DVec               getYError(const Index j);
-    // get total fit variance matrix and its pseudo-inverse
+    // get total fit variance & correlation matrices and their pseudo-inverse
     const DMat & getFitVarMat(void);
     const DMat & getFitVarMatPInv(void);
+    const DMat & getFitCorrMat(void);
+    const DMat & getFitCorrMatPInv(void);
     // set data to a particular sample
     void setDataToSample(const Index s);
     // get internal XYStatData
@@ -111,6 +115,7 @@ public:
                         const DoubleModel &model, const Ts... models);
     // residuals
     XYSampleData getResiduals(const SampleFitResult &fit);
+    XYSampleData getNormalisedResiduals(const SampleFitResult &fit);
     XYSampleData getPartialResiduals(const SampleFitResult &fit, const DVec &x,
                                      const Index i);
 private:

lib/Statistics/XYStatData.cpp

@@ -60,6 +60,11 @@ double FitResult::getCcdf(void) const
     return Math::chi2Ccdf(getChi2(), getNDof());;
 }
 
+double FitResult::getCorrRangeDb(void) const
+{
+    return corrRangeDb_;
+}
+
 const DoubleFunction & FitResult::getModel(const Index j) const
 {
     return model_[j];
@@ -75,9 +80,11 @@ void FitResult::print(const bool printXsi, ostream &out) const
         getChi2(), static_cast<int>(getNDof()), getChi2PerDof(), getCcdf(), 
         getPValue());
     out << buf << endl;
+    sprintf(buf, "correlation dynamic range= %.1f dB", getCorrRangeDb());
+    out << buf << endl;
     for (Index p = 0; p < pMax; ++p)
     {
-        sprintf(buf, "%8s= %e", parName_[p].c_str(), (*this)(p));
+        sprintf(buf, "%12s= %e", parName_[p].c_str(), (*this)(p));
         out << buf << endl;
     }
 }
@@ -216,7 +223,7 @@ DVec XYStatData::getXError(const Index i) const
 
 DVec XYStatData::getYError(const Index j) const
 {
-    checkXDim(j);
+    checkYDim(j);
     
     return yyVar_(j, j).diagonal().cwiseSqrt();
 }
@@ -259,6 +266,20 @@ const DMat & XYStatData::getFitVarMatPInv(void)
     return fitVarInv_;
 }
 
+const DMat & XYStatData::getFitCorrMat(void)
+{
+    updateFitVarMat();
+    
+    return fitCorr_;
+}
+
+const DMat & XYStatData::getFitCorrMatPInv(void)
+{
+    updateFitVarMat();
+    
+    return fitCorrInv_;
+}
+
 // fit /////////////////////////////////////////////////////////////////////////
 FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
                           const vector<const DoubleModel *> &v)
@@ -337,9 +358,10 @@ FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
         result    = (*m)(chi2);
         totalInit = result;
     }
-    result.chi2_ = chi2(result);
-    result.nPar_ = nPar;
-    result.nDof_ = layout.totalYSize - nPar;
+    result.corrRangeDb_ = Math::cdr(getFitCorrMat());
+    result.chi2_       = chi2(result);
+    result.nPar_       = nPar;
+    result.nDof_       = layout.totalYSize - nPar;
     result.model_.resize(v.size());
     for (unsigned int j = 0; j < v.size(); ++j)
     {
@@ -379,6 +401,27 @@ XYStatData XYStatData::getResiduals(const FitResult &fit)
     return res;
 }
 
+XYStatData XYStatData::getNormalisedResiduals(const FitResult &fit)
+{
+    XYStatData res(*this);
+    
+    for (Index j = 0; j < getNYDim(); ++j)
+    {
+        const DoubleFunction &f  = fit.getModel(j);
+        const DVec           err = getYError(j);
+        Index                row = 0;
+        
+        for (auto &p: yData_[j])
+        {
+            res.y(p.first, j) -= f(x(p.first));
+            res.y(p.first, j) /= err(row);
+            row++;
+        }
+    }
+    
+    return res;
+}
+
 XYStatData XYStatData::getPartialResiduals(const FitResult &fit,
                                            const DVec &ref, const Index i)
 {
@@ -530,8 +573,11 @@ void XYStatData::updateFitVarMat(void)
         chi2DataVec_.resize(layout.totalSize);
         chi2ModVec_.resize(layout.totalSize);
         chi2Vec_.resize(layout.totalSize);
-        fitVar_    = fitVar_.cwiseProduct(makeCorrFilter());
-        fitVarInv_ = fitVar_.pInverse(getSvdTolerance());
+        fitVar_     = fitVar_.cwiseProduct(makeCorrFilter());
+        fitCorr_    = Math::varToCorr(fitVar_);
+        fitCorrInv_ = fitCorr_.pInverse(getSvdTolerance());
+        fitVarInv_  = Math::corrToVar(fitCorrInv_, fitVar_.diagonal().cwiseInverse());
+
         scheduleFitVarMatInit(false);
     }
 }

lib/Statistics/XYStatData.hpp

@@ -48,12 +48,13 @@ public:
     Index                  getNPar(void) const;
     double                 getPValue(void) const;
     double                 getCcdf(void) const;
+    double                 getCorrRangeDb(void) const;
     const DoubleFunction & getModel(const Index j = 0) const;
     // IO
     void print(const bool printXsi = false,
                std::ostream &out = std::cout) const;
 private:
-    double                      chi2_{0.};
+    double                      chi2_{0.}, corrRangeDb_{0.};
     Index                       nDof_{0}, nPar_{0};
     std::vector<DoubleFunction> model_;
     std::vector<std::string>    parName_;
@@ -88,9 +89,11 @@ public:
     DVec           getXError(const Index i) const;
     DVec           getYError(const Index j) const;
     DMat           getTable(const Index i, const Index j) const;
-    // get total fit variance matrix and its pseudo-inverse
+    // get total fit variance & correlation matrices and their pseudo-inverse
     const DMat & getFitVarMat(void);
     const DMat & getFitVarMatPInv(void);
+    const DMat & getFitCorrMat(void);
+    const DMat & getFitCorrMatPInv(void);
     // fit
     FitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
                   const std::vector<const DoubleModel *> &v);
@@ -104,6 +107,7 @@ public:
                   const DoubleModel &model, const Ts... models);
     // residuals
     XYStatData getResiduals(const FitResult &fit);
+    XYStatData getNormalisedResiduals(const FitResult &fit);
     XYStatData getPartialResiduals(const FitResult &fit, const DVec &ref,
                                    const Index i);
 protected:
@@ -130,7 +134,7 @@ private:
     std::vector<DVec>                    xData_;
     std::vector<DVec>                    xMap_;
     Mat<DMat>                            xxVar_, yyVar_, xyVar_;
-    DMat                                 fitVar_, fitVarInv_;
+    DMat                                 fitVar_, fitVarInv_, fitCorr_, fitCorrInv_;
     DVec                                 chi2DataVec_, chi2ModVec_, chi2Vec_;
     DVec                                 xBuf_;
     bool                                 initXMap_{true};

physics/2pt-fit.cpp

@@ -3,6 +3,7 @@
 #include <LatAnalyze/Core/Plot.hpp>
 #include <LatAnalyze/Functional/CompiledModel.hpp>
 #include <LatAnalyze/Io/Io.hpp>
+#include <LatAnalyze/Numerical/DWT.hpp>
 #include <LatAnalyze/Numerical/MinuitMinimizer.hpp>
 #include <LatAnalyze/Numerical/NloptMinimizer.hpp>
 #include <LatAnalyze/Physics/CorrelatorFitter.hpp>
@@ -23,7 +24,7 @@ int main(int argc, char *argv[])
 {
     // parse arguments /////////////////////////////////////////////////////////
     OptParser            opt;
-    bool                 parsed, doPlot, doHeatmap, doCorr, fold, doScan;
+    bool                 parsed, doLaplace, doPlot, doHeatmap, doCorr, fold, doScan, noGuess;
     string               corrFileName, model, outFileName, outFmt, savePlot;
     Index                ti, tf, shift, nPar, thinning;
     double               svdTol;
@@ -52,10 +53,14 @@ int main(int argc, char *argv[])
                   "only do the uncorrelated fit");
     opt.addOption("" , "fold"   , OptParser::OptType::trigger, true,
                   "fold the correlator");
+    opt.addOption("l" , "laplace"   , OptParser::OptType::trigger, true,
+                  "apply Laplace filter to the correlator");
     opt.addOption("p", "plot"     , OptParser::OptType::trigger, true,
                   "show the fit plot");
     opt.addOption("h", "heatmap"  , OptParser::OptType::trigger, true,
                   "show the fit correlation heatmap");
+    opt.addOption("", "no-guess"  , OptParser::OptType::trigger, true,
+                  "do not try to guess fit parameters");
     opt.addOption("", "save-plot", OptParser::OptType::value, true,
                     "saves the source and .pdf", "");
     opt.addOption("", "scan", OptParser::OptType::trigger, true,
@@ -74,6 +79,7 @@ int main(int argc, char *argv[])
     ti           = opt.optionValue<Index>("ti");
     tf           = opt.optionValue<Index>("tf");
     thinning     = opt.optionValue<Index>("t");
+    doLaplace    = opt.gotOption("l");
     shift        = opt.optionValue<Index>("s");
     model        = opt.optionValue("m");
     nPar         = opt.optionValue<Index>("nPar");
@@ -83,6 +89,7 @@ int main(int argc, char *argv[])
     fold         = opt.gotOption("fold");
     doPlot       = opt.gotOption("p");
     doHeatmap    = opt.gotOption("h");
+    noGuess      = opt.gotOption("no-guess");
     savePlot     = opt.optionValue("save-plot");
     doScan       = opt.gotOption("scan");
     switch (opt.optionValue<unsigned int>("v"))
@@ -110,6 +117,17 @@ int main(int argc, char *argv[])
     nt      = corr[central].rows();
     corr    = corr.block(0, 0, nt, 1);
     corr    = CorrelatorUtils::shift(corr, shift);
+    if (doLaplace)
+    {
+        vector<double> filter = {1., -2., 1.};
+        DVec           buf;
+
+        FOR_STAT_ARRAY(corr, s)
+        {
+            DWT::filterConvolution(buf, corr[s], filter, 1);
+            corr[s] = buf;
+        }
+    }
     if (fold)
     {
         corr = CorrelatorUtils::fold(corr);
@@ -152,13 +170,14 @@ int main(int argc, char *argv[])
     fitter.setThinning(thinning);
 
     // set initial values ******************************************************
-    if (modelPar.type != CorrelatorType::undefined)
+    if ((modelPar.type != CorrelatorType::undefined) and !noGuess)
     {
         init = CorrelatorModels::parameterGuess(corr, modelPar);
     }
     else
     {
-        init.fill(0.1);
+        init.fill(1.);
+        init(0) = 0.2;
     }
 
     // set limits for minimisers ***********************************************
@@ -282,7 +301,7 @@ int main(int argc, char *argv[])
                 DMat  id = DMat::Identity(n, n),
                       var = fitter.data().getFitVarMat();
                 
-                p << PlotMatrix(Math::varToCorr(var));
+                p << PlotCorrMatrix(Math::varToCorr(var));
                 p << Caption("correlation matrix");
                 p.display();
                 if (svdTol > 0.)

utils/Makefile.am

@@ -9,9 +9,11 @@ endif
 bin_PROGRAMS =            \
     latan-plot            \
     latan-sample-combine  \
+    latan-sample-dwt      \
     latan-sample-element  \
     latan-sample-fake     \
     latan-sample-ft       \
+    latan-sample-merge    \
     latan-sample-plot     \
     latan-sample-plot-corr\
     latan-sample-read     \
@@ -25,6 +27,10 @@ latan_sample_combine_SOURCES  = sample-combine.cpp
 latan_sample_combine_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_combine_LDFLAGS  = -L../lib/.libs -lLatAnalyze
 
+latan_sample_dwt_SOURCES  = sample-dwt.cpp
+latan_sample_dwt_CXXFLAGS = $(COM_CXXFLAGS)
+latan_sample_dwt_LDFLAGS  = -L../lib/.libs -lLatAnalyze
+
 latan_sample_element_SOURCES  = sample-element.cpp
 latan_sample_element_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_element_LDFLAGS  = -L../lib/.libs -lLatAnalyze
@@ -37,6 +43,10 @@ latan_sample_ft_SOURCES  = sample-ft.cpp
 latan_sample_ft_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_ft_LDFLAGS  = -L../lib/.libs -lLatAnalyze
 
+latan_sample_merge_SOURCES  = sample-merge.cpp
+latan_sample_merge_CXXFLAGS = $(COM_CXXFLAGS)
+latan_sample_merge_LDFLAGS  = -L../lib/.libs -lLatAnalyze
+
 latan_sample_plot_corr_SOURCES  = sample-plot-corr.cpp
 latan_sample_plot_corr_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_plot_corr_LDFLAGS  = -L../lib/.libs -lLatAnalyze

utils/sample-dwt.cpp

@@ -0,0 +1,199 @@
+/*
+ * sample-dwt.cpp, part of LatAnalyze 3
+ *
+ * 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
+ * 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/Math.hpp>
+#include <LatAnalyze/Core/OptParser.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
+#include <LatAnalyze/Io/Io.hpp>
+#include <LatAnalyze/Numerical/DWT.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char *argv[])
+{
+    // argument parsing ////////////////////////////////////////////////////////
+    OptParser    opt;
+    bool         parsed, doPlot, ss, saveAll;
+    unsigned int level;
+    string       inFilename, outFilename, filterName;
+    
+    opt.addOption("l", "level", OptParser::OptType::value, true,
+                  "number of levels", "1");
+    opt.addOption("f", "filter", OptParser::OptType::value, true,
+                  "filter name", "haar");
+    opt.addOption("s", "ss", OptParser::OptType::trigger, true,
+                  "super-sampling (inverse DWT on data)");
+    opt.addOption("a", "all", OptParser::OptType::trigger, true,
+                  "save all-levels");
+    opt.addOption("o", "output", OptParser::OptType::value, true,
+                  "output file name, or directory name if --all is used (default: result not saved)", "");
+    opt.addOption("p", "plot", OptParser::OptType::trigger, true,
+                  "show plot");
+    opt.addOption("" , "help"  , OptParser::OptType::trigger, true,
+                  "show this help message and exit");
+    parsed = opt.parse(argc, argv);
+    if (!parsed or (opt.getArgs().size() != 1) or opt.gotOption("help"))
+    {
+        cerr << "usage: " << argv[0];
+        cerr << " <options> <input file>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
+        cerr << "Available DWT filters:" << endl;
+        for (auto &fv: DWTFilters::fromName)
+        {
+            cerr << fv.first << " ";
+        }
+        cerr << endl << endl;
+
+        return EXIT_FAILURE;
+    }
+    inFilename  = opt.getArgs()[0];
+    level       = opt.optionValue<unsigned int>("l");
+    filterName  = opt.optionValue("f");
+    ss          = opt.gotOption("s");
+    saveAll     = opt.gotOption("a");
+    outFilename = opt.optionValue("o");
+    doPlot      = opt.gotOption("p");
+
+    // DWT /////////////////////////////////////////////////////////////////////
+    DMatSample            in = Io::load<DMatSample>(inFilename), res;
+    Index                 nSample = in.size(), n = in[central].rows();
+    vector<DMatSample>    out(ss ? 1 : level, DMatSample(nSample)), 
+                          outh(ss ? 0 : level, DMatSample(nSample)),
+                          concath(ss ? 0 : level, DMatSample(nSample));
+    DMatSample            concat(nSample, n, 1);
+    DWT                   dwt(*DWTFilters::fromName.at(filterName));
+    vector<DWT::DWTLevel> dataDWT(level);
+
+    FOR_STAT_ARRAY(in, s)
+    {
+        in[s].conservativeResize(n, 1);
+    }
+    if (!ss)
+    {
+        DMatSample buf(nSample);
+
+        cout << "-- compute discrete wavelet transform" << endl;
+        cout << "filter '" << filterName << "' / " << level << " level(s)" << endl;
+        FOR_STAT_ARRAY(in, s)
+        {
+            dataDWT = dwt.forward(in[s], level);
+            for (unsigned int l = 0; l < level; ++l)
+            {
+                out[l][s]     = dataDWT[l].first;
+                outh[l][s]    = dataDWT[l].second;
+                concath[l][s] = DWT::concat(dataDWT, l, true);
+            }
+            concat[s] = DWT::concat(dataDWT);
+        }
+        cout << "Data CDR " << Math::cdr(in.correlationMatrix()) << " dB" << endl;
+        cout << "DWT CDR " << Math::cdr(concat.correlationMatrix()) << " dB" << endl;
+        for (unsigned int l = 0; l < level; ++l)
+        {
+            cout << "DWT level " << l << " CDR: L= ";
+            cout << Math::cdr(out[l].correlationMatrix()) << " dB / H= ";
+            cout << Math::cdr(outh[l].correlationMatrix()) << " dB" << endl;
+        }
+        for (unsigned int l = 0; l < level; ++l)
+        {
+            cout << "DWT detail level " << l << " CDR: ";
+            cout << Math::cdr(concath[l].correlationMatrix()) << " dB" << endl;
+        }
+    }
+    else
+    {
+        Index ssn = n;
+
+        cout << "-- compute inverse discrete wavelet transform" << endl;
+        cout << "filter '" << filterName << "' / " << level << " level(s)" << endl;
+        for (int l = level - 1; l >= 0; --l)
+        {
+            dataDWT[l].first.resize(ssn);
+            dataDWT[l].second.resize(ssn);
+            dataDWT[l].first.fill(0.);
+            dataDWT[l].second.fill(0.);
+            ssn *= 2;
+        }
+        FOR_STAT_ARRAY(in, s)
+        {
+            dataDWT.back().first = in[s];
+            out[0][s] = dwt.backward(dataDWT);
+        }
+    }
+    if (!outFilename.empty())
+    {
+        if (!ss and saveAll)
+        {
+            Latan::mkdir(outFilename);
+            for (unsigned int l = 0; l < level; ++l)
+            {
+                Io::save<DMatSample>(out[l], outFilename + "/L" + strFrom(l) + ".h5");
+                Io::save<DMatSample>(outh[l], outFilename + "/H" + strFrom(l) + ".h5");
+                Io::save<DMatSample>(concath[l], outFilename + "/concatH" + strFrom(l) + ".h5");
+            }
+            Io::save<DMatSample>(concat, outFilename + "/concat.h5");
+        }
+        else
+        {
+            Io::save<DMatSample>(out.back(), outFilename);
+        }
+    }
+
+    // plot ////////////////////////////////////////////////////////////////////
+    if (doPlot)
+    {
+        Plot p;
+        DVec x;
+
+        x.setLinSpaced(n, 0., n - 1.);
+        p << PlotRange(Axis::x, 0., n);
+        p << Title("original") << PlotData(x, in);
+        if (!ss)
+        {
+            Index  ln = n, step = 1;
+
+            for (unsigned int l = 0; l < level; ++l)
+            {
+                ln   /= 2;
+                step *= 2;
+                x.setLinSpaced(ln, 0., n - step);
+                p << Title("level " + strFrom(l + 1) + " L") << PlotData(x, out[l]);
+                p << Title("level " + strFrom(l + 1) + " H") << PlotData(x, outh[l]);
+            }
+            p.display();
+        }
+        else
+        {
+            double step = 1.;
+            DVec   err;
+
+            for (unsigned int l = 0; l < level; ++l)
+            {
+                step /= 2.;
+            }
+            x.setLinSpaced(out[0][central].size(), 0., n - step);
+            err = out[0].variance().cwiseSqrt();
+            p << Title("supersampled") << Color("3") << PlotPredBand(x, out[0][central], err);
+            p << Color("3") << PlotLine(x, out[0][central]);
+            p.display();
+        }
+    }
+
+    return EXIT_SUCCESS;
+}

utils/sample-merge.cpp

@@ -0,0 +1,94 @@
+/*
+ * sample-merge.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/Core/OptParser.hpp>
+#include <LatAnalyze/Io/Io.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char *argv[])
+{
+    // argument parsing ////////////////////////////////////////////////////////
+    OptParser      opt;
+    bool           parsed;
+    string         outFileName = "";
+    vector<string> fileName;
+    unsigned int   n = 0;
+
+    opt.addOption("o", "output", OptParser::OptType::value  , true,
+                  "output file name (default: result not saved)");
+    opt.addOption("" , "help"  , OptParser::OptType::trigger, true,
+                  "show this help message and exit");
+    parsed = opt.parse(argc, argv);
+    if (!parsed or (opt.getArgs().size() < 1) or opt.gotOption("help"))
+    {
+        cerr << "usage: " << argv[0];
+        cerr << " <sample 1> ... <sample n>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
+        
+        return EXIT_FAILURE;
+    }
+    n           = opt.getArgs().size();
+    outFileName = opt.optionValue("o");
+    for (unsigned int i = 0; i < n; ++i)
+    {
+        fileName.push_back(opt.getArgs()[i]);
+    }
+
+    // figure out dimensions ///////////////////////////////////////////////////
+    Index      nCol, nSample, totRows;
+    DMatSample buf;
+
+    buf     = Io::load<DMatSample>(fileName[0]);
+    nSample = buf.size();
+    totRows = buf[central].rows();
+    nCol    = buf[central].cols();
+    for (unsigned int i = 1; i < n; ++i)
+    {
+        buf = Io::load<DMatSample>(fileName[i]);
+        if (buf.size() != nSample)
+        {
+            LATAN_ERROR(Size, "sample size mismatch");
+        }
+        if (buf[central].cols() != nCol)
+        {
+            LATAN_ERROR(Size, "column size mismatch");
+        }
+        totRows += buf[central].rows();
+    }
+
+    // merge along rows ////////////////////////////////////////////////////////
+    DMatSample out(nSample, totRows, nCol);
+    Index      rowo = 0, r;
+
+    for (unsigned int i = 0; i < n; ++i)
+    {
+        buf = Io::load<DMatSample>(fileName[i]);
+        r   = buf[central].rows();
+        out.block(rowo, 0, r, nCol) = buf;
+        rowo += r;
+    }
+    if (!outFileName.empty())
+    {
+        Io::save<DMatSample>(out, outFileName);
+    }
+
+    return EXIT_SUCCESS;
+}

utils/sample-plot-corr.cpp

@@ -67,7 +67,9 @@ int main(int argc, char *argv[])
     sample = sample.block(0, 0, sample[central].rows(), 1);
     var    = sample.varianceMatrix();
     corr   = sample.correlationMatrix();
-    p << PlotMatrix(corr);
+
+    cout << "dynamic range " << Math::cdr(corr) << " dB" << endl;
+    p << PlotCorrMatrix(corr);
     p.display();
     if (!outVarName.empty())
     {