About alternatives to Matlab

[Mark Morss]

I doubt that anyone would dispute that even as boosted by Numpy/Scipy,
Python will almost certainly be notably slower than moderately
well-written code in a compiled language.  The reason Numpy exists,
however, is not to deliver the best possible speed, but to deliver
enough speed to make it possible to solve large numerical problems with
the powerful and flexible Python language.  As observed by Hans
Latangen in _Python Scripting for Computational Science_, 2nd ed.,
Springer 2005, scientific computing is more than number crunching:

"Very often programming is about shuffling data in and out of different
tools, converting one data format to another, extracting numerical data
from a text, and administering numerical experiments involving a large
number of data files and directories.  Such tasks are much faster to
accomplish in a language like Python than in Fortran, C, C++ or Java."

He might well have mentioned the importance of developing nice-looking
reports once the analysis is complete, and that development is simpler
and more flexible in Python than a compiled language.

In principle, I agree that heavy-duty number-crunching, at least if it
has to be repeated again and again, should be accomplished by means of
a compiled language.  However, if someone has to solve many different
problems just one time, or just a few times (for example if you are an
engineering consultant), there is an excellent argument for using
Python + Numpy.  Unless it affects feasibility, I opine, computational
speed is important primarily in context of regular production, e.g.,
computing the daily "value at risk" in a commodity trading portfolio,
or making daily weather predictions.

