Merge 269d0c338e into 4823426d55

.gitignore

@@ -25,9 +25,10 @@ lib/*Lexer.cpp
 lib/*Parser.cpp
 lib/*Parser.hpp
 
-# Eigen headers
+# Eigen headers and archives
-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.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

lib/Core/MathInterpreter.cpp

@@ -217,146 +217,25 @@ void RunContext::reset(void)
 #define CODE_WIDTH 6
 #define CODE_MOD   setw(CODE_WIDTH) << left
 
-// Instruction operator ////////////////////////////////////////////////////////
+auto readConstant(Program::const_iterator ip)
-ostream &Latan::operator<<(ostream& out, const Instruction& ins)
+    -> std::tuple<double, Program::const_iterator>
 {
-    ins.print(out);
+    double value = 0.0;
+    std::copy(ip, ip + sizeof(double), reinterpret_cast<std::uint8_t*>(&value));
 
-    return out;
+    return std::make_tuple(value, ip + sizeof(double));
 }
 
-// Push constructors ///////////////////////////////////////////////////////////
+auto readAddress(Program::const_iterator ip)
-Push::Push(const double val)
+    -> std::tuple<unsigned int, Program::const_iterator>
-: type_(ArgType::Constant)
-, val_(val)
-, address_(0)
-, name_("")
-{}
-
-Push::Push(const unsigned int address, const string &name)
-: type_(ArgType::Variable)
-, val_(0.0)
-, address_(address)
-, name_(name)
-{}
-
-// Push execution //////////////////////////////////////////////////////////////
-void Push::operator()(RunContext &context) const
 {
-    if (type_ == ArgType::Constant)
+    unsigned int address = 0.0;
-    {
+    const auto end = ip + sizeof(unsigned int);
-        context.stack().push(val_);
+    std::copy(ip, end, reinterpret_cast<std::uint8_t*>(&address));
-    }
+
-    else
+    return std::make_tuple(address, end);
-    {
-        context.stack().push(context.getVariable(address_));
-    }
-    context.incrementInsIndex();
 }
 
-// Push print //////////////////////////////////////////////////////////////////
-void Push::print(ostream &out) const
-{
-    out << CODE_MOD << "push";
-    if (type_ == ArgType::Constant)
-    {
-        out << CODE_MOD << val_;
-    }
-    else
-    {
-        out << CODE_MOD << name_ << " @v" << address_;
-    }
-}
-
-// Pop constructor /////////////////////////////////////////////////////////////
-Pop::Pop(const unsigned int address, const string &name)
-: address_(address)
-, name_(name)
-{}
-
-// Pop execution ///////////////////////////////////////////////////////////////
-void Pop::operator()(RunContext &context) const
-{
-    if (!name_.empty())
-    {
-        context.setVariable(address_, context.stack().top());
-    }
-    context.stack().pop();
-    context.incrementInsIndex();
-}
-
-// Pop print ///////////////////////////////////////////////////////////////////
-void Pop::print(ostream &out) const
-{
-    out << CODE_MOD << "pop" << CODE_MOD << name_ << " @v" << address_;
-}
-
-// Store constructor ///////////////////////////////////////////////////////////
-Store::Store(const unsigned int address, const string &name)
-: address_(address)
-, name_(name)
-{}
-
-// Store execution /////////////////////////////////////////////////////////////
-void Store::operator()(RunContext &context) const
-{
-    if (!name_.empty())
-    {
-        context.setVariable(address_, context.stack().top());
-    }
-    context.incrementInsIndex();
-}
-
-// Store print /////////////////////////////////////////////////////////////////
-void Store::print(ostream &out) const
-{
-    out << CODE_MOD << "store" << CODE_MOD << name_ << " @v" << address_;
-}
-
-// Call constructor ////////////////////////////////////////////////////////////
-Call::Call(const unsigned int address, const string &name)
-: address_(address)
-, name_(name)
-{}
-
-// Call execution //////////////////////////////////////////////////////////////
-void Call::operator()(RunContext &context) const
-{
-    context.stack().push((*context.getFunction(address_))(context.stack()));
-    context.incrementInsIndex();
-}
-
-// Call print //////////////////////////////////////////////////////////////////
-void Call::print(ostream &out) const
-{
-    out << CODE_MOD << "call" << CODE_MOD << name_ << " @f" << address_;
-}
-
-// Math operations /////////////////////////////////////////////////////////////
-#define DEF_OP(name, nArg, exp, insName)\
-void name::operator()(RunContext &context) const\
-{\
-    double x[nArg];\
-    for (int i = 0; i < nArg; ++i)\
-    {\
-        x[nArg-1-i] = context.stack().top();\
-        context.stack().pop();\
-    }\
-    context.stack().push(exp);\
-    context.incrementInsIndex();\
-}\
-void name::print(ostream &out) const\
-{\
-    out << CODE_MOD << insName;\
-}
-
-DEF_OP(Neg, 1, -x[0],          "neg")
-DEF_OP(Add, 2, x[0] + x[1],    "add")
-DEF_OP(Sub, 2, x[0] - x[1],    "sub")
-DEF_OP(Mul, 2, x[0]*x[1],      "mul")
-DEF_OP(Div, 2, x[0]/x[1],      "div")
-DEF_OP(Pow, 2, pow(x[0],x[1]), "pow")
-
 /******************************************************************************
  *                        ExprNode implementation                             *
  ******************************************************************************/
