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:

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/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/MathInterpreter.cpp

@@ -217,25 +217,146 @@ void RunContext::reset(void)
 #define CODE_WIDTH 6
 #define CODE_MOD   setw(CODE_WIDTH) << left
 
-auto readConstant(Program::const_iterator ip)
-    -> std::tuple<double, Program::const_iterator>
+// Instruction operator ////////////////////////////////////////////////////////
+ostream &Latan::operator<<(ostream& out, const Instruction& ins)
 {
-    double value = 0.0;
-    std::copy(ip, ip + sizeof(double), reinterpret_cast<std::uint8_t*>(&value));
-
-    return std::make_tuple(value, ip + sizeof(double));
+    ins.print(out);
+    
+    return out;
 }
 
-auto readAddress(Program::const_iterator ip)
-    -> std::tuple<unsigned int, Program::const_iterator>
-{
-    unsigned int address = 0.0;
-    const auto end = ip + sizeof(unsigned int);
-    std::copy(ip, end, reinterpret_cast<std::uint8_t*>(&address));
+// Push constructors ///////////////////////////////////////////////////////////
+Push::Push(const double val)
+: type_(ArgType::Constant)
+, val_(val)
+, address_(0)
+, name_("")
+{}
 
-    return std::make_tuple(address, end);
+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)
+    {
+        context.stack().push(val_);
+    }
+    else
+    {
+        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                             *
  ******************************************************************************/
@@ -321,35 +442,29 @@ ostream &Latan::operator<<(ostream &out, const ExprNode &n)
     return out;
 }
 
-// Bytecode helper functions ///////////////////////////////////////////////////
-void pushInstruction(Program &program, const Instruction instruction) {
-    program.push_back(static_cast<std::uint8_t>(instruction));
-}
-
-void pushAddress(Program &program, const unsigned int address) {
-    const auto address_ptr = reinterpret_cast<const std::uint8_t*>(&address);
-    const auto size = sizeof(unsigned int);
-    program.insert(program.end(), address_ptr, address_ptr + size);
-}
-
-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);
-}
+#define PUSH_INS(program, type, ...)\
+program.push_back(unique_ptr<type>(new type(__VA_ARGS__)))
+#define GET_ADDRESS(address, table, name)\
+try\
+{\
+    address = (table).at(name);\
+}\
+catch (out_of_range)\
+{\
+    address         = (table).size();\
+    (table)[(name)] = address;\
+}\
 
 // VarNode compile /////////////////////////////////////////////////////////////
 void VarNode::compile(Program &program, RunContext &context) const
 {
-    pushInstruction(program, Instruction::LOAD);
-    pushAddress(program, context.getVariableAddress(getName()));
+    PUSH_INS(program, Push, context.getVariableAddress(getName()), getName());
 }
 
 // CstNode compile /////////////////////////////////////////////////////////////
 void CstNode::compile(Program &program, RunContext &context __dumb) const
 {
-    pushInstruction(program, Instruction::CONST);
-    pushConstant(program, strTo<double>(getName()));
+    PUSH_INS(program, Push, strTo<double>(getName()));
 }
 
 // SemicolonNode compile ///////////////////////////////////////////////////////
@@ -367,10 +482,6 @@ 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);
-        }
     }
 }
 
@@ -381,10 +492,19 @@ 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);
-        const unsigned int address = context.addVariable(n[0].getName());
-        pushInstruction(program, Instruction::STORE);
-        pushAddress(program, address);
+        address = context.addVariable(n[0].getName());
+        if (hasSemicolonParent)
+        {
+            PUSH_INS(program, Pop, address, n[0].getName());
+        }
+        else
+        {
+            PUSH_INS(program, Store, address, n[0].getName());
+        }
     }
     else
     {
@@ -399,37 +519,19 @@ 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);
     }
