From 9dbd32d9f73bbc91c65c2c6a05095f084c6fdb16 Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
Date: Tue, 25 Apr 2017 23:15:41 +0200
Subject: [PATCH] commented a bit Bruno's file

---
 Etude_En_Cours_R_E_bis.R | 171 +++++++++++++++++++++------------------
 1 file changed, 92 insertions(+), 79 deletions(-)

diff --git a/Etude_En_Cours_R_E_bis.R b/Etude_En_Cours_R_E_bis.R
index 69b1da0..40c7954 100644
--- a/Etude_En_Cours_R_E_bis.R
+++ b/Etude_En_Cours_R_E_bis.R
@@ -1,29 +1,37 @@
-setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno")
+#setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno")
 rm(list=ls())
 
-# Lecture des données
-pm  =  read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=";",dec=".",header=T) #,row.names=1)
-n     = dim(pm)[1]
-datedebut = "10/12/2008"
+# Lecture des données: pm = dataframe 2 colonnes, date-time puis PM10 horaire
+pm  =  read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=";",dec=".",header=T)
+#n = dim(pm)[1]
+#datedebut = "10/12/2008"
 
-# Chargement des données météo et indicateurs
+# Chargement des données météo et indicateurs: VarExp matrice des données météo,
+# première colonne = date, première rangée = second jour
 VarExp <- read.table("DATA/meteo_extra_jourMois.csv",sep=";",dec=".",header=T)
 VarExp <- VarExp[-1,]
 
-# Lecture des dates
-dates =  read.table("DATA/Dates_jours.csv",header=F,as.is=T)
+# Lecture des dates: ??? sans doute liste des jours de 10/12/2008 à 06/2016...
+dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T)
 dates = dates[,1]
 
-# Contruction des matrices de données
+# Contruction des matrices de données: pm.h = matrice des séries en ligne
 pm.h   <- matrix(pm[,2],ncol=24,byrow=TRUE)
-
 Nlignes = nrow(pm.h)
+
+# Data = matrice des couples de jours (séries de 48 PM10),
+# première ligne = jour 1/2, dernière = jour N-1/N
 Data = cbind(pm.h[1:(Nlignes -1), ], pm.h[2:Nlignes, ])
+
+# dates2 = dates du 2eme jour au dernier
 dates2 = dates[2:Nlignes]
+
 rownames(Data) = dates2
 # df contient l'ensemble des données.
 #df <- cbind(Data,varexp[,-1])
 df <- Data
+
+# Complétion des variables exogènes (du 2eme jour) par les PM10 moyen ce même jour
 PMjour <- apply(df[,25:48],1,mean,na.rm=T)
 dfexp <- cbind(VarExp,PMjour)
 
@@ -46,79 +54,84 @@ Err24 = NULL
 ErrPrev = NULL
 Kvois = NULL
 
