
## 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/>.


no.scalefn <- function(x, ...){
  return(x)
}
no.rescalefn <- function(x, zs, ...){
  if ( is.null(zs) )
    return(matrix(NA, nrow = 0, ncol = 0))
  return(zs)
}


make.dummyfn <- function(huge=200) {

  bp.dummyfn <- function(x) {
    y = rbind(x, rep(huge, ncol(x))) 
    
    attr(y, '.Meta') = attr(x, '.Meta')
    attr(y, '.Meta')$dummyrow = nrow(y)
  
    return(y)
  }

  return(bp.dummyfn)
}







bisquaref <- function(resid,PO,CT, ...) {
 
  resid0 <- resid < sqrt(.Machine$double.eps)

  #c <- 6 * median(resid[!resid0])
	c <- CT * quantile(resid[!resid0],probs=PO)
  v <- resid / c
	co=c^2/6

  ifelse(resid <= c, co*(1-(1 - v^2)^3), co)
}



intpm<-function(ri,PM){
(t(ri)%*%PM%*%ri)}



norm2.RFFnormfn2 <- function(m, PM,PO,CT) {

di=dim(m)
r1=apply(m[1:(di[1]/2),],2,intpm, PM=PM)
r2=apply(m[((di[1]/2 + 1): di[1]),],2,intpm, PM=PM)
r=sqrt(r1+r2)
return(sum(bisquaref(r,PO,CT)))


}




library(FNN) #for the function knn.


RFstepArchetypoidsB <- function(i,nearest,data,ArchObj,sequ,aux,PM,PO,CT){


#datao=data
  x11 <- t(data[,,1])
    x12 <- t(data[,,2])
    x1=rbind(x11,x12)
    data=t(x1)
  
  N = dim(data)[1]
  
  if(sequ){
   ai <- bestModel(ArchObj[[i]])
  }else{
   ai <- bestModel(ArchObj[[i-aux]])  
  }
  
  if(is.null(parameters(ai))){
    stop("No archetypes computed")  
  }
  
    
  if(nearest){
  
    	dime=dim(parameters(ai))
  	ar=parameters(ai)
  	R=matrix(0,nrow=dime[1], ncol=N)
  	for(di in 1:dime[1]){ #dime[1] is the number of archetypes
  
  		
  Raux=matrix(ar[di,],byrow=F,ncol=N,nrow=dime[2])-t(data)
  
  
   dii=dim(Raux)
  
  		
  R[di,]= apply(Raux[1:(dii[1]/2),],2,intpm, PM=PM)+apply(Raux[((dii[1]/2 + 1): dii[1]),],2,intpm, PM=PM)
  
  
  

    	}

  ini_arch <- apply(R, 1, which.min)

 
    
    if( any(duplicated(ini_arch)) == TRUE){ #for very rare ocasions, it can be improved
     # ini_arch=unique(ini_arch)      

	 k=1
       neig=apply(R, 1, rank)     
      indices1 <- neig[1,]
      ini_arch <- indices1[,k]
      
      while(any(duplicated(ini_arch))){
        k=k+1  
       
        indicesk <- neig[k,]
               
        dupl <- anyDuplicated(indices1)
        ini_arch <- c(indices1[-dupl],indicesk[dupl])
       
      }
    }
    
  }else{
    ini_arch <- apply(coef(ai, "alphas"), 2, which.max) 
  }
  
  res <- RFrealArchetypesB(data, huge = 200, i, ini_arch,PM,PO,CT) #as.vector(ini_arch)
  cat("Done!") 
  return(list(res[[1]],res[[2]],ini_arch))
}





RstepFArchetypoidsnB <- function (data, k, norep=3,PM,PO,CT){
  mycall <- match.call()
  as <- list()
  for (i in 1:length(k)) {
    as[[i]] <- list()
    class(as[[i]]) <- "repArchetypoids"
    lass <- RFstepArchetypesModB(data,k[i],verbose=FALSE,nrep=norep,PM=PM,saveHistory=F,PO,CT) 
          as[[i]]<- RFstepArchetypoidsB(k[i],TRUE,data,lass,FALSE,k[i]-1,PM,PO,CT)
     
      
    
  }
  return(structure(as, class = "stepArchetypoids", call = mycall))
}

