139
$\begingroup$

... or are they unnecessary in such a high-level language?

I've been thinking about programming style, coding standards and the like quite a bit lately, the result of my current work on a mixed .Net/Mathematica project.

In object-oriented programming, people talk a lot about "design patterns", essentially, standard ways of approaching particular problems (There's even a conference). The idea has been extended to other domains as well, e.g. for personal productivity.

But some people argue that identified "patterns" are really just a sign that the language in question isn't abstracted enough.

So what, if any, are the patterns of Mathematica programming? Or are they not necessary?

$\endgroup$
5
  • 5
    $\begingroup$ Great question, +1. $\endgroup$ May 10, 2013 at 13:45
  • $\begingroup$ Stating the obvious, but I'm expecting a great answer from Leonid :) $\endgroup$ May 10, 2013 at 13:54
  • $\begingroup$ One thing to mention is that an idiom is a pattern that has gained some acceptance in the community and is commonly used (thus it's somewhat a cultural thing, just like in a natural language). We can identify idioms by looking at what code people typically use. A design pattern can also gain acceptance in the community and become widely known and tested (tested in the sense that it's been proven to be a robust and convenient design by practical usage). But theoretically it's possible to invent completely new design patterns that are not commonly used. The inventor, if s/he posts here, $\endgroup$
    – Szabolcs
    May 12, 2013 at 21:08
  • $\begingroup$ would have to make a case for why his design pattern is effective and why it will work well in a broad range of situation, how it will solve common problems, ect. This is going to be more difficult unless the inventor can present some convincing case studies. My point is that a pattern can be idiomatic (commonly used); it can also be an effective design. The are two independent things. $\endgroup$
    – Szabolcs
    May 12, 2013 at 21:11
  • 3
    $\begingroup$ @JacobAkkerboom I have spent several hours on this, and identified a few dozens of patterns so far. Some of them have been mentioned by others (particularly recent post of @Szabolcs). However, to my mind, a huge amount of work is needed to write posts which really show people how to use those patterns and give enough links, rather than just stating a very brief description (which is what is being currently done for most of them, and most of which is understandable by experts only). As time allows, I will add descriptions for some of them, which may overlap with those already posted by others). $\endgroup$ May 13, 2013 at 13:44

8 Answers 8

96
$\begingroup$

Leaving the general discussion on design patterns for some later time, I will try to collect a comprehensive information on some specific design patterns / idioms, one pattern per post. This time it will be memoization

Basic memoization

What is memoization

Memoization is a technique which enables a program to reuse parts of computations it has done already, to avoid unnecessary repeated work.

Often (although not always) memoization goes hand in hand with recursion. One reason for that is that many recursive functions have very high (often exponential) complexity when implemented straightforwardly, due to the phenomena of overlapping sub-problems. In many such cases, memoization is able to drastically reduce the complexity.

Types of memoization

This is largely a matter of definitions, but, at least in Mathematica, I will distinguish 3 types of memoization:

  • "Simple" memoization. This is basically a trick used to remember a number of past values on which some computationally demanding function has been computing. This is typically done as

    f[x_]:=f[x]= r.h.s.
    

    what is important here is that enabling memoization leads to a speed-up, but not to a drastic change of the complexity of the resulting algorithm. Put in other words, this is a technique of remembering past values for functions without the overlapping sub-problems. A typical example here would be to memoize values of some expensive integral computed on some fixed grid in its parameter space.

  • Algorithmic memoization: this may be technically realized in the same way, but the distinctive feature of this type of problems is that they have overlapping sub-problems, so here memoization leads to a drastic change in algorithm's complexity. Perhaps the simplest well-known example here is Fibonacci numbers:

    fib[0] = fib[1] = 1;
    fib[n_Integer?Positive]:= fib[n] = fib[n-1] + fib[n-2]
    

    where the complexity of the resulting algorithm is reduced from exponential to linear, as a result of memoization.

  • Caching. The main difference between caching and memoization is that memoization is a technique of run-time caching during one and the same function's evaluation (single evaluation process), while general caching is a technique to memorize some computed values for possible future uses, which may be happening for different evaluations. So, we can say that caching is a form of "persistent memoization", whose scope extends to more than one evaluation process.

Why / how it works in Mathematica

The technical reason why Mathematica makes it rather straightforward to implement memoization is that pattern-based definitions are global rules, and more specific rules are applied before more general ones. The rules we create during memoization, are more specific than the rules for the original general definitions, and therefore they are reordered by the system to automatically be tried first. For example, for a function

f[x_]:= f[x]= x^2

When we call it as

f /@ {1,2,3}

we end up with

