+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"
+
+# Chargement des données météo et indicateurs
+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)
+dates = dates[,1]
+
+# Contruction des matrices de données
+pm.h <- matrix(pm[,2],ncol=24,byrow=TRUE)
+
+Nlignes = nrow(pm.h)
+Data = cbind(pm.h[1:(Nlignes -1), ], pm.h[2:Nlignes, ])
+dates2 = dates[2:Nlignes]
+rownames(Data) = dates2
+# df contient l'ensemble des données.
+#df <- cbind(Data,varexp[,-1])
+df <- Data
+PMjour <- apply(df[,25:48],1,mean,na.rm=T)
+dfexp <- cbind(VarExp,PMjour)
+
+Dates = c(
+"16/03/2015",
+"19/01/2015",
+"27/04/2015")
+
+Categorie = c("Epandage", "Chauffage", "Non Polluée")
+
+RepFig = "FIGURES_Etude"
+
+ResDates = NULL
+
+nbvois=10
+j=1 # numéro de semaine
+ij=6 # numéro du jour (0 = lundi)
+
+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])
+}
+
+
+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)
+}
+
+# pdf("Erreur_Epandage_PMjour.pdf")
+ymin = min(na.omit(ErrPrev), Err24)
+ymin = min(ymin, Kvois)
+ymax = max(na.omit(ErrPrev), Err24)
+ymax = max(ymax, Kvois)
+plot(5:24,Err24, type = "l", lwd = 2, cex.axis=1.4, cex.lab = 1.5,
+ ylab = "MAE",xlab = "Heures de prévision", ylim= c(ymin,ymax))
+lines(5:24,ErrPrev, lwd=2, lty = 2)
+legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lwd=2, cex=1.5)
+points(5:24, Kvois, cex=1.8, pch = 19)
+# xx = dev.off()
+
+length(D)