Geneascopy

For scripts to aid with computation or simulation in cellular automata.
Post Reply
User avatar
Rhombic
Posts: 1057
Joined: June 1st, 2013, 5:41 pm

Geneascopy

Post by Rhombic » March 17th, 2019, 5:32 pm

[UPDATED SCRIPT AT THE BOTTOM OF THIS POST]

geneasco.py
I had been toying around with the idea, inspired by NMR, to probe "stuff" (I will expand on this) in CA in a systematic manner that can yield useful results. I was surprised there has not been much study in the dynamics of rules in general and CGoL in particular starting from random configurations, and I have devised a script that can allow both this and, in the future, probing specific pattern families.

What is geneasco.py?
It is a Python script that runs ns 16x16 soups in any rule for kmax generations, recording generation by generation their population AND their population change (the average is given by the results divided by ns, obviously).
The idea behind this is that random soups are not random in their evolution and, on top of significant background noise, can show surprisingly important trends in population and population change when averaged across enough examples.

What information this gives us: Population
I greatly suggest you run the script for CGoL with ns=2000 and kmax=2000 (should be fairly quick, less than half a minute) and plot the resulting cumulative population against generation in Excel. The resulting plot is what seems to be a simple Morse potential-like function. The parameters for that Morse potential-like function will be very interesting to understand, as it will give insight into the dynamics of the rule/soup size/soup density/etc. You will see surprisingly well-defined humps within that general function though, and I think we can safely call this fine-structure, which I think has vital information about the rule's dynamics (think X-ray crystallography in chemistry - the fine detail within the diffraction points gives info on the atom positions in the unit cell).
Unexpectedly, if you run it for 5000 generations especially for loads of scans (40000), you see the population falling again after a peak at about 2100. Interestingly, if you zoom in enough with these many scans, you see a VERY well-defined fine structure!*
finestructure.JPG
finestructure.JPG (102.19 KiB) Viewed 2535 times
Which brings us to this. That "fine-structure" I have found to differ extremely significantly between relatively closely-related rules. I think further studying this will be extremely interesting.

What information this gives us: Population Change
Population change converges apparently oscillating around 0 (actually, around the expected end heat of a soup).
Plotting ln(change^2)/2 vs generation shows something unique: after a given generation threshold, the behaviour shows a slow decrease with a lot of noise mostly under a maximum y=mx+c, related to the chaotic heat of the rule. However, before this generation threshold, in the case of CGoL being around 160 generations, the curve followed is not linear and very precise, behaving in an entirely different manner. The complexity of the system seems to be defined by generation 160-ish for CGoL, before which the evolution of most soups has strikingly similar features (in fact, the dip at about 28 ticks is actually where you are most likely to find B-heptominos, r-pentominos, pis etc!).
concertedevolution.JPG
concertedevolution.JPG (59.31 KiB) Viewed 2530 times
*NOTE: I upload the raw results from the 5000-generation 40k scans as a text file.
gscopy_ns_40k.txt
(89.09 KiB) Downloaded 122 times
Future work:
Discrete Fourier transforms will be HUGELY useful for the fine structure (not the Morse potential-like function which is already easily analysed). For some reason, the frequency domain always shows more processable data. I will try to find out how to do this.

Code: Select all

import golly as g

class soup:
    def __init__(self,kmax):
        self.pop=[0]*kmax #Population of pattern
        self.dpop=[0]*kmax #Change of population wrt previous generation
        

kmax=int(g.getstring("Maximum number of generations:","2000"))
ns=int(g.getstring("Number of scans:","20000"))
alpha=soup(kmax)
k=0
g.show("Press q to stop at any point.")
for j in xrange(ns):
    g.new("Spectrum_scan_"+str(j))
    g.select([0,0,16,16]) # These two lines can be substituted to analyse different patterns
    g.randfill(30)
    temp_pop=int(g.getpop())
    alpha.pop[0]+=temp_pop
    alpha.dpop[0]=0 # IMPORTANT: Ignore first data point for FT
    for i in xrange(kmax-1):
        g.run(1)
        alpha.pop[i+1]+=int(g.getpop())
        alpha.dpop[i+1]+=int(g.getpop())-temp_pop
        temp_pop=int(g.getpop())
    try:
        g.select(g.getselrect())
        g.clear(0)
    except: pass
    event=g.getevent()
    if event.startswith("key q none"):
        ns=j
        break
