big update: first fit implementation

examples/Makefile.am

@@ -16,6 +16,7 @@ endif
 
 noinst_PROGRAMS =           \
     exCompiledDoubleFunction\
+    exFit                   \
     exMat                   \
     exMathInterpreter       \
     exMin                   \
@@ -26,6 +27,10 @@ exCompiledDoubleFunction_SOURCES = exCompiledDoubleFunction.cpp
 exCompiledDoubleFunction_CFLAGS  = -g -O2
 exCompiledDoubleFunction_LDFLAGS = -L../latan/.libs -llatan
 
+exFit_SOURCES                    = exFit.cpp
+exFit_CFLAGS                     = -g -O2
+exFit_LDFLAGS                    = -L../latan/.libs -llatan
+
 exMat_SOURCES                    = exMat.cpp
 exMat_CFLAGS                     = -g -O2
 exMat_LDFLAGS                    = -L../latan/.libs -llatan

examples/exFit.cpp

@@ -0,0 +1,51 @@
+#include <iostream>
+#include <cmath>
+#include <latan/MinuitMinimizer.hpp>
+#include <latan/Plot.hpp>
+#include <latan/RandGen.hpp>
+#include <latan/XYStatData.hpp>
+
+using namespace std;
+using namespace Latan;
+
+const Index  nPoint = 20;
+const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast<double>(nPoint);
+
+int main(void)
+{
+    // generate fake data
+    XYStatData data(nPoint, 1, 1);
+    RandGen    rg;
+    double     x_k;
+    auto       f = [](const double x[1], const double p[2])
+                     {return p[1]*exp(-x[0]*p[0]);};
+    
+    for (Index k = 0; k < nPoint; ++k)
+    {
+        x_k                = k*dx;
+        data.x(0, k)(0, 0) = x_k;
+        data.y(0, k)(0, 0) = f(&x_k, exactPar) + rg.gaussian(0.0, 0.1);
+    }
+    data.yyVar(0, 0).diagonal() = DMat::Constant(nPoint, 1, 0.1*0.1);
+    data.assumeXExact(0);
+    
+    // fit
+    DVec init = DVec::Constant(2, 0.5);
+    
+    DoubleModel model(1, 2, f);
+    FitResult p;
+    MinuitMinimizer minimizer;
+    
+    data.fitAllPoints();
+    p = data.fit(model, minimizer, init, true, Minimizer::Verbosity::Debug);
+    cout << "a= " << p(0) << " b= " << p(1)
+         << " chi^2/ndof= " << p.getChi2PerDof() << endl;
+    
+    // plot result
+    Plot plot;
+    
+    plot << LogScale(Axis::y) << PlotData(data);
+    plot.display();
+    
+    return EXIT_SUCCESS;
+}

latan/Chi2Function.cpp

