From 9632c00e92297d6cd93ea325794be9c7a7538b2d Mon Sep 17 00:00:00 2001 From: "kathryn.doering" Date: Mon, 6 Jul 2020 16:16:56 -0700 Subject: [PATCH] fix, test: #22 get readme example to run with new iter_vec input --- R/runSSMSE.R | 22 ++++---- R/utils.R | 6 +-- README.Rmd | 10 ++-- README.md | 81 +++++++++++++++--------------- man/figures/README-plot_SSB-1.png | Bin 9540 -> 9534 bytes 5 files changed, 61 insertions(+), 58 deletions(-) diff --git a/R/runSSMSE.R b/R/runSSMSE.R index fdcac77d..ae91ecfa 100644 --- a/R/runSSMSE.R +++ b/R/runSSMSE.R @@ -151,16 +151,16 @@ run_SSMSE <- function(scen_list = NULL, #Now set the global, scenario, and iteration seeds that will be used as needed seed<-set_MSE_seeds(seed = seed, iter_vec = iter_vec) - # Get directory of base OM files for each scenario as they may be different - rec_stddev<-rep(0,length(scen_list$scen_name_vec)) - n_impl_error_groups <- rep(0,length(scen_list$scen_name_vec)) - rec_autoCorr<-list() - for(i in 1:length(scen_list$scen_name_vec)){ - if(!is.null(scen_list$OM_in_dir_vec)){ - OM_dir <- locate_in_dirs(scen_list$OM_name_vec[i], scen_list$OM_in_dir_vec[i]) + rec_stddev<-rep(0,length(scen_list)) + n_impl_error_groups <- rep(0,length(scen_list)) + rec_autoCorr<-vector(mode = "list", length = length(scen_list)) + for(i in 1:length(scen_list)){ + tmp_scen_list <- scen_list[[i]] + if(is.null(tmp_scen_list[["OM_in_dir"]])){ + OM_dir <- locate_in_dirs( OM_name = tmp_scen_list[["OM_name"]]) }else{ - OM_dir <- locate_in_dirs(scen_list$OM_name_vec[i], scen_list$OM_in_dir_vec) + OM_dir <- locate_in_dirs(OM_in_dir = tmp_scen_list[["OM_in_dir"]]) } # Read in starter file start <- r4ss::SS_readstarter(file.path(OM_dir, "starter.ss"), @@ -188,11 +188,11 @@ run_SSMSE <- function(scen_list = NULL, } } - rec_dev_list <- build_rec_devs(yrs = scen_list$nyrs_vec, scen_list$nyrs_assess_vec, scope, rec_dev_pattern, rec_dev_pars, rec_stddev, length(scen_list$scen_name_vec), scen_list$iter_vec, rec_autoCorr, seed) + rec_dev_list <- build_rec_devs(yrs = nyrs_vec, nyrs_assess_vec, scope, rec_dev_pattern, rec_dev_pars, rec_stddev, length(scen_list), iter_vec, rec_autoCorr, seed) - impl_error <- build_impl_error(scen_list$nyrs_vec, scen_list$nyrs_assess_vec, n_impl_error_groups, scope, impl_error_pattern, impl_error_pars, length(scen_list$scen_name_vec), scen_list$iter_vec, seed) + impl_error <- build_impl_error(nyrs_vec, nyrs_assess_vec, n_impl_error_groups, scope, impl_error_pattern, impl_error_pars, length(scen_list), iter_vec, seed) # pass each scenario to run for (i in seq_along(scen_list)) { @@ -341,7 +341,7 @@ run_SSMSE_scen <- function(scen_name = "scen_1", } for (i in seq_len(iter)) { # TODO: make work in parallel. - iter_seed <- as.vector(mode = "list", length = 3) + iter_seed <- vector(mode = "list", length = 3) names(iter_seed) <- c("global", "scenario", "iter") iter_seed$global <- scen_seed$global iter_seed$scenario <- scen_seed$scenario diff --git a/R/utils.R b/R/utils.R index 4a8c1375..9944dd58 100644 --- a/R/utils.R +++ b/R/utils.R @@ -473,11 +473,11 @@ create_out_dirs <- function(out_dir, niter, OM_name, OM_in_dir, MS = "not_EM", #' Locate the OM model files #' -#' @param OM_name Name of OM model. +#' @param OM_name Name of OM model.Defaults to NULL #' @param OM_in_dir Relative or absolute path to the operating model. NULL if -#' using SSMSE package model. +#' using SSMSE package model. Defaults to NULL. #' @return A list with on comonent, OM_in_dir, which contains the model location -locate_in_dirs <- function(OM_name, OM_in_dir) { +locate_in_dirs <- function(OM_name = NULL, OM_in_dir = NULL) { # checks if (!is.null(OM_name)) assertive.types::assert_is_a_string(OM_name) if (is.null(OM_name) & is.null(OM_in_dir)) { diff --git a/README.Rmd b/README.Rmd index cf39a83a..c3b5b505 100644 --- a/README.Rmd +++ b/README.Rmd @@ -57,10 +57,11 @@ You can read the help files with ## An SSMSE toy example Suppose we want to look at 2 scenarios, one where Steepness (H) is specified correctly and one where it is specified incorrectly in an estimation model (EM): + 1. **H-ctl**: Cod operating model (H = 0.65) with correctly specified cod model EM (fixed H = 0.65) 2. **H-1**: Cod operating model (OM; H = 1) with misspecified cod model EM (fixed H = 0.65) -Note that this is a toy example and not a true MSE, so there is not significantly different structure between the cod EM and OM, except for different steepness. We will assume we want to run the MSE loop for 6 years, with a stock assessment occuring every 3 years. The cod model's last year is 100, so the OM is initially conditioned through year 100. Then, after conditioning the operating model through year 100, assessments will occur in years 100, 103, and 106. Note that the assessment run in year 106 will generate future catch for years 107, 108, and 109, but these are not feed back into the operating model because the MSE loop is specified to only run through year 106). +Note that this is a toy example and not a true MSE, so the OM and EM structures for both scenarios are identical, except for different steepness between the OM and EM in scenario 2. We will assume we want to run the MSE loop for 6 years, with a stock assessment occuring every 3 years. The cod model's last year is 100, so the OM is initially conditioned through year 100. Then, after conditioning the operating model through year 100, assessments will occur in years 100, 103, and 106. Note that the assessment run in year 106 will generate future catch for years 107, 108, and 109, but the future catch values are not input into the operating model because the MSE loop is specified to only run through year 106). First, we will load the `SSMSE` package and create a folder in which to run the example: ```{r, echo=FALSE, results="hide"} @@ -113,7 +114,7 @@ EM_datfile <- system.file("extdata", "models", "cod", "ss3.dat", package = "SSMS sample_struct <- create_sample_struct(dat = EM_datfile, nyrs = 6) # note warning sample_struct ``` -This sample_structure suggest that catch will be added to the estimation model every year (years 101 to 106), but an index of abundance (i.e., CPUE) and age composition (i.e., agecomp) will only be added in year 105. We can modify this sampling strategy however we would like. +The sample structure specifies that catch will be added to the estimation model every year (years 101 to 106), but an index of abundance (i.e., CPUE) and age composition (i.e., agecomp) will only be added in year 105. The user could modify this sampling strategy (for example, maybe age composition should also be sampled from FltSvy 2 in Yr 102; the user could add another line to the dataframe in `sample_struct$agecomp`). Note that length comp (lencomp) includes an `NA` value for year. This is because no consistent pattern was identified, so the user must define their own input. @@ -156,10 +157,13 @@ The function `SSMSE_summary_all` can be used to summarize the model results in a summary <- SSMSE_summary_all(normalizePath(run_res_path)) ``` Plotting and data manipulation can then be done with these summaries. For example, SSB over time by model can be plotted. The models include the Operating Model (cod_OM), Estimation model (EM) for the historical period of years 0-100 (cod_EM_init), the EM run with last year of data in year 103 (cod_EM_103), and the EM run with last year of data in 106 (cod_EM_106). -```{r plot_SSB} + +```{r results="hide"} library(ggplot2) # use install.packages("ggplot2") to install package if needed library(tidyr) # use install.packages("tidyr") to install package if needed library(dplyr) +``` +```{r plot_SSB} summary$ts <- tidyr::separate(summary$ts, col = model_run, into = c(NA, "model_type"), diff --git a/README.md b/README.md index a378f5db..1c6d0490 100644 --- a/README.md +++ b/README.md @@ -63,21 +63,24 @@ You can read the help files with Suppose we want to look at 2 scenarios, one where Steepness (H) is specified correctly and one where it is specified incorrectly in an -estimation model (EM): 1. **H-ctl**: Cod operating model (H = 0.65) with -correctly specified cod model EM (fixed H = 0.65) 2. **H-1**: Cod -operating model (OM; H = 1) with misspecified cod model EM (fixed H = -0.65) - -Note that this is a toy example and not a true MSE, so there is not -significantly different structure between the cod EM and OM, except for -different steepness. We will assume we want to run the MSE loop for 6 -years, with a stock assessment occuring every 3 years. The cod model’s -last year is 100, so the OM is initially conditioned through year 100. -Then, after conditioning the operating model through year 100, -assessments will occur in years 100, 103, and 106. Note that the -assessment run in year 106 will generate future catch for years 107, -108, and 109, but these are not feed back into the operating model -because the MSE loop is specified to only run through year 106). +estimation model (EM): + +1. **H-ctl**: Cod operating model (H = 0.65) with correctly specified + cod model EM (fixed H = 0.65) +2. **H-1**: Cod operating model (OM; H = 1) with misspecified cod model + EM (fixed H = 0.65) + +Note that this is a toy example and not a true MSE, so the OM and EM +structures for both scenarios are identical, except for different +steepness between the OM and EM in scenario 2. We will assume we want to +run the MSE loop for 6 years, with a stock assessment occuring every 3 +years. The cod model’s last year is 100, so the OM is initially +conditioned through year 100. Then, after conditioning the operating +model through year 100, assessments will occur in years 100, 103, and +106. Note that the assessment run in year 106 will generate future catch +for years 107, 108, and 109, but the future catch values are not input +into the operating model because the MSE loop is specified to only run +through year 106). First, we will load the `SSMSE` package and create a folder in which to run the example: @@ -93,7 +96,6 @@ library(r4ss) #install using remotes::install_github("r4ss/r4ss@development) # Create a folder for the output in the working directory. run_SSMSE_dir <- file.path("run_SSMSE-ex") dir.create(run_SSMSE_dir) -## Warning in dir.create(run_SSMSE_dir): 'run_SSMSE-ex' already exists ``` The cod model with H = 0.65 is included as external package data. @@ -106,10 +108,7 @@ cod_mod_path <- system.file("extdata", "models", "cod", package = "SSMSE") file.copy(from = cod_mod_path, to = run_SSMSE_dir, recursive = TRUE) ## [1] TRUE file.rename(from = file.path(run_SSMSE_dir, "cod"), to = file.path(run_SSMSE_dir, "cod-1")) -## Warning in file.rename(from = file.path(run_SSMSE_dir, "cod"), to = -## file.path(run_SSMSE_dir, : cannot rename file 'run_SSMSE-ex/cod' to 'run_SSMSE- -## ex/cod-1', reason 'Access is denied' -## [1] FALSE +## [1] TRUE cod_1_path <- file.path(run_SSMSE_dir, "cod-1") # make model read initial values from control file and not ss.par start <- r4ss::SS_readstarter(file = file.path(cod_1_path, "starter.ss"), verbose = FALSE) @@ -122,14 +121,14 @@ r4ss::SS_changepars(dir = cod_1_path, ctlfile = "control.ss_new", ## [1] "SR_BH_steep" ## These are the ctl file lines as they currently exist: ## LO HI INIT PRIOR PR_SD PR_type PHASE env_var&link dev_link dev_minyr -## 107 0.2 1 1 0.7 0.05 0 -4 0 0 0 +## 107 0.2 1 0.65 0.7 0.05 0 -4 0 0 0 ## dev_maxyr dev_PH Block Block_Fxn Label Linenum ## 107 0 0 0 0 SR_BH_steep 107 ## line numbers in control file (n=1): ## 107 ## wrote new file to control_modified.ss with the following changes: ## oldvals newvals oldphase newphase oldlos newlos oldhis newhis oldprior -## 1 1 1 -4 -4 0.2 0.2 1 1 0.7 +## 1 0.65 1 -4 -4 0.2 0.2 1 1 0.7 ## newprior oldprsd newprsd oldprtype newprtype comment ## 1 0.7 0.05 0.05 0 0 # SR_BH_steep # remove files with M = 0.2 @@ -188,11 +187,13 @@ sample_struct ## 1 105 1 2 0 0 1 -1 -1 500 ``` -This sample\_structure suggest that catch will be added to the +The sample structure specifies that catch will be added to the estimation model every year (years 101 to 106), but an index of abundance (i.e., CPUE) and age composition (i.e., agecomp) will only be -added in year 105. We can modify this sampling strategy however we would -like. +added in year 105. The user could modify this sampling strategy (for +example, maybe age composition should also be sampled from FltSvy 2 in +Yr 102; the user could add another line to the dataframe in +`sample_struct$agecomp`). Note that length comp (lencomp) includes an `NA` value for year. This is because no consistent pattern was identified, so the user must define @@ -217,7 +218,7 @@ dir.create(run_res_path) # run 1 iteration and 1 scenario of SSMSE using an EM. run_SSMSE(scen_name_vec = c("H-ctl", "H-1"), # name of the scenario out_dir_scen_vec = run_res_path, # directory in which to run the scenario - iter_list = list(1:5, 1:5), # run with 5 iterations each + iter_vec = c(5,5), # run with 5 iterations each OM_name_vec = NULL, # specify directories instead OM_in_dir_vec = c(cod_mod_path, normalizePath(cod_1_path)), # OM files EM_name_vec = c("cod", "cod"), # cod is included in package data @@ -246,12 +247,8 @@ should be installed automatically when SSMSE is downloaded. # Summarize 1 iteration of output summary <- SSMSE_summary_all(normalizePath(run_res_path)) ## Extracting results from 2 scenarios -## Warning in ss3sim::get_results_all(directory = dir, user_scenarios = scenarios): -## ss3sim_scalar.csv already exists and overwrite_files = FALSE, so a new file was -## not written -## Warning in ss3sim::get_results_all(directory = dir, user_scenarios = scenarios): -## ss3sim_ts.csv already exists and overwrite_files = FALSE, so a new file was not -## written +## Starting H-1 with 5 iterations +## Starting H-ctl with 5 iterations ``` Plotting and data manipulation can then be done with these summaries. @@ -263,7 +260,6 @@ in year 103 (cod\_EM\_103), and the EM run with last year of data in 106 ``` r library(ggplot2) # use install.packages("ggplot2") to install package if needed -## Warning: package 'ggplot2' was built under R version 4.0.2 library(tidyr) # use install.packages("tidyr") to install package if needed ## ## Attaching package: 'tidyr' @@ -283,6 +279,9 @@ library(dplyr) ## The following objects are masked from 'package:base': ## ## intersect, setdiff, setequal, union +``` + +``` r summary$ts <- tidyr::separate(summary$ts, col = model_run, into = c(NA, "model_type"), @@ -304,14 +303,14 @@ summary$scalar %>% # plot SSB by year and model run - need to correct using code from the # think tank ggplot2::ggplot(data = subset(summary$ts, model_run %in% c("cod_OM", "cod-1_OM", "cod_EM_106")), - aes(x = year, y = SpawnBio)) + - geom_vline(xintercept = 100, color = "gray") + - geom_line(aes(linetype = as.character(iteration), color = model_type))+ - scale_color_manual(values = c("#D65F00", "black")) + - scale_linetype_manual(values = rep("solid", 50))+ - guides(linetype = FALSE)+ - facet_wrap(~scenario)+ - theme_classic() + ggplot2::aes(x = year, y = SpawnBio)) + + ggplot2::geom_vline(xintercept = 100, color = "gray") + + ggplot2::geom_line(ggplot2::aes(linetype = as.character(iteration), color = model_type))+ + ggplot2::scale_color_manual(values = c("#D65F00", "black")) + + ggplot2::scale_linetype_manual(values = rep("solid", 50)) + + ggplot2::guides(linetype = FALSE) + + ggplot2::facet_wrap(. ~ scenario) + + ggplot2::theme_classic() ``` ![](man/figures/README-plot_SSB-1.png) diff --git a/man/figures/README-plot_SSB-1.png b/man/figures/README-plot_SSB-1.png index 351dcd1e69e40098d46420cb6b496ea6ffebe917..bcdf77f1fae27eeaf1aba5c427fdaab431b37ed6 100644 GIT binary patch literal 9534 zcmcI~c|26_+yBT6Vc*&+GZ)`Rn)lo!89kocr9@`?{~|zVCD2*SX&3254r$h7-V{P$)LT zQzy+)s68+W#qb%+h}FFbN=_9-4UZ&-UB8qAG*RSPYcWc~HHSWZk8uyyvl$zn) zn%&`=-Q7AL{<@L$y4~Hz#>Tf4)3>|3{Z6Ri6yk7-J0dc?yE`(1XxjZ+?)7zd7g2!# zcHdHVhf{u|hlh7Z(vdo(0@1X)i=cla5OGA7%i-rR(q}AJPT5^Uq4vG{{by+Lrk+Qk zgi(ek0n5Ow`E1tQ5mm?7Y)}6@Q%2`|sEk5r(>GV;sM&>hr9VjMv508MEaF zr#2%nGEV0W?}pZ{Z2SrfZ++IWq_?4Epw{kxn1$i6aRXztF95J8otcl@Op7#=uguM3P1iYh!XFu@bl|L~7{f}h{6i)yc_vAP)I1E8=hK1h}^}lZ?@1E5%nzd3+#m*Ez59L0-il-ft1e-N1(XW|L)M!*F^ufpDcvfaKxnFo$_Fe@u2ylHW5xW z0~odSq>h4yB`^b=kLRNnTF_A-DS_4Uy%ZTmB{XBCr}4@t%e1jQcYJXIC}|^~>qbNj z8;YmsskTHm#1xre?(qL!jnf%0&-UDb?(6?*Z+ry5EZXqBUU8D?-}ZyEq;nuI7HXk^ zN0e6=_&=>&v2FO(!H@gbKFx>lfthaIe;rJgrhWho1w@RfL641rL~`$LEeg{tTnzh3 z)v@#^nD@NB?zu2$1(KRRppWRak|#jA$OB3070Oi{&cj9&?~A>}GSp8xwZ=`h##^#q zeDLJP`kXe{{qz0uo=mI!^E;<};yzOBz%SNwIfwZ?94=e)SXRH4niB3`6rZlkfih)( zyzJ(wd5lu$rdKVk*W=ux?) zf9X%XBMiDOUEhOOX7m@O?#;?Pui^ZrmABmK z{talK%f#|n`DAe^Rep!zKwTqj$CSA>kKNd5)m zVp9gkuc4dw-pWz{oou*vUU4(=&m8NqMu+)1SxvRQ`q+`k-MpRnJm#Jr&u(i97p z%<+4B6fo8m(G^d>VQvtF`1W4wU0N+5=Nq@ZeLG*T4A-6SJk+!=Vy^(8;Cd1uP4OvL zj{9A*NRFB`HJ6>o152~N^zy~)<%|ha!ucbyZ`Qtavkv*JJuxw-ixl$UpK{h&rbVyz zY6Eec^Ct{Hb@_{t)o$)C?&q|LO zfiS6eGX}4FE$EX~XR$MRZxtDajEU|*mB1`(_IBWFj@Bj8)t-;r6Jxiaedj8;df)r9 zLq?JvBcn1xQd6k@im};vzrLk-kyNF8W$o^n19)E)=MbFdf_aY3<}7m5E8mj_?Gt`K zLJ3n)w+A+_79Q0`aAVjq52eutZHg7WQ-^l(9xHChSQ*d+$`ioX`#X}&ZAr3J7~VdK zr}5AVxGXd9xgT!yUF^CkK_JVKFn&C;vw?^lo@R3UUcX6G5nK>2PUgm1bh~m{x)O#8 zRc^oe(lMEcr=?Fp_6e%XGhZ?spqbj4Qj<&WPcQPqd#=^I{e)4$Yl;GA$a)o5d!@Bw z-f>to83WmM5B6mqa|AS*TNTTHf)OQfm^c^h_`$H(P=fUj#kH=KBoOBzGm+hq_zxGK z7xiJ(%p$Wn9Vq0cl`T*n@{OK*JN;~}Dl%L4$!VZrxvi>zmzMKrLi0$e%o!>2Aw!HT zUNacjE^#E$oyB=ArEvR2pUvIT3)LfF1of{16AG*z0v?i8_3jwVyhINCLKH85@~Sg1 z9E*$@NJp{iwaQ!4PL#mgNTL1r17^kmfdA~5JxX;YYx3NMesYrN+A530L4YZY| zhU;e|GTF^bI^yX`=6)WD*AnF<2ye+J^9H#w0 zPWFNf>YX$vKsSkw`RVLqeIxaFZ3?sxLS~($(N9}V%zPYjDP`o&b{b&jumQc5{UxYD zXPI$C6iz@`KaD`vx`kQIxfjhR%5Mn1W0hYx1B!g?xgtf_|D2ch($aST5*3;*f)frO z-Q=-!am>B}-R-?KSPesv|80<7*Fw7nED7Y{X|a$w{X{Cx{teTk)Z2JG1}6{+xujUr zMUUa}dvQY`Paf;s>-QkImy6Z@zpl+$%U)Kw0wT$N<{&o?mmZAs5b>8G2R%?~yQ{Aa z$a3h4OA&B*Wo;RNf&GmI$cjt28;bLI90Gg^HpKTrP`*SiWDkP0um4pC0Ntsl0a*gU z{)ZGb>BD0<#Rgjht44}lPlfWe0mhTuc<-00m}<@xKMhPX3IN|ankxfKGl4h{UEl=* z^&}8%5hxo5V|)Trod=(5ErC_>_y-4&Dt}pOAv}W20tF1FeEhdi5$sGLM+QPKN0u_6 zArQm~#6nE*k^~22-2N>xodl^U9>i}Hx;cwOzzTezPr&#^W^ot+3NXx&uXh9T*Fns} z<2|#Gq1N^vmm)AB9-jwpDMmuX-x{MKCK%q2w91Lmn?-sdZiE%rZ=4VB!>2(^c`(D+ zZynCF!9HQPJF++GiFsVLQFHrHujKOu{C%CMQ}(-3OWF+_q(zZ<=v^_AT~L>mSm4 z>QgKQeQM)Tx5=4^`S{#Gt!Se%p{_PJeQ{!)SuUGv#WZhWdZNOeWmWb<2lkao%JMv6 zYX4F1?0)t`T1P!rxAFo{ROT%IED)o5kQ_gSyUSWtE3y&4)*U;(z7n(fWi{&d^U9+7 zx(i!LQ0Mr5wADf%ySJr6h8AD5^TTr_xpJn+9wyDE(j6YTYYMa8$1n9?op5TszZc`C z`*|jj??>N-&Vf;_AnJ+w@Dwd6(!8IEn;7)aCR*qjY}GK23l7@kOExX%AngWlRh9+?{ zg6?e{O3YeaMrw&W-L(_7<=)4YuV#mRSKMv8xYJOaZT~bQvD87?80<6Lm{G^o?KKdd zQvQ+95b(~ief6k!<=t01-Z%Samg8d5aaC+(qr-T@m6`_V&SSkAp>);el4?Q z;Cl#%?`ng(7vP<4bT)iTU2bAnt6tNN4h1BCV7lX;=7NX89=(3LCFhzVfymglq-3&fi#TCw5;LEFIV$}7j~_p~BGyKn zwnLx5EXm>|5SFohoUq&MwEHdBp_4k+m3jxh5|`?m79DPTCk(qW=hWDsWI)eXHha%h z1RqnAaIygXwqx%#{JJ%6x#M>evv@94u;^5{3oG$jOWDP&k{Nsd<|fYejmfAtjgb|F zFxXHS&-Zh-<#?*{##r?N@$jA_RCysAj2%5~GsZ_I2U+e>93`Q`LPiO41QRjk8J3E$ z58QVb>qC5r;mRxVmvcv_9!zyCKxl(_bXm)JzL*q(EQp=69Kn*}|c4DLCO;qnv zJ8M~N{~a6dPW7K>H{9n!IdOG#-X-)|_kd43*WZ_Blxu9%yQVlFu6NDGTPlPGJW<{7 zo|vNqB%4WQ1(*Qa@&0MneE1^HxUG34^XOCnGTO19=E~ z7?DI5vu(Hg6VWIZTn|>( zKMFDjSuX?Qdv~xbq!tVC&{IhN5^|3nPfL;@)0{V4T=fKVi3c2j+Y9U%Ptrn+2=%5U z;iM`E%E^nZ3IVQy5wEG^6^X{}fu(CYz?ox^swo}ItloYIeHVttVGJEgm}EcVEkFYO zQ}G1ITD==vri(YEe_huG8sLLQPcSN(un$uvDx8Y%W!3H`eZ~q@vDmX_KBd%ROo!1$ z@calt0i9(1DoRzCWUj-vSizf)sUl>b8$_HvBtn)dK#ZJ%4tqyu zLdA%BY6Gk30DRND|LL>SfD*0@;Ghj&H*b3vD0q?Y{e@^BzEe04VJ>EUTjnSG@}Gvd zFEMX(xXPuVO-Io&T*@T1My{UYoxy|)jXV+E>0c-q0>%chKk{F8tb zSyug`QCf_VZM!$6B|xm694Yw(#%=@{nx(VfDZSq zWS}c%&tDSkOKR7Ih|kZ`3}PC)p7~)ao-#<=)2A^VS>Y}MlvALS``xR-!9WfJqqF6T zA;!uxJ1vPGHHw3jP)KnU}XXj*XXxjik*|WBVpIe9X?&u5qAciB(#@UUV2x%SB z)XN^}(@$tx^^~lJS!M%q0xlHzgx0KRw=%7_s;N*FlO!a}AWYEb#+%ZooqaIXL)jf7 z1v4>McKY9F9N)^oteV{Q!I%*`;P{MdtAlU|w7(G82x6D4iVT`}cNSk!OcC z+DuRD%D73Yn)57H00c6JLXr2mJ8>^q_UTVgR}P%{)?C7?9GjO@_xhto*$9-{AWdGwU+atyMU_J zDlHshA?9&5|tJM@7KQ`Bd_>3J)>3o9?5T~xi^~O zu3ZHR@Ok}}cHVQsS3G_tv3jK44t`=@K7FjU>$NP%Zw!*W+dChbVuO*i%GbO#ClAHI{f4d1`(Kycfu7D1~er{~#` zK$PGVcK1>y5{ZWNE98)%nk|}VFTV!yHZnV2S*YnW(^T7FQ(N=QO6RS)l-aSWtuO1k zR_?qp6RT=+ZE@vtjt*h#bsSH<1;_b&)Z}Nxb90{@Z@)HY_SPk%w5h5|x=)kS*w205 zh-_y77IY3mFizsj@@Cn)2Y-W7XMkBvNBv0PqME^jHTqgbpH#8GieScc0@G^$x>d)* z!x}#_+H>`LL66kHrsNtc`!zK=GyIy|0SU)Ib;Xx0)_@40Vj>W$VfDGdXufuJmHoDjz(|;~UN? zaQEFxm7-{7Q6b_2o+;pV|}lU_E{*jqrR4vumiuRE&ajvjj(TVT+ro zC-;XZ0NgsE+qUCqjxNJ%r3w!@kK5gS1j}lz##?L;+_>+6KLtc3fl6JOTRePWM|R8m3N#31V6gc z-#ESdRu`rG;q=<25!O5wVjq~d`z&aN)3TvU95a9TBewcdpj9VM;pOSg5KY$WSn1hI zlc(`C8SVB^F0YMYXcF1ii&wcV!3eS5r78`MV;=p**!rK5!eqPXAiCx!hRFSqPeB0upGNu zCXbeV6Z?9F=xDnwO~-!;244P!wM)@6%@ zTeqWmET8I8g^oxG@s6Hn#R9XIUA(?)Dyi2%hgL0LC=N-u|W&NecBj88QfTUj z1x62z3r>7^OiDkX)~6Cu$PDknTY}S^w4K8?;}A%@9K}VuEZWK2+p7KT+ryTnRHyBaXiK_u?MkLz zbv7A8NQHxiUIqeg#-J-{tv+)DEmtTYWCNCCWhw0)8O~}98bHf|Az5LDDP_wx9vD$^-Ql{w|fEvPPUnwRDkweM`i zUJ%n4r>ba|o|6g*lz6bU;ud50S5g^(m_k^W-ttL6vtUvn4PLz{Jw6lY+;v@=ym!t2 zI0hq_km8rr6N+T?+_FVnXMWj(%d}%>z%J2zVmUX)7AVit606MLxFJctIF&(c^n!4M zips1vZ{Tx5I?(&7ry~SMe27))3Q;TmFkM`_$X9vd_dJoM-Xq=GM&oH)4g2^kZ6Tu` z5LbtGUlkfT4FW}E%I=C=eAsH_WMfv)_R%@>zQMN)Z%e?$PlCX_;53M~QNzn(r%6?%+sV;gmXc>ZzUdINg4q&8>c-ty+W^K|oNOa5qr+A-_8kJy& zJh0ymxYYy|i{DSj6u|cp5>PdN5vQI1vquS6bm-Nk$p0{s689TLX*0XOFnuG4GoEn% z*Uxe~9U{Sc=Q``B=7V(pLszCoaXf@Of2rWj=?H~E z^;(8Myr+Y&;Nmlmv)hM5i{QrV>VNoA57J0UOHYpScGfX{=0C)!?8$4gPG7Rw6~CEV z|54fdvN5=9J?@88=+GqU4@+8Vs>!vu{`8qW{dK{WzjFHQ&9z!NVHGNW^^4WV;%4nH z;eB8D{lmZ`!QRoGpZNgZ&hGw0x~4Ka^V?JTGx!w+{y=_^ScWWDD@v;@E=d1%mEBca zs=VvorL&snx&C@tB&CyBkU7u2S858HjbxNIa-0$Eea2g%;`84!lAtZF>q5cbe%0zA z1v-5rZ?wk#QGyrt9ncBOyLIs|1soYX<3o1*aq=${9F6Tr&h^Va|Hn3*$R3-Ti~t3O zf6M2jlRu+SjK#md1prAnE`rv#3jCk$hGaf&$oyrP<=fM9&;4Ph-KT$+?D-?~Qbwh@ z6aA+^ydvjR`;T73lo&g-3XRvjj1Uo7Y6->;!b zYb1EB^oH5YsZA(o?`ID4Q*_lhEvp))y*O5AE8TZ&s^3D0LC8>)4uX5wFcSaEiO!id z^6pTlcUwU6lQ|^IUtfR6ZBQ^MXwy5P_B@6#h5htZpfSXEV}W#Z21&4whKFZr?S3Ad zsJLIcHW!!&i-Gz`7QI-*s&aSDZD(+zC0WipyZBSpcnm}lsr=N}^n>Uf@IzE3qHgI} zfbuFuP!0(s(rY7@rY!pwnxEI^js`lU!FXRyBuoEEbE9?s!jFz5^S5e7y=pRm4WeXB z<9VraN5{Kn%Y4UA`Ea=oU_kmA*DHTpTCeLzNz?XL8uSp05!Cw>ce;U9JqiJAnayZ^ oKxj;KdTrT``VFT4&G#*>oX`JY9=mo+qmaW;&+OzwopaFt0hi%eVE_OC literal 9540 zcmd6Nc|4Ts-~Y@Q%Mh|>&z3EzY&Eu<2qO#zQ8;BALPxg9awqD$AMC>-t>ReO=e*^SSQgtTBCuG5H`Y51>)n7&qr(#BqCw{7*;YS3fmT*2R)LjPDQ~=gpm=*@ zvGC9P$EU$NEzP@u2}=Z81qNEB1O~$2%D~Duft6E1a5YWLl)sfztWr{buq#vEFw?LI zxSBLrGYt!Cr!}Od(Hqj}%xrVf?56DOiKgrcW@VsNB@C2eRXH_PImK-7Ml^V18ydVD z=xGgfW&@Mnz+^TBi8f7SH!+#5t*viTtllu0qwWZL8kU~s4U5p3%!vsYi8)mnIK^bb z5U>LCO$w8q_CuRaXHI0pJeUC^F`2OT4-HrxmgT>=L5FRI`PJMx7=hqv|M5e1T_t-W z5Qh;KCPswo50~Sq^RR24)=Gk|GXevT`w1kaHDN6^?AwsgyDi6nH57+Ozp0|cIy%u$ z6V<=vwe7BdA8TW{zlEOqpQCQJZC)yuUgE(|tlV})iFI|NwVrNI-;v%AYn1~++uE3W z%zg~fL|ovtih&PXa-c^X2hpwGsExT=!$CYME?@#!X0~j_e((G6R7ASFv+1>}#sA!~ zL$5k{>ca;f-39xHCiPQ~IIc=b+7$v=qeg2$KQ7Vk;{Uw|sY7^MRdAg0PwQFo}W{w>aj^nwx%eJr7lU)%Nk_Xb`QJtGOFUEuejIQ@ zkPE+nt6@Y$pJlCal5ZmkKdd7vD1Y;*nti9WNQH?bqX8ft8=O?;|!s{NfVr@&r_(Y(A zpxMLAEcm*u5tB8c5;kN?-fp z^VlhU1q?G%OSOCgKmLH)Fk7g|Uhe!BXuFm&cuV-D6oSFmRDg}S1SG6*A^51Rh0t%{ zWu?(m-J*O5Rr^8 z2W;oXSrlL%8}&uLVH~IL10_c0ADyne#ajbOSJ%uKzI5MWrtDzp4+6my&tL zj(`$Ul7t>e(oqS2@7pX>G}ZRt2-6c{p286KoEF*fyWM)z;wY@Ka)C$JdQ+0Gq*@BUv^%0PdG)ygsw~t&3gE6Ap?Syic?bWb(Bqq+4`eH;$gXn%PcOk zyK8l*p-&~~(D!H6Nf-Jd?LtTTZ7{`N`l3_ph6t}$e@6G(PTuTDVWoP3d-S$S+x|{B zhqyH&&oaL=?(%u;NnACfwTk49!k*RtwS~D{*sHQXzIW0`1rW+ngQSJTh+R7vp+9&t z3{`NYHu=kZNXZ(h=Wdhgs^9Ll$MxwbmMph5s&;2Md{_K@l3a+Nv+yx9>E3w5Mu_e= zX*fOB-plKhuAl(%avrx7W*l${yVtrCcogC>X1FbN77a7Z=2mE82C(&~?ksp^YA9rv zeA6bbNBnVCjh!@Svym&RFxC{)dEwm=JZ!D@_H(A<->il%#WK$~nD@%7h?ll&t1wJ1 z%YmGKU9riVP(Wd=+Ckyas3Dqvr%}=oUlujwQsx)8-s4Egk_R^1U!&@hl^92(m%kG< zjWBZ@rcph2R+9Qmx?MYV5kiad$$UtL5<#(BXNvI30=Of@jaMhpQIyZA?9`)*&?7-! zyfo_$H*9mmW)$?AiuAB3VsiyO1C zTH*APGvC7P=FQL{=lLqG%X#wzJeF~^UiX)O#<6a5;TM=}OM|SWCRm`(bDq&Im)09z z-UE45T@)>3n|J|*qWS!}TAJODUOW|Xmz{r@($Jw4iv4Zr9BfI`6L)ebuZynS2eCNQ z5Po>FIud~VyIwcIipli9u8ec8$iq^P-+(hkNpqbD2b%%4b+3ghqr{?n%t6iX@ZmYF%HM|B(QM+OMT_!`KuDZkx;Hstw?fHD

