How and why?
Mathematica Techniques
By Allan Hayes

For this issue, I look at targetting action at parts of expression, pitfalls from delayed evaluation, parallel assignments using the listability of Set and SetDelayed, updating data lists, wire frame and white surfaces, stacked graphics, and the adaptive sampling used by Plot and ParametricPlot.

Targetting where to act, Parts and Levels

In Mathematica programming we often target action at particular parts of an expression.  We can do this by using pattern matching, or, more precisely, by using the addressing, or indexing, of parts that is made possible by the regular structure of expressions. Here I shall look at using indices.

Expressions

All Mathematica expressions are of two kinds:

Atoms		symbol		a,a4b,$ab 
			integer 			0,2,-2 
			rational		2/3,-2/3 
			real			3.141, -12.8 
			complex		2 + 3/4 I 		( I = (Sqrt[-1]))
			string			"strings are in quotes" 

	Normal Expressions	
	head[element1, element2,..] where head and elements are expressions;  
	for example f[a,b][ a, g[1.5,2.4] ]

Because of the special forms allowed for convenience this is not always apparent, but we can always see it by using FullForm:

FullForm[a+b]

Sometimes we have to prevent evaluation

{FullForm[1+2], HoldForm[FullForm[1+2]]}

Indices

Each part of  an expression has at least one index, a sequence of integers specifying how to get to the part from outside the expression. For clarity a sequence like 1, 2, 3 will often be written in parentheses, (1, 2, 3).

h[ a, b, c] 
(indices L to R)		0  1  2  3
(indices R to L)	  	  -3 -2 -1		(-4 is not an index for h)

indices of marked part
h[a,b][c,d[e,a]]	(2), (-1)
h[a,b][c,d[e,a]]   	(0,1),(0,-2)				
h[a,b][c,d[e,a]]   	(2,2),(2,-1),(-1,2),(-1,-1)	
h[a,b][c,d[e,a]]   	()

The index of the whole expression is the empty sequence — no instructions, so stay where you are.

Part specification
used by: Extract, MapAt, ReplacePart, Delete, Insert, ReplacePart (but not by Part)

A part specification allows us to specify a single part or several parts (with repetitions allowed) as the target. To keep indices separate it is often necessary to wrap each one in list brackets.

Part specifications are of three kinds:

({{a list of indices each wrapped in list brackets, {{1},  {3, 2},  {-1, 0}}}, {a single index wrapped in list brackets, {3, 2};  {}}, {a single integer, 3}})

Examples

expr = h[a, b][c,k[d, p[e]]];

MapAt[F,expr, 1]		(*map at part 1*)

MapAt[F,expr, {2,1}]     (*map at part 2,1*)

MapAt[F,expr, {{2},{1}}] (*map at parts 2 and 1*)

MapAt[F,expr, {{0,2},{2,2},{},{}}]
(*map at parts 0,2 and 2,2 and twice at the whole expression*)

Position gives part specifications for the occurences of a pattern.

Position[expr, c|_p] (*c|_p = Alternatives[c,_p]*)

Oddly enough Part does not use part specifications as defined above:

Part[h[a,k[b,c]],{{2},{1}}]

But we do have

Part[h[a,k[b,c]],2,1]

Part[h[a,b,c],{2,2,1}]

Also, Part uses Cartesian products very conveniently for extracting submatrices and resetting parts:

(mat= Table[i+j,{i,4},{j,4}])//MatrixForm

Part[mat,{4,2},{1,3,1}]//MatrixForm

mat[[{2,3},{1,4}]] = {{A,B},{C,D}};
mat//MatrixForm

Levels specifications
used by Map, MapAll, MapIndexed, MemberQ, Scan, FreeQ, Position, Apply, Cases, Count, DeleteCases, Level

A convenient way to specify the target uses the level of a part, that is the length of its index. For example:

Map[F,h[p[a],q[b[c]]]] 
(*map at parts at level 1,excluding heads*)

Map[F,h[p[a],q[b[c]]],{2}] 
(*map at parts at level 2,excluding heads*)

Heads of expressions are rather special, and Mathematica allows us to specify whether we want to include them or not (a part is in a head if its index includes 0)

Map[F,h[p[a],q[b[c]]], Heads -> True]
(*map at parts at level 1,including heads*)

Map[F,h[p[a],q[b[c]]], {1,3}]
(*map at parts at levels 1 to 3, excluding parts of heads*)

Suppose that we want to map F at all atoms in h[p[a],q[b[c]]]. The problem is that they are at different levels; however it can be done, and we can include or exclude heads.

Map[F,h[p[a],q[b[c]]], {-1}]

Map[F,h[p[a],q[b[c]]], {-1}, Heads->True]

