Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

What is in your Mathematica tool bag? [closed]

One of the nice things about the Mathematica notebook interface is that it can evaluate expressions in any language, not just Mathematica. As a simple example, consider creating a new Shell input cell type that passes the contained expression to the operating system shell for evaluation.

First, define a function that delegates evaluation of a textual command to the external shell:

shellEvaluate[cmd_, _] := Import["!"~~cmd, "Text"]

The second argument is needed and ignored for reasons that will become apparent later. Next, we want to create a new style called Shell:

  1. Open a new notebook.
  2. Select the menu item Format/Edit Stylesheet...
  3. In the dialog, beside Enter a style name: type Shell.
  4. Select the cell bracket beside the new style.
  5. Select the menu item Cell/Show Expression
  6. Overwrite the cell expression with the Step 6 Text given below.
  7. Once again, select the menu item Cell/Show Expression
  8. Close the dialog.

Use the following cell expression as the Step 6 Text:

Cell[StyleData["Shell"],
 CellFrame->{{0, 0}, {0.5, 0.5}},
 CellMargins->{{66, 4}, {0, 8}},
 Evaluatable->True,
 StripStyleOnPaste->True,
 CellEvaluationFunction->shellEvaluate,
 CellFrameLabels->{{None, "Shell"}, {None, None}},
 Hyphenation->False,
 AutoQuoteCharacters->{},
 PasteAutoQuoteCharacters->{},
 LanguageCategory->"Formula",
 ScriptLevel->1,
 MenuSortingValue->1800,
 FontFamily->"Courier"]

Most of this expression was copied directly form the built-in Program style. The key changes are these lines:

 Evaluatable->True,
 CellEvaluationFunction->shellEvaluate,
 CellFrameLabels->{{None, "Shell"}, {None, None}},

Evaluatable enables the SHIFT+ENTER functionality for the cell. Evaluation will call the CellEvaluationFunction passing the cell content and content type as arguments (shellEvaluate ignores the latter argument). CellFrameLabels is just a nicety that let's the user identify that this cell is unusual.

With all of this in place, we can now enter and evaluate a shell expression:

  1. In the notebook created in the steps above, create an empty cell and select the cell bracket.
  2. Select the menu item Format/Style/Shell.
  3. Type a valid operating system shell command into the cell (e.g. 'ls' on Unix or 'dir' on Windows).
  4. Press SHIFT+ENTER.

It is best to keep this defined style in a centrally located stylesheet. Furthermore, evaluation functions like shellEvaluate are best defined as stubs using DeclarePackage in init.m. The details of both of these activities are beyond the scope of this response.

With this functionality, one can create notebooks that contain input expressions in any syntax of interest. The evaluation function can be written in pure Mathematica, or delegate any or all parts of the evaluation to an external agency. Be aware that there are other hooks that relate to cell evaluation, like CellEpilog, CellProlog and CellDynamicExpression.

A common pattern involves writing the input expression text to a temporary file, compiling the file in some language, running the program and capturing the output for ultimate display in the output cell. There are plenty of details to address when implementing a full solution of this kind (like capturing error messages properly), but one must appreciate the fact that it is not only possible to do things like this, but practical.

On a personal note, it is features like this that makes the notebook interface the center of my programming universe.

Update

The following helper function is useful for creating such cells:

evaluatableCell[label_String, evaluationFunction_] :=
  ( CellPrint[
      TextCell[
        ""
      , "Program"
      , Evaluatable -> True
      , CellEvaluationFunction -> (evaluationFunction[#]&)
      , CellFrameLabels -> {{None, label}, {None, None}}
      , CellGroupingRules -> "InputGrouping"
      ]
    ]
  ; SelectionMove[EvaluationNotebook[], All, EvaluationCell]
  ; NotebookDelete[]
  ; SelectionMove[EvaluationNotebook[], Next, CellContents]
  )

It is used thus:

shellCell[] := evaluatableCell["shell", Import["!"~~#, "Text"] &]

Now, if shellCell[] is evaluated, the input cell will be deleted and replaced with a new input cell that evaluates its contents as a shell command.


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];

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.


This is not a complete resource, so I'm throwing it here in the answers section, but I have found it very useful when figuring out speed issues (which, unfortunately, is a large part of what Mathematica programming is about).

timeAvg[func_] := Module[
{x = 0, y = 0, timeLimit = 0.1, p, q, iterTimes = Power[10, Range[0, 10]]},
Catch[
 If[(x = First[Timing[(y++; Do[func, {#}]);]]) > timeLimit,
    Throw[{x, y}]
    ] & /@ iterTimes
 ] /. {p_, q_} :> p/iterTimes[[q]]
];
Attributes[timeAvg] = {HoldAll};

Usage is then simply timeAvg@funcYouWantToTest.

EDIT: Mr. Wizard has provided a simpler version that does away with Throw and Catch and is a bit easier to parse:

SetAttributes[timeAvg, HoldFirst]
timeAvg[func_] := Do[If[# > 0.3, Return[#/5^i]] & @@ 
                     Timing @ Do[func, {5^i}]
                     ,{i, 0, 15}]

EDIT: Here's a version from acl (taken from here):

timeIt::usage = "timeIt[expr] gives the time taken to execute expr, \
  repeating as many times as necessary to achieve a total time of 1s";

SetAttributes[timeIt, HoldAll]
timeIt[expr_] := Module[{t = Timing[expr;][[1]], tries = 1},
  While[t < 1., tries *= 2; t = Timing[Do[expr, {tries}];][[1]];]; 
  t/tries]