I currently have a NextJS project using p5js deployed on https://artifical-life.vercel.app/smoothlife/smoothLife.
It is an implementation of https://arxiv.org/pdf/1111.1567.pdf which is essentially a cellular automata generalized on a continuous domain. It currently runs at around 10 to 14 frames per second and I want to increase this.
Here's the source file with the logic: https://github.com/JaylenLuc/ArtificalLife/blob/main/artificial_life/pages/smoothlife/smoothLife.tsx
export default function P5Sketch () {
/********************************************************
* UNIVERSAL LIFE CONSTANTS AND/OR VARS
********************************************************/
const renderRef = useRef(null);
var WIDTH_HEIGHT = 125 //the true number of cells WIDTH_HEIGHT ^ 2
//const HEIGHT = 150
var SIZE = 4
const RGB_MIN_RANGE = 255 //min range
const [strokePolicy, setStrokePolicy] = useState(false)
const [initOption, setInitPolicy] = useState("center")
const [seedUser, _setSeed] = useState(0)
const [colorScheme, _setColorScheme] = useState(0)
const setColorScheme = () => {
_setColorScheme((colorScheme + 1) % 3)
//console.log(colorScheme)
}
/****
* radius checks
* ****/
const ra_DEFAULT = 11
// const ra_DEFAULT = 11
var noLoop = false
const setnoLoop = (loopVal : boolean = !noLoop) => {
noLoop = loopVal
}
const [ra, _setOuterRadius] = useState(ra_DEFAULT)
const [ri, _setInnerRadius] = useState(ra_DEFAULT/3)
var ri_area = Math.PI * (ri*ri)
var ra_area = (Math.PI * (ra*ra)) - (ri_area)
/****
* sigmoid alpha values
* ****/
const alpha_n_DEFAULT = 0.028
const alpha_m_DEFAULT = 0.147
//αn = 0.028, αm = 0.147
const [alpha_n, setAlphaN] = useState(alpha_n_DEFAULT)
const [alpha_m, setAlphaM] = useState(alpha_m_DEFAULT)
/****
* birth and death interval values given by [b1, b2] and [d1, d2]
* ****/
const d1_DEFAULT = 0.267
const d2_DEFAULT = 0.445
const b1_DEFAULT = 0.278
const b2_DEFAULT = 0.365
// const d1_DEFAULT = 0.365
// const d2_DEFAULT = 0.549
// const b1_DEFAULT = 0.257
// const b2_DEFAULT = 0.336
//var d1 = d1_DEFAULT
const [d1, _setd1] = useState(d1_DEFAULT)
const [d2, _setd2] = useState(d2_DEFAULT)
const [b1, _setb1] = useState(b1_DEFAULT)
const [b2, _setb2] = useState(b2_DEFAULT)
/****
* delta time/ time stepping
* ****/
const dt_DEFAULT = 0.25 //time step
const [dt, setDeltaT] = useState(dt_DEFAULT)
const setDefaultParams = () => {
_setd1(d1_DEFAULT)
_setd2(d2_DEFAULT)
_setb1(b1_DEFAULT)
_setb2(b2_DEFAULT)
setRadius(ra_DEFAULT)
}
var cellsArray: number[] = []
const push_cellsAray = (row : number, col : number, val : number ) => {
cellsArray[WIDTH_HEIGHT * row + col] = val;
}
const get_cellsAray = (row : number, col : number ) => {
return cellsArray[WIDTH_HEIGHT * row + col];
}
//randomize the grid
//determine from the randomziex numbers what we render
//One PRESET
// var d1 = 0.365
// var d2 = 0.549
// var b1 = 0.257
// var b2 = 0.336
// var ra = 11 //outer radius
// var ri = ra/3 // inner radius
// var ri_area = Math.PI * (ri*ri)
// var ra_area = (Math.PI * (ra*ra)) - (ri_area)
// const ra_DEFAULT = 11
// var alpha_n = 0.028
// var alpha_m = 0.147
// var dt = 0.05 //time step
/********************************************************
* Event handlers
********************************************************/
const setSeed = (seed : number) => {
_setSeed(seed);
resetGrid();
};
const setRadius = (e : number) =>{
if (e <= 25){
_setOuterRadius(e)
_setInnerRadius(e/3)
}
}
/********************************************************
* GRID FUNCTIONS
********************************************************/
const random_number = (row: number, col : number, seed: number = seedUser) => {
//console.log("random_number func : ", seed)
if (seed > 0){
let random = Math.sin(seed + row * col) * 10000;
// console.log(random - Math.floor(random))
return random - Math.floor(random);
}else{
//console.log("default")
return Math.random();
}
}
const randomizeFullGrid = () => {
// for (let row = 0 ; row < WIDTH_HEIGHT; row ++){
// cellsArray.push([])
// for (let col = 0; col< WIDTH_HEIGHT; col ++){
// cellsArray[row].push( Math.random());
// }
// }
}
const initGridZero = () => {
//this function will likely be called first
for (let row = 0 ; row < WIDTH_HEIGHT; row ++){
for (let col = 0; col< WIDTH_HEIGHT; col ++){
push_cellsAray(row, col, 0)
}
}
}
const randomizeCenterGrid = (centerWidth : number) => {
//setStrokePolicy(false)
let center_grid = (Math.floor(WIDTH_HEIGHT/2))
let center_diff = (Math.floor((WIDTH_HEIGHT * centerWidth)/2))
let center_start = center_grid - center_diff
let center_end = center_grid + center_diff
for (let row = 0 ; row < WIDTH_HEIGHT; row ++){
for (let col = 0; col< WIDTH_HEIGHT; col ++){
if ((col > center_start && col <= center_end) &&
(row > center_start && row <= center_end)) push_cellsAray(row, col, random_number(row,col));
else push_cellsAray(row, col, seedUser != 0? random_number(row,col)/10 : 0);
}
}
}
const resetGrid = () => {
// console.log("in resetGrid: ", ra)
// console.log("in resetGrid: ", ri)
// console.log("in resetGrid: ", ri_area)
if (cellsArray.length > 0){
cellsArray = []
}
setnoLoop(false)
arbitrateMode()
}
const fillGrid = (p : any, myShader : any) => {
//it takes the rand nums in the grid and draws the color based on the numbers in the grid
let xPos = 0;
let yPos = 0;
for (let row = 0; row < WIDTH_HEIGHT; row ++){
xPos += SIZE
for (let col = 0; col< WIDTH_HEIGHT; col ++){
yPos += SIZE
let current_state = get_cellsAray(row, col)
let fill_value = (current_state * RGB_MIN_RANGE);
/*
fill value is the rgb 255 * the decimal which
is a percentage of how far it is from black until white
this can be changed with the fill value determining
only some of the RGB values , play around
*/
//myShader.setUniform('myColor', [fill_value,fill_value,fill_value])
p.fill(fill_value,fill_value,fill_value)
p.circle(yPos, xPos , SIZE);
// p.fill(70,70,70)
// p.circle(yPos, xPos , 5);
}
yPos = 0
}
}
const clamp = function(value : number, lower_b : number, upper_b : number) {
return Math.min(Math.max(value, lower_b), upper_b);
};
const clamp_test = function(value : number, lower_b : number, upper_b : number) {
if (value < lower_b ) return lower_b
else if (value > upper_b ) return upper_b
else return value
};
const generalizeTransitionFunc = ( ) => {
let row = 0
let col = 0
while ( row < WIDTH_HEIGHT){
while(col < WIDTH_HEIGHT){
let m_n : Array<number> = fillingIntegralN_M(row, col)
let new_value = dt * transitionFunc_S(m_n[1], m_n[0]) // [-1,1]
//console.log(new_value)
//smooth time stepping scheme
//console.log("not clamped: ", new_value)
//f(~x, t + dt) = f(~x, t) + dt S[s(n, m)] f(~x, t)
push_cellsAray(row, col, clamp_test(get_cellsAray(row, col) + new_value, 0, 1 ))
col++
}
col = 0
row += 1
}
}
/********************************************************
* INTEGRAL FUNCTIONS
********************************************************/
const emod = (pos : number, size : number ) => {
return ((pos % size) + size ) % size
}
//outer and inner filling
const fillingIntegralN_M = (_row : number = 0, _col : number = 0) => { //return value between 0 -1 normalize
let c_row : number = _row
let c_col : number = _col
let m : number = 0
let n : number = 0
for (let d_row = -(ra - 1) ; d_row < (ra - 1); ++d_row){ //iterate over the outer radius
let real_pos_row = emod(d_row + c_row, WIDTH_HEIGHT)
for (let d_col = -(ra -1); d_col < (ra -1); ++ d_col){
let real_pos_col = emod(d_col + c_col, WIDTH_HEIGHT)
if (d_row*d_row + d_col* d_col <= ri*ri){ //inner
m += get_cellsAray(real_pos_row,real_pos_col)
//M ++
}else if (d_row*d_row + d_col* d_col <= ra*ra) {//outer
n += get_cellsAray(real_pos_row,real_pos_col)
//N ++
}
}
}
m /= ri_area
//m= clamp(m , 0 ,1 )
n /= ra_area
//n = clamp(n , 0 ,1)
return [m,n] // inner, outer
}
/********************************************************
* SIGMOIDS AND THE TRANSITION FUNCTION
********************************************************/
const sigmoid1 = (x : number , a : number, alpha_val : number) => {
return 1/(1 + Math.exp(-(x-a) * 4/alpha_val))
}
const sigmoid2 = (x : number , a : number , b : number) => {
return sigmoid1(x,a, alpha_n) * (1- sigmoid1(x,b, alpha_n))
}
const sigmoidM = (x :number , y : number, m : number) => {
return x * ( 1 - sigmoid1(m, 0.5,alpha_m) ) + (y * sigmoid1(m, 0.5,alpha_m))
}
const transitionFunc_S = (n : number, m : number ) => {
//If the transition function in the discrete time-stepping scheme was sd(n, m) then the smooth one is s(n, m) = 2sd(n, m)−1.
return 2 * sigmoid2(n, sigmoidM(b1,d1,m), sigmoidM(b2,d2,m) ) - 1
}
const arbitrateMode = () => {
ri_area = (Math.PI * (ri*ri))
ra_area = ((Math.PI * (ra*ra)) - (ri_area))
switch(initOption){
case "full":
randomizeFullGrid();
break;
case "center":
randomizeCenterGrid(0.35);
break;
}
}
/******************************************************
********************************************************/
// const vertex = `
// precision highp float;
// uniform mat4 uModelViewMatrix;
// uniform mat4 uProjectionMatrix;
// attribute vec3 aPosition;
// attribute vec2 aTexCoord;
// varying vec2 vTexCoord;
// void main() {
// vTexCoord = aTexCoord;
// vec4 positionVec4 = vec4(aPosition, 1.0);
// gl_Position = uProjectionMatrix * uModelViewMatrix * positionVec4;
// }
// `;
// // the fragment shader is called for each pixel
// const fragment = `
// precision highp float;
// uniform vec2 p;
// uniform float r;
// const int I = 500;
// varying vec2 vTexCoord;
// void main() {
// vec2 c = p + gl_FragCoord.xy * r, z = c;
// float n = 0.0;
// for (int i = I; i > 0; i --) {
// if(z.x*z.x+z.y*z.y > 4.0) {
// n = float(i)/float(I);
// break;
// }
// z = vec2(z.x*z.x-z.y*z.y, 2.0*z.x*z.y) + c;
// }
// gl_FragColor = vec4(0.5-cos(n*17.0)/2.0,0.5-cos(n*13.0)/2.0,0.5-cos(n*23.0)/2.0,1.0);
// }`;
useEffect(() => {
const p5 = require("p5");
p5.disableFriendlyErrors = true;
var myShader: any;
var fps ;
const p5instance = new p5((p : any) => {
// p.preload = () => {
// // load each shader file (don't worry, we will come back to these!)
// myShader = p.createShader(vertex, fragment);
// }
p.setup = () => {
//console.log("HERE:", p5.disableFriendlyErrors )
p.createCanvas( WIDTH_HEIGHT * SIZE, WIDTH_HEIGHT * SIZE).parent(renderRef.current);
//p.createGraphics( WIDTH_HEIGHT + 200,WIDTH_HEIGHT + 200)
if (!strokePolicy) p.noStroke()
let color = ""
switch(colorScheme){
case 0 : p.colorMode(p.RGB); break;
case 1 : p.colorMode(p.HSB); break;
case 2: p.colorMode(p.HSL); break;
default: p.colorMode(p.RGB);
}
arbitrateMode()
}
p.draw = () => {
//fps = p.frameRate();
p.background(0,0,0);
//p.shader(myShader)
// console.time('generalizeTransitionFunc')
if (!noLoop){
generalizeTransitionFunc()
// console.timeEnd('generalizeTransitionFunc')
//console.log( "loooping", cellsArray[13][55])
}
fillGrid(p, myShader);
console.log(p.frameRate());
}
})
return () => {
//console.log("cleaning up...");
// comment this out to get 2 canvases and 2 draw() loops
p5instance.remove();
};
}, [ strokePolicy, seedUser, initOption, b1, b2, d1, d2, dt, ra, ri, ri_area, ra_area, noLoop, colorScheme])
//the entropy of the universe is tending to a maximum
return(
<div className={styles.foreground}>
<Sparkles
color="random"
count={80}
minSize={9}
maxSize={14}
overflowPx={80}
fadeOutSpeed={30}
flicker={true}
/>
<meta name="viewport" content="width=device-height"></meta>
<div className={styles.title}>
The Universe moves to an Entropic Maximum
</div>
<section>
<div className= {styles.life_box} ref={renderRef}></div>
</section>
<div className={styles.buttonlayout}>
<ButtonLayout setStrokePolicy = {setStrokePolicy} strokePolicy = {strokePolicy}/>
<BigBangButton resetHandler={resetGrid}></BigBangButton>
<SeedInput setSeed={setSeed} seedUser = {seedUser}/>
</div>
<br></br>
<div className={styles.buttonlayout}>
<ColorButton changeColor = {setColorScheme}/>
<FourDButton setNoLoop = {setnoLoop} noLoop = {noLoop}/>
</div>
<div className={styles.buttonlayout}>
<ParamNav setd1 = {_setd1} d1 = {d1} setd2= {_setd2}
d2 = {d2} setb1={_setb1} b1={b1} setb2={_setb2} b2={b2} setrad={setRadius} rad = {ra} resetSettings={setDefaultParams}/>
</div>
</div>
)
}
Any ideas on how to increase framerate?
Here's the high level implementation:
It currently has to initialize a 1D array with numbers from 0 to 1:
const push_cellsAray = (row : number, col : number, val : number ) => {
cellsArray[WIDTH_HEIGHT * row + col] = val;
}
const randomizeCenterGrid = (centerWidth : number) => {
//setStrokePolicy(false)
let center_grid = (Math.floor(WIDTH_HEIGHT/2))
let center_diff = (Math.floor((WIDTH_HEIGHT * centerWidth)/2))
let center_start = center_grid - center_diff
let center_end = center_grid + center_diff
for (let row = 0 ; row < WIDTH_HEIGHT; row ++){
for (let col = 0; col< WIDTH_HEIGHT; col ++){
if ((col > center_start && col <= center_end) &&
(row > center_start && row <= center_end)) push_cellsAray(row, col, random_number(row,col));
else push_cellsAray(row, col, seedUser != 0? random_number(row,col)/10 : 0);
}
}
}
loop through the array again, check all its neighbors in the given inner radius ri and outer radius ra and apply a transition function to every value of the array after getting values M and N from the fillingIntegralN_M function to get a new value for that given row and column:
const clamp_test = function(value : number, lower_b : number, upper_b : number) {
if (value < lower_b ) return lower_b
else if (value > upper_b ) return upper_b
else return value
};
const sigmoid1 = (x : number , a : number, alpha_val : number) => {
return 1/(1 + Math.exp(-(x-a) * 4/alpha_val))
}
const sigmoid2 = (x : number , a : number , b : number) => {
return sigmoid1(x,a, alpha_n) * (1- sigmoid1(x,b, alpha_n))
}
const sigmoidM = (x :number , y : number, m : number) => {
return x * ( 1 - sigmoid1(m, 0.5,alpha_m) ) + (y * sigmoid1(m, 0.5,alpha_m))
}
const transitionFunc_S = (n : number, m : number ) => {
//If the transition function in the discrete time-stepping scheme was sd(n, m) then the smooth one is s(n, m) = 2sd(n, m)−1.
return 2 * sigmoid2(n, sigmoidM(b1,d1,m), sigmoidM(b2,d2,m) ) - 1
}
const fillingIntegralN_M = (_row : number = 0, _col : number = 0) => { //return value between 0 -1 normalize
let c_row : number = _row
let c_col : number = _col
let m : number = 0
let n : number = 0
for (let d_row = -(ra - 1) ; d_row < (ra - 1); ++d_row){ //iterate over the outer radius
let real_pos_row = emod(d_row + c_row, WIDTH_HEIGHT)
for (let d_col = -(ra -1); d_col < (ra -1); ++ d_col){
let real_pos_col = emod(d_col + c_col, WIDTH_HEIGHT)
if (d_row*d_row + d_col* d_col <= ri*ri){ //inner
m += get_cellsAray(real_pos_row,real_pos_col)
//M ++
}else if (d_row*d_row + d_col* d_col <= ra*ra) {//outer
n += get_cellsAray(real_pos_row,real_pos_col)
//N ++
}
}
}
m /= ri_area
//m= clamp(m , 0 ,1 )
n /= ra_area
//n = clamp(n , 0 ,1)
return [m,n] // inner, outer
}
const generalizeTransitionFunc = ( ) => {
let row = 0
let col = 0
while ( row < WIDTH_HEIGHT){
while(col < WIDTH_HEIGHT){
let m_n : Array<number> = fillingIntegralN_M(row, col)
let new_value = dt * transitionFunc_S(m_n[1], m_n[0]) // [-1,1]
//console.log(new_value)
//smooth time stepping scheme
//console.log("not clamped: ", new_value)
//f(~x, t + dt) = f(~x, t) + dt S[s(n, m)] f(~x, t)
push_cellsAray(row, col, clamp_test(get_cellsAray(row, col) + new_value, 0, 1 ))
col++
}
col = 0
row += 1
}
}
then fill the grid with a color value normalized from 0 to 255 based on those values:
const fillGrid = (p : any, myShader : any) => {
//it takes the rand nums in the grid and draws the color based on the numbers in the grid
let xPos = 0;
let yPos = 0;
for (let row = 0; row < WIDTH_HEIGHT; row ++){
xPos += SIZE
for (let col = 0; col< WIDTH_HEIGHT; col ++){
yPos += SIZE
let current_state = get_cellsAray(row, col)
let fill_value = (current_state * RGB_MIN_RANGE);
/*
fill value is the rgb 255 * the decimal which
is a percentage of how far it is from black until white
this can be changed with the fill value determining
only some of the RGB values , play around
*/
//myShader.setUniform('myColor', [fill_value,fill_value,fill_value])
p.fill(fill_value,fill_value,fill_value)
p.circle(yPos, xPos , SIZE);
// p.fill(70,70,70)
// p.circle(yPos, xPos , 5);
}
yPos = 0
}
}
It repeats this process indefinitely.
The main hit to performance comes from increasing the radius and increasing the board size.
Is it possible to improve performance without the use of shaders?
The main hit to performance comes from increasing the radius and increasing the board size.
I didn't look at your code in detail, but your fillingIntegralN_M function essentially computes a 2D convolution, doesn't it?
For large kernel sizes (and in practice even a radius of a few grid units might be "large" enough), it's more efficient to do that with a fast Fourier transform (FFT).
Basically, what you would do is precalculate the FFTs of your inner and outer "filling integral" kernels, i.e. of grids where the cells within the inner and outer radii of the origin have value 1 (or 1/M or 1/N, if you want to bake the normalization factors into the kernels as an optimization) and all cells outside the radius have value 0, with "anti-aliasing" applied for cells at the boundary as described in the paper.
(This assumes that your grid "wraps around" at the edges, so that the cell at (width-1, height-1) is diagonally adjacent to (0, 0). The Fourier transform naturally operates on this kind of cyclic domain. This is also what your emod function seems to implement, so you should be good here.)
On each iteration of your simulation loop, you then take the FFT of your simulation grid, multiply it pointwise with the FFTs of each kernel, and take the inverse FFTs of the resulting grids to obtain two grids that contain the m and n values for each point. Then you just use these m and n values to update the grid and repeat.
(As an aside, if your update rule was completely linear, you could do all the calculations in "Fourier space" and only occasionally take inverse FFTs to display the grid to the viewer. But since the SmoothLife update rule is explicitly non-linear, you do need to jump back and forth between "normal space", where non-linear maps are easy but convolutions are hard, and "Fourier space", where non-linear maps are hard but convolutions are easy. But luckily a well optimized FFT implementation makes the transition back and forth pretty quick.)
You probably don't want to write your own JavaScript FFT implementation, at least not at first — it'd be a whole separate project, and you'd have to spend an enormous amount of work to get it not only working correctly but as fast as existing optimized implementations are.
Fortunately, Googling for "javascript fft" turns up several results, including fft.js, jsfft and the math.fft function in math.js. Out of those, at least the last one explicitly claims to support 2-dimensional FFTs, but if you use a library that doesn't, you can always turn a 1D FFT into a 2D FFT by applying the 1D FFT to each row of the array, then transposing the transformed array (i.e. turning rows into columns and vice versa) and applying the 1D FFT to the rows again.
FWIW, apparently p5.js also has an FFT module, but it seems to be designed for audio processing applications and I'm not 100% sure if you can use it to get the "raw" FFT of arbitrary numerical data. (The FFT.analyze function might do that, but the documentation isn't very clear to me.) It also seems to have some limitations, such as a maximum input/output size of 1024 values, meaning that using it would limit you to a 1024 by 1024 point grid.
Addendum: Here's some sample code to demonstrate FFT convolution using the jsfft library and how to apply it to the SmoothLife algorithm. I wrote it mostly based on the paper you linked, so it might not exactly match your code, but hopefully it's close enough that you can adapt it.
import("https://cdn.jsdelivr.net/npm/jsfft@0.0.4/+esm").then(fft => {
// Grid size (powers of 2 are faster!)
const WIDTH = 512, HEIGHT = 256
// Simulation parameters
const RA = 20, RI = RA / 3 // outer and inner radius
const B1 = 0.257, B2 = 0.336, D1 = 0.365, D2 = 0.549, ALPHA_N = 0.028, ALPHA_M = 0.147
const DT = 0.05 // simulation time step
const INIT_SCALE = 0.512
// Create a w * h sized FFT convolution kernel based on a function that maps (x,y) pairs to real
// numbers. The FFT kernel is scaled so that the complex norm of its DC coefficient becomes 1.
function makeKernel(width, height, realFunction) {
const kernel = new fft.ComplexArray(WIDTH * HEIGHT)
for (let y = 0; y < height; y++) {
for (let x = 0; x < width; x++) {
kernel.real[x + width * y] = realFunction(x, y)
}
}
const fftKernel = kernel.FFT()
const scale = 1 / Math.sqrt(fftKernel.real[0]**2 + fftKernel.imag[0]**2)
for (let i = 0; i < kernel.length; i++) {
fftKernel.real[i] *= scale
fftKernel.imag[i] *= scale
}
return fftKernel;
}
// Helper function to clamp a real number to the interval [0, 1]:
function clampTo01(realNumber) {
return Math.max(0, Math.min(1, realNumber))
}
// Predefine two nornalized convolution kernels: one a filled circle of radius RI around the origin,
// the other an annulus of inner radius RI and outer radius RA.
const mKernel = makeKernel(WIDTH, HEIGHT, (x, y) => {
const dx = Math.min(x, WIDTH - x), dy = Math.min(y, HEIGHT - y), dist = Math.sqrt(dx**2 + dy**2)
return clampTo01(RI + 0.5 - dist)
})
const nKernel = makeKernel(WIDTH, HEIGHT, (x, y) => {
const dx = Math.min(x, WIDTH - x), dy = Math.min(y, HEIGHT - y), dist = Math.sqrt(dx**2 + dy**2)
return clampTo01(RA + 0.5 - dist) * (1 - clampTo01(RI + 0.5 - dist))
})
// Convolve a Fourier-transformed grid with a precomputed kernel (e.g. as returned by makeKernel).
function fftConvolve(freqs, kernel) {
const output = new fft.ComplexArray(freqs.length)
for (let i = 0; i < freqs.length; i++) {
// Complex number multiplication: (a + bi)*(c + di) = (ac - bd) + (ad + bc)i
output.real[i] = freqs.real[i] * kernel.real[i] - freqs.imag[i] * kernel.imag[i]
output.imag[i] = freqs.real[i] * kernel.imag[i] + freqs.imag[i] * kernel.real[i]
}
return output.InvFFT().real
}
// Apply the SmoothLife rule to a grid using the given nonlinear s(n, m) function.
function smoothLifeStep(grid, dt, sFunction) {
const freqs = new fft.ComplexArray(grid).FFT()
const m = fftConvolve(freqs, mKernel)
const n = fftConvolve(freqs, nKernel)
for (let i = 0; i < grid.length; i++) {
grid[i] = clampTo01(grid[i] + dt *(2 * sFunction(n[i], m[i]) - 1))
// Alternative (smoother) update rule:
// grid[i] = (1-dt) * grid[i] + dt * sFunction(n[i], m[i])
}
}
// Define the s(n, m) function.
const sigma1 = (x, a, alpha) => 1 / (1 + Math.exp(-4 * (x - a) / alpha))
const sigma2 = (x, a, b) => sigma1(x, a, ALPHA_N) * (1 - sigma1(x, b, ALPHA_N))
const sigmaM = (x, y, m) => x * (1 - sigma1(m, 0.5, ALPHA_M)) + y * sigma1(m, 0.5, ALPHA_M)
const s = (n, m) => sigma2(n, sigmaM(B1, D1, m), sigmaM(B2, D2, m))
// Set up a canvas for drawing.
const canvas = document.createElement("canvas")
canvas.width = WIDTH
canvas.height = HEIGHT
document.getElementById("smoothlife").appendChild(canvas)
// Draw grid to canvas.
function drawToCanvas(grid) {
const ctx = canvas.getContext("2d")
const img = ctx.getImageData(0, 0, WIDTH, HEIGHT)
const data = img.data
for (let i = 0; i < WIDTH * HEIGHT; i++) {
const j = 4*i, grayscale = 255 * grid[i]
data[j] = data[j+1] = data[j+2] = grayscale // RGB
data[j+3] = 255 // fully opaque
}
ctx.putImageData(img, 0, 0)
}
// Debug: uncomment to draw the kernels on the canvas.
//drawToCanvas(new fft.ComplexArray(mKernel).InvFFT().real); return
// Initialized a random grid.
const grid = new Float32Array(WIDTH * HEIGHT)
for (let i = 0; i < grid.length; i++) grid[i] = Math.random() * INIT_SCALE
drawToCanvas(grid)
// Run SmoothLife.
const fpsDiv = document.getElementById("fps")
let startTime = undefined
let simTime = 0
function run(now) {
if (startTime === undefined) startTime = now
smoothLifeStep(grid, DT, s)
drawToCanvas(grid)
simTime += DT
const seconds = (now - startTime) / 1000
const fps = (simTime / DT) / seconds
fpsDiv.textContent = "t = " + simTime.toFixed(3) + " (FPS: " + fps.toFixed(2) + ")"
window.requestAnimationFrame(run)
}
window.requestAnimationFrame(run)
})<div id="smoothlife"></div>
<div id="fps"></div>On my laptop, this example runs at about 20 FPS with the 512 by 256 grid size. It could probably be optimized to improve that. The grid size doesn't have to be a power of 2, but the FFT code seems to run noticeably faster when it is.
One detail worth noting is that I'm not actually using a 2D FFT, but just a simple 1D FFT. For convolution it turns out that it doesn't actually matter what kind of FFT you use, as long as you're using the same kind consistently. And since I'm already using one-dimensional arrays to store two-dimensional grids (just like your code does), just using a one-dimensional FFT simplifies the code and probably improves performance.
(Well, it almost doesn't matter. Using a 1D FFT like this instead of a proper 2D FFT does change the global connectivity of the grid so that e.g. the cell one step to the right of (WIDTH-1, 0) isn't (0, 0) but rather (0, 1). But that's a very subtle difference and basically impossible to notice from the animation.)
Also, I'm still having some trouble reproducing classic SmoothLife behavior like that seen in this video, even with the same parameter values, so my code might still have some subtle bugs. (Or it might differ from the reference implementation in some undocumented details; I took at quick peek at the code on SourceForge and it seems to have a bunch of options for things like different integration schemes and different kinds of sigmoid functions etc.)