Re: About alternatives to Matlab

Carl Banks wrote:
fill a (map3 (fun b c d -> b + c + d) b c d)

which will be much faster because it doesn't generate an intermediate
array.

Ah, but, this wasn't about temporaries when you spoke of "eagerly
allocating arrays", was it?

I had thought that all of the array operations were allocating new arrays at
first but it seems that at least assignment to a slice does not.

Does:

a[:] = b[:] + c[:]

allocate a temporary for b[:] + c[:]?

But yes, you're right, in this example
temporary arrays are created; arrays that Ocaml would not need. (I
don't exactly understand why you'd need a functional language to get
this optimization, though.)

I thought that functional programming could get you both the brevity of
Python's slicing and the performance of C's compilation but I was wrong.
You can get the brevity of the Python approach in C.

You might have guessed that you can get rid of most temporary arrays,
at the expense of readability, with numpy. For example, replacing

d1[:] = odd[:] - C1*even[:]

with

numpy.multiply(-C1,even,d1)
numpy.add(d1,odd,d1)

eliminates the intermediate. No unnecessary array allocations; no
closures necessary.

Right. Performance will probably go from 5x slower to 2x slower because
you're traversing the arrays twice instead of once.

I'm curious whether this shared-slicing isn't, in some ways,
advantageous over the Ocaml way.

That's exactly what I thought to start with. The author of F# is working on
getting slicing into his language but the use of slicing in this Python
benchmark corrupted my fragile little mind. Slicing is a terrible way to
approach this problem if you're using a compiled language like F#.

I first wrote an OCaml translation of the Python and wrote my own
little "slice" implementation. I have since looked up a C++ solution and
translated that into OCaml instead:

let rec d4_aux a n =
let n2 = n lsr 1 in
let tmp = Array.make n 0. in
for i=0 to n2-2 do
tmp.(i) <- a.(i*2)*.h0+.a.(i*2+1)*.h1+.a.(i*2+2)*.h2+.a.(i*2+3)*.h3;
tmp.(i+n2) <- a.(i*2)*.g0+.a.(i*2+1)*.g1+.a.(i*2+2)*.g2+.a.(i*2+3)*.g3;
done;
tmp.(n2-1) <- a.(n-2)*.h0 +. a.(n-1)*.h1 +. a.(0)*.h2 +. a.(1)*.h3;
tmp.(2*n2-1) <- a.(n-2)*.g0 +. a.(n-1)*.g1 +. a.(0)*.g2 +. a.(1)*.g3;
Array.blit tmp 0 a 0 n;
if n > 4 then d4_aux a (n lsr 1)

let d4 a = d4_aux a (Array.length a)

Not only is that shorter than the Python, it is much faster:

0.56s C++ (direct arrays)
0.61s F# (direct arrays)
0.62s OCaml (direct arrays)
1.38s OCaml (slices)
2.38s Python (slices)
10s Mathematica 5.1

Note that all implementations are safe (e.g. C++ uses a.at(i) instead of
a[i]).

I presume that in Ocaml, the way
you'd "share" array data is to create a closure can apply some
operation to selected elements of the array, returning the result.

Yes. That is certainly one way of doing it.

But could this as easily modify the array in-place?

Absolutely. A closure would capture a reference to the array. Arrays are
mutable so you could alter the array from within the closure. In F# you
could even execute the closures concurrently to alter different parts of
the same array at the same time. I just tried that and it is actually
slower to multithread this.

I presume a good
compiler could be able to inline the function and optimize the call to
an in-place operation, but I'm not sure how well this works in
practice.

Surprisingly well it seems: F# and OCaml are almost as fast as C/C++!

Here's a version of D4_Transform that uses no temporary arrays (aside
from the work arrays s1, d1, and d2, which are allocated only once).
It's about 40% faster for me than the one with infix operations. I'd
be curious how it compares to a correct Ocaml version. (I'd still
expect Ocaml to be at least twice as fast.)

I get:

1.57s Python (in-place)

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.)

Writing a JIT to compile this kind of stuff is easy.

