proper cmake package config

.github/workflows/build-macos.yml

@@ -1,6 +1,6 @@
 name: Build macOS
 
-on: [push]
+on: [push, workflow_dispatch]
 
 jobs:
   build:

.gitignore

@@ -12,22 +12,19 @@ sandbox/*
 
 # autotools
 autom4te.cache/*
-*.in
-*.in~
-config.h*
 configure
 configure~
 .buildutils/*
 aclocal.m4
 
 # bison/flex generated sources
-lib/*Lexer.cpp
-lib/*Parser.cpp
-lib/*Parser.hpp
+lib/LatAnalyze/*Lexer.cpp
+lib/LatAnalyze/*Parser.cpp
+lib/LatAnalyze/*Parser.hpp
 
-# Eigen headers
-lib/Eigen/*
-lib/eigen_files.mk
+# Eigen headers and archives
+lib/LatAnalyze/Eigen
+eigen-*.tar.bz2
 
 # CI builds
 ci-scripts/local/*

CMakeLists.txt

@@ -0,0 +1,64 @@
+# package config
+cmake_minimum_required(VERSION 3.11.0)
+project(
+  LatAnalyze
+  VERSION 3.6
+  LANGUAGES C CXX)
+
+# includes
+include(CMakePackageConfigHelpers)
+include(FetchContent)
+include(FindPackageMessage)
+include(GNUInstallDirs)
+
+# C++ standard
+set(CMAKE_CXX_STANDARD 14)
+
+# fetch and create symbolic link to Eigen headers
+set(FETCHCONTENT_BASE_DIR ${CMAKE_BINARY_DIR}/deps)
+FetchContent_Declare(
+  Eigen3
+  GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
+  GIT_TAG 3.4.0)
+FetchContent_GetProperties(Eigen3)
+if(NOT eigen3_POPULATED)
+  FetchContent_Populate(Eigen3)
+  message(STATUS "Eigen3 fetched")
+endif()
+file(CREATE_LINK ${eigen3_SOURCE_DIR}/Eigen ${CMAKE_SOURCE_DIR}/lib/LatAnalyze/Eigen
+     SYMBOLIC)
+
+# dependencies
+find_package(GSL REQUIRED)
+find_package(HDF5 REQUIRED COMPONENTS C CXX)
+find_package(Minuit2 QUIET)
+if(Minuit2_FOUND)
+  set(Minuit2_MSG "yes")
+  find_package_message(
+    Minuit2
+    "Found Minuit2: ${Minuit2_DIR} (found version \"${Minuit2_VERSION}\")" " ")
+else()
+  set(Minuit2_MSG "no")
+  message(STATUS "Minuit2 not found")
+endif()
+find_package(NLopt QUIET)
+if(NLopt_FOUND)
+  set(NLopt_MSG "yes")
+  find_package_message(
+    NLopt "Found NLopt: ${NLopt_DIR} (found version \"${NLopt_VERSION}\")" " ")
+else()
+  set(NLopt_MSG "no")
+  message(STATUS "NLopt not found")
+endif()
+
+# subdirectories
+add_subdirectory(lib)
+
+# summary
+message(STATUS "---------------------------")
+message(STATUS "LatAnalyze v${PROJECT_VERSION} configured")
+message(STATUS "---------------------------")
+message(STATUS "    install prefix: ${CMAKE_INSTALL_PREFIX}")
+message(STATUS "Minuit2 minimisers: ${Minuit2_MSG}")
+message(STATUS "  NLopt minimisers: ${NLopt_MSG}")
+message(STATUS "---------------------------")

bootstrap.sh

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

build.sh

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

examples/exRand.cpp

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

lib/CMakeLists.txt

@@ -0,0 +1,81 @@
+# rpath config
+set(CMAKE_BUILD_WITH_INSTALL_RPATH TRUE)
+set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_LIBDIR}")
+
+# config.h header for package version and name
+configure_file(config.h.in ${CMAKE_CURRENT_BINARY_DIR}/config.h @ONLY)
+
+# find bison & flex
+find_package(BISON REQUIRED)
+find_package(FLEX REQUIRED)
+
+# ASCII parser
+bison_target(
+  AsciiParser LatAnalyze/Io/AsciiParser.ypp
+  ${CMAKE_CURRENT_BINARY_DIR}/AsciiParser.cpp COMPILE_FLAGS "-Wno-deprecated")
+flex_target(AsciiLexer LatAnalyze/Io/AsciiLexer.lpp
+            ${CMAKE_CURRENT_BINARY_DIR}/AsciiLexer.cpp)
+add_flex_bison_dependency(AsciiLexer AsciiParser)
+
+# math parser
+bison_target(
+  MathParser LatAnalyze/Core/MathParser.ypp
+  ${CMAKE_CURRENT_BINARY_DIR}/MathParser.cpp COMPILE_FLAGS "-Wno-deprecated")
+flex_target(MathLexer LatAnalyze/Core/MathLexer.lpp
+            ${CMAKE_CURRENT_BINARY_DIR}/MathLexer.cpp)
+add_flex_bison_dependency(MathLexer MathParser)
+
+# library target
+add_library(LatAnalyze SHARED)
+file(GLOB_RECURSE EIGEN_HEADERS LatAnalyze/Eigen/*)
+target_include_directories(LatAnalyze PRIVATE ${CMAKE_CURRENT_BINARY_DIR})
+include(source-list.cmake)
+target_sources(
+  LatAnalyze
+  PRIVATE ${BISON_AsciiParser_OUTPUTS} ${FLEX_AsciiLexer_OUTPUTS}
+          ${BISON_MathParser_OUTPUTS} ${FLEX_MathLexer_OUTPUTS}
+          ${CMAKE_CURRENT_BINARY_DIR}/config.h ${LATAN_SOURCES}
+  PUBLIC FILE_SET
+         public_headers
+         TYPE
+         HEADERS
+         FILES
+         ${EIGEN_HEADERS}
+         ${LATAN_HEADERS})
+target_link_libraries(LatAnalyze PUBLIC GSL::gsl hdf5::hdf5 hdf5::hdf5_cpp)
+if(Minuit2_FOUND)
+  target_link_libraries(LatAnalyze PUBLIC Minuit2::Minuit2)
+endif()
+if(NLopt_FOUND)
+  target_link_libraries(LatAnalyze PUBLIC NLopt::nlopt)
+endif()
+set_target_properties(LatAnalyze PROPERTIES VERSION ${PROJECT_VERSION})
+target_include_directories(
+  LatAnalyze PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+                    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+
+# installation
+install(
+  TARGETS LatAnalyze
+  EXPORT LatAnalyzeTargets
+  LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+  ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+  RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
+          FILE_SET public_headers
+          DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+install(
+  EXPORT LatAnalyzeTargets
+  FILE LatAnalyzeTargets.cmake
+  NAMESPACE LatAnalyze::
+  DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/LatAnalyze)
+configure_package_config_file(
+  ${CMAKE_CURRENT_SOURCE_DIR}/Config.cmake.in
+  "${CMAKE_CURRENT_BINARY_DIR}/LatAnalyzeConfig.cmake"
+  INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/LatAnalyze)
+write_basic_package_version_file(
+  "${CMAKE_CURRENT_BINARY_DIR}/LatAnalyzeConfigVersion.cmake"
+  VERSION "${PROJECT_VERSION}"
+  COMPATIBILITY AnyNewerVersion)
+install(FILES "${CMAKE_CURRENT_BINARY_DIR}/LatAnalyzeConfig.cmake"
+              "${CMAKE_CURRENT_BINARY_DIR}/LatAnalyzeConfigVersion.cmake"
+        DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/LatAnalyze)

lib/Config.cmake.in

@@ -0,0 +1,28 @@
+@PACKAGE_INIT@
+
+include(CMakeFindDependencyMacro)
+
+find_dependency(GSL REQUIRED)
+find_dependency(HDF5 REQUIRED COMPONENTS C CXX)
+set(HAVE_Minuit2 @Minuit2_FOUND@)
+set(HAVE_NLopt @NLopt_FOUND@)
+if(HAVE_Minuit2)
+  find_dependency(Minuit2 REQUIRED)
+  list(APPEND _LatAnalyze_supported_components MINUIT2)
+  set(LatAnalyze_MINUIT2_FOUND True)
+endif()
+if(HAVE_NLopt)
+  find_dependency(NLopt REQUIRED)
+  list(APPEND _LatAnalyze_supported_components NLOPT)
+  set(LatAnalyze_NLOPT_FOUND True)
+endif()
+foreach(_comp ${LatAnalyze_FIND_COMPONENTS})
+  if (NOT ";${_LatAnalyze_supported_components};" MATCHES ";${_comp};")
+    set(LatAnalyze_FOUND False)
+    set(LatAnalyze_NOT_FOUND_MESSAGE "Unsupported LatAnalyze component: ${_comp}")
+  endif()
+endforeach()
+
+include("${CMAKE_CURRENT_LIST_DIR}/LatAnalyzeTargets.cmake")
+
+check_required_components(LatAnalyze)

lib/LatAnalyze

@@ -1 +0,0 @@
-.

lib/LatAnalyze/Core/Math.cpp

@@ -18,6 +18,7 @@
  */
 
 #include <LatAnalyze/Core/Math.hpp>