@@ -442,29 +321,35 @@ ostream &Latan::operator<<(ostream &out, const ExprNode &n)
     return out;
 }
 
-#define PUSH_INS(program, type, ...)\
+// Bytecode helper functions ///////////////////////////////////////////////////
-program.push_back(unique_ptr<type>(new type(__VA_ARGS__)))
+void pushInstruction(Program &program, const Instruction instruction) {
-#define GET_ADDRESS(address, table, name)\
+    program.push_back(static_cast<std::uint8_t>(instruction));
-try\
+}
-{\
+
-    address = (table).at(name);\
+void pushAddress(Program &program, const unsigned int address) {
-}\
+    const auto address_ptr = reinterpret_cast<const std::uint8_t*>(&address);
-catch (out_of_range)\
+    const auto size = sizeof(unsigned int);
-{\
+    program.insert(program.end(), address_ptr, address_ptr + size);
-    address         = (table).size();\
+}
-    (table)[(name)] = address;\
+
-}\
+void pushConstant(Program &program, const double value) {
+    const auto value_ptr = reinterpret_cast<const std::uint8_t*>(&value);
+    const auto size = sizeof(double);
+    program.insert(program.end(), value_ptr, value_ptr + size);
+}
 
 // VarNode compile /////////////////////////////////////////////////////////////
 void VarNode::compile(Program &program, RunContext &context) const
 {
-    PUSH_INS(program, Push, context.getVariableAddress(getName()), getName());
+    pushInstruction(program, Instruction::LOAD);
+    pushAddress(program, context.getVariableAddress(getName()));
 }
 
 // CstNode compile /////////////////////////////////////////////////////////////
 void CstNode::compile(Program &program, RunContext &context __dumb) const
 {
-    PUSH_INS(program, Push, strTo<double>(getName()));
+    pushInstruction(program, Instruction::CONST);
+    pushConstant(program, strTo<double>(getName()));
 }
 
 // SemicolonNode compile ///////////////////////////////////////////////////////
@@ -482,6 +367,10 @@ void SemicolonNode::compile(Program &program, RunContext &context) const
         {
             n[i].compile(program, context);
         }
+        // Where a variable has just been assigned, pop it off the stack.
+        if (isAssign) {
+            pushInstruction(program, Instruction::POP);
+        }
     }
 }
 