Here we are using the depth of an expression: Depth[expr] = 1 + the maximum index length for the parts of expr outside heads (so, atoms are expressions with depth 1). There is no option Heads ->True or False to let us include parts of heads.

Map[F,h[p[a],q[b[c]]], {-2}] 
(*map at parts with depth 2, excluding parts of heads*)

Map[F,h[p[a],q[b[c]]], {-2,-1}]
(*map at parts with depth between 1 and 2, excluding parts of heads*)

The general level specification allows depth and level to be used together.

A level specification, is of one of the three forms  
	{m}, {m}, {m, n}, where m, n  are integers, Infinity or -Infinity
	(perhaps followed by  Heads->True, or Heads->False). 
These are interpreted as follows:

({{ , i ≥ 0, i < 0}, {{i,    }, level ≥ i, depth ≤ -i}, {{ , i}, level ≤ i, depth ≥ -i }, {{i}, level = i, depth = -i}, {i, 0 < level ≤ i, depth ≥ -i}})

Heads -> False excludes parts of heads (parts with zero in their index)
Heads -> True allows all parts.

The default settings for functions that have the option Heads are:

({{\(FreeQ, \ Position	\), \(Heads\  -> True	\)}, {\(Apply, Cases, Count, DeleteCases, Level, 		
Map, 
              MapAll, MapIndexed, MemberQ, Scan\), \(Heads\  -> False	\)}})

Use all features of functions

In addition to using part and level specifications it is a good idea to look up the full information in the Help browser from time to time for other controls. For example, with Position , if we  want to find the first two positions of integers, in level 2 and not in a head, then we can use

Position[p[3][h[2.0,2],p[4]], _Integer, {2}, 2, Heads->False]

Cases has even more control features, for example

Cases[p[3][h[2.0,2],p[4]], n_Integer->S[n^2], {2}, 2, Heads->True]

Pitfalls from delayed evaluation and iterators

(1) What is going on below?

Table[D[Sin[x],x],{x,2,3}];

The explanation is that Table has the attribute HoldAll, so it is free to evaluate its entries in a special way: it evaluates first,  x = 2; D[Sin[x],x] giving D[Sin[2],2], hence the first message; then x = 3; D[Sin[x],x], which gives D[Sin[3],3] and causes the second message.

The solution is to to force the evaluation of D[Sin[x],x] to occur before the numerical values are assigned:

Table[Evaluate[D[Sin[x],x]],{x,1,2}]

This works because it is equivalent to

Table[Cos,{x,1,2}]

The same problem and solution occur with  Do, NIntegrate, NProduct, NSum, Plot, ParametricPlot, Plot3D, ParametricPlot3D, Product, Sum.

(2) Plot can plot a single expression or a nested list of expressions, so why do we get the following?

le = {x,{x^2}};

Plot[le,{x,0,1}];

The explanation is as follows. 
With Plot[expr, {x,a,b}], the unevaluated expr is inspected:  
• if  expr is not a list then Mathematica tries to construct a line for it on the assumption that for any real number, r, x=r;expr will evaluate to a real number (this was  not the case in the example above);
• if expr is a nested list, then it is flattened, to {(Subscript[expr, 1]),(Subscript[expr, 2]), . .} say; the lines for (Subscript[expr, 1]), (Subscript[expr, 2]) ... are constructed where possible, and, finally,  the lines that have been made are displayed together

The solution for the example is to force le to evaluate to {x,{x^2}} before the inspection is made

Plot[Evaluate[le],{x,0,1}];

The same problem and solution occur with ParametricPlot, and ParametricPlot3D.

Please note that the flattening is done on the unevaluated expression, so we get:

Plot[{Sin[6x],le},{x,0,1}];

Parallel assignment

It is sometimes convenient to use the listability of Set and SetDelayed:

{f[x_/;x>2],{g[x_]/;x>3,h[y_],z}} = {1,{2,3/;y>4,{4}}};

?f g h z

The crucial point is that the list structure on right side must be an extension of that on the left.
Here is an application of this technique within a program.

