## 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("MstepArchetypesMod.R")


X=rbind(c(1,2),c(5,NaN),c(6,6),c(4,7))
plot(X)


norep=20
set.seed(2016)
lass10 <- MstepArchetypesMod(data=X,k=1:2,verbose=FALSE,nrep=norep,saveHistory=F) 
 screeplot(lass10)


numarq=1; 
 
a=bestModel(lass10[[numarq]])
 parameters(a) #the mean

#[1,] 4.003687 5.004524


numarq=2; 
 
a=bestModel(lass10[[numarq]])
 parameters(a) 

#Fig 1 supplement
plot(X,ylab='',xlab='')
points(parameters(a),col=2,pch=8)
points(t(matrix(c(1,2,5,6.46),ncol=2)),col=4,pch=2) #from toy.m

#[1,] 1.051004 2.633046
#[2,] 5.468799 6.227073

set.seed(2016)
lass10 <- MstepArchetypesMod(data=X[-2,],k=1:2,verbose=FALSE,nrep=norep,saveHistory=F) 
 screeplot(lass10)


numarq=2; 
 
a=bestModel(lass10[[numarq]])
 parameters(a) 

         [,1]     [,2]
#[1,] 5.065543 6.464278
#[2,] 1.002196 2.003660


numarq=1; 
 
a=bestModel(lass10[[numarq]])
 parameters(a) 

# 3.669342 5.003415

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

###################
 #require(nnls)
nnls.alphasfn <- function(coefs, C, d, ...) {
 

  n = ncol(d)

  for ( j in 1:n )
    coefs[,j] = coef(nnls(C, d[,j]))

  return(coefs)
}



  bp.dummyfn <- function(x, ...) {
	huge=200
    y = rbind(x, rep(huge, ncol(x)))

    attr(y, '.Meta') = attr(x, '.Meta')
    attr(y, '.Meta')$dummyrow = nrow(y)

    return(y)
  }

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

#Example in a well-known data set

#In shapes the missing landmarks could not be completely random but, for example, extremities

#So, we will make NA in 50\% but in specific variables


library(mlbench)
R=matrix(0,ncol=4,nrow=100)

