lot of fixes, XYSampleData fit now seems to work (more checks are probably needed)

examples/Makefile.am

@@ -1,11 +1,3 @@
-if CC_GNU
-    COM_CFLAGS = -Wall -W -pedantic
-else
-if CC_INTEL
-    COM_CFLAGS = -Wall
-endif
-endif
-
 if CXX_GNU
     COM_CXXFLAGS = -Wall -W -pedantic
 else
@@ -26,55 +18,59 @@ noinst_PROGRAMS =           \
     exRootFinder
 
 if HAVE_MINUIT
-noinst_PROGRAMS += exFit exMin
+noinst_PROGRAMS += exFit exFitSample exMin
 endif
 
 exCompiledDoubleFunction_SOURCES  = exCompiledDoubleFunction.cpp
-exCompiledDoubleFunction_CFLAGS  = -g -O2
+exCompiledDoubleFunction_CXXFLAGS = $(COM_CXXFLAGS)
 exCompiledDoubleFunction_LDFLAGS  = -L../lib/.libs -lLatAnalyze
 
 exDerivative_SOURCES              = exDerivative.cpp
-exDerivative_CFLAGS              = -g -O2
+exDerivative_CXXFLAGS             = $(COM_CXXFLAGS)
 exDerivative_LDFLAGS              = -L../lib/.libs -lLatAnalyze
 
 if HAVE_MINUIT
 exFit_SOURCES                     = exFit.cpp
-exFit_CFLAGS                     = -g -O2
+exFit_CXXFLAGS                    = $(COM_CXXFLAGS)
 exFit_LDFLAGS                     = -L../lib/.libs -lLatAnalyze
+
+exFitSample_SOURCES               = exFitSample.cpp
+exFitSample_CXXFLAGS              = $(COM_CXXFLAGS)
+exFitSample_LDFLAGS               = -L../lib/.libs -lLatAnalyze
 endif
 
 exInterp_SOURCES                  = exInterp.cpp
-exInterp_CFLAGS                  = -g -O2
+exInterp_CXXFLAGS                 = $(COM_CXXFLAGS)
 exInterp_LDFLAGS                  = -L../lib/.libs -lLatAnalyze
 
 exIntegrator_SOURCES              = exIntegrator.cpp
-exIntegrator_CFLAGS              = -g -O2
+exIntegrator_CXXFLAGS             = $(COM_CXXFLAGS)
 exIntegrator_LDFLAGS              = -L../lib/.libs -lLatAnalyze
 
 exMat_SOURCES                     = exMat.cpp
-exMat_CFLAGS                     = -g -O2
+exMat_CXXFLAGS                    = $(COM_CXXFLAGS)
 exMat_LDFLAGS                     = -L../lib/.libs -lLatAnalyze
 
 if HAVE_MINUIT
 exMin_SOURCES                     = exMin.cpp
-exMin_CFLAGS                     = -g -O2
+exMin_CXXFLAGS                    = $(COM_CXXFLAGS)
 exMin_LDFLAGS                     = -L../lib/.libs -lLatAnalyze
 endif
 
 exMathInterpreter_SOURCES         = exMathInterpreter.cpp
-exMathInterpreter_CFLAGS         = -g -O2
+exMathInterpreter_CXXFLAGS        = $(COM_CXXFLAGS)
 exMathInterpreter_LDFLAGS         = -L../lib/.libs -lLatAnalyze
 
 exPlot_SOURCES                    = exPlot.cpp
-exPlot_CFLAGS                    = -g -O2
+exPlot_CXXFLAGS                   = $(COM_CXXFLAGS)
 exPlot_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
 
 exRand_SOURCES                    = exRand.cpp
-exRand_CFLAGS                    = -g -O2
+exRand_CXXFLAGS                   = $(COM_CXXFLAGS)
 exRand_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
 
 exRootFinder_SOURCES              = exRootFinder.cpp
-exRootFinder_CFLAGS              = -g -O2
+exRootFinder_CXXFLAGS             = $(COM_CXXFLAGS)
 exRootFinder_LDFLAGS              = -L../lib/.libs -lLatAnalyze
 
 ACLOCAL_AMFLAGS = -I .buildutils/m4

examples/exFit.cpp

