version 1.0

DESCRIPTION

@@ -0,0 +1,30 @@
+Package: QPot
+Version: 1.0
+Date: 2015-10-22
+Title: Quasi-Potential Analysis for Stochastic Differential Equations
+Authors@R: c(
+	person("Christopher", "Moore", role = c("aut"), email = "life.dispersing@gmail.com"),
+	person("Christopher", "Stieha", role = c("aut", "cre"), email = "stieha@hotmail.com"),
+	person("Ben", "Nolting", role = c("aut"), email = ""),
+	person("Maria", "Cameron", role = c("aut"), email = ""),
+	person("Karen", "Abbott", role = c("aut"), email = ""),
+	person("James", "Gregson", role = ("cph"), email = "", comment = "author of expression_parser library: https://github.com/jamesgregson/expression_parser")
+	)
+Depends: R (>= 3.0.2)
+Imports: MASS
+Maintainer: Christopher Stieha <stieha@hotmail.com>
+Description: Tools to 1) simulate and visualize stochastic differential equations and 2) determine stability of equilibria using the ordered-upwind method to compute the quasi-potential.    
+License: GPL-2
+URL: http://www.r-project.org, https://github.com/bmarkslash7/QPot
+BugReports: https://github.com/bmarkslash7/QPot/issues
+NeedsCompilation: yes
+Packaged: 2015-10-23 02:06:16 UTC; stieha
+Author: Christopher Moore [aut],
+  Christopher Stieha [aut, cre],
+  Ben Nolting [aut],
+  Maria Cameron [aut],
+  Karen Abbott [aut],
+  James Gregson [cph] (author of expression_parser library:
+    https://github.com/jamesgregson/expression_parser)
+Repository: CRAN
+Date/Publication: 2015-10-23 08:50:26

MD5

@@ -0,0 +1,42 @@
+486b290357eebd501b76d4deb55be27d *DESCRIPTION
+f2aab4d8f56b75ecbc1374060aea63e4 *NAMESPACE
+f76207f3f595bc7121e1504a94a81dd6 *NEWS
+e22e22ee53ede1d0f148cf65f98ae393 *R/Model2String.R
+cbae596e7010d35585c4e8a884b33937 *R/QPContour.R
+dc6c74d10d9bdd674f28cdf84b4e6e8e *R/QPGlobal.R
+be81b23c728b092ca37e65db22393f37 *R/QPInterp.R
+6dd0af1f92df143e0ea46f2f611cfebf *R/QPotential.R
+c2f2763e5bb2ec53773e48cad43fe701 *R/TSDensity.R
+ce4db9fadf7df9237abf9d65f6addd37 *R/TSPlot.R
+68d9393221a304a02b148bcfe2b0c9c8 *R/TSTraj.R
+a27d1ba670bb6915dc74d6552b8f7fcd *R/VecDecomAll.R
+68b15decb9f3014597cd9b9a356f78ea *R/VecDecomGrad.R
+b07c8f92a695765e5b7c2f14906be61a *R/VecDecomPlot.R
+5c41a3568232f0d80412f09512ddb306 *R/VecDecomRem.R
+d3dc19290eb5da696bbeaee829c1fabe *R/VecDecomVec.R
+c70c8462933d65667a9f3916716b1da6 *R/test-default.R
+2d3995a70b02de03e92a581a6e76922f *README.md
+f038e36263be3e8501717d54d112627f *demo/00Index
+a052b33da402f04c8ec5c9d59ee34145 *demo/Ex1_test.R
+2c452a8a29d47762245da4ced108d7c3 *demo/Ex2_test.R
+deb551628bece73999ecf5fac8b9b0d5 *demo/Ex3_test.R
+a2af3e77eb4f5236cdab59a348e7687a *man/Model2String.Rd
+0676e8c0ea2cb9feb4380cf847bc326a *man/QPContour.Rd
+7ead564ac89fcf6db4c312192702814f *man/QPGlobal.Rd
+d54e29f68beb1299ae07d0ef18bd1a80 *man/QPInterp.Rd
+6df010710ec10546e271438680bd1023 *man/QPotential.Rd
+7bf470aacdbeabdbe4e1f601fe75b777 *man/TSDensity.Rd
+2405fd02e303b1ecf9a93d9f6af281d2 *man/TSPlot.Rd
+c0f2315ca2cef55386beb10e831033ca *man/TSTraj.Rd
+f1750050e6a266eb30bcf5d6d4c38ada *man/VecDecomAll.Rd
+a5dadf8b1235c8d09ce35b2a3d515d8f *man/VecDecomGrad.Rd
+abb3cf59fbb2a4c7b8c3a9a098d18409 *man/VecDecomPlot.Rd
+5ce6905fb9a82e2b0ee5134f6fae93f1 *man/VecDecomRem.Rd
+af7aa0b7ea5018f8953c68b42efffe51 *man/VecDecomVec.Rd
+3b3f9506d4149ca47bebf09e8bc31ef2 *man/defaultTest.Rd
+4884e5a52bce9acf6488c38b8976ebd7 *man/figures/Example3.png
+f349d3fa85852bfefe11ef0abbb8bd11 *src/Makevars
+6ced1da2b745886b816f9ff206f2edf3 *src/Makevars.win
+7f1c9d34bf1a028f338345cf522c3a34 *src/expression_parser.c
+7022188b5b6811cc0d3512c34ab8ddda *src/expression_parser.h
+b176e7c9d494b36fd4ebc6609fe84d02 *src/upwindorderedv4.c

NAMESPACE

@@ -0,0 +1,21 @@
+useDynLib(QPot, quasipotential)
+
+importFrom("grDevices", "colorRampPalette", "rgb")
+importFrom("graphics", ".filled.contour", "abline", "arrows", "axis",
+			"contour", "lines", "par", "plot", "polygon")
+importFrom("stats", "density", "rnorm")
+importFrom("utils", "read.table")
+
+export(Model2String)
+export(QPContour)
+export(QPGlobal)
+export(QPInterp)
+export(QPotential)
+export(TSDensity)
+export(TSPlot)
+export(TSTraj)
+export(VecDecomAll)
+export(VecDecomGrad)
+export(VecDecomRem)
+export(VecDecomPlot)
+export(VecDecomVec)

NEWS

@@ -0,0 +1,27 @@
+#########################################################################
+Install_instructions_self
+#########################################################################
+R CMD REMOVE QPot
+R CMD build QPot
+#copy tar.gz to /home to install
+#can't install across partitions
+R CMD INSTALL QPot_
+
+#################################
+examples taken from functions in R/
+##################################
+
+QPContour
+
+# #' @examples
+# #' # First, use a surface (example from QPGlobal)
+# #' global.qp <- QPGlobal(list(local.1,local.2),c(0,4),c(0,4),c(-1,5),c(-1,5))
+# #'
+# #' # Second, input that surface into QPContour
+# #' QPContour(surface=global.qp, density=c(100,100), xlim=c(-0.5,20), y.lim=c(-0.5,20), n.filled.contour=20, n.contour.lines=20, col=c("red", "white", "blue"), contour.lines = TRUE)
+
+
+QPGlobal
+# #' @export
+# #' @examples 
+# #' QPGlobal(list(local.1,local.2),c(0,4),c(0,4),c(-1,5),c(-1,5))

R/Model2String.R

@@ -0,0 +1,100 @@
+#' Converts equations with parameters to strings with numbers
+#'
+#' Function that converts differential equations from function-format to string-format.  Specifically, it reads in a function, searches for the differential equations within the function, and returns a list of strings containing the differential equations from the function.   Will also replace parameters of those equations with numerical values.
+#' 
+#' @param model.function function containing the differential equations as given to TSTraj()
+#' @param parms a named vector of paramters and their respective values for the deterministic equations.
+#' @param x.lhs.term string containing the left hand side of the first equation to search for, default is 'dx'
+#' @param y.lhs.term string containing the left hand side of the second equation to search for, default is 'dx'
+#' @param supress.print Default it FALSE, supress output.  TRUE prints out equations from function
+#' @return equations a list with two elements, the first is the x equation, the second is the y equation
+#'
+#' @examples
+#' #deSolve-style function call:
+#' model.parms <- c(alpha=1.54, beta=10.14, delta=1, kappa=1, gamma=0.476, mu=0.112509)
+#' ModelEquations <- function(t, state, parms) {
+#' 	with(as.list(c(state, parms)), {
+#'	dx <- (alpha*x)*(1-(x/beta)) - ((delta*(x^2)*y)/(kappa + (x^2)))
+#'	dy <- ((gamma*(x^2)*y)/(kappa + (x^2))) - mu*(y^2)
+#'	list(c(dx,dy))
+#'	})
+#' }
+#'
+#' Model2String(ModelEquations, parms = model.parms, x.lhs.term = 'dx', y.lhs.term = 'dy') 
+#'
+#' 
+#' #Second example with individual strings
+#' # Note the use of dx and dy
+#' test.eqn.x = "dx = (alpha*x)*(1-(x/beta)) - ((delta*(x^2)*y)/(kappa + (x^2)))"
+#' test.eqn.y = "dy = ((gamma*(x^2)*y)/(kappa + (x^2))) - mu*(y^2)"
+#' model.parms <- c(alpha=1.54, beta=10.14, delta=1, kappa=1, gamma=0.476, mu=0.112509)
+#' equations.as.strings.x <- Model2String(test.eqn.x, parms = model.parms, 
+#'  x.lhs.term = 'dx', y.lhs.term = 'dy')
+#' equations.as.strings.y <- Model2String(test.eqn.y, parms = model.parms, 
+#'  x.lhs.term = 'dx', y.lhs.term = 'dy')
+
+
+Model2String <- function(model.function, parms = 'NULL', x.lhs.term = 'dx', y.lhs.term = 'dy', supress.print = FALSE) {
+	if (!supress.print) {
+		print("Note: This function is supplied as duct tape.  Long equations, equations spanning multiple lines, equations with strange notation, etc, may not work.  Always check the output.")
+	}
+#	if (parms[1] == 'NULL') {stop("Need to define parms, the names and values of the model parameters")} 
+
+	temp <- deparse(model.function, width.cutoff = 500)
+# for some reason, inputing a string causes a " to be added
+	temp <- gsub(pattern = '\"', replacement = '', x = temp)
+
+#Dump function into a list of character strings
+#Go through each string and determine if it contains an equation
+	foundx = 0	#flag for making sure dx is only found once
+	foundy = 0	#flag for making sure dy is only found once
+
+#remove the lhs and return the rhs
+	for (i in 1:length(temp)) {
+	#when searching, first look for the lhs defining whether the derivative is for x or y
+	#once found, look inside the string and use either '<-' or '=' to separate
+	# the lhs from the rhs
+		if ((foundx == 0) && isTRUE(grepl(pattern = x.lhs.term, x = temp[i]))) {
+			foundx = 1
+			if (grepl(pattern = '<-', x = temp[i])) {
+				location <- regexpr(pattern = '<-', text = temp[i])
+				x.equation = substr(temp[i], start = (location+2), stop = nchar(temp[i]))
+			}
+			else if (grepl(pattern = '=', x = temp[i])) {
+				location <- regexpr(pattern = '=', text = temp[i])
+				x.equation = substr(temp[i], start = (location+2), stop = nchar(temp[i]))
+			} else {stop("Equation does not contain = or <-")}
+		}
+
+		if ((foundy == 0) && isTRUE(grepl(pattern = y.lhs.term, x = temp[i]))) {
+			foundy = 1
+			if (grepl(pattern = '<-', x = temp[i])) {
+				location <- regexpr(pattern = '<-', text = temp[i])
+				y.equation = substr(temp[i], start = (location+2), stop = nchar(temp[i]))
+			}
+			else if (grepl(pattern = '=', x = temp[i])) {
+				location <- regexpr(pattern = '=', text = temp[i])
+				y.equation = substr(temp[i], start = (location+2), stop = nchar(temp[i]))
+			} else {stop("Equation does not contain = or <-")}
+		}
+	}
+
+	equations = c()
+	if (foundx == 1) {equations = x.equation}
+	if (foundy == 1) {equations = c(equations, y.equation)}
+
+#if parameters are not declared, then we do not have to replace anything
+	if (!(parms[1] == 'NULL')){
+#replace the parameter names in the equations with their values
+		allnames <-names(parms)
+		for (i in 1:length(parms)) {
+			currname <- allnames[i]
+			value = toString(parms[[i]])
+			equations <- gsub(pattern = currname, replacement = value, x = equations)
+		}
+	} #if (!(parms[1] == 'NULL')){
+
+	if (!supress.print) {print(paste(equations))}
+	return(equations)
+
+}

R/QPContour.R

@@ -0,0 +1,91 @@
+#' Contour plot of quasi-potential surfaces
+#'
+#' This function allows you to plot quasi-potential surfaces
+#' @param surface the surface to be plotted, from \code{\link{QPGlobal}}.
+#' @param dens vector respectively for the number of \code{x} and \code{y} points to be plotted.
+#' @param x.bound two-element vector for the surface's minimum and maximum x values.
+#' @param y.bound two-element vector for the surface's minimum and maximum y values.
+#' @param xlim numeric vectors of length 2, giving the x coordinate range.
+#' @param ylim numeric vectors of length 2, giving the y coordinates range.
+#' @param n.filled.contour numeric value for the nubmber of breaks in the filled contour.
+#' @param n.contour.lines numeric value for the number of breaks in the contour lines.
+#' @param c.parm contour line adjustment (see details).
+#' @param col.contour a vector of colors used for the filled contour region.
+#' @param contour.lines if TRUE, then contour lines plotted over filled contour; vice versa if FALSE.
+#' @param ... passes arguments to both \code{\link{plot}} and \code{\link{contour}}.
+#' @details Because, in general, capturing the topological features can be subtle, we implemented a feature in \code{\link{QPContour}} to keep the filled.contour region while changing the contour lines.  Specifically, filled.contour takes the range of the surface values (\eqn{\phi}), divides by the number of the specified contours (i.e., \code{n.filled.contour}), and creates a contour at each break, which happenes to be equal across the range.  But because visualizing some topology may (i) require looking between contour breaks and (ii) adding contour lines would overload the plot with lines, we use an equation to modify the distribution of contour lines.  Namely, adjusting the \code{c} argument in the \code{\link{QPContour}} function adjusts the \eqn{c} paramter in the following equation: \deqn{max_\phi \times (\frac{x}{n-1})^c.}.  This allows the user to keep the same number of contour lines (i.e., specified with \code{n.contour.lines}), but focus them toward the troughs or peaks of the surfaces. At \eqn{c=1}, the contour lines correspond to the filled.contour breaks.  If \eqn{c > 1}, then the contour lines become more concentrated towards the trough.  Similarly, if \eqn{c < 1}, then the contour lines are more focused at the peaks of the surface.  As an example: \cr \figure{Example3.png}.
+#'
+#' @examples
+#' # First, System of equations
+#' 	equationx <- "1.54*x*(1.0-(x/10.14)) - (y*x*x)/(1.0+x*x)"
+#' 	equationy <- "((0.476*x*x*y)/(1+x*x)) - 0.112590*y*y"
+#' 
+#' # Second, shared parameters for each quasi-potential run
+#' 	xbounds <- c(-0.5, 10.0)
+#' 	ybounds <- c(-0.5, 10.0)
+#' 	xstepnumber <- 150
+#' 	ystepnumber <- 150
+#' 
+#' # Third, first local quasi-potential run
+#' 	xinit1 <- 1.40491
+#' 	yinit1 <- 2.80808
+#' 	storage.eq1 <- QPotential(x.rhs = equationx, x.start = xinit1, 
+#'		x.bound = xbounds, x.num.steps = xstepnumber, y.rhs = equationy,
+#'		y.start = yinit1, y.bound = ybounds, y.num.steps = ystepnumber)
+#' 
+#' # Fourth, second local quasi-potential run
+#' 	xinit2 <- 4.9040
+#' 	yinit2 <- 4.06187
+#' 	storage.eq2 <- QPotential(x.rhs = equationx, x.start = xinit2, 
+#'		x.bound = xbounds, x.num.steps = xstepnumber, y.rhs = equationy,
+#'		y.start = yinit2, y.bound = ybounds, y.num.steps = ystepnumber)
+#' 
+#' # Fifth, determine global quasi-potential 
+#' 	unst.x <- c(0, 4.2008)
+#' 	unst.y <- c(0, 4.0039)
+#' 	ex1.global <- QPGlobal(local.surfaces = list(storage.eq1, storage.eq2), 
+#'		unstable.eq.x = unst.x, unstable.eq.y = unst.y, x.bound = xbounds, 
+#'		y.bound = ybounds)
+#' 
+#' # Sixth, contour of the quasi-potential
+#' 	QPContour(ex1.global, dens = c(100,100), x.bound = xbounds, 
+#'		y.bound = ybounds, c.parm = 5)
+
+QPContour <- function(surface, dens, x.bound, y.bound, xlim = 'NULL', ylim = 'NULL', n.filled.contour=25, n.contour.lines=25, c.parm=1, col.contour, contour.lines = TRUE, ...){
+	x.range <- max(x.bound)-min(x.bound)
+	y.range <- max(y.bound)-min(y.bound)
+
+	row.range <- nrow(surface)-1
+	col.range <- ncol(surface)-1
+
+	row.min <- min(which(surface != 0 , arr.ind = T)[,1])
+	row.max <- max(which(surface != 0 , arr.ind = T)[,1])
+	col.min <- min(which(surface != 0 , arr.ind = T)[,2])
+	col.max <- max(which(surface != 0 , arr.ind = T)[,2])
+
+	x.min <- ((row.min-1)/row.range)*x.range + min(x.bound)
+	x.max <- ((row.max-1)/row.range)*x.range + min(x.bound)
+	y.min <- ((col.min-1)/col.range)*y.range + min(y.bound)
+	y.max <- ((col.max-1)/col.range)*y.range + min(y.bound)
+
+	if(missing(xlim)) {xlim <- c(x.min,x.max)}
+	if(missing(ylim)) {ylim <- c(y.min,y.max)}
+
+ 	sub.x <- seq(row.min, row.max, length.out=dens[1])
+	sub.y <- seq(col.min, col.max, length.out=dens[2])
+
+	sub.x.val <- ((sub.x-1)/row.range)*x.range + min(x.bound)
+	sub.y.val <- ((sub.y-1)/col.range)*y.range + min(y.bound)
+
+	eq.sub <- surface[sub.x, sub.y]
+
+	plot(0 , type = "n" , xlim = xlim , ylim = ylim , las = 1, ...)
+	min.eq.sub <- min(eq.sub , na.rm = T)
+	max.eq.sub <- max(eq.sub , na.rm = T)
+	contour.breaks <- seq(min.eq.sub , max.eq.sub , length = n.filled.contour)
+	eq.max <- max(surface, na.rm = T)
+	line.contour.breaks <- (eq.max)*(((0:n.contour.lines)/(n.contour.lines-1)))^c.parm
+	myRamp <- if(missing(col.contour)){colorRampPalette(c("#FDE725FF","#E3E418FF","#C7E020FF","#ABDC32FF","#8FD744FF","#75D054FF","#5DC963FF","#47C06FFF","#35B779FF","#28AE80FF","#20A486FF","#1F9A8AFF","#21908CFF","#24868EFF","#287C8EFF","#2C728EFF" ,"#31688EFF","#355D8DFF","#3B528BFF","#404688FF","#443A83FF","#472D7BFF","#481F71FF","#471163FF","#440154FF"))(n.filled.contour)}else{colorRampPalette(col.contour)(n.filled.contour)}
+	.filled.contour(sub.x.val , sub.y.val , eq.sub , levels = contour.breaks , col = myRamp)
+	if(contour.lines==TRUE){contour(sub.x.val , sub.y.val , eq.sub , levels = line.contour.breaks, drawlabels = F ,  add = TRUE , col = "black" , ...)}
+	}