?f
Global`f
f[1]=1

f[2]=4

f[3]=9

f[x_]:=f[x]=x^2

What I want to stress here is that the simplicity of memoization in Mathematica is due to the fact that we are reusing powerful core system mechanisms.

Advanced memoization

Here I will mention a few more advanced cases and /or illustrate with some examples

Algorithmic memoization example - longest common subsequence

Here I will show a less trivial example of memoization applied to a longest common subsequence problem. Here is the code:

Clear[lcs];
lcs[{}, _List] := {};
lcs[_List, {}] := {};
lcs[x_List, y_List] /; Last[x] === Last[y] := 
   Append[lcs[Most[x], Most[y]], Last[x]];
lcs[x_List, y_List] :=
   lcs[x, y] = (
     If[Length[#1] > Length[#2], #1, #2] &[
       lcs[Most[x], y], lcs[x, Most[y]]
     ]);

The resulting complexity of this algorithm is quadratic in the length of the lists, as it is also for an imperative version of the algorithm. But here the code is quite simple and self-documenting. As an illustration, here is a longest increasing subsequence for some random list:

lst=RandomInteger[{1,500},100];
lcs[lst,Sort@lst]//AbsoluteTiming

(* {0.129883,{30,89,103,160,221,227,236,250,254,297,307,312,330,354,371,374,374}} *)

The implementation can be improved by using linked lists. It also has an annoying problem that one has to manually clear the memoized values - but this issue will be addressed below.

Memoization for functions of more than one parameter, and using patterns in memoized definitions.

So far, we only looked at memoization of functions depending on a single parameter. Often the situation is more complex. One such example is here. Generally, we may want to memoize functions rather than just values. The general idiom for this is the following (I will show one extra parameter, but this is easy to generalize):

f[x_,y_]:=
   Module[{ylocal},
      f[x,ylocal_] = r.h.s (depends on x and ylocal);
      f[x,y]
   ];

What you see here is that the function is first adding a general definition, which is however computed using a fixed x, and then actually computes the value for f[x,y]. In all calls after the first one, with the same x, the newly added general definition will be used. Sometimes, as in the linked example, this may involve additional complications, but the general scheme will be the same.

There are many more examples of various flavors of this technique. For example, in his very cool solution for selecting minimal subsets, Mr.Wizard was able to use a combination of this technique and Orderless attribute in a very powerful way. The code is so short that I will reproduce it here:

minimal[sets_] :=
  Module[{f},
    f[x__] := (f[x, ___] = Sequence[]; {x});
    SetAttributes[f, Orderless];
    f @@@ Sort @ sets
  ] 

This is a more complex example of memoization, where the newly generated definitions involve patterns and this allows for a much more powerful memoization.

Caching, and selective clearing of memoized definitions

Sometimes, particularly for caching (persistent memoization), one may need to clear all or part of the memoized definitions. The basic idiom for doing so was given in this answer by Simon

ClearCache[f_] :=  DownValues[f] = 
     DeleteCases[DownValues[f], _?(FreeQ[First[#], Pattern] &)]

Sometimes, one may need more complex ways of caching, such as e.g. fixed number of cached results. The two implementation I am aware of, which automate this sort of things, are

In the linked posts there are examples of use.

Additional useful techniques / tricks

Self-blocking and automatic removal of memoized definitions

Self-blocking can be thought of as a separate design pattern. Its application to memoization is for cases when one needs the memoized values to be automatically cleared at the end of a computation. I will show two examples - Fibonacci numbers and longest common subsequence - since both of them were described above. Here is how it may look for the Fibonacci numbers:

ClearAll[fib];
fib[n_Integer]:=
  Block[{fib},         
    fib[0] = fib[1] = 1;
    fib[m_Integer?Positive]:= fib[m] = fib[m-1] + fib[m-2];
    fib[n]
  ]

You can see that the main definition redefines fib in the body of Block, and then calls fib[n] inside Block. This will guarantee that we don't generate new global definitions for fib once we exit Block.

For the Fibonacci numbers, this is less important, but for many memoizing functions this will be crucial. For example, for the lcs function, we can do the same thing:

lcs[fst_List, sec_List]:=
  Block[{lcs},
    definitions-of-lcs-given-before; 
    lcs[fst,sec]
  ]

I have used this technique on a number of occasions and find it generally useful. One related discussion is in my answer here. One particularly nice thing about it is that Block guarantees you the automatic cleanup even in cases when exception or Abort[] was thrown during the execution of its body - without any additional effort from the programmer.

Encapsulation of cached definitions

Often one may need to cache some intermediate result, but not expose it directly on the top-level, because the direct users of that result would be some other functions, not the end user. One can then use another pattern of creation of a mutable state by using Module-generated variables.

The main idiom here is

 Module[{cached},
   cached:=cached = some-computation;
   f[...]:= f-body-involving-cached;
   g[...]:= g-body-involving-cached;
 ]

Some examples involve

Memoization of compiled functions, JIT compilation

I will give here a simple example from my post on meta-programming: create a JIT-compiled version of Select, which would compile Select with a custom predicate:

ClearAll[selectJIT];
selectJIT[pred_, listType_] :=
  selectJIT[pred, Verbatim[listType]] = 
    Block[{lst},
     With[{decl = {Prepend[listType, lst]}},
      Compile @@ 
       Hold[decl, Select[lst, pred], CompilationTarget -> "C", 
          RuntimeOptions -> "Speed"]]];

This is a general method however. More examples can be found in

Symbol's auto-loading / self-uncompression

This technique is logically closely related to memoization, although perhaps is not memoization per se. I will reproduce here one function from this post, which takes a symbol and modifies its definition so that it becomes "self-uncompressing" - meaning that it stores a compressed self-modifying definition of itself:

(* Make a symbol with DownValues / OwnValues self - uncompressing *)
ClearAll[defineCompressed];
SetAttributes[defineCompressed, HoldFirst];
defineCompressed[sym_Symbol, valueType_: DownValues] :=
  With[{newVals = 
     valueType[sym] /.
       Verbatim[RuleDelayed][
         hpt : Verbatim[HoldPattern][HoldPattern[pt_]], rhs_] :>
            With[{eval = Compress@rhs}, hpt :> (pt = Uncompress@ eval)]
        },
        ClearAll[sym];
        sym := (ClearAll[sym]; valueType[sym] = newVals; sym)
];

In the mentioned post there are more explanations on how this works.

Links

Here I will list some memoization-related links which came to my mind (some of them I gave above already, but provide here once again for convenience). This list will surely be incomplete, and I invite the community to add more.

$\endgroup$
6
  • 1
    $\begingroup$ Excellent stuff Leonid, very very useful. $\endgroup$ May 13, 2013 at 22:24
  • $\begingroup$ @DavidMcHarg Thanks! One thing I dislike about how this post ended up is that it is too long. If we make many posts like this here, this pages will become slow to load and hard to navigate. I think I will eventually move some part of this answer to the page on memoization I linked to here, and here will only keep the very minimal description and the links. But ultimately this has to be decided by the community, so I would appreciate comments on this from you and others. $\endgroup$ May 13, 2013 at 22:35
  • $\begingroup$ To augment your discussion of self-blocking code, sometimes you don't need to remove existing definitions, but add to them. So, Internal`InheritedBlock is useful, instead. I recently did this in a recursive descent parser I wrote where I needed to add additional behaviors to the children of some of the tokens. $\endgroup$
    – rcollyer
    May 18, 2013 at 3:58
  • $\begingroup$ @rcollyer Yes, that's a good point. Somehow I didn't trust Internal`InheritedBlock for user-defined symbols (don't remember why - I think something went wrong a few times), so in such cases I use the technique which can be sketched as f[x_]:=With[{dv = DownValues[f]},Block[{f},DownValues[f]=dv; new-defs;code]]. But whatever the implementation, this version of technique is useful. Since I plan to have a separate post on self-blocking, which wasn't the main topic here, I did not provide all these details. $\endgroup$ May 18, 2013 at 15:06
  • 1
    $\begingroup$ @shrx Right-hand side. $\endgroup$ Jun 17, 2013 at 20:47
51
+50
$\begingroup$

Linked lists: general discussion

I would like to describe here a pattern which is IMO quite important yet seems to be generally under-appreciated and under-used in our community (and probably in Mathematica programming in general) - namely, linked lists.

Linked lists are important for several reasons. Here is a (partial) list of benefits they can bring to the table:

  • Idiomatic recursive functional programming. Linked lists allow one to use recursion in Mathematica in ways idiomatic for functional programming. Many people actually consider linked lists and recursion to be one of the corner stones of FP. Since main Mathematica data structures - lists, are actually arrays rather than linked lists, recursion using lists is not a very practical tool in Mathematica, because of performance reasons (we will discuss that in more detail below).

  • Transparent, easy to guess complexity of algorithms implemented using linked lists. This is quite important particularly in Mathematica, where a lot of internal details may be hidden behind concise programming constructs and built-in functions. In general, Mathematica is a two (or multi) scale language in terms of performance, where one performance scale is roughly given by built-in functions, while another one by the best-written top-level code (the real situation is yet more complex, but this schematic view will be good enough for us). Many algorithms which may be concisely implemented may be fast for smaller samples but then change their behavior, since their speed may be due to small time constants associated with built-in functions, rather than due to the asymptotic behavior of the resulting algorithm itself. And in any case, my point is that an reasonably efficient Mathematica solution may require an expert to quickly determine the complexity of the resulting algorithm.

    The algorithms using linked lists typically don't use the "standard" machinery such as vectorization and working with large amounts of data at once, so they are much easier to understand since they use much less of the internal machinery.

  • Predictable and usually quite reasonable speed. While solutions based on linked lists may be slower than vectorized versions, they are still among the fastest solutions one can achieve purely with the top-level code. And sometimes, they can actually lead to some of the fastest solutions for a problem (one example is a Boggle problem).

  • Memory efficiency. The solutions using linked lists are typically quite memory efficient, often much more so that vectorized solutions.

  • Generality. Because solutions based on linked lists don't usually require to have or know all the data at once, they can be used for a more general class of problems than the techniques such as vectorization or structural operations. In particular, there are several classes of problems where linked lists provide by far the most natural and clear solutions.

One final comment here is that, to an extent, declaring linked lists as a design pattern means that here we face the language limitations. Were linked lists so natural to use within the language, and we would not feel the need to describe it as a pattern. To some extent also, it is because the dominant paradigm of practical Mathematica programming is to work with as much data at once as possible and use vectorized and structural operations, an approach which is not quite compatible with linked lists and the way one works with them.

Basic linked lists

Construction

A linked list is a construct of this sort:

llist = {1,{2,{3,{4,{5,{}}}}}}

It can be conveniently constructed from a flat list, using Fold:

toLinkedList[l_List]:=Fold[{#2,#1}&,{},Reverse@l]

so that llist could have been constructed as toLinkedList[Range[5]]. The reverse operation of converting to a flat list can be also done very conveniently using Flatten:

Flatten[llist, Infinity]

(* {1, 2, 3, 4, 5}  *)

When elements are lists themselves, one needs to use generalized linked lists, which I will cover later.

The main structural elements of a linked lists are head (the first element), and tail (a list containing the rest of the elements). Note that tail is also a linked list. Extracting one or the other is very easy - just use First for head and Last or tail:

First@llist

(* 1 *)

Last@llist

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

The defining algorithmic property and an illustration - accumulation of results.

In the algorithmic aspect, the main, defining property of a linked list, which makes it very different from usual Mathematica list, is that it is cheap to copy. The reason is that any linked list is a list of two elements only - head and tail. And because expressions in Mathematica are immutable, the copying of tail does not involve a deep copy, but rather (internally) a pointer to the list is copied. This is in contrast to normal Mathematica lists, whose copying will in most cases have a linear complexity in the size of the list.

Everything else follows from this property. For example, collecting of intermediate results is easy with linked lists. The main idiom here is

list = {new-elem, list}

One practical example of using this idiom can be found here. To illustrate this here, consider a toy implementation of Select based in linked lists:

selectLL[l_List, pred_] :=
   Module[{result  = {}},
      Do[If[pred[el], result = {el, result}], {el, l}];
      Reverse@Flatten[result]];

Here is how it is used (note that it will only work correctly for list elements which are not themselves lists):

selectLL[Range[10], OddQ]

(* {1, 3, 5, 7, 9} *)

Now, the analog we may write using normal lists will be

selectNaive[l_List, pred_] :=
  Module[{result  = {}},
    Do[If[pred[el], AppendTo[result, el]], {el, l}];
    result];

But in this case, entire list is copied every time when a new element is added, which results in much worse performance even for relatively small lists:

selectLL[Range[20000],OddQ];//AbsoluteTiming
selectNaive[Range[20000],OddQ];//AbsoluteTiming

(* 
 {0.037109,Null}
 {0.561523,Null}
*)

Note by the way that the built-in Select is not dramatically faster than our implementation based on linked lists:

Select[Range[20000], OddQ]; // AbsoluteTiming

(* {0.014648, Null} *)

One can argue that Sow and Reap can be used for accumulation of results, and in fact may be more efficient in that than linked lists. This is true, but the Sow - Reap model is more limited because partial lists of results accumulated "so far" are not available to the user.

Generalized linked lists

The general construction for linked lists is similar to the one described above, but the List head is replaced by some (arbitrary) head (say,ll). For certain reasons, it is better in many cases if ll is HoldAllComplete (basically, to prevent search for UpValues, which may be expensive and slow things down). So, we end up with

ClearAll[ll];
SetAttributes[ll,HoldAllComplete];
toLinkedListGen[l_List]:= Fold[ll[#2,#1]&,ll[],Reverse@l]

for example

llgen = toLinkedListGen[Range[5]]

(* ll[1, ll[2, ll[3, ll[4, ll[5, ll[]]]]]] *)

and the code to convert this back to the flat list may look like

List @@ Flatten[llgen, Infinity, ll]

(* {1, 2, 3, 4, 5} *)

This construction is valid for lists with arbitrary elements (except those having themselves the head ll). If you need to use linked lists inside other linked lists, you can define several such heads.

One thing to keep in mind here is that you have to make sure that the element you are adding has been fully evaluated (if that's essential) before writing

llist = ll[elem, llist]

if ll is HoldAllComplete.

Typical use cases and examples

Here I will try to show some examples and give some links for the typical use cases of linked lists, to achieve various goals.

Implementing algorithms requiring complex list iteration, possibly with look-ahead analysis.

A typical such problem is the merge sort algorithm. I will show here the implementation I posted here:

Clear[toLinkedList]; 
toLinkedList[x_List] := Fold[{#2, #1} &, {}, Reverse[x]];
Module[{h,lrev},        
    mergeLinked[x_h,y_h]:=
        Last[
            {x,y,h[]}//. 
                {                        
                    {fst:h[hA_,tA_h],sec:h[hB_,tB_h],e_h}:>
                        If[hA>hB,{tA,sec,h[hA,e]},{fst,tB,h[hB,e]}],
                    {fst:h[hA_,tA_h],h[],e_h}:>{tA,h[],h[hA,e]},
                    {h[],sec:h[hB_,tB_h],e_h}:>{h[],tB,h[hB,e]}
                }
        ];

    lrev[set_]:=
        Last[h[set,h[]]//. 
            h[h[hd_,tl_h],acc_h]:>h[tl,h[hd,acc]]];

    sort[lst_List]:=
        Flatten[                
            (h[#1,h[]]&)/@lst//. 
                x_List:>
                    Flatten[
                        {toLinkedList[x],{}}//. 
                            {{hd1_,{hd2_,tail_List}},accum_List}:>
                                {tail,{accum,lrev[mergeLinked[hd1,hd2]]}}
                    ],
            \[Infinity],
            h
        ]
]

Here, two different linked lists are used. One is for a list of already merged lists, which at every iteration gets half its previous lengths, since already ordered sub-lists are merged pairwise. For this one, a normal List head is used, and this "higher-level" merging procedure is implemented in sort. Another one is used for individual sub-lists which are merged - they are first converted to linked lists with head h. The lower-level merging procedure is implemented in mergeLinked, while lrev reverses such linked list.

This implementation is a bit unusual in that it uses exclusively linked lists, local rules and ReplaceRepeated, instead of recursion. It was originally posted on MathGroup, where it is explained in detail, and the corrected version I use in this post was posted here. The main idea was to illustrate that ReplaceRepeated can be just as efficient as other means, when used with appropriately constructed patterns. My implementations were certainly motivated by a similar discussion in the book of David Wagner.

In any case, the naive implementation using normal lists, which can look e.g. like

Clear[merge];
merge[x_List,y_List]:=
    Block[{merge},
        Flatten[
            merge[x,y]//. 
                {                        
                    merge[{a_,b___},{c_,d___}]:>
                        If[a<c,{a,merge[{b},{c,d}]},{c,merge[{a,b},{d}]}],
                    merge[{},{a__}]:>{a},
                    merge[{a__},{}]:>{a}
                }
        ]
    ]

Will have an unacceptable complexity due to the list copying at every stage.

Notes on recursion

Note that while in the previous example recursion was not used explicitly, this was rather artificial, and normally the solutions involving linked lists are recursive. To finish with the merge sort example and illustrate this point, here is the recursive implementation of merge, taken from here

ClearAll[merge];
merge[a_ll, ll[], s_, _] := List @@ Flatten[ll[s, a], Infinity, ll];
merge[ll[], b_ll, s_, _] := List @@ Flatten[ll[s, b], Infinity, ll];
merge[ll[a1_, atail_], b : ll[b1_, _], s_, f_: LessEqual] /;f[a1, b1] :=
    merge[atail, b, ll[s, a1], f];
merge[a : ll[a1_, _], ll[b1_, brest_], s_, f_: LessEqual] :=
    merge[a, brest, ll[s, b1], f];
merge[a_List, b_List, f_: LessEqual] :=
    merge[toLinkedListGen@a, toLinkedListGen@b, ll[], f];

(where I replaced toLinkedList with toLinkedListGen, the latter defined previously in this post).

Four additional important things to note here: first, note how we can effectively combine pattern-matching on linked lists and function overloading (more than one definition / global rule), to end up with code which is quite expressive and self-documenting.

Second, note that the last definition is for normal lists, and automatically translates them to linked lists first, while the first two definitions automatically translate the result into the flat list again. This is a common technique, it allows the function to be self-contained, in that linked lists are only used internally, while both input and output is in terms of normal lists.

Third, you can see that the list of accumulated results (which is also a linked list), is being passed to the function as its third argument. This is also a common /standard technique - make the function pass all additional information that it changes in the process, as extra parameter(s). Not only does this allow one to avoid introducing mutable state into the function (to hold this changing information), but also, and perhaps more importantly, it allows the function to be tail-recursive in the Mathematica sense. This is an absolute requirement for recursive solutions to be practical. The main reason why the function won't be tail-recursive otherwise is explained in the link I just gave, but basically you will have something like

f[x_]:= (
   Do-something-with-mutable-state;
   f[x-1]
)

And the presence of CompoundExpression around is what makes it non-tail-recursive.

Last thing to note is that, for large enough lists, one may need to increase the setting of $IterationLimit, or set it to infinity altogether. This can always be done on the level of function's definitions. For example, for merge, we could change the last definition to

merge[a_List, b_List, f_: LessEqual] :=
  Block[{$IterationLimit = Infinity},
    merge[toLinkedListGen@a, toLinkedListGen@b, ll[], f]
  ];

Note also that if the system complains about $RecursionLimit being exceeded, this means that your function is not properly tail-recursive (as far as Mathematica is concerned). All properly tail-recursive Mathematica functions can only exceed $IterationLimit, but not $RecursionLimit.

General list processing based on some look-ahead analysis

A typical class of use cases for linked lists consists of problems of splitting lists into various groups, according to their structure. I will show a simple example, takesn from here, where the built-in Split would do, but in general Split will not suffice, since the splitting condition may involve more than just two adjacent elements.

This function will count the lengths of sequence of adjacent elements, using linked lists:

Clear[cValsRec];
cValsRec[l_List] := cValsRec[toLinkedList@l, {{}, 0}];
cValsRec[{}, counts_] := Most@Flatten[counts];
cValsRec[{val_, tail : {val_, _List}}, {prev_List, count_}] :=
    cValsRec[tail, {prev, count + 1}];
cValsRec[{val_, tail_List}, {prev_List, count_}] :=
    cValsRec[tail, {{prev, count + 1}, 0}];

and it can be used e.g. as

test = Sort@RandomInteger[5, 15]

(* {0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5} *)

cValsRec[test]

(* {2, 3, 3, 4, 3}  *)

As before, you can see the same pattern: the first definition is a "starting" one and converts initial list to a linked list. The last definition describes the final stage and conversion the result (which is also accumulated in a linked list) back to normal list. And the results are accumulated in a separate linked list which is passed to the next function call as a separate argument.

Sometimes, instead of splitting the list, one has to do something else, like e.g. deleting or adding some elements, but also based on near-by elements. All these cases are problems for which functional iterators such as Fold are generally ill-suited. I also add such cases to the set of cases I cover in this section, since they are conceptually similar.

Some further links for similar problems solved using linked lists:

I particularly recommend to look into the code of the last two links, since those use a number of interesting techniques in combination.

Some general remarks for this section: first, I'd like to stress one clear advantage of linked lists for this type of problems: they are completely straight-forward to use. You don't have to figure out some clever vectorizing scheme or sequence of list-manipulation steps. You don't have to know all subtleties of various built-in functions involved, and their combinations. This means that linked lists provide a straight-forward implementation path for users who are not Mathematica experts.

Second, the code using linked lists and recursion tends to be readable and self-documenting. This is quite important, since cryptic code is hard to understand and maintain, and it is easier to create bugs when the code is cryptic.

Third, code based on linked lists can be completely generic, since it does not rely on vectorization, packed arrays etc. At the same time, it remains completely idiomatic, since linked lists are immutable, and no mutable state is used. Admittedly, for numerical problems where the structure of the data may allow (much) more efficient vectorized solutions, linked lists may not be the fastest. But in those cases, one can dispatch to such special cases automatically, as it has been done e.g. for the merge function here.

Using linked lists for tree traversals and tree-building

This is a more advanced use of linked lists, and here they can fully reveal their potential. The thing is, problems involving trees are usually naturally recursive, and while for flat list creation / processing it is very often possible to get away by using standard Mathematica arsenal for list manipulation, here itis often not enough.

While I plan to extend this part of the post soon, right now I will illustrate this with one relatively simple example - a solution for a variant of the cutting stock problem in Mathematica, and then provide links to a few more involved cases.

So, the problem there was to list the ways in which one could sum elements from a list (with repetitions), so the sum is less-than-or-equal to some value. So, here are the values of the list and the maximum value, from the original question:

i = {7.25, 7.75, 15, 19, 22};
m = 22;  

The solution using linked lists is fairly simple:

ClearAll[getCombs];
getCombs[l_List,lim_Integer]:=
    Cases[            
        getCombs[{},toLinkedList[Reverse[Sort[l]]],lim],
        getCombs[acc:Except[{}],__]:>Flatten[acc],
        \[Infinity]
    ];
getCombs[accum_,l:{head_,tail_List},lim_?Positive]:=
    {            
        If[lim>=head,getCombs[{accum,head},l,lim-head],Sequence@@{}],
        getCombs[accum,tail,lim]
    };

The first definition is both the "starting" and "ending" one, and it realizes the "wrap-around definition" pattern - that is, preprocess data (in this case, convert to linked list the original list), call the main function, and then post-process the data (Cases here, I will explain what it does in a moment).

The second definition is the one doing the real work. Borrowing the explanation from the original answer, the way it works is to create a tree by going backwards, subtracting a given element from the current maximal for the sum, to get the maximal for the sums of remaining elements. The case of repetitions is accounted for by the line If[lim>=head,...]. The resulting numbers are also accumulated in linked lists (accum parameter), and finally they are found with Cases and converted back to normal lists using Flatten. In this way, we avoid unnecessary array copying.

As you can see, the most important difference here is that a tree is created as a result, rather than a list. That tree is then post-processed, and the results collected from it, using Cases. This construction allows us to take the best from both worlds: linked lists + recursion, and pattern-based expression destructuring. And since the patterns are "syntactic" See my comment to a post here, it is really fast. Because the algorithm is essentially recursive, and the result is a tree rather than a more regular list, the standard list-manipulation approaches are at a loss here.

Here are some more links to applications of linked lists involving trees:

In some of the examples linked, the use of linked lists or associated data structures (and trie can be viewed as a generalization of a linked list) led to some quite dramatic performance boosts. You can read explanations in the linked posts for more details.

Remarks

I tried to outline the usefulness of linked lists and illustrate the main ideas behind their use, and techniques associated with them. One of the main advantages they offer is the possibility of direct and straight-forward implementation of recursive algorithms in Mathematica, with quite reasonable performance.

This post is a work in progress. I plan to add more links, and make a special separate link section where all links will be collected. I also may add a couple of more sections, such as a section of using linked lists to trade speed for memory.

All remarks, suggestions and additions are welcome and encouraged - this is a community post as are all answers here.

$\endgroup$
6
  • 1
    $\begingroup$ Not to scoff (I did give this an upvote) but I think these are not quite true linked lists. If I am correct, then they suffer for certain types of computation in consequence. Specifically, it seems one cannot add a new cell in constant time in the middle of a list (assuming one already has a "pointer" to the parent, say). Best I can tell, one would have to use Part[...] where the ... sequence has as many elements as is needed to reach the cells between which we are adding a new one. $\endgroup$ Jun 4, 2013 at 14:33
  • $\begingroup$ @DanielLichtblau You are right. There are no pointers for these objects in the sense of e.g. C. So, they are more restricted than the "normal" linked lists, but other than being immutable, they still are linked lists. So, this restriction seems to have more to do with the general immutability of Mathematica than with these constructs not really being linked lists. I think. $\endgroup$ Jun 4, 2013 at 15:32
  • $\begingroup$ It's partly philosophy, but...one really ought to be able to add and remove links in a "linked list". (I realize there are many applications that do not require this particular facet of such entities.) Uh oh, this is starting to sound like the ongoing "What should we be allowed to call a List?" tempest on MathGroup. $\endgroup$ Jun 4, 2013 at 15:48
  • $\begingroup$ @DanielLichtblau Well, I would say that this construct is as close to a linked list as one can get in Mathematica if we exclude some top-level emulations of pointer semantics which are almost certainly doomed to be slower or much slower. But I agree with you that this construct (as described in my post above) falls short of the expected linked list's behavior in a few practically important aspects. Also agree with Mathgroup remark :). $\endgroup$ Jun 4, 2013 at 15:56
  • 1
    $\begingroup$ @Daniel and Leonid What would be a good list-like data structure if I want to remove elements from arbitrary positions? Say I have a list, a function selects an element from the list to process it further and removes it from the list. Then rinse and repeat until the list is empty. (I hope that wasn't too vague.) $\endgroup$
    – sebhofer
    Sep 20, 2013 at 8:42
25
$\begingroup$

"Design patterns" are a standard approach to a class of problems. They may become idioms if the community adopts them. They may be simple, but non-trivial to figure out and identify (for oneself) as a pattern worthy of re-use.

Here's one particular pattern I personally like to use:

Example: the "TransformBy" approach

There are functions which transform lists based on the properties of list elements, e.g. Sort, Split, Gather. Often it's advantageous to do the transformation based on a transform of the list elements, so we have SortBy, SplitBy, GatherBy. It's particularly useful when we deal with equivalence classes which are easier to find based on the comparison f[x] == f[y] than a completely custom comparison function, or for ordering-related problems (put < in place of ==).

This approach becomes a "design pattern" once I internalise it enough that I start approaching problems in terms of it.

There's no built-in MaxBy or MinBy (probably because the result of these wouldn't necessarily be unique), but they're often very useful, so I define them

MaxBy[list_, fun_] := list[[First@Ordering[fun /@ list, -1]]]
MinBy[list_, fun_] := list[[First@Ordering[fun /@ list, 1]]]

Then we can use them in so many ways. I commonly do things like MaxBy[dataTable, Last], MinBy[dataList, Variance], MaxBy[lists, Length], MinBy[points, Norm], etc.

What makes this a "design pattern" is that I internalised the approach enough that it comes as naturally to me as using built-ins such as Outer. Or, to put it less subjectively, that the pattern turned out to be general enough that it was possible and useful for me to internalise it.

Other candidates for a ...By pattern (with caveats): Union, Intersection, Tally, DeleteDuplicates etc. All these deal with equivalence classes and some of them speed up operations using sorting. The ...By version would typically be faster than giving a custom comparison function (which precludes using sorting or binary search for speedup).

Some other general patterns:

  • Memoization in the form of f[k_] := f[k] = ...

  • The Villegas-Gayley trick

  • Vectorization: instead of processing a list (or lists) element wise, use operations that work with whole lists and keep them packed. I'm not sure if I'd call this a design pattern or a paradigm, but in Mathematica, the pattern of vectorizing things is rather standard: Use arithmetic operations on lists; use UnitStep for comparisons (>, >=, <, <=); use 1 and 0 for True/False and add/multiply them for and and or operations; use Pick for conditional extraction. Here's a (very simple) example. Doing this usually leads to a significant speedup, but it also creates code that is not easy and quick to modify. Some other languages, such as MATLAB, provide much better support for this style of programming, but it's possible in Mathematica too (with the same performance benefits).

  • ...Block for resource management, analogous to JavaBlock or NETBlock. This is a bit specialized an advanced, and it's useful for package authors. It's a very useful pattern for simplifying resource management when one needs to deal with objects that have to be created and destroyed explicitly. (I extended TetGenLink to do this as well by wrapping TetGenCreate using Villegas-Gayley to keep track of created objects.) See also here.

  • One could argue that the standard package structure is a design pattern as well: BeginPackage[...], mention public symbols (coupling this conveniently with setting usage messages), Begin["`Private`"], definitions, End[]; EndPackage[]. Theoretically it could be done differently. There's nothing forcing me to write packages using this pattern. But this is a time tested pattern for a package structure that makes things easy to manage and robust.

  • Creating closures (i.e. functions with a persistent state) using Module, as shown here.

  • (will update as I come across them)

$\endgroup$
7
  • $\begingroup$ Some ideas: closures, using Reap and Sow, Options, using Block to temporarily have global variables, With to inject code inside held expressions (Villegas-Gayley trick), use of Compile for speed critical code, use of functional programming, HoldAll attributes, Dispatch, any form of struct (the most simple one is struct["name"]=property) $\endgroup$
    – faysou
    May 13, 2013 at 6:00
  • $\begingroup$ Many examples can be found here mathematica.stackexchange.com/a/8617/66 and could be used in this post. $\endgroup$
    – faysou
    May 13, 2013 at 6:01
  • $\begingroup$ Other ideas: nested lists, lazy lists, using UpValues on structures of the form headType[data1,....] $\endgroup$
    – faysou
    May 13, 2013 at 12:24
  • $\begingroup$ @Faysal Some good examples (e.g. the use of closures), but some of these I wouldn't consider a design pattern. Using Sow/Reap or using Compile is nothing more than using a standard language feature. Compare this with implementing closures using Module which involves bringing together several language features (Module, Function, Temporary` variables) in a specific and non-ovbious way to get something new that is not directly given to you by the language. $\endgroup$
    – Szabolcs
    May 13, 2013 at 13:18
  • $\begingroup$ Yes I agree that some examples are standard language features, but they are still methods often used to solve problems, and not common to other languages, idioms like you said. $\endgroup$
    – faysou
    May 13, 2013 at 15:22
