changes for new BacArena 1.7 cran release

- minor errors and fixed to pass cran checks
- HeatMapFeeding() was quite buggy...

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BacArena
 Title: Modeling Framework for Cellular Communities in their Environments
-Version: 1.6
+Version: 1.7
 Authors@R: c( person("Eugen", "Bauer", , "[email protected]", role = c("aut")), 
               person("Johannes", "Zimmermann", , "[email protected]", role = c("aut", "cre")))
 Author: Eugen Bauer [aut],
@@ -11,7 +11,7 @@ Description: Can be used for simulation of organisms living in
     metabolic models determine the uptake and release of compounds. Biological
     processes such as movement, diffusion, chemotaxis and kinetics are available
     along with data analysis techniques.
-URL: https://github.com/euba/BacArena
+URL: https://BacArena.github.io/
 BugReports: https://github.com/euba/BacArena/issues
 Depends:
     R (>= 3.0.0),
@@ -31,7 +31,6 @@ Imports:
     plyr,
     Rcpp
 Suggests:
-    sybilSBML,
     parallel,
     knitr,
     rmarkdown

NEWS

@@ -1,3 +1,17 @@
+BacArena v.1.7 (Release date: 2018-03-21)
+==============
+Changes: 
+- fine-tuning on growth model
+   * several duplications per time step
+   * fixed no growth bug for very high accumulated biomass
+- improved handling of non-biomass objectives
+- better status printing
+- more then one compound in chemotaxis
+- wrong unit of diffusion constant
+- several small updated in plotting methods
+
+
+
 BacArena v.1.6 (Release date: 2017-03-21)
 ==============
 Changes: 

R/Arena.R

