R/ICMg.links.sampler.SeedsAndWeights.R

In ChristophRau/wMICA: weighted Maximal Information Component Analysis

`ICMg.links.sampler.SeedsAndWeights` <-
function(L, C,Seeds=c(), Weight=0, alpha=10, beta=0.01,  B.num=8, B.size=100, S.num=20, S.size=10, C.boost=1) {

  cat("Sampling ICM...\n")
  ## L is Nx3 matrix of link endpoints, with the third column being weights.  Will auto-fix if you have only two rows
  ## C is the number of latent component
  ## Seeds is a Kx2 matrix of nodes and modules.
  ## alpha and beta are the hyperparameters
  ## Weight is the weight (in extra rounds) to give the seeds

  ## B.size is the size of one burnin round, after which convergence and component counts are reported
  ## B.num is the number of these burnin rounds
  ## S size is likewise the size of one sample round, after which samples are taken
  ## S.num is the number of these sample rounds

  Nlinks <- dim(L)[1]
  Nnodes <- max(L)

  cat("nodes:",Nnodes,"links:",Nlinks,"components:",C,"alpha:",alpha,"beta:",beta,"\n")

  out <- list()
  out$z <- vector("double", Nlinks)     # link l is generated by z[l]
  out$n <- vector("double", C) # n[z] links generated by component z
  out$q <- matrix(0, C, Nnodes) # k[z,i] links generated by z end at i

  ## Start from random initialization
  out$z <- sample(1:C, Nlinks, replace=TRUE)
  out$n <- as.vector(table(factor(out$z, levels=1:C)))

  if(ncol(L)==2){
    L=cbind(L,1)
  }

  for (l in 1:Nlinks) {
    out$q[out$z[l], L[l,1]] <- out$q[out$z[l], L[l,1]]+L[l,3]
    out$q[out$z[l], L[l,2]] <- out$q[out$z[l], L[l,2]]+L[l,3]
  }

  for (i in nrow(Seeds)){
    out$q[Seeds[i,2],Seeds[i,1]] <- Weight
  }

  ## Initialize random number generator
  ICMg.randominit()

  ## Initialize some necessary data structures
  samples <- list(z=matrix(0, S.num, Nlinks), conv=vector("numeric", B.num+S.num),
                  counts=matrix(0,B.num+S.num,C))
  conv <- vector("numeric", Nlinks)

  cat("Sampling", B.num*B.size + S.num*S.size,"iterations\n")

  ## Do burnin and monitor convergence
  cat("Burnin iterations:",B.num*B.size,"\n")
  cat("I: 0 n(z):",out$n,"\n")
  for (b in 1:B.num) {

    Lindices <- sample(Nlinks)
    out <- ICMg.links.wrapper.Weights(L, B.size, Nlinks, Nnodes, Lindices, C,
                                      out$z, out$q, out$n, alpha, beta, conv, C.boost)
    samples$conv[b] <- mean(log(out$conv))
    samples$counts[b,] <- out$n
    cat("I:",b*B.size,"conv:",samples$conv[b],"n(z):",samples$counts[b,],"\n")
    #print(out$q)
  }

  ## After burnin take samples
  cat("Sample iterations:",S.num*S.size,"\n")
  for (s in 1:S.num) {

    Lindices <- sample(Nlinks)
    out <- ICMg.links.wrapper.Weights(L, S.size, Nlinks, Nnodes, Lindices, C,
                                      out$z, out$q, out$n, alpha, beta, conv, C.boost)
    samples$conv[B.num +s] <- mean(log(out$conv))
    samples$counts[B.num+s,] <- out$n
    cat("I:",B.num*B.size + s*S.size,"conv:",samples$conv[B.num +s],
        "n(z):",samples$counts[B.num+s,],"\n")
    samples$z[s,] <- out$z
    #print(out$q)
  }

  ## Save parameters
  samples$B.num <- B.num
  samples$B.size <- B.size
  samples$S.num <- S.num
  samples$S.size <- S.size
  samples$alpha <- alpha
  samples$beta <- beta

  cat("Sampling finished\n")
  return(samples)
}