On the Goodness of Programming Languages

Since my last post on C I have been thinking a lot about what makes a programming language good. This is also motivated by my recent emersion into JavaScript that has me wondering how good of a language is JavaScript. I've had similar thoughts about C++, Java, Lisp and others in the past.

I think that the only sensible conclusion is there is no way to measure the goodness of a language in any objective way. I am certain that I can can come up with a metric that would declare C as the best language ever developed, and another that would make the winner Lisp, or Mathematica or JavaScript. For example if, to simplify things a bit, I decide that the metric for a language's goodness is that:

1. It be functional (first class functions, closures, etc.)
2. It be ubiquitous (available on almost every computing device)

Then it is clear that JavaScript is the best language ever. I can also easily derive a metric that would declare Mathematica as the best language ever (and on some days I am convinced of the absolute truth of this!).

But the truth is that there can never be the best language ever. This is not because programming language technology is constantly (or at least in fits and starts) evolving. The reason there can never be a best programming language begins with the idea of Turing Completeness but it does not end there. Any language worth learning is TC and by that metric is no better than any other. Yet, it is often argued that some languages are more expressive than others (e.g. Lisp is more expressive than Fortran) and therefore one can solve a wider class of problems faster in Lisp than Fortran. But that does not mean there are not valid metrics by which Fortran can be declared the world champion of programming languages.

Humans spend a way too much time debating what the best is. What is the best car, best restaurant, best sport, best OS and best programming language? We do this because we want to be associated with the best. I drive an X, X is the best, therefore I have superior taste. For programming languages such arguments are silly not because we each have our own tastes. They are silly because if you are a programmer and don't sometimes see the superiority of Language A over Language B (and visa versa) then my guess is you are a rank amateur (or really boring person to sit next to at the office). So next time you find yourself in one of those "my language is better than yours" debates (and I've been there too) stop yourself, pick up the best editor :-) and write some code instead!

The Elegance of Rule-Based Programming in Mathematica

Not too long ago I had to take a programming test as part of the application process for a development position at a hedge fund. The test had two problems - one easy and one not so easy. The solutions to both problems had to be implemented in your choice of Java, C or C++. You had two hours to finish.


I have never been good at presure cooker coding tests like this. I personally don't think these kind of tests bring out the best in most developers but there are some good devlopers who work well under these conditions. I just am not one of them.


In any case I recently thought of the hard problem on this test and how trivial it would be to solve in Mathematica. Indeed I whipped up and tested this solution on my train ride home which is only about 45 mins!


The solution is very short and I am guessing is much shorter than the average solution you can find in Java or C in the equivalent amount of time. The brevity come form exploiting mathematica's powerful rule-based approach. Most of the code is comments!

Problem

A collection of particles is contained in a linear chamber. They all have the same speed, but some are headed toward the right and others are headed toward the left. These particles can pass through each other without disturbing the motion of the particles, so all the particles will leave the chamber relatively quickly.


You will be given the initial conditions by a string containing at each position an 'L' for a leftward moving particle, an 'R' for a rightward moving particle, or a '.' for an empty location. Initially, no location in the chamber contains two particles passing through each other.


Create an animation of the process. At each unit of time, you want a string showing occupied locations with an 'X' and unoccupied locations with a '.'. Create code that for a function animate that is given an integer speed and a string giving the initial conditions. The speed is the number of positions each particle moves in one time unit.


The function will return a list of strings in which each successive element shows the occupied locations at the next time unit. The first element of the return should show the occupied locations at the initial instant (at time = 0) in the 'X', '.' format. The last element in the return should show the empty chamber at the first time that it becomes empty.

Solution

ClearAll[simulation, step, L, R, r, l, rule1, rule2, rule3, animation];

(* 
This function does most of the work. It simulates a single step in the animation.
It uses 3 rules which are initialized in the calling function and are visble
here via the dynamic scoping of Block.

There are 3 transformations performed within a Do. The Do's role is to run the
transformations velocity times as a means for making the particles move
that many steps.

The first transform takes the input and adds an empty cell to the start and end of the chamber.
This is done to avoid having special rules to deal with the 
boundary conditions at the end of the chamber. 

The second transformation uses rules that are applied repeatedly using ReplaceAll (//.). 
I discuss the rules below. 

Delete is used to remove the dummy empty cells added in the first step and 
a third transformation maps lower case r and l back to upper case before the Do
repeats.

*)
step[chamber_List, velocity_Integer] := Module[{work = chamber},
  Do[
   work = 
    Delete[work  /. {P__} :> {{}, P, {}} //. {rule1, rule2, 
        rule3}, {{-1}, {1}}] /. {r -> R, l -> L}, {velocity}];
  work
  ]

(*
The simulation sets up 3 rules. It is not necessary to define the rules here in a
Block. The code is written this way largely because I incrementally developed the
rules in the global scope and grafted them into a program later on.

rule1 is responsible for moving a R partcle to the Right. NOTE: When a particle moves
I change its case from upper to lower to make it invisble to subsequent rules.

rule2 is responsible for handling an intermediate condition where a right moving
particle lands in a cell of another right moving particle that has not moved yet.

rule3 deals with left moving particles being careful to also handle the case where a
previously moved partcle Z may be in the same cell as the L. 

The hardest aspect of this solution was distilling the transformations down to these
three rules. It took some trial and error. It is possible I missed some corner case 
so let me know if you see a bug!!

The idea here is that FixedPointList drives the simulation until it reaches a steady
state and Most removes the repetitive last entry. 

A final transformation maps particles to "X" as dicated by the specifications.
*)

