


source("FstepArchetypesModB.R")
source("FstepArchetypoidsB.R")
source("FFunctions_to_calculate_real_archetypes_with_swap2B.R")

library(archetypes)
library(fda)


data1 <- read.table('Fertilitynonames.txt',header=T,sep='\t',dec=',',na.strings='..',row.names=1)


codnoc=c(7,35,37,55,56,64,65,66,67,70,86,87,88,89,112,113,115,121,122,123,137,138,139,154,156,
158,159,181,185,191,192,213,220)

data1=data1[-codnoc,]


colnames(data1)=1960:2013


data2 <- read.table('totalexpnonames.txt',header=T,sep='\t',dec=',',na.strings='..',row.names=1)



colnames(data2)=1960:2013

data2=data2[-codnoc,]

data=array(0,dim=c(dim(data1),2))
data[,,1]=as.matrix(data1)
data[,,2]=as.matrix(data2)

colnames(data)=1960:2013
rownames(data)=rownames(data1)

#  ---------------  argument vectors  --------------------

years <- seq(1960, 2013, 1)    #  


#  -------------------  set up the  basis  ----------------------

#    

yearbasisfd <- create.bspline.basis( c(1960,2013), 32)

#  --------------  make fd object  ------------

tdata=array(0,dim=c(dim(t(data1)),2))
tdata[,,1]=as.matrix(t(data1))
tdata[,,2]=as.matrix(t(data2))

rownames(tdata)=1960:2013
colnames(tdata)=rownames(data1)

tempfd <- Data2fd(years,tdata, yearbasisfd)


#tempfdLUX <- Data2fd(years[c(-2,-4)],as.matrix(t(data[104,c(-2,-4)])), yearbasisfd)
tempfdLUX =Data2fd(years[c(-2,-4)],tdata[c(-2,-4),104,], yearbasisfd)


#Singapore 1 : 2000 year
#tempfdSGP <- Data2fd(years[c(-41)],as.matrix(t(data[178,c(-41)])), yearbasisfd)
#tempfdSGP <- Data2fd(years[c(-41)],as.matrix(t(data[151,c(-41)])), yearbasisfd)

tempfdSGP=Data2fd(years[c(-41)],tdata[c(-41),151,], yearbasisfd)

#Israel 4 missing values: 1960 and 1963, 1964, 1965

#tempfdISR <- Data2fd(years[c(-1,-4,-5,-6)],tdata[c(-1,-4,-5,-6),110,], yearbasisfd)
tempfdISR <- Data2fd(years[c(-1,-4,-5,-6)],tdata[c(-1,-4,-5,-6),85,], yearbasisfd)


#plot(tempfd[c(-124,-178),]) #mirar ajuste
# plotfit.fd(t(data[c(-124,-178),]), years, tempfd[c(-124,-178),])



#dim(tempfd$coefs)[1]  32 223   2
# dim(tempfdISR$coefs)[1] 32  2





tempfd$coefs[,85,2]=tempfdISR$coefs[,2]
tempfd$coefs[,104,1]=tempfdLUX$coefs[,1]
tempfd$coefs[,151,1]=tempfdSGP$coefs[,1]


me=mean.fd(tempfd)
#va=var.fd(tempfd)
s=sd.fd(tempfd)


X=array(0,dim=c(dim(t(tempfd$coef[,,1])),2))
X[,,1]=t(tempfd$coef[,,1])
X[,,2]=t(tempfd$coef[,,2])


X[85,,2]=tempfdISR$coefs[,2]
X[104,,1]=tempfdLUX$coefs[,1]
X[151,,1]=tempfdSGP$coefs[,1]

Xs=X
Xs[,,1]=scale(X[,,1])
Xs[,,2]=scale(X[,,2])

PM=eval.penalty(yearbasisfd, rng=c(1960,2013))

K=5 #
set.seed(2015)

norep=20
lass10 <- FstepArchetypesModB(data=Xs,k=K,verbose=FALSE,nrep=norep,PM=PM,saveHistory=F) 
ai <- bestModel(lass10[[1]]) #archetypes
 
an4 <-  FstepArchetypoidsB(K,TRUE,Xs,lass10,FALSE,K-1,PM)

aw4 <-  FstepArchetypoidsB(K,FALSE,Xs,lass10,FALSE,K-1,PM) #which_alpha archetypoids

ini_arch=c()
 for (j in 1:K){
 ini_arch[j]=which.max(ai$betas[j,])
 }
 
  x11 <- t(Xs[,,1])
    x12 <- t(Xs[,,2])
    x1=rbind(x11,x12)
    data=t(x1)
 
abeta <- FrealArchetypesB(data,huge = 200, K, ini_arch,PM) #which beta archetypoid

