inst/doc/QPot_Example2.R

In QPot: Quasi-Potential Analysis for Stochastic Differential Equations

###########################################################################
### This 'tangle' R script was created from an RSP document.
### RSP source document: 'QPot_Example2.md.rsp'
### Metadata 'keywords': ''
###########################################################################

##### preparation #####
	# 0.0.1 write a function to create a legend for contour plots
	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
	}
	# 0.0.2 preset color palettes for controus
	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")
require(QPot)
	require(R.devices)
	require(R.rsp)
	require(viridis)
	var.eqn.x <- "-(y-beta) + mu*(x-alpha)*(1-(x-alpha)^2-(y-beta)^2) "
	var.eqn.y <- "(x-alpha) + mu*(y-beta)*(1-(x-alpha)^2-(y-beta)^2)"
	model.parms <- c(alpha = 4, beta = 5, mu = 0.2)
	parms.eqn.x <- Model2String(var.eqn.x, parms = model.parms, supress.print = TRUE)
	parms.eqn.y <- Model2String(var.eqn.y, parms = model.parms, supress.print = TRUE)
require(phaseR)
model.ex2 <- 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 <- -(y-beta) + mu*(x-alpha)*(1-(x-alpha)^2-(y-beta)^2)
	dy <- (x-alpha) + mu*(y-beta)*(1-(x-alpha)^2-(y-beta)^2)
	list(c(dx,dy))
}
toPNG("EX2vectordiagram", {
	# draws the vector field
	flowField(deriv = model.ex2, xlim = c(2, 7), ylim = c(2, 7), 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.ex2, xlim = c(2, 7), ylim = c(2, 7), 
		parameters = model.parms, col = c("blue", "red"), points = 250)
	# draw and label the equilibria
	# open circles are unstable, black are stable
	points(4,5 , pch = 21 , col = "black" , bg = "white" , cex = 1.5)
	text(4,5 , labels = expression(italic(e[0])) , adj = c(0,-.25) , cex = 1.5)
	traj <- trajectory(model.ex2,y0=c(runif(1,2,7),runif(1,2,7)) , tlim = c(2,7), xlim = c(2,7) , ylim = c(2,7) , parameters = model.parms , t.end = 250, lwd = 1.5 , pch = 16)
	traj <- trajectory(model.ex2,y0=c(runif(1,2,7),runif(1,2,7)) , tlim = c(2,7), xlim = c(2,7) , ylim = c(2,7) , parameters = model.parms , t.end = 250, lwd = 1.5 , pch = 16)
	traj <- trajectory(model.ex2,y0=c(runif(1,2,7),runif(1,2,7)) , tlim = c(2,7), xlim = c(2,7) , ylim = c(2,7) , parameters = model.parms , t.end = 250, lwd = 1.5 , pch = 16)
	# label the x and y axis
	mtext(expression(italic(x)), side = 1, line = 2.5)
	mtext(expression(italic(y)), side = 2, line = 2.5)
})
model.state <- c(x = 3, y = 3)
	model.sigma <- 0.1
	model.time <- 500 #2500
	model.deltat <- 0.005
	rcat("\n")
ts.ex2 <- 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")
toPNG("EX2timeseriesFig1", {TSPlot(ts.ex2, deltat = model.deltat)})
toPNG("EX2timeseriesFig2", {TSPlot(ts.ex2, deltat = model.deltat, dim = 2)})
toPNG("EX2timeseriesFig3", {TSDensity(ts.ex2, dim = 1)})
k2 <- MASS::kde2d(ts.ex2[,2], ts.ex2[,3], n = 200)
	k2dns <- k2$z/sum(k2$z)
	k2cut <- cut(k2dns, 100, label = FALSE)
	crramp <- colorRampPalette(tsdens.col)
	colr <- crramp(100)
toPNG("EX2timeseriesFig4", {
		holdthese <- par(no.readonly = TRUE)
		par(mar = c(4, 4, 5 , 5))
		TSDensity(ts.ex2, dim = 2, xlab  ="", ylab = "", xlim = c(2.5, 6.5), ylim = c(2.5, 6.5), contour.levels = 25, contour.lines = T, las = 1, col2d = tsdens.col, contour.lwd = 0.25, kde2d.n = 200, xaxs = "i", yaxs = "i");
		legend.col(col = colr, lev = k2dns, xadj = 0.1)
		par(holdthese)
	})
bounds.x = c(-0.5, 11.5)
	bounds.y = c(-0.5, 11.5)
	step.number.x = 2000 # 4000
	step.number.y = 2000 # 4000
	xinit = 4.15611
	yinit = 5.987774
eq1.qp <- QPotential(x.rhs = parms.eqn.x, x.start = xinit, x.bound = bounds.x, 
					x.num.steps = step.number.x, y.rhs = parms.eqn.y, y.start = yinit, 
					y.bound = bounds.y, y.num.steps = step.number.y)
k2dns <- seq(min(eq1.qp, na.rm = T), max(eq1.qp, na.rm = T), 0.01)
		k2cut <- cut(k2dns, 100, label = FALSE)
		crramp <- colorRampPalette(qp.col)
		colr <- crramp(100)
toPNG("EX1globalvisualizationFig5", {
		holdthese <- par(no.readonly = TRUE)
		par(oma = c(0,1,0,2))
		QPContour(surface = eq1.qp , dens = c(1000, 1000), x.bound = bounds.x, y.bound = bounds.y, 
					c.parm = 1, xlab = expression(italic("x")), ylab = "");
		legend.col(col = colr, lev = k2dns, xadj = 0.1)
		par(holdthese)
	})
require(plot3D)
	
	#first, we have to subset the x and y axis because it is scaled from 0 to 1
	frac.x <- c(0.2,0.5)
    frac.y <- c(0.3,0.6)
    
    #then we reduce the global quasi-potential matrix to contain only these values
    ex1.global <- eq1.qp	
    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])]

	#regular data, can see the valley for the limit cycle
    dens.sub <- c(200, 200)	#pull only 200 rows and columns to speed up graphing
    global.sub <- global.sub[round(seq(1,nrow(global.sub),length.out=dens.sub[1])) , round(seq(1,ncol(global.sub),length.out=dens.sub[2]))]
    global.sub[global.sub > 0.01] <- NA # limit the z axis to give the best shot
toPNG("EX2globalvisualization3D", {
    persp3D(z = global.sub, theta = 25, phi = 55, col = viridis(100, option = "C"), shade = 0.1, 
        colkey = list(side = 4, length = 0.85), contour = list(levels = c(0.01, 0.001, 0.0001, 0.00005 ) ), zlim = c(-0.01, .011), 
        zlab = intToUtf8(0x03A6))
    #spend lots of time playing with theta, phi, zlim, contour levels, etc. to produce a decent graph
})
#natural log transformation
    global.sub <- log(global.sub)
toPNG("EX2visualization3Dln", {
    persp3D(z = global.sub, theta = 25, phi = 60, col = viridis(100, option = "C"), shade = 0.1, 
        colkey = list(side = 4, length = 0.85), contour = list(levels = seq(-10,-4,1) ), zlim = c(-15, -4), 
        zlab = intToUtf8(0x03A6)) # zlim = c(-10,-4)
})
VDAll <- VecDecomAll(surface = eq1.qp, x.rhs = parms.eqn.x, y.rhs = parms.eqn.y, x.bound = bounds.x, y.bound = bounds.y)
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, 11), arrow.type = "equal", 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)
})
toPNG("EX1vectorfielddecompREMprop", {
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("EX1vectorfielddecompREMequal", {
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.35, head.length = 0.025)
})
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)
})