t?I^p!deXX8uD|8568Ck{tP@9_=PO^`-wfQ{UZ2ylv^YeLHf zg&H2gKcTL6_#MQurE+K^KNdPwR|- z^k50}@DC0ekWue6LkX_h@gwWQ6(PJVw6Tm9#G!{L^FiagAgx&duf&JH5su4j$ z10)lmYGPB5C-Kxl9I&*ksm1u*8s1uzk!h^zzGvTKao}G7o7r-w)+o6 z>^~4TKOmO&FiII5dh`RrLJx6SbmPnDq^B) z;Sr1{3mrls1#40ip+j*9SfutHqT@UE$J?=EHj}A3pUHlG98Hxg>TN0%RdqSjiksxq_$> zw*gLg3iq2j&h-ard*{DaLAcHcVu&@rziQliIiwf6yG<~@_4_S;{2^A>Z=v0UlHOEl ztvu%JhqKEyBK`CM8@2mG0dEC9nn!B_tDQNrh>ALzI%?il2WOctgIp3Y5lTIEP>c_+ zY3zfD++H0B6l%%t{bZ9qHJR@<-uJ>pxp_6*znmG-@7~ddQ;o_YI)tJ z7&S+YLhf-WWy~NVD}lwLsmJEmR{8^HrkTBj&5iZ>;F#5yT&lNPv*rriCjB;*V!fPi zAM%!k@VTsWzI(Hu&1xOh_uq_hF^z>3$r&~p4`Q0C*SY}@9#Eft)efi0u zv>j#Ste(w{sN{|es{3x!s`kfY7%Z@-bj#wMXBuwMLKZsC$wRJq767PUWqw+uuu|-s zEl=Ev9KJhvc~+6ip)W=iPo=4}MzySx_tH0Z!+i`J_Xt`C0(IJI{~ld?PXsOOUSNf_ zUAhFU^%e#k_L2-HZ11X?XSru|#!dvCd5WC@ks ziH=F&j4jf1LE#dLHj;P$zD*%ryhi%rpV67KPMee!j=6MMC}!JoWn}P;Y7| z=|MLdv(Q}R7c^h<+Vh}uE!l@!d;i-=I9d~a_ngkY3^4-=^)oJJ4bq{$+3Mt+GNg{( zkGbx+x^s$L(buazmsi@pzyH)fJ&R74+UG&IX4yQska`JO|IN-}iPt)J)LVncMbko@ z7|t?RByu0*NMNPBfx_fMYLgqwjDTzn1BaK$s@v*U^_pHd@CfAFG}-xPSe`D6*vafa z48F4RNxax`Rh&u}vdflc_9uN`lAbK!CJ#5)X^wruC?#gGeA|%L?@?!NDFN<2LzUZy zbtXy48hlh=epjsz_2B-3@4mEQ9xu_h!=`28k%vd=6QpYsCUGR$oy@)SGO3qt$eaj| zWvq{bVPldxfgaV4h`FLWbZy9rZLh*x6(lp7i#C5duvPZ0qCT7 z9k)v&c>6vZu`#ASy=Cj}XF{O(eN*7s0iss&I^|-o^}fy}_pvDI7MBhpMFU?6gFmjnMC~ju z_BYoJeR)JeYl7g9;k48oSSdcQcR-rjF{UjC1&)Jmf1azY<)0KK!%C-Lm0p!CflP?ET*w@S&Q%$#imxP~)}Jkjb{~0{GN4M1oj3ykefT zG&MM<-0jKt{4AM;+H!I6HLK=$yJs*82TWW$ci#Hs+jwsB`;ibn)VwB&n;3t^B9HA` zksr4I`O5o?J2Rbfj;lB{Vd!af_L&6*Cyiy@Ftuon(AGCNFk9~2O8*fX z@*TF}kqaW#?(YYY$r}5nfNPv6*_{plE@J333R!OoOgIX&4&SZtxkcb27K=#3EpT_|LpAU~Xi)odlA-HLlxxbW3Q;8UY_w{%|9H8%Qhn)Xj zg(M1dQ6*$xrzJHf4YaNuWFd~A`-(hmk#=41%n=>AgRzw#@;`__voT1vsJ8;ja)15a z%0^hrNsuBj@^1r&DO|)?!LLx7_@Z3@O9!duWB0E8c2)_=c<_={bF{-657-y*3qm(< zOZXg`TfWDJ|A&b|&pq6W278ifp5{Vp?kot>Oq1_vWL4j=rLu&x)Jw@79)}%@C4}S&0(is6B>ULkjN3i0?IS*9gyX3La$L6tB0vP;hSpyO3=Uyu1 zu!()U**W=xE%K;&tq>C(13+~C$ZMl~3r;}5>2YV1@5mKw04w#N*L#{$|K~Ivt*Cgj zZa;QwLy`r%W~JMZuN-(=+1Peu@m|+0Y&sbu_$a$r-7>R$%oD{D5-6vM$M~h6B4RXU zq0@jIfkjgex=m^XkKw@pRI-?K{8@4X87E9Do%jdbJUUwTzRi`wd8Y_RtJ3naV9RM(`l@_1Ncki2O;8`WTuW z`KH%JIk8I@c!n#{xhN_nnHpvRVsBo1P~M;^iQnan-1D2ljd_KYop`6`Xo2#3|GV5_ zfLc@)gcuSGwa$;)Rnua~0>r5~CBXjVJA}*F*ToKLsT-5Gl>Ql_H7v90u@A!xP)*ux=2vbfJc0ig#s_Y>l zsL(*X$Ex$puZM@KH%*1vAmb{!C17P=hm`zS*IX`J)AX5qi}S0v9un8S=sD_--a}IHBNN`${XK#;G#8S*?nBs=Nc@!tQaW?jzU z*PIu8F34Ojjwo3awl+oD98xO?)u>gVH1{+dk5$>=A-}G9q&^mR9)#Z*lQ6pyy9+7c zVc$Z#f^Xgu(r>@>it(w|Sy=%}7bh9l8!fbd`SH>pq}nFqMZCR)7sV&FOZ%E*yu!VP z63n!WKTEV)z%J~=dNFFjtDs*TYY0P0_0PPh6S>0p?HWxn=+*Tq$AV#>$Pz&%h(&qW z916%*>g~!R#%e&rmnc=Nf^~c}*DuB_nx_jF&;%REWp;*WQ1vt`{;#`Z93wRQP)>cs z$k#M5tSPfW3`lMb*Oe}IzGLxxXEr*EQSIL|_Xcag4+R9{GzqVMtWYV$_))DFtI-z} z$JOj_*!iqT3^bIbNl*_Z<87H-qjs=2e_V&qrF-L(wE^9AUJv`CHNPO{ez}-m5|S9) z_{>+myYhvS71iK*^>Y40JqmoTJ=M62i`o3L&IkKOgXTcr)ja+sQSV&T#HB4o%ZiH6 zj}R`>>NQMPiQMhqsGwWyo_yQRkm5||W~%gzx}oZEwZG5BZ*@!9@Zm^?(^pcDAc{JT z;d5!$aC#^O`fW))pNw++9O^-6Vj3-FHrTU~_~_b+K&)UX{p za-BI3F(WWa7pqQ-on(i^Vi48*uY5LpAC@B=JR86mg3@tl1lzkW=%RBZ*>&wiU$62j zWQfJNKS9*U*FCT`Hp>0l`6t^m#lw3#;R-pmiA+z@0W7|2A)wSdRgcJUQ#k=;>MmXO zvy?PFcKBZT?j-_WH5+w-ebcLEZmQVRb^!8H{8G6hi&P)T(o)N--uO1Sjyre3{V0N%lqY}}?^I)HA}YA6)ykFdCBF9Iy(m-|G0qk}6!N|o z*S6RpVnNV8#}2y#iqYElU2!Qru^0neglmJ43T2mlp=2@sJO)ztsMDQp2s`1EosIf=9=uwZvp($4 zvvcntlRFgIAO}fx;RVGU)$_w1Z0v^(NXSq=mBj!{1<}L68QxXDt~mtlNFa*lfj`vZ z)Ka`?mE;h|0s85bNoSfwH*;P|&W!vK7ord+Gbe>J^MQq;1MUyHkyY26ZwQP;i6mV0 zxs+gyObnJW=wcoA_tWNrJE6HHnC^H8^y3xJS$Y&3|ETUam(pa(U%p}vkjGa~lo4t0 z9`=9;_S=vIkV4Lby43-fSZIV0%&|P)EI~rTkCVG6(|k=U#SKq+V(5M@yeX|}F zt(Q)GJpbyLblQURCLTHRNc#8-;PprwG@xS(wlgr>-R zY*#VW?|ECq)hMq&wP~#ikcWT8)VOO(|8v^Yi$<3b75kY|>Be3*Z(7p0k@}t90uLMf z>u-Faa2@4>j{=E)Z;l23u5zVvc&76i90Go)ozTFgHf>z!)qD{?WLEHDS@7HDXU3*2 z(|umi@XAXShFK7@Z(*d3lAzDHnXS;vCqV6%(&?@%*cjQv89o$5N_r6-mSmrO6{I0` zc!(Fe#%xA&9vUZ{xYJ+pZmA@6*bIg3uFy6$dF=3_7G? zB$8EoU3IIa<3bmf12q&}9RX>Z+)vXEc6@_JkuFSR&nJ)9y2$gST+5r?dD?;D#=Ec;(8 zU9AenM9KVeF&kF%17I~(T;nMeRXc#zt9ku^QDX`>Im}k$Sivv0JtL1E&7Ylm-!eQ@ zn0LEmT$Yh@^O~dtwKL4$;As)T-7N2VVLf?eiN8@ySGz(U6-k+YU{XD43FoGCIPy`% z$^IfO2p`zCxM?M~w@ve-g=9fE{SE zlJb}JtX)i%?)R`F!O=>4d^PNsKt%*g87V3Y`RH>&m#e)J4voiS z{u2JFuO)WhI8QDA2su{z%unUvH|8YTHC38Ajk|yQ@a3&fC;pYKSnZ(3Z-zpDPSoX}gP_1pV$}6}Zb)dUjjJpDGmf>xVIalX>YC^u zN3UExmD?iE_0u%twHTqC2C*~2ANY{mL>6njZEeX8V9c0KIN^8oJ2B;7Q^~#i-gfl+ z+nvrB9ov7+Z5w>Wwr<#Nz3-itUe-UO8d}1~r>(4lK(IXiaTmbvG4H6`hWyXI^Z3S# z1q^Ns`49a>m`b~B8{{{cD$nNh?2FLP7BvR*_->q|Q4b{SE= z0kyZz1epE25%KQ@iN6pnF~45h9C@t&;o-xmov4PJo2fPi5hue}DsOa2xcIgenBk8D z9}Ht2w0-yQFK}|#Q{Q}e<>Pq!_r3|I^^MjT{GwL)tjjJ#?{*wAa<4Dg`Qtia@$?YT zNMCIKH(eMs+v+~xp>zo9`J8h@gqi3g_IzipuqQEaa?-EnbANPOK^I&}Qjd2-rPF3& zT!HV$`N{TL%IuPwYLiWYMDX&B4ZYjI$+_>bZ@Wt3iYmX3gj;CL@RSktRN&@YuvXCH!L+QU4u|~5xNroai>)}4m8@I|L zVj@&~YjwA$m)g>$%UMkT#ktU%f5Fkba35B(i1dr=yW762JGaWRp9Pj0)WKmWYrMr2 z99>r^o8Dzc+xj