Merge 269d0c338e into 05138baa08

.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
 

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.**

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

@@ -29,7 +29,8 @@ using namespace Latan;
  ******************************************************************************/
 DMat MATH_NAMESPACE::varToCorr(const DMat &var)
 {
-    DMat res = var, invDiag = res.diagonal();
+    DMat res = var;
+    DVec invDiag = res.diagonal();
     
     invDiag = invDiag.cwiseInverse().cwiseSqrt();
     res     = (invDiag*invDiag.transpose()).cwiseProduct(res);
@@ -37,6 +38,28 @@ DMat MATH_NAMESPACE::varToCorr(const DMat &var)
     return res;
 }
 
+DMat MATH_NAMESPACE::corrToVar(const DMat &corr, const DVec &varDiag)
+{
+    DMat res = corr;
+    DVec varSqrtDiag = varDiag.cwiseSqrt();
+    
+    res = (varSqrtDiag*varSqrtDiag.transpose()).cwiseProduct(res);
+
+    return res;
+}
+
+double MATH_NAMESPACE::svdDynamicRange(const DMat &mat)
+{
+    DVec s = mat.singularValues();
+
+    return s.maxCoeff()/s.minCoeff();
+}
+
+double MATH_NAMESPACE::svdDynamicRangeDb(const DMat &mat)
+{
+    return 10.*log10(svdDynamicRange(mat));
+}
+
 /******************************************************************************
  *                             Standard C functions                           *
  ******************************************************************************/

lib/Core/Math.hpp

@@ -70,6 +70,11 @@ namespace MATH_NAMESPACE
     
     // convert variance matrix to correlation matrix
     DMat varToCorr(const DMat &var);
+    DMat corrToVar(const DMat &corr, const DVec &varDiag);
+
+    // matrix SVD dynamic range
+    double svdDynamicRange(const DMat &mat);
+    double svdDynamicRangeDb(const DMat &mat);
     
     // Constants
     constexpr double pi  = 3.1415926535897932384626433832795028841970;

lib/Core/MathInterpreter.cpp

@@ -217,146 +217,25 @@ void RunContext::reset(void)
 #define CODE_WIDTH 6
 #define CODE_MOD   setw(CODE_WIDTH) << left
 
-// Instruction operator ////////////////////////////////////////////////////////
-ostream &Latan::operator<<(ostream& out, const Instruction& ins)
+auto readConstant(Program::const_iterator ip)
+    -> std::tuple<double, Program::const_iterator>
 {
-    ins.print(out);
-    
-    return out;
+    double value = 0.0;
+    std::copy(ip, ip + sizeof(double), reinterpret_cast<std::uint8_t*>(&value));
+
+    return std::make_tuple(value, ip + sizeof(double));
 }
 
-// Push constructors ///////////////////////////////////////////////////////////
-Push::Push(const double val)
-: type_(ArgType::Constant)
-, val_(val)
-, address_(0)
-, name_("")
-{}
-
-Push::Push(const unsigned int address, const string &name)
-: type_(ArgType::Variable)
-, val_(0.0)
-, address_(address)
-, name_(name)
-{}
-
-// Push execution //////////////////////////////////////////////////////////////
-void Push::operator()(RunContext &context) const
+auto readAddress(Program::const_iterator ip)
+    -> std::tuple<unsigned int, Program::const_iterator>
 {
-    if (type_ == ArgType::Constant)
-    {
-        context.stack().push(val_);
-    }
-    else
-    {
-        context.stack().push(context.getVariable(address_));
-    }
-    context.incrementInsIndex();
+    unsigned int address = 0.0;
+    const auto end = ip + sizeof(unsigned int);
+    std::copy(ip, end, reinterpret_cast<std::uint8_t*>(&address));
+
+    return std::make_tuple(address, end);
 }
 
-// Push print //////////////////////////////////////////////////////////////////
-void Push::print(ostream &out) const
-{
-    out << CODE_MOD << "push";
-    if (type_ == ArgType::Constant)
-    {
-        out << CODE_MOD << val_;
-    }
-    else
-    {
-        out << CODE_MOD << name_ << " @v" << address_;
-    }
-}
-
-// Pop constructor /////////////////////////////////////////////////////////////
-Pop::Pop(const unsigned int address, const string &name)
-: address_(address)
-, name_(name)
-{}
-
-// Pop execution ///////////////////////////////////////////////////////////////
-void Pop::operator()(RunContext &context) const
-{
-    if (!name_.empty())
-    {
-        context.setVariable(address_, context.stack().top());
-    }
-    context.stack().pop();
-    context.incrementInsIndex();
-}
-
-// Pop print ///////////////////////////////////////////////////////////////////
-void Pop::print(ostream &out) const
-{
-    out << CODE_MOD << "pop" << CODE_MOD << name_ << " @v" << address_;
-}
-
-// Store constructor ///////////////////////////////////////////////////////////
-Store::Store(const unsigned int address, const string &name)
-: address_(address)
-, name_(name)
-{}
-
-// Store execution /////////////////////////////////////////////////////////////
-void Store::operator()(RunContext &context) const
-{
-    if (!name_.empty())
-    {
-        context.setVariable(address_, context.stack().top());
-    }
-    context.incrementInsIndex();
-}
-
-// Store print /////////////////////////////////////////////////////////////////
-void Store::print(ostream &out) const
-{
-    out << CODE_MOD << "store" << CODE_MOD << name_ << " @v" << address_;
-}
-
-// Call constructor ////////////////////////////////////////////////////////////
-Call::Call(const unsigned int address, const string &name)
-: address_(address)
-, name_(name)
-{}
-
-// Call execution //////////////////////////////////////////////////////////////
-void Call::operator()(RunContext &context) const
-{
-    context.stack().push((*context.getFunction(address_))(context.stack()));
-    context.incrementInsIndex();
-}
-
-// Call print //////////////////////////////////////////////////////////////////
-void Call::print(ostream &out) const
-{
-    out << CODE_MOD << "call" << CODE_MOD << name_ << " @f" << address_;
-}
-
-// Math operations /////////////////////////////////////////////////////////////
-#define DEF_OP(name, nArg, exp, insName)\
-void name::operator()(RunContext &context) const\
-{\
-    double x[nArg];\
-    for (int i = 0; i < nArg; ++i)\
-    {\
-        x[nArg-1-i] = context.stack().top();\
-        context.stack().pop();\
-    }\
-    context.stack().push(exp);\
-    context.incrementInsIndex();\
-}\
-void name::print(ostream &out) const\
-{\
-    out << CODE_MOD << insName;\
-}
-
-DEF_OP(Neg, 1, -x[0],          "neg")
-DEF_OP(Add, 2, x[0] + x[1],    "add")
-DEF_OP(Sub, 2, x[0] - x[1],    "sub")
-DEF_OP(Mul, 2, x[0]*x[1],      "mul")
-DEF_OP(Div, 2, x[0]/x[1],      "div")
-DEF_OP(Pow, 2, pow(x[0],x[1]), "pow")
-
 /******************************************************************************
  *                        ExprNode implementation                             *
  ******************************************************************************/
