First evaluate in place all subexpressions matching pattern before substituting result

Question

ReplaceAll[expr, {patt1:>rhs1, patt2:>rhs2, …}] works by looking for and making all possible replacements of subexpressions in expr with rhs1,rhs2… that match patt1,patt2…. Then, the main evaluator evaluates all the rhs1, rhs2… appearing in the modified expr.

This scheme of replacement works well especially if rhs1,rhs2… cause side effects (like i++). However, if rhs1,rhs2… are computationally intensive functions without side effects, and expr is a large expression containing many instances matching patt1,patt2…, the main evaluator ends up evaluating the big functions rhs1,rhs2… many times.

The canonical solution to this problem in Mathematica is memoization where rhs is supposed to be defined like rhs := rhs = ... so that even though many instances of rhs1,rhs2… appear in the modified expr, the main evaluator will recognize precomputed values based on an ever growing list of DownValues.

However, I believe another possible solution to this problem, especially suitable if rhs1,rhs2… are not recursive functions, is to follow a different scheme: (1) scout out and note the positions of subexpressions in expr uniquely matching the patt1, patt2, (2) evaluate in place the corresponding rhs1,rhs2… exactly once for each unique match, and finally (3) substitute all results into expr at their noted positions. I believe that apart from saving some memory, this scheme has the advantage of (a) not rebuilding a hash table of memoized DownValues, (b) not looking for a possible match in the hash table, and (c) not scanning the entire expr each time rhs is evaluated in the modified expr.

So my question is: is there a built-in replacement function that implements the replacement scheme described above? I would find it useful. Otherwise, is there a method to write a top-level implementation of the scheme?

Here is a first-attempt at an implementation:

