--- a
+++ b/R/LRAcluster.R
@@ -0,0 +1,628 @@
+#' @name LRAcluster
+#' @aliases LRAcluster
+#' @title integrated analysis of cancer omics data by low-rank approximation
+#' @description The LRAcluster function is the main interface of this package, it gets a list of matrix as input and outputs the coordinates of the samples in the reduced space and the explained potential.
+#' @param data a list of data matrix,please keep the columns are the same order of samples
+#' @param types a list of data types, can be binary, gaussian, or poisson
+#' @param dimension the reduced dimension
+#' @param names data names
+#' @return A list contains the following component:
+#'
+#'         \code{coordinate} A matrix of the coordinates of all the samples in the reduced space
+#'
+#'         \code{potential}  ratio of explained variance
+#' @keywords internal
+#' @author Dingming Wu, Dongfang Wang
+#' @references Wu D, Wang D, Zhang MQ, Gu J (2015). Fast dimension reduction and integrative clustering of multi-omics data using low-rank approximation: application to cancer molecular classification. BMC Genomics, 16(1):1022.
+LRAcluster <- function(data, types, dimension = 2, names = as.character(1:length(data)))
+{
+  #--------#
+  # binary #
+  #--------#
+  epsilon.binary<-2.0
+  check.binary.row<-function(arr)
+  {
+    if (sum(!is.na(arr))==0)
+    {
+      return (F)
+    }
+    else
+    {
+      idx<-!is.na(arr)
+      if (sum(arr[idx])==0 || sum(arr[idx])==sum(idx))
+      {
+        return (F)
+      }
+      else
+      {
+        return (T)
+      }
+    }
+  }
+  check.binary<-function(mat,name)
+  {
+    index<-apply(mat,1,check.binary.row)
+    n<-sum(!index)
+    if (n>0)
+    {
+      w<-paste("Warning: ",name," have ",as.character(n)," invalid lines!",sep="")
+      warning(w)
+    }
+    mat_c<-mat[index,]
+    rownames(mat_c)<-rownames(mat)[index]
+    colnames(mat_c)<-colnames(mat)
+    return (mat_c)
+  }
+
+  base.binary.row<-function(arr)
+  {
+    idx<-!is.na(arr)
+    n<-sum(idx)
+    m<-sum(arr[idx])
+    return (log(m/(n-m)))
+  }
+
+  base.binary<-function(mat)
+  {
+    mat_b<-matrix(0,nrow(mat),ncol(mat))
+    ar_b<-apply(mat,1,base.binary.row)
+    mat_b[1:nrow(mat_b),]<-ar_b
+    return (mat_b)
+  }
+
+  update.binary<-function(mat,mat_b,mat_now,eps)
+  {
+    mat_p<-mat_b+mat_now
+    mat_u<-matrix(0,nrow(mat),ncol(mat))
+    idx1<-!is.na(mat) & mat==1
+    idx0<-!is.na(mat) & mat==0
+    index<-is.na(mat)
+    arr<-exp(mat_p)
+    mat_u[index]<-mat_now[index]
+    mat_u[idx1]<-mat_now[idx1]+eps*epsilon.binary/(1.0+arr[idx1])
+    mat_u[idx0]<-mat_now[idx0]-eps*epsilon.binary*arr[idx0]/(1.0+arr[idx0])
+    return (mat_u)
+  }
+
+  stop.binary<-function(mat,mat_b,mat_now,mat_u)
+  {
+    index<-!is.na(mat)
+    mn<-mat_b+mat_now
+    mu<-mat_b+mat_u
+    arn<-exp(mn)
+    aru<-exp(mu)
+    idx1<-!is.na(mat) & mat==1
+    idx0<-!is.na(mat) & mat==0
+    lgn<-sum(log(arn[idx1]/(1+arn[idx1])))+sum(log(1/(1+arn[idx0])))
+    lgu<-sum(log(aru[idx1]/(1+aru[idx1])))+sum(log(1/(1+aru[idx0])))
+    return (lgu-lgn)
+  }
+
+  LL.binary<-function(mat,mat_b,mat_u)
+  {
+    index<-!is.na(mat)
+    mu<-mat_b+mat_u
+    aru<-exp(mu)
+    idx1<-!is.na(mat) & mat==1
+    idx0<-!is.na(mat) & mat==0
+    lgu<-sum(log(aru[idx1]/(1+aru[idx1])))+sum(log(1/(1+aru[idx0])))
+    return (lgu)
+  }
+
+  LLmax.binary<-function(mat)
+  {
+    return (0)
+  }
+
+  LLmin.binary<-function(mat,mat_b)
+  {
+    index<-!is.na(mat)
+    aru<-exp(mat_b)
+    idx1<-!is.na(mat) & mat==1
+    idx0<-!is.na(mat) & mat==0
+    lgu<-sum(log(aru[idx1]/(1+aru[idx1])))+sum(log(1/(1+aru[idx0])))
+    return (lgu)
+  }
+
+  binary_type_base <- function( data,dimension=2 ,name="test")
+  {
+    data<-check.binary(data,name)
+    data_b<-base.binary(data)
+    data_now<-matrix(0,nrow(data),ncol(data))
+    data_u<-update.binary(data,data_b,data_now)
+    data_u<-nuclear_approximation(data_u,dimension)
+    while (T)
+    {
+      thr<-stop.binary(data,data_b,data_now,data_u)
+      message(thr)
+      if (thr<0.2)
+      {
+        break
+      }
+      data_now<-data_u
+      data_u<-update.binary(data,data_b,data_now)
+      data_u<-nuclear_approximation(data_u,dimension)
+    }
+    return (data_now)
+  }
+
+  #----------#
+  # gaussian #
+  #----------#
+
+  epsilon.gaussian=0.5
+
+  check.gaussian.row<-function(arr)
+  {
+    if (sum(!is.na(arr))==0)
+    {
+      return (F)
+    }
+    else
+    {
+      return (T)
+    }
+  }
+  check.gaussian<-function(mat,name)
+  {
+    index<-array(T,nrow(mat))
+    for(i in 1:nrow(mat))
+    {
+      if (sum(is.na(mat[i,])==ncol(mat)))
+      {
+        war<-paste("Warning: ",name,"'s ",as.character(i)," line is all NA. Delete this line",sep="")
+        warning(war)
+        index[i]<-F
+      }
+    }
+    mat_c<-mat[index,]
+    rownames(mat_c)<-rownames(mat)[index]
+    colnames(mat_c)<-colnames(mat)
+    return (mat_c)
+  }
+
+  base.gaussian.row<-function(arr)
+  {
+    idx<-!is.na(arr)
+    return (mean(arr[idx]))
+  }
+
+  base.gaussian<-function(mat)
+  {
+    mat_b<-matrix(0,nrow(mat),ncol(mat))
+    ar_b<-apply(mat,1,base.gaussian.row)
+    mat_b[1:nrow(mat_b),]<-ar_b
+    return (mat_b)
+  }
+
+  update.gaussian<-function(mat,mat_b,mat_now,eps)
+  {
+    mat_p<-mat_b+mat_now
+    mat_u<-matrix(0,nrow(mat),ncol(mat))
+    index<-!is.na(mat)
+    mat_u[index]<-mat_now[index]+eps*epsilon.gaussian*(mat[index]-mat_p[index])
+    index<-is.na(mat)
+    mat_u[index]<-mat_now[index]
+    return (mat_u)
+  }
+
+  stop.gaussian<-function(mat,mat_b,mat_now,mat_u)
+  {
+    index<-!is.na(mat)
+    mn<-mat_b+mat_now
+    mu<-mat_b+mat_u
+    ren<-mat[index]-mn[index]
+    reu<-mat[index]-mu[index]
+    lgn<- -0.5*sum(ren*ren)
+    lgu<- -0.5*sum(reu*reu)
+    return (lgu-lgn)
+  }
+
+  LL.gaussian<-function(mat,mat_b,mat_u)
+  {
+    index<-!is.na(mat)
+    mu<-mat_b+mat_u
+    reu<-mat[index]-mu[index]
+    lgu<- -0.5*sum(reu*reu)
+    return (lgu)
+  }
+
+  LLmax.gaussian<-function(mat)
+  {
+    return (0.0)
+  }
+
+  LLmin.gaussian<-function(mat,mat_b)
+  {
+    index<-!is.na(mat)
+    reu<-mat[index]-mat_b[index]
+    lgu<- -0.5*sum(reu*reu)
+    return (lgu)
+  }
+
+  gaussian_base<-function(data,dimension=2,name="test")
+  {
+    data<-check.gaussian(data,name)
+    data_b<-base.gaussian(data)
+    data_now<-matrix(0,nrow(data),ncol(data))
+    data_u<-update.gaussian(data,data_b,data_now)
+    data_u<-nuclear_approximation(data_u,dimension)
+    while(T)
+    {
+      thr<-stop.gaussian(data,data_b,data_now,data_u)
+      message(thr)
+      if (thr<0.2)
+      {
+        break
+      }
+      data_now<-data_u
+      data_u<-update.gaussian(data,data_b,data_now)
+      data_u<-nuclear_approximation(data_u,dimension)
+    }
+    return (data_now)
+  }
+
+  #---------#
+  # poisson #
+  #---------#
+
+  epsilon.poisson<-0.5
+
+  check.poisson.row<-function(arr)
+  {
+    if (sum(!is.na(arr))==0)
+    {
+      return (F)
+    }
+    else
+    {
+      idx<-!is.na(arr)
+      if (sum(arr[idx]<0)>0)
+      {
+        return (F)
+      }
+      else
+      {
+        return (T)
+      }
+    }
+  }
+
+  check.poisson<-function(mat,name)
+  {
+    w<-paste(name," is poisson type. Add 1 to all counts",sep="")
+    message(w)
+    index<-apply(mat,1,check.poisson.row)
+    n<-sum(!index)
+    if (n>0)
+    {
+      w<-paste("Warning: ",name," have ",as.character(n)," invalid lines!",sep="")
+      warning(w)
+    }
+    mat_c<-mat[index,]+1
+    rownames(mat_c)<-rownames(mat)[index]
+    colnames(mat_c)<-colnames(mat)
+    return (mat_c)
+  }
+
+  base.poisson.row<-function(arr)
+  {
+    idx<-!is.na(arr)
+    m<-sum(log(arr[idx]))
+    n<-sum(idx)
+    return(m/n)
+  }
+
+  base.poisson<-function(mat)
+  {
+    mat_b<-matrix(0,nrow(mat),ncol(mat))
+    ar_b<-apply(mat,1,base.poisson.row)
+    mat_b[1:nrow(mat_b),]<-ar_b
+    return (mat_b)
+  }
+
+  update.poisson<-function(mat,mat_b,mat_now,eps)
+  {
+    mat_p<-mat_b+mat_now
+    mat_u<-matrix(0,nrow(mat),ncol(mat))
+    index<-!is.na(mat)
+    mat_u[index]<-mat_now[index]+eps*epsilon.poisson*(log(mat[index])-mat_p[index])
+    index<-is.na(mat)
+    mat_u[index]<-mat_now[index]
+    return (mat_u)
+  }
+
+  stop.poisson<-function(mat,mat_b,mat_now,mat_u)
+  {
+    index<-!is.na(mat)
+    mn<-mat_b+mat_now
+    mu<-mat_b+mat_u
+    lgn<-sum(mat[index]*mn[index]-exp(mn[index]))
+    lgu<-sum(mat[index]*mu[index]-exp(mu[index]))
+    return (lgu-lgn)
+  }
+
+  LL.poisson<-function(mat,mat_b,mat_u)
+  {
+    index<-!is.na(mat)
+    mu<-mat_b+mat_u
+    lgu<-sum(mat[index]*mu[index]-exp(mu[index]))
+    return (lgu)
+  }
+
+  LLmax.poisson<-function(mat)
+  {
+    index<-!is.na(mat)
+    lgu<-sum(mat[index]*log(mat[index])-mat[index])
+    return (lgu)
+  }
+
+  LLmin.poisson<-function(mat,mat_b)
+  {
+    index<-!is.na(mat)
+    lgu<-sum(mat[index]*mat_b[index]-exp(mat_b[index]))
+    return (lgu)
+  }
+
+  poisson_type_base<-function(data,dimension=2,name="test")
+  {
+    data<-check.poisson(data,name)
+    data_b<-base.poisson(data)
+    data_now<-matrix(0,nrow(data),ncol(data))
+    data_u<-update.poisson(data,data_b,data_now)
+    data_u<-nuclear_approximation(data_u,dimension)
+    while(T)
+    {
+      thr<-stop.poisson(data,data_b,data_now,data_u)
+      message(thr)
+      if (thr<0.2)
+      {
+        break
+      }
+      data_now<-data_u
+      data_u<-update.poisson(data,data_b,data_now)
+      data_u<-nuclear_approximation(data_u,dimension)
+    }
+    return (data_now)
+  }
+
+  #----#
+  # na #
+  #----#
+
+  nuclear_approximation<-function(mat,dimension)
+  {
+    svd<-svd(mat,nu=0,nv=0)
+    if (dimension<length(svd$d))
+    {
+      lambda<-svd$d[dimension+1]
+      svd<-svd(mat,nu=dimension,nv=dimension)
+      indexh<-svd$d>lambda
+      indexm<-svd$d<lambda
+      dia<-array(svd$d,length(svd$d))
+      dia[indexh]<-dia[indexh]-lambda
+      dia[indexm]<-0
+      mat_low<-svd$u%*%diag(c(dia[1:dimension],0))[1:dimension,1:dimension]%*%t(svd$v)
+    }
+    else
+    {
+      mat_low<-mat
+    }
+    return (mat_low)
+  }
+
+  #------------#
+  # LRAcluster #
+  #------------#
+  check.matrix.element<-function(x)
+  {
+    if (!is.matrix(x))
+    {
+      return (T)
+    }
+    else
+    {
+      return (F)
+    }
+  }
+
+  ncol.element<-function(x)
+  {
+    return (ncol(x))
+  }
+
+  nrow.element<-function(x)
+  {
+    return (nrow(x))
+  }
+
+  check<-function(mat,type,name)
+  {
+    if (type=="binary")
+    {
+      return (check.binary(mat,name))
+    }
+    else if (type=="gaussian")
+    {
+      return (check.gaussian(mat,name))
+    }
+    else if (type=="poisson")
+    {
+      return (check.poisson(mat,name))
+    }
+    else
+    {
+      e<-paste("unknown type ",type,sep="")
+      stop(e)
+    }
+  }
+
+  eps<-0.0
+  if (!is.list(data))
+  {
+    stop("the input data must be a list!")
+  }
+  c<-sapply(data,check.matrix.element)
+  if (sum(c)>0)
+  {
+    stop("each element of input list must be a matrix!")
+  }
+  c<-sapply(data,ncol.element)
+  if (length(levels(factor(c)))>1)
+  {
+    stop("each element of input list must have the same column number!")
+  }
+  if (length(data)!=length(types))
+  {
+    stop("data and types must be the same length!")
+  }
+  nSample<-c[1]
+  loglmin<-0
+  loglmax<-0
+  loglu<-0.0
+  nData<-length(data)
+  for (i in 1:nData)
+  {
+    data[[i]]<-check(data[[i]],types[[i]],names[[i]])
+  }
+  nGeneArr<-sapply(data,nrow.element)
+  nGene<-sum(nGeneArr)
+  indexData<-list()
+  k=1
+  for(i in 1:nData)
+  {
+    indexData[[i]]<- (k):(k+nGeneArr[i]-1)
+    k<-k+nGeneArr[i]
+  }
+  base<-matrix(0,nGene,nSample)
+  now<-matrix(0,nGene,nSample)
+  update<-matrix(0,nGene,nSample)
+  thr<-array(0,nData)
+  for (i in 1:nData)
+  {
+    if (types[[i]]=="binary")
+    {
+      base[indexData[[i]],]<-base.binary(data[[i]])
+      loglmin<-loglmin+LLmin.binary(data[[i]],base[indexData[[i]],])
+      loglmax<-loglmax+LLmax.binary(data[[i]])
+    }
+    else if (types[[i]]=="gaussian")
+    {
+      base[indexData[[i]],]<-base.gaussian(data[[i]])
+      loglmin<-loglmin+LLmin.gaussian(data[[i]],base[indexData[[i]],])
+      loglmax<-loglmax+LLmax.gaussian(data[[i]])
+    }
+    else if (types[[i]]=="poisson")
+    {
+      base[indexData[[i]],]<-base.poisson(data[[i]])
+      loglmin<-loglmin+LLmin.poisson(data[[i]],base[indexData[[i]],])
+      loglmax<-loglmax+LLmax.poisson(data[[i]])
+    }
+  }
+  for (i in 1:nData)
+  {
+    if (types[[i]]=="binary")
+    {
+      update[indexData[[i]],]<-update.binary(data[[i]],base[indexData[[i]],],now[indexData[[i]],],exp(eps))
+    }
+    else if (types[[i]]=="gaussian")
+    {
+      update[indexData[[i]],]<-update.gaussian(data[[i]],base[indexData[[i]],],now[indexData[[i]],],exp(eps))
+    }
+    else if (types[[i]]=="poisson")
+    {
+      update[indexData[[i]],]<-update.poisson(data[[i]],base[indexData[[i]],],now[indexData[[i]],],exp(eps))
+    }
+  }
+  update<-nuclear_approximation(update,dimension)
+  nIter<-0
+  thres<-array(Inf,3)
+  epsN<-array(Inf,2)
+  while(T)
+  {
+    for (i in 1:nData)
+    {
+      if (types[[i]]=="binary")
+      {
+        thr[i]<-stop.binary(data[[i]],base[indexData[[i]],],now[indexData[[i]],],update[indexData[[i]],])
+      }
+      else if (types[[i]]=="gaussian")
+      {
+        thr[i]<-stop.gaussian(data[[i]],base[indexData[[i]],],now[indexData[[i]],],update[indexData[[i]],])
+      }
+      else if (types[[i]]=="poisson")
+      {
+        thr[i]<-stop.poisson(data[[i]],base[indexData[[i]],],now[indexData[[i]],],update[indexData[[i]],])
+      }
+    }
+    nIter<-nIter+1
+    thres[1]<-thres[2]
+    thres[2]<-thres[3]
+    thres[3]<-sum(thr)
+    epsN[1]<-epsN[2]
+    epsN[2]<-eps
+    if (nIter>5)
+    {
+      if (runif(1)<thres[1]*thres[3]/(thres[2]*thres[2]+thres[1]*thres[3]))
+      {
+        eps<-epsN[1]+0.05*runif(1)-0.025
+      }
+      else
+      {
+        eps<-epsN[2]+0.05*runif(1)-0.025
+      }
+      if (eps< -0.7)
+      {
+        eps<- 0
+        epsN<-c(0,0)
+      }
+      if (eps > 1.4)
+      {
+        eps<-0
+        epsN<-c(0,0)
+      }
+    }
+    if (sum(thr)<nData*0.2)
+    {
+      break
+    }
+    now<-update
+    for (i in 1:nData)
+    {
+      if (types[[i]]=="binary")
+      {
+        update[indexData[[i]],]<-update.binary(data[[i]],base[indexData[[i]],],now[indexData[[i]],],exp(eps))
+      }
+      else if (types[[i]]=="gaussian")
+      {
+        update[indexData[[i]],]<-update.gaussian(data[[i]],base[indexData[[i]],],now[indexData[[i]],],exp(eps))
+      }
+      else if (types[[i]]=="poisson")
+      {
+        update[indexData[[i]],]<-update.poisson(data[[i]],base[indexData[[i]],],now[indexData[[i]],],exp(eps))
+      }
+    }
+    update<-nuclear_approximation(update,dimension)
+  }
+  for (i in 1:nData)
+  {
+    if (types[[i]]=="binary")
+    {
+      loglu<-loglu+LL.binary(data[[i]],base[indexData[[i]],],update[indexData[[i]],])
+    }
+    else if (types[[i]]=="gaussian")
+    {
+      loglu<-loglu+LL.gaussian(data[[i]],base[indexData[[i]],],update[indexData[[i]],])
+    }
+    else if (types[[i]]=="poisson")
+    {
+      loglu<-loglu+LL.poisson(data[[i]],base[indexData[[i]],],update[indexData[[i]],])
+    }
+  }
+  sv<-svd(update,nu=0,nv=dimension)
+  coordinate<-diag(c(sv$d[1:dimension],0))[1:dimension,1:dimension]%*%t(sv$v)
+  colnames(coordinate)<-colnames(data[[1]])
+  rownames(coordinate)<-paste("PC ",as.character(1:dimension),sep="")
+  ratio<-(loglu-loglmin)/(loglmax-loglmin)
+  return (list("coordinate"=coordinate,"potential"=ratio))
+}