Mathematica 7 compares to other languages

[Xah Lee]

alright, here's my improved code, pasted near the bottom.

let me say a few things about Jon's code.

If we rate that piece of mathematica code on the level of: Beginner
Mathematica programer, Intermediate, Advanced, where Beginner is
someone who just learned tried to program Mathematica no more than 6
months, then that piece of code is Beginner level.

Here's some basic analysis and explanation.

The program has these main functions:

• RaySphere
• Intersect
• RayTrace
• Create
• Main

The Main calls Create then feed it to RayTrace.
Create calls itself recursively, and basically returns a long list of
a repeating element, each of the element differ in their parameter.

RayTrace calls Intersect 2 times. Intersect has 2 forms, one of them
calls itself recursively. Both forms calls RaySphere once.

So, the core loop is with the Intersect function and RaySphere. Some
99.99% of time are spent there.

------------------

I didn't realize until after a hour, that if Jon simply give numerical
arguments to Main and Create, the result timing by a factor of 0.3 of
original. What a incredible sloppiness! and he intended this to show
Mathematica speed with this code?

The Main[] function calls Create. The create has 3 parameters: level,
c, and r. The level is a integer for the recursive level of
raytracing . The c is a vector for sphere center i presume. The r is
radius of the sphere. His input has c and r as integers, and this in
Mathematica means computation with exact arithmetics (and automatic
kicks into infinite precision if necessary). Changing c and r to float
immediately reduced the timing to 0.3 of original.

------------------
now, back to the core loop.

The RaySphere function contain codes that does symbolic computation by
calling Im, which is the imaginary part of a complex number!! and if
so, it returns the symbol Infinity! The possible result of Infinity is
significant because it is used in Intersect to do a numerical
comparison in a If statement. So, here in these deep loops,
Mathematica's symbolic computation is used for numerical purposes!

So, first optimization at the superficial code form level is to get
rid of this symbolic computation.

Instead of checking whethere his “disc = Sqrt[b^2 - v.v + r^2]” has
imaginary part, one simply check whether the argument to sqrt is
negative.

after getting rid of the symbolic computation, i made the RaySphere
function to be a Compiled function.

I stopped my optimization at this step.

The above are some _fundamental_ things any dummy who claims to code
Mathematica for speed should know. Jon has written a time series
Mathematica package that he's selling commercially. So, either he got
very sloppy with this Mathematica code, or he intentionally made it
look bad, or that his Mathematica skill is truely beginner level. Yet
he dares to talk bullshit in this thread.

Besides the above basic things, there are several aspects that his
code can improve in speed. For example, he used pattern matching to do
core loops.
e.g. Intersect[o_, d_][{lambda_, n_}, Bound[c_, r_, s_]]

any Mathematica expert knows that this is something you don't want to
do if it is used in a core loop. Instead of pattern matching, one can
change the form to Function and it'll speed up.

Also, he used “Block”, which is designed for local variables and the
scope is dynamic scope. However the local vars used in this are local
constants. A proper code would use “With” instead. (in lisp, this is
various let, let*. Lispers here can imagine how lousy the code is
now.)

Here's a improved code. The timing of this code is about 0.2 of the
original. Also, optimization is purely based on code doodling. That
is, i do not know what his code is doing, i do not have experience in
writing a ray tracer. All i did is eyeballing his code flow, and
improved the form.

