Merge 269d0c338e into c73b609ac5

Changed linear model's parameters

.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
 
@@ -31,7 +32,6 @@ lib/eigen_files.mk
 # CI builds
 ci-scripts/local/*
 
-# CLion stuff
+# VS Code Studio stuff
-.idea/*
+.vscode
-CMakeLists.txt
+*.code-workspace
-cmake-build-debug/*

.travis.yml

@@ -1,9 +1,5 @@
 language: cpp
   
-notifications:
-  email: false
-  slack: ukqcd:mQLXCtz8D2cg89xT8j1a4wku
-  
 cache:
   directories:
     - ci-scripts/local
@@ -12,56 +8,32 @@ cache:
 matrix:
   include:
     - os:        osx
-      osx_image: xcode7.2
+      osx_image: xcode10.1
       compiler: clang
     - os:        osx
-      osx_image: xcode7.2
+      osx_image: xcode10.1
       compiler: gcc
-      env: VERSION=-5
+      env: VERSION=-7
     - compiler: gcc
       addons:
         apt:
           sources:
             - ubuntu-toolchain-r-test
           packages:
-            - g++-4.9
+            - g++-7
-            - libgsl0-dev
             - flex
             - bison
-      env: VERSION=-4.9
+      env: VERSION=-7
-    - compiler: gcc
-      addons:
-        apt:
-          sources:
-            - ubuntu-toolchain-r-test
-          packages:
-            - g++-5
-            - libgsl0-dev
-            - flex
-            - bison
-      env: VERSION=-5
     - compiler: clang
       addons:
         apt:
           sources:
             - ubuntu-toolchain-r-test
           packages:
-            - g++-4.8
+            - g++-7
-            - libgsl0-dev
             - flex
             - bison
-      env: CLANG_LINK=http://llvm.org/releases/3.6.0/clang+llvm-3.6.0-x86_64-linux-gnu-ubuntu-14.04.tar.xz
+      env: CLANG_LINK=http://releases.llvm.org/7.0.1/clang+llvm-7.0.1-x86_64-linux-gnu-ubuntu-14.04.tar.xz
-    - compiler: clang
-      addons:
-        apt:
-          sources:
-            - ubuntu-toolchain-r-test
-          packages:
-            - g++-4.8
-            - libgsl0-dev
-            - flex
-            - bison
-      env: CLANG_LINK=http://llvm.org/releases/3.7.0/clang+llvm-3.7.0-x86_64-linux-gnu-ubuntu-14.04.tar.xz
 
 before_install:
   - export LATDIR=`pwd`
@@ -70,9 +42,8 @@ before_install:
   - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export PATH="${LATDIR}/ci-scripts/clang/bin:${PATH}"; fi
   - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export LD_LIBRARY_PATH="${LATDIR}/ci-scripts/clang/lib:${LD_LIBRARY_PATH}"; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew update; fi
-  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install gsl; fi
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install bison flex cmake; export PATH="/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:$PATH"; fi
-  - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install bison flex; export PATH="/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:$PATH"; fi
+  - if [[ "$TRAVIS_OS_NAME" == "osx" ]] && [[ "$CC" == "gcc" ]]; then brew install gcc@${VERSION#-}; fi
-  - if [[ "$TRAVIS_OS_NAME" == "osx" ]] && [[ "$CC" == "gcc" ]]; then brew install gcc5; fi
   
 install:
   - export CC=$CC$VERSION
@@ -85,7 +56,8 @@ install:
   - ./install-deps.sh `pwd`/local
   - cd ..
   - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then export LD_LIBRARY_PATH=${LATDIR}/ci-scripts/local/lib:$LD_LIBRARY_PATH; fi
+  - if [[ "$CC" == "gcc-7" ]]; then export CXXFLAGS='-Wno-int-in-bool-context'; fi
 
 script:
   - cd ci-scripts
-  - ./install-latan.sh `pwd`/local $TRAVIS_OS_NAME
+  - ./install-latan.sh `pwd`/local

Makefile.am

@@ -1,3 +1,5 @@
 SUBDIRS = lib utils physics examples
 
+bin_SCRIPTS=latan-config
+
 ACLOCAL_AMFLAGS = -I .buildutils/m4

Readme.md

@@ -1,22 +1,6 @@
 # LatAnalyze 
-Contributors: Antonin Portelli, Matt Spraggs
 
-License: GNU General Public License v3
+[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![DOI](https://zenodo.org/badge/10201777.svg)](https://zenodo.org/badge/latestdoi/10201777) [![Build Ubuntu](https://github.com/aportelli/LatAnalyze/actions/workflows/build-ubuntu.yml/badge.svg)](https://github.com/aportelli/LatAnalyze/actions/workflows/build-ubuntu.yml) [![Build macOS](https://github.com/aportelli/LatAnalyze/actions/workflows/build-macos.yml/badge.svg)](https://github.com/aportelli/LatAnalyze/actions/workflows/build-macos.yml)
-
-<table>
-<tr>
-    <td>Last stable release</td>
-    <td><a href="https://travis-ci.org/aportelli/LatAnalyze">
-    <img src="https://travis-ci.org/aportelli/LatAnalyze.svg?branch=master"></a>
-    </td>
-</tr>
-<tr>
-    <td>Development branch</td>
-    <td><a href="https://travis-ci.org/aportelli/LatAnalyze">
-    <img src="https://travis-ci.org/aportelli/LatAnalyze.svg?branch=develop"></a>
-    </td>
-</tr>
-</table>
 
 ## Description
 LatAnalyze is a C++11 library for statistical data analysis based on bootstrap
@@ -41,22 +25,93 @@ The main features are the following:
 The head of the `master` branch always points to the latest stable release. The `develop` branch is the main unstable branch of LatAnalyze.
 
 LatAnalyze is written in C++11 and requires a rather recent C++ compiler to be built. It has been successfully built on various Linux and OS X platforms using clang (from 3.7), GCC (from 4.9) and the Intel C++ compiler (2016).
-The strict dependencies are the [GSL](http://www.gnu.org/software/gsl/) and [LatCore](https://github.com/aportelli/LatCore).
+The only strict dependencies is the [GSL](http://www.gnu.org/software/gsl/).
 Additionally, autoconf, automake (from 1.11), libtool, bison (from 3.0) and flex are necessary to build the library. Unless you use a very exotic system, these tools are standard on any Unix platform and should be already present or easy to install through a package manager.
 Optional dependencies are [HDF5](https://www.hdfgroup.org/HDF5/) (built with C++ support), [Minuit](http://seal.web.cern.ch/seal/snapshot/work-packages/mathlibs/minuit/) and [NLopt](http://ab-initio.mit.edu/wiki/index.php/NLopt).
 
+Below are instructions for a quick installation. For a more customised installation, one first needs to generate the build system by running `./bootstrap.sh` in the root directory. Then the library can be built and installed through the usual GNU mantra `./configure <options> && make && make install`. Use `./configure --help` to obtain a list of possible options for `./configure`. Because Eigen expressions rely a lot on inlining and compiler optimisations it is strongly recommended to set the `CXXFLAGS` variable to `-O3 -march=native -mtune=native`.
 
+### General quick installation
 For a quick installation with all possible extensions execute: 
 ```
-./install-latan.sh <prefix> {osx|linux}
+./install-latan.sh <prefix>
 ``` 
-in the `ci-scripts` directory where `<prefix>` is where you want LatAnalyze (and its dependencies) to be installed and `{osx|linux}` should be the name of your OS. This script will automatically download, build and install LatCore, HDF5, Minuit and NLopt. It is assumed that the GSL can be found in a standard location (_e.g._ it has been installed through a package manager for Linux and is present in `/usr/local` for OS X).
+in the `ci-scripts` directory where `<prefix>` is where you want LatAnalyze (and its dependencies) to be installed. This script will automatically download, build and install GSL, HDF5, Minuit, and NLopt.
 
-For a more customised installation, one first needs to generate the build system by running `./bootstrap.sh` in the root directory. Then the library can be built and installed through the usual GNU mantra `./configure <options> && make && make install`. Use `./configure --help` to obtain a list of possible options for `./configure`. Because Eigen expressions rely a lot on inlining and compiler optimisations it is strongly recommended to set the `CXXFLAGS` variable to `-O3 -march=native -mtune=native`.
+### Quick installation on macOS with dependencies through Homebrew
+All the dependencies of LatAnalyze can be installed through the [Homebrew](https://brew.sh) package manager.
+```
+brew install automake autoconf libtool bison flex gsl minuit2 nlopt hdf5
+```
+Then generate the build system in LatAnalyze main directory by running the `./bootstrap.sh` script. Finally, build the library
+```
+mkdir build
+cd build
+../configure --prefix=<prefix> --with-minuit=/usr/local --with-nlopt=/usr/local \
+             --with-hdf5=/usr/local --with-gsl=/usr/local                       \
+             CXXFLAGS='-g -O3 -march=native -mtune=native'
+make -j <n>
+make install
+```
+where `<prefix>` should be replaced by the desired prefix for LatAnalyze installation, and `<n>` is the number of parallel build processes (typically twice your number of cores).
 
 ## History
+#### v3.5.1
+Various fixes and cleaning of outdated code.
+
+#### v3.5
+Additions:
+* 'Impulse' & line type plots
+* Plot line width & dash modifiers
+* Plot palettes (`category10` by default)
+* Multivariate Gaussian RNG
+* 2-pt fitter 'scan' mode over all possible fit ranges
+* Command line utility for plotting data
+
+Changes:
+* Complete overhaul of the header structure
+* Integration of LatCore in LatAnalyze
+* p-value is now a 2-sided chi^2 test, 1-sided value kept as 'chi^2 CCDF'
+
+Fixes:
+* Matrix plot data now saving correctly
+* Many compatibility fixes
+
+#### v3.4
+Additions:
+* `latan-config` utility to easily compile LatAnalyze-based programs
+* Linear and constant models to the 2-point fitter
+
+Changes:
+* HDF5 is now a compulsory dependency
+
+Fixes:
+* Variance matrix computation fix.
+
+#### v3.3
+Additions:
+* Sample plot CL utility.
+* Infinity as a math constant.
+* Option to dump bootstrap sequence while resampling.
+* FFT through the GSL.
+
+Changes:
+* GSL integrator accepts infinite bounds.
+* `latan-sample-combine` accepts mixes of `DSample` and `DMatSample`.
+* More general `latan-sample-element` command.
+
+#### v3.2.2
+Additions:
+* The math interpreter supports `inf` for infinity.
+
+Changes:
+* Vector version of `setUnidimData`.
+
+Fixes:
+* Variance matrix computation fix.
+
 #### v3.2.1
-Fix:
+Fixes:
 * Wrong argument number check in `latan-resample`
 
 #### v3.2 (needs LatCore 1.1)

bootstrap.sh

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

ci-scripts/install-deps.sh

@@ -1,16 +1,16 @@
 #!/usr/bin/env bash
 
-if (( $# != 1 )); then
+if (( $# != 2 )); then
-  echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
   exit 1
 fi
 PREFIX=$1
+NTASKS=$2
 
 set -ex
 mkdir -p local/build
-for d in nlopt minuit hdf5; do
+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
-./install-latcore.sh ${PREFIX}

ci-scripts/install-gsl.sh

@@ -0,0 +1,27 @@
+#!/usr/bin/env bash
+
+NAME='gsl-2.6'
+
+if (( $# != 2 )); then
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
+  exit 1
+fi
+PREFIX=$1
+NTASKS=$2
+
+set -ex
+INITDIR=$(pwd -P)
+mkdir -p ${PREFIX}
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}/local/build
+wget http://ftpmirror.gnu.org/gsl/${NAME}.tar.gz
+tar -xzvf ${NAME}.tar.gz
+mkdir -p ${NAME}/build
+cd ${NAME}/build
+../configure --prefix=${PREFIX}
+make -j${NTASKS}
+make install
+cd ${INITDIR}/local
+touch .built.gsl
+cd ${INITDIR}

ci-scripts/install-hdf5.sh

@@ -1,22 +1,26 @@
 #!/usr/bin/env bash
 
-NAME='hdf5-1.8.16'
+NAME='hdf5-1.10.8'
 
-if (( $# != 1 )); then
+if (( $# != 2 )); then
-  echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
   exit 1
 fi
 PREFIX=$1
+NTASKS=$2
 
 set -ex
-INITDIR=`pwd`
+INITDIR=$(pwd -P)
-cd local/build
+mkdir -p ${PREFIX}
-wget http://www.hdfgroup.org/ftp/HDF5/releases/${NAME}/src/${NAME}.tar.gz
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}/local/build
+wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/${NAME}/src/${NAME}.tar.gz
 tar -xzvf ${NAME}.tar.gz
 mkdir ${NAME}/build
 cd ${NAME}/build
 ../configure --prefix=${PREFIX} --enable-cxx
-make -j4 
+make -j${NTASKS}
 make install
 cd ${INITDIR}/local
 touch .built.hdf5

ci-scripts/install-latan.sh

@@ -1,19 +1,23 @@
 #!/usr/bin/env bash
 
 if (( $# != 2 )); then
-  echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
   exit 1
 fi
 PREFIX=$1
-OS=$2
+NTASKS=$2
 
 set -ex
-./install-deps.sh ${PREFIX}
+INITDIR=$(pwd -P)
+mkdir -p ${PREFIX}
+cd ${PREFIX}
+PREFIX=$(pwd -P)
+cd ${INITDIR}
+./install-deps.sh ${PREFIX} ${NTASKS}
 cd ..
 ./bootstrap.sh
 mkdir -p build
 cd build
-if [[ "$OS" == "osx" ]]; then EXTRA_FLAGS='--with-gsl=/usr/local'; fi
+../configure --prefix=${PREFIX} --with-minuit=${PREFIX} --with-nlopt=${PREFIX} --with-hdf5=${PREFIX} --with-gsl=${PREFIX} CXXFLAGS="${CXXFLAGS} -O3 -march=native -mtune=native"
-../configure --prefix=$PREFIX --with-minuit=$PREFIX --with-nlopt=$PREFIX --with-latcore=$PREFIX --with-hdf5=$PREFIX $EXTRA_FLAGS CXXFLAGS='-O3 -march=native -mtune=native'
+make -j${NTASKS}
-make -j4
 make install

ci-scripts/install-latcore.sh

@@ -1,25 +0,0 @@
-#!/usr/bin/env bash
-
-if (( $# != 1 )); then
-  echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
-  exit 1
-fi
-PREFIX=$1
-
-set -ex
-INITDIR=`pwd`
-cd local/build
-if [ -d LatCore ]; then
-  cd LatCore
-  git pull origin master
-else  
-  git clone https://github.com/aportelli/LatCore.git
-  mkdir LatCore/build
-  cd LatCore
-  ./bootstrap.sh
-fi
-cd build
-../configure --prefix=${PREFIX}
-make -j4 
-make install
-cd ${INITDIR}

ci-scripts/install-minuit.sh

@@ -1,22 +1,25 @@
 #!/usr/bin/env bash
 
-NAME='Minuit2-5.34.14'
 
-if (( $# != 1 )); then
+if (( $# != 2 )); then
-  echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
   exit 1
 fi
 PREFIX=$1
+NTASKS=$2
 
 set -ex
-INITDIR=`pwd`
+INITDIR=$(pwd -P)
-cd local/build
+mkdir -p ${PREFIX}
-wget http://www.cern.ch/mathlibs/sw/5_34_14/Minuit2/${NAME}.tar.gz
+cd ${PREFIX}
-tar -xzvf ${NAME}.tar.gz
+PREFIX=$(pwd -P)
-mkdir -p ${NAME}/build
+cd ${INITDIR}/local/build
-cd ${NAME}/build
+rm -rf root
-../configure --prefix=${PREFIX} --disable-openmp
+git clone https://github.com/root-project/root.git
-make -j4 
+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

@@ -1,22 +1,27 @@
 #!/usr/bin/env bash
 
-NAME='nlopt-2.4.2'
+NAME='2.6.1'
 
-if (( $# != 1 )); then
+if (( $# != 2 )); then
-  echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
+  echo "usage: `basename $0` <prefix> <ntasks>" 1>&2
   exit 1
 fi
 PREFIX=$1
+NTASKS=$2
 
 set -ex
-INITDIR=`pwd`
+INITDIR=$(pwd -P)
-cd local/build
+mkdir -p ${PREFIX}
-wget http://ab-initio.mit.edu/nlopt/${NAME}.tar.gz
+cd ${PREFIX}
-tar -xzvf ${NAME}.tar.gz
+PREFIX=$(pwd -P)
+cd ${INITDIR}/local/build
+wget https://github.com/stevengj/nlopt/archive/v${NAME}.tar.gz
+tar -xzvf v${NAME}.tar.gz
+NAME=nlopt-${NAME}
 mkdir -p ${NAME}/build
 cd ${NAME}/build
-../configure --prefix=${PREFIX} --with-cxx --without-guile --without-python --without-octave --without-matlab --with-pic
+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,13 +2,13 @@
 
 # Initialization
 AC_PREREQ([2.63])
-AC_INIT([LatAnalyze],[3.2.1-dev],[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])
 AC_CONFIG_SRCDIR([examples/exMathInterpreter.cpp])
 AC_CONFIG_MACRO_DIR([.buildutils/m4])
-AM_INIT_AUTOMAKE([1.11 -Wall -Werror foreign])
+AM_INIT_AUTOMAKE([1.11 -Wall -Werror foreign subdir-objects])
 AM_SILENT_RULES([yes])
 AC_CONFIG_HEADERS([config.h])
 AM_CONDITIONAL([HAVE_AM_MINOR_LE_11],
@@ -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],
@@ -48,14 +48,6 @@ AC_ARG_WITH([hdf5],
     [try this for a non-standard install prefix of the HDF5 library])],
     [AM_CXXFLAGS="$AM_CXXFLAGS -I$with_hdf5/include"]
     [AM_LDFLAGS="$AM_LDFLAGS -L$with_hdf5/lib"])
-AC_ARG_WITH([latcore],
-    [AS_HELP_STRING([--with-latcore=prefix],
-    [use this option for a non-standard install prefix of the LatCore library])],
-    [AM_CXXFLAGS="$AM_CXXFLAGS -I$with_latcore/include"]
-    [AM_LDFLAGS="$AM_LDFLAGS -L$with_latcore/lib"])
-CFLAGS="$AM_CFLAGS $CFLAGS"
-CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS"
-LDFLAGS="$AM_LDFLAGS $LDFLAGS"
 
 # Get compilers informations
 AX_COMPILER_VENDOR
@@ -78,30 +70,31 @@ AC_DEFINE_UNQUOTED([GXX_VERSION],["$GXX_VERSION"],
 			[version of g++ that will compile the code])
 
 # Checks for libraries.
+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])])
 AC_CHECK_LIB([gsl],[gsl_blas_dgemm],[],[AC_MSG_ERROR([GSL library not found])])
-AC_CHECK_LIB([nlopt_cxx],[nlopt_create],
+AC_CHECK_LIB([nlopt],[nlopt_create],
     [AC_DEFINE([HAVE_NLOPT],
     [1],
-    [Define to 1 if you have the `NLopt' library (-lnlopt_cxx).])]
+    [Define to 1 if you have the `NLopt' library (-lnlopt).])]
     [have_nlopt=true]
-    [LIBS="$LIBS -lnlopt_cxx"],[])
+    [LIBS="$LIBS -lnlopt"],[])
 AM_CONDITIONAL([HAVE_NLOPT], [test x$have_nlopt = xtrue])
 AC_CHECK_LIB([hdf5_cpp],[H5Fopen],
-    [AC_DEFINE([HAVE_HDF5],
+    [LIBS="$LIBS -lhdf5_cpp -lhdf5"],
-    [1],
+    [AC_MSG_ERROR([HDF5 library not found])], [-lhdf5])
-    [Define to 1 if you have the `HDF5' library (-lhdf5_cpp).])]
-    [have_hdf5=true]
-    [LIBS="$LIBS -lhdf5_cpp -lhdf5"],[],[-lhdf5])
-AM_CONDITIONAL([HAVE_HDF5], [test x$have_hdf5 = xtrue])
 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>],
+	[AC_LANG_PROGRAM([#include <Minuit2/VariableMetricMinimizer.h>],
-	[ROOT::Minuit2::BasicMinimumError dummy(0)])],
+	[ROOT::Minuit2::VariableMetricMinimizer dummy()])],
 	[LIBS="$LIBS -lMinuit2"]
 	[AC_DEFINE([HAVE_MINUIT2],
 				[1],
@@ -111,36 +104,59 @@ AC_LINK_IFELSE(
 	[have_minuit=false]
     [AC_MSG_RESULT([no])])
 AM_CONDITIONAL([HAVE_MINUIT], [test x$have_minuit = xtrue])
-LDFLAGS=$SAVED_LDFLAGS
+LDFLAGS="$LDFLAGS -lMinuit2Math"
-SAVED_LDFLAGS=$LDFLAGS
+AC_MSG_CHECKING([for ROOT::Math::MinimizerOptions in -lMinuit2Math]);
-LDFLAGS="$LDFLAGS -lLatCore"
-AC_MSG_CHECKING([for LatCore::XmlReader in -lLatCore]);
 AC_LINK_IFELSE(
-	[AC_LANG_PROGRAM([#include <LatCore/XmlReader.hpp>],
+	[AC_LANG_PROGRAM([#include <Minuit2/Math/MinimizerOptions.h>],
-    [LatCore::XmlReader dummy()])],
+	[ROOT::Math::MinimizerOptions dummy()])],
-	[LIBS="$LIBS -lLatCore"]
+	[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])],
-    [AC_MSG_RESULT([no])]
+	[have_minuitmath=false]
-    [AC_MSG_ERROR([LatCore library not found])])
+    [AC_MSG_RESULT([no])])
+AM_CONDITIONAL([HAVE_MINUITMATH], [test x$have_minuit = xtrue])
 LDFLAGS=$SAVED_LDFLAGS
+CXXFLAGS=$CXXFLAGS_CPY
+LDFLAGS=$LDFLAGS_CPY
 
 # Checks for header files.
 AC_HEADER_STDC
 
+cwd=`pwd -P`; cd ${srcdir}; abs_srcdir=`pwd -P`; cd ${cwd}
+LATAN_CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS"
+LATAN_LDFLAGS="$AM_LDFLAGS $LDFLAGS"
+LATAN_LIBS=$LIBS
+LATAN_SHORT_SHA=`git rev-parse --short HEAD`
+LATAN_SHA=`git rev-parse HEAD`
+LATAN_BRANCH=`git rev-parse --abbrev-ref HEAD`
+AM_CXXFLAGS="-I${abs_srcdir}/lib $AM_CXXFLAGS"
+AM_CFLAGS="-I${abs_srcdir}/lib $AM_CFLAGS"
+CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS"
+LDFLAGS="$AM_LDFLAGS $LDFLAGS"
 AC_SUBST([LIBS])
-AC_SUBST([AM_CFLAGS])
+AC_SUBST([AM_CXXFLAGS])
 AC_SUBST([AM_LDFLAGS])
+AC_SUBST([LATAN_CXXFLAGS])
+AC_SUBST([LATAN_LDFLAGS])
+AC_SUBST([LATAN_LIBS])
+AC_SUBST([LATAN_SHA])
+AC_SUBST([LATAN_BRANCH])
 
-AC_CONFIG_FILES([Makefile lib/Makefile utils/Makefile physics/Makefile 
+AC_CONFIG_FILES([latan-config], [chmod +x latan-config])
-                 examples/Makefile])
+AC_CONFIG_FILES([Makefile])
+AC_CONFIG_FILES([lib/Makefile])
+AC_CONFIG_FILES([utils/Makefile])
+AC_CONFIG_FILES([physics/Makefile])
+AC_CONFIG_FILES([examples/Makefile])
 AC_OUTPUT
 
 echo "*********************************************"
 echo "* ${PACKAGE_NAME} v${VERSION}" build options
 echo "*********************************************"
 echo "* C++ compiler     : $CXX"
-echo "* HDF5 support     : `test x$HAVE_HDF5_TRUE = 'x'   && echo yes \
-                                                          || echo no`"
 echo "* Minuit minimizers: `test x$HAVE_MINUIT_TRUE = 'x' && echo yes \
                                                           || echo no`"
 echo "* NLopt minimizers : `test x$HAVE_NLOPT_TRUE = 'x'  && echo yes \

examples/Makefile.am

@@ -9,6 +9,7 @@ endif
 noinst_PROGRAMS =           \
     exCompiledDoubleFunction\
     exDerivative            \
+    exDWT                   \
     exFit                   \
     exFitSample             \
     exIntegrator            \
@@ -17,8 +18,10 @@ noinst_PROGRAMS =           \
     exMathInterpreter       \
     exMin                   \
     exPlot                  \
+    exPValue                \
     exRand                  \
-    exRootFinder
+    exRootFinder            \
+    exThreadPool
 
 exCompiledDoubleFunction_SOURCES  = exCompiledDoubleFunction.cpp
 exCompiledDoubleFunction_CXXFLAGS = $(COM_CXXFLAGS)
@@ -28,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
@@ -60,6 +67,10 @@ exPlot_SOURCES                    = exPlot.cpp
 exPlot_CXXFLAGS                   = $(COM_CXXFLAGS)
 exPlot_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
 
+exPValue_SOURCES                  = exPValue.cpp
+exPValue_CXXFLAGS                 = $(COM_CXXFLAGS)
+exPValue_LDFLAGS                  = -L../lib/.libs -lLatAnalyze
+
 exRand_SOURCES                    = exRand.cpp
 exRand_CXXFLAGS                   = $(COM_CXXFLAGS)
 exRand_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
@@ -68,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/exCompiledDoubleFunction.cpp

@@ -1,4 +1,4 @@
-#include <LatAnalyze/CompiledFunction.hpp>
+#include <LatAnalyze/Functional/CompiledFunction.hpp>
 
 using namespace std;
 using namespace Latan;

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

@@ -1,6 +1,6 @@
-#include <LatAnalyze/Derivative.hpp>
+#include <LatAnalyze/Numerical/Derivative.hpp>
-#include <LatAnalyze/CompiledFunction.hpp>
+#include <LatAnalyze/Functional/CompiledFunction.hpp>
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
 
 using namespace std;
 using namespace Latan;

examples/exFit.cpp

@@ -1,10 +1,10 @@
 #include <iostream>
 #include <cmath>
-#include <LatAnalyze/CompiledModel.hpp>
+#include <LatAnalyze/Functional/CompiledModel.hpp>
-#include <LatAnalyze/Io.hpp>
+#include <LatAnalyze/Io/Io.hpp>
-#include <LatAnalyze/GslMinimizer.hpp>
+#include <LatAnalyze/Numerical/GslMinimizer.hpp>
-#include <LatAnalyze/Plot.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
-#include <LatAnalyze/XYStatData.hpp>
+#include <LatAnalyze/Statistics/XYStatData.hpp>
 
 using namespace std;
 using namespace Latan;

examples/exFitSample.cpp

@@ -1,7 +1,7 @@
-#include <LatAnalyze/CompiledModel.hpp>
+#include <LatAnalyze/Functional/CompiledModel.hpp>
-#include <LatAnalyze/GslMinimizer.hpp>
+#include <LatAnalyze/Numerical/GslMinimizer.hpp>
-#include <LatAnalyze/Plot.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
-#include <LatAnalyze/XYSampleData.hpp>
+#include <LatAnalyze/Statistics/XYSampleData.hpp>
 
 using namespace std;
 using namespace Latan;

examples/exIntegrator.cpp

@@ -1,5 +1,5 @@
-#include <LatAnalyze/CompiledFunction.hpp>
+#include <LatAnalyze/Functional/CompiledFunction.hpp>
-#include <LatAnalyze/GslQagsIntegrator.hpp>
+#include <LatAnalyze/Numerical/GslQagsIntegrator.hpp>
 
 using namespace std;
 using namespace Latan;

examples/exInterp.cpp

@@ -1,4 +1,4 @@
-#include <LatAnalyze/TabFunction.hpp>
+#include <LatAnalyze/Functional/TabFunction.hpp>
 
 int main(void)
 {

examples/exMat.cpp

@@ -1,6 +1,6 @@
-#include <LatAnalyze/Io.hpp>
+#include <LatAnalyze/Io/Io.hpp>
-#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
 
 using namespace std;
 using namespace Latan;

examples/exMathInterpreter.cpp

@@ -1,5 +1,5 @@
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
-#include <LatAnalyze/MathInterpreter.hpp>
+#include <LatAnalyze/Core/MathInterpreter.hpp>
 
 using namespace std;
 using namespace Latan;

examples/exMin.cpp

@@ -1,5 +1,5 @@
-#include <LatAnalyze/CompiledFunction.hpp>
+#include <LatAnalyze/Functional/CompiledFunction.hpp>
-#include <LatAnalyze/GslMinimizer.hpp>
+#include <LatAnalyze/Numerical/GslMinimizer.hpp>
 
 using namespace std;
 using namespace Latan;

examples/exPValue.cpp

@@ -0,0 +1,23 @@
+#include <LatAnalyze/Core/Math.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char* argv[])
+{
+    double chi2, ndof;
+
+    if (argc != 3)
+    {
+        cerr << "usage: " << argv[0] << " <chi^2> <ndof>" << endl;
+        
+        return EXIT_FAILURE;
+    }
+    chi2 = strTo<double>(argv[1]);
+    ndof = strTo<double>(argv[2]);
+
+    cout << "Two-sided p-value: " << Math::chi2PValue(chi2, ndof) << endl;
+    cout << "chi^2 CCDF       : " << Math::chi2Ccdf(chi2, ndof) << endl;
+
+    return EXIT_SUCCESS;
+}

examples/exPlot.cpp

@@ -1,7 +1,7 @@
-#include <LatAnalyze/CompiledFunction.hpp>
+#include <LatAnalyze/Functional/CompiledFunction.hpp>
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
-#include <LatAnalyze/Plot.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
-#include <LatAnalyze/TabFunction.hpp>
+#include <LatAnalyze/Functional/TabFunction.hpp>
 
 using namespace std;
 using namespace Latan;

examples/exRand.cpp

@@ -1,11 +1,15 @@
-#include <LatAnalyze/Io.hpp>
+#include <LatAnalyze/Io/Io.hpp>
-#include <LatAnalyze/CompiledFunction.hpp>
+#include <LatAnalyze/Functional/CompiledFunction.hpp>
-#include <LatAnalyze/Plot.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
+#include <LatAnalyze/Statistics/Random.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
 
 using namespace std;
 using namespace Latan;
 
+constexpr Index  n             = 8;
 constexpr Index  nDraw         = 20000;
+constexpr Index  nSample       = 2000;
 const     string stateFileName = "exRand.seed";
 
 int main(void)
@@ -36,5 +40,24 @@ int main(void)
     p << PlotFunction(compile("return exp(-x_0^2/2)/sqrt(2*pi);", 1), -5., 5.);
     p.display();
 
+    DMat       var(n, n);
+    DVec       mean(n);
+    DMatSample sample(nSample, n, 1);
+
+    cout << "-- generating " << nSample << " Gaussian random vectors..." << endl;
+    var   = DMat::Random(n, n);
+    var  *= var.adjoint();
+    mean  = DVec::Random(n); 
+    RandomNormal mgauss(mean, var, rd());
+    sample[central] = mgauss();
+    FOR_STAT_ARRAY(sample, s)
+    {
+        sample[s] = mgauss();
+    }
+    cout << "* original variance matrix:\n" << var << endl;
+    cout << "* measured variance matrix:\n" << sample.varianceMatrix() << endl;
+    cout << "* original mean:\n" << mean << endl;
+    cout << "* measured mean:\n" << sample.mean() << endl;
+
     return EXIT_SUCCESS;
 }

examples/exRootFinder.cpp

@@ -1,5 +1,5 @@
-#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
-#include <LatAnalyze/GslHybridRootFinder.hpp>
+#include <LatAnalyze/Numerical/GslHybridRootFinder.hpp>
 
 using namespace std;
 using namespace Latan;
@@ -7,8 +7,8 @@ using namespace Latan;
 int main(void)
 {
     constexpr double a = 1., b = 10.;
-    DoubleFunction f1([a](const double *x){return a*(1.-x[0]);}, 2);
+    DoubleFunction f1([](const double *x){return a*(1.-x[0]);}, 2);
-    DoubleFunction f2([b](const double *x){return b*(x[1]-x[0]*x[0]);}, 2);
+    DoubleFunction f2([](const double *x){return b*(x[1]-x[0]*x[0]);}, 2);
     vector<DoubleFunction *> system = {&f1, &f2};
     GslHybridRootFinder solve;
     DVec init(2), x;

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

latan-config.in

@@ -0,0 +1,79 @@
+#! /bin/sh
+
+prefix=@prefix@
+exec_prefix=@exec_prefix@
+includedir=@includedir@
+
+usage()
+{
+  cat <<EOF
+Usage: latan-config [OPTION]
+
+Known values for OPTION are:
+
+  --prefix     show LatAnalyze installation prefix
+  --cxxflags   print pre-processor and compiler flags
+  --ldflags    print library linking flags
+  --libs       print library linking information
+  --help       display this help and exit
+  --version    output version information
+  --git        print git revision
+
+EOF
+  
+  exit $1
+}
+
+if test $# -eq 0; then
+  usage 1
+fi
+
+cflags=false
+libs=false
+
+while test $# -gt 0; do
+  case "$1" in
+    -*=*) optarg=`echo "$1" | sed 's/[-_a-zA-Z0-9]*=//'` ;;
+    *) optarg= ;;
+  esac
+  
+  case "$1" in
+    --prefix)
+      echo $prefix
+    ;;
+    
+    --version)
+      echo @VERSION@
+      exit 0
+    ;;
+    
+    --git)
+      echo "@LATAN_BRANCH@ @LATAN_SHA@"
+      exit 0
+    ;;
+    
+    --help)
+      usage 0
+    ;;
+    
+    --cxxflags)
+      echo @LATAN_CXXFLAGS@
+    ;;
+    
+    --ldflags)
+      echo @LATAN_LDFLAGS@
+    ;;
+    
+    --libs)
+      echo @LIBS@
+    ;;
+    
+    *)
+      usage
+      exit 1
+    ;;
+  esac
+  shift
+done
+
+exit 0

lib/Core/Eigen.hpp

@@ -0,0 +1,123 @@
+/*
+ * Eigen.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+// Eigen inclusion
+#define EIGEN_DONT_PARALLELIZE
+#define EIGEN_MATRIXBASE_PLUGIN <LatAnalyze/Core/EigenPlugin.hpp>
+#include <LatAnalyze/Eigen/Dense>
+
+// copy/assignement from Eigen expression
+#define EIGEN_EXPR_CTOR(ctorName, Class, Base, ExprType) \
+template <typename Derived>\
+ctorName(const ExprType<Derived> &m): Base(m) {}\
+template<typename Derived>\
+Class & operator=(const ExprType<Derived> &m)\
+{\
+    this->Base::operator=(m);\
+    return *this;\
+}
+
+#define FOR_MAT(mat, i, j) \
+for (Latan::Index j = 0; j < mat.cols(); ++j)\
+for (Latan::Index i = 0; i < mat.rows(); ++i)
+
+BEGIN_LATAN_NAMESPACE
+
+const int dynamic = Eigen::Dynamic;
+
+// array types
+template <typename Derived>
+using ArrayExpr = Eigen::ArrayBase<Derived>;
+
+template <typename T, int nRow = dynamic, int nCol = dynamic>
+using Array = Eigen::Array<T, nRow, nCol>;
+
+// matrix types
+template <typename Derived>
+using MatExpr = Eigen::MatrixBase<Derived>;
+
+template <typename T, int nRow = dynamic, int nCol = dynamic>
+using MatBase = Eigen::Matrix<T, nRow, nCol>;
+
+template <int nRow, int nCol>
+using SFMat = Eigen::Matrix<float, nRow, nCol>;
+
+template <int nRow, int nCol>
+using SDMat = Eigen::Matrix<double, nRow, nCol>;
+
+template <int nRow, int nCol>
+using SCMat = Eigen::Matrix<std::complex<double>, nRow, nCol>;
+
+// vector types
+template <typename T, int size = dynamic>
+using Vec = MatBase<T, size, 1>;
+
+template <int size>
+using SIVec = Vec<int, size>;
+
+template <int size>
+using SUVec = Vec<unsigned int, size>;
+
+template <int size>
+using SFVec = Vec<float, size>;
+
+template <int size>
+using SDVec = Vec<double, size>;
+
+template <int size>
+using SCVec = Vec<std::complex<double>, size>;
+
+typedef SIVec<dynamic> IVec;
+typedef SUVec<dynamic> UVec;
+typedef SDVec<dynamic> DVec;
+typedef SCVec<dynamic> CVec;
+
+// block types
+template <typename Derived>
+using Block      = Eigen::Block<Derived>;
+template <typename Derived>
+using ConstBlock = const Eigen::Block<const Derived>;
+
+template <typename Derived>
+using Row      = typename Derived::RowXpr;
+template <typename Derived>
+using ConstRow = typename Derived::ConstRowXpr;
+
+template <typename Derived>
+using Col = typename Derived::ColXpr;
+template <typename Derived>
+using ConstCol = typename Derived::ConstColXpr;
+
+// map type
+template <int stride>
+using InnerStride = Eigen::InnerStride<stride>;
+template <int rowStride, int colStride>
+using Stride = Eigen::Stride<rowStride, colStride>;
+template <typename Derived, typename StrideType = Stride<0, 0>>
+using Map = Eigen::Map<Derived, Eigen::Unaligned, StrideType>;
+template <typename Derived, typename StrideType = Stride<0, 0>>
+using ConstMap = Eigen::Map<const Derived, Eigen::Unaligned, StrideType>;
+
+// Index type //////////////////////////////////////////////////////////////////
+typedef MatBase<int>::Index Index;
+
+#define FOR_VEC(vec, i)  for (Latan::Index i = 0; i < (vec).size(); ++i)
+#define FOR_ARRAY(ar, i) FOR_VEC(ar, i)
+
+END_LATAN_NAMESPACE

lib/Core/EigenPlugin.hpp

@@ -0,0 +1,60 @@
+/*
+ * EigenPlugin.hpp, part of LatAnalyze
+ *
+ * Copyright (C) 2015 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/>.
+ */
+
+Derived pInverse(const double tolerance = 1.0e-10) const
+{
+    auto         svd   = jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV);
+    const auto   u     = svd.matrixU();
+    const auto   v     = svd.matrixV();
+    auto         s     = svd.singularValues();
+    double       maxsv = 0.;
+    unsigned int elim  = 0;
+    
+    for (Index i = 0; i < s.rows(); ++i)
+    {
+        if (fabs(s(i)) > maxsv) maxsv = fabs(s(i));
+    }
+    for (Index i = 0; i < s.rows(); ++i)
+    {
+        if (fabs(s(i)) > maxsv*tolerance)
+        {
+            s(i) = 1./s(i);
+        }
+        else
+        {
+            elim++;
+            s(i) = 0.;
+        }
+    }
+    if (elim)
+    {
+        std::cerr << "warning: pseudoinverse: " << elim << "/";
+        std::cerr << s.rows() << " singular value(s) eliminated (tolerance= ";
+        std::cerr << tolerance << ")" << std::endl;
+    }
+    
+    return v*s.asDiagonal()*u.transpose();
+}
+
+Derived singularValues(void) const
+{
+    auto svd = jacobiSvd();
+    
+    return svd.singularValues();
+}

