What is in your Mathematica tool bag? [closed]

Question

We all know that Mathematica is great, but it also often lacks critical functionality. What kind of external packages / tools / resources do you use with Mathematica?

I'll edit (and invite anyone else to do so too) this main post to include resources which are focused on general applicability in scientific research and which as many people as possible will find useful. Feel free to contribute anything, even small code snippets (as I did below for a timing routine).

Also, undocumented and useful features in Mathematica 7 and beyond you found yourself, or dug up from some paper/site are most welcome.

Please include a short description or comment on why something is great or what utility it provides. If you link to books on Amazon with affiliate links please mention it, e.g., by putting your name after the link.

---

Packages:

1. LevelScheme is a package that greatly expands Mathematica's capability to produce good looking plots. I use it if not for anything else then for the much, much improved control over frame/axes ticks. Its newest version is called SciDraw, and it will be released sometime this year.
2. David Park's Presentation Package (US$50 - no charge for updates)
3. Jeremy Michelson's grassmannOps package provides resources for doing algebra and calculus with Grassmann variables and operators that have non trivial commutation relations.
4. John Brown's GrassmannAlgebra package and book for working with Grassmann and Clifford algebras.
5. RISC (Research Institute for Symbolic Computation) has a variety of packages for Mathematica (and other languages) available for download. In particular, there is Theorema for automated theorem proving, and the multitude of packages for symbolic summation, difference equations, etc. at the Algorithmic Combinatorics group's software page.

Tools:

1. MASH is Daniel Reeves's excellent Perl script essentially providing scripting support for Mathematica v7. (Now built in as of Mathematica 8 with the -script option.)
2. An alternate Mathematica shell with a GNU readline input (using python, *nix only)
3. ColourMaths package allows you to visually select parts of an expression and manipulate them. http://www.dbaileyconsultancy.co.uk/colour_maths/colour_maths.html

Resources:

1. Wolfram's own repository MathSource has a lot of useful if narrow notebooks for various applications. Also check out the other sections such as

   • Current Documentation,
   • Courseware for lectures,
   • and Demos for, well, demos.

2. The Mathematica Wikibook.

Books:

1. Mathematica programming: an advanced introduction by Leonid Shifrin (web, pdf) is a must read if you want to do anything more than For loops in Mathematica. We have the pleasure of having Leonid himself answering questions here.
2. Quantum Methods with Mathematica by James F. Feagin (amazon)
3. The Mathematica Book by Stephen Wolfram (amazon) (web)
4. Schaum's Outline (amazon)
5. Mathematica in Action by Stan Wagon (amazon) - 600 pages of neat examples and goes up to Mathematica version 7. Visualization techniques are especially good, you can see some of them on the author's Demonstrations Page.
6. Mathematica Programming Fundamentals by Richard Gaylord (pdf) - A good concise introduction to most of what you need to know about Mathematica programming.
7. Mathematica Cookbook by Sal Mangano published by O'Reilly 2010 832 pages. - Written in the well known O'Reilly Cookbook style: Problem - Solution. For intermediates.
8. Differential Equations with Mathematica, 3rd Ed. Elsevier 2004 Amsterdam by Martha L. Abell, James P. Braselton - 893 pages For beginners, learn solving DEs and Mathematica at the same time.

Undocumented (or scarcely documented) features:

7. How to customize Mathematica keyboard shortcuts. See this question.
8. How to inspect patterns and functions used by Mathematica's own functions. See this answer
9. How to achieve consistent size for GraphPlots in Mathematica? See this question.
10. How to produce documents and presentations with Mathematica. See this question.

Answer 1

Todd Gayley (Wolfram Research) just send me a nice hack which allows to "wrap" built-in functions with arbitrary code. I feel that I have to share this useful instrument. The following is Todd's answer on my question.

