FitInterface data layout complete, first version of the new XYStatData

examples/exFit.cpp

@@ -1,35 +1,46 @@
 #include <iostream>
-#include <LatAnalyze/FitInterface.hpp>
+#include <LatAnalyze/XYStatData.hpp>
 
 using namespace std;
 using namespace Latan;
 
 int main(void)
 {
-    FitInterface f;
+    XYStatData f;
     
     f.addYDim("q1");
     f.addYDim("q2");
     f.addXDim("x1", 6);
     f.addXDim("x2", 5);
     f.addXDim("x3", 5);
-    f.registerDataPoint(f.dataIndex(0,0,0), 0);
+    f.y(f.dataIndex(0,0,0), 0) = 2;
-    f.registerDataPoint(f.dataIndex(1,1,1), 0);
+    f.y(f.dataIndex(1,1,1), 0) = 4;
-    f.registerDataPoint(f.dataIndex(2,2,2), 0);
+    f.y(f.dataIndex(2,2,2), 0) = 5;
-    f.registerDataPoint(f.dataIndex(2,3,3), 0);
+    f.y(f.dataIndex(2,3,3), 0) = 10;
-    f.registerDataPoint(f.dataIndex(0,0,0), 1);
+    f.y(f.dataIndex(0,0,0), 1) = 1;
-    f.registerDataPoint(f.dataIndex(1,1,1), 1);
+    f.y(f.dataIndex(1,1,1), 1) = 2;
-    f.registerDataPoint(f.dataIndex(2,2,3), 1);
+    f.y(f.dataIndex(2,2,3), 1) = 4;
+    f.fitPoint(false, f.dataIndex(2,2,2), 0);
     f.fitPoint(false, f.dataIndex(1,1,1), 1);
-    f.assumeXXCorrelated(true, 0, 0, 0, 1);
+    f.assumeXXCorrelated(true, 0, 0, 1, 0);
-    f.assumeXXCorrelated(true, 1, 1, 0, 1);
+    f.assumeXXCorrelated(true, 0, 1, 1, 1);
-    f.assumeXXCorrelated(true, 2, 2, 0, 1);
+    f.assumeXXCorrelated(true, 0, 2, 1, 2);
-    f.assumeYYCorrelated(true, 0, 0, f.dataIndex(0,0,0), f.dataIndex(1,1,1));
+    f.assumeXXCorrelated(true, 0, 0, 1, 2);
-    f.assumeYYCorrelated(true, 1, 1, f.dataIndex(0,0,0), f.dataIndex(2,2,3));
+    f.assumeXXCorrelated(true, 3, 2, 4, 2);
-    f.assumeXYCorrelated(true, 0, 0, 0, f.dataIndex(1,1,1));
+    f.assumeYYCorrelated(true, f.dataIndex(0,0,0), 0, f.dataIndex(2,3,3), 0);
+    f.assumeYYCorrelated(true, f.dataIndex(0,0,0), 1, f.dataIndex(2,2,3), 1);
+    f.assumeXYCorrelated(true, 0, 0, f.dataIndex(1,1,1), 0);
+    f.assumeXExact(true, 0);
+    f.assumeXExact(true, 1);
+    f.assumeXExact(true, 2);
     cout << f << endl;
-    f.updateLayout();
+    f.setXXVar(0, 0, DMat::Identity(6, 6));
+    f.setXXVar(0, 2, DMat::Identity(6, 5));
+    f.setXXVar(1, 1, DMat::Identity(5, 5));
+    f.setXXVar(2, 2, DMat::Identity(5, 5));
     DEBUG_MAT(f.makeCorrFilter());
+    DEBUG_MAT(f.getFitVar());
+    DEBUG_MAT(f.getFitVar().cwiseProduct(f.makeCorrFilter()));
     
     return EXIT_SUCCESS;
 }

lib/FitInterface.cpp

@@ -23,50 +23,6 @@
 using namespace std;
 using namespace Latan;
 
-// error checks ////////////////////////////////////////////////////////////////
-#define checkXDim(i)\
-if ((i) >= getNXDim())\
-{\
-    LATAN_ERROR(Range, "X dimension " + strFrom(i) + " out of range");\
-}
-
-#define checkXIndex(vi, i)\
-if ((vi) >= getXSize(i))\
-{\
-    LATAN_ERROR(Range, "index " + strFrom(vi) + " in X dimension "\
-                + strFrom(i) + " out of range");\
-}
-
-#define checkYDim(j)\
-if ((j) >= getNYDim())\
-{\
-    LATAN_ERROR(Range, "Y dimension " + strFrom(j) + " out of range");\
-}
-
-#define checkDataIndex(k)\
-if ((k) >= maxDataIndex_)\
-{\
-    LATAN_ERROR(Range, "data point index " + strFrom(k) + " invalid");\
-}
-
-#define checkDataCoord(v)\
-if (static_cast<Index>((v).size()) != getNXDim())\
-{\
-    LATAN_ERROR(Size, "number of coordinates and number of X dimensions "\
-                "mismatch");\
-}\
-for (unsigned int i_ = 0; i_ < (v).size(); ++i_)\
-{\
-    checkXIndex((v)[i_], i_);\
-}
-
-#define checkPoint(k, j)\
-if (!isXUsed(k, j))\
-{\
-    LATAN_ERROR(Range, "no data point in Y dimension " + strFrom(j)\
-                + " with index " + strFrom(k));\
-}
-
 /******************************************************************************
  *                     FitInterface implementation                            *
  ******************************************************************************/