@@ -9,7 +9,7 @@ using namespace std;
 using namespace Latan;
 
 const Index  nPoint = 30;
-const double xErr = .01, yErr   = .1;
+const double xErr = .2, yErr   = .1;
 const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast<double>(nPoint);
 
 int main(void)
@@ -17,7 +17,7 @@ int main(void)
     // generate fake data
     XYStatData  data;
     RandGen     rg;
-    double      x_k, y_k;
+    double      x1_k, x2_k;
     DoubleModel f([](const double *x, const double *p)
                   {return p[1]*exp(-x[0]*p[0]);}, 1, 2);
     
@@ -25,11 +25,11 @@ int main(void)
     data.addYDim("y");
     for (Index k = 0; k < nPoint; ++k)
     {
-        x_k          = k*dx + rg.gaussian(0.0, xErr);
+        x1_k      = rg.gaussian(k*dx, xErr);
-        y_k          = f(&x_k, exactPar) + rg.gaussian(0.0, yErr);
+        x2_k      = rg.gaussian(k*dx, xErr);
-        printf("% 8e % 8e % 8e % 8e\n", x_k, xErr, y_k, yErr);
+        data.x(k) = x1_k;
-        data.x(k) = x_k;
+        data.y(k) = rg.gaussian(f(&x2_k, exactPar), yErr);
-        data.y(k) = y_k;
+        printf("% 8e % 8e % 8e % 8e\n", data.x(k), data.y(k), xErr, yErr);
     }
     cout << endl;
     data.setXError(0, DVec::Constant(nPoint, xErr));
@@ -42,9 +42,9 @@ int main(void)
     MinuitMinimizer minimizer;
     
     p = data.fit(minimizer, init, f);
-    cout << "a= " << p(0) << " b= " << p(1);
+    cout << "a= " << p(0) << " b= " << p(1) << endl;
-    cout << " chi^2/ndof= " << p.getChi2PerDof();
+    cout << "chi^2/ndof= " << p.getChi2PerDof() << endl;
-    cout << " p-value= " << p.getPValue() <<endl;
+    cout << "p-value= " << p.getPValue() <<endl;
     
     return EXIT_SUCCESS;
 }

examples/exFitSample.cpp

@@ -0,0 +1,54 @@
+#include <iostream>
+#include <cmath>
+#include <LatAnalyze/CompiledModel.hpp>
+#include <LatAnalyze/MinuitMinimizer.hpp>
+#include <LatAnalyze/RandGen.hpp>
+#include <LatAnalyze/XYSampleData.hpp>
+
+using namespace std;
+using namespace Latan;
+
+const Index  nPoint  = 30;
+const Index  nSample = 1000;
+const double xErr = .2, yErr   = .1;
+const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast<double>(nPoint);
+
+int main(void)
+{
+    // generate fake data
+    DMatSample   x(nSample, nPoint, 1), y(nSample, nPoint, 1);
+    XYSampleData data(nSample);
+    RandGen     rg;
+    double      x1_k, x2_k;
+    DoubleModel f([](const double *t, const double *p)
+                  {return p[1]*exp(-t[0]*p[0]);}, 1, 2);
+    
+    data.addXDim("x", nPoint);
+    data.addYDim("y");
+    FOR_STAT_ARRAY(x, s)
+    {
+        for (Index k = 0; k < nPoint; ++k)
+        {
+            x1_k         = rg.gaussian(k*dx, xErr);
+            x2_k         = rg.gaussian(k*dx, xErr);
+            data.x(k)[s] = x1_k;
+            data.y(k)[s] = rg.gaussian(f(&x2_k, exactPar), yErr);
+        }
+    }
+    cout << data << endl;
+    
+    // fit
+    DVec init = DVec::Constant(2, 0.5);
+    DMat err;
+    SampleFitResult p;
+    MinuitMinimizer minimizer;
+    
+    p   = data.fit(minimizer, init, f);
+    err = p.variance().cwiseSqrt();
+    cout << "a= " << p[central](0) << " +/- " << err(0) << endl;
+    cout << "b= " << p[central](1) << " +/- " << err(1) << endl;
+    cout << "chi^2/ndof= " << p.getChi2PerDof() << endl;
+    cout << "p-value= " << p.getPValue() << endl;
+    
+    return EXIT_SUCCESS;
+}

