About alternatives to Matlab

[Filip Wasilewski]

sturlamolden wrote:
> Boris wrote:
> > Hi, is there any alternative software for Matlab? Although Matlab is
> > powerful & popular among mathematical & engineering guys, it still
> > costs too much & not publicly open. So I wonder if there's similar
> > software/lang that is open & with comparable functionality, at least
> > for numerical aspect. Thanks!
>
> I have used Matlab for years, and has recently changed to Python. In
> addition one needs NumPy and Matplotlib, and perhaps SciPy. Other
> useful packages are PyGTK for GUI, PyGame for multimedia, PyOpenGL for
> 3D graphics, Mpi4Py for parallel computation, etc. You will find python
> packages for nearly any conceivable task. Unlike Matlab, Python is a
> general purpose programming language, and a distant cousin of Lisp. The
> Python language is fare more expressive and productive than Matlab, yet
> even more easy to use.

Huh, so you don't think it's just damn sexy to do OOP, networking and
databases by multiplying matrices?;-)


> The NumPy package is the core requirement for numerical work in Python.
> It is quite different form Matlab, but I think it is more powerful.
> Particularly, arrays are passed by reference (not by value), and
> indexing creates view matrices.

For those interested in drilling down on this topic there is a Numpy
for Matlab Users guide at http://scipy.com/NumPy_for_Matlab_Users


> To compare Matlab with NumPy we can e.g. use the D4 discrete wavelet
> transform. I have here coded it in Matlab and Python/NumPy using Tim
> Swelden's lifting scheme.
[...]

Actually you have not. The algorithm you presented gives completely
wrong results. Have a look at quick&dirty(TM) implementation bellow.

import math
import numpy

def d2_lwt(x, axis=-1, level=None, copy=False):
    """
    Daubechies2 (4 point support) Lifting Wavelet Transform in NumPy
    """
    C1 = 1.7320508075688772    # sqrt(3)
    C2 = 0.4330127018922193    # sqrt(3)/4
    C3 = -0.066987298107780702 # (sqrt(3)-2)/4)
    C4 = 0.51763809020504137   # (sqrt(3)-1)/sqrt(2)
    C5 = 1.9318516525781364    # (sqrt(3)+1)/sqrt(2)

    if not isinstance(x, numpy.ndarray) or x.dtype.kind != 'f':
        x = numpy.array(x, dtype=numpy.float64)
    elif copy:
        x = x.copy()

    max_level = int(math.floor(math.log(x.shape[axis],2)))
    if level is None:
        level = max_level
    else:
        assert level <= max_level, "Level param too high"

    coeffs = []
    cA = x.swapaxes(0, axis)
    while level > 0:
        level -= 1

        # lazy
        even = cA[::2]
        odd = cA[1::2]

        # dual
        # using `even = even + C1*odd` may speed up things on
        # some machines (only for not in-place transform).
        even += C1*odd

        # primal
        odd[0] -= C2*even[0] + C3*even[-1]
        odd[1:] -= C2*even[1:] + C3*even[:-1]

        # dual
        even[:-1] -= odd[1:]
        even[-1] -= odd[0]

        # scale
        even *= C4
        odd *= C5

        cA, cD = even, odd
        coeffs.append(cD.swapaxes(axis, 0))

    coeffs.append(cA.swapaxes(axis, 0))
    coeffs.reverse()
    return coeffs

if __name__ == "__main__":
    d = [1,5,3,2,4,8,5,2]
    data = [
        numpy.array(d, dtype=numpy.float32),
        numpy.array(d, dtype=numpy.float64),
        numpy.array([d]*2, dtype=numpy.float64),
        numpy.array([d]*2, dtype=numpy.float32),
        numpy.array(d, dtype=numpy.int32),
        [d]*2,
        [[d,d]]*3,
    ]

    for i,x in enumerate(data):
        print "Case %d:" % (i+1)
        print "x in:\n", x
        coeffs = d2_lwt(x, axis=-1, level=None, copy=False)
        print "coeffs:"
        for c in coeffs:
            print c, c.dtype
        print "x out:\n", x
        print


Please excuse me not including Matlab version here, I would like to
have a free weekend.

As far as the speed comparison is concerned I totally agree that NumPy
can easily outperform Matlab in most cases. Of course one can use
compiled low-level extensions to speed up specific computations in
Matlab, but it's a lot easier and/or cheaper to find very good tools
for Python.


> If anyone wonders why I think Travis Oliphant and the NumPy team should
> be knighted, then this is the answer.

Would love to see the ceremony:)

cheers,
fw