Dynamic Programming with delayed evaluation

Question

By using dynamical programming, we can save intermediate steps for recursive relations, as in

f[n_]:= f[n] = f[n-1] + f[n-2]

However, this only stores definitions for explicit values of n. But I need to be able to store function definitions, where the number of parameters of the function depends on each step. More specifically, I have something like $$ t^{a_1 \cdots a_n} = A \, t^{a_1 \cdots a_{n-2}} + B \, t^{a_2 \cdots a_n} \\ t^a = C^a \\ t^{ab} = D^{ab} $$

Where all $a_i$ are indices with unspecified values. This recursive definition allows me to express any $t^{a_1 \cdots a_n}$ as a combination of $C$'s and $D$'s. It's easy to implement this in mathematica:

t[a_]:= c[a]
t[a_,b_]:= d[a,b]
t[a__]:= A t@@Most[Most[{a}]] + B t@@Rest[{a}]

but for high $n$, calculation times become very long (and especially because my function is a bit more complicated than this). This is where we normally would use dynamic programming, something like:

t[a__]:= t[a] = A t@@Most[Most[{a}]] + B t@@Rest[{a}]

but this is not what we need, because if I execute

t[a,b,c]

it only makes a definition for a, b and c literally, meaning that t[b,a,c] will have to be recalculated recursively. This is because it is stored as Set, and not a SetDelayed:

t[a,b,c] = A c[a] + B d[b,c]

What I would need is some kind of magic dynamic programming that stores function definitions as in

t[a__]:= Magic.... =A t@@Most[Most[{a}]] + B t@@Rest[{a}]

such that when I execute

t[a,b,c]

it stores

t[a_,b_,c_] := A c[a] + B d[b,c]

etcetera.

I've looked into Leonid's answer for this question, but I can't easily see how to adapt it for a growing number of arguments.

Any ideas? Many thanks in advance!

EDIT

Thanks to Leonid's answer, I now understand a bit how this can be solved. Unfortunately, the way I'd need it, is a bit more complicated. I have several function definitions, with Condition, PatternTest, Optional and next to the recursive ('growing') variables I also have some that aren't.

At first I thought I'd figure it out myself based on Leonid's answer (hence the accept), but I didn't. The idea is to make a wrapper memoize, which I wrap around my function definition such that it memorises the function definitions. An example in pseudocode:

memoize[
    t[dot[a__], b_?EvenQ, c_:1]/;b!=c := .... some recursive function of t ...
]

The first thing I tried is to retrieve all patternames, and replace them with Unique[] ones:

ClearAll@memoize
SetAttributes[memoize, HoldAllComplete]
memoize[expr_SetDelayed] :=
    (* first we retrieve the lhs and rhs *)
    With[{funcLHS = Hold[expr] /. Hold[SetDelayed[x_, y_]] :> Hold[x] , 
          funcRHS = Hold[expr] /. Hold[SetDelayed[x_, y_]] :> Hold[y]},
        (* next we retrieve the names of the patterns used *)
        With[{patternNames = First /@ Cases[funcLHS, _Pattern, Infinity]},
            (* we make some locally scoped patterns *)
            With[{varsExt = Table[Unique[], {Length[patternNames]}]},
                (* and we express the lhs and rhs in function of the former locally scoped patterns *)
                With[{lhs = funcLHS /. Rule @@@ Transpose[{patternNames, varsExt}], 
                      rhs = funcRHS /. Rule @@@ Transpose[{patternNames, varsExt}]},

                   {lhs,rhs}/.{Hold[x_],Hold[y_]}:>SetDelayed[x,y]
    ]]]]

And this works, but it doesn't do anything useful. Now I would need to replace the central part with Leonid's answer somehow, which means I have to inject vars into lhs and rhs, with the additional difficulty that Length[vars] won't equal Length[varsExt]. So this is what I tried:

ClearAll@memoize
SetAttributes[memoize, HoldAllComplete]
memoize[expr_SetDelayed] :=
    With[{funcLHS = Hold[expr] /. Hold[SetDelayed[x_, y_]] :> Hold[x] , 
          funcRHS = Hold[expr] /. Hold[SetDelayed[x_, y_]] :> Hold[y]},
        With[{patternNames = First /@ Cases[funcLHS, _Pattern, Infinity]},
            With[{varsExt = Table[Unique[], {Length[patternNames]}]},
                With[{lhs = funcLHS /. Rule @@@ Transpose[{patternNames, varsExt}], 
                      rhs = funcRHS /. Rule @@@ Transpose[{patternNames, varsExt}]},

                    {lhs, Hold[
                        With[{vars = Table[Unique[], {Length[{#}]}] & /@ varsExt},
                            With[{pts = (Pattern[#1, _] &) /@ vars},
                                {lhs/.MapThread[ Pattern[#1, _] -> Sequence @@ (Function[x, Pattern[x, _]] /@ #2) & , {varsExt, vars} ]
                                ,
                                rhs/.MapThread[ #1 -> Sequence @@ #2 & , {varsExt, vars} ]
                                }/.{Hold[x_],Hold[y_]}:>Set[x,y] ;
                                lhs/.x_Pattern:>x[[1]]
                        ]]]
                     } /. {Hold[x_], Hold[y_]} :> SetDelayed[x, y] 
    ]]]]

But the pattern replacing on lhs doesn't work, and it doesn't take care of Condition, Optionaland PatternTest. So I'm totally stuck, although I have the feeling that it's only a few changes to make it work..

Answer 1

I will offer a rather cryptic solution using nested version of the injector pattern, but it should be possible to also rewrite it using more readable methods, if really needed.

Solution

Here is the code:

ClearAll[t];    
t[a_]:=  c[a]
t[a_,b_]:= d[a,b]
t[a__]:=
    With[{vars=Table[Unique[],{Length[{a}]}]},
      With[{pts=(Pattern[#1,_]&)/@vars},            
        pts/. 
          {x__}:>
            (
              Most[Most[vars]]/. 
                {arhs___}:>(Rest[vars]/. 
                  {brhs___}:>Set@@Hold[t[x],A t[arhs]+B t[brhs]])
            );
        t[a]
      ]
    ];

We get:

t[a, b, c, d]

(* A d[a, b] + B (A c[b] + B d[c, d])  *)

and now we have the following definitions generated in the process:

?t
Global`t
t[a_]:=c[a]

t[a_,b_]:=d[a,b]

t[$21_,$22_,$23_]=A c[$21]+B d[$22,$23]

t[$17_,$18_,$19_,$20_]=A d[$17,$18]+B (A c[$18]+B d[$19,$20])

t[a__]:=With[{vars=Table[Unique[],{Length[{a}]}]},
   With[{pts=(Pattern[#1,_]&)/@vars},pts/. {x__}:>(Most[Most[vars]]/. {arhs___}:>
      (Rest[vars]/. {brhs___}:>Set@@Hold[t[x],A t[arhs]+B t[brhs]]));t[a]]]

How it works

Basically, dummy symbols were used to programmatically create a pattern for the l.h.s., and then to create the corresponding memoized definition.

In more details, here is what happens: whenever t[x,y,...] is called, the general definition t[a__]:=... is only invoked if there isn't yet a definition for the exact number of input arguments Length[{a}]. In the process of evaluation of t[a], such definition is created, so all subsequent times calls with the same number of arguments will use the new specialized definition, not the general one.

The general idea here can be expressed in the following pseudo-code:

t[a__]:=
  Module[{helper vars},
    t[var1_,...,varn_] = rhs[var1,...varn];
    t[a]
  ]

where we first construct a specialized definition for a fixed number of input arguments, and then use it right away for our specific arguments a.

Now, how it works: let's say we start with the call

t[m,n,k]

so that

Clear[m, n, k];
a = Sequence[m, n, k];

models what a becomes. We first generate a set of dummy variables:

vars = Table[Unique[], {Length[{a}]}]

(* {$3, $4, $5} *)

and the corresponding patterns:

pts = (Pattern[#1, _] &) /@ vars

(* {$3_, $4_, $5_} *)

What we would like now to do is to construct this code:

t[$3_, $4_, $5_] = A * t[$3] + B * t[$4, $5]

In principle, we could just try to use

Evaluate[t @@ pts] = A t @@ Most[Most[vars]] + B t @@ Rest[vars]

to create the definition. However, we must be more careful. Direct evaluation like Evaluate[t @@ pts] will lead to the general definition t[a__] applied now for t[$3_, $4_, $5_] (since patterns are also arguments), and this is something we don't want. In other words, we want to prevent unwanted evaluation during assignment.

This is where injector pattern comes handy. Nested injector pattern is used to collect more than one part of an expression within a single nested rule, and is actually an overkill here. The definition can be constructed with this simpler code:

pts /. {x__} :> Set @@ Hold[t[x], A t @@ Most[Most[vars]] + B t @@ Rest[vars]]; 

The construct pts /. {x__} :> ... is used to inject x (in this case, a sequence of paterns $3_, $4_, $5_ ) into the right-hand side before thr r.h.s. is allowed to evaluate. Here we use it to programmatically construct t[$3_, $4_, $5_] - the l.h.s. of the new definition. The construct Set@@Hold[lhs, rhs] is used to fool the built-in variable renaming mechanism, which would otherwise rename our variables in an attempt to protect them from name collisions (which is exactly what we don't want here).

So, the final code can take a simpler form and look like

t[a__]:=
    With[{vars=Table[Unique[],{Length[{a}]}]},
        With[{pts=(Pattern[#1,_]&)/@vars},                
            pts/. 
               {x__}:>Set@@Hold[t[x],A * t @@ Most[Most[vars]]+B * t @@ Rest[vars]];
            t[a]
        ]
    ];

where, as explained, the first line of code in the inner With creates a specialized definition, and then this definition is immediately used on the actual arguments a. All subsequent times, when called with the same number of arguments as in a, t will use the specialized defintion automatically, and not the general t[a__].

Answer 2

You can also try this code. It uses only explicit variables a1,a2,a3... in the definitions. The auxiliary function ComputeDefinition recursively constructs f for given length n of the sequence.

Clear[f,ComputeDefinition]

f[{a1_}] = c[a1];
f[{a1_, a2_}] = d[a1, a2];
ComputeDefinition[1] = 1;
ComputeDefinition[2] = 1;
ComputeDefinition[10] = 1;

ComputeDefinition[n_] := (
   f[Table[ToExpression["a" <> ToString[i] <> "_"], {i, 1, n}]] =
   A f[Table[ToExpression["a" <> ToString[i]], {i, 1, n - 2}]] +
   B f[Table[ToExpression["a" <> ToString[i]], {i, 2, n}]];
   ComputeDefinition[n + 1])

ComputeDefinition[3]
(* 1 *)

?f
(*
Global`f

f[{a1_}]=c[a1]

f[{a1_,a2_}]=d[a1,a2]

f[{a1_,a2_,a3_}]=A c[a1]+B d[a2,a3]

f[{a1_,a2_,a3_,a4_}]=A d[a1,a2]+B (A c[a2]+B d[a3,a4])

f[{a1_,a2_,a3_,a4_,a5_}]=A (A c[a1]+B d[a2,a3])+B (A d[a2,a3]+B (A c[a3]+B d[a4,a5]))
........
........
........
*)


f[{t1, t2, t3, t4, t5, t6}]

(*
A (A d[t1, t2] + B (A c[t2] + B d[t3, t4])) + B (A (A c[t2] + B d[t3, t4]) +
B (A d[t3, t4] + B (A c[t4] + B d[t5, t6])))
*)