<%@meta language="R-vignette" content="--------------------------------
  DIRECTIVES FOR R:

  %\VignetteIndexEntry{Example 1: A consumer-resource model with alternative stable states}
  %\VignetteAuthor{Christopher Stieha, Christopher Moore}
  %\VignetteEngine{R.rsp::rsp}
--------------------------------------------------------------------"%>

<%
		holdthese <- par(no.readonly = TRUE)
		tsdens.col <- c("lightsteelblue", "white", "indianred")
		qp.col <- 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")
		legend.col <- function(col, lev, xadj){ 
				opar <- par
				n <- length(col)
				bx <- par("usr")
				box.cx <- c(bx[2] + (bx[2] - bx[1]) / 1000 + 0, bx[2] + (bx[2] - bx[1]) / 1000 + (bx[2] - bx[1]) / 50)
				box.cy <- c(bx[3], bx[3])
				box.sy <- (bx[4] - bx[3]) / n
				xx <- rep(box.cx, each = 2) + xadj
				par(xpd = TRUE)
				for(i in 1:n){
					yy <- c(box.cy[1] + (box.sy * (i - 1)),
					box.cy[1] + (box.sy * (i)),
					box.cy[1] + (box.sy * (i)),
					box.cy[1] + (box.sy * (i - 1)))
					polygon(xx, yy, col = col[i], border = col[i])
				}
				par(new = TRUE)
				plot(0, 0, type = "n", ylim = c(min(lev), max(lev)), yaxt = "n", ylab = "",
				xaxt = "n", xlab = "", frame.plot = FALSE)
				axis(side = 4, las = 2, tick = FALSE, line = (.25 + xadj))
				par <- opar
			}
%>


In this vignette, we will extend the example from the R article and include the code to produce the figures.

### Example 1: A consumer-resource model with alternative stable states
<br>
```R
	var.eqn.x <- "(alpha*x)*(1-(x/beta)) - ((delta*(x^2)*y)/(kappa+(x^2)))"
	var.eqn.y <- "((gamma*(x^2)*y)/(kappa+(x^2))) - mu*(y^2)"
	model.parms <- c(alpha = 1.54, beta = 10.14, delta = 1, gamma = 0.476, kappa = 1, mu = 0.112509)
	parms.eqn.x <- Model2String(var.eqn.x, parms = model.parms)
	parms.eqn.y <- Model2String(var.eqn.y, parms = model.parms) 
```
	
<%
	require(QPot)
	require(R.devices)
	require(R.rsp)
	var.eqn.x <- "(alpha*x)*(1-(x/beta)) - ((delta*(x^2)*y)/(kappa+(x^2)))"
	var.eqn.y <- "((gamma*(x^2)*y)/(kappa+(x^2))) - mu*(y^2)"
	model.parms <- c(alpha = 1.54, beta = 10.14, delta = 1, gamma = 0.476, kappa = 1, mu = 0.112509)
	parms.eqn.x <- Model2String(var.eqn.x, parms = model.parms)
	rcat("\n")
	parms.eqn.y <- Model2String(var.eqn.y, parms = model.parms) 
	rcat("\n")
%>
<br>
Here, we confirm that Model2String() correctly combined the equations and the parameters.  We can suppress this output from Model2String() using suppress.print = TRUE.

### Step 1: Analyze the deterministic skeleton
<br>
To analyze our system of differential equations, we use the package phaseR, which requires the system of equations to be implemented as a function.  Our function Model2String() was made to help users of QPot transition from their use of previous packages requiring this function format.

```R
	require(phaseR)
	model.ex1 <- function(t, y, parameters){
		x <- y[1]
		y <- y[2]
		alpha <- parameters["alpha"]
		beta <- parameters["beta"]
		delta <- parameters["delta"]
		kappa <- parameters["kappa"]
		gamma <- parameters["gamma"]
		mu <- parameters["mu"]
		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))
	}
```
<%
	require(phaseR)
	model.ex1 <- function(t, y, parameters){
		x <- y[1]
		y <- y[2]
		alpha <- parameters["alpha"]
		beta <- parameters["beta"]
		delta <- parameters["delta"]
		kappa <- parameters["kappa"]
		gamma <- parameters["gamma"]
		mu <- parameters["mu"]
		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))
	}
