
#As model 1 from Arribas and Romo (2014)

source('LFunctions_to_calculate_real_archetypes_with_swap2.R')
source('LstepArchetypesMod.R')
source('LstepArchetypoids.R')
source('fhdr.R')
source('robout.R')


library(archetypes)
library(Anthropometry)
library( roahd )
library(rainbow)
library(mrfDepth)
library(fda)
library(ks)
library(cluster)
library(fda.usc)


# Auxiliary rep.col and rep.row functions (http://www.r-bloggers.com/a-quick-way-to-do-row-repeat-and-col-repeat-rep-row-rep-col/)
rep.row<-function(x,n){
  matrix(rep(x,each=n),nrow=n)
}
rep.col<-function(x,n){
  matrix(rep(x,each=n), ncol=n, byrow=TRUE)
}

# Covariance functions      
cov.fun=function(d,k,c,mu){ #Modelo Hydmann y Shang 2010
        k*exp(-c*d^mu)
}
cov.fun2=function(d,k,c){   #Modelo Fraiman y Svarc 2013
        k*exp(-c*d)
}

# Contamination models
f.out1=function(t,m){            
	u<-runif(1,0,1) 
	if (u<0.5){
		f<-4*t+exp(-(t-m)^2/0.02)/sqrt(2*pi*0.01)-1.8
	}
	else{f<-4*t-exp(-(t-m)^2/0.02)/sqrt(2*pi*0.01)+1.8}
	return(f)
	}
f.out2=function(t,th){4*t+2*sin(4*(t+th)*pi)}
f.out5=function(t){30*(1-t)*t^(3/2)}  





##############################################
#Figures example: Sect. 2.4.3 Functional outlier detection by archetype analysis
#############################################



#Figure 1

ite=1
K=2
n=100
out.rate=.02


p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)


set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    


p=length(t)
n=ncol(y)


 X=t(y)
 set.seed(2018)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=2, norep)
#lass10d=LstepArchetypesMod(data=X, k=2, norep)



a2a=t(y[,lass10d[[1]][[1]]])
#a2=bestModel(lass10d)
#a2a=a2$archetypes



matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, a2a[1,], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #archet
lines(t, a2a[2,], type = "l",lty=2,xlab="",ylab="",col=3,lwd=3) #archet


alfa=seq(0,1,0.01)
windows()

matplot(t, y, type = "n",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
for(i in 1:101){
lines(t, alfa[i]*a2a[1,]+(1-alfa[i])*a2a[2,], type = "l",lty=1,xlab="",ylab="",col=5,lwd=3) 
}
lines(t, a2a[1,], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #archetypoids
lines(t, a2a[2,], type = "l",lty=2,xlab="",ylab="",col=3,lwd=3) #archetypoids



#Figure 2

ite=1
K=3
n=100
out.rate=.02


p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)


set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers
trueout<-sort(sample(1:n,n*out.rate))
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

for (i in trueout){
	e=rnorm(p)
    y[,i]= f.out5(t)+t(L)%*%e                     #Outliers for Model 1                    
}

p=length(t)
n=ncol(y)


 X=t(y)
 set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')


 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, y[,trueout[1]], type = "l",lty=1,xlab="",ylab="",col=1) #outliers
lines(t, y[,trueout[2]], type = "l",lty=1,xlab="",ylab="",col=1) #outliers


matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, a3[1,], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #archetypoids
lines(t, a3[2,], type = "l",lty=2,xlab="",ylab="",col=3,lwd=3) #archetypoids
lines(t, a3[3,], type = "l",lty=3,xlab="",ylab="",col=4,lwd=3) #archetypoids

  

#############################
#with out.rate=0-> to see ternaryplot without holes

#Figure 3


K=3
n=100
out.rate=0
ite=1

set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers
trueout<-sort(sample(1:n,n*out.rate))
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

for (i in trueout){
	e=rnorm(p)
    y[,i]= f.out5(t)+t(L)%*%e                     #Outliers for Model 1
                         
}

p=length(t)
n=ncol(y)


 X=t(y)
 set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')



matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, a3[1,], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #archetypoids
lines(t, a3[2,], type = "l",lty=2,xlab="",ylab="",col=3,lwd=3) #archetypoids
lines(t, a3[3,], type = "l",lty=3,xlab="",ylab="",col=4,lwd=3) #archetypoids


####Figure 4


th=0.5
ite=1
K=3
n=100
out.rate=.02

no=round(n*out.rate)


p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)


set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers as Arribas (same variance, different mean)
trueout1<-sort(sample(1:n,n*out.rate))
#trueout1=1:no
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

for (i in trueout1){
	e=rnorm(p)
    y[,i]= f.out5(t)+t(L)%*%e                     #Outliers for Model 1                    
}


### Generate other kind of Outliers

#same mean, different variance
trueout2=sort(sample(1:n,n*out.rate))
trueout=c(trueout1,trueout2) # 7 67 72 92


t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

e=matrix(rnorm(no*p),p,no)

