Get the code for this example The zip function in Python accepts a number of sequences (such as lists or arrays) and returns a list of tuples.  The first element of each tuple comes from the first sequence, the second from the second sequence, and so on.  For example: [code language=“Python”]a = range(0, 5) b = range(5, 10) c = range(10, 15) sequence_of_tuples = zip(a, b, c) print(sequence_of_tuples)[/code]

[(0, 5, 10), (1, 6, 11), (2, 7, 12), (3, 8, 13), (4, 9, 14)]

I’ve just released a self-contained Fortran module that solves a system of linear equations using the LU decomposition. Download the Fortran linear solver from github This module is based on code that was implemented and released on the Web by Jean-Pierre Moreau.  His implementation was based on one of the Numerical Recipes books.  I updated his code to a more strict Fortran 90 standard and added the necessary comments so that it can be built as a Python module using f2py.  I replaced Jean-Pierre’s Fortran test program with a simpler, self-contained program.  I also included a Python script that implements the same test case. I created this module because sometimes a self-contained routine is more appropriate than a full library.  I am compiling a library that implements a custom boundary condition for a proprietary computational fluid dynamics solver (CFD-ACE+).  The library has to be written in Fortran, and it has to be built using a proprietary set of build scripts.  I could either try to reverse-engineer the build process and to modify it to link to a shared library like LAPACK, or I could implement a self-contained solver.  Since Pierre had already implemented the solver, I was able to slightly modify his code and get it working relatively quickly.

I’ve recently gained a lot of experience with Sage, an open-source alternative to MATLAB, Mathematica, Maple, MuPAD, and Magma.  Here are a couple of links to check out:

Public notebook servers--try it online right now! Interactive examples with Sage

• Mathematical publications citing Sage

Sage vs. Mathematica

My experience with mathematical software started my freshman year at the University of Illinois when I signed up for a calculus class that was taught almost entirely with Mathematica.  I grew to love Mathematica’s symbolic computation and plotting capabilities, although I found its programming language to be cumbersome.  Once I was no longer a student, Mathematica was no longer an option due to the hefty licensing fees.  With the caveat that I haven’t used Mathematica in many years,  I will say that Sage compares very favorably to my experience with Mathematica.  Sage has a notebook interface that allows you to integrate code, results, text, typeset equations, and graphics in an interactive document that can be viewed with any standards-compliant web browser.  I don’t know how Sage stacks up against Mathematica in areas of advanced mathematics.  I have heard that Sage is far ahead of Mathematica in number theory, because the lead developer of Sage is a number theorist.

Linear time-invariant (LTI) systems are widely used in the field of signal processing.  Scipy contains powerful tools for simulating LTI systems in the scipy.signal package, but they are not well documented.  I will provide a simple example that demonstrates how to use a few of the core classes and functions in scipy.signal for simulating LTI systems with Python.

Define an LTI system

You will need to have Scipy installed, and you will need to have Matplotlib as well to make plots. We will start with an example of a first-order LTI system, which is characterized by a single parameter known as the time constant.  scipy.signal defines the class lti to represent a linear system in Python.

Scipy contains a good least-squares fitting routine, leastsq(), which implements a modified Levenberg-Marquardt algorithm.  I just learned that it also has a constrained least-squared routine called fmin_slsqp().   I am using simple upper and lower bound constraints, but it’s also possible to specify more complex functional constraints. What I did not realize, at first, is that fmin_slsqprequires a different type of objective function thanleastsq.   leastsqrequires you to write a function that returns a vector of residuals, andleastsqautomatically squares and sums the residuals.fmin_slsqp is actually more flexible, in that it can use any objective function that returns a single scalar value.  To implement least-squares curve fitting, your objective function will need to find the residual at each data point, square the values, and sum them up.  Hopefully this tip will save you some time. Check out the scipy optimization tutorial for more examples.  Here is the original paper by Dieter Kraft which introduces the algorithm used by fmin_slsqp.

In a previous post, I described how Python’s Pickle module is fast and convenient for storing all sorts of data on disk. More recently, I showed how to profile the memory usage of Python code.  In recent weeks, I’ve uncovered a serious limitation in the Pickle module when storing large amounts of data: Pickle requires a large amount of memory to save a data structure to disk. Fortunately, there is an open standard called HDF, which defines a binary file format that is designed to efficiently store large scientific data sets. I will demonstrate both approaches, and profile them to see how much memory is required. I am writing the HDF file using the PyTables interface. Here’s the little test program I’ve been using: [code language=“Python”]#!/usr/bin/env python from numpy import array

I  have recently started using f2py to call Fortran from Python.  I have found this useful for two reasons: speeding up Python scripts by calling compiled Fortran code, and using Python as a unit testing framework for Fortran modules.   Unfortunately, the documentation for f2py is rather sparse, and may not be completely up to date.   In this note, I will hopefully prevent you from wasting a lot of time figuring out how to pass array arguments, and return array results.

In a previous post, I explained how to use the Python profiler.  The profile is great for finding out which parts of the code run the slowest, or are called most often.  However, the profiler doesn’t give any information about how much RAM is being consumed, or where it’s being consumed.  If your program needs so much memory that it starts swapping to disk, its speed can be reduced by orders of magnitude.  On the positive side, your code may run much faster if it fits entirely in the processor cache.  In this post, I will introduce two tools that can help you understand the RAM usage of your Python code.

Lookup tables and spline fitting are widely used by scientific programmers.  A particular function may not have an analytic solution–in other words, it can’t be expressed as an equation of elementary functions.  This might happen if the function were empirically determined from experimental data, or if the equation can’t be algebraically solved for one variable.  If the analytical function is available, but takes a long time to evaluate, a lookup table or spline approximation can be considerably faster.  In a previous post, I showed how to use the function interp1d from scipy.interpolate as a lookup table.  In a later post, I showed that interp1d is actually rather slow, and scipy.interpolate.UnivariateSpline is much faster.  Now, I will show some benchmark results, and explain a potential pitfall when using UnivariateSpline. The two functions I will use for this demonstration are not very complicated.  Each function is a piecewise approximation of a more complex function.  Here is one of the functions.  For small z, the function uses one approximation, and for large z, it uses a far-field approximation.  For intermediate values, a polynomial is used to smoothly match the two solutions.