A bit of interesting (?) history: That style of hack for "wrapping" a built-in function was invented around 1994 by Robby Villegas and I, ironically for the function Message, in a package called ErrorHelp that I wrote for the Mathematica Journal back then. It has been used many times, by many people, since then. It's a bit of an insider's trick, but I think it's fair to say that it has become the canonical way of injecting your own code into the definition of a built-in function. It gets the job done nicely. You can, of course, put the $inMsg variable into any private context you wish.

Unprotect[Message];

Message[args___] := Block[{$inMsg = True, result},
   "some code here";
   result = Message[args];
   "some code here";
   result] /; ! TrueQ[$inMsg]

Protect[Message];

Answer 2

Caching expressions

I find these functions very helpful to cache any expression. The interesting thing here for these two functions is that the held expression itself is used as a key of the hashtable/symbol Cache or CacheIndex, compared to the well-known memoization in mathematica where you can only cache result if the function is defined like f[x_] := f[x] = ... So you can cache any part of a code, this is useful if a function is to be called several times but just some parts of the code must not be recomputed.

To cache an expression independently of its arguments.

SetAttributes[Cache, HoldFirst];
c:Cache[expr_] := c = expr;

Ex: Cache[Pause[5]; 6]
Cache[Pause[5]; 6]

The second time the expression returns 6 without waiting.

To cache an expression using an alias expression that can depend on an argument of the cached expression.

SetAttributes[CacheIndex, HoldRest];
c:CacheIndex[index_,expr_] := c = expr;

Ex: CacheIndex[{"f",2},x=2;y=4;x+y]

If expr takes some time to compute, it is much faster to evaluate {"f",2} for example to retrieve the cached result.

For a variation of these functions in order to have a localized cache (ie. the cache memory is automatically released outside the Block construct) see this post Avoid repeated calls to Interpolation

Deleting cached values

To delete cached values when you don't know the number of definitions of a function. I consider that definitions have a Blank somewhere in their arguments.