I am aware of the power and flexibility of the OCaml language, and
particularly that an OCaml user can easily switch back and forth
between interpreted and compiled implementation.  I'm attacted to OCaml
and, indeed, I'm in the process of reading Smith's (unfortunately not
very well-written) _Practical OCaml_.  However, I also understand that
OCaml supports only double-precision implementation of real numbers;
that its implementation of arrays is a little clunky compared to
Fortran 95 or Numpy (and I suspect not as fast as Fortran's); and that
the available libraries, while powerful, are by no means as
comprehensive as those available for Python.  For example, I am unaware
of the existance of an HDF5 interface for OCaml.

In summary, I think that OCaml is an excellent language, but I think
that the question of whether to use it in preference to Python+Numpy
for general-purpose numerical analysis must rest on much more than its
computational speed.

Jon Harrop wrote:
> Filip Wasilewski wrote:
> > Besides of that this code is irrelevant to the original one and your
> > further conclusions may not be perfectly correct. Please learn first
> > about the topic of your benchmark and different variants of wavelet
> > transform, namely difference between lifting scheme and dwt, and try
> > posting some relevant examples or use other topic for your benchmarks.
>
> Your lifting implementation is slightly longer and slower than the
> non-lifting one but its characteristics are identical. For n=2^25 I get:
>
> 1.88s C++ (816 non-whitespace bytes)
> 2.00s OCaml (741 b)
> 2.33s F# (741 b)
> 9.83s Your Python (1,002 b)
>
> The original python was 789 bytes.
>
> > Neglecting this issues, I wonder if it is equally easy to juggle
> > arbitrary multidimensional arrays in a statically typed language and do
> > operations on selected axis directly from the interactive interpreter
> > like in the NumPy example from my other post -
> > http://groups.google.com/group/comp.lang.python/msg/a71a5bf00a88ad57 ?
> > I don't know much OCaml so it would be great if you could elaborate on
> > this.
>
> It is probably just as easy. Instead of dynamic typing you have parametric
> polymorphism. If you want to make your function generic over arithmetic
> type then you can pass in the arithmetic operators.
>
> >> 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]).
> >
> > So why not use C++ instead of all others if the speed really matters?
> > What is your point here?
>
> 1. Benchmarks should not just include two inappropriate languages.
> 2. Speed aside, the other languages are more concise.
>
> > Could you please share full benchmark code, so we could make
> > conclusions too?
>
> I'll paste the whole programs at the end...
>
> > I would like to get to know about your benchmark
> > methodology. I wonder if you are allocating the input data really
> > dynamically or whether it's size is a priori knowledge for the
> > compiler.
>
> The knowledge is there for the compiler to use but I don't believe any of
> them exploit it.
>
> >> In this specific context (discrete wavelet transform benchmark), I'd have
> >> said that speed was the most important thing after correctness.
> >
> > Not necessarily, if you are doing research with different kinds of
> > wavelets and need a general, flexible and easily adaptable tool or just
> > the computation speed is not considered a bottleneck.
>
> Brevity is probably next.
>
> > Language speed is a great advantage, but if it always was the most
> > important, people would talk directly to the CPUs or knit DSP chips in
> > theirs backyards whenever they need to solve a calculation problem.
>
> Sure.
>
> > Please don't make this discussion heading towards `what is the fastest
> > way of doing d4 transform and why OCaml rules` instead of `what is the
> > optimal way of doing things under given circumstances and still have a
> > free afternoon ;-)`.
>
> Comments like that seem to be based on the fundamental misconception that
> writing C++ like this is hard.
>
> > Different tasks need different, well-balanced
> > measures. Neither Python nor OCalm nor any other language is a panacea
> > for every single problem.
>
> Absolutely. OCaml is as good as the next (compiled) language in this case.
> Python and Matlab seem to be particularly bad at this task.
>
> Here's my complete OCaml:
>
> let a = sqrt 3. and b = 4. *. sqrt 2.
> let h0, h1, h2, h3 =
>   (1. +. a) /. b, (3. +. a) /. b, (3. -. a) /. b, (1. -. a) /. b
> let g0, g1, g2, g3 = h3, -.h2, h1, -.h0
>
> let rec d4_aux tmp a n =
>   let n2 = n lsr 1 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 tmp a (n lsr 1)
>
> let d4 a =
>   let tmp = Array.make (Array.length a) 0. in
>   d4_aux tmp a (Array.length a)
>
> let () =
>   print_endline "Generating test data...";
>   let x = Array.init (1 lsl 25) (fun _ -> Random.float 1.) in
>   print_endline "Benchmarking...";
>   let t1 = Sys.time() in
>   ignore(d4 x);
>   Printf.printf "Elapsed time is %.6f seconds\n" (Sys.time() -. t1)
>
> and C++:
>
> #include <iostream>
> #include <vector>
> #include <ctime>
> #include <cmath>
>
> using namespace std;
>
> double a = sqrt(3), b = 4 * sqrt(2);
> double h0 = (1 + a) / b, h1 = (3 + a) / b, h2 = (3 - a) / b, h3 = (1 - a) /
> b;
> double g0 = h3, g1 = -h2, g2 = h1, g3 = -h0;
>
> void d4(vector<double> &a) {
>   vector<double> tmp(a.size());
>   for (int n = a.size(); n >= 4; n >>= 1) {
>     int half = n >> 1;
>
>     for (int i=0; i<n-2; i+=2) {
>       tmp.at(i/2) = a.at(i)*h0 + a.at(i+1)*h1 + a.at(i+2)*h2 + a.at(i+3)*h3;
>       tmp.at(i/2+half) = a.at(i)*h3 - a.at(i+1)*h2 + a.at(i+2)*h1 -
> a.at(i+3)*h0
> ;
>     }
>
>     tmp.at(half-1) = a.at(n-2)*h0 + a.at(n-1)*h1 + a.at(0)*h2 + a.at(1)*h3;
>     tmp.at(2*half-1) = a.at(n-2)*h3 - a.at(n-1)*h2 + a.at(0)*h1 -
> a.at(1)*h0;
>
>     copy(tmp.begin(), tmp.begin() + n, a.begin());
>   }
> }
>
> int main() {
>   cout << "Generating data...\n";
>   vector<double> a(1 << 25);
>   for (int i=0; i<a.size(); ++i) a.at(i) = rand();
>   cout << "Benchmarking...\n";
>   double t1 = clock();
>   d4(a);
>   cout << "Elapsed time " << (clock() - t1) / CLOCKS_PER_SEC << "s\n";
> }
>
> --
> Dr Jon D Harrop, Flying Frog Consultancy
> Objective CAML for Scientists
> http://www.ffconsultancy.com/products/ocaml_for_scientists/index.html?usenet