From 09cf9c19b5c6a04bc23c58b7ac8a4bbae2c6827d Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
Date: Mon, 13 Feb 2017 18:15:20 +0100
Subject: [PATCH] fixes and improvements

---
 R/D_Neighbors.R                               |  12 +-
 R/D_Persistence.R                             |  17 +-
 R/Data.R                                      |   6 +
 R/Forecast.R                                  |   2 +-
 R/S_Average.R                                 |  24 +-
 R/S_Neighbors.R                               |  30 +-
 R/S_Persistence.R                             |  23 +-
 R/getData.R                                   |  39 +--
 R/getForecast.R                               |   7 +-
 R/plot.R                                      |  19 +-
 R/utils.R                                     |  36 ++-
 ...t_13-01-2017.rnw => report_2017-01-13.rnw} |   0
 ...t_02-02-2017.Rnw => report_2017-02-02.Rnw} |   0
 ...-02-2017.ipynb => report_2017-02-02.ipynb} |  13 +-
 reports/report_2017-03-01.ipynb               | 281 ++++++++++++++++++
 15 files changed, 432 insertions(+), 77 deletions(-)
 rename reports/{report_13-01-2017.rnw => report_2017-01-13.rnw} (100%)
 rename reports/{report_02-02-2017.Rnw => report_2017-02-02.Rnw} (100%)
 rename reports/{report_02-02-2017.ipynb => report_2017-02-02.ipynb} (98%)
 create mode 100644 reports/report_2017-03-01.ipynb