#y[,trueout2]=mu+t(L)%*%e +e
y[,trueout2]=mu +e


 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, y[,trueout1[1]], type = "l",lty=1,xlab="",ylab="",col=1) #outliers 1
lines(t, y[,trueout1[2]], type = "l",lty=1,xlab="",ylab="",col=1) #outliers 1


lines(t, y[,trueout2[1]], type = "l",lty=1,xlab="",ylab="",col=2) #outliers 2
lines(t, y[,trueout2[2]], type = "l",lty=1,xlab="",ylab="",col=2) #outliers 2


p=nrow(y)
n=ncol(y)


 X=t(y)
 set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#out: 7, 67




#Remove out, and repeat the procedure

y=y[,-out]

p=length(t)
n=ncol(y)


 X=t(y)
set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#no more outliers in shape, but due to variance?

yd=diff(y)
y=yd


p=nrow(y)
n=ncol(y)

t=seq(0,1,len=p)

 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data

lines(t, y[,trueout2[1]-2], type = "l",lty=1,xlab="",ylab="",col=2) #outliers 2
lines(t, y[,trueout2[2]-2], type = "l",lty=1,xlab="",ylab="",col=2) #outliers 2


 X=t(y)
set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#out 90 (originally 92)

#remove outlier
y=y[,-out]
p=nrow(y)
n=ncol(y)


 X=t(y)
set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#out 70 (originally 72)


#remove outlier
y=y[,-out]
p=nrow(y)
n=ncol(y)


 X=t(y)
set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#out: 7 points->no more outliers



#Figure 5
################
#amplitude outlier
################



th=0.6
ite=1
K=3
n=100
out.rate=.02

no=round(n*out.rate)


p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)


set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers as Arribas (same variance, different mean)
trueout1<-sort(sample(1:n,n*out.rate))
#trueout1=1:no
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

for (i in trueout1){
	e=rnorm(p)
    y[,i]= f.out5(t)+t(L)%*%e                     #Outliers for Model 1                    
}


### Generate other kind of Outliers

#amplitude
trueout2=sort(sample(1:n,n*out.rate))
trueout=c(trueout1,trueout2) # 7 67 72 92


for (i in trueout2){
	
    y[,i]= y[,i]+2                                    
}



 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, y[,trueout1[1]], type = "l",lty=1,xlab="",ylab="",col=1) #outliers 1
lines(t, y[,trueout1[2]], type = "l",lty=1,xlab="",ylab="",col=1) #outliers 1


lines(t, y[,trueout2[1]], type = "l",lty=1,xlab="",ylab="",col=2) #outliers 2
lines(t, y[,trueout2[2]], type = "l",lty=1,xlab="",ylab="",col=2) #outliers 2





p=length(t)
n=ncol(y)


 X=t(y)
 set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#out 72,92 , amplitude outliers


#Remove out, and repeat the procedure

y=y[,-out]

p=length(t)
n=ncol(y)


 X=t(y)
set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#7, 67, shape outliers

#Remove out, and repeat the procedure

y=y[,-out]

p=length(t)
n=ncol(y)


 X=t(y)
set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#no more outliers 



################
#isolated outlier: different in some interval
################

th=0.6
ite=1
K=3
n=100
out.rate=.02

no=round(n*out.rate)


p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)


set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers as Arribas (same variance, different mean)
trueout1<-sort(sample(1:n,n*out.rate))
#trueout1=1:no
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

for (i in trueout1){
	e=rnorm(p)
    y[,i]= f.out5(t)+t(L)%*%e                     #Outliers for Model 1                    
}


### Generate other kind of Outliers
out.rate=.05
trueout2=sort(sample(1:n,n*out.rate))
trueout=c(trueout1,trueout2) #7 67 22 72 77 90 92


for (i in trueout2){
	e=rnorm(p)
    y[,i]=  y[,i]+10*c(dnorm(seq(-2,2,.3)),rep(0,36))
                               
}



 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, y[,trueout1[1]], type = "l",lty=1,xlab="",ylab="",col=1) #outliers 1
lines(t, y[,trueout1[2]], type = "l",lty=1,xlab="",ylab="",col=1) #outliers 1

for(i in trueout2){
lines(t, y[,i], type = "l",lty=1,xlab="",ylab="",col=2) #outliers 2
}



p=length(t)
n=ncol(y)


 X=t(y)
 set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#out 7,67 , shape outliers


#Remove out, and repeat the procedure

y=y[,-out]

p=length(t)
n=ncol(y)


 X=t(y)
set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])


# 21 70 75 88 90: isolated outliers


#Remove out, and repeat the procedure

y=y[,-out]

p=length(t)
n=ncol(y)


 X=t(y)
set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])


#no more outliers

################
#shift outlier: same shape but translated in time
################


th=0.6
ite=1
K=3
n=100
out.rate=.02

no=round(n*out.rate)


p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)


set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers as Arribas (same variance, different mean)
trueout1<-sort(sample(1:n,n*out.rate))
#trueout1=1:no
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