x={x0, x1, x2, x3};
e={e1,e2,e3};
Evaluate[Partition[x,2,1,p[x_/;#1<=x<=#2]&]] = e;

?p

Updating data: Hash and Dispatch

The following problem was posted to com.soft-sys.math.mathematica by Joel Cannon.
Given two lists of triplets, {x,y,z}, called rough and fine;  then for each pair {x,y,r} in rough and {x,y,f}in fine replace {x,y,r}with {x,y,f}. There is some ambiguity, for example from

rough = {{1,2,3},{3,1,6},{2,2,5},{1,2,1}};
fine = {{2,2,4},{1,2,0},{0,0,3},{1,2,9}};

we might get

{{1, 2, 0}, {3, 1, 6}, {2, 2, 4}, {1, 2, 0}}

or

{{1, 2, 9}, {3, 1, 6}, {2, 2, 4}, {1, 2, 9}}

Carl Woll gave the following solution, which from each triple {x,y,z) in fine gets an assignment g[x,y] = z. He noted that g[x,y]=z can be retrieved using hashing methods, which is not true for definitions that involve patterns like g[x,y,_]:={x,y,z}.

repair[1][rough_, fine_] := 
Module[{f,g,h},
	f[{x_,y_,z_}] := g[x,y] = z;
	f /@ fine;
	_g := Sequence[]; (*default for g*)
	h := {#1,#2,First[{g[#1,#2],#3}]}&;
	Apply[h,rough,{1}]
]

repair[1][rough, fine]

Since assignments are converted to replacement rules internally, and since there is a function Dispatch to speed up the application of rules, it seems reasonable to try to convert this code into one using replacement rules (the new head k in the following code makes sure that replacement takes place just where required).

repair[2][rough_,fine_]:=
Module[{r,k,R,F},
	{R ,F} = Apply[k,{rough,fine},{2}];
	r = Dispatch[hh=F/.k[x___,z_]:> (k[x]:>z)];
	R/. k[x___,z_]:>{x, k[x]/.r/._k:>z
	}
]

hh

repair[2][rough, fine]

Notice that repair[2] gives {1,2,0} where repair[1] gave {1,2,9}. This is because when repair[1] stores assignments, later assignments to the same list replace earlier ones; but when repair[2] builds a list of replacement rules, all the rules are kept and, for each replacement, only the first applicable one is used.

Here is the same code without Dispatch.

repair[3][rough_,fine_]:=
Module[{r,k,R,F},
	{R ,F} = Apply[k,{rough,fine},{2}];
	r = F/.k[x___,z_]:> (k[x]:>z);
	R/. k[x___,z_]:>{x, k[x]/.r/._k:>z
	}
]

To compare speeds I'll  use the test supplied by Carl.

Do[
R[n]= Flatten[Table[{i,j,Random[Real,89]},{i,n},{j,n}],1];
F[n]= Flatten[Table[{i,j,Random[Real,99]},{i,20},{j,20}],1],
{n,2,20,2}
];

TableForm[
Table[
( Do[repair[i][R[n],F[n]],{10}];//Timing//First)/Second,
{n,2,20,2},
{i,1,3}
],
TableHeadings->{Table[n,{n,2,20,2}], Table[repair[i],{i,1,3}]}
]

Wire frames and white surfaces with grid lines

These forms are sometimes preferable for printing

Plot3D as a white surface with grid lines

Plot3D[Sin[x+Sin[y]],{x,0,2Pi},{y,0,2Pi}, Shading->False];

Plot3D, wire frame

Plot3D[Sin[x+Sin[y]],{x,0,2Pi},{y,0,2Pi},HiddenSurface ->False];

ParametricPlot3D, white surface with grid lines

ParametricPlot3D[{x,y,Sin[x+Sin[y]]},{x,0,2Pi},{y,0,2Pi}, Shading->False];

ParametricPlot3D, wire frame

<<Graphics`Shapes`
ParametricPlot3D[{x,y,Sin[x+Sin[y]]},{x,0,2Pi},{y,0,2Pi},
DisplayFunction->Identity];
Show[WireFrame[%],DisplayFunction -> $DisplayFunction];

Stacked Graphics

Here is a rather unusual use of ParametricPlot3D:

f[x_, t_] := Exp[-(x - Cos[Pi t])^2]

gtab = 
Table[{
	ParametricPlot3D[
		{x, t, f[x, t] s, {Hue[t/5, 0.5, 1], EdgeForm[]}}, 
		{x, -5, 5}, {s, 0, 1}, 
		PlotPoints -> {50, 2}, (*only end points in s direction*)
      	Lighting -> False, DisplayFunction -> Identity], 
 	ParametricPlot3D[
 		{{x, t, f[x, t]},{x,t,0}}, {x, -5, 5}, 
      	PlotPoints -> 50, DisplayFunction -> Identity]},
    {t, 0, 3, 0.25}
];

Show[gtab,  DisplayFunction -> $DisplayFunction, 
      BoxRatios -> {1, 1, 1}, PlotRange  -> All];

Plot: adaptive sampling

Paul Abbot sent me a note on "No more aliasing" in Volume 7, No.1,  pointing  out that the answer is "No". Here is an adaptation of the example that he gave.

Plot[{Sin[x], x+Sin[x]},{x,0,24Pi}, PlotRange ->All];

This turns out to be because the bending of the graph is used in deciding whether to add extra sample points; and as measured for Plot,  x + Sin[x] bends less than Sin[x]. In trying to understand this I have arrived at the following account of how Plot evaluates (it seems to be substantially correct but please see the example). ParametricPlot evaluates simililarly.

Suppose we are evaluating Plot[f[x],{x,a,b}] with options 
	MaxBend -> mb 		(default 10)
	PlotPoints  ->  pp 	(default 25), 
	PlotDivision -> pd 	(default 30) .

To start with, a sample of pp numbers in (a,b) is constructed: the first number is slightly larger then a, the last one is slightly smaller than b, and they are almost equally spaced. 
Then a process that may add extra numbers to the sample begins: this examines triples of successive sample numbers starting at the first triple.
Suppose that a triple p, q, r is being examined.
A crucial measure used in the examination is the normalised bend for p, q, r. This is the bend, measured in degrees,  in the line through {{p, f[p]},{q,f[q]},{r,f[r]}} after it has been scaled so as to just fit into a square. For example

Show[Graphics[{Line[{{0,0}, {4,2},{6,1}}]}], AspectRatio ->Automatic, PlotRange ->All, Frame->True, FrameTicks -> False];

is scaled to

Show[Graphics[{Line[{{0,0}, {4,2},{6,1}}]}], AspectRatio ->1, PlotRange ->All, Frame->True, FrameTicks -> False];

Adding extra sample points is limited by the requirement that they may not be added between consecutive sample points with distance apart less than 
	min = (b-a)/((pp-1)pd).

In detail, when the triple p, q, r is being examined, the action is given by (adapting the Mathematica function Which)

({{Which, Action}, {(q - p <min  and  r - q< min
 or  normalised  bend < mb)  
and s is the next sample after r , leave p, q, r ;
examine q, r, s.}, {q - p ≥ r-q and q-p ≥ rd, 
, insert m = (p + q)/2 in the sample;
 examine p, m, q.}, {q - p < r-q  and r-q ≥ rd
, insert m = (q + r)/2 in the sample;
examine p, q, m}})

Here is some code that displays and lets us animate the steps; the output is the list of corresponding graphics objects,

PlotSteps[expr_, {x_,a_, b_}, opts___]:=
Module[{xs = {},ps, gr,pr},
	Needs["Utilities`FilterOptions`"];
	Print[Definition[Utilities`FilterOptions]];
	xs={};
	gr=Plot[AppendTo[xs,{x, expr}];expr,{x,0,2}, 		
	PlotRange -> {All,All},
	DisplayFunction -> Identity, opts];
	pr = FullOptions[gr, PlotRange];
	Table[
		ListPlot[
			ps=Sort[Take[xs, t]],	
			
			PlotJoined ->True,PlotRange->pr,
			Epilog-> {PointSize[.02],Point/@ps}
		],
	{t,2,Length[xs]}
	]
]

I'll show the steps in a graphics array, using the following function to convert the list to a table. Please note that replacement must be used to insert the option Ticks -> None into each graphic object in  ls; we cannot simply use it as an option to GraphicsArray. Also, the list lst has to be padded to make its length a multiple of n, because Partition discards any remainder that is not of length n .

ToGraphicsArray[lst_, n_]:=
GraphicsArray[Partition[Join[lst/.
Graphics[g_,ops___]-> Graphics[g, Ticks->None,ops]
, Table[{Graphics[{}]}, {Mod[Length[lst],n]}]],n]];

First, extract the output from PlotSteps, using Block[{$DisplayFunction=Identity} , ...] to temporarily turn off the display.

steps = 
Block[{$DisplayFunction=Identity},
PlotSteps[Sin[Pi x], {x,0, 2}, PlotPoints ->6,
PlotDivision ->2]];

Now show the array

Show[ToGraphicsArray[steps,4]];

Does this help explain the difference between the plots of  Sin[x] and x + Sin[x] that we started with? The clue is the effect on the normalized  bend. Adding x to Sin[x] corresponds to changing this

Show[Graphics[{Line[{{0,0}, {4,2},{6,1}}]}], AspectRatio ->1, PlotRange ->All, Frame->True, FrameTicks -> False];

into to this

Show[Graphics[{Line[{{0,0}, {4,6},{6,7}}]}], AspectRatio ->1, PlotRange ->All, Frame->True, FrameTicks -> False];

Giving a marked reduction in the normalized bend, which would in turn reduce the likelihood of adding extra sample points. Perhaps someone can come up with an algebraic analysis.

Note: there is still a little difference between what happens and what is predicted by the preceding account.
We would expect the following to make a plot with five sample points, but it gives seven.

Plot[x^2, {x, 0, 2}, PlotPoints ->5, PlotRange->All, PlotDivision->1, MaxBend -> 200];

Cases[%, Line[pts_]:> pts, Infinity][[1]]