From 7a01b27cf4d3450e9abc8a93fe8f9c7ae5f844b8 Mon Sep 17 00:00:00 2001 From: pjentsch <pjentsch@uwaterloo.ca> Date: Tue, 16 Mar 2021 16:34:13 -0400 Subject: [PATCH] documentation! --- IntervalsModel/simulation_data/hh.dat | Bin 32467 -> 32467 bytes IntervalsModel/simulation_data/ws.dat | Bin 31333 -> 31253 bytes IntervalsModel/src/IntervalsModel.jl | 16 ++++++- IntervalsModel/src/hh_durations_model.jl | 27 +++++++----- .../src/interval_overlap_sampling.jl | 31 +++++++------ IntervalsModel/src/rest_durations_model.jl | 19 +++++++- IntervalsModel/src/utils.jl | 41 ++++++++++++++++-- IntervalsModel/src/ws_durations_model.jl | 32 +++++++++++--- 8 files changed, 129 insertions(+), 37 deletions(-) diff --git a/IntervalsModel/simulation_data/hh.dat b/IntervalsModel/simulation_data/hh.dat index eb57911038e6bf2d99523665eb316df80e6f8ec2..fea7ea8ac0d3dd37a0ab30ff5ce899e8d56515ce 100644 GIT binary patch literal 32467 zcmb82X;_x!w#PwnRs;oeL=+rQQ50oThT8!o1rZfR96-Rt0TmEM6W;fE9B>}eU20ic znUzyHHf)<LGt+E0Sem729<$OWo73~O_y2cJ*Lz*(!|_A^Ypr|T>t4fut><Nb+0l0o zZRF$S<<-Z}yTP!Mg2Ix*f`LW(Wo;Xbzg6_uwWn{ad2vF?giB*yf8M)^jt9(~KD{u% ztY~I&>A)J!O|K??-u2V-W)v2ToHJ{BVcYaA|5gXo^!Cn}=w0KTgMK6a($#b<{W3W9 z^RAbXS5j7#KfSQDZG*S|toU%obKgc>^RIXN?GsDuc{fX&SzJ~)Ag^Tl%(TM1(m5rC zGYX5#^j_VbYrVa^F8;4GFW-w^zJXq~Z~hg{-~Tyg|F3_6UQI4G@%Hxe^RAOIv#7Lm zW^vmFJz~!6DysPENa;xTHox>&dwMq-SX5e8QZ#AKe@gJz1A*Oqye{^3{u9`(uGh8R z&VK^E{u5=j-j5b03sZ#Mh0TSdg)M}OgaN`SLhf~uoFQx~%oJijR&q1pUBdZ7Kj8vl ze<5)YANP27g6xPdM{-MHUttpu{XLTLPrMyG>^NU4eYmiz&<N9nA;L*Qx3IO4c-u<O z7LFHA7kb7|TvK%(Cu}49i`{?y8~h*eqRzaBT|1o<XOf3K^&KNUarcu<AIP_^^g+Tt z!frz1D3*M;u!oSmdU~7_FZZK#j$i6SU#oRqDTJTxC9f1RkztbS3CWB4$4d`?;4$@R zrE}t`Et$CB3H-r-Z`qL-Jm@HWu8{sDNDdW}4}Re>aiOP=^naG_4;01-;XQl?c@KZY zPoBi>lAp1{@j~uJN)8kb@Hi(<>dC(0y;RxL5A^s4;Z>IGiGw<_UzX^cKKV*6_0aP! zc_!%`zGGiUdiEavK>Wb76|xT&vL|9Cj}uZC&Z!T5oF+T;#Kk>$nIpS<h1AbSa<s6a z$2oSZrN?h0$z{T1A^Ecp;WPcAFYt?gaL#*k<&XH_VTAP9kvIE>{K#jq>~n?er)<f5 zFPKO6Vke#V5|TH3;T#?hk=;xo`GEbTCx5<M?3oIkGuN#p6DRQ%OFu$LT+9VBbx)Mt zP$7GOxOk6#!BgJD9{(}&L;lzoc*MgTk_X=j;$vTPP90{*4|Qi>43<7!NdG~;ANU(C zJMv`DW=QWZBoFFEo^5nKP{=#%0rE)Fd0!##!IMJisV99S9xzsR#E<_l>E{chgv<ea zW^QZgJW`1NB+10bcXyHW<UK|*@z<8joH4iX0{`%bzL5{#^Gx}v64HO($(6pdkp04( zBuif=j1!VKabnLtpgwbCA0xy*K{9zz?|#x-Lgs`$6fQk|f*;g}K65`>_Vk5�Ni; zWJmo{BojA0qaXO;{=KqC=DSDV$rE|J?B@!J6Mj)=Pd}q&Uo9j)>PFn%b^iCiMbNQp zePu%}s1a{cTew2*$Z(EWiKV0b!UcA9fb`UWnIjJD*$uoye15vmP9;9#?yqxdKwNNx z9`G*rLEfW2+@sg*AnsES@`=^^*z*m5JIp@!ntId&J+nX`s59@t5q1Z&7$|>9LheOL zhMUY1?~y<82g{y3m~H%f#zEcaAAN?C$n17zm>oYz?-56|Wc+fUyy+|P(|6*g&*UBB z;SWCX9(mQ19e&A^-*&{~`OSkD)QdVGGY`B6QV;Cln~&n)ea^W@9OOyf)RlR{o_e8Y zzTg4(=`()eIew6d3#4A~jryYZ%yWqH!!O9Lr=BsgBVYU>(+7|_Ab#pje~FhpLS6Y* z^6ldt^2Lt%B_H@r{0$T@e8rx;*emn_Kj`TnaZ-oI^4ncVT-1Sj!*}B4UG@m}(H?$+ zB_{}pn>^S@@S1wUSL%%aZt_DuX_BcI`Sq3_z6VQYe-j7!#_1g1(kJRq{6==@L3j^O zsY8tHd5`Z2b>kiS#cw6<lQ+Cel0W!Fy!a)5>PdX)v9Iru2mOQxd_T}L2k<6V_u&KY z!x!GAe<1nLcTc|ZF8ARBd6NhIqAvXv4}Q5%e8}wi!LnyQ$P@pmI!A^lG18L<h=1&n z;VXTne)K<C{-_(g#xF>o_{Sf7<eYcm6LqIP@R)qzt!JIcx3=PAuM$tN^z?^(;1_kq z4|?7qU-Wz*Jm29x>Ox<z2YDZR-esPcha|;Ko#;RN$Wu@K=p(#=r__PI5Dy3+qx3F$ z^A7ikn{(<HqkEk54#<7-U~k6h9`A6DANCYHfd}N7uKW1I4|!29>JFb0bq_u9agR7S z=RM-V5AU&$(DQwzPVAk(|1GjeF~Hq{l9`cA$xDUV!YUzLW5y;*KUIiddRr+y@6t2w z!+HMq31-Ge_uyz9$#aD>g|={&u#vD(I88|W<0a1%h6)=CYYC|XyPfyRced=P8(bh? zdI}e^WMAO%j*-kg{0@-5uW+r$y)}}lH}4V82%Y1HnPfJoZw=YOF=nHo^oc_1?k#zl zkZ%-q9w<G1a7i95geOBJyM@HxKr-<%d-Q?#>dP)($gB>R%#2`{AwBP~D;G--ui*nc zh}HQZ;TR$Pr2g=jI*yYab&i${@0n5f2v3+<?yt~2;_2<-m$?3~IEV{gEZ2Sb4j<q* zb%dA9EB#{EP+$5@eArPZ`UuZ>r-$NTXEFb29{KXlO6lPhec~RwoBC2u_{zO{@`FA7 z@w6W)JK`j6;wjO2f{?m#UL^f$VX!b>$R0_SOh4L5rVi|1{C3cJx{!J4Dw)2K7x^#H zxhGG;q^BOkB=;2J=WfaH2p%$@><!+7S4(t{K2MOGFJ!;L)1lHM^P7TqxQ`qod*Y-X z%q{t`7pBXeI<Oyzk31@5N1w?DfAHFpU4d|vaFUR^hfC)BHcc}8nkSjMkPm#JKljRR zwvc^D-{Bp<>Db@sd5`+mlOOV1D0!@Kyl{n(cc@zv>C1)i6#wMQ`|J<q7QT$o{Ujm$ zt|>WBi2n%5(L(A?U-{jzSm#rO9feWCTp@AN_aV}E6;dbc*%#R7$PWL+(@lE#Q6`ys z(!X}nGhZNmUZZpR(Ni*YCQrUUygN&F#Dib@)=%fTLgr(PWb$kvdAKlJxKv2q%q8`> zN9Xf}mXLU}B)1h3Km2Pi{c0gRz+WBd;X8YNl=OU0*{9@39jG7ug(vK{Opo`wOQv7! z4fZl|xMVlT<39Xk4yhA7h8M$hkG}DHAz1nlA$?dWxxcVnSXW42$d~#o(>Z&B{tzej z%oE=M?v2#_MZ#LbC?R#Gj(`7Kgf~2Q9dCk0zH)`^2F@ExkL<Z4=s~dTIET}L(z8>r z8zMb^=`~y>{x-6MyAvhz4*uXG{<_GnwTFF*<Q76^0DoY-$360B<#8WB#7X@4VfP|~ z$?`)F>3yv98A9^uD4E^FZV8sYr7%p$H<6hk4&H<FO=Ztc!!I)NlIIxNqsJe0^vplO zBOdbLef*Lq{TZlx%vKZ0yhA?h26jx8?6QRPAHKmW;^tlEq^s`L5K>o={<9+n$PRw6 zyRfHDyaQ(I9)7SR-e8@R5A#I6*i#Scn5cWiNu8;uCr^n7f5eTRodYk37eB<seDfZ8 z;5SZj5$90J{vP^h51IbwOV7I??=dgDm+Rq&Z&w}Zi4Q-a()Sg@(^$#y!82awmb$~! z!MaC%c%QmbXOQ0l#81D#NcrPFe523oZR*Va;~abHNgtR4WPS?~2YSBo%r|;|Go&g$ z{AEiXA|!tNj+UM}!)x-3);WCRTgy3h<2~XbK6uM}>_xtJ)C1m7C%(DxAxm)#5Ka)X zZ{QR0@Y{p@*rV(PWFPtAeP7AsgFoUVzqPWfFGP<$b%fXW<2O%?$35o9PkQPLa=)g| z;SY7?UG@;aWjHU?y|F_21K+88ADzQL{NbPakUxB+{_u~uLGr^d^F-ZwpSm!&s}&FP z&3jSO<3CF>`O;_l>#cL*pD3C9;2FGwXZWN3>=ojSke^Ipf)GCNK7Ui7hadEVJqXe- z;-l~IiF&~U<`UlbR~-0fPUvHAozoZUNM56KK0(O$8ouXCkDnTnLxf|5>=o)bO?qVV znJqo%^ozWi8~Q-r#1HS_4fm-(^@6AL4Zl<MK6T*T0Umnl2p_X`PX4?HugRbK8QIeZ z`o=x@O&sKzsC$!y_~Tv9xd;CJZ;@!-4HL!-J>Osl#Oj>drbtc_;txBx2U271r|TZy zAZibi12fP|_U(lIgmFUd5f||c)OmLye(?v#=^Zssl0Eg|KJhRU{N{zj{dEtnP)~Lg z@laPdKwf=ypZAHET>^)QAHV!|<J*pZb{Kl{p?<_k+{A&NT{TE?QD5po-QYO!k_UW9 z)_ryaJfQ!)M_%+bSoeAgxyLzq5FdP_&cx4s@&}2VzK81_>M=$#{@5uKrRV-=$<&AU zs1FEl=nMYgEA_(;|J0B7*_C`d$Oj&DRy@SvnO}&`J?~RTc#U6V;-YWk<;RnsU8Sc_ z+^25L9r?iP9NmKlV3PFhh42dAfV@lG5wfR$#7(`?6Cd{U2YxZ9<iq>?c4WTj6Lq8y z)S3911LB51yoZ141g}$-5Ajep&dHBF$cwnBC;ejY;}_nehbNp94|@lD_BnOKKl#8L z>O$S|hkd5<gU8qtH+*BR$%{S^2fU)b!SYKy$i$C-;=(U`3_sMF{tz#7OT5$%fApWY z$eVZB8}y6*Vh>+Jl?U~t5A>aQiHkjr9dVF%d-)@7_(Z;*=gbXpw$**U8`K5t;^CjZ zGN;T9e5GE*O<#$R{NNdN!j8DP4^QY%lHv&#Qg{3_ue?KC>`(d%kH~MB{P7<3#y&{r z)V+^n>cMx1Ipcl&z(;sMeepxw@QL4~)EEEb^e*p!A<|QS{8KmjMI0bJrtaADF7Yy# z#7BPA7yeQY@*<9G#fLsgGW8@r_6YUHj{VB_kT~Ek{batV3wrz#H+{)a9O&UeuJp_c z^}#Mo=kS(#!5{2+pZmniUce9ihR58euJD8YQU~}?UAaeJ@z4A0KjzN!9{DmCNy@X6 z5PtmqZxPONWK&B>%+#cV^hv^I!aBlWVWluh*h0wt5Xrg1RzkRhUrYKxA-U2MzD49j z9!au?OT5ochMRZEu3A`Ac#n`euzTo5L!HkORtV`Wen(3`Nr*qr*$ESL-cmSL2p5`4 zE*H)ak`J>%|8jLs-o(rM6LpT?HIkVr@}=(FhX?Eu9}hopmi!v%oVbuVr@zES-|?TX z`}2gn=Q&gGueIzrpDuZ^u(^=>GQ;qkeCac@P9GNPKKfM2<TF!psW4B-yX?p^>ET6P z$@HH-qlZ__Jo7{R^>p7`*jPwE>02Y|tA+SsZi1yJuO^b&ZG6+GN)K=08FePEEZKDy zQir;dnX?qh^t+AZWkSAT)Hz#v`cNpjNXT9YmW-bPlBx3&$?$HB<NzUa!TBob;a8&M z4np2-C7HcYBDt4vq>#Aa9s4a*=kT5RrVb-?PTtfZSbFlDD0#kcg^>7$NakC?Zwl_y zmu|B271FmJl9?xXPrvy#2g#0nQY9029m%|dU8VHxgyhTn#nRUn;)gyk|D5yPVD8vI z^aH+=S10+UF3i&+>FFEqz|UNr({K0+zsP&8>_!R6(=C}gk{5GWOXqwCs9UV`#4%7Z z^G08p!*ZRo@1rGC|JjmxpE@xQyxT%{%mr};NWWT0{m3(3dhQ2H<~M1gWc=_haWF6R zi@w7<&YAb6@<(5KONI}lB~vHffxpx>M0Q}3WcWb5HKflHrg+$456|J>LfJJGb`rA3 z=?~vK>POwDc)XJ!nfj3jeW1_&vTG{L^YB9*=?gr<AG{_H_&!&DMhT;Z@P|H75Bk^9 z<NZ*{ywAJbkJEXskUi%onfTHrvlrkGyy&BI@?-B&zbKuD3x^7Mhq&0Q@CZF~!u=eN z_lRS%^z=DKGX2Ic^NAh%mHKwqJ>Fwqvq#vAgJegY10=(9;-p^a=_ma%x(9E$S5tcS z8hrly-y+OPl5F8Fza8NavmPuvYBNwWGs3-M>8pgfLb$`X7;cp7oY^dpOmBiEhY5p( z<PO*105c7D=o#-5FT03c&rDM@W|X+*=v`!Hv_g8g3D@B=-+t;vFJg7CnQ*ZXZjpb4 z^s_vWe8``A(|b5a|JXHfXRZ8^59jlxr{3^@`oQ5p*^wuHd2g=Hsncx9#7F+aq-TcG zJp9&^%q)`^bt}<%ijbKtlFV-59pa?V@Qru5&raq({<GwlZ%z}*)Pvo_y+t~wzT{b7 z`tHJ1Au|vErb$o!jU+Q~agwP6b?+}de5c;<gt&R19fY1;vP^#33CWN86DK@K_po0o zd9ug-1j!AB@R&NbmmZmVGUr`&o+N}P%prcM6Lp6_<k3d=Gd=u}hqv_PNj&hJJeXhn z5-;Zy<cB_INTxsNnfG{|j~4PS{+T2C0Ds{db%&<~@<ScwNoL;RBmEz#bK+nx4VJ!C z$os@MM0)tno`Qc<bY3h>6w=RDl7ocodG<#g>Dep917DbL=AAlY&;EjM%rpK=<R4zF zl1zWu57eKzSRgy@Q_q>wGr!T2`wCq`>M}_(aj_5JQ*E7-7jZGyxjLtRya&%#>YRA8 zC5H>?Q-x&eLOm*_r*GV2zNtg5?BH9d<f+1GLgIk`@Q1v)$9@^Ad(>fo<Z>bPhsXF0 z);WAGlH6Sg&)^;PV{XWgJ~5Y_bst_0l1x1G6aUQJTG_=22MT$gK8%;1`*S2$3(3<@ zGWF;$IY2m3NdC-oN9h|2nG^ENk)FK+-$zKlMo664kq7njksbXHlMJ6CB-0Q0?v|c- zb0xD^h%-)l`q4)+b*?3uxC$lTEu_ErA+J$74-pdoFv--N{j*y7C?PyT=6;gS*+cLJ zKIG|~{;=Q2NRL0xr%RtD3>MOt4wBik)R{S)pmX|^BpK`}IaJ6T;=hgbyboWoYoK%L z`S-s?*iq!d%yiKme<8i&4SF(G=g8%f>19pHnL^%SW^L*5Kfpr{l+15H;)OfxKx#+4 z)U<`}gEJ)aeoM*BAT!%adcJwVlG)7(lA8*tTQkY*KH_3FqjbJlNL^M)t}mo+<UdFH zmBJ7qJDgd-Ki@dG!#V!Eb)P(mpE_aZD?4@_?;=xAc1AnddvcqdvrOmsuPvGS#7bTu zq#yLNoAm5}O3B=#Ze68^AG}AL%nm%_H=A4cs2lyHE+sl|BxJYoO@>cpI_FyfZ;2yU z=j20uMoCZHbtSWV$eVt4(m6aMPTrfU^T|TyV}Rt|!f`_G-y@m$s3&}}be`sMk3L|B zzjE0zPw*DrW#~Ll$hU<0)sTLkFh@wg*(=EK7JgCpMY`8VNPO@)LVEg1Kd3MK36))e z5MJ|5=_@_;fG2zlsKad8WeR%;X9|g*d1bDcmt@(|AN+EU-wMQCDtq!F556Jr4So<e z_SAK${NxMSTgxRAFZ*hS^u)n;gSg>gq3oD@?5R(M&goNQ$*ID{Lgp}AGWCMz@R9z* z$9A$`B<vy#7N!e{k9+v*q;vQ`Uo!boCwR_x2xQ-k);;F9sbu;Jf03E%Rk9-w>PP>G zuSj;xH+H$w6EF7c8{)%$sO+gfbAuoH43Fps^_!ymD}=;JzQNM7uX%_4!28T^uI#%D z8wlzDAP*T{F?YlTujoJbN_9U)7$k&8$nb*u>@Vhud4vD(l029j>W$r8y@wt9j{JJ* zoOiKf&(JsOOh4A@9(wjV`M@*u?Eg9*`NMni<NM6KW_i2=zsRGf&f!I@WcW<|*t5h> z+`O|;_uwOSWc~*0yq=If-(E6tQLiZJ^M&~1e3|s*0YBg|{V4Ic-%&E}hD+vNYsu83 zp=9E7OAZyn1M<9EdcH^0mpw_{!emFk%Ot}m;^e!;oG{n)oqj~?K65%uGId~I{rztd zcIzD3#t8chM+@0~Ig*ow><D6=CViFg9w9q}Sc{}5R%Vf!q0f*VwHYEgQ@B6~=N3sG zDI^c>Q4?wh$3kQe7luoQJJTgoAL3$15kI>ZdpI;k_nE1=l8JMwWa>%W3DT2)s$|}o z?QzaM>I=7w?C?h)=s#SfuFNLqaFQ9JPattqH{KbccP9(s2l4cio_FIUQ&(myM|$c^ z-zub^BqSc{N}lYfiL&DzW-Ceh8AABPEYolDA&$Yar!UkiSbBJX9wfiEvV$M=6B+vg z*%1%#f$V7NS|oew1TT1(x+TeuI*=!EMe2N(kR8iTM$desFO@y}gE**ryv|z)=@<U6 zgEwKag9r49{nAY5xkB=wkGwZg=fq2VY0|SFh#xsZ=WT?{2mG8MJ@ZB%=^y9Jd93Wo zFIqBvS|E9`kbZDZe#{N~yNB$Hh4^JIvZQBT@Jk$pI)_K-@lT(a1NK{<?lIrgyIT73 zLiht8;0yl7%8vNRkG(>C*ui_=rBC#eI?a+_`o{MF|I9_D?BEx?>Lxw>p}+W{&(w>$ zMCu;>?<1MHg=g#u?h`-1G2kh4fxj;DM;)jy{pWsX*|A5c1NRbj&U_7&+|ol&AK0s$ z^Dg-$>mGZIe8`W!kT-LTALgo;{Ll~hk3D+<d;GHhxKG{O@<%_ZD}1GIC9<Qw%_TFJ z@D={@4*P<74%I#S1<#ll=FGD{@QgX&J?h_4eyKD4B!2opy@;cO?h)T&$uot-1F!i$ z_x7;MmE2uOKX`|Dd53uB$i9<surNVb=z-{&H{zsT)RVeY>OMS%55z}Y=;0^(10FD+ z<jp&Lf0_4ky_YG(FMKDT5S=FpxlcWqL-N2+2ify3@y?Q-yr?7hnS0Lh$6i3@UHH>S zeyJ1vBu@4={Gd+UgSX^IePZMnJ~DTKUcUeSe`8F)eZ0liQnzX4w`KB?O>Q&&(Mf*< zL|Z1k!TesIzTq-{{t3}*-mEl76MFZ5>L-`k5b@de=B2jzDl7M)^;K?DdMfXuQ(N3- zd1A()oo!vlzL_=h>b<s^xWD}SG+)cKDl&(APjs7M|9s<czsKF?u6q-H{=J>cT&Q&{ zcm7n{6l{I|@!?x+v${vW4LgojneFd>)B5HXm$`hp$<%#2Y||$>;^32y+9u&_mlrk_ z+oo(?-Dk2pxJ{Egmk%dTw9HRu#+=%JpUZ5m6(8_RE!%9koHFH;&n>epcIBO6NwyhX zGUw>o@s>%uH2-2qC(AswbKS;g>)YnLS8ERR|IKZt?P`2DYqDk9CVf$HeuLZi{&LG~ z9&VcttBdkKzG|7>71x70KU!%rPqu2)aE9B2r`=x=?Cmxa7MD!9_mFMkcG{Eo^{~x- zkMI2I*lEjTEpOT9YbV~lu`5%LJyB_XzO{GzV@utpw|AEj`Q>i2ebw@%F&%6(W^S8H zEuw96WBJ~)6O%0S*-I_gZaUyLbz??;+0M!PVz<2f>az<?m(EGU#{B3u5p`~yxIW7= z-#?et;nRb*+1NS#lSh|Yrc>^Q<<kdQ=B;V1Gk!enGC#dq^X0e=wwc<cQM=;qZsXVZ zmyMTiS!PR@4ZBASvrUWp2jzd>#cc}u-E(zFp4%j3hrbZm*)pX&>-X*C>-6dN@ewz? z+-BX_(MPY(ahtaDV%v<*aGB!W&#vhAfy;bzYj3xG`L=m6>*M)ZcUz|L*n!xOesGzf z2fulEugh(^rmnyC?SU#2Rd=21z6{$uwAHK+xLIX()Hv_cr$wdt^W4sF&wc1NuQ%*^ z=tz6FS+ORy;lVem%*>ccnGFgpv+JI<BbE)e&6)ptuQGJ3ZD!6o6tZBwZKmxR^=jBe zhj**|2Zx`iG{avEX|r>&+sx_t`@nx3w9K%RttuMDxJ-+dU)`+TzS11|xN>Gnxn;g- zSNv7P8mEszi&yk|);2>Hz4OGWcP!)kLz?LoV4GeGFJ(VJ!Flhlpsg?MwM_pu@Bdu- znq`ddC&#UeRpw6o!((%bs?4UH>sHhoZkeX-Z>9bEs>}T1v-9~MocLP}e=zj(#<ux% z{d1Et=US%rg!B&@zvedE8m~HUhgqiY?yj$VH_l~xjPJE#%VMYRJ<IxMJM~^aFXh|s zoqBz7V(R_#{cQ8+<A-|P-03oJjthOe=I<`E@#wdIobI^L>}he~fyA}832PL)*XIEz zpS1mz9s0V=rXKy?{Nmz5v!v_3ua<|~ro3v;oy*2DFFdmJg0;bA@_h5RuRT|3+UIWA z5i!S^n<=}pGaqu9OE(Le-h19<YQOpPw?{Iac`3X*;ffR2H|v@We8ri=@;MVa-mGt# z!;kma-r^It`RrQx{%@l!^F+|%lNE2d%+i#xf7JfS_18RA|GdL(YHjk#`siY%dFjty z5$`(v&R=>YZRoRACb#vGhrBMh&A4+-pT2$6GTqmHx_7}BE^{#DQp-_J{q5-F<V8+j z*ALpVCTE>x;@iBxG^~qlCa?JCvUwLQGyBH=-{-8dO~bQ~URhb<HcQUjedWh0=iNQ0 zm!18=ndi6Wo$K1#HlMxx&!i^KKIym9ZdvZc@4x%8_Q!8p=F5~RyJvN(GAFK;pZB|H znV!wAM!uTvGB5wvd#(4+u*}q1hjIe<yUo$75$Wq%SZ2lQ9didCt~CD|e>^Cx$k`u% ztV!80!kM2%gMVEW<Tk~p@_ZLK^LqcFEg=gE-R8v_=e_>*Rh4-r=F5@4ePfw-ecBx< z9AcTEk{jojuCa|PU{Cg?-!0SPcGB#_PQ9#Advo63=){%ZWY4~2+dTeWWPpE&ZAR@l z7N35_ZDM22|6CI8GV!O*g~wmAOvkBL2DF&(G8NeeW?wmDnK$O0uAaTtZSp>Co_Ha{ zZGP>WePY>fRpw00thv=&-Nt?9?#-PV*k+k)_#gSBZFAzcGgogo`{cp9?(g$yj%7~& zbaLvYK+8mhUABL5cw5#kw^wVI%Z&4XeO=@t%alcosoTAsWj2;RZ2oa<q4{y!GnEtD zR+{qo`lg+;kKglZ-TcePEHmZC{zu=*sxqI&%sO?;nWw1D-DgCnS*Bu9)=w|LQ)S9) zjEtD=)GxL7hwtu8u+71<F<GrvTPEb>hPV1RvrR>O{dHeEyc@Zoz5mfg&OC3K+}P_W z+w3Z@{=LDyZu3-+e%5;{oxD%`ME<tdHZ_L~Z*b2WZu4HPW63#dEVJGB`nif~+Y}}r z&N=w8^L=UfkMy(7ylzXmR1!JD`F_-jZ#d>vx4D^9_mRn6T;^8DfkQ20oIUo)_ZMzF z<M_KWYVPiRmN_sgs$grRZB`^^?97~No18aFZ&asNnrYpeJbvEk%eA0}^8@~Ln;TEQ z(7VPt%gnj(M|FcR+g!Ro@8M9VZoQt1I}yIiGUFS)JbaB$l?k~1$duYawmDqd(dWq! z%Ph}Kdm!Ko%TzafEOJ0sx9Rj;T-%e*{o8NUz1DEQWuCrZAFTeD%PdIPz3_L-Hg@sP zM<z6Po5vqpxh8bA+kE$Nt?fSVStfhW=_mSZw9SF4J{<>r?=pb}n|$kDbD2BEwa&f# zhh<{IE=PY;V3}7F&pw!(>Ne->SFUtC>N1Dki|kx~iEW0C>O6k>o0jRbx%=6$&JN!T z-}M>d?=r!O{Tnq-u}$^fU*>s*xXh3DjT_zbY0Hc}^V*$h)7<8nf2|!7aL6*J-hR9G zq?68kK6PUM>Rihl^_|^!-M)p!^-AX8JDXi5<4XRq#GSTz=}5P7Yl>ygywvjg=;;ps znzs1qhp{d*^~U}+ZF20tzAwAI3#&|{DO*4P;Fhxwzwhb&U7sqG8gCEou+%m;)ArXs z@UU&38QwE}&;-kLh>kG}ojKhR_T9yoo%|Blhql`Am1Vl7$DKM{V4HCdmMsZz_|@sa z`rMOl>#y&`6Sc=VzZdeNBOi&d&AQ_cH3)FtZJzV>@SitZX7lELLH^e*^ZeC_Lw;4Z zxf6HjNc~@{OuMklUEg%RJG1W}v3%Vym&ty$I3UdVuKU0EVCQsauisvM;L5W*EpxeN zy^C+%vCNW%H$Qs)yxVm6;L1ajCOi9UtFiVt{CR2nsm~ktvQ5mf@S+RO-rUlC)4Q)< zcbm8T->YwLtTb<Z+A-nn``zYF-QAs6546n5=id2n+jCCbPww8m$@xy-**?<$qaBth zIqlOl?sJ!^pa0YA!<~69NItykk+yF0XvXTU-*0o79lqBazMg8C$n?1Hjyv;k@#6HZ zm#4bT=l9$d5Pz)7Tx`(zNQVfwc{BLrv#q^s({o;I<J%5@7T?-d<>&0VktbVy@nw{4 zzMlTly!SJmyg$j_73Ii2-#>LVe6VGjUw>q)32^p>6_k8m8_TpzT|dkJN#}P><nvGT zKH)OgYD}`{B)QCT|DiswIN#&RpS`fKcA3k(Gh*%!9i4p6Y<|sp**RbHTk4)&VGGUp z52k&1vA6U4#`(UcINzDricg{wtKH^-4c`T|zvAqZkqvieZ?nwGH-0<!XNhH|pKNux zVrr$SUhw4XmCpReby+fHSFziidA{V<5{KW1CoYK|>Gb3J$$f{;Is2pg!9KrstYMoS z6E2@Qe$!<h3{Bnc7iyW?;XbD$-?B`>yOUS%A7-0v_4<E4s+DDihtE8=%K6>!<0}8S z>{!dxy_nJJrt|yt>ZF;GNe7(Yfd!k=lb!F<Lp}P<3-0Ue)pwrAezAeubY9apW#TT| z__h4%t|F(NN8;B-wyNniK3%s5{PuK}@!I^vq8CS6X2kZFW`;QX^Wn1o#gFZAzO(Ib zExp~<HuJiTnDu<9Z5}-v(`@fW_y64IclU+*(;JU#>HL0+dghPy9bD$?MFlB49RJ@w z*7VJN&i#~EPCnUpyOU?*RqpzyDvdqy&i6sip1(P|=1=L)p3nWf&FtSwT;{3#6+;4s zIQz8sm6zYX%QiKpd0l_tYbTFicg&dgwA;K;)#szV&hOERf9!kV_ycZ}7F#)Vz_XR+ zlL=Ed{?m!~-7j05UenZVcD7D@wsoM(m^SZ+HShGle;@RGZTXSMkGRdCIcH6HP?cF# zG%R$|yOrjbZh5;;Ty>fI0*B}RvD-HMU9n}_#_%nc+r;1Kc0SViUF7+9h41mek8008 K?0m>SOZ#t!cRUpU literal 32467 zcmb82XL#0S()JSwy(ILmfe=DMM<59`L+B-;CZPvN2%(cg2`#yE*8l+#MMV%46~(rI zqJWBGLs3>-R1^`c6kT*##0u~2v**vp<NqA{;pKzpF>}o|*OYVSf9I3&!$#B!3keD7 z9uZn2H7_eWFFR{M?#zN)YfQP8yX(T9e{MKDHE-(qjMol^*41(H!g=$uXBOlxT#!GY zLhxot-H6cYDbwd?XN@adG%x$sl&$}xu3s@Ubkwxa3MZb482g8=V$d>T6sHlP)kaOv zE6ANWFFXI%8gKo+{Db+=f8Y5+ooY9ZAM;fUZ8&n_f`aVi>3Q=Oj?A8(UznFYKYKxe z-mBVWV`xaox&L(*Qu|y;?FJ#0uly0sU;hQo{=fb;2&sFnZfIyoL}-;!3v=`H7cRK9 zMyG_oKc2h%@^ks)N_Ix1ywWAK)_~mng1p=rh5sqR9}hHW9~N@1Tkt;(+E)#^&@K3% z1|k26GDYuq5e^p)5=IJF3G;>J!lA-c;ap)G;aXt_VN>BOVKZTCA@478oKKM4T-aMU zMA%5!O}IqZLTJKj!WiLT;jO|N!p_1b!p6d?!eSxuW=JNE`I1)(x!*u?4`G}zTSy<M zb9d>fQ%lL2!iB>1!f;^~A@!;!d5SPnI8R7F`$#5_o|1W&_a;kUEaV<>c9kBP{OLnG zo%a<IC;g0;ew~o}AoZ=UbLvKa;iF&Yg+lVmkqrOo3-$hE8~pKa*?%B}_^1QC@#vhm zsLxvIR|qExsf$Z;L*W7;{)me@P)B6)h}Av#Qs(g6U-D?-G$FhlAvsY<eB=plh%-rc z@QpZnNk2xoQpmj`$?(iAIUrmk#15V`l|E9qT*x_nh1Xu4Q{SPI=_CGlf27V6guGAQ z^qqY1mnM7aO&;)>`W48IeOFmBeIc&44n6%2l|EZYziLaKD<m)Sp+5ALeM+44ZL02* z*Feeiqk?4gS(0Og<d-g)eaZX7q=%nlB@Ylz7ZML~QpZ6$Um|2bjF((j*i$%3$a~~P zAL&<~?C{6@_LiPGo-CRByvsRql<ZmyR|)AC=QE_I@8r+-r&#CoX`E#C7(5M=p1p^j zypnWIJmg0|R_lC@5dN^&=m&kLp3Hf%?vWSsM!o14al-HUy2m;1W=WqctRp0D=AOFZ zk9fJaLHF37xsp2w;mJJ7)Sr5=7vU-VAwG}p@%w=92=mQ3b8wsNCkV5J<dq~j)1fEd zhSI}J`p3M{Ki-d&J$3CU8NO2gH0jBUIY^S8{6Nmz=^Wnk{jDiI{U9#t$zHlecJM1x zGWwO0=@)a?S9<mV^Fdtffmqq`E_2mhdgAj+rqB50oi;ir4$kQtJY@gzT^pr)@QgiA zJ;<B;@T-sRG0*Hz`iDLBsVDmcA@Pwn_LFta_n&-vJMKkFb_)xI&U^UhJ?Fhex<@^z zXR7pbh4ksKe~U0P^qk&Oi>kWIw-JsrBkZCW*}=Ia$31#L?YNIW?$PVEx{rVKDbf>n z1IhSpEtzi;Il~R^wUAvmA)M(Uxt$RIRU{K{Cx?vR%F?s5nE~qCRp-=`-oh=oN&dtM zr{Z*<*(E;eNu9}y{E3_QnaNi2NBjdM6EAt<7dz@h-zw?ecwxM-jgXx{AMp=Y$p`=B zL;rc7``AUxKfEPg?$b|ZI!yMph14ZUay4OhA$`IReI*|Hgdg-XbicQdc!?9g`0FD( z>Oo(y<C{R;sVn@zj`v2$AKwVx?<YMx$dF83`L<9W_ysSx-&6PK5BI4P{A0H!$)5bM z$Bz7oi+&Lgb3t9$ZR{QD!aT7n=?DE`KO`s~ct<?cnS0dP`7Zv*FIoPE3FC!nLi~^i zzfqk0Ns&GA5HIl%N2u(m!$8M9{6|RNQb=9l2}nNh2ELFV=fpA1@jmtqrKcan%e(ND zxup;AfO@bmkg<dBnR=(Wka~f{%l=~zGY{lJJnRkfzz=(pe7Fbl4)L@9xZhv#)e_>L zIp`%lc_m4%A*3G(lBrt<$=oLn;(|9BvZK$`llX{dvh0vKCvWuBoA=-cb;XYF8GJww z*42CXgP-X6&W&=!MV*L`eS<ykat?A1AJBvB3-%oI%zmMN@TIciLT3IV9r2Pcbz;8R zYole)J$O&uc%ORny~Qv5qVLt@kA9OEb#nSA|BAX-L6{^=5W+X!r4R6f`1sw(UPTX* zH}xbQ{KH$~s-pMVpZF(Vc!QtnvZs&iLEhyaI9B#fzs|mLPCT5`CwR-g!$0+7|5Vj` z+~@mAJ>ebkJKrM@>Or50FVx{5UNUF&fw;IwT+aSF?+tO(i+hpM<Da@9!(;f%yMO&# zgty>gPx-7Q%n_~>^6g{h+dAwfNoMxSC9_-U5wRxe964VyTr7}WCX5#{gY150h@Hbb z<i~F+YRx<B_Kos8M+gVVXM*FL8q*u(EZNN$BCnS`RY-l|AhWwz=j25^^nu?p?63v0 zPZJIo(jR_XQGfcv`_zj#vUQ)CpuhB;d)P(GzK+l%oGGL)y(Ci?;^f=loH6)7TxGgX zT=dtOAG0_{_N#=v!#Qy=)5Kp{_Ff_Vqrc1tyAOXL=dE<VhcHPvSy){dCUgneX{C~h zhyHSJuFmNTzgegYd3BXtZy|jkALasHF?U6>Cm!aAeqh&KcJzz7=1RXpNZv~%^V@Wd z<UAqo&6XS?^a<Aq$sd34c&yIh$sozpZ=vKqLih}ih-0<RsUvx0NZ(9IznRw=($gpE zN`1J;P9}f+(?8;;Z}IZW+@wpU5BQxW{d6I9q<>wdCvWnlZ)G|sk4VYYgzOdS&`NrE zv{W+t%$YYb?}LMNFI&i-tSY%g$UCu;mkHqsbATOtoP7^3sWab$4f0c4Nc`+=>QYnZ zyic7cIquP4ctxLxCsX#f3CWxOESDZ$P=EMM-ABoe`691*(z9RKZ{$ax*zfFL`a^#7 zCm?_1kuJHfkUfh(=RDw-ImUlw-Jd3eUv(ve%o+ai<CYz9dL&aH_8587*J9Zv35mCc zWcWkB;XzZKQxEcDPPso*b_0YRg!GB|Vs0aKK1B#G<0TUZ^<v-g9fJQGWse<sfYgn9 z)nva~$ajo$@~fruB|`W_e(;(;bdz0*u%i%rc+^AsC?W3>2mKwWbNn+0$l*F)BZN2b zn*7-(^oe)y17c5o&}Ym4N?|V{=LwSO8+>L@<mlX4N8Y2Z6=lc1=RK1?Ovt++=h!!r z9XzWbxwYdyaV?cTN$9-C9P%FbsQYN$Bkvf=*iDiQKmPi+2sLOfn=V3T4i3Q`e)~{! zW`3FO;g1?HGwpOvto){d8{>7(40e*dP)JV9G`}U-HG^fBCFC1`|Jl;B19=DTGYbW> zqbDOJ=L_jgOUdOz@*O0Zofj>cUQxG1haH?FKjMY6^ryY<!CC4<J*flzqWAb8ru)>1 zT@Wul^~4{3oV$&Dnbl<7?<eFO9<dW{(>Zz3=V{VYuM)`<9GEVdc;UrB>7D!|5BeJ^ zJ7!^pWa`>bGI_&)^5Gop<*?@-`NA*aC7xBfN5A;3125>C$KjuP`lX*DB!21(ze{z_ zd=U?IBCck#t0IJN>m_sUmJCnX7fI3+AN5%$eX4M#koc%;hV;BkT=<)-bNWm_`7K6W z$*+>^*9z$i@$i0;&dHBmO<c2dP8|4Q{z7$*J-;Q1v#sNtIZTkAegzzM?3Xpt_YhJC zcmhA+a~s(ar?c<WmvfUndy6?=AwBiQZ$Ig&AAOiC{c<7uk3Eb%^?^UswNUrMgv3Rh zyaUhi&wG)&NB^lm_o;6U*->}m;yvQwe2DBbglWQjA$;K-c*NY04|76X({;b5aK4aw z;SZ#4OJvtyNZ)5k#;&hq@=K9SeBC9t6!Kl956lVmzz=rWx)&~F57U>8((`?SUsI*; zC8VG9m%8WZoOkIfbC#!b_9%5??~K$rd9%0RG4b@2U6v3(@SZ)sNauX7$g7(4*fFOC z(hn697w^xIp1kP`eWj17vLhbyrw{a}o$Tfa`OcFMeRb{+{6*>>Jf)BDtF_M2(~p4k z#7DoV>v)~Bcj!OAZz}3MS4f@Mo6da5mp#PyDp~j8<zUI=k3Ds+r}Mf({PG<^Zsj<S z{x5!$q^ExHf%@_;b3$C~qX^xnuhhG`<2~vRkC-E8y!eOb;qpTts0a7>j__TBht7NC z=e$RL@E}3&B?^}b8wiQ-uYZdWLwDJbBfaDsgAaNR=hAequdt^uLpWRrN8lDa9(#H` zK=u=aa57#pyOUZltJFMQcFYhtV@ICor^%if;2RMsJzPzaOuhppvpXXtlNY@tpIddF zBcz_xZK(8cmHX6%I>E)!vL|mi&+etK=*f4c?(rVy%ntRVz98R9{B)8Z5I(>U^1^?T z?Bj%dE8zB2=`)3q!YCo{j+4xL)Umhp^bwx$9=n?QCCDDW(|7udpBb{FudO6gcl;(w zKS@Yk#!EIK^MuUKAs+mYC;foW^n>~m7xf|jc)bTt(2tOQwvb&$9f_ZL!7q8U7x2$q z(+~PcT<GCN2fdpm#6NYMAbk%Z{D&vh2VPMp@&cIy<_UZ5^R1_j^n<zUqxb1Ec@Q7F zm-m=&_=cYM$RGdY4R6^e#LG7cKYSy3m-ndCFvW|%T*<^uUwb+7;+s5Idh#P~_E&qI zGbij1;+v&&?9&}`L&@;5yJYGgDVceMAG|}pQ)GvJxa15WeIsAa;cb1{!Bggfyoftf zcEn4+=wDl%!z=pByvOStzO%on6Mdro6J>vkkbL0-{b2s7GrS;9@*+<1X(hkRCG!r# z3*w@_^q0QyTLXP-`Gsf99lT@@5(n=PC;7sM#qvjf%q8EmK{|);v67h|`avJyb2Hg3 z62eFN=Jemv;g|2#80qO3{OTe-c@ii4tGUjp8~2>`U@yZ<`j0=p1LWITesdh}CrGAF z>`C?>a|CZ@%AWqkOXgkXriS$Mh0IOTf3e38JmB}zG}%W9iI4vEm7aRiKjyfd&Z$Ei z$zz4^f<AGd^IX}{cl^?4?lXs+XX+kv1OJ&z<_?)Yum{MSdFdyAoa2x8;0f;$7j>jg z*y9g-{NW#-g4~;|IEb6MbKaxQ&U@@f;vqlw1^S+PXQJ@0e~Z8su#Rl9gd2qTruW!I z=sZ$bQ8-mtK{!uHZHNVq!xaz?Wb57lAv0YlIY~&oIg;VnVuu_mnO(^)LVt_STMFUs zP|3_nhGgQHE}2<n=Z%)0-zdbxJLJ<#cJzT=MBVA(eA&Tq>bhQf>QG&BQ(;vh9A)Rg zA?&-!jy}Xlt|Ckl@=fDg$h*{|M0WV6->s!5U;K}h-XmNn%o8$`)VZ$o++QM@c-S@A zRnqx7A^yuGlRt5>v$?lXcJw1(GW_f=d6kg*(SLRnd~YkeaA89s{i4rvrOy<?zXHi| zLgHlRL#2-tt`s`s<-EPl;cp$u)R*_*5Ao8MVX|K?r0(@3Ge7X1JUOQyePj<`+DUFH zEE3idvU{0Fm-NgB_2M4$-dlF`IZQJ75D)u|yvVzi?3q8mWabc_aNbGh#K$}(N}nnu zUUv9e>EZoA$(cgpf?w1Pzwnm&U>~jf;X=M~@Sk_fb<RF#PU+J!onyy*HI|-vkCdDw zgh%WV>dgN0%I-EHb2LIS^&oEc{dk>|FZU}*PhH6)NqV=Cxf~~%`oS~m#=F$TC;JE? zc9kS235lQjuqOuV93IfeZqma~;%g#(XJMX@yrLw-*Br_0scDj_cU8&EDLjU^3v^B$ zdq^HAga;{-X9?*ub5H)<W4@883w+``$#;|aBp>`U$Jvg!@+Gs+21q8}F_P)qJjqRk zU4+ybq^`t2S9b7{e!&y!NS|xSev**5*oV{u9#TK*!2Tj0>RL(uItX)xnL_d!B$>Uz zcciZL><Rk8-ef;e7xoNsWa<7CA@9JmTGCS=_|12VeL%j{13u%QxRT_L`?V$W&IHNy ziMV+e8NbXq{@@?yKKWz*(1%GsR#;O=ABms)#7`e@l|6i6eq7S?9(`v24$wJq6gupg zQ~E?7c!zr%bT3j!zn%B$=-hdaK8@5lcJK+m6?L8|{OjK$<U-70vhfSUg_VWuEM|DE z^voD`<lIQ-><s*;Nk3XxMc78jPU9Xsf*FpMU7oOokQ!2pQt9E|I?1bra35}<XJ+9< zs_e^z3xseV&cYFT*j;u3VTN#ykZ%!oL#2-tP7v}9ST31&=SpU$5Py5eJ?iC=zPb=@ zFe6Q*r*G5+|M;cf)So)?jU%p_^5+&N2y=z(Hu_vFeT1;SaE*{}0DY(=eP7{vA^e3G zMUMO9BojCFD3E@lkUWbeGh6VS{!P_6ybYC1{`4C=@}kc8ZLNFERh;A~A-p7Sc#fW3 z&_VX>vW}AB6ZOWPUARDYw+X2?@scn6V?GLHPaN=x`#Cy?|L{3odfv^F+(WodNL_d* zUi$Gu^6Dvhk+8dvZ{;Y-_$A)i(r*+}Kl-vz`W3=72Xb${^vi_gzeF-SzgRN+kG#4` zPrau|Ccg%f=_htArRN>`#6E4W^Kv15<39Dro;|>R;9cegKh$Tk{IcgJN%jhvuPT!1 zJ9`QL^nq`BFWF}b$=ffPy;3BZZ=y+NFEB4-q{sha$*mknKC7fhCSU4C9pMLg&(yuH z!eK)C4Nu|aBAs&&nS6PVI(3sh_t^vFRjTuB$GiOYA<pqSrw+{3IO&sx?DY|n8wu$T zbDJqW@*v6NOJ3{&=8x|JbB<k%?&k~HpX7}n>P$ZLnfT!y{@IK0Z=U>b6!s7%3YiD^ z01w!kg|cJMSCd>-2oK?9ZRwd;>Pf!jMPJ!h1-h3aWM1J-xb*O)MDlXsEy86&?orQB z>EThNWO&D%P-p7OIdj=v_sBCxa${kdkh$UaKwIet2rCPz5B{U1XRa1VrZ3Dp{b;W9 z2}1H4E*Zb@w4L<Si*vqP&bZQLKUfG~*}L?e`cgOMql@k}6M7vo@5M?_J>V(v(MRkX z%AW5$eSsg`gE!bm{Fi%S(mU_bXOQ`Wx73BWi2JX9i!d|v1dej2tKOzY@sc^G_ryq# z;n)V*vjYN>*(vmJw!<%R!Kqn}^NEs&2<avE)RdkRFL930J$^gnNv2oSf!$qI=WT_= zUtKcs50y->`7Jy{dg?%3YB~H>l<XGH6OspU#z;^7mP#gnW{<e(6aMi>T<lch#m^G? zWrpYjd|`L@ksbBr9)5Ti9;M2jxH2UZ_jJkB2|wh^yVNsI_RJ#u=_)-tg!o2F&+fs0 z8|kx!<XtA2cpFKkUi2N_kuUd|ZQj8j-w5~upXeLAo45z+UG9P8H{WsIPjWq>Ur3#2 zN+xgm#yRm9Is7o^i>0UD10+uo_7_ePCJA{D|9s1+Ut8JL6!J}2Avs0p6Pg1#?;$;X zDwdoqr2lm!Hx-T%!mq)S;cu>F;-~%r>FF<gB##80Q!nCRpAaYer&#u_gz#jo<ds6= zr%t?+q;vR6J?I~C!87J6U-#TX`d3+U4I%Y_@9=?d(q!53ezauf8hxzv>^1H&FYto* zBV|wid~2Cc@}V#AlDO#GMEO}F<UNo)*(30yqwGQEYqa$3gycs&^r@T9<Am*n>}B$$ zu9-ULejkTD{UNUmozrLLD@=O6JM8@g=@W(2pM69g%voRA!Q<N`)1NHK<lS2`{G`6{ z5kKs`rLyO}E)M@*$w@-`gdg%;rt?x^jxb$FU5JN!Q94JauH;95(bGTt!INnDDG{<q z!X-xv$$PD2{5FxCBBT$*1;QJX-D)BCnYT{Tv!B@~yqBeOXFXa<&m3?bBmE-By@2Fc zp<BovYa@A;Fjq+Z3MErN=A8Pwbk2NW-$8of<URU7&T(&oWc)IZGo^<ooWoc49(|_H z^k=c|6E|}ZFFp5|uThTo;4|L?-h(&9&D;|oeZ-!9!5$<Z_B{8f7xA;#ocEX~=RNq) zU3qZM+`?1#5q$aU-y&^v7p~E>B<VW~>8V$8nsAM9y^tEiMQUa`hwIfOPZN$2))o>M zGehmD6*JAf1-jQu*k4E-MUsaKlO6I1$#aF%h0HFqlp#HNf}FE!(Zgxt;k>T=xP)+- zdSJ(nTPizxPJHYfW^je<;)TP7_)Cz?w}U$K4t1fP%VpnOI8WG9NMER1l=Sq4elU}L zb<Q`iTym<A*;^>Nw-9-vWOg8N!pjXhA0nJB%ooz<PLlBt|KU44i<KR5F>8~gr;o*w z@jF0rGvR0lQU~~mzjVjF){>_Rd5=Dj4}8m(U8-=9u&qN6FA}Baexzjh3V-3>LY>oZ z;$!#b>YUw+p1gcIXCJJPJVQ83NPe92t;y7RA0hb_N#@(de9&L`*Hd=z0l$9f*_G_K z`O@RRR5HA0|Ivr`I)_JmyU?@CnMZgMse9}--otLK&e;#}r>peTo83=esncZHO%c*Z z`aqxIIp?copCzR4e6z`)y0XXkjnhx}m@oR<PWo_Rv5@z%#~*RCSEwWQZFGO6kaxHb z?|SKck#MRIem0N{FPJCz*ih&A8zq^2>XnS1I>1Bt#Jkk7p6>BmX}M(LX(XBan7a<r zQ%~~Bk)Ap+hfN)RQYEv8*r&`9yrK@o4?ppTf9wXzuSZB9;2HI$KHQ^k>>u8z{_u^x zP$2)*2mf=V&lAE!z6-;oA0(vSybIrXpZF#@?r}~Xc$a$87vk@x`-wv49Uk+&qL0xI zJN8K%>DlYV&-nzMBd1A@ap;#w9wB7E#Y@I-tmKYD;_M?C9#LoJp`p&Z3-JTbxZh0Y z<-#UH&e?ZK(o>gPBy*l9nLX~bYp?TVLiTTxWcDHP(@)NO%C41gy>PUUdasu37gA65 zk29|7vSYs{Nv6NN*Iat|#(UVaFR*7%@g92O;$CO@jT06No%i5JlFqqDTpsDugiD3g ziMWv)gw+0@{~P0zKb)-6t%}PQM|?CZc~FV<bItg^ytm(`{c~cQ*>0D;8WUc;rqpZi zy;%FF1{VTWSQ1ui(r+GXdU9c!`zf!jT)cHy*P*3$t$b`*&EtMczT|!_dYa!Z)Tq~Z zMSY*W^759jzwdTg?P_JqTJ7~%;!pkeytgZ0Hyh`l`NsjT{gAr8Tjj>9Ec{%Ds5|!g z?4y#a*^MT;?P~HbX_t<9Y<bZ`GkbRRS{u*YqvJNYtxdl{<2U^pu&47b42($k+R(M% z=Im+gweI_eTt4}X%YG>z`}o5hz1A*Zlf7sCc3|6;Bab!n+l*Bm$NKkrEqcL&SDsqF z#y-5SPOsy!ZhP$hd(XEU;<gJdS6zPiYoC=rf4T3FuLAaMt4R%C3G-XmjXT!gxIJJ~ zu4g?xWQgBBxpVJ>PrXoV&+mC+bd_?CJ@G`wz=TQxJC*irlRlANYkzr0Nv{JQJ9GQw zf4_65*LLqeu)a-2mp!^7Z}D%d{MPgJGfzL>uFNJ)^Zh()QJKZ>dGV1AbKLgM(N})T zN$}d_FQU&r`;gmyzPV=0`Oo~8F@NCVo~!+q*R0mzTCqNRtbOvEp?PJt`i{oa_gyZw z^`4KT&({vv`_o%)J+e1oJL^1jZ@qg;Ejjnf%JjM=_SCj3&xEJ?EFtIB&z64Vx4Xs` zoNv&?Z~4vE*W1<9XYo<P9y?k!V402LJ9pj|uu9hp2X8nS%>Qyyh0Sl3+SwX+RJgIk zV^I^cidS?lvoF`}uDx}I*D{8-m=NtNvD~{y+<GRM$M=;#8@-~)Z~e~~t&hDkU~5<Z zx~=2AB^Es~E9U7t-PY?+%tKvf7hA!Ho0eQHFR`6*jh~qMgwI~T)^<w!=^k76`QW*= zPI+uX%#H(tulX!?``KX^-wxQ2dgs@s_VwGin8<30RsELKeC|uD&zD-t-N(arzUsGz z<%K)S5<T{QuXp<$Ug@*A4tvY%&-GiqX4`jO^p@JnqUhocS4ypC;N^RoepqT9wvKP` zQWLL@X+HOvEuGvJ`(vXWoofWFB-`(bjr806S+klSk0`Y+0|zg=F(+UrJ;&bJ>hW9C zr~h^Mhd8%Qj(vM_gYBhu^OY~pogH6njqlrcId7!j&ZO3Q_w<Zn%dS~rcF`iQO*{I2 z(huonc5%+0e)F&UY}muqo4h@%*nXV&%GKHF9;=;^_|G@*_Sx^Zj!C%xnA_@H8#Zd- zFFw2D$|Ji^K3W>gf9i&YfAibpYbWhHdNg2t3M*uc`Dl&xiY~lVXNS)`|9blH*W0`{ z?rQnrZOdF1c(TX8F2$EwhXoJrE8Om}tSOVLJTo?!?}#>KtvmVbrQ6Tccy_VR4prK8 zPx{_cyZg3Q(Z}BN*p!c(y>=qmXKxHY@L7wNe!H)@%_mdqdF;{V9}gY#_hLJ8d_=Q8 z4Fk5MS*^xZzVO;BclWCO?=qh~yY0%nzM;W7e0MqOi-|rvk@oHExz)V(z{c=Yd(2}k zPQUwOd=;-%nE%18HHVd1#X-@pz59~Kc7J-ui8{eJhu**ERF|b*Th%yhTi`3dJ=H$> zuEWdy)^0)6nde8jt>9eqw8U{G_U_>GMQ;stS=opN^^PS6d3*iK(yMQ}?cwEvukEZ~ zX0PsiEhg&&k3E=@^yTIuE*r6-|B~Gq0lT!WZH3sgZW}Z>FXz;EUbCE6uRXgmVCA1& zNPjG-Z@=`3ukIREYNu~>^^Dlzv0a-Eo(he1*$My2UISh)v5)Jw7=B~-8vEqzxtAk$ zm)iMZ6P_#F<+q`+Z{K_Ne!soAa9X4A3V!=OVc^BQ`W_2^rp=?D=KC%CwHtRXOAOBI zn$x@H7nfSpqjz+E<1N3fKK(+@wCoZ)bZz(SAK&p>WMPHX)uRITX{SR?dew2+C*$mn zQ9=H`l=xukkzd^Q*2`OV9th_1_L9D9AFNbnjpw{m_DgWiV;4kSEF0ywgV%Qd6u!%6 zvyXpTDYMLPOFFci-~YVFT-)FIz3Si+t3GAEEAS7m^^SWq@WoZXm7Sd0cXV)`x9&>$ zYU6g7{cYQoFVF2Rwg;k~-`ZwunN@%Ho&~>l^_#2V>F-Wm56;_~(^Kvb?t_D$^{wz- zl*fL4{^KPh<4Vl^ucz<Zmg}{>o4?v~`eZQQ68BfhF#(I-QSI<;%gbzHv-P9iIv=o% z-x}|`S?IB$uTLH|xn{BTxcvF5qdfwal~k{Ba%PEDxo_V$5A7(i&YRwTWX)3^`?}w* zCJzSZrCF_WzrFgX$G$Axo>%eXfbBYQs>2UYyX?ELZuJhe@mk$C4u+f!DYfW}2VVXC z$AHzE{XlWb0>3S6y(PO^@H_Ot{sYh4Y*1{KrrlkA(-DuA9T^t#{ZGO9*zuEVTZrF| zT)O(w$~pn7aIV7(Pd4=0Ya6<J@Iz#YjeEImQRAPz_Exi%GX}orwr|FManFHZ{dV7T zPwv^RK1-?oVdq9ZpFI+H_li5Od(7J2xHNrPnLV|>@^`DA_1GuXE<|ViRAOJ=cYDH_ zOMbiMn{^MJit*Z$FU=cv*LkmnTrX_&O?t6aU2`e@pK*Ts?t!{-_XqjEywmZQKfE{i zJ-o-?VgL9td*jC5b(L59?B(Of^PVsES=+WZ`m~xA+)sPg{dV&kuPyGe=dq8!^IMgX zVbf1;c3Hog6{^%r_S)ZG`gX(q;GCb|-h1uwpF9?`<!0Q49mV!U&b|-A3;dQjWpdz$ zO~LQzzrU^7CHOu5Ds}y`alt;f9r5sUqo#U;-`)JN=N5ad#knWaPu29;>DHBJU4GSP zVTUVKuh=g*=a$-G$2hkowExGoz0FGO?SCGdx35)jZ%ldnuD}y+s}Or;%+(hH_Un#s z(zezrv8#3>+0)NwT@J;p9n>!PJD}UbmK~3JENA$epG=$Lwyd2W41O=fXP+PFG`sU` zuYLV_kL&Xe7u(j)Ui|p+;P=Cq|N6zn!8vWdquRS`R(fsgp}Y21T~unF_qGUsq)KqE zHb(B;GR9}M>uye+`=;L_4!oBYw#8+AyPtifNtVxc4}3Pc!jq-;!I4CFbhzJ693OFI zc#AUIAN73QUncqOT5{1(`@Rj>E8pI7@M>C#-Aw-F^vK|M_4`L&-_yBmsXg9&Vv~O6 zv(+VGPmSH?vZ{|aKT#_fZ%X>4cg7zJz7sm~ktV_YJfP#vPrCeGYOfX4`>yawpRJmg z)8NK&w;j4G;-f?5K_2}YyM1_q+ur~C1NCod>9JO&M@}cWz4q0#AD6D|?za1uHvdQ6 zT(|YCp8n{zG_Mu>9RA6Xbf5L8^K3+qt!^9sa9O*okl=oM)b-gN3xajOboIwAnI1d( z@Xv4DbG^jgn7=>A^{d;)zO^^{Xpi8WA2>4XZ16kMb?U+Ie;nbm#7%E^J-WHnPJGob zxl&QU-st+OfBzPr+12uK{o_0qIC=K4Z=lQ0Up|r%^Sjq#AK2f$d8pSe&OCl>N|RFC zSpMaEg*Qqp``+sXFCFq(lXb%mE`2dzSL$!;>E7+L!|(5LjX&eJe>WVG^sgU_?V*p$ zLJC*9tj_kSX;~izEcDc}W8=Q^SpSYUQ!Y&^u~ioz-j}c>nD^=x{YD0P+2_LP9<#@k z+LnY5My;LTw<$Y6h})A?W>NM-%}sk;mQ>UAWS^&eHr%~_?$5UbzpMIPn%rPrcVF+| z{N%s>;MUV6R_F7K9c4#ccJGx(&W_&fwHdL)cV;&9*h9Hfhn&c8+u!Pqx{>;(&-$*N zw=bow%Q~(+zq{lKzug{p_vTrLJb(Ia*`xQx+^5!9-G+1CACy>Pk8L`bTr>E+=G_+C zr&jy*{r~+voYtrM=Fn4qd!g@#2NQFAf9&bJ*}*ytuac7c_2(Yj{A1<Uil>&^(C%lS zJ3Pr{16Dm*{cLdFW<FSw6EQisC&JQNw0D=-^_w^Q*1O`iu#CiB2NK;jVE>RAD=L;+ z{Ue|BxzXGG=kM=1i<|l0treV?W63Z6R^+m{gq)c80j@ve3-a^sE-(51w7>s_m)m}u z9kBc38n4ZdcUkhL;p?9Bm)NIwcHO-Cd-tFB_?;j8o-fEK@>$55IkQjpbpM}w|LyOC zU7yzb_p`zIfB4EH{fdJ7kiQQuzF7Y5^2V<J^Y=mN()%ht6yvwr4^)5s%(KA{$7i37 F{D0KrmInX; diff --git a/IntervalsModel/simulation_data/ws.dat b/IntervalsModel/simulation_data/ws.dat index 687a6b2c3245be85104c1a36ed856e2dcbce8fbe..ee07dd956d8f31c12c7c6005b403b532cdc55b5b 100644 GIT binary patch literal 31253 zcmb82cYN3NoyQYEMu2a^2(k$yVedT?APgZ1ksT070s#^rK-eIAgA6_GUC%pDJv}{Z zr?cwS)^QK**(#6PR_j2m?V+O5b6inQ@B0nUC-V5VzXz#*^u_!AdGGc8<*yE%lGnIS zQc}|3wB&}<s+N>hl`Y9FU%b9U!v&v}-*x6YADwt&Vb#L7=RW&ba#NniuU@&bZ1MW? z)s^dVYlWIgP1BN73QAU$Et$Dt&C0S41-DkSo7YZGE?ksc>&1uCW<<Q&p=4U2uW88* z3QMZimoHvfwys0Npa0*M->mxaxqfGwH2CU;=MOeWZjryba(&tOlB$)f^UF%sZKx_+ zRaUtkfA#yHOioIA>vCgtrk9nhjBLM^6r?BB{WuEi^8Z5lzv;j9q^55*O-@cqORiVA zx_sTb)s-C@_8s=p{pDLecy!&&-FKuFJk>wBac=p#^;PAi8!px&a**DuPSRV0!v9I{ zRX^#>pzwdvlP-qa9{+>E>EPAiNiYNK1MURN!13TPuno8e><R7#`+@pt2d%xK(1*Zf zU@veTxD{*(js~^AAG$5L3v2-{0kwAmdJtF+W`f4q8+t4_2{evY&?`XmKL+i)@@wI5 z0H=XVLF=G=5&St|SFkg96tuqT&491HHPFW247v%}0kl5Lp;JNoWc>PF&Aahf|84Mh zfcA3)^a8L8m<<}o5a_;OW6*loH{-I7)-@hS0s1B2!Wd>i8;A2HJ&N}eU@2%{th@U) zocHyh{gIy!e<WzUJ)rx8wL#-|zT8jys{fA2Uk@s89>!t)HY2w^Ccg=~JD3CNCmp&Q zXkOZL&bIJ=5L^YCe=F#Tp!4V6S|9V*kMlGXJ@?i5(Vlg5&g5&q3Vr7>1KPN}Qv=~! zZ|yq&IlM0hXMpp;6`=F)zBuQ`XCC&^{uiP@8uU)pgSM~6X<atP{2IS~DvIg1hBkid zWIYD(uAKLvH~cZ6b7x!);2#CGe;u^*WPjZ=>vJt~xuE`)KM8*}XdT?ge(<fI{g&@s z=OL%vHqh?B_Uxy6)<s%0PQSm}Gma*)cmrtdd#|)-Uhajo*a`g`K)>_qkB9GE$~UiJ zF+KZkAGGgW6e6Di8i(;)U(vZYzPacbhj+qyYR~#P2XT4tr}4SB=I<P8-+b-I2I3fz zARhZ@J&KX50Ik~|Xzy!J=viO}=-gWe@2qvwj`5#F&pk8#UhtK(F3S7vT~yC|Q;5EL z#-m<Z%<oQU=f}9?cH;ii(aQnnf%eI`w0EU^jMutfNv|LNQox$}bMBj==l4$iL~$Bl z++PlM66sm@uE_Zvu<r4CX~+4F*Ij+{uc^HL+`D+baeo<Ci1+fpMf`RhMaLWH&DVBa z-lu?bz^UL&Z~$npyscNkSAPO@5ok@0L+=8uk>3F0@LRkVIpd!OoeX-TtX(?%ZJ;-* z7<wJ(w{0<W6{x<qt0{c*s{?I5{w7%l-<^*C=5lw;KMi^Hv~NAVefyCsi`m}+oekPg z>yQWE-^A9dIehD@9s6znY9VL;y~)Pw?r7H=<L(;2^&5&kzqym3mACI(;9mn;H}mWZ z-*2Mxvl;$4a0BSAG9T@#zYw_+(Ehs@=G&8Z`<)7{9pl%pw71E>4efJB^y7T{VqEq` z+TTW%=o^RjobN2&t-tx3mve7?_F)Ekok07mTp@h-z_=#EZv>j>Z0Nx;y>`%|b3PaT z8ZZa+4mh9v;SU1~K=X56{B}0x-EX>nXTi6R#=Q!@^P-%4kixrnQU4R*SAn}h`(gh^ z!cPay+d5cJ?^qk;j8nO8@U`dO8J~3@f!q?XJ!l=RTLpaYnDaLtzI$t(^y}Pu*A^ge z-JKufIn29x6hS*zW1;=+HV}FP=seg*?T_YtAJ_qOkGnt@gWfUkj(cm~osn~oq^*y8 zX#Je`<LFHQofG%d{Ec@Oa`~WjI0mhp{gu{^^W(lYK~Md?(9VzFQ|D1R^-dyhf4!^M z!dFiJZQ(oT%6s3er~B(3>DT)8$BzD;Tm8l5Gm-Z@FdW)@7MGVEAB(pSwC~2F9p~M? zJI4divmf>?2flsrKG;9)rz3X+)Sme&=X~2Y_fWZb-0rvaQ7=(F?Voz?f%~_K_|}5< z(fIZ6d}+tM8HV0Y(D;o@dx`RKUQ*C^U#CJR(lb8esHs1{ANJk(P(P7B`|lnqUw(x+ zt(W@NV@`s4>Hh@p&WU!^*MCjr^%tL;MEcH&{?wmJT$leXqU=20{7teJdKzf27C_Gj z4}#mkdZ4r9j`<t6FYkN6^`JF&<{QDcR=c3B<xptldqR6_>O;Gu##six2k4EE9tGc7 zNr};ZyR1hM??b@Gptt=vbZyYOdh1HzTW9lG4u28ow{}0YJJ224+h<*kb2aa~K;_NP z{+Rz#<W_;#fek>v0p1Sn`%T%2oHyFI`@;`F^I8|<XG2$j_S+ko2jAN>3)(p`9`kdq zq_yWwl;0aW()P6keml@UI#(0n=Y#gwJ+PnZY2Vx9Ui3jf6*Rsn(9WUzIt{+Gb{oQX zK9)hZiOE}+$?)B$CD5(F>EKv!Ht26kZ~032_N585b+m8Bw}toNp!4p0IB&*je1*tY zfjOZ2q<`nid!{|-&N%9$zYlbOjmNs$5Bu$%=!~B8ZJ({9`qv=W9qb1hM@MMo<(I&> z&kLb<gVLR#uLie*_T4<Kr*$t!PXFqAXD0BT1I_@wTmFXj9vc5h<h<j~wRtq--M#bs z;oKUxb@bjGK+isnfVMyOe;)iJpmSh6*4KM5A35(`S7`0oU+d^Ty1&}BE*0o&$9TO9 z-U;K@o^fwNU%Tp^fUmrJpq%kaJ0I#hPu4?y`(S?Zt&`t<_55zTSA&SlKHEq8C4VDw z-VN*R+#Tm#yT(xf-}u~P@34Cpm$&Z5y%;-*<gLH;D?{FV64|%D+TR|_!}#4lzYEvI z;%N#!6>J5X?^bB@nFej1#&3Sg8MpG*Upe*M2lrULJp8C<zvjSq9tJ|Y=gMivytLzd zO+arGxEC}o?Ip^`J?ii`>e&bNYwFMXxc|ndexf+t1M^qjc~7LLyz?-VJiS}-dTGZy z6YsnF%dt~adHp5IH}21Q)}Q(Z@q78-B1_?U1C(t8-`umHt$il6vz!6#PWr92XQI1k zPo(o>_6|Wi%gSqiI`5l7cYPbQJ+@{i;ZFjc!2{5K^Ze#mD{tK(<jg+}dKK6Mw0?`B z-7)iV_pGZs>P<Qtvn%b5XdScbH^M$jd%LWkesa(^uOra*!TOr7`Dte`^6rLyt!okQ z_S+jd4!%2YKE~xPG(m19=x-M5JqEtNiL#*0$2k3V^x=I6Xy2W8zirYb$PECSf}KI* zbPnpmH!u66zVR<Z&N!8`j{2F7oZm8UOgVh>(5|<4H1C_h-JpGNFPsm*wUx;EZatl! zS-d;1#xWUwJJ5Kn!xZ>^!BWt=jE7zi>URyadt*PH3;S;!dm?YYS3sMWa`q{Mck85_ zevRL6uXn(`&|gK&j(g%B>31S>ZNMR5FVOt0V|V!1gX_SC;I*LlWd^i!<UG%TZ~yl| zJ9qZgyXAd#PNyL+KNH&cCqdh%&Cupw2)!NjzS^e@`0YXCSq*&*w0`zwDSYq4PH6Mi zzx6i1<;WT5Jm}h>_U-Rf_!*%4p`Qluy?4^ixBKFpdw1-k{&KP7owPq$@Qv3!?*`v_ z)sJ#<+Br%?&%EbDD}Myqd2;Wz!1rA}<J`#mFwi{bLfZ%XBH#GkFXI@E-gTgS<MI9& zk9Vi&3i8fxcFewcX#W83#-ab>1opiP?y-F}-@53zhu*PG@CSqX_wE{(v~_a-?St{@ z&+phm?8fDkcMq+ndKKuYue_+<3gkM20caiU^APyzTOaLsH?6mK%)K-&?Ip^`I(t9- zPRz!CB0az7+mNfNzY)k+f>km7ME;xu>#V;v=qJ*%UK@~W31))vdTB?y)?Yj7pNRRr zlDz(wqG$Yx{AD4hKlP3G^1ns4AZyR2LHk?f0JJx>2lQUhUN(jHHl#pjfY#6(zW~19 z2zPNDe0NW~%iue!gQ2bQc4%|7=cVwq<8I`@_gknxcio<uzuyAiy%l}2vkml{G!yy| zs6FfJO|<^j*Lr$eyzRzk?=!G#o=2cZf_1@RU=FDNTxjd+?hb`-d_~a8w}Z|HCxPx} zA87R^LJtS80Vji5pmSk<`tf}!a@y?(T?v}M_Upkf0PT-?7?<-r5IOtWDn`$NE&-i~ z0D3p*w|o(_-(2UsKK${Zd*uAKfp7f}Lpu+~=?!1ayZKrlf9o#h-8nHn_iQfj&Zlu~ zg+Cgc2O8fF=*^&ZilME~0ciWDzWq{f6LQYMPH6kv3EI7H4!r`jUiK#ozIC=f2jMRO zXT{LpKwID&=N@SJ?&%Wv#%;g)!@mLS2zraPr<`%Q7p2HMueG7AxBOx7jl;eygFgzi zPQ}pP0q4Q}ui)K&TYu%9<8tJr-4E}xb5n-g1kn0<*Q}d+yAHWYp!xbevd;^7KLNT& z?V+9L2GHIs``#D6b@Th+eYaocJsf%Gv<|d)-Z;C#PXQZ)`#|+lp}T_hK<BCpwDFpc zdu%@1SI+s9o)WWvHMDuDKN`OKAm2LY#`LVK`(nT9Bj??<j!oeA1f2`>FhA?!To}h8 z^xO;eoHzZ&<?Tye%)k0^dFRl0HlSyoae4dX?+WkM1oWDL8KCz3&2ByHtM_U@^2R?G z+B(}m@2hq5cbs#koN+5J?L4dJ{5lWH?ZCgZcgs8LJu|+#$d||Tl+#W>-t|`qZT!Zi zy+rvq2ij>Dvy({A{`5w^rv8k7AbRGdej<PFf%Q|q1$GkY>8~$x13>Ryyk6SzZpQnr z{(_j_E6MB6J8t}m{54DvZ!P$j|1DAluM_BRG5PBB;N6@DLwl3N>tk~A)l1{OHt23< z#`K&$ceN?+-qrxxIEta2No#LB&Ya&^f6I7l{e~UHj<?G=N5XfOa-dT|Z?5?j!(R(} zYy3tTkGIQjU1?0-n|LzD_nUGw#_tDh{Q7lgq*ozle17BC!Ji14xAD3&&a(BiU&>hz z>pBX1?%-kQA)q^Ezs$opmm}xA+yHGn`pJUt+~q?L0*`~vwQ=u;e;sJux5nh1W9LHJ zxig;C=#_)ky=Bb4H^z5wX$Ry)Z@hJ1$@^l^-$Z+${gyeOx$w2`{9gm#x-WvB1^V6y z+J23Jw%^X9d+OYHPpr57v){&LJ+)(<MCaW6?N1ltu#VnHQF}{~b6&SY7lHP%652fN zr*X{Y-TLas{)xSi+XgyU1EKR`eCN|WG`^0=#ph@;{CM2n#m>C10C#|Gz!GpL=$_fH zvGARz0%+&lxOT_n)<Dk%eP00W{!2SY-rYjv>}xA%-~CPMyw>O4x${mOgKyo9TfJ_) zn}>C<pY?dR9?H2l?w@(;-}trT{5hY_zjwC;KUJV}z69EM?cXf;&ZT<FdC#QBA@3Yq z2W?*J4~Cxs%GZu_tKLTBc7gWGd2a{by)#eyzmNCcpmiyR){k-2fxi{hPARl=5|>}c zyLr2p?t$;xH=ca-oG<CPeeXaz^4>M?lkqm--8-ip=g{wk@5U3C)2{XLuDVa1@KXdT zryb>ew~y*gMz1TVo^z$1cDo?w{Yr+O2ab!`F@EiMSN%>or}{H4?djLJ%x5NkoV!)f z^T9-VzPCcIrvALg?qNktKaoH8A?|PD72>pB*4O!PZ{zh!<X=1LXVlQI^7^y>)+>>} zEB)PY`QIYm)C0tH0`!JvLOU~~pmSq#?zXqf9aC>P@_wr}LmvU>f&D<geb(MxxA$9- z^G2xeEF0%=<gArF)6QJprN_qPji)#K>9Kde`Nv~++zoen1Mk+s-y-ATTR(4--`;V& zw*{Tqk<fj>5#S!sJ}Bob*%x<b7V^gBY+3*Hybl1?n-SCV#yCgQc-M~oGw;E?7lZcA zZ?t(Czx{KE?VE9L#ZC%nzuH6RfPF#hW1sv@<Xk#;+R^W5^h>}}a24nrwS+c4^J)Tr zFX(R<>(mUs-%K$BzH@Ed`f0?wH^TijF6Y==ZT<Dv0sVZ?xh#WT2x?baIp@+o8`n|v zoL}|Le?RZu0q67({4CJAFMu|F<I~?#-nC<&E8#nz&b9q=-`qccn_4&belvFb8_nD9 z{#p;`V+r!kvvE789eFn|<8WS#)4g+koonq*#f~@^+B$E6o(ygSt+VxW9((fcJQ<(y zdvCOFyc^L|PJQ=t9`7f?HK6?|fbIc~1<lvrj{X+7j(6|m3h36L`{#VA?;Luk-0uPC z4FrAf9JA+~ZHMnZ7>D&VKI`CoId@~wUkMtQ^JKivMJ95GLG38#zDqmDi;;JJPsH?x zKp&6sjazxKIdZ<+Kl9LzeXt(Z)xAB29sAS=dS5J#3~2oshjXkQ>!G~!unYaUpz_vR zdF{Gy&fgC7;`Y4{`Mh@l8-v#0{8zyr0eVNAXYZl&=H7Y7hN5Sj&VzOKz72@U8MpG% z)<-?}E?$54Lp|?^`8xO3N&n`fUKVjlPl9$X%v*oPr9Jaj-+a=sGZ7pKUP(_ovyrc< zzii~Eg6_ZiiQ?>yd|u37JLp7u#%Vse$cyoMX~+4tPTEm#Z_MwN<kfS}<N3z*oqP2} z<u3nQq&~4M0c(TKhCMcBYdHkDjbL3c8T9t&Lwmcd*(~_;Kx^!74~IVpbk?o0zv=q( z?k&3t+MRIM&BNKBj+}AIuLnN`bT^klcLTjq&g2&O#;f0s@SU~p&_$s2t_5v<toIoB z*3}tSekkwW*Z|sFcM#fq%+va$@?HpfV{)LK5BpUB-#8CJmw?()uNVAc&^WGzb}qU= zHwFEcHiPz?-w4`owex2^r}6Gw42Et1>Tf-?b(MCGt)Klrf_xTeetvW9_YB@=g6@_2 z_G1z6&h=#IhM;}D0opj-;XUvxz&SBF<rl;ETeT0mEogn6GxyuNl_O_-+O^KJd0z?o z4RAmG7JDzOr}^7Y>(>=K_S3l@27fE)H*`F7G3frAXDNLBO@wxTtba%N($2YaX+NDe z>*fABkN$SBZ%47KU+d*vY|Oj)&V#nj&haw%*3UhcPT^hu9iZJa<1~-HF+c8?ajI`! za*?-g*4h3U?-Jyk7v-hR-#g*^SoEBK{nvr-{yMk%JH~q{xD@og8?=7pPlazE^y7OD z@AldKa?afg^-7TU+iTyPSK|_UA)gIeH~Z8bzVR%F&H(Fzem9jfKIuy2wC^3)kNWPj zbsB`8cHB4h{7y^fBJccZM>*?|g`9P>jtk+NxASFu&O<tK)^Q88am3}F&lAWw|HfyZ zv~QiOkT-uZZhuG2zxj2J={sk}naX=f4CR|=E8gv6ybsFlLe4&CLhH}CmEXyG{g~bq zXy;fx_f9?IY=XRU*4H{{-*@fC?f4yPkG)N?xU`ojAMcm;?WgmUNY8ujch&k`$=_(= zO9NZS;!WglJ#yAhe<R>0(z9Ofh5Z<UUc6q~8OVFQ@9H~$HI>(&eKr0>`g`%OKlQEq z<$sIRMfNDz9CY_qLAL;{xxH%#-(6`4?d%<f_C~4a&X48YZ>!&Qf0KAi8X$K)I3BbH zb)el1bMFgZ|JHsK`~q+&XuQ^JJAD23hSsn0?vVA^fZPO7zs~Mv_>(~O?1OdoW}ZOK z{**ymkB-pJw6kv9=06d+WYFC({{8T+yYls7ddH!yb5rQjnA~vaPGA*i{j5VD_|91# z^cv9lu)ogj9NtCiVZ7zMd+W@{`dSC;YhM1gt3cmfJ^;N3>;+mce~UX$+Oz+;$Xn+e z=((Wz`c0S#e;nvMSqI}VPxqh}^3H+vyav8?c5VvbJMW#LXMx&10&P9!L0d;}R3Uu( zF(_uoZ>oJMjlK7X*|nbD6>oDTa`wslvk(3f&~N%=Xy-?L`(Zp|kh8AN(+v2VK)-#~ zLw)Nw8o5rO@deP<)qU6nKLxZe#_ygwXWfwVPG1FWfBQk(pR1w$Z6yC%_&vdspnlB% z2KbFZ>tQ{O&pBv-oO@?{_SgE_SMTc-^yF(tJ^gJ(&iZ+e?8{2tn}X)yT{50_yo>Ia zdINcPFZ5&m?5FjWe-J(Umj-R${mrDF_uYJr*SMCVZ~ppG&N-8IPt>nKUq9*_r+&R3 z*5`NvJ>{qJ?%ZiV1-|vukMD<hH|}ZB*2DNV!8e|x(Ar5P@0{8n=fr)D+jqX5j~y{T zae3`K@7B$`*9d*{Tm)@A+o7Eg@5d_m#(M;MK4{*`^@nd=ddK99TY2}&IZ-bkJ^QHM z8u;qj2m9KF_j0fg82`@3{ad#w$ZOYmH7@O$w{dB2F#3M)yjxe&%SNvu__z9V{_Usw ziQ>$_kNc~9G5kb&*4=uF&V9UIHMQfs*Hm7A@%R(zJ4cD)P4CqxsrvtyF}Q8Q*_1m5 zpUn=MeVu>H2QP&GsCCU<ZQrVYF8lQ9;DP4%q~_E)mmU1B=;g8PQzG6s{{>&AZvFRP zReTUJgSGDu+&}0CA7=+|_UqWZ_(vaP2V<tBJP?Ftg4_0;yD$HZkE_M{ymhqr)aRkt zsqr7StacLIGc@Jk$9<};gm#*}_2$|R7o9h+Q$B2Yt<XU4iMi9>{qscwb^rDI^=pUB z(C6vk(5ah#@xW6bWuFdu{qemv_lALn%ySPm7}GTjHMptb+efy475chw?)lq04-7p7 z8DEX-wc$^prQnmBw~qdEQYiJ_AO85=f+6o`2g`F`J2d~}=<S8h*;n2CY2+_hQ2YDq z?~QC^&Uj?{sV?C^Mx@?Vdf`;~k1w`1y837DhKhG>xM#{8{lbiYa@U=wQmU>0@8etE zxllK>{`#*<@A=QX53_^quRQQ@%czo@UYmKW)x}7U{rIWpZvAao@#C}4v`P)dPX`bG z-LTS!t6f(octh_+_n&?2Om>h|)$jejUsi|o_AjS5|Cg&mzbntyz2_HWqJWRKT0iHF zt3u;{YP93+*Q0-Y)&1IK6%V~1r7-@FUA{NAL1g^+gQK3PI2_fX;eESbD+wbF@|In_ zdF)T$i{w6=c<T2nzK&dkMDoKweB=7_k?pdXSzBKV)r0Oym0!FZRcy>{Pd{|i|5XbH znQgP%rF<Cr3G&{!WBcQ`T;vYuzkOlr+M#H0^w-ZES@uzOaG>q^=2OGKLet}}{qcpD z9u9dE`<&bKT&J)$b4L~Z;%A}d$n)u7LYH@5eEq%mvV$|JyDs#<5SFJ`r?qeVsJdjQ zgXd12X?rFbibFrmz3I+srJ(PaRd+ABs2_amXMd<ylpMaldiC&kHbjFIJUKo0wgFKz z`#+jA;z$2{Ch`+xZ<#Wy>7uZTk2Sri{O_Yy-G6d<i=xlMTF)Od`}XoE)*$C(+RtB! zoW1kN_<=={(D0Hk8|=B*y~!h9zJBw?0Ny@-HuZ~o;qwn~|K!OP+oFMdzx9LXU;iR1 zV7*>%{Lkv@KsrtwH01rvvmx(0Y4^`b+ZN{d)bQLjCnC*}CqG^^eQns0PtT{He>a?f zps-}n^4%|m{}}uErh50>8`f#c@!Rq;zYGhK*Yfs9hea8U%i1;Yzr)Q4hHd!G`Rl4D zXv(7>EPA(f*wLWth)-Xh7VT{C>-STCUK!<)zvfWy&#I%kP`~{%&wmm{+O2K78~$ry z^sj<8D^4Xp9knp{-HYjKj%>Qv%}35Z_1S~f3A{4s=6lCS#l7z*H&1=()(=CSH(z~Y z$DJW7s_7#g-aa_vdto$hHrja4mZ)6W&wlakSGtCwA073;vrj}~!S3yQZ~VvmLZLA! z?O$sg#s2Ek553qSvikYCZ+)3rEmYa#t^*UJvOL}3(cRsmp$vZYeCq32pGE(=@a%=h z9<LXc@pSNVi+?!xLiO{vhyFhE2WLZTxxJpus|afo1{D4m-1W@IPqYjrF7bmGo^DvY zEULlj;NXpW_ElD|hGHSJBDK}TyQ5A0vTMQKo?#G?osj=Y**}$j9SNWQ#@5D{THkT* zop6fI{owBHE7nx+(M!Mou=qrH?*39L)a?A~jknI48y5R?pl<lj+kV@u$Cm@b(L2}T zZ0-Mu8f6xjq=FA``@G-KOVOMTzB}*hhwDU_B6#_illtrnZ`A4FmeohcwFo;E{ATHs zS$%g#{kbf$O9%eGcOO65F=~&Gm)wQ1%<S=L-q(+XDJgzBm|w8y#*d<_aw)1th1<qf z-4WF=)X4kRq)y|m3Y)xc&`VE!UVZ-*jGoJO)PLsIsGgs-p76z#>Sk&BlEz<)Bvc)E z`lq)vjq>@+K!ZA8ez(b+)mvH<YeA>f`EAEUV-a+@^T^$=McG`^`%CMg?@O;Uo*(va z4_*lS`B%kCA03zURQ0%hqoH7wz3;(u!*f3QG}_s3KCj=rh9gj;&#=})!#A$U$-Jrh Hi=6-e!!fj! literal 31333 zcmb7{d3;ypm4_3+u;!Am?}R0+0TN{2681IhAUg>pKuAI$2|Gx_zJz5)t4^zoR@>Uz znOduDU1+PV&U9ha>S%4NKdbE&cdFw!)~a)Fcpl2<&h(G*55GC*J@0$A?>YC!$XRnz z8zm(r4NpsMKCiN<xU#rtO3C`_?#)+!SaQ#$``$hK*qX{Uzgza?Bgw5eo?KB@R=mEt zq+(0elm_9=q*iIk&E^)C7Z)wwzO}5l``j<ptJ^e8PM*Iuxxv#9r7end4MWSc`JSdF zH=SQtSzWTethlOs^PhdX`-Sp{FK1s$Z~ECUetN8Fa=ZMBE!D-73oFYi@{0?rwpSLH z7jLP?Uz41=<fNomuYIi1!s5cR=<QdNg0@L1??qu<`(J4P|N5_OQma>6B_}7PB{!a5 zQBqY^v88+SA)}u8X36fizF)Pt_O7(KPvj)0PARFXt}I!%{Yot&2W_(&CA~T<{9oIw zCP|luh5u`tbS2!K`0oe~2G4;bz#?!3*a6%CmV>?{wvX9QK%Na&fg8c@;BK%ITnugj zXMyX%6JT3#0k|3L2b$kJ<St+*a0WOOoEek%BaZ~<f>S{A?u>j2JOoY#TY&ArfnXtM z-K_U@=<UNAWc!efTn8=%8-v!#JgtkecII(1_KtqUO`MzWYUJslb=BTJp5}ZGXq{Rk zo45IykME>lr(8ejq0amE&Afa^f7UaVcdWxA<O)#U+kf+vU-nIX8ux~Q*1rr{9+<!U z$>ZF*k49bs%46rG8+v(SUyZAl^BQmuXkV<$0rY3Ve6Sj{UJH?(JNs<C%u8NdkG|Y< zJ`N*WZ~Jas+Bd>3Zg0P=`%&yBfYxmkvhR;UmRIYL+ko2lKsN9G$nJwu<guW6<{~c# zt=FoUo%8K`=3yM;u(#jak#~UB-+Zi3OU|vI{TJ7Aem&S4l&{W<c{Jr*f8&vzH|ux; zy}Yz8);*8&0<Z}<5!B9lyHDyk?;W#qUt1UBv5xlRAosS$?srADF525C<7k1M^ROG) zJ+c_tIyleE(A!7vTlaF#oxf(t?yup<WuSH6h-|(+k&P>!k9O*fOCDJ-?X1Hd-gWQm z*K^}{e!Rbed+q`A8-iXQ`HuCKFZOv6_BEjXt+VxXE{m~qkGTiNptn!PYu?U71$O2! z8(F(*<msUGQNJI(^JE>opT~J6m;t)4mF0{5v`z=HcRsr#p9Xh;`ZEvfK7ez1ZQRaX z0q4s=^VhB=dgFI*nP(@?wb!qFFXem|Xq?(-pg#ic1?`LW`ZvG6*qsBjK;v@1^g}-{ zX0M&~_Fdyrwy*l}o_rX|yWZ1}^P~N8?0i>y>*f0cv6HXf_nxx%<;nb5e)6RVePPVM z`_6gp!nu8N-t-&yH#p{R1hVfLr}wY@w}zZt#vNzT-?oF$%YAo(y{h6o7u*Nx+k0|i zHRo+Wb8CoP0P5Es`diAHxWg7>uOD};HQd3uGd>R4-)!y-?dAAN?5smJa%0e$I)GdP z%28+BUYBuhy*eSAznpYu>fdj)9BRZpx$5qiiQb)|e`nR%Xo1~Qa1SUCGLenf`18@5 zul;pajMI1I=yC4p-~7$PZ@9l5)w?U4Df3b`fBpMy&cL7f&p_@4wgNYS?v5<vW?))O zPDQp}^1^uKoqhG(o$HP?ZfCIp@1}!CK=~*?%*XhhBj26Ry*f}nIM4RgIa-LFarzCH zhxR=gyRqOt(D(eepJ*Q0*x5JZGH?CMBkNYcy^Wyf@~Jy|^U`k%^wwn{vUQXn&Xx7D zKNGR{x2*d~evIMVcg7)GkFLloK;zkkY(2`5(?IKPAN037W_JX6D3}Gd1pPOPaWz5j zd*Y4g?ek=0^R-^?A-@YZU}t^vk<C93**x5j$6|KYcN+RBp!GY0Y<|v{{BZxvgPzzI zgYs3Ld%uoz`%!@G{*jl?jrM+*JU`Ao_ow;RpqKyh(mY3VZamJldGF-hxsZ?A4dc8c zXnv=W?UQvIh2FZj*BYaD4>{KXdgm+$IS(|>GUTb?eo+1lN0wLpkgb<}%S10v!~*oz zy)Ux)$Zz8`59?+>oGa0}bniFi-9Dgk8K-l$j&tkn_pJwd`EUN(IY0JeFZS}mxSdOR zXI+P3KLXTlBXTb2el@T0=;hB*WaHk7Y(MpLEN16-+PRY-bFkB&@9Jj}=f-;yd0Xs# z>*4(Rj(yeM`W14&612~gkd4dw`2Lufy>`~gxzoSjg(b0j@=HF*C+}IW2HewLjQd#? zd*42+MsFRxw;sLzoOkz^^>A*jqxa(P$e-=Jy8t{0TKBlWwVdm35b|Nrxcz=z`)>_< zXpX(O<IKrfbExIqZ{-o>so*Kl7!DzuYr~kVzCHR9a4zWk_I5aW<LrR!x6hsxpsxUZ zzX|dIus7(pw+(U+Q2QOod7!(b9kOUII-@tg>yiD|WFXty*2t-#92L8um&?wwvcIjo zZ{OT4a_21Xq=4=|{hyBcu`Y71JLmSpS>A<So}5IU49d++<W69p7$zedpZPUM@3+=( zh<(xj2<+ngu>Q`dvtYjR-yK@WJFP(Tsza^?-SsVzjm!AWTmIS4rr0~9*0(SEW1x1{ zT|UZF>t#GMxToJ{$nw#;7*99OwR0EC<Mo_dFYDn>_uFdx&VhZ~!hP+gAiE2#vwaxE z`AN{Yog3%9E9WJk^>&AQzL|5smF887-g;ZddFUsC=fL>fXy2Q2dEj0!zge7H7xOc& zC7cfgdxG-(5VC#IuXAIa^|umx_p`i}cX^x-0j;O!@@FyUXTWu!eH@D15R{kp$$Op` zV`rV@wf*vZKX&Hjd-6kGTbC)=F9Iin=H>mR=-mVEFL|+nbKkdL)=@ioW4zWkoBKUL z^D98MKgMH!&CB?FcO>`B&-c^O+y4W|^12Q3F0cTs0PTnGSr7Mtb0;72xOW3s1nvRd zA2rA=!FgagXrHauSoHGDepsh|oR0(7fXl$a;Bas>==;XuJox^N*kyuaK>KRm##g|( z@fe3Zvu@6f^KPHy`*hy%{b|VZ%I}YTThI9fP+oRJwx0IgI(vT*cKbl%?u~38-G4PP zyZ$lRI%lAFE_~O0Ig@ky<2+j@?cIAN*!%t_WbF%)jaNU;i{}~GodxY*KC<@O$uoIk zT=pS1cCQQa8PGlLJ^N^1wLgM=+>i5XT=F*q``Xw&^U|MjHp9+$<iGw;{f&2~qOSqn zw}+5-f^mPn(Cg2Bj6`qT^5fcnYjnq^1L)3nM&<l?&dqrq@)pqFvO|!aVP`<!_Rw8C z0DC#H9Jw2)pBs=Hf!4;F``gEyjaRJVUKwb<$;hj~wxGMi{05@eejIXp@G$7E?}}{P z=Hq+LviaFVcc5|DyUD!Q64d`*<P2~==<JGn(7W?mAvXrCt9_KK#<3x0*9qBr+=%SV zn3r*{;M{MU_pIMO&YOeAm4ZA8><n5zd8FR_oF(J14)R7m$n#aaD^?);Ev!R!w*6)p zzq0i*5BYArr|`})a1=Na+y=^<6Ug$({Pg2oiSo?46>(qw%7aw&_R$?b1ig0Zt&{J| zQ{&yqz0IKUwLqQ?HU!Og7P9frL2eED-f?7kCLf${{n@V`*qfiSJeDu^TVA&3p8YYt z8uZtJT|nz$obt>1SZC`vn0v;rKlxG0xqYzz9ntp#joUruZyVz&#!eo{pN;6Jg2TZk zpgi$?{Tiq7TL<@60r#zMfGofD*9X0E%|xC78i%+ZeM?aOxUYtxmp|K)_3M52l01=L z>YH=VdK^UdyJY-vp5$Ry3A+Ey&-dj;A$FTU`zN1vq8}Q2&${}a_m*Ji{VHVnum)Kk zuS1rv&WpbT%y&6<r@%qrTF`#%KsKNK$nm-Ci2eX*eMcb6FZ-5<-Z<SeW6>M8`_S*| zUe1q!#@7zn{o!7A-ZDAQ2IZge<)C+un4kU7UcULgFt0rB`;K<jL4Mdr_n`axDEG&K zzN@`)m2*B3tN=ZqfxHX!UF)Zv?{CFU9{IlY*WP%2w}yMpg?!T9JyMOGd&zohZ+`K8 z;`hk!lyRNpy%nJEYp>m5&gX%~)f0I&=soA}IC|&KdsEOGxAy1I$NgA8-!H*VzMhTU za}M<9{98xgv5&rg`ft3Gh2HzKkgbpVH12N^cKVaY)6mD`zV_c5p0&YtHaH9{1oi1R z##rZbZVdJy8GRqn*|ZnND3?Wd{!#870p-9*WWUYs)}846ZLt8^+WSr0j{XMFxa3ev z^!7dvxi#pn?T_rPFmHE>`8z9eQU4RTUkCODt*?DCuHKxxD~#KnX<Y79&*hx;Tg5xZ z-3-|sY~3cKcjrt*HlJ+duAqHrj649$0IR`npz(P>6}{g!`|nO@z`1exEnbIy8K~Yq z&P8v1+=++KF9qd^{Le>kyrYrVfacQ+IUPI)de3@wjOkO5{U&-o4!!eaT*J}JAM0Yj z#&W(JTn);*lgRScd9xn&c^Y>5X^U)sji(U3_vRzZEBn<Qz5E@5+!+i&`S09JMPCFu z7n_jfSx;nn)DXD^Xn&d_`=0(rqMr=z1f5IkF8|#x#^=sAjy&$WyUkbG_)4+U?>6Kt zP~PlEwh!`9zt%zi<Y2E~_lJEqZ|C6v_A|ltpuDob*3-F{r|ui$mUr^dddWlWt%p20 zfnWEmbuvzWJ2b*B53~>3%a0zMd*A*{KyUrzVFCIQ&^_RH#=Rsz<y&{`M}X#SJ&m`P z^Fq*lWqi(iI_FKnf#4+YdeHj&zJ2ohZ~!~wan6d+8=rMJgkHPm$o6$QvOLQ~?gh#h z`)S|z#qJ$M?gtJ4OTi|f@4JU@KySay+wX?wy|FWI>+Sr=|7O_v{&mRK*E;P$FAv=B z*3Y^R#%>iT&+WT&>E2hCNBWt^ee2*n4MHyud{;knIiCocZ!)rRPeyjH8Asf119sNa z{#bwS`>uVQ!o5CVBhY=M{TTH6wU5TF{b}sfk45(Psr4}~`Db2!&$QRhy=h&IOMWki z`SG56#<}!fI`*SM?VaPeANQYe$s6C1kKQvL{k6f5b@rVt=#9_4?K=xOci+ex=fU^J zV;A=~8@>Jtk$vB|-5b~bTVp%6##V;x48?QT-d*SJmXo`<zZt9p%}s8}!R4H*7k8r1 z0V`s5smO9>E3&)aUaUhu2OJB^Nq22_OrL_>AM_pX8>gIFf}Op#9>!D4x%Y}<@3%zG z2aV63wm@&5?s)r^#<@FQ4#}+%oR@&cVZ4*jy8~*FjoW+9ZaL?BL3d*zvO7&b_TRdA zzbE!H!Tq4$gc@Y)ehT>z==;9ox5B!6eh~YW;9yWc@<@AUziG^m=iAV`d(>MGWzk(} zUCduJ9{-KvZaImcTrkdu8R$2FQ$c^r_$@k$J|B$F(;4*cfmGyUpz+7^ozM9?a5(5L zb$`e!`(i!?*tZ7lzxA+R0q51Ad@&w(rTfC)9$T=tU%pd}-gE2io>;@Vb#q?3psxZm zLHjI^YSEj&^J{;na6Sk$9`~l-tWBKz-U{SKU}tb?Om98p!C21ScgDK|{Vq_weK4+m zoacduK<n-M?w{kFj{<vv?hos}9{m|`BRBx`+i(Ar<(+=bOP(2r_SVyLd2auv;x7+0 zPH``K`{o?_8`O8?cYo~7TRwJ0?|tw29d#eL$K2o6SsqQ~on@f*ewTNnm(TJ+{`#JE zFwYw9d2XEYFpqP2>u<vA&>L4nWcQ4A*4cNxSBJf_{`#Pomkp5fz!b0tI2)9A3y}MQ zYr%9dGltGpA^LM*XRs1H3Z4S(gM9NH=cpDt_0~u9UEk54_U3J#?kDHcxk<*4^HYuN z`=^oJx9)ZMqFoVo?rHhs_fWsu=VNc2Es=fK{%dc4^>3b~+|%Cul!4y3{C%2*-njN5 zyAR#>+G|&ao%_SMhN7PfdQYC2kNo#u1@_Ln_6PrlAN~8jaT)LS*gfY)f6lEu^d0Np zhWi_0?@U4N9XofQI){rmkNb<q?Ht)}<8}{S`)`fIylG9`MeZ&+ua7M3XM*mG<;c50 zcdoN$t$Zf~yAEIr(AjY3EJVK*)W12lL*Ee`3Z{WMpmDhC{T7|%Tzw;C-_Jm{hH_9l zXM72E?ha?<2K45+7uo#!AU6Rs!2_UmwomrQU1FcCyYHUhzWqCfJOaECl%sOjI?Uo+ zUMxe-0Npj#btQW1Z+%;#_kH6zg<hVlLYA{5k*h%KxF9BvMwTa=k)0j;l7~J6lrL+M z{T9qXE&{FhdSrKeDRM5@8nhnvOaJo8ytOwk<5|r+qrePs9cZ5N#J;%Oyx$vpWqDGL z-ks!qd9a6bfBUE(i{5X5^|GJMIPVI|H~ZBAy>Ztd+kg3U2EF|D+;^;}@yX{&-0KUv z`{y7RgXX^*c`;Z5)`9j{JHPSP*?ySsA?~U7H?4ckxU63e_VU|#^W6UUO`eLq`@p`r zkL0oT8Q4z)kATK|4%zz16Z^1*b9p15oh$Xmxg2}XXCqGp{Z4F0-Us@Q`HVxq3tR*q z1nrajQtv)%i=B1w{!#SyQU3bd!1~$uQP{i3<d1W(mvedQyYf-~Ip4l-Jl(i&o!cSz z2i>FkbDr}!?*L|l_BR9B_!=V{kN%B!4d=$;d}N}}1v9|$;7-uJm4z&?+@ssj8(%WA zb2Sgy{up0(^zr)Vqqj~2k>!o~Oh@lL$rJOoPu6iH_S#v$rs%H+%|m{=m*k6mGjI2T zby~?g*4OvkyZtyH47z{i{R#9n;8xK5(~zyBdFrPP=QF`7a5Si$@3^Oo-#MEav+si3 z8axhK*U2&a&dA0gZ=3_|opbrooO|}ec{eWkT8W)J(q233>Yg%h_ndin?*Q*j1l`x( zv+voQYwv!D`>|fWZyx)&?>u==-s?}EHpAX`x*@Ls-S@uZzF*F{^^iaN(Z~C%zj)mC z(fwlF_VwC-Yq*Oh;6qL=Mb@{qcUGJwb9JUxanId&7TH;9fo!h><c8oRP>wm%CFtE@ zJ&`wn=2nMnyseOZuK?NF%5m-G>KW|hKm+6qa3GimT65#P4!!=J@#*OIfpWpR^+7MX zi><qLGA?(W@2ugz{c~4~)@wa>a>#GMcJz6m`V?e$xBfhL*Op=@kBn2EIs5i69eclZ zEn|M=O(uG0#a(RO+$~+OGmcc`WuX08f$T2!eP=~K&hBRH&AS3w{@U+_=&frZvh}M* zJ^`9{QS6@Q@w}Q~HxrZ}*4cg(aBkhLvpm*rD0b>acaZm6VW<Bc$mZpHqIo(O-s{J` zv7mcr5wi2%AGsK`kM_%aeRnH%^2>Mq*2Vd-8v8Wx3~1k{A!mU4sY13tCCD|Pb#903 z985-T3cCNSul=|F&9QT?jw0I^_m=UB`qkb(%;LWBZbEkNjYsyKG04_=9I}3`cNTj0 zl6CZ5{VAWtelDnAdE}hQFYOy)e-N}E#+ifOymumdzXx&v4hF54eQ{6tzV)zvdEA=; zTF2SQ8$jdifb4tLN4tHT%Y)&_%|Pwk1I|}B&i!7A>(HCeL}cr79NGE~L)OlDSb*Mn zDo0)f&IIL0A#ygD0y_V`uV3%!&+pq7?v;Vo*ZFb2uIF6)!^rNbBgoE=JaHe$%Z=FC zAN~6Ne9pCVKCMSB=gzbBAB4UQI2Wt}<<(&10JM)&k!wKva0<C4I1jY$Cy@_<*3Gyl zp*KExGy{Fy&OI$pQ?YZOYVUXC0O$6>zE{NFmv_#;eX<U5`+V-pyQQ%>jLUv`-+E}T z-3{D#uFS*lMjq$!dkfLa6Yq6JUjXW79dg`{yzqVbE5Drw@A<n%f5z#ah`%Eromby+ z9_5Gs=5QX*67L=`?$7U;{`O<n5;X2@$k+Z`qZPV3&~J*nCIh{(x{I~j!g(?{2{gyu z$j+2IwI_OKbrrI?J8OC9M}WTL?6pJhZngH#tULAucG_FR-sqjl70A}bS=)i$-vHK9 zZg=9`p1I46<2ufrC3n8|a&IVh#_Kmszt+pR<c{@U%YEy!6xn|H4ah(*x8$Dr7?1H; zkA~cHmgGqm`T?NdSZ87c`W~Qln2T)v`XcWI%R%Fm+xE?K?dM@X1(c7Qk<H6*iF{bd zxp|4ZV)|3a<~tEtf7T-nz4e#Rr_on~2S8_7Y>j>mX#b7N{Cjck+{lB2=#_^c%Mbnf zo_^$M4fg8uk=;SgtLJ?<H%{l{Bzo&3?^mL?jy1^gQ@&cKZk*ds=hr#04~Manr`B;X zdi(FU#&hHLH^(OIt;ZnbT+ld+kjueha5m^()n0zefBS7-`uDeu{1hAUZXq}Vw9nQ# z554bMNB2up&VA4Pm!Piz?c-_W4A40&MfRP6$kxfeSr30x7?*p}{o?*u#5?^#zn${e zIQwyKot?kcvHM4m<&FDp4to3H_sRR_V|*>JKM4B0u>Ld9?*!dD^4Wa4a_*kYM%Isc zZi(sTue@2tx%ta`>uf$}u#?B`^~va)fuq5rpx*=gY28n8Ze7cekAudaiEMn<!FjdL z?n~|EO)u{60G$i@@7y`x>aC;k9^=0Cme(1v_pOU_CBN2VCr?t5Mf2%_KE995a|7q& z!C9btXb!Ub)jsE=m%sWwhu(fUXO-wnLE|<rc_UwqOFR9^2YF>3j9+{Asr)R!ul#qv zY43hHgk8M9+NX0~3l@XUc}ryD65SWZuf6A6xfjnv{(9~`_f!G*?6>*Ick8dcb&UHl zKi|*bUHfi-ych4I^VE)ezT>>hKfklSV_eDH8xBqbcY}jL^V6U6tG^8Fwu8p4UE8b{ zNjc|2Z;4>>O=U0M^62}yi6FDx{a<)C`EqVBx7|yZ^HRb;miD-`<F4kHa}$X~@Mvzz z)0v;Y7uvnI>$mR@Xc8Xp`S8~J4km?ve6s1E|4?{c_{R@!Y54nj(aj%hOv&8-?7O+a zn>EXKpNO8X|1WsAb?v0@910Eo^3yhp-o5x<ZgAJB%|Cl`);qbu#e-v)*M*kBm#TmI zg$?x)?HpauWmKg4>E+BbBZJ6!Qln03`y&r&+rDw&(+1%m^IpHF%`?#-hxb4Bmz(PE zEgQY=sqD`}^F;8_7ni@D)a&i==A2i4J-*=G=#A>bf4*zLo00DR3qRg}vFMG+@QK2V zl$}uy`3qiNGHzk`M|$>@n@?xF6&@yni|^mE`IeW%aPMhQx_CpU@br_Mhkrjesz)OD z^3Qhs^4K@u3hR~jVBUZ~gnyK_Te&LbMC5D1s_&IQ(k8T@-v0L5KlBJ|8B7>5V{E@h zQB{Law+(#uwT@Q`T=3H4ImKZfi6HyL!<#;e3Vio(s`?LD8|MG<u?I39hzgN<XvN$i z;ZY*Er{A_0KZ%Mz?1g{3{*$qhD&?kk7QNIZ`l~_eva;5fqu>(3pa-gMU(qzumf!cM z4{m)YOd|iwuRnL^H^V>P=soYBUk_WJ2v!V#sqvo7@VWl~y7FMVPr|H&yL-HS*Zn^V zty=AQGwZq6qHNCIGw6;xBl|?qeAokDyD@A;BDnMYkI(-kj3f~}_~S=!uc!>WUft!c zXEuhB1mEx3W5dS!c;|ioXx8p4J?Ylx-v?eSemgg~{TGvVfAs#F;mx0X>o*0Bqw4Q_ z=4-uD>%vTHn}2WF!i-RJ=<*%&*E|-6VKu^}LjNCBetXQ*E8Y$7{@u5~-sR_?h6Y#f z-Zo_J?#ENZyEor|uII)2(C$pl-d&z^?JL27uK)J*<i=6lzx+X)+k4eVnRWT0w)aGZ ze)Qqa8!Hxs_V*9@(S_S0&EDF|%X^Q8n!c5thJFwgG!e{Cf9;Q7t*6`<es|=JyHmr~ ze}BlP3y1z8wC|Xi)4uh8h4vFBeydHZW|5y3=bOCorN4x#Ww%{QUm5i}eB;xvq%3F= z4q1?Xa?J;gA|gMr|NRRE(H|z12+Gb+Ib85z)Pm;<_pbj@SjXUt+lpqt_Ve(Mzg#?Y zeoK^X$4`sDdiS94k4-~Pb#K)>oT(3f{Q1Z)2Uq-jW7hrw$Jd7TiNEXq{)1r@k@M?6 z|M`x8in6n~SA7Nxm#p0Ud_=H{w?Db({$Wx4Cks+;&HB&qiTx|OZY_>-`(dlczkYN; zXw-IJvuB_CQ|Pkil1>-y{Wv^&qWIxs|JtzLcOp3W(!CG=>XR!zQ{Mi(>7pM-<2Ng% zZQom?jni_~w5p@wFa<4-d^7L5KZk#O_t%$t&!}fw80CczM!k5pp2GLu*>mO5>EVsC zpOoJH#rn$s_RNQi9^86mOgp_WbmXuz(K7Qqc;VuGPh>_sx|-qdi@*P5ZZyK*JN@U{ zhW{NpD_(f%I~RTx{?XvWhaYbmQAy*lMxp<s58XWGhdE)DZhG*6(p#5?%>8}(N5h7` zUY}1Q_^id++duk6MB}pZq{gKorGoo%J5;<<KUgpSZ28PKEkf1j|I+WfqeF@(u09?N zsvG*q<sIR?Jek^Q(-SX-e(iYBdFfA1Rxb{ZE@VwBD~}RQ1npia`O{<XMNWgh&p-X! z*RrDG26eSLH@=Y;mGh2A=X4!$#dJ-g{eLV9y(WU8eY$>qaQo1Id3vV-JwJ}@9w~b2 z4@08P1nW~jE#EURJnZw?1B0hWYjyCEVb9(Ftt*XMdDrR*N9&W=)bQJ1ZIT`yC%yjB zy-Pwa2WQ^PYBD1t&QlMyZg#$&-Ie>lG^6(PJK?jpJpF}XkB$!0nVS0eq2rO&)p`aO zGCsTG*AYJ+Y&{{d^HRu|IUhgu!0(<8HR;*=`u?Dv$UQohjb0Ka-|Q=oJY4Z=Xmf0M zdFlLULKDH}^6#Y9MBHmsaIsBxG$y~?dU0s_P2m;~PR!nxcVEQW+n&1Z-nQ3=f9xv! z&Bv)Bp%UTEe;HT&-t%E@!JO_*PW?JE`*ZQnx179EF~dm&?_|H-VcOxa$*Eu5{p`HR z)1nRC=9j)04!1_ZKRvYL`LER1x6gvJZD-GWBRBZ(mtJ1~M7J>3TfR1RRkJ^a_pe&# z?;Uq;Y{c{N8{SIX5w8AVSIhIgCPb`w;`1YuE<GN`_gCvT$2NX_aQ&<VUHa60e!g3M z+TY&v{DhizQPZw|r~h=m`{Zb`Ty>oYR*!jlYsv0tD_*rSu0(M2hyR!}D_XKw1ADU5 z@4oTJdcGch=9M?1(qFyy^;agovNUSlukIcGa{E?qhM&Hq@9g<*RPU?Tzf$yE=CP1S zSKAwSD-k@l?Y=#CMEAA+Ye<P;#8VyrbLwxRtgo8A6+E@8^jiD>>z9e+k+|mK>~-hc z75^e6i@#8U<JF%Q9FNNVt#yyw{`Y;uKmPYG6{GsUyj=bLQ}E=wFN`h-8J!4L9lZHe Jluja%|36u~a(n;) diff --git a/IntervalsModel/src/IntervalsModel.jl b/IntervalsModel/src/IntervalsModel.jl index e79d487..5d73060 100644 --- a/IntervalsModel/src/IntervalsModel.jl +++ b/IntervalsModel/src/IntervalsModel.jl @@ -22,10 +22,14 @@ using Plots const rng = Xoroshiro128Plus() -#lets us use nicer names for the indices +#consts that let give us nicer names for the indices const YOUNG, MIDDLE,OLD = 1,2,3 const durmax = 144 -const comparison_samples = 100 + +""" +Number of start times to sample for a given set of distribution parameters. +""" +const comparison_samples = 70 const sparam = [60, 12] # Swap age brackets for numbers const swap_dict = OrderedDict("Y" => YOUNG, "M" => MIDDLE, "O" => OLD) @@ -98,6 +102,8 @@ function bayesian_estimation(fname, err_func, priors_list, dists, particles; alp end """ + do_hh(particles::Integer) + Fit HH durations. This function specifies a list of distributions to fit, and loads the priors for "hh", and then calls `bayesian_estimation` with `err_hh`. This function passes parameters to `bayesian_estimation` depending on what produces the best fit for the HH data. """ function do_hh(particles) @@ -115,6 +121,9 @@ end """ + + do_ws(particles::Integer) + Fit WS durations. This function specifies a list of distributions to fit, and loads the priors for "workschool", and then calls `bayesian_estimation` with `err_ws`. This function passes parameters to `bayesian_estimation` depending on what produces the best fit for the WS data. """ function do_ws(particles) @@ -130,6 +139,9 @@ function do_ws(particles) end """ + + do_rest(particles::Integer) + Fit rest durations. This function specifies a list of distributions to fit, and loads the priors for "rest", and then calls `bayesian_estimation` with `err_rest`. This function passes parameters to `bayesian_estimation` depending on what produces the best fit for the Rest data. """ function do_rest(particles) diff --git a/IntervalsModel/src/hh_durations_model.jl b/IntervalsModel/src/hh_durations_model.jl index 2f2eebc..fe80396 100644 --- a/IntervalsModel/src/hh_durations_model.jl +++ b/IntervalsModel/src/hh_durations_model.jl @@ -1,9 +1,6 @@ const HHYMO = DataFrame(CSV.File("$PACKAGE_FOLDER/network-data/Timeuse/HH/HHYMO.csv")) -const cnst_hh = ( - subsize = size(HHYMO)[1], - Wghttotal = sum(HHYMO[:,"WGHT_PER"]), -) +const Wghttotal = sum(HHYMO[:,"WGHT_PER"]) function make_dat_array() durs = hcat( @@ -27,28 +24,34 @@ function make_dat_array() end const dat = make_dat_array() #assign a constant data array +""" -#error function + err_hh(p,dist) + +Calculate error between observed total contact durations in HH data and simulated total contact durations based on `dist(p[i]...)`. This one works a bit differently than the WS or rest error function, since we do not have a distribution of total durations to compare to, but rather one average total duration, and we do have the explicit samples of numbers of contacts. +""" function err_hh(p,dist) p_matrix = zeros(3,3) for i in 1:9 p_matrix[i] = p[i] end age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD] - duration_subarray = dat.durs - num_contacts_subarray = dat.nums + duration_subarray = dat.durs #durations, these are call subarray because we used to only sample a subset + num_contacts_subarray = dat.nums #numbers of contacts for each duration # display(p_matrix) - AGERESP = dat.AGERESP + AGERESP = dat.AGERESP #age of the respondent errsum = 0 - @inbounds for i = 1:cnst_hh.subsize + @inbounds for i = 1:length(duration_subarray) #loop through entire file age_sample = AGERESP[i] @inbounds for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages + + #get num_contacts from the HH data durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],num_contacts_subarray[i,age_j])) .% durmax + + #this returns the MEAN of the distribution of total durations, given we have durations given by durs, and we sample comparison_samples from that distribution. expdur = tot_dur_sample(comparison_samples,durs) - errsum += (expdur/comparison_samples - duration_subarray[i,age_j])^2 #compute total + errsum += (expdur - duration_subarray[i,age_j])^2 #compute total end end - # display("2") - return errsum end \ No newline at end of file diff --git a/IntervalsModel/src/interval_overlap_sampling.jl b/IntervalsModel/src/interval_overlap_sampling.jl index 8918004..c64d83f 100644 --- a/IntervalsModel/src/interval_overlap_sampling.jl +++ b/IntervalsModel/src/interval_overlap_sampling.jl @@ -1,4 +1,5 @@ # Set the underlying parameters for the intervals model +# This is the distribution of start_times const startdist = Normal(sparam...) @@ -10,7 +11,14 @@ function coverage!(cov,S_j,E_j) push!(cov,Interval(S_j,E_j)) end end -#compute the total duration of a sample of intervals + +""" + + tot_dur_sample(n::Integer,durlist::Vector{Integer}) + +Returns the mean of the distribution of total durations, given a set of contact durations `durlist`, and a number of samples `n`. +Note that this is different from `tot_dur_sample!` in that it returns only a mean and does not mutate it's arguments. This could be reducible to that one but not without losing some performance I think. Used in `err_hh` where we only need a mean. +""" function tot_dur_sample(n,durlist) if isempty(durlist) return 0 @@ -20,7 +28,6 @@ function tot_dur_sample(n,durlist) int_list = Vector{Interval{Int,Closed,Closed}}() sizehint!(int_list,numcontact*2) - start_matrix = trunc.(Int,rand(rng,startdist,(numcontact,n))) @inbounds for i in 1:n empty!(int_list) @@ -32,10 +39,18 @@ function tot_dur_sample(n,durlist) union!(int_list) total_dur += mapreduce(Intervals.span,+,int_list) end - return total_dur + return total_dur/n end +""" + + tot_dur_sample(sample_list::Vector{Integer},durlist::Vector{Integer}) + +Mutates `sample_list` to contain samples from the distribution of total durations, given a set of contact durations `durlist`. + +This is different from `tot_dur_sample` because it modifies it's arguments. Used in `err_ws` and `err_rest` where we need a distribution. +""" function tot_dur_sample!(sample_list,durlist) if isempty(durlist) sample_list .= 0.0 @@ -45,7 +60,6 @@ function tot_dur_sample!(sample_list,durlist) n = length(sample_list) int_list = Vector{Interval{Int,Closed,Closed}}() sizehint!(int_list,numcontact*2) - # @show rand(rng,startdist,(numcontact,n)) start_matrix = trunc.(Int,rand(rng,startdist,(numcontact,n))) for i in 1:n empty!(int_list) @@ -58,12 +72,3 @@ function tot_dur_sample!(sample_list,durlist) sample_list[i] = mapreduce(Intervals.span,+,int_list) end end - - -function as_symmetric_matrix(l) #turn a vector of length 6, l, into a symmetric 3x3 matrix, probably a nicer way to do this exists - return [ - l[1] l[2] l[3] - l[2] l[4] l[5] - l[3] l[5] l[6] - ] -end \ No newline at end of file diff --git a/IntervalsModel/src/rest_durations_model.jl b/IntervalsModel/src/rest_durations_model.jl index e40e217..53678df 100644 --- a/IntervalsModel/src/rest_durations_model.jl +++ b/IntervalsModel/src/rest_durations_model.jl @@ -6,25 +6,42 @@ const kerneldensity_data_rest = ( M = kde(rest_data.M.durs,weights = rest_data.M.weights), O = kde(rest_data.O.durs,weights = rest_data.O.weights), ) +""" + err_rest(p::Vector{Number},dist::Distribution) -#error function for rest +Calculate error between observed contact durations in rest data and simulated contact durations based on `dist(p[i]...)`. +""" function err_rest(p,dist) + #we get the vector p as a 9 element vector from KissABC, but more convenient as a 3x3 matrix p_matrix = zeros(3,3) for i in 1:9 p_matrix[i] = p[i] end + #sample a matrix of neighourhood sizes for distributions in rest_distributions + #this is actually not correct? neighbourhood sizes should be symmetric! neighourhoods = rand.(rng,rest_distributions) + + #create a matrix of distributions of type dist, from the parameters in p_matrix age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD] + + #initialize the list of samples and the total error sample_list = zeros(comparison_samples) errsum = 0 for age_sample in YOUNG:OLD for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages if neighourhoods[age_sample,age_j] > 0 + #get durations from our candidate distribtions for each of the contacts in neighbourhood durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],neighourhoods[age_sample,age_j])) .% durmax + + #this MODIFIES sample_list to contain samples from the distribution of total_durations, given intervals of length dur. tot_dur_sample!(sample_list,durs) + + #compute a kernel density from the list of samples kde_est = kde(sample_list) + + #compute L1 distance between observed distribution of total durations from survey, and observed distribution from simulation errsum += mapreduce(+,0:1:durmax) do i return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i)) end diff --git a/IntervalsModel/src/utils.jl b/IntervalsModel/src/utils.jl index 5215ce8..a02d85d 100644 --- a/IntervalsModel/src/utils.jl +++ b/IntervalsModel/src/utils.jl @@ -1,12 +1,32 @@ import Pandas: read_csv using DataFrames +""" + get_params(params::Vector{Number}) + +This function reorganizes a Vector of 6n parameters into a 3x3 symmetric matrix of n-tuples of parameters. +""" function get_params(params) p_list = [as_symmetric_matrix(params[i:i+5]) for i = 1:6:length(params)] return zip(p_list...) |> collect end + +""" +Turn a vector of length 6, l, into a symmetric 3x3 matrix, probably a nicer way to do this exists but meh. +""" +function as_symmetric_matrix(l) + return [ + l[1] l[2] l[3] + l[2] l[4] l[5] + l[3] l[5] l[6] + ] +end + +""" +Load rest data from `network-data/Timeuse/Rest/RData`. +""" function get_rest_data() path = "$PACKAGE_FOLDER/network-data/Timeuse/Rest/RData" data_table_by_age = map(collect(keys(swap_dict))) do age @@ -18,6 +38,9 @@ function get_rest_data() return (;data_table_by_age...) end +""" +Load WS data from `network-data/Timeuse/WS/WorkschoolData`. +""" function get_ws_data() path = "$PACKAGE_FOLDER/network-data/Timeuse/WS/WorkschoolData" data_table_by_age = map(collect(keys(swap_dict))) do age @@ -29,6 +52,9 @@ function get_ws_data() return (;data_table_by_age...) end +""" +Load priors data from `IntervalsModel/network-data/POLYMOD/AALPoisPriors.csv` and format into a DataFrame with some additional columns for the index of the age class (per `swap_dict`). +""" function get_priors() df = DataFrame(read_csv("IntervalsModel/network-data/POLYMOD/AALPoisPriors.csv")) df_w_indicies = transform( @@ -41,12 +67,21 @@ function get_priors() return df_w_indicies end -function filter_priors(key) + +""" + + filter_priors(location_key::String) + +Filter priors dataframe based on location_key, which should be one of "rest", "workschool", or "home". Returns an array of Categorical distributions. +""" +function filter_priors(location_key) priors = get_priors() - if !(key in priors[!,:location]) + + if !(location_key in priors[!,:location]) throw(ArgumentError("location key not in priors file!")) end - priors_probs = filter(row-> row["location"] == key, eachrow(priors)) |> + + priors_probs = filter(row-> row["location"] == location_key, eachrow(priors)) |> df -> map(r -> (r[:index_out],r[:index_in]) => convert(Vector,r[string.(0:143)]),df) |> Dict |> df -> Dict(map(v -> v => Categorical(df[v]), collect(keys(df)))) diff --git a/IntervalsModel/src/ws_durations_model.jl b/IntervalsModel/src/ws_durations_model.jl index cafdb6c..91aadca 100644 --- a/IntervalsModel/src/ws_durations_model.jl +++ b/IntervalsModel/src/ws_durations_model.jl @@ -1,29 +1,49 @@ -const ws_distributions = CovidAlertVaccinationModel.initial_workschool_mixing_matrix const ws_data = get_ws_data() + +""" +Kernel density functions based on the durs in ws_data and their corresponding weights. This defines an interpolated density function that we can compare to the density obtained from the simulation. +""" const kerneldensity_data_ws = ( Y = kde(ws_data.Y.durs;weights = ws_data.Y.weights), M = kde(ws_data.M.durs;weights = ws_data.M.weights), O = kde(ws_data.O.durs;weights = ws_data.O.weights), ) +""" + err_ws(p::Vector{Number},dist::Distribution) -#error function for ws +Calculate error between observed contact durations in workschool data and simulated contact durations based on `dist(p[i]...)`. For these workschool simulations, we use a different sampling procedure for neighbourhood sizes based on `ws_sample` in CovidAlertVaccinationModel. +""" function err_ws(p,dist) - p_matrix = zeros(3,3) + #we get the vector p as a 9 element vector from KissABC, but more convenient as a 3x3 matrix + p_matrix = zeros(3,3) for i in 1:9 - p_matrix[i] = p[i] + p_matrix[i] = p[i] #this will fill p_matrix with the elements of p in column-major order end - age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD] + + #create a matrix of distributions of type dist, from the parameters in p_matrix + age_dists = [dist(p_matrix[i,j]...) for i in YOUNG:OLD, j in YOUNG:OLD] + + #initialize the list of samples and the total error sample_list = zeros(comparison_samples) errsum = 0 + for age_sample in YOUNG:OLD + #Sample neighbourhood sizes. Note that this includes the weekend/weekday coinflip, implemented in CovidAlertVaccinationModel. neighourhoods = CovidAlertVaccinationModel.ws_sample(age_sample) - for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages + + for age_j in YOUNG:OLD #for a given age_sample loop over possible contact ages if neighourhoods[age_j] > 0 + #get durations from our candidate distribtions for each of the contacts in neighbourhood durs = trunc.(Int,rand(rng,age_dists[age_sample,age_j],neighourhoods[age_j])) .% durmax + + #this MODIFIES sample_list to contain samples from the distribution of total_durations, given intervals of length dur. tot_dur_sample!(sample_list,durs) + + #compute a kernel density from the list of samples kde_est = kde(sample_list) + #compute L1 distance between observed distribution of total durations from survey, and observed distribution from simulation errsum += mapreduce(+,0:1:durmax) do i return abs(pdf(kde_est,i) - pdf(kerneldensity_data_ws[age_sample],i)) end -- GitLab