Home

The purpose of this markdown is to explore Wolfram’s elementary cellular automata, in terms of rule space, and detmerine if I can find the so-called edge of chaos.

The first thing we have to do is make a function that generatees the rules of interest.

library(pheatmap)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.2.0
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
# Define a rule
#state <- c("111", "110", "101", "100", "011", "010", "001", "000")
state <- lapply(seq(0, 7), function(i) {
    result <- intToBits(i)[1:3] %>% as.numeric() %>% rev() %>% paste(collapse = "")
    return(result)
}) %>% as.character() 

outcome <- lapply(seq(0, 255), function(i) {
    result <- intToBits(i)[1:8] %>% as.numeric()
    return(result)
}) %>% do.call(cbind, .) %>% as_tibble()
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
names(outcome) <- paste("rule", seq(0, 255), sep = "_")

rule_set <- bind_cols(state = state, outcome = outcome)

# The procedure
ApplyRule <- function(v, rule_set, rule_num) {
    # Set the rule
    rule <- tibble(state = state, outcome = rule_set[[rule_num + 2]])
    
    # Iterate through the vector
    result <- lapply(seq_along(v), function(i) {
        # Edge case
        if(i == 1 | i == length(v)) {
            return(0)
        }
        
        # The rule
        curr <- paste(v[i - 1], v[i], v[i + 1], sep = "") 
        out <- rule[rule$state == curr,]$outcome
        return(out)
    }) %>% unlist()
}

Now let’s make Rule 30.

# Define a vector 
v <- c(rep(0, 30), 1, rep(0, 30))

# Define an output matrix
top <- v
mat <- lapply(seq(30), function(i) {
    curr <- ApplyRule(v = v, rule_set = rule_set, rule = 30)
    v <<- curr
    return(curr)
}) %>% do.call(rbind, .)

mat <- rbind(top, mat)


# Visualize
pheatmap::pheatmap(mat, cluster_rows = FALSE, cluster_cols = FALSE)

And now with a less fixed initialization.

# Define a vector 
v <- c(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0), 10))

# Define an output matrix
top <- v
mat <- lapply(seq(30), function(i) {
    curr <- ApplyRule(v = v, rule_set = rule_set, rule = 30)
    v <<- curr
    return(curr)
}) %>% do.call(rbind, .)

mat <- rbind(top, mat)


# Visualize
pheatmap::pheatmap(mat, cluster_rows = FALSE, cluster_cols = FALSE)

Great. Now let’s try to get at Rule Space. We have this rule_set object that we just made. If we transpose it we might be able to get something. 8 dimensions of rules.

rule_mat <- as.data.frame(rule_set)
rownames(rule_mat) <- rule_mat$state
rule_mat$state <- NULL
rule_mat <- t(rule_mat) %>% as.data.frame()
head(rule_mat)
##        000 001 010 011 100 101 110 111
## rule_0   0   0   0   0   0   0   0   0
## rule_1   1   0   0   0   0   0   0   0
## rule_2   0   1   0   0   0   0   0   0
## rule_3   1   1   0   0   0   0   0   0
## rule_4   0   0   1   0   0   0   0   0
## rule_5   1   0   1   0   0   0   0   0

Perfect. Now let’s add a column that corresponds to the Wolfram Class of each rule. Note that we’re only dealing with the first instance of an equivalence. Some rules will output the same thing, so we’re simplifying the data accordingly.

# Wolfram's Class I
class_1 <- c(0, 8, 32, 40, 128, 136, 160, 168)

# Wolfram's Class II
class_2 <- c(1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 19, 23, 24, 25, 26, 27, 28, 29, 33, 34, 35, 36, 37, 38, 41, 42, 43, 44, 46, 50, 51, 56, 57, 58, 62, 72, 73, 74, 76, 77, 78, 94, 104, 108, 130, 132, 134, 138, 140, 142, 152, 154, 156, 162, 164, 170, 172, 178, 184, 200, 204, 232)

# Wolfram's Class III
class_3 <- c(18, 22, 30, 45, 60, 90, 105, 122, 126, 146, 150)

# Wolfram's Class IV
class_4 <- c(54, 106, 110)

Now we prune the rule set.

rownames(rule_mat) <- lapply(rownames(rule_mat), function(i) {
    str_split(i, "_")[[1]][2] 
}) %>% unlist() %>% as.numeric()
rule_mat <- rule_mat[rownames(rule_mat) %in% c(class_1, class_2, class_3, class_4),]

Ok, now let’s add the class.

rule_mat$class <- lapply(rownames(rule_mat), function(i) {
    if(i %in% class_1) return(1)
    if(i %in% class_2) return(2)
    if(i %in% class_3) return(3)
    if(i %in% class_4) return(4)
}) %>% unlist()

Ok, now let’s add Shannon entropy. Low means we’re dealing with mainly or only 0 or 1. High means we’re dealing with even 0 and 1.

ShannonEntropy <- function(x) {
  # Calculate the probability of each unique element in the vector
  prob <- table(x) / length(x)
  # Compute the entropy
  entropy <- -sum(prob * log2(prob))
  return(entropy)
}

rule_mat$entropy <- apply(rule_mat, 1, function(i) {
    result <- i[1:8] # get out the binary
    result <- ShannonEntropy(result)
    return(result)
}) %>% unlist()