lib/Core/Exceptions.cpp

@@ -1,7 +1,7 @@
 /*
  * Exceptions.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Exceptions.hpp>
+#include <LatAnalyze/Core/Exceptions.hpp>
 #include <LatAnalyze/includes.hpp>
 
 #ifndef ERR_SUFF

lib/Core/Exceptions.hpp

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

lib/Core/Mat.cpp

@@ -1,7 +1,7 @@
 /*
  * Mat.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
 #include <LatAnalyze/includes.hpp>
 
 using namespace std;

lib/Core/Mat.hpp

@@ -1,7 +1,7 @@
 /*
  * Mat.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,11 +21,7 @@
 #define Latan_Mat_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/IoObject.hpp>
+#include <LatAnalyze/Io/IoObject.hpp>
-
-#define FOR_MAT(mat, i, j) \
-for (Latan::Index j = 0; j < mat.cols(); ++j)\
-for (Latan::Index i = 0; i < mat.rows(); ++i)
 
 BEGIN_LATAN_NAMESPACE
 

lib/Core/Math.cpp

@@ -1,7 +1,7 @@
 /*
  * Math.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
 #include <LatAnalyze/includes.hpp>
 #include <gsl/gsl_cdf.h>
 
@@ -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                           *
  ******************************************************************************/
