Merge branch 'release/3.3'

# Conflicts:
#	configure.ac

.gitignore

@@ -3,7 +3,7 @@ build/*
 build_debug/*
 sandbox/*
 
-# Apple stuffs
+# Apple stuff
 *.xcodeproj
 .DS_Store
 
@@ -30,3 +30,12 @@ lib/eigen_files.mk
 
 # CI builds
 ci-scripts/local/*
+
+# CLion stuff
+.idea/*
+CMakeLists.txt
+cmake-build-debug/*
+
+# VS Code Studio stuff
+.vscode
+*.code-workspace

.travis.yml

@@ -7,6 +7,7 @@ notifications:
 cache:
   directories:
     - ci-scripts/local
+    - ci-scripts/clang
 
 matrix:
   include:
@@ -44,37 +45,43 @@ matrix:
         apt:
           sources:
             - ubuntu-toolchain-r-test
-            - llvm-toolchain-precise-3.6
           packages:
-            - clang-3.6
+            - g++-4.8
             - libgsl0-dev
             - flex
             - bison
-      env: VERSION=-3.6
+      env: CLANG_LINK=http://llvm.org/releases/3.6.0/clang+llvm-3.6.0-x86_64-linux-gnu-ubuntu-14.04.tar.xz
     - compiler: clang
       addons:
         apt:
           sources:
             - ubuntu-toolchain-r-test
-            - llvm-toolchain-precise-3.7
           packages:
-            - clang-3.7
+            - g++-4.8
             - libgsl0-dev
             - flex
             - bison
-      env: VERSION=-3.7
+      env: CLANG_LINK=http://llvm.org/releases/3.7.0/clang+llvm-3.7.0-x86_64-linux-gnu-ubuntu-14.04.tar.xz
 
 before_install:
+  - export LATDIR=`pwd`
+  - cd ci-scripts
+  - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]] && [ ! -e clang/bin ]; then wget $CLANG_LINK; tar -xf `basename $CLANG_LINK`; mkdir clang; mv clang+*/* clang/; fi
+  - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export PATH="${LATDIR}/ci-scripts/clang/bin:${PATH}"; fi
+  - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CC" == "clang" ]]; then export LD_LIBRARY_PATH="${LATDIR}/ci-scripts/clang/lib:${LD_LIBRARY_PATH}"; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew update; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install gsl; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install bison flex; export PATH="/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:$PATH"; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]] && [[ "$CC" == "gcc" ]]; then brew install gcc5; fi
   
 install:
-  - export LATDIR=`pwd`
   - export CC=$CC$VERSION
   - export CXX=$CXX$VERSION
-  - cd ci-scripts
+  - echo $PATH
+  - which $CC
+  - $CC  --version
+  - which $CXX
+  - $CXX --version
   - ./install-deps.sh `pwd`/local
   - cd ..
   - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then export LD_LIBRARY_PATH=${LATDIR}/ci-scripts/local/lib:$LD_LIBRARY_PATH; fi

Makefile.am

@@ -1,3 +1,3 @@
-SUBDIRS = lib utils examples
+SUBDIRS = lib utils physics examples
 
 ACLOCAL_AMFLAGS = -I .buildutils/m4

Readme.md