-    if (n.getName() == "-" and n.getNArg() == 1) {
-        pushInstruction(program, Instruction::NEG);
-        return;
-    }
-
-    if (getNArg() != 2) {
-        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
+    IFNODE("-", 1)   PUSH_INS(program, Neg,);
+    ELIFNODE("+", 2) PUSH_INS(program, Add,);
+    ELIFNODE("-", 2) PUSH_INS(program, Sub,);
+    ELIFNODE("*", 2) PUSH_INS(program, Mul,);
+    ELIFNODE("/", 2) PUSH_INS(program, Div,);
+    ELIFNODE("^", 2) PUSH_INS(program, Pow,);
+    ELSE LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
 }
 
 // FuncNode compile ////////////////////////////////////////////////////////////
@@ -441,8 +543,7 @@ void FuncNode::compile(Program &program, RunContext &context) const
     {
         n[i].compile(program, context);
     }
-    pushInstruction(program, Instruction::CALL);
-    pushAddress(program, context.getFunctionAddress(getName()));
+    PUSH_INS(program, Call, context.getFunctionAddress(getName()), getName());
 }
 
 // ReturnNode compile ////////////////////////////////////////////////////////////
@@ -451,7 +552,7 @@ void ReturnNode::compile(Program &program, RunContext &context) const
     auto &n = *this;
     
     n[0].compile(program, context);
-    pushInstruction(program, Instruction::RET);
+    program.push_back(nullptr);
 }
 
 /******************************************************************************
@@ -481,7 +582,7 @@ MathInterpreter::MathInterpreter(const std::string &code)
 // access //////////////////////////////////////////////////////////////////////
 const Instruction * MathInterpreter::operator[](const Index i) const
 {
-    return reinterpret_cast<const Instruction*>(&program_[i]);
+    return program_[i].get();
 }
 
 const ExprNode * MathInterpreter::getAST(void) const
@@ -491,7 +592,7 @@ const ExprNode * MathInterpreter::getAST(void) const
 
 void MathInterpreter::push(const Instruction *i)
 {
-    pushInstruction(program_, *i);
+    program_.push_back(unique_ptr<const Instruction>(i));
 }
 
 // initialization //////////////////////////////////////////////////////////////
@@ -594,7 +695,7 @@ void MathInterpreter::compile(RunContext &context)
             root_->compile(program_, context);
             for (unsigned int i = 0; i < program_.size(); ++i)
             {
-                if (static_cast<Instruction>(program_[i]) == Instruction::RET)
+                if (!program_[i])
                 {
                     gotReturn = true;
                     program_.resize(i);
@@ -625,145 +726,20 @@ void MathInterpreter::operator()(RunContext &context)
 
 void MathInterpreter::execute(RunContext &context) const
 {
-#define BINARY_OP_CASE(instruction, expr) \
-        case Instruction::instruction: { \
-            const auto second = context.stack().top(); \
-            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;
-        }
+    context.setInsIndex(0);
+    while (context.getInsIndex() != program_.size())
+    {
+        (*(program_[context.getInsIndex()]))(context);
     }
-
-#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)
 {
-    return programToStream(out, program.program_);
+    for (unsigned int i = 0; i < program.program_.size(); ++i)
+    {
+        out << *(program.program_[i]) << endl;
+    }
+    
+    return out;
 }

lib/Core/MathInterpreter.hpp

@@ -20,8 +20,6 @@
 #ifndef Latan_MathInterpreter_hpp_
 #define	Latan_MathInterpreter_hpp_
 
-#include <cstdint>
-
 #include <LatAnalyze/Functional/Function.hpp>
 #include <LatAnalyze/Global.hpp>
 #include <LatAnalyze/Core/ParserState.hpp>
@@ -81,26 +79,108 @@ private:
  *                         Instruction classes                                *
  ******************************************************************************/
 // Abstract base