@@ -90,7 +46,8 @@ void FitInterface::addXDim(const string name, const Index nData,
         xIsExact_.push_back(isExact);
         xDimIndex_[name]  = xDimName_.size();
         maxDataIndex_    *= nData;
-        updateDataSize();
+        createXData(nData);
+        scheduleLayoutInit();
     }
 }
 
@@ -99,6 +56,8 @@ void FitInterface::addYDim(const string name)
     yDimName_.push_back(name);
     yDataIndex_.push_back(map<Index, bool>());
     yDimIndex_[name] = yDimName_.size();
+    createYData();
+    scheduleLayoutInit();
 }
 
 // size access /////////////////////////////////////////////////////////////////
@@ -209,6 +168,11 @@ Index FitInterface::getYFitSize(const Index j) const
     return size;
 }
 
+Index FitInterface::getMaxDataIndex(void) const
+{
+    return maxDataIndex_;
+}
+
 // Y dimension index helper ////////////////////////////////////////////////////
 Index FitInterface::dataIndex(const vector<Index> &v) const
 {
@@ -248,6 +212,7 @@ void FitInterface::fitPoint(const bool isFitPoint, const Index k, const Index j)
 {
     checkPoint(k, j);
     yDataIndex_[j][k] = isFitPoint;
+    scheduleLayoutInit();
 }
 
 // variance interface //////////////////////////////////////////////////////////
@@ -255,6 +220,7 @@ void FitInterface::assumeXExact(const bool isExact, const Index i)
 {
     checkXDim(i);
     xIsExact_[i] = isExact;
+    scheduleLayoutInit();
 }
 
 void FitInterface::addCorr(set<array<Index, 4>> &s, const bool isCorr,
@@ -275,25 +241,26 @@ void FitInterface::addCorr(set<array<Index, 4>> &s, const bool isCorr,
     }
 }
 
-void FitInterface::assumeXXCorrelated(const bool isCorr, const Index i1,
+void FitInterface::assumeXXCorrelated(const bool isCorr, const Index r1,
-                                      const Index i2, const Index vi1,
+                                      const Index i1, const Index r2,
-                                      const Index vi2)
+                                      const Index i2)
 {
-    array<Index, 4> c{{i1, i2, vi1, vi2}};
+    array<Index, 4> c{{r1, i1, r2, i2}};
     
-    checkXIndex(vi1, i1);
+    checkXIndex(r1, i1);
-    checkXIndex(vi2, i2);
+    checkXIndex(r2, i2);
-    if ((i1 != i2) or (vi1 != vi2))
+    if ((i1 != i2) or (r1 != r2))
     {
         addCorr(xxCorr_, isCorr, c);
     }
+    scheduleLayoutInit();
 }
 
-void FitInterface::assumeYYCorrelated(const bool isCorr, const Index j1,
+void FitInterface::assumeYYCorrelated(const bool isCorr, const Index k1,
-                                      const Index j2, const Index k1,
+                                      const Index j1, const Index k2,
-                                      const Index k2)
+                                      const Index j2)
 {
-    array<Index, 4> c{{j1, j2, k1, k2}};
+    array<Index, 4> c{{k1, j1, k2, j2}};
     
     checkPoint(k1, j1);
     checkPoint(k2, j2);
@@ -301,39 +268,64 @@ void FitInterface::assumeYYCorrelated(const bool isCorr, const Index j1,
     {
         addCorr(yyCorr_, isCorr, c);
     }
+    scheduleLayoutInit();
 }
 
-void FitInterface::assumeXYCorrelated(const bool isCorr, const Index i,
+void FitInterface::assumeXYCorrelated(const bool isCorr, const Index r,
-                                      const Index j, const Index vi,
+                                      const Index i, const Index k,
-                                      const Index k)
+                                      const Index j)
 {
-    array<Index, 4> c{{i, j, vi, k}};
+    array<Index, 4> c{{r, i, k, j}};
     
-    checkXIndex(vi, i);
+    checkXIndex(r, i);
     checkPoint(k, j);
     addCorr(xyCorr_, isCorr, c);
+    scheduleLayoutInit();
 }
 
 // tests ///////////////////////////////////////////////////////////////////////
-bool FitInterface::isXUsed(const Index k) const
+bool FitInterface::pointExists(const Index k) const
 {
     bool isUsed = false;
     
     for (Index j = 0; j < getNYDim(); ++j)
     {
-        isUsed = isUsed or isXUsed(k, j);
+        isUsed = isUsed or pointExists(k, j);
     }
     
     return isUsed;
 }
 
-bool FitInterface::isXUsed(const Index k, const Index j) const
+bool FitInterface::pointExists(const Index k, const Index j) const
 {
     checkYDim(j);
     
     return !(yDataIndex_[j].find(k) == yDataIndex_[j].end());
 }
 
+bool FitInterface::isXUsed(const Index r, const Index i, const bool inFit) const
+{
+    vector<Index> v;
+    
+    checkXDim(i);
+    for (Index j = 0; j < getNYDim(); ++j)
+    {
+        for (auto &p: yDataIndex_[j])
+        {
+            if (p.second or !inFit)
+            {
+                v      = dataCoord(p.first);
+                if (v[i] == r)
+                {
+                    return true;
+                }
+            }
+        }
+    }
+    
+    return false;
+}
+
 bool FitInterface::isFitPoint(const Index k, const Index j) const
 {
     checkPoint(k, j);
@@ -341,115 +333,192 @@ bool FitInterface::isFitPoint(const Index k, const Index j) const
     return yDataIndex_[j].at(k);
 }
 
+// make correlation filter for fit variance matrix /////////////////////////////
+DMat FitInterface::makeCorrFilter(void)
+{
+    updateLayout();
+    
+    DMat  f = DMat::Identity(layout.totalSize, layout.totalSize);
+    Index row, col;
+    
+    for (auto &c: xxCorr_)
+    {
+        row = indX(c[0], c[1]);
+        col = indX(c[2], c[3]);
+        if ((row != -1) and (col != -1))
+        {
+            f(row, col) = 1.;
+            f(col, row) = 1.;
+        }
+    }
+    for (auto &c: yyCorr_)
+    {
+        row = indY(c[0], c[1]);
+        col = indY(c[2], c[3]);
+        if ((row != -1) and (col != -1))
+        {
+            f(row, col) = 1.;
+            f(col, row) = 1.;
+        }
+    }
+    for (auto &c: xyCorr_)
+    {
+        row = indX(c[0], c[1]);
+        col = indY(c[2], c[3]);
+        if ((row != -1) and (col != -1))
+        {
+            f(row, col) = 1.;
+            f(col, row) = 1.;
+        }
+    }
+    
+    return f;
+}
+
 // register a data point ///////////////////////////////////////////////////////
 void FitInterface::registerDataPoint(const Index k, const Index j)
 {
     checkYDim(j);
     yDataIndex_[j][k] = true;
+    scheduleLayoutInit();
 }
 
 // global layout management ////////////////////////////////////////////////////
+void FitInterface::scheduleLayoutInit(void)
+{
+    initLayout_ = true;
+}
+
 void FitInterface::updateLayout(void)
 {
-    Index size;
+    if (initLayout_)
+    {
+        Index size, ifit;
         
-    layout_.totalSize  = 0;
+        layout.totalSize  = 0;
-    layout_.totalXSize = 0;
+        layout.totalXSize = 0;
-    layout_.totalYSize = 0;
+        layout.totalYSize = 0;
-    layout_.xSize.clear();
+        layout.xSize.clear();
-    layout_.ySize.clear();
+        layout.ySize.clear();
-    layout_.dataIndex.clear();
+        layout.xDim.clear();
-    layout_.xTrans.clear();
+        layout.yDim.clear();
-    layout_.dataTrans.clear();
+        layout.xFitDim.clear();
-    for (Index i = 0; i < getNXDim(); ++i)
+        layout.yFitDim.clear();
-    {
+        layout.x.clear();
-        if (!xIsExact_[i])
+        layout.y.clear();
+        layout.xFit.clear();
+        layout.yFit.clear();
+        ifit = 0;
+        for (Index i = 0; i < getNXDim(); ++i)
         {
-            size = getXFitSize(i);
+            if (!xIsExact_[i])
-            layout_.xSize.push_back(size);
-            layout_.totalXSize += size;
-            layout_.xTrans[i] = static_cast<Index>(layout_.xSize.size() - 1);
-        }
-    }
-    for (Index j = 0; j < getNYDim(); ++j)
-    {
-        size = getYFitSize(j);
-        layout_.ySize.push_back(size);
-        layout_.totalYSize += size;
-    }
-    layout_.totalSize = layout_.totalXSize + layout_.totalYSize;
-    layout_.dataIndex.resize(layout_.ySize.size());
-    for (Index j = 0; j < getNYDim(); ++j)
-    {
-        for (auto &p: yDataIndex_[j])
-        {
-            if (p.second)
             {
-                layout_.dataIndex[j].push_back(p.first);
+                size = getXFitSize(i);
-                layout_.dataTrans[p.first] = layout_.dataIndex[j].size() - 1;
+                layout.xSize.push_back(size);
+                layout.totalXSize += size;
+                layout.xDim.push_back(i);
+                layout.xFitDim.push_back(layout.xDim.size() - 1);
+                layout.x.push_back(vector<Index>());
+                layout.xFit.push_back(vector<Index>());
+                for (Index r = 0; r < getXSize(i); ++r)
+                {
+                    if (isXUsed(r, i))
+                    {
+                        layout.x[ifit].push_back(r);
+                        layout.xFit[i].push_back(layout.x[ifit].size() - 1);
+                    }
+                    else
+                    {
+                        layout.xFit[i].push_back(-1);
+                    }
+                }
+                ifit++;
+            }
+            else
+            {
+                layout.xFitDim.push_back(-1);
+                layout.xFit.push_back(vector<Index>());
+                for (Index r = 0; r < getXSize(i); ++r)
+                {
+                    layout.xFit[i].push_back(-1);
+                }
             }
         }
+        for (Index j = 0; j < getNYDim(); ++j)
+        {
+            Index s = 0;
+            
+            size = getYFitSize(j);
+            layout.ySize.push_back(size);
+            layout.totalYSize += size;
+            layout.yDim.push_back(j);
+            layout.yFitDim.push_back(layout.yDim.size() - 1);
+            layout.y.push_back(vector<Index>());
+            layout.yFit.push_back(vector<Index>());
+            layout.data.push_back(vector<Index>());
+            layout.yFitFromData.push_back(map<Index, Index>());
+            for (auto &p: yDataIndex_[j])
+            {
+                if (p.second)
+                {
+                    layout.y[j].push_back(s);
+                    layout.yFit[j].push_back(layout.y[j].size() - 1);
+                    layout.data[j].push_back(p.first);
+                    layout.yFitFromData[j][p.first] = layout.y[j].size() - 1;
+                }
+                else
+                {
+                    layout.yFit[j].push_back(-1);
+                    layout.yFitFromData[j][p.first] = -1;
+                }
+                s++;
+            }
+        }
+        layout.totalSize = layout.totalXSize + layout.totalYSize;
+        layout.nXFitDim  = static_cast<Index>(layout.xSize.size());
+        layout.nYFitDim  = static_cast<Index>(layout.ySize.size());
+        initLayout_      = false;
     }
 }
 
-// WARNING: NO INDEX CHECKS IN INDEXING FUNCTIONS
+Index FitInterface::indX(const Index r, const Index i) const
-Index FitInterface::indX(const Index vi, const Index i) const
 {
-    Index ind, iTrans;
+    Index ind = -1;
     
-    iTrans = layout_.xTrans.at(i);
+    if (layout.xFit[i][r] != -1)
-    ind    = layout_.totalYSize;
-    for (Index a = 0; a < iTrans; ++a)
     {
-        ind += layout_.xSize[a];
+        Index ifit = layout.xFitDim[i], rfit = layout.xFit[i][r];
+        
+        ind = layout.totalYSize;
+        for (Index a = 0; a < ifit; ++a)
+        {
+            ind += layout.xSize[a];
+        }
+        ind += rfit;
     }
-    ind += vi;
     
     return ind;
 }
 
 Index FitInterface::indY(const Index k, const Index j) const
 {
-    Index ind = 0;
+    Index ind = -1;
     
-    for (Index a = 0; a < j; ++a)
+    if (layout.yFitFromData[j].at(k) != -1)
     {
-        ind += layout_.ySize[a];
+        Index jfit = layout.yFitDim[j], sfit = layout.yFitFromData[j].at(k);
+        
+        ind = 0;
+        for (Index b = 0; b < jfit; ++b)
+        {
+            ind += layout.ySize[b];
+        }
+        ind += sfit;
     }
-    ind += layout_.dataTrans.at(k);
     
     return ind;
 }
 
-DMat FitInterface::makeCorrFilter(void) const
-{
-    DMat f = DMat::Identity(layout_.totalSize, layout_.totalSize);
-    
-    for (auto &c: xxCorr_)
-    {
-        f(indX(c[2], c[0]), indX(c[3], c[1])) = 1.;
-        f(indX(c[3], c[1]), indX(c[2], c[0])) = 1.;
-    }
-    for (auto &c: yyCorr_)
-    {
-        if (isFitPoint(c[2], c[0]) and isFitPoint(c[3], c[1]))
-        {
-            f(indY(c[2], c[0]), indY(c[3], c[1])) = 1.;
-            f(indY(c[3], c[1]), indY(c[2], c[0])) = 1.;
-        }
-    }
-    for (auto &c: xyCorr_)
-    {
-        if (isFitPoint(c[3], c[1]))
-        {
-            f(indX(c[2], c[0]), indY(c[3], c[1])) = 1.;
-            f(indY(c[3], c[1]), indX(c[2], c[0])) = 1.;
-        }
-    }
-    
-    return f;
-}
-
 // IO //////////////////////////////////////////////////////////////////////////
 ostream & Latan::operator<<(ostream &out, FitInterface &f)
 {
@@ -457,7 +526,12 @@ ostream & Latan::operator<<(ostream &out, FitInterface &f)
     for (Index i = 0; i < f.getNXDim(); ++i)
     {
         out << "  * " << i << " \"" << f.xDimName_[i] << "\": ";
-        out << f.getXSize(i) << " value(s)" << endl;
+        out << f.getXSize(i) << " value(s)";
+        if (f.xIsExact_[i])
+        {
+            out << " (assumed exact)";
+        }
+        out << endl;
     }
     out << "Y dimensions: " << f.getNYDim() << endl;
     for (Index j = 0; j < f.getNYDim(); ++j)
@@ -474,7 +548,7 @@ ostream & Latan::operator<<(ostream &out, FitInterface &f)
             out << "\b) fit: " << (p.second ? "true" : "false") << endl;
         }
     }
-    out << "X/X correlations (i1 i2 vi1 vi2): ";
+    out << "X/X correlations (r1 i1 r2 i2): ";
     if (f.xxCorr_.empty())
     {
         out << "no" << endl;
@@ -492,7 +566,7 @@ ostream & Latan::operator<<(ostream &out, FitInterface &f)
             out << endl;
         }
     }
-    out << "Y/Y correlations (j1 j2 k1 k2): ";
+    out << "Y/Y correlations (k1 j1 k2 j2): ";
     if (f.yyCorr_.empty())
     {
         out << "no" << endl;
@@ -510,7 +584,7 @@ ostream & Latan::operator<<(ostream &out, FitInterface &f)
             out << endl;
         }
     }
