Cellular Automata Rule Generating Polynomials

Rodrigo A. Obando, Columbus State University

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

Cellular Automata rules are represented by integers where we encode the output of the function without knowing the details on how it might be implemented. The CellularAutomaton function in Mathematica only requires these integers, along with the values of r and k, to evolve rules for a given input. The BooleanFunction function in Mathematica provides a polynomial of the function given its representing integer which we can use to get the functions when k = 2. Unfortunately, there is no Mathematica function that does the same for other values of k. We present a uniform way to generate these functions using certain kind of polynomials. Every function in a rule space can be generated by changing the coefficients of its corresponding polynomial. We find that when the value of k is prime the polynomial is relatively simple and straightforward. When the value of k is composite we need a composition of polynomials, each based on a prime factor of k. A possible use of these polynomials is the creation of Probabilistic Cellular Automata by defining one or more of the coefficients as random variables. For every step, for every cell, the random coefficients will take on different values and hence generate different rules. This may allow for a better analysis of these stochastic processes.

Introduction

We have used integers to encode the rules in evolving cellular automata. Our interest has been in looking at how particular rules create behaviors whether simple, cyclic, complex or random. In particular, when using these systems we are able to find the functions that generate these encoded rules for k = 2, the binary case. We have published these functions as Stephen Wolfram did in many of his publications and in particular in his book A New Kind of Science. We can use the Wolfram Language function BooleanFunction to find the function that generates any rule for the case of k = 2.
What about the case for k > 3? Can we find those functions? Even though the functions could be of any type, I propose the use a particular kind of polynomials to find these functions.

Using polynomials to find 1-variable modulo k functions

We will restrict our polynomials to be of k-1 degree at the most. Therefore, for modulo 2 functions we will use polynomials of degree 1 at the most.

A 1-variable Boolean polynomial where x is the input variable, y is the output variable and (Subscript[c, 0]) and (Subscript[c, 1]) are coefficients.

(Subscript[c, 0])+(Subscript[c, 1])x==y

Number of maps for modulo 2 with one variable.

NumberOfMaps[2,1]

There are 4 possible 1-variable functions, and we may use the polynomial above to find the coefficients by solving a set of equations for each map.

MatrixForm[{Solve[(Subscript[c, 0 ])+(Subscript[c, 1])0==0&&(Subscript[c, 0])+(Subscript[c, 1])1==0,{(Subscript[c, 0]),(Subscript[c, 1])},Modulus->2]->0,
Solve[(Subscript[c, 0 ])+(Subscript[c, 1])0==1&&(Subscript[c, 0])+(Subscript[c, 1])1==0,{(Subscript[c, 0]),(Subscript[c, 1])},Modulus->2]->1+x,
Solve[(Subscript[c, 0 ])+(Subscript[c, 1])0==0&&(Subscript[c, 0])+(Subscript[c, 1])1==1,{(Subscript[c, 0]),(Subscript[c, 1])},Modulus->2]->x,
Solve[(Subscript[c, 0 ])+(Subscript[c, 1])0==1&&(Subscript[c, 0])+(Subscript[c, 1])1==1,{(Subscript[c, 0]),(Subscript[c, 1])},Modulus->2]->1},TableHeadings->{{"{0,0}","{0,1}","{1,0}","{1,1}"},{"->"}}]

Solution for output {0,1}

(Subscript[c, 0])+(Subscript[c, 1])x/.Flatten[Solve[(Subscript[c, 0 ])+(Subscript[c, 1])0==1&&(Subscript[c, 0])+(Subscript[c, 1])1==0,{(Subscript[c, 0]),(Subscript[c, 1])},Modulus->2]]

We can generalize solving these equations for all modulo polynomials as follows.

Computes a Vandermonde matrix for given modulo m.

