Re: About alternatives to Matlab

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

.

Relevant Pages

  • Re: traumatized by pointer casting
    ... The array might not be word-aligned. ... It might have an odd address. ... Instead it will generate a hardware trap. ... The implementation determines the boundaries. ...
    (comp.lang.c)
  • Re: Triangular Matrix
    ... in which zero "0" elements do not have representation. ... R odd and C> N-R: if either test fails, ... How many non-zeroes are in the rows before row R? ... How big must the linear array be? ...
    (comp.lang.c)
  • Re: Triangular Matrix
    ... in which zero "0" elements do not have representation. ... R odd and C> N-R: if either test fails, ... How many non-zeroes are in the rows before row R? ... How big must the linear array be? ...
    (comp.lang.c)
  • Re: Way for computing random primes in standard C.
    ... There are a number of algorithms for finding primes, but, ... If you don't want an array, you could generate random even number from 0 to ... Make them all odd by adding one. ... the number through a naive prime number proof algorithm to begin with, ...
    (comp.lang.c)
  • Re: mXCreateNumericArray (Crashing program)
    ... > from matlab it crashes with "Out of memory. ... > just a 2D array and it still gives the same error.. ... > Below is my fortran file which i compile from matlab with mex cmsl4Interated. ... That way when the mex file returns MATLAB knows that there is nothing in those arrays. ...
    (comp.soft-sys.matlab)