-    out << "X/Y correlations (i j vi k): ";
+    out << "X/Y correlations (r i k j): ";
     if (f.xyCorr_.empty())
     {
         out << "no";

lib/FitInterface.hpp

@@ -33,10 +33,21 @@ class FitInterface
 private:
     typedef struct
     {
+        Index                           nXFitDim, nYFitDim;
+        // X/Y block sizes
         Index                           totalSize, totalXSize, totalYSize;
+        // size of each X/Y dimension
         std::vector<Index>              xSize, ySize;
-        std::vector<std::vector<Index>> dataIndex;
+        // lookup tables
-        std::map<Index, Index>          xTrans, dataTrans;
+        // xDim        : x fit dim ifit -> x dim i
+        // x           : x fit point ifit,rfit -> x point r
+        // xFitDim     : x dim i -> x fit dim ifit (-1 if empty)
+        // xFit        : x point i,r -> x fit point rfit (-1 if empty)
+        // data        : y fit point jfit,sfit -> y point index k
+        // yFitFromData: y point indec k,j -> y fit point sfit (-1 if empty)
+        std::vector<Index>                  xDim, yDim, xFitDim, yFitDim;
+        std::vector<std::vector<Index>>     x, y, data, xFit, yFit;
+        std::vector<std::map<Index, Index>> yFitFromData;
     } Layout;
 public:
     // constructor
@@ -58,6 +69,7 @@ public:
     Index getXFitSize(const Index i) const;
     Index getYFitSize(void) const;
     Index getYFitSize(const Index j) const;
+    Index getMaxDataIndex(void) const;
     // Y dimension index helper
     template <typename... Ts>
     Index              dataIndex(const Ts... is) const;
@@ -67,32 +79,37 @@ public:
     void fitPoint(const bool isFitPoint, const Index k, const Index j = 0);
     // variance interface
     void assumeXExact(const bool isExact, const Index i);
-    void assumeXXCorrelated(const bool isCorr, const Index i1, const Index i2,
+    void assumeXXCorrelated(const bool isCorr, const Index r1, const Index i1,
-                            const Index vi1, const Index vi2);
+                            const Index r2, const Index i2);
-    void assumeYYCorrelated(const bool isCorr, const Index j1, const Index j2,
+    void assumeYYCorrelated(const bool isCorr, const Index k1, const Index j1,
-                            const Index k1, const Index k2);
+                            const Index k2, const Index j2);
-    void assumeXYCorrelated(const bool isCorr, const Index i, const Index j,
+    void assumeXYCorrelated(const bool isCorr, const Index r, const Index i,
-                            const Index vi, const Index k);
+                            const Index k, const Index j);
     // tests
-    bool isXUsed(const Index k) const;
+    bool pointExists(const Index k) const;
-    bool isXUsed(const Index k, const Index j) const;
+    bool pointExists(const Index k, const Index j) const;
+    bool isXUsed(const Index r, const Index i, const bool inFit = true) const;
     bool isFitPoint(const Index k, const Index j) const;
+    // make correlation filter for fit variance matrix
+    DMat makeCorrFilter(void);
     // IO
     friend std::ostream & operator<<(std::ostream &out, FitInterface &f);
 protected:
-public:
     // register a data point
     void registerDataPoint(const Index k, const Index j = 0);
     // add correlation to a set
     static void addCorr(std::set<std::array<Index, 4>> &s, const bool isCorr,
                         const std::array<Index, 4> &c);
-    // abstract method to update data container size
+    // abstract methods to create data containers
-    virtual void updateDataSize(void) {};
+    virtual void createXData(const Index nData) = 0;
+    virtual void createYData(void) = 0;
     // global layout management
+    void  scheduleLayoutInit(void);
     void  updateLayout(void);
-    Index indX(const Index vi, const Index i) const;
+    Index indX(const Index r, const Index i) const;
     Index indY(const Index k, const Index j) const;
-    DMat  makeCorrFilter(void) const;
+protected:
+    Layout layout;
 private:
     std::vector<std::string>           xDimName_, yDimName_;
     std::map<std::string, Index>       xDimIndex_, yDimIndex_;
@@ -101,7 +118,7 @@ private:
     std::vector<std::map<Index, bool>> yDataIndex_;
     std::set<std::array<Index, 4>>     xxCorr_, yyCorr_, xyCorr_;
     Index                              maxDataIndex_{1};
-    Layout                             layout_;
+    bool                               initLayout_{true};
 };
 
 std::ostream & operator<<(std::ostream &out, FitInterface &f);