simulation[chamber_List, velocity_Integer] := 
 Block[
  {rule1 = {X___, {R}, {P___}, Y___} :> {X, {}, {P, r}, Y} ,
   rule2 = {X___,  {R | r, r}, {P___}, Y___} :> {X, {r}, {P, r}, Y},
   rule3 = {X___, {P___}, {L, Z___}, Y___} :> {X, {l, P}, {Z}, Y} },
  Most[FixedPointList[
     step[#, velocity] &, chamber]] //. {{(L | R) ..} -> "X", {} -> 
     "."}
  ]

(* 
The main function does little but convert from the string encoding specified by the problem to a more convenent symbolic form where each position in the chamber is one of {L}, {R} or {} for empty and the chamber itself is a list rather than a string.
*)

animation[chamber_String, velocity_Integer] := 
 simulation[
  (Characters[chamber] //. {"." -> {}, a_String :> {Symbol[a]}}), 
  velocity]

The Test Cases

1) The single particle starts at the 3rd position, moves to the 5th, then 7th, and then out of the chamber.

In[360]:= animation["..R....",2]//Grid
Out[360]=
..X....
....X..
......X
.......

2) At time 1, there are actually 4 particles in the chamber, but two are passing through each other at the 4th position.

In[361]:= animation["RR..LRL",3]//Grid
Out[361]= 
XX..XXX
.X.XX..
X.....X
.......

3) At time 0 there are 8 particles. At time 1, there are still 6 particles, but only 4 positions are occupied since particles are passing through each other.

In[362]:= animation["LRLR.LRLR",2]//Grid
Out[362]=
XXXX.XXXX
X..X.X..X
.X.X.X.X.
.X.....X.
.........

4) These particles are moving so fast that they all exit the chamber by time 1.

In[363]:= animation["RLRLRLRLRL",10]//Grid
Out[363]= 
XXXXXXXXXX
..........

5) The empty chamber test

In[364]:= animation["...",1]//Grid
Out[364]=
...

6) A complicated test

In[365]:= animation["LRRL.LR.LRR.R.LRRL.",1]//Grid
Out[365]= 
XXXX.XX.XXX.X.XXXX.
..XXX..X..XX.X..XX.
.X.XX.X.X..XX.XX.XX
X.X.XX...X.XXXXX..X
.X..XXX...X..XX.X..
X..X..XX.X.XX.XX.X.
..X....XX..XX..XX.X
.X.....XXXX..X..XX.
X.....X..XX...X..XX
.....X..X.XX...X..X
....X..X...XX...X..
...X..X.....XX...X.
..X..X.......XX...X
.X..X.........XX...
X..X...........XX..
..X.............XX.
.X...............XX
X.................X
...................

I'll post the Mathematica Notebook on MathematicaCookbook.com

Hold Everything!

In Mathematica Cookbook recipe 2.2 I discuss a "technique" for holding arbitrary arguments. Well truth be told, this recipe is pretty lame and was not supposed to make it into the final book. This post is intended to make of for that lameness and hopefully clarify some of the bewildering multitude of ways Mathematica offers for Holding.

The Hold Family of Attributes

Mathematica's default behavior is to evaluate every expression it sees. Here is an example.

In[1]:= a=1;b=2;c=3; d= c; e:= d; f :=e;
{a,b,c, d,e, f}
Out[2]= {1,2,3,3,3,3}

Here we associate symbols a, b, c with integers 1, 2,3. We then associate d with the value of symbol c and we associate e with symbol d telling Mathematica not to evaluate d just yet by using SetDelayed (:=). We also associate f with symbol d likewise delaying evaluation. Later, when we evaluate the list containing these symbols Mathematica keeps evaluating until there is nothing left to do and we get the result as integers. By using Trace we can see the steps Mathematica goes through.

In[3]:= Trace[{a,b,c, d,e, f}]
Out[3]= {{a,1},{b,2},{c,3},{d,3},{e,d,3},{f,e,d,3},{1,2,3,3,3,3}}

Okay, simple enough. However, occasionally you want to write functions that act on expressions before they are evaluated. In fact, even if you never had a reason for doing so, Mathematica itself needs this capability. For example, Mathematica could not implement the function SetDelayed if it did not have a way of saying "don't evaluate". Rather then creating certain functions with special no-evaluating behavior, Mathematica takes a general approach via the concept of attributes. You can inspect the attributes of a symbol using the command Attributes.

In[4]:= Attributes[SetDelayed]
Out[4]= {HoldAll,Protected,SequenceHold}

Here we see SetDelayed has two attributes in the hold-family: HoldAll and SequenceHold. Let's explore these using our own symbols and SetAttributes. Here I use symbols f1, f2 and so on without associating them with values because that is not necessary to illustrate the behavior of the attributes.

In[5]:= (* Make sure these symbols have no values or attributes*)
ClearAll[f1,f2,f3,f4,f5,f6,f7] ;

In[6]:= (* Here you can see that our list of symbols evalauates as before when we wrap f1 around it*)f1[{a,b,c,d,e,f}]
Out[6]= f1[{1,2,3,3,3,3}]

In[7]:= (*Here we use the HoldAll and associate it with f2. This says that no arguement of f2 should be evaluated*)
SetAttributes[f2, HoldAll]
In[8]:= f2[{a,b,c,d,e,f}]
Out[8]= f2[{a,b,c,d,e,f}]

In[9]:= (* It does not matter how many seperate arguments are passed, they are all held*)f2[a,b,c]
Out[9]= f2[a,b,c]

