From 3793aca27c882d29707c7473a381c6a4d6894c76 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 16 Mar 2016 20:03:32 +0000
Subject: [PATCH] lot of fixes, XYSampleData fit now seems to work (more checks
 are probably needed)

---
 examples/Makefile.am     | 80 +++++++++++++++++++---------------------
 examples/exFit.cpp       | 20 +++++-----
 examples/exFitSample.cpp | 54 +++++++++++++++++++++++++++
 lib/XYSampleData.cpp     | 31 +++++++++++-----
 utils/Makefile.am        |  8 ----
 5 files changed, 123 insertions(+), 70 deletions(-)
 create mode 100644 examples/exFitSample.cpp

diff --git a/examples/Makefile.am b/examples/Makefile.am
index fce581f..0f2fa07 100644
--- a/examples/Makefile.am
+++ b/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_LDFLAGS = -L../lib/.libs -lLatAnalyze
+exCompiledDoubleFunction_SOURCES  = exCompiledDoubleFunction.cpp
+exCompiledDoubleFunction_CXXFLAGS = $(COM_CXXFLAGS)
+exCompiledDoubleFunction_LDFLAGS  = -L../lib/.libs -lLatAnalyze
 
-exDerivative_SOURCES             = exDerivative.cpp
-exDerivative_CFLAGS              = -g -O2
-exDerivative_LDFLAGS             = -L../lib/.libs -lLatAnalyze
+exDerivative_SOURCES              = exDerivative.cpp
+exDerivative_CXXFLAGS             = $(COM_CXXFLAGS)
+exDerivative_LDFLAGS              = -L../lib/.libs -lLatAnalyze
 
 if HAVE_MINUIT
-exFit_SOURCES                    = exFit.cpp
-exFit_CFLAGS                     = -g -O2
-exFit_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
+exFit_SOURCES                     = exFit.cpp
+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_LDFLAGS                 = -L../lib/.libs -lLatAnalyze
+exInterp_SOURCES                  = exInterp.cpp
+exInterp_CXXFLAGS                 = $(COM_CXXFLAGS)
+exInterp_LDFLAGS                  = -L../lib/.libs -lLatAnalyze
 
-exIntegrator_SOURCES             = exIntegrator.cpp
-exIntegrator_CFLAGS              = -g -O2
-exIntegrator_LDFLAGS             = -L../lib/.libs -lLatAnalyze
+exIntegrator_SOURCES              = exIntegrator.cpp
+exIntegrator_CXXFLAGS             = $(COM_CXXFLAGS)
+exIntegrator_LDFLAGS              = -L../lib/.libs -lLatAnalyze
 
-exMat_SOURCES                    = exMat.cpp
-exMat_CFLAGS                     = -g -O2
-exMat_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
+exMat_SOURCES                     = exMat.cpp
+exMat_CXXFLAGS                    = $(COM_CXXFLAGS)
+exMat_LDFLAGS                     = -L../lib/.libs -lLatAnalyze
 
 if HAVE_MINUIT
-exMin_SOURCES                    = exMin.cpp
-exMin_CFLAGS                     = -g -O2
-exMin_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
+exMin_SOURCES                     = exMin.cpp
+exMin_CXXFLAGS                    = $(COM_CXXFLAGS)
+exMin_LDFLAGS                     = -L../lib/.libs -lLatAnalyze
 endif
 
-exMathInterpreter_SOURCES        = exMathInterpreter.cpp
-exMathInterpreter_CFLAGS         = -g -O2
-exMathInterpreter_LDFLAGS        = -L../lib/.libs -lLatAnalyze
+exMathInterpreter_SOURCES         = exMathInterpreter.cpp
+exMathInterpreter_CXXFLAGS        = $(COM_CXXFLAGS)
+exMathInterpreter_LDFLAGS         = -L../lib/.libs -lLatAnalyze
 
-exPlot_SOURCES                   = exPlot.cpp
-exPlot_CFLAGS                    = -g -O2
-exPlot_LDFLAGS                   = -L../lib/.libs -lLatAnalyze
+exPlot_SOURCES                    = exPlot.cpp
+exPlot_CXXFLAGS                   = $(COM_CXXFLAGS)
+exPlot_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
 
-exRand_SOURCES                   = exRand.cpp
-exRand_CFLAGS                    = -g -O2
-exRand_LDFLAGS                   = -L../lib/.libs -lLatAnalyze
+exRand_SOURCES                    = exRand.cpp
+exRand_CXXFLAGS                   = $(COM_CXXFLAGS)
+exRand_LDFLAGS                    = -L../lib/.libs -lLatAnalyze
 
-exRootFinder_SOURCES             = exRootFinder.cpp
-exRootFinder_CFLAGS              = -g -O2
-exRootFinder_LDFLAGS             = -L../lib/.libs -lLatAnalyze
+exRootFinder_SOURCES              = exRootFinder.cpp
+exRootFinder_CXXFLAGS             = $(COM_CXXFLAGS)
+exRootFinder_LDFLAGS              = -L../lib/.libs -lLatAnalyze
 
 ACLOCAL_AMFLAGS = -I .buildutils/m4
diff --git a/examples/exFit.cpp b/examples/exFit.cpp
index 1da5ed7..850a142 100644
--- a/examples/exFit.cpp
+++ b/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);
-        y_k          = f(&x_k, exactPar) + rg.gaussian(0.0, yErr);
-        printf("% 8e % 8e % 8e % 8e\n", x_k, xErr, y_k, yErr);
-        data.x(k) = x_k;
-        data.y(k) = y_k;
+        x1_k      = rg.gaussian(k*dx, xErr);
+        x2_k      = rg.gaussian(k*dx, xErr);
+        data.x(k) = x1_k;
+        data.y(k) = rg.gaussian(f(&x2_k, exactPar), yErr);
+        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 << " chi^2/ndof= " << p.getChi2PerDof();
-    cout << " p-value= " << p.getPValue() <<endl;
+    cout << "a= " << p(0) << " b= " << p(1) << endl;
+    cout << "chi^2/ndof= " << p.getChi2PerDof() << endl;
+    cout << "p-value= " << p.getPValue() <<endl;
     
     return EXIT_SUCCESS;
 }
diff --git a/examples/exFitSample.cpp b/examples/exFitSample.cpp
new file mode 100644
index 0000000..7c6ddd3
--- /dev/null
+++ b/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;
+}
diff --git a/lib/XYSampleData.cpp b/lib/XYSampleData.cpp
index 9a22d28..1f09bfe 100644
--- a/lib/XYSampleData.cpp
+++ b/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;
     }
@@ -232,7 +240,7 @@ SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init,
     FOR_STAT_ARRAY(result, s)
     {
         setDataToSample(s);
-        sampleResult = data_.fit(minimizer, init, v);
+        sampleResult    = data_.fit(minimizer, init, v);
         result[s]       = sampleResult;
         result.chi2_[s] = sampleResult.getChi2();
         result.nDof_    = sampleResult.getNDof();
@@ -265,7 +273,7 @@ void XYSampleData::computeVarMat(void)
     {
         // initialize data if necessary
         setDataToSample(central);
-        
+
         // compute relevant sizes
         Index size = 0, ySize = 0;
         
@@ -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>());
 }
diff --git a/utils/Makefile.am b/utils/Makefile.am
index 75b2a77..83648cc 100644
--- a/utils/Makefile.am
+++ b/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