@@ -121,6 +138,52 @@ Index FitInterface::dataIndex(const Ts... coords) const
     return dataIndex(coord);
 }
 
+/******************************************************************************
+ *                       error check macros                                   *
+ ******************************************************************************/
+#define checkXDim(i)\
+if ((i) >= getNXDim())\
+{\
+    LATAN_ERROR(Range, "X dimension " + strFrom(i) + " out of range");\
+}
+
+#define checkXIndex(vi, i)\
+if ((vi) >= getXSize(i))\
+{\
+    LATAN_ERROR(Range, "index " + strFrom(vi) + " in X dimension "\
+                + strFrom(i) + " out of range");\
+}
+
+#define checkYDim(j)\
+if ((j) >= getNYDim())\
+{\
+    LATAN_ERROR(Range, "Y dimension " + strFrom(j) + " out of range");\
+}
+
+#define checkDataIndex(k)\
+if ((k) >= getMaxDataIndex())\
+{\
+    LATAN_ERROR(Range, "data point index " + strFrom(k) + " invalid");\
+}
+
+#define checkDataCoord(v)\
+if (static_cast<Index>((v).size()) != getNXDim())\
+{\
+    LATAN_ERROR(Size, "number of coordinates and number of X dimensions "\
+                "mismatch");\
+}\
+for (unsigned int i_ = 0; i_ < (v).size(); ++i_)\
+{\
+    checkXIndex((v)[i_], i_);\
+}
+
+#define checkPoint(k, j)\
+if (!pointExists(k, j))\
+{\
+    LATAN_ERROR(Range, "no data point in Y dimension " + strFrom(j)\
+                + " with index " + strFrom(k));\
+}
+
 END_LATAN_NAMESPACE
 
 #endif // Latan_FitInterface_hpp_

