Merge pull request #19 from AndrewYongZhenNing/develop

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

.github/workflows/build-macos.yml

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

.gitignore

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

CMakeLists.txt

@@ -1,64 +0,0 @@
-# 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.4.0.tar.bz2
+./update_eigen.sh eigen-3.3.8.tar.bz2
 autoreconf -fvi

build.sh

@@ -0,0 +1,24 @@
+#!/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  n             = 8;
+constexpr Index  size          = 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(n, n);
-    DVec       mean(n);
-    DMatSample sample(nSample, n, 1);
+    DMat       var(size, size);
+    DVec       mean(size);
+    DMatSample sample(nSample, size, 1);
 
     cout << "-- generating " << nSample << " Gaussian random vectors..." << endl;
-    var   = DMat::Random(n, n);
+    var   = DMat::Random(size, size);
     var  *= var.adjoint();
-    mean  = DVec::Random(n); 
+    mean  = DVec::Random(size); 
     RandomNormal mgauss(mean, var, rd());
     sample[central] = mgauss();
     FOR_STAT_ARRAY(sample, s)

lib/CMakeLists.txt

@@ -1,81 +0,0 @@
-# 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

@@ -1,28 +0,0 @@
-@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/Core/Math.cpp

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

lib/Core/Math.hpp

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

lib/Core/Plot.cpp

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

lib/Core/Plot.hpp

@@ -227,7 +227,6 @@ 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;
@@ -315,14 +314,13 @@ class LogScale: public PlotModifier
 {
 public:
     // constructor
-    explicit LogScale(const Axis axis, const double basis = 10);
+    explicit LogScale(const Axis axis);
     // destructor
     virtual ~LogScale(void) = default;
     // modifier
     virtual void operator()(PlotOptions &option) const;
 private:
     const Axis axis_;
-    const double basis_;
 };
 
 class PlotRange: public PlotModifier

lib/Core/Utilities.hpp

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

lib/Global.cpp

@@ -25,10 +25,11 @@ using namespace Latan;
 
 PlaceHolder Latan::_;
 
-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 + "] ";
+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) + "] ";
 
 void Env::function(void)
 {}

lib/LatAnalyze

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

lib/LatAnalyze/Physics/DataFilter.cpp

@@ -1,83 +0,0 @@
-/*
- * 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

@@ -1,139 +0,0 @@
-/*
- * 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/Makefile.am

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

lib/Numerical/DWT.cpp

@@ -32,91 +32,46 @@ DWT::DWT(const DWTFilter &filter)
 {}
 
 // convolution primitive ///////////////////////////////////////////////////////
-template <typename MatType>
-void filterConvolution(MatType &out, const MatType &data, 
-                       const std::vector<double> &filter, const Index offset)
+void DWT::filterConvolution(DVec &out, const DVec &data, 
+                            const std::vector<double> &filter, const Index offset)
 {
-    Index n = data.rows(), nf = n*filter.size();
+    Index n = data.size(), nf = n*filter.size();
 
-    out.resizeLike(data);
+    out.resize(n);
     out.fill(0.);
     for (unsigned int i = 0; i < filter.size(); ++i)
     {
-        FOR_MAT(out, j, k)
+        FOR_VEC(out, j)
         {
-            out(j, k) += filter[i]*data((j + i + nf - offset) % n, k);
+            out(j) += filter[i]*data((j + i + nf - offset) % n);
         }
     }
 }
 
-void DWT::filterConvolution(DVec &out, const DVec &data, 
-                            const std::vector<double> &filter, const Index offset)
-{
-    ::filterConvolution(out, data, filter, offset);
-}
-
-void DWT::filterConvolution(DMat &out, const DMat &data, 
-                            const std::vector<double> &filter, const Index offset)
-{
-    ::filterConvolution(out, data, filter, offset);
-}
-
 // downsampling/upsampling primitives //////////////////////////////////////////
-template <typename MatType>
-void downsample(MatType &out, const MatType &in)
-{
-    if (out.rows() < in.rows()/2)
-    {
-        LATAN_ERROR(Size, "output rows smaller than half the input vector rows");
-    }
-    if (out.cols() != in.cols())
-    {
-        LATAN_ERROR(Size, "output and input number of columns mismatch");
-    }
-    for (Index j = 0; j < in.cols(); j++)
-    for (Index i = 0; i < in.rows(); i += 2)
-    {
-        out(i/2, j) = in(i, j);
-    }
-}
-
 void DWT::downsample(DVec &out, const DVec &in)
 {
-    ::downsample(out, in);
-}
-
-void DWT::downsample(DMat &out, const DMat &in)
-{
-    ::downsample(out, in);
-}
-
-template <typename MatType>
-void upsample(MatType &out, const MatType &in)
-{
-    if (out.size() < 2*in.size())
+    if (out.size() < in.size()/2)
     {
-        LATAN_ERROR(Size, "output rows smaller than twice the input rows");
+        LATAN_ERROR(Size, "output vector smaller than half the input vector size");
     }
-    if (out.cols() != in.cols())
+    for (Index i = 0; i < in.size(); i += 2)
     {
-        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);
+        out(i/2) = in(i);
     }
 }
 
 void DWT::upsample(DVec &out, const DVec &in)
 {
-    ::upsample(out, in);
-}
-
-void DWT::upsample(DMat &out, const DMat &in)
-{
-    ::upsample(out, 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 /////////////////////////////////////////////////////////////////////////
@@ -180,26 +135,3 @@ 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/Numerical/DWT.hpp

@@ -22,7 +22,6 @@
 
 #include <LatAnalyze/Global.hpp>
 #include <LatAnalyze/Numerical/DWTFilters.hpp>
-#include <LatAnalyze/Core/Mat.hpp>
 
 BEGIN_LATAN_NAMESPACE
 
@@ -41,18 +40,12 @@ 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/Numerical/MinuitMinimizer.cpp

@@ -30,12 +30,8 @@ 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/Physics/CorrelatorFitter.cpp

@@ -253,16 +253,39 @@ DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts)
     }
 }
 
-DMatSample CorrelatorUtils::fold(const DMatSample &c)
+DMatSample CorrelatorUtils::fold(const DMatSample &c, const CorrelatorModels::ModelPar &par)
 {
     const Index nt  = c[central].rows();
     DMatSample  buf = c;
-
-    FOR_STAT_ARRAY(buf, s)
+    int         sign;
+    bool        fold = false;
+    
+    switch (par.type)
     {
-        for (Index t = 0; t < nt; ++t)
+        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)
         {
-            buf[s](t) = 0.5*(c[s](t) + c[s]((nt - t) % nt));
+            for (Index t = 0; t < nt; ++t)
+            {
+                buf[s](t) = 0.5*(c[s](t) + sign*c[s]((nt - t) % nt));
+            }
         }
     }
 

lib/Physics/CorrelatorFitter.hpp

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

lib/Statistics/Histogram.cpp

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

lib/Statistics/Histogram.hpp

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

lib/Statistics/MatSample.hpp

@@ -103,6 +103,10 @@ public:
                      const Index nCol);
     // resize all matrices
     void resizeMat(const Index nRow, const Index nCol);
+    // covariance matrix
+    Mat<T> covarianceMatrix(const MatSample<T> &sample) const;
+    Mat<T> varianceMatrix(void) const;
+    Mat<T> correlationMatrix(void) const;
 };
 
 // non-member operators
@@ -379,6 +383,79 @@ 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_