R/ctmcProbabilistic.R

In spedygiorgio/markovchain: Easy Handling Discrete Time Markov Chains

#' @title rctmc
#'
#' @description The function generates random CTMC transitions as per the
#'   provided generator matrix.
#' @usage rctmc(n, ctmc, initDist = numeric(), T = 0, include.T0 = TRUE,
#'   out.type = "list")
#'
#' @param n The number of samples to generate.
#' @param ctmc The CTMC S4 object.
#' @param initDist The initial distribution of states.
#' @param T The time up to which the simulation runs (all transitions after time
#'   T are not returned).
#' @param include.T0 Flag to determine if start state is to be included.
#' @param out.type "list" or "df"
#'
#' @details In order to use the T0 argument, set n to Inf.
#' @return Based on out.type, a list or a data frame is returned. The returned
#'   list has two elements - a character vector (states) and a numeric vector
#'   (indicating time of transitions). The data frame is similarly structured.
#' @references 
#' Introduction to Stochastic Processes with Applications in the Biosciences 
#' (2013), David F. Anderson, University of Wisconsin at Madison
#' 
#' @author Sai Bhargav Yalamanchi
#' @seealso \code{\link{generatorToTransitionMatrix}},\code{\link{ctmc-class}}
#' @examples
#' energyStates <- c("sigma", "sigma_star")
#' byRow <- TRUE
#' gen <- matrix(data = c(-3, 3, 1, -1), nrow = 2,
#'              byrow = byRow, dimnames = list(energyStates, energyStates))
#' molecularCTMC <- new("ctmc", states = energyStates, 
#'                      byrow = byRow, generator = gen, 
#'                      name = "Molecular Transition Model")   
#'
#' statesDist <- c(0.8, 0.2)
#' rctmc(n = Inf, ctmc = molecularCTMC, T = 1)
#' rctmc(n = 5, ctmc = molecularCTMC, initDist = statesDist, include.T0 = FALSE)
#' @export
rctmc <- function(n, ctmc, initDist = numeric(), T = 0, include.T0 = TRUE, 
                  out.type = "list") {
  if (identical(initDist, numeric()))
    state <- sample(ctmc@states, 1) # sample state randomly
  else if (length(initDist) != dim(ctmc) | round(sum(initDist), 5) != 1)
    stop("Error! Provide a valid initial state probability distribution")
  else 
    state <- sample(ctmc@states, 1, prob = initDist) # if valid probability distribution,
  # sample accordingly
  
  # obtain transition probability matrix from the generator matrix
  trans <- generatorToTransitionMatrix(ctmc@generator)
  
  
  states <- c()
  time <- c()
  if (include.T0 == TRUE){
    states <- c(states, state)
    time <- c(time, 0)
  }
  
  t <- 0
  i <- 1
  while (i <= n){
    idx <- which(ctmc@states == state)
    t <- t + rexp(1, -ctmc@generator[idx, idx])
    state <- ctmc@states[sample(1:dim(ctmc), 1, prob = trans[idx, ])]
    
    if(T > 0 & t > T)
      break
    
    states <- c(states, state)
    time <- c(time, t)
    i <- i + 1
  }
  
  out <- list(states, time)
  if (out.type == "list")
    return(out)
  else if(out.type == "df"){
    df <- data.frame(matrix(unlist(out), nrow = length(states)))
    names(df) <- c("states", "time")
    return(df)
  }
  else
    stop("Not a valid output type")
}

#' @title Return the generator matrix for a corresponding transition matrix
#' 
#' @description Calculate the generator matrix for a 
#'              corresponding transition matrix
#' 
#' @param P transition matrix between time 0 and t
#' @param t time of observation
#' @param method "logarithm" returns the Matrix logarithm of the transition matrix
#'
#' @return A matrix that represent the generator of P
#' @export
#'
#' @examples
#' mymatr <- matrix(c(.4, .6, .1, .9), nrow = 2, byrow = TRUE)
#' Q <- transition2Generator(P = mymatr)
#' expm::expm(Q)
#'  
#' @seealso \code{\link{rctmc}}
transition2Generator<-function(P, t = 1,method = "logarithm") {
  if (method == "logarithm") {
    Q = logm(P)/t
  } #else 
  return(Q)
}


