library(stats)

# input x = (x_0, x_1, ... x_{n+1}), and simulate the process on this,
# begun at w, up to time T, 
# killed upon hitting x_0 or x_{n+1}.
simul <- function(w, x, T=1000) {
	x <- sort(x);
	n <- length(x);
	dx <- diff(x);
	dir <- runif(T);
	loc <- c(which(w==x),rep(NA,T));
	if(is.na(loc[1])) {error("initial point not in state space");}
	i<-1;
	while (i<T && loc[i] != 0 && loc[i] != n) {
		loc[i+1] <- loc[i]+2*( (dir[i] < (dx[loc[i]-1]/(dx[loc[i]]+dx[loc[i]-1])) ) - 0.5 );
		i <- i+1;
		if (is.na(loc[i])) {browser()}
	}
	loc <- loc[1:i]; 
	dt <- rexp(i-1,1/(dx[loc[1:(i-1)]]*dx[loc[1:(i-1)]-1]));
	return(list(x=x, loc=loc, dt=dt));
}

# use this to plot output of simul.
plotsim <- function(x,loc=NA,dt=NA,...) {
	if (is.na(dt)) { dt <- x$dt; loc <- x$loc; x <- x$x; } 
	plot(c(0,0,sum(dt),sum(dt)),c(min(x[loc]),max(x[loc]),min(x[loc]),max(x[loc])),type="n",xlab="",ylab="",xaxt="n",yaxt="n");
	lines(rep(diffinv(dt),each=2)[-1],rep(x[loc],each=2)[-(2*length(loc))],lty=rep(c(0,1),length(loc))[-1],...);
	points(rep(0,length(x)),x,pch=20,cex=0.5);
}

#returns (m,n)-th approx to q-timescale
tq <- function(q,m=-15,n=5){
	xq <- q^(m:n);
	return(sort(c(xq,-xq,0)));
}

#returns n-th approx to cantor set
cantor <- function(n) {
	ccc <- c(0,2);
	for(i in 1:n) {
		ccc <- rep(ccc,each=2) + (-1)^(rep(c(0,0,1,1),length(ccc)/2))*rep(c(0,2),2)/3^i;
		ccc <- sort(ccc);
	}
	return(ccc);
}