DeleteCachedValues[f_] := 
       DownValues[f] = Select[DownValues[f], !FreeQ[Hold@#,Pattern]&];

To delete cached values when you know the number of definitions of a function (goes slightly faster).

DeleteCachedValues[f_,nrules_] := 
       DownValues[f] = Extract[DownValues[f], List /@ Range[-nrules, -1]];

This uses the fact that definitions of a function are at the end of their DownValues list, cached values are before.

Using symbols to store data and object-like functions

Also here are interesting functions to use symbols like objects.

It is already well known that you can store data in symbols and quickly access them using DownValues

mysymbol["property"]=2;

You can access the list of keys (or properties) of a symbol using these functions based on what dreeves submitted in a post on this site:

SetAttributes[RemoveHead, {HoldAll}];
RemoveHead[h_[args___]] := {args};
NKeys[symbol_] := RemoveHead @@@ DownValues[symbol(*,Sort->False*)][[All,1]];
Keys[symbol_] := NKeys[symbol] /. {x_} :> x;

I use this function a lot to display all infos contained in the DownValues of a symbol:

PrintSymbol[symbol_] :=
  Module[{symbolKeys},
    symbolKeys = Keys[symbol];
    TableForm@Transpose[{symbolKeys, symbol /@ symbolKeys}]
  ];

Finally here is a simple way to create a symbol that behaves like an object in object oriented programming (it just reproduces the most basic behaviour of OOP but I find the syntax elegant) :

Options[NewObject]={y->2};
NewObject[OptionsPattern[]]:=
  Module[{newObject},
    newObject["y"]=OptionValue[y];

    function[newObject,x_] ^:= newObject["y"]+x;
    newObject /: newObject.function2[x_] := 2 newObject["y"]+x;

    newObject
  ];

Properties are stored as DownValues and methods as delayed Upvalues in the symbol created by Module that is returned. I found the syntax for function2 that is the usual OO-syntax for functions in Tree data structure in Mathematica.

For a list of existing types of values each symbol has, see http://reference.wolfram.com/mathematica/tutorial/PatternsAndTransformationRules.html and http://www.verbeia.com/mathematica/tips/HTMLLinks/Tricks_Misc_4.html.

For example try this

x = NewObject[y -> 3];
function[x, 4]
x.function2[5]

You can go further if you want to emulate object inheritance using a package called InheritRules available here http://library.wolfram.com/infocenter/MathSource/671/

You could also store the function definition not in newObject but in a type symbol, so if NewObject returned type[newObject] instead of newObject you could define function and function2 like this outside of NewObject (and not inside) and have the same usage as before.

function[type[object_], x_] ^:= object["y"] + x;
type /: type[object_].function2[x_] := 2 object["y"]+x;

Use UpValues[type] to see that function and function2 are defined in the type symbol.

Further ideas about this last syntax are introduced here https://mathematica.stackexchange.com/a/999/66.

Improved version of SelectEquivalents

@rcollyer: Many thanks for bringing SelectEquivalents to the surface, it's an amazing function. Here is an improved version of SelectEquivalents listed above with more possibilities and using options, this makes it easier to use.

Options[SelectEquivalents] = 
   {
      TagElement->Identity,
      TransformElement->Identity,
      TransformResults->(#2&) (*#1=tag,#2 list of elements corresponding to tag*),
      MapLevel->1,
      TagPattern->_,
      FinalFunction->Identity
   };

SelectEquivalents[x_List,OptionsPattern[]] := 
   With[
      {
         tagElement=OptionValue@TagElement,
         transformElement=OptionValue@TransformElement,
         transformResults=OptionValue@TransformResults,
         mapLevel=OptionValue@MapLevel,
         tagPattern=OptionValue@TagPattern,
         finalFunction=OptionValue@FinalFunction
      }
      ,
      finalFunction[
         Reap[
            Map[
               Sow[
                  transformElement@#
                  ,
                  {tagElement@#}
               ]&
               , 
               x
               , 
               {mapLevel}
            ] 
            , 
            tagPattern
            , 
            transformResults
         ][[2]]
      ]
   ];

Here are examples of how this version can be used:

Using Mathematica Gather/Collect properly

How would you do a PivotTable function in Mathematica?

Mathematica fast 2D binning algorithm

Internal`Bag

Daniel Lichtblau describes here an interesting internal data structure for growing lists.

Implementing a Quadtree in Mathematica

Debugging functions

These two posts point to useful functions for debugging:

How to debug when writting small or big codes using Mathematica? workbench? mma debugger? or something else? (ShowIt)

https://stackoverflow.com/questions/5459735/the-clearest-way-to-represent-mathematicas-evaluation-sequence/5527117#5527117 (TraceView)

Here's another function based on Reap and Sow to extract expressions from different parts of a program and store them in a symbol.

SetAttributes[ReapTags,HoldFirst];
ReapTags[expr_]:=
   Module[{elements},
      Reap[expr,_,(elements[#1]=#2/.{x_}:>x)&];
      elements
   ];

Here's an example

ftest[]:=((*some code*)Sow[1,"x"];(*some code*)Sow[2,"x"];(*some code*)Sow[3,"y"]);
s=ReapTags[ftest[]];
Keys[s]
s["x"]
PrintSymbol[s] (*Keys and PrintSymbol are defined above*)

Other resources

Here's a list of interesting links for learning purpose:

A collection of Mathematica learning resources

Updated here: https://mathematica.stackexchange.com/a/259/66

Answer 3

One of the things that bothers me about the built-in scoping constructs is that they evaluate all of the local variable definitions at once, so you can't write for example

With[{a = 5, b = 2 * a},
    ...
]

So a while ago I came up with a macro called WithNest that allows you to do this. I find it handy, since it lets you keep variable bindings local without having to do something like

Module[{a = 5,b},
    b = 2 * a;
    ...
]

In the end, the best way I could find to do this was by using a special symbol to make it easier to recurse over the list of bindings, and I put the definition into its own package to keep this symbol hidden. Maybe someone has a simpler solution to this problem?

If you want to try it out, put the following into a file called Scoping.m:

BeginPackage["Scoping`"];

WithNest::usage=
"WithNest[{var1=val1,var2=val2,...},body] works just like With, except that
values are evaluated in order and later values have access to earlier ones.
For example, val2 can use var1 in its definition.";

Begin["`Private`"];

(* Set up a custom symbol that works just like Hold. *)
SetAttributes[WithNestHold,HoldAll];

(* The user-facing call.  Give a list of bindings and a body that's not
our custom symbol, and we start a recursive call by using the custom
symbol. *)
WithNest[bindings_List,body:Except[_WithNestHold]]:=
WithNest[bindings,WithNestHold[body]];

(* Base case of recursive definition *)
WithNest[{},WithNestHold[body_]]:=body;

WithNest[{bindings___,a_},WithNestHold[body_]]:=
WithNest[
{bindings},
WithNestHold[With[List@a,body]]];

SyntaxInformation[WithNest]={"ArgumentsPattern"->{{__},_}};
SetAttributes[WithNest,{HoldAll,Protected}];

End[];

EndPackage[];

Answer 4

This is recipe from Stan Wagon's book...use it when built-in Plot behaves erratically due to lack of precision

Options[PrecisePlot] = {PrecisionGoal -> 6};
PrecisePlot[f_, {x_, a_, b_}, opts___] := Module[{g, pg},
   pg = PrecisionGoal /. {opts} /. Options[PrecisePlot];
   SetAttributes[g, NumericFunction];
   g[z_?InexactNumberQ] := Evaluate[f /. x -> z];
   Plot[N[g[SetPrecision[y, \[Infinity]]], pg], {y, a, b},
    Evaluate[Sequence @@ FilterRules[{opts}, Options[Plot]]]]];

I often use the following trick from Kristjan Kannike's when I need "dictionary-like" behavior from Mathematica's downvalues

index[downvalue_, 
   dict_] := (downvalue[[1]] /. HoldPattern[dict[x_]] -> x) // 
   ReleaseHold;
value[downvalue_] := downvalue[[-1]];
indices[dict_] := 
  Map[#[[1]] /. {HoldPattern[dict[x_]] -> x} &, DownValues[dict]] // 
   ReleaseHold;
values[dict_] := Map[#[[-1]] &, DownValues[dict]];
items[dict_] := Map[{index[#, dict], value[#]} &, DownValues[dict]];
indexQ[dict_, index_] := 
  If[MatchQ[dict[index], HoldPattern[dict[index]]], False, True];

(* Usage example: *)
(* Count number of times each subexpression occurs in an expression *)
expr = Cos[x + Cos[Cos[x] + Sin[x]]] + Cos[Cos[x] + Sin[x]]
Map[(counts[#] = If[indexQ[counts, #], counts[#] + 1, 1]; #) &, expr, Infinity];
items[counts]

When evaluation results are confusing, sometimes it helps to dump evaluation steps into a text file

SetAttributes[recordSteps, HoldAll];
recordSteps[expr_] :=
 Block[{$Output = List@OpenWrite["~/temp/msgStream.m"]}, 
  TracePrint[Unevaluated[expr], _?(FreeQ[#, Off] &), 
   TraceInternal -> True];
  Close /@ $Output;
  Thread[Union@
    Cases[ReadList["~/temp/msgStream.m", HoldComplete[Expression]], 
     symb_Symbol /; 
       AtomQ@Unevaluated@symb && 
        Context@Unevaluated@symb === "System`" :> 
      HoldComplete@symb, {0, Infinity}, Heads -> True], HoldComplete]
  ]

(* Usage example: *)
(* puts steps of evaluation of 1+2+Sin[5]) into ~/temp/msgStream.m *)
recordSteps[1+2+Sin[5]]

Answer 5

Remember that The Mathematica Book is also available online at http://reference.wolfram.com/legacy/v5_2/ - though it's superseded by the current documentation at http://reference.wolfram.com

Answer 6

I've mentioned this before, but the tool I find most useful is an application of Reap and Sow which mimics/extends the behavior of GatherBy:

SelectEquivalents[x_List,f_:Identity, g_:Identity, h_:(#2&)]:=
   Reap[Sow[g[#],{f[#]}]&/@x, _, h][[2]];

This allows me to group lists by any criteria and transform them in the process. The way it works is that a criteria function (f) tags each item in the list, each item is then transformed by a second supplied function (g), and the specific output is controlled by a third function (h). The function h accepts two arguments: a tag and a list of the collected items that have that tag. The items retain their original order, so if you set h = #1& then you get an unsorted Union, like in the examples for Reap. But, it can be used for secondary processing.

As an example of its utility, I've been working with Wannier90 which outputs the spatially dependent Hamiltonian into a file where each line is a different element in the matrix, as follows

rx ry rz i j Re[Hij] Im[Hij]

To turn that list into a set of matrices, I gathered up all sublists that contain the same coordinate, turned the element information into a rule (i.e. {i,j}-> Re[Hij]+I Im[Hij]), and then turned the collected rules into a SparseArray all with the one liner:

SelectEquivalents[hamlst, 
      #[[;; 3]] &, 
      #[[{4, 5}]] -> (Complex @@ #[[6 ;;]]) &, 
      {#1, SparseArray[#2]} &]

Honestly, this is my Swiss Army Knife, and it makes complex things very simple. Most of my other tools are somewhat domain specific, so I'll probably not post them. However, most, if not all, of them reference SelectEquivalents.

Edit: it doesn't completely mimic GatherBy in that it cannot group multiple levels of the expression as simply as GatherBy can. However, Map works just fine for most of what I need.

Example: @Yaroslav Bulatov has asked for a self-contained example. Here's one from my research that has been greatly simplified. So, let's say we have a set of points in a plane

In[1] := pts = {{-1, -1, 0}, {-1, 0, 0}, {-1, 1, 0}, {0, -1, 0}, {0, 0, 0}, 
 {0, 1, 0}, {1, -1, 0}, {1, 0, 0}, {1, 1, 0}}

and we'd like to reduce the number of points by a set of symmetry operations. (For the curious, we are generating the little group of each point.) For this example, let's use a four fold rotation axis about the z-axis

In[2] := rots = RotationTransform[#, {0, 0, 1}] & /@ (Pi/2 Range[0, 3]);

Using SelectEquivalents we can group the points that produce the same set of images under these operations, i.e. they're equivalent, using the following

In[3] := SelectEquivalents[ pts, Union[Through[rots[#] ] ]& ] (*<-- Note Union*)
Out[3]:= {{{-1, -1, 0}, {-1, 1, 0}, {1, -1, 0}, {1, 1, 0}},
          {{-1, 0, 0}, {0, -1, 0}, {0, 1, 0}, {1, 0, 0}},
          {{0,0,0}}}

which produces 3 sublists containing the equivalent points. (Note, Union is absolutely vital here as it ensures that the same image is produced by each point. Originally, I used Sort, but if a point lies on a symmetry axis, it is invariant under the rotation about that axis giving an extra image of itself. So, Union eliminates these extra images. Also, GatherBy would produce the same result.) In this case, the points are already in a form that I will use, but I only need a representative point from each grouping and I'd like a count of the equivalent points. Since, I don't need to transform each point, I use the Identity function in the second position. For the third function, we need to be careful. The first argument passed to it will be the images of the points under the rotations which for the point {0,0,0} is a list of four identical elements, and using it would throw off the count. However, the second argument is just a list of all the elements that have that tag, so it will only contain {0,0,0}. In code,

In[4] := SelectEquivalents[pts,  
             Union[Through[rots[#]]]&, #&, {#2[[1]], Length[#2]}& ]
Out[4]:= {{{-1, -1, 0}, 4}, {{-1, 0, 0}, 4}, {{0, 0, 0}, 1}}

Note, this last step can just as easily be accomplished by

In[5] := {#[[1]], Length[#]}& /@ Out[3]

But, it is easy with this and the less complete example above to see how very complex transformations are possible with a minimum of code.