Skip to content

Commit f139c32

Browse files
bsouhaibbsouhaib
committedApr 13, 2015
first commit
0 parents  commit f139c32

8 files changed

+1552
-0
lines changed
 

‎README.md

Whitespace-only changes.

‎config.R

+11
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
library(nnet)
2+
library(mboost)
3+
library(gtools)
4+
5+
source("multistep.R")
6+
source("methods.R")
7+
source("utils.R")
8+
9+
# At the moment Nearest neighbors models are not available
10+
# source("mykknn.R")
11+

‎main-package.R

+51
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
rm(list=ls())
2+
source("config.R")
3+
source("strategies.R")
4+
5+
library(forecast)
6+
7+
myts <- window(sunspot.year,start=1700,end=1988-11)
8+
future <- as.numeric(window(sunspot.year,start=1988-11,1988))
9+
10+
# Forecasting horizon
11+
H <- 12
12+
13+
# Cross validation procedures
14+
control <- strategy_control("ts-cv",train.percentage = 0.7, n.fold=5)
15+
16+
# Learners
17+
bst1 <- strategy_learner(name = "BST", interactions=1, nu = 0.3, max.mstop = 500)
18+
bst2 <- strategy_learner(name = "BST", interactions=2, nu = 0.3, max.mstop = 500)
19+
20+
# KNN NOT AVAILABLE AT THIS STAGE
21+
# knn <- strategy_learner(name = "KNN")
22+
23+
mlp <- strategy_learner(name = "MLP", set.hidden = c(0,1,2,3,5), set.decay = c(0.005,0.01,0.05,0.1,0.2,0.3), nb.runs = 1, maxiter = 100, drop.linbias = F)
24+
lin <- strategy_learner(name = "MLP", set.hidden = 0, set.decay = 0, nb.runs = 1, drop.linbias = F)
25+
26+
# Recursive forecasts
27+
resultsREC <- recursive(time.series = myts, H = H, set.embedding = c(1,2,3,4,5), objectives = details$objectives, learner = mlp, control = control)
28+
frec <- as.numeric(resultsREC$forecasts)
29+
30+
# Direct forecasts
31+
resultsDIR <- direct(time.series = myts, H = H, set.embedding = c(1,2,3,4,5), objectives = details$objectives, learner = mlp, control = control)
32+
fdir <- as.numeric(resultsDIR$forecasts)
33+
34+
# Boosted AR forecasts
35+
resultsRECT <- rectify(time.series = myts, H = H, embeddings.base = c(1,2,3,4,5), embeddings.rect = c(1,2,3,4,5), objectives = details$objectives, learner.bmodel = lin, learner = bst1, control = control)
36+
fbase <- as.numeric(resultsRECT$base.forecasts)
37+
frect <- as.numeric(resultsRECT$forecasts)
38+
39+
# Auto.arima forecasts
40+
farima <- as.numeric(forecast(auto.arima(myts), h = 12)$mean)
41+
42+
43+
# Base AR forecasts in red
44+
# Boosting AR forecasts in blue
45+
# Real future in green
46+
matplot(cbind(fbase, frect, future), col = c("red", "blue", "green"), type = 'l', lty = 1)
47+
48+
# Comparing recursive, direct, arima and boosted ar forecasts
49+
mat <- cbind(future, frec, fdir, farima, frect)
50+
mycol <- rainbow(ncol(mat))
51+
matplot(mat, col = mycol, type = 'l', lty = 1)

‎methods.R