@@ -107,8 +130,14 @@ DEF_STD_FUNC_1ARG(fabs)
 
 // p-value
 auto chi2PValueVecFunc = [](const double arg[2])
+{
+    return 2.*min(gsl_cdf_chisq_P(arg[0], arg[1]), gsl_cdf_chisq_Q(arg[0], arg[1]));
+};
+
+auto chi2CcdfVecFunc = [](const double arg[2])
 {
     return gsl_cdf_chisq_Q(arg[0], arg[1]);
 };
 
 DoubleFunction MATH_NAMESPACE::chi2PValue(chi2PValueVecFunc, 2);
+DoubleFunction MATH_NAMESPACE::chi2Ccdf(chi2CcdfVecFunc, 2);

lib/Core/Math.hpp

@@ -1,7 +1,7 @@
 /*
  * Math.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,8 +21,8 @@
 #define	Latan_Math_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
-#include <LatAnalyze/MathInterpreter.hpp>
+#include <LatAnalyze/Core/MathInterpreter.hpp>
 
 BEGIN_LATAN_NAMESPACE
 
@@ -70,10 +70,17 @@ 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
-    const double pi = 3.1415926535897932384626433832795028841970;
+    constexpr double pi  = 3.1415926535897932384626433832795028841970;
-    const double e  = 2.7182818284590452353602874713526624977572;
+    constexpr double e   = 2.7182818284590452353602874713526624977572;
+    constexpr double inf = std::numeric_limits<double>::infinity();
+    constexpr double nan = std::numeric_limits<double>::quiet_NaN();
 }
 
 
@@ -151,6 +158,7 @@ DECL_STD_FUNC(fabs)
 namespace MATH_NAMESPACE
 {
     extern DoubleFunction chi2PValue;
+    extern DoubleFunction chi2Ccdf;
 }
 
 END_LATAN_NAMESPACE

lib/Core/MathInterpreter.cpp

@@ -1,7 +1,7 @@
 /*
  * MathInterpreter.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,9 +17,9 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/MathInterpreter.hpp>
+#include <LatAnalyze/Core/MathInterpreter.hpp>
 #include <LatAnalyze/includes.hpp>
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
 
 using namespace std;
 using namespace Latan;
@@ -217,145 +217,24 @@ void RunContext::reset(void)
 #define CODE_WIDTH 6
 #define CODE_MOD   setw(CODE_WIDTH) << left
 
-// Instruction operator ////////////////////////////////////////////////////////
+auto readConstant(Program::const_iterator ip)
-ostream &Latan::operator<<(ostream& out, const Instruction& ins)
+    -> std::tuple<double, Program::const_iterator>
 {
-    ins.print(out);
+    double value = 0.0;
+    std::copy(ip, ip + sizeof(double), reinterpret_cast<std::uint8_t*>(&value));
 
-    return out;
+    return std::make_tuple(value, ip + sizeof(double));
 }
 
-// Push constructors ///////////////////////////////////////////////////////////
+auto readAddress(Program::const_iterator ip)
-Push::Push(const double val)
+    -> std::tuple<unsigned int, Program::const_iterator>
-: 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
 {
-    if (type_ == ArgType::Constant)
+    unsigned int address = 0.0;
-    {
+    const auto end = ip + sizeof(unsigned int);
-        context.stack().push(val_);
+    std::copy(ip, end, reinterpret_cast<std::uint8_t*>(&address));
-    }
-    else
-    {
-        context.stack().push(context.getVariable(address_));
-    }
-    context.incrementInsIndex();
-}
 
-// Push print //////////////////////////////////////////////////////////////////
+    return std::make_tuple(address, end);
-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, ...)\
+// Bytecode helper functions ///////////////////////////////////////////////////
-program.push_back(unique_ptr<type>(new type(__VA_ARGS__)))
+void pushInstruction(Program &program, const Instruction instruction) {
-#define GET_ADDRESS(address, table, name)\
+    program.push_back(static_cast<std::uint8_t>(instruction));
-try\
+}
-{\
+
-    address = (table).at(name);\
+void pushAddress(Program &program, const unsigned int address) {
-}\
+    const auto address_ptr = reinterpret_cast<const std::uint8_t*>(&address);
-catch (out_of_range)\
+    const auto size = sizeof(unsigned int);
-{\
+    program.insert(program.end(), address_ptr, address_ptr + size);
-    address         = (table).size();\
+}
-    (table)[(name)] = address;\
+
-}\
+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());
+        const unsigned int address = context.addVariable(n[0].getName());
-        if (hasSemicolonParent)
+        pushInstruction(program, Instruction::STORE);
-        {
+        pushAddress(program, address);
-            PUSH_INS(program, Pop, address, n[0].getName());
-        }
-        else
-        {
-            PUSH_INS(program, Store, address, n[0].getName());
-        }
     }
     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,);
+    if (n.getName() == "-" and n.getNArg() == 1) {
-    ELIFNODE("+", 2) PUSH_INS(program, Add,);
+        pushInstruction(program, Instruction::NEG);
-    ELIFNODE("-", 2) PUSH_INS(program, Sub,);
+        return;
-    ELIFNODE("*", 2) PUSH_INS(program, Mul,);
+    }
-    ELIFNODE("/", 2) PUSH_INS(program, Div,);
+
-    ELIFNODE("^", 2) PUSH_INS(program, Pow,);
+    if (getNArg() != 2) {
-    ELSE LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+        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 //////////////////////////////////////////////////////////////
@@ -690,12 +589,12 @@ void MathInterpreter::compile(RunContext &context)
         if (root_)
         {
             context.addVariable("pi",  Math::pi);
-            context.addVariable("inf", HUGE_VAL);
+            context.addVariable("inf", Math::inf);
             ADD_STDMATH_FUNCS(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);
+#define BINARY_OP_CASE(instruction, expr) \
-    while (context.getInsIndex() != program_.size())
+        case Instruction::instruction: { \
-    {
+            const auto second = context.stack().top(); \
-        (*(program_[context.getInsIndex()]))(context);
+            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)
+    return programToStream(out, program.program_);
-    {
-        out << *(program.program_[i]) << endl;
-    }
-    
-    return out;
 }

lib/Core/MathInterpreter.hpp

@@ -1,7 +1,7 @@
 /*
  * MathInterpreter.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -20,9 +20,11 @@
 #ifndef Latan_MathInterpreter_hpp_
 #define	Latan_MathInterpreter_hpp_
 
-#include <LatAnalyze/Function.hpp>
+#include <cstdint>
+
+#include <LatAnalyze/Functional/Function.hpp>
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/ParserState.hpp>
+#include <LatAnalyze/Core/ParserState.hpp>
 
 #define MAXIDLENGTH 256
 
@@ -79,108 +81,26 @@ private:
  *                         Instruction classes                                *
  ******************************************************************************/
 // Abstract base