#' Returns a generator matrix corresponding to frequency matrix
#' 
#' @description The function provides interface to calculate generator matrix corresponding to 
#' a frequency matrix and time taken
#' 
#' @param P relative frequency matrix
#' @param t (default value = 1)
#' @param method one among "QO"(Quasi optimaisation), "WA"(weighted adjustment), "DA"(diagonal adjustment)
#' @param logmethod method for computation of matrx algorithm (by default : Eigen)
#' 
#' @return returns a generator matix with same dimnames
#' 
#' @references E. Kreinin and M. Sidelnikova: Regularization Algorithms for
#' Transition Matrices. Algo Research Quarterly 4(1):23-40, 2001
#' 
#' @examples
#' sample <- matrix(c(150,2,1,1,1,200,2,1,2,1,175,1,1,1,1,150),nrow = 4,byrow = TRUE)
#' sample_rel = rbind((sample/rowSums(sample))[1:dim(sample)[1]-1,],c(rep(0,dim(sample)[1]-1),1)) 
#' freq2Generator(sample_rel,1)
#' 
#' data(tm_abs)
#' tm_rel=rbind((tm_abs/rowSums(tm_abs))[1:7,],c(rep(0,7),1))
#' ## Derive quasi optimization generator matrix estimate
#' freq2Generator(tm_rel,1)
#' 
#' @export
#' 
freq2Generator <- function(P,t = 1,method = "QO",logmethod = "Eigen"){
  if(requireNamespace('ctmcd', quietly = TRUE)) {
  if(method == "QO"){
    out <- ctmcd::gmQO(P, t, logmethod)
  } else if(method == "WA") {
    out <- ctmcd::gmWA(P, t, logmethod)
  } else if(method == "DA") {
    out <- ctmcd::gmDA(P, t, logmethod)
  }
  } else {
    warning("package ctmcd is not installed")
    out = NULL
  }
  return(out)
}

#' @title Returns expected hitting time from state i to state j
#' 
#' @description Returns expected hitting time from state i to state j
#' 
#' @usage ExpectedTime(C,i,j,useRCpp)
#' 
#' @param C A CTMC S4 object
#' @param i Initial state i
#' @param j Final state j
#' @param useRCpp logical whether to use Rcpp
#' 
#' @details According to the theorem, holding times for all states except j should be greater than 0.
#' 
#' @return A numerical value that returns expected hitting times from i to j
#' 
#' @references Markovchains, J. R. Norris, Cambridge University Press
#' 
#' @author Vandit Jain
#' 
#' @examples
#' states <- c("a","b","c","d")
#' byRow <- TRUE
#' gen <- matrix(data = c(-1, 1/2, 1/2, 0, 1/4, -1/2, 0, 1/4, 1/6, 0, -1/3, 1/6, 0, 0, 0, 0),
#' nrow = 4,byrow = byRow, dimnames = list(states,states))
#' ctmc <- new("ctmc",states = states, byrow = byRow, generator = gen, name = "testctmc")
#' ExpectedTime(ctmc,1,4,TRUE)
#' 
#' @export
ExpectedTime <- function(C,i,j,useRCpp = TRUE){
  # take generator from ctmc-class object
  Q <- C@generator
  
  # in case where generator is written column wise
  if(C@byrow==FALSE){
    Q <- t(Q)
  }
  NoofStates <- dim(C)
  
  Exceptj <- c(1:NoofStates)
  # create vector with all values from 1:NoofStates except j 
  Exceptj <- which(Exceptj!=j)
  
  # build matrix with vlaues from Q such that row!=j or column!=j
  Q_Exceptj <- Q[Exceptj,Exceptj]
  
  # check for positivity of holding times except for state j
  if(!all(diag(Q_Exceptj)!=0)){
    stop("Holding times for all states except j should be greater than 0")
  }
  
  # get b for solving the system of linear equation Ax = b where A is Q_Exceptj
  b <- rep(-1,dim(Q_Exceptj)[1])
  
  # use solve function from base packge to solve Ax = b
  if(useRCpp == TRUE){
    out <- .ExpectedTimeRCpp(Q_Exceptj,b)
  } else {
    out <- solve(Q_Exceptj,b)
  }
  
  # out will be of size NoofStates-1, hence the adjustment for different cases of i>=<j
  if(i<j){
    return(out[[i]])
  } else if(i == j) {
    return(0)
  } else {
    return(out[[i-1]])
  }
}

