From 65571bdd3eec6bec56600299001c86723e80fa87 Mon Sep 17 00:00:00 2001
From: dramanica It is usually advisable to plot the locations directly on the raster
that will be used to extract climatic variables, to see how the
locations fall within the discrete space of the raster. For this
@@ -179,9 +179,9 @@ SDMs with
tidymodels
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag() masks stats::lag()
#> ✖ recipes::step() masks stats::step()
-#> • Search for functions across packages at https://www.tidymodels.org/find/
+#> • Use tidymodels_prefer() to resolve common conflicts.
#> Loading required package: spatialsample
Accessing the data for this vignette: how to use
@@ -162,13 +162,13 @@
Preparing your data
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
lacerta <- st_as_sf(lacerta, coords = c("longitude", "latitude"))
-st_crs(lacerta) <- 4326
Preparing your data
library(pastclim)
download_dataset(dataset = "WorldClim_2.1_10m")
@@ -197,7 +197,7 @@
Preparing your data -9.8 43.3))"
)
-crs(iberia_poly) <- "lonlat"
+crs(iberia_poly) <- "+proj=longlat"
# crop the extent
land_mask <- crop(land_mask, iberia_poly)
# and mask to the polygon
@@ -231,32 +231,81 @@
Preparing your data
Now, we thin the observations to have one per cell in the raster (it -would be better if we had an equal area projection…):
+Before we start thinning the data we need to make sure that all our +data (points and rasters) have the same geographic coordinate reference +system (CRS) by projecting them. In some of the pipeline steps +(e.g. thinning data, measuring areas) using an equal area projection may +make a significant difference, especially for large-scale projects.
+You can use the website projectionwizard.org
(https://link.springer.com/chapter/10.1007/978-3-319-51835-0_9)
+to find an appropriate equal area projection for any region.
To define our projection within the code, we will use a proj4 string,
+which provides information on the type of projection, its parameters and
+the units of distance in which the new coordinates will be expressed (if
+you are using projectionwizard.org
it will provide yo with
+the string as well).
In this case, we will use a Albers Equal Area Conic projection +centred on the Iberian peninsula, with km as units. The proj4 string +is:
+iberia_proj4 <- "+proj=aea +lon_0=-4.0 +lat_1=36.8 +lat_2=42.6 +lat_0=39.7 +datum=WGS84 +units=m +no_defs"
For rasters (maps), we use the terra
function
+project
to change the CRS. We pass the raster object and
+the proj4 string as arguments:
+land_mask <- terra::project(land_mask, y = iberia_proj4)
Now we need to project the data points to the same CRS as the raster.
+We will do so using the appropriate sf
function:
+lacerta <- st_transform(lacerta, iberia_proj4)
Plotting the data, we will see that the shape of the land mask has +slightly changed following the new projection.
+
+
+ggplot() +
+ geom_spatraster(data = land_mask, aes(fill = land_mask_1985)) +
+ geom_sf(data = lacerta) +
+ guides(fill="none")
Now, we thin the observations to have one per cell in the raster +(given our project, each cell is approximately the same size):
+
set.seed(1234567)
lacerta <- thin_by_cell(lacerta, raster = land_mask)
nrow(lacerta)
-#> [1] 226
+#> [1] 233
+pres_data <- terra::extract(land_mask, lacerta)
+summary(pres_data)
+#> ID land_mask_1985
+#> Min. : 1 land:233
+#> 1st Qu.: 59
+#> Median :117
+#> Mean :117
+#> 3rd Qu.:175
+#> Max. :233
ggplot() +
geom_spatraster(data = land_mask, aes(fill = land_mask_1985)) +
geom_sf(data = lacerta) +
guides(fill="none")
Now, we thin further to remove points that are closer than 20km.
-However, note that the standard map units for a ‘lonlat’ projection are
-meters. tidysdm
provides a convening conversion function,
-km2m()
, to avoid having to write lots of zeroes):
+Now, we thin further to remove points that are closer than 20km. As +the units of our projection are m (the default for most projections), we +use a a convenient conversion function,
+km2m()
, to avoid +having to write lots of zeroes:+#> [1] 112set.seed(1234567) lacerta_thin <- thin_by_dist(lacerta, dist_min = km2m(20)) nrow(lacerta_thin) -#> [1] 111
Let’s see what we have left of our points:
-+ggplot() + geom_spatraster(data = land_mask, aes(fill = land_mask_1985)) + geom_sf(data = lacerta_thin) + @@ -275,7 +324,7 @@
Thinning stephttps://doi.org/10.15468/dl.53js5z: -
+++library(rgbif) occ_download_get(key = "0121761-240321170329656", path = tempdir()) library(readr) @@ -284,33 +333,53 @@
Thinning step# keep the necessary columns lacertidae_background <- backg_distrib %>% select(gbifID, decimalLatitude, decimalLongitude) %>% - rename(ID = gbifID, latitude = decimalLatitude, longitude = decimalLongitude) + rename(ID = gbifID, latitude = decimalLatitude, longitude = decimalLongitude)
In this case as well, we need to use the appropriate projection (the +same defined before) for the background. If the projections do not +correspond the analyses will stop giving an error message.
++ +st_crs(lacertidae_background) <- "+proj=longlat" +lacertidae_background <- st_transform(lacertidae_background, crs = iberia_proj4)+# convert to an sf object lacertidae_background <- st_as_sf(lacertidae_background, coords = c("longitude", "latitude")) -st_crs(lacertidae_background) <- 4326
#> Warning: st_crs<- : replacing crs does not reproject data; use st_transform for +#> that
We need to convert these observations into a raster whose values are the number of records (which will be later used to determine how likely -each cell is to be used as a background point):
-+-lacertidae_background_raster <- rasterize(lacertidae_background, land_mask, +each cell is to be used as a background point). We will also mask the +resulting background raster to match the land mask of interest. +
+lacertidae_background_raster <- mask(lacertidae_background_raster, + land_mask) +ggplot() + + geom_spatraster(data = lacertidae_background_raster, aes(fill = count)) + + scale_fill_viridis_b( na.value = "transparent")+lacertidae_background_raster <- rasterize(lacertidae_background, + land_mask, fun = "count") - -plot(lacertidae_background_raster)
+guides(fill="none") +#> <Guides[1] ggproto object> +#> +#> fill : "none"
We can see that the sampling is far from random, with certain locations having very large number of records. We can now sample the background, using the ‘bias’ method to represent this heterogeneity in sampling effort:
-++ return_pres = TRUE) +#> Warning in sample_background(data = lacerta_thin, raster = lacertidae_background_raster, : There are fewer available cells for raster 'NA' (112 presences) than the requested 336 background points. Only 9 will be returned.set.seed(1234567) lacerta_thin <- sample_background(data = lacerta_thin, raster = lacertidae_background_raster, n = 3 * nrow(lacerta_thin), method = "bias", class_label = "background", - return_pres = TRUE)
Let’s see our presences and background:
-+ggplot() + geom_spatraster(data = land_mask, aes(fill = land_mask_1985)) + geom_sf(data = lacerta_thin, aes(col = class)) + @@ -321,7 +390,7 @@
Thinning step
+download_dataset("WorldClim_2.1_10m") climate_vars <- get_vars_for_dataset("WorldClim_2.1_10m") climate_present <- pastclim::region_slice( @@ -339,44 +408,102 @@
Thinning steparticle on additional features of tidymodels with tidysdm. +
And now we project the climate variables in the same way as we did +for all previous spatial data:
++climate_present <- terra::project(climate_present, y = iberia_proj4)
Next, we extract climate for all presences and background points:
-++lacerta_thin <- lacerta_thin %>% bind_cols(terra::extract(climate_present, lacerta_thin, ID = FALSE))
Before going forward with the analysis, we should make sure that +there are no missing values in the climate that we extracted:
+++summary(lacerta_thin) +#> class geometry bio01 bio02 +#> presence :112 POINT :121 Min. : 7.222 Min. : 5.920 +#> background: 9 epsg:NA : 0 1st Qu.:11.556 1st Qu.: 8.822 +#> +proj=aea ...: 0 Median :13.040 Median : 9.280 +#> Mean :12.931 Mean : 9.301 +#> 3rd Qu.:14.297 3rd Qu.: 9.881 +#> Max. :16.672 Max. :12.124 +#> bio03 bio04 bio05 bio06 +#> Min. :32.96 Min. :320.6 Min. :21.22 Min. :-2.665 +#> 1st Qu.:38.80 1st Qu.:424.3 1st Qu.:23.94 1st Qu.: 1.409 +#> Median :40.77 Median :493.5 Median :25.78 Median : 3.184 +#> Mean :40.89 Mean :491.7 Mean :25.95 Mean : 3.034 +#> 3rd Qu.:43.49 3rd Qu.:558.5 3rd Qu.:27.54 3rd Qu.: 5.030 +#> Max. :47.59 Max. :667.4 Max. :32.72 Max. : 8.221 +#> bio07 bio08 bio09 bio10 +#> Min. :14.74 Min. : 1.817 Min. :14.56 Min. :14.57 +#> 1st Qu.:20.78 1st Qu.: 6.940 1st Qu.:17.96 1st Qu.:18.25 +#> Median :23.13 Median : 8.485 Median :18.96 Median :19.02 +#> Mean :22.91 Mean : 8.616 Mean :19.13 Mean :19.29 +#> 3rd Qu.:25.03 3rd Qu.:10.440 3rd Qu.:20.25 3rd Qu.:20.39 +#> Max. :30.30 Max. :17.639 Max. :23.72 Max. :23.87 +#> bio11 bio12 bio13 bio14 +#> Min. : 0.8458 Min. : 397.5 Min. : 48.67 Min. : 2.00 +#> 1st Qu.: 5.8650 1st Qu.: 769.1 1st Qu.:106.65 1st Qu.:13.18 +#> Median : 7.5043 Median :1056.4 Median :143.07 Median :18.88 +#> Mean : 7.3158 Mean :1049.5 Mean :146.48 Mean :22.34 +#> 3rd Qu.: 9.2647 3rd Qu.:1347.5 3rd Qu.:186.54 3rd Qu.:30.63 +#> Max. :11.7415 Max. :1732.6 Max. :249.79 Max. :51.97 +#> bio15 bio16 bio17 bio18 +#> Min. :24.92 Min. :128.6 Min. : 19.38 Min. : 25.38 +#> 1st Qu.:40.18 1st Qu.:301.1 1st Qu.: 65.13 1st Qu.: 68.60 +#> Median :51.05 Median :401.9 Median : 90.20 Median : 99.58 +#> Mean :47.88 Mean :412.4 Mean : 95.30 Mean :105.46 +#> 3rd Qu.:56.32 3rd Qu.:528.8 3rd Qu.:130.62 3rd Qu.:153.24 +#> Max. :69.06 Max. :695.0 Max. :178.73 Max. :195.92 +#> bio19 altitude +#> Min. :106.2 Min. : 22.58 +#> 1st Qu.:283.8 1st Qu.: 192.58 +#> Median :387.5 Median : 512.13 +#> Mean :397.9 Mean : 538.88 +#> 3rd Qu.:528.8 3rd Qu.: 733.15 +#> Max. :695.0 Max. :1477.87
We can see that there are no missing values in any of the extracted +climate variables. If that was not the case, we would have to go back to +the climate raster and homogenise the NAs across layers +(i.e. variables). This can be achieved either by setting the same cells +to NA in all layers (including the land mask that we used to thin the +data), or by interpolating the layers with less information to fill the +gaps (e.g. cloud cover in some remote sensed data). interpolate the +missing
Based on this paper (https://doi.org/10.1007/s10531-010-9865-2), we are interested in these variables: “bio06”, “bio05”, “bio13”, “bio14”, “bio15”. We can visualise the differences between presences and the background using violin plots:
-+-lacerta_thin %>% plot_pres_vs_bg(class)
+
We can see that all the variables of interest do seem to have a different distribution between presences and the background. We can formally quantify the mismatch between the two by computing the overlap:
-++#> bio19 bio16 bio12 bio13 bio02 bio15 bio07 bio04 +#> 0.9959022 0.9362035 0.9268597 0.9140456 0.7975679 0.7452740 0.7092859 0.6983330 +#> bio05 bio17 bio14 bio06 bio03 bio18 bio11 bio08 +#> 0.6646394 0.5635174 0.5631364 0.4638880 0.4473745 0.4445754 0.4186998 0.3892218 +#> bio09 bio10 bio01 altitude +#> 0.2898833 0.2812704 0.2706408 0.2580447lacerta_thin %>% dist_pres_vs_bg(class) -#> bio09 bio12 bio16 bio19 bio13 bio05 bio10 -#> 0.43907819 0.41888524 0.41487381 0.40742724 0.40492411 0.38854703 0.38610145 -#> bio02 bio07 bio04 bio08 bio17 bio15 bio18 -#> 0.35191109 0.35036167 0.32450555 0.31879785 0.28143659 0.27152095 0.25007068 -#> bio01 bio14 bio03 bio11 altitude bio06 -#> 0.24589097 0.24294699 0.18414624 0.11169528 0.07271380 0.06742951
Again, we can see that the variables of interest seem good candidates with a clear signal. Let us then focus on those variables:
-+suggested_vars <- c("bio06", "bio05", "bio13", "bio14", "bio15")
Environmental variables are often highly correlated, and collinearity is an issue for several types of models. We can inspect the correlation among variables with:
-+- +pairs(climate_present[[suggested_vars]])
We can see that some variables have rather high correlation (e.g. bio05 vs bio14). We can subset to variables below a certain threshold correlation (e.g. 0.7) with:
-+climate_present <- climate_present[[suggested_vars]] vars_uncor <- filter_collinear(climate_present, cutoff = 0.7, method = "cor_caret") @@ -388,7 +515,7 @@
Thinning step
+lacerta_thin <- lacerta_thin %>% select(all_of(c(vars_uncor, "class"))) climate_present <- climate_present[[vars_uncor]] names(climate_present) # added to highlight which variables are retained in the end @@ -405,7 +532,7 @@
is automatically replaced byFit the model by cross-validationgeometry
X
andY
columns which are assigned a role ofcoords
, and thus not used as predictors): -+lacerta_rec <- recipe(lacerta_thin, formula = class ~ .) lacerta_rec #> @@ -420,7 +547,7 @@
Fit the model by cross-validation -
+lacerta_thin %>% check_sdm_presence(class) #> [1] TRUE
We now build a
workflow_set
of different models, @@ -437,7 +564,7 @@Fit the model by cross-validationgam_formula() due to the non-standard formula notation of GAMs (see the help of
sdm_spec_gam()
for an example of how to do this). -+lacerta_models <- # create the workflow_set workflow_set( @@ -466,7 +593,7 @@
object suitable toFit the model by cross-validationrsample
tisysdm
with the functionblockcv2rsample
. -+library(tidysdm) set.seed(100) #lacerta_cv <- spatial_block_cv(lacerta_thin, v = 5) @@ -477,7 +604,7 @@
Fit the model by cross-validationWe can now use the block CV folds to tune and assess the models (to keep computations fast, we will only explore 3 combination of hyperparameters per model; this is far too little in real life!): -
++#> ✔ 4 of 4 tuning: default_maxent (831ms)set.seed(1234567) lacerta_models <- lacerta_models %>% @@ -487,19 +614,23 @@
Fit the model by cross-validation ) #> i No tuning parameters. `fit_resamples()` will be attempted #> i 1 of 4 resampling: default_glm -#> ✔ 1 of 4 resampling: default_glm (190ms) +#> → A | warning: glm.fit: algorithm did not converge, glm.fit: fitted probabilities numerically 0 or 1 occurred +#> There were issues with some computations A: x1 +#> There were issues with some computations A: x3 +#> +#> ✔ 1 of 4 resampling: default_glm (260ms) #> i 2 of 4 tuning: default_rf #> i Creating pre-processing data to finalize unknown parameter: mtry -#> ✔ 2 of 4 tuning: default_rf (848ms) +#> ✔ 2 of 4 tuning: default_rf (465ms) #> i 3 of 4 tuning: default_gbm #> i Creating pre-processing data to finalize unknown parameter: mtry -#> ✔ 3 of 4 tuning: default_gbm (4s) +#> ✔ 3 of 4 tuning: default_gbm (1.7s) #> i 4 of 4 tuning: default_maxent -#> ✔ 4 of 4 tuning: default_maxent (1.2s)
Note that
-workflow_set
correctly detects that we have no tuning parameters for glm. We can have a look at the performance of our models with:+autoplot(lacerta_models)
Now let’s create an ensemble, selecting the best set of parameters @@ -509,9 +640,11 @@
Fit the model by cross-validation -
+lacerta_ensemble <- simple_ensemble() %>% add_member(lacerta_models, metric = "boyce_cont") +#> Warning: glm.fit: algorithm did not converge +#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred lacerta_ensemble #> A simple_ensemble of models #> @@ -529,34 +662,34 @@
Fit the model by cross-validation#> Metric used to tune workflows: #> • boyce_cont
And visualise it
-+autoplot(lacerta_ensemble)
A tabular form of the model metrics can be obtained with:
-++#> 1 default_glm boyce_cont 0.789 0.211 2 +#> 2 default_glm roc_auc 1 0 3 +#> 3 default_glm tss_max 0.969 0.00942 3 +#> 4 default_rf boyce_cont 0.845 0.0835 3 +#> 5 default_rf roc_auc 0.978 0.0222 3 +#> 6 default_rf tss_max 0.943 0.0469 3 +#> 7 default_gbm boyce_cont 0.890 0.0651 3 +#> 8 default_gbm roc_auc 0.986 0.0139 3 +#> 9 default_gbm tss_max 0.931 0.0458 3 +#> 10 default_maxent boyce_cont 0.684 0.193 3 +#> 11 default_maxent roc_auc 0.978 0.0222 3 +#> 12 default_maxent tss_max 0.943 0.0469 3lacerta_ensemble %>% collect_metrics() #> # A tibble: 12 × 5 #> wflow_id .metric mean std_err n #> <chr> <chr> <dbl> <dbl> <int> -#> 1 default_glm boyce_cont 0.573 0.115 3 -#> 2 default_glm roc_auc 0.775 0.0138 3 -#> 3 default_glm tss_max 0.486 0.0337 3 -#> 4 default_rf boyce_cont 0.709 0.0856 3 -#> 5 default_rf roc_auc 0.794 0.00648 3 -#> 6 default_rf tss_max 0.537 0.0363 3 -#> 7 default_gbm boyce_cont 0.659 0.0472 3 -#> 8 default_gbm roc_auc 0.789 0.00707 3 -#> 9 default_gbm tss_max 0.524 0.0152 3 -#> 10 default_maxent boyce_cont 0.651 0.157 3 -#> 11 default_maxent roc_auc 0.804 0.00653 3 -#> 12 default_maxent tss_max 0.572 0.0111 3
Projecting to the present
We can now make predictions with this ensemble (using the default option of taking the mean of the predictions from each model).
-+prediction_present <- predict_raster(lacerta_ensemble, climate_present) ggplot() + geom_spatraster(data = prediction_present, aes(fill = mean)) + @@ -565,14 +698,15 @@
Projecting to the present geom_sf(data = lacerta_thin %>% filter(class == "presence"))
We can subset the ensemble to only use the best models, based on the -Boyce continuous index, by setting a minimum threshold of 0.7 for that -metric. We will also take the median of the available model predictions -(instead of the mean, which is the default). The plot does not change -much (the models are quite consistent).
-+Boyce continuous index, by setting a minimum threshold of 0.6 for that +metric (this is somewhat low, for a real analysis we would recommend a +higher value of 0.7 or higher). We will also take the median of the +available model predictions (instead of the mean, which is the default). +The plot does not change much (the models are quite consistent). +prediction_present_boyce <- predict_raster(lacerta_ensemble, climate_present, - metric_thresh = c("boyce_cont", 0.7), + metric_thresh = c("boyce_cont", 0.6), fun = "median" ) ggplot() + @@ -584,24 +718,24 @@
Projecting to the present -
+lacerta_ensemble <- calib_class_thresh(lacerta_ensemble, class_thresh = "tss_max", - metric_thresh = c("boyce_cont", 0.7) + metric_thresh = c("boyce_cont", 0.6) )
And now we can predict for the whole continent:
-+- +prediction_present_binary <- predict_raster(lacerta_ensemble, climate_present, type = "class", class_thresh = c("tss_max"), - metric_thresh = c("boyce_cont", 0.7) + metric_thresh = c("boyce_cont", 0.6) ) ggplot() + geom_spatraster(data = prediction_present_binary, aes(fill = binary_mean)) + geom_sf(data = lacerta_thin %>% filter(class == "presence")) + scale_fill_discrete(na.value = "transparent")
Projecting to the future @@ -612,16 +746,16 @@
Projecting to the future -
+download_dataset("WorldClim_2.1_HadGEM3-GC31-LL_ssp245_10m")
Let’s see what times are available:
-+get_time_ce_steps("WorldClim_2.1_HadGEM3-GC31-LL_ssp245_10m")
#> [1] 2030 2050 2070 2090
We will predict for 2090, the further prediction in the future that is available.
Let’s now check the available variables:
-+get_vars_for_dataset("WorldClim_2.1_HadGEM3-GC31-LL_ssp245_10m")
#> [1] "bio01" "bio02" "bio03" "bio04" "bio05" "bio06" "bio07" "bio08" "bio09" #> [10] "bio10" "bio11" "bio12" "bio13" "bio14" "bio15" "bio16" "bio17" "bio18" @@ -631,15 +765,19 @@
Projecting to the future -
++climate_future <- pastclim::region_slice( time_ce = 2090, bio_variables = vars_uncor, data = "WorldClim_2.1_HadGEM3-GC31-LL_ssp245_10m", crop = iberia_poly )
Project the climatic raster with the same projection that we have +been using for the analysis:
++climate_future <- terra::project(climate_future, y = iberia_proj4)
And predict using the ensemble:
-+prediction_future <- predict_raster(lacerta_ensemble, climate_future) ggplot() + @@ -657,7 +795,7 @@
Dealing with extrapolation
tidysdm
offers a couple of approaches to deal with this problem. The simplest one is that we can clamp the environmental variables to stay within the limits observed in the calibration set: -+- +climate_future_clamped <- clamp_predictors(climate_future, training = lacerta_thin, .col= class) @@ -667,7 +805,7 @@
Dealing with extrapolationggplot() + geom_spatraster(data = prediction_future_clamped, aes(fill = mean)) + scale_fill_terrain_c()
The predictions seem to have changed very little.
An alternative is to allow values to exceed the ranges of the calibration set, but compute the Multivariate environmental similarity @@ -675,14 +813,14 @@
Dealing with extrapolation
We estimate the MESS for the same future time slice used above:
-+- +lacerta_mess_future <- extrapol_mess(x = climate_future, training = lacerta_thin, .col = "class") ggplot() + geom_spatraster(data = lacerta_mess_future) + scale_fill_viridis_b(na.value = "transparent")
Extrapolation occurs in areas where MESS values are negative, with the magnitude of the negative values indicating how extreme is in the interpolation. From this plot, we can see that the area of extrapolation @@ -690,7 +828,7 @@
Dealing with extrapolation
We can now overlay MESS values with current prediction to visualize areas characterized by spatial extrapolation.
-+- +# subset mess lacerta_mess_future_subset <- lacerta_mess_future lacerta_mess_future_subset[lacerta_mess_future_subset >= 0] <- NA @@ -704,7 +842,7 @@
Dealing with extrapolation scale_fill_viridis_b(na.value = "transparent") + geom_sf(data = lacerta_mess_future_subset, fill= "lightgray", alpha = 0.5, linewidth = 0.5)
Note that clamping and MESS are not only useful when making predictions into the future, but also into the past and present (in the latter case, it allows us to make sure that the @@ -728,7 +866,7 @@
Visualising the co marginal prediction. We can then plot the predictions against the values of the variable of interest. For example, to investigate the contribution of
bio05
, we would: -+- +bio05_prof <- lacerta_rec %>% step_profile(-bio05, profile = vars(bio05)) %>% prep(training = lacerta_thin) @@ -742,7 +880,7 @@
Visualising the co ggplot(bio05_data, aes(x = bio05, y = pred)) + geom_point(alpha = .5, cex = 1)
It is also possible to use DALEX,to explore
@@ -757,11 +895,12 @@tidysdm
models; see more details in the tidymodels additions article.Repeated ensembles
+will just use two fast models to speed up the process, and use +pseudo-absences instead of background. ++#> ✔ 2 of 2 tuning: default_maxent (1.2s)# empty object to store the simple ensembles that we will create ensemble_list <- list() -set.seed(123) # make sure you set the seed OUTSIDE the loop +set.seed(1234) # make sure you set the seed OUTSIDE the loop for (i_repeat in 1:3) { # thin the data lacerta_thin_rep <- thin_by_cell(lacerta, raster = climate_present) @@ -776,7 +915,7 @@
Repeated ensembleslacerta_thin_rep <- lacerta_thin_rep %>% bind_cols(terra::extract(climate_present, lacerta_thin_rep, ID = FALSE)) # create folds - lacerta_thin_rep_cv <- spatial_block_cv(lacerta_thin_rep, v = 5) + lacerta_thin_rep_cv <- spatial_block_cv(lacerta_thin_rep, v = 3, n = 5) # create a recipe lacerta_thin_rep_rec <- recipe(lacerta_thin_rep, formula = class ~ .) # create a workflow_set @@ -800,7 +939,7 @@
Repeated ensembleslacerta_thin_rep_models <- lacerta_thin_rep_models %>% workflow_map("tune_grid", - resamples = lacerta_thin_rep_cv, grid = 10, + resamples = lacerta_thin_rep_cv, grid = 3, metrics = sdm_metric_set(), verbose = TRUE ) # make an simple ensemble and add it to the list @@ -809,21 +948,21 @@
Repeated ensembles} #> i No tuning parameters. `fit_resamples()` will be attempted #> i 1 of 2 resampling: default_glm -#> ✔ 1 of 2 resampling: default_glm (226ms) +#> ✔ 1 of 2 resampling: default_glm (161ms) #> i 2 of 2 tuning: default_maxent -#> ✔ 2 of 2 tuning: default_maxent (7.1s) +#> ✔ 2 of 2 tuning: default_maxent (1.2s) #> i No tuning parameters. `fit_resamples()` will be attempted #> i 1 of 2 resampling: default_glm -#> ✔ 1 of 2 resampling: default_glm (226ms) +#> ✔ 1 of 2 resampling: default_glm (170ms) #> i 2 of 2 tuning: default_maxent -#> ✔ 2 of 2 tuning: default_maxent (6.8s) +#> ✔ 2 of 2 tuning: default_maxent (1.1s) #> i No tuning parameters. `fit_resamples()` will be attempted #> i 1 of 2 resampling: default_glm -#> ✔ 1 of 2 resampling: default_glm (225ms) +#> ✔ 1 of 2 resampling: default_glm (171ms) #> i 2 of 2 tuning: default_maxent -#> ✔ 2 of 2 tuning: default_maxent (7.4s)
Now we can create a
-repeat_ensemble
from the list:+lacerta_rep_ens <- repeat_ensemble() %>% add_repeat(ensemble_list) lacerta_rep_ens #> A repeat_ensemble of models @@ -845,16 +984,21 @@
function forRepeated ensemblesautoplot()
repeated_ensemble
objects. -We can then predict in the usual way (we will take the mean and -median of all models):
-diff --git a/dev/articles/a0_tidysdm_overview_files/figure-html/autoplot_ens-1.png b/dev/articles/a0_tidysdm_overview_files/figure-html/autoplot_ens-1.png index 9e96e4a4aa6e5af36c598ad33b3150629a7be30e..eb1b28bf57071cad86cd7626ddcb6b45dfae4fd3 100644 GIT binary patch literal 56001 zcmd431z1&k^ewtEPzg~1X;9qK-5{k>8x)ZS2?3Fk?#2Q{0oh1*NrQBQih?xKf`D{O zgY=sV&$;!z_wIfF8{gaCIfBSud#zv0Ip!E+{$435$`ayH;-OF|LX4c0G75EC5`{X6 zI*SWGnF?yNgfF;;3bIo06<#~HtGwaMIU6}GdlZVs5BVP_cxu%cg}Q>mNZnL*iCr9a zcOjV?K3eTqv~nss_u}m-oUhOPlbHfYq^5HGl`h`AD4o4=C0BdTBF!|dCz+We can then predict in the usual way. We will take the mean and +median of all models, without filtering by performance, and plot the +results:
+- + +lacerta_rep_ens <- predict_raster(lacerta_rep_ens, climate_present, fun = c("mean", "median") ) ggplot() + geom_spatraster(data = lacerta_rep_ens, aes(fill = median)) + scale_fill_terrain_c()
Note that the predictions are quite similar to the ones we obtained +before, but the predicted suitable range is somewhat larger, probably +because we included models that are not very good (as we did not filter +by performance) in the ensemble.
*|-}+cy^f~eO+Yjvjp^N6JlKB{I{T|L9A0JnZz`1g| zU3_m*HNr-DsK&iaUS3{V7uSd6*Pp0J-72T1_pa?7%=hKSSk-K*$Jw_FbML_v z=V}$n?Ck99E*5o=yprlkQy?Rx61d~JyAe`v%eKtP%G&+GpnlS+Dl8#^Gre1DLZmVs9yUd}(R5zNB z-1Ss*n25)p4 ^y2Y!XV(bG0W1DekUqffujO#wud9> z!E&@G?zWMk=Y8OaJlo#ZQujo#MBeY=DusaMf%{sbjo&NJzPz>jXCk=7W5>?&XZigz zM6{}}FSGa5xI625A8d-OxAIynlAZVCcUgIe>2hDL(J+eH-}ktiAg)U!w(-_&`}eTv z(5a7Z^Y(hgcrT0LDm4dNL%ZXzIY|g`drK`Jz*sBX_9eIvAD+FDYrVId ?W95M?PuaO^hS^GaTCzmebI*yv86S8N&* z?>5h@>oWFocRk^Vt)453jgy0;uT5|?r?8NxEml~?awQrH?%8wnyBf1uV5ojFtBr z&okE=sd U;*^ZMltx+j3q3bOxyQRTpzhO_MDMXmVZQ%LUPZttN>?_fE z{Yi><|1E56VzMABu^Y1QRafCL(GWOTZks~uzQ_$ToIO}MA%;r?yQL+Zy&DFt VL z3-`QtL*qTSv>ca5bmK%Ev${;+R**rlv$NwfZYGI?8^@+v*Y3y=iLUt6?yo7_jpC73 zd4Gp{?+vxEGBRft`8DcBC~Zy63>a+H+y`ZS6eh9vjeo=n+w*xJdQ|9e48uB}AFe9n z;o%`U?^oiytakt_D#xj}dw*x0-+DyUakld&jCH<<<9upyG2f1dxaZ2JGj8Lr7&jL5 zkF;UlBl;Ht-H*PiJP_T>FKndq)~WX+>JRqXxCW!ddf}|;k3!>UxP3j3RcT%Kr4sYL zoGA2MtbP0UGjtw_H-4Aj*7w@EZh(4d+x`A-_uH$t`K*2ll|1=T*xW1hiu%U5uzIiT zP^l%q>)In-n@^|BzbC(>KJVuto*=$6NsK%(@8edbI AAQJzRI9OuA=vSZgkSjyXeQP%1Xx}7k!oqWialDB!q35pPa4Hf>+M@YY z%*- y#^nGTu|aazSZ%q{gvd*Ewot8b fA8-jN(U{&my|?d*L!D5B3wxtS zwZE&D#n&UL!@d< $Mjpj2`4zpiJW)c>=x~**iz<>vqaj&SK^YZd;e!CVdhh?}EnJ9TB zk>-Ze?9NPr#^cBD8Scfmx}s5uBgk$n8})Jz)c8g~!6OBC*K2HSJiQ-CB`A;OsN1(I zG4IXNuA6R)i7D&I(JE?)wI=s%2%=%dxNQ$RXK$SIt9OH_h{nLfW)wbw0mP;v8EE>C zbmJam7E(wak6X!IMdK-lHWvEtDSur$Tn+d;u8pJPw4^jtoa#1sA9uYrQs{dB_&ug? z3~~ulYjS`Am(MC!lh9JGxr@qScT=0)S{kBC^`y6|UTrX6he72i|8)Y>Mdbabp3<90 zszY~sc~n1tJG+Iu{6U!&CBOMQJUqOO*_2@O&o9nlQ!`c4S$KJi^p?+{e7>-+d`O(I zwhoq8P*7 z8eJpV1*PXcHf6#c?6&A1Ca1@{$!DwQb^#c1XcwMg;^fS4;5=Ly62Dl`K9Q`E(@1H- z%gU;N&C-cf(n@apg+v>R{=D8o;}*1<=))}do^kOcN|FsQRP38KZz5aZzKT8fwpgBO z!{R|_O259$y(fLfU&7;$P1w%d+#IEdL&k{Lo~7lFLTTMfhY(fnhPP;B6ngvAR8`TK zP0F&h&ZOH|Rb!hS&G{~da2IohM@Q0MMJ~>4p?^(Agomq7>`MjFtWsysFL)nqdv9*n z9UZ)|G W^iZr5>Zk*5YpmMbwi$&9j zxv46*d1s8m4ZH0vmGMLR75nq+ISVuMcNmkNLSu#Xr?g?xbxT2q1rNx+FS$aj@2&Gb z64}nqhbM73!ue39>HUkrZqa{6Q;lxab35+vaJ#N$$!X!~8QM aCA;INX_G6%!kdJvvzS)>r DT1De HQyy#w027CC35s0$p>r`$r zD9s1-GyKDqj-lvv_2oxp)=waMUS|dr7wo;C5=1MOm2j|VBC?QEv;m=B4>D%=YtGt? z{JKMtIsTKV#b6Xkc#-3Ry4c#cYtxo50AViI9n~V}2tzu`pb;J#YUCRRAU~qidNj{Q zC`S&RY+ZMtj{M|S(p$H!#Z0HA!QMad-x+HpAhQh?nFPCS&Kr-kyHr1Ji=JbY*m;#2 z@ac!d{>trjDNYi+%)~@VUX%75JE})TkJ~ni5)STQH@XuJRRMp6cFgcbz{_H$JuaHx zA{>^SPX0ldsb$lE_o~16^_@895zM@h2Z~bmZ>Zh0g)kzAMPL^a#rBnR%gZ%6t5=>) zw?^52K7?!=hDIKYg^B5{D!2GLgyFU&5xrwbnY~m!-K(MB3y3QHU~fC7rGD%K5_P{2 z-;)#Bdg&^ #_2Qa!AK>vsX=T)S>$>3kyrv&kFl{nk@i@W{7K@ zEypw2zaNkVXXjPQBjKktikIJh`i^j%BBtvdvkKNNjoHUH^zS;-jq^U-sf0C+;10-+ zd4p>JxB3x4-g=0@KLD5zMxBP?gVe$ipb9BSHOzJ{nbHnfV3srdqh5QtdTw)Rkq7r~ zEb;L!hN>=asbjx!y8zg7+h1**a;n;de5FniG>>4#+_I6ZC|;BKf=1fz$)-@MCsZiZ z)@>ACO9d31P?y|;b*Qk=mkTql0)?34LT=eSV5R(m;Gw0y?t2osVXPVe9i?|6+~&fn zv+oeoKmg9irzhKpziZ}I{T4w6E-EUDdz7(x%*GT(no``O9H~&D^v(G2LCMwCb^d#b zUtYy@wE4Xg0tC06M_ki=iRki>AhtAEg2XgY+F(VA0|x}MAM9>suYOY4gq;YhZL>G& z_8GB*jc#NaiXUvKtQ# q(cGBN%FgTZzj&?K^n0?*C Z>K3XB>NjHzAvBnp>+yfcxA>@;fN^~I}}AKuisZ)cU43-ucFX+dS$TTdjB^KXTg z)pC 1(9C_ukMb8Uyt( z#+Lst-6MC|Vnh&!FBuB2o%I G|by7RA_LdkN(;Xp1izEeh!0YPfO^l2ku) zdMtL{?w^ZFiJoe@<#P5?rvvv{Ps8JHZ1?XF5Zycx%@LG g-Rdy1@zZa~wL6DxvwRJFCUuvjdKMvns@yqPvQy?x2v-rktp z7b;;qe0-`8g}W^+Eh-TE5#$GNT-XFqsDgFan7zMqa7d!zaHqKR $YqF|PNzB{bp`f6cJ}l9;g#82s zPWYbKx8>R0g0*E+)QA5)qxoNOia`G%6MKcfS}&iVKEqiey*o1}CkHo!iCZsUgEpA4 z*EkXov#{mYL4x`ZKWE4fT#F4YvG3yJ^NNd$RU-_t_IbnO<2lNP;y# ?39AM&O%PmmzfJiNLyPN8lVm-8JXZCv&~Dnxo{1r8ZtQ=--U!Ey?H}`Db)FQ z9CoR%GYbmBMrXS8N;9vtx>|JSzzZ(+zLF%ge?bju5i}5zEz8^e71nxsdbh()qr&U| z|3gaZWOATP**@7s?LYm1X1;)zT26_|yc8to_xvh_Z<38r?7ZDWm8*3mwsIl DQyh4{w-enLhN|JL)oKf%9^HzP6}soJ2!gf+%s~%G*dXbYu_m`UCh7v zEXAvg0$867zi6Qz i5d8C7@3_9hB7^y{( G!WKJ%>BHh_;F_LY1Bi`{~%3_^EZh7 zpXAs7*WKyg7#JSTP&TAHv*l|Ld1MXL@9i*BGN39}U4X5CxWv3<8ZEr!Ai`eKk-V1) z^e=N>YF2YXLa({K5H*CZqMn&q8|pnp2z(l@gvZ07{2bMDi=f{B23LW0ry;q_t)cD< z$xm(Hky;cf-)$^1v{XN&DAGpc7EmtKP9 o_K50;K4HU$!}^uI>#f z7zpyU9gj?ouc~@_R*HVNt|~0?1~;YvD9;9T2wA*$DzdPsXv$8FnjuI10pL#3v!@`6 zepih!w5u9!Zf?e+!@|NcQ&UshJ=fRPD6`m(KRKF@ km-+@t5aAbbp zEMY_&+^AdvYv$LaET}X>bJFw|rRhV;cYW)}&>ErTa%jWMN`iW>^jz(jFmsXkUh~(7 z?giSqP-9{@y&B-QZPu 9{ZA6qeI#-qO@`UrQ_MZ)#JkoiD1cuGTBjnhSUYm)9gV`X5(+ zAVZ^2lJXRLcqF^?uLH>1q7^D~l^ pXdRtS2!~{n=GGJ$h6ZFt@DFrYnyAN(vAE zQL+3qE&LL*G0CEsK#lV~JXZLr56!V&awUykxFYO*#Q_ccM$=dJ75E~T8(VGEUEx|V zgoDbm$ j{IPO;-8njs?#uV0TQ9r{EGQ5nzv*J)z^p*@Sl zd*tcUphinhrnelfR0B1qKT9nekch)~X|g#WoNj+!;DlNN=;D&cZIm=34st+NcR_)g zd3dk~T%~}fK3qff3i51{p``}(xCr^#JMC9L_M*QA>fvGUqXV_R64&*|mLt_#ugN(3 zdNNg}R$Z_oNFHyr!54)u(R`4_M{3+H_W!IR)iw|bk?%HXh3%w>85LE9sQWVRr|USHSe5i)$aUZTy}M-XosDx6H^XtE573NUv)lUmIw@qG zGHZQAgF&PZ19J23_kl>+3KdYUUX3=WHVWqd`AR6X7W&fnXQ%Pem=xLIyi(C?*KQw| zhq;;s^k$y0o|I^gCC|5yq@|@ zpu4YQhRUq(X&--HqslL2d7@c4IKmy+GICz^FNC4PA|hBx@S+S%=sklqfZ(UHuI8S0 zJ6wyX%K%u;IDj>h2mWP{i(qPUvL{pD2%+710N&~PHgOTy4+bCsB8PhJL9Rfi;n{%M zItkXZ>1oWnoZSE)ho9{cpmKbUSD##= q-spS`Tm{n>u9<$Lr zFFS07XD>ysrr%$nZxrC`rv9jRFz1i3j8=g_wr5a?d yp*a;D=!M_B ztiSV{H>K!9+ l0iKr)}bSUyFhDMS>>V4jt&ph^7?XepytN* zKn<0;i7a}@ m=vl)D=*S}|C7@tO> 0zxi zu#1W5f&gv1NyBIQ79<9T&jK)Kl!8{VAgO6AZ!&fv!GRQNN%#V<|9RKI@jWA#`j1%= z!dZJyqgBwC?%%&3Qff7v1+un9Z niDrVT}`*Kc2{a!wVgpdqU}Npv>0D*AoFP1)=lr#) u0)&NUA?$o_YmEcvHTb0p<)T1 zIO2QU7IIlhL^E@8a!!YiK;&e_41#Es(t90KBeByOJK|==XuOVg`}-k@cqclwNyx|) zL@wL)eOcI)$Aa1lk>3zE4eJ0BAK;hb24gT@s3nMdj%O7C%xZShnwf9_qzOr3Xi3qL znZdcKKwF0+0wAA?=3O;k1G%g})VZl^9fik;6!2L1KK1pciB%?ccJ`OOe?RyYs#yK0 zBve3=7Zs-7GEtvT;6fb$LFEWAkurit+(R?Z?L{PoqMt|72OZy}H^qDvY)+>$>c`&J z8^qne%+e1iwr0!GsKD^ 08>8p{h z-*y#>E}QA-shOEZoy!EcZQQa8tRr58X~r$#NEYju0VKH)^eT<^;}PAPlr@Ch%hy z@x?l`V8Es$Q0m$2H*P%m^pI dGhDF8g?kF-4l2Ss{aInVMe5G zzXHazk$SY~-^i)m9}yNtYX1W$1x1&%xmTLZkpY0p!*f`h B8kmEG3?0f^*{m_I-qP)=8jbC?R# z$nriqEb-b2mDAU(zj);Ri2OVz>eyU^Le0NK6<-GWALf&Vld~UmD2Ls-bPbClt6>pL z*YyaUu7N3>57=f6+vzF(-689N4>!RgU{kU)W+OD_Xbt)haQfb1jtHoU5wz|DT94mx zz9-V?I>fHN)!);v?1Z}bW^?LNE01GX23dUPef)YZhc^CxIV@~)4NF&$qDR~o*f9Zn z=YFFWuVRPS1{_*qWxT!rXTXqVqxbnE1-j*3earpn>!c*s@{YpDN>G9l^r*yPP8r+^ z%AkDbl|K36f02jaHezI0{86N)r^kU2vH9i0a5p*?WOR@lGVaC*1nLT)A}gJj^F|N1 z$|$)W`e>I~<${)$nj {2+Dpj8i8W7>8m8iO@PYdy0! z2>t3oZ~+Rcw{Iv}8XqwOgfiUif_vq6-&$m80SJ}9_q&QSD&_<#{K6}#Vz6Rhx z)yCyv*J^sF=~lhdy*;-zHg`u}kf OGISafhB zanaqWP9Ey>G5N%R0BhW{^ saFFFzDCwp2#CLV;WJ{Z!E*Sec* VgHe9Y20r z=@mYai0tuxDF?QW=ohydy*OC>F2CByeKoWcRM9$S^M${*D7Ac#|D7pQSaW5q@Jf>S zkXU}7Q-YV^1F{6$h}wvus*71!QKVZMZ9%8q)U(?H_HyjlL*!DSaAp1mL+is2&*q=$ zkYIqD7;?e;p@1|3OC83uITaLGX?-W`DzbDwMahj$Ak~C+3Cr}___!f}(%pbE{o74{ zY1j&D$_GN1ln)!7`zt)*&xQ9^3XjM$D-1i1nb~f;gp_wCidwm)0BnIb-0`Rdf}xp- zRb}Js1p-_()nL!?hzJe$YcLFaW J% e{8} zC!nY-uK!}bo>CU2&;< (E%5O3;PWDsFAPBLwHS}Kx-{gyMh4(kV8rSMZ!CrV^; zLooeb$M4Y3Fngwhe+|8fIH|(C7O_x#1pY7d?fy?9@c)Aj WNGgPC>OX-q6Q#YE)xuS=MYCrQh2jJA1PB%_UO@o(kaEZj&;UpMqjDX* zBWDqI_BQN)F49du&r|x)j^^bZP=*?ihvz61>REo*&Qa}N&C!lgTfpY%QlU~W9oNhM zNxaClL>nr~(q9DZS)SCKCHPm+N$5xi9#l2LD687KrMa1{qytC;AY<6nR9=eOH=CCr znQnPI`!CPO^Zcg0P8ARY`eIvJAXqE~qCPVG|Dz^L#E~MwIG%%*{m)9ivPJ4`u1Lpn zAo4N{7IyJpbF?Hj?DXrBBJdJynz-vh{ibLj4a;Bm*h3B=pE~Fdk6C_y%%u4xBXVCb z)Uc^V6b+au-k?6O&^v9Xwy(0`RudqVW1s!bfPA{RaU!JtmMmnhNxKTRO<@t;MyI#E z2yZGmp}oCCjQ7H0dLg-sX4i03$I;AF1uy%Ak_H4G#H*>mS(6G3(LA6M?xA+DO+gi6 z`DH%G04nGCp3G3R8R*rOzn -DQ?=%#_EQ4+Y ?^e{Knk!X>phmTvP6SR^F}-N3ArKuK{{EMDdhRYS>}1N3 zq^|tUrr%@$NOs!l=idq5u`($F+>n?D)tlGKbGq51!{wv-9&r*pejDjN@P&gh3tN~6 zW;*x}BBZ;~mGhSacz %?k? z>-scsEyX>gi6L@%sR7xNL^L93;T-!R_5X@Om*-?x8}v?p{wXnFew7zeTtAmIc-&>m z4EN|rf>e*LjcC-~1JLunMPy(6v>HxEHe&Z|A7k_hJqcv8WZ}E;=unx=dxc7oP&iUH zw%=ALNwLASPEk2aWY-^BQ|4bmZEj7+SVwo$d2WOt1>q~HtbMZOuL@3GP)gc+qwV%? zbdp2pyK+-TE2`R>Q+4ojMdM?y-B&z%o|yKNd&(N*5FAyLY$RLC68_A$E0yc8mfo{U z4M2UF0CUV+^d{I4K{!{}*48c&SQW23DQUjKc^P%-;i>0q?QzWgCNZfK9Df{s4_{lh z@D4SSxC)W|?x(u#+MHXXecD5wo+}BvrKYaeR8L)@g&l}Zg}j?yehJld3j+H&n|%H2 zGb3W$5mZxCYaCjwOTbwV#7byp@ZX@zx8Iu@^46aaN c@GN8IEHf=G~!*yFX!NQu(Jp{S}yMC>XQNx@2V9@IUbWTVot^+0^ zw@^ >|dR0!jh?Rt;rEYhg^<}-m z`sU}C_}~!D4(zW0`v*oD$Pdt@M;>k-Zlp`}jvgIu6PKPlb)s_2mkMFJ5z`@fYn}`} z2&J+rzX@a+cZ4?>XAsxQzL)nk^_TDcSz*a4Ui7=o^QN-Ian%zj8K)_(SdsiTqTZ6F z0hXhZCXvf}#T^t%mOEOm^(vbDHfv3-t*ZTbr21$!%@2u#YhXW;$6D1LIDw?i3}$X1 zSkwkzz}*J40OR~U)*F0dQTO=jz^H-7Ksd_T+Hi)3kaX?-94(dW`>{O60pO8J2M7f* zBENQ5zF*Hhw_+dYuq3c?beRqfk*nuu8u%taP=dCRScDdXd4rkoTx(*i{4@T#1U!(y zzgOP3xZRDAgR6EzBqID %A zX^wmh)Lu~ +@yQeS~Kh^*&yog&ZC|t?a)TUOQ(l4eyq 6aO|M_(iql}cg! so*0AwTF=#0lGv&Dg?zBHIEMfY7%FxA-+y#J-LT6NWNQUV-htj7PiwkrwsA(PD zWT{ydkfvOfLh$AXL_r143FfRL!Ma_)@DaHB3iz&P!d^BXD$PaDt#?R7AjXOsAi%LK zR#etZd)v!`=^fxa)E?3YVrSX`5K)*hWBpcN>YDg$XHxySAxhsQe2((Or>wVB5PcJK zL3I7e>3nj-+8`a9kY>)>J0T+HvFG9Q6RnabbY&9UI4}U?!_m;>(({CjMq(YRcc`R? zO78-3f$ah3t4wy8K%+#-feYir7OjMyndK!4YtjsSS_J}7Pb-G4Du#cT?+RG-u{HIg z+b33gcFv^)CfTvm2Fr3pI4aU5muUVKVtoz1Ejjh?HmV!hy-#Cyt54{jEB+~d5~bM! zr2z*9qTBB7LO#+abO5%>a5TsRZBGQB3~(XcFo !v1NnI>%C G9z;i+AKn^`yQtqIErTXJ|_1Ce@sbUv5inGif{LkCJJRZ2kvd* znkPSa*2l!b8>#~9pdtsTVQgx!#G$5mLWhTBXGy=_P@=v-g43Yc^aiXZQ@u}G#L%$? z-9Yl#N~b053D@x6vN8s6ZD$D|AjEft@}>(gc3~Ye7H O$?R2UZo* zScjQ7=je6NdAfgXitWy(OmBc{!|%E402vFDRX*hoI1pjkp^?i0L+bqICkU0?{U)1B zxlPxCMGY0VK~>7HUJC$+E=%S(&Ru G!^+ 7LsJ5ui-F B%$mnd$LepjBrFPHBb0XIYG|x`2_A9<)PMkObnO9zwif zV7pte_TGy@+KyanM&9-Lzklj?1?PBHLlrleXG*k6%++f8_7FK`uwlv$EYm$^o`@I@ zdNxA+y7k&t`wot{Af%__-bXG1#GMHs C(qn&b7mcUyuQKM-I#RXE7$f8zz )IW}8~dgtz;mXIrfq8AR6HQl4e4CGoF%mydk&5F;EYU6#voIocH zR^0Heyp^9>zD%Kwm*hrVUe-Gxbc>$eJnKu(ZJV8)IyIt^)|W$qREbNQXXWX xQihVK`UpQ$v$ zgdrcsmUgMKUDv|$y~V*5^u2uCv&dWhqj}&HzoB2chT6m5*wjTHGNGNrn%1$Q#Ld29 zy*|iW;UZs=Sw}N1(cr+9)UbYxJ)HNRaU>d-*t>CBzpMU}_1*sX4#WRhPm8xkf0E5> z)F-(=`QUM9{yY9xm5~>}p9*Z@zd8O?lE2US&Fd7-3V!)sQ9fPb KdjG; zkX=AtR|DuvpZqZFwS`R__@l@q1*kk5d;8~10RtZlo~zN-U%Dpq#)m23TT@e}Qy$1m z7o;!t=j)f d_$2V-kP`4_zRNb{;{pUxL~-*lbA*3|RLXZISyHO>+ezLZb9rKw4-o~QE)I+crd z6GUAv0n;lBn%uAHX jkKcbP>7$g?>^&8M1zUi8lmUfdw8V2M93ZyCa zV9o$~#DTFu&i~z%PqW0F6vQhI{q3oih?Jb1*2O`jWmDS4MTpkr7mkaYn<3~=ataE( zKWgm`I`PRk9zqMQY4>|VQ86)YV2P L6+(qLcS^WWRiOcZI*-+2xn9tCgagmHnV+w9hq^rC8A{Eu^$;n4A zUcPKvTF%uq%L$)?=CKdZHp9{gx$?z}7pKr@^h0QW@DB HTUwrt$6)WzkE&jOtH`u)zp|8L^w`O;_(g)oe1<|12L;5$Gc@Kc`dqYpe=Z=RfkI<^pfDI1 zaq#i*KEUw>8LBKW(r2QgN;wZ!p?Opiixz$w6Bozc$cnin6G*8xQ^0sH{tc|(6!2Wh zfkA*cw5h479H4Gs|L_$ao=Y%!GEPnc8*XAf=KKh+g8B$oT!k*_$i*9$gM6*A!o+1( z!_V~cWep6jLi6zX#ALpk;62C+tN{pu`Up#4snQ_q+Bks{5}fE!Q0EWSdfxEbTN0d{ znz||^M1D_n^%OW{PNgYprE_1$+|t&jgs#`u!1c{;x;ZS5&_G!(1!Te E-_ZlsD2_3CdGx3UfmBh`yFAYR53=>pYv!?m*#tZ%1t6-F7YO9Hi2| z!pL|A1lSkKnO~sw;yX0+2at!!6VQs1gK OTN?uwz0*M3kdhZ~}ZJU*VmdVDftcV*Mf87`$M|Lp z-WTXOxky5yn$XgjYJ+>6vWB;z4&Z{KAF%C7Sn@#CQ^0mB0A$-5&EK*}X%S||N7))L zMg?)<3s@D9B?f!9r$@ClMvw@si~_
+2>Qz z=_v5nPrjk<9q&7-DJAtHc0Y%Wk54uJz%(acFTW9l>}<3ID~vm}Pk`}{zBQ;m2vHxL z^Nre4OJ0ms*6mEks9_Uu<8C<1KuENpF@Y^E6Iuw`?FkgMkPQlYaAfw{qxrK3xlK$= z2q^i^qAm~<8*RK+t+#O7S)0 63ux}=yMWzwtxUuo^XE5g11OzTa zq2&Mibv;x#Y?RJ$%)z&B-+b!oB%sUvHMHy9yl_gV#*Gv HccBamzg&bQY}fz%3fn&XfV~UC$j}VMxjtCERU<=; F8`vPaGZ|W9CA86bn}|Z}H{4!lrxHY9YGiVdV>L?5 zdb7g5*Fx+;JuD_d{n_0;fr8@E;?0s5f*avA+z
P^pZ)x)bL!NoNdCuOUgBPdJ0_@)aC*kCnHk0ecZl!>u%fD=%~|fwoi8B8|N8yg z%7NN~ed-#$xp{Ndw$%M}tION=qZ_A5%B*^Re(T1UCa%XJHkQ9*=Ts^bKS=ImCMdAi z{qgOEs(v8*QL1RPMfMpj*NH3ktqBrzAcGdoM+SYd8TRz>xKqkK!~4|F%cI4tPKX z(MfRh>(@kNcA!eww*}SoX1p5jsUo1qVHAr#2lOlo^@dKOcF@(?+1V9$8tRoc+HGj@ zL-_!yri$@Wr!9lNX%vgoyX*AyEpqaVa#mKCN1Jh@<1uDR$FR$}rO~ 4_7Pz)IdNH5;zvN0IB6%1KF;-lzRZVSL*+CU!V;PJ(>UN{mv7r5nzF@&EcE>tiG| zq>jw40{e6l<+MCv*&7LMbf72KdI%$ppMUoD>YJJ{@vf$kVrF-ye!K!ULS~3OJYgtw ztdfzDG3H4TU_b3f?sr_}5i%(u&OTuFF{YtMEBrv?Xe%qzkLLX7qb_B7nBQy8ZPY-( zrc-_a7;+hK%kxnZ=IK;CgZ*j !oYQ;@+vAI^+ zO$9Uw{P;l{V;d{x-uEbm&U^nNNX;Da-iNN$hr0{N4J*b9eIk78?>~Q&xbMuoM*1zP zLHtUe%bh!AlCInFOz~W&jSp}KZ&b4Y{c!@NS@f6?oR`^gN5HP#^|~ntO%AA!$YCrK zb^PHxt>z8_1(Qdi_r$q>Xk0i2xB Q~0_x{!u=wuTWTrlHS(p;&LA@S2BLPLd!@^@d!kTgU?95X)GYV+qYWOVr zT9d*rsmISv`WFa~FKF&q;L5%x%wW |Y@f71+-x5&%2E(WXpQ?@R&9}gzcaf43{7z4SF&l%!iyDx#R*a6CLK#72 zp@ }c_3Sg9k9wJ^&8HXH4I+nkpb!2k zA0G+i^9!eNV45PK-*kWJrGXzmf@{pF3A`Uwl$S?7Q4t=~8e+JF>*eJ&`=f}_Ct%o_ zQ@0Y`uT|x=L}_SfIH?ue5RR&>s+xeF;LmXM(^C+59zuinK(QGye2AETvOn|$^a#~A z0ABl9>DV%PH^2XOl2W8+b}QD4(Doz?)@>^&Wi8k4M9+PrB4_FbJ=xkuHgB4ROH}-{ z04rt(3a&w#{|^3SCK5cRFnOSk&cKl)ji3lv7B0iU1XJ)nPAo4chm$Sn0pm<@O8)%$ zvorUBAOaW7O9mHRcDI* I 0kri<$&_?o0frkI-RoC4fFQHWvPbvc+gs6Ib8+fFJrB zd`7*GXl~xTx#qYG=63{YQBa&qoF>GNV1M8Xw|ni%71XmUAHRM*?R~U2ilh1jjU0f| zBP?F{rq_Ex&C6a4(r8I35g}pNa=!^hMK0ce0v3Hle-;mqXRtwutXW$n_c-nv2pyCd zNZ&Qavh;_Pq0KJP&; 8^vV=9oTsd1}_2&MLderp;%JJz8ZsD~}yUB8DegXH! zhi&~>%3v6(F@V4TJ#NDcC#kYrymaXq@X6zFgxFdto6)J8AK~1Jw^32Apf4;*Hkek4 zC2VD7MG{I(5Ht}WMnI#zH?0)b_!H_6IBDYf+qV~BQb7!)1eEj{EbbO@1Tz2xTU+n6 z#@!$T>i?y3CIXi~*iFepC3i(gXaFVysVktp!vO5ME32zd71N)qXXp0257*b91nBY^ zZens_Bt@3!OGCprFye1Df&qI(Q%^wE)CfCt723{eYHMquD7`1XgZ4Vu)bN0_DI9|W z*-SRCuhbnKR(XuJv}E^f`DXT`CQp2DG_b t0hoB}H0yq@1pC&*8 zESwO4;1|2Q_v`Ik?V*2P=FlTLB;>x)x%W^M^_SbqKnB-2QzCGE8?D{C)PLdBO`rUH z&VYb`dDZn%b-2XTpjFW282PkkuU{(!#)J*{C3tj%j2KCNZMA}q-;x8?oSSLA%|+3W z?0PQi@JdM34N`^F+t}LPgz5m$c~)GEh$G8;#W-!FBj{Z>1nyI_(jikTcpG^0o=8y# zJG;-Y x0^ulups#6k%V*T@`AmbO zWxH+1w|ec}L9-{jF(f3(uRI4(x)}mQMxteBfsrhPfs9^r=v*Qaa#|#R_WZfUe&0)k zN}>)PYgqz_9+|W3da&z=k_0M)jaENX)ezjX {!wT%l1IXRaYf5i8u&wwOqzvii8jcqnG zGMd_hwVAm-F<9tx05QTpGLrmVSeOB<|AA7=m}=sgVzX{2U0=Yo0@Cdl_CgY4VX#Sp z=jCTtet}$qWH{JWJ^>QZkO-CZ^nxHFq?8DF@6Iv4p}BDx>?`+6?4jQimepqfnUi;A z;XI=cV0_rGgZ}F)XHP&-)HVKk8bv}vas>)4GB~Sm?MG>GF$ps>^BUa)1Y{w!mxUf4 zeDtmYSRZV_C@Li*lMHSoB_@)U&CRb6dJEv3xW$p0p(KV$rM4XLgatT5491-mO088m z_wN@R20}LCvEF7*w}|=$r T^vs!7Uu|y4L50(WSp@oZr8Xe0Lb0t%i zMHG(5=rrY%larJ5@bD 9-3`FLT>JehWmorHiYp!Pw`Q{;ytfR~AS -BM=ex0sTQV_ot2ZcWv?zP`RsPzV%T zB?L%_sNH>W8r8P${I~#*6(=21Y3`+*68xjG$B(;4Mj8fi)~PM%*{PGo7#zAZDrv*l zYj@fdkeB>Mr+7IF%Xy*Craw^tKNYh@Ocd#7I`h&eW(!~}k$|$0mdwi&`^KhgCJQPg zzV1(2cY(_^Rm1&+t6D@=XR!D1uce2#%B%!s?QnMCL;Gxh>Je8_$hj?>Lw!d&D(B%d z2)3l8q;v^@=KO^V^`X{z6`Yc8p@L_&{a=1VkmWfZy+1%(@IgB2Q~A(gblr5=?FoO= z_uk%nfN(Z(b1%R_vEOA{466aEQnCbsg3QN3;9Mq?A7n_hvi)ozQAgGhGJ-5{bSPI~ zR9`VOD)`wM+X1QJp=^#S+yv%D9=;N1*~&UcJiXaNV~R((e-GJ*P7bwY8LWQP7X0*X zr(80&Wyf?&clTZFn;Y)#eoclz;#v2aJ*l)>e3Kp_QW~lA4o@<`!YX=&A9`)3mbzqX z#S~xAw9bXBCSv2sN4~+M_Rv|oG(Kz$a3}VQmm$yDWz_B4w` 4;({&v89qPM_$z#pZVz)Q#eMjwbw&)_Wq+s zHvz`CJl1tSc;UiU{0&K`t-Gp=OuP*5(*HuXJ4hP?=+O$$=Q49Sc1KS``ZAy|3rSV~ z#iaaZT>Gg`7WV8Hne_u6jYqz(j`;kiO@^;x!o&O`{5E#rUw{3zeD!LBs=_2}F%axF zVZ~vD{nnj3yX#%RqFPhRAMQ5PeRF16n|-UF;FO~P!_HrgOAbEow56N7^dm YAONzOSXV6|mt!E^h8ke0)oAHexdZk6sw)I4BcaX29KDBrwVRbYG6J-Sl`s z^MeO|w-R>fnB_&`vCO`rLWEw<=Q%q?|BL$C)8i3x*!>*18nuiM*j%|f2ZlaB(P-e- zJ?*+NH*yoNpRm*{qA`05?|1oA;x=7?Yw0PJqO9UNYabPEj`5PT+P5XW<7#~G9pO65 z?}I0XP1PM<{G*A;_E=<_6sLX1CZkjhDvI_q&M`}%%Ug_lcST3v*Set}{e#oFu`}*{ zJ9Mu_ng%mJ+ib%6Pa##6$lV*I-uC|h$ET<1)mn_LRdkPZe90D-_GG^e9^Ap@H#8>t z@R%TT!XqFcEhR-k1@<0*z)+<2Hcl~r?r5qBP4??;rJZj>jyP2H{1VA#SXl0IBl4BA z-mBoKB4w>?$5Ho^wY@73OG`7LA`GePZqp~HxV(HJselnSGr&U=iU5ucF|=g{1_mkz z0)2h&AiW*#ZCDEii`yYZr#*u4YIm`Mw(#-A77y^LCwv8k@6$`GZC|doGLt+Zlti#7 z WBLz?*s5I$i{+L3!_f_wJGe-CWO4pH3T`Kfg&pU@0V| zl8hQ_V%4{7*)k7!+!Ik-0VpIbJh=SQNS7+Vb|wiq*r@X`$$dJ$AJisnkp+ORNEVJ3 z!%c7RRXEF&97Zjm1-gKMw~!k00F0o~tweea`XhBgv`KH;ty^)?Ws54125MUSM{OG0 z*R9z5By%#xvAG~sSi^qU;{E*&iQv(nH?r~?OJ2R^bOBZf%z$;*S&Aqrja>!VGznus zC!g%1!zB@zKuXH>jQ)@I_Evm528#d8i}_@GzJy-sa# q%5y5dLVl(s}zuXU?;Y z8SV^UIep58FZmEwi-D;Yz5Eda?(shHZ!3TCP{Bu@(G{1h?@2<_2e~pu;a(?R`O?(1 zL0)`hHNdNTI8@F-&GW)NhrY%Gv&m-+AY-ovuvfUYZEM9W$6D^aT9}Y9kNN`;Om0 zKX_2yDq=24O$ z+PZ8Nd)(Z89uwr88}(WOdPTQtGZ7(><)1_hxFk-Rs%d z?qFOs1ZhAmhwrWTNe2$l_4TprTheV-s}7_TLz@%@BtWqKN98d-S~38yG$Snve8Kpl zz0l%hpisGl6 2yNwauTW4siQKE>j<<|!>8fSKCw@XIXK3c|P7HZl@S44{2*RVPK^ ztU{hwfH`c^HA3N+tWS=8yk6X>d7Nk4Q9*~=4cYdj;&L|UEz38DPROAP 1yrTjG18uqtWfKw-kPW5Tvc4$pXW6%J->p!iLZg`7 zZ+_fLe0J^G+pUcBT3T9PzklEV=>B{|!;YEEL>1l*eZB<1R&2d6BA-((F%joa?cnFX z0eQPr&c6V%ht!zRQ3AsJ*@^ydK6Z~>;Fg-2oU8!4#4dV;sj{*X`OWL=%On`74Cg5k zh)XZ-Oq3HWma7( 7QqIY``)U-slwXtu1Y8b4-TWasr^IB2g z>P=l`-U@xIlsd&kmEuJTHxJQ;&(qnn`9+}S_Pz|O;+o$@xxeZ^KBCsv?n`vmp&U5` zuGJmpPFn(?(w9AUn-rwE|H4s0<$+WC6Qyxe{j!5cj*Ra5@O^hRo$$Nolk8IQdj$CI zdV0QtT;6reJ7s|ds9FXXj|4Q*xQjI%ku-3g95_ajNQ5MmklPXbaL-5w2CSN|P< {$=0?I0LK6Z+)I*0 znM>N3k^3Y;4r8u)(A}+h{zGJtGaFi#$g#*$LFHk%n*m@Rv=3{)+o`Loqnr11%AYsZ zR8n(bVBnyP%p$x;?}j@c0N~)5p?JH}Yabe2G$Km}^RpRsD>5_%@+;3lQ05TiRFKeV zl@x!^9|nEk4pUUX-&}y3YqXLbWUZDzXkEgpcEh#su1BEg#7@cM_IX9jMr)m`ep4S@ zR<3w(GcJ%SEh>5q?DsaPsGHI9xpwW^Rh(>skYi4P?rp8TvT#izS~{&vMQ{IIz00Dj za^}RpuXaH~UPx^e5=to^KwE`ADr%m_fw426^}WR1(1-^1u#;pXXZC3 fk}A3LZ9uIuI@x#`>UJXX|x9W#5UxB14$?sV#Svifn!EM-w|8)g0aq|W13TwEGg z#98U}*vfm(h|{sLvhvi0?jcJ6>)J$Z#zB%=wL-BsD?1x>4OC1R)xNiNiL$aX9R(Tm zr25t_YY`IBsftA1T;6u3nKjiZ>jJ;|Z@PVSE@T^+|5$qs-_@O)kv`Dw$2v+11y|}2 zaOF?QhRa; h0F!$R?yG*6-cBn29M)Q1K5nNkv^P zh)$zQAR9-#DE2Nt#4kOuNC^Z8c!&Ln;u|1#q>m5gGwI{sE-D7gb{i?B5mcT?x;Q@~ zg&m1bz_!2#WuV@dmmfQKZZm=a^TFl2K$b9BbP_Y8);QKs7Zf0m8~~Y@TpU0#%%Z LC1{HV(-yA;gh5Z TQy-;N>_3S<6qx{$POu oI=LiK`wEuU5zOZuT%3XT-^JB%` zxIqBC78b^(m31|#ep(P?<;Tzv7tqZy8@%Z4rlK5Yz=;Th6+ZKj<^z;>g0!gp31JpmKz}+~ug%8Gg9A<8FQJ-O!v)eSD zK2z}}r26SB?I+bl+X<-~s0X0zq}S_zz6)$It=;3`9Op-YTr(k92@3*2WL&oF*ww2U zpPir&;igyMwoV-6`U+4shmk6t8rUx?s*m)VRm9pTV0@MEiN>l2x1vtXU$C*0Zo!OF zID>+U&-E19k`VVbn~(vgYc6zf|GrB}S^0sV-(9Sa)>r*G*H46>D|I%hj+-toMxQ@; zBKA}TP-wGSH1Th}SL~dL*}CcTMcRBsqIY6Tgk2A>U3(~2eIB(maE7tmdZNB;#Q{&8 z*i;_+qu-sVTWot``{_Bcw1L~6L!G_1_b@iu_o(kmxBVAke|*1xd9l~m@@~zxpz@~r zyLatmu6RUUz%d48;Dbg+JSw~_`T6+-js;rsv}OkKpHC#5M1X)(h@7Ol#FR=Q) i|a>Mgjh4%yEE(kq&1jY{r?Dp9@9OVzu zjm2;DfDByQ4n;lnlgdeHT+kO_`};IUsi-j>hX(A9+@J7OEbo0bV&aJ>IG72@;bd zRC|HzLhf47Aqx)?Zv%McE2cRJM+qM~PhMUgH5^Bw7w6-KFF8)@IURpqb98dY0|ZaS zUOZ@G-P$$FS15YhhNdtraETY9s43mg(Ms7IxocwETAfdDXJ+^&7|$Lh_mGDVl)QXm zQn)N|f?6iwNG*Pyq{NL>2M610bltz}0|kGJLD>PAF6aS*+&XcBb_00lW&weO;sLN? z!y+R|M;|aTP_@ JT?i=ZIa83vm=7cRgiGICbRlxkIA)e`(l%mlLNMF0!nIR@yo_T z+~rJ;@$~q{hq0Av2Xcf)BLcRdhrM~TX36_0hPuAV#QFivuB!bHhK?32R3SKk|9(q) z6yN*p2DvcLh2|B()gXkbE^DEdh8f&9Z{LQXjdol{h6GyVXb|cpIGK?U8(~S{OHr dESg>=p7&r_!Nh>oN) zkf>hPMI=MSOPFEpkaPY(&(*sVo!lD|zH 1|gj(x?%dek5^X{ s3o8@*_$^Mp z@Ao2+p_5J&_H7HKH3;5>^sj4re*fW-6RsD|pI7xerDYd0fQ2cSGlnN@(U@@!^q=BF zNcQRkT_Uj!uMP=?)AsAyEoO{HENy6O)(*MEC^xI=cxfD9+qgYCRKWNjten+)dV2fe z)>4@RXpqc)( j&wC<1F*_>V5E^n2@sfMM>s2&uU+eoCAx>sB*dOu>J_-^`n?RY$=^t`!zur7 ztjA `i8;$tuU>N#U~ zJ4^2aC}DRCtyjC*$3j+O_vas&8i})_PeBeM?KLqmF%Jr Y>N1c`$l)esk9J*qEv2% J#!j+ z;JQ+@P(8>C8LwUy)YYx!;NTcpL}-{%ueWP+@iY#*J}5nE!+qyYD}C@K;pv6$AKzHw zAc+aM>%CUXKW3-g-nJ9di5Ip>(tXv?V;ihCT>jm(xhai#QAY7ZWMlUaxz_b|Y~zp- zCJdIpYR(mZ`0xfn@@rSu9kiB0!1pEX!P$0ig$X#Q@F5UA)KYz^cW`J(F6|}s1rFal ztftV<{c}`__t9l~A>_#t1*>pd5AgSY;r`jDeQvmMOnxeuZ^8Tp6%vzI<{gv95jra3 z)T)^>++nU5v}d1f==FZ@HKp{aX^>z+Wek#%*hEjVbOSj_rM*V`v sHoyqcYabx>{L( zc$OlL3#<9ux zHxb)-0 zIvmr0rErw#=eg2tov?Rs&>v`dOSz)bJO^F|pes+u95g&2V7_h5P8|<4b)RnV)_XjP z2J*Mr&MjSYrp*%6%fI2+VcH{!Jb@xKOv}!X+!Tw=2`^4bPc*@R@TP5t@6&SXEAXH= z (@j^+5nnPQNf%gj=etKLaWz>bb?pMvk35yAe%7 z^t(ae3{i^y32AxX4R>FIg2h0cetwqjh~f>va41f2_xL$VzGBv#s3-jd(6NPw=Xg#Z z&~-P20wnizPrf|A#eTO`hZ`&s^nijf9EP%(SoaWCrFroP(8v`yu)(;5YE1(#`{vD= zUh#I%S)3PHY-`q<9>4$b#Ohdf-HMn#(F3z<@ADxvGB35fSwx~WWXo|sz$ku=#L5y) zDL}$cdiWUWfe6kv?|p@rh}#e2%b*9a5D6nHA~TS_3D|t0a0?^-ul{B)= drLjv1|wl0tRe(s zW=>Ad4w4Yno2<9hL?HYpBmjwt@$uR*!KiPHN={#xx591b96BmnyEmygzeg8I4+I*@ z73eZ5zB8i!(ybWk;lWVoE2qRg?Uvyu=h=qSuYGkTMbE&;^h?$O-~l_v*A&&3we{xP zZi|Fh;QYd`-_(3uVhtb-DFEI>9==xB$ddZ{1dCrsi$?Bju6g|+o$h99&*`4l)_L)2 ze04#V^Fw)xh3&8o6|}U>uIQp-)Vyh jcaX{xVfXSOhWmDdGc z+A@?nJlh3xUV4gy)_X8qi7qNC7W(2Q?YH`Nzqw#Cs_gJ!-$5n%gBK5|nR))cm@@H{ ziO1G#0d@R(kW}Rsr_#QX=AaSpprDJtxNfYJIn*$Gl9Eciz&jbKy 8sQvfUCFI>0qVYOH1_UHk_Sg@}+;P>leD1e-S|vs(G3+Ea$;BY!jaBwnQINg@blD zfIb?1)ae>C-sYp8Z@V~4T3cJ&CM+#2Q`cNKhfb#>sc~wTl6q_Q)pr;JNDucp?QUPQ zV#TJ>n(h&@ISJa;)5Cv}9gKSUO&}VfhxP_=nSP<+F`wu5^us`$X#MH}Z}=GB)w!Yx zR5W_AF+;ES@6-pt*GJ6d6Lb*M6@>J(dhr6Vw*{|Yxd0q_W!wLwr(F6?^slWxm`MN2 zu^G!CG{;kAwg2B2# $PNEE@8XyEB-7AgVUrt>wGm*prsu{R&QLY{jx`d)Bh9 zcyRsli?*+zWNgft9yg8C&E-Hconb%13mgQI^ctib^Gi!heU{A ZT;~U7&9s&-Te%yL$C$zNJrrUDp);xyL3w!wEn7pM}`krt5Lx$7GS* zp>KpNs`kzTYNHHIwGe j-ll<+e#Q(-O 4}BtH9DJ;T(F(3*_58DnJ$RP&dAtbO@|8Bs`SlCS!^I0u zxO+<3xu2bwf9%Nj;5Zt!LqOnN?#$H1EciKI1q_9NW^hmMPsgHz!}p6 uy7#(PjEf) z1NV ?}u zIcD3 lI{^Q%IsZuet7$8m#@hm~(vvY4%r@Tv=Cui@&L(uBx;AlzvY$G| zlDl>BiTw-yZUM6eF>LQr2RR1#>RAkL0iMtw^p1LP$>PfB*5`9m0Qh>shxW=DgIK*r zP*Cva;9x$y<(`|?2~hV ;Jv=lS09O(~b0VGrdwhen~=hEX&cpNKo`C9`L< zjhkE0OgTLw$Pec&Mo>?|)HG=a`*WG%YD>_;M2|?n+h1>b4kSMS-Jm=y1QEy*9D}h8 zq$wB5{BM(|2&lc*`P0^eAqrBf-`1~ZQ>=4N+e|2-;fLRbD>RmlF8S-LLee^)r*}gc zWoSL@rP_a8vR!)>GaaJRN&y=gN5);=tzk;sLb*H4kXmAeVpUjOH%>IPlH>PAAUnps z2E?b?o78hOw&qZ54VyUbki{VTXZL5>Th(M-0q?MjSTmwjC`7KQ+13c6jw{CH7Xntg zrc|V@VzT+pxd88zSZ`{c5!0kQ1q`Qum%+K$O$HBvV;t7mgQlE=`&rr+ot=uSH|#Cy zc%MWr1qB3`urLc+H>_e;cM*zmlu8WmDuY{WekgOxgCxmF0QfcRh%SQzh?AHJVIb|T z!`rOc =P(qL1$e9Da;{$_(e@Q ;0gheuX_! zkX>xsDZnn61|S_^iG%xQ6x-5BBq^1wX8zlyg@uJ`dUMz1+R!G_!<(hA)(|P?7jjei zde8da`%EPt3h%AoHgoTA^^=-d1|#g8B*?KqGqwO-S7f1Y(PBruyoTkwd(3zfFrli> z4>_tAWOCgy=vb6ndp8_wYsfAVbtqIGyKq|iZ2x|f%nQFf+rDUOPX`%C)g=k^yUcEm z7xDgsmp>HeSk(4vJ#FQ0TEM$?fHbqPn-ro|2&jgPIrijmbP))o1OX?=KGIY|K!9k5 z7houL<8;Gv4y-!RpjS}K*R9*0XrUafqT2E~{URqpl+dz8y0&T8uH_h-dxuzg@*S3h zsrHU^;3@b$Y~tGSEcB?sE(QurM iznXkZ1>=Sl z%%PK>+K;nnA96Q40Pws )B6O!z{C z8{1aTfdQ_COWCE&%=kgJP!11@*CDq?;3qZ;sZ9lekJ~r@eMh25(;3qs*S>x|OZZt} zt`ROBh=3sM$8x9fYG&IDBPSv(b