diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index d8c2982685..8c3d28b814 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -12,5 +12,10 @@ #include "dataset.h" namespace NNPDF{ + matrix ComputeCovMat(CommonData const& cd, std::vector const& t0); + matrix ComputeSqrtMat(matrix const& inmatrix); + void ComputeChi2_basic(int const nDat, int const nMem, + const double* data, matrix const& L, + real *const& theory, real *chi2); template void ComputeChi2(const T*, int const&, real *const&, real *); } diff --git a/libnnpdf/src/NNPDF/thpredictions.h b/libnnpdf/src/NNPDF/thpredictions.h index ee7e322ede..2aeadf0f97 100644 --- a/libnnpdf/src/NNPDF/thpredictions.h +++ b/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 diff --git a/libnnpdf/src/NNPDF/utils.h b/libnnpdf/src/NNPDF/utils.h index 77888a5a75..cd27f287e0 100644 --- a/libnnpdf/src/NNPDF/utils.h +++ b/libnnpdf/src/NNPDF/utils.h @@ -149,9 +149,10 @@ std::string joinpath(const std::initializer_list &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; //!< the dimension pair @@ -181,9 +182,6 @@ std::string joinpath(const std::initializer_list &list); real ComputeMom(int const& n, const real *x, int const& m);//!< Compute mth moment of distribution void Compute68cl(std::vector const& x, real &up, real &dn);//!< Compute the 68% c.l. void Compute95cl(std::vector const& x, real &up, real &dn);//!< Compute the 95% c.l. - - void CholeskyDecomposition(matrix const& inmatrix, matrix & sqrtmat); - } /*! @} */ diff --git a/libnnpdf/src/chisquared.cc b/libnnpdf/src/chisquared.cc index 3fc9c6062b..021fcaa917 100644 --- a/libnnpdf/src/chisquared.cc +++ b/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 ComputeCovMat(CommonData const& cd, std::vector const& t0) + { + const int ndat = cd.GetNData(); + const int nsys = cd.GetNSys(); - template - 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(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 ComputeSqrtMat(matrix const& inmatrix) { - matrix 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 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 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 + void ComputeChi2(const T* set, int const& nMem, real *const& theory, real *chi2) + { + matrix 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(const Experiment* set, int const& nMem, real *const& theory, real *chi2); template void ComputeChi2(const DataSet* set, int const& nMem, real *const& theory, real *chi2); } diff --git a/libnnpdf/src/dataset.cc b/libnnpdf/src/dataset.cc index 0f8b99e511..c53f9dfbbb 100644 --- a/libnnpdf/src/dataset.cc +++ b/libnnpdf/src/dataset.cc @@ -20,6 +20,7 @@ #include #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) diff --git a/libnnpdf/src/experiments.cc b/libnnpdf/src/experiments.cc index a685b66f28..165f5177ba 100644 --- a/libnnpdf/src/experiments.cc +++ b/libnnpdf/src/experiments.cc @@ -17,6 +17,7 @@ #include #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) diff --git a/libnnpdf/src/thpredictions.cc b/libnnpdf/src/thpredictions.cc index 0b5c1a5f08..f78bce33dd 100644 --- a/libnnpdf/src/thpredictions.cc +++ b/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) { diff --git a/libnnpdf/src/utils.cc b/libnnpdf/src/utils.cc index bc5764a028..64e624ca6b 100644 --- a/libnnpdf/src/utils.cc +++ b/libnnpdf/src/utils.cc @@ -433,23 +433,4 @@ std::string joinpath(const std::initializer_list &list) return sum; } - // *********** Cholesky decomposition of a matrix *************** - - void CholeskyDecomposition(const matrix &inmatrix, matrix &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); - } - } diff --git a/validphys2/src/validphys/reweighting.py b/validphys2/src/validphys/reweighting.py index 46fd51a814..bcb98785c7 100644 --- a/validphys2/src/validphys/reweighting.py +++ b/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):