@@ -492,19 +381,10 @@ void AssignNode::compile(Program &program, RunContext &context) const
     
     if (isDerivedFrom<VarNode>(&n[0]))
     {
-        bool hasSemicolonParent = isDerivedFrom<SemicolonNode>(getParent());
-        unsigned int address;
-        
         n[1].compile(program, context);
-        address = context.addVariable(n[0].getName());
+        const unsigned int address = context.addVariable(n[0].getName());
-        if (hasSemicolonParent)
+        pushInstruction(program, Instruction::STORE);
-        {
+        pushAddress(program, address);
-            PUSH_INS(program, Pop, address, n[0].getName());
-        }
-        else
-        {
-            PUSH_INS(program, Store, address, n[0].getName());
-        }
     }
     else
     {
@@ -519,19 +399,37 @@ void AssignNode::compile(Program &program, RunContext &context) const
 
 void MathOpNode::compile(Program &program, RunContext &context) const
 {
+#define PUSH_BINARY_OP(op, instruction) \
+    case op: \
+        pushInstruction(program, Instruction::instruction); \
+        break;
+
     auto &n = *this;
     
     for (Index i = 0; i < n.getNArg(); ++i)
     {
         n[i].compile(program, context);
     }
-    IFNODE("-", 1)   PUSH_INS(program, Neg,);
+    if (n.getName() == "-" and n.getNArg() == 1) {
-    ELIFNODE("+", 2) PUSH_INS(program, Add,);
+        pushInstruction(program, Instruction::NEG);
-    ELIFNODE("-", 2) PUSH_INS(program, Sub,);
+        return;
-    ELIFNODE("*", 2) PUSH_INS(program, Mul,);
+    }
-    ELIFNODE("/", 2) PUSH_INS(program, Div,);
+
-    ELIFNODE("^", 2) PUSH_INS(program, Pow,);
+    if (getNArg() != 2) {
-    ELSE LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+        LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+    }
+
+    switch (getName()[0]) {
+    PUSH_BINARY_OP('+', ADD)
+    PUSH_BINARY_OP('-', SUB)
+    PUSH_BINARY_OP('*', MUL)
+    PUSH_BINARY_OP('/', DIV)
+    PUSH_BINARY_OP('^', POW)
+    default:
+        LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
+    }
+
+#undef PUSH_BINARY_OP
 }
 
 // FuncNode compile ////////////////////////////////////////////////////////////
@@ -543,7 +441,8 @@ void FuncNode::compile(Program &program, RunContext &context) const
     {
         n[i].compile(program, context);
     }
-    PUSH_INS(program, Call, context.getFunctionAddress(getName()), getName());
+    pushInstruction(program, Instruction::CALL);
+    pushAddress(program, context.getFunctionAddress(getName()));
 }
 
 // ReturnNode compile ////////////////////////////////////////////////////////////
@@ -552,7 +451,7 @@ void ReturnNode::compile(Program &program, RunContext &context) const
     auto &n = *this;
     
     n[0].compile(program, context);