+236
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,236 @@
1+
2+
strategy_learner <- function(name=c("KNN","MLP", "MLP", "BST"),...) {
3+
4+
name <- match.arg(name)
5+
6+
parameters <- list(...)
7+
par.names <- names(parameters)
8+
9+
if(name=="KNN"){
10+
11+
if(!is.element("inc",par.names)){
12+
13+
parameters$inc <- 2
14+
}
15+
16+
}else if(name=="MLP"){
17+
18+
if(!is.element("set.hidden",par.names)){
19+
20+
parameters$set.hidden <- c(0, 1, 2, 3, 4)
21+
}
22+
23+
if(!is.element("set.decay",par.names)){
24+
25+
parameters$set.decay <- c(0, 0.001, 0.1, 0.5)
26+
}
27+
28+
if(!is.element("n.runs",par.names)){
29+
30+
parameters$n.runs <- 1
31+
}
32+
33+
if(!is.element("maxiter",par.names)){
34+
parameters$maxiter <- 100
35+
}
36+
37+
##### Checking #####
38+
39+
if(parameters$n.runs<1)
40+
stop("n.runs must be > 0 !")
41+
42+
if( any(parameters$set.decay<0) || any(parameters$set.hidden<0) )
43+
stop(" All elements of set.decay and set.hidden must be >= 0 ")
44+
45+
}else if (name=="BST"){
46+
47+
if(!is.element("interactions",par.names)){
48+
49+
parameters$interactions <- 1
50+
}
51+
52+
if(!is.element("nu",par.names)){
53+
54+
parameters$nu <- 0.2
55+
}
56+
57+
if(!is.element("max.mstop",par.names)){
58+
59+
parameters$max.mstop <- 100
60+
}
61+
62+
if(!is.element("save.all",par.names)){
63+
64+
parameters$save.all <- FALSE
65+
}
66+
67+
##### Checking #####
68+
if(!(parameters$nu>=0 && parameters$nu<=1))
69+
stop("nu must be in [0,1] !")
70+
71+
72+
}
73+
74+
RET <- c(list(name = name),parameters)
75+
class(RET) <- c("learner")
76+
RET
77+
}
78+
79+
80+
strategy_procedure <- function(type=c("recursive","direct"), objectives = NULL) {
81+
82+
type <- match.arg(type)
83+
84+
RET <- list(type = type, objectives = objectives)
85+
class(RET) <- c("procedure")
86+
RET
87+
}
88+
89+
90+
strategy_control <- function(method = c("ts-cv","cv"),...) {
91+
92+
method <- match.arg(method)
93+
parameters <- list(...)
94+
par.names <- names(parameters)
95+
96+
if(method == "ts-cv"){
97+
98+
if(!is.element("train.percentage",par.names)){
99+
100+
parameters$train.percentage <- 0.7
101+
}
102+
}
103+
104+
if(!is.element("n.fold",par.names)){
105+
106+
parameters$n.fold <- 5
107+
}
108+
109+
validation <- c(list(method = method), parameters)
110+
111+
RET <- list(validation = validation)
112+
class(RET) <- c("control")
113+
RET
114+
}
115+
116+
117+
strategy <- function(time.series , H, set.embedding = seq(1,5), data.set=NULL, procedure=strategy_procedure(), learner = strategy_learner(), control=strategy_control()){
118+
119+
####### TESTING INTEGRITY OF STRATEGY ########
120+
stopifnot(length(time.series) > 0)
121+
stopifnot(H>0)
122+
stopifnot(max(set.embedding) < length(time.series))
123+
124+
if(any(set.embedding<1)){
125+
stop(" All elements of set.embedding must be > 0 ")
126+
}
127+
128+
if( procedure$type=="recursive" && learner$name!="KNN" && learner$name!="MLP" && learner$name!="BST" && !is.null(procedure$objectives) ){
129+
stop(" Only KNN can be used with a recursive procedure and objectives != NULL ")
130+
}
131+
132+
if( learner$name=="BST" && learner$interactions==2 && max(set.embedding) == 1 ){
133+
stop(" To allow two interactions (interactions=2), all elements of set.embedding must be > 1 ! ")
134+
}
135+
136+
137+
#if(learner$name=="BST" && !is.null(procedure$objectives)){
138+
#
139+
# if( !identical( as.list(procedure$objectives), lapply(seq(H),identity)) ){
140+
# stop(" With BST, objectives must be seq(H) or be equal to NULL ! ")
141+
# }
142+
#
143+
# }
144+
############################################
145+
146+
147+
max.embedding <- max(set.embedding)
148+
y.label <- "Y"
149+
x.label <- "X"
150+
151+
all.x.variables <- paste(x.label,seq(max.embedding),sep="")
152+
all.y.variables <- paste(y.label,seq(H,1),sep="")
153+
154+
if(is.null(data.set)){
155+
data.set <- EmbedTimeSeries(time.series, c(seq(H-1, 0), -seq(max.embedding)), y.label, x.label)
156+
}
157+
158+
# CHECK HERE THAT h is in the objective for direct and ALL the objectives are all in seq(H)
159+
160+
161+
if( procedure$type=="direct" && is.null(procedure$objectives) ){
162+
163+
procedure$objectives <- lapply(seq(H),identity)
164+
}
165+
procedure$unique.objectives <- unique(procedure$objectives)
166+
167+
168+
169+
170+
RET <- list(time.series=time.series,
171+
H=H,
172+
data.set=data.set,
173+
n.init.obs=nrow(data.set),
174+
set.embedding=set.embedding,
175+
n.embedding=length(set.embedding),
176+
max.embedding=max(set.embedding),
177+
all.x.variables=all.x.variables,
178+
all.y.variables=all.y.variables,
179+
procedure=procedure,
180+
learner=learner,
181+
control=control)
182+
183+
class(RET) <- c("strategy")
184+
RET
185+
}
186+
187+
188+
189+
type.strategy <- function(x){
190+
x$procedure$type
191+
}
192+
193+
multistep.matrix.strategy <- function(x){
194+
195+
if(x$learner$name!="KNN")
196+
stop(" This function is only available for KNN !")
197+
198+
error.matrix <- ifelse(type(x)=="recursive",
199+
recursive.matrix(x),
200+
direct.matrix(x))
201+
202+
error.matrix
203+
}
204+
205+
206+
print.learner <- function(x){
207+
cat("Learner : ",x$name,"\n")
208+
invisible(x)
209+
}
210+
211+
print.procedure <- function(x){
212+
cat("Type : ",x$type,"\n")
213+
cat("Criterion : ",x$criterion,"\n")
214+
invisible(x)
215+
}
216+
217+
218+
print.control <- function(x){
219+
cat("Validation : ",x$validation,"\n")
220+
}
221+
222+
223+
print.strategy <- function(x){
224+
225+
cat("Length : ",length(x$time.series),"\n")
226+
cat("H : ",x$H,"\n")
227+
cat("Set.embedding : ",x$set.embedding,"\n")
228+
229+
print(x$procedure)
230+
cat("--- \n")
231+
print(x$learner)
232+
cat("--- \n")
233+
print(x$control)
234+
235+
invisible(x)
236+
}

