R/simulation_functions.R

In WarrenDavidAnderson/dynamicNetworkID:

#' Initial conditions
#'
#' This function extracts the initial conditions for the dynamic variables corresponding to measurements.
#' @param datHMF The data matrix with scaled measurements measurements. Scaled center measurements are provided for 
#' each measurement/annotation (e.g., gene/organ) combination. This format is produced by scale_zeroOne(). Defaults to datHMF.
#' @return An initial condition vector.
#' @export
#' @examples ic = get_init_condits(datHMF)
get_init_condits = function(datHMF=datHMF){
  ic = rep(0, length(unique(datHMF[,2]))*(ncol(datHMF)-2))
  x = 1
  for(ii in 1:length(unique(datHMF[,2]))){
    ic_ind_ann = which(datHMF[,2] == unique(datHMF[,2])[ii])
    for(jj in 1:(ncol(datHMF)-2)){
      ic_ind_feat = which(names(datHMF) == names(datHMF)[3:ncol(datHMF)][jj])
      ic[x] = datHMF[ic_ind_ann,ic_ind_feat][1]
      x = x + 1
    }
  }
  return(ic)
} # end get_init_condits


#' Dynamic model
#'
#' This function generates the mathematical formulation for the dynamic model.
#' @param time Vector of simulation times.
#' @param y Dynamic variables.
#' @param params Parameter matrix, the Jacobian.
#' @return Model object for input into the ode() function from the deSolve library.
#' @examples out <- ode(y=ic, times=time, parms=J, func=odemodel, method=lsoda)
#' @references \url{https://www.jstatsoft.org/article/view/v033i09}
odemodel <- function(time, y, params) {
  dy <- apply(params, 1, function(x) sum(y*x))
  list(c(dy))
}


#' Simulate the dynamic model and record the agreement with the data.
#'
#' This function generates simulations and corresponding errors.
#' @param datHMF The data matrix with scaled measurements measurements. Scaled center measurements are provided for 
#' each measurement/annotation (e.g., gene/organ) combination. This format is produced by scale_zeroOne(). Defaults to datHMF.
#' @param Nr Number of annotations (e.g., organs).
#' @param tlen Duration of simulation time.
#' @param rateConst Rate constant matrix.
#' @param ic Initial conditions.
#' @param dt Time step.
#' @param epsilon Term for computing error.
#' @param numPar Number of parameters.
#' @return Simulation errors.
#' @examples MLAM[indLam,5] <- sim_error(datHMF,Nr,tlen,rateConst,ic,dt,epsilon,numPar)
#' @references \url{https://www.jstatsoft.org/article/view/v033i09}
# model simulation - get errors
sim_error = function(datHMF=NULL,Nr=NULL,tlen=NULL,rateConst=NULL,ic=NULL,dt=NULL,epsilon=NULL,numPar=NULL){
  
  errors = rep(0,ncol(rateConst))
  for(sims in 1:ncol(rateConst)){
    
    par = rateConst[,sims]
    
    # Jacobian matrix, J
    param_num = length(par)
    J = matrix(c(0),sqrt(param_num),sqrt(param_num))
    for (ii in 1:sqrt(numPar)) {
      ind <- ((ii-1)*sqrt(param_num)+1):(ii*sqrt(param_num))
      J[ii,] <- par[ind]
    } # ii (parameter loop)
    
    # implement similation 
    time <- seq(0,tlen,by=dt)
    out <- deSolve::ode(y=ic, times=time, parms=J, func=odemodel)
    
    # loop through annotation/feature combinations to compute errors
    sim_time = out[,1] + 4
    err = c()
    for(ii in 1:Nr){
      indorg = which(datHMF[,2] == unique(datHMF[,2])[ii])
      orgcols = ((ii-1)*22+1):(22*ii) + 1
      expr_time = as.numeric(datHMF[indorg,1])
      ind_sim_time = unlist( sapply(expr_time,function(x){which(sim_time==x)}) )
      sum_sq_err = apply( ( datHMF[indorg,3:ncol(datHMF)] - out[ind_sim_time,orgcols] )^2, 2, sum )
      var_term = apply(out[,orgcols],2,var) + epsilon
      err = c( err, sum ( sum_sq_err / var_term ) )
    }
    errors[sims] = sum(err)
  } # sims loop
  return(errors)
} # end sim_error