ClearAll[ReplaceAllOnce];
ReplaceAllOnce[expr_, rules_] := 
  ReplacePart[expr, (Position[expr, #, {0, Infinity}, Heads -> True] -> (Replace[#, rules])) & /@ 
    DeleteDuplicates[Cases[expr, Alternatives @@ rules[[All, 1]], {0, Infinity}, Heads -> True]]]

However, I don't like it much because it has to traverse through expr multiple times.

---

As requested, here is an example. Suppose I want to apply to this expression,

expr := (a parity[1] + b parity[2] + c parity[3])*
            (d parity[1] + e parity[2] + f parity[3]);

The following set of rules

rules = {parity[x_Integer?OddQ] :> (Pause[1]; 1), 
         parity[x_Integer?EvenQ] :> (Pause[1]; 2)};

A straightforward ReplaceAll takes 6 seconds (because there are 6 instances of parity).

ReplaceAll[expr, rules] // AbsoluteTiming
(*  {6.00317, (a + 2 b + c) (d + 2 e + f)}  *)

Then using the search-and-replace once method using ReplaceAllOnce proposed above takes only 3 seconds, and uses 5680 bytes of memory according to MemoryInUse[]:

ReplaceAllOnce[expr, rules] // AbsoluteTiming
(*  {3.00413, (a + 2 b + c) (d + 2 e + f)}  *)

The timing is similar to using memoization, but this uses 12800 bytes (larger, to store the memoized cases):

Clear[parity];
parity[x_Integer?OddQ]  := parity[x] = (Pause[1]; 1);
parity[x_Integer?EvenQ] := parity[x] = (Pause[1]; 2);

expr // AbsoluteTiming
(* {3.00486, (a + 2 b + c) (d + 2 e + f)} *)

Answer 1

First, let me discourage you from solutions using Position to find the places where the rules should be applied. At least, if you want to reproduce the behavior of ReplaceAll: ReplaceAll goes from the trunk to the leaves. Once a substitution is made, no more substitutions are done in the parts.

Consider for example

a + b f[x] /. {b y_f :> r[1], f[y_] :> r[2]}
(*a + r[1]*)

If the substitution started at the leaves, the result should have been a + b r[2]. That means that code similar to

Position[a + b f[x], Alternatives[f[y_], b y_f]]
(*{{2, 2}, {2}}*)

would have to check for nested positions and keep those closer to the trunk. That is not trivial and better left to be done by Mathematica.

Since we need to find the places for replacement and we cannot use Position, we should use ReplaceAll. Here, the answer by Coolwater almost hits the right spot. If one checks carefully, that answer traverses expr twice.

Memoization should be your best approach, but your comment on Carl Woll's answer leads me to believe that you might be worried about polluting the namespace. [Edit: Actually, I do not understand why you have the constraint against memoization. See appended paragraphs]

The following code is not what I consider the best approach:

ReplaceAllOnceH[expr_, rules_List] := 
 With[{t = System`Utilities`HashTable[]}, 
  ReplaceAll[expr, 
   Replace[rules, 
    HoldPattern[
      a_ :> b_] :> (Pattern[vv$, a] :> 
       If[System`Utilities`HashTableContainsQ[t, vv$], 
        System`Utilities`HashTableGet[t, vv$], 
        With[{v = b}, System`Utilities`HashTableAdd[t, vv$, v]; 
         v]]), {1}]]]

It is basically a poor man's memoization. However, let us review how it stands against your claims on a method based on Position:

(a) not rebuilding a hash table of memoized DownValues, (b) not looking for a possible match in the hash table, and (c) not scanning the entire expr each time rhs is evaluated in the modified expr.

(a) Rebuilding a hash table is expensive only when going over the size of the bucket. (b) Looking up in a hash table is an O(1) operation (on average). (c) Done by ReplaceAll.

---

PS. In my computer, ReplaceAllOnceH is barely faster than MemoizedReplaceAll in some cases and barely slower in others.

---

Could you elaborate on your requirement against memoization?

If it is about memory consumption, you could use Clear to remove the memoized values. Or you could use something similar to the solution above where the hash table is destroyed automatically at the end of the replacement.

If it is speed, unless you have a set of expressions where you know that Position will not lead to nested replacements, then you would have to check for nested position. Otherwise, you could be spending resources on computing a replacement that you will not use. For example,

(1 + a b + (a + (1 + a b) c) d)/(b + (1 + b c) d)
Position[%, 1 + _]
(*the above yields 3 positions*)
ReplacePart[%%, Thread[% -> Extract[%%, %, Framed]]]
(*there are only two boxes in the result*)

Having only two boxes in the result emulates the behavior of ReplaceAll but there were 3 calculations done.

Now, manipulating the results of Position to remove the nested replacements is an operation that grows like the square of the number of replacements. Memoization grows linearly.

Answer 2

In case you do have some knowledge about the expressions to be replaced and you can confidently predict that Position will not lead to nested replacements, you might still want to use a hash table to hold the intermediate results.

ReplaceAllOnceH2[expr_, rules_List] := 
 With[{p = Position[expr, Alternatives @@ rules[[All, 1]]], 
   t = System`Utilities`HashTable[]}, 
  ReplacePart[expr, 
   Thread[p -> 
     Extract[expr, p, 
      If[System`Utilities`HashTableContainsQ[t, #], 
        System`Utilities`HashTableGet[t, #], 
        With[{w = Replace[#, rules]}, 
         System`Utilities`HashTableAdd[t, #, w]; w]] &]]]]

Answer 3

Edit: Carl Woll's answer made me realize it should be done differently. The code below Sows the replacement values (wrapped by eval and Hold), because these may be non-injective wrt. the variables in the patterns.

eval prevents the second ReplaceAll from replacing parts that didn't match the pattern in the first place. Hold prevents evaluation until the association is constructed (whereby duplicates are deleted). First is applied to release Hold.

replaceOnce[expr_, rules_] := Module[{eval, expr2, reap},
 {expr2, {reap}} = Reap[expr /. Fold[MapAt[#2, #, {All, 2}] &, rules, {Hold, eval, Sow}]];
 expr2 /. First @@@ AssociationThread[reap -> reap]]

replaceOnce[expr, rules] // AbsoluteTiming

{2.0040347, (a + 2 b + c) (d + 2 e + f)}

Old answer:

You can use ReplaceAll to find the cases and Sow them as non-delayed rules. By putting the non-delayed rules in an Association, duplicates are automatically deleted. The association holds the evaluation, however, so Map[Identity] is applied.

Even if your rules is short, there may be many distinct matches, and if that is case replacing with an Association is usually preferable.

evaluateRules[expr_, rules_List] :=
     Module[{eval, rules2 = rules},
        rules2[[All, 1]] = eval /@ rules2[[All, 1]];
        Identity /@ (Association[Reap[expr /. #;][[2, 1]]] /. rules2) &[
            Pattern[a, #] :> Sow[a -> eval[a]] & /@ rules[[All, 1]]]]

expr /. evaluateRules[expr, rules] // AbsoluteTiming

{3.0014318, (a + 2 b + c) (d + 2 e + f)}

Answer 4

1.05174 seconds

This is very simple, first get the Position, then evaluate only the rules that are been used (non empty Position list ) using ParallelTable to avoid sequential evaluation. Using Table would take 2 seconds.

expr := (a parity[1] + b parity[2] + c parity[3])*(d parity[1] + 
     e parity[2] + f parity[3]);

rules = {
   parity[x_Integer?OddQ] :> (Pause[1]; 1),
   parity[x_Integer?EvenQ] :> (Pause[1]; 2),
   nonexistent[_] :> (Pause[1]; 0)
   };

LaunchKernels[]
ClearAll[replacerh];
replacerh[expr_, rules_] := With[
  {pos = Position[expr, #] & /@ rules[[All, 1]]},
  ReplacePart[expr,
   ParallelTable[
    If[
     Length[pos[[k]]] > 0,
     pos[[k]] -> rules[[k, 2]],
     Nothing
     ]
    , {k, Length[pos]}
    ]
   ]
  ]

replacerh[expr, rules] // AbsoluteTiming
(* {1.05174, (a + 2 b + c) (d + 2 e + f)} *)

Answer 5

Here is an implementation of the memoization approach:

SetAttributes[memoize,HoldAll];

memoize[a_] := memoize[a] = a

MemoizedReplaceAll[expr_, lhs_:>rhs_] := Internal`InheritedBlock[{memoize},
    ReplaceAll[expr, lhs :> memoize[rhs]]
]

MemoizedReplaceAll[expr_,rules_List] := Internal`InheritedBlock[{memoize},
    With[{new=Replace[rules, RuleDelayed[lhs_, rhs_] :> lhs:>memoize[rhs], {1}]},
        ReplaceAll[expr, new]
    ]
]

An example:

SeedRandom[1]
data = RandomInteger[10, {10, 10}]

{{1, 4, 0, 7, 0, 0, 8, 6, 0, 4}, {1, 8, 5, 1, 1, 1, 3, 2, 10, 1}, {6, 0, 2, 6, 4, 5, 4, 3, 0, 1}, {3, 5, 3, 0, 3, 2, 3, 9, 5, 1}, {5, 2, 3, 9, 1, 0, 4, 4, 1, 5}, {2, 7, 9, 9, 8, 10, 0, 10, 10, 7}, {4, 9, 2, 6, 3, 2, 1, 1, 6, 1}, {1, 6, 8, 6, 5, 6, 0, 10, 7, 9}, {1, 4, 4, 10, 3, 5, 2, 3, 1, 2}, {5, 10, 8, 3, 6, 3, 10, 3, 10, 3}}

Two functions that are slow:

f[x_] := (Pause[.5]; x^2)
g[x_] := (Pause[.5]; x^3)

Using ReplaceAll vs MemoizedReplaceAll:

r1 = ReplaceAll[data, x_?PrimeQ :> f[x]]; //AbsoluteTiming
r2 = MemoizedReplaceAll[data, x_?PrimeQ :> f[x]]; //AbsoluteTiming

r1===r2

{18.0379, Null}

{2.00427, Null}

True