diff --git a/R/D_Neighbors.R b/R/D_Neighbors.R
index f7bb5a8..dba5f37 100644
--- a/R/D_Neighbors.R
+++ b/R/D_Neighbors.R
@@ -4,14 +4,18 @@
 #' @inheritParams getZeroDeltaForecast
 getNeighborsDeltaForecast = function(data, today, memory, horizon, shape_params, ...)
 {
-	if (any(is.na(shape_params$weights) | is.na(shape_params$indices)))
+	first_day = max(1, today-memory)
+	filter = shape_params$indices >= first_day
+	indices = shape_params$indices[filter]
+	weights = shape_params$weights[filter]
+	if (any(is.na(weights) | is.na(indices)))
 		return (NA)
 
-	gaps = sapply(shape_params$indices, function(i) {
+	gaps = sapply(indices, function(i) {
 		data$getSerie(i+1)[1] - tail(data$getSerie(i), 1)
 	})
 
-	scal_product = shape_params$weights * gaps
-	norm_fact = sum( shape_params$weights[!is.na(scal_product)] )
+	scal_product = weights * gaps
+	norm_fact = sum( weights[!is.na(scal_product)] )
 	sum(scal_product, na.rm=TRUE) / norm_fact
 }
diff --git a/R/D_Persistence.R b/R/D_Persistence.R
index 2f7935a..979bf05 100644
--- a/R/D_Persistence.R
+++ b/R/D_Persistence.R
@@ -4,6 +4,19 @@
 #' @inheritParams getZeroDeltaForecast
 getPersistenceDeltaForecast = function(data, today, memory, horizon, shape_params, ...)
 {
-	# Return level of last similar "tomorrow" (thus -6) computed on "horizon"
-	data$getSerie(today-6)[1] - tail(data$getSerie(today-7), 1)
+	#return gap between end of similar day curve and first day of tomorrow (in the past)
+	first_day = max(1, today-memory)
+	index = today-7
+	repeat
+	{
+		{
+			last_similar_serie_end = tail( data$getCenteredSerie(index), 1)
+			last_similar_tomorrow_begin = data$getCenteredSerie(index+1)[1]
+			index = index - 7
+		};
+		if (!is.na(last_similar_serie_end) && !is.na(last_similar_tomorrow_begin))
+			return (last_similar_tomorrow_begin - last_similar_serie_end);
+		if (index < first_day)
+			return (NA)
+	}
 }
diff --git a/R/Data.R b/R/Data.R
index 311453b..a17e262 100644
--- a/R/Data.R
+++ b/R/Data.R
@@ -42,36 +42,42 @@ Data = setRefClass(
 		{
 			"Get time values at specified index"
 
+			index = dateIndexToInteger(index, .self)
 			data[[index]]$time
 		},
 		getCenteredSerie = function(index)
 		{
 			"Get serie values (centered) at specified index"
 
+			index = dateIndexToInteger(index, .self)
 			data[[index]]$serie
 		},
 		getLevel = function(index)
 		{
 			"Get level for the serie at specified index"
 
+			index = dateIndexToInteger(index, .self)
 			data[[index]]$level
 		},
 		getSerie = function(index)
 		{
 			"Get serie values (centered+level) at specified index"
 
+			index = dateIndexToInteger(index, .self)
 			data[[index]]$serie + data[[index]]$level
 		},
 		getExoHat = function(index)
 		{
 			"Get exogeous predictions at specified index"
 
+			index = dateIndexToInteger(index, .self)
 			data[[index]]$exo_hat
 		},
 		getExoDm1 = function(index)
 		{
 			"Get exogenous measures the day before specified index"
 
+			index = dateIndexToInteger(index, .self)
 			data[[index]]$exo_Dm1
 		}
 	)
diff --git a/R/Forecast.R b/R/Forecast.R
index 1a3505b..57817e7 100644
--- a/R/Forecast.R
+++ b/R/Forecast.R
@@ -49,7 +49,7 @@ Forecast = setRefClass(
 		},
 		getIndexInData = function(index)
 		{
-			"Get (day)index at specified index"
+			"Get (day) index in data where prediction took place"
 
 			pred[[index]]$index
 		}
diff --git a/R/S_Average.R b/R/S_Average.R
index 9694693..819531d 100644
--- a/R/S_Average.R
+++ b/R/S_Average.R
@@ -2,8 +2,8 @@
 #'
 #' @title Average Shape Forecaster
 #'
-#' @description Return the (pointwise) average of the all the centered day curves in the past.
-#'   Inherits \code{\link{ShapeForecaster}}
+#' @description Return the (pointwise) average of the all the (similar) centered day curves
+#'   in the past. Inherits \code{\link{ShapeForecaster}}
 AverageShapeForecaster = setRefClass(
 	Class = "AverageShapeForecaster",
 	contains = "ShapeForecaster",
@@ -16,12 +16,24 @@ AverageShapeForecaster = setRefClass(
 		predict = function(today, memory, horizon, ...)
 		{
 			avg = rep(0., horizon)
-			for (i in (today-memory):today)
+			first_day = max(1, today-memory)
+			index = today-7 + 1
+			nb_no_na_series = 0
+			repeat
 			{
-				if (!any(is.na(data$getSerie(i))))
-					avg = avg + data$getCenteredSerie(i)
+				{
+					serie_on_horizon = data$getCenteredSerie(index)[1:horizon]
+					index = index - 7
+				};
+				if (!any(is.na(serie_on_horizon)))
+				{
+					avg = avg + serie_on_horizon
+					nb_no_na_series = nb_no_na_series + 1
+				};
+				if (index < first_day)
+					break
 			}
-			avg
+			avg / nb_no_na_series
 		}
 	)
 )
diff --git a/R/S_Neighbors.R b/R/S_Neighbors.R
index c44bd9b..b8e32cc 100644
--- a/R/S_Neighbors.R
+++ b/R/S_Neighbors.R
@@ -23,9 +23,33 @@ NeighborsShapeForecaster = setRefClass(
 			if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2)))
 				first_day = 2
 
