From cf5370ae0362688888bcf161e8a9146bc3a48c4f Mon Sep 17 00:00:00 2001 From: im359 Date: Fri, 20 Mar 2026 13:46:51 -0400 Subject: [PATCH] feat: add maximin analysis script, update core libraries to use native len(), and generate final comparison strip plots --- analysis/sf_c_44_maximin_plot.pdf | Bin 0 -> 18724 bytes analysis/sf_c_44_plot.pdf | Bin 0 -> 19366 bytes analysis/sf_c_44_stats.txt | 6 + analysis/sf_e_110_maximin_plot.pdf | Bin 0 -> 37440 bytes analysis/sf_e_110_plot.pdf | Bin 0 -> 38128 bytes analysis/sf_e_110_stats.txt | 6 + analyze.py | 168 +++++++++ analyze_maximin.py | 157 ++++++++ .../committee_generation/leximin.py | 6 +- .../committee_generation/maximin.py | 2 +- .../committee_generation/newleximin.py | 268 ++++++++++++++ .../committee_generation/og_leximin.py | 334 ++++++++++++++++++ .../committee_generation/op1_leximin.py | 250 +++++++++++++ .../committee_generation/op2_leximin.py | 273 ++++++++++++++ .../committee_generation/op4_leximin.py | 213 +++++++++++ 15 files changed, 1679 insertions(+), 4 deletions(-) create mode 100644 analysis/sf_c_44_maximin_plot.pdf create mode 100644 analysis/sf_c_44_plot.pdf create mode 100644 analysis/sf_c_44_stats.txt create mode 100644 analysis/sf_e_110_maximin_plot.pdf create mode 100644 analysis/sf_e_110_plot.pdf create mode 100644 analysis/sf_e_110_stats.txt create mode 100644 analyze.py create mode 100644 analyze_maximin.py create mode 100644 src/sortition_algorithms/committee_generation/newleximin.py create mode 100644 src/sortition_algorithms/committee_generation/og_leximin.py create mode 100644 src/sortition_algorithms/committee_generation/op1_leximin.py create mode 100644 src/sortition_algorithms/committee_generation/op2_leximin.py create mode 100644 src/sortition_algorithms/committee_generation/op4_leximin.py diff --git a/analysis/sf_c_44_maximin_plot.pdf b/analysis/sf_c_44_maximin_plot.pdf new file mode 100644 index 0000000000000000000000000000000000000000..8d905c834da4516641a57200193508795a6aafa1 GIT binary patch literal 18724 zcmbun2RxVG`#+wQO?Dz)cJ^!Y%E%tsBQs^M?AMD>%E;a$D>9Om8QF;>vydX8gp3k0 zLW$pfQ}1`f=lkyQ|Gzwr``p)kopY{po$H+IdCs}d&7-ZLBn%TllJS&`LGCn=!5~nG ztED5Eq$EUC-^a}cA}Vizvv6{?gNSNd*x6tqaDYJ-Z*7FoCVGXg4`!)Ti|du9xf0RcuOX#XK!Kc;9>_se}C%XYNcm` zgO~!vDkuS>*m&b0qN>h-2y)-A^53tj5VN1ef&EzlP;O9q4;w(beep$gZLqGM9#%ko zK>304HEgUMEaY6hffAv>9~ur55rd1v;3xl=O#jxGn$OW;!+4YR@XsR(Mxw+}({v?b8-?5uAUr(|TBomri1 z;jN6T85tb(#1v?*_l^`W?KYo$;QeVkZEMxAXXNeZR%81O_RX`$?5<3#zV=?t=zBN0 z`PrAhb@W&|bz1rGXgLS{2h%%g78ln(X7#Lm^Z1;k_}YD7)uyX`=9^;Ws|?O6qmcq? zrSjTat%{@FS2imidf~?A>~0BVVG+&rEI~e(EOu;Pr&bkMEOgYo4C)I;ZBAn6y1tsM zaGpO(iL!BheQq*{eX}X+yhTiz_ovG~6TP)5J$>G1w?AlvJzYsT=cg%rES9LN@tFs| zy3|$Yr_F9NgIPmc5?Se`Qp8J_)!G*VtN1Ok_EO64@RuSrgx!KB5sN+OjmZmz_xP)< z7HU$SrE3lFRDDRL*Nq^Gz0l~8?6*ccqKRlUa&j&bOUV>`wf1PnsZ~3z$cXc;pVNuTtyk-)24PHo2EKHdKcX~k%vy%dW))e)u- zt*Zj~r}c|o&xg;c<6n%PqgjQ!-`J>`N-`pU@IJ{=kWpeOFunP1qBS|`#l zL5#`ikBMgzN(Vd>5L&FJPE7P(aXA=aiCjb1t6iRWw5?xn-ZIXhiIffy_zOD+iH90O)GI{ zxHs0AT=U5rGkZ5ZQrP$Lx49!Ayvo{65dKA%}Jg+lEa}( z8#Ct^qk7ClPKZ6yKKjd0TA6B|Aitc6ZpkW;cK3Qz^sy3Out;*Z2xqq|p^; zuR*n8#bhMyG!-P9x+3w}((n-XhQ7fd7wri>zf9;SNOk47cW|&~%rM?&FqS59m;wT5RmjQvA;iX@Aj#M7vA0yF%RyXh()Q7N`uOyIhJ z4$JgDEPUPXG~zYy`4hWwc$mUFHITcFYB9GEd^U}y7RMsgvc{E<<&y=H+;gE~pg-Rb zr90?p=)@g!0ZN!TEcw;6hr{PplQDtJ2G=Gj$z*tb$Z5eBtdbkh!L}o0N92t)szvB| z>lm%oSCr9_3GOk-(}|gou!a+G4wCEW_=LiH>YNXIZ^_FAOyG&}(`M6IzRy)pd^sN% zj=V~j*Ke02C97L{D?cvrdNav#aUn*La}`f3M!!9L!fGQ;8&_v~(dx57x)NvE$`|67 z9k-CP6`rZG_$Ld(J}3FcCZ4@rcl+#s!x~;OOV^OODbp4Iz8PnxZv6aEOWuz$NfsSW=WZN6yjX?(GgnnvdF!|WS?9!!Z-Q&j zoBCs$>DKsJc`v!!;D2$Fh-eK;6R+fb`SJ+4B2|d(R%qNK{^xJt4jy^uFOJ};71*B0 z=H8hzIL3olbo?WiiE)8qk|5sgMmj&cyt<7^#kHb`#bGL5IB9(b9rrP^OF^v4Mf~xZ z>vdenCiY;R?3*0?XPA6g;$!hf!u>~d=r}rFy@Q-yf*|rAM?L7*Q_MS*UPXrYzw7`0G(>0q`BO(_gbPPVNNtmn{#-?^o93X-@WXu&lQbTebVJl z|Lo*PmQt1IYZ>|+Ed0!JZa+9?7Fvr3OX(SSj{9UA&%>#uDlaY8@$M-i+$`t9Vo0r{VqwIZih4UvuwkZ$cl0hugmKAAZ6H*#IoBet54qhUKP5BOghW+ z{IIsd$qPwx`m2E>t*wC%CpahVwfyYQEUC2L3;6Oj*dVfcNFPaPIk?_<;>84u2DA3D zS)=GEkCRvFkZ4Yxz-RP%1sl0rXEGKeNXYmpwC&jY4o!+L7kyJ{zFdFVfo=N`e>Xv% zxwP1dJ&6y3jS5euFPUGm2a8>?7Pb+jP`c}S^OCbtZQtGG^2Y{GZ{~lqs_OY@){AAX zzbb1v@4}WR+~@ZA>kU4Q#`}|;YF4pHt(@!dIHEdBpU)ncujE&~X5MW~6}{;(>=~Zn zdN~ok)#$UewKCH4W>_n`JAp8-YpCSOoYgnnQsBie{GCN&W0{1U)SvXe1Z~rkca}|6 zT(%mrU8}s(aIsb58$Gs`V>&FGFoAkmZ!2g!$+7Lwyi|>?vl9DOsW z!$Wm=1lj5mn-(9ZIUVTK4QLduy8N1g&*VsElc?|25C+}GMA~JR6 z;T7y#-mNy?EwK^Tm;f@Tz_pZ>s?a-4SFm5JN(%2ccyCntYw0Yidf|7fWFrTF&%J0%H)zq)64a zlK#|o@yU3MQjf*kC(DkV>z|e#*FSsR3DRfJ>KQfbOa6Lp^rYXt_ca2gRVj_|IVJGS>5@jl4^fUHOKqE**YD<;7>$k)%^)d|me4&UpF8`HLsdQP$V*w&t^ytnze+ z_qFbLE?=&*k-Vv@jQ!@j`=t0>;{(&B2O~w2TptE9k9r;*e|Z&I*WhU8`Qc3Y*46QO zv@Y=tH0;de1KauH2$y^P28YnK_u)17c4)Kk!;=U_8!?nYP56(3OdR~>hWF;NnR(xzK+YFmM?uR!NQJz*DP};7b z6@$~FJQo`Kah$?xF+R!fLFM+-L$-IVOWc(y%3oiOx?lEgNBqK4{4vwEMAqatXii0z zQ)oYlog2fl1!4b9eo6d^U1^;p&XEVvlr%Y&SvbB- z^H4GucO!~w3D(lh=X$Arb8NH;+SqXgBLvT@YkB9t^g!EEjXG<}#SF_lnibH@}XV-)#}0dz_F{QJ!!ZgdaE8$1w-7z_0W#TNzohjR$AAD z&%VPQ)6`>*jcl|v%kHs&54tv`*SF|f#J^a`{c^M^XSKvOwEPSLc14=h*|I|);~ zkCwgIz5d*{Y!yngv`X+sR@(Z2D*rsw|86Y^}394Ox;85A^2nZyLYu@yY1Q z+r;jo@g*39D&+8=;S;BMhjSF`lx^~8LC^6mU%W3Foz>A8XR@rH$6kkap894~#Y>i~ zBosW-(ag>>?d`bDf-lA7@GweUiHe)R{lU;|pKl_6=o^2evnD%_Ds?VVqiI`=C^=lk zSGO%8bGxF9$bcz_6lE_1J#KdzODT#$G zf*cv`CJ$lqTE;Wh1+t|2r5T41CBw~7%|Rv7i;fNV1wXH^_Hsv!vU(Yu&Tqv~E!rev z0c%2?j9a()X6&Rl?_8?L?P%s5EZ2hH@U6ExzX!1+O9sVi?o^AQ-l9iLkcC+4m z-CEt=A&xq6h3kOU0arVIGe`tf>@WIEf|?5f3@}LPV38v&Uq)`O=ACFreSLM5k7quD zb%asTF}XqASdT3ZnRxB${IXkR=8A+N30yqX_dk>J0(jxfBK=lN!ay^iUhKYE7jeXNEb(6J9geBE=qkW z!Eamf0jDeC!W7EE*11x+{O;Dv3vzv%6@ol#DlbD29(;@Tk7u5DPR8xLjwV!#EwkrX zohWx15;1SmJE>@Ct?5zX5P3#(Lz{lEzU@B6k##3Ud46pCi}&Z-o|0ckc4G8C{e_6_ z);O;kmBGkpq2^MlqqHs3WS-v^x)%RbDx*|AfQjNh#g2hzcif+2>u{H1{YY`I505I$OmocTS3! zC-O|L`Hy8^x4rDEt4ae6x?Iv_-;0m#Pi(YfdYrDtrWj+0Os$YmMJ@YpOx84Spxt9v znKpwjzS$L|@*BpNzIsS)$L#O{$^+LYe^(v?{@=<&>8QFqf>JboW!P>xG_=^n*Pb2v zd3RX1o1R35r>2g_(uxaT6t;;88{14k-QvAs=V7=lq!q0EuyjI|J|N}oIj_JTmd7@y zs1WQ;DF#pT9w$5jn%lo}cdBwGrkb~VUF#+|Eu_nU{m zQfqSP(;T2Ma0T}_g^3~lTVDt?+7PM@MH0p1zi6AO%^!fZ=94CCV`NHqQ@s#8hYmFsLX?cf*SP8t*Z?QX_|4%99@HS z$J4yixVS28J4Q{cA+{5G>g^9cS5wm!Z%nXpCPk84gv7fIl}x$ZjA7dKy546*yoBD2 zFg|zWyIXaqyVRlo%HxM0^yWq2gf6PO6*KB-~O|fx~2f{+iIHY5_nDOq~ z`TJftwtFXNv2y9MB6J*2<^}SPq-&Q{-OAS0aVwH8m@CP;`tJN$;REz_00e|1G5_6o zQrZ$R3LF$4Fu|hB+$&seY)W^8ZI|AZ5v22Lt6!->aQD*$C{w_a%BOfrv&#Y=M9O_~ znfws+5OPOe7xL7XIt;gvkB+B1c{~Gg%0kR4rbQ8UIAb1{FGb9Y`TCIT8kWUo!9ayC zyqGiQJtBtxrog2vZ$_byooG|O1gUArTV37ouG?8j-KU(plgLim4iE~_)=yZ|X^8hU zs_I5DC&Ayd<$H1Zl}L8pIyq;p@w9|1ZDdeuG;NdBt@ivYYt8E&IJVx+jOwRB7>0(8 zscg$+6f0a4w?<0hF`JZZ0wADul)n__pB`Hl7N7V%D~3FiSCJOHaAFsQ%K zpst!ry9~+bEdAN+XUjMGe00x@H+^-zAg@~_&=d}P4E5`dtAD3vM(ZBjReP??R5Jgv zh^a}MPgrAQG>>%6Sjr+tjTug)QvA5`3OPJeK%wwS@j%X|TyijRoaS;4`9nX<&XlNs z>u?2S-CHTRYz6lWzuaRN9p-deSsIJkc!t6}?JVPmZ*&$vY`}>eU88Q3Td}T@RujoN zt31MgFLOv?AxXb9F2;_V?g7KNoh7UJ?7&^TGG={+RNo#3*IgzX` zZ1C7ckNNy|Qq=RyI3j@}mc=O5G-=Y%A`U?JVt?V!tDrUDh=zpf zL6Ks*P%*SF6ovvvG8ch9T6^co~*>+*Nf zIUK+f>K)4Rc3_nHjUt-CZYtxvF25y5s+L<8RF5^9|Kf^OtiOAbqfmjp_c^tv8f{F% z7wLp=m!>2#yVfnHY$#uMj;B=Y6DF4%c1qZS2M}LAnyTzM<#4-W{IgT}=j}sJosg;r zROkS_`d0&whK_0k0pRv$8NSGHoJ(vfkWZ7J-5u_x3dRtUmM0`8ry7w{!>+LlNo~e+ zV8y<%T*)12H1C2&6mdd$6Yg~E%6GrG(8n&63J-Ch# zk%MRDrwXqxoVm_Hr%B<-N@bGwuo$qTPc6#9MoVsDUQ9z3A7 zZsdJFSUR;R&1*41&iAc`1ZOWBSsq+MSn%GZNAFz%zMSgOUH;a!v3b0l`pNd$i~}@y z0P4j2g^?;ox6-sgNCi17EY;#QBJvs^C3d3gyyqB5xak2Zvnmu$xNFd2WyS()ad@On8BCKZZD-FQ_L{?$Agfnj$3sSxKQ z?GQ$V%CvB1h2t4`@jVl#2@dg$GIkvx&jIL3{4daxR4<0sA4)>he^YOp>slx0hC6Gd z_pSkS;4s@A-)s)?poZ4cx8~gqIJzahsTYqkXh{pbPirJMko!KbZhqcuKhgMMM(suT zsV99*nbeKPna2+h?EoBw{EHdffvQSAVj#k#E=1gmKalK3#(A^#`r$+Mc9eVv6nhY4 zLj47PQ;b*hAcT-8Y}^dw^5348OgwD^tKUUaKD9=Br4)Xagv`hl?~a;5-h9JM7ML>Q9h*`rHfP>;Seo8J>~X~kdb9#zG%PecLd!sojwn! zL8%D8$dZj=2#Y(A9mM8ybYCUc$m2HHhL2oGe(ox3baG%r1 zE7OHG+{n?*htg@WN2eCdFdIphIw1rGa* z9xSC@A593O=rIdEX4K;QvB-`wQZ0XXO=asd=K!Okt^g`6fnU@!iN{#)^~uZ0kEnb? z-nm$#c?1+i6DxQ&`pG`|h0oA$jK4^(l1HZ4W3G@R1d3UX52~Tl-2~_wvd7vGUnre% z1O}mOm8OvXJoDiih(`zYQHs_yx=2on(;A0&%RFMU&PB3Fw|$bQk58bQLG2oSTftJ6 z>&5cLySO$_p35krBA{VSNN?gYGpZD!rEx&XYV@_32>UnIp4(xm=4zmj`5gYpJzX~@ z=fMTH+kK+%`1D@>X|2FIDr+c1-5NEsGb;?+K&V^L1C3sFQ9kluU+vz{8h2 zW_y5W2OuFB>Ti~kwJ=m|P>RuSWLF^mo<=pV6FKMNcCMDrWgjP`$@cENd&`Q53KpQ+ z!5W@Y&AyF{mJtd{=kk?nhYGo&RRcn04O9)KO*7Ov?h&0+*&L#mWYMq=yx9?Uye_e8 zj$ldGw9<19@8gj(+-6l2R+Xg^ER7kHg?$T?vl+6aG#hh*3~S{Na7)IGbSNdgp_FZu)wzls*vycDSL zhGue}x|?oaaK3@n`vX6*+{>}+jlLWjB#Y69t$KydUM#wsd5BdVmlhsn)I`}a$U#AQ zF8y$c+fZ~Wq>_c}HuX7T=D~`Kifo>>M?W>2K2BSigd#zt~i-?DmbLa7qtPPAB-M;kp(({!{ zISu_6=7dfsD6}#r9M>+e3X?dYqOM|)Up9u+-WqyT8{I1t{I>i~Z*ZT=Xy?hWp18=m1mWOI3PRJcN0N^Mm+pLG zW7{g}?W!VvEIM+LTJZDpYc?EH?aZ$%7&*ci-alf33+?8&}3KI%n zR_gaZKb=(`GqIJx=&+IMJXq#ip)OIJSHbF}g+0f6-o&o{BwHVG*XOdAG`fxW^uhIr zr`WeA1`U*|8NJXdj%S@*W4;-k?hkf3s}#W%msP|mnVr|w%S$#Ho5AKGKw$Acy$t8U za>l*jW-sbSWGzC~_CvEuW!gg@Gx4e&3O(C;cVk z095>&{X+x-tq07Dp!l@;`5|O8O38ckBI&T}Jw$XpcvZbZy#Yd8H#l?}$z#xkED&o2T~?{mF_B{3hG<7iZ?oAF_q=M;Ve=6|y<4n9^O(J6<})JOX{)1J$$4 zd0RpBe8$oyaPi#~s{CQtN1a5oHn+v(Cl^Tu>~~lgyCAzWZ+v0Hto5}hrwf_9&&_5dj@cVz@6x3gI*1g63 zLvm(nTtrro(Xmclj`p~o&gX=^d+x-|_-N3XcW%rw4yn1UWfWz|Ltn2^9k7uF_^m%m z{~H{+=ga)t=;eW&kIxY39#F4?Kn~(BT$N(18WkZ2Fv-lwjcMs%pLK`N2V-|X`e(Y( zSxE)!2h<1VNQqn^$YfyWV6OR8Ta+B5Xe_Pot7|~rH4$3!xQAmbTWC2@RUm{zJ2LfH zys(W+t^X*R^6uU1dbeXLB5y%sKUX_3u$|8}<@8KA74%3zbZn~Oc<1e*#@D%#%~$jA zV5A4s<^W8C08sdUI$o46#svV}AHI?KLHu1tmqviO&-Me#q7ixmcbsmLgo(aaPiV0I zy#?G>_4QY@V$&CGf*+A~-J$)c8=_L!WWe!aDGgp;i?So^)V=+1I&a#&EAUj??E>|?mOf^~6cCJ2v;WX2QgblnR^|^rg-SSWOnumwkkAL)kd_dzI0Lb8f zgUEz+R6Ph}NgDeE*RI(Xys=6V44vK`Cg~9-v^R06g2J|hdOM1LY_9z+jM=Vs*yj0>1+#Lg zG}*3_vyXV_nABfbaM@9vN`s43W2M9(vr0xO7I7{H#}KqJ$Y2$Y;Nly3PYGpdqso+` z(H`_)Z;Qs_Le4}N!WOW?&gc6y#4q$O7?rF{C$X8ii|{dtmawj^qY0RF97 z->iUZ#Hmur0}B}&84sBXuOV7yB_9+w*e^XOl_&NL*q(;CNlNy+bBo}t=b~Z+Z?I0a zR6(pJ2gG-~>ZBjRWXE%bjZ(ZiMtW(n58*KU8l|og$=nv|jj5v80gHlvx`)4c14t^i6_}|4?`3U~?aJi6fj;u<+)=@sV z)TcMRTX3hTXN&0mMUkBY)P4Z=5&Mgs6{#o*DnRX^sVm!G;U<;~3wV9?Hh*>Gw>9q3 z_5pK~c>LOyu4(C_L61IeiOWZBa7-JrnA|JtewI!VnI83k^3Ax$FnPW>Q_G(+0TYow zG!2LS#ik-k2SWucBra?)?8tCOiF#}0o_MnJ;7Qbg0bA?Ykm@9kv=!aNLXLNO(!OaT z84Y8@y^hG)4D1*Ft)%!WfOPcRF~J}FY9V|4BVPIWv0)^T2^d>rr8yZt!%Xux4S198^jx#8vO zQd+zNI0pFhJE)$cq;!g$l3OwPFX%p)AfLf$pFgH4aU^Jx_y*NHbcoKq2j4pDbW+vy zBiG8cs>6m~cJQ*4^^^b95A5fTx_@kL?<}0CoCVhAdrU(|N5Mc)!N$?Tz*En{1q*s- z4=hgJ-ogV=uBf`j_c$DdA`=DEDA-`FJRIC`t{xCT)!_LuQ9Vyf+}@!)unYjLXj(V} zN899n&w4+--2QPkVHgw+g9symbzhhmuptA312_jb$_5spVS)2-*gNwk0)>KL3;6pd zW=19o3a1Sm#RE>bfr9~XPEOGUIIspD)cf~<;13d+sDgv7tqpL@4m?+93ISi5!Ma&k z0Z0CT`niGU@0@IGao?jJ4tDl92>4c#m8-L}1q2D0xSNfKgR3==%E_ z;(*hA0Fwi70Ulm}h&qvpIzvQ(0$m|67?6oOL==z*XiuPOo`8CR+yS9{0QrGeTY-~( zqJ|FEID0^aAQ1TZkl5`Cf=(88SU}Bx)!edJEAady1`dqi zUnmA7zH+1DMEOw{N{6gXW8ocnaJl6A3jvH=7a)x+618-NiF3-7%eL0||Z z{O6;WKY)%Y1Puj_Q2zf9WdC3DKmfA@Bw)};VBiu*h(pi_G#Lh%%ZLF@4+l9p6$8?Pd13%3j27Rg zf!DptAjHVPD<}k5DKMY`;{?;9fQs+u0q_Fl?!|!>7r2=CUOI3`9|#PX4hca4S8+gD;B_wm2V_BIfJy+bVEKEE1}YGc0Z=FzSSAWcv!|3F z!hzR4`vFTrf;Ui5MWEPV)xawd_Syq1Zal|4%NUP*5n@-m~u+V9*%I zj{q9Nz72qO@N3x9HDDK@P5c!0>;bfiA7Rf{zS{(74?hP`uDyF8570J#1hD34C<3d!T4mx5W4{N}p_hU8?`1i8*$2eg8+vC`G3wv+Y0B+3|u$6uKo?ij9xyQ8UY5;e+ zcjEv7D!CWC0Q95%m=j?j66E{M|1C&T=nXcL?%(@gAVt+@oP3z}Ea;`h9N**vvj1aFD?JbNlf< zx4O5a1C+e)n)jAifIA$}Y4+3ads)=?qCv+COgZ+_djcl&J@&811O~6Yhu*;0wjc9> z>;ZE?6TW~A?9)J34jAR$&A;Cc0L%E{|CE4Ux1ZI&&%ttkJeP(23_1S!nt-UhEWjjc zWdSTGg1%N3{=eBo0P#UiC>YseM*JVlfMNZVE%JYI0ubZh(oz4D4Jh^&x+p%9zQj zIO%J-`yGGbx9tzvy;}@-5|%^15C}Scfl8~M_A_xli4+#dH?^^M7`>zV@!77SPeK0? z&(~q*^_=%PEVVSBvn+E-NaZ!!x31Y2NYbSF2Pof+u4bS47$4jA5MQ0j!UAG|i;}0x z6*IdNr*Xi31=uTo=yG3G``^J3RdcWgQ7$mpg38_dHiRdDsrJn9`^hK3ll+@l9{9ck zh~-ZJpxe#K73bt&3GqURz@Q@HLJ)f#&JBA~RMh#OXCkg1cKl?b0Q$D}v;w{t@pCda zYg>q=h1CV1m>*dHta}AGxVk8SU)eaJa1suM126&(1>7tOB`yp_@-%7;dQqug&aaZq8u z)4=iecN#bX{Xzq7(7%@p1r73dK7<$)7&(4R3x@(D;;%FW8kht8N`s-r|EwDl13I5y z(!v15@P}TI;5PYhd?*YscmI_JMTmhe?3etY2ytL;@*52a1%1b_d`LLpxqqjjF@KbS z0^augoezb;{9X?z3Ph*B<%bggQ>Q4*AG(5~;Q%7}EiD@ThaCXG?$5f3!T-=HR1Eot z9RMjo?*sn+u|KF7280T~(!`T6jFa(IQf9DhX zvmOX>^dIE{UF)ZI!+BUZIN5mYea%SE!Pf>D>mZ`KuCCxPy@y&cJ V1GhIU0;|hV;G;-9idsr!{}0{ddU5~& literal 0 HcmV?d00001 diff --git a/analysis/sf_c_44_plot.pdf b/analysis/sf_c_44_plot.pdf new file mode 100644 index 0000000000000000000000000000000000000000..86bce9e0f0f9a9e2970fe8614aaebb97f3a8ecaa GIT binary patch literal 19366 zcmbun1zc3k_dkwwBP~iSC6e0*^H!lh8O!j%({ie_3k1P1YO z^q`fKgFsD#eQ^+|iUYyH)5jG8HFR*r;UNfsz!;*SK#TKs1_@<;c|qOBn*c#>6F@Bu z96fMO1c>y`ty&O4!-U{Kz(G(u1%?g;0?ywXf(CDCp(bt)&hFl>5X_IK{yt77I0D25 zXjWAnU`^7j01PtsEO4?q24o z=L;5VSIM1YPT$@)62T@*9#4G^w!2PW$<^0yPY@Ryh*dW zq3-oI()XWiup<|lt`?mXy})SMvV^$<`0Te&iCcKbzAn}`hXck z`%w8@woCfRFm=&ppx88FAeeFLl#KLQ^w|7d*-A;)r4yoK%h@>W)?0(@Wn0hk$1!sk zH+{n%`)q8Y&UA9iy_+F~Ob;?~K0J9!f;XTiQ|1esQ=f3L8#jllD_##P;(Bi?=fcA1 zsZH-AXJ@T*>TBAr*=oDYjvunGuaTs9xreip+w(l;iS7A!?)s7~Gdg8=m-oJ(aFl#( zezA4Dp@|Y&()Ii-D`e3lq6e$;}yIkEjT>6;HDTn?VMo zMqML1K+~$ch*If4MsK)BVV@)k`+7gcA^}|RumLJ%sHx(y)(4DDeI^%75FW&Kb zUo1Sz@=G!+|A>BP5-vr(Ak7UOQ-6-Bp6-6I+~Vw$`ieN~lnTd``01NQ8d%|K9TN`*oOrPY$>yo=7kF3s!ei-kdB<@IWpF6AA#=u{R~ zd9Oa9BXjWX=g#a@L)$wv3MUpj1e+}RqDyJ;J!KM>0h&RKbk4?s`3#$4L&?X_-@M`R zrIuDCMebfI_NrMy+DPuZ>ksFH`g(e}ovsq5rHfD6Xh$|ZzR5lG{Lwgvp?lPMjjXN% zGzA_b@fmK{jcrXJ)}~afxu!C`j6s$|^rvD@)cfdmv|DY27HWzOb$fPektaQp;Ng%{ zuXPCBGoYj3{N&K0Tzwc!NuRtXYIO$gg3(bsFZMce{{q`|*pu)yD?j@RmgK@1_7g_e z$?T0cW6~Uy<47^bLS{0l{gO*UBCAYEb&8$LFh-FK^jR!xd6zTAGDGE<$@p9%;?Eh$ z6flU9SZZIdSwfxVVP+*C)Xj=H^6HcOd?~xd3kPOVw=^~r9W@P`b_i3bO)<)}J=dng z4_>O=9;6?ptLnSB?gPpmKp)lkOF3AGq%@u$c|?ibf0 zN$T+Mr_~a1-e1SzB;ORtp_BM^ME~I-zm5)tX^LT%pz}Qa%a7~HinuEyX2ti9s*vZ6 z;nD|q2xKL(l?2j?@E{$IfG4I+jbUtI9JW&L^~D6NhuPmNR-IV9E^ueYaV9_hys%>F z%Jk#3$@^ChjaBcDX47iOR@F12^IJpaMog$Y5A+ZhO%{Zg2cPesa$;g+v|AVIWbL;` z8#pmH?L){pU%ly5P3;SL`M6uwszou1C9ho0Plq+_vRr8}XTY@ANDu{Y#CSlzSO^x)t>XU?_rXu2G} zs8s03U@g|PPoq+=ZSQM0D;2Y|<}7)|2mLvV%AlmzX&qxgxV3h3Sq5PnS z$|MeX6=G#x;+^qtac+85{TdG{;;l*U7Ee=EbL-~ki?iV~WA7;Thm2ikc%WTWujivq zK9<2!esF}b!j0>=R-ufH*WN=9{FL`TRXZIxDQ^pD(F2=Xd87f0-qu7qP5k8a*Lc zjrbzDHftm#68T0W%U3{!R;*O1OYjYaUcJCjX~kivQtEXYl{=MzM=yH4j%rj~QA9=@ zH5TI$VN8F&28*HOG^(c7v)j7*%q?iYTDarWGRBlx6~~e5-nEZN&&AAK{E)FjE*LU15C%OArJu!7BbtpOKje+2 zv`rqQlc~xnPLV>87UjiN&Fh?%(skLknaKU69bbG$z7X{w{5 z^ErrOFb99WUvRI&esOC>NzwbJ)btHd>tjQ?S2xd4rmMY8HeEVNZuNZ8 zaqj)#(@N>-wP0#V2dx*OQ|=UZh1Lq=tTuF5Zeg0-Qw~5@CD&!i*!QYh*2=cN@n0HC zpR3sr5tfuAxFO~B32fi=v4_r*l807ExhL`3x1m-E8EQrL&yZqtw*=rCX^Jnxl1AK$ z8t;4UsZOVielo{=ukcBic?XqjGfQOOyQ0-o3=x;+KIt4Lajq@@h)MPJRK#b=cP1y+ zjVE6W5qY^7zxZ`&sRAXaoOFmLqb7t$RIm8z=2C@S3|)-kNZ3e6RN@3lhg{@7_(6g& zXHLO1r7wRzKh>y@N3^!>yTYf=sR&ulC4!>7t4H~(SMSnfqs8RwOiVU! z3ijL+XwR=W+4qPrYih?BeZ|)g^@}uZk$b=Q zfUL>pJzyjsl71*XALf%D-fH%K{u+;sEz9%_=qj1m^1CB4umhu*UxX&-$Se&Rv+ z7IOS^Bj4Pd`$U^ki{*-=K5p(iMQqXJEsKknUz|5+88f)}=2n{9DQEOs@AqOp-YhFi z;{!``24CyZ(6h_Xt)p$1+mK22nZ5Hq^DoC9T0qb0KsV=Jl9uqcirP{h8PiyoI8hRF zDn2!_;6f4kg@ORkQ?6X-yf2TWJ(8SvU*9yIOSBikqfSr&$ZQ01hgKhVGL;d`~{ln1Fbo0ma(=m zdFtirrB<3kgMmutq2sZ0elN;RtzVOB2im!27zu0f5)vDZ^yYrnn%g?+eWuZM+-!xj zWPVlgy#TA^e#Tc*hz%Rw7|NOX2*x+gTV$8XuhV5DX59MrF8-4qJ>Asz;`yrYdJK5J z4Xt|~^Viyg%)|q;9$TR1XA}e3+_}Fvd96PY`}j@lqwxBb+oaK?EXUl?FOz}^e!<^f z-Q^Dd2zhrNIzyFanbCPe+z`*hAceqTBw=dy6$_MGN z4rfyuKSPUFGfK&|R2cQ+nce8<;h|Wjs@bP&XNu3k?Q>`utTa0O<)h^@J=u#H9U1G@ zb_}EA=(^F$8x9270PM(_MNjPJndS3v^LTv)6YFRv-7m2%b{YAt);`jWM?O9F@bsAq ziMz$;ee#{I!a@jUa=mEu-I!_1=la3RZ$E3le({=R@lm5j<$H6HIkASzRwS=0Y#0XY zrToVys9uV*eBWPxAu%EUaLnkp^p|>9Gc$b<)@QX&aW4I8mB=V=1vg}xw62dEKZ(6M zu_E&B`O{pSZ`~LBh8sV@zBD?0cjDu8O}+Hp*X2Q`CzBY?d>r6>BFGs} zC4m176}(FJ5p)`VxH>=F$8Em|^p%hZEasWW#9YRk&a#>W4bb3p{C{1Aeb^-%)e%%#<$R+z2T5mr~`UY9$373hy0R?$8f&C2s{PwICb;8>^ni%BhuH;VNHp34W?0)m!cttUPCKuj zK`R<LCljV%lcHy|4bFhr{fso&nYQC(V_9do@hlVJ@K?wDr&BDw z8y1swSh^teD(kt*8YEkx9~vGmg-@Kd8OSC?usIDtRV@h1=`p=MX$k9YP6m;Zwl4|$ z^-Z{v;#yqn@&+6b4}Ds5n%YerQXVfCeBRfZ|Dn<)rp62jkCe`4*>5j)!k_8m#}b)0 z&r*{;EgVKJ?Jul{4W@jR-)n*&!JkX{DtOHQxWh>HJSV4Q1y>bL_x`d&xu+!ks|2MJ zcRxM~sr~?CUim=wTuH&1NS8sM|BpdIqIRJgXk*z27nP`A-!@q!y=*l#5S-gzxPjlo zkjk=4)GCucAf?IU1m*X%dT`h&C?>Jm}o zy={ELb3qT6A~;jFP`6gg0$^nbrLGgp%7YPT`vn~cJFaos1#;Ga2tG_>8`o) z&UzmH@mjfJ52l0LQv~^fO^tkpGfha^k5iAQ`~ILMzGcMFc&@D_=E8-e>zNly+r{~F zyRAFnDh5{b&Lv7Trd7El$jY%cnEpd`niC$)x5Phv{4gjK|BO4(T)Vgf%lH;|0UPdI zBYV)Zb9V2CDf#e}Ea}^ElHuZ+*fHnNb$n_EmLg_8=T?Z;vkQM!Ntdge*}rw=`IpTP zU%ykuA37sQWOTrdsy{~uCcTR}ld9uQ1_u;UAxdg*`{(h?9|{jOXHA`bCL+9ao_m}_ z%_E~(*UE%98FeA@?9!@l-KF=kp2IFK)7Ax3$XIa`?ug z5tfq*qgQ^VLZWC=Du{9H;ANXn%tu;B@^~MirD&uH7FThiD5vr-XCIvW8jqOv(6#e4 zt~kjjcU0AEagjNmJ~^&DyiqGIyr%gi8tE_c)~#p$QQu7R_o)POouq0v{tweN z-lJ0XttN-n9G&(3E8XMF^w$m99yWE|qTl<`Q%yw_pYr(Cv95b`$1^-Rg0w$V@Lruf zsKaPJ{;X7grAlsmWS!4?jKD`k_TVPqAF@ZFQ2#A^wG=&n3V!Oy^~(j&)$%zD`Qbui zE8VxZ8HQSg#I2*|$sUrjYOJ!TEXSs^@*z$^7ZIm!AY@k$k7xTD-1N5X!gV1XY^G(! z24s(_oM<#Yaw)%jWj}UG${0<+IbR9M&!HPkK3C0bMG(&nP<)v7*7(JKPn&SIkJ^~|X^#f)5NapmmE(cx7~v_NIm%4(rkL)T{-+tx9DNgp^j zqE0;D5@$R;Myha@MCZF5B@y9)Tc&>s4~h71;h~MRyzjv1TfVS=Z6_Ih+bYtX7xQUr z%=j)FwW4t2O<_kTK~gAu1AA^_BNcu1;2Bqci?0#}Q5v09(^_m{nJ`Pba_f%ooXTJ#;0)Q^Dh$ zYUUx=>DZ2rLnDmKO>bOh%=q2xCQH1WO4YS{AKqrr(@Uk^RHbokM7`zses=g)+t?Q- zeSTABA_4<9mH!Z!H1faYg~VV03k^e&r>Lt_(iSV4pXcY)Y#-@U<_9ElvETouy*0l5FyxwH71ox^)}6$! zfW?&vFGO$^0=XL_Wq?g3=cgYCj#7Y`?R%zpyQHKM^JMdEb85DJ@+yD-Lsp*bplm_G zT9@8u*3J-@X%pS<+n*YkSj*R^c?Htq=p3R`d`By1y)P$nZUv?fTT-oHHqKie-TS0d zchrKzRq7dqlY8zRY`=HB9cx`eKvju}r3NgYK9bPbl?K6^7OwRtBIXxmQbaDtt$cbW z5r{p*_?ixrz|X$0Ar=yji6~eOULoD@ue4=1xxuT}vM3BWK*KR=W$UOD^nn5`rcc@cqo2o(0e2T$HmHc^$I z{xvoVx+?Tu(3e-eLegdBc{N#%sG;teMx@XPbC?D_JgsI{xGJwY>~@^;JMWp-5uK0< z6=TS~5TtUW;FO*yTRC1WVJjGpnA=K#P`fJXadx z6>o8%!0Af4zN?29o%QZj3TBT#GqF*i&cy^j-q*;6`r4y$CFO^ZINohDtBt z>EiM_tGgJ~xP{01(wP93^vQh-`&cqv&vHF?-rA)49%C)=w*-+yuqqOZ-h~8>b-cS3 zsh=&d+2%c1y*M0fY&O~Y#pk$+ahX_aEW8JH`fhU5OC38Fzo`BjN2_h*icd+|SZ4>H zYl%w`R%o2ae9PZxN06$M;n8?chqxrBS~^sIKYv3xBZ?|ne>I=3^ECGREHt!Ztd`;C z3weZ6tqQvlgqv5%I=lMx2oPzFwpyaC_`7$yvN&w z>)6+{_(!J*6k=tl_i67(%<4zzEMtRNK4HIVZY+#?@<$)`9(^DCk`LL9V7X_ngSF4k zoYk0qW);9(&C_)A9wCe`ebiUC;q{rmj|Sff&#Z?6TK%aD6kqBot+p^$tJmGRgHX7X zAR1C2QM*5EsRWfooV}7WvfWeA#H5sWFI_?J-fEVP@IDvcqCF%= zZCo5gI3fbR;YiGXYnfh>5+fO$nrtOF>WU&nF=*J_jOp3dF`XxY59xUOLp_Zs!-Qjk zVh+5x|BUIm8iw69?2&mzmdX^m;QF( zOu=}IRY+?8_*T|}fQiwDajBG}vfa4@!}-js^7+SEtm6XmA1#B9 zMJZ&JWe2`ZRStQfCo3?>OIwJLl@z~z^3E&ou+K*Zj90()uW#_wFb#dR%_Sm1BB&F) z3nf)c=wNP#P?K+&D$&r=J~cRa0f0>jdl9GX7Z~40QqVO9tMgH$A3Vx#Nu}*HVVIxI8${{#h-ugKY_Bc{0-4^OzKXqkyl6tV5HHo5iW_qyIpyW zYQXZS3a5T{+r&3a^F?OzhlX;{EwQ<%`*>)Y#@;h%z@*suZf%pWpwL*F!D#`N%Hppm ze<6Y54J{gq!1Scv<2scI%PFO7SDwC`ty@zs9%r}y`1XF$P=V3V-9BvJ-u%n#&DTyo zVjN!$vDIu|NgLwl9@sJrlfX^-(+VrujGS0ze$K>G9WEa3P#xlwhl&;(h}ANiaxHVZ!w~7BxcQVTAf$lhj<4TcwRMI+T2G7-%kRWI!C zHWOG^OlBYVT3 ziHa&?uelVbt_0N}>HQ2}T7k69c*t_&hXnk5F%(7Qm_YAHJYRcp+JeO_5({j>u+#nThnyRiPGHSt4^XhF9m?uATdX6kkTtoJ zeLU9VwG(dc=0jP^?gRJiuw~AESK|Xr1H^oFhtoeAO>PESrJdGb#L+&f#B+O2-}UO+ zC!&6_>Mi4XEW|*>M5Fr+WPFuMB;j;=g}Cn}&rkPtV2tG7lv%Up;uK%ghxH7Ztj)13 z&&~!egxUK9Gi#3T?Mva`#}giS$#!dIbo17hU$pG-b*^gem6;P+qk+B098=lxd*UWU zOBaqyS&wlYp1OeSoh$wleXq$rvn4`zW%J1B@aOMjyrM3j{}v^R2vUP%cA;Hr$+9eD z{PZd(_45Tc1j!VJT@a#}bswiEp2qDDEg=yttl^KIhQ}_3$21C0Qmj;wyE=Vc!g!-u z=Ews=%?oGU(u_|xb-Sk%-nn7$;&fH4JT8nOX;-=O_8w-RdQ&jxU;Osori$;!V;^rj zPGT>{hV|d&p&oi=GO>w~Eog4;|6bTTQ8{)-eA&#w-AZSG=EkKK7e$&06YSFp`6GO^ zhd&YND-i^PMC@W_qm5d;x?%LoXbMgkf0f3wMG#>ssg8L`3IJast|>%MfDHj*3`Jl(r>>Q(Vj zqjuFi?|17=%~mXm``Ny-E#_2YSUQq2#6?w&v(K1Midx>R*a#GAqN4GKTP>O}s+||T z(SC|cVm!q<_qD$j?Uwgx()|umGRyUGs?)*kWV%A%okKl`6>OGYy;8#jM$0UXgvNzl z91_~2WFk|caMR(7N6Y*|_8}1&>}X8^yp2JBUdMRKqN#T{PiblcDd&13(U396ka}&B z6cRJpt{$myeYjnUtWm4wgp>AWuCr(xIfnLAsR_^BRmbH)bDt+)dl&Dz7bB`qNu`j&YqNvTuqa0^JbsUCdy4!=?msIQ@WK)PBqMkT+$n z9C13uTWgw|p9Jo6s2Pi)vr|Q(0cpZkCR2w`W!zy5j(+Lwj1d-7gI=f=ULT=-cRF^S zZGG}_M!gCu(+zuu4k=cCfajqOCdXHdwK;F13;CJBi$G=`!&_$q87Z_MYlQgsGVP=9 z$YzZbpx4%;+^Y6Zx^gt`fI`G2{;YuoVv}hSJ?Ze$fUe24^6hHH?m8G@-n?dw4GnrS_X{H zR>Lsy((+Y@$8E336FwFje^_BkyS!XE6a9J*sq%goBCHWXS#b1jHu?>)j9oDLXWwYg zKtcm78>cP^EGBuAbhhuBzEr z-hmw59l}2U5#|(|s&U0k=k*>6){O_7XWx3S1fJew-woS7efs^t3Pl8-!O^?fFG(6; z0J^8Lk`qoxJ~Dhw&ft7r3}NfIa0vx{7OY=noFgJ7UPae_Im_vVu&F_nf{@^iW-aXA zGDB%{12oPVPn!w(1eZrR2UzD=(BP*N|NF)skg){m)SRh0=5`9AjUl z$)tYY-Ej6S&PV>whi_K-=d@ktyBLZ-C(>)|+ zdv6K_qe@keR?J9=RR(vp@$g}LNqsw!++{aq%_&UYP7YpZJ}N-=ru*YRwYlU6o~IUinj zJoZ&uMRxI^iG<$BhbgP-`F+zT&7>baUpse{G{q>~%pr_1S@x#=#`y7ro%KHGdi}@4VisSZCAcSIwL!2;0!OapxRo4h<_>x-N~Qc{BE2Pq zUyu520;SWSgzbs4YnMp4bqU$A@s_O&y$|{68II;qR{D-6WI^f<2o^CNrQ&*6dqRyj z;KsgpZ8kk=Esuv!8=aoFd&5WFaX@6_4s0C$%he=A%> z(>8p)a^QN$=dCiTQH_EnQWSSHhu+-~HaCT(I`#bKk>iUJ+Nr7?9H~5pB~IsL4{7RZ zI;Cnh@VDYt1cHOXz&>YL8^B+FGRUN`t4!$0Y}}SYH(rul@|5~W=&`vgHHp)ksT}U>SzZsTLuz$p z%L{9{Jq_?j4<562Z92?5Ox6FX`U$gf3n^Pv6Y?JZ#i57h>J1!$7)=jbPoKnZmf9mx z-nQ!J1(UCo3CQIY_75JUok_~&^%o;^c$HI4@IPSYS8{m}eKGC^Qp@Fan`T{hXRw`2 z{dam3m+QeI!~O-2%BWZBTy;i9u+^7uWmvfpF0!^3Mh^Gr8|b%LAYb!i*=fogtt~p- z*|HsJp0OM=F*3;s;rjmN{o8l4XTI-|z@NhZJuJWi|JRoGZccZRNQ?>K7r{t@Q}F2D zPQfWy2lmttN(_cc2+}2Ia7(GsIl;UfJqzt!94VcxAaF=CcpRRU!%9Yk8=_!t7qbBk z_|O!94?Qe;Q?Y2(;NFMw^0+o&i;?vfyDrRrRupRWK&S1I@Ux<>IUBxp;X381q&GQ( z#@*{2^=!kEH5$XrHzaAvjacO^v!Uz)kA~M4510CCUkS6IzQJ~&HA0Id!efJqw+IO_ zTu|(3aYPuCg$q5wyZcec`X7uvTiAN}xky!E{+hg1_Q3NCIz=Nqdisldy6<*t<=JU% z=4EGTtEOhYK3%807TumIr=S7v4_(H@dKFL{T;MGd9=yYLxVD$H)g|Y#*`j?XZ;WWX z1zmk9ujhLk*7QQ2s#&ga*wg^b#4-OxEybgGM_l;Zm$T@a&U0^!F4%SXzRegqL4Dut z`vH!-(@yNEp&BK1EDY{4?dRGrAC>7o`hgY8*nZCGQjP?@liR_Pki^*bk<0Z(4n4iR z)*nAGnJ%4*AuNwdYc0PWt#&m&6)bA zV=8+F_$TrtR>QT#qWKNuvi7G);=FH!KEp6vyOwTJlvo>g6_)g=!IPc$Sb>c|K|4TK^tSe0l@wG8|`UGsQ0s#alrTa zdYhqaoQ=$nV4NmvZ5lWb6J>gRnXuWA{**;}?gTFC4o!ar%Nyfp&C*tL{>LlXh?*N{ zSMomNqRzR(Ilun!BV9!$biJF6bu};hc#?Q6*#+}e+Y#M&7WZe|dt2*bb0)*RJt-(M zSA0~m*OJ2cz@*HZ^XvJ+lGiBD-OA25Ps?K-%W+2OU7kk5=}-GiEg|>caBKO-#~P-D z28KP=q{p6=2sGkbLmq;&E#0(FNH?=0P`9E-(zid-*>wUhd+e z(f%yQb~{O!l_#z5IU_^<*<-v{<8v`nr1J$wyx6H2XX?9-HifZka~@A!-&3Z$7`C)k z^X_`v*cczro6sI2gF^(6A$EhvB#pHE$&{#DhQ-$+T}qxiWs1klZH-airJ`hrQzTf= zt`;mqjypVYHUsrIP7u zyy(PY_X{|e>y0toqrJ;^HS*a?{gt+NgjqRtA3F%TG9JlBNHyT)r6CLImYEL8-sby} zEQzQnP5!9zi-q^dm004d)e|uOY+s6?NrY&#gi`o2UefE>u%68Ek!8!u_j75yc79SK zoX|?{wT~DwP9yiq_llbB-@KVV;;&tJ9;CBykK^I&$)W!Dou>BP7stF6F^^}Sx1O8* z5?`~hu=$njn7z$8A~Ga`EaAI3pU_LzVpIV(GS+iDIZGpv?F(|Z%bVR+Zda*L1%!Q_ zgZRqHjra*k5u6v}6U8rb&$ib?oM!IJZ1vw%xC2+3EReL!4D1~rWWjeL;G~V}{o`^4 z?M%J4)GW@Kd~NLz0lTD4^AA=R)`oA$P9`-sxEzx&Xpdc{?ylWq5&EsUH%zB<^XS95 zOfg|6)|mHUuLO)SviWR-TVh)9YXkMc|Nwbxy*C_Iv&2l6vL;ni%U8>`Ad%|BNK>0Tcpj?sH>6pa6XcC^N zGPtq}v3Au+td{4pd*Ovy$E|d))+-IHRlV_v;B7uKnaEB1H-i8?{tL6fck_Lgd{PJC zsRCckZPIQj7CsTpKV?wj5%YGdMo_1PE`-!6K5As-KzUAAwEe}|?gzrA$2Zuz6@{hy zvh5j*3Lo%jq=Y!i?6c)X4CDrBhNRM)}uv_TKN)z{qVn5b9JZOqnqOWqV zw(f1ZQs(*ul{$&4PsQJ4xP5ZiWNfRXqIF=C;?@bN??ie<1p7$u;@nF$xdjrRUIg|dhK!$PvB_Fm+lvp8UVz54Eh9QwGN z_}dK6C-uhYie)(4iSR@O1|Z(|?Lqp$h55?a%=k2hK zZb;+Les6rCl>eoPLP)k$Zu7+0pa*Ipqd^iMdb6<=T9w)!>vArl!#(Tfo*cP7Dr5Ic zjZ`>!M;;ZK9IBllS2%yen`zSNdRS&7|3WOp@er4U(omTiD~H#xNFn+*u#Yp9`DU`u zq^KujUl`>ca~bZj#j(|Nd4q%Z3G9eR73hI|G^}a@a+`_8k6B+^qaGkw9`!I+dXTls zeuL?gkg)m$NdF_Tg-7^ru^zkhOev`~;aLm-qBh>-+WU?9#?%^L`L1Oq+)4KV##L0{!y^<2OBVE`%Rhe|J|m0tEbK%*n^g%K?G{RNNQm@9yIa zyyfYDcY{DfaQ;5DP#cgHg+TQnG5|-05OBxM1Ok4#5a8{M^T#{+_~Rg0 z;E2Nk0(AuHfItCez!o7;S0G9hpmPT<(zMXy5U3{v>IH!U&H6y#aG)MP2o&H77*(Ky z0f3N!1_0KAfnEULKm!q?Pz!fwf*T-65PAL@Y`Q(RKR+h=S9g9H{q3RpA43g{qo;!_ z9uW4f23-m71V+(f5x_bQwgBEpBEVyH6$f9<9SOo{p_V_hK#VRBJ?!qJYcQb!9>M1zJrYov?JH6mr2oGX?a${(fEG+(9X}Hq_yh|M z13NE+0uoRWRt6#sc+?nRTA+c(rNM+oKoD3oEs!7>GzNl31C;>tR2l`)$^Z<&fIN5w zpCf?;L`O>l?}K$>0aJ{T*(n3B+nqs5(}Gu!39wUO!T`z%zKaGrzEcN43$(kP2libM z(lXl#Esfp2VgXo%0m?CO1kk~qc2PjPXkZaRV!`)N5HxU=0fYr!w-XR(4GIGk0(b@6 z-ySqjfB+AGKxx4?(ZD;~Lix!Yc-__?uq70D0|P|_vJG|(yaH)^Jiz955-s}2nEg}L zP6l9brwRnnW8gJ3$Q8gj_&TTtKZon*72uZkCnMW6pripJAfN$cAJl_?>IbR7!fj;& zY69>7sRd*T2BY16_M-wU1`GL_Ktc-Cm_8bF417isA3hDrngp~ot1k?*)4f-|x)QO*0P$c;2K~r^Mhbwc+XOz4?GrPAY0%aI8|K?p{cA7*Be_ktUGl?j{%HA!#RQD!HiaJq zIH2Fo`9lDG*{%l<0lseBF59t_+dF+gTcGI#+|ZpT+Xfeyp*wfmMtghn1)w890GDIC zrXBl>{?RmOl7TgA`;7o#N&Aue*N_5d_S+AEW_~ab3~1y|RUv>5?K}mIJOs1T^1oK( z*Nq=_U?%Rg^Y3$zqn|cY3H~b>`R8{kpejlL5!A^6*lGmruM*;alSu)c0x4l&W}6uK ze-Hx-_e-^?|49jekAGW7|4%ZY(O(HM{|6y(hPutA48Y_sHWwU`AJX_&M4-t2Dg;OA zr`+_wwHLVm4UGIM_#H6v^M^qG1)u6Z{$8}upYT#z3XXsy5ileQm;{o)-Tm(y0Du0B zAwnVnAF)*ig0<_mm0r{P{|kFoqR#<7T;4>V?L9K zdwhl|Xlb9e_BxD`6f3p?g0Y@eXuWx`Wf5F5vNvGroPCqPEq+G>{YM8@1!d(6TiiO< z+)Ct_vqQr)E+;hb&A&-W>gpubWpr?Wm=oev7z?EBDw6exoMeDH`%{)XqT2a|6{wE8 zGl+kI1sD|W_V11a061$~4S$5q0Y?7cmsNn@4gpdAApnm1dioGN-5nu;NGUi>N=5?W zMj-g&4@03||2&iO@pl!ag#x(TIlu|{Es|ef^L2KCI662T2b%d=1wgvpfV+>kD){>? zhg1(EUN88|D! zKZy5_GPpEoiT*5ugM#>j4i2tzf0Us>ZTzzgTx|a=16QEm%77cp-(}!OcYn|!rNNc_ zk9QF;V4?cG3@Hs*wBO5s?WVu`1{lvj+l8YMfJOh~U6eFvmVRpk4v6zlIv7$KTy=k= zg8_i<&v$|EZ9$XrI~@ub&%er0pzZjBPUi1+(Fo93{Z5Al{{AYH0T=M!=>Vwp7aax; z9<%>JhXU=;A7#=2V)~;@8u54Eq>%s|{F4q1IJ1A2VgJ%em<$X+Ie*Z}ps;_*R0abY zXYlVIeTK)IIp0P%THY=IegUMhV@L4>2Yl^Z})w+XE625iu$ooIDZq61OAP=CyKUqoXe<&ay4l z=bs|`zk`>PgO!hypM#SJ$jQme%E8ai$Ma8zP2yjK|1RSW6!vi7h6a{Q09w5`Lxz=3Si4*#+t39@ju1pSBK-NOxJ<_Pe~>&xG0s9SN( zzgp3_FNZG!-vzK2p)hzM5wK%IU;@k#Llx(^L2*&kKjd%k2~DO)+g?rT);B#&km;^k3oe2(lZ47?BU5ge8yot^Vcun|tKFxF*m=g)^J=i#U zIjYSn6M26GrN3o^4(15BzyMz}ETrsLk3w ziMR^*hrOmZh@!g+Tv#)Y-Z4-Lw0BbfBL_rNdt?N!v3GbYZ+jQ!esC z=9?*y^=+qtZLNE8UZ2!2^=E4Pi^a^agpuvmrRSGJJOA>6`Ro4MEtBI14u6IE{np-x z?l41g+w$~BuK~5c$K|KnXN9|4k^Lnn&8M`|M!rseKR;UBulq4kJ?$uv)vok=u?)>D zzN~gT*Jn@o(+Ew=`f=zHT`ExMbTeyV|Z@TuXgo^*MM*v4=>VDY$3a zD3mp9gKN)0NxXM{ZGYolKga|3mnIo~Z(|JE}uF9MIB>tKJeB+x%dd|C!C>c=2*+6n(Dg$aJyRHgfg5`DrTfY}_}g_lq5y zl6Uq){8h$XwYRbS$*-jcncTk9Zh*9+J!BTNd|vPc_4d4{P0xc=@WtJ_Z&NGlwRJ&B zrt4?~1%bCdUuSX>M0z~I$3bt7i=LlwF|XPfJqh{Rvav{!pTYboow>b|GyrgPL3 zDmC0(dh&*V#t^IkGP_Pf`6Uco3+>9w^>!AJ@5eKKunnYv<91f0U53x~hzJto9xHL{ z&>5i_d^ui<_igJ)oR|$?HQA|jSUwlboMdvO;ns%A4KHDQ{>YH2@@xresrq0@4&CS} z6ZSCh_QfaUGM4jZ`HTUrW4t_X8!%qxB?Z`O(|bCz3$o=U%*m6v?kB!U^>VptAqr*D zlB9$**7t=6YLT^IKJu*j3%;Q1S-czMkiOHc=DIDeZwEe|_;=nWo6T+5$1&Y=MxQr$ z$$6H3ef50EpKWc_gOFV+?>$`zm=?l17d@4IZM@Cx#&Et;?gKM`x9-EcY%p*fJuk?f zI)xiEN43YJwHSYmexkGeHm*EwwS7OWR~!})YM*$uvgGGY9A42hh<0bGZHauSe>)B5 zs$aiSmUSY?A%zPd_96>xIeNIq{kpQ}#3H0>R5>3S8#)s7X-#j74u8)&UziIl4+INcQICP)QDdR^GGzh89x~cIrBlVQ%?qPi5 z`bI9{FsiAM)Oy=uITuaIpHn@^0rAw#DeLC9TSxP|b)zfD(Du=J*y(Zd%5_{M&V1O| z3f7ujwnTrpJs|2jEk_FFTQKUW4KLlHU^{S+yr@Oe9$3tEED$*f0fUxd`tPs;vAFz zYyX5XM=RY`r2a^4JmYkyo6usInZbH?Yl2!@tmM3g`EW8@?5W{lp35tPNr7s-lcj`B z&$=E#i;8AmAv7}8|A~jLM<5?sAW!H`lG*IA>iC$S&Rb@n@%T70WqC5|CF;-Ci?j=S zOP&n@+PL$?-{cjxl~Hw4;k9plDxTnF&|LijTiUBtP)jNouMBS&L2CR9!|C6S+pQUkmVd)OM2-)1pV_SHh z|E)uB@>r@PA0_{Lv++>wqc=X4U;aq^Q=O-Tta!50icvlW#fNZ zcPgYys2#kgc7NH28-Ng-UGh^vg|V)qM51?7)sWYlF-sdx{y?8rYL4qcU|akQ(nFG( z7#sFwLj@LGt2@YtiC3@FCq~$8utH2KC&wOwj$NCVN6XTeCl^yAc&nJJ_0_g6kKGI#E z)gakmso%0Cr=COre8*`v(|>01MO(q9=@EKwIe##AB5@Om`bxx_0&<~kh2Z9Bi<^1> zHcvgkk=KqXc98t6j{F305)UQM__kd0(Px}h;5D?`(^g#kc@Z)Tu1swtQOPN8zwn|G z&4*F8Dlg$RKBNA6F8mrT-Y9>Q`O9PB{r6Ps5$+czqU*rI{zg>;MlF;obcR8ii@SjU z>;}eY9lMf*apWbg<=b5{Yx@A4I2iCalMi^|tsW%Z&_yWO=2swhv!R4otByrdix3(Y zxUm!9z~sU*DQ%B%(q&xN%g?4@tX_M73^oiY{MuxDKi_hIm-!{Vbze_?hXjQuUU7Mk z>}-L(T01*;1I_<>9xQ1{C^h3G{G%h844)8D=#+kLFoJ$`o{_IcPV=ezUB-1(Mtm?e zLyz^vC98l_l=;Y!}r&y;^cSB;LYuoUm8Re^wPoB2K`BMnn-&q`YL0R(xXvs z7h|zQ73$JKpCl7T~27@1C>N#aH!m_Jf$|CnxS!g5WV-*Tcd6dhi12^sY{{ zhxSwWKU$(VHzXr+9YOh%01a;LdnwYc$C@m}P5I0~s5n>h zs$sdvL-}o_C*Oj&bKJ>^S_ILk340xRBrj?0VWRVidhM1UAph10IS+<$IwEy>a9Cj4 zL=KeJaig9RDuLY6nG5k-rti<`H4#@AUXmglZA)J=#cDhnrMX(-+awdyuWR4*QkV2lB zc44brum{@v!TJ5SbPj=Et< z$=QU-KQ_MeMN0YSjGrZFMP9x1QNN|ma=y+fjxwDBqdf2udlMrW7Pn}amR#LS)@kd_ zarz_^cL9IJa~%g~3im^Cn+PILQzYp-^QND&eLKa|)&K4HTvaZl1dnNHer1`w(tbgo zs%~Up5_7>YQvkO1$iP{s?VYHsK|?GZEtL7}jfC~*2gRgqcUuQF4s)6BPJ)Q6lhc26Fmdg6@l8sbe8)1C- z0|!r6DI_qv8DhO#-xxgOH*D;Xs9C*9isKl5=`8um-Y-~)ZF?R{IYaKFIswbll?*W5GmI}TO8_L&wNHLWLjt>tOFMNa`$GQiN%Wc+D=DYr8>xlqj zF~Cv!`70qS1y;Jug?&%ckmXKUJ+m6ypraO1zigxJg=F(u@Iz$Tcy~#6MZNb%XWe1< z+N)=uF}_TC-;vHlx0mh7eZMmtw&lQRN8Z!fsiYLsN&Ktswk6aVr-8h{vt#f`?zfZH zr6_~LK$d0;K3s~fa!AkntYKV}`@>v_9G#IBa6;tIlI1#I%Urs#oHf)WFZ{{kVk+xbMe^IwQJ`S#^nQ0)^}WNHCbU+Y>sZO&d#@+*-^B!UEF5l`V*V?qU*JcO>w zWk;fOHMD?F=`4HP*7~u6+#jslqA1cQ>*S;8l&T@(H^NPp8UlRf{A78B>9)T^YO*Hw zA924h^v}b;cBc2HuhAP*lGOvv;|%p(Exci4!_;-X8RmhIDBUThA$8 zn;xY2LnDb@SJX*A(cT@2)SkHkIt97?xSR?#obd+tE7sHl*&kFnX~5;*>w9^H@TUv# zfQ)PqhW6t>N$AX9koM^%)rL;Ahu~SOgmy=`8Xzolb5`jFu~Af$K6#Y6!{t>>lb;qA zH1~IwxtmWVBSe)AT#$xVxdJ&smnAPeuvA{Mm0#cbe-p_g>h^aVBnsFt3+avid0~7q zQVU-W_Gp_t8}Cdrquhiz8PIBmT2wk1vLvb?$RT$>yLAF(bC6`8@5_dyZJGYcK>KkWL?#d zq=aBs3fiIi61tPfdBL(i9z?MQkiMlJ^o=$c2_nUxqph67MC_q)aKvHD6+xyJyhY_oVc40jg?_t=9Tf!djg`|Z^ zDQfiyAmxVC%-4?to9e|Oo5>_jm{+D>Q5r@N>DN&frPW2V{WX_qPHYr7b{v9sqXh>zwb^2CG1pj(4A~vC%2A~Y4`**DzSfxIS6?Aq)hbon> zV+$?l+(Hzfoo|>{ua@mc#?RSKH5+0r#MQf2sZIFDo0K!yJBr!^jZ`L@h*-^}lRZW` z6m%ra@y5RG+@6l9y^Oli1~Dy(=-FE-x7

