R/dsmodel.R

In dsmodels: A Language to Facilitate Simulation and Visualization of Two-Dimensional Dynamical Systems

#' Defines a model object encapsulating a dynamical system
#'
#' To begin, define a two dimensional function that outputs
#' a two-dimensional list. For instance:
#' \code{fun <- function(x,y) list(x,y)}. Then
#' \code{dsmodel(fun)} will initialize the model.
#' \strong{Make sure that the function defining your model indeed outputs a \code{list}.}
#' The model is used to hold the data for graphics.
#' To display a desired graphic, add the corresponding
#' feature of type "dsproto" to your model.
#'
#' Models are constructed
#' incrementally using the + operator to add features
#' to the existing dsmodel object. A \code{\link{dsrange}} must be one of the objects added to a model.
#'
#' @section Methods:
#' \code{dsmodel} objects support the following methods, which may be helpful for advanced users.
#'
#' \code{dsmodel$points(format="list", filter="all")} returns a list of the points in the model.
#' \describe{
#' \item{\code{filter}}{Valid values are \code{"all"}, \code{"fixed"}, \code{"attractor"}, and \code{"sim"}.
#'  \code{"sim"} returns only points generated by \code{\link{simattractors}}. }
#' \item{\code{format}}{If \code{"list"}, return a list with \code{x}, \code{y}, \code{col}, and \code{ind} fields holding
#'  the coordinates, color, and index. Other formats are \code{"objects"}, returning a vector of \code{dspoint} objects,
#'  and \code{"pairs"}, returning a list of pairs of coordinates.}
#' }
#'
#' \code{dsmodel$display()} forces the model to re-render the plot from scratch. Primarily useful if \code{display=TRUE} was set.
#'
#' \code{dsmodel$basins()} returns a list of which fixed points have a basin. This requires \code{simbasins()} to have been composed
#' with the model, and is primarily useful when testing if a dynamical system is globally stable. In that case, the method
#' will return a list of length 1. The list will contain the indices of the fixed points, as given in
#' \code{dsmodel$points(format="list", filter="attractor")}. An index of 0 means that some points never moved within
#' epsilon of an attractor.
#'
#' \code{dsmodel$sim.is.stable()} attempts to determine if the system is stable by simulation. If no attractors have
#' been composed with the model, \code{simattractors()} is composed with defaults. If \code{simbasins} has not
#' been composed with the model, it is be composed with defaults. If every point is drawn to a single attractor, the
#' system has been deemed stable. Note that boundary points on the range will not be tested.
#'
#' @family Foundation
#' @param fun Function with two inputs and two outputs which defines the dynamical system. The output should be a list, preferably with field names x and y.
#' @param title A string title for the graph. Text can be input in the form of pseudo-LaTeX code within quotes.
#'  See \code{\link[latex2exp]{TeX}} for more details.
#' @param display If set to \code{FALSE}, the model will be drawn only when the user calls \code{`MODELNAME`$display()}. Otherwise,
#'  the model will be drawn with every \code{dsmodels} object added after the \code{dsrange} is added.
#' @import grDevices latex2exp
#' @include dsproto.R
#' @seealso \code{\link{dsrange}}
#' @export
#' @examples
#' library(dsmodels)
#'
#' fun <- function(X,Y) {
#'   list(
#'     x = X/exp(Y),
#'     y = Y/exp(X)
#'   )
#' }
#' # Add dsRange to see the actual range.
#' model <- dsmodel(fun)
#'
#' dsmodel(function(x,y) {
#'   list(
#'     x = x^2,
#'     y = x/(y+1)
#'   )
#' }, title = "Another function showing $f(x)=x^{\\alpha}$!")
dsmodel <- function(fun, title="", display = TRUE) {
  texTitle <- TeX(title)
  if(length(formals(fun)) != 2)
    stop("dsmodel: Please make sure your function has 2 distinct, variable inputs.")
  dsproto(
    `_class` = "model",
    `_inherit` = NULL,
    fun = fun,
    title = texTitle,
    dim = 2,
    range = NULL,
    background = c(),
    feature = c(),
    visualization = c(),
    dots = c(),
    autoDisplay = display,
    print = function(self, ...) {
      invisible(self)
    },
    apply = function(self, x, y, iters=1, accumulate=TRUE, crop = TRUE) {
      if(is.null(x) || is.null(y))
        stop("dsmodel: Please make sure your x and y values are defined in latest object created.")
      tmp = self$fun(x[1], y[1])
      if(length(tmp) != 2)
        stop("dsmodel: Please make sure your function outputs a list with 2 values. (for example: list(x+1,y^2)")
      if(is.element("x", names(tmp)) && is.element("y", names(tmp)))
        properNames = TRUE
      else {
        properNames = FALSE
        if(!is.null(names(tmp)))
          warning("dsmodel function has outputs names that are not \"x\" and \"y\". Assuming first output is x, and second is y.")
      }
      if(accumulate) {
        iterAccum = vector("list",iters+1)
        startvals = list(x=x, y=y)
        if(crop){
          startvals <- self$cropframe(startvals)
        }
        iterAccum[[1]] <- startvals
        iterSeq = 1:iters
        if(iters==0)
          iterSeq <- NULL
        for(i in iterSeq){
          tmp=self$fun(x,y)
          if(any(is.nan(tmp[[1]])) || any(is.nan(tmp[[2]])))
          {
            warning("dsmodel: model undefined, NaN computed. (Division by zero). Removing points with NaN value and continuing procedure..")
            tmp = NaNRemove(tmp)
          }
          if(crop){
            tmp = self$cropframe(tmp)
          }
          if(!properNames)
            names(tmp) <- c("x","y")
          x = tmp$x
          y = tmp$y
          iterAccum[[i+1]] <- tmp
        }
        iterAccum
      } else {
        for(i in 1:iters){
          tmp=self$fun(x,y)
          if(crop){
            tmp <- self$cropframe(tmp)
          }
          if(!properNames)
            names(tmp) <- c("x","y")
          x = tmp$x
          y = tmp$y
          x = tmp[[1]]
          y = tmp[[2]]
        }
        if(!properNames)
          names(tmp) <- c("x","y")
        tmp
      }
    },
    cropframe = function(self, vals) {
      xin <- vals[[1]]
      yin <- vals[[2]]
      tmp1 <- self$range$xlim
      tmp2 <- self$range$ylim
      filterfunx <- function(x) (x <= max(tmp1) && x >= min(tmp1))
      filterfuny <- function(x) (x <= max(tmp2) && x >= min(tmp2))
      xrid = boolpos(filterfunx, xin)
      yrid = boolpos(filterfuny, yin)
      rid <- c(xrid,yrid)
      if(!is.null(rid)) {
        list(
          x = xin[-rid],
          y = yin[-rid]
        )
      }
      else {
        list(
          x = xin,
          y = yin
        )
      }
    },
    render = function(self, obj = NULL) {
      if(self$autoDisplay){
        rerender = FALSE
        if(is.null(obj))
          rerender = TRUE
        if(is.range(obj)) {
          if(!is.range(self$range))
            stop("Range not added properly: severe error. Please notify developers.")
          if(!self$range$rendered)
            rerender = TRUE
        }
        else if (!is.null(self$range)) {
          if(is.feature(obj) ||
            (is.visualization(obj) && is.null(self$feature)) ||
            (is.background(obj) && is.null(self$visualization)  && is.null(self$feature))) {
            obj$render(model = self)
          } else {
            rerender = TRUE
          }
        }
        if(rerender) {
			self$display()
		}
      }
    },
    recalculate = function(self) {
      if(!is.null(self$range)) {
        if(is.range(self$range)) {
          self$range$render(model = self)
        }
        if(!is.null(self$background)) {
          for(ba in self$background)
            ba$recalculate(model = self)
        }
        if(!is.null(self$visualization)) {
          for(vi in self$visualization)
            vi$recalculate(model = self)
        }
        if(!is.null(self$feature)) {
          for(fe in self$feature)
            fe$recalculate(model = self)
        }
      }
    },
    points = function(self, format="list", filter="all") {
      points <- Filter(is.dspoint, self$feature)
      if(filter == "attractor")
        points <- Filter(function(p) {p$attractor}, points)
      else if (filter == "sim")
        points <- Filter(function(p) {p$artificial}, points)
      else if (filter == "fixed" || filter == "fixedpoint")
        points <- Filter(function(p) {p$fixed}, points)
      else if(filter != "all")
        stop("dsmodel: valid filters for dsmodel$points method are: \"all\", \"attractor\", \"fixed\", \"sim\".")
      if(format == "objects")
        points
      else if (format == "list") {
        res <- pointsToList(points)
        res$inds = c(res$inds,recursive=TRUE)
        res
      }
      else if (format == "pairs")
        Map(function(pnt) c(pnt$x, pnt$y), points)
      else
        stop("dsmodel: valid formats for dsmodel$points method are: \"objects\",  \"list\", \"pairs\"")
    },
    display = function(self){
    	self$range$render(model = self)
        for(bg in self$background)
          bg$render(model = self)
        for(vi in self$visualization)
          vi$render(model = self)
        for(fe in self$feature)
          fe$render(model = self)
    },
		basins = function(self){
      basins <- Filter(is.dsimage,c(self$background,self$feature))
      if(length(basins) == 0)
        stop("dsmodel$basins requires simbasins() to have been composed with the model.")
      if(length(basins) > 1)
        stop("dsmodel$basins can't handle models with multiple simbasins().")
      basin <- basins[[1]]
      if(!basin$regionsCalculated)
        basin$recalculate(self)
      res <- unique(c(basin$colMatrix))
      if(is.element(0,res))
        warning("dsmodel$basins: simbasins did not find an attractor for every point.")
      res
		},
		sim.is.stable = function(self) {
		  attractors <- Filter(
		    function(x) {
		      if(is.dspoint(x))
		        x$attractor
		      else
		        FALSE
		    }, self$feature)
		  if(length(attractors) == 0)
		    self + simattractors()
		  basins <- Filter(is.dsimage,c(self$background,self$feature))
		  if(length(basins) == 0)
		    self$simbasins()
		  if(length(basins) > 1)
		    stop("dsmodel$basins can't handle models with multiple simbasins().")
		  basin <- basins[[1]]
		  if(!basin$regionsCalculated)
		    basin$recalculate(self)
		  res <- unique(c(basin$colMatrix))
		  (length(res) == 1) && !(is.element(0,res))
		}
  )
}
#' Checks if object is a mode.
#'
#' @param x Object to check
#' @export
#' @keywords internal
#' @usage NULL
is.model <- function(x) inherits(x,"model")

