From 80dde324b1119d546a131223d9b17f3d0553c164 Mon Sep 17 00:00:00 2001
From: Henrik Stovring <hersto@rm.dk>
Date: Mon, 24 Jun 2024 16:50:36 +0200
Subject: [PATCH 1/6] reverted try_sandwich.R (data read) added vectorized
 example for derivatives

---
 sandpit/sand_test_HS.R | 94 ++++++++++++++++++++++++++++++++++++++++++
 sandpit/try_sandwich.R |  1 +
 2 files changed, 95 insertions(+)
 create mode 100644 sandpit/sand_test_HS.R

diff --git a/sandpit/sand_test_HS.R b/sandpit/sand_test_HS.R
new file mode 100644
index 0000000..7d87217
--- /dev/null
+++ b/sandpit/sand_test_HS.R
@@ -0,0 +1,94 @@
+# How to compute derivatives vectorized using the numericDeriv() function
+
+
+library(bbmle)
+library(Rwtdttt)
+library(haven)
+library(data.table)
+library(readxl)
+
+# load data
+# test data for sandwich from Henrik
+
+tmp <- haven::read_dta(system.file("extdata", "score_ex.dta", package="Rwtdttt"))
+
+tmp$packcat <- factor(tmp$packsize)
+
+cl <- tmp ## [order(tmp$pid),]
+
+fit <- wtdttt(data = cl,
+               rxtime ~ dlnorm(logitp, mu, lnsigma),
+               parameters = list(logitp ~ packcat, mu ~ packcat, lnsigma ~ packcat),
+               #id = "pid", # do this to stop wtdttt() from filtering the data
+               start = 0, end = 1
+)
+parm_form <- unlist(strsplit(gsub(" ", "", unlist(strsplit(fit@formula, ":", fixed=T))[2]), ",", fixed=T))
+
+#####
+# based on code from bbmle::calc_mle2_function
+# FIXME: painfully slow
+# FIXME: redo without special case code for each density type
+# FIXME: redo to handle non-standard parameter naming
+
+parnames <- grep("delta", names(fit@call$start), value=T, fixed=T, invert=T)
+
+if (typeof(parm_form)=="character" && length(parm_form)>0) {
+  ## linear model specified for some parameters
+  vars <- sapply(strsplit(parm_form, "~", fixed=T),"[",1)
+  models <-  paste0("~", sapply(strsplit(parm_form, "~", fixed=T),"[",2))
+} else {parm_form <- c(); vars <- c(); models <- c()}
+
+# FIXME for any missing components, substitute "~1"
+
+for (i in seq(along=parnames)) {
+  if (!(parnames[i] %in% vars)) {
+    vars <- c(vars, parnames[i])
+    models <- c(models, "~1")
+    parm_form <- c(parm_form, paste0(parnames[i], "~1"))
+  }
+}
+
+vpos <- list()
+for (i in seq(along=parm_form)) {
+  vname <- vars[i]      ## name of variable
+  vpos[[vname]] <- which(parnames==vname)
+}
+
+
+
+mm1 <- model.matrix(formula(models[vpos[["logitp"]]]), data=as.data.frame(fit@data))
+mm2 <- model.matrix(formula(models[vpos[["mu"]]]), data=as.data.frame(fit@data), drop=F)
+mm3 <- model.matrix(formula(models[vpos[["lnsigma"]]]), data=as.data.frame(fit@data), drop=F)
+
+mm <- cbind(mm1, mm2, mm3)
+
+logitp_est <- mm1 %*% fit@coef[1:2]
+mu_est <- mm2 %*% fit@coef[3:4]
+lnsigma_est <- mm3%*% fit@coef[5:6]
+
+
+
+myenv <- new.env()
+myenv$logitp <- logitp_est
+myenv$mu <- mu_est
+myenv$lnsigma <- lnsigma_est
+myenv$x <- as.numeric(getElement(fit@data, fit@depvar))
+
+nD <- numericDeriv(quote(dlnorm(x, logitp, mu, lnsigma, log=T)), c("logitp", "mu", "lnsigma"), myenv)
+grad_mat <- attr(nD, "gradient")
+
+# also possible to just compute one derivative at a time
+nD2 <- numericDeriv(quote(dlnorm(x, logitp, mu, lnsigma, log=T)), c("logitp"), myenv)
+grad_mat2 <- diag(attr(nD2, "gradient"))
+
+# The u-vectors in the Stata output (ml score)
+grad_vecs <- cbind(as.numeric(getElement(fit@data, fit@depvar)),
+                   diag(grad_mat[1:nrow(grad_mat), 1:nrow(grad_mat)]),
+                   diag(grad_mat[1:nrow(grad_mat),
+                                 (nrow(grad_mat) + 1):(2*nrow(grad_mat))]),
+                   diag(grad_mat[1:nrow(grad_mat),
+                                 (2*nrow(grad_mat) + 1):(3*nrow(grad_mat))]))
+grad_vecs[1:10,] # This is the same as Stata's computation!! HOORRAY!!
+
+# Each column should now be replicated and multiplied with corresponding column
+# of model matrix
diff --git a/sandpit/try_sandwich.R b/sandpit/try_sandwich.R
index 19f9767..6cc4d74 100644
--- a/sandpit/try_sandwich.R
+++ b/sandpit/try_sandwich.R
@@ -8,6 +8,7 @@ library(readxl)
 
 # load data
 # test data for sandwich from Henrik
+
 tmp <- haven::read_dta(system.file("extdata", "score_ex.dta", package="Rwtdttt"))
 
 tmp$packcat <- factor(tmp$packsize)

From a693f52d3d49da7b6bd4597db9cc687a415745a7 Mon Sep 17 00:00:00 2001
From: Malcolm Gillies <74748374+mbg-unsw@users.noreply.github.com>
Date: Tue, 25 Jun 2024 08:48:12 +1000
Subject: [PATCH 2/6] Update README.md

Added brief example code
---
 README.md | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/README.md b/README.md
index d025daa..4790ce9 100644
--- a/README.md
+++ b/README.md
@@ -19,6 +19,23 @@ The parametric Waiting Time Distribution has been developed to provide a data-dr
 
 This implementation of estimation procedures mimics that found in the corresponding Stata package (wtdttt). The model is fit using maximum likelihood estimation.
 
+## Example
+
+```R
+library(haven)
+df <- read_dta(system.file("extdata", "ranwtddat_discdates.dta", package="Rwtdttt"))
+
+fit_r <- ranwtdttt(data = df,
+                   rxdate ~ dlnorm(logitp, mu, lnsigma),
+                   id = "pid",
+                   start = as.Date('2014-01-01'),
+                   end = as.Date('2014-12-31'),
+                   reverse = T, robust = F
+)
+
+summary(fit_r)
+```
+
 ## More
 
 ### Known bugs

From b8d8fba522d4a99415200613693b3b2f06baa999 Mon Sep 17 00:00:00 2001
From: Malcolm Gillies <malcolm.gillies@unsw.edu.au>
Date: Tue, 25 Jun 2024 09:16:12 +1000
Subject: [PATCH 3/6] Removed outdated comment

---
 sandpit/try2.R | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/sandpit/try2.R b/sandpit/try2.R
index 758d822..b78aeab 100644
--- a/sandpit/try2.R
+++ b/sandpit/try2.R
@@ -1,7 +1,7 @@
 # Quick test
 
 library(bbmle)
-library(Rwtdttt)
+#library(Rwtdttt)
 library(haven)
 library(data.table)
 library(tidyverse)
@@ -71,7 +71,6 @@ summary(fit1)
 df <- haven::read_dta(system.file("extdata", "wtddat_covar.dta", package="Rwtdttt"))
 
 ### fit waiting time distribution (exp)