1|g&@FWA9b$>sbo{OByv>k)3E*y!uG zq^dN4eou!@+$ly)Lz8BFQ+_y8(oW`A53+1$bN0na&aom zOdD4|1%G>W!Fd%X&9FkVaC2Fs^{poISYm}Z-iJw3QH82~2lg;i4d-Ag{ zEnZVPT#{!u@Yi$fGrd?m0GHL@7whY9*6o$%e*?L-y1wiB&ibk$jiUB{8CRX!cs@O_ zI43+#dCFEGOKhnvQE174`u;_`g6QkT*6A=T|dxzQZ=gcWULlH#F*UoCI(4yI2-=3rMY`Mt?t z2NLs(@aK2{o(km3h%L&P+fl!x1Oq=*oV0jF76LF&?0x;_q|0H93)bSXwLkPCThSnzia`5}K_1T!stfL} z2ukxQY!r=2e+SKdUsNmyhC^*K<&8G!!ho}DJxb)zF%~Uw1;$dS1t{taK~Ord$V>ct zN>0Hc&}BG^nK*7q7gyQ_QD!Pf3jYonUIBQp=KdF;=jn^Vbwt)-Q?BV*>$Po@Q8T{C;QN zP5V{Z#Kh@Ly;yyyZ|zup0&+hOLSOdn4FDjWC*V z9+{OT_(gK~AA%*o+>fl@NOQP7ukETc;p)6Ik)YjVj`^A|vrS=}vzc3J*@Ua{XD@c~ z7>epgP>OK}rK`-?&70!r{I+`zk!*euTkKhJ@o6Sdy(U=$;v|+OUYc6F>GVb`+E>EM zgSlaxZ?_4;+ap`1v;Y^VYq&SgX#^@p+)pc9sq+cirrq>-(Wkza1f`qray+rdzDj~l z)Yd|a5^AK@2po|E+EV9iSJ0l7mP7TIKMhWpFEH;4C5NH#iDoF8txs*bo7%5M7 zO#N0PFAc!ft;+fM1>0PgMv}uFBako|czYt2lO)KGI52f^IsV?XY=BT7SlYbA6b`t2 zFk1fQ$SIpdz};MU$fU>vfX(43DwHpvrO0yj6dj>v`UEDIM@g1NH`Uf{6$}9@DJY&4 z!yaO!F#I+n=(U9Rf5J)l{{2!tYP1Z^eVrFbOj1(0lvx3NryfGS$cLOA%E!2NlLm>b zKiy$2aWEdBsI5n}K`71h+l(WDW<^kIA4b{=D7GO15z053&KCcK#Z<>|AZiP_$}ZF< zwunGAfYsZ}#7c*=1i^L>g=2s&M_xEC#X1`hbYS*E6ry@eT55bU)j$qHW=3{ z>GrEAo3w0yc-m|%R!qSae+bcvva2$TLmfsX{sV{T2Hvuh|9EE)nXUpozA`}>cEICL z$9+Y^E6ICvdVTNC8!3ki9|pxx`_X4Ql<7O|e~(&zwnve{7x*a|X?2#4V8J1}a+ZzA zgW1P9>W{my{7zAxCa5e0=SIqMdwYw=3nU!W|GAu;;b>(WDAG&b^a&knnBA(xkB8tM zxvnj>ZJ7JK0eZ3nN~s)8R-Dpo%?^JpB(*Ocf;v=dQ8iWEm_^>@J`s?d_TdPD+O%dsGn zk5N3GxP!c(@rxh#AG`D;RX_Wt7sqv#rD*H^R}@j!`i=g|GiN*_iE?fL!rg`ffeQ1&WDVBS)))*SO(@YP-Q6>Nh{6)pbfzW5Hly$d zda4&b#*lY0AHpQ z#_g>&C*kclU{iozgTa4^l?M;d?lQ_LdW`QQa(wU3GdoDI&evF-1xyE z*{CaF{=5m5U|?1_6M!jN$+Z?pR_N3R8M}^15I%5e6u09&%vY-bHb1CB?6L zAQ;9;E5U_4wj#b#D6UEkHF?uG;(sAfrSq-gDw4QoPI_X%)XUY6(8lR2{i7gfY|S`c zc-B5-av@Ugw15a`k~@VV&UfI1a_*`ge);{cUhwPkT{V!a8c`!O-hlh#K(i@-#w8nt zuqA!LP_YdZ3-guLt0YHrrD)a9$UdRw)#VjbA}8O;uQOK!cM{lb4Xl2li?_nzCoBwm zR9Wi?CQ|1NQr!~hl&bt!0ac3 zO_)V?#7E&(FaH#Eqh7@3VZq}#zUnb-wjS4uur*-5dR4sDf>noFi;v-dQiAi%QcGs# z9fCOCyM<5Uugfeg!u2Mn_zS1>L22|YX(CSo->z0h^FwKmNi&0W{dx6uthnPc;CHRo z8PG&L`x~q$6_It0+9ZA1Hd$I4B<(2>+C|q>ptX87j0IpE$Kv~{V<{MEB?3Z%K=U5% zu|?8Jkw?}kuIG)s;VNo~XNCPejy z2jnugv)$lg=*&2kT9;rD25F=k+Q5Q;M@@<$e5;}m^#>K4X+_wc#!eAJvRw>Y_=Ezn zNTgv-2#nj8Ejv*d7bC|u?yuzDXspTIWL}sdU=gDCWDMtc)oCRJSOsBi?j55}WAaT}f7b3(xygB(L`biqH>4VopEHbTS$ znIA#jn@IQX+r?fz6)jnKij?YaT-{t!@pDckyo>|X!I&A1mZH1f=Sv+iW7eevNff(# zh#%TDw3<1(;YD%|n)0O_*JgL|(_>2k%2IomL!D#7Mpuh5T(v6u<>~mt5}ept&Zdh72+peCapKWi~ia~#kKw_BdU434__3dms|o7)W|d$x!hVKO-~Wd$$`i9nqz;+xg78RW-qcM z1>rP{l7bOnpg)4z-a56{Wq-#n4##Y@32cdFpdrJ;6JviNm`Y_*Rt{Qrr}JP**Xg-L z-UC&I;v_z6l`q7bX$m|oRMDf>ubj}CPE3v`f423Gy39eB*7Q*4syEV|9?FB*LMn)246}EX7UwW=42JPR3FeSw zRjOuM(ejW;3Ol}eeHFps?6S!)?1kW%dly+t2fmzC2!(IYe72Tb{6YxBN_xXL%oW4y zIu)%B>6gN=HX;l*ZkIVR9Xv)I~J?CaDyqUDuX zKGHbEDiEc5#9cf+DE$0%LV_F@AxyaqO(j?yB&g}9x|N5qkStSdWUIf*80b(A4|iK; z9&>1?;1MnP{tR45S0rh+nW-u)z6Zbk!;^R`KG!>G+;zztt^`_uglo2l^k|tjr{AA; z6O_PkY#zhy4KKvc<_@1c*VAjxQFYNpLVVfSGef($O4)iC54k+^(V|rkm7^b#NA8Bk1b@`ETp--GwF;wx2RiTodiyEvCbYdn3cV+fI z1E6?Q!9kfaVaxb!x^X)=l(<|+!td|?1FRxS62pqy6{a*XKaxwi`9lJz_Y*%!%q;AE zu>X;$GArJHuK8`P?#IHuu`()29VbWpJy1BI%k?V_tFxZJfom(; z+w=J3*$cor+V~V*3OgaIc(|84Ag9`CTC2k4F=Z;1Ba(sP^*S0+2DA{#%Yv)vc!+u& zK|eb*Kig$kTa}KA2b78UMvCR+mm~n!?pcw^mr-t7B}A%}3_a)t$IjrV-t-yL))Q$Q zb-q(ZQ7RaLt>1^Lc@Q3kqnZ9?r^#MAU~@Y1*{KxhKN0UMG#{JO6-R2bHRmd=4!4_} z27;+$D=+i;LH)e#xUTRVB|WCzBiczsGI)M?+%&_D{hxWjd+{<@gB_gsk9T% zaH3jV5ozjO&+qpeY#-WYI&|)+rdX1w=E^~IfkABNqZ*u(B)U6+DkXh0q=T#ZXPtO@ znMf5L2mC1$JWCo|IlKAFqFf?UYSaw0w^0&&D5kOTFk5V^a#vV0?6zM&2D+^?1Xdfe zE<%x8hjT-WSoBPl4Zz34d7@1YiC2zPj*=Wtz-!+kmc?RM9`CFO89*iQkK(JvG6RQ& zEtLgosUz@YTrtKiX_O{3`jPtACF)_JC zjis>94U$=5XXsK=9bsSf8mQ3NsDnX<6MQofoN-LH?d~7gG`l#cFh*Gmo!^FD$9jwZ z!XX@bsI@;u|9wIZGuv3Jd=g2w78^TxvOLtRiaCNLST=%L$l{JREtbjX-RF3YBPExJ z_{TwiX3&i#Q;7sTnNBEQD(ZN$%RYl_*SV6pUbMJlWI6N&1s*(qk!z33X79E1y6*u> zx`t*X5T$7{kw1WStD0{tN(H0&*BH*DG1el&(IeDN@XaDY8Dm?m-6ZI2H8Qx!#yESSs zO9}7ARJFp~3>~_8edxwhDhfGPGC0~W6MItog$$D`IBb1cn(b5Nu0m%2lZi{B#XdH7 z_1sCck2{ym(SkksmbEOpttGyFg;EzsVdB3V>omo;d4gKxITj4Zou?ng28s?OaKV!p zu$Dt1oth*)mXP)R6Axp2JVoE zF#}}5PRYmb^YMPz@HF;q2)&(8XSWM3Stn1{<(Jk-VX9=yo&yBxiZ^nXNM9Rkm!23% zKuwgg{)L$a<8lrvA_7o06s>Xw%loN>T$wkLZY;M0Dre2c3V6J#h^@XZTLID+tahyW zSA?;XO92y171}pN+@x`Gv*Ib!v59}w3bVowWDv^nS711lWP$_4Rd8xPbD1fE?1@kC zTCPLyFWwQdh0;WYv0wSs8~VjdE9Hte%3$Be)n-JT_s4=-*=lbcl2%+B1C5fQX--MVW-dP09*MP%0DuyNhTNy{eK#Jyb}hxqb_a z#WPM)`cjRiaUw^dNn}-K3d>lc(Z9eOzz;{!r&D-!fN6ym6R3o?*vLlDV@MB{6g(w4 zK#pdf zT7va*O7Bo2k#6Q?*XHy58QyQt31DGrP?onr9=ik*X!6p7qy9Qh7hNA8aR{K~rWi-lRvd5$$sKpJ*3Y}NXRB=VtSs6{3aDoOi|`Qsqd z8g_}sWMw?Kr!XqBNU8(bUUn&7zro2tX&a=CB%-*L!3y;4+uLxO5oafhKc<;7bUB2A z&w}Jh9!5cYUVt36FPVk0A<$PnASthFMA$LNoLNsDZxp|Tu$PJPN6^@zhqGaFp#gy_ zd|OG%^`jd9Oo3P|^KsKLP{;ToqmD5mhyFX$>u9dEkR zQ7&AOVl7dI#yIm}6w;QWe()>a)!I@rk}|bJYT6?Ou`#WUQZH;T&w#H<%u{t{>G<=+ zjVK|X2`(!!Hnx^q3?6NEeS4q*5~tOJ6(MyBT{i0V53x3j6Iw*<@w1Gr@}aq?DQisy zn9;aiU9`1bS#qLB|CgJyGqGYziBFt6B|H(za@zVg)_@E1Y#a8`UZ}4hP_`j#*Pj#F zRn^BU0=yBLPQwH}UYFkE{2VVM9y5?}FnN*Q{ydK4QKv3c*}`^A&>rAV{7C)7r`uQb z4Fk)T#Jc42Y~4#F8(IUhvW&CrP86blAZR14DiTD&p4DjwKDw@I>#KXUPJP+l7Uqjo!&>9)e<3R-yUa6~gm{Ia0RC+pUIy~Boj<+N)%C-Y*ghrY+{_!I4_1krfvSS= zvcKqUN`)maIwW580um=8&*KaOOFNJv(M)?v$ghuXi4qWq^smVYOWi54udV6kZC@}Q zr#_j>5sx#grha1JH`^=96oH9MvHV21*;W1QFjna^oNWhIL;Ox_i(kyNCRIt+zZ;7B z7S4-5Zew|b5rn%%8^c-ZW*S<9Oh52EZ%i2Zv%gMLnqd8OrZ&%YI#9D#W@nY{Y^ZaP zzix#{E45qgLSmRWuVAbQ2G~o2y6yJCWH!NlO=<_5zbX`y$_k+XC)iI$1{4$7U^ZG} zzuu$tSYj@U@3`>Jw>S|^O?W3u>PuAi%pY2BoXT+O`VByWjFd-6mE>eOOGL#cu&P|T zE9T~+SL<@P5_~E>?yH2TsEP6ryHCW@c){~8-1FxJO5 z0iqYlzQlf|qWi`4rSMGzPj{D5rY*3z2y~T}(p{YN=D`RL>_=!vBZ<=WRVW={W$5J|wVi`dSR^%-OB1pRJFN-d zjj}1Jpw2$1g!P2$zTRn(l*|QKzp)(Sn=FIB6~TS7QF=gxNr*F(OXxn)%I#EE73yQA zvkluL2OPvWt;o*fc&Yg&X6xcp2~|)LUbo{9wjrIyyI_9Ggy!WmYNMA8T-~=;o@PT} ze}`4=89yzQqbA66OFG{G3Hg-7pXQ2mSbAof!G%LdzIx;ACUEMPRhvt^Dzd@%J@lY` zD>aZwKdZQmN`sS)ZSjGd)osG(oh|(c%}#eRt0j3iZNrW`7R-P~xNITQSA}g^#@qXb zdZ09coGSKvp}qLDrvyZkD~Vo7;HizGH)D#09_1R!BB+!k*dOO$9=_{ zOcD5TswjPjFxMX?>O&_*-%2lDT1kQ$B6*c2`yb|D3H)N>!xK(tf>NW75~u;%uTJNe z$;UbaQUSpoWa6Dx$3y2&gRRpQLpI1QD93{NGLBoA4cD!CZ;Rr)M9K8f{2IxtYQL#1 z)k5@|NP|j~$hc1LGD<|vcv~IFbz$W2LL5u+c;(0ihNdF2z2`Iu+%kmbDJq&!RN!yg z{jD4rajG7St(ReXGmK5yCF2=>%zBibE!=5@bfN8Fv13C0DsbaID3e&W=alLw%D!nO zYsI*1gS(sb=9-0FE?bCp^?u`~y`d2K#2{LnfKa>i_&M|`JgtS=mzMm=rw+EnYn@)r zN2i$2&L!w2kW%4`d%*TE_RGE15$=foC7t84b$hFuq~yBhquQRI$N7POEK*cxJ^<+N^G zo{#uC8$+6kal?EF9JHM~7SD!#e8Fwv2`^7EDLvOi4`Na4iNfGQ>Ucpw)iFjsL->{T zAtCCntNzpN?oVYzbjI)zvUb>+hJnwZns?O3IlC(+I@u{-~7Z(HO*b0Sr!dMsA zp8={1)A`<)4Qfyjp1$Y{twLkG7V+K|;FIcN(OiZ~dnAwdTp8ROUbT4zykw0k&=4c- zMwFH9SDJ5wZAl7r84dy*_`|#vb(i2cF6`HWr!OoG%8}LG8gc9HTHv}Wdx^~9f!UHx zV6liWlD&u3sPTEk8UyKz^O?SA-j6rN+}4I8U&Y_pRh%5j$j-VpL*71L z=kg<>VHiPYWw@=0X>9EF-1@bo08XZ;BBdt%qB(T~wAg<_y6(jiI%MtAH}K>3sq9qtHe?anz^$r)3$@G8dzyQ~W*f8XYtr zyax>eSPeNZl1E760hjfYHRM_c7d$MQr{SDQ6sFREK__Th6Zf<4#Mxt=KGcS1mQj{k zVYEWCsv2VxhI@K%>bo57hrL7{|#}1Q30I(ju zdSr)Da<`JHrp8-2F{Fj6nLqUnWmizlEb|y3f80x@$zcYA>*TPLg%5Rq(a~3(ODwJyf>~>IKxIQKP#pDVEpuA&*(*j%zXe z=`!qXmwb;2PazixMVNd;E=qivrcXHzZ;IlZSet}A?{@MXp_<=NiK)hbpL4X&WgCxQ z-%qzZh1%#N<-=hLf_D*2Lc+bVKe3E?frlKzMK4C&z?7kfUt-wM3Kj)1e4oGZTjAS$ zKEoxyTchY*{tFMi=qbdFiaB1nZ=cfy20ovTfslJue+-B!6O@yvKm-sCV;|$$-atAA zyVZG`JQ;7Z@8O<8Fl@psmEy%LP$t&=($VR&aEF*kfbb*j5Tj4)$Wp3gmY)R>oA0-v zRAGISNye27c7Qz$Ntpri?$Ss-k2(UVd|T$gY%}<3hkHP>8!R+6K(JoY^?pvTeOC)ot`GM0=w0t&?u3; z22(4N0W{#o(T*`8~cpp-k9=JjlJuE;MD zR^5$Qi2kJsJWM2TW`(o(&CE-317uCZR@YV9Srb~yqT+la4mbyOO)NDlYZbZ=M&-Q#Vz(vHeh=x6j3j<$@<8h zu9+|9V*{9@52t)gp@!<+_O!^d=9&p@^{+20%W<>tbf&>JA!;ph!u~>t2%n&5S82~? z|5Qi8g^G-*G~@Xpza^*8mBS&s40ec(_q3xbaIAhHG&eGTp=U=mMw*IRyC-Nz`*iH^ z0xEnpQxdZ6#;|@RlK`T^WpFb-6SuNp^F-BZCJmayN2;;=Bd#J)xf5oCUutlHe*vT^ z+Z1`KWJLft2BDnXIt+x-qddf$MWB+9UAZD$k?!6PNq$4HytyIVA-m+Y{;bvLrB*`( zT!xw2!eWzMoV7{uwGDP0^V#jZ61N`tlYF-KdJeX-N;APLiLSttPg7@8$3aQL42It_ zDN;z&*1p6oCuW*c>Pw!Nq9NZo>L*dB#}sIHMnZWYrI{(xaa9bZ8*UlSZ1fP1+UqV& zgK6iuY}o9ZZ6Bq(tRFOb#GBtv*$$k&uu+y-Lh22Q3JP6o;HxnS>eCY&1Q~{+o0Tiq zF6?_D0wxYz|3g@I`?nQJi})$xw9DBlDSZP8=q%)j@Pnv6#RdjM)UB)~ADlHN^Mx3u zZYZhnT=G5g?GqrYoWD-HaJEP-+|^V%3sVhns|U3T+i>sBPfr-WF^&PY&X^NR;$30i zd=uR6+9bwO$HIidbn&Zd6mB&3>{j{%afj&#en@5}v?-O==SfFm{WVdtqggn{=A-E0 zOIi=Lai!s$OJ}P-xNG~#+2hQ_3Y&mu@LiCR(!S_+3~S(b{`(?Aq3UQU{1{RCg7!b| zh>9*$&uFy$dAGUKMAn4?w;xVI_n|Repx_u9J&5!i>vg-R@XeyQWNVSJXvS*hq%R^rCo~(ZhJH=F1^dbuC*kGA|Wx4K5aF@FpFINTWSyiT_ktZW;u8_8 zqV7sG6D&|b7u>|!e0Py`g90+SRj5N;e*r)+HwnZewMJ~>@P96=GGmqQ7Nh-H@##x7&hCwDm8o5L-25L@dYObg zFotSBE1g2v${Z`JsB+3L^Qt%#W&t4U!DN+uIs3l+I08c>YY^`vz(vy%9V5W9wB1>q$&Wre6PK;a|Z z{r8P}?Ltw~#F||RE?>>`v?9NILOqID`$(+(spC&75x@MT{A9BZCTRamqJGazeeeT@ z8@G~LY+go)%WbS5y^+PIGK>|)1#vYcy6sa!rz*?GXP;mVqpJIYkvV&Rnz^C z!*q8^a{x(+QyjWexQn~~b>5}dcBt$v|K>-Pok`|;9P$UEek(B!HgZe!8t$crv z*L%HqjeGXmYu2n;v-Ygd%-(y>cg*a!ln^E;Wrp-~g~Zh58>D0=A>~fiVc2nD+3hoB zS-7R#S6D_`_?T=diYQubw!U*J71yw$T2YVTf*COJNWaEH#ZtSIm$N3XlN4b*!W^v6U0rQ@rSJ$7eP3iSWR>AV zvzFtbj~}L&@-wRbrD?=ULE5-qM5VMhxyi1p-E_Lis-}`~nYq3b3T`wXUW%F{`=M=v zJ6cF|g3#EzA)K`^@+-L>7L5_8Pj9A#*_2!NHy??lj}3=ur`FVwvKSCA;x8>Z%6}&x zCmZbuTxL5UBUg-%L@B38M|orxqIt7vGx}>=-YVpy&WlhcDF~#wk~Ct`-j7Y~>2nJv zhaFmq54Qzt51GF176m;l)@gY5}P5X}Rqsv{BMi496Ckn-T${=Fd&}jAh zXrD@D+MkFk>RzBdKpbUrz*4! zdNb|$!r}S^{_rjad!g9E1o|O8PWH5DNL(xS;F;6z^C0kNJPacTJhD>!dmOYA(q&}U zS1AI5IfYo8`I6`eAL{mmnx4f6*G?K;Jy7pg=YcImFi9lxAj8~d&%P;vWj-`2IJ1wg zZ7MyItjlQAeCt82BzAOj+jb~_Q?MG5@Z^iELLa-x-S5Jc^<5ITWve%M;O^Vk{5IifC0FHYMMW~6X}EfS>l*L|9Jd=sF!|W*YafEvm~r`yCV%y3#eyXV$gS=Jj{&Cg&PA zI$|%y*la!xOfWBxx=YV<N>YM&j2v)^@D zQ4fi_wPNj?F&#Ob_UMAhX7+uy8qwJAOT1~lW9(FHV}j~6cd0$}E;cF8^p&b(%4n!u zCSQtHXIINaCsmo@dc>HStwj3Av@$h?26OWTj0J|kWYMlm86~dkb7OB42HaWu$zj2H zSmO+10?9U?gn|bgTDE!Agsdf2w?Xl&wV4B`3h$q5l$q52 zU(3LW87c^^VUbc?L@DWJK^Qg`lT_(-4fZ@jc+wjfi?T)BR^`NBd{RcxXjLYH$(n00 zsg^zj<@W@pds=-;zf^lI&q;~@OgfjAxS<$rF)R_R4^5Wn*9q&}a>I@#{B-v_8{$C| zmdlMjxCO37ZP5PLDVx*;a?R(6?n5wJe?0S1BP{|N`gV&+?kV9^Ip&cbMqap8PUJ_a zCDa85&dPj2Pdp8e(QkWPS5<1-^OU7~usWJtA=rM8&_V+1S-1BDld2!aOJkuD6;XY8 zlCRp3?aTG)G2TmbMxMs>J~{6vs84*uVK60{p*As{*~`rm3CrtJJgnoZG&a6m#1Qc?=mqVhEGtHSa_xztvi1B)KcY;dtFyM0|?lEKo3G zJM%LK-RHrG(^#p=V$e`UMPUHPjF_)=^o=#U$vx|5jjSnhN7qE7>z&}xqL|{xWL;=9TcbukWl#o$Pcu)-ee6 z$a&vd%f?J@BUC0NI6-Xi;i0tY2Hn4}=AvmP?in4boQ1Y173qL>rzes+uj2}FpB1id zHt?A~@wD0p@>lEuEHp2bz^NgWkA-$*v7cB^4F)wFNY7J{eReHmg2DXQ^3?=0Y`b16 ziY%B_AJD`ky)IZl~?`-1mL8?`Q1oRaWr2jF7vC&lpfj~y3WrGJz{!xEYhxN9&a88`ZnBa;}_0@ z6pm4iKk`(V$@AW&BnEc82s6Az-JKhyu@unqsIrqo&sUNqenV4jxKTRRIcX97agYs8 zit-zgV?u|Xmb{gayt00Ee%0z$>XZUkIY{GxLr~OaL3O21;0o?1#wmW;<10ZHB&?Wv zWTOJ=))^;xppyC5&@$4y1C!z+?i;*zBjYc{n0+6$5Oqu6NE<}4&LK?-8X?lfNQrK| z35i`VGKknU-}9UY;y=AGpH>7FhdwO&{x!>Y;SZck|?U(hX4>_V_ z_^7V;YWTh9RD|g(v%|V$r#ftuw~FNfrzrgLYa%ZV9=R!~;)OsmD|=P7n5CC$y%J(z z$0fDv9p9hCWCXqMGN;?NL%pYEM*bF&Gp@0m8l@Jf;IKu*Ys2u?u92=u9s&}WFGqDv zp#0o5Cp+U5T*e>~zk1K4@>ShSN>s15Bd4vUr1%Yb()1@>?vDP)XHHMhD6ZzITp1Q8 z7>wXw_Ex%KtyL1?b>&&$VVENLzCGjC9LZ8sHET0*X@^L1HEX6F`_g`i8LNXrZ1A0w)FE$J=GokYUx{wjioFrKI<}Ft| zQ5V%B8Ieq~lMC(o&4*NbXbWT|p{n?9m&$)uZ*UDm0-`S&C5vB%nKrMzKPNk_chE46`zsEWxQ*0bXLtW3M(u#XiuS+V^h{ zW2JtzH)ki}U}WZLT#y}lr(`(m@R2Qr9*1?Hf>Uv=t_}^BazqbyPAZ+TwqQzYTuQ%r1@~YM^f6}IO9zzswdukT*@LszK&w9T?!H~o|}fx zmBNyp0@Xoc&2F6Gu<^Tk@2wv51-t||LmWR)=~2yU%csS63WcFrg2-cJ3xl)0CX#TlkmHHS?z zNjf+Hrs>`fj@6`I8FicFxzw5yEs`*}O*9^kC^v%cy|Whn9uN_I zNUwlXTlCb<^$ml-BC++y`#yeHh3me>V7o7x9E8di8H|LZ4lOvCpO0^H4vr8lE*WlJ z9Y?w9VdAVQ)D>4|8WL_wg3BsmN}$rKuRX6J!68U_$Afxb85>WQDL)vtry>o;r#obm zlj%ux(?{X;?Q^UQ!L9dAH_&nj4H@mOG&OKvN>B}J+abTdjA+4@5T}+nyVoq0^2qkR zXbwlXvhc!|1c8rf@tJ3x`Sn`In;l6Rf@q@(IT$kWlwLWv0!6sU=L}0TbIo=vU%}p% z7H40pi>PCw&}kQUX4;u0Z++=qyvE76?Hzqnh{?xRyG>Y3dqgZd9Rzs!l7xS?Fw4 zSv!?}_BDCUK7#;Txl5Ej>-^Mu7&uEC2Bav6mt-h|$;D{jOKza4TNYmlv`+UPF+Y|& zF2nwYicreK34HerLG?X51||$I*XWCL;=oiYzt3dVA6d*%aYcoJvzM$~p~FI>FM~y< z$B(b|J)sIX>x51|cl1fguy>kN%W>nwMeb}rkq9juSiDXGdIhTGK8^ibs-j1v_gzWJ zsYRR*irHfeoJDK+OHr*F9EHKIAj~m2Gu|gCErk2MEyA=AmNxQ( zaTFUn@seO(9(hF6kilU~QKgBNpS&jFiU}pT3hwu!t?CuC__{CVj!n;gxyo6=SFvPL zMLcZe-4X?EdgFsrtla;jK*Jr2GN`KBN$tvDTH%w_(@ z_`MMLeM&BTDz_&{XTmhR+3cxv;hm_L4aUYMYt>o2A0)9lrUVJzb}06$KM3raOyD9J zHvbEh5+rQ?{QvW{E*Y}gt&tkFV*k5HtX63X^7Mm^qVa6Szx zl6Lw?x!Or|n{^>{VRDpD@_Rv7X%$Gctyrcln(s4w+VB0G^k@OQpycpqEu7(${6y^I z6HGsQaFN{_s%knIk{%^@)4D>*fM~H4^%IINq17;YYRZ}h+6{vSR!H}f7K#vE41}F_ zrzl!U0}IQvl{}H=-3I&>VU$%P^;{8ttuEDjDEi50lB`>$o~cf(1zm0S(??YnnM2bY z_r1}1J+COZrsKpHYfI(yOEvCGR2*SxY?4;_dfak{n;MK#OL=@z`%2O z<0?k2$`hkgm-IJqdFMiPsSsOlH*n;D!At9Vw*ncnU7k_po-ojDNE$nA$?e`-a>L8E zpUxWJ;a&9<|5jpym*~G|TGTtq&%Q=ZF{zqXfHoax<0H}ze51qQd;jLbE%=L*!hoj& zB4q8X*eNwuF+49CTAp%KONQS(&PAi844F=bu=AmuP6-L#c{Or+I2L zrq7&;3x8JYS;$``(H5G8wjO~_5G&ctR3U!*oJw3pQEOkhK8Yg`p%{&i~3IfkE_mNX2%*& zNUb?g8iT{$_pLMxoxbn$Kl>iV@Bj61$l3IG(VE(gfn;iC;&g9KcB&Oy2mNGkps#O! zI`r;dZI^SwjKNT=>rs}*z18(k_Bb;QC)bA--^i(Cb<6QV;Y@=6`@LNf}}a?Z9t z)}AR`_dfb4`C@C~-r;&$YwJnIUZ;lB*@wNpR^P)u)6*4+e6fsbJ@)50&Y#vsnU3P` zZQ)CPPhWCi8XIY=ZqV>i|ERH(^;)h|loVHDZhW*za!>eV;nKYw;glLxYD&86}$zE?QYJnbYON`iBxHHA8FIty_rATK#HTa&wuC>HG5AZg z`9kg!9Wn@UO^|9vwX8U|9xwOtpMGBVKiS)x_qVVRLiu#Iu`ijIV0QM+^z5s_*}<-$ z>DkPg?J*EC{&_6Izl)(K%)r(bk7!%`Hv!88O*Rfjby$Lppio4_`$yA5aUj^_N) zX^(SSPsCL%ljnC#*7nwozqd>%!!_ORtXN6h{j_mkqdD=^i!IJ+;^mp3wMgai$B72= zh2!@pZ1$EQSqNYEW=Hth^GmPN=`Jfi_S&32l_Vln-0eX1-G4gY8aOZEcko38TX<i|ATKcOs8wWUa?2vl1F<*@U_#v&!+cEqc@oD<)iOz z&jo(pu>If=R<+a7UF|oRbK{0lgSh~oAn-ZMp=sk8e4glRrh9(%gS>y?wcd@UxD*2M zs>dg3E4!m7jXP0kD{=eYcdFCkCRkx3hR#Eh4HN6{S5)pCZAJCPH2WP7sN7PD<4obw zU=rNgVH$5ao&GrbbulVO;$BmNb0t&j-0rub7vi9O8vn1`Lxuyu9nt@>?Uey^l=p2L zUWB)mP%y#s>Yg6tPP7A%7Zvv+4c3MHxA$*yL~nbgkg#&9+`vHVuu=t~o^3&xj0C17 zW*#J%`9+L%tzXKuw$jjG8dN^Z83M4Ax5O-K2U&u(O6gj7>zT)klaoz}*@#cf3T}74 z>-()9pAI^=o_Rg~zVUfqa(0wYiAvw*qvY(_Ot|T!Wnf$3)4u8By(5kJ%op<-0Heoa z@Rt#>mUm;}8aZFiwijm}`5uk06lK0x0b~q0f9enyf1$H-()#&LL)&VT%c?Ne^|eE1 zqqS{C@1y-2M<&fPlGn%OP(KWsYN zEpu);Z8beL^4d3*INICODQK|*_iwzDNT2%61)$i_$`Jy<_ zrZl^~IE~w@}e9m!iwO(Sn-5^QgC^qj+r^?a1%+hgm?7Kdca~8hI zQ2+WH8`E5yKT(+e+4GQc{xjnHFRv9N_mTd&_wmmirO!7?MsA)Ag$g2sf9*4UzQr<< zY78>@5AC0uPygIzSqK7bt^7Z>SXLWoB5VVm)DJ)85<^=!lh8_LnQDz8B+aYI@!-fd z3&(SEH6*ALrLQ;`(@FPRWh6_~bf1$mgt1MB6qtVR19zkw_S9Di8dSQqP|eauUC@LJ z_J0-z|8K>_WBUl(z&Ns-s_i(7I%r3_dc?h_6Jd9K+lh^Pp5amvFg=N+%WQB5``xqar)QRn8&=WV*?32!&@zf(8oeJi&KuhN4-LIg9&XwA59yAkj1 z772ZxNsM!)y8S{dCa1e8-e^m=faw{+@xzH>EjwPS z0#E2vcMB8SlK0JHIyCW1_AjGV6^JgQyS9$64EZK;h0h1T9o3oGt`g^CHW_!s^5a95 ze6>3gvyZDQG4(Iy;ksL@U1lxPFI+yHgkA3F?Qu1X%*Gp;)og@8*?oGp*+Dk1>;Ymq1_S2Gz z+paKYbsD{d$ZA|$wkndw)vd@vft63RfYnD8aBtsid~s!O_gnwv=oxx1edWS-VWKbA zNy0%^Rif+;FF%lcdnX>0kS^F7#T#@bO?blUdkvF3-CFSb?^$JBwUlhfvMFLU?`h9+ z=8q1)9iL)HbKYaUpm&hlqyOfRFtFfX+%t(PPUsN8AtgfjNZP(nJ^WU{*_i&$e};o? z?LPe!mHf@rM%8ONjPdZK5dXC;mzwNtQHLQLn+2o%w8)RiE5no)U-Iux^(&ZM2J>l&q(~d#>Td59I-pb!G z86VBDs}Si=%@nI|R)|dv?=u-qcraaHd*q(|zu6uNhyQQe%O|M0VKU={>_5!s-zr@) z5FaYgzNY%6m2jZ>@)e`dW%N-ra>Xq&*^S5)awh0q{#EGR2B_$k$W(@l#xtku9o8K% zbK?b3?mp2QvbXECCA0HNH))06@o5WqSX<@z=4Imd$469>Uh}w;<|#Fr{6%|_*11cJi&ep5z);HA}h!%tiwhD*H`;0DLD9Ud9PhR8e%x9$@WrLrV68ctN}c|TJ);GR>1!I`j&l! zYuvX>2cfs;&#n*wdmKslqo|ykU|nE5#EZ$FhSE1<#8Lz_?Kpm z&fM_5vkC3j6gX0B_0QNWELqX`AqT<{vj>R+kJ#_ox)~hvXoM=htXQ~85tz1k!z-wd zuFG1I2*%W$raxHFl|;}Zj}R9lg-)`Oh=7~0JST0YBU_jkYgx4i6jDbg$adavW4zX> zpI-|ib1eJrepHU{cr;C}?sezW)sN&6E2-se)*mR7It{VJwVta${Ng`3ajmf*}wFU`;R_BSAp{3$IX0|3332rMBT7Ip#F{%QzMv z1Fm67(dG+bSF(p+LGQ=14=u*70ey<$jmhIC<%w{>GN45%>bX>PfkMb1o%l^5gwWq8 zRHq3jBxo=!0Eg>=p-@0MsFo6B3AlH#VeZe(f`rXdBT?hcaW1#hx;-;9a&Hp>J<+TX zzD{%99>q1E%=Q9ZkiVDln=Vk;-vlWQ;Y18(fiWlp#ax zO~k?dYd1(HU#gB9P}%a$U|QN|y%O$qiZ&syiSeu`(lJy7=MjW>)ORF<+;t0fdSju> zt3nAJ52H5sW_Y}W?-8xzBVw2-R}Q#+gAl>_8$O$8v~JR8CUXaj@=dF3AUa&CaV?84 zSsz18UwC;iKIJ5Hm&ugoBWE61<1Qr0)GV)kl&h`fQX)~bTAt(o(eyg+1@^iC0z%=! z|GV?VHAQ3PmUgCwg(^+rm zD{(}Yvc!IZ#d1C5zLVok#q+rbVa$;zz64s-){g9boRj=o((!dNo`}r|1`pZu1uJqj zp}wZ8+R-%0&`*qoUMzm)Vm*&UR;|vhR74rOZ*UP?M$E`qRu?ogHGIh8wY5V8IKpDoj}5 zFD$67;?yaHGqXZ*J@?JlgCQSny}9Ni=UcMcCEU%CkS?&_>-gs%RZPfSLwg%;R2qvF z-sLkk%J7M3ii%;AsGm*y!d!3S!B->1pty|>&E}RX9w;5jJCI2Y#g13s%EN!@Cw%&W zKcIc0n(*18I8?g&@-ml9my;GvCYdy;`APz1ae-#ewU_f+U)J}bm^YoHALCon@8H&9 z$~Y>$Li;3zBDayOTM-{?dzrkIa?aL*-fU&$32G&cu3YM!jqZq&!Q$vUg=>_9&5qFT3puV+`|)MLhIz_Zy6CtvELn)Fr^QKdr7 zt5;Bo>=-WJe4c9B%w`2xyEzHz_pCbud&3F8>Phl6Jd%?em=e6K^ww{c+F(J1HW(s+9MRBl1Jj>h2e;FM`>2a~uz8})q^)3W4y-yl#f&8Q`Q0RyWbHvc zHZG4*QM6iUs4mdx0`LuiA^x|Qsl`bXp+j)cH+@2Lq(D;ML;8BeGiRnMlU}3vjJ*L4 zT3CT>;ojkNiz73{^YREv+Yec$+FTaQ=^8FMU>*7xuG`y|aRIK$H+hQeylU2fV5B86{=k^?`_Kr~o9pG0l zXwU_CRroJVEFaTO+6KbGIMbEJ#lxG{U{3;qehdlpuGtFFmk><+dLIQQ%oO$;l!n;K zmb031QZE+Wczb6tiCsASr-ya3Qb%vwCtTIznk^vnIU}-Vk1zHavCQ7bw>)AMyFkke zke%RPU=tX$0dRzBI4W5f6+7Ot4Hh0(*mbFvlE~*We3z66UcE-GB?zcv%LZBVq5iYTEeH87f>&fk+I+U`j zHFxNBE!j_z$Qp05+@0QG;k9E2HHX53iZ4X_ksu8Y9mUSips6ivju^j`vMVmx4tpai zU?PkY8S?a{D5-UVz^(zEoh7o34qa8G^jBub;giKYVfK8YVjQ zl%|q?^Zo7gaj$Mui$qFy)Tmjm;+0!`MiVq5?~-8MONB>agU`*hR8Sh5lt?M z2<5xa`AcqvUj`3lU3N`(TP@~`#xj8MsIa|>i7a`Ag6vt849CvrFS6_HQaM<*);8@D`vQZ?KDFy2?c$YqNCxwomc zvhXxJ;d}Tn?ov92Q58g{3paUV?raC_JE5Zox_&rgjWKArz-*!(qeS ztcG<~=@yrBBZ^L<&qlXJ`a8(Aep=%nRYC~;Cya7!OH%{s5SQ`nX_99=oW@^Xgwcc`hJjDNzLuc{=6I}po{${L^j|Nr; z!8mev3e7O?vF`z^gD>B;0$Y!@Yzs`W0!Q}+&ki~jPz!PTT8n}}R#=owp8JoTF z$oy@~CybcgLhWyx#K{jg%HM~rqoT>s+We^+{DSxl5%`MJCyy*Ki6DVjyXhm->X?q@Mbh|#H81wHJJV0f0)yNbREY#->kiu##E@3KiPfn`laql>sXhCrt; zYfX4pdWV1lf#B|mh{Pu$GyCC-eP4i%g^9dsy zjGd?uCY72HDTAECqP*(=7=*k0&*^LM6+)E7I>Ph|s=gqgL*aih#;IvFRRGIm>*AqQ zLcl6niw~c-cc=9kq96;SL{66CYHx>5qt$D3KV=1_SUzG$S`~*zF+mxf@cTg}w<13! zmt_?8>+q-ze|vYWa+!MJj-KGy{7%FTv;?gnJ@Y`Kc+qEO2UEA$9klRSAYq!Y0Q;AC z?$5{|7&OsOZK;#F3Dn;XaT^@*7gto^@*fgWnx8)f<8>qjw z7@3f%!JcU$8m56^);QzjSFmd0r-F*{^G8aj~6BbFerxmiy#iA>p}bF zeNWrJpOs*bE9S4E!RZ^R)L#2i*h#F_DC9K`-&*BSPLyk>N@UP1vWyVrR8m#4OjN38 z3}($VtTjo`(B-bxUaYh4eAzl>^tMMX#wF%Bat%xOSmV?6sg&kloiLH4wdcQcUPfU_mATry5d<=_JzIzSS@M!#1LrlL^=wek_f9R0X zOpi!JUwqWFMBdOlay(0rS7NV%Hp{*;G9H%q_ts)}@lV|*zOpwKV$J-alV;kSvfG?; z$*|g!x%7^=PC-Y><;f3GXA+(6ycoQ96@~*$mvX9N7Y-At?Dx|hM=O1+RYgk+s_7jx z+;6a(8reP-xZb5H9JKp=aEJinObn{P{l^$+%dagwe`vo3EHNdXgthXrDWW4k-5vn~U(6M>y z!!hKRKURXXS!1g*JS<%K@I~r1JC6tC%>^MtUFsU@Ee5c4W?@R)5(^`Pc6*8p3)~qp zQyndxOkbMQqwOzWMem)W^0?o1zo1(#fO>!P$OQ&N=m20Gj7C;i7>2j3ka`ZrB_dM# zFv($c6FVBSyuf>IT1_JNE4+8dRYF}9h<9ZSYc`)oBb z@8UjZ_G@?UQ`J%o@m47gkv8z+VQ7(y8)ooRvWyMwtcn!7DCYzk;51OsH3wfs3BGxN z&G;Aw(p-`1YO;W8qX%7{bhmfKiF9L+^ednr$-Q>aCi+ zWDMtuHo&hfW^~v#CQm6~sQ5rL1%B5D*0IQ2ti~K$wy+NR^6`T})ys&_T1h4yE?-gy zZsUyDozhXgUa+K03{WhpAtSUGYKv%lctfc3#y4_dqP7Uj>`Wd4OFQ-=-`L2u;fJ-4 z&AYl8jds5g>#p4mzn~*7fN}zVakHMUhoZ=ssIX#Mf@WrWw3$2O`+CMO`p?6O%Qfm? zy^k(5QugWj^pymaaFnAEq6^mkfwi(fWA$I)$n#hh&_N*!NM=|@*Z#9z0G|0%GVCux zm3*8E5e73bv&qVjZR=)QafPmhx}SXx$aW&P6c5!6d>)i1&UXtvo05r{rv7U~NouV8 zH3?l`ZGGb2h4AvOKIYk6o~@v(++oa`QR%b^yw*+)0W%1~Cr?sz9>-QkJp#w=)j3cy zn&um`cqU2)zvAYf{m{tJ^LV`JU4B%He*r24_kz}3fN5X=j{A$-Mc4@d-0LTJexLxS znawGXQo?a7VaXH)x~qqFvZ#@+S6_Ij?$ZsA!@87dGQp+W)}gO(d&|f^YlkTnH|sOM z-OPYiH3-;Z^k_eRxm2*^+8ZR<@wf=T`>?*I>SGT>9HSv6Yo1&iwDZ*}ZR+5cY8&Hw zsrGGg-XU36#(BF9F}3ODgg)Qj&-3AZf)(*1BlSKWgMK8{J?XExiZOnB)VhYCk^ymbn>&(fozHC!kvvCJM#y_HMysg(Sofa+-t^0 z*DT*i(0xws;(@OD8t!;}8}AOrhi~Q|ue>E7gImWWc41dc3y0-rXz6=5$R-gEK)RW~ z5k2gn!Bze)M$sC_M$VAjGfzaA=WojRC^`%A4sAJK(vcFI=zVR+jpu=s%9m~>?xU8d zt_H51RegQhGBLr#@HwFCg3h@BkU{?jk@0F>bwiiNX&Sn+6Jk>|Z<%%_eCccg=QTDK zVU(1I(fi!!91xE|^T&iTyN-DnBxoKp{N9xbN5ds4UC>oaO@TbGuf98sldkEw{Ym51 z1UBWw$I)D~tDTDw8ry~m>oJlIlPd8H>E80|uh_^hslGL5wI!0wfb!M3iwlBQ6b#eM zBg(;+pXTf-ZM<=aci zj3%yp9GCdZ>34P!=$EwY%eSSJ+D@E^cio>P&3mh?3{s7LxHHhZ{ZiMg^TC9Z6yojs z`R0fZN6}R)D~HGEre?+w7ufIuWC{73MMJgtt3{3ghJ5j z6?&({^4o~Juj7zer7^X%fmlr94)woXOWGM~5S@!_tg|r{&u@#|!0D_;H3&Fq><(0U zd3a-VDUF-Wl01An@DqzRLNt%$sqThb>HXTd2UlDucaBG7zkJ4g$|^R0L9;Iag?}~r z3IGeGn!`_ZhXtSrJ+J`A6T$&iN^OGY3meyQ$S>z~VCvwS2Nt-KaJxT424QzC;N%vZv2QJ*1}(<$Gauyh(}X<5 zC+0p<@|ct(w57aE-ciM~RhmWV=*~tybxrkAt#sc`#b@VW#zx7xSbDq*6ubblK>p^X z3GujgV5SPZ%5sQzCRH%Wm3LR8=w|qrvnp1VCVXEs%jnSIVY<@Hjxe(aA3EQ#>E1e^ z=#*j;?8z`AdR*{^K{3JCLWuf0V;x=W`H_Q|_!o=O;?bf+ zfZ36*u4s4fGPX=a;5+}vTy;?=JC|oVN6d^8&>Gr$mn2F?-G(lU-X(d!yktOU^tAHz zn@obJ%;;9a`8l--{6e8iZ5QZt0Stgb{^B`|fR->3uzIzzPkAbJIhx;FBcF5Nv~?hQ zM4z$!dRSdDbH=uIQZe&K9SPqIzO2UCiT<1LmDD<3_kd^h&HNRKZILz+!R_|x&rmbP zP-Q1ZinU}fF%FMC*5Rz4#gMq);6yxU`7|)Co_Qq_WWhu32Hr7R2 zeiOY}^aQMqfNHzOMj8@_~qd^oh+^00h!464BTDJEv-RFDB@xcTyU_q0gQ511W?=5FciW(^Vsmeb8a{1$)!5I>*;vLX-!3^3UP2jJo@5WfS6-x0*`1n2~U zKmZ}GAbvMM13)uC6ku3DI-rpcpgHg!po9b-zk$7#haF%%!E@IAQc!dB3z(5{(L|%sh2?7BL0_Bl1g@N#k5c)ww9?u&C6U0LvkxC#Ng-i&*Igw=r zfQJ8&0q_F#p67wJ7O0@mc@hv5K0gWrPznLi5fCWQz#sL(fqDgifdeLtEC&Y(07oIf zSjgjff+G14=SDytk@cT<8q$D(4uC=NkaY?GWzLQAOF886+<%ZY!I39mq=}GfBb$ai z0_nUzkk$Q2cmhBB>_4J@WB>(!h=2ku28szFwE{Fp7Du|lukQMF1iB3Gmm=p9;DP`X z@K!HUeWV}!M?R7ZNj!HZKoau&KT?oNfx&p^*M15R2F96R3F!zwd;sYO|4!$&2K)l) z6aPr({($s}U+LUee)9v*S`$;6DEl3;#bk0U+_;+6DekHlWgfG9vyTj7X4z)KUmg@*mbhn)jzQ{%azn$^IZB zjsDASYRJ(PIRB^J&k4U6e!vSu4U3!BmHep0m7Qw$M-? zquXK?6)$MAYu~Xe5+ltB2vmF+Q^&OYIU%m&C7LRcxj9JRBU+XyU(lp1UhRVC5lG+v zx2ZJEZ2|nsuVQb7M6tk-iZm|JR9RMa z*js?SV0;jupLsxb9v&|4BK-W0|GCEJ>}Jb_#}A-uD^E-7e;0GHvH@9`Tiyby`6U8i zJ+Hvt*+~w0UvSEaK*3M|EI@=HFabCe!V89TfWaJG|8OX8H)|U_WH`kGK*`TP0CNZc zqnHioKQ!cBiTndO{iMMJgpm(Pf1^PWfPel@13>fdGy%i~G$CMQ_^m7g3>fcUG~ncq zI>1PW{EZI^xapr|q2S-^0A5W6;^x2WfC(V}Y#U4vfcd}i!QeoQ{*49|Mvm~vzyI0_ zE{GiU{zU`yLx!5)XfPn|{!Rm4&HaOhjO)Mg!J){x*}vt%1%Qu4Ujf9Q-5__J*YL1gg#cO3{}WUTx*O%U>@eu7YBMEi{o z4*o+I080O113@8VH2Zg1A@HAV6M_NQ_jf+PFn`blks=+&e;t^$l}v4bX09Y$|1~8%U8U4qOhkKvxpIf3o6||4@LQt|ZJd_Wuw> z|ECrEpH_y%;D6z;{{IR5hxYh3A3c#zgma`%^b~v|Dktrbq1Q)13Ytk^S0`0#@zC5*VG!_K{oxdO( zqA`NQnEu!I0N;<5wSbS~Imf@xx4V6wf6MFr|E4?kFYfJaFJ5cCf1tGcey9cn{7!#V z>i5U@64<0M_(Kx4{`}Iux83)u(^&0!l>DJx zY^AaIB=k4i@tjr7sCwD(`PyUetS+Fm|1av_sI6nA&+80FA*bWpKCh43d{(9Jdnewp zY|e9xwN6uALFr4{ezEzEMqlsGKaI+_j6SDMkI&ti-%qM?7;IlIm}tKiUdEfJJ9D8u z{3iTsbaku6<`#4oyXwpmIBoMfw{P6|Pi%Zj^B3>>9)^kMKQenuHQ79ljYR7bT2y6y zx{WIPiF?oXE^iNZ!uyKuNA&RBjyu0@*Ao_OFXJo;^rd}j3_mtoq59P>%`RtKO{s26 zAHNezDDh>eK91b#{9TtLG;kd3uK9=%N=)esu=SWu?`a>DNPi!1Jx7^V&2E;2ExV9; z@Njtk!S1O2QR{j8qLPz8rLFTo2eq@`rnsDA~8PZ z)0*q!epTG)%9;A%&R%Ao57X6_*2``W!j-mqqWcs1y)m=ftdG^K-$WgGo3+%+DYG5O z7*guY2Ka^Ddg4kZrz_saA5BCDOxnrWJBKBN%O0{`k07av-$EX%1%_^3O8h&nxtIF9 z=qm!4o#VNywXij?e%0^gt7|ykXS>Z8C_9SLKfAQ)cNsT8c{sg@kD~X>z)3o?TDNHd$hOCx&}KRw*4v zFRKr&zw1v@>CHBN-{TGZ5Cv5Ho}5u6j|-2Rvj5y_vPG0UCS5)o+4!5sH{V1Pjw)Gu zZSmXw)xk7i)bnp>CTrfhn>9PDv~-N?lc1N)Zf-r+U5S07;r*k><8ruFeU}cr_I1n!wU%X>$ z!AJG3dkvKZzovvWZMd3b(Cadkf|wZ}=$OHR&08mMzQCCe!n#>D*Z6WtV zSM!^bU)eS}wgFR~n1LRTS}k11wXMU)R%HR3*`mRHHmfH1yc1?h0wJ1TIEP*DirSk% zk4d6Gnx=n$K)RgwKr4#e%48V11-nu2Vy;HGS0R%!OB-poAVs(NtaiKA-!~b$>n#wE zR(fhl>x*cY-r`}B=V7vslXqyo`5V-7b?A-#bm026G`u`DHlo?~;NY#{T~NoqgnhZ+ zBKw);+EGpRpu4?O8e_-H%2z7(0H4POVPLskhsgHdd2@0)GK(m1|A? zss=hJY6!V)Zr;{mfoK_Xie!TO!=(PI9C0@}ufn!*FZ^?3v2LdCDXscem;ObYHKXZ5 zY_1*q#6!0 zK;`&pc~h^8Cuz~>(VWNfvMjS~;p8#t*2!OED0%@n+osOp(a@+5daPt3%x+n?7Kqj{ zOn3Qf`?T2fFh7txHL)4R7~Q;+XuWTi{23&Q%O?(iYnk$J;%eLdTxPXLh)|Q>X$`YA z7Kg8lRGZkMg25C^Nk)t(u5uBIPF7u9wZHj!p+I}0qp`JV+MvK4dH)dbw@aT-#3?>s z%b79wIr%{lcUFEp9Rt!1Xgzq}YoyuCoK|MLbg3a8#eJ%Vtw*gyAsyn@6Rjxh;lW<6 zexF;@)pxY|JH%BR|K~^3LxuNjF>NmO$-LPN61M7g^-?9?2}6{H6YycD+Zs~qF9y>C zWvQ2lK=$7VjiU|cJQBH8U5?+8qwY3S*b5BO|=$CL5NAXOoaz$uJW1{_?xb#B& zD(0~b`~t4eeoO!cUDc+z@;5iSyr2)TWJ)87^aWUS)NZd6V0COYWNX0|t~^x)jZ(4R z4<_bI>6i1M#@fpmq~cX9bhR+dTVHo-ji`UQ?KJTpQ^EY)yIi&kIF)rc7TlR~&JX2G zM6S(Yz0k;%na^1~4o{Tz= z@yBHcRvy&_SliM!4vb#a03n0)Js@_AOS!wcTfGMzMsmBXdJtV>jOk%it+V4zSN2S8 zx7rXjETPWz^f3*;mj9uy{jryT{e^6kt#_BVe@b>1el9&`I&peL<*{a4U04$cLV{~% zcG(a!nl{#Ej$f(HW74@(D4#zzeI_r#>vw zg$_xEmvhG%9WfgD-g)>(DTx}DR#r|v{u^>D6!h^D*8+t{w4z%P9MZZ9WN|e$rZ!w> zQRz4NyvD-_w@?gz$J<7Ablz(>+$ejAesK&P*NacNq~0Z2Mn}LkJWFc$TG4mG1}wI} zW-4Xb)m;<$v0kK!#mW*Roa1g7SKT1n1wD)3XIe71uG-l>!9xoDJ(?K16KtAH9_!E4 z7NEPst|IeOnQF|Cn`24R7qc#!BXk{@2+=M5IVt!l_jZRt5!x47XX_sDpUvK#^j+{( z9pB0g9Qh|hJ-;~xE$=1D*(Bhg7ST3Faa-N0R;$cYy`O_7ji$J4@fUR%0H+dQ zM^--`CjYyg(B}p(SMe83dBt7%u4_ZpLyY+H*)__aR?|&7e2&@?%5S+-*&;S|^A%}T zr=39iR$ZaxBgjwv(2_KL5QZ1uXEd7_^d`B)O{4P90sZ&W_eX`Qv@(o6SfEz#@X#*3 zxZGfSBH<<2`7hhYOjczB8LU%=iV@qrO{2|>!G&Vn%rtGLcf&sVq?vp*9i>B>!-q1{ ztO35By7b-^QQSWU`Up9S-Nj4gi}!$JjyDLMoYqIqw3-R%Z-2X??d@mO?_E>glX~HF z>6;k*tLu<~A6?4V6woUhvH^PM9Y(eS6X+8d-4k4R)G#?3jR_9C-mMQvwe(H9lvNyh znu=`wwg_~G&Q5HKd9PYIx#_|S{=Xq(FP(WE?-1P$!gZm1v0J`KQSs){ZAe zx)o+iJ=3bHpVRYKykoOcd&Tc}B95a}C-5Fx>Z%2Z%uK>kBVJhglaLAw!a1*w6~mmI z7i8~)J-w8#6EE@U71=O%;;Ph(KaP)^DGoM)qxyfVLb{06B_k+Xk~Ap4<7?i#I_k7n zTpUo{<#wqhqmzc^%nkk%`8*ad^s>QEh&{;^GPyE znjrB7l9b&MS?T)tchU1#cn?>yV5uSbr_Z1?V-PvqN5%tAGYR(st^SusQCqQ4TlnbeLzFk zT~T_SN(99>EfGwb_}njbevAo`2n`CZ{Bb#amh_g7wqB8SP?z*ruK4-n#VS%*=RNn! z+XW%6$)5vA@4w)&CyTs2d;HauO_e0d{j)X$r+OqnHdFlt5Hv8oXohLe9?ic+qX>d&iPCS~kMN%%LZtAShx&VBPwX zBfD$f!qY9mMS&0$x0wHOnF?rR0ISoCR>}=HmK*dZZ68W%lT;Gtp~SQ%RBL}FpW=k0 zz(_l9iwwjZv$f^mYNgIEtFwN`IU9X9inqy_m{I}paiA% z={Z@CmH6^LvQ@ZumcRmM-2y)p?Q73iQSRs3TN$F+W zR4GOL{Ug4_+E4DS5hWC%tW}F5H0rT(25@?nbDe82jUNW~FESjs)JF~|zCVV5M=3wt zyXE%OcA2Yz>_H9Rb=4e5>eLmP>9dv^H{o!-4pCqKoHXS_q_?>`|G>c!Ru@X{D_)*0>+{?5-^3VhgLFn9({3+wK5Qk#%7JBxOR{^~D>tGHWYoS-*lG$2{{(}j$Je;$ z3<9@H9i@j!9zunYjvu%Gq;^?X*a|!nh>oH~r)iVU!jG=zSpT+xxQ6v>@~8@9QT`pW zDHow75<>muo_d(|TGmNl$8l(QMi(L%1IeoiB=KvzqExl&fe`|dE6O*Wo1%(Fe zI&k&%UTk3w!#`12>qN3sFN-0lQ$g4z9r7G$oYC51oiGWl;S=G3X zHJNHIYf?z3Y7F41f(JV7gyb82v*JqFleV)F1Vxl0!!dfrL?Oa1K85;!JoN?ELYwYC$I3etZ`m$8 zx8ca$Zz`cOuP!gd3F*&0KS=h>=P!a^RKSREDOoVQ+vB*on`Bm})s#gJ_})RZYh;J$ z#$d~bwUxTvDMs7hAW1?K8{-!)Pa2N0&%Wb`_ z$k_lWxkMuiTdoyVVEg?LWFQW1i4m+d#kJZS8RjiWZm8(Jao(-;5F&1b{S&Gg+_ct9$RaKU)mOs ziAywT1R+kb#4Vv%wG;t;d_zc-0;5VIFkM$8ib_QM9fEQmvP_+VZ_n$eXg4Oq>AVnh zbD{;u$~sm3K3i&v+`Y>ctnwB`6hXqfb0T$9l?`Wi3CI}R2;Oj|Kql}@?COnR6m^;_ za<;lD#^D8vR3&(T{D8OO3}T{%)f3}imFayf))Y)R0N4^1fk?D)jyEhzBh9Co^2qF& z>%wkFSdnAHEO34c^bFen;o+y&QjKXgIY?{UMQZ|)SvfM=YiCjTY-P=%%%Pbi4e_tUiS# zAkR!kq4zcw>7Y+Bg;v%vvY9|q4ihavFC0VXM*7k7ZC*LigmxO%QhFMnp%MWni|l*v z_+8`pvvtZ;rbT;9df?A<;pF~zW5wV>=-k}rY$0d~YY`Mz_ON1YlHCt)KX0m)mk#Rc z671pzMw{Zqs8mstUKnH=o!Xq)bffN;-g<$e0Z_Ne(V)?I^hF{}QrTMrGKt-Lhg!aG zlgeC61;s12DpC;}Feupk>km8^BQth~EI*|{k%4bWz-biC!5m@$j9y^b$1F^ItZp`m zVP*%OAY435A(lGZmBka9>aZxy!LEL#*(4>kLEPh=OszH zt1*lxuEr7Qv4df>>o043+F1=8Tx21Kkc3I|R}*YCm{K}AZbrDS{L17&sEr9Xo`cdM zB`=fyU_Edf+M$DE$fm-M_5rCk=$a|#oZEF+t6#hMcGVENSv8aqt6A;IT84${+fDyo zIC$IRgi**fwmBYH#nV9O07SKCOhDY@Phm5uFHaZ%w};e5dEzIN1E$xj#ihHX=~J*W zC>D!X8YLT&=9l2@uEW!t#vx?aA+iCRq8)jYo%9K=Jyh0!{$Yj60JZ zzDF#N#72{mR>JDg?l*y4#83hU98|Tl^R#=rsT1XTk@@lNPasd!l{taGJ8{TRYK;Ch z#G5h$Y=&QQL6rlAM?X{hg7bR}B822lUd_JwdT1B^G|;YhcT&~PeRGS#sP^?+>&vCH zEx2tv(;mxl$x#D+a7L2)iIkCH6gp5!$hc35E{`KbZ!U(qUTDUHnM6}%LM-I#iIvI*rWY+S^l0Er{NBux;U_^4%EaM0Bu%f`7e&JOeh%{YYy7m#3W+Lwm z>y5QkyAQg9$u;5WI8H%}2n{7k6``>pdZ)Y zr6=@mq5L12RCL3nB=;d^shJr1-aP)a*RCT?1o84FFHWZXaQmy#Xiae<49=(d!E|51 zx85o4#pp7*%i4DTZM>YQONZLI9yT@5;Ny#hg(@o%I|b|{C{j#Dg{5}`#Z|vv`HyNP zfXe2*`)vIhSvovg?wb2WpOZc7#=~fpE08cnf?`=Yj^HU8&TsJ@brn<{*FauYbO1|%84?$7xk;szPmw@%h$EGH`yO9!fcCCu znoOgWQ?rI>nIw85Tb`ebyt6_cde1 zRUJz~_@}N2P8$y$R;5m56)uXfOG<#fZ8G(4T?xuy z^>=a?Gyf$iQtdHiWD!O7qB5QurUCD5vbUYZygHVI?5IoDT7a_s32NfMU;#7~kp;prof$o8D@-_8`apq`2Lojd* zxK~v;vG9QkH?kpwvHe?!0&fMN-#?tdM?&_bEQ%oWDb-^aqJ` zw~|IQofb|8$GN!PZ_^DaCpBd+vEL}er6<)-c@Gh{jELU{y|K_b*HzL+ucNzyF`Z>= zQb}rjX2g>pd>T;?b$veC>+;Xp;sitaFNn7Utk-ujS1|azt$;spx|5@Ng-aKQ>O#S}3|LJ&z4`$^jDK4Pk_X2OC{ZUvyeiH(tkb;9 zZ&e@`ZQldEiVej@oIgnIES{^Dl7>tnO*VM1L*6ijB0Fl@Sq-}GecpFRN3XIDjVZwM z6?MPw5PMPkl6tNIQOGM9qMaLm--NA?!-C|)&20yhU{j=#hy36Q%}v0EU%wenRMD)! zk-8^htAWm(2cYFM+JC3%e~>&oIed>3}r1=fD4?zlZ=6;!3jxCYRrsVKi)+AyLs z7?g?@n0NTDIqN7z=){+5Z!z-&08yz1F0Nme;=;nciVIzPi(J}^oM*cRm7CVTwcKw#FI*Jwz z_QX=kUe)eFsM-#Q`jMZq1$Q(UQ>aQWO;RI*Ss*^hWojBn_<3%+L2>a?Las~&mCfk^ zxL9T3Kws?ctfYxTS$!5r3XzB3!sEzmIT|nb@8D>NejJ4k!=tQ#4btzte^?>V%F0Wi zjLk{_3+IT^SoU;|sJ;W%UbphXzke+%$2xZ1c=wkB3iHp28aU}6?J$)2)|rUbSsuK% z=PZLEW#0w(w~oAJg34oM>D3H@TMw;5Kf_YDLLi`zVMEuKBG?0d3Y&arICusoFZ#e) zRuF-~QIR~cj11*#y?XjS5}xVkEEYu;DcL}zktxO)egWOWKH4E!svS(zGU97zVZ2gx zl&F{h24h(2+w}op?m|`u=8hq|5IEoST+U$FF+F%Z3UogF4G?o85fv+^gj&UAx|B(V zgwb<8qtH!pcMr-Ji2nq}zZEUbPDLpT%hK*UQn|*EG~Yd@BMQ#!{Bc)q{jgD4h9`j34>a|_WgHk@5j4_%C62uqHA zw_w~Q0o6prkewoKm=ko~NTLQPeJnkbnIf4Ts}`r0C6-9AU%|ZEW!U|o1B#p=xq>yE z*^nLvt}TDrXaAl6p~yOL0tC1_rA)k_SpQ+Ti#cWAx#L>cr+SCOoWa3)4tgodRT+LJeRl-USi@?ssoO#)pJj!B+xm< z_GYf42sn!k;S$Q{TG~w2qQ;ZoTC!-u^~8*AZVArawE1hQwjhIb-s2wj>hwCGk@ zOB}e3nlnM4Pj^lDi9$rg{dh$*)ZkMf)d-jONkC>W0ZYL<`SAd)@ebFzj^$ZFObpG4 zkrr1n$8jp=MbZF`=x?MnZ+JsvGV%G{Ak*3~XIwCEe0FJnjZxJ8pD=hQ$gsn7rHT-z z1e9lKNoxRt1!_(~@_Nn_aB?5vPa#kg)P*HZ9sKK6!!2KRT7vs{QrrV%ep%Sol*ysq zv(joC1f>_mn;d-|3bbh&+Yy>4_?;IJ??ZL}e(?I2X;4glFu}1M;g0cLBQ65T{wG^( zCQkA!DDNzr?vXZfKZ#L%7yC@B6mBFu!z2ov$ui0&f7NP9??5F=|St7qLYYfa%bi)BtDC1TD> z&=*fbh<ea(xZYu_3M&F@K)YRsv7t7&y;Dz@HlFX~Mt5|1|Yk;o_jdPO&2S+W2 z2H;ttcB|N#1@>o{{C3k+P-PJ;taSHV0&*OaR>uU#y(CEQ@laO(3F4wSNYg@~cHYwp z(a49XmHh)p_T08D$)seNBTl>uhO4#~xj4Zwhi9H~wIdmjg@Nsn5 zTgelLwU2X$?vMu#Fp$i`-rxaHHIZeS4^j*%rhVB4p9#y?Ej@#B%PV+;Os5Cv(7NLU zmB@fFjE0Z@n0y&;3xdiH^LInRiI>*%AslMdAphOX;K#^$i#;gR`dDD>q+GYwH8c_M zLV~B^xC!f~I!`5f2PLm}asrKlyUd-4AHFtDj`?hmquvK#2d#42j1Z2S@ZrcpMp+2r zs)Ia9s~-9o+BH+;J3lucp-DJ3OrQ^=zdrbvGXkO9bRknO+`1NDSTyk{l$dRoIb9 zb+Zui*QRYYLUYPbr~M!c$Y35zD{wXO3^D3O0qIC4k@FdQ&R;Fb_kT&*+O%^L=Nqu9 zbiD;uw7@1LBT{iQwdW4HCA_gxhmN16pDj4EFCp$YQbM5N@x%>1>r4uZXF4=gRqnY6L_I3^L$+@c@+?Z+v*A)BfqYvE zc|7gfrV2BGQ9e3$|3tueh_=Bo(hG0NA0@)?8M$p3c&}T0!&< zZ$eiqSD)<_TY;6yL{8?os-h%I{qfnTXMX@YcyYyYD^|SCJ!s!ORnSU(xQ-mETsX7d z9aMidnhj+KB`f;EW1y*58_DwxgDPZPmKwzgVU##{96&DdqFd>xEVYkgT%I1Uv2Hda zN!%BB32Vz@C_9qX?fi`?GTKEF&sq;xl~3TZ=H0J;$|hctSv<2(ASr5H@@kPblMyZ{ zC>{3jE90kBY~mFsUTrUaH5JzOsrgUpVSy+T5)Ow`Gi@Sz&_$nzJxz5#*G6mZiplL=Gk zQgu|FIb0SL9XBYGB?n%?E>DSaFam&$A8KPJwUa((Utnu~k3Ck7&5l(C)!N ziiw%+vJQC}`4p^JkTzc!;I*f)IIGNJz!6Dy4$)VFF0Y0GOJzL+7BdnhuNmzi04DSkgx>oVjERR=3+Ov<7t(94)|_* zr}$Oq0gu$|M-ZIE!Ur?vSmAkvg#pUvMa+=3ixYvs!~xFiDSw?=7|UQ=N$?^Cy6^LL zynPO8TQ+i5jHv6~AaoM}(C=w&f-AC3pND&JojEsw<3gFBbXuytH^w8e0q0D45(uu< ziKwCXjN``}JN`%29A+MewJ}RjY)KuZcXbJ@w zU_!d*S5uG6k*no~ft_qFYP`dM%-}S&nmdg|VxBtmW≺gHaAbe&R3@G3w zXWmn%g79pJR+Z4)T{hAfsH@BGi*x11s^DVxWlS~x*Gaj-cx)6PH@uw~jBWMp^St%G z@Jlt35vrqnMnaCtN=|e*=9bWYeRrTOtMUgsxE%X{-a5BD7XkWk2_sdW%Hn{$Ig1EN z-gNM-`?K0lDFk?_TX%ltXepT#VdD)RgOyAZzmvPCz)7QcoC}AX7v?BKRWOmVJrO6R zLG$%j5%6gXBtJVE^vI)YT+iC%6=ZkRzT?8~=ca6FtlBhauir|jhoVM32aLX~O9^!` zLph7uva}1Yc47H0>~5znKh)%K?;hA=fX+IS$r^qaZ69u(_eawZNl677I6CO4P8T?R z^MoiyfAKo6$Qg3JN%UUDByS7VaGAcAv@3D2M?1v8KI~~;Gl@sTn4fF@- zpG|n0bM7jxjJbC@td9#op^Rv%NLICov~(zfY$$!`>W)ouc3CxF!8?7h53?ChvJgrA z_~6z#t!?rV`p-Y@cNa&4I7sc)TxW&JY^u=$Zrpgh;ozh+1FSCOMNxkaPZ2#mSI?rH z8#C$GkakinVT!n&0z1>3qApZgBMWtWirW_>@-$h%*>kGui2eJS$L{JCkiQWK9+!5x zxPJTkv)X*FZ;a^BhlTZ}dxC@u_MHYHw4LCD=Ez3}sgT$UL*f^jTWY1esM%4BAa!uy zD(CMmWfu@OG2o$(DSWp%B#bl)X}9|)U-2*_Iy#5E=x;lyX>bh{lS`bj$e9r3F-ZRX zS}IZrQwnfw3*zA}X#@tqOnVaetj%L6OsdtiPX~xUKU}ihwCE;#R!vuwI}tD6aa)j- z_K|S1Ah={F6__-u*EfpGsJOv=S9q+mV;B)W&bY%R4WyxdFsBrYJGQ_o-W840uIsWOBO2vKtO zF{_k#3@WxNq~2EsE+2x@q2^OAY`>d??@Ni6toL|nt{V8PlMs#LYKazU+D!|xLslLh zRFU4;H(;%QMJ2$^hd(SMOXT`wuX<*p_pov-Jy4rC4E}~4gln>Y@~27bf?ps1eAnQ6 z-5@m?>Sr<5Qex8uJM<46 zb7dxu%SLSL-Ivk!Xa1U9eiJPyk0nA@(9))`7SuA@80~{z$5uB(GhyJ7G0U}Im$67z z+?|aeHOSbDD-r#7`HlnLf-@AIDCf~CaFWDP6t~=-(44`q*8OJQ*f}H8pDoZi?z(xk z)FUKC!x#WZ8*E^vJK7WCHBQkx4fZI~L#oe%kS__8moYKP7jf1TZoL6i4qx1Id+CkmMs&0=EU z#fz$>KIzl5^w-a5vUEkP@H%q%B z8_9eAzRvsF!6@R^0xEt3dsirIflyoYfdo}a+2^;?omB`oatS{vgcMb&aFRvW+>Z{> zAKaA$uzHBp8Zq`$+4R*q2YT%lDfpa+&QmFEpCQ;<}^aXz^?N^Ja4^WpGV?-)&I`vDdq&u+S z+0E}rRVeye1f%WsXMIqk;d}k)uK;eEYoJkB#eVQ!H0#u{<9+!U8PJ9nTYBw4orvGpSV}!=R-E=UjGKU%Aa;s*2TgDk zh)b>Z<^aq-pdxAvv}q1>U4&q(lo0(+%*`f{UWrl72YhymYy9wAIeBDY`;i&1EXZd4 z6aktdjvjJsgFGe3)K>Rz)M$p{k*>pP<`|qLR{$WCBZsH$I3PpggJ@c^IGxbsxn@GF zM&$jvzYN9)e||iYiy${{d(RY%Qycsqk-&s;pHHy#Te#X*EFaO2zf6%^vD+f~K}CKK zYE&JhBxda3V0r0`xBABOv8gjabO26;{cEIQ<(;LGQ*=J(QqDK1(F15s23?pQGA?qlf+P($tTZM)P@ae zo>$IVXg;uvca14)&4qnKA7 z!*phZ9MY0|(iWv?bmd#Za5mY$Wy`c`CvEG{V0%b-LL%M|gVAw3XI zjUeD`*qM%dnB+a^$DCSN9~;=MWY;rAgAS^ookvH+eKq(#}Y z+$<<(jpVK6#vdC7+5U{H_O+h2x=Wx82L?Wk=XYJ z6(tGN8%YDy9rcT_bAIt55z6TxuWkld{!|s93;mv^39cj)PD6LVFJ#9OzrKHeAnwBC zH^E##cF)um+en387PBJZ!2dw+y%z?Pgcp}o_IfQ2pL}h3A_@Ku#iMSoO#e9GW0XD& z&x;9v`LJE?owzz}g&sB*TFdJf)#%<{Bq2*_3kmBaa1Vfxg7?73O`UqTcePN?g;I}; z!^syg5|u{qWFthk%F`6rvH$(8Rb}?2(W4^aHmoPsw|~ZM_;u0uYUitqfF-2OLPVO4 zK;NRGg-+I%5obM2nLta6S5#y^rO<~o8Ba6P#OzlM zXj%f_16VrOFy+4{UZ<|jaXCU}R62sdcBsdCYQnQ)4!#?RaR26igri@|SSk1owH}VR zC=cIgZRQ)H-Tz*D@DcI!W02!0dw}xYBO|D`_>2!M`FBlQn}aZ^r_StXJ5C48w9IPL zOe8xQ;>eM7cH}seI^Q_8Xa--nLKzl}$P{t|g_@SCY@5@0NJVC?s%-PF`^q7d+7vjSLgnQB;@Le6TijMv2CNbw?6CH#MjxIzSJmo-&$cU>G> zNQi)Wea7D2Q83;pqk+uTxNpoPeQr#} z&V-NncPk405P$tBoqqVZ7Wv8EBrL%LiUKi$OcCZ8V`JEROhi0QX@rA>=c9#Lc!ya8 zayI(Iv4$TeMx$o!NMdDop^^nPQ^13p$a&-Bk%qg_j8kg#dfi&BUcLhHoK19BQ2dWM z(oPkNsd#&DkwI#xp8YkHK^GPN$B_a0=_ug|lEXr>%V&yQp^&~+?IHeGa8)Xz5o32A5BtVKC{SP$$jc}wFwEj#7ftrj)%R^{AW%Sw)6PU}m zNf)xHNAuye*6;mQ-cY}?-lTrg+F>$%n<&o0OhjJlGkt`Mejt`6pYIo69IH1!w0)E< zYm`YM9f}PJk0ZX+Pv5#@U5#Ko)wJLZQxFhaF0_kxDure&-rJ<;%n8&ETx1Q%+7~~o zW6h}DgOAgvd&H7-<8--?2)Z3%XK_Vd*Rn<*7{(Lj%Hk8NU|f*wB;=6@6P-ygy~*kk zzQ=85T+!W|w7@zQEIbjx6`j-R4xNRA=mWCnZ0VPqYud2o@F58-euzH$bkkxoN|XpC ztL$p3y|A6-pgo2uX`qka|KS!*1bg>?BMq3W*dFkwc zO+3rA^1_iPJB2XfEQnP(tn*a``xqfd`by)2*c({mctQBNq{a#(!CYr%@-W@F;A*#U zYu-5qdWc>%HtP9yzFBhngnie+oTW(G95qRr^J0p9fZEtKmM2&$O?ePFFn9$5*Q1zT zT3QN6GCmg9&@*7tL=7=4dGSZGw;*mE3s}`3eXJfl`8iam+=)ZkzYQIoFs1ets3I~6 z&0CK`;F=%P@WiXdUuj7i&oTqnK++w=?a}396NN&@Q@Pqw$Tk(-)dD}3<82IdU!phr zUay9_!<0WT&H)0-ZQ;K&+KQ8YvhBfag+xt52erbtLDIYD;p1Z-U+Bxf$HZ#zB$7oM zBmZe4@Oz2t%CFA~w`Aot**|&S#-rO1kR7UX zIKXyVfp~*t4sjMu9V0;;W+1+u4N$d$WMZ#4N#0l89{oHkh#!US(~Meq%#-A{Yy-s%PJ}AfKjYgr?Sm{65+6JX2F>COrP+5L|YX5>_`5Xz(Pj zGR(1SW#)5=d}`GaHJYYlSI7uX^`w}m&kG$7#Je>o^fV}jLf&m?ih`)3CpALKN}kx4 zED|MeRV;AbG&hxUDFvfVsM5NS1W#KwF;) z*5};q>@$JWtU>GErf;! z(VXoVq`Lq*(DMLBzv+p2QwPtant8I|^X4-hu-yBZsA1F%0qmCFxgx^np?;yti&%)w zV{a~!g4h_}YAQG)i*C=YJf4bm4)t{8w2voM#Z3+WP1wPS%elc?pM%X;=SUO5BbMVh z1#oxL+qFneu*Vs})g9tmCPrg}OdHPA0@IW4K%x*pv+{<)!B{jnO8xYl72DO8-~TqI z^knjCUU#Jm$6fm9ava^uvp`=M4k+3QWB*9(-ovA;jsA{Aatx1d$ZB-`Z$tsa`?g0H z*rPo>-C8yyG~>`4mb^%tSPd*y(SVzXBxrbpu3wB*c?T5pbo3mH?0gh-``yX6fSc8g z@F)75rh*q`^QOWXE32t(c>_Q+tVx$q@cf5;L>BgC z9QB@wcw@#LDO+XQ)A(?xW)-K6nn`Yo%#d8vp^O6MA@F+DY8{k<)ceC-TS)=1wR9BN z2yKy%cp@p6YJDiU2cGV%=@QU)dNn)n+^?=r2eLgxSFOL#cFoES;3gQF?l?KwdhdBGsq9&#pO;i~vs;isc6p^?p*#t~&oVy)m1tUT`)CyAu?mJXIQDd1T;bS}t}N5=CbudYDu%bRhy&*t^68BD;k>fvP$Y=i z;cA6Llc({SVf5p{I;(8x)8Jo41#8;t_~lE$ zRhe-eL4PRUS#A`>F2@=2lA!X_!xZJ`1UC=Jrow3MxKx@Km%#y=({6s&=VD)XH!EO-Q+H`Xiv#r9UUh_Mf(d?Xnk z*`x^n+m?*Uu!aL35=tthvRX`l*>h*8P=}A)Z*>Q^_-28QjP?^5?E5mvTNuIFB4l2V zS?YD+D?M&#kRf7%#JeFH1+?<}T#YM8fpw?KUAltg1e9^}@tLu5gbb#aN+#8i3WyUE z!6Y|^?<00X_+F|?Gp|xueO3}4RL7+2J#UhlyRB2UbWCLm8BRQF`rUQ{(Wvsme^WWc zPymvUv}wR17kmH4_=9B*8o&;FBxtNAQ-*=iLCKrgklM&bnJKEfU3iStWtnlDt9fRC zGmS`9bifpD>7`3*(2d}lvZpLA&ccq3QY#=v6asI-QV9-De|VH2?yV(T0UW-UlrQ6Sy0+T_TbP)FS2aJzUm}AV5!(l46O{6o@yg>NvyRIt< z%e=9hmvn9N@05gExmKuf3yMR_1vMBR(M7cOuj96uDy_ZiVzWAuDIfCsCPXL13B%;djXKN#aiNW1u}WKssY3kwrL= zFp{3q_GREt511G%vNr8NnNCvn+(VqQrIw%$ygFh9c}nY-N!n+jLJ)x#*Sr=h=+}wT zf!;k4e+v=;*(j@mvMFr+5MDHXt6qV`nZpKip6pQ~WNzsoEiXqXf$@0?67vi4vcLt1 z+I5ky7?`~l=d`+r({3$Ut| z?th%_PLVnYNW%#@bazQ9-Hp;EA`KGKok~h5DF_mhQX(ZNA|N0j9U`5-eeinUdzJV5 z`uv|iKI5J}dt%L6Yu5VA?AdeTg3+(zz7)A`JbSChGyJ75#{5-9-3)(OWn1FPhU9p; zz}?ejrC)R8losT+!g(rtWIA`U{to5SJ zoZi8GfXnf&^v({d2%kbw@0J{F=zwecc7R5eLueniVLk3QfuXQK`3cDiZ&B8$j>g)7 z02*_JaB4K4dzH*fTWlNImUqeV=%Tm>%$aKOyG)Y;T7`<*G<*-PSz0{Qvi~w_>h7$! zzBUo)wn?jKapNL?-X>8aqfY?r2w z>Xl3hd-ot*7dif~lLU`4_7_-7sPAA1Dw@Tqwdte2Co}9lzAZ(@WK7ag(?&BPb)e_y zBR)(_6|+jVwdSW=fmW&WDtC>;(JYw%G}X%eo|z(w%Z)T1#`|MY>9-}A=y$wkjXAd* z9*n~V@b)F|iiZco%vUjt_tp3tb#%b-`mU&GQYILq>*MYp>*9{^>y3Fl&e zdOrEu5z@I|4GAV${pgn*D$@MMetQI-5@Rd9*I8f2mzI9~=x~bziXEbqCyr^u_99BT zJ4AN@<(clsc)tUqU2XLM@JmXjv~3F7q-1|cgHn!sF&%0;s9cL|HWkZM|Ax>mDKkk> zP+nNbnj*O8X4El}tzv+hWEDz?7s1ZlV5sb}U$3i!3e zW*Y%b@vK)LNVpf4aT*Y(1y?1 zlZZ`tP}5zlRq=HPj4fA2HZUiX+uf)Y^id4C(Ve7KhJsv@bHCRX|XF1*18P^QCG)fq&NL?v)(>WpwKmfG-@Z=M0@qUR6-Y)L~Ei# zbzBhL+d9cRO2=&cQmyOPBYNZDPzqii6qU{z@IXRkpTW}Qffc_W)+^?O-82A$v1dQqgT`IsPGlwZy)CB*B_`!#GjhgC#zySnPn}} znE~A?8`)SKLp#MCS57iP;j>u}l)~6evXL`uSXkFeq}bA#IgK=ncCMVdl72d-9k(KJ zO0^sA4MShWW-FjymZ*^2k7?gvdz5F$o=<7ZmOqZK(i=KpSJOwEDIddLMh$&}H!JZm z*qWRCHB9s78R-+W1PD2Du~*u)rLV;4XcE$14(Ke(=f-WCQq-d%a3aYB9sJ5%G*4fW z)pp8pQ!F!X*s+X6XRV*n`_Z;i-Vp~PuF5qkRA~HJ(i7e+jPE;2ioTTkmcy?6J6=VE z{R*fdmGRW3&SnrBU@athpun*3^IFGdBE;RsVgO&kiUa}zR&qu}MnZQ_btJI9lwZh7e6G36ca%(>4}kysYC+FM3QPD&keZO^joW&-Hw zdmqlxH( zTWQ_B;h7_0cg)20+Bx0t_R49Gunx)CtkY-Ic?}dTcQrBzdv9S(bNMTl-UBN01k|$y zkn!}ihgT=U)rZxuq@eU(ha|=aRpK2yAIRWT7sHW8rx>6bausjXDiAS9C9lz(*(hu= z^p%W3?SLoPr`6@VI9SnY??Udz-q4x$%dnavf9vK;K@#Eo9cSf~YJj7QWvEL=AXAY8 zr(AO2J!?7m6pv}rSRzn0q7PRE6^AnkB}Lp?nHLJex>H7TXc0U^xKy3Dw9sI@KTeWG z9uM!$X;}`RJ_O6&5IUdI$yDxAuNufD zbiX$13tE|A8J73K%8L)cCL(%@XU-n>-6iv-vZ|!T18>S6WhPp>CJ|jq6gJ{QPAb+4 zBqQ7WILnJ`o+_NP9pM-Qx`C)A7Gfr8Z<-$=t$HZ ze89pOe6t$(vib!lq0xJ^>cspc%yHUJ?VXsfC!&JZbj?G@^{a`XiI%{u;kggWk)`SA zrrl?3Sa`_nL-$0(>o_$Q>2{*?21H+pYU@*qMiitUGCpF1v{oN|Q}A_(BXknPk)s>l zVxH<}eo12fseST=8EZ~~(Gi=2V+)R+?I+Zp_rir8@~$JOuPmq|AZMuan0DKiJI<$R z*%_3hxg)Q;Gwd?)>#69&!}muR^X2WHBonyM1mE2tZsKmCrbOx4CC1GW4^Y?WjM@6Y z`j|Ks#TKv20cvq;i5PvQbbbLZm?Y&8y2<+1jYvQPGdiSVXc>ECHgy+EyT)8PSF%Q? z(R4QaK6tE^zuh3GvX9TvS&bPL3gy2+?}S=IE7q{&XVh$Uh2LAKV|``#)$sH+-1<&0)GO{r7CXSujNQW~ydtj5U`*06QB2z` z4OPruWw6GP#j)EP6FtrS)QFz3R5+YU`0yESA~MN4{%tfWobN-4e2}ZU0}?gl*<*N_ z=2oiVgg9zZ@Ep$t*KAf7b(Cay(R<>IZ|U+Rsip&w&7YoY3kl(K;5* z!Dr0PpT3!T6S+qQGh4U!O*AT7mfkju2G--esn36oXKVB6+h zTB>2!uzH!2eu-GH_E!Js>TL_nyv;8heB;2Lrkm^m^r(YDd1PAmDAq~i`f4y9+E+wI z+aNn@;*PE$RkLV4Q}PS@AmEZ7g85wOo~q5mUe>W(MhC@UItz)Fp)1}dVyc09Od;dL zCJg$t$>U`!si6#_A7>b^cFKHVp$qbBEqeZLpj$ma4|Y$SNk24JJcRrf1vCEB7oqfW z)tSgyt>^Z)=^mbPSXtJrnWIoSY>4+R?Q*!+`z?jruSlW50GvhwCx3hF_!C1SemK96L(v z<{WG_)68S2M7%i;Hf(>I;@!2!8C2o`5*s3GtruTv=4N&6F#T$ENUId5ND$uvF||ll zDd=Fp>3?}kasBi;`^cS~Ja@DArzEW3!d}=GE>pLMu#Pp~pSe=LYSy9tt|;3i_o>KB zB=_Pd)qcW`4`$$b;rV5TiXgB4%}t}!)ko6n8a}R!J5q}4gaoVML1aQVnoZ|V8|Y&2 zUX@$t-~C+eaM+d=OIOaLj%*-gt}pi;MQ={G_%Se+-1Qu{1mqH)nM5s)@{xQ%|L~E9 z+v;qd)I+>jOh>CzST6cG4_Y&ECYB|eOr{VkRI<*;-Ot1X67D=p1YjwO5n zSy!{B27VuJwbdSql1ne<{N%=$?S88}T>^UbskW_ZiL{g0(;X!HcV!Hn^K$W@$dxjR zbF5|_$E6~7nchEWdmpwt!aIJaLhSaV{;G$3ulB`=_!c-Z-aULO&9lmedBgoa%EnS- zzY8BbbIEZAUC~#ysypqTR3&&svnhB>?K*z&=cG?qotq#Ts(7u<5ntj{-*4^k)ZCK0 znd}&q2f0EC-qkmGcSV>Lr?b7X7?t?$QC1d4DSxi!mD!X0bfHOCE>Y&Ft5FXtL(%2G zTM16R?`##lEl7MqQuw5I=l-q(s?7b!(3^>7(0A&u`-2?UJ@^7To(4Yy-#%>Evh7sQ z3i`O{tblZ5K!a;GFbL`Ca}nHy!ORcAgOoPtJ_8?lg(Ml#WT{Fq-B3*lXGt?{JE+A) znRM+A+wu2Tw;7kQgHq}p&U;kT7_C2HVqS;cNQEA&wynAf`wL3PeHH7Bc|48rFj428 zwxJNJWMJu`hp2^)(3KuLc^MDXGIGJ;6&OyoosKUpC?cFOt~kfkHjL2iiw|_TS@M=! zM&P15xM|qV&&0Pbp0EoQ^@aXY>2*Cmd-qCl5NneDU7i#SjE)GayVl=kQRKxyScPfU zmfY4hcha@5zHE9%r(PgU85KSy;G~qwnVNm#%E5OV?cN8W(Y_%#dpvdl3URIV5o!>B z<{(E>0IO}oAaxMric?940WQV!7bwA**~e{McnSCNP-9jglS+PXo-}H%1$woK(Vz;Q zhDbVM6IKZ^^%F%K6BbvyaEWPj1_$acu#I=UzDgb=b4Sef7W=-0`STztdtD!|_@t@_ z(N%_5h%vM7ITeOQZq*r>8g(tpEA^>Rx8As#$JE%YOR3=SY6Vp+^o$(EwZ}7GPDF^o zgBttAdg5zvwA5BkM3g$4YG7~5t%e+vhpp$E%J)Lv3QOQE4(< zM2_q!rp`+T0L@d!>sKA{NDca>H9P}!(Q`X0S>C@1o}-U)T@G)2RF9nEex;VxW;N-y@Vgq2F($i%ZGe zLGSw4kYpc9(%00?__l{jf=eIi+skzx(e&#~)u04x7Pe)*wWj8OQ{cS#fKftba|ZHI zysW#xT<+^N2=}dSMSnGX7X z{aI^~#X^k*s66^BeQ{=ct+Y7_ZqrpBn(EL%u;z8l2ieVN*w;^Lo;QcyuIp{;lF%9m zY;}`{k`xeIsCS^SD|vf^#E2R^T7!b^gvPVmEuFiGj$ z0MjW4PNi!N`Meswg8HMj$ziRswNWREVMJ*O^!K!uY_SR@Xae;<^(n!Shcb|y!OsW6 zt05>?<4>?KO{wQp)eC{|Nu$JNYVOdPeU;3`Sr1vz$*)-L?azJBBcxE#cm8$V$+6{? zN1ofmeA3(eqCz+ysk(Lkd1mJ;-|xo;^{2_rCpxaTR`)Q;vOT2~!k(X{I+>hwj-Ag< zR5YJW@c8>3uk1VJQ%66)7CUHiyqZc#9wHL%cf8(MSJzVAcksm~dRn2*$=H8xk2(9~ zo>Iga^yQSf=;kPOPtz_|XILgCxsiJZC3Z$SRA-lf8=VB(gIxQ1$ZZZlST* zsYUj?Sif09u_sSO<;Sj!+8kpO+MPTTVR~Zu(qaeCcw>2CM0}6obcyn8Mmj;u`$_@M zkj){sn%(I$ktcqwL|qNv{Iqzo&#Kgx%LzXUO#_3pm&-NjfjpIM~Q1GDAD zL9tSg=G|leljqx~xL=%9C%Z?WmW`QCfHLsm+Ki?le@jq35Bz*(> zNcX)&nE!TvUCZIZ-MZD1SLHKn?X!E!6?KF12NM0y6!)DR z&vyrt<4@fq8y_0DOKnE%fE<7EG{ zHDTY-jSLruSN@YpFn|4nljQx_uD*Ry@nfIy=~q5`hkM5UU%tB?oN1k`thbmlwH)X{ z7R~4BEg%nVUVj{HDKsk|z~(&OYgFT2^4&FP(qG(V@vks;pI#&x^8x2KRj#3aMspW8 z7W*olX(MoFzxI3_V`RXmhOP_aq>IJjDnx{Sz~MYm4Zq|W;eLLG`9T=r*rNE^iGBVP z8=Y@Qk6N|GPxc0~{f|$t-t{|Q-=7#)_oO;I9`W@*o9DPs@piurdivb?)%o7>-S!IK z({-_FL{{Mi=Eu{k>rMU#e$%OgVxl>zn?Cnavi*)v(_0MZ`Wo-mKiVHVc(s3ePuVZ# zE^zomYOAB)nT2-9iNEI@J?(bb%UV$$w`<21@%Tbn@^`r&J$IZ{ZGH}*NNL~Y1g)KF0W4|&~>i}zc-kf6Q5oX$Ztu_ zJUJQf6K(b>G-wJb_V=ovIM}t5b~p|?2t7ZVm_A9zr;I)6o|=sf!OwQ89kxGD(%sqg zJ#cNTpg)_Q;yFoNr1qN$^*HVZ(2bDIChdCjW*jrgq~JhY+5Y$7N#=VBU3 z*T!3-@riS*vfWu8nTFvwy6EE72-B1ujB2I%F42)BsJl(c=)#ykg%p?^^?^ImbbD$m z_;t!%n<*w~A}{H}C5K7F|9~(UAD<3by|M<0q_9v9=c`-_4qK?9gjjfTpOvj0f%QUw z9z~Y~h5}WV8B#%8hMTFk>r3d%cSwEYL~nm50*W}{`j-=}1;PKDJEAdt_^n_p=}nb3 ztXVCTJ(JP;$doJj7ly=L1j5VAxs&mWKGeM zalM4}9h1vvnOtF2rALG{BNHtSNC&xrN7{Fq*jDF&j$%N=eD}jH=8^Wq8f(sXZ6oO0?wNH@lGDtIMx53DSC?VDf1nzCn!@igC9zgpw3bjF>wWIKm(^Qv;C2p25hmxXltOuj zf|>3W!SGb0{D4^vY|OHC_mxO=c6Vc({?^R``j>Erxl`Th>o_TLoDasko9UV7y>6dS zqli-4wMVJQ5ip^;w0u|`^iE_8n-1W0xIxdXNSKG-XwVVEg9}yg*6c`lbW&N4u1%SX z?Pj6M#8{+V`1SB3%&?2pf~hzIooShLza)s*uFGPvl$j7pqqh5{v!d>s{#VHE^QRJB z@HQHI7T7GAEZP^!Bqn|rN>$p<;kUYOjYI01^UaN6iHXMhsfop{91K~V`t1;DHNCHv zMUvPz%d(JQrDM(D8zXYqcW*bm;@I2${+20fg2q!@sjy9uV8bd=Fvzk(__}@j9MSjB zqCxR#{4J4OK^&=qW0pr%^s>~;!C#KDp0QPvGoMH&i&TBNdY(OfeE9w33^R(=kMR=e zAkM-6Wyrw5{C_ZJ5|kZLA%H@PJ>Vv4JsN-fy@0hL?X&*`3-j_rnsEx*+bInydRlaG zyon+H%UjM>kG6&F2d%AV^z%|97Lpc+$jvwMe8%6(8R@Zfu7p?%)h%fVQewhYpJqSf zd@{xDMKE^#vB4g(c+*f0-6$V7Hov=0z7-p<#nWT|clS@Cpfk5sjO;a^-KQ5ZmbtaG zgef0!%_5a@aZVC@S*hJtKRo%aEUM(jcFmjn*c%g{Z}&XMEJ+&>X0xtbb=Tx}qzq%A z>)9^eTFC!+M{4k(l7o3i;bSPwm1V=W>+5LGm$&aGN%13CWsU z-f%6oi@bGXU!82EzT*`h(XPF$G@D!e`xTRp0o*$&_7q-9N9c6yNq9U&4FDZw&mxdnhmO|K`1Hys9fY16IiX<2;_Nr}H|Zg9Vy;DjO~M{Y^|9 z`VYRMj-ZgrZ;?o^MI@8bLw$Ibpgwg_;Vq%@bZ51fj)omp9WXP48DaK5VPomLwVL9O za-VKq75vPt$>(lmneCmMf%`Ttyqs9iog>vlVkBupbM~seK@i!lQYp)$9Ns13H(Pg; z?)Z3+s#R3Y-|T15?=K_0`moa8OLDZAiTDM6Buk1LgT~=g>SsdSrisj7b^<1IlC6Eb zH5G}$!9JzEwr^43Ly3(xlwBFhbh0rzylIufihNrE`(J9C_u($F-zg6s+?_t>An+YS z5%WhU zdSv0EBE-;-%tS)qMhuUS)}JDp8D^^)HT&gKM#e~X-m#;-(W;$Z4kB@QcI0+kitBLv ziS&l&{T0Op((uKU(pIZE^2APE3{j1GWr%OwihV9iZWC*Eu)WyraG6wAo0*8VHr646 zwfc3NFSi(Mji!nmEQ;lnx<^{@RaFzP4rQ=yYI!#(&in^oHIE$=-e9;%dz=Eg4s!3F+=pc-|6Aby=8K3?8iU?>zI z4yvI5nFq!WCd{qgG)T}iB?39l3=8;3$MvPDzFVsh=(&2Cpdpo^9g<5PiOnUnT!O2h zus?B9YJv%940wI_9*BW&U73(*DJrUke?0VWNJzgCx5dyqLQ0$NmCndmY27`cZwa!V z(NgJb*{dNWeY!tGcP%Lr*DN&N`9tZP}Y!+!N>) zc3H0ldmW>UNUNef%8Imf<-xglA?~#uNg%hI1v|Yl(637Z@hp!cH}@tuJq7&;R&n9c z4CIRk?A}4};JmeanMDM5N#Gr+>GBMrd1Y#k+tkzD2jPq@VFf(#Esf&OSruTRR#JQ6*Z%%lAtSeg`U@ZOGSF} zg_bN;2TB>!$4AsA(hq2y>r6gb-bn6tr+a&lRWlGQNZznN7XimqjsCQ>wohuu&ZF5# ztN+N)gCu$C%HkE0R2zS)Y0L9N%(F24Yrk_4b_rI63G)2`2Q`%)J0-9t7Rd~A-fd+L z-qXA_)pYE1M_RLly(t3H1@?UtSHGZaMB?(Gx6ZiSK%~%z+dx14UU*|n{$HZJrLyG(Ka)8X6$sgwtpaW}Sdaoc?b&*pdn+QusJU(Sj`B`cY}vPpG0 zYEWg6ND`YZ#*-HpsAuc7Pit(f?nBXUJ4HRgwV>I-u0fY_Pl>2Ml;UxpbQB5+h#ne=fSW{HEw7VR&s7y|hqDFn((b-FQ zx0AMqZdDKA>PFln*%LfbI9U8(&k!o)_+dL@fgaWgB^fYP7BtOGos*xL(DNWJr>%cE z;2uby{J~kJX4S7}SMAJwLciam$rUS4VnJ1MtC65wuIlw`sMw=uHt#&n%BvYoa=qt1jDw~~xA38|e?KdBN$bAlju`>jH zY6!i(D&bk({N1@WTO78O=WvNkFG2VGf8fun!L%Z=21`DBN9 zeAz)RT&XepghOZ`rq7i$d5<}}C*4zf_7!_4NCWn~ikEcg61*z-2P~G2ZX<35VWFMh zl*Got`J{F|5pep^L{J~oMS_hGbllo|a`foQ*Nw?_F_kPP)ng@I&AM{+PNEY!v-p1+ zx|x+Sa@Q_i@fO=;0m;2{0-Nh`#rK9S9_`~=95aesLggjM?vLg?vcTY$fnZIvXgLQe z%%;YDuX&09)dra=@B3%4P3^td;Ye!qi3Fd{;lOh8YCD2*tIpB*?_ir_%l+9uhEihH1MKS zLPJmDr02b7+h<6xf-}WejF?x>BAw~Z$^u;|5)WT}x#GRW{P8P~PI-977Pd!&B~MfP z<4<1~r4q63mLU-sTXvK^%P&q@Uls}4vypX`iyDm-5MoHsecv=O69eX*Oy4=66nbP2gnfytpkQ>B`-2bD8+T0uhtkBYZc8U5%p5s zlnQvS-ym1x1=KKS1%=xLq=6iWIzbdNRp}8_GPGI9@qIH&s7TBc6upGmWcXmET8*_YYD&B5b6kn!1FvOE4PmAIuIfWmO81HaDdTC+1s} z;S^__YmeKKF_7wQ@L4Wt_A)re_XpTeHeT5k4TL4L|2T*-;ACDTQOOEYe-6hVu!MW2 z7N3SFF&2_=E|T2?`}S`4Th=b$71S1{MkkYHP3yh2RE6?PSnGBAorv437FP2wM}#pt zsRxV%ODtXTqdabUusf>^ChuxY9eU~|`N|Vm;e0H0qp_cPjfOaZHnl>LPv5FVi)vvWPPqH>Q=eZbY>U57-VVtNZ zp5Y2@kmn=A^DiF`U!A*z3J<=ZDyP}}ayRXRXSa!Y0=XM<7H@)D#5f&YPa$;Jtjpfcb|-@lQ|c)*A%HfRlH zgYP?8e6-hW9nl2^HES9UH7QbS^|4gOYqY(QE_@4p}pVx%6kc>fb82Xm)KB5Di@ybC>-ux8>W!eWampyVB_&fV%F38Eaa2&n&4jOf}(!8m*f_}I#_NUaeBg*g$XLiztoXHN@5x#>vs&*J% zl>um|z;vt@E?0GMX5TbgQ^Yi0$|!RR>%3xdIyNjL)cR7#r{;_99-pWN60Zw%l=k zGRU*QqV@3-7F_~S{>0YEMga@RU@V#Yg{EkJ*Y|)0BCC6a|S&)cMiWY*p}JSp~{7$;-`FB6VJZWuxei;(Dg^9C&4@`b9$Xi`K|* z9Pzj%7AI3#5%iR5P2K2u6GWayIqEt+__EA3mz1uSqa#naGu$a_l7Z-v#@R4R~@Y617`Vxl@ z1^E5HUyEoomH~^HtD+Am1b{Wn=6l>;UY%Cw@PaI~5*cYS#oi9B2Fo{QzDjZmG2Da> zSC#D>MEIq&!X5^d+=*C8dX`@JR*O@0==*2A@~>BB?%(1co!$vIMv2!5x@8tf5GVZ7 z^kDqXb$bomYmiWNSb$wSj@wHT2pUz?3!5uR?07fc53=hV^AuTqWUurBKe{IVJj1rg zq=Ci@*fb>daWc8ln?V(8BO1fvEvMn#lIM?*XjI(OBcgPh@ViGC@bHZ@FiM?2M5lqO zs2QIS8e>w8RNj@P^QgP>t=XU}sqy`wuZHheqjh?$HfolG`^ayojf^kh-6h}y0{??i zCmSb70?ZK5++GE-C>RFqm-fAAJ324H{2-sVjKWLPK%x4^o6J^hxk@g#Vd%~hr&5AU z8$|-GdXYuAFsp)!f<=NtEnP5UhHkY{TKY})YR%aiyUzBOasBr_GSSY_ClSjSH&4`7 z43oB0gDVmsN5B`n!og9B2gb6K`jGWA8E_TV9=4P0b@Y&ks*>^MLDSH;`#56)-QV$5 zg7{Z2>F6a;1_t|ssRNg$vTG*@Z^8}7C}cIx-R;y#Dm$2XKx=2DopeUh4;lpzD(eHa9@UV9`ah`bT;#b}ryCan8d(^+I%mn7w=DABi;^QOF+D!v@@DSx2bHpW$g)zqrRT=$bcLg0)#?qdYVuSpN7`ik_`|CTUYa$Y9qEwHwVk zuvG>@a_kaweVsNtvUGFo2@(@64Xq4s;MMQ;#y4TVGh|LTAGb@yatYM?lNVwz7+ebg z<6sov#bM~bUL2y6_90il<$N2+$%q@5Lc=YMYXNpJw=XcYHpj5k13^Fyh&PC78KkI} zP;&`df&s4;|EDG4114_hz{G8k?WM$%EwzE~PoG9M1FJHmugF!vrgLmOdhe8*N0}#{ zbj%yj?=x3PeU4qvc&pjDPf<-a$WZ zx+z!@DfspQCfySlNPSVFtI-^)i5kT8(ap{UE5h}9gnvQP<7t+%g4`WZz4X56MCB(# zw5m6jkUQUWD&`m|9_FN{DajLYYjRiW5zmp0Qw1!(e%HYJau)SX7 z9TU+y^tk$oSywlm{_c0eo6A07zY7BJ`@by?^ZmhKy;!YAk}^_eM7IDWSU-*3>5Koe3-L>SLJ9;oJ&RWDjLfk4(E7&t>IU7cc?GfBrmK>q$&y zWIi}{ug0F7&Lq#^nn!|o@N0IS$+-sFo+lp~Kj%d@`xhWXurKM&C71>VpzuE!E`p8# z;9fn&@dX7qPHc{Y6yi@>@Jq(YP+iIFU77k)+;27;{#W>F1jZ5l)(DK{oH$8&oRPZrKdc^q1BF{@RR(OBOmR% zcjcxoXn37Q(K#wYJwCuQ@RI_b;C)U{4f-7zD};3E6SeH5WFT`wnqX&x`_A;iP*u*^ z$0+_5A@*eheLahJV$|zNU7XNmZ`~dD?;pB@ad|iMP8Q$ek?>kYCv;(!eG&}KOV`l$ zs*_H{?}v0Vd?$F^L4~b!Btq5_%S=j})H6+hpX+Zzmmie{|BUiAPuzhVli*8rhjD!% zxf11_gni@^m8HPt^NMdTn#acIY1ad~E)mWpfDHO4h>S}^(G^t^t8tKHC&apF+9H)B zZ2o)<>kTFbekAZa2w!rdvO%0WO$+hQY&)i5kf3SwFh7nl2i4B_sq)^8(S|@g9m-wf zCdv<5l+#T$i*wYz3L}Z(eW1Yb;Av*T0GcF8RJmL<+?DM32~Vth=&k5t$eJ6MgUO(( zz@4Er-O}y(BswD(ZWcUS+wC_^ zJ2S@|CE)MBOgDwk9Y<9xE*_qsnwT1dU&6ynkR{|#mULBt-SyJI?0G+{ow7It*19Ot z^0dKrv!zTL(@`GkDvwJy)w6s2 zEs0w@6oOJK*E=qf*GkxJh(%(VO5fZHx@Hu6sQq0pac8hjcq+D`#@Yn9mx$>B`Cb22Qwv5x+NQDpj( zZeIcl|H%6C0n6CxLoaR)@j>CYzi^Uefv{NKN1m-YxN|))NPS5olh%HCCaCuv;VK)L<=yH|0h)@f*xL^TBuNq#BrQ1MZh<*-aL@ zkiMqnLtkNA$;xjcs?RR5`N*3v#ao@=_7cD4{-D;_dC{|`zC-j^ce&3l>D49JhyM>= z6^KR&698;SsJgtZ9}}ijc;ILMCu}v5r#nm&ox`U3@hEkzz4KxvBd&u?!ahWq4D&kF z`Y+1gyvx9g%!q2ipPo`3!z~n`ykut_FcJC95knz=@Rp2ELy!Piy;|ERKa*gJ;_*_; zW9>g{>5m%LrfV|{tw~}?-_}elW?0Y?^G@f^YM2~*dz*JLrH0Eb;AL$SPgz21gmrju zn_b#V9OHr4QWccq0uK=td)~!Vh`&a38ySx1g6$9EJlFLS*bWB`%7d=DI5k_ z3iw)pb)3(s>r*{iPr6dB1Y#g-3($-$jcp~1s$Cy;Cx?zci{TXcuwmZPX^=Q}7KWBvO|?6w7R-r0V(CPE|YA7&=F1$l3BYQ+}5jL94XiDnWs(+N1l5+9_8GyJTU|_o& z7ca2b2;m3b(LkU8&N2EypsJa>s~sW}m>Ud6fGyzP-!UT`9t3gfu1*$iKv4t*I6ShB z7EYFSjy52~HnJZ*{h5TrBV%W6ZRKj^Xkq0B)QQ*y;^u5-VFf}!5oa^tfxVSAz?6%q ztDTLlI|#AO#=^>h{x)pI>g$vg_FILBM7nS(aOcc%pL>>OqMK&ACAML0^(5x z2>`UHgAi+NS|G%GUJpl0D_1uQCs!+wAh4ot2I4UX>HzToWFXoC@z{VMU?86zZ~>m* z0rA-5@Hl{YfHs{#5C~9@3y8-RpaY-}s10B*&;UT!J%D~-!>E`T4v&tVrMoS_LIe=d z1OX{O_nZFHtNll7{vr1ZJ@_v|1|(>2X5$91`H#dc>1JUCXrv$%U!L>>2n_fN_#fT*33LoV za4`7)kAU<4wH_ESOW*|*ninu!0x$s(90tb`1m-gQfb^mKd>|OG9TgGy06G!R5&2<& z%3NGw{D}PjSHSuC90ue?1VkM_10S$o1p-DeK!6trKt+NAAbwyj4F?7UAJ8~IBJe>$ zP(eN%Ab{X}a1b9KPzf+X`FVl70ssMEAda{qp2L6ufq?jd@`ySG0VfO>_>qRVUUUY= zkAt`(h(L4-5#WGwBFgds9sf}WkQZq8A`Yy$K=}nO0w2HN#Z?dhRB#|24uJw4{LwBi z&@LZfK45}~a=aiu;3@zx3vs;&P((iX1rrcgMEe&)Log7a0bnQ&qD?-a%mq_^5{I~6 z=ntYLUc?O;!6F3Nh^`^7K)8?xqPZUdhwm?${aw|M2teSEDxg4*fnt0JssPFn#Sv=o zQ(QlqiXX zX8>P$abp7ldUz4D0|wEL*d0Ju5HTDcdk|m*E}lAo0I6KW9DzaeBjyCC&tEaXmjM9v z;`x6ZB49-?YWb&M{Hxi&JR;yGFUq)pc>ikK72z;1(%eA6o8XJ`KU^W8EkDu$?+9!K z`w_ozs27Vl0E2!w<%`7?;LaTwEh5a57;G#tS7@%_~_!ifXwc#+{hz5>z4&pn57zD}wibapONMLa^+&M1&xI@=XfW=oBFQDyQ|GN9{Hvqo;`B?=F2H=zb@m-8V z90V~Tlf8SR4ZRutrsJa*78UW#Pm_)-T%HAR&33lad$sDEC2WNqp>o_l0*fgs`RY3e z@p7>j8UWv?OgvLue|N5PuV{$K<8!!a{k2yN=4v-askaz~MGG2j+jeY=M2OP^0_7h^ z*U*1mkB{wWM^Pa#GXrV6M@bXp@f$shQ@!L(1!7eEL$os**sVqCDw2tmP606;(l zAbemxUI-VMmjw)FVf%+ddAV9y;{etliURc+f^T7e>{woazhyLC-UKqmu|2rQq1mS4@od)AWSh#_wO`d<^}$h zJ}8X$xBfz52%z(8Sp;VMHBAt}3csZB@`8Vr5tNtzx4r>Okl^3w1GMEgp7QY_3_arS zZ#{tW!3BRy1Ay0WGy~A)H=5!6zxNF;_?sR;`60jQ5R@PGJAM582pI5h`uGI^9QI3^ z02l!r{+&+%g4ho5OByeLgMUrqM|=hN?|gz6Gv;5?AiwiU5C-7uU&|u!{jX{Kzv~qk z*i!#XJ{TA=f%$h?7#I%VxnI)+e(Ns`u;jmvEf|2pf29Qmu;uqOI0DZg{{Gejz$^Tw zlQ1X@Sgras9~=SSekse3KqtSX3I4`s7(n!|@_@mhz(&De%EICQtk2!m%ntY|WE{jY zg_fPS6=3T?Jep2U2s3?wS`;0v0bKIq25_ow?q;s;7iJL#7XS<3Ff+@l$>IDzYI=!h literal 0 HcmV?d00001 diff --git a/analysis/sf_e_110_stats.txt b/analysis/sf_e_110_stats.txt new file mode 100644 index 0000000..05a20c7 --- /dev/null +++ b/analysis/sf_e_110_stats.txt @@ -0,0 +1,6 @@ +Validation: PASS +Sum of Probs: 110.0000 (Target k=110) +Average Time: 169.60s +Gini Fairness: 51.6610% +Min Prob: 0.026144 +Max Prob: 1.000000 diff --git a/analyze.py b/analyze.py new file mode 100644 index 0000000..1b35083 --- /dev/null +++ b/analyze.py @@ -0,0 +1,168 @@ +# coding: utf-8 +import csv +import argparse +import logging # <-- NEW: Imported the logging module! +from argparse import ArgumentParser +from dataclasses import dataclass +from math import ceil, isclose +from pathlib import Path +import time +from typing import Dict, List, Any, NewType, Union, Tuple +import sys +import os + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from scipy.stats.mstats import gmean + +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), 'src'))) + +try: + from sortition_algorithms.committee_generation.leximin import find_distribution_leximin +except ImportError as e: + print(f"CRITICAL ERROR: Could not locate leximin.py. Details: {e}") + exit(1) + +AgentId = NewType("AgentId", Any) +FeatureCategory = NewType("FeatureCategory", str) +Feature = NewType("Feature", str) + +@dataclass +class FeatureInfo: + min: int + max: int + selected: int = 0 + remaining: int = 0 + +ProbAllocation = NewType("ProbAllocation", Dict[AgentId, float]) + +@dataclass +class Instance: + k: int + categories: Dict[FeatureCategory, Dict[Feature, FeatureInfo]] + agents: Dict[AgentId, Dict[FeatureCategory, Feature]] + +def read_instance(feature_file: Union[str, Path], pool_file: Union[str, Path], k: int) -> Instance: + feature_info = {} + with open(feature_file, "r", encoding="utf-8") as file: + reader = csv.DictReader(file) + for line in reader: + category = FeatureCategory(line["category"]) + feature = Feature(line["name"]) + if category not in feature_info: + feature_info[category] = {} + feature_info[category][feature] = FeatureInfo(min=int(line["min"]), max=int(line["max"])) + + categories = list(feature_info) + agents = {} + with open(pool_file, "r", encoding="utf-8") as file: + reader = csv.DictReader(file) + for i, line in enumerate(reader): + agent_id = AgentId(str(i)) + agents[agent_id] = {category: Feature(line[category]) for category in categories} + for category in categories: + feature_info[category][Feature(line[category])].remaining += 1 + return Instance(k=k, categories=feature_info, agents=agents) + +def leximin_probabilities(instance: Instance, backend: str = "mip") -> ProbAllocation: + portfolio, output_probs, _ = find_distribution_leximin( + instance.categories, + instance.agents, + instance.k, + [], + solver_backend=backend + ) + selection_probs = {agent_id: 0. for agent_id in instance.agents} + for panel, probability in zip(portfolio, output_probs): + for agent_id in panel: + selection_probs[agent_id] += probability + return ProbAllocation(selection_probs) + +@dataclass +class ProbAllocationStats: + gini: float + geometric_mean: float + min: float + max: float + +def compute_prob_allocation_stats(alloc: ProbAllocation) -> ProbAllocationStats: + n = len(alloc) + k = round(sum(alloc.values())) + sorted_probs = sorted(alloc.values()) + gini = sum((2 * i - n + 1) * prob for i, prob in enumerate(sorted_probs)) / (n * k) + geometric_mean = gmean([max(p, 1e-10) for p in sorted_probs]) + mini = min(sorted_probs) + maxi = max(sorted_probs) + return ProbAllocationStats(gini=gini, geometric_mean=geometric_mean, min=mini, max=maxi) + +def analyze_instance(instance_name: str, instance: Instance, backend: str, num_runs: int = 3): + Path("analysis").mkdir(exist_ok=True) + print(f"Running analysis with {backend} backend ({num_runs} runs)...") + + timings = [] + alloc = None + for i in range(num_runs): + start_time = time.time() + alloc = leximin_probabilities(instance, backend) + timings.append(time.time() - start_time) + print(f" Run {i + 1}/{num_runs} complete.") + + total_prob = sum(alloc.values()) + is_valid = isclose(total_prob, instance.k, rel_tol=1e-5) + + stats = compute_prob_allocation_stats(alloc) + + log_path = Path("analysis") / f"{instance_name}_{instance.k}_stats.txt" + with open(log_path, "w") as f: + f.write(f"Validation: {'PASS' if is_valid else 'FAIL'}\n") + f.write(f"Sum of Probs: {total_prob:.4f} (Target k={instance.k})\n") + f.write(f"Average Time: {sum(timings)/num_runs:.2f}s\n") + f.write(f"Gini Fairness: {stats.gini:.4%}\n") + f.write(f"Min Prob: {stats.min:.6f}\n") + f.write(f"Max Prob: {stats.max:.6f}\n") + + n = len(alloc) + probs = list(alloc.values()) + + plt.figure(figsize=(10, 3)) + + rng = np.random.default_rng(42) + jitter = 0.22 + y_vals = rng.uniform(-jitter, jitter, size=n) + + plt.scatter(probs, y_vals, s=18, alpha=0.55, color="tab:blue", label="Leximin") + + plt.axvline(x=instance.k/n, color='black', linestyle='--', linewidth=1, label=f'Equalized (k/n = {instance.k/n:.4f})') + plt.axvline(x=stats.min, color='red', linestyle=':', alpha=0.5, label=f'Min: {stats.min:.4f}') + plt.axvline(x=stats.max, color='orange', linestyle=':', alpha=0.5, label=f'Max: {stats.max:.4f}') + + plt.yticks([0], ["Leximin"]) + plt.grid(axis="x", alpha=0.25) + plt.ylim(-0.6, 0.6) + plt.title(f"Marginal Selection Probabilities: {instance_name} (n={n}, k={instance.k})") + plt.xlabel("Selection Probability") + plt.legend(loc="upper right") + + output_path = Path("analysis") / f"{instance_name}_{instance.k}_plot.pdf" + plt.savefig(output_path, bbox_inches="tight") + plt.close() + + print(f"VALIDATION: {'✅ PASS' if is_valid else '❌ FAIL'} (Sum={total_prob:.2f})") + print(f"Gini Score: {stats.gini:.4%}") + print(f"Min/Max Probs: {stats.min:.6f} / {stats.max:.6f}") + +if __name__ == '__main__': + # --- NEW: Turn on the "chatty" debug logs! --- + # To make it less spammy later, change logging.DEBUG to logging.INFO + logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s') + + parser = ArgumentParser() + parser.add_argument('instance_name') + parser.add_argument('panel_size', type=int) + parser.add_argument('--backend', default='mip') + args = parser.parse_args() + + data_path = Path("data") / f"{args.instance_name}_{args.panel_size}" + instance = read_instance(data_path / "categories.csv", data_path / "respondents.csv", args.panel_size) + analyze_instance(args.instance_name, instance, args.backend) \ No newline at end of file diff --git a/analyze_maximin.py b/analyze_maximin.py new file mode 100644 index 0000000..3025c5b --- /dev/null +++ b/analyze_maximin.py @@ -0,0 +1,157 @@ +# coding: utf-8 +import csv +import argparse +from argparse import ArgumentParser +from dataclasses import dataclass +from math import ceil, isclose +from pathlib import Path +import time +from typing import Dict, List, Any, NewType, Union, Tuple +import sys +import os + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from scipy.stats.mstats import gmean + +sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), 'src'))) + +# IMPORT SWAPPED TO MAXIMIN +try: + from sortition_algorithms.committee_generation.maximin import find_distribution_maximin +except ImportError as e: + print(f"CRITICAL ERROR: Could not locate maximin.py. Details: {e}") + exit(1) + +AgentId = NewType("AgentId", Any) +FeatureCategory = NewType("FeatureCategory", str) +Feature = NewType("Feature", str) + +@dataclass +class FeatureInfo: + min: int + max: int + selected: int = 0 + remaining: int = 0 + +ProbAllocation = NewType("ProbAllocation", Dict[AgentId, float]) + +@dataclass +class Instance: + k: int + categories: Dict[FeatureCategory, Dict[Feature, FeatureInfo]] + agents: Dict[AgentId, Dict[FeatureCategory, Feature]] + +def read_instance(feature_file: Union[str, Path], pool_file: Union[str, Path], k: int) -> Instance: + feature_info = {} + with open(feature_file, "r", encoding="utf-8") as file: + reader = csv.DictReader(file) + for line in reader: + category = FeatureCategory(line["category"]) + feature = Feature(line["name"]) + if category not in feature_info: + feature_info[category] = {} + feature_info[category][feature] = FeatureInfo(min=int(line["min"]), max=int(line["max"])) + + categories = list(feature_info) + agents = {} + with open(pool_file, "r", encoding="utf-8") as file: + reader = csv.DictReader(file) + for i, line in enumerate(reader): + agent_id = AgentId(str(i)) + agents[agent_id] = {category: Feature(line[category]) for category in categories} + for category in categories: + feature_info[category][Feature(line[category])].remaining += 1 + return Instance(k=k, categories=feature_info, agents=agents) + +def maximin_probabilities(instance: Instance, backend: str = "mip") -> ProbAllocation: + portfolio, output_probs, _ = find_distribution_maximin( + instance.categories, + instance.agents, + instance.k, + [], + solver_backend=backend + ) + selection_probs = {agent_id: 0. for agent_id in instance.agents} + for panel, probability in zip(portfolio, output_probs): + for agent_id in panel: + selection_probs[agent_id] += probability + return ProbAllocation(selection_probs) + +@dataclass +class ProbAllocationStats: + gini: float + geometric_mean: float + min: float + max: float + +def compute_prob_allocation_stats(alloc: ProbAllocation) -> ProbAllocationStats: + n = len(alloc) + k = round(sum(alloc.values())) + sorted_probs = sorted(alloc.values()) + gini = sum((2 * i - n + 1) * prob for i, prob in enumerate(sorted_probs)) / (n * k) + geometric_mean = gmean([max(p, 1e-10) for p in sorted_probs]) + mini = min(sorted_probs) + maxi = max(sorted_probs) + return ProbAllocationStats(gini=gini, geometric_mean=geometric_mean, min=mini, max=maxi) + +def analyze_instance(instance_name: str, instance: Instance, backend: str, num_runs: int = 1): + Path("analysis").mkdir(exist_ok=True) + print(f"Running MAXIMIN analysis with {backend} backend ({num_runs} runs)...") + + timings = [] + alloc = None + for i in range(num_runs): + start_time = time.time() + alloc = maximin_probabilities(instance, backend) + timings.append(time.time() - start_time) + print(f" Run {i + 1}/{num_runs} complete.") + + total_prob = sum(alloc.values()) + is_valid = isclose(total_prob, instance.k, rel_tol=1e-5) + + stats = compute_prob_allocation_stats(alloc) + + # Plotting (Scatter Strip Plot Aesthetic) + n = len(alloc) + probs = list(alloc.values()) + + plt.figure(figsize=(10, 3)) + + rng = np.random.default_rng(42) + jitter = 0.22 + y_vals = rng.uniform(-jitter, jitter, size=n) + + # Changed dot color to Red for Maximin + plt.scatter(probs, y_vals, s=18, alpha=0.55, color="tab:red", label="Maximin") + + plt.axvline(x=instance.k/n, color='black', linestyle='--', linewidth=1, label=f'Equalized (k/n = {instance.k/n:.4f})') + plt.axvline(x=stats.min, color='red', linestyle=':', alpha=0.5, label=f'Min: {stats.min:.4f}') + plt.axvline(x=stats.max, color='orange', linestyle=':', alpha=0.5, label=f'Max: {stats.max:.4f}') + + plt.yticks([0], ["Maximin"]) + plt.grid(axis="x", alpha=0.25) + plt.ylim(-0.6, 0.6) + plt.title(f"Maximin Selection Probabilities: {instance_name} (n={n}, k={instance.k})") + plt.xlabel("Selection Probability") + plt.legend(loc="upper right") + + output_path = Path("analysis") / f"{instance_name}_{instance.k}_maximin_plot.pdf" + plt.savefig(output_path, bbox_inches="tight") + plt.close() + + print(f"VALIDATION: {'✅ PASS' if is_valid else '❌ FAIL'} (Sum={total_prob:.2f})") + print(f"Gini Score: {stats.gini:.4%}") + print(f"Min/Max Probs: {stats.min:.6f} / {stats.max:.6f}") + +if __name__ == '__main__': + parser = ArgumentParser() + parser.add_argument('instance_name') + parser.add_argument('panel_size', type=int) + parser.add_argument('--backend', default='mip') + args = parser.parse_args() + + data_path = Path("data") / f"{args.instance_name}_{args.panel_size}" + instance = read_instance(data_path / "categories.csv", data_path / "respondents.csv", args.panel_size) + analyze_instance(args.instance_name, instance, args.backend) \ No newline at end of file diff --git a/src/sortition_algorithms/committee_generation/leximin.py b/src/sortition_algorithms/committee_generation/leximin.py index 1a17c89..628a3a8 100644 --- a/src/sortition_algorithms/committee_generation/leximin.py +++ b/src/sortition_algorithms/committee_generation/leximin.py @@ -249,8 +249,8 @@ def _run_leximin_main_loop( fixed_probabilities: dict[str, float] = {} reduction_counter = 0 - while len(fixed_probabilities) < people.count: - logger.debug(f"Fixed {len(fixed_probabilities)}/{people.count} probabilities.") + while len(fixed_probabilities) < len(people): + logger.debug(f"Fixed {len(fixed_probabilities)}/{len(people)} probabilities.") dual_model, dual_agent_vars, dual_cap_var = _dual_leximin_stage( people, @@ -321,7 +321,7 @@ def find_distribution_leximin( # Find initial committees that cover every possible agent committees, covered_agents, initial_report = generate_initial_committees( - new_committee_solver, agent_vars, 3 * people.count + new_committee_solver, agent_vars, 3 * len(people) ) report.add_report(initial_report) diff --git a/src/sortition_algorithms/committee_generation/maximin.py b/src/sortition_algorithms/committee_generation/maximin.py index d5265ce..09ca690 100644 --- a/src/sortition_algorithms/committee_generation/maximin.py +++ b/src/sortition_algorithms/committee_generation/maximin.py @@ -297,7 +297,7 @@ def find_distribution_maximin( # Find initial committees that cover every possible agent committees, covered_agents, init_report = generate_initial_committees( - new_committee_solver, agent_vars, people.count + new_committee_solver, agent_vars, len(people) ) report.add_report(init_report) diff --git a/src/sortition_algorithms/committee_generation/newleximin.py b/src/sortition_algorithms/committee_generation/newleximin.py new file mode 100644 index 0000000..55f8467 --- /dev/null +++ b/src/sortition_algorithms/committee_generation/newleximin.py @@ -0,0 +1,268 @@ +# ABOUTME: Leximin algorithm for committee generation. +# ABOUTME: Maximizes minimum probability, breaking ties by second-lowest, third-lowest, etc. + +import logging +from typing import Any + +import numpy + +try: + import gurobipy as grb + + GUROBI_AVAILABLE = True +except ImportError: + GUROBI_AVAILABLE = False + grb = None + +from sortition_algorithms.committee_generation.common import ( + EPS, + generate_initial_committees, + ilp_results_to_committee, + setup_committee_generation, +) +from sortition_algorithms.committee_generation.solver import Solver, SolverSense, solver_sum +from sortition_algorithms.features import FeatureCollection +from sortition_algorithms.people import People +from sortition_algorithms.utils import RunReport, logger + + +def _dual_leximin_stage( + people: People, + committees: set[frozenset[str]], + fixed_probabilities: dict[str, float], +) -> tuple: + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + assert len(committees) != 0 + + model = grb.Model() + agent_vars = {agent_id: model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for agent_id in people} + cap_var = model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) + model.addConstr( + grb.quicksum(agent_vars[agent_id] for agent_id in people if agent_id not in fixed_probabilities) == 1 + ) + for committee in committees: + model.addConstr(grb.quicksum(agent_vars[agent_id] for agent_id in committee) <= cap_var) + model.setObjective( + cap_var + - grb.quicksum(fixed_probabilities[agent_id] * agent_vars[agent_id] for agent_id in fixed_probabilities), + grb.GRB.MINIMIZE, + ) + + model.setParam("Method", 2) # optimize via barrier only + model.setParam("Crossover", 0) # deactivate cross-over + + return model, agent_vars, cap_var + + +def _add_report_update( + report: RunReport, dual_obj: float, value: float, upper: float, committees: set[frozenset[str]] +) -> None: + at_most = f"{dual_obj - upper + value:.2%}" + dual_obj_str = f"{dual_obj:.2%}" + gap_str = f"{value - upper:.2%}" + report.add_message_and_log( + "leximin_is_at_most", + log_level=logging.DEBUG, + at_most=at_most, + dual_obj_str=dual_obj_str, + num_committees=len(committees), + gap_str=gap_str, + ) + + +def _run_leximin_column_generation_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + dual_model: Any, + dual_agent_vars: dict[str, Any], + dual_cap_var: Any, + committees: set[frozenset[str]], + fixed_probabilities: dict[str, float], + people: People, + reduction_counter: int, + report: RunReport, +) -> tuple[bool, int]: + while True: + dual_model.optimize() + if dual_model.status != grb.GRB.OPTIMAL: + for agent_id in fixed_probabilities: + fixed_probabilities[agent_id] = max(0.0, fixed_probabilities[agent_id] - 0.0001) + dual_model, dual_agent_vars, dual_cap_var = _dual_leximin_stage( + people, + committees, + fixed_probabilities, + ) + logger.debug(f"Status: {dual_model.status} REDUCE PROBS for {reduction_counter}th time.") + reduction_counter += 1 + continue + + agent_weights = {agent_id: agent_var.x for agent_id, agent_var in dual_agent_vars.items()} + new_committee_solver.set_objective( + solver_sum(agent_weights[agent_id] * agent_vars[agent_id] for agent_id in people), SolverSense.MAXIMIZE + ) + new_committee_solver.optimize() + new_set = ilp_results_to_committee(new_committee_solver, agent_vars) + value = new_committee_solver.get_objective_value() + + upper = dual_cap_var.x + dual_obj = dual_model.objVal + + _add_report_update(report, dual_obj, value, upper, committees) + if value <= upper + EPS: + for agent_id, agent_weight in agent_weights.items(): + if agent_weight > EPS and agent_id not in fixed_probabilities: + fixed_probabilities[agent_id] = max(0, dual_obj) + return True, reduction_counter + + assert new_set not in committees + committees.add(new_set) + dual_model.addConstr(grb.quicksum(dual_agent_vars[agent_id] for agent_id in new_set) <= dual_cap_var) + + +def _solve_leximin_primal_for_final_probabilities( + committees: set[frozenset[str]], fixed_probabilities: dict[str, float] +) -> list[float]: + primal = grb.Model() + committee_vars = [primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for _ in committees] + eps = primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) + primal.addConstr(grb.quicksum(committee_vars) == 1) + for agent_id, prob in fixed_probabilities.items(): + agent_probability = grb.quicksum( + comm_var for committee, comm_var in zip(committees, committee_vars, strict=False) if agent_id in committee + ) + primal.addConstr(agent_probability >= prob - eps) + primal.setObjective(eps, grb.GRB.MINIMIZE) + primal.optimize() + + probabilities = numpy.array([comm_var.x for comm_var in committee_vars]).clip(0, 1) + return list(probabilities / sum(probabilities)) + + +def _run_leximin_main_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + committees: set[frozenset[str]], + people: People, + report: RunReport, +) -> dict[str, float]: + fixed_probabilities: dict[str, float] = {} + reduction_counter = 0 + + while len(fixed_probabilities) < people.count: + logger.debug(f"Fixed {len(fixed_probabilities)}/{people.count} probabilities.") + + dual_model, dual_agent_vars, dual_cap_var = _dual_leximin_stage( + people, + committees, + fixed_probabilities, + ) + + should_break, reduction_counter = _run_leximin_column_generation_loop( + new_committee_solver, + agent_vars, + dual_model, + dual_agent_vars, + dual_cap_var, + committees, + fixed_probabilities, + people, + reduction_counter, + report, + ) + if should_break: + break + + return fixed_probabilities + + +def find_distribution_leximin( + features: FeatureCollection, + people: People, + number_people_wanted: int, + check_same_address_columns: list[str], + solver_backend: str = "highspy", +) -> tuple[list[frozenset[str]], list[float], RunReport]: + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + report = RunReport() + report.add_message_and_log("using_leximin_algorithm", logging.INFO) + grb.setParam("OutputFlag", 0) + + new_committee_solver, agent_vars = setup_committee_generation( + features, people, number_people_wanted, check_same_address_columns, solver_backend + ) + + # cap initial committees at 100 to speedup + committees, covered_agents, initial_report = generate_initial_committees( + new_committee_solver, agent_vars, 100 + ) + report.add_report(initial_report) + + fixed_probabilities = _run_leximin_main_loop(new_committee_solver, agent_vars, committees, people, report) + + probabilities_normalised = _solve_leximin_primal_for_final_probabilities(committees, fixed_probabilities) + + # ========================================== + # WHITEBOARD VALIDATION CHECKS + # ========================================== + total_prob = sum(probabilities_normalised) + min_prob = min(fixed_probabilities.values()) if fixed_probabilities else 0.0 + + print("\n" + "="*50) + print(f"WHITEBOARD CHECKS (Dataset size: {people.count})") + print("="*50) + + # Check 1: Do probabilities add to 1? + print(f"1. Do probabilities add to 1? -> {total_prob:.6f} ", end="") + if abs(total_prob - 1.0) < 1e-5: + print("(✅ PASS)") + else: + print("(❌ FAIL)") + + # Check 2: The exact objective (lowest selection probability) + print(f"2. Objective (Lowest Selection Prob) -> {min_prob:.6f}") + print("="*50) + + # ========================================== + # STRIP PLOTS & CSV GENERATION + # ========================================== + try: + import csv + import matplotlib.pyplot as plt + + dataset_size = people.count + print(f"\nGenerating CSV and Strip Plot for dataset size {dataset_size}...") + + # save to CSV + csv_filename = f"individual_probs_{dataset_size}.csv" + with open(csv_filename, mode='w', newline='') as file: + writer = csv.writer(file) + writer.writerow(["Agent_ID", "Selection_Probability"]) + for agent_id, prob in fixed_probabilities.items(): + writer.writerow([agent_id, f"{prob:.6f}"]) + + # generate Strip Plot + prob_values = list(fixed_probabilities.values()) + plt.figure(figsize=(10, 3)) + plt.scatter(prob_values, [0] * len(prob_values), alpha=0.3, color='blue', edgecolors='none') + plt.title(f'Strip Plot of Selection Probabilities (N={dataset_size})', fontsize=14) + plt.xlabel('Probability of being selected', fontsize=12) + plt.yticks([]) + plt.gca().spines['left'].set_visible(False) + plt.gca().spines['right'].set_visible(False) + plt.gca().spines['top'].set_visible(False) + plt.tight_layout() + + img_filename = f"strip_plot_{dataset_size}.png" + plt.savefig(img_filename, dpi=300) + plt.close() # Close figure to save memory + print(f"✅ Saved {csv_filename} and {img_filename}!\n") + except Exception as e: + print(f"⚠️ Could not generate plot/csv: {e}\n") + + return list(committees), probabilities_normalised, report \ No newline at end of file diff --git a/src/sortition_algorithms/committee_generation/og_leximin.py b/src/sortition_algorithms/committee_generation/og_leximin.py new file mode 100644 index 0000000..628a3a8 --- /dev/null +++ b/src/sortition_algorithms/committee_generation/og_leximin.py @@ -0,0 +1,334 @@ +# ABOUTME: Leximin algorithm for committee generation. +# ABOUTME: Maximizes minimum probability, breaking ties by second-lowest, third-lowest, etc. + +import logging +from typing import Any + +import numpy + +try: + import gurobipy as grb + + GUROBI_AVAILABLE = True +except ImportError: + GUROBI_AVAILABLE = False + grb = None + +from sortition_algorithms.committee_generation.common import ( + EPS, + generate_initial_committees, + ilp_results_to_committee, + setup_committee_generation, +) +from sortition_algorithms.committee_generation.solver import Solver, SolverSense, solver_sum +from sortition_algorithms.features import FeatureCollection +from sortition_algorithms.people import People +from sortition_algorithms.utils import RunReport, logger + + +def _dual_leximin_stage( + people: People, + committees: set[frozenset[str]], + fixed_probabilities: dict[str, float], +) -> tuple: + """Implements the dual LP described in `find_distribution_leximin`, but where P only ranges over the panels + in `committees` rather than over all feasible panels: + minimize ŷ - Σ_{i in fixed_probabilities} fixed_probabilities[i] * yᵢ + s.t. Σ_{i ∈ P} yᵢ ≤ ŷ ∀ P + Σ_{i not in fixed_probabilities} yᵢ = 1 + ŷ, yᵢ ≥ 0 ∀ i + + Args: + people: People object with all agents + committees: set of feasible committees + fixed_probabilities: dict mapping agent_id to fixed probability + + Returns: + tuple of (grb.Model, dict[str, grb.Var], grb.Var) - (model, agent_vars, cap_var) + + Raises: + RuntimeError: If Gurobi is not available + """ + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + assert len(committees) != 0 + + model = grb.Model() + agent_vars = {agent_id: model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for agent_id in people} # yᵢ + cap_var = model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) # ŷ + model.addConstr( + grb.quicksum(agent_vars[agent_id] for agent_id in people if agent_id not in fixed_probabilities) == 1 + ) + for committee in committees: + model.addConstr(grb.quicksum(agent_vars[agent_id] for agent_id in committee) <= cap_var) + model.setObjective( + cap_var + - grb.quicksum(fixed_probabilities[agent_id] * agent_vars[agent_id] for agent_id in fixed_probabilities), + grb.GRB.MINIMIZE, + ) + + # Change Gurobi configuration to encourage strictly complementary ("inner") solutions. These solutions will + # typically allow to fix more probabilities per outer loop of the leximin algorithm. + model.setParam("Method", 2) # optimize via barrier only + model.setParam("Crossover", 0) # deactivate cross-over + + return model, agent_vars, cap_var + + +def _add_report_update( + report: RunReport, dual_obj: float, value: float, upper: float, committees: set[frozenset[str]] +) -> None: + at_most = f"{dual_obj - upper + value:.2%}" + dual_obj_str = f"{dual_obj:.2%}" + gap_str = f"{value - upper:.2%}" + report.add_message_and_log( + "leximin_is_at_most", + log_level=logging.DEBUG, + at_most=at_most, + dual_obj_str=dual_obj_str, + num_committees=len(committees), + gap_str=gap_str, + ) + + +def _run_leximin_column_generation_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + dual_model: Any, # grb.Model + dual_agent_vars: dict[str, Any], # dict[str, grb.Var] + dual_cap_var: Any, # grb.Var + committees: set[frozenset[str]], + fixed_probabilities: dict[str, float], + people: People, + reduction_counter: int, + report: RunReport, +) -> tuple[bool, int]: + """Run the column generation inner loop for leximin optimization. + + The primal LP being solved by column generation, with a variable x_P for each feasible panel P: + + maximize z + s.t. Σ_{P : i ∈ P} x_P ≥ z ∀ i not in fixed_probabilities + Σ_{P : i ∈ P} x_P ≥ fixed_probabilities[i] ∀ i in fixed_probabilities + Σ_P x_P ≤ 1 (This should be thought of as equality, and wlog. + optimal solutions have equality, but simplifies dual) + x_P ≥ 0 ∀ P + + We instead solve its dual linear program: + minimize ŷ - Σ_{i in fixed_probabilities} fixed_probabilities[i] * yᵢ + s.t. Σ_{i ∈ P} yᵢ ≤ ŷ ∀ P + Σ_{i not in fixed_probabilities} yᵢ = 1 + ŷ, yᵢ ≥ 0 ∀ i + + Args: + new_committee_solver: Solver for finding committees + agent_vars: agent variables in the committee solver + dual_model: Gurobi model for dual LP + dual_agent_vars: agent variables in dual model + dual_cap_var: capacity variable in dual model + committees: set of committees (modified in-place) + fixed_probabilities: probabilities that have been fixed (modified in-place) + people: People object with all agents + reduction_counter: counter for probability reductions + output_lines: list of output messages (modified in-place) + + Returns: + tuple of (should_break_outer_loop, updated_reduction_counter) + """ + while True: + dual_model.optimize() + if dual_model.status != grb.GRB.OPTIMAL: + # In theory, the LP is feasible in the first iterations, and we only add constraints (by fixing + # probabilities) that preserve feasibility. Due to floating-point issues, however, it may happen that + # Gurobi still cannot satisfy all the fixed probabilities in the primal (meaning that the dual will be + # unbounded). In this case, we slightly relax the LP by slightly reducing all fixed probabilities. + for agent_id in fixed_probabilities: + # Relax all fixed probabilities by a small constant + fixed_probabilities[agent_id] = max(0.0, fixed_probabilities[agent_id] - 0.0001) + dual_model, dual_agent_vars, dual_cap_var = _dual_leximin_stage( + people, + committees, + fixed_probabilities, + ) + logger.debug(f"Status: {dual_model.status} REDUCE PROBS for {reduction_counter}th time.") + reduction_counter += 1 + continue + + # Find the panel P for which Σ_{i ∈ P} yᵢ is largest, i.e., for which Σ_{i ∈ P} yᵢ ≤ ŷ is tightest + agent_weights = {agent_id: agent_var.x for agent_id, agent_var in dual_agent_vars.items()} + new_committee_solver.set_objective( + solver_sum(agent_weights[agent_id] * agent_vars[agent_id] for agent_id in people), SolverSense.MAXIMIZE + ) + new_committee_solver.optimize() + new_set = ilp_results_to_committee(new_committee_solver, agent_vars) # panel P + value = new_committee_solver.get_objective_value() # Σ_{i ∈ P} yᵢ + + upper = dual_cap_var.x # ŷ + dual_obj = dual_model.objVal # ŷ - Σ_{i in fixed_probabilities} fixed_probabilities[i] * yᵢ + + _add_report_update(report, dual_obj, value, upper, committees) + if value <= upper + EPS: + # Within numeric tolerance, the panels in `committees` are enough to constrain the dual, i.e., they are + # enough to support an optimal primal solution. + for agent_id, agent_weight in agent_weights.items(): + if agent_weight > EPS and agent_id not in fixed_probabilities: + # `agent_weight` is the dual variable yᵢ of the constraint "Σ_{P : i ∈ P} x_P ≥ z" for + # i = `agent_id` in the primal LP. If yᵢ is positive, this means that the constraint must be + # binding in all optimal solutions [1], and we can fix `agent_id`'s probability to the + # optimal value of the primal/dual LP. + # [1] Theorem 3.3 in: Renato Pelessoni. Some remarks on the use of the strict complementarity in + # checking coherence and extending coherent probabilities. 1998. + fixed_probabilities[agent_id] = max(0, dual_obj) + return True, reduction_counter # Break outer loop + + # Given that Σ_{i ∈ P} yᵢ > ŷ, the current solution to `dual_model` is not yet a solution to the dual. + # Thus, add the constraint for panel P and recurse. + assert new_set not in committees + committees.add(new_set) + dual_model.addConstr(grb.quicksum(dual_agent_vars[agent_id] for agent_id in new_set) <= dual_cap_var) + + +def _solve_leximin_primal_for_final_probabilities( + committees: set[frozenset[str]], fixed_probabilities: dict[str, float] +) -> list[float]: + """Solve the final primal problem to get committee probabilities from fixed agent probabilities. + + The previous algorithm computed the leximin selection probabilities of each agent and a set of panels such that + the selection probabilities can be obtained by randomizing over these panels. Here, such a randomization is found. + + Args: + committees: set of committees + fixed_probabilities: fixed probabilities for each agent + + Returns: + list of normalized probabilities for each committee + """ + primal = grb.Model() + # Variables for the output probabilities of the different panels + committee_vars = [primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for _ in committees] + # To avoid numerical problems, we formally minimize the largest downward deviation from the fixed probabilities. + eps = primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) + primal.addConstr(grb.quicksum(committee_vars) == 1) # Probabilities add up to 1 + for agent_id, prob in fixed_probabilities.items(): + agent_probability = grb.quicksum( + comm_var for committee, comm_var in zip(committees, committee_vars, strict=False) if agent_id in committee + ) + primal.addConstr(agent_probability >= prob - eps) + primal.setObjective(eps, grb.GRB.MINIMIZE) + primal.optimize() + + # Bound variables between 0 and 1 and renormalize, because numpy.random.choice is sensitive to small deviations here + probabilities = numpy.array([comm_var.x for comm_var in committee_vars]).clip(0, 1) + return list(probabilities / sum(probabilities)) + + +def _run_leximin_main_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + committees: set[frozenset[str]], + people: People, + report: RunReport, +) -> dict[str, float]: + """Run the main leximin optimization loop that fixes probabilities iteratively. + + The outer loop maximizes the minimum of all unfixed probabilities while satisfying the fixed probabilities. + In each iteration, at least one more probability is fixed, but often more than one. + + Args: + new_committee_solver: Solver for finding committees + agent_vars: agent variables in the committee solver + committees: set of committees (modified in-place) + people: People object with all agents + output_lines: list of output messages (modified in-place) + + Returns: + dict mapping agent_id to fixed probability + """ + fixed_probabilities: dict[str, float] = {} + reduction_counter = 0 + + while len(fixed_probabilities) < len(people): + logger.debug(f"Fixed {len(fixed_probabilities)}/{len(people)} probabilities.") + + dual_model, dual_agent_vars, dual_cap_var = _dual_leximin_stage( + people, + committees, + fixed_probabilities, + ) + + # Run column generation inner loop + should_break, reduction_counter = _run_leximin_column_generation_loop( + new_committee_solver, + agent_vars, + dual_model, + dual_agent_vars, + dual_cap_var, + committees, + fixed_probabilities, + people, + reduction_counter, + report, + ) + if should_break: + break + + return fixed_probabilities + + +def find_distribution_leximin( + features: FeatureCollection, + people: People, + number_people_wanted: int, + check_same_address_columns: list[str], + solver_backend: str = "highspy", +) -> tuple[list[frozenset[str]], list[float], RunReport]: + """Find a distribution over feasible committees that maximizes the minimum probability of an agent being selected + (just like maximin), but breaks ties to maximize the second-lowest probability, breaks further ties to maximize the + third-lowest probability and so forth. + + Args: + features: FeatureCollection with min/max quotas + people: People object with pool members + number_people_wanted: desired size of the panel + check_same_address_columns: Address columns for household identification, or empty + if no address checking to be done. + solver_backend: solver backend to use ("highspy" or "mip") for committee discovery. + Note: The dual LP still uses Gurobi. + + Returns: + tuple of (committees, probabilities, output_lines) + - committees: list of feasible committees (frozenset of agent IDs) + - probabilities: list of probabilities for each committee + - output_lines: list of debug strings + + Raises: + RuntimeError: If Gurobi is not available + """ + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + report = RunReport() + report.add_message_and_log("using_leximin_algorithm", logging.INFO) + grb.setParam("OutputFlag", 0) + + # Set up an ILP that can be used for discovering new feasible committees + new_committee_solver, agent_vars = setup_committee_generation( + features, people, number_people_wanted, check_same_address_columns, solver_backend + ) + + # Find initial committees that cover every possible agent + committees, covered_agents, initial_report = generate_initial_committees( + new_committee_solver, agent_vars, 3 * len(people) + ) + report.add_report(initial_report) + + # Run the main leximin optimization loop to fix agent probabilities + fixed_probabilities = _run_leximin_main_loop(new_committee_solver, agent_vars, committees, people, report) + + # Convert fixed agent probabilities to committee probabilities + probabilities_normalised = _solve_leximin_primal_for_final_probabilities(committees, fixed_probabilities) + + return list(committees), probabilities_normalised, report diff --git a/src/sortition_algorithms/committee_generation/op1_leximin.py b/src/sortition_algorithms/committee_generation/op1_leximin.py new file mode 100644 index 0000000..3e1933a --- /dev/null +++ b/src/sortition_algorithms/committee_generation/op1_leximin.py @@ -0,0 +1,250 @@ +# ABOUTME: Leximin algorithm for committee generation. +# ABOUTME: Maximizes minimum probability, breaking ties by second-lowest, third-lowest, etc. +# Option 1: The "Rebuilding the World" Bottleneck + +import logging +from typing import Any + +import numpy + +try: + import gurobipy as grb + + GUROBI_AVAILABLE = True +except ImportError: + GUROBI_AVAILABLE = False + grb = None + +from sortition_algorithms.committee_generation.common import ( + EPS, + generate_initial_committees, + ilp_results_to_committee, + setup_committee_generation, +) +from sortition_algorithms.committee_generation.solver import Solver, SolverSense, solver_sum +from sortition_algorithms.features import FeatureCollection +from sortition_algorithms.people import People +from sortition_algorithms.utils import RunReport, logger + + +def _setup_dual_model( + people: People, + committees: set[frozenset[str]], +) -> tuple: + """Sets up the initial Gurobi model for the dual LP exactly once.""" + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + assert len(committees) != 0 + + model = grb.Model() + + # Variables: yᵢ for each agent, and ŷ (cap_var) + agent_vars = {agent_id: model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for agent_id in people} + cap_var = model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) + + # Initial constraint: Σ yᵢ = 1 + # (Since no probabilities are fixed yet, all agents have coefficient 1) + sum_constr = model.addConstr( + grb.quicksum(agent_vars[agent_id] for agent_id in people) == 1, + name="sum_unfixed" + ) + + # Initial committee constraints + for committee in committees: + model.addConstr(grb.quicksum(agent_vars[agent_id] for agent_id in committee) <= cap_var) + + # Initial objective: minimize ŷ (cap_var) + model.setObjective(cap_var, grb.GRB.MINIMIZE) + + # Change Gurobi configuration to encourage strictly complementary ("inner") solutions. + model.setParam("Method", 2) # optimize via barrier only + model.setParam("Crossover", 0) # deactivate cross-over + + return model, agent_vars, cap_var, sum_constr + + +def _update_dual_model( + model: Any, + agent_vars: dict[str, Any], + cap_var: Any, + sum_constr: Any, + fixed_probabilities: dict[str, float], + people: People +): + """Updates the persistent dual LP in-place without rebuilding it from scratch.""" + # 1. Update the sum constraint coefficients. + for agent_id in people: + if agent_id in fixed_probabilities: + model.chgCoeff(sum_constr, agent_vars[agent_id], 0.0) + else: + model.chgCoeff(sum_constr, agent_vars[agent_id], 1.0) + + # 2. Update the objective function + model.setObjective( + cap_var - grb.quicksum( + fixed_probabilities[agent_id] * agent_vars[agent_id] + for agent_id in fixed_probabilities + ), + grb.GRB.MINIMIZE + ) + + +def _add_report_update( + report: RunReport, dual_obj: float, value: float, upper: float, committees: set[frozenset[str]] +) -> None: + at_most = f"{dual_obj - upper + value:.2%}" + dual_obj_str = f"{dual_obj:.2%}" + gap_str = f"{value - upper:.2%}" + report.add_message_and_log( + "leximin_is_at_most", + log_level=logging.DEBUG, + at_most=at_most, + dual_obj_str=dual_obj_str, + num_committees=len(committees), + gap_str=gap_str, + ) + + +def _run_leximin_column_generation_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + dual_model: Any, + dual_agent_vars: dict[str, Any], + dual_cap_var: Any, + sum_constr: Any, + committees: set[frozenset[str]], + fixed_probabilities: dict[str, float], + people: People, + reduction_counter: int, + report: RunReport, +) -> tuple[bool, int]: + """Run the column generation inner loop for leximin optimization.""" + while True: + dual_model.optimize() + if dual_model.status != grb.GRB.OPTIMAL: + for agent_id in fixed_probabilities: + # Relax all fixed probabilities by a small constant + fixed_probabilities[agent_id] = max(0.0, fixed_probabilities[agent_id] - 0.0001) + + # Update persistent model instead of rebuilding + _update_dual_model(dual_model, dual_agent_vars, dual_cap_var, sum_constr, fixed_probabilities, people) + + logger.debug(f"Status: {dual_model.status} REDUCE PROBS for {reduction_counter}th time.") + reduction_counter += 1 + continue + + agent_weights = {agent_id: agent_var.x for agent_id, agent_var in dual_agent_vars.items()} + new_committee_solver.set_objective( + solver_sum(agent_weights[agent_id] * agent_vars[agent_id] for agent_id in people), SolverSense.MAXIMIZE + ) + new_committee_solver.optimize() + new_set = ilp_results_to_committee(new_committee_solver, agent_vars) + value = new_committee_solver.get_objective_value() + + upper = dual_cap_var.x + dual_obj = dual_model.objVal + + _add_report_update(report, dual_obj, value, upper, committees) + if value <= upper + EPS: + for agent_id, agent_weight in agent_weights.items(): + if agent_weight > EPS and agent_id not in fixed_probabilities: + fixed_probabilities[agent_id] = max(0, dual_obj) + return True, reduction_counter + + assert new_set not in committees + committees.add(new_set) + # Add the new constraint directly to the persistent model + dual_model.addConstr(grb.quicksum(dual_agent_vars[agent_id] for agent_id in new_set) <= dual_cap_var) + + +def _solve_leximin_primal_for_final_probabilities( + committees: set[frozenset[str]], fixed_probabilities: dict[str, float] +) -> list[float]: + """Solve the final primal problem to get committee probabilities from fixed agent probabilities.""" + primal = grb.Model() + committee_vars = [primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for _ in committees] + eps = primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) + primal.addConstr(grb.quicksum(committee_vars) == 1) + for agent_id, prob in fixed_probabilities.items(): + agent_probability = grb.quicksum( + comm_var for committee, comm_var in zip(committees, committee_vars, strict=False) if agent_id in committee + ) + primal.addConstr(agent_probability >= prob - eps) + primal.setObjective(eps, grb.GRB.MINIMIZE) + primal.optimize() + + probabilities = numpy.array([comm_var.x for comm_var in committee_vars]).clip(0, 1) + return list(probabilities / sum(probabilities)) + + +def _run_leximin_main_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + committees: set[frozenset[str]], + people: People, + report: RunReport, +) -> dict[str, float]: + """Run the main leximin optimization loop that fixes probabilities iteratively.""" + fixed_probabilities: dict[str, float] = {} + reduction_counter = 0 + + # Build the model exactly ONCE + dual_model, dual_agent_vars, dual_cap_var, sum_constr = _setup_dual_model(people, committees) + + while len(fixed_probabilities) < people.count: + logger.debug(f"Fixed {len(fixed_probabilities)}/{people.count} probabilities.") + + # Update the math before running the inner loop + _update_dual_model(dual_model, dual_agent_vars, dual_cap_var, sum_constr, fixed_probabilities, people) + + should_break, reduction_counter = _run_leximin_column_generation_loop( + new_committee_solver, + agent_vars, + dual_model, + dual_agent_vars, + dual_cap_var, + sum_constr, + committees, + fixed_probabilities, + people, + reduction_counter, + report, + ) + if should_break: + break + + return fixed_probabilities + + +def find_distribution_leximin( + features: FeatureCollection, + people: People, + number_people_wanted: int, + check_same_address_columns: list[str], + solver_backend: str = "highspy", +) -> tuple[list[frozenset[str]], list[float], RunReport]: + """Find a distribution over feasible committees that maximizes the minimum probability of an agent being selected.""" + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + report = RunReport() + report.add_message_and_log("using_leximin_algorithm", logging.INFO) + grb.setParam("OutputFlag", 0) + + new_committee_solver, agent_vars = setup_committee_generation( + features, people, number_people_wanted, check_same_address_columns, solver_backend + ) + + committees, covered_agents, initial_report = generate_initial_committees( + new_committee_solver, agent_vars, 3 * people.count + ) + report.add_report(initial_report) + + fixed_probabilities = _run_leximin_main_loop(new_committee_solver, agent_vars, committees, people, report) + + probabilities_normalised = _solve_leximin_primal_for_final_probabilities(committees, fixed_probabilities) + + return list(committees), probabilities_normalised, report \ No newline at end of file diff --git a/src/sortition_algorithms/committee_generation/op2_leximin.py b/src/sortition_algorithms/committee_generation/op2_leximin.py new file mode 100644 index 0000000..5fa2acd --- /dev/null +++ b/src/sortition_algorithms/committee_generation/op2_leximin.py @@ -0,0 +1,273 @@ +# ABOUTME: Leximin algorithm for committee generation. +# ABOUTME: Maximizes minimum probability, breaking ties by second-lowest, third-lowest, etc. +# Option 2: Fixing the "Batch" Heuristic + +import logging +from typing import Any + +import numpy + +try: + import gurobipy as grb + + GUROBI_AVAILABLE = True +except ImportError: + GUROBI_AVAILABLE = False + grb = None + +from sortition_algorithms.committee_generation.common import ( + EPS, + generate_initial_committees, + ilp_results_to_committee, + setup_committee_generation, +) +from sortition_algorithms.committee_generation.solver import Solver, SolverSense, solver_sum +from sortition_algorithms.features import FeatureCollection +from sortition_algorithms.people import People +from sortition_algorithms.utils import RunReport, logger + + +def _dual_leximin_stage( + people: People, + committees: set[frozenset[str]], + fixed_probabilities: dict[str, float], +) -> tuple: + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + assert len(committees) != 0 + + model = grb.Model() + agent_vars = {agent_id: model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for agent_id in people} + cap_var = model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) + model.addConstr( + grb.quicksum(agent_vars[agent_id] for agent_id in people if agent_id not in fixed_probabilities) == 1 + ) + for committee in committees: + model.addConstr(grb.quicksum(agent_vars[agent_id] for agent_id in committee) <= cap_var) + model.setObjective( + cap_var + - grb.quicksum(fixed_probabilities[agent_id] * agent_vars[agent_id] for agent_id in fixed_probabilities), + grb.GRB.MINIMIZE, + ) + + model.setParam("Method", 2) + model.setParam("Crossover", 0) + + return model, agent_vars, cap_var + + +def _add_report_update( + report: RunReport, dual_obj: float, value: float, upper: float, committees: set[frozenset[str]] +) -> None: + at_most = f"{dual_obj - upper + value:.2%}" + dual_obj_str = f"{dual_obj:.2%}" + gap_str = f"{value - upper:.2%}" + report.add_message_and_log( + "leximin_is_at_most", + log_level=logging.DEBUG, + at_most=at_most, + dual_obj_str=dual_obj_str, + num_committees=len(committees), + gap_str=gap_str, + ) + + +def _run_leximin_heuristic_for_additional_committees( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + dual_model: Any, + dual_agent_vars: dict[str, Any], + dual_cap_var: Any, + committees: set[frozenset[str]], + people: People, + agent_weights: dict[str, float], + upper: float, + value: float, +) -> int: + """Run heuristic to find additional committees to reduce outer loop calls.""" + counter = 0 + new_set = None + + for _ in range(10): # Try to harvest up to 10 extra committees + if new_set is not None: + for agent_id in new_set: + # Penalize agents in the last found committee so it finds different ones + agent_weights[agent_id] *= (upper / value) if value > EPS else 0.5 + + sum_weights = sum(agent_weights.values()) + if sum_weights < EPS: + break + for agent_id in agent_weights: + agent_weights[agent_id] /= sum_weights + upper /= sum_weights + + new_committee_solver.set_objective( + solver_sum(agent_weights[agent_id] * agent_vars[agent_id] for agent_id in people), + SolverSense.MAXIMIZE, + ) + new_committee_solver.optimize() + new_set = ilp_results_to_committee(new_committee_solver, agent_vars) + value = sum(agent_weights[agent_id] for agent_id in new_set) + + if value <= upper + EPS or new_set in committees: + break + + committees.add(new_set) + dual_model.addConstr(grb.quicksum(dual_agent_vars[agent_id] for agent_id in new_set) <= dual_cap_var) + counter += 1 + + return counter + + +def _run_leximin_column_generation_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + dual_model: Any, + dual_agent_vars: dict[str, Any], + dual_cap_var: Any, + committees: set[frozenset[str]], + fixed_probabilities: dict[str, float], + people: People, + reduction_counter: int, + report: RunReport, +) -> tuple[bool, int]: + while True: + dual_model.optimize() + if dual_model.status != grb.GRB.OPTIMAL: + for agent_id in fixed_probabilities: + fixed_probabilities[agent_id] = max(0.0, fixed_probabilities[agent_id] - 0.0001) + dual_model, dual_agent_vars, dual_cap_var = _dual_leximin_stage( + people, + committees, + fixed_probabilities, + ) + logger.debug(f"Status: {dual_model.status} REDUCE PROBS for {reduction_counter}th time.") + reduction_counter += 1 + continue + + agent_weights = {agent_id: agent_var.x for agent_id, agent_var in dual_agent_vars.items()} + new_committee_solver.set_objective( + solver_sum(agent_weights[agent_id] * agent_vars[agent_id] for agent_id in people), SolverSense.MAXIMIZE + ) + new_committee_solver.optimize() + new_set = ilp_results_to_committee(new_committee_solver, agent_vars) + value = new_committee_solver.get_objective_value() + + upper = dual_cap_var.x + dual_obj = dual_model.objVal + + _add_report_update(report, dual_obj, value, upper, committees) + if value <= upper + EPS: + for agent_id, agent_weight in agent_weights.items(): + if agent_weight > EPS and agent_id not in fixed_probabilities: + fixed_probabilities[agent_id] = max(0, dual_obj) + return True, reduction_counter + + assert new_set not in committees + committees.add(new_set) + dual_model.addConstr(grb.quicksum(dual_agent_vars[agent_id] for agent_id in new_set) <= dual_cap_var) + + # --- NEW HEURISTIC TRICK INJECTED HERE --- + counter = _run_leximin_heuristic_for_additional_committees( + new_committee_solver, + agent_vars, + dual_model, + dual_agent_vars, + dual_cap_var, + committees, + people, + agent_weights, + upper, + value + ) + if counter > 0: + logger.debug(f"Heuristic added {counter} extra committees.") + + +def _solve_leximin_primal_for_final_probabilities( + committees: set[frozenset[str]], fixed_probabilities: dict[str, float] +) -> list[float]: + primal = grb.Model() + committee_vars = [primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for _ in committees] + eps = primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) + primal.addConstr(grb.quicksum(committee_vars) == 1) + for agent_id, prob in fixed_probabilities.items(): + agent_probability = grb.quicksum( + comm_var for committee, comm_var in zip(committees, committee_vars, strict=False) if agent_id in committee + ) + primal.addConstr(agent_probability >= prob - eps) + primal.setObjective(eps, grb.GRB.MINIMIZE) + primal.optimize() + + probabilities = numpy.array([comm_var.x for comm_var in committee_vars]).clip(0, 1) + return list(probabilities / sum(probabilities)) + + +def _run_leximin_main_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + committees: set[frozenset[str]], + people: People, + report: RunReport, +) -> dict[str, float]: + fixed_probabilities: dict[str, float] = {} + reduction_counter = 0 + + while len(fixed_probabilities) < people.count: + logger.debug(f"Fixed {len(fixed_probabilities)}/{people.count} probabilities.") + + dual_model, dual_agent_vars, dual_cap_var = _dual_leximin_stage( + people, + committees, + fixed_probabilities, + ) + + should_break, reduction_counter = _run_leximin_column_generation_loop( + new_committee_solver, + agent_vars, + dual_model, + dual_agent_vars, + dual_cap_var, + committees, + fixed_probabilities, + people, + reduction_counter, + report, + ) + if should_break: + break + + return fixed_probabilities + + +def find_distribution_leximin( + features: FeatureCollection, + people: People, + number_people_wanted: int, + check_same_address_columns: list[str], + solver_backend: str = "highspy", +) -> tuple[list[frozenset[str]], list[float], RunReport]: + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + report = RunReport() + report.add_message_and_log("using_leximin_algorithm", logging.INFO) + grb.setParam("OutputFlag", 0) + + new_committee_solver, agent_vars = setup_committee_generation( + features, people, number_people_wanted, check_same_address_columns, solver_backend + ) + + committees, covered_agents, initial_report = generate_initial_committees( + new_committee_solver, agent_vars, 3 * people.count + ) + report.add_report(initial_report) + + fixed_probabilities = _run_leximin_main_loop(new_committee_solver, agent_vars, committees, people, report) + + probabilities_normalised = _solve_leximin_primal_for_final_probabilities(committees, fixed_probabilities) + + return list(committees), probabilities_normalised, report \ No newline at end of file diff --git a/src/sortition_algorithms/committee_generation/op4_leximin.py b/src/sortition_algorithms/committee_generation/op4_leximin.py new file mode 100644 index 0000000..6c119e0 --- /dev/null +++ b/src/sortition_algorithms/committee_generation/op4_leximin.py @@ -0,0 +1,213 @@ +# ABOUTME: Leximin algorithm for committee generation. +# ABOUTME: Maximizes minimum probability, breaking ties by second-lowest, third-lowest, etc. +# Option 4: Heuristic Pricing (The "Good Enough" Approach) + +import logging +from typing import Any + +import numpy + +try: + import gurobipy as grb + + GUROBI_AVAILABLE = True +except ImportError: + GUROBI_AVAILABLE = False + grb = None + +from sortition_algorithms.committee_generation.common import ( + EPS, + generate_initial_committees, + ilp_results_to_committee, + setup_committee_generation, +) +from sortition_algorithms.committee_generation.solver import Solver, SolverSense, solver_sum +from sortition_algorithms.features import FeatureCollection +from sortition_algorithms.people import People +from sortition_algorithms.utils import RunReport, logger + + +def _dual_leximin_stage( + people: People, + committees: set[frozenset[str]], + fixed_probabilities: dict[str, float], +) -> tuple: + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + assert len(committees) != 0 + + model = grb.Model() + agent_vars = {agent_id: model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for agent_id in people} + cap_var = model.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) + model.addConstr( + grb.quicksum(agent_vars[agent_id] for agent_id in people if agent_id not in fixed_probabilities) == 1 + ) + for committee in committees: + model.addConstr(grb.quicksum(agent_vars[agent_id] for agent_id in committee) <= cap_var) + model.setObjective( + cap_var + - grb.quicksum(fixed_probabilities[agent_id] * agent_vars[agent_id] for agent_id in fixed_probabilities), + grb.GRB.MINIMIZE, + ) + + # ---> OPTION 4 SPEEDUP: Let Gurobi use its default Simplex method <--- + # We are commenting out Paul's forced Barrier method. + # model.setParam("Method", 2) # optimize via barrier only + # model.setParam("Crossover", 0) # deactivate cross-over + + return model, agent_vars, cap_var + + +def _add_report_update( + report: RunReport, dual_obj: float, value: float, upper: float, committees: set[frozenset[str]] +) -> None: + at_most = f"{dual_obj - upper + value:.2%}" + dual_obj_str = f"{dual_obj:.2%}" + gap_str = f"{value - upper:.2%}" + report.add_message_and_log( + "leximin_is_at_most", + log_level=logging.DEBUG, + at_most=at_most, + dual_obj_str=dual_obj_str, + num_committees=len(committees), + gap_str=gap_str, + ) + + +def _run_leximin_column_generation_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + dual_model: Any, + dual_agent_vars: dict[str, Any], + dual_cap_var: Any, + committees: set[frozenset[str]], + fixed_probabilities: dict[str, float], + people: People, + reduction_counter: int, + report: RunReport, +) -> tuple[bool, int]: + while True: + dual_model.optimize() + if dual_model.status != grb.GRB.OPTIMAL: + for agent_id in fixed_probabilities: + fixed_probabilities[agent_id] = max(0.0, fixed_probabilities[agent_id] - 0.0001) + dual_model, dual_agent_vars, dual_cap_var = _dual_leximin_stage( + people, + committees, + fixed_probabilities, + ) + logger.debug(f"Status: {dual_model.status} REDUCE PROBS for {reduction_counter}th time.") + reduction_counter += 1 + continue + + agent_weights = {agent_id: agent_var.x for agent_id, agent_var in dual_agent_vars.items()} + new_committee_solver.set_objective( + solver_sum(agent_weights[agent_id] * agent_vars[agent_id] for agent_id in people), SolverSense.MAXIMIZE + ) + new_committee_solver.optimize() + new_set = ilp_results_to_committee(new_committee_solver, agent_vars) + value = new_committee_solver.get_objective_value() + + upper = dual_cap_var.x + dual_obj = dual_model.objVal + + _add_report_update(report, dual_obj, value, upper, committees) + if value <= upper + EPS: + for agent_id, agent_weight in agent_weights.items(): + if agent_weight > EPS and agent_id not in fixed_probabilities: + fixed_probabilities[agent_id] = max(0, dual_obj) + return True, reduction_counter + + assert new_set not in committees + committees.add(new_set) + dual_model.addConstr(grb.quicksum(dual_agent_vars[agent_id] for agent_id in new_set) <= dual_cap_var) + + +def _solve_leximin_primal_for_final_probabilities( + committees: set[frozenset[str]], fixed_probabilities: dict[str, float] +) -> list[float]: + primal = grb.Model() + committee_vars = [primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) for _ in committees] + eps = primal.addVar(vtype=grb.GRB.CONTINUOUS, lb=0.0) + primal.addConstr(grb.quicksum(committee_vars) == 1) + for agent_id, prob in fixed_probabilities.items(): + agent_probability = grb.quicksum( + comm_var for committee, comm_var in zip(committees, committee_vars, strict=False) if agent_id in committee + ) + primal.addConstr(agent_probability >= prob - eps) + primal.setObjective(eps, grb.GRB.MINIMIZE) + primal.optimize() + + probabilities = numpy.array([comm_var.x for comm_var in committee_vars]).clip(0, 1) + return list(probabilities / sum(probabilities)) + + +def _run_leximin_main_loop( + new_committee_solver: Solver, + agent_vars: dict[str, Any], + committees: set[frozenset[str]], + people: People, + report: RunReport, +) -> dict[str, float]: + fixed_probabilities: dict[str, float] = {} + reduction_counter = 0 + + while len(fixed_probabilities) < people.count: + logger.debug(f"Fixed {len(fixed_probabilities)}/{people.count} probabilities.") + + dual_model, dual_agent_vars, dual_cap_var = _dual_leximin_stage( + people, + committees, + fixed_probabilities, + ) + + should_break, reduction_counter = _run_leximin_column_generation_loop( + new_committee_solver, + agent_vars, + dual_model, + dual_agent_vars, + dual_cap_var, + committees, + fixed_probabilities, + people, + reduction_counter, + report, + ) + if should_break: + break + + return fixed_probabilities + + +def find_distribution_leximin( + features: FeatureCollection, + people: People, + number_people_wanted: int, + check_same_address_columns: list[str], + solver_backend: str = "highspy", +) -> tuple[list[frozenset[str]], list[float], RunReport]: + if not GUROBI_AVAILABLE: + msg = "Leximin algorithm requires Gurobi solver which is not available" + raise RuntimeError(msg, "gurobi_not_available", {}) + + report = RunReport() + report.add_message_and_log("using_leximin_algorithm", logging.INFO) + grb.setParam("OutputFlag", 0) + + new_committee_solver, agent_vars = setup_committee_generation( + features, people, number_people_wanted, check_same_address_columns, solver_backend + ) + + # We keep our Option 3 win here (capped at 100) + committees, covered_agents, initial_report = generate_initial_committees( + new_committee_solver, agent_vars, 100 + ) + report.add_report(initial_report) + + fixed_probabilities = _run_leximin_main_loop(new_committee_solver, agent_vars, committees, people, report) + + probabilities_normalised = _solve_leximin_primal_for_final_probabilities(committees, fixed_probabilities) + + return list(committees), probabilities_normalised, report \ No newline at end of file