@@ -0,0 +1,282 @@
+/*
+ * Chi2Function.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2014 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 <latan/Chi2Function.hpp>
+#include <latan/includes.hpp>
+#include <latan/XYStatData.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                         Chi2Function implementation                        *
+ ******************************************************************************/
+// constructors ////////////////////////////////////////////////////////////////
+Chi2Function::Chi2Function(const XYStatData &data)
+: DoubleFunction(0)
+, data_(data)
+, buffer_(new Chi2FunctionBuffer)
+{
+    resizeBuffer();
+}
+
+Chi2Function::Chi2Function(const XYStatData &data,
+                           const vector<const DoubleModel *> &modelVector)
+: Chi2Function(data)
+{
+    setModel(modelVector);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+Index Chi2Function::getNArg(void) const
+{
+    if (nPar_ < 0)
+    {
+        LATAN_ERROR(Memory, "no model set");
+    }
+    
+    return nPar_ + data_.getStatXDim()*data_.getNFitPoint();
+}
+
+Index Chi2Function::getNDof(void) const
+{
+    if (nPar_ < 0)
+    {
+        LATAN_ERROR(Memory, "no model set");
+    }
+    
+    return data_.getYDim()*data_.getNFitPoint() - nPar_;
+};
+
+Index Chi2Function::getNPar(void) const
+{
+    if (nPar_ < 0)
+    {
+        LATAN_ERROR(Memory, "no model set");
+    }
+    
+    return nPar_;
+};
+
+void Chi2Function::setModel(const DoubleModel &model, const Index j)
+{
+    if (static_cast<Index>(model_.size()) != data_.getYDim())
+    {
+         model_.resize(static_cast<unsigned int>(data_.getYDim()));
+    }
+    if (model.getNArg() != data_.getXDim())
+    {
+        LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
+    }
+    for (unsigned int l = 0; l < data_.getYDim(); ++l)
+    {
+        if (model_[l]&&(l != j))
+        {
+            if (model_[l]->getNPar() != model.getNPar())
+            {
+                LATAN_ERROR(Size, "model number of parameter mismatch");
+            }
+        }
+    }
+    model_[static_cast<unsigned int>(j)] = &model;
+    nPar_                                = model.getNPar();
+}
+
+void Chi2Function::setModel(const vector<const DoubleModel *> &modelVector)
+{
+    if (static_cast<Index>(model_.size()) != data_.getYDim())
+    {
+        model_.resize(static_cast<unsigned int>(data_.getYDim()));
+    }
+    if (modelVector.size() != static_cast<unsigned int>(data_.getYDim()))
+    {
+        LATAN_ERROR(Size, "number of models and y-dimension mismatch");
+    }
+    for (unsigned int j = 0; j < model_.size(); ++j)
+    {
+        if (!modelVector[j])
+        {
+            LATAN_ERROR(Memory, "trying to set a null model");
+        }
+        if (modelVector[j]->getNArg() != data_.getXDim())
+        {
+            LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
+        }
+        model_[static_cast<unsigned int>(j)] = modelVector[j];
+        if (modelVector[j]->getNPar() != modelVector[0]->getNPar())
+        {
+            LATAN_ERROR(Size, "model number of parameter mismatch");
+        }
+    }
+    nPar_ = modelVector[0]->getNPar();
+}
+
+void Chi2Function::requestInit(void) const
+{
+    buffer_->isInit = false;
+}
+
+void Chi2Function::resizeBuffer(void) const
+{
+    Index size;
+    
+    size = (data_.getYDim() + data_.getStatXDim())*data_.getNFitPoint();
+    buffer_->v.conservativeResize(size);
+    buffer_->x.resize(data_.getXDim());
+    buffer_->invVar.resize(size, size);
+    buffer_->xInd.resize(data_.getStatXDim());
+    buffer_->dInd.resize(data_.getNFitPoint());
+}
+
+// compute variance matrix inverse /////////////////////////////////////////////
+void Chi2Function::setVarianceBlock(const Index l1, const Index l2,
+                                    ConstBlock<DMatBase> m) const
+{
+    Index nPoint = data_.getNFitPoint();
+    
+    FOR_VEC(buffer_->dInd, k2)
+    FOR_VEC(buffer_->dInd, k1)
+    {
+        buffer_->invVar(l1*nPoint + k1, l2*nPoint + k2) = m(buffer_->dInd(k1),
+                                                            buffer_->dInd(k2));
+    }
+}
+
+void Chi2Function::initBuffer(void) const
+{
+    const Index xDim  = data_.getXDim();
+    const Index yDim  = data_.getYDim();
+    const Index nData = data_.getNData();
+    Index       is, kf;
+    
+    // resize buffer
+    resizeBuffer();
+    
+    // build index tables
+    is = 0;
+    for (Index i = 0; i < xDim; ++i)
+    {
+        if (!data_.isXExact(i))
+        {
+            buffer_->xInd(is) = i;
+            is++;
+        }
+    }
+    kf = 0;
+    for (Index k = 0; k < nData; ++k)
+    {
+        if (data_.isFitPoint(k))
+        {
+            buffer_->dInd(kf) = k;
+            kf++;
+        }
+    }
+    
+    // set y/y variance matrix
+    for (Index j2 = 0; j2 < yDim; ++j2)
+    for (Index j1 = 0; j1 < yDim; ++j1)
+    {
+        setVarianceBlock(j1, j2, data_.yyVar(j1, j2));
+    }
+
+    // set x/x variance matrix
+    FOR_VEC(buffer_->xInd, i2)
+    FOR_VEC(buffer_->xInd, i1)
+    {
+        setVarianceBlock(i1, i2, data_.xxVar(buffer_->xInd(i1),
+                                             buffer_->xInd(i2)));
+    }
+
+    // set y/x variance matrix
+    FOR_VEC(buffer_->xInd, i)
+    for (Index j = 0; j < yDim; ++j)
+    {
+        setVarianceBlock(j, i, data_.yxVar(j, buffer_->xInd(i)));
+    }
+
+    // inversion
+    buffer_->invVar = buffer_->invVar.inverse().eval();
+    buffer_->isInit = true;
+}
+
+// function call ///////////////////////////////////////////////////////////////
+double Chi2Function::operator()(const double *arg) const
+{
+    if (!model_[0])
+    {
+        LATAN_ERROR(Memory, "null model");
+    }
+    
+    const Index    xDim   = data_.getXDim();
+    const Index    yDim   = data_.getYDim();
+    const Index    nPoint = data_.getNFitPoint();
+    Index          is;
+    ConstMap<DVec> p(arg, nPar_, 1);
+    ConstMap<DVec> xi(arg + nPar_, getNArg() - nPar_, 1);
+    double         res;
+    
+    // initialize buffers if necessary
+    if (!buffer_->isInit)
+    {
+        initBuffer();
+    }
+    
+    // set the upper part of v: f_j(xi, p) - y_j
+    is = 0;
+    for (Index j = 0; j < yDim; ++j)
+    FOR_VEC(buffer_->dInd, k)
+    {
+        const DoubleModel *f = model_[static_cast<unsigned int>(j)];
+        double            f_jk, y_jk = data_.y(j, buffer_->dInd(k))(0, 0);
+        
+        if (!f)
+        {
+            LATAN_ERROR(Memory, "null model");
+        }
+        for (Index i = 0; i < xDim; ++i)
+        {
+            if (data_.isXExact(i))
+            {
+                buffer_->x(i) = data_.x(i, buffer_->dInd(k))(0, 0);
+            }
+            else
+            {
+                buffer_->x(i) = xi(is*nPoint + k);
+                is++;
+            }
+        }
+        f_jk = (*f)(buffer_->x, p);
+        buffer_->v(j*nPoint + k) = f_jk - y_jk;
+    }
+    
+    // set the down part of v: xi_ik - x_ik
+    FOR_VEC(buffer_->xInd, i)
+    FOR_VEC(buffer_->dInd, k)
+    {
+        double x_ik  = data_.x(buffer_->xInd(i), buffer_->dInd(k))(0, 0);
+        double xi_ik = xi(i*nPoint + k);
+        
+        buffer_->v(i*nPoint + k) = xi_ik - x_ik;
+    }
+    
+    // compute result
+    res = buffer_->v.dot(buffer_->invVar*buffer_->v);
+    
+    return res;
+}

latan/Chi2Function.hpp

@@ -0,0 +1,82 @@
+/*
+ * Chi2Function.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2014 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_Chi2Function_hpp_
+#define Latan_Chi2Function_hpp_
+
+#include <latan/Global.hpp>
+#include <latan/Function.hpp>
+#include <latan/Model.hpp>
+#include <memory>
+
+BEGIN_NAMESPACE
+
+/******************************************************************************
+ *                                 <Class>                                    *
+ ******************************************************************************/
+// forward declaration of XYStatData
+class XYStatData;
+
+class Chi2Function: public DoubleFunction
+{
+private:
+    struct Chi2FunctionBuffer
+    {
+        DVec       v, x;
+        DMat       invVar;
+        bool       isInit{false};
+        Vec<Index> xInd, dInd;
+    };
+public:
+    // constructor
+    Chi2Function(const XYStatData &data);
+    Chi2Function(const XYStatData &data,
+                 const std::vector<const DoubleModel *> &modelVector);
+    // destructor
+    virtual ~Chi2Function(void) = default;
+    // access
+    virtual Index getNArg(void) const;
+            Index getNDof(void) const;
+            Index getNPar(void) const;
+            void  setModel(const DoubleModel &model, const Index j = 0);
+            void  setModel(const std::vector<const DoubleModel *> &modelVector);
+            void  requestInit(void) const;
+    // function call
+    using DoubleFunction::operator();
+protected:
+    // function call
+    virtual double operator()(const double *arg) const;
+private:
+    // access
+    void resizeBuffer(void) const;
+    
+    // compute variance matrix inverse
+    void setVarianceBlock(const Index l1, const Index l2,
+                          ConstBlock<DMatBase> m) const;
+    void initBuffer(void) const;
+private:
+    const XYStatData                    &data_;
+    std::shared_ptr<Chi2FunctionBuffer> buffer_;
+    std::vector<const DoubleModel *>    model_;
+    Index                               nPar_{-1};
+};
+
+END_NAMESPACE
+
+#endif // Latan_Chi2Function_hpp_

latan/CompiledFunction.cpp

@@ -29,12 +29,12 @@ using namespace Latan;
  ******************************************************************************/
 // constructors ////////////////////////////////////////////////////////////////
 CompiledDoubleFunction::CompiledDoubleFunction(const unsigned nArg)
-: DoubleFunction(nArg)
+: DoubleFunction(nArg, nullptr)
 {}
 
 CompiledDoubleFunction::CompiledDoubleFunction(const unsigned nArg,
                                                const string &code)
-: DoubleFunction(nArg)
+: DoubleFunction(nArg, nullptr)
 {
     setCode(code);
 }

latan/CompiledFunction.hpp

@@ -37,9 +37,8 @@ class CompiledDoubleFunction: public DoubleFunction
 {
 public:
     // constructors
-    explicit CompiledDoubleFunction(const unsigned nArg);
+    CompiledDoubleFunction(const unsigned nArg);
-    explicit CompiledDoubleFunction(const unsigned nArg,
+    CompiledDoubleFunction(const unsigned nArg, const std::string &code);
-                                    const std::string &code);
     // destructor
     virtual ~CompiledDoubleFunction(void) = default;
     // access

latan/Exceptions.cpp

@@ -43,6 +43,7 @@ CONST_EXC(Size, Logic("size error: " + msg, loc))
 CONST_EXC(Runtime, runtime_error(Env::msgPrefix + msg + ERR_SUFF))
 CONST_EXC(Compilation, Runtime("compilation error: " + msg, loc))
 CONST_EXC(Io, Runtime("IO error: " + msg, loc))
+CONST_EXC(Memory , Runtime("memory error: " + msg, loc))
 CONST_EXC(Parsing, Runtime(msg, loc))
 CONST_EXC(Program, Runtime(msg, loc))
 CONST_EXC(Syntax, Runtime("syntax error: " + msg, loc))

latan/Exceptions.hpp

@@ -53,6 +53,7 @@ namespace Exceptions
     DECL_EXC(Runtime, std::runtime_error);
     DECL_EXC(Compilation, Runtime);
     DECL_EXC(Io, Runtime);
+    DECL_EXC(Memory, Runtime);
     DECL_EXC(Parsing, Runtime);
     DECL_EXC(Program, Runtime);
     DECL_EXC(Syntax, Runtime);

latan/Function.cpp

@@ -48,9 +48,10 @@ const DoubleFunction::vecFunc DoubleFunction::nullFunction_ = nullptr;
 
 // constructor /////////////////////////////////////////////////////////////////
 DoubleFunction::DoubleFunction(const Index nArg, const vecFunc &f)
-: buffer_(new DVec(nArg))
+: buffer_(new DVec)
-, f_(f)
+{
-{}
+    setFunction(f, nArg);
+}
 
 // access //////////////////////////////////////////////////////////////////////
 Index DoubleFunction::getNArg(void) const
@@ -58,14 +59,21 @@ Index DoubleFunction::getNArg(void) const
     return buffer_->size();
 }
 
-void DoubleFunction::setFunction(const vecFunc &f)
+void DoubleFunction::setFunction(const vecFunc &f, const Index nArg)
 {
+    buffer_->resize(nArg);
     f_ = f;
 }
 
-void DoubleFunction::setNArg(const Index nArg)
+// error checking //////////////////////////////////////////////////////////////
+void DoubleFunction::checkSize(const Index nPar) const
 {
-    buffer_->resize(nArg);
+    if (nPar != getNArg())
+    {
+        LATAN_ERROR(Size, "function argument vector has a wrong size (expected "
+                    + strFrom(getNArg()) + ", got " + strFrom(nPar)
+                    + ")");
+    }
 }
 
 // function call ///////////////////////////////////////////////////////////////
@@ -76,38 +84,23 @@ double DoubleFunction::operator()(const double *arg) const
 
 double DoubleFunction::operator()(const DVec &arg) const
 {
-    if (arg.size() != getNArg())
+    checkSize(arg.size());
-    {
-        LATAN_ERROR(Size, "function argument vector has a wrong size (expected "
-                    + strFrom(getNArg()) + ", got " + strFrom(arg.size())
-                    + ")");
-    }
     
     return (*this)(arg.data());
 }
 
 double DoubleFunction::operator()(const std::vector<double> &arg) const
 {
-    if (arg.size() != static_cast<unsigned int>(getNArg()))
+    checkSize(static_cast<Index>(arg.size()));
-    {
-        LATAN_ERROR(Size, "function argument vector has a wrong size (expected "
-                    + strFrom(getNArg()) + ", got " + strFrom(arg.size())
-                    + ")");
-    }
     
     return (*this)(arg.data());
 }
 
 double DoubleFunction::operator()(std::stack<double> &arg) const
 {
+    checkSize(static_cast<Index>(arg.size()));
     for (Index i = 0; i < getNArg(); ++i)
     {
-        if (arg.empty())
-        {
-            LATAN_ERROR(Size, "function argument stack has a wrong size (expected "
-                        + strFrom(getNArg()) + ", got " + strFrom(i)
-                        + ")");
-        }
         (*buffer_)(getNArg() - i - 1) = arg.top();
         arg.pop();
     }

latan/Function.hpp

@@ -42,9 +42,8 @@ public:
     // destructor
     virtual ~DoubleFunction(void) = default;
     // access
-    Index getNArg(void) const;
+    virtual Index getNArg(void) const;
-    void  setFunction(const vecFunc &f);
+    void    setFunction(const vecFunc &f, const Index nArg);
-    void  setNArg(const Index nArg);
     // function call
     double operator()(const DVec &arg) const;
     double operator()(const std::vector<double> &arg) const;
@@ -54,7 +53,10 @@ protected:
     // function call
     virtual double operator()(const double *arg) const;
 private:
-    std::shared_ptr<DVec> buffer_;
+    // error checking
+    void checkSize(const Index nPar) const;
+private:
+    std::shared_ptr<DVec> buffer_{nullptr};
     vecFunc               f_;
     static const vecFunc  nullFunction_;
 };

latan/Global.hpp

@@ -76,6 +76,9 @@ using Vec = Mat<T, dynamic, 1>;
 typedef Vec<int>    IVec;
 typedef Vec<double> DVec;
 
+#define FOR_VEC(vec, i) \
+for (Latan::Index i = 0; i < (vec).size(); ++i)
+
 // block types
 template <typename Derived>
 using Block      = Eigen::Block<Derived>;
@@ -129,7 +132,7 @@ inline T strTo(const std::string &str)
     return buf;
 }
 