-			# Predict only on non-NAs days (TODO:)
-			if (any(is.na(data$getSerie(today))))
-				return (NA)
+			# Predict only on (almost) non-NAs days
+			nas_in_serie = is.na(data$getSerie(today))
+			if (any(nas_in_serie))
+			{
+				#TODO: better define "repairing" conditions (and method)
+				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
 			fdays_indices = c()
diff --git a/R/S_Persistence.R b/R/S_Persistence.R
index 017402f..08647e3 100644
--- a/R/S_Persistence.R
+++ b/R/S_Persistence.R
@@ -2,8 +2,8 @@
 #'
 #' @title Persistence Shape Forecaster
 #'
-#' @description Return the last centered last (similar) day curve (may be a lot of NAs,
-#'   but we cannot do better). Inherits \code{\link{ShapeForecaster}}
+#' @description Return the last centered last (similar) day curve.
+#'   Inherits \code{\link{ShapeForecaster}}
 PersistenceShapeForecaster = setRefClass(
 	Class = "PersistenceShapeForecaster",
 	contains = "ShapeForecaster",
@@ -15,11 +15,20 @@ PersistenceShapeForecaster = setRefClass(
 		},
 		predict = function(today, memory, horizon, ...)
 		{
-			#return centered last (similar) day curve (may be a lot of NAs, but we cannot do better)
-			last_similar_serie = data$getCenteredSerie(today-6)[1:horizon]
-			if (any(is.na(last_similar_serie))) #TODO:
-				return (NA)
-			last_similar_serie
+			#return centered last (similar) day curve, avoiding NAs until memory is run
+			first_day = max(1, today-memory)
+			index = today-7 + 1
+			repeat
+			{
+				{
+					last_similar_serie = data$getCenteredSerie(index)[1:horizon]
+					index = index - 7
+				};
+				if (!any(is.na(last_similar_serie)))
+					return (last_similar_serie);
+				if (index < first_day)
+					return (NA)
+			}
 		}
 	)
 )
diff --git a/R/getData.R b/R/getData.R
index c2e6842..5139a4d 100644
--- a/R/getData.R
+++ b/R/getData.R
@@ -13,14 +13,13 @@
 #' @param date_format How date/time are stored (e.g. year/month/day hour:minutes;
 #'   see \code{strptime})
 #' @param working_tz Timezone to work with ("GMT" or e.g. "Europe/Paris")
