brew install imagemagick –with-fontconfig –with-librsvg –with-fftw— title: “Lenia” output: html.notebook —
SIZE <- 2^8
MID <- SIZE / 2
kernel.core <- function(r, kernel.type)
{
rm <- pmin(r, 1)
if (kernel.type == 0)
return ( (4 * rm * (1-rm))^4 )
else
return ( exp(4 - 1 / (rm * (1-rm))) )
}
kernel.shell <- function(r, peaks, kernel.type)
{
k <- length(peaks)
kr <- k * r
peak = peaks[pmin(floor(kr), k-1)+1]
return ( (r<1) * kernel.core(kr %% 1, kernel.type) * peak )
}
delta.func <- function(n, mu, sigma, delta.type)
{
if (delta.type == 0)
return ( pmax(0, 1 - (n - mu)^2 / (sigma^2 * 9) )^4 * 2 - 1 )
else
return ( exp( - (n - mu)^2 / (sigma^2 * 2) ) * 2 - 1 )
}
R <- 52
peaks <- c(1/2, 2/3, 1)
mu <- 0.3
sigma <- 0.03
dt <- 0.1
kernel.type <- 0
delta.type <- 0
x <- seq(0, 1.2, 0.01)
plot(x, kernel.core(x, kernel.type), type="l", ylim=c(0, 1), col="cyan", xlab="r", ylab="kernel")
lines(x, kernel.shell(x, peaks, kernel.type), col="blue")
plot(x, delta.func(x, mu, sigma, delta.type), type="l", ylim=c(-1, 1), col="magenta", xlab="n", ylab="delta")
jet <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))(100)
show.world <- function(A, vmin=0, vmax=1, is.display=TRUE)
{
par(mar=rep(0, 4))
image(1:dim(A)[1], 1:dim(A)[1], A, #pmax(pmin(A,vmax),vmin),
zlim=c(vmin, vmax), asp=1,
col=jet, axes=FALSE, xlab="", ylab="")
}
fft.shift <- function(A)
{
r1 <- 1; r4 <- dim(A)[1]; r2 <- ceiling(r4/2) - 1; r3 <- r2 + 1
c1 <- 1; c4 <- dim(A)[2]; c2 <- ceiling(c4/2) - 1; c3 <- c2 + 1
A <- rbind(A[r3:r4, c1:c4], A[r1:r2, c1:c4])
A <- cbind(A[r1:r4, c3:c4], A[r1:r4, c1:c2])
return ( A )
}
calc.kernel <- function(R)
{
I <- matrix(1:SIZE, SIZE, SIZE)
X <- (I-MID) / R
Y <- t(X)
D <- sqrt(X^2 + Y^2)
kernel <- kernel.shell(D, peaks, kernel.type)
kernel.sum <- sum(kernel)
kernel.norm <- kernel / kernel.sum
kernel.FFT <- fft(kernel.norm)
return (list("kernel"=kernel, "kernel.FFT"=kernel.FFT, "D"=D))
}
output <- calc.kernel(R)
kernel <- output$kernel; kernel.FFT <- output$kernel.FFT; D <- output$D
#show.world(D, 0, 5)
show.world(kernel)
#show.world(fft.shift(abs(kernel.FFT)))
load.cells <- function(id)
{
if (id==0) {
name <<- 'Orbium bicaudatus';
R <<- 13; peaks <<- c(1); mu <<- 0.15; sigma <<- 0.014; dt <<- 0.1; cells <<- matrix(c(0,0,0,0,0,0,0.1,0.14,0.1,0,0,0.03,0.03,0,0,0.3,0,0,0,0, 0,0,0,0,0,0.08,0.24,0.3,0.3,0.18,0.14,0.15,0.16,0.15,0.09,0.2,0,0,0,0, 0,0,0,0,0,0.15,0.34,0.44,0.46,0.38,0.18,0.14,0.11,0.13,0.19,0.18,0.45,0,0,0, 0,0,0,0,0.06,0.13,0.39,0.5,0.5,0.37,0.06,0,0,0,0.02,0.16,0.68,0,0,0, 0,0,0,0.11,0.17,0.17,0.33,0.4,0.38,0.28,0.14,0,0,0,0,0,0.18,0.42,0,0, 0,0,0.09,0.18,0.13,0.06,0.08,0.26,0.32,0.32,0.27,0,0,0,0,0,0,0.82,0,0, 0.27,0,0.16,0.12,0,0,0,0.25,0.38,0.44,0.45,0.34,0,0,0,0,0,0.22,0.17,0, 0,0.07,0.2,0.02,0,0,0,0.31,0.48,0.57,0.6,0.57,0,0,0,0,0,0,0.49,0, 0,0.59,0.19,0,0,0,0,0.2,0.57,0.69,0.76,0.76,0.49,0,0,0,0,0,0.36,0, 0,0.58,0.19,0,0,0,0,0,0.67,0.83,0.9,0.92,0.87,0.12,0,0,0,0,0.22,0.07, 0,0,0.46,0,0,0,0,0,0.7,0.93,1,1,1,0.61,0,0,0,0,0.18,0.11, 0,0,0.82,0,0,0,0,0,0.47,1,1,0.98,1,0.96,0.27,0,0,0,0.19,0.1, 0,0,0.46,0,0,0,0,0,0.25,1,1,0.84,0.92,0.97,0.54,0.14,0.04,0.1,0.21,0.05, 0,0,0,0.4,0,0,0,0,0.09,0.8,1,0.82,0.8,0.85,0.63,0.31,0.18,0.19,0.2,0.01, 0,0,0,0.36,0.1,0,0,0,0.05,0.54,0.86,0.79,0.74,0.72,0.6,0.39,0.28,0.24,0.13,0, 0,0,0,0.01,0.3,0.07,0,0,0.08,0.36,0.64,0.7,0.64,0.6,0.51,0.39,0.29,0.19,0.04,0, 0,0,0,0,0.1,0.24,0.14,0.1,0.15,0.29,0.45,0.53,0.52,0.46,0.4,0.31,0.21,0.08,0,0, 0,0,0,0,0,0.08,0.21,0.21,0.22,0.29,0.36,0.39,0.37,0.33,0.26,0.18,0.09,0,0,0, 0,0,0,0,0,0,0.03,0.13,0.19,0.22,0.24,0.24,0.23,0.18,0.13,0.05,0,0,0,0, 0,0,0,0,0,0,0,0,0.02,0.06,0.08,0.09,0.07,0.05,0.01,0,0,0,0,0), ncol=20)
}
}
clear.world <- function()
{
world <<- matrix(0L, SIZE, SIZE)
}
random.world <- function()
{
world <<- matrix(runif(SIZE*SIZE) * 0.6, SIZE, SIZE)
}
add.cells <- function()
{
w <- dim(cells)[2]
h <- dim(cells)[1]
i <- MID - floor(w / 2)
j <- MID - floor(h / 2)
world[j:(j+h-1), i:(i+w-1)] <<- cells
}
multiply.cells <- function(n=2)
{
w <- dim(cells)[2] * n
h <- dim(cells)[1] * n
cells2 <- matrix(0L, h, w)
for (i in 1:n)
for (j in 1:n)
cells2[seq(i,h,n), seq(j,w,n)] <- cells
cells <<- cells2
R <<- R * n
}
ifft <- function(A) { fft(A, inverse=TRUE) / length(A) }
calc.once <- function()
{
world.FFT <<- fft(world)
potential <<- fft.shift(Re(ifft(kernel.FFT * world.FFT)))
delta <<- delta.func(potential, mu, sigma, delta.type)
delta <<- matrix(delta, SIZE, SIZE)
world <<- pmax(0, pmin(1, world + delta * dt))
world <<- matrix(world, SIZE, SIZE)
}
clear.world()
load.cells(0)
multiply.cells(4)
output <- calc.kernel(R)
kernel <- output$kernel; kernel.FFT <- output$kernel.FFT; D <- output$D
add.cells()
#show.world(world)
library(animation)
saveGIF({
for (gen in 1:50) {
calc.once()
show.world(world)
}
}, movie.name="orbium.gif", interval=0.1)
## Executing:
## convert -loop 0 -delay 10 Rplot1.png Rplot2.png Rplot3.png
## Rplot4.png Rplot5.png Rplot6.png Rplot7.png Rplot8.png
## Rplot9.png Rplot10.png Rplot11.png Rplot12.png Rplot13.png
## Rplot14.png Rplot15.png Rplot16.png Rplot17.png Rplot18.png
## Rplot19.png Rplot20.png Rplot21.png Rplot22.png Rplot23.png
## Rplot24.png Rplot25.png Rplot26.png Rplot27.png Rplot28.png
## Rplot29.png Rplot30.png Rplot31.png Rplot32.png Rplot33.png
## Rplot34.png Rplot35.png Rplot36.png Rplot37.png Rplot38.png
## Rplot39.png Rplot40.png Rplot41.png Rplot42.png Rplot43.png
## Rplot44.png Rplot45.png Rplot46.png Rplot47.png Rplot48.png
## Rplot49.png Rplot50.png 'orbium.gif'
## Output at: orbium.gif
## [1] TRUE
#print(dim(world))
#print(dim(ifft))
#print(dim(potential))
#print(dim(delta))
#show.world(world)
#show.world(Re(fft.shift(a)) / SIZE^2)
#show.world(potential)