File=open("geneascopy_ns_%d.txt"%ns,"w+")
File.write("Gen\tPop\tdPop\tln(dPop^2)/2\tln((Pop-minpop)+1)\tabs(dPop)\t\tNS=%d\n"%ns)
minPop=min(alpha.pop)
for m in xrange(kmax):
    File.write(str(m)+"\t"+str(alpha.pop[m])+"\t"+str(alpha.dpop[m])+"\t=ln(C%d^2)/2\t"%(m+2)+"=ln(B%d-"%(m+2)+"%d+1)\t"%minPop+"=abs(C%d)\n"%(m+2))
#
File.write("\nPop{")
for n in xrange(kmax):
    File.write("%f,"%(float(alpha.pop[n])/float(ns)))
File.write("}\n\ndPop{")
for op in xrange(kmax):
    File.write("%f,"%(float(alpha.dpop[op])/float(ns)))
File.write("}")
#
File.close()
g.exit()
Last edited by Rhombic on March 20th, 2019, 9:01 am, edited 5 times in total.
SoL : FreeElectronics : DeadlyEnemies : 6a-ite : Rule X3VI
what is “sesame oil”?

User avatar
Rhombic
Posts: 1057
Joined: June 1st, 2013, 5:41 pm

Re: Geneascopy

Post by Rhombic » March 18th, 2019, 12:15 pm

I have run 120000 scans through 8000 generations in CGoL and I can confirm that, for a larger number of scans, the fine structure becomes identifiable for more generations (even down to generation 5000).

I have noticed some interesting behaviour by plotting ln((pop-minpop)+1) vs ln(changepop^2)/2, with CGoL showing areas that correspond to artificially-random patterns (i.e. initial generations of a soup) and areas that correspond to natural-like behaviour (such as 60 generations after the interaction between a Pi and a constellation).
This addresses the issue of obtaining useful "natural-like" soups for synthesis and its premise could be used in design of parents! (the check is computationally really cheap)

More things, I have found something so far unique to CGoL (vs tlife, Snowflakes, etc) which is that after the initial dip bottleneck, the population actually increases again. Also, tlife, Snowflakes and even HighLife, rules that seem very close to CGoL otherwise, do NOT have the same concerted evolution behaviour when plotting ln(changepop^2)/2 vs generation, implying that it has a less complex dynamic behaviour.

B37/S23 (DryLife) shows an increase in population after an initial dip, as expected. However, as opposed to CGoL, after the initial dip in ln(changepop^2)/2 vs generation, it shows a much earlier onset of dispersion in change in population, which keeps steadily increasing instead of decreasing.

So far, the only rule to show similar behaviour to CGoL is B37e/S23, which, beyond the collective [linear growths]*P(linear growths) after gen 5470, shows a similar dip followed by a rapid increase with loads of fine structure, although obviously reaching higher populations. A less pronounced concerted evolution region than with CGoL.
SoL : FreeElectronics : DeadlyEnemies : 6a-ite : Rule X3VI
what is “sesame oil”?

User avatar
Rhombic
Posts: 1057
Joined: June 1st, 2013, 5:41 pm

Re: Geneascopy

Post by Rhombic » March 20th, 2019, 8:59 am

I am using this as a brute force discrete Fourier Transform script:

Code: Select all

//
// dft.c - Simple brute force DFT
// Written by Ted Burke
// Last updated 7-12-2013
//
// To compile:
//    gcc dft.c -o dft.exe
// 
// To run:
//    dft.exe
//
// Modified by Rhombic, 20-03-2019
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#define N 3999 //Number of generations. Remove first data point and reduce this number by one when working with dPop.
#define PI2 6.2832

