## Copyright 2017 Irene Epifanio

## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.


###################
#Loading functions 
###################

source("LstepArchetypesMod.R")
source("LstepArchetypoids.R")
source("LFunctions_to_calculate_real_archetypes_with_swap2.R")


  CalculaExponencialDeFormas = function( logba, base ) {
  #### cat( "\n" )
  #### cat( "dim of logba: ", dim( logba ), "\n" )
  #### cat( "dim of base:  ", dim( base ), "\n" )
  #### logba está en forma de vector
  #### base en forma de matriz
  m = dim( base )[ 2 ]
  k = dim( base )[ 1 ]+1
  #### print(m)
  #### print(k)
  gamma = as.vector( base )
  dentro = ( sqrt( 1 - t( logba ) %*% logba ) ) * gamma + logba
  #### print(dim(dentro))
  Z = matrix( dentro, nrow = k-1 )
  return( Z )
}

my_defh = function( nrow ) {

    #### # Set timer.
    #### t1 = proc.time()[ "elapsed" ]

    k <- nrow
    h <- matrix( 0, k, k + 1 )
    if( nrow > 0 ) {
      for( j in seq( 1, k ) ) {
          val = -1 / sqrt( j * ( j + 1 ) )
          h[ j, seq( 1, j ) ] = val
          h[ j, j+1 ]         = - j * val
      }
    }

    #### # Print elapsed time.
    #### t2 = proc.time()[ "elapsed" ]
    #### cat( "Elapsed time in my_defh (seconds): ", t2-t1, "\n" )

    h
}


library(Anthropometry)
library(shapes)
library(archetypes)

###################
#Data set
###################
d3=digit3.dat
d3o=d3[,,-1] #removing the first digit, outlier

###################
#Computing tangent coordinates and mean
###################
formas=procGPA(d3o,tangentcoords = "partial", pcaoutput=TRUE) 
  
tangentes=formas$tan

media=formas$mshape

X=t(tangentes)

###################
#PCA 
###################
shapepca(formas,type="r",mag=1)

#figures PCA
source('dibujosaxes.R')
 shapepcar(formas,pcno=1:4,type="r",mag=1,joinline=1:13)

#distance to the mean
formas$rho

cumsum(formas$percent)
formas$percent
rank(formas$rho)

#######################
#Sparse PCA
#######################
library(elasticnet)
out2<-spca(X,K=4,type="predictor",sparse="varnum",trace=TRUE,para=c(6,6,6,6)) #6 Num. of non-zero loadings


######Plot sparse pca

formas2=formas

formas2$pcar=out2$loadings
formas2$pcasd=sqrt(out2$var.all*out2$pev)
formas2$percent=out2$pev*100

shapepca(formas2,type="r",mag=.3)
shapepca(formas2,type="v",mag=.3)

source('dibujosaxes.R')
 shapepcar(formas2,pcno=1:4,type="r",mag=.2,joinline=1:13)


###############
#k-means landmarks -> Vinue,Simo,Alemany 2014
###################

x=d3o
Nsteps=10
Nclusters=4
n=dim(x)[3]
k=dim(x)[1]
d=dim(x)[2]

set.seed(478392)

meanshapes=0
asig=0
iniciales=sample(c(1:n),Nclusters,replace=F)
meanshapes=x[,,iniciales]
dist=matrix(0,n,Nclusters)
asig=0

for (step in 1:Nsteps)
{

    for(j in 1:Nclusters){
       dist[,j]=apply(x[,,1:n],3,riemdist,y=meanshapes[,,j])
    }
    asig=max.col(-dist)
    numero=as.vector(table(asig));

    for(j in 1:Nclusters){
        if(numero[j]>1){
            meanshapes[,,j]=procGPA(x[,,asig==j], distances=F, pcaoutput=F)$mshape
        }
       if(numero[j]==1){
            meanshapes[,,j]=x[,,asig==j]
        }
    }
}


#################
#figures k-means

#The full Procrustes coordinates 
cm=centroid.size(meanshapes)

par(mfrow=c(2,2))
 plotshapes(meanshapes[,,1]/cm[1],joinline=c(1:13))
 plotshapes(meanshapes[,,2]/cm[2],joinline=c(1:13))
 plotshapes(meanshapes[,,3]/cm[3],joinline=c(1:13))
 plotshapes(meanshapes[,,4]/cm[4],joinline=c(1:13))

par(mfrow=c(2,2))
 plotshapesa3(meanshapes[,,1]/cm[1],joinline=c(1:13))
 plotshapesa3(meanshapes[,,2]/cm[2],joinline=c(1:13))
 plotshapesa3(meanshapes[,,3]/cm[3],joinline=c(1:13))
 plotshapesa3(meanshapes[,,4]/cm[4],joinline=c(1:13))