#' Calculating probability from a ctmc object
#' 
#' @description 
#' This function returns the probability of every state at time t under different conditions
#' 
#' @usage probabilityatT(C,t,x0,useRCpp)
#' 
#' @param C A CTMC S4 object
#' @param t final time t
#' @param x0 initial state
#' @param useRCpp logical whether to use RCpp implementation
#' 
#' @details The initial state is not mandatory, In case it is not provided, 
#' function returns a matrix of transition function at time \code{t} else it returns
#' vector of probaabilities of transition to different states if initial state was \code{x0}
#' 
#' @return returns a vector or a matrix in case \code{x0} is provided or not respectively.
#' 
#' @references INTRODUCTION TO STOCHASTIC PROCESSES WITH R, ROBERT P. DOBROW, Wiley
#' 
#' @author Vandit Jain
#' 
#' @examples
#' states <- c("a","b","c","d")
#' byRow <- TRUE
#' gen <- matrix(data = c(-1, 1/2, 1/2, 0, 1/4, -1/2, 0, 1/4, 1/6, 0, -1/3, 1/6, 0, 0, 0, 0),
#' nrow = 4,byrow = byRow, dimnames = list(states,states))
#' ctmc <- new("ctmc",states = states, byrow = byRow, generator = gen, name = "testctmc")
#' probabilityatT(ctmc,1,useRCpp = TRUE)
#' 
#' @export
probabilityatT <- function(C, t, x0, useRCpp = TRUE){
  
  if(!is(C,"ctmc")){
    stop("Provided object is not a ctmc object")
  }
  if(t < 0){
    stop("Time provided should be greater than equal to 0")
  }
  # take generator from ctmc-class object
  Q <- C@generator
  
  # in case where generator is written column wise
  if(C@byrow==FALSE){
    Q <- t(Q)
  }
  NoofStates <- dim(C)
  
  
  # calculate transition functoin at time t using Kolmogorov backward equation
  if(useRCpp == TRUE){
    P <- .probabilityatTRCpp(t*Q);
  }
  else{
  P <- expm(t*Q)
  }
  
  
  # return the row vector according to state at time t=0 if x0 is provided
  # if x0 not provided returns the whole matrix
  # hence returns probability of being at state j at time t if state at t =0 is x0
  # where j is from 1:noofStates
  if(missing(x0)){
    P <- matrix(P, nrow = NoofStates, dimnames = list(C@states,C@states))
    return(P)
  } else {
    if(x0 > NoofStates || x0 < 1){
      stop("Initial state provided is not correct")
    }
    return(P[x0,])
  }
}




