From 327e3011c1aea6e8affdc5ad2f1c9816e832c6ab Mon Sep 17 00:00:00 2001
From: Alex Phillips <A.M.Phillips@liverpool.ac.uk>
Date: Thu, 17 Aug 2023 12:01:09 +0100
Subject: [PATCH 1/2] test: add tests for resampling

---
 test/generate_test_data.py                    |  13 +++
 .../eight_schools_resampled_samples.npy       | Bin 0 -> 32128 bytes
 test/test_resample.py                         | 108 ++++++++++++++++++
 3 files changed, 121 insertions(+)
 create mode 100644 test/test_models/eight_schools/eight_schools_resampled_samples.npy
 create mode 100644 test/test_resample.py

diff --git a/test/generate_test_data.py b/test/generate_test_data.py
index db4a535..59e1682 100644
--- a/test/generate_test_data.py
+++ b/test/generate_test_data.py
@@ -52,3 +52,16 @@
 log_weights = new_logProbs - logProbs
 log_weights = log_weights - logsumexp(log_weights)
 np.save(os.path.join(model_path, "eight_schools_log_weights.npy"), log_weights)
+
+nsamples = samples.shape[0]*samples.shape[1]
+tmp_samples = samples.reshape((nsamples, 1, samples.shape[2]))
+tmp_log_weights = log_weights.reshape((nsamples))
+
+rng = np.random.default_rng(seed=0)
+resampled_iterations = rng.choice(
+        nsamples,
+        size=nsamples,
+        p=np.exp(tmp_log_weights))
+
+resampled_samples = tmp_samples[resampled_iterations, :, :]
+np.save(os.path.join(model_path, "eight_schools_resampled_samples.npy"), resampled_samples)
diff --git a/test/test_models/eight_schools/eight_schools_resampled_samples.npy b/test/test_models/eight_schools/eight_schools_resampled_samples.npy
new file mode 100644
index 0000000000000000000000000000000000000000..3d07bd5184e90d6586c66a681fe6b51a593ccfe1
GIT binary patch
literal 32128
zcmb_l1#=Zi*N(g08~3Ym7a-U?xCeK4cU|1wA-KCc1iyp3JHahjump*}r@Qn0fv+aD
zRa;~?A(`p!bIxPED79(51})9ahMA2D?b3B{=Ru)mHKFCZl?jc~gm&vcXvm<B{o3~*
z)TJx?eD#if26v^;5AM-%Kv(*6SW;YElqR13i3^X?jM4l*J{%Ufjx4paq8&=LC{ru@
z(hF(mT=OSi#;9fE^+Sb`kwvBjD_`4;Zl^}Wr~acmCK%=ACzi~2%`1{Bjk&feW~Vn$
zFS_-YSwE%1X*n5>7HE+F@qXH^Mn3q~s&d7TDIusidrj+)wI4_sGYh-^eHVt^8+X`m
zDAHozn<l+GdVMjy2-x!==+7r<6Mg!dUV|8`l(3Fl)Ts5S?&&96RPg$9=v;euU-VWr
ztoy?3q3P|}3G<_55Bb;KdY#-J=bL<+{jC0>p_yF4vGRYxmp(}2(=qYQG`)P|I(_ts
z!8KLF{O(ILrPqgERcqf-hy7>H)><K5maGQn{5v&AE&t4^-EwM~P#j-(x@)R;fc$cc
zl)s>|yWA?Ge~mKFwP+d_vE6;Xt^D_0$Mp|$qOslmcGbrX;;}ogP2US+L$P3Ne%S%3
zpG@Z~@9z0R6^Qar%XEr9Z7y3Cx4trWj=%hA%-i0#cKwn@zH_x%lmFKwHE+G&^fgbK
z*La9wX|Wqt7|nb2dSWF9<xZNA`NAH)R?S?Z?}p&atkV;d^v>vBwP*hy$>m|&epy}5
zRf+f-QaQ;=wO^W;`Qn^$nj5~nDdi`RGU8I5(laW(ibm^uC*OLkj=|ZU^c^Z|5YT>9
zw^|xEY}?+o_n4pl=+feX|MIPBR82ZQ>2Zo0&%ZbD_}E#8Wxl1Z=G{1H8gu(#%ik&P
zc$U*$YP888z1PmY<aFN(hvx1GnZCjs2FGRc?hHRnlCQ3>_RtGuqDS0m{=*ISWhO@?
zxAwqn`msC5Yq7HF+MT<7e6jlYua;{r1!3Hi3jejRQv-iS{=GEH6Q^c(I=1;ip|twZ
zW8;qWzf!|)lj<g@edVbdTjTZne(1LLQrNANAEmc*uItCm*I<Ixg?aJawczj7$sN;a
zxlOod!0v{wcRsoV<64F))v;Xwq~vzTYt2p&e_xmvj?!|t^tcrK<kp%>3!>yj$)B};
zh3*(|bNbyEaq+Uf*R({_XT6L~FHgHnu7qT#U{&KzR(R9GJN?e9a-i?#WDyyFm732L
z&93;P*QU+ChMo6BsT<`|@4VLI*c_W5USWZ-nUGdx)~+C&|B`2Y;&m|Mzg%dLx>$p<
zh^sNLz6K3FKDz(6Fh^=%egCt&z1-2{qUqLAH#PXaS?>!6_8=T|fP2St_jVdFIOamr
zuE8GgpD?1WX}mXjk1g0-ZOj8x&K8~LMK343d)grWW_vp<cTTu+xur8Ar+GcOnOO>+
zZwu;$RrNxH%QrI~e|;~-FLuw0IHW0gPTnwoVPQWrxY)$kZA>^SO0(U9oK@IL%U*Xt
zk6EP?<`q;5z``~~Rr0bD@H=|&@u+S-xH;#OeVvmoa^lu&ojSP(W9H7Qrz0okNyT4l
z%)IOt1YLu|Ni}WWNY}%bv~HE<hkWZ39WyKf6kjig`|OTYxAMjB^^Tug+o(Z@7wsrN
zSe&$0d#v>Uol3K2INM;<^y<UiMrjarf5nHu?rymKtNhV&U!3q`>%}^yQ#6Q+l`EDV
zS|F{PvG4V#w9BSS_s5tuI^iap^{15<8-j`-pX|yDFyPC!UK6*qR-uLY>5P*DUz%JG
zSyXA39*ivdwf9B%;eO3J@3Sszp{3=q{YDVZHvIKDK<|b!|IT*V5myHM*i4JT$9mNw
zX4mtUJM6TYb?%7|vUjE+E>kP3y3gy_zF@A&-}zPae6tV~PJF(vugbRM>mo6))Z7-c
zZduB+Cbxe3I5^ETq2hwjTEjiSfBen1@-X%VSvK+{kMdL6Y=144d)Mu#`AZ-9`R9~R
zV=U=;MDLz@zJj~Vj=OrrDzQ)fz0tnXiv2yxx+BV^``F~tI;3?NSSEX?8k^2Ot~`IR
z4oCH?yiSJcaAToa^5I=B$k_Te-ZCu+$GiQyQEI6hLV_BUdVf`gjrwB+({jA9vvq~g
zQ+plc%gb1SD3UHa(F(8c1wU3e&O{=#kL3II+*xVN$BB7MYU!ZeH>k&g=p<C!<4jA|
zORSf(-q7A?k+izd4d%)8rf&tKa)@t=<4R9BeVkIxyLzUncjTX?&+plzFD-+o3Fhef
zF6wK~X)ZY4+dg+)jR>r4b>yLUcp!$@?oXU@(-nOK_Is_h_+^^*;?3&{)4oWF=l93-
z+gN0(+HhU{cQfqdrOm3m3992L_uA}s**N#1w63)0_P8VeO!O8l52UN*jK1D0z8v~(
z;`{rn8B)`~x)tx+E|nfQ$W%q;j_8Wk@`w;`&{MRv3o(=P_t+JsFEod)U#ig(7o_hW
zw!0k5^_KtI&;g$8iRNEpPlV3<VJf;m;!{Y&2vHAR0#Nn2t>$n~ckHWh(13b{CQYTE
zov*k3E`{t!w>dJ?7Kdjn$`5of6Z^8S9m=QPS<ogZ6hT`*){KeO<LT?YSNwjvq8t0L
ziI9`of?nt@p8teUu;X`Inof3^+{N7eXpZS_;Q5;#dG>NS^Lmq==NXVWJ1on2Zah{c
zZa7qWS*)xoRsK|?VFo08*kgRRu?)Jt^gL61M2!Vmqow52k?5C1mG*8N&UM+*p~>AS
zB+-}sJ8Qt)p#%M$VU(8~a2mMFQ9>g7=v__J&3#cQ{4?~`s(~T+Y?Y)cRS*i>dvT+y
z*0lgTP8U2g;sCuJ^OJefQOVoc*ELI`syDlKuWYGd?-j#Wy-9>yi<Co`my{NNUne7q
zLOgV@ZrjP-c!4&Tf4j~3ruma=YJQ-z-py4d7cXqGWtrtSX^3-T(i*!brWutwyuaPa
zMn3&`)X!@9Z%tdZJAQN@@zzx5Mo<sysfAMGk0&q24AkIywb7y7+InNk;<Zik7rW!Z
zt(7yIlu=>vA9|6ibofHs>uotL`0?{H9PQZuxH164O4lEHxvvKTHOn?$idCcX*ShPo
zSLjjhW&Jrz`x)_#9$&i^QMg)v$C$G1g0b-0xSFQIFxb#>x#6oXzI~BX^gE(4t50VB
zc9jvoD&6Y#AT$Qaqt<!){BlME`YvkA1o6Jslm>e}liEdqPLFl}8XDzwg*8>ae*&O-
zw{G%sFMpZW>0t(0oEtu9NoCD!O(Mdtxm`M5GXatAW4EvM_LAqYa(vqYR2({g%Zn2A
zeQJuR$A^ANv}lHwcD2W}nFEVLn&@P0a893QyUU6?s&SH^Rks}B*d+n;3eL!b>j%U3
zNbP-P%KL&=N2Ad3k)m$bEen@_1#y)}{4;&^ny~-$gok1u_EcjQYY!VHfU5lPXsrQH
z_N%5px}bq%w!Pk}XL_V`XufFcuCG#VJpIO1L*X_wc>Mh>R&w6hUMGh%j>n5(>yM?L
zbH|h~4UYzOjX~oj+POpTgn&QxU%q(RXK?3{odUtC+82!hfmhy3o1|DI?0Q2TMi6L~
z(UJ@o>l7MvoilQJ_2zzXD@QBB#z*uMmuVl3ZsQ+OJs7tiH2WHxp~0g(-TG)-53sf?
z@R$arT5b&5FwhTczEh!j6o=~Upxx$#zb*1^TOW4DqonU;s@aBum6tYwaY)SVw$lGz
zARbL+4>l0P>DSN)YO%N%HC_c8?Araf{jmN<EL?ImLv>RF=YhLsWwg+U^)uE3tXs}{
z=?^Jzp6|AkS`7TqeDHYcmvntbjNH8IucW1R>X5S@)Kwl>UEhEk?^gEhcikV|PWfD&
z@lXd=Hl+|Fc;7ongKg}XpP<FOml-+TV|B=1lG4`epf|GTQp4(<D-E2t*|qI_d-RX(
zIs4(Aw^FMvw6g2PfIaSk8IDM3RimkU)jY|r0_y<1Wqy843I%Jwer`+w{Zvs^R(XN<
z=`}*omzMd@1EFAT?}1YWG@E$g_4V)mu;|)&&yRjixY(Btn%jxsef~gq(Jz@B5nUzV
z%bR8q_>mYBV0TQ9b8CkE>#fs@e%vJlydS^giKEv)9auJ9kJr`cMQC*h`(^mF;z}@F
zPRA{e7!i&!tYMrHCHm(9I*ghWTX%~?41V@_JN@uc6`Iw|t#x>>24e%ubg%j;LG*p`
zaVS18>eoIuJseK=%=>DVDP>;{b}B!-K&qyu&NkN`ON-}^&97Y+2kC!nuu~5Y>Eg4)
zRx#+z2BVsOxcPE&WvgdSSXOU*$K34qQnuCE)cw6)NgXtH^VfHNAzjLuYz&&O73(YW
zrPOKij~;J(XkanE^Q|@4?n#fk&;l)`0YA@l|0GpGe`C2)`BEY6u0tws^y^o*ev3Ef
zOv%fpzSdiKB7aEvZtD&>$h^PZkYmbCy*@JD!v$7s5DN6g^wr@*+EwvVeEs#lUnS|d
zrMH;pKIBXF{%s!nVzHeZLv3Zq9~(@ZAJxHPxEB29KT;rNIn#@o>5f;cRGUv}BSf8h
z5so2j;QSVZR_b+Imh=t6l>JjQw$<o;KiIuy(Ha$k*#Ni44f!;d75y@Z=Q2@^%!~9O
zPwU10{pf@TRCeQQ8ekbR;7h-DZsNVmu~_*baO+3ANbI5QR)0?t0%+Q8x-S-)^y!*&
zW5v0+Ly!375k6^$5=0$(7>D+@H+>sEO~ip7Wy|>nxPbSM^D-oUKC})4o35MtYvuXu
zSZi-7_U|^^l0Mx}w)T$f73cg@^4L@3^~fd;psD%Et{oaVeN~--O^^9uR_hZR54YGR
z-ML1M_0#Xt``fyV+@c(*QP#bm%O==i^8Ku@k*=>yPc3)L^PA>Kyo`ONx}TI&e?a^{
z)A*XEF@Ym2P?yHT4x7!uzwD_7gCD7>v>5T0#>dw46H#|y(<j4|_2|J0{O@S!yLYT#
zj?R<G>TXYqhI-&Sm8mgpW5L?5MVYZ^M&sS%mx1^g*S_8O)&98f;Vn%hbndVNqLB7?
zBs;+aaDI9F8|H8PvE}xEbK2^BM17~eT0AdbbUhpUJ3GV|4L5glJv8x;^uhXZjjb`6
zra|XhjdmDsDPQ|LyK0ZgS*B+(NsoHIb19kEji}p)I>)pEiKgQzSDJ*&^WOamzFR*Q
z9W2tneMr$`!R;l3!}n+v-*1zh>+fCs-W=403_EQ~#FpjXhBj)SVXALa<I2v#1=7`X
z1O%Lp#=BZ~lNWt2GIhB8x!J8x!FaJ~+LTs_VWQvJ<pr7^W|p?lWATg9r~h2{1s^k)
zCL*M<^~QxG)$$x_$4>7vBJI9?58do|aJ=KFLEir7(}?SGFbcN+yjw4!idg3hOM{i|
zR<&zk_CMO#9krM`Z+iFwgBC{qB7O0;vDxhI4K!kX-txn>V;;@ZZfL;k`d(j9H4RB6
zTxQG4^O}8)4A}m{?{P$>9Eta(GpxbJnGQ+r81tX`pF2-<I8diaf9LW#%p1_^$A(yc
zlv(d)K5c%mm|Od6O6q+m`17(Z4L0^Hchw;E*x3_ps(Q#%o1R?opOvTFj}8866+DJD
z%#klWm+YGM`h&<{QpdSLru0vR(w!m063Wc9!VNmc9~T6{f(^Q<*`^B`8X)3BQJ)XY
z?r6>cmpydO-?^yj-B5$NjjhU@n)O*aMP>Pl(HBSeT4v7v=z_pq#e<j8c`oJ>BXays
zRQ+8l5Vdj=X1*`wEq?z4A)qdn(3E5(nO0a@4YtBz&-`oO>_g!GptWV|u0FWipLTxF
zc=+{`GpEt>j<?$t(5h+(My%SPKX$SdrpIjSzdX?&L2Q7R@}!uQZ;M?X=16Tm3~5=n
zPqFkU_OWw@SB^x(iz%dcu}RLG;o;F>g{cA?Slim--Zc6K|M`Ntr*$)aIHAS1<QwH~
zT7kN?(w;r_=v9v@^7&*Cc+{%BL5EN2^b78ahBZCr_BUM7Gq7nrx8H?QdM83&j)fw5
z?CUf`{WnrqUqUv9hoLwvEoAdLFU99erROca9u<5aO7}wtY0Q~bM+Y`W+lJ{e+mkBp
z1W$aB7@%Z?Sf_4i7xSiE{7Ma0(e@e?7B2eFadD`}8rh^^aW0OwM^)DN_Y6gC)`(oB
zK0J>aux4h^vrn8X3tDAe$`yNqPzSrM7c8M)AD!8#yc>$BeeU==Mtq<0$ryfmtzT?l
zh^Pa1sBbKEZ?o%i5JKkMDZAsi9&;RRD<41r))?9*=1%p;qOvpy?heOgTWZ{j1EHT>
z!@q`85EgIh9rCh|HD*6-_`s^6Gji8;Klu5cKbCoS2wpwX27G?p8!jhN8+K*l52;0_
z&-_y^)ECRm2wL+&t}Cs~bAEVts#n*u{Z;a_<14q8J@19m43JwN3O*M<2mmX89u*V7
z;La(ZXf9*`^v-+}|NPugOc_nU<WQANUDV{!w?pAUFzk=2iRjFL#m5ogeQ0zXnlb>x
z%>xp_B#*9pLbc=m*9pViVApDbc4vMt+B?0kdcq+Hacuk1gwFSE7jN7(BVg>R5$g#C
zURuopY`kyS!5&`;&g)pt3ybT~!S&1ynXJsEs?h#>(TKexRdN#s<2LY;y*ByI+H%kn
z$#)z5SMKa7)4dn392fb<gB=T9Mnyv&VdGqw9Rz9;jC)gyq<U=Nm~6zYb?e7XbhMU#
z9v^$A-!eONOI>>5d^a<BG_{SrTd8G^qrGv#3|smMCR)In#^&;4d{BSkp5>!0H3(({
z&l;-Zml*I8i5Z5Gz88KKNP&YApN%ipU??4*wcGil90P=3c;Xp@Nv(A7U(+u9T3j@q
zx4Y}Jd!sA1H9Bb9$~+he9Po}rIBUd?#p5=Oy*{J#xI^W8LxcwV2`xJiWk6AfSI7Pi
zutS^S*-fU=oZ_@-`^1Lf(RfcVdD1O<U(1h->sn0>daLvbqjB*2%va-PX)*Frc-e)s
z!qMe=+d7v^$Kyn(tK;G7kznod<dc3V%Ld|NHF#M+6A$m24GMf3C!q5xll_!1Ptdf}
zzU6jHIjCsGu7>lRz{*YoBjGKH-bc$nJS$cBWovosRJKIlb<2O5F|y10+BbqOd*Qi(
z5FE=ST%od?FiwLfH~%(%TP+d)8o%AVvZ?`B3tg&BC~`&uwQ0eNJxcm~pOX1D9J*I-
zJMLc(hQtVnlY}3tpAL#?szuSyc83!xX;2gyS+6$rV;mQbFORDU1T>{afa5m~ewfK%
z#)YxsTu4nsAYZVNI6|xWV9Q9v(bQ2rtsLqQ8}ZOOMa*wb)c3N&vnLD#SYce}A?DtY
z090^zv~fUSg4}uOHmcM=OcoPUMi#VlMVnpQ&Ll;-gY!D=DuQkMZD#_+93B^^cz?k0
z_-6@tLTz&p;WEQm<9;I^Y<!Dp?pd<_!@!oC+oedjvqC#XRf5m^;6G}|?cP<zg(4cL
z%qrsdv8OW&t7*WxYX3fZv?e&uzMc<O?joT7p9b?T5wKVsfD4vax3rlMfZ>zxP4qtH
zk6i?()^H8R;rf?4+guIAAa?M?W=k*TQiHNk3x}K;<!hA<#<NEg23KDk1_sN&9OWZ$
z^S=fhJbb!lkDURy+T_f(P1Hxcr*^WxNdt}}z0}~y`cA7%Gm52REkm#jVx6uHLLCNV
zT`UJ5`@Shl-~2Ff9CV%*O=wx2A8(WgjLuqaJmV`<_i*LhL?D^Q{~CjgGRNZ{#GpC@
zltx5hgrgyB?Cdb{b#esSGN7V&0ysW>IHF{Zq`rj4{~oX6OYpM<oN7dq_l0D!Zl>5O
zj<+j)-uEGxw)>*!lKD9tw2EmdhakMoytBVA8%lmZAM{Rd+hb2-Z-FOjA_V`qFdQ?g
zl$sw?>W&mju%Fv!f5rJcrTGA*@lKj<4~$9(Lt@|JtF`sihcmpiO$<yfgxubb1B3OB
zR)_<Gg^Sk*iuqO(hOY#RP05JE^v-RcRo@tbp#%pcoTh%5okR2DFn0E_-uA~5;q^4-
zhDG1+66f2R8!(C;-%q1N{VeA$&P$InqL1m3An>*aiQ+t39fWN@nf)twbjImvb>@~c
z{xFI2*&Z{?*E-&<Q7~-kzIyn4v<)7Fq<$OTBM=KUZOo^P{2(psr=BFwG>ZAsHVo-Y
zld6BY?5{W;p){XIt2V67OD$fVrWarl3TMB((4&|2q955ECFbndFr53LiZ`nNNJc`R
zO1CyI!8`t#ay+SV2=SAv{eCDnBYv?qjg9@61dI9@8VO=q9KH`KkT@Pv>8B~1j-~dQ
zVd%~Y!{4R>m9-WwbiTVX;6B|T=G|@4=&;ZKNOWy~Y+hr%{-kTBsbT7&o%0Xa%l}<R
z*Mj^4(~IEmO(WgknOgR2<Z!&u0;yqDWo?s(Ced%$%A9|Au9Drm_7BgT^+uxUo0v}<
zFu1XOPo0<}gY`IZlnBrzMs(#sZy0Lv##D>Y&FeDilnsDO_IvBc&52J6qDJMTPVm{A
zgOTodW2oP!0*TjQqe}2TV~9^7HNoU9R99I&t>lSj_dEDaz2$`AYE@CS1|Lkm7Y5bY
z*)afKYycu$vgl#z;>qW~NVihi0Im@@+vFnC{YaX0e4@dQi5&)w%<EIBc=;u@W!op_
zN;AS2RE?eFji0NR{qD5QQhv=oeTL5Qin=XT2WjP+tZ~VVM~x2)B%h<ss7$nh&fiS%
z!~y?Itr<`g^~uEX=YtUlJv()4a|b5~{^+mZhwkJ_oR2lQf%ArcU7>XSP+A}H_y-9h
zx`m7Sb1f3o-Oe97Dij=lYL#t@&2zXrdrE+qKRsg6=d1cb>y5;(l5F6ZLz0-Uzgvp+
zI6MSXh98^${hbFkr_!K)Fb;2qGQ!ddY<wMJ6((>gV<=W!AJO_#a+X9z@L2EdMsRp4
zBnJFFgc=1NSjkUuzVQf^{cHdEpzZ6nF>hA-q5Ff(UXg#iaN!Cygv8^Y+|z!a|DPc6
zKH{uD&fTMhwILWx&M~UL8gq&MuQ@=4HI=FqTQxUVydEjtpK{*rjRvGjiL#$>A^5A<
zQfPcT%Q)wbKm5BaTyRjM!9$`qPSl8ympv!o>Ue~hkNgtE^VyPsF?)`58uuv<G+eCi
z{6i1V{#}A{-|E2rHR((V-c^O>28LT6*Wg=)HiMpR55Nx^10uI;;LSF^rB2L0D|KM7
zOzd|v3>f~f--!kuSWIx)y*6QR43`M<`eEX{9PttxuiGr|V&#oBW!D{P9=B2AcyB@o
z7z>?|Uh;G0MuPXnM?Ilsf}DC9G0$H2LC=dDK7Bp!gF6f!yh?wj<L}fvKO8N`Iu6n$
zm06Yhpp6DZ3I`uY5Nz!}eEfM$iz~y}K%Fab{Irz@rv_x6tr>gO#Cg-x4khb22J8P3
zK{YfUJ`Bhs9S`SqF1ibSNI^W<F>&=*l$b+e0|d`KPmM1`H$Q6ag)vJq&pcnQ7WmWr
z05G}DiCk+8V}MkaqrCiH*ZBP@mN-E0xA8jhJBwzonli=?F+~5xcXPt@<F4rC{z~Hd
zj&B+~ZORCYpAw^gzIS(!J)5VEX@1^I;BJG|C3*xE63S}q^Z&dNb7Or6oFCrl;jGsI
zm^h074Vz+#^9&P>V&0GSg-cA};5GX<N~t4X{(G8eA#>jEJDsE9?`l_Yd}wNQBCGP?
zKtG{}n3OFgvW4aDCwQe|GjLtmSYLsIzth3X=fLp?*;evWqMO_^v|!)0j)gT+I}zhM
z#GwSgr#Xx@*kiQFIzj~`&mAwyKhTzb>mfeBS_3wQ)*9nqvff;Tj$$a`a@BsA(!QF@
zTqpFx4?Dk}AlAea%e#<(V7nLSp@sfEp~j;AM~=lu>q>CTSTvl}lS+aa%3T`&_4l_N
zNyi$x(Gi$(v-el^Z>nePyjWrYgX#9E;=$KxyqlP_hiFAz+!`lvjsH@x!{wN7zgH!^
za*kRi_HvQWyC^i-{={;(Mkn^8DNf+KZ=wZ1btMK|FD8XUX+A}1y@mXBP5G!CBbXej
zPkL!=9kQ?JmevEDw>8kb8@tZ?TbsMm29oc5NKzr{<TYa2RWir<7RKQt$=WW>(c$>t
z0apfQ5-&vU?40Vv<1+qsq6OG?-PBQoi8RIi=@tcIFE2iOU_=GBf$6<7HsH)MeTi<N
zgcrX<?}_(IPDmfub#tw7t>A6^j0k1|wRV9B9A4Y@%S%rLku0L@UgGiDe%tMX1_SnQ
zusP^of=h&fP-{B+bwdASVWqe}jPqK_8iB9n1>-^5tzYF^sIlw9Rx=;c!|*=2b_n7o
zv?k=mADj<59fG6pjy5jq<Bz>xyuVJ34!|VuG1fgV7D!6J4}Si%)S&#|pqP|YKfJlV
z_k-DMKd`6z@lP-aU80s&kLX*IP8?j~j*SZmSlkqZaFT^N&#)Bp|7I09e|gaYLzbRz
zsY?E0;<#Te51H$^?iirq0K2)&^*^I>O$<)(9_WJdn<xDqH`f-$6&VoYCeB&45As$v
zm_9yH3$~rCyy$$&@akIXaxfa}X0)rmK!<xzuDl96rUI#y`j1|%!{<yAWPNqRRX<XU
zxcCY@$J$TmBL2nV)%vj8p^bcz5LmXt_v)V|j#FF@z{73RdJM6(6m@L5g`8nJaCei(
zZxiRws{J)3&V0Y5*(6)AHdXsGPjX`i%W!YFj;lvss>Y+|bYP{M3q50c2uviGymMcJ
z)tyNhnQku!R8Ct_VV_>;+qxA<$(wgZe4cp3G~@H`3A@$jB~o9N-f~HU>5qOp4x{?H
ztvMl010zAdbu*XUM%>%vF=Ujzh5Y>r5mB|=L2P@#7hM!490+ed;*p2IQ6EGJK0iGa
zZx$cRt)7!5H6_?1Ro3DD{R4!A_(4T7nqj51Aa-Z3rmq+Fs7WbS^jpgHFnjBp6bM4)
z20tR6=p%JYO9MSHjmlcW<Pg-_u(V;+BQ1Ioeg3bB3r>*y(>6U2oK~M2kHM@VtgpxN
z30MF9^|cW?xplOk7(SNp*Tm<>3%)Y3mSwsdV)5tOyuUj<UBx_7!3zVzMjvsR`$k&h
zG<RY1OzJZq7N=aN=P|ACtOWmeZY8=t7rEBz@gZkztgz+5{G*o|<(t^DamoHI@i~n2
zS-S>1U2gZ?S8@J`%T?CuKxm$E=mvT|Oy^P0fa^;MFmJ0BeMVRSR!^yr^<$}4;C@53
z0(Tu;GH>S@aCYMy{lhnZC8hg&rRU$qQvZlX=KhGdn7gvUM-`e*d*)QSgDX0Rg+6U+
zXC`pJ8fF4NyW$QPqEEYCCZ2pxmrK<fXyLSPUTWKvK%sATQvN*Ueb_G_9J_Rv7&Uj9
z<9T}n1#aO<Iu3^CjSq!Q#M_tMPpai9M~HZTkT3D|Y?X*O22c;BIg;@W9wAVg?-BS=
z6lm&HpK>C+lhomB$~$9cn!V}%I2{h&AmOc}r|2KQ(j2ycgt<>k=9E$~;P}CoU`#QR
zAagV6h^airxB6$Q$_UnD)-u<9tA0wH7yfC3H#-9_q<p?_5;_7`nd3WoKc%#v-$_a1
zFZ%qMxu%+lUiG?%+JVp2D`_5Q{L&d8sC)kmH|yqw9!y|l495Ug5L9~jKaM!_X1*T9
z_A$-wYdq0<@e|FW@><w4psJS++w;vfl^R$mks5vciac+@Pg{ASZ!?-CHpJuRSrYI>
zddaWHG2wtM=9YJBxoO5ZQ{Opu7td^dWNMkTdbLwOJqSF_btRtRgYLBZ!53PrIh+xb
z<gLL;Qk{<cP)6*#v*CF9frKB{K_z;j7{Q|#Xu;_9vaiFybrC*Bp-1c(FLd`7I@qv=
z;-pcWr^F{QnP@GwIYJpxv?3Brt#Ph<5SouzI_qCQcW`~^eNUOxBIFbZ7yU?@0neAM
z*?0b)M&P#!>Wjl1+qeJH=%eJ$iFpkUejZE=!Ee(zvLW337=n~Ldn#XbH3M~-a;9y*
z#Py-$2v_qVCBe2h!Q<aIqIM;Qum#~4d+l~Vr85_1`qiycWRkBB|KV=?R^ohKWgE1L
z@n4;N@0}^RGcoZ89OVSfs8ji!9fbZ%Wsc63uN#|Y+?Tixd5%ipNJHrTGJe0RM(Bkf
z5pGEB+}0cyndARm_3-)PzwBAqc~k7h7X5lA=SZYZO}>#I4iei>^RiFEnScso``_~e
zTZNMrX@sw#Z#?exs`;#QEnhf}Vv1K!S?o{xy&qoan*<xV=k}9<uX8+P*3K0d_=!IL
zLMWKb>~3G9=+mhm;ChL>$%1Dk9ZybaLMktjkAUH)b7D*M0s*3bIUk2JZw7p5zDSGu
z`-;n!&dQSJKRth?uwMZF=&N;UQ7upUw}UoplU%6`$>+}?SMp;543!PI9(9<tzzc5K
z$ev~Q?yu9t9MqiAAWF~ud~mhZ$L^VWNR6saO*2zmpQyASM`?UqX}u#q51Z)(E;Bzw
z=)n#+iTP_rxGd_w9^7Y=WL2{MtxI_C6w!}_RTA@;MWo{MxzhY9pMQqNiaxk<d9X5-
z(<>HDW@ndPZ$|5t>tmwO*|CyK?G`pN=TR^FVK2#s*PMz+_Ao}E`G9So@qgn)9k+7<
zxt{1-6$qaCP%_x@(WA0a?y`6IxRXnrWaCha6?UWI#NQX;D~B<{m2@WDe{>^BX0p&Z
z-IGN<*%1bYHjfRr7y9SC7`fBG>C4pH17yx`x<nyk+R68yFJ+h*{dpm}6xI%WjhZ%8
zKiEcn8V*u>mhxK{rMTZ}X#kxRJ5@ptu-OU0k^lNm|M@}UxJjBCSv1bx|Eve0)jhYD
zcNBO8^#Q+`ZETo7HjE>L<hO^+;CuUKK~S2XQJR0_ev}9C0w*UOBiB6-Ndlt>es$L4
zFsb;ve2yxabBe^=^VLSYPvV1^?63Wj7@^lc8KSuUPs0SVG|!$Pd)bl3Aux=K`|<3A
zSqU#>5Ij?Kg3!J6(P7Q!`72+CxM50rn&=+r1g@MD4-&tccg=Q3NU30g;POo55?mlX
zc9!?tBg4*|K$%={HHuz2$k_{vue?0<*~D?TmcJyf^X>UtdPVZwfV1}E@7?Mlmu*p;
z85e8=rPs?6E9%f@1&jHpx)IBI<~MJa=LRw(?CRIlAB1|GRp0Ue8@InD1>x8Z0<Mm_
z3*N;`4{b}?s$ZdA=;f0wg)W2eQ*t8?oay2NE_b;?`WB{FdtxbYmF}+K{=_jle5dg$
zZHAw~=_jkO;3#X9gK(rf{bb=5Aa@)I;QmQj1ed++r9$4(QcY&N`<LL7_5ydjsKO9V
zWLwB?BRZJ9uJTqQRmd}^;Vu}tqe{Q2<lo_VLBBvKJs)wMx`KEj@1v$i3E#@1P@z*I
z9rBKvgv`*q!twLU20XYrx&3Ebe+b>K0kO=s@;*-RO^bElyxJH~)c@3>aP2t@`N?ns
zAcuP5U#Y?;^M?2f{K_U$W~xyq((?=7WunaFDm~tN$lPE0B3hj52|+TqhZ3Hyv>!q$
zJM{g*Hg2*bjjI-UA*jX(pyfunE)x<}vj_KiRQSI-8Y**qdx{avepGZWOw7^!|C%_z
z=4DV^@2EI`uCU)g_>c>w8f1Slj3_1evW!rWSnjX+WgQlfygmJ!Bc`?@hsfgJQdP2_
z$jifI&XbUCi0d&Qkq?OZEpLP=uFsjv7_1*2V0H(e-w`k{qGTocs##y1t&@c=)e4_n
zI(2slER?t|%<Y~e{6jjSms>%2`-EmBWmXCN=x`V!<GKvp+1(e1_H7L;TA_i475E}=
zfxmT26n%8xSTGye5^afYa8bCJuXZIM;~5P|RmuO_%f;>b0!QJyneTy3m%jJS8*Gau
zq)yGa_QS-_-Is49OE#d+4TSTIQBh;Ixq{gPW7S%mUtDv;(mTFbbN5W1RKr>LvV*+9
z=#URVS}@yE{lgAm@VZRrB5PloXIX>UWI9m)?!Izf?<yBu!N%H)ivqDfmku~5wa^uf
zDUg)bvnqWaFnYblB@dzRGK&>>^&$_U)94ZkX18c!$T4w!a9uZiO&z{wMCyB~A{)4(
z4e(_~fHMZ+iw*Y!xAARqMq`o#PUz@}8Fxv@nCOFN6E`0pzx0%e)Xu82=SbIi`e>8Q
z(H~6Q@AoPI%x+sQI}lv&)X^OgjXL(xURDVo+C?486(qOZ*v=f>FFe2xD@j$I^Un_Y
zjNGxk3S6=BN@-dWPUubJ-KTT864!^mju1F-Of+~bfU_5tHH=C>^wa|u45%y$6#n+{
z*0|bdey!>^NiW^J;<fA!=4ieAN1eja4lqSdEZ%)yBXb|?S4%Kke{>lSdF#PNXBu|W
zD6a1uKW511tJB<J$r{%^MkG||IK$S%A1!CqH;<dEMK5OT&#{*Iocmu4Slq(%U>`7B
z<FK)%1%EN!AoTt|VNlxdth9e{^WQmpE6{w-`?i=w;g|WQ6Z^WMFU$fzTyxx4Y~u5G
zJ#(4+BZ7=DlI*w}okPqm_I#Qhn4Q?c)kkst8pk;<{*%ZZH+$ULd<oAU?7aN(saXHy
z6X)l_ZgVu_M#vvhNAGH`qn$ELOy+NVY$=obc+}I5#U}1=Yi*`@Kg9J>UDe>ai$7%`
z`c8lOmlFwz{&~op7nJ{+grC+$Zp4hqYHu(b&a|$ipJu+o&03Ip>gsy>GtV^`Zvc~%
z8XuHGc;(N-{-oWLh}GY>wNNL&BRj<VV!h}uDn`nTUOm*q2=1T#lPr9C=aRwXV7cR?
z6wjmJb*4g)(C-Y^gUPd-!AJP*9ur=7a(m4w@p_?8s~9Z!PJJMHCL6x(UmSp2B-`EG
z+JIKEQwMEmn=RTQ1@U&J@xKRT-|2HS6x?1~tFq|l2F2pyc`{7B(%=q5-d6nL^5@ay
z_u)AP;|#+0SwE}<A59SF`S}RN&%bPnMd<BaGUt{}#O{8-SLM9$k(q6*%vN7{!LLQV
z*4}iH)5r!rq`*lQx`Uq*_jS`8FnxTLZH*cSVhPc`nXLna?pVP*85WDM?w?Wk-1{ft
zyFUe*TyX?DRp#~clvz8K;+_J%8!?}xL+AGX*WqG5_liNm#Ybye-X%RQYXG2!J~}O3
zf6^b^S9Db^_qA}kwzmI$Q|C_E1N%<*l8-k3u=}`6x{2p<n1>2oWG$`yg~6<KA_R`{
zNsYwUOY{Cq)G3}XpmaZ^wEkIX{k`xv$AQ1cgQZ2ET$OZ1%y2fpjPSd6HiF56!^(RL
zKh!?*VKcos%^M?`U?|^Fu0$~}%T6RH?w91g$={{HViM|)131p#GezjWKcom<V;R!b
zl3eLzl~6GkhwH#>`jgC(z~|uf0Q~EHpn8Q(8Yn$qa(z@!PsA`Hi|_)Dn|Y~4J*M0y
z{_}>rBXt}h)x@`OpTk&laqi0gLT^HTqCt#k3^qgX(@PFx{ok2*t`PZNteGKycd_Zf
zqZ87LFJ3b5+pWyy9@h<@t%Kj0y3X44{NvkqrXK$VJzJfo!QtF-4;JX(n|?eSb@W{s
zN5%2`xqB%fgU-z@)Yfh~TWBIwrhPlv8~nVQ7B2X@rm>huaYe`HnF(BfsRqP4lJq%R
z<ocKy#Pd2Ft2o}t=OzXFU8VKAO6!%C&c9N6{s?A3+FTFeKa16r=vgDh{87dlFItjN
zjdTIrCz%qC{*RA@9eho=+xQN1gVuV3#gKj25+HI;4tOYzUn#A3RGQycnjhqNW!r2M
z=V1o<$Xu7TIS5Q%opH<{b6pVSYm98Xz@hUs8=1?=-+Rj3hMQy|^Z8C=CAVm{V)4Wn
zf5rW3O7nS2=Y#N^H3jpVxbIOxen#nei$q)oefNBc<5t&7i+q?$VJO!uyu!UYH>I-7
z7`rzdo-EL1vQcL6&_~kmGuyzt>3UFlR(<n*UIMSQG$`I5D1AQj^Hjn75~cG1l-8ee
z-<_it&zY^aag@+0)C@&+ieEal%u?j|4E!tdF@_jG%ADL^8~jWB@CKwXLC7rwE{|e}
zyf>8ABPpHdHdku+WZ#ynrh7SKYj!z5Q~F!s^WNr<Uz;ctpyHkzmNJWZnSQuPVm8#A
zK6#S44>^R3^Q9|`m=W7TzS)>Ldvm0kTlb7!(!D_F_GsScb5a=ztWDiZJ}}N#tn(B8
z%9HUzr+LOtHtyf_Y|PwmChjv$i$f;Gllk_uD(PRtME*pc54evuS|fO-Bgw)ieLh6y
zc>K>O5KFaU@mXJaeKrLsq^SizF-a}_lkMVUjx%kHgVK1|cS}-U%y7lJbEE$=Z>B{>
z7ME5}h5yK2?H=F-MoasTjS+b8v?~&YB~A>uua@l@u<aOulLv=oOxa>ra*o6*&JQZB
zZ&JFRE6t~bkc{9)gbEJhe>(Xc4-&dCYlGk$b3-9=E`CW*hCTID&3<Kan9%0UWy>=X
ze;?Okgb#fm@wzPcB05><naV{9AMEX9flJpa@mEEpl<?~Y)M5jVmr?kvUBZMeY(}K;
zcWkl*_vM>u6wjOExh}D`D93=d(?MWUXs>#aV6l<&eiC2J{8A@EOZ+T0;Qi;f9257=
z*NqlAg<GS*b%x0S*u`Qr#}Pivj3^a!BIn_{g~;7nWG(k#0<WeCSTX`n7Jt#evP1bv
zQ*V-QhWp*jgub+Ik;HMDidw8)M1eA1ewdxn!>oUN0=Q4DF6p6}5p9GPTwZuR5J597
z#_0N&__8L2;myLgsh=oklv4IK=<VY$Mn?$)|9LCk?=I9);+nH1u9F?Qa^F(mdh%SY
z@Ktqkl9^3q>j!^X<Rtu)Sd86~_9~Iv^aH{dV^#(|S<uq^WX1J$oX@>q!rKvF&2fb~
zT9MOHm2%WrW6@rRYD~CR%R}hN+v-8CkqRAsH5kT%j!x-Bj_ig4$+-_DYfjgq2NxcY
z4mP{+{)p`%V05u{qte1x5FCa<Gv~hhwKfETm(U45FpGTiw!^z;{ZWJWGuwOxPp2S$
zIy;&%D}KVyYiC3Ziwk{djd9a9U~QNWcy4hyoxm^ego4E`c@qD_j@{@J-X%U`oydD2
zAE(mwP3e3yp0i->0cM-3;7qy1o2EPM^a&C=ic@N=ETTZ$kHp6gEpxo~@GBD6w^APD
z{M5F&QLZ-VJ#kyUSr-dbW5FrMT%h!N=002n^^;s5s^EF5v>&JLg^MTWmU$}i+_iR;
z3wmoTy%BmZZ2eunVd~a<%!1oq&Vt`_vA~s<kMFri_A=KeH`jyNJN5@y%S9C5QE>si
zM}_kamChHPPVSNY2eVBqCh|_9pS+J^sls-LqAnYV-<y@>h5s<|9Kk2<qCS{ZB|2+Q
zY-RD!8@$0}HU9g8!1K$Og^S#gj)^?RdR4ZgI8RJLGSAyw=7Vb_r^;UKDf9U(#9e0T
zM629G#MkXF4W6GiFG}cZKk0<;wnn7LwX7A11LN(V9iOEEi(%bh8H#$W@gu)QyE>E@
z<R6ax4d}nLi5B^#wY|~xQguzcbhYC7?L5D8yb;_FLq67~KCSHr#RY)H@XSw*!&;K-
zt=M8IbKgaDBDnm3cvz8JlOX(4R}C`P>5)HE>GMzN_seLnoc)ydH6!iL<?in!p<@mM
zmp}XEn|MB9JJKC72C1sI@JV-!0H2G$hhQnihNe1d1n)%8OKH4L>3rbT6yJU0lAqv>
z8~7le;swX#*dpxv+a7bHUYR(KyV?_6UrPQdW@ohu_7b>nvJsTNdt-p58|su_*XQa3
zN5%PbrR$N>{hHEwT1ubyO7p89*S_jhEz2MLeu>GHlN@gSydCA(>dr;Hxctmh{+vcZ
z)O{%5mgPR2wUC+Jq*Jl0n3w-W;j=xhlwI~R&&wn|*5;=abndPd{pQiKU^b7`KEBw}
zkp#;Nt&kGAAtZ0~W7ED5W@iH0dI%lIlPIjFxY9F!6LHDC(yD3hA5A>ppo>v)J(AM;
zd!Cb3BM7FShr2K7qQ#=-6yR2<67vJ)DDpfzZLpYgvO@5s*{xdN%LU>Dvzcdl3BSo&
zEk@C;4i9F!DxQbT^QCraz{<?Zq2zC5xoX>0BJbQtdKv~q<k0)#Ji8Gr_Ux9c7L@L~
z=Chez=qV~|F_+wQjmO&y-;08J5aIhYlNYJ4zj1AF*fe-`b-FslQy%1Vy}fqu1<A;Q
z<6CQFH8aMh`O3n_StK#l&-b(d*_s)&3Uei%Z)$J{rO1Ey*ijFjC!FpBejm%JlKikn
zLD;xgbECx}ACZfur~Zctu=<98#X6=hh!r`NwJ86E1#XhxRcZYYUnf=$xI;0tV{1nv
zr<axE@0h+4_x;Ik=$CgeIQ%T-Au=Irm>oRW;OgNj=7)x%!j~MPh84wky&mHwa<>%3
z|13Y>*p?gd$|Q2Vy@fBD^3nJ`D%O7_k*n$|^LuFe1tExC3}f;Xxw|VuU{;o#nzBmf
z`E^Er#r<VU=UFQ4ci}#)6CvU}OQ#&BLw@SNWxT-c-3h@ZzWdUmzZ~Z$^ZPVfCW?Hn
zL6HJSd>be7sFDnd>n&JZ@Ro_5C45|XNq!pnnAv?Th!s3=XcRWu(8U6z!{&TRhL`X!
znkRwDbi>+3!;F1&2vwX9;5k7nt;9VVRf2G;ECoj-yMg;!>V=_*;`y7!t8i~E1@-3!
zU@66rM2zyqWs<!f47SE|qRo#UbQ8D^`RP`go%@xU@k!!&h#kYRVPx+PJ9D%$m)lcb
z<cK3R`X8P5*u;4;*-P$4vj4d(&vJX@+V@W=&-N%81}6H0`@1!vQ2OUn+Fz@5f6w);
z6Cy<Z^1)Kz_aS)2g4x&A3yoq<Se<~&mB^`TqWm)Er=Oq$zn^14gg9qAr-1t{e`^)L
z|AWIqlrzEer&JoK8l@lHTmFmbYoE(f;`Aq`P3QMyWi5Ga;y8D-D`@rntR7-kqT2{B
z!S~2l&hsB!gr8hNe1pYuTfFrZbK=1G5?^eV#Pzt*X0qVllHjo68zJ*QB%W7xz*6Mj
zI}rX)uy9_c1Nc1KG7zgMrlCqr9q5*$$gclU4nMn(<ewX;YrE?+Cm1X?%<pYr$^8mO
z;ji%U0EL02FFoXkr(~z^OXqT<v6O(WQj6T3DU@f-^7wDl&tbZMYbE6Za{SQ6Oy>7I
ztc@0V&#&#kw#Pfl;pBSguN6>(2W!}r<WU+ijs>#VmBvGg4Y)ha9y~{3Mud2-#*~sA
zkZ6(XOZAprW0UnPROa_{IaNkwx~?O0MHP{c@GDvIyfoI<mj9J6@p~J#gyP(ZHWbxU
zqNBX6$L-rBD6DN);v>(Kcpi0BoGk7WaK?>|QFFFiACk`1Amw~z^1-wCQ3dt<O7}xb
z-%m1?;$gNKewu_2^Q(!C*Xw6_$b%Z3t(zfv%ij0XuYXD)zXtO!Z*nfldsE94?z!VY
ze+!viv-P*Jleo{In+ufoZ*hMS<p^>-!=)UkYc-o~>jQRs*pMtwnb(2!Q6i@{$x-CQ
z|L~DHpSPLv)R^wKX*^iW?zQ@H;Qk59+hDP6)$Xamby^R+v7Ta_2Orjn``u`N@xGFB
zm^eSA3NPU?;zSP3RUJ4T`OsPM`va8bw|K5cf33K$z&QadrtHb{AmOuF8VM(gy&iA>
z$i(wke_P62k9Ro%EWX-;?$hADzz(v+{ZJz#1>Zc;O5FQaRflJF{5B3<>?(3<t;1zz
zt6dgui&iYSwt`XSId@gUWUik(5mu60O7}c5KVDG0$epkb1&dq3y$I2-WR~PMpQW6!
zfs~j;yb`m|E_De7lkIgmsxQeah!MQ^w+NB*Lb^Vs{Sf@#$2yb;&gZ~jnZ=AQTId5t
zSE!rBgXdLs2@&~fE$n4}-_(Z!>8v*eo+Ty=UanNM@XbvyD6V(pzBmQ_JxcR=O7kC_
zk2(~E(X7C)2ok&z-EYEi8Jj4<_Z#BCVvVx4>rjwC+tF67FY$9z9b>+3${*u>Z^bah
z{R>Lxl`8E|-m``X(Tc8^MAzyy97cJ$L@!?K?Js<G3n-Uw%_!HUZJY$Jwwm-Nl-}NR
z0^NhfE^-{D6Za3vq&M2JLi3sKnNj-wOMY*sg7?p)(lvw8UHnkfpE?`L7h>aWd#$>}
zUl)#JWZ%smZjKI|xYmQ;5A>!`qV!7dNV?B(Wrf=gPPH_U)4o2o>!OCZuTmv)XO6pr
z>rVda#d;d+F7ur6vGHKHSlqATB+j{S@uI)#6D{ft<)L(B*9$%gR6PG)>3k%m{dp{=
zy~0cbHkP9smPr5oiDH?xZ4Kc0JoZ`)>tpq1Xax=0vB1g*Emlx${^nZ2LJu(}1QS9i
zF_P}<AW_gikMV_)(*6ad^)yP?bEWZiuA^$NmecRhg%6*+WPYDvQK-zXp}Y|w`u814
z=y1b)V#;xUndg6=aFzMJGevPYL@~Edua+;#e~5<Ce9yE7G;w**eVVuEg1Nn6=*TXB
zvHv7-eN$a6#*)3TJK=hK9^CCG@`@K0nMB@w9Jm~1@Gptify#PuKg}H#_`P5Qox!fX
zdl_UX!2=B7*JE}lHu>Crt)0F0w77q3hovm;KP}NI958~{_jBYYRrtOHpSJ7<=+7p)
zt>IGmNLzWua{BQ7Vd}yHxwg8>eaJRZb)z{_58^}>hi@i+A4O|Balh$CTMQxiBtOPO
z&gwE?@yW8qlG67lE4`k@{nk#RpIV$O?u%}x!lfQ`Vb3~k$$4rJe$mrD;(ku@fssr4
l&f7=6VAq3un4}W-DEiPm%&z5DFwaA2{jSpYohY40@jrJ<MF;=@

literal 0
HcmV?d00001

diff --git a/test/test_resample.py b/test/test_resample.py
new file mode 100644
index 0000000..d9a5893
--- /dev/null
+++ b/test/test_resample.py
@@ -0,0 +1,108 @@
+from pathlib import Path
+import os
+import pytest
+import cmdstanpy
+import numpy as np
+import json
+
+from retrospectr.resampling import resample
+
+TEST_MODELS_PATH = os.path.join(Path(__file__).parent, 'test_models')
+
+
+@pytest.fixture
+def eight_schools_model_file():
+    return os.path.join(
+        TEST_MODELS_PATH, 'eight_schools', 'eight_schools.stan'
+    )
+
+
+@pytest.fixture
+def eight_schools_data_file():
+    return os.path.join(
+        TEST_MODELS_PATH, 'eight_schools', 'eight_schools.data.json'
+    )
+
+
+@pytest.fixture
+def eight_schools_data_json(eight_schools_data_file):
+    with open(eight_schools_data_file) as f:
+        json_data = f.read()
+    return json_data
+
+
+@pytest.fixture
+def eight_schools_data_dict(eight_schools_data_file):
+    with open(eight_schools_data_file) as f:
+        json_dict = json.load(f)
+    return json_dict
+
+
+@pytest.fixture
+def eight_schools_samples():
+    return np.load(os.path.join(
+        TEST_MODELS_PATH, 'eight_schools', 'eight_schools_samples.npy'
+    ))
+
+
+@pytest.fixture
+def eight_schools_log_weights():
+    return np.load(os.path.join(
+        TEST_MODELS_PATH, 'eight_schools', 'eight_schools_log_weights.npy'
+    ))
+
+
+@pytest.fixture
+def eight_schools_resampled_samples():
+    return np.load(os.path.join(
+        TEST_MODELS_PATH, 'eight_schools', 'eight_schools_resampled_samples.npy'
+    ))
+
+
+@pytest.fixture
+def seven_schools_data_file():
+    return os.path.join(
+        TEST_MODELS_PATH, 'eight_schools', 'seven_schools.data.json'
+    )
+
+
+@pytest.fixture
+def seven_schools_samples():
+    return np.load(os.path.join(
+        TEST_MODELS_PATH, 'eight_schools', 'seven_schools_samples.npy'
+    ))
+
+
+@pytest.fixture
+def eight_schools_bad_data_file():
+    return os.path.join(
+        TEST_MODELS_PATH, 'eight_schools', 'eight_schools.bad_data.json'
+    )
+
+
+@pytest.fixture
+def eight_schools_cmdstanpy_fit(eight_schools_model_file, eight_schools_data_file):
+    model = cmdstanpy.CmdStanModel(stan_file=eight_schools_model_file)
+    fit = model.sample(data=eight_schools_data_file, chains=2, iter_sampling=200, iter_warmup=200, seed=0)
+    return fit
+
+
+class TestResampled:
+    def test_good_sample_array(self, eight_schools_samples, eight_schools_log_weights, eight_schools_resampled_samples):
+        resampled_samples = resample(eight_schools_samples, eight_schools_log_weights, seed=0)
+        np.testing.assert_equal(resampled_samples, eight_schools_resampled_samples)
+
+    def test_good_cmdstanpy_fit(self, eight_schools_cmdstanpy_fit, eight_schools_log_weights, eight_schools_resampled_samples):
+        resampled_samples = resample(eight_schools_cmdstanpy_fit, eight_schools_log_weights, seed=0)
+        np.testing.assert_equal(resampled_samples, eight_schools_resampled_samples)
+
+    def test_bad_weights_wrong_len(self, eight_schools_samples, eight_schools_log_weights):
+        tmp_log_weights = eight_schools_log_weights[0:(len(eight_schools_log_weights)-1)]
+        with np.testing.assert_raises(ValueError):
+            resample(eight_schools_samples, tmp_log_weights, seed=0)
+
+    def test_bad_weights_not_sum_to_one(self, eight_schools_samples, eight_schools_log_weights):
+        tmp_log_weights = eight_schools_log_weights
+        tmp_log_weights[0] = 1.0
+        with np.testing.assert_raises(ValueError):
+            resample(eight_schools_samples, tmp_log_weights, seed=0)

From f5398eb143b33c52d799e3f502b2de18cf83d53d Mon Sep 17 00:00:00 2001
From: Alex Phillips <A.M.Phillips@liverpool.ac.uk>
Date: Thu, 17 Aug 2023 12:01:34 +0100
Subject: [PATCH 2/2] feat: add resample method Add resample method to resample
 from log_weights

---
 retrospectr/__init__.py   |  2 +-
 retrospectr/resampling.py | 26 ++++++++++++++++++++++++++
 2 files changed, 27 insertions(+), 1 deletion(-)
 create mode 100644 retrospectr/resampling.py

diff --git a/retrospectr/__init__.py b/retrospectr/__init__.py
index 863d295..b53fce8 100644
--- a/retrospectr/__init__.py
+++ b/retrospectr/__init__.py
@@ -1 +1 @@
-all = ["importance_weights"]
+all = ["importance_weights", "resampling"]
diff --git a/retrospectr/resampling.py b/retrospectr/resampling.py
new file mode 100644
index 0000000..917948d
--- /dev/null
+++ b/retrospectr/resampling.py
@@ -0,0 +1,26 @@
+import numpy as np
+import cmdstanpy
+from retrospectr.importance_weights import extract_samples
+
+
+def resample(samples, log_weights, seed=0):
+
+    if isinstance(samples, cmdstanpy.CmdStanMCMC):
+        samples = extract_samples(samples)
+
+    rng = np.random.default_rng(seed=seed)
+    niters = log_weights.shape[0]
+    nchains = log_weights.shape[1]
+    nparams = samples.shape[2]
+
+    nsamples = niters*nchains
+    flat_log_weights = log_weights.reshape((nsamples))
+
+    resampled_iterations = rng.choice(
+        nsamples,
+        size=nsamples,
+        p=np.exp(flat_log_weights))
+
+    flat_samples = samples.reshape(nsamples, 1, nparams)
+    resampled_samples = flat_samples[resampled_iterations, :]
+    return resampled_samples