22
$\begingroup$

This is just a beginning of some notions of what the patterns of Mathematica programming might be. Others should feel free to add or correct anything here. It might well be that I've simply misunderstood what patterns in software development really are.

Pattern-based overloading to create polymorphous functions

It's pretty common to use Mathematica's pattern-matching functionality to create functions that can do vastly different things depending on the structure of the input. A trivial example:

f[x_Integer]:= x^2
f[x_String]:= StringReverse[x]

Head substitution (Apply @@, @@@) to transform input

Before Total came on the scene in version 5, we all got used to doing things like

Plus @@@ Table[RandomReal[], {50}, {50}]

What this is really doing is changing something with head List to something with head Plus, which just happens in this case to create a summation. But it can be used to transform inputs in non-trival ways. In this example, the symbolic function f is removed entirely.

m = f /@ {{a, b, c}, {b, c, d}};
Accumulate @@@ m
(* {{a, a + b, a + b + c}, {b, b + c, b + c + d}} *)

Replacement rules to modify existing results

We all know that replacement rules (/.) can be used for substituting numerical values for symbolic parameters in algebraic expressions. But Mathematica's pattern-matching functionality, together with replacement rules, can be used more generally to transform intermediate results in quite sophisticated ways.

data = Table[RandomInteger[{10, 30}], {50}, {50}];
data2 = data //. {p___, a_, b_?PrimeQ, q___} -> {p, a, a, q};
$\endgroup$
1
  • 6
    $\begingroup$ Not to criticize, since I think this is an important question, but in my opinion these are really examples of idioms, rather than design patterns as commonly understood. Design patterns used in Mathematica might be e.g. the use of Condition to preprocess arguments into the form required inside a function as well as checking their validity, or the declaration of symbols in packages by giving their usage messages before beginning the Private` context. $\endgroup$ May 10, 2013 at 16:33
10
$\begingroup$

Although some of the techniques below use the core language, they could nonetheless be new to a beginner or intermediate Mathematica user.

Many illustrations can be found in my answer here: Where can I find examples of good Mathematica programming practice? that collects interesting questions and answers, I'm just trying to summarize here some techniques I often see used.

  • Closures: a Module can return a pure function that depends on the Module's created symbols. The pure function will work even outside of the Module.

    newIdFunction = Module[{id=0},id++&];
    newIdFunction[]
    newIdFunction[]
    
  • Use Block for temporary global variables. If you use packages be careful to use the correct variable name in the right context.

    f[]:=a+1
    Block[{a=0},f[]]
    
  • Block a function from executing

    myFunction[x_] := x
    myFunction[x] // Hold@# &
    Block[{myFunction},myFunction[x] // Hold@# &]
    
    • An interesting application of blocking a function from executing is to create a VerboseBlock which can be used to toggle on/off Print statements imbedded in your code.

      ClearAll[VerboseBlock];
      SetAttributes[VeboseBlock, HoldRest];
      VerboseBlock[isVerbose : True | False, block_] := 
       With[{verbose = If[isVerbose, {}, {Print}]}, Block[verbose, block]]
      

      This is usually accompanied by an option passed to your code that toggles the printing on/off:

      Options[mfcn] = {Verbose -> False};
      mfcn[x_, OptionsPattern[]] := (Print[x]; x^2)
      mfcn[4]
      (* 16 *)
      
      mfcn[4, Verbose -> True]
      (* 4
         16
      *)
      

      This could be adapted to a more robust scheme with tag based logging, but I leave that as an exercise to the reader.

  • Replace inside a Held expression or a function with a HoldAll attribute using With.

    With[{x=2},Hold@x]
    

    If you also need to evaluate what you need to replace use the Trott-Strzebonski in-place evaluation trick

  • Any form of struct (the most simple one is struct["name"]=property)

  • Using UpValues on structures of the form headType[data1,....], this allows to do some basic object oriented programming in Mathematica. More on this here.

    x = myHead[2,3];  
    y = myHead2[2,3];  
    f[myHead[x_,y_]] ^:= x+y;  
    f[myHead2[x_,y_]] ^:= x y;  
    f[x]  
    f[y]  
    
  • Temporarily modify the attributes or the definition of a built-in function

    Internal`InheritedBlock[{Map},Unprotect@Map;Map[x___]:={x};Map[1,2,3]]
    

    Comment from rcollyer:
    Internal`InheritedBlock retains the current definition so it can be augmented. Better examples are in BlockStream or withTaggedMsg both of which make use of the original behavior after running the modifications

  • Use OptionsPattern[], for optional arguments

    Options[f]={"PrintInput"->False};
    f[x_,OptionsPattern[]]:=(If[OptionValue@"PrintInput",Print["x=",x];];x)
    f[2,PrintInput->True]  
    
  • Using Reap and Sow

  • Use of Compile for speed critical code
  • HoldAll attributes and defered evaluation
  • Dispatch
$\endgroup$
1
  • $\begingroup$ Internal`InheritedBlock retains the current definition so it can be augmented. As written, you do a complete override. Better examples are in BlockStream or withTaggedMsg both of which make use of the original behavior after running the modifications. $\endgroup$
    – rcollyer
    Jun 8, 2013 at 2:58