-for (Hc in 5:24) {
-
-H=24+Hc
-fen=H
-#fen=12 ou 24 ou H
-L = (H-fen+1):H
-
-# Premier conditionnement : mois
-indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4 | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10
-data = df[indcond,]
-varexp = dfexp[indcond,]
-
-nl = (1:nrow(data))[rownames(data)==Dates[j]]
-dateJPrev = rownames(data)[nl+ij]
-dataj = as.numeric(data[nl+ij, 1:48])
-data = data[-(nl + ij), ]
-varexp = varexp[-(nl + ij), ]
-indNA = attr(na.omit(data[, 1:48]),"na.action")
-data = data[-indNA,]
-varexp = varexp[-indNA,]
-
-# Conditionnement : les jours avec PMjour +/- large
-large = 1
-bornes = mean(dataj[25:48])+c(-large,large)
-indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
-data = data[indcond,]
-varexp = varexp[indcond,]
-
-D = rep(0,nrow(data))
-for (k in 1:nrow(data)) {
-       #D[k] = sqrt(sum((1:H)*(dataj[L] -data[k,L])^2))
-       D[k] = sqrt(sum((dataj[L] -data[k,L])^2))
-}
-ind = order(D)[1:nbvois]
-w = 1/(D[ind]^2)
-w = w/sum(w)
-W = w%o%rep(1,48)
-JourMoy = apply(data[ind, 1:48], 2, mean)
-#JourMoy = apply(W*data[ind, 1:48], 2, sum)
-NomFile = paste("Voisins_Epandage_PMjour_Hc_",Hc,".png",sep="")
-Titre = paste("Jour à prévoir : ",dateJPrev," - ", length(ind)," voisins",sep="")
-#erreur = sqrt(sum((dataj[25:48] - JourMoy[25:48])^2))
-if(Hc==24){erreurPrev=NA}else{erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48]))}
-erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48]))
-png(NomFile)
-matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind), 
-	cex.axis=1.4, cex.main = 1.7, cex.lab=1.5, 
-	xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1)," / ",round(erreurPrev,1)), main=Titre)
-legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind))
-lines(1:48, dataj, lwd=2.5)
-lines(JourMoy, lty = 2, lwd=2)
-abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2)
-xx=dev.off()
-
-
-
-  Err24 = c(Err24, erreur24)
-  ErrPrev = c(ErrPrev, erreurPrev)
-ResDates = cbind(ResDates, rownames(data)[ind])
+for (Hc in 5:24)
+{
+	H=24+Hc #Hc: dernier PM10 connu avant prédiction des 24-Hc restants, le jour 2
+	L = 1:H
+
+	# Premier conditionnement : mois
+	indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4
+		| dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10
+	data = df[indcond,] #restriction aux couples dont le mois du 2eme jour est 2,3,4,9 ou 10
+	varexp = dfexp[indcond,] #de même sur les exogènes (+ PM10 moyens)
+
+	nl = (1:nrow(data))[rownames(data)==Dates[j]] #numéro de ligne où 2eme jour == 16/03/2015
+	dateJPrev = rownames(data)[nl+ij] # ??? la date du 1er jour du couple de jours pile une semaine après le couple 'nl'
+	dataj = as.numeric(data[nl+ij, 1:48]) #extraction de cette série dans dataj
+	data = data[-(nl + ij), ] #suppression de cette série dans data (????)
+	varexp = varexp[-(nl + ij), ] #idem dans les variables exogènes (date 16/03/2015 + 6 jours
+	indNA = attr(na.omit(data[, 1:48]),"na.action")
+	data = data[-indNA,]
+	varexp = varexp[-indNA,]
+
+	# Conditionnement : les jours avec PMjour +/- large
+	large = 1
+	bornes = mean(dataj[25:48])+c(-large,large)
+	indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
+	data = data[indcond,]
+	varexp = varexp[indcond,]
+
+	D = rep(0,nrow(data))
+	for (k in 1:nrow(data))
+	{
+		#D[k] = sqrt(sum((1:H)*(dataj[L] - data[k,L])^2))
+		D[k] = sqrt(sum((dataj[L] - data[k,L])^2))
+	}
+	ind = order(D)[1:nbvois]
+	w = 1/(D[ind]^2)
+	w = w/sum(w)
+	W = w %o% rep(1,48)
+	JourMoy = apply(data[ind, 1:48], 2, mean)
+	#JourMoy = apply(W*data[ind, 1:48], 2, sum)
+	NomFile = paste("Voisins_Epandage_PMjour_Hc_",Hc,".png",sep="")
+	Titre = paste("Jour à prévoir : ",dateJPrev," - ", length(ind)," voisins",sep="")
+	#erreur = sqrt(sum((dataj[25:48] - JourMoy[25:48])^2))
+	if(Hc==24)
+		erreurPrev=NA
+	else
+		erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48]))
+	erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48]))
+	#png(NomFile)
+	matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind), 
+		cex.axis=1.4, cex.main = 1.7, cex.lab=1.5, 
+		xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1),
+			" / ",round(erreurPrev,1)), main=Titre)
+	legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind))
+	lines(1:48, dataj, lwd=2.5)
+	lines(JourMoy, lty = 2, lwd=2)
+	abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2)
+	#xx=dev.off()
+
+	Err24 = c(Err24, erreur24)
+	ErrPrev = c(ErrPrev, erreurPrev)
+	ResDates = cbind(ResDates, rownames(data)[ind])
 }
 
-
 rownames(ResDates) = 1:10
 
 Kvois = NULL
-for (Col in ncol(ResDates):1) {
-K = 0
-for (I in 1:10){
-    for (J in 1:10){
-       if (ResDates[I,1] == ResDates[J,Col]) K = K +1
-     }
-}
-Kvois = c(Kvois, K)
+for (Col in ncol(ResDates):1)
+{
+	K = 0
+	for (I in 1:10)
+	{
+		for (J in 1:10)
+		{
+			if (ResDates[I,1] == ResDates[J,Col])
+				K = K +1
+		}
+	}
+	Kvois = c(Kvois, K)
 }
 
 # pdf("Erreur_Epandage_PMjour.pdf")
@@ -133,4 +146,4 @@ legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lw
 points(5:24, Kvois, cex=1.8, pch = 19)
 # xx = dev.off()
 
-length(D)
+#length(D)
-- 
2.44.0