Replacing 2D Patterns in 2D Lists | Markov Programming Language

Question

Recently I came across with a probabilistic programming language where programs are combinations of rewrite rules - Markov Junior

I think it is quite interesting to try out implementing some features of this language using Mathematica, since it has an extremely powerful pattern matching engine

The first thing is to perform 1D/2D/3D patterns replacement.

For example, we have a rule

rule = {1,0} -> {1,1}

and we want to apply it to the 2D field of 1,0s... randomly

field = {
 {0, 0, 0, 0, 0, 0, 0},
 {0, 0, 0, 0, 1, 0, 0},
 {0, 0, 0, 0, 0, 0, 0}
}

now we have 4 possibilities for the resulting field

 {0, 0, 0, 0, 0, 0, 0},
 {0, 0, 0, 0, 1, 1, 0},
 {0, 0, 0, 0, 0, 0, 0}

 {0, 0, 0, 0, 1, 0, 0},
 {0, 0, 0, 0, 1, 0, 0},
 {0, 0, 0, 0, 0, 0, 0}

 {0, 0, 0, 0, 0, 0, 0},
 {0, 0, 0, 1, 1, 0, 0},
 {0, 0, 0, 0, 0, 0, 0}

 {0, 0, 0, 0, 0, 0, 0},
 {0, 0, 0, 0, 1, 0, 0},
 {0, 0, 0, 0, 1, 0, 0}

According to the idea of Markov language, one need to pick up a random variant and repeat the pattern matching with it.

The corresponding routine looks like by steps

1. Find all matches with a pattern
2. Pick up a random match
3. Replace it according to the rule
4. repeat

Possible implementation

1. This is a following up question to this one

Since our rule is 1D, but the field array is 2D, then, I decided to expand it

RulePermute[x_List] := Module[{array, dump = {}},
  array = Table[_, {i, 2Length[x]-1}, {j, 2Length[x]-1}];
  array[[Length[x] ;; , Length[x] ]] = x;
  Do[
    AppendTo[dump, array];
    array = Reverse@(Transpose@ array);
  , {i,1,4}];

  dump
];

what it does: it expands our pattern {1,0} to the set of 4

{{_,_,_,}
 {_,1,0},
 {_,_,_}}

{{_,0,_,}
 {_,1,_},
 {_,_,_}}

{{_,_,_,}
 {0,1,_},
 {_,_,_}}

{{_,_,_,}
 {_,1,_},
 {_,0,_}}

and etc... Then we can find all matches using this approach

field = {
 {0, 0, 0, 0, 0, 0, 0},
 {0, 0, 0, 0, 1, 0, 0},
 {0, 0, 0, 0, 0, 0, 0}
};

RulePermute[x_List] := Module[{array, dump = {}},
  array = Table[_, {i, 2Length[x]-1}, {j, 2Length[x]-1}];
  array[[Length[x] ;; , Length[x] ]] = x;
  Do[
    AppendTo[dump, array];
    array = Reverse@(Transpose@ array);
  , {i,1,4}];

  dump
];

(* our pattern to match *)
pattern = RulePermute[{1,0}];
(* our replacement for the matched pattern *)
replacement = RulePermute[{1,1}];

(* sketchy way to generate a function, that checks all 4 variations of the pattern *)
seq = {Table[
  {MatchQS[x, pattern[[r]]], r}
, {r, pattern//Length}], {True, False}} //Flatten;

testfunction = FunctionS[x, Which @@ seq ] /. {FunctionS -> Function, MatchQS -> MatchQ};

(* form an array with matched positions *)
filtered = ArrayFilter[
  testfunction,
  field,
  {{1,1,1},{1,1,1},{1,1,1}}
];

(* extract the coordinates of the matched block in the field array *)
pos = Position[filtered, _?NumberQ]//RandomChoice;

Then, we know the exact position, which is for this case {2,5}.

The actual question

2. I have no idea how to replace it properly utilizing some features of Wolfram Language

In general the desired function should act like ReplaceAll2D (if it was there)

ReplaceAll2D[expr, 

{{_,_,_,}
 {_,1,_},
 {_,0,_}} ->

{{_,_,_,}
 {_,1,_},
 {_,1,_}}
]

I believe there is a way to construct 2D pattern replacing function acting on 2D/3D lists without involving the straightforward procedural programming style

Thanks!

Answer 1

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