From 6b12dd621b81fd84046eb8daf561399440bdd3a6 Mon Sep 17 00:00:00 2001 From: jafiorucci Date: Fri, 22 Feb 2019 22:15:13 -0300 Subject: [PATCH] Including the fifth place R code This is the R code to reproduce the results of the GROEC method in the M4 competition, which was ranked fifth for point forecast and interval forecast. --- GROEC - M4_Code.R | 373 +++++++++++++++++++++++++++++++++++++ M4-Method-Description.docx | Bin 0 -> 13646 bytes 2 files changed, 373 insertions(+) create mode 100644 GROEC - M4_Code.R create mode 100644 M4-Method-Description.docx diff --git a/GROEC - M4_Code.R b/GROEC - M4_Code.R new file mode 100644 index 0000000..0dd89c6 --- /dev/null +++ b/GROEC - M4_Code.R @@ -0,0 +1,373 @@ + +# install.packages("devtools") +## devtools::install_github("carlanetto/M4comp2018") + +start_time <- Sys.time() + + +library(parallel) +library(snow) +library(M4comp2018) +library(forecast) +library(forecTheta) +library(dplyr) + +data(M4) + + +########################################################################################## +################ Basic functions ################################################### +########################################################################################## + +ets_forec <- function(x,h,level=0, model = "ZZZ"){ + forecast(ets(x, model=model),h=h, level=level) +} + +arima_forec <- function(x,h,level=0, model=NULL){ + if(is.null(model)) + forecast(auto.arima(x),h=h,level=level) + else + forecast(Arima(y=x,model=model),h=h,level=level) +} + + +SeasonalityTest <- function(input, ppy){ + #Used to determine whether a time series is seasonal + tcrit <- 1.645 + if (length(input)<3*ppy){ + test_seasonal <- FALSE + }else{ + xacf <- acf(input, plot = FALSE)$acf[-1, 1, 1] + clim <- tcrit/sqrt(length(input)) * sqrt(cumsum(c(1, 2 * xacf^2))) + test_seasonal <- ( abs(xacf[ppy]) > clim[ppy] ) + + if (is.na(test_seasonal)==TRUE){ test_seasonal <- FALSE } + } + + return(test_seasonal) +} + +naive2 <- function(y,h){ + input = y + fh = h + ppy <- frequency(input) ; ST <- F + if (ppy>1){ ST <- SeasonalityTest(input,ppy) } + if (ST==T){ + Dec <- decompose(input,type="multiplicative") + des_input <- input/Dec$seasonal + SIout <- head(rep(Dec$seasonal[(length(Dec$seasonal)-ppy+1):length(Dec$seasonal)], fh), fh) + }else{ + des_input <- input ; SIout <- rep(1, fh) + } + + list( mean = naive(des_input, h=fh)$mean*SIout ) +} +########################################################################################## +########################################################################################## + + +########################################################################################## +############## compute groe for OWA metric ########################################### +########################################################################################## +groe_owa <- function(y, forecFunction, n1=length(y)-10, m=5, H=length(y)-n1, p=1+floor((length(y)-n1)/m), ...){ + if(n1>=length(y)){ stop("Error in groe function: n1>=length(y)") } + if(n1<4){ stop("Error in groe function: n1<4") } + if(m<1){ stop("Error in groe function: m<1") } + if(H<1){ stop("Error in groe function: H<1") } + if(p > 1+floor((length(y)-n1)/m)){ stop("ERROR in groe function: p > 1+floor((length(y)-n1)/m)") } + if(p <= 0){ stop("ERROR in groe function: p <= 0") } + #if(!any( g==c("AE","SE","APE","sAPE") )) stop("Error in lossFunction: this g function has not been implemented.") + + n <- length(y) + fq <- frequency(y) + time_y <- time(y) + + predictionErrors <- function(i){ + ni = n1+(i-1)*m + n_pred = min(H,n-ni) + train <- window(y,start=time_y[1],end=time_y[ni]) + test <- window(y,start=time_y[ni+1],end=time_y[ni+n_pred]) + + mean_diff <- mean(abs(diff(train, lag = fq))) # used to compute MASE metric + + if(mean_diff > 0){ + predic_naive2 <- naive2( train, h=n_pred )$mean + errors_sAPE_naive2 <- errorMetric(obs=test, forec=predic_naive2, type="sAPE", statistic="N") + errors_AE_naive2 <- errorMetric(obs=test, forec=predic_naive2, type="AE", statistic="N") + errors_ASE_naive2 <- errors_AE_naive2 / mean_diff + sMAPE_naive2 <- mean(errors_sAPE_naive2) + MASE_naive2 <- mean(errors_ASE_naive2) + + if(sMAPE_naive2 > 0 && MASE_naive2 > 0){ + prediction <- forecFunction( train, h=n_pred)$mean + errors_sAPE <- errorMetric(obs=test, forec=prediction, type="sAPE", statistic="N") + errors_AE <- errorMetric(obs=test, forec=prediction, type="AE", statistic="N") + errors_ASE <- errors_AE / mean_diff + + errors <- 0.5*errors_sAPE/sMAPE_naive2 + 0.5*errors_ASE/MASE_naive2 + }else{ + errors <- 0 + } + }else{ + errors <- 0 + } + + if( i < p && n1+i*m < n) + errors <- c( errors, predictionErrors(i+1) ) + + return(errors) + } + + errors <- predictionErrors(i=1) + + errors <- replace(errors,errors==Inf,0) + + return( sum( errors ) ) +} +########################################################################################## +########################################################################################## + + +########################################################################################## +####### Run the cross-validation process for all models ################################# +########################################################################################## +runModels <- function(M_i){ + h <- M_i$h + p <- 6 + m <- trunc(h/p) + n1 <- M_i$n - M_i$h + if(n1 < 5){n1 <- 5;} + + out_dotm <- dotm(M_i$x, M_i$h, level=c(95,95)) + out_otm <- otm(M_i$x, M_i$h, level=c(95,95)) + out_ets <- ets(M_i$x) + out_arima <- auto.arima(M_i$x) + arima_order <- arimaorder(out_arima) + arima_seasonal <- c(0,0,0) + if(length(arima_order) > 3){arima_seasonal <- arima_order[4:6]; arima_order <- arima_order[1:3]} + + criterions <- matrix(NA,4,4,dimnames=list(c("DOTM","OTM", "ETS", "ARIMA"), c("AIC","AICc","BIC","GROE"))) + + criterions[1,4] <- groe_owa( y=M_i$x, forecFunction=dotm, n1=n1, m=m, H=h, p=p, level=NULL ) + criterions[2,4] <- groe_owa( y=M_i$x, forecFunction=otm, n1=n1, m=m, H=h, p=p, level=NULL ) + criterions[3,4] <- groe_owa( y=M_i$x, forecFunction=ets_forec, n1=n1, m=m, H=h, p=p, model=paste0(out_ets$components[1:3], collapse='') ) + criterions[4,4] <- groe_owa( y=M_i$x, forecFunction=arima_forec, n1=n1, m=m, H=h, p=p, model=out_arima ) + + criterions[1,1:3] = summary(out_dotm)$informationCriterions + criterions[2,1:3] = summary(out_otm)$informationCriterions + criterions[3,1:3] = unlist(out_ets[c('aic','aicc','bic')]) + criterions[4,1:3] = unlist(out_arima[c('aic','aicc','bic')]) + + out_ets = forecast(out_ets, M_i$h, level=95) + out_arima = forecast(out_arima, M_i$h, level=95) + + Rank = apply(criterions,FUN=rank,MARGIN=2) + + list(DOTM=out_dotm, OTM=out_otm, ETS=out_ets, ARIMA=out_arima, Criterions = criterions, Rank=Rank ) +} + +tryRunModels <- function(M_i){ + return(tryCatch(runModels(M_i), error=function(e) NULL)); +} +########################################################################################## +########################################################################################## + + +########################################################################################## +############ run all models for all time series from M4 ################################# +########################################################################################## + + +yearly_M4 <- Filter(function(l) l$period == "Yearly", M4) +quarterly_M4 <- Filter(function(l) l$period == "Quarterly", M4) +monthly_M4 <- Filter(function(l) l$period == "Monthly", M4) +weekly_M4 <- Filter(function(l) l$period == "Weekly", M4) +daily_M4 <- Filter(function(l) l$period == "Daily", M4) +hourly_M4 <- Filter(function(l) l$period == "Hourly", M4) + +cl<-makeCluster(detectCores()-1) + +clusterEvalQ(cl, library(M4comp2018)) +clusterEvalQ(cl, library(forecast)) +clusterEvalQ(cl, library(forecTheta)) +clusterEvalQ(cl, library(Mcomp)) + +clusterExport(cl, "runModels") +clusterExport(cl, "arima_forec") +clusterExport(cl, "groe_owa") +clusterExport(cl, "SeasonalityTest") +clusterExport(cl, "naive2") +clusterExport(cl, "ets_forec") + +clusterExport(cl, "yearly_M4") +clusterExport(cl, "quarterly_M4") +clusterExport(cl, "monthly_M4") +clusterExport(cl, "weekly_M4") +clusterExport(cl, "daily_M4") +clusterExport(cl, "hourly_M4") + + +start_time_yearly <- Sys.time() +yearly_M4.out = parLapply( cl = cl, x=yearly_M4, fun=tryRunModels ) +end_time_yearly <- Sys.time() +save.image("M4-yearly.Rdata") +rm(yearly_M4.out) + +start_time_quarterly <- Sys.time() +quarterly_M4.out = parLapply( cl = cl, x=quarterly_M4, fun=tryRunModels ) +end_time_quarterly <- Sys.time() +save.image("M4-quarterly.Rdata") +rm(quarterly_M4.out) + +start_time_monthly <- Sys.time() +monthly_M4.out = parLapply( cl = cl, x=monthly_M4, fun=tryRunModels ) +end_time_monthly <- Sys.time() +save.image("M4-monthly.Rdata") +rm(monthly_M4.out) + + +start_time_weekly <- Sys.time() +weekly_M4.out = parLapply( cl = cl, x=weekly_M4, fun=tryRunModels ) +end_time_weekly <- Sys.time() +save.image("M4-weekly.Rdata") +rm(weekly_M4.out) + +start_time_daily <- Sys.time() +daily_M4.out = parLapply( cl = cl, x=daily_M4, fun=tryRunModels ) +end_time_daily <- Sys.time() +save.image("M4-daily.Rdata") +rm(daily_M4.out) + +start_time_hourly <- Sys.time() +hourly_M4.out = parLapply( cl = cl, x=hourly_M4, fun=tryRunModels ) +end_time_hourly <- Sys.time() +save.image("M4-hourly.Rdata") +rm(hourly_M4.out) + +stopCluster(cl) + +########################################################################################## +########################################################################################## + + +########################################################################################## +############# models combination and CSV files ###################################### +########################################################################################## +modelsCombination <- function(series, outModels){ + + nn <- length(series) + if(nn != length(outModels)){ + stop("nn != length(outModels)") + return() + } + + nomeLinhas <- rep(NA, nn) + forec <- matrix(NA, nrow=nn, ncol=48) + lower95 <- matrix(NA, nrow=nn, ncol=48) + upper95 <- matrix(NA, nrow=nn, ncol=48) + + weight <- matrix(NA, nrow=nn, ncol=4) + colnames(weight) <- c('DOTM','OTM','ETS','ARIMA') + + for(i in 1:nn){ + models <- outModels[[i]] + nomeLinhas[i] <- series[[i]]$st + h <- series[[i]]$h + GROE <- models$Criterions[,'GROE'] + if( !is.null(models) && all(!is.na(GROE)) && all(GROE>0) ){ + if(!all( near(outModels[[i]]$DOTM$y , series[[i]]$x, 0.01) )){ + stop(paste0("\nERRO na série i = ",i,"\n")) + } + + Score <- 1/GROE + pesos <- Score / sum(Score) + forec[i,1:h] = pesos["DOTM"]*models[["DOTM"]]$mean + pesos["OTM"]*models[["OTM"]]$mean + pesos["ETS"]*models[["ETS"]]$mean + pesos["ARIMA"]*models[["ARIMA"]]$mean + lower95[i,1:h] = pesos["DOTM"]*models[["DOTM"]]$lower[,1] + pesos["OTM"]*models[["OTM"]]$lower[,1] + pesos["ETS"]*models[["ETS"]]$lower[,1] + pesos["ARIMA"]*models[["ARIMA"]]$lower[,1] + upper95[i,1:h] = pesos["DOTM"]*models[["DOTM"]]$upper[,1] + pesos["OTM"]*models[["OTM"]]$upper[,1] + pesos["ETS"]*models[["ETS"]]$upper[,1] + pesos["ARIMA"]*models[["ARIMA"]]$upper[,1] + weight[i,] = pesos + }else{ + cat("\nSerie = ",series[[i]]$st," , ", "i = ",i,"\n") + + out_naive <- snaive(y=series[[i]]$x, h=series[[i]]$h, level=95) + + forec[i,1:h] = out_naive$mean + lower95[i,1:h] = out_naive$lower[,1] + upper95[i,1:h] = out_naive$upper[,1] + } + } + + forec[forec < 0] <- 0 + lower95[lower95 < 0] <- 0 + upper95[upper95 < 0] <- 0 + + colnames(forec) <- colnames(lower95) <- colnames(upper95) <- c(paste0("F",1:48)) + rownames(forec) <- rownames(lower95) <- rownames(upper95) <- rownames(weight) <- nomeLinhas + + write.table(round(forec,3), file = "forec_4groec.csv", append = TRUE, quote = FALSE, sep = ",", + eol = "\n", na = "NA", dec = ".", row.names = TRUE, + col.names = FALSE) + + write.table(round(lower95,3), file = "lower95_4groec.csv", append = TRUE, quote = FALSE, sep = ",", + eol = "\n", na = "NA", dec = ".", row.names = TRUE, + col.names = FALSE) + + write.table(round(upper95,3), file = "upper95_4groec.csv", append = TRUE, quote = FALSE, sep = ",", + eol = "\n", na = "NA", dec = ".", row.names = TRUE, + col.names = FALSE) + + write.table(round(weight,3), file = "weight_4groec.csv", append = TRUE, quote = FALSE, sep = ",", + eol = "\n", na = "NA", dec = ".", row.names = TRUE, + col.names = FALSE) +} + + +load("M4-yearly.Rdata") +modelsCombination(series=yearly_M4, outModels=yearly_M4.out) +rm(yearly_M4.out) + +load("M4-quarterly.Rdata") +modelsCombination(series=quarterly_M4, outModels=quarterly_M4.out) +rm(quarterly_M4.out) + +load("M4-monthly.Rdata") +modelsCombination(series=monthly_M4, outModels=monthly_M4.out) +rm(monthly_M4.out) + +load("M4-weekly.Rdata") +modelsCombination(series=weekly_M4, outModels=weekly_M4.out) +rm(weekly_M4.out) + +load("M4-daily.Rdata") +modelsCombination(series=daily_M4, outModels=daily_M4.out) +rm(daily_M4.out) + +load("M4-hourly.Rdata") +modelsCombination(series=hourly_M4, outModels=hourly_M4.out) +rm(hourly_M4.out) +########################################################################################## +########################################################################################## + +end_time <- Sys.time() + + + + +total_time <- (end_time_yearly - start_time_yearly) + +(end_time_quarterly - start_time_quarterly) + +(end_time_monthly - start_time_monthly) + +(end_time_weekly - start_time_weekly) + +(end_time_daily - start_time_daily) + +(end_time_hourly - start_time_hourly) + +### Computer: +# DELL OPTIPLEX 7050, Windows 10 64bits, equipped with processor i7-7700 and 16GB of RAM Memory, R version 3.4.3 + +total_time/3600 ## hours +# Time difference of 61.63763 + +total_time/(3600*24) ## days +#Time difference of 2.568235 + +total_time / 100000 ## mean time per series +#Time difference of 2.218955 secs + diff --git a/M4-Method-Description.docx b/M4-Method-Description.docx new file mode 100644 index 0000000000000000000000000000000000000000..436fbc688f1aaa05cd2548e722b32f5c454ad029 GIT binary patch literal 13646 zcmeHug;yO*w)f%S?yd*-;I1K9aCdiimq3sNcM0wi+})kv?(UG_ZXb8v%)9qy=KBlY z)LN%k*Q);Q)3tl=s%<4N1rC7$fC9h(000nR^(leV6bt|$fdl~10Wd%<5j$IF6I*9J zWeA|VmrfgNKqlTn> z$EyRoX$)EeVr@vOCeI!$nHC`=a9vuItGNx`H~5CQ)nnj@I~DRTCgJktxweG|<(R>J zmycH+4hcW^5Lh8_5bz08U_sUFv-_;khUlaHAZmcNDvu382~Pw{(&+JA2T-gtW6Uh3 zNjecuAr9$3k%;JU%T^wHky+)aMn{PDsu@T^>D?DKUN{<91n1!I`N`_3Xue<8%MCi9 z$or6?rBsVaSG#;H_?>LCU+Cczrq&9h=M;CO0J*FecPjDGD4jgF1}aHTu&c{}IP{0= z6SSVw=$3vvMZb8rjnSG<_0WhTlEb0CAtA|lT&bsO7n%ogq^e9l9m$XDoh^a2_i+L* zUE$9z3U0&_gm%9gHnHC%>-7~JApbY@B#OanxO^LtebYYpH}%wWG_iJi$N0zfziRn^ zSQq~?_3+r`kKn=x!KaYx27^AF%W$Mqi)#|uZluTHf!|KR+>l1{mfG8LNt)%Cb@#6J zf5^?H2BgmL7V<$=7--C9Eu`v;uUom>G=xFRYYb|od^STAAFH&mwVQj#7Zab&7tUOg zvsDu;OfEob%37JL9W>ak+#F3H39A$pGlG8|Ueeq?aHR&G@~vP=1JT0lRH$x1okMR9 zNrw~97tZu_;YSGOAhdzqBBLR#wPe3zgeA>2J)SkwI7NDlVJRHPhDFZo^@mrzu*e{K z*lgI4@k?))G;mE_Waq^kkAVtaVvq(smXFt?NI=aybAh*w>fd!_&bUV5;f*28-~a#t z00zv>&e8bY|0s&Foso;p+iv#9Huf*20ejoa-st?_Zl$rKmOV@;fu}+30SjIXywBNz zgH7>5mDzQR&@Eu<($f=c} zPn{6Y7!sW|q?V*2pfnL8TM&cF3<;!#;SY^2c=?3DcGUDRKyl0+3G-4Pr5I_&F>qE% zF7lHe`$3y0`X&rl30I~3vrtTLkX5TtvJ?Va;5#H1)%qwI-^FQjIJ56iJo-}11xy}t3CvXwHX#{ zSyAff{`YUBpT`qD&ig>JFumVeKsgF-9H(!MZ?rq#el*M~St8wy@ z=ehfFxenQ75F({ftUYpyMD-22HlagC{qZ$DPE1v%PZo{j#+F}o)%YeI{@j4Tqf;yq zlylkx-;uaOZ|3o+*mf;tYIaf`tkhCoD)jQ*--yOi=Nt-&Xn(g7|kV>udlB{A?=le6Q5 zMlP%KqI1o4M_^URB6>-*{!KA8?xO;W6;xAh_}$0Gi>wOq-?TLK#;!)*dHyZ6XP!!!mZ1FTvtBEctJJI8;txS=(n8lbXSc9kU~=} z*{+b*KH1l+#&(gq)Ycx;_0%|TSda@o;*2Hm_36|)uBN$*e`ao^qPjv_sF-ZDB+yIC z+=E#B<+LlcRp~#WrrUgL4OrO9(Z4M>qQD^u~zXA z6hqq(2}gbL*cfN_haG}En>+=t{M?KDjTjANVL$%5@qj@Izxz$Fuygn&aC&ifJpQFv z)4aABk;c{YU8XvQ;kd}%$)tMiJdX{BC|O@-)V@p+)xETJol1Z!o1_3Yj{xszNG~yQ z*wko%_lZA&Cy9815;MN~T~ftEA}3y>YBU76Fa~qj7O0|i%nq-Aw_Y%aH2U_J8^d?iuz~qoNs}XWfO(FW33EL66yNw{b2QG*$_`4w& z9kW(iZM;GYzzPqQN)_fWKJ*=^>0^Hf3?-YKc!#W)fFsXH< zk$&0V%Ipj@&|3_{a4l|h{03((Eu6O(OB@GQabuk6C_6BlVKQm2T&LAt7)Y-2&{r3G zS?A;=8kgj&p!;1@q`5VomExjcBjy|VyU!Ip8f>0FP{z<4S3b~huYG!ErqL_oj;6Qd zf(GH)bHjiIs!^*iX1|BI4>0XLmJN;zpLXnbW{S0%m(f9RM3Cl=F#rM2LQuD}8N`G9aE)YL2$$uMM~P+HQ=vtq;CyrALd9RUY}b^bS> z78BlbajFBFq{~{?E!8Ul!igo|Qq;(&kS)hV?IG_tq={c5_Q!^&j{_)}PRH-X7$Ps& z7VyQ_+9?$>zU+)c>mpAF3|($wb6mw8P*9x$fYSz77(wXkJ2jZ^16}FYBjMXZBwFCx z47hA&sUPa!+d26!qPV)DdOx0{sEW*8TstJD)XW%vXKJ9?8-q}GpXg#2$$W1tRTJME zZGMf5O_DY+V}nnM6Q$IM_MUQNL)+dt5zQq%aGAowV{iyQ1lrUr_$Ex4N_LC>fE-co z3B|*T43iFCQa|<*qcT{rA|cuD1%C(?e!N0Yx>z})qWBtH#Lt0xn8`%3=nHeJw9EQh zxISdb$Wke2UI6XW$!e?vvE|}igE1tHD(lc>=jAk8B(Q@Rd-otU23VqGlv1|TO~1F9 z2VJ1D6AzMBK^myDe1|&!&HD(GcwIy%q_SLAPEM0Sp{n5hW)Q@(XpM_vYIfVmWSC(R zw0-c##|*$=m=4OQtBS_fR=5=g8;Pyj$8$AIs(dE4VvehBt14jL;pW5U{c{!oQ)13M82sRG%mme<*SBUOfR( zVn*a9pkZ9%B|&-(z2_^4Kc|?*+txZj)Ub>irrL>Z9EmO^{;9<{yu{M8v_I;mA+mV} zsp8L8DHZy3TU6fdAF0v&&=WAmOIwb@Q3al62Ode!=>EvU2+fIP>+JedorZa#L|YAn z(8xY+lBuX5Te@DO{!R{%V4P2Xzv_CS;QEyY_M|i}qLZJ$)|NFbG=<%8yw&3oSEWS|nRF5QMC9 z^nePfry;sO^rm5A)pvjsMFM-gGYr4hr^_1qLe#SM>M2Vb*bU4ON`x~ATZ)`|p2)ve zca9+mMC21R4qm0Ac2~Fi{ivJG5X)r0u{SNGsBvncsc1OE7n~jnBTmWp*SGIb?fweo zQRO%h*|8T6N;s&_wNKoh)rG1c)a<&q74-bPW_YH#UC4kVt_gM?OG{jKs)1tF5-e$^ z7cjr1%th9V3reaB%H{5KL`kLDu+EIQQM6Vv*{XTxthPR0V5QJz1R%w<>m-iBPKNhUsiJmuXMp7rQY>MSS?Ar~wDV z%36V(@fo~#7wgoS(GVqt67sQk*3an$ixDLhvcVFF&95hb@dtB?V4Jko0&cl7L0&bn=B(0D&vkO<1 zfy>gkw8GMlng&5nNsxqr+Fu49HbLLFl+DS0kT1Z$eD@F^Q+xbUN)EGgayu5;d{0!c zI34*%lP=0Slps_q`R2gcgygFjt+*<*zeF5j8G_({sSJBNUa*i?yl&ubFao;)56 zS40(IyWaw;Nm6_3uw>cz$Mr+^4{5fxb}&UZBC&Xfijy)>N?G`?yglNj>0o{C_yhQW zVm3vD0YmZO<9*2^2TR+og-}^@jy2PT^x!WQV8bJEGzbb}$R@;3@nQW#he6bRE(X3* z1}~?gRI8}L6~jY(>jCqt)YpN!cYd-s?VK{W{%gs>xk##D-?cnC9heTABK35F+O?J5f`{`RZOHmr9BJm(7>@q$lL{_%EI@<5z;fi#Z?*H zRk1YUdz8&(xy%wuw;GyTn%jFXTdiw*d~^8KDUjPorbc#}bDBO4|LASiFn}R$XT|$f z?$c`AVXDH#5+F6AA6rbzD^z31kGI&!v=2zCpV_|loa`;I?1xw*QLUXe*n#5q*`)v` zJg@VeNQDdcfm;J;m4HBGP-*w$-Xa_K!lw^lThc9}NVm-s`8?XHGY0(iHAwNW_d~2p zp*~Zmx)HZgZEnlBw@|nW)!a!vgjM5ntu(PM?siZVv=!fC1GKk(c6m?1rSNfo1&3L{ zDP>{-jUmvfUkFIgQ{27YlL-1DLtko?G;zXmMh1zAk1+UgyE_(dhC;qe<2DEZ+WNy1 z<{||YV+u`uDY-GJSf&j-2gRM4&@1C@a8wEJ@F&iJa*Gz*-A9ZGEYY+11cC4gx}lI= z$81`DA8~#G%j0( zn$9Y8qhtR9AJQXCaEVsOvJ7$-sSpp=U-ImirT=<2LG&eYx#^x zmt|{a>#LAaWx-g}B4kx>JPWZse55%?4E}EY8>U5({K}wxpL;m5)igEllseYAKGk#z z@#d+IEme~3ha=$^EIe{Q(0AQIA0Koj03*~HXVq@~|A6Mo-obMA2IXiq&j1IDM zLNTv*Rp?YQIXo&S3y8HAxeE+O4z0X5vaBsNDY#VUh;`!RQLXygosthh=s_PY1OKyK8;sa(77`I(cTuIhP*_v`t%pwAxk`+Mn8kC(n} zy-We$mtCq3pT`OCbqsg;pj(UgHBG;8nzd&xg7j>(`vA=b+O8in`qf?2Mz}X zn{d<|yl{x`_H4IjLhQ`agxB)~5`HFgnf_u)qva!|QIP1^ehhTKYnCeRIYvT+O!6+8 z&>)*AUO;v~`aC4Z-!a{n^XDh8+v+V5qVNY1dqv`qiNo zo!e-ZIVE&4TZx+U#LAlwIyf_qd%cpwO6piWaymmy{sH)~*2d8rMz%pP;M&{9j9qMIzL}NaPfIH~o6_mbVWf(1%j-5n1yH0`77JHm0u$9XMk)$}dyq`Y zr&^bBFb$X$?!4Z|d)W|%IDD}Fnzjj->Vn=v8f5x3{cDkcn>dxiN@(d&EV(v0Q_zgO zNuG%U7eY=IU8OMI%Wyk$ABsgP>tI?>cg`>Hs&;ZMJ z6_E-SLHll5Nk^8#`lf@Dp;tSru_GoEL^Gi} zad8QGby;P4n8l!a$_FF$gk+q2z&!tWe4d+kKbqcmQRP`aqrY}exJTct8lc86;Jv=hojc~6LQze{y zSJ09trhqxN)li=(cb!6+W&&1TT32s0>fYL@jRUQac8vq|<5?O?O-SMu{FoRMR6|U6 z^?7OEYeEav+pF>kGll8mDDg6bL`%yx z(LFQn#1Va0c7J9|G&Y^}#NK9LPZ2SOUn7T!jeZP7$+SnOb}fQIyb>T}*Ol;^yeOx{ z5a7D(}Hv8F>*;TZ3c z5XvMkcDI44jholmNs77yV^Xu6_+NbDm{kCHkTHpMt7caZU?I zye_!uOm4F7hK1UoOpXp_<9DUOJZeIqGdUrHbJfPL^wgcL3sU1(@6i9#o~parUH|Sa zoNfEIul{dX(#gcx*}~S$>5u(1O-05jl?kQIt>D#vQ+q()sL~=$+B()(ddhX{YAO}5 zXr5CMV~l<~?*>WTtCTWy6?nDwYueJSE=^B>`G~M`>cd$+aRy}j7Z1?KZa1md{hm?t zUtUWjx}>yRx)0)|&p*O$Psfgxjf^^_y-2!TXcJ#_&y}$;O(mjmXyPXft;B-#ix+Gq z{ad2eJ&N-6aLK068=ZL2(VXO$;rgN{&EEg)FEgVITD)eY#tNVyxPV!)9H6Wgyls5eSUJXe zpawi8CVAeaA0f-=;zUh;z5jAw8%lo@WY;xchKa+m!D0sN7E@Rph!QT&h2?-1zEE;y z^+1_HLtLX8M=@1gl8;zf3!=u|fq$n)V#}I;P}{_TD=OG~E^N_~%sRD{^@c5oK?EfC?}tehZDA7dVqDFek^Vw zR^nx)SqQ0u&&zX87>kwxOx|g+KUAtNkv3n?^w$ZQF{&AIi;GNUBRHX-x0W}sgPL3I z#;cX8Ji_AVa9OAwApIiB&*8LCI$-#v4=Z3k5)&ap97n;{;=C+Nd8YB{q%;2&dxlb2 zO~SCKB%NkLpVY?BcF_DZdW$Mqy^$R108}{vVozs$sg)96TdFo&4kf^D&zen|F)V1iqru`wCSq* z0GvQ~m-crRz|lPa!}5{hM)~(PoZm6C`3G)BnIfMrz>L5^I`X zwO`oRmn?eFRVt)G#&Sa&JzB7DHh-O0*>CjW7*STc|Kn+iB!Q7G9O6Pv0uu&Z>Ga0o ze3_9J1`DFO~cFrnyHtB^%vx+6ztd_zE z%E#4n8&X-e$#j;SXXNsETMhp`^>+qBy~}(1f@^O-jDLMeHxompzbuvEW5RH6WeboK zZ5xE$CpO>VV+@xf+`+B^%<1+1#s$VXv&;As;EhjDgs%{L-|IxV_{Nv2(okWP$J3Zfv3Sbv)Rrb*3}z{!t&tcr(FvXgsg z(?uRQtv%>1;va@R)#km@nHIE(26xQ-GwYvn8PXoOt3+ZP$%8Af%cv%PzOC>7KC7}` z`0p!k^11XjXOy>u=372W-qFt9>79YS{U6f#mihm067%hArzv&WE;AvtAukBNMVTC@ z^bdeyebI&o23p{~HsdMM1Kq}w(aO1t^HJ5>54a9~IEQ$0$Hu!dpqCLfB>a^on$|#&e}$KwT6ByIdVMX@dh_-J=@_ znU83qdF+Ksb0u!NtSZ@BbCEBOzWDvS0F!rwIi7X-%1UJ#k8SW>hYUVlT>F`FCIRf8 zywPdx*TD`7ZoCC)A)DNZx}%ti`QCBVlRi{%Jzy$;PitNN$pO;(bLDX$HcPCRw?XGt4kOzg4dJRE(i(cB*M=N!=l0~TY?F}fi)|94qNy+Z?t@Ps|6;RU^y}1>|?UT zSbjAPPL(UjQ@5n3m~*$lJ&YpJUX{UR&Rc@!+Ur0Lz{DA!=$DKg&uPD zxVSEleqtTBdbjR0^SkiD4dTto*&9ZkA9l#^hOB>%)4D=a88^qxUiX0+({N3s>xp8u>Y$PZ~E}oykqs`0+=xR)@7cM@=hnZ3uJ-xFGCsHhNNQKTGMM) zVwx%S;XjwPqBpT3qB!?g^v0TO>Q{d7EcjRDC~$|8eFn;xZoXIH_GD1o1VYX!^jWM@ z4^D#7Yhoo~<4)69A@HS&1vmh$HzbD(r1#B_(`1?Jt>I&r26^MLIU{0b?kAPim2M>I zI$=Tc>YsDdctXAOi>OHPV$x(JR#Q9y_QBauu)`Vz;><7SJWNwBO-$5OvNS8gkFrK9 zN*reETRUr7%(yM(Va!Dz@505i5&2@siVt@`o3^$uu&3OpA z$>hWkgxC5jEhyKVPZ6ZDP(Rqr(hw09P5pQ_{Kw;cIk{-qb}vKKI`9}p(ECp?_w7u3oWQ*mTb(XN~TxMG;)k& zj3?=any!wqAt+K3rQs`QEI556MC15K07UYUk-F`^+Wl*75;YQ(`)ys9zT@bUBt6tk z{HH!09=9j8`zB49fxh(U4r^X$8mHGKg0JrPXQV$e4P!i*gk#$z#=XM$^xC+1hCc_X zKP);mp}~4{BpptTZ?d%U$6k71e5J!JO(GqN3VwdMqsFYt1b46hs;W_u{L`OhIY7Ec zg%&)hu=r&wF?S?(rNb0*djlfXrbV`Ik@cn|3KU+P*LZp9@Cpra+oHPVi=^N2 zQif4uywZQT;78vqcs}MGSvKq!Dvu)=A(msR+oJb=#QB;qtko&?L?QG^5uC|`1s)nv zQykyQSp*Uj|CIp*E$Z!>O&?fsQ}E^H6tY{TLJJ%*A@>Awqm^Du{AVGTtze(zBMp_MAk~2r^5_X|9|^vbP(W$+sJ@w;M1! z0);TU?XQ|#w~;X3YjN=-V3uT2o>Nt;#_5$BS`c*^d zFy2XZ-oxG;Z*^zt^(7m2Fp5Ht$yMrMmxJ#QJS`CGaZGF9AEt{h>7w7PwJG0vUAk!h z@)z?=u5_#P&tmXz3mJVPkA3T%VflDs3j^G2T zHjZPZ-PvkmjP?~KPdtYu3tWSuVx)i4l7JWK&=>=_FM50&S=nQdrrNn!Di><4zvI~I zv#`o7j9B}*Q)J~fOLBeUcCFyWK5~ik>qe4qvqh6%m@J-QQR0MEPoI?mCgxf;8W}s$ znkk$#TY%$UP?i?qUr+*EMAjg(YzI>TGa0R^!p&rMK&(`JguLWo)E0$?_?1@(EXZ0 zc0hTr3tOne$W~Z&Yg?u?%Q#N3uQJF`gOqka6xU7r9@&e2o!9B1)IZ)my%m=Ft3r8<)JC^fCtvmD|Bst$+>?pVuR_NzY3vM3o1*LvfQVHi$4CE*md{Ouq_Txr=2-rnB4I;t#}$ zBQ=Vd=xQ{l!5eMv)*I=q&-yT#?F~{3f7OwRDtl#)n)nq{Bu}yS;xFaKE z09_=z27%$~OdDwV$FoTCe=@?ru>K!YGr(DV4T7G76GyggC>79K z01a;Z2U|uFkGQqmb>pAhk2+%@X7&+=3-Jd1M8uk=OOlSC&(+-QLan4?<~TV$(pckz)FA8=%Ds9|jd^gK`N00kUKN|uj9`0~h$ zj-dQ6??>?#z6q8IRoBpF(yP_@TS|oV!NH=qJ^?;wy`Rg4qDx11{;;Y#8UD#0c&@Ep zfBRTC5O;pqh{C?fhS&*pM5#C6=m$rJZ%KA9T>B0Aeg3$Zvk87|=k!cI$hGB%5QY}b z@vlI2pU+(|S&Zev*$JDoV%_ZMS-apBgtzB*DldH4Sr7P~HG2yR&FxWA=x5+% z2Z|(Fog#8L{3{O3ERR)_J%=1Pfo+?d>j5H#(D2Wz?`}6F1Zme!gfTcz&Ns}z}I|eYGj$^E$aYIk($d z1aYBr*yox}8a?v!_GeB-QP3fTB`me*zEL{~nXaMg2s86vrp^$3mFHg`HZv+5uk`i5 zxl@}@0Zp3WY^RKq(o1W#_{7DHX24_F3&3Oz|Kuj3mQ3+;>is#4tfX;SDAE(`Y6B6! zJ4an>oodJD+`%9U$>zTs`8}C#b&kkNcplR6@dI6g3UmE9?nCcJ_KlR?x8|U%e6iFN zX*OhGqf>9J<0a(R`{<$?G?*s`)&nQ?yMtqk%}XiN0pe$6+(Oag$nqdhLIJ|7^SwHj z1k`HY{-jXX0f-HwWgMSktW8(nYTdC9a9bz*V9B8-2$p0WiPR*PN(})v&%f$zk4sy8 zc1m$`?ZBCwJHo+2fe=&_>F0g$h?ZXl;|Rq`R-xLt4p541Z)PN0_pr#Qo*DR9r=acB zck8VPoT1YUkk-!91R1*+yPC@u&4LX%xC8Uh1}Ro1gL3xRchJiyR(Bg1K6WFy%+c(Q zEZ;J#Jx_jXpJZWNSSiIgO`X{7d2WYY*)_KH1h*sbSIf6xKUc z?1B2HP51w)9{!!u@71G!(t3S6UGs0Xq`$*|&z$}Ve~-ZB70K}8M`Fp=b zK7L34-7owrx`6U8=zn>Mzr%mG=l>+YK=Uu|{jYT|F9rD~9)BFefCqHF`JrwsfAszz Dx=Ub( literal 0 HcmV?d00001