int main()
{
    // time and frequency domain data arrays
    int n, k;             // indices for time and frequency domains
    const float x[]={};   // discrete-time signal, x, ordered values of Pop or dPop
    float Xre[N], Xim[N]; // DFT of x (real and imaginary parts)
    float P[N];           // power spectrum of x
    
     
    // Calculate DFT of x using brute force
    for (k=0 ; k<N ; ++k)
    {
        // Real part of X[k]
        Xre[k] = 0;
        for (n=0 ; n<N ; ++n) Xre[k] += x[n] * cos(n * k * PI2 / N);
         
        // Imaginary part of X[k]
        Xim[k] = 0;
        for (n=0 ; n<N ; ++n) Xim[k] -= x[n] * sin(n * k * PI2 / N);
         
        // Power at kth frequency bin
        P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
    }
     
    // Output results to MATLAB / Octave M-file for plotting
    FILE *f = fopen("dftplots.m", "w");
    fprintf(f, "n = [0:%d];\n", N-1);
    fprintf(f, "x = [ ");
    for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
    fprintf(f, "];\n");
    fprintf(f, "Xre = [ ");
    for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xre[k]);
    fprintf(f, "];\n");
    fprintf(f, "Xim = [ ");
    for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xim[k]);
    fprintf(f, "];\n");
    fprintf(f, "P = [ ");
    for (k=0 ; k<N ; ++k) fprintf(f, "%f ", P[k]);
    fprintf(f, "];\n");
    fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
    fprintf(f, "xlim([0 %d])\n", N-1);
    fprintf(f, "subplot(3,1,2)\nplot(n,Xre,n,Xim)\n");
    fprintf(f, "xlim([0 %d])\n", N-1);
    fprintf(f, "subplot(3,1,3)\nstem(n,P)\n");
    fprintf(f, "xlim([0 %d])\n", N-1);
    fclose(f);
	// Log values
	for (n=0 ; n<N ; ++n) Xre[n] = log(Xre[n]*Xre[n])/2;
	for (n=0 ; n<N ; ++n) Xim[n] = log(Xim[n]*Xim[n])/2;
	for (n=0 ; n<N ; ++n) P[n] = log(P[n]);
	FILE *fg	= fopen("dftplots_log.m", "w");
    fprintf(fg, "n = [0:%d];\n", N-1);
    fprintf(fg, "x = [ ");
    for (n=0 ; n<N ; ++n) fprintf(fg, "%f ", x[n]);
    fprintf(fg, "];\n");
    fprintf(fg, "Xre = [ ");
    for (k=0 ; k<N ; ++k) fprintf(fg, "%f ", Xre[k]);
    fprintf(fg, "];\n");
    fprintf(fg, "Xim = [ ");
    for (k=0 ; k<N ; ++k) fprintf(fg, "%f ", Xim[k]);
    fprintf(fg, "];\n");
    fprintf(fg, "P = [ ");
    for (k=0 ; k<N ; ++k) fprintf(fg, "%f ", P[k]);
    fprintf(fg, "];\n");
    fprintf(fg, "subplot(3,1,1)\nplot(n,x)\n");
    fprintf(fg, "xlim([0 %d])\n", N-1);
    fprintf(fg, "subplot(3,1,2)\nplot(n,Xre,n,Xim)\n");
    fprintf(fg, "xlim([0 %d])\n", N-1);
    fprintf(fg, "subplot(3,1,3)\nstem(n,P)\n");
    fprintf(fg, "xlim([0 %d])\n", N-1);
    fclose(fg);
     
    // exit normally
    return 0;
}
The results can be viewed with Octave. So far, getting very interesting results for dPop, but population needs more scans. Will keep you posted.
EDIT:

Preliminary DFT results of dPop for 30% density 16x16 soups in CGoL show signs of periodicity:
dPopDFT.JPG
dPopDFT.JPG (100.97 KiB) Viewed 2428 times
Broad peaks at ca. 21 gen^-1 and 59 gen^-1
SoL : FreeElectronics : DeadlyEnemies : 6a-ite : Rule X3VI
what is “sesame oil”?

User avatar
Andrew
Moderator
Posts: 766
Joined: June 2nd, 2009, 2:08 am
Location: Melbourne, Australia
Contact:

Re: Geneascopy

Post by Andrew » March 20th, 2019, 6:45 pm

I've written a Lua version of Rhombic's geneasco.py script:

Code: Select all

