Merge 269d0c338e into 6739019c83

.gitignore

@@ -25,10 +25,9 @@ lib/*Lexer.cpp
 lib/*Parser.cpp
 lib/*Parser.hpp
 
-# Eigen headers and archives
+# Eigen headers
-lib/Eigen
+lib/Eigen/*
 lib/eigen_files.mk
-eigen-*.tar.bz2
 
 # CI builds
 ci-scripts/local/*

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

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 DWT::filterConvolution(DVec &out, const DVec &data, 
-void filterConvolution(MatType &out, const MatType &data, 
+                            const std::vector<double> &filter, const Index offset)
-                       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);
+    if (out.size() < in.size()/2)
-}
-
-void DWT::downsample(DMat &out, const DMat &in)
-{
-    ::downsample(out, in);
-}
-
-template <typename MatType>
-void upsample(MatType &out, const MatType &in)
-{
-    if (out.size() < 2*in.size())
     {
-        LATAN_ERROR(Size, "output rows smaller than twice the input rows");
+        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(i/2) = in(i);
-    }
-    out.block(0, 0, 2*in.size(), out.cols()).fill(0.);
-    for (Index j = 0; j < in.cols(); j++)
-    for (Index i = 0; i < in.size(); i ++)
-    {
-        out(2*i, j) = in(i, j);
     }
 }
 
 void DWT::upsample(DVec &out, const DVec &in)
 {
-    ::upsample(out, in);
+    if (out.size() < 2*in.size())
-}
+    {
-
+        LATAN_ERROR(Size, "output vector smaller than twice the input vector size");
-void DWT::upsample(DMat &out, const DMat &in)
+    }
-{
+    out.segment(0, 2*in.size()).fill(0.);
-    ::upsample(out, in);
+    for (Index i = 0; i < in.size(); i ++)
+    {
+        out(2*i) = in(i);
+    }
 }
 
 // DWT /////////////////////////////////////////////////////////////////////////

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,13 +40,9 @@ 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;

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

lib/Statistics/StatArray.hpp

@@ -52,13 +52,10 @@ public:
     // statistics
     void bin(Index binSize);
     T    sum(const Index pos = 0, const Index n = -1) const;
+    T    meanOld(const Index pos = 0, const Index n = -1) const;
     T    mean(const Index pos = 0, const Index n = -1) const;
     T    covariance(const StatArray<T, os> &array) const;
     T    variance(void) const;
-    T    covarianceMatrix(const StatArray<T, os> &data) const;
-    T    varianceMatrix(void) const;
-    T    correlationMatrix(void) const;
-
     // IO type
     virtual IoType getType(void) const;
 public:
@@ -195,79 +192,6 @@ T StatArray<T, os>::variance(void) const
     return covariance(*this);
 }
 
-template <typename MatType, Index os>
-MatType StatArray<MatType, os>::covarianceMatrix(
-    const StatArray<MatType, os> &data) const
-{
-    if (((*this)[central].cols() != 1) or (data[central].cols() != 1))
-    {
-        LATAN_ERROR(Size, "samples have more than one column");
-    }
-
-    Index   n1 = (*this)[central].rows(), n2 = data[central].rows();
-    Index   nSample = this->size();
-    MatType tmp1(n1, nSample), tmp2(n2, nSample), res(n1, n2);
-    MatType s1(n1, 1), s2(n2, 1), one(nSample, 1);
-
-    one.fill(1.);
-    s1.fill(0.);
-    s2.fill(0.);
-    for (unsigned int s = 0; s < nSample; ++s)
-    {
-        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          += data[s];
-        tmp2.col(s)  = data[s];
-    }
-    tmp2 -= s2*one.transpose()/static_cast<double>(nSample);
-    res   = tmp1*tmp2.transpose()/static_cast<double>(nSample - 1);
-
-    return res;
-}
-
-template <typename MatType, Index os>
-MatType StatArray<MatType, os>::varianceMatrix(void) const
-{
-    if ((*this)[0].cols() != 1)
-    {
-        LATAN_ERROR(Size, "samples have more than one column");
-    }
-
-    Index   n1 = (*this)[0].rows();
-    Index   nSample = this->size();
-    MatType tmp1(n1, nSample), res(n1, n1);
-    MatType s1(n1, 1), one(nSample, 1);
-
-    one.fill(1.);
-    s1.fill(0.);
-    for (unsigned int s = 0; s < nSample; ++s)
-    {
-        s1          += (*this)[s];
-        tmp1.col(s)  = (*this)[s];
-    }
-    tmp1 -= s1*one.transpose()/static_cast<double>(nSample);
-    res   = tmp1*tmp1.transpose()/static_cast<double>(nSample - 1);
-
-    return res;
-}
-
-template <typename MatType, Index os>
-MatType StatArray<MatType, os>::correlationMatrix(void) const
-{
-    MatType res = varianceMatrix();
-    MatType invDiag(res.rows(), 1);
-
-    invDiag = res.diagonal();
-    invDiag = invDiag.cwiseInverse().cwiseSqrt();
-    res     = (invDiag*invDiag.transpose()).cwiseProduct(res);
-
-    return res;
-}
-
 // reduction operations ////////////////////////////////////////////////////////
 namespace StatOp
 {