-//// optimized specializations
+// optimized specializations
 template <>
 inline float strTo<float>(const std::string &str)
 {

latan/Makefile.am

@@ -28,6 +28,7 @@ liblatan_la_SOURCES =   \
     AsciiFile.cpp       \
     AsciiParser.ypp     \
     AsciiLexer.lpp      \
+    Chi2Function.cpp    \
     CompiledFunction.cpp\
     Exceptions.cpp      \
     Function.cpp        \
@@ -40,13 +41,16 @@ liblatan_la_SOURCES =   \
     MathParser.ypp      \
     MathLexer.lpp       \
     Minimizer.cpp       \
+    Model.cpp           \
     Plot.cpp            \
     RandGen.cpp         \
     Sample.cpp          \
+    XYStatData.cpp      \
 	../config.h
 liblatan_ladir = $(pkgincludedir)
 liblatan_la_HEADERS =   \
     AsciiFile.hpp       \
+    Chi2Function.hpp    \
     CompiledFunction.hpp\
     Dataset.hpp         \
     Function.hpp        \
@@ -57,12 +61,13 @@ liblatan_la_HEADERS =   \
     Math.hpp            \
     MathInterpreter.hpp \
     Minimizer.hpp       \
+    Model.hpp           \
     ParserState.hpp     \
     Plot.hpp            \
     RandGen.hpp         \
     Sample.hpp          \
-    StatArray.hpp
+    StatArray.hpp       \
-    
+    XYStatData.hpp
 if HAVE_MINUIT
 liblatan_la_SOURCES += MinuitMinimizer.cpp
 liblatan_la_HEADERS += MinuitMinimizer.hpp

latan/Mat.hpp

@@ -24,8 +24,8 @@
 #include <latan/IOObject.hpp>
 
 #define FOR_MAT(mat, i, j) \
-for (Index j = 0; j < mat.cols(); ++j)\
+for (Latan::Index j = 0; j < mat.cols(); ++j)\
-for (Index i = 0; i < mat.rows(); ++i)
+for (Latan::Index i = 0; i < mat.rows(); ++i)
 
 BEGIN_NAMESPACE
 

latan/MinuitMinimizer.cpp

@@ -76,8 +76,8 @@ const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
     }
     
     // pre-minimization