#' Format error
#'
#' This function generates reformats the simulation error data from a list to a data.frame.
#' @param simulation_error Simulation errors obtained from storing the results of sim_error() in MLAM.
#' @param alph Vector of alpha terms (L1 regularization parameter).
#' @return Simulation error data.frame.
#' @examples simulate_error = format_error(simulation_error,alph)
format_error = function(simulation_error=NULL,alph=alph){
  mat = c()
  for(ii in 1:length(simulation_error)){
    mat0 = simulation_error[[ii]]
    aa = rep(alph[ii], nrow(mat0))
    mat = rbind(mat, cbind(aa,mat0))
  }
  colnames(mat)[1] = "alpha_value"
  return(as.data.frame(mat))
}


#' Best simulation
#'
#' This function gets the parameters (lambda, alpha, m-values) for the best simulation based on 
#' the minumum of the simulation error.
#' @param simulate_error simulation error data frame generated by format_error().
#' @return Best simulation parameters.
#' @examples best_params = select_best(simulate_error)
select_best = function(simulate_error=NULL){
  ind = which(simulate_error$error == min(simulate_error$error))
  out = simulate_error[ind,]
  rownames(out)=""
  return(out)
}



#' Implement simulation
#'
#' This function performs a simulation for a single set of simulation parameters.
#' @param mRange_best A set of m-values.
#' @param tsim simulation time vector.
#' @param numPar Number of simulation parameters.
#' @param XY Object generated by getXY().
#' @param datHMF The data matrix with scaled measurements measurements. Scaled center measurements are provided for 
#' each measurement/annotation (e.g., gene/organ) combination. This format is produced by scale_zeroOne(). 
#' @param ic Initial conditions.
#' @return Simulation results
#' @export
#' @examples 
#' best_sims = sim_best(best_params,mRange_best,tsim,numPar,XY,datHMF,ic)
#' best_sims$simulations[,1] = best_sims$simulations[,1] + min(datHMF[,1])
sim_best = function(best_params=NULL,mRange_best=NULL,tsim=NULL,numPar=NULL,XY=NULL,datHMF=NULL,ic=NULL){
  
  XY = XY
  numPar = numPar
  Yreg = get_Yreg(datHMF,XY$Y,numPar,mRange_best)
  Xreg = get_Xreg(datHMF,XY$X,numPar,mRange_best)
  fit <- glmnet::glmnet( data.matrix(Xreg), data.matrix(Yreg), family="gaussian",intercept=F,standardize=F,
                 alpha = best_params$alpha_value, lambda = best_params$lambda_value)
  tsim = tsim
  par = fit$beta
  #numPar = length(par)
  
  # Jacobian matrix, J
  param_num = length(par)
  J = matrix(c(0),sqrt(param_num),sqrt(param_num))
  for (ii in 1:sqrt(numPar)) {
    ind <- ((ii-1)*sqrt(param_num)+1):(ii*sqrt(param_num))
    J[ii,] <- par[ind]
  } # ii (parameter loop)
  
  # implement similation 
  sim <- deSolve::ode(y=ic, times=tsim, parms=J, func=odemodel)
  
  # annotate simulation and parameter data
  anns = get_annotations(sim,J)
  
  return(list(simulations=anns[[1]], param_matrix=anns[[2]]))
}



#' Format simulation annotations
#'
#' This function formats the final simulation data. This function is utilized by sim_best().
#' @param sim simulations.
#' @param J Jacobian matrix of parameters
#' @return list(sim,J)
#' @examples best_params = select_best(simulate_error)
get_annotations = function(sim=NULL,J=NULL){
  orgs = unique(datHMF[,2])
  assays = colnames(datHMF)[3:ncol(datHMF)]
  anns = as.vector( sapply(orgs,function(x){paste0(x,"_",assays)}) )
  colnames(sim)[2:ncol(sim)] = anns
  colnames(J) = anns
  rownames(J) = anns
  return(list(sim,J))
}