-
-enum class Instruction : std::uint8_t {
-    ADD,
-    SUB,
-    MUL,
-    DIV,
-    POW,
-    NEG,
-    CONST,
-    POP,
-    LOAD,
-    STORE,
-    CALL,
-    RET,
+class Instruction
+{
+public:
+    // destructor
+    virtual ~Instruction(void) = default;
+    // instruction execution
+    virtual void operator()(RunContext &context) const = 0;
+    friend std::ostream & operator<<(std::ostream &out, const Instruction &ins);
+private:
+    virtual void print(std::ostream &out) const = 0;
 };
 
 std::ostream & operator<<(std::ostream &out, const Instruction &ins);
 
 // Instruction container
-typedef std::vector<std::uint8_t> Program;
+typedef std::vector<std::unique_ptr<const Instruction>> 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/Core/Plot.cpp

@@ -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 //////////////////////////////////////////////////////
@@ -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/Numerical/DWT.cpp

@@ -135,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

@@ -46,8 +46,6 @@ public:
     // 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/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/XYSampleData.cpp

@@ -300,6 +300,67 @@ const XYStatData & XYSampleData::getData(void)
 }
 
 // fit /////////////////////////////////////////////////////////////////////////
+
+void XYSampleData::fitSample(std::vector<Minimizer *> &minimizer,
+                             const std::vector<const DoubleModel *> &v,
+                             SampleFitResult &result,
+                             DVec &init,
+                             Index s)
+{
+    result.resize(nSample_);
+    result.chi2_.resize(nSample_);
+    result.model_.resize(v.size());
+    FitResult sampleResult;
+    setDataToSample(s);
+    if (s == central)
+    {
+        sampleResult        = data_.fit(minimizer, init, v);
+        init                = sampleResult.segment(0, init.size());
+        result.nPar_        = sampleResult.getNPar();
+        result.nDof_        = sampleResult.nDof_;
+        result.parName_     = sampleResult.parName_;
+        result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
+    }
+    else
+    {
+        sampleResult = data_.fit(*(minimizer.back()), init, v);
+    }
+    result[s]       = sampleResult;
+    result.chi2_[s] = sampleResult.getChi2();
+    for (unsigned int j = 0; j < v.size(); ++j)
+    {
+        result.model_[j].resize(nSample_);
+        result.model_[j][s] = sampleResult.getModel(j);
+    }
+
+    
+}
+
+SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
+                                  const DVec &init,
+                                  const std::vector<const DoubleModel *> &v,
+                                  Index s)
+{
+    computeVarMat();
+    
+    SampleFitResult result;
+    DVec      initCopy = init;
+
+    fitSample(minimizer, v, result, initCopy, s);
+    
+    return result;
+}
+
+SampleFitResult XYSampleData::fit(Minimizer &minimizer,
+                                  const DVec &init,
+                                  const std::vector<const DoubleModel *> &v,
+                                  Index s)
+{
+    vector<Minimizer *> mv{&minimizer};
+    
+    return fit(mv, init, v, s);
+}
+
 SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
                                   const DVec &init,
                                   const std::vector<const DoubleModel *> &v)
@@ -307,43 +368,14 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
     computeVarMat();
     
     SampleFitResult      result;
-    FitResult            sampleResult;
     DVec                 initCopy = init;
     Minimizer::Verbosity verbCopy = minimizer.back()->getVerbosity();
     
-    result.resize(nSample_);
-    result.chi2_.resize(nSample_);
-    result.model_.resize(v.size());
     FOR_STAT_ARRAY(result, s)
     {
-        setDataToSample(s);
-        if (s == central)
-        {
-            sampleResult = data_.fit(minimizer, initCopy, v);
-            initCopy     = sampleResult.segment(0, initCopy.size());
-            if (verbCopy != Minimizer::Verbosity::Debug)
-            {
-                minimizer.back()->setVerbosity(Minimizer::Verbosity::Silent);
-            }
-        }
-        else
-        {
-            
-            sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
-        }
-        result[s]       = sampleResult;
-        result.chi2_[s] = sampleResult.getChi2();
-        for (unsigned int j = 0; j < v.size(); ++j)
-        {
-            result.model_[j].resize(nSample_);
-            result.model_[j][s] = sampleResult.getModel(j);
-        }
+        fitSample(minimizer, v, result, initCopy, s);
     }
     minimizer.back()->setVerbosity(verbCopy);