+#include <LatAnalyze/Numerical/GslFFT.hpp>
 #include <LatAnalyze/includes.hpp>
 #include <gsl/gsl_cdf.h>
 
@@ -48,16 +49,42 @@ DMat MATH_NAMESPACE::corrToVar(const DMat &corr, const DVec &varDiag)
     return res;
 }
 
-double MATH_NAMESPACE::svdDynamicRange(const DMat &mat)
+double MATH_NAMESPACE::conditionNumber(const DMat &mat)
 {
     DVec s = mat.singularValues();
 
     return s.maxCoeff()/s.minCoeff();
 }
 
-double MATH_NAMESPACE::svdDynamicRangeDb(const DMat &mat)
+double MATH_NAMESPACE::cdr(const DMat &mat)
 {
-    return 10.*log10(svdDynamicRange(mat));
+    return 10.*log10(conditionNumber(mat));
+}
+
+template <typename FFT>
+double nsdr(const DMat &m)
+{
+    Index  n = m.rows();
+    FFT    fft(n);
+    CMat   buf(n, 1);
+
+    FOR_VEC(buf, i)
+    {
+        buf(i) = 0.;
+        for (Index j = 0; j < n; ++j)
+        {
+            buf(i) += m(j, (i+j) % n);
+        }
+        buf(i) /= n;
+    }
+    fft(buf, FFT::Forward);
+
+    return 10.*log10(buf.real().maxCoeff()/buf.real().minCoeff());
+}
+
+double MATH_NAMESPACE::nsdr(const DMat &mat)
+{
+    return ::nsdr<GslFFT>(mat);
 }
 
 /******************************************************************************

lib/LatAnalyze/Core/Math.hpp

@@ -73,8 +73,9 @@ namespace MATH_NAMESPACE
     DMat corrToVar(const DMat &corr, const DVec &varDiag);
 
     // matrix SVD dynamic range
-    double svdDynamicRange(const DMat &mat);
-    double svdDynamicRangeDb(const DMat &mat);
+    double conditionNumber(const DMat &mat);
+    double cdr(const DMat &mat);
+    double nsdr(const DMat &mat);
     
     // Constants
     constexpr double pi  = 3.1415926535897932384626433832795028841970;

lib/LatAnalyze/Core/Plot.cpp

@@ -72,7 +72,7 @@ string PlotObject::dumpToTmpFile(const DMat &m)
     for (Index j = 0; j < m.cols(); ++j)
     {
     }
-    sprintf(tmpFileName, "%s/latan_plot_tmp.XXXXXX.dat", P_tmpdir);
+    snprintf(tmpFileName, sizeof(tmpFileName), "%s/latan_plot_tmp.XXXXXX.dat", P_tmpdir);
     fd = mkstemps(tmpFileName, 4);
     if (fd == -1)
     {
@@ -515,14 +515,16 @@ void Dash::operator()(PlotOptions &option) const
 }
 
 // LogScale constructor ////////////////////////////////////////////////////////
-LogScale::LogScale(const Axis axis)
+LogScale::LogScale(const Axis axis, const double basis)
 : axis_(axis)
+, basis_(basis)
 {}
 
 // Logscale modifier ///////////////////////////////////////////////////////////
 void LogScale::operator()(PlotOptions &option) const
 {
-    option.scaleMode[static_cast<int>(axis_)] |= Plot::Scale::log;
+    option.scaleMode[static_cast<int>(axis_)]     |= Plot::Scale::log;
+    option.logScaleBasis[static_cast<int>(axis_)]  = basis_;
 }
 
 // PlotRange constructors //////////////////////////////////////////////////////
@@ -697,7 +699,7 @@ Plot & Plot::operator<<(PlotModifier &&modifier)
 
 // find gnuplot ////////////////////////////////////////////////////////////////
 #define SEARCH_DIR(dir) \
-sprintf(buf, "%s/%s", dir, gnuplotBin_.c_str());\
+snprintf(buf, sizeof(buf), "%s/%s", dir, gnuplotBin_.c_str());\
 if (access(buf, X_OK) == 0)\
 {\
     return dir;\
@@ -915,11 +917,11 @@ ostream & Latan::operator<<(ostream &out, const Plot &plot)
     out << "unset log" << endl;
     if (plot.options_.scaleMode[x] & Plot::Scale::log)
     {
-        out << "set log x" << endl;
+        out << "set log x " << plot.options_.logScaleBasis[x] << endl;;
     }
     if (plot.options_.scaleMode[y] & Plot::Scale::log)
     {
-        out << "set log y" << endl;
+        out << "set log y " << plot.options_.logScaleBasis[y] << endl;
     }
     if (!plot.options_.label[x].empty())
     {

lib/LatAnalyze/Core/Plot.hpp

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

lib/LatAnalyze/Core/Utilities.hpp

@@ -108,23 +108,6 @@ inline std::string strFrom(const T x)
 }
 
 // specialization for vectors
-template<>
-inline std::vector<Index> strTo<std::vector<Index>>(const std::string &str)
-{
-    std::vector<Index>  res;
-    std::vector<double> vbuf;
-    double              buf;
-    std::istringstream  stream(str);
-    
-    while (!stream.eof())
-    {
-        stream >> buf;
-        res.push_back(buf);
-    }
-    
-    return res;
-}
-
 template<>
 inline DVec strTo<DVec>(const std::string &str)
 {

lib/LatAnalyze/Global.cpp

@@ -25,11 +25,10 @@ using namespace Latan;
 
 PlaceHolder Latan::_;
 
-const string Env::fullName  = PACKAGE_STRING;
-const string Env::name      = PACKAGE_NAME;
-const string Env::version   = PACKAGE_VERSION;
-const string Env::msgPrefix = "[" + strFrom(PACKAGE_NAME) + " v"
-                              + strFrom(PACKAGE_VERSION) + "] ";
+const string Env::fullName  = strFrom(PROJECT_NAME) + " v" + strFrom(PROJECT_VERSION);
+const string Env::name      = PROJECT_NAME;
+const string Env::version   = PROJECT_VERSION;
+const string Env::msgPrefix = "[" + Env::fullName + "] ";
 
 void Env::function(void)
 {}

lib/LatAnalyze/Makefile.am

@@ -58,6 +58,7 @@ libLatAnalyze_la_SOURCES =           \
     Numerical/RootFinder.cpp         \
     Numerical/Solver.cpp             \
     Physics/CorrelatorFitter.cpp     \
+    Physics/DataFilter.cpp           \
     Physics/EffectiveMass.cpp        \
     Statistics/FitInterface.cpp      \
     Statistics/Histogram.cpp         \
@@ -106,6 +107,7 @@ HPPFILES =                           \
     Numerical/RootFinder.hpp         \
     Numerical/Solver.hpp             \
     Physics/CorrelatorFitter.hpp     \
+    Physics/DataFilter.hpp           \
     Physics/EffectiveMass.hpp        \
     Statistics/Dataset.hpp           \
     Statistics/FitInterface.hpp      \

lib/LatAnalyze/Numerical/DWT.cpp

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

lib/LatAnalyze/Numerical/DWT.hpp

@@ -22,6 +22,7 @@
 
 #include <LatAnalyze/Global.hpp>
 #include <LatAnalyze/Numerical/DWTFilters.hpp>
+#include <LatAnalyze/Core/Mat.hpp>
 
 BEGIN_LATAN_NAMESPACE
 
@@ -40,12 +41,18 @@ public:
     // convolution primitive
     static void filterConvolution(DVec &out, const DVec &data, 
                                   const std::vector<double> &filter, const Index offset);
+    static void filterConvolution(DMat &out, const DMat &data, 
+                                  const std::vector<double> &filter, const Index offset);
     // downsampling/upsampling primitives
     static void downsample(DVec &out, const DVec &in);
+    static void downsample(DMat &out, const DMat &in);
     static void upsample(DVec &out, const DVec &in);
+    static void upsample(DMat &out, const DMat &in);
     // DWT
     std::vector<DWTLevel> forward(const DVec &data, const unsigned int level) const;
     DVec                  backward(const std::vector<DWTLevel>& dwt) const;
+    // concatenate levels
+    static DVec concat(const std::vector<DWTLevel>& dwt, const int maxLevel = -1, const bool dropLow = false);
 private:
     DWTFilter filter_;
 };

lib/LatAnalyze/Numerical/MinuitMinimizer.cpp

@@ -30,8 +30,12 @@ namespace ROOT
 };
 
 // macros necessary in the ROOT-based version of Minuit2
+#ifndef ROOT_Math_VecTypes
 #define ROOT_Math_VecTypes
+#endif
+#ifndef MATHCORE_STANDALONE
 #define MATHCORE_STANDALONE
+#endif
 
 #include <Minuit2/Minuit2Minimizer.h>
 #include <Math/Functor.h>

lib/LatAnalyze/Physics/CorrelatorFitter.cpp

@@ -253,39 +253,16 @@ DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts)
     }
 }
 
-DMatSample CorrelatorUtils::fold(const DMatSample &c, const CorrelatorModels::ModelPar &par)
+DMatSample CorrelatorUtils::fold(const DMatSample &c)
 {
     const Index nt  = c[central].rows();
     DMatSample  buf = c;
-    int         sign;
-    bool        fold = false;
-    
-    switch (par.type)
-    {
-        case CorrelatorType::cosh:
-        case CorrelatorType::cst:
-            sign = 1;
-            fold = true;
-            break;
-        case CorrelatorType::sinh:
-            sign = -1;
-            fold = true;
-            break;
-        case CorrelatorType::linear:
-            cout << "Linear model is asymmetric: will not fold." << endl;
-            break;
-        default:
-            break;
-    }
 
-    if (fold)
+    FOR_STAT_ARRAY(buf, s)
     {
-        FOR_STAT_ARRAY(buf, s)
+        for (Index t = 0; t < nt; ++t)
         {
-            for (Index t = 0; t < nt; ++t)
-            {
-                buf[s](t) = 0.5*(c[s](t) + sign*c[s]((nt - t) % nt));
-            }
+            buf[s](t) = 0.5*(c[s](t) + c[s]((nt - t) % nt));
         }
     }
 

lib/LatAnalyze/Physics/CorrelatorFitter.hpp

@@ -56,7 +56,7 @@ namespace CorrelatorModels
 namespace CorrelatorUtils
 {
     DMatSample shift(const DMatSample &c, const Index ts);
-    DMatSample fold(const DMatSample &c, const CorrelatorModels::ModelPar &par);
+    DMatSample fold(const DMatSample &c);
     DMatSample fourierTransform(const DMatSample &c, FFT &fft, 
                                 const unsigned int dir = FFT::Forward);
 };

lib/LatAnalyze/Physics/DataFilter.cpp

@@ -0,0 +1,83 @@
+/*
+ * DataFilter.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Physics/DataFilter.hpp>
+#include <LatAnalyze/includes.hpp>
+#include <LatAnalyze/Numerical/DWT.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                        DataFilter implementation                           *
+ ******************************************************************************/
+// constructor ////////////////////////////////////////////////////////////////
+DataFilter::DataFilter(const vector<double> &filter, const bool downsample)
+: filter_(filter), downsample_(downsample)
+{}
+
+// filtering //////////////////////////////////////////////////////////////////
+template <typename MatType>
+void filter(MatType &out, const MatType &in, const vector<double> &filter, 
+            const bool downsample, MatType &buf)
+{
+    if (!downsample)
+    {
+        out.resizeLike(in);
+        DWT::filterConvolution(out, in, filter, filter.size()/2);
+    }
+    else
+    {
+        out.resize(in.rows()/2, in.cols());
+        buf.resizeLike(in);
+        DWT::filterConvolution(buf, in, filter, filter.size()/2);
+        DWT::downsample(out, buf);
+    }
+}
+
+void DataFilter::operator()(DVec &out, const DVec &in)
+{
+    filter(out, in, filter_, downsample_, vBuf_);
+}
+
+void DataFilter::operator()(DMat &out, const DMat &in)
+{
+    filter(out, in, filter_, downsample_, mBuf_);
+}
+
+/******************************************************************************
+ *                     LaplaceDataFilter implementation                       *
+ ******************************************************************************/
+// constructor ////////////////////////////////////////////////////////////////
+LaplaceDataFilter::LaplaceDataFilter(const bool downsample)
+: DataFilter({1., -2. , 1.}, downsample)
+{}
+
+// filtering //////////////////////////////////////////////////////////////////
+void LaplaceDataFilter::operator()(DVec &out, const DVec &in, const double lambda)
+{
+    filter_[1] = -2. - lambda;
+    DataFilter::operator()(out, in);
+}
+
+void LaplaceDataFilter::operator()(DMat &out, const DMat &in, const double lambda)
+{
+    filter_[1] = -2. - lambda;
+    DataFilter::operator()(out, in);
+}

