[Speculative] Modifications used for reweighting

libnnpdf/src/NNPDF/chisquared.h

@@ -12,5 +12,10 @@
 #include "dataset.h"
 
 namespace NNPDF{
+  matrix<double> ComputeCovMat(CommonData const& cd, std::vector<double> const& t0);
+  matrix<double> ComputeSqrtMat(matrix<double> const& inmatrix);
+  void ComputeChi2_basic(int const nDat, int const nMem,
+                   const double* data, matrix<double> const& L,
+                   real *const& theory, real *chi2);
   template<class T> void ComputeChi2(const T*, int const&, real *const&, real *);
 }

libnnpdf/src/NNPDF/thpredictions.h

@@ -46,6 +46,13 @@ namespace NNPDF
 
         ThPredictions(const PDFSet*, const PDFSet*, const FKTable*);    //!< Different-beam constructor
 
+        // Empty constructor
+        ThPredictions(std::string pdfname,
+                      std::string setname,
+                      int nPDF,
+                      int nDat,
+                      PDFSet::erType);
+
         ThPredictions(const ThPredictions&);  //!< Copy-constructor
         friend void swap(ThPredictions&, ThPredictions&);
         ThPredictions& operator=(ThPredictions); //!< Copy-assignment

libnnpdf/src/NNPDF/utils.h

@@ -149,9 +149,10 @@ std::string joinpath(const std::initializer_list<std::string> &list);
     size_t const& size(size_t dim) const { return _size[dim]; } //!< Returns the (row,col) size pair.
     T&       operator()(size_t i, size_t j)       { return _data[i*_size[1]+j]; }
     T const& operator()(size_t i, size_t j) const { return _data[i*_size[1]+j]; }
-    //TODO: Does this have to be const? In any case there
-    //should be a const version.
-    T const * data () const {return _data.data();} //!< Return the underlying buffer.
+
+    // Data access
+    T *       data ()       {return _data.data();}  //!< Return the underlying buffer.
+    T const * data () const {return _data.data();}  //!< Return the underlying buffer (const version).
 
   private:
     std::array<size_t, 2> _size; //!< the dimension pair
@@ -181,9 +182,6 @@ std::string joinpath(const std::initializer_list<std::string> &list);
   real ComputeMom(int const& n, const real *x, int const& m);//!< Compute mth moment of distribution
   void Compute68cl(std::vector<real> const& x, real &up, real &dn);//!< Compute the 68% c.l.
   void Compute95cl(std::vector<real> const& x, real &up, real &dn);//!< Compute the 95% c.l.
-
-  void CholeskyDecomposition(matrix<double> const& inmatrix, matrix<double> & sqrtmat);
-
 }
 
 /*! @} */

libnnpdf/src/chisquared.cc

@@ -5,18 +5,91 @@
 // Authors: Nathan Hartland,  n.p.hartland@ed.ac.uk
 //          Stefano Carrazza, stefano.carrazza@mi.infn.it
 
+#include "NNPDF/exceptions.h"
 #include "NNPDF/chisquared.h"