-    result.nPar_       = sampleResult.getNPar();
-    result.nDof_       = sampleResult.nDof_;
-    result.parName_    = sampleResult.parName_;
-    result.corrRangeDb_ = Math::cdr(getFitCorrMat());
     
     return result;
 }

lib/Statistics/XYSampleData.hpp

@@ -103,9 +103,16 @@ public:
     // get internal XYStatData
     const XYStatData & getData(void);
     // fit
-    SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
+    void            fitSample(std::vector<Minimizer *> &minimizer, 
+                              const std::vector<const DoubleModel *> &v, 
+                              SampleFitResult &sampleResult, DVec &init, Index s);
+    SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init, 
+                        const std::vector<const DoubleModel *> &v, Index s);
+    SampleFitResult fit(Minimizer &minimizer, const DVec &init, 
+                        const std::vector<const DoubleModel *> &v, Index s);
+    SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init, 
                         const std::vector<const DoubleModel *> &v);
-    SampleFitResult fit(Minimizer &minimizer, const DVec &init,
+    SampleFitResult fit(Minimizer &minimizer, const DVec &init, 
                         const std::vector<const DoubleModel *> &v);
     template <typename... Ts>
     SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,

lib/Statistics/XYStatData.cpp

@@ -358,7 +358,7 @@ FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
         result    = (*m)(chi2);
         totalInit = result;
     }
-    result.corrRangeDb_ = Math::cdr(getFitCorrMat());
+    result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
     result.chi2_       = chi2(result);
     result.nPar_       = nPar;
     result.nDof_       = layout.totalYSize - nPar;

physics/2pt-fit.cpp

@@ -24,7 +24,7 @@ int main(int argc, char *argv[])
 {
     // parse arguments /////////////////////////////////////////////////////////
     OptParser            opt;
-    bool                 parsed, doLaplace, doPlot, doHeatmap, doCorr, fold, doScan, noGuess;
+    bool                 parsed, doLaplace, doPlot, doHeatmap, doCorr, fold, doScan;
     string               corrFileName, model, outFileName, outFmt, savePlot;
     Index                ti, tf, shift, nPar, thinning;
     double               svdTol;
@@ -59,8 +59,6 @@ int main(int argc, char *argv[])
                   "show the fit plot");
     opt.addOption("h", "heatmap"  , OptParser::OptType::trigger, true,
                   "show the fit correlation heatmap");
-    opt.addOption("", "no-guess"  , OptParser::OptType::trigger, true,
-                  "do not try to guess fit parameters");
     opt.addOption("", "save-plot", OptParser::OptType::value, true,
                     "saves the source and .pdf", "");
     opt.addOption("", "scan", OptParser::OptType::trigger, true,
@@ -89,7 +87,6 @@ int main(int argc, char *argv[])
     fold         = opt.gotOption("fold");
     doPlot       = opt.gotOption("p");
     doHeatmap    = opt.gotOption("h");
-    noGuess      = opt.gotOption("no-guess");
     savePlot     = opt.optionValue("save-plot");
     doScan       = opt.gotOption("scan");
     switch (opt.optionValue<unsigned int>("v"))
@@ -117,6 +114,7 @@ int main(int argc, char *argv[])
     nt      = corr[central].rows();
     corr    = corr.block(0, 0, nt, 1);
     corr    = CorrelatorUtils::shift(corr, shift);
+
     if (doLaplace)
     {
         vector<double> filter = {1., -2., 1.};
@@ -158,6 +156,11 @@ int main(int argc, char *argv[])
         }
     }
     
+    if (fold)
+    {
+        corr = CorrelatorUtils::fold(corr,modelPar);
+    }
+
     // fit /////////////////////////////////////////////////////////////////////
     DVec                init(nPar);
     NloptMinimizer      globMin(NloptMinimizer::Algorithm::GN_CRS2_LM);
@@ -170,14 +173,13 @@ int main(int argc, char *argv[])
     fitter.setThinning(thinning);
 
     // set initial values ******************************************************
