MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Mathematica 7 is now available

  • To: mathgroup at smc.vnet.net
  • Subject: [mg94326] Re: Mathematica 7 is now available
  • From: Jon Harrop <jon at ffconsultancy.com>
  • Date: Wed, 10 Dec 2008 04:49:31 -0500 (EST)
  • References: <gg62p3$56g$1@smc.vnet.net> <gh8hkr$r1d$1@smc.vnet.net> <ghavsp$ock$1@smc.vnet.net> <200812081121.GAA15802@smc.vnet.net> <ghlmbp$jof$1@smc.vnet.net>

Daniel Lichtblau wrote:
> Jon Harrop wrote:
>> This is not specific to OCaml. Here are equivalent implementations of
>> that program in C++, Java, OCaml, SML, Scheme, Lisp and Haskell and all
>> well over 100,000x faster than the Mathematica I gave (and many are more
>> concise as well!):
>> 
>>   http://www.ffconsultancy.com/languages/ray_tracer/
> 
> I do not understand this claim that the speed difference is 5 orders of
> magnitude. Some experiments indicate 3, perhaps.

The code I gave is over 5 orders of magnitude slower by my measurements.

> Some modest improvements to the original Mathematica code, though, bring
> it to maybe 2 orders of magnitude of the OCaml speed.

That is a huge improvement. Thank you.

> In the variant 
> below, I make sure everyone is a machine real number or machine integer,
> whichever is appropriate. I use Compile on RaySphere. There are a few
> other tweaks for general efficiency (some might apply to the OCaml
> implementation, I'm not sure).

I have a few questions about this. :-)

> delta = Sqrt[$MachineEpsilon];
> inf = 100000.;

This is nasty. Why is infinity not represented internally by the floating
point number infinity in Mathematica?

> neglight = Normalize[{1., 3., -2.}];
> nohit = {inf, {0., 0., 0.}};
> 
> Clear[RaySphere, Intersect, Create, RayTrace, Main];
> 
> RaySphere = With[{inf = inf}, Compile[

What is the purpose of "With[{inf = inf}, ..]" here?

>      {{c, _Real, 1}, {o, _Real, 1}, {d, _Real, 1}, {r, _Real}},
>      Module[{v = c - o, b, disc, t},
>       b = v.d;
>       disc = b^2 - v.v + r;
>       If[disc < 0., inf,
>        disc = Sqrt[disc];
>        If[b > disc, b - disc,
>         t = b + disc;
>         If[t < 0., inf, t]]]]]];
> 
> Intersect[o_, d_, ln : {lambda_, n_}, List[c_, r_]] := Block[
>    {lambda2 = RaySphere[c, o, d, r]},
>    If[lambda2 >= lambda, ln, {lambda2, c}]]

Defer calculation of the normal. Fair enough but it is only done once per
pixel.

> Intersect[o_, d_, ln : {lambda_, n_}, List[c_, r_, s_]] := Block[
>    {lambda2 = RaySphere[c, o, d, r]},
>    If[lambda2 >= lambda, ln,
>     Fold[Intersect[o, d, #1, #2] &, ln, s]
>     ]
>    ]

No currying and a named subpattern, ok.

> RayTrace[o_, d_, scene_] := Block[
>    {lambda, n, g, p},
>    {lambda, n} = Intersect[o, d, nohit, scene];
>    If[lambda >= inf, 0.,
>     n = Normalize[o + d*lambda - n];
>     g = n.neglight;
>     If[g <= 0., 0.,
>      {lambda, n} =
>       Intersect[o + lambda d + delta n, neglight, nohit, scene];
>      If[lambda < inf, 0., g]]]]
> 
> Create[level_, c_, r_] :=
>   Block[{obj = List[c, r^2]},
>    If[level == 1, obj,
>     Block[{a = 3*r/Sqrt[12], Aux},
>      Aux[x1_, z1_] := Create[level - 1, c + {x1, a, z1}, 0.5 r];
>      List[c,
>       9*r^2, {obj, Aux[-a, -a], Aux[a, -a], Aux[-a, a], Aux[a, a]}]]]]
> 
> prepareTable = Compile[{{n, _Integer}, {s, _Integer}},
>     With[{nss = N[s]},
>      Flatten[
>       Table[Normalize[{(x + s1/nss)/n - .5, (y + s2/nss)/n - .5,
>          1.}], {s1, 0, s - 1}, {s2, 0, s - 1}, {y, 0, n - 1}, {x, 0,
>         n - 1}], 1]
>      ]];

So you precompute a matrix of all the ray directions using a compiled
function.

> Main[level_, n_, ss_] :=
>   Block[{scene = Create[level, {0., -1., 4.}, 1.], nss = N[ss]},
>    scene = MapAll[Developer`ToPackedArray, scene];
>    Total[Map[RayTrace[{0., 0., 0.}, #, scene] &,
>       prepareTable[n, ss], {3}], 1]/nss^2]
> 
> In[22]:= Timing[ff512 = Main[9, 512, 4];]
> Out[22]= {1542.27, Null}

That is orders of magnitude faster than my version but I am not sure why. I
suspect the main performance improvement comes from the use of packed
arrays but I haven't been able to reproduce it in my own code yet.

> At this point, it might make sense to use ParallelMap instead of Map, in
> Main.

Replacing Map with ParallelMap slows it down, even with
Method -> "CoarsestGrained". How can I get a speedup from my 8 cores?! :-)

Many thanks,
-- 
Dr Jon D Harrop, Flying Frog Consultancy Ltd.
http://www.ffconsultancy.com/?u


  • Prev by Date: Re: Re: A 3D Plot Query
  • Next by Date: Re: Drop elements from list
  • Previous by thread: Re: Re: Mathematica 7 is now available
  • Next by thread: Re: Mathematica 7 is now available