##########################################################
# Unsupervised classification of the altimetric wave forms 
# of the satellite Topex/Poseidon
# R/S+ commandlines
##########################################################
##
# Entering radar dataset
########################
SATDAT <- as.matrix(read.table("npfda-sat.dat"))
attributes(SATDAT)$dimnames <- NULL
##
# Automatical classification: the following commandline 
# allows to get directly the groups
#######################################################
sat.classif <- classif.npfda(SATDAT, 1:ncol(SATDAT), kind.of.kernel="quadratic", semimetric="hshift", threshold=0.05, nb.bw=50, nss=0, mspg=10, centrality="median")
##
# Conclusion: the classification procedure retains five groups:
#              Group11, Group12, Group21, Group22, Group23
###############################################################
Group11 <- sat.classif$Partition[[1]]
Group12 <- sat.classif$Partition[[2]]
Group21 <- sat.classif$Partition[[3]]
Group22 <- sat.classif$Partition[[4]]
Group23 <- sat.classif$Partition[[5]]
########################################################
# Plotting a sample of size 16 of each terminal leaves:# 
# Group11, Group12, Group21, Group22 and Group23       # 
########################################################
s11 <- sample(Group11,16)
s12 <- sample(Group12,16)
s21 <- sample(Group21,16)
s22 <- sample(Group22,16)
s23 <- sample(Group23,16)
x <- 1:ncol(SATDAT)
close.screen(all = TRUE)
screen1 <- c(0.10,0.42,0.55,0.88)
screen2 <- c(0.58,0.90,0.55,0.88)
screen3 <- c(0,0.32,0.11,0.44)
screen4 <- c(0.34,0.66,0.11,0.44)
screen5 <- c(0.68,1,0.11,0.44)
split.screen(rbind(screen1, screen2, screen3, screen4, screen5))
screen(1)
par(mar=c(0,0,0,0)+0.2, oma=c(0,0,2,0), xaxt="n", yaxt="n")
plot(x,x,type="n",bty="n")
mtext(c("GROUP 11","GROUP 12"), side=3, outer=T, line=-4, at=c(0.28,.72))
##
# GROUP11
#########
nbscreen <- split.screen(c(4,4),1)
screen(nbscreen[1])
plot(x, SATDAT[s11[1],],ylim=range(SATDAT), type="l", xlab="",ylab="")
count <- 0
for(i in s11[2:16]){
        count <- count + 1
        screen(nbscreen[count+1])
        plot(x,SATDAT[i,],ylim=range(SATDAT), type="l", xlab="",ylab="")
}
mtext(c("GROUP 21","GROUP 22", "GROUP 23"), side=3, outer=T, line=-22, at=c(0.16,0.5,0.84))
##
# GROUP12
#########
nbscreen <- split.screen(c(4,4),2)
screen(nbscreen[1])
plot(x, SATDAT[s12[1],],ylim=range(SATDAT), type="l", xlab="",ylab="")
count <- 0
for(i in s12[2:16]){
        count <- count + 1
        screen(nbscreen[count+1])
        plot(x,SATDAT[i,],ylim=range(SATDAT), type="l", xlab="",ylab="")
}
##
# GROUP21
#########
nbscreen <- split.screen(c(4,4),3)
screen(nbscreen[1])
plot(x, SATDAT[s21[1],],ylim=range(SATDAT), type="l", xlab="",ylab="")
count <- 0
for(i in s21[2:16]){
        count <- count + 1
        screen(nbscreen[count+1])
        plot(x,SATDAT[i,],ylim=range(SATDAT), type="l", xlab="",ylab="")
}
##
# GROUP22
#########
nbscreen <- split.screen(c(4,4),4)
screen(nbscreen[1])
plot(x, SATDAT[s22[1],],ylim=range(SATDAT), type="l", xlab="",ylab="")
count <- 0
for(i in s22[2:16]){
        count <- count + 1
        screen(nbscreen[count+1])
        plot(x,SATDAT[i,],ylim=range(SATDAT), type="l", xlab="",ylab="")
}
##
# GROUP23
#########
nbscreen <- split.screen(c(4,4),5)
screen(nbscreen[1])
plot(x, SATDAT[s23[1],],ylim=range(SATDAT), type="l", xlab="",ylab="")
count <- 0
for(i in s23[2:16]){
        count <- count + 1
        screen(nbscreen[count+1])
        plot(x,SATDAT[i,],ylim=range(SATDAT), type="l", xlab="",ylab="")
}
##
# Loading mean, median and modal curves for each group:
#######################################################
MEAN.CURVES <- sat.classif$MEANS
MEDIAN.CURVES <- sat.classif$MEDIANS
MODAL.CURVES <- sat.classif$MODES
###
# Displaying the corresponding modal curves
###########################################
par(mfrow=c(1,6), mar=c(1,2,2,0)+.2, oma=c(1,1,0,0), pty="s")
Labels <- sat.classif$Labels
for(j in 1:2)
    plot(x, MODAL.CURVES[j,],ylim=range(SATDAT), type="l", xlab="",ylab="", main=paste("GROUP",Labels[j],sep=""), cex.main=2)