lib/XYSampleData.cpp

@@ -100,6 +100,10 @@ DSample & XYSampleData::x(const Index r, const Index i)
     checkXIndex(r, i);
     scheduleDataInit();
     scheduleComputeVarMat();
+    if (xData_[i][r].size() == 0)
+    {
+        xData_[i][r].resize(nSample_);
+    }
     
     return xData_[i][r];
 }
@@ -120,6 +124,10 @@ DSample & XYSampleData::y(const Index k, const Index j)
     }
     scheduleDataInit();
     scheduleComputeVarMat();
+    if (yData_[j][k].size() == 0)
+    {
+        yData_[j][k].resize(nSample_);
+    }
     
     return yData_[j][k];
 }
@@ -194,7 +202,6 @@ void XYSampleData::setDataToSample(const Index s)
 {
     if (initData_ or (s != dataSample_))
     {
-        data_.copyInterface(*this);
         for (Index i = 0; i < getNXDim(); ++i)
         for (Index r = 0; r < getXSize(i); ++r)
         {
@@ -205,6 +212,7 @@ void XYSampleData::setDataToSample(const Index s)
         {
             data_.y(p.first, j) = p.second[s];
         }
+        data_.copyInterface(*this);
         dataSample_ = s;
         initData_   = false;
     }
@@ -282,10 +290,11 @@ void XYSampleData::computeVarMat(void)
         // compute total matrix
         DMatSample z(nSample_, size, 1);
         DMat       var;
-        Index      a = 0;
+        Index      a;
         
         FOR_STAT_ARRAY(z, s)
         {
+            a = 0;
             for (Index j = 0; j < getNYDim(); ++j)
             for (auto &p: yData_[j])
             {
@@ -311,7 +320,7 @@ void XYSampleData::computeVarMat(void)
             for (Index i2 = 0; i2 < getNXDim(); ++i2)
             {
                 data_.setXXVar(i1, i2,
-                               var.block(a1, getXSize(i1), a2, getXSize(i2)));
+                               var.block(a1, a2, getXSize(i1), getXSize(i2)));
                 a2 += getXSize(i2);
             }
             a1 += getXSize(i1);
@@ -323,7 +332,7 @@ void XYSampleData::computeVarMat(void)
             for (Index j2 = 0; j2 < getNYDim(); ++j2)
             {
                 data_.setYYVar(j1, j2,
-                               var.block(a1, getYSize(j1), a2, getYSize(j2)));
+                               var.block(a1, a2, getYSize(j1), getYSize(j2)));
                 a2 += getYSize(j2);
             }
             a1 += getYSize(j1);
@@ -332,10 +341,10 @@ void XYSampleData::computeVarMat(void)
         for (Index i = 0; i < getNXDim(); ++i)
         {
             a2 = 0;
-            for (Index j = 0; j < getNXDim(); ++j)
+            for (Index j = 0; j < getNYDim(); ++j)
             {
                 data_.setXYVar(i, j,
-                               var.block(a1, getXSize(i), a2, getYSize(j)));
+                               var.block(a1, a2, getXSize(i), getYSize(j)));
                 a2 += getYSize(j);
             }
             a1 += getXSize(i);
@@ -350,12 +359,14 @@ void XYSampleData::computeVarMat(void)
 }
 
 // create data /////////////////////////////////////////////////////////////////
-void XYSampleData::createXData(const string name __dumb, const Index nData)
+void XYSampleData::createXData(const string name, const Index nData)
 {
+    data_.addXDim(name, nData);
     xData_.push_back(vector<DSample>(nData));
 }
 
-void XYSampleData::createYData(const string name __dumb)
+void XYSampleData::createYData(const string name)
 {
+    data_.addYDim(name);
     yData_.push_back(map<Index, DSample>());
 }

utils/Makefile.am

@@ -1,11 +1,3 @@
-if CC_GNU
-    COM_CFLAGS = -Wall -W -pedantic
-else
-if CC_INTEL
-    COM_CFLAGS = -Wall
-endif
-endif
-
 if CXX_GNU
     COM_CXXFLAGS = -Wall -W -pedantic
 else