%>


```R
	# draws the vector field
	flowField(deriv = model.ex1, xlim = c(-0.1, 6), ylim = c(-0.1, 6), parameters = model.parms, 
		add = F, points = 30, col = "grey70", ann = F, arrow.head = 0.025, frac = 1.1, 
		xaxs = "i", yaxs = "i", las = 1)
	# draw the nullclines and suppress the output by assigning it to a variable
	supp.print <- nullclines(deriv = model.ex1, xlim = c(-0.1, 6), ylim = c(-0.1, 6), 
		parameters = model.parms, col = c("blue", "red"), points = 250)
	# draw and label the equilibria
	# open circles are unstable, black are stable
	points(0,0 , pch = 21 , col = "black" , bg = "white" , cex = 1.5)
	text(0,0 , labels = expression(italic(e[u1])) , adj = c(-0.1,-.25) , cex = 1.5)
	points(1.4049, 2.8081 , pch = 21 , col = "black" , bg = "grey30" , cex = 1.5)
	text(1.4049, 2.8081 , labels = expression(italic(e[s1])) , adj = c(-0.1,1.25) , cex = 1.5)
	points(4.2008, 4.0039 , pch = 21 , col = "black" , bg = "white" , cex = 1.5)
	text(4.2008, 4.0039 , labels = expression(italic(e[u2])) , adj = c(-0.1,1.25) , cex = 1.5)
	points(4.9040, 4.0619 , pch = 21 , col = "black" , bg = "grey30" , cex = 1.5)
	text(4.9040, 4.0619 , labels = expression(italic(e[s2])) , adj = c(.4,-.35) , cex = 1.5)
	# label the x and y axis
	mtext(expression(italic(x)), side = 1, line = 2.5)
	mtext(expression(italic(y)), side = 2, line = 2.5)
```
<%toPNG("EX1vectordiagram", {
	# draws the vector field
	flowField(deriv = model.ex1, xlim = c(-0.1, 6), ylim = c(-0.1, 6), parameters = model.parms, add = F, points = 30, col = "grey70", ann = F, arrow.head = 0.025, frac = 1.1, xaxs = "i", yaxs = "i", las = 1)
	# draw the nullclines and suppress the output by assigning it to a variable
	supp.print <- nullclines(deriv = model.ex1, xlim = c(-0.1, 6), ylim = c(-0.1, 6), parameters = model.parms, col = c("blue", "red"), points = 250)
	# draw and label the equilibria
	# open circles are unstable, black are stable
	points(0,0 , pch = 21 , col = "black" , bg = "white" , cex = 1.5)
	text(0,0 , labels = expression(italic(e[u1])) , adj = c(-0.1,-.25) , cex = 1.5)
	points(1.4049, 2.8081 , pch = 21 , col = "black" , bg = "grey30" , cex = 1.5)
	text(1.4049, 2.8081 , labels = expression(italic(e[s1])) , adj = c(-0.1,1.25) , cex = 1.5)
	points(4.2008, 4.0039 , pch = 21 , col = "black" , bg = "white" , cex = 1.5)
	text(4.2008, 4.0039 , labels = expression(italic(e[u2])) , adj = c(-0.1,1.25) , cex = 1.5)
	points(4.9040, 4.0619 , pch = 21 , col = "black" , bg = "grey30" , cex = 1.5)
	text(4.9040, 4.0619 , labels = expression(italic(e[s2])) , adj = c(.4,-.35) , cex = 1.5)
	mtext(expression(italic(x)), side = 1, line = 2.5)
	mtext(expression(italic(y)), side = 2, line = 2.5)
})
%>

![Null clines of the current system of equations](figures/EX1vectordiagram.png)

### Step 2: Stochastic simulation
<br>
	In order to simulate the stochastic equations, we need:
	
