Numerical Python
For the past few months, I've been covering different software packages for scientific computations. For my next several articles, I'm going to be focusing on using Python to come up with your own algorithms for your scientific problems. Python seems to be completely taking over the scientific communities for developing new code, so it is a good idea to have a solid working knowledge of the fundamentals so you can build solutions to your own problems.
In this article, I start with NumPY. My next article will cover SciPy, and then I'll look at some of the more complicated modules available in the following article.
So, let's start with the root module from which almost all other scientific modules are built, NumPY. Out of the box, Python supports real numbers and integers. You also can create complicated data structures with lists, sets and so on. This makes it very easy to write algorithms to solve scientific problems. But, just diving in naively, without paying attention to what is happening under the hood, will lead to inefficient code. This is true with all programming languages, not just Python. Most scientific code needs to squeeze every last available cycle out of the hardware. One of the things to remember about Python is that it is a dynamic language where almost all functions and operators are polymorphic. This means that Python doesn't really know what needs to be done, at a hardware level, until it hits that operation. Unfortunately, this rules out any optimizations that can be made by rearranging operations to take advantage of how they are stored in memory and cache.
One property of Python that causes a much larger problem is polymorphism. In this case, Python needs to check the operands of any operator or function to see what type it is, decide whether this particular operand or function can handle these data types, then use the correct form of the operand or function to do the actual operation. In most cases, this is not really an issue because modern computers have become so fast. But in many scientific algorithms, you end up applying the same operations to thousands, or millions, of data points. A simple example is just taking the square of the first 100,000 numbers:
import time
a = range(100000)
c = []
starttime = time.clock()
for b in a:
c.append(b*b)
endtime = time.clock()
print "Total time for loop: ", (endtime  starttime)
This little program uses the range
function to create a list of the
first 100,000 integers. You need to import the time
module to get access
to the system clock to do timing measurements. Running this on my desktop
takes approximately 0.037775 seconds (remember always to make several
measurements, and take the average). It takes this much time because
for every iteration of the loop, Python needs to check the type of the
b variable before trying to use the multiplication operator. What can
you do if you have even more complicated algorithms?
This is where NumPY comes in. The key element that NumPY introduces is an Ndimensional array object. The great flexibility of Python lists, allowing all sorts of different types of elements, comes at a computational cost. NumPY arrays deal with this cost by introducing restrictions. Arrays can be multidimensional, and they must all be the same data type. Once this is done, you can start to take some shortcuts by taking advantage of the fact that you already know what the type of the elements is. It adds extra overloading functions for the common operators and functions to help optimize uses of arrays.
All of the
normal arithmetic operators work on NumPY arrays in an elementwise
fashion. So, to get the squares of the elements in an array, it would
be as simple as array1 * array1
.
NumPY also uses external standard, optimized libraries written in C or FORTRAN to handle many of the actual manipulations on these array data types. This is handled by libraries like BLAS or lapack. Python simply does an initial check of each of the arrays in question, then hands them as a single object to the external library. The external library does all of the hard work, then hands back a single object containing the result. This removes the need for Python to check each element when using the loop notation above. Using NumPY, the earlier example becomes:
import numpy
import time
a = numpy.arange(1000000)
starttime = time.clock()
c = a * a
endtime = time.clock()
print "Total time used: ", (endtime  starttime)
Running this toy code results in an average run time of 0.011167 seconds for this squaring operation. So the time is cut by one third, and the readability of the code is simplified by getting rid of the loop construct.
I've dealt only with onedimensional arrays so far, but NumPY supports multidimensional arrays just as easily. If you want to define a twodimensional array, or a matrix, you can set up a small one with something like this:
a = numpy.array([[1,2,3,4], [2,3,4,5]])
Basically, you are creating a list of lists, where each of the sublists
is each of the rows of your matrix. This will work only if each of the
sublists is the same size—that is, each of the rows has the same number
of columns. You can get the total number of elements in this matrix,
with the property a.size
. The dimensions of the matrix are stored in
the property a.shape
. In this case, the size is 8, and the shape is
(2, 4), or two rows and four columns. What shape did the array in the
earlier example have? If you go ahead and check, you should see that
the shape is (1000000). The other key properties of these arrays are:

ndim: the number of dimensions.

dtype: the data type of the elements.

itemsize: the size in bytes of each element.