for (i in trueout1){
	e=rnorm(p)
    y[,i]= f.out5(t)+t(L)%*%e                     #Outliers for Model 1                    
}


### Generate other kind of Outliers
out.rate= 0.03
trueout2=sort(sample(1:n,n*out.rate)) #  
trueout=c(trueout1,trueout2) #  7 67 72 90 92



ts=seq(-0.1,0.9,len=p)

for (i in trueout2){
	e=rnorm(p)
    y[,i]=  30*ts*(1-ts)^(3/2) + t(L)%*%e   
                               
}



 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data

for(i in trueout1){
lines(t, y[,i], type = "l",lty=1,xlab="",ylab="",col=1) #outliers 1
}


for(i in trueout2){
lines(t, y[,i], type = "l",lty=1,xlab="",ylab="",col=2) #outliers 2
}

lines(t, y[,72], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #outliers, not detected

p=length(t)
n=ncol(y)


 X=t(y)
 set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#out  7 67 90 92



#Remove out, and repeat the procedure

y=y[,-out]

p=length(t)
n=ncol(y)


 X=t(y)
set.seed(1234)
 norep=20
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)



a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

library(vcd)
ternaryplot(t(alphascce),grid=T,dimnames=c(1,2,3),col=1,cex=.6,main='')




wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])






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

########################################################
#Table 1 simulated data in Section 3.2 as Arribas et al.
########################################################


#Table 1

out.rate=.05 #0, 0.15


VP<-c()
FP<-c()
VPr<-c()
FPr<-c()

p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)



th=0.5
K=3
n=100
nrun=100
 norep=20


VP.or<-c()
FP.or<-c()
VPr.or<-c()
FPr.or<-c()
  
for (ite in 1:nrun) {
		cat(ite, '\n')
		set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers
trueout<-sort(sample(1:n,n*out.rate))
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

for (i in trueout){
	e=rnorm(p)
    y[,i]= f.out5(t)+t(L)%*%e                     #Outliers for Model 1               
}

p=length(t)
n=ncol(y)

#################################### ADA

 X=t(y)
 set.seed(1234)
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)

a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#outliers by threshold
aa.faltan<-length(setdiff(trueout,out))
aa.sobran<-length(setdiff(out,trueout))


#outliers by multivariate: Mah
sco=t(alphascce)
rownames(sco) = as.numeric(1:n)
s = cbind(sco, rep(1, n))

if( (out.rate == 0.15) & (ite == 35)){outr = robout(s, 1, "mve") # error with mcd Error in solve.default(cov, ...) 
}else{
outr = robout(s, 1, "mcd")}



aa.faltanr<-length(setdiff(trueout,outr$outliers))
aa.sobranr<-length(setdiff(outr$outliers,trueout))



################ OutlierGram (Arribas and Romo, 2014)

fD = fData( t, t(y) )
oo=outliergram( fD, display = F )
sout.OG.faltan<-length(setdiff(trueout,oo$ID_outliers))
sout.OG.sobran<-length(setdiff(oo$ID_outliers,trueout))


################ Robust Mahalanobis Distance (See Hydmann and Shang 2010)  

yf=fds(t,y)
 colnames(yf$y)=1:n

FO.robMah=foutliers(yf,method = "robMah") 
sout.robMah.faltan<-length(setdiff(trueout,FO.robMah$outliers))
sout.robMah.sobran<-length(setdiff(FO.robMah$outliers,trueout))

################ Integrated Square Error method (Hyndman and Ullah (2007); See Hydmann and Shang 2010)        
FO.HU<-foutliers(data=yf, method = "HUoutliers",order = 2, lambda = 3.29)
soutHU.faltan<-length(setdiff(trueout,FO.HU$outliers))
soutHU.sobran<-length(setdiff(FO.HU$outliers,trueout))


#################################################################
x=t(y)
rownames(x)<-as.character(1:n)

################ Package fda.usc          
Datosf<-fdata(x,t)

############################LRT (see Febrero et al. 2007)
flo=outliers.lrt(Datosf,nb=200,smo=0.05,dfunc=depth.mode,trim=out.rate)
#flo=foutliers(yf,method = "lrt",trim=out.rate) 
soutFLO.faltan<-length(setdiff(trueout,flo$outliers))
soutFLO.sobran<-length(setdiff(flo$outliers,trueout))


################ Depth Based Trimming (Febrero, Galeano, Gonzalez-Manteiga, 2008)  
#dt=foutliers(yf,method ="depth.trim",trim=out.rate)
    
dt<-outliers.depth.trim(Datosf,nb=200,smo=0.05,trim=out.rate,dfunc=depth.mode)
soutDT.faltan<-length(setdiff(trueout,dt$outliers))
soutDT.sobran<-length(setdiff(dt$outliers,trueout))


################ Depth Based Weighting (Febrero, Galeano, Gonzalez-Manteiga, 2008)  
#dp=foutliers(yf,method = "depth.pond",trim=out.rate)
    