```R
	model.state <- c(x = 1, y = 2) 	# initial condition
	model.sigma <- 0.05				# the level of noise
	model.time <- 1000 #12500				# the length of the simulation
	model.deltat <- 0.025			# the time step
```
<%
	model.state <- c(x = 1, y = 2) 	# initial condition
	model.sigma <- 0.05				# the level of noise
	model.time <- 1000 #12500				# the length of the simulation
	model.deltat <- 0.025			# the time step
	rcat("\n")
%>
<br>
	Then we can simulate the equations using TSTraj():
	
```R
	ts.ex1 <- TSTraj(y0 = model.state, time = model.time, 
		deltat = model.deltat, x.rhs = var.eqn.x, y.rhs = var.eqn.y, 
		parms = model.parms, sigma = model.sigma)
```
<%	
	ts.ex1 <- TSTraj(y0 = model.state, time = model.time, 
		deltat = model.deltat, x.rhs = var.eqn.x, y.rhs = var.eqn.y, 
		parms = model.parms, sigma = model.sigma)
	rcat("\n")
%>

The function TSPlot allows us to easily view the output from TSTraj().
The default plot is the time series plotted through time with a histogram:
```R
	TSPlot(ts.ex1, deltat = model.deltat)
```

<%toPNG("EX1timeseriesFig1", {TSPlot(ts.ex1, deltat = model.deltat)})%>
![Stochastic simulation of a time series](figures/EX1timeseriesFig1.png)

We can view the times series plotted along on the state variables by setting dim = 2: 
```R
	TSPlot(ts.ex1, deltat = model.deltat, dim = 2)
```
<%toPNG("EX1timeseriesFig2", {TSPlot(ts.ex1, deltat = model.deltat, dim = 2)})%>
<img src="figures/EX1timeseriesFig2.png">
	
If we only want to see the histograms, we can see them in one dimensional space
```R
	TSDensity(ts.ex1, dim = 1)
```
<%toPNG("EX1timeseriesFig3", {TSDensity(ts.ex1, dim = 1)})%>
<img src="figures/EX1timeseriesFig3.png">

or in two-dimensional space with 
```R
	TSDensity(ts.ex1, dim = 2)
```
<%
		k2 <- MASS::kde2d(ts.ex1[,2], ts.ex1[,3], n = 200)
		k2dns <- k2$z/sum(k2$z)
		k2cut <- cut(k2dns, 100, label = FALSE)
		crramp <- colorRampPalette(tsdens.col)
		colr <- crramp(100)
%>

<%toPNG("EX1timeseriesFig4", {
	par(mar = c(4, 4, 4 , 4))
	TSDensity(ts.ex1, dim = 2)
	legend.col(col = colr, lev = k2dns, xadj = 0.1)
	par(holdthese)
	})
%>
<img src="figures/EX1timeseriesFig4.png">


### Step 3: Local quasi-potential calculations
<br>
To calculate the quasi-potential for a system of equations, we first need to define some initial conditions and parameter values.

```R	
	bounds.x = c(-0.5, 20.0) 	# upper and lower limit of X
	bounds.y = c(-0.5, 20.0) 	# upper and lower limit of Y
	step.number.x = 1000 #4100 	# number of division between upper and lower limit
	step.number.y = 1000 #4100 	# number of division between upper and lower limit
	eq1.x = 1.40491				# x value of first equilibrium
	eq1.y = 2.80808				# y value of first equilibrium
	eq2.x = 4.9040				# x value of second equilibrium
	eq2.y = 4.06187				# y value of second equilibrium
```
<%
	bounds.x = c(-0.5, 20.0)
	bounds.y = c(-0.5, 20.0)
	step.number.x = 1000 # 4100
	step.number.y = 1000 # 4100
	eq1.x = 1.40491
	eq1.y = 2.80808
	eq2.x = 4.9040
	eq2.y = 4.06187
%>