#alphas of an4, aw4 or abeta


x11 <- t(Xs[,,1])
    x12 <- t(Xs[,,2])
    x1=rbind(x11,x12)
    data=t(x1)
    
huge=200
 archs_ini <- an4[[1]]
 data=data
  
  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]))
  }

colnames(alphas)=rownames(data1)


matplot(1960:2013,t(as.matrix(data1[c(an4[[1]]),])))
windows()
matplot(1960:2013,t(as.matrix(data2[c(an4[[1]]),])))



#figures
matplot(1960:2013,t(as.matrix(data1[c(an4[[1]]),])),cex.lab=1.5,cex.axis=1.5,xlab="Years",ylab="TFR",cex=1.2)
windows()
matplot(1960:2013,t(as.matrix(data2[c(an4[[1]]),])),cex.lab=1.5,cex.axis=1.5,xlab="Years",ylab="LEB",cex=1.2)




x11 <- t(Xs[,,1])
    x12 <- t(Xs[,,2])
    x1=rbind(x11,x12)
    data=t(x1)
    
huge=200 
 archs_ini <- an4[[1]]
 data=data
  
  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]))
  }

colnames(alphas)=rownames(data1)



alphas[alphas>1]=1
alphas[alphas<0]=0

countnames <- read.csv2('countries2.txt',header=T)

countnames=countnames[-codnoc,]

countnames[an4[[1]]]

 data1[an4[[1]],]
  data2[an4[[1]],]

 t(data1[an4[[1]],])
  t(data2[an4[[1]],])


library(rworldmap)



da=data.frame(code=rownames(data1),cn=countnames,a=t(alphas))


#da=da[-codnoc,] #190 countries

 sPDF2 <- joinCountryData2Map(da, joinCode = "NAME", nameJoinColumn = "cn",verbose=T) #falla channel island porque su libreria no le pone codigo
 
 
  mapCountryData( sPDF2, nameColumnToPlot="a.1",catMethod=seq(0,1,0.1),missingCountryCol=3,mapTitle = 'Archetypoid 1') #
 mapCountryData( sPDF2, nameColumnToPlot="a.1",catMethod=seq(0,1,0.1),missingCountryCol=3,mapTitle = '',aspect='variable') #
 windows()
  mapCountryData( sPDF2, nameColumnToPlot="a.2",catMethod=seq(0,1,0.1),missingCountryCol=3,mapTitle = 'Archetypoid 2') # Japon tiene alpha >1 no represntaria, debrķa cambiar huge en programa!
  mapCountryData( sPDF2, nameColumnToPlot="a.2",catMethod=seq(0,1,0.1),missingCountryCol=3,mapTitle = '',aspect='variable') # 
 
 windows()
  mapCountryData( sPDF2, nameColumnToPlot="a.3",catMethod=seq(0,1,0.1),missingCountryCol=3,mapTitle = 'Archetypoid 3') #
   mapCountryData( sPDF2, nameColumnToPlot="a.3",catMethod=seq(0,1,0.1),missingCountryCol=3,mapTitle = '',aspect='variable') #
  
 windows()
  mapCountryData( sPDF2, nameColumnToPlot="a.4",catMethod=seq(0,1,0.1),missingCountryCol=3,mapTitle = 'Archetypoid 4') #
   mapCountryData( sPDF2, nameColumnToPlot="a.4",catMethod=seq(0,1,0.1),missingCountryCol=3,mapTitle = '',aspect='variable') #
  
 windows()
  mapCountryData( sPDF2, nameColumnToPlot="a.5",catMethod=seq(0,1.,0.1),missingCountryCol=3,mapTitle = 'Archetypoid 5') #
   mapCountryData( sPDF2, nameColumnToPlot="a.5",catMethod=seq(0,1,0.1),missingCountryCol=3,mapTitle = '',aspect='variable') #
 
 
 
 

########################
#removing countries  with conflict: Sierra Leone, Cambodia and Rwanda



Xs=X[c(-31,-150,-145),,]
Xs[,,1]=scale(Xs[,,1])
Xs[,,2]=scale(Xs[,,2])

K=5
set.seed(2015)

norep=20
lass10 <- FstepArchetypesModB(data=Xs,k=K,verbose=FALSE,nrep=norep,PM=PM,saveHistory=F) 
ai <- bestModel(lass10[[1]]) #archetypes
 
an4 <-  FstepArchetypoidsB(K,TRUE,Xs,lass10,FALSE,K-1,PM)

aw4 <-  FstepArchetypoidsB(K,FALSE,Xs,lass10,FALSE,K-1,PM) #which_alpha archetypoids