dp<-outliers.depth.pond(Datosf,nb=200,smo=0.05,dfunc=depth.mode)
soutDP.faltan<-length(setdiff(trueout,dp$outliers))
soutDP.sobran<-length(setdiff(dp$outliers,trueout))


######################FOM (see Hubert et. al 2015)

yt=array(y,dim=c(dim(y),1))
Result <- fOutl(yt,  type = "fAO", diagnostic = TRUE);

# The user may opt to draw a cut off line separating the outliers.
# which will be plotted in red
ffom=fom(Result, cutoff = TRUE) 

soutFOM.faltan<-length(setdiff(trueout,which(ffom$data[,4]=='red')))
soutFOM.sobran<-length(setdiff(which(ffom$data[,4]=='red'),trueout))

######################Functional Boxplots (Sun and Genton 2011)
out.fb=fda::fbplot(y,plot=F)$outpoint 
sout.fb.faltan<-length(setdiff(trueout,out.fb))
sout.fb.sobran<-length(setdiff(out.fb,trueout))

#########################High density regions (Hydmann and Shang 2010) 
Datos=fds(t, y)
colnames(Datos$y)<-seq(1,n,1)
#fboxplot(data = Datos, plot.type = "functional", type = "hdr", alpha = c(out.rate,0.5), projmethod="PCAproj",ylab="x(t)",xlab="",xpd=NA)
hdrout<-fhdr(Datos, alpha = c(out.rate,0.5),projmethod="PCAproj",xlab='',ylab='',plotlegend=F)
sout.hdr.faltan<-length(setdiff(trueout,hdrout))
sout.hdr.sobran<-length(setdiff(hdrout,trueout))



################ SUMARY
VPr.or<-rbind(VPr.or,1-c(aa.faltan,aa.faltanr,sout.OG.faltan,sout.robMah.faltan,soutHU.faltan,soutFLO.faltan,soutDT.faltan,soutDP.faltan,soutFOM.faltan,sout.fb.faltan,sout.hdr.faltan)/length(trueout))
FPr.or<-rbind(FPr.or,c(aa.sobran,aa.sobranr,sout.OG.sobran,sout.robMah.sobran,soutHU.sobran,soutFLO.sobran,soutDT.sobran,soutDP.sobran,soutFOM.sobran,sout.fb.sobran,sout.hdr.sobran)/(n-length(trueout)))
VP.or<-rbind(VP.or,length(trueout)-c(aa.faltan,aa.faltanr,sout.OG.faltan,sout.robMah.faltan,soutHU.faltan,soutFLO.faltan,soutDT.faltan,soutDP.faltan,soutFOM.faltan,sout.fb.faltan,sout.hdr.faltan))
FP.or<-rbind(FP.or,c(aa.sobran,aa.sobranr,sout.OG.sobran,sout.robMah.sobran,soutHU.sobran,soutFLO.sobran,soutDT.sobran,soutDP.sobran,soutFOM.sobran,sout.fb.sobran,sout.hdr.sobran))

}


apply(VPr.or,2,mean)
apply(VPr.or,2,sd)
apply(FPr.or,2,mean)
apply(FPr.or,2,sd)



########################################################
#Table 2 simulated data in Section 3.2 for amplitude outliers
########################################################

out.rate=.05 


VP<-c()
FP<-c()
VPr<-c()
FPr<-c()

p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)



th=0.5
K=3
n=100
nrun=100
 norep=20


VP.or<-c()
FP.or<-c()
VPr.or<-c()
FPr.or<-c()
  
for (ite in 1:nrun) {
		cat(ite, '\n')
		set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers
trueout<-sort(sample(1:n,n*out.rate))
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

for (i in trueout){
	    y[,i]= y[,i]+2                   
}

p=length(t)
n=ncol(y)

#################################### ADA

 X=t(y)
 set.seed(1234)
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)

a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#outliers by threshold
aa.faltan<-length(setdiff(trueout,out))
aa.sobran<-length(setdiff(out,trueout))


#outliers by multivariate: Mah
sco=t(alphascce)
rownames(sco) = as.numeric(1:n)
s = cbind(sco, rep(1, n))

if( (out.rate == 0.15) & (ite == 35)){outr = robout(s, 1, "mve") # error with mcd Error in solve.default(cov, ...) 
}else{
outr = robout(s, 1, "mcd")}



aa.faltanr<-length(setdiff(trueout,outr$outliers))
aa.sobranr<-length(setdiff(outr$outliers,trueout))



################ OutlierGram (Arribas and Romo, 2014)

fD = fData( t, t(y) )
oo=outliergram( fD, display = F )
sout.OG.faltan<-length(setdiff(trueout,oo$ID_outliers))
sout.OG.sobran<-length(setdiff(oo$ID_outliers,trueout))


################ Robust Mahalanobis Distance (See Hydmann and Shang 2010)  

yf=fds(t,y)
 colnames(yf$y)=1:n

FO.robMah=foutliers(yf,method = "robMah") 
sout.robMah.faltan<-length(setdiff(trueout,FO.robMah$outliers))
sout.robMah.sobran<-length(setdiff(FO.robMah$outliers,trueout))