-class Instruction
+
-{
+enum class Instruction : std::uint8_t {
-public:
+    ADD,
-    // destructor
+    SUB,
-    virtual ~Instruction(void) = default;
+    MUL,
-    // instruction execution
+    DIV,
-    virtual void operator()(RunContext &context) const = 0;
+    POW,
-    friend std::ostream & operator<<(std::ostream &out, const Instruction &ins);
+    NEG,
-private:
+    CONST,
-    virtual void print(std::ostream &out) const = 0;
+    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;
+typedef std::vector<std::uint8_t> 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);
 
 /******************************************************************************
  *                            Expression node classes                         *

lib/Core/MathLexer.lpp

@@ -1,7 +1,7 @@
 /*
  * MathLexer.lpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -25,7 +25,7 @@
 %option yylineno
 
 %{
-    #include <LatAnalyze/MathInterpreter.hpp>
+    #include <LatAnalyze/Core/MathInterpreter.hpp>
     #include "MathParser.hpp"
     
     using namespace std;

lib/Core/MathParser.ypp

@@ -1,7 +1,7 @@
 /*
  * MathParser.ypp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -19,7 +19,7 @@
 
 %{
     #include <LatAnalyze/Global.hpp>
-	#include <LatAnalyze/MathInterpreter.hpp>
+	#include <LatAnalyze/Core/MathInterpreter.hpp>
 
     using namespace std;
     using namespace Latan;

lib/Core/OptParser.cpp

@@ -0,0 +1,299 @@
+/*
+ * OptParser.cpp, part of LatAnalyze
+ *
+ * Copyright (C) 2016 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/Core/OptParser.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+static char optRegex[] = "(-([a-zA-Z])(.+)?)|(--([a-zA-Z_-]+)=?(.+)?)";
+
+/******************************************************************************
+ *                         OptParser implementation                           *
+ ******************************************************************************/
+// regular expressions /////////////////////////////////////////////////////////
+const regex OptParser::optRegex_(optRegex);
+
+// access //////////////////////////////////////////////////////////////////////
+void OptParser::addOption(const std::string shortName,
+                          const std::string longName,
+                          const OptType type, const bool optional,
+                          const std::string helpMessage,
+                          const std::string defaultVal)
+{
+    OptPar par;
+    
+    par.shortName   = shortName;
+    par.longName    = longName;
+    par.defaultVal  = defaultVal;
+    par.helpMessage = helpMessage;
+    par.type        = type;
+    par.optional    = optional;
+    auto it = std::find_if(opt_.begin(), opt_.end(), [&par](const OptPar & p)
+    {
+        bool match = false;
+
+        match |= (par.shortName == p.shortName) and !par.shortName.empty();
+        match |= (par.longName == p.longName) and !par.longName.empty();
+
+        return match;
+    });
+    if (it != opt_.end())
+    {
+        string opt;
+
+        if (!it->shortName.empty())
+        {
+            opt += "-" + it->shortName;
+        }
+        if (!opt.empty())
+        {
+            opt += "/";
+        }
+        if (!it->longName.empty())
+        {
+            opt += "--" + it->longName;
+        }
+        throw(logic_error("duplicate option " + opt + " (in the code, not in the command line)"));
+    }
+    opt_.push_back(par);
+}
+
+bool OptParser::gotOption(const std::string name) const
+{
+    int i = optIndex(name);
+    
+    if (result_.size() != opt_.size())
+    {
+        throw(runtime_error("options not parsed"));
+    }
+    if (i >= 0)
+    {
+        return result_[i].present;
+    }
+    else
+    {
+        throw(out_of_range("no option with name '" + name + "'"));
+    }
+}
+
+const vector<string> & OptParser::getArgs(void) const
+{
+    return arg_;
+}
+
+// parse ///////////////////////////////////////////////////////////////////////
+bool OptParser::parse(int argc, char *argv[])
+{
+    smatch        sm;
+    queue<string> arg;
+    int           expectVal = -1;
+    bool          isCorrect = true;
+    
+    for (int i = 1; i < argc; ++i)
+    {
+        arg.push(argv[i]);
+    }
+    result_.clear();
+    result_.resize(opt_.size());
+    arg_.clear();
+    for (unsigned int i = 0; i < opt_.size(); ++i)
+    {
+        result_[i].value = opt_[i].defaultVal;
+    }
+    while (!arg.empty())
+    {
+        // option
+        if (regex_match(arg.front(), sm, optRegex_))
+        {
+            // should it be a value?
+            if (expectVal >= 0)
+            {
+                cerr << "warning: expected value for option ";
+                cerr << optName(opt_[expectVal]);
+                cerr << ", got option '" << arg.front() << "' instead" << endl;
+                expectVal = -1;
+                isCorrect = false;
+            }
+            // short option
+            if (sm[1].matched)
+            {
+                string optName = sm[2].str();
+                
+                // find option
+                auto it = find_if(opt_.begin(), opt_.end(),
+                                  [&optName](const OptPar &p)
+                                  {
+                                      return (p.shortName == optName);
+                                  });
+                
+                // parse if found
+                if (it != opt_.end())
+                {
+                    unsigned int i =  it - opt_.begin();
+                    
+                    result_[i].present = true;
+                    if (opt_[i].type == OptType::value)
+                    {
+                        if (sm[3].matched)
+                        {
+                            result_[i].value = sm[3].str();
+                        }
+                        else
+                        {
+                            expectVal = i;
+                        }
+                    }
+                }
+                // warning if not found
+                else
+                {
+                    cerr << "warning: unknown option '" << arg.front() << "'";
+                    cerr << endl;
+                }
+            }
+            // long option
+            else if (sm[4].matched)
+            {
+                string optName = sm[5].str();
+                
+                // find option
+                auto it = find_if(opt_.begin(), opt_.end(),
+                                  [&optName](const OptPar &p)
+                                  {
+                                      return (p.longName == optName);
+                                  });
+                
+                // parse if found
+                if (it != opt_.end())
+                {
+                    unsigned int i =  it - opt_.begin();
+                    
+                    result_[i].present = true;
+                    if (opt_[i].type == OptType::value)
+                    {
+                        if (sm[6].matched)
+                        {
+                            result_[i].value = sm[6].str();
+                        }
+                        else
+                        {
+                            expectVal = i;
+                        }
+                    }
+                }
+                // warning if not found
+                else
+                {
+                    cerr << "warning: unknown option '" << arg.front() << "'";
+                    cerr << endl;
+                }
+            }
+        }
+        else if (expectVal >= 0)
+        {
+            result_[expectVal].value = arg.front();
+            expectVal                = -1;
+        }
+        else
+        {
+            arg_.push_back(arg.front());
+        }
+        arg.pop();
+    }
+    if (expectVal >= 0)
+    {
+        cerr << "warning: expected value for option ";
+        cerr << optName(opt_[expectVal]) << endl;
+        expectVal = -1;
+        isCorrect = false;
+    }
+    for (unsigned int i = 0; i < opt_.size(); ++i)
+    {
+        if (!opt_[i].optional and !result_[i].present)
+        {
+            cerr << "warning: mandatory option " << optName(opt_[i]);
+            cerr << " is missing" << endl;
+            isCorrect = false;
+        }
+    }
+    
+    return isCorrect;
+}
+
+// find option index ///////////////////////////////////////////////////////////
+int OptParser::optIndex(const string name) const
+{
+    auto it = find_if(opt_.begin(), opt_.end(), [&name](const OptPar &p)
+                      {
+                          return (p.shortName == name) or (p.longName == name);
+                      });
+    
+    if (it != opt_.end())
+    {
+        return static_cast<int>(it - opt_.begin());
+    }
+    else
+    {
+        return -1;
+    }
+}
+
+// option name for messages ////////////////////////////////////////////////////
+std::string OptParser::optName(const OptPar &opt)
+{
+    std::string res = "";
+    
+    if (!opt.shortName.empty())
+    {
+        res += "-" + opt.shortName;
+        if (!opt.longName.empty())
+        {
+            res += "/";
+        }
+    }
+    if (!opt.longName.empty())
+    {
+        res += "--" + opt.longName;
+        if (opt.type == OptParser::OptType::value)
+        {
+            res += "=";
+        }
+    }
+    
+    return res;
+}
+
+// print option list ///////////////////////////////////////////////////////////
+std::ostream & Latan::operator<<(std::ostream &out, const OptParser &parser)
+{
+    for (auto &o: parser.opt_)
+    {
+        out << setw(20) << OptParser::optName(o);
+        out << ": " << o.helpMessage;
+        if (!o.defaultVal.empty())
+        {
+            out << " (default: " << o.defaultVal << ")";
+        }
+        out << endl;
+    }
+    
+    return out;
+}

lib/Core/OptParser.hpp

@@ -0,0 +1,103 @@
+/*
+ * OptParser.hpp, part of LatAnalyze
+ *
+ * Copyright (C) 2016 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 LatAnalyze_OptParser_hpp_
+#define LatAnalyze_OptParser_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                       command-line option parser                           *
+ ******************************************************************************/
+class OptParser
+{
+public:
+    enum class OptType {value, trigger};
+private:
+    struct OptPar
+    {
+        std::string shortName, longName, defaultVal, helpMessage;
+        OptType     type;
+        bool        optional;
+    };
+    struct OptRes
+    {
+        std::string value;
+        bool        present;
+    };
+public:
+    // constructor
+    OptParser(void) = default;
+    // destructor
+    virtual ~OptParser(void) = default;
+    // access
+    void addOption(const std::string shortName, const std::string longName,
+                   const OptType type, const bool optional = false,
+                   const std::string helpMessage = "",
+                   const std::string defaultVal = "");
+    bool gotOption(const std::string name) const;
+    template <typename T = std::string>
+    T    optionValue(const std::string name) const;
+    const std::vector<std::string> & getArgs(void) const;
+    // parse
+    bool parse(int argc, char *argv[]);
+    // print option list
+    friend std::ostream & operator<<(std::ostream &out,
+                                     const OptParser &parser);
+private:
+    // find option index
+    int optIndex(const std::string name) const;
+    // option name for messages
+    static std::string optName(const OptPar &opt);
+private:
+    std::vector<OptPar>      opt_;
+    std::vector<OptRes>      result_;
+    std::vector<std::string> arg_;
+    static const std::regex  optRegex_;
+};
+
+std::ostream & operator<<(std::ostream &out, const OptParser &parser);
+
+/******************************************************************************
+ *                     OptParser template implementation                      *
+ ******************************************************************************/
+template <typename T>
+T OptParser::optionValue(const std::string name) const
+{
+    int i = optIndex(name);
+    
+    if (result_.size() != opt_.size())
+    {
+        throw(std::runtime_error("options not parsed"));
+    }
+    if (i >= 0)
+    {
+        return strTo<T>(result_[i].value);
+    }
+    else
+    {
+        throw(std::out_of_range("no option with name '" + name + "'"));
+    }
+}
+
+END_LATAN_NAMESPACE
+
+#endif // LatAnalyze_OptParser_hpp_

lib/Core/ParserState.hpp

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

lib/Core/Plot.cpp

@@ -1,7 +1,7 @@
 /*
  * Plot.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,9 +17,9 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Plot.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
 #include <LatAnalyze/includes.hpp>
-#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
 
 using namespace std;
 using namespace Latan;
@@ -112,26 +112,33 @@ 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())
     {
-        LATAN_ERROR(Size, "x and y vector does not have the same size");
+        LATAN_ERROR(Size, "x and y vectors do not have the same size");
     }
 
     DMat d(x[central].rows(), 4);
     string usingCmd, tmpFileName;
 
-    d.col(0)    = x[central];
+    d.col(0)    = x[central].col(0);
-    d.col(2)    = y[central];
+    d.col(2)    = y[central].col(0);
-    d.col(1)    = x.variance().cwiseSqrt();
+    d.col(1)    = x.variance().cwiseSqrt().col(0);
-    d.col(3)    = y.variance().cwiseSqrt();
+    d.col(3)    = y.variance().cwiseSqrt().col(0);
     tmpFileName = dumpToTmpFile(d);
     pushTmpFile(tmpFileName);
+    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,39 +149,151 @@ PlotData::PlotData(const DVec &x, const DMatSample &y)
     string usingCmd, tmpFileName;
 
     d.col(0)    = x;
-    d.col(1)    = y[central];
+    d.col(1)    = y[central].col(0);
-    d.col(2)    = y.variance().cwiseSqrt();
+    d.col(2)    = y.variance().cwiseSqrt().col(0);
     tmpFileName = dumpToTmpFile(d);
     pushTmpFile(tmpFileName);
+    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())
     {
-        LATAN_ERROR(Size, "x and y vector does not have the same size");
+        LATAN_ERROR(Size, "x and y vectors do not have the same size");
     }
 
     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";
+    d(0, 0)     = x;
-    tmpFileName = dumpToTmpFile(data.getTable(i, j));
+    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)
+{
+    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 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 ///////////////////////////////////////////////////////
@@ -199,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;
@@ -212,10 +332,59 @@ PlotFunction::PlotFunction(const DoubleFunction &function, const double xMin,
     }
     tmpFileName = dumpToTmpFile(d);
     pushTmpFile(tmpFileName);
+    if (!abs)
+    {
         setCommand("'" + tmpFileName + "' u 1:2 w lines");
     }
+    else
+    {
+        setCommand("'" + tmpFileName + "' u 1:(abs($2)) w lines");
+    }
+}
 
 // PlotPredBand constructor ////////////////////////////////////////////////////
+void PlotPredBand::makePredBand(const DMat &low, const DMat &high, const double opacity)
+{
+    string lowFileName, highFileName, contFileName;
+    DMat   contour(low.rows() + high.rows() + 1, 2);
+
+    FOR_MAT(low, i, j)
+    {
+        contour(i, j) = low(i, j);
+    }
+    FOR_MAT(high, i, j)
+    {
+        contour(low.rows() + i, j) = high(high.rows() - i - 1, j);
+    }
+    contour.row(low.rows() + high.rows()) = low.row(0);
+    contFileName = dumpToTmpFile(contour);
+    pushTmpFile(contFileName);
+    setCommand("'" + contFileName + "' u 1:2 w filledcurves closed" +
+               " fs solid " + strFrom(opacity) + " noborder");
+}
+
+PlotPredBand::PlotPredBand(const DVec &x, const DVec &y, const DVec &yerr,
+                           const double opacity)
+{
+    if (x.size() != y.size())
+    {
+        LATAN_ERROR(Size, "x and y vectors do not have the same size");
+    }
+    if (y.size() != yerr.size())
+    {
+        LATAN_ERROR(Size, "y and y error vectors do not have the same size");
+    }
+
+    Index nPoint = x.size();
+    DMat  dLow(nPoint, 2), dHigh(nPoint, 2);
+
+    dLow.col(0)  = x;
+    dLow.col(1)  = y - yerr;
+    dHigh.col(0) = x;
+    dHigh.col(1) = y + yerr;
+    makePredBand(dLow, dHigh, opacity);
+}
+
 PlotPredBand::PlotPredBand(const DoubleFunctionSample &function,
                            const double xMin, const double xMax,
                            const unsigned int nPoint, const double opacity)
@@ -236,15 +405,10 @@ PlotPredBand::PlotPredBand(const DoubleFunctionSample &function,
         dHigh(i, 0) = x;
         dHigh(i, 1) = pred[central] + err;
     }
-    lowFileName  = dumpToTmpFile(dLow);
+    makePredBand(dLow, dHigh, opacity);
-    highFileName = dumpToTmpFile(dHigh);
-    pushTmpFile(lowFileName);
-    pushTmpFile(highFileName);
-    setCommand("'< (cat " + lowFileName + "; tac " + highFileName +
-               "; head -n1 " + lowFileName + ")' u 1:2 w filledcurves closed" +
-               " fs solid " + strFrom(opacity) + " noborder");
 }
 
+
 // PlotHistogram constructor ///////////////////////////////////////////////////
 PlotHistogram::PlotHistogram(const Histogram &h)
 {
@@ -261,11 +425,33 @@ PlotHistogram::PlotHistogram(const Histogram &h)
     setCommand("'" + tmpFileName + "' u 1:2 w steps");
 }
 
+// PlotImpulses constructor ////////////////////////////////////////////////////
+PlotImpulses::PlotImpulses(const DVec &x, const DVec &y)
+{
+    if (x.rows() != y.rows())
+    {
+        LATAN_ERROR(Size, "x and y vector does not have the same size");
+    }
+
+    DMat   d(x.rows(), 2);
+    string tmpFileName;
+
+    for (Index i = 0; i < x.rows(); ++i)
+    {
+        d(i, 0) = x(i);
+        d(i, 1) = y(i);
+    }
+    tmpFileName = dumpToTmpFile(d);
+    pushTmpFile(tmpFileName);
+    setCommand("'" + tmpFileName + "' u 1:2 w impulses");
+}
+
 // PlotMatrixNoRange constructor ///////////////////////////////////////////////
 PlotMatrixNoRange::PlotMatrixNoRange(const DMat &m)
 {
     string tmpFileName = dumpToTmpFile(m);
     
+    pushTmpFile(tmpFileName);
     setCommand("'" + tmpFileName + "' matrix w image");
 }
 
@@ -306,6 +492,28 @@ void Color::operator()(PlotOptions &option) const
     option.lineColor = color_;
 }
 
+// LineWidth constructor ///////////////////////////////////////////////////////
+LineWidth::LineWidth(const unsigned int width)
+: width_(width)
+{}
+
+// LineWidth modifier //////////////////////////////////////////////////////////
+void LineWidth::operator()(PlotOptions &option) const
+{
+    option.lineWidth = static_cast<int>(width_);
+}
+
+// Dash constructor ///////////////////////////////////////////////////////////
+Dash::Dash(const string &dash)
+: dash_(dash)
+{}
+
+// Dash modifier //////////////////////////////////////////////////////////////
+void Dash::operator()(PlotOptions &option) const
+{
+    option.dashType = dash_;
+}
+
 // LogScale constructor ////////////////////////////////////////////////////////
 LogScale::LogScale(const Axis axis)
 : axis_(axis)
@@ -371,6 +579,31 @@ void Title::operator()(PlotOptions &option) const
     option.title = title_;
 }
 
+// Palette constructor /////////////////////////////////////////////////////////
+Palette::Palette(const std::vector<std::string> &palette)
+: palette_(palette)
+{}
+
+// Palette modifier ////////////////////////////////////////////////////////////
+void Palette::operator()(PlotOptions &option) const
+{
+    option.palette = palette_;
+}
+
+// category10 palette //////////////////////////////////////////////////////////
+const std::vector<std::string> Palette::category10 =
+{
+    "rgb '#1f77b4'",
+    "rgb '#ff7f0e'",
+    "rgb '#2ca02c'",
+    "rgb '#d62728'",
+    "rgb '#9467bd'",
+    "rgb '#8c564b'",
+    "rgb '#e377c2'",
+    "rgb '#7f7f7f'",
+    "rgb '#bcbd22'"
+};
+
 /******************************************************************************
  *                          Plot implementation                               *
  ******************************************************************************/
@@ -383,7 +616,7 @@ Plot::Plot(void)
 // default options /////////////////////////////////////////////////////////////
 void Plot::initOptions(void)
 {
-    options_.terminal     = "qt";
+    options_.terminal     = "qt noenhanced font 'Arial,12'";
     options_.output       = "";
     options_.caption      = "";
     options_.title        = "";
@@ -394,6 +627,9 @@ void Plot::initOptions(void)
     options_.label[0]     = "";
     options_.label[1]     = "";
     options_.lineColor    = "";
+    options_.lineWidth    = -1;
+    options_.dashType     = "";
+    options_.palette      = Palette::category10;
 }
 
 // plot reset //////////////////////////////////////////////////////////////////
@@ -423,6 +659,16 @@ Plot & Plot::operator<<(PlotObject &&command)
             commandStr         += " lc " + options_.lineColor;
             options_.lineColor  = "";
         }
+        if (options_.lineWidth > 0)
+        {
+            commandStr         += " lw " + strFrom(options_.lineWidth);
+            options_.lineWidth  = -1;
+        }
+        if (!options_.dashType.empty())
+        {
+            commandStr        += " dt " + options_.dashType;
+            options_.dashType  = "";
+        }
         if (options_.title.empty())
         {
             commandStr += " notitle";
@@ -450,22 +696,20 @@ Plot & Plot::operator<<(PlotModifier &&modifier)
 }
 
 // find gnuplot ////////////////////////////////////////////////////////////////