###############################
#Archetypes with landmarks
################################
norep=20
#Archetypes screeplot for determining K
set.seed(2016)
lass10 <- LstepArchetypesMod(data=X,k=1:10,verbose=FALSE,nrep=norep,saveHistory=F) 
 screeplot(lass10)


#Obtaining archetypes for a certain number (lass10 executed for 1:10), in array arquetipos

numarq=4; #elbow
 
a=bestModel(lass10[[numarq]])


#explained
1 -bestModel(lass10[[numarq]])$rss/(norm2.Frnormfn(X)/dim(X)[1])


#Archetypes in arquetipos array
q0=parameters(a) #(a$betas)%*%t(tangentes)
k=13
m=2

 arquetipos=array(rep(0,k*m*numarq),dim=c(k,m,numarq))
 for( i in 1:numarq)
 {
   
   q1 = CalculaExponencialDeFormas( q0[i,],preshape(media) )
   h <- my_defh( k - 1 )
   q1 <- t( h ) %*% q1
   arquetipos[,,i]=q1
}


###################
#Archetypoids with landmarks
###################
#finding RSS-archetypoids for different K, beginning for three initializations and using the minimum RSS of each three. Screeplot.

set.seed(2016)
norep=20
lass10d <- stepLArchetypoids3(data=X, k=1:10, norep)

screeplot(lass10d)

#elbow
numarq=4; #elbow

#narquetipoides contains the individual numbers which are archetypoids
narquetipoides=lass10d[[numarq]][[1]]
# 2 28 15  7



#Plots archetypes and archetypoids
plotshapes(arquetipos,joinline=c(1:13))

par(mfrow=c(2,2))
 plotshapes(arquetipos[,,1],joinline=c(1:13))
 plotshapes(arquetipos[,,2],joinline=c(1:13))
 plotshapes(arquetipos[,,3],joinline=c(1:13))
 plotshapes(arquetipos[,,4],joinline=c(1:13))


par(mfrow=c(2,2))
 plotshapesa3(arquetipos[,,1],joinline=c(1:13))
 plotshapesa3(arquetipos[,,2],joinline=c(1:13))
 plotshapesa3(arquetipos[,,3],joinline=c(1:13))
 plotshapesa3(arquetipos[,,4],joinline=c(1:13))



plotshapes(d3o[,,narquetipoides],joinline=c(1:13))

windows()
par(mfrow=c(2,2))

plotshapes(d3o[,,narquetipoides[1]],joinline=c(1:13))
plotshapes(d3o[,,narquetipoides[2]],joinline=c(1:13))
plotshapes(d3o[,,narquetipoides[3]],joinline=c(1:13))
plotshapes(d3o[,,narquetipoides[4]],joinline=c(1:13))




par(mfrow=c(2,2))

plotshapesa3(formas$rot[,,narquetipoides[1]],joinline=c(1:13))
plotshapesa3(formas$rot[,,narquetipoides[2]],joinline=c(1:13))
plotshapesa3(formas$rot[,,narquetipoides[3]],joinline=c(1:13))
plotshapesa3(formas$rot[,,narquetipoides[4]],joinline=c(1:13))


#################
#alphas archetypes, barplot with purity order

al=a$alphas
 
 gini <- function(x, unbiased = TRUE, na.rm = FALSE){
    if (!is.numeric(x)){
        warning("'x' is not numeric; returning NA")
        return(NA)
    }
    if (!na.rm && any(na.ind <- is.na(x)))
        stop("'x' contain NAs")
    if (na.rm)
        x <- x[!na.ind]
    n <- length(x)
    mu <- mean(x)
    N <- if (unbiased) n * (n - 1) else n * n
    ox <- x[order(x)]
    dsum <- drop(crossprod(2 * 1:n - n - 1,  ox))
    dsum / (mu * N)
}

ag=apply(al,1,FUN=gini)

aog=al[order(ag,decreasing=T),]
barplot(t(aog),col=rainbow(4))

order(ag,decreasing=T)

###################
#alphas for the mean

numarq=4; #elbow
 
a=bestModel(lass10[[numarq]])
q0=parameters(a) #(a$betas)%*%t(tangentes)


K=numarq
huge=200
# archs_ini <- lass10d[[numarq]][[1]]
 data=matrix(0,ncol=24,nrow=1) #mean
  
  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
zs=  rbind(t(q0), rep(huge, n))

#  zs <- x_gvv[, archs_ini] 
  zs <- as.matrix(zs)
  
  alphas <- matrix(0, nrow = K, ncol = n)
  for (j in 1 : n){
    alphas[, j] = coef(nnls(zs, x_gvv[,j]))
  }

 alphas
###############


##########################

#Paralelel coordinates:
matplot(t(a$alphas),type='l',ylab='alphas')

#Simplexplot
simplexplot(a)

