##########################################################
# Unsupervised classification of the spectrometric curves
# R/S+ commandlines
##########################################################
##
# Entering spectrometric dataset
##################"#############
SPECDAT <- as.matrix(read.table("npfda-spectrometric.dat"))
attributes(SPECDAT)$dimnames <- NULL
SPECURVES <- SPECDAT[,1:100]
##
# Automatical classification: the following commandline 
# allows to get directly the groups
#######################################################
spec.classif <- classif.npfda(SPECURVES, 2, 20, c(0,1), semimetric="deriv", threshold=0.05, nb.bw=100, nss=0, mspg=10, centrality="mean")
##
# Conclusion: the classification procedure retains three groups:
#              Group11, Group21, Group22
################################################################
Group1 <- spec.classif$Partition[[1]]
Group21 <- spec.classif$Partition[[2]]
Group22 <- spec.classif$Partition[[3]]
#########################################################
# Plotting terminal leaves: Group1, Group21 and Group22 # 
#########################################################
###
# Computing the second derivatives of spectrometric curves
##########################################################
SPECD2 <- t(approx.spline.deriv(SPECURVES, 2, 20, c(0,1))$APPROX)
###
# Loading the modal curve for each group
########################################
MODAL.CURVES <- spec.classif$MODES
###
# Computing the second derivatives of modal curves
##################################################
MODAL.CURVESD2 <- t(approx.spline.deriv(MODAL.CURVES, 2, 20, c(0,1))$APPROX)
###
# Displaying groups: original data, second derivatives 
# and corresponding modal curves
######################################################
x <- seq(850, 1050, length=100)
par(mfrow=c(3,3))
##
# GROUP1
########
plot(x,SPECURVES[Group1[1],],ylim=range(SPECURVES), type="l", xlab="Wavelenghts",ylab="Absorbances", main="GROUP 1")
for(i in Group1[-1]){
        par(new=T)
        plot(x,SPECURVES[i,],ylim=range(SPECURVES), type="l", xlab="",ylab="", axes=F)
}
plot(x,SPECD2[Group1[1],],ylim=range(SPECD2), type="l", xlab="Wavelenghts",ylab="Absorbances", main="GROUP 1: second derivatives")
for(i in Group1[-1]){
        par(new=T)
        plot(x,SPECD2[i,],ylim=range(SPECD2), type="l", xlab="",ylab="", axes=F)
}
plot(x, MODAL.CURVESD2[1,], ylim=range(SPECD2), type="l", xlab="Wavelenghts", ylab="Absorbances", main="GROUP 1: modal curve")
##
# GROUP21
#########
plot(x,SPECURVES[Group21[1],],ylim=range(SPECURVES), type="l", xlab="Wavelenghts",ylab="Absorbances", main="GROUP 21")
for(i in Group21[-1]){
        par(new=T)
        plot(x,SPECURVES[i,],ylim=range(SPECURVES), type="l", xlab="",ylab="", axes=F)
}
plot(x,SPECD2[Group21[1],],ylim=range(SPECD2), type="l", xlab="Wavelenghts",ylab="Absorbances", main="GROUP 21: second derivatives")
for(i in Group21[-1]){
        par(new=T)
        plot(x,SPECD2[i,],ylim=range(SPECD2), type="l", xlab="",ylab="", axes=F)
}
plot(x, MODAL.CURVESD2[2,], ylim=range(SPECD2), type="l", xlab="Wavelenghts", ylab="Absorbances", main="GROUP 21: modal curve")
##
# GROUP22
#########
plot(x,SPECURVES[Group22[1],],ylim=range(SPECURVES), type="l", xlab="Wavelenghts",ylab="Absorbances", main="GROUP 22")
for(i in Group22[-1]){
        par(new=T)
        plot(x,SPECURVES[i,],ylim=range(SPECURVES), type="l", xlab="",ylab="", axes=F)
}
plot(x,SPECD2[Group22[1],],ylim=range(SPECD2), type="l", xlab="Wavelenghts",ylab="Absorbances", main="GROUP 22: second derivatives")
for(i in Group22[-1]){
        par(new=T)
        plot(x,SPECD2[i,],ylim=range(SPECD2), type="l", xlab="",ylab="", axes=F)
}
plot(x, MODAL.CURVESD2[3,], ylim=range(SPECD2), type="l", xlab="Wavelenghts", ylab="Absorbances", main="GROUP 22: modal curve")
###
# Displaying the splitting score behavior along the procedure
#############################################################
Splitting.score <- spec.classif$Ssc
Split.names <- c(paste("GROUPS","\n","1 & 2"), paste("GROUPS", "\n","21 & 22"))
par(mfrow=c(1,1))
barplot(Splitting.score, names=Split.names, main="Splitting scores", ylab="Percentages")
abline(h=0.05, lwd=3)