data: the buffer that stores the actual data.
You also can reshape arrays. So if you wanted to take the earlier example of the first 100,000 integers and turn it into a threedimensional array, you could do something like this:
old_array = numpy.arange(100000)
new_array = old_array.reshape(10,100,100)
This will give you a new 3D array laid out into a 10x100x100 element cube.
Now, let's look at some of the other functions available to apply to
arrays. If you remember from earlier, all of the standard arithmetic
operations are overloaded to operate on arrays one element at a time. But,
most matrix programming languages use the multiplication element to mean
matrix multiplication. This is something to keep in mind when you start
using Python. To get a true matrix multiplication, you need to use the
dot()
function.
If you have two matrices, A and B, you can multiply
them with numpy.dot(A, B)
.
Many of the standard mathematical functions,
like cosine, sine, square root and so on, are provided by NumPY and are
called universal functions. Just like with the arithmetic operators,
these universal functions are applied elementwise across the array. In
science, several common functions are used. You can get
the transpose of a matrix with the object function
a.transpose()
. If
you need to get the inverse of a matrix, there is the module function
inv()
, so you would execute
numpy.inv(a)
. The trace is also a module
function, given by numpy.trace(a)
.
Even more complicated
functions are available. You can solve systems of equations with NumPY. If
you have a matrix of coefficients given by a, and a vector of numbers
representing the righthand side of your equations given by y, you can
solve this system with numpy.solve(a,y)
. In many
situations, you may
be interested in finding the eigenvalues and eigenfunctions of a given
system. If so, you can use numpy.eig(array1)
to get
those values.
The last thing I want to look at here is a class that provides even more shortcuts, at the cost of more restrictions. Matrices (2D arrays) are so prevalent that NumPY provides a special class to optimize operations using them as much as possible. To create a new matrix, say a 2x2 matrix, you would write:
A = numpy.matrix('1.0, 2.0; 3.0, 4.0')
Now, you can get the transpose with just A.T
. Similarly, the inverse
is found with A.I
. The multiplication operation will do actual matrix
multiplication when you use matrix objects. So, given two matrix objects
A and B, you can do matrix multiplication with A*B
. The solve function
still works as expected on systems of equations that are defined using
matrix objects. Lots of tips and tricks available on the NumPY
Web site, which is well worth a look, especially as you start out.
This short introduction should get you started in thinking of Python as a viable possibility in "real" numerical computations. The NumPY module provides a very strong foundation to build up complex scientific workflows. Next month, I'll look at one of the available modules, SciPY. Until then, play with all of the raw numbercrunching possibilities provided by NumPY.
Joey Bernard has a background in both physics and computer science. This serves him well in his day job as a computational research consultant at the University of New Brunswick. He also teaches computational physics and parallel programming.
Linux has become a key foundation for supporting today's rapidly growing IT environments. Linux is being used to deploy business applications and databases, trading on its reputation as a lowcost operating environment. For many IT organizations, Linux is a mainstay for deploying Web servers and has evolved from handling basic file, print, and utility workloads to running missioncritical applications and databases, physically, virtually, and in the cloud. As Linux grows in importance in terms of value to the business, managing Linux environments to high standards of service quality — availability, security, and performance — becomes an essential requirement for business success.
Sponsored by Red Hat
If you already use virtualized infrastructure, you are well on your way to leveraging the power of the cloud. Virtualization offers the promise of limitless resources, but how do you manage that scalability when your DevOps team doesn’t scale? In today’s hypercompetitive markets, fast results can make a difference between leading the pack vs. obsolescence. Organizations need more benefits from cloud computing than just raw resources. They need agility, flexibility, convenience, ROI, and control.
Stackato private PlatformasaService technology from ActiveState extends your private cloud infrastructure by creating a private PaaS to provide ondemand availability, flexibility, control, and ultimately, faster timetomarket for your enterprise.
Sponsored by ActiveState
Trending Topics
diff u: What's New in Kernel Development  Jul 23, 2014 
Great Scott! It's Version 13!  Jul 21, 2014 
Adminer—Better Than Awesome!  Jul 17, 2014 
It Actually Is Rocket Science  Jul 16, 2014 
Android Candy: Repix, Not Just Another Photo App  Jul 14, 2014 
Wanted: Your Embedded Linux Projects  Jul 10, 2014 
 diff u: What's New in Kernel Development
 Divx# Watch The Other Woman Full HD Online Streaming Viooz
 Numerical Python
 Use Linux as a SAN Provider
 NSA: Linux Journal is an "extremist forum" and its readers get flagged for extra surveillance
 Great Scott! It's Version 13!
 RSS Feeds
 Tech Tip: Really Simple HTTP Server with Python
 Adminer—Better Than Awesome!
 <Watch> HD! Watch Walking On Sunshine Online Full Movie Streaming
Featured Jobs
Linux Systems Administrator  Houston and Austin, Texas  Host Gator 
Senior Perl Developer  Austin, Texas  Host Gator 
Technical Support Rep  Houston and Austin, Texas  Host Gator 
UX Designer  Austin, Texas  Host Gator 
Web & UI Developer (JavaScript & j Query)  Austin, Texas  Host Gator 