ini_arch=c()
 for (j in 1:K){
 ini_arch[j]=which.max(ai$betas[j,])
 }
 
  x11 <- t(Xs[,,1])
    x12 <- t(Xs[,,2])
    x1=rbind(x11,x12)
    data=t(x1)
 
abeta <- FrealArchetypesB(data,huge = 200, K, ini_arch,PM) #which beta archetypoid


countnames <- read.csv2('countries2.txt',header=T)

countnames=countnames[-codnoc,]

countnames=countnames[c(-31,-150,-145)]

countnames[an4[[1]]]


d1=data1[c(-31,-150,-145),]

d1[an4[[1]],]


d2=data2[c(-31,-150,-145),]

d2[an4[[1]],]


###############################################
#some extra results, without erasing those three countries, with k=5

#Residuals

 n <- ncol(t(data))
zs <- x_gvv[, archs_ini] 
resid <- zs[1: (nrow(zs)-1),] %*% alphas - x_gvv[1: (nrow(x_gvv)-1),]

rssm<-c()
di=dim(resid)

for(i in 1:n){
ri1=resid[1:(di[1]/2),i]
ri2=resid[((di[1]/2 + 1): di[1]),i]
rssm[i]=(t(ri1)%*%PM%*%ri1) + (t(ri2)%*%PM%*%ri2) 
}

countnames[order(rssm)]
#####################################

#how many countires for each archetypoid with high alphas > 0.7

countnames[alphas[1,]>0.7]
countnames[alphas[2,]>0.7]
countnames[alphas[3,]>0.7]
countnames[alphas[4,]>0.7]
countnames[alphas[5,]>0.7]

#####################################
#Fig. 3 with the rest of countries in gray

matplot(1960:2013,t(as.matrix(data1[c(an4[[1]]),])),cex.lab=1.5,cex.axis=1.5,xlab="Years",ylab="TFR",cex=1.2,ylim=c(1,9),lty=1:5,type='l',lwd=2)
wo=1:dim(data1)[1]
wo=wo[-an4[[1]]]
for (i in wo){
lines(1960:2013,t(as.matrix(data1[i,])),col=8,cex.lab=1.5,cex.axis=1.5,xlab="Years",ylab="TFR",cex=1.2)}
windows()
matplot(1960:2013,t(as.matrix(data2[c(an4[[1]]),])),cex.lab=1.5,cex.axis=1.5,xlab="Years",ylab="LEB",cex=1.2,ylim=c(20,85),lty=1:5,type='l',lwd=2)
for (i in wo){
lines(1960:2013,t(as.matrix(data2[i,])),col=8,cex.lab=1.5,cex.axis=1.5,xlab="Years",ylab="TFR",cex=1.2)}

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

#Panorama plot modified

tempfds=tempfd #standardized functions
tempfds$coefs[,,1]=t(Xs[,,1])
tempfds$coefs[,,2]=t(Xs[,,2])

 et=eval.fd(years,tempfds) #evaluate fd in years
 
 library(fda.usc)
 #compute the lp distance between bivariate functions et and archetypoids
 
 
 d=metric.lp(t(et[,,1]),t(et[,an4[[1]],1])) + metric.lp(t(et[,,2]),t(et[,an4[[1]],2]))
 
 #archetypoid 1
 i=1 #change the number of archetypoid each time
   ix <- order(d[, i])
 plot(1:dim(d)[1], d[ix, i],xlab='',ylab='',col=8)
 
 
 atypes.col=(seq(length = K) + 1)
 atypes.pch = (seq(length = K) - 1)
 
 atypes.cex = rep(1.5, K)
 
points(match(an4[[1]],ix), d[an4[[1]], i],col=atypes.col,pch= atypes.pch,cex=atypes.cex,lwd=3)

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

#Pentagon (simplex) plot