#' Calculating full conditional probability using lower rate transition matrix
#' 
#' This function calculates full conditional probability at given 
#' time s using lower rate transition matrix
#' 
#' @usage impreciseProbabilityatT(C,i,t,s,error,useRCpp)
#' 
#' @param C a ictmc class object
#' @param i initial state at time t
#' @param t initial time t. Default value = 0
#' @param s final time
#' @param error error rate. Default value = 0.001
#' @param useRCpp logical whether to use RCpp implementation; by default TRUE
#' 
#' @references Imprecise Continuous-Time Markov Chains, Thomas Krak et al., 2016
#' 
#' @author Vandit Jain
#' 
#' @examples
#' states <- c("n","y")
#' Q <- matrix(c(-1,1,1,-1),nrow = 2,byrow = TRUE,dimnames = list(states,states))
#' range <- matrix(c(1/52,3/52,1/2,2),nrow = 2,byrow = 2)
#' name <- "testictmc"
#' ictmc <- new("ictmc",states = states,Q = Q,range = range,name = name)
#' impreciseProbabilityatT(ictmc,2,0,1,10^-3,TRUE)
#' @export
impreciseProbabilityatT <- function(C, i, t=0, s, error = 10^-3, useRCpp = TRUE){
  ##  input validity checking
  if(s <= t){
    stop("Please provide time points such that initial time is greater than or equal to end point")
  }
  
  if(!is(C,'ictmc')){
    stop("Please provide a valid ictmc-class object")
  }
  noOfstates <-length(C@states)
  
  if(i <= 0 || i > noOfstates){
    stop("Please provide a valid initial state")
  }
  ### validity checking ends
  if(useRCpp == TRUE) {
    Qgx <- .impreciseProbabilityatTRCpp(C,i,t,s,error)
  } else {
  ## extract values from ictmc object 
  Q <- C@Q
  range <- C@range
  
  ### calculate ||QI_i||
  
  #initialise Q norm value
  QNorm <- -1
  
  for(i in 1:noOfstates){
    sum <- 0
    for(j in 1:noOfstates){
      sum <- sum + abs(Q[i,j])
    }
    QNorm <- max(sum*range[i,2],QNorm)
  }
  
  ### calculate no of iterations
  # The 1 is for norm of I_s i.e. ||I_s|| which equals 1
  n <- max((s-t)*QNorm, (s-t)*(s-t)*QNorm*QNorm*1/(2*error))
  
  ### calculate delta
  delta <- (s-t)/n
  
  ### build I_i vector
  Ii <- rep(0, noOfstates)
  Ii[i] <- 1
  
  
  ### calculate value of lower operator _QI_i(x) for all x belongs to no ofStates
  values <- Q%*%Ii
  Qgx <- rep(0, noOfstates)
  for(i in 1:noOfstates){
    Qgx[i] <- min(values[i]*range[i,1], values[i]*range[i,2])
  }
  Qgx <- delta*Qgx
  Qgx <- Ii + Qgx
  for(iter in 1:n-1){
    temp <- Qgx
    values <- Q%*%Qgx
    for(i in 1:noOfstates){
      Qgx[i] <- min(values[i]*range[i,1], values[i]*range[i,2])
    }
    Qgx <- delta*Qgx
    Qgx <- temp + Qgx
  }
  }
  return(Qgx)
}



#' Check if CTMC is irreducible
#' 
#' @description 
#' This function verifies whether a CTMC object is irreducible
#' 
#' @usage is.CTMCirreducible(ctmc)
#' 
#' @param ctmc a ctmc-class object
#' 
#' @references 
#' Continuous-Time Markov Chains, Karl Sigman, Columbia University
#' 
#' @author Vandit Jain
#' 
#' @return a boolean value as described above.
#' 
#' @examples 
#' energyStates <- c("sigma", "sigma_star")
#' byRow <- TRUE
#' gen <- matrix(data = c(-3, 3,
#'                        1, -1), nrow = 2,
#'               byrow = byRow, dimnames = list(energyStates, energyStates))
#' molecularCTMC <- new("ctmc", states = energyStates, 
#'                      byrow = byRow, generator = gen, 
#'                      name = "Molecular Transition Model")
#' is.CTMCirreducible(molecularCTMC)
#' 
#' @export
is.CTMCirreducible <- function(ctmc) {
  
  if(!is(ctmc, 'ctmc') ) {
    stop("please provide a valid ctmc class object")
  }
  
  ## gets the embeded chain matrix
  embeddedChainMatrix <- generatorToTransitionMatrix(ctmc@generator)
  
  
  ## forms a markovchain object related to embedded transition matrix
  markovchainObject <- new("markovchain",states = ctmc@states,
                           transitionMatrix = embeddedChainMatrix)
  
  ## returns result using is.irreducible function on embedded chain transition matrix
  return(is.irreducible(markovchainObject))
  
}



