About alternatives to Matlab

[Carl Banks]

Jon Harrop wrote:
> I don't know Python but this benchmark caught my eye.
>
> >>def D4_Transform(x, s1=None, d1=None, d2=None):
> >>   """
> >>   D4 Wavelet transform in NumPy
> >>   (C) Sturla Molden
> >>   """
> >>   C1 = 1.7320508075688772
> >>   C2 = 0.4330127018922193
> >>   C3 = -0.066987298107780702
> >>   C4 = 0.51763809020504137
> >>   C5 = 1.9318516525781364
> >>   if d1 == None:
> >>      d1 = numpy.zeros(x.size/2)
> >>      s1 = numpy.zeros(x.size/2)
> >>      d2 = numpy.zeros(x.size/2)
>
> Are these definitions ever used? It looks like s1, d1 and d2 are all
> redefined below without reference to these previous values.

No, they're never redefined (even in the recursive version).  Slices of
them are reassigned.  That's a big difference.

> >>   odd = x[1::2]
> >>   even = x[:-1:2]
> >>   d1[:] = odd[:] - C1*even[:]
> >>   s1[0] = even[0] + C2*d1[0] + C3*d1[-1]
> >>   s1[1:] = even[1:] + C2*d1[1:] + C3*d1[:-1]
> >>   d2[0] = d1[0] + s1[-1]
> >>   d2[1:] = d1[1:] + s1[:-1]
> >>   even[:] = C4 * s1[:]
> >>   odd[:] = C5 * d2[:]
>
> Does that line create an array that is never used? If so, is C5 also never
> used?

No it doesn't create a new array.  But it does have an effect.
Explained below.

(Actually, it does create a temporary array to store the intermediate
result, but that array does not get bound to odd.)


> >>   if x.size > 2:
> >>
> >>D4_Transform(even,s1[0:even.size/2],d1[0:even.size/2],d2[0:even.size/2])
>
> What is the result of this function?

It modifies the array even in place.  This also has an effect.


> I'm interested in translating this function into other languages, like
> OCaml, to see how good Python's performance is (I am amazed it can beat
> Matlab, not that I've used Matlab).

Matlab has a few *cough* limitations when it comes to hand-optimizing.
When writing naive code, Matlab often is faster than Python with numpy
because it has many commerical man-year of optimizing behind it.
However, Matlab helps v


> In particular, I think you are eagerly
> allocating arrays when, in a functional language, you could just as easily
> compose closures.

You are completely wrong.

What you seem to be missing is that slice assignment shares data.  For
example, the statement "odd = x[1::2]" does NOT create a new array.  It
takes a slice of x's data (in this case, the odd elements of x: slices
need not be contiguous), and binds it to the symbol odd.  The array
elements are shared; changing elements in odd also changes elements of
x.

So later on, the statement "odd[:] = C5*d2[:]" reassigns the elements
of odd, which, because the array elements are shared, also reassigns
the odd elements in x.

(Notice the difference between "odd[:]=" and "odd=".  "odd[:]=" is
slice assignment; it changes the elements in the given slice of odd, in
this case the whole array.  "odd=" rebinds the symbol odd to a
different object, which could be a new array, a slice of an existing
array, or any other Python object.)


This data sharing more or less accomplishes the same thing that the
"closures" you speak of accomplish (in this case), only without the
functional.


> For example, this program is 3x faster than the Python on my machine:
>
> let rec d4 s1 d1 d2 x =
>   let c1 = 1.7320508075688772 in
>   let c2 = 0.4330127018922193 in
>   let c3 = -0.066987298107780702 in
>   let c4 = 0.51763809020504137 in
>   let c5 = 1.9318516525781364 in
>   let n = Array.length x in
>   let odd i = x.(2*i) and even i = x.(2*i + 1) in
>   let d1 i = odd i -. c1 *. even i in
>   let f = function -1 -> n/2 - 1 | i -> i in
>   let s1 i = even i +. c2 *. d1 i +. c3 *. d1 (f(i-1)) in
>   let d2 i = d1 i +. s1 (f(i-1)) in
>   let even = Array.init (n/2) (fun i -> c4 *. s1 i) in
>   if n > 2 then d4 s1 d1 d2 even else s1, d1, d2
>
> but I'm not sure it is correct!

It's not correct, but what you left out is probably low cost.

OCaml is compiled to machine code, right?  And types can be inferred at
compile time, correct?  Well then of course it's faster.  It seems to
me a big help is the ability to fold multiple array operations into a
single loop, which is optimization a dynamically-typed language like
Python can't easily make.  (It'd require are really smart JIT compiler
or some concessions in dynamicity.)


Carl Banks