‎multistep.R

+717
Large diffs are not rendered by default.

‎shared-main.R

+162
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
2+
# This main file is for running multiple strategies and models
3+
4+
knn.rec.matrix <- knn.dir.matrix <- NULL
5+
mlp.rec.matrix <- mlp.dir.matrix <- NULL
6+
lin.rec.matrix <- lin.dir.matrix <- NULL
7+
bst1.rec.matrix <- bst1.dir.matrix <- NULL
8+
bst2.rec.matrix <- bst2.dir.matrix <- NULL
9+
10+
11+
rec.matrix <- dir.matrix <- NULL
12+
base.results <- NULL
13+
14+
init.set.embedding <- set.embedding
15+
init.embeddings.base <- embeddings.base
16+
init.embeddings.rect <- embeddings.rect
17+
18+
for(i in seq_along(strategies)){
19+
20+
istrategy <- strategies[i]
21+
print(istrategy)
22+
23+
if(istrategy!="MEAN" && istrategy!="NAIVE"){
24+
25+
26+
###################
27+
# REC-LIN will have the same forecasts as the base forecasts of all RFY
28+
details <- GetDetails(ifelse(istrategy=="REC-LIN", "RFY-NONE", istrategy), H = H)
29+
30+
31+
learner <- switch(details$learner,
32+
KNN = knn,
33+
MLP = mlp,
34+
MLPmo = mlp,
35+
BST1 = bst1,
36+
BST2 = bst2,
37+
LIN = lin,
38+
LINMIS = lin,
39+
NONE = NULL)
40+
41+
learner.bmodel <- lin
42+
43+
set.embedding <- init.set.embedding
44+
embeddings.base <- init.embeddings.base
45+
embeddings.rect <- init.embeddings.rect
46+
47+
if(exists("DGP") && DGP == "SUNSPOT"){
48+
49+
if(details$learner == "LINMIS"){
50+
set.embedding <- embeddings.base <- embeddings.rect <- c(1, 2, 3)
51+
}
52+
if(details$procedure == "RFYMIS"){
53+
embeddings.base <- c(1, 2, 3)
54+
}
55+
56+
}
57+
#print(set.embedding)
58+
#print(embeddings.base)
59+
#print(embeddings.rect)
60+
61+
###################
62+
63+
64+
if(details$learner == "LIN")
65+
{
66+
rec.matrix <- lin.rec.matrix
67+
dir.matrix <- lin.dir.matrix
68+
69+
}else if(details$learner == "MLP")
70+
{
71+
rec.matrix <- mlp.rec.matrix
72+
dir.matrix <- mlp.dir.matrix
73+
74+
}else if(details$learner == "KNN"){
75+
76+
rec.matrix <- knn.rec.matrix
77+
dir.matrix <- knn.dir.matrix
78+
79+
}else if(details$learner == "BST1"){
80+
81+
rec.matrix <- bst1.rec.matrix
82+
dir.matrix <- bst1.dir.matrix
83+
84+
}else if(details$learner == "BST2"){
85+
86+
rec.matrix <- bst2.rec.matrix
87+
dir.matrix <- bst2.dir.matrix
88+
89+
}else{
90+
rec.matrix <- dir.matrix <- NULL
91+
}
92+
93+
# To avoid problems with RFYJT
94+
if(grepl("RFY",details$procedure)){
95+
rec.matrix <- dir.matrix <- NULL
96+
}
97+
98+
if(details$proc.type == "recursive"){
99+
100+
results <- recursive(time.series = trainset, H=H, set.embedding = set.embedding, objectives = details$objectives, Xtest = Xtest, error.matrix = rec.matrix,
101+
learner = learner, control = control)
102+
forecasts <- results$forecasts
103+
104+
if(details$learner == "LIN" && is.null(lin.rec.matrix) ){
105+
lin.rec.matrix <- results$error.matrix
106+
}else if(details$learner == "MLP" && is.null(mlp.rec.matrix) ){
107+
mlp.rec.matrix <- results$error.matrix
108+
}else if(details$learner == "KNN" && is.null(knn.rec.matrix)){
109+
knn.rec.matrix <- results$error.matrix
110+
}else if(details$learner == "BST1" && is.null(bst1.rec.matrix)){
111+
bst1.rec.matrix <- results$error.matrix
112+
}else if(details$learner == "BST2" && is.null(bst2.rec.matrix)){
113+
bst2.rec.matrix <- results$error.matrix
114+
}
115+
116+
}else if(details$proc.type == "direct"){
117+
118+
results <- direct(time.series = trainset, H=H, set.embedding = set.embedding, objectives = details$objectives, Xtest = Xtest, error.matrix = dir.matrix,
119+
learner = learner, control = control)
120+
forecasts <- results$forecasts
121+
122+
if(details$learner == "LIN" && is.null(lin.dir.matrix) ){
123+
lin.dir.matrix <- results$error.matrix
124+
}else if(details$learner == "MLP" && is.null(mlp.dir.matrix) ){
125+
mlp.dir.matrix <- results$error.matrix
126+
}else if(details$learner == "KNN" && is.null(knn.dir.matrix)){
127+
knn.dir.matrix <- results$error.matrix
128+
}else if(details$learner == "BST1" && is.null(bst1.dir.matrix)){
129+
bst1.dir.matrix <- results$error.matrix
130+
}else if(details$learner == "BST2" && is.null(bst2.dir.matrix)){
131+
bst2.dir.matrix <- results$error.matrix
132+
}
133+
134+
}else if(details$proc.type == "rectify"){
135+
136+
137+
results <- rectify(time.series = trainset, H=H, embeddings.base = embeddings.base, embeddings.rect = embeddings.rect, objectives = details$objectives, Xtest = Xtest, base.results = base.results, learner.bmodel = learner.bmodel,
138+
learner = learner, control = control)
139+
140+
#if(is.null(base.results)){
141+
# base.results <- results$base.results
142+
#}
143+
forecasts <- results$forecasts
144+
}
145+
146+
}else if(istrategy == "MEAN"){
147+
148+
mforecasts<-rep(mean(trainset),H)
149+
forecasts <- matrix(rep(mforecasts,Nbtest),nrow=Nbtest,ncol=H,byrow=T)
150+
}else if(istrategy == "NAIVE"){
151+
naiveforecasts<-rep(tail(trainset, 1),H)
152+
forecasts <- matrix(rep(naiveforecasts,Nbtest),nrow=Nbtest,ncol=H,byrow=T)
153+
}
154+
155+
all.forecasts[[i]]$forecasts[, seq(H), id.run] <- as.matrix(forecasts)
156+
157+
if(grepl("RFY",istrategy)){
158+
all.forecasts[[i]]$comp1[, seq(H), id.run] <- as.matrix(results$base.forecasts)
159+
all.forecasts[[i]]$comp2[, seq(H), id.run] <- as.matrix(results$rectifications)
160+
}
161+
162+
}