+#include "gsl/gsl_matrix.h"
+#include "gsl/gsl_linalg.h"
 
 namespace NNPDF
 {
+  /**
+   * Generate covariance matrix from CommonData and a t0 vector
+   */
+  matrix<double> ComputeCovMat(CommonData const& cd, std::vector<double> const& t0)
+  {
+    const int ndat = cd.GetNData();
+    const int nsys = cd.GetNSys();
 
-  template<class T>
-  void ComputeChi2(const T* set, int const& nMem, real *const& theory, real *chi2)
+    if (t0.size() != ndat)
+      throw LengthError("ComputeCovMat","invalid number of points in t0 vector!");
+
+    auto CovMat = NNPDF::matrix<double>(ndat, ndat);
+    for (int i = 0; i < ndat; i++)
+    {
+      for (int j = 0; j < ndat; j++)
+      {
+        double sig    = 0.0;
+        double signor = 0.0;
+
+        if (i == j)
+          sig += pow(cd.GetStat(i),2); // stat error
+
+        for (int l = 0; l < nsys; l++)
+        {
+          sysError const& isys = cd.GetSys(i,l);
+          sysError const& jsys = cd.GetSys(j,l);
+          if (isys.name != jsys.name)
+              throw RuntimeException("ComputeCovMat", "Inconsistent naming of systematics");
+          if (isys.name == "SKIP")
+              continue;
+          const bool is_correlated = ( isys.name != "UNCORR" && isys.name !="THEORYUNCORR");
+          if (i == j || is_correlated)
+            switch (isys.type)
+            {
+              case ADD:   sig    += isys.add *jsys.add;  break;
+              case MULT:  signor += isys.mult*jsys.mult; break;
+              case UNSET: throw RuntimeException("ComputeCovMat", "UNSET systype encountered");
+            }
+        }
+
+        CovMat(i, j) = sig + signor*t0[i]*t0[j]*1e-4;
+      }
+    }
+    return CovMat;
+  }
+
+  matrix<double> ComputeSqrtMat(matrix<double> const& inmatrix)
   {
-    matrix<double> const& L = set->GetSqrtCov();
-    const double* data  = set->GetData();
-    const int nDat      = set->GetNData();
+    const size_t n = inmatrix.size(0);
+    if (n <= 0)
+      throw LengthError("CholeskyDecomposition","attempting a decomposition of an empty matrix!");
+
+    gsl_matrix_const_view inmatrix_view = gsl_matrix_const_view_array(inmatrix.data(), n, n);
+    const gsl_matrix *inmatrix_gsl = &(inmatrix_view.matrix);
 
+    matrix<double> sqrtmat(n,n);
+    gsl_matrix_view sqrtmat_view = gsl_matrix_view_array(sqrtmat.data(), n, n);
+    gsl_matrix *sqrtmat_gsl = &(sqrtmat_view.matrix);
+
+    // Copy and decompose inmatrix
+    const int copy = gsl_matrix_memcpy (sqrtmat_gsl, inmatrix_gsl);
+    if (copy != 0 ) throw RuntimeException("CholeskyDecomposition", "Error encountered in gsl matrix copy");
+    const int decomp = gsl_linalg_cholesky_decomp(sqrtmat_gsl);
+    if (decomp != 0 ) throw RuntimeException("CholeskyDecomposition", "Error encountered in gsl decomposition");
+
+    // Zero upper-diagonal part of matrix left by gsl (probably unneccesary)
+    for (int i = 0; i < n; i++)
+      for (int j = 0; j > i; j++)
+        sqrtmat(i, j) = 0;
+
+    return sqrtmat;
+  }
+
+  // TODO to sort this out, need to make data and theory vectors
+  void ComputeChi2_basic(int const nDat, int const nMem,
+                   const double* data, matrix<double> const& L,
+                   real *const& theory, real *chi2)
+  {
     // Forward solve Lx = diffs
     double x[nDat];
     for (int n = 0; n < nMem; n++)
@@ -31,6 +104,17 @@ namespace NNPDF
    return;
   }
 
+  template<class T>
+  void ComputeChi2(const T* set, int const& nMem, real *const& theory, real *chi2)
+  {
+    matrix<double> const& L = set->GetSqrtCov();
+    const double* data  = set->GetData();
+    const int nDat      = set->GetNData();
+
+    ComputeChi2_basic(nDat, nMem, data, L, theory, chi2);
+    return;
+  }
+
   template void ComputeChi2<Experiment>(const Experiment* set, int const& nMem, real *const& theory, real *chi2);
   template void ComputeChi2<DataSet>(const DataSet* set, int const& nMem, real *const& theory, real *chi2);
 }

libnnpdf/src/dataset.cc

@@ -20,6 +20,7 @@
 #include <memory>
 
 #include "NNPDF/dataset.h"
+#include "NNPDF/chisquared.h"
 #include "NNPDF/fastkernel.h"
 #include "NNPDF/thpredictions.h"
 #include "NNPDF/randomgenerator.h"