-    MnMigrad migrad0(minuitF, parameters, 0);
+    MnMigrad        migrad1(minuitF, parameters, 1);
-    FunctionMinimum min = migrad0();
+    FunctionMinimum min = migrad1();
     
     if (verbosity >= Verbosity::Debug)
     {

latan/Model.cpp

@@ -0,0 +1,91 @@
+/*
+ * Model.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2014 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 <latan/Model.hpp>
+#include <latan/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                           Model implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+DoubleModel::DoubleModel(const Index nArg, const Index nPar, const vecFunc &f)
+: size_(new ModelSize)
+{
+    setFunction(nArg, nPar, f);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+Index DoubleModel::getNArg(void) const
+{
+    return size_->nArg;
+}
+
+Index DoubleModel::getNPar(void) const
+{
+    return size_->nPar;
+}
+
+void DoubleModel::setFunction(const Index nArg, const Index nPar,
+                              const vecFunc &f)
+{
+    size_->nArg = nArg;
+    size_->nPar = nPar;
+    f_          = f;
+}
+
+// error checking //////////////////////////////////////////////////////////////
+void DoubleModel::checkSize(const Index nArg, const Index nPar) const
+{
+    if (nArg != getNArg())
+    {
+        LATAN_ERROR(Size, "model argument vector has a wrong size (expected "
+                    + strFrom(getNArg()) + ", got " + strFrom(nArg)
+                    + ")");
+    }
+    if (nPar != getNPar())
+    {
+        LATAN_ERROR(Size, "model parameter vector has a wrong size (expected "
+                    + strFrom(getNPar()) + ", got " + strFrom(nPar)
+                    + ")");
+    }
+}
+
+// function call ///////////////////////////////////////////////////////////////
+double DoubleModel::operator()(const DVec &arg, const DVec &par) const
+{
+    checkSize(arg.size(), par.size());
+    
+    return f_(arg.data(), par.data());
+}
+
+double DoubleModel::operator()(const vector<double> &arg,
+                               const vector<double> &par) const
+{
+    checkSize(static_cast<Index>(arg.size()), static_cast<Index>(par.size()));
+    
+    return f_(arg.data(), par.data());
+}
+
+double DoubleModel::operator()(const double *data, const double *par) const
+{
+    return f_(data, par);
+}

latan/Model.hpp

@@ -0,0 +1,67 @@
+/*
+ * Model.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2014 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_Model_hpp_
+#define Latan_Model_hpp_
+
+#include <latan/Global.hpp>
+#include <latan/Mat.hpp>
+#include <memory>
+#include <vector>
+
+BEGIN_NAMESPACE
+
+/******************************************************************************
+ *                           Double model class                               *
+ ******************************************************************************/
+class DoubleModel
+{
+private:
+    typedef std::function<double(const double *, const double *)> vecFunc;
+    struct ModelSize {Index nArg, nPar;};
+public:
+    // constructor
+    DoubleModel(const Index nArg = 0, const Index nPar = 0,
+                const vecFunc &f = nullFunction_);
+    // destructor
+    virtual ~DoubleModel(void) = default;
+    // access
+    virtual Index getNArg(void) const;
+    virtual Index getNPar(void) const;
+    void  setFunction(const Index nArg = 0, const Index nPar = 0,
+                      const vecFunc &f = nullFunction_);
+    // function call
+    double operator()(const DVec &data, const DVec &par) const;
+    double operator()(const std::vector<double> &data,
+                      const std::vector<double> &par) const;
+protected:
+    // function call
+    virtual double operator()(const double *data, const double *par) const;
+private:
+    // error checking
+    void checkSize(const Index nArg, const Index nPar) const;
+private:
+    std::shared_ptr<ModelSize> size_;
+    vecFunc                    f_;
+    static const vecFunc       nullFunction_;
+};
+
+END_NAMESPACE
+
+#endif // Latan_Model_hpp_

latan/Plot.cpp

@@ -17,29 +17,127 @@
  * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#define _POSIX_C_SOURCE 199209L
-
 #include <latan/Plot.hpp>
 #include <latan/includes.hpp>
+#include <latan/Mat.hpp>
 
 using namespace std;
 using namespace Latan;
 
 /******************************************************************************
- *                         PlotCommand implementation                         *
+ *                              Plot objects                                  *
  ******************************************************************************/
-// constructors ////////////////////////////////////////////////////////////////
+// PlotObject access ///////////////////////////////////////////////////////////
-PlotCommand::PlotCommand(void)
+const string & PlotObject::getCommand(void) const
-{}
+{
+    return command_;
+}
 
+void PlotObject::setCommand(const string &command)
+{
+    command_ = command;
+}
+
+string PlotObject::popTmpFile(void)
+{
+    string res = tmpFileName_.top();
+    
+    tmpFileName_.pop();
+    
+    return res;
+}
+
+void PlotObject::pushTmpFile(const std::string &fileName)
+{
+    tmpFileName_.push(fileName);
+}
+
+
+// PlotObject test /////////////////////////////////////////////////////////////
+bool PlotObject::gotTmpFile(void) const
+{
+    return !tmpFileName_.empty();
+}
+
+// PlotCommand constructor /////////////////////////////////////////////////////
 PlotCommand::PlotCommand(const string &command)
 : command_(command)
 {}
 
-// access //////////////////////////////////////////////////////////////////////
+// PlotData constructor ////////////////////////////////////////////////////////
-const std::string & PlotCommand::getCommand(void) const
+PlotData::PlotData(const XYStatData &data, const int i, const int j)
+: data_(data)
+, i_(i)
+, j_(j)
 {
-    return command_;
+    DMat d(data_.getNData(), 4);
+    char tmpFileName[MAX_PATH_LENGTH];
+    int fd;
+    FILE *tmpFile;
+    string usingCmd;
+    
+    d.col(0) = data_.x(i_, _);
+    d.col(2) = data_.y(j_, _);
+    d.col(1) = data_.xxVar(i_, i_).diagonal().array().sqrt();
+    d.col(3) = data_.yyVar(i_, i_).diagonal().array().sqrt();
+    usingCmd = (data_.isXExact(i_)) ? "u 1:3:4 w yerr" : "u 1:2:3:4 w xyerr";
+    strcpy(tmpFileName, "./latan_plot_tmp.XXXXXX.dat");
+    fd = mkstemps(tmpFileName, 4);
+    if (fd == -1)
+    {
+        LATAN_ERROR(System, "impossible to create a temporary file from template "
+                    + strFrom(tmpFileName));
+    }
+    tmpFile = fdopen(fd, "w");
+    for (Index r = 0; r < data_.getNData(); ++r)
+    {
+        fprintf(tmpFile, "%e %e %e %e\n", d(r, 0), d(r, 1), d(r, 2), d(r, 3));
+    }
+    fclose(tmpFile);
+    pushTmpFile(tmpFileName);
+    setCommand("'" + string(tmpFileName) + "' " + usingCmd);
+}
+
+/******************************************************************************
+ *                             Plot modifiers                                 *
+ ******************************************************************************/
+// Color constructor ///////////////////////////////////////////////////////////
+Color::Color(const string &color)
+: color_(color)
+{}
+
+// Color modifier //////////////////////////////////////////////////////////////
+void Color::operator()(PlotOptions &option) const
+{
+    option.lineColor = color_;
+}
+
+// LogScale constructor ////////////////////////////////////////////////////////
+LogScale::LogScale(const Axis axis)
+: axis_(axis)
+{}
+
+// Logscale modifier ///////////////////////////////////////////////////////////
+void LogScale::operator()(PlotOptions &option) const
+{
+    option.scaleMode[static_cast<int>(axis_)] |= Plot::Scale::log;
+}
+
+// PlotRange constructor ////////////////////////////////////////////////////////
+PlotRange::PlotRange(const Axis axis, const double min, const double max)
+: axis_(axis)
+, min_(min)
+, max_(max)
+{}
+
+// PlotRange modifier ///////////////////////////////////////////////////////////
+void PlotRange::operator()(PlotOptions &option) const
+{
+    int a = static_cast<int>(axis_);
+    
+    option.scaleMode[a] |= Plot::Scale::manual;
+    option.scale[a].min  = min_;
+    option.scale[a].max  = max_;
 }
 
 /******************************************************************************
@@ -47,21 +145,7 @@ const std::string & PlotCommand::getCommand(void) const
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
 Plot::Plot(void)
-: gnuplotBin_(GNUPLOT_BIN)
+{}
-, gnuplotArgs_(GNUPLOT_ARGS)
-, gnuplotPath_("")
-, terminal_("")
-, title_("")
-{
-    scaleMode_[Axis::x] = 0;
-    scaleMode_[Axis::y] = 0;
-    scale_[Axis::x].min = 0.0;
-    scale_[Axis::x].max = 0.0;
-    scale_[Axis::y].min = 0.0;
-    scale_[Axis::y].max = 0.0;
-    label_[Axis::x]     = "";
-    label_[Axis::y]     = "";
-}
 
 // destructor //////////////////////////////////////////////////////////////////
 Plot::~Plot(void)
@@ -78,9 +162,28 @@ Plot::~Plot(void)
 }
 
 // plot objects ////////////////////////////////////////////////////////////////
-Plot & Plot::operator<<(const PlotCommand &command)
+Plot & Plot::operator<<(PlotObject &&command)
 {
-    plotCommand_.push_back(command.getCommand());
+    string commandStr;
+    
+    while (command.gotTmpFile())
+    {
+        tmpFileName_.push(command.popTmpFile());
+    }
+    commandStr = command.getCommand();
+    if (!options_.lineColor.empty())
+    {
+        commandStr += " lc " + options_.lineColor;
+        options_.lineColor = "";
+    }
+    plotCommand_.push_back(commandStr);
+    
+    return *this;
+}
+
+Plot & Plot::operator<<(PlotModifier &&modifier)
+{
+    modifier(options_);
     
     return *this;
 }
@@ -178,52 +281,53 @@ void Plot::display(void)
 ostream & Latan::operator<<(ostream &out, const Plot &plot)
 {
     std::string begin, end;
+    int x = static_cast<int>(Axis::x), y = static_cast<int>(Axis::y);
     
-    if (!plot.terminal_.empty())
+    if (!plot.options_.terminal.empty())
     {
-        out << "set term " << plot.terminal_ << endl;
+        out << "set term " << plot.options_.terminal << endl;
     }
-    if (!plot.output_.empty())
+    if (!plot.options_.output.empty())
     {
-        out << "set output '" << plot.terminal_ << "'" << endl;
+        out << "set output '" << plot.options_.terminal << "'" << endl;
     }
-    if (plot.scaleMode_[Plot::Axis::x] & Plot::Scale::manual)
+    if (plot.options_.scaleMode[x] & Plot::Scale::manual)
     {
-        out << "xMin = " << plot.scale_[Plot::Axis::x].min << endl;
+        out << "xMin = " << plot.options_.scale[x].min << endl;
-        out << "xMax = " << plot.scale_[Plot::Axis::x].max << endl;
+        out << "xMax = " << plot.options_.scale[x].max << endl;
     }
-    if (plot.scaleMode_[Plot::Axis::y] & Plot::Scale::manual)
+    if (plot.options_.scaleMode[y] & Plot::Scale::manual)
     {
-        out << "yMin = " << plot.scale_[Plot::Axis::y].min << endl;
+        out << "yMin = " << plot.options_.scale[y].min << endl;
-        out << "yMax = " << plot.scale_[Plot::Axis::y].max << endl;
+        out << "yMax = " << plot.options_.scale[y].max << endl;
     }
-    if (!plot.title_.empty())
+    if (!plot.options_.title.empty())
     {
-        out << "set title '" << plot.title_ << "'" << endl;
+        out << "set title '" << plot.options_.title << "'" << endl;
     }
-    if (plot.scaleMode_[Plot::Axis::x] & Plot::Scale::manual)
+    if (plot.options_.scaleMode[x] & Plot::Scale::manual)
     {
         out << "set xrange [xMin:xMax]" << endl;
     }
-    if (plot.scaleMode_[Plot::Axis::y] & Plot::Scale::manual)
+    if (plot.options_.scaleMode[y] & Plot::Scale::manual)
     {
         out << "set yrange [yMin:yMax]" << endl;
     }
-    if (plot.scaleMode_[Plot::Axis::x] & Plot::Scale::log)
+    if (plot.options_.scaleMode[x] & Plot::Scale::log)
     {
         out << "set log x" << endl;
     }
-    if (plot.scaleMode_[Plot::Axis::y] & Plot::Scale::log)
+    if (plot.options_.scaleMode[y] & Plot::Scale::log)
     {
         out << "set log y" << endl;
     }
-    if (!plot.label_[Plot::Axis::x].empty())
+    if (!plot.options_.label[x].empty())
     {
-        out << "set xlabel '" << plot.label_[Plot::Axis::x] << "'" << endl;
+        out << "set xlabel '" << plot.options_.label[x] << "'" << endl;
     }
-    if (!plot.label_[Plot::Axis::y].empty())
+    if (!plot.options_.label[y].empty())
     {
-        out << "set ylabel '" << plot.label_[Plot::Axis::y] << "'" << endl;
+        out << "set ylabel '" << plot.options_.label[y] << "'" << endl;
     }
     for (unsigned int i = 0; i < plot.headCommand_.size(); ++i)
     {

latan/Plot.hpp

@@ -24,6 +24,7 @@
 #include <stack>
 #include <vector>
 #include <latan/Global.hpp>
+#include <latan/XYStatData.hpp>
 
 // gnuplot default parameters
 #ifndef GNUPLOT_BIN
@@ -37,26 +38,125 @@
 BEGIN_NAMESPACE
 
 /******************************************************************************
- *                             Plot commands                                  *
+ *                              Plot objects                                  *
  ******************************************************************************/
-class PlotCommand
+class PlotObject
+{
+public:
+    // destructor
+    virtual ~PlotObject(void) = default;
+    // access
+    std::string                 popTmpFile(void);
+    const std::string &         getCommand(void) const;
+    // test
+    bool                        gotTmpFile(void) const;
+protected:
+    // access
+    void pushTmpFile(const std::string &fileName);
+    void setCommand(const std::string &command);
+private:
+    // plot command
+    std::string command_;
+    // stack of created temporary files
+    std::stack<std::string>  tmpFileName_;
+};
+
+class PlotCommand: public PlotObject
 {
 public:
     // constructors
-    PlotCommand(void);
+    explicit PlotCommand(const std::string &command);
-    PlotCommand(const std::string &command);
     // destructor
     virtual ~PlotCommand(void) = default;
-    // access
+private:
-    virtual const std::string & getCommand(void) const;
-protected:
     std::string command_;
 };
 
+class PlotData: public PlotObject
+{
+public:
+    // constructors
+    explicit PlotData(const XYStatData &data, const int i = 0, const int j = 0);
+    // destructor
+    virtual ~PlotData(void) = default;
+private:
+    const XYStatData &data_;
+    const int i_, j_;
+};
+
+/******************************************************************************
+ *                             Plot modifiers                                 *
+ ******************************************************************************/
+enum class Axis {x = 0, y = 1};
+
+struct Range
+{
+    double min, max;
+};
+
+struct PlotOptions
+{
+    std::string  terminal     {""};
+    std::string  output       {""};
+    std::string  title        {""};
+    unsigned int scaleMode[2] {0,0};
+    Range        scale[2]     {{0.0,0.0},{0.0,0.0}};
+    std::string  label[2]     {"",""};
+    std::string  lineColor    {""};
+};
+
+class PlotModifier
+{
+public:
+    // destructor
+    virtual ~PlotModifier(void) = default;
+    // modifier
+    virtual void operator()(PlotOptions &option) const = 0;
+};
+
+class Color: public PlotModifier
+{
+public:
+    // constructor
+    explicit Color(const std::string &color);
+    // destructor
+    virtual ~Color(void) = default;
+    // modifier
+    virtual void operator()(PlotOptions &option) const;
+private:
+    const std::string color_;
+};
+
+class LogScale: public PlotModifier
+{
+public:
+    // constructor
+    explicit LogScale(const Axis axis);
+    // destructor
+    virtual ~LogScale(void) = default;
+    // modifier
+    virtual void operator()(PlotOptions &option) const;
+private:
+    const Axis axis_;
+};
+
+class PlotRange: public PlotModifier
+{
+public:
+    // constructor
+    explicit PlotRange(const Axis axis, const double min, const double max);
+    // destructor
+    virtual ~PlotRange(void) = default;
+    // modifier
+    virtual void operator()(PlotOptions &option) const;
+private:
+    const Axis   axis_;
+    const double min_, max_;
+};
+
 /******************************************************************************
  *                               Plot class                                   *
  ******************************************************************************/
-
 class Plot
 {
 public:
@@ -70,26 +170,13 @@ public:
             log    = 1 << 1
         };
     };
