From 81761f7784e74be4f8857221b179f7ea3302bac1 Mon Sep 17 00:00:00 2001 From: Collin Erickson Date: Sun, 11 Feb 2024 11:44:08 -0500 Subject: [PATCH] Added kernel stuff to readme --- R/GauPro_S3.R | 32 --------- R/kernel_Matern52.R | 5 ++ R/kernel_model.R | 9 ++- R/kernel_periodic.R | 2 + README.Rmd | 59 +++++++++++++++- README.md | 64 ++++++++++++++++-- man/summary.GauPro.Rd | 64 +----------------- scratch/ToDo.md | 12 +--- .../README-diamond_construct_kernel_fit-1.png | Bin 0 -> 12452 bytes 9 files changed, 136 insertions(+), 111 deletions(-) create mode 100644 tools/README-diamond_construct_kernel_fit-1.png diff --git a/R/GauPro_S3.R b/R/GauPro_S3.R index 6ed385f..0b31b23 100644 --- a/R/GauPro_S3.R +++ b/R/GauPro_S3.R @@ -23,38 +23,6 @@ predict.GauPro <- function(object, XX, se.fit=F, covmat=F, split_speed=T, ...) { object$predict(XX=XX, se.fit=se.fit, covmat=covmat, split_speed=split_speed) } -#' if (F) { -#' # Plot is automatically dispatched, same with print and format -#' #' Plot for class GauPro -#' #' -#' #' @param x Object of class GauPro -#' #' @param ... Additional parameters -#' #' -#' #' @return Nothing -#' #' @export -#' #' -#' #' @examples -#' #' n <- 12 -#' #' x <- matrix(seq(0,1,length.out = n), ncol=1) -#' #' y <- sin(2*pi*x) + rnorm(n,0,1e-1) -#' #' gp <- GauPro(X=x, Z=y, parallel=FALSE) -#' #' if (requireNamespace("MASS", quietly = TRUE)) { -#' #' plot(gp) -#' #' } -#' #' -#' plot.GauPro <- function(x, ...) { -#' x$plot(...) -#' # if (x$D == 1) { -#' # x$cool1Dplot(...) -#' # } else if (x$D == 2) { -#' # x$plot2D(...) -#' # } else { -#' # # stop("No plot method for higher than 2 dimension") -#' # x$plotmarginal() -#' # } -#' } -#' } - #' Summary for GauPro object #' #' @param object GauPro R6 object diff --git a/R/kernel_Matern52.R b/R/kernel_Matern52.R index 515ccf5..d67f1ec 100644 --- a/R/kernel_Matern52.R +++ b/R/kernel_Matern52.R @@ -1,5 +1,10 @@ #' Matern 5/2 Kernel R6 class #' +#' +#' \eqn{k(x, y) = s2 * (1 + t1 + t1^2 / 3) * exp(-t1) } +#' where +#' \eqn{t1 = sqrt(5) * sqrt(sum(theta * (x-y)^2))} +#' #' @docType class #' @importFrom R6 R6Class #' @export diff --git a/R/kernel_model.R b/R/kernel_model.R index d18398f..b29dee8 100644 --- a/R/kernel_model.R +++ b/R/kernel_model.R @@ -495,14 +495,21 @@ GauPro_kernel_model <- R6::R6Class( # Update K, Kinv, mu_hat, and s2_hat, maybe nugget too self$K <- self$kernel$k(self$X) + diag(self$kernel$s2 * self$nug, self$N) + nugincreased <- FALSE while(T) { try.chol <- try(self$Kchol <- chol(self$K), silent = T) if (!inherits(try.chol, "try-error")) {break} - warning("Can't Cholesky, increasing nugget #7819553") + nugincreased <- TRUE + # warning("Can't Cholesky, increasing nugget #7819553") oldnug <- self$nug self$nug <- max(1e-8, 2 * self$nug) self$K <- self$K + diag(self$kernel$s2 * (self$nug - oldnug), self$N) + # cat("Increasing nugget to get invertibility from ", oldnug, ' to ', + # self$nug, "\n") + } + if (nugincreased) { + warning("Can't Cholesky, increasing nugget #7819553") cat("Increasing nugget to get invertibility from ", oldnug, ' to ', self$nug, "\n") } diff --git a/R/kernel_periodic.R b/R/kernel_periodic.R index fae6aa0..eba610e 100644 --- a/R/kernel_periodic.R +++ b/R/kernel_periodic.R @@ -4,6 +4,8 @@ #' #' \eqn{k(x, y) = s2 * exp(-sum(alpha*sin(p * (x-y))^2))} #' +#' \eqn{$k(x, y) = \sigma^2 * \exp(-\sum(\alpha_i*sin(p * (x_i-y_i))^2))} +#' #' @docType class #' @importFrom R6 R6Class #' @export diff --git a/README.Rmd b/README.Rmd index e8b57e8..e20f51f 100644 --- a/README.Rmd +++ b/README.Rmd @@ -107,6 +107,39 @@ output. plot(dm) ``` +### Constructing a kernel + +In this case, the kernel was chosen automatically by looking at which dimensions +were continuous and which were discrete, and then using a Matern 5/2 on the +continuous dimensions (1,5), and separate ordered factor kernels on the other +dimensions since those columns in the data frame are all ordinal factors. +We can construct our own kernel using products and sums of any kernels, +making sure that the continuous kernels ignore the factor dimensions. + +Suppose we want to construct a kernel for this example that uses the power +exponential kernel for the two continuous dimensions, ordered kernels on +`cut` and `color`, and a Gower kernel on `clarity`. +First we construct the power exponential kernel that ignores the 3 factor +dimensions. +Then we construct + +```{r diamond_construct_kernel} +cts_kernel <- IgnoreIndsKernel$new(k=PowerExp$new(D=2), ignoreinds = c(2,3,4)) +factor_kernel2 <- OrderedFactorKernel$new(D=5, xindex=2, nlevels=nlevels(diamonds_subset[[2]])) +factor_kernel3 <- OrderedFactorKernel$new(D=5, xindex=3, nlevels=nlevels(diamonds_subset[[3]])) +factor_kernel4 <- GowerFactorKernel$new(D=5, xindex=4, nlevels=nlevels(diamonds_subset[[4]])) + +# Multiply them +diamond_kernel <- cts_kernel * factor_kernel2 * factor_kernel3 * factor_kernel4 +``` + +Now we can pass this kernel into `gpkm` and it will use it. + +```{r diamond_construct_kernel_fit} +dm <- gpkm(price ~ carat + cut + color + clarity + depth, + diamonds_subset, kernel=diamond_kernel) +dm$plotkernel() +``` @@ -139,6 +172,27 @@ proper types, then the default kernel will properly handle it by applying the numeric kernel to the numeric inputs and the factor kernel to the factor and character inputs. +| Kernel | Code | Continuous/
discrete | Equation | Notes | +| --- | --- | --- | --- | --- | +| Gaussian | `Gaussian` | cts | | Often causes issues since it assumes infinite differentiability. Experts don't recommend using it. | +| Matern 3/2 | `Matern32` | cts | | Assumes one time differentiability. This is often too low of an assumption. | +| Matern 5/2 | `Matern52` | cts | | Assumes two time differentiability. Generally the best. | +| Exponential | `Exponential` | cts | | Equivalent to Matern 1/2. Assumes no differentiability. | +| Triangle | `Triangle` | cts | | | +| Power exponential | `PowerExp` | cts | | | +| Periodic | `Periodic` | cts | $k(x, y) = \sigma^2 * \exp(-\sum(\alpha_i*sin(p * (x_i-y_i))^2))$ | The only kernel that takes advantage of periodic data. But there is often incoherance far apart, so you will likely want to multiply by one of the standard kernels. | +| Cubic | `Cubic` | cts | | | +| Rational quadratic | `RatQuad` | cts | | | +| Latent factor kernel | `LatentFactorKernel` | factor | | This embeds each discrete value into a low dimensional space and calculates the distances in that space. This works well when there are many discrete values. | +| Ordered factor kernel | `OrderedFactorKernel` | factor | | This maintains the order of the discrete values. E.g., if there are 3 levels, it will ensure that 1 and 2 have a higher correlation than 1 and 3. This is similar to embedding into a latent space with 1 dimension and requiring the values to be kept in numerical order. | +| Factor kernel | `FactorKernel` | factor | | This fits a parameter for every pair of possible values. E.g., if there are 4 discrete values, it will fit 6 (4 choose 2) values. This doesn't scale well. When there are many discrete values, use any of the other factor kernels. | +| Gower factor kernel | `GowerFactorKernel` | factor | $k(x,y) = \begin{cases} 1, & \text{if } x=y \\ p, & \text{if } x \neq y \end{cases}$ | This is a very simple factor kernel. For the relevant dimension, the correlation will either be 1 if the value are the same, or $p$ if they are different. | +| Ignore indices | `IgnoreInds` | N/A | | Use this to create a kernel that ignores certain dimensions. Useful when you want to fit different kernel types to different dimensions or when there is a mix of continuous and discrete dimensions. | + +Factor kernels: note that these all only work on a single dimension. If there +are multiple factor dimensions in your input, then they each will be given a +different factor kernel. + ## Combining kernels Kernels can be combined by multiplying or adding them directly. @@ -146,7 +200,7 @@ Kernels can be combined by multiplying or adding them directly. The following example uses the product of a periodic and a Matern 5/2 kernel to fit periodic data. -```{r combine seed} +```{r combine seed, include=F} set.seed(99) ``` @@ -158,4 +212,7 @@ gp <- gpkm(x, y, kernel=Periodic$new(D=1)*Matern52$new(D=1), nug.min=1e-6) gp$plot() ``` +For an example of a more complex kernel being constructed, +see the diamonds section above. + diff --git a/README.md b/README.md index daf1b1e..9deef4a 100644 --- a/README.md +++ b/README.md @@ -113,6 +113,42 @@ plot(dm) ![](tools/README-plot_dm-1.png) +### Constructing a kernel + +In this case, the kernel was chosen automatically by looking at which +dimensions were continuous and which were discrete, and then using a +Matern 5/2 on the continuous dimensions (1,5), and separate ordered +factor kernels on the other dimensions since those columns in the data +frame are all ordinal factors. We can construct our own kernel using +products and sums of any kernels, making sure that the continuous +kernels ignore the factor dimensions. + +Suppose we want to construct a kernel for this example that uses the +power exponential kernel for the two continuous dimensions, ordered +kernels on `cut` and `color`, and a Gower kernel on `clarity`. First we +construct the power exponential kernel that ignores the 3 factor +dimensions. Then we construct + +``` r +cts_kernel <- IgnoreIndsKernel$new(k=PowerExp$new(D=2), ignoreinds = c(2,3,4)) +factor_kernel2 <- OrderedFactorKernel$new(D=5, xindex=2, nlevels=nlevels(diamonds_subset[[2]])) +factor_kernel3 <- OrderedFactorKernel$new(D=5, xindex=3, nlevels=nlevels(diamonds_subset[[3]])) +factor_kernel4 <- GowerFactorKernel$new(D=5, xindex=4, nlevels=nlevels(diamonds_subset[[4]])) + +# Multiply them +diamond_kernel <- cts_kernel * factor_kernel2 * factor_kernel3 * factor_kernel4 +``` + +Now we can pass this kernel into `gpkm` and it will use it. + +``` r +dm <- gpkm(price ~ carat + cut + color + clarity + depth, + diamonds_subset, kernel=diamond_kernel) +dm$plotkernel() +``` + +![](tools/README-diamond_construct_kernel_fit-1.png) + ## Using kernels A key modeling decision for Gaussian process models is the choice of @@ -134,6 +170,27 @@ default kernel will properly handle it by applying the numeric kernel to the numeric inputs and the factor kernel to the factor and character inputs. +| Kernel | Code | Continuous/
discrete | Equation | Notes | +|-----------------------|-----------------------|---------------------------|--------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------| +| Gaussian | `Gaussian` | cts | | Often causes issues since it assumes infinite differentiability. Experts don’t recommend using it. | +| Matern 3/2 | `Matern32` | cts | | Assumes one time differentiability. This is often too low of an assumption. | +| Matern 5/2 | `Matern52` | cts | | Assumes two time differentiability. Generally the best. | +| Exponential | `Exponential` | cts | | Equivalent to Matern 1/2. Assumes no differentiability. | +| Triangle | `Triangle` | cts | | | +| Power exponential | `PowerExp` | cts | | | +| Periodic | `Periodic` | cts | $k(x, y) = \sigma^2 * \exp(-\sum(\alpha_i*sin(p * (x_i-y_i))^2))$ | The only kernel that takes advantage of periodic data. But there is often incoherance far apart, so you will likely want to multiply by one of the standard kernels. | +| Cubic | `Cubic` | cts | | | +| Rational quadratic | `RatQuad` | cts | | | +| Latent factor kernel | `LatentFactorKernel` | factor | | This embeds each discrete value into a low dimensional space and calculates the distances in that space. This works well when there are many discrete values. | +| Ordered factor kernel | `OrderedFactorKernel` | factor | | This maintains the order of the discrete values. E.g., if there are 3 levels, it will ensure that 1 and 2 have a higher correlation than 1 and 3. This is similar to embedding into a latent space with 1 dimension and requiring the values to be kept in numerical order. | +| Factor kernel | `FactorKernel` | factor | | This fits a parameter for every pair of possible values. E.g., if there are 4 discrete values, it will fit 6 (4 choose 2) values. This doesn’t scale well. When there are many discrete values, use any of the other factor kernels. | +| Gower factor kernel | `GowerFactorKernel` | factor | $k(x,y) = \begin{cases} 1, & \text{if } x=y \\ p, & \text{if } x \neq y \end{cases}$ | This is a very simple factor kernel. For the relevant dimension, the correlation will either be 1 if the value are the same, or $p$ if they are different. | +| Ignore indices | `IgnoreInds` | N/A | | Use this to create a kernel that ignores certain dimensions. Useful when you want to fit different kernel types to different dimensions or when there is a mix of continuous and discrete dimensions. | + +Factor kernels: note that these all only work on a single dimension. If +there are multiple factor dimensions in your input, then they each will +be given a different factor kernel. + ## Combining kernels Kernels can be combined by multiplying or adding them directly. @@ -141,10 +198,6 @@ Kernels can be combined by multiplying or adding them directly. The following example uses the product of a periodic and a Matern 5/2 kernel to fit periodic data. -``` r -set.seed(99) -``` - ``` r x <- 1:20 y <- sin(x) + .1*x^1.3 @@ -154,3 +207,6 @@ gp$plot() ``` ![](tools/README-combine_periodic-1.png) + +For an example of a more complex kernel being constructed, see the +diamonds section above. diff --git a/man/summary.GauPro.Rd b/man/summary.GauPro.Rd index 9aef4a0..6f687e8 100644 --- a/man/summary.GauPro.Rd +++ b/man/summary.GauPro.Rd @@ -2,38 +2,7 @@ % Please edit documentation in R/GauPro_S3.R \name{summary.GauPro} \alias{summary.GauPro} -\title{if (F) { - # Plot is automatically dispatched, same with print and format - #' Plot for class GauPro - #' - #' @param x Object of class GauPro - #' @param ... Additional parameters - #' - #' @return Nothing - #' @export - #' - #' @examples - #' n <- 12 - #' x <- matrix(seq(0,1,length.out = n), ncol=1) - #' y <- sin(2*pi*x) + rnorm(n,0,1e-1) - #' gp <- GauPro(X=x, Z=y, parallel=FALSE) - #' if (requireNamespace("MASS", quietly = TRUE)) { - #' plot(gp) - #' } - #' - plot.GauPro <- function(x, ...) { - x$plot(...) - # if (x$D == 1) { - # x$cool1Dplot(...) - # } else if (x$D == 2) { - # x$plot2D(...) - # } else { - # # stop("No plot method for higher than 2 dimension") - # x$plotmarginal() - # } - } -} -Summary for GauPro object} +\title{Summary for GauPro object} \usage{ \method{summary}{GauPro}(object, ...) } @@ -46,36 +15,5 @@ Summary for GauPro object} Summary } \description{ -if (F) { - # Plot is automatically dispatched, same with print and format - #' Plot for class GauPro - #' - #' @param x Object of class GauPro - #' @param ... Additional parameters - #' - #' @return Nothing - #' @export - #' - #' @examples - #' n <- 12 - #' x <- matrix(seq(0,1,length.out = n), ncol=1) - #' y <- sin(2*pi*x) + rnorm(n,0,1e-1) - #' gp <- GauPro(X=x, Z=y, parallel=FALSE) - #' if (requireNamespace("MASS", quietly = TRUE)) { - #' plot(gp) - #' } - #' - plot.GauPro <- function(x, ...) { - x$plot(...) - # if (x$D == 1) { - # x$cool1Dplot(...) - # } else if (x$D == 2) { - # x$plot2D(...) - # } else { - # # stop("No plot method for higher than 2 dimension") - # x$plotmarginal() - # } - } -} Summary for GauPro object } diff --git a/scratch/ToDo.md b/scratch/ToDo.md index 77c69fa..06c6065 100644 --- a/scratch/ToDo.md +++ b/scratch/ToDo.md @@ -1,14 +1,12 @@ # GauPro to do -* readme should have list of all kernels - -* readme example of combining kernels - * EI/CorEI/AugEI/qEI: nopt, test, doc * Use t-dist * Add documentation for kernels, esp. factor ones +* Add readme/documentation for trends + * Speed up triangle, ratquad, periodic, powerexp k/grad * Knowledge gradient: multiple starts for optim @@ -40,10 +38,4 @@ diff(range(Z))^2. Or else give message to normalize. * Standardize for X. -* plotmarginal w/ factors needs level names - * factor kernel: clean up jitter/start - -* Fix GHA/Linux bug with cubic - -* increasing nugget to get invertibility, only print last diff --git a/tools/README-diamond_construct_kernel_fit-1.png b/tools/README-diamond_construct_kernel_fit-1.png new file mode 100644 index 0000000000000000000000000000000000000000..400df4d4f431058787fb021b9c124e54d0142bbc GIT binary patch literal 12452 zcmbVzc|27A_poa0TlOuKeMv|tYeUGsMwBf=l4Q?XBZMrG>1+1}@z_bp6MM~#wpCDi=5{HV0Cr%(DMaAKmd?k1n7wXHZtHECBRDq2+;vI z7y)r+K#~o(!vV-%1r)gfC=YOtA5aqnG=u@|8-R`|pnnrEx&;_Z0j9TsN3wvq9AKpY zSlj>xz(EyoR0o_iz~uohx`3-b;BE+b8UfxWfR8EQ_Xr3u2Ldetm^JXu1_*u% zgxUjPjzGk7Ai@cVbOxec0v}z0ICtQaC-B)DNbmuY{eTpIATl$F}&vgn*eqP zIM@S@_JN}V;P?008mf>UPPdL007ur=_CgAi`3=5sT%;Wmj~R}q;;Vrl)FH1>Xq$!kHGJTRSQY+Ccg$LC4W8r zb!REnO~?LDdf%OtUi%ZRrF)|a#4<~D6}?x!0Y$49$?}!v9t`;g6T{#5Ol?ZH@{wh3 zuh?8OGW$xH$=RFX@6zate8o@xxs;3K1^9~I=qZtqK)pV3VFS@N+Q z<{56@(Zn^vA}%$u;`@QoLPV#fG{Z;voG+8AkW|1(x5-uz^i`A-4K z3W6UqI$1;xE#;hi#MxAJ+)q`zs%2`B_wBi5%q3X{sJu5sz83q$6sGWcu+GT3Ft9Rc zL&oqks$gx5X}(!UMUOgB^R#`p^lL`_HZRVmVmgU(-=`>+qlvFAs;9Zrg3)TE$`w0t z|BU}L5D~0M(oN#CpWLE_`+{9xY04cUI}UQkWpWv&b<0Mg(%3u9Lk$vsVQ7sYytEHs zQcY>`NmRo?+o>tM(saMx6~`nBdMh#Kb{e*41y3o030@<=#U+< zu`QQ%o$x;%Vjjy{?F+WXumA|f`c;1v&v%Na4DNZ#Ld#uuR|e9;KRnQZ1@M4f4MrZ$ zmWy`?DdRuNP8$9Eh%V67b(B?_#`%B)=)nQyZ=E1Z~{a1__NQ z<2QvqH6B&4fLnMBl8mOW4RL6nSn+6Rl$yZK$7;ibN!0ow3hT1!M2KT%9)srcI-F?J z9?}ZgjUZ4t(nJ1XFNbS;)RC7t$(7OSFMHOFCyna)#9zu}R%h>VOM~-`eh?EL-Ug>0 zKiBWGl#VWSSt?Cd8>wUa#jvfm3a)G^ZH1zv}4ck-!BM75rU7B0hG5iQ5)xgV) zYzA$J?zdLmz7E+x=ER2A6m$3l!Om8&(=G1@g9_OBOLwY6Rtp9T2!~rk6$wK;V;X%u zO18KTV^xr(WXA=dSoS8@y@MGMacg^nb=y=C@UjdZs7TqQwMxK2VCN1{*2daDJ80z~ zyWC`XT382nR>ZMwijk~rNmN17;^3KafID-nodR@qF{9Hr+_zmf|H3=PqB`=GZHcn- zm%pYvR0@Bc6K&rNKd_u;tQ@7Zz$VJp5(m`j8JEhCt$e>;_Pi-`W^-ZB5JE|QoFs<9 zG1U;haEBVCEtFV#cWdz3VhiKoX;)SS^3i`GA^t+@4>M%QNzrVqMV6v*vzOe9K5ZJ8 zPQF-u&+7eNixB6;Y_>KI8=IJp$Y3h1*3{Q9onH}+9+ldH(`l@`22Km#&A}tQgGV^6 zfvl*lc`EhO{UJ&Qy(EuuCh3XOIU?242|U*u&_2zon!1l4T{(5VcJD z3!sv$Ebr@!fv@O3GIy=AiI0P45Q9VCZVEPT5LsIMc_8kC+pj^iJ|UzWCc0)YH@V}R z5%cCnb(U8ymkCbtQPBbqW~tSvJK-e4b#I}QydAo|yR7y|mb8&rmzyK8mVHxYQ9y2L^sIuRM4>}FCP5L{k-#3;!hLjz8lM$*9 zJ=s{%5H-#6&F*to*TWk$&W@ zlovK`x#z+Ze)tl<2Q5#iG5>9v>sf#(lLcdbmb+UwG2a=h%LkvG=?{s<+cidx&zEM2 zD)DUEgqHtQM~a5@7p!DikeZ!d=OcU;pG8U@zGCuDD+kLY$QdI;>Ufe=1h>5S#c-W= zuQI(T*Mn1L{yJw{zp@%V{gn%-rgdZQLLw>jIi%qJf_-_6wLIsY^)WjBeQB3Vkjzn%lRxB#>^CK8#BGHzCe6qv(hj6#OqtVFf41Yd3Jz&n62ysB zJB0Y5O|&7|qe#PRN+C+@>{(R~i3h2EFf*&Vn#=~}l+&iNdTRBp$4z7*i=)jj zrI3s0+0Mi4?AeeqOCQ5DG?$GeC?bAfNbFYMJ*3F^7vxE`$~j~rZ~_vRmHY52#hkVV zC3W9**iEwI@8CEhCVGi>B}*6)am|5F-1HP93x3%$%*pnx&@wipvUt(80i-G%kGv#qns;t2gA*jPe+0m?^9i)^DAgncYHs{5^Lx$ zN^-jEedyJ|U@}CLN%i@=8R&TFwCdrjyAvtVt1Ehg&vXtXJGA&Be*L~jNcnqiMwjr( zBjR894Wfx9Hwa#M*66%pTg$#m!}j~)Y%}GHagUwTLC)#uW^-c8&+wjyco#c$QJH*p zSb1O}U?Q@cVjQmKRHo*NF^2i@6(EMyJnwrv*bU`N-p;Gy@_lu}c)92f>B^lW*vsHB z`brFqfx%>uR7j)#L1qu3`OV#I>-;p0^j&#AhIG|}p2f8S#WfPK)-CM*Zmh!q;Hc!EPuB?CW@X4al#k9 zS@aU8YN1BHfB47mNt_wkes!hYtFtbg_t`3QM2S^gRe{c_9@0YT825aw+zTbWAC9s# z-bsE)3o{$^YmfTKoZ+bId9U}~yl!OVmrL$U(v@BgaqKTcF{KfWOsaXwZ`t;|r6c>> z96J~)(-(yUK8+rt9GxDoY2C-;Dw|1z!4oTv8;Lx_S&mhsS%;amzuNr?*nZ(&><2v9 zda)-y>WirEHbQQO_7~V~iOD-A$|omU&(j+9Wucy4Jvu}^IrKgK7g!Knk>HDAHQDdh zON+FXns!nTJyOHD6&|oBh|$F=GO!3cc`seXP`6-E#gBShgLpQfx}V_Yc|j4{4|e^+ z;5gD_Wigo{3XWVg$_iDw3T8{a5K>8QLtO{Do2OO1F zor)0#rLY@E7Mzvb*|lAjPl72ouRw4#ph*dmtYn0hs6i7xY*vra$b|6i`~@HYvc8{C zWKvbNZ+A%0V>$2=S|1r9U4O?78b|?}ZWn9=bDUc9>{D&d%4<1KZ`KP=olD2TUUFaU zCL>8U=^>Qcb7GZuA?P5Xz52LXUq8|a{#v%?RU0|bAD3{e9>P96?y&>8ew{)J!kZtkJzO8~ z7}iq3^VNm>otIXyTIS|$fls1cT|-WbS22AED`h>{GE*rIH8JPmn(%v%VkVg7TD?u- z^O#NNIfdhYpvmE(B?fE{)yH|7Z0f+IHV=~v-qWSF@_RY)eus37{}0E#%~qcS zj;~#a(tUEV%3ZcEeBFxgTwUI6=@gt@bCjisivIAxsV41%fCJ`h08NNz05=J1?WYpI z*3e@wb}dRueKtja@}_b5jAv{aya}Q#+cZO{BKP6zY4>))%r(a^c{R!%LQ>;8w2XTR z)(%V3mGEYAHLFKUTbYTGX$6|UN#Kj*;c#|B90T*%5#@COPd}B;t<0K8_JTj#tFwz2 z)c%=mf}i~rez?tGvX0n1t)t4F?JRvB?H_tu>|gXmAo^>B9+ma(PH5KCsBf=?zk9qs zcu`+5rx{#55fwcl)n8#|50Wy12R7bJDT`uc5w05M{>AIzPPCKo!n#`fw-IvcO+lAnl=pkn`K=D|o`Ut=S+f2w@}?HrufFC{M3!K_MutL-HMVndME-)tKyhSImc$b+9X z9paQ+L7EovU%w}9?aI`_31yL^cyka5e@BSp2Gz@2re|4qquaKw;Xc=+N~dq%FtyJO z94l}Q#Pg2&Jl8Q9&B$Kr!S#nAv?^lfhg1udRJaLaJ#zCa1`ydYa&*eMlq0O0a>_%Y=hvQGmq)q zLk37~EZn7|+j|M75;V}0(VB!y+FHn_;LLT~8XY$!=Tml&Krkg|3Zzy(&$`jnD#C6nEicJ{XOM?YXLdVaqHTyq*u3*P0 zo@BCBD)!(TONu(k@$2}BeG-=tF;ydR)3=h3`l}*kT`H2Uy?)Mh*mPkp0U8g-`u)!Q z($VhHa3;;1=Gb93kn#IsiUOgeNNP$P%DVr~h?9UeT56-YaKd=)?J}Hj$Wpkz6N=EgPt6sFdjU>@0&JFU*LL95th$G_RuacN}-qNBRcg7R5 zsUpUZj-u?}MR0g7{kb%VAT`rI2sNX7{qyhuzW19vwpl*^YRx%(I}no4hkNo;gJ3a$ zFjImnKjzjy*Dmc|!{(I}HlFq?E#kt^AvqeQp6e~wV7&>Z?`&LQ7891&QRv5#O>?u0 z$EQE2&E#>G`mtZe1U}3kR8?YV(!prAO|Wb|U>r$U|O}!awif+y+lu z7S0T*SF{{XV@0_fE*e2iiXV;S)gqn1ch@jo0_;`zSz6&uw`y7=z(DT$f+ho-q~z)l zF6o_(BW$=diAq{Jjk7jT?=$UjLQkKuZ%b;KB?8wuBA=*on~y&~hUuTB|L?msqi7W! zj5nPt9|~Defjm|H`%AdzP82U%&`hgKgVK`}W$?BN%Oi2rYv9S9QsZKOI4cNaHDV(d zynA%s{^xNt75a+_X?a~P*YAtFNh2??VyO;MxeB9tGbbDd zC%K2b;Nr(iL;sdUh9Ihx2gXH0dhE;zw1(Q;Tuc;9N@C`LpIxe(I|}3rCHWai@-zIO zUu^++RK*AR)TCK4L#XNOGg1w53&~OXMB`7KgOP5mC+!iN2&i zY5UCd9Pio-VANXLICLyOroUcfc|Ma0p+wO3?elW5fK?6F*)LX4yb^Ho0-ACQ9H^~+ zvGUMg7JB$~heW>(lYhJq`-W$T-*v6EnQu&=@LZmmvwA|E&;Y-_O%J8yq;5TCPaR#4 z+IiFaktAJMQ!ae+5O4MhT*prw11JrS!YFI=OaanMJ<4V{ePPct+Jv#_Q!ex)efqt zQ{zD*z7l#ojs;_UL+h;QJdOVKLK}NeP&Svr9*BoIzGo1m*9*6$v)XKWKH;BcYp&@1 zdFKWD5^#-!jl{t>)ZJ)rwb?hft%Jlg9J}S>+iB_RdBAi_@2(;p-^#?g3+ajEsS zl^5v>pDC31|1^@fzMlRHc^a07;$cys=G(Bc_q7roC86OH2IZk#oKzpVSWo4%j-%uR z7%A7stB?E=g^pYU^`jV1y)$nM1SsH0CKDVV`<|NNH1X6*7{Pk+v^Ei_=2=Wky(>$X zT6d>Iw>Do9_8{NVA_Os!I36m2MS-KU=;5WlWyP&uBk~U|{jID5-h6|o47b;*8{=6= z4^8pYza+bGinF*e`w)J%C!HWLb?*UF^~(lSl>IxBZZetsa?hdc(%3ts$L(S=56s^c zzF}01P`~b%<6!Dh-)jttcX;zZ19-tW;w(b3c*%HsGSpx!s5a>LUy z(D}ZZM74S}_!ySnuo%Cl2~Hk^DM@4T@<>2nP(>%wMR~nDO4(;?rIRb(J=u;5u1&l_M@4)9fe+}{Fn@mC@ zxAUY1M4;N_eHA_#Y0n21B9Ipg1y8&up0&pf#-Q-h+Y7*#y^9Nw!EXnq+yC0kn35INL)VaIU28hOtdK)`V%( zjv1<4)>dJN=cVT>i=O&@pCfrr#LH82F@5_|zU+L>el5I`cnU!qviZSWrI9UHqlCDw zQ1!=6@+(_=5z$vEL7xb^_(m{5HVZP3PpeM(VhOdq0ooPF;R7%hmO#GY$jA{t9+E1I z9>Jn_1A5wMb?bC#_a=L7tllh>79O4_@r5eUu0zbYLc*zLN}B9z;CxW+H`m8Qsl?@> zO`Kpw08z3hV6u-fakxzx_GWN=`p37n3R)C)GKw)EWtsH&l9-IgD)s4Bvn6bd!U85^ z&EOv6Z49fOM7PJK!>vh$=tu^c zFYA|l?Jp%L;f(U8Ua{2!7stYS+`#~d4bts}S{K?%{E6pS_<6$NPN#|GX=2y)TyUVk z0|OTu-Fz<%7bo5D>ogW`5|y2Huu1esw_H0P*qMM9YSdB=us)_$R5h=sGsd%B;+7QM ziOb`8wzh>5yfMGSH!kksOe<&d*&8gXOea&1 zO*-ZBWBQugLg8F#v?z7F&J8;EK!*I<7Q5*Qai)ofmyQ1AP;zuWgB}ec1<)O-_a^SiWE?IvVS_snN1TD3K5v9LN6E91+Ov)-#Wcem6>Mkr&C3O!k z&*Ji=W~N{IcMQekp@B4fQUo5i-SBFOmCygO$jxW30a7K3Xd}AY=bUb2L&Gl)3jsq0 z9DT+K?^%L1SKkyp{Kx0N8xhJcW;X}VS2C^nqY*FoD4sNih?z<17x0!2H#OirHn z%(u(_OfZoUweStHTp6*N?L}_zrQoHrRtmx142Ne{9D5|zNK`sQ&q?9RM~EBQ#ycIm zr)?H~6e|=jy-$I;Lbq+gbG}ob=c7)-lC^JGDZgCVB||_o^rX2vK`neB8T$L8kUH7% zBhUqZ5OyhS@$rUt)+@!FYl5cUr@r{?MgNJTJjUqzIl zPGme}OFQc?>VF|jkCYm`!q43~=$;$|QXb{Or49uy^Q}OJumXs;Y+THD(sL#09!skA zvRna~q9&%#FTg8ccpkIxi@mbOulu)bMS3b-HAN@lJtzp64NS90 zc%HMx{Rc95v_|^Y_Js{TP}<3w-0~8$2`uhzIXAb1|9YnWujlChde;5tITp^uxpl*S zV7=DFKczx`kdSHiIMdB=$a3mOvAhInW#*DF6wC^rGSNTocMln0QdJ|Fm~m_GrIW7X zmn8b^Z2*yX3sf8$H~rpEsDCH3?nBF!Iu)56LAyS!i!<6zyP{&B%rCpB46}UmWk?D~a&KPnv zq-;o+b?u;#UI!oLBGR~>BhyjL4`;hQ!YqmXwm!i}MC)5$hXuNhxTQ%+ z#KlCD&NBZp0u|`p9)HDED&1RYneM5v!tXY=j-D*tCiwv}mOleNLcMu=@P@PRJ(wpeCdCF%noDJC)jP}k|0(}*#H67{io?+v0Y;El zyLGPn?z>AR$!0pAUO|d)$P0O|H1%v^J14q3TzLwuoCD6lbhVG`y&1A)f=Zi6)$R&Z z7pwC4d)dI-N`!ojx9$ z$=tPgO1el}n%5&yu_RxCuy{EH?|k{_+g$FZvMAf0|E&xB`RTJ??+hl@-7Ziej7XcYu-TA*`KWWJ$A#xx354eObE&_lPYq3A!hDgMCeqIX&Mxnn^bb%haeeS;HZ2d z_#mm6i>4`5?^(n{vRcwow6D%nGF3fE|JU2}f&4pyjn;zDq^cB@{B6&-f0WmA<#-}B z^z^XN3-S&*0%oO=54F5vMAU+Owf7?+S@o>xt)#SYU?FKCB zXW4_asHNQEcvH?-ZyJ1t37utBkoqi{Hzon6sb88|tPNgEyFRD=wJ`6t0F3z{O zv++8u(un<{856FaN_62jJwnX*IUi=h@Z@z}Mz{v%hsCJigEZ2cPu%K4vR2ny;_7V% z1{M(xqMxR7Qj@rtRbS}86xtxpSDk9PqPK!Fv|A(kr6T_Gpd;>EWzB%^>WHA`9XvJ> z2M@K|FgSnY!M6FqFHn%yPEc}L>2?eWj4?P5*USUs>&l0*%3j&r^rc79r_YQWB>qJ} zU%AcGKr|V-FW|5{wh`zb{RpoI*4Lq#eD3Q{DmkHLiq z_l`^45hh+eGDTl=J2K{rnGrTc`)HQ2@_&Du@DaI!etj|$Cp1p{lKJZGSLBXfL81W` z^^-SLx`VYmmhUA5`e8NqjCb~qb1v)%X}a)YOeaiB8mV~~s&QLTWJF#7srPEwpH`}E z>cT0fpG@QOnEo`!XP@AYJXG^3Ml|b&MG@=suY9i#iDHvAS>}AhOY0B`Or8%uSTv=b zB^TI8Q`sL}{yp%tIt8MHGXPBE`4IW{u2sZk3=mE>ppV{_- z;peY)pbZRgyqNzGJ2S&GgP*5|T*#GV4tNtbZhZsVsytZ!>6sOw!0mme=cH!S`2=_+ zZx$AK^64Sp)3I<*)|9A9`g+77zEx15a3D=PDL3<1)c1|-9K9w2{9>G7Pdatx4LeXb zL}aPc@0*j=|6uZYtFGsl##eMRILj|S`BP4D6i!LB3=9MvsjgyF^%pg?e-13NS0l-7 zs5yo`mcM0i!1vDeIakEN)UV{v0b4nqdqEj_sL#V1daA+LR=O0+lauW8HJob;>r zUV$D(EZ@pIVO7)SL&sgYYf^(yW1cX(x~H_uvF952O3zW|K)HP!sUBZAw{u!p2Nrdu z>aVUH!>Wk&+BU75$@K7oDIh0n_>{MHf#Z@1B#lCyoGe;Fli^w7F)zMz4L5r9xfOg1vj_sj{CgR)VKsH zP|SO{T0-$XLw>x^NAay*>ABKHL(baD)h+aKu zln^p;8GOZ*lfS{yw7ANo*bLtx$<9)EIu4VnMHErF)%w=!M2Xa|ne2S5yT60=$6I*~ z?OykUW6NRw(d?s$gkS5O`}&0Pf>+k}JuC5k#Hi^o0*#%c34%`x=adtQeD-`TJN8m) zfjrl;e%}y(BwWteEGjp6nf|LSEmVzgP}HIUFp5?K&$HD(5&n8L=^gu7wUPZ?nl;JU zD^(Kg4Cb0C2Vb}_`4l{uZ`D%i{vzfyp#jr)!)}73b=SHJFpteXXM+RwL>x|%1XLL0 zHOoaa&i#ad!^um_B5$s07pA{@hqy2r3u(q3ICeuSm&fW;KxoBff&v37jD9-@@v*YtU z={;QQ|LG0~%_2iJy-RpKr(ww`*9!mo`ywck^QeQktLQqyJ$W$0AW=(dA}46xIK&G4 zxJB3Q!bmeE?D9qM1Clcob2n52xlX7gdM_Rx{wJ0FQFh~Js@O*?Tn8Wx#NEXd2=E&h zFdY52FUG_a#!_>^kJ?e$=SkLAx$ho?o)67X5#Zt1dwA+9+1lS<9yxURRWM$i{qLC7 zo!&@S!Rr-T3%ay;Y%Q)-vFoo`wlP-tn6kKNsa3hyw<`3o1?#0JVjcVnBj9i+uO+VT zV7SwISxPtqEKu6j*#osdv&-%;)DHQ1CwjzV4>!w>BW%Fu=`pLN4PpkEw-M0J$;Rw83fuTk+`e^L~%ZYHD2!C6VEfgO#v(c~Pa;=XAM% z8{?xb5R=*N%(eZ98cUH8@N}50P;^iE9H-E)R_b;Y{i5(MvNik;YrGbv5!)4JOdcuL zjX|H_bNx8^{D5NSfM3}h1t?i=l&@aH4GpC0grkBSDu69+k+3NU#Agu4!vmgS{a;rP zJ)yk+vHtT>->Dx8f%^y&OI1ul5=r4aJ#^_kIKBr3X&ASZ!y?$Ar5mD*DV$4v0b5VD z3lsXhm?Bczm)a43ru3|yK7m zk0+#Q85KnjsC&Jr!&iL;9gbeFZ3?*D6o*O&FRHdIamZWHk6-sQ$}&EHb76ib7Q$DS z%uWPgz2}q1q=7Vr2R;uihBJJ5V7>8Sq&`0al@)njgd-z!hLi_#FEc>xyW1cy>g{~M zy(RhC?%gJZ$6$eYsS5>}Q1VNa{pYlbj>#zT-