lib/LatAnalyze/Physics/DataFilter.hpp

@@ -0,0 +1,139 @@
+/*
+ * DataFilter.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_DataFilter_hpp_
+#define Latan_DataFilter_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Core/Math.hpp>
+#include <LatAnalyze/Statistics/StatArray.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
+#include <LatAnalyze/Numerical/Minimizer.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                    Generic convolution filter class                        *
+ ******************************************************************************/
+class DataFilter
+{
+public:
+    // constructor
+    DataFilter(const std::vector<double> &filter, const bool downsample = false);
+    // filtering
+    void operator()(DVec &out, const DVec &in);
+    void operator()(DMat &out, const DMat &in);
+    template <typename MatType, Index o>
+    void operator()(StatArray<MatType, o> &out, const StatArray<MatType, o> &in);
+protected:
+    std::vector<double> filter_;
+private:
+    bool                downsample_;
+    DVec                vBuf_;
+    DMat                mBuf_;
+};
+
+/******************************************************************************
+ *                         Laplacian filter class                             *
+ ******************************************************************************/
+class LaplaceDataFilter: public DataFilter
+{
+public:
+    // constructor
+    LaplaceDataFilter(const bool downsample = false);
+    // filtering
+    void operator()(DVec &out, const DVec &in, const double lambda = 0.);
+    void operator()(DMat &out, const DMat &in, const double lambda = 0.);
+    template <typename MatType, Index o>
+    void operator()(StatArray<MatType, o> &out, const StatArray<MatType, o> &in, 
+                    const double lambda = 0.);
+    // correlation optimisation
+    template <typename MatType, Index o>
+    double optimiseCdr(const StatArray<MatType, o> &data, Minimizer &min,
+                       const unsigned int nPass = 3);
+};
+
+/******************************************************************************
+ *                 DataFilter class template implementation                   *
+ ******************************************************************************/
+// filtering //////////////////////////////////////////////////////////////////
+template <typename MatType, Index o>
+void DataFilter::operator()(StatArray<MatType, o> &out, const StatArray<MatType, o> &in)
+{
+    FOR_STAT_ARRAY(in, s)
+    {
+        (*this)(out[s], in[s]);
+    }
+}
+
+/******************************************************************************
+ *            LaplaceDataFilter class template implementation                 *
+ ******************************************************************************/
+// filtering //////////////////////////////////////////////////////////////////
+template <typename MatType, Index o>
+void LaplaceDataFilter::operator()(StatArray<MatType, o> &out, 
+                                   const StatArray<MatType, o> &in, const double lambda)
+{
+    FOR_STAT_ARRAY(in, s)
+    {
+        (*this)(out[s], in[s], lambda);
+    }
+}
+
+// correlation optimisation ///////////////////////////////////////////////////
+template <typename MatType, Index o>
+double LaplaceDataFilter::optimiseCdr(const StatArray<MatType, o> &data, 
+                                      Minimizer &min, const unsigned int nPass)
+{
+    StatArray<MatType, o> fdata(data.size());
+    DVec init(1);
+    double reg, prec;
+    DoubleFunction cdr([&data, &fdata, this](const double *x)
+    {
+        double res;
+        (*this)(fdata, data, x[0]);
+        res = Math::cdr(fdata.correlationMatrix());
+        
+        return res;
+    }, 1);
+
+    min.setLowLimit(0., -0.1);
+    min.setHighLimit(0., 100000.);
+    init(0) = 0.1;
+    min.setInit(init);
+    prec = 0.1;
+    min.setPrecision(prec);
+    reg = min(cdr)(0);
+    for (unsigned int pass = 0; pass < nPass; pass++)
+    {
+        min.setLowLimit(0., (1.-10.*prec)*reg);
+        min.setHighLimit(0., (1.+10.*prec)*reg);
+        init(0) = reg;
+        min.setInit(init);
+        prec *= 0.1;
+        min.setPrecision(prec);
+        reg = min(cdr)(0);
+    }
+
+    return reg;
+}
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_DataFilter_hpp_

lib/LatAnalyze/Statistics/Histogram.cpp

@@ -146,16 +146,6 @@ double Histogram::getX(const Index i) const
     return x_(i);
 }
 
-double Histogram::getXMin(void) const
-{
-    return xMin_;
-}
-
-double Histogram::getXMax(void) const
-{
-    return xMax_;
-}
-
 double Histogram::operator[](const Index i) const
 {
     return bin_(i)*(isNormalized() ? norm_ : 1.);

lib/LatAnalyze/Statistics/Histogram.hpp

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

lib/LatAnalyze/Statistics/MatSample.hpp

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