OneVarModuloMatrix[modulo_]:=Append[#,1]&/@(Mod[((#)^(Range[modulo-1,1,-1])),modulo]&/@Range[modulo-1,0,-1]);

Computes the coefficients for a polynomial for a given output vector.

OneVarModuloCoeff[outpuVector_]:=Module[{L=Length@outpuVector},
LinearSolve[OneVarModuloMatrix[L],outpuVector,Modulus->L]];

Computes a one-variable modulo m polynomial for a given output vector.

OneVarModuloFunction[outpuVector_]:=Module[{L=Length@outpuVector},
Clear@x;
OneVarModuloCoeff[outpuVector].(x^(Range[L-1,0,-1]))];

Table for all 1-variable polynomials modulo 2 for all possible outputs

MatrixForm[{#,OneVarModuloFunction@#}&/@Table[IntegerDigits[n,2,2],{n,0,3}],TableHeadings->{None,{"output","function"}}]

Table for all 1-variable polynomials modulo 3 for all possible outputs

MatrixForm[{#,OneVarModuloFunction@#}&/@Table[IntegerDigits[n,3,3],{n,0,26}],TableHeadings->{None,{"output","function"}}]

When we try this for modulo 4 we find that there are no solutions for all possible output vectors. In fact, whenever we try for modulo k, k a composite number, not all output vectors have solutions.
It seems that these 1-variable modulo k functions can be realized as polynomials of degree k-1 or less only if k is prime.

Theorem

All possible modulo k 1-variable functions can be realized as polynomials of degree k-1 if and only if k is prime.

Proof

In order to prove the theorem is true, we need to show that all output vectors are obtainable given an input vector (vector of coefficients modulo k.) We are creating a map from modulo k vectors of dimension k (output vectors) to modulo k vectors of dimension k (coefficient vectors.)

Let V be a Vandermonde matrix over Z(k) built from a vector of values from 0 to k-1. Let coefficient vector c be a modulo k vector of k dimensions and the output vector o be a modulo k vector of k dimensions where the entries of  o and c come from the additive group Z(k) then

o=V.c

We need to show that transformation V is bijective. Therefore, we need to show that V is both injective and surjective.

We first show that V is injective, so let's assume

o=V.a
o=V.b
V.a=V.b
a=(V^(-1))V.b
a=b

This implies that the transformation V must invertible for the transformation to be injective and therefore its determinant should be nonzero. 
Given how the Vandermonde matrix V is constructed, we may compute the determinant by the following function.

Function to compute the determinant of the special version of the Vandermonde matrix used.

VDet[k_]:=Mod[Times@@(Mod[Factorial[#],k]&/@Range[1,k-1]),k];

Since we are operating on a ring modulo k, we have that if k is prime the ring has no zero divisors but if k is composite it does have zero divisors. Therefore, if k is composite V is not invertible because its determinant is 0 under modulo k since at least one zero divisor will be a term in the computation of the determinant.

Table showing modulo k and the determinant of the V matrix for that modulo using both VDet and Det functions.

MatrixForm[{#,VDet[#],Det[OneVarModuloMatrix@#,Modulus->#]}&/@Range[2,23],TableHeadings->{None,{"modulo","VDet","Det"}}]

We second show that V is surjective, let's assume

o=V.a
a=(V^(-1)).o

The coefficient vector a exists in the k dimensional vector space modulo k only if the transformation V has an inverse; it follows from the result above that this is only true when k is prime.
QED

Polynomials of higher degree

Continuing this result the following functions were integrated in a Paclet named Obando`CAGenetics which is under development. This package is introduced in another talk, ‘Introduction to a Unified Model of Cellular Automata,’ at 1:30 pm.

In the case of binary cellular automata we may use BooleanFunction as follows.

PrintableEquation[BooleanFunction[110,{(Subscript[a, -1]),(Subscript[a, 0]),(Subscript[a, 1])}]]
PrintableEquation[BooleanFunction[30,{(Subscript[a, -1]),(Subscript[a, 0]),(Subscript[a, 1])}]]
PrintableEquation[BooleanFunction[90,{(Subscript[a, -1]),(Subscript[a, 0]),(Subscript[a, 1])}]]
PrintableEquation[BooleanFunction[37,{(Subscript[a, -1]),(Subscript[a, 0]),(Subscript[a, 1])}]]

For the case where k > 2, the function in the Obando`CAGenetics paclet is ModulusPFunction and it takes an integer as the rule, the k value and the variables we want to use. The number of variables determines the value of r. In the following example we use various rules in k=3, r=1.

ModulusPFunction[3456,3,{(Subscript[a, -1]),(Subscript[a, 0]),(Subscript[a, 1])}]
ModulusPFunction[5999036270526,3,{(Subscript[a, -1]),(Subscript[a, 0]),(Subscript[a, 1])}]
ModulusPFunction[1059299700530,3,{(Subscript[a, -1]),(Subscript[a, 0]),(Subscript[a, 1])}]
ModulusPFunction[4932442728002,3,{(Subscript[a, -1]),(Subscript[a, 0]),(Subscript[a, 1])}]

RuleEquation wraps ModulusPFunction so that you don’t specifically add the list of variables, instead it uses the value of r to provide the list of variables.

RuleEquation[3456,3,1]
RuleEquation[5999036270526,3,1]
RuleEquation[1059299700530,3,1]
RuleEquation[4932442728002,3,1]

The function ModulusPFunction makes use of the following functions.

Function to create the list of variables for a given value of r.

RuleVariables[1]

Function to create the list of all powers in the polynomial.

ModulusPCore[3,RuleVariables[1]]

Function to create the transformation matrix to map the output of the rule to the coefficients in the polynomial.

ModulusPMatrix[3,RuleVariables[1]]//MatrixForm

Function to create the list of the coefficients to be applied to the Core computed above using LinearSolve.

ModulusPCoeff[5999036270526,3,RuleVariables[1]]

There are 2 other functions where we can convert from rule output to coefficients and vice-versa.

Function to convert a rule output to coefficients. The coefficients are returned encoded in an integer, similar to what we use for rule output.

OutputToCoeff[5999036270526,3,1]
RuleExpand[OutputToCoeff[5999036270526,3,1],3,1]

Function to convert a list of coefficients to rule output.

CoeffToOutput[1030840625691,3,1]
RuleExpand[1030840625691,3,1]

Polynomials for Composite k

Since these polynomials do not work for composite k, I am working on a polynomial of mix modulus. We can find the integers that make up a composite number by using IntegerDigits. This is still work in progress.

Probabilistic Cellular Automata

A possible use of these polynomials is to explore probabilistic cellular automata. The variables in the polynomials are the inputs from to the rule from one evolution step to the next. If we make one or more of the coefficients to be random variables, we can look at how those variables generate different rules and how that affects the overall behavior of the cellular automata evolution. At every evaluation of every cell in the cellular automata we can sample those variables and change the actual rule used in the evolution. In effect, we would take the output of the ModulusPCoeff function and make some of the elements of the list to be randomly generated then used the list to compute the rule using the ModulusPMatrix function.

Simple Example

Function to convert a list of coefficients to a cellular automata rule as a list of rules instead of an integer. You may define which coefficients are random.

PRule[coeff_,randomCoeff_List,k_,r_]:=MapThread[#1:>#2[[RandomInteger[{1,k^Length[randomCoeff]}]]]&,{Tuples[Reverse[Range[0,k-1]],1/Log[Length[RuleExpand[coeff,k,r]],k]],Transpose[Map[RuleExpand[CoeffToOutput[FromDigits[ReplacePart[RuleExpand[coeff,k,r],#],k],k,r],k,r]&,MapThread[#1->#2&,{randomCoeff,#}]&/@Tuples[Range[0,k-1],Length[RuleExpand[randomCoeff]]]]]}]

Manipulate to show evolution of rules, one deterministic, one probabilistic, k = 2, r = 1 (ECA.)

Manipulate[Grid[{{ModulusPCore[2,RuleVariables[1]]},{GraphicsGrid[{{ArrayPlot[CellularAutomaton[{theRule,2,1},{{1},0},100],PlotLabel->theRule],ArrayPlot[CellularAutomaton[PRule[OutputToCoeff[theRule,2,1],coeff,2,1],{{1},0},100],PlotLabel->"Probabilistic"]}}]}}],{{theRule,90},0,255,1},{{coeff,{3,5}},Range[8],ControlType->TogglerBar}]

Manipulate to show evolution of rules, one deterministic, one probabilistic, k = 3, r = 1.

Manipulate[GraphicsGrid[{{ArrayPlot[CellularAutomaton[{theRule,3,1},{{1},0},100],PlotLabel->theRule],ArrayPlot[CellularAutomaton[PRule[OutputToCoeff[theRule,3,1],{3,5},3,1],{{1},0},100],PlotLabel->"Probabilistic"]}}],{theRule,2727936019226}]

Conclusion

Finding the polynomial that defines the rule function of a cellular automata allows for other kinds of analysis of cellular automata behavior. It provides a new avenue to explore the many types of cellular automata from the point of view of the polynomials and not just from the rule output standpoint.
This also plays an important role in the Unified Model of Cellular Automata I will present today at 1:30 pm. It makes the creation of equation trees of rules and analysis based on the polynomials that make up the gene functions in the rules.

Code for the talk (snipet from Obando`CAGenetics)

RuleExpand::usage = "RuleExpand[rule, k, r] expands the rule into integer digits of modulo k at the appropriate size for the value of r.";

RuleExpand[rule_, k_:2, r_:1] := IntegerDigits[rule, k, k ^ (2 r + 1)
    ];
RuleVariables::usage = "RuleVariables[r] returns the variables for a rule in rule space with r.";

RuleVariables[r_:1] := Flatten[{Table[Subscript[Global`a, x], {x, -Ceiling[
    r], -1}], {Subscript[Global`a, 0]}, Table[Subscript[Global`a, x], {x,
     1, IntegerPart[r]}]}];

GeneVariables::usage = "GeneVariables[r] returns the variable for the genes in rule space with r.";

GeneVariables[r_:1] := Flatten[{Table[Subscript[Global`a, x], {x, -Ceiling[
    r], -1}], Table[Subscript[Global`a, x], {x, 1, IntegerPart[r]}]}];

RuleEquation::usage = "RuleEquation[rule, k, r] returns the rule equation for the rule in rule space with k and r.";

RuleEquation[rule_, k_:2, r_:1] := If[k == 2,
    BooleanFunction[rule, RuleVariables[r]]
    ,
    If[PrimeQ[k],
        ModulusPFunction[rule, k, RuleVariables[r]]
        ,
        RuleVariables[r]
    ]
];

GeneEquation::usage = "GeneEquation[gene, k, r] returns the gene equation for the gene in rule space with k and r.";

GeneEquation[gene_, k_:2, r_:1] := If[k == 2,
    BooleanFunction[gene, GeneVariables[r]]
    ,
    If[PrimeQ[k],
        ModulusPFunction[gene, k, GeneVariables[r]]
        ,
        GeneVariables[r]
    ]
];

CoeffToOutput::usage = "CoeffToOutput[coeff, k, r] returns the output signature given a list of coefficients.";

CoeffToOutput[coeff_, k_, r_] := FromDigits[Reverse[Mod[ModulusPMatrix[
    k, RuleVariables[r]] . Reverse[RuleExpand[coeff, k, r]], k]], k];

OutputToCoeff::usage = "OutputToCoeff[output, k, r] returns the coefficient list given an output signature.";

OutputToCoeff[output_, k_, r_] := FromDigits[ModulusPCoeff[output, k,
     RuleVariables[r]], k];

ModulusPCore::usage = "ModulusPCore[modulus, vars] returns the list of variables and their powers in the given modulus.";

ModulusPCore[modulus_, vars_] := Module[{cores, finalCore},
    If[Length[vars] == 1,
        finalCore = Times @@ #& /@ ({Prepend[Table[vars[[1]], modulus
             - 1], 1][[#]] ^ (# - 1)}& /@ Range[modulus])
        ,
        cores = ModulusPCore[modulus, {#}]& /@ vars; finalCore = Fold[
            Flatten[Outer[Times, #1, #2]]&, cores]
    ];
    finalCore
];

ModulusPMatrix::usage = "ModulusPMatrix[modulus, vars] returns the transformation matrix given the modulus and list of variables.";

ModulusPMatrix[modulus_, vars_] := Mod[ModulusPCore[modulus, vars], modulus
    ] /. Inner[#1 -> #2&, vars, #, List]& /@ (IntegerDigits[#, modulus, Length[
    vars]]& /@ Table[n, {n, 0, modulus ^ Length[vars] - 1}]);

ModulusPCoeff::usage = "ModulusPCoeff[output, modulus, vars] returns the coefficients for the polynomial that generates the given output for the modulus and variables.";

ModulusPCoeff[output_, modulus_, vars_] := Reverse[LinearSolve[ModulusPMatrix[
    modulus, vars], Reverse[RuleExpand[output, modulus, (Length[vars] - 1
    ) / 2]], Modulus -> modulus]];

ModulusPFunction::usage = "ModulusPFunction[output, modulus, vars] returns the function of the polynomial that produces the given output for modulus and list of variables.";

ModulusPFunction[output_, modulus_, vars_] := Reverse[ModulusPCoeff[output,
     modulus, vars]] . ModulusPCore[modulus, vars];

ModulusPFunctionFromCoeff::usage = "ModulusPFunctionFromCoeff[coeff, modulus, vars] returns the function of the polynomials for a given list of coefficients with modulus and list of variables.";

ModulusPFunctionFromCoeff[coeff_, modulus_, vars_] := Reverse[coeff] . ModulusPCore[
    modulus, vars];

ModulusKDCoeff::usage = "ModulusKDCoeff[modulus, var, value] returns the Kronecker Delta coefficients for the given modulus and variable and value.";

ModulusKDCoeff[modulus_, var_, value_] := ModulusPCoeff[FromDigits[Table[
    KroneckerDelta[n, value + 1], {n, modulus, 1, -1}], modulus], modulus,
     {var}];

ModulusKDFunction::usage = "ModulusKDFunction[modulus, var, value] returns the Kronecker Delta function for the given modulus and variables and value.";

ModulusKDFunction[modulus_, var_, value_] := Reverse[ModulusKDCoeff[modulus,
     var, value]] . ModulusPCore[modulus, {var}];

PrintableEquation::usage = "PrintableEquation[equation] takes an equation and replaces the boolean operators with common operators.";

PrintableEquation[equation_] := ToExpression[ToString[FullForm[equation
    ] /. Or -> Plus /. And -> Times /. Not -> OverBar /. Xor -> CirclePlus
     /. False -> 0 /. True -> 1]];