8
$\begingroup$

One of my favorite Mathematica paradigms is what we used to call "self modifying code." Back in assembly language days (before the dinosaurs) this was considered a VERY BAD THING and was discouraged -- many modern languages forbid it. But it's a normal and common part of the Mathematica syntax. For example, say you need to implement a filter on an image, but you might want to do many kinds of filters depending on what kind of image it is. In this case it is very common to write the line

 filt[image]

and then in the preceeding code to define filt=GaussianFilter or filt=BandPassFilter or filt=RangeFilter. This "self modifying code" strategy is endemic throughout built-in functions as well. Consider something as simple as ListConvolve. Normally (in most languages I have ever seen) the convolution function combines two sets of data using multiplication and summation. Mathematica's ListConvolve takes two optional arguments that let you specify what "multiplication" means and what function "addition" is going to be. This is a very powerful construct.

$\endgroup$
8
  • 3
    $\begingroup$ It's powerful, yes, but I don't see how it is "self-modifying." Self-modifying code in MMA can be written but rarely is: it would consist of editing string expressions that would then be interpreted as MMA code and executed. What you discuss here appears merely to be the ability to pass almost any object as parameters, including objects that are thought of as executable. That has not historically been considered a "very bad thing" (although the dangers of such a thorough lack of typing are well known). $\endgroup$
    – whuber
    May 10, 2013 at 17:59
  • 1
    $\begingroup$ It's self modifying because by reading the line of code filt[image] you have no idea what operation will take place (until the program actually runs). This is similar to what one used to do in assembly language: build up a string or a number and then execute that number as a command. I was trying to give a very simple example here where it seems perfectly natural. As you say, it is a mixing of typing: in this case, a variable that is set in the code (filt=GaussianFilter) becomes a function. $\endgroup$
    – bill s
    May 10, 2013 at 18:03
  • $\begingroup$ Not knowing what operation will take place is similar to not knowing (at compile time) what values the parameters will have, either. "Self-modifying" does not mean "unknown at compile time." Your example would be implemented in Assembly or C or their ilk using a dispatch table. Although the exact procedure will not be known until runtime, there is control over which procedures can be executed and in particular any procedure that might get executed has already been created before the execution started. During runtime, no code is actually modified. $\endgroup$
    – whuber
    May 10, 2013 at 18:11
  • 3
    $\begingroup$ @whuber I agree with your points except that I think it should be noted that there is no distinction between code and data in Mathematica, so conversion to a string and back is not necessary. Self-modifying code could be implemented by transforming a held expression, or even without any explicit transformation--memoization is, I think, a good example, as the need to distinguish between computing a value and retrieving it from a LUT is avoided by defining new downvalues. Here, admittedly, the line between the code and the process of evaluating it is somewhat blurred. $\endgroup$ May 10, 2013 at 19:47
  • 5
    $\begingroup$ Memoization by itself is a very important design pattern in Matheamtica: used in the right way I could get amazing speed-ups for real-world consulting projects over the years. Lately I used transformation of held expressions to transfrom f \[RightTeeArrow] (2 x) (which displays nicely in a notebook) to f[x_]:=f[x]=(2 x), where the independent variables need to be explicitly specified. Now, if only I could control the precedences of the mysterioius FrontEnd's parsing ... (but this is not so serious, just put () around the right hand side). $\endgroup$ May 10, 2013 at 20:33