screeplot.stepArchetypoids <- function(as,nar=NULL){
xx=c()
l=length(as)

for (i in 1:l){
xx[i]=as[[i]][[2]]
}
if (is.null(nar)){
plot(1:l,xx,type='b',xlab='Archetypoids', ylab='RSS')}
else{
plot(nar,xx,type='b',xlab='Archetypoids', ylab='RSS')}
}


RstepFArchetypoidswB <- function (data, k, norep=3,PM,PO,CT){
  mycall <- match.call()
  as <- list()
  for (i in 1:length(k)) {
    as[[i]] <- list()
    class(as[[i]]) <- "repArchetypoids"
    lass <- RFstepArchetypesModB(data,k[i],verbose=FALSE,nrep=norep,PM=PM,saveHistory=F,PO,CT) 
          as[[i]]<- RFstepArchetypoidsB(k[i],FALSE,data,lass,FALSE,k[i]-1,PM,PO,CT)
     
      
    
  }
  return(structure(as, class = "stepArchetypoids", call = mycall))
}


RstepFArchetypoidsbB <- function (data, k, norep=3,PM,PO,CT){
  mycall <- match.call()
  as <- list()
  for (i in 1:length(k)) {
    as[[i]] <- list()
    class(as[[i]]) <- "repArchetypoids"
    lass <- RFstepArchetypesModB(data,k[i],verbose=FALSE,nrep=norep,PM=PM,saveHistory=F,PO,CT) 
    ai <- bestModel(lass[[1]]) #archetypes
    ini_arch=c()
 for (j in 1:k[i]){
 ini_arch[j]=which.max(ai$betas[j,])
 }
 
as[[i]] <- RFrealArchetypesB(data,huge = 200, k[i], ini_arch,PM,PO,CT) #which beta archetypoid
      
  }
  return(structure(as, class = "stepArchetypoids", call = mycall))
}





stepRFArchetypoids3B <- function (data, k, norep=3,PM,PO=0.5,CT=6){
  mycall <- match.call()
  as <- list()
   asn=list()
   asw=list()
   asb=list()
  for (i in 1:length(k)) {
    asb[[i]] <- list()
    asn[[i]]=list()
    asw[[i]]=list()
    as[[i]]=list()
    class(asn[[i]]) <- "repArchetypoids"
    lass <- RFstepArchetypesModB(data,k[i],verbose=FALSE,nrep=norep,saveHistory=F,PM=PM,PO,CT) 
    ai <- bestModel(lass[[1]]) #archetypes
    ini_arch=c()
 for (j in 1:k[i]){
 ini_arch[j]=which.max(ai$betas[j,])
 }
 
Xs=data

  x11 <- t(Xs[,,1])
    x12 <- t(Xs[,,2])
    x1=rbind(x11,x12)
    data1=t(x1)

#from cand_beta
asb[[i]] <- RFrealArchetypesB(data1,huge = 200, k[i], ini_arch,PM=PM,PO,CT) 

#from nearest
asn[[i]]<- RFstepArchetypoidsB(k[i],TRUE,data,lass,FALSE,k[i]-1,PM=PM,PO,CT)
#from cand_alpha
asw[[i]]<- RFstepArchetypoidsB(k[i],FALSE,data,lass,FALSE,k[i]-1,PM=PM,PO,CT)

aux=c(asn[[i]][[2]],asw[[i]][[2]],asb[[i]][[2]])
wmi=which.min(aux)
auxa=list()
auxa[[1]]=asn[[i]]
auxa[[2]]=asw[[i]]
auxa[[3]]=asb[[i]]
as[[i]]= auxa[[wmi]]
      
  }
  return(structure(as, class = "stepArchetypoids", call = mycall))
}

