In martinEcarnot/signe: Useful functions to handle spectra

Load libraries

library(MASS) # For ginv function
# library(FactoMineR)
library(signal)
# library(plyr)
library(caret)
# library(dplyr)
library(prospectr)
library(rnirs)
# library(ggplot2)
# library(plotly)
library(nirsextra)

source("Script_R_2020/sp2dfclo.R")
source('Script_R_2020/SIGNE_maha0.R')

Load globalmatix

# brb3="/home/ecarnot/Documents/INRA/Projets/SIGNE/2019/globalmatrix"
brb3="globalmatrix"
load(file=brb3)
# sp=globalmatrix

Set parameters

### FIXATION DES PARAMETRES UTILISES:
## Nombre de repetitions de la boucle de PLSDA:
repet= 2
## Parametres du Savitsky-Golay (p=degre du polynome, n= taille de la fenetre, m=ordre de derivation)
p=3 #2
n=11#11 #Faire avec un n plus gros ?
m=1 #1
## Nombre de VL max autorisees
ncmax=30 #35

Pretreatments

# p=rbind(list('adj',''),list('red',c(1000,1,1)),list('snv',''),list('sder',c(2,3,w)),list('red',c(w/2,w/2,1)))
pr=rbind(list('adj',c(551,1352)),list('adj',c(552,1352)),list('snv',''),list('sder',c(m,p,n)),list('red',c(n/2,n/2,1)))
# pr=rbind(list('adj',c(551,1352)),list('adj',c(552,1352)),list('snv',''),list('sder',c(m,p,n)))
sp_pre=pre(globalmatrix,pr)

Constitution d' jeu de calibration

sp=sp_pre
cep=as.factor(substr(rownames(sp),9,9))
clo=as.factor(substr(rownames(sp),9,13))
sp=sp2dfclo(sp,cep,clo)
sp$annee=as.factor(substr(rownames(sp),1,4))
sp$parcelle=as.factor(substr(rownames(sp),18,18))
sp$date=as.factor(substr(rownames(sp),1,8))
## On enleve 20190628, jour a + de 43°C ?
sp=sp[sp$date != "20190628" & sp$parcelle=="G",]
# sp=sp[sp$parcelle=="G",]
# sp=sp[sp$parcelle=="G" & sp$date != "20190730",]

annees=c("2017","2018","2019")
for (i in 1:3) {
  idcal= which(sp$annee !=annees[i])
spcal= droplevels(sp[idcal,])
spval=droplevels(sp[-idcal,])
predmF=as.data.frame(matrix(nrow = nrow(spval), ncol = ncmax))
distances=as.data.frame(matrix(nrow = nrow(spval), ncol = ncmax))

rplsda=caret::plsda(spcal$x, spcal$cep,ncomp=ncmax)
sccal=rplsda$scores
spval_c=scale(spval$x,center=rplsda$Xmeans,scale = F)
scval=spval_c%*%rplsda$projection  # score_val=predict(rplsda,sc_val,type="scores") : ne marche pas

for (ii in 2:ncmax) {
  predmF[,ii]=SIGNE_maha0(sccal[,1:ii], spcal$cep, scval[,1:ii])$class
  distances[,ii]=SIGNE_maha0(sccal[,1:ii], spcal$cep, scval[,1:ii])$dist
}

tsm=lapply(as.list(predmF), spval$cep, FUN = table)
diagsm=lapply(tsm, FUN = diag)
perokm =100*unlist(lapply(diagsm, FUN = sum))/nrow(spval)
perokmT=perokm
print(perokm)
}

Pour justifier PLSDA externe seuelemnt CV pour montrer

# On moyenne par feuille pour ne pas avoir un spectre d'une mm feuille en cal et en val
pos_feuille=floor((as.numeric(substr(rownames(sp),15,16))+2)/3)
date_cep_feuill=paste(substr(rownames(sp),1,13),pos_feuille, sep = '-')
spm=aggregate(sp[,3],list(date_cep_feuill),mean)
# spm=aggregate(sp[,3],list(date_cep),head,1)

rownames(spm)=spm[,1]
cep=as.factor(substr(rownames(spm),9,9))
clo=as.factor(substr(rownames(spm),9,13))
spm=sp2dfclo(as.matrix(spm[,-1]),cep,clo)
spm$annee=as.factor(substr(rownames(spm),1,4))
spm$date=as.factor(substr(rownames(spm),1,8))

idcal= which(spm$annee !=annees[3])
spcal= droplevels(spm[idcal,])
spval=droplevels(spm[-idcal,])


nblock=5
seg=segmkf(nrow(spcal),K=nblock,type = "random", nrep=20)
r=1
b=1
perokmT=NULL
for (j in 1:length(seg)) {
spcal1=spcal[-unlist(seg[[r]][b]),]
spval1=spcal[unlist(seg[[r]][b]),]
rplsda=caret::plsda(spcal1$x, spcal1$cep,ncomp=ncmax)
sccal=rplsda$scores
spval_c=scale(spval1$x,center=rplsda$Xmeans,scale = F)
scval=spval_c%*%rplsda$projection  # score_val=predict(rplsda,sc_val,type="scores") : ne marche pas

predmF=as.data.frame(matrix(nrow = nrow(spval1), ncol = ncmax))
for (ii in 2:ncmax) {
  predmF[,ii]=SIGNE_maha0(sccal[,1:ii], spcal1$cep, scval[,1:ii])$class
  # distances[,ii]=SIGNE_maha0(sccal[,1:ii], spcal1$cep, scval[,1:ii])$dist
}

tsm=lapply(as.list(predmF), spval1$cep, FUN = table)
diagsm=lapply(tsm, FUN = diag)
perokm =100*unlist(lapply(diagsm, FUN = sum))/nrow(spval1)
perokmT=rbind(perokmT,perokm)
b=b+1
if (b==nblock) {b=1;r=r+1}
}
# print(colMeans(perokmT))

# Validation
rplsda=caret::plsda(spcal$x, spcal$cep,ncomp=ncmax)
sccal=rplsda$scores
spval_c=scale(spval$x,center=rplsda$Xmeans,scale = F)
scval=spval_c%*%rplsda$projection  # score_val=predict(rplsda,sc_val,type="scores") : ne marche pas
predmF=as.data.frame(matrix(nrow = nrow(spval), ncol = ncmax))
for (ii in 2:ncmax) {
  predmF[,ii]=SIGNE_maha0(sccal[,1:ii], spcal$cep, scval[,1:ii])$class
}
tsm=lapply(as.list(predmF), spval$cep, FUN = table)
diagsm=lapply(tsm, FUN = diag)
perokm =100*unlist(lapply(diagsm, FUN = sum))/nrow(spval)
# print(perokm)

plot(colMeans(perokmT))
points(perokm)