Mathematica 7 compares to other languages

[jason-sage at creativetrax.com]

Xah Lee wrote:
> 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.
>
>   

For the interested, with MMA 6, on a Pentium 4 3.8Ghz:

The code that Jon posted:

Timing[Export["image-jon.pgm", Graphics at Raster@Main[2, 100, 4]]]
{80.565, "image-jon.pgm"}


The code that Xah posted:

Timing[Export["image-xah.pgm", Graphics at Raster@Main[2, 100, 4.]]]
{42.3186, "image-xah.pgm"}

So Xah's code is about twice as fast as Jon's code, on my computer.

The resulting files were identical (and both looked like pure white 
images; I thought they'd be interesting!).

-J