-# XXXX need to change wtdttt(... subset= ...) to use non-standard evaluation
 fit3 <- wtdttt(data = df,
                last_rxtime ~ dexp(logitp, lnbeta),
                start = 0, end = 1, reverse = T, subset = packsize==200

From 19617ffefd96545f33e2bdf8758ea37535c79b34 Mon Sep 17 00:00:00 2001
From: Malcolm Gillies <Malcolm.Gillies@unsw.edu.au>
Date: Tue, 25 Jun 2024 17:26:13 +1000
Subject: [PATCH 4/6] Removed old sandpit files

---
 inst/extdata/score_ex.xlsx                 | Bin 0 -> 35456 bytes
 sandpit/d1functions.R                      |  60 -------------
 sandpit/{try2.R => examples.R}             |   0
 sandpit/good_practises.R                   |  31 -------
 sandpit/mlwtdttt_lnorm.R                   |   6 --
 sandpit/sand_test_HS.R                     |  94 --------------------
 sandpit/sandwich_estimation_code_malcolm.R |  80 -----------------
 sandpit/try.R                              |  97 ---------------------
 sandpit/try_sandwich.R                     |  44 ----------
 sandpit/wtd_bbmle_try.R                    |  62 -------------
 sandpit/wtd_class_test.R                   |  45 ----------
 sandpit/wtdttt.R                           |  50 -----------
 sandpit/wtdtttpreddur.R                    |  11 ---
 sandpit/wtdtttpredprob.R                   |  17 ----
 14 files changed, 597 deletions(-)
 create mode 100644 inst/extdata/score_ex.xlsx
 delete mode 100644 sandpit/d1functions.R
 rename sandpit/{try2.R => examples.R} (100%)
 delete mode 100644 sandpit/good_practises.R
 delete mode 100644 sandpit/mlwtdttt_lnorm.R
 delete mode 100644 sandpit/sand_test_HS.R
 delete mode 100644 sandpit/sandwich_estimation_code_malcolm.R
 delete mode 100644 sandpit/try.R
 delete mode 100644 sandpit/try_sandwich.R
 delete mode 100644 sandpit/wtd_bbmle_try.R
 delete mode 100644 sandpit/wtd_class_test.R
 delete mode 100644 sandpit/wtdttt.R
 delete mode 100644 sandpit/wtdtttpreddur.R
 delete mode 100644 sandpit/wtdtttpredprob.R

diff --git a/inst/extdata/score_ex.xlsx b/inst/extdata/score_ex.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..9190808f4042f05069df8203f5eee8f2a79eefc7
GIT binary patch
literal 35456
zcmY(qWmp}-5-o~Lkl^m_?iyTzyF+ky*8~W*ao6B(!Civ8yM_%rcwplec$;(Hz3-j+
z!F(`1-%M9muUf0BQTq%FhXVx#g#?u%3(!^RA4^AvhJqr4g@VF{!htfBa&h*wa`rUU
z@^iKFFlO_0a)gG$_|I*~Gf?kC->oHiQe}_>LndVS3U63NpCUZjQ{4<%KJ_EV3{3+t
z<IGOsjck+=OY%5WMmKm%*uz>@(}})#98*$jX2o+*QL$ffm1%SJ#^TLrYtc#bT$PMU
zL5J8uh~f!xrkJEZYAhEqCv^s%nnP4-wkNUp=-X~Oh%{q{N{_{LMk40~Q^;*j{zutO
zj~VQd{J1mcOilKOi<V@Vl$nIqe^QJyaOqS{^ym|#IhF;1X&lTM(1-1-Jdp(k54yB<
zr>Ws9<6#PvEKg8Haj;NC(^uoMYJ97*JJS95Uo>>Y>|Q|w564?$b)%yx5$z);si#SY
zJ6HcS_Vr^$$a9!g%rKxgd{FHQR}sORmzs#PD7%5MYX9TX`QskA5{QrRAihF@xYX3$
z%F%<J?fvJwE8qY8Ztp|hEiOR<J(L4Sx_;rM2gjJeg26Es*)R;T0;EY~8?vY~yt!FQ
zfX_zHeNCsy3&>lCmEt*oa|B=!2+8!g8T7XZyP(f4EyLSnxl_fD$fYy0eIU@hlA9?e
znxD{1l*`Eb7MjL5S<sPTpG@>v$8s+c$YnNdnTJZ0^402(>W9TRTP7}4z5A<&Ssb)5
zbzds4wxflF;d;2At@s!l!6@Bz4<}~-5dL?9a$?PZDae;o!$3h1K_=+y$o}6s`nb3|
zn7g<*yhqLdo$vdocYA;0-D;D5#SCy@M3~Y}1?A{S;Boq6{;J*#`kJVcW1bV9HNV`p
zJ0_{FwiZ@Ssr2Vr+*&7RQlNqx>z<RKA_fZw0Y68~6QTEXdz9XBr5YAp=j!8Kka?OM
z4DltET_vw@LNSPn4z*)ZmZ<iC-b+Z|r<wFbU!RR<oBMZCrTu5p80hrVuSS;W6;Iih
z;*1SW0g;<WD%1#mSY>Y&1$S}K8~nIO_+as*_%lq}*)7Fe=M+Y><gA|5FD}SA?1h1V
z`}zYivGFMvuBW|?Q`7(T%aP*&=`Cc2CXgB9{MRoYwqM+>EHypd?VN2q{?BH5zyEIU
zW8W<<c^0ynP$NxCC-O=+F(frNNi*-@T<eOKjM1ja&GDIj{E)V?wGUY1%nx}y(_@S;
zK5k-13ZF?g?;4ju(CI=rFU4`GKa8|-tNB1RV?L+BC}bUynCyutJ3r$0vtgmr3}@7$
zUH%GMONtaK$K>zTivrPi9i<4V_R+l!Jl@<~GE}`EG|>LTZb~O}i_3g`K9AW2z56M3
z*ITi^*-TW6&Ob?XFBuwt0S}?r-zeCiszaIUZxHuB;v`o~((!>gxU%7I*QV6}iP2}~
zdI&6S(r&QM2gl2=3bEO`312PQ#Ff8D{>w1kXY5<y2f4JznwCAm|F7#tz6Iugg=~x>
z2n?b9zpnH2bF})uK6}6MZtsKMEiXk`72>mqH^o$8%RD&AigWY=1*3K-tkoK?P5DO+
z;J`@zs9kGaO5&&77tB6Yl2umzfZ}z;@QU18U$-70s^h+a)pJZ=k8}OshcY&5GM;W!
zsIiKjk(iHBuO?AzsAiVk-efF!o@@lsI1wh&<eyqOqhs+-mrryEDp_X3pD*zvPKB2H
ziydTP@wjQerps@SNId!$oo|YRHz>r9t;&P(SjGNen`GoQ5$NyB0(cTFLS-+9#Y)PQ
ze#sac*inpr)TlOAV`Ov(SyCORu5jB<dla$tNQaxX;ruLq^S@McJ7uw)(z^fbMSf7I
zx61;TQW)UF)_fz8!GEY_Y=4vpCvhv)c+@Q)88wr*$h0oZ>?QdFXZ;8FkjtmjjAo3B
z?cv-z(x&2dqCk6Qm4MU{SNF=M;v7yv`88VBCj7PAzMqb}(=RAq;I2-A>Z}_R@*i?d
zVL$WID}V&DnVUa^e?0YN4(D3gX)?X$-uVC8iBl|Y6rcK*34G~W-u>U@bIY%Mc!aFV
z1Y~tS{&)HQ|5)%2d+)d2?S0t0#U-jbKo)SgFJ?Z-i|B|?GMP*V;rp)-r`|qJrJT3<
z6vE2luesrKWC6)9H#hs^pgU`{E-LtQWB5~T1k}76#V!x{zR}qXpS2ax?_>9?*Ma2C
z?;kb@x<aoh44bKH0wdX!iBxrlZc4a6k`}5J%EYc*icFex;BNd<VmKm<a3NxKJG63{
zA10Kvl~&T1s_yzZ=X}EWgNn1HD_!Ixrq4Cp-wd*-RnXdt9m#+(dA`zEMZGUh`3Ej}
zu~)NP-_KSmdNsdo&Q|;Boaz_EVT&cpBkq{;ktVZoWXYo9thoN<v<1S}9LGJ%Z?rvM
zQ5`4w=ia3o-HQLO=fHG($U9#lI9dk-MezT5&coKq%F~1WKNrsbgB$N--tB$-yQTIx
zrErhBu|DYHyws`p^a(ibwxeWbX30|X@p<11v(C2;q#qR=oOZnhU*Ax>cdB^P&ZUA1
zwgtZ{{r1%UQyKX9{08oOyCLs;G!=h+$q;|NIShH2>Fc-%zTVFdd0QlZeK-vH7kTlT
z|MvX6+j{XjD*pN}<n89V?+vsY`~t4;!vu@JfyjHG+xlK_#ozw=^gT~ryq%AV|GUg+
z%YVHd6`xEAc|1OZ{I2^b9`d@n`ws-UfFs2nLvG679*kY<UjC|zzkp|VgADgBf`hH%
z@?SxwuMf}U_2U0t`(AJM^WUDFlJY}d63AbfCSF0iuMZFqv}$?#7n!m0JUA=<(nfFx
z5EI6&-~I8@xRW3JHsK!fw09WtR{!>Lei8Dt7!srS_A>hRuy6X-=L9lczB)$iS;&XH
zB0uD5QuXa>mH%<FZ@r?g?zQW+i#+&k+vWA<x$pU5RU9(Wk8dyJ;;(}~n=!lp_IKa*
z60YKY^FL3|ioZ^Z2i}YS+Z8Rk-u?HN{4LlJ`1Wr<;q70>cHbL<_;1L=i+%fV4-bRl
zfrk`N!nDu-the7DuU%d*_xoNhP2a$)Z>QyN!EfM%w<nhnT-yUv-M-h$Rq|KxMab)M
zUkLSXC%#i~LdawOYu9enmg&LUtLfX(>#NV(XR+V+Apv(GXZerdb*nd)*x9;g8@Iqa
z)n^;Kz`K91`egarce4`@XF+c@@v~xgvl|a*9&a@Xv+2LQpUnRbN<CSx>`?~sChYC4
z`dkY?8VCP#8>kmodHBm8@_1)4o^keYxGwJHvdy2bp}^UxVNvk0_w2jE{^={>o5{#-
zP+!QM;r)u@ePIh$(JUp+Q=0RjqDx`A)<l_OE#=#Tp<Lg?GsCX?%*^+dFqc9t);G5r
z<&?KF6X788{PL0vscnxtQ99jbw)>Q;)1+D|gQgd|)vJV65TfO<X%g-|3WMb^w0<K$
z-EQN{PsW?n16>No^I<>sD1_FC1WLy^riYX2z{bEZ+V2ZTGQGp#(KA7ujB1}N73p0<
zJsQ$(qzD#-GCT__JyND*gTI2`We4&jWF{MuB}~?|f+$n|jwqfsKts9Ku`LGN3z*pT
z%RJ%Yj}OTZxvBfRR+bsqC$-TntOO`r=fB#EJeOPcn=kce>;z4VJyVGDhD5gU1go};
zkQh6Owd*MslfLDZSk!4>{tPOu6Beo#L&tTiJv3qNqj-O%*7!o+i}0@jJv&jpf)D2K
zINyTvM!|}(ELZhE$LLtDE{#18jX$>_ixj0~gsU~+m$BbwqiI>f;%4;?($(uYwA`<l
z7g?P0Rh<SVu~!KtDXYIGevw#kZ=9F~cIWE54}mA|+-jm<&;|Wx6Q6k;4d@IeB&g^q
zF6?|PN^jRf)im!rE)^5@(-4eOP2JR$5KCh(<Wp0qmo|E^g{gXKGI}_UJy{Xh?q=5{
z2ytC|9ukwg(M50ARN(fCV*Bza=*PZoRJ~}cyr}f0G$#zhLnmdViRM#|`t#@pq-5lA
z-1sT_3oaR$bgT8K#yNK9Cz~WV;jk;H<czsY>}wxGf4<da(e_BFyb+StZyxQ<6aMb=
zI%$&K`_qB9Un@@8dg(WQ0?h-1vsCGKXRjBw)<H2Zuvf8f=8t=V<Q?>zRLLl1KB7^|
z=|U>*grT3df;t+nn*8PA1r{~UDyZ_V3g!h)uWsE7P=GsTmh6t<=DjYA5*cLCz(l9E
zCF$K^D|V$7<D73X){dOTFj&)LZ)1D9#>qhjel(D9_2rE@<A8}>T0&SHxzO4vS(6I=
z%sfy694#c8p3K4whz|Uln`{!F8ZTjDx1KG3{EWQrWnSRea+QV?Z${1z`1t9h52iBb
z`7^mG(>j{f>bc;GmVG(vK_h-y{=)_ug_lJDCH#;b(y$FO+R(U;oHltMifNI6FXaUC
z{Q90h{d)MNAphAo6gj!}sQp29^+P}0tLoqntFj{O+-Uqy?VLIW^axR=%XwG*WTVVz
znU{Zc)yDb^DhS8`pIPHO&g`+oEb+75zaZH81L_8)hGPU_2?qzytEU*Pec`XJ1rQFM
zc#doR!tHv;k@LIFXJprDdcTfkATcpEx(F{zrOW*sILrRB=w3tPh&7tso7dH{;O`hp
z0)6#y@0L^F7H%}2_?A!L1qvM>EStu*mP%U85$RdDOhC318rk-nvx@Td)3?v`fB8h@
z(1Q`|NS+j62-gK@Y~$uRo33Vn?gP&BhH0k8x#;jo)W)H3tiY0J_CX}hDFeC~9l;6>
z&Aq{G1BW~V%3!Y<#~X{r6{HmQ^p#DXT5`V)iapgc;A|Te9wVDmqP6qm32W7gp%JUy
z9+PbxgPe!E%$=sU)=0}vQGMgLuc@{@=|dEhT)4S9Q`!Sc372i$^6}cx^H`}i*YDA}
zV#t#u-(4L93cnQM=2HC15FEN;Z-Ps%emgwppHp-@^!bVHhos_zLGmj7;?K;^`UZIe
zmxksV8N3jdL9lWuhp)>I-_P{y1??1A9b<MWIMvTHSCvu4XWc_MqSigPB-I+(fuGn(
z@bq<9HDv%p28|bfBS|tVWQLSG)w#IlXkF|nEwZVB;1;p<fRFR3FbOadI<GRnGYM#6
zpC?mm+XObjs)V#d#FIa$gHum|a?n3N&`Vp=u}g@trP*iL3RiDO4q5SIx6zEQHPxqu
z%_c;1xtmR@7tBn~dH^SteB5+IZTV$q7{OlV_K6-{iX{W#Q<FgI&R$znm+RhOOFYBd
zhB40bVJV-u!Jf%>0?oc1U~rGjtN(C^j-Vkm3H&d-d@+6vgPn64w~`^u7OGq`*QDkG
zz6Tb$b^#3SI{}2%2DZkZB;lObJzuR>xioq9Vknye{MuFw5~f}C>b%@^2nE;^`(FX@
zcf%3m^0pHIMqv+>!bt^MQ=vgiz1S=>7xiEJh@46v`cNCO6)UJs!*lOqM>i`*as(Q8
zqFvL2gfGYMR0cXg=S*Xtm$ov0{>x-YHXlQpUqY}XXd$s@%#lB@<o@Q~relzxEN=$4
zYzWS*;Q{Ug)^GV*PNU~=)Z|}~N4IqGm-)C@Mh@!Q){I@|yJ}=)3FyiE+4M~k`!6X*
zG`s6|j6Yz>-%{mNe`r<^)6O{P3l2cxMd4Pb>%<7{HEBLGu+F~20Y$l;N$l1RuS65F
z@j>V9r@V0LI}DBEFQ4YbVN(Two-u4NXbN5wPwxDTVvkK+z&((usYHpSwNaO6nh(lm
zF9mHlE$!Ean-oO>D^LSJescI+z)T$Va`n4uIsT1v5$=+?mvLp4%Rx}1NqD{3Sk0%%
z$D!2AJiJC*!zqc<y?Y5Mi5+WolckRSsgbHRjnsrAS=gmE!;O>9y=VBqA24CX+MOW5
z=i@iu1DHzLmM2<EoYw8fF<(Za<XZZoj=T#Fh`s>v+WFJ~x?gy$tT)71o6@6;aQC#W
z#7~`BTV|#bE$HQc`?><VKb=Of8}P3+I}_L$uXhD*|73exk~E#LGK%%vHUoG50(!xM
zE3r%w1kvol9m6HQ+v;dr4bd}LqsvA@raV0<97tk*D6|ymANK+p5~8VOf-#V>AJhId
z#~Q@9q~`Osj|JRp8$fbjp&vVfrI;NZbKJwF>BAa#AYxXUv;~!estFrL7L!rWoB`@b
zz$z(O)tO>Z#_)Z(gRWw{p!ch`^EUw%F#ab2(AvR)+LudVz%xnEJ=@p@aIf+0ic?*-
zpT@-Px9*?o+y}1tp6@K^+P_85gH;0aGGjbbSqhnUxqDf}LLD+%h|i`Yk)~;S_3Y5^
z>ka?l4?37-XrI(&brIKW+h`5yBBzB3l(?ez>{ke4<gJ5*lP;(|Ke|L*(C_M#$myd9
zEM8)9e|zCvlreXt!RKm6ZA#B7!irKz9OXHGM=^cMI_HFV(uh&ZT1+}qIG3n9GHHeO
z<!Ex+LU^tXbIsnVb|UWXU*^7%k@omhyazSxzxx^muwZa{`H^x#Oh9;kGEVAiMN-YF
z=5&k$G-u9<XS!Waivkjd6u=7`n_KI_p&;0*D9|3b`8@Dt^<VCTP!q2_#*u+~DNOxj
z0@5>8WSBB~WYR|L*UZ5Lt)QwYZR3&QOcHcFK@t{s0Stv>e|z9}uS81VY*m~s`;Q$O
zfuvndSwS2N#gjy0u?JR4X5KM`=5K&P-gzmN>eK^f^)(^-J`$$CAWJ-(grJyb!oAzI
z308@ss^{<c-{s_kI0q;<6N&srFtSN<DZaZcR<+iAVWkTXlR^6U5w|GnOxj)BjZj6*
z$gIpA0@VHv6tbWpBxBCY3(D)y(K$wz;S@fF9_hEEN*TGoLK-(;{Zs-4b$UK4{nly@
z5U1{*IfGL)KkO<Tg#cMyFT;VE4Hlvzt#_qv^wNPv(cN#heRCyBM1Ul^j*y#LR)^pO
zVBf0g1V~%;H7?2|yxpDaH=J(Pw+j}Qxog)^(9-Ce7Kn@VG+A=PePqemOEEC3f($GE
z3MlXhgt)SWSKgPFIN80%B%OTkcKB1I{UAWo-eIc!$MCkBbD**i+X$-XPlBVp!ccsW
zfDiugeypL6>EUD#@hM4Ig)ujmPmd7w4lHIH*2?~B#sz2b*|%R_HPJQ15_;9P7paiY
z6#0ChBdJ63i}2XX+NR#D=dnDnn?(%f6wrw7?e9*|p_UL+zf&bJ+dloLy55bmtR_rG
zACXxut2yyE?c|pQVS&>UMH!bua0XVfz!WWf$<4LZo-lx#Q{(RCq&j>m@39e*6qJo-
z@U$m({u3zytYa2&FId~N(e>DxptFpy3VzV9J|?lEQ@}Yo-zI6)W*QRoxHyn=Y;O|%
zuT>)Yls6*e4&LUp0gmqGAEU1eYTq%xW$u2pWP@9B&7($?xL}Friqx<=RgL;X2p}Kr
z^epkDogsK?t}(Foe~M^!UAcCyZBdC$Bv<$}g-e_;u3<rWfxk@ZR3#j*)#l3#^Gjy{
zA*mvn2f;03j@Q7*r@7r(L&`fc)AgRs-x7Y=kVet~JvAbb7_2L$oZEIY11z-kLp&hf
zsY?=+vnYxQr2QRbd#*ZJWmQyuMsm5P-I~y{BE(rV_8Qx4uD$#Wvq$<5U-3-qk9@e#
zXMw6AewG6BHE5y|YVMtR^X}rb@hO9L3NRbn66q1ej)H8;PVPgtSRqb_->_-l;)1;D
z6~k2`+WPV8?|nUN=?MjNzwemL36FZyEbSRl0-TaHEPtZW#v})Cg-0zR**rY3VP)RB
zmbzgD&F8`|fWEEMBVw%oGq(dN>xz&;4R7UMn<9p1ZX$lWYfiW!4GexxDE7q{=12y8
zUIYp-4c}ejWV&j;-{zyh$}oQFw<=(w>ST>ok*$Qge~l4;ASAQW>--)bYs~nGBBJjQ
zK3?Y*K5g(Z;n{HawX!k>lNWJ&9LV(oP5*3EJZvZjXp4_rVh>Gbv8_G39)V@k^(I4x
zbBBW?U_{(<Jkc`}Axv|6;Z&!?S=#m1`7c(A((+7J<QDsi)A!CODLYT@YaT9C1F8Fa
z<xg!>R+6sVl1z`2<nkKo-QlmTFcPK=<Nz)4cn;I{j^arvQJA{5>a@_qCB4^L{h&z~
zrg%mhzn|@a(O3xx>)Q61YVLM_M;rLu7YdU3D>6?Y*sN*7XC!E^O*<n9_Of<hM@|yB
z2!g0Rx$q#jwfPzQcxp(n=MK?fHRuXDz}k42AlU5~)z4Ek6sIju8Ax9w_g7Q8)Lf5y
zQ6I-rr98yh;+Tew7sxw0A=Uo=(-t+zLO$3sjIBN7bmyMvzMgNj>|R5!u>-OY4DL1j
zT+2-zxXtFf#Cb2N=<xz9E}oW?E7N3d{A-y3xFS!R@gp31+yE1?N(vc-H)}|dMC4-J
zry{j|S0CVUax0`QAquj9A*5643dO*ctDvG6s37^45Y1;|^k<F;5b^`tly1Uiu3zvv
zZs;AA8(NUoh;yOJyiGF5G4Es)<b)Z8u+EPU934VABDvVhyb-?iMz+fWlHzq#c$n2)
z582I*35{8By4z90D-6RP41h)zBKtn+JDJ|QzfM-twc|G$m@w1)wGYW5^b%Y2^G=M}
zvT+S0^c!FlsDQjxFbCrX)Uv&6AX|EAImIa1Mlzw={yoJ&k6^-*tgAj*fN{Kzq^IKx
zla{6PL%bSNP#OW9MVZ^bkxbv44Y%Cg8<hcQ{tpZ%)nhtph|3Lh5<8@i&beY$HQ`6m
z*)9J1QQQkPX)991t~$TbKRwQ}A4rX|M+l&o)X5@&3k|DLjdI2)O9xUXDW>B$cnwg&
zmujS%#dwX2qls;Q1ZxYg$za2e7U%+$;m~Ny9VgGe^)$7&bM$s~xBPMxN^vVv>Tmr(
zpT_B~7$~?R+RT(w_IKI!KQWnR@r4=l3ri1^7$58lFg=BPJ}pWiEq{#m(b;_f7oVOB
zeW6>}I7No_<3~c{68KKCC>LO^FLn<YJ@B)16W==iTspRhnYNwR1<T(gXEe4_AlBE#
zMY<|^CA~S$aKAwabC$kcpc3+^LQg&%TFol$?MNloWfd{OBX0op{j@4KD^EI8TI*F8
z-#z&0Lh@#Wa6A<2*+_6~2{NOf;TgfCWv>+Ab^pnk?@GMtKPjIQTTo)vq^8^_zZL{2
z^FJ%))7}<SFUimu+B1{<4dEy575t4*U2qyx>iL)^ALDSujV`yfV9o&p)J44nXO`S+
z(MKC+iBG?cKgDp6V`oTYGr!i4Q~gYlF>Asb10?|cV&v4cpVNpHY$B1ckJm9taURK7
z+<M(>L&d<dhRKW4yh>s^GB3xO##&w#>DrTpKaTR4Oj-_EcDh=Y{2iq~Wp_WA<cXgA
zCtQwBEf-g#vsPE;hv23?U9+a!O2Y*BcRJ|+-70u~wMz{0P&h6Pt>-Y=>O}@9`>7ea
zy)M&ig3GbP^Pye%M#G&^!p^;@!Yi^syT{i{Yc4gFXm)H6Z7f@*|DTMDcYoFv?m5dC
zql>@qgVxgm$Lm$Q1Ct=P;k}q^m)l*`pOhc<c-8s37z}T`%n|fjDVmq2uu^cm7PF(=
zaX+?DASE-1o|rgau#ojf4Qb7H^+?MQD|;rfehw^{@QKVn9Lk9P+d<DW07AR}58d0i
z@}{ozKYeF=rvONi^T~_oSX<zUbsT{OJwLqmgOL=oby;1Cj4?g)n-C`VX+Q>p>fHC(
zPow0X#LK<_!aQW`C(W_=`uKLv=#KH6thtGi@zpWY`%!7kGhB2DJnGX_$O`)0USSnR
zqU@i9U-~;*N>w_25dwjWE~6^gi-gyEgt3yEm3w1j4tRAnz-LgV2yTBa>&B@%M^(+`
z)n9!Wk@)8_HRwm^imbnh;8^=6%(NPnzsvxUdlEtSbwF#<#4T)P#Yy0y2hmWu)slOS
zEjD|nP&C}~xJ}DNrXF+~Jg?6bm+}gc#MQ7ieNMGXuGl%sYfh57wiX$XV*1+*ufC;M
z<Cy&l%vo_bnnE)Y(^gMXGu){Bj%1Ga*6%8($Q)ao$0VxGwVSV4GS*b|oj5%9I@5RW
z;=n>CX(K)U53HW$?Nm&6se=E^Vu%tRvHI1p4t-9oN-9cTE3@2B9^jx`xEEoGpji(c
zpjHLXu(ml$MMR$c8>cyk7hl{$BlB1Wu6i%!Zst{hKo-HhB^ka0tv@f)o<&h~cRLK0
zGvLOKu!W?$<2ExjLk<2CPhXT1CJ!M-Q+GT85|+}vM<<qx%h6*isqWvzPiXn>s@Ml9
z@-pjG@Qm|VI<s3I6iKnjmAxjEq#Iw81;fSEmHBxs!KPOhm%KhDV;qX{X<4P{dGClk
zs}VU9UC(X!>qGb$J$T;=>?S{sZ=83ob)@9!Ac2;Z+cq{dVo4dgN-<+H!fP1l)v8gb
zfM&|E0xe3YS(7K+xWxso(GFHn(j5oVSCpZ%bf-BW;aPmw9h=K$SeV<Vs-V+MSDtK0
zcbBE)Z>nb2cBu6-hij}#784<d(1X|q_Sm&%4>G$$3Sm>(u8Y*FD%n2Nfh_$)$A+18
z^E?G)dUR`<&Gg^H^l`ae?*Jm@_Nqe`X7*!$?_Oz&1EHl*rO7HbZ2sCYUOVYtj8=^e
zHCYol+3=?Uz+bz8use%XkIfzMX$_RQ^-d%QP5=QjeD}GG0YbVsmy?`%p$`+?$B+ON
z^Qd2^%5RbpOCbW5eD1~!lkz|a(oT<##>H0<Nok)ROfC-JCdd0|I7-5k&B{Qn)-_wm
z=5Cp;rWHL}ZnXlD`eCygN_X!PM{j<Mk3*ux??(;qzyTn5hVNgl+{uy)x>22c3}Y3l
z7semBQis;DL&v09eQ?|O7nazU2nv?9w{fL0p+s1qYXNtD6)4EYEosD3I+E+=>s-Va
z_nJ7DD+kYdw%(a{_-OLo%iJh!Ezks>POX@9UjtQ|U~f8KqUzbN1UIjrxsqlX+;4f8
zPv-FSl<hY^H(9t(<2{N5W5C8XN3GewzCtI>UJ`v@T$w@)jF)=-z7fXri>QZWi;A_r
zzEeBTR%4?rOD~qyA$GU$$JI(fLersE2Lu|0I3-CibaVge@<NE#m~-lGihi(&k4mDk
zzwLIROMIs<`A`S%J+C2&d+S2JLY!%tlZ+>_*7mM*R_3Nl6R(sD>-cyX-ILoNjz+ek
z?h$z?&#BNrl4`D;HY@!Wl2bxsnQvJW`n`q@wb;x<Sr~l`93hvPHM~U!0K6IR*XzrN
zG_L%;j|4KSj+le%6JKY9qSu(N$C;SdSgyZ->8-@Dt*jfK-o@WDi_U*CV_xdxyLLDO
zXw2}H=E|*LW20Egv65}JD?)dpF4Uaz8m79pGJurA8POG+PX^&y>--OrW>kWeN|r4!
z^050Wnqmaa8(0J`yBI5{8l&m3{<Z%h$!E`p<P7QAaES($X^n2k*N)1LspGbOL5t1C
zM*_XbFxYf;$dc*LhSI?tT^rJ<pbM34e7=^M<~0nzx>P&=B%R5p<$1foQAQLuPC&CX
zJ0N^x?*asZkdnhB27tW-<)27#uXeQ6yegDDG6%zXkR!$!!B%C$BRSQ`<M)~ir&DSP
z@#UQg0}#R2+M9Xg$UUEhRZhw#YL!!URYl{KA|~$<?m3NzoySTi<iNM;lOXG18$6RU
z#L6wo92rewpSkw8$2v6AIZX=sS}t*dJ_&J^V8&Bov=FGhbMx*$|DjCLXyw$mK*hkk
ztPOt@(l)c$x(v{0$C#;>X9IURih8|*iB&4T`5XI-jwym&S!dm_c@ZE@&{*J20HqPj
zntGHZZ8wlvdJ1m|Dl}uuc~ZOD9q*pN?5A;hIa8g;zC)4Ha8iD?y{)%OVAjzO?kY+x
z@~meF<2AQ;qb%xAVm!?lmh*gM_gvo=`Oxt5iYqe0Nvaj~;3HS}9q_x|nfd(MK0Rl5
zMRY_rf8lcNxUJ5`RxJiRF*6$n>C-c<XO)F#lJ7vp0Z`^=l+DcIrdYfeBN{2oQefD-
zcqp^IVIbSjxQD1PGQ)4dvBn}MJHOw-C^Il=VVJHIC24_-dyVG9$e2lV!($@xkYI2-
z9w3zskI@(hvwL{;N=)u1!#dM7#3p<$)11sLuDraNI}gC1zom)XDU+KLKv$(C88lf1
z(o;2}u$f}hB90$cASa@9?)$vCU79`k3PW$)D{yS?vG85MPOwLP_b7;&^NP0LD)CRN
z(5gA)8>mWAt~;+*bV@ROEWH>dD=<<Fs=IT46%S$WYh|WdyCpm#Ldz)0C+-P-hp5@i
zK)St;t3^MOP3W>Gd$SI-G>hHd)vQ#5l{sF%<hxpF?R_{hr<cN?0O8gTj<252h2YC{
zpw{&Wl!qAsWaROgiL;j<zDFI8odY1P66T@VU<)0C2c6t$Op?zb7K#U24r39z;`5_Q
zS;*0zkTpC0VHtT^l?#4n6zehbY$~RIMIihI`fk`7KWCba+J_l{zg6k@Z(Y0LeOh5J
z1}t1%dY0alRd|cemiv-4lrxqZ*9XqU+AgKue$^<>4IFQo(0m(fZFj*f!jVg?f1$yg
z%kW+~w7$z*%Z-kLI$AI)kG6k9-QEkEwVw)&rBy#8%%C`WKbkiI@(zrmqF75O6wc?N
z@DOX<!|#AfFs6WED}3hhtfl_nB!c$O>H7slP&*$uf8wt7<3c6B3Y&WMH2ErF)^lhC
z9Xu~{2+qMT0%HQ{qcsV|A!$CVh`=Thcia=AFKu)v3T`?Tfd462<$5*<9uL^f%-i5(
zP4bx}E1Hp(=r=huK$g>Op$6WcjLyw&U?oT9oF-M2<v9*cv%&91Ra+SO5OEUof`but
zFqO6yh`JKZM5sxl#h(vR>5wS#yq3_Ly-9!|Y*P_MTQA^#mcs8_Dw}DEXG2piiq);N
zwzhln<5yu=l`o_vkF?$I!KjnY?Eq$tl-|KMngf^WQR3-o#Zw|TPO8@feDFrs09%?h
zmuuNXpL}W*KT;M=qgV`0xc@kgZpCPD`?B|*r3jL+3W_G;$=;K2VAA9;T?R_h<dM~1
zE?YP?UL$r+`Rntk7QF-ahB43psuPpEf4?gBr$tW1DC}u58ON_Zc+oaET~k14tI5WT
z-|l5JD%MY#4Ae=a0;IgX5Piorzg^r^)BbpZBt#1l&gZrWr3NzM5X~^B-PaDXFBjA`
zmz}5h$og=So|5e6x@DDU!*piA>fK@BH&Zgo&Pd|t5M1|{jmfQOSA+UDeL;ED$>c8}
zF)dcBaEQ|SLDWCJ48jJBtQI94G=grmWut_ycclk3lrt`8+opGj69AhEH0T2rSvVtN
z!%u>$O^jSXv#7u{eV~2A(6V{70$TI_Urc1LMr1*6LFY6iP4IlzOC}au!Ah+DxDOs8
zBK?KrsLYS}IJzn`s03;Iz)@%Y7Nr{;xrRx@3a%iCu8hV+TgOnDcJ@<RAyU7R!P3p4
zsiUSPV#>zmow|*ZC$UaED-|q!Fkp>uHz0ANMOaFUtz6I9HNIIZ-LzdIZcP3_q%+B1
zmJHreSfeg)qE1b#*O#y?Lb9%-_tY^@)FS~GrvDE2M0?#fuoZ1BcuvNpe$8az8Q2Gq
z*6|4&$+>e!y;iJAZrH1zU$BMoJQuFoHw-SDM=GF1N+Pd+%G+Id5-{qCBfr73yx+<0
zmc7Z-G8}97mc>^93t@Mc%*aTQM`q)c4F2|=qn~C3gVfw&bFbacDpL;MIlI-T4j+`M
z0}+!nGKpIm6x|%Wy4<pEZFT<}zbOBhdy46_cF6k%md(F2b5|Q;0#G*5T>tz$(utm|
zA#kGHsXNOs)W8{x849{yE>vm0jWYKdS#IEBF!E^Mwld1YGL84TuW-UcjE1d-WV9i*
z=<hPbGFW2on(}rSNjC?xE_IiLV;_YzDkl(&-UM}k1jnRW21yryQ72R?u!2rcNAYTI
zSy}TXe}3G%EBv&|i5-9fC(PC~bSO*D+IyXQM8yz|sENw*Ue1?KlgWZMWz%!0<z}Yg
z`o20CehJykG^XLMEsHs)W5MPlrh{5{2T8C+YK<9nQ;52Nk99iE+;-ojXB&t4TYxbM
z0bdgwz7XMfAH<6TsveDKsc>0~v9phCVlzS8Cl+Oj(&W<1>yN27?lnGOLuLu|>37;m
z-N0CkOSS6W!y;C60uZmeOZ}zjX(v(T$(Tt`sughy0mj8Y2CT2Xz7=Z2gHOL3v=8@W
z@U@qHs)x!VmNW+AjcAdMp7bohJ6N)hq1of{UBSUhX~<Y{RbBJ%{VZu8oa`@S#Jj2x
zB2QmdWMeend<-G*n8v0<9xc&_iXzNH9M2}Cjc81xoDh{P>}y;tBccoNrDHXM9RZOX
zA#)RV6rMXZBS*U$X0Wcdai+(z+WQ#Cz@8@`fTq4ZQw1sUNC@dGYI~u!)fmJtqZ0f$
zTUL^AfON|dqG=CINv_i&#Y<|*UfV<a#hEnu#q_zXk1ZKkI7jF#VFDT3rq#p3x&R+*
z-qwa=(LEV?Tl;<tby0c^cC_MLLo8{=#3MQ^V-l@LFm0=%3|gDOimZ1VmLbE-!h5YK
zK247Mvn3Yu=u_;CP*d~x2$tvJ{^;oyZ0`m>c*_ykfeALEV!EpelsWoTT9$_n3_Nvn
z!n@dF+3Db1J>?PtQ`aTer_8p=;;`^Nv3Be{8!VmQ2m?N^<#bL;#HOuA-?t@L7)!b@
z7aLN3o8)rHVu%rCkXshKd#}ybd_7E9ez<dhv|6#>;IVUj$jJIMgauwCEdF%pPkA(D
z9+n-MoQ4*07C)cD%S;kbj(_^HaF>Dd;l@g1YG$;7tqJABdzJF&EOg*q5O8)CK!_|1
z2;lF`ce}+PPIzr1#g3HBxtM00Es^UIRzD5=kuxmHpN7L?YJ+Zj`W;m`0{QTttP*g2
z;)m<}=r17v^(JF6y6OuKSXTf&=SU;P!}XSc4npcRo2>=uh4DGJYg<+<m@PS#pgz8Q
z7Nv8D<9IA8$WfNy|D!ffP1c~DI!F%=yqlz$QD9oNcAJ|97U;*(r;WvFIOJ<lH9e##
z9z=6B(JZ`&)xAFL)qC~vqy*~`5>6(jmYz)k);whz@Fo>Jn9K*->fBlgZoLkDe1JGI
z@#wwbpyJ`f{%soXq5))$Hd&S`-=}T<gdy}tUGLVN7GIgqAXZA+qZLtXV;YQV6=l-r
z<|zZ~m?Dd;>Z;6<clUfdCc}cRUJ?Bd1p`k^|4AbTPaNs3lI-L!B<&7ztOGpwrdsr;
z+?z5%Wouik8fgzqqh;>?MA)q}a-Pc*1f}?bvcWcP$iC(Sxxbjy6+tlf3DInn#c;8j
zJpb+JF&ssWDO+XgO<hfWXg2jZ`4j>{?JrluEd0XD?Qu)Czq}LncZ$QK+3O)t8t`q^
znsah`g2y>G*-%5WPryyHe}%hANpZ4GYd%TyYbds4VYp_bl+1jMKo5aY<y~?#(~+*p
zvh$lABq7}lvHe+oN&Amr8ld#y`-%C!OJsX}f!2qFZAl;{l~$}@ow?RPH|BaD_N`1h
zXFTPt8lHVwPb<j}u1pO7T*WHZLy8q#=6}JbJanWD$Q=le)46dfO3e~27rGZo*Eb+3
zEsM}j7J&2yBMcG?{<!Ni*G4|!(+98!iE;3ljIlYJ(zrAkn($EZI*f>LJqlsTD=)_f
zVI0Q1IyC9+>=j;lDa1{iqH-39I;*p?6M;60jy-Vgnx7|sxP0B!-}>U~8dp#ubP}0$
zl*bU3?5_mS?EMSL5PE)sbX;^F?0@HCncB4S-L#BWXvi*3UO|fC0>G8?9-9OHwz3bl
z{LDG3EP)(G5MfYIH(eAp`q4ZPa}si5D0J^9?k1cIL)&~Ae;~GUutuaTRPLF~=gCcK
z@j%>sY@Ti+JOOo@;`n~r9%w>Y)=`79bgU081ago=NJ9}{HQB1&e~)XxetN<m-*ifD
z!!Fi`xk`hh8ap}sXM<rqT(@21+a&{T#G(wyO?DDKUS5}8D7A=igI}a%BQN*r^#RHz
zHINB!cvx=vVtt~^Hv4U>JDn@eTAXgdTTSAJOFf4~v7&K0gn)r!cyB8?2FVeWuPBs!
zeVcMjO(;T_yO|_Jg4eA56x&%_=Xf2PtnO{t9{Ne9O}xivEmcmhvWmh>|Mq%fGBd~F
z#)LlytW$EH6{WG&iKK91YhMPK`6F2j0ljVPV<H(a>-M35&aAC7HDhD`Y7)zkLmg4c
z?CxzB3S#Hup5|O(-F$w;H-Cy0aXV2L0<pc%@!Et%Eaf7aKdDnW4-E02xo=cT=;is&
ztR%iMyAg-dEqfp34U0yIe96o`TUG)21(7gAdhCvN`nL>%v^`2Za#lvkO<}M{9*~NG
zTBIfGId&h5nE5FH^DiIu)_x_aE&H_iM2lg#Ya7y@%20_TI+op#0ie6?`)2zNq${Zw
zyd}Ee&SAX<Ik3ii@B4zabI>fI&ni3R`-uGyksSN`fe#w=z9FGm`A~+`zoddAxKx9j
zXNAYxIiu*wruHwq$Kt-Z4`G-#rc17E*Qo&-0eT}p6{Cl9vms4|k`An!4OZ)RL%)eU
zWc#gCKuaOw$A~<)^rh;!;O}0i)wzPc0XA-Z(-faPT-FCeJAE@5G?T)Pw1Z~eJ9uSX
zF4F*=3AuNksqBT0O+jc{ACqwOWS-hJ&2XP2mL^5maxV&Z(A{ta$$~kiCS5WNJ^;19
zL+LBE1Pk9b(<{OZ_qT*RdNE<!Ib>DQ6|6rf2x#Hp`zd`I6pyX3kU6lb-kj#6o-9g9
zXi5EKfPs}I-34iuMDkHwITn$d3M%OkL0=8i`Urpa^L<gI8h$`jqW1x){M9LPXMe%V
z(aq0e*rnIaXL&t1@hF-)JNZuwKFIqP(#Xxr66$yg)DMJP!^hlqqUz+vWR>7BU}&bo
zM2%E~d~H^pLtEyNGv~Z<Qa@>y>i8;5Txx<?V-xg?Gs^H6g#6BGiZ3)NPfk}A0>V=^
zK~p$x$5HYRkTV77EvckjhzHu-F_84P1{bM8vh4B1sc$zlich?Kp1MDy#q`&JIdn8G
zS|k3>BsaRlp>6HR8AOsobY;(*&?){%N6Vzy_0cC-9(l<td6jO$ugqKkZxmFzo6Q`}
zUDIi-<e9XKr|>!EV`FkWeGCkhZDI8%^HNfBgKq`63e^%V43LBeseWnx<tyKcmR2TI
z?u2J`S1D``{nChkC^_uVCJ%B>nkYf&+`v&X@In{0=va3;G%|8MtMi@m|Fuj%D9!8T
z+OC|K(doMly4cslCkLI%%YfYUtyxE7ijeEV?ef^;E>;Xqmp+tqX|s6;LGav+%Jl)4
z{=o-@Vc(;~ec6bM3L3d4_e9*25cEjj9iKqReDehkPsfTGa7y^hlRbHjyUOrddsBNC
z5nob}L@<1zXtE;_t@NA!&)L#ZRzftSY1&&?nM%CaYE{P48L(V1vb8<Eyb{fL^@eW9
zFDDjiU5#n~g3imgD<XzyKlI(@3axZ!4)@3YCH<-K=phjw-sx1QqhEJ>1J;<p_LMiy
zTY5E0GSNP^x_0%K;%hgLkt<5=UxY5)++I*TeoX@Kqs>lPjGJmv*&?P9?{BJPg#Q-_
zxbm{8#&jG&=i`GG6_XPSI^FCSSPv`K%wv^xvo1Iu#Ia^U_YZtRnU4TR0G;KOC7XaM
z=p1vw<85l33AtJ2_r-nk4U}lc-N;6$$N&E5(Ko9|xQ^rJ(}VCQ5!fTDM+WW}LmG&5
zic$RRnU+tdOfmif(S5+<-n^I7?BVD8rW_JZV*kp)l{<AsJLG?8F^V6qiI8aDN|4Np
zW@<hRynnQCbr^@Pp-tHo6D7ncjtl5~-b7aqfuEhF9;?@2Cc;T0CK+ZduOGB5MiixW
zVsH-fh87H@x5oEY7~&16cl|$ZQ|80g2IE=_0NQBv`xG&>IjQ8x-6Bs4<S>bF9cITS
zWPdQ8!v5j%Rr%PFQ(eogN#iW7P5!K(ux=lo)jb$HO_4Bbr8xpOW#YQ6Q>7@|84FaR
zOi<;uX2J3c>@@yH&F5m6L3|LvG<S(3NrO76vWrzVNB$r~!#?R9`5031RYz(~-1^B8
zjduv%CV8&eIm?c@57N@DC8Y+a{a9b%?2#(6Qn|>gm91bDURrCudLUUfXw)wvT9|yT
z44i_@u7oQnQ!-A5xrs`J1eazI=SFERO}AJe5F^|~o_i=0RH0a(+C56@bML-k27_lh
z7d*9u>6S{x?0#7C!Zp7uMOd<MGJTl(BGl4WJX8Er4AS*Qw)aBh%`qd(AKZ({e8Sx1
zXVcP*+9&526%kuDYtpmUvAXzR1!CQjb#wnNqKSE-_)<ylSuwqZ-PryotIk7p444C1
z6463@0h@~rbRKW?n+FnI=SF@iLC-FHkT+_Edal0B69+sYGm3v;%ricD(oi@ot3u$k
z@Ql+-*n^ITD(#pqJpx>>cc6{1)uMcbL<Qk+nuM?<x6^^JLZ(UiCy&1$(=BePL%lrB
zbQpvvdsJd&W~W$)=}b#{SyUKKg;Uo{ZXZZ;Oh<X(JrUM?xxP12;u;jz<S;bZI+V))
zH56<hp*pG??!IM+U)AK(h<!-EtQqquR+UJ`r*DajB1XXXYP(6O6smKSLHb~c=Oe1E
zvvN}D&i~xiu_44u8`VZ(!oKT1r9;fQ@AJ|UGGO|gR4TJjd)5T3Z^xqM?}daxOb^))
zrYYZ)9y%G<aZI<eB9Cg#>r-dTHR9|V+;v-5g>5YZ>iQqQO%-nFp*9gNInCqKlJaT3
zz`!+W=a<<{N^A2V+Ex)yZ5uh{5dja>wBl~bnud|sgCjHOr&{R$bnXL90t*<%+&W1q
z>&r3rm4`Bicq>C!8VJ|>ZI|llK#=9}TW-uHf%#gy(NdX{?L*tQo%b9+e`VfK^&z)1
z481_CavrP9X;MiadB938R!FQ#2eOs=nd^8PPwz7!5x0)H`iks!&geyi7UjUuxNaHS
zrCUHaW;1LQkkWBS=b@Y};1-Yw%?^zeOappNY;T7Kh7}3_(+M2xAS`}ZDV3cIRmr3E
z<=wm9S0g<5r=y!RRmka~7@TFBPje!PnXuzG6IHKi0Xj=GO<Vj&RzHb7<1O71BB}M9
z6ekne6Pv246q>}KIlr{;A5(rur+TjZBU=khJi#%b-kLLDI^%?c*h3vqCACUxY3*LS
z+=Z1<kNMpM{)fD$(ee%CJo`~gD~%}t2+gvXURgKu<wG{)q(EKkj1s0Urqz{FxF`fN
zY`u9}+BAzs4-Emb>-ZqG>7x-=^(2_-PXWwD8#&S~ZC6DtOkeB5RnXhg(zH9`t_&N=
zCQe*2_NlT!yJ_9n<t>Ul_=$82=w==@31@*;IXwlg?=w8UpcQbhuvM2AqU}>TS8nJC
z>yS*q4*BG}slN4!Jba!7lx1#e=M9$gF9_ZInG0<18Fvdv*}u!yR2i#|GMgr8<LYSm
z<J@KKnTfX9s7tN1lNP8%h`rt&-YbS1>?Ia|Yd7;vq{9DCgWWdQV(ZsZDT(he8RKrM
zj!9_7TMWxo51fW^W=;F;J=+l$foz*JNYB5M0;pg6Itd7MHA&j|1^<GJaqKfK;gg^n
z#i?6Tt6^PACLsS;iWF;edIYLqph;FSEaez6aBO4zz{zdgn#0rjyx8oRWlPj{j6rg<
zJq@Wuk)q98Sf>BxGz82w7FYhm&fdUQW_2!V2N>pWlnIj=+ltC)2BehT(J3oe6~L+t
zYDKR2YYDEmsRyD22^p*a@y}bXold%vRS7R-rmoD*6ZRV5JiM6PK{sR;S4-Gd^BbFp
zqGgjtPI=i2>UrDR0ZS0n%_dbnKR-JH3$7K!ow8@*e+-fQX04}!xx7PD)ssdsoh{&L
z{x|a=E}bt6k$6i!tu1_E>MKIn)f3nL!%VdeA`Fyv(ORk})4ntzoH4S;1+H}s&2d&$
z=wE`zx*O-Z78TTS!SQVW_9I75n5lU8J5>f7)#!!x=vWyQpJ=6Gd$z#BTb$!F<%RW}
z67EXid?BGQn|1Yf*w;k*KjUR<n9;32<8ZVRzm(x^`d%XQZ>fMpgO!%RnRV-0h*_%d
zIo&a7v32aYp4u<P^G6mH81G55+HiT?p4xsO-V<F_{_qoxHKQTA)F{Q!NR|)2fl9;7
zFU!2JaUsMBIN8tFPQjs${z7ucBiR<=tV}LN-CmCrk6ueJ(@xWmAY6q~9bT|Y8LBt@
zFPb;DcN|gDPDLK)h9`7rxLr#GY7BG@_Y7M*V4A(4bHJhnLJ<IUu>aD#W8!a-Tfxt}
zf8d}y77@|L^f!AI4p_4Fal2sR>cugO@@IvPCTd&;xM+B>`_lGJhoBIKxpK~HTKrKX
zS7PzXbg51rpBZQ8{-<8rtSqF%0SWt=v)W#Hy_%OLNcjY=l~kYfjBIBAx>kC=aQTg|
zzNJzHD??q3p;ai;@%G`%2d6C&eX8@?zjWpkfTA#%a-QPFbg4Mkj0<_c85XuLB}1di
z1rc3AV%lQX8x({E%o`VrJWg51*}rOw>62m7m#6)#KOig*^~;S|v=V%=mJ>nj3L1SN
z$z*KbKteWh{3x_Vw5mDe!4ZA9SGiM&E=+xh!@1gZB2>_?Fq+8TzJ%Q^jNY-ECZg#3
z`IGM)Uq<*Tip>XTtHO$3$`Hnqkp&IRqYhD?o)Tp2Q$dr}v|7K>uFr&Q=&zeeD|^Ux
zq_Sn&!PTk8aKD!Qfj(f#NyZtW%v>)H!P8Hvl^4MK05CQDeO%k%CIjSXXS>WY5|KI{
z@8flz*b`a<?}nVsRNg+xRzen2qiGS__oPQZjbO-e8-;G|loQ%OaoQdg7Whd9_Vi9T
zKlSX(7V1i=;86Ax{%>`O$}4!0IOohby%8WB1yZGCf8Q#w%#edrdG^-`FPR;lAgj+;
zxX4k<lY~(z%4k94!)=BNsC<(9rX_g{;iSy4e0qUGDxyy!jD1R}v)U>as1E8F!qmMg
zPH_xdn%m|kc=V~~0_T@BF2ZU%d>fM;G4TeZ_;4Sr2Hj(p(CI2kr-sh^vu4L00jyuT
zdvqICol-}&hetZcI1PXOa`hhTO5R2ITI@t@ZJ<3v^qP5m?ZEF2->i8xsr5$SaAeQ6
z{jcPmc=S6=i7cLgJ4L=Z$weaur0+eDSQyA$#DK@)tv)@Ck)&(J`f?D0_bL&HQ?5&e
zbR|Up;0c#6@92qr&W*0NyPU!@s$@?6#N!^jsfX=pc*oGIp(m_^h=%X48zdsGI>C;8
z+vGU@kkk_5EgP7pPYgPgDa=xo^3(s3*jhY4nBsXchh%3MBo&3Ltn5Ygw>`ILK1`t)
zk$X3XXFXTy8`KD5ab)=2%n{6|P%fo`Edv<eUOa#DVX!fsAQc8+a%ANL^|eKrN}Qo7
z5uM%p`?zn1J7POAF5Vke`Z;?XF{bsiv$`N;JK<MRk-c45@P(eJFTHwmcY;uZGsYq8
zn_>d)g3^;H$1WNdX?ZybdK_;@98f?N1Ip9a6Dxu?USl^c5bii)UG->;EVJz?vufXA
z=ztNPCnZ>1M~5ToM0k1Z+_|n*MfulxhayGD6do}rw$E8W;6SgOqVr5$zr%~pdWFcr
zB*W}^N&QWx1h(>1p4PAKwD|bDF%MRuHySL7_}83ixtr-Ip|mIC@#1p4-`X7Ia=c0p
znU&aW2UKF|;>w%pS&dg$a#IIt1r6_1;eZVsu3M{5-n5Lw8yuNUAHAEz@DQPp<?bj6
ztDWdYI9y!xyIA3u9QQ>Iuf#nLYD}`SVY%n_*iU}H$1s6f!3jKlL5-lp;UsR#c`dm6
z_#d=P-+hW=&rn)OELrr^S+TW2Z+Y-vncjvV<vOAT;R#qTzXeNP2rSDcsM_AhJG$K=
z<9^`WFOr~VK?8Op8;_0>anJ!~vZbYsD89_)#ioS1u>6>h>D3~kP<Sm?DI<88mTu@q
zY1fxvXZW6J;dFqSm>8GS5tLoVifufcbjv9lG;YJ%s=|u=U>eucG3v-0Iu$RIpsIHx
zKjk8Sv`b=1j)yOF37WS`J&fnjeG>Rim?C4N{euUI6cqB@ELD*lvOwioHd`Zir0{4}
z_Ez&5a`m#mC{aB+y$Ww7{Z5vzcH~wvC;C^C(_t@1ysW5MA~)cJ(eWbAn$ic+lV-aE
z0x#+A!fV_q$VSS@myGu7vYv~sy-+kV_o7*~QAz|fp4<K;J=Jt0(q7hix`t5rMN^y!
zL#8O0O-9I@T)kw_DQAht_lL+aw!CBg1C#Gaqt!>FaA<h^kyp9y?K)zCMwecs<PPmw
zcDpsv-HrsxyYVL*&Mez=iwI+2`50YTd^xwkG|aFIT9R{j{}|@lU*0unb%_HK$n91e
zT<Ul|t#eioEvpBBh@or;vo-DVjdgptEK*Uj{JH(4BR5aN(M!uI7>(XjE1W~DV&?yc
ziK_~zvTM462m;a|t<ojkE#2KI-QE4t-5|{&4$>jr-QC@t2c%2ze~#bHf5jDS_OtiQ
zthHwC@y!T~mWq=oOl7973JdK<QM2W4RPsCBjORO$%w18Bv;I%x>dP4VY+l{>4bkj-
zr*cK?8-GjV*lXzKSigjPh+39w%}M@V1Rpi1GNqz|u~+pMX4`f!evnFPHv!Md^F@(Z
zM;Q^Xs`NJp+mQMPPd5B3i-YJ|N=%kMqNyS{RDaSTDL${P5S^0j;tBk>7ho>Dv_M)-
z&qjh38af8c{0|*2Ca$!_>G-=?YSZ}|e*_@r>4>|mBFgTo$?EU&H82NuV~=AMFbBw(
zi}5R{m4d>BP2~$a0iqY63tTv6vmOufnt4NUO3b`Q8lzaxGQ#w`P&)?;${Lgv#kb#<
zr$Cv`?6=Vo_u1cL&_puBJ4$t*ZY&8JWnCbF_}{tvcL31wNr$xCv-BtJ;0v=(^p{$2
z`R+H;&)-7Rv^ryW5??C9IlJLTO|j!5Yho!)`T6F<2oerqlS470hn9sEd_(Zgb?8y$
zoUek+K3Htew7FBnvB<kLp)Dm5^!_z2M9YYcVy&>3_*b6!TQ(bAx;|zlelNPqKf(VK
z8n{q)MYc2(dn!KxuR|)%`FOU4=9O~=vjp$RK-fG}#Jm(*gk)0?v07Ng9X&+zvpD4n
zI6k@>HyloD48z|rj>A_S5%%%>*phf26@Q8%cr1ul`&(7*?zw>g@H-tYn$R%4hsu9T
zb^9?6UaPyNjy>s{+h<aNUIfbv;VX(`w_k_yRq=a+T$J}bScWf64vMHo;$p9)e~^*{
zxfhQpfo6+G#A}E0QW6|1jax6ln|*qY_`x`CaFBcj)$Sur?4Yu`K<BGbQCDusmZmaV
zwUg)@FPjn%wRc)01mCNqu-e>oy}WvUQprmxtP(R)W-iWfWgf7MI3LEIK4`ujBcDx@
zhE=GrzyGEq!5hLTf4w1{V*1y~(zy}4Pc<o?l2B9Vf(iGQ91_8UnKlxa8^uiz<r6h1
z&0~!T?#|V&X73AKJtL6k4Y{;pm}<#pXR%oLV?O*EEIh<I87NgD4rTyzxmz{O6wwT&
z>3S*ay7b^J`7dY1t_QACze0AfohOhm&3zY3T))Qor4%~~LmE7^$LRl)y%iF+ROuxb
zfYy_@m|lEMW1A3Yy=2$ym{Ly4czNc!M1X4~&IEye9b4}C2MyLa7QFY%)GyxRAeW%y
z;7_;lVp(>X=~H8CPE5jEhY4E2wyRsov{|7(60P~@#Ya=;IPswVH9&{9TX^>-Ca#2o
zAsk$`vg!j${%*D+^PWyW$LV-dvPwlKiCS58eSO!}GjKUpOhMslW>)eysp_}afZ*@?
zoaF7%J6RzvKNVD%Sg&WWYI~bV<?*>^r1G8CcP!`jQOGQ@S~4{pzIka73h44lI|l4c
z{at8M4>;4rg=uX%w2&=Pbc53QedGhjpM|T_*$Mhog`uJ1u~y>eWixQq#HSz-o3@lF
z4lZKL6eg1U`#fiASN#2Dx%SvN_JJ{Vx`Id#PJQ8g+Sc*u@(U#*B~^rMwLQpt`F{Tk
zX}-wLm<~aTJa&rf&!1Usnf{bK1;tOF)w4{iWno8aaxFEF+0kLCIKg^Xj5t9H85{2j
zAAjBimSZ12P3g|%HhRQat**FZ>`3LVaSVPEkGb48E{?Fka!uMzx;0H*BK_q=z9$Iz
z8zI9`9ZN9lPtj@aE%RA+@Zd@=orsESdQPD!K;`374k=2NMh=Wl$Ku1a?W!!NhK=CT
zh@C+oro$3u19Bocp+Hi&wv$Rz^h2go9)Sga=5*%;^3cHPm8<hPABWGfI}hi~`slog
z=XAdXB6q!0aIrk&xf-)a`d-p3g9`O({#wm+?TyE|;^bh<(iC~DjrT94Qy6jS9=)x4
z#yQ}Fq+Q<xb1Ot%tUr-2D3i7L7f1gP0wsN_ym3TgrF5-SOSn+?`+Vu2No(C2v!l}x
zkJapda-<xU;aIlU=~z&&o^@qm3s!V4S3n%c(M_aY0y2mhJlC3TL3yiRP7;uvi*!o|
zG4i`^uRce6jPhRb2})s4ih^1z&=H9H6!mi(iiNTep#j%TpXs?%E+tf=)LoH41T&_x
zA%3~}3em>2<dD~Yp77NNKzR6eWbDc&k$|fbs}dDBM{l2R_#bmBZ^Pl=@8blbuHMxO
zDGlxtR<$f;*b0uI^X)hmC7+70i#-K5V^xcfh}?xzkcYae2S>?o8|CeE)?Bi-n&f_5
zVPn>MsRTYR4)85PjX@bt@rJ~uy^E`}q-Mxlh^u3_Y<8$qZjct=;a+OQ6aHhV6OT1s
zHpcg&g{K<!>BS~5H2^S@^<QHt^I_32Fu=Cdbcd(h#vdqw8n=M0iHwu?v(`>g6JOTQ
z)zR|)+{g(hJWu<)tvxL1FwXPnzaJz4v$<e{wLYjSy<W{DVCU(9WKKyXpj7YX5~rcF
zyACj2f|-Q=x%s*SsJ=){1ImftNte9TWKIxW(<q=kal{ScI>gWe@Pdh~G9PP67U6jT
zKdFBulLO&C2BS<<w0O9Mh0baBo#>T*Tc7R1J1<g7h>+Y(M7}b2TB|qktjaQHdYx0^
z+{ROorgNfv9aG0@doKsd(CUdvKsmCVW6Aolu!4JZ-l~nq&#6LLW#YiaM`3sVLt#lu
zg{XYYG*wfSpS4b?1wXFxEsdVBW8~QDa315knC^;<WZ(!&NUgEyp7rt*J@v<77Tfe<
zf+_s!-)>Hvi`rD7Y&Rh>VvDBoAIZATiXnf<h4{>$%-V=Xeo@Ymw}Lma?yiW6zW*z!
zp<E?3^|z^d1aebX-FUsSHY}lLYN{>9d}3R%(ihSK@0Nr6s}OOV9bGuqMTsL8_lg!G
z;8LJcmdkg+v9o$8IyKodtnZ-_g(|ff#5&nId*S9rW&6m>p{zxrbTsSoe({Ghadf^}
z!ZGRO6Qz7*=_6UJF9l|i`#JNnayRJ<ve`-B&7Ksj3nzxFRJ_gERTxLHR@~U!{e+zl
z69VZ`IxpM~L+ZZ$P~Z%&?OD96hR>98Hs7(QbDE<#!<2wwZn|rZ2*2|q^qKaK<$>Jl
zMU1wtWqf7)CEC&+xKhp0bccjlVJY23-a^$f;x_X(bzlR~RMJU~x-L{}TL*2Jg~wtK
zCltK$O!5tM5dI5Rw<u;%^Iw*^#%rYWO@obIMxwy-jS9xxRFtG7iP<gxOGp>np@(Cu
z?>pEh$p0)gI?AhLv<5%3{}%ID^D}w>5Odju_&*HLE;ICgpSK9;lewDNAEjv*m`AqL
z%1f?J?>-{(gCb^s)ULIjytAK<*vu~}YUd043I)3CNuSKs8j|414tv#2^MiQBaCdLx
zmsB~HKuxW^t%<4G`lLr7W?Ld(P(6Zcxgx_B!m`XGA619!d==L^>2(Ekj>O#B_l&rb
zved=D2EjVZ)w}yrR-Sp$M+eguxM~vJ#_FGYZQ~EIW9fW}dMnbJ_sfy|DkCH0(>4pW
z(bvQu?Kl4(c-ix_<qLe!@;jZ{VP@`)a8GqrKTyph{aZXTJehPgI`+nqmAPK#UY(n&
zlK@^ub!m=$4~m7i{-8Q2@chCyTIdJc!Myb&<Z&t`zZPnvHfv8CdwS49vM)ef!HdMo
zOpqgX5kbMc5E$`HER_KJ;#_vtHuWBj+S<S#AW#)4%tY|u4b%C>8exB3O^K-C9$JXG
zGe*)f!%9h~!OR(+EE-B4wBd#wk0obRki6RV{O@#C(Be9nH1GGW{jF@@Hlqh)6HvD7
zWHZ_<^q*1#ZZUtCV5aq~UE3M>s_~w+CjQ4!{i>OL+~+k4C`D&3M(y=0JglZX)qXVu
z+h!8UsZ^q#B!0N{Gc<c&s@|Ss+}PEo9{QC1yT-%GwX|GBn08;ZaVHZRgI1|KZSRJ0
zvR>t|4<TE65xl>TsZ)jIRaD#S1Sy}`axi>omA9jBSZ3SRqGP49dfb7GARqH1oLh>m
zJ0|-x)vN!n2IG;@q`6CBS<m=fSya<UI8WCb{*5ma0k1y<<b=j$iK*limGB!|A0o#{
zd=oUNxA8Z;PTt#7R66ohh3q$iSg)+8h(5?|;!oZn4%HyGy&ghok>Ef~|IW=&t21yT
z^vR_~)ne{E^gvw=#>=BvDecj^HC-6J<5{B9N9%>OCGd_DGabMeTary|XWB#)D9mfY
ziYUvfawQwx5bip&8Osu436efH<@KHneB=3(rSP`&&824gR}i+<TbEAD(Yk$2wFpxs
z_TEo}P$eBt+4%?8L@i0&GPM47>w^k{aODwjt@&EYaf7M6=6b$NCLg)_lw%kd2I2nj
z*V{_Thtu2d+z_{HE3f7;5o9N~a<lWis@x+J?=IBeD*yDsnqO7?E+~a5QB?^k#dK6w
zW6zx;cpizxjN@S}F$H#{^oHbKn!{sxdT$?5PM?dA$dW+K2=_A2C8hhzxal*HIK-oe
z6KhEWQv$odk%C3xDTGmEV61W+X_K;(!6$*pUYUNW^uRx8S3IgYm&T^jUoBdD?5~O}
zmsynR)of6X-!M^139!PscM1GP*)?t;fH50POr=`h3aJJ8GD_N@VeXk9`~@KBA@l!o
zD+iY~sf732ed2Ej&UU+U5sjhe`Ku&B&;Nzskyc}!pS;!V{8e9}6wkdtIW$V{7%uQ*
z%;tAQW!gUCjJ(jm7<`TvW`f=D;Xa_Vk(dp-rBd0r<I`x%PGNMy^gHCgSuD^8;F5UN
zjU^H@<ucGBaG0`_c`rO6uz9(rssARH!Fv-^C=oQT9X=W}G@Y_kMbZnd&b(n(k$6Ma
zA+aOd4~=6{Uk5D&IgNARp;AdbJ_MDxD3>2F`X;Emi@N(J!s|z5t<TIdCAR*o&?>z0
z4;9?{Gz5Ie9&L84ZmCp8_+i7>rpQU!>kdMipalqq#sV#K$bRQE!<9C8J}Dnlkn%(p
zpBxCV{Ie{f8g-qb!l8{|yN+n^uu(8Ur^q-~Lkr98@w7#i=x5_eHFoy)=>HUPzU#?#
zv}E|WnjWkj5SI!5At`Jl8x`Coj+;*h5A2mujZ7#RV=iE+6t?Y}um_IER#hrIobWi+
z<n6+BzpEpn9KmP>aI%dDz`du^Ok63D77!dWn33M~rM66VCYajOlfwzG>OPj*J4g`_
zt4;5VZ^%R7#@VEHIC1ztE&;O$)MrPRqYSjc{U#-7F~9v}o>kfrsAMoy-LPGKg_RZ)
zw8YFpC4v~EBUPRM1CE=@8443rWkWOA9GgbP$HRS4EIyX;$VZp84{wul9p{EbjcEIR
z{&{&Vaalg~P+PQyGx*FrE8g&tT;WtW2}-EE3x<rI3WTQ+k6o>-dpPmA?~ehEuG?l%
zev#_TMGbNENWSs1+bP2_?+Fk@yiS2E)RQiyP)obY?s?NhlNuE5ez;2g7qYSn#x3k!
zR80@}*0@4dkr7$hQ5-v*Q_Q4>OBCZ15rc(RT6j^knuXH<`NO@%=I))pbg;XZHCzF<
z|67ZK;N#=wGDU&#7g6;`ry#z~4Jh@tq2!k<wb74!V8z$@(`^Hij+i25EoK5)wkGcI
z0b%%=$#GudW^-Md5oZo{UdWE`itDD(>}nS1Ed1^=4lnKbF|4Q{+pXvdM^WnLUPF*)
zCZV)Q#l7WPV8m>JdY`)bYk7$I0FW`TrvMoPHYOcbJ|jFX#xak7yC5-MH7B^&yuIHR
zHps{KH<v3hQB`ab5?1e4pF2D1xp8TX4U!Ni{4-Qpty4IvQon(_W(6KmsH<w#Pc?nJ
z{@;hRTBLf*MpuLhHp=I&BeAm*zXVnbtc8ju>D=Fv*c>HMq%EjPUn6bg3C!kP2xD%_
z&~(_<cZ?5xu{T?WCBKGD4;EW1yAyt9b<IcKig#sw0uD~r_c226XiRzyu;a5@D``f;
zeBwFqhGJK6?<p2f8>R|j1EQ-%51{JmIN<klCp`m<Lz4B_VL`MfvtkIZ*{;-hMJ|GO
z&Xvv4X{uQK`;0ne0@+g|RdlehX=#g4x`m3Ov*(mY7Xa{%wYQzv%;GPg)Fv^w!Du@d
z)vAU-$ZRldk;=z|LO^#aFu^#_W-wx8fVl553{VsR7VP7pHav@Wuchcs$n>OR7+Rah
zoFuVP@IjmUB$y&)@PM!%T%{LFUBFE<CyPA-Q(MX%bnrhpXwlX!6-41H@0U%mMi9A@
zX}56m1P~5uUiCJ9!+yf*9ZaO$Nw<~4+G4HbR3#8=Vb!kjyC$>5ofs{P%f^0JM0UG%
z4?O#rG{Z+_@)-%?Qq9}}-3~fWM;m7wYA4YsY5^NTvy<?0qDi74-;4%&EK`a?I-Meq
zA2LQ2fVyf4%U=X<2yGlf4HycVmjmToc3qigAMk06VdY%YS^kKr)OoO4UYgD9K#1r!
z<Y#ui{*C@9sBSJDrPDyqYFbp+^sk^xM=7Y2t({q1-UWZOC!>!|gniUUr*e1MP)Vf;
z;%lXWmW7z%s3h8>tSu{x;Ru39gGQ1oZQZCZes*q{2VfrEywJVr$eGawk_0~oeNJak
zu}Xr&K@<APEYugF<KR4xN`RzUR-vLouqJd1%Lo0J-9C7F_r!A=xsgTwoZeRHkVRa&
z(FgDYVaq@KAeB%{+q&|5Lp3uTzz$U?j22}jhe5S}s_M6kVIuStsM4s}=uri!I>$bw
zeGL0Y`X(cLzo$?&LDH`73U;N)+%zWsU@D8*8)6U8)Is{wyIX<Ey)qU{RS?nig*aQ&
zIVo4v776J8(N&CNvK+#~SA4GHGxYrPn%1`kt31ht@?@MZk6a?<l5a^3#Xsdl-Fb$4
z!LTXIHzr2XV&EIh+vr|^?`070u7j91V({%<VrZNc8~uf8AQBh*k7ZD#U?n@yN~|;n
zn;Q74IF<yk%g%(M9o#MPZ+>cG`ZZkICpIrV;4_7ta=IxZ>mcuFr#;m4;CoRG$8E>c
zp|d|w5|&XTw%WGzuaPw9-bh8WY`e#^Umj<H9D?_N&f@L(dH^Lw_^f2DN`;Q=+LG_W
z^x#o^am-Zh9d-K+RRaf;;c!lq+Hmi?1GF6FC*<WTo&goT%!W<USUHzDsfKDq{Lr+F
z`{>0CdZ`dGF^SmFv%ZSwE>|x?q~zTQwUivy^@KT}@~<Qr!Cxi%{(*$_g!OR#{&m<V
z5u4im_D{Cl8MmSK0uh7cjFsoLcM~(yxKx4~DgJQ2<=DEht_kbk<pY9PY)5PMHD@CW
zX^Qbx44AQM;vjFb!!bCJ=&A~g8;Sp}f!OFRWn+v8J|^0UqNtI(<M%(dLzSwgG(8zf
ztfuUpM)n&6(N?V=A#>qj5LidqZbvSD7a@WuW1Hf>l*;g^3UN)JvHD};-5hNFiEu{N
zbyGAEqgs%H2Q)<KywRmUHi-;DA^4nIo9n4U!_Nd8ly_9Dp-J!9Wamj%^s_G3px?6I
z?p$^(V<i0(GaFK}j=+~niP?>Yyd;YcF_jcgS1ASUO^GP{S;#5s>GiCpM8Zg8NNSp_
zHacddcE&3%1QGQMUN(|qw6F^*cs%;naI%{ilP`#uYs43F)neVU>lucuz{G-%X26Vv
zqu(OMTT5zM=HRsDnM4rbP1wK9W;8aq5N1WEtfqrCYA<Rah9JjHRog+BS*&s?tY${_
zzs?69I`MWTlTkU;A64_Ct4!m0!i2pideiU#5wa#>g|o`|N|=xit{i0qduU`ib8;Wr
zIZb4z&1rk5+5-cXQH2%N3S242=tlG?EQG1N%=9*XG$|a#NAIyarnmD<Qmz;nxhN6$
z2+9HKVrkQ|Fm*)(zIC&!Rivm*xTKoX3efl%k99tOi|F%)pr3nm1KHcFRKNKXkmbO?
zpz%BY&7W>Ap7zV2gNr}-2m#!{!(4@DIsXoFWI$prsrHM+6NYrf{o~l_c`>nC7q@yb
zr9&V{A>e<s>3ieC`w=X_%P2oiPyJzmW(aCla9B;4e?in*Cxa+0FU~k1Zm3MLx*#}q
zh(m8kUAnF-RVnsYV&#<SUq$4^b8$H#M0<)2kIix%B@esyo}tUn54e}Mgc<}1OnA&f
zVt@MwWvMQWnRH`BTA_>f!B^wH#I-CO`8_c%7?%b0Wp1%#yhaVCcYfJ|Se(9+celo!
zm3x4Pja7yEZs8hH0m)t0*AdVeyI8MVd5bc@FCo?#+f<a_Xt!e4t0)+T@kL^I;VCV{
zq?Ov8%xr2gxA;%W8;otO2LE_f*{MQx)WH>8+!J-LiY+%`uk0114t}mX52){C%uUq&
zKVg^qgvoVj>U|%fipW1bk#$6U%@<1m889@vTFX9Kjg)G}=;a|4ucwv|V3Ro1d=w(1
zPIrj^DYLV~!zp{wubCVqbc1C|9{sg~$4=?DWk$o8YRpmMVYB(^x~@$&exNYT=41Xo
zpzx#Df}5X{>dv<xW!hH<Iu&KzIn&SW%-q~bgL2sK&qc)qt<T4=CNsPhZqhe4LH~4!
z3;W4`I>db=phIi`sr<j5cz2_r&X_l6S2${(-Z)YM;jE6n0KEK_PVQ-LMDQbah|Kv&
zH(!S=;eZ%dM<+-s`R!HAB(E=q1y{;+HZE&h#spJ3brSQ>x=*^LnuYW#YRr(lQz54X
zqtm{4N*LsfH`MRjg_G`y<u-B(EZuk9=~l$lQkixI#8Al?z<e6v`v0^V5J$;Q0_BWv
z&&|(Vv-x##$1ENPI_9_d_YNcdKyRnIFNQ=SKtT=v-DfP<s|9X8X5Smjm4>loTg+$)
zPboEqOlU!T_b#~jX|C9lLZ$Tr@8UK1@;q7ptD}fD9ww+=qPtFbTpy$a-P>GRGG|F3
zsCpjk4nM$t!j~Xlz*Z`;HD>^8ZkILpT)Bq@&uP4o+*i+j$<zT@DKRa`+7qLS1pf6d
zE9aL=PlWt9H2}ol)sOHV$4Wo9-L&=rhaE6)vhO+*{f9zqrs%Mtcdy&)jK9wva>yfK
zHa8W9VWo*Pb|58Z@MUFg>=+HGUu5{(+~$t2?=JM~;T5l#Q^Scqiz2H1Q1+3XL)_vM
z*%=-ApVOm<&7vE%_RWFS=~yyR*I=muV3vtd+;uOmj<~qwmYl<5mJ%^CX=^N9-%Yx$
z+;v7FjCwND^OI^E@S9b8y0OLnOA*8!ah|s$oy|!|(Ya4uMAT&_I-hxMoImk(eFa+T
z>He1J4Ack4i1{kWSvEb=l{PN`ZA-bO?t#J%z}7Sk!vGSE*y%6X()7hP-&yr}QaG5s
z7%b|X@{oc>di;S3FBN@-$r5Q#^IhjCzs*C4$Hd26CA)Wuz-L;DPx3>oS`)YWFhVi`
z!n_^Od~mV}^f0S1oM_|N2Etao_o(6W;KXeoFmImipt@ZS+x0)L(1Q{2d#n+!ZZOOp
z)-(X0U#o>3@JK%Ns$#3CxwuM4{E{Byt#SJ_EtTY-iDp4S&I_)HqSv;wU(E|#%j3aY
zK>2W7G4c!8+v(dG5@d7kwSS(}o$S=oL=aS?XSEE&=RNlFiFNTmU=Av9k1$4>MFmWS
z!`4qsA^3EM+BBryJe5sMh-7f*l^%rsjr<*C<I-z8Y128XhM{}=xzf}c?<hjAghJAd
zT^myTr0R><ZFS5>QX(yn5-xsRK#vQXd#NvTA~Fk9A)YqOJGP{H^Puatr#~g}9`Wqc
zODVS5-OiZkY#7a$rtHmw%Qc7NwgxgsdqLu=raeNA>1pk!{%ks0D~+q2KerX84Htp$
z+{^%AX#-Qb2$N>H#xl$CNN=*Uw+{6a){60!C$^OCdyc7-O(hds*l++_lb4aYGvMy)
zPiJU(?xmgp7EuivFv{1$2=m)Yb4}lKp{0F!<`~)45Zlj^bVs}{Geh{%)|N=YQ830<
z*E_>IE}EmR*yG^9mkkK@PgK|x2KP!_|IBFpF^H7TG$!q|3YulS<;%cdCK4n0^!P!`
zC`Vt}=-@ZY2rHIcPGOfB^K|)eBW<*9*HXTux?xPyPHcuJxrb;x140o2CrI3ipGB3W
z<p*ooC9{&mE9Xf2vil<&^fT^Mz$cXRlSD;<(gWT6F-v~0uL80}Kq2ItF_mrH<aUH%
z@-HF^CDi&Z^BYraa|SP;QfA|Qn_5!TcJ5N@pa2>r+j&lH2#W2Cmxz0|s;bV(y^@!V
zw4H$4m`F+vnj8#{S?qj7o*n&J)1Nb-NvzZfe%D4j2c2ztI4}@Tv>_5wn)=t@f5H7*
zE#lcFcgL4G8Ll`IMEa6@pzX){AF#Q2v&+Tuh$@{F1I*Mn4zpAw2~9Y9qO9WLJNa2b
z8$biKgAFeqk4Z~qU+zEIIhYDYbNZxjIsd$p(epo;0y5X3cQT!3hf=CfHE|<bf4k}>
zw)0Nm!CJ7U4sR{2@LDWZWsQ@`4BoL{KkoBD+9jnoBqHg@X1g&5yzWc}xz#=<`BNQd
z0V#LtC`Wb>iHMn-N6QX3SSdhy=b6bR!JQ`G9^ZJN#x~uO0Eu#d4s!;uKxmOpBcAuF
z83<}czu9&7_Y~Ii&aX^>%o+8&<kH6plbGo$2n}S$GpAcms#aScz1?f@<B=WQ5O-wq
zbMoJtagt3_A6v8n;*FG9+Q=vxEq+Kd#B%Ti&QXYBKv|?(7V7V)0fUv$5qAIAVYwPI
zxmI2sGhr7Vk8QqFNV68aO71X;E*Dp~yw7G01XKp+Y7MXQZ^p`tAQjNtfr~F-7Nfh`
zt3bZXzMbpV;GQ>sDt}xv0UaIj0rmnjY;wO}*00R7R6KP=5Za(S;OyqV$Ko?~3Z_P9
zC4WD>Gun_;UWV4=B{v_yRjc8>W2!WY<zw$L&r3C+Yr!q1>7DfJj?!>E!(iZTC>B_o
z#A)?ml{(fLA!csY?J|Uh)`}_kDTo|KTQnm44SiY1+-FwbDO){V<aV+2=v!Hxy5l=}
z_QdAQ(F#_bn$(L@(++AMAnKORsIBd(bmJyDVabs#_V{Ql&<9jt7{PqFk$T6~JuyEV
zE7UN|WgC3Qxck-W1y)ppEXOFq!&yypl$(D^Zv@cx*0*Ib>d93-r=29ppz7(t_(FTA
zPY4C+TWN~_^MLlVUC<J~En){fyzN<yJ-<^pvQuzzk1%vcENLt;<nDL(PP(|G+=SqC
z>t1fKUurW|O1!fv;C8~^o>1>Z`Lf5kuQwbevq`q{!Td6|En`+-nu>J0q4eQ$md!4m
zHp)||-5NBbZJV7QsJpwM7MB1BoEgOH5XdL8U1yVF+Y%MOGOtK0FBJNs2Rr+#^<8tC
zso!iOs=1?#A>!sVEY)t{Nv~eq(9j5ZQaS!abs5ga{bSx_ZD*B?bcUl>%EN>X_1dpG
zvFx@-YUs;q?U^zaSggX|Y6T0GKG)X2r|0ci5hN++zpVNF%`sCQnF(xMZ}|VtYshh{
z;fz@|#_5d{VLvBJXZ^4GrLQFF9aml&R~9-kr08M7*7Kw<0R8(mZ5DM_=Lc|aXVc8V
z_G;4d3TUML7U}CfjQI<c&u_J4<?Y!5`~s<O-lXdtg_Fe}t8Es?jM5COWJcFF))WmE
zafKvsUobzg=gCKuBoX>9=^U>XkDWQ#c4c(d=3S`|od=3<lr(o8lYCaOHF?|tmx?~H
z#puMk_Lt|M%){NaHUn?E3+&EPPSRh}=*8ogl_p;iZ_gqOU8$w$=}848o67)Orqbj_
z0jOPyzm~FJ>IfTIdyoRhcT?O^v{lmcTc8!)&cpAYJ*|ln`lJlaPSVXMHgV}3_-h;&
zo&QXk8fbBK4ykFa%}YG2<oOqGy@0{0ZNhwE#Fs=M^hv`#W=Hh8FN4PRI;ZHNV@+8`
zprBzsgd|VYzB@`s3Cxzoe6?f(HV5FyWR&vJOe{U)IAbwHX7HzNPz+nM4*1IB(ixPx
z2;ov_g&V)>P-d7F%)+VBngQ<CGs)khU@KpCov5uHcL#iyG_)<@I_G!$(U}#x4sXGD
zqS^}E{RP<-MVoYL4T=0r4y$n1er5LHKbz9Og})tJ-gI2qKUMB!H}NrZv!S5^7GNBU
zSez??q1V+pt{ni@kNIQoDP)|`{0D6(pnNn(h`1~*xP|s$+@f87>FC5welKw%{u{_R
zRCN=8<fQIw?e$bn@kIS9v+~?V(E>QAFn#kU*d#WA+(*M}9Rax7CS=p=MRO{|o0bzM
zyU%f11pa<RHh|+RW~>^~XQ;~}%?Xkc&{60#m($%0q(A#YBL!FbP2`WYC=%KcN`y~8
ztJ<RR*UsB>1`9aC#>*QdHsE<j(|6LWa^yguAExp$y<El8=nIT)k*^r9sd&fk4OKZw
zYzH;MsC-=Uq}18+<Bqs=z>MkaEmKaBm8g_@w)`(3b|GmD?Pf5qnaldMx+me^FCF`u
zR|5tlVxsjwxA#p`C2I6U&17PE@NElGC`}PO5$dsA0h>d_(<dzJ)(kplX|HU&sFBie
z8}X1@+<rte-VB~SDkks-PKiug)T#pF87%7|aLcCWx(iBx)2O9H<;FL7^1%h;b3g@@
zczaX{=@$+g344n*0fKa|cfPtt4q0oD%(L(@ROD7utIPPnNEN=0`@lt=pflik2m41G
zODP`uD^pL{mDm+j1(xQwrcPve6eDO+lp7JOI>3j_ZwQjI{r4fazWG?wDm&qZNq)}G
z^3oN76tjpLm%Yse!cR$ElB4g<+FLj8eEBz0YMBJduF8JtrQi3<OKF*J{FME&msjIu
ziO?GAfL&SgYV$2&kpVbV$Xi$eV|jd!Ya8+kxMQ{X((0OSZ@d@!wj(i2m5a+xr8e8%
ziLo!zKa>#HqN*VXQN)F-I4r|b!j-Ax0a<)cGlJe@zT*H|GdkJArjcm9ud<-tQO_nM
zE4#jb7US$spuPXX#0`4?tj1!;e41ZB3`tRO@d5^@pB~?rBf^7_2fFm#sHP(*8Nz>~
zVyu~#Jvix3FSdPU?4Ul;+LMWGYMeIQs}|UShp)>~x^Dwq041-)Wcv@_nsczen=g1Z
zx}r=MJdZ+~<tU#)Q@;aPx~6k6*_rSvW8s)X>hmNbcD93tdl8bt(6XE%%oViFTctL9
zCmoFF#OsZBo|XByzs_HQ=8-8Q78kax>b=*uV^PAJX&KOCvBGSPX?|(ylFXl}A-^4{
z2t_!9ep))(F>FL;vJ)E$7#<(EW`3hS<`(Oo5pvA?&nSUvKUdv2en_{t`du)>Se|%x
zt+n>e)fFp(cGCmwKf=9cdzwL+vYAGL*hNn8I}N)xExnkPtQlxrmedqaX|!wr9+46X
z%=&zF1#Z<O!205ckO8RQ2`t%UUOC3)J`+l&Sl`7ji<=%SoGzzjVa(+)3GfU-SzPa>
zI2lPCqdgD}>As3Y{tW$4B-g0jAcGL?D4BJ2&pa)D=Z?L65=l8)XPQrk`kk9OZKy31
z@C(Y2!@EkBDEmzV4U(}~`8b868%fJioOC0L7xL9sxZM*Af<11jQ|Zq|OvZ`F3C`1*
z4}VBx3=Uba{%cesXtuDswxt-8!i%j`rqlYQrjKH194A?+`srSj&=FMjMa7<#GS#pA
zt6gCqop~6<jsfs-mRK}@auM_&hMFmU0Orf%1Uv)9+L~U^&x%<^=fM1}$}p&8Xd+9S
z&F(p4>XQ5yGRpv#ew9j{P8X<U07ISMum657mb51^fd4@g8!w-pLJ4`TJ3^B6QI23O
zPA9^><vgexX+@StuO#qU<?M?Q+sYG7U5Skan^830_fTCUHNoL4Jna~p>{xk6=A)W`
z2w^EcRM}iL3&3qXHkS%f*uYzVI7oirVf_eUAQ<sTj(U0ejK1WS1$a20h$B?&SKn7h
z((GW`S`~RSPUw#V4Oqi**%vdptWB0=bZ;rq3WeTnFw@K1Ajx*Phx|sP98%|vG_>$f
zhU~YwH9rN-Mw<LDKw+1D8Bkv>nKFKD<=+NiifQPV%Tl)@0zs%iY~Z%$;nNc%k?+Es
z_xA07=#}s0YpwyK-1^kNj%ea;`G{LeFhji11kvI#(2?)VMrrztikCW%Mnt`?<V7fk
zR7JM`cETxt>f^dQ`_nU2vS-f?=WzN{?Ee=CR#v<v+NScE*o9H3Z(eEdQtH5i+(?56
zR~8x5V&uL6!7Jv9=tTE;7kg`L?sCY+*L?_tbVNj$P70WzGtshCO2TE4!!3W<^<NR*
zVl+<1^Gmeq53abJS%IeiLgf}0u?-KAYvXtHtj<1W!|J=0g-rtg9{H%c`y_IIk_hBL
z!tU2el`9=OE;rkc-Z@)<T5rw9#x5mI%CagdZF#XLu(A>O{1!6d9xN0yIw9Oox`aO)
z-EJ@uL>;(WGA#3D{{EDz4vbFc@GG|z*<2Z8zi@X3M$1V|+1wj9pc|SFstpF{2y`Ji
z&J){P*k37wzh=)K0&2-2HB>G6Ixu@tPg|@G$J)6+OC;((_69c&o5o{8Cp0b)87)TC
zOI~O&9<n367_5LDN&Q!~9tJu|vE&K%F42W>7MQ!Xw&d^9KVb6=D5}v0GQ26g+=Mm}
zO}aPYw#A6u;$?R~g)&|O04tKvxYrXu9qIEo85l@|(3vYGz1`PM?Ox)J`=Q@T`+gVT
zhsM*fkBtJ!mGynAy@RsZd%Fu>N+gghfKF}ZVQ_6QgsuSd>aA=Ove7T!lR<#0k+BH_
zk;?EQkkYxlpAtu?cYF649N&Pe5Q(oq$!f>{;^z}76PQdZj1wI*a^|`ziB&)c?>XA7
zKU$LSQ+NrU4i?|i$+Au(&HctPU4T6|#zOb+nS(n2FtWb0k^EmjXhozY(GW+ehF!9r
zfvX~p0-TLU5kp4NV^K;*yKYmpPb`tn>o>sF?;a(e)=K)?s%*`lUt(vbyw4HYivZgg
zLwR--TqKPOJeDf`(%Gle6QZ~8(h!%<;}mV;xUCw+2s9N+%N4ji)vTQW70A3WIPh^s
zE2bED?3TMceRPZ%1EMC=)18Qn6aQgs;`W;*^1A*MxcfrJGe;4Y7Hd)AQ8P?|g^sqj
zw#gD6<W#~a^jo09Qc6kS3Gae@(<c4TFUQ}&X8C9x>6lGdr_lR9D;U`nCFK4J4C|hN
z{Kzi%hCcy(l(4I@;tGI;oz1Be?>2+KB`5-4Z*lCyj*qj?KLPKnmvXS5@bNDZ4WWzX
z$P<nc`~CI_+okA<UX?{}>=zhHo}ZD4_;E(!F%8BcxpJ~nMR{3C-wWZn3aohNf1UZ=
zp)I8Mm)JapdA%7Jy$Wt5P6B~Uii#4_asFd=jW$NP;wcT#N>$>cq2y1OZMnU9zNG1p
z3yZaN($1^@hLkR~vI<DIkX79*Sk%1Y6t)W+%P3U4_Ehwl=UlC@W}cmE_&R!Z58g-I
z(Y`l`ax5ACbWPdu`e@@~XjL<1yKq&5Wi!FG#4tIa8$Bseg3(Y#SmDyTkkz$>9hO%Z
z`+kilz7V$BYoz0=asnpi&+TRKH(7o!m3{a9{%SbAKC*S`2U1#v^sIVC{r(G3?aY#S
z>mLaP%%_?+q9%@_$hvbNIDU#l7HQ`!=bvaQ@yy&V>Tt%j<W61fk$<$k%nAA%J^glA
z#Hj`|^Vc>dE@M0xXgt~IO0E)*0{0_0(rhkj?4k_Pft-QT@lc?Qc>b?}sf)}b@af9|
z4iZ0v`sOBBSRt7{mT^sLjxe5X;wK>>F%kFu&(NM%ODWkJp=E?{Mqow28s!emKF2|1
z<?$k>9dcmBS+qh4zaWDbalLy7JnkizwfvAp@&{`9uKV$T!{6<apIpDh(@636?+?Av
z%4<F+pL`z-ipz%Kfu8?tMz46sjQrnd!3}+gUsOny#|&4q8_DR_T`$ZSq_5E`ul$tj
znEd$S<-M!4L^O5t2q%{lu|<TL{1T8N+dKAar-J>kbA86I+;H~YtT_zWyzVnFB%W^o
z&73bZRbsWJh}f^n;sZ%QuY+TK5-_z>qS~d(A@=IwEN(v>$H{f3;kW0s-|Zg-#keqm
z_~03e=kaXDQYOdsFCI_#-z|FSKmES(=pl#;E_%C}6&9@b-;7vk^})VxUYKKDKZAJ|
zaX0Z9lkmNtXaqWf4mljxQngo?7LwD?+Q`OVFVXsK4Fy_0M*Jx;kXqi9mfU;hocB5F
z%@v^9calHbrC~*G#aG)(L{%4;<>~S3Sqb`Ox%$Qs**--O5l7L>_ST(li9$#h(<A~0
zjByDre4lW=VNGFjN8L;|8|g}dB85Sq-`4Qo1e2&4*Zcj(urgTC8c?-YhVLWLCUY$P
zc^<oXx3BYGP9EdPluO?(W_B>F*iZNj1~Qnq4Pm{GEZ$+;Ams9a-t3PbelKZc+uYKw
z3Hd~QxzImC7BEY7STE^+neK(Wnp4u~<2K0u>-M0*YoZ0AVZ$Pm6X(9`AX?nfxocz>
zQT-DlUD!vhvFrEq|G{Y4Rd9nbWuB|y<#6fZgqK4rMJgKda&i3aY-PJa@k(beZ(zEl
zr<lMkZ{%y}r#5=a@1M{nKfjM$mL4LnKNNeOzy1#?Wo{mgLV59w)4d1d<#SJ0aMFDU
zqhb^B)Ro5!1uLeu9dZ*y_^rM;MsZX9f)7H&R4@5Tiyq5z!tfQ<sJHdn@~6bQB+r^;
zIuXLvlq}{;$K8#Liz;^P^`Li7fewuWD-!mR*+0!z{n_SP<X8lsk;J{m(yOC2+<%Q{
zOT9~_bGV&aC5$0G@OK}?0!lXJ6&y*d>Rn-&wyTdd1m9He$PtrO?}V#4(6=u|l2-~$
zg3{?PSi!Bo)V%OLF@R!*_RaW#vA-sl#{dqTfLy+Sr?JFu8~nmt#L4qH^+An2pC*li
zvn^m2y>k&ukav&jO(Q2>{}@bbEY;qn<WjwDp#@<Lu23eOm%g=^36cwO>y+!oY-it(
zTH(@%0GXH05Y?$a7?wB<f^-E=e;koZDwnH3f9*)-*@4?zr6GA+&o8zLa~(^S=5BMT
z`D>$Bm1L7ebbDq)+u6+{uSn{btI);J*$G#mLH41d+3L`{+gMphLl$rDx{4fDXQ6w8
zcvuZKi39gCw?dj*EF%UfC03{L<awX^q>@f%K4#?i^==OPPx_Wota$ySFqg4ZPt?^R
zMvypLnIs4@^<xPw=SG?{7d>Q<fKkskS(1yv;ss2fNDWJ5&Y<yRI5{230t_#YfBf{I
z;I8qEIEIb>VSvZyi1FefNPwhv#FLM4r7HICF8UPFOOVM9u%}E$;rT#)v7pX|A_Ly7
zy-PT%xZrcj_}^IYpUYnR8W%vDCqG{_!4=5KwRD|2Qm`ZlDDBuGn4cPKQNR>a76<5p
zS-*odQ09*fiSC{RhSffUL{h$J{!tZwpMwk0u!mI<+ZI+U_>=K+AjdWLUTHF6dl9&A
zbug)@lAb(<AS&TthzaKowjgS&8}Fcw@x;v!V7Z-&`{VCth+a7vWTxidMDVh9=`jq5
zY}FdDGWTaxfXU0PX*&AMLXyZ?20DWwIakZXyaN6)l!Qmnh-|iRUB#YyOb0`NMJE<F
z#8<<YbMEi&u&uYzt*L6hMxeXt@Ica(9Ny~GaeEqP!x<>(`M3Kf2NQ4yiqueuHO174
zgykZdfSfD6uS?#K#WY5ETfcfUvTRM54~YHlnx93R%;bxKA=aaqYAOAi<IAq#$fY`=
zV*jgxHAx%Za`XO9930NwwRa3hsX&kYfnq-Gi$jnAPS6G+=Uto0*Bh>gbOuR~K<<o+
zRD-`sK#!mH>G!2&A_N_j;TGxH>X-NMM<1+@oak1LJycj%+$<1kHs;h9mp!|FMmB7n
z#WZ-Tz{m?2{+FEl2v4Sf5Y3luuon@(&H|H%5Jc#X>h-d^<UHL6WNo3@cf-p#1@Dgp
z%`VW{HsItXKM;U@-(|=cXUj@onm)t0P~L1)Hw_Zjp>EfT-02$=n+Zd%GagS`1oNo9
zm6f)v?=pGywkrEOXoIgVI=^kZbq&oTg84A;KG46HnLM#qtA<;h&M<bg6kISk#T>!9
zu1?C931wTtH<c9DdA5YO7b3L4Ew{X`j>`n2Ep2dI@U1K(#RMU6Gw>7r_r?R>V$>4e
zbr`4+jo140Gc5<M<37&jI#ZUi?7N$|9p9?mzAerowt&Lr1SH-orZ`F~We5RdA48Kq
zdZqPg<LW5%^rL#Cgv;SlsF32liYPoMDVCZ<@H05FU>TnCDIGbCeWChh#a!|$`yw{G
zcZ@<lE@%C?A}70_1{hy0{5nr;?ObSt|6$_pb3?EF)U6A2HC8@ma-u0S3}yUo`O6Q8
zo7&#;p}1aIqH?Kt@|BvU*;5_86W^x@KL-2}HTXrmG8%pN<(+9+o5#kyKumA+^QSCR
z(x?CdKKGGmGim~6J)@@GQ?8Y5n-q?-s&k7>UMw9;l3^MjIde9ojJJ{?`EOUE@6R=b
zTz}yr#xtq30{Q)KJTkwWQ6X|-Emba$HWK@48OlIz--k@_XHnW#rkv!NO&#Yyc@PEe
zJ{Y-EGh5}W-)+Mr`Pt!Jd+j9^penLDZSqzxSA^B9^reQNCK{hBZ*Q>D5Rt=`OV0;|
zzW?B*7C4H2A#2e7SIK<1X9y&$A(!;+vthP_m-B9pc&(@EKfgRhSA}^!m-AK4Ed8DA
z$M3J~M<=YsDR|y@lgi&OXBu9#sCpiHyA?I?beOr!3tfFBncWZZW<`V%m2ZD=6$i(2
z=dK^aLFD_88vh!9QX=GERZ4Jw)UlqJs{mqKq!8JSsyk9596tpS`?l-9rJRLN8u*Sy
zEJ6(Z3(^qht~<jl5m;>EyUJ20;=>2y|59MHVb-%4plMh{yexDGhq)>WxR{$wZ}YQ+
z+dCWWgXllUpZoM6PQpzECJN@`FKMb2kH}h{HK%yz4La8CGtgz}!eA@fxQh`QKq{VR
zGcIi65E<KE+8Vx%>n`}>dM96tunteS1{IK^4HBw~FIR_0FhZ-zq{e%sxuLGPGTx;F
z`Mt&p=G5Gd_GxOf9uhNo^ix?`Gjt{7jVq)_(n9w&?dPv)uIEMz%~4wB@RM|D9x4pR
zFiH|ep}DIs5@L<&yq#Oh@<!tCnx><OFF;}xwP>F~ybnJi6+c9aIUDDl-8Og^{Ejpi
z2yL$mji@+(*`($C0rE&OckaR?T>S(+7xIor1frwy?9SrtLdG>?^VBT4MN6(L<yVfg
z!5zMf$1wep`nV?*@mg48gQ1|h@6EOBCcLn93J!|ulu5ftCosav3<J<YGdSyiE3~Jo
zzh2E@{;7dzwK~9ha$`din;xNwC#P$^xnFa7X*s|C1G~=9FUdIo0?s}I$z*y)tc(@>
z(}MZ^wEX*KRLo{U2=0O4kyuiwHp{tL%~lws_b;3c#Br!!UO}vt+-qX#!sb-ZgmU5f
z1a~2BV|5_$E5gVAl|_nE$h9tHE4DIb&HwDBg1vtU20TH9=Smsv^~;xE-G+#_SIy|m
zVnX8>z&Qh&xCIliSFj3c4vWENJQ>eSR@M@q*nEixUPATc&28e&`7nO6AHiD)eLZQ5
z74Zx{$Sy8`mbeW#Bvq}$fX&23F`#TSt;8xMxnz5Ne}L^@g(;e#;&P?@*pwZv$dS#u
z6ijOG1pp6TeE@hk=oMt>A|i&_|1n};CHSfHg=Bb&h!`>{xIg@nW+36hDkqdhOp+iU
zJzHfi>~I2PCRUl@>2?S}H}&U6L$iuaD9Xu<s6aa!E^of$<Z>0*ntKd>2hw44tJ`P%
z!G)s(Rt%vbab`PQGy;-l{BdH}Fc#4m3-+F)x~z1q4T8w8)nObxE!){S55r-!RghWT
zR)9F@Yd_i~?z6BY`bCm|vq`6)o4~wL>){`g%MOEfPv?qabB2I<@aE>KN<6Ux#Ail}
zm30rPCiA6G1aEK4*SdoHMQa6k-PRJrK8!JCP2Vj{_1#1S(+KW7kTY}IiUy_kYw3o2
zHQJWomUR9~D^XoQ=<YbvWHG^3OV9r>V~rdZoRVG%!?OK`&mfmbMzGe>#RY_B2Gzks
zAfT7Ecl)mQpuHKD12dqP(0Sqs*DE7f8yD4Ffs)!{Fe{ga6n9h5Cbl0p3;w-Ye_x37
z{ell;iafPrA3xKPOvd7WHv{_izva2Fa@I>KVylkO@f54b4Ch07v(x+dv|fK`ao#~<
z#7Ug}d#@3RK!~WHVYd|6%K9JG-MhftB5lAm36;e~os`Z79ctg#&p2`Ozh}Yms)96(
z!9X^C)~@NDA)Jw<G`!*Yw)Fnf>oN&%j!25sRta=#*>7IcMPXiCr~R%eqaP=Z+{(Le
zZ*f;V=nzEPvcBwupgU&4PfJ_De<7=ytN;MCACOPFelU>g!v!_mVvf&G<ICH7B7?_1
zavKLcoo_UVReoRcTWeLiFwM>7MAqDIG9H|>Dq<R*@pjhXH<##>O>5SD_y<jm{-byw
zAf>z?Ev|%xs^lCPf!ouk)S+VaEr^(ne72Z-1-H9T1{aG_Pc}<eD5%e5nW2Tj`g=ys
z0#T|aHCwjgb|^1z`?2Y>rD&C~BnRpdB8^t^t~`MvPI--Edvir)Z3@cdLHL<%`!t6;
zZ4&)a3J%i9F|e3(SKlZirGEIZZF+ZG<9{&oDlK2vfg*?e5nAMMu^Gp*KSgG@jUy9f
z+}>57JD_wVn!YU_3gjtXEBN)o*2gTAL#D;|9%FaAnbovQs_llq6}q+E*g*J$+O`(3
z?YZg~w_KEk%?%4|b_AQYEQ8&UhruSM9J6Hg#W!8j`3Dw2Gy(ZUY4+XV@{hk|A3iHm
z2~x+aCa-?*5?HF?&C%&OFh-W5_IwkJpVsmCo5)F~w`v&(=&snZ21i5=hjBk=D!xEf
zgUoXA#l`Qr{=9GY_eFmfpttafTw8v7jU%Mg53{S_59&nbmxryS^AxVsM1_U73Gt$S
z>XvF_wzCOCaN+N!5RM-|VYuu1_4cG%7_s3EV&*fL*YuaIa}VK{>4kZ|%PcwSVhsFX
z4fne7Y@{2NtMWzq3FwxZTC=*=L$=S+o$;&pRwo3j-vStUIy2-GmC|fTYt+`LLGvm+
z)A9gmesKu%_d2BFFYj|{9=;rMR4W{*BpV|&DIApDhV7y};l(Bu8}qfQ0AG(PJVa7&
z8d_z>`!+uQNivLT<Rkz%%Q*v8-PsImJvwh4eZ>M3xVhw40)g{KA)O;jOcp(LH<Bgx
zRc!*=T*+h>x}+RFv`*bMe`F4wi2aq+afrPx1YO_Jlc##)r1u)$gyT$krt23_X>RA8
zj<|k1y`G>paI;o_{b&n!n?H^pEqw4sKNU;a%~ENRd&Rr1@oP9fX(=X)gmY-zM<AYM
zQ<L@H3&y-d`iwpl@G0UXBk_;|SCS?l4tWQ&XT(lVq!e~MCR(Sji>M~-H2P(*IBKSr
zw|?FcCpsQU`KKuG>o1<9#eB8VV&rn!)zk*4crM%t3^dX8!e6IieJnIMPU&roCt4SP
zn*lTXgMmaL)QT!EtuUBUFsQ+pHUt+Y>$~`GlB;rh<z*sZViQSByfLy<HUNH>O+`ga
zmAn4||2Q3?)!_Q7wSTlm(xZXYXQ{2DzDSYRVDTrYT#kMVsq~#)j8Z-Zma{zowW5x`
z-moFoujpl_NoPsa8IRzc1$6~ie)d|(38cKvMTKOelh0@yZ@uFd^G;3TAogpj^zY({
z=$ftNt)6#e>l@~HBUbEfY@2gv*c}Zb`Pz4VH(&++`6>R;E09m4M<c;!IpI?(a>~o7
z{6gCGv}azCW8ECX?=etx?f8N`TSG*^bZ@utL2(Vk$8C#6#O!aVsnQ0v;b4f&yP6pe
zH5^2jHYr861h>(D%yyAY{w_>DgYh1IwJja)CJ;VloMbmzoCFJ%#2J=%@uO}9^$%h8
zd+ezck&kHuZAd4&t;Noe{D~jSBPjtWQ$9pp9DDWWVf>%n4y+{=$$r%dW0^HH-ncU&
zi_HF1@@on^L~}zO%~WG*=4_q6U3gs*N#Wdq=Bh&@Dvshl<N;tYL+7sSQw*m+@vv&b
z8WOA5g`4WRTfLG?Y%OcRzrb*|_I$&5XTmrQLxVb#I8%37#Ey`dsaVCE$6&TtnNplk
z82>&H!$3Y;)c#?+$rvU&^3!qNb(=+}*DKh-X|RAI7AP_2K4EO|7S~cHzSar=11LL1
z)Yj7pu2;Wae`OzUeIMjMAH!|G1qBD%xaE~tf0?cE5<Isgl2pdfV(yBx&(%qmHBVKQ
zAD%m@G3Q?CU*{}>!bpU!=Q_Zf9v%W>ymMeAnbZ=QyVFB%I2TZ(t0J+N-b>G&=fu#l
znM>rXeOr}JP|vx))*_aZJB0~{=}i<U(iYdtqpAL&yr!}4v-Ohe@cNJD>D}ms5{pV>
zsXz|Aa=5iXMr?-)b%%-`lm6BvFsB8`F+qO{PGw7CFoy~4EvmZHM5|dRz=WC|eg(y?
zHR#oDW<AQs>(Y>l|9M@>5{A5SS;We1ZSY{AWrJD59kLazJZY;x@EQ3Igu0bpRx74o
zCfHUi1+0mbGVQ1=920LL=^m0{>r8ZA?xkyj>Li9g)k7l}lpN2`XAiIYv=T+F@fWBA
zpeha%@#-8=HQb|Dz9f9~H%*6_BAZ)lxBkwHSlaFs!a)0R-wMIzri@m%-xq*EpY#E!
zVAj%R=e94q){2Ski+(STrPr)%{2bKO2+<7G;pYomqg%7SE8~RGSbqi>w)~RH{p$*-
z2jZWwh!x$|@Vlbjrztm#tyhr>=Y!3d)};YJGU{BJucUvmt<)|WcIkbMRUxpYSC9U`
zF|7W^>Y&UlT=Fib_jFIgyKiFIt=S&(^Yn#x-8;HSYVzKhYv<eQZWC;`ED+eYD>?U^
z%`-{$<f}Uik1Z;&RCw~WWoay@qP*kz8Q;5`Whz$g-TF)?au={!9J5O&_nS)6zY0UQ
z2P>MQ!rGR2_xygkG*zqkX~>Ivx0IQd?LFx<<5<fru7#P_?-VTc<OM$6yAyUi(z=#4
zfGzUGOJJ)1yv*p>+`FeEuK>rYnw~Y3UKd`o+i2$wCg2e;J6XR^nN=wN_Gc)l$!Hq8
zVRj6sfQ0XX@MVHGdtG=LFO|(Yl&f)N!|_F_Dmi&K4W8^jZp-iGbE~@I9>e1qOXS|(
zoO|olFRoWJ11^9Dd&<>LnzdQ1E?8n5C2Ty6^I!d$2m8M#E&D6u5OFHwW3si5#@9nB
z*>8Wkr7t`lc;>=kjkFt6W39`6=x%A=C4ST=dtdO~z^^??TdfxTn%U&`#=U!O!10XT
z&V5Q&A3ogSjF`g}yc04kx=o!=aMo%O)4QN1fv|{O^Oq%GfE^6Z_7yLi@4bEQF)z2A
zfAQSN3+ESf{X228g1c5P@pALto8o``4zy{XYU<kH67@NjgYE7#^%Dh-QEeZ+H@fR~
zEEU##qIp}ew>bU0tU}B;e!Id_P;>a!y_A{%&WP)SI=$wi)t3W>byBkPx7^?G^LqWC
z-}nFjJO4EQ|DV}^zXMtSfA#;Z|9|^GLx49U69a<?0|NsG!(_MaaK8g5rtvW{Fvu}8
zFbD%hL4c7#0ZcP60MD1pD@n~Oi4UnPNG*=ltH{kkH{jySy(~+B%Ef?&a6`303<w5_
z#22OJ6zl7ONOX<vc3Zm705vEB-6Msf(F7<`k)sdMTb^H(os^%S4YC`nIWKaa$bJRt
zkAwOPWE8{!K0rSK%_+`EEJ{rYE-A{)OE1Reg=391^BaJUT?ACmi{g@Y7KjNYl{u+c
zwHD4_`Qs;0%PF9CeiW^4Y+$X(&PNY{k7^G&?==8@zmADP0!4EJJ6v;dMrvwFG0+o0
X$`ImuR-n-g4D3L72x#q!4iFCj6<7Gi

literal 0
HcmV?d00001

diff --git a/sandpit/d1functions.R b/sandpit/d1functions.R
deleted file mode 100644
index 5defc42..0000000
--- a/sandpit/d1functions.R
+++ /dev/null
@@ -1,60 +0,0 @@
-#' Derivative of WTD log-likelihood wrt omega (=logitp)
-
-d1domega <- function(dfun, x, logitp, param1, param2 = 0, delta = 1) {
-  prob <- exp(logitp) / (1 + exp(logitp))
-
-
-  if (dfun == dlnorm) {
-    mu <- param1
-    sigma <- exp(param2)
-    mean <- exp(mu + exp(2 * lnsigma)/2)
-
-    gfun <- pnorm(log(x), mu, sigma, lower.tail = FALSE) / mean
-  }
-
-  if (dfun == dweib) {
-    beta <- exp(param1)
-    alpha <- exp(param2)
-
-    gfun <- exp(-((x*exp(lnbeta))^exp(lnalpha)) -
-                  lgamma(1+1/exp(lnalpha))) * exp(lnbeta)
-  }
-
-  if (dfun == dexp) {
-    beta <- exp(param1)
-    gfun <- exp(-((x*exp(lnbeta))) ) * exp(lnbeta)
-  }
-
-
-  d1 <- prob * (delta * gfun - 1) / (delta * gfun * exp(logitp) + 1)
-
-  return(d1lomega = d1)
-}
-
-#' Derivative of WTD log-likelihood wrt param1 (and param2)
-
-d1dtheta <- function(dfun, logitp, param1, param2 = 0, delta = 1) {
-  prob <- exp(logitp) / (1 + exp(logitp))
-
-  if (dfun == dlnorm) {
-    mu <- param1
-    # s is also known as lnsigma
-    s <- param2
-    sigma <- exp(param2)
-    mean <- exp(mu + exp(2 * lnsigma)/2)
-
-    gfun <- pnorm(log(x), mu, sigma, lower.tail = FALSE) / mean
-
-
-    y <- exp(-s) * (log(x) - mu)
-
-    d1FRDdmu <-(exp(-s) * dnorm(-y) - pnorm(-y)) / mean
-    d1FRDds <- (exp(-s) * y * dnorm(-y) - exp(2 * s) * pnorm(-y)) / mean
-
-    d1ldmu <- d1FRDmu / (gfun + exp(-logitp) + 1)
-    d1lds < d1FRDs / (gfun + exp(-logitp) + 1)
-    return(d1lmu = d1ldmu, d1lnsigma = d1lds)
-  }
-
-
-}
diff --git a/sandpit/try2.R b/sandpit/examples.R
similarity index 100%
rename from sandpit/try2.R
rename to sandpit/examples.R
diff --git a/sandpit/good_practises.R b/sandpit/good_practises.R
deleted file mode 100644
index 3211826..0000000
--- a/sandpit/good_practises.R
+++ /dev/null
@@ -1,31 +0,0 @@
-#######################################################
-#                                                     #
-# Implementing good practises for creating a package  #
-#                                                     #
-#######################################################
-
-# why we moved all .R files out of the R folder?
-
-# open a new .R file within the R folder, where the new function will be stored
-use_r("dlnorm")
-# makes functions available, although it does not exist in the global environment
-load_all()
-# check if function exists in the Global Environment
-exists("dlnorm", where = globalenv(), inherits = FALSE)
-
-# If you see TRUE instead of FALSE, that indicates you’re still using a script-oriented workflow and sourcing your functions. Here’s how to get back on track:
-
-# 1) Clean out the global environment and restart R.
-# 2) Re-attach devtools with library(devtools) and re-load regexcite with load_all().
-# 3) Redefine the test input x and call strsplit1(x, split = ",") again. This should work!
-# 4) Run exists("strsplit1", where = globalenv(), inherits = FALSE) again and you should see FALSE
-
-# gold standard for checking that an R package is in full working order
-check()
-
-# converting dlnorm()’s special comment into man/dlnorm.Rd
-# updates the NAMESPACE file, based on @export tags found in roxygen comments
-document()
-
-# preview of the help file
-help(dlnorm)
diff --git a/sandpit/mlwtdttt_lnorm.R b/sandpit/mlwtdttt_lnorm.R
deleted file mode 100644
index f5e0d3d..0000000
--- a/sandpit/mlwtdttt_lnorm.R
+++ /dev/null
@@ -1,6 +0,0 @@
-mlwtdttt_lnorm <- function(logitp, mu, lnsigma) {
-  
-  a <- -sum(log(plogis(logitp) * pnorm(-(log(obstime)-mu)/exp(lnsigma))/exp(mu + exp(2*lnsigma)/2)+(plogis(-logitp))/delta))
-  return(a)
-  
-}
diff --git a/sandpit/sand_test_HS.R b/sandpit/sand_test_HS.R
deleted file mode 100644
index 7d87217..0000000
--- a/sandpit/sand_test_HS.R
+++ /dev/null
@@ -1,94 +0,0 @@
-# How to compute derivatives vectorized using the numericDeriv() function
-
-
-library(bbmle)
-library(Rwtdttt)
-library(haven)
-library(data.table)
-library(readxl)
-
-# load data
-# test data for sandwich from Henrik
-
-tmp <- haven::read_dta(system.file("extdata", "score_ex.dta", package="Rwtdttt"))
-
-tmp$packcat <- factor(tmp$packsize)
-
-cl <- tmp ## [order(tmp$pid),]
-
-fit <- wtdttt(data = cl,
-               rxtime ~ dlnorm(logitp, mu, lnsigma),
-               parameters = list(logitp ~ packcat, mu ~ packcat, lnsigma ~ packcat),
-               #id = "pid", # do this to stop wtdttt() from filtering the data
-               start = 0, end = 1
-)
-parm_form <- unlist(strsplit(gsub(" ", "", unlist(strsplit(fit@formula, ":", fixed=T))[2]), ",", fixed=T))
-
-#####
-# based on code from bbmle::calc_mle2_function
-# FIXME: painfully slow
-# FIXME: redo without special case code for each density type
-# FIXME: redo to handle non-standard parameter naming
-
-parnames <- grep("delta", names(fit@call$start), value=T, fixed=T, invert=T)
-
-if (typeof(parm_form)=="character" && length(parm_form)>0) {
-  ## linear model specified for some parameters
-  vars <- sapply(strsplit(parm_form, "~", fixed=T),"[",1)
-  models <-  paste0("~", sapply(strsplit(parm_form, "~", fixed=T),"[",2))
-} else {parm_form <- c(); vars <- c(); models <- c()}
-
-# FIXME for any missing components, substitute "~1"
-
-for (i in seq(along=parnames)) {
-  if (!(parnames[i] %in% vars)) {
-    vars <- c(vars, parnames[i])
-    models <- c(models, "~1")
-    parm_form <- c(parm_form, paste0(parnames[i], "~1"))
-  }
-}
-
-vpos <- list()
-for (i in seq(along=parm_form)) {
-  vname <- vars[i]      ## name of variable
-  vpos[[vname]] <- which(parnames==vname)
-}
-
-
-
-mm1 <- model.matrix(formula(models[vpos[["logitp"]]]), data=as.data.frame(fit@data))
-mm2 <- model.matrix(formula(models[vpos[["mu"]]]), data=as.data.frame(fit@data), drop=F)
-mm3 <- model.matrix(formula(models[vpos[["lnsigma"]]]), data=as.data.frame(fit@data), drop=F)
-
-mm <- cbind(mm1, mm2, mm3)
-
-logitp_est <- mm1 %*% fit@coef[1:2]
-mu_est <- mm2 %*% fit@coef[3:4]
-lnsigma_est <- mm3%*% fit@coef[5:6]
-
-
-
-myenv <- new.env()
-myenv$logitp <- logitp_est
-myenv$mu <- mu_est
-myenv$lnsigma <- lnsigma_est
-myenv$x <- as.numeric(getElement(fit@data, fit@depvar))
-
-nD <- numericDeriv(quote(dlnorm(x, logitp, mu, lnsigma, log=T)), c("logitp", "mu", "lnsigma"), myenv)
-grad_mat <- attr(nD, "gradient")
-
-# also possible to just compute one derivative at a time
-nD2 <- numericDeriv(quote(dlnorm(x, logitp, mu, lnsigma, log=T)), c("logitp"), myenv)
-grad_mat2 <- diag(attr(nD2, "gradient"))
-
-# The u-vectors in the Stata output (ml score)
-grad_vecs <- cbind(as.numeric(getElement(fit@data, fit@depvar)),
-                   diag(grad_mat[1:nrow(grad_mat), 1:nrow(grad_mat)]),
-                   diag(grad_mat[1:nrow(grad_mat),
-                                 (nrow(grad_mat) + 1):(2*nrow(grad_mat))]),
-                   diag(grad_mat[1:nrow(grad_mat),
-                                 (2*nrow(grad_mat) + 1):(3*nrow(grad_mat))]))
-grad_vecs[1:10,] # This is the same as Stata's computation!! HOORRAY!!
-
-# Each column should now be replicated and multiplied with corresponding column
-# of model matrix
diff --git a/sandpit/sandwich_estimation_code_malcolm.R b/sandpit/sandwich_estimation_code_malcolm.R
deleted file mode 100644
index c78d645..0000000
--- a/sandpit/sandwich_estimation_code_malcolm.R
+++ /dev/null
@@ -1,80 +0,0 @@
-# proof of concept sandwich estimation
-
-library(haven)
-library(bbmle)
-library(numDeriv)
-
-
-tmp <- read_dta("inst/extdata/score_ex.dta")
-
-tmp$packcat <- factor(tmp$packsize)
-
-delta <- 1
-
-dlnorm2 <- function(x, logitp, mu, lnsigma, log = TRUE) {
-
-  prob <- exp(logitp) / (1 + exp(logitp))
-  sigma <- exp(lnsigma)
-  mean <- exp(mu + exp(2 * lnsigma)/2)
-
-  if (log == FALSE) {
-
-    prob * pnorm(log(x), mu, sigma, lower.tail = FALSE) / mean + (1 - prob) / delta
-
-  } else {
-
-    log(prob * pnorm(log(x), mu, sigma, lower.tail = FALSE) / mean + (1 - prob) / delta)
-
-  }
-
-}
-
-
-fit1 <- mle2(rxtime ~ dlnorm2(logitp, mu, lnsigma),
-             parameters = list(logitp ~ packcat, mu ~ packcat, lnsigma ~ packcat),
-             start = list(logitp = 0, mu = log(delta/5), lnsigma = 0),
-             data = tmp)
-
-summary(fit1)
-
-score_mat <- t(sapply(1:nrow(tmp),
-
-                      function(i) {
-
-                        logl <- function(x) {
-
-                          mm <- model.matrix(~packcat, data=tmp)[i,]
-
-                          dlnorm2(as.numeric(tmp[i, "rxtime"]),
-
-                                  logitp=t(x[1:2]) %*% mm,
-
-                                  mu=t(x[3:4]) %*% mm,
-
-                                  lnsigma=t(x[5:6]) %*% mm,
-
-                                  log=T)
-
-                        }
-
-                        grad(logl, fit1@coef)
-
-                      }
-
-))
-
-colnames(score_mat) <- names(fit1@coef)
-
-bread <- fit1@vcov
-
-applyTapplySum <- function(X,index) apply(X, 2, function(col) tapply(col, index, sum))
-
-psi <- applyTapplySum(as.matrix(score_mat), tmp$pid) # sum by cluster
-
-meat <- crossprod(as.matrix(psi))
-
-rownames(meat) <- colnames(meat) <- colnames(psi)
-
-sand <- (NROW(unique(tmp$pid))/(NROW(unique(tmp$pid))-1)) * (bread %*% meat %*% bread)
-
-sand
diff --git a/sandpit/try.R b/sandpit/try.R
deleted file mode 100644
index 12b5149..0000000
--- a/sandpit/try.R
+++ /dev/null
@@ -1,97 +0,0 @@
-rm(list = ls())
-
-library(data.table)
-library(tidyverse)
-library(haven)
-library(bbmle)
-
-# load data
-df <- haven::read_dta("data/wtddat_dates.dta")
-df <- as.data.table(df)
-
-# create time variable
-df_2 <- create_time(data = df, event.date.colname = rx1time, start = as.Date('2014-01-01'))
-
-
-# fit waiting time distribution
-fit1 <- wtdttt(data = df_2,
-               obstime ~ dlnorm(logitp, mu, lnsigma),
-               event.date.colname = rx1time,
-               event.time.colname = obstime,
-               start = as.Date('2014-01-01'), end = as.Date('2014-12-31'),
-               init = list(mu = 4.23, lnsigma = -1.26),
-               )
-
-###### ERROR
-# if the user provides the init argument, the following error is thrown:
-# Error in minuslogl(logitp = 1.65, logitp = 1.61280168302535, mu = 4.23,  :
-# formal argument "logitp" matched by multiple actual arguments
-#
-# SOLVED removing lpinit as element of the list argument 'init'
-#
-# N.B. fit1@fullcoef includes the estimates of three parameters (logitp, mu, lnsigma) if init values are NOT provided by the user
-#      fit2@fullcoef includes the estimates of two parameters (mu, lnsigma) if init values ARE provided by the user
-#
-# commented Henrik rows in wtdttt on redefinition of density functions
-# run again external dlnorm function in dfunctions.R (with delta=365 in the definition of the function)
-
-#############
-
-summary(fit1)
-
-predict(fit1, family = "lnorm")
-predict(fit1, family = "lnorm", type = "prob", distrx=c(2,100))
-
-plot(fit1)
-
-########################### RANDOM INEDX DATE
-
-# by hand
-
-# create sampling window and define a random index date for each individual
-set.seed(84)
-
-df_r <- df %>%
-       mutate(r_index_date = sample(as.Date(as.Date("2014-01-01"):as.Date("2014-12-31")), nrow(df), replace=T),
-              obstime = as.numeric(rx1time - r_index_date))
-
-# remove dispensations not within the observation window
-df_sel <- df_r %>%
-       filter(obstime>=0)
-
-########## ERROR
-# i) if init values provided:
-# Error in minuslogl(logitp = 1.76, logitp = 1.76026080216868, mu = 4.23,  :
-# formal argument "logitp" matched by multiple actual arguments
-#
-# --> i) SOLVED: revoming lpinit from the init values provided
-#
-# ii) if not provided:
-# Error in optim(par = c(logitp = 1.76026080216868, mu = -Inf, lnsigma = NaN :
-#                          valore non finito passato da optim
-###########
-
-# calling the created function
-df_sel <- create_time_random(data = df, event.date.colname = rx1time)
-
-###########
-
-# fit waiting time distribution
-fit2 <- wtdttt(data = df_sel,
-               obstime ~ dlnorm(logitp, mu, lnsigma),
-               event.date.colname = rx1time,
-               event.time.colname = obstime,
-               start = as.Date('2014-01-01'), end = as.Date('2014-12-31'),
-               init = list(mu = 4.23, lnsigma = -1.26)
-)
-
-summary(fit2)
-
-predict(fit2, family = "lnorm")
-predict(fit2, family = "lnorm", type = "prob", distrx=c(2,100))
-
-plot(fit2)
-
-################
-
-
diff --git a/sandpit/try_sandwich.R b/sandpit/try_sandwich.R
deleted file mode 100644
index 6cc4d74..0000000
--- a/sandpit/try_sandwich.R
+++ /dev/null
@@ -1,44 +0,0 @@
-# proof of concept sandwich estimation
-
-library(bbmle)
-library(Rwtdttt)
-library(haven)
-library(data.table)
-library(readxl)
-
-# load data
-# test data for sandwich from Henrik
-
-tmp <- haven::read_dta(system.file("extdata", "score_ex.dta", package="Rwtdttt"))
-
-tmp$packcat <- factor(tmp$packsize)
-
-cl <- tmp[order(tmp$pid),]
-
-fit1 <- wtdttt(data = cl,
-               rxtime ~ dlnorm(logitp, mu, lnsigma),
-               parameters = list(logitp ~ packcat, mu ~ packcat, lnsigma ~ packcat),
-               #id = "pid", # do this to stop wtdttt() from filtering the data
-               start = 0, end = 1
-)
-
-summary(fit1)
-
-fit1@idvar <- "pid" # XXXX nasty hack to add the ids back in again for the sandwich
-
-Rwtdttt:::sand_vcov(fit1)
-
-# Equivalent Stata output
-
-# . matrix list Dcl
-#
-# symmetric Dcl[6,6]
-# logitp:     logitp:         mu:         mu:    lnsigma:    lnsigma:
-#   ipack       _cons       ipack       _cons       ipack       _cons
-# logitp:ipack   .61015454
-# logitp:_cons   -.2340662    .2340662
-# mu:ipack   .01835992  -.00759435   .04215239
-# mu:_cons  -.00759435   .00759435    -.032557     .032557
-# lnsigma:ipack    .0696015    -.045143  -.03004152   .02833014   .12762532
-# lnsigma:_cons    -.045143     .045143   .02833014  -.02833014  -.07785109   .07785109
-
diff --git a/sandpit/wtd_bbmle_try.R b/sandpit/wtd_bbmle_try.R
deleted file mode 100644
index 617857c..0000000
--- a/sandpit/wtd_bbmle_try.R
+++ /dev/null
@@ -1,62 +0,0 @@
-library(tidyverse)
-library(bbmle)
-library(haven)
-
-
-wtddat <- read_dta("data/wtddat_covar.dta")
-
-wtddat <- mutate(wtddat,
-                 wtdtimes = 1 - last_rxtime,
-                 packcat = factor(packsize))
-
-delta <- 1
-
-dwtdttt <- function(x, logitp, par1, par2, log = FALSE, family = "lnorm") {
-
-  if (family=="lnorm") {
-    prob <- exp(logitp) / (1 + exp(logitp))
-    sigma <- exp(par2)
-    mean <- exp(par1 + exp(2 * par2)/2)
-
-    if (log == FALSE) {
-      prob * pnorm(log(x), par1, sigma, lower.tail = FALSE) / mean + (1 - prob) / delta
-    } else {
-      log(prob * pnorm(log(x), par1, sigma, lower.tail = FALSE) / mean + (1 - prob) / delta)
-    }
-
-  } else if (family=="weibull") {
-    prob <- exp(logitp) / (1 + exp(logitp))
-    alpha <- exp(par1)
-    beta <- exp(par2)
-
-    if (log == FALSE) {
-      prob * exp(-((x*exp(par2))^exp(par1)) - lgamma(1+1/exp(par1))) * exp(par2) + (1-prob)/delta
-    } else {
-      log(prob * exp(-((x*exp(par2))^exp(par1)) - lgamma(1+1/exp(par1))) * exp(par2) + (1-prob)/delta)
-    }
-  } else if (family=="exponential") {
-    prob <- exp(logitp) / (1 + exp(logitp))
-    alpha <- par1^0
-    beta <- exp(par2)
-
-    if (log == FALSE) {
-      prob * exp(-((x*exp(par2))^alpha) - lgamma(1+1/alpha)) * exp(par2) + (1-prob)/delta
-    } else {
-      log(prob * exp(-((x*exp(par2))^alpha) - lgamma(1+1/alpha)) * exp(par2) + (1-prob)/delta)
-    }
-  }
-}
-
-## Testing the density function - does it make sense?
-wtddat$wtddens <- dwtdttt(wtddat$wtdtimes,
-                          logitp = 1, par1 = log(.15), par2 = -0.4, log = FALSE, family = "exponential")
-ggplot(data = wtddat) +
-  geom_point(mapping = aes(x=wtdtimes, y=wtddens))
-
-
-## Covariates can now be included (SIC!!!)
-fit1 <- mle2(wtdtimes ~ dwtdttt(logitp, par1, par2, family = "exponential"),
-             parameters = list(logitp ~ packcat, par1 ~ packcat, par2 ~ 1),
-             start = list(logitp = 0, par1 = log(delta/5), par2 = 0),
-             data = wtddat)
-summary(fit1)
diff --git a/sandpit/wtd_class_test.R b/sandpit/wtd_class_test.R
deleted file mode 100644
index 0312f69..0000000
--- a/sandpit/wtd_class_test.R
+++ /dev/null
@@ -1,45 +0,0 @@
-# Rwtdttt
-# Combine Henrik, Malcolm & Sabrina's progress so far
-# Create a "wtd" model class and use it to set up a 'predict' method
-
-# Test it
-
-library(tidyverse)
-library(bbmle)
-library(haven)
-
-x <- read_dta("data/wtddat_dates.dta")
-
-# XXXX a bunch of preprocessing that should actually be
-#      implemented in wtdttt()
-
-start <- as.Date('2014-01-01'); end<-as.Date('2014-12-31')
-x$obstime <- 0.5 + as.double(x$rx1time - start, units="days")
-
-delta <- as.double(end - start, units="days") + 1
-ntot <- nrow(x)
-nonprevend <- sum(x$obstime > (delta * 2/3))
-prp <- 1 - 3 * nonprevend / ntot
-lpinit <- qlogis(prp)
-
-muinit <- mean(log(x$obstime[x$obstime < 0.5 * delta]))
-lnsinit <- log(sd(log(x$obstime[x$obstime < 0.5 * delta])))
-
-# ###############
-
-x.w <- wtdttt(obstime ~ dlnorm(logitp, mu, lnsigma), data=x,
-	start=as.Date('2014-01-01'), end=as.Date('2014-12-31'), reverse=F,
-	init=list(logitp=lpinit, mu = muinit, lnsigma = lnsinit))
-
-summary(x.w)
-
-# Predict duration at default quantile=0.8
-predict(x.w)
-
-# Predict probability
-predict(x.w, type="prob", distrx=c(10,100))
-
-# Diagnostic plot
-plot(x.w)
-
-
diff --git a/sandpit/wtdttt.R b/sandpit/wtdttt.R
deleted file mode 100644
index 53a677f..0000000
--- a/sandpit/wtdttt.R
+++ /dev/null
@@ -1,50 +0,0 @@
-if (!require("zoo")) install.packages("zoo")
-library(zoo)
-if (!require("haven")) install.packages("haven")
-library(haven)
-if (!require("tidyverse")) install.packages("tidyverse")
-library(tidyverse)
-
-wtdttt <- function(time, start, end) {
-  
-  # continuity correction for discrete data
-  tstart <- as.numeric(start) - 0.5
-  tend <- as.numeric(end) + 0.5
-  delta <- tend - tstart
-  
-  # if(reverse==TRUE) {
-  #   obstime <- as.numeric(as.Date(tend)-as.Date(time)) }
-  #   else {
-  #   obstime <- as.numeric(as.Date(time)-as.Date(tstart))
-    # }
-  
-    time <- as.numeric(time)
-    obstime <- as.numeric(as.Date(time)-as.Date(tstart))
-  
-  # initial values
-  ntot <- length(time)
-  nonprevend <- sum(obstime > (delta*2/3))
-  prp <- 1 - 3*nonprevend/ntot
-  lpinit <- log((prp)/(1-prp)) # logit
-  
-  # df <- df %>%
-  #        mutate(logtime = case_when(
-  #          obstime<delta*0.5 ~ log(obstime),
-  #          TRUE ~ obstime))
-  
-  muinit <- mean(log(obstime[which(obstime<delta*0.5)]))
-  lnsinit <- log(sd(log(obstime[which(obstime<delta*0.5)])))
-  
-  # start: initial parameters to be used as a starting point for the optimization
-  mle2(mlwtdttt_lnorm, start = list(mu=muinit, lnsigma=lnsinit, logitp=lpinit))
-  
-}
-
-########### 
-
-data_wtd <- read_dta("G:/Il mio Drive/Dottorato unipi/R functions/Stata files/wtddat_dates.dta")
-
-m <- wtdttt(time = data_wtd$rx1time, start = as.Date('2014-01-01'), end = as.Date('2014-12-31'))
-
-summary(m)
-
diff --git a/sandpit/wtdtttpreddur.R b/sandpit/wtdtttpreddur.R
deleted file mode 100644
index 06d292c..0000000
--- a/sandpit/wtdtttpreddur.R
+++ /dev/null
@@ -1,11 +0,0 @@
-# Case with iadmean==FALSE
-#           disttype==lnorm
-
-wtdtttpreddur <- function(iadpercentile) {
-  
-   mu <- m@fullcoef[2]
-   lnsigma <- m@fullcoef[3]
-   
-   rxdur <- exp(qnorm(iadpercentile)*exp(lnsigma)+mu)
-  
-}
\ No newline at end of file
diff --git a/sandpit/wtdtttpredprob.R b/sandpit/wtdtttpredprob.R
deleted file mode 100644
index 2985f1d..0000000
--- a/sandpit/wtdtttpredprob.R
+++ /dev/null
@@ -1,17 +0,0 @@
-wtdtttpredprob <- function(distrx) {
-  
-  mu <- m@fullcoef[2]
-  lnsigma <- m@fullcoef[3]
-  
-  prttt <- pnorm(-(log(distrx)-mu)/exp(lnsigma))
-  return(prttt)
-}
-
-
-################# PROVA
-
-data_predprob <- read_dta("G:/Il mio Drive/Dottorato unipi/R functions/Stata files/lastRx_index.dta")
-
-wtdtttpredprob(distrx = data_predprob$distlast)
-
-plot(data_predprob$distlast, prttt, type = "l")

From 34ee6b02405f4ee75149917466e2493213f28033 Mon Sep 17 00:00:00 2001
From: Malcolm Gillies <74748374+mbg-unsw@users.noreply.github.com>
Date: Tue, 25 Jun 2024 17:34:18 +1000
Subject: [PATCH 5/6] Update README.md

Added Acknowledgment of Jesper Hallas
Added link to examples.R
---
 README.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/README.md b/README.md
index 4790ce9..b318bc8 100644
--- a/README.md
+++ b/README.md
@@ -35,6 +35,7 @@ fit_r <- ranwtdttt(data = df,
 
 summary(fit_r)
 ```
+Please see [examples.R](sandpit/examples.R) for more.
 
 ## More
 
@@ -57,6 +58,8 @@ The package was developed by:
 * Malcolm Gillies
 * Olga Paoletti
 
+Thank you to Jesper Hallas for kindly allowing us to reproduce the synthetic dispensing data in the `drugpakud.dta` dataset.
+
 ### License
 
 Rwtdttt is licensed under the GNU General Public License, Version 3.

From d60209e8488dea40ac2ead214805795da659fa71 Mon Sep 17 00:00:00 2001
From: Malcolm Gillies <74748374+mbg-unsw@users.noreply.github.com>
Date: Tue, 25 Jun 2024 17:39:16 +1000
Subject: [PATCH 6/6] Update README.md

Added link to TODO
---
 README.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/README.md b/README.md
index b318bc8..891ef63 100644
--- a/README.md
+++ b/README.md
@@ -43,6 +43,8 @@ Please see [examples.R](sandpit/examples.R) for more.
 
 Robust variance calculation, used by default in the randwtdttt() function, is very slow. We are working on improving this.
 
+Other planned improvements are listed in [TODO](sandpit/TODO)
+
 ### Bug reports, feature requests etc.
 
 This is a young project with some rough edges. Please submit any bug reports, feature request and comments as a Github issue.