@@ -337,13 +337,13 @@ setMethod("addSubs", "Arena", function(object, smax=0, mediac=object@mediac, dif
     if( !template )
       diffmat <- conv * diffmat
     else
-      smax    <- smax * (arena@n * arena@m) / sum(diffmat!=0)
+      smax    <- smax * (object@n * object@m) / sum(diffmat!=0)
   }
   if( is(diffmat, "list") ){
     if( !template )  
       diffmat <- lapply(diffmat, function(x){conv*x})
     else
-      smax    <- sapply(seq_along(smax), function(i){smax[i] * (arena@n * arena@m) / sum(diffmat[[i]]) })
+      smax    <- sapply(seq_along(smax), function(i){smax[i] * (object@n * object@m) / sum(diffmat[[i]]) })
   }
 
 
@@ -497,7 +497,7 @@ setMethod("addEssentialMed", "Arena", function(object, org, only_return=FALSE, l
   min_val <- -ex@lowbnd[which(ex_max<0)] # get lower bounds also for -INF cases
   min_val[min_val==Inf] <- 1000
 
-  predefined <- sapply(arena@media, function(x){sum(x@diffmat)}) # only use compounds which are not already in present in the media
+  predefined <- sapply(object@media, function(x){sum(x@diffmat)}) # only use compounds which are not already in present in the media
   idx <- which(min_id %in% names(predefined)[which(predefined==0)])
 
   if(length(idx)==0) return(object)
@@ -1473,6 +1473,7 @@ setMethod(show, "Arena", function(object){
 #' @slot mfluxlist A list of containing highly used metabolic reactions per time step. 
 #' @slot shadowlist A list of containing shadow prices per time step. 
 #' @slot subchange A vector of all substrates with numbers indicating the degree of change in the overall simulation.
+#' @slot exchangeslist A list of containing exchanges per time step. 
 setClass("Eval",
          contains="Arena",
          representation(
@@ -2850,6 +2851,7 @@ setMethod("checkCorr", "Eval", function(object, corr=NULL, tocheck=list()){
 #' @param spec_nr Number of the specie
 #' @param sub_nr Maximal number of substances to be show
 #' @param cutoff Shadow costs should be smaller than cutoff
+#' @param noplot Do not plot
 #' @details Returns ggplot objects
 setGeneric("plotShadowCost", function(object, spec_nr=1, sub_nr=10, cutoff=-1, noplot=FALSE){standardGeneric("plotShadowCost")})
 #' @export
@@ -3089,46 +3091,42 @@ setMethod("plotSubDist2", "Eval", function(object, sub, times=NULL){
 #' @return Heatmap (ggplot2)
 #' 
 #' @examples 
-#'  Do Not Run
-#'  sim : simulation object
-#'  speciesA = 1 # E. Coli (wild): sequence number in names(sim@specs) => 1
-#'  speciesB = 3 # E. Coli (mutant): sequence number in names(sim@specs) => 3
-#'  var_nr = 30 # (manually selected)
-#'  HeatMapFeeding(object = sim, speciesA = 1, speciesB = 3, var_nr = 30)
+#' sim <- sihumi_test
+#' HeatMapFeeding(object = sim, speciesA = 1, speciesB = 2, var_nr = 78)
 #' 
 setGeneric("HeatMapFeeding", function(object, speciesA, speciesB, var_nr){standardGeneric("HeatMapFeeding")})
 #' @export
 #' @rdname HeatMapFeeding
 setMethod("HeatMapFeeding", "Eval", function(object, speciesA, speciesB, var_nr){
-  A<-BacArena::plotSpecActivity(object,spec_list = speciesA, var_nr = var_nr, ret_data = T)
-  z<-data.frame()
-  chronos <- (seq_along(object@simlist)-1)
-  for (t in chronos) 
-  { a<-data.frame()
-  a <-BacArena::findFeeding3(object, time = t, mets = A$sub[1:var_nr], plot = F)
-  if (nrow(a)!=0) 
-  {z<-rbind.data.frame(z,a)
-  }
+  A <- BacArena::plotSpecActivity(object,spec_list = speciesA, var_nr = var_nr, ret_data = T)
+  z <- data.frame()
+  chronos <- 1:(length(object@simlist)-1)
+  for (t in chronos) { 
+    a<-data.frame()
+    a <-BacArena::findFeeding3(object, time = t, mets = A$sub[1:var_nr], plot = F)
+    if (nrow(a)!=0)
+      z <- rbind.data.frame(z,a)
   }
   d <- data.frame()
   # prod:speciesA cons:speciesB   -> 1                       
   # prod:speciesB cons:speciesA   -> -1
   # NA -> 0
-  for (h in 1:nrow(z))
-  {if (z[h,1] == names(object@specs)[speciesA]) (d[h,1]=1)
-    if (z[h,1] == names(object@specs)[speciesB]) (d[h,1]=-1)
+  for (h in 1:nrow(z)){
+    if (z[h,1] == names(object@specs)[speciesA]) 
+      d[h,1] <- 1
+    else if (z[h,1] == names(object@specs)[speciesB]) 
+      d[h,1] <- -1
+    else
+      d[h,1] <- NA
   }
   colnames(d)[1] <- "status"
   q <- cbind(z,d)
-  p <- reshape::cast(q, met ~ sim_step, value = "status")
-  p[is.na(p)] <- 0
-  # GGPLOT METHOD based on https://stackoverflow.com/questions/8406394/how-to-produce-a-heatmap-with-ggplot2
-  p2 <- p
-  p.m <- reshape::melt(p2)
-  p3 <- ggplot2::ggplot(p.m, ggplot2::aes(sim_step,met)) +
-    ggplot2::geom_tile(ggplot2::aes(fill = value), colour = "white") +
+  q$status[is.na(q$status)] <- 0
+  p3 <- ggplot2::ggplot(q, ggplot2::aes_string("sim_step","met")) +
+    ggplot2::geom_tile(ggplot2::aes_string(fill = "status"), colour = "white") +
     ggplot2::scale_fill_gradient2(name = "Producer", low = "red", mid = "green", high = "blue", breaks=seq(-1,1,by=1),
                                   labels = c(names(object@specs)[speciesB], "no feeding", names(object@specs)[speciesA])) +
     ggplot2::xlab("Simulation Step") + ggplot2::ylab("Exchange Reactions") 
-  return((p3)) })
+  return((p3)) 
+})
 

R/Organism.R

@@ -120,7 +120,7 @@ Organism <- function(model, algo="fba", ex="EX_", ex_comp=NA, csuffix="\\[c\\]",
     if(length(pot_biomass)==0) stop("No objection function set in model")
     print("No objective function set, try set automatically...")
     print(paste("found:", pot_biomass))
-    print(paste("set new biomass function:", pot_biomas[1]))
+    print(paste("set new biomass function:", pot_biomass[1]))
     sybil::obj_coef(model)[which(sybil::react_id(model)==pot_biomass[1])] <- 1
     rbiomass <- pot_biomass[1]
   }else{
@@ -392,9 +392,10 @@ setMethod("optimizeLP", "Organism", function(object, lpob=object@lpobj, lb=objec
     shadow <- glpkAPI::getRowsDualGLPK(lpob@problem@oobj)[ex@met_pos]
     names(shadow) <- ex@met_id
   }else if(lpob@problem@solver=="cplexAPI" & solve_ok){ # TODO: cplex shadow costs needs update of optimzeProb call in optimizeLP() which is having side effects ...
-    ex <- findExchReact(object@model)
-    shadow <- cplexAPI::getPiCPLEX(lpob@problem@oobj@env, lpob@problem@oobj@lp, 0, lpob@nr-1)[ex@met_pos]
-    names(shadow) <- ex@met_id
+    warning("cplex shadow costs are not supported yet")
+    #ex <- findExchReact(object@model)
+    #shadow <- cplexAPI::getPiCPLEX(lpob@problem@oobj@env, lpob@problem@oobj@lp, 0, lpob@nr-1)[ex@met_pos]
+    #names(shadow) <- ex@met_id
   }
 
   names(fbasl$fluxes) <- names(object@lbnd)
@@ -822,7 +823,7 @@ setMethod("growth_par", "Bac", function(object, population, j, fbasol, tstep){
 #' arena <- Arena(n=20,m=20) #initialize the environment
 #' arena <- addOrg(arena,bac,amount=10) #add 10 organisms
 #' arena <- addSubs(arena,40) #add all possible substances
-#' chemotaxis(bac,arena,1)
+#' chemotaxis(bac,arena,1, "EX_glc(e)")
 setGeneric("chemotaxis", function(object, population, j, chemo){standardGeneric("chemotaxis")})
 #' @export
 #' @rdname chemotaxis
@@ -924,13 +925,14 @@ setMethod("simBac", "Bac", function(object, arena, j, sublb, bacnum, sec_obj="no
 #' @param lpobject linar programming object (copy of organism@lpobj) that have to be a deep copy in parallel due to pointer use in sybil.
 #' @param sec_obj character giving the secondary objective for a bi-level LP if wanted.
 #' @param cutoff value used to define numeric accuracy
+#' @param with_shadow True if shadow cost should be stores (default off).
 #' @return Returns the updated enivironment of the \code{population} parameter with all new positions of individuals on the grid and all new substrate concentrations.
 #'
-setGeneric("simBac_par", function(object, arena, j, sublb, bacnum, lpobject, sec_obj="none", cutoff=1e-6){standardGeneric("simBac_par")})
+setGeneric("simBac_par", function(object, arena, j, sublb, bacnum, lpobject, sec_obj="none", cutoff=1e-6, with_shadow=FALSE){standardGeneric("simBac_par")})
 #' @export
 #' @rdname simBac_par
-setMethod("simBac_par", "Bac", function(object, arena, j, sublb, bacnum, lpobject, sec_obj="none", cutoff=1e-6){
-  const <- constrain(object, object@medium, lb=-sublb[j,object@medium]/bacnum, ub=object@ubnd*bacnum, #scale to population size
+setMethod("simBac_par", "Bac", function(object, arena, j, sublb, bacnum, lpobject, sec_obj="none", cutoff=1e-6, with_shadow=FALSE){
+  const <- constrain(object, object@medium, lb=-sublb[j,object@medium]/bacnum, #scale to population size
                      dryweight=arena@orgdat[j,"biomass"], tstep=arena@tstep, scale=arena@scale)
   lobnd <- const[[1]]; upbnd <- const[[2]]
   fbasol <- optimizeLP(object, lb=lobnd, ub=upbnd, j=j, sec_obj=sec_obj, cutoff=cutoff, with_shadow=with_shadow)

R/Stuff.R

@@ -172,6 +172,7 @@ lsd <- function(y){lb=mean(y)-stats::sd(y); ifelse(lb<0,0,lb)}
 #' @param ret_data Set true if data should be returned
 #' @param num_var Number of varying substances to be shown (if mediac is not specified)
 #' @param unit Unit for the substances which should be used for plotting (default: mmol)
+#' @param useNames Use substance names instead of ids
 #' 
 #' @return list of three ggplot object for further formating
 #'