From 764b7c703186f0bf5ad77a137f79870bd66138dc Mon Sep 17 00:00:00 2001
From: Michael <m.langguth@fz-juelich.de>
Date: Mon, 18 Jul 2022 13:40:40 +0200
Subject: [PATCH] Integrate Jupyter Notebooks for evaluating IFS forecasts and
evaluation on shared subdomain.
---
Jupyter_Notebooks/first_cond_quantile.png | Bin 36035 -> 0 bytes
Jupyter_Notebooks/get_era5_forecasts.sh | 61 ++
.../get_era5_metrics_netcdf.ipynb | 394 ++++++++++
Jupyter_Notebooks/get_metrics_joint_dom.ipynb | 338 ++++++++
.../juwels_juwelsbooster_compare_old.ipynb | 684 -----------------
Jupyter_Notebooks/performance_check.ipynb | 724 ------------------
6 files changed, 793 insertions(+), 1408 deletions(-)
delete mode 100644 Jupyter_Notebooks/first_cond_quantile.png
create mode 100644 Jupyter_Notebooks/get_era5_forecasts.sh
create mode 100644 Jupyter_Notebooks/get_era5_metrics_netcdf.ipynb
create mode 100644 Jupyter_Notebooks/get_metrics_joint_dom.ipynb
delete mode 100644 Jupyter_Notebooks/juwels_juwelsbooster_compare_old.ipynb
delete mode 100644 Jupyter_Notebooks/performance_check.ipynb
diff --git a/Jupyter_Notebooks/first_cond_quantile.png b/Jupyter_Notebooks/first_cond_quantile.png
deleted file mode 100644
index 6ff3a7a8a081c4a874d2e2a8c8d3d0d2e47d1fb5..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001
literal 36035
zcmeAS@N?(olHy`uVBq!ia0y~yVCrCCVBEmL#=yX!)iLEc0|SF)iEBhjaDG}zd16s2
zgKuI<K~8>2PG*uqS!z*nW`3Trp0S>xjzUIBNkOrdzJ4xTfnI)5y8hj}e?BuXFmM)l
zL>4nJa0`PlBg3pY5)2GI2A(dCAr*7p+%2uSdiC+~kJlx8dZ*~UoFt*3DXVn1x8;kQ
zm${d>W^c)2A*CjbCax6?mkf99o;|Vq^Y6!fzmr=YgqP{(ynp_<P=5FK$<Jd1?^(|m
zJ{Hy8(b3UyLRs#eUIU1=Cr(7b)y2g{g@@HeP*6~i)0r{Q#l^+tP=JGoprGJHgHdVP
z2aEEm-P{6#Dn2tVEc2Xv?CwzkLBWaP3AeZ9hOgP$c_)^yJjkJj_0kg0!`IeE-;L!F
z6cnsv?J-op8_QOHSvidq_jOERyeM(`v_x{Rp|ITTv>wC1y@t<w48QjnzV0!+c_-2j
zWLSAe3FAeHV1^H@W%e6aTVG(!;D5rb$GBf2IcfK{79}MmrFWt$nHCv7(|w@5K(gTP
z#<y%QICwT+)7yCM9itTU#>QWH8mmBw>$c0P1`Z$AeOs@st4*qoQqbJ+>mpnBTiXZx
z3nskoQdd$^dQ>+tZKZ9pkGBHzg%6UuC%66E^W(_3?gO{a@0_RK=YDmG*VI>{Q-fSw
zT#o!za&<5{wuiZfO=ROWwtZdMQQ6BMWlU7bzqjXGXiTt+i_5?HTc@$}%N4x7rki(v
z-`Tgfx6gf_e6-8;<D;W@V~=oX2eG&;P?#Wa<MoR?4hNA1GYprBe#qETB=q1)@^e*2
zy&bD2C@Cp@`hTon-oKqscHYkO63zSf+uz)jd1&GXfyiG9b^Iyp_iEp46^#)1^-rKr
zA>!9X&SMh!3p6%dY<(&sDA;*L&aUP~&HYuO+D7T;e41D|J!YHfW?ouy@bF>5*dD_r
z9up769_{&{<m#~H+Pw|e-c|l#y~DV!wQk!1i3RP3%i0W&vCrUTVTj+o_JNX;(nObp
z1ck%x{OZnaJTW_qSQ-5N{lC4xe}1B}yGF<g0d>DQ3Mwir4fE@Mb#A{`b$YIA-vuRC
zhZ;}ihKmPWWY%_{Y%scZ>vo>*)Q%1hmhFX)-HM){^F2S$*0T1OiT5-eMwzuOtHahR
zB_Hdl`Fb_nrT5hf$?Kcb{ZCEPeSJ>Jy<(!7LR7YQMD}glZ#Sc|e@ACej?CV>cI#2+
zXen107ZsV9eKj+u=|%_bC`kPF=BD$ml9v%%vqbCa>gL$j&vWaQ3S1XsnRjQ0<H06Y
zPA;xR2?v{aWGpV+b?yJ4lnX9fzDMq<Fq~<Wn)K_-%ZjS1OS{YSK?XH4vA(#yT|e{k
zvebWnem*%hwY&KFxs6Fjx$J(uP~P|NSGJ9<?cLZv9QvDjlEefBCtmmX^yK8u!pCf>
zr=}#nxv_DH_w>G#lhutA4lt~X-R(BVqHu<FxnAk3D~0on)6dD6PqQcr<X+U#(b2Md
zbNcyDPp8LQ*^BGN%+QVA_TjR>{l?_ue4yBC=ab!0``b)6YRiL5OTDM)$Im<1%-$|}
zQ)k^mF;Gd#X|AZK=rhZtbDnLrkbpo#e?R}xE>UI$=Qf^@)nUD&+F=6PVLJCJpUaBr
z$N7AFbCZ*gZ<*(0HJykJ4a)9)2G!qkzK5wy3Iqkk#Dc;P4;-VmWE}kU^)=&)u(h*_
zD{t9We+$@KRoW+O?bah<C}B~cAYcE-FmRE}!Ta~;KRY{nbLwfaq$3@U>*My`iR|_8
z?4RN4;__rF7Z+F3j}H%@oS8ZK>+9>2RXl~J=|mp7wA9;kvKnu<n66U&zdv&<3Y%;y
zKW!*~A7@qi$|WE`;P0=m%H{9w<lNlUy8V9L?5C%ve?DjZe#dUFiG56<u-X&k+{Uvo
z{ro(eiVqEHzOy#u-?z)UvcfSSKtR^IY)RhTT`a7uACHR1Z^*o?1}drcR)0S++gzWS
zoliwWqvP7z=*{`}?YyRHJ$?7e#jEOwj*^nn`YXrf>(89G|L^nv-`}~yQStj~W^PJ7
zy$KWnX=hH%GR;1uGgaKSs-#Oy_f#`G{|vj@U7P1G+R+{bs<EEf-`riU&ns>A<X|&9
z!;x<B{$6SGplvykYqp-5Y20p-d1*uPaXzcEHzyW4w|h)fVhvv(_w-1o@VVDL5{6B&
zyUR3X_OV@Q(oj-T`ak)>gM?!}lFsw(YAZfGV2s?9ve19NT|@HBJ;}%W0#}FWmc6;*
z*u=sav9oBY@t@c|6%+aGel#>RF#P-Rm_KNlkLQ|*jdvntD|;X=i?=R&BcbX&ZAtd^
zb)fKic6RpVeYMpLa&|Qy@9*tBIYm?0y8PXe+TY(|b`%_x<U2pd(s^Uj(Hoo7`S0Dk
zS5cwSw<izcGX1krr?+HYesO1~arwJDik_2HeCFHD^_r^1_~G~a{huGU%NPCoQ)!fb
z?#!pBr%%o>T>SYD$i=HdH2vmSIG&qh8M(L0v{P8UEWGrIrN~WC0$KX^(`o%5KYtcJ
zKPQ`Yecjo4w$&>_RthmPGG199KVL}ItKjRa&}04b=l}iveYoW?Gdo|vmW)8ZITi=2
z-|szsYiqVBUwm+uqllp3#Cj1Sp+_}5F}q4KrOa|XTDe45gs+!7CUI?TwENjvraQ~t
zN~wBHNw~Qw^+>0%`We%8aeKWsCLP^T^V8^PmuT|QF3}@>vbSZgPMp2~oTVz)&9DDA
zGw0@}jrsTO4jBGOxwogX=EuW!P#7NTl@8vN;`#3G?&D3Y+!4FW-UdGuP;hKwDS2_B
z(J1xQhN7olt%tAA0jE}*r+0Q1Us)SH{nOLa#SacJ9_x`5PW#`+EA27UNVWL+xuq#5
zC&lckn7FI-b(`h0FE20u`S<(%x%0KZzC1iGU$5fUBe8Bp(;s%O4p0qM^8H;b8?Thh
zGM||dTQUR<la3@D=@9gpV<FgZ;DEy{vs|x<O0GMKpYz?X`OKSrZ4KvTB?%#+LsK+^
zPfXL5-kN>=*{P}8=V}8NuPVPi`+(~IsXd0v|FvJ%3Ah-w`+4YecNZ6zdEcJTuMgT^
zS9@-bC9}iUtg8~nX*?@~mzTX;AGdecX8VVimU?%}+t<y=ma!~qnWh`PDetZogWZn@
z%uY^DF;V<K-^YkP4rZ`toF4zFf`5UcfbJbH6(uF5#3f!+yR^gCP0$Y4bL*9I?G#dF
zIKHQiSK2J)<D;W4{c>}Eetv#=fn&45^R|{2gS<NyRo~y8-CO<LV~&NPu$s?;u(eSV
z$#(*zf9l8FJZF6&YWL4nK9gfEtWvGTl`bwWN?R{4_ctzl<YM>xjj@f5jmKOoQ)%-&
zmGpCS&IO)c9kzCn-(0IZyGpx_Qcr#OdOd!!f@9P2bywF!8b3QTliA^L8?U5s+8N{i
z$>z<ge<{>4&8fTDymqVO!sl~%w+nT4bhrdLIy&~q+0I%Mx!I`ZN5Qjmb6?+?t{tvt
zV`H<U{=eN+t<Xg=JBwOcS~gUCObS~Y_4Uq$>uVy7C#m_K+LU^FhI#(HLl33Tsp@uI
zyQR3!0u<E$ZazCZ`{U=&$<NQtWw^2;usLX{SF2^($45s^GB364E`J}i)Jrt&{5;=1
z6(0}D6hAz~8n!0Fv0F@+K_hI9#PL2^<*F|)^j5aF=uhuWGCP)Her(B)_P?6FhBsqF
z<pc#Mo_O-??CdSs*ZY`Qxi%C(_X}F;b(G)!&w)v*-k`+CaD9FJ@`8tluB;4h&$zf~
zW6jT^X}ZzhLM;A#yPf~><Hv=mr>Fh-bXx!1TdsD&c5kLdhRdWj@OpgVZ&zj5CY7wU
zW04Q2WpMJ%&CN3{3YFyK<VxP(n_K()o8@(XHVOMW8_!89iI0!<n&jP?5xm?laA(m|
zgX=r${#x-$nJg%Kdn;yd)znK%y(g=Bi-m`~#vBN34w=OiADvxp^PVC5tzYt)a0hQ0
z<$0ij!X6aD_wL<e*phj<Ep~U=#g)P8o6`UPXD?ItpSQ%RmFvyz?fp+qPTrhxQHe`j
z@6VO3Ve8|3H>aJosQ#t{$|X-uOx%)pH%fVmvPQUVQXT(~hpK(ij3OM*H+6J$xMc9l
z+ohbCpt!_)dSAD=KGTDjFHfGCX*~Dgg{af&^J|=TmAtelc)$?5tK{L0jmhVE+Bc`4
zUzT%o(;VCCYd<@7`OY@uY8RZbI!fR%vv~ZY9o`F%cXf1}`2XPh{rcsludltizhA!L
z{r>-Q_V)V?j=#CTe}2)^Q-!atXd0!RIZ=MUw!OdqeDV6orCy>z%l&%i&zHZutMqb~
z#>9rzQ3C6kJ>wtkSaP8EgKk$xhsTlPudhO9nP$6n3aQ@MUvDpIoc82;d_Au)-`cHP
za&J$|xw)zE+nY!>ez`OM{{EiN)$Zu%=rhBhk)2<TLBPJQhEK-AVSU_Q5k6(k@W0{r
zc&@PRt36m)>eMePC^&KYm&g6~mzH=6uZ!8~v?b%>i@UqcZL7W%e6we{u_={%s#fTu
zOG~{$wfWI*@yUkC$Cj@v{q)51)%Eq~kM&AhR)5oRZs%j=DAw)#5p%WahibrWiOU}=
z*(H*#7)={H^|U5;ba=E(P;g`d6}~4Yseb(aeR=%;x^oWVGBQ5RY`hik_kQ0|_?RvE
zSWn{bZ*R}@w0Dc^tEHctlXGXs#3v^w2QTxH3=emW6J9q<qk5@f8jry#iOV9#8>UWP
zxI376Dx-Wt&(?40;Zr<aT%P!fN=TeIHC20M&{8f2RxXi*!)?4G8UYGt4onPP9TvDf
zFE;Pq9?yeKtUsU4&i9yW^)>kV?&9ZsbFE5~Zf(hISbY9rkEF5341<N)S}N;hr~6LH
zUHg`&hQ)=kiaFRRes3Ym)&4bGV-G#PcWlPCV>9j^n-RTm{r(<9QTf~XXdRhDi$M*H
z|3{yloqck)xqj&Cu+*=wuAY0#$}LvXcDVZcyTENZk*B6=yC)<lY)(6ybhwT8JWu;{
zy;!CEdwYCloAK)F>#O<AIr05|{rUollRLL8U%9o+&fvbpWF8O6T|c*M{Z?af{};~%
zW{VzW5p%YE=Et6Fy!MW57jtZb+5yK0!4Jff&sbl(b^BqAnV{fA@m+<F*=m1%NxZbg
zGa}n_wwdn5MXuGacC6p`t4mBb%Hn#^JiA&e&q*qk?@EPLy<C2MdHM0fhlGPotclM&
zrkpq(k$pFoZM(s-80IR*;Kt$uE$*-FrIQyuo>8=X_13ae-Rq{B=f~|XeSPhN&N@(`
zc`ZszFUDiJ-`pJq519fNyDcq#e(uZb>-`djNf!OPtG;T9YK0_xes&hrXo}pN7PvL*
zYI<Dp;Wl27z?#U-Onst3d-m)(6r*osd@N`4wP&T}mUm80%{JVh+p{vSr}CeR_c~2?
zkz^lR1+fi+9^J`Rvgy`)w|<lRz#ecyw^c+?kkb@YDg621tQWJxVWCs&jeWJ&huiu2
zg`XYkll9(|a`MKWN@MxjFM`MW<?GWd_7py5TN|}?k;g=(Eg2V^CQJ}ucyzRTd4Wda
zjx610uMXrbVBBy^;&Z3=yWcL7o1~MO8rTn{EZDX6+Wi0X=YKrexkvvKf3dk#a!s#c
zbgyBwe*8Wey;Lzl!OSZw1aEFiO}@1y^Uu%c^AqJ5=iAlJ%DK7e;>KimAr+4Ux3*@_
zG*0Ks-}^OeMs`H@&u_Q$Pfk+h4qoPS@bcx!K})>^MMR#wdUZnjP^fN4Jj2|zTemiR
za<EVCSsm_U$+4$9`R0Dl19Sdgnw7bBZr0kf*=x_v)(Bh_9^xV>*!iaH-5t-jx3=yq
ze=o<7zxS)y?(+A`T)V|8D*pVpzYS_|=ib@^DnqZWjV^w2LeOi9#>XpP{pQ(tUR@Qs
zGIsa2&7wue`($~fOkT{am2f>!+8pwTN!;mv{1HK)WajHJCG~5z@;53kl;19y{E2aH
zgW|&5lNgrES*vpXk53a6oLKPn!$W6K#oWdt=@bwkU|IZZN&fwP42ew-4>q%FhOd(W
z^&mr6h17gLYc8Vo<&9k6_jh-d>;L_U*;UeMQTWIp=?F*Ewj9Z4RZ8a{t+H1PZuod?
z#<l}Y33Io8d%}H=GhePjyg<fXK<95pyWufujjqN0i#j?ycHG*UE&lWK^YjA^j4Z6I
z7x&dxTa~^Nc_xvq`}X#B_3UeF=GatDS{J*!=);49&)b|EUR+#!d2e<3t1BxHH!`z-
z`20E8V#@SKt2#GTOD2~wwl=;m{K8lwnS93Jf%^lKzG<g_IdCpmwDnq>-=n*?C5jm=
zSiZ0_HLQ;Xja_ucFflX#e6#tygiVEkqod;z-`Ue<8mHTRT~YD*+1ZelL9N~5`jZTk
z-9Y`wPtR(rCcHg6+x+F_<@}aKPYztaF5Yl%u66rNo9HL0J96Y)vfnbvtZjccL1c#f
zug_|IjK)vbZe8v4y`N#x<A`VaavL}`7F!Cqx~SOv{qeZ}%lrHC4cl^W8|B<ENIKf}
z^p3oUh|A?=zCT~BUT<*h!^`FKFD>_<pZ-lYW@nM=+1cjL|NQ){;xl7{Td$O$ZB<FD
z^UwJ%w<-nVRnInSC@KBt&ejF>tfIDNP0@)oI<qgn?q{o2>8par$9QXMYJR+0z5c{R
z<>SHEcT|7ZySuwQ{n3%mhUA-jHG`Mge2%EJoW*{z^zxmWlpfCy_H%^w_ncBuQhJmH
zs<~@^etL6zJ3m7*PxA3TUSYoX_x5@p?Gn|D-loII$jBpWHAU9C%*Oab?)`mpZ7M$*
z6h3mXv9*=7tFd@}eSNsal>cW!`)8yxvNo*NFWz9N=JmqzXT$4<5qh8`m%cXj+#JiA
zKOc|Z*;gx_%wzxe%jMwXx)B>1($3B*e1FgO)6>((@9Zo-ImhyI@%o#KT)9D&<Fz%B
zGwti`b{0M5dcVp;<>9J^j0Uk|Gd}58?rDqY?C@wQczH?n&!0a({{4Q>xFT+^Rq(Qw
z*rQ)xU0ofrE@tNb|No>x#kJ22hu7EE{x3YT(XcIYbDClHHJvs-S+6fIE*{=2C3JGf
z68=SoYESBIOxJHM`>n0(V-X<lA}ClnO}_rm#I&=sET1dRF-T<Ml{VusE=akvqp;-d
zt*HkNIDlH&%lzm2y}Gh;{;q3>+xaC;vqXM>f1iGCj^(;FZc7&1b*!8H%uTmRNl9pB
z%)XkLph)`uF4k*`#=$w3#U4|&rWXG{Iaz&Wz(S_q-`}g3y}5B_j^$)FUMZX6D<5~2
zzRtM3jQ98V_v{OjkM~v0(Rt`~&=)kkqT)5BV`K91Pd}f}|9m@tzvbbDTc^#lt)668
zTlMYjZTGcNTP<pSd<fR<W)RnpTN1d~ZHi9hq`$wvPc}?;3tJzzwm{?Ixzx`JT^$`R
zDUObgpm-J047#u-GngT2ONQXH8M4-80V{)4&&{=V-;{E4hFz_ds8&eDv#-ni=6Zd9
zcXwyyXSI9x?tOWAxjpahu7^*bioPoq1dmqy?B|oS@tCaU`{Vcf{hhMbW>;5*D#z{l
zKYyjJfx(2dv$H;44UcE+b3f9|#>-`_`EB#M%FoYyPfgc9-@+;EF<mcq$93l?f~+or
zf)oD-78DpPcJE)7e}5mSm-po4WaIR6a||AScs{>=S=`>LJ9{c8=iJ<6n0`*?)02~*
zuWap;w?DVD__;>R4g-eyHJ><LyTz7fU0oH$>l))%3mO;r$Ii|tv%q7bQpxLUz3lvQ
z7nXPme|vZLG`oBagM)(uL+Yt13!PfI?(8n_KX3nk&RpwqJ5fGnCMM9BQ{<+Ug^tZ^
zcXpTaH|!{U9JDsdwAHzgdn+gu|9D@Enq`)I>VExy-%n3Y&NR!7nz65CbK2QOzO&8l
z?5_uPoei?D>44gD&!iL;6+we)U*6oDT>ShTW55y*!FTueuD0Yo`u2{TK(C=EL-yLe
zPu_uA!heLbb&pB(82Ze!nfdYY@y|b>&lfh&i#cHUY6mE*-q_e|QTWK9_E*WP>+AoA
zXg~P%_4UguD~0=?T{t`2TtE4E-_n4EPV3f9I$^vsOZS=9u^HRv-3YrcQ7mM^KF9o+
zjQO!EWv-yk+P~iF@9!=y_m{W(^TGM~wmFu?Mzz06(k*0cDhxD(mLweSlfAM!ynnHK
zzg=<M;ROp6K;v!O^6q+pQio}_*le@hAd4y2yH<w<ak6nmWt+Fxt>e;R-Nw}I%)gp_
z5BHR<*ZxVpYprVCGO4S><H*<3`ui{JD$VYbv*kLlE_U~ciOTK8|L^UsUKz4d2vqY#
zZb~^b&vv#^>Zz8U$Ii|+XRLUw>mGgb(o*l6JB!ut?k-nPJ2OLXRp%a^wM!oDa@6ww
zw&5-38fLDDZ1<>acc!0>y$<E#Q$g*{oqBP5Jet{fHG-G%xcAE)J$}5MSK7?taH59?
z$JyED>eb)g%&{u%(u>_y@Z`iqgNFth8XcaK)g~%Bv;F%1{{5Y}y;Y@bJQ52EG}h~u
ze|0>=dS&}HtLqutuCYlct1-pzS_$fNF@SO{$cV?sco(_%`|Yp$3mTT0H*X%pnz+5Q
zc9p&^dVGwx=-HW*dn!Mlm}gr(BYREMR*<EjVK#<8fB$mw@u_)D(a>Ape`5dBb(wih
zjos&O)G=OIAFQPGD6IP1n}^4GrKjjdPqVH5c44tQzm2UeCkMxZl9!i2g6#Zq2AP*s
zw&mZKQ<HgjfB*b0QSG30F_z!n-Sxh;B@@(i`T2Buya?aZ-iN2ej{G||Lwn(Sd8W)p
zUVoX%LOY%vGV1K;`0(MOGiYRJna@lH36l(gR&McSZoN|H9Mz1|&n@wtZ3fDOTA`~9
zN?(OcQt=deX64>5r|RA(!@|lcXkC_L_w$Lc-M=5ny8Y_;GU2R8f*qzEx4X~qiShG%
zi*-yJuDhlR2?{1&Ugo<c=jJ47^E{B(laBZKo|$2|GIFz;K-sHVEu6xea&8)dMwxB?
zU%X<@@bdEV<ylu(fvUnrX7))Ig-P45cs|)*bW%AjJ69?@_m+B`?7YGkMhp72L6z9V
zXa@%dP|cWdu<6I|-@>|4CVunnWS`AA*v$U<+wJ_#$;bJkwq_me6jndf;#^-}4=RWs
z8%OFyZRt2USzTM^l=R!JQCa7sb~Ah^FSm4YQAzs!uwDMdEK}{|<9$aDA8vkkclYM%
z?|GZj&T?5yyXZY#PZm_rtPEP3bb6Yu&n%OZmal%Cn``~$#l_}3JBvRbm#+s^?pEvH
zKKvN<_p;U}`^wyZBGSoTQpsr#B~~aYDGBLWmAr76sugNc`^%*A)04!{&(8A5+s%<M
zOcD?ha++<H3$lH=-`qv<`|D(EDkkiy{9o_2=jrL`lTEY3+<GJ&3knQ!Z*6(_dj0;g
z<X)FJ+wFQ=(|b;G)UcQwGIV&#_+hT~X%!_U!Q1TeH3xbmjX_ZhYHMYM9Jg<MeSQ7%
zs;{qNb{4hX*;%}~;GvUI%83Uy?;joQ2K8Q5y{9ou&<@vI<}=f2U(L@uLHrXxJo(fh
zvFk~k?E=FOa($j#bPZ%A^49Ks`O2uX!{f--{QGwQ{(ScT^XCu4ikO{7$NS~g0~fh4
zfRawt*H<6^e!u_t?c1|cG=np4ZhC4w{lxF5)8jitwZ)p*`IqJ1-gd4)Q|RK3C;r<E
zk1=cXD~p7&?P=J2fFr@O?`2+Hn1zVo#PFvN54ZD3oAn%QW*3%`@!41N^Tb?h^9ExX
zy|_Iymif*udVWs!)HL1Gb1aJw`NS4KI>HGWgq@-hSn>Pq_Q(whj5jx>Du=EL$y+(O
zWA%<z{EsYU95>i}&<zmtNL~Au^-Xc0q>GEnBv2=$^3#)r@%!so*w`-ZsVu&<%=h)3
z`oF(QyF@gd=3140`0^#?=BCs?e?Ff-S0EuGa^%EBWl#s<@$vrhdExtNDqEd9MPft+
zgdbN-*(3QwE+T#HR=r)@J33s1%HQ2tSp58)h)%==+v;x<6rI_wt&2T<Wo2-uxV{{N
zu$qrT=A|V$H#Ri7^-6)pC02#3e6;7-xw+Pr@ArOxac{5n_xJaoA82HLx2mt>_6{ut
zF~((#I~!dOuqL0GeryJl-md)}9WG3<yUS9)ytv44Ze_OD6b;4$FE1|#jUv_l{-)wF
zp@CQ0%%J>TOw#c_-UEMse}8#-Ie+EnXQ{uxy#);nPSp-Ks{E8v9_oH1H^wyZn_Ko;
z`3WK!4U1Jbs`l+`w=U!f>U`d{2{Zz7b$fn%m#Fr!Wxlf|>}o77E_Od2<$Zr&Eoe^R
z_qVr;y{GHR*w@W*YUQ##Z(aP1XJ^sVN1M;vDXXZoe0zI)bNTzYQ&Y8{AG_7DqrdR_
z<hW#?AeqazYgkg4of=9H_+-D8HeDEUZ%^gh$D3btcXVX0iP<?x(l~8H!9%A-E}fsQ
z`~%H}EOcV^n{Vg4%xC70ce~%8n5KI>_&6&&d+-tu!FBQb&-F+eYsBo>pxvdvUNiZO
zS$6I{ae2wYow<v0vfr|1zxBw@EmP3{pgKW##`ASXogF7;Km7IeHE7I_OH^w^%1NPT
z=jY3()q8kwtc~8T=H4d*n#D>y+_o|Q{yu|;AFjvOKRq+kn8Bq}s5O88-)VEL%e7>>
zKFnIPBx&v4ytRKNJlL+(T)l37?2FN{9K&OGjE=>a9s85e^Re=gL~_r|{2t5no}KAE
zk?B2^zpncT3Kp)GulvymZQ1PEvq!=(iG^3nq~P8E-|zQ(&N9)QtnRNCzAi?_s-#2C
zw#p(}BLDtA-_>DjHDY&}tc~9O?Be7t+wa#&=jZ2xO3TvEQ7kSD$v(ab+#6J6^x~P^
zxl<J%)GBcNxGEw1e`7LeVB(vD-2<Ty-G*vq7E|0^T>dRx9lrk3ia_O6AuA94{QO)q
zWQD-)vbT$(x97R^N=;p`K;hNZ)$RwISa+7b7K^X@sT#dKkJotb*;%I9C+68&L%YaL
ztlWz{O<m?a+tiTJ^D?F9W_k~&!Lc_+$3pn-#+{Ma%zgaDGZu%{;kEhMxslm#b+X^e
zWv{jG+t=UOaA;wwfrER(jICvzTg&+SJohFue~mA#ba8odKYB-j;^*h*)vLa|0FD0I
z|NGHgey@_fA#`=vM1w>p(CEg-q@xy9Uo_T6Z!cTdCL|!>u)FN7M(8RLb^m!Ppoy&8
z+j{Tqt$zIQ<E<@IJNB5xY@Os4F0*=>3`-b8W`pXn8F3cVZtU(;Uf{d&nq7(G#H7a@
z$FDovFF$H{sOsxt<~3Wt_ZZec$S9rjs0ox<gqxb1=k~Y#{qx77<b}Z7TU(3QO|1F;
zE_PM;`g5<YuI?06X8Zi?>|)1ew!~>Cch~>7+gbel*@uURE2^qKfjWMW+Wc;i`b6`_
z)pfa9U5i@{R&U)lZ|gRO3;gaGcNkN)UF$0kv|y3wPOh4}-^pH`_iyFL&c}Lo+w^&o
zd)6DKtu{=1{q*$o!%HV9DNUTZbz1(upKU%fjanr=jMC0Xq+86fueakAQc3vw>ME#H
zB&_atVMF5K<<|nfefzd7?d+^Qzu)Z!&F<~1{jDX_wWnb<Q&LsJH{sRGeGYHjbxI-o
zZK(KE=HxRPGQ#E8w_IcUzyA8hYkpg=J==QiS?!i(%R4&$ls`W|A5@Wmriqe|_ZjBi
z+G6l~eaud$e}8_023ic0jy!0WuS@v(>8Z`X@{i^BYmc9rsy%o9@9*!=FY}#!Vy5x(
z)2|disTGGh^svQLC*59re0O?m2wRWg;zs_)?`yXnm%Cl6_$7Ukx}WE#8+L+%obK16
zo}HWPeQix-<hC5i!be9Ee|>rR<^BEpm1kB2F7}wFBN@FdXQAh0HPC#f*HkUb@AIw8
z`=m^>0yd}lMs3MB2x?V+`1DCBYD)&Ffh>JmwN}8nQX*N7e@TCuiK&0i@hIs}tMc=D
zB)PA!dPQZ=kBRtf(Wek0b=)JX$MBv*;DXz$^}0JcT1wyDvHbJrPsGL~RtC`E*{La-
zOFSkrH3+Nw1#HWS<dd~>`Sj#u<c<PG*DjI5dGbO+hrYbL{PNmb>Fs%UmBQ9Ue0g{G
z_RgMp9e=FX$?V$7$WimARKMi+GWBiyKh<>{o6%|U%kdifI_cy|ZTpfKTkL+xF@G^Q
zc7<cEa-o?3!=Av^yE;1dbe@}UfBw`|?M`WPz2awQ5|4I?f*g^p8@szqwf5H+8M~Sp
z>tc5oeR*-w;9<f2eYGO`aWeMy_IY=Axi0mZdSiROeD3XSS2<EA{-3xyN`QUG3%m1|
z)xJE>_s@}M_`+Xn&(y)^b7Z657qtaFTfc47=F(6Sx*4;(toO%9*~=H7ot)ghv-o+?
zA{Wj{>i&EOw&mUi)&9@V&sQ&gb_UeCdUJDg`%X6x4-V6;D-RwVY*z7`Go$eFvB2ei
zbC<KN>Nv6g@wTlql!ewF)lvHNY0cI%vrN1F=2{t+zl&KFwpJ>wJ|IBAca{mLk)xJ%
zWkt=mo9QL5u4tZ#I<2Uvc<$Ok(A?ls@9Ci7l}AUrpI_SLab(HX8S#s!hPb$>JhHC+
zW%Bg&^yQVGpMfTvKvU&s{(-c#a0-9A8Xh0Gz=07IP1o1O&VBgd)KqPgygM`e=2}^v
z2aSXUFY^hs5UFg95fbP%yf2w-WV3PYW<#}$!nZd`CztI8O}ad?x)zmoW=3OwKR;;I
zMBJW=GxP1|YX-}_7Pzywdb&~SsepYol}yacGcAkNq|I^^d}o=Q`<#Aumg}=KGiO?t
z>s?;%|NPQY@4Ph=J&rtEC3YmMq4;icvX4lD%7jmykCPZ~#)PN|3KlkYi|YsNtEmK~
z-?h=(S43^q`u*)~@wy75^m8(xF-!gZf1E&V$-+l0zO&6rzbW6@o-eQJHDy8M<}}ds
zPS8@Xjm6LXgw*}sJe$?&60>g3vuB>K9Ceh0IFED)HZ!yHfhHwDGraBcbsdIn-qZC0
z7do+S%f0Oen&RG+!YQg1QZP?kNa&FI{F*~2CMt9C@Te3%I?^L$+V$(}>uiB84~t;%
z*vV5qc{?8m2L^_X$;Z=9PEvhyclUM4cWr#Kr*7x(Kf5jWcE-g;tm=MqIE48aZES2n
z1$JX&BWUvZ!GXr**CNFxepuL`p`^6l2ei_qjaRw)+nX;hFCX7&>p#y%#dp>e=XSn{
zcC}Si-`|1y+&_Q6-#^nVS88L@QKkoPZf<^ad;9sydoM07zPTZhd9u3yH3bnt@uF0R
znrk&yvy$CiRGMyVOlE(-_q*KQs;`d@HnW3L#sq=mhawae6+v@!XJ#4~e|qBi>)YGY
z6P4XTQy~n;ub3qt>rpIzcE%**Lc@Xu3JjA}z1<c%wcZI+pLpWb)LMahgJU*3y-xWb
zn{nW{UG_pB!HMQ){{H?hX^_C6?B1thVKL*%%3$NvQzECP>2e=9)+@a-Y^~IE{dmx%
z3MiZS%rvq*|Nr&+{TX+6neHxo%LGa(clOmr>$WN7Jz3RYaoqL0pr9al(twpq1eB8`
zl3!d{$lT1%uNJl@!lvfOgkN7@Z%#QWv}ezrAD>R^XI@=3bu<6Nt=ZRK+}x~wd6_S>
z%)x>`KR$-63OPA9R7fyJL}2Y!-`GNVLBYbcI|?69GEVpFF$DEzn^?J-Hq`tqk}}Bv
zjg2hznws?Q&rcQKSyw9K3JMHBLFe8l0~!m@ySocCc93+W<Kenh9Ve3SSf1_e=-AT;
z>IvMddaYaZ<OJh^DVo6_zI{7(`}XYH+w+5$dWnKo0?e@}{PFE}{>sSB&w{Ua2&r<-
zwJuM)xX88U=hNvaCnp`fdR1Uz{er+VN=l!Cn3$L{?(efb+9kTU@bNJgHmh9|j`vE7
z=jZ2xhTs;tbVh8?lP!IH?dg@3!RHEoe0b;#nssa86rOYa^XJdWXJ?smwJWVx1P$1n
z;9eWIS1LX}9yFqUWkn!pQtFzO+M*p5AC;b-np(UrZ&SiSCa-BaFKx>I{`&gk*RM_W
z|Lc;D^?3gM^|e#p-cC*~Yq3CO?92atf`XOnvaheZxFs|A)HL1GCnu|)oT@Dzy)7s4
zKm(&ow^;9ayWc$`nn4?~uj}2~l6kmOSbfgJ5660?fBgIzxFsX-+xz?e$9g0q_tjW_
ze0;qAV65wryWsU5=S5^>d=56TT9&*J0FAZ1zdv6uc9(!@RtP9#{rvRQvg(UQ=A|VM
zZ|Cn<^`B>Ru6Y0Nch<MIW*`6g`T55$UzUWgkNfiKDtCT<zWFr|m$+%5^(Rke%JHgb
zhOLnRwdUOW<<8~KG){kZq*J)&*URNHRwWuf7CfAsiAOpFK}}qzRxVI!0U7@I#Aew4
zYr>8{$Br#w31ag9+ncyH_VF*FzwbXiJDa^d@9wL2NgW-3)IsI3CFkEif9Bjiuu<>R
zwYAYJV|E&e>%|;6e!LyjIxfFgdAu^l!GWRh(GkXk^Yd&$>nlKW!L__kdy{)^rma2u
z`d!hE$*1KFrIRl)#xwqu&*l1?>3lz0p=<SNRVAUBZ}09Fx3{-9%emnI8aoJEE9E&!
z<>2AN%?A!RNSS0DXkz6C_07J#xCk0-KKEdI{r|c@fBsawUb{VFUybE4iJF?44<A1k
z7O!=AQvbBiuH$lJh=cnBh7V6ZRVyhyQu_Y>e)@?CilAQjfddYp^#VHzAA?psFldLZ
zSy1=)*PVT}vlqE`f4UkTKhY@FivhIoWP+kIsGZ)}*r?(;sipY&xsTs&=kMO^rSyqw
zhqA&m=2;Qh(|0_)r=%p5_4d|QP%-!8L*jvt5tTPKBp&INo*unDFK}yCC@53Bxv>#6
z3A!SB`?}(pdq89R*VaZCe|Zu3e5Z3e-_>6yl>}R?nI1`A$<k1I6m~7@*O!-vA0BQ$
zIm=W#c$v>akBLeQpz*>!S?f;^4l=*Gw)XUslao88&DVY2v&?rksH$MtlzRHum6gGo
zAuA^AUZte;@PI7iN8<(UhH1T9zx7GK(^nGGjLhcb;7}-gb3-C|hEXcl>1n!)tG~Yk
zrRdAc{eu@ev4Ym)U0UJ^>QW!=7UwVKP3G~OtakMF?b&sIe_h;JoW97l``I$7PLGPh
z(;8{XX9Nz2|B&AC??_(u+T1svoc^*)9o11<znhtz@4~83?WE&<Z)H3_b!?rsul6^n
z9$&My<kgkV^LD>^^!4=_zP-JD{AjoMT;YWKd#k5&i|b9$i?!N*#Z_fbrpV@MH$!vA
zD#n8jdGqS-yqQ;7?&6}f;P?0UpKqqmH_W+Vz{Je#IYmP-{ro&t*Deu;H#asmpPH&&
z{OwI7Xm!*iRd3J`i$Ixx*CLlr&?J80qa%uLJrZ2)N`iB&)s{;!=&%Y!X5Zg-ZCUN3
zrTa@=T;du{v#(81cIN}F{Ygkr04)SO(8#>RV`5V~zx<^gg~>(F&z*gGdb*@-l}YvY
zccpEI<Lm!UeRFfOamon+CT8Z9K})-4XDchMd(vmu!P%g{cB}cFf_%Y=9Ot9FPfydG
zX_(A*Ah<aoY;Bb5-m0%V3Lmo_?G}GNJAa>I<faq`&`iR)xz@(@|7sX&YHDT}rFO9`
z>~Jx9m*E)0pHg_Z*40IcHJOJ~Naeuk)7=M~*+Ju8piv742L@0(C}w973oGl!ijPUx
z)<!?)X}=bw6|&;MpP!#Q<!r4ExA7+b{`QuuU1{Qk)9foHlJ7B!F)BCAudJHDb?uhG
zbg@_VmlQ+=3%4I`=P&;Er?N{_o6Ul4(j43BX+4t0Oa}G;YHDh0K@*!>v#ur`>ybQ{
z`cZCs)>SRL|9^^?`OI8Y*Cv#iIiahX>Cxg3-}gvIXRo!EtNnjoM8Ne(EE6-cq-Bu`
zXl0;d6U(_qm41+4y{2kyto)oNq~>$t(a~<rz(p(!ptV1;*3-1Z*BO+&2zYLJ?Uupf
zIj&05tS;ZbILdUmX?1t}Nsr#1x3u!}vp4tlN;ep#o)VCe;jwVKw<2)yi>s@}9UUFF
zWL@pLxY+&i+qb-+QD{(uBzk+^Q3>bGX=mrG?d57$di1FwW9{4bud2H`TvW>6-&<P$
z|6k3=qv9<sEuaJlnhseRw3NZ3?2W|J)6<_{Sm?aOZ|<z6-qR;(hwCxO*;E{OaIjg@
zEJp$~-{>=QlILVK%fkyhICi8@=$dS(HqB6NsiE35?suVQBsSlXD843<?EE~rre?q4
zu`lMw<`^HlV|;9l@v#){_eGr@f1*QIg(N;XF_EF<-JO|tcb5mR3ej|OasrJ9NgAg;
zd2w-Z&8L&<GL}V0B>Bv9Z@Ij?vs2OrJk<#r%iCN1J!oT+Yk8QvOPq7w%h%~Wn%q2$
z_gHxndJS(k+FiR9f9U15jSaR3?p(M&88jw<x>zAZO>koU+U)D=Kx^Y9jngLB*Vi#@
z%fIjU^wd<#@^><Jc6K_kyQc6+8nsyFUs(|dn!NUzX=IpujOY0lb0sDJqT>ovHy!9d
zAoxJ~0n-GY72GKdS}n;QH6d*b3+8SuD|?a7Jcn5=BD=jMubVH1X;)PC@;{fqns~Un
z9J>7J>1oi~6#w~l=W-7oJO~<SIdH%sYFp09IhMsAzI=J|`Mmw|+TY*K)wO(of4}^j
z^`454O2X=XpgPc{U+%7MsbHmW1~|>dTP)kQ^2?7Bvtttc39JIwqomhO&a+ikSjODP
z^mFajXgTY5?hj6NUuz3<5SjR);Nv4#(3rNEZd8k5+uPgQKfl}k-s=1M$jxn0TeAuw
zgMo(@xpvPfmb$m6ax!Rmxb$_{zkk2;K`W6%xF?$TGcBtArSObJ&Cc?A;*00{g*Ihn
zTbC@f*SFjF_fhiiLk7owNF{HRO5PT?uf{TlTU}{l_tt4Ui=R(Y^%fJ;jdD3N!%)&V
zjpyd(bat8M5UucaJ}ZNkR{Z#&7+?EUwEo}E<;F64@%!e0R>76N3i<Q*@69clldV@P
zE9EPcvm9YRaC!4)t0R1O|KHjS>$Ce#(U`cqTUkjEwDM9Te4R|MwE44NUtdqr3~sZn
z{`TSHad}}?FA<RCpjOq!B-YSXArHU2ynJPSynMkHi)(A6r>_cKUG(mbC6|cCgBVv=
zmnZiN+#O=rd5n%NsjXRTs3v+rKH_zoQ7p&bnNK?3ddl5aR}y*$njioz1N-yiBWPiT
z{r^A4nU|In+PSTZ-mVuH7Y7>IOFcbJP(RLwK};vY0o1pLE~<QUW8<3m{qu5enn*Tx
z{HcB?arBPx@v1ZIy@vJ7&kJpi=`!45sM@j6=XfE@!Al3k7?OL;14LdtYx~sUq4Fu?
z0Ox@d6O}u~bfxCkev|Z^WKows*CbO&*1GJ;!^7>M&R(amx{z_&8OU<D>+50{`^+=~
zji8uhUfQusN=c}3N6iGT217N;KF!oFrjL!;2i_mZ$bS2?U{|`1%;F`CW=w}0o*ppG
zeyf@6!@t1%f(U5gK^LR`oM&#$AudnqKYn=V+{P<??9^24pP$d$cS@PQdH*K;{Jf<#
zKR^BW{d=-aWzoLB-=ew1bT}M5JUBpA#y4*sPEOF0(<MPmy`nf>pMZ*8ar-2*V=l(W
zdf0*@vUf*i-`{*Kjpy%8+dYR8SL`{^cOkqa!QsX4#)ZG#85TEo9te{Oem~`+Oox(?
zoLT-opLsTwFRrhbumAtoeE+{+tBdDvop!L9eX><)7N{;??l;$>?9B{g89A$x4k1;q
zfDH+bese4s8`|61r)r0XS%^#wKeV)e#@em?y@quMg2L~af4TR%z~|K8J%7dA<xl)u
zr#Z31#jX6^orjbC?Le#TL8YXkGaIOR{_)XKMvIys2C1i~J^k_V@e=>}^WN|K-6vyN
z#2_Fdas*T+Utcf({r!D)&q*qCtV+H5C%ZnGT6|n#s;onq`JHf=BWtg%iFE$<<|e30
zWw^VmH2M3xyP!dgY+YHak_U6k?>)S-GMGogpn;j4&*JdIcXxMdMr=^9t@`rd-{0S$
zIIjKuEpSbQVR@LRi;_k;PtZ%In&ci!&OaG<FLDOC9O31YvpMnc@$r*0jn$3Q&M5fL
zvuSu-SUgoL6x1uQEPC<))X+{o)^qN2+POKM^XJP$%k-7O%Y)X%SeA!6yC^Z1vm7~k
zAo_rj!Uc(9$qE~lsO;P<^C{jgPo{3VR`lmb;iaYC><#nh%U@jVzTB~yt)i+*NJ7Hn
z@v&Y|i?2gWSL*Sx-o?4Mw^dZM)b5z9A8%)4V*^?QaJY?Ex$Mo2cPsijT$J8%bnF*2
z6X@!2u>hr^*xhBIWub5G?w)>SW$@#N4;?3~`7-?Y`F#G$@bz-$`S+G&TwGN1_v>}g
zh)_fI_jel$9y&>x=k*x2y}7v=+_0<@ocQ6{q=t;1lZ-N(uHD;y&8+aKt*gr&MHQ8n
zlatjiZ%+5myT8x(&5eyWH>dNj47Qf*0S#=ty1IIE)z_?~V?8JT{QSH!e*eDWFCXRW
z|A5xSf@0e+>Bs`VxmFCInlxyg)RIliCm0iYEYo^UzFhRXt7FHRdwZ)jqqb<|-rn}~
z+uPeAYon%uVrgTNE5oz1v)w^M4tuLiL8Z@;PGOUb3k=CTAn$~(4inLjn{&9G|MBC;
z&bPK?#uvB?PW<q2mHV0R1!rtICa`@dx>)P#5|i}o%uG-v4{gS`a0;&oT+DV%LN9)w
zPb-(`ja{YMA0Ho2KRHP?WJSQi%6H%I*ZbezmJ6!S=FOY;<=x%c_x4tUW)VR>C?(@?
z)+4qDo*$bb?`*%_;FwK5$GolA_=*B|y12x3*8Kd$*w?55il@9gI~Z#I{e1rA)z$6=
z3lwIV=bxKlnA{<%J?-<J%FoX}UUZjd3|Qd6xXge4y04R6kGzf9bY-65G6n<T6D`S`
zUe(Wfbw4V5aZB~WI7Qv9>Owz3-3?Z5u?ve_xzqn2yu93hGN{whCvWcu>IGJOPyo%)
zUtHvRZfbyB`S*9RT;h5@o72u_Twc~&?3*T7D6hNq`69cHYq$2V-5S}LVR0xi@xW(G
zecNEciPvwf3SIr-<Hv*7uTRh4|M%KX+tgE2BC=oJ+^jBXl=9%i!^5B@#Eaegtqwmt
zIa&SZi^csrN?(gHfR^yy-j=&1|9;%V0GBxFb#YnxM?h;oJDw$k--~v6vVT|ETPaXy
z{@ff(h9~Fd&IXNaeE<FXH>l@uZjPn%vokY6TYT*7>>83wK3`cG3>pw+WMn+jFF)Tf
z`51$Z>HCEp9lLk_$yd#OE0g^;H2i+J%ai@!>77}o+QrY#9K3yd_QS*Nn=>z~6+JuS
z>EXel8?{9t{oI_G9R-aGo!d91pO;Iwu&MiFvDAC|veMVrSXfv<gTGZ@Up;*$)!AV)
z-&*aIm|*4NGjpw{AM2H74%m_rxXge4yBy1y{dKliSA{-)b#?WVlat$7xy3+pbO$bG
zS!7;PsrvHb;9To+j~NDvpsA0o+1JaWV_jV0+U8|<bnJ0h8@2VsRBdr#HJ^mLyGm0|
zP3b&vz#;3}8qbT1Tsb*75*{7t%(=JcrX(MuS>7Fw-DPh<+XG5pT}iyZua>J_>62Rk
z_o9vy!g}%h&KziDmNZJ?fObc6?(VvJGrZvEr>CINr-lZGE^+;He|~<Rt35X&d+y>{
z-6EP!%l+npR==2LUn_~Ob#eKp3mWaP0S&*UoSxRpDXey3xxf7LGJbykcXxIk{{4Qx
z|J7BY8?PDV-?O>4F7`D~d-eBsOA8(z%DJ`WWF@GjEW+pdr1$Z*hK!z%2FHF#BzH+C
zC+%7-qonj8!*osbcD=1xS08<PdU}dRVAHv|*4E$8&opjl;S}0X_BIMMceFA2c*czl
zipL~C4SSmk1JK;x`u+bNUG}#xP4;zh`DeYI#U<G%`T^4h?g-{R!n^)fC!aa~PGfS1
zOI&{xuQ;Ur16s|#|L-^H%*)H3etLTP$=TV{)qH0O*jAP7`~R<cnr^h)qa&S+pz*w)
zpPxUU>~E)NU?315A1|gG^#l~tf6VpxJC4nmb|B(G^#Z-b@X(*a7dA@>7S3E7y}d2d
z^vn0}%YA2?Ra8`fRy9si_df?3oxHLFv;Y;f`T#t8y|4E7B>i|h(7H{~ii)?}?}K{9
z`Fp=ko1Lwy^y&VN<OyEs54f{)`IFB~KFub>uxsZ+A0;8rr)OtRzqU4dvTn3l)b_mI
z_ICE;{qoCmZf*jN$ful~WLW>N=3f2(zXs1^cbD~o+95(hPVIcM5=JQ;hRLh8_jI`M
zyyNJ|Yz#ervz)!-^{z!79x0$j=b&W+?R>H)rt8ZyfYRXuiN-s7tIPNOc*MPC>z4fc
z^AZoY2`ans9N3swv?b$Whm@(-$w{io*VaUWR*xNO<=(y4PD$zELD^b?^$eNwHoupB
zu;z7_lG4PJn^I4McBg>amv3)vJskD=&i3TveD~^p=eqUDc!GTY`~7}N^Sn96GXFjv
zmj|sszqr`l`uxhUwX+T_>~QgT$I+qP;NNTbQtozIuOVlnpT#e=u8uuMpoUfA;Wp3~
zIM6`$nh3*v|9)ln7_zXlU*1ue{OQTb!%MxV|M+s*|K|35`JbPkvp@Lt_4SWmzk>Eu
z6!ytj9(r`No0IQa&?hCKjXRYUzOl+NR5!#ka`!NQV?KDm!F-CRi_!{F?J$dI-o@_y
z=Wc9Fj@(_Q`||Sg=iBesDJLFk0p)ZV%OcQN92<{>!<LMT8?J52xTv%#<>a9)nU_C&
z`LZN_e;uf$sS!Bujh(KN(8e8m_!n*c)@7)se?d@#TZUKnI7`Bot>31tKCPl8c)3^F
z9F)F6l~TdOL!c!~Jd#EQ|NqryU0ZYV)Ku+{U%sTAnPIpk`}#S_KcKmmW;R~XMA3r?
z!PHk*R)U5&=Ux>QER>wi=927Vl3?*+>((*>gA-G~i%&f9@M2MU;p1bekB)SLR^9g7
z|C^x^xM)MdK_*a@a(8z*Xdtin*%?LOStg0+`bD+HuCI?*_nTwUBVpJyO*h(V`b}ZM
zPM$N5=WdBDa!!87^mNVE`W;`6J=Dm*x2Np&wYBf`+B}w&|NT{JRrbcCPsZ}arKR38
zYC*A-bhPWw-|zR2^hi$r`1ttcC7#01cb&8@d!rDyr{dg(iZ?e5_w3#4InQS2>eY%$
zd1;7IJL@{tdyFaNhihF`n7~VCrs+x-KR=gxVuIqGeYLlD=6`u{5j1$e|Np<)udc2J
zH9lsU<#rjiH8e2nDtY;64*$oC`e$EWS{iI2GSTC->}Q2C#(jpz=CIUlS(kT6WFp7Y
zTU)bB-rSg2`S}@Rz={CHw6n8Z-`v<}S@?+M>8YuYK~qa!Q#wE^*Ec4)KHq8o=L2)@
z?QKuLyu4gdQ8B@?IIXFv>BE;VCE=i99B$Bv36BcHG{)~ehRb^mm&;YpEpeXd@`N{e
zZR*EIM?r;V<>zM)-@ZM&B5?7K$NlyvW*92p-Bqd#no_f=m=Lx$s_@H;z;p|n>Tf!x
z+1HkMP1TxXS33)|n5Ob`np?lz+5!!wrs^nx`oo4>dpl12KNEF&N8#fg1rM2AyF?P7
zotc?(X2wL_=xqT@y+lC^;&+w2oMB&Y2ilA<Ph3F2;oKa{omF48HYOcS`taZ&s3X53
z@b1U!3p-A{J}G<oj!EHH+yA`DKGHI;9o3n?G0CpoTEA<1h!T^4jLez0x3`1l1Lfr8
zeCAqBy<h)dHa<SSskwP&@N&NUb-#537r7)}Sl|e%;xjHTYTfC!CSv1<&!3a;@2l;R
zHt(CF8N6$ynUd0a{ptNlb#o<?Uom}dG&}I?!hGMxHP>zhpSSj3Jj3<S{Qdjw_wCyU
z8pNDm_v@tPFNuVgmzLIix#$k6a5p44f|_!_zPx1IQS;O2?(T9>;df<qIDY|)uWy;I
z()q=KXYw)}O^(%kS(LnHYjJGgwDcD|XSQAQ<NkZ+o<uNT2Je@m;F&H;TS2SbDn2B%
z@k+Y|1PGj*tp5DX&CQ@OH)->{fSpCD*Vf1TuZh?Q+R@7^ZFWS$*~5e5@-pAW;p^jK
z_SelVeSHnIU2~CJ@2Ua~rB6)nUX&RBw!AH$f2yR?foJX!YyJ;Eg?EQBx5mWE2{!Ud
zn}Hf!d@>dcpryl-RwWvsb^4%<Z{E}OK=UlOwq_e=UQ!WK@n8Ur9)SYx@-knBDSEN9
zI)&9&&0pAY!n^4B{b`x=H|jD5aIWAtDakFIBg7Qg(CweIpQGl(lY_@3;tz-`fZDfl
z786e_D}D}Islt}6JJUEl?bnx=B6=}1KpVAI1So=*^-fap{PFAc`pW-*zcYeXcF(W-
zr3spkvapx|sz}n#NMv1KC;M#TmX3};#?$?alh@k*?73c+{np`uv`l|xeD$?kE3e(!
zbM4liP9fEz@<JCS@G3J^FA>lPa-Xa<Xy<tCuP=(Ay$ARA&j%GakB|3TA69HIY<s`|
z|GWtk1a585Ki|mAz9M*e-=T#aC-#5YWt;y~_&e(nZ3lx1pss!c_d5MPk0(VhFRAv)
zSUN3q>XVbL{rzoe+}^4mKYl!T|K7iqOVpzH8IPoKTFS*muAsFC44^D?V`DOCY3m~Q
z{&S!~oll>X^6z=v-5=xf#J=)w*nE8kmK~3eUAxtOz-+-6;oV-$(y{SM6HgqBI{oJ6
zX3(-L(8Thuudf*@s;fae*k+mM`+@RK!9ynJcD~g6`)WZpv9Pln*Zry3ly=tZ{k^?&
zuWQHb*Z}fJ{gb}HGxP0|i%kx&ui)@0$n88N!FtB{SPq*R!`g-?2U5(Zd8!<;Dtf|E
zS68<tax+^2kJ`@a?|L5}9ZkNwtMtm6$jRB)*G<%qw-eKg@wm06Pw&jXfB&}R+`MGT
z0cx_$+U4f*M6~F*^0R;EZ+v60;N)R!+qK1{f8#YhlVd6jzaGudsOAqi7`6K8f0ww0
z2M#!Zs&O_R3DC+c&j|{Qpk>gBhuc8OPAhcPhr|5#58l2#du3&C=EX&<$0Xj}-+z9A
zV>4(1!qL%j&);viedbt9Tot-{*KRK*CBZK1u%ODja@p;PYyS%UxODfaK<9~9S63e|
zzhBGlu-L6v!Zb^SQ&=tK<Rn#4N7=gkUBEIQ$!BL~I$v7i$+)BPvl^&fbzz~iL^7x^
z>fFw!>Nm&Y-I}frkCrEW?q`lRyuEg7`7bG{<Rr$OYq!?#*cGGn$SU{Nmc-}h=K9RD
znR#icck$z6yr8|>zrMci6jlcpKc?B&c#6+O=kImhSMw9pKDoTix0nyKJth47Q}4Hu
z$xV$-*KVCXe`7lPjhwZ$i>-t*SK3aj|Ms}wepBh|uwP$aKffa;A<?mXxjLu;G}o&1
z(Zj>-pmoZglhp*x^I{mn*2Q=}I?{P^vO52R4A(E;z8yP#`n2p-fr%bY-fwp<=YM4H
z>+!5Ld6lGKCj+SPZRe8(ZQ}(k<4}m-SJRLjHM{8Psg3F9=Pg&;vk%mxu_{e^cV{Q4
zUdX+@?cvj>N2BxiKHavfqr*i_kF(?03=@Z4$IYJd*3S8}<AsCqgDIe8sSnjvR9YmB
z(=P0-E(aZY0NSw+X$vzgh~Hla8ZrgttgYF{KRi4PS}w%kJzZ~U>FaATdnyF;^Ydq!
zWG=dRQb|b<v@+nBgT2ga-RK0H#Au=17ta+ItnEIfqVzD~{k^?Ee*WB;c$kgh)Ku-`
zOFSogOjhFsHRqjLx#n1vc9p!mv@zqN((}iUPfypsye`%{Jv(`#lIsk^WVV?`sYReJ
zjFNUZ>yeuWtYj1mc>4r)yzfkwV`yKq^*U(r0Vlh%dtb`oHeS&5*2Tr{!RzB}*MR1+
z40CUp*x1;B29|1neK~mfvT(zL2MJuFT1O@-yLU($s~v9RT?}qs@3}TvuA{>xLXWd!
z;b)F{EC$D>uxzu<7TMgfXC`!Xff1C13m+d_8o$5pNVj;u{r^9cLCu7(uR@ESo;upZ
z$~`Yjv`0vlYj5@Uw10no`pmbxYpJ-cqhm+=X?H_+7pbE<6R);%i-T5Q|NHkZ=l;IC
zH<dU*E1N;nVBhc8_fMZLzSzD0*!lDP{QUfDwu07aSeL&8%|IR4_(ZP3@V`l)uZzkb
z$gEOwk0$dLwtI|w=HJL&yY*#kyp~{P<iUdnK}S2-)&2q<CITw2w&mW|2w1@IKq3RQ
zW)9l3I=`p#v&TFeONR9G^NwD>K7D@OFV6h@e65fb4y!^|itxC)#3}2&e~{d>l3@)?
z62r_!W(zmH<t$DNzLRu6N!{L4Y5eo^bI_nRXkh|q9b)<WI4%*518bwVYXmN0G0nPi
z;PBz*R&Mc6xAXT;w5cpwv-MnNcHOU+>W_}~o(9GLOrxc{SF0!~`4_nJfADYUbZEVB
zzIN?a)rK4g5t)a9E^*F#tG_REXk^N{zi;jo&EQQXFM~iEns05%1WiD73adW@ZI}&P
z8wFZJWmWo$q0jx|)m5RO>37iH3(&Mo`TICf@|@-@IPv-fYo?^M2kbIyTc=$(e&^aP
zXNP(j#vqqC=a-k4t9wn+0JSsz{{9}kHp)~lb{A-R0#tf!Ogj2#`TV-1kB^QX>6gD>
z8P_eQ8?mp(611f2+L}n6_<eVNo``gDiL-zAqFI88B{F-h;WMKNToqg1+_QA}%kLTo
zIxr#S{5)HeoEriS%Y0`Ct%)!+N;|`0;imKb<>hwJ5R7T|u^EQRa~98(H_y7FvB;(K
zPz$H<7mrAnC%g~(#EvvO%)0PBy0P!rjFe+DTpYw7v`C*)QR?4W_t)y}t*y%b^K3wC
z)`ivm7#=))=y<G064czBX_UIiy<ZNrk_og3QPMcg;_$;RQSFRtYi1Uw6oHNe*p@RB
zRJ%wYl4bm;t>FFQcb&yAxp|Cw4DXre?^@`i)Rf%0=jG+)lQo0Y)codj7`BD2jS|$4
zvpHZW3ObC%bH1Ied!Nk0<l~@$(w0T8-Ao0~&&ft@&+~nHYU)hWY%%*kADXRKD=R7e
zSDtR4wC$Q-&6~ew$v$NZWcmWn?guSeOLTRKNqBv2Eoh^*ZuB;WfISt3YCbbSN0Ko7
z`}Z$oML?r*`nd^ax!`@F`)Yn3nsa}FBQu|@Rm<|_>Y(xRsal~$-mxw&PxcpPaxXGe
zn`OAH+3?u&i)PtxC9~g_WWTj&br$^o^4?zS=xsTQm7kt~_6)9x-7O|BFYn~!^x^Yo
z<@9rN=2(~a-P>Eex&D7$x&;p(U)qTYiXu7@6F{qYbfe9-<=k|-zAkoa$h?jd(x;Y!
zHs!Y;n~~@c8GbL^WnS7^RRx6uH}eY~9q9xu4fD7EJ0)Uc661!%!)&WUS09_I9e!ex
zDmUnm1<T@RJjHW8C#zlDo*&NuI^yHPLg$a)>WcnV3kp_F-tqYQ@4pUKpb=#!Sv%0w
z>8XI#VY<&9Ch0_O0(JkBj&>bA+|CbbBD;2rJw4Vd%_C<sBX)P0ttcO3)V3VYpP!zB
z4$fHY-Va)j-pO!BP{5wmj`MA`hLT>%|9`c=zQ6D9@8{qD@0WJ;_PnQ0PEJnj`}y&>
zyypx9MbP%?@O3die*eA<8i%nuy*=;lqTuCzG5hOe8T9q_LAlp!y58GQzdAZReiR(H
z>9B2xIKY>DhVy|<ALnNIjn{5H%F~$c(XuV=tQ2S?!^6Yvj2V}fa4Ng?C<HEYX-H-%
z`ughX&#%|xL9-p8RzvOYZ!VoeM=N_lJ+5D0Ul;%RQMfK{@2NXGi|+=B3kq_^zvJjQ
z*ub{9V@JnUvH1A-W&ZQ|EaX)EXPI<9I{N?r_YI&9(CgQy|Nj228MK50)ZdKV_5VGm
zkc`Zk3k#iJ-rsM3@ZdqvQ9SqdR`1@Ur=;{q>fMZI%UV{y*A$v5nG9Nr85b7^nwkXd
zy2`k?h_T@5DN(J^RYyRlEqwcy_UXyVAHROJMauuYvNG5t>q>|HzaPvEda=6<Qcei$
z*|X<Pki4KEXZ^E4_R^pCpPrq4J@4+Wt9OqI6yDmNFaP}9+{exQb__DL{O9M{djI;;
zH&tit)*gBLd3kqtfu>!&rfMBMGt<~}s#fTX*PwoM`Mt_^Y4f}ZCYeE?CeBQw)R;JV
zLBVh%#CB%g$x|eg@0lIzF*=rGdQ7I?M1PUtvc-naY|E#3swjcl-g~RRgVt=B<lMMm
zIpN6d+qWh8{O4M^etUEC<&BNX7Z<sL23kO)^8T&P4M#eKFE4Op1`UY4ySrPI&(%eV
zbGkj~gb3z;M=pwg)K?G}&@Cy?Qm_xWfAy(Q<>kLW9`|p_zTPLS?w8QBFluX7qS`7x
zPEOF0BGA%D7FO24jY+Pc*=Nw0NzC*0`~N9<P0@IAe!hIe#l`Nzrdd}E`scd3s3d`A
zQ1lDVv_^>S;J2yCx$l;+L;94;r(;&Nzf5GUN<cj|2G9V`$;s-$3mh1q8-mJ&bulx8
zm-~IZ7M*XHc!=fI_4V><UszdL4O31Cu<=T{tO{9qV`s5?;iDske)-BuN)y|ynUVw_
zIA`ZJ7^<BT{;kL4c{f09qWPq))3UCv+F1L$43w!3G%{!2-)9@R*y`VPJ|3PYkNfS@
zZf(hAVP^+5C64t<Kc92!*~G=}{h&JZ->=vDJ%(4-#YPLqb$7TVz2oRm{J%k8vEY0t
zt6*g?D?2-AhO);HG}7TYS<Nu*j6~29kB2tT&fEWA;?^q_v$Lob)F;ZlWug_fM&cQ0
zedgI&uI+rX92^`E-fq92_U_KkJ^%mxUbk*iN5>wwJBA*q#uvU={(i+E9TOoo(friU
z&(AZjtPljPA<xqeT_tkA?)O?_nR|Pyr$=qg3S8vExh`g>)9Y(%EejqT2!1GV;QICH
z&(6*k77%clYgKw<Z}s*=3p+Y?s4DO=X-OxCG4~d%S@fyVFpV!jR^m~cQRkj(O-)Up
zlU&~K`ONqE`T6vJe}2a7srb0^@ie_yuSqJN5t~vtUtU`J_`Lo9C#UuItN6{aIKL<F
z=clKjO_iV#m4i*JNBZUM_pDY@Qp#tT&hC=^_Pv9tg06wE!2yE>8UbgccE1e0e_7&j
zc28t>&(6f2jmKt~oVQS$XP7q6@LAf6?0JT2?^f$|pV+^(?r&Am^K-t-{pNyBwVJ9G
znsk1it%~QQBb8^qzP=9XZ>|hoEe2Zt;6LB)$M4_4OFSkTT;EmxUJf*k__?@O%G9Fl
zjYO}M=_QR+LBYbJ(>5K~Zpp9RD%t2`!L~MK<F#d*udU-*9({}<jQQ{XhcnhU-f`G@
zAv`Z3&SI+jlh&r@=FCe=K-)Es-oHP8OXlTGH9w2$emrD<Zd7qCYMJkBw@*(_f`;-e
zi=Kc^zT=Uz`LWVb4KyB=e{T=y_^sUA+bp-M`^>no^sA8I#Ow2|*;ZOaX7@9=#}v#r
zWH>e>{ajaG&&eib5wVTe{)z57%m3rD#OJ=P-x4aO2R`uxotqWaJ$<^kZPk|r>F4Kx
zmJICte`cmJXtn&^-Q~ijSs`uw@_x_H&7G?}L!Oh5FYQQ&An4SD0|y-LRX(3<@L9R5
zqr>K8toF(jlVfKLkI5JvyTkd{^0vg}gR)N<xtJI&o+Z|#Fc+|EJeelsxnxV--zu}*
zU;pF8K#l76_xAEg7#xt~b9Qo)uqx4bd1-0#^K)}WbRq=m>+5%wzyB9vVN>=-;%K+{
z^78lhK--b_)mDSnUaX7R855;1C|J39N3!Zvp#%OF%QpWMJbtI$v*DCOY_iYJ3w2%V
zrikpwk&m8nD=NF!aC*hLJw3-H;tz;dJQGS)`hV<bw>YTgE`ELvv@`PU?d^|Wz65od
zJm*+!{JaBnz|+p6)IM43w(Z-^laKW*3|$>|K1cA+%jNS|L~qw)WMowFpEoCRb6Vy~
z8$rQJ&374&f0%ACF4=aiE@kc8w6%XZJec3`-2Lx#Y{tQb%zpn~hAPc@)v&c})mF3h
zTg~22p3K+1XU&2I3dee--Iw{y1P$H9#l?Y+-kC6A!{t^sHnvMky~Q`BpFcNI*&Vcn
z<M+3>pw7hk+L$+wA2+wRvv1A5e(e7J_=f>5E^)K0e<`G~zKY6bmb;&^c5iNvBwx+v
zN!KNr&(t3@Tpsr+O1#14n8fE5q8BaNO4Ru7s((85_viEZppy3P?sDdUg-)!XBPiN<
z604(f?(LbGb#+zY$49Q9>CGD(lQS<ZIce-JyFPk*U((U8hYue*f(pW|+2SjMmtW&Z
z6%_m#e&=FS1zX3rIkguXHtQD)X1^89e%mO%=iR=Png5uQ9=_=Q;o+#DmiSw>YmX2p
zJGu2ry}7q{Ht0OEoSQ~rYa#^Iz(owG5xb)x@!WiS|D#=^mIV(Od}o{0YJIr7r*bm%
zuq#l3d1s!QuT-nEpy1EzcP=(r$Vepg+{{>O`#oc^;j`QY@*i%l-nY<DO|x%b%S<Qp
zRNlFdKTO&BZPwPZ{f5ud6TC9+J2hw?;Ob477OtwH(J^6yz#`Xf&~b7<{(L@fS@z~e
z@Z!mgVQZs~w#(O@I6K=MG?e)Gc>nXm{Ps0zH-5j}em~>Z7R}}U^VO=pyzr0p25rtM
za(B?-JHpr}nf&kBfA*Ck9JTR}HtuV_7IpiF=$*}?eII_F=wE0*n`yK2{p;+<m{&X5
z>$Bh5+f?0cA+LUYWpKM)?JtH4J3d~n-?Bb_zu)C$zCT|spU=U=v!v|ptsZ%MyE98}
z^n-S1$L*~$EPE5d@ag&a^LKU@S3aLxerUUql9JHFnaA1Q?RwUBOy<9Mdvtd0t^(~d
z+LIthpoEo+JI6VT>%}an{{HUF?B_FzZL_bfd2?&)>y<Av)@;>@*uW4UAOG&&UhgIr
z&KV||LZ;c*YMzyzpJVAf(<l|RFfeej+tbVb_DjRp$N9&GySV(b1??t#Yy8}=PbcQz
z`p&IpYgfP5e3S(}<frcMufo^Y)*4KXwU>MZ+DqTe#{1*-`u(7D<anh_65ia{sNyl<
z!OrXcvrIHW)#<vkbIkLX{HykHafvfDHSgGyB(5K~q~ztLAHRPy2W(7o1tm|=V4VzT
z1^yJx;1{>HYJ(O{6h1lvTCSoMy2|GBF3<+*sxL1<9rYD~i#Mg56nb>D+r3A^@X%gA
zB_*MqJ9HJKl7pBu8Q-gakIcTlW1oxCdbdR`ogMP_cA%3L-`v<3vM%Q5%4E=<h-Nlk
z(B=?OYOVYGt4G@W-OROZ*P@=CpAVWT2hBi$hS-lufToASL)=|lloZR?ESbdjgX_fq
zjNljk(!cI81;)hcP5l4p`TTmd;AK8byr=W=^Yicf`AquJvEJWTIvyQ3;NaFLb8>>B
zGbm(0V`}g2>^%QAyW;Ds(5UTsXCEEyW(2MBy0^DF%Tz>gV!_isu_G@Xv=%5AFjTx~
zGwRqQ6uvHIVeoQ4P<Ni+?nlGb)#1rKptX|zvb%F{n@O4FocQ_q`P|P6ezQz8&jeNV
zM{ms%^_^|@^v}=FN4iA23m+e|EpN&Z6r31-QkKzB{=$){-A8!8i&dtA#`8dvv_((3
z?%lfw8r%hC6c^6h40Ye%#oGP-61>c3=AmYGe$Xz~_}Z^mi=UfrxprxxGyAcAd47k7
zhgw1VJ6>Me`;GaRtBXogfxClEa}@7W-m<T+LT_!$J$<N^n^7ZjlS=S%zuK%$n+DVD
zYoKEZK*`X-fdRDS@6(f$=W`Nkjz~E7NEm`vHa|GfSn~SX+m#nIm6U|$-7(}(Qa-@1
zk-c{B<4d#Kh0iuiRD$vdXamEmE1ICvOFX`&asB%D!Mcs0BgVX^>(zX}TMkk)!!X&X
z?oS2i^nf#Apnd&)^7ekmdL$V^5z)>kt5p5%jo!-X9UXt{&o;l9S1IsZGTDe(-6=e1
z=f5wC>r=jpObpkn`tkzQJq}$RW?230&GKdI7(ui28CO?{%E`%r1`1Y%t_E!;ouuab
z>rA$sby-gvuXNF;C!T-){sk?jkBf`rYFAQvWcDt@F^B04OH-q_&EFkAdn0V07G^p2
z+%U_z0Uq^;-)Hk6<E`;Ao3A&%-`QCVnhWWbHZOX7j2G0t1TF3eUw^Oox74}0*3--0
z-VziLH~`ve(j}^0@%?UjY@n{7U?tyHg{F4tca{thTq>Kc?aN*3`}@f2Uaoeb!n-?*
z+t1Cl4qg*sC}o;;<Y>1zV~3P!){KITwW+tZWd8W^12k<28n^%z2D{7O7qc${O)Kmw
zeGNL+fr*J}MeJ^~-R18=-NK32o2=PZsxa)^a_wHh;k~YLr`GTP*R_7V{_%eK<A)D7
zPt%PSGS9oSywAGm35Ra<wkKCs2D7lT7QVcs3ToANi|Nj@mOTqvNz=^6Yf=2HXY<sL
zi~H>^eKA#1`gE(no%>OhLbvog?N3gi2#wsF21=NZkM~bDPWLm*zjp@I=l%4lNKRZt
z<j9AIhfmJ8w>Qtf_vFULWYEy%T<h{nn^I5L^Hi6=zX#g7Dxwk40216$m<*~;D%b8<
z%({_z!wHFE*&j2Pv#(58%Uf8v(xoBm>Z*^&<?A0jed-E2;^5X+?RoR&CHCFiSG#+2
z>xI*&x$EoeK|O@?^K3h%zTW)reRo%Bw?*M2hJbZ3me9`a-#d?8U0m**l-2F{-T3A}
zj7)O~)BK9eK9|Ck6DIusUz>7afuo3aSP$q_xz*w84IX}YcXu~vKmDfE(?Zr|IiPmW
zogIZgU$5W);`a9U#fR^K2HK`-iN3zRzC2Gme%~HY%X`nZJ0JD0ILN%#?cQ;9WAbs(
z@Z<xILI%*;l6QBN`rCer-Cd@;Z{I!@?`b^#{{Eo3qwnwTzP!KxfADKv`=Zy^bV0i{
z!sBaMLC5oH1TFEfxaQ&F^5iP0SN=D))_>yh*c}Cp2b<ZAOJ9XBF*ASs^eO4diHYYb
zx#aoVL9tT({oS80m;LvAK4-lqb~orGxQ#hCjXpg+{rygL&aEw-pw9T48-__oIu8E$
z_&D>{mdqX3eO+9hgcclEPE$|z(Mk5nO73~y-7oetA9O(5e7jnY`F6I#>V7I38Xf(1
zzdHW@{%&6^^!ok2-+c1&^5^DQIu{fe)c^U&Ui<r7>g{d0=XpWL&GAT?bSSxY6}-7&
zc(3+*EGYaor=2wcoriKyBJU+jYR}4+$AQJ?F8KazczbNd%L7u`Z?C`8nC=p@J^#L1
z?XNF)_EZWtfbwv~hlHA%nt4{=K)EmbI$wQ#J?OZhpU>xm7R_mguTu$M7xTYje)jcs
zhIw}^nwpwkTv>TJc=3E!7ngf1dijnf$3mE=F@`tD9e8|J>o}wCX7%-;QQ|qw0ic5*
zwq{)gtwxnL&wKLw{r-9_nW*fUcC}XT_x;w3i;KInv$!2pHm;AiH_N~G=T1S!T2%`R
z2~g+SxcpsA-QTa_ZoN`V3pA9J@>R-pwr*#>k+wEBv8S^u`ASh-_S>>#pHP{}8T;Q}
za{B8YmG%DXrqb78pf#DG&8U}_dasPys`a4a&+k1G4}-QfUB5nkwt4;}tI{k`hiqHU
zO;Fo(XYup8cY8rK`R25<poKQ`Yd(4IxZ>{O65|b?P)HI>NPO{I{qoCF8SS+lf3y`9
z6?NkG&3SaR8?=^hnoeZX^5yF8eKLh|pVxrT4CoM41`Qq0G)jGRYO1zQ^tPUPw$+<*
zZ*P0v=A5m|Ev|Ru`0;j7$y)fxg^7vj(AqjBCCILU?F}vGZ%$|1@n+W|4;Rp)fL~u<
zA7AL)&d?#O&bM#hKG5!oCnqP{f8Qf(U3OtlW%0eL*Seso0np*2Q#6I8&2knvG%}sH
zeSc}GH)wPabl5{*=lb+>b3jMd#^oM$adEkGL{_(BKf~T0!$Nt`T8=me>x8A!r&RJc
z7^j~rcyU1ybcE%Kz{M7|zf6?f`)cwm_LRS$SA6I-D8aqBpa`nutx8`V>65jdC(S3n
zRVwZ5tf{Kr(*(4`bkzLko$<H->jGMQvAIi0Nl6b>t=k<)d!X<^xbOHG+eVoK(GI5O
zQ#|h}fDS-8Gt+pAb~s<dG~MV)MyXy|*Vg<zbk!5IVdnHS-IH@HjemZ6n*8U-$0sKy
zHoxEh|DXJG)0&i9TQZ^h9t~@Mm0Vil`B`g~ijvZ&sKTELq72DCTnP*TVjes(M>pzy
zUViP?Wj*~1^-CPJ)`3og2CdgwzFhtG_I&l)Utjb~@AXQXzqqzmT04B*lA@=lKr1Ri
z6LO&4E=#?qFTWODcO>fc-QDG#lE!Kb>i+YdJfB~G?E3Za&r)3-E&)<Y9xdzKw03Ly
znyqKoZq;0~_2QbXC)aHCy?>-YV;$)DtEep*poN<|%HPXfT@_mW?8&wBpb`0+ni|jw
zG0+hspPrn2a(1?OeSLk^w>LlU6jc2ERSG(=M=Na2i#dAgrLV61`1w;)1~e`A?34SM
zV>3)0q_f{ProTDO;<N4AzW1++PQHG9VWIPtb+NPW?k)$N#BrokIP=z)m)y>`K<y*p
zxbXFHb6;Lw4w5=E6SO@0r2V}8X=kOR%=6C7G)@Q2=YS5Xir$v9epi-~lG1wB@+}I}
z7#`WA$LB2A@w!WCJ?N~#iRSrnW$*9#9_x_=?aPecku!JeG|;92<8;4CDxOY<S~zDK
zBr?g@d|<5p{_gMbg3AB@YOTuOftI#|=7vG%Lq0z@7t|^cTDc>2rR}_V|36Ip(x5h3
zc)5z_Bmn^dhqt%3f(nPiM>g-vx!c`^<L>V&?FJQ!ox<vm-o8D1XJ>KdjSUa&*+4M@
z8hryTah{?P2wGkA^6u{M-0b^0IyziT^fDbyj+vMq`;*x7liA&3&aW4rlh^AWka^%A
zKE*T6etZ7?XQ0gmRbOA(NZCwP^%kr9{nlJeFXl)8!w(-h*Wa&xZ@XvDo*$38_2;Z@
z>y<LCdZ_&J>FMd4)6dHprJOj>#w!i75j5tf5>#|`f1RR3!Heq8N1nT{W18IXaQ&xa
z$7cMnaPyhY7W6oRxr+6|^{DEv8y`PneB7Y5aP~Cjw`;eW$E8+2$-KTkKK;v!i=eYZ
zUteF(yrJ;18<)5qkHg7Hs{cR!0cG(;uH2w@{_AUNe?IEg|M2Nk)ylBL?EG>T*S+gZ
z5|-xP-Ztw?l#7d!>GUlMQQ6NUvc03Sy_xni^UBrMy$L-Rb^B&i^+(bBc9tK^F36tX
zGhyjsyl+!@=ek4Qf|wWD<_7OYr@P7(e}5MX8VcD}@)C5A%yr0#wQC|aHaWF&?fviv
zbc___oOn>(&&k7+^5Mb3DLRp#TGzHIDk_4;WqTxziyj_gJ+M(EzF*E(gwNGQ<x=r!
zn~ryBJu4d)iyE%m`fZw_8utc%jij}A1*&`mW%_p;M`!c@yZqIJ|IvPN8JV8`etvd7
z3%1SEbs{G*GP4D2$+$S(>OW|*D)scVJ-^@W2CX;+9T21HB?9U)JwDd^f8*5;phd*l
z*Lt>QUuQBXdlO;z@5kf+ahigHodPEOi)KG_tl>SvD92FlY%hMq&?4<|^6?yocMSS@
zrzTBLTB~Yt<KkuouPGWeKOVM&M!Q2-hkd-AzkhFt1?Yszef##A<lLC>?(Xj5x3^5I
zzP@sOc&PQ}=JfOSe(}j`!E3$F&zrk7`}!n}K&5@Z-&up&SrgS=b~Q-sV*8?yFl(!s
zMoDqdgH5$-EKR@sU3#JR{H{KvqaAi$$=?S$lx;^r;xeC^kM3OCobC@ge`}^uYFyT%
zm+r!G)fE*Nb{3~wmA^Z)%y;$-t5Pk{1`EyLWi{`1f}+W-SL(;lpMnAcpcSe&c9-j4
zT^(Mp`c+s^u+!t-jfRXK$(){z`V~_v9?gk3z<(f8rdH^GI`WE<U;5L#E`Zjs9q*F`
z9h_SK_pA8jWxmCI0exb+QgU+o|CB2p9%9Y9y2|zZJX_E<u691zMa9q0frh2)>+7So
z<@`L<u(SGm-ZcGqKhVzKU$56kZpje*^<GU$=~GVrCmFfqhM)Qt&#o*sU-<CD@(cY_
z*KGYRZ~cyG!U>6Fxvff4?vFfVcKQEP&DQM}(*;eO%(1VR<LBoGP4hG`GRIZTw93{6
zomLgTKF%cn-W(=YF3=*LMQ**PzP!8~Z;|-zF=)~3mdwkbDb$yjmj2cc1qHwdPynpm
z+mw82mblTeBWz)hW+W}2v-R7_lRN%H#%)%}CqYJTPI>Prddl_n_4Vh|<LekfH7IB@
zxcq+Y_oeGWOTj>U-zq<=1uyqYeR5(V=%Cy#QEgD~CuoVs$NuN>yUTjh&d$2HEjQY0
zx*o5?=Crf_o?C&AvsTLktshun!4%<j;P=VD*KW;R{~>N;<5df`Wg*N-Qprhe@Ajpw
zt<6|_H^1kn<cy@AkGyxyStC9!5<K&1%l3O!-Q41O0{U?_pxpWS+1dD4fkqZ;d*p1Z
zOb)m6FR%Xoj)jdaXnmY5s6Pl<tpATYn%6w<j>q1rub}#yolhp={k^>|okCATLqNw2
z1{G8~m>he;HjPQMV9N{JlT7P%6qqU-m=?b0x81m2WrN<A+S&x04O(A57cID}t>eDr
z;memNn^?I$rf6K;|LXbn;^%&#RXCvSJ|>w$%<O!iNm|fS8BSrf2_`R%k8@kvR(xpq
z_xHDP(h&|PCnwM*`-mL{3qcL|6Z@a^i5)p<v24mpHW>!Fw&!QpY<*o(CQ$on$7>6<
zi_`ued~<X2&&U1tGt6?O0vEX`W?x&QU;A!v^>@&12*zn=Kuaq?W7W(3WP_IZoc#3k
z^val>pSag=%)YJ%I=!y!{k>kpHaXiWgMtSRi(I?E#jjFTQWDbKsSFxJjL5F9dN}cu
z+YH-kv$L~I4@Y_TN|~+*S$WBO{_l)KuR#aoWn54IHMyRio(?*4kC~mXW@V^vQ&SUY
zi^1_eS!L%oo-;FzpDztd5)>4im|lF`#Js(L=ida52OkS8N?u%;&ez`WUvYP3u==|@
zJDpdDt@W6w#0uIteq}|V{@y3CS5^ehw5c@8y}d0JS{Ew)pLlxuE#qTTnA@1AuHCvk
zuH?^bg-v=}zo%Hdxw-kd_uTJ)zuyNf?Km|>)A`Gbi=gum;^N{!6PRUhZ|&W1KEcC-
z<MFZH$ET)hgT~N9LtI^4j(B>%-T6>(^VV<AFL18mmrmwmyguLVh*){)&L{V08l`fr
zjo$t*{NaZ~t=ynh#j(4~0(Y0?K07=6_}1*}8Ta;BHnZ_oJyZ~B2wv_7+68}YO=PoS
z8>lAsnQ64QORKY^!$Si++x6wh#;%jR55xri9ML{<HFZ+``}h0**S#$Gk#Xp?jAfC4
zX;ug*r_L};1|1@FdwYKRtu2}QpzW!k(WX8b%ZSZsyq!X-NiQ!g1vL#o+m@R{tRBrs
z-Y|Eo*$v+B)9W+8zq<>nB0=ZhZOsa00Nr4*G3jW=tt~I@=T(1s5%}lNpB*(njX=8<
zBQ`XAd3kwr$;+TsVQYT{X)gpd6vNg=eSEw9J`-p;CTQZ}(UH#o)lXeqT%KsY<LEGK
zILg(&e)-<&?~mTUKfkm1xkdRqnc3#~plK%-&|wVc8fGZF^A$Zg@$t@U7FJf!yi^;n
zwA)mz&>MTJ&Dr_oQht1RSn}cmqrU$7>3YH6-`xeB&DSZaEe1NPG4t{={k=!#b##EX
zxN&r77p_^f=@{>KF`3gBa_{b%dUm$?=E~1$$NFSJhx##sPIvhE>9qb~8F}q6omnQC
zPP@zW&GUbpn`;eQSMm4z{q$>VB0)7Gv<tfT!=Jyuzk~Me-mCv#%h1Lrd+O`!>!1e1
zoTnm!f)hVH?rX4Mn<0CP>5H*Hf6&3G>N7`EyFQrZ-a694DSTr~rf@sIJm2Q`hxzSG
z*w%l)SKaSDT@SSSWRkkS->Q(6D<U>36+Svr`Ox;={r%?;G%|x4X}VEc5)L-8f{t1K
z^6F~$zrVl#KVEWmd*0ooi;G-AM@{Vi_p2MUaQejs#U4Y@f{OS1e!tsa%Lux1=h7X+
zvu7VS2h6K1&VDPs)$G=5yR&C!8mD_q*OP5*Yy_|Ovo6p3^!&X4zM7vE|NnlE&tmy#
z|MTzn`yao4RjvN^1~f&{DXcz8%~uL^WYEJyt@EqzfllQxOkz3Q&j0+~-Q5gF`ebKU
zeSP(J^-fT~I8M~J`0+fkBWX735?((*`EssRsgi}o3^rb=2ag{&FJJy#d*?UM>A?5)
zTHF8op=?|ACE?kbnL3f1TEujt46?3hBpvVjd#9LR-p=ReC-6}X44~<o$jxb>0&s?*
z>GWM49UZ%O>|r%js=6B`@c#b(^9_v5plT`g^fXqP%SZOSYLx9)+p{iUAyek1C5aCX
zG#=>`Zhv`sdGJ!Nsp9kgKRC$z{E9hf;oH&f>|fr3f`#^56?UbswB>Y<-j<{2KhK6o
z+KeYYK3>kQhU0_R;fJl<;wxfzn}N;)zOo`v#<psTUF|Q>2}R%D-28k8d|(D>)#19h
zy<82`r;BTcuS<cIJL?6(>(e=teVCJd-aUU)e_P`7O6}7k6F~<OE^_VOlz-pODD~8l
zCRXl$g$XD2fBn&N9dvYK`1-gXzkY41{G66{X2!wq_v`=P*~wY{{+=&rC*rNG+Iy?N
zm*x3S*So7-eMeVGNojo>=vD)@X@+UaAEXWVZzT2nyeTpHPp5M*!{3;|WgRQBA(y7i
zGR@vp{M?UAOy|d`b+un!1ls+0zzjN*{msqImNh>NKqdC^<L#g!q&1P7)3*ghWlNgp
z$$$=?K0nVk<-~->sI6IhHx+eubevfIe3jUddk$g?I46i^a6e%(lT6;UYXy(styfo9
zM{Z1F1sz}W=t!qazue!cG3FW?9iWq8K?y%-saN36qSUObtB#(UstwxYz5o9|>z9|7
zRzKSb8n#b5+6B6Qgq>gR!Ryzj&&{>gjNSD`TXS|tM@I{PcyH1=#xI+%o#Xynb3o!V
z^9G(Dhb4mVMeRQ7n%Y^gW^2joYrGAhJtr!jlfsMcx(ml0?-XutWM(%?I>G@uwPdbU
zsf1;b%F3XnRqu9&uZvN1?~?(oZus%>v2pFMl2zgB|9PEKRZ>!Vcsf+KV`YPEfBuGR
z^H{r>{~cI=Y=)4YmUl~X59o3x&}Je&S*t(KU*FwdKmXp|YS5|vt3p@1P0<L<xV=pm
z<o~y~w(7*}_#l5+(cw@FC+N%?Rqtsj-{0K@7sL4i(|wAQdn_5hY`RwWxu^Et0igvf
zA6PHA3#aZr4;sS+_0K1%a(hlzOa1lbB`B*oHnUw^6{@|+t+sCQ=KcHa8yg#+oO>(L
z3|c^AV{6+f^tAWac}*pyiF20iER5xT6t>{i=CX+F_wDa~3%+>0=#z_tzuxTc`*&4-
zR$Cjj_0eH|`z2XdSLNK@HB~cs8EE*NPu9w1ciG!{(%?x|<20V1pPnwxySs~pg@vI{
zG#b?FaJgr*j;rH7!^&Pm&UsreM@h3>VO_?sJSJqB`MrhC?I)(`N>9}eSJTjVQNQmT
zXb1YokDy}#Q%+85ov7?S$sp0mEbq>VFE1~9PSLnneBQSFoy^lyQx^v>_XC}FIZ4%f
z-zUzUpbB<tjFHFDsNF2blb_UE9b;M>8#u26bd(+o8`~s{!bkNrUw+$&>&JmtM7+MX
zw&ue@cF+<k&`}#3laF6o<~v(FuK4AprJzCf)#2+WDL68votfcyb=CPC4KYE%iRK+)
zL6y#bzunFUjWtDVOfsza@SwhQ-PUP0H>XdwtuC{%wFTW7@a@gbko9qMnc4Y3TfFxh
z$khFMxg0dtx;lJ4=yIe#fB%B|sV+(!;l9Pn>jMf34CY#uGA$^7e^14G+MA{4ZOh(B
z)Ya92j${B$vaH|#Z`PHS!JDhTzN$C1Pu!Ms^TVf4Nk=*aw`5=E<LCdoPEkooDDXx6
zSMv*BUS9tA{X1xc;@pBCN(T=f1f3W4=jZ3iiq33IOiVwX&CUm{EdKZJ-+u`{{rG)8
zprg2ob(ugT6??0{fBJY_URwrqT(yDT%x7wgil3U8i|fa!6hAw2XGbAu5oFMEzq#V`
zKu5DiZOLfdzTF&DX5QOV3A#Z<(l{+(kxOU&p6fdbA4^!5>6E{_bMW-(ZqT|ngTzDO
z5z2ys6R&qK+Zo9HNE&noH|S^{P|m)-?r;C`d;9C>@2mYCv?@dsbb#+n<8;ve4A9z{
z@9*xO=V@0|RQ&Vn^?Hy?MKpsx+)SSjx;;rNY|V*dz0#2IWZbH-kLleE4QmyZC-U-d
z|NVSE-z5K@OnyFST=&$@;^&~WHd)PgQO3nZ=WXBL*qFS;f4&?8KR>@#_&T5S^K3N(
z797xC-`vr0;<blvht!f!$2d4R6f7)eY`<5<{b1+A57*X4gJ!XpdQTS;5;_Fh8Wgw3
zB6ztUXn(o?+~@D#gI1$T*wt7Zleo0h+kLW{ug7dNUC=@EZ@1s)`@RcQ5WaS0{KyMh
zH7BMYcV<uJXHZo-zw%6Yd~NIe`SRc2-DO^oc6L_D>ubF=KR?-)zS)+08#DxPdYUft
zf~v2tK+6>xA~&bC2J3>FOTG&lez1WKQGI-@_weDv%@rRX*>3%RsFfR(;LgrAPhPt;
zbahzC%S){#FE81azB$$_4LYnsJA56B%;|?efBl-I<~!@dy0(svjy;M>pxga?nD4Pq
zVrD$xwqW_zvj54)k1!u&40STsWxo}joxfw-mPA%Iwo98*z1M6ld3R@K%+4aq?~&P{
zqZ4eazk&AogF5XqjnhB9+x_0|bChY;74YIG&=#j_YojkO^%kF`=KHE^73c;`X{Aj)
zN%y!+7&kGbGBW?k>iJmn2C-;mUUv7-ErE;OHdcSnyR^*L8+1e2?S2_cr)z5>U*6iP
zeM|y0oB8bQ?D^7sjGIzUI=#8E@#Wp!=IVZP5}urx2re@}P4Teox+yW4`37sswrg@G
z$5I%-F}kkVD(%#-$@J&viy8I1;wpNR?5e(K*xA{Ejt~dklCsQ4vZ=Wlv^zKXSdXHH
zM#m&oZ=>8>CeKeicF()B<KP5E=Mxi^*?Xl-nI1elJG-T&#pbi&){c%2iwu!smC7T3
zERubE7HAf5SqSdhsqp{*i=PvD56CVM_|VU2Xk+uryW;=9TF`A>GmTPze7zq3^TT0&
z(82rpd%ub;cI#b~d3l*c@|@?aY-~oEmsHqzBo1uNz7AUBT>AQ&YSfmDbAlH@7eM}3
zxdhtODV@L~Ang7?J9TZXz%Faef}3{N=47qidTfUG!r${sZDBiK0^@d-ygURNTD`oS
z-?dvT_3y8*paVWY#}}PcpAR~Z?)k^+Z*LM0HnG-xy&68pqL8WK(b4YDN5$g>bfdO}
z$Jq%A3U=<e2Oh)V`BHf>QFZOs*G~8AS^qGpFz;gO?lo*(yS4oP^xx6h%b9`?e(d;J
zzw=l(sF(#U(EIV@hk%TX&s42Y3EQeK!HXxG=ifV1{eJKBKXVL|-9S4Yr)q^B+M?G5
zN{Kz@W(rKnutADzw`L#keXzmsS^fj70O8{fQ@9<Jl6xwzbY;k#)GB;oH=Ub<Lm_ro
z38=sXEicW!rUSa;U|a6(8Rq$N*Vn}+zqzrI0koZUj%D#A`}(?9S691#d~~!^NVUtl
zOIJxrN$`NQw7}Y};SJjVAHP)l@>KNKPkjZkzI`42zZAuee6*-~9Gt$}VNJwF(8NhA
zm#C8891GCa576O;r>E;P3qUsaI6pel3912(_sM4If|~i^k9v7K^iNAnW_0&{xc)$@
z0_z3g3WhZd8?L=GKlX>mhx2VThXG`8;o(ZZxmKW^`4c7xJl}R@WpMDmno6yZ6#{7&
zIzb)vjy+m}b_$O;np7L~9AX|QeP9o`5cT`V<Hx@aR3wDHc>dJEdivB)lAoTQZePA!
z-L*?(Va?A^pfy`Mk()GlKzAtZ==#8TF{)Zh<G;}HSazwl`p&qxNh+Q@%HB#DCLJj#
z4|H*HQ4#9Y(O$V{Be>!|a#aM>N_SD=VFlf`NbP%Rd$wM?n_BT$GFgvVJ+@?b%K=gG
zyRS=^zCYl2zj+_?g>r^pd=B>yDzAT4?UB3I_Y2p$&&PrsY7R>zZ;Q(Qylc12yCvvK
zf62wk+JY{mM7}fki23=Da`6YI3$Umb{OTwI8Z|Ju0=|S$DNzG-%|M5T1Xo8sNAeky
z<THWEXKa(tq$Zy+O7@wi-q)GlGc&zsr{u3$={=pVmnrQJIyU3Su^Ata&0t=5|6peF
z8I@$8CiTAV4fh!v43Axx*nCZ*_%P_E*UWi^J}v+F7<vqIdkl~F8h-CJyxn8ydN+32
zhm;=8w4TVk9!u$8_xGL@{@rW1Y`@{Nm>2)$U-uZ^yc2RK{G7z)(-NDHN)&5<dB0!2
z$I#ASMtI?VW`<;+=;Sl+F1&v@&G4Civd_`u-`sz^e^9uvUT=+I+M(pzN>|GT=7|X0
zmx*pva}Zmglh6>~YdEi8_Os~(X$|pTr+Cj4=NcaNEqai<z<Ppg28RgqzRE{CI{wHq
z9D%MN<oGKqk^G7&x$%5M%mKHArG{!#qq)WQB{5bv@EzEj{^mVj3F9s1pAFsz^pDNB
zp{JerLyNJ7Q-tXrdsD%mzvBNlCM*3aXDDEM(fIiP%RL9A9#kj$^j-Kn@u4-t0SV~s
zw%V6!*dLyYT76ym2U8J)ScUqR;%OJox5`^DJMs9q%NFqfMh?bDF5BKS9Jl8a^Z0z^
z8^?cHhS?37{yEk>Jf-C~V?zEjEI64wpSM8n#pfH_UL43haQ^O7-sSHV!~}Rk`qtb(
znEZ7Q+Yi@e?;HMSzd!SMjn+Dwe}5H~x^@eeh`1hokuk?`X8+<85kr>Lo@1<wdGc0d
z^a%I%J)Wl_-x)0=YCNmsVn^Ve3!T3Yc^p}g5MwB%EWzUFIK@?PimT1%<144{d;hKa
z?)ur+ZC;e`n}7aq$^BiWtNrG#s(xQtn$_6bu<}6N0gDH}s`A-$3tjjhNp&pIJHUB>
zC!vb#U2DCCTFN3jcAaDOrQ6-jSf(+}U;4S~isQVpwalJ#8|n{M2X5|JC_82S)ExF4
zM>qAfZ<Zh2Z>{$;J7Al@eW~*9ONUEQlh_^|`^56fo86r8`^=k{cYdh4<TzWH?U9s=
zRQXm8bB2GiFS66TFF)07=-2t1!sPM4e#W}3+LvbRO04j|B(Aw*T}eROw{^|y8L}CR
ziyjtd$Un{IZ7sS!*>4@Q8l&*#%DXQdF6F(n-XO9=(BjOQxaEE&3NKR@sDF4YH2t*`
z`|2(Of$WYY^OrGh^OdbN+TL^Mc+6XFA4WItsP$}%!h6GBFnGRZ=4XEY__6P`jHt!U
z3M;}@Vw?GAOC+ReXX*O$&fs~&_lEP2)DB(y#N~UHGSzF|d{i)wj?1pRWi6Yn>oC_?
zX2;^07RFrvXWUdjpnl-?q<P2hYBPV}<dA2W?H#o6yT_!4buVwdI#m$0pgBdxZkivz
z(Pq_|6{QQC(;|XP12l_2ZcM*h^zbg!W95r?%-v!=*JGEM1jn27wD&1`S9|%+<Vb$l
zUh-?x65EE&4X1-5s;@cf-|7A!UcZ*{UEmAdmmv!bqIcT9PM9jlmc(ngdFK!I^`;uK
zua|u8{1GuF&ZNiwHEYSj3Y%M5uHGdR(iYnGP2atLrtT@%%<on!8cc)OrzWMHR@uw*
zh1o4#?<@1W^}Qjc=~DkbzGmOUp|H|;1Al=Z@6%OUlh!f!9qK*2y$BSrsV^Gp*<U%b
z^FPq8O4%toe@U6;yU6VKRuQQ!jXT&BR;u3%{vq-~aD!=qogLd#K?dXBtCy%Z2u||z
z+}_6^XSn&~>bO_0CiNUP{BdN8bE977CcXnT#?jaAl`<$~a@Xjb4G+Gcmud5i$9Q%;
zzZJte#=PlYnCuvzhtx<ek<Qi%m7Y{}u6o&Jz5`~H%8uXu$K!dpky~|dYrRT+Pk-pX
z&1{dtWIfFpjz2p8ujz{ngO1T=*R-^^vFr?!*0IR3raclaWqat&P`g~NBq2?nSw8)+
zite%b&E}VW7Tns%7pA#&o|~-uXZ8!S?=61|M6~UfabvW8{wTL;x0d*_$9xHoqfW=G
z%6>UGJxjCn%Sv|HC7*b87dIdLx~Mti*yRtBf4m!~H{RPiZDsH2=yiYleXdWvbu1@=
z^8?R^tx>0AkA9k7{jELi(Uxr4B}RQL_k7>}bKZ0LK||FgN6xuD%DTnD>#Z_YJ9{p4
z6kq`%h9f%96gZkd$X%iXBqG4#D43|g24T9kG)jQg35W=?B|)?VxS%LpEjo$q`8S1s
zA(yfz)oHBs6)&59Mfs_vg|J!Br7vR3F1ud}@C>cp7c|RvowsVHv!}TCEQQQ@OZqSW
zZx0KeW>lF!C*_Y@^z*ID6kY~+{twNz^iq83<Y}x@tC$(I_Niw3lDC(OfAlR`I>~N|
zN#LaoYVT%Fn)lS%C^6gnm&c`K-#gdMfAe|#OuO$c@A+%-l7DAU+;L9*7BqcJ-;&vv
zcfa-Z{5@qy{G`5VCT6zBPpj@#f0^g`JN41sIVT_AT$aMsG);&tX~V8hwl&>57!PnC
z_|G7I-Jvh>K%(xR<@^PlTPyc%J19P}j?IGY{;V9X57nOkZyu}6@0o7RbA*kLHO!hh
z@Y4QA)!*;&emPdI8Tod?7TtXUzup~KdVul3%OBe&x=-2peD;C9mtS@`?U<>1y!`Tw
z_@Fxde4~%z%->TF*&Vv?y0>BS$JQD0ldo)VQu<?Z_q&(c8=VL9KWr=HE||aMtI?<R
ztarHYwCS_e3tbSm;D5o>Gj&T~S6!vo|EpXnoK0ea(QjE-YfnD2W`pA28C&i?=+{VA
z@rsiEacI?^PrGOOu1%|2-H>_UTiUlswKtxgcR%u<;oT~@y7ofY+cm95YgkX5*%#57
z^3U+^-iv4MWj@{$);9UW-lM-~-Q45z!*<fCw`;3yKCpgZf1uzwd)CcAC2!3h@YHng
z%HC|gbgeZ%dtLs~-+T0pH$ON2_>bYtf&1>svnSSPeY<z$-T?z+>2B%j<`2IrEMxaG
znAi0iZ2t4_(PEXi{5tF&rQiFOv>y3xv%=_Nv&0(BJsde2E3=Q)7SCET-*Z-LzU@}S
zh9!qpZS`Mnev%}sn>jltY0|aHUsgZa7pd}6v8F$}=%dum8LFidOO`HKdGhJk#wEX1
zODCI%ymw`eh%5O2<kMcSOD3MDFWv9{^e)3Qe!5BAB_YrE-oH5S6<<_I+h&n2BUZtA
zf4#-*;})Lhy<X{OPFiwab+1;Yzo)EftX}4}B|km+W4)F{d#q;9o5S=+?eRCMIi10t
zdwr&tv@c0k_4PCAoxUac((jIt^kS9RwP!Xk@65ih!t3!#Z^D*<OB$ZDpMH9&@si0i
zf2q}+sh3{Q(GA{Nrk43x<@A)0c{3-yo3bVJ*}sW9Hf{=alfQg<>!E-yIm?|=yviOQ
z^(L3>Tax*sZ1(P@*}qTSoiN`kO6R4gXS{Z}|I2xWKc?j}1PA|e);{rC;%7m>>ZE60
zyJ9a*ySi$>QBGI>^lyd#nr7BUF3OqPsQj_^c0FGS(?8|n^V$cbr@S+?l`(hA&-PL<
zj_&w=>U)k>oUjjT-^495?|j|EysuVnk`-e*bN@60KJK6WnxA{v9ZOe!jjJy_GSRl-
z;19Xgw#5$@e~fGV+90QS?)L%a1C9q259~Uy^5w1kM&?5C*YSmqQ=5INwm+8N-`(){
zWB*AWsrT77LNn%BzHiyRU5weD@z=xs-|Cs)^hnP%UtaQ$_XhV3?Yp-HCHT)UpF48b
zP5hz!fz`&9HLMX6p4f4fyV|zDKbX7r%4?Sgy)ic|9xy*BO-Q)(c=mxgA@iaOr{9}d
z(SP*(ERTBGJEHkJZfCDA$(gq<JWBs^u&s~JQ#+v<tbZ)8e$}6Ju6^kjzteYatGAx@
z+4TMR-HG#={GUJBYkl?M(F5Cq_w+aZHZH6R*f;krOH4a|<BPdpz9&uEU%acH{Xt;O
z`8!Q^`{rdow=F8-J$t>-t?i@c{(lbU53GL}e@vM4T3}DE9>Xs4!0_EmtTu(DS2ISN
zPrk1C^48sz|Js%m`7&Nx@Ah}<#IHL<-tSxTymVG)eMncm$UR0Q;UB?z=jxYR?0$CW
z_>v#8Md$xl@2rshSX-xVqUgTxEstK|TIml_yY`7lOxpSFtVo>~_x$&T2R3W}{@1JA
z_a?BI@ll6zb&A~_>!so!)ed}b6g|MK>G$yIZ?1R$wC?P-`@s3(cFe5X505>V`LbsD
zlzTsSR_s5w|4L?#!n_%eYPNS8@%R>6Jv-Fz`*3gZLvhugyOw=EeW~)h&m6r2F;R&M
zmp&BQTZ!B8ng7{#zNOauXswp`_xi@nvrmFoOjq3bq3e}ut1tUAX6e*J=R98Co3s1B
zs8O3dFWV#GofYoOKkR#1Gry!RW@g21^Q-4hO{vvgP~BA@vF`7laHW3<QE3Y2n68~|
z{m&~goBOuG*2xdgp0Hlx^RM@lSL$r%WaaAjBL91`uEa}PvfkmS@oaW9m3yE1I!5VV
z@{Wb^2A9@*e1AKqbF*q*MEq3F-|T<p7P#u|op_-7$L<x%|CkG0ygr{%d!n~x?!nY`
z<p+*?=1t+&m)OB>vUIx8%U)%kBNnS_TDmPGVk^1L<DN0^XL5UZEOHwAnKMSu_kQ4i
zvf_!;n;9_?#=%v4gZJ2e6#ue!=a+Zu9M=79<UjLf;WQDoi}K~Fz9&wavqeb1KdCw8
zai#I>$Q>8Yh%xhTH#z*|&@a^w*2dDIN~c7$<(j{R{ra7vcQ5+DwMnx|i>EKRRHbQf
z@r>QdfMxGgKb^k*I^)iP7gNvue%$u5R`pe*(5EN2B^f5wO^8|7UC^BDwrfi;=ghEa
z>Sx}3_PJ;M$aYHQ+`6{!_nw}<@M_P8PbEj^AKagucK*@Z%KLZzYDd@Kv^jIe@960&
zs{;FXxBjwzXj)RrU8tJdXg`1P*&I(Uy&c{C=a*PYuYbMxs==1W5`X3foZ7ob^q*c-
z-}al173V+e&z`h*W>IZsTK@AdA}pOf^(#9H$FfBH}JL*BF2`Ke{^_GJJ19Wv8I
ze98XS!r<bf>z|)5y}?s{aBJJT1KXA`K1;pCHvN6C&~@#J@85geVy@tjUt_<-XAA4j
z&#AsG<wvR~M_DB<yi}Fh>oYggXX3mwIWw2A{%GH);IDq?XiV4HqUS%$+Ll}{WPT|1
zV2(;O`-7iGMr(DR-uk8c@`UkXmx%1ox3>6L{19mMEl$*quJ1M2oO$5ufv=xBji)m|
ze>6E_o&UwYgSVeIPcl1q!kkn7bNNKmp6ETi{~dqm?unJ--jk5VnrwEueTVh^4-XZM
z&*r$Q7Hw}j)PE@Vw#HwD6*CM3ny*eZUy~?SZ24)gRcmBlbL0Lg9se_bMgO=^%DicB
z`qOFqP8f-26bU&9pE-Q|j%?(!UthkImq|uyP2^|r-2KqjY;ptRGw+%!U0$!3Mv1)C
z<o$EMh3#;ioZY+47k?FZM(+IVwzvJqjKdq#%F5K{zl~Bze#q{>Jf`BFg|51@<8twU
zTOZ=Wl-2}&suWFQ$PZL{CuVYHU##kkn_QmLAF@5to9W#ecc{qrEo-ts&7EiJ2R2T6
z_VD}{{?m81-#HT(Ww}!Ond%+aI)j4CTBnRTyoA`2F3g+JAT4cscFoF}9rcG|``nqr
z(+;tRO*k%*uwY5^*)@0CGctee);oIayzq|$|6Nn7w77g&&P}|s_d(jU^9pu5n0|7-
zJO2N0@BQYYcZasAM8DISdd+*Qh0uJS;0M(e@9!Nc`Ny*F$$7VD-o0DezBaV(+#u|A
zB-$&v)?wDxE@`#k1L<Lp^)7kM==s0k#?>?Ha@5|PGLn2hRn}wK_Mm<D9&9^sHFxU9
zCAyQnEzX~N`TmW_d+8l|?iVUnc&P7D|64JsF(b`e<!<xt?&%NMzQ%^?J8*mMy&KMa
zScP|kknv2_ugvo{&p0Ww#V}EU?Jb*ev_`!2TQ<+x%>9q|vP3_8yyV!-n{(65A}`H7
zur!nBQumMcp9i=$?PQME`L5>~{@~a%wjVk7AD;YBZG6-=`qJYAU00Y5xy+T0E$cs!
z{o{Bxdri@ag)f)Z>+ZYpb;g!ICxY(u-dwX#eWiAuOxdC7w{NPN=X?k?j*fbDPs@P)
z{Y@#2m*;nIm(QAccJ-3DRT{_FP0BmU+q(YI<ho;`J&wlG>M!+bcpjPlaC#H}NSK?C
z<(%<m&9hHFvJ}i)!hQMg{Y766oy}*O?e&ZCp7oX^XX1|Ez7z7``AS~Bqpxe`dU0i>
zsm1^JbMT6DmHxfmcX{V?S|2{MZl%gvu1`^lFOw2lmRNsCPtcB8pC6pl6F$3pl46zl
zJ?VWi@hQJz_N=`->yd2#_nNq=cP_emR54iiT>5ow!?B0CKUmgH{<xL@`SbUy1;ZKs
z_pJYFA@p8wZomHo$p-!7n?27i;eOR-AkbXaI<t3ziLm(!wuiZUHsl;%;&bIabN2^T
zeXEC$zpyd-=7p!Goqo!FHLj^XM(?0~$oY#8x~f^`M-)rh&4~S>bK~hGK8ADt-}Bm~
z1A{+44Vbpn=Gm$8%pH9HJI)?XFDhH|^KMGn!o)w9lvUT|yEkj^UNz(9^hs?G9v_}}
z=&3EE#LcFvm2CEhc;8ijFwEqCY@E4SrrE#BLU6w8K9TFECx@6_`y5nr_w42mg(c@7
z*b9B?_PgKXZ@0_lXubG+%Rj9((mSL-%-r<3@!=2tK5s3}B{OZ+R+dltpSEZDC*~TN
zC$>|U1Zt+0O<r>T!>?+A`Z77@bH2hWZgt)MxUER0K+$e?z1ctIj2Sccs>koC5wKI<
zHL-BU4%YdKH7uFD9>%j3Ugp-AeSCjw1COqf<P-Ntn_u}Kc^WI9KmBoIZE}(K2Cab3
zla-DFEZH)d$G%;5us7B=S5)J;C1a{1Q!kWPvrfZQcpul#k^}AjJ*?M^Ht$_M*F9kW
z5hHU@f9rzTou=JKKYKHVcFzv_rM-DgpuiNBm@^=$fY>y~>XdWp+XL+me|V@K@b@Z1
zla)F|4%B0Nn8XL}#Bp?SFdIVpa~euKM?n2Lj-~}21_I3xrqUq?9-1lqQ6HF`=AROK
S`yvAa1B0ilpUXO@geCwFN#}Y1
diff --git a/Jupyter_Notebooks/get_era5_forecasts.sh b/Jupyter_Notebooks/get_era5_forecasts.sh
new file mode 100644
index 00000000..f01bc770
--- /dev/null
+++ b/Jupyter_Notebooks/get_era5_forecasts.sh
@@ -0,0 +1,61 @@
+#!/usr/bin/env bash
+
+yr=$1
+
+indir=/p/fastdata/slmet/slmet111/met_data/ecmwf/era5/grib/${yr}
+indir_ref=/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/savp/20210901T090059_gong1_savp_cv12/
+outdir=/p/scratch/deepacf/deeprain/video_prediction_shared_folder/era5_forecast_ref/${yr}
+
+if [[ ! -d "${outdir}" ]]; then
+ mkdir ${outdir}
+fi
+
+hh_s=6
+hh_e=12
+
+declare -a hh_list=(18)
+
+for hh in ${hh_list[@]}; do
+ hh0=$(printf "%02d" ${hh})
+
+ # slice and retrieve 2m temperature from forecast-files
+ for mm in {01..12};
+ do sf_files=(`ls ${indir}/${mm}/fc_${hh0}/*_{6..12}_sf_fc.grb`);
+ for sf_file in ${sf_files[*]};
+ do newfile=`basename "${sf_file}"`
+ newfile="${outdir}/${newfile/.grb/.nc}"
+ echo "Processing file '${newfile}'"
+ cdo --eccodes -f nc copy -selname,2t -sellonlatbox,0.,27.3,38.4,54.9 ${sf_file} ${newfile}
+ done
+ done
+
+ date1=`date -d "${yr}0101" '+%Y%m%d'`
+ date2=`date -d "${yr}0101 ${hh0}" '+%Y-%m-%d %H:%M:00'`
+ date_end=`date -d "{yr}1231 ${hh0}" '+%Y%m%d'`
+
+ while [[ "$date1" -le "$date_end" ]]; do
+ outfile=${outdir}/${date1}${hh0}_allfc.nc
+ echo "Merging forecats for run at ${date2}..."
+ cdo mergetime ${outdir}/${date1}_${hh}00_*.nc ${outfile}
+ echo "Manipulate attributes and dimensions..."
+ ncrename -O -v 2t,2t_era5_fcst ${outfile} ${outfile}
+ ncap2 -O -s init_time=0 ${outfile} ${outfile}
+ ncatted -O -a calendar,init_time,a,c,"proleptic_gregorian" -a units,init_time,a,c,"hours since ${date2}" ${outfile} ${outfile}
+ ncrename -O -d time,fcst_hour -v time,fcst_hour ${outfile} ${outfile}
+ ncap2 -O -s 'fcst_hour=int(fcst_hour)' ${outfile} ${outfile}
+
+ # add reference data as possible
+ patt=${indir_ref}/vfp_date_${date1}${hh0}_sample_ind*.nc
+ if ls $patt 1>/dev/null 2<&1; then
+ file_src=`ls $patt`
+ ncks -d fcst_hour,5,11 -v 2t_ref,2t_persistence_fcst -A ${file_src} ${outfile}
+ else
+ echo "Could not find reference data for ${date2}..."
+ mv ${outfile} "${oufile/.nc/_ref_miss.nc}"
+ fi
+
+ date1="$(date -u --date="$date1 tomorrow" '+%Y%m%d')"
+ date2="$(date -u --date="$date2 tomorrow" '+%Y-%m-%d %H:%M:00')"
+ done
+done
+
diff --git a/Jupyter_Notebooks/get_era5_metrics_netcdf.ipynb b/Jupyter_Notebooks/get_era5_metrics_netcdf.ipynb
new file mode 100644
index 00000000..9a07a3a2
--- /dev/null
+++ b/Jupyter_Notebooks/get_era5_metrics_netcdf.ipynb
@@ -0,0 +1,394 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5272c6d8-2d18-4de1-a437-5fe6461ca743",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os, sys\n",
+ "sys.path.append(\"../utils/\")\n",
+ "import xarray as xr\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "\n",
+ "from statistical_evaluation import perform_block_bootstrap_metric, avg_metrics, calculate_cond_quantiles, Scores"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "61e6b683-b9c6-4691-9c1f-130324ddb7a8",
+ "metadata": {},
+ "source": [
+ "# Evaluation of the ERA5 short-range forecasts\n",
+ "\n",
+ "Define the path to the file and load the data."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "453ad6f2-5655-47bf-8f39-7a1d73b059ab",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "indir = \"/p/home/jusers/langguth1/juwels/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/era5_forecast\"\n",
+ "yr=2019\n",
+ "\n",
+ "era5_fcst_file = os.path.join(indir, \"{0}_era5_short_range_fcst.nc\".format(yr))\n",
+ "\n",
+ "if not os.path.isfile(era5_fcst_file):\n",
+ " raise FileNotFoundError(\"Could not find file with all ERA5 forecasts '{0}'\".format(era5_fcst_file))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "459fa9bf-dbd0-48ee-8184-ccfec9ddbd59",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "era5_fcst = xr.open_dataset(era5_fcst_file)\n",
+ "print(era5_fcst)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "07b1324f-c6df-45c8-b77a-e334c2bb2966",
+ "metadata": {},
+ "source": [
+ "Next we initialize the function for calculating the MSE and call it to evaluate the ERA5 and persistence forecasts. <br>\n",
+ "If you require further evaluation metrics, just expand the cell accordingly, e.g. add the following lines <br>\n",
+ "```\n",
+ "ssim_func = Scores(\"ssim\", [\"lat\", \"lon\"]).score_func \n",
+ "\n",
+ "ssim_era5_all = ssim_func(data_fcst=era5_fcst[varname_fcst], data_ref=era5_fcst[varname_ref])\n",
+ "ssim_per_all = (data_fcst=era5_fcst[varname_per], data_ref=era5_fcst[varname_ref])\n",
+ "```\n",
+ "in case you want to evaluate the SSIM as well."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e066baeb-3190-4012-a0c1-ace1100be3aa",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fdata_clim = os.path.join(\"/p/project/deepacf/deeprain/video_prediction_shared_folder/preprocessedData/T2monthly/\", \"climatology_t2m_1991-2020.nc\")\n",
+ "\n",
+ "data = xr.open_dataset(fdata_clim)\n",
+ "\n",
+ "# copied from load_climdata in main_visualize_postprocess.py\n",
+ "var=\"var167\"\n",
+ "dt_clim = data[var]\n",
+ "\n",
+ "lon_dom = np.arange(0., 27.5, 0.3)\n",
+ "lat_dom = np.arange(54.9, 38.1, -0.3)\n",
+ "\n",
+ "# get the coordinates of the data after running CDO\n",
+ "coords = dt_clim.coords\n",
+ "nlat, nlon = len(coords[\"lat\"]), len(coords[\"lon\"])\n",
+ "# modify it our needs\n",
+ "coords_new = dict(coords)\n",
+ "coords_new.pop(\"time\")\n",
+ "coords_new[\"month\"] = np.arange(1, 13)\n",
+ "coords_new[\"hour\"] = np.arange(0, 24)\n",
+ "# initialize a new data array with explicit dimensions for month and hour\n",
+ "data_clim_new = xr.DataArray(\n",
+ " np.full((12, 24, nlat, nlon), np.nan),\n",
+ " coords=coords_new,\n",
+ " dims=[\"month\", \"hour\", \"lat\", \"lon\"],\n",
+ ")\n",
+ "# do the reorganization\n",
+ "for month in np.arange(1, 13):\n",
+ " data_clim_new.loc[dict(month=month)] = dt_clim.sel(\n",
+ " time=dt_clim[\"time.month\"] == month\n",
+ " )\n",
+ "\n",
+ "data_clim = data_clim_new.sel(lon=lons, lat=lats, tolerance=0.01, method=\"nearest\")\n",
+ "print(data_clim[\"lat\"])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0dd70ea1-9341-4b7a-9086-405ee75c6a64",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from tqdm import tqdm\n",
+ "from skimage.metrics import structural_similarity as ssim\n",
+ "\n",
+ "# overwrite some score-functions which inherently expect a separate fcst_hour-dimension next to batch_size \n",
+ "\n",
+ "def calc_acc_batch(data_fcst, data_ref, **kwargs):\n",
+ " \"\"\"\n",
+ " Calculate acc ealuation metric of forecast data w.r.t reference data\n",
+ " :param data_fcst: forecasted data (xarray with dimensions [batch, fore_hours, lat, lon])\n",
+ " :param data_ref: reference data (xarray with dimensions [batch, fore_hours, lat, lon])\n",
+ " :param data_clim: climatology data (xarray with dimensions [monthly, hourly, lat, lon])\n",
+ " :return: averaged acc for each batch example [batch, fore_hours]\n",
+ " \"\"\"\n",
+ " \n",
+ " print(\"Start calculating ACC\")\n",
+ " if \"data_clim\" in kwargs:\n",
+ " data_clim = kwargs[\"data_clim\"]\n",
+ " else:\n",
+ " raise KeyError(\"%{0}: climatological data must be parsed to calculate the ACC.\".format(method)) \n",
+ "\n",
+ " batch_size = data_fcst.shape[0]\n",
+ " acc = np.ones([batch_size])*np.nan\n",
+ " for i in tqdm(range(batch_size)):\n",
+ " img_fcst = data_fcst[i, ...]\n",
+ " img_ref = data_ref[i, ...]\n",
+ " # get the forecast time\n",
+ " img_month = img_fcst[\"fcst_hour\"].dt.month.values\n",
+ " img_hour = img_fcst[\"fcst_hour\"].dt.hour.values\n",
+ " img_clim = data_clim.sel(month=img_month, hour=img_hour) \n",
+ "\n",
+ " img1_ = img_ref - img_clim\n",
+ " img2_ = img_fcst - img_clim\n",
+ " cor1 = np.sum(img1_*img2_)\n",
+ " cor2 = np.sqrt(np.sum(img1_**2)*np.sum(img2_**2))\n",
+ " acc[i] = cor1/cor2\n",
+ " \n",
+ " # convert to data array \n",
+ " acc = xr.DataArray(acc, coords={\"fcst_hour\": data_fcst[\"fcst_hour\"]}, dims=[\"fcst_hour\"])\n",
+ " return acc\n",
+ "\n",
+ "def calc_ssim_batch(data_fcst, data_ref, **kwargs):\n",
+ " \"\"\"\n",
+ " Calculate ssim ealuation metric of forecast data w.r.t reference data\n",
+ " :param data_fcst: forecasted data (xarray with dimensions [batch, fore_hours, lat, lon])\n",
+ " :param data_ref: reference data (xarray with dimensions [batch, fore_hours, lat, lon])\n",
+ " :return: averaged ssim for each batch example, shape is [batch,fore_hours]\n",
+ " \"\"\"\n",
+ " method = Scores.calc_ssim_batch.__name__\n",
+ " batch_size = np.array(data_ref).shape[0]\n",
+ " ssim_pred = []\n",
+ " for i in tqdm(range(batch_size)):\n",
+ " ssim_pred.append(ssim(data_ref[i, ...],data_fcst[i,...]))\n",
+ " \n",
+ " # convert to data array \n",
+ " ssim_pred = xr.DataArray(np.asarray(ssim_pred), coords={\"fcst_hour\": data_fcst[\"fcst_hour\"]}, dims=[\"fcst_hour\"]) \n",
+ " return ssim_pred"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e76c4db9-f4da-4664-8054-fdd99cf5f64b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# to get and configure the score-functions\n",
+ "score_dims = [\"lat\", \"lon\"]\n",
+ "scores = [\"mse\", \"ssim\", \"acc\", \"texture\"]\n",
+ "# the reference data\n",
+ "varname_ref, varname_fcst, varname_per = \"2t_ref\", \"2t_era5_fcst\", \"2t_persistence_fcst\"\n",
+ "\n",
+ "# initialize empty dictionaries to store the score functions and to save the corresponding results\n",
+ "score_data = {}\n",
+ "for score in scores:\n",
+ " # overwrite functions that are incompatble here\n",
+ " if score == \"acc\":\n",
+ " score_func = calc_acc_batch\n",
+ " elif score == \"ssim\":\n",
+ " score_func = calc_ssim_batch\n",
+ " else:\n",
+ " score_func = Scores(score, score_dims).score_func\n",
+ " # parsing the persistence forecast as float32 increases throughput by a factor of 10!\n",
+ " score_data[score] = {\"era5_all\": score_func(data_fcst=era5_fcst[varname_fcst], data_ref=era5_fcst[varname_ref], data_clim=data_clim,),\n",
+ " \"per_all\": score_func(data_fcst=era5_fcst[varname_per].astype(\"float32\"), data_ref=era5_fcst[varname_ref], data_clim=data_clim)} "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c612a06f-e193-4d4a-85a4-1adc11fde81e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(score_data[\"mse\"][\"era5_all\"])\n",
+ "print(score_data[\"ssim\"][\"era5_all\"])\n",
+ "print(score_data[\"texture\"][\"era5_all\"])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b91695dc-7d3f-47de-8e67-03e5675aeac1",
+ "metadata": {},
+ "source": [
+ "Next, we initialize the data arrays to store the metrics for each forecast hour. <br>\n",
+ "Note that the ERA5 short-range forecasts only start twice a day at 06 and 18 UTC, respectively. Besides, the have only data starting from lead time 6 hours, but for consistency with the video prediction models, the data arrays cover all lead times between forecast hour 1 and 12. The unavailable values will be set to None."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "95611d9c-7551-466d-84b6-ca24be2d3977",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "init_times = [6, 18]\n",
+ "fcst_hours = np.arange(1, 13)\n",
+ "nhours = len(fcst_hours)\n",
+ "nboots=1000\n",
+ "\n",
+ "# MSE\n",
+ "mse_era5_fcst = xr.DataArray(np.empty(nhours, dtype=object), coords={\"fcst_hour\": fcst_hours}, dims=[\"fcst_hour\"])\n",
+ "mse_era5_fcst_boot = xr.DataArray(np.empty((nhours, nboots), dtype=object),\n",
+ " coords={\"fcst_hour\": fcst_hours, \"iboot\": np.arange(nboots)},\n",
+ " dims=[\"fcst_hour\", \"iboot\"])\n",
+ "\n",
+ "mse_per_fcst = xr.DataArray(np.empty(nhours, dtype=object), coords={\"fcst_hour\": fcst_hours}, dims=[\"fcst_hour\"])\n",
+ "mse_per_fcst_boot = xr.DataArray(np.empty((nhours, nboots), dtype=object),\n",
+ " coords={\"fcst_hour\": fcst_hours, \"iboot\": np.arange(nboots)},\n",
+ " dims=[\"fcst_hour\", \"iboot\"])\n",
+ "\n",
+ "# SSMI\n",
+ "ssim_era5_fcst, ssim_per_fcst = mse_era5_fcst.copy(), mse_per_fcst.copy()\n",
+ "ssim_era5_fcst_boot, ssim_per_fcst_boot = mse_per_fcst_boot.copy(), mse_per_fcst_boot.copy()\n",
+ "\n",
+ "# ACC\n",
+ "acc_era5_fcst, acc_per_fcst = mse_era5_fcst.copy(), mse_per_fcst.copy()\n",
+ "acc_era5_fcst_boot, acc_per_fcst_boot = mse_per_fcst_boot.copy(), mse_per_fcst_boot.copy()\n",
+ "\n",
+ "# ACC\n",
+ "txtr_era5_fcst, txtr_per_fcst = mse_era5_fcst.copy(), mse_per_fcst.copy()\n",
+ "txtr_era5_fcst_boot, txtr_per_fcst_boot = mse_per_fcst_boot.copy(), mse_per_fcst_boot.copy()\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e5700456-043b-4656-8ac4-759f42fc03c4",
+ "metadata": {},
+ "source": [
+ "Finally, we populate the initialized data arrays by looping over the forecast hours for which data is available. <br>\n",
+ "Additionally, we perform block bootstrapping to estimate the uncertainty of our evaluation metrics."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f3ea77c8-60c1-48bd-8285-95d6af66217a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def handle_scores(scores1_all, scores_ref_all, fhh):\n",
+ " \n",
+ " scores1 = scores1_all.sel(fcst_hour=(scores1_all.fcst_hour.dt.hour.isin(fhh)))\n",
+ " scores_ref = scores_ref_all.sel(fcst_hour=(scores_ref_all.fcst_hour.dt.hour.isin(fhh)))\n",
+ " score1_mean, score_ref_mean = scores1.mean(), scores_ref.mean()\n",
+ " # two runs per day -> 2*7 correpsonds to a block length of one week\n",
+ " score1_boot = perform_block_bootstrap_metric(scores1, \"fcst_hour\", 2*7)\n",
+ " score_ref_boot = perform_block_bootstrap_metric(scores_ref, \"fcst_hour\", 2*7)\n",
+ " \n",
+ " return score1_mean, score_ref_mean, score1_boot, score_ref_boot\n",
+ " \n",
+ "\n",
+ "for fh in fcst_hours[5::]:\n",
+ " print(\"Handling scores for forecast hour '{0:0d}'\".format(fh))\n",
+ " fh_curr = (init_times + fh)%24\n",
+ " # MSE\n",
+ " mse_era5_fcst[fh-1], mse_per_fcst[fh-1], mse_era5_fcst_boot[fh-1, :], mse_per_fcst_boot[fh-1, :] = handle_scores(score_data[\"mse\"][\"era5_all\"],\n",
+ " score_data[\"mse\"][\"per_all\"], fh_curr)\n",
+ " # SSIM\n",
+ " ssim_era5_fcst[fh-1], ssim_per_fcst[fh-1], ssim_era5_fcst_boot[fh-1, :], ssim_per_fcst_boot[fh-1, :] = handle_scores(score_data[\"ssim\"][\"era5_all\"],\n",
+ " score_data[\"ssim\"][\"per_all\"], fh_curr)\n",
+ " # ACC\n",
+ " acc_era5_fcst[fh-1], acc_per_fcst[fh-1], acc_era5_fcst_boot[fh-1, :], acc_per_fcst_boot[fh-1, :] = handle_scores(score_data[\"acc\"][\"era5_all\"],\n",
+ " score_data[\"acc\"][\"per_all\"], fh_curr) \n",
+ " # TEXTURE\n",
+ " txtr_era5_fcst[fh-1], txtr_per_fcst[fh-1], txtr_era5_fcst_boot[fh-1, :], txtr_per_fcst_boot[fh-1, :] = handle_scores(score_data[\"texture\"][\"era5_all\"],\n",
+ " score_data[\"texture\"][\"per_all\"], fh_curr)\n",
+ " \n",
+ " #mse_era5_curr = mse_era5_all.sel(fcst_hour=(mse_era5_all.fcst_hour.dt.hour.isin(fh_curr)))\n",
+ " #mse_per_curr = mse_per_all.sel(fcst_hour=(mse_per_all.fcst_hour.dt.hour.isin(fh_curr)))\n",
+ " #mse_era5_fcst[fh-1], mse_per_fcst[fh-1] = mse_era5_curr.mean(), mse_per_curr.mean()\n",
+ " ## two runs per day -> 2*7 correpsonds to a block length of one week\n",
+ " #mse_era5_fcst_boot[fh-1, :] = perform_block_bootstrap_metric(mse_era5_curr, \"fcst_hour\", 2*7)\n",
+ " #mse_per_fcst_boot[fh-1, :] = perform_block_bootstrap_metric(mse_per_curr, \"fcst_hour\", 2*7)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5ec002e6-cb1b-4c66-9ed0-9f5e2bd3fae9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(mse_era5_fcst)\n",
+ "print(ssim_era5_fcst)\n",
+ "print(acc_era5_fcst)\n",
+ "print(txtr_era5_fcst)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "276f4f4c-e926-4009-9198-c4e43271e87d",
+ "metadata": {},
+ "source": [
+ "Finally, we put the data arrays into a joint dataset and save the results into the netCDF-file."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "88ccbdce-5d6a-4c40-971c-ebe9a05a8c88",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# create Dataset and save to netCDF-file\n",
+ "ds_mse = xr.Dataset({\"2t_era5_mse_avg\": mse_era5_fcst, \"2t_era5_mse_bootstrapped\": mse_era5_fcst_boot, \n",
+ " \"2t_persistence_mse_avg\": mse_per_fcst, \"2t_persistence_mse_bootstrapped\": mse_per_fcst_boot,\n",
+ " # SSIM\n",
+ " \"2t_era5_ssim_avg\": ssim_era5_fcst, \"2t_era5_mse_bootstrapped\": ssim_era5_fcst_boot, \n",
+ " \"2t_persistence_ssim_avg\": ssim_per_fcst, \"2t_persistence_mse_bootstrapped\": ssim_per_fcst_boot,\n",
+ " # ACC\n",
+ " \"2t_era5_acc_avg\": acc_era5_fcst, \"2t_era5_acc_bootstrapped\": acc_era5_fcst_boot, \n",
+ " \"2t_persistence_acc_avg\": acc_per_fcst, \"2t_persistence_acc_bootstrapped\": acc_per_fcst_boot,\n",
+ " # TEXTURE \n",
+ " \"2t_era5_texture_avg\": txtr_era5_fcst, \"2t_era5_texture_bootstrapped\": txtr_era5_fcst_boot, \n",
+ " \"2t_persistence_texture_avg\": txtr_per_fcst, \"2t_persistence_texture_bootstrapped\": txtr_per_fcst_boot, \n",
+ " })\n",
+ "\n",
+ "outfile = os.path.join(indir, \"evaluation_metrics.nc\")\n",
+ "\n",
+ "print(\"Save evaluation metrics to '{0}'\".format(outfile))\n",
+ "ds_mse.to_netcdf(outfile)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a7678cc4-7333-402a-bcf4-1bc15712255d",
+ "metadata": {},
+ "source": [
+ "## DONE!"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "PyDeepLearning-1.1",
+ "language": "python",
+ "name": "pydeeplearning"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.8.5"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/Jupyter_Notebooks/get_metrics_joint_dom.ipynb b/Jupyter_Notebooks/get_metrics_joint_dom.ipynb
new file mode 100644
index 00000000..e54043e1
--- /dev/null
+++ b/Jupyter_Notebooks/get_metrics_joint_dom.ipynb
@@ -0,0 +1,338 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d636ee12-e299-485f-b84d-6f35c05fa766",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os, sys\n",
+ "import glob\n",
+ "sys.path.append(\"../utils/\")\n",
+ "import xarray as xr\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "\n",
+ "from statistical_evaluation import perform_block_bootstrap_metric, avg_metrics, calculate_cond_quantiles, Scores"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8510684d-9374-4e40-bc1c-4d69181c925c",
+ "metadata": {},
+ "source": [
+ "# Evaluation over a smaller (joint) domain\n",
+ "\n",
+ "The following cells will first merge all forecast files under `indir` into a single netCDF-file.<br>\n",
+ "Then the data is sliced to the domain defined by `lonlatbox` and all subsequent evaluation is performed on this smaller domain.<br>\n",
+ "The evaluation metrics are then saved to a file under `indir` named `evaluation_metrics_<nlon>x<nlat>.nc` where `nlat` and `nlon` denote the number of grid points/pixels in latitude and longitude direction of the smaller domain, respectively. <br>\n",
+ "\n",
+ "Thus, first let's define the basic parameters:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "440b15fa-ecd4-4bb4-9100-ede5abb2b04f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "indir = \"/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/savp/20210901T090059_gong1_savp_cv12/\"\n",
+ "model = \"savp\"\n",
+ "# define domain. [3., 24.3, 40.2, 53.1] corresponds to the smallest domain tested in the GMD paper\n",
+ "lonlatbox = [3., 24.3, 40.2, 53.1]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fd759f01-2561-4615-8056-036bdee6e2c7",
+ "metadata": {},
+ "source": [
+ "Next, we perform a first merging step. For computational efficiency, we merge max. 1000 files in the first step.<br>\n",
+ "Since the data is not sorted by the dimension `init_time` when querying along the sample index, we sort it before saving to intermediate files.<br>\n",
+ "\n",
+ "Given that the merging step has already been performed, no further processing is required.<br>\n",
+ "If this is not the case, we start with the sample indices between 0 and 999:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e6726da3-d774-4eda-89d6-e315a865bb99",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def get_fname_ndigits(indir, prefix, suffix, n, patt=\"[0-9]\"):\n",
+ " flist = []\n",
+ " for i in range(1, n+1):\n",
+ " fn_search = os.path.join(indir, \"{0}{1}{2}\".format(prefix, i*patt, suffix))\n",
+ " flist = flist + glob.glob(fn_search)\n",
+ " \n",
+ " if len(flist) == 0:\n",
+ " raise FileNotFoundError(\"Could not find any file under '{0}' with prefix '{1}' and suffix '{2}' containing digits.\".format(indir, prefix, suffix))\n",
+ " return flist\n",
+ "\n",
+ "# get list of files with sample index between 0 and 999.\n",
+ "vfp_list = get_fname_ndigits(indir, \"vfp_date_*sample_ind_\", \".nc\", 3)\n",
+ "outfile = os.path.join(indir, \"vfp_{0}_forecasts_sample_ind_0_999.nc\".format(model))\n",
+ "\n",
+ "if not os.path.isfile(outfile):\n",
+ " print(\"File '{0}' does not exist. \\n Start reading data with sample index between 0 and 999 from '{1}'...\".format(outfile, indir))\n",
+ " data_all = xr.open_mfdataset(vfp_list, concat_dim=\"init_time\", combine=\"nested\", decode_cf=True).load()\n",
+ " data_all = data_all.sortby(\"init_time\")\n",
+ " print(\"Data loaded successfully. Save merged data to '{0}'.\".format(outfile))\n",
+ " data_all.to_netcdf(outfile, encoding={'init_time':{'units': \"seconds since 1900-01-01 00:00:00\"}})"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1c0222c0-386d-44f4-9532-4e824b14828c",
+ "metadata": {},
+ "source": [
+ "Then, we proceed with the rest. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "54f4aa3e-3a39-496e-ae97-65f79d9cd598",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i in np.arange(1, 9):\n",
+ " outfile = os.path.join(indir, \"vfp_{0}_forecasts_sample_ind_{1:d}000_{1:d}999.nc\".format(model, i))\n",
+ " if not os.path.isfile(outfile):\n",
+ " print(\"File '{0}' does not exist. Start reading data with sample index between {1:d}000 and {1:d}999 from '{2}'...\".format(outfile, i, indir))\n",
+ " data_all = xr.open_mfdataset(os.path.join(indir, \"vfp_date_*sample_ind_{0}???.nc\".format(i)), concat_dim=\"init_time\", combine=\"nested\", decode_cf=True).load()\n",
+ " data_all = data_all.sortby(\"init_time\")\n",
+ " print(\"Data loaded successfully. Save merged data to '{0}'.\".format(outfile))\n",
+ " data_all.to_netcdf(outfile, encoding={'init_time':{'units': \"seconds since 1900-01-01 00:00:00\"}})"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "bdf16158-0ce5-40a3-848d-f574a1b9d622",
+ "metadata": {},
+ "source": [
+ "Still, xarray's `open_mfdataset`-method would not be able to concatenate all data since the `init_time`-dimension is not montonically increasing/decreasing when looping through the files. <br>\n",
+ "Thus, we have to merge the data manually.\n",
+ "The merged dataset is then saved to separate datafile for later computation.\n",
+ "\n",
+ "If the data has already been merged, we simply read the data from the corresponding netCDF-file."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "92f15edf-c23f-4803-b3c5-618305194de5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "outfile_all = os.path.join(indir, \"vfp_{0}_forecasts_all.nc\".format(model))\n",
+ "\n",
+ "if not os.path.isfile(outfile_all):\n",
+ " \n",
+ " print(\"netCDF-file with all forecasts '{0}' does not exist yet. Start merging and sorting all precursor files.\".format(outfile))\n",
+ " all_files = sorted(glob.glob(os.path.join(indir, \"vfp_{0}_forecasts_sample_ind_*.nc\".format(model))))\n",
+ " \n",
+ " if len(all_files) == 0:\n",
+ " raise FileNotFoundError(\"Could not find any precursor files.\")\n",
+ "\n",
+ " for i, f in enumerate(all_files):\n",
+ " print(\"Processing file '{0}'\".format(f))\n",
+ " tmp = xr.open_dataset(os.path.join(indir, f)).load()\n",
+ " if i == 0:\n",
+ " all_fcst = tmp.copy()\n",
+ " else:\n",
+ " print(\"Start merging\")\n",
+ " all_fcst = xr.merge([all_fcst, tmp]) \n",
+ "\n",
+ " # sort by init_time-dimension...\n",
+ " all_fcst = all_fcst.sortby(\"init_time\")\n",
+ " # ... and save to file\n",
+ " print(\"Finally, write all merged and sorted data to '{0}'.\".format(outfile_all))\n",
+ " all_fcst.to_netcdf(outfile_all)\n",
+ "else:\n",
+ " all_fcst = xr.open_dataset(outfile_all).load()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0fcf1cb1-ba0d-4262-8e23-12ba44b6e2d0",
+ "metadata": {},
+ "source": [
+ "Now, we slice the dataset to the domain of interest (defined by `lonlatbox`)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ede23e56-5be8-48be-b584-0eb8741acbf3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "all_fcst_sl = all_fcst.sel({\"lon\": slice(lonlatbox[0], lonlatbox[1]), \"lat\": slice(lonlatbox[3], lonlatbox[2])}) \n",
+ "print(all_fcst_sl)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e21b89c8-57ab-4070-9b4c-ec0fe24c37b9",
+ "metadata": {},
+ "source": [
+ "Next we initialize the function for calculating the MSE and call it to evaluate the ERA5 and persistence forecasts. <br>\n",
+ "If you require further evaluation metrics, just expand the cell accordingly, e.g. add the following lines <br>\n",
+ "```\n",
+ "ssim_func = Scores(\"ssim\", [\"lat\", \"lon\"]).score_func \n",
+ "\n",
+ "ssim_era5_all = ssim_func(data_fcst=era5_fcst[varname_fcst], data_ref=era5_fcst[varname_ref])\n",
+ "ssim_per_all = (data_fcst=era5_fcst[varname_per], data_ref=era5_fcst[varname_ref])\n",
+ "```\n",
+ "in case you want to evaluate the SSIM as well."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c2b70b80-6b86-4674-b051-6a23aaa821ea",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mse_func = Scores(\"mse\", [\"lat\", \"lon\"]).score_func\n",
+ "varname_ref, varname_fcst, varname_per = \"2t_ref\", \"2t_{0}_fcst\".format(model), \"2t_persistence_fcst\"\n",
+ "\n",
+ "mse_model_all = mse_func(data_fcst=all_fcst_sl[varname_fcst], data_ref=all_fcst_sl[varname_ref])\n",
+ "mse_per_all = mse_func(data_fcst=all_fcst_sl[varname_per], data_ref=all_fcst_sl[varname_ref])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7745356d-ad44-47b6-9655-8d6db3433b1a",
+ "metadata": {},
+ "source": [
+ "Then, we initialize the data arrays to store the desired evaluation metrics..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b49db031-126c-44b1-b649-4f70587fac89",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fcst_hours = all_fcst_sl[\"fcst_hour\"]\n",
+ "nhours = len(fcst_hours)\n",
+ "nboots=1000\n",
+ "\n",
+ "mse_model_fcst = xr.DataArray(np.empty(nhours, dtype=object), coords={\"fcst_hour\": fcst_hours}, dims=[\"fcst_hour\"])\n",
+ "mse_model_fcst_boot = xr.DataArray(np.empty((nhours, nboots), dtype=object),\n",
+ " coords={\"fcst_hour\": fcst_hours, \"iboot\": np.arange(nboots)},\n",
+ " dims=[\"fcst_hour\", \"iboot\"])\n",
+ "mse_per_fcst = xr.DataArray(np.empty(nhours, dtype=object), coords={\"fcst_hour\": fcst_hours}, dims=[\"fcst_hour\"])\n",
+ "mse_per_fcst_boot = xr.DataArray(np.empty((nhours, nboots), dtype=object),\n",
+ " coords={\"fcst_hour\": fcst_hours, \"iboot\": np.arange(nboots)},\n",
+ " dims=[\"fcst_hour\", \"iboot\"])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "55967405-02d1-46e8-b3c3-8952d0e28bd2",
+ "metadata": {},
+ "source": [
+ "... and populate them by looping over all forecast hours."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5090e71c-f20f-43e6-94f6-71cbd0b6006d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i, fh in enumerate(fcst_hours):\n",
+ " mse_model_curr = mse_model_all.sel(fcst_hour=fh)\n",
+ " mse_per_curr = mse_per_all.sel(fcst_hour=fh)\n",
+ " mse_model_fcst[fh-1], mse_per_fcst[fh-1] = mse_model_curr.mean(), mse_per_curr.mean()\n",
+ "\n",
+ " mse_model_fcst_boot[i, :] = perform_block_bootstrap_metric(mse_model_curr, \"init_time\", 24*7)\n",
+ " mse_per_fcst_boot[i, :] = perform_block_bootstrap_metric(mse_per_curr, \"init_time\", 24*7)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d526324d-5d19-4193-8208-e609d9c65205",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(mse_model_fcst)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b42cd738-b966-4b24-ad13-351d9b88f9e8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(mse_model_fcst)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b13c7287-7a8c-4133-bccc-f250bf25dad7",
+ "metadata": {},
+ "source": [
+ "Finally, we put the data arrays into a joint dataset and save the results into the netCDF-file."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7cd455ee-4749-46dd-8095-9e43744a1563",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# create Dataset and save to netCDF-file\n",
+ "ds_mse = xr.Dataset({\"2t_{0}_mse_avg\".format(model): mse_model_fcst, \"2t_{0}_mse_bootstrapped\".format(model): mse_model_fcst_boot, \n",
+ " \"2t_persistence_mse_avg\": mse_per_fcst, \"2t_persistence_mse_bootstrapped\": mse_per_fcst_boot})\n",
+ "\n",
+ "outfile = os.path.join(indir, \"evaluation_metrics_{0:d}x{1:d}.nc\".format(len(all_fcst_sl[\"lon\"]), len(all_fcst_sl[\"lat\"])))\n",
+ "\n",
+ "print(\"Save evaluation metrics to '{0}'\".format(outfile))\n",
+ "print(ds_mse)\n",
+ "ds_mse.to_netcdf(outfile)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "69bb9464-6fdb-489b-ba59-0170040144ee",
+ "metadata": {},
+ "source": [
+ "## Done!"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "PyDeepLearning-1.1",
+ "language": "python",
+ "name": "pydeeplearning"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.8.5"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/Jupyter_Notebooks/juwels_juwelsbooster_compare_old.ipynb b/Jupyter_Notebooks/juwels_juwelsbooster_compare_old.ipynb
deleted file mode 100644
index d788742d..00000000
--- a/Jupyter_Notebooks/juwels_juwelsbooster_compare_old.ipynb
+++ /dev/null
@@ -1,684 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": 16,
- "metadata": {},
- "outputs": [],
- "source": [
- "import os, glob\n",
- "import math\n",
- "import pickle\n",
- "import numpy as np\n",
- "import xarray as xr\n",
- "import matplotlib\n",
- "matplotlib.use('Agg')\n",
- "from matplotlib.transforms import Affine2D\n",
- "from matplotlib.patches import Polygon\n",
- "import matplotlib.pyplot as plt\n",
- "%matplotlib inline\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 2,
- "metadata": {},
- "outputs": [],
- "source": [
- "base = \"/p/project/deepacf/deeprain/video_prediction_shared_folder/models/\"+ \\\n",
- " \"era5-Y2010toY2222M01to12-160x128-2970N1500W-T2_MSL_gph500/convLSTM/\"\n",
- "fname_timing_train = \"/timing_training_time.pkl\"\n",
- "fname_timing_total = \"/timing_total_time.pkl\"\n",
- "\n",
- "fname_timing_iter = \"timing_per_iteration_time.pkl\""
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 3,
- "metadata": {},
- "outputs": [],
- "source": [
- "# some auxiliary functions\n",
- "def orderOfMagnitude(number):\n",
- " return np.floor(np.log(number, 10))\n",
- "\n",
- "def total_times(infile):\n",
- " with open(infile,'rb') as tfile:\n",
- " #print(\"Opening pickle time: '{0}'\".format(infile))\n",
- " total_time_sec = pickle.load(tfile)\n",
- " return np.asarray(total_time_sec/60)\n",
- "\n",
- "def log_total_times(infile):\n",
- " total_time_min = total_times(infile)\n",
- " return np.log(total_time_min)\n",
- "\n",
- "\n",
- "def get_time_dict(base, wildcardspec, tfilename, gpu_id_str=\"gpu\", llog = False):\n",
- " time_dict = {}\n",
- " flist_hpc = sorted(glob.glob(base + wildcardspec))\n",
- " wrapper = total_times\n",
- " if llog: wrapper = log_total_times\n",
- " for tfile in flist_hpc: \n",
- " ngpus = get_ngpus(tfile, gpu_id_str)\n",
- " time_dict[\"{0:d} GPU(s)\".format(ngpus)] = wrapper(tfile + tfilename)\n",
- " return time_dict\n",
- "\n",
- "def get_ngpus(fname, search_str, max_order=3):\n",
- " \"\"\"\n",
- " Tries to get numbers in the vicinty of search_str which is supposed to be a substring in fname.\n",
- " First seaches for numbers right before the occurence of search_str, then afterwards.\n",
- " :param fname: file name from which number should be inferred\n",
- " :param search_str: seach string for which number identification is considered to be possible\n",
- " :param max_order: maximum order of retrieved number (default: 3 -> maximum number is 999 then)\n",
- " :return num_int: integer of number in the vicintity of search string. \n",
- " \"\"\"\n",
- " \n",
- " ind_gpu_info = fname.lower().find(search_str)\n",
- " if ind_gpu_info == -1:\n",
- " raise ValueError(\"Unable to find search string '{0}' in file name '{1}'\".format(search_str, fname))\n",
- " \n",
- " # init loops\n",
- " fname_len = len(fname)\n",
- " success, flag = False, True\n",
- " indm = 1\n",
- " ind_sm, ind_sp = 0, 0\n",
- "\n",
- " # check occurence of numbers in front of search string\n",
- " while indm < max_order and flag:\n",
- " if ind_gpu_info - indm > 0:\n",
- " if fname[ind_gpu_info - indm].isnumeric():\n",
- " ind_sm += 1\n",
- " success = True\n",
- " else:\n",
- " flag = False\n",
- " else:\n",
- " flag = False\n",
- " indm += 1\n",
- " \n",
- "\n",
- " if not success: # check occurence of numbers after search string\n",
- " ind_gpu_info = ind_gpu_info + len(search_str)\n",
- " flag = True\n",
- " indm = 0\n",
- " while indm < max_order and flag: \n",
- " if ind_gpu_info + indm < fname_len:\n",
- " if fname[ind_gpu_info + indm].isnumeric():\n",
- " ind_sp += 1\n",
- " success = True\n",
- " else:\n",
- " flag = False\n",
- " else:\n",
- " flag = False\n",
- " indm += 1\n",
- " \n",
- " if success:\n",
- " return(int(fname[ind_gpu_info:ind_gpu_info+ind_sp]))\n",
- " else:\n",
- " raise ValueError(\"Search string found in fname, but unable to infer number of GPUs.\")\n",
- "\n",
- " else:\n",
- " return(int(fname[ind_gpu_info-ind_sm:ind_gpu_info]))\n",
- " \n",
- " \n",
- " "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 4,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Total computation with 16 GPU(s): 152.50984706878663\n",
- "Total computation with 32 GPU(s): 81.80640578667322\n",
- "Total computation with 4 GPU(s): 554.5182513117791\n",
- "Total computation with 64 GPU(s): 45.01537701288859\n",
- "Total computation with 8 GPU(s): 287.91878341039023\n"
- ]
- }
- ],
- "source": [
- "# Juwels\n",
- "wildcard_juwels = '20210115T135325_langguth1_test_venv_juwels_container*old'\n",
- "total_time_min_juwels = get_time_dict(base, wildcard_juwels, fname_timing_total, \"gpus\")\n",
- "training_time_min_juwels = get_time_dict(base, wildcard_juwels, fname_timing_train, \"gpus\")\n",
- "for key in training_time_min_juwels.keys():\n",
- " print(\"Total computation with {0}: {1}\".format(key, training_time_min_juwels[key]))\n",
- "\n",
- "overhead_time_juwels = {}\n",
- "for key in training_time_min_juwels.keys() & total_time_min_juwels.keys():\n",
- " overhead_time_juwels[key] = total_time_min_juwels[key] - training_time_min_juwels[key]\n",
- " \n",
- "#print('Juwels total time in minutes', get_time_d)\n",
- "#print('Juwels total training time in minutes', training_time_min_juwels)\n",
- "#overhead_time_juwels = np.array(total_time_min_juwels) - np.array(training_time_min_juwels)\n",
- "#print('Juwels overhead time in minutes', overhead_time_juwels)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 5,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Total computation with 1 GPU(s): 566.7376739541689\n",
- "Total computation with 4 GPU(s): 159.4931242307027\n",
- "Total computation with 8 GPU(s): 92.15467914342881\n",
- "Total computation with 16 GPU(s): 46.11619712909063\n",
- "Total computation with 32 GPU(s): 33.09077355464299\n",
- "Total computation with 64 GPU(s): 23.24405464331309\n"
- ]
- }
- ],
- "source": [
- "# Juwels booster\n",
- "wildcard_booster = '2020*gong1_booster_gpu*'\n",
- "total_time_min_booster = get_time_dict(base, wildcard_booster, fname_timing_total)\n",
- "training_time_min_booster = get_time_dict(base, wildcard_booster, fname_timing_train)\n",
- "for key in training_time_min_booster.keys():\n",
- " print(\"Total computation with {0}: {1}\".format(key, training_time_min_booster[key]))\n",
- "\n",
- "#print('Juwels Booster total time in minutes', list_times(base, wildcard_booster, filename_timing_total))\n",
- "#print('Juwels Booster total training time in minutes', list_times(base, wildcard_booster, filename_timing_train))\n",
- "overhead_time_booster = {}\n",
- "for key in training_time_min_booster.keys() & total_time_min_booster.keys():\n",
- " overhead_time_booster[key] = total_time_min_booster[key] - training_time_min_booster[key]\n",
- "#print('Juwels overhead time in minutes', overhead_time_booster)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 6,
- "metadata": {},
- "outputs": [],
- "source": [
- "def time_per_iteration_mean_std(infile):\n",
- " with open(infile, 'rb') as tfile:\n",
- " time_per_iteration_list = pickle.load(tfile) \n",
- " \n",
- " time_per_iteration = np.array(time_per_iteration_list)\n",
- " return np.mean(time_per_iteration), np.std(time_per_iteration)\n",
- "\n",
- "def iter_stat(base, wildcardspec, gpu_id_str=\"gpu\"):\n",
- " stat_iter_dict = {}\n",
- " flist_hpc = sorted(glob.glob(base + wildcardspec))\n",
- " for tdir in flist_hpc: \n",
- " ngpus = get_ngpus(tdir, gpu_id_str)\n",
- " ftname = os.path.join(tdir, fname_timing_iter)\n",
- " mean_loc, std_loc = time_per_iteration_mean_std(ftname)\n",
- " stat_iter_dict[\"{0:d} GPU(s)\".format(ngpus)] = {\"mean\": mean_loc , \"std\": std_loc}\n",
- " return stat_iter_dict\n",
- "\n",
- "def time_per_iteration_all(infile):\n",
- " with open(infile,'rb') as tfile:\n",
- " time_per_iteration_list = pickle.load(tfile)\n",
- " return np.asarray(time_per_iteration_list)\n",
- "\n",
- "def all_iter(base, wildcardspec, gpu_id_str=\"gpu\"):\n",
- " iter_dict = {}\n",
- " flist_hpc = sorted(glob.glob(base + wildcardspec))\n",
- " for tdir in flist_hpc: \n",
- " ngpus = get_ngpus(tdir, gpu_id_str)\n",
- " ftname = os.path.join(tdir, fname_timing_iter)\n",
- " iter_dict[\"{0:d} GPU(s)\".format(ngpus)] = time_per_iteration_all(ftname)\n",
- " return iter_dict \n",
- " "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 30,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "JUWELS (0.6151515198034729, 0.20104178037750603)\n",
- "Booster (0.3521572324468615, 0.3656996619706779)\n"
- ]
- }
- ],
- "source": [
- "# Juwels\n",
- "print('JUWELS', time_per_iteration_mean_std('/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2010toY2222M01to12-160x128-2970N1500W-T2_MSL_gph500/convLSTM/20201210T140958_stadtler1_comparison_1node_1gpu/timing_per_iteration_time.pkl'))\n",
- "# Booster\n",
- "print('Booster', time_per_iteration_mean_std('/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2010toY2222M01to12-160x128-2970N1500W-T2_MSL_gph500/convLSTM/20201210T141910_gong1_booster_gpu1/timing_per_iteration_time.pkl'))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 31,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Juwels mean and standart deviation {'16 GPU(s)': {'mean': 0.8209993402058342, 'std': 0.2627643291319852}, '32 GPU(s)': {'mean': 0.8590118098249986, 'std': 0.4078450977768068}, '4 GPU(s)': {'mean': 0.7445914211655112, 'std': 0.13789611351045}, '64 GPU(s)': {'mean': 0.9353915504630987, 'std': 0.6640973670265782}, '8 GPU(s)': {'mean': 0.7804724221628322, 'std': 0.21824334555299446}}\n"
- ]
- }
- ],
- "source": [
- "# Juwels\n",
- "print('Juwels mean and standart deviation',iter_stat(base, wildcard_juwels))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 9,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Booster mean and standart deviation {'1 GPU(s)': {'mean': 0.3521572324468615, 'std': 0.3656996619706779}, '4 GPU(s)': {'mean': 0.41844419631014446, 'std': 0.5273198599590724}, '8 GPU(s)': {'mean': 0.48867375665101026, 'std': 0.4378652997442439}, '16 GPU(s)': {'mean': 0.4786909431320202, 'std': 0.49638173862734053}, '32 GPU(s)': {'mean': 0.6439339113469129, 'std': 1.4395666886291258}, '64 GPU(s)': {'mean': 0.8176603168024377, 'std': 2.1044189535471185}}\n"
- ]
- }
- ],
- "source": [
- "# Booster\n",
- "print('Booster mean and standart deviation',iter_stat(base, wildcard_booster))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 34,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Plotting \n",
- "# Bar plot of total time and training time --> overhead time\n",
- "\n",
- "# dictionaries with the total times\n",
- "tot_time_juwels_dict = get_time_dict(base, wildcard_juwels, fname_timing_total)\n",
- "tot_time_booster_dict= get_time_dict(base, wildcard_booster, fname_timing_total)\n",
- "\n",
- "# dictionaries with the training times\n",
- "train_time_juwels_dict = get_time_dict(base, wildcard_juwels, fname_timing_train)\n",
- "train_time_booster_dict = get_time_dict(base, wildcard_booster, fname_timing_train)\n",
- "\n",
- "# get sorted arrays\n",
- "# Note: The times for Juwels are divided by 2, since the experiments have been performed with an epoch number of 20\n",
- "# instead of 10 (as Bing and Scarlet did)\n",
- "ngpus_sort = sorted([int(ngpu.split()[0]) for ngpu in tot_time_juwels_dict.keys()])\n",
- "nexps = len(ngpus_sort)\n",
- "tot_time_juwels = np.array([tot_time_juwels_dict[\"{0:d} GPU(s)\".format(key)] for key in ngpus_sort])/2.\n",
- "tot_time_booster = np.array([tot_time_booster_dict[\"{0:d} GPU(s)\".format(key)] for key in ngpus_sort])\n",
- "\n",
- "train_time_juwels = np.array([train_time_juwels_dict[\"{0:d} GPU(s)\".format(key)] for key in ngpus_sort])/2.\n",
- "train_time_booster = np.array([train_time_booster_dict[\"{0:d} GPU(s)\".format(key)] for key in ngpus_sort])\n",
- "\n",
- "overhead_juwels = tot_time_juwels - train_time_juwels \n",
- "overhead_booster= tot_time_booster - train_time_booster\n",
- "\n",
- "names = [\"Juwels\", \"Juwels Booster\"]"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 31,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "400.0\n",
- "278.0\n",
- "100.0\n",
- "2.0\n"
- ]
- }
- ],
- "source": [
- "plot_computation_times(tot_time_juwels, tot_time_booster, labels, [\"Juwels\", \"Juwels Booster\"], \\\n",
- " \"./total_computation_time\", log_yvals=False)\n",
- "\n",
- "plot_computation_times(overhead_juwels, overhead_booster, labels, [\"Juwels\", \"Juwels Booster\"], \\\n",
- " \"./overhead_time\")"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 12,
- "metadata": {},
- "outputs": [],
- "source": [
- "#print(labels)\n",
- "#raise ValueError(\"Stop!\")\n",
- "#x = np.arange(len(labels)) # the label locations\n",
- "#width = 0.35 # the width of the bars\n",
- "\n",
- "#fig, ax = plt.subplots()\n",
- "#rects1 = ax.bar(x - width/2, np.round(tot_time_juwels, 2), width, label='Juwels')\n",
- "#rects2 = ax.bar(x + width/2, np.round(tot_time_booster, 2), width, label='Booster')\n",
- "\n",
- "def plot_computation_times(times1, times2, ngpus, names, plt_fname, log_yvals = False):\n",
- " \n",
- " nlabels = len(ngpus)\n",
- " x_pos = np.arange(nlabels)\n",
- " \n",
- " bar_width = 0.35\n",
- " ytitle = \"Time\"\n",
- " ymax = np.ceil(np.maximum(np.max(times1)/100. + 0.5, np.max(times2)/100. + 0.5))*100.\n",
- " print(ymax) \n",
- " if log_yvals: \n",
- " times1, times2 = np.log(times1), np.log(times2)\n",
- " ytitle = \"LOG(Time) [min]\"\n",
- " ymax = np.ceil(np.maximum(np.max(times1)+0.5, np.max(times2) + 0.5))\n",
- " \n",
- " # create plot object\n",
- " fig, ax = plt.subplots()\n",
- " # create data bars\n",
- " rects1 = ax.bar(x_pos - bar_width/2, np.round(times1, 2), bar_width, label=names[0])\n",
- " rects2 = ax.bar(x_pos + bar_width/2, np.round(times2, 2), bar_width, label=names[1])\n",
- " # customize plot appearance\n",
- " # Add some text for labels, title and custom x-axis tick labels, etc.\n",
- " ax.set_ylabel(ytitle)\n",
- " ax.set_title('Comparison {0} and {1} with convLSTM model'.format(*names))\n",
- " ax.set_xticks(x_pos)\n",
- " ax.set_xticklabels(labels)\n",
- " ax.set_xlabel('# GPUs')\n",
- " print(np.ceil(np.maximum(np.max(times1)+0.5, np.max(times2) + 0.5)))\n",
- " ax.set_ylim(0., ymax)\n",
- " ax.legend()\n",
- " \n",
- " # add labels\n",
- " autolabel(ax, rects1)\n",
- " autolabel(ax, rects2)\n",
- " plt.savefig(plt_fname+\".png\")\n",
- " plt.close()\n",
- " \n",
- "\n",
- "def autolabel(ax, rects):\n",
- " \"\"\"Attach a text label above each bar in *rects*, displaying its height.\"\"\"\n",
- " for rect in rects:\n",
- " height = rect.get_height()\n",
- " ax.annotate('{}'.format(height),\n",
- " xy=(rect.get_x() + rect.get_width() / 2, height),\n",
- " xytext=(0, 3), # 3 points vertical offset\n",
- " textcoords=\"offset points\",\n",
- " ha='center', va='bottom')\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 13,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Plot mean + std \n",
- "# Juwels\n",
- "dict_stat_juwels = iter_stat(base, wildcard_juwels, gpu_id_str=\"gpu\")\n",
- "#print(dict_stat_juwels)\n",
- "iter_mean_juwels = np.array([dict_stat_juwels[\"{0:d} GPU(s)\".format(key)][\"mean\"] for key in labels])\n",
- "iter_std_juwels = np.array([dict_stat_juwels[\"{0:d} GPU(s)\".format(key)][\"std\"] for key in labels])\n",
- "\n",
- "dict_stat_booster = iter_stat(base, wildcard_booster, gpu_id_str=\"gpu\")\n",
- "iter_mean_booster = np.array([dict_stat_booster[\"{0:d} GPU(s)\".format(key)][\"mean\"] for key in labels])\n",
- "iter_std_booster = np.array([dict_stat_booster[\"{0:d} GPU(s)\".format(key)][\"std\"] for key in labels])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 29,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "(21225,)\n"
- ]
- }
- ],
- "source": [
- "iter_time_juwels = all_iter(base, wildcard_juwels)\n",
- "iter_time_booster= all_iter(base, wildcard_booster)\n",
- "\n",
- "max_iter_juwels = np.shape(iter_time_booster[\"{0:d} GPU(s)\".format(labels[0])])[0]\n",
- "max_iter_booster = np.shape(iter_time_booster[\"{0:d} GPU(s)\".format(labels[0])])[0]\n",
- "\n",
- "arr_iter_juwels = np.full((nexps, max_iter_juwels), np.nan)\n",
- "arr_iter_booster= np.full((nexps, max_iter_booster), np.nan)\n",
- "\n",
- "\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 37,
- "metadata": {},
- "outputs": [],
- "source": [
- "# box plot instead of errorbar plot\n",
- "# Juwels\n",
- "#data_juwels = list_time_per_iteration_all_runs(base, wildcard_juwels)\n",
- "data_juwels = all_iter(base, wildcard_juwels, gpu_id_str=\"gpu\")\n",
- "# Booster\n",
- "#data_booster = list_time_per_iteration_all_runs(base, wildcard_booster)\n",
- "data_booster = all_iter(base, wildcard_booster, gpu_id_str=\"gpu\")\n",
- "def simple_boxplot(time_per_iteration_data, title):\n",
- " # Multiple box plots on one Axes\n",
- " fig, ax = plt.subplots()\n",
- " ax.set_title(title)\n",
- " ax.boxplot(time_per_iteration_data, showfliers=False) # Outliers for initialization are disturbing \n",
- " plt.xticks([1, 2, 3, 4, 5 ,6], ['1', '4', '8', '16', '32', '64'])\n",
- " #plt.savefig('boxplot_'+title)\n",
- " #plt.close()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 86,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "886\n",
- "64.08639097213745\n",
- "31.232596397399902\n",
- "(1326,)\n",
- "***********\n",
- "2100\n",
- "4.405388832092285\n",
- "29.095214366912842\n",
- "(2653,)\n",
- "***********\n",
- "36981\n",
- "7.751298189163208\n",
- "26.409477949142456\n",
- "(42450,)\n",
- "***********\n",
- "3843\n",
- "66.00082683563232\n",
- "29.385547637939453\n",
- "(21225,)\n"
- ]
- }
- ],
- "source": [
- "print(np.argmax(data_booster[\"64 GPU(s)\"]))\n",
- "print(np.max(data_booster[\"64 GPU(s)\"]))\n",
- "print(data_booster[\"64 GPU(s)\"][0])\n",
- "print(np.shape(data_booster[\"64 GPU(s)\"]))\n",
- "print(\"***********\")\n",
- "\n",
- "print(np.argmax(data_juwels[\"64 GPU(s)\"][1::]))\n",
- "print(np.max(data_juwels[\"64 GPU(s)\"][1::]))\n",
- "print(data_juwels[\"64 GPU(s)\"][0])\n",
- "print(np.shape(data_juwels[\"64 GPU(s)\"]))\n",
- "print(\"***********\")\n",
- "\n",
- "print(np.argmax(data_juwels[\"4 GPU(s)\"][1::]))\n",
- "print(np.max(data_juwels[\"4 GPU(s)\"][1::]))\n",
- "print(data_juwels[\"4 GPU(s)\"][0])\n",
- "print(np.shape(data_juwels[\"4 GPU(s)\"]))\n",
- " \n",
- "print(\"***********\")\n",
- "print(np.argmax(data_booster[\"4 GPU(s)\"][1::]))\n",
- "print(np.max(data_booster[\"4 GPU(s)\"][1::]))\n",
- "print(data_booster[\"4 GPU(s)\"][0])\n",
- "print(np.shape(data_booster[\"4 GPU(s)\"]))\n",
- "\n",
- "#simple_boxplot(data_juwels, 'Juwels')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "simple_boxplot(data_booster, 'Booster')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 81,
- "metadata": {},
- "outputs": [],
- "source": [
- "# Try more fancy box plot \n",
- "def more_fancy_boxplot(time_per_iteration_data1, time_per_iteration_data2, ngpu_list, title):\n",
- " nexps = len(ngpu_list)\n",
- " # Shuffle data: EXPECT JUWELS FIRST FOR THE LEGEND! NOT GENERIC!\n",
- " data = []\n",
- " for i in np.arange(nexps):\n",
- " data.append(time_per_iteration_data1[\"{0} GPU(s)\".format(ngpu_list[i])])\n",
- " data.append(time_per_iteration_data2[\"{0} GPU(s)\".format(ngpu_list[i])])\n",
- " \n",
- " # trick to get list with duplicated entries\n",
- " xlabels = [val for val in ngpu_list for _ in (0, 1)]\n",
- "\n",
- " # Multiple box plots on one Axes\n",
- " #fig, ax = plt.subplots()\n",
- " fig = plt.figure(figsize=(6,4))\n",
- " ax = plt.axes([0.1, 0.15, 0.75, 0.75]) \n",
- " \n",
- " ax.set_title(title)\n",
- " bp = ax.boxplot(data, notch=0, sym='+', vert=1, whis=1.5, showfliers=False) # Outliers for initialization are disturbing\n",
- " plt.xticks(np.arange(1, nexps*2 +1), xlabels)\n",
- " ax.set_xlabel('# GPUs')\n",
- " ax.set_ylabel('Seconds')\n",
- " \n",
- " # Reference: https://matplotlib.org/3.1.1/gallery/statistics/boxplot_demo.html \n",
- " box_colors = ['darkkhaki', 'royalblue']\n",
- " num_boxes = len(data)\n",
- " medians = np.empty(num_boxes)\n",
- " for i in range(num_boxes):\n",
- " box = bp['boxes'][i]\n",
- " boxX = []\n",
- " boxY = []\n",
- " for j in range(5):\n",
- " boxX.append(box.get_xdata()[j])\n",
- " boxY.append(box.get_ydata()[j])\n",
- " box_coords = np.column_stack([boxX, boxY])\n",
- " # Alternate between Dark Khaki and Royal Blue\n",
- " ax.add_patch(Polygon(box_coords, facecolor=box_colors[i % 2]))\n",
- " # Now draw the median lines back over what we just filled in\n",
- " med = bp['medians'][i]\n",
- " medianX = []\n",
- " medianY = []\n",
- " for j in range(2):\n",
- " medianX.append(med.get_xdata()[j])\n",
- " medianY.append(med.get_ydata()[j])\n",
- " ax.plot(medianX, medianY, 'k')\n",
- " medians[i] = medianY[0]\n",
- " # Finally, overplot the sample averages, with horizontal alignment\n",
- " # in the center of each box\n",
- " ax.plot(np.average(med.get_xdata()), np.average(data[i]),\n",
- " color='w', marker='*', markeredgecolor='k')\n",
- " \n",
- " # Finally, add a basic legend\n",
- " fig.text(0.9, 0.15, 'Juwels',\n",
- " backgroundcolor=box_colors[0], color='black', weight='roman',\n",
- " size='small')\n",
- " fig.text(0.9, 0.09, 'Booster',\n",
- " backgroundcolor=box_colors[1],\n",
- " color='white', weight='roman', size='small')\n",
- " #fig.text(0.90, 0.015, '*', color='white', backgroundcolor='silver',\n",
- " # weight='roman', size='medium')\n",
- " fig.text(0.9, 0.03, '* Mean', color='white', backgroundcolor='silver',\n",
- " weight='roman', size='small')\n",
- "\n",
- " \n",
- " plt.savefig('fancy_boxplot_'+title.replace(' ', '_'))\n",
- " plt.close()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 82,
- "metadata": {},
- "outputs": [],
- "source": [
- "more_fancy_boxplot(data_juwels, data_booster, ngpus_sort, 'Time needed to iterate one step')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "flist_hpc1 = sorted(glob.glob(base + wildcard_juwels))\n",
- "flist_hpc2 = sorted(glob.glob(base + wildcard_booster))\n",
- "\n",
- "\n",
- " \n",
- "\n",
- "print(get_ngpus(flist_hpc1[2], \"gpu\"))\n",
- "print(get_ngpus(flist_hpc1[0], \"gpu\"))\n",
- "\n",
- "print(get_ngpus(flist_hpc2[2], \"gpu\"))\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.8.5"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
diff --git a/Jupyter_Notebooks/performance_check.ipynb b/Jupyter_Notebooks/performance_check.ipynb
deleted file mode 100644
index 3caf9018..00000000
--- a/Jupyter_Notebooks/performance_check.ipynb
+++ /dev/null
@@ -1,724 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": 108,
- "metadata": {},
- "outputs": [],
- "source": [
- "## import all required modules\n",
- "import os, glob\n",
- "import numpy as np\n",
- "import pickle\n",
- "# for plotting\n",
- "import matplotlib\n",
- "matplotlib.use('Agg')\n",
- "from matplotlib.transforms import Affine2D\n",
- "from matplotlib.patches import Polygon\n",
- "import matplotlib.pyplot as plt\n",
- "%matplotlib inline"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 144,
- "metadata": {},
- "outputs": [],
- "source": [
- "## some auxiliary functions\n",
- "#\n",
- "#colors = ['darkkhaki', 'royalblue']\n",
- "colors = [\"midnightblue\", \"darkorange\"]\n",
- "\n",
- "def val_order(number):\n",
- " return int(np.floor(np.log10(number)))\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "def get_ngpus(fname, search_str, max_order=3):\n",
- " \"\"\"\n",
- " Tries to get numbers in the vicinty of search_str which is supposed to be a substring in fname.\n",
- " First seaches for numbers right before the occurence of search_str, then afterwards.\n",
- " :param fname: file name from which number should be inferred\n",
- " :param search_str: seach string for which number identification is considered to be possible\n",
- " :param max_order: maximum order of retrieved number (default: 3 -> maximum number is 999 then)\n",
- " :return num_int: integer of number in the vicintity of search string. \n",
- " \"\"\"\n",
- " \n",
- " ind_gpu_info = fname.lower().find(search_str)\n",
- " if ind_gpu_info == -1:\n",
- " raise ValueError(\"Unable to find search string '{0}' in file name '{1}'\".format(search_str, fname))\n",
- " \n",
- " # init loops\n",
- " fname_len = len(fname)\n",
- " success, flag = False, True\n",
- " indm = 1\n",
- " ind_sm, ind_sp = 0, 0\n",
- " # check occurence of numbers in front of search string\n",
- " while indm < max_order and flag:\n",
- " if ind_gpu_info - indm > 0:\n",
- " if fname[ind_gpu_info - indm].isnumeric():\n",
- " ind_sm += 1\n",
- " success = True\n",
- " else:\n",
- " flag = False\n",
- " else:\n",
- " flag = False\n",
- " indm += 1\n",
- " # end while-loop\n",
- " if not success: # check occurence of numbers after search string\n",
- " ind_gpu_info = ind_gpu_info + len(search_str)\n",
- " flag = True\n",
- " indm = 0\n",
- " while indm < max_order and flag: \n",
- " if ind_gpu_info + indm < fname_len:\n",
- " if fname[ind_gpu_info + indm].isnumeric():\n",
- " ind_sp += 1\n",
- " success = True\n",
- " else:\n",
- " flag = False\n",
- " else:\n",
- " flag = False\n",
- " indm += 1\n",
- " # end while-loop \n",
- " if success:\n",
- " return(int(fname[ind_gpu_info:ind_gpu_info+ind_sp]))\n",
- " else:\n",
- " raise ValueError(\"Search string found in fname, but unable to infer number of GPUs.\")\n",
- "\n",
- " else:\n",
- " return(int(fname[ind_gpu_info-ind_sm:ind_gpu_info]))\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "# functions for computing time\n",
- "def compute_time_tot(infile):\n",
- " with open(infile,'rb') as tfile:\n",
- " #print(\"Opening pickle time: '{0}'\".format(infile))\n",
- " total_time_sec = pickle.load(tfile)\n",
- " return np.asarray(total_time_sec/60)\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "def compute_time_tot_log(infile):\n",
- " total_time_min = compute_time_tot(infile)\n",
- " return np.log(total_time_min)\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "def get_time_dict(base, wildcardspec, tfilename, gpu_id_str=\"gpu\", llog = False):\n",
- " time_dict = {}\n",
- " flist_hpc = sorted(glob.glob(base + wildcardspec))\n",
- " print(flist_hpc)\n",
- " wrapper = compute_time_tot\n",
- " if llog: wrapper = compute_time_tot_log\n",
- " for tfile in flist_hpc: \n",
- " ngpus = get_ngpus(tfile, gpu_id_str)\n",
- " time_dict[\"{0:d} GPU(s)\".format(ngpus)] = wrapper(tfile + tfilename)\n",
- " return time_dict\n",
- "#\n",
- "def calc_speedup(comp_time, ngpus, l_ideal= False):\n",
- " nn = np.shape(ngpus)[0]\n",
- " if l_ideal:\n",
- " spd_data = np.array(ngpus, dtype=float)\n",
- " else:\n",
- " spd_data = comp_time\n",
- "\n",
- " spd_up = spd_data[0:nn-1]/spd_data[1::]\n",
- " \n",
- " if l_ideal: spd_up = 1./spd_up\n",
- "\n",
- " return spd_up\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "# functions for iteration time data \n",
- "def iter_time_mean_std(infile):\n",
- " with open(infile, 'rb') as tfile:\n",
- " time_per_iteration_list = pickle.load(tfile) \n",
- " \n",
- " time_per_iteration = np.array(time_per_iteration_list)\n",
- " return np.mean(time_per_iteration), np.std(time_per_iteration)\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "def iter_stat(base, wildcardspec, gpu_id_str=\"gpu\"):\n",
- " stat_iter_dict = {}\n",
- " flist_hpc = sorted(glob.glob(base + wildcardspec))\n",
- " for tdir in flist_hpc: \n",
- " ngpus = get_ngpus(tdir, gpu_id_str)\n",
- " ftname = os.path.join(tdir, fname_timing_iter)\n",
- " mean_loc, std_loc = iter_time_mean_std(ftname)\n",
- " stat_iter_dict[\"{0:d} GPU(s)\".format(ngpus)] = {\"mean\": mean_loc , \"std\": std_loc}\n",
- " return stat_iter_dict\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "def read_iter_time(infile):\n",
- " with open(infile,'rb') as tfile:\n",
- " time_per_iteration_list = pickle.load(tfile)\n",
- " return np.asarray(time_per_iteration_list)\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "def get_iter_time_all(base, wildcardspec, gpu_id_str=\"gpu\"):\n",
- " iter_dict = {}\n",
- " flist_hpc = sorted(glob.glob(base + wildcardspec))\n",
- " for tdir in flist_hpc: \n",
- " ngpus = get_ngpus(tdir, gpu_id_str)\n",
- " ftname = os.path.join(tdir, fname_timing_iter)\n",
- " iter_dict[\"{0:d} GPU(s)\".format(ngpus)] = read_iter_time(ftname)\n",
- " return iter_dict \n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "# functions for plotting\n",
- "def autolabel(ax, rects, rot=45):\n",
- " \"\"\"Attach a text label above each bar in *rects*, displaying its height.\"\"\"\n",
- " scal = 1\n",
- " if rot <0.:\n",
- " scal = -1\n",
- " for rect in rects:\n",
- " height = rect.get_height()\n",
- " ax.annotate('{}'.format(height),\n",
- " xy=(rect.get_x() + rect.get_width()*scal, height),\n",
- " xytext=(0, 3), # 3 points vertical offset\n",
- " textcoords=\"offset points\",\n",
- " ha='center', va='bottom', rotation=rot)\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "def plot_computation_time(times1, times2, ngpus, names, plt_fname, log_yvals = False):\n",
- " \n",
- " nlabels = len(ngpus)\n",
- " x_pos = np.arange(nlabels)\n",
- " \n",
- " bar_width = 0.35\n",
- " ytitle = \"Time [min]\"\n",
- " max_time = np.maximum(np.max(times1), np.max(times2))\n",
- " time_order = val_order(max_time)\n",
- " ymax = np.ceil(max_time/(10**time_order) + 0.5)*(10**time_order) + 10**time_order\n",
- " # np.ceil(np.maximum(np.max(times1)/100. + 0.5, np.max(times2)/100. + 0.5))*100.\n",
- " if log_yvals: \n",
- " times1, times2 = np.log(times1), np.log(times2)\n",
- " ytitle = \"LOG(Time) [min]\"\n",
- " ymax = np.ceil(np.maximum(np.max(times1)+0.5, np.max(times2) + 0.5))\n",
- " \n",
- " # create plot object\n",
- " fig, ax = plt.subplots()\n",
- " # create data bars\n",
- " rects1 = ax.bar(x_pos - bar_width/2, np.round(times1, 2), bar_width, label=names[0], color=colors[0])\n",
- " rects2 = ax.bar(x_pos + bar_width/2, np.round(times2, 2), bar_width, label=names[1], color=colors[1])\n",
- " # customize plot appearance\n",
- " # Add some text for labels, title and custom x-axis tick labels, etc.\n",
- " ax.set_ylabel(ytitle)\n",
- " ax.set_title('Comparison {0} and {1} with convLSTM model'.format(*names))\n",
- " ax.set_xticks(x_pos)\n",
- " ax.set_xticklabels(ngpus)\n",
- " ax.set_xlabel('# GPUs')\n",
- " ax.set_ylim(0., ymax)\n",
- " ax.legend()\n",
- " \n",
- " # add labels\n",
- " autolabel(ax, rects1)\n",
- " autolabel(ax, rects2)\n",
- " print(\"Saving plot in file: {0}.png ...\".format(plt_fname))\n",
- " plt.savefig(plt_fname+\".png\")\n",
- " plt.close()\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "def plot_speedup(comp_time_hpc1, comp_time_hpc2, ngpus, names):\n",
- " fig = plt.figure(figsize=(6,4))\n",
- " ax = plt.axes([0.1, 0.15, 0.75, 0.75]) \n",
- " \n",
- " spd_up1 = calc_speedup(comp_time_hpc1, ngpus)\n",
- " spd_up2 = calc_speedup(comp_time_hpc2, ngpus)\n",
- " spd_ideal= calc_speedup(comp_time_hpc2, ngpus, l_ideal=True)\n",
- " \n",
- " plt.plot(spd_up1/spd_ideal, label= names[0], c=colors[0], lw=1.5)\n",
- " plt.plot(spd_up2/spd_ideal, label= names[1], c=colors[1], lw=1.5)\n",
- " plt.plot(spd_ideal/spd_ideal, label= \"Ideal\", c=\"r\", lw=3.)\n",
- " \n",
- " xlabels = []\n",
- " for i in np.arange(len(ngpus)-1):\n",
- " xlabels.append(\"{0} -> {1}\".format(ngpus[i], ngpus[i+1]))\n",
- " plt.xticks(np.arange(0, len(ngpus)-1), xlabels)\n",
- " ax.set_xlim(-0.5, len(ngpus)-1.5)\n",
- " ax.set_ylim(0.5, 1.5)\n",
- " legend = ax.legend(loc='upper left')\n",
- " ax.set_xlabel('GPU usage')\n",
- " ax.set_ylabel('Ratio Speedup factor') \n",
- " \n",
- " plt_fname = \"speed_up_{0}_vs_{1}.png\".format(*names)\n",
- " print(\"Saving plot in file: {0}.png ...\".format(plt_fname))\n",
- " plt.savefig(\"speed_up_{0}_vs_{1}.png\".format(*names))\n",
- "#\n",
- "# ****************************************************************************************************\n",
- "#\n",
- "def boxplot_iter_time(time_per_iteration_data1, time_per_iteration_data2, ngpu_list, names):\n",
- " nexps = len(ngpu_list)\n",
- " # create data lists for boxplot-routine\n",
- " data = []\n",
- " for i in np.arange(nexps):\n",
- " data.append(time_per_iteration_data1[\"{0} GPU(s)\".format(ngpu_list[i])])\n",
- " data.append(time_per_iteration_data2[\"{0} GPU(s)\".format(ngpu_list[i])])\n",
- " \n",
- " # trick to get list with duplicated entries\n",
- " xlabels = [val for val in ngpu_list for _ in (0, 1)]\n",
- "\n",
- " # Multiple box plots on one Axes\n",
- " #fig, ax = plt.subplots()\n",
- " fig = plt.figure(figsize=(6,4))\n",
- " ax = plt.axes([0.1, 0.15, 0.75, 0.75]) \n",
- " \n",
- " ax.set_title(\"Time per iteration step\")\n",
- " bp = ax.boxplot(data, notch=0, sym='+', vert=1, whis=1.5, showfliers=False) # Outliers for initialization are disturbing\n",
- " plt.xticks(np.arange(1, nexps*2 +1), xlabels)\n",
- " ax.set_xlabel('# GPUs')\n",
- " ax.set_ylabel('Time [s]')\n",
- " \n",
- " # Reference: https://matplotlib.org/3.1.1/gallery/statistics/boxplot_demo.html \n",
- " box_colors = colors\n",
- " num_boxes = len(data)\n",
- " medians = np.empty(num_boxes)\n",
- " for i in range(num_boxes):\n",
- " box = bp['boxes'][i]\n",
- " boxX = []\n",
- " boxY = []\n",
- " for j in range(5):\n",
- " boxX.append(box.get_xdata()[j])\n",
- " boxY.append(box.get_ydata()[j])\n",
- " box_coords = np.column_stack([boxX, boxY])\n",
- " # Alternate between Dark Khaki and Royal Blue\n",
- " ax.add_patch(Polygon(box_coords, facecolor=box_colors[i % 2]))\n",
- " # Now draw the median lines back over what we just filled in\n",
- " med = bp['medians'][i]\n",
- " medianX = []\n",
- " medianY = []\n",
- " for j in range(2):\n",
- " medianX.append(med.get_xdata()[j])\n",
- " medianY.append(med.get_ydata()[j])\n",
- " ax.plot(medianX, medianY, 'k')\n",
- " medians[i] = medianY[0]\n",
- " # Finally, overplot the sample averages, with horizontal alignment\n",
- " # in the center of each box\n",
- " ax.plot(np.average(med.get_xdata()), np.average(data[i]),\n",
- " color='w', marker='*', markeredgecolor='k', markersize=10)\n",
- " \n",
- " # Finally, add a basic legend\n",
- " fig.text(0.86, 0.15, names[0],\n",
- " backgroundcolor=box_colors[0], color='white', weight='roman',\n",
- " size='small')\n",
- " fig.text(0.86, 0.09, names[1],\n",
- " backgroundcolor=box_colors[1],\n",
- " color='white', weight='roman', size='small')\n",
- " #fig.text(0.90, 0.015, '*', color='white', backgroundcolor='silver',\n",
- " # weight='roman', size='medium')\n",
- " #fig_transform = ax.figure.transFigure #+ ax.transAxes.inverted() #+ ax.figure.transFigure.inverted()\n",
- " #ax.plot(0.1, 0.03, marker='*', markersize=30, color=\"w\", markeredgecolor=\"k\", transform=fig_transform)\n",
- " fig.text(0.86, 0.03, '* Mean', color='black', backgroundcolor='white', \n",
- " weight='roman', size='small', bbox=dict(facecolor='none', edgecolor='k'))\n",
- "\n",
- " plt_fname = \"boxplot_iter_time_{0}_vs_{1}\".format(*names)\n",
- " print(\"Saving plot in file: {0}.png ...\".format(plt_fname))\n",
- " plt.savefig(plt_fname+\".png\")\n",
- " plt.close()\n",
- " \n",
- " "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 110,
- "metadata": {},
- "outputs": [],
- "source": [
- "## some basic settings\n",
- "base_dir = \"/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/\"\n",
- "\n",
- "wildcard_hpc1 = '20210325T095504_langguth1_juwels_container_[1-9]*gpu*' # search pattern for finding the experiments\n",
- "wildcard_hpc2 = '20210325T095504_langguth1_jwb_container_[1-9]*gpu*'\n",
- "\n",
- "gpu_id_str = [\"gpu\", \"gpu\"] # search substring to get the number of GPUs used in the experiments,\n",
- " # e.g. \"gpu\" if '64gpu' is a substring in the experiment directory\n",
- " # or \"ngpu\" if 'ngpu64' is a substring in the experiment directory\n",
- " # -> see wilcard-variables above\n",
- "names_hpc = [\"Juwels\", \"Booster\"]\n",
- "\n",
- "# name of pickle files tracking computing time\n",
- "fname_timing_train = \"/timing_training_time.pkl\"\n",
- "fname_timing_total = \"/timing_total_time.pkl\"\n",
- "\n",
- "fname_timing_iter = \"timing_per_iteration_time.pkl\"\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 111,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "['/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_16gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_1gpu', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_32gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_4gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_64gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_8gpus']\n",
- "['/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_16gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_1gpu', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_32gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_4gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_64gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_8gpus']\n",
- "{'16 GPU(s)': array(53.40843068), '1 GPU(s)': array(930.4968381), '32 GPU(s)': array(45.96871045), '4 GPU(s)': array(217.45655225), '64 GPU(s)': array(35.7369519), '8 GPU(s)': array(106.4218419)}\n",
- "{'16 GPU(s)': array(34.26928383), '1 GPU(s)': array(492.70926997), '32 GPU(s)': array(35.05492661), '4 GPU(s)': array(100.99109779), '64 GPU(s)': array(30.98471271), '8 GPU(s)': array(49.63896298)}\n",
- "['/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_16gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_1gpu', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_32gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_4gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_64gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_juwels_container_8gpus']\n",
- "['/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_16gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_1gpu', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_32gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_4gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_64gpus', '/p/project/deepacf/deeprain/video_prediction_shared_folder/models/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM_container/20210325T095504_langguth1_jwb_container_8gpus']\n"
- ]
- }
- ],
- "source": [
- "## evaluate computing time\n",
- "# dictionaries with the total times\n",
- "tot_time_hpc1_dict = get_time_dict(base_dir, wildcard_hpc1, fname_timing_total, gpu_id_str=gpu_id_str[0])\n",
- "tot_time_hpc2_dict= get_time_dict(base_dir, wildcard_hpc2, fname_timing_total, gpu_id_str=gpu_id_str[1])\n",
- "\n",
- "print(tot_time_hpc1_dict)\n",
- "print(tot_time_hpc2_dict)\n",
- "\n",
- "# dictionaries with the training times\n",
- "train_time_hpc1_dict = get_time_dict(base_dir, wildcard_hpc1, fname_timing_train, gpu_id_str=gpu_id_str[0])\n",
- "train_time_hpc2_dict = get_time_dict(base_dir, wildcard_hpc2, fname_timing_train, gpu_id_str=gpu_id_str[1])\n",
- "\n",
- "# get sorted arrays\n",
- "# Note: The times for Juwels are divided by 2, since the experiments have been performed with an epoch number of 20\n",
- "# instead of 10 (as Bing and Scarlet did)\n",
- "ngpus_sort = sorted([int(ngpu.split()[0]) for ngpu in tot_time_hpc1_dict.keys()])\n",
- "nexps = len(ngpus_sort)\n",
- "tot_time_hpc1 = np.array([tot_time_hpc1_dict[\"{0:d} GPU(s)\".format(key)] for key in ngpus_sort])\n",
- "tot_time_hpc1[0] = tot_time_hpc1[0]#*2.\n",
- "tot_time_hpc2 = np.array([tot_time_hpc2_dict[\"{0:d} GPU(s)\".format(key)] for key in ngpus_sort])\n",
- "\n",
- "train_time_hpc1 = np.array([train_time_hpc1_dict[\"{0:d} GPU(s)\".format(key)] for key in ngpus_sort])\n",
- "train_time_hpc1[0] = train_time_hpc1[0]#*2.\n",
- "train_time_hpc2 = np.array([train_time_hpc2_dict[\"{0:d} GPU(s)\".format(key)] for key in ngpus_sort])\n",
- "\n",
- "overhead_hpc1 = tot_time_hpc1 - train_time_hpc1\n",
- "overhead_hpc2= tot_time_hpc2 - train_time_hpc2"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 112,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "[492.70926997 100.99109779 49.63896298 34.26928383 35.05492661\n",
- " 30.98471271]\n",
- "Saving plot in file: ./total_computation_time_Juwels_vs_Booster.png ...\n",
- "Saving plot in file: ./overhead_time_Juwels_vs_Booster.png ...\n",
- "Saving plot in file: speed_up_Juwels_vs_Booster.png.png ...\n"
- ]
- },
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEFCAYAAAAYKqc0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA010lEQVR4nO3dd3hUZfbA8e9JCAQIzSQgEJqKSJFEWmyAgAIiihQRFYG1YENl7e4urmuvu66u5afrigoiiqKABZUiiIgQJNIFUSGClAhIqCnn98d7EyOkTMhMbjI5n+e5DzN37tx77iScefPe955XVBVjjDHhJcLvAIwxxgSfJXdjjAlDltyNMSYMWXI3xpgwZMndGGPCkCV3Y4wJQ1X8DqCk4uLitHnz5n6HYYwxvktJSdmhqvEFvVbhknvz5s1ZsmSJ32EYY4zvROSnwl6zbhljjAlDltyNMSYMWXI3xpgwVOH63AuSmZlJWloaBw4c8DuUcis6OpqEhASioqL8DsUYUwbCIrmnpaVRq1Ytmjdvjoj4HU65o6qkp6eTlpZGixYt/A7HGFMGwqJb5sCBA8TGxlpiL4SIEBsba3/ZGFOJhEVyByyxF8M+H2Mql7BJ7uVBTExMSPZ777338sQTT4Rk38aY8GTJ3RhjwpAl9yCbO3cu/fv3z3s+ZswYxo8fz+LFixk0aBAA77//PtWrV+fQoUMcOHCA4447DoDvv/+evn370rFjR7p27cqaNWuO2P/TTz9NmzZtaN++PcOGDSubkzLGVDhhMVomv3vumcrKlT8HdZ9t2zbmvvsGlmofp5xyCsuWLQNg/vz5tGvXjsWLF5OVlUVycjIAo0eP5oUXXqBly5YsWrSI66+/ntmzZ/9hP4888gg//PAD1apVY9euXaWKyRgTvsIuuZdXVapU4fjjj2f16tV8/fXX3HLLLcybN4/s7Gy6du1KRkYGX375JRdddFHeew4ePHjEftq3b89ll13GhRdeyIUXXliGZ2CMqUjCLrmXtoVdWlWqVCEnJyfvef7hh926deOjjz4iKiqKs88+m1GjRpGdnc3jjz9OTk4OdevWzWvdF+aDDz5g3rx5TJ8+nQcffJDly5dTpUrY/RiNMaVkfe5B1qxZM1atWsXBgwfZtWsXs2bNynuta9euPPXUU5x22mnEx8eTnp7O2rVradeuHbVr16ZFixa8/fbbgLvxKDU19Q/7zsnJYdOmTfTo0YNHH32U3bt3k5GRUabnZ4ypGKzJFyRZWVlUq1aNJk2aMHToUNq1a0eLFi045ZRT8rZJTk5m69atdOvWDXBdLL/88kveGPSJEydy3XXX8cADD5CZmcmwYcNITEzMe392djbDhw9n9+7dqCo33XQTdevWLdPzNMZUDKKqfsdQIp06ddLD67mvXr2a1q1b+xSRk5qaytVXX83XX3/taxxFKQ+fkzEmeEQkRVU7FfRayLplROR/IrJNRFYUs11nEckSkSGhiiXUXnjhBS655BIeeOABv0MxxhggtH3u44G+RW0gIpHAo8AnIYwj5K699lpWrVpF7969/Q7FGGOAECZ3VZ0H/FrMZjcC7wDbQhWHMcZURr6NlhGRxsBA4PkAth0tIktEZMn27dtDH5wxxlRwfg6FfAq4U1VzittQVV9U1U6q2ik+vsCJvo0xxuTj51DITsCb3jDAOKCfiGSp6ns+xmSMMWHBt5a7qrZQ1eaq2hyYAlxfkRN7ZGQkSUlJJCYm0qFDB7788sug7v+hhx4K6v6MMeEtlEMhJwELgVYikiYiV4rItSJybaiO6afq1auzbNkyUlNTefjhh7n77ruDuv+jSe7Z2dlBjcEYU3GEcrTMJaraUFWjVDVBVV9W1RdU9YUCth2lqlNCFUtZ++2336hXrx7gygjcfvvttGvXjpNPPpnJkycXuX7Lli1069aNpKQk2rVrx/z587nrrrvYv38/SUlJXHbZZQBMmDCBLl26kJSUxDXXXJOXyGNiYrj11ltJTExk4cKFPpy9MaY8CL/yA3PGwrZlwd1n/STo8VSRm+Qm3wMHDrBly5a8Ur3vvvtuXot+x44ddO7cmW7duvHll18WuP6NN96gT58+/PWvfyU7O5t9+/bRtWtX/vOf/+QVFVu9ejWTJ09mwYIFREVFcf311zNx4kRGjBjB3r17SU5O5sknnwzuZ2CMqVDCL7n7JLdbBmDhwoWMGDGCFStW8MUXX3DJJZcQGRlJgwYN6N69O4sXLy50fefOnbniiivIzMzkwgsvJCkp6YhjzZo1i5SUFDp37gy4L5b69esDru9/8ODBZXXaxphyKvySezEt7LJw2mmnsWPHDo5mTH63bt2YN28eH3zwAaNGjeKWW25hxIgRf9hGVRk5ciQPP/zwEe+Pjo4mMjLyqGM3xoQHK/kbAmvWrCE7O5vY2Fi6du3K5MmTyc7OZvv27cybN48uXboUuv6nn36iQYMGXH311Vx11VUsXboUgKioKDIzMwHo1asXU6ZMYds2d2Pvr7/+yk8//eTb+Rpjyp/wa7n7JLfPHVzL+tVXXyUyMpKBAweycOFCEhMTEREee+wxjj322ELXv/rqqzz++ONERUURExPDa6+9Brgp+Nq3b0+HDh2YOHEiDzzwAL179yYnJ4eoqCieffZZmjVr5uMnYIwpT6zkbyVin5Mx4cWXkr/GGGP8Y8ndGGPCkCV3Y4wJQ5bcjTEmDFlyN8aYMGTJ3RhjwpAl9yCJiYkpcP2oUaOYMuXoaqLde++9PPHEE6UJyxhTSVlyN8aYMGTJPchUlTFjxtCqVSvOPvvsvBIBACkpKXTv3p2OHTvSp08ftmzZAsBLL71E586dSUxMZPDgwezbt8+v8I0xYSL8krtI6JYATJ06lbVr17Jq1Spee+21vBmZMjMzufHGG5kyZQopKSlcccUV/PWvfwVg0KBBLF68mNTUVFq3bs3LL78cso/HGFM5WG2ZIJs3b15eKd9GjRrRs2dPANauXcuKFSs455xzADdLUsOGDQFYsWIFf/vb39i1axcZGRn06dPHt/iNMeHBknsZUVXatm1b4OxIo0aN4r333iMxMZHx48czd+7csg/QGBNWwq9bRjV0SwC6deuWV8p3y5YtzJkzB4BWrVqxffv2vOSemZnJypUrAdizZw8NGzYkMzOTiRMnhuZzMcZUKtZyD7KBAwcye/Zs2rRpQ9OmTTnttNMAqFq1KlOmTOGmm25i9+7dZGVlMXbsWNq2bcv9999PcnIy8fHxJCcns2fPHp/PwhhT0VnJ30rEPidjwouV/DXGmErGkrsxxoQhS+7GGBOGwia5V7RrB2XNPh9jKpewSO7R0dGkp6dbAiuEqpKenk50dLTfoRhjykhYDIVMSEggLS2N7du3+x1KuRUdHU1CQoLfYRhjykhYJPeoqChatGjhdxjGGFNuhEW3jDHGmD+y5G6MMWGoyOQuIpEiYsVOjDGmgikyuatqNtBMRKqWUTzGGGOCIJALqhuABSIyDdibu1JV/xmyqMorVVjyJJx8FUTX9TsaY4wpVCB97t8DM7xta+VbKp+tS2D+XfBqO/hxpt/RGGNMoYptuavqPwBEJMZ7nhHIjkXkf0B/YJuqtivg9cuAOwEB9gDXqWpq4KH74NjOcOlC+GgkvNMXTr4aznoSqlbO7zpjTPlVbMtdRNqJyDfASmCliKSISNsA9j0e6FvE6z8A3VX1ZOB+4MUA9um/YzvD5Uuh8x2w4mV49WTYONvvqIwx5g8C6ZZ5EbhFVZupajPgVuCl4t6kqvOAX4t4/UtV3ek9/QqoOLdPVomGbo/CsC8gsiq83QtmjYHMvcW/1xhjykAgyb2mqs7JfaKqc4GaQY7jSuCjIO8z9BqdBpcvgw5jYdlz8FoipH3hd1TGGBNQct8gIuNEpLm3/A03giYoRKQHLrnfWcQ2o0VkiYgsKXf1Y6JqQI9/wcVzQXNgcjeYeytk7vc7MmNMJRZIcr8CiAfeBd4B4oA/BePgItIe+C8wQFXTC9tOVV9U1U6q2ik+Pj4Yhw6+hG4w4ltIvBZS/gmvnwJbFvkdlTGmkgokuZ+tqjepagdV7aiqY4FzSntgEWmK+8K4XFW/K+3+yoWqMXD2czDkU8jaD5NOh/l3Q9ZBvyMzxlQygST3uwNc9wciMglYCLQSkTQRuVJErhWRa71N7gFigedEZJmILCl0ZxVNs7Nh5HJo+yf4+hGY0BG2pvgdlTGmEpHCJrgQkXOBfsBQYHK+l2oDbVS1S+jDO1KnTp10yZIK9D3ww0fwyVWwdysk/xVO/asbYWOMMaUkIimq2qmg14pquW8GlgAHgJR8yzSgT7CDDFstzoWRK6D1pfDVfTAxGbZ/63dUxpgwV2jLPW8DkdrAXq+IGCISCVRT1X1lEN8RKlzLPb/178Ono+HATjjt79DlTogIi/lSjDE+ONqWe65PgOr5nlcHPgtGYJXOCQNg5EpoOQgW/M1dcE1f7XdUxpgwFEhyj85fT8Z7XCN0IYW5GnHQ/03o/xbs2uCGTC5+HHKy/Y7MGBNGAknue0WkQ+4TEekI2B06pdXqIhi10vXJz7vD3fy0c53fURljwkQgyX0s8LaIzBeRL3AjZ8aENKrKomYDuOBd6DcB0le58gVL/+3udDXGmFIIpOTvYhE5CWjlrVqrqpmhDasSEYHWl0GTHvDJ1TBnLKybCn3+B3WP8zs6Y0wFFegE2a2ANkAH4BIRGRG6kCqpmEYwcIZL6tu+gdfaQ+oLbvYnY4wpoUDquf8deMZbegCPAReEOK7KSQTa/cnd3drodPjsOpjSG37b6HdkxpgKJpCW+xCgF/CLqv4JSATqhDSqyq52Uxg8E85+AbYsdBOCLP+fteKNMQELJLnvV9UcIMu7oWkb0CS0YRlEIPEa14qvfwp8ciVM7Q8Zm/2OzBhTAQSS3JeISF3c7EspwFJcQTBTFuq0gKGzoce/YdMcGN8WVk2wVrwxpkiFJncROcN7+GdV3aWqL+BK/Y70umdMWZEI6HATjEiF2Dbw0eUwbZArRmaMMQUoquX+tPdvXitdVX9UVat65Zd6LeHiedDtcVdtcnxbWPuW31EZY8qhosa5Z4rIi0CCiDx9+IuqelPowjKFioiEzrfBcf3g41Ew42L47h3o9awrbWCMMRTdcu8PzMaVGkgpYDF+im0Dl3wJZz4I66fCq21h3Xt+R2WMKScKbbmr6g7gTRFZraqpZRiTCVREFUj+CxzXHz4aCdMGQuvh0PNpiK7nd3TGGB8VO1rGEnsFEN8eLvva1Yhf+ya82g42fOh3VMYYHwVafsCUd5FRcPq9cOkiqFYPpp4HM6+Eg7v9jswY4wNL7uGmQQcYngJd7oaV493drT/Z3CrGVDaB1JaJFZFnRGSpiKSIyL9FJLYsgjNHqUo16PqQu+BapQZMOcfVqTmUUfx7jTFhIZCW+5u4kgODcXVmtuNqulc6OTk5ZGQc8DuMwDVMhsu/gY63QOr/uUqTmz73OypjTBkIZILsFara7rB1y1X15JBGVgg/J8jetOlXkpPvp1GjurRs2eCIJTY2xpe4ApL2BcwcBbu+hw43w5kPQZTNlmhMRVbUBNnFTtYBfCIiw4DcWyGHADODFVxFEh0dxV139WPduq2sW7eVN974in37DuW9Xq9eTVq2rH9E0m/UqC4RET5f3kg405UvmHeXm+3phw+hz3hofLq/cRljQiKQlvseoCaQO/dbBLDXe6yqWjt04R3Jz5b74XJycti8eTfr12/NS/i5y6+/7s3brnr1qpxwwpFJv3nzOKKiIss+8I2zYeYVsGcTdLwVzrgPqkSXfRzGmFIpquVebHIvb8pTci9KenrGEQl/3bqtbN68K2+bKlUiaN487oikf8IJ9alRo1poAzy0Bz6/Db59EY5pDee+Csd2Du0xjTFBVarkLiLdClqvqvOCEFuJVZTkXpi9ew/ma+lvy0v6P/64g+zs3yfGbty4npfs6+dL+iHo1/9xphsPv/cX6HIXnDrOjbYxxpR7pU3u0/M9jQa6ACmq2jN4IQauoif3whw6lMWPP+44oqW/fv02Dhz4fT7yY46pScuWDTjxxN8Tfm6/vogc3cEP7IK5f3bj4uNOhnNfg/pJwTgtY0wIBbVbRkSaAE+p6uBgBFdS4ZrcC5OTk8PPP+8qIOlvZefOfXnb1axZ7Q/9+rlJv3nzWKpUCbBf//vp8Olo2L/DteC73O3ufDXGlEvBTu4CrFTVNsEIrqQqW3IvjKoW0q+/jS1bduVtFxUVSYsWcUck/eOPr0+NGlWP3PH+dJh9E6x5Axp0hL7jIa7dkdsZY3xX2m6ZZ4DcjSKAJOBHVR0ezCADZcm9eBkZB1i/ftsRif+nn9Lz+vVFhISEeoclfdfyr1evpqsR/9l1cGg3nPYPV0M+IpCRs8aYslLa5D4y39MsXGJfEMT4SsSS+9E7eDCLH3/cntfC/+67X1i3bisbNmz/Q79+XFwMLVs2IOnEaIYfO4EWWZ9zKLYTUee/jsSe5OMZGGPys6GQpkjZ2Tmkpf16xAie9eu3snv3Pga0+o4He82hRlQWr2/ox/Kqg4mLr0tcXAyxsTHExcXkPY6NjaF69QK6e4wxQXdUyV1ElvN7d8wRVLV9cMIrGUvuZUdV2b59D+vWbSVt7Urab32Q1lWXsHRbM0ZP78/mXQVfbK1ZsxqxsTWJi6uVl/AP/wJwr9UkNjaGqlWtu8eYo3G0yb2Z9/AG79/XvX+H4+5MvauYg/4PN1XftsNr03ivC/BvoB+wDxilqkuLORdL7n5ShVWvwWfXonVOYN95H7BjXw3S0/eyY8ceduzIID3dLbmP86/LzMwucLe1a0cX+EVwzDG5j2vlfVnUq1cj8NE/xoS50va5f6Oqpxy2bqmqdijmfd2ADOC1QpJ7P+BGXHJPBv6tqslFBoMl93Jh42yYej7UbgoXzYKYRsW+RVX57bcD7Nixp8gvgPyPc3KO/N0UEerWrfGHLqGC/jrI/bKoW7e6/3V9jAmR0hYOExE5I/ciqoicTmDT880TkeZFbDIAl/gV+EpE6opIQ1XdEkBMxk9Ne8Lgj+HdfjC5O1w0G2o3KfItIkKdOtWpU6c6xx9fv9hD5OTksHPnPn799Y9/FRz+RbBmzRbS0/eyc+feAvcTGRnBMcfUPOJL4MgvBvdlULt29NHfDGZMORJIcr8S+J+I1PGe7wKuCMKxGwOb8j1P89aV7+Ru//EPsx6uahr0vUYAsd7SsrQ7Syt1OMYUrRwOTCk2uatqCpCYm9xVtcwn5RSR0cBogKZNg59IjDEm3AQyzV4DEXkZeFNVd4tIGxG5MgjH/hnI/7d8grfuCKr6oqp2UtVO8fHxQTi0McaEt0CuNI3HTc6Re9XsO2BsEI49DRghzqnA7grR365qS0HLtlR4Ng6ebwg7VvkfT4iXd99ZQqOGY3ni8Y98j8WWcrCUQ4Ek9zhVfQtvsg5VzQIKHtOWj4hMAhYCrUQkTUSuFJFrReRab5MPgQ3AeuAl4PqjOQFTTsS3h6FzQXPgrbNgxwq/IwqpQYM6MmRIJ5566hMWLdrgdzjGHCGQoZBzcZNjf6qqHbxW9qOq2r0M4juCDYUs535dC2/3hKyDcNFnYV06OCPjAL17P0lmZjaffnobdevanLSmbBU1FDKQlvstuC6U40VkAfAabny6MUc6phUM/dxNvv12T9ia4ndEIRMTE82zzw5n69bd3HnnW1S0Uh4mvAUyXn0p0B04HbgGaKuq34Y6MFOB1TsBLv4cqtaBt3vBlkV+RxQyp5zSjNtvP5fp01OZNCl8z9NUPIGMlqkB3AWMVdUVQHMR6R/yyEzFVqeFS/DV42DKOZD2hd8RhcwNN/TkzDNbMm7cVNat2+p3OMYAgXXLvAIcAk7znv8MPBCyiEz4qN3UddHUbAjv9oVNc/2OKCQiIiJ4+unLiI6O4vrrX+fgwSy/QzImoOR+vKo+BmQCqOo+wG7TNIGp1di14Gs3c+UKfvrM74hC4thj6/Cvfw1j5cqfefjhGX6HY0xAyf2QiFQHV/5XRI4HDoY0KhNeah4LQ+dA3RNgan/44WO/IwqJ3r3bMWrUmbz44ufMmbPa73BMJRdIcv878DHQREQmArOAO0IalQk/Neq7BB/bBt4f4CbjDkPjxp3PSSc15Oab32D79j1+h2MqsUBGy3wKDAJGAZOATqo6N7RhmbBUPdaVCI5PhGmDYN27fkcUdNWrV+X55y8nI+MgY8e+QU5Ojt8hmUoq0ELX3YFeQA+ga+jCMWEvuh4M+RQadIbpQ2HNZL8jCrpWrRry978PYM6cNfz3v/P8DsdUUoEMhXwOuBZYDqwArhGRZ0MdmAlj1erAkJnQ6DT48FJYNcHviIJuxIjT6dOnHQ8+OIPly63msCl7gbTcewJ9VPUVVX0FN3NSz9CGZcJe1Vpuwo+E7vDRCFgx3u+IgkpEeOKJi4mNjeH6619n3z4bg2DKViDJfT2Qv4h6E2+dMaUTVRMGzoBm58DMP8G3L/odUVDFxsbwzDOXsWHDdsaNm+p3OKaSCSS51wJWi8hcr4jYKqC2iEwTkWkhjc6Ev6gacOH70KIffHoNfBNePX5nnNGSMWN6MWnSIqZPX+Z3OKYSCWSavXtCHoWp3KpEwwXvwoyLYfYYyDkEHf/sd1RBc9ttfVmwYB233z6ZU05pSkLCMX6HZCqBQIZCfq6qn+Muph4DZOSu89YbU3pVqsH5b8OJQ2DuLfD1o35HFDRRUZE8++xwcnKUG26YQFZWsdMhGFNqhSZ3EZkhIu28xw1xyf0K4HURGVs24ZlKJTIKzpsEJ10C8++Chff7HVHQNGsWxyOPXMTixT/w1FOf+h2OqQSKarm38KpAAvwJN1nH+UAyLskbE3wRVeDc16HN5fDlPbBgXLmdxqykbPYmU5aKSu6Z+R73wk2Lh6ruwZtyz5iQiIiEPq9Auyvhqwdg/t1hk+AfemgwTZocw5gxE9i1a5/f4ZgwVlRy3yQiN4rIQKADrr4MXhGxqLIIzlRiEZHQ+0VIvBYWPwqf3xoWCT4mJprnnrvcZm8yIVdUcr8SaIurKXOxqu7y1p+Kq/FuTGhJBPR6Dk65CVL+BbNvdBNwV3A2e5MpC4UOhVTVbbiyA4evnwPMCWVQxuQRgR5PQWRVWPIE5GTC2c+7xF+B3XBDT+bP/45x46bSuXMLWrZs4HdIJsxU7P8hpnIQgW6PQfJf3F2sM6+EnIo9nNBmbzKhZsndVAwicMYDcNq9sHI8fDwScip2QrTZm0woWXI3FYcInP53OPNBWD0RPrgMsjOLf185ZrM3mVAJpORvgohMFZHtIrJNRN4RkYSyCM6YAiX/Bbo9Dt+95UoWZB/yO6JSsdmbTCgE0nJ/BZgGNAQaAdOx0TLGb51vgx7/hvVTYdpgyKq4JXXzz9508802e5MJjkCSe7xXyz3LW8YD8SGOy5jidbjJDZXcMAPevxAy9/sd0VHLnb1p7lybvckERyDJPV1EhotIpLcMB9JDHZgxAUm6Dnr/F36cCe+dD5kV965Pm73JBFMgyf0KYCjwC7AFGIKrNWNM+XDyldB3PGyaA+/2g0MZfkd0VGz2JhNMgZT8/UlVL1DVeFWtr6oXqurGsgjOmIC1HQHnToCfv4B3+sLB3/yO6KjY7E0mWAq9Q1VE7lDVx0TkGeCIAhiqelNIIzOmpFpf4soGf3AJvNMbBn0M0XX9jqrEcmdveuaZzzjrrJM4//wkv0MyFVBRMzHlDrpdUhaBGBMUJw6BiCiYfhFMORsGfwLVK97MRzZ7kymtQrtlVHW693Cfqr6afwEq7lUrE/5OGAAD3oMdK+DtnrBvu98RlZjN3mRKK5ALqncHuM6Y8uO4fnDhNNi5Ft7qAXu3+h1RidnsTaY0ippm71yvv72xiDydbxkPBFTUQ0T6ishaEVkvIncV8HpTEZkjIt+IyLci0u+oz8SYwzXvDQM/gN0/wFtnQcZmvyMqsT/O3vS93+GYCqSolvtmXH/7ASAl3zIN6FPcjkUkEngWOBdoA1wiIm0O2+xvwFuqegowDHiupCdgTJGa9oTBH8OeNJjcHX7b5HdEJfbgg272phtusNmbTOCK6nNP9frXTzisz/1dVd0ZwL67AOtVdYOqHgLeBAYcfhigtve4Du4LxZjgSugKg2fCvm3wVnfY/aPfEZVIrVpu9qZt236z2ZtMwALpc28uIlNEZJWIbMhdAnhfYyB/MynNW5ffvcBwEUnDzdF6YyBBG1NijU+HIZ/CgZ2uBb+rYk1QbbM3mZIKtHDY87h+9h7Aa8CEIB3/EmC8qiYA/YDXRY6cYkdERovIEhFZsn17xRv5YMqJhl3golmQmQGTu8HOdX5HVCI33NCTM89sybhxU1m3ruJdIDZlK5DkXl1VZwHi3a16L3BeAO/7GWiS73mCty6/K4G3AFR1IRANxB2+I1V9UVU7qWqn+HirWWZKoUEHGDrHlQme3B3SK04NdZu9yZREIMn9oNeaXiciY0RkIBATwPsWAy1FpIWIVMVdMJ122DYbgV4AItIal9ytaW5CK749XDzXTbb91lluPHwFYbM3mUAFktxvBmoANwEdgcuBEcW9SVWzgDHATNzdrm+p6koRuU9ELvA2uxW4WkRSgUnAKLWrRaYsxLaBiz+HiCow+SzYtszviAJmszeZQEhJc6k3xHGYqk4MTUhF69Spky5ZYhURTJDsXO/uYs3McKUKju3kd0QB2b//EOed9xQ7duxh1qw7iI+v5XdIxgcikqKqBf7SFnUTU20RuVtE/iMivcUZA6zHlQA2puKrd4JrwVetA2/3gs1f+R1RQKpXr8pzz9nsTaZwRXXLvA60ApYDVwFzgIuAgap6+Hh1YyquOi1cgq8R76pJpn3hd0QBOemkhtxzzwU2e5MpUFHJ/ThVHaWq/4cbstgG6KOqy8okMmPKUu2mMPRzqNkQ3u0Lm+b6HVFARo48w2ZvMgUqKrln5j5Q1WwgTVUPhD4kY3xSq7Frwddu5mZ0+ukzvyMqls3eZApTVHJPFJHfvGUP0D73sYhUzGlujClOzWPdOPi6J8DU/vDDx35HVCybvckUpKjaMpGqWttbaqlqlXyPaxf2PmMqvBr1XYKPbQPvD4Dvpxf/Hp/lzt40adIipk9f5nc4phwIZJy7MZVP9VhXqiA+EaYNgnXv+h1RsW67rS8dOjTj9tsnk5b2q9/hGJ9ZcjemMNH1XLGxBp1h+lBYM9nviIpkszeZ/Cy5G1OUanVgyExodDp8eCmsClbNvNCw2ZtMLkvuxhSnai0Y/BEkdIePRsCKV/yOqEg2e5MBS+7GBCaqJgycAc3OgZlXQOoLfkdUJJu9yVhyNyZQUTXgwvfhuPPgs+tg1o2udHA5ZLM3GUvuxpRElWi4YCp0vAWW/Qfe6lFuJ9622ZsqN0vuxpRUZBSc9SSc9yZsT4UJHSFtvt9RFchmb6q8LLkbc7ROuhguXeQuuL7dE5b+G8pZ94fN3lR5WXI3pjTi2sJli6HFeTBnLHw4HDL3+h3VH9jsTZWTJXdjSqtaHRjwLpz5IKyZBG+cWu4m37bZmyofS+7GBINEQPJfYPDH7gLrxM7lribNuHHnc9JJDbn55jfYvn2P3+GYELPkbkwwNe8Nw1OgzvHw3gWwYBzklI8yADZ7U+Viyd2YYKvTHIZ9AW3/BF894EoH7y8fhbzyz9700ks2e1M4s+RuTChEVYc+L8PZL8DGWW645NZv/I4K+H32pocestmbwpkld2NCRQQSr4Fh8yEnC948HVa+6ndUNntTJWHJ3ZhQa5gMl6dAw9Pg41Hw2fW+ly2w2ZvCnyV3Y8pCjfow5BPodDukPg+Tu8Oen30NyWZvCm+W3I0pKxFVoPtjcP7bsGMFTOgAmz73NSSbvSl8WXI3pqydOAQu+xqq1YO3e8GSf/pWtsBmbwpfltyN8UNsa5fgTxgAn98KM4bBoQxfQrHZm8KTJXdj/FKtNpw/Bbo+AuumwBvJ8OtaX0IZNKgjgwd3tNmbSqg818mX8hxcQTp16qRLlizxOwxjguunWfDBMMg+CH1fg5YXlnkIe/YcoHfvJ8jMzOazz26nbt0aZR5DeXbwYBarV29m2bKNLFu2idTUjQwc2IGbbjrHt5hEJEVVOxX0WpWyDsYYU4BmvVzZgmmDYdpA6HI3nHE/RESWWQi5szcNGPA0d9zxFv/3fyMRkTI7fnmSnZ3D+vVbWbZsk5fMN7J69WYOHXLXJGJjY0hKakLTprE+R1o4S+7GlBe1m7obnmbfCF8/DFuXQL83oEZcmYWQO3vTww9/wKRJi7j00lPL7Nh+UVU2bkzPS+SpqZtYvjyNvXvdzV0xMdVITGzCVVd1JympCUlJTWncuF65/+KzbhljyqNv/wuzb4CaDeGCd6BBxzI7dE5ODsOGvUBKyk98/PEttGzZoMyOXRa2bfstr2slN5nv3Olq8FerVoW2bRuTmOiSeFJSU44/Pp6IiPJ5ebKobhlL7saUV78shmlDYN9W6PUcnHxF2R36l9306vU4jRrVZcaMsVSrVjH/yN+9ez+pqS6B5yb0LVt2ARARIZx0UsO8RJ6Y2ISTTmpI1aoV51wtuRtTUe3bAR9cAhs/g/ajocfTUKVamRz6k09WMGrUy4we3Z17772wTI5ZGvv2HWLlyp9JTd3IN9+4hL5hw/a811u0iCMxsSlJSU1ITGzKySc3pkaNsvksQ8UuqBpTUdWIcxOALPgbfP0IbFvmhk/WbhLyQ+efval791b06NE65McMVGZmNmvWbCE19ffulbVrfyE729Wob9iwDomJTRk6tDOJia5VXtlG/4S05S4ifYF/A5HAf1X1kQK2GQrcCyiQqqqXFrVPa7mbSmvdVPh4JERGQ/83oWnPkB9y//5DnHfeU+zYsYdZs+4gPr5WyI95uJycHDZs2JGvRb6RlSs3c+BAJgB169b4Q9dKUlJTjj22TpnH6QdfumVEJBL4DjgHSAMWA5eo6qp827QE3gJ6qupOEamvqtuK2q8ld1Oppa+BaYNg51p381On21xp4RBas2YL/fr9i1NPPZ4JE64O6cVFVeXnn3d5LXLXtZKauok9ew4Abjap9u0T/pDImzWLLfcjV0LFr26ZLsB6Vd3gBfEmMABYlW+bq4FnVXUnQHGJ3ZhKL/YkuGwRzLwS5t0BWxZB31egauha1LmzN/3lL+/w0kvzuOaas4K27/T0jHwXO92yY4crwxAVFUnr1g0ZOLBDXl95y5YNqFKl7Mb+V2ShTO6NgU35nqcByYdtcyKAiCzAdd3cq6ofH74jERkNjAZo2rRpSII1psKoWgv6T4aUZJh3J0xcCRdMdYk/REaOPIPPP1/LQw/N4PTTT+DkkxNKvI+MjAN8+22a1yJ3feWbNrlKlCJCy5b16dmzdV4ib926EdHRUcE+lUojlN0yQ4C+qnqV9/xyIFlVx+TbZgaQCQwFEoB5wMmququw/Vq3jDH5bJwDMy6GrP3QdzycODhkh0pPz+Ccc56gZs1qzJx5S5EjTQ4ezGLVqp/zbtP/5puNrF+/La8WS5Mmx+SNWklKakr79gnExESHLPZw5Ve3zM9A/kv6Cd66/NKARaqaCfwgIt8BLXH988aY4jTtAZcvdWULpg+BznfAmQ+62vFBljt709ChzzNu3FSefHIY4G7V/+67X/ISubtVfwuZme5W/fj4WiQmNmHAgFPy+spjY2OCHp/5o1Am98VASxFpgUvqw4DDR8K8B1wCvCIicbhumg0hjMmY8FMrAS6eB3NuhsWPubIF570JNeKDfqjc2ZueeeYzMjOz2bTpV779No39+920gbVqRZOY2ITRo7vn3eHZqFHdSnvB00+hHgrZD3gK15/+P1V9UETuA5ao6jRxP/Engb5ANvCgqr5Z1D6tW8aYIqx4BT67zk3rd8E7cGznoB8iMzObIUOeZfnyNNq2bZxXbyUxsSnHHRdXbm/VD0d2h6oxlcnWpW645N4t0PM/0P7qoB8iOzsHVbWRKz4rKrnbV6wx4aZBB1c+OOEs+HQ0zLwKsg4E9RCRkRGW2PduhYzNfkdRKEvuxoSj6rEw6ENI/iuseBne7Aq//eR3VBXb/nRY9y7MGgPj28ILx0LKv/yOqlBWW8aYcBURCWc+4PrdPxoBr3d0ZQuane13ZBXDwd/g5/mwcbZbtqcCClVqQEJXaDMSjuvnd5SFsuRuTLg7YQBcttj1w7/TB854ELrcGfKyBRVO5j74eQFsmgObZsMvS0CzIbIaNDodzrgPmvRwX5aRVf2OtliW3I2pDI45ES79Cj65Cr64G3752t30VK2235H5J+sg/LLo95b5lq8gJ9PdI3BsF0i+G5r0hEanQZWKd4OVJXdjKouqMXDeJGh4Knx+G0zsDAOmQmwbvyMrGzlZrjW+aba7s3fzAndnr0RA/Q7Q8c+uZd74TPdZVXCW3I2pTESg41g3omb6UJjYBfq8Aq0u8juy4MvJdv3km+a4lvnP8+HQHvdafHs3+UmTnpDQDaLr+hpqKFhyN6YySujmhktOvwhmDIUtt0K3R0JStqDMqEL6KpfIN82BtLlwYKd7rV4raD3ctcybnBWSu3fLmwr8kzTGlEqtxnDxXJh7C6Q8CdtSXLXJGvX9jiwwqrBr/e8t801zYJ9XNbxOCzhhoJvQpEkPiGnkb6w+sORuTGUWWRV6/cddQPzsGni9g5vGr9GpfkdWsN82/p7IN86GjDS3PqYRNOvtEnnTHi65V3KW3I0x0HaE64eeNggmd4OeT0P7a/wfLrn3F3fxc5M3omW3V1ewepyXyHu6fvN6Lf2PtZyx5G6MceonwWVL4KPhrvjYlq+g1/MQVb3sYtifDpvm/t4y/3W1W1+tjiun0OFml9Tj2rpRLqZQltyNMb+rfgwMnAFf/gO+ug+2f+uqS4aqm+Pgbkib//vwxNy7QKNqQuOu0O5PrnUen+TuuDUBs+RujPkjiYAz/uGVLRgOEzrBeW9A8z6l33fm3t/vAt0429We1xx3F2jjM+CM+/PdBWpT7JWGJXdjTMGO7++6aaYNgnfOdYk3+e6SdYdkHXTdOxtnu9b5lkW/3wXa8FRX2KxpT/e4At4FWp5ZcjfGFK7eCXDpQvhkNCz4m1e24NXCb/rJznSt8dwRLZsXuHLDEgENOkLHW9xolkZnhMVdoOVZhZusQ0S2A37XLo0Ddvgcg58q+/mDfQZgnwH4/xk0U9UC78iqcMm9PBCRJYXNflIZVPbzB/sMwD4DKN+fgY0lMsaYMGTJ3RhjwpAl96Pzot8B+Kyynz/YZwD2GUA5/gysz90YY8KQtdyNMSYMVdrkLiL/E5FtIrIixMcZLCIqIuXyirqIRIrINyIyIwT7ThKRr0RkmYgsEZEuwT7G0RCRP4vIShFZISKTRCSod8+IyEXe/nMO/7mLSHsRWei9vjzYxy4ipkJ/30XkRhFZ48X0WJCP20xElnq/AytF5FpvfQ0R+SDfcR8J5nELiCNaRL4WkVTveP/I99pEEVnr/T78T0SCfmtscT93EZkW9FykqpVyAboBHYAVJXhPvRIeoxYwD/gK6OT3ORcS4y3AG8CMYH8GwCfAud7jfsDccnC+jYEfgOre87eAUUE+79ZAK2Bu/p877qbBb4FE73ksEFlG513g7zvQA/gMqOY9rx/kz6Jqvn3HAD8CjYAaQI9828zP/V0J0fkLEOM9jgIWAafm+90Ub5kEXBfkz6DInzswyPs/GHAuCmSptC13VZ0H/FrCt73nfcNeICKB3N17P/AocKDEAZYBEUkAzgP+W4K3LfFaOj1Fiq2xqkDuDMx1gM1HEWYoVAGqez/DGgQWV8DnraqrVXVtAS/1Br5V1VRvu3RVzS5p8EejiN/364BHVPWgt922AHZ3sdfKvVVEipzSSFUP5e4bqIbXW6Cq+1R1Tu42wFIgIbCzKTl1MrynUd6i3msfeq8r8HWAcdzu/SVwjYgUN8t4oT93EYnBNbAeKPlZFa3SJvejdBbwT2AIsFpEHhKREwraUEQ6AE1U9YMyjK+kngLuAHJK8J4Tca2bMcAqEfmLiBQ2zc1Y4HER2QQ8Adx99KEGh6r+7MWyEdgC7FbVTwJ4a0nOu6h9qIjM9Loq7ijh+0PhRKCriCwSkc9FpHNxb1DVF4BzcV+M80Rkioj0FSm46IyINBGRb4FNwKOquvmw1+sC5wOzSnkuRfK6IJcB24BPVXXRYa9HAZcDHxe3L1X9i7ftccBSEXlFRM4sZPOifu73A08C+0p8QgEEWWkXoDlH+acQrkX6IJAFDD7stQjcn+TNvedzKWfdMkB/4Dnv8VkE2C1z2D7ica3+LKBLAa8/nfvZAEOBz8rBedcDZnuxRwHvAcODed75tvvDzx24DdclFIdLjAuBXmV47kf8vgMrgGdwXRJdvPikBPsUXLfGZmBaMds2wrWMG+RbVwX4CBhbhp9DXWAO0O6w9S8BTx3F/iKB4cBu4OkCXi/w5w4k5X5mpclFhS3Wci9E7re8t9yXb311EbkUeBfoA9wMfHrY22sB7YC5IvIjcCowrZxdVD0DuMCL702gp4hMyL+B1+LK/Qyuzbe+johcA0wDWgJX4PoUDzcS9zkBvI1LHn47G/hBVberaiYuvtPzbxCE8y5MGjBPVXeo6j7gQ1w/uJ/SgHfV+Rr3V1xc/g1E5MHcz+Ow9V2A53Bf4m9RzF9m6lrsK4Cu+Va/CKxT1adKeR4BU9VduOTeN3ediPwd96V9S0Hv8Vrmy0Tkw3zrRER6Aq8C9+A+hycLeHthP/fTgE7e/8EvgBNFZG6pTzBXWX1blseFEn5bAo/hvoH/A5xSgvfNpZy13A+L7ywCv6A6AfgeeARoWcy2q4GzvMe9gJRycK7JwEpcC0pw/zFvDOZ5F/Zzx/3VsNQ7dhXchczzyvDcj/h9B64F7vMen4jrOimy5Y7Xh4y7YD4UqFrEtgn8fvG6HvAdcLL3/AHgHSCiDM49HqjrPa6Ou4Db33t+FfBlbpwB7u8yYK0Xfz+KuDAeyM+9pLkooBjL6hervC24/tMtQCbum/XKAN7TD4g+imOFU3K/AKgS4LZnAilAKm50Qke/z9WL6x/AGlwr8nW80RxBPO+B3u/UQWArMDPfa8O9L5cVwGNleM4F/r7jRqpM8OJZCvQMYF8dcdUIAznuOd4XQar372hvfQLuguZqYJm3XBXC828PfOPFsAK4J99rWbgv7tw47glgf2cC8SU4fpE/91Akd7tD1RhjwpD1uRtjTBiy5G6MMWHIkrsxxoQhS+7GGBOGLLkbY0wYsuRuwoqINBCRN0Rkg4ikeJX4BnqvnSUiu72bUVZ7N64gIqNE5D+H7WduObvpzJgSseRuwoZX0Os93N2Ax6lqR2AYfywENV9Vk4BOwHCvBpAxYceSuwknPYFD6gpbAaCqP6nqM4dvqKp7cTdYFVj4rTAi8qOIxHmPO+XeLi4i3fOVLPhGRGqJSIyIzPKKRS0XkQH59jNOXA3xL8TVlL/NW3+8iHzs/dUxX0ROOpoPwphAytYaU1G0xd1lWSwRicXV/LkfKLYSYgBuA25Q1QVeGdfcMs8DVfU37wvhKxGZhvurYTCQiCtethT3RQOu1sq1qrpORJJxtVt6BiE+U8lYcjdhS0Sexd0mfkhVcxN4VxH5Blcg6xFVXVlE33pJbt9eAPxTRCbiCnGleSVkHxKRbt7xGgMNcEXb3lfVA8ABEZnuxRuDK2L2tvxeMr5aCWIwJo8ldxNOVuJaxACo6g1ei3lJvm3mq2r/w96XjivulN8xwI4CjpHF792ZeVOlqeojIvIBrv7QAhHpg/vLIB5XUyfTq/5X1LR6EcAu75qAMaVife4mnMwGokXkunzragTwvsXAGSJyLLi+dFyLeVMB2/6IK5wF+b5IROR4VV2uqo96+zsJN/vUNi+x9wCaeZsvAM4XN69nDK62Pqr6G/CDiFzk7VNEJDGA+I05grXcTdhQVRWRC4F/ebPdbAf2AncW876tInIz8KE3m1AGcImqFjRD1T+Al0Xkfly1z1xjvQSeg/sL4iNcXf/pIrIc99fDGu94i72+929xVSOX4yZ6AFdK9nkR+RuuP/5NXEVFY0rEqkIa4wMRiVHVDBGpgZtEfbSqBnQx2JhAWMvdGH+8KCJtcH3wr1piN8FmLXdjjAlDdkHVGGPCkCV3Y4wJQ5bcjTEmDFlyN8aYMGTJ3RhjwpAld2OMCUP/D4LzvOpD9rk/AAAAAElFTkSuQmCC\n",
- "text/plain": [
- "<Figure size 432x288 with 1 Axes>"
- ]
- },
- "metadata": {
- "needs_background": "light"
- },
- "output_type": "display_data"
- }
- ],
- "source": [
- "# plot the computing time\n",
- "print(tot_time_hpc2)\n",
- "plot_computation_time(tot_time_hpc1, tot_time_hpc2, ngpus_sort, names_hpc, \\\n",
- " \"./total_computation_time_{0}_vs_{1}\".format(*names_hpc), log_yvals=False)\n",
- "\n",
- "plot_computation_time(overhead_hpc1, overhead_hpc2, ngpus_sort, names_hpc, \\\n",
- " \"./overhead_time_{0}_vs_{1}\".format(*names_hpc))\n",
- "# plot speed-up factors\n",
- "plot_speedup(tot_time_hpc1, tot_time_hpc2, ngpus_sort, names_hpc)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 113,
- "metadata": {},
- "outputs": [],
- "source": [
- "## evaluate iteration time\n",
- "# get iteration times\n",
- "iter_data_hpc1 = get_iter_time_all(base_dir, wildcard_hpc1, gpu_id_str=gpu_id_str[0])\n",
- "iter_data_hpc2 = get_iter_time_all(base_dir, wildcard_hpc2, gpu_id_str=gpu_id_str[1])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 114,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Saving plot in file: boxplot_iter_time_Juwels_vs_Booster.png ...\n"
- ]
- }
- ],
- "source": [
- "# plot the iteration time in box plots\n",
- "boxplot_iter_time(iter_data_hpc1, iter_data_hpc2, ngpus_sort, names_hpc)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 115,
- "metadata": {},
- "outputs": [],
- "source": [
- "def get_slowiter(iter_time, threshold):\n",
- " inds_slow = np.where(iter_time > threshold)[0]\n",
- " return iter_time[inds_slow], np.shape(inds_slow)[0]\n",
- "\n",
- "def ana_slowiter(itertime1, itertime2, thres, names):\n",
- " slowt1, nslow1 = get_slowiter(itertime1, thres)\n",
- " slowt2, nslow2 = get_slowiter(itertime2, thres)\n",
- " \n",
- " if nslow1 > 0:\n",
- " print(\"{0:d} slow iteration steps on {1} with averaged time of {2:5.2f}s (max: {3:5.2f}s)\"\\\n",
- " .format(nslow1, names[0], np.mean(slowt1), np.max(slowt1)))\n",
- " else: \n",
- " print(\"No slow iterations on {0}\".format(names[0]))\n",
- " \n",
- " if nslow2 > 0:\n",
- " print(\"{0:d} slow iteration steps on {1} with averaged time of {2:5.2f}s (max: {3:5.2f}s)\"\\\n",
- " .format(nslow2, names[1], np.mean(slowt2), np.max(slowt2)))\n",
- " else: \n",
- " print(\"No slow iterations on {0}\".format(names[1]))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 116,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "***** Analyse single GPUs experiments *****\n",
- "1 slow iteration steps on Juwels with averaged time of 5.18s (max: 5.18s)\n",
- "No slow iterations on Booster\n",
- "***** Analyse 4 GPUs experiments *****\n",
- "No slow iterations on Juwels\n",
- "No slow iterations on Booster\n",
- "***** Analyse 8 GPUs experiments *****\n",
- "No slow iterations on Juwels\n",
- "No slow iterations on Booster\n",
- "***** Analyse 32 GPUs experiments *****\n",
- "No slow iterations on Juwels\n",
- "No slow iterations on Booster\n",
- "***** Analyse 32 GPUs experiments *****\n",
- "No slow iterations on Juwels\n",
- "No slow iterations on Booster\n",
- "***** Analyse 64 GPUs experiments *****\n",
- "No slow iterations on Juwels\n",
- "No slow iterations on Booster\n"
- ]
- }
- ],
- "source": [
- " \n",
- "## settings\n",
- "names = [\"Juwels\", \"Booster\"]\n",
- "slowiter_time = 5. # arbitrary threshold for slow iteration steps\n",
- "\n",
- "# analyze single GPU experiments\n",
- "print(\"***** Analyse single GPUs experiments *****\")\n",
- "itertime_juwels = iter_data_hpc1[\"1 GPU(s)\"]\n",
- "itertime_booster = iter_data_hpc2[\"1 GPU(s)\"]\n",
- "\n",
- "ana_slowiter(itertime_juwels[1:], itertime_booster[1:], slowiter_time, names)\n",
- "\n",
- "# analyze 4 GPUs experiments\n",
- "print(\"***** Analyse 4 GPUs experiments *****\")\n",
- "itertime_juwels = iter_data_hpc1[\"4 GPU(s)\"]\n",
- "itertime_booster = iter_data_hpc2[\"4 GPU(s)\"]\n",
- "\n",
- "ana_slowiter(itertime_juwels[1:], itertime_booster[1:], slowiter_time, names)\n",
- "\n",
- "# analyze 8 GPUs experiments\n",
- "print(\"***** Analyse 8 GPUs experiments *****\")\n",
- "itertime_juwels = iter_data_hpc1[\"8 GPU(s)\"]\n",
- "itertime_booster = iter_data_hpc2[\"8 GPU(s)\"]\n",
- "\n",
- "ana_slowiter(itertime_juwels[1:], itertime_booster[1:], slowiter_time, names)\n",
- "\n",
- "# analyze 16 GPUs experiments\n",
- "print(\"***** Analyse 32 GPUs experiments *****\")\n",
- "itertime_juwels = iter_data_hpc1[\"16 GPU(s)\"]\n",
- "itertime_booster = iter_data_hpc2[\"16 GPU(s)\"]\n",
- "\n",
- "ana_slowiter(itertime_juwels[1:], itertime_booster[1:], slowiter_time, names)\n",
- "\n",
- "# analyze 32 GPUs experiments\n",
- "print(\"***** Analyse 32 GPUs experiments *****\")\n",
- "itertime_juwels = iter_data_hpc1[\"32 GPU(s)\"]\n",
- "itertime_booster = iter_data_hpc2[\"32 GPU(s)\"]\n",
- "\n",
- "ana_slowiter(itertime_juwels[1:], itertime_booster[1:], slowiter_time, names)\n",
- "\n",
- "# analyze 64 GPUs experiments\n",
- "print(\"***** Analyse 64 GPUs experiments *****\")\n",
- "itertime_juwels = iter_data_hpc1[\"64 GPU(s)\"]\n",
- "itertime_booster = iter_data_hpc2[\"64 GPU(s)\"]\n",
- "\n",
- "ana_slowiter(itertime_juwels[1:], itertime_booster[1:], slowiter_time, names)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Summary\n",
- "- Occasionally, a few iteration steps are slow\n",
- "- However, performance degradation seems to be much worser on Booster than on Juwels\n",
- "- Higher chance for slow iteration steps on Booster in general"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 157,
- "metadata": {},
- "outputs": [],
- "source": [
- "def boxplot_iter_total_time(iteration_time, total_time, ngpu_list, name, log_yvals=False):\n",
- " nexps = len(ngpu_list)\n",
- " bar_width = 0.35\n",
- " # create data lists for boxplot-routine\n",
- " iter_time_all = []\n",
- " for i in np.arange(nexps):\n",
- " iter_time_all.append(iteration_time[\"{0} GPU(s)\".format(ngpu_list[i])])\n",
- " \n",
- " # trick to get list with duplicated entries\n",
- " xlabels = [val for val in ngpu_list for _ in (0, 1)]\n",
- " nlabels = len(xlabels)\n",
- "\n",
- " # Multiple box plots on one Axes\n",
- " #fig, ax = plt.subplots()\n",
- " fig = plt.figure(figsize=(6,4))\n",
- " ax = plt.axes([0.1, 0.15, 0.75, 0.75]) \n",
- " \n",
- " bp = ax.boxplot(iter_time_all, positions=np.arange(0, nlabels, 2), notch=0, sym='+', vert=1, showfliers=False, widths=bar_width) # Outliers for initialization are disturbing\n",
- " ax.set_xlabel('# GPUs')\n",
- " ax.set_ylabel('Time [s]')\n",
- " \n",
- " # Reference: https://matplotlib.org/3.1.1/gallery/statistics/boxplot_demo.html \n",
- " num_boxes = len(iter_time_all)\n",
- " medians = np.empty(num_boxes)\n",
- " for i in range(num_boxes):\n",
- " box = bp['boxes'][i]\n",
- " boxX = []\n",
- " boxY = []\n",
- " for j in range(5):\n",
- " boxX.append(box.get_xdata()[j])\n",
- " boxY.append(box.get_ydata()[j])\n",
- " box_coords = np.column_stack([boxX, boxY])\n",
- " ax.add_patch(Polygon(box_coords, facecolor=colors[1]))\n",
- " # Now draw the median lines back over what we just filled in\n",
- " med = bp['medians'][i]\n",
- " medianX = []\n",
- " medianY = []\n",
- " for j in range(2):\n",
- " medianX.append(med.get_xdata()[j])\n",
- " medianY.append(med.get_ydata()[j])\n",
- " ax.plot(medianX, medianY, 'k')\n",
- " medians[i] = medianY[0]\n",
- " # Finally, overplot the sample averages, with horizontal alignment\n",
- " # in the center of each box\n",
- " ax.plot(np.average(med.get_xdata()), np.average(iter_time_all[i]),\n",
- " color='w', marker='*', markeredgecolor='k', markersize=10)\n",
- " \n",
- " ax2 = ax.twinx()\n",
- " x_pos = np.arange(1, nlabels+1 ,2)\n",
- " \n",
- " ytitle = \"Time [min]\"\n",
- " max_time = np.max(total_time)\n",
- " time_order = val_order(max_time)\n",
- " ymax = np.ceil(max_time/(10**time_order) + 0.5)*(10**time_order) + 10**time_order\n",
- " # np.ceil(np.maximum(np.max(times1)/100. + 0.5, np.max(times2)/100. + 0.5))*100.\n",
- " if log_yvals: \n",
- " total_time = np.log(total_time)\n",
- " ytitle = \"LOG(Time) [min]\"\n",
- " ymax = np.ceil(np.max(total_time) + 0.5)\n",
- " \n",
- " # create data bars\n",
- " rects = ax2.bar(x_pos, np.round(total_time, 2), bar_width, label=names, color=colors[0])\n",
- " # customize plot appearance\n",
- " # Add some text for labels, title and custom x-axis tick labels, etc.\n",
- " ax2.set_ylabel(ytitle)\n",
- " ax2.set_xticks(np.arange(0, nlabels))\n",
- " ax2.set_xticklabels(xlabels)\n",
- " ax2.set_xlabel('# GPUs')\n",
- " ax2.set_ylim(0., ymax)\n",
- " \n",
- " # add labels\n",
- " autolabel(ax2, rects, rot=45) \n",
- "\n",
- " plt_fname = \"iter+tot_time_{0}_vs_{1}\".format(*names)\n",
- " print(\"Saving plot in file: {0}.png ...\".format(plt_fname))\n",
- " #plt.show()\n",
- " plt.savefig(plt_fname+\".png\")\n",
- " plt.close()\n",
- " \n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 158,
- "metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "Saving plot in file: iter+tot_time_Juwels_vs_Booster.png ...\n"
- ]
- }
- ],
- "source": [
- "boxplot_iter_total_time(iter_data_hpc2, tot_time_hpc2, ngpus_sort, names_hpc[1])"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.8.5"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
--
GitLab