# Example code to explore (re)construction of point patterns using simulated annealing. # Ap - Cuxhaven, 04.01.2013. Modified to accelerate computations by using C++ on 17.01.2013. # Functions rm(list = ls()) set.seed(-1) library(Rcpp) sourceCpp("/Users/afs80b/Dropbox/Rccp/numericUtilities.cpp") # sourceCpp("/Users/arnepommerening/Dropbox/Rccp/numericUtilities.cpp") # Initialisation density <- 0.025 xmax <- 100 ymax <- 100 nPoints <- rpois(1, density * xmax * ymax) x <- runif(nPoints, min = 0, max = xmax) y <- runif(nPoints, min = 0, max = ymax) Number <- seq(1, nPoints, by = 1) oldData <- data.frame(cbind(Number, x, y)) CEtarget <- 1.60 length(oldData$x) # Simulated annealing CE <- calcCE(xmax, ymax, oldData$x, oldData$y) E0 <- calcEnergy(CE, CEtarget) energyAim <- 5E-15 maxSimSteps <- 200000 simData <- oldData coolingFactor <- 0.9 Te <- 0 #0.00005 j <- 1 while (E0 > energyAim && j < maxSimSteps + 1) { xn <- round(runif(1, min = 1, max = length(simData$Number)), 0) tz <- round(runif(1, min = 1, max = 3), 0) if(tz == 1) { cofioX <- simData$x[xn] cofioY <- simData$y[xn] simData <- simData[-xn,] simData$Number <- seq(1, length(simData$Number), by = 1) } else if(tz == 2) { Number <- length(simData$Number) + 1 x <- runif(1, min = 0, max = xmax) y <- runif(1, min = 0, max = ymax) simData <- rbind(simData, cbind(Number, x, y)) } else { cofioX <- simData$x[xn] cofioY <- simData$y[xn] simData$x[xn] <- runif(1, min = 0, max = xmax) simData$y[xn] <- runif(1, min = 0, max = xmax) } CE <- calcCE(xmax, ymax, simData$x, simData$y) E1 <- calcEnergy(CE, CEtarget) Accepted <- T if(E1 > E0) { if(Te > 0) { u <- runif(1, min = 0, max = 1) p <- exp((E0 - E1) / Te) if(u < p) # Accept E0 <- E1 else { # Reject Accepted <- F if(tz == 1) { x <- cofioX y <- cofioY Number <- xn simData <- rbind(simData, cbind(Number, x, y)) } else if(tz == 2) { simData <- simData[-length(simData$Number),] } else { simData$x[xn] <- cofioX simData$y[xn] <- cofioY } } } else { # If T == 0 accept only improvements, i.e. reject in this case Accepted <- F if(tz == 1) { x <- cofioX y <- cofioY Number <- xn simData <- rbind(simData, cbind(Number, x, y)) } else if(tz == 2) { simData <- simData[-length(simData$Number),] } else { simData$x[xn] <- cofioX simData$y[xn] <- cofioY } } } else # If E1 < E0 always accept E0 <- E1 Te <- Te * coolingFactor cat("Iteration: ", j, "Point #: ", xn, " E0: ", E0, " E1: ", E1, " CE: ", CE, "T: ", Te, "tz: ", tz, "Number of points: ", length(simData$x), "Accepted: ", Accepted, "\n") j <- j + 1 } library(spatstat) oldDataP <- ppp(oldData$x, oldData$y, window = owin(c(0, xmax), c(0, ymax)), unitname = c("metres", "metres")) simDataP <- ppp(simData$x, simData$y, window = owin(c(0, xmax), c(0, ymax)), unitname = c("metres", "metres")) par(mar = c(0, 0, 0, 0), mfrow=c(1, 2)) plot(oldDataP, main = "", cex = 0.8, bg = "black", fg = "black") plot(simDataP, main = "", cex = 0.8, bg = "black", fg = "black")