Notice the difference between evaluating f1 which has no attributes and f2 which has attribute HoldAll. In essence, f2 acts the same as the built-in Mathematica command Hold.

In[10]:= Attributes[Hold]
Out[10]= {HoldAll,Protected}

Sometimes you want only the first argument to a function to be held unevaluated. For that you use attribute HoldFirst.

In[11]:= SetAttributes[f3, HoldFirst]
f3[a,b,c]
Out[12]= f3[a,2,3]

Alternatively you may want all arguments but the first to be held. Here you use HoldRest.

In[13]:= SetAttributes[f4, HoldRest]
f4[a,b,c]
Out[14]= f4[1,b,c]

So far I think the mechanics of HoldAll, HoldFirst and HoldRest should be pretty clear. Don't worry yet if you don't get why you would want to use these in your own code, I'll get to that later. Just make sure you are comfortable with the idea of using attributes to suppress Mathematica's desire to evaluate before reading on.


Okay, now I want to point out some important exceptions. Well, not really exceptions but rather clarifications. HoldAll, HoldFirst and HoldRest do not suppress every action that Mathematica takes when it sees an expression. To illustrate this, please recall that Mathematica has a long hand way of specifying a sequence of things.

In[15]:= Sequence[1,2,3]
Out[15]= Sequence[1,2,3]

A sequence is different form a list in that Mathematica will magically flatten out sequences and splice the result into any function call.

In[16]:= f1[Sequence[1,2,3]]
Out[16]= f1[1,2,3]
In[17]:= f1[Sequence[1,2,Sequence[3,4,5]],7,8,9]
Out[17]= f1[1,2,3,4,5,7,8,9]

Now, recall that f2 has attributes HoldAll. What do you think Mathematica does if we give f2 a sequence?

In[18]:= f2[Sequence[1,2,Sequence[3,4,5]],7,8,9]
Out[18]= f2[1,2,3,4,5,7,8,9]

Ah! The automatic flatting is not suppressed by HoldAll nor is it suppressed by HoldFirst or HoldRest. The flattening out of sequences is something you usually don't want to suppress. This is why it is an exception to the rule for the standard hold family of attributes. However, you may recall that SetDelayed had another hold-like attribute called SequenceHold and you can probably guess what it does!

In[19]:= SetAttributes[f5, SequenceHold]
f5[Sequence[1,2,Sequence[3,4,5]],7,8,9]
Out[20]= f5[Sequence[1,2,3,4,5],7,8,9]

Notice how the outer sequence remains. Why did the inner sequence get flattened? Simply because it is evaluated within Sequence which naturally does not have the SequenceHold attribute.


Okay, lets review. HoldAll, HoldFirst and HoldRest are attributes that suppress evaluation of specific arguments but don't suppress sequence flattening. For that you use SequenceHold. You can use one of the hold attributes together with SequenceHold to get both behaviors.

In[21]:= SetAttributes[f6, {HoldAll,SequenceHold}]
f5[Sequence[a,b,c]]
f6[Sequence[a,b,c]]
Out[22]= f5[Sequence[1,2,3]]
Out[23]= f6[Sequence[a,b,c]]

But we are not done yet! There is even a stronger form of holding that suppresses normal evaluation, sequence flattening and more! First, recall that sometimes you want to tell Mathematica that you want something to evaluate despite the presence of HoldAll and friends. To see how there can be a stronger form of holding we must first consider Evaluate.

In[24]:= f2[Evaluate[a,b,c]]
Out[24]= f2[1,2,3]

Evaluate is a way of saying to Mathematica that you know what you are doing and you want evaluation to take place despite the presence of HoldAll, etc.. This typically arises when you want to plot a function that you obtain by integration (or other function generating operations).

In[25]:= Attributes[Plot]
Out[25]= {HoldAll,Protected}

In[26]:= Plot[Evaluate[Integrate[Sin[x],x]], {x, 0, 2Pi}]
Out[26]:=

The attribute HoldAllComplete has super-powers because it can even shield evaluation by Evaluate!

In[27]:= SetAttributes[f7,HoldAllComplete]
In[28]:= f6[Evaluate[{a,b,c}]]
f7[Evaluate[{a,b,c}]]
Out[28]= f6[{1,2,3}]
Out[29]= f7[Evaluate[{a,b,c}]]

The built in command HoldComplete is the counterpart to the Hold command.

In[30]:= Attributes[Hold]
Attributes[HoldComplete]
Out[30]= {HoldAll,Protected}
Out[31]= {HoldAllComplete,Protected}
In[32]:= HoldComplete[Evaluate[{a,b,c}]]
Out[32]= HoldComplete[Evaluate[{a,b,c}]]

Functions with attribute HoldAllComplete also suppress upvalue evaluation. I am not going to discuss up upvalues here but you can refer to the Mathematica documentation or my cookbook.

HoldForm, Unevaluated, Defer oh my!

At this point you might think we have exhausted all the ways you can suppress evaluation but no. Mathematica has some more subtle ways you can keep its evaluation engine in check. Perhaps the easiest to understand is HoldForm. This command suppresses evaluation just like Hold except the command gets hidden when display in output form. The following should make the difference clear.

In[33]:= Hold[1+2]
Out[33]= Hold[1+2]

In[34]:= ReleaseHold[%]
Out[34]= 3

In[35]:= HoldForm[1+2]
Out[35]= 1+2

In[36]:= ReleaseHold[%]
Out[36]= 3

Unevaluated can be thought of as a temporary or one-shot Hold. It suppresses evaluation the first time the Mathematica evaluator sees it but not subsequent times. The following examples are a good illustration of the difference.

