From 80441182b01b48a5ce8acf30cd87124a8ec1e96a Mon Sep 17 00:00:00 2001 From: Nicholas Clark Date: Mon, 13 Nov 2023 10:50:46 +1000 Subject: [PATCH] ignore random imputation when generating default priors --- R/get_mvgam_priors.R | 83 ++++++++++++++++++-------------------- R/mvgam.R | 7 ++-- man/mvgam.Rd | 3 +- man/sim_mvgam.Rd | 2 +- src/mvgam.dll | Bin 1061888 -> 1061888 bytes tests/testthat/Rplots.pdf | Bin 22901 -> 22886 bytes 6 files changed, 46 insertions(+), 49 deletions(-) diff --git a/R/get_mvgam_priors.R b/R/get_mvgam_priors.R index c058e90d..c72e884e 100644 --- a/R/get_mvgam_priors.R +++ b/R/get_mvgam_priors.R @@ -433,25 +433,6 @@ get_mvgam_priors = function(formula, # Fill in missing observations in data_train so the size of the dataset is correct when # building the initial JAGS model. - replace_nas = function(var){ - if(all(is.na(var))){ - # Sampling from uniform[0.1,0.99] will allow all the gam models - # to work, even though the Poisson / Negative Binomial will issue - # warnings. This is ok as we just need to produce the linear predictor matrix - # and store the coefficient names - var <- runif(length(var), 0.1, 0.99) - } else { - # If there are some non-missing observations, - # sample from the observed values to ensure - # distributional assumptions are met without warnings - var[which(is.na(var))] <- - sample(var[which(!is.na(var))], - length(which(is.na(var))), - replace = TRUE) - } - var - } - data_train[[terms(formula(formula))[[2]]]] <- replace_nas(data_train[[terms(formula(formula))[[2]]]]) @@ -509,7 +490,7 @@ get_mvgam_priors = function(formula, if(int_included){ all_paras <- c('(Intercept)', all_paras) # Compute default intercept prior using brms - def_int <- make_default_int(response = data_train[[terms(formula(formula))[[2]]]], + def_int <- make_default_int(response = data_train$y, family = family) para_priors <- c(paste0(def_int$class, ' ~ ', def_int$prior, ';'), para_priors) @@ -1240,12 +1221,18 @@ make_default_scales = function(response, family){ #' @noRd make_default_int = function(response, family){ - int_prior <- get_prior(y ~ 1, - data = data.frame(y = response), - family = family_to_brmsfam(family)) - int_prior$prior[which(int_prior$class == 'Intercept')] - prior_string(int_prior$prior[which(int_prior$class == 'Intercept')], - class = '(Intercept)') + if(all(is.na(response))){ + out <- prior_string("student_t(3, 0, 3.5)", + class = '(Intercept)') + } else { + resp_dat <- data.frame(y = response[!is.na(response)]) + int_prior <- get_prior(y ~ 1, + data = resp_dat, + family = family_to_brmsfam(family)) + out <- prior_string(int_prior$prior[which(int_prior$class == 'Intercept')], + class = '(Intercept)') + } + return(out) } #' @noRd @@ -1266,25 +1253,33 @@ update_default_scales = function(response, call. = FALSE)) } - y <- response - link <- family$link - if(link %in% c("log", "inverse", "1/mu^2")) { - # avoid Inf in link(y) - y <- ifelse(y == 0, y + 0.1, y) - } + if(all(is.na(response))){ + out <- paste0("student_t(", + paste0(as.character(c(df, '0', '3')), + collapse = ", "), ")") + } else { + y <- response[!is.na(response)] + link <- family$link + if(link %in% c("log", "inverse", "1/mu^2")) { + # avoid Inf in link(y) + y <- ifelse(y == 0, y + 0.1, y) + } - y_link <- suppressWarnings(linkfun(y, link = link)) - scale_y <- round(mad(y_link), 1) - if(is.finite(scale_y)) { - scale <- max(scale, scale_y) - } - if(!center){ - location_y <- round(median(y_link), 1) - if (is.finite(location_y)) { - location <- location_y + y_link <- suppressWarnings(linkfun(y, link = link)) + scale_y <- round(mad(y_link), 1) + if(is.finite(scale_y)) { + scale <- max(scale, scale_y) } + if(!center){ + location_y <- round(median(y_link), 1) + if (is.finite(location_y)) { + location <- location_y + } + } + out <- paste0("student_t(", + paste0(as.character(c(df, location, scale)), + collapse = ", "), ")") } - paste0("student_t(", - paste0(as.character(c(df, location, scale)), - collapse = ", "), ")") + + return(out) } diff --git a/R/mvgam.R b/R/mvgam.R index e6562c0a..b2a9d923 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -183,7 +183,8 @@ #'\cr #'*Priors*: A \code{\link[mgcv]{jagam}} model file is generated from \code{formula} and #'modified to include any latent -#'temporal processes. Prior distributions for most important model parameters can be +#'temporal processes. Default priors for intercepts and any scale parameters are generated +#'using the same practice as \pkg{brms}. Prior distributions for most important model parameters can be #'altered by the user to inspect model #'sensitivities to given priors (see \code{\link{get_mvgam_priors}} for details). #'Note that latent trends are estimated on the link scale so choose priors @@ -730,9 +731,9 @@ mvgam = function(formula, replace_nas(data_train[[terms(formula(formula))[[2]]]]) # Compute default priors - def_priors <- adapt_brms_priors(c(make_default_scales(data_train[[terms(formula(formula))[[2]]]], + def_priors <- adapt_brms_priors(c(make_default_scales(orig_y, family), - make_default_int(data_train[[terms(formula(formula))[[2]]]], + make_default_int(orig_y, family)), formula = orig_formula, trend_formula = trend_formula, diff --git a/man/mvgam.Rd b/man/mvgam.Rd index 3d4da0e7..f6735ee7 100644 --- a/man/mvgam.Rd +++ b/man/mvgam.Rd @@ -266,7 +266,8 @@ in a Bayesian framework using Markov Chain Monte Carlo. \cr \emph{Priors}: A \code{\link[mgcv]{jagam}} model file is generated from \code{formula} and modified to include any latent -temporal processes. Prior distributions for most important model parameters can be +temporal processes. Default priors for intercepts and any scale parameters are generated +using the same practice as \pkg{brms}. Prior distributions for most important model parameters can be altered by the user to inspect model sensitivities to given priors (see \code{\link{get_mvgam_priors}} for details). Note that latent trends are estimated on the link scale so choose priors diff --git a/man/sim_mvgam.Rd b/man/sim_mvgam.Rd index 3632d758..23390fb4 100644 --- a/man/sim_mvgam.Rd +++ b/man/sim_mvgam.Rd @@ -84,7 +84,7 @@ be replicated \code{n_series} times. Defaults to small random values between \co \item{prop_missing}{\code{numeric} stating proportion of observations that are missing. Should be between \code{0} and \code{0.8}, inclusive} -\item{prop_train}{\code{numeric} stating the proportion of data to use for training. Should be between \code{0.25} and \code{0.75}} +\item{prop_train}{\code{numeric} stating the proportion of data to use for training. Should be between \code{0.2} and \code{1}} } \value{ A \code{list} object containing outputs needed for \code{\link{mvgam}}, including 'data_train' and 'data_test', diff --git a/src/mvgam.dll b/src/mvgam.dll index 2fc4de66739d5c838c4de5bfd34dd9071474c817..662ebc286498efabb17ff3339f98a8d5b37ba5ab 100644 GIT binary patch delta 75 zcmWN|y$wJ>0EXc>oPYK(#l--6Ynb2^Um_|-QC!U8oxlt#p8oNV^N-tCbI!*$njtKR fNEi`|M5MxsOyr^vrKm(L8qo&2{{H7>>E7uBR`w&M delta 75 zcmWN|yAeP@07cQUSbzG^fK&yW%N4SMH)BL(8Zl~dE7%%Doc#DO`Y`UR8Rlbi#+?fy f5=O)#5vi~u6S*ivDJoHmMzj}Q|NQf|bno;5E8`-B diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index de4dd13c455c14b809ddc5b02f4470406e09bbd1..192c4208b118db3a1a4c94c77a516d72d24e7403 100644 GIT binary patch delta 9623 zcmZu%2{@E(+r}hCrOlGH5|u33cZ%>9l7xm(_GOsKZY&RZQ6^M0vWqBtW-QreNDL-4 zvM)0v3}Z=SH~1gb_kQ31{eOoekH>R8_kCU0d7jsK&tV0_JcVI?rnXOBQBMB37Vod` zWmf&7%nAvr%;x}Wqf)6hBQnKyNh*Q$yfvkF`1F;SeJoDrxDS-QGCXZ^^XRE_JDRH# z1gSNZ6EYEZ8zrX5538*wXT3IZW@pyW{iQaOE4BT2U<3RQ-?Htu((eYZ$0KLpg}5Pc zcn)r7eV!(41AK|bF}{`J3QB>E1)?-=YYV?QN$bJ`n_l(S{VVN2=*l!8QtV#0LnmtO ze5E0`Hw^-}Z|rQ1CS`Nh?QEiA@7xF+rmb{_2NjDiZ}rm(jX7H>SEtxkw-p!YVb5p6 zRr{?jEitt%U!hIRyX{>tRs3YbyudT~;PypyP29#p5SidlVT=={cQEa#JDW-_$d?5G zu#>zUiMAp41;6BI0=9QH2Qw=&@wlDQB>XyG6XWvC;nCC9K)FwPzm=_+b!fRy0)Fc7 z3{9V7%?rQ6GR~U}$FEfd;mCN0Va9fn)UJrv%}`b%Ls(M3TQvW!eo!sC?$#)zpsybv zvb(BpXs}~ZF0~e%n$*5K(Y3Q04^VbCt0;_Bt(|#JN_{AuZc%HawN<5s!)>j`Yw=t(QfcA$oaiLN^Ek-H=@rO^Up5-22`vU8wI>KVu8fMctisfvyHs zjb*;CvLst_mu%}FpllU8j4J{k5D5x+EP|qS?*mGsNJ8CyaW^aAXFajfIYGwUieZ}M zVe6*QJMw8OslMk>;j-u+{(QNT+MyAsN?!17+vcE07s=&&L zfN3jq2?(_=p|sH5wega`gSycy+n4=yS1DBEsPHvYrJT(+`O8`%$A&1!u^+C$m;F%> z3xg#gn@<)6hwyxvIITcQkMI_;5^FNgcm5IWm#$X3qvknJ;sJR+qKSinXDeq3A(($N zRQ~;6K&8u#A+N`M0rkFvUabX#6DXn#%$HUzIg z^O*eeqvWhIfz3^Kvl=@UAG0W2XCF-E(-fLiDc~`mD#BDs982FC zNwGb(aFvP%9$&>f63>+4wL`sCVkCj#N`cB4ZAzuUq->4ty0G;^?W0OCx?q*Uf!0}p zk68>0Y1`@zBOhh)&vklHONh@h z7G*PQi~%Dw7p-w=2-f=PUZz+rMl>YVkam-90^sU7Q{f3y49%er<=C zMhx1Y@PD-Ox?b_&yL!Q(g>ORbnycjR3ME@}Z0W1kouRHTBq|@SCd|e!Ra5TZK4XPb z0hs^rPNi*Oz?p&bYu@~k!<5P(LELtSpgh!XxNA!4&iu=f>-}jnnhSZ6wX*&_pMau} z5wjC(9$h>o(Vvn7N33aI=TaBmyZ-nPR||97zx~0~I;t|hw)ki0RLDq~*oS`gg&NnM zR6_ci`DZ+7ROtyj;_V?eH!YD!jtczCkHB^QbVf_VJ>C)tf!?5jP1$>9>3Uk@hSL$- z9})eZLZ{WoZ*|;|6Vo#Y(LiKQp|^D551yNQ@JdhE+TyR9SgL>P+jrk)$CfTLoQ#seDD2M>>I8 zQ?hwZo7?hMk@-tZNnS7P#v*peOp$b>vD#cunzQq))8;D$vt+0Ii^>Kic8$G7Wicxv zLtiXcwj0y}5$I9nAok90@=v?QBvO~#8=+)mTO%|A81mS}J&+ViYaa5~s5LWiKgToD zp_(x-rZ*j$tVE|4k<@gj#mVP`5&-d1)>>i?_S!^x4)&bRR0nT5T#o~Z75lod6;O+m zfxcEMvfH1VOD%#rX^`)j>fTkWFD?-#N9#{VeUz)Kz^w)e!At1TYfIFk>ZzL~=qJ6n z&hGh9Qy*^ZBY)@ee7n4XQK9r^+)}8z>Y~#sjRaNX!m!47Vk9JWbtG2+f-31YmU%!f zCtob5_zhKJbMa+1cALwG|AOzthL{sGe2P@{1P?SN$f$B|`nSl1%%hDi zlZ_U%T6C?BE~EZDUEinA0hvLC`7%g) z+#X$TJjT}n%~V?G)LG0E{Z6xzjSK08Funj&usdyfp~D!}fl(aa%W6z1+AqQUM0Bi% zu5%!{D6;P~<3$cu$ScvosCmWNhdVUAjjy?N`UVLa4=Hu%;K6dlciN3=Bz0nHC5F7J zTD|(v*;of@jYY+^JP-;gLofZ@bJu=SiwXx#0Ux2Vty71i0xzv`|4eZ=CPCek(zCHt zPWTv^DsE$SvcL4rP<5LudTyTG+bq5lS(NJT+=eTve694-l7;-SzPLy8u_E}^t@%4@ z?&XBLCh2ZYi^ zgFZ2W>)+~w*cn8N4im?qVZ#Y?fUSehuSFEE1SNXJ&9=;&WuqJ+<<1XZDz~weeT4-Wc8gxjH>-SKRan;JXR-zn~*X6+c;hyiAtP!S>$Zo zn};(mgpGX(A_}kyk2YrV*zmq^!Pt}_ibB)_h&Zd8aKDg$r zg`7ax#ex;9AUgr+CB;%+*RjWF$nR&}zN@IBszzr#@)detkn`qtUQ#NRBcJDBdkR+l zsx96Ap1EJ!Hmp-l@9m425J+DcQ-z}R*wm`C6!QIfFRIe=Zk$&!OUp=Cmr*q?&2{ZK zE?Jsf6vAaCOgWsicw_h9BGzUT`&6 z6mFHAiu73>9Vc7Xq0-;Ts)7!%g5E)$Awkg<->Jd$?8Y1{vNk0f+xJ36l_-x;o+K4~ z3DrmeUal7DHm5$&LCPog`>dRlRKr5zJ5l@tIasQkdraqT)X2;ms|44c6rps-ft1a! zmJZ*k62n==Z{n{_kY3i0#1WgWk9QuR+eP}g7Sc9Ebt zjx|rdB+yejwdGZE8$DD{RZPhUhVK z*mY#GhPRfB{EfjPr{VdVTlr>v?IZA?I~L@4*~4iX0X}HYp!Q0VZAtp1LZT0*c(gLD zqmX~FEt4v5^-~F{7$^#yTg7_$a4WgUZzntOD5DU-`K>k#+onz%c?qYk>b{!d4mbzy zVnNxbm$^OOxx3o`QtGUNRheb#POyg>hYuj#aVZ@}8$_A;3#D`9T$*^2vHkfZ+%HJa zc}ciWYvN^HUNA<`IY(&G;Z-mu)>2JSfmm=>K;~84up1=GQj*~u(#ymTO%LFqV>NAUm!0emjUxw>jm&j! zT^Fn8s5)w7ljjSDhipv7dn9Wo%~%!ul-Vi@4$E`feO1ge2@+ocgk=kwXzg;+0v`c+d zifR4c2i8t@>R853Ob|JML)3{gCVgtQH)&f#DsOJY)O9<`-gBsd^*!>+|WgH8u*gd1^^yJKx7Vwe4_&a_B9~r-KP;2Q#=(k z6?ROv&HbEry1pH7)21UOneH#o71G+Mt_`-xtGDMe-qbOUa(x1neo~{{y`|o-CwjBg ze=|`u)cA=|8CgM=o?bQ)K5@9@@OPn_Q2M0Okn}>dLOM4p)cEj~kQ=ZX_9KdS&XboO zh26C3l1@XR{E|e%ZI&+ZbIU8LjGFZZDX_#3z50HeKY>`2vZ>^Cegv!Vt}7}s$594v zGv@h8Abt_Lno<4{z(t_2dP#C0k)N>Da<8sW;{CF38U8+~Ex~MSUR_b2B zw{!VAyDp6#$)V07=%ORX^SQdXu7pIqO%j@61>E1ULvNSakk1~o?V~vI+ctRylN8VDAp83`>@yAUc%bmEz+ZCr10mN^fbGZH*9+V=fTP*4z zy(NGOiNL&7QQYW3VW99#jyRiR0pinkl7`|OC9>sg2QwcE?PH(F*RPEzwKe+uOw=)C z_$!RLve|K^kU#GZ2x zjlqQ#f&-=*p%1if^<7+!eOsSPjXMk-vAageDuTxuX|80$`cl208(s2y>;Ku)9?W8#f^q$Y)Rooh0sWX0 zYm~)mR_AVA8}>A0|7}%IgR2kJ17!{z{H}!JP85yA&AV1+9p!|k?boMd&AG;T1osmi zguzX}h|M2W=y3&u*Pb+&H1q@H;#24C>Re7shD^aliedjm# zI~oJtu%G!nj#zRUvAk55tHN=KWVP(}VF0Ro^{`l;VAV+0W%9KpAx3_#T`S+OS~u3` z65YaVNmVAR<+1a@ZlF9WO=DCU%TyU2cw@ag^n>OV4J7{C5&xEZhc-_!V%NtpZ;+px zt@y~rM04tFNdGHT{6D1v4&AKkX5rh=sDc?ab14?SzJ6OJp9>;sj8>)J>zU434rMEwD{Yd|hI3|C)3a zgjJxxR@8RqcBTRJ?H8qtRqp%R*5=(MLDwjHP>NzD{ddU>!fX5Jylv?eZ-cK-ZcX zWX&(t;cuBDm)t7uWip)FG#!7_qT{ONvu2af<-_e6Qv)sj^_FQqYhRV7Za)Y$>i$p( z(#or;{L<+f%aaBs!v3XMFm>3;Y1-jWhyH9ljjsmQ*a4_&-Y=Q9SNm1Aj%t_gk%#o* zYTLR-K1c}`*eDtp*-`P1PICGqqLv7z_lK-09`4&YS%@M;1 zO3#y;Q>_;YpV#N!B`KMB(^OvZMdwJN`dh=EjYh{KH63CB%pRNWHvTn8ym5{2RHF2G z$ew<*@9CDuN^v`%YuxZXo;f& zB&7>MDgUr|Oaxkt-~+4lW2HFbZhwAW{8!6o5(MFP++fOXBJ=oIMDYF{-eA%s3*D=ylc6kGXky9G6&udx6bszyrL%@J2~eo96P=J@#3-?K2!1*4K#3mq$1* zJU?hSs6O;s@sQ?h%QMm5ct05Pl(<7UB@S1(Si+kcJgd9M|J(D{LrWi$Kxf>1;o)X5AkqN>N27eS3%sKqU44O^?IHV%~cEsxz=x=%`J@{sPPGbrZGC>><3Mu)QiEepQz|Ftxh1Y0L1Zz zPZph#l3$SnKb8dF;Tn7%eGGs|5jUDb)nhO>e@!N%(hTj&A-G=IBA;gWp7tcq>A*W3 z$K)o?&M3J%z)7aPF=q|vkE47tzBA9__KF-M72tuVYFabQ6ta%oP^ewy+w-1_%WY|iqK8a`Qd z9314he+}}^NxkBOKAPaVhvv%oeG9&+S1ylHxD5dx-SSCYr)M^{c+ZAsjypn)0}wg!?jkw zNSGca3nnGNFZr6sArm=#?}v`6;Qc2LmP9Ro@H1^=P*mhFm2vu1qN{U1Fwj}W@0s*l zwZ_(tX>&8nZwlE=H~*8eov$}@+O$TMyaF4mwQHStLrd}uXr4A$8#vsp`u*+s-15S;46x`&>gSz%0kDr0lMt}T=0cn&A zleADm{VjNZLOT}FYL|T$^BtMAbMDKI@HdnT%EL~!xIg@hB~;+WmA(r!Fk6nNqmF6~ zSC$TQUI-d$@zYeuh4l?j+p$3er6x713^gWZJL1??e}pb|bbH|kYXwc8QRWXsU~V49 z7eC?NB3(-IKt=<{LzKWD_x3Ed^(EDNMj7%&`iCcKizHiKV(#PP{2)7rPMe&mhAo;; zsR%otIAgHKCs@9scRLHa(EsvksB~snST@<4qf>43JsqO9ipZ)q;e`)uM9}A9A04lA z;lrPDGo<2(CH!*EpjG>-sJH?AQJof_#DiTM_06FdU&l`dyz1UQuXuehplaelqLwtF z>7B$y8ozN6$8G5(E^pH4XHY+IzxBhZ-dAcPVrX$I+_<59`SWFF!a!KJ-Gk|4_x;1{ z)dv&&^WL_9Kl?B{?ZheNsHprkMXBh_i^>gUry69QKH1;DAMoO|hA8)Gu@GvQ3Exja z6x=^r$svmNBZ`W~Jc+$Vai5su%9tDDyjJDE114Z3$Lktl_CwvzWe~D4h z4pDAFJ2@{kV9g+TY#%PgXtDzNW@`ARF>7=WL^%r0!jwS#*0i%vfCg@Nc>EXVIT40% zgUkjQ57!G5%iP(P0Q(iMu>mfoT;w=#^F{EsN2`c|ID{~OVV2J-`hIzvc`rB*bPFjd zb2^IHITFmB$N5kjfO*iOXRa54oNDP930F~YB)G; z&51blXAo*Ol&NU45WaNFe-;lg@bJ5j5*r&(Zjh-6{HeyS2M%s#K;RPk-TXs`r3Y7+ ztyy4`e**p9o6+fbrns`^)BKyMfgPcSBFu!SXbo1N*mvD@O{et*k_mk6vWt{MA^Pz_ z*XfjwcE#|$h_7HEpxZEKo^ zUD-qw8XHgtQS`^nh1;*Y%`gjK(fr8)Ww!NZ_g^9NAOdtHjqvw}G;zzE&>UF=Y#nNl zd8XFn)w!_&pRrea?im2>a}|Uz_edSAzw8`JD){@8E}C#dnOB%SpGX@1nAW|$zAAn$ zrPlfm#;k>g+2f=%e~$sMGqed2u=MM%)mh=LE|%-b?fV2Sfvr=Pn$Ed@2znbV40qZ> zWKjqKWW$l*EV8svKXoQ-0G^)QDhD=@B?8GKUg4F4ULI4&FW z&mwAvTAOhMX0L~TjtQDP!J_M6UncuR{>trS*AGe`gc{0dOaP&$RFxw27&KEXab>I&^+<_O z7b#{!z%8zSV`%%Jaujz2sh21l!nkD}R_hdo{I^L4=>juhdQ4O_sEvUn4RrCO#RBMpw7*2;b*e8j;1oxZFe(wmBLxfGa1rTZcH*n`c z;CeKa*&T2mnCq4nhm5eYPaN)R`Hof2+vC{Nt?SHDQtoT5WdAjTH@vBgZ(F>-&d3I- z`vg#SY_e5=oOc&wfhBTD3{^f^b7O^O<&w2DcsdRNw|oNe zGz-YE3vj6zv}1h@A(Sw28u!dK>T3vN^fbclo3ph&&bVi#|B@0=4e)C3*f7vPa&YZw zzO`rJ{3?RM8*abNMCiGh2KUV&{-Wiht~#B{22F6i`TjpAwT|kMD@jgcud2pQ9NrF( zC%TomsN9d-DbFPRCU_dku5PXT!HwELL14@A1{rsfkA3)1 zM`JWd*dqu-`q61)w)1gZJtve~iuUif!z48xIf8h~Bdu9UPj!}j)`1irySoc+egAB= z$FNIDf!1YE@^AB&X8m!pU|~qvRepe!8gthV+*^wtnP0Hk$N~if?1N1Cp^o5VKKu(~ z1dx445U`>q{ZU@b8zJyF@^!9B&Bk1Ml#VVbyFa`*BA^qp$`l7heH<*@7Rwejm}eGn z2XA`t2KU6CalZq?&=In61G@d3dI>@Zu$kOb*V~FTV7{U3wci@ZC`UEEJNPEE!a=r? z3JojRD%@zFFwgTddy2!qU!aiQn*9ib2gv__;{_Btnb6uI^P$3O?=-UL7;mY`^T)+a zi$YGP_v%>_osSfTGSux^;$k5Jh4Vjtn_f-?3}A`#sAU=9ACm-DmR|66_%>L7?c%w0 zuV{>uA^wJ_x&jWb=F|F8!xG~Tln#2|-jaIADq`*@ ze58Wa_T!6qRO>T^y*Xm{{#3g;<1Yjg0d32=34W%sXLv!;5$iD?8~C2Q_>p8;CF=wt z)zStKnF#vgyN=$h-9va0Z@`;nk>T-?y2S+ruYXDFG-eJJwoky`K<+fZo2L;_q6Q&I zl@PCWEHE5-o+pA+M1UXt$nESDCTVK^;9j@>WB%-b8Yp*nrI6_6xvjhgHib>uyl>oC zD(#RiX z`KUNNmpVT-aP}SH4<$hb*gAm#V8bL!8WY^L?uDw^?o1~{V1kXCeWyEMI=FE9kEO_d zY1R5o|H}ykOlXX^<@ZM~V6=DezPK`7SqWCSZ!H80(sE+djO?C(pLcnJxhgt#w=DCV z$&YLzcZW_l=Q+6vut1C8pO3_HBGm@k7tnHu&z4c!F@PmDXolJNRhaz*6QNHv7-CXe zWCEh(^4NgQnL@3)DiEz|@5zMU?skibuz7xrH=JS5RA)B-<3{lm6CtJP?VL?2!&9F( z(f5~#(8HeTZh$WYA9hn?ge?^Q<6&H8ekce|vvus8dB)_=`bn?R)*cYl>=e`fWu3TB zZk|D@5tfbJ5+MRjfBe57%+Z6ATIxe&k+z{`k`W)`><+DVzG)O+ugN656lK(ZJS2C}ah zwx?WVMZkh4r|i|VJJKBqU;9QZ%za?R5Zsg%EZS42}6JFlwaEEaNRY+Ypv*% zt(t~Pf8QsI+Sb|Hu9en*G;>jEbzH%Ky11FZ2K0 zQ;<#Yx@&h^Ug^Icmz9;3m;FEYUUB*^ImJ&*2D~x@aGG;J1sZe1QMP&;`*^P`WV@Zs4 zklir$eK%(8{0?>B-}`%ifBtbEhxd8Euj_hU&+B#01eTtPrGE`QEF~?a_(uKo!S8&F zZs9`;3eiyd^MGNb;uQ}vUcF}V2O91*TbPZzm{#!Ei-I#+OI zSo2k~ zi=@H5>RG>?Q48N4)7_on>?wxaJ(`2~J7e;iUqIflbtd!9GIevJkzoZ|w_RC2NC*&+ zl2N@q>E|gdaQ1_MQ$U%G`2NKw(_IbImE4LNr*QAu(Ne)cdzHZT1w+q@joiE1o20?y zsk=bN?zZSwaFjK<+y5O)GqAP1*=bW~Q&GNCO#(K!o44hr*#^#90XRfLH`?ZcRUi%# zZ9@rJr4|IEiwDs`@b67+i#t1sxsgR~ae)8qir^cs6+W|x3RmSQ3EdjyzO6!g4R5xu z_Xl@(N4JvN#5Z!CtQI!ZB;W0A+3ZBB1HfKn#oof&-lWa$YV|CtBCv3lnh9)BWxm>%sUsGT=huUB*-uIV99q~?Yif=QU| z3fZQQnhruz*8;#C2=-_(iR@=*>9+FbYk_r+c+@g7_Wh`F(pp}v%A`6VasQXGP9s*m zcFkio zyK|btGdA10`3f=1INmQuR~1KK*Skoo!qS*QcuTiG)) zz0&S7Tx>15jZ~6P;LWL`ah^UB+veRf6KtE>0l!(xl1k2W8kPj@PHQ*?=r@fV!W#{W z*JiI2vQAwrM>*}d_$vp>h`9&V{*z->%d*1ek*rsgz3WlX_DhEsh?u^FJ+C;%G-+4P z!7s=rW}I~MW0V4aihM|sL>&7Ed9?o{<&b8N0R?^*+iQakRJ;+MuiZ!(VR|=IB8}$W z-m51~{oKbVm!=t7VfsZhP16xJ1;;tGpHnCb$4xsmmJF8M4E<^aUca`sRU+k0|A=8T zjEiSKo%=hU{g?#M>e!2T(J0Z1@=}&^(IX-1*57Ztx1Xzq-rWm)sQgOY{)@ZvW(QmX zQH6NxXE(aiC3+8DuuFBJoSS!Fifaz#vIv-`O4SZ)}VhsBZXa~SS@}{ZtO77r&$sGoU!I^qgfZ{75XOq zOg3a%06#ZuzxxKsFMA|-DDg)k>5}Cxu77rE+;ibm{dba|0xH=hckEWl!iR-XkNZDF zy>xEQxzCU0bg0ewdO6-q(MTc;$mj263ZQ!Ad@Ko8xG>sPR7QI%?NQO)S|-{}R2H9Q zkL-Z_oU;ImY$3g;jT$_ExiUcp@)aM+CDFG4@Ir;+zqi$UvNi<|H~cihhc-)A-|h5R zJxTnoeh)gF^HM#Y@?2A1G)pinMp*@(G!gS~CnH}wQO%^ytGZ8 z0e2g#d{GS~H@OcHfl0j?BdRVRZ&ELKC-vrFUzjoz3uXL#;oQ>gge+aW_csxxM(diV zndWwX8J>aMo06{_t=RtGr~-cY`e~d(>~SrRg}UC7yV2# z)4F;r;qtC%+{e=iTQG7{YJv75S3;%_kTiEQgO$ZJR>$#WEmX&`m-U099fEsD?+Q6f zw`fa^oMYXe)yiio{p`0fqRYxbm3*h@RUXpMB0)VoRA(~emuy)f`<7~?V(ut*d3j3g z>fEh1jp>-wf+o-AHuH}USMxVsOw~+fku$zbiP4UEj<=a}-)(90m5AerB+9VCfOsc! zVcwsCXKv<`@4q6Bh)?RZnWyatgnFdxS+!8y;cH`U#fp5AkwnhyaL>)J_$tzE8&#)= z-*HZPF&jLzq6R6PoB$OWa3;aKC@Q7E5ACCaK^4IGD9PK z6HZJhWr!q>Smi@ksYuz6k5paca8%-y*yj9)l!EmF>2|+ZhoM{*<=KU7ejvGYG;v$E zeZXc#F~gIlqg*bpVE#GXhlIH9^EY48@nI zn*w*h%RJlqGpPlekz5J01MZ5m%Nld-Q6k!tRdbc81r9tYFhzjari?O&XJngiYlUp$ zbRNm$>9U8%g0`5z%ytzZtvow{)kc`kOT~o~WoqGxHYo+~!}~VITDTItcl%xvQ{>Yk z+j>{lQzH!e5_JL6Pa0kkYvD7eWq%^F>|XJNkatc58XZZ*)@cj4f8Za}Q{#Jwh9!w- zUueF-&WCZL^J@(W-n?fQ+rGjV7uxht z93nt>#{hFmFc?dxanAjJK1wTxov{?@?%c2v>2r}Q10@+p8)`H6ihc`Klum_ep}mi$ z36f5F78^@MMq@$2(NJkqV%u+!LM|LVY@0q}3xO)ijhhu-2-5;$Td{b$IPa8EJE&>Y zq^OpobX>`Ae?tLzUqlSd9E%UXiUz)@)^8SWo=FhyGncr9ehc%Uwp$CAG`*WMuicud zC&P-`j2FTk7Pt~%Wv;qy=5w8oVL6IaEHs;Ns?;$qa;7b+&AfjG;*6IZWNQ_yld^GE z%q=$6-#uJBb!#`x3pLVK%oc`O$o&;~4@0yyo!~f~ASobTqW*TdFx5!g5jbN?>0V85 zEV!aQ`O|V>CqGgC^Hf%EL=E32C9c$JZH;qr0g~YT{_#SVPPS;x$v1>`Uj(2_evnDIqTG_!R$Q+@ z2uzzQgPWdiGcO5+*H2uy7YVoup+{rZIF;*#scF_N_!qps7~+UKidVca?f!44SBhr7 z-)ls|Aqg;V($r&I$I*9vr$L_bx7Hq<#AF>mTq{@svV+B>)ox#8A_^9 zUt2_CqJH~F;rrII=3AE))XRA#Z1^q3Yf{}Kk3ps%Sa^v-eMek6^ntN>Xwd@GfP!GD zDs4o|D+v)5@{AG><+xvjJ{g%Nkf0?{Y<%5O{oCiFjPHSK4j)}($L8%%=b^nkFH!7@ zB~uH0V)4<6ou7V36=nCw`5|bX!{qyAVF{mQlxO7;7@{At-53!=%j;g75*rWzRlVYx zW|=qhu;X3Io3RCMAmViPuO^*_g#5U=;;+`M@3`2!Ka4qvw9(C@5!TIfy{QG)+Y?jr zCo;I3&m@S>gUXgD@W*EJoRcq1FAl&EW#+)T{;mLtd&jPnp4Oc<9v&LrQs>Af6X+lr zPAow~M3ftSvgHIb_HHSj;2pUCvBgxS?{PV-m*nB{YZLIlKDj`bYADvC z+Q%0F#q{k8D)GBcB=NUkW5`q_Ofld~*o}c~vz(iNL}`UubVcmQyIVM_SS=PBNepg# zw|pHnC|eUZkl1Fv9olI_^tduf3{;p|zV)&GBB`#h-Ke1Ha?^(@1iWz0405|Rm2&Y# z&FeRuUVIVir)-T7%6tF(dP?;;j6XW;7ia`T-BKw>yo2fG_J5h+JI*m&#J{)x`p(C% z9OmP&uN>bz-cy!Qt?`__$4A#2)3YephIGWzdPBNDkudPxVD^!DIV&fR`btzr*fE*dcpnVZ^F5!WS3o@= zQF6!T_llWeta%T7ve6kIo>31&co_A6ak`dhh4I10NTz$1_@m56dM$3{eW?&FRWWL` zl(MoNjY7qWP@q`-PKJ;u%bbhv7ouv!v~T6X>WdKO=3DkR9FBI+K8??_BAs**ShNbf z(jeYBv@D{DaWfJi43-4z#M`B7-F3fuDZx`n24-`Z<)-_o&jcbjl&qb~fO zEX%C*sPm{z@NEB*Tz|P| zzy8v&ZM=YE4;H)8C_&jvwVF_Lb~MSuH>i7SQvAf}9w0`0sP?N|L$R!{%}jwpo4^;$ zBT!#Nb>S5rGyjGwLHCz0vB2%e9{-r7-oHAPO1{8s$ycX}*jboI-IrVl^@fUb@4R95 zkU#hG=Dn(_keOdl>*KA?vm3!)kJ=UL*37?L> z^RY7@EW0{8BzYuM9z}H|18=z^N;VTC%&^Kg+_HUg=w_fQ_|+RyVD|3TMZMZs;e5F+ zvnN54iXZl~YbbM{jI)-AqdM3yAGIo@V#|Gw5sSG>z92=^^bF1s4K$Jg{2?d(r8h~> z#Y8&X4GTVhEqp!)%eBk*xJ{EIT}bjXQkb0&F!vG79NDD_y@Q*W!`}s+Y~6W9IeXj5 zP_>VL`_(h!SYZ&3?PBz}8kp#Un$TDfw%m+FOi#P8)fyM@?MDC7bKkF*iscPXJ<1rU z)lP#FOEf+pmnN!g%$OwredhhOAN*+3SmWvPRFBj-XIIijM;O-JlUAdk$omPEU+(kL zeHS+QNkR1O*!f3>hneT)%X}uc10Sm1?h@&UNXLJkdc|hmZ+m^(PqCkG7pXRLo*y|E z);{!=@2zSxE@IK|!OpT6YTjsE>VxvM{sxiB_NJUD>qnIP}P>0 zW6}4xd6i4!YhXSd6W8gs#(53;(3yMnl26o&vrNA78KRF#B}-)LbO_BTG|b>|3E=Vr zY4;Z<66(ksLa}lKZ6ZDfJ1KqRK1Y70p=*;JduQ$?L^fu>7QU-G{etss*)`hl>#8%1 zFv9}(m2W(-0y033CE=)*1phJ4o-fx}^4gh{6_>a^wOVK6x-VArIwbuKd*ZBkovBuz z&OzMn{@s~1x+rz;RLN5l>5IXIuz=+#zBDWPw$9IHQJX)|M_}_0la*|6uo*9rS+1Gd zcI&4mnz47suKIo%WmuqsZD}I8Uof5lKn~Tvjm+VZ91b`p^GH=FT1@%2%_oD>r<{G4<`qLMiM+&hAPoS8A)Z&V4?A`?7O? zTWh!Pn0=7Oi;}Co?QWpF<@p^Ds}45WFFPMD3+n;!sxLvz3M~F`2t&{{^s{2=4Fik_ z)RqeWjQ@;29-Q^X-uh3mY;1#@4E@kHF#&fN6gJ-I$GBI~3te`IT~2H@r^m*&YE1ZR zoVDuuR&c27u+r}9`H#=84-6C}IFZ(KV(Znw9WSgnoJ|$ep+(v14qhi_Pd;#gU`hb7 zL?MQ2NI|ugjo0#ERnyIFj|ssf^-u82VpVi#V3BW%>=zj~?#TBwA=2{q37evEBzB{R z7IvD;;hx`eGIGSfksMvAH|Hx}!T|)w9vhmt{l&Bn?;IecrWTACX z_nE+{)K+O~CJ|U%DwfdCsF?vrSgD+JKrExrZ<$_#e;{BU|=63d|&p| zm&_-y6-e(0Fh@R!*km`f@R-`)8M%oG>Am&VjcOo)9u|LBNiWfY1Pgp`uux;4o%t=3 z{Xc&0)Cm^_7dCGHLpu@Z?_Z;)x!5v%ufFK{Zvdt%uFn60^E3~P%GP3Ky}WDFBS}J* z&%&@AU_(FMHc+>(DhRvuK_CQ8FDw%z5~l+l3#rl}pT>UPY{LjwS_X?Pk&IB(&Qk~d zw4bYtfV8heRvVFbj2{&~iS1`m(9bO6EcY@heKq#adL5r;*QHESqQfzjAB{6100r~E z5fZX(Fqv_Gms8jwmaUkdd>A(x6L2M%#{0N>$aCv6J4`p(1MRnU5rv4J#cE2LrToXB zXf9tJv5C%T2O*>Mw2Q3;R{Y*&@&|#mcIHP@83*9=etPHjS5l&YK-wK>EI;$T8}A#g z0fxiw(KafqSNC%j6is)b94jgI!L&Lga#4eNB<&Iis2QKVgx_ybP+ z>B9I%(@N#CKa&R7=ylQN&SluGS67Eh6jokEPS&-v`pq9Dd`b@BwONkm%wA>q~sVy zV)hQzgm~%dd5&KFg3Pa5D<(vzKVwF53k!Tqf|VH3IQ_3`@hi1F=$zz~a(wHoxWmf5 zJ?@zI+A5H7791qZ+Zm14NFk+d)~i+dWFrrht*n^35Pp z%LtX~Q0O~J`hu%xJSY`wu>Nzm?(P^+BV9bWh^ST;tRX&IatAv*_A2*jz!+y)@Has4 zw^3YnQ6Z5h>goq1GGgApbQo!F49diz{ufh+B{#j?j6mX~7HRAcF=K6gz3Ag7LDtJN z9kFsQyIfgYvvmpNb$PgCFfC?CP}@|8>210WDc^ zOC~MVz=yafO#AX&EoStg-kPCX6qe>esJ;YPV?vG0jZGC3X`F6&OVhz)xkg1srbUJG z{w{i;b)pGZ+}^SiTt94}-@6?yX2qPZOOvj_}~WN8aTCrGg#o_(CN`Z_o7o{?09DH zD?$R5(p*ECnr}}gAQ>_9bQtHjF=N%}MxaDssAu~SBo23b^_PYxGcWfwfR3<9kC}G< zK~W@fTg7xRD5VLrS@6f+ZLg4pH(gbABW-GXy7RG!$0IU~6_R4n8l_&q88 zXWC>b_I@D+%1Vamsy?y1Btj=uEM`4tSeYOA#IBV8po6yP2&=f&W03Z!5INwM=edl| zbh`6Q7%TeQ8S5o!6j5uY$L7Xolxm~cKGDNO)vDe-vD={|eBXC@(E(EdV5BSuks7Lr zk^xga!)!>uHiH_u=9axfO0}HN z+@=&M9|H`8=m{aAsz-oq?{$-I&AKQ*I>PPWt4P=(Ijuj|&?sfV?ndrd}!i9W+>t_JV$b zLp{5PAS6bJiKJ-UbLT)KolQCd#7D>ilh3|a>jxUzgGk4@u>96EI>H1!Ci$=S*&)>y z`W!>=-}U#roSi0g(y^MfAfDK%{XTiLacC5(e$Mp{BZf`gbsCVhSLBR742ffh6#O+C z$f5UEuN}cNh2YzQ7ZF3$O(u92WaxeE(Miu=9I?VwIt&{{R2e5>A-%w)G(p+7(fdwR ztukLtje&6cDaW{okG^L6&(OD=2Mb4j#791tw$~NjkF$gayPfkG3}bz2e;shJP@3v) z2>wydt1uTmec?u#+mlk;Rd_4>iEz5eZgFO8@%sY}{%=Rx1@sd9J zOY{)CW8~YWdOIoQo+=zaW@)JspR`Sv8prQ#bSr}u4dHwyTRit%bEzb~#rFFjh=UT` zTxW&MA_2bJ7jM1nkYs02g$(t0Kd}Q1vnS6l3H{6T1hc28#B1Dw-1?Q$0(bsJU=Fw& z0OtAE0|XVDpr&?*XDH{)CxIP-api4?xPXRzA(*Am9glU#&0l7J{*4_xYpwc6Yicx^ z9TKgzD25QlE2Vwq5ONzk&sZ+xWYx9!|*xg zrKRgc0sk?|IoFPpN~U6$JGiloec-b}lJwyyrO*7%9*Az|h?4fN_kH8Ieh&^pO9EFb1ABg|k9&?Rs2xmr&?=fqjt=03aRZN%_$579B zWj=GojDM*K!rT}$#+hRG-YIKp?Zw*Pgn@2) z05fgfy#boim^+6sS})0c#dC*Ms%H}ujdTulGagN5Qo@0$xh;d34EscaGpP$4Aj^(l zZKd}0elsS2e)!pnJLz(GNnDcE^}n6LuN1O9w|~yyf_qyRQru(u7?ipoN+JDH2M+iv zE5xpzapfHr*nBf3%0Te39(RW(eKZDLXk~sfa*Ym?Zx7|YGgysZHtM58xetmfY+0PF zB52tkpe}@-Fr$L48#)Qk68ulP?y$ha7bP;^(OecF=_oLVu6Nvibq%99WM2z|6d$O* zgO3zbYkB{~Zj_GT0^J^Y2pq+lFET(#Hc#wGIdp4`KKcP=ZV?m8LpjrvEH@-fm*k62 zpWEj$*u{`yw*t0{vmWK15&I$qKd+F2>;^E%lkg_JzVSh+&VcK(J~<2o34mawzjQ%t z7!*ZYHy(O#37O-h<2R_2FMWZo1& zfnKfiumq3+A?{bbg(k4E3#5~9H^_hI$T4r`UhH`P=HmqbEDJk4E2}ywkddW#+A)tF z5*Nq_T@Y%l@IF;3y;e^~Zi|a0R$;dG0iv0+0aoVMPz0ZE#Z0adg|I`_o8YEklem*N7=aYfZE_5m zNe`VKGI^VwVR#60#2)jP=(O@4w#SU%hV1E}(dAy0J_%>#k{FaMEk*`&eHGT2oV^1-&aPN4OKlOj25-NX!_1Z>>(r@>Dvo0faAQ+L z2BiXeV_+aMCTCgKlgXlY01VZvHd}*w}((_|V3hZ#V>}tj1g^j#mU$Ze_=1*Za zhw-3`AO%s$ILbE>@B<2yQrAAX#snFLvRy*&xyXYzAZoC+zI}HKn%jMl#4|k?`%O3u zAu%bLl7^WOHn4xqc1=JxGz=D+@3Mn*