-    if ((modelPar.type != CorrelatorType::undefined) and !noGuess)
+    if (modelPar.type != CorrelatorType::undefined)
     {
         init = CorrelatorModels::parameterGuess(corr, modelPar);
     }
     else
     {
-        init.fill(1.);
-        init(0) = 0.2;
+        init.fill(0.1);
     }
 
     // set limits for minimisers ***********************************************

utils/sample-dwt.cpp

@@ -17,7 +17,6 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include <LatAnalyze/Core/Math.hpp>
 #include <LatAnalyze/Core/OptParser.hpp>
 #include <LatAnalyze/Core/Plot.hpp>
 #include <LatAnalyze/Io/Io.hpp>
@@ -54,13 +53,7 @@ int main(int argc, char *argv[])
         cerr << "usage: " << argv[0];
         cerr << " <options> <input file>" << endl;
         cerr << endl << "Possible options:" << endl << opt << endl;
-        cerr << "Available DWT filters:" << endl;
-        for (auto &fv: DWTFilters::fromName)
-        {
-            cerr << fv.first << " ";
-        }
-        cerr << endl << endl;
-
+        
         return EXIT_FAILURE;
     }
     inFilename  = opt.getArgs()[0];
@@ -75,45 +68,22 @@ int main(int argc, char *argv[])
     DMatSample            in = Io::load<DMatSample>(inFilename), res;
     Index                 nSample = in.size(), n = in[central].rows();
     vector<DMatSample>    out(ss ? 1 : level, DMatSample(nSample)), 
-                          outh(ss ? 0 : level, DMatSample(nSample)),
-                          concath(ss ? 0 : level, DMatSample(nSample));
-    DMatSample            concat(nSample, n, 1);
+                          outh(ss ? 0 : level, DMatSample(nSample));
     DWT                   dwt(*DWTFilters::fromName.at(filterName));
     vector<DWT::DWTLevel> dataDWT(level);
 
-    FOR_STAT_ARRAY(in, s)
-    {
-        in[s].conservativeResize(n, 1);
-    }
     if (!ss)
     {
-        DMatSample buf(nSample);
-
         cout << "-- compute discrete wavelet transform" << endl;
         cout << "filter '" << filterName << "' / " << level << " level(s)" << endl;
         FOR_STAT_ARRAY(in, s)
         {
-            dataDWT = dwt.forward(in[s], level);
+            dataDWT = dwt.forward(in[s].col(0), level);
             for (unsigned int l = 0; l < level; ++l)
             {
-                out[l][s]     = dataDWT[l].first;
-                outh[l][s]    = dataDWT[l].second;
-                concath[l][s] = DWT::concat(dataDWT, l, true);
+                out[l][s]  = dataDWT[l].first;
+                outh[l][s] = dataDWT[l].second;
             }
-            concat[s] = DWT::concat(dataDWT);
-        }
-        cout << "Data CDR " << Math::cdr(in.correlationMatrix()) << " dB" << endl;
-        cout << "DWT CDR " << Math::cdr(concat.correlationMatrix()) << " dB" << endl;
-        for (unsigned int l = 0; l < level; ++l)
-        {
-            cout << "DWT level " << l << " CDR: L= ";
-            cout << Math::cdr(out[l].correlationMatrix()) << " dB / H= ";
-            cout << Math::cdr(outh[l].correlationMatrix()) << " dB" << endl;
-        }
-        for (unsigned int l = 0; l < level; ++l)
-        {
-            cout << "DWT detail level " << l << " CDR: ";
-            cout << Math::cdr(concath[l].correlationMatrix()) << " dB" << endl;
         }
     }
     else
@@ -132,7 +102,7 @@ int main(int argc, char *argv[])
         }
         FOR_STAT_ARRAY(in, s)
         {
-            dataDWT.back().first = in[s];
+            dataDWT.back().first = in[s].col(0);
             out[0][s] = dwt.backward(dataDWT);
         }
     }
@@ -145,9 +115,7 @@ int main(int argc, char *argv[])
             {
                 Io::save<DMatSample>(out[l], outFilename + "/L" + strFrom(l) + ".h5");
                 Io::save<DMatSample>(outh[l], outFilename + "/H" + strFrom(l) + ".h5");
-                Io::save<DMatSample>(concath[l], outFilename + "/concatH" + strFrom(l) + ".h5");
             }