-    program.push_back(nullptr);
+    pushInstruction(program, Instruction::RET);
 }
 
 /******************************************************************************
@@ -582,7 +481,7 @@ MathInterpreter::MathInterpreter(const std::string &code)
 // access //////////////////////////////////////////////////////////////////////
 const Instruction * MathInterpreter::operator[](const Index i) const
 {
-    return program_[i].get();
+    return reinterpret_cast<const Instruction*>(&program_[i]);
 }
 
 const ExprNode * MathInterpreter::getAST(void) const
@@ -592,7 +491,7 @@ const ExprNode * MathInterpreter::getAST(void) const
 
 void MathInterpreter::push(const Instruction *i)
 {
-    program_.push_back(unique_ptr<const Instruction>(i));
+    pushInstruction(program_, *i);
 }
 
 // initialization //////////////////////////////////////////////////////////////
@@ -695,7 +594,7 @@ void MathInterpreter::compile(RunContext &context)
             root_->compile(program_, context);
             for (unsigned int i = 0; i < program_.size(); ++i)
             {
-                if (!program_[i])
+                if (static_cast<Instruction>(program_[i]) == Instruction::RET)
                 {
                     gotReturn = true;
                     program_.resize(i);
@@ -726,20 +625,145 @@ void MathInterpreter::operator()(RunContext &context)
 
 void MathInterpreter::execute(RunContext &context) const
 {
-    context.setInsIndex(0);
+#define BINARY_OP_CASE(instruction, expr) \
-    while (context.getInsIndex() != program_.size())
+        case Instruction::instruction: { \
-    {
+            const auto second = context.stack().top(); \
-        (*(program_[context.getInsIndex()]))(context);
+            context.stack().pop(); \
+            const auto first = context.stack().top(); \
+            context.stack().pop(); \
+            context.stack().push(expr); \
+            break; \
+        }
+    auto ip = program_.begin();
+
+    while (ip != program_.end()) {
+        const auto instruction = static_cast<Instruction>(*ip);
+        ip++;
+
+        switch (instruction) {
+        BINARY_OP_CASE(ADD, first + second)
+        BINARY_OP_CASE(SUB, first - second)
+        BINARY_OP_CASE(MUL, first * second)
+        BINARY_OP_CASE(DIV, first / second)
+        BINARY_OP_CASE(POW, std::pow(first, second))
+        case Instruction::NEG: {
+            const auto operand = context.stack().top();
+            context.stack().pop();
+            context.stack().push(-operand);
+            break;
+        }
+        case Instruction::CONST: {
+            double value = 0.0;
+            std::tie(value, ip) = readConstant(ip);
+            context.stack().push(value);
+            break;
+        }
+        case Instruction::POP:
+            context.stack().pop();
+            break;
+        case Instruction::LOAD: {
+            unsigned int address = 0;
+            std::tie(address, ip) = readAddress(ip);
+            context.stack().push(context.getVariable(address));
+            break;
+        }
+        case Instruction::STORE: {
+            unsigned int address = 0;
+            std::tie(address, ip) = readAddress(ip);
+            context.setVariable(address, context.stack().top());
+            break;
+        }
+        case Instruction::CALL: {
+            unsigned int address = 0;
+            std::tie(address, ip) = readAddress(ip);
+            auto& stack = context.stack();
+            stack.push((*context.getFunction(address))(stack));
+            break;
+        }
+        case Instruction::RET:
+            break;
+        }
     }
+
+#undef BINARY_OP_CASE
+}
+
+Program::const_iterator instructionToStream(
+    ostream &out, Program::const_iterator ip)
+{
+    const auto instruction = static_cast<Instruction>(*ip);
+    ip++;
+
+    switch (instruction) {
+    case Instruction::ADD:
+        out << "ADD";
+        break;
+    case Instruction::SUB:
+        out << "SUB";
+        break;
+    case Instruction::MUL:
+        out << "MUL";
+        break;
+    case Instruction::DIV:
+        out << "DIV";
+        break;
+    case Instruction::POW:
+        out << "POW";
+        break;
+    case Instruction::NEG:
+        out << "NEG";
+        break;
+    case Instruction::CONST: {
+        double value = 0.0;
+        std::tie(value, ip) = readConstant(ip);
+        out << CODE_MOD << setfill(' ') << "CONST" << value;
+        break;
+    }
+    case Instruction::POP:
+        out << "POP";
+        break;
+    case Instruction::LOAD: {
+        unsigned int address = 0;
+        std::tie(address, ip) = readAddress(ip);
+        out << CODE_MOD << setfill(' ') << "LOAD" << address;
+        break;
+    }
+    case Instruction::STORE: {
+        unsigned int address = 0;
+        std::tie(address, ip) = readAddress(ip);
+        out << CODE_MOD << setfill(' ') << "STORE" << address;
+        break;
+    }
+    case Instruction::CALL: {
+        unsigned int address = 0;
+        std::tie(address, ip) = readAddress(ip);
+        out << CODE_MOD << setfill(' ') << "CALL" << address;
+        break;
+    }
+    case Instruction::RET:
+        out << "RET";
+        break;
+    }
+
+    return ip;
+}
+
+ostream &programToStream(ostream &out, const Program &program)
+{
+    auto ip = program.begin();
+
+    while (ip != program.end()) {
+        const auto i = std::distance(program.begin(), ip);
+        cout << setw(4) << setfill('0') << right << i << "  ";
+        ip = instructionToStream(out, ip);
+        out << '\n';
+    }
+
+    return out;
 }
 
 // IO //////////////////////////////////////////////////////////////////////////
 ostream &Latan::operator<<(ostream &out, const MathInterpreter &program)
 {
-    for (unsigned int i = 0; i < program.program_.size(); ++i)
+    return programToStream(out, program.program_);
-    {
-        out << *(program.program_[i]) << endl;
-    }
-    
-    return out;
 }

lib/Core/MathInterpreter.hpp

@@ -20,6 +20,8 @@
 #ifndef Latan_MathInterpreter_hpp_
 #define	Latan_MathInterpreter_hpp_
 
+#include <cstdint>
+
 #include <LatAnalyze/Functional/Function.hpp>
 #include <LatAnalyze/Global.hpp>
 #include <LatAnalyze/Core/ParserState.hpp>
@@ -79,108 +81,26 @@ private:
  *                         Instruction classes                                *
  ******************************************************************************/
 // Abstract base
