diff --git a/CMakeLists.txt b/CMakeLists.txt deleted file mode 100644 index 2dd7b4e4ea..0000000000 --- a/CMakeLists.txt +++ /dev/null @@ -1,109 +0,0 @@ -cmake_minimum_required (VERSION 3.0.2) - -# Use rpaths for now, previously there were issues with osx -SET(CMAKE_SKIP_BUILD_RPATH FALSE) - -# when building, don't use the install RPATH already -# (but later on when installing) -SET(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE) - -SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") - -# Set a default build type if none was specified -if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) - message(STATUS "Setting build type to 'Release' as none was specified.") - set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE) - # Set the possible values of build type for cmake-gui - set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" - "MinSizeRel" "RelWithDebInfo") -endif() - -project(nnpdf) - -set(CMAKE_CXX_STANDARD 14) -set(CMAKE_CXX_STANDARD_REQUIRED ON) -set(CMAKE_CXX_EXTENSIONS OFF) - -set(VERSION "\"4.0\"") -set(nnpdfcpp_VERSION 4.0) - -# check for dependencies -find_package(PythonInterp 3 REQUIRED) -find_package(PkgConfig REQUIRED) -pkg_search_module(LIBARCHIVE REQUIRED libarchive) -pkg_search_module(SQLITE3 REQUIRED sqlite3) -pkg_search_module(GSL REQUIRED gsl) -pkg_search_module(YAML REQUIRED yaml-cpp) - - -option(NNPDF_DEV "n3fit and validphys in developer mode" ON) -set(PROFILE_PREFIX "" CACHE STRING "Where you store the 'data' folder. Default empty uses CMAKE_INSTALL_PREFIX/share/NNPDF.") - -if (PROFILE_PREFIX) - set(PROFILE_PREFIX "${PROFILE_PREFIX}") -else (PROFILE_PREFIX) - set(PROFILE_PREFIX "${CMAKE_INSTALL_PREFIX}/share/NNPDF") -endif() - -# LHAPDF -find_program(LHAPDF_CONFIG lhapdf-config REQUIRED) -if (LHAPDF_CONFIG) - exec_program(${LHAPDF_CONFIG} - ARGS --cflags - OUTPUT_VARIABLE LHAPDF_CXX_FLAGS - ) - set(LHAPDF_CXX_FLAGS ${LHAPDF_CXX_FLAGS} CACHE STRING INTERNAL) - exec_program(${LHAPDF_CONFIG} - ARGS --libs - OUTPUT_VARIABLE LHAPDF_LIBRARIES - ) - set(LHAPDF_LIBRARIES ${LHAPDF_LIBRARIES} CACHE STRING INTERNAL) - #This is to stop the LHAPDF specific warnings from spamming the output - exec_program(${LHAPDF_CONFIG} - ARGS --incdir - OUTPUT_VARIABLE LHAPDF_INCLUDES - ) - set(LHAPDF_INCLUDES ${LHAPDF_INCLUDES} CACHE STRING INTERNAL) - include_directories(SYSTEM ${LHAPDF_INCLUDES}) -endif(LHAPDF_CONFIG) - -# APFEL -find_program(APFEL_CONFIG apfel-config REQUIRED) -if (APFEL_CONFIG) - exec_program(${APFEL_CONFIG} - ARGS --cppflags - OUTPUT_VARIABLE APFEL_CXX_FLAGS - ) - set(APFEL_CXX_FLAGS ${APFEL_CXX_FLAGS} CACHE STRING INTERNAL) - exec_program(${APFEL_CONFIG} - ARGS --ldflags - OUTPUT_VARIABLE APFEL_LIBRARIES - ) - set(APFEL_LIBRARIES ${APFEL_LIBRARIES} CACHE STRING INTERNAL) -endif(APFEL_CONFIG) - -set(DEFAULT_CXX_OPTIONS "-Wall -Wextra -fvisibility-inlines-hidden -fmessage-length=0 -ftree-vectorize -fPIC -fstack-protector-strong -O2 -pipe") - -#strip linker flags to avoid duplication of asan flags -string(REPLACE "-fsanitize=address" "" CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") -string(REPLACE "-fsanitize=address" "" CMAKE_EXE_LINKER_FLAGS_DEBUG "${CMAKE_EXE_LINKER_FLAGS_DEBUG}") - -set(CMAKE_ALL_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${CMAKE_MODULE_LINKER_FLAGS} ${CMAKE_SHARED_LINKER_FLAGS}") - -set(CMAKE_CXX_FLAGS "${DEFAULT_CXX_OPTIONS} ${LHAPDF_CXX_FLAGS} ${APFEL_CXX_FLAGS} ${YAML_CFLAGS} ${SQLITE3_CFLAGS} ${GSL_CFLAGS} ${LIBARCHIVE_CFLAGS}" CACHE STRING "compile flags" FORCE) -set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -g" CACHE STRING "debug compile flags" FORCE) -set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}" CACHE STRING "linker flags" FORCE) -set(CMAKE_EXE_LINKER_FLAGS_DEBUG "${CMAKE_EXE_LINKER_FLAGS_DEBUG}" CACHE STRING "debug linker flags" FORCE) - -# libnnpdf configuration -add_subdirectory(libnnpdf) - -# evolven3fit -add_subdirectory(evolven3fit) -install(FILES ${PROJECT_SOURCE_DIR}/validphys2/src/validphys/datafiles/theory.db DESTINATION ${PROFILE_PREFIX}/) - -if(NNPDF_DEV) - install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} -m pip install --no-deps -e ${PROJECT_SOURCE_DIR})") -else(NNPDF_DEV) - install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} -m pip install --no-deps --ignore-installed ${PROJECT_SOURCE_DIR})") -endif(NNPDF_DEV) diff --git a/conda-recipe/build.sh b/conda-recipe/build.sh deleted file mode 100644 index a8f21fde68..0000000000 --- a/conda-recipe/build.sh +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/bash - -mkdir build -cd build -cmake .. -DCMAKE_INSTALL_PREFIX=${PREFIX} -DNNPDF_DEV=OFF -make -j${CPU_COUNT} -make install diff --git a/conda-recipe/conda_build_config.yaml b/conda-recipe/conda_build_config.yaml index a44b901cc3..95a95fa5b0 100644 --- a/conda-recipe/conda_build_config.yaml +++ b/conda-recipe/conda_build_config.yaml @@ -5,7 +5,3 @@ numpy: pin_run_as_build: lhapdf: x.x.x - apfel: x.x.x.x - gsl: x.x.x - yaml-cpp: x.x.x - libarchive: x.x diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 3e4a615a98..1d2487686d 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -5,29 +5,23 @@ package: source: git_url: ../ +build: + noarch: python + script: {{ PYTHON }} -m pip install . -vv + number: 0 + detect_binary_files_with_prefix: True + # import CONDA_BUILD_SYSROOT from env variable for osx + script_env: + - CONDA_BUILD_SYSROOT # [osx] + requirements: - build: - - {{ compiler("cxx") }} - - {{ compiler("c") }} - - swig - - cmake - - pkg-config - # See https://github.com/NNPDF/nnpdf/pull/1280 - - sysroot_linux-64==2.17 # [linux] - - clang==14 # [osx] - - cmake==3.26 # [osx] host: - - lhapdf # with python wrapper - - sqlite - - gsl # Gsl doesn't link openblas on old debian 7 - - libarchive - - yaml-cpp >=0.6.3 - - apfel >=3 # see https://github.com/scarrazza/apfel - - python - - numpy + - python >=3.9,<3.12 - poetry-core >=1.0.0 - - poetry-dynamic-versioning >=1.1.0 + - poetry-dynamic-versioning + - pip run: + - python >=3.9,<3.12 - tensorflow >=2.10 - psutil # to ensure n3fit affinity is with the right processors - blas==1.0 *mkl* # [osx] # Host's blas is mkl, force also runtime blas to be @@ -35,14 +29,10 @@ requirements: - seaborn - lhapdf - sqlite - - gsl - numpy - - libarchive - - yaml-cpp >=0.6.3 - - apfel - pkg-config - - python # requirements for validphys - reportengine ==0.30.28 # see https://github.com/NNPDF/reportengine + - curio >=1.0 # reportengine uses it but it's not in its dependencies - matplotlib >=3.3.0,<3.8 # see https://github.com/NNPDF/nnpdf/pull/1809 - blessings >=1.7 - scipy >=0.19.1 @@ -50,15 +40,14 @@ requirements: - requests - prompt_toolkit - validobj - - sphinx >=5.0.2,<6 # documentation. Needs pinning temporarily due to markdown - - recommonmark - - sphinx_rtd_theme >0.5 - - sphinxcontrib-bibtex - pineappl >=0.6.2 - eko >=0.14.1 - fiatlux - frozendict # needed for caching of data loading - - curio >=1.0 # reportengine uses it but it's not in its dependencies + - sphinx >=5.0.2,<6 # documentation. Needs pinning temporarily due to markdown + - recommonmark + - sphinx_rtd_theme >0.5 + - sphinxcontrib-bibtex test: requires: @@ -70,14 +59,8 @@ test: source_files: - "*" -build: - number: 0 - detect_binary_files_with_prefix: True - # import CONDA_BUILD_SYSROOT from env variable for osx - script_env: - - CONDA_BUILD_SYSROOT # [osx] about: home: https://nnpdf.mi.infn.it/ - license: SECRET + license: GPLv3.0 summary: "NNPDF analysis framework" diff --git a/conda-recipe/run_test.sh b/conda-recipe/run_test.sh index 762dffc7d3..5be3fb78d7 100644 --- a/conda-recipe/run_test.sh +++ b/conda-recipe/run_test.sh @@ -12,5 +12,5 @@ platformstr=`uname` if [[ "$platformstr" != "Darwin" ]]; then pytest --pyargs n3fit else - echo "Skipping tests on Mac" + echo "Skipping n3fit tests on Mac" fi diff --git a/evolven3fit/CMakeLists.txt b/evolven3fit/CMakeLists.txt deleted file mode 100644 index f4e8d09787..0000000000 --- a/evolven3fit/CMakeLists.txt +++ /dev/null @@ -1,25 +0,0 @@ -# Include files (should this information not be known at this point?) -include_directories(${PROJECT_SOURCE_DIR}/evolven3fit) -include_directories(${PROJECT_SOURCE_DIR}/libnnpdf/src/) -set(EXECUTABLE_OUTPUT_PATH ${CMAKE_BINARY_DIR}/binaries) - -configure_file( - "${PROJECT_SOURCE_DIR}/evolven3fit/evolven3fit.cc.in" - "${PROJECT_SOURCE_DIR}/evolven3fit/evolven3fit.cc" -) - -# Add files to the make -add_executable(evolven3fit ${PROJECT_SOURCE_DIR}/evolven3fit/evolven3fit.cc - ${PROJECT_SOURCE_DIR}/evolven3fit/exportgrid.cc - ${PROJECT_SOURCE_DIR}/evolven3fit/evolgrid.cc ) - -# Set all flags -set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${NNPDF_LDFLAGS} ${GSL_LDFLAGS} ${APFEL_LIBRARIES} ${YAML_LDFLAGS}") - -# I am pretty sure this should not be a thing -string(REPLACE ";" " " CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}") -target_link_libraries(evolven3fit nnpdf ${YAML_LDFLAGS} ${APFEL_LIBRARIES} ${GSL_LDFLAGS}) - -install(TARGETS evolven3fit DESTINATION bin - PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ -GROUP_EXECUTE WORLD_READ WORLD_EXECUTE) diff --git a/evolven3fit/evolgrid.cc b/evolven3fit/evolgrid.cc deleted file mode 100644 index a25b455e6b..0000000000 --- a/evolven3fit/evolgrid.cc +++ /dev/null @@ -1,310 +0,0 @@ -// $Id$ -// -// NNPDF++ 2012-2015 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -#include -#include -#include -#include -#include "evolgrid.h" -#include "exportgrid.h" -#include -#include -#include -#include -#include - -using namespace NNPDF; -using namespace std; - -//______________________________________________________ -vector generate_q2grid(const int nq2, - const double q2min, - const double q2max, - const double lambda2) -{ - vector q2grid; - const double tau_min = log( log( q2min / lambda2 )); - const double tau_max = log( log( q2max / lambda2 )); - const double delta_tau = (tau_max - tau_min) /( (double) nq2 - 1); - - for (int iq2 = 0; iq2 < nq2; iq2++) - q2grid.push_back(lambda2 * exp(exp( tau_min + iq2*delta_tau))); - return q2grid; -} - -//______________________________________________________ -vector< vector > generate_q2subgrids( const int nq2, - const double q2min, const double q2max) -{ - const double eps = 1E-4; - const double lambda2 = 0.0625e0; - const int nfpdf = APFEL::GetMaxFlavourPDFs(); - - const std::array hqth = {{ - pow(APFEL::GetThreshold(4),2), - pow(APFEL::GetThreshold(5),2), - pow(APFEL::GetThreshold(6),2), - }}; - - // Determine subgrid edges - vector subgrid_edges={q2min}; - for (int i=0; i<(nfpdf-3); i++) - if (hqth[i] > q2min) subgrid_edges.push_back(hqth[i]); - subgrid_edges.push_back(q2max); - - // Determine point distribution across subgrids and generate subgrids - const std::vector point_dist = generate_q2grid(nq2, q2min, q2max, lambda2); - std::vector> q2grid_subgrids; - for (size_t i=0; i -eps && (q2pt - max_edge) < eps;}; - const int npoints_sub = std::count_if(point_dist.begin(), point_dist.end(), in_sub); - - if (npoints_sub < 2) - throw std::runtime_error("Too few points to sufficiently populate subgrids. More Q points required"); - - const vector< double > subgrid = generate_q2grid(npoints_sub, min_edge, max_edge, lambda2); - q2grid_subgrids.push_back(subgrid); - } - - return q2grid_subgrids; -} - -//______________________________________________________ -vector compute_pids(int maxnf, bool qed) -{ - vector pids; - for (int i = -maxnf; i <= maxnf; i++) - { - int fl = i; - if (i == 0) fl = 21; - pids.push_back(fl); - } - if (qed) - pids.push_back(22); - return pids; -} - -//______________________________________________________ -vector convert_pids_to_indexes(vector const& pids) -{ - vector indexes; - for (auto fl: pids) - { - if (fl == 21) - indexes.push_back(7); - else if (fl == 22) - indexes.push_back(0); - else - indexes.push_back(fl+7); // follows APFEL ordering from -7:6 where -7 is the photon. - } - return indexes; -} - -//______________________________________________________ -EvolveGrid::EvolveGrid(vector const& initialscale_grid, - map const& theory): - fnq2(50), - fnf(max(stoi(theory.at("MaxNfPdf")), stoi(theory.at("MaxNfAs")))), - fqed(stoi(theory.at("QED"))), - fq2min(pow(stof(theory.at("Q0")), 2)), - fq2max(1E5*1E5), - finitialscale_grid(initialscale_grid) -{ - const auto xg = initialscale_grid[0].GetXGrid(); - cout << "- Initialising evolution with " + std::to_string(xg.size()) + " x-points" << endl; - - // Initialize APFEL - APFEL::SetParam(theory); - - // reset ren/fac to 1 if EScaleVar is 0 - if (!stoi(theory.at("EScaleVar"))) - APFEL::SetRenFacRatio(1.0); - - // Fetch initial-scale x-grid - APFEL::SetNumberOfGrids(1); - APFEL::SetExternalGrid(1, xg.size()-1, 3, (double*) xg.data()); - - // General settings - APFEL::SetQLimits(sqrt(fq2min), sqrt(fq2max)); - APFEL::SetFastEvolution(false); - APFEL::EnableEvolutionOperator(true); - - APFEL::InitializeAPFEL(); -} - -//______________________________________________________ -void EvolveGrid::WriteInfoFile(string const& infofile, int nrep) const -{ - // skip if file exists - struct stat s; - stringstream infodata; - if(stat(infofile.c_str(), &s) == 0 && S_ISREG(s.st_mode)) - { - cout << "- LHAPDF info already exported: " << infofile << endl; - return; - } - else - cout << "- Exporting LHAPDF info file: " << infofile << endl; - - // LHAPDF6 HEADER - const auto xg = finitialscale_grid[0].GetXGrid(); - const auto q2subgrids = generate_q2subgrids(fnq2, fq2min, fq2max); - - infodata << "SetDesc: \"NNPDF x.x\"" << endl; - infodata << "SetIndex: " << endl; - infodata << "Authors: NNPDF Collaboration." << endl; - infodata << "Reference: arXiv:xxxx.xxxxxx" << endl; - infodata << "Format: lhagrid1" << endl; - infodata << "DataVersion: 1" << endl; - if (nrep == -1) - infodata << "NumMembers: REPLACE_NREP" << endl; - else - infodata << "NumMembers: " << nrep << endl; - infodata << "Particle: 2212" << endl; - infodata << "Flavors: ["; - for (auto fl: compute_pids(fnf, fqed)) - infodata << fl << (fl == 22 ? "" : ", "); - infodata << "]" << endl; - infodata << "OrderQCD: " << APFEL::GetPerturbativeOrder() << endl; - infodata << "FlavorScheme: variable" << endl; - infodata << "NumFlavors: " << std::max(APFEL::GetMaxFlavourPDFs(), APFEL::GetMaxFlavourAlpha()) << endl; - infodata << "ErrorType: replicas" << endl; - - infodata.precision(7); - infodata << scientific; - infodata << "XMin: "<< xg[0] << endl; - infodata << "XMax: "<< xg.back() << endl; - infodata << "QMin: "<< sqrt(fq2min) << endl; - infodata << "QMax: "<< sqrt(fq2max) << endl; - infodata << "MZ: " << APFEL::GetZMass() << endl; - infodata << "MUp: 0\nMDown: 0\nMStrange: 0" << std::endl; - infodata << "MCharm: " << APFEL::GetThreshold(4) << endl; - infodata << "MBottom: " << APFEL::GetThreshold(5) << endl; - infodata << "MTop: " << APFEL::GetThreshold(6) << endl; - infodata << fixed << "AlphaS_MZ: " << APFEL::AlphaQCD(APFEL::GetZMass()) << endl; - infodata << scientific; - infodata << "AlphaS_OrderQCD: " << APFEL::GetPerturbativeOrder() << endl; - infodata << "AlphaS_Type: ipol" << endl; - infodata << "AlphaS_Qs: ["; - - for (int s = 0; s < (int) q2subgrids.size(); s++) - for (int iq = 0; iq < (int) q2subgrids[s].size(); iq++) - infodata << sqrt(q2subgrids[s][iq]) << ((s == (int) q2subgrids.size()-1 && iq == (int) q2subgrids[s].size()-1) ? "]\n" : ", "); - - infodata << "AlphaS_Vals: ["; - for (int s = 0; s < (int) q2subgrids.size(); s++) - for (int iq = 0; iq < (int) q2subgrids[s].size(); iq++) - infodata << APFEL::AlphaQCD(sqrt(q2subgrids[s][iq])) << ((s == (int) q2subgrids.size()-1 && iq == (int) q2subgrids[s].size()-1) ? "]\n" : ", "); - - infodata << "AlphaS_Lambda4: 0.342207" << endl; - infodata << "AlphaS_Lambda5: 0.239" << endl; - write_to_file(infofile, infodata.str()); -} - -//______________________________________________________ -vector EvolveGrid::WriteLHAFile() const -{ - const auto xg = finitialscale_grid[0].GetXGrid(); - vector outstream(finitialscale_grid.size()); - const vector pids = compute_pids(fnf, fqed); - - for (size_t i = 0; i < finitialscale_grid.size(); i++) - { - outstream[i] << scientific << setprecision(7); - outstream[i] << "PdfType: replica\nFormat: lhagrid1\nFromMCReplica: " << finitialscale_grid[i].GetReplica() << "\n---" << std::endl; - } - - // compute q2 subgrids - const auto q2subgrids = generate_q2subgrids(fnq2, fq2min, fq2max); - for (auto subgrid: q2subgrids) - { - for (auto& stream: outstream) - { - // Print out x-grid - for ( auto x : xg ) - stream << x << " "; - stream << std::endl; - - // Print out q2-grid - for ( auto q2val : subgrid ) - stream << sqrt(q2val) << " "; - stream << std::endl; - - // Print out final-state PIDs - for (auto fl: pids) - stream << fl << " "; - stream << std::endl; - } - - // Compute Evolution Operators - // Ordered evol_op[q][ox][ofl][ix][ifl] - // This is a big array - const int evol_op_size = subgrid.size()*xg.size()*xg.size()*14*14; - vector evol_op(evol_op_size, std::numeric_limits::signaling_NaN()); - - for ( size_t iq = 0; iq < subgrid.size(); iq++ ) - { - const double qi = sqrt(fq2min); - const double qt = sqrt(subgrid[iq]); - - std::cout << " Qi = " << qi << " -> Qt = " << qt << " GeV" << std::endl; - std::cout << "-------------------------"< -#include -#include -#include -class ExportGrid; - -/** - * @brief The EvolveGrid class. - * Allocates APFEL evolution once for all replicas. - * Computes the evolution operator and writes info and LHA files. - */ -class EvolveGrid -{ -public: - /** - * @brief Class constructor. - * @param initialscale_grids vector with initial grids per replica - * @param theory the theory mapping for APFEL - */ - EvolveGrid(std::vector const& initialscale_grid, - std::map const& theory); - - /** - * @brief Export LHAPDF info file, if infofile does not exist. - * @param infofile the desired file path - * @int nrep if not set (ie. set to -1) will print REPLACE_NREP in the info file - * if specified (ie. != -1) will print nrep in the info file. - */ - void WriteInfoFile(std::string const& infofile, int nrep = -1) const; - - /** - * @brief Export the LHA evolved grid file. - * @param replica_file the desired file path - * @param rep replica number, by default =0 (for nnfit) while can take - * larger numbers for revolve/evolvefit programs. - */ - std::vector WriteLHAFile() const; - -private: - const int fnq2; - const int fnf; - const bool fqed; - const double fq2min; - const double fq2max; - std::vector const& finitialscale_grid; -}; - -/** - * @brief Generates a list of nq2 Q2 points between q2min and q2max. - * Distributed as linear in tau = log( log( Q2 / lambda2 ) ) - * @param nq2 number of q2 points (nodes). - * @param q2min minimum q2 value for the grid. - * @param q2max maximum q2 value for the grid. - * @param lambda2 coefficient for the distribution density. - * @return a vector with the q2grid nodes - */ -std::vector generate_q2grid(const int nq2, - const double q2min, - const double q2max, - const double lambda2); - -/** - * @brief Compute Q2 grid. - * This is suitable for PDFs not AlphaS (see the nfpdf variable) - * @param nq2 the number of q2 points for the grid. - * @param q2min minimum q2 value for the grid. - * @param q2max maximum q2 value for the grid. - * @return the q2 subgrid vector - */ -std::vector> generate_q2subgrids(const int nq2, - const double q2min, - const double q2max); diff --git a/evolven3fit/evolven3fit.cc.in b/evolven3fit/evolven3fit.cc.in deleted file mode 100644 index e01eee48d1..0000000000 --- a/evolven3fit/evolven3fit.cc.in +++ /dev/null @@ -1,145 +0,0 @@ -// $Id$ -// -// NNPDF++ 2016 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it - -#include -#include -#include -#include - -#include -#include -#include -#include - -#include "exportgrid.h" -#include "evolgrid.h" - -using namespace NNPDF; -using std::cout; -using std::endl; -using std::cerr; -using std::string; -using std::stringstream; -using std::stoi; - -#define DBPATH "@PROFILE_PREFIX@/theory.db" - -// Check if folder exists -bool CheckConsistency(string const& folder, string const& exportfile) -{ - bool status1 = false, status2 = false; - struct stat s, t; - if (stat(folder.c_str(), &s) == 0) - if (s.st_mode & S_IFDIR) - status1 = true; - if (stat(exportfile.c_str(), &t) == 0) - if (t.st_mode) - status2 = true; - if (status1 == true && status2 == true) return true; - else return false; -} - -// return theoryid from runcard -int get_theory_id_from_runcard(string const& filteryaml_path) -{ - YAML::Node runcard; - try { - runcard = YAML::LoadFile(filteryaml_path); - } - catch(YAML::BadFile &e) - { - throw FileError("evolven3fit", "runcard not found: " + filteryaml_path); - } - return runcard["theory"]["theoryid"].as(); -} - -/** - * This program: - * - takes as input a fit folder and a theoryID, - * - loads a vector of ExportGrid for all replicas generated by nnfit, - * - computes the DGLAP evolution operators for the theoryID - * - applies the evolution operators to the ExportGrid objects - * - outputs the evolved PDFs in the LHAPDF format to the fit folder. - */ -int main(int argc, char **argv) -{ - // Read configuration filename from arguments - if (argc != 3 && argc != 5) - { - cerr << "\nusage: evolven3fit [configuration folder] [max_replicas] [--theory_id ID]\n" << endl; - exit(EXIT_FAILURE); - } - - const string fit_path = argv[1]; - const int maxreplica = stoi(argv[2]); - const string filteryaml_path = fit_path + "/filter.yml"; - int theory_id; - - // Get theory id - if (argc == 5) - { - const string flag = argv[3]; - if (flag == "--theory_id") - theory_id = stoi(argv[4]); - else - { - cerr << "\nNot supported flag: " << flag << endl; - exit(EXIT_FAILURE); - } - } - else - theory_id = get_theory_id_from_runcard(filteryaml_path); - cout << "Theory ID = " << theory_id << endl; - - // load theory from db - std::map theory_map; - NNPDF::IndexDB db(DBPATH, "theoryIndex"); - auto keys = APFEL::kValues; - keys.push_back("EScaleVar"); - db.ExtractMap(theory_id, keys, theory_map); - - // load grids - vector initialscale_grids; - vector replicas; - for (int nrep = 1; nrep <= maxreplica; nrep++) - { - const string folder = fit_path + "/nnfit/replica_" + std::to_string(nrep); - const string path = folder + "/" + fit_path + ".exportgrid"; - bool status = CheckConsistency(folder, path); - if (status) - { - initialscale_grids.emplace_back(path); - replicas.push_back(nrep); - } - else - { - cout << "Skipping exportgrid (missing file): " << path << endl; - } - } - - if (initialscale_grids.size() == 0) - throw NNPDF::RuntimeException("main", "nrep = 0, check replica folder/files."); - - string infofile = fit_path + "/nnfit/" + fit_path + ".info"; - auto dglapg = EvolveGrid(initialscale_grids, theory_map); - dglapg.WriteInfoFile(infofile); - - const auto outstream = dglapg.WriteLHAFile(); - for (size_t i = 0; i < outstream.size(); i++) - { - stringstream replica_file; - replica_file << fit_path - << "/nnfit/replica_" - << replicas[i] - << "/" - << fit_path - << ".dat"; - write_to_file(replica_file.str(), outstream[i].str()); - } - - return 0; -} diff --git a/evolven3fit/exportgrid.cc b/evolven3fit/exportgrid.cc deleted file mode 100644 index 6e2aaedabb..0000000000 --- a/evolven3fit/exportgrid.cc +++ /dev/null @@ -1,203 +0,0 @@ -#include "exportgrid.h" - -#include -#include - -#include - -#include -#include - -using namespace std; -using NNPDF::PDFSet; -using NNPDF::RuntimeException; -using NNPDF::write_to_file; - -// Grid upon which fitted PDFs are exported -// Defined via the 'applgrid' prescription as -// linear in y(x) = ln(1/x) + a(1-x) with a = 30 -// Limits: [1E-9, 1] -// 197 points (maximum that APFEL can handle) -const std::vector export_xgrid = -{{ - 1.000000000000000e-09, 1.297084823439570e-09, 1.682429034742569e-09, - 2.182253154205826e-09, 2.830567417398188e-09, 3.671485978929410e-09, - 4.762228629353150e-09, 6.177014273761803e-09, 8.012111098984379e-09, - 1.039238706072454e-08, 1.347980640738050e-08, 1.748445036917782e-08, - 2.267881188811028e-08, 2.941633703008346e-08, 3.815547465958784e-08, - 4.949087072321288e-08, 6.419382957083709e-08, 8.326479519868589e-08, - 1.080014229938285e-07, 1.400868730811297e-07, 1.817043317937722e-07, - 2.356855515453774e-07, 3.057035125953233e-07, 3.965223098417466e-07, - 5.143212572365697e-07, 6.671152451366758e-07, 8.652999229731433e-07, - 1.122358752414873e-06, 1.455779955476825e-06, 1.888245605146129e-06, - 2.449173524549460e-06, 3.176716500287171e-06, 4.120354152327973e-06, - 5.344252657520903e-06, 6.931618978063155e-06, 8.990342582381449e-06, - 1.166030301122581e-05, 1.512283122887690e-05, 1.961295293492122e-05, - 2.543522071345024e-05, 3.298416834359921e-05, 4.277070539720159e-05, - 5.545612481058487e-05, 7.189583136325140e-05, 9.319542279796139e-05, - 1.207823677313300e-04, 1.564972094665545e-04, 2.027089363284954e-04, - 2.624597993319508e-04, 3.396452441689850e-04, 4.392344430004219e-04, - 5.675356601045333e-04, 7.325076157255367e-04, 9.441121054524513e-04, - 1.214693176869783e-03, 1.559353061182245e-03, 1.996274511413378e-03, - 2.546914937365516e-03, 3.235975102131256e-03, 4.091034365095647e-03, - 5.141759770839620e-03, 6.418650960623169e-03, 7.951379403063506e-03, - 9.766899996240997e-03, 1.188761392513640e-02, 1.432989476439189e-02, - 1.710322794602712e-02, 2.021007339250794e-02, 2.364639713695425e-02, - 2.740269157283572e-02, 3.146525061324443e-02, 3.581748292824286e-02, - 4.044110601633167e-02, 4.531713439738071e-02, 5.042663479500688e-02, - 5.575126100843393e-02, 6.127360193905193e-02, 6.697738294982548e-02, - 7.284755899865170e-02, 7.887033222927267e-02, 8.503311978014517e-02, - 9.132449102786790e-02, 9.773408797837715e-02, 1.042525382086388e-01, - 1.108713665472371e-01, 1.175829093728782e-01, 1.243802338015993e-01, - 1.312570629450312e-01, 1.382077077072888e-01, 1.452270051356506e-01, - 1.523102630659852e-01, 1.594532106521559e-01, 1.666519542939869e-01, - 1.739029384555782e-01, 1.812029108733327e-01, 1.885488916790972e-01, - 1.959381459991933e-01, 2.033681596297647e-01, 2.108366174291031e-01, - 2.183413841065613e-01, 2.258804871240649e-01, 2.334521014595030e-01, - 2.410545360116810e-01, 2.486862214527622e-01, 2.563456993587234e-01, - 2.640316124686842e-01, 2.717426959427826e-01, 2.794777695041488e-01, - 2.872357303648326e-01, 2.950155468476644e-01, 3.028162526268661e-01, - 3.106369415195031e-01, 3.184767627680818e-01, 3.263349167616716e-01, - 3.342106511491565e-01, 3.421032573036267e-01, 3.500120671016855e-01, - 3.579364499855710e-01, 3.658758102796432e-01, 3.738295847359622e-01, - 3.817972402864939e-01, 3.897782719819471e-01, 3.977722010992863e-01, - 4.057785734023404e-01, 4.137969575406706e-01, 4.218269435745480e-01, - 4.298681416141745e-01, 4.379201805632053e-01, 4.459827069569899e-01, - 4.540553838875619e-01, 4.621378900076507e-01, 4.702299186071416e-01, - 4.783311767556753e-01, 4.864413845060586e-01, 4.945602741533477e-01, - 5.026875895451769e-01, 5.108230854390865e-01, 5.189665269032351e-01, - 5.271176887569979e-01, 5.352763550484283e-01, 5.434423185656607e-01, - 5.516153803797675e-01, 5.597953494166407e-01, 5.679820420558005e-01, - 5.761752817540883e-01, 5.843748986924983e-01, 5.925807294444404e-01, - 6.007926166639503e-01, 6.090104087923975e-01, 6.172339597824495e-01, - 6.254631288380691e-01, 6.336977801694852e-01, 6.419377827620891e-01, - 6.501830101583613e-01, 6.584333402519444e-01, 6.666886550930888e-01, - 6.749488407047076e-01, 6.832137869083856e-01, 6.914833871596969e-01, - 6.997575383922505e-01, 7.080361408699164e-01, 7.163190980467328e-01, - 7.246063164340254e-01, 7.328977054742707e-01, 7.411931774214037e-01, - 7.494926472270083e-01, 7.577960324322238e-01, 7.661032530649272e-01, - 7.744142315419215e-01, 7.827288925758362e-01, 7.910471630864785e-01, - 7.993689721163776e-01, 8.076942507502913e-01, 8.160229320384573e-01, - 8.243549509233821e-01, 8.326902441699869e-01, 8.410287502988437e-01, - 8.493704095226000e-01, 8.577151636849855e-01, 8.660629562026835e-01, - 8.744137320097212e-01, 8.827674375042057e-01, 8.911240204974589e-01, - 8.994834301652264e-01, 9.078456170010206e-01, 9.162105327713991e-01, - 9.245781304731123e-01, 9.329483642920292e-01, 9.413211895637338e-01, - 9.496965627357548e-01, 9.580744413312983e-01, 9.664547839144387e-01, - 9.748375500567046e-01, 9.832227003049778e-01, 9.916101961506623e-01, - 1.000000000000000e+00 - }}; - - -ExportGrid::ExportGrid(PDFSet const& pdf, - const int imem, - const int irep, - const double Q20): - fRep(irep), - fQ20(Q20), - fXgrid(export_xgrid), - fLabels(pdf.fl_labels()), - fPDFgrid() -{ - // Populate initial scale x-grid - for (auto x : export_xgrid) - { - array evl{}; - array lha{}; - pdf.GetPDF(x, fQ20, imem, evl.data()); - PDFSet::EVLN2LHA(evl.data(), lha.data()); - - // Cast to double for export - array lha_double{}; - std::copy(lha.begin(), lha.end(), lha_double.begin()); - - fPDFgrid.push_back(lha_double); - } -} - -ExportGrid::ExportGrid(string filename): - ExportGrid(YAML::LoadFile(filename)) -{ - std::cout << "Read ExportGrid from: " +filename <()), - fQ20(input["q20"].as()), - fXgrid(input["xgrid"].as>()), - fLabels(input["labels"].as>()), - fPDFgrid() -{ - for (auto xfx : input["pdfgrid"]) - { - std::array input = xfx.as>(); - if (fLabels != PDFSet::fl_labels()) - { - std::array output; - for (size_t fl1 = 0; fl1 < PDFSet::fl_labels().size(); fl1++) - for (size_t fl2 = 0; fl2 < fLabels.size(); fl2++) - if (PDFSet::fl_labels()[fl1] == fLabels[fl2]) - { - output[fl1] = input[fl2]; - break; - } - fPDFgrid.push_back(output); - } - else - fPDFgrid.push_back(input); - } - if (fPDFgrid.size() != fXgrid.size()) - throw RuntimeException("ExportGrid::ExportGrid", "Mismatch in x-grid and number of PDFgrid entries"); -} - -// Because yaml-cpp doesn't like scientific notation -std::string ExportGrid::Format(double in) -{ - std::stringstream cvrt; - cvrt << scientific << setprecision(14) << in; - return cvrt.str(); -} - -void ExportGrid::Write(const std::string filename) const -{ - YAML::Emitter data; - data.SetOutputCharset(YAML::EscapeNonAscii); - - // Convert x-grid to formatted strings - std::vector xg_str; - std::transform(fXgrid.begin(), fXgrid.end(), std::back_inserter(xg_str), - [](double in) -> std::string { return ExportGrid::Format(in); }); - - data << YAML::BeginMap; - data << YAML::Key << "replica" << YAML::Value << fRep; - data << YAML::Key << "q20" << YAML::Value << Format(fQ20); - data << YAML::Key << "xgrid" << YAML::Value; - data << YAML::Flow << xg_str; - - // Yaml-cpp doesn't like arrays - const std::vector labels(fLabels.begin(), fLabels.end()); - data << YAML::Key << "labels"; - data << YAML::Value; - data << YAML::Flow; - data << labels; - - // Print pdf grid - data << YAML::Key << "pdfgrid" << YAML::Value; - data << YAML::BeginSeq; - for (auto pdf: fPDFgrid) - { - // Convert xf-grid to formatted strings - std::vector xfg_str; - std::transform(pdf.begin(), pdf.end(), std::back_inserter(xfg_str), - [](double in) -> std::string { return ExportGrid::Format(in); }); - - data << YAML::Flow << xfg_str; - } - - data << YAML::EndSeq; - - data << YAML::EndMap; - cout << "- Printing grid to file: " << filename << endl; - write_to_file(filename, data.c_str()); -} - diff --git a/evolven3fit/exportgrid.h b/evolven3fit/exportgrid.h deleted file mode 100644 index 4adbc7ff29..0000000000 --- a/evolven3fit/exportgrid.h +++ /dev/null @@ -1,55 +0,0 @@ -// exportgrid.h -// Definition of initial scale export PDF grid -#include -#include -#include -#include - -using std::vector; -using std::array; -using NNPDF::real; - -/** - * This class writes to a YAML file the values of a PDFSet - * (for a given replica number and input scale) in a grid of x points - * for all 14 flavours enumerated in the PDFSet class. - */ -class ExportGrid -{ -public: - // Generate an ExportGrid from a PDFset - ExportGrid(NNPDF::PDFSet const& pdf, // The PDFset generating the ExportGrid - const int imem, // The member of `pdf` to be exported - const int irep, // The replica ID of `imem` from `pdf` - const double Q0); // The scale at which to export - - // Read an ExportGrid from file - ExportGrid(std::string filename); - ExportGrid(YAML::Node input); - - // Print exportgrid to file - void Write(const std::string filename) const; - - // Return the initial-scale x-grid - vector GetXGrid() const {return fXgrid;} - - // Return the initial-scale x-grid - double GetPDF(int ix, int ifl) const {return fPDFgrid[ix][ifl];} - - // Get the replica ID - int GetReplica() const {return fRep;} - - // Get PDF grid - vector> const& GetPDFgrid() const { return fPDFgrid; } - - // Set PDF grid - void SetPDFgrid(vector> const& pgrid) { fPDFgrid = pgrid; } - -private: - static std::string Format(double in); - const int fRep; - const double fQ20; - const vector fXgrid; - const array fLabels; - vector> fPDFgrid; -}; diff --git a/evolven3fit/varflavors.py b/evolven3fit/varflavors.py deleted file mode 100644 index 41e786907b..0000000000 --- a/evolven3fit/varflavors.py +++ /dev/null @@ -1,51 +0,0 @@ -""" -varflavors.py - -When producing a PDF with different maximal flavor number than the nf during the fit, we first run -evolven3fit by manually selecting a theory id of the theory correspodning to the desired nf -evolution. - -After running evolven3fit, but before running posfit, this script should be run to replace the -``AlphaS_MZ'' and ``MZ'' values in the .info file, with the ``alphas'' and ``Qref'' valuse from the -theory database. -""" - -from argparse import ArgumentParser -from pathlib import Path - -from validphys.loader import Loader -from validphys.theorydbutils import fetch_theory - - -ll = Loader() -path_db = ll.theorydb_file - -if __name__ == "__main__": - parser = ArgumentParser() - parser.add_argument( - "--fit_folder", required=True, help=f"Path to the folder containing the evolved fit" - ) - parser.add_argument( - "--theory_id", required=True, help=f"ID of the theory used to evolve the fit" - ) - args = parser.parse_args() - - path_fit_folder = Path(args.fit_folder) - fit_name = path_fit_folder.name - path_info_file = path_fit_folder / f"nnfit/{fit_name}.info" - - theory_dict = fetch_theory(path_db, args.theory_id) - Qref = theory_dict["Qref"] - alphas = theory_dict["alphas"] - - path_temp_info = path_info_file.parent / "temp.info" - with open(path_info_file, "r+") as f: - with open(path_temp_info, "w") as tempfile: - for line in f: - if line.startswith("AlphaS_MZ:"): - line = f"AlphaS_MZ: {alphas}\n" - if line.startswith("MZ"): - line = f"MZ: {Qref}\n" - tempfile.write(line) - - path_temp_info.rename(path_info_file) diff --git a/libnnpdf/CMakeLists.txt b/libnnpdf/CMakeLists.txt deleted file mode 100644 index 9fbf64e4ee..0000000000 --- a/libnnpdf/CMakeLists.txt +++ /dev/null @@ -1,58 +0,0 @@ -set(prefix ${CMAKE_INSTALL_PREFIX}) -set(exec_prefix "${prefix}") -set(includedir "${prefix}/include") -set(libdir "${prefix}/lib") - -configure_file( - "${PROJECT_SOURCE_DIR}/libnnpdf/src/NNPDF/config.h.in" - "${PROJECT_SOURCE_DIR}/libnnpdf/src/NNPDF/config.h" - ) - -configure_file( - "${PROJECT_SOURCE_DIR}/libnnpdf/scripts/nnpdf.pc.in" - "${PROJECT_SOURCE_DIR}/libnnpdf/scripts/nnpdf.pc" - ) - -configure_file( - "${PROJECT_SOURCE_DIR}/libnnpdf/src/NNPDF/common.h.in" - "${PROJECT_SOURCE_DIR}/libnnpdf/src/NNPDF/common.h" - ) - -# add preprocessor flag -add_definitions(-DDEFAULT_NNPDF_PROFILE_PATH="${prefix}/share/NNPDF/nnprofile.yaml") - -# Note that BEFORE is important here: Otherwise we might be reading -# the headers from a previous installation. -include_directories(BEFORE src/NNPDF src) -include_directories(${GSL_INCLUDE_DIRS} ${YAML_INCLUDE_DIRS}) -FILE(GLOB_RECURSE Headers "src/NNPDF/*.h") -add_library(nnpdf SHARED src/common.cc - src/commondata.cc - src/chisquared.cc - src/dataset.cc - src/experiments.cc - src/fastkernel.cc - src/fkgenerator.cc - src/fkset.cc - src/lhapdfset.cc - src/logger.cc - src/nnmpi.cc - src/nnpdfdb.cc - src/parametrisation.cc - src/pdfset.cc - src/positivity.cc - src/positivity.cc - src/randomgenerator.cc - src/thpredictions.cc - src/utils.cc - src/pathlib.cc - ${Headers} -) - -target_link_libraries(nnpdf ${LHAPDF_LIBRARIES} ${GSL_LDFLAGS} ${SQLITE3_LDFLAGS} ${LIBARCHIVE_LDFLAGS} ${YAML_LDFLAGS}) - -install(FILES ${PROJECT_SOURCE_DIR}/libnnpdf/scripts/nnpdf.pc DESTINATION lib/pkgconfig) -install(DIRECTORY src/NNPDF DESTINATION include) -file(WRITE ${PROJECT_SOURCE_DIR}/libnnpdf/REAMDE.md "The share folder of NNPDF has been moved, see: https://github.com/NNPDF/nnpdf/pull/1861\n") -install(FILES ${PROJECT_SOURCE_DIR}/libnnpdf/REAMDE.md DESTINATION share/NNPDF) -install(TARGETS nnpdf DESTINATION lib) diff --git a/libnnpdf/nnprofile.yaml.in b/libnnpdf/nnprofile.yaml.in deleted file mode 100644 index e3f5486208..0000000000 --- a/libnnpdf/nnprofile.yaml.in +++ /dev/null @@ -1,36 +0,0 @@ -# Local resource locations -data_path: '@PROFILE_PREFIX@' - -# Remote resource locations -fit_urls: - - 'https://data.nnpdf.science/fits/' - - 'https://nnpdf.web.cern.ch/nnpdf/fits/' - -fit_index: 'fitdata.json' - -hyperscan_urls: - - 'https://data.nnpdf.science/hyperscans/' - -hyperscan_index: 'hyperscandata.json' - -theory_urls: - - 'https://nnpdf.web.cern.ch/nnpdf/tables/' - -theory_index: 'theorydata.json' - -lhapdf_urls: - - 'http://lhapdfsets.web.cern.ch/lhapdfsets/current/' -nnpdf_pdfs_urls: - - 'https://data.nnpdf.science/pdfs/' -nnpdf_pdfs_index: 'pdfdata.json' - -#Server side uploading locations -upload_host: 'nnpdf@vp.nnpdf.science' -reports_target_dir: "validphys-reports/" -reports_root_url: 'https://vp.nnpdf.science/' -fits_target_dir: "WEB/fits/" -fits_root_url: 'https://data.nnpdf.science/fits/' -pdfs_target_dir: "WEB/pdfs/" -pdfs_root_url: 'https://data.nnpdf.science/pdfs/' -hyperscan_target_dir: "WEB/hyperscans/" -hyperscan_root_url: 'https://data.nnpdf.science/hyperscans/' diff --git a/libnnpdf/scripts/nnpdf.pc.in b/libnnpdf/scripts/nnpdf.pc.in deleted file mode 100644 index ef494719d0..0000000000 --- a/libnnpdf/scripts/nnpdf.pc.in +++ /dev/null @@ -1,10 +0,0 @@ -prefix=@prefix@ -exec_prefix=@exec_prefix@ -includedir=@includedir@ -libdir=@libdir@ - -Name: nnpdf -Description: The nnpdf library -Version: @VERSION@ -Cflags: -I@includedir@ -std=c++14 -Libs: -L@libdir@ -lnnpdf diff --git a/libnnpdf/src/NNPDF/chisquared.h b/libnnpdf/src/NNPDF/chisquared.h deleted file mode 100644 index 2e23a7acc6..0000000000 --- a/libnnpdf/src/NNPDF/chisquared.h +++ /dev/null @@ -1,39 +0,0 @@ -// $Id -// -// NNPDF++ 2012-2015 -// -// Authors: Nathan Hartland, n.p.hartland@vu.nl -// Stefano Carrazza, stefano.carrazza@mi.infn.it - -#pragma once - -#include "common.h" -#include "experiments.h" -#include "dataset.h" - -namespace NNPDF{ - matrix ComputeCovMat_basic(int const nDat, - int const nSys, - std::vector const& sqrt_weights, - std::vector const& central_values, - std::vector const& stat_error, - sysError** const systematic_errors, - bool const mult_errors, // account for multiplicative uncertainties in building the CovMat - bool const use_theory_errors, // account for MC uncertainties in building the CovMat - bool const th_cov_matrix, // account for theoretical uncertainties in building the CovMat - std::string filename, - std::vector bmask); - - matrix ComputeCovMat(CommonData const& cd, std::vector const& t0, - const bool th_cov_matrix = false, - std::string filename = "", - std::vector bmask = {}, - double weight=1.); - matrix ComputeSqrtMat(matrix const& inmatrix); - matrix read_theory_covmat(int ndata, const std::string filename, std::vector bmask); - - 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/common.h.in b/libnnpdf/src/NNPDF/common.h.in deleted file mode 100644 index 74e958216f..0000000000 --- a/libnnpdf/src/NNPDF/common.h.in +++ /dev/null @@ -1,35 +0,0 @@ -#pragma once - -// $Id: common.h 2701 2015-03-18 14:54:20Z s1044006 $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -@LIBNNPDF_HAVE_SSE@ -@LIBNNPDF_HAVE_MPI@ - -/* Universal includes */ -#include -#include - -namespace NNPDF -{ - /*! Returns the version string for libnnpdf */ - std::string getVersion(); - - // ***** Convolution alignment targets **** - typedef float real; - #ifdef SSE_CONV - const int convoluteAlign = 16/sizeof(real); - #else - const int convoluteAlign = 1; - #endif - //Make extern to avoid multiple declaration - void SetVerbosity(int); - int GetVerbosity(); - std::ostream& get_logger(); -} - diff --git a/libnnpdf/src/NNPDF/commondata.h b/libnnpdf/src/NNPDF/commondata.h deleted file mode 100644 index 293f165e05..0000000000 --- a/libnnpdf/src/NNPDF/commondata.h +++ /dev/null @@ -1,252 +0,0 @@ -// $Id: commondata.h 3004 2015-06-11 16:13:23Z s1044006 $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once - -#include "common.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -/** @defgroup commondata CommonData - * \brief Classes and functions related to the handling of NNPDF standard CommonData files. - */ - -/*! \addtogroup commondata - * @{ - */ - -namespace NNPDF -{ - /*! - * \brief Consted information on a CommonData file - precursor to full CommonData class. - * - * This struct is used to communicate the basic required information to build a CommonData instance. - * In this way, the general metadata is read from the file outside the CommonData constructor, - * and therefore can be consted in the CommonData instance. This should therefore be used when - * reading CommonData FILES. - */ - struct dataInfo - { - /*! - * dataInfo constructor (required by SWIG) - */ - dataInfo(int nD, int nS, std::string sN, std::string tF, std::string sF): - nData(nD), - nSys(nS), - SetName(sN), - targetFile(tF), - systypeFile(sF) - {}; - - const int nData; //!< Number of datapoints in prototype CommonData - const int nSys; //!< Number of systematic uncertainties in prototype CommonData - - const std::string SetName; //!< Prototype CommonData set name - - const std::string targetFile; //!< Target filename for prototype CommonData - const std::string systypeFile; //!< Target filename for prototype SYSTYPE - }; - - /*! - * \brief Minimal Consted information on CommonData - precursor to full CommonData class. - * - * This struct is used to communicate the basic required information to build a CommonData instance, - * in the event where an input file is not used (for example, in buildmaster when constructing an empty - * commondata) - */ - struct dataInfoRaw - { - /*! - * dataInfoRaw constructor (required by SWIG) - */ - dataInfoRaw(int nD, int nS, std::string sN, std::string pT): - nData(nD), - nSys(nS), - SetName(sN), - ProcType(pT) - {}; - - const int nData; //!< Number of datapoints in prototype CommonData - const int nSys; //!< Number of systematic uncertainties in prototype CommonData - const std::string SetName; //!< Prototype CommonData set name - const std::string ProcType; //!< Process type for prototype CommonData - }; - - /*! Specifies which type of systematic error */ - enum sysType { - ADD, //!< Additive systematic - MULT, //!< Multiplicative systematic - UNSET //!< Unset systematic - }; - - /*! - * \brief Systematics struct - contains information on a single source of systematic error - * - * This struct is used to detail how individual sources of systematic error should be treated, - * specifying whether they are additive or multiplicative, or whether their treatment should be - * randomised. Furthermore the actual value of the uncertainty in the two cases is also held. - */ - struct sysError - { - double add; //!< Value of the systematic error when additive - double mult; //!< Value of the systematic error when multiplicative - sysType type; //!< Type of the systematic error (ADD/MULT) - std::string name; //!< Name of the systematic error (for correlation purposes) - bool isRAND; //!< Tag specifying whether the treatment should be randomised - - // Normal constructor - sets quiet_NaNs as sentinel values - sysError(): - add(std::numeric_limits::quiet_NaN()), - mult(std::numeric_limits::quiet_NaN()), - type(UNSET), - name("CORR"), - isRAND(false) - {}; - - // Copy constructor - sysError(sysError const& o): - add(o.add), - mult(o.mult), - type(o.type), - name(o.name), - isRAND(o.isRAND) - {}; - }; - - /*! - * \brief Generate a dataInfo struct from file - * - * This function returns a new dataInfo struct from the provided filenames. - * The return value is used to generate a full CommonData instance. - */ - dataInfo genInfoStruct(std::string const& targetfile, std::string const& sysfile); - - /*! - * \brief Extract systematics suffix number - * - * This function returns the systematic ID present in the filename provided. - */ - std::string extractSysID(std::string const& sysfile); - - /** - * \class CommonData - * \brief Class to handle NNPDF CommonData file format - */ - class CommonData { - - private: - CommonData(); //!< Disable default constructor - - static void VerifyProc(std::string const& proc); //!< Verify process types - - protected: - CommonData(dataInfo const&); //!< Constructor from file - CommonData(dataInfoRaw const&); //!< Constructor for empty CommonData - - std::string fSetName; //!< Dataset name - int fNData; //!< Number of datapoints - double *fData; //!< Data value array - - // Kinematical variables - std::string *fProc; //!< Process (determines kinematics) - double *fKin1; //!< First kinematic variable array (x/pT) - double *fKin2; //!< Second kinematic variable array (Q/y) - double *fKin3; //!< Third kinematic variable array (inelasticity) - - // Uncertainties - int fNSys; //!< Number of systematic errors - std::string fSysId; //!< SYSTYPE ID - double *fStat; //!< Statistical error array - sysError **fSys; //!< Systematic errors - - /*! - * Read data from provided CommonData filenames - */ - void ReadData(std::string const& targetfile, std::string const& sysfile); - - public: - // ******************************* CommonData Public Constructors ************************************* - - CommonData(const CommonData& set); //!< copy constructor - CommonData(CommonData&& set); //!< Move constructor - CommonData& operator=(CommonData); //!< Copy-assignment - CommonData& operator=(CommonData&&); //!< Move-assignment - CommonData(const CommonData& set, std::vector const& mask); //!< Masked copy constructor - friend void swap(CommonData&, CommonData& ); - - virtual ~CommonData(); //!< The destructor. - - // ******************************* Process types ******************************************** - - typedef std::map > kinMap; - static const kinMap kinLabel_latex; - static const kinMap kinLabel; - // ******************************* CommonData Get Methods ************************************* - - std::string const& GetSetName() const {return fSetName; }; //!< Returns set name - - int const& GetNData() const { return fNData;}; //!< Returns number of data points - int const& GetNSys() const { return fNSys;}; //!< Returns number of systematic uncertainties - - double const& GetData(int i) const { return fData[i];} //!< Return data value for point i - double const& GetStat(int i) const { return fStat[i];} //!< Return statistical uncertanty for point i - double GetUncE(int i) const; //!< Return total uncorrelated error for point i - double GetCorE(int i) const; //!< Return total correlated error for point i - - sysError const& GetSys(int i, int l) const { return fSys[i][l];} //!< Return lth systematic for point i - sysError** GetSysErrors() const { return fSys; } //!< Return full systematic matrix - - double* GetData() const { return fData; }; //!< Return whole data array - - // Process/Kinematics - std::string const& GetProc(int i) const { return fProc[i]; }; //!< Return the process type for point i - double const& GetKinematics(int const& idat, int const& ikin) const //!< Return the ikinth kinematic value for point idat - { - if (idat < 0 || idat >= fNData) - { - throw std::out_of_range("CommonData::GetKinematics: Data point out of range"); - } - - switch (ikin) - { - case 0: - return fKin1[idat]; - break; - - case 1: - return fKin2[idat]; - break; - - case 2: - return fKin3[idat]; - break; - - default: - throw std::out_of_range("CommonData::GetKinematics: Kinematical variable out of range"); - } - }; - - // ********************************* CommonData File IO ***************************************** - - static CommonData ReadFile(std::string const& filename, std::string const& sysfile); //!< Returns a new CommonData read from file - void Verify() const; //!< Verifies the current CommonData, checking that all fields are set - void Export(std::string const& targetdir) const; //!< Writes the current CommonData instance to file - }; - -} - /*! @} */ diff --git a/libnnpdf/src/NNPDF/config.h.in b/libnnpdf/src/NNPDF/config.h.in deleted file mode 100644 index 07f5548113..0000000000 --- a/libnnpdf/src/NNPDF/config.h.in +++ /dev/null @@ -1 +0,0 @@ -#define VERSION @VERSION@ diff --git a/libnnpdf/src/NNPDF/dataset.h b/libnnpdf/src/NNPDF/dataset.h deleted file mode 100644 index f28eee54f2..0000000000 --- a/libnnpdf/src/NNPDF/dataset.h +++ /dev/null @@ -1,85 +0,0 @@ -// $Id: dataset.h 2990 2015-06-08 19:06:21Z s1044006 $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once - -#include "common.h" -#include -#include -#include - -#include "fkset.h" -#include "pdfset.h" -#include "commondata.h" -#include "utils.h" - -namespace NNPDF -{ - class ThPredictions; - - /** - * \class DataSet - * \brief Class for datasets loading and computing of observables - */ - class DataSet : public CommonData, public FKSet - { - private: - bool fIsArtificial; //!< Flag to determine if data is artifical - bool fIsT0; //!< Flag to determine if covmat is T0 - - // Data information - std::vector fT0Pred; //!< The t0 predictions - defaults to data in case of no t0-std::vector - mutable matrix fCovMat; //!< The covariance matrix - mutable matrix fSqrtCov; //!< The Cholesky decomposition of the covariance matrix - - // private methods for constructor - void GenCovMat() const; //!< Generate covariance matrix - double fWeight = 1; //!< Factor that divides the covariance matrix and increases - //!the weight of this dataset in the fit proportionally. - - DataSet(); //!< Disable default constructor - - public: - DataSet(CommonData const&, FKSet const&, double weight=1.); //!< Constructor - DataSet(const DataSet&, std::vector const&); //!< Masked Copy constructor - DataSet(const DataSet&) = default; //!< Masked Copy constructor - virtual ~DataSet(); //!< The destructor. - - // **************** DataSet T0 Methods ********************** - - void SetT0( ThPredictions const&); //!< Set T0 predictions - void SetT0(PDFSet const&); //!< Set T0 predictions - bool const& IsT0 () const { return fIsT0; } //!< Return t0 covmat flag - - // ************************ Data Get Methods ****************************** - - double const& GetT0Pred(int i) const { return fT0Pred[i];} //!< Return t0 prediction - double GetWeight() const {return fWeight;} - - matrix const& GetCovMat() const; //!< Return fCovMat - matrix const& GetSqrtCov() const; //!< Return the Cholesky decomposition of the covariance matrix - double GetSqrtCov(int i, int j) const; //!< Returns an element of the Cholesky decomposition - - bool const& IsArtificial() const { return fIsArtificial; } //!< Returns the artificial flag - - // **************** Update data values ******************************** - - void UpdateData(double* newdat); //!< Update data - void UpdateData(double* newdat, double* norm); //!< Update with a rescaling - also rescales additive uncertainties - - void SetArtificial( bool const& artificial) { fIsArtificial = artificial; } - - // **************** Rescale Uncertainties **************************** - - void RescaleErrors(const double mult); //!< Rescale uncertainties by a multiplicative factor - - // ************************************************************************** - - }; - -} diff --git a/libnnpdf/src/NNPDF/exceptions.h b/libnnpdf/src/NNPDF/exceptions.h deleted file mode 100644 index 6dbfd805ae..0000000000 --- a/libnnpdf/src/NNPDF/exceptions.h +++ /dev/null @@ -1,80 +0,0 @@ -// $Id$ -// -// NNPDF++ 2015 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it - -#pragma once - -#include -#include -#include - -namespace NNPDF -{ - /** - * - * Simple and generic exception handlers - * - implements runtime and logic errors - * - */ - - /// Error to be thrown when runtime error is detected - class RuntimeException: public std::runtime_error - { - public: - RuntimeException(const std::string& tag, const std::string& what) : std::runtime_error(std::string("[") + std::string(tag) + std::string("] error: ") + std::string(what)) {} - }; - - /// Error to be thrown when logic error is detected - class LogicException: public std::logic_error - { - public: - LogicException(const std::string& tag, const std::string& what) : std::logic_error(std::string("[") + std::string(tag) + std::string("] error: ") + std::string(what)) {} - }; - - //_______________________________________________________________________ - class FileError: public RuntimeException - { - public: - FileError(const std::string& tag, const std::string& what) : RuntimeException(tag,what) {} - }; - - class EvaluationError: public RuntimeException - { - public: - EvaluationError(const std::string& tag, const std::string& what) : RuntimeException(tag,what) {} - }; - - class InitError: public RuntimeException - { - public: - InitError(const std::string& tag, const std::string& what) : RuntimeException(tag,what) {} - }; - - class RangeError: public LogicException - { - public: - RangeError(const std::string& tag, const std::string& what) : LogicException(tag,what) {} - }; - - class LengthError: public LogicException - { - public: - LengthError(const std::string& tag, const std::string& what) : LogicException(tag,what) {} - }; - - class LogError: public LogicException - { - public: - LogError(const std::string& tag, const std::string& what) : LogicException(tag,what) {} - }; - - class UserError: public LogicException - { - public: - UserError(const std::string& tag, const std::string& what) : LogicException(tag,what) {} - }; - -} diff --git a/libnnpdf/src/NNPDF/experiments.h b/libnnpdf/src/NNPDF/experiments.h deleted file mode 100644 index 0eba7826a8..0000000000 --- a/libnnpdf/src/NNPDF/experiments.h +++ /dev/null @@ -1,112 +0,0 @@ -// $Id: experiments.h 2990 2015-06-08 19:06:21Z s1044006 $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once - -#include -#include - -#include "commondata.h" -#include "dataset.h" -#include "pdfset.h" -#include "utils.h" - -namespace NNPDF -{ - /** - * \class Experiment - * \brief Class for handling experiments - */ - class Experiment - { - public: - Experiment(std::vector const& sets, std::string const& expname); //!< Experiment constructor - Experiment(Experiment const&, std::vector const& sets); //!< Copy constructor with dataset replace - Experiment(Experiment const&); //!< Copy constructor - ~Experiment(); //!< Experiment destructor - - void MakeReplica(); //!< Shift the exp data and produces art. data - void MakePerPointReplica(int pointindex); - void MakeClosure(const std::vector& predictions, bool const& noise); - void MakeClosure(PDFSet* pdf, bool const& noise); //!< Make fake exp data using theory predictions - - int GetNSet() const {return (int) fSets.size();} //!< Return the number of sets in the experiment - DataSet const& GetSet(int i) const {return fSets[i];} //!< Return the ith DataSet - std::vector const& DataSets() const {return fSets;} //!< Return a reference to the vector of DataSets. - - std::string const& GetExpName() const {return fExpName;} //!< Return the experiment name - std::string const& GetSetName(int i) const { return fSets[i].GetSetName(); } //!< Return the dataset name - - int GetNData() const { return fNData; } //!< Return the number of points in the experiment - const double* GetData() const { return fData.data(); } - - bool IsArtificial() const { return fIsArtificial; } //!< Return the artificial flag - bool IsClosure() const { return fIsClosure; } //!< Return the artificial flag - bool IsT0() const { return fIsT0; } //!< Return t0 covmat flag - - matrix const& GetCovMat() const { return fCovMat; } //!< Return fCovMat - matrix const& GetSqrtCov() const { return fSqrtCov; } //!< Return the Cholesky decomposition of the covariance matrix - - void ExportCovMat(std::string); //!< Export covariance matrix - void ExportSqrtCov(std::string); //!< Export Cholesky decomposition - - void SetT0(const PDFSet&); //! bmask = {}); //!< Read in covmat for rep gen, generate covmat and sqrt - void LoadFitCovMat(std::string filename, bool ThUnc, std::vector bmask = {}); //!< Read in covmat for fitting, generate covmat and sqrt - matrix GetSqrtFitCovMat(std::string filename, std::vector bmask = {}); //temporary external function for calc sqrtcovmat with th error - - private: - - Experiment(); //disable default constructor - Experiment& operator=(const Experiment&) ; //disable copy-assignment - - std::string fExpName; //!< Stores the name of the experiment - - std::vector fSets; //!< Vector of contained datasets - - int fNData; //!< Number of data points - int fNSys; //!< Number of additive systematics correlations - - std::vector fData; //!< The experimental data - std::vector fT0Pred; //!< The t0 predictions - std::vector fSqrtWeights; //!< The weights - - matrix fCovMat; //!< The covariance matrix - matrix fSqrtCov; //!< The Cholesky decomposition of the covariance matrix - matrix fSamplingMatrix; //!< The Cholesky decomposition of the total replica generation covmat - - std::vector fStat; //!< The statistical errors - sysError **fSys; //!< The syscor - int **fSetSysMap; //!< Map for ordering of systematics in datasets - - bool fIsArtificial; //!< The artificial flag - bool fIsClosure; //!< If data is closure data - bool fIsT0; //!< Flag to determine if covmat is T0 - - - void ClearLocalData(); //!< Clear data pulled from datasets - - void PullData(); //!< Pull experimental data from datasets - void GenCovMat(); //!< Generate covmat and inverse - - }; - - // Auxiliary tools for replica generation - - /** - * @brief pseudodata Returns the artificial data replica - * used by nnfit for a given data seed and replica number. - * - * Note that the user is responsable to destroy the objects. - */ - std::vector pseudodata(std::vector const& exps, - unsigned long int dataseed, - int replica_id); - -} diff --git a/libnnpdf/src/NNPDF/fastkernel.h b/libnnpdf/src/NNPDF/fastkernel.h deleted file mode 100644 index 45ae4b2a40..0000000000 --- a/libnnpdf/src/NNPDF/fastkernel.h +++ /dev/null @@ -1,202 +0,0 @@ -// $Id$ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -// libnnpdf 2014 nph - -#pragma once - -#include -#include -#include -#include -#include - -#include - -#include "common.h" - -namespace NNPDF -{ - - // Helper for gcc - template - std::string ToString(T val) - { - std::stringstream stream; - stream << val; - return stream.str(); - } - - /** - * \class FKHeader - * \brief Class for parsing the header of FK tables - */ - class FKHeader - { - public: - FKHeader(); //!< Prototype constructor - FKHeader(std::istream&); //!< Construct from istream - FKHeader(std::string const& filename); //!< Construct from file - FKHeader(FKHeader const&); //!< Copy-construct - ~FKHeader(); //!< Destructor - - void Read(std::istream&); //!< Read FKTable header from ostream - void Print(std::ostream&) const; //!< Print FKTable header to ostream - void ResetFlavourMap(); //!< Resets the flavourmap to the maximal version - - typedef std::map keyMap; - typedef enum {VERSIONS, GRIDINFO, THEORYINFO, BLOB} section; - - // ********************************* Tag Manip ********************************* - - template - void AddTag( section sec, std::string const& key, T const& value) - { AddTag(sec, key,ToString(value)); } - void AddTag( section sec, std::string const& key, const char* value) - { AddTag(sec, key,ToString(value));} - void AddTag( section sec, std::string const& key, std::string const& value); - - // ********************************* Tag Getters ********************************* - - bool HasTag( section sec, std::string const& key ) const; - std::string GetTag( section sec, std::string const& key) const; - template T GetTag( section sec, std::string const& key) const - { return static_cast(atof(GetTag(sec,key).c_str())); } - - protected: - void RemTag( section sec, std::string const& key ); //!< Remove existing tags - - // Printing helper functions - std::string SectionHeader(const char* title, section) const; - - void PrintKeyValue(std::ostream&, section) const; - void ReadKeyValue( std::istream& is, section sec ); - - void PrintBlob(std::ostream& os, std::string blobname) const; - void ReadBlob(std::istream& is, std::string blobname); - - // Map fetchers helper functions - section GetSection(std::string const&) const; //!< Fetch the appropriate section for a title - const keyMap* GetMap(section const& sec) const; //!< Fetch the appropriate map for a section - keyMap* GetMap(section const& sec) //!< Fetch the appropriate map for a section - { return const_cast(const_cast(this)->GetMap(sec)); }; - - // ********************************* Attributes ********************************* - - // Key-value pairs - keyMap fVersions; - keyMap fGridInfo; - keyMap fTheoryInfo; - keyMap fBlobString; - - friend class FKTable; - }; - - /** - * \class FKTable - * \brief Class for holding FastKernel tables - */ - class FKTable - { - public: - FKTable( std::istream& , - std::vector const& cFactors = std::vector() - ); // Stream constructor - FKTable( std::string const& filename, - std::vector const& cfactors = std::vector() - ); // Filename constructor - FKTable(FKTable const&); //!< Copy constructor - FKTable(FKTable const&, std::vector const&); //!< Masked copy constructor - - virtual ~FKTable(); //!< Destructor - void Print(std::ostream&); //!< Print FKTable header to ostream - void Print(std::string const&, bool const& compress = true); - - // ******************** FK Get Methods *************************** - - std::string const& GetDataName() const {return fDataName;} - - double const& GetQ20() const { return fQ20; } - double const* GetCFactors() const { return fcFactors; } - double const* GetCUncerts() const { return fcUncerts; } - - int const& GetNData() const { return fNData;} //!< Return fNData - int const& GetNx() const { return fNx; } //!< Return fNx - int const& GetTx() const { return fTx; } //!< Return fTx - int const& GetDSz() const { return fDSz; } //!< Return fDSz - - double * GetXGrid() const { return fXgrid; } //!< Return fXGrid - real * GetSigma() const { return fSigma; } //!< Return fSigma - - int * GetFlmap() const { return fFlmap; } //!< Return fFlmap - int const& GetNonZero() const { return fNonZero; } //!< Return fNonZero - - bool const& IsHadronic() const { return fHadronic;} //!< Return fHadronic - - std::string GetTag(FKHeader::section sec, std::string const& key) const { return fFKHeader.GetTag(sec, key); } - - // GetISig returns a position in the FK table - int GetISig( int const& d, // Datapoint index - int const& ix1, // First x-index - int const& ix2, // Second x-index - int const& ifl1, // First flavour index - int const& ifl2 // Second flavour index - ) const; - - // DIS version of GetISig - int GetISig( int const& d, // Datapoint index - int const& ix, // x-index - int const& ifl // flavour index - ) const; - - protected: - void ReadCFactors(std::string const& filename); //!< Read C-factors from file - bool OptimalFlavourmap(std::string& flmap) const; //!< Determine and return the optimal flavour map - - - // FKHeader - FKHeader fFKHeader; - - // Metadata - std::string fDataName; - std::string fDescription; - int fNData; - - // Process information - double fQ20; - bool fHadronic; - int fNonZero; - int * fFlmap; - - // x-grid information - int fNx; - int fTx; - int fRmr; - int fDSz; - - // X-arrays - double * fXgrid; - - // FK table - real * fSigma; - - // Cfactor information and uncertainties - bool fHasCFactors; - double * fcFactors; - double * fcUncerts; // the *squared* uncorrelated uncertainty of the C-factors - - private: - FKTable(); //!< Disable default constructor - FKTable& operator=(const FKTable&); //!< Disable copy-assignment - - void InitialiseFromStream(std::istream&, std::vector const& cFactors); //!< Initialise the FK table from an input stream - - int parseNonZero(); // Parse flavourmap information into fNonZero - }; - -} diff --git a/libnnpdf/src/NNPDF/fkgenerator.h b/libnnpdf/src/NNPDF/fkgenerator.h deleted file mode 100644 index 6e243a4bce..0000000000 --- a/libnnpdf/src/NNPDF/fkgenerator.h +++ /dev/null @@ -1,57 +0,0 @@ -// $Id$ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -// libnnpdf 2014 nph - -#pragma once - -#include -#include -#include -#include "NNPDF/common.h" -#include "NNPDF/fastkernel.h" - -namespace NNPDF -{ - /** - * \class FKGenerator - * \brief Class for filling FastKernel tables - */ - class FKGenerator : public FKTable - { - private: - FKGenerator(); //!< Disable default constructor - FKGenerator( const FKGenerator& ); //!< Disable default copy-constructor - FKGenerator& operator=(const FKGenerator&); //!< Disable copy-assignment - - std::unique_ptr fAccumulator; - - public: - FKGenerator( std::istream& ); // Constructor - virtual ~FKGenerator() = default; - - void Finalise(); //!< Copy Accumulator to fSigma - - // Fill the FKTable - void Fill( int const& d, // Datapoint index - int const& ix1, // First x-index - int const& ix2, // Second x-index - int const& ifl1, // First flavour index - int const& ifl2, // Second flavour index - real const& isig // FK Value - ); - - // DIS version of Fill - void Fill( int const& d, // Datapoint index - int const& ix, // x-index - int const& ifl, // flavour index - real const& isig // FK Value - ); - }; - -} diff --git a/libnnpdf/src/NNPDF/fkset.h b/libnnpdf/src/NNPDF/fkset.h deleted file mode 100644 index ff08bbc8b5..0000000000 --- a/libnnpdf/src/NNPDF/fkset.h +++ /dev/null @@ -1,67 +0,0 @@ -// $Id$ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -// libnnpdf 2014 nph - -#pragma once - -#include -#include - -#include "common.h" -#include "fastkernel.h" - -namespace NNPDF -{ - class FKTable; - - // Operator for compound FK tables - typedef void (*SigmaOp) (int const&, std::vector const&, real* ); - - /** - * \class FKSet - * \brief Class for holding sets of FK tables, alongside information on compound observable computation - */ - class FKSet - { - public: - FKSet(SigmaOp, std::vector const&); //!< Constructor - virtual ~FKSet(); //!< Destructor - - FKSet(FKSet const&); //!< Copy-constructor - FKSet(FKSet &&); //Move constructor - FKSet(FKSet const&, std::vector const&); //!< Masked copy-constructor - - SigmaOp GetOperator() const { return *fOperator;}; - - int const& GetNSigma() const { return fNSigma; }; - int const& GetNDataFK() const { return fNDataFK; }; - bool const& IsHadronic() const { return fHadronic; }; - - std::string const& GetDataName() const {return fDataName;}; - const FKTable* GetFK(size_t const& i) const {return fFK[i];}; //!< Fetch FK tables - - // parse operator string - static SigmaOp parseOperator(std::string const& op); - - friend void swap(FKSet &, FKSet &); - FKSet& operator=(FKSet); //copy-assignment - - private: - FKSet(); //disable default constructor - - SigmaOp fOperator; //!< Operator acting on combination observables - int fNSigma; //!< Number of FastKernel Grids for this dataset - int fNDataFK; //!< Number of datapoints - bool fHadronic; //!< Hadronic Observables - - std::string fDataName; //!< Name of dataset - FKTable **fFK; //!< FastKernel tables - }; - -} diff --git a/libnnpdf/src/NNPDF/lhapdfset.h b/libnnpdf/src/NNPDF/lhapdfset.h deleted file mode 100644 index a4ca343bd2..0000000000 --- a/libnnpdf/src/NNPDF/lhapdfset.h +++ /dev/null @@ -1,47 +0,0 @@ -// $Id: lhapdfset.h 3177 2015-08-18 14:43:31Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once -#include "pdfset.h" -#include "LHAPDF/LHAPDF.h" - -#include - -namespace NNPDF -{ - - /** - * \class LHAPDFSet - * \brief Class for handle LHAPDF methods - */ - class LHAPDFSet : private LHAPDF::PDFSet, public NNPDF::PDFSet - { - public: - //!< Constructor - LHAPDFSet(std::string const&, erType); - LHAPDFSet(std::string const&, int const&); - LHAPDFSet(const LHAPDFSet&) = delete; - LHAPDFSet& operator=(const LHAPDFSet&) = delete; - - virtual ~LHAPDFSet(); //!< Destructor - - //!< Get PDF array in evolution basis - virtual void GetPDF(real const& x, real const& Q2, int const& n, real* pdf) const; - - //!< Check if flavor exists - bool hasFlavor(int pdgid) const; - - //!< Return single flavor - real xfxQ(real const& x, real const& Q, int const& n, int const& fl) const; - - protected: - // Member PDFs - std::vector fMemberPDFs; - - }; -} diff --git a/libnnpdf/src/NNPDF/logger.h b/libnnpdf/src/NNPDF/logger.h deleted file mode 100644 index b9460c8d7b..0000000000 --- a/libnnpdf/src/NNPDF/logger.h +++ /dev/null @@ -1,79 +0,0 @@ -// $Id -// -// NNPDF++ 2013 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once - -#include "common.h" - -#include -#include -#include -#include - -namespace NNPDF -{ - - class Logger; - - /** - * \class LogManager - * \brief Container class for individual loggers - */ - class LogManager - { - public: - static void InitPath(std::string const& path) {GetLM()->fBasePath = path;}; //!< Initialise the Manager's base path - - static void AddLogger(std::string const&, std::string const&); //!< Add a new log file to the manager - static Logger& GetLogger(std::string const&); //!< Get a log from the manager - - static void AddLogEntry(std::string const& log, std::string ent); //!< Add entry 'ent' to log ID 'log' - - static void ExportLogs(); //!< Export all active loggers - - private: - LogManager(); - ~LogManager(); - - // LogManager get/init method - static LogManager* GetLM(){ - if (logInstance == 0) - logInstance = new LogManager(); - - return logInstance; - }; - - // Attributes - static LogManager* logInstance; //!< Singleton log instance - - typedef std::map LogMap; //!< LogMap typedef - LogMap fMap; //!< hash map of available loggers - std::string fBasePath; //!< Base path of logger - - }; - - /** - * \class Logger - * \brief Basic logging class. All methods private - can only be handled by logmanager - */ - class Logger - { - Logger(std::string const&, std::string const&); - - void AddLogEntry(std::string log) {fLogEntries.push_back(log);}; - void Export(); - - const std::string fLogname; - const std::string fFilename; - - std::vector fLogEntries; - - friend class LogManager; - }; - -} diff --git a/libnnpdf/src/NNPDF/nnmpi.h b/libnnpdf/src/NNPDF/nnmpi.h deleted file mode 100644 index 4ede9b3fe1..0000000000 --- a/libnnpdf/src/NNPDF/nnmpi.h +++ /dev/null @@ -1,49 +0,0 @@ -// $Id$ -// -// NNPDF++ 2015 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it - -#pragma once - -#include "common.h" - -namespace NNPDF -{ - /** - * @brief The MPI class: handler for openMPI world - */ - class MPI - { - public: - static void Init(); - static void Finalize(); - static int TaskID() { return getInstance()->ftaskid; } - static int NumTasks() { return getInstance()->fnumbertasks; } - static void DecomposeDomain(int domain_size, int taskid, - int* subdomain_start, - int* subdomain_size); - static void SendJobs(int, int, int, real*, real*, real*); - static void RecvJobs(int, int, real *); - static void ResetWorld(); - - private: - MPI(); - ~MPI(); - void SlaveConvoluteLoop(); - - static MPI* getInstance() - { - if (!mpiInstance) - mpiInstance = new MPI(); - return mpiInstance; - } - - // Members - int fnumbertasks; - int ftaskid; - int fworld; - static MPI* mpiInstance; - }; -} diff --git a/libnnpdf/src/NNPDF/nnpdfdb.h b/libnnpdf/src/NNPDF/nnpdfdb.h deleted file mode 100644 index 4ed9581283..0000000000 --- a/libnnpdf/src/NNPDF/nnpdfdb.h +++ /dev/null @@ -1,92 +0,0 @@ -#pragma once -/* - * nnpdfdb.h - * Database access routines for nnpdf family applications - * * nathan.hartland@physics.ox.ac.uk 02/15 - */ - -#include -#include -#include -#include -#include -#include - -#include "exceptions.h" - -namespace NNPDF -{ - - class IndexDB; - template - T dbquery(IndexDB const& db, int const& id, std::string const& field); - - // Index-based databases - class IndexDB - { - public: - IndexDB(std::string const& _path, std::string const& _table); - ~IndexDB(); - - const int& GetNEntries() const {return fNEntries;}; - void ExtractMap( const int& id, std::vector const& keys, std::map& map); - std::vector ExtractString(const int &id, const std::vector &keys); - - private: - sqlite3* fDB; - - const std::string fPath; - const std::string fTable; - - int fNEntries; - - // Returns the value of a key for a database entry ID - template - friend T dbquery(IndexDB const& db, int const& id, std::string const& field); - - // Returns a list of database entry IDs with field matching value - friend std::vector dbmatch(IndexDB const& db, std::string const& field, std::string const& value); - }; - - // *************************************************************************************************************** - - // Return a queried value - template - T dbquery(IndexDB const& db, int const& id, std::string const& field) - { - std::stringstream query; - query << "select " << field << " from "<> retval; - - return retval; - } else - { - std::stringstream err; err << "SQLITE Error: " << sqlite3_errmsg(db.fDB); - throw RuntimeException("dbquery", err.str()); - } - } - - // Needs to be specialised to avoid splitting the string - template<> - std::string dbquery(IndexDB const& db, int const& id, std::string const& field); - std::vector dbmatch(IndexDB const& db, std::string const& field, std::string const& value); - -} diff --git a/libnnpdf/src/NNPDF/parametrisation.h b/libnnpdf/src/NNPDF/parametrisation.h deleted file mode 100644 index 594e7cffd0..0000000000 --- a/libnnpdf/src/NNPDF/parametrisation.h +++ /dev/null @@ -1,128 +0,0 @@ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once - -#include -#include -#include -#include -using actfunction = std::function; - -#include "common.h" - -namespace NNPDF -{ - /** - * \class Parametrisation - * \brief Abstract parametrisation base class - */ - class Parametrisation - { - public: - Parametrisation(std::string name, int nParameters); - Parametrisation(Parametrisation const& ); - virtual ~Parametrisation(); - - virtual void Compute(real*,real*) const = 0; - virtual Parametrisation* Duplicate() = 0; - void CopyPars(Parametrisation const* ); - - void ReadPars(std::ifstream&); - void WritePars(std::ofstream&); - void SetPars(std::vector const& param); - - real* GetParameters() {return fParameters;} - int const& GetNParameters() const {return fNParameters;} - - std::string const& GetParamName() const {return fParamName;}; - - virtual void InitParameters() = 0; - - protected: - const int fNParameters; //!< Total number of parameters - real* const fParameters; //!< Parameters - - private: - const std::string fParamName; - }; - - /** - * \class MultiLayerPerceptron - * \brief General MLP class - */ - class MultiLayerPerceptron : public Parametrisation - { - public: - MultiLayerPerceptron(std::vector const& arch); //!< Network constructor - MultiLayerPerceptron(MultiLayerPerceptron const&); //!< Network copy constructor - ~MultiLayerPerceptron(); //!< Network destructor - - void InitParameters(); //!< Initialize (or reinitialize) parameters - void Compute(real*,real*) const; //!< Returns a fArch[fNLayers-1] long array of output for a given input array - Parametrisation* Duplicate(); //!< Returns a parametrisation based on an MLP - - void SetActivationFunction(actfunction const& f) { fActFunction = f; } //!< set the activation function - - // Nodal GA - const int* GetArch() const {return fArch;} - //!< Returns the number of parameters that belong to a specific node (including biases). - int GetNumNodeParams(int const& layer) const; - //!< Returns a pointer to the fParameters coordinate representing the parameters for a specific node - real* GetNodeParams (int const& layer, int const& node); - - protected: - const int fNLayers; //!< Number of layers - int* fArch; //!< Network architecture - - real** fWeightMatrix; //!< Weights/Thresholds - mutable real** fOutputMatrix; //!< Neuron Activation - - actfunction fActFunction; //!< Activation function - }; - - /*! - * \class SingleLayerPerceptron - * \brief An implementation of a perceptron with a single hidden layer - * This class can take an arbitrary number of extra parameters, for use - * in derived classes which are, for example, fitting preprocessing exponents. - */ - class SingleLayerPerceptron : public Parametrisation - { - public: - SingleLayerPerceptron(std::vector const& arch, unsigned int extra_pars = 0); - virtual Parametrisation *Duplicate(); - virtual void Compute(real*,real*) const; - virtual void InitParameters(); - protected: - // Number of hidden nodes - const int fNHidden; - }; - - /*! - * \class SingleLayerPerceptronPreproc - * \brief A Single-hidden-layer perceptron with preprocessing - * Preprocessing is performed to the output of the network as - * output*x^{|alpha|}(1-x)^{|beta|} - * where alpha and beta are the last and second-to-last parameters - * respectively. - */ - class SingleLayerPerceptronPreproc : public SingleLayerPerceptron - { - public: - SingleLayerPerceptronPreproc(std::vector const& arch): - SingleLayerPerceptron(arch, 2){}; - virtual Parametrisation *Duplicate(); - virtual void InitParameters(); - void Compute(real*,real*) const; - protected: - // This attribute scales the parameters corresponding to preprocessing - // exponents in order to balance the training between them and NN params. - const real fScaleFac = 0.2; - }; - -} diff --git a/libnnpdf/src/NNPDF/pathlib.h b/libnnpdf/src/NNPDF/pathlib.h deleted file mode 100644 index ff7ae7bc38..0000000000 --- a/libnnpdf/src/NNPDF/pathlib.h +++ /dev/null @@ -1,42 +0,0 @@ -#pragma once -#include - -/** - * pthlib.h - * - * Locate the nnprofile configuration file and extract basic - * locations needed for other C++ projecs. - */ - -namespace NNPDF -{ -/** - * @brief get_profile_path. - * This is detemined by the environment variable - * NNPDF_PROFILE_PATH if present, and otherwise by the compile time - * constant DEFAULT_NNPDF_PROFILE_PATH. - * @return The location of the NNPDF configuration file. - */ -std::string get_profile_path(); - -/** - * @brief get_data_path - * @return The location of the nnpdfcpp data path, - * containing experimental data and theory settings. - */ -std::string get_data_path(); - -/** - * @brief get_results_path - * @return The location of the nnpdfcpp results path, containing - * completed fits. - */ -std::string get_results_path(); - -/** - * @brief get_config_path - * @return The location of the nnpdfcpp config path, containing - * default configurations for plotting and fiatlux - */ -std::string get_config_path(); -} // namespace NNPDF diff --git a/libnnpdf/src/NNPDF/pdfset.h b/libnnpdf/src/NNPDF/pdfset.h deleted file mode 100644 index 9115352f31..0000000000 --- a/libnnpdf/src/NNPDF/pdfset.h +++ /dev/null @@ -1,105 +0,0 @@ -// $Id: pdfset.h 3177 2015-08-18 14:43:31Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once - -#include -#include -#include "common.h" - -#include - -namespace NNPDF -{ - /** - * \class PDFSet - * \brief An abstract class that provides the basic interface to PDF sets - */ - - class PDFSet - { - public: - //! Error types (Etype) - enum class erType - { ER_NONE, //!< PDF set is not an error set - ER_MC, //!< 1sigma error for NNPDF - ER_MC68, //!< 68cl error for NNPDF - ER_MCT0, //!< NNPDF set use only replica 0 for T0 - ER_EIG, //!< 1sigma error for CTEQ & MSTW - ER_EIG90, //!< 90cl error for CTEQ & MSTW - ER_SYMEIG //!< Symmetric eigenvectors for PDF4LHC - }; - - // LHA-style flavour basis - enum {PHT,TBAR,BBAR,CBAR,SBAR,UBAR,DBAR,GLUON,D,U,S,C,B,T}; - static const std::array fl_labels() { return {{"PHT","TBAR","BBAR","CBAR","SBAR","UBAR","DBAR","GLUON","D","U","S","C","B","T"}}; } - - // NNPDF-style EVLN basis - enum evlnBasis { EVLN_GAM, EVLN_SNG, EVLN_GLU, EVLN_VAL, EVLN_V3, - EVLN_V8, EVLN_V15, EVLN_V24, EVLN_V35, EVLN_T3, - EVLN_T8, EVLN_T15, EVLN_T24, EVLN_T35 }; - - protected: - PDFSet(std::string const&, int const& fMembers, erType const& etype ); //!< Constructor - - const std::string fSetName; //!< Name of the PDF Set - - int fMembers; //!< Number of members in the PDF set - const erType fEtype; //!< Type of error handling for the set - - public: - virtual ~PDFSet(); - - // Verbosity control - static bool Verbose; - - //! Configuration of PDFs - virtual real GetPDF(real const& x, real const& Q2, int const& n, int const& fl) const; //!< Get single PDF - virtual void GetPDF(real const& x, real const& Q2, int const& n, real* pdf) const = 0; //!< Get PDF array - - // Get methods - const std::string& GetSetName() const { return fSetName; } //!< Return the set name - const erType& GetEtype() const { return fEtype; } //!< Fetch the error type of the PDF set - const int& GetMembers() const { return fMembers; } //!< Return number of memeber PDFs - - enum intType { XFX = true, FX = false }; // Integration type helper - // Integrate the PDF - real IntegratePDF( int const& mem, int const& fl, real const& Q2, - PDFSet::intType xfx, bool& gslerror, - gsl_integration_workspace *, - real xmin = 0.0,real xmax = 1.0) const; - - static std::string errString(erType const& type) - { - switch (type) - { - case erType::ER_NONE: - return "No Errors"; - case erType::ER_MC: - return "Monte Carlo"; - case erType::ER_MC68: - return "Monte Carlo 68pc"; - case erType::ER_MCT0: - return "Monte Carlo t0"; - case erType::ER_EIG: - return "Eigenvector 68pc"; - case erType::ER_EIG90: - return "Eigenvector 90pc"; - case erType::ER_SYMEIG: - return "Symmetric eigenvector"; - } - - return "Unrecognised Error Type"; - } - - // Transformations - static void LHA2EVLN(const real *LHA, real *EVLN); - static void EVLN2LHA(const real* EVL, real* LHA); - }; - -} diff --git a/libnnpdf/src/NNPDF/positivity.h b/libnnpdf/src/NNPDF/positivity.h deleted file mode 100644 index bf9569225f..0000000000 --- a/libnnpdf/src/NNPDF/positivity.h +++ /dev/null @@ -1,55 +0,0 @@ -// $Id: positivity.h 3187 2015-08-23 11:08:42Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2013 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk -#pragma once - -#include -#include - -#include "common.h" -#include "fastkernel.h" -#include "commondata.h" - -namespace NNPDF -{ - - class PDFSet; - - /** - * \class PositivitySet - * \brief Positivity observables class manager - */ - - class PositivitySet : public CommonData, public FKTable - { - public: - PositivitySet(CommonData const&, FKTable const&, real const& lambda); //!< Constructor - PositivitySet(const PositivitySet&); //!< Copy constructor - virtual ~PositivitySet(); //!< Destructor. - - void ComputeErf (const PDFSet*, real*) const; //!< Compute error function for supplied PDF set - - void ComputeNViolated( const PDFSet*, int* res) const; //!< Number of violated points - void ComputeNUnacceptable(const PDFSet*, int* res) const; //!< Number of unacceptable points (less than bounds) - - void SetBounds(const PDFSet*); //!< Set bounds - - void GetPredictions(const PDFSet&, real **result, int* ndata, int* npdf); - - private: - - PositivitySet(); // disable default constructor - PositivitySet& operator=(const PositivitySet&); // disable copy-assignment - - void ComputePoints(const PDFSet*, real*) const; //!< Compute the positivity points - - const real fLambda; //!< Lagrange Multiplier - - std::vector fBounds; - }; - -} diff --git a/libnnpdf/src/NNPDF/randomgenerator.h b/libnnpdf/src/NNPDF/randomgenerator.h deleted file mode 100644 index e8cd122986..0000000000 --- a/libnnpdf/src/NNPDF/randomgenerator.h +++ /dev/null @@ -1,54 +0,0 @@ -// $Id: randomgenerator.h 2120 2014-12-02 17:14:52Z s1044006 $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once - -#include "common.h" -#include "gsl/gsl_rng.h" -#include "gsl/gsl_randist.h" - -#include - -namespace NNPDF -{ - /** - * @brief The RandomGenerator class - * This class manages the random number generators - */ - class RandomGenerator - { - private: - RandomGenerator(int,unsigned long int); - ~RandomGenerator(); - - gsl_rng *fR; - - public: - - // Initialise and fetch RNG - static void InitRNG(int method,unsigned long int seed); - static RandomGenerator* GetRNG(); - - // Verbosity control - static bool Verbose; - - unsigned long int GetRandomInt(); - unsigned long int GetRandomUniform(unsigned long int); - double GetRandomUniform(); - double GetRandomUniform(double, double); - double GetRandomUniformPos(); - double GetRandomGausDev(const double); - - template - void ShuffleVector( std::vector& vec ) - {return gsl_ran_shuffle (fR, &vec[0], vec.size(), sizeof(T));} - - // Set Methods - void SetSeed(unsigned long int); - }; -} diff --git a/libnnpdf/src/NNPDF/thpredictions.h b/libnnpdf/src/NNPDF/thpredictions.h deleted file mode 100644 index 20b73b88da..0000000000 --- a/libnnpdf/src/NNPDF/thpredictions.h +++ /dev/null @@ -1,103 +0,0 @@ -// $Id: thpredictions.h 3044 2015-06-18 19:06:34Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once - -#include "NNPDF/common.h" -#include "NNPDF/experiments.h" -#include "NNPDF/pdfset.h" -#include "NNPDF/fastkernel.h" -#include "NNPDF/fkset.h" - -namespace NNPDF -{ - /** - * \class ThPredictions - * \brief Computes the predictions and statistics for a PDFSet and DataSet - */ - class ThPredictions { - - private: - ThPredictions(); //!< Disable default constructor - - real *fObs; //!< Theory predictions for Observables - real fTconv; //!< time required for the convolution in msecs - - int fNpdf; //!< Number of PDFs used in the computation - int fNData; //!< Number of data points - - std::string fPDFName; //!< Name of PDF used to do the calculation - std::string fSetName; //!< Name of dataset in the calculation - - PDFSet::erType fEtype; //!< Uncertainty type - - static void GetNZPDF(const PDFSet*, const FKTable*, real* pdf); //!< Form up PDF array - static void GetNZPDF(const PDFSet*, const PDFSet*, const FKTable*, real* pdf); //!< Form up PDF array (differing beams) - - public: - ThPredictions(const PDFSet*, const Experiment*); //!< Constructor - ThPredictions(const PDFSet*, const FKTable*); //!< Constructor - ThPredictions(const PDFSet*, const FKSet*); //!< Constructor - - 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 - ThPredictions(ThPredictions &&); //!< Move constructor - ~ThPredictions(); //!< Destructor - - // Operators - ThPredictions& operator+=(const ThPredictions&); - ThPredictions& operator-=(const ThPredictions&); - ThPredictions& operator*=(const ThPredictions&); - ThPredictions& operator/=(const ThPredictions&); - - const ThPredictions operator+(const ThPredictions &) const; - const ThPredictions operator-(const ThPredictions &) const; - const ThPredictions operator/(const ThPredictions &) const; - const ThPredictions operator*(const ThPredictions &) const; - - // Fast static convolution functions - 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];}; - real GetObsCV (int const& idat) const; //!< Get Central Value - real GetObsError (int const& idat) const; //!< Get Error value (fEtype) - real GetTConv() const { return fTconv; } //!< Get convolution time - - int GetNPdf() const { return fNpdf; } //!< Return the number of pdfs - int GetNData() const { return fNData; } //!< Return the number of datapoints - - std::string GetSetName() const { return fSetName; } //!< Return the set name - std::string GetPDFName() const { return fPDFName; } //!< Return the set name - - // Set Methods - void SetPDFName( std::string const& newname ) {fPDFName = newname;}; - void SetDataName( std::string const& newname ) {fSetName = newname;}; - - // Print thpredictions to file - void Print(std::ostream&, bool latex = false) const; - - // Verbosity - static bool Verbose; - }; - - void swap(ThPredictions& lhs, ThPredictions& rhs); -} diff --git a/libnnpdf/src/NNPDF/timer.h b/libnnpdf/src/NNPDF/timer.h deleted file mode 100644 index 160255859f..0000000000 --- a/libnnpdf/src/NNPDF/timer.h +++ /dev/null @@ -1,52 +0,0 @@ -// $Id: timer.h 2079 2014-11-12 11:42:33Z s1044006 $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include -#include -#include - -namespace NNPDF -{ - - /** - * \class Timer - * \brief Computes the calculation time. - */ - class Timer { - private: - - timeval startTime; - - public: - - void start(){ - gettimeofday(&startTime, NULL); - } - - double stop(){ - timeval endTime; - long seconds, useconds; - double duration; - - gettimeofday(&endTime, NULL); - - seconds = endTime.tv_sec - startTime.tv_sec; - useconds = endTime.tv_usec - startTime.tv_usec; - - duration = seconds + useconds/1E6f; - - return duration; - } - - static void printTime(double duration){ - printf("%5.6f seconds\n", duration); - } - - }; - -} diff --git a/libnnpdf/src/NNPDF/utils.h b/libnnpdf/src/NNPDF/utils.h deleted file mode 100644 index 00a5151946..0000000000 --- a/libnnpdf/src/NNPDF/utils.h +++ /dev/null @@ -1,217 +0,0 @@ -// $Id: utils.h 2825 2015-05-03 09:14:46Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#pragma once - -#include "common.h" -#include "NNPDF/exceptions.h" -#include -#include -#include -#include -#include -#include - -/** @defgroup utils Utils - * \brief libnnpdf utility functions - */ - -/*! \addtogroup utils - * @{ - */ - -namespace NNPDF -{ - - // This is all so we can call both joinpath({"a", "b", "c"}) which - // requires a special treatment and auto v = vector{"a", "b", "c"}; joinpath(v). See - // https://stackoverflow.com/questions/4757614/why-doesnt-my-template-accept-an-initializer-list#4763493 - namespace - { - template inline std::string joinpath_inner(const T &list) - { - // This is implemented following - // https://github.com/python/cpython/blob/05d68a8bd84cb141be9f9335f5b3540f15a989c4/Lib/posixpath.py#L75 - const auto sep = '/'; - auto path = std::string{""}; - - for (const auto &it : list) { - //TODO: Use this in C++17 - //if (std::empty(it)) { - if (it.empty()) { - continue; - } - if (*std::cbegin(it) == sep) { - path = it; - } else if (path.empty() || *std::crbegin(path) == sep) { - path += it; - } else { - path += sep + it; - } - } - return path; - } - } - - /** - * @brief Produce a UNIX path from a container of string-like objects. - * @param A iterable of objects that can be joined with strings and iterated over. - * @return A string representing a joined path. - * - * The resuts should be the same as os.path.join in Python (on POSIX), except - * that an empty string is returned from an empty iterable. - */ - template std::string joinpath(const T &list) - { - return joinpath_inner(list); - } - std::string joinpath(const std::initializer_list &list); - - /** - * @brief write_to_file Writes string to file handling failures. - * - * This function creates and writes data using a POSIX file - * descriptor. The data buffer is written to the filesystem through - * write and fsync methods. The implementation also throws - * exceptions at multiple layers (opening, writing, flushing and - * closing the file). This is quite useful on clusters like - * LXPLUS/EOS where IO issues are non negligible. The downside of - * this method is the requirement of loading in memory the entire - * buffer as a string (and obviously the portability to non POSIX - * systems). - * - * @param data the data string - * @param filename the file name - */ - void write_to_file(std::string const& filename, const std::string &data); - - /** - * @brief untargz Decompress tar.gz files in istream. - * @param filename the input filename - * @return the istream object - */ - std::string untargz(std::string const& filename); - - /** - * @brief targz Store to disk the compressed data - * @param filename the output file - * @param data the stream buffer - */ - void targz(std::string const& filename, std::string const& data); - - // ******************* SWIG helpers ***************************** - - /*! - * \class istream_proxy - * \brief Wrapper converting a filename to an std::ifstream for use with the SWIG - */ - class istream_proxy - { - private: - std::ifstream fStream; //!< Internal ifstream object - public: - istream_proxy(std::string const& filename): //!< Constructor takes a filename and initialises the internal ifstream - fStream(filename.c_str()) {} - - std::istream& stream() {return fStream;} //!< Returns the internal ifstream - }; - - /*! - * \class ostream_proxy - * \brief Wrapper converting a filename to an std::ofstream for use with the SWIG - */ - class ostream_proxy - { - private: - std::ofstream fStream; //!< Internal ofstream object - public: - ostream_proxy(std::string const& filename): //!< Constructor takes a filename and initialises the internal ofstream - fStream(filename.c_str()) {} - - std::ostream& stream() {return fStream;} //!< Returns the internal stream - }; - - /** - * @brief The matrix class for the covmat related objects - */ - template - class matrix - { - public: - //!< matrix constructor - matrix(size_t row = 0, size_t col = 0): _size{{row,col}} - { - if (row*col != 0) - _data.resize(row*col); - } - - //!< resize matrix and fill with v - void resize(size_t row, size_t col, T v) - { - _size = {{row,col}}; - _data.resize(row*col, v); - } - - //!< clear matrix size and content - void clear() - { - _size = {0,0}; - _data.clear(); - } - - //operators - 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]; } - - // There is no doubt a better way to do this - std::vector operator*(std::vector in) const { - if (_size[0] != in.size()) - throw RangeError("matrix-vector product", "Mismatch of matrix and input vector dimension"); - std::vector out(_size[0],0); - for (size_t i=0; i<_size[0]; i++) - for (size_t j=0; j<_size[1]; j++) - out[i] += _data[i*_size[1]+j]*in[j]; - return out; - } - - // 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 - std::vector _data; //!< the data array - }; - - // ******************* Numerical ***************************** - double integrate(double data[], size_t npoints, double h); //!< Basic simpson rule integrator - - // ******************* std::string Tokenisers ********************** - std::vector split(std::string const& input); //!< Split std::strings, return std::string std::vector. - void split(std::vector& results, std::string const& input); - - std::vector rsplit(std::string const& input); //!< Split std::strings, return real std::vector. - void rsplit(std::vector& results, std::string const& input); - - std::vector isplit(std::string const& input); //!< Split std::strings, return int std::vector. - void isplit(std::vector& results, std::string const& input); - - // ******************* Basic Stat Functions ********************* - real ComputeAVG(int const& n, const real *x); //!< Compute average from x points - real ComputeAVG(std::vector const& x); //!< Compute average from std::vector x - real ComputeStdDev(int const& n, const real *x); //!< Compute the std deviation - real ComputeStdDev(std::vector const& x); //!< Compute the std deviation - real ComputeEigErr(int const& p, const real *x); //!< Compute error in the Hessian method for PDFs - real ComputeSymEigErr(int const& p, const real *x); //!< Compute error in the symmetric Hessian method for PDFs - 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. -} - -/*! @} */ diff --git a/libnnpdf/src/chisquared.cc b/libnnpdf/src/chisquared.cc deleted file mode 100644 index 78ce0e4aa7..0000000000 --- a/libnnpdf/src/chisquared.cc +++ /dev/null @@ -1,276 +0,0 @@ -// chisquared.cc -// -// NNPDF++ 2012-2015 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it - -#include -#include -#include -#include -#include -#include -#include - -#include "NNPDF/chisquared.h" -#include "NNPDF/pdfset.h" -#include "NNPDF/thpredictions.h" -#include "NNPDF/utils.h" -#include "NNPDF/randomgenerator.h" -#include "NNPDF/exceptions.h" - - -#include "gsl/gsl_matrix.h" -#include "gsl/gsl_linalg.h" - -using namespace std; -namespace NNPDF -{ - - /** - * Generate covariance matrix from basic quantities. - */ - matrix ComputeCovMat_basic(int const ndat, - int const nsys, - std::vector const& sqrt_weights, - std::vector const& central_values, - std::vector const& stat_error, - sysError** const systematic_errors, - bool const mult_errors, - bool const use_theory_errors, - bool const th_cov_matrix, - std::string filename, - std::vector bmask) - { - if (central_values.size() != stat_error.size()) - throw LengthError("ComputeCovMat_basic","mismatch in points between central_values and stat_error!"); - - 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; - - // Statistical error - if (i == j) - sig += pow(stat_error[i],2); - - for (int l = 0; l < nsys; l++) - { - sysError const& isys = systematic_errors[i][l]; - sysError const& jsys = systematic_errors[j][l]; - if (isys.name != jsys.name) - throw RuntimeException("ComputeCovMat", "Inconsistent naming of systematics"); - if (isys.name == "SKIP") - continue; // Continue if systype is skipped - if ((isys.name == "THEORYCORR" || isys.name == "THEORYUNCORR") && !use_theory_errors) - continue; // Continue if systype is theoretical and use_theory_errors == false - 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: if (mult_errors) { signor += isys.mult*jsys.mult; break; } - else { continue; } - case UNSET: throw RuntimeException("ComputeCovMat", "UNSET systype encountered"); - } - } - - // Covariance matrix entry - CovMat(i, j) = (sig + signor*central_values[i]*central_values[j]*1e-4); - // Covariance matrix weight - CovMat(i, j) /= sqrt_weights[i]*sqrt_weights[j]; - } - } - - // Add theory uncertainty - if(th_cov_matrix) - { - auto ThCovMat = read_theory_covmat(ndat, filename, bmask); - for (int i = 0; i < ndat; i++) - { - for (int j = 0; j < ndat; j++) - { - CovMat(i, j) += ThCovMat(i, j); - } - } - } - - return CovMat; - } - - /** - * Generate covariance matrix from CommonData and a t0 vector - * This should be deprecated in favour of a version of `Experiment` that does not contain an FK table. - * CommonData should be considered only as an I/O class with all logic belonging to `Experiment`. - * - * The covariance matris is created accounting for multiplicative and MC uncertainty by default. - */ - matrix ComputeCovMat(CommonData const& cd, std::vector const& t0, - const bool th_cov_matrix, // Specify whether or not the theory error should be used - std::string filename, - std::vector bmask, - double weight) - { - const int ndat = cd.GetNData(); - const int nsys = cd.GetNSys(); - - std::vector sqrt_weights(ndat, sqrt(weight)); - std::vector stat_error(ndat, 0); - for (int i=0; i ComputeSqrtMat(matrix const& inmatrix) - { - const size_t n = inmatrix.size(0); - if (n <= 0) - throw LengthError("CholeskyDecomposition","attempting a decomposition of an empty matrix!"); - - // Originally the decomposition of the covariance matrix was found directly. - // However, when the covariance matrix is near singular it is preferable to - // find the decomposition of the corresponding correlation matrix and then - // to rescale this, so this is what is done here. - // See https://www.gnu.org/software/gsl/doc/html/linalg.html#cholesky-decomposition - // for a more detailed discussion of this - std::vector sqrt_diags; - matrix inmatrix_corr(n,n); - - // Find sqrts of diagonal entries of input covariance matrix - for (int i = 0; i < (int) n; i++) - sqrt_diags.push_back(sqrt(inmatrix(i, i))); - - // Find correlation matrix from covariance matrix - for (int i = 0; i < (int) n; i++) - for (int j = 0; j < (int) n; j++) - inmatrix_corr(i, j) = inmatrix(i, j) / (sqrt_diags[i] * sqrt_diags[j]); - - gsl_matrix_const_view inmatrix_view = gsl_matrix_const_view_array(inmatrix_corr.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 - for (int i = 0; i < (int) n; i++) - for (int j = i + 1; j < (int) n; j++) - sqrtmat(i, j) = 0; - - // Rescale the decomposition of correlation matrix by multiplying by sqrts - // of diagonals of the covariance matrix. This then gives the decomposition - // of the covariance matrix - for (int i = 0; i < (int) n; i++) - for (int j = 0; j < (int) n; j++) - sqrtmat(i, j) = sqrtmat(i, j) * sqrt_diags[i]; - - 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++) - for (int i = 0; i < nDat; i++) - { - x[i] = (data[i] - theory[n+nMem*i]); - for (int j = 0; j < i; j++) - x[i] -= L(i, j) * x[j]; - x[i] /= L(i, i); - chi2[n] += x[i]*x[i]; - } - 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); - - -/* -* Reads in covariance matrix for an experiment from pandas dataframe -*/ - matrix read_theory_covmat(int ndata, const std::string filename, std::vector bmask = {}) - { - // open file - ifstream file(filename.c_str()); - - if (!file.good()) - throw FileError("experiments", "Cannot read covariance matrix file from " + filename); - - string line; - const int first_lines_to_skip = 4; //experiment, dataset, id, header - const int first_columns_to_skip = 3; //experiment, dataset, id - - int file_matrix_size = ndata; - if (!bmask.empty()) - { - if(std::accumulate(bmask.begin(), bmask.end(), 0) != ndata) - throw RuntimeException("read_theory_covmat", "wrong mask size."); - file_matrix_size = bmask.size(); - } - - matrix covmat(ndata, ndata); - - for (int i = 0; i < first_lines_to_skip; ++i) - std::getline(file, line); - - int ix = 0; - for (int idat = 0; idat < file_matrix_size; idat++) - { - getline(file, line); - if (bmask.empty()) - { - auto row = split(line); - for (int jdat = 0; jdat < file_matrix_size; jdat++) - covmat(idat, jdat) = stod(row[first_columns_to_skip + jdat]); - } - else - { - if (bmask[idat] == 1) - { - int jx = 0; - auto row = split(line); - for (int jdat = 0; jdat < file_matrix_size; jdat++) - { - if (bmask[jdat] == 1) - { - covmat(ix, jx) = stod(row[first_columns_to_skip + jdat]); - jx++; - } - } - ix++; - } - } - } - - return covmat; - } - -} diff --git a/libnnpdf/src/common.cc b/libnnpdf/src/common.cc deleted file mode 100644 index 07cef7bc46..0000000000 --- a/libnnpdf/src/common.cc +++ /dev/null @@ -1,55 +0,0 @@ -// $Id: common.cc 3319 2015-09-25 12:54:27Z s1044006 $ -// -// NNPDF++ 2013 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include "NNPDF/common.h" -#include "NNPDF/config.h" - -namespace{ - -class NullBuffer : public std::streambuf -{ -public: - int overflow(int c) { return c; } -}; - -std::ostream& get_null_stream() -{ - static NullBuffer nullbuffer; - static std::ostream nullstream(&nullbuffer); - return nullstream; -} - -int Verbosity = 1; - -} - -namespace NNPDF -{ - std::string getVersion() - { - return std::string(VERSION); - } - - void SetVerbosity(int level) - { - Verbosity = level; - } - - int GetVerbosity() - { - return Verbosity; - } - - std::ostream& get_logger() - { - if (Verbosity){ - return std::cout; - } - return get_null_stream(); - } -} diff --git a/libnnpdf/src/commondata.cc b/libnnpdf/src/commondata.cc deleted file mode 100644 index 872f735e2f..0000000000 --- a/libnnpdf/src/commondata.cc +++ /dev/null @@ -1,601 +0,0 @@ -// $Id: commondata.cc 639 2013-03-25 19:08:38Z s1044006 $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "NNPDF/common.h" -#include "NNPDF/commondata.h" -#include "NNPDF/randomgenerator.h" - -namespace NNPDF -{ - - // Kinematics type labels - const CommonData::kinMap CommonData::kinLabel_latex = { - { "DIS", {"$x$","$Q^2 (GeV^2)$","$y$"}}, - { "DYP", {"$y$","$M^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "JET", {"$\\eta$","$p_T^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "DIJET", {"\\eta","$\\m_{1,2} (GeV)","$\\sqrt{s} (GeV)"}}, - { "PHT", {"$\\eta_\\gamma$","$E_{T,\\gamma}^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "INC", {"$0$","$\\mu^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWK_RAP", {"$\\eta/y$","$M^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWK_PT", {"$p_T$ (GeV)","$M^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWK_PTRAP", {"$\\eta/y$","$p_T^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWK_MLL", {"$M_{ll} (GeV)$","$M_{ll}^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWJ_RAP", {"$\\eta/y$","$M^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWJ_PT", {"$p_T (GeV)$","$M^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWJ_PTRAP", {"$\\eta/y$","$p_T^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWJ_JRAP", {"$\\eta/y$","$M^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWJ_JPT", {"$p_T (GeV)$","$M^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "EWJ_MLL", {"$M_{ll} (GeV)$","$M_{ll}^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "HQP_YQQ", {"$y^{QQ}$","$\\mu^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "HQP_MQQ", {"$M^{QQ} (GeV)$","$\\mu^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "HQP_PTQQ", {"$p_T^{QQ} (GeV)$","$\\mu^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "HQP_YQ", {"$y^Q$","$\\mu^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "HQP_PTQ", {"$p_T^Q (GeV)$","$\\mu^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "HIG_RAP", {"$y$","$M_H^2 (GeV^2)$","$\\sqrt{s} (GeV)$"}}, - { "SIA" , {"$z$", "$Q^2 (GeV^2)$", "$y$"}} - }; - - const CommonData::kinMap CommonData::kinLabel = { - { "DIS", {"x","Q2","y"}}, - { "DYP", {"y","M2","sqrts"}}, - { "JET", {"eta","p_T2","sqrts"}}, - { "DIJET", {"eta","m_12","sqrts"}}, - { "PHT", {"eta_gamma","E_{T,gamma}2","sqrts"}}, - { "INC", {"0","mu2","sqrts"}}, - { "EWK_RAP", {"etay","M2","sqrts"}}, - { "EWK_PT", {"p_T","M2","sqrts"}}, - { "EWK_PTRAP", {"etay","p_T2","sqrts"}}, - { "EWK_MLL", {"M_ll","M_ll2","sqrts"}}, - { "EWJ_RAP", {"etay","M2","sqrts"}}, - { "EWJ_PT", {"p_T","M2","sqrt(s)"}}, - { "EWJ_PTRAP", {"etay","p_T2","sqrts"}}, - { "EWJ_JRAP", {"etay","M2","sqrts"}}, - { "EWJ_JPT", {"p_T","M2","sqrts"}}, - { "EWJ_MLL", {"M_ll","M_ll2","sqrts"}}, - { "HQP_YQQ", {"yQQ","mu2","sqrts"}}, - { "HQP_MQQ", {"MQQ","mu2","sqrts"}}, - { "HQP_PTQQ", {"p_TQQ","mu2","sqrts"}}, - { "HQP_YQ", {"yQ","mu2","sqrts"}}, - { "HQP_PTQ", {"p_TQ","mu2","sqrts"}}, - { "HIG_RAP", {"y","M_H2","sqrts"}}, - { "SIA" , {"z", "Q2", "y"}} - }; - - - // Generate a dataInfo struct given a target filename - dataInfo genInfoStruct(std::string const& targetfile, std::string const& sysfile) - { - // Open commondata file - std::ifstream datafile; - datafile.open(targetfile.c_str()); - - if (!datafile.good()) - throw std::invalid_argument("genInfoStruct: Cannot read commondata file from: " + targetfile); - - // Read metadata - std::string SetName; - int nSys, nData; - - datafile >> SetName - >> nSys - >> nData; - - datafile.close(); - - // Setup info struct - const dataInfo infoStruct = { - nData, - nSys, - SetName, - targetfile, - sysfile - }; - - return infoStruct; - } - - std::string extractSysID(std::string const& sysfile) - { - const std::string sysid = sysfile.substr(sysfile.find_last_of("_")+1, sysfile.length()-4-sysfile.find_last_of("_")-1).c_str(); - return sysid; - } - - // Helper function to parse strings into SYSTYPES - sysType parseSYS(std::string const& str) - { - if (str.compare("ADD") == 0) return ADD; - else if (str.compare("MULT") == 0) return MULT; - else - throw std::invalid_argument("parseSYS: Unrecognised systematic: " + str); - - return MULT; - } - - // CommonData base class constructor - CommonData::CommonData(dataInfo const& info): - fSetName(info.SetName), - fNData(info.nData), - fData(new double[fNData]), - fProc(new std::string[fNData]), - fKin1(new double[fNData]), - fKin2(new double[fNData]), - fKin3(new double[fNData]), - fNSys(info.nSys), - fSysId(extractSysID(info.systypeFile)), - fStat(new double[fNData]), - fSys(new sysError*[fNData]) - { - get_logger() << std::endl << "-- Reading COMMONDATA for Dataset: " << fSetName<::quiet_NaN(); - fKin1[i] = std::numeric_limits::quiet_NaN(); - fKin2[i] = std::numeric_limits::quiet_NaN(); - fKin3[i] = std::numeric_limits::quiet_NaN(); - fStat[i] = std::numeric_limits::quiet_NaN(); - fProc[i] = ""; - fSys[i] = new sysError[fNSys]; - } - - // Read commondata - ReadData(info.targetFile, info.systypeFile); - Verify(); - } - - // CommonData base class constructor - CommonData::CommonData(dataInfoRaw const& info): - fSetName(info.SetName), - fNData(info.nData), - fData(new double[fNData]), - fProc(new std::string[fNData]), - fKin1(new double[fNData]), - fKin2(new double[fNData]), - fKin3(new double[fNData]), - fNSys(info.nSys), - fSysId("DEFAULT"), - fStat(new double[fNData]), - fSys(new sysError*[fNData]) - { - get_logger() << std::endl << "-- Reading COMMONDATA for Dataset: " << fSetName<::quiet_NaN(); - fKin1[i] = std::numeric_limits::quiet_NaN(); - fKin2[i] = std::numeric_limits::quiet_NaN(); - fKin3[i] = std::numeric_limits::quiet_NaN(); - fStat[i] = std::numeric_limits::quiet_NaN(); - fProc[i] = info.ProcType; - fSys[i] = new sysError[fNSys]; - } - } - - // CommonData copy constructor - CommonData::CommonData(const CommonData& set): - fSetName(set.fSetName), - fNData(set.fNData), - fData(new double[fNData]), - fProc(new std::string[fNData]), - fKin1(new double[fNData]), - fKin2(new double[fNData]), - fKin3(new double[fNData]), - fNSys(set.fNSys), - fSysId(set.fSysId), - fStat(new double[fNData]), - fSys(new sysError*[fNData]) - { - // Setup masked data - for (int i = 0; i < fNData; i++) - { - fData[i] = set.fData[i]; - fStat[i] = set.fStat[i]; - - fSys[i] = new sysError[fNSys]; - - fKin1[i] = set.fKin1[i]; - fKin2[i] = set.fKin2[i]; - fKin3[i] = set.fKin3[i]; - - fProc[i] = set.fProc[i]; - VerifyProc(fProc[i]); - - for (int l = 0; l < fNSys; l++) - fSys[i][l] = set.fSys[i][l]; - } - - Verify(); - } - - - void swap(CommonData & lhs, CommonData & rhs) { - using std::swap; - swap(lhs.fSetName, rhs.fSetName); - swap(lhs.fNData, rhs.fNData); - swap(lhs.fData, rhs.fData); - swap(lhs.fProc, rhs.fProc); - swap(lhs.fKin1, rhs.fKin1); - swap(lhs.fKin2, rhs.fKin2); - swap(lhs.fKin3, rhs.fKin3); - swap(lhs.fNSys, rhs.fNSys); - swap(lhs.fSysId, rhs.fSysId); - swap(lhs.fStat, rhs.fStat); - swap(lhs.fSys, rhs.fSys); - } - - CommonData & CommonData::operator =(CommonData other){ - using std::swap; - swap(*this, other); - return *this; - } - - - CommonData::CommonData(CommonData&& other): - fSetName(std::string()), - fNData(0), //Set this to maintain the class invariant - fData(nullptr), - fProc(nullptr), - fKin1(nullptr), - fKin2(nullptr), - fKin3(nullptr), - fNSys(0), - fSysId("-1"), - fStat(nullptr), - fSys(nullptr) - { - using std::swap; - swap(*this, other); - Verify(); - } - - // CommonData masked copy constructor - CommonData::CommonData(const CommonData& set, std::vector const& mask): - fSetName(set.fSetName), - fNData(mask.size()), - fData(new double[fNData]), - fProc(new std::string[fNData]), - fKin1(new double[fNData]), - fKin2(new double[fNData]), - fKin3(new double[fNData]), - fNSys(set.fNSys), - fSysId(set.fSysId), - fStat(new double[fNData]), - fSys(new sysError*[fNData]) - { - // Setup masked data - for (int i = 0; i < fNData; i++) - { - fData[i] = set.fData[mask[i]]; - fStat[i] = set.fStat[mask[i]]; - - fSys[i] = new sysError[fNSys]; - - fKin1[i] = set.fKin1[mask[i]]; - fKin2[i] = set.fKin2[mask[i]]; - fKin3[i] = set.fKin3[mask[i]]; - - fProc[i] = set.fProc[mask[i]]; - VerifyProc(fProc[i]); - - for (int l = 0; l < fNSys; l++) - fSys[i][l] = set.fSys[mask[i]][l]; - } - - Verify(); - } - - // Destructor - CommonData::~CommonData() - { - for (int i=0; i> setname - >> nsys - >> ndata; - - // Verification - if (setname != fSetName) - throw std::runtime_error("CommonData::ReadData: Setname Mismatch."); - - if (nsys != fNSys) - throw std::runtime_error("CommonData::ReadData: N_Uncertainty Mismatch."); - - if (ndata != fNData) - throw std::runtime_error("CommonData::ReadData: NData Mismatch."); - - // Read data - for (int i=0; i> idat - >> fProc[i] - >> fKin1[i] - >> fKin2[i] - >> fKin3[i] - >> fData[i] - >> fStat[i]; - - VerifyProc(fProc[i]); - - if (idat != i+1) - throw std::runtime_error("CommonData::ReadData: Datapoint Mismatch."); - - // Read systematics - for (int l = 0; l < fNSys; l++) - datafile >> fSys[i][l].add >> fSys[i][l].mult; - - } - - datafile.close(); - - // ************ Reading Systypes ******************* - - std::ifstream h; - h.open(systypefilename.c_str()); - - if (h.fail()) - throw std::runtime_error("CommonData::ReadData: Error opening systype file "+systypefilename); - - // Verify number of systematics in SYSTYPE file adds up - h >> nsys; - if (nsys != fNSys) - throw std::runtime_error("CommonData::ReadData: Number of systematics for " + fSetName + " doesn't match"); - - // Read types and names for systematics - for (int l = 0; l < fNSys; l++) - { - int sysnum; std::string tystr,sysname; - h >> sysnum >> tystr >> sysname; - - const bool random = (tystr.compare("RAND") == 0); - const sysType type = random ? MULT:parseSYS(tystr); - - for (int i = 0; i < fNData; i++) - { - fSys[i][l].name = sysname; - fSys[i][l].isRAND = random; - fSys[i][l].type = type; - } - } - - h.close(); - - Verify(); - get_logger() << "-- COMMONDATA Files for "< -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "NNPDF/dataset.h" -#include "NNPDF/chisquared.h" -#include "NNPDF/fastkernel.h" -#include "NNPDF/thpredictions.h" -#include "NNPDF/randomgenerator.h" -#include "NNPDF/exceptions.h" - -using namespace NNPDF; - -/** - * Principal dataset constructor - * \param data NNPDF::CommonData containing the experimental data - * \param set NNPDF::FKSet containing the corresponding theory calculation - * \param weight the factor by which the importance of the dataset in the fit chi² - * is increased. - */ -DataSet::DataSet(CommonData const& data, FKSet const& set, double weight): - CommonData(data), - FKSet(set), - fIsArtificial(false), - fIsT0(false), - fWeight(weight) -{ - fT0Pred.reserve(fNData); - - // Init T0 Vector - for (int i = 0; i < fNData; i++) - fT0Pred.push_back(fData[i]); // Default T0 to data -} - -DataSet::DataSet(const DataSet& set, std::vector const& mask): - CommonData(set, mask), - FKSet(set, mask), - fIsArtificial(set.fIsArtificial), - fIsT0(set.fIsT0), - fWeight(set.fWeight) -{ - fT0Pred.reserve(fNData); - - // Building Obs array - for (int i = 0; i < fNData; i++) - fT0Pred.push_back(set.fT0Pred[mask[i]]); -} - -/** - * Cleanup memory - */ -DataSet::~DataSet() -{ -} - -/** - * Generate covariance matrix and inverse - */ -void DataSet::GenCovMat() const -{ - fCovMat = ComputeCovMat(*this, fT0Pred, false, " ", {}, fWeight); - fSqrtCov = ComputeSqrtMat(fCovMat); -} - -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; - } - } -} - -void DataSet::SetT0(ThPredictions const& t0pred) -{ - // switch flag - fIsT0 = true; - - if (fNData != t0pred.GetNData()) - throw RangeError("DataSet::SetT0","number of datapoints in set and predictions do not match!"); - - for (int i=0; i &DataSet::GetCovMat() const -{ - if (!fCovMat.size(0)) GenCovMat(); - return fCovMat; -} - -/** - * @brief DataSet::GetSqrtCov - * @return - */ -matrix const& DataSet::GetSqrtCov() const -{ - if (!fSqrtCov.size(0)) GenCovMat(); - return fSqrtCov; -} - -/** - * @brief DataSet::GetSqrtCov - * @param i - * @param j - * @return - */ -double DataSet::GetSqrtCov(int i, int j) const -{ - if (!fSqrtCov.size(0)) GenCovMat(); - return fSqrtCov(i, j); -} - -/** - * Update data values - used by MakeArtificial and MakeReplica - */ -void DataSet::UpdateData(double *newdata) // Update data only -{ - for (int i = 0; i < fNData; i++) - fData[i] = newdata[i]; // Set new data values -} - -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]; - } -} diff --git a/libnnpdf/src/experiments.cc b/libnnpdf/src/experiments.cc deleted file mode 100644 index 1ded3486bb..0000000000 --- a/libnnpdf/src/experiments.cc +++ /dev/null @@ -1,659 +0,0 @@ -// $Id: experiments.cc 3222 2015-09-09 16:52:06Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "NNPDF/experiments.h" -#include "NNPDF/chisquared.h" -#include "NNPDF/pdfset.h" -#include "NNPDF/dataset.h" -#include "NNPDF/thpredictions.h" -#include "NNPDF/utils.h" -#include "NNPDF/randomgenerator.h" -#include "NNPDF/exceptions.h" - -using namespace std; -namespace NNPDF{ - -/** - * Constructor - */ -Experiment::Experiment(std::vector const& sets, std::string const& expname): -fExpName(expname), -fNData(0), -fNSys(0), -fSys(NULL), -fSetSysMap(NULL), -fIsArtificial(false), -fIsClosure(false), -fIsT0(false) -{ - // Check - if (sets.size() == 0) - throw LengthError("Experiment::Experiment","No DataSets in constructor!"); - - fIsT0 = sets[0].IsT0(); - - // Create datasets - for (size_t i=0; i< sets.size(); i++) - { - fSets.push_back(sets[i]); - if (sets[i].IsT0() != fIsT0) - throw UserError("Experiment::Experiment", "Supplied datasets must all be either T0 or EXP"); - } - - // Pull data from datasets - PullData(); - GenCovMat(); -} - -/** - * Copy-Constructor - */ -Experiment::Experiment(Experiment const& exp): -fExpName(exp.fExpName), -fNData(exp.fNData), -fNSys(exp.fNSys), -fCovMat(exp.fCovMat), -fSqrtCov(exp.fSqrtCov), -fSamplingMatrix(exp.fSamplingMatrix), -fSys(NULL), -fSetSysMap(NULL), -fIsArtificial(exp.fIsArtificial), -fIsClosure(exp.fIsClosure), -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(); -} - -/** - * Copy-Constructor - */ -Experiment::Experiment(Experiment const& exp, std::vector const& sets): -fExpName(exp.fExpName), -fNData(exp.fNData), -fNSys(exp.fNSys), -fSys(NULL), -fSetSysMap(NULL), -fIsArtificial(exp.fIsArtificial), -fIsClosure(exp.fIsClosure), -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(); -} - -/** - * Destructor - */ -Experiment::~Experiment() -{ - ClearLocalData(); -} - -/* - * Clears data pulled from DataSets - */ -void Experiment::ClearLocalData() -{ - // Already clear - if (fData.size() == 0) - return; - fData.clear(); - - 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; -} - -/** - * @brief Experiment::MakeReplica - * Produces the shifts by using a MC sampling - * Modifies the experiment and correlated datasets - */ -void Experiment::MakeReplica() -{ - cout << "-- Generating replica data for " << fExpName << endl; - - if (fNData == 0) - throw RangeError("Experiment::MakeReplica","you must run ReadData before making a replica."); - - RandomGenerator* rng = RandomGenerator::GetRNG(); - - double *rand = new double[fNSys]; - double *xnor = new double[fNData]; - sysType rST[2] = {ADD,MULT}; - - // Compute the sampling covariance matrix with data CVs, no multiplicative error and no theory errors - if (fSamplingMatrix.size(0) == 0) - { - matrix SM = ComputeCovMat_basic(fNData, fNSys, fSqrtWeights, fData, fStat, fSys, false, false, false, "", {}); - fSamplingMatrix = ComputeSqrtMat(SM); // Take the sqrt of the sampling matrix - } - - // generate procType array for ease of checking - std::vector proctype; - for (int s = 0; s < GetNSet(); s++) - for (int i=0; i< GetSet(s).GetNData(); i++) - proctype.push_back(GetSet(s).GetProc(i)); - - bool isArtNegative = true; - while (isArtNegative) - { - isArtNegative = false; - - for (int l = 0; l < fNSys; l++) - { - rand[l] = rng->GetRandomGausDev(1.0); - if (fSys[0][l].isRAND) - fSys[0][l].type = rST[rng->GetRandomUniform(2)]; - for (int i = 1; i < fNData; i++) - fSys[i][l].type = fSys[0][l].type; - } - - // Generate normal deviates - vector deviates(fNData, std::numeric_limits::quiet_NaN()); - generate(deviates.begin(), deviates.end(), - []()->double {return RandomGenerator::GetRNG()->GetRandomGausDev(1); } ); - const vector correlated_deviates = fSamplingMatrix*deviates; - - // Generate additive theory noise directly from the covariance matrix - vector artdata(fData); - for (int i=0; iGetRandomGausDev(1.0)*fSys[i][l].mult*1e-2); break; - case UNSET: throw RuntimeException("Experiment::MakeReplica", "UNSET systype encountered"); - } - } - else // Noise from correlated systematics - { - switch (fSys[i][l].type) - { - case ADD: break; - case MULT: xnor[i] *= (1.0 + rand[l]*fSys[i][l].mult*1e-2); break; - case UNSET: throw RuntimeException("Experiment::MakeReplica", "UNSET systype encountered"); - } - } - } - - artdata[i] = xnor[i] * artdata[i]; - - } - - // If it's not a closure test, check for positivity of artifical data - if (!fIsClosure) - for (int i=0; i SM = ComputeCovMat_basic(fNData, fNSys, fSqrtWeights, fData, fStat, fSys, false, false, false, "", {}); - fSamplingMatrix = ComputeSqrtMat(SM); // Take the sqrt of the sampling matrix - } - - // generate procType array for ease of checking - std::vector proctype; - for (int s = 0; s < GetNSet(); s++) - for (int i=0; i< GetSet(s).GetNData(); i++) - proctype.push_back(GetSet(s).GetProc(i)); - - bool isArtNegative = true; - while (isArtNegative) - { - isArtNegative = false; - - for (int l = 0; l < fNSys; l++) - { - rand[l] = rng->GetRandomGausDev(1.0); - if (fSys[0][l].isRAND) - fSys[0][l].type = rST[rng->GetRandomUniform(2)]; - for (int i = 1; i < fNData; i++) - fSys[i][l].type = fSys[0][l].type; - } - - // Generate normal deviates - vector deviates(fNData, std::numeric_limits::quiet_NaN()); - generate(deviates.begin(), deviates.end(), - []()->double {return RandomGenerator::GetRNG()->GetRandomGausDev(1); } ); - const vector correlated_deviates = fSamplingMatrix*deviates; - - // Generate additive theory noise directly from the covariance matrix - vector artdata(fData); - - artdata[pointindex] += correlated_deviates[pointindex]; - - // Generation of the experimental noise - xnor[pointindex] = 1.0; - - for (int l = 0; l < fNSys; l++) - { - if (fSys[pointindex][l].name.compare("THEORYCORR")==0) continue; // Skip theoretical uncertainties - if (fSys[pointindex][l].name.compare("THEORYUNCORR")==0) continue; // Skip theoretical uncertainties - if (fSys[pointindex][l].name.compare("SKIP")==0) continue; // Skip uncertainties - if (fSys[pointindex][l].name.compare("UNCORR")==0) // Noise from uncorrelated systematics - { - switch (fSys[pointindex][l].type) - { - case ADD: break; - case MULT: xnor[pointindex] *= (1.0 + rng->GetRandomGausDev(1.0)*fSys[pointindex][l].mult*1e-2); break; - case UNSET: throw RuntimeException("Experiment::MakeReplica", "UNSET systype encountered"); - } - } - else // Noise from correlated systematics - { - switch (fSys[pointindex][l].type) - { - case ADD: break; - case MULT: xnor[pointindex] *= (1.0 + rand[l]*fSys[pointindex][l].mult*1e-2); break; - case UNSET: throw RuntimeException("Experiment::MakeReplica", "UNSET systype encountered"); - } - } - } - - artdata[pointindex] = xnor[pointindex] * artdata[pointindex]; - - - - // If it's not a closure test, check for positivity of artifical data - if (!fIsClosure) - if (artdata[pointindex] < 0 ) - { - // If it's negative and not an asymmetry - const bool is_asymmetry = proctype[pointindex].find("ASY") != std::string::npos; - if (!is_asymmetry) - { - cout << "Datapoint " << pointindex << " has been set to 0" << endl; - isArtNegative = true; - break; - } - } - - // Don't ever discard the replica now - if (isArtNegative) continue; - - // Update data in set - int index = 0; - for (int s = 0; s < GetNSet(); s++) - { - double* data_pointer = artdata.data() + index; - fSets[s].UpdateData(data_pointer); - fSets[s].SetArtificial(true); - index+=fSets[s].GetNData(); - } - } - - // Update local data - PullData(); - - // Now the fData is artificial - fIsArtificial = true; -} - - -void Experiment::MakeClosure(const vector& predictions, bool const& noise) -{ - cout << "-- Generating closure data for " << fExpName << endl; - - // Set closure flag - fIsClosure = true; - - for (size_t s = 0; s < predictions.size(); s++) - { - auto & theory = predictions[s]; - auto & set = fSets[s]; - auto newdata = vector(set.GetNData()); - - for (int i = 0; i < set.GetNData(); i++) - newdata[i] = theory.GetObsCV(i); - - set.UpdateData(newdata.data()); // MakeClosure treated as shifts rather than normalisations - - } - - // Rebuild uncertainty breakdowns - PullData(); - - // Add fluctations in data according to uncertainties. - if (noise) - MakeReplica(); - - // Rebuild covariance matrix - GenCovMat(); -} - -void Experiment::MakeClosure(PDFSet* pdf, bool const& noise) -{ - vector predictions; - for (auto& set: fSets) - { - predictions.emplace_back(pdf,&set); - - } - MakeClosure(predictions, noise); -} - -/** - * Pulls data from held datasets and allocates covariance matrices - */ -void Experiment::PullData() -{ - // Clear local arrays - ClearLocalData(); - fNData = 0; fNSys = 0; - - // List of types that do not need to be correlated across the experiment - const std::array special_types = { - "UNCORR", - "CORR", - "THEORYUNCORR", - "THEORYCORR", - "SKIP" - }; - - // We need to correlate all subset systematics that share a name. - // Therefore we now loop over the subsets, and form a total systematics vector, - // taking into account these correlations. Futhermore we build a map (fSetSysMap) - // which maps each dataset systematic to a correponding systematic in the experiment - // (i.e an index of total_systematics) - vector total_systematics; - fSetSysMap = new int*[GetNSet()]; - for ( int i = 0; i< GetNSet(); i++ ) - { - DataSet& subset = fSets[i]; - fNData += subset.GetNData(); // Count number of datapoints - fSetSysMap[i] = new int[subset.GetNSys()]; // Initialise map - for (int l = 0; l < subset.GetNSys(); l++) - { - const sysError& testsys = subset.GetSys(0,l); - // Check for systematics that are uncorrelated across the experiment - 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(); - total_systematics.emplace_back(testsys); - } - 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 ) - { - // This already exists in total_systematics, it's not a new systematic - found_systematic = true; - 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"); - break; - } - } - // If the systematic doesn't already exist in the list, add it - if ( found_systematic == false ) - { - fSetSysMap[i][l] = total_systematics.size(); - total_systematics.emplace_back(testsys); - } - } - } - } - - // Initialise data arrays, use quiet NaNs to indicate mis-filling - fData = vector(fNData, std::numeric_limits::quiet_NaN()); - fStat = vector(fNData, std::numeric_limits::quiet_NaN()); - fT0Pred = vector(fNData, std::numeric_limits::quiet_NaN()); - fSqrtWeights = vector(); - fSqrtWeights.reserve(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++) - for (int i = 0; i < fSets[s].GetNData(); i++) - { - DataSet& subset = fSets[s]; - fData[idat] = subset.GetData(i); - fStat[idat] = subset.GetStat(i); - fT0Pred[idat] = subset.GetT0Pred(i); - //The cast is so that it does not complain about an ambiguous resolution, - //since apparently there is sqrt(float) and sqrt(long double) but not sqrt(double) - fSqrtWeights.push_back(std::sqrt((long double)subset.GetWeight())); - - // Loop over experimental systematics, find if there is a corresponding dataset systematic - // and set it accordingly. - fSys[idat] = new sysError[fNSys]; - for (int l = 0; l < fNSys; l++) - { - fSys[idat][l].name = total_systematics[l].name; - fSys[idat][l].type = total_systematics[l].type; - fSys[idat][l].isRAND = total_systematics[l].isRAND; - // Find the dataset systematic corresponding to experimental systematic 'l' - // If it doesn't exist (dataset 'subset' doesn't have this systematic) this - // will be subset.GetNSys() - int* map = fSetSysMap[s]; - const int dataset_systematic = find(map, map+subset.GetNSys(), l) - map; - if (dataset_systematic == subset.GetNSys()) - { - // This subset doesn't have this systematic, set it to zero - fSys[idat][l].mult = 0; - fSys[idat][l].add = 0; - } else { - // Copy the systematic - fSys[idat][l].mult = subset.GetSys(i, dataset_systematic).mult; - fSys[idat][l].add = subset.GetSys(i, dataset_systematic).add; - } - } - idat++; - } - - return; -} - -/** - * Generate covariance matrix and inverse - */ - -void Experiment::GenCovMat() -{ - // Compute the Covariance matrix with t0 and NNPDF3.1 theory errors - fCovMat = ComputeCovMat_basic(fNData, fNSys, fSqrtWeights, fT0Pred, fStat, fSys, true, true, false, " ", {}); - fSqrtCov = ComputeSqrtMat(fCovMat); -} - - -/** -* Read in covariance matrix for replica generation from file, and generate covariance matrix and its square root -*/ -void Experiment::LoadRepCovMat(string filename, bool ThUnc, std::vector bmask) -{ - fSamplingMatrix = ComputeCovMat_basic(fNData, fNSys, fSqrtWeights, fT0Pred, fStat, fSys, false, false, ThUnc, filename, bmask); - fSamplingMatrix = ComputeSqrtMat(fSamplingMatrix); -} - -/** -* Read in covariance matrix to be used in fit from file, and generate covariance matrix and its square root -*/ -void Experiment::LoadFitCovMat(string filename, bool ThUnc, std::vector bmask) -{ - fCovMat = ComputeCovMat_basic(fNData, fNSys, fSqrtWeights, fT0Pred, fStat, fSys, true, true, ThUnc, filename, bmask); - fSqrtCov = ComputeSqrtMat(fCovMat); -} - -matrix Experiment::GetSqrtFitCovMat(string filename, std::vector bmask) -{ - matrix const CovMat = ComputeCovMat_basic(fNData, fNSys, fSqrtWeights, fT0Pred, fStat, fSys, true, true, true, filename, bmask); - matrix SqrtCovMat = ComputeSqrtMat(CovMat); - return SqrtCovMat; -} - -void Experiment::ExportCovMat(string filename) -{ - ofstream outCovMat(filename.c_str()); - outCovMat << setprecision(5); - outCovMat << scientific; - for (int i=0; i pseudodata(vector const &exps, - unsigned long int dataseed, int replica) -{ - // make a copy of the experiments - vector output; - output.reserve(exps.size()); - for (auto &e : exps) { - output.emplace_back(new Experiment(*e)); - } - - // select the appropriate random seed, using dataseed - // as initial condition and replica for the filtering selection - unsigned long int seed = 0; - RandomGenerator::GetRNG()->SetSeed(dataseed); - for (int i = 0; i < replica; i++) { - seed = RandomGenerator::GetRNG()->GetRandomInt(); - } - RandomGenerator::GetRNG()->SetSeed(seed); - for (auto e : output) - { - // take exps and MakeReplica - e->MakeReplica(); - - // keep rng flow in sync with nnfit by calling tr/val random seed shuffle - for (auto &set : e->DataSets()) { - // Creating Masks - vector mask(set.GetNData()); - std::iota(mask.begin(), mask.end(), 0); - RandomGenerator::GetRNG()->ShuffleVector(mask); - } - } - - return output; -} - -} // namespace: NNPDF diff --git a/libnnpdf/src/fastkernel.cc b/libnnpdf/src/fastkernel.cc deleted file mode 100644 index 3377eb8cf0..0000000000 --- a/libnnpdf/src/fastkernel.cc +++ /dev/null @@ -1,899 +0,0 @@ -// $Id$ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "NNPDF/fastkernel.h" -#include "NNPDF/utils.h" -#include "NNPDF/exceptions.h" - -namespace NNPDF -{ - - // ********************************************************* - - // Section delineators for FK headers - static const int FK_DELIN_SEC = std::char_traits::to_int_type('_'); - static const int FK_DELIN_BLB = std::char_traits::to_int_type('{'); - static const int FK_DELIN_KEY = std::char_traits::to_int_type('*'); - // ********************************************************* - - /** - * @brief Prototype Constructor for FKParser class - to be used only for constructing new FK tables - * @param filename The FK table filename - */ - FKHeader::FKHeader() - { - - } - - /** - * @brief Constructor for FK header parsing class - * @param str The input stream - */ - FKHeader::FKHeader(std::istream& str) - { - Read(str); - } - - /** - * @brief Constructor for FK header parsing class - * @param filename The FK table filename - */ - FKHeader::FKHeader(std::string const& filename) - { - std::istringstream is(untargz(filename)); - Read(is); - } - - FKHeader::FKHeader(FKHeader const& ref): - fVersions(ref.fVersions), - fGridInfo(ref.fGridInfo), - fTheoryInfo(ref.fTheoryInfo) - { - - } - - FKHeader::~FKHeader() - { - - } - - void FKHeader::Read(std::istream& is) - { - if (!is.good()) - throw FileError("FKHeader::FKHeader","cannot open FK grid file "); - - int peekval = (is >> std::ws).peek(); - while ( peekval == FK_DELIN_SEC || - peekval == FK_DELIN_BLB ) - { - std::string sectionTitle; getline( is, sectionTitle ); - // Remove delineating character and trailing underscores - sectionTitle = sectionTitle.substr(1, sectionTitle.size()); - sectionTitle.erase(std::remove(sectionTitle.begin(), sectionTitle.end(), '_'), sectionTitle.end()); - - // Terminating case - if (sectionTitle.compare("FastKernel") == 0) - return; - - // Blob read - if (peekval == FK_DELIN_BLB) - ReadBlob(is, sectionTitle); - else - ReadKeyValue(is, GetSection(sectionTitle)); - - // Re-peek - peekval = (is >> std::ws).peek(); - } - } - - void FKHeader::Print(std::ostream& out) const - { - out << SectionHeader("GridDesc", BLOB)<(FKHeader::GRIDINFO, "HADRONIC")) - { - for (int i=0; i<14; i++) - { - for (int i=0; i<14; i++) - generalMap << "1 "; - generalMap<(key,value.substr(0, trimpos))); - } - - bool FKHeader::HasTag( section sec, std::string const& key) const - { - const keyMap *tMap = GetMap(sec); - keyMap::const_iterator iMap = (*tMap).find(key); - - if (iMap != (*tMap).end()) - return true; - return false; - } - - std::string FKHeader::GetTag( section sec, std::string const& key) const - { - const keyMap *tMap = GetMap(sec); - keyMap::const_iterator iMap = (*tMap).find(key); - - if (iMap != (*tMap).end()) - return (*iMap).second; - else - throw FileError("FKHeader::GetTag","key " + key + " not found in header!"); - - return std::string(); - } - - FKHeader::section FKHeader::GetSection(std::string const& title) const - { - if (title.compare("VersionInfo") == 0) return VERSIONS; - if (title.compare("GridInfo") == 0) return GRIDINFO; - if (title.compare("TheoryInfo") == 0) return THEORYINFO; - throw FileError("FKHeader::GetSection","Unrecognised section title: " + title); - return BLOB; - } - - const FKHeader::keyMap* FKHeader::GetMap( section const& sec ) const - { - switch (sec) - { - case VERSIONS: return &fVersions; - case GRIDINFO: return &fGridInfo; - case THEORYINFO: return &fTheoryInfo; - case BLOB: return &fBlobString; - } - - return NULL; - } - - void FKHeader::RemTag( section sec, std::string const& key ) - { - keyMap *tMap = GetMap(sec); - keyMap::iterator iTag = (*tMap).find(key); - if (iTag != (*tMap).end()) - (*tMap).erase(iTag); - else - throw FileError("FKHeader::RemTag", "key " + key + " not found in header!"); - } - - - void FKHeader::PrintKeyValue( std::ostream& os, section sec ) const - { - const keyMap *tMap = GetMap(sec); - for (keyMap::const_iterator iPair = (*tMap).begin(); - iPair != (*tMap).end(); iPair++) - os << "*" << (*iPair).first <<": "<<(*iPair).second<> std::ws).peek() == FK_DELIN_KEY ) - { - // Read Key-Value pair - std::string keypair; getline( is, keypair ); - const std::string key = keypair.substr(1, keypair.find(':')-1); - std::string val = keypair.substr(keypair.find(':')+1, keypair.size()); - - // Remove leading spaces from value - const size_t startpos = val.find_first_not_of(' '); - val = val.substr(startpos, val.size()); - - // Add the tag - AddTag(sec, key, val); - } - - return; - - } - - void FKHeader::PrintBlob(std::ostream& os, std::string blobname) const - { - keyMap::const_iterator iBlob = fBlobString.find(blobname); - if (iBlob != fBlobString.end()) - { os << (*iBlob).second << std::endl; } - else - throw InitError("FKHeader::PrintBlob","Blob " + blobname + " not initialised"); - } - - void FKHeader::ReadBlob(std::istream& is, std::string blobname) - { - // While the next char is not a delineator - std::stringstream blobstream; - while ((is >> std::ws).peek() != FK_DELIN_SEC && - (is >> std::ws).peek() != FK_DELIN_BLB ) - { - std::string blobline; getline( is, blobline ); - blobstream << blobline < const& cFactors) - { - std::istringstream is(untargz(filename)); - fFKHeader = FKHeader(is); - fDataName = fFKHeader.GetTag(fFKHeader.GRIDINFO, "SETNAME"); - fDescription = fFKHeader.GetTag(fFKHeader.BLOB, "GridDesc"); - fNData = fFKHeader.GetTag(fFKHeader.GRIDINFO, "NDATA"); - fQ20 = pow(fFKHeader.GetTag(fFKHeader.THEORYINFO, "Q0"), 2); - fHadronic = fFKHeader.GetTag(fFKHeader.GRIDINFO, "HADRONIC"); - fNonZero = parseNonZero(); - fFlmap = fHadronic ? new int[2*fNonZero]:new int[fNonZero]; - fNx = fFKHeader.GetTag(fFKHeader.GRIDINFO, "NX"); - fTx = fHadronic ? fNx*fNx:fNx; - fRmr = fTx*fNonZero % convoluteAlign; - fDSz = fTx*fNonZero + ((convoluteAlign - fRmr)%convoluteAlign) ; - fXgrid = new double[fNx]; - fSigma = new real[fDSz*fNData]; - fHasCFactors = cFactors.size(); - fcFactors = new double[fNData]; - fcUncerts = new double[fNData]; - - InitialiseFromStream(is, cFactors); - } - - /** - * @brief Constructor for FK Table - * @param is An input stream for reading from - * @param cFactors A vector of filenames for potential C-factors - */ - FKTable::FKTable( std::istream& is, std::vector const& cFactors ): - fFKHeader(is), - fDataName( fFKHeader.GetTag (fFKHeader.GRIDINFO, "SETNAME")), - fDescription( fFKHeader.GetTag (fFKHeader.BLOB, "GridDesc")), - fNData( fFKHeader.GetTag (fFKHeader.GRIDINFO, "NDATA")), - fQ20( pow(fFKHeader.GetTag(fFKHeader.THEORYINFO, "Q0"),2)), - fHadronic( fFKHeader.GetTag (fFKHeader.GRIDINFO, "HADRONIC")), - fNonZero(parseNonZero()), // All flavours - fFlmap(fHadronic ? new int[2*fNonZero]:new int[fNonZero]), - fNx( fFKHeader.GetTag (fFKHeader.GRIDINFO, "NX")), - fTx(fHadronic ? fNx*fNx:fNx), - fRmr(fTx*fNonZero % convoluteAlign), - fDSz( fTx*fNonZero + (convoluteAlign - fRmr)%convoluteAlign), - fXgrid(new double[fNx]), - fSigma( new real[fDSz*fNData]), - fHasCFactors(cFactors.size()), - fcFactors(new double[fNData]), - fcUncerts(new double[fNData]) - { - InitialiseFromStream(is, cFactors); - } - - /** - * @brief Method for initialisation from stream - * @param is the input stream after reading the FK header - */ - void FKTable::InitialiseFromStream( std::istream& is, std::vector const& cFactors ) - { - - get_logger() << fFKHeader.SectionHeader(fDataName.c_str(),fFKHeader.GRIDINFO)<> maxMap[i]; - - // Build FLMap - int index = 0; - for (int i=0; i> maxMap[i][j]; - - // Build FLMap - int index = 0; - for (int i=0; i> fXgrid[i]; - - // Sanity checks - if ( fNData <= 0 ) - throw RangeError("FKTable::FKTable","Number of datapoints is set to: " + std::to_string(fNData) ); - - if ( fNx <= 0 ) - throw RangeError("FKTable::FKTable","Number of x-points is set to: " + std::to_string(fNx) ); - - if ( fNonZero <= 0 ) - throw RangeError("FKTable::FKTable","Number of nonzero flavours is set to: " + std::to_string(fNonZero) ); - - get_logger() << fNData << " Data Points " - << fNx << " X points " - << fNonZero << " active flavours"< also zeros pad quantities - for (int i=0; i datasplit; - rsplit(datasplit,line); - const int d = datasplit[0]; // datapoint - const int a = datasplit[1]; // x1 index - const int b = datasplit[2]; // x2 index - - for (int j=0; j datasplit; - rsplit(datasplit,line); - - const int d = datasplit[0]; - const int a = datasplit[1]; - - for (int j=0; j const& mask): - fFKHeader(set.fFKHeader), - fDataName(set.fDataName), - fNData(mask.size()), - fQ20(set.fQ20), - fHadronic(set.fHadronic), - fNonZero(set.fNonZero), - fFlmap(fHadronic ? (new int[2*fNonZero]):(new int[fNonZero])), - fNx(set.fNx), - fTx(set.fTx), - fRmr(set.fRmr), - fDSz(set.fDSz), - fXgrid(new double[fNx]), - fSigma(new real[fDSz*fNData]), - fHasCFactors(set.fHasCFactors), - fcFactors(new double[fNData]), - fcUncerts(new double[fNData]) - { - if (fNData == 0) - throw RangeError("FKTable::FKTable", fDataName + " has 0 points, check your filter mask and your TR/VAL split!"); - - // Copy X grid - for (int i = 0; i < fNx; i++) - fXgrid[i] = set.fXgrid[i]; - - // Copy flavour map - if (IsHadronic()) // Hadronic - { - for (int fl = 0; fl < fNonZero; fl++) - { - fFlmap[2*fl] = set.fFlmap[2*fl]; - fFlmap[2*fl+1] = set.fFlmap[2*fl+1]; - } - } else // DIS - { - for (int fl = 0; fl < fNonZero; fl++) - fFlmap[fl] = set.fFlmap[fl]; - } - - // Copy reduced FK table - for (int i = 0; i < fNData; i++) - { - for (int j = 0; j < fDSz; j++) - fSigma[i*fDSz + j] = set.fSigma[mask[i]*fDSz + j]; - fcFactors[i] = set.GetCFactors()[mask[i]]; - fcUncerts[i] = set.GetCUncerts()[mask[i]]; - } - } - - /** - * @brief FKTable destructor - */ - FKTable::~FKTable() - { - delete[] fSigma; - delete[] fFlmap; - delete[] fXgrid; - delete[] fcFactors; - delete[] fcUncerts; - } - - - /** - * @brief FKTable print to ostream - */ - void FKTable::Print(std::ostream& os) - { - get_logger() << "****** Exporting FKTable: "<> std::ws).peek(); - if (peekval == FK_DELIN_KEY) - nDelin++; - getline(g,line); - } - - // Read C-factors - for (int i = 0; i < fNData; i++) - { - real tmp, tmp_error; - if (! (g >> tmp >> tmp_error ) ) - throw FileError("FKTable::FKTable","Error reading cfactor file: " + cfilename); - fcFactors[i] *= tmp; - fcUncerts[i] += tmp_error*tmp_error; - } - - g.close(); - } - - // GetFKValue returns the appropriate point of the FK table - int FKTable::GetISig( int const& d, // Datapoint index - int const& a, // First x-index - int const& b, // Second x-index - int const& ifl1, // First flavour index - int const& ifl2 // Second flavour index - ) const - { - - if (!fHadronic) - throw EvaluationError("FKTable::GetISig","Hadronic call for DIS table!"); - - // Identify which nonZero flavour ifl1 and ifl2 correspond to - int j = -1; - for (int i=0; i < fNonZero; i++) - if ( ifl1 == fFlmap[2*i] && ifl2 == fFlmap[2*i+1]) - { - j = i; - break; - } - - // Not in FLmap - if (j == -1) - return -1; - - // Not in dataset - if (d >= fNData) - return -1; - - // Not in x-grid - if (a >= fNx || b >= fNx) - return -1; - - // Return pointer to FKTable segment - return d*fDSz+j*fTx+a*fNx+b ; - } - - // DIS version of GetFKValue - int FKTable::GetISig( int const& d, // Datapoint index - int const& a, // x-index - int const& ifl // flavour index - ) const - { - - if (fHadronic) - throw EvaluationError("FKTable::GetISig","DIS call for Hadronic table!"); - - // Identify which nonZero flavour ifl corresponds to - int j = -1; - for (int i=0; i < fNonZero; i++) - if ( ifl == fFlmap[i] ) - { - j = i; - break; - } - - // Not in FLmap - if (j == -1) - return -1; - - // Not in dataset - if (d >= fNData) - return -1; - - // Not in x-grid - if (a >= fNx) - return -1; - - return d*fDSz+j*fNx+a; - } - - bool FKTable::OptimalFlavourmap(std::string& flmap) const //!< Determine and return the optimal flavour map - { - bool nonZero[fNonZero]; - for (int j=0; j < fNonZero; j++) - { - nonZero[j] = false; - for (int d=0; d < fNData; d++) - { - if (nonZero[j]) break; - - for (int a=0; a> iNonZero; - nNonZero += iNonZero; - } - else - { - for (int j=0; j<14; j++) - if (fmBlob.good()) - { - bool iNonZero = false; fmBlob >> iNonZero; - nNonZero += iNonZero; - } - else - throw RangeError("FKTable::parseNonZero","FlavourMap formatting error!"); - } - } - - return nNonZero; - } - -} diff --git a/libnnpdf/src/fkgenerator.cc b/libnnpdf/src/fkgenerator.cc deleted file mode 100644 index 64e559b4fa..0000000000 --- a/libnnpdf/src/fkgenerator.cc +++ /dev/null @@ -1,99 +0,0 @@ -// $Id$ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include -#include -#include - -#include "NNPDF/fkgenerator.h" -#include "NNPDF/utils.h" -#include "NNPDF/exceptions.h" - -namespace NNPDF -{ - - FKGenerator::FKGenerator( std::istream& is ): - FKTable(is), - fAccumulator(new double[fNData*fDSz]()) - { - } - - // Fill the FKTable - beware there is no error checking performed! - void FKGenerator::Fill( int const& d, // Datapoint index - int const& ix1, // First x-index - int const& ix2, // Second x-index - int const& ifl1, // First flavour index - int const& ifl2, // Second flavour index - real const& fk // FK Value - ) - { - - if (fk==0) - return; - - if (d >= fNData) - throw RangeError("FKGenerator::Fill","datapoint " + std::to_string(d) + "out of bounds."); - - if (ix1 >= fNx) - throw RangeError("FKGenerator::Fill","xpoint " + std::to_string(ix1) + " out of bounds."); - - if (ix2 >= fNx) - throw RangeError("FKGenerator::Fill","xpoint " + std::to_string(ix2) + " out of bounds."); - - if (ifl1 >= 14) - throw RangeError("FKGenerator::Fill","flavour " + std::to_string(ifl1) + " out of bounds."); - - if (ifl2 >= 14) - throw RangeError("FKGenerator::Fill","flavour " + std::to_string(ifl2) + " out of bounds."); - - // pointer to FKTable segment - const size_t jFL = 14*ifl1 + ifl2; - const size_t iSig = d*fDSz+jFL*fTx+ix1*fNx+ix2 ; - - // Assign FK Table - fAccumulator[iSig] += fk; - - return; - }; - - // DIS version of Fill - void FKGenerator::Fill( int const& d, // Datapoint index - int const& ix, // x-index - int const& ifl, // flavour index - real const& fk // FK Value - ) - { - - if (fk==0) - return; - - if (d >= fNData) - throw RangeError("FKGenerator::Fill","datapoint " + std::to_string(d) + " out of bounds."); - - if (ix >= fNx) - throw RangeError("FKGenerator::Fill","xpoint " + std::to_string(ix) + " out of bounds."); - - if (ifl >= 14) - throw RangeError("FKGenerator::Fill","flavour " + std::to_string(ifl) + " out of bounds."); - - // pointer to FKTable segment - const size_t iSig = d*fDSz+ifl*fNx+ix; - - // Assign FK Table - fAccumulator[iSig] += fk; - - return; - }; - - void FKGenerator::Finalise() - { - for (int i=0; i -#include -#include -#include -#include -#include -#include - -#include "NNPDF/fkset.h" -#include "NNPDF/fastkernel.h" -#include "NNPDF/utils.h" -#include "NNPDF/exceptions.h" - -namespace NNPDF -{ - - /** - * DataSet Operators - * Ratio -> for W asymmetries etc - */ - static void OpNull(int const& Nvals, std::vector const& obs, real*out) - { - for (int i=0; i const& obs, real*out) - { - for (int i=0; i const& obs, real*out) - { - if (obs.size()!=2) - throw LengthError("OpRatio","number of FK grids is incorrect"); - - for (int i=0; i const& obs, real*out) - { - if (obs.size()!=2) - throw LengthError("OpAsy","number of FK grids is incorrect"); - - for (int i=0; i const& obs, real*out) - { - if (obs.size()!=4) - throw LengthError("OpSmn","number of FK grids is incorrect"); - - for (int i=0; i const& obs, real*out) - { - if (obs.size()!=20) - throw LengthError("OpCom","number of FK grids is incorrect"); - - for (int i=0; i const& obs, real*out) - { - if (obs.size()!=10) - throw LengthError("OpSmt","number of FK grids is incorrect"); - - for (int i=0; i const& fktabs): - fOperator(op), - fNSigma(fktabs.size()), - fNDataFK(fktabs[0]->GetNData()), - fHadronic(fktabs[0]->IsHadronic()), - fDataName(fktabs[0]->GetDataName()), - fFK(new FKTable*[fNSigma]) - { - // Copy FKTables - for (size_t i=0; iIsHadronic() != fHadronic) - throw EvaluationError("FKSet::FKSet", "Hadronic status mismatch!"); - - if (fFK[i]->GetNData() != fNDataFK) - throw EvaluationError("FKSet::FKSet", "NData mismatch!"); - } - }; - - - // FKSet copy-constructor - FKSet::FKSet(FKSet const& set): - fOperator(set.fOperator), - fNSigma(set.fNSigma), - fNDataFK(set.fNDataFK), - fHadronic(set.fHadronic), - fDataName(set.fDataName), - fFK(new FKTable*[fNSigma]) - { - // Copy FKTables - for (int i=0; iGetNData() != fNDataFK) - throw RangeError("FKSet::FKSet", "NData mismatch!"); - }; - - void swap(FKSet & lhs, FKSet & rhs) - { - using std::swap; - swap(lhs.fOperator, rhs.fOperator); - swap(lhs.fNSigma, rhs.fNSigma); - swap(lhs.fNDataFK, rhs.fNDataFK); - swap(lhs.fHadronic, rhs.fHadronic); - swap(lhs.fDataName, rhs.fDataName); - swap(lhs.fFK, rhs.fFK); - } - - FKSet & FKSet::operator =(FKSet other){ - using std::swap; - swap(*this, other); - return *this; - } - - - FKSet::FKSet(FKSet && other): - fOperator(OpNull), - fNSigma(0), - fNDataFK(0), - fHadronic(false), - fDataName(std::string()), - fFK(nullptr) - { - using std::swap; - swap(*this, other); - } - - - // FKSet masked copy-constructor - FKSet::FKSet(FKSet const& set, std::vector const& mask): - fOperator(set.fOperator), - fNSigma(set.fNSigma), - fNDataFK( mask.size() ), - fHadronic(set.fHadronic), - fDataName(set.fDataName), - fFK(new FKTable*[fNSigma]) - { - // Copy FKTables - for (int i=0; iGetNData() != fNDataFK) - throw RangeError("FKSet::FKSet","NData mismatch!"); - }; - - FKSet::~FKSet() - { - for (int i=0; i -#include -#include - -#include "LHAPDF/LHAPDF.h" -#include "NNPDF/lhapdfset.h" -#include "NNPDF/pdfset.h" -#include "NNPDF/utils.h" -#include "NNPDF/exceptions.h" - -using namespace NNPDF; - -/** - * The constructor - */ -LHAPDFSet::LHAPDFSet(std::string const& pdfname, erType etype): -LHAPDF::PDFSet(pdfname), -NNPDF::PDFSet(pdfname, size(), etype) -{ - // Verify erType - const std::string LHError = errorType(); - switch (etype) - { - case erType::ER_NONE: - case erType::ER_MCT0: - fMembers = 1; - break; - - case erType::ER_MC: - case erType::ER_MC68: - if (LHError.compare("replicas") != 0 ) - throw RuntimeException("LHAPDFSet::LHAPDFSet", "Error: ErrorSet Types do not match: ER_MC(68) and " + LHError); - fMembers-=1; // Remove replica zero - break; - - case erType::ER_EIG: - case erType::ER_EIG90: - if (LHError.compare("hessian") != 0 ) - throw RuntimeException("LHAPDFSet::LHAPDFSet", "Error: ErrorSet Types do not match: ER_EIG(90) and " + LHError); - break; - - case erType::ER_SYMEIG: - if (LHError.compare("symmhessian") != 0 ) - throw RuntimeException("LHAPDFSet::LHAPDFSet", "Error: ErrorSet Types do not match: ER_SYMEIG and " + LHError); - break; - } - - // Load member PDFs - if (fMembers == 1) - fMemberPDFs.push_back(mkPDF(0)); - else - mkPDFs( fMemberPDFs ); - - get_logger() << pdfname<< " Initialised with " << fMembers<<" members and errorType "<hasFlavor(pdgid); -} - -void LHAPDFSet::GetPDF(real const& x, real const& Q2, int const& n, real* pdf) const -{ - // Skip replica 0 for MC like ensembles - const int iMem = (fEtype == erType::ER_MC || fEtype == erType::ER_MC68) ? n+1:n; - - if (!fMemberPDFs[iMem]->inPhysicalRangeXQ2(x, Q2)) - throw RangeError("LHAPDFSet::GetPDF","kinematics out of set range: x="+ std::to_string(x) +" Q2="+ std::to_string(Q2)); - - //** PDG CODES **// - // Let's be super explicit here - std::array LHA; - - LHA[TBAR] = fMemberPDFs[iMem]->xfxQ2(-6, x, Q2); - LHA[BBAR] = fMemberPDFs[iMem]->xfxQ2(-5, x, Q2); - LHA[CBAR] = fMemberPDFs[iMem]->xfxQ2(-4, x, Q2); - LHA[SBAR] = fMemberPDFs[iMem]->xfxQ2(-3, x, Q2); - LHA[UBAR] = fMemberPDFs[iMem]->xfxQ2(-2, x, Q2); - LHA[DBAR] = fMemberPDFs[iMem]->xfxQ2(-1, x, Q2); - - LHA[T] = fMemberPDFs[iMem]->xfxQ2(6, x, Q2); - LHA[B] = fMemberPDFs[iMem]->xfxQ2(5, x, Q2); - LHA[C] = fMemberPDFs[iMem]->xfxQ2(4, x, Q2); - LHA[S] = fMemberPDFs[iMem]->xfxQ2(3, x, Q2); - LHA[U] = fMemberPDFs[iMem]->xfxQ2(2, x, Q2); - LHA[D] = fMemberPDFs[iMem]->xfxQ2(1, x, Q2); - - LHA[GLUON] = fMemberPDFs[iMem]->xfxQ2(21, x, Q2); - LHA[PHT] = fMemberPDFs[iMem]->xfxQ2(22, x, Q2); - - // Transform - LHA2EVLN(LHA.data(), pdf); -} - -real LHAPDFSet::xfxQ(real const& x, real const& Q, int const& n, int const& fl) const -{ - // Skip replica 0 for MC like ensembles - const int iMem = (fEtype == erType::ER_MC || fEtype == erType::ER_MC68) ? n+1:n; - - if (!fMemberPDFs[iMem]->inPhysicalRangeXQ(x, Q)) - throw RangeError("LHAPDFSet::GetPDF", "kinematics out of set range: x=" + std::to_string(x) + " Q=" + std::to_string(Q)); - - return fMemberPDFs[iMem]->xfxQ(fl, x, Q); -} - - - - diff --git a/libnnpdf/src/logger.cc b/libnnpdf/src/logger.cc deleted file mode 100644 index f1b4f4af4a..0000000000 --- a/libnnpdf/src/logger.cc +++ /dev/null @@ -1,162 +0,0 @@ -// $Id -// -// NNPDF++ 2013 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include "NNPDF/common.h" -#include "NNPDF/logger.h" -#include "NNPDF/exceptions.h" -#include "NNPDF/utils.h" - -#include -#include -#include -#include -#include - -using std::cout; -using std::cerr; -using std::endl; -using std::string; - -namespace NNPDF -{ - - // Symbol for logmanager singleton - LogManager *LogManager::logInstance = 0; - - /** - * @brief LogManager::LogManager() - * Private constructor for the Log file manager class - */ - LogManager::LogManager(): - fMap(), - fBasePath("") - { - - } - - /** - * @brief LogManager::~LogManager() - * LogManager singleton destructor - */ - LogManager::~LogManager() - { - fMap.clear(); - } - - /** - * @brief LogManager::AddLogger - * Adds a new logger class to the manager. - * Logger classes are kept in a hash map, indexed by the hash of logname - * @param logname The identifier for the log, e.g NNFIT or GAMinimizer - * @param filename The target filename to which the log is written out. Relative to base path in LogManager - */ - void LogManager::AddLogger(string const& logname, string const& filename) - { - LogManager* LM = GetLM(); - - std::hash str_hash; - const size_t hashval = str_hash(logname); - - // Collision - if (LM->fMap.find(hashval)!=LM->fMap.end()) - { - cerr << "LogManager::AddLogger Error - hash collision for new log: \'" <fBasePath + "/" + filename; - LM->fMap.insert(std::pair(hashval,Logger(logname, targetfile ))); - } - - /** - * @brief LogManager::GetLogger - * Returns a reference to specific logger, identified by name - * @param logname The identifier for the log, e.g NNFIT or GAMinimizer - */ - Logger& LogManager::GetLogger(string const& logname) - { - LogManager* LM = GetLM(); - std::hash str_hash; - const size_t hashval = str_hash(logname); - - // Collision - LogMap::iterator iLog = LM->fMap.find(hashval); - if ( iLog == LM->fMap.end()) - { - cerr << "LogManager::GetLogger Error - log: \'" <fMap.begin(); - for (; iLog != GetLM()->fMap.end(); iLog++) - (*iLog).second.Export(); - cout << "-- LogManager:: Log Export Completed"< -#ifdef SSE_CONV -#define NNMPI_REAL MPI_FLOAT -#else -#define NNMPI_REAL MPI_DOUBLE -#endif -#endif - -#define UNUSED(expr) (void)(expr) - -namespace NNPDF { - -#ifdef SSE_CONV - #include - - static inline void convolute(const real* __restrict__ x, const real* __restrict__ y, real& retval, int const& n) - { - __m128 acc = _mm_setzero_ps(); - - const __m128* a = (const __m128*) x; - const __m128* b = (const __m128*) y; - - for (int i=0; iTaskID() == 0) - { - MPI::ResetWorld(); - int v = -1; - for (int d = 1; d < getInstance()->NumTasks(); d++) - MPI_Send(&v, 1, MPI_INT, d, 0, MPI_COMM_WORLD); - } - MPI_Finalize(); -#endif - } - - void MPI::Init() - { -#ifdef OPENMPI - MPI_Init(NULL, NULL); - MPI_Comm_size(MPI_COMM_WORLD,&getInstance()->fworld); - MPI_Comm_rank(MPI_COMM_WORLD,&getInstance()->ftaskid); - getInstance()->fnumbertasks = getInstance()->fworld; - if (getInstance()->ftaskid > 0) getInstance()->SlaveConvoluteLoop(); -#endif - } - - void MPI::DecomposeDomain(int domain_size, int taskid, int *subdomain_start, int *subdomain_size) - { -#ifdef OPENMPI - if (taskid == 0) - if (getInstance()->fnumbertasks > domain_size) - getInstance()->fnumbertasks = domain_size; - - *subdomain_start = domain_size / getInstance()->fnumbertasks * taskid; - *subdomain_size = domain_size / getInstance()->fnumbertasks; - if (taskid == getInstance()->fnumbertasks - 1) - *subdomain_size += domain_size % getInstance()->fnumbertasks; -#else - UNUSED(domain_size); - UNUSED(taskid); - UNUSED(subdomain_start); - UNUSED(subdomain_size); -#endif - } - - void MPI::SlaveConvoluteLoop() - { -#ifdef OPENMPI - MPI_Status status; - int offset, chunksize, Npdf, Dsz, NData; - real *theory; - real *pdf; - real *sigma; - - while(true) - { - // Perform the receiver - MPI_Recv(&offset, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status); - if (offset == -1) break; - - MPI_Recv(&chunksize, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &status); - MPI_Recv(&Npdf, 1, MPI_INT, 0, 2, MPI_COMM_WORLD, &status); - MPI_Recv(&Dsz, 1, MPI_INT, 0, 3, MPI_COMM_WORLD, &status); - MPI_Recv(&NData, 1, MPI_INT, 0, 4, MPI_COMM_WORLD, &status); - - theory = new real[NData*Npdf]; - MPI_Recv(&theory[offset*Npdf], chunksize*Npdf, NNMPI_REAL, 0, 5, MPI_COMM_WORLD, &status); - - pdf = new real[Dsz*Npdf]; - MPI_Recv(pdf, Dsz*Npdf, NNMPI_REAL, 0, 6, MPI_COMM_WORLD, &status); - - sigma = new real[Dsz*NData]; - MPI_Recv(sigma, Dsz*NData, NNMPI_REAL, 0, 7, MPI_COMM_WORLD, &status); - - for (int i = offset; i < (offset+chunksize); i++) - for (int n = 0; n < Npdf; n++) - { - theory[i*Npdf + n] = 0; - convolute(pdf+Dsz*n,sigma+Dsz*i,theory[i*Npdf + n],Dsz); - } - - MPI_Send(&theory[offset*Npdf], chunksize*Npdf, NNMPI_REAL, 0, 0, MPI_COMM_WORLD); - - delete[] theory; - delete[] pdf; - delete[] sigma; - } -#endif - } - - void MPI::SendJobs(int NData, int Npdf, int Dsz, real *theory, real *pdf, real *sigma) - { -#ifdef OPENMPI - int offset, chunksize; - for (int d = 1; d < getInstance()->NumTasks(); d++) - { - MPI::DecomposeDomain(NData, d, &offset, &chunksize); - MPI_Send(&offset, 1, MPI_INT, d, 0, MPI_COMM_WORLD); - MPI_Send(&chunksize, 1, MPI_INT, d, 1, MPI_COMM_WORLD); - - MPI_Send(&Npdf, 1, MPI_INT, d, 2, MPI_COMM_WORLD); - MPI_Send(&Dsz, 1, MPI_INT, d, 3, MPI_COMM_WORLD); - MPI_Send(&NData, 1, MPI_INT, d, 4, MPI_COMM_WORLD); - - MPI_Send(&theory[offset*Npdf],chunksize*Npdf,NNMPI_REAL, d, 5, MPI_COMM_WORLD); - - MPI_Send(pdf, Dsz*Npdf, NNMPI_REAL, d, 6, MPI_COMM_WORLD); - MPI_Send(sigma, Dsz*NData, NNMPI_REAL, d, 7,MPI_COMM_WORLD); - - //printf("Sent %d elements to task %d offset= %d\n",chunksize,d,offset); - } -#else - UNUSED(NData); - UNUSED(Npdf); - UNUSED(Dsz); - UNUSED(theory); - UNUSED(pdf); - UNUSED(sigma); -#endif - } - - void MPI::RecvJobs(int NData, int Npdf, real *theory) - { -#ifdef OPENMPI - MPI_Status status; - int offset, chunksize; - for (int d = 1; d < getInstance()->NumTasks(); d++) - { - MPI::DecomposeDomain(NData, d, &offset, &chunksize); - MPI_Recv(&theory[offset*Npdf], chunksize*Npdf, NNMPI_REAL, d, 0, MPI_COMM_WORLD, &status); - } -#else - UNUSED(NData); - UNUSED(Npdf); - UNUSED(theory); -#endif - } - - void MPI::ResetWorld() - { -#ifdef OPENMPI - getInstance()->fnumbertasks = getInstance()->fworld; -#endif - } - -} diff --git a/libnnpdf/src/nnpdfdb.cc b/libnnpdf/src/nnpdfdb.cc deleted file mode 100644 index a8d4fe47e5..0000000000 --- a/libnnpdf/src/nnpdfdb.cc +++ /dev/null @@ -1,122 +0,0 @@ -/* - * nnpdfdb.cc - * Database access routines for nnpdf applications - * nathan.hartland@physics.ox.ac.uk 02/15 - */ - -#include "NNPDF/nnpdfdb.h" - -namespace NNPDF -{ - - IndexDB::IndexDB(std::string const& _path, std::string const& _table): - fDB(0), - fPath(_path), - fTable(_table), - fNEntries(0) - { - - // Setup db connection - if( sqlite3_open(fPath.c_str(), &fDB) ) - { - std::stringstream err; err << "Can't open database: " << sqlite3_errmsg(fDB); - throw FileError("IndexDB::IndexDB", err.str()); - } - - // Determine number of entries - sqlite3_stmt *statement; - const std::string query = "select * from "+fTable; - const int retcode = sqlite3_prepare_v2(fDB, query.c_str(), -1, &statement, 0 ); - if ( retcode != SQLITE_OK ) - { - std::stringstream err; err << "SQLITE Error: " << sqlite3_errmsg(fDB); - throw RuntimeException("IndexDB::IndexDB", err.str()); - } - - while (sqlite3_step(statement) == SQLITE_ROW) - fNEntries++; - - } - - IndexDB::~IndexDB() - { - sqlite3_close(fDB); - } - - void IndexDB::ExtractMap( const int& id, std::vector const& keys, std::map& map) - { - for (size_t i=0; i(*this, id, keys[i])); - } - - std::vector IndexDB::ExtractString( const int& id, std::vector const& keys) - { - std::map map; - ExtractMap(id,keys,map); - - std::vector v; - for( std::map::iterator it = map.begin(); it != map.end(); ++it ) - v.push_back( it->second ); - - return v; - } - - // Needs to be specialised to avoid splitting the string and returning NULL strings - template<> - std::string dbquery(IndexDB const& db, int const& id, std::string const& field) - { - std::stringstream query; - query << "select " << field << " from "< dbmatch(IndexDB const& db, std::string const& field, std::string const& value) - { - std::stringstream query; - query << "select id from "< indices; - while ( sqlite3_step(statement) == SQLITE_ROW ) - { - std::stringstream s; - s <<(char*)sqlite3_column_text(statement, 0); - - int ID; s >> ID; - indices.push_back(ID); - } - return indices; - } - -} diff --git a/libnnpdf/src/parametrisation.cc b/libnnpdf/src/parametrisation.cc deleted file mode 100644 index e88715837e..0000000000 --- a/libnnpdf/src/parametrisation.cc +++ /dev/null @@ -1,425 +0,0 @@ -// $Id: parametrisation.cc 3195 2015-08-27 10:29:43Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include -#include -#include - -#include - -#include "NNPDF/parametrisation.h" -#include "NNPDF/randomgenerator.h" -#include "NNPDF/exceptions.h" - -using namespace NNPDF; - -// ******************** Base ****************************** -/** - * @brief Parametrisation base class constructor - * This is an abstract class handling parameter arrays and I/O - * @param name a string specifying the type of parametrisation - * @param nParameters the number of parameters requested by the derived class - */ -Parametrisation::Parametrisation( std::string name, int nParameters ): -fNParameters(nParameters), -fParameters(new real[fNParameters]()), -fParamName(name) -{ -} - -/** - * @brief Parametrisation base class copy constructor - * @param o the parametrisation to be copied - * - * This constructor copies the model parameters from another parametrisation. - */ -Parametrisation::Parametrisation( Parametrisation const& o): -fNParameters(o.fNParameters), -fParameters(new real[fNParameters]()), -fParamName(o.fParamName) -{ - CopyPars(&o); - return; -} - -/** - * @brief Parametrisation base class destructor - * - * Frees parametrisation array - */ -Parametrisation::~Parametrisation() -{ - delete[] fParameters; -} - -/** - * @brief Parametrisation::CopyPars - * @param t - */ -void Parametrisation::CopyPars(Parametrisation const* t) -{ - if (fNParameters != t->GetNParameters()) - throw EvaluationError("Parametrisation::CopyPars", "number of parameters does not match: " + std::to_string(fNParameters) + " vs " + std::to_string(t->GetNParameters())); - - for (int i=0; ifParameters[i]; - - return; -} - -/** - * @brief Parametrisation::ReadPars - * @param instream - input for parameters - */ -void Parametrisation::ReadPars(std::ifstream& instream) -{ - for (int i=0; i> fParameters[i]; -} - -/** - * @brief Parametrisation::WritePars - * @param outstream - output for parameters - */ -void Parametrisation::WritePars(std::ofstream& outstream) -{ - for (int i=0; i const& param) -{ - for (int i=0; i const& arch) -{ - int nparam = 0; - for (size_t i=1; i const& arch): -Parametrisation(std::string("MultiLayerPerceptron"), perceptron_count_parameters(arch)), -fNLayers(arch.size()), -fArch(0), -fWeightMatrix(0), -fOutputMatrix(0), -fActFunction(sigmoid) -{ - int pSize = 0; // size of parameter array - int* pMap = new int[fNLayers-1]; // map for weight matrix - - // Architecture information - fArch = new int[fNLayers]; - for (int i=0; i 0) - { - pMap[i-1] = pSize; // startpoint of this layer - pSize+=fArch[i]*(1+fArch[i-1]); // size of this layer - } - } - - // Alloc parameter array - if (fNParameters != pSize) - throw EvaluationError("MultiLayerPerceptron::Constructor", "Error in the computation of nParameters"); - - // Alloc WeightMatrix (map for parameters) - fWeightMatrix = new real*[fNLayers - 1]; - - // Alloc activation function (output) matrix - fOutputMatrix = new real*[fNLayers]; - fOutputMatrix[0] = new real[fArch[0]+1]; - - for (int i=1; i 0) - { - pMap[i-1] = pSize; // startpoint of this layer - pSize+=fArch[i]*(1+fArch[i-1]); // size of this layer - } - } - // Alloc WeightMatrix (map for parameters) - fWeightMatrix = new real*[fNLayers - 1]; - - // Alloc activation function (output) matrix - fOutputMatrix = new real*[fNLayers]; - fOutputMatrix[0] = new real[fArch[0]+1]; - - for (int i=1; iGetRandomGausDev(1.0); -} - -/** - * @brief MultiLayerPerceptron::Compute - * @param in - * @param out - */ -void MultiLayerPerceptron::Compute(real* in,real* out) const -{ - // setup input - for (int i=0; i= fNLayers) - throw RangeError("MultiLayerPerceptron::GetNumNodeParams", "layer requested (" + std::to_string(layer) + ") is out of bounds!"); - - return fArch[layer-1] + 1; -} - -real* MultiLayerPerceptron::GetNodeParams(int const& layer, int const& node) -{ - if (layer <=0 || layer >= fNLayers) - throw RangeError("MultiLayerPerceptron::GetNodeParams","layer requested (" + std::to_string(layer) + ") is out of bounds!"); - - if (node < 0 || node >= fArch[layer] ) - throw RangeError("MultiLayerPerceptron::GetNodeParams","node requested (" + std::to_string(node) + ") is out of bounds!"); - - return &fWeightMatrix[layer-1][node]; -} - -// ******************** Single (hidden) layer network ************************** - -/** - * @brief SingleLayerPerceptron::SingleLayerPerceptron - * @details The constructor for a single (hidden) layer perceptron class - * @param arch a vector specifying the NN architecture (must be of the form 2-N-1) - * @param extra_pars number of additional parameters required (for derived classes) - */ -SingleLayerPerceptron::SingleLayerPerceptron(std::vector const& arch, unsigned int extra_pars): -Parametrisation(std::string("SingleLayerPerceptron"), perceptron_count_parameters(arch) + extra_pars), -fNHidden(arch[1]) -{ - if (arch.size() != 3) - throw EvaluationError("SingleLayerPerceptron::Constructor", "Requested architecture invalid: must have three layers" ); - if (arch[0] != 2) - throw EvaluationError("SingleLayerPerceptron::Constructor", "Requested architecture invalid: it must have two nodes in the first layer" ); - if (arch[2] != 1) - throw EvaluationError("SingleLayerPerceptron::Constructor", "Requested architecture invalid: it must have a single node in the third (and final) layer" ); - - InitParameters(); -} - -/** - * @brief SingleLayerPerceptron::InitParameters - * @details Initialises the parameters of the neural network, normally distributed. - */ -void SingleLayerPerceptron::InitParameters() -{ - for(int i=0; iGetRandomGausDev(1.0); -} - -/** - * @copydoc MultiLayerPerceptron::Duplicate() - */ -Parametrisation* SingleLayerPerceptron::Duplicate() -{ - const std::vector arch = {2, fNHidden, 1}; - return new SingleLayerPerceptron(arch); -} - -/** - * @brief Single layer neural network evaluation - * @details Computes the output of the neural network for a given input - * @param in, a 2-element array containing {x,log(x)} - * @param out, a pointer to the location where the output should be written - */ -void SingleLayerPerceptron::Compute(real* in,real* out) const -{ - out[0] = 0; - // Note no explicit use of fNParameters here, as it includes 'extra_pars' - const real* final_node = fParameters + 3*fNHidden; - for (int i=0; i arch = {2, fNHidden, 1}; - return new SingleLayerPerceptronPreproc(arch); -} - -/** - * @brief Initialise the parameters of the preprocessed single layer perceptron - * Values are initialised as uniform between 0 and 2/5 for the low/high-x exponents - * respectively. The lower limit is set for integrability, the higher according to - * maximum typically observed effective exponents. - * The final two parameters are used for preprocessing exponents - * small-x: fNParameters - 1 - * large-x: fNParameters - 2 - */ -void SingleLayerPerceptronPreproc::InitParameters() -{ - SingleLayerPerceptron::InitParameters(); - fParameters[fNParameters - 2] = RandomGenerator::GetRNG()->GetRandomUniform(0, 2.0)/fScaleFac; - fParameters[fNParameters - 1] = RandomGenerator::GetRNG()->GetRandomUniform(0, 5.0)/fScaleFac; -} - -/** - * @brief Initialise the parameters of the preprocessed single layer perceptron - * Preprocessing is performed on the output of SingleLayerPerceptron::Compute - */ -void SingleLayerPerceptronPreproc::Compute(real* in,real* out) const -{ - // Implemented here as alpha, beta > 0. To set lower limits use normal preproc ranges. - // Use last two parameters as preprocessing exponents - SingleLayerPerceptron::Compute(in,out); - out[0] *= pow(in[0], fScaleFac*fabs(fParameters[fNParameters -2])); - out[0] *= pow(1.0-in[0], fScaleFac*fabs(fParameters[fNParameters -1])); -} diff --git a/libnnpdf/src/pathlib.cc b/libnnpdf/src/pathlib.cc deleted file mode 100644 index 552cdfb43e..0000000000 --- a/libnnpdf/src/pathlib.cc +++ /dev/null @@ -1,60 +0,0 @@ -#include "NNPDF/exceptions.h" -#include "NNPDF/pathlib.h" -#include -#include - -namespace NNPDF -{ - -namespace -{ - -const YAML::Node &load() -{ - auto p = get_profile_path(); - try { - const static auto res = YAML::LoadFile(p); - return res; - } catch (YAML::Exception &e) { - throw FileError("pathlib", "Could not process profile file '" + p - + "': " + e.what()); - } -} - -template T get_key(std::string param) -{ - auto node = load()[param]; - // yaml-cpp error messages when a key is not found are - // pretty horrible, so we write our own. - if (!node) { - throw RuntimeException( - "pathlib", "'" + param + "' " + "key not found in profile file '" - + get_profile_path() + "'"); - } - return node.as(); -} - -} // namespace - -std::string get_profile_path() { - char *env_path = std::getenv("NNPDF_PROFILE_PATH"); - if (env_path) { - return env_path; - } - // This is defined at compile time - // It has to be volatile because tools may want - // to substitute the path in the binary, and - // the compiler needs to actually compute e.g. the - // length each time. - volatile auto p = DEFAULT_NNPDF_PROFILE_PATH; - return p; -} - -std::string get_data_path() { return get_key("data_path"); } - -std::string get_results_path() { return get_key("results_path"); } - -std::string get_config_path() { return get_key("config_path"); } - -} // namespace NNPDF - diff --git a/libnnpdf/src/pdfset.cc b/libnnpdf/src/pdfset.cc deleted file mode 100644 index bc75841a89..0000000000 --- a/libnnpdf/src/pdfset.cc +++ /dev/null @@ -1,203 +0,0 @@ -// $Id: pdfset.cc 3177 2015-08-18 14:43:31Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include - -#include "NNPDF/pdfset.h" - -#include - -namespace NNPDF -{ - - // Default verbosity level - bool PDFSet::Verbose = true; - - /** - * The constructor - */ - PDFSet::PDFSet(std::string const& pdfname, int const& members, erType const& etype): - fSetName(pdfname), - fMembers(members), - fEtype(etype) - { - get_logger() << "PDF: " << pdfname<<" ErrorType: "<< errString(etype) << " booked"<pdf->GetPDF(x, params->Q2, params->mem, params->fl)/x; - } - - static double int_xf (double x, void * p) { - struct int_param * params - = (struct int_param *)p; - return params->pdf->GetPDF(x, params->Q2, params->mem, params->fl); - } - - real PDFSet::GetPDF(real const& x, real const& Q2, int const& n, int const& fl) const - { - real* EVLN = new real[14]; - GetPDF(x,Q2,n,EVLN); - - const real retVal = EVLN[fl]; - delete[] EVLN; - - return retVal; - }; - - real PDFSet::IntegratePDF( int const& mem, int const&fl, real const& Q2, - intType xfx, bool& gslerror, gsl_integration_workspace *gsl, - real xmin, real xmax) const - { - double int_res = 0; - double int_err = 0; - - // gsl parameters - struct int_param gslparam = { this, mem, Q2, fl }; - - // GSL function - gsl_function F; - F.function = xfx ? &int_xf:&int_f; - - F.params = &gslparam; - - // Integration - int status = gsl_integration_qags (&F, xmin, xmax, 0, 1E-4, 10000, gsl, &int_res, &int_err); - - if (status == GSL_EDIVERGE || status == GSL_ESING || status == GSL_EROUND) - gslerror = true; - - return int_res; - } - - /** - * Rotate flavour basis PDFs to evolution basis - * \param LHA the les houches 13 pdfs - * \return evln the lha in the evln basis - */ - void PDFSet::LHA2EVLN(const real *LHA, real *EVLN) - { - const real uplus = LHA[U] + LHA[UBAR]; - const real uminus = LHA[U] - LHA[UBAR]; - - const real dplus = LHA[D] + LHA[DBAR]; - const real dminus = LHA[D] - LHA[DBAR]; - - const real cplus = LHA[C] + LHA[CBAR]; - const real cminus = LHA[C] - LHA[CBAR]; - - const real splus = LHA[S] + LHA[SBAR]; - const real sminus = LHA[S] - LHA[SBAR]; - - const real tplus = LHA[T] + LHA[TBAR]; - const real tminus = LHA[T] - LHA[TBAR]; - - const real bplus = LHA[B] + LHA[BBAR]; - const real bminus = LHA[B] - LHA[BBAR]; - - EVLN[0]= LHA[PHT]; // photon - EVLN[1]=(uplus + dplus + cplus + splus + tplus + bplus); //Singlet - EVLN[2]=(LHA[GLUON]); // Gluon - - EVLN[3]=( uminus + dminus + sminus + cminus + bminus + tminus ); //V - EVLN[4]=( uminus - dminus ); // V3 - EVLN[5]=( uminus + dminus - 2*sminus); // V8 - EVLN[6]=( uminus + dminus + sminus - 3*cminus); //V15 - EVLN[7]=( uminus + dminus + sminus + cminus - 4*bminus ); //V24 - EVLN[8]=( uminus + dminus + sminus + cminus + bminus - 5*tminus); // V35 - - EVLN[9]=( uplus - dplus ); // T3 - EVLN[10]=( uplus + dplus - 2*splus ); // T8 - EVLN[11]=( uplus + dplus + splus - 3*cplus ); //T15 - EVLN[12]=( uplus + dplus + splus + cplus - 4*bplus ); //T24 - EVLN[13]=( uplus + dplus + splus + cplus + bplus - 5*tplus ); // T35 - - } - - void PDFSet::EVLN2LHA(const real* EVL, real* LHA) - { - // Basis {"PHT","SNG","GLU","VAL","V03","V08","V15","V24","V35","T03","T08","T15","T24","T35"}; - - LHA[PHT] = EVL[0]; - - LHA[GLUON] = EVL[2]; - - LHA[U] = ( 10*EVL[1] - + 30*EVL[9] + 10*EVL[10] + 5*EVL[11] + 3*EVL[12] + 2*EVL[13] - + 10*EVL[3] + 30*EVL[4] + 10*EVL[5] + 5*EVL[6] + 3*EVL[7] + 2*EVL[8] ) - / 120; - - LHA[UBAR] = ( 10*EVL[1] - + 30*EVL[9] + 10*EVL[10] + 5*EVL[11] + 3*EVL[12] + 2*EVL[13] - - 10*EVL[3] - 30*EVL[4] - 10*EVL[5] - 5*EVL[6] - 3*EVL[7] - 2*EVL[8] ) - / 120; - - - LHA[D] = ( 10*EVL[1] - - 30*EVL[9] + 10*EVL[10] + 5*EVL[11] + 3*EVL[12] + 2*EVL[13] - + 10*EVL[3] - 30*EVL[4] + 10*EVL[5] + 5*EVL[6] + 3*EVL[7] + 2*EVL[8] ) - / 120; - - LHA[DBAR] = ( 10*EVL[1] - - 30*EVL[9] + 10*EVL[10] + 5*EVL[11] + 3*EVL[12] + 2*EVL[13] - - 10*EVL[3] + 30*EVL[4] - 10*EVL[5] - 5*EVL[6] - 3*EVL[7] - 2*EVL[8] ) - / 120; - - LHA[S] = ( 10*EVL[1] - - 20*EVL[10] + 5*EVL[11] + 3*EVL[12] + 2*EVL[13] - + 10*EVL[3] - 20*EVL[5] + 5*EVL[6] + 3*EVL[7] + 2*EVL[8] ) - / 120; - - LHA[SBAR] = ( 10*EVL[1] - - 20*EVL[10] + 5*EVL[11] + 3*EVL[12] + 2*EVL[13] - - 10*EVL[3] + 20*EVL[5] - 5*EVL[6] - 3*EVL[7] - 2*EVL[8] ) - / 120; - - LHA[C] = ( 10*EVL[1] - - 15*EVL[11] + 3*EVL[12] + 2*EVL[13] - + 10*EVL[3] - 15*EVL[6] + 3*EVL[7] + 2*EVL[8] ) - / 120; - - LHA[CBAR] = ( 10*EVL[1] - - 15*EVL[11] + 3*EVL[12] + 2*EVL[13] - - 10*EVL[3] + 15*EVL[6] - 3*EVL[7] - 2*EVL[8] ) - / 120; - - LHA[B] = ( 5*EVL[1] - - 6*EVL[12] + EVL[13] - + 5*EVL[3] - 6*EVL[7] + EVL[8] ) - / 60; - - LHA[BBAR] = ( 5*EVL[1] - - 6*EVL[12] + EVL[13] - - 5*EVL[3] + 6*EVL[7] - EVL[8] ) - / 60; - - LHA[T] = ( EVL[1] - - EVL[13] - + EVL[3] - EVL[8] ) - / 12; - - LHA[TBAR] = ( EVL[1] - - EVL[13] - - EVL[3] + EVL[8] ) - / 12; - } -} diff --git a/libnnpdf/src/positivity.cc b/libnnpdf/src/positivity.cc deleted file mode 100644 index 086add0520..0000000000 --- a/libnnpdf/src/positivity.cc +++ /dev/null @@ -1,182 +0,0 @@ -// $Id: positivity.cc 3187 2015-08-23 11:08:42Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include "NNPDF/positivity.h" - -#include "NNPDF/utils.h" -#include "NNPDF/fastkernel.h" -#include "NNPDF/pdfset.h" -#include "NNPDF/thpredictions.h" -#include "NNPDF/exceptions.h" - -#include -#include -#include -#include - -namespace NNPDF -{ - -/** - * @brief Positivity Constraint Set - * @param The base commondata structure - * @param The FKTable corresponding to the commondata - * @param The Lagrange multiplier - */ -PositivitySet::PositivitySet(CommonData const& data, FKTable const& table, real const& lambda): -CommonData(data), -FKTable(table), -fLambda(lambda) -{ - get_logger() << "Positivity Lagrange Multiplier: "<GetMembers() > 1) - throw LengthError("Positivity::SetBounds","bound PDF contains more than 1 replica"); - - const int Ndat = CommonData::GetNData(); - real *tmp = new real[Ndat]; - ComputePoints(pdf,tmp); - - fBounds.clear(); - for (int i = 0; i < Ndat; i++) - fBounds.push_back(-0.25*fabs(tmp[i])); - - delete[] tmp; -} - -/** - * @brief Compute contribution to error function from Positivity violation - * @param pdf The input PDF set which will be used to compute the error function - * @param res The error function for each pdf member (res[i]) - * if the theoretical prediction for such observable is < 0 we penalize the error function - */ -void PositivitySet::ComputeErf(const NNPDF::PDFSet* pdf, real* res) const -{ - const int Npdf = pdf->GetMembers(); - const int Ndat = CommonData::GetNData(); - - // Compute positivity points - real* tmp = new real[Npdf*Ndat]; - ComputePoints(pdf,tmp); - - // Contribution to Error Function - for (int i = 0; i< Ndat; i++) - for (int j = 0; j< Npdf; j++) - if (tmp[i*Npdf+j] < 0) - res[j] -= fLambda*tmp[i*Npdf+j]; - - delete[] tmp; -} - -/** - * @brief Computes the theoretical predictions for Positivity observables - * @param pdf The PDF set used in the convolution - * @param res The observable, in vector of members dimention res[i] - */ -void PositivitySet::ComputePoints(const PDFSet* pdf, real* res) const -{ - for (int i=0; iGetMembers(); i++) - res[i] = 0; - - ThPredictions::Convolute(pdf,this,res); -} - -void PositivitySet::GetPredictions(const PDFSet &pdf, real **result, int *ndata, int *npdf) -{ - - *npdf = pdf.GetMembers(); - *ndata = CommonData::GetNData(); - auto space = new real[(*npdf)*(*ndata)]; - ComputePoints(&pdf,space); - *result = space; -} - -/** - * @brief Computes the total number of violated data points (negative observables) - * @param pdf The PDF set used to compute the theoretical predictions - * @param res The total number of points which violates the positivity observable per PDF member - */ -void PositivitySet::ComputeNViolated( const PDFSet* pdf, int* res) const -{ - const int Npdf = pdf->GetMembers(); - const int Ndat = CommonData::GetNData(); - - // Compute positivity points - real* tmp = new real[Npdf*Ndat]; - ComputePoints(pdf,tmp); - - // Contribution to Error Function - for (int j = 0; j< pdf->GetMembers(); j++) - { - res[j] = 0; - for (int i = 0; i< Ndat; i++) - if (tmp[i*Npdf+j] < 0) - res[j] += 1; - } - - delete[] tmp; -} - - - -/** - * @brief Computes the total number of unacceptable data points (observables less than bounds). Ignores first and last points. - * @param pdf The PDF set used to compute the theoretical predictions - * @param pdf The PDF set used to compute the theoretical predictions - * @param res The total number of points which are unacceptable per PDF member - */ -void PositivitySet::ComputeNUnacceptable( const PDFSet* pdf, int* res) const -{ - const int Npdf = pdf->GetMembers(); - const int Ndat = CommonData::GetNData(); - - // Compute positivity points - real* tmp = new real[Npdf*Ndat]; - ComputePoints(pdf,tmp); - - // Contribution to Error Function - for (int j = 0; j< pdf->GetMembers(); j++) - { - res[j] = 0; - for (int i = 0; i< Ndat; i++) - if (tmp[i*Npdf+j] < fBounds[i]) - res[j] += 1; - } - - delete[] tmp; -} - -} diff --git a/libnnpdf/src/randomgenerator.cc b/libnnpdf/src/randomgenerator.cc deleted file mode 100644 index 4b5c8a392c..0000000000 --- a/libnnpdf/src/randomgenerator.cc +++ /dev/null @@ -1,150 +0,0 @@ -// $Id: randomgenerator.cc 2120 2014-12-02 17:14:52Z s1044006 $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include -#include "NNPDF/randomgenerator.h" -#include "NNPDF/exceptions.h" - -namespace NNPDF -{ - // Singleton - static RandomGenerator* rngInstance; - - // See GSL documentation - const gsl_rng_type* rngs[] = { gsl_rng_ranlux, - gsl_rng_cmrg, - gsl_rng_mrg, - gsl_rng_mt19937, - gsl_rng_gfsr4, - gsl_rng_ran0, - gsl_rng_ran1, - gsl_rng_ran2, - gsl_rng_ran3, - gsl_rng_rand, - gsl_rng_ranlxd1, - gsl_rng_ranlxd2, - gsl_rng_ranlxs0, - gsl_rng_ranlxs1, - gsl_rng_ranlxs2, - gsl_rng_taus, - gsl_rng_taus2 - }; - - // Default Verbosity - bool RandomGenerator::Verbose = true; - - - // Management methods - void RandomGenerator::InitRNG(int method,unsigned long int seed) - { - if (rngInstance != 0) - delete rngInstance; - - rngInstance = new RandomGenerator(method, seed); - } - - RandomGenerator* RandomGenerator::GetRNG() - { - if (rngInstance == 0) - throw InitError("GetRNG()","RNG has not been initialised!"); - - return rngInstance; - } - - /** - * @brief RandomGenerator::RandomGenerator - */ - RandomGenerator::RandomGenerator(int method, unsigned long int seed) - { - // just export GSL_RNG_TYPE=mrg to change the generator - //gsl_rng_env_setup(); - const gsl_rng_type *T = rngs[method]; - - fR = gsl_rng_alloc(T); - if (Verbose) - std::cout << "- Random Generator allocated: " << gsl_rng_name(fR) << std::endl; - gsl_rng_set(fR, seed); - } - - /** - * @brief RandomGenerator::~RandomGenerator - */ - RandomGenerator::~RandomGenerator() - { - gsl_rng_free(fR); - } - - /** - * @brief RandomGenerator::SetSeed - * @param s - */ - void RandomGenerator::SetSeed(unsigned long int s) - { - gsl_rng_set(fR, s); - } - - /** - * @brief RandomGenerator::GetRandomInt - * @return a random integer from the generator. - * Min max are equally likely and depends on the - * algorithm used. - */ - unsigned long int RandomGenerator::GetRandomInt() - { - return gsl_rng_get(fR); - } - - /** - * @brief RandomGenerator::GetRandomUniform - * @param n - * @return an integer from 0 to n-1 - */ - unsigned long int RandomGenerator::GetRandomUniform(unsigned long int n) - { - return gsl_rng_uniform_int(fR, n); - } - - /** - * @brief RandomGenerator::GetRandomUniform - * @return a double from [0,1) - */ - double RandomGenerator::GetRandomUniform() - { - return gsl_rng_uniform(fR); - } - - /** - * @brief RandomGenerator::GetRandomUniform - * @return a double from [a,b) - */ - double RandomGenerator::GetRandomUniform(double a, double b) - { - return (b-a)*gsl_rng_uniform(fR) + a; - } - - /** - * @brief RandomGenerator::GetRandomUniformPos - * @return a double from (0,1) - */ - double RandomGenerator::GetRandomUniformPos() - { - return gsl_rng_uniform_pos(fR); - } - - /** - * @brief Gaussian Random number generator - * using Box-Muller algorithm - * @param x - * @return - */ - double RandomGenerator::GetRandomGausDev(const double sigma) - { - return gsl_ran_gaussian(fR, sigma); - } - -} diff --git a/libnnpdf/src/thpredictions.cc b/libnnpdf/src/thpredictions.cc deleted file mode 100644 index 62d0eb2782..0000000000 --- a/libnnpdf/src/thpredictions.cc +++ /dev/null @@ -1,717 +0,0 @@ -// $Id: thpredictions.cc 3184 2015-08-21 19:54:00Z stefano.carrazza@mi.infn.it $ -// -// NNPDF++ 2012 -// -// Authors: Nathan Hartland, n.p.hartland@ed.ac.uk -// Stefano Carrazza, stefano.carrazza@mi.infn.it -// Luigi Del Debbio, luigi.del.debbio@ed.ac.uk - -#include -#include -#include -#include -#include -#include -#include -#include - -#include "NNPDF/common.h" -#include "NNPDF/experiments.h" -#include "NNPDF/thpredictions.h" -#include "NNPDF/dataset.h" -#include "NNPDF/fkset.h" -#include "NNPDF/fastkernel.h" -#include "NNPDF/pdfset.h" -#include "NNPDF/utils.h" -#include "NNPDF/timer.h" -#include "NNPDF/nnmpi.h" -#include "NNPDF/exceptions.h" -using namespace NNPDF; - -// Default verbosity level -bool ThPredictions::Verbose = true; - -#ifdef SSE_CONV - #include - - static inline void convolute(const real* __restrict__ x, const real* __restrict__ y, real& retval, int const& n) - { - __m128 acc = _mm_setzero_ps(); - - const __m128* a = (const __m128*) x; - const __m128* b = (const __m128*) y; - - for (int i=0; iGetMembers()), -fNData(exp->GetNData()), -fPDFName(pdfset->GetSetName()), -fSetName(exp->GetExpName()), -fEtype(pdfset->GetEtype()) -{ - // New theory array - fObs = new real[fNData*fNpdf]; - - Timer timer=Timer(); - timer.start(); - - int index = 0; - for (int i=0; iGetNSet(); i++) - { - ThPredictions::Convolute(pdfset,&(exp->GetSet(i)),fObs+index); - index+=pdfset->GetMembers()*exp->GetSet(i).GetNData(); - } - - fTconv=timer.stop(); -} - -/** - * Constructor for observables only - */ -ThPredictions::ThPredictions(const PDFSet *pdfset, const FKTable *fktab): -fObs(NULL), -fNpdf(pdfset->GetMembers()), -fNData(fktab->GetNData()), -fPDFName(pdfset->GetSetName()), -fSetName(fktab->GetDataName()), -fEtype(pdfset->GetEtype()) -{ - // New theory array - fObs = new real[fNData*fNpdf]; - - Timer timer=Timer(); - timer.start(); - - // Calculate observables - ThPredictions::Convolute(pdfset,fktab,fObs); - - fTconv=timer.stop(); -} - -/** - * Constructor for observables only - */ -ThPredictions::ThPredictions(const PDFSet *pdfset, const FKSet *fkset): -fObs(NULL), -fNpdf(pdfset->GetMembers()), -fNData(fkset->GetNDataFK()), -fPDFName(pdfset->GetSetName()), -fSetName(fkset->GetDataName()), -fEtype(pdfset->GetEtype()) -{ - // New theory array - fObs = new real[fNData*fNpdf]; - - Timer timer=Timer(); - timer.start(); - - // Calculate observables - ThPredictions::Convolute(pdfset,fkset,fObs); - - fTconv=timer.stop(); -} - -/** - * Constructor for predictions with different beam PDFs - */ -ThPredictions::ThPredictions(const PDFSet *pdf1, const PDFSet *pdf2, const FKTable* fktab): -fObs(NULL), -fNpdf(pdf1->GetMembers()), -fNData(fktab->GetNData()), -fPDFName(pdf1->GetSetName()+"_"+pdf2->GetSetName()), -fSetName(fktab->GetDataName()), -fEtype(pdf1->GetEtype()) -{ - - if (pdf1->GetMembers() != pdf2->GetMembers()) - throw RangeError("ThPredictions::ThPredictions","PDFs for different beams must have the same number of members!"); - - if (pdf1->GetEtype() != pdf2->GetEtype()) - throw RangeError("ThPredictions::ThPredictions","PDFs for different beams must have the same error type!"); - - // New theory array - fObs = new real[fNData*fNpdf]; - - Timer timer=Timer(); - timer.start(); - - // Calculate observables - ThPredictions::Convolute(pdf1,pdf2,fktab,fObs); - - fTconv=timer.stop(); -} - - -/** - * Copy Constructor - */ -ThPredictions::ThPredictions(const ThPredictions& o): -fObs(new real[o.fNData*o.fNpdf]), -fTconv(o.fTconv), -fNpdf(o.fNpdf), -fNData(o.fNData), -fPDFName(o.fPDFName), -fSetName(o.fSetName), -fEtype(o.fEtype) -{ -// Copy predictions - for (int i=0; iGetMembers(); - const int NData = fk->GetNData(); - const int Dsz = fk->GetDSz(); - const int Psz = sizeof(real)*Dsz*Npdf; - - real* sigma = fk->GetSigma(); - real *pdf = 0; - int err = posix_memalign(reinterpret_cast(&pdf), 16, Psz); - std::memset(pdf,0,Psz); - if (err != 0) - throw RangeError("ThPredictions::Convolute","ThPredictions posix_memalign " + std::to_string(err)); - - // Fetch PDF array - GetNZPDF(pdfset, fk, pdf); - - // Switcher between standard sequencial (+optional openmp) convolution - // or openmpi point-to-point communication implementation -#ifndef OPENMPI - // Calculate observables - #pragma omp parallel for - for (int i = 0; i < NData; i++) - for (int n = 0; n < Npdf; n++) - { - theory[i*Npdf + n] = 0; - convolute(pdf+Dsz*n,sigma+Dsz*i,theory[i*Npdf + n],Dsz); - } -#else - int offset, chunksize; - MPI::DecomposeDomain(NData, 0, &offset, &chunksize); // compule data for master - MPI::SendJobs(NData, Npdf, Dsz, theory, pdf, sigma); // prepare slaves - for (int i = 0; i < chunksize; i++) - for (int n = 0; n < Npdf; n++) - { - theory[i*Npdf + n] = 0; - convolute(pdf+Dsz*n,sigma+Dsz*i,theory[i*Npdf + n],Dsz); - } - MPI::RecvJobs(NData, Npdf, theory); -#endif - - // Delete pdfs - free(reinterpret_cast(pdf)); - return; -} - - /** - * Static FK level convolution function - differing beams - NO CHECK ON MATCHING HERE - */ -void ThPredictions::Convolute(const PDFSet *pdf1, const PDFSet *pdf2, const FKTable *fk, real* theory) -{ - const int Npdf = pdf1->GetMembers(); - const int NData = fk->GetNData(); - const int Dsz = fk->GetDSz(); - const int Psz = sizeof(real)*Dsz*Npdf; - - real* sigma = fk->GetSigma(); - real *pdf = 0; - int err = posix_memalign(reinterpret_cast(&pdf), 16, Psz); - std::memset(pdf,0,Psz); - - if (err != 0) - throw RangeError("ThPredictions::Convolute","ThPredictions posix_memalign " + std::to_string(err)); - - // Fetch PDF array - GetNZPDF(pdf1, pdf2, fk, pdf); - - // Switcher between standard sequencial (+optional openmp) convolution - // or openmpi point-to-point communication implementation -#ifndef OPENMPI - // Calculate observables - #pragma omp parallel for - for (int i = 0; i < NData; i++) - for (int n = 0; n < Npdf; n++) - { - theory[i*Npdf + n] = 0; - convolute(pdf+Dsz*n,sigma+Dsz*i,theory[i*Npdf + n],Dsz); - } -#else - int offset, chunksize; - MPI::ResetWorld(); - MPI::DecomposeDomain(NData, 0, &offset, &chunksize); // compule data for master - MPI::SendJobs(NData, Npdf, Dsz, theory, pdf, sigma); // prepare slaves - for (int i = 0; i < chunksize; i++) - for (int n = 0; n < Npdf; n++) - { - theory[i*Npdf + n] = 0; - convolute(pdf+Dsz*n,sigma+Dsz*i,theory[i*Npdf + n],Dsz); - } - MPI::RecvJobs(NData, Npdf, theory); -#endif - - // Delete pdfs - free(reinterpret_cast(pdf)); - return; -} - -/** - * Static convolution function - required for a fast fit - */ -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 - { - 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; iGetMembers(); - - const double* xgrid = fk->GetXGrid(); - const int* NZFlmap = fk->GetFlmap(); - const double Q20 = fk->GetQ20(); - - const int Nfl = 14; - const int NonZero = fk->GetNonZero(); - const int Nx = fk->GetNx(); - const int Tx = fk->GetTx(); - const int DSz = fk->GetDSz(); - - if ( fk->IsHadronic() ) - { - // Hadronic process - real* EVLN = new real[Nx*Nfl]; - - for (int n = 0; n < NPDF; n++) - { - for (int i = 0; i < Nx; i++) - pdfset->GetPDF(xgrid[i],Q20,n,&EVLN[i*Nfl]); - - for (int fl=0; flGetPDF(xgrid[i],Q20,n,EVLN); - for (int fl=0; flIsHadronic()) - throw EvaluationError("ThPredictions::GetNZPDF","Dual-beam convolution requires an hadronic FK table!"); - - // Get PDF set data - const int NPDF = pdf1->GetMembers(); - - const double* xgrid = fk->GetXGrid(); - const int* NZFlmap = fk->GetFlmap(); - const double Q20 = fk->GetQ20(); - - const int Nfl = 14; - const int NonZero = fk->GetNonZero(); - const int Nx = fk->GetNx(); - const int DSz = fk->GetDSz(); - - // Hadronic process - real* EVLN1 = new real[Nx*Nfl]; - real* EVLN2 = new real[Nx*Nfl]; - - for (int n = 0; n < NPDF; n++) - { - for (int i = 0; i < Nx; i++) - { - pdf1->GetPDF(xgrid[i],Q20,n,&EVLN1[i*Nfl]); - pdf2->GetPDF(xgrid[i],Q20,n,&EVLN2[i*Nfl]); - } - - for (int fl=0; fl -#include -#include -#include -#include -#include -#include -#include -#include - -#include "NNPDF/utils.h" -#include "NNPDF/exceptions.h" - -#include "gsl/gsl_matrix.h" -#include "gsl/gsl_linalg.h" - -#include -#include - -namespace NNPDF -{ -std::string joinpath(const std::initializer_list &list) -{ - return joinpath_inner(list); -} - //__________________________________________________________________ - class archive_wrapper_read - { - public: - archive_wrapper_read(): a(archive_read_new()) {} - operator struct archive*() {return a;} - ~archive_wrapper_read() - { - archive_read_free(a); - } - private: - struct archive * a; - }; - - //__________________________________________________________________ - class archive_wrapper_write - { - public: - archive_wrapper_write(): a(archive_write_new()) {} - operator struct archive*() {return a;} - ~archive_wrapper_write() - { - archive_write_finish_entry(a); - archive_write_free(a); - } - private: - struct archive * a; - }; - - //__________________________________________________________________ - class archive_entry_wrapper - { - public: - archive_entry_wrapper(): entry(archive_entry_new()) {} - operator struct archive_entry*() {return entry;} - ~archive_entry_wrapper() - { - archive_entry_free(entry); - } - - private: - struct archive_entry * entry; - }; - - //__________________________________________________________________ - void write_to_file(std::string const& filename, std::string const& data) - { - int file = creat(filename.c_str(), S_IRUSR | S_IWUSR); - if (file == -1) - throw FileError("write_to_file::creat", "Error creating file " + - filename + ". Errno = " + std::string(strerror(errno))); - - int wr_err = write(file, data.c_str(), data.size()); - if (wr_err == -1) - throw FileError("write_to_file::write", "Error writing to file " + - filename + ". Errno = " + std::string(strerror(errno))); - - int fs_err = fsync(file); - if (fs_err == -1) - throw FileError("write_to_file::write", "Error flushing file " + - filename + ". Errno = " + std::string(strerror(errno))); - - int close_err = close(file); - if (close_err == -1) - throw FileError("write_to_file::write", "Error closing file " + - filename + ". Errno = " + std::string(strerror(errno))); - } - - //__________________________________________________________________ - std::string buf_from_file(std::string const & filename){ - std::string buf{}; - std::ifstream is(filename.c_str()); - if (is.fail()){ - throw RuntimeException("utils", "Could not open " + filename); - } - - is.seekg(0, std::ios_base::end); - const auto fileSize = static_cast(is.tellg()); - buf.resize(fileSize); - - is.seekg(0, std::ios_base::beg); - is.read(&buf[0], fileSize); - - is.close(); - return buf; - } - - std::string untargz(std::string const& filename) - { - auto a = archive_wrapper_read{}; - - //The lifetime of this is managed the archive struct - struct archive_entry *entry; - std::string buf{}; - - auto filterr = archive_read_support_filter_gzip(a); - - if (!(filterr == ARCHIVE_OK || filterr == ARCHIVE_WARN)) - { - throw RuntimeException("untargz", "Failed to setup support for gzip format"); - } - - auto readerr = archive_read_support_format_tar(a); - if (readerr != ARCHIVE_OK){ - throw RuntimeException("untargz", "Failed to setup support for tar format"); - } - - // if these operations fail, most likely you have a plain txt file. - // Apparently libarchive, when supporting all the formats can, at random, either fail to - // recognize the format or fail to open the header. In fact, this is different for different - // FKTables. - // Could not understand the format - if (archive_read_open_filename(a, filename.c_str(), 10240) != ARCHIVE_OK){ - buf = buf_from_file(filename); - // Could not get the header - } else if (archive_read_next_header(a, &entry) != ARCHIVE_OK) - { - buf = buf_from_file(filename); - } - else - { - // get the entry size - auto entry_size = archive_entry_size(entry); - if (entry_size == 0) - throw RuntimeException("untargz", "Could not obtain the uncompressed size from the header in " + filename); - - // read buffer - buf.resize(entry_size); - char throwaway; - if ((archive_read_data(a, &buf[0], entry_size) != entry_size) || - (archive_read_data(a, &throwaway, 1) != 0)) - throw RuntimeException("untargz", archive_error_string(a)); - } - - return buf; - } - - //____________________________________________________________________ - void targz(std::string const& filename, std::string const& data) - { - auto a = archive_wrapper_write{}; - if (a == NULL) - throw RuntimeException("targz", "Empty archive write."); - - // Allocate the compressors and file - if ((archive_write_add_filter_gzip(a) != ARCHIVE_OK) || - (archive_write_set_format_ustar(a) != ARCHIVE_OK) || - (archive_write_open_filename(a, filename.c_str()) != ARCHIVE_OK)) - throw RuntimeException("targz", "Cannot allocate compressed file " + filename); - - // Prepare entry - auto entry = archive_entry_wrapper{}; - if (entry == NULL) - throw RuntimeException("targz", "Empty archive entry"); - - // Some options - archive_entry_set_pathname(entry, filename.c_str()); - archive_entry_set_perm(entry, 0644); - archive_entry_set_filetype(entry, AE_IFREG); - archive_entry_set_size(entry, data.size()); - - // Write header and data - if (archive_write_header(a, entry) != ARCHIVE_OK) - throw RuntimeException("targz", "Cannot write header."); - - if (archive_write_data(a, data.c_str(), data.size()) != (int) data.size()) - throw RuntimeException("targz", "Written length does not match data length"); - } - - // /very/ basic integrator - double integrate(double data[], size_t npoints, double h) - { - double integral=0; - - integral+=data[0]+data[npoints-1]; - - for ( size_t j=1; j<(npoints)/2 ; j++ ) - integral+=2*data[2*j -1]; - - for (size_t j=1; j<(npoints)/2 + 1; j++) - integral+=4*data[2*j - 2]; - - return integral*h/3.0; - } - - /** - * Split string into string vector - */ - std::vector split(std::string const &input) - { - std::stringstream strstr(input); - std::istream_iterator it(strstr); - std::istream_iterator end; - std::vector results(it, end); - return results; - } - - void split(std::vector& results, std::string const& input) - { - std::stringstream strstr(input); - std::istream_iterator it(strstr); - std::istream_iterator end; - - results.assign(it, end); - return; - } - - /** - * Split std::string into real std::vector - */ - std::vector rsplit(std::string const& input) - { - std::vector results; - char *buffer = new char[input.size() + 1]; - sprintf(buffer, "%s", input.c_str()); - char *token = strtok(buffer, " \t"); - while (token) - { - results.push_back(atof(token)); - token = strtok(NULL, " \t"); - } - delete[] buffer; - return results; - } - - void rsplit(std::vector& results, std::string const& input) - { - results.clear(); - char *buffer = new char[input.size() + 1]; - sprintf(buffer, "%s", input.c_str()); - char *token = strtok(buffer, " \t"); - while (token) - { - results.push_back(atof(token)); - token = strtok(NULL, " \t"); - } - delete[] buffer; - return; - } - - /** - * Split std::string into integer std::vector - */ - - std::vector isplit(std::string const& input) - { - std::stringstream strstr(input); - std::istream_iterator it(strstr); - std::istream_iterator end; - std::vector results(it, end); - return results; - } - - void isplit(std::vector& results, std::string const& input) - { - std::stringstream strstr(input); - std::istream_iterator it(strstr); - std::istream_iterator end; - - results.assign(it, end); - return; - } - - /** - * Compute average - * \param n number of points - * \param x array with values - * \return the average as real - */ - real ComputeAVG(int const& n, const real *x) - { - real sum = 0.0; - for (int i = 0; i < n; i++) - { - sum += x[i]; - } - - return sum / n; - } - - /** - * Compute average - * \param x std::vector with values - * \return the average as real - */ - real ComputeAVG(std::vector const& x) - { - if (x.size() != 0) - { - int n = (int) x.size(); - real sum = 0.0; - for (int i = 0; i < n; i++) - sum += x[i]; - - return sum / n; - } - - return 0; - } - - /** - * Compute the standard deviation - * \param n number of points - * \param x array with values - * \return the std dev as real - */ - real ComputeStdDev(int const& n, const real *x) - { - real sum = 0.0; - real avg = ComputeAVG(n, x); - for (int i = 0; i < n; i++) - sum += (x[i]-avg)*(x[i]-avg); - - sum /= n-1; - - return sqrt(sum); - } - - /** - * Compute the standard deviation - * \param x std::vector with values - * \return the std dev as real - */ - real ComputeStdDev(std::vector const& x) - { - if (x.size() != 0) - { - real sum = 0.0; - int n = (int) x.size(); - real avg = ComputeAVG(x); - for (int i = 0; i < n; i++) - sum += (x[i]-avg)*(x[i]-avg); - - sum /= n-1; - - return sqrt(sum); - } - - return 0; - } - - /** - * Compute the 68% c.l. - */ - void Compute68cl(std::vector const& x, real &up, real &dn) - { - up = 0; - dn = 0; - if (x.size() > 0) - { - int esc = (int) (x.size()*(1-0.68)/2); - std::vector xval(x); - std::sort(xval.begin(),xval.end()); - up = xval[xval.size()-1-esc]; - dn = xval[esc]; - } - } - - /** - * Compute the 95% c.l. - */ - void Compute95cl(std::vector const& x, real &up, real &dn) - { - up = 0; - dn = 0; - if (x.size() > 0) - { - int esc = (int) (x.size()*(1-0.95)/2); - std::vector xval(x); - std::sort(xval.begin(),xval.end()); - up = xval[xval.size()-1-esc]; - dn = xval[esc]; - } - } - - /** - * Compute the errors for the Hessian method - * \param p number of pdfs - * \param x array with values - * \return the error for the Hessian method as real - */ - real ComputeEigErr(int const& p, const real *x) - { - real err = 0; - const int nvec = (p-1)/2.0; - - for (int i = 0; i < nvec; i++) - err += pow(x[2*i+1]-x[2*i+2], 2); // Eigenstd::vector - - return sqrt(err)/2.0; - } - - /** - * Compute the errors for the symmetric Hessian method - * \param p number of pdfs - * \param x array with values - * \return the error for the Hessian method as real - */ - real ComputeSymEigErr(int const& p, const real *x) - { - real err = 0; - for (int i = 1; i < p; i++) - err += pow(x[i]-x[0], 2); // Eigenstd::vector - - return sqrt(err); - } - - - /** - * Compute the mth moment of the distribution - * \param n number of points - * \param x array with values - * \return the std dev as real - */ - real ComputeMom(int const& n, const real *x, int const& m) - { - real sum = 0.0; - real avg = ComputeAVG(n, x); - for (int i = 0; i < n; i++) - sum += pow(x[i]-avg,m); - - sum /= n; - - return sum; - } - -} diff --git a/pyproject.toml b/pyproject.toml index d64fb1216a..77ad372966 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,6 +43,8 @@ include = [ n3fit = "n3fit.scripts.n3fit_exec:main" validphys = "validphys.scripts.main:main" # Fitting scripts +evolven3fit = "n3fit.scripts.evolven3fit_new:main" +# Keep the `new` suffix so people's scripts don't break evolven3fit_new = "n3fit.scripts.evolven3fit_new:main" vp-setupfit = "n3fit.scripts.vp_setupfit:main" varflavors = "n3fit.scripts.varflavors:main"