From 3acfd712232bf36b91f3ec3906e545ced2c2cb6a Mon Sep 17 00:00:00 2001 From: Thibaut Jombart Date: Thu, 7 Mar 2019 11:12:36 +0200 Subject: [PATCH 1/8] fixing the NaN likelihoods occuring when first case count is 0 --- R/get_R.R | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/R/get_R.R b/R/get_R.R index d58a8d8..b36f480 100644 --- a/R/get_R.R +++ b/R/get_R.R @@ -167,16 +167,29 @@ get_R.integer <- function(x, disease = NULL, si = NULL, x_with_tail <- rep(0, last_day) x_with_tail[dates] <- x - all_lambdas <- EpiEstim::overall_infectivity(x_with_tail, si)[-1] - dates_lambdas <- seq_along(all_lambdas) + 1 - x_lambdas <- all_lambdas[dates] + + + ## Important note: EpiEstim::overall_infectivity assumes the first case is + ## imported, so we need to tweak the data provided to ensure the data passed + ## on to the function is the first non-zero case count. If data are all + ## 'zero', then only the first data point is ignored. + + if (sum(x_with_tail) == 0) { + dates_to_ignore <- 1 + } else { + first_cases <- min(which(x_with_tail > 0)) + dates_to_ignore <- seq_len(first_cases) + } + + all_lambdas <- EpiEstim::overall_infectivity(x_with_tail, si) + dates_lambdas <- seq_along(all_lambdas) + x_lambdas <- all_lambdas[-dates_to_ignore] loglike <- function(R) { if (R <= 0 ) return(-Inf) - sum(stats::dpois(x[-1], lambda = R * x_lambdas, log = TRUE)) + sum(stats::dpois(x[-dates_to_ignore], lambda = R * x_lambdas, log = TRUE)) } - GRID_SIZE <- 1000 R_grid <- seq(0, max_R, length.out = GRID_SIZE) From 1041ad008646639b7b08be037155c1388a4318bc Mon Sep 17 00:00:00 2001 From: Thibaut Jombart Date: Thu, 7 Mar 2019 13:47:00 +0200 Subject: [PATCH 2/8] adding test for no cases --- R/get_R.R | 40 ++++++++++------ R/internals.R | 8 +++- tests/testthat/rds/R_1.rds | Bin 14969 -> 14998 bytes tests/testthat/rds/print_1.rds | Bin 330 -> 327 bytes tests/testthat/test_get_R.R | 82 ++++++++++++++++++++++++++++----- tests/testthat/test_sample_R.R | 27 +++++------ 6 files changed, 117 insertions(+), 40 deletions(-) diff --git a/R/get_R.R b/R/get_R.R index b36f480..54f109f 100644 --- a/R/get_R.R +++ b/R/get_R.R @@ -114,6 +114,11 @@ get_R.default <- function(x, ...) { get_R.integer <- function(x, disease = NULL, si = NULL, si_mean = NULL, si_sd = NULL, max_R = 10, days = 30, ...) { + + if (sum(x, na.rm = TRUE) == 0) { + stop("Cannot estimate R with no cases.") + } + dates <- seq_along(x) last_day <- max(dates) + days @@ -169,25 +174,34 @@ get_R.integer <- function(x, disease = NULL, si = NULL, x_with_tail[dates] <- x - ## Important note: EpiEstim::overall_infectivity assumes the first case is - ## imported, so we need to tweak the data provided to ensure the data passed - ## on to the function is the first non-zero case count. If data are all - ## 'zero', then only the first data point is ignored. - - if (sum(x_with_tail) == 0) { - dates_to_ignore <- 1 - } else { - first_cases <- min(which(x_with_tail > 0)) - dates_to_ignore <- seq_len(first_cases) + ## Important note: we cannot compute the likelihood of the first day with + ## case(s), as the force of infection on that day is by default 0: the + ## likelihood is in fact conditional on the first day with cases. The + ## resulting log-likelihood for the first day with cases will be -Inf. We keep + ## track of this day and ignore this data point when summing over + ## log-likelihood. In the general case, this will be the first element of x + ## (x[1]), but we also implement the general case where the incidence curve + ## could have started before the first cases. Note that the force of infection + ## for day 1 is also by definition not known (NA), so the first day of data + ## (x[1]) is always ignored in any case. + + ignored_likelihood <- 1 # fist day always ignored + if (! sum(x_with_tail) == 0) { + days_with_cases <- which(x_with_tail > 0) + first_day_with_cases <- min(days_with_cases, na.rm = TRUE) + ignored_likelihood <- union(1, first_day_with_cases) } all_lambdas <- EpiEstim::overall_infectivity(x_with_tail, si) dates_lambdas <- seq_along(all_lambdas) - x_lambdas <- all_lambdas[-dates_to_ignore] + lambdas_for_x <- all_lambdas[dates] loglike <- function(R) { - if (R <= 0 ) return(-Inf) - sum(stats::dpois(x[-dates_to_ignore], lambda = R * x_lambdas, log = TRUE)) + if (R < 0 ) return(-Inf) + all_likelihoods <- stats::dpois(x, + lambda = R * lambdas_for_x, + log = TRUE) + sum(all_likelihoods[-ignored_likelihood]) } GRID_SIZE <- 1000 diff --git a/R/internals.R b/R/internals.R index 2e50d32..764aeee 100644 --- a/R/internals.R +++ b/R/internals.R @@ -54,10 +54,14 @@ ## Sanitize very low log-likelihood values and put them back on their original -## scale, summing to 1. +## scale, summing to 1. We need to consider the max of log-like could be 0. loglike_to_density <- function(x) { - out <- x / abs(max(x, na.rm = TRUE)) + out <- x + x_max <- max(x, na.rm = TRUE) + if (x_max != 0) { + out <- out / abs(x_max) + } out <- exp(out) out <- out / sum(out) out diff --git a/tests/testthat/rds/R_1.rds b/tests/testthat/rds/R_1.rds index 2a968ce1f6d918768e59bce82b30da03226260ff..76cd5c7d6b4f5314a9871d30e97c4a1e71e345ab 100644 GIT binary patch literal 14998 zcmW-nRZv@98-;Oq4Nh>UxRc^g+@WZ3cXxNE6o(IKaV@Sv3dM^%3B`-M!=L$cv!7Y- zv)9>k=3FG1Lz9RG_x}&#X)vk*(b(E;8XuO(5E2m&9%F3#Z;{b2@?&(0C_LF%3&r?N_)a+B_{Ac9-}zedBDY=Mk0SoB0C5BH#b;|MRtRg5Bm|7p zlSb1lCG{5+a^^0T0O;fv1m+qjZCfB#+-=ui0F%S3JFOoscrV%+n&LuT&tj{iwMsK~ zc;i8p8;(O_mSxRxg+18YISXT~hKNR`hUl`7EtOxxPouupb|jrtGW1LdxqTG3@wnP8 zudIza=_pYhVUV-N;TLJt<^VD>0R_1jna-Z0(fQv6NBl)vK1m9NunZqhktVBz_;bUF zUEXcg1rR;ozPRC?gRTsB48972C`U)twY>A%zF!a78`m$mCHsVIB?{YBxi}kB$US&(s;t}!cRCuw*mX}Qa$A9V zCVIW&tvjMxRtdEualP~4dtbl8Ht)vJ=_dDrc(C>iJT=sZX2H(P%r=mW`tKdh(H|NH zBw;Z9Mrt*KP2pbD#0B8~;PUbd&9IFym~f+k*~~kzIEv;OsR?shOA6H^(TBZ;t#E1q-c=!Gb zShi{#QVufV7Rn@EfI;vrPM%GT7eFJm;b))H#D0DVKR@)?*`W3mVB17Re&O9u*S({G z+i-E)-=%)fcWd*Jo;r5KN&febKX!J*JF~Jr%noAw+wS(F4TPEERn(?Cm;O`|9NU_L zh^9O44kXk?Hu;3uQ+_>piQYebIk#}8{5G~8l5!Yz4RBpA9(wH;^Tqrz3Ag^b2XzPR zmL3FJLOH2teD%t{iP5aXpI>zJ)ll_5fjd$BQb+HD&vH}~gmNyZQ?FyUpTul7Q^Rrn zdO%+(PMm93;b|^ALREOgPc|3eXf7W@Gx*p5eC6{2&_;%x|E+8R)6f3&F)A$_H#7z@zTuNY2 z0dT3du;jx3d;SydKHC)Cv{bMN3u&l&Fk{q`GSw9J_c@zJ{={F-btV|$ds%FF^)0P3 zW@Nyd06-?+hU5)w>%+K`%IaUfJZn`@M{^}r)E{i1N4q4|xOjZ_G%BIksd*+f(T=6u^!o9hh(jR8C8?d-k4>O?Vi z6gQyVOimu;-kXDN#uL)tJkn@uYy`uU{xb^gz{(du&@aEYCsk1O_B7%{jQi`9e!u49 z7mA>vhqtyV^U?!te^~lfHwE1+AZ7X!d%1Psl;d%tOMCNWs;pv(I_-KbldQcdZ>V&u*HJpp?Wn(ve{1{#GCcwJiG;E^ z|2(7wWVe4|mNL{c%Qpiw{oEs2*q3Lz*=(X&dLE>{&fQw?FplF#lPbi`bx5@siBf9% zT}C}9dTBxEBu$NhlGoV7ztt3idLE>M(ka*)*yHIXO$&g1aun8dx=S$b0rr(qyap^P z!1N7ZpDl$xkwcLb%$IqIpomYu`}oo$Ec>#{|DlL3#8jLvgy43F03XE7h5ufl>e&D9 z4weSdy-DO3)l+u@F0cO*T|($7^kiOo{N!GG^GCBNi0I)>H|XG9{0xBbrZ05x4O8C1 zmJ{i4+$rf`D_-y*jEd6X3^%Sg3823O2w2veYX)sZVJt-O(6e`ogr^YEBdoAt-(;4G z~<}6@|SS;oh5ZN?XY@tK%FF-UX!{QsF_iB1)Fih{rF{sT=-f=4yjHW7zC!N?y5*%DI(A5iscPUWCUq2AU zlV%nN!>0L)FH)C*6^KhPE>xZl^37RQj2Nv@j96j!eWA*LUb1E6G&gv^ zkCDOR2E_Cf} zdNE92ESrJA;s$#A^q~w>Z@x7CR7g61_7Wou4$}FLFB(LWDPVk_&Ob{LVVKUJdY{fO zQs6p6XHWf3e-x__5hUqHXKzaVaSoDF>+p7iHFj|YpRqmYLy)n(hs>}O_|@tY!_Wk- zqRqeAi;Vgp7!$@R;QmvL$>vVOntg6Jy9bgL2wC|dAZxn6{pTUJ5mhOis?piCO#8zdV{vSrfPf~4tgS{f4YLo~#Yru6R% zgL2wx0`^w_J@IHrx7@2n8!q??69$6^|1M6cAblg=?6+T2^fZf%s1MmZhnF(=S~(HU9!z~s<*5i5+ZJc%Hc+O8ppNYXv8zum zwbR$6i%)l7SpjJ3WPvVJS*d)%AD9^~KNPS@aD8i8T zoF+a!a}|vr-mcb7>t9f@*Ll;01uNuKpLE$+WrJ;K6Cbb`IcR}-kSkSu-PZ$Vm^Zo? zXy&DBEF4GVCaF(rQy7O;9(LV@AauG4L@t03b9BIm5xEuV)BXK+bC3``>RuBt`sV6gl47b7%Pk0}x9_i+?{ls|8ibp@upgKC6X931D1x z80KL7Qo7}`(LG3~TOv&mnY~0%kAfjtLi8=>ncg$4T&)n%y$2z!)^L|hy16Iaa2nkb zO+r)$Ep0MU6OJas@2_EaTEhXd>CA`e1WSycc(EC@(g>EAKN;t4M|^Em(`>dgkbjB& z+GwhIwND_aHP|H=AK3X>5&f9-`w75fAGV3-Q=-`zlK_qECxaOl=WPN`m?R%h>%Z;m z1JfwJQ#9q`Zp|osGB{vK-+tZ})EW#>i2vw$#mrb9<(#BLZ62 zMt&XRUwp!6wUj7V`L7^{tUSp(-NTNa{kx^0*2t4$ZU@Wp^syq@F~c4umnNNcZT466 z8_+*=sc{1#*7H3AJeZ{P7_`y8z^w|iSoWSb!1Ch;>a2vCz@_ZuuTAw@Z-dx?@`&3k z3~vMx4@Kl-4&pQGM_S_AcAejL_y);3jb~by`)VJw2A-6oF^jih#)mAYGcb_2c2ws# z3jwDTn?ZJ!)3y&DC_WQ00#m!QfnHh>2@Ts~}aUX7Su)!vMdD*gUEc1R6oqPo_o@sJRc9eslg6z$E0QEow-!}*jAfA@HY2U}=mw!y zSp-);355Svv{uAEmKFF%QOZPC9&=ktPRdjQOXLN#1U|IW3e?gp1`(jM^vOFIC{()- z&B9fEk*JcFZ#^R_kMaKg`EA`^Ti*DigMnGKyS+M6)tAC5`8bMI%|2L#ssQhrLW7ip zfp2yCk@X}dti6K#J|>Y%Ye-b55RJ=rNT)pIwo-HqL7q8Xbm{xQaxA_xg?|Rc)ne(x z{}L9O)O2wJeeLs(8WgFHb-gpFCaLMYYmFAsXBk~fk{6mnbj?%=q{OvGMfD0#ek=lE zWwoMX--)Oo-yE@88WN*G%>wZl?@ShQEHbZXJbr+J%-c{(y&*@A{2YE#FrKOSk&pw%PdGc$xAS>VeM4I3o`r3w#HEhLPV+gc~rgJ`_XfZVY`*Cdr zWS2Rg|4MFeNmEidfvpJJ6D2rSJ^zE&`GmdY-aR~}rK$f!3wzf=Ol$`6!=dp?KZr~W9(LD z{wVr&2FdZCINDj=S{vbGKU0wB-fkAaE=sx9eDw1h#7>_7#MjPohU3)3EW&Q({!a(* zDBf)2drA72^Gg`umf}Z>T`lRz4u!<|98DG?aAV63^~L!l5>=AI>gdUz(B4a zk^nn`GJi%&@Etw*ohttswhc&j=q-MvcY7t}a7AGj-G@BXm$TFoB$o}MaGQdl%r-TY zX0}01sY)$FNrKO9>|3j3$z{k4t$a-Gt`-(7vXG^dKsZPAr277ksEG$ANcgL>{hKD6{$5i!8VYR%O<(aI{ z_=L|kHq>MnA+Ut+TK~dWtd)wVDb=Q)9?_oD!U8n4yT(KTtI1OukJ}O*9)8PHsy-p} zhn+_7npfnK9*^+Z*#?-m!!km+t`5{Gx&6P9=s zLa5OmwgbJo+pg!SG2aN#Bcu~Sss+1fb0if`K~+K@2U_mgY(G;~1`s<%@6P@wAlc}0 z?rHsBF+h8ic@D=GM2}wncD#3i=Eco;4_X#O8)qo=7XmWsP z{4`$M336FgmM@8}a^$5X%~3`yno7~D$Fo+fa8AOn$`(}l$ki2GdbT4&b4k;o4BZ_} z5`UVJu@1wCy{*d~5+2*g#kyVR@B_VwTQPy}IS2mc?Bp^(rwtH%RNA8@>^b0$uWFUU zal+(V&L=#(5DpJt@3gPK4{V86{Aa0vwPJRr?HGuUhI3~6yJZT8(nY;1os|_m1%Bt- z<2fnZs6x`}b?EK9s)PNcYJ_H7g@% z;+0GXg1rfFshwZtOsj+j%3cj3oHgnEM*Jax^xb6QG+RGO2-%Ri75?~Sir%-fUpbagAyYQ*f@#4EqosoevJ9jUQY`}nsj`9ddB28Dd!mcw zuFYoEYl??Jvy8m`O7u5A)g3{(?~$RGhnV5W2LtbyXQZ)CHZoz*vzJGIPCRq|fPvu) zJ2Zf8pUrnQ0wj(ekpl^901~Am-r=Me7JRF}YSv&h9!tPsVf8y98h=VK?X9^uQnRnT ze3Cd$xPO(DPaYLs$9q8ad%h>$L!?^+beK}M(|NUI-lP$x9G2)DJqQB`ry=!SHh8tq+XJ(O#lRKMk&$c< z9W=F7%0PtRC}>1ljK;G8nq=E?eUk+MW%7Q)w#dH%b0zC(1pv~R0jWi@wRZQC$WLNFNR zT=Xe6L93+@L2f4{g|D&sm+{gp2>bTy8vNue5QRCJWS=x2>-JqP{-p$nSATHkaa9V4 zGRo!>`MHl76)2<=k|2q+jZZ4bKFY31ctXi0yTN$mZ7#LZ)s1NE)q0 zFzX`h!AtdLrmcU88*lf14Po9PuFEgD|H!3l+?+YY23ez>Xh?Jbu8>!nfPY=TQzM9N znzD`lL|FN`wYAaMkF=7A$W21&h}Mzg|XNB(hyyV15twJD+S%~wPB$Fon+3mVOb z#%fgIjD6AN{GIT!9HW3h8iW+RFXm<;m~LK~bjxQQxRf&Sjrb|-VV3D3EcX$#ZhyX2 z2ys-vk2<&dI>}4Ay+l+@)GnZ^9bP}T_xr&!4?^WDR0Q3a{Pg}M#33-Rt_WSIP`dF_ zOc$#c;}jgFCkCXVz%$!6EGgv5LnlG^cP{sqP>xjsmk%I(VfLQlo?9zK+o`GN{Ui$T zDnqEBQ9r+1TWWQdXe&Utd#|3d;Sy!Kk7AF$q8M@d3yq%Ns2G$Ia%nYICJ#kPwpo z?~KTe+Dh`u3s`@yX=ZTqV+icetLAg?MGPu#LhWDS%I4;mqd>qfrYhN)WtcaZg*ll!Izw#TY!+Umj+Ru!}mF#~FioVC4tmR?xKI5Oa-=fH8)qkK< zNXLg9`1Z2X*0=9VGLWd!)z|P*{@|jm{~1^pTg%~AUA<{4Vgq-uy~cM-SOkC2*Qr{f z0&qJk(C7~c=HM~)=>EiQ!Z0$Uz8pf&gJV5*z8pt*YOLB^{F9DSI9q()**OLNplA4& zcuQ8o7H8M&p9b8t68$9tyoF>+o$24^^o}B#*4o!(I`ArYfZiPgk4qp?w+SMB*iUU? zr)*&!ZlfYzxiy{epC#t=kq8#lOG)W2CSGyi(YK$en{kf@^1*!9ZftHOO%6=PzHn+= z5zJx+2;khG;ii4X;KiFPGnUCzaBhL-!df{_0sr;GTC^%yskS&XGy<$_{4E@^DT&Z< zu0G~WcR_PTe5M1!6qaI-gp8HLj~>jex3=aYO&3Smc;P>IRtDHyI7+xVO;hQAdqftt zQMQkGP9?nDo{Z0a!l3N^LARg+&C7N9-X5Rhn><2$&a4Av5h4GPEDleZ8`LBeNFPog zpw{cRg$cdTb!}aH#$M_5zMp>RM(jR)%YQm`0EdL;08fljPe!UncUhRjO~jsNwS-6- zYuLEzLLF3xeOCskLs`wUH-j160>~S=^O&fqTrf_wuM?NVp5UJ^7O#u*D8hRW=JWkg zS-@LNf>;@gg~8{#lDSsg*}sZ zrzxDzOO>^xgO0g^aP4QwG4vInZR+y<7{a=Of-;jkbz`CJ)xoKx{KaZkR!Q_knDbTi z49_UxN*RDJF|`L?aA@4?mO0{Ghib+4!51#yp)SiJ6tF@k`m8>FCH2rAqWz2%0a~BP zw=g`ydvh63cp{)jYY@tptaT!d0#%9n^G)Flo~s10uzy3?Ufv27af^T+&E;IUGOyx= zv`8yoMS`Nv7gddjJJ1K)Rf-lt6%>Q|S>>a4e6-^4r@P2UTu8+e)UF*qfWe9s8|Nd2 zHPxoipAFh+1}&X`1TN}-^y_Hq{tq;OEkJ%Wd{|wBN`k;2Ftrj!q-~ntB!~^hqoR)TIi*L8U-iR zR|B~*4cOC{6xOB9Ig?g3Ihx1}e%c({D<&c6jWpb!Qmdiry{stVYyt>b1+#=I(oACE->o8yQ8xkZlB z0AW?>`5U>ELA%sS`4TbSU0*$WN&N47?rXvzJt7%{tPtz1aVk8$kcnA353KYVdPo$p zg7iw0dTR~$c}$~(_9^oXWmIEwZ%j~7DzMi4Dp?U!y3ylHUtmM<7#iT$JbU*EO}W!U z{8XAL>epelwsjqlga6!Erj0hrf2m-)`x1UIZQwVp0pb@++mxb!SBa?Pj;mTOI`TpO z=*T|{?08JaezId?D1%WppHAESF-s(3_b;){5K3qqk^{IQ(SHIo5i=WnXOW@5BD8xyak6DDX!e14>y^_alAh2zv8ZnWi>mdAq% z=y+NtMZWtkoKZSw1|q%Uc+y;|G~N2Pf>8w*xr@3Q@W4QOf7jA|>E9zWJ2W~((!Z~+ zTRuNC<7Kp0ro8{Qm3~L2P7escm(EzYcQQ;z0OE5OCucZEr-A9zNTWW+IYvb4=`%zb zq;44FU7<_gK&G<^{cwRvy&2y!XQblhkSfO7;_>LF-SBSPBA~7_x(sP}NbR-tYHMsb z1Ic}kMk5T;Zmo)g22_Y>^};i&h05X3M8j(giC5Z(4%3P?6N}#NyHa;$79>{#*iY>k znut*E)iQ_{@<6AAYi)g9EcE1ckKQLw#?1oBcGn;fjhZ`f1#tnm@;<&(zT}8^BHp=` zT7^6Ct&98m8a_fnUx-9;PwE3_ziX&iNh0xyQ*^Q%&|vp^d&($}_HZI4>f?F<+7nTR zNoGB=6t=}wJCy~n^54d3D+MR8QhX?AU@!R)q&NSk#6?=W-J0K3$`v2emg?0dhdu!nT`Y7%!ibR2Q@^CJ2jhfG)ZOCl zveRA+kp&0DG$?sS2dyY0GEkN-#$y*eB!|S@;Er>ga1Uq9>&jH1M#pq1n*%-4 z%@+Fa;Gbv+lBTtpH*YBLx1$hUq5Y2^7yBr$Z5+{9Ae~}fHC#vq5YaJTXJ#bIP2=w& zwiHN2&i1P}S13re5Y4u+eM%$__7t}3=CH60OCj7-e|Y#+Qn~XkmGH1(F?tf`Zd~{s zGjcvjf)I}XwEH)ja1e~Ydt*+3WWrjv8dqH{fne|=3#%t5(g$Xg<{{3dW@WpV_K>v*znmDh0cJw{$9Uz~(mGkMM z;{EC2I0+HjAS*&Av`PK)i}`2ELmjx`TxPMeGZoyxIqNoZ7!19Ca@P8;gZHplUaQd} zh8yB5me`hv6b5PC7(#u5gh9wUGZYhLFT@PI>LUm-Lc~tkKG$@CU8m-SNW1?IxQ<~7 z7)Apwp5?@TJpAFi$T!y8ZrnnBA^2&)M`Sd>Z2@R4u|<49*7LoU%s~lpVJDVXszm}B z^$)KPdO&SCMuo<6?OY#r@~R%Q!RoLD;IU%HC;VLeG}czoC&INlOsD5v=K}lL%D+?% zH;BqqZl+QMuCMGkyMkyevSBVCV}Du*s(iFq%5?%??0YqFb-;~KtY9V&Ns!R2N^wk+Fo&r;{t*k|Iw+^0W zvsJp{&KQ+Os^1l`gbI&k-bbwTz|7C+KnAz_9f z{2olfT9V5I`14$;8ucd#{JEbsX328)z!vrkZu;|p5tr{b zSLrN6W1JY9;i6{7;yR$=#vR=f7K?*OqCFzU2 z`6~TgnXmcf4drv$sHVWG^HSfI$n4kENPJhJ$YsE4j>6h6==TgY)%aX&Qv`K&6TbjV zanPfzNWUWDOa<|rdyE~m43&d+DlYr45(ZyP2Al|% z91;4-lc?O7T;YrNM%;gzIzUy&q_T&Udt6wBY`Xm?f-G)rigX(2fJSY`=JbhV1~qkA zUpiK|Dtt7urrUXN73#iSUvkzq3TyT9tKziy@8;Jp1uY;hO%>~{2!X%YeDy61Glwg@*0wH=E?CxG7d-Dh{d^@pu% z&(1x?Dh$*mNxDB5(z%J9@k4-qJRQyx8rA*nxEtSyYZ}Pw&|LVo8Z3gpqP8PV_}`#Q zLP?rmHtP~Txp_KqF6TF5ZY0s!^#l9!CC-gT=h`RHV2^&m$q*t>fN~Y%28?(UwzH;g-Y0ew_$(%=WN ze&dtw(2T!n5B9ehZrjFgsvN3_r0Ut&D9^Az-1ag%supE~>6?A5Iq%^}t#93I(JeT2 zc(Sz~Pzu&b9jI3O00k65^06(4Vb{TRe`B2I!H*U*qG(!(2QsaNT+{W?p2xdi`Gncw z*K-sm6IU-_f15o~P{Tb`yvpX_QREd$xmCqmQ;aiBUF}T2RfKG-EJRp57Y!AsJ8;)q z5OdS()2hMar7V!MDiq3o3e&cY*ln@*a&zdhmc0~3P)nUm-RB$ z6*Rp-s-UvnY){k)pa1lZpUBjYT(R-(Vp0;{ZG7b-C~wKwE#Xc~XF(TJ4W8Hc!O}G{ zGq`mRK?Ko_ckNxu*PBszmbAa#cmqMLIk&5pa~T}Z?Q))y(F;H$-=CW>;b}9?IUK6ZTK4|@6ZAg~otv&L9LlIXr$Mc@@ZjtL`z-W>fj<9whxgu0b2F3C zdUj=$XQO&PlHI59;G?QB+d0mGxBud6vexMvAs$z+)^a2(q#eU?e$+_M8*8QW>g)rO z><4~DQ>=}}9|F>cBZ$puKC3_TixHk|wySbNxX|yG;+2z6%LYP{oc+cuFdLWkf=&BV z>=eGyxn`cj-Mtid6Si%C=AS*^uD*p8-PE@B%GA1Q?%~d{!xi#Ta;9$~ z9pW87=1Mz$T>pln9r_u~Zf~kD=m*|^j|=Gds#IY=jligU2V3MZU3gQ08< zm&9==?iK)o^WLeq)#^20?sRs^4nkd7EU3|n32XSvQWmb$8JQoaX!(ltNoZ|Mf2Clx{4><`n2OS;9MMh=MT)x z@<4ZLuqgyFn4`8m|NSZHcZ$4DtkkBSPzd%}Gd7et$7xL7SoZp$f#zhA|4EiQ1jz_7 z^8{DA0jy__Y%;lu?)3C}AlRQ6yx5YmO1cNzlgEw{8nv z(4U28+kkq7p^DpQ@ z?|k?fSI?28V9llvz1G8pS!9|3rDk}0!Eu0ZfX!|72YUNxL{?p;!A1|dxj<0#!SYXb zaDr$IA!#=aM__ewHHRFy`O9Kf=d$xaGq1eYaY3rNKY_RBSy%X;V5D4w!X(N>(Cdy; zGl99tDwpn9GcgDK2u-@FF~YaA{vR{^D2@D`-C_;pXdJBoYE#$)&zg6zELnlh@mVa@ zr9=`g@nD`wh1WX#7xq{VD{q^*`^}#v(p*XNCRX`N>2_)Z4$6=3rTNzO2Nx@IhMeUon*6VX z$ouWpkSnuF)a%SPL)dHV=6$bja$Pf;AIxm70tkY+a@4y`pS+{Uaz+`2bG^nQJNpAj zUNd-d*2R5;0JiA+^f)B3>NF_(82v7NlYI#LUnz+NJmG%ci)h)8dxB5K8ZysY{lZok zBVD#)Ng`Lb^T;@%ST^|!vg1O#VXLL+w+%J+=7Nt2=9QOR9}<^RXR{zUHj;s&%!jg( zxviJzHjTtNM_Zo@U-qeUSglBq4x%H;AitAcv4z5qVrOEVwH^m@7ApQt&Sz&@~h{31j$K>*%$9*Qrd7X_EE z*=~Q>l>ETTjv4Hmh{MDvNgi;n3bwKhj_!%0$f+$5>KoY>&e_e|Vw++b|K6I|E3}Tm zgil)4O)Y^!SnQn|YmF2hf66Tt7@$00<>jN6mn)G|sea8}hlfgdDMxx1Sc|+ z+aKnBjR;w|Uo(G^<+b>tObEzE?lJTT!UALy^>(?JlP4w0Zho5-h99Zm0coADW_ zmYr9F+S@PYlc6h>UgF3}pW8w1lU(p+Pp6)0N=MvLGqz1_MT9)5d3XeJQ-CesRJEv2F$ZR;l)9<1M*@f_ z2yv_qAkojOYtEJpf$bg!?Xh#0!HR}vwnK}Z<_ub+TDpQLLuZ_G=UKYo3)Rw7Y7>j+ zs$f-_9(ceO`Q7-54g=sDoBx`rDQXPr?k~OjxUg*-?C&vC&gKs_nihs~U4K)T~y;l~!KcC;RN)v5~A^fH}%5+7*8v^xZgF(aDgWh9M8t?=` zB4tzBISDO<=s!h{Ps61u9J>^VI`MMJ1`F}gFS-3pO?LQ+f=z$oqSZLCEWQvV;Gu9r-4o~5h#SxFXhOXNWr|rSS>&|Fm7=DJ}&)5t@{?>pVla25n^e_Cfh>pn$ICZg>amq?_~c8D2{S{t&ob< zxL(l|$Q({gIMA`Wqr)72=8?YTdz$&tljy(@@&+NHxv3YGHlP7^*k) zulxG+Gt1jgQ_=PH-hX3X1oVB^b(f9*dHJ0=`5wqwXASichALbcE?P3BOd9FJSJ;ygAY+Q*zCbpv~dElElyfsiO z-L>a=|MKz?YH9OsYFq65z#;Df+3VjxX1O#?OG7n+wL1DjFRoKtL;yA$O~Z$OO(CHO zS}#YT5iIUKoZ(~CxeP*TPP;-#iy8+JYAIbNR&xjkW zUqyJ>PKCB6lj1)Oy81r`Ha#O6Vc6POGYng4z5zDpzCsl4;hM)?!hdffq$H4BzqDBFk5)PddPTW?|Vxf)+gUP*SA3AwtN0*@T8!# z@*m(6Yzx-dn<*BySEk*kvwm?Kt9R~SlOH8AI*LM}*Pu}+4C{}O;J3P^5(DCj7al!S z%dfky{e)#?UqQo_FL)|!FWHt_Z;*KGOUiM`^JDBk8e@4DNj4D?hq$6_K+rt2d-ZuN zpSA+)@V4@6PJ9&^;(4@6;&|LArXS*#(5#H0CWy$jkbA^CKv+2~NOe!0%l1zp8FkyE zKS@TQnC+o6Wyxf1CB)Q*)vDCB+I3P)#wYFdK@vlJ=c*kQiLbO{3OoKAklC_~Ufhnk z7~Z-u_32c_quAs99`oT8>9Tx*)t&6gP!nD0>$xG4jKW{rT)7xj(gZO{NrZL>)om$$ ztYZ4)RP?`=R>RFKkP;Q2&4R&hM5Ss$f`1l%s!KKM+ncKl2a%H&mTuj`h&^}3gN zJ>OV~wy1Q0kcLjC6HzYL~j1cBH^f*#T`UOf-Vq?Zibc?F$Dh>E)>yaZ@$*GcOQF<34Q`#;!{`pVOlsnbld)3_ z5z(uXRCwj#Bv9^Gh~7aC5cyZAp4;wu($$Q-^@a073Ht#Qi`O}WI%CO~dp*9|5 zZ*awojcH@^{l}}a(d%k}n`?Q<90f3Vw6eS$xJf3=@(+-nLbyT0`sF($e_d+bL8gr{ zgqv=SMVyOv8NqKVby?bvQ}MB$B(C`e{+4!pecvJpw_W{=d4N^*(MTnIbQpz(A_{nw zgJycoOCV?aUmrY1h1Of{5uF3MphhNxs}bqF0ZBvF#b@jWMwIa#idkcmouZ*`jF9pB zhK0tWCujWfp!YmgyFVe@Qi-q^ZuMOe&HAx}kcBP;Ho*GsT9M<4IOHN^F#8aNq_da< zag;c;84cI)D+{FmFU3Qiy0r#Kmle`~NAd8__0JXVQBcShDNi{G?8Rg=o3uGi@*5#i z>wbE>j!OEcI?=E@P&$`Gq}-U9^d5vw=Z&dCizT!`JGPQws$?|QF1>#iNrvIKhU!JA zs^^~=Y8pC$Lzz_ES<)of^f$)-c`_t4Jns+LTxw`>5a6!C`N;Mdl!i4FduN{ZN%5H~ z>%&dawHe!Xu`W8jTdMV&7W%t;)`iK)wM(npsPjc9qHk14K>lIi#)}Vgoc_f5Yd(Vk zX<8}*+W6mG($bgRPXyS?YJ(@@MH52f3Jdhu%+8z4qPik?fV2sr}ZOj{!>09*7`tMhzf|pyjsZdSh z7gXx{ZAmoV8bvyY5sexo1jCC?t4e9>(GbCza&F% zg0yE4Fkw>?V#O&^cO=og@ffC;SMDCJREDd*O-2EO$v@Bxx@(dt%_E@vQ>F`0^f`Nn NS248V%lCzY`yWLXC=>ty literal 14969 zcmW-nQ+Op!+l6D>HYT=hd%}rrI}>MO+nU%mc5J?}jmeJfZ~uQEbltU{)vKO9=!3et zNMqo@{{Mo!^oP~KQrLJ*!f_#YH~sIf;z9vLkh%bZ{Pj%k<@UT@lNTbig=@Pn5_O@URDq z{H{B!1gjWrXgOS&7$C7PPA{RQWpg=@?aqGfzoaAONOBRqT)MF}kXqvf^AMk%-%-7>qq5^~ zz@G2t6@0$6v!-s)o}VjZA0L1`y7`nDh5AARtmLgpm_()1mkHZreQ7` zGxq67=Ked(CJieS5a%M%e&o5V=bB=XKqlGcuoC(bDpkHgO1OVbPFmDZF$l)>?qR7y zvcAGrR*e}mC?6CC^YH$7gr`|?>gN@M?28XaqqgMS#j^&i{x@G;D?yqKZ(t9Db z?7iBqSLgNqH70D?^YOd*%b`(IfxdI_>bF@W;B|*X({kWr&$ROEXE!a0-tF(cX8gv% z<26t3h|@=d?ePQ1Q^Ns2!T#&hGpFRlHTbOT!$4QLKk9 z{(Dl5!mGq;+4r!uqF0GcV-NQ1F}zGk6RAr^4f5=9yyW@0=u}}YDqz@;B9&Gr>`U&x z=K;*DkdOWN;lB3>B$MGdKnF>esj4}xqZN?+w{M=Km2TUOi6P28vw34n_plC1g_!kT zxTH$+SkR)lhR=1X8wW1T`15~O`|VH-euv7CyPFW9h3fk3KL7YV>iq(oNj!%4g(_fU z7IreJX+^pYn!_V4;KM~SVpZZdCJbT0<1xzS@`Nj~6y+IJX>Toy zcd#UPHzV?3S=q|qKK|iuCg`!+gp1|NTIV}%k@p|C5T22c&&CM_Vre1lVNJREI~F5X zFiC8CJ@5o_k70Csl_>ANd&L}}9ot{u?oOQVqz`SxhO$1@@rd_eD)>}(*2#BS6J{7Y zwS(yE<@)3y&N6CMa0FN2FH=&4dYGzH6Y)r^2%^?{d-H~eEYY)}j!zoaCf|n<#h38g zcP4!i*SA*98|1^6P2tZjiS~p6p4d~TPZ*c_D{XO+QIG!m0lU{K<+HJut6f{(o%<}# z2eaH~FM_@t11AJXylY5i>B_Ne2otcw0OL3{zxz!^JirU7x1Ayd7K{rK(gwlvizi8_ z7!hiyf>PPsWo4Kc5pL)s4J7Tm0Qw(H(#Yb5s{`CR2+93ncnVq%=z5G|=L2Xpni^V5 zR~8V$sA<&}cxjHqsPW#0{;*ozjYRK^`;p5lp1Mj1wu9FtyaUG})tIhEMnBhELSKgG?UDY1oLPB`yzZ0`G3MtK zIxg)O8qc+~X_*ij4=oBhPbX5WRb&okLS&9)?+MiKOq9Z~kl%kBP8+=-wII^(xy)+B z7aGLnk?-W|jTi?f!wd&k+eH4G&MXJ}!%PPz5fR8BU^%$m*kU#IB)oM#z`IqXfn9UZ~oVIig>9 zKZ0Ir+B|ikaL9yFRLDXgr%A(pkPD+2k@bZ*@L>)OTPZg4Racc{%ozPUk+hoRrg(BZ2F1!8;*s&!iz1`rWo*zagkcG(kvm=L452@)RO z#Tx})_c~s@dWDi4WtY5&%|3{N9Oat4Xr%k1c_kerev=m+9u4||R8)$h4WVKDxg@j* z)Ff6vaasx!6c!4TWXrZkauXC0ihWks!k8r12wHTPTDhg|N!iSy3|ViM*Y$`nzMiNs zhUY#nVdI!EzK^J|^)N5kC}(f1sD^Srj40<9tSGcydE&qbN#34lNzfMy6yYNV$=Fd1 z`-~xyyn->3j+X^n&YqAM`!`G&6k+&80=Z&~!VqX5b57`@KM&ZBw4AlZs)!*yQlYR44BHj?^FxkIRFNq7!YW@kx#tt3I)&S=z{+uR(_*6Fb?9uvC z0A%-|Vs_S9sy5-rIK2_Tu==kH4g`W2s3JXr&yO%vv-^Hj%CX_bS|IYLQdYLQzz0zv zm7xdK_-QCD#_16_)lQ(f$5i>LY)A@fA{1d8W)c)F3mN?GUn&g>c;x_542^i{+%v;i zDqbH5moqVZIp;g*WHRa9HwV&WfR987r79ZsDJs1mzzH;>)92tIfksSvF@i3)W_0XR z?5b0aN|x7C5aQA^3_eW$4v-x23D16Hv>t&Vl^l5u@A_ft`Y>n;QX}N0hl=?@TErwt zA5QUspk(hC@O+uNGBWxe{n1BqZnBgX6i)rG{4|FvE@FVKxG~|&yG2nWZpTse1G(@4QEujgWt4r(`8>Yd{j~Q`A2ri&Sr^Jvw_LI zYTAzW>ro}nB&isI3qOC}q08<6{xMQ!f=Kn}zBrXy&b>JiYc@gMaG}n>bTtWUHbuQ= z!-dBynvyP;r8IKgsFwjsZoh}FHIgcGMi${zg^pe5$|n?*mu5ON0tO1Bm#Uep{agpt zH0fM@QF3L}dHdyCN4&>K4?s94f4_o0{w%i0yDN`by&PDACwJJxpcH<)z)ux-U*BRj z9m+Mw6Qm_BKqCiY$%L3Nk&owuuQToelY+&PO#OL(!7Brk#95~jw;k6csee_`TB~HD)cdU8qpJ07Ulvft1ZLl@H&PPzXsQw zF9dADbel{e{BdYT7Uu$1H<_EXSFA}#W}KfpSyr-1M`b*`36m{(q)Xh>^l$Y@%00-1 zxW2`9=w0=i#JGarBo&89vh+fxTf<}i+<-|2a!RBdyHN}$GOmy|g>{gT#hEvvBvFs2 z2B9U7_)4Zzo+S~jP`hH(%a)Mw)u4v{HBWVh#`7+ zT}~ugVR2o*J({HCt1&x5$fip9g+jqsV|CQl5QQ^gVn8XADEKYlM-}lHw>1+sJ^6eO zO6PMme+lS_IBM8QqMXmYlS7#>*``!M8ZBdjnK1jLR9Q5CSEfqvNEFN};}@pzsfjp( zPsL;iJ!)i5M)fDofAhgf9($9`@A$FFvYI5gB6rtmUo3g2!!m(+7WlyhHDR_*)gPfG zDob)f`+6VsY2)vdD7d0_m#`Kbchs1lOj|5&O-49jWO@2CuMVa782o$fUGS^Q({?x@Hw?{0ze{1gpVOgRHqBj6famrj(@b{!gbQ(j&JR z;};&eRLv>YL`BM52v%dAhPy4fOf*(wp9T_PKc|AHu=+5+3R@PFY$C;qU$i+V{W!?} zjx_iPohlJyB}U#pGmL_|hA>Tb8L$=Eu zd4o?D{zATVq*lYtof(wClD5`D8&Q@A)M%s)?HFX+ z_7TyU^;Ri1g^p;H#ClY-lXsQWE+w(A)K893p~p&0>WtB+47Kour8ee+!iL%i>C6nO z6`|hG)PhnQ`!p?wOu*Gg`L$m+{X&YaDY0w3G-bJ65@qOho~s&$m~kqNKdVh8lpmN& zLOdD}$YCtwwkJW7D7_9#bt3Q( zKl*2bqgLl*lJPtfsAh%DZs5lQTIH`Mzwq1K1D7Pdj!JcM%~3|^&&E0}YKx?A=#uEi zR{CF02xOWy&6;^47||=Zm8QM5@=cmI1SQcP?eu)PiM}9H`@dXHg+HtgAg#mpmNsO6 zF_C`duuJo*-?h6v}idxPXXD zd@#6`W`TvulhY4`vc&ID9-i$Opb^?5S*Te&iT?@t-@g5U=)X0Qij=9G?bxNp1l58e z6xova&oKUpDh;e25fCwSuG zYAD`Wg6zILK_Yfu3p6Sr8u2>=zkt#S zCI6GIJGZ?%XwKBNjW1rbE9DOCKC3H1%MLqR3taegfG%MB5T!TdJH$C>-dQo@a?wP%= zP59Z`*n(WPzgFDY+T4OlK6R!yOtG1fOurVo=wld5VD^5uikoH!IhW5Kz0K`~!=vhTN%4jd#b%s{4I|@YY1wr4yR^vl9Q{t;w)Q?s0k7fDULt)U;sQyqrO7csL+SZc%PFV0gGMD@z5Yht$;c zds-hQ->%%9oaSj^+QXczR!R|Ao@Kvj9+A1G%^`=suo;xOX2cc-NL2GD{IIhYs7Y z>U&cU52Kwy57%&&l;<7{4sK^96^GRv8dO(5qFirJKJS(72z$r@b(NmkINn4G0qcn1 z)VT{O7A#}gb)rsYZS}aGVnL6#bSPGN5glgCV+(oBAppnR7Hz*Kgcw5T6zJ@8h)VvJ zknXUg4@u5vvbrVdO?EW6)~iSSPL^8GzWt_88s@vvRtvBGXa0CmU=Ha>>}*IASH`R7 z4{Ieu4v@c=oYC0nX6DTXt!?r2BR>%x&GXtkn(i1W46!|b-%bZ*pPAIjLK|811^DaL zB8OJurE)9J0vnOMN9opGLKRBVx>lRHClH>Z_yu;c@}1mN%JeQdehkf}bjx3`7lWb? z4pu9A`;Y7;+@=3#5*ofQn(IXzxohwZE?+lIB;%Cc4ey1JBzHja0;_Y#k0{+7ubmB4 zwYH6=7Wii}`rF>trWfHzy|0@9VihdOJm37ZO(Do)!)THxApw+P=&`_apn(ht!RwQDz zE`f-A1ZW%^3O}XJ*@t5fm*2((q1f)B?|G=MgEJS$1zx$T$xhfyt^ZA1MpnM_tC@LH z4?lEnPo1bl=yZCT1)rmmt#n*mJTx+ptuUA?$Ei~c$7`zEWS_t`HeX6|cu$gf&Mc46 zhF(hmi_F?rX)MKn$Vmy8&?acV{EW&ReK>|aQQnYRx*=ko1B{OregTr<3@p>Mz!*n7 zf@`OagGTgwle~P<8XkIw1Qe%(EL?2LVq>5i3jfxgtdlit5%8?&qk*S?Plm082mqoB z&-A?9A?;%KAH8WcjFk$fW>cs(6zY9RP1TWWDh;bvgmQw@y2)s7Od%R@7(4o|{t!2$ z`)G5Eh#=P#TSfg8P%?I4ZFpxDLDl&sl*(af32v7$(AiuS0@(O2R`ys@GA2gm<}geK zuSv0M)FmI+^p^FcLMH=!mwI_+ z&2by@1qshvldQkj)%)Z}DU6rZRQ?(Gq1I$XwH{vEgK7Pze!i{;0kF;A-S<5qXwVC- zh#7r?xcIapNs=H9dbAP$*{HzVu#rHD$JvckGh6jaQ!f?#`05$HefSFa$9(l{RtA1i zJ+uF=Vgz3O9H#i=Pe+}TI4mZy;me&$G4e>w_AItZKJCl>CxJn zv(by7lkx1Bl1zE>36U^m+Q@92=Q6QlZYaJ-E@UKk&d?Cg#CvbEU}#LN$E)T9?y}Q2 zyJSzgB7{4bJM=UnZ-YWY3a ziTp!oQzN+C0A^2kIq((z(CR_VNS|!%+||~}Xuv(xIUY`Md(F`u4Y{sysZ<87;WDMO z2)(L~4SWNJR9z9eU@xmpWhQPH3Z%ntFkUnR z)SMYauymTz()wg(qIiXQ8d<~f+UL8ZHtnJUyhq10nD0dGom6?|gGZj9ncT8{!$@JYE>@r<#sP5or&NP< z*!kRWt(Fh5cr#~^x{Wf@QNPn+lYl3(jx#{_k)O13Ll z09CS{yMI|%;pV%i6+7MX7)_{?iBne%qIhiM^2?PTx>&-}g^4~8VxewK!QFg{Y@>KT zT5HB4a^+;lVO|?);SNgY2t~XlXtocM5=dV4YgdidTu74Z0xNZH7g|jvjd^P&DY$~L zXL9m6Lu@IkY%I;g7OnL@cNS{{Ba+_m--56hQSk|^71Oer_|G=reef3!MZ&fK-{=*d zB=qUV`U$^yYRpc$4gUn1g{O+=zhN)IA^DTUh|H&S7fCeNadUM(ZB z3@fmwJCM+guJ1J4c*TZ)EzI~#%b=-U7U?kUK%w*8DcH`;iVsIMtNkO4ml%$Jn*Q$^ zMReE}ISP<=hn7T&C0?)jCvuvv(h0ws+;nJ6?rZT4>pAEDBDRq8zaV7rpj=5hyY9an|U(QA|p^FyK-F-A0#|+1k z;Z>#G7pHtaYD4`-CRZzpi8*>tAxmE>7fYD5&@PiY3uKYuiD!Y-l~{#C5#>L+P}SaegzJPjv?{Dc z8%L#3)uMZz>#bf`a-+vld%KLF#B`4Y=o}ge+1AV~)qJY3t#P9?y;Uf3TDp&nt}1f# zQ($er77`rY$wYeFLP(-d%xaG-O}vVhs{5;AOG<-``+x>29DGxOUy~JW7_VmiT3$b* zxT0IG$RDC0w6et{o+%p9$fe_SNF)KcedUYg{1S)(;s=&c>?gF3rnmo?@dM!kbe6E~ zlfecu#LfT2Y$kyg0mzhV1Lz1%HoQsg52-QaZ5l7zJKp-+tm5Q^bJr6e4L!c?5cbIB!fr~jT zG5s>S6H^U+!~56c4(Dr3Q;=2yCWh;O#jF_!Cl>KvdfRCIpTgt=mutjma_f<_M(ut& zan-Ag1K^@QoIKA#0J|Yn$ibpAGxG+j*p~%vCHN}&c+#);3PmocEGU1bBT{Iw*ih@S zR0OzW99-W0e*6&croqEzZZYFi1xX&H2tZlk#mZ*9nDKJ&95&^X#2(k~KN&AYvIgM! z&p)LhXdN7l1Ls37aEwg1{x_n9|vW^h0LvRvdI9-^Us?=H?+cQ7jxfJSU6RN zjWBy83J{$;)9R_ms zn(mQ$Gtv5g5#D%ZUAO=pNnU<+J=TA!!hUtRN?G6g-e~)>czAEo-{AIg{cUOm|A0Xe z&Z=~P3veLJfiNbiVP`a^wStU^1e{z3{U@+eS;>v{2s-SXV_ zMf(7rlw8d0L5hTiTNlYB3?y9ff9db4-=qi|_GHe#2YmQ4QDe>k+J*L-io|(egiUUH z-2A}$tSswLUFN`dll^1dAM6Iab~Y#5gaY;i7Hrc$9zORj+ROH&YB68*Zk{9h9UzR3 zUzdFbWd<@73=NpLk$&&w*3_(Ez!r2L_7kvI4V+f&Os@Y<7_gdYpMgZ>4xO9!J)zEF z*u(seG^8O!6~0PXz2=>V5e|3KULy#AA?0()7Y<;$`8rvulfsQ^1>R8~)ep4UzYDr|@&GC5kt(zxh*}39@(XJMS$-HCoSIYO7Pv z2ZZ6_(5OmP4H(OA-W!E){GP=`j>qYZ{f+umOF!lggyC;uXox@m4crgq?y*gIhMt;@ zre9t-4BXovcjPvr4m7t|@0Q%?+~~7=LVRl+yfe1fULb4_a7?bSl-k02^KV)acX5yk z%D_Zg3Qj|*WZAD%(&hjTH+L-hej?sjzok;==)#g<&uXjeEJGwT^>UEbX+R~+w8YJ% z=k2itX1wX-5$+{f&h=7)UL~y0)AuqDnZ(*GboVXbW?re{-#a5BiTdBA9jqR(XW^MU z!*bNUosu?<)M}V5h882iKswYbTB2g}a^KMR^e9ddO_&arts?3oUYM=_M5qn^ZNXe6 zUCn5Fb%%D}-^MjE{o8Xi=6%*`5bl`DsJAj~fbgo@J8d@84%rl7F1KQj+SBrQ>U#^1 z+w%nyDpgJJ@Em-xzPIt0U~a%;bwKK`15VHd zeP_faUGE;h-FgxW^Iw2R2$7p^Jz)n7X<;jn9$GsBpdMEPCILUq0R-f07GQKb^Bz zBVY=WS%d4Zt4wE*!hw0!T#>-ufsxb#@D2>uPiE$E=6OKsdXTd^r5qr&(=4DekE`#U z9ka{u7~=M;C{Q5-ig(p_!Om?g^xi=?TFZxbdcTi=*Qvfo;LI9omcfzuXIeJo zEXVcCx{)73E-fd2?@=eztVvI8y`x7^bIeFdIY)hHv)ofJt=%F(_fs$M^!|O1aTM5= z#YLFwt4$ig5x&;g`6T6>48d}hYN&GR zuMxAWCV)-$CIdTvO!g{hZ9Xv%-5D-}~QfQhg%;N>#Io}`AK(M}42 z^1;<$^OZm3n3K=F9c<;^nlFdx+lr)fjauI8x*sg>&ba^k^5 zj?T_MQjHYO#9__K-35-doD1DxnYN@KWO}{c04re5g+Rk!Gwryln|E56E(C?)AH?-K z|AI)S`{^qu5zkB;4IJ9x3x~#_Xgf$?eyUe%mA;|P%V*K0FKB=r|A`^#I7a#5@03=5 zMACe3AR1pliQ$Ghr)!bb!9BUoHGVD;_ z3`&FSvgPD8%OjH3TTJqi76yC0*DO|0_FFq0&z~e1jv5_{h@79Vy#jrbu2L^kM$*Zx}4vkX07@F{5ACxz-!fA zf0hg4i=EueW?h-Uk*21IHt-gXb-hS`^0@^6D{W|`hm{5Mw*uD8XE6p)cVRqBe}U6b z-4`&l@`bdPcS*9`c}}*4c-cnx5+=ts=VfU6)7Z_h z-6TUGNU1Uj@jj1z4a@tqC|j~2E0ssZXMnH)zB8sO)6|$lzTe&E+@!(e>W7@o1T3Km zN1PI_*+5rq_!ToA%!*3Pi4*f2#KpZ?9;qv3!_NL8^+h@Iib79%D$N}9h2l!aU5Je_ zaijM}cBF*yrL(qk`ak)G6^oi544cR+OUHtUDswQ^ByxcnK2gS+{Z=8ptGEqa=(d4& z=JYzm>RmxKEWw&Cr-*gS3;rbvp#)!zs^2g_L(+iC8`7VL8b{uAEFy|8{3)`5l z_o$&LSxAS5k*i;(6&$n~YW8F?4l0#8ch{2pp;NW1z0O&S&iH2-Qi=#_px7){vg8%J z!8>4c;ar_kN8imqDW@D72ex3x?}z~!$Ic1?8O9Mjh7YHip4u@~ioHHtPn}B#s=rRy z#}`@$aL@E#)z*Ocvqw+&pjboVeiC*zH_0^Yv?tKe7iaoc*)0Rs4(z>5+QwJs+`tFX z!<)iJ$ZK6pN~&Xbu#d_7lE&bjvG8F|VHy=~W`Y2nug*`RX~+|$Mom)alfR?uPi5F4 zoyw`*HNu$FrI=c|8#>^YkexjU4_?MkPsVDKXQ#$m%a)7&v>4OO{32(&0r=BTmkGU; zRS1@LEhIz+#ZZ=R@l(-i1OP@Z*%{w4G1i#9`t2Ha!fE|0{_E){D2_w*Q@b`yfTG9x z=k&CJF-^-$)Pk5c>u1k^S>3z>tG|zO>*g$&BEtr;%|&{UP3Q5I2Yofb{S<+=F4dZq z33yo4vJ`R-%(Y_f{x`sx)h`&>N01qQsVAhVjiYElQM_|w1f_xLaaC+C!m9Z$Nd2Q9 zPEpDh^>IuQA|=0^Fj=5=;OMbB2~kEFCPiMpLCFq3Q=bznlj1juB1%ctFO?p6nhO@; z%uihmMI82PD90ZVDHC?b^|1;NDc;Kl_2n=@QM)TY6ogN?#BV=NE`s%g2Q{M&8JAr^@n zI}26C*Te6br)nZ-VKBi{dn)QpgDBd6s2E_afBuE##KMY@HkX*Y z4=Y`Fk*^!p8I&{EbI0g!E`kiIiTz?+zj_NP`^I&BVQEkXh`Xf9hEi;Cdw_!=Iq zlY@o(DO<$r@**JN5_ca1pWu;F=q*cJT0@FpjVw1VaafmP7=5Wrpuz{Aoy;wuf(N=g zo_wpg@V3%`pL<6zFz**mN}X-kNa;H+`ObUjcqLlUxpG#)!A{rrv6g|XtZn^mo0RZK zC2bGWorKp6V9^ z((@6F|NTNiFH{k+By9%tL&ad{(x)Lmc_2uZIT?AYok_mz54=y!F5Zd89* z*e%38Xyb)-##}se=-Ln7Rp~#OZb(r`_z>0qGlh5VB&QMBmQbi*b(9de3EKtXQVV^K zDt0Q->Ba4Zd-gL<$IL|^(tUC&zQ<-3I+h}UE5=ex&3N_KS#wqs0dFGZo5pNkNe*>4pZO{=2s>PxF52PTER z-DZ03Gl=tsjmO;DWpvgmaurBl?8Gg)C_W@PK(Wo4kKYcj$n)H;<7#{s&>INw)m#-h zVmVEJw6}9$Tl6+iIyERnFHF1*FVY$dJ-RIQTh!RB9deH7RKKxWaH{EsUGZD=Q3&A-1t9fc%1PR*F2R^yZR`A7Dho%b+vHv+DCX^vS=3iPz_Qb< zNU(PW;j>Stb8)*jLr-R>RJ`BvMBPR@6fKyE=W|EVUTp73I1dF}3(np{%Df(HNyab> zZFD`pU$B17m-J_)WeZ?)N{ukr3nhUQ7zyH#Xv0Cpg~y$K5DAGMqSUv`VqkJk-f<#? zQj7kT+^C?M#N$Mm{0sY+NzBaksSCP%Hx$P10W$krZvJNJ+dCi}%Ub7*GyXULmcYHb zwX15l&-!Qo?28ltXegw*_JEWbSC4b~f`%~d#~g#2*#&<8kcD!kK?B&@@6!H|7GTxx z+2QZ{i;>B(V#tBb3p;(wnD6aEYz(vBr6Sb@ko3Lh^q~cJD+!=$h6dLd=o~l6BWCsa zFIW8Hzu*N=+xVpz#6y7R))TZF(YoeAFC%F zYgqwiVb(0kTf~#ZcRc2`-AC{DN#Ggj86^QTe#dZsxr|=N(m*K*v&_LFScY^5=V&QB zvyzKVxPqI$-rgcchE_+*XelzYgNseXlv+p281)_oa$0BZBKWnyp1(}b_4s=Eo=*sU z3?FS-Z9)Ws#eWT`42b)+k8yd=BP3$9kM@4sA3D~(0-sT#>BJAbhST%&=;B*;VdnJM zB}Ljg2Zi56;MTkQ#J$5`Y@FqKwQStSpm8X^vxK*QZOY?<^JrB4AG6 zNqGCx!`f<3$HPeD%E(uMke$`d;Qh-s1_c%ZOTahe`#pBXmIRBOtnB1;^qlOpJdtl& zNMmqP@Oi^CgERHV7SZ+I4T45w|66#!83msM(Q)~+WP{NVnZkaLn>dowk&BPM3^wWC z{ep?2>cf3A(dLV!`suU5_m2*P$pnKj+Sn1LeF=zAk?X*ETVvxO)Xyh=c+)+Y6y9Yh z*g!n~h1PkpJ}5RpM%Nd15u~^Xji&dD;ayA*pS4G~{?#EuJ~KgcgMX-yg*3`psY#hx zS*hvS=lR)*yvGrh?mmXGY_%Mq;xnEH&umzhw}Pot``8)DO!((qjVaA_D#u_*Hz9UmZ;mk@B4L z>`ccdD2zH~zxp-4WB$6cvb3NYu~PfO+dfF$m3c;Mnn=dY_?=pSuH1KYnt;Z9_lr4F zsqh2J;mbso6vcnQ~|>hoKJ{MqY);pzI!tAf~NYgXj;)e3Xii#S>N zb;ug?y%~7hk7bndBM{FBeEQAl+l~7u)^PN(`D-64+fh~wJ2j%tj7tBZlzdDkFGMcC zV{u6-DY3J{2EZO|GSn}3d}%tcPT^#xx9 zTxn!+lBl@v@)^%r<(?7;&A7(KEW(LNO(_|EGqcoDMobmCVFlpP>e3%1_oNef}omN34}{q0QCW625%&rFr^bg6$Karbb2`6?wAcl3Kv4@g*}$ zL?37$|9pC~1QS7g`u^JS#;UCK()$$gd_VL$Ikb8D{n{~ka|qF&joWR44IQYKMC+({ z(oiMdZ>tfgBeid&x%8f=`?Xy2;Z5k_4fwVl4XCzzjy#<|IYH~(t&~7fJ;H^=8jqTi z{Ul`l7|ztoO2Ifx8C5Q|N#QJ0QkW{VX`YbJ;kIvYBPw6{uMJHP6%ms{F1ZAJfw7*Kh>f{~}!&P>McPN@5lhijP=| z?wqjc15-5UQblGfhki!=yEz_-~v?PQvzhDk@`KE1EziHnSOXx>P8`nTZ^TJHS5SrAsBkN$%7Bzqsf{^);synFuKeHeH0 zKtDe15&@$2v)AR4s+1rJ+`iY1CeR{p_r7({nn&CLp6k9-=FH%b|K(;un75X~arXxScCLwUZl z1HPk+3Pir`k>tL_kUo`qJ`fAqlipHzUoZBvzl%FwPzCje0?^D~kTdtaKX3&#DDP=` zKQ6$tNjF|WcM<=flXh}MB+a+4g5d)-=;8P&S^1Z~j9gBaGdw3z`^kBWj|b8X1gf#p rSh+~TjpJEG0`4;M`+)+yj3hxWPs1<}wOSZ}^06Q!rW-o65>0BSY0_>T*nlAu5<=*j z7;0fBm7O5?_nfP=6~(f3&+qO%pYQpUG1g(-;|}X|39>6X*X)R$_)8ndVHmA9oVTLh3x5R>h*;1*FxH7p~tCo$QbDKT~DAC9{e4sJ7C)a3+Ca+yn7eLFc Za81QLht{(@eA*{q3{F(005|;m<|8{ literal 330 zcmV-Q0k!@giwFP!0000017%P>Ps1<}wc0K~`B+dTrW-o65=|PrP0|iX3~aDuLP7{# z6GJWRq_PtP|DJP|wxU?F{QP|Hz4JYvGR8Wr+v~7Smms^Mea(98#9!KT8ipY}pw&Mh zU9Ak(5~_Ewa$~Ea$!k@WP&Q^IEsU};gkh^25QJQ#b4!q3p5200u0FD#__8PbR;DqN zNA)SST31!^Co6-|Fncq4ylh4Nm*^FuU@oG$;75kLv+s?E8qKP}o$d}`8jvXENkoYA zL@anrC`nT>6`VYy1wDHz70E1)W Date: Tue, 12 Mar 2019 09:42:58 +0000 Subject: [PATCH 3/8] add skip to test that needs to be rewritten --- tests/testthat/test_get_R.R | 49 +++++++++++++++++++------------------ 1 file changed, 25 insertions(+), 24 deletions(-) diff --git a/tests/testthat/test_get_R.R b/tests/testthat/test_get_R.R index 1c39d82..97e1344 100644 --- a/tests/testthat/test_get_R.R +++ b/tests/testthat/test_get_R.R @@ -43,34 +43,35 @@ test_that("Estimation robust to heading zeros", { test_that("Test against simple example", { skip_on_cran() + skip("REWRITE for changes from 3acfd712232bf36b91f3ec3906e545ced2c2cb6a") - ## ## simulate basic epicurve - ## dat <- c(1, 2, 4) - ## i <- incidence(dat) - - ## ## example with a function for SI - ## si <- distcrete("gamma", interval = 1L, - ## shape = 1.5, - ## scale = 2, w = 0) + ## simulate basic epicurve + dat <- c(1, 2, 4) + i <- incidence(dat) + + ## example with a function for SI + si <- distcrete("gamma", interval = 1L, + shape = 1.5, + scale = 2, w = 0) - ## w <- si$d(0:100) - ## w[1] <- 0 - ## w <- w / sum(w) + w <- si$d(0:100) + w[1] <- 0 + w <- w / sum(w) - ## ## manual computations of expected forces of infection - ## expected_lambdas <- c( - ## day_1 = NA, - ## day_2 = w[1 + 1], - ## day_3 = sum(w[1:2 + 1]), - ## day_4 = sum(w[2:3 + 1]), - ## day_5 = sum(w[c(4,3,1) + 1]) - ## ) - - ## R_est <- get_R(i, si = si) - ## expect_equal_to_reference(R_1, file = "rds/R_1.rds") - - ## expect_identical(i, R_1$incidence) + ## manual computations of expected forces of infection + expected_lambdas <- c( + day_1 = NA, + day_2 = w[1 + 1], + day_3 = sum(w[1:2 + 1]), + day_4 = sum(w[2:3 + 1]), + day_5 = sum(w[c(4,3,1) + 1]) + ) + + R_est <- get_R(i, si = si) + expect_equal_to_reference(R_1, file = "rds/R_1.rds") + + expect_identical(i, R_1$incidence) }) From 19c50d8fb60cbacecad336d0fd92542a1a2296b4 Mon Sep 17 00:00:00 2001 From: "Zhian N. Kamvar" Date: Tue, 28 May 2019 10:59:08 +0100 Subject: [PATCH 4/8] add correct version of EpiEstim to DESCRIPTION --- DESCRIPTION | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 994cdf4..44a275a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,16 +2,26 @@ Package: earlyR Title: Estimation of Transmissibility in the Early Stages of a Disease Outbreak Version: 0.0.1 Authors@R: c(person("Thibaut", "Jombart", email = "thibautjombart@gmail.com", role = c("aut", "cre")), - person("Anne", "Cori", email = "a.cori@imperial.ac.uk", role = c("aut")), - person("Pierre", "Nouvellet", email = "p.nouvellet@imperial.ac.uk", role = c("aut")), - person("Janetta", "Skarp", email = "janetta.skarp13@imperial.ac.uk", role = c("aut"))) + person("Anne", "Cori", email = "a.cori@imperial.ac.uk", role = c("aut")), + person("Pierre", "Nouvellet", email = "p.nouvellet@imperial.ac.uk", role = c("aut")), + person("Janetta", "Skarp", email = "janetta.skarp13@imperial.ac.uk", role = c("aut"))) Description: Implements a simple, likelihood-based estimation of the reproduction number (R0) using a branching process with a Poisson likelihood. This model requires knowledge of the serial interval distribution, and dates of symptom onsets. Infectiousness is determined by weighting R0 by the probability mass function of the serial interval on the corresponding day. It is a simplified version of the model introduced by Cori et al. (2013) . Depends: R (>= 3.3.0) License: MIT + file LICENSE Encoding: UTF-8 LazyData: true -Imports: stats, distcrete, EpiEstim, epitrix, ggplot2 -Suggests: testthat, vdiffr, roxygen2, incidence, knitr +Imports: stats, + distcrete, + EpiEstim, + epitrix, + ggplot2 +Remotes: + annecori/EpiEstim@new-version +Suggests: testthat, + vdiffr, + roxygen2, + incidence, + knitr RoxygenNote: 6.1.1 URL: http://www.repidemicsconsortium.org/earlyR BugReports: https://github.com/reconhub/earlyR/issues From eb5f151d42e1facc95497d1198fda90c788312cb Mon Sep 17 00:00:00 2001 From: "Zhian N. Kamvar" Date: Tue, 28 May 2019 11:06:18 +0100 Subject: [PATCH 5/8] uncomment test and add skip for us to work on --- tests/testthat/test_sample_R.R | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/tests/testthat/test_sample_R.R b/tests/testthat/test_sample_R.R index ffb27c2..9389581 100644 --- a/tests/testthat/test_sample_R.R +++ b/tests/testthat/test_sample_R.R @@ -3,22 +3,25 @@ context("Test sample_R") test_that("Test against reference results", { skip_on_cran() - ## ## simulate basic epicurve - ## dat <- c(0, 2, 2, 3, 3, 5, 5, 5, 6, 6, 6, 6) - ## i <- incidence(dat) + ## simulate basic epicurve + dat <- c(0, 2, 2, 3, 3, 5, 5, 5, 6, 6, 6, 6) + i <- incidence(dat) + + + ## example with a function for SI + si <- distcrete("gamma", interval = 1L, + shape = 1.5, + scale = 2, w = 0) + R_1 <- get_R(i, si = si) + samp_1 <- sample_R(R_1, n = 1e4) - ## ## example with a function for SI - ## si <- distcrete("gamma", interval = 1L, - ## shape = 1.5, - ## scale = 2, w = 0) - ## R_1 <- get_R(i, si = si) + expect_length(samp_1, 1e4) - ## samp_1 <- sample_R(R_1, n = 1e4) + skip("\n!!! TODO (2019-05-28): These need to be validated with expected values\n") - ## expect_length(samp_1, 1e4) - ## expect_true(abs(round(mean(samp_1), 1) - 3.3) < 0.1) - ## expect_true(abs(round(sd(samp_1), 1) - 2.2) < 0.1) + expect_true(abs(round(sd(samp_1), 1) - 2.2) < 0.1) + expect_true(abs(round(mean(samp_1), 1) - 3.3) < 0.1) }) From 6555fd9ae3b4792ad96da9d3c4f3a1b5a5cc4200 Mon Sep 17 00:00:00 2001 From: "Zhian N. Kamvar" Date: Tue, 28 May 2019 11:07:58 +0100 Subject: [PATCH 6/8] update skip warning for test --- tests/testthat/test_get_R.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test_get_R.R b/tests/testthat/test_get_R.R index 97e1344..0c98365 100644 --- a/tests/testthat/test_get_R.R +++ b/tests/testthat/test_get_R.R @@ -43,7 +43,7 @@ test_that("Estimation robust to heading zeros", { test_that("Test against simple example", { skip_on_cran() - skip("REWRITE for changes from 3acfd712232bf36b91f3ec3906e545ced2c2cb6a") + skip("\n!!! TODO (2019-05-28): REWRITE for changes from 3acfd712232bf36b91f3ec3906e545ced2c2cb6a\n") ## simulate basic epicurve From 6371f29b6a18a0aed64147ead70901dfc59116f6 Mon Sep 17 00:00:00 2001 From: "Zhian N. Kamvar" Date: Tue, 28 May 2019 13:20:40 +0100 Subject: [PATCH 7/8] update DESCRIPTION and travis yaml --- .travis.yml | 11 +++++++++-- DESCRIPTION | 1 + 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index ba12623..5064665 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,6 +1,14 @@ language: r cache: packages +addons: + apt: + sources: + - sourceline: 'ppa:chris-lea/libsodium' + packages: + - libsodium-dev + - libicu-dev + matrix: include: - os: linux @@ -13,8 +21,7 @@ matrix: r: oldrel - os: osx osx_image: xcode8.3 - - + r_github_packages: - jimhester/covr diff --git a/DESCRIPTION b/DESCRIPTION index 44a275a..0b7ef07 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -17,6 +17,7 @@ Imports: stats, ggplot2 Remotes: annecori/EpiEstim@new-version + reconhub/epitrix Suggests: testthat, vdiffr, roxygen2, From 467c65d3d38e8020b3fdf6e478d1fe20807e2f9e Mon Sep 17 00:00:00 2001 From: "Zhian N. Kamvar" Date: Tue, 28 May 2019 13:38:02 +0100 Subject: [PATCH 8/8] damn comma --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0b7ef07..1f5bd0e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -16,7 +16,7 @@ Imports: stats, epitrix, ggplot2 Remotes: - annecori/EpiEstim@new-version + annecori/EpiEstim@new-version, reconhub/epitrix Suggests: testthat, vdiffr,