Merge branch 'release/3.3'

.gitignore

@@ -35,3 +35,7 @@ ci-scripts/local/*
 .idea/*
 CMakeLists.txt
 cmake-build-debug/*
+
+# VS Code Studio stuff
+.vscode
+*.code-workspace

Readme.md

@@ -55,6 +55,28 @@ in the `ci-scripts` directory where `<prefix>` is where you want LatAnalyze (and
 For a more customised installation, one first needs to generate the build system by running `./bootstrap.sh` in the root directory. Then the library can be built and installed through the usual GNU mantra `./configure <options> && make && make install`. Use `./configure --help` to obtain a list of possible options for `./configure`. Because Eigen expressions rely a lot on inlining and compiler optimisations it is strongly recommended to set the `CXXFLAGS` variable to `-O3 -march=native -mtune=native`.
 
 ## History
+#### v3.3
+Additions:
+* Sample plot CL utility.
+* Infinity as a math constant.
+* Option to dump bootstrap sequence while resampling.
+* FFT through the GSL.
+
+Changes:
+* GSL integrator accepts infinite bounds.
+* `latan-sample-combine` accepts mixes of `DSample` and `DMatSample`.
+* More general `latan-sample-element` command.
+
+#### v3.2.2
+Additions:
+* The math interpreter supports `inf` for infinity.
+
+Changes:
+* Vector version of `setUnidimData`.
+
+Fixes:
+* Variance matrix computation fix.
+
 #### v3.2.1
 Fix:
 * Wrong argument number check in `latan-resample`

ci-scripts/install-hdf5.sh

@@ -1,6 +1,6 @@
 #!/usr/bin/env bash
 
-NAME='hdf5-1.8.16'
+NAME='hdf5-1.10.1'
 
 if (( $# != 1 )); then
   echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
@@ -11,7 +11,7 @@ PREFIX=$1
 set -ex
 INITDIR=`pwd`
 cd local/build
-wget http://www.hdfgroup.org/ftp/HDF5/releases/${NAME}/src/${NAME}.tar.gz
+wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/${NAME}/src/${NAME}.tar.gz
 tar -xzvf ${NAME}.tar.gz
 mkdir ${NAME}/build
 cd ${NAME}/build

configure.ac

@@ -2,7 +2,7 @@
 
 # Initialization
 AC_PREREQ([2.63])
-AC_INIT([LatAnalyze],[3.2.1-dev],[antonin.portelli@me.com],[LatAnalyze])
+AC_INIT([LatAnalyze],[3.3],[antonin.portelli@me.com],[LatAnalyze])
 AC_CONFIG_AUX_DIR([.buildutils])
 AC_CONFIG_SRCDIR([lib/Global.cpp])
 AC_CONFIG_SRCDIR([utils/sample_read.cpp])

lib/Dataset.hpp

@@ -47,6 +47,8 @@ public:
     // resampling
     Sample<T> bootstrapMean(const Index nSample, const SeedType seed);
     Sample<T> bootstrapMean(const Index nSample);
+    void      dumpBootstrapSeq(std::ostream &out, const Index nSample,
+                               const SeedType seed);
 private:
     // mean from pointer vector for resampling
     void ptVectorMean(T &m, const std::vector<const T *> &v);
@@ -114,6 +116,23 @@ Sample<T> Dataset<T>::bootstrapMean(const Index nSample)
     return bootstrapMean(nSample, rd());
 }
 
+template <typename T>
+void Dataset<T>::dumpBootstrapSeq(std::ostream &out, const Index nSample,
+                                  const SeedType seed)
+{
+    std::mt19937                         gen(seed);
+    std::uniform_int_distribution<Index> dis(0, this->size() - 1);
+
+    for (Index i = 0; i < nSample; ++i)
+    {
+        for (unsigned int j = 0; j < this->size(); ++j)
+        {
+            out << dis(gen) << " " << std::endl;
+        }
+        out << std::endl;
+    }
+}
+
 template <typename T>
 void Dataset<T>::ptVectorMean(T &m, const std::vector<const T *> &v)
 {

lib/FFT.hpp

@@ -0,0 +1,53 @@
+/*
+ * FFT.hpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2017 Antonin Portelli
+ *
+ * LatAnalyze 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 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.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_FFT_hpp_
+#define Latan_FFT_hpp_
+
+#include <LatAnalyze/Global.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                             FFT abstract class                             *
+ ******************************************************************************/
+
+class FFT
+{
+public:
+    enum
+    {
+      Forward  = 0,
+      Backward = 1
+    };
+public:
+  // constructor
+  FFT(void) = default;
+  FFT(const Index size);
+  // destructor
+  virtual ~FFT(void) = default;
+  // size
+  virtual void resize(const Index size) = 0;
+  // FFT
+  virtual void operator()(CMat &x, const unsigned int dir = FFT::Forward) = 0;
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_FFT_hpp_

lib/GslFFT.cpp

@@ -0,0 +1,89 @@
+/*
+ * GslFFT.cpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2017 Antonin Portelli
+ *
+ * LatAnalyze 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 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.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/GslFFT.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                           GslFFT implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+GslFFT::GslFFT(const Index size)
+{
+    resize(size);
+}
+
+// destructor //////////////////////////////////////////////////////////////////
+GslFFT::~GslFFT(void)
+{
+    clear();
+}
+
+// size ////////////////////////////////////////////////////////////////////////
+void GslFFT::resize(const Index size)
+{
+    if (size_ != size)
+    {
+        size_      = size;
+        wavetable_ = gsl_fft_complex_wavetable_alloc(size_);
+        workspace_ = gsl_fft_complex_workspace_alloc(size_);
+    }
+}
+
+// fft /////////////////////////////////////////////////////////////////////////
+void GslFFT::operator()(CMat &x, const unsigned int dir)
+{
+    if (x.size() != size_)
+    {
+        LATAN_ERROR(Size, "wrong input vector size");
+    }
+    else
+    {
+        switch (dir)
+        {
+            case FFT::Forward:
+                gsl_fft_complex_forward((double *)x.data(), 1, size_,
+                                        wavetable_, workspace_);
+                break;
+            case FFT::Backward:
+                gsl_fft_complex_backward((double *)x.data(), 1, size_,
+                                         wavetable_, workspace_);
+                break;
+            default:
+                LATAN_ERROR(Argument, "invalid FT direction");
+                break;
+        }
+    }
+}
+
+// destroy GSL objects /////////////////////////////////////////////////////////
+void GslFFT::clear(void)
+{
+    if (!wavetable_)
+    {
+        gsl_fft_complex_wavetable_free(wavetable_);
+    }
+    if (!workspace_)
+    {
+        gsl_fft_complex_workspace_free(workspace_);
+    }
+}

lib/GslFFT.hpp

@@ -0,0 +1,57 @@
+/*
+ * GslFFT.hpp, part of LatAnalyze
+ *
+ * Copyright (C) 2013 - 2017 Antonin Portelli
+ *
+ * LatAnalyze 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 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.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_GslFFT_hpp_
+#define Latan_GslFFT_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/Mat.hpp>
+#include <LatAnalyze/FFT.hpp>
+#include <gsl/gsl_fft_complex.h>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                                 GSL FFT                                    *
+ ******************************************************************************/
+
+class GslFFT: public FFT
+{
+public:
+    // constructors
+    GslFFT(void) = default;
+    GslFFT(const Index size);
+    // destructor
+    virtual ~GslFFT(void);
+    // size
+    void resize(const Index size);
+    // fft
+    virtual void operator()(CMat &x, const unsigned int dir = FFT::Forward);
+private:
+    // destroy GSL objects
+    void clear(void);
+private:
+    Index                     size_{0};
+    gsl_fft_complex_wavetable *wavetable_{nullptr};
+    gsl_fft_complex_workspace *workspace_{nullptr};
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_GslFFT_hpp_

lib/GslQagsIntegrator.cpp

@@ -19,6 +19,7 @@
 
 #include <LatAnalyze/GslQagsIntegrator.hpp>
 #include <LatAnalyze/includes.hpp>
+#include <LatAnalyze/Math.hpp>
 
 using namespace std;
 using namespace Latan;
@@ -55,9 +56,26 @@ double GslQagsIntegrator::operator()(const DoubleFunction &f, const double xMin,
     
     gslF.function = fWrap;
     gslF.params   = reinterpret_cast<void *>(&const_cast<DoubleFunction &>(f));
-    
+    if ((xMin > -Math::inf) and (xMax < Math::inf))
-    gsl_integration_qags(&gslF, xMin, xMax, 0.0, precision_, limit_, workspace_,
+    {
-                         &result, &error_);
+        gsl_integration_qags(&gslF, xMin, xMax, 0.0, precision_, limit_,
+                             workspace_, &result, &error_);
+    }
+    else if (xMax < Math::inf)
+    {
+        gsl_integration_qagil(&gslF, xMax, 0.0, precision_, limit_,
+                              workspace_, &result, &error_);
+    }
+    else if (xMin > -Math::inf)
+    {
+        gsl_integration_qagiu(&gslF, xMin, 0.0, precision_, limit_,
+                              workspace_, &result, &error_);
+    }
+    else
+    {
+        gsl_integration_qagi(&gslF, 0.0, precision_, limit_,
+                             workspace_, &result, &error_);
+    }
     
     return result;
 }

lib/Makefile.am

@@ -31,6 +31,7 @@ libLatAnalyze_la_SOURCES = \
     FitInterface.cpp       \
     Function.cpp           \
     Global.cpp             \
+    GslFFT.cpp             \
     GslHybridRootFinder.cpp\
     GslMinimizer.cpp       \
     GslQagsIntegrator.cpp  \
@@ -59,10 +60,12 @@ libLatAnalyze_la_HEADERS = \
     Dataset.hpp            \
     Derivative.hpp         \
     Exceptions.hpp         \
+    FFT.hpp                \
     Function.hpp           \
     File.hpp               \
     FitInterface.hpp       \
     Global.hpp             \
+    GslFFT.hpp             \
     GslHybridRootFinder.hpp\
     GslMinimizer.hpp       \
     GslQagsIntegrator.hpp  \

lib/Math.hpp

@@ -72,8 +72,9 @@ namespace MATH_NAMESPACE
     DMat varToCorr(const DMat &var);
     
     // Constants
-    const double pi = 3.1415926535897932384626433832795028841970;
+    constexpr double pi  = 3.1415926535897932384626433832795028841970;
-    const double e  = 2.7182818284590452353602874713526624977572;
+    constexpr double e   = 2.7182818284590452353602874713526624977572;
+    constexpr double inf = std::numeric_limits<double>::infinity();
 }
 
 

lib/MathInterpreter.cpp

@@ -690,7 +690,7 @@ void MathInterpreter::compile(RunContext &context)
         if (root_)
         {
             context.addVariable("pi",  Math::pi);
-            context.addVariable("inf", HUGE_VAL);
+            context.addVariable("inf", Math::inf);
             ADD_STDMATH_FUNCS(context);
             root_->compile(program_, context);
             for (unsigned int i = 0; i < program_.size(); ++i)

lib/TabFunction.cpp

@@ -29,18 +29,11 @@ using namespace Latan;
 // constructors ////////////////////////////////////////////////////////////////
 TabFunction::TabFunction(const DVec &x, const DVec &y,
                          const InterpType interpType)
-: interpType_(interpType)
 {
     setData(x, y);
+    setInterpolationType(interpType);
 }
 
-//TabFunction::TabFunction(const XYStatData &data, const Index i, const Index j,
-//                         const InterpType interpType)
-//: interpType_(interpType)
-//{
-//    setData(data, i, j);
-//}
-
 // access //////////////////////////////////////////////////////////////////////
 void TabFunction::setData(const DVec &x, const DVec &y)
 {
@@ -54,10 +47,10 @@ void TabFunction::setData(const DVec &x, const DVec &y)
     }
 }
 
-//void TabFunction::setData(const XYStatData &data, const Index i, const Index j)
+void TabFunction::setInterpolationType(const InterpType interpType)
-//{
+{
-//    setData(data.x(i), data.y(j));
+    interpType_ = interpType;
-//}
+}
 
 // function call ///////////////////////////////////////////////////////////////
 double TabFunction::operator()(const double *arg) const
@@ -154,12 +147,6 @@ DoubleFunction Latan::interpolate(const DVec &x, const DVec &y,
     return TabFunction(x, y, interpType).makeFunction();
 }
 
-//DoubleFunction Latan::interpolate(const XYStatData &data, const Index i,
-//                                  const Index j, const InterpType interpType)
-//{
-//    return TabFunction(data, i, j, interpType).makeFunction();
-//}
-
 map<double, double>::const_iterator TabFunction::nearest(const double x) const
 {
     map<double, double>::const_iterator ret;

lib/TabFunction.hpp

@@ -45,13 +45,11 @@ public:
     TabFunction(void) = default;
     TabFunction(const DVec &x, const DVec &y,
                 const InterpType interpType = InterpType::LINEAR);
-    //TabFunction(const XYStatData &data, const Index i = 0, const Index j = 0,
-    //            const InterpType interpType = InterpType::LINEAR);
     // destructor
     virtual ~TabFunction(void) = default;
     // access
     void setData(const DVec &x, const DVec &y);
-    // void setData(const XYStatData &data, const Index i = 0, const Index j = 0);
+    void setInterpolationType(const InterpType interpType);
     // function call
     double operator()(const double *arg) const;
     // factory
@@ -65,9 +63,6 @@ private:
 
 DoubleFunction interpolate(const DVec &x, const DVec &y,
                            const InterpType interpType = InterpType::LINEAR);
-//DoubleFunction interpolate(const XYStatData &data, const Index i = 0,
-//                           const Index j = 0,
-//                           const InterpType interpType = InterpType::LINEAR);
 
 END_LATAN_NAMESPACE
 

physics/Makefile.am

@@ -6,14 +6,10 @@ if CXX_INTEL
 endif
 endif
 
-bin_PROGRAMS = latan-fit-2pt latan-fit-phys
+bin_PROGRAMS = latan-2pt-fit
 
-latan_fit_2pt_SOURCES  = fit-2pt.cpp
+latan_2pt_fit_SOURCES  = 2pt-fit.cpp
-latan_fit_2pt_CXXFLAGS = $(COM_CXXFLAGS)
+latan_2pt_fit_CXXFLAGS = $(COM_CXXFLAGS)
-latan_fit_2pt_LDFLAGS  = -L../lib/.libs -lLatAnalyze
+latan_2pt_fit_LDFLAGS  = -L../lib/.libs -lLatAnalyze
-
-latan_fit_phys_SOURCES  = fit-phys.cpp fit-phys-env.cpp
-latan_fit_phys_CXXFLAGS = $(COM_CXXFLAGS)
-latan_fit_phys_LDFLAGS  = -L../lib/.libs -lLatAnalyze
 
 ACLOCAL_AMFLAGS = -I .buildutils/m4

physics/fit-phys-env.cpp

@@ -1,509 +0,0 @@
-#include "fit-phys-env.hpp"
-#include <LatCore/XmlReader.hpp>
-#include <LatAnalyze/CompiledModel.hpp>
-#include <LatAnalyze/Io.hpp>
-
-#define DRATIO(a,b) static_cast<double>(a)/static_cast<double>(b)
-
-using namespace std;
-using namespace Latan;
-
-void FitEnv::reset(void)
-{
-    nT_.clear();
-    nL_.clear();
-    variable_.clear();
-    varData_.clear();
-    varName_.clear();
-    varScalePow_.clear();
-    quantity_.clear();
-    quData_.clear();
-    quName_.clear();
-    ensemble_.clear();
-    point_.clear();
-    macro_.clear();
-    scaleVar_ = nullptr;
-}
-
-Index FitEnv::getVarIndex(const string name)
-{
-    if (name == "nT")
-    {
-        return 0;
-    }
-    else if (name == "nL")
-    {
-        return 1;
-    }
-    else
-    {
-        auto it = variable_.find(name);
-        
-        if (it != variable_.end())
-        {
-            return it->second.index;
-        }
-        else
-        {
-            LATAN_ERROR(Range, "no variable with name '" + name + "'");
-        }
-    }
-}
-
-string FitEnv::getVarName(const Index i)
-{
-    if (i < static_cast<Index>(varName_.size()))
-    {
-        return varName_[i];
-    }
-    else
-    {
-        LATAN_ERROR(Range, "no variable with index " + strFrom(i));
-    }
-}
-
-Index FitEnv::getQuIndex(const string name)
-{
-    auto it = quantity_.find(name);
-    
-    if (it != quantity_.end())
-    {
-        return it->second.index;
-    }
-    else
-    {
-        LATAN_ERROR(Range, "no quantity with name '" + name + "'");
-    }
-}
-
-string FitEnv::getQuName(const Index i)
-{
-    if (i < static_cast<Index>(quName_.size()))
-    {
-        return quName_[i];
-    }
-    else
-    {
-        LATAN_ERROR(Range, "no variable with index " + strFrom(i));
-    }
-}
-
-DVec FitEnv::getPhyPt(void)
-{
-    DVec phyPt(varName_.size());
-    
-    for (unsigned int i = 0; i < varName_.size(); ++i)
-    {
-        phyPt(i) = variable_[varName_[i]].physVal;
-    }
-    
-    return phyPt;
-}
-
-vector<const DoubleModel *> FitEnv::getModels(void)
-{
-    vector<const DoubleModel *> m;
-    
-    for (auto &q: quantity_)
-    {
-        m.push_back(&q.second.model);
-    }
-    
-    return m;
-}
-
-#define XGFV(type, ...) XmlReader::getFirstValue<type>(node, __VA_ARGS__)
-
-void FitEnv::parseXml(const string paramFileName)
-{
-    XmlReader                paramFile(paramFileName);
-    const XmlNode            *node = nullptr;
-    set<unsigned int>        nTs, nLs;
-    map<string, set<string>> varFileNames, quFileNames;
-    
-    reset();
-    nSample_ = paramFile.getFirstValue<Index>("nSample");
-    scale_   = paramFile.getFirstValue<string>("scale");
-    
-    // macros
-    if (paramFile.hasNode("macros", "macro"))
-    {
-        node = paramFile.getFirstNode("macros", "macro");
-        while (node)
-        {
-            macro_[XGFV(string, "symbol")] = XGFV(string, "value");
-            node                           = paramFile.getNextSameNode(node);
-        }
-    }
-    
-    // ensembles
-    node = paramFile.getFirstNode("ensembles", "ensemble");
-    while (node)
-    {
-        string   name, spacing;
-        Ensemble ens;
-        
-        name            = XGFV(string, "name");
-        ens.nT          = XGFV(unsigned int, "nT");
-        ens.nL          = XGFV(unsigned int, "nL");
-        ensemble_[name] = ens;
-        node            = paramFile.getNextSameNode(node);
-        nTs.insert(ens.nT);
-        nLs.insert(ens.nL);
-    }
-    
-    // fit variables
-    {
-        string  name;
-        VarInfo var;
-        
-        name            = "T";
-        var.physVal     = HUGE_VAL;
-        var.dim         = -1;
-        variable_[name] = var;
-    }
-    {
-        string  name;
-        VarInfo var;
-        
-        name            = "L";
-        var.physVal     = HUGE_VAL;
-        var.dim         = -1;
-        variable_[name] = var;
-    }
-    node = paramFile.getFirstNode("variables", "variable");
-    while (node)
-    {
-        string  name;
-        VarInfo var;
-        
-        name            = XGFV(string, "name");
-        var.physVal     = XGFV(double, "physical");
-        var.dim         = XGFV(int, "dim");
-        variable_[name] = var;
-        if (name == scale_)
-        {
-            scaleVar_ = &(variable_[name]);
-        }
-        node = paramFile.getNextSameNode(node);
-    }
-    if (!scaleVar_)
-    {
-        LATAN_ERROR(Definition, "scaling variable '" + scale_
-                    + "' not defined");
-    }
-    for (auto &v: variable_)
-    {
-        v.second.index = varName_.size();
-        varName_.push_back(v.first);
-    }
-    for (auto &v: variable_)
-    {
-        varScalePow_.push_back(DRATIO(v.second.dim, scaleVar_->dim));
-    }
-    
-    // fitted quantities
-    node = paramFile.getFirstNode("quantities", "quantity");
-    while (node)
-    {
-        string           name, code = "";
-        Index            nPar;
-        QuInfo           q;
-        DoubleModel      m;
-        shared_ptr<DVec> buf(new DVec(varName_.size()));
-        
-        name            = XGFV(string, "name");
-        nPar            = XGFV(Index, "model", "nPar");
-        q.dim           = XGFV(int, "dim");
-        for (auto &v: variable_)
-        {
-            code += v.first + " = x_" + strFrom(v.second.index) + "; ";
-            code += v.first + "_phy = " + strFrom(v.second.physVal) + "; ";
-        }
-        code           += XGFV(string, "model", "code");
-        DEBUG_VAR(code);
-        m               = compile(code, variable_.size(), nPar);
-        auto wrap       = [m, buf, this, q](const double *x, const double *p)
-        {
-            double s = x[scaleVar_->index];
-            
-            for (unsigned int i = 0; i < varScalePow_.size(); ++i)
-            {
-                if (i != scaleVar_->index)
-                {
-                    (*buf)(i) = x[i]*pow(s, varScalePow_[i]);
-                }
-                else
-                {
-                    (*buf)(i) = x[i];
-                }
-            }
-            
-            return pow(s, -DRATIO(q.dim, scaleVar_->dim))*m(buf->data(), p);
-        };
-        q.model.setFunction(wrap, m.getNArg(), m.getNPar());
-        quantity_[name] = q;
-        node            = paramFile.getNextSameNode(node);
-    }
-    for (auto &q: quantity_)
-    {
-        q.second.index = quName_.size();
-        quName_.push_back(q.first);
-    }
-    
-    // data points
-    node = paramFile.getFirstNode("points", "point");
-    while (node)
-    {
-        string ensemble, fileName;
-        Point  point;
-        
-        ensemble = XGFV(string, "ensemble");
-        auto it  = ensemble_.find(ensemble);
-        if (it == ensemble_.end())
-        {
-            LATAN_ERROR(Parsing, "unknown ensemble '" + ensemble + "'");
-        }
-        macro_["_ensemble_"] = ensemble;
-        point.isActive       = XGFV(bool, "active");
-        point.ensemble       = &(it->second);
-        for (auto &v: variable_)
-        {
-            if (v.first == "T")
-            {
-                fileName = strFrom(point.ensemble->nT);
-            }
-            else if (v.first == "L")
-            {
-                fileName = strFrom(point.ensemble->nL);
-            }
-            else
-            {
-                fileName = macroSubst(XGFV(string, v.first));
-            }
-            point.fileName[v.first] = fileName;
-            varFileNames[v.first].insert(fileName);
-        }
-        for (auto &q: quantity_)
-        {
-            fileName                = macroSubst(XGFV(string, q.first));
-            point.fileName[q.first] = fileName;
-            quFileNames[q.first].insert(fileName);
-        }
-        point_.push_back(point);
-        node = paramFile.getNextSameNode(node);
-    }
-    macro_.erase("_ensemble_");
-    
-    // compute data indices
-    for (auto &v: varFileNames)
-    {
-        varData_.push_back(vector<Data>());
-        for (auto &f: v.second)
-        {
-            Data d;
-            
-            d.fileName            = f;
-            varIndex_[v.first][f] = varData_.back().size();
-            varData_.back().push_back(d);
-        }
-    }
-    for (auto &q: quFileNames)
-    {
-        quData_.push_back(vector<Data>());
-        for (auto &f: q.second)
-        {
-            Data d;
-            
-            d.fileName           = f;
-            quIndex_[q.first][f] = quData_.back().size();
-            quData_.back().push_back(d);
-        }
-    }
-    
-    // compute point coordinates
-    for (auto &p: point_)
-    {
-        p.coord.resize(varName_.size());
-        for (unsigned int i = 0; i < varName_.size(); ++i)
-        {
-            p.coord[i] = varIndex_[varName_[i]][p.fileName[varName_[i]]];
-        }
-    }
-}
-
-#undef XGFV
-
-std::string FitEnv::macroSubst(const std::string str) const
-{
-    std::string res = str, pat;
-    
-    for (auto &m: macro_)
-    {
-        pat      = "@" + m.first + "@";
-        auto pos = res.find(pat);
-        
-        if (pos != string::npos)
-        {
-            res.replace(pos, pat.size(), m.second);
-        }
-    }
-    
-    return res;
-}
-
-void FitEnv::load(void)
-{
-    for (unsigned int i = 0; i < varName_.size(); ++i)
-    {
-        auto &v = varData_[i];
-        
-        if ((varName_[i] == "T") or (varName_[i] == "L"))
-        {
-            for (auto &d: v)
-            {
-                d.value.resize(nSample_);
-                d.value.fill(strTo<double>(d.fileName));
-            }
-        }
-        else
-        {
-            for (auto &d: v)
-            {
-                d.value = Io::load<DSample>(d.fileName);
-                if (d.value.size() != nSample_)
-                {
-                    LATAN_ERROR(Size, "sample loaded from file '" + d.fileName
-                                + "' has a wrong number of element (expected "
-                                + strFrom(nSample_) + ", got "
-                                + strFrom(d.value.size()) + ")");
-                }
-            }
-        }
-        
-    }
-    for (auto &q: quData_)
-    {
-        for (auto &d: q)
-        {
-            d.value = Io::load<DSample>(d.fileName);
-            if (d.value.size() != nSample_)
-            {
-                LATAN_ERROR(Size, "sample loaded from file '" + d.fileName
-                            + "' has a wrong number of element (expected "
-                            + strFrom(nSample_) + ", got "
-                            + strFrom(d.value.size()) + ")");
-            }
-        }
-    }
-
-}
-
-XYSampleData FitEnv::generateData(const bool phyUnits, const bool corr)
-{
-    XYSampleData  data(nSample_);
-    Index         k, k1, k2, ind;
-    const Index   sInd = getVarIndex(scale_);
-    DSample       scale, tmp;
-    int           dim;
-    const int     sDim = scaleVar_->dim;
-    
-    // add dimensions
-    for (unsigned int i = 0; i < varName_.size(); ++i)
-    {
-        data.addXDim(varData_[i].size(), varName_[i],
-                     ((varName_[i] == "T") or (varName_[i] == "L")));
-    }
-    for (auto &q: quName_)
-    {
-        data.addYDim(q);
-    }
-    // add data
-    for (auto &p: point_)
-    {
-        k     = data.dataIndex(p.coord);
-        scale = varData_[sInd][varIndex_[scale_][p.fileName[scale_]]].value;
-        for (unsigned int i = 0; i < varName_.size(); ++i)
-        {
-            ind = varIndex_[varName_[i]][p.fileName[varName_[i]]];
-            dim = variable_[varName_[i]].dim;
-            tmp = varData_[i][ind].value;
-            if (phyUnits and (varName_[i] != scale_))
-            {
-                FOR_STAT_ARRAY(tmp, s)
-                {
-                    tmp[s] *= pow(scale[s], DRATIO(dim, sDim));
-                }
-            }
-            data.x(p.coord[i], i) = tmp;
-        }
-        for (unsigned int j = 0; j < quName_.size(); ++j)
-        {
-            ind = quIndex_[quName_[j]][p.fileName[quName_[j]]];
-            dim = quantity_[quName_[j]].dim;
-            tmp = quData_[j][ind].value;
-            if (phyUnits)
-            {
-                FOR_STAT_ARRAY(tmp, s)
-                {
-                    tmp[s] *= pow(scale[s], DRATIO(dim, sDim));
-                }
-            }
-            data.y(k, j) = tmp;
-        }
-    }
-    // add correlations
-    if (corr)
-    {
-        for (unsigned int p1 = 0; p1 < point_.size(); ++p1)
-        for (unsigned int p2 = p1; p2 < point_.size(); ++p2)
-        {
-            if (point_[p1].ensemble == point_[p2].ensemble)
-            {
-                k1 = data.dataIndex(point_[p1].coord);
-                k2 = data.dataIndex(point_[p2].coord);
-                for (unsigned int i1 = 0; i1 < varName_.size(); ++i1)
-                for (unsigned int i2 = i1; i2 < varName_.size(); ++i2)
-                {
-                    data.assumeXXCorrelated(true, point_[p1].coord[i1], i1,
-                                            point_[p2].coord[i2], i2);
-                }
-                for (unsigned int j1 = 0; j1 < quName_.size(); ++j1)
-                for (unsigned int j2 = j1; j2 < quName_.size(); ++j2)
-                {
-                    data.assumeYYCorrelated(true, k1, j1, k2, j2);
-                }
-                for (unsigned int i = 0; i < varName_.size(); ++i)
-                for (unsigned int j = 0; j < quName_.size(); ++j)
-                {
-                    data.assumeXYCorrelated(true, point_[p1].coord[i], i, k2, j);
-                    data.assumeXYCorrelated(true, point_[p2].coord[i], i, k1, j);
-                }
-            }
-        }
-    }
-    
-    return data;
-}
-
-ostream & operator<<(ostream &out, FitEnv &f)
-{
-    for (unsigned int i = 0; i < f.varName_.size(); ++i)
-    {
-        out << f.varName_[i] << ":" << endl;
-        for (auto &d: f.varData_[i])
-        {
-            out << "  * " << d.fileName << endl;
-        }
-    }
-    for (unsigned int i = 0; i < f.quName_.size(); ++i)
-    {
-        out << f.quName_[i] << ":" << endl;
-        for (auto &d: f.quData_[i])
-        {
-            out << "  * " << d.fileName << endl;
-        }
-    }
-    
-    return out;
-}

physics/fit-phys-env.hpp

@@ -1,82 +0,0 @@
-#ifndef fit_phys_env_hpp_
-#define fit_phys_env_hpp_
-
-#include <LatAnalyze/MatSample.hpp>
-#include <LatAnalyze/Model.hpp>
-#include <LatAnalyze/XYSampleData.hpp>
-
-class FitEnv
-{
-public:
-    // fit variable info
-    struct VarInfo
-    {
-        double       physVal;
-        int          dim;
-        Latan::Index index;
-    };
-    // fitted quantity info
-    struct QuInfo
-    {
-        Latan::DoubleModel model;
-        int                dim;
-        Latan::Index       index;
-    };
-    // ensemble
-    struct Ensemble
-    {
-        unsigned int  nT, nL;
-    };
-    // point
-    struct Point
-    {
-        bool                               isActive{true};
-        const Ensemble                     *ensemble{nullptr};
-        std::map<std::string, std::string> fileName;
-        std::vector<Latan::Index>          coord;
-    };
-    // data container
-    struct Data
-    {
-        std::string    fileName;
-        Latan::DSample value;
-    };
-    // table types
-    typedef std::vector<std::vector<Data>> DataTable;
-    typedef std::map<std::string, std::map<std::string, unsigned int>>
-            IndexTable;
-public:
-    FitEnv(void)          = default;
-    virtual ~FitEnv(void) = default;
-    void                reset(void);
-    Latan::Index        getVarIndex(const std::string name);
-    std::string         getVarName(const Latan::Index i);
-    Latan::Index        getQuIndex(const std::string name);
-    std::string         getQuName(const Latan::Index i);
-    Latan::DVec         getPhyPt(void);
-    std::vector<const Latan::DoubleModel *> getModels(void);
-    void                parseXml(const std::string paramFileName);
-    std::string         macroSubst(const std::string str) const;
-    void                load(void);
-    Latan::XYSampleData generateData(const bool phyUnits, const bool corr);
-    friend std::ostream & operator<<(std::ostream &out, FitEnv &f);
-private:
-    Latan::Index                       nSample_;
-    std::string                        scale_;
-    std::vector<unsigned int>          nT_, nL_;
-    DataTable                          varData_, quData_;
-    IndexTable                         varIndex_, quIndex_;
-    std::map<std::string, VarInfo>     variable_;
-    VarInfo                            *scaleVar_{nullptr};
-    std::vector<std::string>           varName_;
-    std::vector<double>                varScalePow_;
-    std::map<std::string, QuInfo>      quantity_;
-    std::vector<std::string>           quName_;
-    std::map<std::string, Ensemble>    ensemble_;
-    std::vector<Point>                 point_;
-    std::map<std::string, std::string> macro_;
-};
-
-std::ostream & operator<<(std::ostream &out, FitEnv &f);
-
-#endif // fit_phys_env_hpp_

physics/fit-phys.cpp

@@ -1,77 +0,0 @@
-#include <LatAnalyze/Io.hpp>
-#include <LatAnalyze/MinuitMinimizer.hpp>
-#include <LatAnalyze/NloptMinimizer.hpp>
-#include <LatAnalyze/Plot.hpp>
-#include "fit-phys-env.hpp"
-
-using namespace std;
-using namespace Latan;
-
-int main(int argc, char *argv[])
-{
-    // parse arguments /////////////////////////////////////////////////////////
-    string paramFileName;
-    
-    if (argc != 2)
-    {
-        cerr << "usage: " << argv[0] << " <parameter file>" << endl;
-        
-        return EXIT_FAILURE;
-    }
-    paramFileName = argv[1];
-    
-    // parse XML & load data ///////////////////////////////////////////////////
-    FitEnv env;
-    
-    env.parseXml(paramFileName);
-    env.load();
-    
-    XYSampleData uncorrData = env.generateData(false, false);
-    XYSampleData corrData   = env.generateData(false, true);
-    
-    cout << "DATA SUMMARY" << endl;
-    cout << "============" << endl;
-    cout << env << uncorrData << endl;
-    
-    // fit /////////////////////////////////////////////////////////////////////
-    auto v = env.getModels();
-    SampleFitResult fit;
-    MinuitMinimizer min1, min2;
-    vector<Minimizer *> min{&min1, &min2};
-    DVec            init(v[0]->getNPar());
-    
-    min1.setVerbosity(Minimizer::Verbosity::Normal);
-    min2.setVerbosity(Minimizer::Verbosity::Normal);
-    min1.setMaxIteration(1000000);
-    min1.setPrecision(1.0e-3);
-    min2.setMaxIteration(1000000);
-    min2.setPrecision(1.0e-5);
-    init.fill(1.0);
-    fit = uncorrData.fit(min, init, v);
-    fit.print();
-    init = fit[central].block(0, 0, init.size(), 1);
-    fit = corrData.fit(min2, init, v);
-    fit.print();
-    
-//    init = fit[central].block(0, 0, v[0]->getNPar(), 1);
-//    min1.setVerbosity(Minimizer::Verbosity::Normal);
-//    fit = corrData.fit(min1, init, v);
-    
-    // plot ////////////////////////////////////////////////////////////////////
-//    Plot p;
-//    DVec phyPt = env.getPhyPt();
-//    phyPt(env.getVarIndex("a")) = 1.;
-//    XYSampleData projData = uncorrData.getPartialResiduals(fit, phyPt, env.getVarIndex("M_Ds"));
-//    
-//    p << PlotPredBand(fit.getModel(_).bind(env.getVarIndex("M_Ds"), phyPt), 0., 3.);
-//    p << PlotData(projData.getData(), env.getVarIndex("M_Ds"), 0);
-//    p.display();
-//    p.reset();
-//    projData = uncorrData.getPartialResiduals(fit, phyPt, env.getVarIndex("a"));
-//    p << PlotPredBand(fit.getModel(_).bind(env.getVarIndex("a"), phyPt), 0., 1.);
-//    p << PlotData(projData.getData(), env.getVarIndex("a"), 0);
-//    p.display();
-//    p.reset();
-    
-    return EXIT_SUCCESS;
-}

utils/Makefile.am

@@ -10,6 +10,8 @@ bin_PROGRAMS =            \
     latan-sample-combine  \
     latan-sample-element  \
     latan-sample-fake     \
+    latan-sample-ft       \
+    latan-sample-plot     \
     latan-sample-plot-corr\
     latan-sample-read     \
     latan-resample
@@ -26,10 +28,18 @@ latan_sample_fake_SOURCES  = sample-fake.cpp
 latan_sample_fake_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_fake_LDFLAGS  = -L../lib/.libs -lLatAnalyze
 
+latan_sample_ft_SOURCES  = sample-ft.cpp
+latan_sample_ft_CXXFLAGS = $(COM_CXXFLAGS)
+latan_sample_ft_LDFLAGS  = -L../lib/.libs -lLatAnalyze
+
 latan_sample_plot_corr_SOURCES  = sample-plot-corr.cpp
 latan_sample_plot_corr_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_plot_corr_LDFLAGS  = -L../lib/.libs -lLatAnalyze
 
+latan_sample_plot_SOURCES  = sample-plot.cpp
+latan_sample_plot_CXXFLAGS = $(COM_CXXFLAGS)
+latan_sample_plot_LDFLAGS  = -L../lib/.libs -lLatAnalyze
+
 latan_sample_read_SOURCES  = sample-read.cpp
 latan_sample_read_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_read_LDFLAGS  = -L../lib/.libs -lLatAnalyze

utils/resample.cpp

@@ -39,7 +39,7 @@ int main(int argc, char *argv[])
 {
     // argument parsing ////////////////////////////////////////////////////////
     OptParser     opt;
-    bool          parsed;
+    bool          parsed, dumpBoot;
     random_device rd;
     SeedType      seed = rd();
     string        manFileName, nameFileName, outDirName;
@@ -56,6 +56,8 @@ int main(int argc, char *argv[])
                   "output directory", ".");
     opt.addOption("f", "format"    , OptParser::OptType::value,   true,
                   "output file format", DEF_FMT);
+    opt.addOption("d", "dump-boot" , OptParser::OptType::trigger, true,
+                  "dump bootstrap sequence");
     opt.addOption("" , "help"      , OptParser::OptType::trigger, true,
                   "show this help message and exit");
     parsed = opt.parse(argc, argv);
@@ -75,6 +77,7 @@ int main(int argc, char *argv[])
     }
     ext          = opt.optionValue("f");
     outDirName   = opt.optionValue("o");
+    dumpBoot     = opt.gotOption("d");
     manFileName  = opt.getArgs()[0];
     nameFileName = opt.getArgs()[1];
     
@@ -124,6 +127,15 @@ int main(int argc, char *argv[])
         
         cout << '\r' << ProgressBar(i + 1, name.size());
         data[name[i]].bin(binSize);
+        if ((i == 0) and dumpBoot)
+        {
+            ofstream file(outDirName + "/" + manFileName + ".bootseq");
+
+            file << "# bootstrap sequences" << endl;
+            file << "# manifest file: " << manFileName << endl;
+            file << "#      bin size: " << binSize << endl;
+            data[name[i]].dumpBootstrapSeq(file, nSample, seed);
+        }
         s = data[name[i]].bootstrapMean(nSample, seed);
         Io::save<DMatSample>(s, outDirName + "/" + outFileName,
                              File::Mode::write, outFileName);

utils/sample-combine.cpp

@@ -25,14 +25,21 @@ using namespace std;
 using namespace Latan;
 
 template <typename T>
-static void loadAndCheck(vector<T> &sample, const vector<string> &fileName)
+static void loadAndCheck(vector<T> &sample __dumb,
+                         const vector<string> &fileName __dumb)
 {
-    const unsigned int n = sample.size();
+    abort();
+}
+
+template <>
+void loadAndCheck(vector<DSample> &sample, const vector<string> &fileName)
+{
+    const unsigned int n = static_cast<unsigned int>(sample.size());
     Index              nSample = 0;
     
     for (unsigned int i = 0; i < n; ++i)
     {
-        sample[i] = Io::load<T>(fileName[i]);
+        sample[i] = Io::load<DSample>(fileName[i]);
         if (i == 0)
         {
             nSample = sample[i].size();
@@ -46,6 +53,62 @@ static void loadAndCheck(vector<T> &sample, const vector<string> &fileName)
     }
 }
 
+template <>
+void loadAndCheck(vector<DMatSample> &sample, const vector<string> &fileName)
+{
+    const unsigned int n = static_cast<unsigned int>(sample.size());
+    Index              nSample = 0;
+    set<unsigned int>  failed;
+    bool               gotSize = false;
+    Index              nRow = 0, nCol = 0;
+    
+    for (unsigned int i = 0; i < n; ++i)
+    {
+        try
+        {
+            sample[i] = Io::load<DMatSample>(fileName[i]);
+            if (!gotSize)
+            {
+                nSample = sample[i].size();
+                nRow    = sample[i][central].rows();
+                nCol    = sample[i][central].cols();
+                gotSize = true;
+            }
+        }
+        catch (Exceptions::Definition)
+        {
+            failed.insert(i);
+        }
+    }
+    for (unsigned int i: failed)
+    {
+        DSample buf;
+        
+        buf = Io::load<DSample>(fileName[i]);
+        sample[i].resize(nSample);
+        FOR_STAT_ARRAY(sample[i], s)
+        {
+            sample[i][s] = DMat::Constant(nRow, nCol, buf[s]);
+        }
+    }
+    for (unsigned int i = 0; i < n; ++i)
+    {
+        if (sample[i].size() != nSample)
+        {
+            cerr << "error: number of sample mismatch (between '";
+            cerr << fileName[0] << "' and '" << fileName[i] << "')" << endl;
+            abort();
+        }
+        if ((sample[i][central].rows() != nRow) and
+            (sample[i][central].cols() != nCol))
+        {
+            cerr << "error: matrix size mismatch (between '";
+            cerr << fileName[0] << "' and '" << fileName[i] << "')" << endl;
+            abort();
+        }
+    }
+}
+
 template <typename T>
 static void combine(const string &outFileName __dumb,
                     const vector<T> &sample __dumb, const string &code __dumb)
@@ -57,7 +120,7 @@ template <>
 void combine(const string &outFileName, const vector<DSample> &sample,
              const string &code)
 {
-    const unsigned int n = sample.size();
+    const unsigned int n = static_cast<unsigned int>(sample.size());
     DoubleFunction     f = compile(code, n);
     DSample            result(sample[0]);
     DVec               buf(n);
@@ -87,7 +150,7 @@ template <>
 void combine(const string &outFileName, const vector<DMatSample> &sample,
              const string &code)
 {
-    const unsigned int n = sample.size();
+    const unsigned int n = static_cast<unsigned int>(sample.size());
     DoubleFunction     f = compile(code, n);
     DVec               buf(n);
     DMatSample         result(sample[0]);

utils/sample-element.cpp

@@ -1,39 +1,103 @@
-#include <iostream>
+/*
-#include <string>
+ * sample-element.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
+ *
+ * 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 <LatCore/OptParser.hpp>
 #include <LatAnalyze/Io.hpp>
 #include <LatAnalyze/MatSample.hpp>
 
-
-int main(int argc, char* argv[])
-{
 using namespace std;
 using namespace Latan;
 
-    if (argc != 4 and argc != 5) {
+int main(int argc, char* argv[])
-        cout << "Usage: " << argv[0] << " <input filename> <row> <column> ";
+{
-        cout << "[output filename]" << endl;
+    // argument parsing ////////////////////////////////////////////////////////
-        return -1;
+    OptParser opt;
+    bool      parsed;
+    string    inFilename, outFilename;
+    Index     r, c, nr, nc;
+    
+    opt.addOption("r", "row",    OptParser::OptType::value  , false,
+                  "row");
+    opt.addOption("c", "col",    OptParser::OptType::value  , false,
+                  "column");
+    opt.addOption("" , "nrow",   OptParser::OptType::value  , true,
+                  "number of rows (default: 1)", "1");
+    opt.addOption("" , "ncol",   OptParser::OptType::value  , true,
+                  "number of columns (default: 1)", "1");
+    opt.addOption("o", "output", OptParser::OptType::value  , true,
+                  "output file name (default: result not saved)", "");
+    opt.addOption("" , "help"  , OptParser::OptType::trigger, true,
+                  "show this help message and exit");
+    parsed = opt.parse(argc, argv);
+    if (!parsed or (opt.getArgs().size() != 1) or opt.gotOption("help"))
+    {
+        cerr << "usage: " << argv[0];
+        cerr << " <options> <input file>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
+        
+        return EXIT_FAILURE;
     }
+    inFilename  = opt.getArgs()[0];
+    outFilename = opt.optionValue("o");
+    r           = opt.optionValue<Index>("r");
+    c           = opt.optionValue<Index>("c");
+    nr          = opt.optionValue<Index>("nrow");
+    nc          = opt.optionValue<Index>("ncol");
     
-    string inFileName = argv[1];
+    // Data extraction /////////////////////////////////////////////////////////
-    auto row = strTo<Index>(argv[2]);
+    auto inputData = Io::load<DMatSample>(inFilename);
-    auto col = strTo<Index>(argv[3]);
-    string outFileName = (argc == 5) ? argv[4] : "";
     
-    auto inputData = Io::load<DMatSample>(inFileName);
+    if ((nr == 1) and (nc == 1))
-    cout << scientific;
-    cout << "central value:\n" << inputData[central](row, col) << endl;
-    cout << "standard deviation:\n";
-    cout << inputData.variance().cwiseSqrt()(row, col) << endl;
-
-    if (not outFileName.empty())
     {
         DSample outputData(inputData.size());
-        FOR_STAT_ARRAY(inputData, s) {
+        
-            outputData[s] = inputData[s](row, col);
+        FOR_STAT_ARRAY(inputData, s)
+        {
+            outputData[s] = inputData[s](r, c);
         }
         
-        Io::save(outputData, outFileName);
+        cout << scientific;
+        cout << "central value:\n" << outputData[central] << endl;
+        cout << "standard deviation:\n";
+        cout << sqrt(outputData.variance()) << endl;
+        if (not outFilename.empty())
+        {
+            Io::save(outputData, outFilename);
         }
     }
+    else
+    {
+        DMatSample outputData(inputData.size(), nr, nc);
+        
+        FOR_STAT_ARRAY(inputData, s)
+        {
+            outputData[s] = inputData[s].block(r, c, nr, nc);
+        }
+        cout << scientific;
+        cout << "central value:\n" << outputData[central] << endl;
+        cout << "standard deviation:\n";
+        cout << outputData.variance().cwiseSqrt() << endl;
+        if (not outFilename.empty())
+        {
+            Io::save(outputData, outFilename);
+        }
+    }
+
+    return EXIT_SUCCESS;
+}

utils/sample-ft.cpp

@@ -0,0 +1,94 @@
+/*
+ * sample-ft.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2016 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 <LatCore/OptParser.hpp>
+#include <LatAnalyze/GslFFT.hpp>
+#include <LatAnalyze/Io.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char *argv[])
+{
+    // argument parsing ////////////////////////////////////////////////////////
+    OptParser    opt;
+    bool         parsed;
+    string       inFilename, outFilename;
+    unsigned int dir = FFT::Forward;
+    
+    opt.addOption("o", "output", OptParser::OptType::value  , true,
+                  "output file name (default: result not saved)", "");
+    opt.addOption("b", "backward", OptParser::OptType::trigger, true,
+                  "backward Fourier transform (forward by default)");
+    opt.addOption("" , "help"  , OptParser::OptType::trigger, true,
+                  "show this help message and exit");
+    parsed = opt.parse(argc, argv);
+    if (!parsed or (opt.getArgs().size() != 1) or opt.gotOption("help"))
+    {
+        cerr << "usage: " << argv[0];
+        cerr << " <options> <input file>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
+        
+        return EXIT_FAILURE;
+    }
+    inFilename  = opt.getArgs()[0];
+    outFilename = opt.optionValue("o");
+    if (opt.gotOption("b"))
+    {
+        dir = FFT::Backward;
+    }
+    
+    // Fourier transform ///////////////////////////////////////////////////////
+    DMatSample in          = Io::load<DMatSample>(inFilename);
+    Index      nSample     = in.size(), l = in[central].rows();
+    bool       isInComplex = (in[central].cols() > 1);
+    CMatSample res(nSample, l, 1);
+    DMatSample out(nSample, l, 2);
+    GslFFT     ft(l);
+    
+    cout << "-- computing Fourier transform..." << endl;
+    FOR_STAT_ARRAY(in, s)
+    {
+        res[s].real() = in[s].col(0);
+        if (isInComplex)
+        {
+            res[s].imag() = in[s].col(1);
+        }
+        else
+        {
+            res[s].imag() = DVec::Constant(l, 0.);
+        }
+        ft(res[s], dir);
+        out[s].col(0) = res[s].real();
+        out[s].col(1) = res[s].imag();
+    }
+    
+    // output /////////////////////////////////////////////////////////////////
+    cout << scientific;
+    cout << "central value:\n"      << out[central];
+    cout << endl;
+    cout << "standard deviation:\n" << out.variance().cwiseSqrt();
+    cout << endl;
+    if (!outFilename.empty())
+    {
+        Io::save<DMatSample>(out, outFilename);
+    }
+    
+    return EXIT_SUCCESS;
+}

utils/sample-plot.cpp

@@ -0,0 +1,51 @@
+/*
+ * sample-plot.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2016 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/Io.hpp>
+#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/Plot.hpp>
+
+using namespace std;
+using namespace Latan;
+using namespace Math;
+
+int main(int argc, char *argv[])
+{
+    string xFileName, yFileName;
+    
+    if (argc != 3)
+    {
+        cerr << "usage: " << argv[0] << " <x sample> <y sample>" << endl;
+        
+        return EXIT_FAILURE;
+    }
+    xFileName = argv[1];
+    yFileName = argv[2];
+    
+    Plot p;
+    DMatSample x, y;
+    
+    x = Io::load<DMatSample>(xFileName);
+    y = Io::load<DMatSample>(yFileName);
+    
+    p << PlotData(x, y);
+    p.display();
+    
+    return EXIT_SUCCESS;
+}