From 06dd5a274962f74f29520bbcd5a01398a62a82a3 Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Wed, 14 Feb 2024 16:57:46 -0800 Subject: [PATCH 01/11] add debug --- src/spey/backends/distributions.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/spey/backends/distributions.py b/src/spey/backends/distributions.py index 8588b5f..265eeee 100644 --- a/src/spey/backends/distributions.py +++ b/src/spey/backends/distributions.py @@ -1,5 +1,6 @@ """Autograd based differentiable distribution classes""" +import logging import warnings from typing import Any, Callable, Dict, List, Text, Union @@ -8,7 +9,9 @@ from autograd.scipy.stats.poisson import logpmf from scipy.stats import multivariate_normal, norm, poisson -# pylint: disable=E1101 +# pylint: disable=E1101, W1203 + +log = logging.getLogger("Spey") __all__ = ["Poisson", "Normal", "MultivariateNormal", "MainModel", "ConstraintModel"] @@ -239,11 +242,13 @@ def __init__(self, pdf_descriptions: List[Dict[Text, Any]]): self._pdfs = [] distributions = {"normal": Normal, "multivariatenormal": MultivariateNormal} + log.debug("Adding constraint terms:") for desc in pdf_descriptions: assert desc["distribution_type"].lower() in [ "normal", "multivariatenormal", ], f"Unknown distribution type: {desc['distribution_type']}" + log.debug(f"{desc}") self._pdfs.append( distributions[desc["distribution_type"]]( *desc.get("args", []), **desc.get("kwargs", {}) From 06e54e5689c54e9d454f730f69ceda542d4d2b33 Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Wed, 14 Feb 2024 16:57:54 -0800 Subject: [PATCH 02/11] remove weight --- src/spey/backends/default_pdf/uncertainty_synthesizer.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/spey/backends/default_pdf/uncertainty_synthesizer.py b/src/spey/backends/default_pdf/uncertainty_synthesizer.py index 8b54a9a..2f59ba7 100644 --- a/src/spey/backends/default_pdf/uncertainty_synthesizer.py +++ b/src/spey/backends/default_pdf/uncertainty_synthesizer.py @@ -18,7 +18,7 @@ def constraint_from_corr( { "distribution_type": "multivariatenormal", "args": [np.zeros(size), corr], - "kwargs": {"domain": domain, "weight": lambda pars: pars[0]}, + "kwargs": {"domain": domain}, } ] else: @@ -26,7 +26,7 @@ def constraint_from_corr( { "distribution_type": "normal", "args": [np.zeros(size), np.ones(size)], - "kwargs": {"domain": domain, "weight": lambda pars: pars[0]}, + "kwargs": {"domain": domain}, } ] @@ -127,7 +127,7 @@ def lam_signal(pars: np.ndarray) -> np.ndarray: { "distribution_type": "multivariatenormal", "args": [np.zeros(len(signal_yields)), corr], - "kwargs": {"domain": domain, "weight": lambda pars: pars[0]}, + "kwargs": {"domain": domain}, } ] From 471736f1b887bda733c339d01e8ad6e694fb046d Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Wed, 14 Feb 2024 17:32:53 -0800 Subject: [PATCH 03/11] update tutp --- docs/figs/sig_unc_chi2.png | Bin 0 -> 21974 bytes docs/tutorials/signal_unc.md | 58 ++++++++++++++++++++++++++++------- 2 files changed, 47 insertions(+), 11 deletions(-) create mode 100644 docs/figs/sig_unc_chi2.png diff --git a/docs/figs/sig_unc_chi2.png b/docs/figs/sig_unc_chi2.png new file mode 100644 index 0000000000000000000000000000000000000000..475039a50cf2957245645bb99fd2bc4fff8760a6 GIT binary patch literal 21974 zcma&ObySpZ)HXUpBi#+sp!ASZigZedAfSW{-QDm5>5z~{BqRhx7-9%vBn5`u7rNStLY-@LNMi6NV8x`6_#{&P1+I$`G?{aW( z7~ILAjsre;U#hOGtc>%-fD8B#i;tg@f`S74?<yTH%#|(Y%&+ zV}EOg{6Wa~LfK<{l^aYV|O;_y=3n7 zpPGtl84msq5$ew_3SQGqY8=~C=YZ^h1wetBhMR)Ne;72m$4^b~=>U(n7&;0W&CpK37Twf;@kdb=w{~uXBxA~YeEC!=+z`qq&5(Gv8 zl}r%9H}~zg&oAD!&*z|^NbM3Sp&A=5XuT{nTe}s%`nfc9wo{)E)H%bL{I9sVB%t)} zjGq}ruuPkYmZ|KON00yeEey@exOU1 zmX|^_7wU|NAWTb{3)xMfP#uexz}PeovE!C+$7lf!2t9EB2I^}Ok3XJWa7^GpTCw~4 zr+X1_>h6DfqW&NmFR*7rcIlk2*91@&1IypPdxIfWqOdxsr0k~@mscy%7;WSS_WGRZ z-1+sn@OhiZ3#oB?Ns_LTuuzW`q-ghQ<5u4My7>#WB?TVW)3bf`g^HXxyX^VIp1_%X zzojr|ue!dYSKkxnp^{}^PrgDW73<_U-+2}1$UMPL)DGt9_6#Uw3|gLV?-Tg6GKHC+ z*l{_L<>8nYEnzuWQy!Oa=uImvRlids^`lQ5y8>GkY#44uLYc{k^0BmBva{cWGcT0B zN{f7ic;jF2Xqbh(Ru(&(^iN8EBDXiuGG|}f@eYWTM;ygMRVz7YbVd{v7P2nXJN3ZtWC|nmpi|( z(VD$U6L{%1vLZ_1Fr^ShN(M4q3g92MH1_nndiDEp#j4kPwZ+km$*)_PTCmfpcBmrelb98;vr~J0E_%IYpcGmz#m>YconJ#?OmR|Gv*SHJywNv~!~JO+ z7?Z{K!EaA_V@`tvMA#}EHgnUN0^IO;UTUodk`gB{zQH9rEj89JYNlDVV_~RQRjweWqzwt<#6l&~}vN$>6VRtjh zxKfMEBQIJoNa)DRtoaqBYD1jBR8IYr4m%sxdcUN$g0p6GzEgj5Z%e`<2%AF)E% z>hw}6rbi%bJ!EViJh6jfhfKdS%L-Tw^4z#kXxxO^8%fq?9efF2>a^iVv+R?M%JmPQ z4PPmzAGw!7{bB8@0sEW2!fL@JvI+cqC3--k-Wl|xY?jZj&<46jAce7lO5$g)qf2rh z@KDSm@N65S`Z;v=R;v9CSy)ERwEU+Gopt(h{F!BHa@HUkKfaw!PMstU%*Ql|B{MFZ zC^Vj#uD|b{-9Z~ZKWCvzeagK^X|XpFKgHqjP(%OD(lZ_~P+8Ze$Kq55=kF5ic*O~D!SHvu6jrY)g1%u_-pNCq9$S;*kq+t_sii{S$?YgvKZuJnW=|Rfy`9VZSoX+c zGdp@i87!@&5crP8c~e(xv#8(*_>>Gax|t+ zJDlKpI^Bm*aT*);m2&4#od%M%%_M-5gOS@fiLq7p!*>jb%k=|;c|KSZ6k)EMm0X7V zqoOLDX&pnLRINp);LRL{DAh$n>-x;9oZrDY8|+c4>n?BO zZ2fmRUTFCjK6m0P!kq6wCF2i5_tou%A@L&hyJK+qqZ}l=w<4Qrh~J+~ZQiL;;^Zyf zrmy@{n2fz+rpbw&nj#-ySA@~z4;aC6s1YXqz6ri*Y7LE8!f!lLOzZSre~7)pJ50a3 zmyx@!=QQ10Bb2&uru@jLWR@xV+fG14wuxn6`i=iGGT;}%bgxURLtk#aMikQ@Vx`qI zjTfG}->rv%i9CW-DZaxZ5PW`YvzCSYBqt2HkM&^MtCH{DprDMe=+;99e_Xsw>)mk~ z^3yK;m}7--h%2s|h${`C#T;l6A*=C=^iA}DP`Kw`x@Rq1&44)2<{%$(A{gk;j$Tw& z=JJE?;$2F78Uguwm-7*-3!SN*VX3 zsckhRRNnH6_TopWN$033VH2cyIucqm9LcaI$ErUSwc+ftx%Q3n*q8L(N?5gcDWURO zd>RSOnUBr%duu@J6m>yKS^PMhx6-W^f;*>7m>R2|SD+oz;qkbRk^{nI1%-uInzEJT z7VlP+0l0!o3R+Ftu45m6fcv>fmW^5pZ~AZ4yGSnWaM>*A)s6N1@E3m^jwJeuXKxAA znAHFK;Et8Bu(v%D9Q#RL<%Ns_@mP$ z{=$?rw{o7$S`WJb$JF(Q|Ar+ogb;ig0V#l8HAM6C|4kLli02(QB-I8J8MgE+`nQFy z#N;T7YA`cgve(_StW;eyGlv8MWAKo~IG>(Mmz!$z7^%n2nk^oFpRRmU_>4buP@vGSdgFBPf`v&`#}A;$^i`%0DdN{ zppg+WXa(4#bCIs~Pk@05CC$b2{|>|$|LSSf37MJp3`zC0?h`_+to>394@H;0du5M? zRC=!w!WOhPv@kGT%mv5BdPC3Q^oKX(8=8NQ*$~ITM7}@;4hAiyTT9mf`@b`Y?8}|c zx_%x@62&XYSCzC8xpr}nByQN*TJ9N%h?wiq89jfcUpv|vlII4s@rrGpLGtj9$af=E zd77ja;DB524T1C%C;T3ff(>{Se)afL)^fIbBI?oBAEq{eH{1gZg>ucn;zjFj<~COo zw+7nDV8hPZk$aKD&Y{);s&aO)LZ7C6eVba@yxHtwznYxWh7Yq#yMV%LCxI6FH<*@9 z^PwU&z?I)|;Rw$mK~nc@Jot6<@!Pj*J%?5ElJIMf00VP0nYR7sl|=t*{om%@p75Wq zD5Gv8_2v6AKQgJW{!4%8N2!&}3LOxvqND^Qo)C;A-OAS}ymp(_#(pdt#_#4tTvhaI zNnyT;J_0|!)h5wD(KS1%s2+8!~&u0GXp z49(^r0h`UpsbjQq<_PCT=2|gH-=g3vypzRFx5nnm;`LE56GF9WZ@bx2v?Qk!n+olBKK)HKEW4Hrh zDgI;&-Mx09P1I-+vizmaN!LM0!H*3;Q48$`mBr$$)cbs=W}co<3iBPo6F2H)5P8;z z1rOqBK_EAu6#n5#0^_w{?h$$z7W%yVUc&RzhIFMdIyYq41R_5R6?z>#Jv)WmLbT%TDpSz%Nj7!@j zcO*KbSrNDVoHI>bX=kA&O3KMamefb-G>V6nHUDu_=`NBmlS`s7gJQ{|mK_r3C&pTnDPXfi;SRl$wE!o=Mxx%SSBb`o zCO;aW$)j9RoNIc9ix3_=wuM1yNBhN~jD0KO=fT5un+*~xqNxTF1ao3SH}eH_6U1dT z0G23vn^HVDM{P2nSrL6u;gI|~dM0OUED|)qO4hd=(+`g(cWvTa+EbUD%!HY-YUAI18HEb#mPRw8p3))1>3x2l_)* zja2IVw$Vfv2p&aRF=1VBe2#6}hlmXI#rQLN-t=!tT;ix27if47I4bWVdoHq)lzt)_ z$Bm=t#8SNsxZuV^bHU1@>pCmr)h)u$oH}99oKr{8af;eeW%j;)Q+?=g&FwWaLptv4 z(mY!MvFZxvUOw4}E$K+WOreq-dk+V-%Phu0#vw*x>2Cy%Ih)n^0t)xms^l zph?oOA}ngg*>5ZACjxOsV5=4`*Em2=keK44ASSh>JzBqceORco9seCao6br97r~dj z=SNvLNdFGxA$f;f@S{&+9>gij-A-|P?av+figx}&fU&osYF&0R@naBY`y~b07fz(? zU{=+TXMZxqk4@=T|jro%C!8t1443brIXsF%I~%h4zVyU~-CQ z1R}K6Rrninu-x-OLj6xuw?oRvudQ<_g5KG(5|P}M8>uW&4xEmQT_13W+o&>#FS=`F z$Tfjir&=LgWHH2!NS|t+&r!<8d%Rq&Qp%BC2RS%~UvM0F$NyL_aEw+lD{MVwt}RLAhn` z<)3wjpb>s4aIn?>ZrPa=HK;uBWNB;bm{EprsNW8nnPX4Y z_3i6?Hi5P#k1`0Yz)TN?Vjs~L1CN?Qy6vB6Wz2qM?6!feU@hAW2*MRJn~$3RYsosb z%!o4`#Fsce6PUR-jCUa^xuFz%X}p7u%x%$~f$he}#%9Pu@W6UK(_x^i!+h$3hat1q zEr=Po@VCYDXD@g#7(de)>c%ejaOmdyToU5MlsSONDgE52+=!^9;_j#6U>|(i$fa^c z0DCLNoi0Ga7t9BU;=^#FMa1&jultCr0J#8J6x6cOM2z&5H&>?E{v$X4M_dJ0Xpn`q zB0yv~W>s^yP8Dnf$ej)#n%55V5T*fgetAY(22**PPVr)dhB{w40gi~@3hGR@OZ998 z3yoYGkir(y$zYDebhQ_BCJi=eg(f@JFV*S|wT*rGs=odbj-iuL19&KdTvgu8u#J(4 z_WbE;d3T~QiA?ht%w&Mf>AkGFk%wE#`ln?lE)~xm;}0E#ifDp1?YxiL5&)s#LnQH| zd#1?7d2ui59*EJf*Sf>N!9EPGXHHZve*ht zNdv5yfdg`mTw3hYjJ{&@6EnRfFg)o(?fJ7G*IDiPZ4`AHLxXn`@N{0}J{3$E1>*FX zco|P!A#nK;4-zg0)0~zEp5Y=)5xO4dPf>OO58S>#)iMxI-()MC#%d5Rs=ie)v9+FhUUl2l7fbr1`5{OuMD-K*vuUvs?n=Q)@ zytVPO9oLXhLmx{0E0^|!35(d(#EhCZ@{4xL2%tckN$*1i3eA=S%7^X57g^5&%*}ND zotxmpTY@59qbrxLz1Gmkxb$|)&HGjC4H<>3Lp-e*xN340`xL=6MeE8bE>9b#RNzz2 zu|ON3ZySijX#zu+KkNo_<`%gyWF`>ah=~@En}kv_f-yfh6g#n|j)p zNE*d2FD5k5?J$bp6h9;kqT2%2`J0`TUJ$c;OJZfOj!a1w3_IIdzA`tv!~7aEAz1*y z{mbHgjAhsh_o7sxlW26yLuR^OZ@vbiUIFZT;ID+wq`%t^@IeK2tUEKdediHlOkgAZ$!UeWnta*ON330Li)TT zLNL5H&C8>3x@_A3!2ZBp1#&G?Hr8bF{SjIe7rJ4A7n5w=BmDPhI71t`V?5|Z2@T{i zmu{1+yQ8_8HgJ#aBtPEb{5}3OyQ!9UrEHCsHPV<@PI8Do2~*`lSZ(w-d#ZGYIuq^E zJzK#GosFKfz;kBKy1_D^jP`m<&8^as=Lm%%os}0N04M(H=Wfqo?UL!lM_4uqj&7n; zLb#3l%kjRsJOALpw@^z>9fCY=9fJFSHhfFeiFigp!_J1q)CQ}d;r#_}$8jdt((Svt z2<|oI%Oy7kE=^2TKdiH}{wblpy9{u9hKM8=lIke{V`zN@xa(9;n?5{*^Lx%8tWBx1 zM9)2nzl5*!H|j;WxQB7-Ss87=qVm-OytW^4E8{$BgNehS2e$GMf? z3*=9G7#k?2zwv>k3HPa_6Tkxjz6MvMM7~P*GS^Hnb)FHS6PIS_5hg`e^p@oB+1^@- ze0DGo9W2=Ll)HY8x(`f*@;kp%ZT?p!_`LwuTDN$+m|4%A+}++hk1YR*PRg*@U#Vqu zVn?{F9UdVjP1zipc&bMH!^F`08^R5Te>3mW#GLVnr%utKKXio}ye%xvR}Kgv zUy$(pCx$rg>9V}g0Mw7E6r(HfXBsryhGN%5>d=*$hJNO)9Mt(mY_L%h?;!^=pDSer zNA;miqz7*Ap9#nioMMP*>?*%)T-Ws5-FnrHt*9sh)PevN2X>!9A^EZe2>|?z3^TZI zQJMMzat*76S^LBzlkzWeGb$sh?DV_2Hxy+rkNTxN zaNwx%e^iV(dpzZWUq-=lTDiP3l9(QUswZvb$MqBT)he97YBt#28@EBUd@i0~lZ#Xa(VGR&nI)?bKX+fo_-3v;NT3Dxw4aF9;`P7B{{x81 zOAdYUg@Q~~bJjb=@k3sfg={+Ye@kSg;v~3|aYO0L0j=|2h9&WY^fH)Jf-D|!R6Y57 zY|TW5D)1Ih{D{MkaNnqYXjOmO3hKJQOZi#LKQ(@57}2|7<*Wmu`C&JYuKvLx6!=k6S|3tTD{@%&F z<7+{4?$V_5{a^UL-YEMeSPqwQ?!N=#WRGW!k+7i*QL9j2G{F+7cRY92+V$*zcqX|2 zju1^0PU>-f$4|i-y5bcc>Tzysyw@RI^ZzvjVjw(RihLTenSMmigkEAhGc`9OWsbwF z0y9M;ouhJwvC5!9+A*1asd&ByS;0Uh5taZ-7@}*NJPYZoCUR#wKHX1ncXJ>v`M7?Y zobk4*qBRw4PtY@{q#;1B5FMLEmDpxe0oWM%C6;ol5ET@xu{?D8ZeLN-b94rGz$ zOvBs^3kb_;U8*J)e@{?qr-w2F&{o&O#Fx+}#vgU2n-q(q9$|H3U!Ugv*Ti2hV&kS` zUQp5`4XhQxuxHYkxoz&3DN%>{fGG1gmqE?GPON#pA;>QV+lz=m+_eO;Gt}c@q@Ume zr5-SVZLZGa9|5s9kr>Z6b_f!ok%YEw!gu#H$K)T(IR|_#sZw$Y4TeG=5ASv&l$B^@ z@CGazAe9W`$_Nj^5HgZ*8J(SI4_{~IFI0yhf-0b-tiOvY-)Z@ zb)Zazu>d>MXBA?75_UG39$u0^jYPU)ZMI>y+wp89-;C9VX0+XK7H+v#^+@B(zij|q zg;UNrjX#NrFzw#2krk$SkC&)acoX>?-lc-LMOUymZa(KuC%>J2Am`T0 zjqp$wM^@^^Zel-vOHQMvkMAZ9(?RaG^|qG*A&p-IfX@sbJ`CwKTZbcM`40q zr|7=_rlG%EG&6+c$)4x#ZEpg`4vv706Nc87a$tkT5WsQ1Oc$qsH`|YE>;x7{1S!E% z2uC$|sqpOcSJ7_b^(#Dx&%vjb_o*9iQOAs!jTTBp3HEk0q66JNYgt|rVYa$}k1&RZ z8TE+VHxN4ASDKe{hfaC*OnYSskYJKdCBlpU&+j=)-qMmbW=#O|6 zHgf4>AlwNpLym_bA%HzJi+s8!gTMbuY6hn3mB3H#yP7P<1`aXvp_uH{{uoVIpg9`>vES*+x_mg+lfUjNn!ygMGr4LM+|~l) zsOuiBtd$95*p9)h7UP+!07ML(bAGY1NY3wgJc_>inxu2mc~70vG$rT%c3B9D+OVK# z-pdw=S*rO#fqmf1>G>PLnCEa$errK4;xPajEB0mu3M{C2top}>3HLHJ?u=auL3Q*;g`f#Pt^DywaictOTn0k zC^i1e3CX$`_zi$_&H2|;qhs#PvSBY;dUs8{!{=heB>Bf3T zKhXGjfY-_`(Fxaxil8sIo0PrZ#K5RdUu92aS>3ijx8M&<6Td5h%JVpvKn?aZ!NT2o zH<}41b=2dAw6ZRU0A%=FdsBc%Pp>W3duXSc2=#HvR# zR8gsZdLc!-w>PrC))yM}Q~m)Z9`ARSQ$;CXR{nM%cMZrPm-H=Bl%+tBM&VRszQH*L zZ&FYe$sm}0zR&@1_5|{as#uMZP+tnAdz<7`QWhWTdqzrS-n&J11mO&3AYn?%LgDb5 zr5qNi7i9U=eg%ZV7A3~zRn_=z0$dO#C^UW8`Sa_q4I;1wxWc&7+j780ORU#^RaMjy zUifgZFT`9#ZF`~Ixs^A#g;Pcc)k%t$Cm(`fNL29@`0-K*^qxL&AuV%cbN7WhNzAdr z=bwe=QDF)qFk$gbxeZT&ODyiEM^O8Pan{w?ElDme-=vnG-LaA$%NqOs)5bEoEk#R; z3-BR6Q17BgWDq}vgfOcc+A~doFlu_U-LulYdfoU{^_lL9!nWRF=hC*9q*j({kUXqJ zQP|q@7bB-49(a(&!XseY^D+qs_!3bA7QB5MR)`Ldb8M9Klx4Gk?R42UWHeSJL5c54 ziXtx-Z|?ly_=C0BS;GpweBdqz-!U*@ZS_oVM4b;a!vswlPqw9(6!`+f)9B`BHi_{r z`b&A&a4pUA57$)v6h$-d0}WBBnqD$v5<)MxMmj}z)VD(m`8z32)ePfw4>=aU+to^` z#~Q~Qh9bpbutPfIA33I-b4)t&x2%c=boWsOGw6(5T2ZJ_&)*Zv>TZ91E7@R{Z65Ws zMDzG^KXz^2{BM`!iUYjddG&B@?43ECFw3Kunko5_>CE$==AYG7l)0Wly>UKiP#5g! zvOrgGRQu~~qtK@b@vQHI<8pj-zn4d#=UMYy48&y&sa`>|#4|79m=FXp>W#DTyZ>ZB z=up!#Ovt4l!4CHg%I<)1wS5NFKo-uhR}4=h4L27KP6ubw7#x zZi)%m<73q+a_VRg?u|u`z3{d3E3@*&KPZ7Oz6{?-b;x}X`4Dl2!+TC$F3T=DxUyvD zF)>5S>*{kV6BRHzIQ+K`)i*nq@j7?&v2y7>b0zgykZUsGcoBd&3PEao?Dqugimh#J zRjl7{(mtIAxc!QRn~y&lW7Xe6Z0Jzg@K%sCK8z%siK?w23X`WorhA8HUn4bVnOdzw zh#x9BgM^l|!~l$N2YMg1%Wq5G~ztowOhzdbGewyXylk%a9xZJuQt zHU65TC4PV2A1Lwrt?95ar$8p+0SqL zNe9A=v^`E=3_Gt-R5%sM-;!4%kigWObqVR-I+yT-)mCh^;-AI$48ATqxFnXN7=zqrxBOG`DT~jN=@on%PWTI#M-Fg((?uKlHXXh=@07ItnOAn=<WzSSGbl~q>QAOIW)hRwqJ%ENLj2^@lmFhgYSWhY+w+d@KsJ|09FoLcQ<{D?&vyMBFkx_tJqDfk0M2EMsrWpHi6< z4%DB+mNrfi3e-E}@#{}Drl@(BsXv=NVdHNswU#LO_RHeXERVd`rV!_yvxge6#!|s!?bz@v+A)|sK}quI);- z4@^LXtoMT9v*1Jg%xk9y*uGd7JuWGrXQedR835_h85oTZEKHm{84G>A?hRyS}F%)bDSb~-g9bYb|6T# zYC2fPyAf|Iae28lH1y-#w;WI(;R_6KPEs)*4i1@k7E@Q;1?;8Kh|yEaA15A1c?eNQ zH7*i<{(FL3AmRRv(5aMib`Fjf88)RP%-1E%N=U#iQ=J@yHMQ?`y0HG^9RV7r+>|-3 zRhcdm#-Rv{Wtf|cFo(e&a*pAQY0UNUck>5%9SILUx7#fR|g zP57=-C1&yLaj_XmJvK037L&okvy-5>sBQLNdHk`#acRX9-h^**iUNpG6Vx_jj32}OIRHQNkN^|1@V3{&<~$NH{qr^VhO&k}VS-@@^W%yW zZ0WAVNCRRG<yT>?2e5p&)PWe$#?kTo+zC5ven4(Z9e;IyAECm?s1 z_#%xQ0|O(TrE_xnT)x)9M4gFvHDxXxXtipvvoB%Ed`@M~!qKe0EQbR(w~B2~fm--V z-}H3EBGP^Z324FFUk2|_FmgRP=BfNc4V)jf1n==X0gI0^M6uF(&|T8} zb_zH5`)a&W1#y@^M*pS7zb~yo!jA~usozvGc6>{$F}#(KJ<@~){((~q5c#%>1#~h5 z;ZuprL-R;OpTTo$6tUsaNaY^GQ1DlpCVe2LmxJ@W$iYV+w?O@sEU67|(~1v$gbw!M z?egv~0*hi}soyc&;O93el6m4igT`l(FHi{DU#O?<^6}p@L<`_0-xU|YkzP!!EvHY5 zRsWh@C+fb{*FGZ1oWA7EgN97RY2HgpEB?Rj5=#MKX~j=X*sdwdMk>npq8Tj0HCRk1 zsvgM}YedO#sC4xFZEm=`GZ3zmVhOW>$`g3KIrcu)n@HyEP9#OV47aczCjt7lJ_A+t z1co(4%ysmdOmZHi7huf-op}TSzg&}BEVeZwhEliG7kDs`D1R##E;pn!9{q*QCtV)p z1mG_Cfo^IVM8Rx$Ay>SsWQpIc%B)8+$bzJTQxrXA?0BSD%6*j$2sGpG^(P^yJOV7l zU&O>_hypy~smsWH3-kv00*nLp(|T!NQ&hat5*VkPFMeWTGLQ9J1VXAI%mGjrK##P$ z8p)AVK-5)*)rMhx3_*8lEOk^`@txHXS%gMhF%uH#Lb!WnFji2 zU&p14>tcORP`SY;*1BDPH@lAA{q#xSe5H9F|BE}!u^1*m6@b?w^l?4ySl!Kwx9>nC zb?yTekeNNskbNe2J+*c5H0HT1#c5|Lb*2#mAfmy76>OU|M=8{Z^u21pSX-@aUab67 zAT8z`60-%Wv5Vr7X}^j~E9whyk;y}{6pSybLQ-Qhs3D;Z4~X;vv-e?mM`FEUucm*` zD&GOJ6I`z?Lral^0D%X1(lTFSm6eTc`y}^FCDFR1c1$k@sXOm5$7AZDq%;Gv!0ICf zVTgN%nlByt{0sp%{E8bAh;{n%TZ={KC}+;<153rFor!p@k6CHO_fR8W5ZUkr)@+~? zpFf(5ud zti$ts{DuZbDod5vC-}5|#q5afDJ;s=k2ZaQ%V)dq@XsE#`S>duUfQ8vEy?$eDQ8f> z#euom%oNdQd2ZRF-0vO@ZDm1s4T0*Nx(}ozz0SJJO$L1*I~yF{&4hRa7uJxFrT}`W zMxDtIs1NXw7u0KMXgRD)IV!Mw7K0h>+5(5cw2*8}f}eaw5ZELl1X@%+ug2fcjTK(m zBqj@;!8QT?0gUI+d%_jLY*}L>79;ZWN9U2}HB>kWOcfDnJ$@fCf)?QW)2LEKKScBhSw|?s505ON+n`gSd;&@s1za=31fLVc-PTXK@IIn&nScC#L zedJr9Dm=V{_!{+S5~P~v;m%&m&oTX#c&E>&GmRxrMWlp-betmK*p04!wZ4sjE z++)#;7tr#%F7b>l1R|7h9OGoRKOCPG(Y18zx%Y*K9T0Cld?-Fw$=(R8OQ@9AH#TpS zr`T=N^m`J4w0<6hs%0K`Uwz4VoPMuZQe|Us*@F#xm*|xO+$G#0nOQ0G?OM`8_-d|23NvWzVpe<^5{DLo>EC1AqW3=oD;ODrkJPh+sceb^)bR`W<%- zW@&IXZGt2RGIgeoe3mM8(6(=3WveJg$Yy;oNE<(sEdrA$zDDT)I0S=@Yc?h8jGS7a zeS<`V6*ls^g)1R~4KZj1?Zh^J5$hDEZYx$|6IQ0^DXVtkIC=+^Iewi5*Hf1fPo)Lc zw?+FPLZ?RE$=sO_#Gy`fn`!Nj>D%lb1RnGz+T8IFq^Oi*I^2AF>b@QRt_~Mr=N(jy zJv6E%Bam7R_T?bQ8OLh$pFA1g zo-!n!C3{t86JrAD*ASJ#qPSNo8u~5noe5(uVXgwEWM@QB#e2w&GZOd5=6g6xZhft< z;hC`+ozr`hY|CAu&vw5Wf`V0>UoW3~lcNhiFz2ez0QQqUrtSeca2EJUyc98CD{{X_ zeX5)ST8n2A^M>AsAK?R32z;7^PdehLTyMZ+Sf8VkwVv6Gzk_$?09w_eZ(!uQbb}vE zM_R) zojF>7IG?m6C zXvksfzP%Lui*R8YFhfmIX@&Xg2v#=`9+<6nLMLKNbWy4^YZX_iwmI?s@R_bcEj!0X zZ{k9s=QzyHQ^hWvhCX52SKIDpHKM_W%0Uo7FC$kBmTD$4+-iO(E9A()*25I^lFk)Q z@?Jz7)H%QQC_(*P#AYuwG2mJ*bX?xBDg2!x8cpQpgqBV<0K6I{yc4{Or=7@+XO6QJ zIqw+IJ|Tl(TlHk?GP-*%*BP534RO@{ky~ja9ke0-kgdKi;7jneghL*H?Jiep*Z+`++3we!JaxK==~RaZ zr3M3{*HDrzNeI7+cwOfm)b%i>ZEK5iN<#XE9JL#qZ>Fr7#lVt+Iw`NU>&eXdDf``O zd2!~)6~us4hs!O4`YFe0OoT=FO<$TsIIpNI6HGE}6WvbCg`d>^Es?F^xE*M1LfZR@ zoGLuX{2}WT4*`7{z96gVZ|h%0^T+U}NdcRGixIHoMa7v+kiv2N31Ht{6y?dV6;rD4 zmYn(&T;ywGGiV{Z57C%qH|gsji#6@XZm??zlgQ~^ofY=Csv7#-UWH`4B$rs=@Kgzk=KzELR^-*rlPC!6`B`~<210(Nkp z>F`8q4OGuTqo$$X=zhcb8T$=GU?-J3i#r{3IZy6W?_@o(G*&isLMtdGyDn%ak|wob z&mm0iN1*DZwaYa=VzK(>1P{4PaXM-yZB**a&;OAJ?S|ck{ftMG47bdx$!AAzNDLK2 zZWkH#2~h$KE@|ANPW?k%Kk)d|k0~ap_(|4|?p8Jz#L}L}%ek&14qF{yBOWjR!!Sg} zy-FBU^l2Zp&|#U#ey^?ughQJ_|JKqDfs{=zhoSjcd0A!>z?SP;d&FNH||a zU3&oDfI9#GjSClgt`y!Y9#2Mk_VbTDa>n>n$i8WHbV=ootdjw4Va%H`?0K7RKG!Rog@1}J#+Buo5 zQ!4nk1+-$JB`g{91Z97BGTbxB5)WF#_-4Fdl(6K5B1=LbR{HQvS z00crfjuj{n7MbMtOo?k6XuiZjH-lv#JFAmZ=3qqIm|@8+g~=+2P{d1XePPRZhQ)PG zbBci!m%Wz~ov5uF)^F)m;NU^WYf7{I$rZcHFCT1F05r`}xpmDU0?K~FCy(uq2Rrss zWGMt=Kw$|gA>q8kXV1+Z)1)>nBEF+NS?Q1qS^a@vfijrl=tvBDbM#ttZ%%mvC`udC zv@WsU$;a|Y zo4e|fzZ@VGWVvrpnL#~})gp_Y1+RolKhelGI!VbcNHN7M>Gty?-hHkwI({D2V?P=E zT7Y97(yhDMTB*2>?g0E@Kim9K{+CNK8$hsJI>6hM4rs%OP{dkYPxw*zicJ&LwY-@0 zjJ3#)3DJs}gidX88&eD{NJlFKh&S)OF%8TA685E5F;F8{G%f>+TV+D-G$0CZatoBh zP=`;!TE3ULSa~MRim?`WUigOIC8OM4rIZ0WXTo&7{Ur3LKzqJq*J~&A`kUTFqE;!; zb(NOcQj55v?1DeON%*h>m}Oz=P-#6nn0gd)pVmcPJyG(T-8YWHI+%A*|nL+2nbK;mVqEX&eLZFq2Ulwl~yT4?tN|t?1LiYw=*vF1VdCL6V(Ib!S8L zIYnIb!r6=d-;1+kYwzBsAX#npc~3X$unmZ{O4osHEWta~Dt>A-DN2LHGz0kGaJ=UhwLr#uK>5TeQtT>FaFa)w_L#>6nYO)rT$2`EKsS!&fw#@mzlBA$P@? zu^IbN0ShvJESQt%PUfJmX{{}br@B1e%CI{t zWH(}~g0jFC);~Q_9H=a~u*SJ4ov4V#Duer(goI?ed!)=k{*vY%~= zz46>}(8A5Z(s`7WB|bLe+LrGzgKFMDO)YOXsx4(pD#q7R1|(d6{NL7`Ev87}I?0=p zN+Tw%7^qplvV__Uq~r|qa%manN*inFOLnc-ZGpD1a#CyQuTiO zNmiK*e+1xXZ3M^LAC`GzXb@}gNzBq0t+HKH-;fR|Te>6{8D^S0ho*E03erQ9eFU~|y&sL{l z)}qI&P2s6i{HHDI6Gly33RNO245+~cAMGBNwv@(~<%B6*UhfVXX`wjCc8&#h~E|XRe`2vehSM>lZ6y4tQ zk<;Fd_;NTgl8ElJdVtKV(B*~4<;hzM{|-D$(SnvpP9fG+Dl#1)2-X_`kwZ{?F}+f# zu5djala7##>UUtk#}c+`a+{3#t*4asK%UkI8!fSB0DHUlR{7(^YhN@wujks!&cUz; z*umPuodujJvV(3Ihg(gn0=5`98LXqvWA%~WpKJQt7t&MCh_9dPUB;dxIJ1-btAo~J zMCs0W+|7#psSmOlO!Wjg3^atAA!>z#5~tNeL3^O@){!4hJ-;EW6wsRh#NxM`&Y+%t z)$+2Cw#iQl#u`^R;*wa^3SYr@CXvhJlV?5Sg|66NlG9@MgB0<^T)0kf05bI{xP*pZ zMrxk0Q~s@!5L&}t_`44&7aZ|bTcV)(7Dtn<%iv#*n$v{S)Nq z_NHH!ut!*-h{`t+6wIWR@W2k)H>O_vGr7%+q^Ll@4D9<76#yX>P@3YXe6TXoEW3AL zr)#-e$6|_ySUJ0B3)?hfN;%kh>F~on)ln?>>Hefw-pZF3z&m@l?buTW&HKj@@8Sxz za#t*X_c7op_)WO)l0M)i;*?o8xva(fCN3%%C#tiTKq){%t&E?p*i~|H?$x`>itWAa zZrNGPsaO2GrTD^nCrS4YlfB`KVFW4QmC2)}<~al6$2(ilruCV6-*YyTgKrbldA* z%d;Ki5j!jTa%R(MrQm7`*h>b7-wg*KuX@(@hTODbnJHd{q)Lrl-rn^N=2Ouq#`knK z4r(X^lwMSX6cck72g)5Ka2KL^(g7%xndH*TQsx3Xx<9h{vNvt_!}hTca0}o`oVE0T zBEEmF`kjWq7*HMP#N`_7d3agb`5O-R8#=4+4WbV>>8@R)z-*(gb2EHRkXGh7BIA*i zCR~x}#H)DfM)gPL)+MKk?^d*z<4r9!#T(GApsF88IJxR>?KSr(6 z!q_9bFjUAO5)s9aiR{_NGRDZyl0A`q&HgJRO9nG$-aC5#=Y7w4xAUGe=gge%_dEA{ z@AEwOx%WQL=kwr6J#&eqwrO*t+b|$n&CW*W%o$bdOe7HPpV$9_fD^4U(u#SQ=?Txf zt_ZnDM<3Xp49H~%i%9rAB1cm>NltG-0LobV@m`zxR$NQHlJ1rtDyST6`xJiM$)H$U z3MzIMZ3@uYtRNhFC`!+b{RB_WGajoc78LVD6uw=|Wr=x7W_wrAJA%?Fc3(pGL7w$g zxR%)NK8fnuV>n7FuBIw)UtYRIUv*;051qotAMxTX4|0g}B!Bn|2cdn|ORVJZ=*bcB z+9!J*5G^*YfTa$JYjS3t$MBB_yfj{nKm9d3XDGvzg(O=zuU9{kr0&%W#NWgkgTp?E zayp%k4-I22g3DS#CJw)eqlCg5^pCC=e47WTjPB<70)Mif>R?CIC^te(pgIL;Omf`7$x70o@g=?>Cm>32mQFWo5*2f zDaGy&215P|e@z7ZGe-4tTK!!4W3NrV5!Zf8yv#L9=ZN8NdBg1WpAT?&qU+D%3|WjW zhO&7>a#S6_e{2SGDYo_f5VoqC(L?DXS%)P;m6ee423&blhXpkI~V84ff?+q&_Bf#;Onw z-8xh%)Qn9Naj#KL^weFRBFf*NK3OZilMAADs&-ISt?D=w%>cDM{iHTL*>F^AgcP3m zVVdh|=(mE0ccqQ@wszzV9)_B0cMmcDab}R~EP>^2qj30s*YwD$v0Y*d&MC5Ooan8G zAr|{(PZDW%(l;)rgy38R+4br1_YIMsy|`v_W^2n3mvxc~^DZ73Cv*DhpDyZKt87+_ z>UZ(HCbpQ~1VJWYX&m9^$>S@5q-CVV0_q-duj5vmZN=rX>b?|0i}@z-{DkAbLh|y2 zUdI3w@(~bdW}LxQn)W znk395l#I7e1Sw|*ztBY0oWh6_@LH7`3nz4?Ks1F3t@?|4+?;Uz>T=DFve|q(F<4ujf^dU=?BhjkV2zIWC7Z&D2T!iLD zf5q}&H1^%Fa3dI;yTBPK8NnU7Jz$lMJr{0jVyfGHk-3*Cn%O$}^%QcSxF0iJLUwM< zeBe+ktY*q->5r! za-~(W99tz`fj4VAdeP>+7B1Mj!F=X-uD`%bSzumZ0F1}r5976dIw z0+X7wszM$tS*Qi5L-I#^`xYc9`c{Za2SGc|rlSbhx7Clns}A8!EYo8#HRkE_WQB`6 zKb#sLt_n2x(AsFxV-aZ66Tb)Ta}Ni%FV4-4nBr7b%+>7spicHZ4YCMp<)r*3of(0d zGSadZG7AKLP!3vRb)+k#=~DZeJ&O@)?(pb*Ght`WKO@%=d;YhQzsMp*BfPl9iuNPe ztoLBFw&6Mq3-|ai&1_-G-@CGeg7BTaa@(&en*nn zK9By#T+w*+-(?7^3}`ZnJhXB;8&^JQV;4eNx*|5i4v#C-T&cqI2Q_8twjRY%JzTs> zDz%P_f=Es%y{Ms~r)I~uCIJm;OB+a7ITp0;Tc6c{*UDYFq*DkoT);(31@k{F^)x^{ z=Q@PD?@?nh@`iHx`O^`4Wx)@>v;mU9%wFHSlVwkS zo%(rq+wIha5Z(~N5=vtuS0QL2yC|+CX!z-IET2~%D^JZCx`zEbqe7#I?&;|iAKVVgt zxE=-YdTWaB{{FaXxjwmF>@q!a-(dNZ$KHnQ>rDaSKSC4zr%-AB&V*)xzMIZX(`o$= zQbS;0WzX8dLP`cu8b02(Nbw3NeKMHc#&T0GbY!09zjc=;`?@b(50mfq9{=mR2D#zI zX>L$|hOxPMF_`|9i6Wofo>;%J+hioa`n<;;Q*o-S;1X$CBpIBo7LbuN^o?s zFWh;!fW*01AlSk5dTMr@x)RaYO|w_ul+u(=(aDrzb zZOg(xS1>6RY@`_Vsf@ep74j#uHkUUuTm42GVYXYVt<-9@@~Zv}Y+mT6Q{6cMDiTk*XQiQyNWIi;Qp<}eQOlts zUIRC^iYh&DWf|m?)GPd|jpSE`&}XgHH+9p%&xxqk%=!$as|g< z{J`bbt(Es;#w)4YVL?>rL|VTK{aqlUh1i|FI4h z6h2_+8S0MAjJRQg?1F4ASGu#OdX8sC372&|bw@=VSFZ^Q6)S1liT@hPJ()^ujTY~5 z+9xXH3|PP3n0~rPiFZ=KXzNIY{I*I{`O2^KmYThA!!<8|M^`#Yexl*{=8$z_V112q zU6|1@KW4yx_g=lcwPD&j+t0gkXZ`l>ab&`6_1kp~R(3{>)ll_swvU27CNEUket5I( zpU}g;SnbfNXKhLCuWz7y1lAj_0N=UG`BoPTilQ;P50FC}$zA&tgE_#c1;{zc|G?@T za1Rl4|Du|N#fxE_^b#3R-_Ara3?L80%pXwRVnF5`T>MvsyMso#xqq|C4}Tav``!zpuY9PlTU^^mlVoIv3 zmq72aWIN#o=BC)mV9SEI*8Fn}00a_~kHm8zWX)ntjbMj=Ev70u!${0$UGFyWcJVahR zc6wL5znhJj_fg^_8p$Tj$WKDsFb2uQ006%Yl?U}qa=*mhkg|BU7h;jvhb>IMt#E7L zry_44hWQ6%kcXBWi}jRHf~KEN_fuD9e(kutRZb^%Qut*=lCE*Gdy_BU@W9-tMGx{RNF)u( zjZWi1=yojZ*x-pmcU^t6JZ?9W^BrFznx9I1L7t_ASP($R>L{mSMV^a_KfA-n(r${smV`-&PKLr_?j$~ z$Mtdgv-r5HCY(A_6vzCVXP%sZR#pAJqa$MKag4Jq5*Tdrxu|Kpn~9XuI|YKfkH}#U z38lM|tOn@!lK8Q=!BL!7vDwCS#^5@m)1N`E;c6WQ+UH2Attj4pK5^qa8SajpIX&De z(B^!a^w^}vw%%G5Fw}|pzZm;&QyjK|KBAn|6Ak#pZi56CQ2B zfqHVp5oPqjOhTBQuY7cchzt-iE4I|ZssNH`lG}xIH1&70Y_KAToA9`_#7hIfJT)!) zcqJ3x=ixZsw6iuZ(IzPYKo|h=ArqI8ef1j=kd+H0fcL?ksy5|Sd44NSB_FXKgWzEF z=Ii^qZJp{XoJu*bAZWxMn7@b%Zl%)82Gt1Vzj%_4tl~^L ze92V~C%B;*m;Bq?n@-vtKC@E5vRFrPQvNIBWY=OX{u)l8fQQjbx0VddDBvoHX93Qt z@tv(5L!CnOwL*a;>4Ud^nco4(GY0t6OmZOW(So@V<{x!jpvEfGbzLK!(m(AU F{~M&HxyS$j literal 0 HcmV?d00001 diff --git a/docs/tutorials/signal_unc.md b/docs/tutorials/signal_unc.md index 312e5ca..ef0ec2c 100644 --- a/docs/tutorials/signal_unc.md +++ b/docs/tutorials/signal_unc.md @@ -30,35 +30,71 @@ Note that theoretical uncertainties have different interpretations, we can inter ```{code-cell} ipython3 :tags: [hide-cell] import spey +import numpy as np +import matplotlib.pyplot as plt ``` We can add uncorrelated signal uncertainties just like in uncorrelated background case $$ - \mathcal{L}(\mu, \theta) = \prod_{i\in{\rm bins}}{\rm Poiss}(n^i|\mu n_s^i + n_b^i + \theta^i\sigma_b^i + \theta_i^{(s)}\sigma_s^i) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta^j|0, 1) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta_{(s)}^j|0, 1)\ , + \mathcal{L}(\mu, \theta) = \prod_{i\in{\rm bins}}{\rm Poiss}(n_i|\mu (n^{(s)}_i + \theta_i^{(s)}\sigma^{(s)}_i) + n^{(b)}_i + \theta_i^{(b)}\sigma^{(b)}_i) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta_j^{(b)}|0, 1) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta^{(s)}_j|0, 1)\ , $$ +where $(s)$ superscript indicates signal and $(b)$ indicates background. + ```{code-cell} ipython3 pdf_wrapper = spey.get_backend("default_pdf.uncorrelated_background") -statistical_model = pdf_wrapper( +statistical_model_sigunc = pdf_wrapper( signal_yields=[12.0, 15.0], - background_yields=[50.0,48.0], + background_yields=[50.0, 48.0], data=[36, 33], - absolute_uncertainties=[12.0,16.0], + absolute_uncertainties=[12.0, 16.0], signal_uncertainty_configuration={"absolute_uncertainties": [3.0, 4.0]}, - analysis="example", - xsection=0.123, ) ``` Similarly, we can construct signal uncertainties using ``"absolute_uncertainty_envelops"`` keyword which accepts upper and lower uncertainties as ``[(upper, lower)]``. We can also add a correlation matrix with ``"correlation_matrix"`` keyword and third moments with ``"third_moments"`` keyword. Notice that these are completely independent of background. Now we can simply compute the limits as follows ```{code-cell} ipython3 -print(f"1 - CLs: {statistical_model.exclusion_confidence_level()[0]:.5f}") -print(f"POI upper limit: {statistical_model.poi_upper_limit():.5f}") +print(f"1 - CLs: {statistical_model_sigunc.exclusion_confidence_level()[0]:.5f}") +print(f"POI upper limit: {statistical_model_sigunc.poi_upper_limit():.5f}") ``` ```python -1 - CLs: 0.99563 -POI upper limit: 0.51504 -``` \ No newline at end of file +1 - CLs: 0.96419 +POI upper limit: 0.89459 +``` + +Let us also check the $\chi^2$ distribution with respect to POI which we expect the distribution should get wider with signal uncertainties. For this comparison we first need to define the model without signal uncertainties: + +```{code-cell} ipython3 +statistical_model = pdf_wrapper( + signal_yields=[12.0, 15.0], + background_yields=[50.0, 48.0], + data=[36, 33], + absolute_uncertainties=[12.0, 16.0], +) +``` + +Using ``statistical_model`` and ``statistical_model_sigunc`` we can compute the $\chi^2$ distribution + +```{code-cell} ipython3 +:tags: [hide-cell] +poi = np.linspace(-3,2,20) +plt.plot(poi, [statistical_model.chi2(p) for p in poi], color="b", label="no signal uncertainties") +plt.plot(poi, [statistical_model_sigunc.chi2(p) for p in poi], color="r", label="with signal uncertainties") +plt.legend() +plt.xlabel("$\mu$") +plt.ylabel("$\chi^2(\mu)$") +plt.show() +``` + +```{figure} ../figs/sig_unc_chi2.png +--- +width: 60% +figclass: caption +alt: chi-square distribution +name: fig1 +--- +$\chi^2(\mu)$ distribution comparisson for statistical model with and without signal uncertainties. +``` From 96423c6b051b16f6a5f1de03c992ec4c7bbb103e Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Wed, 14 Feb 2024 17:34:30 -0800 Subject: [PATCH 04/11] bump the version --- .zenodo.json | 6 +++--- src/spey/_version.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.zenodo.json b/.zenodo.json index 5d2c7d4..cd0af80 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -1,8 +1,8 @@ { "description": "smooth inference for reinterpretation studies", "license": "MIT", - "title": "SpeysideHEP/spey: v0.1.7", - "version": "v0.1.7", + "title": "SpeysideHEP/spey: v0.1.8", + "version": "v0.1.8", "upload_type": "software", "creators": [ { @@ -29,7 +29,7 @@ }, { "scheme": "url", - "identifier": "https://github.com/SpeysideHEP/spey/tree/v0.1.7", + "identifier": "https://github.com/SpeysideHEP/spey/tree/v0.1.8", "relation": "isSupplementTo" }, { diff --git a/src/spey/_version.py b/src/spey/_version.py index 9a13fa2..e83e13c 100644 --- a/src/spey/_version.py +++ b/src/spey/_version.py @@ -1,3 +1,3 @@ """Version number (major.minor.patch[-label])""" -__version__ = "0.1.7" +__version__ = "0.1.8" From d3939b9e668f300172c67a023c746dffab197417 Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Wed, 14 Feb 2024 17:36:43 -0800 Subject: [PATCH 05/11] update changelog --- docs/releases/changelog-v0.1.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/releases/changelog-v0.1.md b/docs/releases/changelog-v0.1.md index a46edc2..26ae0e3 100644 --- a/docs/releases/changelog-v0.1.md +++ b/docs/releases/changelog-v0.1.md @@ -66,6 +66,9 @@ Specific upgrades for the latest release can be found [here](https://github.com/ * Bug fix in signal uncertainty synthesizer ([#34](https://github.com/SpeysideHEP/spey/pull/34)) +* Signal uncertainties were causing a narrower $\chi^2$ distribution due to the weight of the constraint term. + ([#38](https://github.com/SpeysideHEP/spey/pull/38)) + ## Contributors This release contains contributions from (in alphabetical order): From c87a3276a9e796f79df5d171cd3c5c7a67ecc328 Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Thu, 15 Feb 2024 08:16:01 -0800 Subject: [PATCH 06/11] use version instead of pkg utils --- src/spey/system/webutils.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/spey/system/webutils.py b/src/spey/system/webutils.py index b67c5ac..51b77b2 100644 --- a/src/spey/system/webutils.py +++ b/src/spey/system/webutils.py @@ -2,9 +2,10 @@ from typing import Literal, Text import requests -from pkg_resources import get_distribution from semantic_version import Version +from spey._version import __version__ + log = logging.getLogger("Spey") __all__ = ["get_bibtex", "check_updates", "ConnectionError"] @@ -101,7 +102,7 @@ def check_updates() -> None: response.encoding = "utf-8" pypi_info = response.json() pypi_version = pypi_info.get("info", {}).get("version", False) - version = get_distribution("spey").version + version = __version__ if pypi_version: log.debug(f"Curernt version {version}, latest version {pypi_version}.") if Version(version) < Version(pypi_version): From 25cb52cd91fe202b2a0a2ae01884c40bf993931f Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Thu, 15 Feb 2024 13:07:38 -0800 Subject: [PATCH 07/11] couple fixes in layout --- docs/tutorials/intro_spey.ipynb | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/docs/tutorials/intro_spey.ipynb b/docs/tutorials/intro_spey.ipynb index bf47651..21016fa 100644 --- a/docs/tutorials/intro_spey.ipynb +++ b/docs/tutorials/intro_spey.ipynb @@ -257,11 +257,10 @@ "id": "99dcc3b1", "metadata": {}, "source": [ - "````{margin}\n", "```{attention}\n", + ":class: dropdown\n", "Notice that we started POI scan from -0.3, this is because this function is not defined for $\\mu<-0.3$ otherwise the Poisson will get negative values. This value can be checked via ```corr_background_model.backend.config().minimum_poi``` attribute.\n", - "```\n", - "````" + "```" ] }, { @@ -542,8 +541,8 @@ "id": "e0f22242", "metadata": {}, "source": [ - "````{margin}\n", "```{admonition} **Question:** \n", + ":class: dropdown\n", "\n", "***Experimental collaboration provided asymmetric uncertainties but not third moments, can spey compute third moments?***\n", "\n", @@ -554,8 +553,7 @@ "$$\n", "\n", "But this is an assumption, if collaboration provides exact third moments, please always use those. \n", - "```\n", - "````" + "```" ] }, { @@ -767,7 +765,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.16" + "version": "3.9.18" } }, "nbformat": 4, From 8df7c62e6fda0a4bb9622dd4f39b09d445223322 Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Thu, 15 Feb 2024 13:16:18 -0800 Subject: [PATCH 08/11] update fig --- docs/figs/sig_unc_chi2.png | Bin 21974 -> 22288 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/docs/figs/sig_unc_chi2.png b/docs/figs/sig_unc_chi2.png index 475039a50cf2957245645bb99fd2bc4fff8760a6..227a7e21e741e8528bc2d5766863a096fa1ef702 100644 GIT binary patch literal 22288 zcmY(q1yqz>)HXablypjm3@y#jsW@~=N;d-%(jd|b(jX}*0;7Z=B_%b2zyl26Al-~)NeYPx+;Q>q^YWDy z7WVj`2ZTI*oP~Ke=`yPIEdnNsOv5XdA{1AgD=Y5vZ9V7Affbl>mY zb^#R%Qwj&E&v1vNcw0}5ijKI>7b*wwxaZX__G1x@!A7_aaV+AZuJ$w!zl-bFa7)_i zm)PS8OySY$%LaFb=R0TD4`$`e7swcf6>jCajRXznKF^>XX22)G>jfLm zwGX2!BEhIsluKYCrt3S>HCB>aszJpt(#mJAxy8s`^n_E>-S>!LO|!Or$qV-!vK`NkzYl_4f#QV~tXN8z}0ukoJD$p=Mmih=HrmQNi$esozJMQ1Re)^Jyai5}Fg>i(mWyw;L-|p6D+rnHcQEp5%ZwLI20&wcd<*Rd|@kzC)1pE z=032koOVTcJPYyZxox%M7-y&H(lWyD5$p()9c;aYOj8-RCN%&ao}yxVeCV2 zEP>^ejI_t-Eo2MV_A^-&Txbv(x+t@dsHM)n#A~zDt6| zSs%zFhUU!{YDZGfj#_DQN8n=Pg)TxBS*lP?yZ&+K_rRle($RiGJ`dD*;L`>=`Y(+P zIa$grUE`GvPqvY+4ryUwl#ky`Z>~x`18YU<%$-%P6sk$6Kn>`hfRhDQ*6C8$SoAU| zz!!XUg!+wtkE^9%CA5N&Dfvh$1bJ2S=6Q*J|* zl8uGG&sm1qCJic=y$yDI1Yk_xH%=*F78J7|NmFDn=zN!q_5Y@v#RT$B%F8DU`rXh= zXL-R&7QyX;VSR8;hZu^Z9h~eR!-E?5=>{1O_7e&izDFLVCvy<9wN$_PbH}3m7rASs z28y?MEzu4kIXK^^>?aCY108hBP$ZF6NYrbW>XhNV$9d81OGY&7kk|U)>4LefJO=lT z8r9Je1770ru|G(S!heYc2d(RoH!X`*jGr4^yZlJHoKo$Es;o$OMgD5R@z|Ksp8rwc%f9=+F|p+4SP2EYI$quj~q8tCw`L3 zD&2^}%`0$TteLVQSoYh8%am2$FYU2Ud}n`mkP`3CHmk0ZX(#XK;UXlwug%{|zF23V zV6B&ZVOrtf2^zK z2fBK(f={zoO@rXoyX4qF`;tqBsO4`(6HrBV=Q&>6x*XEIEVqPBgY~>=-oR?N<+)}s z{Q9hKh+Dz%nZr;=GOaew%WN}aTJMHOZr)@Zx&7_^PrIBAa1u#p`MyWW!ajJP@YU|u zCLPa0XHUun3sus^g4&U8IdzxGaq169|4_9U7MgM<@$zX#G_~_pWH*R!%O8b_>+1e+ z!if`lrsH!mjgWk^VXuTnBXH8&j^1Za!|A%4U5gQK*M>dyo6Fv?ur}q=X1(n%%F2zy zj~aR`IH#0JlT6aImRzbmEdOcT+X&NfQqEB-=5lldu!-ais6-sBmelM5*}a0FFIV1Y zsd`6Y6kVuu`FgPlBQxZ4vJAhJGtg>F)b`JK@;m}p24vc?@8^W6xbQHfn6V6l z8m0*hPcN9C&uvcf^s_c0)C70x33e>}?cPKBn}qw~<}jTG{`c%OC?9S)<6E!}Gt;+? zeA0ZPT7`RwP> zCAN)4E`WYd2=UtXII~TSwoFg|_l~ogTz4`RWTSIr{(4v zVWC`@GS+kK{`O<0d7&A+``aXVEr$VK&{+C(Hfu4Pl2y1tA)|~sDEOT{1MAvwIrp=| znN4^xAH`bsG;i@qFok}2jy|%KH=Fug8$m6>Ldx;XD#x$#jjmC;7^l2q5xMS1;2;)u z@>R7j+?_dQq2>!l(0iFzZsbng=hz&aBsi7W!}ZuGZC>-8xmBE7j^>X@W6gaeASLUi z9wsyYpBHZs?+2)>i5bNcYNo~N==E~w=T}2SBTW$-uPe%n=6=_uB4aY7I&{#-{uEsN z<@5%0lXlB(6%(x0nUm7B7`G4Xo97w+`L!qRmI8pQF(qCXxV28wf9M}YZ&KJq2&& zv44MF!^M->Re|?1xN^nI1y2d~Z#4uuxP zUGa|4@kvMJwr}%W@(K_9jjcJTDEUMYA(>?(6!3!uomiH+WIyFYA-7?rtVf%MoU86c z`kkpWHBU?`QqcE`?MnRN1)eH!{^ zS%1!7gt*Ss6N5fO<^s@B1fr6(P3^GBpP{B)t}yBa6D8^ChW416r{SDj$U?FEN5JW% z)zRe+5ZV2>N11AH;4nH*fi*9nPzk#akked7LKLq>1>e9c*Bgx~=B$dhcelogcLq~H z*ZpMyE}^av=lV*rTY+cmc|}RX4dd?EjuK)VK2R6`CsG`BYbn{13!ulp(zM0P-2FUv z+yY@`-^bi6W&>_e)kp=J`0>ITxpFC4iwX7EVp*=(YwV`3!Lm|3quzh6G6c!LPu)1!n zzKE9!@%g{Yyo~JbUg$#=1SRkiYsMpA9TDJzLet4O-%~iGt*!1Kj0WE4>`k1HuqPOg() z;>7Ir;6VZnH4MvbZppO7gS8vWF;iNyvK3s?vZbD*`zb6Wc)S!%jv8jT<2Yk0yy`@_l>k=`u>v^LiI9lhevw7X45o+8Ba>f@p@wpQbePJCy}Y4V*-lpJ?kBN#;wv*?>wbJ}B^;V7YZsBRI)xaI2{ zS{HaA4Cg{s)oES3YpJ_i9MFkVTN}hi(MVeojI$C$_Wi*^|KNoZK)AeLY`_-_88Qp+ zgA>iw-YcvLZIZ*ppCAe%VD{^`=Nyy%X&7*OEO+rU(rm&+mxGlrmIT7ZvvUC*fn3U;Bd?nNcb@KwFr zYAjS5SD&O7_~3LeGuYwRag$LNu8w-7EC-AGUscax|G7{3S)IxZpjCs`IHvF^F4x(l zm|UJ*y>0v|MM`gvFgE)-&RW$YO+r6wDNw{+j!*AKx}KCvHrC@t_Tv94oU*FLbZ)2p zT#ZR_r*nXLj9x@mbseW;mv>uoFkW@NZIe_j=cFjoEY{G8yId+7F2ROIwYl0xd1;@$ zjHs<5-cv|eB$CpRLM7B& zr&w?!f68*W3UiY%qv^MG-IlM)qr5mtCuLnj@ockUhZ|U>BBLM`9Ax?>oWYZ0~hOpICz5q#br|@k33sqD>->> zBDo5OAk@A~}wR8P2-Al0kC2{A=<4WbyTxh>citPMED%Mn%Jws2l{%M^!U+o~y2G)Z!$`fsTnE;Ivy|G@q<2G-G7ct$4YI<+b@`EB(Yq?>BGhZyu9l-1j+4{277zlQqIibT$dxWE z9g599W&~h@r1SVTya`=^`$Ae`4#SCmH?Hy_*PF6mEj#ECxGGZ*+=(FDOxMSmqxutq z6l!*+mN;oREW8#dF51`v9ViQGvWUD(x;`7ZmGXLqbs+aUog@G4zC967U5TQ?M@$8` z(BK=1WUGjA7ER>kQ`gZG=|dw@S!I=^^I;k2KlI`!X6ru`2G-88Ipt~6Rr|RSRAmK&z8Se-{&b`=GU`0cDJohHGE6ZfASq=TdRWwgpu+UDWbt}y za0#37JXx-v$>&XFQ8^472bI`|IvcEOKlU4rGfNbA!`5T>PIw2`Jq_@yTNHq~J-1bL z9<4W$^Rn9B6OI_RJLjJ3{Lk5a=N+CM%+Z{}B$q6@OTZ-0&aFG6j zCZfyu!Z?b#8>+m+6^cZd9_pm3zlmZjx`Gp8Qc)bp^KEUgvr?vZhMR%TU(xa@T|;lx zw_Z-G>$pjs2d+nF0L=BDEP=uR@sN&d3ZA2p%irj#X3_Yob3U#0p=9lXpUbg-udPyVY_9%nIRid!xnZB7NWWo7I03Cg6oybp+oneP_oM2w@B zi!L0NhA(0!YX51KbaWYLIPSI;?W|&fAFDSDR_lP)H2I%$S_`3@;&xXel_!h zvH12=tFeD08ZsVo#g9|5kRKAyo-AGIT9*v4bE!<{xguag?(oq_K-2wv&Jg^~>#8WY zJ5m?UQfWvOP^SuGvO~Zk%NZ8c+zOpk-3sGk)mubX3Q_EKC#BdqS~3za_sf00W^&y- zZIv!NFd$R;m>L$qSwPWFa3RAi!l?U0TCzHFHS({$;^g=n6($N+?Sk1xS*mS7?|;I| zLjKtclwHf^1)sXaJq_LVTH#XA{CJ)8w~Yn^p6OJ_OaT`t+e7QK#*US(^3{Vug}b0d&NADUZ{v!*9d`b8J8r3pY4vHdmG zKt6*FSeK+VA#$;YtWfW>hKB&w<|nK%F;3T&=#Uc{dxO;y)GZQP}$B0~HU}|^j zHnlQ}ZzM2Mo`u+*rtx}hco47#0SW1TaOuLV`^I@YMC!5E+F< z_pWOcdJ4+#t(X_5MCJhw%IDuVl(4%IzV<4=&5MccXPu1yzpjRgPyu%5W=q~6bR)sb zsQr{k?=v!i$N!cQ#h}`Kdso=Mulx9{c7JfND(lw2&eA<0tx-GIU~-V{G3jUnoK62_ z1N`mZY3F}+i&K>{)Dp)(;R;GDx8>ytrH;S0r;b;#l|-{i=Tw<+!PM48(3~xz zJt8`1u5YgWs6J{cZ^q)rhJeU!rk2ZnqtG|u7G zHFKKy)ylED#*q;^?G|;wtM0x=mN8~8cK!Z56$ETh®^4>gxGd$v;@A_2qinuSN4 zANG{Fs1-Vu3;&vV`yNU4zJi4HLupRW$?KL14iz{j7~xl9x9@9G#x;X?8o(P&*DSbL z%hiMa46$-hQGm63wv5pjJ3U>HWypuUF6%bg36;q>dVu9h4 z2DtP0Jy`{*?{ZmMI{Ioe^rkekj|72}^s|6LA=0KVK8%HIn<#`oqaEi$v^tW@Q{2ng zy_1K5TXi*!W<1r*1Jgz_Y=&e%1Aqt$BLtWxu=PNBDmg`(zsVe`{5os=Hf_@9wohf2 zxg+MAj->j9946eT&9Wv4Q*rrgn1gx_ds1(gddpw zEA$lAsm9c0{#e7HJKNIRhgXuoX-PLwk2k1)_^82W;i0Yi%1m(${kEw@0pKhsPf5S= zx!Wuc%cJi`@B`;05yGPom_#-*L@3sPqU>QmT-}3;N5O7*_*({`Uw+Cz>0bQZSRHC#uxgV7r8m{{{ z?v1lQF+3Hke(4+dvAdFKk8yAOVHgFLoQKLbasW7cL_^L%Z(Ku(p&C5!T~PVND?{KP zaz1ufw!h?whbG}L{|sG)1ZLSyq#TN|9igiCSuNv$gb84PDKfkO@>y7-x@!cP8h*LH|n>LT2s7=It;9Eg{fX+syFo8sPxaFUNx@kF&%NV|+JPr+o6Ikgb@V^SPXt zW{-+H;C!8<=y4XUzeEZaEg=`VyFG!k{&L4KVyNqwW~t~x)#T4#)}qNuMbgQ1&KEoq zwTSWa`3Y};o0oog@vl<~ggfX0x-CQ?L^8WhrIq0o&R!jvN%;R)P3rwJbT>lm=WuKG z5rM&EK9~4*gy4`!Jtaai-6gCPjV+mJLOa$Ed<0 z^}IE$SisTu(&R>?s@pZk9tw;6 zFr&rG)F)Q(9F5}oVQ^Et7;KbjE1i-0^<1LXp?n>m^Un_H6 zL`@^!AzvTAk_T?8p)ffNdI7=>vBBY{CwuSHGnM#53HdhA_m79Cq42vRW?;#Wxj;J>IWr>Gf;%mSbcMdzN_V6o8)adj*T7f=0B?ow4j@E(vV|kSEqGdu?@iYhb2%#F@p73Imu7_eb+T>m|3j})H>B%w z=f3rhZ|$7_8kC3AJ$JcHMtVo>uSD{2>Ya>y9WjiedfVw>6*(F45hb>?)!dnGnemSA zS8j2AETT#gLQjq+q849pLBvEa88oPl4=E)Ar{FhnLdR7y5al(G7LZ43#R22HWl(Sl zI|KY)DTPU6Q-cst2%e`ttNXm?`?-}ZBON50eIt&Lg6LD$y~&z$R(+ddkI3>Och+x( z@`>9*Y8sO?;3X{ad}vxLx&KMWz3?T@CB@=sEH}kC_l*+dgiSy`Cf_c=PEj%Hjy+AK zD(-sKlsth|?3?%WJ}H7fAg5b%=->c=hFd`B&f3U^hsEn4-liT{!0qLtcrDIeT{2=Y z7b+<-j2_nhaW?z7B4GGVKzH$3W#aykXt^4*TeSp`4y_`^gRqp6e zMwRRL8E2L+^Xu8uvQ>p6&YPQe8ItM)6W?8IDzu z>BZYOC>e^GwpLb@6xsj+=P`mA5YPciU?EP%2tT{+Pq6-Vl7Z0*%P|@jJQC0xdea5b zQT(dqXynTQAOV0}_Atu7wzz)IELN0pgKhxZd)acK#v|C@|0)C~2a2$}^<{3;yq=a~8B=uS=! zm~tHXXLF135CMj|&^Xch%Q#sFU<+tIC&UdY{T!x7{Xzv4Rnj)I$Dugo%mP1=v6S?I zUOS06_#cCFtjI#%Cr<#xyn#N&-;>wp>b76nO7_RM22VgXBdxrub%X}?v%V)!6pf4k z%<~1!Oh3QylNVP$-zD~nH!%E@LJALXl-+>4JcVB-ueF$i*-m-%m2tgZP!<2ctYMt7 zVNp3dxQV0d5|Ugp4|6HQUlO6k_@?ya%n!U#12Sq&m1N0xpvV zHo=|S@~r!zUzx)-*g7ebv?tVn>bHjZ!a5(3j;nJIWGapqT*0mRd;+{i$)#K$InDSC zs9RClo63hjGR&q_JU!dAGaJfBRpC(T9DBlDs5Zhi_I7_)4wIKF$l|y^PWB6r4 z7KPwyQTP?{J^djUz>gU?eq;v-DQs(sdC6>)<)g5PvQa+K6ID0hNGm(_79r&fM!0f6 z{dTmqq+t@p7wMT_OlH|)Ngf*Lxwr|VG{j1n_F7FbImt)(LP?ng zQGnc2eA8fPHF@|4NzJQEEb;vEKjC=$5nHS zVIVB!Jp6zy%P%7_4d6LSdXjwjr=|H})HV1TptbkI30$fZzl&%S%Mp^_(W;y*H%Xcp zhCI-E)Vwijghh^?o^e1CN0o~*PkiXFC$V+yJ?VEm%*j!}Wu{Df)7!@lI+(k>?z&n= z;Rfzzc>h$nS26;J54>hLq0o^e@hWg|raWSs5y)VPIcL%M$}I>)Z3T}Hrt{xZL1Ph( z_=xQ3oatyjz<8(VQ|@0;8-`0E#W@vH#7-$VfolpPzr1p|Rs026{7x()WDo($t-uW_ zB_9RC2r}6lKOErBvGm9;9B(h&WwbNx)nfy38DnZkL`j}z*xaS{TUzNloE4y+Nn^Ah zNH_$~L|dM!Eg?ornuPl<8#p29$*+GWP9Oq6m_I=2P#ku_ao+(k4o~Ph-l^dwj`#yO zgbX@Z2c!*rNT$*%()kJS@PPiLs=F0QH)n%l(I$4w~Q)J-8u00rD6 zBbn!-Q=gHZNs&qW7Z{_l46Ps~($+#iB0PF^ReCMltfvc=>@NV^D|&lP`ZmhWgA1EiH5yy3pYz1#X!xcX3W!sfpD!jYgmJ8yx6MCzzx?*!@q zKIBNYy)z|x&lNSY=ye^YysRO28?K|T9tG4OD!Art&yl{%b^wGITYGCTSl+#{plA3ZGxNJ1%X4>ad@7-r;c0he(2^WYQItvDO64l!)3lOW2gP zJFLKX`1V5xD`ku)%Q>i7dRtD!l2AXWzeirhux>}3B@18X8(R0j`|^ORNlM~W(BJ#O zbg+O27o<~Ao*rsedadr5mTf6sSfFmCax@EEj@~O0{0+YGuY9)=?&M>Tg1FqFfI#<0 zOEX>rF>qh1aW3j2KP3G1nHO+S5$BXnUW*z45A~^X0j?9_hyR2J@Ql5($6O+k0W&~_ zTV`p#z`-)scgA~Bji~lLFXEq=sVW4V7~YzIZE4NrfJ%Ni!yeCD3*J`RlNhQvZP|>R zB$4_Y)L7zfdbN{cyp>Si2zrCe&)@mJnj#SWAaiPQFd3qm@PtR=!KU} z+QZm1oH)rsjQ}(~T%7ap1o2_Lhk`Pd%u$uD7=tL+pL{LXD7p8crjuD(6|mCa?2L0A zH|#SwoM3{WgN&JYF4g5B6#DKLvu$6^<$kBkhO8~)FAc`KY;r&P&@ht#wlVO7m$_o_Nn< z>Xeql0>v;I>;?u#x6qOS+Y+Po zlz!z-pTwHjUmd);IGIFkE&q3~7()e6P5VKlNB1#tVeN2d#9~w^L2^&19W+(r3y`-W zSp_JGtb*0mDfESY$ckS-1~yUTTYye0ApR0MJZEFqRtjwBR!Nc zLCIBWpUd6xYvizlfLq?z14CtU-UG(x^#nn3}n_JJrNz8)Wh0W5w-z z&EYx@7GR^A&3t1cEK!JX{=V=Mg&&4r&Q;OvWezEFg21679=^0% ziZbnQu=(a7hB-U{<$*c}gDruy-Lm~COL?i_`*V3fF*>94GAT6N2JZ4fDk2zdp*1Kp zb2i=_1W-~8guVFK6Yxrx^$Dv$hufn!EV`e#8ByMdi#s_>K#nV&AsFNwHNQ&deRXDK z7UsX5QMi1i4=K_}<-5T!6X6*6GSbyn@!=MB1=j5OByVI?Z>+#&bf&4RW_zRUXK<^oCG1h7UkzNx>*O=? zg-h1CJvtexnfMu^WTR2KYgD&cjXH@8LkS{6Mhf=(T~y$?0fg-AkT(|UyA8nQnJ1)_ zhboiiosXL$KN(h-@>~TAzx0^&9tcA56m7_%eiOoh-4Q`g5uZ9uSsn8vx4cg??j$?^ z5W%*5^>E_QK3i6-`aa+~=J4rJFWav39GbTb zMkgB*yb5!TLJaZ;ZZgDRJ8^R#j$v~u|KmRldttoO;i3v1x`%3@w>*1l_4}73D#%JJ z{18ERr32PHC2bFxn|WM%#0MHwMzI~y_jrY@?(YkaU@vN_oF>R3wmUoN9>*Q*P zS&#JBERdF!`-GlF6so7zZs}zQi^CIONn!6fg0oeGeW-n9L6ZCCptOObIp-Qc_ezhO0lL z>MNleJvX*ROf3>h3$J?vk>Y=eFJc9;f2d}aXtiP$abNPYQw7@4cItkwftILf#*YG3 z=*l}!FC|3;8fMdlKQi=Z`$ihIp2;A#vu=)d3!?W?WUIls(o)BwNB4)$_yQlxjTpXB zgi|^o1f*eknSRVWKf=fX5jW{T;#@G9K4N&8bT>k-OT-vhj-BeG;?y?giW`qZ< ztHt&3Ku5_CM&R#}S0F3XfI+d(9%RU?UIikO-gd>DTklMP#G=p#mzv ztR7R|&ur0bW5=ZeLTo{sySw{_=#b2rd3pFNvWuGXaY9W9<`Nb7K!3;J zv*LW_YU#^My-0Agvpip6hA@|k%)?h2l=m`zDrg%OwC&S;Y(2%-)fmFcl20uZkW(k< z!8x@M`6|zL>du5~=QnnquEy7NLZm?_*PPS^G`+MMBt)Z^IGG|LMvP?W*jXuls+citaXOdHC z=5t9eWrN7LSjPe2)~%D`qmJDEg)yBEhWD+O^eT6k@&dHwEwYTJx%!Sdz9S`SCM4*m zjr68oE$pLYNyDlZ98To`2n|BUNxY3ObmSwpD8(qnbQCrG;8Wbs#I$0AWyOR;OIGF{ zg>C}TBW%qqPBTB6h#zVY`!UL;99nW!kYeT;`HsguB*vv^@^+jGkJEh(eLn9_;%1KR zs!Q=;cTBXCZ`ajF6zQf7Cp5Y@kuA-7!~e`Eprx81A%)~TA#!*5>IOqH`UXbNbioUd zGRN^@;o8ged#;U>g)BsS?45RWY40or3 z;&UAW@b_kGPB2+*VY_+mGEjg%Nc`3d9G_FC`q1U>^MGT(hjv%jiX}FaL>F7j8J*ym z$GxtPg@+6|L630Aqj!#M96kuj)n6>vERvsoFgsk#iB*T_@^HcnOwCRt^Fjf4VGT|; zR^|lKF=a+7;o8YKir_L7!AR+@Wyw{QS8pd>^=~h!8EmES)-7wicLC+7`OA zY1+=4x`eFcygj%w zVuey6=+l+q5U}D*yY>Gv{9C28w@&0Dw7T9DpmjQ#yT0Zbo}+BgQ(^#DBNG$B#}Ly( za;2wHIFgd4{LJSZULFRK8oxyUCjVTa{uUZEAob|}bz_%G(X?kglC?H$Is?mh4MCE*@P z#Fv+_X+*cRK!`+s?d?4tSSSx1>5x+^u@~bq67)Q@Fo6rsMnT#Rl*`0@$E(S27RKu! zM}q-8k?q#~A}BCzxs(5!WvMu-m>S4k_s+nn#XpjI5w7dbC#@e7#k>YD_FZ-(e^3K2 zvKT#?Qv#E5RcIom3fcDLp(?EMuJTepZ8XOsE9&F=TMGpM}_ zv(&wyIu8wulINTu4mdR6I&ad&vr|h;y_Fv+}A86>y_;%toLr`zO@19M+wR(!p^z!qyeLlT;y$R$f?R%fWIfm|j z2}<kPzX&M0?0!mMu!6!KMhavF_A0&TXj&#*t5^AK!s&jWK83f`m zK`eWI`3BgTyXw@dz_p$VxL~iK8h~oS7E6qV;c62Ac2%ve+2NT#^8vsp4hzdS4S^cC z3*jB((06^s>-?EJ^q>Vt>FeeA8kRRroJ$3Q1r9r87i?7PmG_wBS-%OnxerlFhCPvQ1WP`u?fzm?9@KKHj&s*}aCfN3JEVFp ztxjct_#X70Fm34SkAF_3 z^4YElelgEg>vBq`1^|UZg1+G?yrSkLdxBT(U+7ieAxh>6Pupd!6EX;HqEn7P^Xu^5 zk%!dq!M-#RPY50E2xU>zNOkPJtp1u-ID4wGYJI`e5ZuDBr~?(C>~DX_6hgeG5W%i( z`An|Qzv)oor}>h=kP2rHP?O}i0;DSL@oe&64}H=4b3hvbrXG@(?gw1gN_5l)(bCh1 z-My2{6EDYi_e;G7q?g3G;08p%)~^*SJu9V1k}V=gn0)p;l+L&o^p;q^Zz~sQ8DA$gjz;88$0BhwqRxqIc}1vgITSa}=oc?O z!3NpQId*mt7@}#Y!5msTwFWV*DE|eAxvrS{s61-pN%5-OCJg5?ov51lG% z?!g>x>U!>Z242hmdP?Qo$gsz}cRPY-`uqjsd-vw&LKI0xa#ITPztrLJ2@Ekf)DiLX zLAia?on2UV*N1 z2paJQ)s5KBw~d0PvBQ81ld5MsRXP3X+w~zH9ok%QJ;-8r=;MQwI|g$%?OT+{m~VxU z)ldr`2E{8V?e-`~ET9}=0onqJl#eqR@+buNW;}@(mNbsd7V<`b;bEwvW?@|+<`r!J)um)2h@8ZSDN*3ItWAXYe=g5vl8CexR|sjdUFfhvjR8UF~zs&@ei z4e_?HZ}=gomF*^1>b3W4?9XnOeF+V_Y5f{#{S*s}qbg|ujib%B7y@h>33nltY8mt% zset92MiO;Vg^@TDtJJJc*MD5skXhKL#DKiiyD|JVV7MpEo#1VpHqt%gg`Gd*!2*kk zG9V@jF4ggVJYIv2u4f1-NX&^7Yl>HdXU|so)|Xk-RdaYYKnsLG1w_k=Vhpe7Gh$nZ zh24FU*6OL}QpHwJ(TPK8#1E0~RuGac9Tsxmx1J+|v;*r@T<~2k`aE^}n^9H5El=1t z3sHYU=Fa$zEn7&dLGGGH4xcohQ^seI!qAEJk-ifM%IhB%+vr3iBW;Jz6im%HBob603@2pX5zQ;_u3WF z>-Rf$B;?a%FJ!##V6U}q1_o*j_W5&hv3HgTBi@AA0XEqkTme~74~5JNw}4bMV;F33 z-rxt!r%mPz2pui{{B6WmEFAH|8@R7hPVKsBEo?5dG>dYKU=?T#_pM?3?*V9*l1nv| zI`QER>oe(BZBLTf5*beSboN$`EnetFd}&rnw~T^D174An!`yVg6{Q6JT#LRM4wT6V z2W{W`6zg6Hz?07$5t35YU%99;_q9_ddm-da1^Y9?DZCTsyl~vqeI!ynZKZ=oI=>8u z*UD@-l(~0V8$u6~cU0ZIX)KFK1Yvi-4m5mRHG1;H@vd2>>;<{4-rS=;mo;;9z^@_BrSZMKP^ESAEDULJ zkASdOsbbMB6)NXmUDS))mFCw(VLiVrS(Gk>mv|4GG&w(PHMzzk_ux5ysv{2!Rl1v_ zou>+15#W5{3muN{NXr~=iakd?MVss0kS1?sKhixGE;w`=aA?u}^WN8_Ud%PPwmjCh z_&BGECpgTSWBHB7&-LShIyEIc)~Q4pH>Dv~=B^jz5BMXOo39O6DpY2G(j^*6K+9W2 zN`amo$ej4ncdnA$HT&&vTfwlXA+DM#FVc-}rZ~NPbJ6tOD8Nq6T7$*<-lvz-r#D6f z1NPr*>BJ?$xJk$feFf z&`}{^>8C0sPXZusa5vDJAAgWwSG}Iq?%i5?46C3{(nsD#91^ffZB}<=AAFwc5%qI0 z+uU1$bD*u=j=^x`JaF0`dtFmVgGLqf$ydT$r*l}IE;9kPbfrqrLoyTvxZv(}Rz3CK z<(ZB-XyyzX5JcTC>Nkb^0Kz$g>?E_f5HEzT4Cj=w{%O|yDMQn-#3q~Air*0|P1&v& zf%-Wt(0L^f)FngY>ge>|kxk-sa^}8!_J1ll^KdA?_m5jbp^vQDG8rT!%UH4v71_y> zB-@a&7BNghBgEMEEE!`hp@`>+SCoFL>m#r3g?kr(%st=;qsG@6BlLHyKL|IQ&d)!XZ({a@WM%5q^t^ zCVj$B(mss4snz-v9^=C0MMeIWbAzDu3yjv3a}=efY)zjUDlgUdhtpc53_`p6EvD=o z7d}N2ra--4uMpj~7l$v~;|OR=Te8&ZKUSwkju++3*Qo*SX<=s&4_PSKF<6TF-JH{N zSVFGB;E8mBaCr0&>!LHt0uQ!st3H^J@aXKRiDR_5d8~_ye)16)S)Yx`ydXwQEU9nq z(P`=9=_fuXOQkm@T6ezOtF1sq@i=FiWOyW)5?e_97`rs%dpj`o<&Cou4;gEMI@VX!ls`LL%yVL2fzy4YJQVV2RMkk#Hg;yc&bPO29Dy>-g< zG;DyDJQAaX*;Ztr)Pzo39{F6X!HR=RFik=2<4wmECi1_w^JgnuVx?q6(w&EvIDwqV^$Z=e(Z-K^ z4xy64>S2!r&0!Id{LR2O!he4->koSk%Qpsnn4R7`5@Q06Nckwd5yJ{K`^x&p!8LESTG7SjF0vB+=et_>1K|rORW;IJZyJH}$bw*!cBEe^A z0fkeMV2d%8y%SCys^~zjQJ;KQZ#>&Zf2N4qqHkeO&tmVNYC1}-@DxYlGN3c4!EXFQ zWtWA?eSgu#_?Xki&%=Fc?))+Lt3O`-4JXAOb=}JTnuzorQDhc|P0=&x0=j}rpD(L2 zggD$?T1QjOf)=TPYH8}4`~&8XMaIgztnP1qh}Il)6SX$rixHIEI@Fa=Ei~I>jH-o7 zaY@c?R-C2V9fGP}58RSjvWmymW?gjla~ygY%Yq@@tbUfi$uRCumTB<6tA-<^ zopsrApEAJ|#*)5aKx_9=+|z`>C&qmIb|fo^uj9fX<1S?$ZC!BZ-^M$Ky4J1(-~A@b-uFr z;quzn?=!ve{wnhDx)42qsj~u8Ij6y6Z#ML_C6Dr*j7*{@Svsr1c>}59`QNBXt?UKu z4|_`3dU|!U)k@xOv3mc={4WxVh6Vt+}i2oE3RH~uCeZO=}LfR&Hy zODjYllw+6H!rT|70+5;EVh{ypC*u&)}rclcC&;0jORrTnJ z)|2$IxNf=elk;i{kLsgp?I=H@8>bq?<0qqZ9$z=@XNEmvCb@`j$yCy5cJ!qn$t7Up z(q8JIgrfvSlx5clLYT?TQ{Az~{;_52+~ZklX|Z;JX@)QMBI+Rl5nZ<8sk+dv`C1lZ z{%Qte$hxoX_!8f)=Z7;7dX;;3AC}z+*)_cm(+M3pH)%wgK6<@>?dEa4i@pMCZpTWSN5gu#3 z?WL49?NYlP#u50wG=th)cV1vb_DPsB2EW5HnlTL7a2~K7S?{|x@>4IlOn!z|!BRTj z_xT2LvcS-Pb-dbrOqA7CL)V0zpQvgs+(cHJjMBXai? z;`Z9W^Ng_-vu&AYi%zSRy)W<;m`sYd4MN5cU$#J*w^dsXQxj?7Fg4zi?q2Qjtz3l(i`p%|GwUEWq?Rs!c`|bY z{MmF6vu!NLDq#{2QAjPy>kCU`PiDT?me<`@eCh{Eb$_=Uu8FE8#%|^e3A=D|k`U`X zaECotn_CL5r@kB%b1rBJ17(Ixdt@cmXgvnj?jGFHk#DDp7H+O@;H!Io zdQ=N_`|jh7g$;fJxh9>w^SU(&2WoaugWGyK&ibvkWHL2+eq`X(u!uB&{G5fkEF%r$etasMX&hrBQ6Kfi58Dg z){gh^hKcpiNjuH=U)}D`5ful$C%kKy$8oKp4=dM}>?V^>Do)nh;%*_L7aH8ezFSBZ zjt~n`x8G!4>`eROG$dPl@z&)zAP!?)7t(E4TJS6S25LKy4yRuBBCq?6a{DX~ypgF> zeW|O(xS{y~J`z@AhgvU`Ah7PesBmQS49BKaR6b@J_HOsqi&7@GeYWR_cZ)12T3N1= zK!-a=-gCEftPlpIiUM!7_`C6Bc?r2pb0b}oGZ|mGy4+_rp&Lg`i9xP@=OSbDH5*&2 zl%MBQ***^!r8Z1D+r&)+1;0cZP@>nPFVfho*kKOjkLzTB6JQs5Nk#k(TDkY62#~6o z@($Ebrk1Xx1|#?gJz^Cz1IOJpJIr|nbDxk#JXhZSy)m~!H!38zDWR=O3n4T3#!dM! zThVStP8fnNgH_@l%OEYZ<#rp?rH8(YituKEXA2jXCPA@mLiP{-D!n1Z-p|6s9FYTG zNf()yrm!6HLHNRcju!OrsM^c9i{eA26%tQ$@wq8%hu42Bs^Ge&i=P*1q=RU!<0p&% zCbXAaP*CK>F`WiPB`#iT-GoaZU5#WKdRn<*9G;_*VuZMXGda!*(Vb4YQ1hRh$qvyq z&Wp}*@DJ8v?%Rw$@`D=MZqUC zLCbY1yeY515%9x_rN58C*DR92))QpEL5q$QYK|X)q?*V6!ovTfmndhPc$^BU5bTmW z(a=Pd^AGnNc{({pUl<4-23LL{tLFSQ7*~o8ALrD#O)NaQ!WsF2YN`>7w4gOOHcP8Q z2(Xdujnuqy#hi)Rm%{Vr*L@8fAzn{hA<=_#aok28uX_syV8r5|m~#N;xI|TSZd-&? zcUaTjR8^41i!#|%?sb_}?7@>=^iAO56DBp1Z&XcA!RU%kG=79e9v^SywYP@7FrWfJ zL~EY?0a;uu$@?ZP>4YJ>q3+&I_*wNyteoqmU_#Hm+06@i(Kvlf7J;+<8{hTXLb10Dz;+Q|CDOuI4KtK;E|^JuS z`wJc|mH?cA4L%uvu{Ui4~m zW0B?nM<$z_tDN-^09bv=ujvql=Q)5`!5Th5^5A+U6)QdReP>Qc061MAU6^~%SNdyU z>?@wfB{-OV!k|+RDPL?3{AdvWTg3B-?8g50v3$gF{deup|5q{7Hc}~&?@{paNY9Y` zgf@UC@q>wE`jk{k&;7YYEdLj0G?s(@@X~Ag2)r3D1UPItwJ#l^rnTS|bZU5M`fGM# z&c%HV!zZtE$76*-5C@xrV|ELPV^+UmCTx00oardJy*gJEUc}|mriPz+t8@fS*b>=*wwfkUuT;JFH`2)Udh{LS)s`)&>C{oB zFI~wRMXY+UdL>M)+O-5_>*y4)1%aC0GnKzApF42e5_jq(Lz|+wnX`MMEP|Yk0?3#@ z`%`i{OAK{5n~0iklPVg*$5t4b>+z1@js!P>ptdIQTG3VpT@0EssoO3DeWC?`NK+#@ zvns6IAsLvGjGP0amnQ2O*pygv3>3YJM_quLJ~(8ToCFojHVVS7f~FOMWMnZuq(V|r xvae0w8=$Fq1V}5z~{BqRhx7-9%vBn5`u7rNStLY-@LNMi6NV8x`6_#{&P1+I$`G?{aW( z7~ILAjsre;U#hOGtc>%-fD8B#i;tg@f`S74?<yTH%#|(Y%&+ zV}EOg{6Wa~LfK<{l^aYV|O;_y=3n7 zpPGtl84msq5$ew_3SQGqY8=~C=YZ^h1wetBhMR)Ne;72m$4^b~=>U(n7&;0W&CpK37Twf;@kdb=w{~uXBxA~YeEC!=+z`qq&5(Gv8 zl}r%9H}~zg&oAD!&*z|^NbM3Sp&A=5XuT{nTe}s%`nfc9wo{)E)H%bL{I9sVB%t)} zjGq}ruuPkYmZ|KON00yeEey@exOU1 zmX|^_7wU|NAWTb{3)xMfP#uexz}PeovE!C+$7lf!2t9EB2I^}Ok3XJWa7^GpTCw~4 zr+X1_>h6DfqW&NmFR*7rcIlk2*91@&1IypPdxIfWqOdxsr0k~@mscy%7;WSS_WGRZ z-1+sn@OhiZ3#oB?Ns_LTuuzW`q-ghQ<5u4My7>#WB?TVW)3bf`g^HXxyX^VIp1_%X zzojr|ue!dYSKkxnp^{}^PrgDW73<_U-+2}1$UMPL)DGt9_6#Uw3|gLV?-Tg6GKHC+ z*l{_L<>8nYEnzuWQy!Oa=uImvRlids^`lQ5y8>GkY#44uLYc{k^0BmBva{cWGcT0B zN{f7ic;jF2Xqbh(Ru(&(^iN8EBDXiuGG|}f@eYWTM;ygMRVz7YbVd{v7P2nXJN3ZtWC|nmpi|( z(VD$U6L{%1vLZ_1Fr^ShN(M4q3g92MH1_nndiDEp#j4kPwZ+km$*)_PTCmfpcBmrelb98;vr~J0E_%IYpcGmz#m>YconJ#?OmR|Gv*SHJywNv~!~JO+ z7?Z{K!EaA_V@`tvMA#}EHgnUN0^IO;UTUodk`gB{zQH9rEj89JYNlDVV_~RQRjweWqzwt<#6l&~}vN$>6VRtjh zxKfMEBQIJoNa)DRtoaqBYD1jBR8IYr4m%sxdcUN$g0p6GzEgj5Z%e`<2%AF)E% z>hw}6rbi%bJ!EViJh6jfhfKdS%L-Tw^4z#kXxxO^8%fq?9efF2>a^iVv+R?M%JmPQ z4PPmzAGw!7{bB8@0sEW2!fL@JvI+cqC3--k-Wl|xY?jZj&<46jAce7lO5$g)qf2rh z@KDSm@N65S`Z;v=R;v9CSy)ERwEU+Gopt(h{F!BHa@HUkKfaw!PMstU%*Ql|B{MFZ zC^Vj#uD|b{-9Z~ZKWCvzeagK^X|XpFKgHqjP(%OD(lZ_~P+8Ze$Kq55=kF5ic*O~D!SHvu6jrY)g1%u_-pNCq9$S;*kq+t_sii{S$?YgvKZuJnW=|Rfy`9VZSoX+c zGdp@i87!@&5crP8c~e(xv#8(*_>>Gax|t+ zJDlKpI^Bm*aT*);m2&4#od%M%%_M-5gOS@fiLq7p!*>jb%k=|;c|KSZ6k)EMm0X7V zqoOLDX&pnLRINp);LRL{DAh$n>-x;9oZrDY8|+c4>n?BO zZ2fmRUTFCjK6m0P!kq6wCF2i5_tou%A@L&hyJK+qqZ}l=w<4Qrh~J+~ZQiL;;^Zyf zrmy@{n2fz+rpbw&nj#-ySA@~z4;aC6s1YXqz6ri*Y7LE8!f!lLOzZSre~7)pJ50a3 zmyx@!=QQ10Bb2&uru@jLWR@xV+fG14wuxn6`i=iGGT;}%bgxURLtk#aMikQ@Vx`qI zjTfG}->rv%i9CW-DZaxZ5PW`YvzCSYBqt2HkM&^MtCH{DprDMe=+;99e_Xsw>)mk~ z^3yK;m}7--h%2s|h${`C#T;l6A*=C=^iA}DP`Kw`x@Rq1&44)2<{%$(A{gk;j$Tw& z=JJE?;$2F78Uguwm-7*-3!SN*VX3 zsckhRRNnH6_TopWN$033VH2cyIucqm9LcaI$ErUSwc+ftx%Q3n*q8L(N?5gcDWURO zd>RSOnUBr%duu@J6m>yKS^PMhx6-W^f;*>7m>R2|SD+oz;qkbRk^{nI1%-uInzEJT z7VlP+0l0!o3R+Ftu45m6fcv>fmW^5pZ~AZ4yGSnWaM>*A)s6N1@E3m^jwJeuXKxAA znAHFK;Et8Bu(v%D9Q#RL<%Ns_@mP$ z{=$?rw{o7$S`WJb$JF(Q|Ar+ogb;ig0V#l8HAM6C|4kLli02(QB-I8J8MgE+`nQFy z#N;T7YA`cgve(_StW;eyGlv8MWAKo~IG>(Mmz!$z7^%n2nk^oFpRRmU_>4buP@vGSdgFBPf`v&`#}A;$^i`%0DdN{ zppg+WXa(4#bCIs~Pk@05CC$b2{|>|$|LSSf37MJp3`zC0?h`_+to>394@H;0du5M? zRC=!w!WOhPv@kGT%mv5BdPC3Q^oKX(8=8NQ*$~ITM7}@;4hAiyTT9mf`@b`Y?8}|c zx_%x@62&XYSCzC8xpr}nByQN*TJ9N%h?wiq89jfcUpv|vlII4s@rrGpLGtj9$af=E zd77ja;DB524T1C%C;T3ff(>{Se)afL)^fIbBI?oBAEq{eH{1gZg>ucn;zjFj<~COo zw+7nDV8hPZk$aKD&Y{);s&aO)LZ7C6eVba@yxHtwznYxWh7Yq#yMV%LCxI6FH<*@9 z^PwU&z?I)|;Rw$mK~nc@Jot6<@!Pj*J%?5ElJIMf00VP0nYR7sl|=t*{om%@p75Wq zD5Gv8_2v6AKQgJW{!4%8N2!&}3LOxvqND^Qo)C;A-OAS}ymp(_#(pdt#_#4tTvhaI zNnyT;J_0|!)h5wD(KS1%s2+8!~&u0GXp z49(^r0h`UpsbjQq<_PCT=2|gH-=g3vypzRFx5nnm;`LE56GF9WZ@bx2v?Qk!n+olBKK)HKEW4Hrh zDgI;&-Mx09P1I-+vizmaN!LM0!H*3;Q48$`mBr$$)cbs=W}co<3iBPo6F2H)5P8;z z1rOqBK_EAu6#n5#0^_w{?h$$z7W%yVUc&RzhIFMdIyYq41R_5R6?z>#Jv)WmLbT%TDpSz%Nj7!@j zcO*KbSrNDVoHI>bX=kA&O3KMamefb-G>V6nHUDu_=`NBmlS`s7gJQ{|mK_r3C&pTnDPXfi;SRl$wE!o=Mxx%SSBb`o zCO;aW$)j9RoNIc9ix3_=wuM1yNBhN~jD0KO=fT5un+*~xqNxTF1ao3SH}eH_6U1dT z0G23vn^HVDM{P2nSrL6u;gI|~dM0OUED|)qO4hd=(+`g(cWvTa+EbUD%!HY-YUAI18HEb#mPRw8p3))1>3x2l_)* zja2IVw$Vfv2p&aRF=1VBe2#6}hlmXI#rQLN-t=!tT;ix27if47I4bWVdoHq)lzt)_ z$Bm=t#8SNsxZuV^bHU1@>pCmr)h)u$oH}99oKr{8af;eeW%j;)Q+?=g&FwWaLptv4 z(mY!MvFZxvUOw4}E$K+WOreq-dk+V-%Phu0#vw*x>2Cy%Ih)n^0t)xms^l zph?oOA}ngg*>5ZACjxOsV5=4`*Em2=keK44ASSh>JzBqceORco9seCao6br97r~dj z=SNvLNdFGxA$f;f@S{&+9>gij-A-|P?av+figx}&fU&osYF&0R@naBY`y~b07fz(? zU{=+TXMZxqk4@=T|jro%C!8t1443brIXsF%I~%h4zVyU~-CQ z1R}K6Rrninu-x-OLj6xuw?oRvudQ<_g5KG(5|P}M8>uW&4xEmQT_13W+o&>#FS=`F z$Tfjir&=LgWHH2!NS|t+&r!<8d%Rq&Qp%BC2RS%~UvM0F$NyL_aEw+lD{MVwt}RLAhn` z<)3wjpb>s4aIn?>ZrPa=HK;uBWNB;bm{EprsNW8nnPX4Y z_3i6?Hi5P#k1`0Yz)TN?Vjs~L1CN?Qy6vB6Wz2qM?6!feU@hAW2*MRJn~$3RYsosb z%!o4`#Fsce6PUR-jCUa^xuFz%X}p7u%x%$~f$he}#%9Pu@W6UK(_x^i!+h$3hat1q zEr=Po@VCYDXD@g#7(de)>c%ejaOmdyToU5MlsSONDgE52+=!^9;_j#6U>|(i$fa^c z0DCLNoi0Ga7t9BU;=^#FMa1&jultCr0J#8J6x6cOM2z&5H&>?E{v$X4M_dJ0Xpn`q zB0yv~W>s^yP8Dnf$ej)#n%55V5T*fgetAY(22**PPVr)dhB{w40gi~@3hGR@OZ998 z3yoYGkir(y$zYDebhQ_BCJi=eg(f@JFV*S|wT*rGs=odbj-iuL19&KdTvgu8u#J(4 z_WbE;d3T~QiA?ht%w&Mf>AkGFk%wE#`ln?lE)~xm;}0E#ifDp1?YxiL5&)s#LnQH| zd#1?7d2ui59*EJf*Sf>N!9EPGXHHZve*ht zNdv5yfdg`mTw3hYjJ{&@6EnRfFg)o(?fJ7G*IDiPZ4`AHLxXn`@N{0}J{3$E1>*FX zco|P!A#nK;4-zg0)0~zEp5Y=)5xO4dPf>OO58S>#)iMxI-()MC#%d5Rs=ie)v9+FhUUl2l7fbr1`5{OuMD-K*vuUvs?n=Q)@ zytVPO9oLXhLmx{0E0^|!35(d(#EhCZ@{4xL2%tckN$*1i3eA=S%7^X57g^5&%*}ND zotxmpTY@59qbrxLz1Gmkxb$|)&HGjC4H<>3Lp-e*xN340`xL=6MeE8bE>9b#RNzz2 zu|ON3ZySijX#zu+KkNo_<`%gyWF`>ah=~@En}kv_f-yfh6g#n|j)p zNE*d2FD5k5?J$bp6h9;kqT2%2`J0`TUJ$c;OJZfOj!a1w3_IIdzA`tv!~7aEAz1*y z{mbHgjAhsh_o7sxlW26yLuR^OZ@vbiUIFZT;ID+wq`%t^@IeK2tUEKdediHlOkgAZ$!UeWnta*ON330Li)TT zLNL5H&C8>3x@_A3!2ZBp1#&G?Hr8bF{SjIe7rJ4A7n5w=BmDPhI71t`V?5|Z2@T{i zmu{1+yQ8_8HgJ#aBtPEb{5}3OyQ!9UrEHCsHPV<@PI8Do2~*`lSZ(w-d#ZGYIuq^E zJzK#GosFKfz;kBKy1_D^jP`m<&8^as=Lm%%os}0N04M(H=Wfqo?UL!lM_4uqj&7n; zLb#3l%kjRsJOALpw@^z>9fCY=9fJFSHhfFeiFigp!_J1q)CQ}d;r#_}$8jdt((Svt z2<|oI%Oy7kE=^2TKdiH}{wblpy9{u9hKM8=lIke{V`zN@xa(9;n?5{*^Lx%8tWBx1 zM9)2nzl5*!H|j;WxQB7-Ss87=qVm-OytW^4E8{$BgNehS2e$GMf? z3*=9G7#k?2zwv>k3HPa_6Tkxjz6MvMM7~P*GS^Hnb)FHS6PIS_5hg`e^p@oB+1^@- ze0DGo9W2=Ll)HY8x(`f*@;kp%ZT?p!_`LwuTDN$+m|4%A+}++hk1YR*PRg*@U#Vqu zVn?{F9UdVjP1zipc&bMH!^F`08^R5Te>3mW#GLVnr%utKKXio}ye%xvR}Kgv zUy$(pCx$rg>9V}g0Mw7E6r(HfXBsryhGN%5>d=*$hJNO)9Mt(mY_L%h?;!^=pDSer zNA;miqz7*Ap9#nioMMP*>?*%)T-Ws5-FnrHt*9sh)PevN2X>!9A^EZe2>|?z3^TZI zQJMMzat*76S^LBzlkzWeGb$sh?DV_2Hxy+rkNTxN zaNwx%e^iV(dpzZWUq-=lTDiP3l9(QUswZvb$MqBT)he97YBt#28@EBUd@i0~lZ#Xa(VGR&nI)?bKX+fo_-3v;NT3Dxw4aF9;`P7B{{x81 zOAdYUg@Q~~bJjb=@k3sfg={+Ye@kSg;v~3|aYO0L0j=|2h9&WY^fH)Jf-D|!R6Y57 zY|TW5D)1Ih{D{MkaNnqYXjOmO3hKJQOZi#LKQ(@57}2|7<*Wmu`C&JYuKvLx6!=k6S|3tTD{@%&F z<7+{4?$V_5{a^UL-YEMeSPqwQ?!N=#WRGW!k+7i*QL9j2G{F+7cRY92+V$*zcqX|2 zju1^0PU>-f$4|i-y5bcc>Tzysyw@RI^ZzvjVjw(RihLTenSMmigkEAhGc`9OWsbwF z0y9M;ouhJwvC5!9+A*1asd&ByS;0Uh5taZ-7@}*NJPYZoCUR#wKHX1ncXJ>v`M7?Y zobk4*qBRw4PtY@{q#;1B5FMLEmDpxe0oWM%C6;ol5ET@xu{?D8ZeLN-b94rGz$ zOvBs^3kb_;U8*J)e@{?qr-w2F&{o&O#Fx+}#vgU2n-q(q9$|H3U!Ugv*Ti2hV&kS` zUQp5`4XhQxuxHYkxoz&3DN%>{fGG1gmqE?GPON#pA;>QV+lz=m+_eO;Gt}c@q@Ume zr5-SVZLZGa9|5s9kr>Z6b_f!ok%YEw!gu#H$K)T(IR|_#sZw$Y4TeG=5ASv&l$B^@ z@CGazAe9W`$_Nj^5HgZ*8J(SI4_{~IFI0yhf-0b-tiOvY-)Z@ zb)Zazu>d>MXBA?75_UG39$u0^jYPU)ZMI>y+wp89-;C9VX0+XK7H+v#^+@B(zij|q zg;UNrjX#NrFzw#2krk$SkC&)acoX>?-lc-LMOUymZa(KuC%>J2Am`T0 zjqp$wM^@^^Zel-vOHQMvkMAZ9(?RaG^|qG*A&p-IfX@sbJ`CwKTZbcM`40q zr|7=_rlG%EG&6+c$)4x#ZEpg`4vv706Nc87a$tkT5WsQ1Oc$qsH`|YE>;x7{1S!E% z2uC$|sqpOcSJ7_b^(#Dx&%vjb_o*9iQOAs!jTTBp3HEk0q66JNYgt|rVYa$}k1&RZ z8TE+VHxN4ASDKe{hfaC*OnYSskYJKdCBlpU&+j=)-qMmbW=#O|6 zHgf4>AlwNpLym_bA%HzJi+s8!gTMbuY6hn3mB3H#yP7P<1`aXvp_uH{{uoVIpg9`>vES*+x_mg+lfUjNn!ygMGr4LM+|~l) zsOuiBtd$95*p9)h7UP+!07ML(bAGY1NY3wgJc_>inxu2mc~70vG$rT%c3B9D+OVK# z-pdw=S*rO#fqmf1>G>PLnCEa$errK4;xPajEB0mu3M{C2top}>3HLHJ?u=auL3Q*;g`f#Pt^DywaictOTn0k zC^i1e3CX$`_zi$_&H2|;qhs#PvSBY;dUs8{!{=heB>Bf3T zKhXGjfY-_`(Fxaxil8sIo0PrZ#K5RdUu92aS>3ijx8M&<6Td5h%JVpvKn?aZ!NT2o zH<}41b=2dAw6ZRU0A%=FdsBc%Pp>W3duXSc2=#HvR# zR8gsZdLc!-w>PrC))yM}Q~m)Z9`ARSQ$;CXR{nM%cMZrPm-H=Bl%+tBM&VRszQH*L zZ&FYe$sm}0zR&@1_5|{as#uMZP+tnAdz<7`QWhWTdqzrS-n&J11mO&3AYn?%LgDb5 zr5qNi7i9U=eg%ZV7A3~zRn_=z0$dO#C^UW8`Sa_q4I;1wxWc&7+j780ORU#^RaMjy zUifgZFT`9#ZF`~Ixs^A#g;Pcc)k%t$Cm(`fNL29@`0-K*^qxL&AuV%cbN7WhNzAdr z=bwe=QDF)qFk$gbxeZT&ODyiEM^O8Pan{w?ElDme-=vnG-LaA$%NqOs)5bEoEk#R; z3-BR6Q17BgWDq}vgfOcc+A~doFlu_U-LulYdfoU{^_lL9!nWRF=hC*9q*j({kUXqJ zQP|q@7bB-49(a(&!XseY^D+qs_!3bA7QB5MR)`Ldb8M9Klx4Gk?R42UWHeSJL5c54 ziXtx-Z|?ly_=C0BS;GpweBdqz-!U*@ZS_oVM4b;a!vswlPqw9(6!`+f)9B`BHi_{r z`b&A&a4pUA57$)v6h$-d0}WBBnqD$v5<)MxMmj}z)VD(m`8z32)ePfw4>=aU+to^` z#~Q~Qh9bpbutPfIA33I-b4)t&x2%c=boWsOGw6(5T2ZJ_&)*Zv>TZ91E7@R{Z65Ws zMDzG^KXz^2{BM`!iUYjddG&B@?43ECFw3Kunko5_>CE$==AYG7l)0Wly>UKiP#5g! zvOrgGRQu~~qtK@b@vQHI<8pj-zn4d#=UMYy48&y&sa`>|#4|79m=FXp>W#DTyZ>ZB z=up!#Ovt4l!4CHg%I<)1wS5NFKo-uhR}4=h4L27KP6ubw7#x zZi)%m<73q+a_VRg?u|u`z3{d3E3@*&KPZ7Oz6{?-b;x}X`4Dl2!+TC$F3T=DxUyvD zF)>5S>*{kV6BRHzIQ+K`)i*nq@j7?&v2y7>b0zgykZUsGcoBd&3PEao?Dqugimh#J zRjl7{(mtIAxc!QRn~y&lW7Xe6Z0Jzg@K%sCK8z%siK?w23X`WorhA8HUn4bVnOdzw zh#x9BgM^l|!~l$N2YMg1%Wq5G~ztowOhzdbGewyXylk%a9xZJuQt zHU65TC4PV2A1Lwrt?95ar$8p+0SqL zNe9A=v^`E=3_Gt-R5%sM-;!4%kigWObqVR-I+yT-)mCh^;-AI$48ATqxFnXN7=zqrxBOG`DT~jN=@on%PWTI#M-Fg((?uKlHXXh=@07ItnOAn=<WzSSGbl~q>QAOIW)hRwqJ%ENLj2^@lmFhgYSWhY+w+d@KsJ|09FoLcQ<{D?&vyMBFkx_tJqDfk0M2EMsrWpHi6< z4%DB+mNrfi3e-E}@#{}Drl@(BsXv=NVdHNswU#LO_RHeXERVd`rV!_yvxge6#!|s!?bz@v+A)|sK}quI);- z4@^LXtoMT9v*1Jg%xk9y*uGd7JuWGrXQedR835_h85oTZEKHm{84G>A?hRyS}F%)bDSb~-g9bYb|6T# zYC2fPyAf|Iae28lH1y-#w;WI(;R_6KPEs)*4i1@k7E@Q;1?;8Kh|yEaA15A1c?eNQ zH7*i<{(FL3AmRRv(5aMib`Fjf88)RP%-1E%N=U#iQ=J@yHMQ?`y0HG^9RV7r+>|-3 zRhcdm#-Rv{Wtf|cFo(e&a*pAQY0UNUck>5%9SILUx7#fR|g zP57=-C1&yLaj_XmJvK037L&okvy-5>sBQLNdHk`#acRX9-h^**iUNpG6Vx_jj32}OIRHQNkN^|1@V3{&<~$NH{qr^VhO&k}VS-@@^W%yW zZ0WAVNCRRG<yT>?2e5p&)PWe$#?kTo+zC5ven4(Z9e;IyAECm?s1 z_#%xQ0|O(TrE_xnT)x)9M4gFvHDxXxXtipvvoB%Ed`@M~!qKe0EQbR(w~B2~fm--V z-}H3EBGP^Z324FFUk2|_FmgRP=BfNc4V)jf1n==X0gI0^M6uF(&|T8} zb_zH5`)a&W1#y@^M*pS7zb~yo!jA~usozvGc6>{$F}#(KJ<@~){((~q5c#%>1#~h5 z;ZuprL-R;OpTTo$6tUsaNaY^GQ1DlpCVe2LmxJ@W$iYV+w?O@sEU67|(~1v$gbw!M z?egv~0*hi}soyc&;O93el6m4igT`l(FHi{DU#O?<^6}p@L<`_0-xU|YkzP!!EvHY5 zRsWh@C+fb{*FGZ1oWA7EgN97RY2HgpEB?Rj5=#MKX~j=X*sdwdMk>npq8Tj0HCRk1 zsvgM}YedO#sC4xFZEm=`GZ3zmVhOW>$`g3KIrcu)n@HyEP9#OV47aczCjt7lJ_A+t z1co(4%ysmdOmZHi7huf-op}TSzg&}BEVeZwhEliG7kDs`D1R##E;pn!9{q*QCtV)p z1mG_Cfo^IVM8Rx$Ay>SsWQpIc%B)8+$bzJTQxrXA?0BSD%6*j$2sGpG^(P^yJOV7l zU&O>_hypy~smsWH3-kv00*nLp(|T!NQ&hat5*VkPFMeWTGLQ9J1VXAI%mGjrK##P$ z8p)AVK-5)*)rMhx3_*8lEOk^`@txHXS%gMhF%uH#Lb!WnFji2 zU&p14>tcORP`SY;*1BDPH@lAA{q#xSe5H9F|BE}!u^1*m6@b?w^l?4ySl!Kwx9>nC zb?yTekeNNskbNe2J+*c5H0HT1#c5|Lb*2#mAfmy76>OU|M=8{Z^u21pSX-@aUab67 zAT8z`60-%Wv5Vr7X}^j~E9whyk;y}{6pSybLQ-Qhs3D;Z4~X;vv-e?mM`FEUucm*` zD&GOJ6I`z?Lral^0D%X1(lTFSm6eTc`y}^FCDFR1c1$k@sXOm5$7AZDq%;Gv!0ICf zVTgN%nlByt{0sp%{E8bAh;{n%TZ={KC}+;<153rFor!p@k6CHO_fR8W5ZUkr)@+~? zpFf(5ud zti$ts{DuZbDod5vC-}5|#q5afDJ;s=k2ZaQ%V)dq@XsE#`S>duUfQ8vEy?$eDQ8f> z#euom%oNdQd2ZRF-0vO@ZDm1s4T0*Nx(}ozz0SJJO$L1*I~yF{&4hRa7uJxFrT}`W zMxDtIs1NXw7u0KMXgRD)IV!Mw7K0h>+5(5cw2*8}f}eaw5ZELl1X@%+ug2fcjTK(m zBqj@;!8QT?0gUI+d%_jLY*}L>79;ZWN9U2}HB>kWOcfDnJ$@fCf)?QW)2LEKKScBhSw|?s505ON+n`gSd;&@s1za=31fLVc-PTXK@IIn&nScC#L zedJr9Dm=V{_!{+S5~P~v;m%&m&oTX#c&E>&GmRxrMWlp-betmK*p04!wZ4sjE z++)#;7tr#%F7b>l1R|7h9OGoRKOCPG(Y18zx%Y*K9T0Cld?-Fw$=(R8OQ@9AH#TpS zr`T=N^m`J4w0<6hs%0K`Uwz4VoPMuZQe|Us*@F#xm*|xO+$G#0nOQ0G?OM`8_-d|23NvWzVpe<^5{DLo>EC1AqW3=oD;ODrkJPh+sceb^)bR`W<%- zW@&IXZGt2RGIgeoe3mM8(6(=3WveJg$Yy;oNE<(sEdrA$zDDT)I0S=@Yc?h8jGS7a zeS<`V6*ls^g)1R~4KZj1?Zh^J5$hDEZYx$|6IQ0^DXVtkIC=+^Iewi5*Hf1fPo)Lc zw?+FPLZ?RE$=sO_#Gy`fn`!Nj>D%lb1RnGz+T8IFq^Oi*I^2AF>b@QRt_~Mr=N(jy zJv6E%Bam7R_T?bQ8OLh$pFA1g zo-!n!C3{t86JrAD*ASJ#qPSNo8u~5noe5(uVXgwEWM@QB#e2w&GZOd5=6g6xZhft< z;hC`+ozr`hY|CAu&vw5Wf`V0>UoW3~lcNhiFz2ez0QQqUrtSeca2EJUyc98CD{{X_ zeX5)ST8n2A^M>AsAK?R32z;7^PdehLTyMZ+Sf8VkwVv6Gzk_$?09w_eZ(!uQbb}vE zM_R) zojF>7IG?m6C zXvksfzP%Lui*R8YFhfmIX@&Xg2v#=`9+<6nLMLKNbWy4^YZX_iwmI?s@R_bcEj!0X zZ{k9s=QzyHQ^hWvhCX52SKIDpHKM_W%0Uo7FC$kBmTD$4+-iO(E9A()*25I^lFk)Q z@?Jz7)H%QQC_(*P#AYuwG2mJ*bX?xBDg2!x8cpQpgqBV<0K6I{yc4{Or=7@+XO6QJ zIqw+IJ|Tl(TlHk?GP-*%*BP534RO@{ky~ja9ke0-kgdKi;7jneghL*H?Jiep*Z+`++3we!JaxK==~RaZ zr3M3{*HDrzNeI7+cwOfm)b%i>ZEK5iN<#XE9JL#qZ>Fr7#lVt+Iw`NU>&eXdDf``O zd2!~)6~us4hs!O4`YFe0OoT=FO<$TsIIpNI6HGE}6WvbCg`d>^Es?F^xE*M1LfZR@ zoGLuX{2}WT4*`7{z96gVZ|h%0^T+U}NdcRGixIHoMa7v+kiv2N31Ht{6y?dV6;rD4 zmYn(&T;ywGGiV{Z57C%qH|gsji#6@XZm??zlgQ~^ofY=Csv7#-UWH`4B$rs=@Kgzk=KzELR^-*rlPC!6`B`~<210(Nkp z>F`8q4OGuTqo$$X=zhcb8T$=GU?-J3i#r{3IZy6W?_@o(G*&isLMtdGyDn%ak|wob z&mm0iN1*DZwaYa=VzK(>1P{4PaXM-yZB**a&;OAJ?S|ck{ftMG47bdx$!AAzNDLK2 zZWkH#2~h$KE@|ANPW?k%Kk)d|k0~ap_(|4|?p8Jz#L}L}%ek&14qF{yBOWjR!!Sg} zy-FBU^l2Zp&|#U#ey^?ughQJ_|JKqDfs{=zhoSjcd0A!>z?SP;d&FNH||a zU3&oDfI9#GjSClgt`y!Y9#2Mk_VbTDa>n>n$i8WHbV=ootdjw4Va%H`?0K7RKG!Rog@1}J#+Buo5 zQ!4nk1+-$JB`g{91Z97BGTbxB5)WF#_-4Fdl(6K5B1=LbR{HQvS z00crfjuj{n7MbMtOo?k6XuiZjH-lv#JFAmZ=3qqIm|@8+g~=+2P{d1XePPRZhQ)PG zbBci!m%Wz~ov5uF)^F)m;NU^WYf7{I$rZcHFCT1F05r`}xpmDU0?K~FCy(uq2Rrss zWGMt=Kw$|gA>q8kXV1+Z)1)>nBEF+NS?Q1qS^a@vfijrl=tvBDbM#ttZ%%mvC`udC zv@WsU$;a|Y zo4e|fzZ@VGWVvrpnL#~})gp_Y1+RolKhelGI!VbcNHN7M>Gty?-hHkwI({D2V?P=E zT7Y97(yhDMTB*2>?g0E@Kim9K{+CNK8$hsJI>6hM4rs%OP{dkYPxw*zicJ&LwY-@0 zjJ3#)3DJs}gidX88&eD{NJlFKh&S)OF%8TA685E5F;F8{G%f>+TV+D-G$0CZatoBh zP=`;!TE3ULSa~MRim?`WUigOIC8OM4rIZ0WXTo&7{Ur3LKzqJq*J~&A`kUTFqE;!; zb(NOcQj55v?1DeON%*h>m}Oz=P-#6nn0gd)pVmcPJyG(T-8YWHI+%A*|nL+2nbK;mVqEX&eLZFq2Ulwl~yT4?tN|t?1LiYw=*vF1VdCL6V(Ib!S8L zIYnIb!r6=d-;1+kYwzBsAX#npc~3X$unmZ{O4osHEWta~Dt>A-DN2LHGz0kGaJ=UhwLr#uK>5TeQtT>FaFa)w_L#>6nYO)rT$2`EKsS!&fw#@mzlBA$P@? zu^IbN0ShvJESQt%PUfJmX{{}br@B1e%CI{t zWH(}~g0jFC);~Q_9H=a~u*SJ4ov4V#Duer(goI?ed!)=k{*vY%~= zz46>}(8A5Z(s`7WB|bLe+LrGzgKFMDO)YOXsx4(pD#q7R1|(d6{NL7`Ev87}I?0=p zN+Tw%7^qplvV__Uq~r|qa%manN*inFOLnc-ZGpD1a#CyQuTiO zNmiK*e+1xXZ3M^LAC`GzXb@}gNzBq0t+HKH-;fR|Te>6{8D^S0ho*E03erQ9eFU~|y&sL{l z)}qI&P2s6i{HHDI6Gly33RNO245+~cAMGBNwv@(~<%B6*UhfVXX`wjCc8&#h~E|XRe`2vehSM>lZ6y4tQ zk<;Fd_;NTgl8ElJdVtKV(B*~4<;hzM{|-D$(SnvpP9fG+Dl#1)2-X_`kwZ{?F}+f# zu5djala7##>UUtk#}c+`a+{3#t*4asK%UkI8!fSB0DHUlR{7(^YhN@wujks!&cUz; z*umPuodujJvV(3Ihg(gn0=5`98LXqvWA%~WpKJQt7t&MCh_9dPUB;dxIJ1-btAo~J zMCs0W+|7#psSmOlO!Wjg3^atAA!>z#5~tNeL3^O@){!4hJ-;EW6wsRh#NxM`&Y+%t z)$+2Cw#iQl#u`^R;*wa^3SYr@CXvhJlV?5Sg|66NlG9@MgB0<^T)0kf05bI{xP*pZ zMrxk0Q~s@!5L&}t_`44&7aZ|bTcV)(7Dtn<%iv#*n$v{S)Nq z_NHH!ut!*-h{`t+6wIWR@W2k)H>O_vGr7%+q^Ll@4D9<76#yX>P@3YXe6TXoEW3AL zr)#-e$6|_ySUJ0B3)?hfN;%kh>F~on)ln?>>Hefw-pZF3z&m@l?buTW&HKj@@8Sxz za#t*X_c7op_)WO)l0M)i;*?o8xva(fCN3%%C#tiTKq){%t&E?p*i~|H?$x`>itWAa zZrNGPsaO2GrTD^nCrS4YlfB`KVFW4QmC2)}<~al6$2(ilruCV6-*YyTgKrbldA* z%d;Ki5j!jTa%R(MrQm7`*h>b7-wg*KuX@(@hTODbnJHd{q)Lrl-rn^N=2Ouq#`knK z4r(X^lwMSX6cck72g)5Ka2KL^(g7%xndH*TQsx3Xx<9h{vNvt_!}hTca0}o`oVE0T zBEEmF`kjWq7*HMP#N`_7d3agb`5O-R8#=4+4WbV>>8@R)z-*(gb2EHRkXGh7BIA*i zCR~x}#H)DfM)gPL)+MKk?^d*z<4r9!#T(GApsF88IJxR>?KSr(6 z!q_9bFjUAO5)s9aiR{_NGRDZyl0A`q&HgJRO9nG$-aC5#=Y7w4xAUGe=gge%_dEA{ z@AEwOx%WQL=kwr6J#&eqwrO*t+b|$n&CW*W%o$bdOe7HPpV$9_fD^4U(u#SQ=?Txf zt_ZnDM<3Xp49H~%i%9rAB1cm>NltG-0LobV@m`zxR$NQHlJ1rtDyST6`xJiM$)H$U z3MzIMZ3@uYtRNhFC`!+b{RB_WGajoc78LVD6uw=|Wr=x7W_wrAJA%?Fc3(pGL7w$g zxR%)NK8fnuV>n7FuBIw)UtYRIUv*;051qotAMxTX4|0g}B!Bn|2cdn|ORVJZ=*bcB z+9!J*5G^*YfTa$JYjS3t$MBB_yfj{nKm9d3XDGvzg(O=zuU9{kr0&%W#NWgkgTp?E zayp%k4-I22g3DS#CJw)eqlCg5^pCC=e47WTjPB<70)Mif>R?CIC^te(pgIL;Omf`7$x70o@g=?>Cm>32mQFWo5*2f zDaGy&215P|e@z7ZGe-4tTK!!4W3NrV5!Zf8yv#L9=ZN8NdBg1WpAT?&qU+D%3|WjW zhO&7>a#S6_e{2SGDYo_f5VoqC(L?DXS%)P;m6ee423&blhXpkI~V84ff?+q&_Bf#;Onw z-8xh%)Qn9Naj#KL^weFRBFf*NK3OZilMAADs&-ISt?D=w%>cDM{iHTL*>F^AgcP3m zVVdh|=(mE0ccqQ@wszzV9)_B0cMmcDab}R~EP>^2qj30s*YwD$v0Y*d&MC5Ooan8G zAr|{(PZDW%(l;)rgy38R+4br1_YIMsy|`v_W^2n3mvxc~^DZ73Cv*DhpDyZKt87+_ z>UZ(HCbpQ~1VJWYX&m9^$>S@5q-CVV0_q-duj5vmZN=rX>b?|0i}@z-{DkAbLh|y2 zUdI3w@(~bdW}LxQn)W znk395l#I7e1Sw|*ztBY0oWh6_@LH7`3nz4?Ks1F3t@?|4+?;Uz>T=DFve|q(F<4ujf^dU=?BhjkV2zIWC7Z&D2T!iLD zf5q}&H1^%Fa3dI;yTBPK8NnU7Jz$lMJr{0jVyfGHk-3*Cn%O$}^%QcSxF0iJLUwM< zeBe+ktY*q->5r! za-~(W99tz`fj4VAdeP>+7B1Mj!F=X-uD`%bSzumZ0F1}r5976dIw z0+X7wszM$tS*Qi5L-I#^`xYc9`c{Za2SGc|rlSbhx7Clns}A8!EYo8#HRkE_WQB`6 zKb#sLt_n2x(AsFxV-aZ66Tb)Ta}Ni%FV4-4nBr7b%+>7spicHZ4YCMp<)r*3of(0d zGSadZG7AKLP!3vRb)+k#=~DZeJ&O@)?(pb*Ght`WKO@%=d;YhQzsMp*BfPl9iuNPe ztoLBFw&6Mq3-|ai&1_-G-@CGeg7BTaa@(&en*nn zK9By#T+w*+-(?7^3}`ZnJhXB;8&^JQV;4eNx*|5i4v#C-T&cqI2Q_8twjRY%JzTs> zDz%P_f=Es%y{Ms~r)I~uCIJm;OB+a7ITp0;Tc6c{*UDYFq*DkoT);(31@k{F^)x^{ z=Q@PD?@?nh@`iHx`O^`4Wx)@>v;mU9%wFHSlVwkS zo%(rq+wIha5Z(~N5=vtuS0QL2yC|+CX!z-IET2~%D^JZCx`zEbqe7#I?&;|iAKVVgt zxE=-YdTWaB{{FaXxjwmF>@q!a-(dNZ$KHnQ>rDaSKSC4zr%-AB&V*)xzMIZX(`o$= zQbS;0WzX8dLP`cu8b02(Nbw3NeKMHc#&T0GbY!09zjc=;`?@b(50mfq9{=mR2D#zI zX>L$|hOxPMF_`|9i6Wofo>;%J+hioa`n<;;Q*o-S;1X$CBpIBo7LbuN^o?s zFWh;!fW*01AlSk5dTMr@x)RaYO|w_ul+u(=(aDrzb zZOg(xS1>6RY@`_Vsf@ep74j#uHkUUuTm42GVYXYVt<-9@@~Zv}Y+mT6Q{6cMDiTk*XQiQyNWIi;Qp<}eQOlts zUIRC^iYh&DWf|m?)GPd|jpSE`&}XgHH+9p%&xxqk%=!$as|g< z{J`bbt(Es;#w)4YVL?>rL|VTK{aqlUh1i|FI4h z6h2_+8S0MAjJRQg?1F4ASGu#OdX8sC372&|bw@=VSFZ^Q6)S1liT@hPJ()^ujTY~5 z+9xXH3|PP3n0~rPiFZ=KXzNIY{I*I{`O2^KmYThA!!<8|M^`#Yexl*{=8$z_V112q zU6|1@KW4yx_g=lcwPD&j+t0gkXZ`l>ab&`6_1kp~R(3{>)ll_swvU27CNEUket5I( zpU}g;SnbfNXKhLCuWz7y1lAj_0N=UG`BoPTilQ;P50FC}$zA&tgE_#c1;{zc|G?@T za1Rl4|Du|N#fxE_^b#3R-_Ara3?L80%pXwRVnF5`T>MvsyMso#xqq|C4}Tav``!zpuY9PlTU^^mlVoIv3 zmq72aWIN#o=BC)mV9SEI*8Fn}00a_~kHm8zWX)ntjbMj=Ev70u!${0$UGFyWcJVahR zc6wL5znhJj_fg^_8p$Tj$WKDsFb2uQ006%Yl?U}qa=*mhkg|BU7h;jvhb>IMt#E7L zry_44hWQ6%kcXBWi}jRHf~KEN_fuD9e(kutRZb^%Qut*=lCE*Gdy_BU@W9-tMGx{RNF)u( zjZWi1=yojZ*x-pmcU^t6JZ?9W^BrFznx9I1L7t_ASP($R>L{mSMV^a_KfA-n(r${smV`-&PKLr_?j$~ z$Mtdgv-r5HCY(A_6vzCVXP%sZR#pAJqa$MKag4Jq5*Tdrxu|Kpn~9XuI|YKfkH}#U z38lM|tOn@!lK8Q=!BL!7vDwCS#^5@m)1N`E;c6WQ+UH2Attj4pK5^qa8SajpIX&De z(B^!a^w^}vw%%G5Fw}|pzZm;&QyjK|KBAn|6Ak#pZi56CQ2B zfqHVp5oPqjOhTBQuY7cchzt-iE4I|ZssNH`lG}xIH1&70Y_KAToA9`_#7hIfJT)!) zcqJ3x=ixZsw6iuZ(IzPYKo|h=ArqI8ef1j=kd+H0fcL?ksy5|Sd44NSB_FXKgWzEF z=Ii^qZJp{XoJu*bAZWxMn7@b%Zl%)82Gt1Vzj%_4tl~^L ze92V~C%B;*m;Bq?n@-vtKC@E5vRFrPQvNIBWY=OX{u)l8fQQjbx0VddDBvoHX93Qt z@tv(5L!CnOwL*a;>4Ud^nco4(GY0t6OmZOW(So@V<{x!jpvEfGbzLK!(m(AU F{~M&HxyS$j From ab40eb702561ddf79d8fc3490337768d3526e72b Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Thu, 15 Feb 2024 13:16:50 -0800 Subject: [PATCH 09/11] update results --- docs/tutorials/signal_unc.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/tutorials/signal_unc.md b/docs/tutorials/signal_unc.md index ef0ec2c..68224eb 100644 --- a/docs/tutorials/signal_unc.md +++ b/docs/tutorials/signal_unc.md @@ -61,8 +61,8 @@ print(f"POI upper limit: {statistical_model_sigunc.poi_upper_limit():.5f}") ``` ```python -1 - CLs: 0.96419 -POI upper limit: 0.89459 +1 - CLs: 0.96607 +POI upper limit: 0.88808 ``` Let us also check the $\chi^2$ distribution with respect to POI which we expect the distribution should get wider with signal uncertainties. For this comparison we first need to define the model without signal uncertainties: From 3a04654ae8ca124ffdf04a26bc9a5892b0892e00 Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Thu, 15 Feb 2024 13:17:03 -0800 Subject: [PATCH 10/11] further fixes --- .../default_pdf/uncertainty_synthesizer.py | 20 ++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/spey/backends/default_pdf/uncertainty_synthesizer.py b/src/spey/backends/default_pdf/uncertainty_synthesizer.py index 2f59ba7..99520a9 100644 --- a/src/spey/backends/default_pdf/uncertainty_synthesizer.py +++ b/src/spey/backends/default_pdf/uncertainty_synthesizer.py @@ -12,6 +12,17 @@ def constraint_from_corr( correlation_matrix: List[List[float]], size: int, domain: slice ) -> List[Dict[Text, Any]]: + """ + Derive constraints from inputs + + Args: + correlation_matrix (``List[List[float]]``): correlation matrix + size (``int``): size of the signal vector + domain (``slice``): domain of the nuisances + + Returns: + ``List[Dict[Text, Any]]``: + """ if correlation_matrix is not None: corr = np.array(correlation_matrix) constraint_term = [ @@ -67,7 +78,7 @@ def signal_uncertainty_synthesizer( absolute_uncertainties = np.array(absolute_uncertainties) def lam_signal(pars: np.ndarray) -> np.ndarray: - return pars[0] * absolute_uncertainties * pars[domain] + return absolute_uncertainties * pars[domain] constraint_term = constraint_from_corr( correlation_matrix, len(absolute_uncertainties), domain @@ -86,12 +97,12 @@ def effective_sigma(pars: np.ndarray) -> np.ndarray: """Compute effective sigma""" return np.sqrt( sigma_plus * sigma_minus - + (sigma_plus - sigma_minus) * (pars - signal_yields) + + (sigma_plus - sigma_minus) * (pars[domain] - signal_yields) ) def lam_signal(pars: np.ndarray) -> np.ndarray: """Compute lambda for Main model""" - return pars[0] * effective_sigma(pars[domain]) * pars[domain] + return effective_sigma(pars) * pars[domain] constraint_term = constraint_from_corr( correlation_matrix, len(sigma_plus), domain @@ -120,8 +131,7 @@ def lam_signal(pars: np.ndarray) -> np.ndarray: expectation value of the poisson distribution with respect to nuisance parameters. """ - nI = A + B * pars[domain] + C * np.square(pars[domain]) - return pars[0] * nI + return A + B * pars[domain] + C * np.square(pars[domain]) constraint_term = [ { From 45b925fbf1c26e7e22e4a2672e87f54edcf774d2 Mon Sep 17 00:00:00 2001 From: "Jack Y. Araz" Date: Thu, 15 Feb 2024 13:18:43 -0800 Subject: [PATCH 11/11] update plotting code --- docs/tutorials/signal_unc.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/tutorials/signal_unc.md b/docs/tutorials/signal_unc.md index 68224eb..6cc3bea 100644 --- a/docs/tutorials/signal_unc.md +++ b/docs/tutorials/signal_unc.md @@ -37,7 +37,7 @@ import matplotlib.pyplot as plt We can add uncorrelated signal uncertainties just like in uncorrelated background case $$ - \mathcal{L}(\mu, \theta) = \prod_{i\in{\rm bins}}{\rm Poiss}(n_i|\mu (n^{(s)}_i + \theta_i^{(s)}\sigma^{(s)}_i) + n^{(b)}_i + \theta_i^{(b)}\sigma^{(b)}_i) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta_j^{(b)}|0, 1) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta^{(s)}_j|0, 1)\ , + \mathcal{L}(\mu, \theta) = \prod_{i\in{\rm bins}}{\rm Poiss}(n_i|\mu n^{(s)}_i + \theta_i^{(s)}\sigma^{(s)}_i + n^{(b)}_i + \theta_i^{(b)}\sigma^{(b)}_i) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta_j^{(b)}|0, 1) \cdot \prod_{j\in{\rm nui}}\mathcal{N}(\theta^{(s)}_j|0, 1)\ , $$ where $(s)$ superscript indicates signal and $(b)$ indicates background. @@ -81,8 +81,8 @@ Using ``statistical_model`` and ``statistical_model_sigunc`` we can compute the ```{code-cell} ipython3 :tags: [hide-cell] poi = np.linspace(-3,2,20) -plt.plot(poi, [statistical_model.chi2(p) for p in poi], color="b", label="no signal uncertainties") -plt.plot(poi, [statistical_model_sigunc.chi2(p) for p in poi], color="r", label="with signal uncertainties") +plt.plot(poi, [statistical_model.chi2(p, allow_negative_signal=True) for p in poi], color="b", label="no signal uncertainties") +plt.plot(poi, [statistical_model_sigunc.chi2(p, allow_negative_signal=True) for p in poi], color="r", label="with signal uncertainties") plt.legend() plt.xlabel("$\mu$") plt.ylabel("$\chi^2(\mu)$")