@@ -70,38 +71,8 @@ DataSet::~DataSet()
  */
 void DataSet::GenCovMat() const
 {
-  fCovMat.clear();
-  fSqrtCov.clear();
-  fCovMat.resize(fNData, fNData, 0);
-  fSqrtCov.resize(fNData, fNData, 0);
-
-  if (fNData <= 0)
-    throw LengthError("DataSet::GenCovMat","invalid number of datapoints!");
-
-  for (int i = 0; i < fNData; i++)
-    for (int j = 0; j < fNData; j++)
-    {
-      double sig    = 0.0;
-      double signor = 0.0;
-
-      if (i == j)
-        sig += fStat[i]*fStat[i]; // stat error
-
-      for (int l = 0; l < fNSys; l++)
-        if (fSys[i][l].name.compare("SKIP")!=0)
-          if (i == j || ( fSys[i][l].name.compare("UNCORR")!=0 && fSys[i][l].name.compare("THEORYUNCORR")!=0))
-            switch (fSys[i][l].type)
-            {
-              case ADD: sig += fSys[i][l].add*fSys[j][l].add; break; // additive systematics
-              case MULT: signor += fSys[i][l].mult*fSys[j][l].mult; break; // multiplicative systematics
-              case UNSET: throw RuntimeException("DataSet::GenCovMat", "UNSET systype encountered");
-            }
-
-      fCovMat(i, j) = sig + signor*fT0Pred[i]*fT0Pred[j]*1e-4;
-    }
-
-  // Compute sqrt of covmat
-  CholeskyDecomposition(fCovMat, fSqrtCov);
+  fCovMat = ComputeCovMat(*this, fT0Pred);
+  fSqrtCov = ComputeSqrtMat(fCovMat);
 }
 
 void DataSet::RescaleErrors(const double mult)

libnnpdf/src/experiments.cc

@@ -17,6 +17,7 @@
 #include <numeric>
 
 #include "NNPDF/experiments.h"
+#include "NNPDF/chisquared.h"
 #include "NNPDF/pdfset.h"
 #include "NNPDF/dataset.h"
 #include "NNPDF/thpredictions.h"
@@ -427,9 +428,7 @@ void Experiment::PullData()
 void Experiment::GenCovMat()
 {
   fCovMat.clear();
-  fSqrtCov.clear();
   fCovMat.resize(fNData, fNData, 0);
-  fSqrtCov.resize(fNData, fNData, 0);
 
   for (int i = 0; i < fNData; i++) {
     // Diagonal case
@@ -484,7 +483,7 @@ void Experiment::GenCovMat()
     }
   }
 
-  CholeskyDecomposition(fCovMat, fSqrtCov);
+  fSqrtCov = ComputeSqrtMat(fCovMat);
 }
 
 void Experiment::ExportCovMat(string filename)

libnnpdf/src/thpredictions.cc

@@ -208,6 +208,20 @@ fEtype(o.fEtype)
   fObs[i] = o.fObs[i];
 }
 
+/**
+ * Empty Constructor
+ */
+ThPredictions::ThPredictions(std::string pdfname, std::string setname, int nPDF, int nDat, PDFSet::erType erty):
+fObs(new real[nPDF*nDat]()),
+fTconv(0),
+fNpdf(nPDF),
+fNData(nDat),
+fPDFName(pdfname),
+fSetName(setname),
+fEtype(erty)
+{
+}
+
 
 void NNPDF::swap(ThPredictions& lhs, ThPredictions& rhs)
 {

libnnpdf/src/utils.cc

@@ -433,23 +433,4 @@ std::string joinpath(const std::initializer_list<std::string> &list)
     return sum;
   }
 
-  // *********** Cholesky decomposition of a matrix ***************
-
-  void CholeskyDecomposition(const matrix<double> &inmatrix, matrix<double> &sqrtmat)
-  {
-    const size_t n = inmatrix.size(0);
-    if (n <= 0)
-      throw LengthError("CholeskyDecomposition","attempting a decomposition of an empty matrix!");
-    gsl_matrix* mat = gsl_matrix_calloc(n, n);
-    for (int i = 0; i < n; i++)
-      for (int j = 0; j < n; j++)
-        gsl_matrix_set(mat, i, j, inmatrix(i, j));
-
-    const int decomp = gsl_linalg_cholesky_decomp(mat);
-    for (int i = 0; i < n; i++)
-      for (int j = 0; j <= i; j++)
-        sqrtmat(i, j) =  gsl_matrix_get(mat, i, j);
-    gsl_matrix_free (mat);
-  }
-
 }

validphys2/src/validphys/reweighting.py

@@ -73,7 +73,7 @@ def nnpdf_weights(chi2_data_for_reweighting_experiments):
 def effective_number_of_replicas(w):
     N = len(w)
     w = w*N/np.sum(w)
-    return np.exp(np.nansum(w*np.log(N/w))/N)
+    return np.exp(np.nansum(w*np.log(N)-w*np.log(w))/N)
 
 @table
 def reweighting_stats(pdf, nnpdf_weights, p_alpha_study):