-void Plot::getProgramPath(void)
+#define SEARCH_DIR(dir) \
+sprintf(buf, "%s/%s", dir, gnuplotBin_.c_str());\
+if (access(buf, X_OK) == 0)\
+{\
+    return dir;\
+}
+
+std::string Plot::getProgramPath(void)
 {
     int         i, j, lg;
     char        *path;
     static char buf[MAX_PATH_LENGTH];
     
-    // trivial case: try in CWD
-    sprintf(buf,"./%s", gnuplotBin_.c_str()) ;
-    if (access(buf, X_OK) == 0)
-    {
-        sprintf(buf,".");
-        gnuplotPath_ = buf;
-    }
     // try out in all paths given in the PATH variable
-    else
-    {
     buf[0] = 0;
     path = getenv("PATH") ;
     if (path)
@@ -491,16 +735,9 @@ void Plot::getProgramPath(void)
             if (path[i] == ':') i++ ;
         }
     }
-        else
-        {
-            LATAN_ERROR(System, "PATH variable not set");
-        }
     // if the buffer is still empty, the command was not found
-        if (buf[0] == 0)
+    if (buf[0] != 0)
     {
-            LATAN_ERROR(System, "cannot find gnuplot in $PATH");
-        }
-        // otherwise truncate the command name to yield path only
         lg = (int)(strlen(buf) - 1);
         while (buf[lg]!='/')
         {
@@ -509,7 +746,19 @@ void Plot::getProgramPath(void)
         }
         buf[lg] = 0;
         gnuplotPath_ = buf;
+
+        return gnuplotPath_;
     }
+
+    // try in CWD, /usr/bin & /usr/local/bin
+    SEARCH_DIR(".");
+    SEARCH_DIR("/usr/bin");
+    SEARCH_DIR("/usr/local/bin");
+
+    // if this code is reached nothing was found
+    LATAN_ERROR(System, "cannot find gnuplot");
+
+    return "";
 }
 
 // plot parsing and output /////////////////////////////////////////////////////
@@ -523,12 +772,8 @@ void Plot::display(void)
         string        command;
         ostringstream scriptBuf;
 
-        if (!getenv("DISPLAY"))
-        {
-            LATAN_ERROR(System, "cannot find DISPLAY variable: is it set ?");
-        }
         getProgramPath();
-        command     = gnuplotPath_ + "/" + gnuplotBin_;
+        command     = gnuplotPath_ + "/" + gnuplotBin_ + " 2>/dev/null";
         gnuplotPipe = popen(command.c_str(), "w");
         if (!gnuplotPipe)
         {
@@ -553,7 +798,7 @@ void Plot::display(void)
     }
 }
 
-void Plot::save(string dirName)
+void Plot::save(string dirName, bool savePdf)
 {
     vector<string> commandBack;
     string         path, terminalBack, outputBack, gpCommand, scriptName;
@@ -562,23 +807,26 @@ void Plot::save(string dirName)
     
     mode755 = S_IRWXU|S_IRGRP|S_IXGRP|S_IROTH|S_IXOTH;
 
-    // backup I/O parameters
-    terminalBack = options_.terminal;
-    outputBack   = options_.output;
-    commandBack  = plotCommand_;
-
     // generate directory
     if (mkdir(dirName))
     {
         LATAN_ERROR(Io, "impossible to create directory '" + dirName + "'");
     }
     
+    // backup I/O parameters
+    terminalBack = options_.terminal;
+    outputBack   = options_.output;
+    commandBack  = plotCommand_;
+
     // save PDF
+    if (savePdf)
+    {
         options_.terminal = "pdf";
         options_.output   = dirName + "/plot.pdf";
         display();
         options_.terminal = terminalBack;
         options_.output   = outputBack;
+    }
     
     // save script and datafiles
     for (unsigned int i = 0; i < tmpFileName_.size(); ++i)
@@ -646,18 +894,25 @@ ostream & Latan::operator<<(ostream &out, const Plot &plot)
         out << "yMin = " << plot.options_.scale[y].min << endl;
         out << "yMax = " << plot.options_.scale[y].max << endl;
     }
-    if (!plot.options_.title.empty())
+    out << "unset xrange" << endl;
-    {
-        out << "set title '" << plot.options_.title << "'" << endl;
-    }
     if (plot.options_.scaleMode[x] & Plot::Scale::manual)
     {
         out << "set xrange [xMin:xMax]" << endl;
     }
+    else
+    {
+        out << "set xrange [:]" << endl;
+    }
+    out << "unset yrange" << endl;
     if (plot.options_.scaleMode[y] & Plot::Scale::manual)
     {
         out << "set yrange [yMin:yMax]" << endl;
     }
+    else
+    {
+        out << "set yrange [:]" << endl;
+    }
+    out << "unset log" << endl;
     if (plot.options_.scaleMode[x] & Plot::Scale::log)
     {
         out << "set log x" << endl;
@@ -674,6 +929,11 @@ ostream & Latan::operator<<(ostream &out, const Plot &plot)
     {
         out << "set ylabel '" << plot.options_.label[y] << "'" << endl;
     }
+    for (unsigned int i = 0; i < plot.options_.palette.size(); ++i)
+    {
+        out << "set linetype " << i + 1 << " lc " 
+            << plot.options_.palette[i] << endl;
+    }
     for (unsigned int i = 0; i < plot.headCommand_.size(); ++i)
     {
         out << plot.headCommand_[i] << endl;

lib/Core/Plot.hpp

@@ -1,7 +1,7 @@
 /*
  * Plot.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,11 +21,11 @@
 #define Latan_Plot_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
-#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
-#include <LatAnalyze/MatSample.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
-#include <LatAnalyze/Histogram.hpp>
+#include <LatAnalyze/Statistics/Histogram.hpp>
-#include <LatAnalyze/XYStatData.hpp>
+#include <LatAnalyze/Statistics/XYStatData.hpp>
 
 // gnuplot default parameters
 #ifndef GNUPLOT_BIN
@@ -89,14 +89,27 @@ class PlotData: public PlotObject
 {
 public:
     // constructor
-    PlotData(const DMatSample &x, const DMatSample &y);
+    PlotData(const DMatSample &x, const DMatSample &y, const bool abs = false);
-    PlotData(const DVec       &x, const DMatSample &y);
+    PlotData(const DVec       &x, const DMatSample &y, const bool abs = false);
-    PlotData(const DMatSample &x, const DVec       &y);
+    PlotData(const DMatSample &x, const DVec       &y, const bool abs = false);
-    PlotData(const XYStatData &data, const Index i = 0, const Index j = 0);
+    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:
@@ -106,6 +119,24 @@ public:
     virtual ~PlotHLine(void) = default;
 };
 
+class PlotLine: public PlotObject
+{
+public:
+    // constructor
+    PlotLine(const DVec &x, const DVec &y);
+    // destructor
+    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:
@@ -121,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;
 };
@@ -130,11 +162,15 @@ class PlotPredBand: public PlotObject
 {
 public:
     // constructor
+    PlotPredBand(const DVec &x, const DVec &y, const DVec &yerr,
+                 const double opacity = 0.15);
     PlotPredBand(const DoubleFunctionSample &function, const double xMin,
                  const double xMax, const unsigned int nPoint = 1000,
                  const double opacity = 0.15);
     // destructor
     virtual ~PlotPredBand(void) = default;
+private:
+    void makePredBand(const DMat &low, const DMat &high, const double opacity);
 };
 
 class PlotHistogram: public PlotObject
@@ -146,6 +182,15 @@ public:
     virtual ~PlotHistogram(void) = default;
 };
 
+class PlotImpulses: public PlotObject
+{
+public:
+    // constructor
+    PlotImpulses(const DVec &x, const DVec &y);
+    // destructor
+    virtual ~PlotImpulses(void) = default;
+};
+
 class PlotMatrixNoRange: public PlotObject
 {
 public:
@@ -160,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                                 *
  ******************************************************************************/
@@ -180,6 +230,9 @@ struct PlotOptions
     Range                    scale[2];
     std::string              label[2];
     std::string              lineColor;
+    int                      lineWidth;
+    std::string              dashType;
+    std::vector<std::string> palette;
 };
 
 class PlotModifier
@@ -231,6 +284,32 @@ private:
     const std::string color_;
 };
 
+class LineWidth: public PlotModifier
+{
+public:
+    // constructor
+    explicit LineWidth(const unsigned int width);
+    // destructor
+    virtual ~LineWidth(void) = default;
+    // modifier
+    virtual void operator()(PlotOptions &option) const;
+private:
+    const unsigned width_;
+};
+
+class Dash: public PlotModifier
+{
+public:
+    // constructor
+    explicit Dash(const std::string &dash);
+    // destructor
+    virtual ~Dash(void) = default;
+    // modifier
+    virtual void operator()(PlotOptions &option) const;
+private:
+    const std::string dash_;
+};
+
 class LogScale: public PlotModifier
 {
 public:
@@ -286,6 +365,21 @@ private:
     const std::string title_;
 };
 
+class Palette: public PlotModifier
+{
+public:
+    static const std::vector<std::string> category10;
+public:
+    // constructor
+    explicit Palette(const std::vector<std::string> &palette);
+    // destructor
+    virtual ~Palette(void) = default;
+    // modifier
+    virtual void operator()(PlotOptions &option) const;
+private:
+    const std::vector<std::string> palette_;
+};
+
 /******************************************************************************
  *                               Plot class                                   *
  ******************************************************************************/
@@ -311,13 +405,13 @@ public:
     Plot & operator<<(PlotModifier &&modifier);
     // plot parsing and output
     void display(void);
-    void save(std::string dirName);
+    void save(std::string dirName, bool savePdf = true);
     friend std::ostream & operator<<(std::ostream &out, const Plot &plot);
     // plot reset
     void reset(void);
-private:
     // find gnuplot
-    void getProgramPath(void);
+    std::string getProgramPath(void);
+private:
     // default options
     void initOptions(void);
 private:

lib/Core/ThreadPool.cpp