‎strategies.R

+97
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
2+
recursive <- function(time.series, H, set.embedding, objectives = NULL, Xtest = NULL, error.matrix = NULL, ...){
3+
4+
obj <- strategy(time.series, H, set.embedding = set.embedding, procedure=strategy_procedure(type="recursive", objectives = objectives), ...)
5+
learning <- learn(obj, error.matrix)
6+
7+
forecasts <- predict(obj, learning$results, Xtest = Xtest)$forecasts
8+
list(forecasts=forecasts, error.matrix = learning$error.matrix)
9+
}
10+
11+
direct <- function(time.series, H, set.embedding, data.set = NULL , objectives = NULL, Xtest = NULL, error.matrix = NULL, ...){
12+
13+
obj <- strategy(time.series, H, set.embedding = set.embedding, data.set, procedure=strategy_procedure(type="direct", objectives = objectives), ...)
14+
learning <- learn(obj, error.matrix)
15+
16+
forecasts <- predict(obj, learning$results, Xtest = Xtest)$forecasts
17+
list(forecasts=forecasts, error.matrix = learning$error.matrix)
18+
}
19+
20+
rectify <-function(time.series, H, embeddings.base, embeddings.rect, objectives = NULL, Xtest = NULL, base.results = NULL, learner.bmodel, ...){
21+
22+
args <- list(...)
23+
24+
if(is.null(base.results)){
25+
############## Base model ##############
26+
obj <- strategy(time.series, H, set.embedding = embeddings.base,
27+
learner = learner.bmodel,
28+
procedure = strategy_procedure(type = "recursive"),
29+
control = args$control)
30+
learning <- learn(obj)
31+
base.forecasts <- predict(obj, learning$results, Xtest)$forecasts
32+
33+
34+
id.complete <- which(complete.cases(obj$data.set[obj$all.x.variables]))
35+
queries <- obj$data.set[id.complete, obj$all.x.variables, drop=F]
36+
predictions <- predict(obj, learning$results, Xtest = queries)$forecasts
37+
38+
# This code require the same embedding for the base model and the rectification models
39+
# pseudo.data <- obj$data.set
40+
# pseudo.data[obj$all.y.variables] <- NA
41+
# pseudo.data[id.complete, obj$all.y.variables] <- obj$data.set[id.complete, obj$all.y.variables, drop=F] - rev(predictions)
42+
43+
#### This code allow a different embedding for rectifications
44+
max.embedding.rect <- max(embeddings.rect)
45+
x.label <- "X"; y.label <- "Y";
46+
data.set <- EmbedTimeSeries(time.series, c(seq(H-1, 0), -seq(max.embedding.rect)), y.label, x.label)
47+
all.x.variables <- paste(x.label,seq(max.embedding.rect),sep="")
48+
all.y.variables <- paste(y.label,seq(H,1),sep="")
49+
pseudo.data <- data.set
50+
pseudo.data[all.y.variables] <- NA
51+
pseudo.data[id.complete, all.y.variables] <- data.set[id.complete, all.y.variables, drop=F] - rev(predictions)
52+
53+
54+
base.results <- list(base.forecasts = base.forecasts, pseudo.data = pseudo.data)
55+
}
56+
57+
base.forecasts <- base.results$base.forecasts
58+
pseudo.data <- base.results$pseudo.data
59+
60+
if(is.null(Xtest)){
61+
stopifnot(all(dim(base.forecasts)==c(1,H)))
62+
}else{
63+
stopifnot(all(dim(base.forecasts)==c(nrow(Xtest),H)))
64+
}
65+
66+
############## Rectification models ##############
67+
rectifications <- 0
68+
69+
if(!is.null(learner)){
70+
rectifications <- direct(time.series, H, set.embedding = embeddings.rect, data.set=pseudo.data, objectives = objectives, Xtest = Xtest, ...)$forecasts
71+
}
72+
73+
forecasts <- base.forecasts + rectifications
74+
75+
76+
list(forecasts = forecasts, base.forecasts = base.forecasts, rectifications = rectifications, base.results = base.results)
77+
}
78+
79+
80+
rec <- function(time.series, H, ...){
81+
82+
recursive(time.series, H, ...)$forecasts
83+
}
84+
85+
dir <- function(time.series, H, ...){
86+
87+
direct(time.series, H, ...)$forecasts
88+
}
89+
90+
mystrategy <- function(time.series, H, ...){
91+
obj <- strategy(time.series, H, ...)
92+
results <- learn(obj)
93+
forecasts <- predict(obj,results)$forecasts
94+
forecasts
95+
}
96+
97+