################ Integrated Square Error method (Hyndman and Ullah (2007); See Hydmann and Shang 2010)        
FO.HU<-foutliers(data=yf, method = "HUoutliers",order = 2, lambda = 3.29)
soutHU.faltan<-length(setdiff(trueout,FO.HU$outliers))
soutHU.sobran<-length(setdiff(FO.HU$outliers,trueout))


#################################################################
x=t(y)
rownames(x)<-as.character(1:n)

################ Package fda.usc          
Datosf<-fdata(x,t)

############################LRT (see Febrero et al. 2007)
flo=outliers.lrt(Datosf,nb=200,smo=0.05,dfunc=depth.mode,trim=out.rate)
#flo=foutliers(yf,method = "lrt",trim=out.rate) 
soutFLO.faltan<-length(setdiff(trueout,flo$outliers))
soutFLO.sobran<-length(setdiff(flo$outliers,trueout))


################ Depth Based Trimming (Febrero, Galeano, Gonzalez-Manteiga, 2008)  
#dt=foutliers(yf,method ="depth.trim",trim=out.rate)
    
dt<-outliers.depth.trim(Datosf,nb=200,smo=0.05,trim=out.rate,dfunc=depth.mode)
soutDT.faltan<-length(setdiff(trueout,dt$outliers))
soutDT.sobran<-length(setdiff(dt$outliers,trueout))


################ Depth Based Weighting (Febrero, Galeano, Gonzalez-Manteiga, 2008)  
#dp=foutliers(yf,method = "depth.pond",trim=out.rate)
    
dp<-outliers.depth.pond(Datosf,nb=200,smo=0.05,dfunc=depth.mode)
soutDP.faltan<-length(setdiff(trueout,dp$outliers))
soutDP.sobran<-length(setdiff(dp$outliers,trueout))


######################FOM (see Hubert et. al 2015)

yt=array(y,dim=c(dim(y),1))
Result <- fOutl(yt,  type = "fAO", diagnostic = TRUE);

# The user may opt to draw a cut off line separating the outliers.
# which will be plotted in red
ffom=fom(Result, cutoff = TRUE) 

soutFOM.faltan<-length(setdiff(trueout,which(ffom$data[,4]=='red')))
soutFOM.sobran<-length(setdiff(which(ffom$data[,4]=='red'),trueout))

######################Functional Boxplots (Sun and Genton 2011)
out.fb=fda::fbplot(y,plot=F)$outpoint 
sout.fb.faltan<-length(setdiff(trueout,out.fb))
sout.fb.sobran<-length(setdiff(out.fb,trueout))

#########################High density regions (Hydmann and Shang 2010) 
Datos=fds(t, y)
colnames(Datos$y)<-seq(1,n,1)
#fboxplot(data = Datos, plot.type = "functional", type = "hdr", alpha = c(out.rate,0.5), projmethod="PCAproj",ylab="x(t)",xlab="",xpd=NA)
hdrout<-fhdr(Datos, alpha = c(out.rate,0.5),projmethod="PCAproj",xlab='',ylab='',plotlegend=F)
sout.hdr.faltan<-length(setdiff(trueout,hdrout))
sout.hdr.sobran<-length(setdiff(hdrout,trueout))



################ SUMARY
VPr.or<-rbind(VPr.or,1-c(aa.faltan,aa.faltanr,sout.OG.faltan,sout.robMah.faltan,soutHU.faltan,soutFLO.faltan,soutDT.faltan,soutDP.faltan,soutFOM.faltan,sout.fb.faltan,sout.hdr.faltan)/length(trueout))
FPr.or<-rbind(FPr.or,c(aa.sobran,aa.sobranr,sout.OG.sobran,sout.robMah.sobran,soutHU.sobran,soutFLO.sobran,soutDT.sobran,soutDP.sobran,soutFOM.sobran,sout.fb.sobran,sout.hdr.sobran)/(n-length(trueout)))
VP.or<-rbind(VP.or,length(trueout)-c(aa.faltan,aa.faltanr,sout.OG.faltan,sout.robMah.faltan,soutHU.faltan,soutFLO.faltan,soutDT.faltan,soutDP.faltan,soutFOM.faltan,sout.fb.faltan,sout.hdr.faltan))
FP.or<-rbind(FP.or,c(aa.sobran,aa.sobranr,sout.OG.sobran,sout.robMah.sobran,soutHU.sobran,soutFLO.sobran,soutDT.sobran,soutDP.sobran,soutFOM.sobran,sout.fb.sobran,sout.hdr.sobran))

}


apply(VPr.or,2,mean)
apply(VPr.or,2,sd)
apply(FPr.or,2,mean)
apply(FPr.or,2,sd)



########################################################
#Table 2 simulated data in Section 3.2 for shift outliers
########################################################

out.rate=.05 


VP<-c()
FP<-c()
VPr<-c()
FPr<-c()

p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)



