My approach kinda mimics the Markov Algorithm, however it uses Monte-Carlo & Metropolis (like here) approach to sample the patterns instead of a perfect matching.
That make it identical to the original approach (see repo) in the limiting case, when the temperature of the system $T \rightarrow 0$.
States and rules representation
I borrowed several ideas from quantum mechanics, where you can represent a state using vectors. Imagine we have 3 state (black, white, and ..??), then any arbitary state can be expressed as
state = {{1}, {0}, {0}} (* black? *)
Looking back to quantum mechanics this is a wavefunction.
$$
|state\rangle
$$
A pattern in this sence can be represented as a Kroneker product of states at different sites. However, if we do not plan to entangle them, lets assume that they are completely independed (all the time, whatever happend) and place them to 2D array
field = Table[{{1}, {0}, {0}}, {i, 25}, {j, 25}];
In physics we could write them as
$$
|field\rangle = |\psi_{i,j}\rangle...|\psi_{i+1,j}\rangle
$$
Pattern matching
One can think about it if our patter perpresents an energy of interaction inside the corresponding group of state. So we measure the energy for each pattern (rule) at random position on the board and decide if we need to perform a replacment or not. Picking up a random posiiton and changing the states reminds an Ising Model, where Monte-Carlo approach is used to perform local updates.
Then to check (measure), if a cell at i,j position is black of not
Transpose[states[[i,j]]].operator.states[[i,j]] -> "energy"
where operator is a matrix, that must perform a mapping black->1, anything else ->0 without changing an actual state
From the quantum mechanics one can write it as
$$
|black \rangle \langle black|
$$
or in matrix representation (if a black state is {{1},{0},{0}})
black = ( {
{1, 0, 0},
{0, 0, 0},
{0, 0, 0}
} );
To match any state one can ise an Identity operator
indifferent = IdentityMatrix[3];
Then by peforming this operation for each cell in a pattern and multiplying the result, we can estimate the probablity of a match using Boltzman Function
probabillities = (
Exp[-Times @@
Flatten@
Table[
Transpose[part[[i, j]]] . #[[i, j]] . part[[i, j]], {i, 1,
kernelsize}, {j, 1, kernelsize}] / T ]
) & /@ allpatterns;
Using RandomChoice we can choose the desired rule
Pattern replacing
The analogy for this process I would call a transition between states, that is performed individually for each cell in the pattern. Any transisitons can also be performed using operators. For example using
toblack = ( {
{1, 0, 0},
{1, 0, 0},
{1, 0, 0}
} );
on an arbitary state toblack.state we can switch it to black (for sure). Of course it is not even close to the Hermitian operator.
Of it is turned out that the probablity for our pattern is extremely low, then, we might ommit the transition applying Identity operator
nochange = IdentityMatrix[3];
Then, we can apply our operators and make transiitons for the corresponding matched cells and update the board.
Implementation
Basis definition and etc
(* set basis wavefunctions and operators *)
DefineBasis[x_List] :=
Module[{initial, dim, transitions, eoperators, basis, replacements, colors},
dim = Length[x];
basis = Table[Table[{KroneckerDelta[i, j]}, {i, dim}], {j, dim}];
replacements = Association[(Rule @@ #) & /@ Transpose[{x, basis}]];
(* define operators for energies (pattern matching) *)
eoperators =
Join[Table[basis[[i]] . Transpose[basis[[i]]], {i, dim}], {IdentityMatrix[dim]}];
(* define operators fro transitions (pattern replacement) *)
transitions = Table[
Table[basis[[j]] . Transpose[basis[[i]]], {i, dim}] // Total, {j, dim}];
colors =
ColorData[22, "ColorList"][[#]] & /@
Range[1, dim] /. {RGBColor -> List};
<|"dim" -> dim, "basis" -> basis, "colors" -> colors,
"replacements" -> replacements,
"operators" -> <|"energy" -> eoperators,
"transition" -> transitions|>|>
];
(* convert board to vector field *)
Board2States[board_, basis_] := Module[{},
Table[
basis["replacements"][board[[i, j]]], {i, Length[board[[1]]]}, {j,
1, Length[board]}]
];
(* helper function to generalize 1D rule to 2D case *)
RulePermute[x_List, indef_, rules_] :=
Module[{array, dump = {}},
array = Table[indef, {i, Length[x]}, {j, Length[x]}];
array[[1 ;; , 1]] = rules /@ x;
Do[AppendTo[dump, array];
array = Reverse@(Transpose@array);, {i, 1, 4}];
dump];
(* rule -> set of operators generator *)
MakeRule[rule_, basis_, strength_ : -1.01, stationary_ : -1] :=
Module[{rules = rule, strengths = strength, interaction, evolve, diff, listrules, i = 0},
If[Head[rule] === Rule, rules = {rule}; ];
If[Head[strength] =!= List, strengths = Table[strength, {j, 1, Length[rules]}]];
listrules = Table[
i = i + 1;
(* forms a set of operators for detecting a pattern *)
interaction =
strengths[[i]] RulePermute[r[[1]],
basis["operators", "energy"] // Last,
basis["operators",
"energy"][[Position[basis["replacements"] // Keys, #] //
First // First]] &];
(* forms a set of operators for the replacing a pattern *)
(* find a diffenrce between left and right part of a rule and decide which operator to use *)
evolve =
If[#[[1]] === #[[2]], basis["operators", "energy"] // Last,
basis["operators", "transition"][[
Position[basis["replacements"] // Keys, #[[2]]] // First //
First ]]] & /@ Transpose[{r[[1]], r[[2]]}];
evolve =
RulePermute[evolve, basis["operators", "energy"] // Last,
Identity];
Transpose[{interaction, evolve}]
, {r, rules}];
(* if a multiple rules provided, then join them *)
listrules = Join @@ listrules;
(* just some output format *)
<|"perturbation" -> (<|"interaction" -> #[[1]],
"evolution" -> #[[2]]|> & /@ listrules),
"stationary" -> <|
"interaction" ->
stationary Table[
basis["operators", "energy"] // Last, {i,
Length[evolve // First]}, {j, Length[evolve // First]}],
"evolution" -> None|>
|>
]
(* helper to convert assotications to a list of lists for speeding up *)
RulesPreProcess[rules_] := With[{zero = rules["stationary"] },
Join[{#["interaction"], #["evolution"]} & /@
rules["perturbation"], {{zero["interaction"],
zero["evolution"]}}] // Transpose
];
SetAttributes[LocalUpdate, HoldFirst]
(* update function: match and replace *)
LocalUpdate[field_, rules_, T_] :=
Module[{convolved, choice, poses, kernelsize, pos, part, x, y, probabillities, rule, operation},
kernelsize = rules // First // First // First // Length;
(* flip a coin and take a random position at the field *)
{x, y} = List[RandomInteger[{1, Length[field[[1]]] - kernelsize + 1}],
RandomInteger[{1, Length[field] - kernelsize + 1}]];
(* take a part of the field with a size of a pattern *)
part = field[[x ;; x + kernelsize - 1, y ;; y + kernelsize - 1]];
(* convolution, i.e. apply energy operators form all rules and calcululate probabillity *)
probabillities = (
Exp[-Times @@
Flatten@
Table[
(* simple matrix multiplication for each cell *)
Transpose[part[[i, j]]] . #[[i, j]] . part[[i, j]], {i, 1,
kernelsize}, {j, 1, kernelsize}]/T]
) & /@ rules[[1]];
(* pick up a random operation to evolve (or to make a transition to a new set of states) *)
operation =
rules[[2,
RandomChoice[
probabillities -> Range[1, Length[rules // First]]]]];
(* if the probabillity was low and it ended up in stationary state *)
If[operation === None, Return[]];
(* perform the transition by applying operators to the states *)
field[[x ;; x + kernelsize - 1, y ;; y + kernelsize - 1]] =
Table[
operation[[i, j]] . part[[i, j]], {i, 1, kernelsize}, {j, 1,
kernelsize}];
]
(* helper function for coloring *)
GetColor[basis_][{{f_}}] := With[{l = f // Flatten}, (l basis["colors"]) // Total];
Now we are ready to see the show.
Maze Generator
Let us try rule 100->120 - Maze generator.
basis = DefineBasis[{0, 1, 2}];
rules = MakeRule[{1, 0, 0} -> {1, 2, 1}, basis, -1.01, -1];
prules = RulesPreProcess[rules];
board = Table[0, {i, 25}, {j, 25}];
(* spawn a point *)
board[[10, 10]] = 1;
field = Board2States[board, basis];
dfield = field;
Dynamic[BlockMap[GetColor[basis], dfield, {1, 1}] // Raster // Graphics, TrackedSymbols :> {dfield}]
Do[
Do[LocalUpdate[field, prules, 0.01], {i, 1, 10 Length[board]}];
dfield = field;
Pause[0.05];
, {j, 1, 100}]
The result
As one can see there are some extra parameters provided
MakeRule[{1, 0, 0} -> {1, 2, 1}, basis, -1.01, -1]
The values -1.01, -1 stand for the multiplers of energy operators for our rule and for the stationary state, respectively. Our rule has to be bring some energy gain to the system, otherwise there will be no reason for the transistion (unless the temperature of the system is high). -1 as you can see is also negative, this we need to keep the system in the stationary state if there is no matches were found.
Imagine we have -1.01 and 0 as a parameters. Then, if there is no matches our interaction energy is 0, but the energy of the stationary state of the system is also 0! Therefore the probabillity that the rule will be applied is 0.5 no matter what the temperature we have.
Thus we need to make sure, that stationary state is stable, but the rule must have a slightly more stable if we have a match.
To control the temperature one change change the variable in LocalUpdate[field, prules, temperature]. Higher values will bring more fluctiations into the system, that will make sometimes a random transisitons, even if our rule has no matches. Which is also interesting.
Here a couple of examples
Random flip
The simples rule 0->1
basis = DefineBasis[{0, 1}];
rules = MakeRule[{0} -> {1}, basis, -1.01, -1];
Growth Model
Rule 10->11, so we need nucleation or just wait longer (when the temperature is above 0)
basis = DefineBasis[{0, 1}];
rules = MakeRule[{1, 0} -> {1, 1}, basis, -1.01, -1];
Maze-backtracker
basis = DefineBasis[{0, 1, 2, 3}];
rules = MakeRule[{
{1, 0, 0} -> {2, 2, 1},
(* to help unstuck *)
{1, 2, 2} -> {3, 3, 1},
{0, 1, 2} -> {0, 3, 1}
}, basis, {-2.0, -1, -0.9}];
Do[
Do[LocalUpdate[field, prules, 0.03], {i, 1, 10 Length[board]}];
dfield = field;
Pause[0.1];
, {j, 1, 500}]
Comments
There are still a lot issues with bondary conditions. Now it is has a dead region with a thickness of 1 block.
I had a lot of fun and frustration making this. ;)
UPD
- Added support for 2D rules and multiple rules sizes.
- A wildcard Null - applies identity operator (match any)
- new syntax for rules
{Rule[], energy strength}, ...
Updated code (sorry, no comments)
DefineBasis[x_List] :=
Module[{initial, dim, transitions, eoperators, basis, replacements,
colors},
dim = Length[x];
basis = Table[Table[{KroneckerDelta[i, j]}, {i, dim}], {j, dim}];
replacements =
Association[(Rule @@ #) & /@ Transpose[{x, basis}]];
replacements[Null] = Null;
eoperators =
Join[
Table[
basis[[i]] . Transpose[basis[[i]]], {i, dim}], {IdentityMatrix[
dim]}];
transitions = Table[
Table[basis[[j]] . Transpose[basis[[i]]], {i, dim}] // Total
, {j, dim}];
transitions = Join[transitions, {IdentityMatrix[dim]}];
colors =
ColorData[22, "ColorList"][[#]] & /@
Range[1, dim] /. {RGBColor -> List};
<|"dim" -> dim, "basis" -> basis, "colors" -> colors,
"replacements" -> replacements,
"operators" -> <|"energy" -> eoperators,
"transition" -> transitions|>|>
];
Board2States[board_, basis_] := Module[{},
Table[
basis["replacements"][board[[i, j]]], {i, Length[board[[1]]]}, {j,
1, Length[board]}]
];
GetColor[basis_][{{f_}}] :=
With[{l = f // Flatten}, (l basis["colors"]) // Total];
RulePermute[x_List, indef_, rules_] := Module[{array, dump = {}},
array = Table[indef, {i, Length[x]}, {j, Length[x]}];
array[[1 ;; , 1]] = rules /@ x;
Do[AppendTo[dump, array];
array = Reverse@(Transpose@array);, {i, 1, 4}];
dump];
MakeRule[rule_, basis_, stationary_ : -1] :=
Module[{rules = rule, interaction, evolve, diff, listrules, i = 0},
If[Head[rule] === Rule, rules = {rule}; ];
listrules = Table[
If[Depth[r[[1, 1]]] === 3,
interaction = {Function[
y, (basis["operators",
"energy"][[Position[basis["replacements"] // Keys, #] //
First // First]] &) /@ y] /@ Reverse[r[[1, 1]]]};
evolve = {Function[
y, (basis["operators",
"transition"][[Position[
basis["replacements"] // Keys, #] // First //
First]] &) /@ y] /@ Reverse[r[[1, 2]]]};
,
interaction =
RulePermute[r[[1, 1]], basis["operators", "energy"] // Last,
basis["operators",
"energy"][[Position[basis["replacements"] // Keys, #] //
First // First]] &];
evolve =
If[#[[1]] === #[[2]], basis["operators", "energy"] // Last,
basis["operators", "transition"][[
Position[basis["replacements"] // Keys, #[[2]]] // First //
First ]]] & /@ Transpose[{r[[1, 1]], r[[1, 2]]}];
evolve =
RulePermute[evolve, basis["operators", "energy"] // Last,
Identity];
];
Transpose[{interaction, evolve,
Table[r[[2]], {i, Length[interaction]}]}]
, {r, rules}];
listrules = Join @@ listrules;
<|"perturbation" -> (<|"interaction" -> #[[1]],
"evolution" -> #[[2]], "strength" -> #[[3]]|> & /@
listrules),
"stationary" -> <|
"interaction" ->
Table[
basis["operators", "energy"] // Last, {i,
Length[evolve // First]}, {j, Length[evolve // First]}],
"evolution" -> None, "strength" -> stationary|>
|>
]
RulesPreProcess[rules_] := With[{zero = rules["stationary"] },
Join[{#["interaction"], #["evolution"], #["strength"]} & /@
rules["perturbation"], {{zero["interaction"], zero["evolution"],
zero["strength"]}}] // Transpose
];
SetAttributes[LocalUpdate, HoldFirst]
LocalUpdate[field_, rules_, T_] :=
Module[{convolved, choice, poses, kernelsize, pos, part, x, y,
probabillities, rule, operation},
kernelsize = (# // First // Length) & /@ (rules // First);
kernelsize = Max[kernelsize];
{x, y} =
List[RandomInteger[{1, Length[field[[1]]] - kernelsize + 1}],
RandomInteger[{1, Length[field] - kernelsize + 1}]];
part = field[[x ;; x + kernelsize - 1, y ;; y + kernelsize - 1]];
probabillities = MapIndexed[
With[{localsize = Length[#1], r = #1, amp = rules[[3, #2[[1]]]]},
Exp[ -amp Times @@
Flatten@
Table[
Transpose[part[[i, j]]] . r[[i, j]] . part[[i, j]], {i, 1,
localsize}, {j, 1, localsize}]/(T)]
] &, rules[[1]]];
operation =
rules[[2,
RandomChoice[probabillities -> Range[1, Length[rules // First]]]]];
If[operation === None, Return[]];
With[{localsize = Length[operation]},
field[[x ;; x + localsize - 1, y ;; y + localsize - 1]] =
Table[
operation[[i, j]] . part[[i, j]], {i, 1, localsize}, {j, 1,
localsize}]
];
]
And now finally I can recreate a beautiful example with growing flowers, which is a combination of 1D and 2D rules
basis = DefineBasis[{-3, -2, -1, 0, 1, 2, 3, 4}];
rules = MakeRule[{
(* making ground *)
{{-2, -3} -> {-2, -2}, -5},
(* making sky/air *)
{{0, -3} -> {0, 0}, -5},
(* making a soil at the interface ground/sky*)
{{0, -2} -> {0, 1}, -2},
(* plant a seed of a flower in the soil *)
{( {
{0, 0, 0, 0},
{0, 0, 0, 0},
{0, 0, 0, 0},
{1, 1, 1, 1}
} ) -> ( {
{0, 0, 3, 0},
{0, 0, 2, 0},
{0, 0, 2, 0},
{1, 1, 1, 1}
} ), -0.999},
(* grow a plant *)
{{
{0, 0, 0},
{0, 0, 0},
{0, 3, 0}
} -> {
{0, 0, 0},
{0, 3, 0},
{0, 2, 0}
}, -2.5},
(* moves away from the straight line *)
{( {
{0, 0, 0, },
{0, 0, 0, },
{3, 0, 0, },
{ , , 0, }
} ) -> ( {
{0, 3, 0, },
{0, 2, 0, },
{2, 2, 0, },
{ , , 0, }
} ), -2},
(* mirrored version *)
{{{Null, 0, 0, 0}, {Null, 0, 0, 0}, {Null, 0, 0, 3}, {Null, 0,
Null, Null}} -> {{Null, 0, 3, 0}, {Null, 0, 2, 0}, {Null,0,
2, 2}, {Null, 0, Null, Null}}, -2},
(* making branches *)
{( {
{0, 0, 0, 0},
{0, 0, 0, 0},
{0, 0, 3, 0},
{0, , , }
} ) -> ( {
{0, 3, 0, 3},
{0, 2, 0, 2},
{0, 2, 2, 2},
{0, , , }
} ), -2},
(* head of the flower *)
{{
{0, 0, 0},
{0, 3, 0},
{0, 2, 0}
} -> {
{0, 4, 0},
{4, 3, 4},
{0, 4, 0}
}, -1}
}, basis];
prules = RulesPreProcess[rules] // N;
board = Table[-3, {i, 30}, {j, 30}];
board[[10, 10]] = 0;
board[[3, 10]] = -2;
field = Board2States[board, basis] // N;
dfield = field;
basis = basis // N;
Dynamic[BlockMap[GetColor[basis], dfield, {1, 1}] // Raster //
Graphics, TrackedSymbols :> {dfield}]
Do[
Do[LocalUpdate[field, prules, 0.003], {i, 1, Length[board]}];
dfield = field;
Pause[0.1];
, {j, 1, 500}] // Quiet
Here we go
ReplaceAllon it directly? Putting aside the Markov randomness, do these replacement rules work individually for your purposes? $\endgroup${1,0}->{1,1}won't work. I was wondering if it is possible to define the field at a combination of downvalues of symbols or a graph, but did not manage to figure it out. $\endgroup$