--- /dev/null
+ # HACK for test reports: complete some days with a few NAs, for nicer graphics
+ nas_in_serie = is.na(data$getSerie(today))
+ if (any(nas_in_serie))
+ {
+ if (sum(nas_in_serie) >= length(nas_in_serie) / 2)
+ return (NA)
+ for (i in seq_along(nas_in_serie))
+ {
+ if (nas_in_serie[i])
+ {
+ #look left
+ left = i-1
+ while (left>=1 && nas_in_serie[left])
+ left = left-1
+ #look right
+ right = i+1
+ while (right<=length(nas_in_serie) && nas_in_serie[right])
+ right = right+1
+ #HACK: modify by-reference Data object...
+ data$data[[today]]$serie[i] <<-
+ if (left==0) data$data[[today]]$serie[right]
+ else if (right==0) data$data[[today]]$serie[left]
+ else (data$data[[today]]$serie[left] + data$data[[today]]$serie[right]) / 2.
+ }
+ }
+ }
+
Bruno Portier <Bruno.Portier@insa-rouen.fr>, [ctb]
Maintainer: Benjamin Auder <Benjamin.Auder@math.u-psud.fr>
Depends:
- R (>= 3.0)
+ R (>= 3.0),
+ R6
Suggests:
roxygen2,
testthat,
'J_Neighbors.R'
'J_Persistence.R'
'J_Zero.R'
- 'getData.R'
'computeError.R'
'computeForecast.R'
+ 'getData.R'
'plot.R'
'utils.R'
VignetteBuilder: knitr
#' @importFrom R6 R6Class
#'
#' @export
-Data = R6Class("Data",
+Data = R6::R6Class("Data",
private = list(
- .data = "list"
+ .data = list()
),
public = list(
- initialize = function(...)
- initialize(self, private, ...)
- ,
getSize = function()
- getSize(private)
+ getSizeData(private)
,
getStdHorizon = function()
- getStdHorizon(private)
+ getStdHorizonData(private)
,
append = function(new_time, new_centered_serie, new_level, new_exo, new_exo_hat)
- append(private, new_time, new_centered_serie, new_level, new_exo, new_exo_hat)
+ appendData(private, new_time, new_centered_serie, new_level, new_exo, new_exo_hat)
,
getTime = function(index)
- getTime(self, private, index)
+ getTimeData(self, private, index)
,
getCenteredSerie = function(index)
- getCenteredSerie(self, private, index)
+ getCenteredSerieData(self, private, index)
,
getLevel = function(index)
- getLevel(self, private, index)
+ getLevelData(self, private, index)
,
getSerie = function(index)
- getSerie(self, private, index)
+ getSerieData(self, private, index)
,
getExo = function(index)
- getExo(self, private, index)
+ getExoData(self, private, index)
,
getExoHat = function(index)
- getExoHat(self, private, index)
+ getExoHatData(self, private, index)
)
)
-##TODO: @param... @inheritParams...
-
-#' Initialize empty Data object
-initialize = function(o, private, ...)
-{
- private$.data <<- if (hasArg("data")) list(...)$data else list()
- invisible(o)
-}
-
#' Number of series in the dataset
-getSize = function(private)
+#'
+#' @param private List of private members in current object
+getSizeData = function(private)
length(private$.data)
#' 'Standard' horizon, from t+1 to midnight
-getStdHorizon = function(private)
+#'
+#' @inheritParams getSizeData
+getStdHorizonData = function(private)
24 - as.POSIXlt( private$.data[[1]]$time[1] )$hour + 1
#' Acquire a new vector of lists (time, centered_serie, level, exo, exo_hat)
-append = function(private, new_time, new_centered_serie, new_level, new_exo_hat, new_exo)
+#'
+#' @inheritParams getSizeData
+#' @param new_time Time
+#' @param new_centered_serie Centered serie
+#' @param new_level Level
+#' @param new_exo Exogneous variables
+#' @param new_exo_hat Predicted exogenous variables
+appendData = function(private, new_time, new_centered_serie, new_level, new_exo, new_exo_hat)
{
- private$.data[[length(private$.data)+1]] <<- list("time"=new_time,
+ private$.data[[length(private$.data)+1]] <- list("time"=new_time,
"centered_serie"=new_centered_serie,"level"=new_level,"exo"=new_exo,"exo_hat"=new_exo_hat)
}
#' Time values at specified index
-getTime = function(o, private, index)
+#'
+#' @inheritParams getSizeData
+#' @param index Return value at this index
+getTimeData = function(o, private, index)
{
index = dateIndexToInteger(index, o)
private$.data[[index]]$time
}
#' Centered serie values at specified index
-getCenteredSerie = function(o, private, index)
+#'
+#' @inheritParams getTimeData
+getCenteredSerieData = function(o, private, index)
{
index = dateIndexToInteger(index, o)
- private$.data[[index]]$serie
+ private$.data[[index]]$centered_serie
}
#' Level of the serie at specified index
-getCenteredSerie = function(o, private, index)
+#'
+#' @inheritParams getTimeData
+getLevelData = function(o, private, index)
{
index = dateIndexToInteger(index, o)
private$.data[[index]]$level
}
#' Serie values (centered+level) at specified index
-getCenteredSerie = function(o, private, index)
+#'
+#' @inheritParams getTimeData
+getSerieData = function(o, private, index)
{
index = dateIndexToInteger(index, o)
- private$.data[[index]]$serie + data[[index]]$level
+ private$.data[[index]]$centered_serie + data[[index]]$level
}
#' Exogenous measures at specified index
-getCenteredSerie = function(o, private, index)
+#'
+#' @inheritParams getTimeData
+getExoData = function(o, private, index)
{
index = dateIndexToInteger(index, o)
private$.data[[index]]$exo
}
#' Exogeous predictions at specified index
-getCenteredSerie = function(o, private, index)
+#'
+#' @inheritParams getTimeData
+getExoHatData = function(o, private, index)
{
index = dateIndexToInteger(index, o)
private$.data[[index]]$exo_hat
#' @include Forecaster.R
#'
-#' @title Average Forecaster
+#' Average Forecaster
#'
-#' @description Return the (pointwise) average of the all the (similar) centered day curves
-#' in the past. Inherits \code{\link{Forecaster}}
-AverageForecaster = setRefClass(
- Class = "AverageForecaster",
- contains = "Forecaster",
+#' Return the (pointwise) average of the all the (similar) centered day curves
+#' in the past. Inherits \code{\link{Forecaster}}
+AverageForecaster = R6::R6Class("AverageForecaster",
+ inherit = "Forecaster",
- methods = list(
- initialize = function(...)
- {
- callSuper(...)
- },
- predict = function(today, memory, horizon, ...)
- {
- predicted_shape = predictShape(today, memory, horizon, ...)
- #Take care of never passing same_day==FALSE (when pjump == Persistence)
- predicted_delta =
- if (#as.character(substitute(pjump))=="Persistence" && #TODO: doesn't work
- hasArg("same_day") && list(...)$same_day==FALSE)
- {
- args = list(...)
- args$same_day = TRUE
- do.call(pjump, append(list("today"=today,"memory"=memory,"horizon"=horizon), args))
- }
- else
- pjump(data, today, memory, horizon, params, ...)
- predicted_shape + tail(data$getSerie(today),1) - predicted_shape[1] + predicted_delta
- },
+ public = list(
predictShape = function(today, memory, horizon, ...)
{
avg = rep(0., horizon)
#' @include Forecaster.R
#'
-#' @title Neighbors Forecaster
+#' Neighbors Forecaster
#'
-#' @description Predict tomorrow as a weighted combination of "futures of the past" days.
-#' Inherits \code{\link{Forecaster}}
-NeighborsForecaster = setRefClass(
- Class = "NeighborsForecaster",
- contains = "Forecaster",
-
- methods = list(
- initialize = function(...)
- {
- callSuper(...)
- },
+#' Predict tomorrow as a weighted combination of "futures of the past" days.
+#' Inherits \code{\link{Forecaster}}
+NeighborsForecaster = R6::R6Class("NeighborsForecaster",
+ inherit = "Forecaster",
+
+ public = list(
predictShape = function(today, memory, horizon, ...)
{
# (re)initialize computed parameters
if (hasArg(h_window))
return (.predictShapeAux(fdays,today,horizon,list(...)$h_window,kernel,simtype,TRUE))
- # HACK for test reports: complete some days with a few NAs, for nicer graphics
- nas_in_serie = is.na(data$getSerie(today))
- if (any(nas_in_serie))
- {
- if (sum(nas_in_serie) >= length(nas_in_serie) / 2)
- return (NA)
- for (i in seq_along(nas_in_serie))
- {
- if (nas_in_serie[i])
- {
- #look left
- left = i-1
- while (left>=1 && nas_in_serie[left])
- left = left-1
- #look right
- right = i+1
- while (right<=length(nas_in_serie) && nas_in_serie[right])
- right = right+1
- #HACK: modify by-reference Data object...
- data$data[[today]]$serie[i] <<-
- if (left==0) data$data[[today]]$serie[right]
- else if (right==0) data$data[[today]]$serie[left]
- else (data$data[[today]]$serie[left] + data$data[[today]]$serie[right]) / 2.
- }
- }
- }
-
# Determine indices of no-NAs days followed by no-NAs tomorrows
first_day = max(today - memory, 1)
fdays = (first_day:(today-1))[ sapply(first_day:(today-1), function(i) {
h_best_mix = c(h_best_endo,h_best_exo)
return (.predictShapeAux(fdays, today, horizon, h_best_mix, kernel, "mix", TRUE))
}
- },
+ }
+ ),
+ private = list(
# Precondition: "today" is full (no NAs)
.predictShapeAux = function(fdays, today, horizon, h, kernel, simtype, final_call)
{
- dat = data$data #HACK: faster this way...
-
fdays = fdays[ fdays < today ]
# TODO: 3 = magic number
if (length(fdays) < 3)
distances2 = rep(NA, length(fdays))
for (i in seq_along(fdays))
{
- delta = dat[[today]]$serie - dat[[ fdays[i] ]]$serie
+ delta = data$getCenteredSerie(today) - data$getCenteredSerie(fdays[i])
# Require at least half of non-NA common values to compute the distance
if (sum(is.na(delta)) <= 0) #length(delta)/2)
distances2[i] = mean(delta^2) #, na.rm=TRUE)
{
h_exo = ifelse(simtype=="mix", h[2], h)
- M = matrix( nrow=1+length(fdays), ncol=1+length(dat[[today]]$exo) )
- M[1,] = c( dat[[today]]$level, as.double(dat[[today]]$exo) )
+ M = matrix( nrow=1+length(fdays), ncol=1+length(data$getExo(today)) )
+ M[1,] = c( data$getLevel(today), as.double(data$getExo(today)) )
for (i in seq_along(fdays))
- M[i+1,] = c( dat[[ fdays[i] ]]$level, as.double(dat[[ fdays[i] ]]$exo) )
+ M[i+1,] = c( data$getLevel(fdays[i]), as.double(data$getExo(fdays[i])) )
sigma = cov(M) #NOTE: robust covariance is way too slow
sigma_inv = solve(sigma) #TODO: use pseudo-inverse if needed?
prediction = rep(0, horizon)
for (i in seq_along(fdays_indices))
- prediction = prediction + similarities[i] * dat[[ fdays_indices[i]+1 ]]$serie[1:horizon]
+ prediction = prediction + similarities[i] * data$getSerie(fdays_indices[i]+1)[1:horizon]
prediction = prediction / sum(similarities, na.rm=TRUE)
if (final_call)
#' @include Forecaster.R
#'
-#' @title Persistence Forecaster
+#' Persistence Forecaster
#'
-#' @description Return the last centered (similar) day curve.
-#' Inherits \code{\link{Forecaster}}
-PersistenceForecaster = setRefClass(
- Class = "PersistenceForecaster",
- contains = "Forecaster",
+#' Return the last centered (similar) day curve.
+#' Inherits \code{\link{Forecaster}}
+PersistenceForecaster = R6::R6Class("PersistenceForecaster",
+ inherit = "Forecaster",
- methods = list(
- initialize = function(...)
- {
- callSuper(...)
- },
+ public = list(
predictShape = function(today, memory, horizon, ...)
{
# Return centered last (similar) day curve, avoiding NAs until memory is run
#' @include Forecaster.R
#'
-#' @title Zero Forecaster
+#' Zero Forecaster
#'
-#' @description Return 0 (and then adjust). Inherits \code{\link{Forecaster}}
-ZeroForecaster = setRefClass(
- Class = "ZeroForecaster",
- contains = "Forecaster",
+#' Return 0 (and then adjust). Inherits \code{\link{Forecaster}}
+ZeroForecaster = R6::R6Class("ZeroForecaster",
+ inherit = "Forecaster",
- methods = list(
- initialize = function(...)
- {
- callSuper(...)
- },
+ public = list(
predictShape = function(today, memory, horizon, ...)
- {
rep(0., horizon)
- }
)
)
-#' @title Forecast
+#' Forecast
#'
-#' @description Forecast encapsulation
+#' Forecast encapsulation
#'
#' @field pred List with
#' \itemize{
#' \item index: corresponding index in data object
#' }
#'
-#' @exportClass Forecast
-#' @export Forecast
-Forecast = setRefClass(
- Class = "Forecast",
-
- fields = list(
- pred = "list"
+#' @docType class
+#' @importFrom R6 R6Class
+#'
+#' @export
+Forecast = R6::R6Class("Forecast",
+ private = list(
+ .pred = "list",
+ .dates = "Date"
),
-
- methods = list(
- initialize = function(...)
- {
- "Initialize empty Forecast object"
-
- callSuper(...)
- },
+ public = list(
+ initialize = function(dates)
+ initialize(self, private, dates)
+ ,
getSize = function()
- {
- "Number of individual forecasts"
-
- length(pred)
- },
+ getSize(private)
+ ,
append = function(new_serie, new_params, new_index)
- {
- "Obtain a new pair (serie, params)"
-
- pred[[length(pred)+1]] <<- list("serie"=new_serie, "params"=new_params, "index"=new_index)
- },
+ append(private, new_serie, new_params, new_index)
+ ,
+ getDates = function()
+ getDates(private)
+ ,
getSerie = function(index)
- {
- "Serie values at specified index"
-
- pred[[index]]$serie
- },
+ getSerie(private, index)
+ ,
getParams = function(index)
- {
- "Params at specified index"
-
- pred[[index]]$params
- },
+ getParams(private, index)
+ ,
getIndexInData = function(index)
- {
- "(day) Index in data where prediction took place"
-
- pred[[index]]$index
- }
+ getIndexInData(private, index)
)
)
+
+#' Initialize empty Forecast object
+#'
+#' @param o Object of class Forecast
+#' @param private List of private members in o
+#' @param dates vector of dates where forecast occurs
+initialize = function(o, private, dates)
+{
+ private$.dates <<- dates
+ private$.pred <<- list()
+ invisible(o)
+}
+
+#' Number of individual forecasts"
+#'
+#' @inheritParams initialize
+getSize = function(private)
+ length(private$.pred)
+
+#' Obtain a new pair (serie, params)"
+#'
+#' @inheritParams initialize
+append = function(new_serie, new_params, new_index_in_data)
+{
+ private$.pred[[length(private$.pred)+1]] <<-
+ list("serie"=new_serie, "params"=new_params, "index_in_data"=new_index_in_data)
+}
+
+#' Dates where prediction occurs
+#'
+#' inheritParams initialize
+getDates = function(private)
+ private$.dates
+
+#' Serie values at specified index"
+#'
+#' @inheritParams initialize
+#' @param index Return value at this index
+getSerie = function(index)
+{
+ if (is(index,"Date"))
+ index = match(index, private$.dates)
+ private$.pred[[index]]$serie
+}
+
+#' Params at specified index"
+#'
+#' @inheritParams getSerie
+getParams = function(index)
+{
+ if (is(index,"Date"))
+ index = match(index, private$.dates)
+ private$.pred[[index]]$params
+}
+
+#' (day) Index in data where prediction took place"
+#'
+#' @inheritParams getSerie
+getIndexInData = function(index)
+{
+ if (is(index,"Date"))
+ index = match(index, private$.dates)
+ private$.pred[[index]]$index_in_data
+}
-#' @title Forecaster (abstract class)
+#' Forecaster
#'
-#' @description Abstract class to represent a forecaster (they all inherit this)
+#' Forecaster (abstract class, implemented by all forecasters)
#'
#' @field params List of computed parameters, for post-run analysis (dev)
#' @field data Dataset, object of class Data
#' @field pjump Function: how to predict the jump at day interface ?
-Forecaster = setRefClass(
- Class = "Forecaster",
-
- fields = list(
- params = "list",
- data = "Data",
- pjump = "function"
+#'
+#' @docType class
+#' @importFrom R6 R6Class
+Forecaster = R6::R6Class("Forecaster",
+ private = list(
+ .params = "list",
+ .data = "Data",
+ .pjump = "function"
),
-
- methods = list(
- initialize = function(...)
- {
- "Initialize (generic) Forecaster object"
-
- callSuper(...)
- if (!hasArg(data))
- stop("Forecaster must be initialized with a Data object")
- params <<- list()
- },
- predict = function(today, memory, horizon, ...)
- {
- "Obtain a new forecasted time-serie"
-
- # Parameters (potentially) computed during shape prediction stage
- predicted_shape = predictShape(today, memory, horizon, ...)
- predicted_delta = pjump(data, today, memory, horizon, params, ...)
- # Predicted shape is aligned it on the end of current day + jump
- predicted_shape + tail(data$getSerie(today),1) - predicted_shape[1] + predicted_delta
- },
+ public = list(
+ initialize = function(data, pjump)
+ initialize(self, private, data, pjump)
+ ,
+ predictSerie = function(today, memory, horizon, ...)
+ predictSerie(private, today, memory, horizon, ...)
+ ,
predictShape = function(today, memory, horizon, ...)
- {
- "Shape prediction (centered curve)"
-
- #empty default implementation: to implement in inherited classes
- },
+ predictShape(private, today, memory, horizon, ...)
+ ,
getParameters = function()
- {
- "Get parameters list"
-
- params
- }
+ getParameters(private)
)
)
+
+#' Initialize (generic) Forecaster object
+#'
+#' @param o Object of class Forecaster
+#' @param private List of private members in o
+#' @param data Object of class Data
+#' @param pjump Function to predict jump
+initialize = function(o, private, data, pjump)
+{
+ .params <<- list()
+ .data <<- data
+ .pjump <<- pjump
+ invisible(o)
+}
+
+#' Obtain a new forecasted time-serie
+#'
+#' @inheritParams initialize
+#' @param today Index for current prediction
+#' @param memory Depth in data (in days)
+#' @param horizon Number of hours to forecast
+predictSerie = function(private, today, memory, horizon, ...)
+{
+ # Parameters (potentially) computed during shape prediction stage
+ predicted_shape = predictShape(today, memory, horizon, ...)
+ predicted_delta = private$.pjump(private$.data, today, memory, horizon, params, ...)
+ # Predicted shape is aligned it on the end of current day + jump
+ predicted_shape + tail(private$.data$getSerie(today),1) - predicted_shape[1] + predicted_delta
+}
+
+#' Shape prediction (centered curve)
+#'
+#' @inheritParams predictSerie
+predictShape = function(private, today, memory, horizon, ...)
+ #empty default implementation: to implement in inherited classes
+
+#' Get parameters list
+#'
+#' @inheritParams initialize
+getParameters = function(private)
+ private$.params
#' @param horizon Number of time steps to predict
#' @param ... Additional parameters for the forecasting models
#'
-#' @return An object of class Forecast
+#' @return A list with the following items
+#' \itemize{
+#' \item serie: forecasted serie
+#' \item params: corresponding list of parameters (weights, neighbors...)
+#' \item index: corresponding index in data object
+#' }
#'
#' @examples
#' ts_data = system.file("extdata","pm10_mesures_H_loc.csv",package="talweg")
#' #do_something_with_pred
#' }}
#' @export
-computeForecast = function(data, indices, forecaster, pjump=NULL,
+computeForecast = function(data, indices, forecaster, pjump,
memory=Inf, horizon=data$getStdHorizon(), ...)
{
# (basic) Arguments sanity checks
if (!is.character(forecaster))
stop("forecaster (name) should be of class character") #pjump could be NULL
- pred = list()
- forecaster = new(paste(forecaster,"Forecaster",sep=""), data=data,
- pjump =
- if (is.null(pjump))
- function() {}
- else
- getFromNamespace(paste("get",pjump,"JumpPredict",sep=""), "talweg"))
+ pred = Forecast$new()
+ forecaster_class_name = getFromNamespace(paste(forecaster,"Forecaster",sep=""), "talweg")
+ forecaster = forecaster_class_name$new(data=data,
+ pjump = getFromNamespace(paste("get",pjump,"JumpPredict",sep=""), "talweg"))
for (today in indices)
{
- #pred$append(...) is slow; TODO: use R6 class
- pred[[length(pred)+1]] = list(
- "serie" = forecaster$predict(today, memory, horizon, ...),
- "params" = forecaster$getParameters(),
- "index" = today
+ pred$append(
+ new_serie = forecaster$predictSerie(today, memory, horizon, ...),
+ new_params = forecaster$getParameters(),
+ new_index = today
)
}
- new("Forecast",pred=pred)
+ pred
}
#'
#' @description Take in input data frames and/or files containing raw data, and timezones, and
#' output a Data object, roughly corresponding to a list where each cell contains all value
-#' for one day (see \code{?Data}). Current limitation: series (in working_tz) must start at
-#' right after midnight (to keep in sync with exogenous vars)
+#' for one day (see \code{?Data}).
#'
#' @param ts_data Time-series, as a data frame (DB style: 2 columns, first is date/time,
#' second is value) or a CSV file
#' @param predict_at When does the prediction take place ? Integer, in hours. Default: 0
#' @param limit Number of days to extract (default: Inf, for "all")
#'
-#' @return An object of class Data
+#' @return A list where data[[i]] contains
+#' \itemize{
+#' \item time: vector of times
+#' \item centered_serie: centered serie
+#' \item level: corresponding level
+#' \item exo: exogenous variables
+#' \item exo_hat: predicted exogenous variables
+#' }
#'
#' @examples
#' ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc.csv",package="talweg"))
#' exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv",package="talweg"))
-#' getData(ts_data, exo_data, input_tz="Europe/Paris", working_tz="Europe/Paris", limit=150)
+#' data = getData(ts_data, exo_data, limit=120)
#' @export
getData = function(ts_data, exo_data, input_tz="GMT", date_format="%d/%m/%Y %H:%M",
working_tz="GMT", predict_at=0, limit=Inf)
line = 1 #index in PM10 file (24 lines for 1 cell)
nb_lines = nrow(ts_df)
nb_exos = ( ncol(exo_df) - 1 ) / 2
- data = Data$new()
+ data = list()
i = 1 #index of a cell in data
while (line <= nb_lines)
{
exo_hat = as.data.frame( exo_df[i,(1+nb_exos+1):(1+2*nb_exos)] )
level = mean(serie, na.rm=TRUE)
centered_serie = serie - level
- data$append(time, centered_serie, level, exo, exo_hat)
+ data[[i]] = list("time"=time, "centered_serie"=centered_serie, "level"=level,
+ "exo"=exo, "exo_hat"=exo_hat)
if (i >= limit)
break
i = i + 1
}
- if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(1)))
- data$removeFirst()
- if (length(data$getCenteredSerie( data$getSize() )) <
- length(data$getCenteredSerie( data$getSize()-1 )))
+ start = 1
+ end = length(data)
+ if (length(data[[1]]$centered_serie) < length(data[[2]]$centered_serie))
+ start = 2
+ if (length(data[[length(data)]]$centered_serie) <
+ length(data[[length(data)-1]]$centered_serie))
{
- data$removeLast()
+ end = end-1
}
-
+ if (start>1 || end<length(data))
+ data = data[start:end]
data
}
L = length(ref_serie)
# Determine indices of no-NAs days followed by no-NAs tomorrows
- first_day = ifelse(length(data$getCenteredSerie(1))<L, 2, 1)
- fdays_indices = c()
- for (i in first_day:(index-1))
+ fdays = c()
+ for (i in 1:(index-1))
{
if ( !any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))) )
- fdays_indices = c(fdays_indices, i)
+ fdays = c(fdays, i)
}
- distances = sapply(fdays_indices, function(i) {
+ distances = sapply(fdays, function(i) {
sqrt( sum( (ref_serie - data$getCenteredSerie(i))^2 ) / L )
})
indices = sort(distances, index.return=TRUE)$ix[1:min(limit,length(distances))]
yrange = quantile( c(ref_serie, sapply( indices, function(i) {
- ii = fdays_indices[i]
- serie = c(data$getCenteredSerie(ii), data$getCenteredSerie(ii+1))
+ serie = c(data$getCenteredSerie(fdays[i]), data$getCenteredSerie(fdays[i]+1))
if (!all(is.na(serie)))
return (range(serie, na.rm=TRUE))
c()
par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2)
for ( i in seq_len(length(indices)+1) )
{
- ii = ifelse(i<=length(indices), fdays_indices[ indices[plot_order[i]] ], index)
+ ii = ifelse(i<=length(indices), fdays[ indices[plot_order[i]] ], index)
plot(c(data$getCenteredSerie(ii),data$getCenteredSerie(ii+1)),
ylim=yrange, type="l", col=colors[i],
xlab=ifelse(i==1,"Temps (en heures)",""), ylab=ifelse(i==1,"PM10 centré",""))
}
abline(v=24, lty=2, col=colors()[56])
}
- list("indices"=c(fdays_indices[ indices[plot_order] ],index), "colors"=colors)
+ list("indices"=c(fdays[ indices[plot_order] ],index), "colors"=colors)
}
#' @title Plot similarities
if (!requireNamespace("rainbow", quietly=TRUE))
stop("Functional boxplot requires the rainbow package")
- start_index = 1
- end_index = data$getSize()
- if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2)))
- {
- # Shifted start (7am, or 1pm, or...)
- start_index = 2
- end_index = data$getSize() - 1
- }
-
L = length(data$getCenteredSerie(2))
- series_matrix = sapply(start_index:end_index, function(index) {
+ series_matrix = sapply(1:data$getSize(), function(index) {
if (filter(index))
as.matrix(data$getSerie(index))
else
ref_var = apply(ref_series, 2, sd)
# Determine indices of no-NAs days followed by no-NAs tomorrows
- first_day = ifelse(length(data$getCenteredSerie(1))<length(ref_series[1,]), 2, 1)
- fdays_indices = c()
- for (i in first_day:(tail(indices,1)-1))
+ fdays = c()
+ for (i in 1:(tail(indices,1)-1))
{
if ( !any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))) )
- fdays_indices = c(fdays_indices, i)
+ fdays = c(fdays, i)
}
# TODO: 3 == magic number
random_var = matrix(nrow=3, ncol=48)
for (mc in seq_len(nrow(random_var)))
{
- random_indices = sample(fdays_indices, length(indices))
+ random_indices = sample(fdays, length(indices))
random_series = t( sapply(random_indices, function(i) {
c( data$getSerie(i), data$getSerie(i+1) )
}) )
-#' @title dateIndexToInteger
+#' dateIndexToInteger
#'
-#' @description Transform a (potential) date index into an integer (relative to data)
+#' Transform a (potential) date index into an integer (relative to data)
#'
#' @param index Date (or integer) index
#' @param data Object of class \code{Data}
integerIndex <- round( as.numeric(
difftime(indexAsDate, as.Date(data$getTime(1)[1])) ) ) + 1
if (integerIndex >= 1 && integerIndex <= data$getSize())
- {
- #WARNING: if index>=2 && series start at date >0h, result must be shifted
- #Convention: if first date asked, return first matching index (i.e.: 1)
- shift = 0
- if (integerIndex >= 2)
- {
- date1 = as.POSIXlt(data$getTime(1)[1])
- date2 = as.POSIXlt(data$getTime(2)[1])
- shift = (date1$year==date2$year && date1$mon==date2$mon && date1$mday==date2$mday)
- }
- return (integerIndex + shift)
- }
+ return (integerIndex)
stop("Date outside data range")
}
stop("Unrecognized index format")
}
-#' @title integerIndexToDate
+#' integerIndexToDate
#'
-#' @description Transform an integer index to date index (relative to data)
+#' Transform an integer index to date index (relative to data)
#'
#' @param index Date (or integer) index
#' @param data Object of class \code{Data}
as.Date( data$getTime(index)[1] )
}
-#' @title getSimilarDaysIndices
+#' getSimilarDaysIndices
#'
-#' @description Find similar days indices in the past
+#' Find similar days indices in the past
#'
#' @param index Day index (numeric or date)
#' @param limit Maximum number of indices to return
return ( days[1:min(limit,length(days))] )
}
+
+#' getSerie
+#'
+#' Return a time-serie from its centered version + level
+#'
+#' @param data A list as returned by \code{getData}
+#' @param index The index to return
+#'
+#' @export
+getSerie = function(data, index)
+ data[[index]]$centered_serie + data[[index]]$level
"cells": [
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 1,
"metadata": {
"collapsed": false
},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "Loading required package: R6\n"
+ ]
+ }
+ ],
"source": [
"library(talweg)"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 2,
"metadata": {
"collapsed": false
},
- "outputs": [],
+ "outputs": [
+ {
+ "ename": "ERROR",
+ "evalue": "Error in private$.data <<- list(): object 'private' not found\n",
+ "output_type": "error",
+ "traceback": [
+ "Error in private$.data <<- list(): object 'private' not found\nTraceback:\n",
+ "1. getData(ts_data, exo_data, input_tz = \"Europe/Paris\", working_tz = \"Europe/Paris\", \n . predict_at = 7)",
+ "2. Data$new()",
+ "3. .subset2(public_bind_env, \"initialize\")(...)",
+ "4. initializeData(self, private)"
+ ]
+ }
+ ],
"source": [
"ts_data = read.csv(system.file(\"extdata\",\"pm10_mesures_H_loc.csv\",package=\"talweg\"))\n",
"exo_data = read.csv(system.file(\"extdata\",\"meteo_extra_noNAs.csv\",package=\"talweg\"))\n",
projects / talweg.git / commitdiff