th=0.5
K=3
n=100
nrun=100
 norep=20


VP.or<-c()
FP.or<-c()
VPr.or<-c()
FPr.or<-c()
  
for (ite in 1:nrun) {
		cat(ite, '\n')
		set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers
trueout<-sort(sample(1:n,n*out.rate))
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

ts=seq(-0.1,0.9,len=p)

for (i in trueout){
	e=rnorm(p)
    y[,i]=  30*ts*(1-ts)^(3/2) + t(L)%*%e   
                               
}

p=length(t)
n=ncol(y)

#################################### ADA

 X=t(y)
 set.seed(1234)
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)

a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#outliers by threshold
aa.faltan<-length(setdiff(trueout,out))
aa.sobran<-length(setdiff(out,trueout))


#outliers by multivariate: Mah
sco=t(alphascce)
rownames(sco) = as.numeric(1:n)
s = cbind(sco, rep(1, n))

if( (out.rate == 0.15) & (ite == 35)){outr = robout(s, 1, "mve") # error with mcd Error in solve.default(cov, ...) 
}else{
outr = robout(s, 1, "mcd")}



aa.faltanr<-length(setdiff(trueout,outr$outliers))
aa.sobranr<-length(setdiff(outr$outliers,trueout))



################ OutlierGram (Arribas and Romo, 2014)

fD = fData( t, t(y) )
oo=outliergram( fD, display = F )
sout.OG.faltan<-length(setdiff(trueout,oo$ID_outliers))
sout.OG.sobran<-length(setdiff(oo$ID_outliers,trueout))


################ Robust Mahalanobis Distance (See Hydmann and Shang 2010)  

yf=fds(t,y)
 colnames(yf$y)=1:n

FO.robMah=foutliers(yf,method = "robMah") 
sout.robMah.faltan<-length(setdiff(trueout,FO.robMah$outliers))
sout.robMah.sobran<-length(setdiff(FO.robMah$outliers,trueout))

################ Integrated Square Error method (Hyndman and Ullah (2007); See Hydmann and Shang 2010)        
FO.HU<-foutliers(data=yf, method = "HUoutliers",order = 2, lambda = 3.29)
soutHU.faltan<-length(setdiff(trueout,FO.HU$outliers))
soutHU.sobran<-length(setdiff(FO.HU$outliers,trueout))


#################################################################
x=t(y)
rownames(x)<-as.character(1:n)

################ Package fda.usc          
Datosf<-fdata(x,t)

############################LRT (see Febrero et al. 2007)
flo=outliers.lrt(Datosf,nb=200,smo=0.05,dfunc=depth.mode,trim=out.rate)
#flo=foutliers(yf,method = "lrt",trim=out.rate) 
soutFLO.faltan<-length(setdiff(trueout,flo$outliers))
soutFLO.sobran<-length(setdiff(flo$outliers,trueout))


################ Depth Based Trimming (Febrero, Galeano, Gonzalez-Manteiga, 2008)  
#dt=foutliers(yf,method ="depth.trim",trim=out.rate)
    
dt<-outliers.depth.trim(Datosf,nb=200,smo=0.05,trim=out.rate,dfunc=depth.mode)
soutDT.faltan<-length(setdiff(trueout,dt$outliers))
soutDT.sobran<-length(setdiff(dt$outliers,trueout))


################ Depth Based Weighting (Febrero, Galeano, Gonzalez-Manteiga, 2008)  
#dp=foutliers(yf,method = "depth.pond",trim=out.rate)
    
dp<-outliers.depth.pond(Datosf,nb=200,smo=0.05,dfunc=depth.mode)
soutDP.faltan<-length(setdiff(trueout,dp$outliers))
soutDP.sobran<-length(setdiff(dp$outliers,trueout))


######################FOM (see Hubert et. al 2015)

yt=array(y,dim=c(dim(y),1))
Result <- fOutl(yt,  type = "fAO", diagnostic = TRUE);

# The user may opt to draw a cut off line separating the outliers.
# which will be plotted in red
ffom=fom(Result, cutoff = TRUE) 

soutFOM.faltan<-length(setdiff(trueout,which(ffom$data[,4]=='red')))
soutFOM.sobran<-length(setdiff(which(ffom$data[,4]=='red'),trueout))

######################Functional Boxplots (Sun and Genton 2011)
out.fb=fda::fbplot(y,plot=F)$outpoint 
sout.fb.faltan<-length(setdiff(trueout,out.fb))
sout.fb.sobran<-length(setdiff(out.fb,trueout))

#########################High density regions (Hydmann and Shang 2010) 
Datos=fds(t, y)
colnames(Datos$y)<-seq(1,n,1)
#fboxplot(data = Datos, plot.type = "functional", type = "hdr", alpha = c(out.rate,0.5), projmethod="PCAproj",ylab="x(t)",xlab="",xpd=NA)
hdrout<-fhdr(Datos, alpha = c(out.rate,0.5),projmethod="PCAproj",xlab='',ylab='',plotlegend=F)
sout.hdr.faltan<-length(setdiff(trueout,hdrout))
sout.hdr.sobran<-length(setdiff(hdrout,trueout))