-private:
-    struct Range
-    {
-        double min, max;
-    };
-    class Axis
-    {
-    public:
-        enum
-        {
-            x = 0,
-            y = 1
-        };
-    };
 public:
     // constructor/destructor
     Plot(void);
     virtual ~Plot(void);
-    // plot objects
+    // plot operations
-    Plot & operator<<(const PlotCommand &command);
+    Plot & operator<<(PlotObject   &&command);
+    Plot & operator<<(PlotModifier &&modifier);
     // plot parsing and output
     void display(void);
     friend std::ostream & operator<<(std::ostream &out, const Plot &plot);
@@ -98,20 +185,15 @@ private:
     void getProgramPath(void);
 private:
     // gnuplot execution parameters
-    std::string              gnuplotBin_;
+    std::string              gnuplotBin_  {GNUPLOT_BIN};
-    std::string              gnuplotArgs_;
+    std::string              gnuplotArgs_ {GNUPLOT_ARGS};
-    std::string              gnuplotPath_;
+    std::string              gnuplotPath_ {""};
     // string buffer for commands
     std::ostringstream       commandBuffer_;
     // stack of created temporary files
     std::stack<std::string>  tmpFileName_;
     // plot content
-    std::string              terminal_;
+    PlotOptions              options_;
-    std::string              output_;
-    std::string              title_;
-    unsigned int             scaleMode_[2];
-    Range                    scale_[2];
-    std::string              label_[2];
     std::vector<std::string> headCommand_;
     std::vector<std::string> plotCommand_;
 };