-#' @param predict_at When does the prediction take place ? ab[:cd][:ef] where a,b,c,d,e,f
-#'   in (0,9) and define an hour[minute[second]]; time must be present in the file
+#' @param predict_at When does the prediction take place ? Integer, in hours. Default: 0
 #'
 #' @return An object of class Data
 #'
 #' @export
 getData = function(ts_data, exo_data,
-	input_tz="GMT", date_format="%d/%m/%Y %H:%M", working_tz="GMT", predict_at="00")
+	input_tz="GMT", date_format="%d/%m/%Y %H:%M", working_tz="GMT", predict_at=0)
 {
 	# Sanity checks (not full, but sufficient at this stage)
 	if (!is.character(input_tz) || !is.character(working_tz))
@@ -31,10 +30,9 @@ getData = function(ts_data, exo_data,
 		stop("Bad time-series input (data frame or CSV file)")
 	if (is.character(ts_data))
 		ts_data = ts_data[1]
-	pattern_index_in_predict_at = grep("^[0-9]{2}(:[0-9]{2}){0,2}$", predict_at)
-	if (!is.character(predict_at) || length(pattern_index_in_predict_at) == 0)
-		stop("Bad predict_at ( ^[0-9]{2}(:[0-9]{2}){0,2}$ )")
-	predict_at = predict_at[ pattern_index_in_predict_at[1] ]
+	predict_at = as.integer(predict_at)[1]
+	if (predict_at<0 || predict_at>23)
+		stop("Bad predict_at (0-23)")
 	if (!is.character(date_format))
 		stop("Bad date_format (character)")
 	date_format = date_format[1]
@@ -55,11 +53,6 @@ getData = function(ts_data, exo_data,
 	formatted_dates_POSIXlt = strptime(as.character(ts_df[,1]), date_format, tz=input_tz)
 	ts_df[,1] = format(as.POSIXct(formatted_dates_POSIXlt), tz=working_tz, usetz=TRUE)
 
-	if (nchar(predict_at) == 2)
-		predict_at = paste(predict_at,":00",sep="")
-	if (nchar(predict_at) == 5)
-		predict_at = paste(predict_at,":00",sep="")
-
 	line = 1 #index in PM10 file (24 lines for 1 cell)
 	nb_lines = nrow(ts_df)
 	nb_exos = ( ncol(exo_df) - 1 ) / 2
@@ -69,22 +62,20 @@ getData = function(ts_data, exo_data,
 	{
 		time = c()
 		serie = c()
-		repeat {
-		{
-			time = c(time, ts_df[line,1])
-			serie = c(serie, ts_df[line,2])
-			line = line + 1
-		};
-		if (line >= nb_lines + 1
-			# NOTE: always second part of date/time, because it has been formatted
-			|| strsplit(as.character(ts_df[line-1,1])," ")[[1]][2] == predict_at)
+		repeat
 		{
-			break
-		}}
+			{
+				time = c(time, ts_df[line,1])
+				serie = c(serie, ts_df[line,2])
+				line = line + 1
+			};
+			if (line >= nb_lines + 1 || as.POSIXlt(ts_df[line-1,1])$hour == predict_at)
+				break
+		}
 
 		# NOTE: if predict_at does not cut days at midnight,
 		# for the exogenous to be synchronized they need to be shifted
-		if (predict_at != "00:00:00")
+		if (predict_at > 0)
 		{
 			exo_hat = as.data.frame(exo_df[max(1,i-1),(1+nb_exos+1):(1+2*nb_exos)])
 			exo_Dm1 = if (i>=3) as.data.frame(exo_df[i-1,2:(1+nb_exos)]) else NA
diff --git a/R/getForecast.R b/R/getForecast.R
index c3248a9..e126946 100644
--- a/R/getForecast.R
+++ b/R/getForecast.R
@@ -42,7 +42,7 @@ getForecast = function(data, indices, memory, horizon,
 	horizon = as.integer(horizon)[1]
 	if (horizon<=0 || horizon>length(data$getCenteredSerie(2)))
 		stop("Horizon too short or too long")
-	indices = as.integer(indices)
+	indices = sapply( seq_along(indices), function(i) dateIndexToInteger(indices[i], data) )
 	if (any(indices<=0 | indices>data$getSize()))
 		stop("Indices out of range")
 	indices = sapply(indices, dateIndexToInteger, data)
@@ -56,9 +56,6 @@ getForecast = function(data, indices, memory, horizon,
 	pred = list()
 	for (today in indices)
 	{
-		#NOTE: To estimate timing...
-#		print(paste("Predict until index",today))
-
 		#shape always predicted first (on centered series, no scaling taken into account),
 		#with side-effect: optimize some parameters (h, weights, ...)
 		predicted_shape = shape_forecaster$predict(today, memory, horizon, ...)
@@ -69,7 +66,7 @@ getForecast = function(data, indices, memory, horizon,
 		#TODO: this way is faster than a call to append(); why ?
 		pred[[length(pred)+1]] = list(
 			# Predict shape and align it on end of current day
-			serie = predicted_shape + tail( data$getSerie(today), 1 ) - predicted_shape[1],
+			serie = predicted_shape + tail( data$getSerie(today), 1 ) - predicted_shape[1] +
 				predicted_delta, #add predicted jump
 			params = shape_forecaster$getParameters(),
 			index = today )
diff --git a/R/plot.R b/R/plot.R
index e5d4753..a7c9395 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -13,9 +13,9 @@ plotPredReal <- function(data, pred, index)
 	par(mar=c(4.7,5,1,1), cex.axis=2, cex.lab=2, lwd=2)
 	measure = data$getSerie(pred$getIndexInData(index)+1)[1:horizon]
 	yrange = range( pred$getSerie(index), measure )
-	plot(measure, type="l", ylim=yrange, lwd=3)
+	plot(measure, type="l", ylim=yrange, lwd=3, xlab="Temps (en heures)", ylab="PM10")
 	par(new=TRUE)
-	plot(pred$getSerie(index), type="l", col=2, ylim=yrange, lwd=3)
+	plot(pred$getSerie(index), type="l", col="#0000FF", ylim=yrange, lwd=3, xlab="", ylab="")
 }
 
 #' @title Plot filaments
@@ -85,20 +85,23 @@ plotSimils <- function(pred, index)
 #' @description Draw error graphs, potentially from several runs of \code{getForecast}
 #'
 #' @param err Error as returned by \code{getError}
+#' @param cols Colors for each error (default: 1,2,3,...)
 #'
 #' @seealso \code{\link{plotPredReal}}, \code{\link{plotFilaments}}, \code{\link{plotSimils}}
 #'   \code{\link{plotFbox}}
 #'
 #' @export
-plotError <- function(err)
+plotError <- function(err, cols=seq_along(err))
 {
+	if (!is.null(err$abs))
+		err = list(err)
 	par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=2, cex.lab=2, lwd=2)
 	L = length(err)
 	yrange = range( sapply(1:L, function(index) ( err[[index]]$abs$day ) ), na.rm=TRUE )
 	for (i in seq_len(L))
 	{
 		plot(err[[i]]$abs$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""),
-			ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=i)
+			ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i])
 		if (i < L)
 			par(new=TRUE)
 	}
@@ -106,7 +109,7 @@ plotError <- function(err)
 	for (i in seq_len(L))
 	{
 		plot(err[[i]]$abs$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""),
-			ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=i)
+			ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i])
 		if (i < L)
 			par(new=TRUE)
 	}
@@ -114,7 +117,7 @@ plotError <- function(err)
 	for (i in seq_len(L))
 	{
 		plot(err[[i]]$MAPE$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""),
-			ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=i)
+			ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i])
 		if (i < L)
 			par(new=TRUE)
 	}
