[R] quasi-random sequences

[baptiste Auguié]

Hi again,

I've had a go at Prof Ripley's suggestion (Strauss process, code  
below). It works great for my limited purpose (qualitative drawing,  
really), I can just add a few mild concerns,

- ideally the hard core would not be a fixed number, but a  
distribution of sizes (ellipses).

- I could not quite understand the link between the window size,  
intensity, and number of elements returned. Is there a simple  
relation I've missed?

- in general, I just have no clue how rStrauss works, I guess Prof  
Ripley's first reference cited in ?Strauss would be useful for that  
matter as I could not find anything with google but a C code looking  
a bit alike some genetic algorithm (death and birth of objects).

Am I on the right track in thinking that the idea is not so  
dissimilar with the following approach:

1) create a first random distribution (Poisson, apparently)

2) give birth and death to new objects

3) carry on with this evolution until the sample distribution  
respects the given intensity and interaction parameters (e.g, for  
hard core, no overlap is permitted within the radius)

This looks a bit like the idea of an "electrostatic problem",

1) start with an initial configuration

2) define the total potential energy as a cost function to be  
minimized, sum of a boundary (repulsive term), and short distance  
repulsion

3) optimize the positions with respect to this objective function

I tried unsuccessfully this second approach (see below), but I'm sort  
of concerned about the handling of all these parameters by optim().  
rStrauss works beautifully anyway !

Thanks again,


baptiste

##############
#  Code
##############

library(spatstat) # Strauss process


## ellipse function obtained from R mailing list ##
## could use library(car) as an alternative
ellipse <- function(x,y,width,height,theta, npoints=100,plot=F)
                     {
   # x = x coordinate of center
   # y = y coordinate of center
   # width = length of major axis
   # height = length of minor axis
   # theta = rotation
   # npoints = number of points to send to polygon
   # plot = if TRUE, add to current device
   # = if FALSE, return list of components

   a <- width/2
   b <- height/2
   theta<-theta*pi/180
   xcoord <- seq(-a,a,length=npoints)
   ycoord.neg <- sqrt(b^2*(1-(xcoord)^2/a^2))
   ycoord.pos <- -sqrt(b^2*(1-(xcoord)^2/a^2))
   xx <- c(xcoord,xcoord[npoints:1])
   yy <- c(ycoord.neg,ycoord.pos)
   x.theta <- xx*cos(2*pi-theta)+yy*sin(2*pi-theta)+x
   y.theta <- yy*cos(2*pi-theta)-xx*sin(2*pi-theta)+y
   if(plot)
     invisible(polygon(x.theta,y.theta,col="black"))
   else
     invisible(list(coords=data.frame(x=x.theta,y=y.theta),
                    center=c(x,y),
                    theta=theta))
}


getEllipse<-function(dataf,plot=TRUE){
	ellipse(dataf[1],dataf[2],dataf[3],dataf[4],dataf[5],plot=plot)
	}


N <- 200 # initial number (loosely related to the actual number of  
elements)

##############################################
## positions from Strauss hard core process ##
##############################################
positions <- rStrauss(beta = 0.1, gamma = 0, R = 2, W= square(N/2))
# beta: intensity
# gamma: interaction.  0 : hard core (no overlap),  1: complete  
randomness
# R: radius of interaction (size of the core)
# W: window

N2<-length(positions$x)# number of elements

rand.angles <- rnorm(N2,sd=45,mean=45) # ellipse angles
rand.a <- rnorm(N2,sd=0.4,mean=1) # ellipse semi-axes
rand.b <- rnorm(N2,sd=0.4,mean=1) #

par(bty="n")
plot(cbind(positions$x,positions$y) ,cex = rand.a,axes=F,
	xlab="",ylab="",t="n") # just circles: type="p"

ell.pos<-cbind(positions$x,positions$y,rand.a,rand.b,rand.angles)
apply(ell.pos,1,getEllipse) -> b.quiet # draw ellipses

##############################################
##############################################


##############################################
## optimizing an "electrostatic potential" problem (not working)
##############################################

N <-200
p0 <- matrix(rnorm(2*N),ncol=2) # random starting point

# boundary repulsion potential #
x <- seq(-1.2,1.2,l=N)
delta <- 0.1
# 1D example #
y <- -1*( 1 - exp(-abs(x)/delta) - exp(abs(x)/delta))
#plot(x,y,t="l")

# 2D example
xy <- expand.grid(x = x, y=x)
z.x <- -1*( 1 - exp(-abs(xy[,1])/delta) - exp(abs(xy[,1])/delta))
z.y <- -1*( 1 - exp(-abs(xy[,2])/delta) - exp(abs(xy[,2])/delta))

z <- matrix(z.x+z.y, ncol=length(x))
#image(x=x,y=x,z=z)


boundary <- function(xy = p0[1,]){
	z.x <- -1*( 1 - exp(-abs(xy[1])/delta) - exp(abs(xy[1])/delta))
	z.y <- -1*( 1 - exp(-abs(xy[2])/delta) - exp(abs(xy[2])/delta))
	z.x + z.y
	}
	
dist.2d <- function(x = c(1, 0 ),y = c(0,1) , w = 0.1){
	
	(drop((x - y) %*% (x - y)) + w^2 )^(-3/2)
	
	}
	

obj <- function(p = p0){
	
	bound <- sum(apply(p,1,boundary)) # boundary term
	rel.positions <- expand.grid(x=p,y=p) # all relative positions
	bound + sum(apply(rel.positions,1,dist.2d)) # cost
	}
	
optim(p0,obj)	# fails miserably, but this does not sound right anyway

_____________________________

Baptiste Auguié

Physics Department
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK

Phone: +44 1392 264187

http://newton.ex.ac.uk/research/emag
http://projects.ex.ac.uk/atto