Merge branch 'develop' into feature/fit-phys

# Conflicts:
#	configure.ac

.gitignore

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

.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,6 +55,29 @@ 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.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 (needs LatCore 1.0)
+Fixes:
+* Minuit precision fixed
+* Minor fit interface fixes
+
 #### v3.1
 Additions:
 * Wrappers to NLopt and GSL minimisers.

ci-scripts/install-hdf5.sh

@@ -11,7 +11,7 @@ PREFIX=$1
 set -ex
 INITDIR=`pwd`
 cd local/build
-wget http://www.hdfgroup.org/ftp/HDF5/current/src/${NAME}.tar.gz
+wget http://www.hdfgroup.org/ftp/HDF5/releases/${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],[antonin.portelli@me.com],[LatAnalyze])
+AC_INIT([LatAnalyze],[3.2.1-dev],[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/AsciiFile.cpp

@@ -170,7 +170,7 @@ void AsciiFile::open(const string &name, const unsigned int mode)
 
 std::string AsciiFile::load(const string &name)
 {
-    if ((mode_ & Mode::read)&&(isOpen()))
+    if ((mode_ & Mode::read) and (isOpen()))
     {
         if (!isParsed_)
         {

lib/File.cpp

@@ -57,7 +57,7 @@ void File::deleteData(void)
 
 void File::checkWritability(void)
 {
-    if (!((mode_ & Mode::write)||(mode_ & Mode::append))||!isOpen())
+    if (!((mode_ & Mode::write) or (mode_ & Mode::append)) or !isOpen())
     {
         LATAN_ERROR(Io, "file '" + name_ + "' is not writable");
     }

lib/File.hpp

@@ -110,6 +110,12 @@ const IoT& File::getData(const std::string &name) const
         LATAN_ERROR(Definition, "no data with name '" + name + "' in file "
                     + name_);
     }
+    catch(std::bad_cast)
+    {
+        LATAN_ERROR(Definition, "data with name '" + name + "' in file "
+                    + name_ + " does not have type '" + typeid(IoT).name()
+                    + "'");
+    }
 }
 
 END_LATAN_NAMESPACE

lib/FitInterface.cpp

@@ -58,7 +58,7 @@ void FitInterface::addXDim(const Index nData, const string name,
         scheduleDataCoordInit();
         if (!name.empty())
         {
-            xName().setName(getNXDim(), name);
+            xName().setName(getNXDim() - 1, name);
         }
     }
 }
@@ -70,7 +70,7 @@ void FitInterface::addYDim(const string name)
     scheduleLayoutInit();
     if (!name.empty())
     {
-        yName().setName(getNYDim(), name);
+        yName().setName(getNYDim() - 1, name);
     }
 }
 

lib/Function.cpp

@@ -183,7 +183,7 @@ DoubleFunction & DoubleFunction::operator op##=(const DoubleFunction &f)\
     setFunction(res, getNArg());\
     return *this;\
 }\
-DoubleFunction & DoubleFunction::operator op##=(const DoubleFunction &&f)\
+DoubleFunction & DoubleFunction::operator op##=(const DoubleFunction  and f)\
 {\
     *this op##= f;\
     return *this;\

lib/GslHybridRootFinder.cpp

@@ -122,7 +122,7 @@ GslHybridRootFinder::operator()(const vector<DoubleFunction *> &func)
             break;
         }
         status = gsl_multiroot_test_residual(solver_->f, getPrecision());
-    } while ((status == GSL_CONTINUE)&&(iter < getMaxIteration()));
+    } while ((status == GSL_CONTINUE) and (iter < getMaxIteration()));
     if (verbosity >= Verbosity::Debug)
     {
         cout << "--- done" << endl;

lib/Hdf5File.cpp

@@ -207,14 +207,14 @@ string Hdf5File::getFirstGroupName(void)
 {
     string res;
 
-    if ((mode_ & Mode::read)&&(isOpen()))
+    if ((mode_ & Mode::read) and (isOpen()))
     {
         auto firstGroupName = [](hid_t loc_id, const char *name, void *fname)
         {
             H5G_stat_t statbuf;
 
             H5Gget_objinfo(loc_id, name, 0, &statbuf);
-            if ((statbuf.type == H5G_GROUP) && (strlen((char *)fname) == 0))
+            if ((statbuf.type == H5G_GROUP) and (strlen((char *)fname) == 0))
             {
                 strncpy((char *)fname, name, maxGroupNameSize);
             }
@@ -268,7 +268,7 @@ void Hdf5File::load(DSample &ds, const DataSet &d)
 
 string Hdf5File::load(const string &name)
 {
-    if ((mode_ & Mode::read)&&(isOpen()))
+    if ((mode_ & Mode::read) and (isOpen()))
     {
         string           groupName;
         Group            group;

lib/Histogram.cpp

@@ -165,7 +165,7 @@ double Histogram::operator()(const double x) const
 // percentiles & confidence interval ///////////////////////////////////////////
 double Histogram::percentile(const double p) const
 {
-    if ((p < 0.0)||(p > 100.0))
+    if ((p < 0.0) or (p > 100.0))
     {
         LATAN_ERROR(Range, "percentile (" + strFrom(p) + ")"
                     " is outside the [0, 100] range");

lib/Io.cpp

@@ -38,7 +38,7 @@ unique_ptr<File> Io::open(const std::string &fileName, const unsigned int mode)
 {
     string ext = extension(fileName);
     
-    if ((ext == "dat")||(ext == "sample")||(ext == "seed"))
+    if ((ext == "dat") or (ext == "sample") or (ext == "seed"))
     {
         return unique_ptr<File>(new AsciiFile(fileName, mode));
     }

lib/MathInterpreter.cpp

@@ -478,7 +478,7 @@ void SemicolonNode::compile(Program &program, RunContext &context) const
         bool isSemiColumn = isDerivedFrom<SemicolonNode>(&n[i]);
         bool isKeyword    = isDerivedFrom<KeywordNode>(&n[i]);
         
-        if (isAssign||isSemiColumn||isKeyword)
+        if (isAssign or isSemiColumn or isKeyword)
         {
             n[i].compile(program, context);
         }
@@ -513,7 +513,7 @@ void AssignNode::compile(Program &program, RunContext &context) const
 }
 
 // MathOpNode compile //////////////////////////////////////////////////////////
-#define IFNODE(name, nArg) if ((n.getName() == (name))&&(n.getNArg() == nArg))
+#define IFNODE(name, nArg) if ((n.getName() == (name)) and (n.getNArg() == nArg))
 #define ELIFNODE(name, nArg) else IFNODE(name, nArg)
 #define ELSE else
 
@@ -690,6 +690,7 @@ void MathInterpreter::compile(RunContext &context)
         if (root_)
         {
             context.addVariable("pi", Math::pi);
+            context.addVariable("inf", HUGE_VAL);
             ADD_STDMATH_FUNCS(context);
             root_->compile(program_, context);
             for (unsigned int i = 0; i < program_.size(); ++i)
@@ -703,7 +704,7 @@ void MathInterpreter::compile(RunContext &context)
                 }
             }
         }
-        if (!root_||!gotReturn)
+        if (!root_ or !gotReturn)
         {
             LATAN_ERROR(Syntax, "expected 'return' in program '" + codeName_
                         + "'");

lib/MinuitMinimizer.cpp

@@ -58,9 +58,10 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
     using namespace ROOT;
     using namespace Minuit2;
     
-    DVec                        &x = getState();
+    DVec           &x = getState();
-    int                         printLevel = 0;
+    int            printLevel = 0;
-    EMinimizerType              minuitAlg  = kCombined;
+    EMinimizerType minuitAlg  = kCombined;
+    double         prec = getPrecision();
     
     // convert Latan parameters to Minuit parameters
     switch (getVerbosity())
@@ -75,16 +76,21 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
             printLevel = 3;
             break;
     }
+    // The factor of 0.002 here is to compensate the dirty hack in Minuit
+    // source used to match the C++ and F77 versions
+    // (cf. VariableMetricBuilder.cxx)
     switch (getAlgorithm())
     {
         case Algorithm::migrad:
-            minuitAlg = kMigrad;
+            minuitAlg  = kMigrad;
+            prec      /= 0.002;
             break;
         case Algorithm::simplex:
             minuitAlg = kSimplex;
             break;
         case Algorithm::combined:
-            minuitAlg = kCombined;
+            minuitAlg  = kCombined;
+            prec      /= 0.002;
             break;
     }
     
@@ -99,7 +105,7 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
     
     min.SetStrategy(2);
     min.SetMaxFunctionCalls(getMaxIteration());
-    min.SetTolerance(getPrecision());
+    min.SetTolerance(prec);
     min.SetPrintLevel(printLevel);
     
     // set function and variables
@@ -146,7 +152,7 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
         min.Minimize();
         status = min.Status();
         n++;
-    } while (status and (n < getMaxPass()));
+    } while ((status >= 2) and (n < getMaxPass()));
     if (getVerbosity() >= Verbosity::Normal)
     {
         cout << "=================================================" << endl;
@@ -154,7 +160,8 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
     switch (status)
     {
         case 1:
-            LATAN_WARNING("invalid minimum: covariance matrix was made positive");
+            // covariance matrix was made positive, the minimum is still good
+            // it just means that Minuit error analysis is inaccurate
             break;
         case 2:
             LATAN_WARNING("invalid minimum: Hesse analysis is not valid");

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)
@@ -458,7 +472,7 @@ void Plot::getProgramPath(void)
         {
             for (i=0;path[i];)
             {
-                for (j=i;(path[j])&&(path[j]!=':');j++);
+                for (j=i;(path[j]) and (path[j]!=':');j++);
                 lg = j - i;
                 strncpy(buf,path + i,(size_t)(lg));
                 if (lg == 0)

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/StatArray.hpp

@@ -262,7 +262,7 @@ namespace ReducOp
     inline Mat<double> tensProd(const Mat<double>  &v1,
                                 const Mat<double>  &v2)
     {
-        if ((v1.cols() != 1)||(v2.cols() != 1))
+        if ((v1.cols() != 1) or (v2.cols() != 1))
         {
             LATAN_ERROR(Size,
                         "tensorial product is only valid with column vectors");

lib/TabFunction.cpp

@@ -65,7 +65,7 @@ double TabFunction::operator()(const double *arg) const
     double result = 0.0, x = arg[0];
 
 
-    if ((x < value_.begin()->first) || (x >= value_.rbegin()->first)) {
+    if ((x < value_.begin()->first) or (x >= value_.rbegin()->first)) {
         LATAN_ERROR(Range, "tabulated function variable out of range "
                                "(x= " + strFrom(x) + " not in ["
                            + strFrom(value_.begin()->first) + ", "

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)
+    setUnidimData(xData, v);
-    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);
-        }
-    }
 }
 
 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/LatAnalyze

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

physics/Makefile.am

@@ -0,0 +1,19 @@
+if CXX_GNU
+    COM_CXXFLAGS = -Wall -W -pedantic -Wno-deprecated-declarations
+else
+if CXX_INTEL
+    COM_CXXFLAGS = -wd1682 -Wall
+endif
+endif
+
+bin_PROGRAMS = latan-fit-2pt latan-fit-phys
+
+latan_fit_2pt_SOURCES  = fit-2pt.cpp
+latan_fit_2pt_CXXFLAGS = $(COM_CXXFLAGS)
+latan_fit_2pt_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-2pt.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/fit-phys-env.cpp

@@ -0,0 +1,509 @@
+#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

@@ -0,0 +1,82 @@
+#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

@@ -0,0 +1,77 @@
+#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

@@ -7,25 +7,30 @@ endif
 endif
 
 bin_PROGRAMS =            \
-    latan_make_fake_sample\
+    latan-sample-combine  \
-    latan_sample_combine  \
+    latan-sample-element  \
-    latan_sample_plot_corr\
+    latan-sample-fake     \
-    latan_sample_read     \
+    latan-sample-plot-corr\
-    latan_resample
+    latan-sample-read     \
+    latan-resample
 
-latan_make_fake_sample_SOURCES  = make_fake_sample.cpp
+latan_sample_combine_SOURCES  = sample-combine.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_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_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_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,66 @@
  * 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;
     random_device rd;
     SeedType      seed = rd();
-    string        manFileName, nameFileName, cmdName, outDirName = ".";
+    string        manFileName, nameFileName, outDirName;
-    string        ext = "h5";
+    string        ext;
-    Index         binSize = 1, nSample = DEF_NSAMPLE;
+    Index         binSize, nSample;
     
-    opterr = 0;
+    opt.addOption("n", "nsample"   , OptParser::OptType::value,   true,
-    cmdName = basename(argv[0]);
+                  "number of samples", DEF_NSAMPLE);
-    while ((c = getopt(argc, argv, "b:n:r:o:f:")) != -1)
+    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("" , "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)
+        cerr << "usage: " << argv[0];
-        {
+        cerr << " <datafile list> <name list> <options>" << endl;
-            case 'b':
+        cerr << endl << "Possible options:" << endl << opt << endl;
-                binSize = strTo<Index>(optarg);
+        
-                break;
+        return EXIT_FAILURE;
-            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;
-        }
     }
-    if (argc - optind == 2)
+    nSample = opt.optionValue<Index>("n");
+    binSize = opt.optionValue<Index>("b");
+    if (opt.gotOption("r"))
     {
-        manFileName  = argv[optind];
+        seed = opt.optionValue<SeedType>("r");
-        nameFileName = argv[optind+1];
-    }
-    else
-    {
-        usage(cmdName);
     }
+    ext          = opt.optionValue("f");
+    outDirName   = opt.optionValue("o");
+    manFileName  = opt.getArgs()[0];
+    nameFileName = opt.getArgs()[1];
     
     // parameter parsing ///////////////////////////////////////////////////////
     vector<string> dataFileName, name;

utils/sample-combine.cpp

@@ -17,28 +17,19 @@
  * 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)
-{
-    cerr << "usage: " << cmdName;
-    cerr << " [-o <output sample>]";
-    cerr << " <n> <function> <sample 1> ... <sample n>";
-    cerr << endl;
-    exit(EXIT_FAILURE);
-}
-
 template <typename T>
 static void loadAndCheck(vector<T> &sample, const vector<string> &fileName)
 {
     const unsigned int n = sample.size();
     Index              nSample = 0;
     
-    cout << "-- loading data..." << endl;
     for (unsigned int i = 0; i < n; ++i)
     {
         sample[i] = Io::load<T>(fileName[i]);
@@ -137,49 +128,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;
+    opt.addOption("o", "output", OptParser::OptType::value  , true,
-    cmdName = basename(argv[0]);
+                  "output file name (default: result not saved)");
-    while ((c = getopt(argc, argv, "o:")) != -1)
+    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)
+        n = strTo<unsigned int>(opt.getArgs()[0]);
-        {
-            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);
-        }
     }
     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 ////////////////////////////////////////////////////////////
@@ -187,7 +170,7 @@ int main(int argc, char *argv[])
     {
         process<DSample>(outFileName, fileName, code);
     }
-    catch (bad_cast &e)
+    catch (Exceptions::Definition)
     {
         process<DMatSample>(outFileName, fileName, code);
     }

utils/sample-element.cpp

@@ -0,0 +1,39 @@
+#include <iostream>
+#include <string>
+
+#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) {
+        cout << "Usage: " << argv[0] << " <input filename> <row> <column> ";
+        cout << "[output filename]" << endl;
+        return -1;
+    }
+
+    string inFileName = argv[1];
+    auto row = strTo<Index>(argv[2]);
+    auto col = strTo<Index>(argv[3]);
+    string outFileName = (argc == 5) ? argv[4] : "";
+
+    auto inputData = Io::load<DMatSample>(inFileName);
+    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);
+        }
+
+        Io::save(outputData, outFileName);
+    }
+}

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-read.cpp

@@ -24,7 +24,7 @@ using namespace Latan;
 
 int main(int argc, char *argv[])
 {
-    if ((argc < 2) || (argc > 3))
+    if ((argc < 2) or (argc > 3))
     {
         cerr << "usage: " << argv[0] << " <file> [<copy>]" << endl;
 
@@ -46,7 +46,7 @@ int main(int argc, char *argv[])
             Io::save(s, copy, File::Mode::write, name);
         }
     }
-    catch (bad_cast &e)
+    catch (Exceptions::Definition)
     {
         DSample s    = Io::load<DSample>(fileName);
         string  name = Io::getFirstName(fileName);