@@ -442,29 +321,35 @@ ostream &Latan::operator<<(ostream &out, const ExprNode &n)
     return out;
 }
 
-#define PUSH_INS(program, type, ...)\
-program.push_back(unique_ptr<type>(new type(__VA_ARGS__)))
-#define GET_ADDRESS(address, table, name)\
-try\
-{\
-    address = (table).at(name);\
-}\
-catch (out_of_range)\
-{\
-    address         = (table).size();\
-    (table)[(name)] = address;\
-}\
+// Bytecode helper functions ///////////////////////////////////////////////////
+void pushInstruction(Program &program, const Instruction instruction) {
+    program.push_back(static_cast<std::uint8_t>(instruction));
+}
+
+void pushAddress(Program &program, const unsigned int address) {
+    const auto address_ptr = reinterpret_cast<const std::uint8_t*>(&address);
+    const auto size = sizeof(unsigned int);
+    program.insert(program.end(), address_ptr, address_ptr + size);
+}
+
+void pushConstant(Program &program, const double value) {
+    const auto value_ptr = reinterpret_cast<const std::uint8_t*>(&value);
+    const auto size = sizeof(double);
+    program.insert(program.end(), value_ptr, value_ptr + size);
+}
 
 // VarNode compile /////////////////////////////////////////////////////////////
 void VarNode::compile(Program &program, RunContext &context) const
 {
-    PUSH_INS(program, Push, context.getVariableAddress(getName()), getName());
+    pushInstruction(program, Instruction::LOAD);
+    pushAddress(program, context.getVariableAddress(getName()));
 }
 
 // CstNode compile /////////////////////////////////////////////////////////////
 void CstNode::compile(Program &program, RunContext &context __dumb) const
 {
-    PUSH_INS(program, Push, strTo<double>(getName()));
+    pushInstruction(program, Instruction::CONST);
+    pushConstant(program, strTo<double>(getName()));
 }
 
 // SemicolonNode compile ///////////////////////////////////////////////////////
@@ -482,6 +367,10 @@ void SemicolonNode::compile(Program &program, RunContext &context) const
         {
             n[i].compile(program, context);
         }
+        // Where a variable has just been assigned, pop it off the stack.
+        if (isAssign) {
+            pushInstruction(program, Instruction::POP);
+        }
     }
 }
 
@@ -492,19 +381,10 @@ void AssignNode::compile(Program &program, RunContext &context) const
     
     if (isDerivedFrom<VarNode>(&n[0]))
     {
-        bool hasSemicolonParent = isDerivedFrom<SemicolonNode>(getParent());
-        unsigned int address;
-        
         n[1].compile(program, context);
-        address = context.addVariable(n[0].getName());
-        if (hasSemicolonParent)
-        {
-            PUSH_INS(program, Pop, address, n[0].getName());
-        }
-        else
-        {
-            PUSH_INS(program, Store, address, n[0].getName());
-        }
+        const unsigned int address = context.addVariable(n[0].getName());
+        pushInstruction(program, Instruction::STORE);
+        pushAddress(program, address);
     }
     else
     {
@@ -519,19 +399,37 @@ void AssignNode::compile(Program &program, RunContext &context) const
 
 void MathOpNode::compile(Program &program, RunContext &context) const
 {
+#define PUSH_BINARY_OP(op, instruction) \
+    case op: \
+        pushInstruction(program, Instruction::instruction); \
+        break;
+
     auto &n = *this;
     
     for (Index i = 0; i < n.getNArg(); ++i)
     {
         n[i].compile(program, context);
     }
-    IFNODE("-", 1)   PUSH_INS(program, Neg,);
-    ELIFNODE("+", 2) PUSH_INS(program, Add,);
-    ELIFNODE("-", 2) PUSH_INS(program, Sub,);
-    ELIFNODE("*", 2) PUSH_INS(program, Mul,);
-    ELIFNODE("/", 2) PUSH_INS(program, Div,);
-    ELIFNODE("^", 2) PUSH_INS(program, Pow,);
-    ELSE LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+    if (n.getName() == "-" and n.getNArg() == 1) {
+        pushInstruction(program, Instruction::NEG);
+        return;
+    }
+
+    if (getNArg() != 2) {
+        LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+    }
+
+    switch (getName()[0]) {
+    PUSH_BINARY_OP('+', ADD)
+    PUSH_BINARY_OP('-', SUB)
+    PUSH_BINARY_OP('*', MUL)
+    PUSH_BINARY_OP('/', DIV)
+    PUSH_BINARY_OP('^', POW)
+    default:
+        LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+    }
+
+#undef PUSH_BINARY_OP
 }
 
 // FuncNode compile ////////////////////////////////////////////////////////////
@@ -543,7 +441,8 @@ void FuncNode::compile(Program &program, RunContext &context) const
     {
         n[i].compile(program, context);
     }
-    PUSH_INS(program, Call, context.getFunctionAddress(getName()), getName());
+    pushInstruction(program, Instruction::CALL);
+    pushAddress(program, context.getFunctionAddress(getName()));
 }
 
 // ReturnNode compile ////////////////////////////////////////////////////////////