-            Io::save<DMatSample>(concat, outFilename + "/concat.h5");
         }
         else
         {

utils/sample-fake.cpp

@@ -18,30 +18,42 @@
  */
 
 #include <LatAnalyze/Io/Io.hpp>
+#include <LatAnalyze/Core/OptParser.hpp>
+
 
 using namespace std;
 using namespace Latan;
 
 int main(int argc, char *argv[])
 {
+    OptParser            opt;
     Index  nSample;
     double val, err;
     string outFileName;
 
-    if (argc != 5)
+    opt.addOption("r", "seed"      , OptParser::OptType::value,   true,
+                    "random generator seed (default: random)"); 
+    opt.addOption("", "help"      , OptParser::OptType::trigger, true,
+                  "show this help message and exit");
+
+    bool parsed = opt.parse(argc, argv);
+    if (!parsed or (opt.getArgs().size() != 4) or opt.gotOption("help"))
     {
         cerr << "usage: " << argv[0];
         cerr << " <central value> <error> <nSample> <output file>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
 
         return EXIT_FAILURE;
     }
+
     val         = strTo<double>(argv[1]);
     err         = strTo<double>(argv[2]);
     nSample     = strTo<Index>(argv[3]);
     outFileName = argv[4];
 
     random_device         rd;
-    mt19937               gen(rd());
+    SeedType              seed = (opt.gotOption("r")) ? opt.optionValue<SeedType>("r") : rd();
+    mt19937               gen(seed);
     normal_distribution<> dis(val, err);
     DSample               res(nSample);
 
@@ -59,4 +71,4 @@ int main(int argc, char *argv[])
     Io::save<DSample>(res, outFileName);
 
     return EXIT_SUCCESS;
-}
+}

utils/sample-plot-corr.cpp

@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
     var    = sample.varianceMatrix();
     corr   = sample.correlationMatrix();
 
-    cout << "dynamic range " << Math::cdr(corr) << " dB" << endl;
+    cout << "dynamic range " << Math::svdDynamicRangeDb(corr) << " dB" << endl;
     p << PlotCorrMatrix(corr);
     p.display();
     if (!outVarName.empty())

utils/sample-read.cpp

@@ -38,9 +38,23 @@ int main(int argc, char *argv[])
     {
         DMatSample s    = Io::load<DMatSample>(fileName);
         string     name = Io::getFirstName(fileName);
+        Index nRows     = s[central].rows();
+        Index nCols     = s[central].cols();
         cout << scientific;
-        cout << "central value:\n"      << s[central]               << endl;
-        cout << "standard deviation:\n" << s.variance().cwiseSqrt() << endl;
+        cout << "central value +/-  standard deviation\n" << endl;
+        cout << "Re:" << endl;
+        for(Index i = 0; i < nRows; i++)
+        {
+            cout << s[central](i,0) << " +/- " << s.variance().cwiseSqrt()(i,0) << endl;
+        }
+        if(nCols == 2)
+        {
+            cout << "\nIm:" << endl;
+            for(Index i = 0; i < nRows; i++)
+            {
+                cout << s[central](i,1) << " +/- " << s.variance().cwiseSqrt()(i,1) << endl;
+            }   
+        }
         if (!copy.empty())
         {
             Io::save(s, copy, File::Mode::write, name);
@@ -51,8 +65,8 @@ int main(int argc, char *argv[])
         DSample s    = Io::load<DSample>(fileName);
         string  name = Io::getFirstName(fileName);
         cout << scientific;
-        cout << "central value:\n"      << s[central]         << endl;
-        cout << "standard deviation:\n" << sqrt(s.variance()) << endl;
+        cout << "central value +/-  standard deviation\n" << endl;
+        cout << s[central] << " +/- " << sqrt(s.variance()) << endl;
         if (!copy.empty())
         {
             Io::save(s, copy, File::Mode::write, name);