About alternatives to Matlab

[Jon Harrop]

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.

>>   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?

>>   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?

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). In particular, I think you are eagerly
allocating arrays when, in a functional language, you could just as easily
compose closures.

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!

-- 
Dr Jon D Harrop, Flying Frog Consultancy
Objective CAML for Scientists
http://www.ffconsultancy.com/products/ocaml_for_scientists