@@ -552,7 +451,7 @@ void ReturnNode::compile(Program &program, RunContext &context) const
     auto &n = *this;
     
     n[0].compile(program, context);
-    program.push_back(nullptr);
+    pushInstruction(program, Instruction::RET);
 }
 
 /******************************************************************************
@@ -582,7 +481,7 @@ MathInterpreter::MathInterpreter(const std::string &code)
 // access //////////////////////////////////////////////////////////////////////
 const Instruction * MathInterpreter::operator[](const Index i) const
 {
-    return program_[i].get();
+    return reinterpret_cast<const Instruction*>(&program_[i]);
 }
 
 const ExprNode * MathInterpreter::getAST(void) const
@@ -592,7 +491,7 @@ const ExprNode * MathInterpreter::getAST(void) const
 
 void MathInterpreter::push(const Instruction *i)
 {
-    program_.push_back(unique_ptr<const Instruction>(i));
+    pushInstruction(program_, *i);
 }
 
 // initialization //////////////////////////////////////////////////////////////
@@ -695,7 +594,7 @@ void MathInterpreter::compile(RunContext &context)
             root_->compile(program_, context);
             for (unsigned int i = 0; i < program_.size(); ++i)
             {
-                if (!program_[i])
+                if (static_cast<Instruction>(program_[i]) == Instruction::RET)
                 {
                     gotReturn = true;
                     program_.resize(i);
@@ -726,20 +625,145 @@ void MathInterpreter::operator()(RunContext &context)
 
 void MathInterpreter::execute(RunContext &context) const
 {
-    context.setInsIndex(0);
-    while (context.getInsIndex() != program_.size())
-    {
-        (*(program_[context.getInsIndex()]))(context);
+#define BINARY_OP_CASE(instruction, expr) \
+        case Instruction::instruction: { \
+            const auto second = context.stack().top(); \
+            context.stack().pop(); \
+            const auto first = context.stack().top(); \
+            context.stack().pop(); \
+            context.stack().push(expr); \
+            break; \
+        }
+    auto ip = program_.begin();
+
+    while (ip != program_.end()) {
+        const auto instruction = static_cast<Instruction>(*ip);
+        ip++;
+
+        switch (instruction) {
+        BINARY_OP_CASE(ADD, first + second)
+        BINARY_OP_CASE(SUB, first - second)
+        BINARY_OP_CASE(MUL, first * second)
+        BINARY_OP_CASE(DIV, first / second)
+        BINARY_OP_CASE(POW, std::pow(first, second))
+        case Instruction::NEG: {
+            const auto operand = context.stack().top();
+            context.stack().pop();
+            context.stack().push(-operand);
+            break;
+        }
+        case Instruction::CONST: {
+            double value = 0.0;
+            std::tie(value, ip) = readConstant(ip);
+            context.stack().push(value);
+            break;
+        }
+        case Instruction::POP:
+            context.stack().pop();
+            break;
+        case Instruction::LOAD: {
+            unsigned int address = 0;
+            std::tie(address, ip) = readAddress(ip);
+            context.stack().push(context.getVariable(address));
+            break;
+        }
+        case Instruction::STORE: {
+            unsigned int address = 0;
+            std::tie(address, ip) = readAddress(ip);
+            context.setVariable(address, context.stack().top());
+            break;
+        }
+        case Instruction::CALL: {
+            unsigned int address = 0;
+            std::tie(address, ip) = readAddress(ip);
+            auto& stack = context.stack();
+            stack.push((*context.getFunction(address))(stack));
+            break;
+        }
+        case Instruction::RET:
+            break;
+        }
     }
+
+#undef BINARY_OP_CASE
+}
+
+Program::const_iterator instructionToStream(
+    ostream &out, Program::const_iterator ip)
+{
+    const auto instruction = static_cast<Instruction>(*ip);
+    ip++;
+
+    switch (instruction) {
+    case Instruction::ADD:
+        out << "ADD";
+        break;
+    case Instruction::SUB:
+        out << "SUB";
+        break;
+    case Instruction::MUL:
+        out << "MUL";
+        break;
+    case Instruction::DIV:
+        out << "DIV";
+        break;
+    case Instruction::POW:
+        out << "POW";
+        break;
+    case Instruction::NEG:
+        out << "NEG";
+        break;
+    case Instruction::CONST: {
+        double value = 0.0;
+        std::tie(value, ip) = readConstant(ip);
+        out << CODE_MOD << setfill(' ') << "CONST" << value;
+        break;
+    }
+    case Instruction::POP:
+        out << "POP";
+        break;
+    case Instruction::LOAD: {
+        unsigned int address = 0;
+        std::tie(address, ip) = readAddress(ip);
+        out << CODE_MOD << setfill(' ') << "LOAD" << address;
+        break;
+    }
+    case Instruction::STORE: {
+        unsigned int address = 0;
+        std::tie(address, ip) = readAddress(ip);
+        out << CODE_MOD << setfill(' ') << "STORE" << address;
+        break;
+    }
+    case Instruction::CALL: {
+        unsigned int address = 0;
+        std::tie(address, ip) = readAddress(ip);
+        out << CODE_MOD << setfill(' ') << "CALL" << address;
+        break;
+    }
+    case Instruction::RET:
+        out << "RET";
+        break;
+    }
+
+    return ip;
+}
+
+ostream &programToStream(ostream &out, const Program &program)
+{
+    auto ip = program.begin();
+
+    while (ip != program.end()) {
+        const auto i = std::distance(program.begin(), ip);
+        cout << setw(4) << setfill('0') << right << i << "  ";
+        ip = instructionToStream(out, ip);
+        out << '\n';
+    }
+
+    return out;
 }
 
 // IO //////////////////////////////////////////////////////////////////////////
 ostream &Latan::operator<<(ostream &out, const MathInterpreter &program)
 {
-    for (unsigned int i = 0; i < program.program_.size(); ++i)
-    {
-        out << *(program.program_[i]) << endl;
-    }
-    
-    return out;
+    return programToStream(out, program.program_);
 }

lib/Core/MathInterpreter.hpp

@@ -20,6 +20,8 @@
 #ifndef Latan_MathInterpreter_hpp_
 #define	Latan_MathInterpreter_hpp_
 
+#include <cstdint>
+
 #include <LatAnalyze/Functional/Function.hpp>
 #include <LatAnalyze/Global.hpp>
 #include <LatAnalyze/Core/ParserState.hpp>
@@ -79,108 +81,26 @@ private:
  *                         Instruction classes                                *
  ******************************************************************************/
 // Abstract base