head(rule_mat)
##   000 001 010 011 100 101 110 111 class   entropy
## 0   0   0   0   0   0   0   0   0     1 0.0000000
## 1   1   0   0   0   0   0   0   0     2 0.5435644
## 2   0   1   0   0   0   0   0   0     2 0.5435644
## 3   1   1   0   0   0   0   0   0     2 0.8112781
## 4   0   0   1   0   0   0   0   0     2 0.5435644
## 5   1   0   1   0   0   0   0   0     2 0.8112781

Ok, and now for dimension reduction.

library(umap)
set.seed(1234)

dimr <- umap::umap(rule_mat[,1:8])$layout %>% as.data.frame()
names(dimr) <- c("umap1", "umap2")

rule_mat <- dplyr::bind_cols(rule_mat, dimr)

And now we do a few visualizations.

ggplot(rule_mat, aes(x = umap1, y = umap2)) + geom_point(aes(color = as.factor(class)))

Class 3 and 4 show up in pockets. Let’s color by entropy now.

ggplot(rule_mat, aes(x = umap1, y = umap2)) + geom_point(aes(color = entropy))

Not as informative. Let’s try entropy x class.

ggplot(rule_mat, aes(x = entropy, y = class)) + geom_point()

Ok, so half 0 and half 1 allows for the possibility of class 4. But doesn’t guarantee it. What if we filter by Shannon Entropy of 1 and visualize the map again.

rule_mat_trim <- dplyr::filter(rule_mat, entropy > 0.9)

ggplot(rule_mat_trim, aes(x = umap1, y = umap2)) + geom_point(aes(color = as.factor(class)))

That didn’t really do much. I wonder if it has to do with the balance between create and distruct rules. I wonder if this is why the class 3s and 4s cluster together on the map. Without getting too complicated, one thing we can do is just look at the heatmaps of each class.

library(pheatmap)

pheatmap(filter(rule_mat, class == 1)[1:8])

pheatmap(filter(rule_mat, class == 2)[1:8])

pheatmap(filter(rule_mat, class == 3)[1:8])

pheatmap(filter(rule_mat, class == 4)[1:8])

In class 4, we have the 000 and 111 always being destructive. We have most everything else being creative, save either 100, 001, and 010.

Class 3 has a bit more freedom, though we still see 001 and 100 almost always being creative and 000 and 11 almost always being destructive.

Class 2 is all over the place.

Class 1 is creative at 111, 101, and 011, but destructive everywhere else.

Now we’re going to investigate rule space in terms of compression ratios for the output. To do this we’re going to generate the output of a given CA, and see how much we can compress it. The hypothesis is that Class 4 will sit right in the middle, at the edge of chaos.

GenerateEca <- function(rule, num_iter = 200) {
    # Define a vector 
    #v <- c(rep(0, 50), 1, rep(0, 50)) # Center init
    # v <- rep(c(0, 1), 100) # Random init
    v <- c(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0), 10)) # Random init 2

    # Define an output matrix
    top <- v
    mat <- lapply(seq(num_iter), function(i) {
        curr <- ApplyRule(v = v, rule_set = rule_set, rule = rule)
        v <<- curr
        return(curr)
    }) %>% do.call(rbind, .)

    mat <- rbind(top, mat)
    return(mat)
}

# Function to calculate ECA complexity based on compression ratio
# Function to calculate ECA complexity based on compression ratio
CompCx <- function(rule) {
  eca <- GenerateEca(rule)
  
  # Convert the ECA to a binary string
  eca_binary_str <- paste(as.vector(eca), collapse = "") 
  
  # Compress the binary string using gzip compression
  compressed_data <- memCompress(eca_binary_str, type = "gzip")
  
  # Calculate the compression ratio
  original_size <- nchar(eca_binary_str)
  compressed_size <- length(compressed_data)
  compression_ratio <- compressed_size / original_size
  
  # Use the compression ratio as a measure of complexity
  complexity <- compression_ratio
  return(complexity)
}

# Example usage
rule <- 30
complexity <- CompCx(rule)
print(complexity)
## [1] 0.101382

Great! Now let’s run this for all the rules.

rule_mat$complexity <- lapply(seq(nrow(rule_mat)), function(i) {
    curr_rule <- rownames(rule_mat)[i]
    result <- CompCx(rule = as.numeric(curr_rule))
    return(result)
}) %>% unlist()

Now for the visuals.

ggplot(rule_mat, aes(x = class, y = complexity)) + geom_text(aes(label = rownames(rule_mat)))

library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
# Assuming you already have the rule_mat data frame
plot1 <- ggplot(rule_mat, aes(x = umap1, y = umap2)) + geom_point(aes(color = complexity))
plot2 <- ggplot(rule_mat, aes(x = umap1, y = umap2)) + geom_point(aes(color = as.factor(class)))

# Plot the ggplot objects side by side
grid.arrange(plot1, plot2, ncol = 2)

So we find that class 3 cellular automata are on avearge complex, in that they’re harder to compress. I don’t see much of a distinction between class 2 and class 4. Perhaps we need other evaluation metrics here. Class 4 CA can produce packets of information that can move and interact (gliders). They can also in theory produce computation (eg. logic gates). To prove that a rule can do this, we have to literally do it. But perhaps there’s a way that we can actually train a classifier to identify Class 4 behavior, and then produce a score from there? That’s beyond the scope of this markdown, but perhaps that’s where we can go with this.

And finally, to get at the original point of this markdown, if class 3 is the most chaotic behavior and class 1 and class 2 are the most orderly, we do see class 4 showing up next to class 2 and class 3. While we don’t have a ton of data points here (we would need a more complex rule space), this does fit with the interpretation of class 4 sitting at the edge of chaos.