for(j in 3:5)
    plot(x, MODAL.CURVES[j,],ylim=range(SATDAT), type="l", xlab="",ylab="", main=paste("GROUP",Labels[j],sep=""), cex.main=2 )
###
# Displaying the corresponding mean curves
########################################################
par(mfrow=c(1,6), mar=c(1,2,2,0)+.2, oma=c(1,1,0,0), pty="s")
Labels <- sat.classif$Labels
for(j in 1:2)
    plot(x, MEAN.CURVES[j,],ylim=range(SATDAT), type="l", xlab="",ylab="", main=paste("GROUP",Labels[j],sep=""), cex.main=2)
for(j in 3:5)
    plot(x, MEAN.CURVES[j,],ylim=range(SATDAT), type="l", xlab="",ylab="", main=paste("GROUP",Labels[j],sep=""), cex.main=2 )
###
# Computing mean and median curves for groups 1 and 2:
######################################################
nb.groups <- length(sat.classif$Partition)
MEAN.CURVES <- matrix(0, 2, ncol(SATDAT))
MEDIAN.CURVES <- matrix(0, 2, ncol(SATDAT))
MODAL.CURVES <- matrix(0, 2, ncol(SATDAT))
Group1 <- c(sat.classif$Partition[[1]], sat.classif$Partition[[2]])
Group2 <- c(sat.classif$Partition[[3]], sat.classif$Partition[[4]], sat.classif$Partition[[5]])
MEAN.CURVES[1,] <- apply(SATDAT[Group1,], 2, mean)
MEAN.CURVES[2,] <- apply(SATDAT[Group2,], 2, mean)
grid <- 1:ncol(SATDAT)
SM.SAT <- semimetric.hshift(SATDAT, SATDAT, grid)
MEDIAN.CURVES[1,] <- SATDAT[Group1[median.npfda(SM.SAT[Group1, Group1])],]
MEDIAN.CURVES[2,] <- SATDAT[Group2[median.npfda(SM.SAT[Group2, Group2])],]
###
# Computing modal curves for groups 1 and 2
###########################################
Hrange <- range(SM.SAT)
Bw.seq <- seq(Hrange[1], Hrange[2] * 0.5, length = 50 + 1)[-1]
n <- nrow(SM.SAT)
PROBCURVES <- matrix(0, n, 50)
for(i in 1:n){
       	PROBCURVES[i,] <- prob.curve(i, SM.SAT, Bw.seq)
}
res.classif.bw <- classif.bw(PROBCURVES, Bw.seq, 1:n)
Bw.opt <- res.classif.bw$bw
index <- res.classif.bw$index
shi <- classif.hi(SATDAT, SM.SAT, "quadratic", Bw.opt, 1:n, "hshift", "median", grid)
first.split.sat <- classif.part(SATDAT, PROBCURVES, SM.SAT, "quadratic", index, Bw.seq, 1:n, shi, "hshift", 0.05, 0, 10, "median", grid)
Group1 <- first.split.sat$Groups[[1]]
Group2 <- first.split.sat$Groups[[2]]
Band <- first.split.sat$Bw.opt
rank.mode1 <- funopa.mode(Band[1], SM.SAT[Group1,Group1], "quadratic")
MODAL.CURVES[1,] <- SATDAT[Group1[rank.mode1],]
rank.mode2 <- funopa.mode(Band[2], SM.SAT[Group2, Group2], "quadratic")
MODAL.CURVES[2,] <- SATDAT[Group2[rank.mode2],]
###
# Displaying modal, median and mean curves for groups 1 and 2 
#############################################################
par(mfrow=c(2,3), mar=c(3,2,2,0)+.2, oma=c(1,1,0,0), pty="m")
for(j in 1:2){
    plot(1:70, MODAL.CURVES[j,],ylim=range(SATDAT), type="l", xlab="",ylab="", main=paste("Modal curve of GROUP",j,sep=""), cex.main=1.5)
    plot(1:70, MEDIAN.CURVES[j,],ylim=range(SATDAT), type="l", xlab="",ylab="", main=paste("Median curve of GROUP",j,sep=""), cex.main=1.5)
    plot(1:70, MEAN.CURVES[j,],ylim=range(SATDAT), type="l", xlab="",ylab="", main=paste("Mean curve of GROUP",j,sep=""), cex.main=1.5)
}
###
# Displaying the splitting score behavior along the procedure
#############################################################
Splitting.score <- sat.classif$Ssc
Split.names <- c(paste("GROUPS","\n","1 & 2"), paste("GROUPS","\n","11 & 12"), paste("GROUPS", "\n","21, 22 & 23"))
par(mfrow=c(1,1))
barplot(Splitting.score, names=Split.names, main="Splitting scores", ylab="Percentages")
###
# Plotting a sample of wave forms
#################################
ylimits <- range(SATDAT)
sample.sat <- c(18, 10, 12, 44, 133, 9, 24, 96, 7, 22, 15, 18, 33, 132, 21, 1, 2, 3, 4, 5)
par(mfrow=c(4,5), mar=c(1,1,4,2), oma=c(1,2,0,0), mex=.6)
for(i in 1:20) plot(1:70, SATDAT[sample.sat[i],], xlab="",ylab="", type="l", ylim = ylimits, main=paste(sample.sat[i]), cex.main=1.5)