From c16162803a3edf96d2f7a7016b50bee3c7818446 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Mon, 11 Oct 2021 11:23:27 +0200 Subject: [PATCH 01/11] update tests --- .../tests/regressions/test_datasetchi2.csv | 2 +- .../regressions/test_replicachi2data.csv | 182 +++++++++--------- 2 files changed, 92 insertions(+), 92 deletions(-) diff --git a/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv b/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv index 736dd60bbc..8d33741347 100644 --- a/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv +++ b/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv @@ -3,4 +3,4 @@ group dataset NMC NMC 204 1.6064257972017228 ATLAS ATLASTTBARTOT 3 1.9383209541765103 -CMS CMSZDIFF12 28 1.8520436886594902 +CMS CMSZDIFF12 28 1.8520436886594904 diff --git a/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv b/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv index 1f95930af2..6bd5bb4b6a 100644 --- a/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv +++ b/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv @@ -1,103 +1,103 @@ name Total NMC ATLAS CMS npoints 235 204 3 28 -0 1.6399278671240285 1.6064257972017228 1.9383209541765103 1.8520436886594902 +0 1.6399278671240285 1.6064257972017228 1.9383209541765103 1.8520436886594904 1 1.740784425773641 1.734961856769467 1.4798818462768135 1.8111598477501392 -2 1.693875430466506 1.6453399067644068 2.667768113809425 1.9431457442236328 -3 1.7059005594871608 1.681430118214625 1.6187057220806294 1.8935275070520492 +2 1.6938754304665056 1.6453399067644063 2.667768113809425 1.9431457442236326 +3 1.7059005594871604 1.6814301182146243 1.6187057220806294 1.8935275070520494 4 1.660561023557888 1.620679175623579 1.9984478347768313 1.9149266144486803 -5 1.6647552843368754 1.6454811689887392 1.489153791606265 1.823995427523003 -6 1.6605784949948608 1.643989054821457 1.5385064854912143 1.7945235601336216 -7 1.697997010626272 1.694977738806316 1.7572695047949767 1.7136439380821642 -8 1.7293825277390849 1.707333387219617 1.6845966174077818 1.8948247562021334 -9 1.6433795971020966 1.609625930205037 1.949549203659571 1.8564952837923714 -10 1.6666941197175142 1.6357991334540218 2.0934780661867998 1.8460593110869652 -11 1.6749153840441096 1.6565059983987407 2.3477312838963167 1.7369534901904902 -12 1.7249237548419944 1.7074563636898035 1.8750750574822856 1.8360985365250684 -13 1.7235703773540902 1.7113219367645702 2.0250533418572783 1.7805072697381104 -14 1.7498153896799236 1.7329714557322533 1.878439953883101 1.8587542765626104 -15 1.6471200546768532 1.6265082854552781 2.090346249256761 1.7498029953004828 -16 1.6684995117671833 1.6587447787248 1.1778912598601305 1.7921348794945862 -17 1.7120837553108796 1.7008429326771637 2.265997803972079 1.7346332435713958 -18 1.6789513582393922 1.663860589781405 1.808447334974139 1.7750238166402903 -19 1.619247012993369 1.568761881907043 2.385151759118459 1.9050060309660548 -20 1.858487298310556 1.8458849948275036 1.7220975903379312 1.9649172638270043 -21 1.6205701542992663 1.5838155488480128 1.9576559492264805 1.8522373731304838 +5 1.6647552843368754 1.6454811689887392 1.4891537916062647 1.8239954275230033 +6 1.6605784949948605 1.6439890548214566 1.5385064854912145 1.7945235601336222 +7 1.6979970106262727 1.6949777388063163 1.7572695047949767 1.7136439380821642 +8 1.7293825277390846 1.7073333872196168 1.6845966174077818 1.8948247562021334 +9 1.6433795971020968 1.6096259302050373 1.949549203659571 1.8564952837923716 +10 1.6666941197175138 1.6357991334540212 2.0934780661867998 1.8460593110869656 +11 1.6749153840441093 1.6565059983987405 2.3477312838963167 1.73695349019049 +12 1.7249237548419944 1.7074563636898035 1.8750750574822856 1.8360985365250682 +13 1.7235703773540902 1.71132193676457 2.0250533418572787 1.7805072697381108 +14 1.7498153896799242 1.7329714557322538 1.878439953883101 1.8587542765626104 +15 1.647120054676853 1.6265082854552777 2.090346249256761 1.749802995300483 +16 1.6684995117671833 1.6587447787248 1.1778912598601305 1.7921348794945864 +17 1.7120837553108799 1.700842932677164 2.265997803972079 1.734633243571396 +18 1.6789513582393922 1.6638605897814052 1.8084473349741386 1.7750238166402903 +19 1.6192470129933687 1.568761881907043 2.385151759118459 1.9050060309660544 +20 1.8584872983105563 1.8458849948275038 1.7220975903379312 1.9649172638270045 +21 1.6205701542992665 1.583815548848013 1.9576559492264805 1.8522373731304838 22 1.7068035421131555 1.7204394120543445 1.2595331411650834 1.6553783183575017 -23 1.6999475840226124 1.676171232672489 2.159292735740848 1.8239597347465573 -24 1.719875829337278 1.6869050696709849 2.2386608483338835 1.9045072548706352 -25 1.6437752789150344 1.6263906528587058 1.512641089252533 1.7844847890749826 -26 1.783585303712091 1.7700187768796092 1.887318474650766 1.871312873748172 -27 1.631086461979631 1.5943841663389928 2.0472377547374054 1.8539012631373752 -28 1.7395973964458913 1.700741982948145 2.337174973887325 1.9586606686321744 -29 1.5867778883244628 1.546183571004266 1.6087981908225046 1.8801771678182475 -30 1.6392762649343806 1.6125476885135313 1.5123662697433151 1.8476105354853265 -31 1.6143488639176584 1.5880460423569003 1.955047479173907 1.7694802836542982 -32 1.6274534000275678 1.6050499850907358 1.7273855764929598 1.7799712613746228 -33 1.781110111934827 1.779356775560228 1.5153409952320958 1.8223596823107677 -34 1.6027961117997223 1.557894017444925 2.85717562008994 1.795542137639293 -35 1.6702341135354954 1.6377982295016986 1.559850677498459 1.9183794939285554 -36 1.6088178913819464 1.5671384218201183 1.9691278134802441 1.8738779636790195 -37 1.6291577245351962 1.5990897337236325 1.7030370935132746 1.8403088680575101 -38 1.719297828925088 1.697919748678121 1.7615648026268762 1.8705238092563694 -39 1.6255435361041521 1.574392349953094 2.382778273932525 1.9170841704373944 -40 1.7106166667093075 1.6613964971953847 3.048083614153444 1.9259207287988709 -41 1.6005383418463868 1.5724534526025484 1.2877707971216625 1.8386676289862862 -42 1.626491740214993 1.5919383597818288 1.6873946532021895 1.8717124855508462 -43 1.680364529323094 1.6423587908608464 2.055360734589433 1.9170853161266461 -44 1.6841281728512845 1.6564533528283856 2.0764468583665514 1.8437248595700548 -45 1.5545388777715874 1.5186200973376842 1.6086832764775056 1.810431663928677 -46 1.644864339328371 1.5886581009650405 2.8342763328030074 1.9269299338174986 -47 1.7324999729156703 1.7202466235589389 2.0188045653411657 1.7910988833262675 -48 1.7365266333754963 1.698015685124622 2.322484435329752 1.9543252061367686 -49 1.573968516825016 1.5232530425694564 2.292018338198644 1.866533062682633 -50 1.558819007628157 1.5200417997181 2.1444737097891555 1.7785899471698932 +23 1.699947584022612 1.6761712326724887 2.159292735740848 1.8239597347465568 +24 1.7198758293372782 1.686905069670985 2.238660848333884 1.904507254870635 +25 1.6437752789150348 1.6263906528587062 1.512641089252533 1.7844847890749824 +26 1.7835853037120908 1.770018776879609 1.8873184746507663 1.8713128737481715 +27 1.6310864619796315 1.5943841663389933 2.0472377547374054 1.8539012631373752 +28 1.7395973964458915 1.7007419829481452 2.337174973887325 1.9586606686321744 +29 1.5867778883244623 1.5461835710042657 1.6087981908225046 1.8801771678182475 +30 1.6392762649343804 1.612547688513531 1.512366269743315 1.8476105354853265 +31 1.6143488639176589 1.5880460423569007 1.955047479173907 1.7694802836542982 +32 1.6274534000275676 1.6050499850907356 1.7273855764929598 1.7799712613746224 +33 1.781110111934827 1.779356775560228 1.5153409952320958 1.822359682310768 +34 1.6027961117997223 1.557894017444925 2.8571756200899405 1.7955421376392933 +35 1.670234113535496 1.6377982295016988 1.559850677498459 1.918379493928556 +36 1.6088178913819466 1.5671384218201185 1.9691278134802443 1.873877963679019 +37 1.6291577245351962 1.5990897337236325 1.7030370935132748 1.8403088680575101 +38 1.7192978289250882 1.6979197486781215 1.7615648026268762 1.8705238092563692 +39 1.6255435361041526 1.5743923499530943 2.382778273932525 1.9170841704373944 +40 1.7106166667093072 1.6613964971953845 3.048083614153444 1.9259207287988709 +41 1.6005383418463865 1.5724534526025482 1.2877707971216628 1.8386676289862862 +42 1.6264917402149925 1.5919383597818284 1.6873946532021895 1.8717124855508462 +43 1.6803645293230935 1.6423587908608461 2.055360734589433 1.9170853161266457 +44 1.6841281728512845 1.6564533528283856 2.076446858366551 1.8437248595700546 +45 1.5545388777715876 1.5186200973376842 1.6086832764775052 1.8104316639286775 +46 1.6448643393283706 1.5886581009650398 2.834276332803008 1.9269299338174986 +47 1.7324999729156703 1.7202466235589386 2.018804565341166 1.7910988833262678 +48 1.7365266333754965 1.6980156851246222 2.322484435329752 1.9543252061367684 +49 1.5739685168250162 1.5232530425694568 2.292018338198644 1.866533062682633 +50 1.558819007628157 1.5200417997181 2.1444737097891555 1.7785899471698927 51 1.6831904302201117 1.6580087457891148 2.236528645228755 1.8073707508950203 -52 1.690246908332293 1.6697587280826391 1.3798646731169555 1.8727731753528438 -53 1.507289917281574 1.4336699766955976 1.7411565364727977 2.0186066323517697 -54 1.548782514858306 1.4879833340395212 2.088119994191017 1.9339618166095185 -55 1.7027050564029675 1.6654284396198904 2.475172054195395 1.8915275146304842 -56 1.6977592496378975 1.6737669804996218 1.8758939983341882 1.8534742017135901 -57 1.6452258027297535 1.607455762288945 2.184882313797843 1.8625871854697773 -58 1.6602053713898766 1.6227941686211176 2.298092090430885 1.8644277002364409 -59 1.6810998878357646 1.6535803078778093 2.1201859498745335 1.8345547494538574 -60 1.7681944900152904 1.7578708846521391 1.6675960251545732 1.854187736039042 +52 1.6902469083322929 1.669758728082639 1.3798646731169555 1.8727731753528438 +53 1.5072899172815741 1.4336699766955983 1.7411565364727977 2.0186066323517693 +54 1.548782514858306 1.4879833340395212 2.088119994191017 1.9339618166095183 +55 1.7027050564029678 1.6654284396198906 2.475172054195395 1.8915275146304844 +56 1.6977592496378973 1.6737669804996216 1.8758939983341885 1.8534742017135901 +57 1.645225802729753 1.6074557622889447 2.184882313797843 1.8625871854697766 +58 1.6602053713898766 1.6227941686211176 2.298092090430885 1.8644277002364416 +59 1.6810998878357652 1.6535803078778097 2.1201859498745335 1.8345547494538572 +60 1.768194490015291 1.7578708846521396 1.6675960251545734 1.854187736039042 61 1.7184009723525648 1.7094673755827825 1.7357205891842886 1.7816329327290092 -62 1.6831635584337734 1.6445275547369218 1.5992944004935756 1.9736404237187144 -63 1.7723273555002863 1.764419554516597 1.4128024791751541 1.8684618565591435 -64 1.650039422766433 1.624358625245459 1.7322710171068991 1.8283318481684812 +62 1.683163558433773 1.6445275547369214 1.5992944004935756 1.9736404237187146 +63 1.7723273555002867 1.764419554516598 1.412802479175154 1.8684618565591433 +64 1.6500394227664326 1.6243586252454583 1.7322710171068991 1.8283318481684812 65 1.5330371125639766 1.4950577498604019 1.9050615277272447 1.7698841392082425 -66 1.6291577245351962 1.5990897337236325 1.7030370935132746 1.8403088680575101 -67 1.6447876636190484 1.6144363734443947 2.143387630238284 1.812497067039465 -68 1.6691288121166161 1.627326621672826 2.667058687192521 1.8667665701632428 -69 1.730533954918448 1.7108335873411127 2.1399486604643445 1.8301993431019732 -70 1.641417976963368 1.6053015596724276 1.9145294847294074 1.8752899271081434 -71 1.6100814095566427 1.5617962838621486 2.412379726003673 1.8759125057114885 -72 1.6153642244934885 1.590886202757717 1.801845331877016 1.7737239784915861 -73 1.5980786448132416 1.5497699252028407 2.312031929037349 1.8735471786650064 -74 1.6738687503676557 1.6476007887206565 1.7921104521950177 1.8525808600285747 -75 1.781524302436951 1.7741234582593268 1.5042507895587232 1.8651526149680229 +66 1.6291577245351962 1.5990897337236325 1.7030370935132748 1.8403088680575101 +67 1.6447876636190482 1.6144363734443945 2.143387630238284 1.8124970670394653 +68 1.6691288121166161 1.627326621672826 2.667058687192522 1.8667665701632428 +69 1.730533954918448 1.7108335873411127 2.1399486604643445 1.8301993431019734 +70 1.6414179769633674 1.6053015596724272 1.9145294847294074 1.8752899271081431 +71 1.6100814095566427 1.5617962838621486 2.4123797260036732 1.8759125057114885 +72 1.6153642244934883 1.5908862027577169 1.801845331877016 1.7737239784915861 +73 1.5980786448132418 1.5497699252028412 2.3120319290373486 1.8735471786650062 +74 1.6738687503676555 1.6476007887206563 1.7921104521950177 1.8525808600285745 +75 1.7815243024369507 1.7741234582593264 1.5042507895587232 1.8651526149680229 76 1.5391391144414883 1.4801487958594524 2.9456669277985914 1.8182263126794906 -77 1.6273482921436313 1.566560490353848 2.172870146295794 2.0117820779528923 -78 1.6082049689793336 1.5687010433284765 1.650374554293857 1.8915011145804497 -79 1.6252311336565095 1.598573289779018 1.7535733181683792 1.8057016192805346 +77 1.6273482921436313 1.566560490353848 2.1728701462957933 2.011782077952892 +78 1.6082049689793334 1.5687010433284763 1.650374554293857 1.8915011145804492 +79 1.625231133656509 1.5985732897790172 1.7535733181683792 1.8057016192805346 80 1.710942312448346 1.690631118932874 2.0214853202106733 1.8256514000865363 81 1.6222172783375624 1.583147892062274 1.4558364356866889 1.9246921829129704 -82 1.556610403614142 1.5048484827713324 1.8833389150166278 1.898726343532917 -83 1.8035683234886928 1.8012390234686313 1.8247001953326532 1.8182748087944318 -84 1.650848340496093 1.6262304929862843 2.0792882017701513 1.7843026729310496 -85 1.7560914868788136 1.737884446482184 2.069899333479669 1.8551205119184502 -86 1.7214822079683119 1.6841598356444387 2.519004188655497 1.9079535655400475 -87 1.7633224041039264 1.7371573611160664 1.5240617818959412 1.9795884982520489 -88 1.6099494843008484 1.5698306556608805 2.5112292010956794 1.8056781233068828 -89 1.6489349161773976 1.6367559559810332 1.8212787692799088 1.7192019276327861 +82 1.5566104036141417 1.5048484827713322 1.8833389150166273 1.898726343532917 +83 1.8035683234886932 1.8012390234686315 1.8247001953326534 1.8182748087944318 +84 1.650848340496093 1.6262304929862843 2.079288201770151 1.7843026729310498 +85 1.7560914868788136 1.7378844464821845 2.069899333479669 1.8551205119184502 +86 1.7214822079683119 1.6841598356444387 2.519004188655497 1.9079535655400477 +87 1.7633224041039264 1.7371573611160664 1.5240617818959412 1.979588498252049 +88 1.609949484300848 1.56983065566088 2.5112292010956794 1.8056781233068822 +89 1.648934916177398 1.6367559559810334 1.8212787692799093 1.719201927632786 90 1.7091519751224513 1.6985942170700976 1.9883539507359351 1.7561582864024405 -91 1.6194034580546435 1.5500116616442614 2.758220862705391 2.0029561099748467 -92 1.720444971345963 1.678772855568108 1.9085714160795701 2.0038996957917345 -93 1.8319481691933461 1.8257548306283655 1.6774250879884554 1.8936271088673013 -94 1.6395812158946916 1.5883899915574327 2.388944145078526 1.932256965082169 -95 1.6242247583222318 1.5672171736940481 2.2247853271094526 1.975219956814653 -96 1.5995602033165557 1.5635473746690751 1.5047951595056765 1.8720927810136507 -97 1.6064411114535724 1.5505054784799164 2.090761326611409 1.9620807000661546 -98 1.651155114020226 1.6140965055263958 2.5081272917568054 1.829335099717786 -99 1.6189525662792061 1.5767427834943755 1.9080945300041143 1.8955014875981586 -100 1.6527589266758624 1.6294564991751102 1.4758331094786497 1.8414900931667595 +91 1.6194034580546435 1.5500116616442616 2.75822086270539 2.0029561099748467 +92 1.720444971345963 1.678772855568108 1.9085714160795704 2.003899695791734 +93 1.8319481691933468 1.8257548306283662 1.6774250879884551 1.8936271088673016 +94 1.6395812158946916 1.5883899915574327 2.3889441450785265 1.9322569650821693 +95 1.6242247583222316 1.567217173694048 2.224785327109452 1.975219956814653 +96 1.599560203316556 1.5635473746690753 1.5047951595056765 1.8720927810136507 +97 1.6064411114535726 1.5505054784799166 2.090761326611409 1.9620807000661546 +98 1.6511551140202259 1.6140965055263956 2.508127291756806 1.8293350997177857 +99 1.6189525662792057 1.576742783494375 1.9080945300041143 1.8955014875981586 +100 1.6527589266758624 1.6294564991751102 1.4758331094786497 1.8414900931667597 From 8748a8f552f18b13f5c975ca27027e52906a7a12 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 13 Oct 2021 09:41:54 +0200 Subject: [PATCH 02/11] add a lhapdfset module --- validphys2/src/validphys/lhapdfset.py | 161 ++++++++++++++++++++++++++ 1 file changed, 161 insertions(+) create mode 100644 validphys2/src/validphys/lhapdfset.py diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py new file mode 100644 index 0000000000..a0e5ce40b7 --- /dev/null +++ b/validphys2/src/validphys/lhapdfset.py @@ -0,0 +1,161 @@ +""" + Module containing an LHAPDF class compatible with validphys + using the official lhapdf python interface. + + It exposes an interface (almost) compatible with libNNPDF::LHAPDFSet +""" +import logging +from typing import NamedTuple +import numpy as np +import lhapdf + +log = logging.getLogger(__name__) + + +class PDFErrorType(NamedTuple): + """Namedtuple containing the information about the error type of the PDF""" + + name: str + monte_carlo: bool + libNNPDF: int + t0: bool + description: str + + +ER_NONE = PDFErrorType("erType_ER_NONE", False, 0, False, "No error info") +ER_MC = PDFErrorType("erType_ER_MC", True, 1, False, "Monte Carlo") +ER_MC68 = PDFErrorType("erType_ER_MC68", True, 2, False, "Monte Carlo 68pc") +ER_MCT0 = PDFErrorType("erType_ER_MCT0", False, 3, True, "Monte Carlo t0") +ER_EIG = PDFErrorType("erType_ER_EIG", False, 4, False, "Eigenvector 68pc") +ER_EIG90 = PDFErrorType("erType_ER_EIG90", False, 5, False, "Eigenvector 90pc") +ER_SYMEIG = PDFErrorType("erType_ER_SYMEIG", False, 6, False, "Symmetric eigenvectors") +_libNNPDF_errors = [ER_NONE, ER_MC, ER_MC68, ER_MCT0, ER_EIG, ER_EIG90, ER_SYMEIG] + + +def _parse_flavs(f): + # TODO: for now we believe what the user unless they say 0, then they're wrong + if f == 0: + return 21 + return f + + +class LHAPDFSet: + """Wrapper for the lhapdf python interface""" + + def __init__(self, name, error_type, lhapdf_verbosity=0): + if isinstance(error_type, int): + # libNNPDF error types were int + error_type = _libNNPDF_errors[error_type] + self._name = name + self._error_type = error_type + self._flavors = None + self._lhapdf_set = lhapdf.mkPDFs(name) + self._libNNPDF_set = None + self.legacy_interface() + log.info("PDF: %s ErrorType: %s", name, error_type.name) + lhapdf.setVerbosity(lhapdf_verbosity) + + def legacy_interface(self): + """Setting some methods and attribute as per libNNPDF specs""" + self.GetMembers = self.get_members + self.GetSetName = lambda self: self.name + + def get_members(self): + """Get the number of members""" + return len(self.members) + + @property + def is_monte_carlo(self): + """Check whether the error type is MC""" + return self._error_type.monte_carlo + + @property + def is_t0(self): + """Check whether we are in t0 mode""" + return self._error_type.t0 + + @property + def members(self): + """Return the members of the set, this depends on the error type: + t0: returns only member 0 + MC: skip member 0 + Eigen: return all + """ + if self.is_t0: + return self._lhapdf_set[0:1] + if self.is_monte_carlo: + return self._lhapdf_set[1:] + return self._lhapdf_set + + def as_libNNPDF(self): + """Imports libNNPDF to return the PDF set in NNPDF form + needed for classes that still rely in libNNPDF and will need a PDF + (at the moment due to closure tests) + """ + if self._libNNPDF_set is None: + # TODO obvs + from NNPDF import LHAPDFSet + + self._libNNPDF_set = LHAPDFSet(self._name, self._error_type.libNNPDF) + return self._libNNPDF_set + + def xfxQ(self, x, q, member_idx, flavour): + """Return the PDF value for one single point""" + member_pdf = self.members[member_idx] + res = member_pdf.xfxQ(x, q) + return res[_parse_flavs(flavour)] + + @property + def flavors(self): + """Returns the list of accepted flavors by the LHAPDF set""" + if self._flavors is None: + self._flavors = self.members[0].flavors() + return self._flavors + + def grid_values(self, flavors, xgrid, qgrid): + """Reimplementation of libNNPDF grid_values + The return shape is + (members, flavors, xgrid, qgrid) + + Examples + -------- + >>> import numpy as np + >>> from validphys.lhapdfset import LHAPDFSet, ER_MC + >>> pdf = LHAPDFSet("NNPDF40_nnlo_as_01180", ER_MC) + >>> xgrid = np.random.rand(10) + >>> qgrid = np.random.rand(3) + >>> flavs = np.arange(-4,4) + >>> results = valid_pdf.grid_values(flavs, xgrid, qgrid) + """ + # TODO Digest the input flavors in a less confusing way + parsed_flavors = [_parse_flavs(f) for f in flavors] + flavs = [] + non_flavs = [] + for pf in parsed_flavors: + if pf in self.flavors: + flavs.append(pf) + else: + non_flavs.append(pf) + xvals = [] + for x in xgrid: + qvals = [] + for q in qgrid: + mres = [] + for member in self.members: + all_flavs_dict = member.xfxQ(x, q) + mres.append([all_flavs_dict[i] for i in flavs]) + qvals.append(mres) + xvals.append(qvals) + ret = np.array(xvals) + # Add 0s for flavours that were not in the PDF but the user asked for + # Could equally well be done in the loop above but looked like an unnecesary complication + zeros = np.zeros_like(ret[..., 0]) + for f in non_flavs: + idx = parsed_flavors.index(f) + # Do it one by one to keep the index in the right places + ret = np.insert(ret, idx, zeros, -1) + return np.transpose(ret, (2, 3, 0, 1)) + + +for error_type in _libNNPDF_errors: + setattr(LHAPDFSet, error_type.name, error_type) From aeaf1b310ee90d38fccf71e2cf8c27e677181e95 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 13 Oct 2021 09:42:50 +0200 Subject: [PATCH 03/11] use the new module lhapdfset instead of libNNPDF --- validphys2/src/validphys/core.py | 18 ++---------------- validphys2/src/validphys/covmats.py | 4 ++-- validphys2/src/validphys/filters.py | 2 +- validphys2/src/validphys/gridvalues.py | 14 ++++---------- validphys2/src/validphys/sumrules.py | 2 +- 5 files changed, 10 insertions(+), 30 deletions(-) diff --git a/validphys2/src/validphys/core.py b/validphys2/src/validphys/core.py index 261f4b4a89..14cf152025 100644 --- a/validphys2/src/validphys/core.py +++ b/validphys2/src/validphys/core.py @@ -22,7 +22,7 @@ from reportengine.baseexceptions import AsInputError from reportengine.compat import yaml -from NNPDF import (LHAPDFSet, +from NNPDF import ( CommonData, FKTable, FKSet, @@ -37,23 +37,10 @@ from validphys.theorydbutils import fetch_theory from validphys.hyperoptplot import HyperoptTrial from validphys.utils import experiments_to_dataset_inputs +from validphys.lhapdfset import LHAPDFSet log = logging.getLogger(__name__) - -#TODO: Remove this eventually -#Bacward compatibility error type names -#Swig renamed these for no reason whatsoever. -try: - LHAPDFSet.erType_ER_EIG -except AttributeError: - import warnings - warnings.warn("libnnpdf out of date. Setting backwards compatible names") - LHAPDFSet.erType_ER_MC = LHAPDFSet.ER_MC - LHAPDFSet.erType_ER_EIG = LHAPDFSet.ER_EIG - LHAPDFSet.erType_ER_EIG90 = LHAPDFSet.ER_EIG90 - LHAPDFSet.erType_ER_SYMEIG = LHAPDFSet.ER_SYMEIG - class TupleComp: @classmethod @@ -167,7 +154,6 @@ def load_t0(self): """Load the PDF as a t0 set""" return LHAPDFSet(self.name, LHAPDFSet.erType_ER_MCT0) - def __str__(self): return self.label diff --git a/validphys2/src/validphys/covmats.py b/validphys2/src/validphys/covmats.py index f66cfc5ba8..ade32e08e4 100644 --- a/validphys2/src/validphys/covmats.py +++ b/validphys2/src/validphys/covmats.py @@ -489,7 +489,7 @@ def dataset_inputs_covmat( #Copy data to avoid chaos loaded_data = type(loaded_data)(loaded_data) log.debug("Setting T0 predictions for %s" % data) - loaded_data.SetT0(t0set.load_t0()) + loaded_data.SetT0(t0set.load_t0().as_libNNPDF()) covmat = loaded_data.get_covmat() @@ -581,7 +581,7 @@ def covmat( #Copy data to avoid chaos loaded_data = type(loaded_data)(loaded_data) log.debug("Setting T0 predictions for %s" % dataset) - loaded_data.SetT0(t0set.load_t0()) + loaded_data.SetT0(t0set.load_t0().as_libNNPDF()) covmat = loaded_data.get_covmat() if fitthcovmat: diff --git a/validphys2/src/validphys/filters.py b/validphys2/src/validphys/filters.py index d364764b85..bd207a725e 100644 --- a/validphys2/src/validphys/filters.py +++ b/validphys2/src/validphys/filters.py @@ -170,7 +170,7 @@ def _filter_closure_data(filter_path, data, fakepdfset, fakenoise, errorsize): # Load data, don't cache result loaded_data = data.load.__wrapped__(data) # generate level 1 shift if fakenoise - loaded_data.MakeClosure(fakeset, fakenoise) + loaded_data.MakeClosure(fakeset.as_libNNPDF(), fakenoise) for j, dataset in enumerate(data.datasets): path = filter_path / dataset.name nfull, ncut = _write_ds_cut_data(path, dataset) diff --git a/validphys2/src/validphys/gridvalues.py b/validphys2/src/validphys/gridvalues.py index 7a6b7bcd2d..3b77f8bd3a 100644 --- a/validphys2/src/validphys/gridvalues.py +++ b/validphys2/src/validphys/gridvalues.py @@ -10,14 +10,8 @@ import numpy as np -from NNPDF import REALDOUBLE, LHAPDFSet - from validphys.core import PDF - -#Numpy is unhappy about downcasting double to float implicitly, so we have -#to manually make sure all input arrays correspond to NNPDF::real. -FLTYPE = np.int32 -REALTYPE = np.double if REALDOUBLE else np.float32 +from validphys.lhapdfset import LHAPDFSet # Canonical ordering of PDG quark flavour codes QUARK_FLAVOURS = (-6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6) @@ -54,9 +48,9 @@ def _grid_values(lpdf, flmat, xmat, qmat): """Compute lpdf.grid_values with more forgiving argument types""" - flmat = np.atleast_1d(np.asanyarray(flmat, dtype=FLTYPE)) - xmat = np.atleast_1d(np.asarray(xmat, dtype=REALTYPE)) - qmat = np.atleast_1d(np.asarray(qmat, dtype=REALTYPE)) + flmat = np.atleast_1d(np.asanyarray(flmat)) + xmat = np.atleast_1d(np.asarray(xmat)) + qmat = np.atleast_1d(np.asarray(qmat)) return lpdf.grid_values(flmat, xmat, qmat) def grid_values(pdf:PDF, flmat, xmat, qmat): diff --git a/validphys2/src/validphys/sumrules.py b/validphys2/src/validphys/sumrules.py index 7ca3774bcf..5b6db2d17b 100644 --- a/validphys2/src/validphys/sumrules.py +++ b/validphys2/src/validphys/sumrules.py @@ -13,13 +13,13 @@ import pandas as pd import scipy.integrate as integrate -from NNPDF import LHAPDFSet from reportengine.table import table from reportengine.checks import check_positive from reportengine.floatformatting import format_error_value_columns from validphys.core import PDF from validphys.pdfbases import ALL_FLAVOURS, parse_flarr +from validphys.lhapdfset import LHAPDFSet def _momentum_sum_rule_integrand(x, lpdf:LHAPDFSet, irep, Q): From 21175a61d746b2229d4e5c993f03687a30e9b48b Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 13 Oct 2021 14:33:23 +0200 Subject: [PATCH 04/11] add some documentation and simplify; --- validphys2/src/validphys/lhapdfset.py | 79 +++++++++++++++++++-------- 1 file changed, 55 insertions(+), 24 deletions(-) diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py index a0e5ce40b7..d63816cb77 100644 --- a/validphys2/src/validphys/lhapdfset.py +++ b/validphys2/src/validphys/lhapdfset.py @@ -3,6 +3,35 @@ using the official lhapdf python interface. It exposes an interface (almost) compatible with libNNPDF::LHAPDFSet + + The ``.members`` and ``.central_member`` of the ``LHAPDFSet`` are + LHAPDF objects (the typical output from ``mkPDFs``) and can be used normally. + + For MC PDFs the ``central_member`` is the average of all replicas (members 1-100) + while for Hessian PDFs the ``central_member`` is also ``member[0]`` + + Examples + -------- + >>> from validphys.lhapdfset import LHAPDFSet, ER_MC + >>> pdf = LHAPDFSet("NNPDF40_nnlo_as_01180", ER_MC) + >>> pdf.get_members() + 100 + >>> len(pdf.members) + 100 + >>> pdf.central_member.alphasQ(91.19) + 0.11800 + >>> pdf.member[0].xfxQ2(0.5, 15625) + {-5: 6.983360500601136e-05, + -4: 0.0021818063617227604, + -3: 0.00172453472243952, + -2: 0.0010906577230485718, + -1: 0.0022049272225017286, + 1: 0.020051104853608722, + 2: 0.0954139944889494, + 3: 0.004116641378803191, + 4: 0.002180124185625795, + 5: 6.922722705177504e-05, + 21: 0.007604124516892057} """ import logging from typing import NamedTuple @@ -32,13 +61,6 @@ class PDFErrorType(NamedTuple): _libNNPDF_errors = [ER_NONE, ER_MC, ER_MC68, ER_MCT0, ER_EIG, ER_EIG90, ER_SYMEIG] -def _parse_flavs(f): - # TODO: for now we believe what the user unless they say 0, then they're wrong - if f == 0: - return 21 - return f - - class LHAPDFSet: """Wrapper for the lhapdf python interface""" @@ -52,7 +74,7 @@ def __init__(self, name, error_type, lhapdf_verbosity=0): self._lhapdf_set = lhapdf.mkPDFs(name) self._libNNPDF_set = None self.legacy_interface() - log.info("PDF: %s ErrorType: %s", name, error_type.name) + log.info("PDF: %s ErrorType: %s", name, error_type.description) lhapdf.setVerbosity(lhapdf_verbosity) def legacy_interface(self): @@ -87,23 +109,28 @@ def members(self): return self._lhapdf_set[1:] return self._lhapdf_set + @property + def central_member(self): + """Returns a reference to member 0 of the PDF list""" + return self._lhapdf_set[0] + def as_libNNPDF(self): """Imports libNNPDF to return the PDF set in NNPDF form needed for classes that still rely in libNNPDF and will need a PDF (at the moment due to closure tests) """ if self._libNNPDF_set is None: - # TODO obvs + # TODO needed for C++ datatypes from NNPDF import LHAPDFSet self._libNNPDF_set = LHAPDFSet(self._name, self._error_type.libNNPDF) return self._libNNPDF_set def xfxQ(self, x, q, member_idx, flavour): - """Return the PDF value for one single point""" + """Return the PDF value for one single point for one single member""" member_pdf = self.members[member_idx] res = member_pdf.xfxQ(x, q) - return res[_parse_flavs(flavour)] + return res[flavour] @property def flavors(self): @@ -127,15 +154,18 @@ def grid_values(self, flavors, xgrid, qgrid): >>> flavs = np.arange(-4,4) >>> results = valid_pdf.grid_values(flavs, xgrid, qgrid) """ - # TODO Digest the input flavors in a less confusing way - parsed_flavors = [_parse_flavs(f) for f in flavors] - flavs = [] - non_flavs = [] - for pf in parsed_flavors: - if pf in self.flavors: - flavs.append(pf) + # LHAPDFSet.grid_values accepts flavours not included in the the PDF + # keep track of which they are to add them later as 0s + # The LHAPDFSet way is that it is accepted to ask for flavours outside the ones + # accepted by the PDF, they will just be made of 0s + flavs, zero_idxs = [], [] + for idx, f in enumerate(flavors): + if f in self.flavors: + flavs.append(f) else: - non_flavs.append(pf) + zero_idxs.append(idx) + + # Grid values loop xvals = [] for x in xgrid: qvals = [] @@ -147,15 +177,16 @@ def grid_values(self, flavors, xgrid, qgrid): qvals.append(mres) xvals.append(qvals) ret = np.array(xvals) + # Add 0s for flavours that were not in the PDF but the user asked for - # Could equally well be done in the loop above but looked like an unnecesary complication - zeros = np.zeros_like(ret[..., 0]) - for f in non_flavs: - idx = parsed_flavors.index(f) + for idx in zero_idxs: # Do it one by one to keep the index in the right places - ret = np.insert(ret, idx, zeros, -1) + ret = np.insert(ret, idx, 0.0, -1) + + # Finally return in the preferred shape return np.transpose(ret, (2, 3, 0, 1)) +# libNNPDF compatibility for error_type in _libNNPDF_errors: setattr(LHAPDFSet, error_type.name, error_type) From 85dfc9d398d239ee0efef5df2b190dce275a287a Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 14 Oct 2021 09:41:59 +0200 Subject: [PATCH 05/11] fix example --- validphys2/src/validphys/lhapdfset.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py index d63816cb77..4969adeb0f 100644 --- a/validphys2/src/validphys/lhapdfset.py +++ b/validphys2/src/validphys/lhapdfset.py @@ -152,7 +152,8 @@ def grid_values(self, flavors, xgrid, qgrid): >>> xgrid = np.random.rand(10) >>> qgrid = np.random.rand(3) >>> flavs = np.arange(-4,4) - >>> results = valid_pdf.grid_values(flavs, xgrid, qgrid) + >>> flavs[4] = 21 + >>> results = pdf.grid_values(flavs, xgrid, qgrid) """ # LHAPDFSet.grid_values accepts flavours not included in the the PDF # keep track of which they are to add them later as 0s From f701d7f9a76d974bf27daca45ba2166965e4f743 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 14 Oct 2021 12:44:04 +0200 Subject: [PATCH 06/11] add the possibility of getting also the CV --- validphys2/src/validphys/core.py | 51 ++++++++++++++++++++++++--- validphys2/src/validphys/lhapdfset.py | 36 +++++++++++++++---- validphys2/src/validphys/results.py | 36 ++++++++++--------- 3 files changed, 95 insertions(+), 28 deletions(-) diff --git a/validphys2/src/validphys/core.py b/validphys2/src/validphys/core.py index 14cf152025..566f9e8897 100644 --- a/validphys2/src/validphys/core.py +++ b/validphys2/src/validphys/core.py @@ -9,6 +9,7 @@ from __future__ import generator_stop from collections import namedtuple +from contextlib import contextmanager import re import enum import functools @@ -78,10 +79,33 @@ def __dir__(self): PDFSETS = _PDFSETS() class PDF(TupleComp): + """Wrapper class for the validphys PDF object to easily manage + both Monte Carlo and Hessian sets + + Offers a context manager (``enable_central_value``) to include the central value + in Monte Carlo sets. + + Examples + -------- + >>> from validphys.api import API + >>> from validphys.convolution import predictions + >>> args = {"dataset_input":{"dataset": "ATLASTTBARTOT"}, "theoryid":200, "use_cuts":"internal"} + >>> ds = API.dataset(**args) + >>> pdf = API.pdf(pdf="NNPDF40_nnlo_as_01180") + >>> with pdf.enable_central_value(): + >>> preds_with_cv = predictions(ds, pdf) + >>> preds_no_cv = predictions(ds, pdf) + >>> len(preds_with_cv.columns) + 101 + >>> len(preds_no_cv.columns) + 100 + """ def __init__(self, name): self.name = name self._plotname = name + self._lhapdfset = None + self._include_cv = False super().__init__(name) @@ -145,13 +169,14 @@ def rescale_factor(self): else: return 1 - @functools.lru_cache(maxsize=16) def load(self): - return LHAPDFSet(self.name, self.nnpdf_error) + if self._lhapdfset is None: + self._lhapdfset = LHAPDFSet(self.name, self.nnpdf_error) + return self._lhapdfset @functools.lru_cache(maxsize=2) def load_t0(self): - """Load the PDF as a t0 set""" + """Reload the PDF as a t0 set""" return LHAPDFSet(self.name, LHAPDFSet.erType_ER_MCT0) def __str__(self): @@ -164,7 +189,7 @@ def __len__(self): @property def nnpdf_error(self): - """Return the NNPDF error tag, used to build the `LHAPDFSet` objeect""" + """Return the NNPDF error tag, used to build the `LHAPDFSet` object""" error = self.ErrorType if error == "replicas": return LHAPDFSet.erType_ER_MC @@ -203,6 +228,8 @@ def grid_values_index(self): len(pdf))`` for Monte Carlo sets, because replica 0 is not selected and ``range(0, len(pdf))`` for hessian sets. + If ``include_cv`` is set to True, add a central value column for member 0 + for Monte Carlo error sets Returns ------- @@ -215,7 +242,10 @@ def grid_values_index(self): """ err = self.nnpdf_error if err is LHAPDFSet.erType_ER_MC: - return range(1, len(self)) + if self._include_cv: + return ["CV"] + list(range(1, len(self))) + else: + return range(1, len(self)) elif err in (LHAPDFSet.erType_ER_SYMEIG, LHAPDFSet.erType_ER_EIG, LHAPDFSet.erType_ER_EIG90): return range(0, len(self)) else: @@ -228,6 +258,17 @@ def get_members(self): """ return len(self.grid_values_index) + @contextmanager + def enable_central_value(self): + """Context manager within which the central value is included + regardless of the error type of the PDF set""" + # Get a reference to the base PDF set of this class + pdfset = self.load() + pdfset.include_cv = True + self._include_cv = True + yield + self._include_cv = False + pdfset.include_cv = False kinlabels_latex = CommonData.kinLabel_latex.asdict() diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py index 4969adeb0f..aa221e8db0 100644 --- a/validphys2/src/validphys/lhapdfset.py +++ b/validphys2/src/validphys/lhapdfset.py @@ -62,20 +62,40 @@ class PDFErrorType(NamedTuple): class LHAPDFSet: - """Wrapper for the lhapdf python interface""" + """Wrapper for the lhapdf python interface. + + Once instantiated this class will load the PDF set according to whether it is to be + treated as a T0 set (only the CV) or not. + + It is possible to control the LHAPDF verbosity with the flag ``lhapdf_verbosity``. + + For Monte Carlo sets the central value (member 0) is by default not included when taking + the resutls for all members (i.e., when using ``grid_values``). + However, it is possible to add member 0 by changing the ``include_cv`` attribute to True. + + Temporarily: it exposes all libNNPDF attributes that were exposed and used prior to + the introduction of this class + """ def __init__(self, name, error_type, lhapdf_verbosity=0): + log.info("PDF: %s ErrorType: %s", name, error_type.description) if isinstance(error_type, int): # libNNPDF error types were int error_type = _libNNPDF_errors[error_type] self._name = name self._error_type = error_type + # If at this point we already know this is a T0 set, load only the CV + if error_type.t0: + self._lhapdf_set = [lhapdf.mkPDF(name)] + else: + self._lhapdf_set = lhapdf.mkPDFs(name) self._flavors = None - self._lhapdf_set = lhapdf.mkPDFs(name) self._libNNPDF_set = None - self.legacy_interface() - log.info("PDF: %s ErrorType: %s", name, error_type.description) + self.include_cv = False + # Set the verbosity of LHAPDF lhapdf.setVerbosity(lhapdf_verbosity) + # Prepare a Legacy Interface + self.legacy_interface() def legacy_interface(self): """Setting some methods and attribute as per libNNPDF specs""" @@ -105,7 +125,7 @@ def members(self): """ if self.is_t0: return self._lhapdf_set[0:1] - if self.is_monte_carlo: + if self.is_monte_carlo and not self.include_cv: return self._lhapdf_set[1:] return self._lhapdf_set @@ -139,11 +159,15 @@ def flavors(self): self._flavors = self.members[0].flavors() return self._flavors - def grid_values(self, flavors, xgrid, qgrid): + def grid_values(self, flavors: np.ndarray, xgrid: np.ndarray, qgrid: np.ndarray): """Reimplementation of libNNPDF grid_values The return shape is (members, flavors, xgrid, qgrid) + Return + ------ + ndarray of shape (members, flavors, xgrid, qgrid) + Examples -------- >>> import numpy as np diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py index 15f6c52862..0d76628d59 100644 --- a/validphys2/src/validphys/results.py +++ b/validphys2/src/validphys/results.py @@ -48,11 +48,13 @@ class NNPDFDataResult(Result): """A result fills its values from a pandas dataframe For legacy (libNNPDF) compatibility, falls back to libNNPDF attributes""" - def __init__(self, dataobj=None, central_value=None): + def __init__(self, dataobj): # This class is used by both validphys and libNNPDF objects - # when central_value is not explictly passed, fallback to - # libNNPDF object .get_cv() - if central_value is None: + # At this point only the result of the ThPredictions is a vp object + # which includes a special CV column which needs to be pop'd + try: + central_value = dataobj.pop("CV") + except AttributeError: central_value = dataobj.get_cv() self._central_value = np.array(central_value).squeeze() @@ -78,8 +80,8 @@ def std_error(self): class DataResult(NNPDFDataResult): - def __init__(self, dataobj, covmat, sqrtcovmat, central_value=None): - super().__init__(dataobj, central_value=central_value) + def __init__(self, dataobj, covmat, sqrtcovmat): + super().__init__(dataobj) self._covmat = covmat self._sqrtcovmat = sqrtcovmat @@ -105,17 +107,19 @@ class ThPredictionsResult(NNPDFDataResult): """Class holding theory prediction For legacy purposes it still accepts libNNPDF datatypes, but prefers python-pure stuff """ - def __init__(self, dataobj, stats_class, label=None, central_value=None): + def __init__(self, dataobj, stats_class, label=None): + super().__init__(dataobj) self.stats_class = stats_class self.label = label - # Ducktype the input into numpy arrays + # Ducktype the input into numpy arrays for the rawdata + # TODO: once all of them are dataframes, the rawdata could be made of + # dataframes as well try: self._std_error = dataobj.std(axis=1).to_numpy() self._rawdata = dataobj.to_numpy() except AttributeError: self._std_error = dataobj.get_error() self._rawdata = dataobj.get_data() - super().__init__(dataobj, central_value=central_value) @property def std_error(self): @@ -145,20 +149,18 @@ def from_convolution(cls, pdf, dataset): datasets = (dataset,) try: - all_preds = [] - all_centrals = [] - for d in datasets: - all_preds.append(predictions(d, pdf)) - all_centrals.append(central_predictions(d, pdf)) + with pdf.enable_central_value(): + th_predictions = pd.concat(predictions(d, pdf) for d in datasets) except PredictionsRequireCutsError as e: raise PredictionsRequireCutsError("Predictions from FKTables always require cuts, " "if you want to use the fktable intrinsic cuts set `use_cuts: 'internal'`") from e - th_predictions = pd.concat(all_preds) - central_values = pd.concat(all_centrals) + # For Hessian sets + if "CV" not in th_predictions: + th_predictions["CV"] = th_predictions[0] label = cls.make_label(pdf, dataset) - return cls(th_predictions, pdf.stats_class, label, central_value=central_values) + return cls(th_predictions, pdf.stats_class, label) class PositivityResult(StatsResult): From 68cf5a22781e485c13f766c8c29da8de21142605 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Thu, 14 Oct 2021 17:46:57 +0200 Subject: [PATCH 07/11] add docs --- doc/sphinx/source/vp/pydataobjs.rst | 33 +++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/doc/sphinx/source/vp/pydataobjs.rst b/doc/sphinx/source/vp/pydataobjs.rst index 7164322b27..8fbae1eaf3 100644 --- a/doc/sphinx/source/vp/pydataobjs.rst +++ b/doc/sphinx/source/vp/pydataobjs.rst @@ -255,3 +255,36 @@ from the API:: "theoryid": 162 } total_cov = API.dataset_inputs_covmat_from_systematics(**inp) + +Loading LHAPDF PDFs +------------------- + +A wrapper class for LHAPDF PDFs is implemented in the :py:mod:`validphys.lhapdfset` module. +An instance of this module will provide with a handful of useful wrappers to the underlying +LHAPDF python interface. This is also the output of the ``pdf.load()`` method. + +For example, the following will return the values for all 100 members of NNPDF4.0 for +the gluon and the d-quark, at three values of ``x`` at ``Q=91.2``. + +.. code-block:: python + + from validphys.api import API + pdf = API.pdf(pdf="NNPDF40_nnlo_as_01180") + l_pdf = pdf.load() + alpha_s = l_pdf.central_member.alphasQ(91.2) + results = l_pdf.grid_values([21,1], [0.1, 0.2, 0.3], [91.2]) + +For Monte Carlo PDFs the result of the replica 0 (the average of all replicas) is not returned +when using the ``grid_values`` method +(or in general when utilizing any action that computes a value per replica), +for convenience, a ``enable_central_value`` is provided as part of the :py:class:`validphys.core.PDF` +class: + +.. code-block:: python + + from validphys.api import API + pdf = API.pdf(pdf="NNPDF40_nnlo_as_01180") + l_pdf = pdf.load() + with pdf.enable_central_value(): + m101 = l_pdf.get_members() + m100 = l_pdf.get_members() \ No newline at end of file From 4d0c3ff56d7d429f4863118bd38badd70c223c0b Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 20 Oct 2021 09:25:40 +0200 Subject: [PATCH 08/11] moved verbosity to app --- validphys2/src/validphys/app.py | 4 ++-- validphys2/src/validphys/lhapdfset.py | 6 ++++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/validphys2/src/validphys/app.py b/validphys2/src/validphys/app.py index 9dbd70e493..11d1ac800b 100644 --- a/validphys2/src/validphys/app.py +++ b/validphys2/src/validphys/app.py @@ -14,7 +14,7 @@ import logging import contextlib - +import lhapdf from reportengine import app from validphys.config import Config, Environment @@ -129,7 +129,7 @@ def init(self): if not cout: import NNPDF NNPDF.SetVerbosity(0) - NNPDF.SetLHAPDFVerbosity(0) + lhapdf.setVerbosity(0) @staticmethod def upload_context(do_upload, output): diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py index aa221e8db0..5bfb53e2df 100644 --- a/validphys2/src/validphys/lhapdfset.py +++ b/validphys2/src/validphys/lhapdfset.py @@ -68,6 +68,7 @@ class LHAPDFSet: treated as a T0 set (only the CV) or not. It is possible to control the LHAPDF verbosity with the flag ``lhapdf_verbosity``. + This will override the global value set in ``app.py`` For Monte Carlo sets the central value (member 0) is by default not included when taking the resutls for all members (i.e., when using ``grid_values``). @@ -77,7 +78,7 @@ class LHAPDFSet: the introduction of this class """ - def __init__(self, name, error_type, lhapdf_verbosity=0): + def __init__(self, name, error_type, lhapdf_verbosity=None): log.info("PDF: %s ErrorType: %s", name, error_type.description) if isinstance(error_type, int): # libNNPDF error types were int @@ -93,7 +94,8 @@ def __init__(self, name, error_type, lhapdf_verbosity=0): self._libNNPDF_set = None self.include_cv = False # Set the verbosity of LHAPDF - lhapdf.setVerbosity(lhapdf_verbosity) + if lhapdf_verbosity is not None: + lhapdf.setVerbosity(lhapdf_verbosity) # Prepare a Legacy Interface self.legacy_interface() From 039d048066d044d727f480a7758b5185c786a056 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 20 Oct 2021 09:33:06 +0200 Subject: [PATCH 09/11] address review comments --- validphys2/src/validphys/lhapdfset.py | 22 +++++++--------------- validphys2/src/validphys/sumrules.py | 2 +- 2 files changed, 8 insertions(+), 16 deletions(-) diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py index 5bfb53e2df..e407508597 100644 --- a/validphys2/src/validphys/lhapdfset.py +++ b/validphys2/src/validphys/lhapdfset.py @@ -74,7 +74,7 @@ class LHAPDFSet: the resutls for all members (i.e., when using ``grid_values``). However, it is possible to add member 0 by changing the ``include_cv`` attribute to True. - Temporarily: it exposes all libNNPDF attributes that were exposed and used prior to + Temporarily: it exposes all libNNPDF error attributes that were exposed and used prior to the introduction of this class """ @@ -96,17 +96,6 @@ def __init__(self, name, error_type, lhapdf_verbosity=None): # Set the verbosity of LHAPDF if lhapdf_verbosity is not None: lhapdf.setVerbosity(lhapdf_verbosity) - # Prepare a Legacy Interface - self.legacy_interface() - - def legacy_interface(self): - """Setting some methods and attribute as per libNNPDF specs""" - self.GetMembers = self.get_members - self.GetSetName = lambda self: self.name - - def get_members(self): - """Get the number of members""" - return len(self.members) @property def is_monte_carlo(self): @@ -118,6 +107,10 @@ def is_t0(self): """Check whether we are in t0 mode""" return self._error_type.t0 + @property + def n_members(self): + return len(self.members) + @property def members(self): """Return the members of the set, this depends on the error type: @@ -162,7 +155,8 @@ def flavors(self): return self._flavors def grid_values(self, flavors: np.ndarray, xgrid: np.ndarray, qgrid: np.ndarray): - """Reimplementation of libNNPDF grid_values + """Returns the PDF values for every member for the required + flavours, points in x and pointx in q The return shape is (members, flavors, xgrid, qgrid) @@ -183,8 +177,6 @@ def grid_values(self, flavors: np.ndarray, xgrid: np.ndarray, qgrid: np.ndarray) """ # LHAPDFSet.grid_values accepts flavours not included in the the PDF # keep track of which they are to add them later as 0s - # The LHAPDFSet way is that it is accepted to ask for flavours outside the ones - # accepted by the PDF, they will just be made of 0s flavs, zero_idxs = [], [] for idx, f in enumerate(flavors): if f in self.flavors: diff --git a/validphys2/src/validphys/sumrules.py b/validphys2/src/validphys/sumrules.py index 5b6db2d17b..b199af8b02 100644 --- a/validphys2/src/validphys/sumrules.py +++ b/validphys2/src/validphys/sumrules.py @@ -125,7 +125,7 @@ def f(x, lpdf, irep, Q): def _sum_rules(rules_dict, lpdf, Q): """Compute a SumRulesGrid from the loaded PDF, at Q""" - nmembers = lpdf.GetMembers() + nmembers = lpdf.n_members #TODO: Write this in something fast #If nothing else, at least allocate and store the result contiguously res = np.zeros((len(rules_dict), nmembers)) From 226f90fde68dbfeaee74232824e7dc4aff0b8ef0 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Wed, 20 Oct 2021 14:04:49 +0200 Subject: [PATCH 10/11] removed verbosity option --- validphys2/src/validphys/lhapdfset.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py index e407508597..273ba3afaa 100644 --- a/validphys2/src/validphys/lhapdfset.py +++ b/validphys2/src/validphys/lhapdfset.py @@ -67,9 +67,6 @@ class LHAPDFSet: Once instantiated this class will load the PDF set according to whether it is to be treated as a T0 set (only the CV) or not. - It is possible to control the LHAPDF verbosity with the flag ``lhapdf_verbosity``. - This will override the global value set in ``app.py`` - For Monte Carlo sets the central value (member 0) is by default not included when taking the resutls for all members (i.e., when using ``grid_values``). However, it is possible to add member 0 by changing the ``include_cv`` attribute to True. @@ -78,7 +75,7 @@ class LHAPDFSet: the introduction of this class """ - def __init__(self, name, error_type, lhapdf_verbosity=None): + def __init__(self, name, error_type): log.info("PDF: %s ErrorType: %s", name, error_type.description) if isinstance(error_type, int): # libNNPDF error types were int @@ -93,9 +90,6 @@ def __init__(self, name, error_type, lhapdf_verbosity=None): self._flavors = None self._libNNPDF_set = None self.include_cv = False - # Set the verbosity of LHAPDF - if lhapdf_verbosity is not None: - lhapdf.setVerbosity(lhapdf_verbosity) @property def is_monte_carlo(self): @@ -109,6 +103,7 @@ def is_t0(self): @property def n_members(self): + """Return the number of active members in the PDF set""" return len(self.members) @property From b872fda985e3aa52e17e7b595ff216d6b29d00a7 Mon Sep 17 00:00:00 2001 From: juacrumar Date: Fri, 22 Oct 2021 12:18:29 +0200 Subject: [PATCH 11/11] use the correct call to xfxq --- validphys2/src/validphys/lhapdfset.py | 37 +++++++-------------------- 1 file changed, 9 insertions(+), 28 deletions(-) diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py index 273ba3afaa..c572515231 100644 --- a/validphys2/src/validphys/lhapdfset.py +++ b/validphys2/src/validphys/lhapdfset.py @@ -170,35 +170,16 @@ def grid_values(self, flavors: np.ndarray, xgrid: np.ndarray, qgrid: np.ndarray) >>> flavs[4] = 21 >>> results = pdf.grid_values(flavs, xgrid, qgrid) """ - # LHAPDFSet.grid_values accepts flavours not included in the the PDF - # keep track of which they are to add them later as 0s - flavs, zero_idxs = [], [] - for idx, f in enumerate(flavors): - if f in self.flavors: - flavs.append(f) - else: - zero_idxs.append(idx) - # Grid values loop - xvals = [] - for x in xgrid: - qvals = [] - for q in qgrid: - mres = [] - for member in self.members: - all_flavs_dict = member.xfxQ(x, q) - mres.append([all_flavs_dict[i] for i in flavs]) - qvals.append(mres) - xvals.append(qvals) - ret = np.array(xvals) - - # Add 0s for flavours that were not in the PDF but the user asked for - for idx in zero_idxs: - # Do it one by one to keep the index in the right places - ret = np.insert(ret, idx, 0.0, -1) - - # Finally return in the preferred shape - return np.transpose(ret, (2, 3, 0, 1)) + ret = np.array( + [ + [[member.xfxQ(flavors, x, q) for x in xgrid] for q in qgrid] + for member in self.members + ] + ) + + # Finally return in the documented shape + return np.transpose(ret, (0, 3, 2, 1)) # libNNPDF compatibility