#' Set Color Vector From Parameters
#'
#' Outputs a color vector.
#' If the combined length of \code{col} and \code{image} is  2, and \code{iters} is greather than 1,
#' this function creates a gradient from first specified color to the second. The length of said
#' vector will be equal to \code{iters}.
#' If the combined length of \code{col} and \code{image} is greater than two, then the output will
#' be a vector of whatever was input into \code{col} and \code{image} with output length being
#' \code{min(length(col) + length(image), iters)}.
#' @param col A vector or list of colors
#' @param image A vector or list of colors
#' @param iters A number indicating the number of iterations desired.
#'  If set to greater than 1, \code{iters} becomes the length of the output.
#' @keywords internal
#' @export
colorVector <- function(col, image, iters) {
  if(is.null(image))
    image <- ""
  if((length(col) == 2 && iters > 1 && image == "") # GRADIENT
     || (length(col) == 1 && length(image) == 1 && iters > 1)){
    if(!is.null(image)) {
      if(image != "")
        col <- append(col,image)
    }
    colf <- colorRamp(col)
    colptr <- seq(0,1, length.out = iters+1)
    vect <- vector(mode = "list", iters+1)
    for(i in 1:(iters+1)) {
      tmp <- colf(colptr[[i]])
      vect[[i]] = rgb(red = tmp[1],
                green = tmp[2],
                blue = tmp[3], maxColorValue = 255)
    }
    vect
  }
  else{ # NOT GRADIENT
    if(image != "" || is.null(image))
      col <- append(col, image)
    lengthOfCol <- length(col)
    len <- min(lengthOfCol, iters)
    if( iters <= 1 ) len <- lengthOfCol
    vect <- vector(mode = "list", len)
    for(i in 1:len){
      if(col[[i]]=="")
        col[[i]] = "NA"
      vect[[i]] <- col[[i]]
    }
    vect
  }
}