@@ -0,0 +1,117 @@
+/*
+ * ThreadPool.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2021 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Core/ThreadPool.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                         ThreadPool implementation                          *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+ThreadPool::ThreadPool(void)
+: ThreadPool(std::thread::hardware_concurrency())
+{}
+
+ThreadPool::ThreadPool(const unsigned int nThreads)
+: nThreads_(nThreads)
+{
+    for (unsigned int t = 0; t < nThreads_; ++t)
+    {
+        threads_.push_back(thread(&ThreadPool::workerLoop, this));
+    }
+}
+
+// destructor //////////////////////////////////////////////////////////////////
+ThreadPool::~ThreadPool(void)
+{
+    terminate();
+}
+
+// get the number of threads ///////////////////////////////////////////////////
+unsigned int ThreadPool::getThreadNum(void) const
+{
+    return nThreads_;
+}
+
+// get the pool mutex for synchronisation //////////////////////////////////////
+std::mutex & ThreadPool::getMutex(void)
+{
+    return mutex_;
+}
+
+// worker loop /////////////////////////////////////////////////////////////////
+void ThreadPool::workerLoop(void)
+{
+    while (true)
+    {
+        Job job;
+        {
+            unique_lock<mutex> lock(mutex_);
+
+            condition_.wait(lock, [this](){
+                return !queue_.empty() || terminatePool_;
+            });
+            if (terminatePool_ and queue_.empty())
+            {
+                return;
+            }
+            job = queue_.front();
+            queue_.pop();
+        }
+        job();
+    }
+}
+
+// add jobs ////////////////////////////////////////////////////////////////////
+void ThreadPool::addJob(Job newJob)
+{
+    {
+        unique_lock<mutex> lock(mutex_);
+
+        queue_.push(newJob);
+    }
+    condition_.notify_one();
+}
+
+// critical section ////////////////////////////////////////////////////////////
+void ThreadPool::critical(Job fn)
+{
+    unique_lock<mutex> lock(mutex_);
+
+    fn();
+}
+
+// wait for completion /////////////////////////////////////////////////////////
+void ThreadPool::terminate(void)
+{
+    {
+        unique_lock<mutex> lock(mutex_);
+
+        terminatePool_ = true;
+    }
+    condition_.notify_all();
+    for (auto &thread: threads_)
+    {
+        thread.join();
+    }
+    threads_.clear();
+}

lib/Core/ThreadPool.hpp

@@ -0,0 +1,56 @@
+/*
+ * ThreadPool.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2021 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_ThreadPool_hpp_
+#define Latan_ThreadPool_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+class ThreadPool
+{
+public:
+    typedef std::function<void(void)> Job;
+public:
+    // constructors/destructor
+    ThreadPool(void);
+    ThreadPool(const unsigned int nThreads);
+    virtual ~ThreadPool(void);
+    // get the number of threads
+    unsigned int getThreadNum(void) const;
+    // get the pool mutex for synchronisation
+    std::mutex & getMutex(void);
+    // add jobs
+    void addJob(Job newJob);
+    // critical section
+    void critical(Job fn);
+    // wait for completion and terminate
+    void terminate(void);
+private:
+    // worker loop
+    void workerLoop(void);
+private:
+    unsigned int             nThreads_;
+    std::condition_variable  condition_;
+    std::vector<std::thread> threads_;
+    bool                     terminatePool_{false};
+    std::queue<Job>          queue_;
+    std::mutex               mutex_;
+};
+
+#endif

lib/Core/Utilities.cpp

@@ -0,0 +1,137 @@
+/*
+ * Utilities.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/Core/Utilities.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+void Latan::testFunction(void)
+{}
+
+ostream & Latan::operator<<(ostream &out, const ProgressBar &&bar)
+{
+    const Index nTick = bar.nCol_*bar.current_/bar.total_;
+    
+    out << "[";
+    for (Index i = 0; i < nTick; ++i)
+    {
+        out << "=";
+    }
+    for (Index i = nTick; i < bar.nCol_; ++i)
+    {
+        out << " ";
+    }
+    out << "] " << bar.current_ << "/" << bar.total_;
+    out.flush();
+    
+    return out;
+}
+
+int Latan::mkdir(const std::string dirName)
+{
+    if (access(dirName.c_str(), R_OK|W_OK|X_OK))
+    {
+        mode_t mode755;
+        char   tmp[MAX_PATH_LENGTH];
+        char   *p = NULL;
+        size_t len;
+
+        mode755 = S_IRWXU|S_IRGRP|S_IXGRP|S_IROTH|S_IXOTH;
+
+        snprintf(tmp, sizeof(tmp), "%s", dirName.c_str());
+        len = strlen(tmp);
+        if(tmp[len - 1] == '/')
+        {
+            tmp[len - 1] = 0;
+        }
+        for(p = tmp + 1; *p; p++)
+        {
+            if(*p == '/')
+            {
+                *p = 0;
+                ::mkdir(tmp, mode755);
+                *p = '/';
+            }
+        }
+
+        return ::mkdir(tmp, mode755);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+string Latan::basename(const string &s)
+{
+    constexpr char sep = '/';
+    size_t         i   = s.rfind(sep, s.length());
+    
+    if (i != string::npos)
+    {
+        return s.substr(i+1, s.length() - i);
+    }
+    else
+    {
+        return s;
+    }
+}
+
+std::string Latan::dirname(const std::string &s)
+{
+    constexpr char sep = '/';
+    size_t         i   = s.rfind(sep, s.length());
+    
+    if (i != std::string::npos)
+    {
+        return s.substr(0, i);
+    }
+    else
+    {
+        return "";
+    }
+}
+
+VarName::VarName(const string defName)
+: defName_(defName)
+{}
+
+string VarName::getName(const Index i) const
+{
+    if (hasName(i))
+    {
+        return name_.at(i);
+    }
+    else
+    {
+        return defName_ + "_" + strFrom(i);
+    }
+}
+
+void VarName::setName(const Index i, const string name)
+{
+    name_[i] = name;
+}
+
+bool VarName::hasName(const Index i) const
+{
+    return (name_.find(i) != name_.end());
+}

lib/Core/Utilities.hpp

@@ -0,0 +1,250 @@
+/*
+ * Utilities.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 LatAnalyze_Utilities_hpp_
+#define	LatAnalyze_Utilities_hpp_
+
+#ifndef LATAN_GLOBAL_HPP_
+#include <LatAnalyze/Global.hpp>
+#endif
+
+BEGIN_LATAN_NAMESPACE
+
+// Random seed type ////////////////////////////////////////////////////////////
+typedef std::random_device::result_type SeedType;
+
+// Type utilities //////////////////////////////////////////////////////////////
+// pointer type test
+template <typename Derived, typename Base>
+inline bool isDerivedFrom(const Base *pt)
+{
+    return (dynamic_cast<const Derived *>(pt) != nullptr);
+}
+
+// static logical or
+template <bool... b>
+struct static_or;
+
+template <bool... tail>
+struct static_or<true, tail...> : static_or<tail...> {};
+
+template <bool... tail>
+struct static_or<false, tail...> : std::false_type {};
+
+template <>
+struct static_or<> : std::true_type {};
+
+// Environment /////////////////////////////////////////////////////////////////
+void testFunction(void);
+
+// String utilities ////////////////////////////////////////////////////////////
+inline std::string extension(const std::string fileName)
+{
+    return fileName.substr(fileName.find_last_of(".") + 1);
+}
+
+template <typename T>
+inline T strTo(const std::string &str)
+{
+    T buf;
+    std::istringstream stream(str);
+    
+    stream >> buf;
+    
+    return buf;
+}
+
+// optimized specializations
+template <>
+inline float strTo<float>(const std::string &str)
+{
+    return strtof(str.c_str(), (char **)NULL);
+}
+template <>
+inline double strTo<double>(const std::string &str)
+{
+    return strtod(str.c_str(), (char **)NULL);
+}
+template <>
+inline int strTo<int>(const std::string &str)
+{
+    return (int)(strtol(str.c_str(), (char **)NULL, 10));
+}
+template <>
+inline long strTo<long>(const std::string &str)
+{
+    return strtol(str.c_str(), (char **)NULL, 10);
+}
+template <>
+inline std::string strTo<std::string>(const std::string &str)
+{
+    return str;
+}
+
+template <typename T>
+inline std::string strFrom(const T x)
+{
+    std::ostringstream stream;
+    
+    stream << x;
+    
+    return stream.str();
+}
+
+// specialization for vectors
+template<>
+inline DVec strTo<DVec>(const std::string &str)
+{
+    DVec                res;
+    std::vector<double> vbuf;
+    double              buf;
+    std::istringstream  stream(str);
+    
+    while (!stream.eof())
+    {
+        stream >> buf;
+        vbuf.push_back(buf);
+    }
+    res = Map<DVec>(vbuf.data(), static_cast<Index>(vbuf.size()));
+    
+    return res;
+}
+
+template<>
+inline IVec strTo<IVec>(const std::string &str)
+{
+    IVec                res;
+    std::vector<int>    vbuf;
+    int                 buf;
+    std::istringstream  stream(str);
+    
+    while (!stream.eof())
+    {
+        stream >> buf;
+        vbuf.push_back(buf);
+    }
+    res = Map<IVec>(vbuf.data(), static_cast<Index>(vbuf.size()));
+    
+    return res;
+}
+
+template<>
+inline UVec strTo<UVec>(const std::string &str)
+{
+    UVec                      res;
+    std::vector<unsigned int> vbuf;
+    unsigned int              buf;
+    std::istringstream        stream(str);
+
+    while (!stream.eof())
+    {
+        stream >> buf;
+        vbuf.push_back(buf);
+    }
+    res = Map<UVec>(vbuf.data(), static_cast<Index>(vbuf.size()));
+
+    return res;
+}
+
+template <typename T>
+void tokenReplace(std::string &str, const std::string token,
+                  const T &x, const std::string mark = "@")
+{
+    std::string fullToken = mark + token + mark;
+    
+    auto pos = str.find(fullToken);
+    if (pos != std::string::npos)
+    {
+        str.replace(pos, fullToken.size(), strFrom(x));
+    }
+}
+
+// Manifest file reader ////////////////////////////////////////////////////////
+inline std::vector<std::string> readManifest(const std::string manFileName)
+{
+    std::vector<std::string> list;
+    std::ifstream            manFile;
+    char                     buf[MAX_PATH_LENGTH];
+    
+    manFile.open(manFileName);
+    while (!manFile.eof())
+    {
+        manFile.getline(buf, MAX_PATH_LENGTH);
+        if (!std::string(buf).empty())
+        {
+            list.push_back(buf);
+        }
+    }
+    manFile.close();
+    
+    return list;
+}
+
+// Recursive directory creation ////////////////////////////////////////////////
+int mkdir(const std::string dirName);
+
+// C++ version of basename/dirname /////////////////////////////////////////////
+std::string basename(const std::string& s);
+std::string dirname(const std::string& s);
+
+// Progress bar class //////////////////////////////////////////////////////////
+class ProgressBar
+{
+public:
+    // constructor
+    template <typename A, typename B>
+    ProgressBar(const A current, const B total, const Index nCol = 60);
+    // IO
+    friend std::ostream & operator<<(std::ostream &out,
+                                     const ProgressBar &&bar);
+private:
+    Index current_, total_, nCol_;
+};
+
+template <typename A, typename B>
+ProgressBar::ProgressBar(const A current, const B total, const Index nCol)
+: current_(static_cast<Index>(current))
+, total_(static_cast<Index>(total))
+, nCol_(nCol)
+{}
+
+std::ostream & operator<<(std::ostream &out, const ProgressBar &&bar);
+
+// named variable interface ////////////////////////////////////////////////////
+// FIXME: check redundant names and variable number limit
+class VarName
+{
+public:
+    // constructor
+    VarName(const std::string defName);
+    // destructor
+    virtual ~VarName(void) = default;
+    // access
+    std::string getName(const Index i) const;
+    void        setName(const Index i, const std::string name);
+    // test
+    bool hasName(const Index i) const;
+private:
+    std::string                            defName_;
+    std::unordered_map<Index, std::string> name_;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // LatAnalyze_Utilities_hpp_

lib/Core/stdincludes.hpp

@@ -0,0 +1,58 @@
+/*
+ * stdincludes.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_stdincludes_hpp_
+#define	Latan_stdincludes_hpp_
+
+#include <algorithm>
+#include <array>
+#include <chrono>
+#include <complex>
+#include <condition_variable>
+#include <fstream>
+#include <functional>
+#include <iostream>
+#include <iomanip>
+#include <iterator>
+#include <limits>
+#include <list>
+#include <map>
+#include <memory>
+#include <queue>
+#include <random>
+#include <regex>
+#include <set>
+#include <stack>
+#include <string>
+#include <sstream>
+#include <thread>
+#include <type_traits>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+#include <cfloat>
+#include <climits>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <sys/stat.h>
+#include <unistd.h>
+
+#endif // Latan_stdincludes_hpp_

lib/Functional/CompiledFunction.cpp

@@ -1,7 +1,7 @@
 /*
  * CompiledFunction.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,8 +17,8 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/CompiledFunction.hpp>
+#include <LatAnalyze/Functional/CompiledFunction.hpp>
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
 #include <LatAnalyze/includes.hpp>
 
 using namespace std;

lib/Functional/CompiledFunction.hpp

@@ -1,7 +1,7 @@
 /*
  * CompiledFunction.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,8 +21,8 @@
 #define	Latan_CompiledFunction_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
-#include <LatAnalyze/MathInterpreter.hpp>
+#include <LatAnalyze/Core/MathInterpreter.hpp>
 
 BEGIN_LATAN_NAMESPACE
 

lib/Functional/CompiledModel.cpp

@@ -1,7 +1,7 @@
 /*
  * CompiledModel.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,8 +17,8 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/CompiledModel.hpp>
+#include <LatAnalyze/Functional/CompiledModel.hpp>
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
 #include <LatAnalyze/includes.hpp>
 
 using namespace std;

lib/Functional/CompiledModel.hpp

@@ -1,7 +1,7 @@
 /*
  * CompiledModel.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,8 +21,8 @@
 #define Latan_CompiledModel_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/Model.hpp>
+#include <LatAnalyze/Functional/Model.hpp>
-#include <LatAnalyze/MathInterpreter.hpp>
+#include <LatAnalyze/Core/MathInterpreter.hpp>
 
 BEGIN_LATAN_NAMESPACE
 

lib/Functional/Function.cpp

@@ -1,7 +1,7 @@
 /*
  * Function.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
 #include <LatAnalyze/includes.hpp>
 
 using namespace std;
@@ -111,6 +111,20 @@ double DoubleFunction::operator()(void) const
     return (*this)(nullptr);
 }
 
+std::map<double, double> DoubleFunction::operator()(const std::map<double, double> &m) const
+{
+    checkSize(1);
+
+    std::map<double, double> res;
+
+    for (auto &val: m)
+    {
+        res[val.first] = (*this)(val.second);
+    }
+
+    return res;
+}
+
 // bind ////////////////////////////////////////////////////////////////////////
 DoubleFunction DoubleFunction::bind(const Index argIndex,
                                     const double val) const
@@ -165,6 +179,26 @@ DoubleFunction DoubleFunction::bind(const Index argIndex,
     return bindFunc;
 }
 
+// sample //////////////////////////////////////////////////////////////////////
+DVec DoubleFunction::sample(const DMat &x) const
+{
+    if (x.cols() != getNArg())
+    {
+        LATAN_ERROR(Size, "sampling point matrix and number of arguments "
+                          "mismatch (matrix has " + strFrom(x.cols())
+                          + ", number of arguments is " + strFrom(getNArg()) + ")");
+    }
+    
+    DVec res(x.rows());
+
+    for (Index i = 0; i < res.size(); ++i)
+    {
+        res(i) = (*this)(x.row(i));
+    }
+
+    return res;
+}
+
 // arithmetic operators ////////////////////////////////////////////////////////
 DoubleFunction DoubleFunction::operator-(void) const
 {

lib/Functional/Function.hpp

@@ -1,7 +1,7 @@
 /*
  * Function.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,8 +21,8 @@
 #define	Latan_Function_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
-#include <LatAnalyze/MatSample.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
 
 BEGIN_LATAN_NAMESPACE
 
@@ -52,9 +52,12 @@ public:
     double operator()(void) const;
     template <typename... Ts>
     double operator()(const double arg0, const Ts... args) const;
+    std::map<double, double> operator()(const std::map<double, double> &m) const;
     // bind
     DoubleFunction bind(const Index argIndex, const double val) const;
     DoubleFunction bind(const Index argIndex, const DVec &x) const;
+    // sample
+    DVec sample(const DMat &x) const;
     // arithmetic operators
     DoubleFunction   operator-(void) const;
     DoubleFunction & operator+=(const DoubleFunction &f);

lib/Functional/Model.cpp

@@ -1,7 +1,7 @@
 /*
  * Model.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Model.hpp>
+#include <LatAnalyze/Functional/Model.hpp>
 #include <LatAnalyze/includes.hpp>
 
 using namespace std;

lib/Functional/Model.hpp

@@ -1,7 +1,7 @@
 /*
  * Model.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,8 +21,8 @@
 #define Latan_Model_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
-#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
 
 BEGIN_LATAN_NAMESPACE
 

lib/Functional/TabFunction.cpp

@@ -1,7 +1,7 @@
 /*
  * TabFunction.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/TabFunction.hpp>
+#include <LatAnalyze/Functional/TabFunction.hpp>
 #include <LatAnalyze/includes.hpp>
 
 using namespace std;
@@ -29,18 +29,11 @@ using namespace Latan;
 // constructors ////////////////////////////////////////////////////////////////
 TabFunction::TabFunction(const DVec &x, const DVec &y,
                          const InterpType interpType)
-: interpType_(interpType)
 {
     setData(x, y);
+    setInterpolationType(interpType);
 }
 
-//TabFunction::TabFunction(const XYStatData &data, const Index i, const Index j,
-//                         const InterpType interpType)
-//: interpType_(interpType)
-//{
-//    setData(data, i, j);
-//}
-
 // access //////////////////////////////////////////////////////////////////////
 void TabFunction::setData(const DVec &x, const DVec &y)
 {
@@ -54,10 +47,10 @@ void TabFunction::setData(const DVec &x, const DVec &y)
     }
 }
 
-//void TabFunction::setData(const XYStatData &data, const Index i, const Index j)
+void TabFunction::setInterpolationType(const InterpType interpType)
-//{
+{
-//    setData(data.x(i), data.y(j));
+    interpType_ = interpType;
-//}
+}
 
 // function call ///////////////////////////////////////////////////////////////
 double TabFunction::operator()(const double *arg) const
@@ -154,12 +147,6 @@ DoubleFunction Latan::interpolate(const DVec &x, const DVec &y,
     return TabFunction(x, y, interpType).makeFunction();
 }
 
-//DoubleFunction Latan::interpolate(const XYStatData &data, const Index i,
-//                                  const Index j, const InterpType interpType)
-//{
-//    return TabFunction(data, i, j, interpType).makeFunction();
-//}
-
 map<double, double>::const_iterator TabFunction::nearest(const double x) const
 {
     map<double, double>::const_iterator ret;

lib/Functional/TabFunction.hpp

@@ -1,7 +1,7 @@
 /*
  * TabFunction.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,9 +21,9 @@
 #define Latan_TabFunction_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
-#include <LatAnalyze/XYStatData.hpp>
+#include <LatAnalyze/Statistics/XYStatData.hpp>
 
 BEGIN_LATAN_NAMESPACE
 
@@ -45,13 +45,11 @@ public:
     TabFunction(void) = default;
     TabFunction(const DVec &x, const DVec &y,
                 const InterpType interpType = InterpType::LINEAR);
-    //TabFunction(const XYStatData &data, const Index i = 0, const Index j = 0,
-    //            const InterpType interpType = InterpType::LINEAR);
     // destructor
     virtual ~TabFunction(void) = default;
     // access
     void setData(const DVec &x, const DVec &y);
-    // void setData(const XYStatData &data, const Index i = 0, const Index j = 0);
+    void setInterpolationType(const InterpType interpType);
     // function call
     double operator()(const double *arg) const;
     // factory
@@ -65,9 +63,6 @@ private:
 
 DoubleFunction interpolate(const DVec &x, const DVec &y,
                            const InterpType interpType = InterpType::LINEAR);
-//DoubleFunction interpolate(const XYStatData &data, const Index i = 0,
-//                           const Index j = 0,
-//                           const InterpType interpType = InterpType::LINEAR);
 
 END_LATAN_NAMESPACE
 

lib/Global.cpp

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

lib/Global.hpp

@@ -1,7 +1,7 @@
 /*
  * Global.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -20,13 +20,27 @@
 #ifndef Latan_Global_hpp_
 #define	Latan_Global_hpp_
 
-#include <LatCore/LatCore.hpp>
+#include <LatAnalyze/Core/stdincludes.hpp>
 
 #define BEGIN_LATAN_NAMESPACE \
-namespace Latan {\
+namespace Latan {
-using namespace LatCore;
 #define END_LATAN_NAMESPACE }
 
+// macro utilities
+#define unique_arg(...) __VA_ARGS__
+#define DEBUG_VAR(x) std::cout << #x << "= "  << x << std::endl
+#define DEBUG_MAT(m) std::cout << #m << "=\n" << m << std::endl
+
+// attribute to switch off unused warnings with gcc
+#ifdef __GNUC__
+#define __dumb __attribute__((unused))
+#else
+#define __dumb
+#endif
+
+// max length for paths
+#define MAX_PATH_LENGTH 512u
+
 BEGIN_LATAN_NAMESPACE
 
 // Placeholder type ////////////////////////////////////////////////////////////
@@ -47,6 +61,8 @@ namespace Env
 
 END_LATAN_NAMESPACE
 
-#include <LatAnalyze/Exceptions.hpp>
+#include <LatAnalyze/Core/Eigen.hpp>
+#include <LatAnalyze/Core/Exceptions.hpp>
+#include <LatAnalyze/Core/Utilities.hpp>
 
 #endif // Latan_Global_hpp_

lib/Io/AsciiFile.cpp

@@ -1,7 +1,7 @@
 /*
  * AsciiFile.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/AsciiFile.hpp>
+#include <LatAnalyze/Io/AsciiFile.hpp>
 #include <LatAnalyze/includes.hpp>
 
 using namespace std;

lib/Io/AsciiFile.hpp

@@ -1,7 +1,7 @@
 /*
  * AsciiFile.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,10 +21,10 @@
 #define Latan_AsciiFile_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/File.hpp>
+#include <LatAnalyze/Io/File.hpp>
-#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
-#include <LatAnalyze/MatSample.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
-#include <LatAnalyze/ParserState.hpp>
+#include <LatAnalyze/Core/ParserState.hpp>
 
 BEGIN_LATAN_NAMESPACE
 

lib/Io/AsciiLexer.lpp

@@ -1,7 +1,7 @@
 /*
  * AsciiLexer.lpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -25,7 +25,7 @@
 %option yylineno
 
 %{
-	#include <LatAnalyze/AsciiFile.hpp>
+	#include <LatAnalyze/Io/AsciiFile.hpp>
 	#include "AsciiParser.hpp"
     
     using namespace std;

lib/Io/AsciiParser.ypp

@@ -1,7 +1,7 @@
 /*
  * AsciiParser.ypp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -19,9 +19,9 @@
 
 %{
     #include <LatAnalyze/Global.hpp>
-	#include <LatAnalyze/AsciiFile.hpp>
+	#include <LatAnalyze/Io/AsciiFile.hpp>
-    #include <LatAnalyze/Mat.hpp>
+    #include <LatAnalyze/Core/Mat.hpp>
-    #include <LatAnalyze/MatSample.hpp>
+    #include <LatAnalyze/Statistics/MatSample.hpp>
 
     using namespace std;
     using namespace Latan;

lib/Io/BinReader.cpp

@@ -0,0 +1,81 @@
+/*
+ * BinReader.cpp, part of LatAnalyze
+ *
+ * Copyright (C) 2015 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/Io/BinReader.hpp>
+#include <LatAnalyze/includes.hpp>
+
+#if (defined __GNUC__)||(defined __clang__)
+#pragma GCC diagnostic ignored "-Wunreachable-code"
+#endif
+
+using namespace std;
+using namespace Latan;
+
+BinIO::BinIO(string msg, string loc)
+: runtime_error("Binary reader error: " + msg + " (" + loc + ")")
+{}
+
+/******************************************************************************
+ *                         BinReader implementation                           *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+BinReader::BinReader(const string fileName, const uint32_t endianness,
+                     const bool isColMaj)
+{
+    open(fileName, endianness, isColMaj);
+}
+
+// I/O /////////////////////////////////////////////////////////////////////////
+void BinReader::open(const string fileName, const uint32_t endianness,
+                     const bool isColMaj)
+{
+    fileName_   = fileName;
+    endianness_ = endianness;
+    isColMaj_   = isColMaj;
+    file_.reset(new ifstream(fileName_, ios::in|ios::binary|ios::ate));
+    if (!file_->is_open())
+    {
+        LATAN_ERROR(Io, "impossible to open file '" + fileName_ + "'");
+    }
+    size_ = static_cast<size_t>(file_->tellg());
+    file_->seekg(0, ios::beg);
+}
+
+void BinReader::close(void)
+{
+    file_.reset(nullptr);
+}
+
+template <>
+std::string BinReader::read(void)
+{
+    std::string s;
+    char        c = 'a';
+    
+    while (c != '\n')
+    {
+        c = read<char>();
+        if (c != '\n')
+        {
+            s.push_back(c);
+        }
+    }
+    
+    return s;
+}

lib/Io/BinReader.hpp

@@ -0,0 +1,196 @@
+/*
+ * BinReader.hpp, part of LatAnalyze
+ *
+ * Copyright (C) 2015 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 LatAnalyze_BinReader_hpp_
+#define LatAnalyze_BinReader_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+// I/O exception
+class BinIO: public std::runtime_error
+{
+public:
+    BinIO(std::string msg, std::string loc);
+};
+
+/******************************************************************************
+ *                        Byte manipulation utilities                         *
+ ******************************************************************************/
+class Endianness
+{
+public:
+    enum: uint32_t
+    {
+        little  = 0x00000001,
+        big     = 0x01000000,
+        unknown = 0xffffffff
+    };
+};
+
+class ByteManip
+{
+public:
+    static constexpr uint32_t getHostEndianness(void)
+    {
+        return ((0xffffffff & 1) == Endianness::little) ? Endianness::little
+            : (((0xffffffff & 1) == Endianness::big) ? Endianness::big
+               : Endianness::unknown);
+    }
+    template <typename T>
+    static T swapBytes(const T);
+};
+
+/******************************************************************************
+ *                         template implementation                            *
+ ******************************************************************************/
+template <typename T>
+T ByteManip::swapBytes(const T u)
+{
+    static_assert (CHAR_BIT == 8, "CHAR_BIT != 8");
+    
+    union
+    {
+        T             u;
+        unsigned char u8[sizeof(T)];
+    } source, dest;
+    
+    source.u = u;
+    for (size_t k = 0; k < sizeof(T); ++k)
+    {
+        dest.u8[k] = source.u8[sizeof(T) - k - 1];
+    }
+    
+    return dest.u;
+}
+
+/******************************************************************************
+ *                        Utility to read binary files                        *
+ ******************************************************************************/
+class BinReader
+{
+public:
+    // constructor
+    BinReader(void) = default;
+    BinReader(const std::string fileName,
+              const uint32_t endianness = ByteManip::getHostEndianness(),
+              const bool isColMaj = false);
+    // destructor
+    virtual ~BinReader(void) = default;
+    // I/O
+    void       open(const std::string fileName,
+                    const uint32_t endianness = ByteManip::getHostEndianness(),
+                    const bool isColMaj = false);
+    void       close(void);
+    template <typename T>
+    void       read(T *pt, Index size);
+    template <typename T>
+    T          read(void);
+    template <typename T>
+    MatBase<T> read(const Index nRow, const Index nCol);
+private:
+    std::unique_ptr<std::ifstream> file_{nullptr};
+    std::string                    fileName_;
+    size_t                         size_;
+    uint32_t                       endianness_;
+    bool                           isColMaj_;
+};
+
+/******************************************************************************
+ *                         template implementation                            *
+ ******************************************************************************/
+template <typename T>
+void BinReader::read(T *pt, Index n)
+{
+    if (file_ != nullptr)
+    {
+        file_->read(reinterpret_cast<char *>(pt),
+                    static_cast<long>(sizeof(T))*n);
+        if (endianness_ != ByteManip::getHostEndianness())
+        {
+            for (Index i = 0; i < n; ++i)
+            {
+                pt[i] = ByteManip::swapBytes(pt[i]);
+            }
+        }
+        
+    }
+    else
+    {
+        LATAN_ERROR(Io, "file is not opened");
+    }
+}
+
+template <typename T>
+T BinReader::read(void)
+{
+    T x;
+    
+    if (file_ != nullptr)
+    {
+        read(&x, 1);
+    }
+    else
+    {
+        LATAN_ERROR(Io, "file is not opened");
+    }
+    
+    return x;
+}
+
+template <>
+std::string BinReader::read(void);
+
+template <typename T>
+MatBase<T> BinReader::read(const Index nRow, const Index nCol)
+{
+    MatBase<T> m;
+    
+    // Eigen matrices use column-major ordering
+    if (isColMaj_)
+    {
+        m.resize(nRow, nCol);
+    }
+    else
+    {
+        m.resize(nCol, nRow);
+    }
+    if (file_ != nullptr)
+    {
+        read(m.data(), nRow*nCol);
+    }
+    else
+    {
+        LATAN_ERROR(Io, "file is not opened");
+    }
+    if (isColMaj_)
+    {
+        return m;
+    }
+    else
+    {
+        return m.transpose();
+
+    }
+}
+
+END_LATAN_NAMESPACE
+
+#endif // LatAnalyze_BinReader_hpp_