@@ -55,11 +55,47 @@ 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`
+
+#### v3.2 (needs LatCore 1.1)
+Additions:
+* 2-pt function fitter `latan-2pt-fit`
+* Tool to extract one element of a matrix sample `latan-sample-element`
+* Band plotting
+
+Changes:
+* Sample utilities renamed `latan-sample-*`
+* Resample utility renamed `latan-resample`
+
 #### v3.1.2
 Fixes:
 * HDF5 archive URL update in build scripts
 
-#### v3.1.1
+#### v3.1.1 (needs LatCore 1.0)
 Fixes:
 * Minuit precision fixed
 * Minor fit interface fixes

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.1.1],[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])
@@ -131,7 +131,8 @@ AC_SUBST([LIBS])
 AC_SUBST([AM_CFLAGS])
 AC_SUBST([AM_LDFLAGS])
 
-AC_CONFIG_FILES([Makefile lib/Makefile utils/Makefile examples/Makefile])
+AC_CONFIG_FILES([Makefile lib/Makefile utils/Makefile physics/Makefile 
+                 examples/Makefile])
 AC_OUTPUT
 
 echo "*********************************************"

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));
-    
-    gsl_integration_qags(&gslF, xMin, xMax, 0.0, precision_, limit_, workspace_,
-                         &result, &error_);
+    if ((xMin > -Math::inf) and (xMax < Math::inf))
+    {
+        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;
-    const double e  = 2.7182818284590452353602874713526624977572;
+    constexpr double pi  = 3.1415926535897932384626433832795028841970;
+    constexpr double e   = 2.7182818284590452353602874713526624977572;
+    constexpr double inf = std::numeric_limits<double>::infinity();
 }
 
 

lib/MathInterpreter.cpp

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

lib/Plot.cpp

@@ -183,6 +183,20 @@ PlotHLine::PlotHLine(const double y)
     setCommand(strFrom(y));
 }
 
+// PlotHBand constructor ///////////////////////////////////////////////////////
+PlotBand::PlotBand(const double xMin, const double xMax, const double yMin,
+                   const double yMax, const double opacity)
+{
+    setCommand("'< printf \"%e %e\\n%e %e\\n%e %e\\n%e %e\\n%e %e\\n\" "
+               + strFrom(xMin) + " " + strFrom(yMin) + " "
+               + strFrom(xMax) + " " + strFrom(yMin) + " "
+               + strFrom(xMax) + " " + strFrom(yMax) + " "
+               + strFrom(xMin) + " " + strFrom(yMax) + " "
+               + strFrom(xMin) + " " + strFrom(yMin)
+               + "' u 1:2 w filledcurves closed fs solid " + strFrom(opacity)
+               + " noborder");
+}
+
 // PlotFunction constructor ////////////////////////////////////////////////////
 PlotFunction::PlotFunction(const DoubleFunction &function, const double xMin,
                            const double xMax, const unsigned int nPoint)

lib/Plot.hpp

@@ -106,6 +106,16 @@ public:
     virtual ~PlotHLine(void) = default;
 };
 
+class PlotBand: public PlotObject
+{
+public:
+    // constructor
+    PlotBand(const double xMin, const double xMax, const double yMin,
+             const double yMax, const double opacity = 0.15);
+    // destructor
+    virtual ~PlotBand(void) = default;
+};
+
 class PlotFunction: public PlotObject
 {
 public:

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)
-//{
-//    setData(data.x(i), data.y(j));
-//}
+void TabFunction::setInterpolationType(const InterpType interpType)
+{
+    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
 

lib/XYSampleData.cpp

@@ -98,7 +98,7 @@ void SampleFitResult::print(const bool printXsi, ostream &out) const
     Index pMax = printXsi ? size() : nPar_;
     DMat  err = this->variance().cwiseSqrt();
     
-    sprintf(buf, "chi^2/dof= %.1f/%d= %.2f -- p-value= %.2e", getChi2(),
+    sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- p-value= %.2e", getChi2(),
             static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
     out << buf << endl;
     for (Index p = 0; p < pMax; ++p)
@@ -171,6 +171,20 @@ const DSample & XYSampleData::y(const Index k, const Index j) const
     return yData_[j].at(k);
 }
 
+void XYSampleData::setUnidimData(const DMatSample &xData,
+                                 const vector<const DMatSample *> &v)
+{
+    FOR_STAT_ARRAY(xData, s)
+    FOR_VEC(xData[central], r)
+    {
+        x(r, 0)[s] = xData[s](r);
+        for (unsigned int j = 0; j < v.size(); ++j)
+        {
+            y(r, j)[s] = (*(v[j]))[s](r);
+        }
+    }
+}
+
 const DMat & XYSampleData::getXXVar(const Index i1, const Index i2)
 {
     checkXDim(i1);
@@ -253,6 +267,7 @@ void XYSampleData::setDataToSample(const Index s)
 const XYStatData & XYSampleData::getData(void)
 {
     setDataToSample(central);
+    computeVarMat();
     
     return data_;
 }
@@ -277,12 +292,12 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
         if (s == central)
         {
             sampleResult = data_.fit(minimizer, initCopy, v);
+            initCopy     = sampleResult.segment(0, initCopy.size());
         }
         else
         {
             sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
         }
-        initCopy        = sampleResult.segment(0, initCopy.size());
         result[s]       = sampleResult;
         result.chi2_[s] = sampleResult.getChi2();
         for (unsigned int j = 0; j < v.size(); ++j)

lib/XYSampleData.hpp

@@ -80,6 +80,8 @@ public:
     const DMatSample & x(const Index k);
     DSample &          y(const Index k, const Index j);
     const DSample &    y(const Index k, const Index j) const;
+    void               setUnidimData(const DMatSample &xData,
+                                     const std::vector<const DMatSample *> &v);
     template <typename... Ts>
     void               setUnidimData(const DMatSample &xData,
                                      const Ts & ...yDatas);
@@ -141,17 +143,9 @@ void XYSampleData::setUnidimData(const DMatSample &xData, const Ts & ...yDatas)
     static_assert(static_or<std::is_assignable<DMatSample, Ts>::value...>::value,
                   "y data arguments are not compatible with DMatSample");
     
-    std::vector<const DMatSample *> yData{&yDatas...};
+    std::vector<const DMatSample *> v{&yDatas...};
     
-    FOR_STAT_ARRAY(xData, s)
-    FOR_VEC(xData[central], r)
-    {
-        x(r, 0)[s] = xData[s](r);
-        for (unsigned int j = 0; j < yData.size(); ++j)
-        {
-            y(r, j)[s] = (*(yData[j]))[s](r);
-        }
-    }
+    setUnidimData(xData, v);
 }
 
 template <typename... Ts>

lib/XYStatData.cpp

@@ -66,7 +66,7 @@ void FitResult::print(const bool printXsi, ostream &out) const
     char  buf[256];
     Index pMax = printXsi ? size() : nPar_;
     
-    sprintf(buf, "chi^2/dof= %.1f/%d= %.2f -- p-value= %.2e", getChi2(),
+    sprintf(buf, "chi^2/dof= %.1e/%d= %.2e -- p-value= %.2e", getChi2(),
             static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
     out << buf << endl;
     for (Index p = 0; p < pMax; ++p)

physics/2pt-fit.cpp

@@ -0,0 +1,332 @@
+#include <LatCore/OptParser.hpp>
+#include <LatAnalyze/CompiledModel.hpp>
+#include <LatAnalyze/Io.hpp>
+#include <LatAnalyze/MatSample.hpp>
+#include <LatAnalyze/Math.hpp>
+#include <LatAnalyze/MinuitMinimizer.hpp>
+#include <LatAnalyze/NloptMinimizer.hpp>
+#include <LatAnalyze/Plot.hpp>
+#include <LatAnalyze/XYSampleData.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char *argv[])
+{
+    // parse arguments /////////////////////////////////////////////////////////
+    OptParser            opt;
+    bool                 parsed, doPlot, doHeatmap, doCorr, fold;
+    string               corrFileName, model, outFileName, outFmt;
+    Index                ti, tf, shift, nPar, thinning;
+    double               svdTol;
+    Minimizer::Verbosity verbosity;
+    
+    opt.addOption("" , "ti"       , OptParser::OptType::value  , false,
+                  "initial fit time");
+    opt.addOption("" , "tf"       , OptParser::OptType::value  , false,
+                  "final fit time");
+    opt.addOption("t" , "thinning", OptParser::OptType::value  , true,
+                  "thinning of the time interval", "1");
+    opt.addOption("s", "shift"    , OptParser::OptType::value  , true,
+                  "time variable shift", "0");
+    opt.addOption("m", "model"    , OptParser::OptType::value  , true,
+                  "fit model (exp|exp2|exp3|cosh|cosh2|cosh3|<interpreter code>)", "cosh");
+    opt.addOption("" , "nPar"     , OptParser::OptType::value  , true,
+                  "number of model parameters for custom models "
+                  "(-1 if irrelevant)", "-1");
+    opt.addOption("" , "svd"      , OptParser::OptType::value  , true,
+                  "singular value elimination threshold", "0.");
+    opt.addOption("v", "verbosity", OptParser::OptType::value  , true,
+                  "minimizer verbosity level (0|1|2)", "0");
+    opt.addOption("o", "output", OptParser::OptType::value  , true,
+                  "output file", "");
+    opt.addOption("" , "uncorr"   , OptParser::OptType::trigger, true,
+                  "only do the uncorrelated fit");
+    opt.addOption("" , "fold"   , OptParser::OptType::trigger, true,
+                  "fold the correlator");
+    opt.addOption("p", "plot"     , OptParser::OptType::trigger, true,
+                  "show the fit plot");
+    opt.addOption("h", "heatmap"  , OptParser::OptType::trigger, true,
+                  "show the fit correlation heatmap");
+    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] << " <options> <correlator file>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
+        
+        return EXIT_FAILURE;
+    }
+    corrFileName = opt.getArgs().front();
+    ti           = opt.optionValue<Index>("ti");
+    tf           = opt.optionValue<Index>("tf");
+    thinning     = opt.optionValue<Index>("t");
+    shift        = opt.optionValue<Index>("s");
+    model        = opt.optionValue("m");
+    nPar         = opt.optionValue<Index>("nPar");
+    svdTol       = opt.optionValue<double>("svd");
+    outFileName  = opt.optionValue<string>("o");
+    doCorr       = !opt.gotOption("uncorr");
+    fold         = opt.gotOption("fold");
+    doPlot       = opt.gotOption("p");
+    doHeatmap    = opt.gotOption("h");
+    switch (opt.optionValue<unsigned int>("v"))
+    {
+        case 0:
+            verbosity = Minimizer::Verbosity::Silent;
+            break;
+        case 1:
+            verbosity = Minimizer::Verbosity::Normal;
+            break;
+        case 2:
+            verbosity = Minimizer::Verbosity::Debug;
+            break;
+        default:
+            cerr << "error: wrong verbosity level" << endl;
+            return EXIT_FAILURE;
+    }
+    
+    // load correlator /////////////////////////////////////////////////////////
+    DMatSample tmp, corr;
+    Index      nSample, nt;
+    
+    tmp     = Io::load<DMatSample>(corrFileName);
+    nSample = tmp.size();
+    nt      = tmp[central].rows();
+    tmp     = tmp.block(0, 0, nt, 1);
+    corr    = tmp;
+    FOR_STAT_ARRAY(corr, s)
+    {
+        for (Index t = 0; t < nt; ++t)
+        {
+            corr[s]((t - shift + nt)%nt) = tmp[s](t);
+        }
+    }
+    if (fold)
+    {
+        tmp = corr;
+        FOR_STAT_ARRAY(corr, s)
+        {
+            for (Index t = 0; t < nt; ++t)
+            {
+                corr[s](t) = 0.5*(tmp[s](t) + tmp[s]((nt - t) % nt));
+            }
+        }
+    }
+    
+    // make model //////////////////////////////////////////////////////////////
+    DoubleModel mod;
+    bool        coshModel = false;
+    
+    if ((model == "exp") or (model == "exp1"))
+    {
+        nPar = 2;
+        mod.setFunction([](const double *x, const double *p)
+                        {
+                            return p[1]*exp(-p[0]*x[0]);
+                        }, 1, nPar);
+    }
+    else if (model == "exp2")
+    {
+        nPar = 4;
+        mod.setFunction([](const double *x, const double *p)
+                        {
+                            return p[1]*exp(-p[0]*x[0]) + p[3]*exp(-p[2]*x[0]);
+                        }, 1, nPar);
+    }
+    else if (model == "exp3")
+    {
+        nPar = 6;
+        mod.setFunction([](const double *x, const double *p)
+                        {
+                            return p[1]*exp(-p[0]*x[0]) + p[3]*exp(-p[2]*x[0])
+                                   + p[5]*exp(-p[4]*x[0]);
+                        }, 1, nPar);
+    }
+    else if ((model == "cosh") or (model == "cosh1"))
+    {
+        coshModel = true;
+        nPar      = 2;
+        mod.setFunction([nt](const double *x, const double *p)
+                        {
+                            return p[1]*(exp(-p[0]*x[0])+exp(-p[0]*(nt-x[0])));
+                        }, 1, nPar);
+    }
+    else if (model == "cosh2")
+    {
+        coshModel = true;
+        nPar      = 4;
+        mod.setFunction([nt](const double *x, const double *p)
+                        {
+                            return p[1]*(exp(-p[0]*x[0])+exp(-p[0]*(nt-x[0])))
+                                 + p[3]*(exp(-p[2]*x[0])+exp(-p[2]*(nt-x[0])));
+                        }, 1, nPar);
+    }
+    else if (model == "cosh3")
+    {
+        coshModel = true;
+        nPar      = 6;
+        mod.setFunction([nt](const double *x, const double *p)
+                        {
+                            return p[1]*(exp(-p[0]*x[0])+exp(-p[0]*(nt-x[0])))
+                                 + p[3]*(exp(-p[2]*x[0])+exp(-p[2]*(nt-x[0])))
+                                 + p[5]*(exp(-p[2]*x[0])+exp(-p[4]*(nt-x[0])));
+                        }, 1, nPar);
+    }
+    else
+    {
+        if (nPar > 0)
+        {
+            mod = compile(model, 1, nPar);
+        }
+        else
+        {
+            cerr << "error: please specify the number of model parameter"
+                    " using the --nPar function" << endl;
+            
+            return EXIT_FAILURE;
+        }
+    }
+    
+    // fit /////////////////////////////////////////////////////////////////////
+    DMatSample          tvec(nSample);
+    XYSampleData        data(nSample);
+    SampleFitResult     fit;
+    DVec                init(nPar);
+    NloptMinimizer      globMin(NloptMinimizer::Algorithm::GN_CRS2_LM);
+    MinuitMinimizer     locMin;
+    vector<Minimizer *> unCorrMin{&globMin, &locMin};
+
+    FOR_STAT_ARRAY(tvec, s)
+    {
+        tvec[s] = DVec::LinSpaced(nt, 0, nt - 1);
+    }
+    data.addXDim(nt, "t/a", true);
+    data.addYDim("C(t)");
+    data.setUnidimData(tvec, corr);
+    for (Index p = 0; p < nPar; p += 2)
+    {
+        mod.parName().setName(p, "E_" + strFrom(p/2));
+        mod.parName().setName(p + 1, "Z_" + strFrom(p/2));
+    }
+    init(0) = log(data.y(nt/4, 0)[central]/data.y(nt/4 + 1, 0)[central]);
+    init(1) = data.y(nt/4, 0)[central]/(exp(-init(0)*nt/4));
+    for (Index p = 2; p < nPar; p += 2)
+    {
+        init(p)     = 2*init(p - 2);
+        init(p + 1) = init(p - 1)/2.;
+    }
+    for (Index p = 0; p < nPar; p += 2)
+    {
+        globMin.setLowLimit(p, 0.);
+        globMin.setHighLimit(p, 10.*init(p));
+        globMin.setLowLimit(p + 1, -10.*init(p + 1));
+        globMin.setHighLimit(p + 1, 10.*init(p + 1));
+        locMin.setLowLimit(p, 0.);
+    }
+    globMin.setPrecision(0.001);
+    globMin.setMaxIteration(100000);
+    globMin.setVerbosity(verbosity);
+    locMin.setMaxIteration(1000000);
+    locMin.setVerbosity(verbosity);
+    for (Index t = 0; t < nt; ++t)
+    {
+        data.fitPoint((t >= ti) and (t <= tf)
+                      and ((t - ti) % thinning == 0), t);
+    }
+    if (doCorr)
+    {
+        cout << "-- uncorrelated fit..." << endl;
+    }
+    cout << "using model '" << model << "'" << endl;
+    data.setSvdTolerance(svdTol);
+    data.assumeYYCorrelated(false, 0, 0);
+    fit = data.fit(unCorrMin, init, mod);
+    fit.print();
+    if (doCorr)
+    {
+        cout << "-- correlated fit..." << endl;
+        cout << "using model '" << model << "'" << endl;
+        init = fit[central];
+        data.assumeYYCorrelated(true, 0, 0);
+        fit = data.fit(locMin, init, mod);
+        fit.print();
+    }
+    
+    // plots ///////////////////////////////////////////////////////////////////
+    if (doPlot)
+    {
+        Plot       p;
+        DMatSample effMass(nSample);
+        DVec       effMassT, fitErr;
+        Index      maxT = (coshModel) ? (nt - 2) : (nt - 1);
+        double     e0, e0Err;
+        
+        p << PlotRange(Axis::x, 0, nt - 1);
+        p << LogScale(Axis::y);
+        p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1);
+        p << Color("rgb 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1);
+        p << Color("rgb 'red'") << PlotData(data.getData());
+        p.display();
+        effMass.resizeMat(maxT, 1);
+        effMassT.setLinSpaced(maxT, 1, maxT);
+        fitErr = fit.variance().cwiseSqrt();
+        e0     = fit[central](0);
+        e0Err  = fitErr(0);
+        if (coshModel)
+        {
+            FOR_STAT_ARRAY(effMass, s)
+            {
+                for (Index t = 1; t < nt - 1; ++t)
+                {
+                    effMass[s](t - 1) = acosh((corr[s](t-1) + corr[s](t+1))
+                                              /(2.*corr[s](t)));
+                }
+            }
+        }
+        else
+        {
+            FOR_STAT_ARRAY(effMass, s)
+            {
+                for (Index t = 1; t < nt; ++t)
+                {
+                    effMass[s](t - 1) = log(corr[s](t-1)/corr[s](t));
+                }
+            }
+        }
+        p.reset();
+        p << PlotRange(Axis::x, 1, maxT);
+        p << PlotRange(Axis::y, e0 - 20.*e0Err, e0 + 20.*e0Err);
+        p << Color("rgb 'blue'") << PlotBand(0, maxT, e0 - e0Err, e0 + e0Err);
+        p << Color("rgb 'blue'") << PlotHLine(e0);
+        p << Color("rgb 'red'") << PlotData(effMassT, effMass);
+        p.display();
+        p.save("test");
+    }
+    if (doHeatmap)
+    {
+        Plot  p;
+        Index n  = data.getFitVarMat().rows();
+        DMat  id = DMat::Identity(n, n);
+        
+        p << PlotMatrix(Math::varToCorr(data.getFitVarMat()));
+        p << Caption("correlation matrix");
+        p.display();
+        if (svdTol > 0.)
+        {
+            p.reset();
+            p << PlotMatrix(id - data.getFitVarMat()*data.getFitVarMatPInv());
+            p << Caption("singular space projector");
+            p.display();
+        }
+    }
+    
+    // output //////////////////////////////////////////////////////////////////
+    if (!outFileName.empty())
+    {
+        Io::save(fit, outFileName);
+    }
+
+    return EXIT_SUCCESS;
+}

physics/LatAnalyze

@@ -0,0 +1 @@
+../lib

physics/Makefile.am

@@ -0,0 +1,15 @@
+if CXX_GNU
+    COM_CXXFLAGS = -Wall -W -pedantic -Wno-deprecated-declarations
+else
+if CXX_INTEL
+    COM_CXXFLAGS = -wd1682 -Wall
+endif
+endif
+
+bin_PROGRAMS = latan-2pt-fit
+
+latan_2pt_fit_SOURCES  = 2pt-fit.cpp
+latan_2pt_fit_CXXFLAGS = $(COM_CXXFLAGS)
+latan_2pt_fit_LDFLAGS  = -L../lib/.libs -lLatAnalyze
+
+ACLOCAL_AMFLAGS = -I .buildutils/m4

utils/Makefile.am

@@ -7,25 +7,40 @@ endif
 endif
 
 bin_PROGRAMS =            \
-    latan_make_fake_sample\
-    latan_sample_combine  \
-    latan_sample_plot_corr\
-    latan_sample_read     \
-    latan_resample
+    latan-sample-combine  \
+    latan-sample-element  \
+    latan-sample-fake     \
+    latan-sample-ft       \
+    latan-sample-plot     \
+    latan-sample-plot-corr\
+    latan-sample-read     \
+    latan-resample
 
-latan_make_fake_sample_SOURCES  = make_fake_sample.cpp
-latan_make_fake_sample_CXXFLAGS = $(COM_CXXFLAGS)
-latan_make_fake_sample_LDFLAGS  = -L../lib/.libs -lLatAnalyze
-
-latan_sample_combine_SOURCES  = sample_combine.cpp
+latan_sample_combine_SOURCES  = sample-combine.cpp
 latan_sample_combine_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_combine_LDFLAGS  = -L../lib/.libs -lLatAnalyze
 
-latan_sample_plot_corr_SOURCES  = sample_plot_corr.cpp
+latan_sample_element_SOURCES  = sample-element.cpp
+latan_sample_element_CXXFLAGS = $(COM_CXXFLAGS)
+latan_sample_element_LDFLAGS  = -L../lib/.libs -lLatAnalyze
+
+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_read_SOURCES  = sample_read.cpp
+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

@@ -17,74 +17,69 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <LatCore/OptParser.hpp>
 #include <LatAnalyze/Dataset.hpp>
 #include <LatAnalyze/Io.hpp>
+#include <LatAnalyze/includes.hpp>
 
 #ifndef DEF_NSAMPLE
-#define DEF_NSAMPLE 100
+#define DEF_NSAMPLE "100"
+#endif
+
+#ifdef HAVE_HDF5
+#define DEF_FMT "h5"
+#else
+#define DEF_FMT "sample"
 #endif
 
 using namespace std;
 using namespace Latan;
 
-static void usage(const string &cmdName)
-{
-    cerr << "usage: " << cmdName;
-    cerr << " [-n <nsample> -b <bin size> -r <seed> -o <output dir> -f {h5|sample}]";
-    cerr << " <data list> <name list>";
-    cerr << endl;
-    exit(EXIT_FAILURE);
-}
-
 int main(int argc, char *argv[])
 {
     // argument parsing ////////////////////////////////////////////////////////
-    int           c;
+    OptParser     opt;
+    bool          parsed, dumpBoot;
     random_device rd;
     SeedType      seed = rd();
-    string        manFileName, nameFileName, cmdName, outDirName = ".";
-    string        ext = "h5";
-    Index         binSize = 1, nSample = DEF_NSAMPLE;
+    string        manFileName, nameFileName, outDirName;
+    string        ext;
+    Index         binSize, nSample;
     
-    opterr = 0;
-    cmdName = basename(argv[0]);
-    while ((c = getopt(argc, argv, "b:n:r:o:f:")) != -1)
+    opt.addOption("n", "nsample"   , OptParser::OptType::value,   true,
+                  "number of samples", DEF_NSAMPLE);
+    opt.addOption("b", "bin"       , OptParser::OptType::value,   true,
+                  "bin size", "1");
+    opt.addOption("r", "seed"      , OptParser::OptType::value,   true,
+                  "random generator seed (default: random)");
+    opt.addOption("o", "output-dir", OptParser::OptType::value,   true,
+                  "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);
+    if (!parsed or (opt.getArgs().size() != 2) or opt.gotOption("help"))
     {
-        switch (c)
-        {
-            case 'b':
-                binSize = strTo<Index>(optarg);
-                break;
-            case 'n':
-                nSample = strTo<Index>(optarg);
-                break;
-            case 'o':
-                outDirName = optarg;
-                break;
-            case 'r':
-                seed = strTo<SeedType>(optarg);
-                break;
-            case 'f':
-                ext = optarg;
-                break;
-            case '?':
-                cerr << "error parsing option -" << char(optopt) << endl;
-                usage(cmdName);
-                break;
-            default:
-                usage(cmdName);
-                break;
-        }
+        cerr << "usage: " << argv[0];
+        cerr << " <datafile list> <name list> <options>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
+        
+        return EXIT_FAILURE;
     }
-    if (argc - optind == 2)
+    nSample = opt.optionValue<Index>("n");
+    binSize = opt.optionValue<Index>("b");
+    if (opt.gotOption("r"))
     {
-        manFileName  = argv[optind];
-        nameFileName = argv[optind+1];
-    }
-    else
-    {
-        usage(cmdName);
+        seed = opt.optionValue<SeedType>("r");
     }
+    ext          = opt.optionValue("f");
+    outDirName   = opt.optionValue("o");
+    dumpBoot     = opt.gotOption("d");
+    manFileName  = opt.getArgs()[0];
+    nameFileName = opt.getArgs()[1];
     
     // parameter parsing ///////////////////////////////////////////////////////
     vector<string> dataFileName, name;
@@ -132,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

@@ -17,31 +17,29 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <LatCore/OptParser.hpp>
 #include <LatAnalyze/Io.hpp>
 #include <LatAnalyze/CompiledFunction.hpp>
 
 using namespace std;
 using namespace Latan;
 
-static void usage(const string &cmdName)
+template <typename T>
+static void loadAndCheck(vector<T> &sample __dumb,
+                         const vector<string> &fileName __dumb)
 {
-    cerr << "usage: " << cmdName;
-    cerr << " [-o <output sample>]";
-    cerr << " <n> <function> <sample 1> ... <sample n>";
-    cerr << endl;
-    exit(EXIT_FAILURE);
+    abort();
 }
 
-template <typename T>
-static void loadAndCheck(vector<T> &sample, const vector<string> &fileName)
+template <>
+void loadAndCheck(vector<DSample> &sample, const vector<string> &fileName)
 {
-    const unsigned int n = sample.size();
+    const unsigned int n = static_cast<unsigned int>(sample.size());
     Index              nSample = 0;
     
-    cout << "-- loading data..." << endl;
     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();
@@ -55,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)
@@ -66,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);
@@ -96,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]);
@@ -137,49 +191,41 @@ void process(const string &outFileName, const vector<string> &fileName,
 int main(int argc, char *argv[])
 {
     // argument parsing ////////////////////////////////////////////////////////
+    OptParser      opt;
+    bool           parsed;
     string         cmdName, outFileName = "", code;
     vector<string> fileName;
-    int            c;
     unsigned int   n = 0;
 
-    opterr = 0;
-    cmdName = basename(argv[0]);
-    while ((c = getopt(argc, argv, "o:")) != -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 (opt.getArgs().size() >= 1)
     {
-        switch (c)
-        {
-            case 'o':
-                outFileName = optarg;
-                break;
-            case '?':
-                cerr << "error parsing option -" << char(optopt) << endl;
-                usage(cmdName);
-                break;
-            default:
-                usage(cmdName);
-                break;
-        }
-    }
-    if (argc - optind >= 1)
-    {
-        n = strTo<unsigned int>(argv[optind]);
-        if (argc - optind == static_cast<int>(n + 2))
-        {
-            fileName.resize(n);
-            code = argv[optind + 1];
-            for (unsigned int i = 0; i < n; ++i)
-            {
-                fileName[i] = argv[optind + 2 + i];
-            }
-        }
-        else
-        {
-            usage(cmdName);
-        }
+        n = strTo<unsigned int>(opt.getArgs()[0]);
     }
     else
     {
-        usage(cmdName);
+        parsed = false;
+    }
+    if (!parsed or (opt.getArgs().size() != n + 2) or opt.gotOption("help"))
+    {
+        cerr << "usage: " << argv[0];
+        cerr << " <options> <n> <function> <sample 1> ... <sample n>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
+        
+        return EXIT_FAILURE;
+    }
+    if (opt.gotOption("o"))
+    {
+        outFileName = opt.optionValue("o");
+    }
+    code = opt.getArgs()[1];
+    for (unsigned int i = 0; i < n; ++i)
+    {
+        fileName.push_back(opt.getArgs()[2 + i]);
     }
     
     // process data ////////////////////////////////////////////////////////////

utils/sample-element.cpp

@@ -0,0 +1,103 @@
+/*
+ * 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>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char* argv[])
+{
+    // argument parsing ////////////////////////////////////////////////////////
+    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");
+    
+    // Data extraction /////////////////////////////////////////////////////////
+    auto inputData = Io::load<DMatSample>(inFilename);
+    
+    if ((nr == 1) and (nc == 1))
+    {
+        DSample outputData(inputData.size());
+        
+        FOR_STAT_ARRAY(inputData, s)
+        {
+            outputData[s] = inputData[s](r, c);
+        }
+        
+        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-corr.cpp

@@ -42,6 +42,7 @@ int main(int argc, char *argv[])
     cout << "-- computing variance matrix from '" << fileName << "'..." << endl;
     name   = Io::getFirstName(fileName);
     sample = Io::load<DMatSample>(fileName);
+    sample = sample.block(0, 0, sample[central].rows(), 1);
     var    = sample.varianceMatrix();
     p << PlotMatrix(varToCorr(var));
     p.display();

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;
+}