#' Find Positions of False Elements
#'
#' Returns all positions of elements in \code{v} where \code{fun(v)}
#' returns \code{FALSE}.
#' @param fun A \code{function(x)} which outputs a single boolean
#' @param v A vector or list of elements compatible with \code{fun}
#' @keywords internal
#' @export
boolpos <- function(fun, v) {
  holder <- c()
  tracker <- 1
  if(length(v) <= 0)
    return(holder)
  for(i in 1:length(v)) {
    val <- fun(v[[i]])
    if(!is.na(val) && !is.nan(val)) {
      if(!val) {
        holder[[tracker]] <- i
        tracker <- tracker+1
      }
    }
  }
  holder
}

#' NaNRemove
#'
#' Finds and removes all 2-d points which produce NaNs
#' @keywords internal
#' @param twoDList Data of form \code{list(list(),list())}.
#' @export
NaNRemove <- function(twoDList){
  tmp <- twoDList
  if(length(tmp[[1]])!=length(tmp[[2]])){
    stop("dsmodel: X and Y of different length. Internal Error.")
  }
  ridFun <- function(x) !is.nan(x)
  toRidX <- boolpos(ridFun,tmp[[1]])
  toRidY <- boolpos(ridFun,tmp[[2]])
  rid <- c(toRidX,toRidY)
  if(!is.null(rid)) {
    list(
      tmp[[1]] <- tmp[[1]][-rid],
      tmp[[2]] <- tmp[[2]][-rid]
    )
  }
}


#' Abstract a Function which does not Crash Upon Failure
#'
#' This function takes a boolean function and abstracts it. If the given function fails, the program will not crash.
#' @keywords internal
#' @param fun Input function to abstract.
#' @param inp Input to that function.

safe.apply <- function(fun,inp){
  tryCatch({fun(inp)},
           warning=function(w) { FALSE},
           error=function(e) { FALSE })
}