-class Instruction
+
-{
+enum class Instruction : std::uint8_t {
-public:
+    ADD,
-    // destructor
+    SUB,
-    virtual ~Instruction(void) = default;
+    MUL,
-    // instruction execution
+    DIV,
-    virtual void operator()(RunContext &context) const = 0;
+    POW,
-    friend std::ostream & operator<<(std::ostream &out, const Instruction &ins);
+    NEG,
-private:
+    CONST,
-    virtual void print(std::ostream &out) const = 0;
+    POP,
+    LOAD,
+    STORE,
+    CALL,
+    RET,
 };
 
 std::ostream & operator<<(std::ostream &out, const Instruction &ins);
 
 // Instruction container
-typedef std::vector<std::unique_ptr<const Instruction>> Program;
+typedef std::vector<std::uint8_t> Program;
-
-// Push
-class Push: public Instruction
-{
-private:
-    enum class ArgType
-    {
-        Constant = 0,
-        Variable = 1
-    };
-public:
-    //constructors
-    explicit Push(const double val);
-    explicit Push(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    ArgType      type_;
-    double       val_;
-    unsigned int address_;
-    std::string  name_;
-};
-
-// Pop
-class Pop: public Instruction
-{
-public:
-    //constructor
-    explicit Pop(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    unsigned int address_;
-    std::string name_;
-};
-
-// Store
-class Store: public Instruction
-{
-public:
-    //constructor
-    explicit Store(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    unsigned int address_;
-    std::string name_;
-};
-
-// Call function
-class Call: public Instruction
-{
-public:
-    //constructor
-    explicit Call(const unsigned int address, const std::string &name);
-    // instruction execution
-    virtual void operator()(RunContext &context) const;
-private:
-    virtual void print(std::ostream& out) const;
-private:
-    unsigned int address_;
-    std::string name_;
-};
-
-// Floating point operations
-#define DECL_OP(name)\
-class name: public Instruction\
-{\
-public:\
-    virtual void operator()(RunContext &context) const;\
-private:\
-    virtual void print(std::ostream &out) const;\
-}
-
-DECL_OP(Neg);
-DECL_OP(Add);
-DECL_OP(Sub);
-DECL_OP(Mul);
-DECL_OP(Div);
-DECL_OP(Pow);
 
 /******************************************************************************
  *                            Expression node classes                         *

lib/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, 
+template <typename MatType>
-                            const std::vector<double> &filter, const Index offset)
+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())
+    ::upsample(out, in);
-    {
+}
-        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.);
+{
-    for (Index i = 0; i < in.size(); i ++)
+    ::upsample(out, in);
-    {
-        out(2*i) = in(i);
-    }
 }
 
 // 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,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/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/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_

lib/Statistics/StatArray.hpp

@@ -52,10 +52,13 @@ 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:
@@ -192,6 +195,79 @@ 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
 {