#' checks if ctmc object is time reversible
#' 
#' @description 
#' The function returns checks if provided function is time reversible
#' 
#' @usage is.TimeReversible(ctmc)
#' 
#' @param ctmc a ctmc-class object
#' 
#' @return Returns a boolean value stating whether ctmc object is time reversible
#' 
#' @author Vandit Jain
#' 
#' @return a boolean value as described above
#' 
#' @references 
#' INTRODUCTION TO STOCHASTIC PROCESSES WITH R, ROBERT P. DOBROW, Wiley
#' 
#' @examples 
#' energyStates <- c("sigma", "sigma_star")
#' byRow <- TRUE
#' gen <- matrix(data = c(-3, 3,
#'                        1, -1), nrow = 2,
#'               byrow = byRow, dimnames = list(energyStates, energyStates))
#' molecularCTMC <- new("ctmc", states = energyStates, 
#'                      byrow = byRow, generator = gen, 
#'                      name = "Molecular Transition Model")
#' is.TimeReversible(molecularCTMC)
#' 
#' @export
is.TimeReversible <- function(ctmc) {
  
  if(!is(ctmc,"ctmc")) {
    stop("please provide a valid ctmc-class object")
  }
  
  ## get steady state probabilities
  Pi <- steadyStates(ctmc)
  
  ## initialise boolean result
  check <- TRUE
  
  ## no of states
  m <- length(ctmc@states)
  
  ## checks for byrow
  if(ctmc@byrow == FALSE)
    gen <- t(ctmc@generator)
  else
    gen <- ctmc@generator
  
  ## iterates for every state
  for( i in 1:m)
  {
    for(j in 1:m)
    {
      if(Pi[i]*gen[i,j] != Pi[j]*gen[j,i]) {
        check <- FALSE
        break
      }
    }
  }
  
  return(check)
  
}