################ SUMARY
VPr.or<-rbind(VPr.or,1-c(aa.faltan,aa.faltanr,sout.OG.faltan,sout.robMah.faltan,soutHU.faltan,soutFLO.faltan,soutDT.faltan,soutDP.faltan,soutFOM.faltan,sout.fb.faltan,sout.hdr.faltan)/length(trueout))
FPr.or<-rbind(FPr.or,c(aa.sobran,aa.sobranr,sout.OG.sobran,sout.robMah.sobran,soutHU.sobran,soutFLO.sobran,soutDT.sobran,soutDP.sobran,soutFOM.sobran,sout.fb.sobran,sout.hdr.sobran)/(n-length(trueout)))
VP.or<-rbind(VP.or,length(trueout)-c(aa.faltan,aa.faltanr,sout.OG.faltan,sout.robMah.faltan,soutHU.faltan,soutFLO.faltan,soutDT.faltan,soutDP.faltan,soutFOM.faltan,sout.fb.faltan,sout.hdr.faltan))
FP.or<-rbind(FP.or,c(aa.sobran,aa.sobranr,sout.OG.sobran,sout.robMah.sobran,soutHU.sobran,soutFLO.sobran,soutDT.sobran,soutDP.sobran,soutFOM.sobran,sout.fb.sobran,sout.hdr.sobran))

}


apply(VPr.or,2,mean)
apply(VPr.or,2,sd)
apply(FPr.or,2,mean)
apply(FPr.or,2,sd)


########################################################
#Table 2 simulated data in Section 3.2 for isolated outliers
########################################################

out.rate=.05 


VP<-c()
FP<-c()
VPr<-c()
FPr<-c()

p=50
t=seq(0,1,len=p)
d=dist(t,upper=T,diag=T)
d.matrix=as.matrix(d)

#Model 1
t.cov=cov.fun2(d.matrix,0.3,1/0.3)  #covariance function in time
L=chol(t.cov)  # Cholesky Decomposition
mu=30*t*(1-t)^(3/2)



th=0.5
K=3
n=100
nrun=102
 norep=20


VP.or<-c()
FP.or<-c()
VPr.or<-c()
FPr.or<-c()
  