With these values defined, we can now compute the quasi-potential around each equilibrium point.  Note that we don't have to start at the equilibria, but this ensures that we quickly find the lowest quasi-potential and work up from there.
```R
	eq1.local <- QPotential(x.rhs = parms.eqn.x, x.start = eq1.x, x.bound = bounds.x, 
					x.num.steps = step.number.x, y.rhs = parms.eqn.y, y.start = eq1.y,  
					y.bound = bounds.y, y.num.steps = step.number.y)
	eq2.local <- QPotential(x.rhs = parms.eqn.x, x.start = eq2.x, x.bound = bounds.x, 
					x.num.steps = step.number.x, y.rhs = parms.eqn.y, y.start = eq2.y, 
					y.bound = bounds.y, y.num.steps = step.number.y)
```
<%
	eq1.local <- QPotential(x.rhs = parms.eqn.x, x.start = eq1.x, x.bound = bounds.x, x.num.steps = step.number.x, y.rhs = parms.eqn.y, y.start = eq1.y,  y.bound = bounds.y, y.num.steps = step.number.y)
	eq2.local <- QPotential(x.rhs = parms.eqn.x, x.start = eq2.x, x.bound = bounds.x, x.num.steps = step.number.x, y.rhs = parms.eqn.y, y.start = eq2.y, y.bound = bounds.y, y.num.steps = step.number.y)
%>

We can use the function QPContour() to view the compute quasi-potential around each equilibrium.  In some situations, this may be useful.  For this example, it is not and we skip straight to computing the global quasi-potential and visualizing that.

### Step 4: Global quasi-potential calculation
<br>
QPGlobal() takes a list of local.surfaces (in this current example, two) and combines the output of QPotential() into a single global quasi-potential.  We supply two points to help the function (in this case unstable equilibria) to tell QPGlobal() which points should have the same value.  QPGlobal() uses these points to recalibrate the local quasi-potentials and combines these into a global quasi-potential. 

```R
	ex1.global <- QPGlobal(local.surfaces = list(eq1.local, eq2.local), 
					unstable.eq.x = c(0, 4.2008), unstable.eq.y = c(0, 4.0039), 
					x.bound = bounds.x, y.bound = bounds.y)
```

<%
	ex1.global <- QPGlobal(local.surfaces = list(eq1.local, eq2.local), unstable.eq.x = c(0, 4.2008), unstable.eq.y = c(0, 4.0039), x.bound = bounds.x, y.bound = bounds.y)
%>

### Step 5: Global quasi-potential visualization

```R
		QPContour(surface = ex1.global, dens = c(1000, 1000), x.bound = bounds.x, y.bound = bounds.y, c.parm = 1)
		QPContour(surface = ex1.global, dens = c(1000, 1000), x.bound = bounds.x, y.bound = bounds.y, c.parm = 5)
```

<%
	k2dns <- seq(min(ex1.global, na.rm = T), max(ex1.global, na.rm = T), 0.001)
	k2cut <- cut(k2dns, 100, label = FALSE)
	crramp <- colorRampPalette(qp.col)
	colr <- crramp(100)

%>


<%toPNG("EX1globalvisualizationFig5", {
		QPContour(surface = ex1.global, dens = c(1000, 1000), x.bound = bounds.x, y.bound = bounds.y, 
					c.parm = 1, xlab = expression(italic("x")), ylab = "")
		mtext(text = expression(italic("y")), side = 2, line = 2.5, outer = T, at = 0.5)
	})
%>

<%toPNG("EX1globalvisualizationFig6", {
		par(oma = c(0,1,0,2))
		QPContour(surface = ex1.global, dens = c(1000, 1000), x.bound = bounds.x, y.bound = bounds.y, 
					c.parm = 5, xlab = expression(italic("x")), ylab = "")
		legend.col(col = colr, lev = k2dns, xadj = 0.1)
		par(holdthese)
	})
%>


The parameter c.parm modifies the distribution of contour lines throughout the quasi-potential (see the help file for specifics).  But for right now, we can compare c.parm = 1
![Global quasi-potential visualization](figures/EX1globalvisualizationFig5.png)

with c.parm = 5

![Global quasi-potential visualization](figures/EX1globalvisualizationFig6.png)


In case you want to view the quasipotential in three dimensions, we can use the R package plot3D:

```R
	require(plot3D)
	require(viridis)
	frac.x <- c(0.05, 0.3)
	frac.y <- c(0.125, 0.25)
	global.sub <- ex1.global[round(dim(ex1.global)[1]*frac.x[1]):round(dim(ex1.global)[1]*frac.x[2]),
		round(dim(ex1.global)[2]*frac.y[1]):round(dim(ex1.global)[2]*frac.y[2])]
	global.sub[global.sub > 0.01] <- NA
#	par(fin = c(4, 4), oma = rep(0,4), mar = c(3.5, 3.5, 0, 0))
	persp3D(z = global.sub, theta = 25, phi = 25, col = viridis(100, option = "C"), shade = 0.1, 
		colkey = list(side = 4, length = 0.85), contour = T, zlim = c(-0.001, .011), 
		zlab = intToUtf8(0x03A6))
```

<%
	require(plot3D)
	require(viridis)
toPNG( "EX1globalvisualization3D", {
	frac.x <- c(0.05, 0.3)
	frac.y <- c(0.125, 0.25)
	global.sub <- ex1.global[round(dim(ex1.global)[1]*frac.x[1]):round(dim(ex1.global)[1]*frac.x[2]), round(dim(ex1.global)[2]*frac.y[1]):round(dim(ex1.global)[2]*frac.y[2])]
	global.sub[global.sub > 0.01] <- NA
#	par(fin = c(4, 4), oma = rep(0,4), mar = c(3.5, 3.5, 0, 0))
	persp3D(z = global.sub, theta = 25, phi = 25, col = viridis(100, option = "C"), shade = 0.1, colkey = list(side = 4, length = 0.85), contour = T, zlim = c(-0.001, .011), zlab = intToUtf8(0x03A6))
})
%>

![Three dimensional quasi-potential visualization using package plot3D](figures/EX1globalvisualization3D.png)

<%---
```R
	library(rgl)
	dens.sub <- c(4000,4000)
	global.sub <- ex1.global[round(seq(1,nrow(ex1.global),length.out=dens.sub[1])) ,
		round(seq(1,ncol(ex1.global),length.out=dens.sub[2]))]
	# global.sub[global.sub > 0.02] <- NA # to cut off large values for a better view of the basin
	persp3d(x = global.sub, col = "orange", expand = 1.1, xlim = c(0.05, 0.35), ylim = c(0.1, 0.3), 
		zlim = c(0, 0.01), xlab = "X", ylab = "Y", zlab = intToUtf8(0x03A6))
```

<%
	library(rgl)
	dens.sub <- c(4000,4000)
	global.sub <- ex1.global[round(seq(1,nrow(ex1.global),length.out=dens.sub[1])) , round(seq(1,ncol(ex1.global),length.out=dens.sub[2]))]
	# global.sub[global.sub > 0.02] <- NA # to cut off large values for a better view of the basin
%>

![Three dimensional quasi-potential visualization using package rgl](figures/EX1globalvisualizationFig7.png)

---%>


### Step 6: Vector field decomposition
<br>
	Taking the global quasi-potential, we can visualize the vector field, specifically the gradient and the remainder.  First we calculate the vector field decomposition, which gives us the deterministic skeleton, the gradient vector field, and the remainder vector field for both state variables X and Y. 

```R
	VDAll <- VecDecomAll(surface = ex1.global, x.rhs = parms.eqn.x, 
				y.rhs = parms.eqn.y, x.bound = bounds.x, y.bound = bounds.y)
```
<%
VDAll <- VecDecomAll(surface = ex1.global, x.rhs = parms.eqn.x, y.rhs = parms.eqn.y, x.bound = bounds.x, y.bound = bounds.y)
%>


