From 360be1f80eaf28f821ded7426b0b678338c762ff Mon Sep 17 00:00:00 2001 From: Joel Therrien Date: Tue, 23 Jul 2019 11:17:21 -0700 Subject: [PATCH] Add integrated Brier score error measure --- DESCRIPTION | 2 +- NAMESPACE | 3 + R/CR_Response.R | 45 ++++++++ R/cr_integratedBrierScore.R | 114 +++++++++++++++++++ R/createRightContinuousStepFunction.R | 11 ++ R/java_classes_directory.R | 11 +- inst/java/largeRCRF-library-1.0-SNAPSHOT.jar | Bin 115905 -> 116676 bytes man/integratedBrierScore.Rd | 69 +++++++++++ tests/testthat/test_brier_score.R | 38 +++++++ tests/testthat/test_responses_length_1.R | 41 +++++++ 10 files changed, 329 insertions(+), 5 deletions(-) create mode 100644 R/cr_integratedBrierScore.R create mode 100644 R/createRightContinuousStepFunction.R create mode 100644 man/integratedBrierScore.Rd create mode 100644 tests/testthat/test_brier_score.R diff --git a/DESCRIPTION b/DESCRIPTION index 20d8762..0973bd4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: largeRCRF Type: Package Title: Large Random Competing Risks Forests -Version: 1.0.3 +Version: 1.0.4 Authors@R: c( person("Joel", "Therrien", email = "joel_therrien@sfu.ca", role = c("aut", "cre", "cph")), person("Jiguo", "Cao", email = "jiguo_cao@sfu.ca", role = c("aut", "dgs")) diff --git a/NAMESPACE b/NAMESPACE index d26b166..c60de41 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,7 @@ # Generated by roxygen2: do not edit by hand +S3method("[",CompetingRiskResponses) +S3method("[",CompetingRiskResponsesWithCensorTimes) S3method(extractCHF,CompetingRiskFunctions) S3method(extractCHF,CompetingRiskFunctions.List) S3method(extractCIF,CompetingRiskFunctions) @@ -28,6 +30,7 @@ export(extractCHF) export(extractCIF) export(extractMortalities) export(extractSurvivorCurve) +export(integratedBrierScore) export(loadForest) export(naiveConcordance) export(saveForest) diff --git a/R/CR_Response.R b/R/CR_Response.R index 6e29d00..6d8cace 100644 --- a/R/CR_Response.R +++ b/R/CR_Response.R @@ -38,6 +38,51 @@ CR_Response <- function(delta, u, C = NULL){ } +# This function is useful is we ever want to do something like CR_Response(c(1,1,2), c(0.1,0.2,0.3))[1] +#' @export +"[.CompetingRiskResponses" <- function(object, indices){ + newList <- list( + eventIndicator = object$eventIndicator[indices], + eventTime = object$eventTime[indices] + ) + + previous.java.list <- object$javaObject + + new.java.list <- .jcall(.class_RUtils, + makeResponse(.class_List), + "produceSublist", + previous.java.list, + .jarray(as.integer(indices - 1))) + + newList$javaObject <- new.java.list + + class(newList) <- "CompetingRiskResponses" + return(newList) +} + +# This function is useful is we ever want to do something like CR_Response(c(1,1,2), c(0.1,0.2,0.3), c(2,3,4))[1] +#' @export +"[.CompetingRiskResponsesWithCensorTimes" <- function(object, indices){ + newList <- list( + eventIndicator = object$eventIndicator[indices], + eventTime = object$eventTime[indices], + censorTime = object$censorTime[indices] + ) + + previous.java.list <- object$javaObject + + new.java.list <- .jcall(.class_RUtils, + makeResponse(.class_List), + "produceSublist", + previous.java.list, + .jarray(as.integer(indices - 1))) + + newList$javaObject <- new.java.list + + class(newList) <- "CompetingRiskResponsesWithCensorTimes" + return(newList) +} + # Internal function Java_CompetingRiskResponses <- function(delta, u){ diff --git a/R/cr_integratedBrierScore.R b/R/cr_integratedBrierScore.R new file mode 100644 index 0000000..ae825dc --- /dev/null +++ b/R/cr_integratedBrierScore.R @@ -0,0 +1,114 @@ + + +#' Integrated Brier Score +#' +#' Used to calculate the Integrated Brier Score, which for the competing risks +#' setting is the integral of the squared difference between each observed +#' cumulative incidence function (CIF) for each observation and the +#' corresponding predicted CIF. If the survivor function (1 - CDF) of the +#' censoring distribution is provided, weights can be calculated to account for +#' the censoring. +#' +#' @return A numeric vector of the Integrated Brier Score for each prediction. +#' @param responses A list of responses corresponding to the provided +#' mortalities; use \code{\link{CR_Response}}. +#' @param predictions The predictions to be tested against. +#' @param event The event type for the error to be calculated on. +#' @param time \code{time} specifies the upper bound of the integral. +#' @param censoringDistribution Optional; if provided then weights are +#' calculated on the errors. There are three ways to provide it - \itemize{ +#' \item{If you have all the censor times and just want to use a simple +#' empirical estimate of the distribution, just provide a numeric vector of +#' all of the censor times and it will be automatically calculated.} \item{You +#' can directly specify the survivor function by providing a list with two +#' numeric vectors called \code{x} and \code{y}. They should be of the same +#' length and correspond to each point. It is assumed that previous to the +#' first value in \code{y} the \code{y} value is 1.0; and that the function +#' you provide is a right-continuous step function.} \item{You can provide a +#' function from \code{\link[stats]{stepfun}}. Note that this only supports +#' functions where \code{right = FALSE} (default), and that the first y value +#' (corresponding to y before the first x value) will be to set to 1.0 +#' regardless of what is specified.} +#' +#' } +#' @param parallel A logical indicating whether multiple cores should be +#' utilized when calculating the error. Available as an option because it's +#' been observed that using Java's \code{parallelStream} can be unstable on +#' some systems. Default value is \code{TRUE}; only set to \code{FALSE} if you +#' get strange errors while predicting. +#' +#' @export +#' @references Section 4.2 of Ishwaran H, Gerds TA, Kogalur UB, Moore RD, Gange +#' SJ, Lau BM (2014). “Random Survival Forests for Competing Risks.” +#' Biostatistics, 15(4), 757–773. doi:10.1093/ biostatistics/kxu010. +#' +#' @examples +#' data <- data.frame(delta=c(1,1,0,0,2,2), T=1:6, x=1:6) +#' +#' model <- train(CR_Response(delta, T) ~ x, data, ntree=100, numberOfSplits=0, mtry=1, nodeSize=1) +#' +#' newData <- data.frame(delta=c(1,0,2,1,0,2), T=1:6, x=1:6) +#' predictions <- predict(model, newData) +#' +#' scores <- integratedBrierScore(CR_Response(data$delta, data$T), predictions, 1, 6.0) +#' +integratedBrierScore <- function(responses, predictions, event, time, censoringDistribution = NULL, parallel = TRUE){ + if(length(responses$eventTime) != length(predictions)){ + stop("Length of responses and predictions must be equal.") + } + + java.censoringDistribution <- NULL + if(!is.null(censoringDistribution)){ + if(is.numeric(censoringDistribution)){ + # estimate ECDF + censoringTimes <- .jarray(censoringDistribution, "D") + java.censoringDistribution <- .jcall(.class_Utils, makeResponse(.class_RightContinuousStepFunction), "estimateOneMinusECDF", censoringTimes) + + } else if(is.list(censoringDistribution)){ + # First check that censoringDistribution fits the correct format + if(is.null(censoringDistribution$x) | is.null(censoringDistribution$y)){ + stop("If the censoringDistribution is provided as a list, it must have an x and a y item that are numeric.") + } + + if(length(censoringDistribution$x) != length(censoringDistribution$y)){ + stop("x and y in censoringDistribution must have the same length.") + } + + if(!is.numeric(censoringDistribution$x) | !is.numeric(censoringDistribution$y)){ + stop("x and y in censoringDistribution must both be numeric.") + } + + java.censoringDistribution <- createRightContinuousStepFunction(censoringDistribution$x, censoringDistribution$y, defaultY = 1.0) + + } else if("stepfun" %in% class(censoringDistribution)){ + x <- knots(censoringDistribution) + y <- censoringDistribution(x) + + java.censoringDistribution <- createRightContinuousStepFunction(x, y, defaultY = 1.0) + } + else{ + stop("Invalid censoringDistribution") + } + + # Make sure we wrap it in an Optional + java.censoringDistribution <- .object_Optional(java.censoringDistribution) + + } + else{ + java.censoringDistribution <- .object_Optional(NULL) + } + + predictions.java <- lapply(predictions, function(x){return(x$javaObject)}) + predictions.java <- convertRListToJava(predictions.java) + + errors <- .jcall(.class_CompetingRiskUtils, "[D", "calculateIBSError", + responses$javaObject, + predictions.java, + java.censoringDistribution, + as.integer(event), + time, + parallel) + + return(errors) + +} \ No newline at end of file diff --git a/R/createRightContinuousStepFunction.R b/R/createRightContinuousStepFunction.R new file mode 100644 index 0000000..5c200f6 --- /dev/null +++ b/R/createRightContinuousStepFunction.R @@ -0,0 +1,11 @@ +# Internal function +createRightContinuousStepFunction <- function(x, y, defaultY){ + x.java <- .jarray(as.numeric(x)) + y.java <- .jarray(as.numeric(y)) + + # as.numeric is explicitly required in case integers were accidently passed + # in. + newFun <- .jnew(.class_RightContinuousStepFunction, as.numeric(x), as.numeric(y), as.numeric(defaultY)) + return(newFun) + +} \ No newline at end of file diff --git a/R/java_classes_directory.R b/R/java_classes_directory.R index 00d271a..97598c3 100644 --- a/R/java_classes_directory.R +++ b/R/java_classes_directory.R @@ -14,11 +14,14 @@ # Utility Classes .class_DataUtils <- "ca/joeltherrien/randomforest/utils/DataUtils" .class_RUtils <- "ca/joeltherrien/randomforest/utils/RUtils" +.class_Utils <- "ca/joeltherrien/randomforest/utils/Utils" .class_CompetingRiskUtils <- "ca/joeltherrien/randomforest/responses/competingrisk/CompetingRiskUtils" .class_Settings <- "ca/joeltherrien/randomforest/Settings" # Misc. Classes .class_RightContinuousStepFunction <- "ca/joeltherrien/randomforest/utils/RightContinuousStepFunction" +.class_CompetingRiskResponse <- "ca/joeltherrien/randomforest/responses/competingrisk/CompetingRiskResponse" +.class_CompetingRiskResponseWithCensorTime <- "ca/joeltherrien/randomforest/responses/competingrisk/CompetingRiskResponseWithCensorTime" # TreeTrainer & its Builder .class_TreeTrainer <- "ca/joeltherrien/randomforest/tree/TreeTrainer" @@ -50,12 +53,12 @@ .class_LogRankSplitFinder <- "ca/joeltherrien/randomforest/responses/competingrisk/splitfinder/LogRankSplitFinder" .class_WeightedVarianceSplitFinder <- "ca/joeltherrien/randomforest/responses/regression/WeightedVarianceSplitFinder" -.object_Optional <- function(forest=NULL){ - if(is.null(forest)){ +.object_Optional <- function(object=NULL){ + if(is.null(object)){ return(.jcall("java/util/Optional", "Ljava/util/Optional;", "empty")) } else{ - forest <- .jcast(forest, .class_Object) - return(.jcall("java/util/Optional", "Ljava/util/Optional;", "of", forest)) + object <- .jcast(object, .class_Object) + return(.jcall("java/util/Optional", "Ljava/util/Optional;", "of", object)) } } diff --git a/inst/java/largeRCRF-library-1.0-SNAPSHOT.jar b/inst/java/largeRCRF-library-1.0-SNAPSHOT.jar index 28cdf7c109fc0782ae9f7be304332570094042c0..f7ce8af48216751df4011f4df244270d6896c1fe 100644 GIT binary patch delta 20673 zcmZ6y1yoyI&@S9i+})+PyITtscXxM};u?xOMT5IL1PE5#-QC??TU!3U?RUR*|C^Py z_GD(Bk$rY@_MS88iG=HqheJ`8dk>BM4jvxkgIe=E*E4j7>PS2D%N1-i%uAc144T;BiX0QOmuy76Jb zlrq_W8DIci+dpCmdk*jqE2=EX23TF=gSMQTOOVTyd{a@Oc zO%tf?>#V*&@Rseru`w~xXY33LXKKpWV4{*_%Lq)C5(&m&+rgJ5<&$_7PHS#BLfGpbA6zqy; za`t7tT?H(uw1_mtChwwI+8Glt(zb{>ZFx#2iZ(2CHGAb3v?UpGE-$6^EC*{wRwKsSs6e4!z z1Gpo>&_=poz^}sur#%_*9E!|c&9w*vS77K&VdHFF zbTfoTKrJl0scipYNa`ec8Hj=5?w*0UT;nV*x^d2(sMDAF0fkqRkXjRKBP#`mM^oF< zDaaoFv3IujcL_0r{tk_|1l91cB4kyQVu9Y-k}r{y!=^~U?n@O;5+DMj8E@uyOAgi+KD2=`c^k&26j|0fPN2RL;%gd%a}}zj1#F$`1gg5?(sWw$ z+9prid!Q2-sxtF`6CbPw%%tsKG}*(YcE}Y_rCvH({fftP-c=jmpzW`1;Q}5^S-)n9oi8`?<0OoSJG5+jXNSSfycc- z;v6gL8%_*fZg@WbUnl~)^{PXJRSPcYpfl~Sj4(7L2b+lsBOOK^^p=nNAHKING@92V zHdzUNj#MvGmTw`5syCZL7*Qma<*?0{*AmE#32%!r%M)AwUlr6j{Q!n<>s9fxzTdGqs$0ktACC%vFQ!yZMaB zwTAT|_sBZD&$H4K_ZdUs;X1*J&|cwaI3hMu%29X-aELhn<~emQH!s8OK9X{h(kJA13Mv^ zTMZS*=;8}8-<{PZ9(!LJ+XVGCTwP*r=aDzFBUQh&e`D0J{f_M30bNwJ!_V^$E zp^Qj%@JKE&Lpa)1Rc**skhCyrf^`+Z&rB+XZjk1=!)OMX$gB7ikpVP^I~!F*gd4Q(p)RThM+A&fNi zLDo?GBxh#WYkpC+7?LBH8XqRYSv^`mteO+^=LgI?nI!r1!N)9ir&Clfeiu~fMKwD( zuYP_Y&GuONHX@-1HI#O_cu0#9Paz;T%wiO!7#ZvaQcIU8j{EST$?#)^n2}-V$xy|1 z6fg(#%o%pEFo-lOq;iOk-CsTO#G87OPbE@D6~P1CfmwOX9Xc{%W}!$*2^PQMp>GYY zSD^F8y5sv8b&D&DAsdC&E0y?*lXcgCMFtm0^jO&80FRYXcwmKEV5*=VMbzKVJ#!1* zU-urtK8VOpBVk6vR!mkXpU6hQ)g%tZakq>9d``$uOwwN@T*zS5-TmH%MjdVM+I?`; zRNub$;~3RdafLthE{#BWE^6J$_?hGeU-otCp!9p*eUy7s=BesQqO!gl_7Xcfyf-r- zBEpOGu+C{qDD7ATcvkwL`R5J`L4>;FmT0mh0~~v+Q(n=g0e1vFimrvS>!PDCJY$yV zg|TvyO6{G7#xpE?)27jmt3%nHrdR6&{3`%h)Nq&aHBWP=x_arokNMjB!f1eJH`%nz z!0AT2bHc8y_e=x@X$T8Rm7>5|2LJ-h8p-a3iC2p%ztec8v!&Ig@6*z|nmWs`w9tAr zNghR63e%-;%i?&D?yl8_UVr~m{R6X=ly1by{)U}K>2her;b1Nd-N%HBP4C2T8$v>mrpzu^o9bThm$>TSy(dLCfg6qi0GK5Gi}u12&Y z@$ef9xG)!%pV4y7hzQepj9FK>aqKwT7K<}lWiORxKlC50VahB+kMlZ!6j{NFtH;|T zs&3_{td)a+8Cqb;eo=Etb#Vma4CgP65LX#a81J>#p6E$csND6 ziI8VA?WvtGIJyp94|$Q7*Qye_EItQ>r(76_y*@!-kfqrI_s8~>mcM@DJoz2e_e^AH zW%@nz$!m}78U2ll?~PgL8SU&1t^oNB2Yj7h)`!7PMiwoBk0G}L^lc0Me9wkJ9pOmZ zh@4fLk&>&AG`|L&P%vC9!py9UDbzxfYwjdqiW&P!J}EmnINomVt1 zuq<9~4V4zZeRn#XFc+oxp|RiTdL{hy>FD)n-NNg(3oZ7YcFda>M!!@?3ly1#OEDhm zOnL7lX5$V6j+Sr!E|73WpHcE-$?o3nPXt!k>k%9AmQHBiZAn(g-WtrRfeF)HT=EUw zAomspE=c9GKEmjYH;!k-l_QQPB9W)J%2!tSjrQ-$!K^+mXq{ATMz|A(Gkycn?=EqU zjWk?x;eK|c#sAb0%G<&_pQbUAFPnh8pBWS1;w>%A&Gf~`UH zX+8PcVsd9^i9UKv=`z*!wuZ%d{L>4>S;&(=v^|yDMx;wx7l(^h=X|azTSK2|Bj1!1LGQW1KYGw<%i)rT$ zKXfrx*Va;}&rg)Jq{^Pe%scfUDA}$xZhEUm2eJT$UkX^nay|%Ps_~XnEfY6I{jh*x zsCLutNroa}7c1U1q)e69n=AfavPovv!Ecm=I4V_YUI*j_rG`az*wYyotPCr=tU#}G z%{p!zA0!#Ic*r9=jiCXS*-)J^VZL-4E3mLSDA^^DG%|Us%#`hl;{4Pg{;*X-Iw9`% zb=9PCh=tst6n5}1uh`i5dt;-#ed4S(M+KT2{@Fvw%v9*6Jy3$HLc~H;hu?MZV94Qx zu6dfTaR!J({rUX(TPGf5s4YoXzzh>8@p)(#DsItXqecpUN87WySn9@ma8X;86K>6Y zBIdsFv|L>(evt4MT35F%BQu!EpL$SQ)n-B-65-5ON?WL-vSWm|cU;1IDb7TQV{zCN zg4zlaYD*<-+|407CcJWR_VaJ&nwxJKDpAT)d<3@P#|4E-MH-C94UD^#A{y;aW_mJd zRb3oQ!yka5_S|W~ZPr6#f~|#Sbf-r#g=VpAGczL|hiTOV^l?bJdW4Je=vGG|^?H>- zQUj|(a^>+OU@kX zoq=F~^mrw@e9B!E>QCXyo%{__rqbX!L0rMqRpc59od2!Ov}IXl7KrdS33& z*YKzaDNfJVbe$u6)M#~L#C*F&*E{Tfa`--3ygjX935vxE)er@xC;S=pbfT@CAllI= z6@Hkln;jXU;$0iEoE|VkC*!)$WcU>VC6L^2E1YaQ{NpS9r^g-I=My^D$LyoNseX0i z1_%i46W$NC!Oi54Lh(lM;I8`}t>>nwt>qifb*|0b5eAd@eQcOfl#`HI{5t+CPMF!0u!s+rRwP2DpR zLV3lfA<*GSijAgTT9s}LA>dSpnX@Qg9ITkodxnRHrL-dy~G-hII$|7^pFB~!~08zQ7At^Jd) zcQevR4@U2c|6@WWVvTOU61J#UwAXJ-qmkb;zC(voIFWZMvJ$rq6r{i;U0UIt!f}QB z@pfb6u?riOj4A{i;KL60s~K><9$)T}Nbr7DlewO|;)Z6r*#;vuSlW_eOYMf$XVn&8 zkVB{?hb3eRN*gq6JBzWTcE>bfvG<*CRqTfI@|lJ?$*GfvFQ8c_dT6bcKDGRrck41K zN~$N)sU|!W{sa6?E?!j2^h(%D<$>1D5{|$-3$joOCgZcz;DD^jMJph^@&j@}G(yv` zu72Y;Vb&o-6qNkq55$>h-7AE)Nby_(r&$3sE3W3pA7pu!ko!uv-^s4^qU}H%!&`h# zqpKRWm-~ zE%D1@b6$#A66j%N>cn+Nr?!YM9&b zRIrB=XGhS;X08Cq4k^BR@k3j2>Ve3nsUQzpoM_l9c`$qy94QUF4P8`hdkWaHB z2b3f0p{?@c#>sm9vSNK2{8_u?rr*UIlN{x#93H`RAdRIbL9+}YW{nYFit~Qw!$Vz- z$xFA(ZhSSYO7^O3Vd%#|y4Vm0vuUp3GxH33_iKOtmEug+gWRItSYd5eS;w#`(w#46 zP(z##kyylb4mkwQXKm{Ku2r@naZXcUpB(aqpS)EegukG;A0SWf(ednUbxNgf-}N+c z^mt~mN7UO5$ePE`;Eb_NkHwR?wTkRut;6=*+wY*jBD@tWc4YCNa=kUj{j{L8)BRB} zK9ic#iihWNPSGZj;#O!w>H6d zz(A5U|Ueh(;c%K*_);;rsAf zvzNNcu;jQM7UQg$lAPLiClC)Ob#r#jFo)Ln0DDLm$Dyx4=3XgTpDhjLp7G~K3SU{> z6jr*?5zGCwQlIl7$(?%ah(W%Gd-mvy9ZunSs#T3tiDCVm{lvqvYbKwS1()Q{wR%e5 zk6WU^Hy9s(>;8o?H`4ysi8)2#z&X;9Y)G7|Ww@kAVF3-7uZnfE>sZ10h>m*W(eN0I zVa{;b*Rz$)?)7sqYDde2t={)XmK(tf~3 zwSG_Olb2)PP%uAfi+R_(aiVjz+M7=}mV2!Y?7XyG_dOmFZnj|AOXV=PECrJyuiiO- zMjTs&9kF{@(JiS?DM%eK(G^G;5#UZ1wyFRtGP`RExH@VG^dWsw-R!ejZ|VM6u&zXe z_Q#yEPwR4T{r9R#ANryISqas)w8UULE7_^+*}K=#kG{n=wdY{1jK@)IE!`R9VQW|* zc&0-VC1J#wh0wEjEpE_3YOlskv+c{5j}2c=hWf;bP0Zuj8a}U+#E%|p>a}>dnFzMk zkr?n1Wq)JI-%pjyJ%II&ZC5mXuogT}Ao`Xwg;$lzTfB-2!fZx#s_i{8j)H{TnG;Zi`Xt+A_b|g$!;|*ELLb&k#^7P^IqRX5TeTPi(^OY|fcP+FJXfH1?OW**{GMGI?0`7%y7r2lRWu+EzE6AB;MxXVj8 zZx33TyH$ZE26fgM!&)~G^E0`O;9Oav$s8Z}_Ji-za z2;FtR?IOc9YU^zB6qY)eJi`JcERpNPZjpB`wGGB3k?2 zSa3xm z15cN`gNUexVVv3YXpn{~9LO+n8mI;(B{?wns;H~0G*4#(VdhA&QnjGf++?FbXhGqAZu=}BB;fIz z&Xi$)oS3n+-AfI=4MwVj6&Npa9e38))PH!OD_fv-Qx%1Da-mo_HGgiO@!KoFL`&^v z(1xI$lQ`iD@zxXC_gqpCbGj!1rPg@URY3t?` z8W$bjXbWY<_H&C}@oH&oS>yx6Q1C>O6fCWg#qno6k&JhQvb_Yck_8a7pd(>TJ$JfN z@S;qm?6X&&uXPDmRcXFdlDxHniaN2`sE{LszwA2NEnpVQ9<5dgkOH!?RJQ2 z;aI~nw@!j65LGFOFl)*wo9MFVuDjwUW@H)hxuN$Y@d$ntosivV;OX--AO~s3Tc~uq zC>RoeSQ2-<=2%2HX}HyC%hj__3{lo>p;G^)jQk04XXsf59{6d%4<;d)#~-mG+C5e` zSG}ylSc2|F2maF#tC$y{%e<6hH~GYh$rVo-pcBvJ^WWWe&jbAI6$PGtLLGe z4rQ;;mbt>siu)XiYvukHJALPXuiZ4{LbOb8A#-*eeiu{8kl~4bmfpNnuTB>`&nI7FXT5Q2Pe`lG?SC zgU;xk3w(%bAb%U4VX=3@^CMkl!G?hLm0R1y3`eD+$FDrQ{*){%q*5IE)2#`Hc`wo(KJ74#bhJ=`t(aULWDcpVsZ)$Mk#R^r zmZJ2maYK^aFU^^}w-2UrLvQ)!_r{L=ji`9EDXBbSVY8H<4OiMa%L_DV^y%_s34aRwj^6ZL=2| z$L0Wqt84NiH^IgU%LP|2%S|JvxDj#k%B6wBFssHz>0>r-YdNE8vush-(%do9&kB!1 zh_3m!t+}!=4a@oTCKK*lk*{9A&19E(Ys-Z&Ba#IXuw$P4Y^j!8)hsN9_(R%Q6rVfr zPzpL_V!bFUGWkBTY@y%8QB49lh^KB$!+0k8M@(q|3nq+=BRsgJa;?#iV6pVTvPeJW~Ils z#qpDF^0JH41Chvg5_xYU17J4llIlIyuc1`I(p}hQGC-UX_wGU zSay67w&afQ^p6e2faL*9x5ACLJPhq5SUTAi&UnEAylRImq)fmk+5PncTAse?9=PUo zYYNIclQfN;h3Qd=NqWhWd}Y2Kp#Ax{*9hgtIwUonO-C~5O7wH!ZQBg3E?48I=M})8 zvRFM-O1&WD7sdw^Sxu`7m)-QSIt2U4jTAPXJ6y8mMlo7Y7?rb~zDYJA?YctWcp^nB7_F-ty?3o6p3=|{a&$ z3}Lf?4yyoo;E`d#zjmOB5sm*;K;NJM{Gh9502*W-7tr<>213RH@IzPu0K>nk;QU_^ z|1ke`B-PackdlB?a7h-U#Q)PH^uJ6L@aU*F^k1_Ign)l+Fxu+xe=!)d{SOB?MUU_o z3fSjQLR-a4`nkOP_cO{BKj@EJG_czn@SjKYO@sd=V3e!vf1!1n{;`YR0H6fNK4JXr ziUK=tn*H^e1lX4b)iP-aBLjURF{z(D>M2P^H z{cHNq4Eq1ex=&u|KMt7Kv+f@Q<%%U~8UmA~JNf5!@X28D; z@Jaz8+JE3bnX!Pj3jwg8?xH_kU-(H)f71iQtLXlVWXgv#iOmcK)LjBV1Wn~Wk z^#!)dclx(1rhFc7|C(GW`0E|K3L*X{`D8Eik9e#v{|ACthW?B6zlHcG6Fji4Lghb~ zR{5M@57)xKA3zBN0QJ9q|MaC9{Wb2q$pY%`1>k@tdH_U_UqOJtziS6H?g!w5bSDC8 z{{kVR2>|iGJ|TfIfX;tc|Iab|KPV#T10jInKVuM!0KmdO{Z3Z^(qBW6KZa8OA^%)_ zrp!wfo5zq$M{}n>Y^#PcFR~|$_74YkCS`ZI00O(%^f;0fTztD0T z05bA_6zC7752|hlz=NvW0H_cKP5>?9U%?-TY9Qxs06YXj3aI`Un7S~4?B7)VSt$P- zJ}OW~OzpFYy_vhciJQ|G7BhPjSJw&cGgmAP>^Fxb8ft@8)f|(F&e_5f)&|~jT42|9*v_f)gY4{{_NqPeCANpX7ke>Cc%yh&#-Jf%xrHRx8oHU%iRTcCl<4av=UoC=UdppT zhB0%F-+h%QYs6aV$yk1j>W*Nd$sSnpC2tlQ`oUSKLpB`&0CZK6}y}QtNs+ z3}(_I^>?O+&&>`ih%9eos?3twReXfo*4><^6Pvb-C=>q+MpU}B((!`CLZ*&OUBGIsP8A??*2>AOJQuy((3 zB-i27n=1VESy3%H#YY={Ck-Om@m!{o+o89ppNpe%0V8Q9f|m#rGW47?vLgm}2kQYP$_9Ui*LkE47*se$k z=H>=DajHWKSUvosyc9^%4-?>Lo(yjX6Yz2hLq;7t{O18^8jX;fz$oMf(KH68eeCnv z7zv!IZrr50kA!W^V-~JUU;Nty9gDMmyjzDo-(##bMZ{AViavquFZWWU%h9@lGE~38 zWU25Xs4ZKLGgi!tzN( zTKBS|Ra3Tb$(86}(XyPx)W5#2SOoKXW{_|;MD*1Egv2fz&$bcKrqbx~J$2?rneZDQ z%!^MBPBpN3Q+369_Lu)+d1^}?vnW9maA9P9^7+{z5`96g`@HUumLa?`MY@~K-Gorc zmY97x6*V;(_3MiJK1vW6mUkr#Ps1uW zi>3$%G*x?dV^1zl6v zdAzj!g5kM>GhclZjo_H)690*K+ZuPQHbr$~a9ea_sNXhS7fEJIvDZ^2iv8=y5@*tj z4k+;*V^J>ovOw|_V1&v2ymT@>|oJJ2;73bxB?$gnshhj6+KsKI^K*=$`Kk%D&754 z`trrWC*4tCWAa-W(;4w)$Y=F?YBep8n2sOKEw}Kzu zY6}R18wq;ZfT3=Pr7wEZ!WSoD-xy@f@PEMlw|T{jSi-G@efRDd{r}m#D!%)_UMhcR zGZR)@Cksdb1>j#tJLu0A&=9zwi{z(%kn+jOlEagu^=_Docl??I-qeta#+VC$A^ooJ z^S2KSjlFEGRPU40**No!FHSwc9tInWM-0b^0~&5}{>>k@cQuzzH8%{_9Gf)Tn?LHN zct9NOW^Rw){5&l=zFwmTcWq%FPPT6S>b(2-cIC(W(v!d88)=O>-q$v4uTcE>D2(@^1TVygU4-X{Qodm-{io!yc&<|plBHuSz?W%8 z!tXmJ_4ls4+_ONuARjd#OV_tC*3bt%iN?^))G-vJ&#LUS1g5(BOq@4YFaT zfg(c89${h@;^dSN!gR=&f90c5Z?jz0}lr!If3<+FnCj$CDm^aDk z712~}zT<9Oo9Qtni}Ow9SGYTuvu<3e8RP1R^a-2Vd=AZ2y9weAl*I8OsyYK53p?&R zZHDrdnRp7tw5)axmG%~2JZTe8oS8B<*}@~)mD6y-e|mq;TP_%3ndm-ZRa0 z2MvH$cHg>nzD`UIZjGfJ9$(9cLCOoJDm&hmcqVsXmY?HMn1(>E@`0$IXKSvbA0 zAhl4i@j?|)TpY?we9s1qMIC$uddy3+r1IG#i7@19V#`7vGPOhwZNq-AvmkM_sOCx= z<2d5-l|?aR3)z!1c6@Guu+lENbPXredravw7uZFnJkfw70<8F90P%)m#AY6E`H{v) z?LEau_Xf!>Gtz*Jq#G{5qkAu!%&uP|3yvjfoLbeuD9LVClJ}Ff!X_q6c6z$}(2G&~ z&6lL*sD{C4BC{JNlHdI3U16~aaukiojTV;h>cZtVEL7oIFn^^#I0in(yH!m8^SO&^*9Oru6%&Ukgy1+k5m1}vV3h^3eJ z)2+`6^!%bjaWiAO+)MFi{J^n-lbX~l{1X>5VESsh!{CHDSEB7zY@U3qd2=9aT>D<6 zxhJ0`&bV@8qcr9yaH?X!L^3k*c(U~3NV<#;trgjjq8S6(j?7BjNAnu~>?ZHz z;VHF68~Q}{dgrInxa$_P3$;DMw95*7**{Mg3{&0(tcCfsa%F`ur}4cG^RFnnuk*`S z0+ry^OTAi1mI-vN@O`P;qHI>JF1Q2XFCsV1`?GLT`j+bor)HW=gB6d|k`&*hx~wud z-7uz4(C!lyfknRr{N;0x?`Z|{=0sbG3*yrz)qh8a?q$0>ST_sy)t~(NBY~NRqDRLtK+S47pj*u#D{H!-uZoVsj?#>bWwr@%9<$^a;{M<8BHgM=} zg;48Qe?-GeQ<&;zRI@`{uTU++$9+}&&b=-Ahg%yFUQTzDz!G_wA8Mg8%O?Zr64&<^ zR&8F^`SQq_!(5*5XmA2byRPnz`*&%UXfuM$MOr`!ChgYM7~{EnRA_>dHYM+8*ea5r zfsq#-Zzg&IKYtHz+Ki4xA=)!Y*#e>EU3P3Lxg<|DJ}57K{VKPhP$?j+_vec;k3KmS z;A+6Ph$Li?#w#`L*z}U*iTlF2Bo2iAj zgr4S>8A*wOaI5RD9L-;DW}&PpC><<-09TBy5-c9N5>g>iykW#ni%YUxODPMrl!+1?J9Q3MAm ziDxJxy{KZ&GmTuvABOy#g%b7{<0cP8=X?jgd2qc@Fo9%Pt5qXUVxmNYGx1aeDMsCX z=OqXkK7X+!3QpfRJsc4_hzJJ~Tof@)k?PlqayTyLWn_y5@qfj)`5MC>Og)K_WU%51 zC^(&YL@{;A#J+(UR*Tfp-jsxc8@5sDOZRk@{LQyb@l$t(T!2G^wD$TAAA#@XXKvY# zcX<95T%?Mou86dh25cRsI(C)CTk14qm#IxVQory}YphXI$6(8{Lha~)Rqml`c$sJU znUpo1vW&0_4uDhK5>t%ClmrY;UD5{+q+{`=bcY-1#=aP-@!tkT#k=BCF-a#<-&;2Y z?hv03r*2-IO<~sR7;XIWC+Au;N*(S42RYqys7cO$S@&5pZV2W%#SSh}FQ&QQ&9?%W zO}sStg1*_K(_MDWJyY}py$vbqOYjD&AOkLzPLZK%rxOp+F9=yXR@|U3o$=}AE{nad zRg>q@L%MWb*ZfW4sl$7<)j~h$n=7>iZ{2PKacSDoaNOxk5j=dJY2!b| z4=Ph;mj;~ZN;6sE>u0)Xg*y#;iY_SbN{KF8bJQqu^JSH6SLk5@kD)RuQjMhN*pQ-M zY`%QMGecV`^~ofDlq}t;+H6m`tl`LNAGRz6lKX;jXFa}|?OT3vEUZ;0n*2cn`aa{( zVqxl;vNNf`lo$z1&1gvYhk6dTe76_e)#`A>J&n_W3+QK=&wK8`KGne4mKJ8NY*H z;dy6F$2Nr*jqj34eteWhs$UcZ?VJAT0NH6L$-UVZ5N9DCK%FKaTtQ2e8@Bm?rab(c z{pC@uhjzO`d*D^I-yOkC!==dpX(I@KE4<7nGwV?jw3F!dt8vu%U{?d5VI0Q8vuUTo z>C2e3XwpaWGT>ESFHNly#$}$TLcxfW6EZ%)2pMOV2S5%)bB z#Ulvm>oRMyQv&ZA-o2u-5;;tjTyUH*J(Ph1^l9d3(4{fKZ4aIsten3cg~IeC;hMN1 z!DrQBYom8BbcNgQVJb1Mv>(Yyg>hG3}aLxL31(CGt}6wzBxC9vr^dF@GUn~g)6h%LSB6Qh>c6}-Fov?W(V z1Nxp^g@+1;T2q@N%@pZuExeaU;#0ZsetUr>E$}6iyPf+8fx*#(U^2ES{96k^B!_Xo zqJk7>HYjdG^|ifhJ}kQDUYh0iUP>}iAS^EhR_E8%PJ)0=9E?^Jd@ndZ7hcnV@Te_{ zgBfabOQAOC6;3U0YlPS3a-p`6HuV8S91T(&TkeKlgiS2^u=VK=JNElhP_Ha1Y35Vc z2EgRKwWPgvT)#C9Cj0zgE=|f@Sic6p)SyQjbI`3$PvN$H)Xu=QRQjz>x|hYWWBHW< zcEgNgo=_RaZdgLh7+o{*>^VIDBMRwbv1E@XA?cUkY21*gec+gLwf1yr;r(VkjkVrMF~jGn{Gg%f9;G#T z%F2+@I3zkRXzWXJ6zO=BA)+s33AEPog{IWsX_v31a|~g+cPJV=j~zgrS!$^3r$aqa#t5TB*$VgO=CG{AseG-g?x(|88#R`V z^s|k9H%{NP|J3l~$=jLN0WJ1BR0go+vxIY_D{bsqZ5o;e^LqGe2iZ$Ainlje!=)%8 ziaep*RnK8+{_eHP8SeoDI-&Gh!*iEVaM_NCrll;^6?{%2KYR&SVPra4i%n$jo^H1Fmts?gor z0m)>G=S7u)o21%z?)Qjxf%8D887&+(|WaxGXWCdM=X_8vEQpe`zN;(r)l$QQLksQ*1%-^uD z>y2?h#id||-O3Ya&Chh!M9;r5d_?brAEp}{ZJ6dv5#Lz7!BhN|bYR=zeLnf2BR3ziRpTnE<1V>;4 z<@5WlK~LgHcaTn^o$)Tb91)6?XFdvieV=R@`_Yh@GGJGP7^p zc?2B%G`QZ!e4hM{i~E%1DB?ja@|{k-F}fm+XVy$igOjbaL=?acMbuqkz@sOU(9MY`FpU~zqBv00kKl2Tqj3CsLyl6cSTa2nC@UFu%=q5E!{IqL zae_sVlwAj1+`%L#I zYtlgTK1^$oc_*e&T1)I2O+>R2RWEF89;VQE?MI~IBg-uIN9A8QB8OrhigCYO6omz3 zt;fdYwjj3ZDZmH;fbOe-zq}}#*JZd+wN-``=Ev2JzjOwszL#iw7SkP@vlgDptNnsx zFjeh_Wh9>2f)akl-$-(A))Ck_gml8+30>9#55KJV(~I3dg8Qd8{jbQbVe-zzF>ijY zTSUigHeQzyy_Vo3XFS$JI!iYbs!18eYoE&T^sfv z8_dUbj&QdVa7sJ&*FN^exZH65AZM0Ppa9dPeq2)&rxp}Fj+Gaeqe~d8AoC<~5Z8c? zEM4^1JPv&LQDCnK(Wkge|C+3tXN>cYt{A*zTS}h2K~Unh;MW>Be7>lwWTWRix}lqi zTRk#kV-h2&w>-AG^e*|Qz%NgG*@1G_Cz^ED!@st~-0-bt zx<}*}%3$HCs_)4OQ+raCtnUHOj-NDS#{3g8KasBiY=Bdctf1rIeaJo`#giG~-O0YC z@pWbZq(geNqbh=&bXFz4ModuMA)5fnehuttosLMncCp6ZS94m8OcD7zumu<*?D5HC zCca*ZOnE~#eIW8Th9E8d=|^P-7K(whulqkg74B!>mY#oO5`H9=-t?VKgv`7r7|}~` z`a)Vxcnwr77VzYB*+ z_H6ckt#MI(!OjI)>%)wbBZ^km_jXo)H}I&id{a@3?-z;N-9Fn^_EmOl`YhNYvqN=v90A zt5!XrOjy|r?UL&~kQ^M9Fs zYeGb6{?f|wZ?hh+#H*0aFv7~Tuhh|E*4J*|Uwm6I9h{?l@z!;QoW5{(;w ztRPz4XsjUfaO3A=L`fBuBuFeVInh~-#B%Ge0)}PR9Yupr0DG2)DiE@OS2$&0kiCK} zAJYFwH$aBc_)1CrST(4ZSO(#NAKqwcpauqG6$~i`rv_79HE}`SuHm%41!DwZoueQm zbG!~-h+aO|foAGSDAMex0-BzQXz)T6@yEdjY_Lp;gYcjb6%KY(!r&f%r1-^APAkA5 z3i|XQ(1eLsA`@ttLGr@>cakq~USxy=MHqcxAd%>=eaCb(hir~T4{~dWHzM~YMyy{3 zk*$KrW^zcG4&7QlcjAZ#RZr^f8vtd7_f!u3Niy<6-?e2X>{}?F=a3iX05WDI)FcCh zI&C2uQ$~BEG0ee&Mi4Y-gCIa`I;WuJ)~Rn`g0{s>z*;YY5$)0hgF{PbX$y&n993`` zvedrU^-vh85ao2vo5f2}gY{`BF(ETvjCrIXg0LFkb2w;Djs*&R6iTZ8yEG|zrX{98 zz$XDdk;6Y)X2UyE&sHhn3<6Mm@rfzPLp*B%p)E~PA(&WLw`>wuvA0)F*R)B-Pu(Gc-kfke@*K=g-vlwe;8{YxwH$NVv6YkaWZvU|l9Fnlchw!fk3l17ZW~b|!6w*{)gH8rLdVOvXX-al2zjdt zvidpHFF(&hnN2MY?G|sjB6|<){^T~ge*G$J&~Vv3Sykg#pqdM!+YV> zBfrm#ZQ*iA@;-2dUmvq;ZFm>D)4HA*rL`U=as1s!XnrlhM>Wx=I!J*>ooH+if%ctX zoKwtA{}B+JZ;%Ao*As|>#KR!SJPrblICYhhRvngBOi=>}oIax6wIq~v^*GZ{8l*W1 zzf;h$6=wai`{-E@>$c&?i$0LKOwIB|4a;;g zmp5Htp0TzeUrd`}b&E^W{)~(Mm(F9+lSDzdi(tQDmSpfv+TBV5z()uA@GzN#x!nBe=BQFA!Fj(C zggB0P@_2-^KB!=Xwkfo}B>*)8G?Ig=*<{JQh&x`>5L<0ASep5&;o%U)t3@S zdh@7cpcpu1=1Der3`WN2tm}$D6mJb@Pz~$Na(gBQIrH{E62koycne)pYK+~*1fs6V z7G)H@VIrYu?1S*6Fm`D8Em-6OxfExI*ZY#O^odD|-e-1fKbs~UM9=TER!`%V+-Tw{ zDcF19l59i5NCb3qHgTb;ol;TJy9!+-?6m7tg0O;1>0M6gWfDR!cS=RQ_B4>{Q|w5) zPLnwpS;BOk`%YLL{;*5Pp2T7ge<t&A#r zq%YUn9ypXL4iVaYfCSRY3gQwp`g9xs8tFjjux9X;f8CoV^yNjf1oc!mIGNR}J5;s$mTY{|CXn BD!%{# delta 19715 zcmZU51z23m(k+wV?(Q1goj`DR_uvrRHB5p#K?isD;O-VQxVyUtNC>riAIWuz?Q4mrNfB;yop{uH@iL1{Kmxda(thTLF#~n*r z@Lh`Ym4^I8Er3=SZ6Z68Bkp#ZON75m3SlKt!KDTU0Ro4IZU#^pqXl_#P=VW4#3 zOhsE9DNb_NV)6*+_sX360{{K*gKA}aFE+G&R63i(yz7Z=L$3w zm1z}2)y~waVey|b7@I4syy>#vxiEJnFm6pV4Eu{weXk-a>=evb@GP+Qh@$^qJ(`yT z%(l9k{ZYZFsUuJ|YVrQ;E(n*qNp=W?Si`)9w~$KFUoDcKE!x8{Jg#+Bzk#K{`QY;_ z&48}wiDfFQ$#^?Fh-Qv0ErN2#O-FVcpDv7bPVejf;7vuySb4roDk^MYUU(Jmpfg*B z^}^LFnC=)AbCuc8n7tgH8evhRM#MEkQo2jf^o@Vj;&zB9B0xzo8{~K`O5=^{Ozf9 zl^AW?y05J^i@cvR&c|j}nnuUx>VWz08B`~3sp#0@ zu4p=6Zm|{;5wYmEpP2DIA)Zi%Ze4e9NKbOqammD#qTj4$QI>2Yd0^NB4tvc^)I4MH zJ|0p!|GZC~7q`zGd4r55{H6J5LJ}7CAtzD;ITS)s{Z(kmY>8wM_8s)Vjx-QszcMj6 ztFy95+VU+-scxNXK)>$meASb?6$_w1kl@S|{=Qep$6Gw1tnZSccs#<8zvrXl#N5I> zFnFS-^89(C$s*Pf`yDjYUD5DegS9Q5qcK$*_0`E&8wE_( zglrJGm4XsSOAsb&qSvTaTBd=0-)4GlF;DOniWGr_XT^z`rt7M5didiYjzrtjD@mZ( z3WjvEugGtuJ$KmHi?z7evDJ#o4sqLB_Ak1U{+RcC-e(iTl{(-1MLr?Oc_;&Hg;VmQ zs;XGkFVPlRd-_R#&EsnuysfK7@`cj$KeMk5p{9f8TK}jXgtG?E{N)!A1%l@OX^#%s zIPAeyq^YZ7)aDiL8U!x$Ynx@bDtW-Rs3x>WGLd@td^=ophqalNS|AY&b&FZED`DEP-G{_hV3jE6 zfp>39DzZAmA)l8l7b6xhw31xe!0@hj#iGNtGWX^YdRuVSb{*S3W|$0%M=6@1MI@x5 z$+a|X9RKzwvcv@*L6x7CPy56`Zdy8n1U8Q<6U$(GX0BxMJEC@=gPLgWU3mB)kn#hg z&>fDFEDS751Iel~6a)kY0tCdrR&E7|=dYR^q8Qlc#aXZ7WpM_<$s7OZ?&IX*bD1B# zACVU|c8zkb=rsLvMLxs@H;qrheC+LBhV%IS*!}6s!;d*;MTn6NH(dY5*BRA~g{Ay#{?Q=)D$OqAP-!ux@r` zURt1vS8B4FEeWDa?7jGG>d5qfeVVNw%lQLy4Eyj;z8DWxGdccX|GZ}03P(@gJ zK#7uJa{7csbL7RKHF_xUdbanr!xR1NKdiTEx-!AnsG{o|uwj&oLfOlAmF>*58&$2J zTu$nZ@9cicmL&Z}IFR`kZ9O5jf%cw|=1fYOZ>S)jznoLloetPSl{=xeHj^wZj~%`g zC=#~(D~`=fsd5gfl#V4?Kz2=A;whV9O!lVo^x@uaV|Al-(UVyQD!O`dY|bpZaFXHe`X8btHZ7Rn!L|#rFyy%8d9Ak}`lWuvPE+|kevZxhiBJ;g6!!xiupQHNN$XMb)}uz> zOS)j7Velm>(&o1*F#W&FFmi4DhvQ}dx4i|g19l@eG^(?dPlksxu z(yi8LFH@Qr5OR~FR`bJ6OqWMM-6q*Ln-z|Mm2SVedSP z9|jA#x|KPdSMSp$O~0s8pWeCqs#Vv+W#q^s_yG?F7Vn5iN$mPH1Oyk%|EwEvUFv_= ziGv&fS`8L{6yVSD0s2A%_$!yq1)$d)6O=v6<2GRaictyxh`>=TZx^*byr-4pEjzB^ zo^=@wbiFSF)C3`7k+JXSDBtG-RQ%pci8H6BZg=mly52IOO;2k4tSGN6xhqOcEUBYg zO!+innJM@|_g6*v$N3N5UOz#FPim!g9m(D`3pUskS$pw4V<VNb1)7a`yF_l zuXBJbpMqcUPzIrtdWe1TmF;WrFY3<~>HQM6LHfaw=Mx2xju;^anSzcO#TYLG))-T< z#6zaf9rI2Tx4?*-{=ByY{HtE;ZL=TVNliHQqsz2J+t(CuAfdut z2RPV~X4_|>w=S&a>TtDvCv2t=W%50!Rb2&UPA>-dUGl9937;ZB6 z&Ua;*WTV~UG}Wx?(UCbVLG4drXS7MShvm1E=RJP zjn#$+o}V~*-||nJr+76}6`N7!>gVT@kKm!`TX?E zv|RqiEB^l7Im|H--CV^)3W8jD?Xm~4XkO!rh)L|16qYxLVTLfn#8H) z7>$8#xAIOuo+1rt7VzDDmW+||Tr#vj5=4%S>4{@ORf%)6Yb?H1qc>i$T`6BF7uORI z1+L|%maV^)5UDBW2$gNtX7G1fG-F1i`Jhd1Xlz-g#Iq#eY2)Ce9wg`j$b-wOc#{@t zN$K_LjmN4*UBcP6%PgOjIvR<+H_2#x&<*CyLq;ha|D7k=b|e*Rect-EmhwD&O+PaV zJD)`ttQd{>6rlE^mDv>B$8s`-c(y&6CN-s?@1LC29!LX7 zHV;rUx~AY2-u`?m-L>ra(K|wI z`oya6f}1?*Gth77=tu3X^=&K(F-6?Ra8%z=L&r8LUbU*@Bm`U7j(+=op9C5k*IB6_ z`&QKX&QH!?n&c)>eZz5*lUZ5sh?s$$Gm}KhQ}|bsYy$~PBonx@?{6kieVI8joYQW< zq8fCQV2Ty;2?f?X;&x}ATR_d6{ZvBLV`i7lAJ<34qcK1S`Gv1lPQFol67PaKOZOqW zsH8K1+faY7rDNU}qVc36i#{B)#(9be=DxwZ7N_o1g|Y;s0wn(1SLvi(ZM^VwhU8YEKlNCd zXrGFyrThqYGW=-k;~F#d*~#DnQKhq2$iLRx0yVniq24 zf5T+{4wtG!*oDD~%>-jJoQ!@8#1g4-z5Q5UQIxR-)N$>NdNP{T5YQnmFWh_cb$Mu+ zA1&|_1?6k~tVS{xVTiB!M$Ay-w<#~tN|xL@ugm0cGadvB(FX((Gz*3$=X>=51EaH2 zyr4@K7%Hu?S!NnVNs6T}qMxYC-W!G`SZ4GJ1jLF$jf7LP03r6Ich$v{fo}sxIK`7S za|6}H&{aoK0fGQ7!~!vBIdg2GX_Fi2P<(Q7ZV!PTbOW7Xrcq-{Q_c( zC}7LAD1Nh}SUU=0jTaOOhHyC)aS9?ym(_iTSj>gADGF24Z=la$>6%M|(5$XK@gIzU^WG2= zw(~&9!o+pj5xnwdg!GGL_Z>7fT1&OSMzoO={C=yQ39^y!BAcPT8sN7*qV9to1RhN) zXaru|ZgQh;!49YKC2rY|>M&p6y)!#EC1$2jVPOIuowqhRliK?y7ncXpsyf%M}}r#e--G}>Z`+I}3 zHcy|r)Jfb1J3jx?YvQM!i%A@jCK!}vrybp>9ZV+g!h;oC?I$dvV{(6_53_ui-e%u< zODz-Mz#=NH?()pK#>Ee=DS|1aQe!?DYdBhb=dRc^`!dxQ_NP&K38jDt!R%Q0c z>OAmc0?c<^Dz-u3nUzux9sA*ON>z$%);ob)R+{5ddxl6z@_dC?P!x|<_;7Bgs>yic z5n1)&7=mYJJzA4R60K-6@`NNH`*4o|3=2jDV{cy3RP#;KuU|fWU>%gHPVGC8Io&Kx zQqasLAjvin8k&i|IHJ?@|JVq~e9v}b&h}A)dXV^Z&W;-o_$tpWHalKGfi0wDiZV_( zf0>aD?;)Ptk87}L2Qf?vI@2Ln4lt-(DHjwt&7_hQ@Z%(rc0D^GQ4UZk2Pm8+jNc&M z;@C+hEvg_@mr3*`>G_ENfMN0*yMQ0f&WFf(1@^4q!h0XGkfB{dX++g!I)ZCJvh-BK zs=SC+ApjHxgx#b2f}zlFFC1W2Tv1AWtc#EGD{@HVSPLAuf+g}E6Xbx*b5QWGdHd#rsF+IPab!UkCn9gwzO;{h&4AP)~;bR{S--u z(~EQWD~C`&7^THMx6kEWL&r7n&jbCh5>{wK1PllWPQw3rps$HQu6q6782}iNR}O%< z#vI$?S^Qs5_GS41^f-LcmS=S1>sCjqCc-s5e1K6`0r zY5I>oC%2jqHUlW#mfikjKwR9!o(m_R*C^jMQg`M<(Rh=baA#uRDFSK<8Hh}bGKf)3 zUe!m*Umk|CD3IrxutnhW2%GLW%OIn`Y(U2#bjOvpw(>FtOK!8|ns9|k_RN2>Ke(&^ zl!SGTL<8k#Wj!)t(Mopes9ithW>cB^>+8n00~?+Y-TGRpyzVRYkHC79n>`vy#@j01 zlnn)WLoPG(sFVuExjy+0ZlgU&x10QtSqyFZEqSd9Yh8s>YrccyZLQflCj;x%8uq4H zfVWTgQmU)K(l?8%CgJgv!*Bew<70gKz3dc9sY86c;%=6HxkpDkM!6 zOinIZ?6p>yK7=nJ%mWn`bL1!J-{)OcSo3C7zxkpS>gJm0C8Fv?ZUm>chNzY$66*bS zn6s^n9<}K}O1eAOIuXXhtpsyBoqJ`W$;oGBI`siFE>l&$FUz7pftaNbla45Q;#Syu znYFD|_KV_jTK;(c<5P_CQzKRANpv69Ot}?%9QRx4NiDXFcs*bpO;ZUI9CFiatvn~C zx0Us{js0h=N3F_MavQqKYX1yt{;gM~CKrBjcvE=xjK`Bzf+-E`R^bcf9pbHpq_JYz zg|m&Ns8GhzP{kic+?)+EIc?3tg6PB4mhGhtLpA(XPi}PC!*#KE$6fIJE;bV`Hxqh> z+Om3b1B-zX>b=0B+v?hkVMT=$HkDp}Rl4%?_grTMWyHFZ3_tNUR;Yiootp{Rts~p8 zom1J&UIctTTbF*VxX}WwuwnJOaOcvShXPgC>Iag5zXj5Aumx0e2QD0Cr@4OjFTpX~ z`7n9D4EA#Y15~YFw8=5<2By@pfn&JpH;Z_j=jaBBJ@lqP;UwX-1I+ItHR;ela_8QW zAme&qwICjgM!)MpI(h9Wxofi~H%Eq>Ba*x;W{88_M?i$1-903YFWfP zx;WRego*tS!eY;9H#v8=?`$j6VwnlUL>VkZ_jX4o31{^N`6;J~$*QV}aoF?tEp;jO zjOPSctx+Gh*?X&>1W1+K`O76I`f;2wI#7Ik3~T&&fTc0wy%wkO0#s#`Dv zT(wTY*nK+jXs;av9PxZV3+?~TmoCh#`Si{rOV9z&?3GYA7w;6*;_x=*n{zUe^q;fM z)m|ToH-MkLfJPdE4jVrWvEL)wmfvz4iu8#e-qK&g8AH3F3iE$(M8TgB7EAA)>bo#F z7UxE;A|ULx+nMp*kGJjB#O-3sDC$dZ4EOFGwl60dn%WwopG5K>vWP77N*MP&T$7h| zXV2|eUnIG~1wRnBfLdEySVXSq&8AAgd05m=9=y9|q<;I;8P+L#RC?xUY1czc zH}<>|?L7$*|0h^~3VAWw!j!%Md((4M+>#`7JSO=jaV1T(n4-9FRAO2YOC98{%tu4g zC}SX2R5mh77NAG=9S5VIC7qv^D}4wwoZm?`>?Uk(0nb9?wD}f#YJyCI^lh=Xv;+i78i1q}jNLsZpj{HNUu1eqFs}Dg;+4YtA<~IU|QFT(zS6a?x`R3(%*0&X^ zpLwEB7&&7mLl0eT1KDr#i43hISNL`?a4LXL#b(4vI5=c$-+KKGPlQdAf~9wdj0_53 ze;akX+%IIBd=j{hgO+R*ib(Eb;fe{H#5-7*n&=`G{jI*+zMchGqp_S+rg#*6g>a8e z=W)p26F0N1ZPt2!OaQAk)|YB16zo90DiGb@G!*sDnD~hAX4gyKn*SFzvenJc1TTFc zb=%wcB7{+skCKf{+ay*fvrx1@;ksiST0I}O_1X+@Y>rsSB_57b!&774Q4+H9$htx@ zt=r`%?CNZOxHQ3KYd_4JC60jHmXVRc8RD=ic&mn!Q+$P+7pbHx*%;EJwDfqrR22K_ zNlG^@#M*KyM6z2(oRUQo(^YNgx0IPO(1ETuzO~^Kk|-wuGuh99zScC2(OXvZ3rC4V zX_X_|$vdmvJ&q@Nub`7_-ui?+yd4~^66P3T=9(b`iUnXuutVauw?CdwTr$5Tqy5F~ zk_VblA3aZW=EWXSYxwPjGsF5@p8&(n$~&ybDBg$H9bE%=j1?RD8v|wuM|0mIfGzWy zKa-D88SA%5WvS(Kk!qWgwO8fu=V*~sfVPhH?^Y3uoJlUmgk4Wf)T(xD3%EyoiYa;c z%Vdz0t3n28epWd8FI!axsy`5T;2L__(XXOOs7 zHYH%UG3NK=Z|@;lq0v=YKBG%_#-esw4x8-8Aj<7h*_#_1ThGjf$>Ov$uFr;Pe2^`- z3>5p2H>-1`(ZcLDv-knW`u1#VYVv87`^vl;0_~UjcDwW0+gqW%hrMMpkDo%(5fEB_ zKyS3(;I>A15@p9a9OQ}O&N1@ZO?<35uk1}A--I5sgkZ_0$>uLd=CJcnE3w9}u-qHw z%ywNBZS(%j9*!jK zZ<-C!M-3)?7Q!J;arO1pDm!NYcIqbb9frmWpr@RgStn_q9ifF863sMiEz?*!2E_ zd6E~Kdb&d2=C?znoN`qbf6=eRjQN$a`{iH7_?oJ}S|#|^yiS+z`nF!1f4!=HAV!JB zjvhci3u6o^=fa79ASTmwop_?tp+Pi<4%D9+4UJ#TPj&5{&_emGo`UIt{*ww%KVSQ^ ztYY4Xu&<_$p`kN*c50*}WAZd!%9yUot+?UvB^4c07)DtGHI5<=tv+8Yn8Ai( z$tg7bta<|5ud>ohR-eD}T8MIkq*KyX(ilr~ zym!;!^;4`qESpTxDA9PQBT9DVfvaEDPswXE_)FU4TwDYSEcILB9&=Be-56#@ zq3=J-@9`UNc<6a7MbzR=X!dObwaM3BTfpj*h!m-?ZuhVdE3)3;{)G*(lg~D(N`lV5WJ`EBiVT4hq8txA##_2vNnL?b4TV2-0q3Kh*RD@ z<{uBv3@h$5wPz;rrM@nVXruiIUifg#fg3%{7xY*iKQ_2Mr2B>iAFwO$&iNl!gEf_R zdaUpKQPLZwq=0G;ABYt%8dvlCG?VlXh}+#Xkj1VIF1&J7tof2ljyWSAwbUx7>3-FD z=v3x;CwO((A#0p|hY7p@3du8}Cr>B~YxQ;ia>f+*VTly%_e|P_PPrGER0&Fk8rq4O zAImdRS)n2hZ{1QmIB&LXRevqtb#aE6(XH=A=RToWg(8u;M!cIN`hMt^lVws9xo^Q1sSfrZNP$xhAyYf|ppUs-m^#d>?}KrG!2TohO*q60Rvm21 z+Qtl7nZR zXG6>KAd<)Az#t*KRZGPcehiWxnoEx>K0>9fLnv`kYh4h*I19F8Y|(@0w}0ld$n4I6@RuA8DPsg>3!7R4Cvg}V~hMT(o(J~Gi?DVCt%F282~siVi`6m1e|yP8W@Nhso4Pz+YLZk!P)9XqLEKcbUZVXVh8DmbzGuIi za7@zXZu6)Y!$Xp_Pe5;#_Z~w|R72Ak6#Pr9u7u}<@qAsF|WD3k&&HXtZ#1Xo-Jm035Q+vKP_iA`9 zBc{-ng3TnG=W$)vMNOoRz-qs|a^a?-Pm3ZQF_d$3d9XMvDeOfiPFUmjg`kDfUCoU- zEItS`lE7CCrof{_f8VsOSqoaF8#!MDxVm{Ay{eG6%mW4}NZIS8gsCW<3+Nw^qB(dz zr1|cpt!NSb;KyNfRZnONTE`D+fe6lGQXQV5QC6!s4Yoi@iz35hQf0*A-CgdEG4ddq z_efA(a5mj14dc>8>M2}@2-EF}UZ;EyOC_2_Zr=zp$a zkHCX>%hs0i4V;w)c|MFejkB1CWw79s(_Z7u^577QEts<^Gf~t1FiGriut*-}`LxiJ zr1~k+J4>3)h}tXGnT%{uSAyXHi5Od^HC_LXv3IZcfak`H!_=&;2eR#7&N4(^nN1 z@OGp(w{fyu}!1$i+dD>SZ?4;whj)qmfsCHgzTfuswuQLVPf z(yX+KZdIf+TSUKtFuMV27SkbMce-p<=dn`BipU{Jo{Xu}8+Pnc`4H}e{F;g^xl7x< zr^!B%s%7)`y@B~1u?@WLrf=FZY6^RVOmZYh{L&({mAMwveJehV90UK7Ki`h)BeT;b zCbrOw;tsG=kUtvfjgdFx#Z@85F;x(K!VpL$bVo_t_RD2nK%hx|r`$9cY zIpAQWKTr#7UtRaSG5;Jo>P3dPQ#!f$MF!)ZVfzl)Sr%p)IefMlxOvWC;i7z=$e*P$ zPf($S*?O;`!zUsxaotnmU<8M0^ZuVpMVczI9`M%QSWsgerJN=5gMEr&%-IB!O zE?v9e=);Y5)ZG)a%%6h6YqmGPPuprwzQk8_tJS9-DAJVf8r#J(GN2(9k`GNzX>TrR z_3JYLInx91%OVV*mjnYOW8lH!}hT^P1?q_-d%lEk+lE(ivL zzw{QMB_!mv>G;k$%GxGab32Gi@RZ8P$ZmMX4fyAt%X=go1aEjDIAcn}RpduEtLM%u zP=4l<$JDh9n3+}So3eZCFM&DMQjGo9oU-E#L@HS4J8&+p3sZLc#y&@(S~+i3zuSD$ zB%fBphM)WUWIl#B^DO7-RT~$_=AjuMigDMjRUuL5LA2Sp?m$zGC`F6123ATN564oI znd9y@!<%{M?HK{)J+=?)67}4y+$~B`y4u?f7oANB5zDk~6?vx&k7#lTB-ZSYpe<8# zpco;mUxP7EL3^Nxo&65Wl!}zQH=B#T0>sAYb%cBQhEl?gq6dTFE!*ka>utd$Hjtw% z9$U?nHVX8KXCqmlnYI6R?N1xVsaYA7;Y}8*bQy~@{Gl}V1iZpT-r;oZhN1M0KyBK+ z4GufK_aw^1d6B^$WkL*Xdd3^`1{&tAKui^({8zO4BHy~|*z@-k4m|D5GTSC=HGeG7 z%X^!}(wAF8n{xOkTdrmW4Vv;0M7Av3N^aHXI)eCz^|g6%w2QB|l*yx5st0fv$WEgM z%v?$t_4u2si1g>9n9@QOq(CZ8G6%ood%v+ON*pJm?S21g!!adsR5Kw*P(Cif4z#VM zoM6BfB-L(YaCV2JR4n_>fRN!7J?b}L-W50C7o1k4%{w3!D4?1_y~u95t<82xi07G6 zIEWCN5YE8QMPK*d-*_4U*TOP(wE7CtRK|6tLlntF+WR$(N%et(OMX_m|8*-h>DRZ4 zrROc2A1~O8prGafu+GLBUJksYfBf9x7`+B@^z~Ny&=)?l-{KMB&!+;e8=N0c;PNC9 zg8%(gpuG5;VuYy?6N-G!Mb-$RPrT$R#|jjm(~N(r5x&<0pg<%|0Q82B2!PTTAYH!> z)E~f0_VO8#UIRdkqg)+*ZWn?4w;gi>0J;V+E|2g>T%%Xf@Z7$k3JMVTOi9qAQ?kI)<-b3LLVX|BneN?tLLt4Q{FM3kf|6VD`mbvYh_{)x0lscozScb6__x zYASatv7gdYSgNik^gYiK>~_jI%xD;?HN;((icoGww)VqK46e5bhbN_?Z`s z!j4OIcu5tP8*{&)gq`*O06=XRfETW(U|f(;IXEElAD>?One8kUha3VIcN{eLmrJS! zIiw5z*<7~s`ipalpQT^=4juF#CX|{HH~g1w0rWy%e9v~~t*LT+DO`zD$mOlc(lL0^ z)Q>V%$CW>m#8oBfJkyiq(7nV5DH-p@8--kPkWv)@rRHN=_Dc*EazDPPzyvj4G^d~a z0YPE_)Si0=E{-vN_PgAe3G}HQz}_&!^>?Pm{Sm$Q15klhx&g2aLQa74*J$8q{H@P> zsnSyh#qa@mKQt;1d1vQyysnWfd&tr|b=7^ml#mJ=;Nua^~G5klE!kvubmK zWt$ka#@PuaAt-8)3?D0PHQkRYT1sb|n+0{Sz?BKMU9S$@$I_GKNUQdP zZ=cegp?zI9JINqa#baR9@65?`fDXAq_MYztaDI`DNz(u16T9Unad@=A?hqYQbjW{58suqf`$HC-6b4lam5J zy1)W{ewaQHfd?Pmw%}7GFFrSpt9u1&T=@nkYZPP(MvQ}vCl6l2?~yEu6*PG`4o=B# z&Leblh7}{~IAwN&IK+V^^z@U^ykwcVssou`Rr3<+*!^R+sfoKqQ;xd!0^ zj>HruB^lb7Q!O4_F3QaqF<><#VYHHpw&zvzdB1}^ci{|&r0T%|JwYlRK9oRus3{*| zcmX#vDk|Zd!7!5!PgG3Zlw`4$UOwNHOEHd}VGBtysOw|d;`~ofrfpEIO;j~;-_Fee zOEQNi9Na06OmB|a;?gJ7=t>0iP|r>|;x=5!rAe+)x^^pUa4UGE`+*T%9NyHNOdq6! zDUNW0TG(D8=;T;$ZMe9{x-5N*!81Kb=7DkcoT$w+67N_}5>m2D(}v^aRG*6FYT0!i zpg6E37u1KRO(oKA=#Hp3kzx*YE!;Kbk5HngI7;p);cPBp$3HV4!x=gCwiX-lpo_^n zjf~HTs5nR@grWx8Pyultx!5zZ5b>Olgl1WEdlVt<_-_p!#fXm(WW$pcMs4k*B^o)1 zOf>9ys2B<*8^zQOB{qo{UInM=f)^O%iglZ#HlC8h?Nhc6K{$4P(S z(bOQ=+7H+IUB=wVg{_c})6gKnzkmjE3=_0gE zY84aP4YsAa$l}>cFUx?Ohqe*M5@Jiv^k85qyU+*pB zu`(uP#>Le7K?Fbp<$_YB12f->|oqIM0|CG(c%+c(uee)pppRzOEuqIkg-N zjcM^>Q02cGZA^}kP+?p_ODu^Q!B%>5R1}CgMjR;f!;6p?7&jTKRd@PqbFe)pC;ikV z7iyZGmK^Wu^uzU-u;5ZUBaK(TM=>zTgUR{|Gv+C?dQ2J^6A=~2KUW(#kZxCrKTzFR z>caxCWbsr-y|YcF88#o%p1-6`RFuu($79_S-_f=mEFL)ifTu)+!74ofF3b50RTg$O!7;TOEMmd&C&GykjN;b_S=N4+wwdOx2kTd@47AVw$+18ku$ExP*bFHe=6>$VwNO0}MKeK8sF z=s0GSqIBaOVfX>0#_%{?$*jZJznGy&mo7O#oEyPYv z1s3H{R1oLyyvB|Yi-uk;&HGY@NF{;6r9U570N20-ziocU*~F3yPXgBtMY5>NJX!wZTY1dMvRyIbkKFVk_F-3L%bmI(WHL-LM^SWxGVAh(oLxuc zB5SaB2V(C#jnjRb120SS6@N^kd+yd@6$0CBVF!#E5r?mR`Niww;g>MKrde8l5W(Sj zDnrD^_1yKujg(O3n|x+DTv?z zC~(5Y3D4{VFc18gx9~b5{GkKdDB-Jp?XHPm@w;l_oMQf{?#mBK!k&hnIs{F1yb9nq zS3B?Wf%C!8^n@T9?ypJ;E*Egq)2<%iFi$OZ~@gvkX^Wq5`$$Et_E5g5{ zyt814+NIJ>aAc-CxiF!L@)?%-B%>vd4_79)j0xcQB=j;plJNQC6T+3c-6+d>qklvO zp_s4`nm7vfXA@#rnm`pyI65CxQR)L{N7W6p6f5$X+l7UeWfYf3dHEZE{JJC9DQCiT z8ne_3(ebilQk2W2juxqy|HL=lBMLWX#5cfksb~Vd*L3gAJU2IWc1>-*Nsd+WHYspn z=$AK3u)XcWN+h(3+o87*e$%Wl-U9BogsZ!>@o|qv#m>=YfuJjn`RE?TZ`zc=XfjhF z_%$>hA#zZ-`baA79MoL|mgW-YG~?h!Pee(>m|3MOfgN5Yhz`k}$p9+bPwzd78s6zQ zS7P4HvMCy<8yaP$y&KwYjOz-;>HCv900n9CvX&L3lbppO8S*R+K>Jb zzo}y1AA<%9lEvi?!;Xq^>!eV2lVC|MJUs>|pqv9tU8h!E!Ep)`iqic3)&39&;_!;s9eaK0!UIoMrO1WrXUg{WmWf$)?kNyc|fm2PDQ zgibd{_^ahwCY3k*AHpMzi_ehb@g9{vE05sKlv@9O+o!hgPEl_hL|Ctz+ZQzcRYxP6 z!hJd2`%7V1KNgLHr{}7xdfy-)sKo@mL0ZfA z#Mnt{)6v-}w55^=tu2ozhzHw5%i^Uqb4QU^22DJqUrGI z%MGMnAHx&+)yaq5U+N>x%+8=Yc=sUdCgy7gnmf9qD|_uo?kZ|K(>hy|iq*cwC$8BI zV)8CGn_dNj=#KECi?_>iKan?Y+_iGO&q>dIpHd!+14r~+PdCP<-K-Pb* z9rfCn968L&LuaH9h|;SczRJZba@tYgos+;4r}I=(OI1YGqpO zB-p(5P;Cuv=*_&SmtSX04Z5`n2VMS>sL9ZRp zK)&X?Pga!oP&T5hY8Uj`T)CJX5U#(=HZ|p}Pc{`-%#xBTN2n3ZEwavVg-Hpdel4#-qNtw3(A{115!*%HTyS#qo7WrXF(iBJ z9hCObqgkcW0J@`kHut8A9m9Tg&Go=a&x1d>!9kvRc@* z5pQOj=;_j$rF)z~~Xl zmq-qj0ZIhz(^A;9THQs)o+Odo%gqNB!z(rgd7X6LeHng}O31Lwy6bXesHWtyQn3o9 zPD4<2gMnAIIN3UFV0|8(0`q2}3KAdlrak{03NNrUzB{cLlUO(qNp>SHv902IdF<*^ zINngkl!;X31KF*h+V#dZ_@nN6XngW1rowD74A3kGl_G>_D}aPE2lpli6NEaW5z6Kg zLN-t8H%}R~NpFm6huY&9hZ|^X%;AQ9;}+*DC)XJr>N_3+pQB3^QwAR#8bIPSipB~B zM8gUDQiPS;TPPlwL|qWV_o#q_WE`YLY96Y=s_PRq9+E^8Yur9h$%Ji)BfqK9zS8e` z3k;(&rkUUfcW*dfX1>Ua*sCc)bbXCR)(6krCodMqLF>PFU&=5-USPAsn3V`S-F1sS zWJ}s_uV&d8n&e8ghJ7{S#O$3K>Vxo2)rqT@dXp5hQY1Ixzgn+~vk#N3%U`Bp#^m~5lXkg0Q z%y-w0cZ^Lz@!puZLK(}5h`uvtF@|BLOslN@q8Sa68kan3W-Gw7`jrLZXgP%eyI*!= zRi$3^xl-%AgTVLan{&(JnR~>){hGz&wbMMvlMBXfD3bf8${M9;SIK1_7oU_b8{XVg z%dCogYf!p8v7kv7Xs;-H3UVlZq7kJI5cLSbAzS|DOJ!lZlV&UMcBQ=Nq%7_1&EYQr z1c5Kv#C@-UP=9iYUQ5Cifv6A=vDE*QQ*0Qe2OvD>qd*x%;FJn#JAkTzhZDg2oZ@M? zp#jXiWJq)|0EW*Q5^z>(^|?XAI60vGB`*>}1`vB`04CXcX|SFQ5O~hPkfZ?Qo^R|@ z00sb}hS3y&+Y3Y@1E2z&fcIm;g!q2@u@5RU!4@)RT@+JWn5vni@o<#dBz%eb?)Kh? zM_WDvmm>0e<#}gd>!2KYBMpp-zs8P2!UvSmg-?lxq}0d<-UOw9aNaNBB*yhiB2hmg zMVAN36sI#5WlwF!2TQd=fLwPnSZrq*0|hPZ5qGcuL}kse)eXlJoKo@u|0u}WITOJDe1o3_kor$@>eWkb3W5UUob}JyFRg3< zBFUfZ)qkHm75xi*x!i+3X8}n6=l=lEJk#6c0tjBUz<|2)090|~L@1!PTmT~JAG-fO zjwbRiIt|ok_@7g&^v}>M&A46{00hJ}*guYcL7OpOpx}dRM9=o%{@DA1{UQY)TO)oZ zcqx^Asgz~@7eOHI3lw~m4fh!uD)Kf*6MRbUC-{Ud)4!m?l+Vz(PGkgz|B19h=cDYR zJh+tuxD@8!ow=s@(~9~}0WRIMhad969)6{NF$*rgC4QEU{PDYv;YA8Aza@K?(iMQE zi;OQ)n*#8|?Vrf~cg>C1zg$#hd4|SC(~yDM^Zy3=zmxI*?ol_-ANcbWPyFK*k@%m# zspG_ia6!8NG$+&W?0paJ+B7&I692NcDfQwQR1Sdo*|9&`f2I)mzo7MUps)1+N|gVJ z8)&r=Y!+7jS@|d5{}PPvHDKv4r57nUjgAB&t_9FMo5%=leNzUez6bws{>y}x#tZ(> zo&boV9)R&2LDV`hm4)VisGvcbbpZ5dRd_X6b#DBu0w>!^LC4j9XS#P$Oeq<-yGe); z5WN2~QUX#g0nha3&46XBB;{T#UasNL$(D_#Y`7_mNJ(#N7>pv=zXK86I zSlS%(A1M;(z7~M_tkU~4fFrRlD)4|j{|^z6&7ZO80sW@|bX?dX=5y4a8}a_JqM7=v z{%7+8{dv&$2i!t&!4abZPcDXk`Ky!t2hR80{@0)J9Lf2QQ3Q}=4d8`ftLY=YJa|Bl zkRc#g|3z@{^&fk13GV>VprT^HOWIwz0nC=r@Q)TMF8MbiXuR>SwkZp&J#LJ1)Fg>R zfkg+EWc>}TkS;|0BG^(W_?-&;M)hA`VjX|bpzJIF)^n%AViVwtz#X~=PlJE2Nltyw zNboMr%gf>H&nn2;a^AoSpMg=_-}gLp+Cq!yF35i zGA4L#Evy9gM*toI;@>Na`#+74K*_}bg6AF__bxMEfu|4!_Wz=J$N!oj1XWf1b;xou zl(qrv6bX1B|Gkc6%>PAWgNCXAn9tlTn_-=;;3?V%2Yy@ouP8!5d|UVj9rs|24q7Ps zJ0vls;3&LaeO5jN!-CdJ|0>qYzzUOp6wljGB+oW@Vt-%Mfw?Kc9{u}O^=KKG z$}@@ZFI~)t;2iikDLhytI5ks(v*F2vD9&Grd});gr-IclQZK3TL;FMYvl08PO>)zt z=QbFwQ9B+V)xsNjgp|RTXc>%;RqiEoT-BCvKMb{&PhQMQJ7&F{l@243JNdSIU*Kv@ zCKA;15g_0QFlBIBTm;Ma!##(P#DuB6U-JFvTGa-peQD6~lC4fSk{U_BqJ;Eq;5o*0 z3BEdRW{q{>4E`C;;M*j#A!CoeGB;tu$CW}_=Z|z zMpmvG|8XMUV!w-2>!GEWr9(dcVFF*e_Bk%?cdO8va+b&6TO4{#P`HzzolqGPWun^M zxUd_OrY=QMoRTli*HjjoHN(DAmaFx!-<0+ny9qCYND(^V2nYw70toeMPzeMB;_A%? z5pws%tky$CG1Da+$u3f}dR4VMLD7WY5$~1$m~lE31z8CXjLjf)Y7k5_;;mkTg5*bo z3@Uy2rs|v#uH7c_Ulf1#F$>J0m=(Z0865;JtTHcqM((Arpy)@CUD};DK5nuF`X1m5 WXyO%1%W^3$Wua72{P(f$EB^t(68K#J diff --git a/man/integratedBrierScore.Rd b/man/integratedBrierScore.Rd new file mode 100644 index 0000000..7cf2c91 --- /dev/null +++ b/man/integratedBrierScore.Rd @@ -0,0 +1,69 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cr_integratedBrierScore.R +\name{integratedBrierScore} +\alias{integratedBrierScore} +\title{Integrated Brier Score} +\usage{ +integratedBrierScore(responses, predictions, event, time, + censoringDistribution = NULL, parallel = TRUE) +} +\arguments{ +\item{responses}{A list of responses corresponding to the provided +mortalities; use \code{\link{CR_Response}}.} + +\item{predictions}{The predictions to be tested against.} + +\item{event}{The event type for the error to be calculated on.} + +\item{time}{\code{time} specifies the upper bound of the integral.} + +\item{censoringDistribution}{Optional; if provided then weights are + calculated on the errors. There are three ways to provide it - \itemize{ + \item{If you have all the censor times and just want to use a simple + empirical estimate of the distribution, just provide a numeric vector of + all of the censor times and it will be automatically calculated.} \item{You + can directly specify the survivor function by providing a list with two + numeric vectors called \code{x} and \code{y}. They should be of the same + length and correspond to each point. It is assumed that previous to the + first value in \code{y} the \code{y} value is 1.0; and that the function + you provide is a right-continuous step function.} \item{You can provide a + function from \code{\link[stats]{stepfun}}. Note that this only supports + functions where \code{right = FALSE} (default), and that the first y value + (corresponding to y before the first x value) will be to set to 1.0 + regardless of what is specified.} + + }} + +\item{parallel}{A logical indicating whether multiple cores should be +utilized when calculating the error. Available as an option because it's +been observed that using Java's \code{parallelStream} can be unstable on +some systems. Default value is \code{TRUE}; only set to \code{FALSE} if you +get strange errors while predicting.} +} +\value{ +A numeric vector of the Integrated Brier Score for each prediction. +} +\description{ +Used to calculate the Integrated Brier Score, which for the competing risks +setting is the integral of the squared difference between each observed +cumulative incidence function (CIF) for each observation and the +corresponding predicted CIF. If the survivor function (1 - CDF) of the +censoring distribution is provided, weights can be calculated to account for +the censoring. +} +\examples{ +data <- data.frame(delta=c(1,1,0,0,2,2), T=1:6, x=1:6) + +model <- train(CR_Response(delta, T) ~ x, data, ntree=100, numberOfSplits=0, mtry=1, nodeSize=1) + +newData <- data.frame(delta=c(1,0,2,1,0,2), T=1:6, x=1:6) +predictions <- predict(model, newData) + +scores <- integratedBrierScore(CR_Response(data$delta, data$T), predictions, 1, 6.0) + +} +\references{ +Section 4.2 of Ishwaran H, Gerds TA, Kogalur UB, Moore RD, Gange + SJ, Lau BM (2014). “Random Survival Forests for Competing Risks.” + Biostatistics, 15(4), 757–773. doi:10.1093/ biostatistics/kxu010. +} diff --git a/tests/testthat/test_brier_score.R b/tests/testthat/test_brier_score.R new file mode 100644 index 0000000..c8eaa99 --- /dev/null +++ b/tests/testthat/test_brier_score.R @@ -0,0 +1,38 @@ +context("Calculate integrated Brier score without error") + +# This code is more concerned that the code runs without error. The tests in the +# Java code check that the results it returns are accurate. +test_that("Can calculate Integrated Brier Score", { + + sampleData <- data.frame(x=rnorm(100)) + sampleData$T <- sample(0:4, size=100, replace = TRUE) # the censor distribution we provide needs to conform to the data or we can get NaNs + sampleData$delta <- sample(0:2, size = 100, replace = TRUE) + + testData <- sampleData[1:5,] + trainingData <- sampleData[6:100,] + + forest <- train(CR_Response(delta, T) ~ x, trainingData, ntree=50, numberOfSplits=0, mtry=1, nodeSize=5, cores=2, displayProgress=FALSE) + + predictions <- predict(forest, testData) + + scores_test <- integratedBrierScore(CR_Response(testData$delta, testData$T), predictions, event = 1, time = 4, + censoringDistribution = NULL) + # Check that we don't get a crash if we calculate the error for only one observation + scores_one <- integratedBrierScore(CR_Response(testData$delta, testData$T)[1], predictions[1], event = 1, time = 4, + censoringDistribution = NULL) + + # Make sure our error didn't somehow change + expect_equal(scores_one, scores_test[1]) + + # Provide a censoring distribution via censor times + scores_censoring1 <- integratedBrierScore(CR_Response(testData$delta, testData$T), predictions, event = 1, time = 4, + censoringDistribution = c(0,1,1,2,3,4)) + scores_censoring2 <- integratedBrierScore(CR_Response(testData$delta, testData$T), predictions, event = 1, time = 4, + censoringDistribution = list(x = 0:4, y = 1 - c(1/6, 3/6, 4/6, 5/6, 6/6))) + scores_censoring3 <- integratedBrierScore(CR_Response(testData$delta, testData$T), predictions, event = 1, time = 4, + censoringDistribution = stepfun(x=0:4, y=1 - c(0, 1/6, 3/6, 4/6, 5/6, 6/6))) + + expect_equal(scores_censoring1, scores_censoring2) + expect_equal(scores_censoring1, scores_censoring3) + +}) diff --git a/tests/testthat/test_responses_length_1.R b/tests/testthat/test_responses_length_1.R index 8bdde3d..05ac60d 100644 --- a/tests/testthat/test_responses_length_1.R +++ b/tests/testthat/test_responses_length_1.R @@ -14,4 +14,45 @@ test_that("CR_Response of length 1 - no censor times", { expect_true(T) # show Ok if we got this far +}) + +test_that("Can sub-index CR_Response - no censor times", { + + x <- CR_Response(1:5, 1:5) + + index <- 5 + + y <- x[index] + + expect_equal(y$eventTime, index) + expect_equal(y$eventIndicator, index) + + expect_equal(rJava::.jcall(y$javaObject, "I", "size"), 1) + oneJavaItem <- rJava::.jcall(y$javaObject, largeRCRF:::makeResponse(largeRCRF:::.class_Object), "get", 0L) + oneJavaItem <- rJava::.jcast(oneJavaItem, largeRCRF:::.class_CompetingRiskResponse) + delta <- rJava::.jcall(oneJavaItem, "I", "getDelta") + + expect_equal(delta, index) + +}) + +test_that("Can sub-index CR_Response - censor times", { + + x <- CR_Response(1:5, 1:5, 1:5) + + index <- 5 + + y <- x[index] + + expect_equal(y$eventTime, index) + expect_equal(y$eventIndicator, index) + expect_equal(y$censorTime, index) + + expect_equal(rJava::.jcall(y$javaObject, "I", "size"), 1) + oneJavaItem <- rJava::.jcall(y$javaObject, largeRCRF:::makeResponse(largeRCRF:::.class_Object), "get", 0L) + oneJavaItem <- rJava::.jcast(oneJavaItem, largeRCRF:::.class_CompetingRiskResponseWithCensorTime) + delta <- rJava::.jcall(oneJavaItem, "D", "getC") + + expect_equal(delta, index) + }) \ No newline at end of file