‎utils.R

+278
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,278 @@
1+
mynnet <- function(config, ...)
2+
{
3+
# config contient obj$learner
4+
5+
all.param <- list(...)
6+
drop.linbias <- config$drop.linbias
7+
is.linear <- (all.param$size == 0)
8+
skip <- ifelse(is.linear,TRUE,FALSE)
9+
maxit <- config$maxiter
10+
11+
if(is.linear && drop.linbias)
12+
{
13+
n.inputs <- ncol(all.param$x)
14+
rang <- 0.7
15+
mask <- c(0, rep(1,n.inputs))
16+
nwts <- n.inputs + 1
17+
Wts <- runif(nwts, -rang, rang)
18+
Wts[1] <- 0
19+
20+
model <- nnet(..., mask = mask, Wts = Wts,
21+
linout = T,trace = F, maxit = maxit, skip = skip)
22+
23+
24+
}else{
25+
model <- nnet(...,
26+
linout = T,trace = F, maxit = maxit, skip = skip)
27+
}
28+
model
29+
}
30+
31+
32+
33+
predictions.bst <- function(mstop, model, X){
34+
35+
if(mstop > 0){
36+
predictions <- predict(model, X)
37+
}else if(mstop == 0){
38+
predictions <- rep(model$offset, nrow(X))
39+
}else if(mstop == -1){
40+
predictions <- rep(0, nrow(X))
41+
}
42+
predictions
43+
}
44+
45+
46+
#if(best.model$mstop > 0){
47+
# forecasts[,horizon.variables] <- predict(best.model$model, Xtest)
48+
#
49+
# if(learner$save.all){
50+
# details$best.mstop[id.objectif] <- best.model$mstop
51+
# details$all.iterations[,seq(best.model$mstop),id.objectif] <- predict(best.model$model,Xtest,aggregate="none")
52+
# }
53+
#
54+
# }else if(best.model$mstop == 0){
55+
# forecasts[,horizon.variables] <- rep(best.model$model$offset, nrow(Xtest))
56+
#
57+
# if(learner$save.all){
58+
# details$best.mstop[id.objectif] <- 0
59+
# }
60+
# }else if(best.model$mstop == -1){
61+
# forecasts[,horizon.variables] <- rep(0, nrow(Xtest))
62+
# if(learner$save.all){
63+
# details$best.mstop[id.objectif] <- -1
64+
# }
65+
# }
66+
67+
68+
EmbedTimeSeries<-function(x,lags,y.label,x.label){
69+
70+
n.obs<-length(x)
71+
72+
# Remove duplicate lags
73+
lags<-union(lags, lags)
74+
75+
lags<-sort(lags, decreasing=T)
76+
index<-rep(seq(n.obs)) + rep(lags, each=n.obs)
77+
id.out<-which(index <= 0 | index > n.obs)
78+
79+
if(length(id.out) > 0){
80+
index[id.out] <- 1
81+
}
82+
83+
dataset<-x[index]
84+
85+
if(length(id.out) > 0){
86+
dataset[id.out] <- NA
87+
}
88+
89+
dim(dataset)<-c(n.obs, length(lags))
90+
dataset <- data.frame(dataset)
91+
colnames(dataset)<-ifelse(lags >= 0, paste(y.label, lags+1 ,sep=""), paste(x.label, abs(lags), sep=""))
92+
dataset
93+
}
94+
95+
MakeFormula <-function(y.variable, x.variables, interactions=1,nbKnotsUni=20,nbKnotsBiv=5)
96+
{
97+
stopifnot(interactions%in%c(1,2,12))
98+
99+
n.variables<-length(x.variables)
100+
101+
102+
if(interactions==1 || interactions==12)
103+
{
104+
myformula<-uniformula<- paste("bbs(",x.variables,",df=4,center=TRUE,knots=",nbKnotsUni,")",sep="",collapse="+")
105+
}
106+
107+
108+
if(interactions==2 || interactions==12)
109+
{
110+
111+
bivformula<-""
112+
allcomb <- combinations(n.variables,2,v=x.variables)
113+
for(icomb in seq(nrow(allcomb)))
114+
{
115+
term<-paste("bbs(",paste(allcomb[icomb,],collapse=",",sep=""),",df=4,center=TRUE,",
116+
paste("knots=list(",paste(allcomb[icomb,],"=",nbKnotsBiv,collapse=",",sep=""),")",sep=""),")",sep="")
117+
118+
bivformula <- ifelse(icomb==1, term, paste(bivformula,term,sep="+"))
119+
}
120+
myformula <- bivformula
121+
122+
if(interactions==12)
123+
{
124+
myformula<-paste(uniformula,bivformula,sep="+")
125+
}
126+
}
127+
finalformula<-paste(y.variable,"~",myformula,sep="")
128+
finalformula
129+
}
130+
131+
132+
ValidationMatrix <- function(validation, n.obs)
133+
{
134+
135+
if(validation$method == "ts-cv"){
136+
137+
n.train <- floor(n.obs * validation$train.percentage)
138+
n.valid <- n.obs - n.train
139+
140+
if(n.valid < validation$n.fold){
141+
stop(" The validation set is too short to have this number of folds ")
142+
}
143+
size.block <- floor(n.valid/ validation$n.fold)
144+
145+
folds <- unlist(lapply( seq(validation$n.fold), function(i.fold) {c( rep(1, n.train) ,rep( rep(1, size.block), i.fold -1 ) , rep(0, n.valid-(i.fold-1)*size.block ) )} ))
146+
validation.matrix <- matrix(folds, nrow = n.obs)
147+
stopifnot( ncol(validation.matrix) == validation$n.fold )
148+
} else if(validation$method=="cv"){
149+
150+
cvkfold <- function(n, k) {
151+
if (k > n / 2) stop("k > n/2")
152+
fl <- floor(n/k)
153+
folds <- c(rep(c(rep(0, fl), rep(1, n)), k - 1),
154+
rep(0, n * k - (k - 1) * (fl + n)))
155+
matrix(folds, nrow = n)[1:n,, drop = FALSE]
156+
}
157+
158+
validation.matrix <- cvkfold(n.obs, validation$n.fold)
159+
160+
}else{
161+
stop("error in validation$method !")
162+
}
163+
validation.matrix
164+
}
165+
166+
167+
getObjectifJTs <- function(s1, s2, H)
168+
{
169+
RET <- removeOut( lapply(seq(H),"+",seq(-s1,s2)) , H)
170+
RET
171+
}
172+
173+
getS <- function(procedure){
174+
n <- nchar(procedure)
175+
s1 <- as.numeric(substring(procedure, n-1, n-1))
176+
s2 <- as.numeric(substring(procedure, n, n))
177+
list(s1 = s1, s2 = s2)
178+
}
179+
180+
removeOut <- function(mylist,H){
181+
182+
for( i in seq(length(mylist)) ){
183+
vec <- mylist[[i]]
184+
id <- which( !(vec %in% seq(H)))
185+
if(length(id)>0){
186+
RET <- vec[-id]
187+
}else{
188+
RET <- vec
189+
}
190+
mylist[[i]] <- RET
191+
}
192+
mylist
193+
}
194+
195+
GetDetails <- function(strategy,H)
196+
{
197+
res <- unlist(strsplit(strategy,"-"))
198+
procedure <- res[1]
199+
learner <- res[2]
200+
201+
if(procedure == "REC"){
202+
203+
objectives <- NULL
204+
205+
}else if(procedure %in% c("RTI", "DIR", "RFY", "RFYMIS")){
206+
207+
objectives <- seq(H)
208+
209+
}else if( grepl("JTL",procedure) ){
210+
# RJTLab DJTLab
211+
212+
res <- getS(procedure)
213+
objectives <- getObjectifJTs(res$s1, res$s2, H)
214+
215+
}else if( procedure %in% c("RJT", "DJT")){
216+
217+
objectives <- lapply(lapply(rep(H,H),identity),seq)
218+
219+
}else if( grepl("RJT",procedure) )
220+
{
221+
# RJTa
222+
n <- nchar(procedure)
223+
224+
mychar <- substring(procedure, 4, n)
225+
s <- as.numeric(mychar)
226+
objectives <- lapply(rep(s,H),seq)
227+
}
228+
229+
230+
# TO DO : Check if the procedure belongs to the set of possible procedures
231+
232+
if(procedure == "RFY" || procedure == "RFYMIS" || procedure == "RFYJT"){
233+
proc.type <- "rectify"
234+
}else if( grepl("R",procedure) && procedure != "DIR"){
235+
proc.type <- "recursive"
236+
}else if( grepl("D",procedure) ){
237+
proc.type <- "direct"
238+
}else{
239+
stop("procedure not found ! ")
240+
}
241+
242+
243+
details <- list(proc.type = proc.type, procedure = procedure, learner = learner, objectives = objectives)
244+
details
245+
}
246+
247+
248+
################################
249+
# improved list of objects
250+
# from http://stackoverflow.com/questions/1358003/tricks-to-manage-the-available-memory-in-an-r-session
251+
.ls.objects <- function (pos = 1, pattern, order.by,
252+
decreasing=FALSE, head=FALSE, n=5) {
253+
napply <- function(names, fn) sapply(names, function(x)
254+
fn(get(x, pos = pos)))
255+
names <- ls(pos = pos, pattern = pattern)
256+
obj.class <- napply(names, function(x) as.character(class(x))[1])
257+
obj.mode <- napply(names, mode)
258+
obj.type <- ifelse(is.na(obj.class), obj.mode, obj.class)
259+
obj.prettysize <- napply(names, function(x) {
260+
capture.output(print(object.size(x), units = "auto")) })
261+
obj.size <- napply(names, object.size)
262+
obj.dim <- t(napply(names, function(x)
263+
as.numeric(dim(x))[1:2]))
264+
vec <- is.na(obj.dim)[, 1] & (obj.type != "function")
265+
obj.dim[vec, 1] <- napply(names, length)[vec]
266+
out <- data.frame(obj.type, obj.size, obj.prettysize, obj.dim)
267+
names(out) <- c("Type", "Size", "PrettySize", "Rows", "Columns")
268+
if (!missing(order.by))
269+
out <- out[order(out[[order.by]], decreasing=decreasing), ]
270+
if (head)
271+
out <- head(out, n)
272+
out
273+
}
274+
275+
# shorthand
276+
lsos <- function(..., n=10) {
277+
.ls.objects(..., order.by="Size", decreasing=TRUE, head=TRUE, n=n)
278+
}

0 commit comments

Comments
 (0)
Please sign in to comment.