local g = golly()
local ov = g.overlay
local ovt = g.ovtable
local op = require "oplus"
local gp = require "gplus"
local int = gp.int

local kmax = tonumber(g.getstring("Maximum number of generations:","5000"))
local ns = tonumber(g.getstring("Number of scans:","20000"))
local pop = {}
local popchange = {}
for m = 0, kmax-1 do
    pop[m] = 0
    popchange[m] = 0
end

for j = 1, ns do
    g.new("Spectrum_scan_"..j)
    g.select({0,0,16,16})
    g.randfill(30)
    g.update()  -- see progress
    local temp_pop = tonumber(g.getpop())
    pop[0] = pop[0] + temp_pop
    popchange[0] = ""
    for i = 1, kmax-1 do
        g.run(1)
        local p = tonumber(g.getpop())
        pop[i] = pop[i] + p
        popchange[i] = popchange[i] + (p - temp_pop)
        temp_pop = p
    end
end

--[[ enable these lines to save data in a .txt file:

local f = io.open("gscopy_ns_"..ns..".txt","w")
f:write("Gen\tPop\tChangePop\n")
for m = 0, kmax-1 do
    f:write(m.."\t"..pop[m].."\t"..popchange[m].."\n")
end
f:close()

--]]

--------------------------------------------------------------------------------

