From 36cf2b8e12844164f485832d40fa58b4fe8129bc Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Wed, 14 Feb 2018 17:52:46 +0100 Subject: [PATCH 01/16] Separated dataset covariance matrix generation, new function for generation of sqrt matrix --- libnnpdf/src/NNPDF/chisquared.h | 2 + libnnpdf/src/NNPDF/utils.h | 3 -- libnnpdf/src/chisquared.cc | 70 +++++++++++++++++++++++++++++++ libnnpdf/src/dataset.cc | 73 ++++++++++----------------------- libnnpdf/src/experiments.cc | 67 +++++++++++++++--------------- libnnpdf/src/utils.cc | 19 --------- 6 files changed, 127 insertions(+), 107 deletions(-) diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index d8c2982685..1d29472f42 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -12,5 +12,7 @@ #include "dataset.h" namespace NNPDF{ + matrix ComputeCovMat(CommonData const& cd, std::vector const& t0); + matrix ComputeSqrtMat(matrix const& inmatrix); template void ComputeChi2(const T*, int const&, real *const&, real *); } diff --git a/libnnpdf/src/NNPDF/utils.h b/libnnpdf/src/NNPDF/utils.h index fac4618270..5b8d7e3431 100644 --- a/libnnpdf/src/NNPDF/utils.h +++ b/libnnpdf/src/NNPDF/utils.h @@ -134,9 +134,6 @@ namespace NNPDF 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..5555d76f60 100644 --- a/libnnpdf/src/chisquared.cc +++ b/libnnpdf/src/chisquared.cc @@ -5,10 +5,80 @@ // 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(); + if (ndat <= 0) + throw LengthError("ComputeCovMat","invalid number of datapoints!"); + 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) + { + const size_t n = inmatrix.size(0); + if (n <= 0) + throw LengthError("CholeskyDecomposition","attempting a decomposition of an empty matrix!"); + matrix sqrtmat(n,n); + + 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); + if (decomp != 0 ) throw RuntimeException("CholeskyDecomposition", "Error encountered in gsl"); + + 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); + return sqrtmat; + } template void ComputeChi2(const T* set, int const& nMem, real *const& theory, real *chi2) diff --git a/libnnpdf/src/dataset.cc b/libnnpdf/src/dataset.cc index ec7c10c0fe..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" @@ -62,7 +63,7 @@ DataSet::DataSet(const DataSet& set, std::vector const& mask): * Cleanup memory */ DataSet::~DataSet() -{ +{ } /** @@ -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) @@ -109,9 +80,9 @@ void DataSet::RescaleErrors(const double mult) for (int i = 0; i < fNData; i++) { fStat[i] *= mult; // Rescale stat and sys to new data value - + for (int l = 0; l < fNSys; l++) - { + { fSys[i][l].add *= mult; fSys[i][l].mult *= mult; } @@ -178,11 +149,11 @@ void DataSet::MakeArtificial() { std::cout << "-- Generating replica data for " << fSetName << std::endl; - double *rand = new double[fNSys]; + double *rand = new double[fNSys]; double *xnor = new double[fNData]; double *artdata = new double[fNData]; sysType rST[2] = {ADD,MULT}; - + NNPDF::RandomGenerator* rng = NNPDF::RandomGenerator::GetRNG(); // Generate random systematics @@ -193,26 +164,26 @@ void DataSet::MakeArtificial() fSys[0][l].type = rST[rng->GetRandomUniform(2)]; for (int i = 1; i < fNData; i++) fSys[i][l].type = fSys[0][l].type; - } + } int genTries = 0; // Number of attempts at generation bool genRep = false; - while ( genTries < 1E6 ) + while ( genTries < 1E6 ) { - + // Update the data for (int i = 0; i < fNData; i++) { const double xstat = rng->GetRandomGausDev(1.0)*fStat[i]; //Noise from statistical uncertainty - + double xadd = 0; xnor[i] = 1.0; - + for (int l = 0; l < fNSys; l++) { if (fSys[i][l].name.compare("THEORYCORR")==0) continue; // Skip theoretical uncertainties if (fSys[i][l].name.compare("THEORYUNCORR")==0) continue; // Skip theoretical uncertainties - if (fSys[i][l].name.compare("SKIP")==0) continue; // Skip uncertainties + if (fSys[i][l].name.compare("SKIP")==0) continue; // Skip uncertainties if (fSys[i][l].name.compare("UNCORR")==0) // Noise from uncorrelated systematics { switch (fSys[i][l].type) @@ -235,13 +206,13 @@ void DataSet::MakeArtificial() // Generate the artificial data artdata[i] = xnor[i] * ( fData[i] + xadd + xstat ); - + // Only generates positive artificial data (except for closure tests) if ( artdata[i] < 0 and fProc[i].find("ASY") != std::string::npos ) break; } - + // Stop when all the datapoints are positive genRep = true; break; @@ -253,12 +224,12 @@ void DataSet::MakeArtificial() // Update data in set UpdateData(artdata); fIsArtificial = true; - + delete[] rand; delete[] xnor; delete[] artdata; - // Note DO NOT Regenerate covariance matrices + // Note DO NOT Regenerate covariance matrices } /** @@ -273,11 +244,11 @@ void DataSet::UpdateData(double *newdata) // Update data only void DataSet::UpdateData(double *newdata, double* uncnorm) // Update data only and rescale systematics by normalisation { UpdateData(newdata); - + for (int i = 0; i < fNData; i++) { fStat[i] *= uncnorm[i]; // Rescale stat and sys to new data value - + for (int l = 0; l < fNSys; l++) fSys[i][l].add *= uncnorm[i]; } @@ -286,8 +257,8 @@ void DataSet::UpdateData(double *newdata, double* uncnorm) // Update data o void DataSet::UpdateData(double *newdata, sysType *type) // Update data and change systypes { UpdateData(newdata); - + for (int l = 0; l < fNSys; l++) for (int i = 0; i < fNData; i++) - fSys[i][l].type = type[l]; + fSys[i][l].type = type[l]; } diff --git a/libnnpdf/src/experiments.cc b/libnnpdf/src/experiments.cc index fb027619aa..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" @@ -56,7 +57,7 @@ fIsT0(false) if (sets[i].IsT0() != fIsT0) throw UserError("Experiment::Experiment", "Supplied datasets must all be either T0 or EXP"); } - + // Pull data from datasets PullData(); GenCovMat(); @@ -83,7 +84,7 @@ fIsT0(exp.fIsT0) // Copy datasets for (size_t i=0; i< exp.fSets.size(); i++) fSets.push_back(exp.fSets[i]); - + // Pull data from datasets PullData(); } @@ -107,7 +108,7 @@ fIsT0(exp.fIsT0) // Copy datasets for (size_t i=0; i< sets.size(); i++) fSets.push_back(sets[i]); - + // Pull data from datasets PullData(); GenCovMat(); @@ -118,14 +119,14 @@ fIsT0(exp.fIsT0) */ Experiment::~Experiment() { - ClearLocalData(); + ClearLocalData(); } /* * Clears data pulled from DataSets */ void Experiment::ClearLocalData() -{ +{ // Already clear if (!fData) return; @@ -133,17 +134,17 @@ void Experiment::ClearLocalData() delete[] fData; delete[] fStat; delete[] fT0Pred; - + for (int s = 0; s < GetNSet(); s++) delete[] fSetSysMap[s]; - + delete[] fSetSysMap; - + for (int i = 0; i < fNData; i++) delete[] fSys[i]; delete[] fSys; - - fNSys = 0; + + fNSys = 0; } /** @@ -160,7 +161,7 @@ void Experiment::MakeReplica() RandomGenerator* rng = RandomGenerator::GetRNG(); - double *rand = new double[fNSys]; + double *rand = new double[fNSys]; double *xnor = new double[fNData]; double *artdata = new double[fNData]; sysType rST[2] = {ADD,MULT}; @@ -185,20 +186,20 @@ void Experiment::MakeReplica() fSys[0][l].type = rST[rng->GetRandomUniform(2)]; for (int i = 1; i < fNData; i++) fSys[i][l].type = fSys[0][l].type; - } - + } + for (int i = 0; i < fNData; i++) // should rearrange to update set-by-set -- nh { double xstat = rng->GetRandomGausDev(1.0)*fStat[i]; //Noise from statistical uncertainty - + double xadd = 0; xnor[i] = 1.0; - + for (int l = 0; l < fNSys; l++) { if (fSys[i][l].name.compare("THEORYCORR")==0) continue; // Skip theoretical uncertainties if (fSys[i][l].name.compare("THEORYUNCORR")==0) continue; // Skip theoretical uncertainties - if (fSys[i][l].name.compare("SKIP")==0) continue; // Skip uncertainties + if (fSys[i][l].name.compare("SKIP")==0) continue; // Skip uncertainties if (fSys[i][l].name.compare("UNCORR")==0) // Noise from uncorrelated systematics { switch (fSys[i][l].type) @@ -219,7 +220,7 @@ void Experiment::MakeReplica() } } artdata[i] = xnor[i] * ( fData[i] + xadd + xstat ); - + // Only generates positive artificial data (except for closure tests and asymmetry data) if (artdata[i] < 0 && !fIsClosure && proctype[i].find("ASY") == std::string::npos ) { @@ -235,22 +236,22 @@ void Experiment::MakeReplica() { sysType *type = new sysType[fSets[s].GetNSys()]; for (int l = 0; l < fSets[s].GetNSys(); l++) - type[l] = fSys[0][fSetSysMap[s][l]].type; + type[l] = fSys[0][fSetSysMap[s][l]].type; fSets[s].UpdateData(artdata+index, type); fSets[s].SetArtificial(true); index+=fSets[s].GetNData(); delete[] type; } - + // Update local data and recompute covariance matrix PullData(); //GenCovMat(); //Use standard t0 covariance matrix for all replicas - + delete[] rand; delete[] xnor; delete[] artdata; delete[] proctype; - + // Now the fData is artificial fIsArtificial = true; } @@ -343,21 +344,21 @@ void Experiment::PullData() if ( find(special_types.begin(), special_types.end(), testsys.name) != special_types.end() ) { // This systematic is an unnamed/special type, add it to the map and the total systematics vector - fSetSysMap[i][l] = total_systematics.size(); + fSetSysMap[i][l] = total_systematics.size(); total_systematics.emplace_back(testsys); } - else + else { // This is a named systematic, we need to cross-check against existing named systematics bool found_systematic = false; for ( size_t k=0; k < total_systematics.size(); k++ ) { sysError& existing_systematic = total_systematics[k]; - if (testsys.name == existing_systematic.name ) + if (testsys.name == existing_systematic.name ) { // This already exists in total_systematics, it's not a new systematic found_systematic = true; - fSetSysMap[i][l] = k; + fSetSysMap[i][l] = k; // Perform some consistency checks if (testsys.type != existing_systematic.type || testsys.isRAND != existing_systematic.isRAND) throw RangeError("Experiment::PullData","Systematic " + testsys.name + " definition not consistant between datasets"); @@ -365,22 +366,22 @@ void Experiment::PullData() } } // If the systematic doesn't already exist in the list, add it - if ( found_systematic == false ) + if ( found_systematic == false ) { - fSetSysMap[i][l] = total_systematics.size(); + fSetSysMap[i][l] = total_systematics.size(); total_systematics.emplace_back(testsys); } } } } - + // Initialise data arrays fData = new double[fNData]; fStat = new double[fNData]; fT0Pred = new double[fNData]; fNSys = total_systematics.size(); fSys = new sysError*[fNData]; - + // Fill the experiment arrays with the values from its subsets int idat = 0; for (int s = 0; s < GetNSet(); s++) @@ -393,7 +394,7 @@ void Experiment::PullData() // Loop over experimental systematics, find if there is a corresponding dataset systematic // and set it accordingly. - fSys[idat] = new sysError[fNSys]; + fSys[idat] = new sysError[fNSys]; for (int l = 0; l < fNSys; l++) { fSys[idat][l].name = total_systematics[l].name; @@ -417,7 +418,7 @@ void Experiment::PullData() } idat++; } - + return; } @@ -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/utils.cc b/libnnpdf/src/utils.cc index 20beb71f2e..1da34fb1a1 100644 --- a/libnnpdf/src/utils.cc +++ b/libnnpdf/src/utils.cc @@ -429,23 +429,4 @@ namespace NNPDF 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); - } - } From 0e891d0a4ec04a1675971e96b54902d38a70fe10 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Wed, 14 Feb 2018 18:43:14 +0100 Subject: [PATCH 02/16] Add a basic form of the chi2 calculation without using any special types --- libnnpdf/src/NNPDF/chisquared.h | 3 +++ libnnpdf/src/chisquared.cc | 21 +++++++++++++++------ 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index 1d29472f42..addaba1a83 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -14,5 +14,8 @@ 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/chisquared.cc b/libnnpdf/src/chisquared.cc index 5555d76f60..0a672bf740 100644 --- a/libnnpdf/src/chisquared.cc +++ b/libnnpdf/src/chisquared.cc @@ -80,13 +80,11 @@ namespace NNPDF return sqrtmat; } - template - void ComputeChi2(const T* set, int const& nMem, real *const& theory, real *chi2) + // 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) { - matrix const& L = set->GetSqrtCov(); - const double* data = set->GetData(); - const int nDat = set->GetNData(); - // Forward solve Lx = diffs double x[nDat]; for (int n = 0; n < nMem; n++) @@ -101,6 +99,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); } From bf087a08d6beed206eb301aa33311abbea859120 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Thu, 15 Feb 2018 13:02:06 +0100 Subject: [PATCH 03/16] Better handling of very small weights in N_eff:= --- validphys2/src/validphys/reweighting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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): From 4d68741eeccc535d73cb41240ff0687c0a079b4c Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Thu, 15 Feb 2018 13:02:33 +0100 Subject: [PATCH 04/16] Empty ThPredictions constructor --- libnnpdf/src/NNPDF/thpredictions.h | 9 ++- libnnpdf/src/thpredictions.cc | 96 +++++++++++++++++------------- 2 files changed, 63 insertions(+), 42 deletions(-) diff --git a/libnnpdf/src/NNPDF/thpredictions.h b/libnnpdf/src/NNPDF/thpredictions.h index 89050754d3..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 @@ -67,7 +74,7 @@ namespace NNPDF static void Convolute(const PDFSet*,const FKTable*,real*); static void Convolute(const PDFSet*,const FKSet*,real*); static void Convolute(const PDFSet*,const PDFSet*,const FKTable*,real*); - + // Get Methods real* GetObs() const {return fObs;}; //!< Return Obs array real GetObs(int const& idat, int const& imem) const {return fObs[idat*fNpdf + imem];}; diff --git a/libnnpdf/src/thpredictions.cc b/libnnpdf/src/thpredictions.cc index 8be47348b5..f78bce33dd 100644 --- a/libnnpdf/src/thpredictions.cc +++ b/libnnpdf/src/thpredictions.cc @@ -68,7 +68,7 @@ bool ThPredictions::Verbose = true; __m256 t2 = _mm256_hadd_ps(t1,t1); __m128 t3 = _mm256_extractf128_ps(t2,1); __m128 t4 = _mm_add_ss(_mm256_castps256_ps128(t2),t3); - retval = _mm_cvtss_f32(t4); + retval = _mm_cvtss_f32(t4); _mm_empty(); */ @@ -100,7 +100,7 @@ fEtype(pdfset->GetEtype()) { // New theory array fObs = new real[fNData*fNpdf]; - + Timer timer=Timer(); timer.start(); @@ -127,7 +127,7 @@ fEtype(pdfset->GetEtype()) { // New theory array fObs = new real[fNData*fNpdf]; - + Timer timer=Timer(); timer.start(); @@ -150,7 +150,7 @@ fEtype(pdfset->GetEtype()) { // New theory array fObs = new real[fNData*fNpdf]; - + Timer timer=Timer(); timer.start(); @@ -180,7 +180,7 @@ fEtype(pdf1->GetEtype()) // New theory array fObs = new real[fNData*fNpdf]; - + Timer timer=Timer(); timer.start(); @@ -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) { @@ -257,9 +271,9 @@ ThPredictions& ThPredictions::operator+=(const ThPredictions& o) throw EvaluationError("ThPredictions::operator+=","Cannot sum predictions with different numbers of datapoints/PDF members"); if (fEtype != PDFSet::erType::ER_MC && fEtype != PDFSet::erType::ER_MC68 && fEtype != PDFSet::erType::ER_MCT0 && ThPredictions::Verbose == true) - std::cerr << "ThPredictions::operator+= Warning: Observable summation undefined for ErrorTypes other than Monte-Carlo" <(pdf)); return; @@ -423,7 +437,7 @@ void ThPredictions::Convolute(const PDFSet *pdf1, const PDFSet *pdf2, const FKTa if (err != 0) throw RangeError("ThPredictions::Convolute","ThPredictions posix_memalign " + std::to_string(err)); - + // Fetch PDF array GetNZPDF(pdf1, pdf2, fk, pdf); @@ -451,7 +465,7 @@ void ThPredictions::Convolute(const PDFSet *pdf1, const PDFSet *pdf2, const FKTa } MPI::RecvJobs(NData, Npdf, theory); #endif - + // Delete pdfs free(reinterpret_cast(pdf)); return; @@ -463,7 +477,7 @@ void ThPredictions::Convolute(const PDFSet *pdf1, const PDFSet *pdf2, const FKTa void ThPredictions::Convolute(const PDFSet *pdfset, const FKSet *fkset, real* theory) { const int NSigma = fkset->GetNSigma(); - + if (NSigma == 1) ThPredictions::Convolute(pdfset,fkset->GetFK(0), theory ); else @@ -471,22 +485,22 @@ void ThPredictions::Convolute(const PDFSet *pdfset, const FKSet *fkset, real* th SigmaOp op = fkset->GetOperator(); const int NData = fkset->GetNDataFK(); const int Npdf = pdfset->GetMembers(); - + std::vector results; for (int i=0; iGetFK(i), results[i] ); - + op( NData*Npdf, results, theory ); - + for (int i=0; iGetPDF(xgrid[i],Q20,n,&EVLN[i*Nfl]); - + for (int fl=0; flGetPDF(xgrid[i],Q20,n,&EVLN1[i*Nfl]); pdf2->GetPDF(xgrid[i],Q20,n,&EVLN2[i*Nfl]); } - + for (int fl=0; fl Date: Wed, 14 Feb 2018 17:52:46 +0100 Subject: [PATCH 05/16] Separated dataset covariance matrix generation, new function for generation of sqrt matrix --- libnnpdf/src/NNPDF/chisquared.h | 2 + libnnpdf/src/NNPDF/utils.h | 3 -- libnnpdf/src/chisquared.cc | 70 +++++++++++++++++++++++++++++++ libnnpdf/src/dataset.cc | 73 ++++++++++----------------------- libnnpdf/src/experiments.cc | 67 +++++++++++++++--------------- libnnpdf/src/utils.cc | 19 --------- 6 files changed, 127 insertions(+), 107 deletions(-) diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index d8c2982685..1d29472f42 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -12,5 +12,7 @@ #include "dataset.h" namespace NNPDF{ + matrix ComputeCovMat(CommonData const& cd, std::vector const& t0); + matrix ComputeSqrtMat(matrix const& inmatrix); template void ComputeChi2(const T*, int const&, real *const&, real *); } diff --git a/libnnpdf/src/NNPDF/utils.h b/libnnpdf/src/NNPDF/utils.h index 77888a5a75..35c3af7837 100644 --- a/libnnpdf/src/NNPDF/utils.h +++ b/libnnpdf/src/NNPDF/utils.h @@ -181,9 +181,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..5555d76f60 100644 --- a/libnnpdf/src/chisquared.cc +++ b/libnnpdf/src/chisquared.cc @@ -5,10 +5,80 @@ // 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(); + if (ndat <= 0) + throw LengthError("ComputeCovMat","invalid number of datapoints!"); + 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) + { + const size_t n = inmatrix.size(0); + if (n <= 0) + throw LengthError("CholeskyDecomposition","attempting a decomposition of an empty matrix!"); + matrix sqrtmat(n,n); + + 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); + if (decomp != 0 ) throw RuntimeException("CholeskyDecomposition", "Error encountered in gsl"); + + 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); + return sqrtmat; + } template void ComputeChi2(const T* set, int const& nMem, real *const& theory, real *chi2) diff --git a/libnnpdf/src/dataset.cc b/libnnpdf/src/dataset.cc index ec7c10c0fe..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" @@ -62,7 +63,7 @@ DataSet::DataSet(const DataSet& set, std::vector const& mask): * Cleanup memory */ DataSet::~DataSet() -{ +{ } /** @@ -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) @@ -109,9 +80,9 @@ void DataSet::RescaleErrors(const double mult) for (int i = 0; i < fNData; i++) { fStat[i] *= mult; // Rescale stat and sys to new data value - + for (int l = 0; l < fNSys; l++) - { + { fSys[i][l].add *= mult; fSys[i][l].mult *= mult; } @@ -178,11 +149,11 @@ void DataSet::MakeArtificial() { std::cout << "-- Generating replica data for " << fSetName << std::endl; - double *rand = new double[fNSys]; + double *rand = new double[fNSys]; double *xnor = new double[fNData]; double *artdata = new double[fNData]; sysType rST[2] = {ADD,MULT}; - + NNPDF::RandomGenerator* rng = NNPDF::RandomGenerator::GetRNG(); // Generate random systematics @@ -193,26 +164,26 @@ void DataSet::MakeArtificial() fSys[0][l].type = rST[rng->GetRandomUniform(2)]; for (int i = 1; i < fNData; i++) fSys[i][l].type = fSys[0][l].type; - } + } int genTries = 0; // Number of attempts at generation bool genRep = false; - while ( genTries < 1E6 ) + while ( genTries < 1E6 ) { - + // Update the data for (int i = 0; i < fNData; i++) { const double xstat = rng->GetRandomGausDev(1.0)*fStat[i]; //Noise from statistical uncertainty - + double xadd = 0; xnor[i] = 1.0; - + for (int l = 0; l < fNSys; l++) { if (fSys[i][l].name.compare("THEORYCORR")==0) continue; // Skip theoretical uncertainties if (fSys[i][l].name.compare("THEORYUNCORR")==0) continue; // Skip theoretical uncertainties - if (fSys[i][l].name.compare("SKIP")==0) continue; // Skip uncertainties + if (fSys[i][l].name.compare("SKIP")==0) continue; // Skip uncertainties if (fSys[i][l].name.compare("UNCORR")==0) // Noise from uncorrelated systematics { switch (fSys[i][l].type) @@ -235,13 +206,13 @@ void DataSet::MakeArtificial() // Generate the artificial data artdata[i] = xnor[i] * ( fData[i] + xadd + xstat ); - + // Only generates positive artificial data (except for closure tests) if ( artdata[i] < 0 and fProc[i].find("ASY") != std::string::npos ) break; } - + // Stop when all the datapoints are positive genRep = true; break; @@ -253,12 +224,12 @@ void DataSet::MakeArtificial() // Update data in set UpdateData(artdata); fIsArtificial = true; - + delete[] rand; delete[] xnor; delete[] artdata; - // Note DO NOT Regenerate covariance matrices + // Note DO NOT Regenerate covariance matrices } /** @@ -273,11 +244,11 @@ void DataSet::UpdateData(double *newdata) // Update data only void DataSet::UpdateData(double *newdata, double* uncnorm) // Update data only and rescale systematics by normalisation { UpdateData(newdata); - + for (int i = 0; i < fNData; i++) { fStat[i] *= uncnorm[i]; // Rescale stat and sys to new data value - + for (int l = 0; l < fNSys; l++) fSys[i][l].add *= uncnorm[i]; } @@ -286,8 +257,8 @@ void DataSet::UpdateData(double *newdata, double* uncnorm) // Update data o void DataSet::UpdateData(double *newdata, sysType *type) // Update data and change systypes { UpdateData(newdata); - + for (int l = 0; l < fNSys; l++) for (int i = 0; i < fNData; i++) - fSys[i][l].type = type[l]; + fSys[i][l].type = type[l]; } diff --git a/libnnpdf/src/experiments.cc b/libnnpdf/src/experiments.cc index fb027619aa..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" @@ -56,7 +57,7 @@ fIsT0(false) if (sets[i].IsT0() != fIsT0) throw UserError("Experiment::Experiment", "Supplied datasets must all be either T0 or EXP"); } - + // Pull data from datasets PullData(); GenCovMat(); @@ -83,7 +84,7 @@ fIsT0(exp.fIsT0) // Copy datasets for (size_t i=0; i< exp.fSets.size(); i++) fSets.push_back(exp.fSets[i]); - + // Pull data from datasets PullData(); } @@ -107,7 +108,7 @@ fIsT0(exp.fIsT0) // Copy datasets for (size_t i=0; i< sets.size(); i++) fSets.push_back(sets[i]); - + // Pull data from datasets PullData(); GenCovMat(); @@ -118,14 +119,14 @@ fIsT0(exp.fIsT0) */ Experiment::~Experiment() { - ClearLocalData(); + ClearLocalData(); } /* * Clears data pulled from DataSets */ void Experiment::ClearLocalData() -{ +{ // Already clear if (!fData) return; @@ -133,17 +134,17 @@ void Experiment::ClearLocalData() delete[] fData; delete[] fStat; delete[] fT0Pred; - + for (int s = 0; s < GetNSet(); s++) delete[] fSetSysMap[s]; - + delete[] fSetSysMap; - + for (int i = 0; i < fNData; i++) delete[] fSys[i]; delete[] fSys; - - fNSys = 0; + + fNSys = 0; } /** @@ -160,7 +161,7 @@ void Experiment::MakeReplica() RandomGenerator* rng = RandomGenerator::GetRNG(); - double *rand = new double[fNSys]; + double *rand = new double[fNSys]; double *xnor = new double[fNData]; double *artdata = new double[fNData]; sysType rST[2] = {ADD,MULT}; @@ -185,20 +186,20 @@ void Experiment::MakeReplica() fSys[0][l].type = rST[rng->GetRandomUniform(2)]; for (int i = 1; i < fNData; i++) fSys[i][l].type = fSys[0][l].type; - } - + } + for (int i = 0; i < fNData; i++) // should rearrange to update set-by-set -- nh { double xstat = rng->GetRandomGausDev(1.0)*fStat[i]; //Noise from statistical uncertainty - + double xadd = 0; xnor[i] = 1.0; - + for (int l = 0; l < fNSys; l++) { if (fSys[i][l].name.compare("THEORYCORR")==0) continue; // Skip theoretical uncertainties if (fSys[i][l].name.compare("THEORYUNCORR")==0) continue; // Skip theoretical uncertainties - if (fSys[i][l].name.compare("SKIP")==0) continue; // Skip uncertainties + if (fSys[i][l].name.compare("SKIP")==0) continue; // Skip uncertainties if (fSys[i][l].name.compare("UNCORR")==0) // Noise from uncorrelated systematics { switch (fSys[i][l].type) @@ -219,7 +220,7 @@ void Experiment::MakeReplica() } } artdata[i] = xnor[i] * ( fData[i] + xadd + xstat ); - + // Only generates positive artificial data (except for closure tests and asymmetry data) if (artdata[i] < 0 && !fIsClosure && proctype[i].find("ASY") == std::string::npos ) { @@ -235,22 +236,22 @@ void Experiment::MakeReplica() { sysType *type = new sysType[fSets[s].GetNSys()]; for (int l = 0; l < fSets[s].GetNSys(); l++) - type[l] = fSys[0][fSetSysMap[s][l]].type; + type[l] = fSys[0][fSetSysMap[s][l]].type; fSets[s].UpdateData(artdata+index, type); fSets[s].SetArtificial(true); index+=fSets[s].GetNData(); delete[] type; } - + // Update local data and recompute covariance matrix PullData(); //GenCovMat(); //Use standard t0 covariance matrix for all replicas - + delete[] rand; delete[] xnor; delete[] artdata; delete[] proctype; - + // Now the fData is artificial fIsArtificial = true; } @@ -343,21 +344,21 @@ void Experiment::PullData() if ( find(special_types.begin(), special_types.end(), testsys.name) != special_types.end() ) { // This systematic is an unnamed/special type, add it to the map and the total systematics vector - fSetSysMap[i][l] = total_systematics.size(); + fSetSysMap[i][l] = total_systematics.size(); total_systematics.emplace_back(testsys); } - else + else { // This is a named systematic, we need to cross-check against existing named systematics bool found_systematic = false; for ( size_t k=0; k < total_systematics.size(); k++ ) { sysError& existing_systematic = total_systematics[k]; - if (testsys.name == existing_systematic.name ) + if (testsys.name == existing_systematic.name ) { // This already exists in total_systematics, it's not a new systematic found_systematic = true; - fSetSysMap[i][l] = k; + fSetSysMap[i][l] = k; // Perform some consistency checks if (testsys.type != existing_systematic.type || testsys.isRAND != existing_systematic.isRAND) throw RangeError("Experiment::PullData","Systematic " + testsys.name + " definition not consistant between datasets"); @@ -365,22 +366,22 @@ void Experiment::PullData() } } // If the systematic doesn't already exist in the list, add it - if ( found_systematic == false ) + if ( found_systematic == false ) { - fSetSysMap[i][l] = total_systematics.size(); + fSetSysMap[i][l] = total_systematics.size(); total_systematics.emplace_back(testsys); } } } } - + // Initialise data arrays fData = new double[fNData]; fStat = new double[fNData]; fT0Pred = new double[fNData]; fNSys = total_systematics.size(); fSys = new sysError*[fNData]; - + // Fill the experiment arrays with the values from its subsets int idat = 0; for (int s = 0; s < GetNSet(); s++) @@ -393,7 +394,7 @@ void Experiment::PullData() // Loop over experimental systematics, find if there is a corresponding dataset systematic // and set it accordingly. - fSys[idat] = new sysError[fNSys]; + fSys[idat] = new sysError[fNSys]; for (int l = 0; l < fNSys; l++) { fSys[idat][l].name = total_systematics[l].name; @@ -417,7 +418,7 @@ void Experiment::PullData() } idat++; } - + return; } @@ -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/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); - } - } From 0ad5fd9b985b845ef6fdf8b946fdd6ece54297f6 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Wed, 14 Feb 2018 18:43:14 +0100 Subject: [PATCH 06/16] Add a basic form of the chi2 calculation without using any special types --- libnnpdf/src/NNPDF/chisquared.h | 3 +++ libnnpdf/src/chisquared.cc | 21 +++++++++++++++------ 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index 1d29472f42..addaba1a83 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -14,5 +14,8 @@ 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/chisquared.cc b/libnnpdf/src/chisquared.cc index 5555d76f60..0a672bf740 100644 --- a/libnnpdf/src/chisquared.cc +++ b/libnnpdf/src/chisquared.cc @@ -80,13 +80,11 @@ namespace NNPDF return sqrtmat; } - template - void ComputeChi2(const T* set, int const& nMem, real *const& theory, real *chi2) + // 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) { - matrix const& L = set->GetSqrtCov(); - const double* data = set->GetData(); - const int nDat = set->GetNData(); - // Forward solve Lx = diffs double x[nDat]; for (int n = 0; n < nMem; n++) @@ -101,6 +99,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); } From 51eb1375326c74068eff1cb8a217f27920cdb060 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Thu, 15 Feb 2018 13:02:06 +0100 Subject: [PATCH 07/16] Better handling of very small weights in N_eff:= --- validphys2/src/validphys/reweighting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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): From 65e5401cad36e1d4c9ac2bf5f9a956fe0e4dbd64 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Thu, 15 Feb 2018 13:02:33 +0100 Subject: [PATCH 08/16] Empty ThPredictions constructor --- libnnpdf/src/NNPDF/thpredictions.h | 9 ++- libnnpdf/src/thpredictions.cc | 96 +++++++++++++++++------------- 2 files changed, 63 insertions(+), 42 deletions(-) diff --git a/libnnpdf/src/NNPDF/thpredictions.h b/libnnpdf/src/NNPDF/thpredictions.h index 89050754d3..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 @@ -67,7 +74,7 @@ namespace NNPDF static void Convolute(const PDFSet*,const FKTable*,real*); static void Convolute(const PDFSet*,const FKSet*,real*); static void Convolute(const PDFSet*,const PDFSet*,const FKTable*,real*); - + // Get Methods real* GetObs() const {return fObs;}; //!< Return Obs array real GetObs(int const& idat, int const& imem) const {return fObs[idat*fNpdf + imem];}; diff --git a/libnnpdf/src/thpredictions.cc b/libnnpdf/src/thpredictions.cc index 8be47348b5..f78bce33dd 100644 --- a/libnnpdf/src/thpredictions.cc +++ b/libnnpdf/src/thpredictions.cc @@ -68,7 +68,7 @@ bool ThPredictions::Verbose = true; __m256 t2 = _mm256_hadd_ps(t1,t1); __m128 t3 = _mm256_extractf128_ps(t2,1); __m128 t4 = _mm_add_ss(_mm256_castps256_ps128(t2),t3); - retval = _mm_cvtss_f32(t4); + retval = _mm_cvtss_f32(t4); _mm_empty(); */ @@ -100,7 +100,7 @@ fEtype(pdfset->GetEtype()) { // New theory array fObs = new real[fNData*fNpdf]; - + Timer timer=Timer(); timer.start(); @@ -127,7 +127,7 @@ fEtype(pdfset->GetEtype()) { // New theory array fObs = new real[fNData*fNpdf]; - + Timer timer=Timer(); timer.start(); @@ -150,7 +150,7 @@ fEtype(pdfset->GetEtype()) { // New theory array fObs = new real[fNData*fNpdf]; - + Timer timer=Timer(); timer.start(); @@ -180,7 +180,7 @@ fEtype(pdf1->GetEtype()) // New theory array fObs = new real[fNData*fNpdf]; - + Timer timer=Timer(); timer.start(); @@ -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) { @@ -257,9 +271,9 @@ ThPredictions& ThPredictions::operator+=(const ThPredictions& o) throw EvaluationError("ThPredictions::operator+=","Cannot sum predictions with different numbers of datapoints/PDF members"); if (fEtype != PDFSet::erType::ER_MC && fEtype != PDFSet::erType::ER_MC68 && fEtype != PDFSet::erType::ER_MCT0 && ThPredictions::Verbose == true) - std::cerr << "ThPredictions::operator+= Warning: Observable summation undefined for ErrorTypes other than Monte-Carlo" <(pdf)); return; @@ -423,7 +437,7 @@ void ThPredictions::Convolute(const PDFSet *pdf1, const PDFSet *pdf2, const FKTa if (err != 0) throw RangeError("ThPredictions::Convolute","ThPredictions posix_memalign " + std::to_string(err)); - + // Fetch PDF array GetNZPDF(pdf1, pdf2, fk, pdf); @@ -451,7 +465,7 @@ void ThPredictions::Convolute(const PDFSet *pdf1, const PDFSet *pdf2, const FKTa } MPI::RecvJobs(NData, Npdf, theory); #endif - + // Delete pdfs free(reinterpret_cast(pdf)); return; @@ -463,7 +477,7 @@ void ThPredictions::Convolute(const PDFSet *pdf1, const PDFSet *pdf2, const FKTa void ThPredictions::Convolute(const PDFSet *pdfset, const FKSet *fkset, real* theory) { const int NSigma = fkset->GetNSigma(); - + if (NSigma == 1) ThPredictions::Convolute(pdfset,fkset->GetFK(0), theory ); else @@ -471,22 +485,22 @@ void ThPredictions::Convolute(const PDFSet *pdfset, const FKSet *fkset, real* th SigmaOp op = fkset->GetOperator(); const int NData = fkset->GetNDataFK(); const int Npdf = pdfset->GetMembers(); - + std::vector results; for (int i=0; iGetFK(i), results[i] ); - + op( NData*Npdf, results, theory ); - + for (int i=0; iGetPDF(xgrid[i],Q20,n,&EVLN[i*Nfl]); - + for (int fl=0; flGetPDF(xgrid[i],Q20,n,&EVLN1[i*Nfl]); pdf2->GetPDF(xgrid[i],Q20,n,&EVLN2[i*Nfl]); } - + for (int fl=0; fl Date: Wed, 14 Feb 2018 17:52:46 +0100 Subject: [PATCH 09/16] Separated dataset covariance matrix generation, new function for generation of sqrt matrix --- libnnpdf/src/NNPDF/chisquared.h | 2 + libnnpdf/src/NNPDF/utils.h | 3 -- libnnpdf/src/chisquared.cc | 70 +++++++++++++++++++++++++++++++++ libnnpdf/src/dataset.cc | 35 ++--------------- libnnpdf/src/experiments.cc | 5 +-- libnnpdf/src/utils.cc | 19 --------- 6 files changed, 77 insertions(+), 57 deletions(-) diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index d8c2982685..1d29472f42 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -12,5 +12,7 @@ #include "dataset.h" namespace NNPDF{ + matrix ComputeCovMat(CommonData const& cd, std::vector const& t0); + matrix ComputeSqrtMat(matrix const& inmatrix); template void ComputeChi2(const T*, int const&, real *const&, real *); } diff --git a/libnnpdf/src/NNPDF/utils.h b/libnnpdf/src/NNPDF/utils.h index 77888a5a75..35c3af7837 100644 --- a/libnnpdf/src/NNPDF/utils.h +++ b/libnnpdf/src/NNPDF/utils.h @@ -181,9 +181,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..5555d76f60 100644 --- a/libnnpdf/src/chisquared.cc +++ b/libnnpdf/src/chisquared.cc @@ -5,10 +5,80 @@ // 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(); + if (ndat <= 0) + throw LengthError("ComputeCovMat","invalid number of datapoints!"); + 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) + { + const size_t n = inmatrix.size(0); + if (n <= 0) + throw LengthError("CholeskyDecomposition","attempting a decomposition of an empty matrix!"); + matrix sqrtmat(n,n); + + 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); + if (decomp != 0 ) throw RuntimeException("CholeskyDecomposition", "Error encountered in gsl"); + + 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); + return sqrtmat; + } template void ComputeChi2(const T* 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/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); - } - } From e5b03c9c4df243cfc3521640add9f4c29218d4cf Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Wed, 14 Feb 2018 18:43:14 +0100 Subject: [PATCH 10/16] Add a basic form of the chi2 calculation without using any special types --- libnnpdf/src/NNPDF/chisquared.h | 3 +++ libnnpdf/src/chisquared.cc | 21 +++++++++++++++------ 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index 1d29472f42..addaba1a83 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -14,5 +14,8 @@ 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/chisquared.cc b/libnnpdf/src/chisquared.cc index 5555d76f60..0a672bf740 100644 --- a/libnnpdf/src/chisquared.cc +++ b/libnnpdf/src/chisquared.cc @@ -80,13 +80,11 @@ namespace NNPDF return sqrtmat; } - template - void ComputeChi2(const T* set, int const& nMem, real *const& theory, real *chi2) + // 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) { - matrix const& L = set->GetSqrtCov(); - const double* data = set->GetData(); - const int nDat = set->GetNData(); - // Forward solve Lx = diffs double x[nDat]; for (int n = 0; n < nMem; n++) @@ -101,6 +99,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); } From 2662a27abed17178d8149b34d4db82342df2fc43 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Thu, 15 Feb 2018 13:02:06 +0100 Subject: [PATCH 11/16] Better handling of very small weights in N_eff:= --- validphys2/src/validphys/reweighting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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): From 4fa5732a220ce51fc22f6b3ee5df4c4a743645b6 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Thu, 15 Feb 2018 13:02:33 +0100 Subject: [PATCH 12/16] Empty ThPredictions constructor --- libnnpdf/src/NNPDF/thpredictions.h | 7 +++++++ libnnpdf/src/thpredictions.cc | 14 ++++++++++++++ 2 files changed, 21 insertions(+) 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/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) { From 2bcd6d321efab2fb429e5b879eaa1164714057c5 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Wed, 14 Feb 2018 17:52:46 +0100 Subject: [PATCH 13/16] Separated dataset covariance matrix generation, new function for generation of sqrt matrix --- libnnpdf/src/NNPDF/chisquared.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index addaba1a83..1d29472f42 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -14,8 +14,5 @@ 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 *); } From 35fe396c0c857f2cad7b613e95595b77768ebe46 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Wed, 14 Feb 2018 18:43:14 +0100 Subject: [PATCH 14/16] Add a basic form of the chi2 calculation without using any special types --- libnnpdf/src/NNPDF/chisquared.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index 1d29472f42..addaba1a83 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -14,5 +14,8 @@ 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 *); } From 241956aa4ec52279e8efbd5f5780716ded448dac Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Wed, 7 Mar 2018 01:01:51 +0100 Subject: [PATCH 15/16] Cosmetics --- libnnpdf/src/NNPDF/chisquared.h | 2 +- libnnpdf/src/chisquared.cc | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h index addaba1a83..8c3d28b814 100644 --- a/libnnpdf/src/NNPDF/chisquared.h +++ b/libnnpdf/src/NNPDF/chisquared.h @@ -14,7 +14,7 @@ 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, + 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/chisquared.cc b/libnnpdf/src/chisquared.cc index 0a672bf740..bfd4f336b2 100644 --- a/libnnpdf/src/chisquared.cc +++ b/libnnpdf/src/chisquared.cc @@ -19,13 +19,13 @@ namespace NNPDF { const int ndat = cd.GetNData(); const int nsys = cd.GetNSys(); - if (ndat <= 0) - throw LengthError("ComputeCovMat","invalid number of datapoints!"); + 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; @@ -54,6 +54,7 @@ namespace NNPDF CovMat(i, j) = sig + signor*t0[i]*t0[j]*1e-4; } + } return CovMat; } @@ -81,7 +82,7 @@ namespace NNPDF } // TODO to sort this out, need to make data and theory vectors - void ComputeChi2_basic(int const& nDat, int const& nMem, + void ComputeChi2_basic(int const nDat, int const nMem, const double* data, matrix const& L, real *const& theory, real *chi2) { From c483bad4f7b7b919022f9131fe343bea5663e3e9 Mon Sep 17 00:00:00 2001 From: Nathan Hartland Date: Wed, 7 Mar 2018 12:35:54 +0100 Subject: [PATCH 16/16] gsl matrix view version of SqrtCovMat --- libnnpdf/src/NNPDF/utils.h | 7 ++++--- libnnpdf/src/chisquared.cc | 24 ++++++++++++++---------- 2 files changed, 18 insertions(+), 13 deletions(-) diff --git a/libnnpdf/src/NNPDF/utils.h b/libnnpdf/src/NNPDF/utils.h index 35c3af7837..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 diff --git a/libnnpdf/src/chisquared.cc b/libnnpdf/src/chisquared.cc index bfd4f336b2..021fcaa917 100644 --- a/libnnpdf/src/chisquared.cc +++ b/libnnpdf/src/chisquared.cc @@ -63,21 +63,25 @@ namespace NNPDF const size_t n = inmatrix.size(0); if (n <= 0) throw LengthError("CholeskyDecomposition","attempting a decomposition of an empty matrix!"); - matrix sqrtmat(n,n); - 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)); + 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); - const int decomp = gsl_linalg_cholesky_decomp(mat); - if (decomp != 0 ) throw RuntimeException("CholeskyDecomposition", "Error encountered in gsl"); + // 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) = gsl_matrix_get(mat, i, j); + for (int j = 0; j > i; j++) + sqrtmat(i, j) = 0; - gsl_matrix_free (mat); return sqrtmat; }