Eh, getting a JIT to do the optimization I spoke of in Python is not
easy. It would be relatively easy in a statically typed language. In
Python, it'd be tantamount to rewriting a dynamically reconfigurable
version of numpy--you'd have to duplicate numpy's complex
type-dispatching rules.

There aren't any complicated types in the above code. In fact, there are
only two types: float and float array. Type checking is easy in this case,
compilation to C is also easy, then you just dispatch to the compiled C
code when the types are ok. You would want to write your Python code like C
code but that is shorter in this case. You may also want to flag code for
compilation manually in order to avoid the overhead of compiling code
unnecessarily.

My point is that this is fundamentally bad code,

Whoa, there. I realize that for people who prefer functional
programming, with their recursively reductionist way of thinking, it
might seem as if leaving any sort of reducibility in final result is
"fundamentally bad", but that's a pretty strong thing to say.

I was referring to the slicing specifically, nothing to do with functional
programming. My F# and OCaml code are now basically identical to the C++
code.

so why bother trying to write a Python JIT? Why
not just write in a better language for this task? Optimising within a
fundamentally slow language seems silly to me.

Because, for most people, language choice is not a greedy maximizing of
a single issue. Nobody who uses numpy is under the impression that it
can match a statically-typed language that is compiled to machine code
in peak performance. (Matlab is fair game, of course.) But speed is
not the only important thing.

In this specific context (discrete wavelet transform benchmark), I'd have
said that speed was the most important thing after correctness.

I use Python because it's a excellently designed language that fits my
overall needs, and numpy because sometimes I need more speed than
vanilla Python. Only when speed is critical (for example, 3D collision
detection) do I write an extension in a "better language for the task"
(i.e. C).

Have you looked at languages like F#, OCaml, Haskell, SML and so on? They
seem to offer the benefits of both worlds.

This is something that's quite popular in the numerically-intensive
computing community, BTW. Many people use Python to handle boring
stuff like file I/O, memory managment, and non-critical numerical
calculations, and write C or Fortran extensions to do the
numerically-intensive stuff.

Yes, I appreciate that Python is hugely popular, I just don't understand why
it is quite so popular.

I found a series of lectures on VPython, showing some amazing stuff:

http://showmedo.com/videos/series?name=pythonThompsonVPythonSeries

I just found some more interesting stuff here:

http://www.phy.syr.edu/~salgado/software/vpython/

This is really great work but I can't help but wonder why the authors chose
to use Python when other languages seem better suited. I'd like to work on
raising people's awareness of these alternatives, and probably create some
useful tools in the process. So I'm keen to learn what Python programmers
would want/expect from F# and OCaml.

What would it take to make you convert?

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

Relevant Pages

  • Re: Article of interest: Python pros/cons for the enterprise
    ... Your product gets the benefit of the doubt if it was written in Python, but not if written in C? ... In Python, somebody could still define a monolithic array, use different parts of it to represent different things, and then "oops" lose track of their index, thereby overwriting one part of the array with data that should have stayed in the other part. ... you mentioned a "failure to use best engineering practices." ... Even if you could somehow remove every feature open to potential abuse, and thereby dodo-proof the language, the result would be a slow, weak, ugly language that non-dodos would not want to use. ...
    (comp.lang.python)
  • Its still I, Miville
    ... array, where absence of key would mean not absence of element but ... Python is an evolving language, ... sake of legal definition Microsoft is liable to give the object-oriented ...
    (comp.lang.python)
  • (patch for Bash) regex(3) splitting/matching
    ... I usually do this in Python. ... 'help array' will give you more info on other options for 'array' ... int dollarflag, zeropad, compareflag; ... SHELL_VAR *var; ...
    (comp.unix.shell)
  • Re: Finding the instance reference of an object
    ... with a large array will copy the entire array. ... reference model, Python seems rather inconsistent: ... you assume it is some other language. ...
    (comp.lang.python)
  • Re: Loopless syntax for 2d in NumPy (or Numarray)
    ... >> source array and would contain a running list of highest values of ... Learning Python, Python Essential Reference, and Practical Python. ... I have a program that was written in a programming language called SPL ... explicitly declare loops to accomplish a lot of what I'm doing. ...
    (comp.lang.python)