function xyplot(width, height, bgcolor)
    ov("create "..width.." "..height)
    ov(bgcolor or op.white)
    ov("fill")

    -- return a table that makes it easy to create and save XY plots
    local p = {}

    p.tborder = 60  -- border above plot
    p.bborder = 60  -- border below plot
    p.lborder = 60  -- border left of plot
    p.rborder = 60  -- border right of plot

    p.title = "untitled"
    p.xlabel = ""
    p.ylabel = ""
    p.titlegap = 10
    p.xlabelgap = 10
    p.ylabelgap = 10
    p.numgap = 10
    p.lines = false     -- draw connected lines?

    p.axescolor = "rgba 0 0 0 255"      -- black axes
    p.textcolor = "rgba 0 0 0 255"      -- black text
    p.plotcolor = "rgba 0 0 255 255"    -- blue dots
    
    p.textfont = "font 12 default-bold"
    p.numfont = "font 10 default"

    p.borders = function(top, left, bottom, right)
        p.tborder = top
        p.bborder = bottom
        p.lborder = left
        p.rborder = right
    end

    p.xaxis = function(min, max, tick)
        p.minx = min
        p.maxx = max
        if min == max then p.maxx = max + 1 end -- avoid division by zero
        p.tickx = tick
    end

    p.yaxis = function(min, max, tick)
        p.miny = min
        p.maxy = max
        if min == max then p.maxy = max + 1 end -- avoid division by zero
        p.ticky = tick
    end

    p.draw = function(xmin, xmax, xinc, yfunc)
        -- calculate length of axes, origin pixel, and scale on each axis
        local xlen = width - (p.lborder + p.rborder)
        local ylen = height - (p.tborder + p.bborder)
        local ox = p.lborder
        local oy = p.tborder + ylen
        local yscale = (p.maxy - p.miny) / ylen
        local xscale = (p.maxx - p.minx) / xlen

        ov("blend 2")   -- for text and lines
    
        -- draw axes
        ov(p.axescolor)
        op.draw_line(ox, oy, xlen+ox, oy)
        op.draw_line(ox, oy, ox, -ylen+oy)
        
        -- draw title and labels
        ov(p.textcolor)
        ov(p.textfont)
        ov("textoption align center")
        local wd, ht
        if #p.title > 0 then
            wd, ht = op.maketext(p.title)
            op.pastetext((xlen - wd)//2 + ox, -ylen - p.titlegap - ht + oy)
        end
        if #p.xlabel > 0 then
            wd, ht = op.maketext(p.xlabel)
            op.pastetext((xlen - wd)//2 + ox, p.xlabelgap + oy)
        end
        if #p.ylabel > 0 then
            wd, ht = op.maketext(p.ylabel)
            -- rotate this text 90 degrees anticlockwise
            op.pastetext(-p.ylabelgap - ht + ox, -(ylen - wd)//2 + oy, op.racw)
        end
        
        -- draw numbers at ends of axes
        ov(p.numfont)
        wd, ht = op.maketext(""..p.minx)
        op.pastetext(-wd//2 + ox, p.numgap + oy)

        wd, ht = op.maketext(""..p.maxx)
        op.pastetext(xlen - wd//2 + ox, p.numgap + oy)

        wd, ht = op.maketext(""..p.miny)
        op.pastetext(-wd - p.numgap + ox, -ht//2 + oy)

        wd, ht = op.maketext(""..p.maxy)
        op.pastetext(-wd - p.numgap + ox, -ylen - ht//2 + oy)

        -- draw optional numbers and tick marks on axes
        if p.tickx > 0 then
            local x = p.minx - (p.minx % p.tickx)
            while x + p.tickx < p.maxx do
                x = x + p.tickx
                local tickpos = int((x - p.minx) / xscale) + ox
                op.draw_line(tickpos, oy, tickpos, oy+5)
                wd, ht = op.maketext(x)
                op.pastetext(tickpos - wd//2, p.numgap + oy)
            end
        end
        if p.ticky > 0 then
            local y = p.miny - (p.miny % p.ticky)
            while y + p.ticky < p.maxy do
                y = y + p.ticky
                local tickpos = -int((y - p.miny) / yscale) + oy
                op.draw_line(ox, tickpos, ox-5, tickpos)
                wd, ht = op.maketext(y)
                op.pastetext(-wd - p.numgap + ox, tickpos - ht//2)
            end
        end
        
        -- plot the data
        ov(p.plotcolor)
        local oldx =  int((xmin - p.minx) / xscale) + ox
        local oldy = -int((yfunc(xmin) - p.miny) / yscale) + oy
        for xval = xmin, xmax, xinc do
            local x =  int((xval - p.minx) / xscale) + ox
            local y = -int((yfunc(xval) - p.miny) / yscale) + oy
            if p.lines then
                op.draw_line(oldx, oldy, x, y)
                oldx = x
                oldy = y
            else
                -- draw a small "+" mark
                ovt{"set", x, y, x-1, y, x+1, y, x, y-1, x, y+1}
            end
        end
    end

    p.save = function(pngpath)
        if pngpath == nil then
            -- prompt user for png file
            g.update()
            pngpath = g.savedialog("Save plot as PNG file", "PNG (*.png)|*.png", "", "plot.png")
            if #pngpath == 0 then return end    -- user hit Cancel
        end
        -- save plot in given png file
        ov("save 0 0 "..width.." "..height.." "..pngpath)
    end

    p.delete = function()
        ov("delete")
    end

    return p
end

--------------------------------------------------------------------------------

function find_limits(t)
    local minval = math.maxinteger
    local maxval = math.mininteger
    for _,v in pairs(t) do
        if v < minval then minval = v end
        if v > maxval then maxval = v end
    end
    return minval, maxval
end

--------------------------------------------------------------------------------

minpop, maxpop = find_limits(pop)
maxdigits = int(math.log10(maxpop)) + 1
poptick = int(10^(maxdigits-1))
if minpop + 4 * poptick > maxpop then poptick = poptick // 10 end
if minpop % poptick > 0 then minpop = minpop - (minpop % poptick) end
if maxpop % poptick > 0 then maxpop = maxpop - (maxpop % poptick) + poptick end

-- now use overlay to create 2 png files with the desired plots

plot1 = xyplot(1000, 800)
plot1.title = g.getrule()..", 16x16 at 30%\nscans = "..ns
plot1.xlabel = "Generation"
plot1.ylabel = "Cumulative Population"
plot1.xlabelgap = 30
plot1.ylabelgap = 10 * (maxdigits+1)
plot1.lborder = plot1.ylabelgap + 40
plot1.rborder = 40
plot1.xaxis(0, kmax, 1000)
plot1.yaxis(minpop, maxpop, poptick)
-- plot1.lines = true
plot1.draw(0, kmax-1, 1, function(x) return pop[x] end)
plot1.save("plot1_ns="..ns.."_g="..kmax..".png")
plot1.delete()

plot2 = xyplot(1000, 600)
plot2.title = "ln(changepop^2)/2 vs generation"
plot2.xaxis(0, math.min(kmax-1,800), 100)
plot2.yaxis(0, 14, 2)
-- plot2.lines = true
-- note that we must start x at 1 because popchange[0] = ""
plot2.draw(1, math.min(kmax-1,800), 1,
function(x)
    if popchange[x] == 0 then return 0 else return math.log(popchange[x]^2)/2 end
end)
plot2.save("plot2_ns="..ns.."_g="..kmax..".png")
plot2.delete()
Instead of saving the results in a text file it uses overlay commands to create 2 plots that are saved in .png files.

The plots are created using a function called xyplot, so people might like to pull out that code and use it in other projects that need to display results in some sort of plot. The function is quite versatile and easy to use. For example, this piece of code:

Code: Select all

test = xyplot(300, 300, op.black)
test.title = "y = x^2"
test.axescolor = op.white
test.textcolor = op.white
test.plotcolor = op.yellow
test.borders(40, 40, 40, 40)
test.numfont = "font 9 mono-bold"
test.lines = true
test.xaxis(0, 10, 1)
test.yaxis(0, 100, 10)
test.draw(0, 10, 0.1, function(x) return x^2 end)
test.save("test.png")
will create this PNG file:
test.png
test.png (7.73 KiB) Viewed 2411 times

User avatar
Rhombic
Posts: 1057
Joined: June 1st, 2013, 5:41 pm

Re: Geneascopy

Post by Rhombic » November 5th, 2019, 7:07 pm

I apologise for having ignored this project (and CGoL to some extent, actually) for the past few months.
I started my PhD about one and a half months ago and now I am pretty busy. I will attempt to devote some time to this when I have a bit more free time, if that is something I can expect from the coming years!

For those who might wish to continue in this line:
I remember I was unsure whether the discrete Fourier transform script I used actually worked properly for periodic patterns that were out-of-phase wrt the reference wave to which it was assigned, so caveat emptor and keep an eye. And yes, discrete FTs will be the way to go to get any meaningful info out of this, but the dependence on phase coherence with the script I used leads to poor results.
I haven't been able to get to grips with why ln(changepop^2)/2 vs generation shows pseudo-random behaviour ONLY after ca. 160 generations, not before. This seems related to the very nature of the CA rule but I have no idea how to approach it.
SoL : FreeElectronics : DeadlyEnemies : 6a-ite : Rule X3VI
what is “sesame oil”?

User avatar
testitemqlstudop
Posts: 1241
Joined: July 21st, 2016, 11:45 am
Location: in catagolue
Contact:

Re: Geneascopy

Post by testitemqlstudop » November 7th, 2019, 2:19 am

Oh wow you're back?
Rhombic wrote:
March 18th, 2019, 12:15 pm

So far, the only rule to show similar behaviour to CGoL is B37e/S23, which, beyond the collective [linear growths]*P(linear growths) after gen 5470, shows a similar dip followed by a rapid increase with loads of fine structure, although obviously reaching higher populations. A less pronounced concerted evolution region than with CGoL.
Have you tried analyzing rules one isotropic transition from CGOL, like b38/s23 or b3/s238?

User avatar
Hdjensofjfnen
Posts: 1368
Joined: March 15th, 2016, 6:41 pm
Location: r cis θ

Re: Geneascopy

Post by Hdjensofjfnen » November 7th, 2019, 10:37 pm

testitemqlstudop wrote:
November 7th, 2019, 2:19 am
Rhombic wrote:
March 18th, 2019, 12:15 pm
B37e/S23
Have you tried analyzing rules one isotropic transition from CGOL, like b38/s23 or b3/s238?
... B37e/S23 is one isotropic transition from CGOL. Also, Rhombic isn't back per se: they're still quite busy from the looks of it.
EDIT: toroidalet has informed me that the purpose of that post was probably to ask about other Life-like rules. Apologies.
"A man said to the universe:
'Sir, I exist!'
'However,' replied the universe,
'The fact has not created in me
A sense of obligation.'" -Stephen Crane

Code: Select all

x = 7, y = 5, rule = B3/S2-i3-y4i
4b3o$6bo$o3b3o$2o$bo!

Post Reply