4
$\begingroup$

Necessary? No they are not necessary, but they are extremely useful, just as in any other programming system. I would say that the abstraction level of a programming language has no effect on the usefulness of design patterns.

There must be hundreds, perhaps thousands, of Mathematica design patterns.

I recently answered a question concerning table building with what is essentially a design pattern although I did not mention the term "design pattern" in my answer (it would have been better if I had but it slipped my mind). I state it now as an example of a Mathematica design pattern.

The question was (paraphrase): how can a 2D table be built in which the row index i does not advance through an arithmetic progression? The design pattern that answers this question is

indexList = selectorFunction[Range[maxI]];
table = Table[elementFunction[i, j], {i, indexList}, {j, maxJ}]
$\endgroup$
3
$\begingroup$

OOP Design Patterns

Many of the Object-Oriented Programming (OOP) Design Patterns can be (relatively) easily implemented in Mathematica without using preliminary code (packages). This answer in the discussion "General strategies to write big code in Mathematica?" provides more details.

We can say that the Design Patterns implementations discussed in the links above emulate OOP in Mathematica at Design Patterns level instead of at objects level. This allows the emulation to be much more straightforward and it is especially important in relation to (as said in the question):

But some people argue that identified "patterns" are really just a sign that the language in question isn't abstracted enough.

(I.e. if we emulate OOP with objects then we will emulate the shortcomings of OOP languages.)

This document:

"Implementation of Object-Oriented Programming Design Patterns in Mathematica" ,

has examples of OOP Design patterns that are already incorporated in the "standard" Mathematica: see Strategy, Decorator, Observer / MVC.

Monadic programming

This functional parsers implementation in Mathematica, FunctionalParsers.m, introduced and demonstrated in these blog posts, is an application of the so called "monadic programming" style. Recent developments in Mathematica make the use of that style much easier.

The monad for the functional parsers in FunctionalParsers.m is defined for the functions

P: {{_String...}, _ParseTree } -> { {{_String...}, _ParseTree } ... }

This definition allows the chaining of parsers that are of basic, combinator, or transforming types. (This kind of chaining is discussed here.)

For more detailed definitions and explanations see "Functional parsers for an integration requests language grammar".

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.