How and why?
Mathematica Techniques
By Allan Hayes

Some tips on speeding up code; how Plot3D and SurfaceGraphics work; getting packages to load automatically.

Speedups

Faster repeated matrix multiplication

Manipulations such as simplification takes time and memory, but sometimes they are needed to avoid large expressions building up and unnecessary computation. For example, in the last of the following inputs, taken from a posting to the newsgroup sci.math.num-analysis, seventeen 4 × 4 matrices are multiplied together: it takes 209 seconds. However by following every matrix multiplication with Simplify the answer is found in about 16 seconds; and if we use Expand instead of Simplify the answer is given in about one and a half seconds. The technique is generally applicable.

y = {{1/2, 1/2, 1/2, 1/2}, {-(1/2), 1/2, -(1/2), 1/2}, {-(1/2), -(1/2), 1/2, 1/2}, 
    {1/2, -(1/2), -(1/2), 1/2}}; 
y1 = Transpose[y];

t = {{Exp[(I s u)/8], 0, 0, 0}, 
    {0, 1/2 + 1/2 Exp[-I/4 s u], 
     -(1/2) + 1/2 Exp[-I/4 s u], 0}, 
    {0, -(1/2) + 1/2 Exp[-I/4 s u], 
     1/2 + 1/2   Exp[-I/4 s u], 0}, 
    {0, 0, 0, Exp[(I s u)/8]}};

\(t1\  = \ Simplify[t\  . \ t]; 
t2\  = \ Simplify[t1\  . \ t1]; \)

lst
        = {t, y, t2, y, t1, y, t2, y, t1, y1, t2, y1, t1, y1, t2, y1, t};

(e = Dot @@ lst;) // Timing

The problem is that Dot does not automatically simplify its output. This allows very large expressions to build up. When simplification is included we get the answer in about 16 seconds

