



#Arribas and Romo 2014

# 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)}  





out.rate=.15 #0.1



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)


po=0.5 #0.25, 0.75 #change the parameters for obtaining the results of Table 3
ct=6 #ct=1
th=0.5
K=2
n=100
nrun=100
 norep=20


S=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: use functions in Directory Archetypal functions

 X=t(y)
 set.seed(1234)
 lassd <- stepLArchetypoids3(data=X, k=2, norep) # no robust

out=lassd[[1]][[1]]

 set.seed(1234)
 lassdr <- stepRArchetypoids3(data=X, k=2, norep,PO=po,CT=ct) #  robust
# lassdr <- stepRArchetypoids3(data=X, k=2, norep,PO=0.75,CT=1) #  robust
outr=lassdr[[1]][[1]]



#outliers in ADAs solution
anr<-sum(is.element(out,trueout))

ar<-sum(is.element(outr,trueout))


################ SUMARY
S<-rbind(S,c(anr,ar))

}


apply(S,2,mean)


############################### Results
#out.rate=.15

 apply(S,2,mean)
#[1] 0.78 0.32  #for c=6 me

############################ Results
#out.rate=.1


# apply(S,2,mean)
#[1] 0.1 0.0  #for c=6 me

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

#Comparative robust PCA

library(rospca)

#################
#Figure 2
################

ite=17 #7
out.rate=.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)
 lassdr <- stepRArchetypoids3(data=X, k=2, norep,PO=0.5,CT=6) #  robust

outr=lassdr[[1]][[1]]

a2a=t(y[,lassdr[[1]][[1]]])





 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
for (i in 1:10){
lines(t, y[,trueout[i]], 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, a2a[1,], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #robust archetypoids
lines(t, a2a[2,], type = "l",lty=2,xlab="",ylab="",col=3,lwd=3) #robust archetypoids



windows()

 set.seed(1234)
 lassd <- stepLArchetypoids3(data=X, k=2, norep) # no robust

out=lassd[[1]][[1]]


a2a=t(y[,lassd[[1]][[1]]])



 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
for (i in 1:10){
lines(t, y[,trueout[i]], 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, a2a[1,], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #no robust archetypoids
lines(t, a2a[2,], type = "l",lty=2,xlab="",ylab="",col=3,lwd=3) #no robust archetypoids


#PCA

windows()

p=princomp(X)


 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
for (i in 1:10){
lines(t, y[,trueout[i]], 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, apply(X,2,mean)+5*p$loadings[,1], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #pc+ c1
lines(t, apply(X,2,mean)-5*p$loadings[,1], type = "l",lty=2,xlab="",ylab="",col=3,lwd=3) #pc-c1


#robust PCA


set.seed(1234)
resR <- robpca(X, k=2)

windows()


 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
for (i in 1:10){
lines(t, y[,trueout[i]], 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, resR$center+5*resR$loadings[,1], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #pc+ c1
lines(t, resR$center-5*resR$loadings[,1], type = "l",lty=2,xlab="",ylab="",col=6,lwd=3) #pc-c1
lines(t, resR$center+2*resR$loadings[,1], type = "l",lty=3,xlab="",ylab="",col=4,lwd=3) #pc+ c1
lines(t, resR$center-2*resR$loadings[,1], type = "l",lty=4,xlab="",ylab="",col=5,lwd=3) #pc-c1

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

#Table 4

library(EvaluationMeasures)

library(rainbow)
library(rospca)



#Arribas and Romo 2014

# 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)}  





out.rate=.1 #0.15, 0.05, 0



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=2
n=100
nrun=100
 norep=20



S=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)

reference=rep(FALSE,n)
reference[trueout]=TRUE


 X=t(y)

set.seed(1234)
resR <- robpca(X, k=2)

#robust pca
 tprrpca=EvaluationMeasures.TPR(reference,!resR$flag.all)
prpca=EvaluationMeasures.MCC(reference,!resR$flag.all)
fprrpca= EvaluationMeasures.FPR(reference,!resR$flag.all)





#################################### ADA: use functions in Directory Archetypal functions


 set.seed(1234)
 lassdr <- stepRArchetypoids3(data=X, k=2, norep,PO=0.5,CT=1) #  robust
outr=lassdr[[1]][[1]]



#outliers in ADAs solution

ar<-sum(is.element(outr,trueout))
huge=200

data=X
k=2
  archs_ini <- outr
  
  n <- ncol(t(data))
  x_gvv <- rbind(t(data), 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]))
  }
  
  #resid <- zs %*% alphas - x_gvv
  resid <- zs[1: (nrow(zs)-1),] %*% alphas - x_gvv[1: (nrow(x_gvv)-1),]
  
r=apply(resid,2,intpms)

re=r

b=boxplot(re)
ooaa= which(re %in% b$out)

oada=rep(FALSE,n)
oada[ooaa]=TRUE


tprrada=EvaluationMeasures.TPR(reference,oada)
prada=EvaluationMeasures.MCC(reference,oada)
fprrada= EvaluationMeasures.FPR(reference,oada)


################ Integrated Square Error method (Hyndman and Ullah (2007); See Hydmann and Shang 2010)  


yf=fds(t,y)
 colnames(yf$y)=1:n
      
FO.HU<-foutliers(data=yf, method = "HUoutliers",order = 2, lambda = 3.29)


hh=rep(FALSE,n)
hh[FO.HU$outliers]=TRUE


tprrh=EvaluationMeasures.TPR(reference,hh)
prh=EvaluationMeasures.MCC(reference,hh)
fprrh= EvaluationMeasures.FPR(reference,hh)




################ SUMARY
S<-rbind(S,c(tprrpca,prpca,fprrpca,tprrada,prada,fprrada,tprrh,prh,fprrh))

}


apply(S,2,summary)*100
apply(S,2,mean)*100
apply(S,2,sd)*100