In[37]:= (*first I create a simple list *)
list = {1,2,3}
Out[37]= {1,2,3}

In[38]:= (*Now lets see what happens if we try to make list  list self-referential*)
list[[3]] = list; list
Out[38]= {1,2,{1,2,3}}

In[39]:= (* What happens if we doo this again using Hold? *)
list = {1,2,3};
list[[3]] = Hold[list];
list

Out[41]= {1,2,Hold[list]}

In[42]:= ReleaseHold[%]
Out[42]= {1,2,{1,2,Hold[list]}}

In[43]:= ReleaseHold[%]
Out[43]= {1,2,{1,2,{1,2,Hold[list]}}}

Okay, that is somewhat interesting. Each time we invoke ReleaseHold the list expands unveiling another nested version of itself. What do you think happens if we do this experiment with Unevaluated instead?

In[44]:= list = {1, 2, 3};
list[[3]] = Unevaluated[list];
list
During evaluation of In[44]:= $RecursionLimit::reclim: Recursion depth of 256 exceeded. >>
During evaluation of In[44]:= $RecursionLimit::reclim: Recursion depth of 256 exceeded. >>
Out[46]= {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,{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,{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,{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,{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,{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,{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,{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,{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,{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,{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,{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,{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,{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,{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,{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,Hold[{1,2,list}]}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}

Oops! We created an never ending evaluation that eventually blows up. The only way for Mathematica to finally display the "result" is for it to force a Hold into the output to put the breaks on this runaway evaluation train! The reason for this behavior is that Unevaluated allowed us to suppress evaluation of list up until the point where it made the symbol list the new third element of the List associated with the symbol list but after that step, evaluation is no longer suppressed. So we have a symbol which contains a reference to itself and the evaluator can never stop until recursion limit is reached.


Defer is another function that acts like Hold but will allow evaluation each time it is presented to the front-end for evaluation via the user action of hitting Shift-Enter or selection the expression and using Evaluation In Place. We can use the self-referential example to see this. Here each step in the expansion was created by hitting Shift-Enter on the prior output.

In[47]:= list = {1, 2, 3}; 
list[[3]] = Defer[list]; 
list
Out[49]= {1,2,list}

In[50]:= {1,2,list}
Out[50]= {1,2,{1,2,list}}

In[51]:= {1,2,{1,2,list}}
Out[51]= {1,2,{1,2,{1,2,list}}}

In[52]:= {1,2,{1,2,{1,2,list}}}
Out[52]= {1,2,{1,2,{1,2,{1,2,list}}}}

So you can see Mathematica has quite a rich repertoire of functions and attributes whose sole purpose is to keep it for doing what it was designed to do - evaluate! The novice Mathematica user may be mystified by this but these features are essential to the functionality of Mathematica. In other words, the rich sets of features Mathematica provides would not be there if Mathematica did not contain mechanisms for controlling itself. Programming after all is the act where an individual exerts control over a (abstract) machine. But perhaps this explanation is a bit too metaphysical for your tastes. So lets get to some concrete examples. Consider the rich family of Plot functions which take other functions as input. Notice that these all have HoldAll as attributes. Think for a second why this is the case before reading on.

In[53]:= Attributes[{Plot, ParametricPlot, Plot3D}]
Out[53]= {{HoldAll,Protected},{HoldAll,Protected},{HoldAll,Protected}}

Functions that rely on delayed evaluation are typically those that are in essence little evaluation engines themselves. Think about plotting. It must take a function and evaluate it at various points so that it can map the values at those points to the graphics coordinates of points (ultimately pixels on the display device). These evaluations can not be interfered with by the normal evaluation process because in many cases that will change the input before the evaluator can do its thing. If this reasoning is correct than any function in Mathematica which repeatedly evaluates another function should utilize the Hold family of attributes. One of the simplest of these is Table so let's see.

In[54]:= Attributes[Table]
Out[54]= {HoldAll,Protected}

Yep. Various forms of control flow (If, Do, Switch) also must use held arguments for similar reasons. Another class of functions that hold arguments unevaluated are those that must act on symbolic values themselves. Examples are Clear and AddTo (+=). The following little program lists all such symbols in the System` context that have an attribute in the Hold family.

In[55]:= Select[Names["System`*"],Length[Intersection[Attributes[#],{HoldAll,HoldFirst,HoldRest,HoldAllComplete}]]>0&]
Out[55]= {AbortProtect,AbsoluteTiming,AddTo,And,Animate,AppendTo,Arrow3DBox,ArrowBox,Assuming,Attributes,BezierCurve3DBox,BezierCurveBox,Block,BlockRandom,BSplineCurve3DBox,BSplineCurveBox,BSplineSurface3DBox,Button,CancelButton,Catch,Check,CheckAbort,CheckAll,ChoiceButtons,CircleBox,Clear,ClearAll,ClearAttributes,Compile,CompiledFunction,CompoundExpression,Condition,ConeBox,ConsoleMessage,Context,ContinuedFractionK,ContourPlot,ContourPlot3D,Control,ControlActive,ControllerManipulate,CuboidBox,CylinderBox,Debug,DebugTag,Decrement,DefaultButton,DefaultValues,Defer,Definition,DensityPlot,Dialog,DialogInput,DialogReturn,DiskBox,DivideBy,Do,DownValues,DumpSave,Dynamic,DynamicBox,DynamicModule,DynamicModuleBox,DynamicWrapper,DynamicWrapperBox,Exists,FileName,FindArgMax,FindArgMin,FindMaximum,FindMaxValue,FindMinimum,FindMinValue,FindRoot,For,ForAll,FormatValues,FullDefinition,Function,GeometricTransformation3DBox,GeometricTransformationBox,Graphics3DBox,GraphicsBox,GraphicsComplex3DBox,GraphicsComplexBox,GraphicsGroup3DBox,GraphicsGroupBox,Hold,HoldComplete,HoldForm,HoldPattern,If,Increment,Information,InsetBox,Interpretation,InterpretationBox,Line3DBox,LineBox,LineIntegralConvolutionPlot,Literal,MakeBoxes,Manipulate,MatchLocalNameQ,MemoryConstrained,MenuItem,Message,MessageName,MessagePacket,Messages,Module,Monitor,Nand,NCache,NIntegrate,Nor,NProduct,NSum,NValues,Off,On,Or,OwnValues,ParametricPlot,ParametricPlot3D,Parenthesize,Pattern,PatternTest,Piecewise,Play,Plot,Plot3D,Point3DBox,PointBox,Polygon3DBox,PolygonBox,PreDecrement,PreemptProtect,PreIncrement,PrependTo,Product,Protect,Quiet,RasterBox,Reap,RectangleBox,Refresh,RegionPlot,RegionPlot3D,Remove,RuleCondition,RuleDelayed,SampledSoundFunction,Save,Set,SetAttributes,SetDelayed,SphereBox,Stack,StackBegin,StackComplete,StackInhibit,StreamDensityPlot,StreamPlot,SubtractFrom,SubValues,Sum,Switch,SystemException,Table,TagSet,TagSetDelayed,TagUnset,Text3DBox,TextBox,TimeConstrained,TimesBy,Timing,Trace,TraceDialog,TracePrint,TraceScan,TubeBezierCurveBox,TubeBox,TubeBSplineCurveBox,Unevaluated,Unprotect,Unset,UpSet,UpSetDelayed,UpValues,ValueQ,VectorDensityPlot,VectorPlot,VectorPlot3D,WaitUntil,Which,While,With,$ConditionHold,$Failed}

Here you can see that the ultimate form of holding, HoldAllComplete, is more rarely used and when it is it is for rather low-level functions.

In[56]:= Select[Names["System`*"],Length[Intersection[Attributes[#],{HoldAllComplete}]]>0&]
Out[56]= {DebugTag,HoldComplete,InterpretationBox,MakeBoxes,Parenthesize,PreemptProtect,SystemException,Unevaluated}

Recipe 2.2 Reconsidered

In recipe 2.2 of Mathematica Cookbook I consider if it was possible to create functions that Hold other combinations of arguments than provided by HoldAll, HoldFirst and HoldRest. The solution proposed using Hold as a pattern within the function itself. This awkward construct thus required you to use Hold when you invoked the function. To further un-motivate this I presented a rather lame example in the solution.

In[57]:= array1 = Table[0, {10}]; array2 = Table[1, {10}];

arrayAssign[Hold[a_Symbol],aIndex_,Hold[b_Symbol],bIndex_]:= 
Module[{},
a[[aIndex]] = b[[bIndex]];
a[[aIndex]]]

(*Assign elements 2 through 3 in array 2 to array1 *)
arrayAssign[Hold[array1],2;;3,Hold[array2],1];
array1

Out[60]= {0,1,1,0,0,0,0,0,0,0}

There are several reasons this solution is lame. First off, the example is not at all practical. There is little reason to create a function to do this when you can do the same in a one line expression. But that could be forgiven by virtue of being a purely pedagogical example. The real flaw is that this technique requires you to use Hold at the call site which is nothing at all like the behavior of functions with attributes HoldFirst, HoldRest or HoldAll. A more sensibly way to achieve the same effect is to simply use HoldAll and force evaluation where required. So to use this somewhat useless example again...

In[61]:= array1 = Table[0, {10}]; array2 = Table[1, {10}];

SetAttributes[arrayAssign2,HoldAll];arrayAssign2[a_Symbol,aIndex_,b_Symbol,bIndex_]:= 
Module[{aIndex2,bIndex2},
{aIndex2,bIndex2} = Evaluate[{aIndex,bIndex}];
a[[aIndex2]] = b[[bIndex2]];
a[[aIndex2]]]

(*Assign elements 2 through 3 in array 2 to array1 *)
arrayAssign2[array1,2;;3,array2,1];
array1

Out[64]= {0,1,1,0,0,0,0,0,0,0}

Let's Make a Deal - Let Monty Rest!

Some "solved" problems just never go away. Perpetual Motion machines is one example. Another example is the Monty Hall Problem. There is presently a mini-debate in the discussion area of Scientific American's How Randomness Rules Our World and Why We Cannot See It with a significant number of adamant doubters of the standard result that states it is better to switch if Monty shows you a goat.

This problem is amazingly obvious to understand once you analyze it correctly and remove certain ambiguities from the problem statement. Here's my analysis and some Mathematica simulations to add some weight (as if any is needed).

Okay, we all can agree that the probability of NOT picking the DREAM VACATION is 2/3, right? There are two GOATS and one DREAM VACATION.

Now, when Monty shows you the remaining door with a GOAT he just beamed you some very significant information. He's told you that if you picked a GOAT then the probability of getting a DREAM VACATION is 1 if you switch! We already know the probability you picked a GOAT is 2/3 so after he gives you this new info your probability of winning is now 2/3. So switch for GOAT's sake!! If you don't switch, your probability is just 1/3.

Here is a Mathematica program for the non-believers.

GOAT = 0; (* Goat worth zero *)
VACATION = 1; (* Vacation worth one*)

makePrizes[] := Module[{},Switch[RandomInteger[{1,3}],
1,{GOAT,GOAT,VACATION},
2,{GOAT,VACATION,GOAT},
3,{VACATION,GOAT,GOAT}]]

randomPick[doors_List] := Module[{},RandomInteger[{1,Length[doors]}]]

strategy1VS2[trials_Integer] :=
Module[{winnings1=0, winnings2=0, firstPick, secondPick, doors, doors2},
SeedRandom[];
Do[doors = makePrizes[];
firstPick = randomPick[doors];
(*winnings of person who keeps first pick*)
winnings1+= doors[[firstPick]];
(*delete first pick from choices*)
doors2 = Drop[doors,{firstPick}];
(*delete goat from remaining*)
doors2 = Drop[doors2,Position[doors2,GOAT][[1]]];
(*Always pick remaining prize *)
secondPick =doors2[[1]];
(*winnings of person who switches*)
winnings2+= secondPick,{trials}];
{winnings1,winnings2}]

(*Run simulation 10000 times. *)
strategy1VS2[10000]
{3356,6644}

The result {3356,6644} means keeping first choice only paid 3356 over 10000 runs but switching paid 6644!

Now, there are ASSUMPTIONS here (there always are). One assumption is that on each run the position of the prize changes. It turns out that keeping the prize always in any particular door for the entire simulation does not matter (as long as the contestant does not have the information, obviously!)

strategy1VS2A[trials_Integer,init_List] :=
Module[{winnings1=0,winnings2=0,firstPick,secondPick,doors,doors2},
SeedRandom[];
Do[doors = init;
firstPick = randomPick[doors];
(*winnings of person who keeps
first pick*)
winnings1+= doors[[firstPick]];
(*delete first pick from choices*)
doors2 = Drop[doors,{firstPick}];
(*delete goat from remaining*)
doors2 = Drop[doors2,Position[doors2,GOAT][[1]]];
(*Always pick remaining prize *)
secondPick =doors2[[1]];
(*winnings of person who switches*)
winnings2+= secondPick;,{trials}];
{winnings1,winnings2}]

strategy1VS2A[10000,{GOAT,GOAT,VACATION}]
{3316,6684}

strategy1VS2A[10000,{GOAT,VACATION,GOAT}]
{3267,6733}

strategy1VS2A[10000,{GOAT,GOAT,VACATION}]
{3382,6618}

The other assumption is that your not forced to switch before seeing the goat. This IS important!!

strategy1VS2B[trials_Integer] :=
Module[{goatPositions,pos,winnings1=0,winnings2=0,firstPick,secondPick,doors,doors2},
SeedRandom[];
Do[doors = makePrizes[];
firstPick = randomPick[doors];
(*winnings of person who keeps first pick*)
winnings1+= doors[[firstPick]];
(*delete first pick from choices*)
doors2 = Drop[doors,{firstPick}];
(*Randomly choose from remaing*)
secondPick =randomPick[doors2];
(*winnings of person who switches*)
winnings2+= doors2[[secondPick]];,{trials}];
{winnings1,winnings2}]

strategy1VS2B[10000]
{3373,3295}

So information has value, Duh!

So now that you have this information, become a believer, make the switch!

Tan of the Kitchen Sink

In a post on MathForum I made the mistake of doubting Mathematica could compute the Tan of the Kitchen Sink. It is always a mistake to doubt Mathematica's prowess as Daniel Lichtblau of Wolfram pointed out to me:


In[1]:= Tan[Khinchin//Sinc] // N
Out[1]= 0.165514

I stand corrected!


p.s.

It sort of spoils the joke/pun to explain but to non-Mathematca users...

Khinchin's constant is aprox. 2.68545
Sinc[x] = Sin[x]/x
N means "give numeric value"
and // means "use postix" so
this computes N[Tan[Sin[Khinchin]/Khinchin]]

F# For Scientists Misses the Boat On Mathematica Performance

I recently purchased F# For Scientists by Dr. Jon Harrop after the author mentioned it on the Mathematica Mailing List. According to Dr. Harrop,

Mathematica's .NET-Link technology allows Mathematica and .NET programs to interoperate seamlessly. Moreover, Microsoft's new functional programming language F# provides many familiar benefits to Mathematica programmers:
The marriage of Mathematica with F# can greatly improve productivity for a wide variety of tasks.

I am a big fan of Mathematica and functional programming and have been wanting to check out F# for some time so I decided to give the book a shot. It just arrived today so I can't post a full review but I did jump directly to the small section (5 pages) on using F# with Mathematica.

What did I learn? Well this section rightly claims that Mathematica has awesome symbolic math capabilities (it does). But then it goes on to claim that F# can beat the pants off of Mathematica on raw calculation. Thus it suggested F# programmers should call out to Mathematica for symbolic integration but then evaluate the result in F# for speed (to the tune of 3.4 times Mathematica's speed). I was naturally dubious. The explanation of this speed up is give as

The single most important reason for this speed boost is the specialization of the F# code compared to Mathematica's own general purpose term rewriter. ... Moreover, the F# programming language also excels at compiler writing and the JIT-compilation capabilities of the .NET platform make it ideally suited to the construction of custom evaluators that are compiled down to native code before being executed. This approach is typically orders of magnitude faster than evaluation in a standalone generic term rewriting system like Mathematica.

Okay, hold the phone! First off, I did not know the F# language could write compilers. I'll forgive this as poetic use of language. I guess I sort of know what he meant to say. More interesting is that we have gone from 3.7 times to "orders of magnitude". Now, I don't take anything away from the brilliant folks at Microsoft, but the equally brilliant folks at Wolfram have been focusing exclusively on mathematics software for 20 years and you might think they learned a thing or two about computational speed!

Here is the example from the book...

First, he uses Mathematica to integrate a function.


Integrate[Sqrt[Tan[x]],x]

(-2*ArcTan[1 - Sqrt[2]*Sqrt[Tan[x]]] +
2*ArcTan[1 + Sqrt[2]*Sqrt[Tan[x]]] +
Log[-1 + Sqrt[2]*Sqrt[Tan[x]] - Tan[x]] -
Log[1 + Sqrt[2]*Sqrt[Tan[x]] + Tan[x]])/(2*Sqrt[2])

He then goes to show that Mathematica takes 26 secondsto evaluate this function in loop for 360,000 iterations.

He then shows a translator that converts the Mathematica to F# and the F# code does the same work in 7.595 seconds.

So far Dr. Harrop is correct but like some many others who are in a rush to show their new favorite language superior to another's, he forgets to read the manual! Particularly, the section on optimization! If he had he would have found a handy little Mathematica function called Compile. Hmm, sounds promising. And in fact....


cf = Compile[{{x, _Complex}}, Evaluate[Integrate[Sqrt[Tan[x]],x]]]
Timing[Do[cf[x + y I],{x,-3.0,3.0,0.01},{y,-3.0,3.0,0.01}]]

{5.281,Null}

That's 5.281 seconds on my relatively underpowered laptop (Thinkpad X60) !

Some might feel I'm being a bit harsh on Dr. Harrop but after all he made me layout bucks for a book that promised me "many familiar benefits" only to deliver 5 measly pages of half truth. F# programmers may benefit from Mathematica but the jury is still out as to whether the reverse is true.

A Biological Programming Language

I've often said to friends that if I could start my career over again I would go into biology instead of computer science. Now, perhaps, there is a way to have a foot in both worlds.

Little b is a programming language for modeling biological systems. Quoting from the languages site...

The little b project is an effort to provide an open source language which
allows scientists to build mathematical models of complex systems. The
initial focus is systems biology. The goal is to stimulate widespread sharing
and reuse of models. The little b language to allow biologists to build
models quickly and easily from shared parts, and to allow theorists to program
new ways of describing complex systems.

Currently, libraries have been developed for building ODE models of
molecular networks in multi-compartment systems such as cellular epithelia.
Aneil Mallavarapu is the author and inventor of little b, and runs the
project. Little b is based in Common Lisp and contains mechanisms for rule-based
reasoning, symbolic mathematics and object-oriented definitions. The syntax is
designed to be terse and human-readable to facilitate communication. The
environment is both interactive and compilable.

Makes me wonder if Mathematica would be a good enviornment for similar exploration but with more sophisticated tools already built in.

Semantic Vectors Revisited (for 31 bucks!)

Its been a while since I posted (been busy writing my new book) and even longer since I've posted any ideas related to my ideas on using concepts from linear algebra to model intelligence. But I thought I'd share an experience that makes me wonder how anything got done before the WWW!

While doing research on tensors for my book I came across a book called MathTensor: A System for Doing Tensor Analysis by Computer. This book describes software for Tensor Math developed using Mathematica so it instantly caught my interest. This lead me to one of the author's web sites which lead me to an article Tensor Analysis of Matrix Cognition during Medical Decision-Making. Now you can't put the words matrix and cognition next to each other without getting my immediate attention so I jumped to that essay which ultimately lead me to this gem: A scaling method for priorities in hierarchical structures by
Thomas L. Sattay written in the Journal Mathematical Psychology 1977; 15:234-281 (there is no online version but you can buy a PDF copy at ScienceDirect if you are willing to part with $31.

I found this research to be fascinating and gave me much food for thought that I'll try to share when I have more time. For now I'd only like to make the following rather obvious observation. If it was not for the web, there would be close to zero chance that I would have found this article and an even smaller chance I would be reading it within 15 mins of finding it. The only sad part is that it is locked up in some obscure journal that I did not have immediate access to without parting with the cost of a nice dinner. I think publishers of journals need to catch up with the rest of the world and begin opening up their older content to free access. Clearly they can use advertising to subsidize this but perhaps advertisement driven business models have reached a point of saturation. Perhaps its time for a library based approach to become virtualized.

I am sure I could find a library within a reasonable vicinity of my home that had access to this journal but who has the time! Why not offer a version that rather than costing $31 to keep forever, costs me $1 to read for a day and $0.50 for each additional day. DRM technology is certainly good enough to make this work. And I am guessing that the publishers would make more money than by waiting for someone like me who was motivated enough to part with $31. There is a vast amount of lost knowledge hiding in these journals. History has shown that the world benefits greatly when such knowledge is serendipitously rediscovered (think Gregor Mendel and his Bean Plants). Its time to unlock the vaults of knowledge so creativity and discovery can reach new unimagined heights!

Red-Black Tree in 2 hours

Many computer science students have heard of red-black trees. If you use any container class library that has a map (ordered associative container) like STL you probably know that red-black trees are a popular implementation of a map. It may be non-obvious why the constraints associated with red-black trees cause them to remain roughly balanced but what is even more daunting is the implementation of such trees in an imperative language like C!

While working on an implementation for a recipe in my forth coming "Mathematica Cookbook" I found a functional implementation in Haskell (postscript) by Chris Okasaki. I know this has been said a million times before but it never ceases to amaze me how succinct and beautiful the functional approach can be. Using the Haskell solution as a guide, I was able to develop a complete red-black implementation in Mathematica in under 2 hours. This may not sound that impressive but consider that the referenced paper does not show how to implemented a remove operation. Even without the need for a remove, how many programmers could take a red-black tree written in say C and translate it into a completely working implementation in Java in under 2 hrs? I suppose a few can but this is really not a post about bragging rights. The functional approach to software development is just god-damn beautiful at so many levels and it is this and not my hacker abilities which made this exercise possible.

Now, there are caveats. There always are. Any C implementation of a red-black tree is bound to have numerous optimizations and a purely functional solution will not fare well for every application of a map (but surprisingly it is competitive for many. I'll post some C++ comparisons when I have a chance.)

If you'd like to learn more about red-black trees but don't know Haskell I would highly recommend learning the minimum of Haskell you need to understand Okasaki's paper instead of trying to learn about them by digging into a C implementation first.

How do you know when you have mastered a new programming language

This is a minor continuation of my thoughts from The Law of the Excluded Middle does not apply to Programming Languages.

One sure sign that you have mastered a language is when you can create a fairly comprehensive list about what sucks about the language (while simultaneously appreciating why some of this suckiness is a necessary evil).

When you first meet a new language that floats your boat, there is a tendency is to fall in love. This happened to me not too long ago with Erlang. You think, "This language is great. Such and such is so hard to do in Language X and look how easy it is in Language Y.

Well, much like in the real world of human relationships you don't really know what love means until you get married! When you get married to a language (commit to developing a multi-year non-trivial system in it) then your love is surely tested. You learn about the languages warts and its tendency to leave its socks outside the hamper, squeeze the tooth paste from the top and leave the toilet seat in an inconvenient position!

If you still love the language with its warts and all, then you have some of the necessary (but not necessarily sufficient) hallmarks of a master. Either that, or you ave a really good therapist.

p.s. Erlang and I broke up but we are still friends! On a happier note, my long term mistress (Mathematica) and I are really making sparks fly! Hope C++ doesn't catch us.

Interval Math

While doing research for the Numerics Chapter of my forthcoming Mathematica Cookbook I came across a site devoted to research on Interval Math. Interval Math is an approach from the domain of Numerical Analysis that deals with the fact that all measurements are imprecise by abandoning the representation of measured values by numbers. Instead of numbers, it defines all mathematical operations on intervals.

Mathematica (as of version 5) support real (but not complex) interval math where intervals take the form Interval[{min1,max1}...]. All of the typical mathematical operations and functions are defined for intervals.

Interval math is important for computer systems that must act intelligently in the real world. All sensors are approximate. This is true for man-made devices as well as for our own eyes and ears. If a sensor on a robot returns a particular value there is always an inherent error. Rather than deal with errors by sampling and averaging, interval math allows the error to directly be represented in the values that enter downstream computations. This means all intermediate results track the propagation of errors from multiple sources to yield better information. There also seems to be a relationship between interval computation and fuzzy sets but it I have not located any resources except on paid content sites.

It seems that although the study of Interval math began in the US it is largely forgotten while in Germany it is there are conferences and it is part of the qualifying exams for studies in numerical methods.

Some of the less technical resources on the earlier mentioned site are this introduction, an article from American Scientist and even a movie.

Mathematica on LinkedIn and on a Wiki

This is my third post today (penance for not posting for so long!)

I recently started a LinkedIn Group called Mathematica Users Group. If you are a member of LinkedIn you can join the group by clicking here. After creating the group I thought it would be cool to have a Mathematica Wiki and soon discovered that Luc Barthelet thought this was a good idea too but thought of it a few years earlier than I did!

Prof. Ray C. Dougherty's Research

A few weeks ago I attended a Wolfram Research event called Mathematica Publishers Day. The goal of the event was to highlight the capabilities of Mathematica 6 as a platform for technical publishing. I really enjoyed this event but was also pleasantly surprised by a talk that did not quite fit into the overall theme of the event but nevertheless was quite fascinating to me.

The presenter was Prof. Ray C. Dougherty, NYU Linguistics researcher. He used Mathematica to model all possible sine wave based communications systems. The presentation is available via Wolfram. Unfortunately, as with many interesting presentations, you needed to hear the talk to get the most out of it. Here are some interesting excerpts that I remember:

• The Cochlea is computing the second derivative of the auditory input.
• The most mathematically complex communication system is one where the transmitter and receiver have the same anatomy (e.g., wings of insects).
• Bats can hear phase changes because they can rotate their ears. A human can not hear a change in the rotation of a tuning fork but a bat can.
• Prof Dougherty believes he has a Chomsky generative grammar that enumerates all possible animal communication systems.
• He also believes he can map each possible system onto the integers in a natural way.
• From this he concludes that evolution must proceed in jumps.
• He relates this idea to the evolution of all possible Tic Tac Toe Games to illustrate the notion that all such games are not unique and similarly the space of all possible communication systems contains many redundant systems as well.
• He goes on to visualising distributions of the primes to illustrate that there are systems that are not random but whose patterns are too complex for us to model in a simple fashion. Explains how this is related to the ideas in Stephen Wolfram's NKS.

Stephen Wolfram on Software Design and Naming

In this post Stephen Wolfram talks about design reviews and the importance of naming things correctly. I sought of envy him for working in a company where one can afford to spend a total of 10,000 hours in design review and where you sweat over the details of naming every last function. As much as I can see the merits of Agile Development, I have never worked on an Agile project that creating anything nearly as elegant as Mathematica. Sure, Agile focuses on working software and "good enough" software and the economics of most software projects makes this mode of development necessary. But it would be nice to one day work on a project where time to market was less important than the end result.