-class Instruction
-{
-public:
-    // destructor
-    virtual ~Instruction(void) = default;
-    // instruction execution
-    virtual void operator()(RunContext &context) const = 0;
-    friend std::ostream & operator<<(std::ostream &out, const Instruction &ins);
-private:
-    virtual void print(std::ostream &out) const = 0;
+
+enum class Instruction : std::uint8_t {
+    ADD,
+    SUB,
+    MUL,
+    DIV,
+    POW,
+    NEG,
+    CONST,
+    POP,
+    LOAD,
+    STORE,
+    CALL,
+    RET,
 };
 
 std::ostream & operator<<(std::ostream &out, const Instruction &ins);
 
 // Instruction container
-typedef std::vector<std::unique_ptr<const Instruction>> Program;
-
-// Push
-class Push: public Instruction
-{
-private:
-    enum class ArgType
-    {
-        Constant = 0,
-        Variable = 1
-    };
-public:
-    //constructors
-    explicit Push(const double val);
-    explicit Push(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    ArgType      type_;
-    double       val_;
-    unsigned int address_;
-    std::string  name_;
-};
-
-// Pop
-class Pop: public Instruction
-{
-public:
-    //constructor
-    explicit Pop(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    unsigned int address_;
-    std::string name_;
-};
-
-// Store
-class Store: public Instruction
-{
-public:
-    //constructor
-    explicit Store(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    unsigned int address_;
-    std::string name_;
-};
-
-// Call function
-class Call: public Instruction
-{
-public:
-    //constructor
-    explicit Call(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    unsigned int address_;
-    std::string name_;
-};
-
-// Floating point operations
-#define DECL_OP(name)\
-class name: public Instruction\
-{\
-public:\
-    virtual void operator()(RunContext &context) const;\
-private:\
-    virtual void print(std::ostream &out) const;\
-}
-
-DECL_OP(Neg);
-DECL_OP(Add);
-DECL_OP(Sub);
-DECL_OP(Mul);
-DECL_OP(Div);
-DECL_OP(Pow);
+typedef std::vector<std::uint8_t> Program;
 
 /******************************************************************************
  *                            Expression node classes                         *

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 ////////////////////////////////////////////////////

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                                 *
  ******************************************************************************/

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       \
@@ -75,6 +78,7 @@ HPPFILES =                           \
     Core/OptParser.hpp               \
     Core/ParserState.hpp             \
     Core/Plot.hpp                    \
+    Core/ThreadPool.hpp              \
     Core/stdincludes.hpp             \
     Core/Utilities.hpp               \
     Functional/CompiledFunction.hpp  \
@@ -90,6 +94,8 @@ HPPFILES =                           \
     Io/IoObject.hpp                  \
     Io/XmlReader.hpp                 \
     Numerical/Derivative.hpp         \
+    Numerical/DWT.hpp                \
+    Numerical/DWTFilters.hpp         \
     Numerical/FFT.hpp                \
     Numerical/GslFFT.hpp             \
     Numerical/GslHybridRootFinder.hpp\

lib/Numerical/DWT.cpp

@@ -0,0 +1,137 @@
+/*
+ * DWT.cpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Numerical/DWT.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                            DWT implementation                              *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+DWT::DWT(const DWTFilter &filter)
+: filter_(filter)
+{}
+
+// convolution primitive ///////////////////////////////////////////////////////
+void DWT::filterConvolution(DVec &out, const DVec &data, 
+                            const std::vector<double> &filter, const Index offset)
+{
+    Index n = data.size(), nf = n*filter.size();
+
+    out.resize(n);
+    out.fill(0.);
+    for (unsigned int i = 0; i < filter.size(); ++i)
+    {
+        FOR_VEC(out, j)
+        {
+            out(j) += filter[i]*data((j + i + nf - offset) % n);
+        }
+    }
+}
+
+// downsampling/upsampling primitives //////////////////////////////////////////
+void DWT::downsample(DVec &out, const DVec &in)
+{
+    if (out.size() < in.size()/2)
+    {
+        LATAN_ERROR(Size, "output vector smaller than half the input vector size");
+    }
+    for (Index i = 0; i < in.size(); i += 2)
+    {
+        out(i/2) = in(i);
+    }
+}
+
+void DWT::upsample(DVec &out, const DVec &in)
+{
+    if (out.size() < 2*in.size())
+    {
+        LATAN_ERROR(Size, "output vector smaller than twice the input vector size");
+    }
+    out.segment(0, 2*in.size()).fill(0.);
+    for (Index i = 0; i < in.size(); i ++)
+    {
+        out(2*i) = in(i);
+    }
+}
+
+// DWT /////////////////////////////////////////////////////////////////////////
+std::vector<DWT::DWTLevel> 
+DWT::forward(const DVec &data, const unsigned int level) const
+{
+    std::vector<DWTLevel> dwt(level);
+    DVec                  *finePt = const_cast<DVec *>(&data);
+    DVec                  tmp; 
+    Index                 n = data.size(), o = filter_.fwdL.size()/2, minSize;
+
+    minSize = 1;
+    for (unsigned int l = 0; l < level; ++l) minSize *= 2;
+    if (n < minSize)
+    {
+        LATAN_ERROR(Size, "data vector too small for a " + strFrom(level) 
+                          + "-level DWT (data size is " + strFrom(n) + ")");
+    }
+    for (unsigned int l = 0; l < level; ++l)
+    {
+        n /= 2;
+        dwt[l].first.resize(n);
+        dwt[l].second.resize(n);
+        filterConvolution(tmp, *finePt, filter_.fwdL, o);
+        downsample(dwt[l].first, tmp);
+        filterConvolution(tmp, *finePt, filter_.fwdH, o);
+        downsample(dwt[l].second, tmp);
+        finePt = &dwt[l].first;
+    }
+
+    return dwt;
+}
+
+DVec DWT::backward(const std::vector<DWTLevel>& dwt) const
+{
+    unsigned int level = dwt.size();
+    Index        n = dwt.back().second.size(), o = filter_.bwdL.size()/2 - 1;
+    DVec         res, tmp, conv;
+
+    res = dwt.back().first;
+    for (int l = level - 2; l >= 0; --l)
+    {
+        n *= 2;
+        if (dwt[l].second.size() != n)
+        {
+            LATAN_ERROR(Size, "DWT result size mismatch");
+        }
+    }
+    n = dwt.back().second.size();
+    for (int l = level - 1; l >= 0; --l)
+    {
+        n *= 2;
+        tmp.resize(n);
+        upsample(tmp, res);
+        filterConvolution(conv, tmp, filter_.bwdL, o);
+        res = conv;
+        upsample(tmp, dwt[l].second);
+        filterConvolution(conv, tmp, filter_.bwdH, o);
+        res += conv;
+    }
+
+    return res;
+}

lib/Numerical/DWT.hpp

@@ -0,0 +1,55 @@
+/*
+ * DWT.hpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_DWT_hpp_
+#define Latan_DWT_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Numerical/DWTFilters.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                   Discrete wavelet transform class                         *
+ ******************************************************************************/
+class DWT
+{
+public:
+  typedef std::pair<DVec, DVec> DWTLevel;
+public:
+    // constructor
+    DWT(const DWTFilter &filter);
+    // destructor
+    virtual ~DWT(void) = default;
+    // convolution primitive
+    static void filterConvolution(DVec &out, const DVec &data, 
+                                  const std::vector<double> &filter, const Index offset);
+    // downsampling/upsampling primitives
+    static void downsample(DVec &out, const DVec &in);
+    static void upsample(DVec &out, const DVec &in);
+    // DWT
+    std::vector<DWTLevel> forward(const DVec &data, const unsigned int level) const;
+    DVec                  backward(const std::vector<DWTLevel>& dwt) const;
+private:
+    DWTFilter filter_;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_DWT_hpp_

lib/Numerical/DWTFilters.cpp

@@ -0,0 +1,528 @@
+/*
+ * DWTFilters.cpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Numerical/DWTFilters.hpp>
+#include <LatAnalyze/includes.hpp>
+
+// cf. http://wavelets.pybytes.com
+// *here we implement the reverse filters more convenient for convolutions*
+
+using namespace std;
+using namespace Latan;
+
+#define FILTDICT(x) {#x, &DWTFilters::x}
+
+std::map<std::string, const DWTFilter *> DWTFilters::fromName = {
+    FILTDICT(haar),
+    FILTDICT(db2),
+    FILTDICT(db3),
+    FILTDICT(db4),
+    FILTDICT(db5),
+    FILTDICT(db6),
+    FILTDICT(bior13),
+    FILTDICT(bior15),
+    FILTDICT(bior22),
+    FILTDICT(bior24),
+    FILTDICT(bior31),
+    FILTDICT(bior33),
+    FILTDICT(bior35)
+};
+
+DWTFilter DWTFilters::haar = {
+    // fwdL
+    {0.7071067811865476,
+    0.7071067811865476},
+    // fwdH
+    {0.7071067811865476,
+    -0.7071067811865476},
+    // bwdL
+    {0.7071067811865476,
+    0.7071067811865476},
+    // bwdH
+    {-0.7071067811865476,
+    0.7071067811865476}
+};
+
+DWTFilter DWTFilters::db2 = {
+    // fwdL
+    {0.48296291314469025,
+    0.836516303737469,
+    0.22414386804185735,
+    -0.12940952255092145},
+    // fwdH
+    {-0.12940952255092145,
+    -0.22414386804185735,
+    0.836516303737469,
+    -0.48296291314469025},
+    // bwdL
+    {-0.12940952255092145,
+    0.22414386804185735,
+    0.836516303737469,
+    0.48296291314469025},
+    // bwdH
+    {-0.48296291314469025,
+    0.836516303737469,
+    -0.22414386804185735,
+    -0.12940952255092145}
+};
+
+DWTFilter DWTFilters::db3 = {
+    // fwdL
+    {0.3326705529509569,
+    0.8068915093133388,
+    0.4598775021193313,
+    -0.13501102001039084,
+    -0.08544127388224149,
+    0.035226291882100656},
+    // fwdH
+    {0.035226291882100656,
+    0.08544127388224149,
+    -0.13501102001039084,
+    -0.4598775021193313,
+    0.8068915093133388,
+    -0.3326705529509569},
+    // bwdL
+    {0.035226291882100656,
+    -0.08544127388224149,
+    -0.13501102001039084,
+    0.4598775021193313,
+    0.8068915093133388,
+    0.3326705529509569},
+    // bwdH
+    {-0.3326705529509569,
+    0.8068915093133388,
+    -0.4598775021193313,
+    -0.13501102001039084,
+    0.08544127388224149,
+    0.035226291882100656}
+};
+
+DWTFilter DWTFilters::db4 = {
+    // fwdL
+    {0.23037781330885523,
+    0.7148465705525415,
+    0.6308807679295904,
+    -0.02798376941698385,
+    -0.18703481171888114,
+    0.030841381835986965,
+    0.032883011666982945,
+    -0.010597401784997278},
+    // fwdH
+    {-0.010597401784997278,
+    -0.032883011666982945,
+    0.030841381835986965,
+    0.18703481171888114,
+    -0.02798376941698385,
+    -0.6308807679295904,
+    0.7148465705525415,
+    -0.23037781330885523},
+    // bwdL
+    {-0.010597401784997278,
+    0.032883011666982945,
+    0.030841381835986965,
+    -0.18703481171888114,
+    -0.02798376941698385,
+    0.6308807679295904,
+    0.7148465705525415,
+    0.23037781330885523},
+    // bwdH
+    {-0.23037781330885523,
+    0.7148465705525415,
+    -0.6308807679295904,
+    -0.02798376941698385,
+    0.18703481171888114,
+    0.030841381835986965,
+    -0.032883011666982945,
+    -0.010597401784997278}
+};
+
+DWTFilter DWTFilters::db5 = {
+    // fwdL
+    {0.160102397974125,
+    0.6038292697974729,
+    0.7243085284385744,
+    0.13842814590110342,
+    -0.24229488706619015,
+    -0.03224486958502952,
+    0.07757149384006515,
+    -0.006241490213011705,
+    -0.012580751999015526,
+    0.003335725285001549},
+    // fwdH
+    {0.003335725285001549,
+    0.012580751999015526,
+    -0.006241490213011705,
+    -0.07757149384006515,
+    -0.03224486958502952,
+    0.24229488706619015,
+    0.13842814590110342,
+    -0.7243085284385744,
+    0.6038292697974729,
+    -0.160102397974125},
+    // bwdL
+    {0.003335725285001549,
+    -0.012580751999015526,
+    -0.006241490213011705,
+    0.07757149384006515,
+    -0.03224486958502952,
+    -0.24229488706619015,
+    0.13842814590110342,
+    0.7243085284385744,
+    0.6038292697974729,
+    0.160102397974125},
+    // bwdH
+    {-0.160102397974125,
+    0.6038292697974729,
+    -0.7243085284385744,
+    0.13842814590110342,
+    0.24229488706619015,
+    -0.03224486958502952,
+    -0.07757149384006515,
+    -0.006241490213011705,
+    0.012580751999015526,
+    0.003335725285001549}
+};
+
+DWTFilter DWTFilters::db6 = {
+    // fwdL
+    {0.11154074335008017,
+    0.4946238903983854,
+    0.7511339080215775,
+    0.3152503517092432,
+    -0.22626469396516913,
+    -0.12976686756709563,
+    0.09750160558707936,
+    0.02752286553001629,
+    -0.031582039318031156,
+    0.0005538422009938016,
+    0.004777257511010651,
+    -0.00107730108499558},
+    // fwdH
+    {-0.00107730108499558,
+    -0.004777257511010651,
+    0.0005538422009938016,
+    0.031582039318031156,
+    0.02752286553001629,
+    -0.09750160558707936,
+    -0.12976686756709563,
+    0.22626469396516913,
+    0.3152503517092432,
+    -0.7511339080215775,
+    0.4946238903983854,
+    -0.11154074335008017},
+    // bwdL
+    {-0.00107730108499558,
+    0.004777257511010651,
+    0.0005538422009938016,
+    -0.031582039318031156,
+    0.02752286553001629,
+    0.09750160558707936,
+    -0.12976686756709563,
+    -0.22626469396516913,
+    0.3152503517092432,
+    0.7511339080215775,
+    0.4946238903983854,
+    0.11154074335008017},
+    // bwdH
+    {-0.11154074335008017,
+    0.4946238903983854,
+    -0.7511339080215775,
+    0.3152503517092432,
+    0.22626469396516913,
+    -0.12976686756709563,
+    -0.09750160558707936,
+    0.02752286553001629,
+    0.031582039318031156,
+    0.0005538422009938016,
+    -0.004777257511010651,
+    -0.00107730108499558}
+};
+
+DWTFilter DWTFilters::bior13 = {
+    // fwdL
+    {-0.08838834764831845,
+    0.08838834764831845,
+    0.7071067811865476,
+    0.7071067811865476,
+    0.08838834764831845,
+    -0.08838834764831845},
+    // fwdH
+    {0.0,
+    0.0,
+    0.7071067811865476,
+    -0.7071067811865476,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.7071067811865476,
+    0.7071067811865476,
+    0.0,
+    0.0},
+    // bwdH
+    {0.08838834764831845,
+    0.08838834764831845,
+    -0.7071067811865476,
+    0.7071067811865476,
+    -0.08838834764831845,
+    -0.08838834764831845}
+};
+
+DWTFilter DWTFilters::bior15 = {
+    // fwdL
+    {0.01657281518405971,
+    -0.01657281518405971,
+    -0.12153397801643787,
+    0.12153397801643787,
+    0.7071067811865476,
+    0.7071067811865476,
+    0.12153397801643787,
+    -0.12153397801643787,
+    -0.01657281518405971,
+    0.01657281518405971},
+    // fwdH
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.7071067811865476,
+    -0.7071067811865476,
+    0.0,
+    0.0,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.7071067811865476,
+    0.7071067811865476,
+    0.0,
+    0.0,
+    0.0,
+    0.0},
+    // bwdH
+    {-0.01657281518405971,
+    -0.01657281518405971,
+    0.12153397801643787,
+    0.12153397801643787,
+    -0.7071067811865476,
+    0.7071067811865476,
+    -0.12153397801643787,
+    -0.12153397801643787,
+    0.01657281518405971,
+    0.01657281518405971}
+};
+
+DWTFilter DWTFilters::bior22 = {
+    // fwdL
+    {-0.1767766952966369,
+    0.3535533905932738,
+    1.0606601717798214,
+    0.3535533905932738,
+    -0.1767766952966369,
+    0.0},
+    // fwdH
+    {0.0,
+    0.0,
+    0.3535533905932738,
+    -0.7071067811865476,
+    0.3535533905932738,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.3535533905932738,
+    0.7071067811865476,
+    0.3535533905932738,
+    0.0},
+    // bwdH
+    {0.1767766952966369,
+    0.3535533905932738,
+    -1.0606601717798214,
+    0.3535533905932738,
+    0.1767766952966369,
+    0.0}
+};
+
+DWTFilter DWTFilters::bior24 = {
+    // fwdL
+    {0.03314563036811942,
+    -0.06629126073623884,
+    -0.1767766952966369,
+    0.4198446513295126,
+    0.9943689110435825,
+    0.4198446513295126,
+    -0.1767766952966369,
+    -0.06629126073623884,
+    0.03314563036811942,
+    0.0},
+    // fwdH
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.3535533905932738,
+    -0.7071067811865476,
+    0.3535533905932738,
+    0.0,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.3535533905932738,
+    0.7071067811865476,
+    0.3535533905932738,
+    0.0,
+    0.0,
+    0.0},
+    // bwdH
+    {-0.03314563036811942,
+    -0.06629126073623884,
+    0.1767766952966369,
+    0.4198446513295126,
+    -0.9943689110435825,
+    0.4198446513295126,
+    0.1767766952966369,
+    -0.06629126073623884,
+    -0.03314563036811942,
+    0.0}
+};
+
+DWTFilter DWTFilters::bior31 = {
+    // fwdL
+    {-0.3535533905932738,
+    1.0606601717798214,
+    1.0606601717798214,
+    -0.3535533905932738},
+    // fwdH
+    {0.1767766952966369,
+    -0.5303300858899107,
+    0.5303300858899107,
+    -0.1767766952966369},
+    // bwdL
+    {0.1767766952966369,
+    0.5303300858899107,
+    0.5303300858899107,
+    0.1767766952966369},
+    // bwdH
+    {0.3535533905932738,
+    1.0606601717798214,
+    -1.0606601717798214,
+    -0.3535533905932738}
+};
+
+DWTFilter DWTFilters::bior33 = {
+    // fwdL
+    {0.06629126073623884,
+    -0.19887378220871652,
+    -0.15467960838455727,
+    0.9943689110435825,
+    0.9943689110435825,
+    -0.15467960838455727,
+    -0.19887378220871652,
+    0.06629126073623884},
+    // fwdH
+    {0.0,
+    0.0,
+    0.1767766952966369,
+    -0.5303300858899107,
+    0.5303300858899107,
+    -0.1767766952966369,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.1767766952966369,
+    0.5303300858899107,
+    0.5303300858899107,
+    0.1767766952966369,
+    0.0,
+    0.0},
+    // bwdH
+    {-0.06629126073623884,
+    -0.19887378220871652,
+    0.15467960838455727,
+    0.9943689110435825,
+    -0.9943689110435825,
+    -0.15467960838455727,
+    0.19887378220871652,
+    0.06629126073623884}
+};
+
+DWTFilter DWTFilters::bior35 = {
+    // fwdL
+    {-0.013810679320049757,
+    0.04143203796014927,
+    0.052480581416189075,
+    -0.26792717880896527,
+    -0.07181553246425874,
+    0.966747552403483,
+    0.966747552403483,
+    -0.07181553246425874,
+    -0.26792717880896527,
+    0.052480581416189075,
+    0.04143203796014927,
+    -0.013810679320049757},
+    // fwdH
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.1767766952966369,
+    -0.5303300858899107,
+    0.5303300858899107,
+    -0.1767766952966369,
+    0.0,
+    0.0,
+    0.0,
+    0.0},
+    // bwdL
+    {0.0,
+    0.0,
+    0.0,
+    0.0,
+    0.1767766952966369,
+    0.5303300858899107,
+    0.5303300858899107,
+    0.1767766952966369,
+    0.0,
+    0.0,
+    0.0,
+    0.0},
+    // bwdH
+    {0.013810679320049757,
+    0.04143203796014927,
+    -0.052480581416189075,
+    -0.26792717880896527,
+    0.07181553246425874,
+    0.966747552403483,
+    -0.966747552403483,
+    -0.07181553246425874,
+    0.26792717880896527,
+    0.052480581416189075,
+    -0.04143203796014927,
+    -0.013810679320049757}
+};

lib/Numerical/DWTFilters.hpp

@@ -0,0 +1,53 @@
+/*
+ * DWTFilters.hpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_DWTFilters_hpp_
+#define Latan_DWTFilters_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+struct DWTFilter
+{
+    const std::vector<double> fwdL, fwdH, bwdL, bwdH;
+};
+
+namespace DWTFilters
+{
+    extern DWTFilter haar;
+    extern DWTFilter db2;
+    extern DWTFilter db3;
+    extern DWTFilter db4;
+    extern DWTFilter db5;
+    extern DWTFilter db6;
+    extern DWTFilter bior13;
+    extern DWTFilter bior15;
+    extern DWTFilter bior22;
+    extern DWTFilter bior24;
+    extern DWTFilter bior31;
+    extern DWTFilter bior33;
+    extern DWTFilter bior35;
+
+    extern std::map<std::string, const DWTFilter *> fromName;
+}
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_DWTFilters_hpp_

lib/Numerical/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/Statistics/MatSample.hpp

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

lib/Statistics/StatArray.hpp

@@ -51,14 +51,11 @@ 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    meanOld(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;
     // IO type
     virtual IoType getType(void) const;
 public:
@@ -66,7 +63,7 @@ public:
 };
 
 // reduction operations
-namespace ReducOp
+namespace StatOp
 {
     // general templates
     template <typename T>
@@ -148,128 +145,67 @@ 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);
-}
-
-template <typename T, Index os>
-T StatArray<T, os>::covariance(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)
-    {
-        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);
-    }
-    
-    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)
-    {
-        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);
-    }
-    
-    return res/static_cast<double>(m - 1);
-}
-
-template <typename T, Index os>
-T StatArray<T, os>::variance(const Index pos, const Index n) const
-{
-    return covariance(*this, pos, n);
-}
-
-template <typename T, Index os>
-T StatArray<T, os>::varianceMatrix(const Index pos, const Index n) 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);
-
-    invDiag = res.diagonal();
-    invDiag = invDiag.cwiseInverse().cwiseSqrt();
-    res     = (invDiag*invDiag.transpose()).cwiseProduct(res);
+    res -= StatOp::prod<T>(s1, s2)/static_cast<double>(size());
+    res /= static_cast<double>(size() - 1);
 
     return res;
 }
 
-// reduction operations ////////////////////////////////////////////////////////
-namespace ReducOp
+template <typename T, Index os>
+T StatArray<T, os>::variance(void) const
 {
-    template <typename T>
-    inline T sum(const T &a, const T &b)
-    {
-        return a + b;
-    }
+    return covariance(*this);
+}
 
+// reduction operations ////////////////////////////////////////////////////////
+namespace StatOp
+{
     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::svdDynamicRangeDb(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::svdDynamicRangeDb(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;
     string               corrFileName, model, outFileName, outFmt, savePlot;
     Index                ti, tf, shift, nPar, thinning;
     double               svdTol;
@@ -52,6 +53,8 @@ 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,
@@ -74,6 +77,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");
@@ -110,6 +114,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);
@@ -282,7 +297,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,167 @@
+/*
+ * 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/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;
+        
+        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));
+    DWT                   dwt(*DWTFilters::fromName.at(filterName));
+    vector<DWT::DWTLevel> dataDWT(level);
+
+    if (!ss)
+    {
+        cout << "-- compute discrete wavelet transform" << endl;
+        cout << "filter '" << filterName << "' / " << level << " level(s)" << endl;
+        FOR_STAT_ARRAY(in, s)
+        {
+            dataDWT = dwt.forward(in[s].col(0), level);
+            for (unsigned int l = 0; l < level; ++l)
+            {
+                out[l][s]  = dataDWT[l].first;
+                outh[l][s] = dataDWT[l].second;
+            }
+        }
+    }
+    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].col(0);
+            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");
+            }
+        }
+        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::svdDynamicRangeDb(corr) << " dB" << endl;
+    p << PlotCorrMatrix(corr);
     p.display();
     if (!outVarName.empty())
     {