lib/Makefile.am

@@ -45,6 +45,7 @@ libLatAnalyze_la_SOURCES = \
     RootFinder.cpp         \
     Solver.cpp             \
     TabFunction.cpp        \
+    XYStatData.cpp         \
 	../config.h
 libLatAnalyze_ladir = $(pkgincludedir)
 libLatAnalyze_la_HEADERS = \
@@ -77,7 +78,8 @@ libLatAnalyze_la_HEADERS = \
     RootFinder.hpp         \
     TabFunction.hpp        \
     Solver.hpp             \
-    StatArray.hpp
+    StatArray.hpp          \
+    XYStatData.hpp
 if HAVE_MINUIT
     libLatAnalyze_la_SOURCES += MinuitMinimizer.cpp
     libLatAnalyze_la_HEADERS += MinuitMinimizer.hpp

lib/XYStatData.cpp

@@ -0,0 +1,306 @@
+/*
+ * XYStatData.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2016 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/XYStatData.hpp>
+#include <LatAnalyze/includes.hpp>
+
+using namespace std;
+using namespace Latan;
+
+/******************************************************************************
+ *                       XYStatData implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+XYStatData::XYStatData(void)
+{}
+
+// schedule fit var matrix update //////////////////////////////////////////////
+void XYStatData::scheduleFitVarMatInit(void)
+{
+    initVarMat_ = true;
+}
+
+// create data /////////////////////////////////////////////////////////////////
+void XYStatData::createXData(const Index nData)
+{
+    xData_.push_back(DVec::Zero(nData));
+    resizeVarMat();
+}
+
+void XYStatData::createYData(void)
+{
+    yData_.push_back(map<Index, double>());
+    resizeVarMat();
+}
+
+void XYStatData::resizeVarMat(void)
+{
+    xxVar_.conservativeResize(getNXDim(), getNXDim());
+    for (Index i1 = 0; i1 < getNXDim(); ++i1)
+    for (Index i2 = 0; i2 < getNXDim(); ++i2)
+    {
+        xxVar_(i1, i2).conservativeResize(getXSize(i1), getXSize(i2));
+    }
+    yyVar_.conservativeResize(getNYDim(), getNYDim());
+    for (Index j1 = 0; j1 < getNYDim(); ++j1)
+    for (Index j2 = 0; j2 < getNYDim(); ++j2)
+    {
+        yyVar_(j1, j2).conservativeResize(getYSize(j1), getYSize(j2));
+    }
+    xyVar_.conservativeResize(getNXDim(), getNYDim());
+    for (Index i = 0; i < getNXDim(); ++i)
+    for (Index j = 0; j < getNYDim(); ++j)
+    {
+        xyVar_(i, j).conservativeResize(getXSize(i), getYSize(j));
+    }
+}
+
+// data access /////////////////////////////////////////////////////////////////
+double & XYStatData::x(const Index vi, const Index i)
+{
+    checkXIndex(vi, i);
+    scheduleFitVarMatInit();
+    
+    return xData_[i](vi);
+}
+
+const double & XYStatData::x(const Index vi, const Index i) const
+{
+    checkXIndex(vi, i);
+    
+    return xData_[i](vi);
+}
+
+double & XYStatData::y(const Index k, const Index j)
+{
+    checkYDim(j);
+    if (!pointExists(k, j))
+    {
+        registerDataPoint(k, j);
+        resizeVarMat();
+        scheduleFitVarMatInit();
+    }
+    
+    return yData_[j][k];
+}
+
+const double & XYStatData::y(const Index k, const Index j) const
+{
+    checkPoint(k, j);
+    
+    return yData_[j].at(k);
+}
+
+#define checkVarMat(m, var)\
+if (((m).rows() != (var).rows()) or ((m).cols() != (var).cols()))\
+{\
+    LATAN_ERROR(Size, "provided variance matrix has a wrong size"\
+                " (expected " + strFrom((var).rows()) + "x"\
+                + strFrom((var).cols()) + ", got " + strFrom((m).rows())\
+                + "x" + strFrom((m).cols()) + ")");\
+}
+
+void XYStatData::setXXVar(const Index i1, const Index i2, const DMat &m)
+{
+    checkXDim(i1);
+    checkXDim(i2);
+    checkVarMat(m, xxVar_(i1, i2));
+    xxVar_(i1, i2) = m;
+    if (i1 != i2)
+    {
+        xxVar_(i2, i1) = m.transpose();
+    }
+    scheduleFitVarMatInit();
+}
+
+void XYStatData::setYYVar(const Index j1, const Index j2, const DMat &m)
+{
+    checkYDim(j1);
+    checkYDim(j2);
+    checkVarMat(m, yyVar_(j1, j2));
+    yyVar_(j1, j2) = m;
+    if (j1 != j2)
+    {
+        yyVar_(j2, j1) = m.transpose();
+    }
+    scheduleFitVarMatInit();
+}
+
+void XYStatData::setXYVar(const Index i, const Index j, const DMat &m)
+{
+    checkXDim(i);
+    checkYDim(j);
+    checkVarMat(m, xyVar_(i, j));
+    xyVar_(i, j) = m;
+    scheduleFitVarMatInit();
+}
+
+#define checkErrVec(err, var)\
+if ((err).size() != (var).rows())\
+{\
+    LATAN_ERROR(Size, "provided error vector has a wrong size"\
+                " (expected " + strFrom((var).rows()) + ", got "\
+                + strFrom((err).size()) + ")");\
+}
+
+void XYStatData::setXError(const Index i, const DVec &err)
+{
+    checkXDim(i);
+    checkErrVec(err, xxVar_(i, i));
+    xxVar_(i, i).diagonal() = err.cwiseProduct(err);
+    scheduleFitVarMatInit();
+}
+
+void XYStatData::setYError(const Index j, const DVec &err)
+{
+    checkXDim(j);
+    checkErrVec(err, yyVar_(j, j));
+    yyVar_(j, j).diagonal() = err.cwiseProduct(err);
+    scheduleFitVarMatInit();
+}
+
+const DMat & XYStatData::getXXVar(const Index i1, const Index i2) const
+{
+    checkXDim(i1);
+    checkXDim(i2);
+    
+    return xxVar_(i1, i2);
+}
+
+const DMat & XYStatData::getYYVar(const Index j1, const Index j2) const
+{
+    checkYDim(j1);
+    checkYDim(j2);
+    
+    return yyVar_(j1, j2);
+}
+
+const DMat & XYStatData::getXYVar(const Index i, const Index j) const
+{
+    checkXDim(i);
+    checkYDim(j);
+    
+    return xyVar_(i, j);
+}
+
+DVec XYStatData::getXError(const Index i) const
+{
+    checkXDim(i);
+    
+    return xxVar_(i, i).diagonal().cwiseSqrt();
+}
+
+DVec XYStatData::getYError(const Index j) const
+{
+    checkXDim(j);
+    
+    return yyVar_(j, j).diagonal().cwiseSqrt();
+}
+
+// get total fit variance matrix ///////////////////////////////////////////////
+const DMat & XYStatData::getFitVar(void)
+{
+    updateFitVarMat();
+    
+    return fitVar_;
+}
+
+// make total fit variance matrix //////////////////////////////////////////////
+void XYStatData::updateFitVarMat(void)
+{
+    if (initVarMat_)
+    {
+        updateLayout();
+        
+        DMat  &v = fitVar_;
+        Index roffs, coffs;
+        
+        v.resize(layout.totalSize, layout.totalSize);
+        roffs = layout.totalYSize;
+        for (Index ifit1 = 0; ifit1 < layout.nXFitDim; ++ifit1)
+        {
+            coffs = layout.totalYSize;
+            for (Index ifit2 = 0; ifit2 < layout.nXFitDim; ++ifit2)
+            {
+                for (Index rfit1 = 0; rfit1 < layout.xSize[ifit1]; ++rfit1)
+                for (Index rfit2 = 0; rfit2 < layout.xSize[ifit2]; ++rfit2)
+                {
+                    Index i1, i2, r1, r2;
+                    
+                    i1 = layout.xDim[ifit1];
+                    i2 = layout.xDim[ifit2];
+                    r1 = layout.x[ifit1][rfit1];
+                    r2 = layout.x[ifit2][rfit2];
+                    
+                    v(roffs+rfit1, coffs+rfit2) = xxVar_(i1, i2)(r1, r2);
+                    v(coffs+rfit2, roffs+rfit1) = v(roffs+rfit1, coffs+rfit2);
+                }
+                coffs += layout.xSize[ifit2];
+            }
+            roffs += layout.xSize[ifit1];
+        }
+        roffs = 0;
+        for (Index jfit1 = 0; jfit1 < layout.nYFitDim; ++jfit1)
+        {
+            coffs = 0;
+            for (Index jfit2 = 0; jfit2 < layout.nYFitDim; ++jfit2)
+            {
+                for (Index sfit1 = 0; sfit1 < layout.ySize[jfit1]; ++sfit1)
+                for (Index sfit2 = 0; sfit2 < layout.ySize[jfit2]; ++sfit2)
+                {
+                    Index j1, j2, s1, s2;
+                    
+                    j1 = layout.yDim[jfit1];
+                    j2 = layout.yDim[jfit2];
+                    s1 = layout.y[jfit1][sfit1];
+                    s2 = layout.y[jfit2][sfit2];
+                    
+                    v(roffs+sfit1, coffs+sfit2) = yyVar_(j1, j2)(s1, s2);
+                    v(coffs+sfit2, roffs+sfit1) = v(roffs+sfit1, coffs+sfit2);
+                }
+                coffs += layout.ySize[jfit2];
+            }
+            roffs += layout.ySize[jfit1];
+        }
+        roffs = layout.totalYSize;
+        for (Index ifit = 0; ifit < layout.nXFitDim; ++ifit)
+        {
+            coffs = 0;
+            for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
+            {
+                for (Index rfit = 0; rfit < layout.xSize[ifit]; ++rfit)
+                for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
+                {
+                    Index i, j, r, s;
+                    
+                    i = layout.xDim[ifit];
+                    j = layout.yDim[jfit];
+                    r = layout.x[ifit][rfit];
+                    s = layout.y[jfit][sfit];
+                    
+                    v(roffs+rfit, coffs+sfit) = xyVar_(i, j)(r, s);
+                    v(coffs+sfit, roffs+rfit) = v(roffs+rfit, coffs+sfit);
+                }
+                coffs += layout.ySize[jfit];
+            }
+            roffs += layout.xSize[ifit];
+        }
+        initVarMat_ = false;
+    }
+}

lib/XYStatData.hpp

@@ -0,0 +1,75 @@
+/*
+ * XYStatData.hpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2016 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef Latan_XYStatData_hpp_
+#define Latan_XYStatData_hpp_
+
+#include <LatAnalyze/Global.hpp>
+#include <LatAnalyze/FitInterface.hpp>
+
+BEGIN_LATAN_NAMESPACE
+
+/******************************************************************************
+ *                    class for X vs. Y statistical data                      *
+ ******************************************************************************/
+class XYStatData: public FitInterface
+{
+public:
+    // constructor
+    XYStatData(void);
+    // destructor
+    virtual ~XYStatData(void) = default;
+    // schedule fit var matrix update
+    void scheduleFitVarMatInit(void);
+    // data access
+    double &       x(const Index vi, const Index i = 0);
+    const double & x(const Index vi, const Index i = 0) const;
+    double &       y(const Index k, const Index j = 0);
+    const double & y(const Index k, const Index j = 0) const;
+    void           setXXVar(const Index i1, const Index i2, const DMat &m);
+    void           setYYVar(const Index j1, const Index j2, const DMat &m);
+    void           setXYVar(const Index i, const Index j, const DMat &m);
+    void           setXError(const Index i, const DVec &err);
+    void           setYError(const Index j, const DVec &err);
+    const DMat &   getXXVar(const Index i1, const Index i2) const;
+    const DMat &   getYYVar(const Index j1, const Index j2) const;
+    const DMat &   getXYVar(const Index i, const Index j) const;
+    DVec           getXError(const Index i) const;
+    DVec           getYError(const Index j) const;
+    // get total fit variance matrix
+    const DMat &   getFitVar(void);
+protected:
+    // create data
+    virtual void createXData(const Index nData);
+    virtual void createYData(void);
+    void         resizeVarMat(void);
+private:
+    // make total fit variance matrix
+    void updateFitVarMat(void);
+private:
+    std::vector<std::map<Index, double>> yData_;
+    std::vector<DVec>                    xData_;
+    Mat<DMat>                            xxVar_, yyVar_, xyVar_;
+    DMat                                 fitVar_;
+    bool                                 initVarMat_{true};
+};
+
+END_LATAN_NAMESPACE
+
+#endif // Latan_XYStatData_hpp_