Permalink
Browse files

moved warningsignals file into independent R package

  • Loading branch information...
1 parent df10c2c commit 45bb7cd409a6311f173e8e2b94bb55460417ee14 @cboettig committed Dec 2, 2010
View
@@ -0,0 +1,11 @@
+Package: warningsignals
+Title: Detection of Early Warning Signals of a Catastrophic Bifurcations
+Version: 0.0-1
+Date: 2010-12-01
+Author: Carl Boettiger <cboettig@gmail.com>
+Maintainer: Carl Boettiger <cboettig@gmail.com>
+Description:
+Depends: snowfall, NORMT3, sde, Kendall, odesolve
+URL: http://www.carlboettiger.info
+License: GPL (>= 2)
+
File renamed without changes.
File renamed without changes.
File renamed without changes.
@@ -0,0 +1,95 @@
+# saddle_analytics.R
+params <- c(K=1000, e=.5, c=1)
+b <- function(x, params){ x*params["c"]*(params["K"]-x) }
+
+## "Changing Mean" Dynamics
+d <- function(x, e, params){ e*x}
+E<- seq(0,1.5, length=50)
+x <- seq(0, 1300, by=5)
+
+require(rootSolve)
+xhat <- function(e, params, interval=c(0,1000)){
+ F <- function(x) {b(x, params) - d(x, e, params) }
+ uniroot.all(F, interval=interval)
+ }
+
+bifur <- function(){
+ for(e in E){
+ plot(x, b(x,params)-d(x,e, params), type="l", ylim=c(-150, 900), ylab="rate")
+ lines(x, b(x,params), lwd=4, lty=1, col="darkblue")
+ lines(x, d(x,e, params), lwd=4, lty=1, col="darkred")
+ abline(h=0, lwd=4, lty=2, col="darkgrey")
+ crit_pts <- sort(xhat(e, params))
+# print(crit_pts)
+ if(length(crit_pts) == 1){
+ points(crit_pts, 0, pch=19, cex=1.5)
+ } else if(length(crit_pts > 1)){
+ points(min(crit_pts), 0, cex=1.5)
+ points(max(crit_pts), 0, pch=19, cex=1.5)
+ }
+ }
+}
+
+require(animation)
+## GIF animation, loop=0 for infinite loop
+saveMovie(bifur(), interval=.2, loop=0, outdir=getwd())
+
+## shockwave flash animation:
+saveSWF(bifur(), interval=.2, dev="png", swfname="bifur.swf", height=300, width=300)
+
+lambda <- function(x, e, params){ }
+
+
+crit_pts <- sapply(A, function(a){ xhat(a, params) })
+good_xhat <-sapply(1:length(crit_pts), function(i) max(crit_pts[[i]]) )
+alpha <- sapply(1:length(good_xhat), function(i) lambda(good_xhat[i], A[i], params))
+sigma_squared <- sapply(1:length(good_xhat), function(i) b(good_xhat[i], params)+d(good_xhat[i], A[i], params))
+
+png("model1.png", 1000, 1000)
+par(mfrow=c(2,2))
+plot(A, good_xhat, type="l", lwd=4,xlab="bifurcation parameter a", ylab="x hat", cex.lab=2, cex.axis=2 )
+plot(A, alpha, type="l", lwd=4, xlab="bifurcation parameter a", cex.lab=2, cex.axis=2 )
+plot(A, sigma_squared, type="l", lwd=4, xlab="bifurcation parameter a", cex.lab=2, cex.axis=2 )
+plot(A, -sigma_squared/(2*alpha), type="l", lwd=4, xlab="bifurcation parameter a", cex.lab=2, cex.axis=2 )
+dev.off()
+
+
+
+## "Constant Mean" Dynamics
+d <- function(x, a, params) { a*(x-params["K"])+params["e"]*params["K"] }
+A<- seq(1,0,length=100)
+saveSWF(bifur(), interval=.2, dev="png", swfname="bifur_mean_const.swf", height=300, width=300)
+
+
+lambda <- function(x, a, params){ params["p"]*params["e"]*params["K"]*x^(params["p"]-1)/(x^params["p"]+params["h"]^params["p"]) - params["p"]*params["e"]*params["K"]*x^(2*params["p"]-1)/(x^params["p"]+params["h"]^params["p"])^2-a }
+
+crit_pts <- sapply(A, function(a){ xhat(a, params) })
+good_xhat <-sapply(1:length(crit_pts), function(i) max(crit_pts[[i]]) )
+alpha <- sapply(1:length(good_xhat), function(i) lambda(good_xhat[i], A[i], params))
+sigma_squared <- sapply(1:length(good_xhat), function(i) b(good_xhat[i], params)+d(good_xhat[i], A[i], params))
+
+
+
+png("model2.png", 1000, 1000)
+par(mfrow=c(2,2))
+plot(A, good_xhat, type="l", lwd=4,xlab="bifurcation parameter a", ylab="x hat", cex.lab=2, cex.axis=2 )
+plot(A, alpha, type="l", lwd=4, xlab="bifurcation parameter a", cex.lab=2, cex.axis=2 )
+plot(A, sigma_squared, type="l", lwd=4, xlab="bifurcation parameter a", cex.lab=2, cex.axis=2 )
+plot(A, -sigma_squared/(2*alpha), type="l", lwd=4, xlab="bifurcation parameter a", cex.lab=2, cex.axis=2 )
+dev.off()
+
+
+
+
+## Note that updating function definition updates prior functions, so the above doesn't need to redefine bifur():
+f <- function(x){ x }
+g <- function(x) { 2*f(x) }
+f <- function(x) {2*x }
+g(2)
+
+
+
+
+
+
+## Uniroot.all makes the previous code here unnecssary
File renamed without changes.

0 comments on commit 45bb7cd

Please sign in to comment.