From aff4f894cd5983e81ca028772f80a1fdbe71d72c Mon Sep 17 00:00:00 2001 From: Isaac Kinley Date: Sat, 9 Nov 2024 10:10:22 -0500 Subject: [PATCH 1/5] Initial attempt --- R/generics.R | 38 ++++++++++++++++++++++++++------------ R/internal_utils.R | 20 ++++++++++++++++++++ 2 files changed, 46 insertions(+), 12 deletions(-) diff --git a/R/generics.R b/R/generics.R index 12d93f6..c38192e 100644 --- a/R/generics.R +++ b/R/generics.R @@ -490,6 +490,7 @@ deviance.td_ddm <- function(mod) return(-2*logLik.td_ddm(mod)) #' @param del Plots data for a particular delay #' @param val_del Plots data for a particular delayed value #' @param legend Logical: display a legend? Only relevant for \code{type = 'summary'} and \code{type = 'rt'} +#' @param p_lines Numerical vector. When \code{type = 'summary'} the discount curve, where the probability of selecting the immediate reward is 0.5, is plotted. \code{p_lines} allows you to specify other probabilities for which similar lines should be plotted (only applicable for probabilistic models, e.g. \code{td_bcnm}, \code{td_bclm} and \code{td_ddm}). #' @param verbose Whether to print info about, e.g., setting del = ED50 when \code{type = 'endpoints'} #' @param confint When \code{type = 'rt'}, what confidence interval should be plotted for RTs? Default is 0.95 (95\% confidence interval) #' @param ... Additional arguments to \code{plot()} @@ -497,13 +498,14 @@ deviance.td_ddm <- function(mod) return(-2*logLik.td_ddm(mod)) #' \dontrun{ #' data("td_bc_single_ptpt") #' mod <- td_bclm(td_bc_single_ptpt, model = 'hyperbolic.1') -#' plot(mod, type = 'summary') +#' plot(mod, type = 'summary', p_lines = c(0.25, 0.75), log = 'x') #' plot(mod, type = 'endpoints') #' } #' @export plot.td_um <- function(x, type = c('summary', 'endpoints', 'link', 'rt'), legend = T, + p_lines = NULL, verbose = T, del = NULL, val_del = NULL, @@ -541,17 +543,29 @@ plot.td_um <- function(x, # Plot indifference curve lines(pred_indiffs ~ plotting_delays) - # Visualize stochasticity---goal for later. For now, don't know how to do this for td_bclm - # if (x$config$gamma_scale != 'none') { - # if (verbose) { - # cat(sprintf('gamma parameter (steepness of curve) is scaled by val_del.\nThus, the curve will have different steepness for a different value of val_del.\nDefaulting to val_del = %s (mean of val_del from data used to fit model).\nUse the `val_del` argument to specify a custom value.\n\n', val_del)) - # } - # } - # p_range <- args$p_range %def% c(0.4, 0.6) - # lower <- invert_decision_function(x, prob = p_range[1], del = plotting_delays) - # upper <- invert_decision_function(x, prob = p_range[2], del = plotting_delays) - # lines(lower ~ plotting_delays, lty = 'dashed') - # lines(upper ~ plotting_delays, lty = 'dashed') + # Get inverse of function that computes probability of selecting immediate reward + pimm_func <- function(data) predict(x, data, type = 'response') + # pimm_func <- get_pimm_func(x) + inv_pimm_func <- function(p, del) { + df <- data.frame(data.frame(val_imm = NA, val_del = val_del, del = del)) + f <- function(val_imm) { + df$val_imm <- val_imm + p - pimm_func(df) + } + if (f(0) <= 0 | f(val_del) >= 0) { + val <- NA + } else { + sol <- uniroot(f = f, + interval = c(0, val_del), + f.lower = p - 0, tol = 0.01) + val <- sol$root / val_del + } + return(val) + } + for (p in p_lines) { + y <- sapply(plotting_delays, function(del) inv_pimm_func(p, del)) + lines(y ~ plotting_delays, lty = 'dashed') + } if ('indiff' %in% colnames(data)) { # Plot empirical indifference points diff --git a/R/internal_utils.R b/R/internal_utils.R index 15c0aec..bba2c09 100644 --- a/R/internal_utils.R +++ b/R/internal_utils.R @@ -34,3 +34,23 @@ get_candidate_discount_functions <- function(arg) { return(unique(candidates)) } +get_pimm_func <- function(mod) { + # Get a function to compute the probability of selecting the immediate reward + + if (is(mod, 'td_bcnm')) { + frame <- do.call(get_prob_mod_frame, mod$config) + } else if (is(mod, 'td_ddm')) { + linpred_func <- do.call(get_linpred_func_ddm, mod$config) + frame <- function(data, par) { + drift <- linpred_func(data, par) + return(pimm_ddm(drift, par)) + } + } + # Get a function where the parameter values are fixed + func <- function(data) { + frame(data, coef(mod)) + } + + return(func) + +} \ No newline at end of file From e3b810821c12bd5c100ef2c070470a3111fb4598 Mon Sep 17 00:00:00 2001 From: Isaac Kinley Date: Sat, 9 Nov 2024 11:38:04 -0500 Subject: [PATCH 2/5] Faster with grid search optimization --- R/generics.R | 52 +++++++++++++++++++++++++++++++--------------------- 1 file changed, 31 insertions(+), 21 deletions(-) diff --git a/R/generics.R b/R/generics.R index c38192e..0d23b55 100644 --- a/R/generics.R +++ b/R/generics.R @@ -491,6 +491,7 @@ deviance.td_ddm <- function(mod) return(-2*logLik.td_ddm(mod)) #' @param val_del Plots data for a particular delayed value #' @param legend Logical: display a legend? Only relevant for \code{type = 'summary'} and \code{type = 'rt'} #' @param p_lines Numerical vector. When \code{type = 'summary'} the discount curve, where the probability of selecting the immediate reward is 0.5, is plotted. \code{p_lines} allows you to specify other probabilities for which similar lines should be plotted (only applicable for probabilistic models, e.g. \code{td_bcnm}, \code{td_bclm} and \code{td_ddm}). +#' @param p_tol Numerical vector. When \code{type = 'summary'} the discount curve, where the probability of selecting the immediate reward is 0.5, is plotted. \code{p_lines} allows you to specify other probabilities for which similar lines should be plotted (only applicable for probabilistic models, e.g. \code{td_bcnm}, \code{td_bclm} and \code{td_ddm}). #' @param verbose Whether to print info about, e.g., setting del = ED50 when \code{type = 'endpoints'} #' @param confint When \code{type = 'rt'}, what confidence interval should be plotted for RTs? Default is 0.95 (95\% confidence interval) #' @param ... Additional arguments to \code{plot()} @@ -506,6 +507,7 @@ plot.td_um <- function(x, type = c('summary', 'endpoints', 'link', 'rt'), legend = T, p_lines = NULL, + p_tol = 0.001, verbose = T, del = NULL, val_del = NULL, @@ -543,28 +545,36 @@ plot.td_um <- function(x, # Plot indifference curve lines(pred_indiffs ~ plotting_delays) - # Get inverse of function that computes probability of selecting immediate reward - pimm_func <- function(data) predict(x, data, type = 'response') - # pimm_func <- get_pimm_func(x) - inv_pimm_func <- function(p, del) { - df <- data.frame(data.frame(val_imm = NA, val_del = val_del, del = del)) - f <- function(val_imm) { - df$val_imm <- val_imm - p - pimm_func(df) - } - if (f(0) <= 0 | f(val_del) >= 0) { - val <- NA - } else { - sol <- uniroot(f = f, - interval = c(0, val_del), - f.lower = p - 0, tol = 0.01) - val <- sol$root / val_del + if (!is.null(p_lines)) { + # Plot curves for other probabilities + + # Because of the overhead of individual calls to predict(), exhaustive grid search + # is faster than uniroot() or similar to invert predict() + + # Call predict() once on a big grid + val_imm_cands <- seq(0, val_del, length.out = ceiling(1/p_tol + 1)) + grid <- data.frame(val_del = val_del, + del = rep(plotting_delays, each = length(val_imm_cands)), + val_imm = rep(val_imm_cands, times = length(plotting_delays))) + grid$p <- predict(x, grid, type = 'response') + + # Split the grid by delay + # Using split() with a numerical index is faster than calling tapply() or similar + split_idx <- rep(1:length(plotting_delays), each = length(val_imm_cands)) + subgrid_list <- split(newdata, split_idx) + for (p in p_lines) { + # Get the val_imm producing (close to) the desired p at each delay + val_imm <- sapply(subgrid_list, function(subgrid) { + if (max(subgrid$p) < p | min(subgrid$p) > p) { + return(NA) + } else { + return(subgrid$val_imm[which.min((subgrid$p - p)**2)]) + } + }) + # Plot + y <- val_imm / val_del + lines(y ~ plotting_delays, lty = 'dashed') } - return(val) - } - for (p in p_lines) { - y <- sapply(plotting_delays, function(del) inv_pimm_func(p, del)) - lines(y ~ plotting_delays, lty = 'dashed') } if ('indiff' %in% colnames(data)) { From ffbaffd1c3352f539073b7d74f5e6091782ef1f0 Mon Sep 17 00:00:00 2001 From: Isaac Kinley Date: Sat, 9 Nov 2024 11:40:41 -0500 Subject: [PATCH 3/5] Update docs --- R/generics.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/generics.R b/R/generics.R index 0d23b55..0072bda 100644 --- a/R/generics.R +++ b/R/generics.R @@ -490,8 +490,8 @@ deviance.td_ddm <- function(mod) return(-2*logLik.td_ddm(mod)) #' @param del Plots data for a particular delay #' @param val_del Plots data for a particular delayed value #' @param legend Logical: display a legend? Only relevant for \code{type = 'summary'} and \code{type = 'rt'} -#' @param p_lines Numerical vector. When \code{type = 'summary'} the discount curve, where the probability of selecting the immediate reward is 0.5, is plotted. \code{p_lines} allows you to specify other probabilities for which similar lines should be plotted (only applicable for probabilistic models, e.g. \code{td_bcnm}, \code{td_bclm} and \code{td_ddm}). -#' @param p_tol Numerical vector. When \code{type = 'summary'} the discount curve, where the probability of selecting the immediate reward is 0.5, is plotted. \code{p_lines} allows you to specify other probabilities for which similar lines should be plotted (only applicable for probabilistic models, e.g. \code{td_bcnm}, \code{td_bclm} and \code{td_ddm}). +#' @param p_lines Numerical vector. When \code{type = 'summary'} the discount curve, where the probability of selecting the immediate reward is 0.5, is plotted. \code{p_lines} allows you to specify other probabilities for which similar curves should be plotted (only applicable for probabilistic models, e.g. \code{td_bcnm}, \code{td_bclm} and \code{td_ddm}) +#' @param p_tol If \code{p_lines} is not \code{NULL}, what is the maximum distance that estimated probabilities can be from their true values? Smaller values results in slower plot generation #' @param verbose Whether to print info about, e.g., setting del = ED50 when \code{type = 'endpoints'} #' @param confint When \code{type = 'rt'}, what confidence interval should be plotted for RTs? Default is 0.95 (95\% confidence interval) #' @param ... Additional arguments to \code{plot()} From 29d8f4354a1ba72224e966d803a7821dd5897ad1 Mon Sep 17 00:00:00 2001 From: Isaac Kinley Date: Sat, 9 Nov 2024 12:01:13 -0500 Subject: [PATCH 4/5] Update docs --- R/generics.R | 2 +- README.Rmd | 2 +- vignettes/visualizing-models.Rmd | 6 ++++++ 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/R/generics.R b/R/generics.R index 0072bda..66af393 100644 --- a/R/generics.R +++ b/R/generics.R @@ -561,7 +561,7 @@ plot.td_um <- function(x, # Split the grid by delay # Using split() with a numerical index is faster than calling tapply() or similar split_idx <- rep(1:length(plotting_delays), each = length(val_imm_cands)) - subgrid_list <- split(newdata, split_idx) + subgrid_list <- split(grid, split_idx) for (p in p_lines) { # Get the val_imm producing (close to) the desired p at each delay val_imm <- sapply(subgrid_list, function(subgrid) { diff --git a/README.Rmd b/README.Rmd index c99fdbc..ff6c7ef 100644 --- a/README.Rmd +++ b/README.Rmd @@ -103,7 +103,7 @@ We can explore an even wider range of discount functions using nonlinear modelin ```{r} mod <- td_bcnm(td_bc_single_ptpt, discount_function = 'all') -plot(mod, log = 'x', verbose = F) +plot(mod, log = 'x', verbose = F, p_lines = c(0.05, 0.95)) ``` ### Drift diffusion models diff --git a/vignettes/visualizing-models.Rmd b/vignettes/visualizing-models.Rmd index 015a321..5a9b1b1 100644 --- a/vignettes/visualizing-models.Rmd +++ b/vignettes/visualizing-models.Rmd @@ -36,6 +36,12 @@ The plotting function prints some info, telling us it is plotting the discount c plot(mod, type = 'summary', verbose = F, log = 'x') ``` +Additionally, we can plot some information about how stochastic the individual's decision making was. The discount curve shows where the probability of selecting the immediate reward is predicted to be 50%, but we can plot curves for other probabilities as well. For example, we can show where the probability of selecting the immediate reward is 10% and 90% by setting `p_lines = c(0.1, 0.9)`. For more stochastic decision makers, there will be a greater separation between these (note that this is *not* a confidence interval for the discount curve itself): + +```{r} +plot(mod, type = 'summary', verbose = F, log = 'x', p_lines = c(0.1, 0.9)) +``` + For an indifference point model, the discount function is usually plotted alongside the empirical indifference points: ```{r} From c8d257250d1299d71519e54a45a3e6a91daf40c9 Mon Sep 17 00:00:00 2001 From: Isaac Kinley Date: Sat, 9 Nov 2024 12:22:13 -0500 Subject: [PATCH 5/5] Update docs --- NAMESPACE | 2 ++ R/generics.R | 5 +++++ README.md | 2 +- man/figures/README-unnamed-chunk-13-1.png | Bin 7223 -> 7767 bytes man/most_consistent_indiffs.Rd | 14 ++++++++++++++ man/plot.td_um.Rd | 8 +++++++- tests/testthat/test-td_bclm.R | 1 + tests/testthat/test-td_bcnm.R | 1 + tests/testthat/test-td_ddm.R | 1 + tests/testthat/test-td_ipm.R | 1 + 10 files changed, 33 insertions(+), 2 deletions(-) create mode 100644 man/most_consistent_indiffs.Rd diff --git a/NAMESPACE b/NAMESPACE index e44efa2..12ff81a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -32,6 +32,7 @@ export(discount_function) export(invariance_checks) export(kirby_consistency) export(kirby_score) +export(most_consistent_indiffs) export(nonsys) export(td_bclm) export(td_bcnm) @@ -51,6 +52,7 @@ importFrom(stats,BIC) importFrom(stats,approx) importFrom(stats,binomial) importFrom(stats,coef) +importFrom(stats,filter) importFrom(stats,fitted) importFrom(stats,glm) importFrom(stats,integrate) diff --git a/R/generics.R b/R/generics.R index 66af393..fe4220d 100644 --- a/R/generics.R +++ b/R/generics.R @@ -546,6 +546,11 @@ plot.td_um <- function(x, lines(pred_indiffs ~ plotting_delays) if (!is.null(p_lines)) { + classes <- c('td_bcnm', 'td_bclm', 'td_ddm') + if (!inherits(x, classes)) { + stop(sprintf('p_lines can only be used for models of the following classes:\n%s', + paste('-', classes, collapse = '\n'))) + } # Plot curves for other probabilities # Because of the overhead of individual calls to predict(), exhaustive grid search diff --git a/README.md b/README.md index ed00740..b2ad952 100644 --- a/README.md +++ b/README.md @@ -215,7 +215,7 @@ all of the following models are tested and the best-fitting one ``` r mod <- td_bcnm(td_bc_single_ptpt, discount_function = 'all') -plot(mod, log = 'x', verbose = F) +plot(mod, log = 'x', verbose = F, p_lines = c(0.05, 0.95)) ``` diff --git a/man/figures/README-unnamed-chunk-13-1.png b/man/figures/README-unnamed-chunk-13-1.png index 38389d9731a29f476d50e88167a1a40cc3cdf4de..c4dc9922354ca1765e14272dbcc28538233df1aa 100644 GIT binary patch literal 7767 zcmd6Mc_5T;_xBjfU@&&2vhP`veIFCq%2LWs$};v*gsfx9ng+!~ma$Zl?1sixqzDs*QErUuN6yo?|ah}p3hxZWxRk2ImIQaX2>|@JYs~gD`-Cn>q%A1O72MB8G^= z5P?rJ&JBcfQwOAR7(f=_;BZ7D4xp2XI5L@tQ77Wu{>oBlL^7ES7#HPBj*0>_$v8j? z^H-1}3owa0CRV^C7*K{b4?rN64$7PAi(jb=2y{W+NKe-?ICm`%8bnl0kE_gBE}zn> zg*CZPZ|<~~hwtiX)GqtHBaV{Ew@ffadPX9aAVz7{=4rLfEALD3X;LEaXwYsfCsR2y z>kRcM9sMSi7j^7@R6P|U?85{|Q2UNg6f1v~GCToheNNyLfm`X6cl@NZR$&qBr?tkbp0(3tl%}$p={< zRCI9ZI*K}{D}MH9AOAhMC#Ns6YBM{sRFM0>k17PnzljT3I^oV}7IZz@b#r%4ig7y3 z-}SM<-u!l9FYmg&sBcq&L!1w!vY_FY6Cd}|`o7dZw@Y5IJ}A}**2cxZZ;A{f*ub`gxrx zh}J*q>xcQ05Th#q+r}HUw|C@}Z-3HAIwWh>*DGr!@!y`Qp8q|+wpnw7_&dveSGS+y zP5;K6?_9UG1LG_9)mg+3KL!$Risswqw=`=wLVQA*&WktzCPQbl_;X161?2JDqnp2e z+L{mJC#M4A5IHxb zf1NRU&MwQ6S;|$h^DV-6rSiort|`=?b^?muXo5MVFS;gC=4q45v$VKtGA85GisxYo zx3sD}NuK3m^f0a_h}g4bccpFLJg$Pj;GY1W?7lqdnR1psdEX?hWXB({j{0G#upk6E z)V(n@!f{?VRN1|)9uAnQ$Pgs_z5ILE>EY#@+pa}0bxmZdeOk0)u{Ge{XIDcDPfy;r znQye(Kbdhb_z)gs@NL%|Z@0`gEJB*_*cd(U7l?BbtBk>KrXBUH9@?9sYY^kEA}rQ?8og?m?z_I@sQ z6g&{{k>&pU_(*kM{o=M&hfD--Se{(*X6WzjpQ|84VsczHshlz3)U&bsSEg|=%o@v) zXtIxzy#pjgU5oXZ3+wk3QnE>)M{2jFNFDE8))jo&HguQ62@g4Sdbaa1bbX$zW_c4?BTXIXUGL4HFN{*sldY+^%LnJfXvJ`;_g%;1Mxb*qki@|8zYDz%zliZWbeAI(3cb&r;F$UU- zCUxUSgzYrI|8fv8M@|zknIT)XBn18)v5q1~; zCn59qS1TALPbV;E@3MPgWr3KZ&*4+j0$b3oATm&2s=EC%s9A4O7%KF+ylVAY0HZyV z8KqM#6>+Um;D(7}wHie>-VPBtJKBH2*MOa0mXjsR*mY=QS^h)YEuoUaPN0C~99I`; zO>OLA*k!?)$ZYS4`@!L7)~_JVu~X4=Yn;P#%%AQWNew*1xqO&!$*O?&^PgDHv(;Ts z!w45js`&@*9RiZtPDD~Z<4t&uiFCmEBka@<)i3+uFFr0u;QElcdn^8ftEdCUdG)MIm2J-2E7=x0@<~B4 zC9{U7=2+Kqi>Ja?k{Ae{<15_s4BV$_eby-W+OM;@Apj#&dx^`vK^;cy!f6Cf_fJB6KhCVpL;+nKv@~8&A$J9P|*8_k(7LIn>Lbx zEPPc9K4B}J!n~LV%UL!|(ss4svRz~U$;S6jt>L(4eS5dc?`a>cc829f>^usF!SO)v zbr>;gB;ua=uLzokm4n2mmZO;0~(y zR6B*8tVLRy=jhCNzVVuFfFrEWzbl{PPGD7DJTP-o;po4QT+Cy3;mDw2m7Sbw@zx2? zo8;3~+X(v@U9?mpqM}JsX<0#VL=0RBS}&Z}1s z9Y4G-agDn>SkBSk?I;ih8{1~3&VjhH*M4Q_c{H?!r~UZQBag7!ob%K%Hh0NDJ2z^S z>W11|Y4du9o+<6*A?QWA#~MXP-k4^+FF(%s&uUh|&EONK62Gio-)b{=GUgw44$}>I z5^mlpj11y%BqLf5x0ml3T{T?Ibyu6L@aL=?hXlqG-aX`3gh)aVAw#huE>RQ0{A=dg z%lm56@j2#2TbY|`T{4u4lSNw!HM(TT)=HK=RFwPX1uR}7oF(*GH?RZUMd#RMFhaSq zguZyWqcKSTGj@MRuAa_%JLpD=xq+P*uYp%~SSGHB+IQw(uWF|9TG$+XtjqAVY zGjY!JrOzA})3ZPkrV{$s(47a7T5VPDScRj?QfQ(%Njs9QYL9$-G@?I?a^aGoroK^} zBo}n6R`~uPA3}UlGjY3&J?qxEl&tDz3}@lu<)Mlf3XckMt+E=cZ&(T!0Z?a?zU(v+$Fxn$O)ropNnRr76>No*Hw+_v@poQ>O09H-8KB-&}oaznJH*5K}u>)7-eZ z=Uj6R)xPyNrcWgdDo#DWO({SA9v9+o!ivSy0sy~ct+Ra`RqG2;*oO)ACBI;Nzd|QX z;K1N?7>;#6kMTY;D3Ew^qD^B>bW$R-k(U;t; z{20YKP3^H-OqfDWzn|yPda%-=arg!f2ole_F!HGZzUTMVe~FlLeBs6{R+#P=!r(??ZO-X0q{saYt|iWOWg0Gj!E~x%6c!$>0jq_?P#B;v6&ZPqrkpC`f}!SC=GNv zUu=;u-(|;VZVAM^eWQ^RsJW^ic6JUe-%XRizTTXXc5C|^_;MHZy>x%pi(!$no&_a( zU^x||VP`gq#DZFM&W~gRbB;uHX`ul;&-6Fy?g5BQ$kh|kuAle1c*T>EeQsZsk_rN~ zk@p=-kIvM5UVfNcrTI}(4v~}iNtey+-;bf1^b1<|FGkmY1DT2_Z6x82V+k`uu*-AW zGL>S-xN4MWZ>a%XBT-r8N21cMZvDd?6a3-Xw1zuaG5p;p#ryuOQUTk~IFl@MgCa1x zMt5~n%nJ=WLzv-oqKP!s?VmDQR4Zw_up<5x=xlnnRBC>pi_4rRC$^T_nD;UZWrb!! zr$0^s2NoKshGXA6g9!6bL}jT7mS1Di98)GN&nH`oIktN9b4E%MG1y^x<{9tB?i$blC4;BE|oqP1L$M>pKIKS zpi#Oq>W}SFtJwTC_RbIq_0{#0Lb4lVhbiY46W z2D3rhI)sEWaujg$&GN8I>yQFO#6YN;|6sHxt{_ArqV>i zZ{0$bur^zF^ytaco-e-?R;CNsEU}{HutyxFiDBjbH3>h&7_-HaIjh=hXi5`b*~|uL z6pManm|4ZXG_`J}u`Tmp5SK2XK;WgP_vcnBGjCuOEG7Pv$H|iBojrvdm6bK{z(cH` z*;B8L7-QXNH=u&HSTx5zY}0dMdm3l=Ny3ro25+T5rw+@F*p0RqjX&!;Nw&)@96Z;7 z>CzbU-iHEpj4=Qu6R9B|Vwu}ujxDh{K9_dxAF4`#RewN}!lg=Ni-3OPj4<-FP}4&< zi7*7wI*;&Li&?iRQ~-*j3z}~u-c-9)+LpFqSs|M;#fdj{%)D_lfri>scJElC{2el*L5rsec&BB;k=-PVr@X+`7eHU-7AAL^s^% z=G3_iCZX6Dp(^QSyq~Q0)ZqH;0(Q{9SqPCz2YojN;FHLtP*jaCmw@&D?lbp>SuPOJ zC!is#2yaie!nUJDbji9#*%W$c0c0A^zBe1)18~tx5JWOdUMe?0U8w+I_QUvC0=v|Bbr$+jQ>^D z3(kXU$)zi~7d#iSmy$L&YV&f$QHFgAI>S4adn4 zUyTjRJ`tWyZBiD%&?}+^-%_t2zU;2gUs_G1&5^ZL0&wWjlja`WeA#Fyk$T(>OFfg6u&acRX?0oZDuMk)!8-C?aK9QB4^o(*VA2<>;0bHAb8UdXo z#H)`y)%+Wxb|CZTooAL-IB!nZikgBc0w$oOvh$Or(;=-+P3*)h#_R1eV2QHZ)RfhP zn=1Blq3;C>aQwF7cQraVIc893^b|e$M;yckqb7<SPmGjFS7Yc!Lz!AKHPs-yS3TC z_Pl0=E3-eljK`(d-LdecS%ck|ltYtCyE6eTF19{zV!<*I?U$iD2@Q@{?L^waT(NfU z=+xHA(?{L@tTI@bZ5XW>Pq=cbG)vu2`nh(JlT3mK z#}G1+G6UbZ@)F=N1GXC9wkY~*s-%)MA`+;sc11XZ2-~nDKvi$#x+`Ibv<-ZMaA;h- z>L@Sp#=k|05;}H*SJ8E&^_ddGaiM&=vvAKNsbtO-1f&5`SD4y2{|(sV#n}l?SXiYY zDE?x{a+l`ZlcYG*ld9~axo#ffZ;NhU03*B*#YEZ$K!(oO`_pf-;S0a@GqxfU07*(8 z?|++mGP%--03iQAb_?k-0Ij*^5u=?bJgOhK1Uh&B?hK-8=qVdXZc}bS_5Y~;GmZq4 zK>~d-m|1)R9&n|_e1N;?9QRK_b&}AaByWOmDUZ|yBA^4TsV8o)r>;*u!2AcHZGn4!N}^H6NdoXmPE!%!|CLA)~+UlE6_P+k#w>G7k&A{jju z=?_-|+Vx9rDzZ+Q%URaNEj^|UPB3ju2uPwUI=2#+gkznm0>K6tc+yM<{!e4r*a2Aa ze|-XRX>T=XjnPHp=equnGKbZLBmX4WjRoJ+M)l@>b7zA}ibTB~0RCR^e33&oy<+J%N)3#w_NG0n&y1#)?^{`#aLFz;Ua)qi@KKZ*uW5XZxV!e2Ygo_N^zT4EP*949az zcUc0%^hHTR0*8J>YwhUSSIfqI50cud@|Gfz8XJ-kU-yhnH%Sg_9d4?&;X4cAEKn{r zLwp^DZFfe>7o)fwDjHMrGnwlIXmwht;e&#C9(Arv;_B9(uOZ1}xVS>Glp{N{ql)9p zL#;X*G>+@x;Y)jdhHa*kt?&Dg!qUTWJ~l));=pO8PY>ey*}eMRgNkFsN%c_;Gu|1u z#=R6M)6Sb(?I9G-@(1N%H~g^oVIA3vSY)2@S?gU(>ySd%Gn>4GW20HV*A z0q!N|zE2PR`+XK}Ib8%Wj_0g7c5Tn+!I{8|AMoJ7xaEkTT4Yn3EC^h(|BA?}fnwfcYG{maP8PPNL);m14t zn3Oy{IZ)B;=%Ea^Bsbht+jy>iLnTNqUU-Rkvq^l^*XIzt8lMewxF%Y>*uNEwHsTyhkr3$EkAks1P-g}9P!UYkLBE4gwilGSv5di@S zf`B9xLAnqr(gI{pxbK$Ve%uec3n}N!ndf%E2oecFq9Mp62pYVwfC}bF2%1W!{((Uf z5)EEhWD)^f09esS%jmmN=%gfc5*7{USS%U~J}F5^ND{aJjZ8|yg8w8eA&G!ZB7jc{ z77f9ok$@YU1bD#;7E2&t0i8m?QYeHZBms*);-%6E6bc0h7v({TiUL^+* z773erw}BE&UfSlq5D05K^-a_1S>^zN@bc(rqKt#GRc|ftg7xM;FeEJkbUs_G5W&4iS#ix%X9%x7DLr0H2y<~Ij_y>kj z=>6X`SjefcC~-AL!T%?uvh0yaa!q9T$I{A!k%=&C8JlBeS=JF`6zmTTHtO^Dyh{rf z4tPK{`#2qCh)sxGguH9Lw*fK0S~1o=w9*=@)IwBP zSlxUD4}rPSEEilcJzQYT*N5NU z0yaMGpZf&n99t1_KYF8HJi&_2t_S+xv9omCVgiw#^}* z(Y*1U&8x6xIyaiv(?MZw&JUAYSHG;y*yEp_AdAPcG%6(j;mf< zi_Au^P3}KyIxdXKp|Mb7f=?K6o{JLCqjSLTLcqK$Zv%M3UrDHWS~nW?9f(ELaTOeg z(0xP|J>!r6T@qL6iPMjQ;URbFgZ#?h)cY0Ou@P7p4In@pW*?APzCI3#4xXCeS2e^E z6UA=uKQl0f_tI{AnEV+uRc(6uEKDDHD=Ik8K9nea-QMUDOF9FVPoB&@ zc~)2=yNn~9;pW%JH=~1vgtITOr!!d9d=s19XTAx*sT0`t2nBnfU3tztB1l_ z-Gam4ztcZ&dX2J8S;uG1*Y&kp&a|e*629aSwthxVMntqUkysm=Tw_ny_Dtatqbo{^ zPAq>gCrvMZinL51tj=|FZ*t+oC&G868zf{d13i9r6r)t@1RHIgHwK!$xC)W20@i#1 z7u1{JtVKm9g1v*|NW8VX1#$~@hAW7`*oMvd)tVtYdwlWWo4>){ZNir3cp*SwI2H0 zmt>E)k43QRU+kyv(Ezp&B=I8I8U=Xyh{1puCv%zp8ZIO%;Ic+1gxWC@~ zJwwA26v13+Cb%E|$@;dKvSatOa#uJ<+)$yT}poGwckeAChX$sqxUX!yG5&>9gZ??liQu;lcV zsO0c88j)%!Hp{A0DC5DXOdgEWdQpfH@;+;QQj9C2m%X<9Y8LL)H+D?&43*1F3#^Yu z?0_{d*82e*eD!L`?8Ux}S=0F{6v6*mW5+d#7Va<6_%AWEgM%GP>$c7VZX^xKPNk01 z?XZKeDtlM_gu-{ucukdc+W@Yb6lA2vGC8_~R0|mPYD5gl69(YKz zJ)RQ&zAprhy!YdwbL_xAdTsZhPYkvCcIPoupNw=j$@3vSHH5}eWPYA>vk}nN3DidF znYPhUqo=f1W|m!yN*vKsvOu>1GmG<11vr@Or}*9IDiEgGP>Z!tkDM9qU-8sR&XM3S zm%T0DsjM8bQ>Uxxv^0flws^6BXr>%aT(Q+{q;nWh`d%9l)Ur&{(3Wzl6AHG7K?c)S z8#@f}1%;A|Wfo=&1x*A{1PP7>ji83c&~U)iKZ0#7{lG=FJ>HuK6A(J1e?2Wu?DK5n z;k*rgYr~t@{2}V)+noiE+TeEmzj^)0KYE;b&ErtdzV_LsL}4ztZ@>`olkvI58kML} zLYDyiMQC-G{(1_{%eXj#l257$^SJxjZR3vbaf zvTHrfmd=30c82~oPe3&#IxEMYv=JCxU&r@NYbE9mo#C!_Wn-c4@bn%!Kw>%yq|@O(nwlA37CjiMN#)IEP~ECh`Twh3u`5HIhe(l_trs)Vo+NYC=vt*OH> zn#jzqab7B;$uGyy-Wt5b9P9l{S~8ZSi7qGxc3++v*;+5%JPkHpmG#?aAF{D9x(o89 zzesdA*_9{lWQP3uQeFd-`quEP-#;rkiL&P1p76%~tgH6;jKS-tTg~^^oj6oI8=ZfK zU9Al`&KV?DxS?{~af~Dlu-e}symcG9!3a)E9W7EvG52wEN|eQA^iOlvJ13ItGlB zljv`_zb+tkfoH;Kc0C@0Q&Ei>B zWAnX&SYdJ8O~L(m$hfBptXQ#D7xRK$ERZagDEWpfV^#AS?4F`AOc{gUD#gt#)y{D) z`Ep~BU7VV(d>9@9R8!{fEzYHl^bSGXwDjK4<#+`?Oryraa0vVbF946MZ?#m$UhWso z%!@=E0FX<5Q-MJP?E6s@({d=WD=SEpJOPjc_t3FXC|wC%BCp@ZkDkNbFSl;~MU>>M zqBlkWB2Vp~`8q>RixA;0J%9$^VJ=47!`yRkkJTz<`4GB9feDORsTXfBCxsW#YMA@)1>WqVTa9^)n;eOvz+F$0jCb^X!w;~#OmI@M8`$^w0Aprx*TCn z-(!1V7r?JCQeik6BxxuX{u z;sg(Y=J)kSXBt}9)sg2FC%-CGJ|U$Cl}&~@uJNbez9!f={Y0qt&qe&W(`3Uer?T81 z3qB6SxUB>Gz*K|`)Qv_iWylr-uW&HV~WnSPu_?|7u7?gpOCk#8{CqbvLB~?Y{E@+i$kVe0hy?>Rl+JyvgGnm2y9Cu|~Pp0eAfru9VRbw=;!6{5r^`tmBbTNe_#ks-$B z2p7|;{k<~BWEgIUP8qZ|jP8-lS+2>oKA1cIvPUVd;hc^b<_?RvlyV1C0?JURcK3=` zh0b5g8E0R2bBiVO;l&Oa)h%+NhM9dClR9?F&m8DSp}I#kHW-KtrD|zOh5Aaz1!>tWlk!F(&yLxrQVEHYd(1vP?IS-713BJYzOThsNEXb}!`0&|awRC#d1Y(SAP zuC8qaHDXm2##AC}+P2P31$CQ_SG@le*1(aA2g#^@q1tdS8y6(^rsi|lJ8lsBL{~Tu z>|AZRH{3-b9+5vTIp8>*?>E~!&7c0Y{*(E>RQg{#dreP|8CdXDJYSXa835k$2YB;o zTng?IW~*Y>Q!7pjL$(AR`~@0)CO#P_lE*>?)~87Y5Fq+YLu{M;ptFsT62{BGl=;+M z1u#f350Mh~8p&-mDR*A6*j>o-&UAK&06Sg4480;5xD&}7MHh%>n>>AZ9H?8SN2Q3^a8M) zDm^d`DG(H^)M(8T#kWi0$rJQ!)h0YczpH2Z$@8R^TFO`MrN9owzGj473`qSi`lc<# z$9uL7Ww=s0q$9Z=F;#eTsKTdd#|u|8L&>+oKwY$XE_BMvcNx9@?A)5%Ar`2mOI525 z#6#7!{~R6X-&vRH_vmB?UdSG^gJX%xd`a+!ZCjiK>I-T!6B{hDpRs~qO$82X%5am` zTJYWAT_{-O)n*=Ilcm}QuE~I1j&X*+fR93l=?iG(;&Va5G6k@5aPaxN{prg=s%Hrn zGrk2{Di6)yTfgV64)eYNi)0J=Y0stJ8ZXedX{yR-0p8d;o?eV)?93ofokjh!Pj_q~>zs z8{f_E{q=rVK`MR|O@&pj;;x2F?G=lpS@$-Z6}i$pfvt_V&U5zo zqN2iY_s};y0=uV^Uy0V~Y2Ws5TNf52)H$MbM9f&V&tr0SM!%#u7%@-_{JGLI*)AKh zWxTNxBl&Qi6e^zdhsV(?4srR-`916M0An^ztsmu(q zfG|)C^*7v`(p|pQ4g_=qx9Mk6i3sVgvX-)r42_VxwbgOjq_QfkJKb8!KGH9C&dQ@S zZtuiQ1ov3D-7ss~8FI`$?#=rnl4gzV50X$?W|=6*46%Dsu)c{r9!ySfr(}A*)}ffn z_N6lPn)1Y&UWP)dPi{v%QMxlRx>)_L)s(u~@Vg7xwmsaw-3x^?t^5~s{=-@ono{-6x(rUwSCf37T^X98> zi%t|nfk4O`B2oyq8`Y1hz}XEedpzSA8-e9aKi=Aq)h4CaS z-n|<55Sq+8R;-(vM3kas>E^E}P}HF5W?8u@tf%df9{unq>+ON;{HUJ~170QzpH=8p z)9{IXz3|0I%}1ljxcjeUUUc56hISJ_#|uIY`IfI7n7UbH@JYlUB8@U9lZbc;lbZWK zA96csuX?eMb2I;Uaqd6*pC842Y{_kPy)8tbQ84^bJ|?G2I>BQ#D4ycOo3T{>s|rh7 zbL$VHi&a-ngg7!*8oDki1nJG1Z12%%(wwN+RbBl+$*p{H{nzk8$Z1~I-}@{1{`#-# zYTCzGy9|aS#T74nEAz5;&0Sq=jed4dNWR@K;kE84ZcO3- zhLu#bq&)wHGd}Zy_6}m{=n3X;$<7qdaiMcQ7iNafUv~*=>%mxD^Er{ZbTdGoaO`{c zUgbT85+1h0h-KUH!AP1@Hp9OnlEL9E>zMdN?zIkeSaXW~k5Gq2q7WU=DsNkP>i6}N}>a@Ra3=If0gKP#%MIln~eVby;w zufy}rvit1CJ-ZI?{Qbq-R~iF@I72&t5KaWq{}5u;Q@Nt)caJ`K=KHC2AJWFSH|62t z$YX*+U55L|;*c|nZ&spyIUJ1CwBJl%+J3jSD6y}8;Jj0l&}fTY;aKn&pXOn$c9srL z{Px>4z=&|_pBoqRu&+9=d{aK){H>CHrS?*@1`Wj5h3m-%GP|GhR)+ajV$(>J-O%hG7*ODxHzj1UD