for (i in 1:100){
pe=0.5 #
n=150
p=21

#set.seed(2017)
set.seed(i)

datag=mlbench.waveform(150)


wn1=sample(1:n,pe*n) #150 specimens

data <- datag$x
data[wn1,5]=NaN #NaNs in variable number 5

wn2=sample(1:n,pe*n) #150 specimens
data[wn2,10]=NaN #NaNs in variable number 10

wn3=sample(1:n,pe*n) #150 specimens
data[wn3,15]=NaN #NaNs in variable number 15

x1 <- bp.dummyfn(t(datag$x))







#AA: 3 archetypes: original
norep=20
set.seed(2017)
lass10 <- LstepArchetypesMod(datag$x,k=1:3,verbose=FALSE,nrep=norep,saveHistory=F) 
 screeplot(lass10)


numarq=3; 
 
a=bestModel(lass10[[numarq]])
zso=parameters(a) 

ao= matrix(0, nrow = 3, ncol = n)  #alphas estimated
ao <- nnls.alphasfn(ao, bp.dummyfn(t(zso)), x1)

resid <- t(zso) %*% ao  - t(datag$x)

rsso=norm2.Frnormfn(resid)/n



##AA missing
norep=20
 set.seed(1234)

lass10 <- MstepArchetypesMod(data,k=1:3,verbose=FALSE,nrep=norep,saveHistory=F) 
 screeplot(lass10)


numarq=3; 
 
a=bestModel(lass10[[numarq]])
zs=parameters(a) 

ae= matrix(0, nrow = 3, ncol = n)  #alphas estimated
ae <- nnls.alphasfn(ae, bp.dummyfn(t(zs)), x1)

resid <- t(zs) %*% ae  - t(datag$x)

rsse=norm2.Frnormfn(resid)/n


############
#If only work with complete cases


norep=20
sum(complete.cases(data))
set.seed(2017)
lass10 <- LstepArchetypesMod(data[complete.cases(data),],k=1:3,verbose=FALSE,nrep=norep,saveHistory=F) 
 screeplot(lass10)


numarq=3; 
 
a=bestModel(lass10[[numarq]])
zsco=parameters(a) 

aco= matrix(0, nrow = 3, ncol = n)  #alphas estimated
aco <- nnls.alphasfn(aco, bp.dummyfn(t(zsco)), x1)

resid <- t(zsco) %*% aco  - t(datag$x)

rssco=norm2.Frnormfn(resid)/n

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

library(Hmisc)

dataf=data.frame(data)

if(i==28 | i==46 | i==51 | i==54 | i==56 | i==71 | i==88 | i== 94| i ==100){ #it gives error with seed(2017) 
set.seed(1)
 dimpute <- aregImpute(~ X1+X2+X3+X4+X5+X6+X7+X8+X9+X10+X11+X12+X13+X14+X15+X16+X17+X18+X19+X20+X21, data = dataf)

}else{
set.seed(2017)
 dimpute <- aregImpute(~ X1+X2+X3+X4+X5+X6+X7+X8+X9+X10+X11+X12+X13+X14+X15+X16+X17+X18+X19+X20+X21, data = dataf)
}

#combine de data #https://stats.stackexchange.com/questions/63313/how-to-combine-multiple-imputed-datasets

d1=dataf
d1[which(is.nan(dataf$X5)),5]=dimpute$imputed$X5[1]
d1[which(is.nan(dataf$X10)),10]=dimpute$imputed$X10[1]
d1[which(is.nan(dataf$X15)),15]=dimpute$imputed$X15[1]

cdata=d1

d2=dataf
d2[which(is.nan(dataf$X5)),5]=dimpute$imputed$X5[2]
d2[which(is.nan(dataf$X10)),10]=dimpute$imputed$X10[2]
d2[which(is.nan(dataf$X15)),15]=dimpute$imputed$X15[2]

cdata=rbind(cdata,d2)


d3=dataf
d3[which(is.nan(dataf$X5)),5]=dimpute$imputed$X5[3]
d3[which(is.nan(dataf$X10)),10]=dimpute$imputed$X10[3]
d3[which(is.nan(dataf$X15)),15]=dimpute$imputed$X15[3]

cdata=rbind(cdata,d3)


d4=dataf
d4[which(is.nan(dataf$X5)),5]=dimpute$imputed$X5[4]
d4[which(is.nan(dataf$X10)),10]=dimpute$imputed$X10[4]
d4[which(is.nan(dataf$X15)),15]=dimpute$imputed$X15[4]

cdata=rbind(cdata,d4)


d5=dataf
d5[which(is.nan(dataf$X5)),5]=dimpute$imputed$X5[5]
d5[which(is.nan(dataf$X10)),10]=dimpute$imputed$X10[5]
d5[which(is.nan(dataf$X15)),15]=dimpute$imputed$X15[5]

cdata=rbind(cdata,d5)


#AA: 3 archetypes:imputed Hmisc
norep=20
set.seed(2017)
lass10 <- LstepArchetypesMod(cdata,k=1:3,verbose=FALSE,nrep=norep,saveHistory=F) 
 screeplot(lass10)


numarq=3; 
 
a=bestModel(lass10[[numarq]])
zsi=parameters(a) 

ai= matrix(0, nrow = 3, ncol = n)  #alphas estimated
ai <- nnls.alphasfn(ai, bp.dummyfn(t(zsi)), x1)

resid <- t(zsi) %*% ai  - t(datag$x)

rssi=norm2.Frnormfn(resid)/n


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

R[i,]=c(rsso,rsse,rssi,rssco)

cat(i)
}

#save(R,file='Rmatrix.Rda')