@@ -122,7 +125,7 @@ plotError <- function(err)
 	for (i in seq_len(L))
 	{
 		plot(err[[i]]$MAPE$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""),
-			ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=i)
+			ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i])
 		if (i < L)
 			par(new=TRUE)
 	}
@@ -136,7 +139,7 @@ plotError <- function(err)
 #' @param fiter Optional filter: return TRUE on indices to process
 #'
 #' @export
-plotFbox <- function(data, filter=function(index) (TRUE))
+plotFbox <- function(data, filter=function(index) TRUE)
 {
 	if (!requireNamespace("rainbow", quietly=TRUE))
 		stop("Functional boxplot requires the rainbow package")
diff --git a/R/utils.R b/R/utils.R
index 879a7d7..3777f2d 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -8,23 +8,47 @@
 #' @export
 dateIndexToInteger = function(index, data)
 {
+	index = index[1]
 	if (is.numeric(index))
 		index = as.integer(index)
-	#Undocumented "private" method to translate index --> date (is needed)
 	if (is.integer(index))
-		return (index[1])
-	if (is.character(index))
+		return (index)
+	if (inherits(index, "Date") || is.character(index))
 	{
-		tryCatch(dt <- as.POSIXct(index[1]), finally=print("Unrecognized index format"))
+		tryCatch(dt <- as.POSIXct(index), error=function(e) stop("Unrecognized index format"))
 		#TODO: tz arg to difftime ?
-		integerIndex <- as.integer( difftime(dt, data$getTime(1)) ) + 1
+		integerIndex <- round( (as.numeric( difftime(dt, data$getTime(1)) ))[1] ) + 1
 		if (integerIndex > 0 && integerIndex <= data$getSize())
-			return (integerIndex)
+		{
+			#WARNING: if series start at date >0h, result must be shifted
+			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 + ifelse(shift,1,0))
+		}
 		stop("Date outside data range")
 	}
 	stop("Unrecognized index format")
 }
 
+#' @title integerIndexToDate
+#'
+#' @description Transform an integer index to date index (relative to data)
+#'
+#' @param index Date (or integer) index
+#' @param data Object of class \code{Data}
+#'
+#' @export
+integerIndexToDate = function(index, data)
+{
+	index = index[1]
+	if (is.numeric(index))
+		index = as.integer(index)
+	if (!is.integer(index))
+		stop("'index' should be an integer")
+	as.Date( data$getTime(index)[1] )
+}
+
 #' @title getSimilarDaysIndices
 #'
 #' @description Find similar days indices in the past
