general interface for filtering data

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

lib/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,9 +41,13 @@ 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

@@ -0,0 +1,87 @@
+/*
+ * 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::DataFilter(const vector<double> &filter, const bool downsample)
+: filter_(filter), downsample_(downsample)
+{}
+
+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_);
+}
+
+void DataFilter::operator()(DMatSample &out, const DMatSample &in)
+{
+    FOR_STAT_ARRAY(in, s)
+    {
+        (*this)(out[s], in[s]);
+    }
+}
+
+LaplaceDataFilter::LaplaceDataFilter(const bool downsample)
+: DataFilter({1., -2. , 1.}, downsample)
+{}
+
+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);
+}
+
+void LaplaceDataFilter::operator()(DMatSample &out, const DMatSample &in, const double lambda)
+{
+    filter_[1] = -2. - lambda;
+    DataFilter::operator()(out, in);
+}

lib/Physics/DataFilter.hpp

@@ -0,0 +1,58 @@
+/*
+ * 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/Statistics/MatSample.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+class DataFilter
+{
+public:
+    // constructor
+    DataFilter(const std::vector<double> &filter, const bool downsample = false);
+    // filtering
+    virtual void operator()(DVec &out, const DVec &in);
+    virtual void operator()(DMat &out, const DMat &in);
+    virtual void operator()(DMatSample &out, const DMatSample &in);
+protected:
+    std::vector<double> filter_;
+private:
+    bool                downsample_;
+    DVec                vBuf_;
+    DMat                mBuf_;
+};
+
+class LaplaceDataFilter: public DataFilter
+{
+public:
+    // constructor
+    LaplaceDataFilter(const bool downsample = false);
+    // filtering
+    virtual void operator()(DVec &out, const DVec &in, const double lambda = 0.);
+    virtual void operator()(DMat &out, const DMat &in, const double lambda = 0.);
+    virtual void operator()(DMatSample &out, const DMatSample &in, const double lambda = 0.);
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_DataFilter_hpp_