simplexplotm=
function (alphas,arch, k=K, radius = 10, order = NULL, labels_cex = 1, 
    labels = NULL, show_labels = TRUE, points_col = "#00000044", 
    points_pch = 19, points_cex = 1, projection = simplex_projection, 
    show_points = TRUE, show_circle = TRUE, circle_col = "lightgray", 
    show_edges = TRUE, edges_col = "lightgray", show_direction = FALSE, 
    direction_length = 1, directions_col = points_col, ...) 
{
    
   
    if (is.null(order)) 
        order <- 1:k
    if (is.null(labels)) 
        labels <- sprintf("A%s", order)
    if (length(points_col) == 1) 
        points_col <- rep(points_col, nrow(t(alphas)))
    if (length(points_cex) == 1) 
        points_cex <- rep(points_cex, nrow(t(alphas)))
    if (length(directions_col) == 1) 
        directions_col <- rep(directions_col, nrow(t(alphas)))
    params <- alphas
    coefs <- t(alphas)
    proj_z <- projection(params, r = radius - 1)
    proj_h <- coefs %*% proj_z
    proj_labels <- proj_z
    t <- cbind(x = acos(proj_z[, "x"]/(radius - 1)), y = asin(proj_z[, 
        "y"]/(radius - 1)))
    proj_labels <- cbind(x = radius * cos(t[, "x"]), y = radius * 
        sin(t[, "y"]))
    proj_circle <- list(center = cbind(x = 0, y = 0), radius = radius - 
        1)
    proj_edges <- proj_z[as.integer(combn(1:k, 2)), ]
    proj_directions <- vector("list", length = ncol(alphas))
    for (j in 1:ncol(alphas)) {
        s <- proj_h[j, , drop = FALSE]
        d <- matrix(NA_real_, ncol = 2, nrow = nrow(alphas))
        for (i in 1:nrow(alphas)) {
            e <- proj_z[i, , drop = FALSE]
            v <- e - s
            m <- sqrt(sum(v^2))
            v <- v/m
            px <- s[1] + v[1] * direction_length * alphas[i, 
                j]
            py <- s[2] + v[2] * direction_length * alphas[i, 
                j]
            d[i, ] <- c(px, py)
        }
        proj_directions[[j]] <- list(s = s, e = d)
    }
    plot(proj_z, type = "n", asp = TRUE, xlim = c(-radius, radius), 
        ylim = c(-radius, radius), axes = FALSE, xlab = "", ylab = "")
    if (show_circle) {
        symbols(proj_circle$center, circles = radius - 1, inches = FALSE, 
            add = TRUE, asp = TRUE, fg = circle_col)
    }
    if (show_edges) {
        lines(proj_edges, col = edges_col)
    }
    if (show_labels) {
        text(proj_labels, labels = labels, cex = labels_cex)
    }
    if (show_direction) {
        for (d in proj_directions) {
            for (i in 1:nrow(d$e)) {
                lines(rbind(d$s, d$e[i, , drop = FALSE]), col = directions_col)
            }
        }
    }
    if (show_points) {
        points(proj_h, col = points_col, pch = points_pch, cex = points_cex)
    }
    ret <- list(proj_z = proj_z, proj_h = proj_h, proj_labels = proj_labels, 
        proj_directions = proj_directions, proj_circle = proj_circle, 
        proj_edges = proj_edges)
    class(ret) <- "simplexplot"
    invisible(ret)
}

 simplex_projection=
function (x, r = 10) 
{
    phi <- seq(-pi, pi, length.out = nrow(x) + 1)
    phi <- phi[-1]
    cbind(x = r * cos(phi), y = r * sin(phi))
}



simplexplotm(alphas,an4[[1]],K,labels=countnames[an4[[1]]])



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

#CMDS with dt

dt=metric.lp(t(et[,,1]),t(et[,,1])) + metric.lp(t(et[,,2]),t(et[,,2]))
fit <- cmdscale(dt,eig=T,k=2)
x <- fit$points[,1]
y <- fit$points[,2]
plot(x, y, xlab="Coordinate 1", ylab="Coordinate 2",
  main="Metric MDS", type="n")
text(x, y, labels = rownames(data1), cex=.7) 


#GOG: standarized: $GOF [1] 0.8532811 0.8922773

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

####################
#FPCA


pt=pca.fd(tempfds)
> pt$var
[1] 0.84552774 0.08058372
> pt$var[1]+pt$var[2]
[1] 0.9261115

pt$scores

ju=pt$scores[,,1]+pt$scores[,,2] #sum int prod

x=ju[,1]
y=ju[,2]
plot(x, y, xlab="PCA 1", ylab="PCA 2",
  main="PCA scores", type="n")
text(x, y, labels = rownames(data1), cex=.7) 


##############################
#outliers: Multivariate functional outlier detection and Rejoinder: Mia Hubert, Peter J. Rousseeuw, Pieter Segaert


library(mrfDepth)


#fom
Resultf = mfd(et, alpha=0, type = "sprojdepth", crossdepth = TRUE, diagnostic = TRUE);
fom(Resultf)

library(matrixStats)

NFunc <- nrow(Resultf$crossdepth)
    NTObs <- ncol(Resultf$crossdepth)
    AOValues <- 1/Resultf$crossdepth - 1
    fao <- AOValues %*% Resultf$weights
    vao=rowSds(AOValues)/(1 + fao)

plot(fao,vao,type='n',xlab=c("fAO") , ylab=c("sd(AO) / (1+fAO)"))
text(fao,vao,labels=rownames(data1))