norm=Function[#/Sqrt@(Plus@@(#^2))];
delta=Sqrt[$MachineEpsilon];
myInfinity=10000.;

Clear[RaySphere];
RaySphere = Compile[{o1, o2, o3, d1, d2, d3, c1, c2, c3, r},
    Block[{v = {c1 - o1, c2 - o2, c3 - o3},
      b = d1*(c1 - o1) + d2*(c2 - o2) + d3*(c3 - o3),
      discriminant = -(c1 - o1)^2 - (c2 - o2)^2 +
        (d1*(c1 - o1) + d2*(c2 - o2) + d3*(c3 - o3))^2 -
        (c3 - o3)^2 + r^2, disc, t1, t2},
     If[discriminant < 0., myInfinity,
      disc = Sqrt[discriminant]; If[(t1 = b - disc) > 0.,
        t1, If[(t2 = b + disc) <= 0., myInfinity, t2]]]]];

Remove[Intersect];
Intersect[{o1_,o2_,o3_},{d1_,d2_,d3_}][{lambda_,n_},Sphere
[{c1_,c2_,c3_},r_]]:=
  Block[{lambda2=RaySphere[o1,o2,o3,d1,d2,d3,c1,c2,c3,r]},
    If[lambda2≥lambda,{lambda,n},{lambda2,
        norm[{o1,o2,o3}+lambda2 *{d1,d2,d3}-{c1,c2,c3}]}]]

Intersect[{o1_,o2_,o3_},{d1_,d2_,d3_}][{lambda_,n_},
    Bound[{c1_,c2_,c3_},r_,s_]]:=
  Block[{lambda2=RaySphere[o1,o2,o3,d1,d2,d3,c1,c2,c3,r]},
    If[lambda2≥lambda,{lambda,n},
      Fold[Intersect[{o1,o2,o3},{d1,d2,d3}],{lambda,n},s]]]

Clear[neglight,nohit]
neglight=N at norm[{1,3,-2}];
nohit={myInfinity,{0.,0.,0.}};

Clear[RayTrace];
RayTrace[o_,d_,scene_]:=
  Block[{lambda,n,g,p},{lambda,n}=Intersect[o,d][nohit,scene];
    If[lambda\[Equal]myInfinity,0,g=n.neglight;
      If[g≤0,
        0,{lambda,n}=Intersect[o+lambda d+delta n,neglight]
[nohit,scene];
        If[lambda<myInfinity,0,g]]]]

Clear[Create];
Create[level_,c_,r_]:=
  Block[{obj=Sphere[c,r]},
    If[level\[Equal]1,obj,
      Block[{a=3*r/Sqrt[12],Aux},
        Aux[x1_,z1_]:=Create[level-1,c+{x1,a,z1},0.5 r];
        Bound[c,3 r,{obj,Aux[-a,-a],Aux[a,-a],Aux[-a,a],Aux[a,a]}]
        ]
      ]]

Main[level_,n_,ss_]:=
  With[{scene=Create[level,{0.,-1.,4.},1.]},
    Table[Sum[
          RayTrace[{0,0,0},
            N at norm[{(x+s/ss/ss)/n-1/2,(y+Mod[s,ss]/ss)/
n-1/2,1}],scene],{s,0,
            ss^2-1}]/ss^2,{y,0,n-1},{x,0,n-1}]]

Timing[Export["image.pgm",Graphics at Raster@Main[2,100,4.]]]


Note to those who have Mathematica.
Mathematica 6 has Normalize, but that's not in Mathematica 4, so i
cooked up above.
Also, Mathematica 6 has AbsoluteTiming, which is intended to be
equivalent if you use stop watch to measure timing. Mathematica 4 has
only Timing, which measures CPU time. My speed improvement is based on
Timing. But the same factor will shown when using Mathematica 6 too.

I'm pretty sure further speed up by 0.5 factor of above's timing is
possible. Within 2 more hours of coding.

Jon wrote:
«The Mathematica code is 700,000x slower so a 50% improvement will be
uninteresting. Can you make my Mathematica code five orders of
magnitude faster or not?»

If anyone pay me $300, i can try to make it whatever the level of F#
or OCaml's speed is as cited in Jon's website. (
http://www.ffconsultancy.com/languages/ray_tracer/index.html ).

Please write out or write to me the detail exactly what speed is
required in some precise terms. If i agree to do it, spec satisfaction
is guaranteed or your money back.

PS Thanks Thomas M Hermann. It was fun.

  Xah
∑ http://xahlee.org/

☄