latan/StatArray.hpp

@@ -45,8 +45,8 @@ public:
     // access
     Index size(void) const;
     // operators
-          T & operator[](const int s);
+          T & operator[](const Index s);
-    const T & operator[](const int s) const;
+    const T & operator[](const Index s) const;
     // statistics
     void bin(Index binSize);
     T    mean(const Index pos, const Index n) const;
@@ -89,13 +89,13 @@ Index StatArray<T, offset>::size(void) const
 
 // operators ///////////////////////////////////////////////////////////////////
 template <typename T, unsigned int offset>
-T & StatArray<T, offset>::operator[](const int s)
+T & StatArray<T, offset>::operator[](const Index s)
 {
     return Base::operator[](s + offset);
 }
 
 template <typename T, unsigned int offset>
-const T & StatArray<T, offset>::operator[](const int s) const
+const T & StatArray<T, offset>::operator[](const Index s) const
 {
     return Base::operator[](s + offset);
 }

latan/XYStatData.cpp

@@ -0,0 +1,273 @@
+/*
+ * XYStatData.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2014 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 <latan/XYStatData.hpp>
+#include <latan/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                          FitResult implementation                          *
+ ******************************************************************************/
+// access //////////////////////////////////////////////////////////////////////
+double FitResult::getChi2(void) const
+{
+    return chi2_;
+}
+
+double FitResult::getChi2PerDof(void) const
+{
+    return chi2_/static_cast<double>(nDof_);
+}
+
+/******************************************************************************
+ *                         XYStatData implementation                          *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+XYStatData::XYStatData(void)
+: chi2_(*this)
+{}
+
+XYStatData::XYStatData(const Index nData, const Index xDim, const Index yDim)
+: XYStatData()
+{
+    resize(nData, xDim, yDim);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+void XYStatData::assumeXExact(const Index i, const bool isExact)
+{
+    isXExact_[i] = (isExact) ? 1 : 0;
+}
+
+void XYStatData::fitPoint(const Index i, const bool isFitPoint)
+{
+    isFitPoint_[i] = (isFitPoint) ? 1 : 0;
+}
+
+void XYStatData::fitPointRange(const Index k1, const Index k2,
+                               const bool isFitPoint)
+{
+    int size = static_cast<int>(k2-k1+1);
+    
+    isFitPoint_.segment(k1, size) = IVec::Constant(size, (isFitPoint) ? 1 : 0);
+}
+
+void XYStatData::fitAllPoints(const bool isFitPoint)
+{
+    fitPointRange(0, getNData()-1, isFitPoint);
+}
+
+Index XYStatData::getNData(void) const
+{
+    return x_.rows();
+}
+
+Index XYStatData::getNFitPoint(void) const
+{
+    return isFitPoint_.sum();
+}
+
+Index XYStatData::getXDim(void) const
+{
+    return x_.cols();
+}
+
+Index XYStatData::getYDim(void) const
+{
+    return y_.cols();
+}
+
+Index XYStatData::getStatXDim(void) const
+{
+    return isXExact_.size() - isXExact_.sum();
+}
+
+void XYStatData::setNData(const Index nData)
+{
+    resize(nData, getXDim(), getYDim());
+}
+
+void XYStatData::setXDim(const Index xDim)
+{
+    resize(getNData(), xDim, getYDim());
+}
+
+void XYStatData::setYDim(const Index yDim)
+{
+    resize(getNData(), getXDim(), yDim);
+}
+
+void XYStatData::resize(const Index nData, const Index xDim, const Index yDim)
+{
+    x_.conservativeResize(nData, xDim);
+    y_.conservativeResize(nData, yDim);
+    isXExact_.conservativeResize(xDim);
+    isFitPoint_.conservativeResize(nData);
+    var_[xx].conservativeResize(xDim, xDim);
+    var_[yy].conservativeResize(yDim, yDim);
+    var_[yx].conservativeResize(yDim, xDim);
+    FOR_MAT(var_[xx], i1, i2)
+    {
+        var_[xx](i1, i2).conservativeResize(nData, nData);
+    }
+    FOR_MAT(var_[yy], j1, j2)
+    {
+        var_[yy](j1, j2).conservativeResize(nData, nData);
+    }
+    FOR_MAT(var_[yx], j, i)
+    {
+        var_[yx](j, i).conservativeResize(nData, nData);
+    }
+}
+
+#define FULL_BLOCK(m) (m).block(0, 0, (m).rows(), (m).cols())
+#define ACCESS_DATA(xy, ij, k) \
+if ((ij >= 0)&&(k >= 0))\
+{\
+    return xy.block(k, ij, 1, 1);\
+}\
+else if ((ij < 0)&&(k >= 0))\
+{\
+    return xy.block(k, 0, 1, getXDim());\
+}\
+else if ((ij >= 0)&&(k < 0))\
+{\
+    return xy.block(0, ij, getNData(), 1);\
+}\
+else\
+{\
+    return xy.block(0, 0, getNData(), getXDim());\
+}
+
+Block<DMatBase> XYStatData::x(const Index i, const Index k)
+{
+    ACCESS_DATA(x_, i, k);
+}
+
+ConstBlock<DMatBase> XYStatData::x(const Index i, const Index k) const
+{
+    ACCESS_DATA(x_, i, k);
+}
+
+Block<DMatBase> XYStatData::y(const Index j, const Index k)
+{
+    ACCESS_DATA(y_, j, k);
+}
+
+ConstBlock<DMatBase> XYStatData::y(const Index j, const Index k) const
+{
+    ACCESS_DATA(y_, j, k);
+}
+
+Block<DMatBase> XYStatData::xxVar(const Index i1, const Index i2)
+{
+    return FULL_BLOCK(var_[xx](i1, i2));
+}
+
+ConstBlock<DMatBase> XYStatData::xxVar(const Index i1, const Index i2) const
+{
+    return FULL_BLOCK(var_[xx](i1, i2));
+}
+
+Block<DMatBase> XYStatData::yyVar(const Index j1, const Index j2)
+{
+    return FULL_BLOCK(var_[yy](j1, j2));
+}
+
+ConstBlock<DMatBase> XYStatData::yyVar(const Index j1, const Index j2) const
+{
+    return FULL_BLOCK(var_[yy](j1, j2));
+}
+
+Block<DMatBase> XYStatData::yxVar(const Index j, const Index i)
+{
+    return FULL_BLOCK(var_[yx](j, i));
+}
+
+ConstBlock<DMatBase> XYStatData::yxVar(const Index j, const Index i) const
+{
+    return FULL_BLOCK(var_[yx](j, i));
+}
+
+// test ////////////////////////////////////////////////////////////////////////
+bool XYStatData::isFitPoint(const Index k) const
+{
+    return (isFitPoint_[k] == 1);
+}
+
+bool XYStatData::isXExact(const Index i) const
+{
+    return (isXExact_[i] == 1);
+}
+
+// fit /////////////////////////////////////////////////////////////////////////
+FitResult XYStatData::fit(const vector<const DoubleModel *> &modelVector,
+                          Minimizer &minimizer, const DVec &init,
+                          const bool reinitChi2,
+                          const FitVerbosity verbosity)
+{
+    // initialization
+    chi2_.setModel(modelVector);
+    if (reinitChi2)
+    {
+        chi2_.requestInit();
+    }
+    minimizer.setVerbosity(verbosity);
+    
+    // initial parameters
+    const Index nPoint = getNFitPoint();
+    DVec        fullInit = init;
+    Index       is = 0, kf = 0;
+    
+    fullInit.conservativeResize(chi2_.getNArg());
+    for (Index i = 0; i < getXDim(); ++i)
+    if (!isXExact(i))
+    {
+        for (Index k = 0; k < getNData(); ++k)
+        if (isFitPoint(k))
+        {
+            fullInit(chi2_.getNPar() + nPoint*is + kf) = x(i, k)(0, 0);
+            kf++;
+        }
+        is++;
+    }
+    minimizer.setInit(fullInit);
+    
+    // fit
+    FitResult result;
+    
+    result       = minimizer(chi2_);
+    result.chi2_ = chi2_(result);
+    result.nDof_ = chi2_.getNDof();
+    
+    return result;
+}
+
+FitResult XYStatData::fit(const DoubleModel &model, Minimizer &minimizer,
+                          const DVec &init, const bool reinitChi2,
+                          const FitVerbosity verbosity)
+{
+    vector<const DoubleModel *> modelVector(1);
+    
+    modelVector[0] = &model;
+    
+    return fit(modelVector, minimizer, init, reinitChi2, verbosity);
+}

latan/XYStatData.hpp

@@ -0,0 +1,123 @@
+/*
+ * XYData.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2014 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_XYData_hpp_
+#define Latan_XYData_hpp_
+
+#include <latan/Global.hpp>
+#include <latan/Chi2Function.hpp>
+#include <latan/Function.hpp>
+#include <latan/Mat.hpp>
+#include <latan/Minimizer.hpp>
+#include <latan/Model.hpp>
+#include <memory>
+
+BEGIN_NAMESPACE
+
+/******************************************************************************
+ *                           object for fit result                            *
+ ******************************************************************************/
+class FitResult: public DVec
+{
+    friend class XYStatData;
+public:
+    // constructors
+    using DVec::DVec;
+    // destructor
+    virtual ~FitResult(void) = default;
+    // access
+    double getChi2(void) const;
+    double getChi2PerDof(void) const;
+private:
+    double chi2_{0.0};
+    Index  nDof_{0};
+};
+
+/******************************************************************************
+ *                    object for X vs. Y statistical data                     *
+ ******************************************************************************/
+const Index _ = -1;
+
+class XYStatData
+{
+public:
+    enum
+    {
+        xx = 0,
+        yy = 1,
+        yx = 2
+    };
+    typedef Minimizer::Verbosity FitVerbosity;
+public:
+    // constructors
+    XYStatData(void);
+    XYStatData(const Index nData, const Index nXDim, const Index nYDim);
+    // destructor
+    virtual ~XYStatData(void) = default;
+    // access
+    void                 assumeXExact(const Index i, const bool isExact = true);
+    void                 fitPoint(const Index k, const bool isFitPoint = true);
+    void                 fitPointRange(const Index k1, const Index k2,
+                                       const bool isFitPoint = true);
+    void                 fitAllPoints(const bool isFitPoint = true);
+    Index                getNData(void)     const;
+    Index                getNFitPoint(void) const;
+    Index                getXDim(void)      const;
+    Index                getYDim(void)      const;
+    Index                getStatXDim(void)  const;
+    void                 setNData(const Index nData);
+    void                 setXDim(const Index xDim);
+    void                 setYDim(const Index yDim);
+    void                 resize(const Index nData, const Index xDim,
+                                const Index yDim);
+    Block<DMatBase>      x(const Index i = _, const Index k = _);
+    ConstBlock<DMatBase> x(const Index i = _, const Index k = _) const;
+    Block<DMatBase>      y(const Index j = _, const Index k = _);
+    ConstBlock<DMatBase> y(const Index j = _, const Index k = _) const;
+    Block<DMatBase>      xxVar(const Index i1, const Index i2);
+    ConstBlock<DMatBase> xxVar(const Index i1, const Index i2) const;
+    Block<DMatBase>      yyVar(const Index j1, const Index j2);
+    ConstBlock<DMatBase> yyVar(const Index j1, const Index j2) const;
+    Block<DMatBase>      yxVar(const Index j,  const Index i);
+    ConstBlock<DMatBase> yxVar(const Index j,  const Index i) const;
+    // test
+    bool isFitPoint(const Index k) const;
+    bool isXExact(const Index i) const;
+    // fit
+    FitResult fit(const std::vector<const DoubleModel *> &modelVector,
+                  Minimizer &minimizer, const DVec &init,
+                  const bool reinitChi2 = true,
+                  const FitVerbosity verbosity = FitVerbosity::Silent);
+    FitResult fit(const DoubleModel &model, Minimizer &minimizer,
+                  const DVec &init, const bool reinitChi2 = true,
+                  const FitVerbosity verbosity = FitVerbosity::Silent);
+private:
+    DMat                                   x_, y_;
+    Mat<DMat>                              var_[3];
+    IVec                                   isXExact_, isFitPoint_;
+    Chi2Function                           chi2_;
+};
+
+/******************************************************************************
+ *                      XYStatData template implementation                    *
+ ******************************************************************************/
+
+END_NAMESPACE
+
+#endif // Latan_XYData_hpp_

utils/sample_read.cpp

@@ -1,5 +1,6 @@
 #include <iostream>
 #include <latan/AsciiFile.hpp>
+#include <latan/Plot.hpp>
 
 using namespace std;
 using namespace Latan;