lib/Io/File.cpp

@@ -1,7 +1,7 @@
 /*
  * File.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/File.hpp>
+#include <LatAnalyze/Io/File.hpp>
 #include <LatAnalyze/includes.hpp>
 
 using namespace std;

lib/Io/File.hpp

@@ -1,7 +1,7 @@
 /*
  * File.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,9 +21,9 @@
 #define	Latan_File_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/IoObject.hpp>
+#include <LatAnalyze/Io/IoObject.hpp>
-#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
-#include <LatAnalyze/MatSample.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
 
 BEGIN_LATAN_NAMESPACE
 

lib/Io/Hdf5File.cpp

@@ -1,7 +1,7 @@
 /*
  * Hdf5File.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
+ * Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,8 +17,8 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Hdf5File.hpp>
+#include <LatAnalyze/Io/Hdf5File.hpp>
-#include <LatAnalyze/IoObject.hpp>
+#include <LatAnalyze/Io/IoObject.hpp>
 #include <LatAnalyze/includes.hpp>
 
 using namespace std;

lib/Io/Hdf5File.hpp

@@ -1,7 +1,7 @@
 /*
  * Hdf5File.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
+ * Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,9 +21,9 @@
 #define Latan_Hdf5File_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/File.hpp>
+#include <LatAnalyze/Io/File.hpp>
-#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
-#include <LatAnalyze/MatSample.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
 #include <H5Cpp.h>
 
 BEGIN_LATAN_NAMESPACE

lib/Io/Io.cpp

@@ -1,7 +1,7 @@
 /*
  * Io.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,12 +17,10 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Io.hpp>
+#include <LatAnalyze/Io/Io.hpp>
 #include <LatAnalyze/includes.hpp>
-#include <LatAnalyze/AsciiFile.hpp>
+#include <LatAnalyze/Io/AsciiFile.hpp>
-#ifdef HAVE_HDF5
+#include <LatAnalyze/Io/Hdf5File.hpp>
-#include <LatAnalyze/Hdf5File.hpp>
-#endif
 
 using namespace std;
 using namespace Latan;
@@ -42,12 +40,10 @@ unique_ptr<File> Io::open(const std::string &fileName, const unsigned int mode)
     {
         return unique_ptr<File>(new AsciiFile(fileName, mode));
     }
-#ifdef HAVE_HDF5
     else if (ext == "h5")
     {
         return unique_ptr<File>(new Hdf5File(fileName, mode));
     }
-#endif
     else
     {
         LATAN_ERROR(Io, "unknown file extension '" + ext + "'");

lib/Io/Io.hpp

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

lib/Io/IoObject.hpp

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

lib/Io/XmlReader.cpp

@@ -0,0 +1,84 @@
+/*
+ * XmlReader.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/Io/XmlReader.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+XmlParsing::XmlParsing(string msg, string loc)
+: runtime_error("XML reader error: " + msg + " (" + loc + ")")
+{}
+
+/******************************************************************************
+ *                         XmlReader implementation                           *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+XmlReader::XmlReader(const string &fileName)
+{
+    open(fileName);
+}
+
+// IO //////////////////////////////////////////////////////////////////////////
+void XmlReader::open(const string &fileName)
+{
+    name_ = fileName;
+    doc_.LoadFile(name_.c_str());
+    if (doc_.Error())
+    {
+        string errMsg;
+        
+        if (doc_.ErrorStr())
+        {
+            errMsg = doc_.ErrorStr();
+        }
+        LATAN_ERROR(Io, "cannot open file " + fileName + " [tinyxml2 code "
+                  + strFrom(doc_.ErrorID()) + ": " + errMsg + "]");
+    }
+    root_ = doc_.RootElement();
+}
+
+// XML structure access ////////////////////////////////////////////////////////
+const XmlNode * XmlReader::getNextNode(const XmlNode *node,
+                                       const string &nodeName)
+{
+    const char *pt = (nodeName.empty()) ? nullptr : nodeName.c_str();
+    
+    if (node)
+    {
+        return node->NextSiblingElement(pt);
+    }
+    else
+    {
+        return nullptr;
+    }
+}
+
+const XmlNode * XmlReader::getNextSameNode(const XmlNode *node)
+{
+    if (node)
+    {
+        return getNextNode(node, node->Name());
+    }
+    else
+    {
+        return nullptr;
+    }
+}

lib/Io/XmlReader.hpp

@@ -0,0 +1,235 @@
+/*
+ * XmlReader.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 LatAnalyze_XmlReader_hpp_
+#define LatAnalyze_XmlReader_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wsign-conversion"
+#include <LatAnalyze/Io/Xml/tinyxml2.hpp>
+#pragma GCC diagnostic pop
+
+BEGIN_LATAN_NAMESPACE
+
+// parsing exception
+class XmlParsing: public std::runtime_error
+{
+public:
+    XmlParsing(std::string msg, std::string loc);
+};
+
+/******************************************************************************
+ *                        XML parameter file reader                           *
+ ******************************************************************************/
+typedef tinyxml2::XMLElement XmlNode;
+
+class XmlReader
+{
+public:
+    // constructor
+    XmlReader(void) = default;
+    explicit XmlReader(const std::string &fileName);
+    // destructor
+    virtual ~XmlReader(void) = default;
+    // IO
+    void open(const std::string &fileName);
+    // XML structure access
+    template <typename... Strs>
+    static const XmlNode * getFirstNode(const XmlNode *startNode,
+                                        const std::string &nodeName,
+                                        Strs... nodeNames);
+    template <typename... Strs>
+    const XmlNode * getFirstNode(const std::string &nodeName,
+                                 Strs... nodeNames) const;
+    static const XmlNode * getNextNode(const XmlNode *node,
+                                       const std::string &nodeName = "");
+    static const XmlNode * getNextSameNode(const XmlNode *node);
+    template <typename T>
+    static T getValue(const XmlNode *node);
+    template <typename T, typename... Strs>
+    static T getFirstValue(const XmlNode *startNode,
+                           const std::string &nodeName, Strs... nodeNames);
+    template <typename T, typename... Strs>
+    T getFirstValue(const std::string &nodeName, Strs... nodeNames) const;
+    template <typename T, typename... Strs>
+    static std::vector<T> getAllValues(const XmlNode *startNode,
+                                       const std::string &nodeName,
+                                       Strs... nodeNames);
+    template <typename T, typename... Strs>
+    std::vector<T> getAllValues(const std::string &nodeName,
+                                Strs... nodeNames) const;
+    // XML structure test
+    template <typename... Strs>
+    static bool hasNode(const XmlNode *startNode, const std::string &nodeName,
+                        Strs... nodeNames);
+    template <typename... Strs>
+    bool hasNode(const std::string &nodeName, Strs... nodeNames) const;
+private:
+    std::string           name_;
+    tinyxml2::XMLDocument doc_;
+    XmlNode               *root_{nullptr};
+};
+
+/******************************************************************************
+ *                     XmlReader template implementation                      *
+ ******************************************************************************/
+// XML structure access ////////////////////////////////////////////////////////
+template <typename... Strs>
+const XmlNode * XmlReader::getFirstNode(const XmlNode *startNode,
+                                        const std::string &nodeName,
+                                        Strs... nodeNames)
+{
+    static_assert(static_or<std::is_assignable<std::string, Strs>::value...>::value,
+                  "getFirstNode arguments are not compatible with std::string");
+    
+    const unsigned int nName  = sizeof...(nodeNames) + 1;
+    const std::string  name[] = {nodeName, nodeNames...};
+    const XmlNode      *node  = startNode;
+    
+    if (!node)
+    {
+        LATAN_ERROR(Io, "root node is null, no XML file opened");
+    }
+    for (unsigned int i = 0; i < nName; ++i)
+    {
+        node = node->FirstChildElement(name[i].c_str());
+        if (!node)
+        {
+            LATAN_ERROR(Io, "XML node " + name[i] + " not found");
+        }
+    }
+    
+    return node;
+}
+
+template <typename... Strs>
+const XmlNode * XmlReader::getFirstNode(const std::string &nodeName,
+                                        Strs... nodeNames) const
+{
+    if (!root_)
+    {
+        LATAN_ERROR(Io, "root node is null, no XML file opened");
+    }
+    
+    return getFirstNode(root_, nodeName, nodeNames...);
+}
+
+template <typename T>
+T XmlReader::getValue(const XmlNode *node)
+{
+    if (node)
+    {
+        if (node->GetText())
+        {
+            return Latan::strTo<T>(node->GetText());
+        }
+        else
+        {
+            return T();
+        }
+    }
+    else
+    {
+        return T();
+    }
+}
+
+template <typename T, typename... Strs>
+T XmlReader::getFirstValue(const XmlNode *startNode,
+                           const std::string &nodeName, Strs... nodeNames)
+{
+    const XmlNode *node = getFirstNode(startNode, nodeName, nodeNames...);
+    
+    return getValue<T>(node);
+}
+
+template <typename T, typename... Strs>
+T XmlReader::getFirstValue(const std::string &nodeName, Strs... nodeNames) const
+{
+    return getFirstValue<T>(root_, nodeName, nodeNames...);
+}
+
+template <typename T, typename... Strs>
+std::vector<T> XmlReader::getAllValues(const XmlNode *startNode,
+                                       const std::string &nodeName,
+                                       Strs... nodeNames)
+{
+    const XmlNode  *node = getFirstNode(startNode, nodeName, nodeNames...);
+    std::vector<T> value;
+    
+    while (node)
+    {
+        value.push_back(getValue<T>(node));
+        node = getNextSameNode(node);
+    }
+    
+    return value;
+}
+
+template <typename T, typename... Strs>
+std::vector<T> XmlReader::getAllValues(const std::string &nodeName,
+                                       Strs... nodeNames) const
+{
+    return getAllValues<T>(root_, nodeName, nodeNames...);
+}
+
+// XML structure test //////////////////////////////////////////////////////////
+template <typename... Strs>
+bool XmlReader::hasNode(const XmlNode *startNode, const std::string &nodeName,
+             Strs... nodeNames)
+{
+    static_assert(static_or<std::is_assignable<std::string, Strs>::value...>::value,
+                  "hasNode arguments are not compatible with std::string");
+
+    const unsigned int nName  = sizeof...(nodeNames) + 1;
+    const std::string  name[] = {nodeName, nodeNames...};
+    const XmlNode      *node  = startNode;
+
+    if (!node)
+    {
+        LATAN_ERROR(Io, "root node is null, no XML file opened");
+    }
+    for (unsigned int i = 0; i < nName; ++i)
+    {
+        node = node->FirstChildElement(name[i].c_str());
+        if (!node)
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+template <typename... Strs>
+bool XmlReader::hasNode(const std::string &nodeName, Strs... nodeNames) const
+{
+    if (!root_)
+    {
+        LATAN_ERROR(Io, "root node is null, no XML file opened");
+    }
+
+    return hasNode(root_, nodeName, nodeNames...);
+}
+
+END_LATAN_NAMESPACE
+
+#endif // LatAnalyze_XmlReader_hpp_

lib/Makefile.am

@@ -7,13 +7,15 @@ if CXX_INTEL
 endif
 endif
 
+include eigen_files.mk
+
 AM_LFLAGS = -olex.yy.c
-AM_YFLAGS = -d
+AM_YFLAGS = -y -d -Wno-yacc -Wno-deprecated
 
 lib_LTLIBRARIES    = libLatAnalyze.la
 noinst_LTLIBRARIES = libLexers.la
 
-libLexers_la_SOURCES  = AsciiLexer.lpp MathLexer.lpp
+libLexers_la_SOURCES  = Io/AsciiLexer.lpp Core/MathLexer.lpp
 if CXX_GNU
     libLexers_la_CXXFLAGS = $(COM_CXXFLAGS) -Wno-unused-parameter -Wno-unused-function -Wno-deprecated-register
 else
@@ -21,100 +23,127 @@ else
 endif
 
 libLatAnalyze_la_SOURCES =           \
-    AsciiFile.cpp          \
-    AsciiParser.ypp        \
-    CompiledFunction.cpp   \
-    CompiledModel.cpp      \
-    Derivative.cpp         \
-    Exceptions.cpp         \
-    File.cpp               \
-    FitInterface.cpp       \
-    Function.cpp           \
-    Global.cpp             \
-    GslHybridRootFinder.cpp\
-    GslMinimizer.cpp       \
-    GslQagsIntegrator.cpp  \
-    Histogram.cpp          \
     includes.hpp                     \
-    Io.cpp                 \
+    Global.cpp                       \
-    Mat.cpp                \
+    Core/Exceptions.cpp              \
-    Math.cpp               \
+    Core/Mat.cpp                     \
-    MathInterpreter.cpp    \
+    Core/Math.cpp                    \
-    MathParser.ypp         \
+    Core/MathInterpreter.cpp         \
-    Minimizer.cpp          \
+    Core/MathParser.ypp              \
-    Model.cpp              \
+    Core/OptParser.cpp               \
-    Plot.cpp               \
+    Core/Plot.cpp                    \
-    RootFinder.cpp         \
+    Core/ThreadPool.cpp              \
-    Solver.cpp             \
+    Core/Utilities.cpp               \
-    StatArray.cpp          \
+    Functional/CompiledFunction.cpp  \
-    TabFunction.cpp        \
+    Functional/CompiledModel.cpp     \
-    XYSampleData.cpp       \
+    Functional/Function.cpp          \
-    XYStatData.cpp         \
+    Functional/Model.cpp             \
+    Functional/TabFunction.cpp       \
+    Io/AsciiFile.cpp                 \
+    Io/AsciiParser.ypp               \
+    Io/BinReader.cpp                 \
+    Io/File.cpp                      \
+    Io/Hdf5File.cpp                  \
+    Io/Io.cpp                        \
+    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       \
+    Numerical/GslQagsIntegrator.cpp  \
+    Numerical/Minimizer.cpp          \
+    Numerical/RootFinder.cpp         \
+    Numerical/Solver.cpp             \
+    Physics/CorrelatorFitter.cpp     \
+    Physics/EffectiveMass.cpp        \
+    Statistics/FitInterface.cpp      \
+    Statistics/Histogram.cpp         \
+    Statistics/Random.cpp            \
+    Statistics/StatArray.cpp         \
+    Statistics/XYSampleData.cpp      \
+    Statistics/XYStatData.cpp        \
 	../config.h
 libLatAnalyze_ladir = $(pkgincludedir)
-libLatAnalyze_la_HEADERS = \
+HPPFILES =                           \
-    AsciiFile.hpp          \
-    CompiledFunction.hpp   \
-    CompiledModel.hpp      \
-    Dataset.hpp            \
-    Derivative.hpp         \
-    Exceptions.hpp         \
-    Function.hpp           \
-    File.hpp               \
-    FitInterface.hpp       \
     Global.hpp                       \
-    GslHybridRootFinder.hpp\
+    Core/Eigen.hpp                   \
-    GslMinimizer.hpp       \
+    Core/EigenPlugin.hpp             \
-    GslQagsIntegrator.hpp  \
+    Core/Exceptions.hpp              \
-    Histogram.hpp          \
+    Core/Mat.hpp                     \
-    Integrator.hpp         \
+    Core/Math.hpp                    \
-    Io.hpp                 \
+    Core/MathInterpreter.hpp         \
-    IoObject.hpp           \
+    Core/OptParser.hpp               \
-    Mat.hpp                \
+    Core/ParserState.hpp             \
-    Math.hpp               \
+    Core/Plot.hpp                    \
-    MathInterpreter.hpp    \
+    Core/ThreadPool.hpp              \
-    MatSample.hpp          \
+    Core/stdincludes.hpp             \
-    Minimizer.hpp          \
+    Core/Utilities.hpp               \
-    Model.hpp              \
+    Functional/CompiledFunction.hpp  \
-    ParserState.hpp        \
+    Functional/CompiledModel.hpp     \
-    Plot.hpp               \
+    Functional/Function.hpp          \
-    RootFinder.hpp         \
+    Functional/Model.hpp             \
-    TabFunction.hpp        \
+    Functional/TabFunction.hpp       \
-    Solver.hpp             \
+    Io/AsciiFile.hpp                 \
-    StatArray.hpp          \
+    Io/BinReader.hpp                 \
-    XYSampleData.hpp       \
+    Io/File.hpp                      \
-    XYStatData.hpp
+    Io/Hdf5File.hpp                  \
-if HAVE_HDF5
+    Io/Io.hpp                        \
-    libLatAnalyze_la_SOURCES += Hdf5File.cpp
+    Io/IoObject.hpp                  \
-    libLatAnalyze_la_HEADERS += Hdf5File.hpp
+    Io/XmlReader.hpp                 \
-endif
+    Numerical/Derivative.hpp         \
+    Numerical/DWT.hpp                \
+    Numerical/DWTFilters.hpp         \
+    Numerical/FFT.hpp                \
+    Numerical/GslFFT.hpp             \
+    Numerical/GslHybridRootFinder.hpp\
+    Numerical/GslMinimizer.hpp       \
+    Numerical/GslQagsIntegrator.hpp  \
+    Numerical/Integrator.hpp         \
+    Numerical/Minimizer.hpp          \
+    Numerical/RootFinder.hpp         \
+    Numerical/Solver.hpp             \
+    Physics/CorrelatorFitter.hpp     \
+    Physics/EffectiveMass.hpp        \
+    Statistics/Dataset.hpp           \
+    Statistics/FitInterface.hpp      \
+    Statistics/Histogram.hpp         \
+    Statistics/MatSample.hpp         \
+    Statistics/Random.hpp            \
+    Statistics/StatArray.hpp         \
+    Statistics/XYSampleData.hpp      \
+    Statistics/XYStatData.hpp
 if HAVE_MINUIT
-    libLatAnalyze_la_SOURCES += MinuitMinimizer.cpp
+    libLatAnalyze_la_SOURCES += Numerical/MinuitMinimizer.cpp
-    libLatAnalyze_la_HEADERS += MinuitMinimizer.hpp
+    HPPFILES += Numerical/MinuitMinimizer.hpp
 endif
 if HAVE_NLOPT
-    libLatAnalyze_la_SOURCES += NloptMinimizer.cpp
+    libLatAnalyze_la_SOURCES += Numerical/NloptMinimizer.cpp
-    libLatAnalyze_la_HEADERS += NloptMinimizer.hpp
+    HPPFILES += Numerical/NloptMinimizer.hpp
 endif
 libLatAnalyze_la_CXXFLAGS = $(COM_CXXFLAGS)
 libLatAnalyze_la_LIBADD   = libLexers.la
 
 if HAVE_AM_MINOR_LE_11
-AsciiParser.hpp: AsciiParser.ypp
+Io/AsciiParser.hpp: Io/AsciiParser.ypp
-	$(AM_V_YACC) $(YACC) -o AsciiParser.cpp --defines=AsciiParser.hpp $<
+	$(AM_V_YACC) $(YACC) -o Io/AsciiParser.cpp --defines=Io/AsciiParser.hpp $<
 
-MathParser.hpp: MathParser.ypp
+Core/MathParser.hpp: Core/MathParser.ypp
-	$(AM_V_YACC) $(YACC) -o MathParser.cpp --defines=MathParser.hpp $<
+	$(AM_V_YACC) $(YACC) -o Core/MathParser.cpp --defines=Core/MathParser.hpp $<
 endif
 
-BUILT_SOURCES = AsciiParser.hpp MathParser.hpp
+BUILT_SOURCES = Io/AsciiParser.hpp Core/MathParser.hpp
 
 CLEANFILES =           \
-    MathLexer.cpp  \
+    Core/MathLexer.cpp \
-    AsciiLexer.cpp \
+    Core/MathParser.cpp\
-    AsciiParser.cpp\
+    Core/MathParser.hpp\
-    AsciiParser.hpp\
+    Io/AsciiLexer.cpp  \
-    MathParser.cpp \
+    Io/AsciiParser.cpp \
-    MathParser.hpp
+    Io/AsciiParser.hpp
+
+nobase_dist_pkginclude_HEADERS = $(HPPFILES) $(eigen_files) Io/Xml/tinyxml2.hpp
 
 ACLOCAL_AMFLAGS = -I .buildutils/m4

lib/Numerical/DWT.cpp

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

lib/Numerical/DWT.hpp

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

lib/Numerical/DWTFilters.cpp

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

lib/Numerical/DWTFilters.hpp

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

lib/Numerical/Derivative.cpp

@@ -1,7 +1,7 @@
 /*
  * Derivative.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,9 +17,9 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Derivative.hpp>
+#include <LatAnalyze/Numerical/Derivative.hpp>
 #include <LatAnalyze/includes.hpp>
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
 
 using namespace std;
 using namespace Latan;

lib/Numerical/Derivative.hpp

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

lib/Numerical/FFT.hpp

@@ -0,0 +1,53 @@
+/*
+ * FFT.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_FFT_hpp_
+#define Latan_FFT_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                             FFT abstract class                             *
+ ******************************************************************************/
+
+class FFT
+{
+public:
+    enum
+    {
+      Forward  = 0,
+      Backward = 1
+    };
+public:
+  // constructor
+  FFT(void) = default;
+  FFT(const Index size);
+  // destructor
+  virtual ~FFT(void) = default;
+  // size
+  virtual void resize(const Index size) = 0;
+  // FFT
+  virtual void operator()(CMat &x, const unsigned int dir = FFT::Forward) = 0;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_FFT_hpp_

lib/Numerical/GslFFT.cpp

@@ -0,0 +1,90 @@
+/*
+ * GslFFT.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/GslFFT.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                           GslFFT implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+GslFFT::GslFFT(const Index size)
+{
+    resize(size);
+}
+
+// destructor //////////////////////////////////////////////////////////////////
+GslFFT::~GslFFT(void)
+{
+    clear();
+}
+
+// size ////////////////////////////////////////////////////////////////////////
+void GslFFT::resize(const Index size)
+{
+    if (size_ != size)
+    {
+        clear();
+        size_      = size;
+        wavetable_ = gsl_fft_complex_wavetable_alloc(size_);
+        workspace_ = gsl_fft_complex_workspace_alloc(size_);
+    }
+}
+
+// fft /////////////////////////////////////////////////////////////////////////
+void GslFFT::operator()(CMat &x, const unsigned int dir)
+{
+    if (x.size() != size_)
+    {
+        LATAN_ERROR(Size, "wrong input vector size");
+    }
+    else
+    {
+        switch (dir)
+        {
+            case FFT::Forward:
+                gsl_fft_complex_forward((double *)x.data(), 1, size_,
+                                        wavetable_, workspace_);
+                break;
+            case FFT::Backward:
+                gsl_fft_complex_backward((double *)x.data(), 1, size_,
+                                         wavetable_, workspace_);
+                break;
+            default:
+                LATAN_ERROR(Argument, "invalid FT direction");
+                break;
+        }
+    }
+}
+
+// destroy GSL objects /////////////////////////////////////////////////////////
+void GslFFT::clear(void)
+{
+    if (!wavetable_)
+    {
+        gsl_fft_complex_wavetable_free(wavetable_);
+    }
+    if (!workspace_)
+    {
+        gsl_fft_complex_workspace_free(workspace_);
+    }
+}

lib/Numerical/GslFFT.hpp

@@ -0,0 +1,57 @@
+/*
+ * GslFFT.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_GslFFT_hpp_
+#define Latan_GslFFT_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
+#include <LatAnalyze/Numerical/FFT.hpp>
+#include <gsl/gsl_fft_complex.h>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                                 GSL FFT                                    *
+ ******************************************************************************/
+
+class GslFFT: public FFT
+{
+public:
+    // constructors
+    GslFFT(void) = default;
+    GslFFT(const Index size);
+    // destructor
+    virtual ~GslFFT(void);
+    // size
+    void resize(const Index size);
+    // fft
+    virtual void operator()(CMat &x, const unsigned int dir = FFT::Forward);
+private:
+    // destroy GSL objects
+    void clear(void);
+private:
+    Index                     size_{0};
+    gsl_fft_complex_wavetable *wavetable_{nullptr};
+    gsl_fft_complex_workspace *workspace_{nullptr};
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_GslFFT_hpp_

lib/Numerical/GslHybridRootFinder.cpp

@@ -1,7 +1,7 @@
 /*
  * GslHybridRootFinder.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,7 +17,7 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/GslHybridRootFinder.hpp>
+#include <LatAnalyze/Numerical/GslHybridRootFinder.hpp>
 #include <LatAnalyze/includes.hpp>
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_multiroots.h>

lib/Numerical/GslHybridRootFinder.hpp

@@ -1,7 +1,7 @@
 /*
  * GslHybridRootFinder.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,7 +21,7 @@
 #define Latan_GslHybridRootFinder_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/RootFinder.hpp>
+#include <LatAnalyze/Numerical/RootFinder.hpp>
 #include <gsl/gsl_multiroots.h>
 
 BEGIN_LATAN_NAMESPACE

lib/Numerical/GslMinimizer.cpp

@@ -1,7 +1,7 @@
 /*
  * GslMinimizer.cpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,9 +17,9 @@
  * along with LatAnalyze.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/GslMinimizer.hpp>
+#include <LatAnalyze/Numerical/GslMinimizer.hpp>
 #include <LatAnalyze/includes.hpp>
-#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Core/Math.hpp>
 #include <gsl/gsl_multimin.h>
 #include <gsl/gsl_blas.h>
 

lib/Numerical/GslMinimizer.hpp

@@ -1,7 +1,7 @@
 /*
  * GslMinimizer.hpp, part of LatAnalyze
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,9 +21,9 @@
 #define Latan_GslMinimizer_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/Derivative.hpp>
+#include <LatAnalyze/Numerical/Derivative.hpp>
-#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
-#include <LatAnalyze/Minimizer.hpp>
+#include <LatAnalyze/Numerical/Minimizer.hpp>
 #include <gsl/gsl_vector.h>
 
 BEGIN_LATAN_NAMESPACE

lib/Numerical/GslQagsIntegrator.cpp

@@ -1,7 +1,7 @@
 /*
  * GslQagsIntegrator.cpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -17,8 +17,9 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/GslQagsIntegrator.hpp>
+#include <LatAnalyze/Numerical/GslQagsIntegrator.hpp>
 #include <LatAnalyze/includes.hpp>
+#include <LatAnalyze/Core/Math.hpp>
 
 using namespace std;
 using namespace Latan;
@@ -55,9 +56,26 @@ double GslQagsIntegrator::operator()(const DoubleFunction &f, const double xMin,
     
     gslF.function = fWrap;
     gslF.params   = reinterpret_cast<void *>(&const_cast<DoubleFunction &>(f));
-    
+    if ((xMin > -Math::inf) and (xMax < Math::inf))
-    gsl_integration_qags(&gslF, xMin, xMax, 0.0, precision_, limit_, workspace_,
+    {
-                         &result, &error_);
+        gsl_integration_qags(&gslF, xMin, xMax, 0.0, precision_, limit_,
+                             workspace_, &result, &error_);
+    }
+    else if (xMax < Math::inf)
+    {
+        gsl_integration_qagil(&gslF, xMax, 0.0, precision_, limit_,
+                              workspace_, &result, &error_);
+    }
+    else if (xMin > -Math::inf)
+    {
+        gsl_integration_qagiu(&gslF, xMin, 0.0, precision_, limit_,
+                              workspace_, &result, &error_);
+    }
+    else
+    {
+        gsl_integration_qagi(&gslF, 0.0, precision_, limit_,
+                             workspace_, &result, &error_);
+    }
     
     return result;
 }

lib/Numerical/GslQagsIntegrator.hpp

@@ -1,7 +1,7 @@
 /*
  * GslQagsIntegrator.hpp, part of LatAnalyze 3
  *
- * Copyright (C) 2013 - 2016 Antonin Portelli
+ * Copyright (C) 2013 - 2020 Antonin Portelli
  *
  * LatAnalyze 3 is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -21,8 +21,8 @@
 #define Latan_GslQagsIntegrator_hpp_
 
 #include <LatAnalyze/Global.hpp>
-#include <LatAnalyze/Function.hpp>
+#include <LatAnalyze/Functional/Function.hpp>
-#include <LatAnalyze/Integrator.hpp>
+#include <LatAnalyze/Numerical/Integrator.hpp>
 #include <gsl/gsl_integration.h>
 
 BEGIN_LATAN_NAMESPACE