diff --git a/reports/report_13-01-2017.rnw b/reports/report_2017-01-13.rnw
similarity index 100%
rename from reports/report_13-01-2017.rnw
rename to reports/report_2017-01-13.rnw
diff --git a/reports/report_02-02-2017.Rnw b/reports/report_2017-02-02.Rnw
similarity index 100%
rename from reports/report_02-02-2017.Rnw
rename to reports/report_2017-02-02.Rnw
diff --git a/reports/report_02-02-2017.ipynb b/reports/report_2017-02-02.ipynb
similarity index 98%
rename from reports/report_02-02-2017.ipynb
rename to reports/report_2017-02-02.ipynb
index 4f49077..00380f7 100644
--- a/reports/report_02-02-2017.ipynb
+++ b/reports/report_2017-02-02.ipynb
@@ -57,20 +57,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {
     "collapsed": false
    },
-   "outputs": [
-    {
-     "ename": "ERROR",
-     "evalue": "Error in eval(expr, envir, enclos): could not find function \"getData\"\n",
-     "output_type": "error",
-     "traceback": [
-      "Error in eval(expr, envir, enclos): could not find function \"getData\"\nTraceback:\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "# Voir ?getData pour les arguments\n",
     "data = getData(ts_data=\"data/pm10_mesures_H_loc.csv\", exo_data=\"data/meteo_extra_noNAs.csv\",\n",
diff --git a/reports/report_2017-03-01.ipynb b/reports/report_2017-03-01.ipynb
new file mode 100644
index 0000000..4eecb50
--- /dev/null
+++ b/reports/report_2017-03-01.ipynb
@@ -0,0 +1,281 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "library(talweg)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "data = getData(ts_data=\"../data/pm10_mesures_H_loc.csv\", exo_data=\"../data/meteo_extra_noNAs.csv\",\n",
+    "               input_tz = \"Europe/Paris\", working_tz=\"Europe/Paris\", predict_at=7)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Pollution par chauffage"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "p_ch_nn = getForecast(data, seq(as.Date(\"2015-01-18\"),as.Date(\"2015-01-24\"),\"days\"), Inf, 17,\n",
+    "                   \"Neighbors\", \"Neighbors\", simtype=\"mix\", same_season=FALSE, mix_strategy=\"mult\")\n",
+    "p_ch_pz = getForecast(data, seq(as.Date(\"2015-01-18\"),as.Date(\"2015-01-24\"),\"days\"), Inf, 17,\n",
+    "                   \"Persistence\", \"Zero\")\n",
+    "p_ch_az = getForecast(data, seq(as.Date(\"2015-01-18\"),as.Date(\"2015-01-24\"),\"days\"), Inf, 17,\n",
+    "                   \"Average\", \"Zero\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "e_ch_nn = getError(data, p_ch_nn, 17)\n",
+    "e_ch_pz = getError(data, p_ch_pz, 17)\n",
+    "e_ch_az = getError(data, p_ch_az, 17)\n",
+    "options(repr.plot.width=9, repr.plot.height=6)\n",
+    "plotError(list(e_ch_nn, e_ch_pz, e_ch_az), cols=c(1,2,colors()[258]))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "par(mfrow=c(1,2))\n",
+    "options(repr.plot.width=9, repr.plot.height=4)\n",
+    "plotPredReal(data, p_ch_nn, 3)\n",
+    "plotPredReal(data, p_ch_nn, 4)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "par(mfrow=c(1,2))\n",
+    "plotFilaments(data, p_ch_nn$getIndexInData(3))\n",
+    "plotFilaments(data, p_ch_nn$getIndexInData(4))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "par(mfrow=c(1,3))\n",
+    "plotSimils(p_ch_nn, 3)\n",
+    "plotSimils(p_ch_nn, 4)\n",
+    "plotSimils(p_ch_nn, 5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Pollution par épandage"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "p_ep_nn = getForecast(data, seq(as.Date(\"2015-03-15\"),as.Date(\"2015-03-21\"),\"days\"), Inf, 17,\n",
+    "                   \"Neighbors\", \"Neighbors\", simtype=\"mix\", same_season=FALSE, mix_strategy=\"mult\")\n",
+    "p_ep_pz = getForecast(data, seq(as.Date(\"2015-03-15\"),as.Date(\"2015-03-21\"),\"days\"), Inf, 17,\n",
+    "                   \"Persistence\", \"Zero\")\n",
+    "p_ep_az = getForecast(data, seq(as.Date(\"2015-03-15\"),as.Date(\"2015-03-21\"),\"days\"), Inf, 17,\n",
+    "                   \"Average\", \"Zero\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "e_ep_nn = getError(data, p_ep_nn, 17)\n",
+    "e_ep_pz = getError(data, p_ep_pz, 17)\n",
+    "e_ep_az = getError(data, p_ep_az, 17)\n",
+    "options(repr.plot.width=9, repr.plot.height=6)\n",
+    "plotError(list(e_ep_nn, e_ep_pz, e_ep_az), cols=c(1,2,colors()[258]))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "par(mfrow=c(1,2))\n",
+    "options(repr.plot.width=9, repr.plot.height=4)\n",
+    "plotPredReal(data, p_ep_nn, 3)\n",
+    "plotPredReal(data, p_ep_nn, 4)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "par(mfrow=c(1,2))\n",
+    "plotFilaments(data, p_ep_nn$getIndexInData(3))\n",
+    "plotFilaments(data, p_ep_nn$getIndexInData(4))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "par(mfrow=c(1,3))\n",
+    "plotSimils(p_ep_nn, 3)\n",
+    "plotSimils(p_ep_nn, 4)\n",
+    "plotSimils(p_ep_nn, 5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Semaine non polluée"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "p_np_nn = getForecast(data, seq(as.Date(\"2015-04-26\"),as.Date(\"2015-05-02\"),\"days\"), Inf, 17,\n",
+    "                   \"Neighbors\", \"Neighbors\", simtype=\"mix\", same_season=FALSE, mix_strategy=\"mult\")\n",
+    "p_np_pz = getForecast(data, seq(as.Date(\"2015-04-26\"),as.Date(\"2015-05-02\"),\"days\"), Inf, 17,\n",
+    "                   \"Persistence\", \"Zero\")\n",
+    "p_np_az = getForecast(data, seq(as.Date(\"2015-04-26\"),as.Date(\"2015-05-02\"),\"days\"), Inf, 17,\n",
+    "                   \"Average\", \"Zero\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "e_np_nn = getError(data, p_np_nn, 17)\n",
+    "e_np_pz = getError(data, p_np_pz, 17)\n",
+    "e_np_az = getError(data, p_np_az, 17)\n",
+    "options(repr.plot.width=9, repr.plot.height=6)\n",
+    "plotError(list(e_np_nn, e_np_pz, e_np_az), cols=c(1,2,colors()[258]))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "par(mfrow=c(1,2))\n",
+    "options(repr.plot.width=9, repr.plot.height=4)\n",
+    "plotPredReal(data, p_np_nn, 3)\n",
+    "plotPredReal(data, p_np_nn, 4)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "par(mfrow=c(1,2))\n",
+    "plotFilaments(data, p_np_nn$getIndexInData(3))\n",
+    "plotFilaments(data, p_np_nn$getIndexInData(4))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "par(mfrow=c(1,3))\n",
+    "plotSimils(p_np_nn, 3)\n",
+    "plotSimils(p_np_nn, 4)\n",
+    "plotSimils(p_np_nn, 5)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "R",
+   "language": "R",
+   "name": "ir"
+  },
+  "language_info": {
+   "codemirror_mode": "r",
+   "file_extension": ".r",
+   "mimetype": "text/x-r-source",
+   "name": "R",
+   "pygments_lexer": "r",
+   "version": "3.3.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
-- 
2.44.0