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.
Trending Topics
Webinar
Practical Task Scheduling Deployment
July 20, 2016 12:00 pm CDT
One of the best things about the UNIX environment (aside from being stable and efficient) is the vast array of software tools available to help you do your job. Traditionally, a UNIX tool does only one thing, but does that one thing very well. For example, grep is very easy to use and can search vast amounts of data quickly. The find tool can find a particular file or files based on all kinds of criteria. It's pretty easy to string these tools together to build even more powerful tools, such as a tool that finds all of the .log files in the /home directory and searches each one for a particular entry. This erectorset mentality allows UNIX system administrators to seem to always have the right tool for the job.
Cron traditionally has been considered another such a tool for job scheduling, but is it enough? This webinar considers that very question. The first part builds on a previous Geek Guide, Beyond Cron, and briefly describes how to know when it might be time to consider upgrading your job scheduling infrastructure. The second part presents an actual planning and implementation framework.
Join Linux Journal's Mike Diehl and Pat Cameron of Help Systems.
Free to Linux Journal readers.
Register Now!SUSE LLC's SUSE Manager  Jul 21, 2016 
My +1 Sword of Productivity  Jul 20, 2016 
NonLinux FOSS: Caffeine!  Jul 19, 2016 
Murat Yener and Onur Dundar's Expert Android Studio (Wrox)  Jul 18, 2016 
Rogue Wave Software's Zend Server  Jul 14, 2016 
Webinar: Practical Task Scheduling Deployment  Jul 14, 2016 
 SUSE LLC's SUSE Manager
 My +1 Sword of Productivity
 Murat Yener and Onur Dundar's Expert Android Studio (Wrox)
 Managing Linux Using Puppet
 NonLinux FOSS: Caffeine!
 Doing for User Space What We Did for Kernel Space
 SuperTuxKart 0.9.2 Released
 Parsing an RSS News Feed with a Bash Script
 Google's SwiftShader Released
 Rogue Wave Software's Zend Server
Geek Guides
With all the industry talk about the benefits of Linux on Power and all the performance advantages offered by its open architecture, you may be considering a move in that direction. If you are thinking about analytics, big data and cloud computing, you would be right to evaluate Power. The idea of using commodity x86 hardware and replacing it every three years is an outdated cost model. It doesn’t consider the total cost of ownership, and it doesn’t consider the advantage of real processing power, highavailability and multithreading like a demon.
This ebook takes a look at some of the practical applications of the Linux on Power platform and ways you might bring all the performance power of this open architecture to bear for your organization. There are no smoke and mirrors here—just hard, cold, empirical evidence provided by independent sources. I also consider some innovative ways Linux on Power will be used in the future.
Get the Guide