R/cecm.R

In evclust: Evidential Clustering

#' Constrained Evidential c-means algorithm
#'
#'\code{cecm} computes a credal partition from a matrix of attribute data and
#'pairwise constraints using the Constrained Evidential c-means (CECM) algorithm.
#'
#'CECM is a version of ECM allowing the user to specify pairwise constraints to guide
#'the clustering process. Pairwise constraints are of two kinds: must-link contraints are
#'pairs of objects that are known to belong to the same class, and cannot-link constraints
#'are pairs of objects that are known to belong to different classes. CECM can also learn
#'a metric for each cluster, like the Gustafson-Kessel algorithm in fuzzy clustering.
#'At each iteration, the algorithm solves a quadratic programming problem using an
#'interior ellipsoidal trust region and barrier function algorithm with dual solution
#'updating technique in the standard QP form (Ye, 1992).
#'
#'
#'If  initial prototypes \code{g0} are provided, the number of trials is automatically set to 1.
#'
#'Remark: Due to the use of the Matrix package, messages may be generated by R's (S4) method
#'dispatch mechanism. They are not error messages, and they can be ignored.
#'
#'
#' @param x input matrix of size n x d, where n is the number of objects and d the number of
#' attributes.
#' @param c Number of  clusters.
#' @param g0 Initial prototypes, matrix of size c x d. If not supplied, the prototypes are
#' initialized randomly.
#' @param type Type of focal sets ("simple": empty set, singletons and Omega;
#' "full": all \eqn{2^c} subsets of \eqn{\Omega}; "pairs": \eqn{\emptyset}, singletons,
#' \eqn{\Omega}, and all or selected pairs).
#' @param pairs Set of pairs to be included in the focal sets; if NULL, all pairs are
#' included. Used only if type="pairs".
#' @param ntrials Number of runs of the optimization algorithm (set to 1 if \code{g0} is  supplied).
#' @param ML Matrix nbML x 2 of must-link constraints. Each row of ML contains the indices
#' of objects that belong to the same class.
#' @param CL Matrix nbCL x 2 of cannot-link constraints. Each row of CL contains the indices
#' of objects that belong to different classes.
#' @param alpha Exponent of the cardinality in the cost function.
#' @param delta Distance to the empty set.
#' @param xi Tradeoff between the objective function Jecm and the constraints:
#' Jcecm=(1-xi)Jecm + xi Jconst.
#' @param distance Type of distance use: 0=Euclidean, 1=Mahalanobis.
#' @param epsi Minimum amount of improvement.
#' @param disp If TRUE (default), intermediate results are displayed.
#'
#' @return The credal partition (an object of class \code{"credpart"}).
#'
#'
#'@references V. Antoine, B. Quost, M.-H. Masson and T. Denoeux. CECM: Constrained
#'Evidential C-Means algorithm. Computational Statistics and Data Analysis, Vol. 56,
#'Issue 4, pages 894--914, 2012. 
#'
#'Y. Ye. On affine-scaling algorithm for nonconvex quadratic programming.
#'Math. Programming 56 (1992) 285--300.
#'
#'@author Thierry Denoeux (from a MATLAB code written by Violaine Antoine).
#'
#' @export
#' @import Matrix
#' @importFrom stats rnorm
#'
#' @seealso \code{\link{create_MLCL}}, \code{\link{makeF}}, \code{\link{extractMass}},
#' \code{\link{ecm}}, \code{\link{recm}}
#'
#' @examples ## Generation of a two-class dataset
#' \dontrun{
#' n<-30
#' x<-cbind(0.2*rnorm(n),rnorm(n))
#' y<-c(rep(1,n/2),rep(2,n/2))
#' x[(n/2+1):n,1]<-x[(n/2+1):n,1]+1
#' plot(x[,1],x[,2],asp=1,pch=y,col=y)
#' ## Generation of 10 constraints
#' const<-create_MLCL(y,nbConst=10)
#' ## Call of cecm
#' clus<-cecm(x=x,c=2,ML=const$M,CL=const$CL,delta=10)
#' plot(x[,1],x[,2],asp=1,pch=clus$y.pl,col=y)
#' }
cecm <- function(x,c,type='full',pairs=NULL,ntrials=1,ML,CL,g0=NULL,alpha=1,
                 delta=10,xi=0.5,distance=0,epsi=1e-3,disp=TRUE){
  bal<-xi
  x<-as.matrix(x)
  n<-nrow(x)
  nbAtt<-ncol(x)
  beta<-2;
  K<-c
  rho2<-delta^2
  F<-makeF(K,type,pairs)
  nbFoc<-nrow(F)
  card<- rowSums(F[2:nbFoc,])

  if((ntrials>1) & !is.null(g0)){
    print('WARNING: ntrials>1 and g0 provided. Parameter ntrials set to 1.')
    ntrials<-1
  }


  #--------------- constraint matrix reformulation --------
  nbML=nrow(ML)
  contraintesML<-sparseMatrix(i=ML[,1],j=ML[,2],x=rep(1,nbML),dims=c(n,n))
  diag(contraintesML)<-0
  contraintesML<-sign(contraintesML+t(contraintesML))
  contraintesML<-contraintesML*lower.tri(contraintesML)
  nbML<-length(which(contraintesML==1))

  nbCL=nrow(CL)
  contraintesCL<-sparseMatrix(i=CL[,1],j=CL[,2],x=rep(1,nbCL),dims=c(n,n))
  diag(contraintesCL)<-0
  contraintesCL<-sign(contraintesCL+t(contraintesCL))
  contraintesCL<-contraintesCL*lower.tri(contraintesCL)
  nbCL<-length(which(contraintesCL==1))

  # -- Setting q vector

  nbContParObjet <- rowSums(contraintesML+t(contraintesML))
  q<-kronecker(nbContParObjet,Matrix(c(1,rep(0,nbFoc-1)),ncol=1)) # Msg 1
  q<-as.vector(q)

  # -- Setting constraints matrix

  if(nbML==0) nbML<-1
  if(nbCL==0) nbCL<-1

  MLMat<-as.numeric(rowSums(F)==1) # getting singletons
  MLMat<- (MLMat %*% t(MLMat)) * (F %*% t(F))
  CLMat<-sign(F %*% t(F))
  MLMat<-MLMat*-sign(bal)/(2*nbML)
  CLMat<-CLMat*sign(bal)/(2*nbCL)

  # Constraints matrix with the respect of the constraints given in parameters

  MLaux<-kronecker(contraintesML,Matrix(1,nbFoc,nbFoc)) # Msg 2-3
  CLaux<-kronecker(contraintesCL,Matrix(1,nbFoc,nbFoc))

  contraintesMat<-kronecker(Matrix(1,n,n),MLMat)*MLaux+kronecker(Matrix(1,n,n),CLMat)*CLaux
  contraintesMat<-contraintesMat+t(contraintesMat)


  Jbest<-Inf
  for(itrial in 1:ntrials){
    #---------------------- initializations --------------------------------------
    # -- centroids initialization
    if(is.null(g0)){
      g <- x[sample(1:n,K),]+0.1*rnorm(K*nbAtt,K,nbAtt)
      #      g<-as.matrix(read.table('g.txt'))
    } else g<-g0

    # centers calculation for all the subsets
    gplus<-matrix(0,nbFoc-1,nbAtt)
    for(i in 2:nbFoc){
      fi <- F[i,]
      truc <- matrix(fi,K,nbAtt)
      gplus[i-1,] <- colSums(g*truc)/sum(fi)
    }

    # compute the Euclidean distance
    D<-matrix(0,n,nbFoc-1)
    for(j in 1:nbFoc-1){
      D[,j]<- rowSums((x-matrix(gplus[j,],n,nbAtt,byrow = TRUE))^2)
    }

    # compute masses (ECM without constraints)
    # Calculation of masses
    m <- matrix(0,n,nbFoc-1)
    for(i in 1:n){
      vect0 <- D[i,]
      for(j in 1:(nbFoc-1)){
        vect1 <- (rep(D[i,j],nbFoc-1)/vect0) ^(1/(beta-1))
        vect2 <-  rep(card[j]^(alpha/(beta-1)),nbFoc-1) /(card^(alpha/(beta-1)))
        vect3 <- vect1 * vect2
        m[i,j]<- 1/(  sum(vect3) + (card[j]^alpha * D[i,j]/rho2)^(1/(beta-1))  )
      }
    }
    m1 <- cbind(1-rowSums(m),m)

    dis<-setDistances(x,F,g,m,alpha,distance)
    D<-dis$D
    S<-dis$Splot
    Smeans<-dis$Smean

    # -- Setting H matrix
    aux=D%*% cbind(matrix(0,nbFoc-1,1), diag(nbFoc-1))+ cbind(matrix(1,n,1)*rho2, matrix(0,n,nbFoc-1))
    vectDist<-as.vector(t(aux))

    Card<-c(1,card)
    Card<-as.vector(matrix(Card^alpha,1,n*nbFoc))

    H <- (1-bal)*Diagonal(x=vectDist*Card/(n*nbFoc)) + bal*contraintesMat

    #------------------------ iterations--------------------------------
    notfinished<-TRUE
    gold<-g
    iter<-0

    #    ui<-rbind(diag(n*(nbFoc-1)),kronecker(diag(n),matrix(rep(-1,nbFoc-1),nrow=1)))
    #    ci<-c(rep(0,n*(nbFoc-1)),rep(-1,n))

    #    Aeq<-kronecker(diag(n),rep(1,nbFoc))
    #   beq<-rep(1,n)
    Aeq<-kronecker(Diagonal(n),rep(1,nbFoc))
    Aeq<-t(Aeq)
    beq<-rep(1,n)


    while(notfinished){
      iter<-iter+1
      mvec0<-as.vector(t(m1)) # masses used as initialization
      #    opt<-constrOptim(mvec0,f_cecm,grad_cecm,ui=ui,ci=ci,H=H,h=q,n=n,outer.eps=1e-7,outer.iterations = 1000)
      #    masses<-opt$par

      qp<-solqp(Q=H,A=Aeq,b=beq,c=q,x=mvec0)
      masses<-qp$x
      m1<-matrix(masses,n,nbFoc,byrow=TRUE)
      m<-m1[,2:nbFoc]

      # Calculation of centers
      g<-setCentersECM(x,m,F,Smeans,alpha,beta)
      dist<-setDistances(x,F,g,m,alpha,distance)
      D<-dist$D
      Smeans<-dist$Smean

      # H matrix
      aux=D%*% cbind(matrix(0,nbFoc-1,1), diag(nbFoc-1))+ cbind(matrix(1,n,1)*rho2, matrix(0,n,nbFoc-1))
      vectDist<-as.vector(t(aux))
      H<-(1-bal)*Diagonal(x=vectDist*Card/(n*nbFoc))+bal*contraintesMat

      J <- as.numeric(t(masses)%*% H %*% masses + bal)
      delta<-max(abs(g-gold))
      if (disp) print(c(iter,J,delta))

      notfinished <- (delta>epsi)
      gold <- g
    } # end while loop


    if(J<Jbest){
      Jbest<-J
      mbest<-m1
      gbest<-g
      Smeansbest<-Smeans
    }
    res<-c(itrial,J,Jbest)
    names(res)<-NULL
    if (disp) print(res)
  } #end for loop iter
  #--------------------------- end of iterations ----------------------------
  clus<-extractMass(mbest,F,g=gbest,S=Smeansbest,method="cecm",crit=Jbest,
                    param=list(alpha=alpha,beta=beta,delta=delta))
  return(clus)

}