# `generator/nextki` <- function(k) {
#   if(k >= 0) return(-1-k)
#   return(-k)
# }
# 
# `generator/nextk` <- function(k, Kmin, Kmax) {
#   if(is.null(k)) {
#     k <- rep(0, length(Kmin))
#     return(list(ans = TRUE, k = k))
#   }
#   
#   if(length(Kmin) == 0) {
#     return(list(ans = FALSE, k = k))
#   }
#   
#   i <- 1
#   kl <- k
#   kl[i] <- `generator/nextki`(kl[i])
#   while (kl[i] > Kmax[i] || kl[i] < Kmin[i]) {
#     kl[i] <- 0
#     i <- i+1
#     if(i > length(kl)) {
#       k <- kl
#       return(list(ans = FALSE, k = k))
#     }
#     kl[i] <- `generator/nextki`(kl[i])
#   }
#   k <- kl
#   return(list(ans = TRUE, k = k))
# }
# 
# `generator/generator` <- function(Po, Di, odi) {
#   P <- Po
#   N <- nrow(P)
#   
#   if(Di > 22) return(NULL) # bad idea
#   options(digits = Di)
#   odigs <- odi
#   
#   rSum <- rowSums(P)
#   if(! all(abs(1-rSum) < 0.001)) {
#     stop("Sum of each rows of Po should be equal to 1")
#   }
#   
#   P <- P/rSum
#   d <- det(P)
#   
#   if(d <= 0) {
#     cat("Matrix has non-positive determinant")
#     return(NULL)
#   }
#   
#   diagP <- 1
#   for(i in 1:nrow(P)) diagP <- diagP * P[i, i]
#   
#   if(d >= diagP) {
#     cat("Determinant exceeds product of diagonal elements\n")
#     return(NULL)
#   }
#   
#   E <- eigen(P)[[1]]
#   B <- eigen(P)[[2]]
#   
#   print("Eigenvalues")
#   print(E)
#   
#   # risky 
#   if(length(unique(E)) != length(E)) {
#     warning("Matrix does not have distinct eigenvalues")
#   }
#   
#   L <- abs(log(d))
#   addigs <- 2 + round(log10(1/Matrix::rcond(B))) + round(L/log(10)) # problem
#   
#   if(options()$digits < odigs + addigs) {
#     if(odigs + addigs > 100) {
#       print("Eigenvector matrix is singular")
#       return(NULL)
#     }
#     
#     cat('Going to', odigs + addigs, "digits")
#     return(`generator/generator`(Po, odigs + addigs, odigs))
#   }
#   
#   Bi <- solve(B)
#   
#   posevs <- NULL
#   negevs <- NULL
#   bestj <- NULL
#   bestQ <- NULL
#   marks <- rep(TRUE, length(E))
#   
#   for(i in 1:length(E)) { 
#     if(marks[i] && !(Re(E[i]) > 0 && Im(E[i]) == 0)) { # invalid comparison of complex number
#       cj <- Conj(E[i])
#       best <- Inf
#       if(i+1 <= length(E)) {
#         for(j in (i+1):length(E)) {
#          if(marks[j]) {
#            score <- abs(cj-E[j])
#            if(score < best) {
#              best <- score
#              bestj <- j
#            }
#          }
#         }
#       }
#         
#       if(best > 10^(3-options()$digits)) {
#         cat("Unpaired non-positive eigenvalue", E[i])
#         return(NULL)
#       }
#       marks[bestj] <- FALSE
#       if(Im(E[i]) >= 0) {
#         posevs <- c(posevs, i)
#         negevs <- c(negevs, bestj)
#         if(Im(E[bestj]) == 0) {
#           E[bestj] <- complex(real = E[bestj], imaginary = 0)
#         }
#       } else {
#         posevs <- c(posevs, bestj)  
#         negevs <- c(negevs, i)
#         if(Im(E[i]) == 0) {
#           E[i] <- complex(real = E[i], imaginary = 0)
#         }
#       }
#     }
#   }
#   
#   npairs <- length(posevs)
#   # display conjugate pairs
#   
#   Kmax <- rep(0, npairs)
#   Kmin <- Kmax
#   
#   for(i in 1:npairs) {
#     a <- Arg(E[posevs[i]])
#     Kmax[i] <- trunc((L-a)/2*pi)
#     Kmin[i] <- trunc((-L-a)/2*pi)
#   }
#   
#   # display K-max
#   # display K-min
#   
#   best <- -0.001
#   DD <- diag(log(E))
#   DK <- matlab::zeros(N)
#   res <- list(); p <- 1
#   k <- NULL
#   while(TRUE) {
#     
#     dlist <- `generator/nextk`(k, Kmin, Kmax)
#     k <- dlist$k
#     
#     if(dlist$ans == FALSE) {break}
#     
#     # display value of k
#     for(i in 1:npairs) {
#       ke <- complex(real = 0, imaginary = 2*pi*k[i])
#       DK[posevs[i], posevs[i]] <- ke
#       DK[negevs[i], negevs[i]] <- -ke
#     }
#     
#     Q <- B %*% (DD + DK) %*% Bi
#     # Q <- fnormal(Re(Q), options()$digits, 5*(10^(-1-odigs))) # define fnormal of maple
#     qmin <- Q[1,2]
#     for(i in 1:N) {
#       for(j in 1:N) {
#         if(i != j) {
#           if(Q[i, j] < qmin) qmin <- Q[i, j]
#         }
#       }
#     }
#     
#     if(EnvAllGenerators == TRUE) {
#       if(qmin > -.001) {
#         cat("Possible generator with qmin =", qmin)
#         res[[p]] <- round(Q, odigs)
#         p <- p + 1  
#       } else {
#         cat("qmin =", qmin)  
#       }
#       
#     } else {
#       if(qmin >= 0) {
#         cat("Found a generator")
#         return(round(Q, odigs))
#       } else {
#         if(qmin > best) {
#           best <- qmin
#           bestQ <- Q
#         }
#         if(qmin > -.001) {
#           cat("Approximate generator with qmin = ", qmin)
#         } else {
#           cat("qmin =", qmin)
#         }
#       }
#     }
#   }
#   
#   if(EnvAllGenerators == TRUE) {
#     return(res)
#   }
#   
#   warning("No completely valid generator found")
#   
#   if(! is.null(bestQ)) {
#     return(round(bestQ, odigs))
#   } else return(NULL)
#   
# }
# 
# generator <- function(Po, digits = 10) {
#   odigs <- digits
#   options(digits = 15)
#   if(is.matrix(Po)) {
#     P <- Po
#   } else {
#     stop("Po must be matrix")
#   }
#   
#   if(nrow(P) != ncol(P)) {
#     print(P)
#     stop('Po must be square matrix')
#   }
#   
#   if(! all(P >= 0)) {
#     print(P)
#     stop('Po must be non negative square matrix')
#   }
#   
#   `generator/generator`(P, options()$digits, odigs)
# }