From a4fd141dc823051658e723a6fc9c311f27d37fda Mon Sep 17 00:00:00 2001
From: Giovanni Santacatterina <santacatterinagiovanni@gmail.com>
Date: Mon, 11 Mar 2024 11:24:02 +0100
Subject: [PATCH] Faster variances

---
 .RData              | Bin 0 -> 35768 bytes
 .Rhistory           | 542 +++++++++++++++++++++++++++++++++++++++++---
 R/RcppExports.R     |   8 +
 R/variance.R        |  45 +---
 src/RcppExports.cpp |  33 +++
 src/variance.cpp    |  33 +++
 6 files changed, 587 insertions(+), 74 deletions(-)
 create mode 100644 .RData

diff --git a/.RData b/.RData
new file mode 100644
index 0000000000000000000000000000000000000000..c78527400725f136d248cc4155f516ab89159b97
GIT binary patch
literal 35768
zcmV(zK<2+6iwFP!000001LXaCIF(QPH;%9UT2zwkDxs2uM5L6=P*f^O5=9|}B(e+b
z2-&1ik))!KkWk5fCi{f!N%s4_WxwaUKF{%dj^lSczvnsb?|t9ByWh|6_Xn=Ej=AQV
zIp>_$`Ffo*bJ04WE6c9SP9PA-1PYl%AX8bF1nMbmwQaiz1PY1ufk>bb*#4`4%|%<*
z@74dKEP=p5AkfPoJwViRrrsM(Zkm1U_ox#U4eq=Z_^Aa_dAs)IGghG=j?ogP35JlW
z9dfpeyB(5mC7sV*vmH}?x947Krlb11DXt55%21bTgqOtCqhK^spflymJuua9=#A&4
zY&5EV>wDZbX-NKNej&*90w$);F3Wqi0utQs+}Jep8B98K-(_;!fa#g_J;Bwt!9?+E
z=N_jtOgTQE{<T98&E+a>8JUs5RPU-USE~jf{j%reGJ|2XG*WWOdb24euYGQ5!tDg`
z6E7aWFeW75)9>KhkcaxL)s<7&BoT9xlVYt{gl2749G-k&gXwOI4ntz`s7mOA;d%}X
ziPGj$Ll>`s+DY|GBRvezzL7h6zDyG29?EdM03oENHhBbiKZdmHq#$Ljm5_Y=pl7AK
zEU2eTe`$UOA>jtw+@TdtkTCRowT{<LNbfoFW++hs%q4$-%ne+i&+1h!F4x4Q*Pr!|
z=H5iL0x$QO?G{7*As-w=@A`uDnM>!kF$rjLtsZ~qIxkdv@_QZINdn3~v${ZVUKsWF
zALqZ^Xa+_nTz%G=1wn$F#Pq2kWz<Q0#Xk0?3{&owihB(uW6}keM+O!(U_hRY%iwA*
zm}ySZ=t^5e9n>@Dk7$fz;u`gillN2*^QLi~kE9ExDh8ifh?NJ-U$1<U-}8d8y)Qtd
zYCI<1{M_|)DHz?_*jWBRnt-|kYL3qy+Kno<-&q?e6b8xaQ6#i93>lszzh@?A(L8;l
zK-AenOv<{TXLisVGuR_GdLP>h7VW04-p@XVKIAZ=OanxVo(q+vC)Q|GqDf2CbQCQC
zV~^7&Yf=A*!JBVphd`_CRddUj2(;iK${%&C57UD8WxeS00KM4*4T~BWbda+&E1wB~
z_NO5?-HL<3)F#DDiEBI1<n;ny!w2n{<mWoKYu`ppa6?;GE$V{lHJ3Pj<)a|!K4=?y
zc^{2Q1@b+}ECzXX<UN<<G{IDlfmKlAZIE1RQ6P1p1rp9!a;$&w0<*E5THJkXGnz@e
zvG>fGHb^ncl+9e#j+Xq7%d)FiW70??y-4K@7+9n1_$G@2DK;CNmU8x@#Tt`wCzTw~
ze<I$o%!C`vT1RoOKh*~*S1xMWmaha4Wb-bnPzac6EbwVUQU^@B?6BZ$=t6Txx32%F
zaT?OKD$)vMe}S)OEHsllnV5D<=A4w43!1H|&3tBc47CS1SetLM#?)gX;h(=BLe&ws
z+a*NIP`{ST$Ctb!Xnf?WOs)|hihW3t`Cy8{sG*nm`#Nz*duS{kKMo+(H;I!O*NtX`
zIq>1#!<e$7A@_8-8=B16vT;gT9TLZsuGz1gz-)9!wgp8i)TXANYc3-O@_0|2lf7(>
z8jW?&@{UYHLaDKIV>A^l-C+Nzwe1ol!=I;mBsQZahgEy$R5USdr(%y*c@y~A6g2uO
zlz|qX-2V2UK?9APd3@kg2??^T-CWXDvl`RSZM0vOt%4c00kdUX)|e<E{Lr8K3YxLK
z$hA*9A2U8XE@vA)3~F{b#Vg76pa$6oFW0Tz4Hk<!Uz#}Y0Zio%nLefinnDjO=A*PR
zk$)l6v+^Dqp*$)BD`}7>ZglW-#~sY**nV1YxfL2L^4(!n<c?+o)`?I&W*{47a@}2=
z2ui|;uI0i6NZoK&-I{8N<^+N%sMs2fyFZNblRkzzJ=R!Nn1_N97rCwLgK3bozOnb5
zu>_cD`qf9SsK>;pb0_;E%}~c=;zfU*X-Eq_5s_DX0Zs1OsJR2Z1MTsX*r7WDEe$@p
zayd;H)4u9ecbW}A@(1(xUEXg{XWSQA8{!i%D!uvR8g4E$HMO)v%en&Ux9&^%VZR%a
zd&<t&O?p7;S@OvhvqGSkJK;%xs6U#j5jz!pL=(-q9pjc{Q$gR7iRiMQHR_BWFk7Z7
z27W%+ICi!=3{^g8Jf4wNh6#g1mpCqMM}v9JPxgd(LE^!oerwVf&}mB736O)JScW8;
z^&kxW>PyLMN@Zet<b!<~_3I#GeUx(Tjvt^u&meE_Z4uNtD}6`qsXiv?hz&8+rNGye
z-T8U}oS=NmAp18iY0&m#D%pN$8j=Q^Z%7sOV!CH*kG;DJ7|Sr3Tvt>I=4$lzUOplP
z8oz3rwhS^*wa~18VLz*U_|+7d)z82}#jH)1dlVS+vn=t*TEtX8YEl3GXf&9d!GGk1
zHKeZ=E+00r1<pmgS1C3Sz|i_1!#ti86mJoqv1$>6iR>FO+KW<<x+ETO;NTAscxk+r
z!`cOmc1!HxxGWE8R|}q)ULQlP#l(vIhzu}yM0S<4G8bwayBtO}I||8W?m8d+wnE}L
zZ|DB^MrcN6r_u5gU6?VJ+Hj6*C1l)b^{OeiM4zq=+C?{F)VN17)7HoZ%x80W5%q&m
zmBU1-=2>e{1LU&Zj<Wpq#eU_A9u-#H^gOh6r4^_QNaa)Ag+b4p+UTY<O*EQhB68Vf
z3Q~`ze!Xcw1!ivsu8e9L0$uOY!&1&N&;;MvzL%!TkV^WV<z&AM6m1sl%id!MW-NkN
z@>*wsdEr2p<*xFOaTFB}y_-YDLiKLfk8?tXQR?O>O+HX}M%uM&x(szZP`bVQMl@Pb
z=-9R_%MXksr+3$48VmQWdR$GCz_99G6@}|^n3y8Rp)bG(2HuvR&@suygca*z)NlSm
zvm0Yid_KyJi6u&Qut5$r>b-muDCz<kXGYg^aZZ2%Q9boLUD}Wga+>xgK1Bu84{F@;
z!e}H<Ud%`5E*j;#wk`apI2hq^4LSL95|iuR8J7!2ppra;Po_~VXh@a<{rcP?UFN~{
zID>XbNlv7UbOwNV!nzOjKlVU2IXgNn>;|UK&TAU*T|+r5mgd|QI5E{TI_W!C2I!-j
zti8=<hp7dg8af`~m^Ap|J9)!5OuAFDGE$lYjft$Evh3>z<5KbC_DT0Z%WBKPcP7!G
zWz3-BewGPn9~I0OSfLN5-|bymv3?N}Jyvb3MDH=d{49@3%n(>G+&Px;;Sr?DZ#K}M
zx()&xb{anJB0xfiMMLT69Z31uxYbtk7V0(;X}!8d4i$U}>;DcDF!j{;b&r?xqn7f8
z>1A&Ss3l*vxUI<q%*4K{(_O`y?|O&1a+P;u5=BbkQ>!)@pRn1HcBTz=?wJ!_XAl7i
zCQSuiG2778g-@!<d+Q<b$a%IAaUv$XRkpg#-ULSUmf40|W`pML+WRBF!XW)gTs?Pp
zKV}o0I~K|Q1msq!loO9aR5Nugb>}=p6{br=TYZHw^^K;mS7;5Y6nJ(`Q~nSbaUJmb
z{HhI7J$~fAT6qxC@7zBm>ZyuK_$lk($!J1G@jU0#rtOewU1liwY&9zX-11BEh&7la
z8ufWKrGd$kfZZF@uVPXh;H)p4MMLf<t@j;21X^6QUv?k43F_4%=+1A&Kvz%Qhm^WP
zv{+eXFxGnoeeHMSVoPA6{?nR9UK$yYAbVYN_TD?t;hE&|nVSJdKi`*1ka+~=$5X8<
z1>zxNtKC)qp6!^HvUAcu=Kv^Zu-Khq%ZbSbHSgUOiZJbT_UFm^U^G;8XOeNJ0nLqG
zn%We82-4lYn0=WJM~!L>@}nU!NYXyJXeDfmN#8H^3`;XG^|9gWu%BU=K|aM6^8!G6
zU8tj+-YL{OBJQ;6L<gEdmf1GNeVDv*Ts!-y6q=*v-~als8`X8YDzBVN!SrYa=klUJ
z)cp9URUV-X%pZSZrRN$1I${oXNk32o<KIG3Lq9rTVzm2J3yv;KJ2cIE%#?zr+h2vK
zOwFUdeR4hr*aRT)TYjpH(Pl`IuyIe@Q2>@onCWbRKA^q|1+Tu_gJx6q@V8m6!i<b9
z?*?yN0HfDe=Cjwm1Y;Jf)TZw300X__y-uR~U_3ut`O@?~^!}w^%bXho(=;Bl#y}}d
z7;8Kd8(@Rx;%tQ1?!S%6vHS~x28Y4i<;|=02W>EMz0F06S~?^SDYc#5e*{#RC_GSe
zaKyxzktcrLH!x8t@!`miF-RE`9&lK78ezXe_M`p;(AAO{Qt>1Oj5z;Dw<*jAbB!Vo
z6T}5k+I2^V@YNWSME3o{Z&sk`bnkMp%PYXd-n@fT*E!Iq#dq7|`(40vcjheJcn=Gw
z4?gC7T7;I$9!&RzUBhI#wJJ(U(V##!>gOw&7SzGMY|PX002mE@Ps=5#V;WDwFs0ZJ
zMHK(qpwI^~Nih1B$fjX1J(ZqH+}{P31RdU=-@pqQ5jhn-?^*G<e%<O-%16QY5t6H+
zCp*ZgJg@H_sR`y3A6(n(F$5V+$L&{i_JX<o)T6hP&Va=L?}#4A;&DTIj|826K-!rr
z7nh66K#I@vv%@CK(CqYLx9t`K)OWhqLuQjTB<YQH@+1K?;!~`hc19d6u3mPX=_JdV
z$BWB;F3h1hamkjP?;;?7ne^5(SD2vUdCGEbD*|Lw61@~ue-e@!Y9BunzJ=+IS9i7F
zzJnG%xXPS~TZQt&yK4IQSv<)fxkXK1#5C@VjWctjn6Wd%S(xJjq+|%2r_OeR_R8xM
zrW%JJ<AYwy-1ha5c%kjrrn2{toaY(4!|Vj4^pVS)mHW};-Y+t0$-JP<O}jBNYZ?>l
z*6?<otw!T~LhJUn{K7PI{v?+n7QYWS?r#+$fr`o2pA5?~(4#letbaxfQvDWt@+4cq
z99QGh8Pz^CENi$;`@jGuS-=1KWOx~9trnL3K-5QLx@i;73$oCn$Tf3BNk;jrU${kg
zU^McyUQcVaB$`@BX4^P-956?>eR(O*0V(I~O1t~`K#fs5ZJt|~)&84}3zN^l)WV)w
zPUc3y%q~wfd9V@F)IB}<BUgjDg}v85FNk82O!3#eJ4ukbyl>^RPB%>Ox}L*!;sQ%&
z%@~~BZHHz$cfUJlaRl}7$8B-8H^Wr!8S7R<8Pw$xxKLyxiE=M&vXD#W1YL5uiQ&2j
z(A1^_dyWpD!xVgUeo%iom>(7;x7hCiv#2Z3K!yW+yOom>{mu~0|5Dkltq_FCJ(Dq-
zsa}AYU#7aUjRYn%JTPw3Va0u&l9GK}F;REtrZwX)Kx?qOjBSG?>N|BfVuBG330rUO
zwb$u|^o>BV=afB~t4JRG#OH#>^Cbi(%$H%dT}^qqSPl|z-aPfq_64MQWC<|N5>fNS
z@gcz&UBnc&s@>gr7i94;uD#jef{7jnKQ1I+Mf3Zd9`BZZ0!A!PygV-@2#V_N$T8;2
zz?gO2<`-UjP|glMF=bhd7Nn3(%v)(N7}Di_SZE#%@xE<jBQYVhUzDxY`3ojkH#qq`
zUx%2I$4?%<76cN(yy%=d0TbEZ+c~M+f<&FA^s(2cP;^k7peli-mmB@|(76lHWNzO{
z=F@yk2`c@5)%7Myi>iIF#-at%pX)b`$Zo>qS#z`Mgf>VAS^XbOWii>*pu0vz8!ag{
z+e@vh0t;;S>85+5Fyp}!?bnW~kb1*YGxzZtG`LE6EcCM%q|bza5^<J}u*f)X@oOs@
zpNvTzUTuR0XG-@<md=Cmmu4YZX~}3ha`M#e4O>9n^h<^0{AtJ#OFvL{cN=0pX*$AR
z<PK(*U1&%%dX0wlq+ZwtAHo!&b;pjL-GwGh&p%^#e~hUW4=hTY7+|{d(cwakb1Xc#
ze||6DLs0f`%4ltkKbk1@ce9GR1eo>ut;)+Qz~aFN5BSrEK>t2Ovx$6y`tP}N|5_4(
zG=atI2kx8(GsEY@4u?JjHNCs$y0uf#)Q@tnmc*T4scVn19eXku61XRAw@Cz3?H=>&
zdL{%JYoBGFPdW_hPJX%bmOwzi)a`GUlX@`mON<LGs|NJ6`=qe#@WO;EqMuxr)X^A+
zT>r{$9WY0T62Z5-P`j4Z)9vf@!6OuGAo7%eI<?LP+c`;NHlAZy+u2Ic%$54v)!nI>
z_{g-`cTo~eXvn{KC71!CS~7Q!{Di2j<F=MrRtuW^ocC4#M<S{XvtDc9YYnNpV~)g!
z3xTd-t2<xYrBR;fsNO9v3@QFkN__kZ&?l4aFQ4!d(3F<K$u}wjnBJ76Q^MR2=~BXP
zggKT#IoHmWM(2c4-&@TETU%W)Z#$A>M^`~b?HzlSO<KSYE8j;X4{1pBbTcM&e*nz%
z?}G&oU!rlN`pug}--8cjU+<oJfYDgkNPBz;36;`#o=YjUVDcx^gfvGLH2a1-mt9H%
zP3*c*SRlfNMt@M;wKq+Gaal=9%rhU*M(kVCJy3-y{(BaK3)8@qsC4g)fIKAaH$S!6
zZ5x`7Kd?DDwGtDyrpB8tE{Aj#>r#&st)M4mlbwnyi@$;D8~STXAw8K-TZEGeiBW5w
zZ*yG(M*Buxma^r+PyJ1o4j$f(8C`LQ*n|VXNa+}t+jlK6Cw?>%nZ`r<ZEN-WYww`Z
z+sDsWXD$cRo<YNF+;q_Vhxi><=j9-E>upoJM{eNf-D%{x)f%)*FZ!QMXT?2TS%XVp
z5rs>Ne%JMaDC3@Mp<I$78a1Z2yeRZV1D~&#^Ypl&B>8=p>mLZS>UUDglhXmcE2$!`
zI~7^+%7|^VaT=!I9UrjrJdgT?ON;~)3qhX3k+<3>3_)_uHrvB8tT-x@9+we%3X?B|
zSLZ}uh4fXWEvE}r(PG_?_!p~aEZ!llBl|21F@KO09@=a~6G>bmyxzu`;bDoN@w~#+
zTE$zPr~;A#?$lgy8ODqX2|eEj7)%Zd{itoLfNV0niD%ty(PFK1W&2MF@X`FjmROD!
zw0L}~|HPRtNH$tMme6Pc8Dp)GJ75Gcvy%fXH(f$~&x<#i$KQlh{;<c>W#2IIeRAZU
z$S6pbVNxYBkATpdAtS$bxnNoy0zLXdXigEV-&6Mm^hF57+&`BB8Y^mM_xpu`-j?RV
z=z=GZk$a`{%u8=bVqYou^6_cZ8kFk8v!4~0vO?>_H@`<coFl0-8~HI^=nd&qMIR);
zdmL`sYyfGaTU4UbP0_?-KNXwMZZyACbz!wy943GKto3Tb3bRSXB${8`4@URIPEV~J
zMM)HD)_|NaSiCXp&|%2|3Dz68anx=`eI@7XevbG<a=EJg*|Hokf53Z>qD2tO=WF;{
zugr<YCd!v?U3>(|>}~Y+`d?sfUv<3w&lt$KWH_32(}{(jTfZ&|tpYQ*uH|nT8il0W
zw0%>8-yyB;f#T~Qo|s}>D4#mL8qL`S4HUn759xV&y9<0|z;JPuhM9i`CM!6HZcWU>
zBzCL7PStqOe{bVwj$cnO-EaMmtYvR8WzM(g<Kvg8+)Z&?DqR;%$g4jd;8#H9(&j}?
z3y09S&9}lZdM(O(GqDFPmqxR;g%WqFgh7X7Z;Fqw38oK<E(V{_1<lJ+kFUE@h{l;?
z6zQ^SsMq*~|3d#f8cQ|E_tOeM?UgaNZpCn;k)G_%*CcV4KhSQfG(t?(E2_%k4?<n%
zXzS8d2w;puJuIh52aH_1-+W1%rF$I}1{w>#pb0+%+Q#)Ys8(#&07M3Zxx-JD8{{ve
zk$s6}r>+SC=6Kb_*GFqvdbiY~wL=Rn_(@J)C?R9onoaLre>tGVa^26png=2Anf}UU
zL2DpYjcsyBkOy_Gq0B|xzJjJ+80GwW?T;p}J{z!%T#lxnZj%~Hz@XM(BxA3G7@AtO
zi9SGi0qQuL<u)C1!F1<?dSL!0nmOBcZCS})FyepZLZ^%x$~rB{Pw2#uE_uRPkQRuT
zT01YkSrrGV9r|%EA3X&#F0rya8a9w%#lO|;=@HP*r~jtLybN?XnBU6Jjs|r;2N&xt
z>oBn>+%_V49!*@x)GHgA#pGA$htZS>7+E~4EWcMBEM?!zi#zfX<ZaRjY8w^?^FLHy
z1yJvyYIobTXoEUP{-`Q-NZk^>=DjWMLxSjYu9F$x1(v;|C(1h$GeFa+Z>A#=!kCbG
zWAmv3bCkyyDf>#0fSB&O&!cWnVm3KnaW?`9bv`LsCfg^CT6LV(sXcyzY5i}DPHkub
zrEfoOd0EXwQLDDSba8>`wcjpc-d%|4QOnwPbUQ)Hj|h2%(E>Ehe!PiGMg{c5`b01_
ze8GI;%eBk3(lA*`E-N;u8q9MQf2_S)im8n$mh$^Jk&9o-rC&6b?)Erz*>lZl5MMBj
z8wXubCpgetzP<;{=qw>wSv5=^RS`WH83g7ZoI1KHTm#dWJqes57eiuW#xL)?C%`n1
zG_JgL6%*6LSB;)$gJj9)LQekAF*U2HdNTSPq<zB=4*c>2lZ5i8H9!iCy|KR?B`FG~
zH1}uB9bwr!`v!Mb_$t)&Eb~mQunSl)3ty5hNCMMdu_5xY&Y*4AMK|jpjOKS)4(#eL
z#)SFWtIETisLMOX*>$@YrWs`SYpxi<B$4}qyJ`+Wa`INrYU&Y4c$se%{cANCo01w{
zzhx9M*3sNdc{<V3?59tu%>`(_`A*MjF?BE>wQgy5k};&Xn;!M-lfjgPghE~h15KYf
zp4GE=Eu?=OUlU#81kzuoaPoH&(5O+lYgSSXT6m*>Z)@^5$oLY+5H)7$jeP3K?W=9k
zaC%vz$*>5EzklSMIzUD337+8%Zi?uKZQ$(EMolob4Qsdb$b(^#541y3p=j2^Wtl+O
z4@h;KubWW#!Sbj3cij|zqQ&>$Y2Cs7plLM1AW@yAt7;q{Zr!Q@x`(&+jB^Q~=~p-W
z6;_{v^eBGMnMx0|RQtrfpT`Z%@VZ}C?o|hq9VXX|KB;43A}^z&ou!kf7QfgetAI+A
zGp744r+{9{mCN&T(xAoQtI_Rz2nH?ejnZ@mK;Jq}=~uW8P50YZM?ZcFMjRUyO{HE#
zdP%8_#m5OyD>PqYv91N=-N~=xeQ6Ep5f@|p%u7Ibd1iaMtUj7{|9rMnrw%RFKdP&j
z_zCIai$}WpHe<@Zj~9$3y1>}Ss~QR4_JEnImyUXC+=Wzii@~G(nV4kL@urq$3`t3X
z%&hbM=+n6m2WC%&VMh1u0$Q*Q8ldEz9MX_OQ^wtv2ZV<)^+#MdiSGs&`rLE%;R=1!
za9NosHkkrxd#Bc|*=Y~%WvqMM7DYg@>zYG9^bt@?^7(4v)&rPeuz243ls5Vq9FjNB
z#)M>>$rMec6JUzNKd!i}2vVmNwV$i?qB3)(ov)8G!Q51zKaElf$)%HKcUo7Y5h_O*
zlPQg=X}x8W?+Vf2=}#HfBTvA*@{xw&B@IZky0*T#={Z_Tf(@}59AISXT2^oEYtUNy
zpqD+}5p~u2@9J_6M$D$P=VE$VkSfIMwBb}NsL!_Idhsg(Q#Nu;uke2lCe6Ma3Vg<a
zY6b3f9SSf+OS?N?m^@B{bdygP#T%}n4xQA4+0nk3UbHWy$?r87(f=?ZCbI`rFPAp5
zY)ydluaAOO*Al_#C7yo0nj|zgw}166bPJM~U83#$nThNY<nQz(5HM9k(`M*{DrgtE
z9(4b#7$&S)#bM8v4CdDgABfSmfm9FmwM!3KIwQcCZ-n;>noA1Z@!YQ;%^B=hQY=+K
zm9DBnTD)SIH0A9V$Iu1LK<@RM7L&oS05;oF+JOoA;p<wT#zWF4zdcD~ji^QPs#<QP
zBVb;tDAlJnvGhF6TE_DYSbTZ0>>g(g>J5CWpn3~JdgSHktV`EHcZQbR+|gW?pPcnO
z=oX76U+?sHJv|88lZBU76LQgv*{!caF?lQ=&h=sAbPA@P5OlQpI*G|h_hGznJ*I|u
zKHl)S1rz7L%5ZqI^xx*n@OzS9Fr!?kevfx9COvslNt<JU4y_a0g+{C)ts6g8<8DSx
zTX$>}32=dI+*T)EdL2V8VR5fU9+yG7&J}*)PwOyYq;iLin=?uoPJZZ>5(c`xGgEh4
zg@TdcyW|TQk1*wp-&^S4j5>cDkC^n0LmBrse7=vvP@9bH+5D@~XiVa~v$R<(YL(k?
zeHEJ`Cixaf(%(D+GY@m4p8D)z;qI%CJB*?+wI#;t>ft^#Ni7ITtks0{=)`B2T=>x3
zHNW&&i+)tj-i8Jav-Fl$<iQiR1WcOV8FH&V7qoo#yIv*j1x6D4mVJ|C*?Hic@b|e5
z;2UogGEdfo<n41;S}z^O#JTLq?R8>c;UqCLZHyQFIv4qE)6gL_J+QmnBGm>k4WBJG
zyv|4Mj|2-gHE4qP%;8Zu;{uk{<u=@Rsz;+@z7`S|BuqOte^PUw3mA;y-+H=P5i~Cp
z708b<z+miCB5j8NnsHtdULvbv!rE;woI30<f$ii>X+$B&xMBUVa+M2O9P>L?>-8KB
z<u-dR8FYYw)}2D1hK{0<SLRN?v|eJ0Q>o*C)K)NX!Mow@5q)q@qG2Y%i-5Y7X}f2l
zgu$G{vkuvUAxP(!Y~_@31;c5s1B&VcXym0>hwX_`v{-HKvoF&X^}kXN^Y9(V)WAx~
zp11y>e}emdf1w?EpA~(1-!&GW=z&j;+rLLkK77$!`&r|m{2}$;DlQbXw!xO?EddNh
z?pJxNsfY%L#?|b%F<JAkLGNhd5M(?Oda(;wp&Y)4Rg*fbak2!r_pP!5ZC@VGzY9G8
zDMh<D%*MWe5w(Qqj87Xe^&@lV)T&#Ud_MS{-??}&SAVN0(ToE%&1nt3i!Xze#;~Ra
z+X*xkzKXExYYQgn?kkjD{tk^OazveLeTnH_86MG3crkH>($M?QU6{1fPVn)X08|;1
zx84IdqKVNfX9B{<LCebaZLOOf(agcYN7j+TXy({QfiuR%ko47a&)`XK%$S#C({*b{
zqfNzS{dYLPpm}K@+z|@;I^W4Wr<kBYaJgX1%`#A}klBbLIZ<(hqnLZ0HRwGyo@{X}
z8ZkZIOIDSiLo=9ti?aR_m|B;9F+;TiEQIE1w=1Y%MjP7I#;*?+dn@XHZp;P?Od0Wv
z#Y<p9_{fkZp$jb852gxS&PPMLGMl1g3em?kLycd}TTuChTvZX@I!HLCceo{G9qL+*
zY^*Y(A<-$eFt~;f^*+{r@J50OiQhJhSf2O*86TbH3ZJP$ddIh#y;pkCFTE61uR<?S
zA6?>|xj77VZnAwy2nj}$hgW_}efkCD!^;_;&l`dnu2q_rE~6kTnQx_ju`mn2$E7O0
z-(Yg6z%9kh49rH)bdB-}f#lg|XG6~%L%n?T(%n2&m>55K(&R@Sny>5LI!t;CX5OjU
z$wuFXghyYz7K87g4la)dbt@J=Sq@molu;nvD2RXQuodX|<hxzMIv7)E$>DQSJJ6_<
z=q)<OLrkADw=f*DMe#r6l{7X(G<nVXVI-p#F|R^_Cz@d_{F>dEJ<6I7IgZWDp$tsF
zeQI9KT>vefnDhB!y8@H?N-M+*!ofsI-}AV>9jJ2qP3>pZ1W2m4-E~>&CnVd)(LP;z
z3~8>EQ;ZiYG0{Hhlxo`s@M$66t=$-cY;;ujrXve?A4p97v}E;@Cd6%l!yx&<zQ!fZ
z?_e?X*2v=|K}g7|%YB#s7&92p?3<4*N8?-K-Jghwfr_D!rqezwU2?g2UtD%HCT@k!
zWyPYHh=fGmgv9_KIZ>vEKLJx}NML!H8kh{(OS-C|1qoxF_jwFtz~}8+`}b{tEFLpe
zDk6wMabRQ75$iB8!Dp%DGEspltvi)xPt<{&6Hgy;^7*6o?Exl-3g;myEjneHN;jA_
z9xwEk;J`GQ6lD)Jgr+Krn@N)!A#p>qz-HHONF4bPpmX^;q)hWJd+leA#%cH5?tQuj
zhVEb8BxKF9r&#)Bp1XEv>`}*6Jy*o?+lINItRl4FWTn)Z(uj#ttF(+%%E2ft>D~iS
z14-}p^R{rl0=<Y+j;U6N!hU{eexm3C=`UZ7+9u0_QAU~g+Wl?dSAx;k%h$ue+k}-n
zM>Q~D>a+7)-mZ?xZAINVAsVP>Q+w;O7c74WZ%DjXbRN?U-|<L?4xsFUc0X)wjS0MF
z_vfd*K_`!|KjB>!Ch@yWTQz?KnY$eN-aI!%rOwZpp2wsiVNc-&YRV1Jmie{aX@4jt
z=L&&qQ`2Dn1=ZiCe-RUSWYpIz`eQneTWOEI92nn`zcwL63*a7xwO*7oV&-cRYY!Hp
zD7kfLjLh;IVgJ&T%BvvbWpi+{?sYKT7-g!_lnp6dxh@x%wxd3u09%2V4QN<xmxYn2
zCg|QOnOgr1qt0-9Dk<C_6adC;!(vXr+&=TD#g_w32Swf6vMmFWc5ZY^`dA5=_Is>P
zJZIrnU1<KBq;r^XBwDv<*%j0p@hFbi_z5grDSj^VqKl=^IM!J5K@j(#;zWzC3npbB
zwzZJf2Xo=wO~gyinC%G1&6i?zC|~u}mT#q;n4<V#yuqFa6Ko2@cD+4>$s=x3KZ``b
znD>B%)3`C{vHQMN!I<?qASvbK`FO<S$lpIQ7Y!E8Z+M;42?CQIq9-P2+Caqq8y9A_
z5+EB-yFy3p0-7ysG9$|?g53G=^uU9hVCiucWw2on8kZX|Rwm4Wo_P?yo>2!{R`}#t
z<S@a2#@Q`ze#T)!(fIC5rtD}q;!@kaSJo)&?$G**1&j$7CdH1UT$F2h^OW>`Ye*4N
zSe<wF2<p&EUvucw6V(6p2ARBV8ceGl<xK1Mfs`9hblZw8!1$qLoBXA0G*EKn{pCyw
zm~X0s^{ekeD*dCAE=M3HylQUwUM`HftjiK?oCnZ2LwWC5KSHC?*Ut6IH=$ac>g$>M
zoS>IF@ivu?F=J6UL1m!}b!Y7g7_1uri!zi?U)btV%2|cLU863bx$Jw<kc2R3J+$oT
zHd{&5;}c`GY+nN=tqwCZlsO5Oaw7+pytpv2e!L}b2M6dnqC;CV@(GM@AdrvTjDn<^
zfZ;=p#~{OF)FktEIrv_o9{9?f2?n(cb%Vn*z?`bX@~nCXFtRtJJ4Ejhrh=jlu?v=9
zGAs3f;gK3hu_T-gl=TCt?Gtaa3s^X4tas#F+&t>{l(zVx;)VLck{&Z6azWLvtY)L#
zDJ;B=-6pYV6)J8wj>^d31oNbb@zgeUP}`S$#w~~Sd9Y2PfffXMcs#ysiMxp@%eN11
zusIKD6Yk;i@wT8ulQ-}}Iw$I@W841P1cG_ry|0#yUc!ubCC@0Ni>QrT^{0wf9s0hg
zoV;4z8jLPJyzFqf33UGQ2!8D7#p36r?^-YUFyrm(&CBcZ!069Q{W>rkjH^^AtMzn&
z=3R_#32A;7?rk{Xws8WCM3zYrBjhpNx=i80=4r5?z~2#8laJ}X_X}x0pF!R8n!JE8
zmY!-?7he3P%<8A(i-FNXNOI}ab?5$trj?8J%B!7GCEHE4OQ)tG$+G2RvK$*&(hSNT
zn#{t)=65p(1Oz~T{m`~sz8#QUNbH)+;XzZJ3Y?@KF3@L^bje*Y6!mL#s8(w5L&_#2
z3HzdpnBkYaZLoa>Vs^HM@V@8(3-(?=){VD8+A=HcI$u{br?p<@c1R(Z5XyQ9_3oiz
zQuz1R;f81`>GToj`y@zRW3JoyNdZ!QZWWB~Oh<Fe&!04A={c(H;57R#8%WO;TpmUE
z4t~C2RD|)caPBlaw_|k>nm5@~oGSAY5+Yykj%h6c<5E+i-z?9d=8@?J_d{WjsQg$h
z`Q#L8lHfDgvbGFUcJ6%F?9&N)`a<?`)XYQr4_lt@q!Lun`eny{Z-3BdUgHOXQc$a8
z(k`j4Do8m}OV6{p2&$a<(u+;Kz-Y)?#enC$DF4{naBn$bkpJ~<(#d_okRId~xoyf6
z4Sb?EX>D=GMAF79!6VlJvp04)vH3O_P|5e&I@SOgTXO5X1G_Q#`<a%Badz}(Pt_sc
zc8HqJnaNvLJA$S-Q@3-QF(w;*d2z?H58*ch&NV_9OxLNbG>X+?=@^X>@!MjUtW$41
z^JE0|y;%9EXT%H>$C$rfj8}p1t%clu{up#*`xl<gz6tsUiq>ik`9i8hXTjB1Kf%J!
znS1#s>cQf{y`T0uvhbAi+W57;b6~9ZncRbbwP@acC?M8r5Hm;<5!=6Ufq{Dj^_s~l
zP}R0@z4m?r>aX=?ubr1epA(Hjb1{q8Yc#}dd!vSGZ<#NgW%0(v;mbnNt=?eF?Xait
zj0z@%m~>0siNWMsL0QA4Y+zO{$X)t%GWunEuYY-^3nT;wC`fj&V#wFs?RC|Mz;x63
zooC#Kz|ZBO{&Js!AuVH9!tM*}(8%SiI*%?+)KcJ{&@rTsdWQQB7R)C=(hPT7^_v{j
z^YV<L2r(4Xf{$H0Xf%sS*Ai7e&t^eJvqHg5sV!*m(%6~!;}}x2#RrnEcB2LDH5%VV
zYe37_eHGd#<xthvF4KC;7O+GSw0hHJjH+*K(hS|qL}3zinY5P>%~Xm!d_MCEP1hOs
zv+vIajkP>)sxF41rS9Fj>#aJ#ShAvW<%J3`KPMPwR=E+=BJz2vN{&Jj&njJox@It)
zt<Qh4t_bDZ%sZB9wSab2ACgA&0Gg0mE>o~n0owS|L`}?=gN}DKS)0XfLTYqwf5ngQ
zkaB|8^yy+dm~yzwwNQ2sWS`|dx34w<(>!E$_FP<p%6SVb?N$8INZx&$*7RabI-ArH
zSv-Qey<fy7^!kIL70(6*Z~3Cu^*NUbqIF;?Hm<>j$puNlozkZwThX}L`V7HUUm*45
zsras;5lAfOSXR>B52+~=qU?v9K;p$r(zG!FRCiu{CpSI~Da)T2KM>uI8vUhyg*WVm
zjDq~3qh4;1ZJ&HcQ9~Z4DK>r>OJl{=w?@Yk-;00(ttwj&RRXXjW`~L&B|u7GpGWFw
z048Zhbno4mg1VTa)n8n5(KLZQ=jmzy=1qGh0$MmwOGCxSp229)-kI9@{R}rI@x6f7
zQ$*BhUiVARKNk~De$%-*^BOd{<(<xLm<G)Ue5AIMT~JpUr%v%ShU9@w2TZq@qtTP>
zFXQNTXpvzo-}zb_670_1w~SQ4Y-@$h0}e4TU7qB5LgXx9o-Qe<TU`X0!!6sGpOQgA
zx6TrEh7;4ouXkTr#=?=OVQL2)N6?sgXI1DRH>zncFOg_;1atg-rpw4`XzBSrYFoYn
zCJ1>QLR?yyzTbgeqPP!Jwd<|CX(o_*TrE)DI}<&*SZETPN&sbgk~FoQb!ZZF@f3%;
zKuUmWCM9+PO_TQ;Fa5XxMx+zPHhtxS1YX;-hWabP;EJ)h>=(X}aE+fwAnPDtR=z&K
z+$xINYU|Itay0-ZpZ1|n;ZZcsX4oim&K1+8wGU{COh9t<JDClhY^bl-qUo8xG^Tf)
z{C<4%OHj0Dpl`s5g<GDq;CY7+VD#B^|J~Fbta12yVa+y2$T)dC^jH@_3A?90Itda$
zUx=5tLQouL>~#5>R1*k^7vvWTHho7;Rl2w0J{6)dnvzm_!Cf@E{;T85@D?yTy>0EL
zc5xODEa%JePlCiefj3P;J0T-&^W$slAAm0T=RD)e+K@zWS{3%D1k&HJsc0BhW5%`(
zTMcK~A?b_OzGdgGfo$L8;3$eA80p(69CdLcq;!=;`<Fk%3}nnpYRCqy(c&7_L!sbD
z`C)awHU60NIQgORcrF?{`h3l)>zB{~+aT>N=EU@X;fEj2x<Pt#!d=?*J}`2>{{8(L
z9aI!>JY#WE6V=ZdUmI492HAX9HkCV+p@o37&q4dL(DcVyjUAg#qlKygH9m_iEc~g@
zwQ|)5OCL?%ZJQxM()&33E1DM|u`~4jPtRDeaCWlH(0UQIuFl*#`=S^vSU*u{FBn0C
zs-bgx+<Q>_mc-)*Il7=@@~p^2ei^1c*YLC*FF;dvl;yJ?OK4ux)8!z437AwbO1NaG
zh2Ht*uGqRBLQ+_GM$V;JP?e~;+hj5U4f;ClWiD}}fsfzbJiK!gHQaq>Z={d{I<9)x
zZ;Xfr72eWA!y6A_QtWD5%Jfc9|Ne;qpMfEYwe^0uLk(i$-t;5A6S*Lrp)_`qiNW}Z
zS0CNSeIV(Z%$9KeYmgzc<<RrfEdJQ(*}m{e1C;e9y?wGCgHGFxgwOV2kbXMF@ynz=
zBtQJT<(cGH&`3L0fYmT0t;$n3?$ra0XNtsniDh7<JEr5)ZZ$~WwN;<m>Bkzkhr(Lg
zQ(!7gJJLKx3sR%vxRwunMdKDeXV+%_024DuqS;mmV#WsUeD+O3kdYDbDz&Wu73_VZ
zW?aU^1dfgRu@%K=_;9IYr|K2Z`#|lRuqFY`MD0y5=8ngtwMS0&T9YtYy3%NIEkNa@
z&<M%xg_yp4gzTgF5=`psKXqaU4blZ~Ur#Lkg=$h>M#UBgqjz@2{BJK|@J{YrQL~c^
znw`(4J6367BCR6ilSV$OiSVwuET4esQ*1Kb{ufc-8G*v%Ck~-DK_&ZBy(XZ2w<ZI%
zi9k}_m$I1dKrrj@;alDKRj}~&`o>pBwa{STaIcaoKc;M47AkPf7cqbQj9EVX2$D7}
z*U{T^AM_0!6U~le=}PGkt7}venil`YA(`6?iFbch7%IF+i>D%X4(*kPwB<`B{68;(
zuAwWN4YHn~<Wwg+{`10U#PTrjm0KpD=VDWDZHz9M_Ab47;@%LXMZ}{!-pXJ$Kk(uw
zYjrSt*0lF-Y8IHd^2%$=%~J4A*0koOECF>k_^o|pCyf@#Ee0RwR-n!_;WX)HM@R_$
zDqYS_#gw}6h88zl&^-?KxsCWjM!=8LTfDb}smNuSOOra7EbCu&Qq%_1oG5;`%v&Hs
z;PQu!>ta!ts$Ipq7=K7%jB!U+AHfX$7zQjV1xvTY_CHj3h-uyEdEAeSnCL<}CYcce
zsViq=3dQcD(cxgH8@c<@sQ<QZqMZ(=W|&vY$v%d(%2Idjh>wtv#CP}InM%+q<^Le!
zk^p2_cej3DWds&2<@Pe{yU=_Q+z{v&hCWr~zvuF30rQ*Be!xw8FwOmHO^d=W&`?+^
z_H0-a4Klxk^_PvKCFkQ-N)JrIlJ#k)K#4t&*f+?b=jDwlEvYijVlHS(L4zoC`YA}{
zSu}QZ76#0(>;_zlN|?kh-zz~Cg;c}mTFV|OfTh+K{`8g#)Ko$3&j|`<@%D_&n*}>D
zg^j#mE~o-#4%g9+xz~ds0oAkn)T&S$*J01G++r|pw)~e$gcfKMd4|&nPavr;ZZKxd
z8&#}bq0F^^HyVRmnw9u&VS*#4kk$DHFy4PC-c+j;OeDvzzf*b@OzcvuPI<wB=7QgD
zx9#RY{bL)REwVdeiuRXlcQ<CBIgtu3w`HnmDXQEpN39Oh#xCycI?#?ruBR6|4f3MV
zmG9;4xQ;@arS?$#gYB3q*>`W!hlCoWT!&uI7NYT{j7N#Z`&s&#ICe~!1G8=Lm2!T|
z4@n+<qOxw?kioC8VygTjXnyP)ZF^%kYCTOKRb4a$ZQ8pci$x%){7Ln^>LZM*?z`(>
z{2Gj=YZjvHlkR}2gE5{Y!5B27ZYy{&ZVMzz4f?KJvlA?kd$-4js$hn5VvB;^J~Y7_
znjn2f1T7qPm#t9U1qq`2_tcbc1|!YaP);DL|Ec+btCzKb2`!!d9lzqyyiua(%8ccZ
zP1$#}zsCYBO6=Th>a-TMUC&rUFbzXJy{W65zeS@F`(sLm?Yv;Y=emeSWiqB(CAq6!
z@Q1|DE*pop{=hW0viF|}exNCBYNxSw3Z|ZVbnM#FD#*C8y6MrEouEfTH;%AV1oS02
zk4D!KFm?OMYl^pTqrL~9Ud_8#p|=%rT@h+7X!_u+GA&*PEVA8Hvc5TlW?rtk|5~dN
zGwz2@Ejm|0T7&vQtMwjWp*NsqUak+dMqq2DIY}^~somhte1zsDMxWP*+<_!vJ!*4D
zENHug2+zG3m`&<K)#aY`ES|0XhGTgfLQ7I!KjOo{Tu9NCk8Ubp)HmqLnx+`Qtf7ZH
zB$R?*?B!<O&Hk7upeCZVFb~?B;)U-q>M-F6=y^oQL{pyRV-WIy@tgFZ9H~*r5Hs1}
zR>{(_sw$iA+=+(t!{Q&yyZ58)>!w?!q&QK!_-ZBKhQS<=k@(ociWAz$%PIZ3C~4Ku
z@;-iyY7dugHzGKK)&rA=Be_k$J@-|UNrfyP$o1Ih*|0e#?HC%%AD9Q7>_1fX@HC6p
z?!N4@cod}`G4ZyxGsKjimixWpl2DVl`n~wcLQD~p*mh(z6;--B3ctPP2&tc*cRA>F
zqjK2+nXEfyV2N9!WvRdzleQmvlqaKxrbnWX+F2>kw|}lZevc%WKXQNRTZl0Tj{kZ-
zshfcMd*GFU8^V~l!un&P(K$%n`(7ce_aMN0w9%MDWhiAGy(4Ha3^Z!Tf5zI;kfGF&
zZJ|^S`l7Zo+)pmYY{b+}rk}Mz`@@Y66GyEf#lkH#sq6<@be9u(-aCt0-nQOJyd#OJ
zzuFXIgP(&DQ^BQ2t4dK{wUxqYwwsvFl^K(!CWgjTJ<3izJp^W>Wp;XfaKMc0{cRh6
z<b%23WSiPYS0P>ez>q}LFVJ*TS@?4X15y=}zeVqB1@q<k1G%N8sI~I4HixnunpQHc
z-G5;NW<=C%Fue!BOtO-JOPxGqgsAG9?Eeh9-}#Tb_ddalseKWaS9DQJnuI`<)DtwK
zYMLpnX@jPDxvF333W0?kpX^7JdC^jGB5#XS4QO-CwX#X^0x`ZZ-UcL$#vGOnTO6=O
z6~`ZZAxo!VMo3=|cxwX2GiLN610G?5obPpuhCVcPW2?88s0R2sQaTn`bqEd3iX1v<
zsfZ@h^rOB%V72e&NjdJ4J}^w#RCKda8%@8>(E3chf-2sr4D05z_{KHaUnk`>QMZcS
z_sx1{m~Ca@Ry;2O<}QDp=55W#H2D>qD&5mrc6u(lgyJC`>IPY8e}EKQ9kq+AN<gJ{
zchdF!!eG8#Z)>t$2$;K4^?KcmJQ&zicGURH0hF8U`>Jip8q)-7ud3$Mp|ZVIk-Msd
z!Jxff)waD=kde1{`&@%7m>>1>Y2W97uyd02i>uOTUT#a2N0tTZy~f)i6gYsc4%mmi
z;U}QJbDI@+pN__4Z@W__9qyo(`)={pekK~N3>}hxF^l@7bT#q^Sh`|G+SmN;W1wfP
z`%0;cj*#S;ll0->Ewt42f*4<#hlzbW;#X7dgMv37b%xwrz=Us&6?ak_s=sqh<CL5i
zYyNwOa2J(fQvA##gSkk+6kkq$aybf3)>3kER<6LbHPw|o-!`L(?v>X+9G3xKMD9)M
zgc*YQ(0$<~84?&Ds9csbo`nXvl*o?Dub>_G)>(MB5WslYwkfCMYAn5V`uyuN8zEWG
z)qUmucc6IH>6%+!hG5cfQ0(39AXeN{b9!^^EUG+y+pOlAH5wJT@cr!#F-)HE;fz<v
z0t+po=b!I=i3$D%$=Xc}FfN?1ZH>MUrn7%LNV7bEQbm_ZkdFvshUXE+%@qewv%~H;
zUger#`bogPn~h2=d{T;c4ba3CoZ6N4=?Yqk9nbysOaM&|+E7TJG{GS6O0HoyNzkVg
zWHAUi5%Uvg<K-M3FwrY+{%dF!QU-U}53IS4M)_aJGGql&vbB@O?bI;D<f83O427Ud
zbZg5}76UAabrili{u0gUHm;Ycw1RZol6_lmW`g15NRBJ?o0wWnFFl#ogE|xs9h2O3
z1=G)6RnFsO!}Q}AUE+2Dy~*Ayt|WZ|%p>wGANlTrigPO^U1s<pJyYS?F;{a)y;x?-
zKk@|=<+aHMuD4P9;}bq-^yI+2>O#PrFbC+_ldQV;l?fO~@7|y8y&O!WmDlqL+`|;g
zCw0~NE|xy?$T$C)2?=w7C$tPZQ5(B>$$ioQm{uFJovhq}X;V6CimPa--^FNP4=ETG
zyUs;^4CaJH@zo_QAQBBYul_Qs&6-bqBDvdq3L({c`eWX<9Lz8q_Sk)IC8W%&irutP
z1&cQ??bmoEg}P>CT<9#^POYv~y1cgq(n^R$n_d&aRB!4mm1`2^n-!dR@fkxR&jHKL
zS)0Mc14X{-rq_`EI_7n}r70M4FQez|>;O|QoVrHKl0mlRhe&rNY0w~)H|Uv@fM$jZ
z0u{c$1wEgw?mz3v#iWZg-UG!Nkb1^N%eCwn7}4NbH%a3{^<k@ZZ;S+E;@hq74EsKy
z83$jvJ;A(S!acS~#pV)V`s`2{9(xGM`<{30QOf~ysTo=xx3bY3QPiTmYzL?u@INI&
zpGKpnOmdZ8o1o-kJN|1r!e~B)e0Pkv8qF*Pm=w>OfR-S+Or82T&=z^)`B#-_goTa<
z0VE9>s)MI5KFx%T6B5QT;pL#c_x4?<e3tGxY}A<D_mhRc*^_Y@d1zjtivL~DQ!tgf
zdHE2z9+RnW1rA8{VVdBnL2|-<G%mYRCvxmDsMxqEdwz`zOAoBK%!MkDet}Eo`t=gf
z9b|n~`I8(P5A|#BP1z47&mA~-DK7}jt?X{>`Fa^s#Kq4#2*d!><zL0#R}ny-m(lLd
zmju8BR=Vm>MA6h1b@2})7tq{Zu8p{d0m_$bq(knrbiNz<rKQw^>7C!S<+F7mBQWH=
z<xD^L$oFVZ&<FV6=g?8@DVd+*oh@KtW634{<?0|^nR4UmoG{2K&heR6bpcGby>=^<
zSaw%QQ5hjEqTY~Y1$_cxptbeQ<p*wZpeE&dWWk*l&^$KOGnE<)mK5_wtPU?iiilUX
ziQjiHnJr@O7_bP58=fwl+BE?dQ*ZBC49o&$TtkIA_AOxK+OkI1Tyexa_--R-WD{bZ
z^tMzvppDu1g)SA1<$=)<g$lR*ta#$meLw7$9wy&6${n?cM@;9YM!mrpOeBmhSfmJ}
z__>#9v&9%KT1D|cnz;l?9}}P5-5!aCQ`f1*I0vISVddLN_RmpnmrAO22@_23I-v5g
zkc7t0yVriW!-XjxW@j>&*?^g2*E8Ra4WTI_e~O^h7^GM!U!w6HMw7J8t*ew?VX9nF
z_^y}(n4-Dm5Y4I=(wuIc$0t{UQm)Hi`$GI7;Y8rKx_$%D;e3J0$3FmRKc}iQcYOde
zt2W6j^kqTDwJirt?!N#X#vui<{Rtp@jpI>7l?FLtj*XLw{$Sx0+l@2VM=@>Ik)haW
z4;BVwzl2?HgY@War;bmwV^X%$?5l!nn0NwKpIfH@x_Un6;*p~ueNE)P_lYe2xznqy
z?7KhusW6=HJjaBDOPe*msab;OEABp)=pjHd%!!RD7DDxfR#E%>${=;{V7GAPK{Of@
zhTq4XzywLPL`ZE#v)lKc{n6urseE@IT+P1#>E|L((4Y80!r=nA#r_~Bf7*V&+aI9j
zX@gv`4r#P-o^xUFVm~UAzRv&kD-(72v1cE!?LZ~lWd#*EA!Mv;*1h(y64Li+K6rPo
z7Yr`Gef_k294y^mIGyk@2t;q%Cmbt90Li;5!n2RD;&tG`Ev}n2z~rGs=9fKckg-_S
zsUN)>GpHiE^z$R2HJSa;J5x^7z23Y<_eeA(bAD?JRWrcExZNkutiOV~qUT=kc`brz
zC%<y4IvS$M(9%@WiXoJ?b9i@pJH~`H%|1!&A2G=%-t%Q~1ZtR!6x=+*@@wX_xrzi6
zjX2kSIu>OENozZ~10NC5%;Cpmo8%R!!6sqyGc_1f?y<dS069R$zscVN=3htL-F+t;
z62zAY_MF;)CaWb5N9R03<32$#s#GD=|3oTpMIg)n4j=Hp00vNpJ(HRhIu5#2i&FF&
zBOv4ZdaH_cshC`#|0&w<EhI+<Sm#Lxpo!=>1*I4AU_|wx&{Kj9>W<rg&gXjuXx^SC
z(wN17BsCY=;^gaK8hkjh+p>m*dnU!6^L~)L<?S9Zg-w|FF8*EouT_xzQ(s&E@q1KR
z|79;>JtrpK+kds;$p9LO+u?d|!(GHY!Ra4gWDUxT4EK}+h}tCX9Z%L72lJOwq!O)Z
ztT?fFVfU_T&>d)4XK~LA^v1u<F{veBqU)+3K=lS{by1gpWW<f>&$X|*-LwI%Mn0Nc
zR5MJEjaA#J5QS-?>q}PD6F_EDx;wu`3mCJ1e^~kZT}+5rd$O`I6xBwrs<?O31$^I<
zND9ak#x$znT;=6@NZfUP@22NtnAEauIA`WEi&t+BvTSq%Q(LM_xQ6zC@vogl8?RTO
zUi-^$qb_8CieRg(*-|k~`Z0C-X@mgk{m#RDRxk}353+|-nhVjW$h6#tS!*;A|3!bd
z{e3h!Ggx&n`7v0!p1S4Jo=2G8^+x-qQZrcm{;A!{+y~NXuHb2xb!cw8=FEqyf?)ie
z@+0S72hetSn~tu~ZcHCaQ*+oqiN*{pp5F}0g=DUvIvtzCVC3R4yJcabpntJqb1_#2
zU_O}STbk+uWePmD<zXq1BDSw&x92fT^n75W>J<&@Iyj~$WfDN&=x*1IK_Za4g6Fl4
zdm<(XKU$X`{0w#f3N}~Y6NkzBH3;pwg{a+OW!}di2ohi2(_Q0&LC)R>S8@Yez~p(K
z*pc0j0CW7>x6(~&XzXFh)ZXRq(AccdQ29z(Ob9+!C>#}q8P$rPD}59&Bd1cF@^mw%
z#(r6=d3`xz7C6TI7*<1r%dIvfh;?Abhlm$vUMzzI=@XM;P2VA5QEMht`5<Ud+ZBHB
zrZ6hZy3%>hUK7o6h4MY?h+*m292>s5DNOd@5a;S_L(^k>wnsl+3&urN0v~kW1@l!Q
z@d5Q}fO)TLyFm&&Cf?QAdVk#!OghG(J+F|0)aLKGrFUj9+b*`lElR>@O4!E6q3||X
zERG$WToT8$ZryRQ%FUo{`E|sm?2akc>JE>-97Qv0pLi^&B!iI`eM8q3N->+@*MJAj
z5$KDl`{$M^0%$#zAY|3zj<At&kGq34nEKNAZe6|(n$ruMTB9-rDb#Fxc9-jzU?_0)
z_^ud8J-Q}qV5<To-n-saeOUq(*-1VWY4Qi-LCQWdffXRh@8o9d$S}aXAfCOvpbJyZ
z${l&NqY=~f+-;ZHKZDfxGZRlMt0CE(T`;@P7By^~|90DN8q5Zjr!E?tK@+R`4b#7W
z1hr<|&!J`^q#b%(!7rbKmUdMVWxQ1|v9+t!Gg21~(o7ud`+Pz7MCv1YOamq_eY!q#
zu?3Ay%>;*e_<)x3;&-%V%fZ59t0;$?vyjj#lzwS46O76PiKZqyqv>n0uawfVP}f4`
z&cyq{U=A@fDR<d1xp&}@?ZYNWuS|D;?N*LD9G`BuEwdc;w{iGJoj!^_FKrU2W)Q&O
z&7z~HkTw{q;yasfNdXo7oQJcLq|t=xLox@K0vL(9+>@jx1%7syd^)_Z1+;aiRb=V9
zqe&x7(s~0x*WH#+x5c^9{J!sBR}PRc-3J!y#wdcp3;E7rpY6b4+>g-0QGGN@8Yt>0
zx{8Kt-gkX(=m6tUm8<Kn+dzg}503$P0Stu5broAEqUn~v^*j28!GddmZQJr>Oc^_{
z!nlNhCX8c^_Gbe$!Cs{wk|c#HR`fL=nTiIbzZT>V*p`9s%eM%I3Atb*^{(ZKv!B4A
z@#);wC;`xLvz22c%Mne7Q(p9swt~48S9`7s$3tq<xc9+d`@vMvp@yUUZD8X1!0EIP
z9H5QQVp<1CLUJWfs~ck*WYnpgDB>DK^J{7RGU<b0eA|A9o1x-pPF8IHsp51r$gObO
zYvln<xb@Xk-oX^Kt<K;!JgtqXmKNVgLhO)KyHj9#vKMvol2Y=Ba;TMi_=N2`FE9~%
zcuQ<M4b_^y=yQF=!ZXh`l8ov)(8YbdeI7zIlyT0hl^hKzr=u@yOI(2o<62z%c|y@t
z6VHXCS#F@tvGnqDF)`3q93`21-2^S5r1enO5>08~;H%$KK<_!N`d2?3QOk|n)8$47
zFild~`Gfg;G_djXjm{NUz$AMr-V#xR*;LPbJVN;ii8{1virR0`#8HnoHL|Iga4+(1
zW9k=7%5B`f@mdG^ITk+Lt0;ykyE=lzG+BDa#!^P>f;9MM)$!}*CXBjnMhtK{$)Vl{
z=`Klkav@<AIz?i|(UybvZ$I8uh&m>B@khkVVNxqQ+pEGV(04pSSH>q4lJtaC&bLgV
z;j>A{&Rw_(>4w)_v<9p&>54~=blo;c%|5l3|3fvVhE`ggBAZ|aK`6T5`!X>4L8&P#
zrXIvMkk%)1K-B%vz09x55z|7dp6ury2idbVL@ic)FIT$!%0Jc`69o@PH4$0<F=+eg
zWs5&xuJpGoTs8?NK5<#bF6Dy}n+$AlM;^1S2z~fAjioaUJX#mDC($Hw;*9Zlhbj6#
zfy{zJNbD(8d3*aarnQ`%bxjY1gwUinZ`CEh_*mO*@6;l&G;>c+C4q!#86LWhU$r6a
zRQ(5j!W&3FF;%yJwG$d$oO4YX$w$-KBfj4QozRlP`?|E#mjQF8bfwhrRnXrdckuRV
zNzlEE+m7Uu0T$AqhzE+tV~VKLDBtskV7$yTSIb-n%-qPh|MiCknqxC@SaQ-pP3t`N
z+wM#NV@D<vLf)`LBHuHw^^fNuV|iZGiiTfkB8#@_WI_}g`|7#CpKSqVI@TzNv^zuU
zzSSIj38zt>!u^eE@|>VD`e;bpLodwe-Su_qDmUoijkf#T`w7i|^o*t4*o?{SY}fbV
zE0A(<Ec*n%CmQB@_A;?p4AUJY@65XSB41(cboK`X&@7wAcEa2dO-95>yx(&LlkeQk
zp*aPDQOlQXM=oe#Hjzv&)vxKG@1nDY%`r1bTHtC|d9{JX`@?31^4L-Ta96eXuVBz`
zG`YRoN)EIdPF-NW2*r%sEAM&sx}qhDPn<^fFc`eD?ERRDHX0e2$V_m#glV_$_AC_w
z%tnm*@@xBLG``|Vbk_SZ%!qq>>BpDrVDg5ouC)*sBv9uwQkPkQb_Tzs%0(|Ua>Qk2
z=f!a}yxi}g`r%?s887ZWX4(g7mjp?34z8HWK@Lpd4nxdGo96Sm3{lDFV~Ni<lz}e~
zKGnu-XaW5peOk9S%Av90x1Tt*!Z0<0EuLSViN>|*F2PFYF(Yg;_p(7Sq}^K*t%+O@
zrk_31w>WSHEY8%)EZY}?Mx)wzRyaNb&o3}fQ~g~qEn<JOV{sU2e{C7-W+j4{eLwp&
zOYTCd=&?Osg|%qXpJyal{|jU+tWU2Gok5e+8P5W<U!nS=8nmRY9!wOkRZnVuj%oaR
zy(%)yA#J_wNOFD-Sa=y@6eXh$x_6u~G>c-*N2@`hH{<+Z(f!S@l@nQLfzXk7*XuAU
zUv?@~Ka7AP6~s@v@?cC83X_WWFol%319#>8DpA{y8jWTScTAoSUA5yE6-~D7;LKTW
zfN2g5Ievt%ATL~7yBcCJ?6`H}Rm42#4I5{-v?o9+GA)q!8IHzZDf*;;bwlGrKhDYU
z_M*Y#Pe<?C%7Mw@10SOjve5MERqbW&qL6anT-=@aE6`|FTHS|97f2DxY2(%RK&AKg
zJK07Pz<^p3SI$K{OggjHUW*WfNwl|zo9G)b?Z~S9FO+Z;ZPdNr;sr(%aB84YHH*dH
z%ceYeJ29m$vukCVFRGp_a4r+*K{Fk<#7{rk3uzHX)%lw|K*!YCI_NWxme|Of@U=Ef
zQ?m6c1YMA_s$J8uA`z1%7pzuHHe<5qp|Y!1FIaxNEsd*U0L&XLN|+bafEn-cTUpv&
znCg?IbnAQ@rp7QgbwfWeWj-si;b#&g?0v<(u<sfeYWmWWxLyu0lWK%=->5*6({~yV
zIS(@=lA5z5GB81q%dhMnADS98?Jz(12#oT){kWP<4N`}A^?0)_G1)CmXYW)68nTJs
z`%<L=4Y7BIuAQDn-72;IUn-?B@xuG}L(5KM;`5oy=kI*Rl&<q0of=$dinC}*tT6!7
z_=d}CDfVEhId<pszz|G%Wjb@GQVCN@FAA=4Sc1h&iP9@)Il#zTp&%mXL(r4a{p78m
zKI-_oJ<GvW8ja<AJu0eJ3rXkOVjaVtVY*3~xcISjFrKx;WyH(@O$!#T)wHTZb2g$q
zKNK=h*NGgPfevj@_j$N!Lj?mRa-ZTD6bM73)fu4r78BH`zr9dFr~_5A(+&rIK`=eG
z({*$M3x}U~ORguVfVr^nh*#fA(Xg1x=b01Ska{PUUb(FuEoj}!5kIVfX6nDmTCdJR
zOR)B((d9%)q7T#_%~*%#ZAPx_ms3X*9Yy%S*L2MIz9lC<#R^OncvUt``JknPQk0aH
znP6%|V3%cC6qtG~{cP071`L^(_UevEVtPWx9fNRFOne|>7uQ^aNna*+sef|@3j*7<
zly1C43$L`e&USAGt(=yp9G#hHXwXQT;}t)e^AM69{LuzxqonPo#RW0xB9DQ2@i)jO
zYqIK;qyZZ0?_QV`K7=}S_THQO3_<v{AHp<ijM_Cv3Lb<TV$ywU018(?A}>Q|aK~X#
zeSP=wZYx$?zq`VWt<nrlek;3{5%vL++WV;Ld`pljtF1#v6R6WdkffrS0Lkw}!lyTM
zg4(y0-Ot0PLD}wI?O|IAF<DVHbe*~v${t^%y<Qfh=7D{Y+_QC<E>>DrBXAuu1l`1|
zxSyfgnKyT~o^gakP4149()+>eNhZzo>1{~o+4pvHN;B&BR1PG}%YmN#uLq@)o}dm=
zpSf3A1nSHQau;`y2Gi(Ktybx6Fu&WmW4Txx%6}Si@|PPE%yov9YFylbX+h;f<z3$}
z?c399A@5beOl8{%_itR7`oy<T$zlYI9Nql$g_AMrpQDuPyo|<F2LeyQNsLB0&6c^&
z`h&WR;O*T-J&>V%Cv>7}CuZ#ZsnM95glhJu<&^D+1|7<54IxQlXvm?Uu=fZ?oh@Z%
ztA$wit;+e?nQRTId$QVMGKbLuJX~_vvKvxU)t+_b9K^&$>ERXZBA}<}J$qDX7@AOe
zrd!FY1c}BYVT)r1U`WBry<lepW*l1-KVV@Bsf~pSlUA;n*fey8#-0S}?=PNMoAVWo
z#CG1-6Vn9?<#|y2XBwFI*NWGRQ9)gz%~|t>n8gcrRf5<`NVc2IxwfYh^{F4^m$SP9
z>5_KB>+<fPnZ%v(*AucaHF^1{V=YF^rj6=HcdP}Tz6K7l0R+&wx5%($TOp<#8@ipJ
zyctY|iaq;oLjw!voi6XI)&bLd{X-ivf>{2)R_DX&K~#m?qu4Xez=*z(N_hG#nmkTP
zkV_hc42=<P`*m(8b<LO8-rNu~*7TYluq0r@+otjTuORxsCYB-c1A|E^;}peebud&?
zy79kw{L-j^O~hw*^vkItCV7kjTD*H+QJOPQvXiVv1_`6y60TM8C1sG*sZ{ex#S-){
zHQiZ1MM2ZZ!~jk<qp9||I|R#qOjVtC$QCvL6WJxL@6N<yl3d)~VU=nyyXB(v9tsB}
zs)tLj6Mu>(V~P7F3s`)^D9Lne*-uOw5Q(WZlSA1so(()qKtGa*28Jopn6Ynnkk`oq
zNZ0PXu%EI5&8pCc6?K9jweNb|nL$~!#KV0se_jAHtgARSRR*B;@gaY@Q60$j{(7ux
z3q<o7Nu;MI)}Xm?#!zwO8PvCXxcVW`25nmm@rSFzXnwly{zU}<2?KCj^mbK9BPC>1
z(+be6#hMqcm+qm8x{#zGu0qTZ^|n)s#*py5V*AYuS5*3@No<>2875C?R-(mh@YCe6
z_xltA8rs54RXKkM(u|Fs@+N(Tl(rhBD<6)d`8hV%Pd*s3t>5-g>&|0ToT2=*T{8jn
zX0Ll^%T)!&reE+rEFht_!uT+&&2?zPHhJ|)e`nA>7f>`3(t;*-3(<2kb3v|r*B;G>
zE?{o^>&WKLb6~W;>_-3bb7+ykOOS4CfmD|F^X*!LCMOInCwW<T<5*~8CH@RDuAP{l
zLQdeD{C2C+>zZI%Lu~J}SqDhC<>Qwiu8bL`_D3t%ETX9rZk;C`x|mD}d30}UHYD;)
zKWr2xq9s{{ubl4eXxd1|?$z0wm~OG>b!dzkV%8Yx@$%k-Y=_tV>VGJKCgj*Q#}Uj}
z^Gkie`Pd3b%1hp}0TIzmlFHNAqB=}wChL%H)j=vW-f^t>EUFigTDzrD7)`GCV5Ejq
zG1*f0p6!nYR97PPg03e9mTp&9J%F##(xx|^ueFOn`FiUkjBH0t3f7NaEgpvE%r@&b
ztFLA8nuCem${Ap&@Dq0_J^<PHhhA+@o&vMnJ*SX07Z`8=b*h#fnBl$N;Vz{fbw9p3
z<P)xmDTzDXhn6m2%1<Z#c_Cd$dG=E||KJ)h7v2#!_Gt)>y*l%F<Gu<oB|bG__azxJ
zT4~p>mF9q-cdi^zyyJpKK6jDnIu8-k_-d%UHy5TK<(Qk0y9#EGT%u}}?*U1-W3HvG
zH3W@O*+cziWnhk5=9r;g4GUk!9ldQh(Bk%LzBaZ)m^kAxvL?41l6tEGy|3<rWJMAF
zy50at&%A3P?C%PR1w*or+D?MOdM>5gG<`6?>uGs&%w5!XuQ_6;XfYUV%qteU77sp;
zipt9I5l|(;`=GfRL{ko*d^u`IQCHLY3@yDDNLqPPQ6cUKYNcA{7bKg3+2mmP(A>+I
zlC3Dio%s>Wa%HQ1e#z=ji+T7l`y%wgS}?gU8KU78MYq+%1<>Fn$27dz5%pKB(PlR+
zLrbP7E}J21u+V)PvmJ0?>C-%$8>v4aX@YOCJdz5ScU`&4erADm`X#d283LI9LcW%&
zB?yV>{7KKwI6-Qs`U!z8cOcQLr(*x%gP8Q}!3A~mIY@mkrAB351pS70KJ(-kvvdd7
zfrnMfs4mIWWT~wOERAivmb0n|OeMZkda9EKrmDi@<^mVN4E0ve&Yh*0ZP$ALtrqEM
z!Z`hLiO_pA@8)Kn{;UpDJ#J^jDOQ2rlhSu3P9H!_%i_f4q79JtDkjb1RxPA6uix{h
z@j;T@6P1fQXV3`wo86h4lVEDMaG>lmCDf~OYT0U6Gf4jI$`$Fl7ZQa&bHCnx8q;)K
zDmT96fowajinXo&1?h_By+MNtXtCqrW0K4_H2zWX%6db6Fmccg+YCE^$!#LR9{g#L
z)_*v4F2@4Q5Aqsk-P(x8bG7@bm$<;|{!h{E9uO1w1ZEbWI6~^Bcd8W;7GO|bXJfB@
zC}>S`b=|Kv4T>#&I5Tgt{OUvZk5U5-ND*CI@Ok|WOm^~JDB0!?Is(C3xr!7pCUCp|
zE^P;x`Ngx%V&NR7s6}nHwda9UyQ8pAf(^34Dyr?3NJyfk81LuU1*zZiJjnYtgUQwd
z7xKz0FtsZ8;P7f9n7<{oh%(#2{EKxGF_W($?X+m)jhuckvBOB+djB~v$ziw6yI2Ym
zY!|-DDK>&i-sl`3J7KW&>C}eTTBR&MHcUGD<pd_I8Fc;0XNalYy_p+)r!oEY=7`H|
zZK&~Ks*3RDU^IGobXV_5ebBY!Mqb&}f(g;3LTgz4nSFQqToqdknl2P5j*H<19cIKc
z7b43r9p{VAym3Yoj-82ILKw}xSRN9yB@s+*t2$TNwgWWkbi0W@>A}>v&Sz&y?I@|S
zTi^l{gQ?GTu{R5|!NAjSDWy*tz=p;*DyL1rL{ECU*(x1O3>-OB9>Qw(R8ZsCq5&B7
zpEP(aT8df%;@bn1xY69v4_zIvO(9)tnZB7y3F>kc%uARvM6;uAwT(5Vm@awmpz1&$
zCVvRZ;yoe;rt}vzD0(5N$+=2aXn!#H`XliE)(8HO@XM3&Rc;kzn09E*33h|IXML1!
zOdT*{^5l5spbe;Nec0l>f*%u{4{JtL?EpVQ57n7>7=q4o1=-CcX*7HN@y7g_$Dnxr
z*!i0;HBqZtfkilV9J6iJ5VU=~4HK;m9M+gfgSjN>9>BCgbHb7#?R~E?Ip0d@vhOl9
zG(9G;w6q+|d<yuPBrysJ34&J-9Q=qG1~oaYoJTPEiiX?8GdyVMsmHdBCDLHVVp37P
zy#}(;N}c3#5B}Htp*KeuECUlox*v$zNvJK+TRcU>7ft6>XZrSQqw>|mYi@ctLb4UX
zF*KtQGh(w}ealWkUFrkM?MXdg&Y>~&XXqYCE>hwY_RfQpv&7{>r!QbSx13yn$}}2_
zly%wnlOHl_Ea!qc3ed#a@(0Z><EYzetV8w(H>429-ao$=2*&nh$3!gS!VG4A@I?<_
zuy}=A?SiHX=%6aex7Lq?GJ&mO$@RN2rR0U|Hpds3;+uY`dsPvp5hY4+-Wf>hU-?}4
zodOy;{$Xc?Pz_`>Tw=iXNXXzl*IXv93i>9slTR4SfxaHX8%FCDOcO|pziw3n>0cf(
zI90kKgYRsf(7R`-&1m~o^%Gub&guqL#py04?iZA|ABaPP9{dZdyO)E)!_O>&8(cs~
z_;Lki9;<&`;&vH4qnLV3n)70x5y}zQ^1Soj1ygRX@!#lc3C4~|sj0iFL(0vtIJIIj
zs>~5pTs_B+iH8i1W(OlM_|7wZQ9=<?<rI{+dR9U@WA}j_dt@*rX{nya=?$hV&5h2S
zWk>B=!Xk@>3@{wfLX#UW20c8^ufX~PXv)6W*R%cxq+M?-Bs;E$6l8IBjCv4~1q{k}
zoCRS1V#rlN@7HMR)@war$Acyr_HJ?9f{-l!RY`X6IQZ@xoB8OeHJC`2yYFk1g$DeD
zS{k;xW5#)@Gm5c$LHD+H-J>@O!KnD?$~P+>p``@+-7T-8(SYjbtiDz|(8+)0gVHHZ
zOy*NwpuG%7GlwZR1l|XumUC}%a;?O`N82ptBSjcZa@H$1Q`*qjW`(#&>lI+aFus1r
zbQ)@nF#cvOPz-9)?>ses&p?Ag?d5Vs5d4@r9dt35#k+VGxKH=lLZYjL+#35IXoy37
zVJF!Uj6X(ir@2ePlI_r0w{@35J#FI9V@?;a5Ri5x%2yrDl9nfW(RmP)v3&T8^lnyM
z;a=Q${1-?ZEcT0wB!KaY&0^?fAEbWh$<6&{j>+t;%Z@!`f&t^k!uzJaAanZ8tfovD
zTHpzkI+w2x<{Sfgy^gI#5k*M~??VYJ|Bk<Qh%AqWT)*Pj@?21Pso{Iiqi9I)ziEEc
z=OU(DUd<^%e+U^wwbM7UJXreh=>g(pJxIU3a4v<7h$bXo%e>tu3EK4Z9*OOc#MGZP
zkHW?gCbb_bW^9i~lhJ<0-z|nQMKVu3ASW2pV?Pz?&h(-Vi*4rVX1g&xWb}C*TN&!(
zkaU8_Gcd7U?namUS~UM+v^D)wGNg4}=@RqVj0Urxe&BG^MN6@3zK?ygMLlC%GA<Qq
zqF;Ig4q``{s6VM&@vTV)m~|@9yD%Mv7KXAK{SIA2HDzA5I|UO^-#1ym$;x0zQ7ygM
z9JK~5%}0;NhQxx7vG<-gj5LwIhvakHeFV_x?Z0v+#08V0FV)1|zXj$+<;M4?_hG7p
z|IW<bR7_uSQtILPbCAY4V_iG64}`KO^sN|T>8Cs^J@yS(KwIiF|8_opkmxL8y5d?G
z=y@qttn;%FOmlfBY(CxxQmiNV>||NIXnu7?h~FKMu&`fy>@>!NLk^9mj{KmyY3qxx
z2dyzN{n}x!Em@!^?0(y>Rfb^ch<nVwIb%%f*Lx`t7KF+EC$)nGzo54C8pWn$ezXuC
zYL!O6&59Rq9?ou2#Ppk4XFB9rdccSgd_wLDnD0J);7U#uBrZQ9cW!C|Q{2BR>Ru7V
z#1}UXIDEa1mfqhlJgMsh89oQXPCv*2Gx`&%c~@${_({>~8?C&MknY6Xde#UHiLAQH
z(Gdqmbd)rIIK@L+>BW4--gjU|yXbnYz#^JCEqQX)cSW#Rt8;(7N(`#>HOZwD%h05p
ze~YT-8chE+T<yM8%;NuNw=Zp91k4|rA6?%hK{9RMp4SZ*(NDjg`RFbt%AvXJ*m#(U
z+TTwo%Y_Pn#T$Y-79WZ*L+PX4Xm&lA5zd$LloJPYuLTNTNpnGhjbsepnw^l|D;PBz
z_yE(47Mb@OMN!w4H9c!}F!;W^Q}SX`8Ax|@TK9w=1{(Uf`$XC~A@OYEPLbOIm`w?u
zv@q!axl7EFP`3myX_2`3MSUq)IMkz;ZBBy}yv|#svK7)-=M*bmuEm7SO0(ad2SM`8
z$7Sh@`!U_=Ye}q#5F}ZiI&2p*0cH=lluAu!gJ0e+W$MB(rsg_bss307Ci&$JrxMjs
zS-*2b`F2h)`XZLseWL^DT3LB|Ytl_HafAC*hp`eEkz-a?-d6#$772XYUpt^VH{_xh
zY=ar^fAl)0NMPy{N>^}Z9_SEiAJ>!eM{j&2?hB~9VDb~feQmKEOsQv9KMt(|^U4Q4
z70viSYSDDg*cxxZ45AJl@`*xKK1Zu{XohHr&-rxoxHjrvZK%J`gQe>ZX1!|ww^8ft
z@z@wP5!BHt9{ogk8d5)U%99l~Af}hg;bO~aFmd-(@Qbh*NV@b&l2oz+E#6|hPqoTr
z>3aW*zS*T{IOFS_(^`K_YbiLdchCY&FAu+Qv~LaQEbDo*{jnihOp__#w=YFAf!i-G
z(hM+}>zID#GF~t=p7CaHjXvlYxv-ooSQ<><TqhekV1oveIYxdpbf8%;&N8(rgmNA!
zCv}PXqj9CN-8ZV+QA?Rq-K|JnOjKdAsjb-#%7(V{C3ab3Qk;O%W8yn7yXsiLAz5S4
zyEK**La{?ki4STv8wi+DV-fyngA`!uzIJ~~Ndk)o`OWw2vmh<G!@~DdHD;VTF8YXQ
zhDkcHU!LFjfElG5Xkva_Q14y~j@NFvpvCFrLF3c|n2nQ>k<Au}i8nZGHk!?Wq5Eoa
zR`H)8HSd+ztu9+kvU7f2_D&g1iAGEoOj~10;Ss(Gt28hpXRR^kzY)_JpHlCQu*S)E
ztNzYKKTM}uC>QaGW9nP3d5Qj`XkOfiZu0UHm@tbU6>HbQ<g2eLnQxuZkdklUBbyY&
z400RNICBrv?w8m~R6Rt)@JHT#$plOan`W>}E{BY%sII3wpP~X^>4O|IOiZlUCfD*Y
z5Ry_98eJXFf(gI+pKT{KFs0^--Q5aSoaXD%*OzSoBg1C-4)oh-v~yZm?Snb6Q#pUH
zT8w}aSEetOkqlXS=bcV|JPRl8ta!l@avPI8b@1%!`(V*`kVA5p6dKq+sV8~a1dT}$
z-*#z5W0KUeZvlMDV9DQZ44=A;2}hFq`Ht;C&BUj_FwGJ5sbA?jyP^eCEDoo%=4he$
zgA#XB_C7;1Yo+HTyeVLz<0SIG*96MP93xM}vG6=I$b;wMMKr!H^xK$l0%*+0aIs#)
zM2mTpy?ua=QW(vuhb<whJN@X+q>3;aY&1D}^1^P&u(Xmo@~asVh0i&ADhz^>J;UqT
zvlAc<En6wOa~Mr+a2t7tmq2H(<KgvZ%+Pq$G6gKk(j5ctd*Z$dqKWuS-ZV}jFu%2f
z@8sP@&}gb&AD3DNW+KbQZ)lwZ(_7{p+Pj@u@d8^5co(7RchA3nJUoj&J0}l%FNC4V
z{cAp|ik(3%9ft$ek8z^WxA*o)pO6PkOW&R_+ga4@MmtoV?hht&ck!<6iNUmkHkac>
z2w>Rn)ZnUPagb*0axwJ&2hdO4ChL~>1e2dh0l`hGAT?sGN=a%MB)TnhA~_+@R(;yE
zNSq(dHAiXma|fUfF8Q=f`vg?&BC@@vR~Sr<t@Mbl?gP^@NatN=6zF_bvBH5642n#2
z4+$o)=1;_nUAIctK?YAz=@UhO$sHOT>Z=_f{Z^k*^C>4ZzF~_;En_8!zM^&Xx+p{=
zZ+Ko)swt?I?IztwJ{L{wGd;RDSq01wRU~N|8l&!Szw#=_Q&1J{zAB-}1@sB8+Ijs7
zi+?tYZz>CML=(rAYbBpvM5CW&%BAxqAR%*eZ~6LiP@AigpLkvq%vY^_{(f;16IA#Q
zu6n=;DyOzA^p2FFkrNx9r_4v7;k7zE+jXrW?VE4*-d=M|&(9tDU=a?g#N=h!&X!>+
zhs@*Z@s%iB%;@p6Fn=&%uC_1Zz6R)dke2DgW6k1&`P;nxt-<&g7rEtO=O8_<S=vSR
z9U7Aj533ZZ!IWRsit*F?FfHTxGeJrZByo4G4C&Ygrfj>bMSb_6zG<L!XUG=~xL2+X
zDT)Te=ZhNF$oWEI?_MBa&x0!UOY(T%C!ijyqlVXZC7|h_(h1RPYEWKK>Whn(nyB}c
zP4v2maj>W~$h+ET7Ssoncbs}1hKUoKw{z^fj~Y{tj>(=n1QxW5lHvnMVEUDD=H8WA
zkfHZ@T~6(OG}-gE>3iWbl<y}z{x)11{jl{d3z8s!iuIAn(j)$u%st9QIV6cmS2=7+
zW@>0W|JYQD84--Kg$sVKdI>49o4ogT8KIH-pXUuSTv+;|b#yGL3^EMA?Q44Y3yg+l
z8nUa^VEPXl{q*HspxjBxciLYU%#WXvU=!^FKXT7T?--MYr1bJx=S={UetoFzzT5)o
zNo<3sYgN(ok+kq-l`W{ho?T6a&kiizdo?Pmeif4P7TFm#0g$@d(PUG_3rHf}PFrGL
zKv-uDr9#dHG2^1N(?5zr!o*rWrOHf9qq};Z5gNygpeIuS&UTQ-v`Rgu^c7V0`x;mZ
zCa`#`G;`zhLre&8ZhGZI#Iz3=ZbwoUF}XQEmNd$Wv(k@lZeG#>=_eodMm=N20rM2^
zt?!L7?Hq+Y;>lMq-QqxyKcoR^*59)`FKoe#+pky4DxSlXcutBx#Se`Fy@aRbm(Wt6
z{lyWR)nF-v@ap3&A4qx|xH&A~7wR#6Q|R<X7>#F0@b!wnf~4=UoJ;8im@r1}_EvR3
zT~EAz*xR~ediYa^OHZsYCGLEHWV|XED;bH`SalD47nNdWogg5l#fqNqDpII5XSn~_
zf;D6a9G_lqI)VDv=!p~dm9gSRyK-e(5m<~;K6l;M2GRshhL&%&!n8fB_Hsl#LbKZP
zOR5HI!1$<_mG<#8)UBd+#HaKo80$Z<@NGv7q}~50BeuI6#q*gd>5jXASqHlZn*{{X
zf|oqm%ZDG5SFdg`_<kG{>;<>WKFfe~4&C#EHGaS(sZvYz00B~jn_jd^<e<s;*FL<d
z{g5ajJ8LI>77}MAN`hU>QD+-{!k$@&8dv`G-D#o=$}YBl>KlictarR5Hhl*eBg-ms
zdfi5|6Nf(gBxRsQPjSPP4mL0>>FnZBrVVBz3=IQvv%nOdr?iB#L&lM#le%ZVLYB73
zTVUS|m<`?@CY#hCrPO(4@A+6vS+UFA$iW>=Dhx)RH1|OZT1C6mi~GRhH)n0xXa|%r
z9Y0*jD~u_<K7&fU0ch;<34@AV2+g`~l(c!pK%Ga}x47=?fpkSVaV=w4uoPJNd<9Me
zQ)KCIKlV~g+{U$k^{@$~kAF2DQaO$)X!mL#=<tJ?+WYXN?<}N9+*r|Y{SBHr;l1JG
z#wo<yfA{Oq$6PQ+wAv^!<&0@Moem`f2SB~`LaVes3y<2P9HOReA%&fnT+Fc#j7(h)
zUsj!jmR41d?=qRfM6PW9n_Ub@_LFm4?PZSM^!5GZ^mjqM9^;$tNQHr)?>+{31`$vL
z=+Mu14@OO{m25X>Qow@Gf%l=~4<Si+CQPs6DOfs_I&s!r3CvxKU;m-P8PX5q#~Et}
zQHw~vS>S9P8oewfdy>V6=!ZEt)f%hN$d1RW4hLNZ^WV;{GhC$uDc#NQEI1CLr99Q5
z;*JndS2EUjcV{$9&(6MKuLY<+)!DG@l@}zt_Prdqvko)3622<ki3c5PpN9-BP#}?;
zL!^xhf$4x|4*itdpd+>NvrG6qYTva{*`2NnrZ#m3Jach?R6(k}%KgKrKb2JRY}sx!
z_V(EKovj&YY1GrUIno)`yt9$j>vRNJ8|L3pe?Y+GBrLrYsX?O~VtErYv%yr~zNbfb
z8H3iZIVUpvsxV_KaDOl;j|qF<r#AVMF(cwgB`(33>=?v5Z{G%%PT3CCy>NgG`N=P?
zeCr`a!(^!AoB^79Aq6!xpQ44=x~?Rp*J#{2;^1w2B`~hOzHFqu3v{afBIvryp=rA0
zkJp<mK}(L7QiW_8W^`pNzp1to)07<MlntJPr2q|;&D2MjYUl9IEPEB|0iDrKufxE%
zJ#!<YE}EE9yf;g1Bn9=?NFDq9(I0g_?&^Lf<_ik89uo14CZNHm%g>U&mtkVc4}<py
zE<s{)7`&Cg0gRuCe004|6pfxZdT!4L4M^Cr{ak%rBU(K2RD0}A3#3>^Y*=~j0hsyT
z&VT3aEJ{%Oxz%4J3=%wQ_#2g7QE$xk?Jrz>F<G^GZtbEeYP=jh_Cl}@k_9<_%I|7N
z%uY{vAt3?K=6$JnStdUwAH4sA>p=vX)iK)8aETuhUfI;?Q+rXI?adwB{w&>mOngCV
zfQ9pQukR4H&7sMe$obWioS=Jr{#*KHYe;&q{kqQec1S(7Y_qZ9X4Jf}A#p-E0kw!u
z-5HuR0g?W$Dy}00G`hN?^G-++q&;_T<cxg|X7ei!8h<K6Q#W5QY_}Z2#IB7~1`ll@
zeZ(pH{kP|!$4hR#i;W8^-&gq|G3OAbx<?e7b8f;^)!qrN5`gl$KVW&p1U(-F62sh5
zKr4GqP13n4P<S-|)Au|VOyTvXP5B!_id{_`kNrAKEwp(T%H|1{wvT)3aVvq@B>N{(
zhBQn%wbNB{=U3GECU@NbLkk-DX)nt@5QEwu-zG>MG6ZuWmQoSf@?iGy9{F!+S(w(7
z#1rym9+KWU*&Zl3jYgIv%y7azOn(2iF8D|ti>D(>gy(%soP8MOZ6ykR8oeYF)KXAE
z>zf_k>YS+T*wc+aHuiw2<qQoYnItr*rr3F?K@ttIZEG~K*^Q<c)YsHam%u=ttNhJ%
z5on=@f9XzWKI$JymEILVfmA4~skUb;n)7f|9d=!b$!4!F+va9thS@Z(G&%#?RF6^O
z0$BA#-k$VPIshh>w`WIaX<_o_!P=j;1k5HK6XCx28M3X@FgB~LMMEpMtlBsv2U?<=
zFK(Wz0z(SDuR8n(P|K0sxALMz05h>X_nP1%NZ+5;Z*G(a=}+apDB0}-6CFyMV&as*
z#Itp}^6EZl<ol`M6vao7a5`Ckyl@pH*oD*{`LPcz%_|x!v`vA*ny={xUQ0r<lyvW*
z>=mHz&QM;^K{J*et3q48m_r8L7HUneN23)hpdQyC>V83$=kGB@RW~}$J`f7SRO#48
zrFuWm=eosGuAG3{EBwy2yp~4ER=#Nu-as_Uc4)2Q4>Q!S9)EfALolXDm5Ez7&Y=06
z<whYZ+dwQq&F5GJ0d?n@YW0b=pr3G8xr0<0XuQU$)7&nNns1ufe%`MOCX9~D-{xS=
zM<;t^uJ9PlrR~PU+Ov>M<ndiOLIoA1`)OS@7|fhoq4fQ%IHsKu@-z~4frRt+VQc7_
zXlir#i^rxjm=bIdzS_qIGz%S>rfK$onG1`OA^K5hA&MZ6l~{U2arRtx_IJ>AJK<yY
z6P7-D`+2}(N(@vt+aG&-!XM38*a(U?&Z3T32V2gc+L)BEPTxG<0Ml?ZMezI`)L?eU
ze)@tnq@LfrtNBqQ>TLA*syE;Q=AShET1wFXGdk}oxA%#n+1p$1j#^fMv2arn>b-oH
z4&R$dB8Z^HaGT@xE3H7euj}`ad5kH(K@Pq;E?~h+`^_$TGNjM;+u1xmjH%y>+l1G>
z1(U;Bo1Wim11;?LFYah&&99NYC#TG3L3wFNNvu{1n12}_c$JWi8iPijuTf^<%9;3<
zbJu(^@!**Y3e*g6L&|>aG$#Sf-aY;7UiBhmbZwK|QPGGQxgWVLubcx@ugXWxoKyiV
z-Obk<Uq`dzeow(p^9PU~+W9l4>;{<UEQ=8Bd<o)oM^=1mAwZh>Jz>uy&5)5Lca#);
z0uwiwoN5m5g9M%Xc4LXQXeMwAp-N){(&fm}OU!mma=Aw{{j>*+RQ8+yvduyrnXjsT
zLKZ(<Y;6`&)C4V~NoHo>S@Zq!E)kjdSuna@eAk!l8_`%g?McT68`R#LvZj=%3;NdS
zhLdsyFdJ8ZPhFHSSPV|0?--~-eTTUAJ0BYWKhLc{Lplw?M7r1X%?)a3@}bA;Da|J!
z(;}a6N`-)F#>HpqMRGA|pWN*)7IL6w1s0TTW6cAxQ_Es<!_h!&mxhFj7bHnvea^u&
z0KMd<YV_S7jA!55?!F+8(Ds&NA}3hmHJLf}_$3b{9Zsyj-h;r<{%sn(kxwx7!87@W
zIa@S)*IL6_`X?mHx%!z$DnJ5(bS}g52#c==NmK4`0Hw|B+MQCY_H=YrJzsMbjHy3j
zAAHX8w=F8#e%Dh_oDyGcS09A*w!Gyg9|u9+{$bhgHiqcaQx3^idkCo`cZB6n%YwEG
zr$3a37GiSLp!w;G+K_H5QsVi7ioTyH?|HTsqq5@RMnAg*Om0daa~(E91={^PJOIn?
zdRL|^&#gnX3DnxNtzKYh%M<a1KpLh5%|#qkA)-cWKeE}gGBip~($}wA4kkXFb$qw7
z0ur`Q_U^nK3kF3#zrUz64#|XT6BSxHXngBSPc8*3NIu7Nq~q2uNUb=WFZDzLQZnaU
zZ-=i2KRj{>tEQQ#>*J;yk$SplBHQXSP(KT4l}#IL*U&*voUdH*$9YT+biR_a*#{Lb
zS#J(9XaR#o%lzcMC}8ZX^@!u>5Li@y%U86!iKT-$hVO8mz=Q_Z6}H(XXngxjH$|7F
zEB0tr3;D5l`qeCbF>J`<OQPEtKi7eFiI}5ne`a9n%lcOlt4%@08tIU>Su;r3lwIbD
zwqx2O#hCst$*4VK)8oN}7Bt|$<we69O-$NwW-7o_1B_?Nx<;8afgW#nu54l%B-Plq
zTe}@Z<MwBFUKlU}{ogZt6dj+SrB9=)AEM=G@rI5@i0EC&xYqwT)9N&sz6D*<Yx2Ql
z+ZN8VnO%_WKsL`ASb`P=dv3kiCkTmYRjaYCJy`nCy(djy3e+*0${W|FfSJc{U!UDM
zi^h7l`$-y9p^1<Ap~UPw%us2$5Pf7b8l$7K)*%DXnKa8|9Jm|P$S+3TzEp%Xd-lh&
z!h12pef(rekp~!ivuJ0^qk<OCVdoqDM<H!GI+N}FckoqU_r<jthNuHMhL^95Mk9FA
z%6!2ZlOONt;&<K!30-YPZ9;iqNpnHkM}`K1bz{Z3LI{w7Qr7hyl0oxsohd2u%P~El
zfA}75!9=$=dudXSFkwXbxY%hOu)x;25XOHO(^lyh)t!n5{nK-|M}h><<Bu;J%AgDA
z(OZ7oW?(lY^B6qtbv=!K*{kTU77PPR&AOe@3H4w;HVA(T<3&Bo6TdF4uY&YfH`D_z
zWr9L$fhgLk1W;@9;PifBF-RMJEfCynfC*if0@#PUKqdR+4Xc-$XeN#85ZKJpO`k1W
zd2gLU^%g{`xv(av{@H&-v<0J<gQJZa>nLbR&)ILI-58`!(W-T&Ju%ynIL(OOH<&(h
zAb?CaW%0DKKnN?q=<^32_qWTV+1KyFpR)~PvP$luw-Xn^lG|i%3LpVw-8Yyp2vPR{
z<BEs$Q8fKNw<T8}VPa_$g<@d`;`IW5-ZdhC3CD;Xai5LBw7|>b&bzx%XALLx4=DyS
zf~sQ8+Ls}HX6qKoBYa@=t@elJA_qvjur^U)c0XjJ{7QW<xC+uwB)g{9Phon4*sHbL
z(<t(wlaba*7cka(YvJ;%yO0pDWpl@R159?G+oH`M&5DnK<mIJCm^dsiHL<x56Q$fE
zi=%I%vF}D{N{in>^QCXom5?9OwT9$5B$q*YhM{bbRX=2uE=zsZ^$T>JqrA)d5`jA1
zzug%^9bh)K%75u+G8#4G|Du|C6f8Z?`PC}ahB}nD6?_Sr1|@nf&`6?;#dn`C)eN)h
zC6#@rFV_R5RW>;_HB3z5e_auNe?J;YY|R|0DaH)>_P{}<vzU!|t1Mt;1!`F#DQr>0
z013hMAJ~@+Q5bL%OL^!5iGu~_JOkX&pwQ**Qdgf~%6cWklu<g1FRm-PO)&!<AC381
zWLSJTvVOy+t>0L7;I&g#{DNk-aHz;dUci)k+g1Apg~3Q-dUiTw>EuFtlcZu5F!Nx~
zO6do$G3iSriN{C^&5SGVlzkElx?wO<=VPL_^93*XkFxx9?y3zh_ZCRaGg#yQVizQ<
zORpogpF-6~s&{YiXW_;5pq&Cc^ue6g<VDkSkI-P;yTj$)325Mxh^y+Q7SMHr>uD#`
z1=Ow?lHBFth^f8T`j6DLL&ozoz3gw_F%kCKAFVJ$Q>aS%U~V92gWQbMPa@Egfa6hp
zt_V!KU}m|@Bo!^jyyaxy=nO_3Y_1M2bHMcJ;$pe_6=<+0E%@?|Pmoj{?z5ox15($m
z9G+6EM5Q(X?RoF4!AQk5>ILc~rg%Bty*)RCemphbwNspl326_wxdSdknq*2duWBQv
z)~Z}rrUqlWrz>ACwGRAvVz@@^k}&FB`Mz_Kq=^RJoTRmG#-K*;3Ky|40nDl?S3Jp^
zM05SoR~QbbK<{NCqXLa5s50SVH03}7Bq!(IuUr_0Y{bxO;c5HP@JnOkmW4P-bKCjm
zq}g877wFlrF(Vi(GPojdaIHnHa@S;bT(<*tmxo?mu`z`7hqc>jhn+C>onG{h<z<*~
z**M_<MH2P0wa)5r4S=4L%MWUui35`v4vD+fZli@>udro+gql4AEQWU`K>D)-DSMNi
zqtYOX5$!Mm6@PxOTKkk$?~?6O@RLoLp{V~Vn)f`KSG#I-<M;xkjrWAwFPulS-4~gW
zN6kT--sIjaEn!R<NOvJb1%uge`+oI7c}P_rw8(3cKqH%y@AFYAFsUGL=N@inOt)3N
zb5oy!dKLGDtoGwZwS-?&U|S35ewJHRTJQ-BdwbpHJbeHxy*j?ka9tW?u)S;AkuZTt
zO>L<^m_cB!ium+GW)h}r|7g5$;TP(u>pGBn#STnt6__zjW|hw``}FeYN6^38=pZ}4
z3Ftj-Oz!y}0eY?y%EX`7p+z_M9+5{Qn9A{m_$cEvDjD%7Yo8HD%$Lp!b{`3#&|I7F
zNRXvlCp6dH@}v_81kV3DhyFkRRXWFUhl@5g7agr_F522&H2Lp;pl)Sv_y2wQ56k}V
z3;Z9xJ$2RO;{VGB|IIG{nWq0!_Wj@1`wyS?zkkHat5$Zxw&tc+%`cf>FuMBxzH|Sm
z!T)N-f7Rtb)t-Of7k@+d|0)0ZM`i!A_WaEz=)dgGzth;icdGr1L-2oi+`o_H|5We(
zefRjSBlGW%;=hOz|GpjmIRN|LH{#zPxj&3+e_25JXAJa@J>_>|;J;Yezwgd}0*C)>
zBmQu}{#T9iPuc2UwBSGY=-<ov{yEtBol4_hohtu+(E9t*<3BqE|1hlo*--xvAMu|X
z<$slq{Z(x1KjzJUf^Ga(rT+&C{D#i`<0}4rocYJk_*;O-Uj}@ChxzpnoAa-_#~&HG
z-x=zDFID<it^F4T@83YieuI5}gHQZ*(Dqjetlyh~{mVw|-}m#sMQZM^GVAY!v)>zC
zf3KVWUUv3(DDWGSs=qbS@Ry}6{zOLMw+8UvVD8_`|9`K{`yEM{e>}YZ{?gE2l`#3M
z*v#LoIlpx>{)@cuH$;>_%SiubzVTPJ?5{HFKZPN`qwxQqtbzQFdGvQ!i22Lz`#UVb
z{stKSi<Yf_!#<yX%_M)?4gRv#>fh>szn5wK;cWlEzy08EwiNVF4a~pFT7LtB|Eo6k
z{Z>@^r~Kfr_m2NyLii6O=07!{e=m~zH{PW5TW8lF@Tq^b?C(^re|Sy$-{1K5pCb)_
zU~2qXNBa9-{!?iCmu(~a4Waap>b&0y41Z+({|`9sznVh(>j&mf%dOu)g8p8c4SuhA
z{xtjivkB<mun7CN$VvU@R(pOci}@42{U0Ove>0f)t=#%Q7U%yk>i=b2_Fpvc-yzrf
zSFM5n-Ht*2S#tCz){6cPmw$Hg>kpXk52QH%fPMdCe)ji0=JW@e{(r}&?0;Vn{moWd
z{uT#t|KV+adX4P2t^xg3jN)%GWB%-(i+{>n|L{uGpNOx2vwi%(VJF7Fzs>OXPL4mj
zso;-inEs;h^*4*Lf0@JoVbK1wn=Ahg3q`-dlz;j-(4ROD^9Q*4U$r>)&yCV=9@F1^
z-*4Do^B3*p`^~WOR{`1|&g}f%-ly;{D~P{-b>a7>GXLyK`(Fl>|FM?<{Ks<Qf0d>F
z#P-_189n)}%<7LYl0Qofe(#9?EmnblYh{1-HAR01^8W`+{hMEM^6%I7hfmZ0UW)Q3
zUKjK)cJwy`=I@2nzXRp@zamNShqJE#qA>amNwdG+P5!X0|6|$k-)hP~ejfH;?ezao
zY>oI28t|`*w|{n{{=W!Bza@q9M^8}v;ibLbd*AP^;%|2OpNmd^UE%*LR(5}{3iz|I
z^`FWX{!?D?f5%+?vrOPmoFVwV=WPGJd)fc89{p26|33yJzt!%4<UL&fOf&zoG|!*V
z5q}jv|NBwqPrShD_Xf+q+18c6%szki%>}=qy?-=|_SeyoKl|p8zwBtv-(V;Fznt9u
zKewGgA+aur|MM&Ag)2t3wuJxjFY^fG|M?}&&iKlClS}`77wgmiyks{qvNMvpWMgD`
z@qej@&Bp4w)K%8^|5eut)(^>r^=oNqTKa!|;C~nT|9l|P@xL#e{`;X9tgf*BEdDQd
zvkJ00_g@!~&IJQMhX49s3po9E0WjimRk3$W<$o!#TuRI6s)?1Ql<h^9UH(T~Nj3z6
z*c;Yg%d)1!TO1B2zbV@BXz!+!JE($dhj~vnm72pF!d&aZ_NYMFa9sT^Q*+(rG;-3q
zpJ&_3+bSg4OeuHcFXc3?VV4IZE6ys0>+!sE*|QPa?zK%3F5cBzu6oJ+ETy%zW-zJw
zCVg=3qq_V1Uwz^WuP$7@!qMp|r>85q>}s0X$L*(U1n&ABdzq7EUidzi$jlycXWo$e
zm7;N4?V3sIoxE$^wCBBNd2Ou)Z`>trX7>s9-d7Rbs3KnZ!6oJXvl5G})h|K{KhHn4
zxZ?OxZ)xFq<sr$PcdNw2l5-pIr>-o;ljAV1>$+ierCu_7to`Z8qyonm6-s9=FS6;s
zknXy;b-|)%#;y3==aGsp^V-U1nWx@gEnfqO?zU@p_HDQo_+&r-vxn(c?+Gv8=RMGh
z<oJ16VA*WOkTmhjUc;NEL>HG;vWMOIHtmhx8$+$E(l|Of-C7ZCG?%o;qQmIhy$Y(`
z;hXo~UNtP)@I|yrHty~^wvkYQy<boM+TWAx_K7X{-u)7za`v$K_-}cm{>^?Z9G-C(
z93C1ANbrxBotp9y)K9#=>a}ga8QaSO9W5*KhrAPtq}v|IUe)94pAP~vRMEQ*8cEfK
zLe~bC^AsJ=^<cg~elht3U0sHroH1EV-$vwlbmaQZYzLJnz1?Z4#<r1a_5m#i)m{N+
z3z_sz_Slg4numK+b0Q8&GcQpiKJJ=XvFSka%A1OoU0cqqO%cD4zsqn<HfYdJ;VceH
z6kqj(JvUmtg5m5?bKxmnX}0nU=lI;i;P&7Y@IH5L<<6_}o1*2;@n%22zT*)2z-naD
z(aP;><+ZzDkErLDT}DDe!Y5K%!J$%xyQ^MFFVfS=@4p||C%W4C<Q~aGV>!h4Vk7M?
zrp;1~*$uLGRa-Z&J#T&Fp2KOe>G84nnlh^2|I?1OfJ2?N@nObg$fc%4lv^t#)~GgO
zXWSy=wnSpBh8bg+iyebuTx%1_CT)>hqEV<M3aO<`LX%u_iCC0Na+#tF$$W<S3vGSR
z?tc3``|R`mpU3k$?|aUBe(yQIIp=*(WZ?jQ`kO1ADVwBlI<|&p{lQ^Vn&J<}``uQh
zu;A>5K&hDz4&}!JZPjAh-wH+#=hxjsw?-<jXmC2?-0_>;E*QJz)T*qV41sfBBV1{8
z?RO`}uqyReZ`HYM9ZG(qEp{^Zsn-i9v$q#4yGo?jLyg?`d)~I%QX1=3M(<VF271uv
zomYQEskGoS-mnCjc01)F{O1;2wi3C%du@qY+G?6~OFlEbIB)cb0<<O|@YP_N^W=!!
zUh*(==EWJGQpI}D9<20G`XJlJOKDorzr3dg<n}hz{+AfP(c=laKG#bkH$O_-j_YM@
zyd-P1iW(qe)M6m(V!~LB%WE02`ekib9O=rWDsD$B<c7VDxu8&joaHgA{bpv(>l@E%
zLX579iim#&SE4nKKHa!SAxxB*q$85l>76#T8n2uAebd!pkM*H<zhgS2_IJaz;lb|>
zpR2v5zO1%>EI($Lr~2@q9YoqnaAmrd@Jrp=ys5G76{d%();epRG)ZVtUtgu0-BK)J
zicwFy_0X2T5oM`%F_6<>feK&Od+Phgn3GMX_jn&UPkG<JnLn#u!#h%Qn+6juSr(HK
z*O8@annN(Y)XnyO;rqsYgU$rfHvdq{*k#m2Ry$68D=J^}{r%mqcA^@OW2#N%V?c3I
zPtsk#&J)hlc_WV?6W67Fcpvtwk&K?=&BF#d-ev=F{-oCi=iy+!C^5CYh=^u!i?YyD
zx2?a?bh#Jzylo(}Ih3w;qd*+gCQE0y6zpn_B5r1C4%F1fJozcWZ}e&PxrvKH`LMLe
zsnDKdoEH+Z<qfA~QwJj5Ch;!PYi=nkcP1kFHXVELNNZG1Q)|+)LY~z>t&EqK$x-!e
z<5#y0`-RNDT_hOy*6b_o$0y&rqEQ~t`3=mTA3_Yq70S2lSa+Z=x7Kcm6;yCuRnI_x
z)YRe3ZcxnluJ*i25#6}MTlaL+@X<QE;<q8tm~(o+!Oh+>BI-v{?OxhPVTRV0^Ff0q
z8anZ1{1Ip=<C`s&rx|BMR#_!PNndI$J?S@8Q>d#?x3k-fwi(xVmVZbwex3}qAT&qU
zc1M?X2m8Zb6j%F%syCe|_*OQiuLH5Yzw({^(VJ@qH(!M1F+_!`^IMYey;D2WtVhi9
z>+am_zMwb!U8~$)n|hQ#I$^c`ZX1yvhbxM?JItGvo)sVMHWxopGhq;*w!+kG+~Z#M
z@f5m3kYl_kD4LdapXeI>GG2MlkK4}eih3LuZPMo1<a^Ux!D25gYa=nM@`Hv`-l()p
z4xS>5+!R{nD!Of;Hc@{Dm1OrMClsDjWAd$TA91%F;t`|VqHedx>%`-E*<muLKzC@K
zNhK?t$KOC2EY6<3U&>)zyG5yo`70J~Da#(9bL9HlPhxH`1u~1k`0VI18CP6+5i{bV
z=acJx1=(uxZ3&0m{4hBT%O9+FhO8@6j_-3m;A;>O9p$i<RmNI1Ek2U;rbE_!b7Rel
zt%w^lQVn2od`fhRS~dfe9|(g}^4rF@U}*Xi@#@;ojwt&kWer}0DT5wd4(^l=JU7%H
zJ9=fmKRU3mgYPs9c}$zsm3?PwjqYUQa7dR*pnZ3Bh6z)?3+3>L`fd`L%Am<Goz{q*
z*GO8EM2fC55w9%u`5~7=6KG|37HFpD@xN+JneLC&*kUQLiJ>rfHS^erXW@q0Jth=K
zEe&<n>!_)_P9YP9nLj}ElpMBQe`!c(VLJ~KV(kOp4eII}i+qpV&<U$GsUSl5LJ;Ad
zS4oL^dRHU+B_c1xGwYK3!+P%Oml#)=&pdlNqWW!?W81YVccYF~6;E`a?&ZBGi5QBt
zLY=-F=h@LPRw$e9ST#rt8Z7n2$8DS%$+G6OcGWwj-et*$TZU663wLceKuz4`v%SN|
zTWcN7u>H=NCt{U}>pP|+(2WmOj1o)xD&PI6EPPu<IHk)b<*KAr0d8X`?ogWZp*{n#
zJ}Dbv69p;6mZBELD50~@AHqo^gtX+$yH<i)`(n$tr;s8mx~Cs2x{Y<!S_PBWJBsXp
zjqYK_`4%*)l*qP@Mh=hru03I&>L1&fEHI+&G0tLEZH?J;;6v>GLdjf}K&h42vabmb
z91sht?;7}~y<;6T*|yY-D7UWxbjj6aHzPpex2&)@NoKJ4@nKzUwmj;kDr{Ij{RDQy
zL)W4XL$TWI%V{vW=fNiXnH(D__Z|b0o5jwdqWwm`@lhFy5%LXek$puX7b;%}XDP4A
zjtAw-1yT>9E>E+R+bkW|4yeh<t~LN0qMn*|M2h_AsjnWfA7x@yZ(cW*6IT2b6bd~@
z{y6~3u$+vo@WHsqhL<POv$HU+{$AL0jWKKSeNy`}w08yxG#orkh)XDIve=Sd-g_eC
z4E35IYq&N!wMR^NuhLlIPpQ{!YaEd)Y3lB;;3eeO$A6s+_hVIreRa4XSfqxzJQ<u#
z@W8P7ByUbaxXJC}4s>p&%2O^j@ayUEEDL!n1{;7UEK@=~@l-d!Ftf<lLI<>;CJ{}g
z;AzX<d64~Q)A^610D~9?%gy)QfKf}6S3i6(*3A>`K_yeTIkd=}ZSGsrPiQv#cPIN{
z7j@y|%{4QR5XxB^xd^W<%9ATNOCu8F>BSf6B{=ibN^(AhvXb-vxfjpPwL!Qx$vGM>
z;S`T5I0DfOO|&KBfhx=z>+VlL1AR>~IIP{gBYwaUKTwysV@V!fZrpm0r-2age#>h;
z(OGJFQ1HG~Jel-AK*9|R%q_JbNKD2V&WZ{U5>vjBAHtrMB^C}^lCzR|%^jlv`Z8)N
zezuzA@s{LI#8U7cUv@q!s)^V?D_NU{E7*V5`CO3*LC_ivL8ZAe2->Azi}pfT2wH$|
zH>2AFLCd&tFft$)01?ie2PBmFIr$k7qVpi!gQ&A5?-x)cq+O|qMI;~+4%i488lOgZ
z&7%N){KP=>rx!pg?!JP7@krtMPJl$PLCe0+EG3n<l;K*cNt}v5y_CSwI<xHWzv$;c
zfsvRCM8f?mv!Dw4By74t0LZ^W8t_N5#K%rZbLG7lZ1rDbB6xw;SY?)w4HxkgN`d{S
z83;I?r!nwbL1G^G^|t%|)xDUz?!Y|afI+|=^N~#8WeoWp-4XzQ{A~!#0WEn7;(S<k
z3@*4NU<0s(ua~s6m$cx#6<Z$E=K=tN{Z|F`IUV1>F{ppyPAwTN2F%NQWaZweJWxbD
z@UI{di@~FTmJpsw(elH30M{%LOLRjM2xLGZi1S!v79hFJ_6M$Wk{b>%h)~H?G-2@(
zAUF%apFsUMDUnP89#SET`ti^D`2Dq+UbfsoUcftUh-iw>AKHLV9&Vul1Uv~lI}pwz
z2iG|NPTP}qFTaCszR<4!y9mHpwuQi%y7=|^Q{;kpBF>*e;Eh?xrEoYeGKJ))=}W>b
z8zAh(3zZis9}rmuIO*rrP({tZKQg{A1;l3u07617o%GsW?3V`+urz>$tt$WLZ53$d
MKT7~B?luJh0B5H7od5s;

literal 0
HcmV?d00001

diff --git a/.Rhistory b/.Rhistory
index e5d50ea..e9a3b27 100644
--- a/.Rhistory
+++ b/.Rhistory
@@ -1,30 +1,512 @@
-library(devil)
-setwd("~/Desktop/dottorato/rdevil_project/devil/")
-rm(list = ls())
-library(tidyverse)
-library(devil)
-library(nebula)
-data = readRDS("../de_analysis/nullpower/test_data/pb.TRUE.bca.n.10.ct.1.fc.0.5.csv")
-metadata <- data$meta
-input_matrix <- Y <- data$count %>% as.matrix()
-design_matrix = model.matrix(~ tx_cell, data = metadata)
-use_size_factors <- FALSE
-glm.fit <- glmGamPoi::glm_gp(Y, design = design_matrix, size_factors = use_size_factors, overdispersion = T, verbose = T)
-devil.fit <- devil::fit_devil(Y, design_matrix, size_factors = use_size_factors, overdispersion = T, verbose = T)
-plot(rowSums(devil.fit$beta), devil.fit$iterations)
-cor.test(glm.fit$Beta, devil.fit$beta)
-cor.test(devil.fit$overdispersion, glm.fit$overdispersions)
-# De analysis ####
-glm.res <- glmGamPoi::test_de(glm.fit, contrast = contrast)
-devil.res <- devil::test_de(devil.fit, contrast = contrast, clusters = metadata$id)
-contrast <- c(0, 1)
-# De analysis ####
-glm.res <- glmGamPoi::test_de(glm.fit, contrast = contrast)
-devil.res <- devil::test_de(devil.fit, contrast = contrast, clusters = metadata$id)
-de_genes <- 1:25
-not_de_genes <- 25:500
-xs <- seq(0, 1, length = length(not_de_genes))
-plot(xs, sort(glm.res$pval[not_de_genes]))
-plot(xs, sort(devil.res$pval[not_de_genes]))
-# install.packages("devtools")
-devtools::install_github("caravagnalab/rRACES")
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+#geom_segment(vline_df, mapping=aes(xend=-Inf, yend=-Inf)) +
+geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+vline_df
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=-Inf, yend=-Inf)) +
+geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=-Inf, yend=-Inf))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=-Inf, yend=-Inf), col="black")
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black")
+vline_df <- dplyr::tibble(x=-3+sqrt(6), y=.25)
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black") +
+geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed") +
+geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-2), col="black", linetype="dashed", ) +
+geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dotted", ) +
+geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.1 ) +
+geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+#geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+vline_df <- dplyr::tibble(x=-3+sqrt(6), y=.25)
+vline_df2 <- dplyr::tibble(x=+3-sqrt(6), y=.25)
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(vline_df2, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+#geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+#geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggsave("gghorns.png", dpi=400, width = 5, height = 5, units = "in")
+imgurl <- "gghorns.png"
+hexSticker::sticker(
+imgurl,
+package = "scDEVIL",
+p_size=18,
+s_x=1,
+s_y=.95,
+s_width=.75, s_height = .75,
+h_fill = "gainsboro",
+h_color = "black",
+h_size = .5,
+p_color = "black"
+)
+res <- lapply(1:2000, function(i) {
+R <- sqrt(10)
+x <- runif(1, -R, R)
+y <- runif(1, -R, R)
+if (x**2 + y**2 <= R**2) {
+if (x**2 + (y-1)**2 <= 7) {
+return(dplyr::tibble(x=x, y=y, col="out"))
+}
+if (((x >= -3 + sqrt(6)) &  (x <= 3 - sqrt(6))) | (y <= -sqrt(6))) {
+return(dplyr::tibble(x=x, y=y, col="not significant"))
+}
+if (x < -3 + sqrt(6)) {
+return(dplyr::tibble(x=x, y=y, col="Underexpressed"))
+}
+if (x > 3 - sqrt(6)) {
+return(dplyr::tibble(x=x, y=y, col="Overrexpressed"))
+}
+return(dplyr::tibble(x=x, y=y, col="significant"))
+} else {
+return(dplyr::tibble(x=x, y=y, col="out"))
+}
+}) %>% do.call("bind_rows", .)
+vline_df <- dplyr::tibble(x=-3+sqrt(6), y=.25)
+vline_df2 <- dplyr::tibble(x=+3-sqrt(6), y=.25)
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(vline_df2, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+#geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+#geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed") +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggsave("gghorns.png", dpi=400, width = 5, height = 5, units = "in")
+imgurl <- "gghorns.png"
+hexSticker::sticker(
+imgurl,
+package = "scDEVIL",
+p_size=18,
+s_x=1,
+s_y=.95,
+s_width=.75, s_height = .75,
+h_fill = "gainsboro",
+h_color = "black",
+h_size = .5,
+p_color = "black"
+)
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(vline_df2, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+#geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+#geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+geom_hline(yintercept = -sqrt(6), linetype = "dashed", linewidth = .5) +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+hline_df <- dplyr::tibble(y=-sqrt(6))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(vline_df2, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(hline_df, mapping=aes(x=-1, y=y, xend=1, yend=y), col="black", linetype="dashed", linewidth=.5 ) +
+#geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+#geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+#geom_hline(yintercept = -sqrt(6), linetype = "dashed", linewidth = .5) +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(vline_df2, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(hline_df, mapping=aes(x=-3, y=y, xend=3, yend=y), col="black", linetype="dashed", linewidth=.5 ) +
+#geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+#geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+#geom_hline(yintercept = -sqrt(6), linetype = "dashed", linewidth = .5) +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(vline_df2, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(hline_df, mapping=aes(x=-2, y=y, xend=2, yend=y), col="black", linetype="dashed", linewidth=.5 ) +
+#geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+#geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+#geom_hline(yintercept = -sqrt(6), linetype = "dashed", linewidth = .5) +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(vline_df2, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(hline_df, mapping=aes(x=-2.4, y=y, xend=2.4, yend=y), col="black", linetype="dashed", linewidth=.5 ) +
+#geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+#geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+#geom_hline(yintercept = -sqrt(6), linetype = "dashed", linewidth = .5) +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggplot2::ggplot(res %>% dplyr::filter(col != "out"), mapping = aes(x=x, y=y, col=col, size=y)) +
+labs(x="", y="", title = "") +
+geom_point() +
+theme_minimal() +
+geom_segment(vline_df, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(vline_df2, mapping=aes(x=x, y=y, xend=x, yend=-Inf), col="black", linetype="dashed", linewidth=.5 ) +
+geom_segment(hline_df, mapping=aes(x=-2.6, y=y, xend=2.6, yend=y), col="black", linetype="dashed", linewidth=.5 ) +
+#geom_vline(xintercept = -3 + sqrt(6), linetype="dashed") +
+#geom_vline(xintercept = +3 - sqrt(6), linetype = "dashed") +
+#geom_hline(yintercept = -sqrt(6), linetype = "dashed", linewidth = .5) +
+theme(
+panel.border = element_blank(),
+panel.grid.major = element_blank(),
+panel.grid.minor = element_blank(),
+axis.line = element_blank(),
+legend.position = 'none',
+axis.text.x = element_blank(),
+axis.text.y = element_blank(),
+panel.background = element_blank()
+) +
+scale_color_manual(c('Underexpressed', 'Overrexpressed', 'out'), values=c("#d1d1d180", "#d1495bcc", "#00798ccb"))
+ggsave("gghorns.png", dpi=400, width = 5, height = 5, units = "in")
+imgurl <- "gghorns.png"
+hexSticker::sticker(
+imgurl,
+package = "scDEVIL",
+p_size=18,
+s_x=1,
+s_y=.95,
+s_width=.75, s_height = .75,
+h_fill = "gainsboro",
+h_color = "black",
+h_size = .5,
+p_color = "black"
+)
+imgurl <- "gghorns.png"
+hexSticker::sticker(
+imgurl,
+package = "scDEVIL", p_family = "helvetica",
+p_size=18,
+s_x=1,
+s_y=.95,
+s_width=.75, s_height = .75,
+h_fill = "gainsboro",
+h_color = "black",
+h_size = .5,
+p_color = "black"
+)
+hexSticker::sticker(
+imgurl,
+package = "scDEVIL", p_family = "Helvetica",
+p_size=18,
+s_x=1,
+s_y=.95,
+s_width=.75, s_height = .75,
+h_fill = "gainsboro",
+h_color = "black",
+h_size = .5,
+p_color = "black"
+)
+hexSticker::sticker(
+imgurl,
+package = "scDEVIL", p_family = "arial",
+p_size=18,
+s_x=1,
+s_y=.95,
+s_width=.75, s_height = .75,
+h_fill = "gainsboro",
+h_color = "black",
+h_size = .5,
+p_color = "black"
+)
+hexSticker::sticker(
+imgurl,
+package = "scDEVIL",
+p_size=18,
+s_x=1,
+s_y=.95,
+s_width=.75, s_height = .75,
+h_fill = "gainsboro",
+h_color = "black",
+h_size = .5,
+p_color = "black"
+)
+hexSticker::sticker(
+imgurl,
+package = "scDEVIL",
+p_size=18,
+s_x=1,
+s_y=.9,
+s_width=.75, s_height = .75,
+h_fill = "gainsboro",
+h_color = "black",
+h_size = .5,
+p_color = "black"
+)
+hexSticker::sticker(
+imgurl,
+package = "scDEVIL",
+p_size=18,
+s_x=1,
+s_y=.92,
+s_width=.75, s_height = .75,
+h_fill = "gainsboro",
+h_color = "black",
+h_size = .5,
+p_color = "black"
+)
+setwd("~/Desktop/dottorato/rdevil_project/devil")
+setwd("~/Desktop/dottorato/rdevil_project/de_analysis")
+setwd("~/Desktop/dottorato/rdevil_project/de_analysis")
+readRDS("benchmarks/devil.bench.20.rds")
+d <- readRDS("benchmarks/devil.bench.20.rds")
+d$expression
+d$time
+d$expression
+d$min
+d$median
+d$mem_alloc
+d$`itr/sec`
+d$total_time
+d$time
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 856c5c9..d23258d 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -37,3 +37,11 @@ compute_hessian <- function(beta, overdispersion, y, design_matrix, size_factors
     .Call('_devil_compute_hessian', PACKAGE = 'devil', beta, overdispersion, y, design_matrix, size_factors)
 }
 
+compute_scores <- function(design_matrix, y, beta, overdispersion, size_factors) {
+    .Call('_devil_compute_scores', PACKAGE = 'devil', design_matrix, y, beta, overdispersion, size_factors)
+}
+
+compute_clustered_meat <- function(design_matrix, y, beta, overdispersion, size_factors, clusters) {
+    .Call('_devil_compute_clustered_meat', PACKAGE = 'devil', design_matrix, y, beta, overdispersion, size_factors, clusters)
+}
+
diff --git a/R/variance.R b/R/variance.R
index c38f807..4e72a47 100644
--- a/R/variance.R
+++ b/R/variance.R
@@ -1,49 +1,6 @@
 
-compute_scores <- function(design_matrix, y, beta, overdispersion, size_factors) {
-  xmat = design_matrix
-  alpha <- 1 / overdispersion
-
-  mu = size_factors * exp(xmat %*% beta)
-  residuals <- (y - mu) / mu
-  weights <- (mu**2) / (mu + mu**2 / alpha)
-  wr <- as.vector(residuals * weights)
-
-  xmat * wr
-}
-
-compute_clustered_meat <- function(design_matrix, y, beta, overdispersion, size_factors, clusters = NULL) {
-  ef = compute_scores(design_matrix, y, beta, overdispersion, size_factors)
-  k <- ncol(ef)
-  n <- nrow(ef)
-
-  rval <- matrix(0, nrow = k, ncol = k)
-
-  cl <- 1
-  sign <- 1
-  g <- length(unique(clusters))
-
-  for (i in seq_along(cl)) {
-    efi <- ef  # Make a copy to avoid modifying the original
-
-    adj <- ifelse(g[i] > 1, g[i] / (g[i] - 1), 1)
-
-    if (g[i] < n) {
-      efi <- matrix(0, nrow = g[i], ncol = k)
-
-      for (j in seq_along(unique(clusters))) {
-        group <- unique(clusters)[j]
-        mask <- clusters == group
-        efi[j,] <- colSums(ef[mask,,drop=FALSE])
-      }
-    }
-
-    rval <- rval + sign * adj * t(efi) %*% efi / n
-  }
-  rval
-}
-
 compute_sandwich <- function(design_matrix, y, beta, overdispersion, size_factors, clusters) {
-  b = devil:::compute_hessian(beta, overdispersion, y, design_matrix, size_factors)
+  b = compute_hessian(beta, overdispersion, y, design_matrix, size_factors)
   m = compute_clustered_meat(design_matrix, y, beta, overdispersion, size_factors, clusters)
   (b %*% m %*% b) * length(y)
 }
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 29765a2..ba2af52 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -144,6 +144,37 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// compute_scores
+Eigen::MatrixXd compute_scores(Eigen::MatrixXd& design_matrix, Eigen::VectorXd& y, Eigen::VectorXd& beta, double overdispersion, Eigen::VectorXd& size_factors);
+RcppExport SEXP _devil_compute_scores(SEXP design_matrixSEXP, SEXP ySEXP, SEXP betaSEXP, SEXP overdispersionSEXP, SEXP size_factorsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Eigen::MatrixXd& >::type design_matrix(design_matrixSEXP);
+    Rcpp::traits::input_parameter< Eigen::VectorXd& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< Eigen::VectorXd& >::type beta(betaSEXP);
+    Rcpp::traits::input_parameter< double >::type overdispersion(overdispersionSEXP);
+    Rcpp::traits::input_parameter< Eigen::VectorXd& >::type size_factors(size_factorsSEXP);
+    rcpp_result_gen = Rcpp::wrap(compute_scores(design_matrix, y, beta, overdispersion, size_factors));
+    return rcpp_result_gen;
+END_RCPP
+}
+// compute_clustered_meat
+Eigen::MatrixXd compute_clustered_meat(Eigen::MatrixXd design_matrix, Eigen::VectorXd y, Eigen::VectorXd beta, double overdispersion, Eigen::VectorXd size_factors, Eigen::VectorXi clusters);
+RcppExport SEXP _devil_compute_clustered_meat(SEXP design_matrixSEXP, SEXP ySEXP, SEXP betaSEXP, SEXP overdispersionSEXP, SEXP size_factorsSEXP, SEXP clustersSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Eigen::MatrixXd >::type design_matrix(design_matrixSEXP);
+    Rcpp::traits::input_parameter< Eigen::VectorXd >::type y(ySEXP);
+    Rcpp::traits::input_parameter< Eigen::VectorXd >::type beta(betaSEXP);
+    Rcpp::traits::input_parameter< double >::type overdispersion(overdispersionSEXP);
+    Rcpp::traits::input_parameter< Eigen::VectorXd >::type size_factors(size_factorsSEXP);
+    Rcpp::traits::input_parameter< Eigen::VectorXi >::type clusters(clustersSEXP);
+    rcpp_result_gen = Rcpp::wrap(compute_clustered_meat(design_matrix, y, beta, overdispersion, size_factors, clusters));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
     {"_devil_init_beta", (DL_FUNC) &_devil_init_beta, 2},
@@ -155,6 +186,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_devil_conventional_score_function_fast", (DL_FUNC) &_devil_conventional_score_function_fast, 7},
     {"_devil_conventional_deriv_score_function_fast", (DL_FUNC) &_devil_conventional_deriv_score_function_fast, 7},
     {"_devil_compute_hessian", (DL_FUNC) &_devil_compute_hessian, 5},
+    {"_devil_compute_scores", (DL_FUNC) &_devil_compute_scores, 5},
+    {"_devil_compute_clustered_meat", (DL_FUNC) &_devil_compute_clustered_meat, 6},
     {NULL, NULL, 0}
 };
 
diff --git a/src/variance.cpp b/src/variance.cpp
index 6a8028f..584fe1a 100644
--- a/src/variance.cpp
+++ b/src/variance.cpp
@@ -28,3 +28,36 @@ Eigen::MatrixXd compute_hessian(Eigen::VectorXd beta, const double overdispersio
 
   return -H.inverse();
 }
+
+// [[Rcpp::export]]
+Eigen::MatrixXd compute_scores(Eigen::MatrixXd& design_matrix, Eigen::VectorXd& y, Eigen::VectorXd& beta, double overdispersion, Eigen::VectorXd& size_factors) {
+  MatrixXd xmat = design_matrix;
+  double alpha = 1.0 / overdispersion;
+
+  VectorXd mu = size_factors.array() * (xmat * beta).array().exp();
+  VectorXd residuals = (y - mu).array() / mu.array();
+  VectorXd weights = mu.array().square() / (mu.array() + mu.array().square() / alpha);
+  VectorXd wr = residuals.array() * weights.array();
+
+  return xmat.array().colwise() * wr.array();
+}
+
+// [[Rcpp::export]]
+Eigen::MatrixXd compute_clustered_meat(Eigen::MatrixXd design_matrix, Eigen::VectorXd y, Eigen::VectorXd beta, double overdispersion, Eigen::VectorXd size_factors, Eigen::VectorXi clusters) {
+
+    Eigen::MatrixXd ef = compute_scores(design_matrix, y, beta, overdispersion, size_factors);
+    int k = design_matrix.cols();
+    int n = design_matrix.rows();
+    int ng = clusters.maxCoeff();
+    double adj = double(ng) / (ng - 1);
+
+    Eigen::MatrixXd rval = Eigen::MatrixXd::Zero(k, k);
+
+    for (int j = 0; j < ng; ++j) {
+      Eigen::VectorXi mask = (clusters.array() == (j + 1)).cast<int>();
+      Eigen::MatrixXd ef_j = (ef.array().colwise() * mask.cast<double>().array()).colwise().sum();
+      rval += adj * ef_j.transpose() * ef_j;
+    }
+
+    return rval / n;
+}