(e1 = Fold[Simplify[#1 . #2]&, First[lst], Rest[lst]]) // Timing

For more flexibility we can define a variant of Dot that allows us to apply any function after Dot; and, using Expand instead of Simplify, the answer is given in less than two seconds

Dot2[e_][x_, y___] := Fold[e[Dot[#1, #2]]&,  x,  {y}]

(e2 = Dot2[Expand]@@lst);  // Timing

e1  == e2 // Simplify

This can also be done by giving Dot an option; and, if we wish, we can use the attributes, Flat and OneIdentity instead of the function Fold.

\(Unprotect[Dot]; 
Options[Dot]\  = \ {PostProcess -> Identity}; 

    Dot[x___, \ opts__?OptionQ] := 

      Module[{d, p = \ \(PostProcess /. {opts}\) /. Options[Dot]}, \ 
        Attributes[d] = \ {Flat, \ OneIdentity}; \ 
        d[u_, v_] := \ p[Dot[u, v]]; 
		d[x]
	]\ ; 
Protect[Dot]; \)

Dot[Sequence@@lst,  PostProcess -> Expand] // Timing

Quick and simple Together

Here is a quick, direct and simple version ot Together; it arose from a posting to comp.soft-sys.math.mathematica about simplifying large expressions with leaf counts of around 10 million. It seems that because Together expands denominators it makes the problem worse.

\(SimpleTogether[expr_Plus] := 
	
      \(\((\((Plus@@\((\ #\ List@@expr)\))\)/\ #)\)&\)[
        Times@@Union[Denominator/@\(List@@expr\)]]\)SimpleTogether[expr_] :=  expr

Ouputs

expr = 
        (a + b + c) (p + q)^3/(c + d/g) + u/p +  v/p  + q;

SimpleTogether[expr]

MapAll (//@) takes SimpleTogether inside expr.

SimpleTogether //@ expr

Together[expr]

Timings

Do[SimpleTogether[expr], {100}] // Timing // First

Do[SimpleTogether //@ expr, {100}] // Timing // First

Do[Together[expr], {100}] // Timing // First

Rough and Ready Numerics

Occasionally, we just want to get a quick estimate. Some functions have an option PrecisionGoal which we can set low so that that we get a more or less guaranteed precision quickly

NIntegrate[Sin[x + y], {x, 0, 7}, {y, 0, 7}] // Timing

NIntegrate[Sin[x + y], {x, 0, 7}, 
        {y, 0, 7}, 
PrecisionGoal -> 3] // Timing

In some cases we can restrict the number of significant digits used in internal calculations.

Block[{$MaxPrecision = 3}, 
        
NIntegrate[Sin[x + y], {x, 0, 7}, {y, 0, 7}]
          
] // Timing

Nested Lists and Flattening better than Nested Append

\(data1\  = \ {}; 
Do[data1 = Append[data1, n], {n, 1000}] // Timing\)

\(\(\((data2\  = \ {}; 
	Do[data2\  = \ {data2, n}, {n, 1000}]; 
          data2 = \ Flatten[data2])\);\) // Timing\)

data1  ==  data2

The main difference is that during construction data1 = Append[data1,n] constructs a new list at each step whereas with data2 = {data2, n}, the first place simply  points to the old data2 and no reconstruction is done until the final flattening.

Transpose versus Mapping

Transpose is another function that seems to do a lot of work reconstructing an object, but in reality just reindexes it. It was asked on comp.soft-sys.math.mathematica if there was a neater and faster way of constructing
	{{x1,y1,z1},{x2,y2,z2},{x3,y3,z3},{x4,y4,z4}} 
from

\(xydata = {{x1, y1}, {x2, y2}, {x3, y3}, {x4, y4}}; 

    zdata = {z1, z2, z3, z4}; \)

than (essentially)

xyzdata = Transpose[Append[Transpose[xydata], zdata]]

I rather unthinkingly posted the following variant:

Flatten /@ Transpose[{xydata, zdata}]

But this is very inefficient: I am grateful to Hartmut Wolf for reminding me of why this is so. Here are some timings.

\(xydata = Table[i + j, {i, 10000}, {j, 5}]; 

    zdata = \ Table[i, {i, 10000}]; \)

Transpose[Append[Transpose[xydata], zdata]]; // Timing

Flatten /@ Transpose[{xydata, zdata}]; // Timing

MapThread[Append, {xydata, zdata}]; // Timing

Plot3D, SurfaceGraphics

ParametricPlot3D, which we looked at in the last issue, creates a Graphics3D object containing lines or polygons, but Plot3D makes a SurfaceGraphics object containing an array of heights for the function being plotted together with appropriate options to guide Show in constructing the display. We shall see that polygons can be colored either by the option ColorFunction ( for Plot3D or SurfaceGraphics) or by an matrix of color directives in SurfaceGraphics.

\(sgr = 

      Plot3D[x^2 + y^2, \ {x, \(-2\), 1}, {y, 0, 2}, 

        PlotPoints\  -> \ {4, \ 3}, \ 
        AxesLabel -> {"x", "y", "z"}]\)

The array of heights is the first entry in ⁃SurfaceGraphics⁃

Short[InputForm[sgr], 4]

The heights are read matrix fashion — by increasing y top to bottom and increasing x left to right: but the proper x and y values have to be reconstructed by scaling using the option MeshRange -> {{-2,1},{0,2}}. 
Please note that MeshRange is an option for SurfaceGraphics that is inserted by Plot3D, it is not an  option for Plot3D.

(arr = sgr[[1]]) // MatrixForm

This orientation is the transpose of what we get with.

Table[x^2 + y^2,  {x, -2, 1}, {y, 0, 2}] // MatrixForm

I'll try to show what is going on by starting from a more interesting array of heights.

f[x_, y_] :=  
      1/( .001 + (y - 1)^6 + x^6) - 
        1/(( .001 + y + 1)^6 + x^6)  f[x_, y_] /; x^2 + y^2 <  .01  :=  "undefined";

hts  =  
      Transpose[Table[f[x, y], {x, -1, 1,   .2}, {y, -1, 1,  .2}]];

sgr1 = SurfaceGraphics[hts];

As menioned above, for an m by n array of z-values the information that x  and y range over {xm1, xm2} and {ym1, ym2} can be given by the option MeshRange -> {{xm1, xm2},{ym1,ym2}}  (the default is MeshRange ->{{1,m},{1,n}})

The evaluation of Show[sgr1, opts] may be modelled as follows:

The plot range  {{x1,x2},{y1,y2},{z1,z2}} is calculated

PlotRange[sgr1]

Certain clipping and scaling functions are defined:
The function clp clips numerical z-values to lie between z1 and z2.
The linear scaling functions xmsc, ymscl, zsc are defined by
	xmsc takes {1,m} to {xm1, xm2}
	ymsc takes {1,n} to {ym1, ym2}
	zsc takes {z1, z2} to {0, 1}.
For each i in {1,m}and j in {1,n} the point pt[i,j] is defined by
	pt[i,j] = {xmsc[i], ymsc[j], clp[h[i,j]]}.
Then the polygons  
Polygon[{pt[i,j],pt[i+1,j],pt[i+1,j+1],pt[i,j+1]}],  i<m, j<n,  
which have real coordinates and lie inside the plot range are constructed (the plot range condition comes into play if we put in our own values, see sgr4 below).

With ColorFunction -> Automatic (the default) coloring is by simulated lighting.

Notice that four polygons are missing because f[0,0] is not numeric.

sgr2 = 
      Show[sgr1, MeshRange  -> {{-1, 1}, {-1, 1}}, Axes -> True];

I'll turn off the messages:

Off[Plot3D::"gval"]

With ClipFill -> {cfbot, cftop} (default RGBColor[1,0,0]), if the heights of athe polygon's vertices had all been increased by clip then its color is cfbot ( if they  had all been decreased then the color is cftop).
With ColorFunction -> cf, the color of the other polygons is cf[average of the (scaled and clipped) z-values of vertices].

\(\(sgr3 = Show[sgr2, 	ClipFill -> \ {Hue[0, 0,  .6], Hue[0, 0, 1]}]; 
    \)\)

sgr4 = Show[sgr3, ColorFunction ->  (Hue[2/3 #]&)];

We can adjust the plot range, style the mesh (or remove it with Mesh -> False). Notice that a polygon is only displayed if all its vertices are inside the plot range.

\(\(Show[sgr4, \ 
      PlotRange\  -> {{\(- .99\),  .99}, All, {\(-300\), 300}}, 

      MeshStyle -> \ {Thickness[ .01], Hue[0]}]; \)\)

The coloring of the polygons can also be set by an  (m-1) by (n-1) shading matrix of color directives.

(shds  =  Table[Hue[(i + j)/5. ], {i, 2}, {j, 3}]) // 
      MatrixForm

Show[SurfaceGraphics[arr,  shds]];

The function Shades defined below constructs a shading matrix using the options ColorFunction and ClipFill.

\(Shades[hts_, \ opts___?OptionQ] := 

      Block[{dm, dn, cf, clpbot, clptop, sc, zmin, zmax, clip, hts4}, 
	
        pp = 	\({pr, cf, \ {clpbot, clptop}} = 
			
            \(\({PlotRange, ColorFunction, ClipFill} /. 

                  				{opts}\) /. {
                  PlotRange :> PlotRange[SurfaceGraphics[hts]]}\) /. 
              Options[SurfaceGraphics]\); 

        		{zmin, zmax} = \(zz = \ Last[pr]\); 

        		{dn, dm} = \ Dimensions[hts]; 
		
        sc = Interpolation[{{zmin, 0}, {zmax, 1}}, \ 
            InterpolationOrder -> 1]; 
		
        clip[z_] := Max[Min[z, zmax], zmin]; 
		
        hts4\  = 
		
          Table[{hts[\([n, m]\)], hts[\([n + 1, m]\)], 
              hts[\([n + 1, m + 1]\)], 
		hts[\([n, m + 1]\)]}, {n, 1, 
              dn - 1}, {m, 1, dm - 1}]; 
		
        hts4 /. {
		
            p : {__?NumericQ} :> \ 
			
              Which[\ \ \ 
				\ \ Max[p] < zmin, clpbot, 
					
                Min[p]\  > zmax, \ clptop, 
					True, 
                cf[ .25 \((Plus@@\((sc/@\(clip/@p\))\))\)]
					], 

            				{_?\((Not[ListQ[#]]&)\), _, _, _} -> \ Hue[0]			

            			}
		]\)

Construct the shading matrix corresponding to sgr4.

shds = 
Shades[hts, 
	ColorFunction -> (Hue[2/3*#1] & ), 
   	ClipFill -> {Hue[0, 0, 0.5], Hue[0, 0, 1]}
];

Try it out.

\(\(Show[SurfaceGraphics[hts, shds], 

      MeshStyle -> \ {Thickness[ .01], Hue[0]}, 

      MeshRange\  -> \ {{\(-1\), 1}, {\(-1\), 1}}, 
Axes -> True]; \)\)

Compare with the corresponding Plot3D display ( I have deleted a message arising from f being non-numeric at {0,0})

Plot3D[f[x, y], {x, -1, 1}, {y, -1, 1}, 
  	ClipFill -> {Hue[0, 0, 0.6], Hue[0, 0, 1]}, 
  	ColorFunction -> (Hue[2/3*#1]&), 
  	MeshStyle -> {Thickness[0.01], Hue[0]}, 
  	PlotPoints -> 11
];

Autoloading packages

The numbered parts of this section describe how to make
1. the package Graphics`FilledPlot`  load immediately the kernel is launched
2. the package Graphics`Arrow` load the first time that one of its functions is used
3. each  Statistics` package to  load the first time that one of its functions is used.

I'll describe firstly how to do these things using the folder Autoload that you will find at

ToFileName[{$TopDirectory, "AddOns",  "Autoload"}]

(the output may give a different location in your system)
Other possible folders to use are mentioned at the end of this section.

The methods work because any file named init.m that is inside either Autoload or a folder within it is automatically evaluated when the kernel is loaded (allowing init.m files  inside folders within  Autoload, permits several such files to be used, if desired)

For definiteness, open, or if necessary create, a file init.m inside a folder named Kernel that is itself inside the folder Autoload, for reference call this the loading file.

1. To have the package Graphics`FilledPlot`  load immediately the kernel is launched:
- type	<<Graphics`FilledPlot into the loading file
- save the file;
When the kernel starts the loading file is evaluated and the package is loaded just as if you had evaluated <<Graphics`FilledPlot` in a notebook. You can add loading instructions <<... for any number of packages to the loading file.

Check the above by quitting the kernel and restarting it (use menu Kernel).

$Packages

2. To have the package Graphics`Arrow` load the first time that one of its functions is used:
- open the file AddOns\StandardPackages\Graphics\Kernel\init
- copy the text
	(* Declarations for file Graphics`Arrow` *)
	DeclarePackage["Graphics`Arrow`",
     	{"Absolute", "Arrow", "HeadCenter", "HeadLength",
     	 "HeadScaling", "HeadShape", "HeadWidth", "Relative", "ZeroShape"}];
    - paste into the loading file;
    - save the file
 
When the kernel starts the expression DeclarePackage[...] is evaluated: this adds the context "Graphics`Arrow`" to $ContextPath and creates the symbols Absolute, Arrow,... in this context (their long names will be "Graphics`Arrow`Absolute" ...) and gives them the attribute Stub. This means that if one of the symbols, Arrow say, is entered then the context "Graphics`Arrow`" is searched; the symbol Graphics`Arrow`Arrow is found; and because it has the attribute Stub, Needs["Graphics`Arrow`"] is evaluated and the package is loaded, if not already loaded (evaluating any other of the symbols Absolute, Arrow, HeadCenter, ... has the same effect but, since Needs is used to load the package it will only be loaded once - a query like ?Arrow will not load the package and will not give the package usage information)

Check the above by quitting the kernel and restarting it (menu Kernel).

Asking for information with ?.. does not trigger loading

?Arrow

But we do get loading after ordinary evaluation:

Arrow;

?Arrow

Warning: if you create your own packages please note that to have mehod2 working in the background we must follow the convention that the package context should give the directions to the package file when ContextToFileName is used on it, thus:

ContextToFileName["Graphics`Arrow`"]

Method 2 will still work if you don't follow this convention but you may get a large number of messages about shadowing
For example: 
- temporarily change the name of the package file Arrow to MyArrow (its context, which is set in the package code by BeginPackage["Graphics`Arrow"], is still "Graphics`Arrow");
- in the text of the loading file, change Graphics`Arrow` to Graphics`MyArrow` (we kept  DeclarePackage["Graphics`Arrow`", {...}] then on entering Arrow, Mathematica would look for the name Graphics\Arrow.m. which is no longer used). 

There is still problem: on starting the kernel, it will be the context "Graphics`MyArrow`" that will be added to $ContextPath and the symbols, MyArrow` Absolute,..., that will be created (with attribute Stub). Now, if Arrow is entered, although the package Graphics`MyArrow` will be loaded,  because the package context is still "Graphics`Arrow`",  a new set of symbols, Graphics`Arrow`Absolute...  will be created, and this will cause warnings to be generated

Check the above by quitting the kernel and restarting it.

Arrow

-----------  many warnings deleted --------------

In spite of these, everything works

Show[Graphics[Arrow[{0, 0}, {1, 1}]]];

The reason that it works is that the context "Graphics`Arrow`" is placed earlier than the context "Graphics`MyArrow`" in $ContextPath:

$ContextPath

so that,  on evaluating Arrow[ ..] the symbol  Graphics`Arrow`Arrow is found and from it the rules that were generated for it when the package was loaded:

?Graphics`Arrow`Arrow

No rules were generated for the symbol Graphics`MyArrow`Arrow:

?Graphics`MyArrow`Arrow

Change Graphics`MyArrow` back to Graphics`Arrow`
 in the package name and in the loading file

3.. To have each  Statistics` package to  load the first time that one of its functions is used:- type into the loading file <<Statistics`

Check by quitting the kernel and restarting it.

?MovingAverage

MovingAverage

Delete the contents of the "loading file"

Alternative folders

Instead of working in the Autoload folder that we used above we could have used a folder at

ToFileName[{$PreferencesDirectory, "AddOns",  "Autoload"}]

or simply (but in this case we must definitely put the init.m file inside a folder called Kernel)

$PreferencesDirectory

With any of the locations your own files may be lost if Mathematica is re-installed (a copy of your files should be kept elsewhere).

It's tempting to have all packages load at startup, but beside taking up RAM and slowing startup, it runs the risk of producing shadowing of symbols - far better to arrange for loading if a symbol is used.