for (ite in 1:nrun) {
		cat(ite, '\n')
		set.seed(ite)
### Generate Data
e=matrix(rnorm(n*p),p,n)
y=mu+t(L)%*%e    

### Generate Outliers
trueout<-sort(sample(1:n,n*out.rate))
t.cov=cov.fun(d.matrix,0.5,1,1)  #covariance function in time
L=chol(t.cov) # Cholesky Decomposition     

for (i in trueout){
	
    y[,i]=  y[,i]+10*c(dnorm(seq(-2,2,.3)),rep(0,36))
                               
}
p=length(t)
n=ncol(y)

#################################### ADA

 X=t(y)
 set.seed(1234)
 lass10d <- stepLArchetypoids3(data=X, k=3, norep)

a3=t(y[,lass10d[[1]][[1]]])
data=t(y)
huge=200
k=3 #

  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  zs <- rbind(t(a3), rep(huge, k))  
  zs <- as.matrix(zs)
  alphascce <- matrix(0, nrow = k, ncol = n)
  
  for (j in 1 : n){
   alphascce[, j] = coef(nnls(zs, x_gvv[,j]))
  }

wm= which.min(c(sum(alphascce[1,]>th[1]),sum(alphascce[2,]>th[1]),sum(alphascce[3,]>th[1])))
out= which(alphascce[wm,]>th[1])

#outliers by threshold
aa.faltan<-length(setdiff(trueout,out))
aa.sobran<-length(setdiff(out,trueout))


#outliers by multivariate: Mah
sco=t(alphascce)
rownames(sco) = as.numeric(1:n)
s = cbind(sco, rep(1, n))

if(  (ite == 22) || (ite==85)){outr=c()
outr$outliers=NA #outr = robout(s, 1, "mve") # error with mcd and mve Error in solve.default(cov, ...) 
}else{
outr = robout(s, 1, "mcd")}



aa.faltanr<-length(setdiff(trueout,outr$outliers))
aa.sobranr<-length(setdiff(outr$outliers,trueout))



################ OutlierGram (Arribas and Romo, 2014)

fD = fData( t, t(y) )
oo=outliergram( fD, display = F )
sout.OG.faltan<-length(setdiff(trueout,oo$ID_outliers))
sout.OG.sobran<-length(setdiff(oo$ID_outliers,trueout))


################ Robust Mahalanobis Distance (See Hydmann and Shang 2010)  

yf=fds(t,y)
 colnames(yf$y)=1:n

FO.robMah=foutliers(yf,method = "robMah") 
sout.robMah.faltan<-length(setdiff(trueout,FO.robMah$outliers))
sout.robMah.sobran<-length(setdiff(FO.robMah$outliers,trueout))

################ Integrated Square Error method (Hyndman and Ullah (2007); See Hydmann and Shang 2010)        
FO.HU<-foutliers(data=yf, method = "HUoutliers",order = 2, lambda = 3.29)
soutHU.faltan<-length(setdiff(trueout,FO.HU$outliers))
soutHU.sobran<-length(setdiff(FO.HU$outliers,trueout))


#################################################################
x=t(y)
rownames(x)<-as.character(1:n)

################ Package fda.usc          
Datosf<-fdata(x,t)

############################LRT (see Febrero et al. 2007)
flo=outliers.lrt(Datosf,nb=200,smo=0.05,dfunc=depth.mode,trim=out.rate)
#flo=foutliers(yf,method = "lrt",trim=out.rate) 
soutFLO.faltan<-length(setdiff(trueout,flo$outliers))
soutFLO.sobran<-length(setdiff(flo$outliers,trueout))


################ Depth Based Trimming (Febrero, Galeano, Gonzalez-Manteiga, 2008)  
#dt=foutliers(yf,method ="depth.trim",trim=out.rate)
    
dt<-outliers.depth.trim(Datosf,nb=200,smo=0.05,trim=out.rate,dfunc=depth.mode)
soutDT.faltan<-length(setdiff(trueout,dt$outliers))
soutDT.sobran<-length(setdiff(dt$outliers,trueout))


################ Depth Based Weighting (Febrero, Galeano, Gonzalez-Manteiga, 2008)  
#dp=foutliers(yf,method = "depth.pond",trim=out.rate)
    
dp<-outliers.depth.pond(Datosf,nb=200,smo=0.05,dfunc=depth.mode)
soutDP.faltan<-length(setdiff(trueout,dp$outliers))
soutDP.sobran<-length(setdiff(dp$outliers,trueout))


######################FOM (see Hubert et. al 2015)

yt=array(y,dim=c(dim(y),1))
Result <- fOutl(yt,  type = "fAO", diagnostic = TRUE);

# The user may opt to draw a cut off line separating the outliers.
# which will be plotted in red
ffom=fom(Result, cutoff = TRUE) 

soutFOM.faltan<-length(setdiff(trueout,which(ffom$data[,4]=='red')))
soutFOM.sobran<-length(setdiff(which(ffom$data[,4]=='red'),trueout))

######################Functional Boxplots (Sun and Genton 2011)
out.fb=fda::fbplot(y,plot=F)$outpoint 
sout.fb.faltan<-length(setdiff(trueout,out.fb))
sout.fb.sobran<-length(setdiff(out.fb,trueout))

#########################High density regions (Hydmann and Shang 2010) 
Datos=fds(t, y)
colnames(Datos$y)<-seq(1,n,1)
#fboxplot(data = Datos, plot.type = "functional", type = "hdr", alpha = c(out.rate,0.5), projmethod="PCAproj",ylab="x(t)",xlab="",xpd=NA)
hdrout<-fhdr(Datos, alpha = c(out.rate,0.5),projmethod="PCAproj",xlab='',ylab='',plotlegend=F)
sout.hdr.faltan<-length(setdiff(trueout,hdrout))
sout.hdr.sobran<-length(setdiff(hdrout,trueout))



################ SUMARY
VPr.or<-rbind(VPr.or,1-c(aa.faltan,aa.faltanr,sout.OG.faltan,sout.robMah.faltan,soutHU.faltan,soutFLO.faltan,soutDT.faltan,soutDP.faltan,soutFOM.faltan,sout.fb.faltan,sout.hdr.faltan)/length(trueout))
FPr.or<-rbind(FPr.or,c(aa.sobran,aa.sobranr,sout.OG.sobran,sout.robMah.sobran,soutHU.sobran,soutFLO.sobran,soutDT.sobran,soutDP.sobran,soutFOM.sobran,sout.fb.sobran,sout.hdr.sobran)/(n-length(trueout)))
VP.or<-rbind(VP.or,length(trueout)-c(aa.faltan,aa.faltanr,sout.OG.faltan,sout.robMah.faltan,soutHU.faltan,soutFLO.faltan,soutDT.faltan,soutDP.faltan,soutFOM.faltan,sout.fb.faltan,sout.hdr.faltan))
FP.or<-rbind(FP.or,c(aa.sobran,aa.sobranr,sout.OG.sobran,sout.robMah.sobran,soutHU.sobran,soutFLO.sobran,soutDT.sobran,soutDP.sobran,soutFOM.sobran,sout.fb.sobran,sout.hdr.sobran))

}


VPr.or=VPr.or[c(-22,-85),]
FPr.or=FPr.or[c(-22,-85),]


apply(VPr.or,2,mean)
apply(VPr.or,2,sd)
apply(FPr.or,2,mean)
apply(FPr.or,2,sd)

#ite 22 and 85 removed,  robout error


apply(VPr.or,2,mean)
apply(VPr.or,2,sd)
apply(FPr.or,2,mean)
apply(FPr.or,2,sd)

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