From 5057f9fae66396ba730f29640c04dc49259f58e0 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Sun, 8 Mar 2020 20:48:42 -0400 Subject: [PATCH 01/25] Fixed issue 46; start new version --- DESCRIPTION | 4 ++-- NEWS.md | 6 ++++++ src/bayesreg.c | 13 +++++++++---- tests/testthat/test-special-functions.R | 1 + 4 files changed, 18 insertions(+), 6 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index f44032b4..f70b7bdf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BAS -Version: 1.5.5 -Date: 2020-1-24 +Version: 1.5.6 +Date: 2020-3-3 Title: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling Authors@R: c(person("Merlise", "Clyde", email="clyde@duke.edu", role=c("aut","cre", "cph"), diff --git a/NEWS.md b/NEWS.md index d3d2e9f5..ece336c5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# BAS 1.5.6 + +## Bug Fixes + +* Fixed NaN results for the Laplace approximation in the 2F1 function (issue #46) and added unit test. + # BAS 1.5.5 ## Changes diff --git a/src/bayesreg.c b/src/bayesreg.c index e1c73013..358873e6 100644 --- a/src/bayesreg.c +++ b/src/bayesreg.c @@ -390,7 +390,7 @@ double logBF_hyperGprior_laplace(double R2, int n, int p, double alpha) n = sample size p = number of rank of X (including intercept) alpha = prior hyperparameter - n and p are adjusted by subtrating off one + n and p are adjusted by substracting off one because of the flat prior on the intercept */ @@ -400,12 +400,17 @@ double logBF_hyperGprior_laplace(double R2, int n, int p, double alpha) dp = (double) p - 1.0; /* Laplace approximation in terms of exp(tau) = g */ /* Agrees with Mathematica but not paper */ - if (p == 1 || dn <= dp) logmarg = 0.0; + if (p == 1 || dn <= dp) { + logmarg = 0.0; + } else { ghat = (-4.+ alpha + dp + (2. - dn)*R2 - sqrt(-8.*(-2. + alpha + dp)*(-1.0 + R2) + (-4. + alpha + dp + (2.-dn)* R2)*(-4. + alpha + dp + (2.-dn)* R2)))/(2.*(-2. + alpha + dp)*(-1. + R2)); - if (ghat <= 0.0) { Rprintf("ERROR: In Laplace approximation to logmarg, ghat = %f R2 = %f p = %d n = %d\n", ghat, R2, p,n);} + if (ghat <= 0.0) { + Rprintf("ERROR: In Laplace approximation to logmarg, ghat = %f R2 = %f p = %d n = %d\n", + ghat, R2, p,n); + } /* Substitute ghat (unconstrained) into sigma, simplify, then plug in ghat @@ -511,7 +516,7 @@ double log_laplace_2F1(double a, double b, double c, double z) // root1 = (-B - sqrt(D))/(2.*A); // root2 = (-B + sqrt(D))/(2.*A); - // ghat1 = (2.*b)/(c + b*(z - 2.0) - a*z + sqrt(c*c + 4*a*b*z - 2.0*(a + b)*c *z + (a -b )*(a-b)*z*z)); + ghat = (2.*b)/(c + b*(z - 2.0) - a*z + sqrt(c*c + 4*a*b*z - 2.0*(a + b)*c *z + (a -b )*(a-b)*z*z)); // ghat2 = -(c +b*(-2. + z) - a*z + sqrt(c*c + 4*a*b*z - 2.0*(a + b)*c *z + (a -b )*(a-b)*z*z))/(2.*(b - c)*(z - 1.)); diff --git a/tests/testthat/test-special-functions.R b/tests/testthat/test-special-functions.R index 22dd7df5..14afb4bb 100644 --- a/tests/testthat/test-special-functions.R +++ b/tests/testthat/test-special-functions.R @@ -27,5 +27,6 @@ expect_equal(TRUE, hypergeometric2F1(1,1,5, 1.0, method="Laplace", log=FALSE)>0) expect_warning(hypergeometric2F1(3,1,1000, .999)) +expect_equal(TRUE, !is.na(hypergeometric2F1(12, 1, 2, .85, method = "Laplace"))) }) From 30e85d9193eb0d53f98289a7cb5226e0b3d18f59 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Thu, 6 Aug 2020 11:10:45 -0400 Subject: [PATCH 02/25] Update of required packages in vignette and conditional tests for suggests in testhat per CRAN requirements --- DESCRIPTION | 4 ++-- NEWS.md | 4 ++++ man/Hald.Rd | 6 ++++-- man/bodyfat.Rd | 6 ++++-- man/climate.Rd | 6 ++++-- man/protein.Rd | 6 ++++-- tests/testthat/test-bas-glm.R | 3 ++- vignettes/BAS-vignette.Rmd | 1 + 8 files changed, 25 insertions(+), 11 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index f70b7bdf..03b58166 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BAS Version: 1.5.6 -Date: 2020-3-3 +Date: 2020-8-5 Title: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling Authors@R: c(person("Merlise", "Clyde", email="clyde@duke.edu", role=c("aut","cre", "cph"), @@ -66,4 +66,4 @@ NeedsCompilation: yes ByteCompile: yes VignetteBuilder: knitr Encoding: UTF-8 -RoxygenNote: 7.0.2 +RoxygenNote: 7.1.0 diff --git a/NEWS.md b/NEWS.md index ece336c5..b8a67a43 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,8 @@ # BAS 1.5.6 +## Changes + +* Updated tests and examples to use Suggested packages conditionally per CRAN policy +"Error: poisson regression (@test-bas-glm.R#118)" ## Bug Fixes diff --git a/man/Hald.Rd b/man/Hald.Rd index 783f2cc3..d44f8277 100644 --- a/man/Hald.Rd +++ b/man/Hald.Rd @@ -5,12 +5,14 @@ \alias{Hald} \alias{hald} \title{Hald Data} -\format{\code{hald} is a dataframe with 13 observations and 5 variables +\format{ +\code{hald} is a dataframe with 13 observations and 5 variables (columns), Y: Heat evolved per gram of cement (in calories) X1: Amount of tricalcium aluminate X2: Amount of tricalcium silicate X3: Amount of tetracalcium -alumino ferrite X4: Amount of dicalcium silicate} +alumino ferrite X4: Amount of dicalcium silicate +} \source{ Wood, H., Steinour, H.H., and Starke, H.R. (1932). "Effect of Composition of Portland cement on Heat Evolved During Hardening", Industrial diff --git a/man/bodyfat.Rd b/man/bodyfat.Rd index 7ae8099b..0ffdc6f6 100644 --- a/man/bodyfat.Rd +++ b/man/bodyfat.Rd @@ -5,7 +5,8 @@ \alias{bodyfat} \alias{Bodyfat} \title{Bodyfat Data} -\format{A data frame with 252 observations on the following 15 variables. +\format{ +A data frame with 252 observations on the following 15 variables. \describe{ \item{Density}{a numeric vector for the density determined from underwater weighing} \item{Bodyfat}{percent body fat from Siri's (1956) equation} \item{Age}{age of individual in years} @@ -18,7 +19,8 @@ circumference (cm)} \item{"Thigh"}{thigh circumference (cm)} \item{"Knee"}{knee circumference (cm)} \item{Ankle}{ankle circumference (cm)} \item{Biceps}{bicep (extended) circumference (cm)} \item{Forearm}{forearm circumference (cm)} -\item{Wrist}{wrist circumference (cm)} }} +\item{Wrist}{wrist circumference (cm)} } +} \source{ These data are used to produce the predictive equations for lean body weight given in the abstract "Generalized body composition prediction diff --git a/man/climate.Rd b/man/climate.Rd index b6655033..b904178e 100644 --- a/man/climate.Rd +++ b/man/climate.Rd @@ -4,7 +4,8 @@ \name{climate} \alias{climate} \title{Climate Data} -\format{Scientists are interested in the Earth's temperature change since the last +\format{ +Scientists are interested in the Earth's temperature change since the last glacial maximum, about 20,000 years ago. The first study to estimate the temperature change was published in 1980, and estimated a change of -1.5 degrees C, +/- 1.2 degrees C in tropical sea surface temperatures. @@ -28,7 +29,8 @@ The proxies are coded as\cr } \item{\emph{T/M}}{, an indicator of whether it was a terrestrial or marine study (T/M), which is coded as 0 for Terrestrial, 1 for Marine;} -\item{ \emph{latitude}}{the latitude where the data were collected.}}} +\item{ \emph{latitude}}{the latitude where the data were collected.}} +} \source{ Data provided originally by Michael Lavine and available at \url{https://stat.duke.edu/sites/stat.duke.edu/files/climate.dat} } diff --git a/man/protein.Rd b/man/protein.Rd index fe8ab68d..8a30bdc0 100644 --- a/man/protein.Rd +++ b/man/protein.Rd @@ -4,7 +4,8 @@ \name{protein} \alias{protein} \title{Protein Activity Data} -\format{\code{protein} is a dataframe with 96 observations and 8 predictor +\format{ +\code{protein} is a dataframe with 96 observations and 8 predictor variables of protein activity: \tabular{llll}{ [,1] \tab buf \tab factor \tab Buffer \cr [,2] \tab pH \tab numeric \tab \cr [,3] \tab NaCl \tab numeric \tab \cr [,4] \tab con \tab numeric \tab protein concentration\cr [,5] \tab @@ -12,7 +13,8 @@ ra \tab factor \tab reducing agent\cr [,6] \tab det \tab factor \tab detergent\cr [,7] \tab MgCl2 \tab numeric\tab \cr [,8] \tab temp \tab numeric\tab (temperature)\cr [,9] \tab prot.act1 \tab numeric\tab \cr [,10] \tab prot.act2 \tab numeric\tab \cr [,11] \tab prot.act3 \tab numeric\tab -\cr [,12] \tab prot.act4 \tab numeric\tab protein activity }} +\cr [,12] \tab prot.act4 \tab numeric\tab protein activity } +} \source{ Clyde, M. A. and Parmigiani, G. (1998), Protein Construct Storage: Bayesian Variable Selection and Prediction with Mixtures, Journal of diff --git a/tests/testthat/test-bas-glm.R b/tests/testthat/test-bas-glm.R index d5d2d067..8d234463 100644 --- a/tests/testthat/test-bas-glm.R +++ b/tests/testthat/test-bas-glm.R @@ -115,6 +115,7 @@ test_that("missing data arg", { }) test_that("poisson regression", { + if (requireNamespace("glmbb", quietly=TRUE)) { data(crabs, package = "glmbb") crabs.bas <- bas.glm(satell ~ color * spine * width + weight, data = crabs, @@ -124,7 +125,7 @@ test_that("poisson regression", { prob.rw = .95 ) expect_null(plot(crabs.bas)) - expect_equal(0, sum(crabs.bas$shrinkage > 1)) + expect_equal(0, sum(crabs.bas$shrinkage > 1))} }) test_that("glm_fit", { diff --git a/vignettes/BAS-vignette.Rmd b/vignettes/BAS-vignette.Rmd index 1f122ae7..13d88048 100644 --- a/vignettes/BAS-vignette.Rmd +++ b/vignettes/BAS-vignette.Rmd @@ -16,6 +16,7 @@ vignette: > # require(knitr) require(MASS) require(dplyr) +require(GGally) ``` The `BAS` package provides easy to use functions to implement Bayesian Model Averaging in linear models and generalized linear models. From 83399a577540768b2683e991cd80156737321aa0 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Thu, 6 Aug 2020 14:48:35 -0400 Subject: [PATCH 03/25] added crab data as package glmbb has been archived to eliminate the dependency --- DESCRIPTION | 5 ++-- R/bas_glm.R | 6 ++-- R/data.R | 47 +++++++++++++++++++++++++++++ data/crabs.rda | Bin 0 -> 1320 bytes man/bas.glm.Rd | 6 ++-- man/crabs.Rd | 55 ++++++++++++++++++++++++++++++++++ tests/testthat/test-bas-glm.R | 5 ++-- 7 files changed, 112 insertions(+), 12 deletions(-) create mode 100644 data/crabs.rda create mode 100644 man/crabs.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 03b58166..a5888158 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BAS Version: 1.5.6 -Date: 2020-8-5 +Date: 2020-8-6 Title: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling Authors@R: c(person("Merlise", "Clyde", email="clyde@duke.edu", role=c("aut","cre", "cph"), @@ -26,7 +26,6 @@ Suggests: rmarkdown, roxygen2, dplyr, - glmbb, pkgdown, testthat, covr @@ -66,4 +65,4 @@ NeedsCompilation: yes ByteCompile: yes VignetteBuilder: knitr Encoding: UTF-8 -RoxygenNote: 7.1.0 +RoxygenNote: 7.1.1 diff --git a/R/bas_glm.R b/R/bas_glm.R index b04b88d4..b613ec23 100644 --- a/R/bas_glm.R +++ b/R/bas_glm.R @@ -215,15 +215,15 @@ normalize.initprobs.glm <- function(initprobs, glm.obj) { #' betaprior=bic.prior(), family=binomial(), #' modelprior=uniform()) #' # Poisson example -#' if(requireNamespace("glmbb", quietly=TRUE)) { -#' data(crabs, package='glmbb') +#' +#' data(crabs) #' #short run for illustration #' crabs.bas = bas.glm(satell ~ color*spine*width + weight, data=crabs, #' family=poisson(), #' betaprior=EB.local(), modelprior=uniform(), #' method='MCMC', n.models=2^10, MCMC.iterations=5000, #' prob.rw=.95) -#' } +#' #' @concept BMA #' @concept variable selection #' @family BMA functions diff --git a/R/data.R b/R/data.R index d4c16527..5baee1b7 100644 --- a/R/data.R +++ b/R/data.R @@ -193,4 +193,51 @@ NULL #' Bayesian Variable Selection and Prediction with Mixtures, Journal of #' Biopharmaceutical Statistics, 8, 431-443 #' @keywords datasets +#' +NULL + +#' @title Horseshoe Crab Data +#' +#' @name crabs +#' @docType data +#' @description Data on horseshoe crabs (\emph{Limulus polyphemus}). +#' Responseis number of males surrounding a breeding female, +#' color (factor), condition (factor), weight (quantitative), +#' and width (quantitative) of the female. +#' @usage crabs +#' @format A data frame with 173 observations on 6 variables. +#'Individuals (rows of the data frame) are female horseshoe crabs. +#'Variables other than \code{satell} refer to these females. The variables are +#' \describe{ +#' \item{color}{color. The colors given in +#' Agresti are \dQuote{light medium}, \dQuote{medium}, \dQuote{dark medium}, +#' and \dQuote{dark}. Here they are abbreviated to \code{light}, +#' \code{medium}, \code{dark}, and \code{darker}, respectively.} +#' \item{spine}{spine condition. The conditions given in Agresti are +#' \dQuote{both good}, \dQuote{one worn or broken}, and +#' \dQuote{both worn or broken}. +#' Here they are abbreviated to \code{good}, \code{middle}, \code{bad}, +#' respectively.} +#' \item{width}{carapace width in centimeters} +#' \item{satell}{number of satellites, which males clustering around the +#' female in addition to the male with which she is breeding.} +#' \item{weight}{weight in grams.} +#' \item{y}{shorthand for \code{as.numeric(satell > 0)}.} +#' } +#' @details +#' Quoting from the abstract of Brockmann (1996). \dQuote{Horseshoe crabs +#' arrive on the beach in pairs and spawn \ldots during \ldots high tides. +#' Unattached males also come to the beach, crowd around the nesting couples +#' and compete with attached males for fertilizations. Satellite males form +#' large groups around some couples while ignoring others, resulting in +#' a nonrandom distribution that cannot be explained by local environmental +#' conditions or habitat selection.} +#' @source Agresti, A. (2013) \emph{Categorical Data Analysis}, +#' Wiley, Hoboken, NJ., Section 4.3.2, +#' \url{http://www.stat.ufl.edu/~aa/cda/data.html} +#' +#' Brockmann, H. J. (1996) +#' Satellite Male Groups in Horseshoe Crabs, \emph{Limulus polyphemus}, +#' \emph{Ethology}, \bold{102}, 1--21. +#' @keywords datasets NULL diff --git a/data/crabs.rda b/data/crabs.rda new file mode 100644 index 0000000000000000000000000000000000000000..dcd3285bcbaa72375bd9aa3657f39320fa70004b GIT binary patch literal 1320 zcmV+@1=so?iwFP!000001EpBKk6T3$Uq9Q=UzZ~iX(CZk(QwVd`ST7VAu5Cr98yGu z6WhrNUyOW)g9=wv&`>25RFp(~HB<-*DhN{e8$gjrl`2)pG4p1A@AgghJ)e}7_V(?} z>~DVavFB&M`O2m3=B1`{uIwrkaOlxa#{8>KLXBs=6O+K8= zsI%SqL@UzDc(iyA7W_O};3CWpu3lS0bvPLxd~o>KU73xJfSCH3j;<`{3wnCw-5#*B zofUbE^CHR9<(-QH&qY|iZB$?v!VSe!v|lIpW3#e#K5UG69P{#AWOJ5jF>DNxCo75} zFW2kxnWwWn&JsgeT0J^((%oh4?(*f)7+2@&d_GR^$H90!o5-Ou4h64-+oXAGIC$+g zjX#GwE5~ElbFy}BiyZql_G;Jpahu1mt!wutk1vNE#p|m5U|aMrx-5SGqS$v;zw5HN z{b@3WZPj;GpYnuX@KU>h%6lsBgYR9z+co^Uz`;2am&Iu-ufF0z-qU#cdTvko4pcu- zd*q{Q^n1!Hjo;*9N6&*N>eg30#NS6;t8AV-*qb^GbOm430eGl4_Na3keW7btd3Mwe ze0xTp=G*c@e$f|*7x_2&LVRu17yM1#P*3y)@DV3)(!O40{h(f!C+dX$?Q1=Hx~BM$ zz?k@dv5O`;z!)3uISIc_8)koo;!v=>VtdIMqbt5uK9nK-Uq~iK0-d!zH>3_ zzi0Xh`9gh^4-LG6)OFv`wc?6LOV>?3i;oB%RhNDF$48sfwd&#ce(FXSdn?Yxx;FG| z?D>LC)syq{Q=Z-*Rqsqm{Wm>5>uO$HvC)CM&)JX3`U}}Q4?pi?H}oRwe^CK1e9ZIM z@chL(BhSvS2VRA#UyG%^+qy@T zs1NH8X}&1$lCagp8uiwgY1FBzBkY7^S+gi~ zfL6*$;|q=HzQn|Q8t)2MURfT~dsVG2=a>)Ua-6)cd?!Bnnf&CW;UimX8_Q_Z*YrgErl8QI`Y5yIm!MKen}k9 zD()SvhxQWt^Xa4M)C*wK)tDg(3c_Z(CMMxIXS)3x> zlYZr~9IU@Qn-KR9*AQ+MJ7;uym9M)p^Z>!ZU7e&vYM zpFZZZdHPey(fbG2C)gYxjF;E&&F^TmoXlp}9Z$U9R7xM`d|hmgN6XRHbU|(nmBwOz eyahUZzoHfWIXO9b%gJ}kNB;w{L!qMT6aWB{pTy1p literal 0 HcmV?d00001 diff --git a/man/bas.glm.Rd b/man/bas.glm.Rd index 6a2f0524..f15153f7 100644 --- a/man/bas.glm.Rd +++ b/man/bas.glm.Rd @@ -228,15 +228,15 @@ pima.BIC = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7, betaprior=bic.prior(), family=binomial(), modelprior=uniform()) # Poisson example -if(requireNamespace("glmbb", quietly=TRUE)) { - data(crabs, package='glmbb') + + data(crabs) #short run for illustration crabs.bas = bas.glm(satell ~ color*spine*width + weight, data=crabs, family=poisson(), betaprior=EB.local(), modelprior=uniform(), method='MCMC', n.models=2^10, MCMC.iterations=5000, prob.rw=.95) -} + } \references{ Li, Y. and Clyde, M. (2019) Mixtures of g-priors in Generalized diff --git a/man/crabs.Rd b/man/crabs.Rd new file mode 100644 index 00000000..5f74256b --- /dev/null +++ b/man/crabs.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{crabs} +\alias{crabs} +\title{Horseshoe Crab Data} +\format{ +A data frame with 173 observations on 6 variables. +Individuals (rows of the data frame) are female horseshoe crabs. +Variables other than \code{satell} refer to these females. The variables are + \describe{ + \item{color}{color. The colors given in + Agresti are \dQuote{light medium}, \dQuote{medium}, \dQuote{dark medium}, + and \dQuote{dark}. Here they are abbreviated to \code{light}, + \code{medium}, \code{dark}, and \code{darker}, respectively.} + \item{spine}{spine condition. The conditions given in Agresti are + \dQuote{both good}, \dQuote{one worn or broken}, and + \dQuote{both worn or broken}. + Here they are abbreviated to \code{good}, \code{middle}, \code{bad}, + respectively.} + \item{width}{carapace width in centimeters} + \item{satell}{number of satellites, which males clustering around the + female in addition to the male with which she is breeding.} + \item{weight}{weight in grams.} + \item{y}{shorthand for \code{as.numeric(satell > 0)}.} + } +} +\source{ +Agresti, A. (2013) \emph{Categorical Data Analysis}, +Wiley, Hoboken, NJ., Section 4.3.2, + \url{http://www.stat.ufl.edu/~aa/cda/data.html} + + Brockmann, H. J. (1996) + Satellite Male Groups in Horseshoe Crabs, \emph{Limulus polyphemus}, + \emph{Ethology}, \bold{102}, 1--21. +} +\usage{ +crabs +} +\description{ +Data on horseshoe crabs (\emph{Limulus polyphemus}). + Responseis number of males surrounding a breeding female, + color (factor), condition (factor), weight (quantitative), + and width (quantitative) of the female. +} +\details{ +Quoting from the abstract of Brockmann (1996). \dQuote{Horseshoe crabs + arrive on the beach in pairs and spawn \ldots during \ldots high tides. + Unattached males also come to the beach, crowd around the nesting couples + and compete with attached males for fertilizations. Satellite males form + large groups around some couples while ignoring others, resulting in + a nonrandom distribution that cannot be explained by local environmental + conditions or habitat selection.} +} +\keyword{datasets} diff --git a/tests/testthat/test-bas-glm.R b/tests/testthat/test-bas-glm.R index 8d234463..c7eab696 100644 --- a/tests/testthat/test-bas-glm.R +++ b/tests/testthat/test-bas-glm.R @@ -115,8 +115,7 @@ test_that("missing data arg", { }) test_that("poisson regression", { - if (requireNamespace("glmbb", quietly=TRUE)) { - data(crabs, package = "glmbb") + data(crabs, package="BAS") crabs.bas <- bas.glm(satell ~ color * spine * width + weight, data = crabs, family = poisson(), @@ -126,7 +125,7 @@ test_that("poisson regression", { ) expect_null(plot(crabs.bas)) expect_equal(0, sum(crabs.bas$shrinkage > 1))} -}) +) test_that("glm_fit", { data(Pima.tr, package = "MASS") From 765512eb87ba598b6fc3a5c489c9fe3c1c588fcf Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Thu, 6 Aug 2020 15:03:00 -0400 Subject: [PATCH 04/25] fixed crabs documentation and loading package in example/tests --- R/bas_glm.R | 8 ++++---- R/data.R | 11 ++++++----- man/bas.glm.Rd | 8 ++++---- man/crabs.Rd | 3 --- 4 files changed, 14 insertions(+), 16 deletions(-) diff --git a/R/bas_glm.R b/R/bas_glm.R index b613ec23..3458f53f 100644 --- a/R/bas_glm.R +++ b/R/bas_glm.R @@ -189,8 +189,7 @@ normalize.initprobs.glm <- function(initprobs, glm.obj) { #' @keywords GLM regression #' @examples #' -#' library(MASS) -#' data(Pima.tr) +#' data(Pima.tr, package="MASS") #' #' #' # enumeration with default method="BAS" @@ -216,9 +215,10 @@ normalize.initprobs.glm <- function(initprobs, glm.obj) { #' modelprior=uniform()) #' # Poisson example #' -#' data(crabs) +#' data(crabs, package="BAS") #' #short run for illustration -#' crabs.bas = bas.glm(satell ~ color*spine*width + weight, data=crabs, +#' crabs.bas = bas.glm(satell ~ color*spine*width + weight, +#' data=crabs, #' family=poisson(), #' betaprior=EB.local(), modelprior=uniform(), #' method='MCMC', n.models=2^10, MCMC.iterations=5000, diff --git a/R/data.R b/R/data.R index 5baee1b7..7f22e447 100644 --- a/R/data.R +++ b/R/data.R @@ -196,15 +196,16 @@ NULL #' NULL -#' @title Horseshoe Crab Data +#' Horseshoe Crab Data #' -#' @name crabs -#' @docType data -#' @description Data on horseshoe crabs (\emph{Limulus polyphemus}). +#' Data on horseshoe crabs (\emph{Limulus polyphemus}). #' Responseis number of males surrounding a breeding female, #' color (factor), condition (factor), weight (quantitative), #' and width (quantitative) of the female. -#' @usage crabs +#' +#' @name crabs +#' @docType data +#' #' @format A data frame with 173 observations on 6 variables. #'Individuals (rows of the data frame) are female horseshoe crabs. #'Variables other than \code{satell} refer to these females. The variables are diff --git a/man/bas.glm.Rd b/man/bas.glm.Rd index f15153f7..a35e08da 100644 --- a/man/bas.glm.Rd +++ b/man/bas.glm.Rd @@ -202,8 +202,7 @@ are mixtures of g-priors that provide approximations to the power prior. } \examples{ -library(MASS) -data(Pima.tr) +data(Pima.tr, package="MASS") # enumeration with default method="BAS" @@ -229,9 +228,10 @@ pima.BIC = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7, modelprior=uniform()) # Poisson example - data(crabs) + data(crabs, package="BAS") #short run for illustration - crabs.bas = bas.glm(satell ~ color*spine*width + weight, data=crabs, + crabs.bas = bas.glm(satell ~ color*spine*width + weight, + data=crabs, family=poisson(), betaprior=EB.local(), modelprior=uniform(), method='MCMC', n.models=2^10, MCMC.iterations=5000, diff --git a/man/crabs.Rd b/man/crabs.Rd index 5f74256b..9af67d1f 100644 --- a/man/crabs.Rd +++ b/man/crabs.Rd @@ -34,9 +34,6 @@ Wiley, Hoboken, NJ., Section 4.3.2, Satellite Male Groups in Horseshoe Crabs, \emph{Limulus polyphemus}, \emph{Ethology}, \bold{102}, 1--21. } -\usage{ -crabs -} \description{ Data on horseshoe crabs (\emph{Limulus polyphemus}). Responseis number of males surrounding a breeding female, From f93f3f6635eeeff2210c2917f2906dc9939b9d0b Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Mon, 24 Aug 2020 14:03:06 -0400 Subject: [PATCH 05/25] update to submit to CRAN --- DESCRIPTION | 2 +- cran-comments.md | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index a5888158..dc5bb39c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BAS Version: 1.5.6 -Date: 2020-8-6 +Date: 2020-8-24 Title: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling Authors@R: c(person("Merlise", "Clyde", email="clyde@duke.edu", role=c("aut","cre", "cph"), diff --git a/cran-comments.md b/cran-comments.md index c750a5a2..a90bb714 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,6 +1,6 @@ -# BAS 1.5.5 Comments to CRAN +# BAS 1.5.6 Comments to CRAN -resubmission to fix errors found in package version 1.5.4: +resubmission to fix errors found in package version 1.5.5: CRAN repository db overrides: X-CRAN-Comment: Archived on 2020-01-20 as check issues were still not From a33732dbe4963233db59e68b198332417a836092 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Tue, 8 Sep 2020 21:37:53 -0400 Subject: [PATCH 06/25] add mean.x to glm objects before passing to predict.bas --- DESCRIPTION | 2 +- R/predict.R | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index dc5bb39c..9d0e0c00 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -12,7 +12,7 @@ Authors@R: c(person("Merlise", "Clyde", email="clyde@duke.edu", person("Yingbo", "Li", role="ctb"), person("Don", "van de Bergh", role="ctb")) Depends: - R (>= 3.0) + R (>= 3.5) Imports: stats, graphics, diff --git a/R/predict.R b/R/predict.R index dc1c3a10..190f4fe6 100644 --- a/R/predict.R +++ b/R/predict.R @@ -75,6 +75,10 @@ predict.basglm <- function(object, top <- 1 } + # predict.bas code assumes intercept is based on centered X's + # add mean.x to object that is all zeros. + + object$mean.x = rep(0, object$n.vars - 1) # get predictions on linear predictor scale pred <- predict.bas( object, From c4eb072b36ff467ac265b166f40566ba30a789d9 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Wed, 9 Sep 2020 16:38:21 -0400 Subject: [PATCH 07/25] Addresses initial issue of run time error with se.fit = TRUE for issue #50, however se.fit values returned are not correct in glms --- R/predict.R | 21 ++++++++++++++++----- tests/testthat/test-predict.R | 27 +++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 5 deletions(-) diff --git a/R/predict.R b/R/predict.R index 190f4fe6..7d8098c4 100644 --- a/R/predict.R +++ b/R/predict.R @@ -75,10 +75,7 @@ predict.basglm <- function(object, top <- 1 } - # predict.bas code assumes intercept is based on centered X's - # add mean.x to object that is all zeros. - object$mean.x = rep(0, object$n.vars - 1) # get predictions on linear predictor scale pred <- predict.bas( object, @@ -571,12 +568,21 @@ fitted.bas <- function(object, df <- object$df[best] + mean.x = object$mean.x # glms don't have centered X for intercept so need t + # to center X and newX to get the right hat values with orthogonal X + if (is.null(mean.x)) { + mean.x =colMeans(object$X[,-1]) + X = sweep(X, 2, mean.x) + } + + + shrinkage <- object$shrinkage[best] if (insample) { xiXTXxiT <- hat(object$X[, model + 1]) - 1 / n } else { X <- cbind(1, X[, model[-1], drop = FALSE]) - oldX <- (sweep(object$X[, -1], 2, object$mean.x))[, model[-1]] + oldX <- (sweep(object$X[, -1], 2, mean.x))[, model[-1]] #center # browser() XRinv <- X %*% solve(qr.R(qr(cbind(1, oldX)))) xiXTXxiT <- apply(XRinv^2, 1, sum) - 1 / n @@ -606,6 +612,11 @@ fitted.bas <- function(object, df <- object$df[best] + mean.x = object$mean.x + if (is.null(mean.x)) { + mean.x =colMeans(object$X[,-1]) + Xnew = sweep(Xnew, 2, mean.x) + } shrinkage <- object$shrinkage[best] if (insample) { @@ -620,7 +631,7 @@ fitted.bas <- function(object, } else { Xnew <- cbind(1, Xnew) - Xold <- cbind(1, sweep(object$X[, -1], 2, object$mean.x)) + Xold <- cbind(1, sweep(object$X[, -1], 2, mean.x)) xiXTXxiT <- sapply( object$which[best], FUN = function(model, Xnew, Xold) { diff --git a/tests/testthat/test-predict.R b/tests/testthat/test-predict.R index c1eb2c54..b037fad7 100644 --- a/tests/testthat/test-predict.R +++ b/tests/testthat/test-predict.R @@ -47,3 +47,30 @@ test_that("predict.bas.glm", { se.fit = TRUE)) #expect_null(plot(confint(pima_pred))) }) + +test_that("se.fit.glm", { + data("Pima.tr", package="MASS") + data("Pima.te", package="MASS") + +pima.bic = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7, + method="BAS", + betaprior=bic.prior(n=200), family=binomial(), + modelprior=beta.binomial(1,1)) + +pred.bic = predict(pima.bic, newdata=Pima.te, se.fit = TRUE, top=1, type="link") + +form = paste("type ~ ", + paste0((pred.bic$best.vars[pred.bic$bestmodel[[1]] + 1])[- 1], + collapse = "+")) + +pima.glm = glm(form, data=Pima.tr, family=binomial()) +pred.glm = predict(pima.glm, newdata=Pima.te, se.fit=TRUE, type='link') + +expect_true(all.equal(pred.glm$fit, pred.bic$fit, check.attributes = FALSE)) + +# issue #50 in github regarding se.fit failing; debugging indicates se.fit is +# incorrect +# Should be expect_true +expect_false(all.equal(pred.glm$se.fit, pred.bic$se.fit, check.attributes = FALSE)) + +}) From c4b3e3477a224d5a5170f77fc22d4af995e15d3b Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Mon, 9 Nov 2020 22:46:48 -0500 Subject: [PATCH 08/25] added check if p-value is 0.0 or less than machine precision in eplogprob and created unit test close #54 --- R/eplogprob.R | 1 + inst/testdata/eplogp-testdata.Rdata | Bin 0 -> 81130 bytes tests/testthat/test-eplogprob.R | 22 ++++++++++++++++++++++ 3 files changed, 23 insertions(+) create mode 100644 inst/testdata/eplogp-testdata.Rdata create mode 100644 tests/testthat/test-eplogprob.R diff --git a/R/eplogprob.R b/R/eplogprob.R index d10fb75a..cf8fc97a 100644 --- a/R/eplogprob.R +++ b/R/eplogprob.R @@ -54,6 +54,7 @@ eplogprob = function(lm.obj, thresh=.5, max = 0.99, int=TRUE) { stop("Full model is not full rank, use `initprobs='marg-eplogp'` instead\n") } else { + pval[pval < .Machine$double.eps] = .Machine$double.eps prob = 1/(1 - exp(1)*pval*log(pval)) prob[pval > 1/exp(1)] = thresh prob[prob > max] = max diff --git a/inst/testdata/eplogp-testdata.Rdata b/inst/testdata/eplogp-testdata.Rdata new file mode 100644 index 0000000000000000000000000000000000000000..79cbb90af657fc1d876d1d8fbc38a15ae5ccdf08 GIT binary patch literal 81130 zcmV(tKk`FP&%x7#^}D&~?*=1f#n zRJ2rEsA#BY={8?d(U};liONz@ZK2sbpx!*excN+Db7k{gY@4t4Y(AJbPt20iGl%U5 zyGICB_psE$d3wUa`<)vjrIg^Wk)(gQe~{3=$M7eVpN6pJ-lm{jD@~|B)RN8Kc8Jv6 z8@!a2*g>jnm&tMUU?F78=)9Cd{t`wb7FnIN{|KSRk=+qjqoH=)uHd~`38{0+P)b-i zk~sd;cp&T&9U&a^-j*&zj#PT!r_GwzPo944-|1PIO`dRN>KurkA(XlAnrZNjlNNlJ zT)2+v62u{H|0u^SQqSfK&${?O;?(o!F;D*UCKNbcUmM7MNf^@5_`iKEMI02_&FRb@ zN}k@l(pkqVu{*nV;F$~tHfNB&suYpBm8G*UhqHZygR`nu6PJhp*^RqULr#Yh-& z;rIIMr(?ef<*kYuzBA`Yl^@r3vGKek^~A$<{+rezWNhUGuetaUdX5H-biq`FpoqmY z`saP5Va#WzZ2CFk$cmCgxVJX3&ywqx(ZlbA=JpvXzL%E>(FT87D zzbHk%1sX!%O8+4GCljX=|CTP(rV;w=mrg89Kf?v*-)}2&=m}$DSSZ|hnpBOrkwG1F zow&Sz&NyvCgD?v>XiX5$B{c&Xef#7Mi4$rr9kDuGK*%(C zD03!4KVcgwklOU>WiKtk|FV0#c7PaRarf5_x>N^JF}?iwu^>glJc>4l?lV1k=xXDv z@9jWBukiCThRXj)deZfInP@bjb1tNHSw@^ZWAcXf=E$M=mt+WaEBwKkEJ^5&@Y$%1nUN<{ww3z-iz00}+KlqAyd~86hAXc>N+M{7Nxsy@ zR|H$fZDynOI8tup^SSH)iIWC#0ya&&Rw|oum=-yRNKdT~gR>FnNZ9htLo_#T<0*8L58u$>XsK4VZgY`G`So z3!(4kI*ueWAqCM7+m{mvRihW}v&pT5sV!Bi-~DYQ?>}+YHx|sq*?au&B_F;gjV3<0 zo^!7zHRSazJY`b|RxZ^i0e?eC-OEm|5)!OPgS^H0UREJe@1^gDT_{2*+z@xZyir^7kaa9VL&in8iAOGbWWc zlH<6){3LYRcZ_t5dXvZeXtZ{R-Xb}hbTmhQE)a&z1)ZjwQh>42uT~a4;Ri)(i5+Sd-ziGx&lDH)JCL~>ConSEgRr?}$A9*@! zQJ?Yr4#J{D@~RTYc|wf))=@^!W`aNWIQv$u--N!hbVT1TFC6#v`P;I#lRS|qIsPAm zE+I`Y>SpFEk$aV+wX{>s35~&=rGpi$q*&>}Uw2yfkVgM$8Mzp85j^L$pY9l-B`zrs z9A#jsB+t9pmD^8Qk_sg=G`RwUgzkmWySF|Tkn*i_Y%({m5qzdgow?qj1TT+~xRaq9 zc_Mi%(AEAcp<;3WZP}Y3Lf4P0T<*^>sd<^ldyrTm#9IY!e|e=#>RkVPCv0Q0e^o=@ zX?#^BmHX2TS?bIQ&Ga(sy<0^I{&Wj#TYeMLe62ce&Bv8c=A#7e|9*_n7Jg@O?gkY> zeEg7UFK$opI@q%=2MCa-xr+||crZ>H$6u2^_a84Q{$F_rTU!SyyT!M<)@z)UT^Cn& z`hJ!)rB{|_GdoJ0Bn>hc=_E2~S|Fvkz^QDW#vA~^Mi%TvLr{6yP zb+w(H(0#^hHP&cL>R+2{8JFWF<;z96`;M`Y z*8YjDCu|=Q;?%K`mh8_+x%GpVvNNiL_W5Dyyhv@57hlK&U<0W;vYm zfRvCuYNWXOk5v7t-Yd2>ia0a+rG@wJ8lii!S>`y^Awue3DEE@i1YyDR`(jYJF=_T% z^_-@KKWVmV=Vbi8nB);1<+6!NAkEdAb^^^HalGKA!BJ~N;+*65J;J)#gwpp6*@@Ct zQl0g_RLVjdc_vPz_doF$1hGzMec-VGVZu@=p%$P3V`rQb?z?i3m%hJy?%qcstv)?6 z;R{|T&!~5YJz6dy%rdVQ?{wfI_I$Yh^e$>sdhg%BaNnn&=fAvpxyqBuqoJLp_LR(^0ir= zaopNT95oEN$Z+-?siD!rA-!>dG@v$4*eA|Hn#Ud(yA3gf5$n5?^RF05vs**O#l$5- zC!y@6nW-A7({Cq~Fm#;M65n0JQqVxkZJ}B1rez_O+kPj>Hr5kH8sRS)0%`~gi?iEi zgN~C56AYiydh7_*|K^L8>ZN#1e3-RN>ZV(dB z=Tn!Ja>+eXt&9ga9}tIj`q|wI=qFB}&KGWa?S`$NFT830P)1sthd({M`bh}daT%LroFydezU4Liz3gePni{6m?dNu6Woxei!=^aT3jo9Ln`zMMJa93 zB2MYb$d`L35E9hy{MV`&3E73eU+5n9kQed!w#Mx5gmTFBgbg7l;yBC9SBV2+B+nH+ zN{vwfp`2OQ8^gs%NIs7Y4RS6ex69g^v4Wok5tu= zX62^!Gh3QSwZncp&->3n)6}?JwV(r`Hl5w~x@~j*IGFKN45kuNqH0g4XI>DJWqzu6 zW~m4vS4OMEb{o=2+Aym@<|v^Q^J@P4(-WjA*K0w|OW#T1*@NEP8VAWER6pbd>9!NH zmdlm;+KPnqf$bX;hU%o$9%-YU67HnzzPc?6CS|1h$$)p~SZGLL>Tef1BbW(eg}+3> z2?J7PLLq|S*+FPMeSM5+wU01(DG~Ffp^_Bob78B`93n2UbO#LXPa_O2Y>oQin@O6~ z<~*2uUP>zWezi3zog~yk>8kbrx|2sLEKBMgdkM|C^M5a$$RU(Vmv)Fda*(FwzjMqT zW64v$o*!n4`c9l#4mN&id6iUkU26yvO(S$#d8nKopCL{^?z#VrZJAJ-)Hro6%a9bx zUd&uC5XC>zN7ZjuYm;ZvG6RRU{~%6zc05uQnjj?C8v1S%t;87{A&U5U3F73w77cTaGlYK64gZQB8B)IQrm10Y7ID1)a>gO;O+NDUv2IO0L0ppP8IpZZ5=zun zyLHqSVf7iWo0nxZsjDM`acgF{Hf_OBdef9Re0K4DTAC55!e6^gH7rjEUQb=U8+DXa zmwSJEubUkqKWl$7F4clK^^dh6!tVxYW%d5jH=SwX>}poEWq=B)yZdtEhIBf~NB((n zs9J|G_|LSKSRW(|bZVmqn1xA&eOFqPzvU7t7rgxs%uEpSk>d#&kE%$G`vLy%3TDY; z^Lq@9U&N4dxBk|^q3fiK-G7~8Y+a;mO^yD*lQQBm-%sP8&uq!#E?)H=oI#}FsYN5^ z>F=bbLP7rBg&LBh)_=~k#+f`lo1Edu(L+e_UzdInZcM6I@i7q~8>-cG2e`+nb1Axx@-9}Nw&{6w%i<&RqW#So&N)ip;L z1W1eDTfMGWyv6~QewShO|H#u@>$wv=-;>7=)AD@vP$hu!EcU{^Fj8z$^Tj&lB*Dn@ zuk?ko{{Q!MT%#3VNs+9@-Zi8Z%Blbgb6yhlFgMN=pW8@@i^a;A^fcRWQ7ABVu}Qb_ z!0E8T_FIA3I5o~pxqRv&4%&LoFlUlDtI1N^W>Nv&0#DPW1tM_mNYCefmI>@x3Ey0Q zDX7*S;{7o84l8tDQ(qCVf?iIpK#JD@Oq$+0F%w3Ijhcr&`~#h#;p*gGUEUh#Vzo=iDJ)(xa8m!u5HYQ<5c|}cD9u;^C`8D>dqyc+S`#=R?mxb z0TXqHKj~nfC|y*Dmp0Bkk~U5^Ho)G!yG1#6Z9wb4;Yc4@0cdh*y44u{8YjZd%g)XZ z6ErH3`ZF4rpz5r}p&KGKxT15dv4f=!Cfyeos}^T)+*;(n`&K15{ye_hDYvX-Q_zRG_3!ZU0x@u*x&NP#KmXRr3tyWsfFk97AYC^#44 zXhH4B2^E?7y%APikatAy42@eE3_Q_h>}vUq!)`>`4Vk9|)xAZDSJ6W_HF7wZURoVT zi71mpjW1wMW?)8n>^Zc@jBzaWUBJQb{_OXJd|}=1cg<@lHW*{r|HPPT0;;cgem9+9 zhECyKPF!d9!klBmhj{KH94PMMK25a)m$GE*$L1E`XKvckfrwFTc_0~4liv-$GM@hE z*!mmhIRz+xt0`Ef8^M-KT*by?)j?Oc*^LPWV62_>sl zFH*MA<1FpPJDR$5m_xQdv8(9Cr3sQ&`?o%<-F~o7ZP^-@dpbKEI$y%1UnWiHzgeJ= zQr9GHJz@6K@sR-6YM2(6XtVtA3kEhmc${u%f?GuPM%Vx%Lz4s zI?fPZN3mU45E1b*!=H(GKIRqP_Vnmj(_hGq*l-mVg$qlON|aTyZ*n?V8i` zv(PYp-#|P2I?i63552R`3)jB${9ZeA31`L+2$;#;z$&JR{vbgbc&#h@aK$Pas?v05 z4?TCmfxkn0YMakuO{0mhK2u07UQSZeMR2Dc z>4vJ86xxFS{yw*FlGIEK~ed@%r*8DW?(!2Pk^uUawV>ryzy)+5xD#ez6nb%J`dP8R1 z@BJ=&#c{cF`*rT6C>%`B*qhhdj+O5|(&o$BV(X>h^Wj`Vo9p0dKb3fkpf*J;# z#(sTpj5n`<#7Kjt*g02B+PVM4Qg94bR#j@Ow?Bj4cDXmkk@qlX$~Dga%XgePu&}US z&H~=rS-sldejNHs%Ad#Ie~u+D?W3BueSnz+?-TruiqQ8|S(zpz6owVw+X>TW;uM4a zT^1EqSPE6%tD3}x1NTq5j*{ecP1n$IVp|M|kzJ>Co0xv@@dO)~=X#x33b2bG|? z&EqW@=Z4MED{tTK(#M7^I_(@D(l}D0yK6Ew45u9r@Tz;|lGN*I`3IBk;YNek*?%2I zFtI)0(Rqd!urxQJB<8ad`u&SsU*unekE|Rt)zOQX@!ibz#FH!>6#aau%X0{u{N2Y4 zetgC`qkWz%!%{FNu~p@FnkNoMtQ?ljQ-;9{N`Dkv_^{lUmUq;d2dWPqZ@6l504Z-- zMXyAKVDsH`A4Es~Ko7fd)Ri+!kn#5MlDcp0DtFRN=nxYZX&F_6?j+-#H)5l( zkCV24qs<75-~3-NBH3@{eM`hh!McTxefj6v1KM+UP#L zN*JX7#^;+FfMZ721_rELb_Vs<;>ez) z)psI_1oh*>ml@2`;4f0RT>Xg=yQu=V1wA|jGqIiayBh?c+QoyD&S4P-Voule9J~j8 z{I^H$zT1s89Fm6{OV+Ud)%y9a*9|cE;F{_$8Yw7F<5;ZLd4oL`61~-af;eq4ytt$N zIj(+<9cC~2fkVBzBPr`zSae#LC24;XQkWEb|61&Y^l(Fa8(#%+waD zjvL%2Y2-1#T|Z{g#SXeW(YWQh3H0)%vN{U-VgKMk+78+ZoPV|@??06TI94^4Y-{`x zvoB^8O!7Q|iM*&}s|FhAZ}hCqiQ|WU)g}$r6GAXPT6S|Z&k{$vZ+=YC-iCjQQiA%p zVqoskPL>T`3z!Z$uGY}Q1PRIAaW;Z)aVqfUTa{HIoHBkH%zAhO`n|@@b*?{#ZlS0p ziT^I(^3_|lcm29HbHDKIk&Br)Z^H1N|4kNB)+T5M`QPK5?m*4>)@f`m4D_{mvYAWo zXmbq`)UohASA0{V3odb{&wU?@CK>3?8Q&j%0d@H^4|XoK!knJouXM@Fo9l9bZ@O?V zbly5Sclhc*C}%2C`6M2Qc^M)ft8C-3CGtVkF6t8KqNNRCRb7hi{|Z*bc}ZQ@p!*rlUbWq*+MTs#eoZR6>+^jE-fh5J1VsZGR| zN7A3>87{-~Cwp5uOuxX;iMN(E9eZ&u;roA-s4JND=T^fj<_R2n8X0jpAQo3O8XjeN zWnz9J+w~C|0~|XcBouJw70hjQl2LF8!l6Af^K~+XuzbP!08`@ToZFMIH+%g#wDFCe zJrO8`fppUr<)%TVq=ollQ1h=V7U> zp51ZaCz94b-y1%-2}!*&7yBP7VfT)*_nuv$PDngcD?H1R&`+Pi@h*R|j6eY6F8e?Lzw%TR<7w=2VH{`9yaDRZv8 zMgXdOxg-XP4nkE`!Vq_T4c5o1I)$*=t*P8|`$(L2QT^t7Ehj?8zhU#%*u>24W;RI3y)6%I-IBD}~+b9&Sqe zK0(v>=6E4_`SGiPkSDEA5lzV-X&nK@!n)c>0+#s^}ejsQ%z9q;}kQwoQkWF zAyrq-O5!p-oBgzHI;?$R6=HfBfz9O=+m5RLg04vsjhN6etPoSS-WI-sEe+pUS9;{& zTgM*PF}r`5H!1P!tlc!O>iX&K4L*!LU)o!|r_`}|@m&99%NO`hq>@IBD-kD(_N;py zih+rSU~0Zo>M$dJWK2@14*tFWL+_wy3+*E;=yJsk`#X~J6H^kfNyw$`bs-&es>Gfz zQF;p(*yK*=#pS?Q*N@Zp-Gy*K7yH!oMX@_`Wh;-fHZ>xv!VtEVp(rvO$AN2mbm~H zKPuF{djW8YirznMz70FBe7n%YSVmC$FD~5>`2&CFT|y0-7-25woK~Ww1a=3QNU{dYr^*FT!0yf2~4K7(QXvK%ZueqGY6MFq393AOFa9>A6*azT~2vv?9azvl9PES7*JF3s0R*OQ^Xq7VXQhv2vId(rIbUpVW2sqFq9 zVeC4z-IPJH5^9V@z1L}~pkc=M+9;DPw#B=%&Mx=jV073i*5B8mfs_nus$YO0_9>lz zb8T4gB=_w}PfLt+sZ>lxR)YHS8fG+1;h@UN3qE;Pu<}?kZC+NCpt?`A*Gc^#x5)dl zTxs#f^jCQ8$$55YY_045Tm1z~wLVDNIEmqaBmc5v?pzoB+JDO85rim^fUWT z|9?-1y5~m@iD|lFN!U>yzWa(WpO}}N+fRklUB~*aWV2yf#bd6SNjDgK|DR{x(+^nD zbA3ZwVg?6R&PHSmD&PR?fm}ky5EoqrZdtoDLIc0i@Umwz&ewe>Lb?l}lTP1U_r*>a zj&nFUSlo`eTs7(?%sZe_XI$F5#uOHo6y}3_1h7)*cc`n(Bj|47y~BLwA$+IEq=>bh zg4P%OwOoU3FzjwtNcd2i@M(5cp3H)VPafubS*#Wfs+)CBz_9yLAqq~=?5 zh+`*?v>$wUw*DnH94+HpuzHFM(d{PO3!>PtV@sXqnGu+bRIui6;=`#8Hiq=aGPr87 zKPKrH8o6P!gWn5QF1(?}Dtxm%Iz?>q z(w050jn}X;FH}v4#S^C^CqBzDeufF~-@+sqh`mn_*{XZ*$7ZShO#`!qNU2S}$k{Ck zJ+b08{A$!Vp-NRmeSySj6P4!zbrU$#uGl?u`XkJ^&YgQ}wv27;zN`r&ia2xVq7ZrX z4HkZE;Zd+!fOe+)-&X`qVArdfx6au`FhBNS|D;O|te@UGdhl5X_V^woj?1gyw3>S8 zp(=5p>~g)$8)gH22YhJ_MmGE5BR}#A+Tu{-=L9pZSs3A27x??(9Sq34iMT{RkF{Bn z(p)uJu<=N-m^<G165e9RADY@N@03$2&gQeRjX;)p(P!j5kbq3f(ledU}p z4#o1Fc0IonhsK6^IGwLR{(-dJ+Yac!OyZPR<;Z0iv3CejY;c362*;a&_vo?4;_=QM zHSe(djsIY9Q!kG4xCXErUm|IT>^?B>Qzod=RpgHS4TXi&u_8+JA()oW=1wYj39Eiw zadqQE&{Id;GSEE*1Hs`Io&pRoR_jPHewYc%kGJqlR2{-aiq`cQ$6oCB&i6mIp$K*J zD+-N?RnRT2$XMJo2kSdHUl+`-Vuw`zt!_CpoIi$#SWA|nZ8E?&Y@8N~FNv9_&8Lx=_>K&sN;vRMP`8I_Rn1;mC7YDi2VTe87g0p zM?qM$?|Wp9>X4>o-kGwi9>#?Phr!_uCcdNJ`ip8OmbA}%zKjSYnF?*E&iQa6S#gF6o}9W8Ceftna6?N`4cZ|${t1s=xtwRmjuZVR+}PW_2! z8-v=DB99Ztu3#P07?a(;Z0Hd=*SjRW8^$wFO()$GCv zg(D@9d5`0+{KXpRi0eEWo(0iIXlk5`5ayq>)i8rEPRzbj)PiPSb z+l?=*tgpd>w79mzfv3=`-L+*yMgr26;!0xoSz@pCpZS9q&coQjWX;LLAEEcAtdspu zcU=6LMNav>!{%38Z>TnA!pxHp<`uJZ_{HR2)xA&KVVH07QQvQG?AUuZ=FUz7g6eW| z5nDnWTq1o3HnAvOY7oLS_yBwcFJ=4`xHV>Z6MrypvdV zJ!#dF-3)tTwCmdU?Zw&$e{z#WgR%AOeigxY=~(|{Y|rJjJJ?M9MNq?KoTM`2)E9Kr zg`x!dQ;HM$Sl=WPC8ulw&7V1Pmn{zxRAuRZ@3IcT$BJ)7&T5q;P0AMj@U2z2D1EE2 zcJE!7vgs>r8X3XaZEQceAN51!a$$b=U_VeMU)`c-nT4v0yM*r5_7Pil%U==Nh=tXg zif(lFVX%7rYGp686n4Du96R`<8!8lM*IjRJ`t;(R;cIut2`Z+yygY;tOo(m-nK&y# z`X6UE<*Rbo`H|nJA$1STi9G1I+h_?(#|9tHpWnbf>tHPEHT&X> zH;$Kf3q&nCLDy8TnPKW1=r%Ss?BzMLnG5t41=81XDu++Zy3+~!CrqO6I6Gr!$Nh8P zkBz|8-K>kJoaAAtI9MX$Y9us%icZv2UV^b*`AQu=gP1xyD8YX91Z376J$7r&6K1X* z?2GM=#io%S&Inm*=oNX+)3~1+Yh5I-X20RVS%Y8Sj)?umLbnG6VJE_&{@BBtCUr&R zR&m9ERzV7}smlNCIVHhN1uK5cY&mkj3WyzK{z7aFs}XWEzD?qci4J?1G}$VI}QOK)am5Oef7m57DZ@K8ScHi^EkFXNYk(F_uaft$qADVIcy~g zYCb)2#MUb-f@vdVST*a!GxXO9raJ8`R*W|1+r@w19?<;2RxM?ZLRNi(MmMysSb++c zpZmsGeP)2F;iSKZ?@Qs#E9R_0+jlt4VKgcza0&;Fzv@g&WZ{fhaLK%R9S(CZY8@7h z!|KFn&obUBg8Cq5&gnlj*!wmmXrEahOz(A7Ww>GpJ<6lGr@gG8-&fqO+A9?vcP41r zZc)X&36GljBR^oE{MlPmu6ywB%-QKSIxQG2o%?l7@;HnbPF_CtC>2&mE-d$4AH^>3 zos7>FJaOvd)~{Ae-=WdD%INz|F8q6lBk`KXHkekwYcc08OVE@xvGJxwK*Jo*`LG8e zFtDe%KGuK%7vIf2>i2fUxsRWIJ$}psHTroDYNKV)Gu3)@SLI2PCVSI4g`3uDw_wbrA;M%G|v0 z>KrWFM_N@{9l^oi#crVxfX2T!PVO&!fI~8jF)@-C2Sb)-jp58<58DVle0ev$O4 z+8$ZkC1Kj8*iV>A9g~ZkGAzP#pz&T-@T2S_nD?mDFt_DAtUk4zGSVHzwv_*4%)5y4g+POM=shLPgb6uC)T zSXMB7o<4sD+Fo(WEU=%&VbY0#rcDv2#~zEp6DDXL(kXf@whk*)EXS%sHo5ikp<{O< zd|McJC^5#+s-G zf0w|Z?e`x@^geM$;M4tbUhLPRK2W`K0f)IXhrX}e#6gYGZYL%NnBPsOtZ(g!l>PN= z^bxf)Hcc_i3pY6Yp>BdMaIMc?9|Q95jXUqGyJ2w5lx)|JfRz&d-wh0*I63Ab z)pYVUv{pUSlfPwyqjmc?axTrl_j;6e@l2{$F2q}j*9wmtK<#zJa{NFl57tpbOF@%^gJ;0py0P{JS*n-|LEJ! zbqNM2V+LMywoul=dp7XlCdUiq^4=)A4Re>>?s&KTM^MiPBpgXq$F)v}!Yr>iXbJm9 zEB8kLMwtHSwrRED{El;5HJ>kFpd7p1*?%_>r`jSN7Y^W1tJB_Nm*a8T-ByRe_$1bw z=)H7Tih|b8xom}()-b@hUeYkNOB8Ou^{1Z?uF3aj9slSpCz z?`T%5J~p?Go_)Km3n@;MT3bSxaKOj9^muR}P*`QZ&e*-gnbda&5^SAu&QJbo2)8Ya zUHf>w%pncBg$%`}`HY}hRZu}#(H#c5mDSsPPeA`st$NqpIXHD?L_TKZCNz3D+Kivv z32n2N6EpiSL7ubuwUM}4nD6>tPaA`uP;y^X~h|1hp{xdM0lgL z&3E?3GDLoLI-h^_7ET4ac=|Y$;JC7MZ0k*J$bS`W_|TymdR>k+Su~Yl3HQ{cDyBr3 z9@2YzdE)G*hmGn>*&o2!0fR{IL|N#^?$O;g#yB+KMYWHz6~0OBIwQ~W3hKAVUscLE zjRg||v*o!-Fmd4E%PdU;=+n6+DSN;N+G)(LY**97j-)y!Ywsa!D)Y;>T)2t-i}u?a z2KQoXW>(>7co_El?kpSV=fIwtlD!;v??X?mHz%)6A}p3*gMcyO?5{o7q@I^SI%DGR z{0|qfWO&xFUA7t5w#47-y?dIVy8Czi&oyTly8OsTm|6#>t~#Gtt$j>TNmq2MSM7xH zSby=D(pNBhWcHJvYAZ~7DzF=V6M11~5Kv z5zuMG2>tgNHRFdvan(yeK%veLDWRK^wb6~;wb2?!mZNa)S2iEz-w!Mr5C8AE@*FHE z2m1ye(Z!m5mR2!S+CUjRcF~N&06yG5a`_lpaOo4@4*|_)T=$?kxpL43%REKdYCTn8I9o#_;0#tDAkZ$lB}U z1`W3i?p2t-fwnyvm7LkRFnLro=OLpYKY==P}J4w|yoX3SHJNvgk=${*IYLsgVKP`dZv z1iD29oEm}tR+=q;WM4zWkqr(uvo`2lQnBj)-U9leSsVn%JAOTQ&q2hZ4&x&gopd&f@~`~z z8>X|^trH#o-OCUrl6rWC^dhnEue28B>KSZp8Mh0meTB>S&)+#(XpRHT1@WtXe>dyi zedmL9UL3IJI?;Cg6)qn4_(w*M;aJiqljhMBj4~`&%wwWpiH4YVVw_)aK*DQ`_FzB z{C3Ufi$r=e_IEyh>TF$%V@Jl#PIGbM@a5VA%F4!=z^{C0);SCtjt8?e#y4Q@&gs$L zGP0PuAs0BkvAM6zIfDLt`|+ol*u$F06Vck4&?d6N`eHnmWq+ z|7v6BR5wl>@6Ne&RD+;izDXUc`X1Z3J`8YgUBqRtF25r@d&#ZmoTv;sRABa(RMqb1 zH?fO$)xM_E1UoqI(5$~;fc1j2hL3l-W14Vcdh^ONf>GDP%p53R z%=@kgzPtxmKiS)3hpEEx7ug;#urwv{B+?q@W4XH(JSt%D`L+l1*L7jxuDlncqzH~a zyS%5%hz&*(qN}^B4Pc;Guj{``E~FTU)mh$9f{9L_-mK&lm^e|H5TW}Wm$vTtJM+mM zS_Z5nK3+oEA4`DPDz4T zqq)qM&vA@oe6oet(Q^%|{7m>K^2G=$wL|Y!#fxz{g0_EPq!h+)eC5>BI02p4)fr=T ztzaa?c1}!S5@V9~{OnfEBXKpl$N=ep+C>!B+wiFdz{KXwo9+9}Ctj$?{k25o!^_)Ypj)?wW{FluIe z0G1ku}* zfYPAz^U9Y%?6!Gn=_9ruCRhYg{&s$YHNO4}HK(Vr;P2V(?#Uz2^RdG3u;p`@Nmk{5 zMztM&Xv&TL)MbU9^GSu*zoftny|08FTQm$gvGrD4=E2B;>wFO=kn(wO zOqnbBbJuXB+j&R8do}Fa`z+eFsSPM&wqgz}#w2x5WcIcX*I=S!8-4YqahTnrS^DVK zR{S#^_*;CN9Avh~+3HF3;IjGMf48r9L&bF24f66PC%oM06Mi@YI$TK~mBv3XeTnAo zh1F&pe6XUjuJa$XzkL#>rlNtA@IzTHo`X=z6Rs`jeH)g{^ecGUB(XW+=D$>-NSN}w z_}xCc5}F@uT&7WA#isJGeNxu*m@7nom4@yst|Zg{D;YS0E${Vs51p5W)q;8HGn)Rm z;5@3Wb|j1h$F7>@ z1#}&L7A3>!2U9wCJ_|NDV5_8`#5}7c%u?+RZHm7SWm&}r;fEe$uOeM_qS`Rd*rY3D z)lOhbnHuG-lLk~XKiiu9WD8EM$2AN;EQ7LwjP>6vo|q$}a85bDk)%HvrLV$z8G7GW z2~k`+u!;Av_Mub|Mh(GEuKFT&$8FZ@tt#l75{?@v9N=6 z!wb-V`UorpIbqSNmYqx(1qbHKS=re=v6p}8l=5Fk?5-4zVE1|iEe$Hia(WtJn9)?+ zd3XsjVi@FI?0&)2NbxVNso z)?2_WXz<#Wn*b{h`*kY}?N8N4oZREU!P2UYZ8bHt z=)VsKzA&6rYh8l;^H)!KL_LT6o}8V0OXW!6R||i6AscGA-3GS_c*DeyZPM+J=dgpn z|7B5o9h6n&@^i6I<6OP#53dhraP@m%O}LN%HZW=Wy6~KVuMc*g41F2~bt*P(DQ^|v zk9z4l$473sbhp!>C_xl9ek!TGc_0L{e|Nojwjhh^7ecqL8QEjT(9NC9!gA1)DE>Nnv9JyV%d&q!zw|;~R>Gd;a}_`FSeQkB|MK>f=zbf^j5N`Th-|3w%W~ zo(uE46l_OqwWj;jnLY$9|2ihawkbg6?#X(VjR5RsIz1a-yzoCykN)v%A2SDu=`0#W z7G5DRGMrZDU#N(AKIzM6&*x&k=hcY*eX=lC)b4+ddKu#6yWfA&N<<3dUZ6Xih|@wz zq4(9BA)lH?>ss{zC~nV_qR|h>xi43Iu3Zs!s7%6U0mo9j^R zAwv2;76wYy?#~Ea(^qFAbOPZtlr!ef^Zt>{Z@lc!@i8kE!qymQ?^-urj zzkySsk`ohR)VORM#O}f&4Qn?qi}yct#L=#KE_dtS^Q*v1)tzU)|YkY^>$G!jiNVdqny4+Er3; zfU{QO5)BjdH}TzcE#D3Ej`tjY<*UHJ3jes)zaALvd44>#nu1BEQl^>yIKikj`;*7} z>ToWmHJtl`6hWo`_^gYmFBV4ZtoEzDh@~nD%ERA8aK?nSIrcy)OjX>ye0r}Mrbq3S z6+dzcXQovKo+_WgS&QS_XMQqcdC_6{k1{OKmluEWi_LpT_Nsl*|K1xKEHCbn6Xd~R z|6|c|KWuQwkgI#!&;yuCJbOQ1!5T>2|Im z)V*qieMKK2_DX1k4ox&vIy6f9Tvf#6>FrCj>$k9Hg-$Dh^uT_T%g+TZcVNQ{Yk|36 zmvG>nS0z6KH!c^HN)L)}!G(%LiiBS!j`OHzU}76 zKhKWidbNwhzwlL-<$bUFKnSj@r3I1 zQ_aWDM8KQ`O}BWZ2h8Yj(CX0{V-5Yq?X26laC`tbK*zs{;X_>LbF6VuPKchE0E*W4 zTj3t}Ve?lo?&Fn&r7a}yfDIY_|d&edZ=0M%nkAU z4xjUPe6shez-hyc&d-6*aN6pHN50HDb_Fg3d^U2x)zqcv_y@N0;vI#|b z`7J<)15+38fvnB%J0ruf^Aap8wG3W+rG~vKy@`(e?{Vst@-#=8B#v3LarX{v!G@-U z7xBW4*f?-igRW!{ql=H8nk?DepT?`RcE6mVEBi~8F+&;*#^2BJtoRRSX%>3fwz~j6 zxq0?bo+C^=E#Gndwhb&@DXMkB1^n7K{@!UU2y2CUq8{s;VI9q{mEXn|*zF{oJ2x(l zV=vqU8z&neFNNWl$V>%}u3fzo`fvo>)x6ktN&STp-S12zwc=P1?a(eKyMmM@_t2Bu zUf@`hS6os37z|9>@z+<>AmxyA!fiWVSdjPdjh_nGtQYgXTlDGJ(*9xH>*^m&YWCv# zs`?Ly2Cn~B%9)3ytOHdOA|^PaArNA~b{3}3deVPvY$dm-9`zP-<%e-6kwJyZd$?fo z_{HTN7h&=4=|IUq0jT2*IK4$_046Tn;k(n^4|Ca4_O~yQ(Er2kIGt!ZE(R+<9yqE@Uj(&oShOA@ z1T(69*+O*@IHW?t;X_&2{Bu5Jfo>5)-oGenSKh!m7v~8Z-B&mjZ`t`Ya|j1ZB8q+P z)j`*#XUe_916XL~yREqB7B=tfS>YO=ff^TE51Cm$r0l+tyes7zF8tYl=lQl0SRa{# zhcC0R-D0)Gz$G1rn7%OYYDwFigAaL5S`xTGWA7p-qK-?EJDv;$ox|~h=x-~_n|xc} zdV9#^Gftit`F(-94cda;u1|Yh0m{~CJ{`r^xatGdH8zJ~p|Y9b^o{~teD`B};W2OQ zTJ^n9H6IH_U2iLxFI~imC3BB0zm2gu^8(wHWIoisR}g? z53q{8rQxW_1awTtf9_Pg4b!Ec`=xFYP@nxn<&AMNF5z!_!~8XZD)0!?oj1pEAl$%Z zXM!>;2FH6`8qLM7ssxetzr8rA^@Y*LuO2^Ef4JaEX2NWSoDlUzd1&a%Qa>X57uqlV z$qv*{gsBe0qFNby7>42i0TpLl5qY)nIj{nTl=clt=X@vWjup^*^a{c9AGbNBH^Vq3 z|Mi28x;~Da4$jFjctlXgf4cklvmf-?={@scoj^*nMZcs)1Fo6oJGy#~L9MN3(nP;5 zcJ1qSq}lI|9UE30cih}y&O~+2&bA7gzTIk>H$6zu{SiH^bkPF3nM-H$uSi4x+dB!s znGfEpNh$VB70^&C$x{94Jk*RuSg&23h7ps37tgz6q2b^g=ZDKTVX;r{-Onc`NV#`W z^7uj=EUTz}u{re|ChXs*mxXP|`Pmhs>yIaNTl^wpoR7k|OW8I09b0kgle#Xiq$P}` zPS_s&dmFpBmJ{V?>2S98d#g{d5q7CCtPMK`;j|4SMcj)UM>X=yzd1F))Q|T|LKcIt zn$l-HANmI=?d}0Q;q37J)p03kNYvtc4*9aEFj;b z4NYBdAH^4GAw{oK`?Z@EPW@V8u!yOIXYI3kDP%=R@u;)ZDFk5nujySyq(9YJ8{I~JiFv@9+H1b zM;PPEaP_s*wXmz^(0|yv@bMdSSU+Sj&m7o`Yfm1gbe{SSvzk{*8T1>WWZ^(s`tWY( zvLA5x!E*!Z*qxU3v?B;A^{nXjY<(QI@V3sm@ex~G9(xOIkB2gK*5po8Mx54R{!F!- z38zmI6q8?6Fv`&PvcTvk^gK&mtyH;;10GfGW8dmp3Uaa+)$jn-bgT;cVlU1BXI9%u99+p}K^9Fn+ zfzRV{jQ;lLm?;ii=VQ$#oPBWo{w|I)7e#RXGDmt+pf6IMTXSylAKc7;@f|-0nqPkT)j0ct)rsN~uE#Qj9AZuFrZRqi{f2$im1KpYncRPfIvG=od zeKS`!ye(P3cS}_n*W-_$4rm&LIzz_23+MA-pskFJdN(~ZPA=X1bdnu*%XdiW99hIipG#Weh&y7CEBT;^`+S386WC7m*=nJzFTT})$idpB(Sm-|v;_BTNj zJG~=+=o_pEP#9{{1hJm>es>w24HQlXm%2TifRdc8i=I}(IFMq_bxU0lDXSc_(R!b7 zg5~4{J)C1jYr!v5Y&yl1={#)_mZ1HA@$VI$yCP zC!T$Rf5KDxnip?j(eH=-7cI}hPpV;)e}-b1v-3Z32R(nR&z|PDzNv(Z}r~j&H*@+x<0?51wGr+OzYuJM&;_n`+ACTQoSz;(O^C z>toE9w#^D*<%RU6aB;)0_1NX^DQy{b68iP`tIzy2gqm&F2XCzpKosTRh46?`Tg>5B*oA z7B#G>q3hlOw(BbwVWucEJMeWOWT))8BdVnY{p^fq#e8XTDx<0JlhIimq;!hx(mju@ zH}vH-w+BNr)$S2)XJHt>_xp+$>!xpM(U=()ea1QtkB?emZ(uInhiV5W4fG$%Er~kg zisN$_!nD9AIad@Q}UwX1x|Q~#3&~TLSEI` zt2M`dz`%Eb3qRyEVElCVRa36}P@-HP_&0ewNtKqLDbQmMy|*Y1=JflqtH+?Qg8e^e zT|5=cIS_`eeWwlv{xrgAs_>Mlf{PIFl!YT|bM7xGj(u(my9}6cYoRd40X`2cv4$-f z;`BGa;@so81ikz7{dVhaFiGpaoUqA3_kE7t(Oa>{Ihg{pug|jn=jokEGPxbTrr56Y z^N^Xz6_{l>I58Y}5?e3-=lkTLIl3<%jAMHCjiCMb*5FLs1tg#C87G?9a5}-J-?VxX z`UG^N1pAM{P|B--Kz}v7!={*aFGc`Yo;V7>=~;&9Sn8{zMp^jpNy&1Eg*BvwIldWt zl?y{1>Jkb)(lAAvuOcmy0rQVPneV!_9nzmQeVNsy#Wmr#RA*Z9pk(}4Y0|SUSQHlf zJhQ=rk;lzd=o_zL?%gc3Om^OBUaY-!Av{Nt_Zm@ z9J7+yf}7TH#&X!UkYNzKR}-?sKJ7442*!ox9Zo$1>@ct5`X9%hN6`E< zwet97VVty3Q#(L2iWHA;K0$0Bp>JgB%dNV-&>JWlATl$BUH1gzPn>Lr`N7Dey_K}s z$uH~|MA3z~X}1y$$#+=0d$2dALjx#lm7jyWenY3`qchN65i+8iX?1qgBnGO)NmCAjsVSP|AHFqfB&0fsW zFEcD^*bB1*uV2^KgyWBK)fTsgc4%1558+{dgWaPnbrpNsU@Yu1>x*qCp=K(YKKY_M zEdC0k8Z5Ac;@m4bdvm0bLQd!T2;PP64_w>CbY^hsOh4;SjbWTL{!sVEN*!C1PpyV| zltKyPf9)ElRADyMkx^z^8>T~z7elyJVc^2b;jJfW@ORp*+2`?ea;xl|pxom{?7PV> zn4`Y{(L6G4!KY%eTSD|EAL$4Gj6*p-(yzi=y2oiX;p5P=-H-bn7iBYF#b_OrkHBny zCC{6?&!EV}T(J7pZ|pXj*Id^x#F(8Yli%DD!SOr2#@(hLVc;Zh{L5pfpg+53YzyTt ze9^i7Hfe2hUg&n7y?I>`#tOuP)PyXsPDGE5gZeBkCMea@8;@h}SJfzngPbrVAQn&{ zlntV+|M4I6dXxyYT7~r$ zTS8#;%^RL^c6u1F;oP3#`V{*gD~rqYi{ZRG)BPQDZ(vaLp6lxGXE+nl67~q|vFA_X zkI0r%NH#C|ns?k2hkI*xPp)-gZTWHWb1_QTa^_0t`;A9fzIH>M;ix%ogmRSadGChl zORKVv)VpB3eDUzJBu;ER61vq4zQTylwmYA6rGXOr^(c*0JT{wXYp+KpLTmBd&EByL zobser%)81CYjg%OJ)^_8B)*VZoYDt<;=Ry37DH0GJrp_8ssV+%(NV?e*_)n|{Nhon zD6HIedvDpUjkzC~cA40du$rBuIH+z0^;F6N`QCX@b?@->#JmkIoYA92{hv7F>zVHr zcMXdBqRRGdzYSAQVtbVG$1(4UJY&E|Dy+-fMLBbp7Ym$=f4{Yy!zP>jSm&fFf=Vr7 zH>qZbeIt)G_iS=;tO#VbZ?#48o!#%`Us@!!%FH$AyfNJH7`prM_dhI_p0Ixo88E%d zdDe+G6^5TrJZfrff?13BZ4Q5@ut_-bYw?>LoSZ7QvA^^fn^H+d9bPJ8OYy`}WIs0mk$lwV67n8&`-o3G};5eF-$%f5d22@NKzQ9jYL z*u7!>oGWhz2JF9-X}#DBJsk8*o57X5M5RDpGLK?eTrvW|Uo zzZGUDUKL0$1w(SB+_nDUA?%Gi$RC$rh_hd#SFiuNiA^D=EjinKaH&N-g-+fZN3tZ! znk56^-@Z2fp0less7m~be7fm(%M;$kPv7Do@0|Bu_C#20t;iSHMT50lsbm`2d0>K< zv5|jw7gmSeIgZ*w^QAay{z?^+>XXoOc3xG2J}K&u!9)m@Vs-n?fZbz)kp#jvWqHpZ(qvaG0R_Ax58gXAchE*XiC9!3Co?hFq3kib40e z=eKr!%ZC0XqM-dNA9ls6%^j%^1PssAJr$vVWx7vYPS`3!AEmP?`M-SZ?K1Cisb7J) zVvb#}%jBWs=%3qf%gb>YbxU;~Ji>`c#|L~d;;_z|`6mCrS2)ir9=qzk1Pii~F|kj2 zv2feBmRS))9QH`|lKpQNOeAzvde#OVOT!g z`}?0x>{SHmNwR$>_AtoPuq0Mm!Ii zWuYk9p>wB3BrF|1))$(62^YQVO1j?4KtHp$?Xxa67&AELs&<6|D-J~M9TfD#rFOI2 z6aBpewZ_%>XNR7^j3>K|^FluMGpWBgAt!|+GJhY$(m%#N-)jH)`b8LXYWzBP{55ur zuIXn;iDR!un4Y`OSLkV*udRHb45=$y=6-8Ggw+9V3mcrpkym?Xi{3fH;_XX=hl&Bm zP8d8_IzB;CGahX1E4+rw*>kTe%O`Q>d(?e%+C^-N>ti`za|F8=f^&GvX0eL;)(wuv z@6dAK{^b!XMkrb@m5RHN4s}~KDh<{2uyvwCTxmxVw9{=^FFsC0^ilBUEgy!u)ld5j z|N27@8DbfwAAvQ^)&oC3oP>c~Oc-ta2X)(|7rm{hU{L?OvDD>hta*5*>G~!Ij_|s@ z>~>{>KjzKJvj@wd-*e6G=L-v1)*9Nu-#@pdoIZSn{6a^M zOF?Y$hmDPJ0}?VbvwiUGC568Da7%Dq`a7K7>nVHHT?M+||IHLS6NF=92Ctv= z{Rh3YT8}CNDF65W=iE0MG*Il4J0`XWpT%0=KSEU}y>Wguy82?55%yj$8hJ#u z4~7_1XAd2u#o>;9>0yb?*gG?A&6U3dvtgkNH%@QCC8`+B7M(~Kp+2|nVB`puO&4Pd zuk*mP@5C#|bw=!Qcp-eJ$O5V^ZWq55NDI^X-|`d(l%Z*t=7-^kDd=T){%j%`+5>LE|x zVg7_s#7(^=9ByPka3W>1-de@PudaCDs%>w1#l?9T8E4s8)@g*=%MS}v8mOV2%5h&D zp9+rFoUwYxw-bjrB`c|^i(t0hrr@wJEeyX5tYPoyhVjS0+`>FXKNuPAh=hr6%T6b}ZD3W+ zrTlMv8FnAvBAT34fD=|3z9K5AknC?1b5%(Z+QY({^sTdTiT+cTszfig*WTsH`{{(s z%!(@<%#sk*@ZW-?p)7oNjd^S*)D9zyS0AP=8)2i9y=!w$A#^0{GvYUSNU%8hat?eu zfm7$V>xY$v;Izgex@V=H&?_iDuG~ou|9|M=WSr>W!KD+acHyj!wzYDv9MscS z-LdRC18qd`@Up;vFg>VsEAg5d_SenlCmmIUstdK@jNAEP%JatXV`T{(_^RV0@@2Em zztp-!N{3?Ywj*vow_e4@xr{n{j;~nm)Bbs<^!&$5Jg{_WDE3D8VVDW2G`v?uA!*;3(41|^!RchN{V!@upvCghEtZqt2&xBc zzkjPNU`w2F1C3@lG>+NqDEF$zWxHp~=TvS$TEplSe}4i}#vSUm^8SK?Pl_yU?qeAD z#IG)ql@@11A3ciyd>9ug>Xb~00|X6=^szUenIT!!QNrQw792F#HuK@(riXtzsupTr z1}pso`#Sbn!8q0T17}{cV!!a8OYxF7Hv5p};{xILKTmHtU}D|((G`3=@2N}$9fsk} z)|7e@k3FIV8L!nUu)uY}cwvhNLB;PN`RVq4DE!!$QfBf77kJkF1+9{?t<)gNQsx@W zi^b;2*(AWE-7&Wv*En#(IMd)|$#I}qvrg6=)WoXYSQ;<5`S~fC>|eCE!)Sf41vlGs zn5mOJt5xq0^9=k|aTl|pMImeCF1-S|wR&tc?UpnSRwx+y2%m#YG3(H+lABr6=iS!t1;b!jj z?AF-KMegTHi_Bd>8NAiECNYlo`kHI%;zT<1!01@ zI<)$p2TraZp0D2QqnX3}-D&UyXRd8s-{E41ZDtd@YbHNoa>aYUGsFIyobp%sT52RV zS~FMIAT9LzXC<1cx?z&#ryU+SepqpM@~Ea|1(LTTi1xw)7`7%Or5kPoWk>PRSMi!K zJ+5`-;PHGcIN5nMHi;K<>F(ED{40i2y`dX|`a6O0#K~j-fnzWe={~trH5EHO*ha6# zc;b@f(@aJUK`3%LawYLy67))FREOKN;HsPUmFzPYV8J2?9x=;6y-(J@tp(Fi2YJSe z4LX?GT2c4&mLHVtV0u>Fq5;kHHSNqB0$9l*Gk*8=ckEOnYVPgwh2a~re0B2=ps{D8 zo3UvLD+7z%_Y;6K$!zcZgI{6&ySs9Q_gNreP508xwn-cb?&#f7cM{6a{Lpi0xC{R> z>0$??|6$kXpr}U!JTQV5#_ZK`IN@@Z>R(GWcIv6cD2_D1c=f7jMuIFSm? zr@!}idyPOcN--Gq4kcs=w$X93~&22`sjEPL471TG5ePa)CiQ9 zls)#xuUqnV-2BlGGvfuO{_j4KTb67^cX;o>(YJn8uitOS`5n?(a%Cm(OvE9n?;k6y z>MUucECp{@VQd=dEL6y3@2~j@!<;YD>@7lZu8ajsh`TV;qZZOB{2nIW-hE#8 z;V||-)RDdYb{q#~YpcE{hTuR)!}phJ;aIb3`lE@58;iEC=x^;%!g-ernmh8tI5e?; z=K6JS=sIi9EN&rH9fEeT$*kO3)5|6K8=Dsn(bw1{D~+J5h7KxC)z( znElK#tiZV|!DI@B8mqbNdW^+4aO%m19=}^5POqt*+WK=2r(eu|Fx(Xnluz^*+F$xW z@BO>lOu7X02Qm)D+kQDdP%OciwhWCqFH*Ku zJ;VwF_bgm{g$?uJFK?*vWB0TEV!h~XKyi}~YhRNlsl{x3tTRNh$2i*Q*Vo%PSd(;~ z_ty!WHGV}*O#g!cuBS%K@t2^syfnX)njU*^h3(|#9)vleB8z_=eYl`PdT+Y0c9*^C!d%K_6B~~xJM&{rN+zm%U^7A3(K0x4en-ws)tju_Yo%y zDqbuyli`Co?%ND38LH5K`0=FrkwmDH64Oy)ro|Pez*8@nPvdClz^-R-Lqi|%a?AL51wAqQKn!ucyP_^ zd5cB!I}85Sm$a&<(9)1 z?I+Ga=Y>1RWed4Kflf|7Ce&7<`BZ%$Cz zc4$CFmlLXYJA|8FU&mPMXE`s@yfD8nHFdPi2?tH(()X?GhVDcEf@8azAg5q$O0_ir z1DyvdR=>7DFB_X*tfm8w#M%rrzjTI~iO&sjA1yX{Ojpu&2Rn?2eWd$Xq6K~Gw%_f| z<}mP5d<##@C7iGQ;+>+}L2ijMlH>XL1bSkx-Eq)d-rQ#ed-G9*jNzd;ouh5o@91_= z>|74^C`LJK?d8X|%K_3$ML%(3E3+yk0>%975_RayryhH}5 z$u%R zEBWvb^7$L}?&I$W+rV*wO>c~NRd@eq29)HCh9teahlLfq`N2&g(5L8Q3w$-u@@3#$<1M=?BncIE>k4enKVba3`sDeyud%qh2lcM&!rWKdL#uT@Kwh=K zH?@-i+Fuk*(;FHQ#GhKW-0OH{%KtB z841wgR)?i0OOvJte4yo=c4qOFP-tvrk4vIH4U|;h1G2OASlxQ&jK=|s&HM6YpLxpw zP4h$d_Wxysj#ztMz*e_fTA z`lS#yB>$_(`0R&M9|qpMO*}%%kTB+Vh%}Y@IBBrwZ;xSx1{lihl;w_9k z5jt>1Itx~-+27KpRA6^Q#>U+X(y%c*Jsp@?hyBmi{kHXdzDtz z*$1NzjhDI{SaG_hMA~`tdi%df<*x5x!UYD62f4p!U~%j1d(qf|-6Fg;0(G1)^dib9 zqTvvXUJeS^S9Qfs+Dpy!c73qnZZL7ocM1D4EE?!*#9*Sqnfdwi3LFvo7{GA!3rQ=; z=A9AWMbfH`vlqTB$03QQd&Q*!aQ=DM>0R|DNV&I*<){@aH0wQP%vhs`nPQf`? z62|n{IU)$wPkr5~)Gdh1dS6mo$JH?JoOz>V$plXR{;2eq^DUIzOCyf=Za~+K)~GGz zg)n_t`6jplK}G#a>IIt$6v(nCsmoTuO6jt|y}uDK{l39c-S;-WmYMin)tCq)vss-Q z4sWsM-#0ggO&%ZEa?|pnb^^BOZN1^T?-iD<@#P(>41~e$Zm2+2js>S#haK<8LEja9 zH}(s^VBklg!$W3uoHM(eAc@yuG)Hv%^F4>4b*lBnv*2Ou6Z=7gd{)Or0k(1Xr~cR( zU0ss;=??Z8hBh4e?E~Wk43+cD*`2#vWWujaI`;>xv?z83E^aqQKvuN|i{V9?6=GR;JojHnSu{nv4P%h@T~h$=#pE?CuSYt z$I4EM)JPZh4^#KOdUO(&ngzZ*R_lRU)@PrQWiuzeL!Wj~c+ks2=yj*E3oOd?{AAWo z#@^fXd0{;Juz=&owrA|OVYuYwA3mLC=u)YTbEmO|ZVOT7!<&9KvvXH-#g9ig!!4Kj z({mIWE}eFiP}<~%g4Lx|fixV5FaQ1Hnl*ecUHGM-(+`D0{vp#vR!|oumnLyE5v%oz zW#7kmLGfU{F89v81XcZ8<7e$-Fvc~Z`>gUabm#Br|1a1G$KzJsvxFBy`cML>yIep@ z#I}>?RjhD2{bx%S!XGpvE&kg7Z8K}LJ`kbUTTWRZ}Qh`qPV{-XfzoAb`D5Ki) z0jxOG1pk`ehfNE#$3GM};z)zm($Z!hq}rZZzMiGKiD+^AKW=ha`;VdPUhe{Mkez#l z>g6m9OU{$j4` z2P}0T#)T5EMHP}CX1(%KnBR0lyW`in^h>cY=u#9@VBU>`yLVR`OCN#$@R{gUv7@l` zXi#y~#uJJ(pP20ZvICcDEh;|AXpmHPO7~W)EV1x3%cXX^3}WjoT_IloK%7`3!Uh!l zpzCN7uj$QVoY%gt?SIc5N>`|!xJ3KHw}Z}t?+@8y?g=^O%kuUxNTYl{y1El5uXU+k zn+_$Z8i=!-R7_CsZI#swOv8Dn{GhD8@i4slb|dkNEY{4t56EnQIpR^YorY5{=qWTnQi0l@%PyGTN`F|zgjmIT zG)Rjg4At)SUV|3)&^IvG^|iML+js8ZJ{$WE>r|{CdjHzYjqt3~&C-R?LswH*JM!Is4edT$?~`O(u4)8_|+2ZhbCRLs%S zvcdtAnM2OJppU@0?-Hwfj7p&A`&n%}kG&X^I)3B4Mg(^LZaXrm4%i;ywz6+;8BU$A zGAmBXg9%rilyx^|94TqMZm(cM&`NMQ|K8C7W3dv=)-_zP*qNd9DIyMrU-oUCy(op_ ze$(c5g-n>j+o}2eUkuhM8jQN+)L>A8?8$V)L$Dan(PEqVkDyw0>8xd2g1Qi|5j!cw z;rG-B9$(CaNn_;WPjkd#r#Sh;-+bXyrBgj|cNxn5dzrVVTorqzMfk6Czk{!D>EFL8 zV1S8ZY*Tyht&`Lx7oEQfacth-m8U&h%Ax76;klTML?{ag2+X_`j^k!qKe@iXiNE8y zJv!b!z!p(2o~^%~V6`8bi*=&HUs-1}r4#&4ax^#(8DR9pU%PmnuKxg4fvQY^vc z=LUl;fj(INvww(3G5~5Bb>*@uWMQO!l{=z07nU?ysLx~%K)<4GOP`@Ij^0W6uV={> z#*d}59Sch*8SMBCw14Vh<*(0Ns%$E_wBUFsr0gxs?c8-Lq7FtFYRoMxOouQx`H*=czR&)z@5SOUd=EOs(x{|RU+1#zo))UZI3n6Xus+lgj4 z-|GLIZqxuNPW0{>X2vjNx##fR9B~*sxIZc2cN)x|&^suvt_7{07sWqQm*KSiO(%;Eee#$tN&+ZijJ=AT>@(e-J3OybZM@0cWk9|AXq4Y#!#YpIF*-@Dc4XekdF{ z_|4sJJ5D?>q!BH;2A!9E=C<6tMA8j*1|FL)!;oL)QrDa>;uP&j%X3Q-S6}QusNkdu zBbp$}IKlRRzjJwK&L9k(_kXSNTHS}4HRa>3R#iBb^5^(Ai3AvV6~9Iu)QMB2JOde- zG%(gHy?;-KYU!kvyiU7o1~jdSQ% zh`R$dcRuk3+*yV8&(5D(&p*SNW#bcvOD%9(`R5s$_`5ielsE0g9}1n;uk;yCe1r+x zYjsWY_cnRuYb)LN70h1v`EU5EF1AT)us!_#8dl^qYr-VvFgioE-t3nnwuCCSxINf` zmFl}aTWWH!!;#BBNLC3Ko^!^CH#TA29$UQ+3L`KaB7F34z;2k+V6JZT8piQx$5|Ko zDV#XQbNszQ61F!s-F(rBOCxT~Z07#A`^~9!# z=>KC%vFrQ;lFn^|ZhpWO79QTFT>4x@(#{Rq8#H`|XUn4fM=Uh4-%?B3<~>4#*lu6m zq+F!@?p1zqWOF?RhVQiZGCBx2QMabjfX7{8!izCYm(k2;7xO7U!oSL-=>MyJ1 z9JwKjeF~*;QRV~uWGJmZxcD7%LhR36F6+j4(S;o4-8}?#Qb=jnHkHl#x$!kJmaU*MS$5FEYW2A90;<jB@LUTn9MBivmvjJH zcrJ=x;HSW^zBAzpgd$87{i%LF;{u}$y`lOG{5Y+q*;*oU6Z-$NHg`?#$JDpej8+|c zVd@8OaZPPM_V`7gkXn2JvmyB(+AszC1PymiUWkVwby{=ocdJDYa3D z&e)k_@p*RgEz2}6y?Fh_Z~GEX-BXGdrO3fZxW`c6&+|x`64d;*5d_7>Vavxk++bcw z(P3z(7E-#$1*;U_!_whIh0M2y;DL#A(d$o2nA3ZXZ(Ca@^jx3F_np&#;tf$f(Q;v2 zS3GFix7>*Jnjy0(9yg)=m*8>Dw>j8Y{Pu9z%q-3yE@fex7lt{fC(05(Jb>czGeIa< z6>Ae2=yn??Z`Seh_Vo)5*y5d+q|#0Y9ZxP?e-o~Tm5oqpW8{w$A)op7R9eE&PKDYYaAergTVh)yG<^uY{Ee0YC%Y%#L_|nH!P(J`|HS=)qU*Xo|MNc_o9UJsTK31x ziqP~$MH}o9&ty^&9EUm@(tJsP3MTjLTQ=@K2KCz0A6=V49PwwyF(dt#9U$pUVEw^h(2qI+PJ>qZxi=uBO4jw$)aH`P~p!c>BM` zI18BNcewFRZynne%Zgf=1fU`0ma3QTVTJROJed>cE1}(;7(crEDR%Znzmxqsg0-)X96MO04^?|kol%_~fQDOQ_abik zLf^yrhnMIRu>0uJKn7DO)V7(rPrK0IkT-d|d}arqv5mrQJ~x| zez^Aj8MdGHZYcIOgI-I;h1GjQP&)RH-y_EcmiGy{%LG@!XqsC8PU;9~8uV)nt-TJ@ zq71Ay4@+^_NcRcL&qZv!_j!-cgWr(Cd!1A1J|8yiq~TiDkiou`_$$Hly*TvWuNnC; zZEW>z;ThZU9frEe|rrF2P=QZpT0OD1YL$- ze=5}+!;!H^=Zdag!tq`%GuC4daHeOl(!GQmS|@#Fu8QkoLD@%*v&q|_|FqTh`u}X9 zzH=b%!MX>mCUv8$qKX+>@ zIv5Cfr+HIv18O!S;^VLVFrV5nbS+sMx-LA(Q<3n3u`>IMBQf7#UF*hh%$=RE$~&uc zoyvujKVP)^IVwc^&rVz?j5J^h@b~BsJSV z+T%Q_kUdNPrGQZu(oMYBeX}>9%i|60?QOa+*)8pN;bt8!yCi}9D{frwbaoI5n%M0B z2>YBTX}FxcuyHF?4*E#$Z*<(>pd}}P|7de0CT`COfM&sxr{?NcOe(#gq6#SZU`_&JOIe zO%Gf=%1+YspI6}++y|6a;{5$TfQ7xiZ8wLH67**Vin;#Wgkmj4Pf_=eNZxKRYn8|+ z>4ui1nUj-os)Db&&FVOg+&)CFwJL^vpIQ|@f9Zk1h`5TcHcQZrzL&l}{SE!pd}NUB zL#!>W3I6ij9B1Ai3V&pG6*@lCEhKgYZ)zxIG1intVq6t_n*hpj?l(JoWq%Q zC6{VwEMfTkWGxGaULI1K@8pE`fnEOc563p^g{L>gg#l(8Zk!((?Z=v2jsrUdKEvoA zC8dj?BZEla+rLN|Hi{^u6%ID zG0IgN!IxDy;KbTRExm*zwL-kAD}p%ppx4{~{7#(F>YLj|9S5rsd69BD`2_WW8{dD% zPeL_QwLJTYgV@4Y^kI8sJ4~NlS+86Pfa!m(8S54auqY!vpHs?Kw$p-??Q>64ltv9&X6wl!Ig3HY`f6MmQFeSHSc{87gAE z27Y|8$0_IDhv!yK0>yO8o-=*yu&TfxcR4GDpvpTlP^h8@Q%vFRw>Ld@+QK2PrGN##ZT|u<3lpib*vLb((W@W<_IBADtRo;7g3md>5y7Qw_su3`BeD zn{iqtR9XC`3rgS2ILo-paq82{ z85(js%rKbH9-L>zj!xG{&O-;FZO`eUN-q=4c*R&RS+WajMegdy$n3&Nr50`pKZviMB2GatrH>id);{Ro<)4cvnobwKHfxHX&o3g=&R(XW>**ImwZMqat2hk%Np`*kf?5salKyd=i`JX&+#02F+GfC z{opoVQhYqmjkXjbX!f@*( zMEhMe*6o{zp+Ae!3ZE))^hU3x=HUdWn~>8Ua`gp@)8K8tP%`vKa+{5c9EbS)e{bW}xg&(5e|Bq|8=b z+Tl_Ho#r)jO^$e2_-(Zmdbb?@x_n4T+NXmQcZTm#k?A}2p1Eyhu({xDe(@OmBz!qrRefhiTaX(yq!!qGOE`*NrTbxGo8evxS;F9p) z3qbpPJwN$91qWj#e~LW4j5EIF!EgAF!o1)knt*s4&c-y;3)(#+FvlxBb1b`t%Nnb5 z$#283Q{v=b?_+21?~2;VJtgwyRm6+6py7@Cscf@Tnnzx!<176)FS&3tX-4- z6T)~JlbWZc8DEEDqv)O8hp3kzKk?u{X2}F7;u3Qe(AbXs3Ofzk7>ytRQ$Vc0YT7i- zfpHnAEw)#){z%|lK(XYqcRfsbwcj|Ecpcjc{uQUDT*C!1$-v5koBQG9GEeY7Htc`% zj#F=o5zzFxsxPDxq4$Eq9Q9@|mcRZmPG|3lE$uS8zQTyb)%CQ83Rj@5s;6P_jOFHD z+go^b+ak1?{E;lRI*ie}hgl0PhoOkNsO}RJ1k~xVie>gyh~YFox`kU1$@G_!ED6y# z=z2=}>nT>4G?3hTJmLaOo)aUJisG?!oPB(Mx&zGqdzPLdvH&Sr^KN2FvPdg+_VT0& zLHgIRAKQ&gko?$IzpC60yW}F*4=l2xp^>LUj>2yY_$3;C-A)>Q#akq9I39xe|CG#x z-_&8>=pKUZ3BZ}kl|y{@r!n-IRePsa4ix=I+IIUxDfSG>Z0VQnz?qrU!jUsnEc~Va zO2_*Kw(#ALGp_v)opoQ-lt|rzi4#wrp4Q`nu_He`wyQ=0SvNP3uGIla=fa+nTnnI6 zS?NoqnKv$I`c%D}V}d{Xlt1(aW&`O8+sw`H>o}yBRyO$B3A$IgGT3A^VX9VmXK2;` z)oJN>L-&0G<1Us50pxu|I?lE6(-!(j+o{){ztcD28&w}_Dp ziAj3f1J%7A;gX(3#X;@UIA3tks>$>e_CI}c^V^#q92R5CIv~`HQ(l*wf>Z=xNXy$u z?%j82iW2|1OHLSOUn=$_zV?Atu4odM@>3Y->V=9!Y~aQ;d$TGJYP;S;yJ-0{#ZBCMdd8Qc2yQu%I?KujQ@)VZy1s`zac)T}rrpxA> zBqdI^b3sqpZgb8d( z3+6@W zlCQBO{IE{${0OF`n7NLXj9}duU6yaH4)jk1izv*yWB(EEQk&h>J->rN3yvlo8WiO0~KRYSHG*mN z%$nd2A7SF(^~KW)vM}X-aP;N$GgzWu8#vVV4nA_2g4$6nj||qYrC9Ssm8;c^2j($>d+eZgvZgswiPZ_utq&zKl2%ptc}6Kexs+Wd_vg{n6uZGXzruv0wlrnO5tP8ncYG}u8=#760PItqb4C*-!H zwjcI*tUp;7j>M%mp6^dPTEpCaBkqXXBbd-6arbff0Mu3TC>pt&;k;^qoYcooEM583zOLB=Re@olK0*zqCSqv{=Iy$^cI9}-|F!%=ye&k2$$_K(Cx9LA!+ zy4|`v$w19Kv>$FH3FHY}2ay*P~9hT%eU88a?LrXC&qgJN@*2h(3&z0u_G5}iFpfaUDxaQXb)k+&36iBCF%&+U3 z+$%jIU7vch_fdYVroQ>+fje~8&RiJVmN)EoolJ*Tc_UH%&+|AH^LC2nmcC;zjc_8IE=MMqP8;~L0y4d>aIgYF}DYpo1hk^MEGk)6Vpi;rbvEt@3w53h( z$u}n9;Npw-O#8HPcr?m~dwUTMNRh;!$)3jH+g&v_`{QB2QQkRe_A`tgjQl8UyMwUB zFza%uz%`g`CR{&1U4p{~bKj{8SvWkht&*etKK@Z4^n6=SMA~T9=s+SFMt(8vQWp3L zYr6`xBrj}o`)`Ad5C4wf{DqZ4@!4U>kaV9q!f^rz98Ux?%!-hZ61> zaNp$U-ZKL!QaJpec}LfWqZnlSYsZss4{&fPb;#g|DlW%onM*P#!}QZ@o;ia*U_`Ks z&Aq-GDe0~)#~u*jCu;^hEzcC^PCRGQZFMIw{wMQ`LG=d~H9}E~pECUR*x{Y0^%DC; zUifkOKgKFasil9rn~3yLV~$n}jZjbDVw@~bnBvU0Z5M0&%fY|{JC#4 z$>P}H>k$8OaT15ED_Cy^mSUIY!gt=c?l7FWTAx;)f>g3Z!na_42ohh4%i+k!xsSoS zo}KxGxL`)3_hp0QLoo1UTzf;|1}tu?8|K-25gMAm9vJL5h0!-=Jqg1P zq26IPm-+)fsJFe~!*lu!QgnV(lMIT9bXu}bGgG(^OaXpb* zW&=aQ8khV&i9)AVQpf3ArAV1iw{p%P;T+SG+t*Z5pnO78;pz6vIDG3t-WU5RXsO}Y z-hJ;Hk{%mVI~P##m40Gfm>)!VZIO z^Sf8daO_~4q-(Dg%)XH@x9LwpQb_E?{cDC`EdM=yS~)_Mc4vXqzkrrQ2KDT ztx@xq?=_r^yqjjW-2_K-s_y-GmWK~r3;3^u1R`yG?dF}hyD(lIAuiv?33I30!o&lz zVR7RrW64wjWVrmi`-o}>J*PG%?_m%UpB02_)J))tLaK&XxdCk4P@TQ=u?$-dK2CFq zX#&zzq#$STX&9X7rpxa)Lb6_d^Cuk^pcVG9DYRXNL8l7lQW+~4VG#Ms|C|KF{vJ=F zZ`xt?*M$lXrC?ao`241Upbf`Q%E zHg@bwSVi^Wd^8Zf$xZ7bQup|YTk9s6J{~ZK70anAn|FJOEMq<2-NiVdku+ud@n|TH za~Gc&C=-Pa?ww&DkN$?~@G6P2oU=qaJ?%H!@+V;C-F9QM`-iY8H`VI>5*cc)j123Z z^~4EFEvu#i30yM&S^417Em&dj$rUlUkFmC*6V&`|SoEL!i)Zv#U_!B)BdJ&v|GxY2 zXx(5dc4!@c;bZ#=dOrI4`}OeR5FHZ^_XtAEJI3a4lbbLi5$e=D%7^(=c1OPQ_``x1 z`CEehF(_sIa<*_r8-4$iElB9Shhr3jnh(o~I4xFM@=WwOmX%4%D*LZOC2P*khqB|? zag{gm?Y3E@cHHzzcsUE*ur8$?p^p90;mVo^J1`-vIPdYLWE}tJGXKxo4C_y3l00(v zVOF84F^TOZkgS_jlU!>tq|j;U#!xCytMhh$I+p^yrrRt6LLMS%!PfXC=jQnup86k| z8N`8lGdA|`H?XS}k5mP?V0%tV`Qf!&Fs%2I=u>$C2jibmM(<`rlG!e4jk97nU2S*u z^Un}$--z9H%qj&N=bbXPB%4F~F^;eCdxmidt{uAm#0%$iwB~oJZNd77yQ{v@a&V^J z@P)UvA^!#^nX(skYc&kkmU&@#y73GJ`t~4$-q!M&83t!DB{_yX9bU_{Hux!hY=2 zFmyOeo`ddO?phXy8!%I}EoJwKM>w6R8~VLo1AlGYElobEgHy+sI(8ZezzS7iNaj=* zG{+T`K;dr4T`&8)!)+H#)(}r&KZRkzpe((J*|uQzuZc0h4Gm)yj`E;qTbu zsEK_IlK%^d`l$8+7DY;>ENO{IRE$|0GctnJgzgpFR7*_%Ay@v!bsA}9R(S;hmN-(^)+r>+}uaALd^%AowIP^jXqP;<0|kPW=wCQ zE5yF1gT1~1`(R3Pa8H~K0jkOt z9XOty0!r~#vuA-?u#&fMK-04WhV#|phdyWnxl(Oj?$2|qx?*}K*wzV$pDk^?EZ*cZ z)jAI&o++H&>K=8=a2lHRdZr{$6DHpMj&!&xi|fZ1YuNa1!&v^b5o3oM3>B>qehcv6 zAdjD?)wNqtxaV0|YVYP;Gx)5y#%;j?IjwTlVqwU((*N%c9T$Od{H17v%NkBAF~*R7 zIbrs%pRw)@BQPNF{`E2c|FGXt<1(|SBMvENnvTq|AWhn_G~7WAn(hY6MsLKxK;GAR zSNb;uy6(LZzLx1QZolKh@31gvNc|<-`aBFWwoO?bGe#_AJ)NB=dl~yhD{Bw5@4`lr z^Iz5Z8*t$CXX6ik+R!KKG9MoNoj`AuY3M)Jf}?Fo5;v|&!hp4KCT|T7bZz(lZ>Dk% z=jv_fpGMz@H%Dm^CzN?N>+xVbG9VE;5?Og0h{FWBLOXpgXAYd=J#j`%?f??6$h<0& zNQ0r5E!O;Z8-co%Zf8Dp9m~2$hm3xm`qBQKWDfFn@^p3W2c=1(;Ck#Btq(;lbMiy`pYQ- zoyvh2`i)~y`TS^fOiee^y8O!-+MW`bmahohkK8~~-$*p;u7fb6`pw+rX%h?)B#!Ua zz6$>WUQTu#H^OOCwbNdE=W$x+yHy>Z2$Au}`a8~#;aCy;s(7P>8>h4^=pQdWgL1*> z6O#wsf$SYT{U%}ykQ_D_w=;l_l^WHO!BU(l zTQLbb)Pcj(k4=B6_&|N#9=}!mgOeAA__jY@M~dc9>Vr5rXtuSt%^KQDUeMC*9w@{ni4a$8^+=QlhuvLi&*Bo7Mkm~ z4kLwp!e{n8hj_E6N1rx)fRQZQoFbj?(C40!dF!_aE-rt}d^dd+N$;de%X}Q5$Gn(r zHkSxw7tYcLHe6;@<%wOMaETXud zK(|yjr&~jZH9_M}k|{TVcAUE=HvbqlvpTz%ex=7*LRDnQk#@QN)#Fhe7#3+eUDz-bK z);kz_rkv8Z^FEGKelq#^TLEPknLTH62>VKpe7-2~jJV~cm-l7yeK^e7vf8$H1{PSK zSlTzV!BU!%j?|nD4$?JWOS*d$`db7q>`3{EB@b(EEh?zt$k{(Pnt5AsHIe1q$&Om+ z>fl(Jy=a2@XCiv*&kqn80;+agY5Il(9wHwSn#3@1pbDq{^G0fn#iyfz&N$Tdc&fFv z4OSj6EGE=s1ND7j=z>uyY&`ue)ceyJXetfTakq(^`_J6;_TKM6di&s%Ufnhru&_Uz z%gO^G#)a=MQ~R)~x5adyj4IMRClmVpA0w61|HMdjJARY6)22E@2MI$32S+AKumiI) z;Cl+TysnA3DOj#_4Ybe= z$5ElrK{if%pvU-V9PZeG728U$mDxYUfQAU?*tjxGv9PY4z8kaIpQl`;1BkF;`GxB` zJq3mw%%^=vY;Y2lwwXnHVc)5nr{2(R;Q*!G;wzOOSIj>~-J4j3Mw@j~4B z-J5+ve)EZI@_rshI9@nDVl0kS+a@7_&}+~ga{0kFw>g-tvt^4vZ-jH!RsR)=nnGU* z+4<3-WE?X)ccAi!ERGAEc@jR=i{*V9PsX#a5uL9-z#j0O& zF+dIFypZZtg|s`r{x#=ag2DAO>bpW-!o<*lD7TehXweQg_}p_1s6Uj#ufEB^*;mfe zHZ|{{BfC=Nf)5eLq{g;O8>B%GL$^T9qa>VLKf!r0qaRl!E#ESpV}gN36V=)X4j5H9 zX|EP>2Sc?=0;7-j5n1?s%dZwJz{ZYW$T^xn{wp64XcA9vg-;z5Sd6D^ZWJ6)js#k6o$KLuxd9*w0w zx5Bj|^PbR@6j*VMKH=3=4b4(vw|lpEBE`}5sJQhBh`K9ID3tC4TEJ-)`t!lq``?@T zqpd>l^?CON2}t%gw-s-sLVd|#=2sKG$O^XT`U{5_9D}Mr%Hpjd`0oa1GH`#qxfOvLK9b@!T_P)kJi*)@hG& zhk@KDVU1z67&hc z%B`6K!&5V{X2yfqx070YkZ+UIztp?BBx=E^ALG}^j6CQ(^)xMddpl6<8`jJjKI8oM zJ9qb6Dhsg9YSK-YqR;X=Dd?K;`C1j+rK5cW@O59>DqQ}iOj^wq2QEyj* zaZ>Q1OJ$)sbc8XQXxO&HX!yz(r#r`hwCbgB{O~^Pnw+>A8!r#Dq$&1kGkut4cG?{q zw2EZfWjBeD<5+QX=8aoyJ`y>fCq`=p;gT9%N6xK07>kjeU8os`iEBccXR|e6Ts2VM z=Zqc>?f()OwKoO}w{cNt7Zq{lj_r{zHxqF>@<`&LmnKlvVWMyIW(jA@F1{Mt>j|W5 zi)OlNybxV{pjhlL9Z-B`eTQvtVN>nD>K!!#&~q>G!Qxr~mMXWt(g@p*<9ff^3MluX zN$dU4q{1s43hyFb*XqE$y#psG-HJHS!M9LTtPFK`Y<>rArH6&7KYwo={Rhj=de<#Y zCZX`ejVp1wjnIC)+IKsZ4_B(lWhXB*;<&wg@XxZXFsyg6+w`a{G*yc0njj4t#Kf6) zW%fdu`)DHP6CRwOBOm%{dlE~U9mPh5hp^%NdAwb41-iEKa2?O^gsDsVYtDg!*xw~K z{qy&GpcdWZewOMC6TW@!|Gw{o359CE{F0AQ%N#r!J);fF0$LUd?}M<+xWWI1SUggc zZ+0wS{Y<34n8Kv$!4GR!e{C7|9)-a&`mLSJTtvEkd;XQS%tB6;2W9Ao3wH3D&w1}I z#y-a#-h#!iuyBX&vT#)vb~3+OaDGjO;NqiBe^fMaq9u&d+`51bqt}kIXTQhEiGd^h zxt1`WrF!7f=W9rI`ogG?bqJ;mYgrQM?!ZuN^UmO%H*m(yEW176IS`*d)31>##mp4J zCJ?&@3;QgMza39RN{jnL-VdzMkK}8`JM&>Xi`au!-uSR)jE9^><$P#CKcof)ge4BF%P4mXC`>pqm=#lyWn5X4j~V09bN9~o=L-rg~#a&PS2rH z#jdg6$PyZ3HijA#uHa;m-+xR$A|X=sv)M~!7Z~;r*mbUX5{VLLKGQzeBjxcwT`}`~ z=yFT;$%yxb$~WhPjdBU-6KnQ)u|f`q{G0!DF5bpyhPQOCt+%l2Z^(zSSQqSGO#5*4 zm<*OAJIZWlYQ^;T2hH#B9YHcH`AEoZIaulb&RDZ|0wx$wE*k8-0`;#l+fCMAV56^+ z%ArkvSaevji8}icrf=QkX+F3MhVu6ri7xF&YO-Wv^6D;Z7I1#UrQwJza+y3V*IA(B z&PNMRUN#utQ}Iss*7wbN9x%->`3ICIyq0f2WkLT)viL)42vjk-Ex3D`Vh>fULArqi z!$mk=Inj%4^?%Et(;s?|7EIPcAvCJqI&IX_gDsrmOl|3cF#PS5*{$g{Y@WCY^FQBX z@z(DP*H1WMcf6XdPon}R@eF4!UOJD|@yi(n$#Xz;t?p7*SHe=8q4LyM6X3k_k*BI^ z4bWs}`LrIqg2L11Lgo)2h0(UaU(W^laF)d_3=Vf9P15GrIet&*lleICi};>M_bWGg z%U?yT<7KygcGLluw0upjE`BA_H$(ULo2)=Kh|Jzm`xNUw^<7rANX41;_}$l8m|&>b z?2vMjBKGZP?#h397JHYHB&GgK2a;pVi?}#D0)6GllF8{aIP&)Mhv`j^Bj-Uw%sU!T z2gGvZ_J6?w!UGOou0uFEwZZ2fT89Ipugz&pf3a_x-+vPMlRdL}^$L2r0Z)YTSe{1_(Jl=>CepO>L;*9V$NJHG` zW*au`oXk6TeuF@FlGEw$H(Mw(O6I&%_<-F5?GSN`akjb;j?Qdln(~fQO zcd6aY-)1l^%sFTDRS7AzMrW$_ts#xoabfb|BN%ddUbQD9sZpYrXsuCPn4fYR3TCjHBV(8+i6i$1&^6hyhgni*ZtBoG* zz^2HLnTprup!Do!+q zE(Mlb`Vtkibzni@@}KHaPw2h9URU?U2`O%QX8-7ZB2D(6|10zV`wk#s%d#9dflg2L zh4rRqRY2R^ufd}*b(j4yGB@E+Htj;a?&kgvkUKm)x(&#eI4h5sNkgoW&D-u>ia>o6 zSzN~Iis?g)YgLl3vEM=Id{Fi-9L&~2kIQL*uNA2&@qs{LxN>dqBOQ@GgKp2ncTJcs z@$Zt$m%+$~ecvv=T*aPM^Y5`)J7It$*5gxI&peto-3j6ZrcW#%$t0_S7ce=!RI! zfpjOB>Zhq6O5uX8V7a?w;TTvgik*tRbB3_hYADIjc?p%Detg!rl@mxd4m}@~s-W@7 zSxH^FVVoYnOgOOrJkAL#^WZq`oq#{u(upq>0S|@pOe9`N+-h;(id>< zm;CXpmO2<3P22DY2!T*#z5fdL@5SL=BY)VsF2i^#(>=R|iv-4tI&Koc1{9YUdwO2C z1I;xea>Lz(Kv%_|_mx2yrs#57vPL;z*n$4aVLKBnP)>K{gutD=S4r*l!XTx=zZPvA+fC_C5Fb<+HH+{Lvi) zhM&+k<_YPwQwEgW3DfHgOU32ABMZ7xT`<5-&(Y9Q3)P3dvm>6o#S{Xg$jQTE{tg7XYk@l3^v>dM1=_BvwOu%ScKf#Z zZE?^EKPtacs&J0mu#r=}8G6o(@hO&n$0j2^=FhEsu+*|dzxCP%472z8mL&YdK4B9p zF@-4TO`2I$KcfUA#4%YFS`bVrtOQ)@(SkZ%F3-}=N?hRLk!Wv7!I?7wFH1h?!@$uQ z&P_CidDVx*U)fLL%8o~O9O@t9$e7EQ0trtT@@A+!%Pj*#m-zw{W2B)w<9jCc)h;64 zoo}UKI)Si4cTt;5osP&HspVrD{Q=28<$66Xc*1bUnN@Y^MjY&Ao-&u&%qKM&HSWAi zxDpnTdq!UYX1Lua&l8HE{Oj+7;lWPWW2I&EPHPDJ#MWcphnpkC^KSVWH3lG2jf55j zoI-DAM9zU-C|*8eH0PcDW8`MFe9^#a1D@MTtG;^ajq} zKFNI3*c(HG&+07feTy`{j`EI++o6!et4@&$hx)0}b?tL~P$;xnE596;}|9!!b0Y=K{&n%nu&=c(EMWk?Vb!Vod2wwqbS^g z*?AK;$LAY>^r$Kv`JOdW@cvc;zYZ!KWleb3Ux06xr4 z9_=W4wCV98>Y^-CCXj3vLfi4u3o?iJw>&yz56$WoN0&M_`Q(D7CbOLzv zafg_$D`6(@`mn`=<6h>SAe_@XdvvjTKfLn)_4`t=B!OA4`0mT3dr0~n)O=p~FV+b) zUZ6()2Tcxatey{dA#K;34q=m@2fMR;;@{@Mq)5x90p=nYIJEJa%6SZD{cp6%DsEtr z|8Pqqw>9>K{7pR_7K&pht)y-l1R?1o>4(Ms2%Oa?vW;KghXaR%*=_~x#zBjpLt(d0 z!RW2)(P3YUprvi4@X^RJP?Q*zy=rEW!nib*>aC7L|86h1KHd$)lR*8r}@Bv8dx?ii#@piAI!1+Q}B-=WBNO@XZkv1m>J$4 zXLOfuGYHw_6_4g1t7$*#~FqlwQ~-SYqFe?H*PY z;aGT!NK_pnA!*}xC!?$~_Wdqu<$OMi|K`QMi#8=fU6o*AAR!PIL_Axh${!&0xy<|R z^8}bqOIcx9W`~9Af2raRC2_3FKDFd2;_Wd8-f@l8I};g?9FgE7A1B+mzDv3Blw z@Wa$iex02TU%!(FBtgdV6FK8hdDCMpc+wPUGh%*YCXJi9du5U-bqHqd)AwGxE`)<2 z;e&2N`>HlW+U@ zAai&`>d41NY;A@31Ef^!FmG7DdTA90!g{C;5(C)1s{}}oxUom)T)o?I8T@%#txa7$ z5=nAlMkzrN#<%Twb){Jp{|2bIx)kmq(Bb{e zi^hX6Q}_0$@Zn=v=2y+ged9D_g|RLOb2`YJ}ZDV>t6CEMhBkCFxg2oS(y?(1%Mrj_I(f(Q&zcLmDU<=l)HR z60uQi?|6NOG&Tmj&*-<6!r7VSi@ZJh*!887|3lU!u6Oq*61S#dAKg2z8<$=|-?Q1% zcRxNS(g#$hn|*7AIhmk|kZ(t@U0eF=yE zO@|#_VuLZ}JyrMZ!w8I3E?1vt-Nm&UT6NBUmtc-ou`>UO}iR1$_@5ibGv2H`sY07#JmR`;+p~XlL=w)MV zM-y+rl<8N$!M$g&>c)aC*K=X4v>Ui9TzUfs9YnVnypY4;6Zbd7damH8FOT7C2{~9c zd+|18?FdY-uXkVS^u*@3MPgTL7ol;>q#*rIT`Up3sIrl84f>y-;MzNP0m>d!a&R}) zL*u~XEy<)!{{A(#wWm@GdS0eVJloAmWN;LDy0u~r2aMR67|xs{GJxB}8S_yf&8o2Y zpT7lzahzrvDWljXrRQmT!G%CWrHidS;Bl1vV>^Mx)cMp!poK~=th6=Y+~T|aw-hFjN|)DA0^kEFcXq@)KF`)f!IGVb3A506J(5{%cjG0CI1e4 zWjq0-zoAbDP|8l>zX z(v_{hyKlV)eF@iJdVT$fUqwC@*vbEZ`Q+nl={q%mw6bPCVWti1_k;NN#53ZA_J~9I zk6A2?(Dm}S=s?=x@}ZXh`+ofOIYR$$8qO)^TgRLig#m?K>q#y>&{S2?8(Y8v;~huO zel&RjQ->uBE8f$vSuKd4PFEb(@3B4Ud0z}OrraFCI`?rj>CV;uPn-Fu`N%ouxg4x5 zc0EXYrw+@aE@3;~k6~ZeB71vZ9uA*B7#q2x7kX&{R)oh6uryY;H($IG3+3knUq5t! z{@f!MWP86s|IEHL2I~PF2|v*3S6BmeCoS|}?|F%_PcB_l&6Kk%ZoL~0%Pdsz1<;Dk|Klh@u61cu6Y6Qszz&3jq=aZ57_DaQ_d_pMe3;^~7_it-M~TfE2T zySE7^Zs z*@8@sP5*YB(ChEq4qaP!6yJKf1(*KId7OD&2bwIdxS2m$!(P)pvTuSNaM51zxv3Zj zxCB__HBy)1`}eUT(OtPPr#0|(+tsZweDHp}bTkVtC^Q@ukrsg&Hrf!~_*ty8NNGP* z_zA1CVs4rZ-^LoBQG?g~l{gnjwEmXf2UDUVz0U%2fqYD@LQrxQS7dg7jn{I83Z{ab zU2JD@oGtQNK|CAI^>3?DT$G1_b3=!{y8q!w(Zyr3I&4rBCz;Ez`ymYWL|l(v+RW|u zA!5_J??KLX^--?Q6yp?SHB7TFo2P+K`V*0T$oOp$w)3s$HFxC@^InbK>hUU5~sfp z%;rVc?$qal!I(!RpKEQfaVqbpP{=ek>OMKoA@=}Q^ z+PAhH55|TjrQbj2KSTe4*Xj;`?FjVzz2S0VKXJ&omHbuk4ED#%*=312gYT_mr(3^@ zF?`nENsNsJ2YHAt1y{E~9bbLuOs*TOWxEWp-pzs;`MIfZOFpdb6R8h4B7}`O4jpF+ zaZo$L@x{d|oyZ_rkUYO|7>1UH_KoW8!}|6VT|yMVSDT)cHK&;arNLsO1S8x z-$gzs5wuu~jEjJ|o@q@BFBc*cNhDf8RvWv9YMGw4ox`l6T{}{WS|H8kyW8XOBq&-d zc*I_Qb<;bav;C0KB`|(Wv5Y=g2P?NqWItWc%z&Am8* zGdXV)V|VZXt#9ksQ)CxJv*q(UL(~X#P23&R=9hrF7|pS`5{gs_rBuqq5hO95J^I06 z7RV>R&Gjp-LCH2b%R#ymNLIP8{^)lak*;X=g6PIJoV7Q3<$OL7Ivq~1>yG9^>-q6O zqft{NH??lN%fk;1pRaaGd8=TF2TRjYyHOatz~aGj>kIzg`y1cOUBNYps_nnDP7>(o zJo%eR(L}nWES_DGirA}tD|%Vh9mo}PZ%d`JaEjr0-k#iApeCpq2hwv8w`lk3eZ8Cx zv@w}46;r&>^?uuGFS|VC(=tzl$vfhl?VLo;?_pff93<#!G{LgiU*&6wDLAv2^~P|I z8`7>D%^4iug7eN^+dNDk!O|5TF|uYLwtfCQdu#j!%oxkK>%U~fD5IoQ7FG!Yqn*Pw zj^e%O=lq&Ktn36NWrVX%-C%|;I)i~|{bWexE*ri0>?@A^4r8SC`{3A7$voM)DqJF; z1F!KB7-;VKGrY+YHH+5?XTH3|g(P`1nXiBi(Uz@(hR)EJsjzB2Y=U)jJ05bJ9l@qv z<6p0}XK=aUiMqsC5tfz}+YW7fL7ZOc-g)ILe7qfQcUDOmrd}+CbKU+AMw8xHIPLyU zpzC0|{`~1VsB7T$ddzpKWC(<0JE;ha>y?PpayO!vToa&J((iXdwYY)zn zC$F?ot8vVKwSAx?9cB+LEE*4=g|7PJ-g?*X;*5vxrA)_q=>JOjeL@JZOQR$pzR(%! zwNzfj9@&oxZd8#g{?DNJNxu;390iMJ?#cYS!-or8PbI1%`eEdZcC60w0+G&VXOT|* zcWm#`XnI+91lv4oUtYJeBhtyKQB0}UIOr0)?ES$7$4g9_c)1c`kyYCyeCRYZUX~uH zo$rLfSLGgl-_>xY^@m+>nIg8yw=o?gkHN;3>Pug5MB~a`R=Yu-O;4Dg*ME6u6tin0 z5A_D_z{QfsaY`+L(4%_xPu@2HAl0_jkywr_Isu$U$udWiJe&@cS2!QVcT)GJmpQ#PVv#x zyNTnL+imPCiO^@)dzdA51{Pk;5^b$IU@|(;FKrmGiE&4BT;pfRUZU;F@&AT>^byT29WzgFyOn@ReofMWjS9%jrQIP8IIl zFShFgj=oO~%b=)jI$rJ9BWo(K;aeM_B_n~|lY8zNdQ3vy&#e=zcP2LX2w~XRBM;`& zce0NaB;ugi>q9NRD?r=1>(}@E>rnooA+_j@1=Mr@JTB|#3A3SPjGt;=LI1|L`N_9Z z*q&(4mnuV#wWO+5eGUSQ8ArEX+LMjc^AZhi1?jk8BgZkM{0@fo{@N-`8{#U{ykaGz zKBmX%+^f-Jgyy}|q>2sM) zM;ef_MWhi($`0IYyOQadb_WRBB!tlqDhv8Ij92S^45qV;h>ket| zDt`40DF2O!`X>=!i1oFT$d+xJd{L)-u^cvY&4>R0K0v|0hhQL#-JAT{NQy%$+uuXK z+4Eo{vdKU_NC}2jqs&74cL7niEv9TY3EQ{I1wVQ50LVEj!RjS%V7A`8cW`tAhbqZs z@h|Qp*}z3Dg-y2gU?6#02-v7l0zn*&j(X#pvvc$I533svrcFhuGS8$fdySS9T(0r4h%KoUY`!7!ep9M1;%v zBVpj(^TLv{Q{O4}mS%YSd-uDyFZ{rBp_IEE$IH6Vh!mv|)2YMVoADdXai1Du0 zVQmTh*f%A{8XY@=q~8i<3rRCL(k;84EO{F!CzexAHZ@?8-qYOr(>tIl;Irp4L(C=z zLdQ0*laQh8b?$zBJ4`D7A;)!Wg;u*0EUU|r(9hxbNvwH0P@Ue!`pBoig3-gKby_7h z2j6lvU8#qTSEl=ZW>IlrXVAi*Gd4J4GIGL~l@-D=vJ8vg_XCOHT41D1HQEP`j~)rQ zfD`;Y`~?veSU0nHA|&V|bZg$lu{cp&kWKlw+{Ojuj!VRfBNvf;!VIks+(+8;=G0O3 z9&F`VH{<&~f~_ain?_%7;PRy8&d@!=Fj>;Cb1`{8{G4oxKHM6GRAqj#wWURzx%ski zeFqcvFSYIz7^}rG)vc?)9vi^mB1eyazZ}fddzxmr9*3oCU%mEz<-*yiS5k}xR?w)B zwRWh(6*{hFIC-g+X$wa7h#GByI^%lc3=y+oiK z|0eOo`8kkQMQ%y-|A6K^?OdfZy-05Mnaxz?fasR6h>-cikaxzjfcrWPR^@B64quIc z$}etq@@G|`Q#;D}93wY$iHSe2J0t~R$Np|FdZ3NvYA1S~;;gabI+`l~xepU+as&%g zUmUVxyJA$<2GbKRKR@`Z;#z09=ia1`Fg|bcSKEA%$apGjQGs*@7MZkV*{`oaTe;CM zo?VZjqqaiw$z4CJrB$&GZ*nY+(YtT%13OTkx(N6G+m4|f!&i?C_TrT7(2ZXYt^vul z_=)40I3%-d?9u(S$yc5a@?Y=ehS?|WOe~7;kYc;&J!Y1Ms{)m(+DF#^_nzLqt?@+D z6IaZUGl<{&O$BD!)0*$7p4_Z&1k3kVXK~1*KE*qp;~zDafdg~cFfayRwj_h zr1oR{KWEG?d3ZdlxC53a9Ax9L8Zs+0%HByD!Q%I}jY1L`=NqaVo!Hs&kHa7T|9g+4 ziWVlk=4Atl`!g+{h-h5Cw3QlU9tW|_j}F~r;=uMVc3G-dxuAf>*zZ)!9-wJ&h*#MX zvHTTJw3LB4%-+>5gl>OW8y=f77?;MMTzl@ezdr~=wrv`G4uY^0ANSO%eGw){P4)LA zY$q}j1gE~w?ZH8%ErcJuL(tckdt3U)aTv^P?EC3D0bQ5OLwDs00C}{~WGN{JRsynu z<2QMx<=7zos>BA&$KG@9KYA9ZH7i5sSqAafz_-?u9-WYIdgOe?*Cre~M5;dBR)*_t zzdQbN6u_M2InAlJO|UfeqWNg08g>W?4>n-+QP85!S zg$IlNgZ~x6h;Xsg&M+2eHVfX*7Wo1jS&i>a6OH=VDOQjPJeCANm6J<%OP&z!I06blXT&Ann+1_9sMY)3><9 zVg-0%;_uAjT2u^K!1 z^Iul^i$g~3DuyMVxHv=k8?1a3M|db^;Z2i7y3n}?aZ6vYomdnpH0FTwr522I{(+d! znZGnvrHzZvwohHwPJqElqb0o+Rb2Y(?sD+pejJszD(8skfWoaS#v$4j&`ED1U*u;9 zLzKMsW`-{?&nUHER2YC2U7vsYI$XoQ-6CU}^$b`<_A@E;JPPC^4+aO7ALE$9r=>TC zWr%biOqZ|N9)$5B3s*zir_h$LE82+87zdc{?w?bhg=HbLRhqRp4A&T(*igSj+!9xR zdpLmu>U>;!DeGHs^2%ky=eGl4nq}LsY&#zqza{&b1s=oj_0L-@|EqxJ%cKc>tcCp) z8>TacTuAfEC2KP7LE0Oy0Ief2K;d}&(r&y52Qs!5tnaPFrLdicXUsW(%-MF;o3t6A zZAw>ON?Jo19{xe{4~KasGs%3FEG)V=`zdJg0gQb<_)-3V2hiw|S^7^ZVEm5*)-DYX% zh@>OC+Qr;0abA;;Mek%SPRDj~)vA=?ypYRz=Y3N+6unnk$2tgS=Qk>)UEN^3gxdK> zkqXtt48i`pxp6ozL4ni809FlGB2z>k!uT9l+#6E@u32wDKL)~J(6aLvHC~vl;#UyW z{{f{!m9K^Pbz!dS6|Le}0^-d%v|MrnJUGl2n)?JRZHk*;T%yCW%MT=MzBR&%sjf&r z-ys;@l6-<|;Uab$S9Q5DKgAw{rq%|M9Z=salP&LFfU!3{#-5Kq;do8a{d4!*aOBGX z^>YRZR(*N3c^=R~LtfiR^tWOd^r16#c{vG%+bD{xv#(+3Co(?TxQt6%?c2f$AEEE= zjq_#K>#&iZU-{0s2Z2G0BO`Jf1v1ngJhb}R=7j~&*iGc3S7Twxn1j4C zL>m_MeSbvkkAzXRHpdyi4kC+ka+ku#1?V6f+j;O0U{`iLDWkj#C*kfQuC zPWs$$B6Bx~AlJoyT>QP9OiTl4Qlk8W$OKk2 zor)|y{0HWg+s`i<*aB_0=h^hF+?!r^ZY1tS5o~-2FO7+{C(?7M`xS_N#L6$+TlP6x zBKfbmJ@0)koYuBUn0@Yt#J4RkHulpKnO*I*)%!}&uZP+7=e0qI`lUiXb!+op$dyzn)m|NhS_Hx+gc7oUrwqOwU$GF;=qtk%MB!3I2s>sZGj%|jq0xZI?x;YNciNr zNtn28xu-c+7@CqjZgk2U66mbetn8GTFyCPB=|z)fSZ{1gKC!x)*UKh_(>hhyBC~yO zV7d{`ZK;XM6&%Ed{!iYDIU_jvMSH_#cPcbRz8Dg08o^508T6xi;E*n^NMEIcS*2V4 zU*=1ou}!DJ{_u0G5Ack9@>T|C^MGzAdka>nr`;0yorTTJDXDKBWZ@Y5^9CZTB#!gY z*Zl0S#nInkpVx!vu~nh$mA%+0?Dzic*i@7aQI)wbIRe?yLztiIVfAUCBq}hR+42iU zO%4SbxZc3VhuIyRA@0!h%~bQn?FbyH7u^%ycN<%dTL#8+{h!Wr{aPZF0TJ7}3zmtd%-`K%V%9{T-yd0GW`;N#ub7xJUl zku)yH!aG`x^VY|&w%_%D2D2Os%M~ij<-VJlcXZ#}N9P1pzf8bz@Pp?gcLlJ0URQ)w z=o=J&qMTZ7+zO*b8n5^7{{a&T8{$6izv944LeiTnen)kSusJ}i^-+uM@l6n)yD{Orxw%*?4 z(^&niL;3`U$DuU%7Xp>$?~3PF=&_er&z%tGhr=QIe`xzWuvHKCJ(<(cx zJp0|S(%^(;?)n^;4aczNI76~jT@?_oMfZ?Te}@&mo43NhxnY4d=Kvu-5jGsHT4^fz z_=JV+#r5oU?8;gk7w`;*adHa%{r9%u$QSUeEW#1Utb2OiI*KFdwFCQsduO0-i|7M- z@+kK7E|9{Gf5QRk*TH6i2(+soFQ{!aBW1Lstecq*+hXsOTFRSaBbc5P{fpT0Co%3r zn=SryAqEVk|G@=1(F_M?B`g|Q%DCKd0cT&Pm@HSv!1zJBac?JPq|9l+*>|Qe(e+kn zN6K658a*u@xGw}FpCzt# zxM9C*+US)y6z=j!9=z@dEjwQFHu=|KXP1@CeYRC>7W$IK>v8}o+a$~NP6y%S(^rPy zB?h2rLn(bv^KtCh_Fpe8U>ZwuITEiio`x~r>=VvOKVZ4j(?n%@4jW8-kC*OIg0>X> zegfy}CfB$`iTZ@&yo!}8m&+>xgGjcy;CcwwOSr~YUev)Ky@TK7Un)X=l*$7u&7II~ zd|WTGToQVJG1pmq6@yN?Ed#30v|uo2#rQq-|Ni&OVfQk#F-)+?3SMQu50f_^r}qUN zLW&M$H~Tw5=n^^n(UV z7tdm2Ui*I+9&tmT=@mBhOBZ2aGURbx^CV8c&CQ@4F2P|2{zu9^uVMDK`7cH0GVD6b z?)zesTZe4T+_)psVQ7zvhDDz)fljo0d(sy@=ukFtBDzUK+l3!#n$N9~@@t##*9-S> z_Wj*#t&WdS-$J%o4OoSLE?TxjCBLzmu$$jv;4r5Ad7qG$WQ?SPNs9~8tw7r4@Y5p$ zku2+1wD;01kSDgxX>3U1*x*JhfohEN#)77_2!G7-_@fg%GXYKB3C?DX^0-DFZtG&x zfwc?IJtbLB0;UDUS`Y zSKZ-<)2>2nNJ#D}XfVbZ+v0;_&r4wHC{N<<7x6Ivw3U6?PaY@^w-N4G=I?-mn(3FSMUuN zw=I564zo|cBt_gB>&$*XSe8Klqr2r$j{}gH&be}NXJNZ~@IJS}OW5@FD&t}55lC6- z`@A=}7CKb#-&XLf#J-moY~nh3ajA{oh3EZC7^Z*H_)qdPQtilOo&g?QXjsZVXt0^* zM{TQn|IXnu=MQQ9v2Qqi?fNUy>B0Yb`fP8$cD%1Kgmf*An$}ssNK(Xg&R=dYozidh z?#UF?#q%lhwY9@ac$StemjZU2e9^3AVvW@D1;flAGe{BIcbfD{92@o5p_mhQgVurJF-D#Nuq&(_$H*4Vo@-lQ7(|Ln~G9t|Y@(5>3pEQX!2 z9fRdwoBZ7%Oy<0O9_w!S?my&SfD7@GkaOlH_U^enH2&HNCl-6Nyw7sN%#GBSgWnY4 zpRC2aDn%8jjJrF@=`OflK6PfM{Rk#7t58`FIpch>P^WKv2Z5g5mwWY)1a6$Yr!EwTN@jgImQP@$nMH{wEU87#yPxfVkrvl-)yGrV zw|qERgkcX*!OpNtu>gi&K9&)@TL-NnO({qIJ|WP#&3vO~cfnGI)I`paMd+0dl+@?7 zhcWFiBisKj!steR+npi;(2}>+A9|pG{k^vv0xor7x?`2B{ktwCnpWoCZZC$(pR=|G z?(eWuD=RVkyf%TUpL9(;@lC&EXd)cGav;otXR?sAl<`wltiS0cYk4AaGWS2I_u9xq()DqvSs`Vafx zG5CKhIDox6&6Hc2+&JMo`N5fy8KVMH6n_^IfoN~I^~t3|B8%jY-w)&*VDRIux^kC$ zIHdoplgq3CCNh7gd_hqhXR+DxNqPl3y3~*G=GI`b-0Av`aS6!J7v6t?{XJ05*giA4 zlLDQ2mxkkaHo(NL=W_R(=b%Dbd8^P_Hf%1m7=JdCfJ+AS42DhRu=sOsJm%gl?6deL zPh7CZqH}w~Tx)A!MDFr<|1me1Q;4u!a`hu_d8u2;ap5|U%@b9rR0kaVz9?~nF$+c& zBGteCJPPApN6q$0=ECHx{#7HzeL!Y+Y&gQn3PS?s>`MJc(3Z&6A^7loI~v`2{llZKhFLDXaW3FA!Nnt6=cK$;f!moAAQ(8W}<&3g^u+I&eKMXVP8_BLK` z5bA-xKB>wtm&K8GEbjH2+Z&iFHJbcNzzv71{r^n<3qi`J%-IIbfzK6%?J~9`TpHD{ zeV_CL$uPvC-bTj$b6Y+M5+C8h;6J4n+8bDcTOpp=J0YE6hgwjfA{Om)-Uw#h1&ddx z+Zr_sfR?d$3-y37k^Vz9r#H`OXyLp2)2SpI`cCe>;KsL#{U7$HN4~3tb*2uEL1kr} zI{k%C^`9Ax4X@noIVZ3=54Rll9Jq-i6Gh>#+?Qc_`ItLT@eKlf-jLNT<#p(Laf~mj zbPLv+^q%Chev8$+o%R0pGr_3GufF2jj2Ln7-7(tZLR<_=wu+Cw3d=9=Cg)u(Lvm7= z28)CP)Oj@huutFgC>{N&M^jIbRC>>4*Rn0nWwAb8{O5%o8wI-^1vk%oA%oqkz795& z|D3)npNqN6G)1ihW^CaXJ2~^R3Q1w5-!s|2!}8vlpjgjtT>kRCEQb9*__KKEYlqZf z=zl90oH59WZHl49B9~|w{^lWNchU#OMe=uqg!#Z2r(%$RHxVX==W?!-ma*#IpFD}7 zgGkP_J4HL;j`bGN_mWQa;$+Waq!>Jg@i1NHE9NoK|3#mC_FV!@`f$%Np5%v7x2Hgo zKZeVfm%9~-m!Pqr-XnK>89PQhwm!|ifZ3UkgIa#m;ot9iB{H^}#4U-!@y>02M7ph( zcRN0QgvrL^0LSBTusZfIy|q1*owKyf6v)Mq#Wx?g2`n)1$E1V(Km-iQb)VGR?Eov4 zx;^E)F5-HfIlH5JFl_Mh>`>fig7ZfS86w)EoB1Q*e3$GEE-%#e-==QEsz>@^bm6;U z=&B+g_2@Gqosk|*NyiJT_qfhJXH?U^k&=lsaNvjd`x>wSCIv*vX#o; zIeZgo9v58cOzffWy|~KB(t|i(@A|K^<{@@m9DdUC&H^bED$CV-{9tp(Sn6i&7%te( z$OU+M5Se-tc&uK@;rd#`$-AB^kUM4>`?*LHquI~gdY6t!sbxs>>N^J+j*MSTw;4fV z+SAmHelH*^JMwU(Z+b*SgD=tF9(wQPu;g@Fz~nnfsGHga`J0BPqHr2UOFae^f|_7b z>1vv_!4+t=KYgXNsBW_#e5`5D7_i~>gPa-F1#G|jT~br~Hm-gpxhUuT#WoYRQHf7& zFhTsUR%lTTr#|{`({j5BEqT_LQ?^%Pm)YTv=AQyErLCH_u{{P`p;Y1PoH};rP{eNE z9KaQuMmBCaAslH@dvIFb1x9a2Icx0k#Qaj9r@{_}*jSgHtXb=g)ZapFf75I?^IG2e zpnx0{qHSmQ&SGe9xX+cQD2wCi@1B@jAA^}&K6`%NmBS)#iG58a+R#rHKX3lx08aZB zWUx;iz=m!~%Ly@S?3dH()bm(?sPy82LL)(>{hKqI=3RuYg5J!Ylj6`v_^dhh_8yi# zy{b@ZeHzI{I*vj2&tl6J{qf_WKd{oq1a~bdl z5~rwj9B(D$?`Bm@F1NvqGj+d+S02FR!NkD56k+K8IxK!WSqnDOs(-yZnTxcinB{H2 zfP>b-ac(^sP*JsirpWXh*5;>=g-#6ONPP!aINNuiycsKOn+Ke5tj?phq~ogAOF^Y6 zj!phPw)JT|4My`v9D+X`hjAgZo^&ZIpd_9)Qj%rE(NmJg_Z^lc(w)0{tD37AI(L;e zd=pxPxnaKL1MZ#BUOgGn^nDY#gf;C(wk`JkbO}w@^?eMLY)Qs8;W$@EEBQ&iD2BD=7%^aqHs9-oDM{DUftumlNG4V zubo1wXh;okR^t(OfZ)m$#%Zen%$86MzRy$*%~v8AQ(ePg=Dy754(U)VuMn`j8F>S{ z0-b58MM2njA>g<}awb-n{FfYNDS=~fNw-|49Dn3iiEOo>hk~LrQU_&NaoFd__sYR( zobhif;rhN0=C;0~F-tRJ8?CABrQJ9TPt!(DPE)a1SNomK-%{w5J-p?KTpP5##xb#7 zX*gf?%76XEJy^J8eRL{NkU(#G`=IW>W~4fkqoTrJ<0nF-Qqa#HByAjQ&2_$mxje;p zchf%M0NaD6>HJ+-Jnt9jKcR%-Cw|&{|E3Y>Psh~9$423t9^)Untv`_ZS7@*C%}4Mz zCx@@A_a`ooRUX{^N(ClTiIv%$n>l=EU+R$X8qQo;iI$I0h4$s8H=YZpq3w5jjh8?^ zmOCVJZtOILE>6B?Yp1)g^x$72>t7ul+$gR49;yQ~;idY)^D&!q)t8?>EQF-r(yHU< zcf(Th>d@;98z@{C)Z}~Xf@2z@@0G9R!XlgY89Sz3u;y~d?rshTv>O);6%F+f7?`n; z^vwuz$S0ESe__F~L%WupFK9zq{fnVTZ}eeyoYoj!;DZa>=(D%=SHom(jkNf#t1yw4 zD5^j>kFj+Bw1;X-p)79ribaMe4(!?56C?TqCsgl6DVNq`)j6GxLySpK9;_~Wr~Arg z9kw$XE@l9#o(is@?!d7>2Z|Jq(`|BELrlz{LF|zJgfI8DLE-bf+Yhh5$G&YE9oI?= z2#oDQCVLY`pilL{>x7PQXb&o%4_(v1>brwu4!lQj*tSFY@QE$(_j6qesXq|Orjeh6 zw55>p!XV+I>Pr|5xKk-U`VDK<_EN!+!=!IV`$ONm*vsA@oA~)Gv`f=j zxfXR`WY_o0MLs*RnV(#oU^}i-OgQ~s?ZkTTM835DZee>vir|sm zwb&62|Qogmuii_W@?yLN!0rh0V^Lwc$Fn{tntzqsjtU88_RQdT~*K(gw z@3U_FI{)9D4|$73x|xL#Ek6mMp?{<4F-{z_jL5wj8;R}RNv#1b>e%tQ@%*S~5l;WA zcKvty4>W|1`f>QVVZFa}aFL!p4DqRaje9c>vzlYC+AIi2l^WIYpps#%F4!rVw-*+j zS)=x}-+)owrB^SKW^kIKdT-Wc5*NOxFRQTI!(`}p<@CN-m=n#up{PQ{#Oo~D$s27r z$aDCi;F=)}IR5P0PPBx9Y`K436VAw ztlD7wRJVR*;ykRHH22=f41(^vBLh9IS7CmizEBF00M+&JE2MuGFt*IXX788+q-&vh zJ{>+-?ZXgQb3A*~Lrzc%wCgY@k+-+&TLVN&owafMFB-eQzIS8J;v>*g=dZSX3xyHU z$L>mkrZ}C#*XKl%JXtCKVx2mFvRe8ReGs`=$hA%2ZB>5(6 zT>UN?G8+heDs@qg_GLJIG}BV1{}eP-ux-=STY{vL+aXOizQJ6})NS8wGSC#Xjr*UC zHq2Vn8T$D{X5f?}`4{ktlw0tiu!M=+Di#&?Vx;F4)Mr zaSRrJ)u#u?RKfVgkQ)m>zu|h5^2P}~j%^kVbc+>QNTCy5{rU3|OpE3Eg_)H>vFifQ zMC?PLN&U!{@qC6;t9^V^v}Pd1Z23ID`#7dpURy}HQUNsT&yTD9EijiQ{vxlt5oiEC zg65|&B49@LqNOs<-v63zA@~5uB}-2ZK1LidD=@wp^%bhazjm*3wgJU(NqXtz0r(X; z-;iYf24|m&U4JmVS&xfv1^3+Vz`{{>vK{AsAgOxh8E7|PO)CpG&si#Vdr3dG&k4cF zcgL!Q|Hwj1;d7B78(AFe8Z*M^<3vX9POZH@S71<0cy#bp5KdXY`kI!T3QgL__NrT7 z0ZKEkNCYg7cso(Fef=0IRelFS!uKm1MEKDru6j9&W$X2gN=N#P=m zDuRt}cG_i9k+5RsA;%`tK%|SW9I4fn#uKDX#+8w00E_2FOa=4kk;18@ z8PS~vYbB~Iu4|_-A<2IDD_1!V3g4!&c!tBwkAx?0H|F4liDKrP#ZaW)2R(uz>V3|9((0hA6{7ol+@K*nd^K>WQsB*4DY+-+; zw|fx`UXHAXE@Ahuq7N$L5)O@8~^k?@EGhs#pjw^b8yQt18GndlHG zEx9;VlBse`WayZO z4Yqg^hMxuxVyj!GhHu|%oVC57s^->>OO0zaMW?Pq{XYI4x#5$W{1DJFcdi}B-UgcQ z(kj3~@1UcKo*$u4(V(Lx%mOIlI<~W4%3-o9L9ACr9|yL*d3wO?G=cHqwlGZ+KMsf< zzWeF5943}0d4Ekf52V-={O6-EQ))JsNbMT1#EyNaSuFOlxG#&l_ zI0V4-UEAdQO58ZcDlK~|hv>2UQMjO|ruW^5$m;Cq$py1X7RlXXTcQ`r^=zB&2qs+)(d z6Y?z&=Dx#-d(Ju8jzOsJjVWH4+zSf?FVoPzJRlw|HTI=6;;86Y&RGRX7!2-Xkhoij z8;K`s-F=*~W8sXk^Rp;ie7kF2sX!FUgIMGOR)wLpxGjk2`3EO8e%_SJe2R5#O1CGN zL!sV;pm1~I3Q!%thBGvfApKPKC-p&3_;dT>{`j@CxaiX+b}H!z(q6PK1g}$o(s)=Q z#-bcr1(`>j7=L2heMReV&aK#$eJZMbSO62>J-ECdnxOPLgZg{BXe?U`mOt>H8%~_N zw$0~)C3f*IJ$;*MjB7tN-lX3AjEg59-9O!?j%$s$@z<#ZCLDO*QtR6{J)-}uuRt2C z6ITTc1D@l8%}{dqUu{U74Gw4W2!s|p&i`cee?jlwmSg!-E;ut2!*$qx6q2mslGpr) zaM~xXq>@(+DKifWE^_%nY;>MpL*9NYWHku;C$R@6t2>H+?l*@B!}23jnIo_EV9}7YI!1Co-7Ai?R20?=Y`D390NG{rn<_AyFeV`N$1!?E0;9 zS$_8w?63RHB&lZy9V0RBGx~v8{oDHe#77fcsOcw%AN&vJ-vw7tWM9FIlF;#vM}E*c z*s6BpWhzjs9u>SdT*LK1Q_YLZ2Z&p))Cu?17!c{@s@o2Gh{AZM$^Cwx!;s;X(5yl1 z#V)4DCl&%uVP{WQCv}wzh9yb<2hOQr*~Q@q`|lh;k=`%ioScZAcL*!SjUm`Kx$tgq z{t!$pTEs_Nauc=~D0q*V_T$K;+e?;Lz2mz$5*Lv(2rShg>*t{OgUw zzwL%$L8+wrzrjzq5X~IEKsj%d*2T_1LI5Qu1lo&0VXahyDaEp$=3^Jnk{`@Uf;m9XgNx$W~RW!qIgnZNJ7N=&+N%!;{v%+eH_g2pu;+NQ$beiwd*fFDrwrXY{^rT`{urpsi% zOT^cLqr)mR9_UCe(u7}W&~6z1JUHh)7Jg{w6q5gel;*D)adNi!xh{fmQ+YcMpVgb= z_*j4>mOf)9Rb!|Vk}{6d)P;Uic{dJoe(Z^r5D;Fo!1?(@ZfnIzy~-zFPPDBwVscX18R!2Mw?5iiN%=VN`uG?^sL|l17qxSyoD* zjaN{QNAW69MD#34-BLvQI@3=5AA6AMUw^$Ichl=qw025fPsfZXU8BdHInWs0q$w<8 z19f}nE{C1h!YTb%#O)@^IB@(I_buTjZ2miZ*yvg_b~&#&RUa^isfALZ?WZ>D8-7o( zhiMgOQ{1KQdhUc(=X;-8mJMKWSWZz@rxwZ%Wcp+flL=esrE*fQ1wnh!Rx_t%E}R}} zC%#Y=!9dN)`jPBLm}^z$dcf9$B&{nMy!8)Zf&YuN^-){sZyI|Tn9hZLnv&8_1M6{^ z-|K$x*)f=wtr=G=aOg@?Bim~}z|ibp@fFcX6xYu^CYas#c~ zuNa#z?pnKcULQMC9=$F7U=R6XoNczR<&iS^Yw%>tA*5_wI`FjeGL+FBG0puNjDs~h z-{f5L__I_j{(DgqYIP}I8`&cCWTz~f>O-1AtOZz(GGnRA(FgtN{(t+k@^8b^`hJ z{)gsn@i_HzD$P2$1!|eTh5)M{ysi8C!ENa^mfJlk@c!6>%@pIG0sToZ_JUI}_WU1M ziWjkH-*vkezXb+t6!YmUr(o4fXnfD3C(u@?)~9|e6UQB!?;oD- zgxhk{~@IXa8vBXUSw(IQu4te)u9R#_kOuRaRmC@)W1R z)|1#YYgo>(PR3G|`t*X3Conv1Ba+~ej*I!rqKr;{n;eoIhH~$4DY&dkXU|?NtLdVi zvx&w6xw*=qBxdXl;JssI7y^U0`FW;ZWJ7)UuuATm!2ka}{ldPzIZrfz`uBH^;hqJU zWGVg`xmpQRXL;=2U($#1xcb;S`5Ulg^ysntpNBA6vNzdI#vhvm+NW6~Td>XT%6`=k zR?yXc{@)PyCzuM|zGb-j3N-WDjf~{wSQA>kk1OyV_Iqa^@Gl93K^9|0 zHKE5)vNce!Of3Qz`nHv%Qhcz^TyNo*q65bBKhRAUdVw|CW$#{34Z%>;>WTR;!MNl% z`$M$y9aOY7E}h`j#&X7JcD*cFpp9=yewMr&h7V;sI=5>=TQ}!_Cd2Z$?kqnR*J((g z+h6&0k!cV~AzfSP)H0A9YqY+w>0u8b zkCQFDxWzmfhV&#S7`Uady;%Cx_dB65RP)vGV^I~5YcFOWN!|*iF>d*9op#Xf!sPC| zgN9Q=Gbz8Vn6S8Ij6Im`6U>@Rwh=x)gPwxR)IUPvNO{VAvaS9x6o{Cp80$ZRx-s$h zuiWckF>g|!WIzg=*4{}P5zCNr-TdDt?1$kUl>D(v7jgB-#ACUdDkMLk%^c=!B+&V& z&jgJzLSwbD!mWdwocm(ezR-)ip=IAI#}K7NY&r9J#z{SYv%f+IT6zD%d_FNt)l>wl zV@^5kS76z!uUF3R8>Z0bEn;iHdjrOkO6#^IvjK_HwfpF-CeluG^8Piw2*Z;%*q8VN zaEcT+mjcX??9WUh8cspQlHN-c25pzTK75zEe&z$y~gSG)nFL-sCvqnwE{>6Ei7N7e34q~ zdm+)`H`F9N70r1N0TV_tir+rpgi8OL71JV<&>PWqc7E+DtZNwVpvHF~HDW`0={_$^ z+ld>>Ee_l=e3}ogq)mi%ENS$yu_^llX`FcvEA_lHa z&QHFk;Iql|6gk`8n=gUFB!_>dzX9p34O^+&cFf*Qqfj+CkaAMSf2@of3f7B)*q1*; z-^wBJ_X~2^k-byvV)Y7;tgM611z*CA(g^0)pYpN3RPXkDkt@)y)oJ})`-<%W>z9+Z zB9Z+&Z&>lud{}C}ns(VU4tnW>))zi;K+@8I`nbGJKc%ZSxGt}PQ}>)HM>E@?j#hK* zXR#3s4~gJN=# z8S|x=V#`3`BC8!&2p$ba5d5=N7{Sa_0Hz$g0&BdC;NA7|&f6YoA~@>$^eSo#RZ{~ONj zej@|D!*r2~?@l8*McP83umYOx{yV<4vKwY11P+K>euVzRIm$QAZ+eNGN3x*80|NcO z-MpMpWiTwl9WWhai`}8mKX^zTMB3t=sDf>O2@GL=Q$xq@BS|QiuXdL@j&@F1uiS8f zX2a(azYZ5btBo?%Jkc4GIMg0=n|Q;t#a-^p^$I}xBRDVW9fXuM!GTrTWJvJ4@R@Z~ z8>)lrn)Ndxpz|`IH;nf`Mi&K}*>7EAebx0Zz+|K7X0Sj+U?!|ewL*reGn`EO1n0j4z zrTFq1<}pxXcJ7RWrIDHf6~_daVSTZoR#l4|EbYMx<3U)>^8Ub)@fd7aq9xG%J`cT< ziFL_u$*^EnYS-VVhyy!5=eN2MVLboALl&wJE(xA%xsqKEV=LSj39p&3m5|kv{XLP$ zVDWIAL*h8J{Texc#AMSu%|A7Z6|&>RUf=y%q9O#k^J8aE`7Oi5r-3fcJ$$&~+9;W#H@r$2P+D*ih*CPBv>hl6$JTIqSLu$6kg^nCG8&t@r(Y*ZTc&_P*9$XRT*F&%U4gx;|^|wKs2Q%10pu{1ZN&8+*Y8 z+HQt<2Y5BWWF{IoS^vVpV^hCB9#Vm!^r!4=vx&HpO7B13CIVRjKdb2W8z6C=Q|$gt z$pX5kLp&oZ2H0X2j8E>Q8E#6HG$V}B*OIx?wDWo z&CuhMBNn}sVtJfG!inb(B1_z+u%){G&EviQ3K-8!H;)=kq>F|b9FcG!-)%AF)#3ox!h4JnfEXl*yDBE z;|%u8Fh={=ka67Y^oUveWgHx=VyhxNz?mq))Ye!RESuZj$W0u88g9Gq%+kU*$n$(7 zwEH0rH0W!l=Pbi8i$y@bhZii68`TqvTcFc-uG&oU8mxcK1jKi&<)D;77-kHiC+{kX(^Enh51oZwXayo&WaZ#8tQm~rNr zMS!880g%m=BIQ|n3AC5h@cTAJAa~trK3F_ozz`&N|9i_u0nM#ZNtdv@IJQfg|7!jJ zeE)1#>Q?Y>s8JPGNNxWOf3JOyAzqM!>C3gr?sB;>c%=PI%bYdzIha_7L``FPiIex1 znH3z_-p$=GJdTv%<4Pk%GcZC_Re0ygOQ8ETKPMyaO@ zi8+*4UypBt{=>mx&xK#(;0*=k!58tyR22Wuo^{aRq zw+~LG-+3wej~=A$Aj_M=yV5!FHO)Z6{VQ!SdfX zFDq?@aAp$l|cy=GwyB@31QY^&*mE%U2v?g$M z{%%UY=^GrhQW|#Z=ftHS>=XGR%1|%l{XDYfDlF*>eXeb5#Bm$axvrqkP;MW>O8<2; zpWZooba&ie>~qR~tnue8Oe9+gvM**pmB=0=BR)S-;Pbcq;6pGbcM$yZ_JeCA}O4?pN^__;SF^kHAg7Kju0jH7haQ+z)qajbG5 z+jzVucBt|TYqoENAU(xXgNo3YiCZjZB}Hdd>becCE4jV%{HJsV`)g@X~hN_M_W z$I%N{El;qWfq^sHjbE$uVeWmSfxKuKjuEZSafkQ7r1l*zi{`gL6>-mJVg8@zsgCps zoN|I$m3`-qn{U=bNx1t6N ze!hB*>2fW8z8e+j_vZzS$x>)GY%QRArid4^Hs>w(F?rzeCm7}0clFc1&oF)Qv5DyL z1Q7L4s8hu6VSlBz=FSgUSaw>O{-RJT{^4~#!|QSlI>MtVbXDR&*)}e?!MwPc<6R96 zGkORN%+d|J2e-m7J-?z~wLjR?l`L$fL?g-2z3ce(E(kRWbJNsNf|zq3_G~BrMDj2D z28OGaSjo`Nosdz7-NGE@8x7GgGMzcgqyGRtcNq!m?!AYj4nh2-tLm`ee?sE0>3+xz zDN22Mqz8VM_V>nWSHq0O#m9-wkAcEt&6B-P7k}SsHjCJJ0n@eBobgRB;OqOJDsK5* z&@A4>yrpppi;jn|m`U{FNM1sbhhrI1*dwpFNtR;OsdG~7KVt|qfm_J@*A#I4o?g?M zIWLlTKbYZAjm1_6jcvPXJ)r-9C7og5!Oh&AR}q&k1m({b*e{4*#8H99F(EEW;7ae&6w3){%=oT2+s$Rps;!iCY@ zLF$clI3{w5?z7xapwKhj^r`s(!&9#e6C<**tmEEFOr#2drfOhY$&*o-Nt}IV@I@77 zl0U_~8ZAph|Apl53lZ6wUDb#!G>i)OBw7vPuvR{r+*ihn64O%!s8{ITc8;%-)g=jx1;;G8VXY?Ym{7W)@J=0O-|i6?isr*rZ9+(&PDxfNI7RLPPD(*A-VcW^);%ugMFpwP(`jh1e&Z$c8O>(v%&`MP^ zottcer9$SJOXX+aEA^pF@=j@JJo|Ez<3|fpN+UUl#r4oI7yfeZZTrpqvRt(<@B`Gx z9y!C#w2ZS^tO1!sQRq~aQWoG%fvBqzgra|fFv3u*yyfsq437}mCzyH!XVz3^C$G_9 zZ3&Zplv_B?-{JZh*y)ZlKAnZ{6}t(vrpBztGw)%Bd(``T--0*!$d^**%L$_{N-Av0 z$8h*(^&W1=dZ5~jvl_e|f`#|0J}aWk*tn==m8kgvmV|dLpFElZWRu$|<)>62BH#D( zf-^si@)PT>yZ6Ba!wyez0aHvnuqE7~-x{eqZalT$Jr6U^C5?#!To^?e;|cq#PGAxW zyBhgYkw9Y;Ke11F1sl>sFLZ3BhxRM2&-@5tFs8lRGvZz*toU!7XVJTb<4%X4b+SkS z<&Zx93SWk46F0W*F;1*83|DK{*@X)h-2Z}4x51EDH#to101|nOoI0fWAp2TD@)O4f zoN)PHR{|8)UN8Tt^5=rEiw5C2qrI%o+H` zJcTQ7u6-aKjE1GxJ{6Cyw*cj?mu=A*2_#KY%&*A^z{r{w>rcW1XkeRvcK@{iPFWNh zewt1KBL9=Wm6i0c@G{jT#gzqHUiPQ@JR`$K)e1}Ky;@ibwLS{z4`6lqx@G(>RVXYa zeFzQefPY!E-$*eWIKFH?d40YYhxxiW>SncY^1AEpgZFK4PV&0ThulF}P8B!Xuj2@v z`=kO6#fRdc%O%NWeQBJPa&`C2bpx{RG&{}bGw`n_U(T za8>%Ag9N7mG?ZO;FlQw z&E5SOzE#gf=N+m;GPg3b_V)!?4h^(w+`&;mQ*!qA^G5?XyAl*=)31;1j)II&f-Y}z zsO*Kvay}gG?N5&DGe^osmy_rhH<-OD$m$oBi2W|jfs`kCID?|j5@jxsUsbHYp!XdI z1utZ~T-_W;dQLdviyw^LmXxVWdjp|9>(4mfpMhp#QLXow9n_hbacT*=<3Ik<=FDz!Lr(jWs;J{+jt#c&HK;-^KwJ^zQ}@mZ{?U=EWKy@uB2 zM4%)U461P4h3*xe)OVs{K=Qd8YL%1#iwe8zhnq<_`emAz+qxO&i(j`--vX@Xczh(H zr-;CoRbZ(!j@V{`8iLlL*!{kMJA(NXbo(qi%bgX5vJ2goX}e~SczDVpn*B3@_Hj|k z_cB=+J^b_Zo`x0d?mNwJWswtdXCuQ;HWgq$vt57O_6dlrjhtb7e-TNYcgJ&HE#hSA z$Fqm{BVoh~`vN8cVdj2GhfVkhE*;O1?c$e1^2IE_`Nn1-^I!N-spE_zKl;?kc}*BI zsa$;Ggdh~XP^jqL%8XN7B1un_7olZmX;APM0+5|shoYHZK(};?=;!ArVCbsN{-BI2 z5Oc;y^yN`CSb5j=p2T@3Va4kK^ipCtnA1Z`Qwpv#44dPW+|4O;SGxq=)oVFK-;hPHDS{$O$$4 zS1Ep7^}QlcW!AXA+KfVP-#)&Y)*qYv!y0H$djyCNMyL`y{9%EX{jF2eMVK`%I=biB zXPgHB}&OI~VE1kcN~=z9gUK5*QJ_>ibj?aZb>QGq_&?$F}tb zn$VcSV&Hz8=d_#rC+d}8QXzm;3mMiY7CJ!X7-4`X*Pv?$W1&vl0~on^$Ch-d8W(m? z&oT(6{m9q(=^+!`b=7E%j^A+ZcW_ba#iQn0)|Hu`vhSX z`mN^j=864bK$fIf# zl}DSw9QgYC{7912`!3kN- z(uc$64?>Wp!-)%7&QMn(A1NHaj-`(U%ynTWQZ>)?zPWM}C`(F=-%8(N4F{DcwEZ8B zULr`9UEPaIvbh28N}bV;xE_?RISp-tEH*Ne@1biTwsLq}2UnUs-0yGt%t~*@MWe%d zFrd#QtHxRj^9#qaW7znyKbm=TM41!F8NM9vjdvp@b!=!kFTa4+Heux8kOnrL-y;!0 zrNado)r5=AoBcp;WUS`#=03W&M|0Zm3`|Juwmcqd4Afjlo7gX!IBh9x_`~Hm{8f3> zX?=PZEVgR%)~k#_!tBGqKLa`iw8o(k@-ImRG{%!Vb+}mZOjp#mHRV^3$dtVWx zq?-0~h>)=@{(*p zYa~}GhT@5FAS&-Xs9-UL%YxZ$Lc$Lq{pO$F%eR7Xt>%K>MKMbp=MFIUb-fEsluI;n zsXbV}ms{(ljwrM}T8O>S)PQMwFOYOIZeVnr!{65~4$$Y-G;6TD$#qvOL!}cpeI)7dJzg;IhhKpn3v$yNjvHACVo;>3@$V&O=!sbj5YhG@BMV41F`Q9hbV3U3v zJUK7Q6yplNGX3t!%m!oamEKc3WX|IF>Hf7&IX&!||5(^@^>G2?dyZg%v`dJlw@<&@ z6NjtG=@P9RK`_oum@t#P4U=M2y(`&y*fw)nKaHyt;>`?gh2tFxXx<4ZJav7Ib&r%x z@*I|-Cg??@zrs9@T}sh%dZU29Rm1!v{9K?WxGR7?))XsvPfG0l{S=3{q#oPb69`Dx z80U5(52llIqVL5$#rbNni}L;#v3dDe>dfe7e~FurJ7{2v6lv9%;6ER6DEV#Q;Gh|l z`iabED@sC?U~A*rwF8)3bnH&84iU3bjl$1sWx`z7QTgk2PB^bRYW_FrF(gP5FU?8q z$6=p71viH}7`gB?)7$Ut-aQVB* z->WAE*FRpv{4Jf8(hfIq;IHXtv*k4yD0dxV3oORrfq%+;eBt}IGV2ixmGG=)$3TiGV}{x@_M zyOd?5@B)=6^JkSc6ze=DDL;!I<9Kz;i}-)?1V-Tn65CrkoKs=mc9U2MEe1wD0;X;- zq^F>w$(o8o#e9dO22-G|XYqxPLOXO(Ba}ytjgctHbSi=VF?I$h6|*I2K%4Zb<_U{% zXy=?e{Zj286kceEHoe-44Of3%WH@gO^T+sO&;Ag=nM&`MU+Oz>l~~jjMG?c=8&Vs`>X|2b_>TQgG>uJ!bwoo}iQ#4ii15;oriI|HJ8xm9)3d-ZzEe z;PS>;&O-zmzQ=qFzi(kz=;5K#oaY3Z**xPvNe^Lp@p1OiwI9&1@MHZ{j|rsIMK4?Y z_CfMtt@dLS7VO*V9NV{V3)HaPRh^gGfh9bzB^INtkQno9PVqw*)b@g-{vUdn%U3^8 zak(4D&1^sLR?Fh*DRXr*No{OqmAlM(@j6uJsV-#6d?e6Fa6Y8|sYBBBK(7pWgpuZ* zN)`vRk*u0yNYP@!0{($-8jiV8vrUy>KJz{FO%N};h>33cq?J{tHy18P`5aj5kivEC z%%qh&!J9m_9VVuRp~hM0B&X(9Toj*~k=#8ArSHyh8WvPy%3~Jc#dnO@;{J2!uUs^c zCj4Xlj%>x1SJw8+v)+*16t!jE&;qCWKTlCvZ^7Eif)-D6IjlKLhx?DLVRlHa%=i=& z)H&a;lFkpsF{2KpS@sejwS75frzf(>0WR&%MI4w}vu{!F{T`gHb51>Kt_NdVZ1+ns zvO}X4n-#m{1f+Z!wLZz~j3vpRKeeTuz=eY0o$X7u(0~2l*AW&0T+t{zC|$db6Mp5{ zO{M0rQQPpo)uj&VMecYCN*uwug7Y87$p~$%31dDsqS%_-`B1=I62^@?MGsvuhY4P- zg`HhyIQ$?f{u#|r=qioSokPT{yT%}uHt5eAy{|Y`ev@}e${N?7K4x?G-VW5_LV*vr{bAt8?Hq$q z39JpBa5rI+z}YR=X&CpJz(^3o4pQD#Bpz_FH{&UUvYn&ute4ARZs1)&I`3&19GGg* zKh%grA3`i?LiDlU>rAxzX8jM6-XzN~ABU2ej*9R5=b-<_fWw-w0g`>2xil9HA%nWy zYMQ`|l$re*wF^@?b0~%5k@Hg+eJ>HEe&92X+3Hwr%<1B6_WpjJx&1g*NjT+B@qh)r zB)jtY8<@8J^K#PiCdZ#2PJX)n6^fh>pYdMvMr!kTU4UdRwAJmVtK6y!-RHjS9~Qg- zz1QFn%jGc`e)D%Rms`W^BpEmevUKu6Tl`hUkgk1%^b*nM2e*5 zgstjT4692A;p6pjn_UT(IO~1li}H94bo4yP&(ukPo-Inr4Asved9kPdaY_+%YTEpW zVg&+?nP+X9HU}1Nxs1Xm)p5NyBei;m0yK~6IDVfNBG4WD75mN88HNULB!mV#BGJ3% zg6Nh_{$IPN@B799i`u(1&#Zbup9S-$qo;U*5)gk`c$5bR?(@qnu6+lp{ELn@AppwL z6n~$aMcDX=V@^=Z3%(I6C6C|SgUy|@tUX^`u+3hart0eoj@YQbVf;=S zu`_GVLE$`*(L8k2gvt+;^de@bX+|J5<_y!fABH3?nc2`GRva1mDweoM3P((TuHTz` z52ZB<33H@gm|hP4biw~H*2g$^jQHuo$a<6L*P2r>@w>XDv_%-k{n;UyRH(Az+%zWMvZOomKS#IiH!<@O}p+?8UV8rq7bMq7G*n8k|@m1^lI5~KkHKk%3 z^h+M11b$A(i6b|Yg_5SAyG4C(w~H+%4F)-;1ja%~XfWUB;Ab!vGuYznN`=Yil=T|+ z09aruT{wJK28aal-{oRIVa(RDNIWbDMuw!fF7vbFxHxA^{k>6WA$_lVtMUyS!-_Ma z8Sg;Oc~*wr_|VfG*mcEhEOCaOJA$Ko+G5M-3;3cU@$MDfgo*Uw`o+ zT4&Nfi@AlAOIA;VgBGBZNgJr#OWjH6s)!KK1+`2fbOECE4IH?aNW)>s!rn+fp#bR=|{X{Fzqw_m+&we$wq|7 zBALqAxGkpETjUYWpUXWXynPo;-h480iDn1RPRVz0JQ2f*H>^xjFW7*RgTEW1?fYC>Xm!Q~7JAYl4AA2*#wn}Pez~XL^U8=92 zVL$h$)dPam0{ZJzn+7H$SYuhO1XvBKXTLn(yV5ob%?@nEu|0kFr5v2r!9gcS4!G8!0 zF}3!CVzpS~UJDOo)PVA6#nn~b8YXP-?e?#ChPL^mOC@$YVO2h?C~-+1Ne75$pO`)@ z*kUkzTx!%FDR(wrZ3$C=c$K1`IZF*N-Yxn|_SH-L8}^+|{d6(@-U(VM;@Qwrp<(%Y2gSigOjq<{G2gC-E$U-05QB%~8s>#yeV^zR z!MdN=nx^J>@8n$|2en@PbNwjJ9#XbB+G~a#ojW==JZ>PhDxi75JR1rUUe#N&j3D)z zSNK*ECc5%rL$a7;_jCUTw zOne04Ze0xc8HwI1X@3LDv8Bs}_Z?v_V$_1@zsU(jL9E)hKVVJi7~6P8nsb zf9}V=#mg7H5=EfnKD+mxvQN?;hX&gW_ouDm`an5%wO2F@?{6Wn$I%qfJmymFdhiQ|b4h0ei}Im1 zuin1JjfQ!TCE*>g$xZ$y-5v;qfgM0u0dbOs7X%i3^vq*`GdA zb3;nbjH5-NEwmjmbBNw!10(0C+2i?JuqH{3mG^4_j4_-(sTie$aq_lntk-WK>DO+x zWs^*tv~4@)(I$Xn_Jk!-yu#gn{((NAB-CMR>_$VJ!B8AJR~%v%&no%O=mloPXhu39kg zF?Z&Hk{A&A9_`khxC3o07ptoUXMt*#`sByRNtn56#jYBuhl4>@4W|XdaP%hk{ylzn zFevS}E*||2Cp^eyD{41kL_BhbZF?7tyB!*nu8J<$!pR#|cOP;1+bUP3%Q%$0NM%pY z$-|nPb=12IVc3;Ug~8jUoO_GJ@z3AS-!@)6$HC5>lUh}3P@0=2 z8mFfOlo38c;^PNULJ7E;LePNaPG>WFq76*wIyYEdjDjvXH@B=^6wLW>imSne6-hU& zw>{yt#Gxd`Mo&X}T#^^mKQY!qpcz+qQFllcaVguP6=5quE@VlhZL$|La&@ ztJ+>9Zmuqu6Vg0#DOa(g_PVbz!-X%|jRb-j_#sbI z*m!G}Hx>r8a*l2KUMuJEsI$EPecyZ8zzeskrG9AXdJ;#L4uh4|U!~G7_>mk=Bes8O z8?4A0J~nYp2TI{>!QXez!;~0}otI<{_EquG{)^g;)MFt%YWtXA_(M3y#{FKXlj3Er z>AL{YOf<-Z@ejF%!hWqJ$F5xm)5jTE#SwKfu={7)PkEClc z`;)B4u}fe|?jO@Efu=pxtMC4Hr0!DL#x_)p-wrh^GTr};acl!-@^pq6)jzk7x#uVn z^=`=qe<9#M7uSC3n{PlGES=xgeiAzjs%Cw{G;m$5732yRV4Ysm@q^+n%=I5ijXI-( zE4hYU1N*h0>2(Kxr-&pF@1C`moo&MuJ^S3V_hXTiXD(Z~YYN*7FDX^-{s=ulLp?De zg_AN(;SMkVZQl0@WlLren7DSmwt~YLn%};hbr*2JX^xm4Cf7OaP5LLD@^}TRRiuAQ zM1MrzqaJp2_YPw!>uJTb&g}&>66(#FtbCjMG5l_G#BLz($<#?aFN?JW#y^e-3*bZu z-D=UpGHku|TZzn8gl+tzh;TW|cHT74Kca{*oT|XN>09 z*Zh{}UWc^kuJQNKd8nEHJv$UcR#;QU|^cx==hWk0|L>rYEte<8Sv&CN-V?C7cH`mCJzhf znj>BLLJ^9^HO$x(rg0%t#q5OW9Mmk@nmyKCg|Q=JVSJ31&{3gK@l@{!&eJN65BxPj zvMH0^Sv6lAt*Salr~e#R5_a+&*~5dxXReIN6~4Ikv98#1ZU9s5_9PpIWx@3B=?ZNp zIw(FsGd#IV5=O2bJ>6OngjM^#eoVh14s)XkdRZiWT~FVBXeJA zG~2;)O?K|`Q4RRJ_~LwCCIT_uE;@*C5qlzP4c_>c!bs@r3|E?PI2xCddwWJ0+qr0i zXQEz0#%i$7rIz2=tnRm8dzU6Iv+R4A^LG@dXXrV+CO^Uu-Gx8x0xavTwb> z{++hvy+1}^*ew3MDf1m@ZL730xbYJyZ@+FZ<2x+Ln{$!=$O(OqO=3Kzw%|mO1NZIi zv$*l})d!}bBd{#+;!m9ABLZDRdH=O%OR%sg%3(EBjif_4REb$1Xb8EXPLuT)NvS5h z^5JDbWXSr_=w<+OXW&e{c{uhT+@HC!eF>IN-)sKn_YDSVkBQN4yvD)jl6hCujIJ{G}h0>2sPXzqlK)dFToMT&QaXNyhIzdDp zDjDw$pZXRABZ2yMhaXX3OkLYbypjVJPdyu5*QkT>vrcvw1GQi#hb(0~Hh?5au~+ZX zWS~bVX=sbG15hQrx0!kT!t`>!^3;_KZ1!yTDPYdW0q(UBsFtxa$v!%bg$>&K!&Dw)=Yr(uBH=TVMt->erHg<$d|LQ) z;u==8d+ZSF7QmE#2I?M%Ye4PV;;hYZ2gtdm7eAYc=+p2N(J! z#7{HAXx(WcA!lct7_XbW`o|Y1u6KAY()dG*Z01QbwSF9|Mh5q*Rk)t^=X_A+6Su@p0i8{0=LpQ^46XNyeIEk_EIs(*C{wLWVX5L$WuIDXeeoES^u!gSqrH>42l& z=<0s@f&ZZzTslB*O+CE`ox`<4aWY4NEI)Ly*5Vx2H;0=&nw^L0Ql5e4haZvJrB5(2 zy@ZR_`|paizr-*uUdB6*OmY63-FkP#Zzw-VHE%X^fMw+uWpa@*SZ7hI@<@jTJ>w%D zJjgJ^KOc>B?W5g+c#{6Repx!K9`^U?pJvCd@@y4@qF|)#&YV8?dl{*}G`3VUT43=+ zUHz)23xOuoaCK@IVoJ*0PaF-uu!8C7ydz^35Z4}Vz59|5dN1$WlXb=rVs6}Iy)M=Y z#f&8rGH)M32ey`;o?wNR0sc#>ePoynA$ogW-U@RBarZA353z6Sfi&&jY#ie57u)|- zAE_p7Hm~&ip)P~)lE=ybXT>6qr0R3T4c1e_{6KKx7;hAGz!)0E&Kn30cm?pG7 zlDBDr`js|wv$y|%6rp=6zx*GP{x(ObF*aiFla9XcwW>h$QX0EsOToqPlgu8+6tTru z!zz&02IfLZQ(R*|u|AWX32NhTvDEZpyM#{x(^yZQ*NGA&69`8`Z5m;%p4A~Ddkxla zux{zDo<*|8MzT&o3sBSBxtW>vz<|Yp?e?~#Flt}7qH5G8L9CS4Gy*%=s56M~fQ724Z0a;tgkiClqD+jd3Sr;Rr$5-ZN zJSPF0qa^u5l{R2>F^kh|)DSzimkn;GgcdN&WJGGzJ0n?5;k(*_Sm^h7dcbkI6P9im z%FF7|0Wr;@?k2~37$M)d#sAR(EXF=n2KkF%Pm`ZP)uJ)x%ay)lK16JKFi}(`1N%d z8|fcgcn>UK^F;UlY*QQT@)9jDJ7|To*O^>A_YYuYaequVg#r2w&59^ZXah;!Hl0^s z7lAJ6%L^^VQ?T@?m1Jxzi0cO)rdXUS#3gzoZK-SLA@=&GnZUh2keGO93(?jaJ7cmE zgWGA~bJx$@8}a!CtOuGi+2V_^_KEjSqRJtlm@wb?lEV$7BJb#;r2q?c_f1Oux&l<= zrFNq8RdhWy-78O%iT!Ppk6R)E%0K_9y>mz!I{MujJJok$id^wet24i#=Vx`I!J#~$ zs>Lel7{qS!uhIOQmCd?k_iVe-B?M$r)zO8bPdKmn@#nMQMTkq+&#Ag`4#_f(KW?J| zta#>@$a)S!5vQy1l>#!9zL(+HXUvC_o?j{)_r8ak9OF|rKD>v8Pp|hzyo!NhHEm6j zR~`KF{_s1L1lSo=U*eD{218-Y2@hvdflT+Oy6TDzE{=t3GCJJ`iubP&GZiXMio`NC z6#sx}9VcF%fQvw8Yshx@eFLpZ;ptr`5ju81QoG5i2hHc#*rkOUp+$l_R>b53)-)SG z`Z~c5q>XRiLeljieexvF>smc16{RFEckRF>6ImL??i~ag%iGcN1}oUkm&j9Eesa@4 znIBr_tKdXJ;)&8x37k~jyC8c)4#HOB<6=1Nf!s9NrX6OBqY0PZj{c*;k_4vO-;BFp z;gxLAgtR|WUYO;ZTQ}m2uy5aN?lxRD+)nRR9|gZlt_KvWCE$EdJM)$0QTVC2qgem$ z0nESgdFsec9q2l~5*3;U6{%l`>TQly zdQYXWzBoXg_JcK0LhtFRE*!(5W@41(bRCc;bV-i)M6puUT-9l94vLiTZ7~cBg1$@k zRj)oh0b={pQxEgRq2V-l<-Wp?Fn9V?_`kb)SlBOl%gWamhrZU?aU?&0<@*`mIxh5L zpLlfY!ebX8+Q%(Cf2#-$O;KUW{YPLbgT*N7*Lk1_G;FOMwuY`J{U>T(y@aNSt9RCv z(*DEgjk~SH0$+Z?I#apu4dyTc?cE2r1GPkJM@KMB#$Gm;O*~!9L zEM$yuZY7n#aEm<$&!wF>p4Cv=q?n6jVPlmZe`X9>zoeP+)f5+I=pHMtGZPp-wl`F5 zzkus!S)_BHTfrDbM!x@TI7FYV75c}?3zHcWakBf3Vff$Jmi@2n23kpVRC^YHZ%->l;aQ%8nZ51?^j=vLp;P2XYQ ziikDifdYjKL$73i;aF;I#=%=>pp8GSyVdLrj^CJy{Pgq(uC5(^=NaM(lmwwA1?C~3 zcxenvo(o3mdTw()`6@2>a=ww>>4=lrBV`rEhafOjIK5|_7XQhrhsau5LzNffYlq*5 z3+Sx6EL$!I!qCcinyYhXNL%t%?|ITY*vFn3TAf|~X^d{0-6<5An~IJMsaB~Lo!N4Wxl z?!utaos%6nx<`q_F28hhpDg$d)iT4tjPKMQb#E9Hd6KjCdkK2JJ2aPf+d=75E6Pxo zIxMZ+y8O}NJq&!_cYf94{1gEYhiOXj0F#c?4ioMVWX6cm$l0`#cFgTLyD%(G7 zDreAp^uq_{%1-=z6iF>$*m@>==NUDikO~HTIMuN3PyMGr>nqs1_h|V0!_T0_+L&Lf zlnv&eOnT-C9K-x0{JKswBiOP__#dNTI%J&Jy6EpA3FSPC!HUm`a9_ow&bX5S23me^ zS5n!Jq^%;v2+L;$3|F(BR0#RPLijV`y&3`7FDPl!Uml6WI}0Q_#`v)7bb82=dqM#X ztt40b2U{4k6ihZe{0JsaLv^#=!*i7<-T`}D)Sl)`pu z?_HcMdTS@+_y$QQRp}4>ZGv$n`qz0D2XXw{u7c}UWiZ9bW_ZHg4$~$&QRA(0I|+NA@^l`_o@gTYC+gbK4Jn_<0{% zLC5@X2M9E&CTYnJVC?vSbZF#4yrKlDl^Mb0YnLEHsvOPlO( zp!15&*yBjHjVC9n$bbh&Qs5?(1*dI#VK~} zTVzT`=pFdNZW)a)LFln<#CjY~C*HjL-XadANTaM#c*TYNt$GV`meT$>#Uv;b<3jC*1e*-5u9CJRT&_H#>)3J(R zZ)j_;w&y6F2C~gi*@xXe&`9H~!lTaw8#~GduKh8AJR^de^3@kO+!5xLnbm;9(p|*C zy?2q+(E75zF%s9#=5<)sOhB`tE6BZB#^qatb3B$5D7ARArd!>M9rT0F_xR~y?t0Co z9P1P$>k-VVlT|T$oaL#Jh%OeepI8wZ9Ri}ruZxR6$WR>`tlsr-5y(Qa=j)PcVBtRL zcJ^c{v?_cxY2hJ3C)bRTwBLK6y0}|P&+NxhsUX*H8y!$v_@lZeVGq2Z7HVw?GMy$)4#EF9REoKpG`@_t2yuk7K%m-?aYS`miOnhLLodpuF#cPqLTCi$Qm%G&^FKoT{mxJru zEF`nGbKfNw5@>jW9x$*MZ~E6Or5MxF&73tV{yW|dTjz=d>9vm{$)ZBRdK(Fe@n(jc zjwSd*K;mna+$EUcPi}tIn}j2W)m$`IUSZF?VUOtbHS9mp&FEoA1oE!@OK*=O%pMze zT-{QDZBMu)mV0NQ?ZV)v??3(89U>_K{ z(X_aCx)&z}=AYb4jj%Qk7QoQ&4REirm-R%B&fTIH&iO0-2B_ok6AG9J_ zdk_cS)Sl~6CPB}eWAP0(OGs{x@z7)AgR-`{U)FI)kdmYM^mg$CP@dRl2shk>p{;JQ zReYDR+{ZY&X>=4RlP9#q-Y;P9?LD;jS?Q5-Hb(y3VGUT%V$F8XV#fKePChL|UO1uY zVPjj37KE1lp^e3C70~#ihJDTW7cRz9?rHL#!A`M7+lSk9G3^Y^Zbh2j z0)_|u%YH7Wpy7;Xj;@d-tU4R4GKoBb0WnRrhX*oo@obGM&#Bvxe{yDrLFINB6tE_E z1s#Rivg0AY#At9P?b#tQc}pB)=-{*6nuZH`p|U)dsYpy_csr=%57Un*2fUK};H_a* zrcXpK4D9DGHCRsqYVh@L-T^vj+Epv_Xe0+oJ05hzRp`L*%hDw67EbIxrTm~ob2rWu ztV}c99fRNOiAJ_OWw1c2UAARL9f_g0T}h{^aGq8vP18IA2Nmk;-#gVqX?}S0LVZ&K z4XNv^;&w}D=UbK8bWj`uP0!i$c~B`FG_>04jm=Gs46kKsa3X7q;q1ggY$Co8r3Mpl zfn9vp(XMJ3)fbQTi!XxJCzL`178Q8&i>MOxX%;&&ip>t`PGD06`)IbVHjrJ+;wKo` zas9f`^KxEKoD9rkJUhM$%}=(ZNs}&O`{ebBcG(9|-Jl*CXT0e<23uY88?-TBZQ%g5 zMg(WKH8e7pF5{d~pI7mECH%eqKI*%#AZ&!#sxzBc!`#ZT68F#t*y5fLb}CN-rk%-} z&y1vDG5#poDQh!#efM1F>gB{P2Yz<8DJvY?SoiwzR2_#RXuntQiiZrt_A^e}KXCc2 z74Hv?n9V+{aA|6X0F+o?e;mxB1YJ$MXRjD@L4&mYfMbVQ0ey%{*i>&btkX>4c-Uub z))hAdbBU^jEskX=4IK)hhe?gSqP?Ud zB*vOJ7^qAZu-dW^%MM9Fxy!*75APk=%g3xF`L+?(v`p8t6t zjbK96mlsA6J~;HGu$mcCa9~UJm5F)_h)8(ldTeVVZV11PtYUeO(=Su2&yJd4{kMsd z1#2k+&74_F>`Nz@ZQo8JwK8D*@hkEdn7+WU=>rbRf)=hQ{AH;=-G<}pw4UXz|6pQh z7mw+RAl~0NzTYo&lh1t)Gt%n_!UT(C+-0xR&_`DE)8nHf&?sLWK5_makm8I^tLQ(# z^)jl2HaQ=c?vK(wVls#EA9aOVC$=Nyu49BtZ3j*+H}6uorwzmH_iy%`8!w;@*Gaur z#01?n?Q|_$>|k2Z=G@L9Wr#1j@wCAq3>rtJr+8NP0x6xdM9M4|EI5fQ1i@{%5IHsB zIg$u{Zig(G#D8O5-_Pks_KaBUpt}EJk^vT_f&9p-3=H2h3HdgD36?Wwj2Hc@F`e9N zaA`;zma^^8m8}j({UcKf{smMGH!e2`9y~J{GUFdHiGRgj z?Bg{gOT>_{3fI%T6<4y4~dqdSf#!;|V8{pr_lQ<-tXtgG>8i*gL;j#NiG8{iuDm$n`jurab*YU;h}AU&$-9 zi&#O*LUb?pcUf$368`GDkqWb%6T4djy)oY=|5XINA-3K3e)L$z7Abx?t0%H$kZhls zXE^)^R}HJ9+@KMLtyPm_G>VaI@^MGs*A^(}?zs8VZWYRZ?Uvar^vyiYy#LQl0uJ1Y zrTF!1V1F{T;?mJJ7>j#OTjnQ%)-@(MCJ#$s^+kp~k>y&4I7yFI1-@Myxif!AD9MAIj1XGtB4NeCs z!oa#r=j^E<7_Z?Jn@?4Q$!tw(_}AOeW*>-sU3aiYaegS|crMm|4XxQ=DMX|B-CNJD zB*XHv1wpa0PXt=74LkOz1?+$IWOvU+YgoP7^~s62iv3njH5F3g&@|8SiRXnZ)=5dw z%E`CE6rZUYomnzW%n8W+sk;Znv+t`kt`*^*Hd>MT@4td`d=ykUwQ1m^7wzH!hhv|N0a^UXtFZ@|ICda8Cw6}%>VBq>pjC~hLR>vHy7pq o`#a0i$F`FHw^$tiGb^Xf!^Xxs`{wC?o>Ex Date: Mon, 9 Nov 2020 23:02:35 -0500 Subject: [PATCH 09/25] updated arguments in `image.bas` so that user could over-ride hard-coded options (names and margins) for better formatting of plots and added unit test close #43 --- R/image.R | 31 +++++++++++++++++++++++-------- man/image.bas.Rd | 15 ++++++++++++--- tests/testthat/test-image.R | 8 ++++++-- 3 files changed, 41 insertions(+), 13 deletions(-) diff --git a/R/image.R b/R/image.R index 8706b5a6..3f30477b 100644 --- a/R/image.R +++ b/R/image.R @@ -28,11 +28,15 @@ #' @param subset indices of variables to include/exclude in plot #' @param drop.always.included logical variable to drop variables that are #' always forced into the model. FALSE by default. +#' @param namesx character vector for predictors. By default, this is NULL and +#' names are extracted from `x$namesx`. Should be the same length as +#' p + 1 (including the intercept). #' @param offset numeric value to add to intensity #' @param digits number of digits in posterior probabilities to keep -#' @param vlas las parameter for placing variable names; see par -#' @param plas las parameter for posterior probability axis -#' @param rlas las parameter for model ranks +#' @param vlas las parameter for placing variable names; see `par()` +#' @param plas las parameter for posterior probability axis; see `par()` +#' @param rlas las parameter for model ranks; see `par()` +#' @param mar mar parameter for margins; see `par()` #' @param ... Other parameters to be passed to the \code{image} and \code{axis} #' functions. #' @note Suggestion to allow area of models be proportional to posterior @@ -55,8 +59,12 @@ #' @family bas plots #' @method image bas #' @export -image.bas <- function(x, top.models = 20, intensity = TRUE, prob = TRUE, log = TRUE, rotate = TRUE, color = "rainbow", subset = NULL, drop.always.included = FALSE, - offset = .75, digits = 3, vlas = 2, plas = 0, rlas = 0, ...) { +image.bas <- function(x, top.models = 20, intensity = TRUE, prob = TRUE, + log = TRUE, rotate = TRUE, color = "rainbow", + subset = NULL, drop.always.included = FALSE, + namesx = NULL, + offset = .75, digits = 3, + vlas = 2, plas = 0, rlas = 0, mar = NULL, ...) { postprob <- x$postprobs top.models <- min(top.models, x$n.models) best <- order(-x$postprobs)[1:top.models] @@ -75,7 +83,12 @@ image.bas <- function(x, top.models = 20, intensity = TRUE, prob = TRUE, log = T which.mat <- which.mat[, subset, drop = FALSE] nvar <- ncol(which.mat) - namesx <- x$namesx[subset] + + if (is.null(namesx)) namesx = x$namesx + if (length(namesx) != x$n.vars) { + stop("length of namesx does not equal number of variables") + } + namesx <- namesx[subset] scale <- postprob prob.lab <- "Posterior Probability" @@ -113,7 +126,8 @@ image.bas <- function(x, top.models = 20, intensity = TRUE, prob = TRUE, log = T if (rotate) { - par(mar = c(6, 6, 3, 5) + .1) + if (is.null(mar)) mar = c(6, 6, 3, 5) + .1 + par(mar = mar) image(0:nvar, mat, t(which.mat[top.models:1, , drop = FALSE]), xaxt = "n", yaxt = "n", ylab = "", @@ -129,7 +143,8 @@ image.bas <- function(x, top.models = 20, intensity = TRUE, prob = TRUE, log = T axis(1, at = (1:nvar - .5), labels = namesx, las = vlas, ...) } else { - par(mar = c(6, 8, 6, 2) + .1) + if (is.null(mar)) mar = c(6, 8, 6, 2) + .1 + par(mar = mar) image(mat, 0:nvar, which.mat[, nvar:1, drop = FALSE], xaxt = "n", yaxt = "n", xlab = "", diff --git a/man/image.bas.Rd b/man/image.bas.Rd index 1075f1fc..df352321 100644 --- a/man/image.bas.Rd +++ b/man/image.bas.Rd @@ -15,11 +15,13 @@ color = "rainbow", subset = NULL, drop.always.included = FALSE, + namesx = NULL, offset = 0.75, digits = 3, vlas = 2, plas = 0, rlas = 0, + mar = NULL, ... ) } @@ -53,15 +55,21 @@ white image (greyscale image)} \item{drop.always.included}{logical variable to drop variables that are always forced into the model. FALSE by default.} +\item{namesx}{character vector for predictors. By default, this is NULL and +names are extracted from `x$namesx`. Should be the same length as +p + 1 (including the intercept).} + \item{offset}{numeric value to add to intensity} \item{digits}{number of digits in posterior probabilities to keep} -\item{vlas}{las parameter for placing variable names; see par} +\item{vlas}{las parameter for placing variable names; see `par()`} + +\item{plas}{las parameter for posterior probability axis; see `par()`} -\item{plas}{las parameter for posterior probability axis} +\item{rlas}{las parameter for model ranks; see `par()`} -\item{rlas}{las parameter for model ranks} +\item{mar}{mar parameter for margins; see `par()`} \item{...}{Other parameters to be passed to the \code{image} and \code{axis} functions.} @@ -101,6 +109,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{plot.confint.bas}()}, diff --git a/tests/testthat/test-image.R b/tests/testthat/test-image.R index ef9bc015..939dee8e 100644 --- a/tests/testthat/test-image.R +++ b/tests/testthat/test-image.R @@ -5,8 +5,12 @@ test_that("test image plots", { hald.ZSprior <- bas.lm(Y ~ ., data = Hald, prior = "ZS-null") expect_null(image(hald.ZSprior, drop.always.included = TRUE)) expect_null(image(hald.ZSprior, drop.always.included = FALSE)) - expect_null(image(hald.ZSprior, rotate = FALSE)) - expect_null(image(hald.ZSprior, prob = FALSE)) expect_null(image(hald.ZSprior, intensity = FALSE)) expect_null(image(hald.ZSprior, intensity = TRUE, color = "blackandwhite")) + # tests related to issue #43 + expect_null(image(hald.ZSprior, rotate = FALSE, mar = c(6, 8, 6, 2) + .1)) + expect_null(image(hald.ZSprior, prob = FALSE, rotate = TRUE, + mar = c(6, 8, 6, 2) + .1)) + expect_error(image(hald.ZSprior, prob = FALSE, rotate = TRUE, + namesx = c("Intercept"))) }) From c47f2d8a1be44814e2d74d8c7173f28c0aa445f0 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Mon, 9 Nov 2020 23:26:10 -0500 Subject: [PATCH 10/25] * missing `mean.x` in GLMs caused `predict` to error when `se.fit = TRUE` [Issue #50](https://github.com/merliseclyde/BAS/issues/50) * Prediction under the MPM failed with GLMs [Issue #51](https://github.com/merliseclyde/BAS/issues/51) close #50 close #51 issue #52 regarding SE still open! --- R/bas_glm.R | 22 ++++++---- R/bas_lm.R | 2 + R/predict.R | 75 ++++++++++++++++++++++++++++------ man/bas.glm.Rd | 22 ++++++---- man/bas.lm.Rd | 2 + man/fitted.Rd | 1 + man/variable.names.pred.bas.Rd | 1 + tests/testthat/test-predict.R | 68 ++++++++++++++++++++++++------ 8 files changed, 154 insertions(+), 39 deletions(-) diff --git a/R/bas_glm.R b/R/bas_glm.R index 3458f53f..0361edd3 100644 --- a/R/bas_glm.R +++ b/R/bas_glm.R @@ -159,18 +159,26 @@ normalize.initprobs.glm <- function(initprobs, glm.obj) { #' \item{priorprobs}{the prior probabilities of the models selected} #' \item{logmarg}{values of the log of the marginal likelihood for the models} #' \item{n.vars}{total number of independent variables in the full model, -#' including the intercept} \item{size}{the number of independent variables in -#' each of the models, includes the intercept} \item{which}{a list of lists +#' including the intercept} +#' \item{size}{the number of independent variables in +#' each of the models, includes the intercept} +#' \item{which}{a list of lists #' with one list per model with variables that are included in the model} #' \item{probne0}{the posterior probability that each variable is non-zero} #' \item{coefficients}{list of lists with one list per model giving the GLM -#' estimate of each (nonzero) coefficient for each model.} \item{se}{list of +#' estimate of each (nonzero) coefficient for each model.} +#' \item{se}{list of #' lists with one list per model giving the GLM standard error of each -#' coefficient for each model} \item{deviance}{the GLM deviance for each model} +#' coefficient for each model} +#' \item{deviance}{the GLM deviance for each model} #' \item{modelprior}{the prior distribution on models that created the BMA -#' object} \item{Q}{the Q statistic for each model used in the marginal -#' likelihood approximation} \item{Y}{response} \item{X}{matrix of predictors} -#' \item{family}{family object from the original call} \item{betaprior}{family +#' object} +#' \item{Q}{the Q statistic for each model used in the marginal +#' likelihood approximation} +#' \item{Y}{response} +#' \item{X}{matrix of predictors} +#' \item{family}{family object from the original call} +#' \item{betaprior}{family #' object for prior on coefficients, including hyperparameters} #' \item{modelprior}{family object for prior on the models} #' \item{include.always}{indices of variables that are forced into the model} diff --git a/R/bas_lm.R b/R/bas_lm.R index 294c856f..e2af5e4c 100644 --- a/R/bas_lm.R +++ b/R/bas_lm.R @@ -302,6 +302,7 @@ is.solaris<-function() { #' \item{X}{matrix of predictors} #' \item{mean.x}{vector of means for each column of X (used in #' \code{\link{predict.bas}})} +#' \link{weights} used in model fitting #' \item{include.always}{indices of variables that are forced into the model} #' #' The function \code{\link{summary.bas}}, is used to print a summary of the @@ -844,6 +845,7 @@ bas.lm <- function(formula, result$xlevels <- .getXlevels(mt, mf) result$terms <- mt result$model <- mf + result$weights <- weights class(result) <- c("bas") if (prior == "EB-global") { diff --git a/R/predict.R b/R/predict.R index 7d8098c4..9ce4e679 100644 --- a/R/predict.R +++ b/R/predict.R @@ -230,10 +230,13 @@ predict.bas <- function(object, estimator = "BMA", na.action = na.pass, ...) { - if (!(estimator %in% c("BMA", "HPM", "MPM", "BPM"))) { + if (!(estimator %in% c("BMA", "HPM", "MPM", "BPM", "MPMold"))) { stop("Estimator must be one of 'BMA', 'BPM', 'HPM', or 'MPM'.") } + if (estimator == "MPM") { + object = extract_MPM(object) + } tt <- terms(object) if (missing(newdata) || is.null(newdata)) { @@ -277,7 +280,7 @@ predict.bas <- function(object, df <- object$df - if (estimator == "MPM") { + if (estimator == "MPMold") { nvar <- object$n.vars - 1 bestmodel <- (0:nvar)[object$probne0 > .5] newX <- cbind(1, newdata) @@ -415,12 +418,11 @@ predict.bas <- function(object, if (se.fit) { if (estimator != "BMA") { - se <- .se.fit(fit, newdata, object, insample) + se <- .se.fit(fit, newdata, object, insample, type) } else { se <- .se.bma( - Ybma, newdata, Ypred, best, object, - insample + Ybma, newdata, Ypred, best, object, insample, type ) } } @@ -521,8 +523,12 @@ fitted.bas <- function(object, if (is.null(top)) { top <- nmodels } + + if (estimator == "MPM") { top = 1 } + + if (estimator == "HPM") { - yhat <- predict( + yhat <- predict( object, newdata = NULL, top = 1, @@ -530,6 +536,7 @@ fitted.bas <- function(object, na.action = na.action )$fit } + if (estimator == "BMA") { yhat <- predict( object, @@ -548,6 +555,15 @@ fitted.bas <- function(object, na.action = na.action )$fit } + if (estimator == "MPMold") { + yhat <- predict( + object, + newdata = NULL, + top = top, + estimator = "MPMold", type = type, + na.action = na.action + )$fit + } if (estimator == "BPM") { yhat <- predict( object, @@ -558,10 +574,19 @@ fitted.bas <- function(object, )$fit } + yhat <- predict( + object, + newdata = NULL, + top = top, + estimator = estimator, + type = type, + na.action = na.action + )$fit + return(as.vector(yhat)) } -.se.fit <- function(yhat, X, object, insample) { +.se.fit <- function(yhat, X, object, insample, type) { n <- object$n model <- attr(yhat, "model") best <- attr(yhat, "best") @@ -569,7 +594,7 @@ fitted.bas <- function(object, df <- object$df[best] mean.x = object$mean.x # glms don't have centered X for intercept so need t - # to center X and newX to get the right hat values with orthogonal X + # to center X and newX to get the right hat values with weights if (is.null(mean.x)) { mean.x =colMeans(object$X[,-1]) X = sweep(X, 2, mean.x) @@ -578,9 +603,33 @@ fitted.bas <- function(object, shrinkage <- object$shrinkage[best] + if (insample) { - xiXTXxiT <- hat(object$X[, model + 1]) - 1 / n - } else { + + if (!is.null(object$family$family)) { + if (type == 'link') { + mu.eta <- object$family$mu.eta(as.vector(yhat)) + weights <- mu.eta^2/object$family$variance(object$family$linkinv(yhat)) + } + else { + mu.eta <- object$family$mu.eta(object$family$link(as.vector(yhat))) + weights <- mu.eta^2/object$family$variance(as.vector(yhat)) + } + } + else { + if (!is.null(object$weights)) { + weights <- object$weights + } + else { + weights = rep(1, object$n) + } + } + # browser() FIX issue #52 + xiXTXxiT <- hat(diag(sqrt(weights)) %*% object$X[, model + 1])/weights - 1 / sum(weights) + } + else { + + #Fix below! FIX issue #52 X <- cbind(1, X[, model[-1], drop = FALSE]) oldX <- (sweep(object$X[, -1], 2, mean.x))[, model[-1]] #center # browser() @@ -589,9 +638,9 @@ fitted.bas <- function(object, } scale_fit <- 1 / n + object$shrinkage[best] * xiXTXxiT if (is.null(object$family)) { - family <- gaussian() + object$family <- gaussian() } - if (eval(family)$family == "gaussian") { + if (object$family$family == "gaussian") { ssy <- var(object$Y) * (n - 1) bayes_mse <- ssy * (1 - shrinkage * object$R2[best]) / df } @@ -607,7 +656,7 @@ fitted.bas <- function(object, )) } -.se.bma <- function(fit, Xnew, Ypred, best, object, insample) { +.se.bma <- function(fit, Xnew, Ypred, best, object, insample, type) { n <- object$n df <- object$df[best] diff --git a/man/bas.glm.Rd b/man/bas.glm.Rd index a35e08da..298411f0 100644 --- a/man/bas.glm.Rd +++ b/man/bas.glm.Rd @@ -164,18 +164,26 @@ components: \item{priorprobs}{the prior probabilities of the models selected} \item{logmarg}{values of the log of the marginal likelihood for the models} \item{n.vars}{total number of independent variables in the full model, -including the intercept} \item{size}{the number of independent variables in -each of the models, includes the intercept} \item{which}{a list of lists +including the intercept} +\item{size}{the number of independent variables in +each of the models, includes the intercept} +\item{which}{a list of lists with one list per model with variables that are included in the model} \item{probne0}{the posterior probability that each variable is non-zero} \item{coefficients}{list of lists with one list per model giving the GLM -estimate of each (nonzero) coefficient for each model.} \item{se}{list of +estimate of each (nonzero) coefficient for each model.} +\item{se}{list of lists with one list per model giving the GLM standard error of each -coefficient for each model} \item{deviance}{the GLM deviance for each model} +coefficient for each model} +\item{deviance}{the GLM deviance for each model} \item{modelprior}{the prior distribution on models that created the BMA -object} \item{Q}{the Q statistic for each model used in the marginal -likelihood approximation} \item{Y}{response} \item{X}{matrix of predictors} -\item{family}{family object from the original call} \item{betaprior}{family +object} +\item{Q}{the Q statistic for each model used in the marginal +likelihood approximation} +\item{Y}{response} +\item{X}{matrix of predictors} +\item{family}{family object from the original call} +\item{betaprior}{family object for prior on coefficients, including hyperparameters} \item{modelprior}{family object for prior on the models} \item{include.always}{indices of variables that are forced into the model} diff --git a/man/bas.lm.Rd b/man/bas.lm.Rd index ceecf8ac..08a98685 100644 --- a/man/bas.lm.Rd +++ b/man/bas.lm.Rd @@ -259,6 +259,7 @@ BMA object} \item{X}{matrix of predictors} \item{mean.x}{vector of means for each column of X (used in \code{\link{predict.bas}})} +\link{weights} used in model fitting \item{include.always}{indices of variables that are forced into the model} The function \code{\link{summary.bas}}, is used to print a summary of the @@ -442,6 +443,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/fitted.Rd b/man/fitted.Rd index 38d8990d..fde3ced8 100644 --- a/man/fitted.Rd +++ b/man/fitted.Rd @@ -84,6 +84,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, \code{\link{plot.confint.bas}()}, diff --git a/man/variable.names.pred.bas.Rd b/man/variable.names.pred.bas.Rd index ca0de94f..6fd743a1 100644 --- a/man/variable.names.pred.bas.Rd +++ b/man/variable.names.pred.bas.Rd @@ -48,6 +48,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/tests/testthat/test-predict.R b/tests/testthat/test-predict.R index b037fad7..f5d62174 100644 --- a/tests/testthat/test-predict.R +++ b/tests/testthat/test-predict.R @@ -36,18 +36,42 @@ test_that("predict.bas.glm", { pima_pred <- predict(pima_gprior, estimator = "HPM", se.fit = FALSE) + pima_top <- predict(pima_gprior, + estimator = "BMA", top=1, + se.fit = TRUE) + + expect_equal(pima_pred$fit, pima_top$fit, check.attributes = FALSE) + +}) + + + +#Fixed Issue #51 +test_that("MPM and predict glm", { + data("Pima.tr", package="MASS") + data("Pima.te", package="MASS") + pima_gprior <- bas.glm(type ~ ., data = Pima.tr, + betaprior = g.prior(g=as.numeric(nrow(Pima.tr))), + family=binomial()) + pima_MPM = extract_MPM(pima_gprior) + + expect_equal(predict(pima_gprior, estimator = "MPM", se.fit = FALSE)$fit, + predict(pima_MPM, se.fit = FALSE)$fit, + check.attributes = FALSE) + + + pima_pred <- predict(pima_gprior, + estimator = "MPM", type = "link", + se.fit = FALSE) + pima_fit <- fitted(pima_gprior, + estimator = "MPM") + + expect_equal(pima_pred$fit, pima_fit, check.attributes = FALSE) - # should not error - expect_error(predict(pima_gprior, - estimator = "HPM", - se.fit = TRUE)) -# expect_null(plot(confint(pima_pred, parm = "mean"))) -# should not error - expect_error( predict(hald_gprior, newdata=Pima.te, estimator = "HPM", - se.fit = TRUE)) - #expect_null(plot(confint(pima_pred))) }) + +# Issue #52 SE's are incorrect for glms and weighted regression test_that("se.fit.glm", { data("Pima.tr", package="MASS") data("Pima.te", package="MASS") @@ -57,6 +81,7 @@ pima.bic = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7, betaprior=bic.prior(n=200), family=binomial(), modelprior=beta.binomial(1,1)) +fit.bic = predict(pima.bic, se.fit = TRUE, top=1, type="link", estimator="HPM") pred.bic = predict(pima.bic, newdata=Pima.te, se.fit = TRUE, top=1, type="link") form = paste("type ~ ", @@ -64,13 +89,32 @@ form = paste("type ~ ", collapse = "+")) pima.glm = glm(form, data=Pima.tr, family=binomial()) +fit.glm = predict(pima.glm, se.fit=TRUE, type='link') pred.glm = predict(pima.glm, newdata=Pima.te, se.fit=TRUE, type='link') -expect_true(all.equal(pred.glm$fit, pred.bic$fit, check.attributes = FALSE)) +expect_true(all.equal(fit.glm$fit, fit.bic$fit, check.attributes = FALSE)) + # issue #50 in github regarding se.fit failing; debugging indicates se.fit is # incorrect -# Should be expect_true -expect_false(all.equal(pred.glm$se.fit, pred.bic$se.fit, check.attributes = FALSE)) +# Should be expect_equal + +expect_equal(fit.glm$se.fit, fit.bic$se.fit, check.attributes = FALSE) +expect_equal(pred.glm$se.fit, pred.bic$se.fit, check.attributes = FALSE) + + +}) + +# Added feature issue #53 +test_that("MPM and predict in lm", { + data(Hald, package="BAS") + hald_bic = bas.lm(Y ~ ., data=Hald, alpha=13, prior="BIC", + modelprior = uniform()) + hald_MPM = extract_MPM(hald_bic) + expect_equal(predict(hald_bic, estimator = "MPM")$fit, + predict(hald_MPM)$fit, check.attributes = FALSE) + expect_equal(predict(hald_bic, estimator = "MPM")$fit, + predict(hald_bic, estimator = "MPMold")$fit, + check.attributes = FALSE) }) From fdf331444e3214c48923e9b486dc07f028a6aa22 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Mon, 9 Nov 2020 23:29:22 -0500 Subject: [PATCH 11/25] drop steps to run pkgdown due to errors in building necessary package `textshaping` required by pkgdown --- .travis.yml | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/.travis.yml b/.travis.yml index 0fb9b17b..3340ad09 100644 --- a/.travis.yml +++ b/.travis.yml @@ -20,13 +20,13 @@ os: # - osx after_success: - Rscript -e 'covr::codecov()' - - Rscript -e 'devtools::install_github("merliseclyde/BAS")' - - Rscript -e 'pkgdown::build_site()' -deploy: - provider: pages - skip-cleanup: true - github-token: $GITHUB_PAT - keep-history: true - local-dir: docs - on: - branch: master +# - Rscript -e 'devtools::install_github("merliseclyde/BAS")' +# - Rscript -e 'pkgdown::build_site()' +#deploy: +# provider: pages +# skip-cleanup: true +# github-token: $GITHUB_PAT +# keep-history: true +# local-dir: docs +# on: +# branch: master From 1333fb0f280ce0e79850ac163a1f2d5c55a6a4ff Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Mon, 9 Nov 2020 23:37:07 -0500 Subject: [PATCH 12/25] created github actions to run pkgdown to create BAS website --- .github/.gitignore | 1 + .github/workflows/pkgdown.yaml | 48 ++++++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+) create mode 100644 .github/.gitignore create mode 100644 .github/workflows/pkgdown.yaml diff --git a/.github/.gitignore b/.github/.gitignore new file mode 100644 index 00000000..2d19fc76 --- /dev/null +++ b/.github/.gitignore @@ -0,0 +1 @@ +*.html diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml new file mode 100644 index 00000000..3c908d39 --- /dev/null +++ b/.github/workflows/pkgdown.yaml @@ -0,0 +1,48 @@ +on: + push: + branches: + - main + - master + +name: pkgdown + +jobs: + pkgdown: + runs-on: macOS-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - uses: actions/checkout@v2 + + - uses: r-lib/actions/setup-r@v1 + + - uses: r-lib/actions/setup-pandoc@v1 + + - name: Query dependencies + run: | + install.packages('remotes') + saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2) + writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version") + shell: Rscript {0} + + - name: Cache R packages + uses: actions/cache@v2 + with: + path: ${{ env.R_LIBS_USER }} + key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }} + restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1- + + - name: Install dependencies + run: | + remotes::install_deps(dependencies = TRUE) + install.packages("pkgdown", type = "binary") + shell: Rscript {0} + + - name: Install package + run: R CMD INSTALL . + + - name: Deploy package + run: | + git config --local user.email "actions@github.com" + git config --local user.name "GitHub Actions" + Rscript -e 'pkgdown::deploy_to_branch(new_process = FALSE)' From 656f3e94facc74c5349c86a2708e24ca47b9d5e4 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Mon, 9 Nov 2020 23:47:19 -0500 Subject: [PATCH 13/25] updated description to drop pkgdown from suggested packages --- DESCRIPTION | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9d0e0c00..6e11959b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BAS -Version: 1.5.6 -Date: 2020-8-24 +Version: 1.6.0 +Date: 2020-11-09 Title: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling Authors@R: c(person("Merlise", "Clyde", email="clyde@duke.edu", role=c("aut","cre", "cph"), @@ -14,10 +14,10 @@ Authors@R: c(person("Merlise", "Clyde", email="clyde@duke.edu", Depends: R (>= 3.5) Imports: - stats, graphics, + grDevices, + stats, utils, - grDevices Suggests: MASS, knitr, @@ -26,23 +26,22 @@ Suggests: rmarkdown, roxygen2, dplyr, - pkgdown, testthat, covr -Description: Package for Bayesian Variable Selection and Model Averaging - in linear models and generalized linear models using stochastic or - deterministic sampling without replacement from posterior - distributions. Prior distributions on coefficients are - from Zellner's g-prior or mixtures of g-priors +Description: Bayesian Variable Selection and Model Averaging + in linear models and generalized linear models implemented using + prior distributions on coefficients based on + Zellner's g-prior or mixtures of g-priors corresponding to the Zellner-Siow Cauchy Priors or the mixture of g-priors from Liang et al (2008) for linear models or mixtures of g-priors from Li and Clyde (2019) in generalized linear models. Other model selection criteria include AIC, BIC and Empirical Bayes - estimates of g. Sampling probabilities may be updated based on the sampled - models using sampling w/out replacement or an efficient MCMC algorithm which - samples models using a tree structure of the model space + estimates of g. Models may be sampled using Markov Chain Monte + Carlo, a deterministic sampler (for enumeration) or + sampling without replacement. Sampling probabilities may be updated based on + the sampled models using sampling w/out replacement using a tree structure of the model space as an efficient hash table. See Clyde, Ghosh and Littman (2010) for details on the sampling algorithms. Uniform priors over all models or beta-binomial prior distributions on From ff8b1465825548da9910eda097e69143304f62d1 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Mon, 9 Nov 2020 23:54:19 -0500 Subject: [PATCH 14/25] created function to extract the MPM related to issue #53 --- R/extract_models.R | 70 ++++++++++++++++++++++++++++ man/extract_MPM.Rd | 55 ++++++++++++++++++++++ tests/testthat/test-extract_models.R | 36 ++++++++++++++ 3 files changed, 161 insertions(+) create mode 100644 R/extract_models.R create mode 100644 man/extract_MPM.Rd create mode 100644 tests/testthat/test-extract_models.R diff --git a/R/extract_models.R b/R/extract_models.R new file mode 100644 index 00000000..f6757ca9 --- /dev/null +++ b/R/extract_models.R @@ -0,0 +1,70 @@ +#' Extract the Median Probability Model +#' @description Extracts the Median Probability Model from a bas object +#' @param object An object of class "bas" or "basglm" +#' @return a new object with of class "bas" or "basglm" with the Median +#' Probability Model +#' @details The Median Probability Model is the model where variables are +#' included if the marginal posterior probabilty of the coefficient being +#' zero is greater than 0.5. As this model may not have been sampled (and even +#' if it has) it is oftern faster to refit the model using bas, rather than +#' search the list of models to see where it was included. +#' @examples +#' data(Hald, package=BAS) +#' hald_bic = bas.lm(Y ~ ., data=Hald, alpha=13, prior="BIC") +#' extract_MPM(hald_bic) +#' +#' data(Pima.tr, package="MASS") +#' Pima_bas = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7, method="BAS", +#' betaprior=CCH(a=1, b=nrow(Pima.tr)/2, s=0), family=binomial(), +#' modelprior=uniform()) +#' extract_MPM(Pima_bas) +#' @family bas methods +#' @export +extract_MPM = function(object) { +# if (!(class(object) %in% c("basglm", "bas"))) { +# stop("requires an object of class 'bas' or 'basglm'") } + nvar <- object$n.vars - 1 + bestmodel <- as.numeric(object$probne0 > .5) + + if (is.null(object$call$weights)) { + object$call$weights = NULL } + + if ( !("basglm" %in% class(object))) { + # call lm + newobject <- bas.lm( + eval(object$call$formula), + data = eval(object$call$data, parent.frame()), + weights = eval(object$call$weights), + n.models = 1, + alpha = object$g, + initprobs = object$probne0, + prior = object$prior, + modelprior = object$modelprior, + update = NULL, + bestmodel = bestmodel + ) + + } +else { + glm_family = eval(object$family, parent.frame())$family + family <- get(glm_family, mode = "function", envir = parent.frame()) + newobject <- bas.glm( + eval(object$call$formula), + data = eval(object$call$data, parent.frame()), + weights = eval(object$call$weights), + family = family, + n.models = 1L, + initprobs = object$probne0, + betaprior = object$betaprior, + modelprior = object$modelprior, + update = NULL, + bestmodel = bestmodel + ) +} + newobject$probne0 = object$probne0 + mf = object$call + mf$n.models = 1 + mf$bestmodel = bestmodel + newobject$call = mf + return(newobject) +} diff --git a/man/extract_MPM.Rd b/man/extract_MPM.Rd new file mode 100644 index 00000000..d64b59e3 --- /dev/null +++ b/man/extract_MPM.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/extract_models.R +\name{extract_MPM} +\alias{extract_MPM} +\title{Extract the Median Probability Model} +\usage{ +extract_MPM(object) +} +\arguments{ +\item{object}{An object of class "bas" or "basglm"} +} +\value{ +a new object with of class "bas" or "basglm" with the Median +Probability Model +} +\description{ +Extracts the Median Probability Model from a bas object +} +\details{ +The Median Probability Model is the model where variables are +included if the marginal posterior probabilty of the coefficient being +zero is greater than 0.5. As this model may not have been sampled (and even +if it has) it is oftern faster to refit the model using bas, rather than +search the list of models to see where it was included. +} +\examples{ +data(Hald, package=BAS) +hald_bic = bas.lm(Y ~ ., data=Hald, alpha=13, prior="BIC") +extract_MPM(hald_bic) + +data(Pima.tr, package="MASS") +Pima_bas = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7, method="BAS", + betaprior=CCH(a=1, b=nrow(Pima.tr)/2, s=0), family=binomial(), + modelprior=uniform()) +extract_MPM(Pima_bas) +} +\seealso{ +Other bas methods: +\code{\link{BAS}}, +\code{\link{bas.lm}()}, +\code{\link{coef.bas}()}, +\code{\link{confint.coef.bas}()}, +\code{\link{confint.pred.bas}()}, +\code{\link{diagnostics}()}, +\code{\link{fitted.bas}()}, +\code{\link{force.heredity.bas}()}, +\code{\link{image.bas}()}, +\code{\link{plot.confint.bas}()}, +\code{\link{predict.basglm}()}, +\code{\link{predict.bas}()}, +\code{\link{summary.bas}()}, +\code{\link{update.bas}()}, +\code{\link{variable.names.pred.bas}()} +} +\concept{bas methods} diff --git a/tests/testthat/test-extract_models.R b/tests/testthat/test-extract_models.R new file mode 100644 index 00000000..ffe6dc13 --- /dev/null +++ b/tests/testthat/test-extract_models.R @@ -0,0 +1,36 @@ +test_that("extract Median Probability Model", { + + data(Hald, package="BAS") + hald_bic = bas.lm(Y ~ ., data=Hald, alpha=13, prior="BIC", + modelprior = uniform()) + hald_MPM_manual = bas.lm(Y ~ ., data=Hald, alpha=13, prior="BIC", + modelprior = uniform(), + n.models = 1L, + bestmodel = as.numeric(hald_bic$probne0 > .5) + ) + hald_MPM = extract_MPM(hald_bic) + expect_equal(hald_bic$n.vars, hald_MPM$n.vars) + expect_equal(as.numeric(hald_bic$probne0 > .5), + as.vector(which.matrix(hald_MPM$which[1], hald_MPM$n.vars))) + expect_equal(predict(hald_bic, estimator="MPM")$fit, + predict(hald_MPM)$fit, + check.attributes = FALSE) + + data(Pima.tr, package="MASS") + Pima_bas = bas.glm(type ~ ., data=Pima.tr, n.models= 2^7, method="BAS", + betaprior=CCH(a=1, b=nrow(Pima.tr)/2, s=0), + family=binomial(), + modelprior=uniform()) + Pima_MPM_man = bas.glm(type ~ ., data=Pima.tr, method="BAS", + betaprior=CCH(a=1, b=nrow(Pima.tr)/2, s=0), + family=binomial(), + modelprior=uniform(), + n.models = 1L, + bestmodel = Pima_bas$probne0 > 0.5) + + + Pima_MPM = extract_MPM(Pima_bas) + expect_equal(as.numeric(Pima_bas$probne0 > .5), + as.vector(which.matrix(Pima_MPM$which[1], Pima_MPM$n.vars))) + expect_equal(coef(Pima_MPM)$coef, coef(Pima_MPM_man)$coef) +}) From c5488eaaac6d7cc070b576288be288cde2f32c96 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Mon, 9 Nov 2020 23:58:47 -0500 Subject: [PATCH 15/25] cleanup documentation --- .Rbuildignore | 1 + NAMESPACE | 2 ++ R/BAS-package.R | 12 ++++++--- R/coefficients.R | 42 +++++++++++++++++++---------- man/BAS.Rd | 10 ++++--- man/coef.Rd | 43 ++++++++++++++++++++---------- man/confint.coef.Rd | 1 + man/confint.pred.Rd | 1 + man/diagnostics.Rd | 1 + man/force.heredity.bas.Rd | 1 + man/plot.confint.Rd | 1 + man/predict.bas.Rd | 1 + man/predict.basglm.Rd | 1 + man/summary.Rd | 1 + man/update.Rd | 1 + tests/testthat/test-coefficients.R | 23 ++++++++++++++++ 16 files changed, 107 insertions(+), 35 deletions(-) diff --git a/.Rbuildignore b/.Rbuildignore index f3398b28..c90d5c61 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -29,3 +29,4 @@ CONTRIBUTING.md ^cran-comments\.md$ ^doc$ ^Meta$ +^\.github$ diff --git a/NAMESPACE b/NAMESPACE index 3d1d8dfc..4756407e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -35,6 +35,7 @@ export(cv.summary.bas) export(diagnostics) export(eplogprob) export(eplogprob.marg) +export(extract_MPM) export(force.heredity.bas) export(g.prior) export(hyper.g) @@ -56,4 +57,5 @@ export(which.matrix) import(grDevices) import(graphics) import(stats) +import(utils) useDynLib(BAS, .registration=TRUE, .fixes="C_") diff --git a/R/BAS-package.R b/R/BAS-package.R index 3c96fffc..2d58dfe4 100644 --- a/R/BAS-package.R +++ b/R/BAS-package.R @@ -9,7 +9,7 @@ #' probabilities may be updated based on the sampled models. #' #' \tabular{ll}{ Package: \tab BAS\cr Depends: \tab R (>= 2.8)\cr License: \tab -#' GPL-3\cr URL: https://www.stat.duke.edu/~clyde\cr } +#' GPL-3\cr URL: https://www2.stat.duke.edu/~clyde\cr } #' #' Index: \preformatted{ } #' @docType package @@ -44,15 +44,19 @@ #' data("Hald") #' hald.gprior = bas.lm(Y ~ ., data=Hald, alpha=13, prior="g-prior") #' -#' # more complete demos +#' # For more examples see the demos #' -#'demo(BAS.hald) #' \dontrun{ +#' demo(BAS.hald) #' demo(BAS.USCrime) #' } -#' @import stats +#' +#' # or package vignette +#' #' @import graphics #' @import grDevices +#' @import stats +#' @import utils NULL diff --git a/R/coefficients.R b/R/coefficients.R index fdf43f73..9fc781db 100644 --- a/R/coefficients.R +++ b/R/coefficients.R @@ -48,28 +48,42 @@ #' @examples #' #' data("Hald") -#' hald.gprior = bas.lm(Y~ ., data=Hald, n.models=2^4, alpha=13, +#' hald.bas = bas.lm(Y~ ., data=Hald, n.models=2^4, alpha=13, #' prior="ZS-null", initprobs="Uniform", update=10) -#' coef.hald.gprior = coefficients(hald.gprior) -#' coef.hald.gprior -#' plot(coef.hald.gprior) -#' confint(coef.hald.gprior) +#' coef.hald.bas = coefficients(hald.bas) +#' coef.hald.bas +#' plot(coef.hald.bas) +#' confint(coef.hald.bas) #' #' #Estimation under Median Probability Model -#' coef.hald.gprior = coefficients(hald.gprior, estimator="MPM") -#' coef.hald.gprior -#' plot(coef.hald.gprior) -#' plot(confint(coef.hald.gprior)) +#' coef.hald.bas = coefficients(hald.bas, estimator="MPM") +#' coef.hald.bas +#' plot(coef.hald.bas) +#' plot(confint(coef.hald.bas)) #' +#' # Highest Probability Model +#' coef.hald.bas = coefficients(hald.bas, estimator="HPM") +#' coef.hald.bas +#' plot(coef.hald.bas) +#' confint(coef.hald.bas) #' -#' coef.hald.gprior = coefficients(hald.gprior, estimator="HPM") -#' coef.hald.gprior -#' plot(coef.hald.gprior) -#' confint(coef.hald.gprior) +#' # Best Predictive Model #' -#' # To add estimation under Best Predictive Model +#'# Note: output from coef is sorted in order of highest probability model. +#'# To extract coefficients of the highest probability model refit the model #' #' +#' hald.pred <- predict(hald.bas, estimator = "BPM") +#' BPM = as.vector(which.matrix(hald.bas$which[hald.pred$best], +#' hald.bas$n.vars)) +#' +#' hald.BPM = bas.lm(Y ~ ., data = Hald, +#' alpha=13, +#' prior="ZS-null", +#' modelprior = uniform(), +#' bestmodel = BPM, +#' n.models = 1) +#' hald.coef <- coef(hald.BPM) #' @rdname coef #' @family bas methods #' @export diff --git a/man/BAS.Rd b/man/BAS.Rd index 1b0ef23a..5291e035 100644 --- a/man/BAS.Rd +++ b/man/BAS.Rd @@ -16,7 +16,7 @@ probabilities may be updated based on the sampled models. } \details{ \tabular{ll}{ Package: \tab BAS\cr Depends: \tab R (>= 2.8)\cr License: \tab -GPL-3\cr URL: https://www.stat.duke.edu/~clyde\cr } +GPL-3\cr URL: https://www2.stat.duke.edu/~clyde\cr } Index: \preformatted{ } } @@ -24,12 +24,15 @@ Index: \preformatted{ } data("Hald") hald.gprior = bas.lm(Y ~ ., data=Hald, alpha=13, prior="g-prior") -# more complete demos +# For more examples see the demos -demo(BAS.hald) \dontrun{ +demo(BAS.hald) demo(BAS.USCrime) } + +# or package vignette + } \references{ Clyde, M. Ghosh, J. and Littman, M. (2010) Bayesian Adaptive @@ -61,6 +64,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/coef.Rd b/man/coef.Rd index 7e141985..de69d9b3 100644 --- a/man/coef.Rd +++ b/man/coef.Rd @@ -65,28 +65,42 @@ model it will be centered at the sample mean of Y. \examples{ data("Hald") -hald.gprior = bas.lm(Y~ ., data=Hald, n.models=2^4, alpha=13, +hald.bas = bas.lm(Y~ ., data=Hald, n.models=2^4, alpha=13, prior="ZS-null", initprobs="Uniform", update=10) -coef.hald.gprior = coefficients(hald.gprior) -coef.hald.gprior -plot(coef.hald.gprior) -confint(coef.hald.gprior) +coef.hald.bas = coefficients(hald.bas) +coef.hald.bas +plot(coef.hald.bas) +confint(coef.hald.bas) #Estimation under Median Probability Model -coef.hald.gprior = coefficients(hald.gprior, estimator="MPM") -coef.hald.gprior -plot(coef.hald.gprior) -plot(confint(coef.hald.gprior)) +coef.hald.bas = coefficients(hald.bas, estimator="MPM") +coef.hald.bas +plot(coef.hald.bas) +plot(confint(coef.hald.bas)) +# Highest Probability Model +coef.hald.bas = coefficients(hald.bas, estimator="HPM") +coef.hald.bas +plot(coef.hald.bas) +confint(coef.hald.bas) -coef.hald.gprior = coefficients(hald.gprior, estimator="HPM") -coef.hald.gprior -plot(coef.hald.gprior) -confint(coef.hald.gprior) +# Best Predictive Model -# To add estimation under Best Predictive Model +# Note: output from coef is sorted in order of highest probability model. +# To extract coefficients of the highest probability model refit the model +hald.pred <- predict(hald.bas, estimator = "BPM") +BPM = as.vector(which.matrix(hald.bas$which[hald.pred$best], + hald.bas$n.vars)) + +hald.BPM = bas.lm(Y ~ ., data = Hald, + alpha=13, + prior="ZS-null", + modelprior = uniform(), + bestmodel = BPM, + n.models = 1) +hald.coef <- coef(hald.BPM) } \references{ Liang, F., Paulo, R., Molina, G., Clyde, M. and Berger, J.O. @@ -103,6 +117,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/confint.coef.Rd b/man/confint.coef.Rd index f098cf9b..a5001f94 100644 --- a/man/confint.coef.Rd +++ b/man/confint.coef.Rd @@ -64,6 +64,7 @@ Other bas methods: \code{\link{coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/confint.pred.Rd b/man/confint.pred.Rd index e3440318..f298f360 100644 --- a/man/confint.pred.Rd +++ b/man/confint.pred.Rd @@ -55,6 +55,7 @@ Other bas methods: \code{\link{coef.bas}()}, \code{\link{confint.coef.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/diagnostics.Rd b/man/diagnostics.Rd index 26fe4278..6b9621c3 100644 --- a/man/diagnostics.Rd +++ b/man/diagnostics.Rd @@ -53,6 +53,7 @@ Other bas methods: \code{\link{coef.bas}()}, \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/force.heredity.bas.Rd b/man/force.heredity.bas.Rd index 9415c581..9496f688 100644 --- a/man/force.heredity.bas.Rd +++ b/man/force.heredity.bas.Rd @@ -57,6 +57,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{image.bas}()}, \code{\link{plot.confint.bas}()}, diff --git a/man/plot.confint.Rd b/man/plot.confint.Rd index 2ac605db..eb06fefd 100644 --- a/man/plot.confint.Rd +++ b/man/plot.confint.Rd @@ -52,6 +52,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/predict.bas.Rd b/man/predict.bas.Rd index 62f8aa6c..d38d492f 100644 --- a/man/predict.bas.Rd +++ b/man/predict.bas.Rd @@ -126,6 +126,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/predict.basglm.Rd b/man/predict.basglm.Rd index db958923..4c4b1bdf 100644 --- a/man/predict.basglm.Rd +++ b/man/predict.basglm.Rd @@ -100,6 +100,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/summary.Rd b/man/summary.Rd index 17935d89..ce9826d3 100644 --- a/man/summary.Rd +++ b/man/summary.Rd @@ -45,6 +45,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/man/update.Rd b/man/update.Rd index 71ba3bda..3c0b4612 100644 --- a/man/update.Rd +++ b/man/update.Rd @@ -56,6 +56,7 @@ Other bas methods: \code{\link{confint.coef.bas}()}, \code{\link{confint.pred.bas}()}, \code{\link{diagnostics}()}, +\code{\link{extract_MPM}()}, \code{\link{fitted.bas}()}, \code{\link{force.heredity.bas}()}, \code{\link{image.bas}()}, diff --git a/tests/testthat/test-coefficients.R b/tests/testthat/test-coefficients.R index 99150f1f..1ad63a3c 100644 --- a/tests/testthat/test-coefficients.R +++ b/tests/testthat/test-coefficients.R @@ -33,3 +33,26 @@ test_that("plot posterior coefficients", { expect_null(plot(crime_coef, ask=FALSE)) }) + +test_that("BPM coefficients", { + # unit test for merliseclyde/BAS#49 + data("Hald") + hald.bas = bas.lm(Y~ ., + data = Hald, + alpha = 1, + modelprior = uniform(), + prior = "BIC", n.models = 2^4) + + hald.pred <- predict(hald.bas, estimator = "BPM") + BPM = as.vector(which.matrix(hald.bas$which[hald.pred$best], + hald.bas$n.vars)) + + hald.BPM = bas.lm(Y ~ ., data = Hald, + alpha = 1, + prior = "BIC", + modelprior = uniform(), + bestmodel = BPM, n.models = 1) + hald.coef <- coef(hald.BPM) + + expect_equal(hald.coef$postmean[BPM == 1], hald.bas$mle[[hald.pred$best]]) +}) From 748076a9c9899af872ff4afabe288222b4250535 Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Tue, 10 Nov 2020 00:25:57 -0500 Subject: [PATCH 16/25] fixed typo in example for extract_MPM --- R/extract_models.R | 2 +- man/extract_MPM.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/extract_models.R b/R/extract_models.R index f6757ca9..a9b47f8e 100644 --- a/R/extract_models.R +++ b/R/extract_models.R @@ -9,7 +9,7 @@ #' if it has) it is oftern faster to refit the model using bas, rather than #' search the list of models to see where it was included. #' @examples -#' data(Hald, package=BAS) +#' data(Hald, package="BAS") #' hald_bic = bas.lm(Y ~ ., data=Hald, alpha=13, prior="BIC") #' extract_MPM(hald_bic) #' diff --git a/man/extract_MPM.Rd b/man/extract_MPM.Rd index d64b59e3..1983df96 100644 --- a/man/extract_MPM.Rd +++ b/man/extract_MPM.Rd @@ -24,7 +24,7 @@ if it has) it is oftern faster to refit the model using bas, rather than search the list of models to see where it was included. } \examples{ -data(Hald, package=BAS) +data(Hald, package="BAS") hald_bic = bas.lm(Y ~ ., data=Hald, alpha=13, prior="BIC") extract_MPM(hald_bic) From af2a488ae01240ceafd5962a99807605a552e99e Mon Sep 17 00:00:00 2001 From: merliseclyde Date: Tue, 10 Nov 2020 00:41:34 -0500 Subject: [PATCH 17/25] Fix typo in details of `extract_MPM` --- R/extract_models.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/extract_models.R b/R/extract_models.R index a9b47f8e..a1d9abd7 100644 --- a/R/extract_models.R +++ b/R/extract_models.R @@ -6,7 +6,7 @@ #' @details The Median Probability Model is the model where variables are #' included if the marginal posterior probabilty of the coefficient being #' zero is greater than 0.5. As this model may not have been sampled (and even -#' if it has) it is oftern faster to refit the model using bas, rather than +#' if it has) it is often faster to refit the model using bas, rather than #' search the list of models to see where it was included. #' @examples #' data(Hald, package="BAS") From ea0ae0da83e6fc20187849f82b84acbedb818d24 Mon Sep 17 00:00:00 2001 From: Merlise Clyde Date: Tue, 9 Nov 2021 22:05:18 -0500 Subject: [PATCH 18/25] update version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6e11959b..0b551c3c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: BAS -Version: 1.6.0 +Version: 1.6.5 Date: 2020-11-09 Title: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling Authors@R: c(person("Merlise", "Clyde", email="clyde@duke.edu", From b4054bf7da63b460868adc288b97801b2bbfa6f2 Mon Sep 17 00:00:00 2001 From: Merlise Clyde Date: Fri, 12 Nov 2021 14:15:51 -0500 Subject: [PATCH 19/25] fixed codecov.yml (removed random insert from yubikey!) --- codecov.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/codecov.yml b/codecov.yml index 7e2d5d34..ca364b8f 100644 --- a/codecov.yml +++ b/codecov.yml @@ -3,6 +3,5 @@ language: R sudo: false cache: packages after_success: - - Rscript -e 'vvenlehbnbdbvkkvlvgvkiverbvrrtnunvghvjkntdlc - covr::codecov()' + - Rscript -e 'covr::codecov()' From 62cdb814d6ae65980e223493377e8b8d93d706b4 Mon Sep 17 00:00:00 2001 From: Merlise Clyde Date: Fri, 12 Nov 2021 14:26:04 -0500 Subject: [PATCH 20/25] submit package --- CRAN-RELEASE | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 CRAN-RELEASE diff --git a/CRAN-RELEASE b/CRAN-RELEASE new file mode 100644 index 00000000..bb8ae62c --- /dev/null +++ b/CRAN-RELEASE @@ -0,0 +1,2 @@ +This package was submitted to CRAN on 2021-11-12. +Once it is accepted, delete this file and tag the release (commit 49d1620). From b7fdf202d0102a3446599a7087df349d194310b7 Mon Sep 17 00:00:00 2001 From: Merlise Clyde Date: Fri, 12 Nov 2021 14:30:47 -0500 Subject: [PATCH 21/25] fix --- codecov.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/codecov.yml b/codecov.yml index 7e2d5d34..ca364b8f 100644 --- a/codecov.yml +++ b/codecov.yml @@ -3,6 +3,5 @@ language: R sudo: false cache: packages after_success: - - Rscript -e 'vvenlehbnbdbvkkvlvgvkiverbvrrtnunvghvjkntdlc - covr::codecov()' + - Rscript -e 'covr::codecov()' From 1aa03e9aea13d171389168e4eca823a751625b16 Mon Sep 17 00:00:00 2001 From: Merlise Clyde Date: Fri, 12 Nov 2021 14:49:29 -0500 Subject: [PATCH 22/25] fixed merge issues --- .github/workflows/pkgdown.yaml | 48 +--------------------------------- DESCRIPTION | 4 ++- 2 files changed, 4 insertions(+), 48 deletions(-) diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index 856a21e5..b72cc8ed 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -1,65 +1,20 @@ -<<<<<<< HEAD -on: - push: - branches: - - main - - master -======= # Workflow derived from https://github.com/r-lib/actions/tree/master/examples # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: branches: [main, master] tags: ['*'] ->>>>>>> master name: pkgdown jobs: pkgdown: -<<<<<<< HEAD - runs-on: macOS-latest -======= runs-on: ubuntu-latest ->>>>>>> master env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - uses: actions/checkout@v2 - -<<<<<<< HEAD - - uses: r-lib/actions/setup-r@v1 - - - uses: r-lib/actions/setup-pandoc@v1 - - - name: Query dependencies - run: | - install.packages('remotes') - saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2) - writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version") - shell: Rscript {0} - - - name: Cache R packages - uses: actions/cache@v2 - with: - path: ${{ env.R_LIBS_USER }} - key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }} - restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1- - - - name: Install dependencies - run: | - remotes::install_deps(dependencies = TRUE) - install.packages("pkgdown", type = "binary") - shell: Rscript {0} - - - name: Install package - run: R CMD INSTALL . - - - name: Deploy package - run: | - git config --local user.email "actions@github.com" - git config --local user.name "GitHub Actions" -======= + - uses: r-lib/actions/setup-pandoc@v1 - uses: r-lib/actions/setup-r@v1 @@ -75,5 +30,4 @@ jobs: run: | git config --local user.name "$GITHUB_ACTOR" git config --local user.email "$GITHUB_ACTOR@users.noreply.github.com" ->>>>>>> master Rscript -e 'pkgdown::deploy_to_branch(new_process = FALSE)' diff --git a/DESCRIPTION b/DESCRIPTION index 32cb0588..7ef5f798 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -25,7 +25,9 @@ Suggests: GGally, rmarkdown, roxygen2, - dplyr + dplyr, + testthat, + codecov Description: Bayesian Variable Selection and Model Averaging in linear models and generalized linear models implemented using prior distributions on coefficients based on From 70f3172d43740d1d38d0e25870046fba941f0207 Mon Sep 17 00:00:00 2001 From: Merlise Clyde Date: Fri, 12 Nov 2021 14:50:30 -0500 Subject: [PATCH 23/25] old release --- CRAN-RELEASE | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 CRAN-RELEASE diff --git a/CRAN-RELEASE b/CRAN-RELEASE new file mode 100644 index 00000000..bb8ae62c --- /dev/null +++ b/CRAN-RELEASE @@ -0,0 +1,2 @@ +This package was submitted to CRAN on 2021-11-12. +Once it is accepted, delete this file and tag the release (commit 49d1620). From 4b08dde3022995ae94c984468835741b2010c25a Mon Sep 17 00:00:00 2001 From: Merlise Clyde Date: Fri, 12 Nov 2021 14:56:37 -0500 Subject: [PATCH 24/25] same as in master --- .github/workflows/pkgdown.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index b72cc8ed..59ae3087 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -14,7 +14,7 @@ jobs: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - uses: actions/checkout@v2 - + - uses: r-lib/actions/setup-pandoc@v1 - uses: r-lib/actions/setup-r@v1 From c3179eb7d63f258368723b7b3fadade5b89598da Mon Sep 17 00:00:00 2001 From: Merlise Clyde Date: Fri, 12 Nov 2021 15:02:16 -0500 Subject: [PATCH 25/25] fixed dependency --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7ef5f798..4d18dca9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -27,7 +27,7 @@ Suggests: roxygen2, dplyr, testthat, - codecov + covr Description: Bayesian Variable Selection and Model Averaging in linear models and generalized linear models implemented using prior distributions on coefficients based on