We can then plot the deterministic skeleton
```R
	VecDecomPlot(x.field = VDAll[,,1], y.field = VDAll[,,2], dens = c(25, 25), 
			x.bound = bounds.x, y.bound = bounds.y, xlim = c(0, 11), ylim = c(0, 6), 
			arrow.type = "equal", tail.length = 0.25, head.length = 0.025)
```
<%toPNG("EX1vectorfielddecompDETSKEL", {
VecDecomPlot(x.field = VDAll[,,1], y.field = VDAll[,,2], dens = c(25, 25), x.bound = bounds.x, y.bound = bounds.y, xlim = c(0, 11), ylim = c(0, 6), arrow.type = "equal", tail.length = 0.25, head.length = 0.025)
})
%>
![Global quasi-potential visualization](figures/EX1vectorfielddecompDETSKEL.png)


followed by the gradient vector field

```R
VecDecomPlot(x.field = VDAll[,,3], y.field = VDAll[,,4], dens = c(25, 25), x.bound = bounds.x, 
		y.bound = bounds.y, arrow.type = "proportional", tail.length = 0.25, head.length = 0.025)
```
<%toPNG("EX1vectorfielddecompGRAD", {
VecDecomPlot(x.field = VDAll[,,3], y.field = VDAll[,,4], dens = c(25, 25), x.bound = bounds.x, y.bound = bounds.y, arrow.type = "proportional", tail.length = 0.25, head.length = 0.025)
})
%>
![Global quasi-potential visualization](figures/EX1vectorfielddecompGRAD.png)


and the remainder vector field, which is the "force that causes trajectories to circulate around level sets of the quasi-potential."
```R
VecDecomPlot(x.field = VDAll[,,5], y.field = VDAll[,,6], dens = c(25, 25), x.bound = bounds.x, 
		y.bound = bounds.y, arrow.type = "proportional", tail.length = 0.35, head.length = 0.025)
```
<%toPNG("EX1vectorfielddecompREM", {
VecDecomPlot(x.field = VDAll[,,5], y.field = VDAll[,,6], dens = c(25, 25), x.bound = bounds.x, y.bound = bounds.y, arrow.type = "proportional", tail.length = 0.35, head.length = 0.025)
})
%>
![Global quasi-potential visualization](figures/EX1vectorfielddecompREM.png)




The force arrows can be drawn equal size or can be drawn proportional to their value by setting arrow.type to "equal" or "proportional".

<%toPNG("EX1vectorfielddecompFig8", {
			holdthese <- par(no.readonly = TRUE)
			par(mfrow = c(2,2), mar = c(2,2,1,1), oma = c(3,3,1,1))
			VecDecomPlot(x.field = VDAll[,,3], y.field = VDAll[,,4], dens = c(25, 25), x.bound = bounds.x, y.bound = bounds.y, arrow.type = "proportional", tail.length = 0.25, head.length = 0.025)
			mtext(side = 3, text = expression(Gradient~field~(-nabla~phi(x, y))))
			mtext(side = 2, text = "Proportional arrow lengths", line = 3.75)
			mtext(side = 2, text = expression(italic(y)), line = 2.25)
			VecDecomPlot(x.field = VDAll[,,5], y.field = VDAll[,,6], dens = c(25, 25), x.bound = bounds.x, y.bound = bounds.y, arrow.type = "proportional", tail.length = 0.35, head.length = 0.025)
			mtext(side = 3, text = expression(Remainder~field~(bold(r)(x, y))))
			VecDecomPlot(x.field = VDAll[,,3], y.field = VDAll[,,4], dens = c(25, 25), x.bound = bounds.x, y.bound = bounds.y, arrow.type = "equal", tail.length = 0.15, head.length = 0.025)
			mtext(side = 2, text = "Equal arrow lengths", line = 3.75)
			mtext(side = 2, text = expression(italic(y)), line = 2.25)
			mtext(side = 1, text = expression(italic(x)), line = 2.25)	
			VecDecomPlot(x.field = VDAll[,,5], y.field = VDAll[,,6], dens = c(25, 25), x.bound = bounds.x, y.bound = bounds.y, arrow.type = "equal", tail.length = 0.15, head.length = 0.025)
			mtext(side = 1, text = expression(italic(x)), line = 2.25)
			par(holdthese)
})
%>
![Global quasi-potential visualization](figures/EX1vectorfielddecompFig8.png)
