Numba and LLVMPy

It's been a busy year so far.    All the time spent on starting a new company, starting new open source projects, and keeping up with the open source projects that I have interest in, has meant that I haven't written nearly as many blog-posts as I planned on.   But, this is probably a good thing at least if you follow the wisdom attributed to Solomon --- which has been paraphrased in this quote attributed to Abraham Lincoln.

One of the things that has been on my mind for the past year is promoting array-oriented computing as a fundamental concept more developers need exposure to.    This is one reason that I am so excited that I've been able to find great people to work on Numba (which intends to be an array-oriented compiler for Python code).      I have given a few talks trying to convey what is meant by array-oriented computing, but the essence is captured by the difference between the life.py example in the Python code-base and a NumPy version of the same code.

I have seen many, many real world examples of very complicated code that could be simplified and sped up (especially on modern hardware) by just thinking about the problem differently using array-oriented concepts. 

One of the goals for Numba is to make it possible to write more vectorized code easily in Python without relying just on the pre-compiled loops that NumPy provides.    In order to write Numba, though, we first needed to resurrect the llvm-py project which provides easy access to the LLVM C++ libraries from Python.   This project is interesting in its own right and in addition to forming a base tool chain for Numba, allows you to do very interesting things (like instrument C-code compiled with Clang to bitcode), build a compiler, or import bitcode directly into Python (a la bitey).

While the documentation for llvm-py left me frustrated early on, I have to admit that llvm-py re-kindled some of the joy I experienced when being first exposed to Python.    Over the past several weeks we have worked to create the llvmpy project from llvm-py.   We now have a domain http://www.llvmpy.org, a GitHub repository, a website served from GitHub, and sphinx-based documents that can be edited via a pull request.    The documentation still needs a lot of improvement (even to get it to the state that the old llvm-py project was in), and contributions are welcome.  

I'm grateful to Fernando Perez, author of IPython, for explaining the 4-repository approach to managing an open source web-site and documentation via github.   We are using the same pattern that IPython uses for both numba and llvmpy.   It took a bit of work to get set-up but it's a nice approach that should make it easier for the community to maintain the documentation and web-site of both of these projects.    The idea is simple.   Use a project page (repo llvmpy.github.com) to be the web-site but generate this repo from another repo (llvmpy-webpage) which contains the actual sources.   I borrowed the scripts from the IPython project to build the page from the sources, check-out the llvmpy.github.com repo, copy the built pages to the repo, and then push the updates back to github which actually updates the site.    The same process (slightly modified) is used for the documentation except the sources for the docs live in the llvmpy repo under the docs directory and the built pages are pushed to the gh-pages branch of the llvmpy-doc repo.    If you are editing sources you only modify llvmpy/docs and llvmpy-webpage files.   The other repos are generated and pushed via scripts.

We are using the same general scheme to host the numba pages (although there I couldn't get the numba.org domain name and so I am using http://numba.pydata.org).   With llvmpy on a relatively solid footing, attention could be shifted to getting a Numba release out.  Today, we finally released Numba 0.1.   It took longer than expected after the SciPy conference mainly because we were hoping that some of the changes (still currently in a devel branch) to use an AST-based code-generator could be merged into the main-line before the release.  

Jon Riehl did the lion's share of the work to transform Numba from my early prototype to a functioning system in 0.1 with funding from Continuum Analytics, Inc.   Thanks to him, I can proudly say that Numba is ready to be tried and used.    It is still early software --- but it is ready for wider testing.   One of the problems you will have with Numba right now is error reporting.  If you make a mistake in the Python code that you are decorating, the error you get will not be informative -- so test the Python code before decorating it with Numba.    But, if you get things right, Numba can speed up your Python code by 200 times or more.    It's is really pretty fun to be able to write image-processing routines in Python.   PyPy can do this too, of course, but with Numba you have full integration with the CPython stack and you don't have to wait for someone to port the library you also want to use to PyPy.

Numba's road-map is being defined right now by the people involved in the project.  On the horizon is support for NumPy index expressions (slices, etc.), merging of the devel branch which uses the AST and Mark Florisson's minivect compiler, improving support for error checking, emitting calls to the Python C-API for code that cannot be type-specialized, and improving complex-number support.  Your suggestions are welcome.

More PyPy discussions

I'm very glad that my co-founder of Continuum Analytics,  Peter Wang, has published his recent follow-up blog-post that hopefully clarifies his perspective on the on-going dialogue about CPython and PyPy.

Peter is a fundamentally good-natured person, and he is a lot of fun to be around --- even when he is disagreeing with you.   I'm very fortunate to be working with him on a daily basis.   He can be opinionated, but his ability to connect deeply to a wide-variety of subjects means that you come away from a dialogue with him having learned something (even if you still remain unconvinced by his views).  

Peter is also one of the smartest people I've ever met.   One of my great memories in life is sitting at dinner with Peter and Eric Weinstein while those two great minds treated me, Wes McKinney, and Adam Klein to the most impressive display of metaphor ping-pong I've ever seen covering a wide-variety of topics from social justice to string theory.  I could keep up with the dialogue, but not enough to really participate meaningfully --- and the other two Ivy-league-educated dinner partners were in the same boat.

I fundamentally agree with Peter's perspective that CPython-the-runtime is and will remain the centerpiece of the Python conversation.    In fact, I would say that even more focus needs to be on CPython-the-runtime.   It is great to see improvements in Python 3.3 like the completion of the memory-view implementation and the fixing of the internal string (Unicode) representation, but there are many other improvements that could be made.

It is a wonderful and inspiring thing to see great developers think out of the box with novel projects like Jython, IronPython, and PyPy.   Nonetheless from my perspective we still have a long way to go to really connect the average developer with ideas of array-oriented computing that could really help the continuing onslaught of parallel-devices-in-search-of-software.   As a result, it feels like those wanting Java, .NET, and machine-code integration would be better served by more attention on JPype, Python.NET, LLVMPy, and even CorePy.   Such efforts would also be better for the entire user-base of Python --- especially a majority of industry uses of Python.

But regardless of my perspective, I'm encouraged by the PyPy developer enthusiasm, and I do want to encourage dialogue regardless of my views.   As a result, I am very happy to report that both NumFOCUS and Continuum Analytics recently joined forces to sponsor Maciej Fijalkowski on a small project to create an embedded version of PyPy --- a "PyPy-in-a-Box."  This is an integration of PyPy to the CPython run-time (so that you can speed-up a particular CPython function by calling out to a library-version of PyPy).   This is proof-of-concept code so it is not appropriate for production --- but it is a good example of what is possible when we all work together to promote the Python ecosystem.

The online project is here:  https://bitbucket.org/fijal/hack2/src/default/pypyembed  and you can get a binary version that works on 64-bit Linux here:  http://baroquesoftware.com/~fijal/pypy-1.9-in-a-box-linux64.tar.bz2.

This approach needs more development to be a viable tool in the CPython ecosystem, but one of my suggestions to the PyPy community is that they focus on "shedding-tools" like this one for the CPython world --- so that everyone can benefit from their innovations.   With an integration effort like embeded PyPy, one can also make better comparisons with tools like Numba --- another dynamic-compilation run-time that uses LLVM and LLVM-py.     Numba has made a lot of progress in the last few months.   In fact, I recently gave a talk on the project at the well-attended SciPy2012 conference in Austin.   You can view my slides that outline and motivate the project online.   An actual release of the project is imminent, but you can already use Numba to very easily write signficant Python code using NumPy arrays that executes at "C-speeds."  But, that is worth another blog-post of its own....

Transition to Continuum

Our lives are punctuated by transformational events:  the birth of a child, finishing school, the passing of a loved one, meeting someone special.   Even without the regular beating of celestial rhythms to provide opportunities for renewal we would have these moments to measure our lives by.   Once in a while, the rhythmic and asynchronous coincide providing a particularly poignant opportunity for change.   Jan 1, 2012 was just such a time for me as I left my position as President of Enthought to start a new venture with Peter Wang (author of Chaco) and others.    Our new company is Continuum Analytics, Inc. (or just Continuum).  Our nascent website initially targeted only to the Python initiate is http://www.continuum.io.  

While I am ecstatic about the new venture, I will definitely miss the team that we've built around the world that has delivered Enthought's second consecutive record year.   This team of exceptional individuals has been very successful at improving and expanding the Python story in a few targeted companies inside of the Fortune 50 as well as making it easy to install Python for the masses.    Those who have taken the time to first install and then learn Traits, TraitsUI, Chaco, MayaVi, and the rest of the Enthought Tool Suite have had their efforts rewarded with increased productivity in the creation of rich client UIs and improved pluggable, scriptable, and component-based architecture.     It has been a highly educational experience to participate with Enthought.  There is much you learn about business, people, and the world when a software consulting company grows from 1 office with fewer than 17 people to 4 offices around the world and nearly 50 people.   I will always be grateful to the Enthought founders, employees (past and present), and customers (past and present) for the relationships, the trust, and the thoughtful times we shared in learning, growing, and serving each other.

My heart, however, has always been and continues to be with NumPy and SciPy which need more support than Enthought can currently provide --- so I must move on.   It took a lot of trust from my wife when (with 3 small children at home) she patiently waited for me while I spent all of 1999 writing Multipack (which in 2001 formed the bulk of SciPy).  It also took trust when in 2005 (with now 5 children at home) she watched me sacrifice my tenure-track position by writing NumPy instead of more papers.   In 2012 (with now 6 children at home), I'm asking her to trust me one more time while I leave a comfortable salary with a good company to put more effort full time in helping take NumPy and SciPy to the next level.

Over the past 4 1/2 years consulting with large companies I have learned a great deal about what NumPy (and SciPy) can and should be.  These and related tools in the Python ecosystem need to become significant pieces to real solutions to the data analytics challenges that face us.   R, Hadoop, and other (proprietary) solutions are already staking their claim on the space that Python should be dominating.   Python has significant traction in science and analysis but too little publicity in the nascent nomenclature of data analytics.    In order to accelerate the processing capabilities of Python and related tools, much progress needs to be made.   My New Year's resolution this year is to begin to contribute more substantially to that progress by organizing a new company that will hopefully allow many people to spend significant time directly on NumPy and SciPy during working hours.  I also hope to assist any public, non-profit efforts towards that mutual goal as well.   I also hope to be able to spend more time myself on NumPy and SciPy.

To realize my hopes long term, the company must succeed.   For the company to succeed it must find customers --- people willing to buy something that it sells.   People are appropriately particular about what they buy.  Making products that delight will require a lot of work from Continuum, but I am excited to help organize and work alongside the best team we can put together to do it.   This may also mean different business models and licensing around some of the NumPy-related code that the company writes.   I recognize this may cause some raised eye-brows.   I deeply value making code freely available.    I'm a Jeffersonian at heart and believe that ideas (including code) should be shared freely.   Six years ago I experimented by selling my "Guide to NumPy" long enough to make sufficient money to justify the effort.   The book ended up in the public domain and contributed substantially to the current NumPy documentation.   This is an illustration of how resources can be allocated to full-time attention and then later made available for all to enjoy.   Of course there are other models that also work to accomplish similar ends and we will be actively exploring a few of them.

Despite my ideals, my wife thanks me that I'm a pragmatist with children to provide for.   In addition, I have watched wearily as it's been difficult to find volunteer labor (including my own) to turn NumPy into the data-management and data-analytics substrate that it should already be.   All of this happens while huge sums of money are wasted at companies large and small inefficiently transforming raw, but inaccessible data into something closer to information that can be used for decisions by the domain experts.    The information available is not what it can be.  The amount of effort it takes to transform the data to actionable information is not where it can be.  The wide-spread understanding about how to program parallel and distributed machines is not where it can be.  We can and must do better in figuring out how to get full-time attention on NumPy and related tools while still making them widely available.

At Continuum, we have a vision for significantly changing how people manipulate, transform, and uncover their data.   We also have customer-driven plans to achieve it, and we are going to put our full energy into it.   So far, the development team consists of Peter Wang, me, Mark Wiebe, Francesc Alted (PyTables), and Bryan Van de Ven.   We will also be getting part-time but important development help from Hugo Shi and Andy Terrel.  In addition, we are building an initial support/business staff to help us build and grow the business.     We plan to continue to collaborate with others in the community both commercial (e.g. Wes McKinney in his new startup: Lambda Foundry) and open (e.g. Fernando Perez, Brian Granger, Min Ragan-Kelley of IPython fame). If you are interested in either joining us or collaborating with us, please send us an email at info@continuum.io.  Also, please follow us on Twitter @ContinuumIO or Like us on Facebook.  

We are actively looking for customer partners, as well.  If you are interested in learning more about where we are heading and how that might help you, please drop us a line, or come see us at PyCon this year.   We will also be at Strata, and afterwords we will be hosting a Python Data Workshop ("PyData") at the Googleplex.   Please sign up for the PyData workshop wait-list at http://pydataworkshop.eventbrite.com/ (we could only find room for 50 people at the Googleplex).  However, given that the event is free of charge, I'm expecting some people who have reserved their spot may not actually be able to attend.  So, signing up on the wait-list is still worthwhile.

This year will be an exciting one for us.  When I get a spare moment, I still hope to finish a few of the blogs that I've started and possibly include some more that describe more of what I've learned over the past several years as a scientist/engineer-turned-software developer, lessons about running a software company, more of where we are headed at Continuum, reflections on open source, and other more technical ramblings.

Thoughts on porting NumPy to PyPy

Last weekend, I attended GitHub's PyCodeConf in Miami Florida and had the opportunity to give a talk on array-oriented computing and Python.  I would like to thank Tom and Chris (GitHub founders) for allowing me to come speak.  I enjoyed my time there, but I have to admit I felt old and a bit out of place.   There were a lot of young people there who understood a lot more about making web-pages and working at cool web start-ups than solving partial differential equations with arrays.  Fortunately, Dave Beazley and Raymond Hettinger were there so I didn't feel completely ancient.   In addition, Wes McKinney and Peter Wang helped me represent the NumPy community.

At the conference I was reminded of PyPy's recent success at showing the speedups some said weren't possible from a dynamic language like Python --- speed-ups which make it possible to achieve C-like speeds using Python constructs.   This is very nice as it illustrates again that high-level languages can be compiled to low-level speeds.    I also became keenly aware of the enthusiasm that has cropped up in porting NumPy to PyPy.   I am happy for this enthusiasm as it illustrates the popularity of NumPy which pleases me.   On the other hand, in every discussion I have heard or read about this effort, I'm not convinced that anyone excited about this porting effort actually understands the complexity of what they are trying to do nor the dangers that it could create in splitting the small community of developers who are regularly contributing to NumPy and SciPy and causing confusion for the user base of Python in Science.

I'm hopeful that I can provide some perspective.   Before I do this, however, I want to congratulate the PyPy team and emphasize that I have the utmost respect for the PyPy developers and what they have achieved.   I am also a true-believer in the ability for high-level languages to achieve faster-than-C speeds.   In fact, I'm not satisfied with a Python JIT.  I want the NumPy constructs such as vectorization, fancy indexing, and reduction to be JIT compiled.    I also think that there are use-cases of NumPy all by itself that makes it somewhat interesting to do a NumPy-in-PyPy port.    I would also welcome the potential things that can be learned about how to improve NumPy that would come out of trying to write a version of NumPy in RPython.

However, to avoid detracting from the overall success of Python in Science, Statistics, and Data Analysis,  I think it is important that 3 things are completely clear to people interested in the NumPy on PyPy idea. 

1. NumPy is just the beginning (SciPy, matplotlib, scikits, and 100s of other packages and legacy C/C++ and Fortran code are all very important)
2. NumPy should be a lot faster than it is currently.
3. NumPy has an ambitious roadmap and will be moving forward rather quickly over the coming years.  

NumPy is just the beginning 

Most of the people who use NumPy use it as an entry-point to the entire ecosystem of Scientific Packages available for Python.    This ecosystem is huge.   There are at least 1 million unique visitors to the http://www.scipy.org site every year and that is just an entry point to the very large and diverse community of technical computing users who rely on Python. 

Most of the scientists and engineers who have come to Python over the past years have done so because it is so easy to integrate their legacy C/C++ and Fortran code into Python.   National laboratories, large oil companies, large banks and many other Fortune 50 companies all must integrate their code into Python in order for Python to be part of their story.   NumPy is part of the answer  that helps them seamlessly view large amounts of data as arrays in Python or as arrays in another compiled language without the non-starter of copying the memory back and forth. 

Once the port of NumPy to PyPy has finished, are you going to port SciPy?  Are you going to port matplotlib?  Are you going to port scikits.learn, or scikits.statsmodels?   What about Sage?   Most of these rely on not just the Python C-API but also the NumPy C-API which you would have to have a story for to make a serious technical user of Python get excited about a NumPy port to PyPy.  

To me it is much easier to think about taking the ideas of PyPy and pulling them into the Scientific Python ecosystem then going the other way around.    It's not to say that there isn't some value in re-writing NumPy in PyPy, it just shouldn't be over-sold and those who fund it should understand what they aren't getting in the transaction.  

C-speed is the wrong target

Several examples including my own previous blog post has shown that vectorized Fortran 90 can be 4-10 times faster than NumPy.   Thus, we know there is room for improvement even on current single-core machines.   This doesn't even take into account the optimizations that should be possible for multiple-cores, GPUs and even FPGAs all of which are in use today but are not being utilized to the degree they should be.   NumPy needs to adapt to make use of this kind of hardware and will adapt in time.

NumPy will be evolving rapidly over the coming years

The pace of NumPy development has leveled off in recent years, but this year has created a firestorm of new ideas that will be coming to fruition over the next 1-2 years and NumPy will be evolving fairly rapidly during that time.    I am committed to making this happen and will be working very hard in 2012 on the code-base itself to realize some of the ideas that have emerged.   Some of this work will require some re-factoring and re-writing as well.   I would honestly rather collaborate with PyPy than compete, but my constraints are that I care very much about backward compatibility and very much about the entire SciPy ecosystem.    I sacrificed a year of my life in 1999 (delaying my PhD graduation by at least 6-12 months) bringing SciPy to life.    I sacrificed my tenure-track position in academia bringing NumPy to life in 2005.   Constraints of keeping my family fed, clothed, and housed seem to keep me on this 6-7 year sabbatical-like cycle for SciPy/NumPy but it looks like next year I will finally be in a position to spend substantial time and take the next steps with NumPy to help it progress to the next stage.

Some of the ideas that will be implemented include:

• integration of non-contiguous memory chunks into the NumPy array structure (generalization of strides)
• addition of labels to axes and dimensions (generalization of shape)
• derived fields, enumerated data-types, reference data-types, and indices for structured arrays
• improvements to the data-type infrastructure to make it easier to add new data-types
• improvements to the calculation infrastructure (iterators and fast general looping constructs)
• fancy-indexing as views
• integration of Pandas group-by features
• missing data bit-patterns
• distributed arrays 

Over conversations with many people this year, more ideas than room to talk about them have emerged and I am excited to start seeing these ideas come to fruition to make NumPy and Python the best solution for data-analysis.    Beginning next year, I will be pushing hard for their introduction into the NumPy/SciPy ecosystem --- with a careful eye on backward compatibility which has long been one of NumPy's strengths.

A way forward

I would love to see more scientific code written at a high-level without sacrificing run-time performance.   The high-level intent allows for the creation of faster machine code than lower-level translations of the intent often does. I know this is possible and I intend to do everything I can professionally to see that happen (but from within the context of the entire SciPy ecosystem). As this work emerges, I will encourage PyPy developers to join us using the hard-won knowledge and tools they have created.

Even if PyPy continues as a separate eco-system, then there are points of collaboration that will benefit both groups. One of these is to continue the effort Microsoft initially funded to separate the C-parts of NumPy away from the CPython interface to NumPy. This work is now in a separate branch that has diverged from the main NumPy branch and needs to be re-visited. If people interested in NumPy on PyPy spent time on improving this refactoring into basically a NumPy C-library, then PyPy can call this independent library using its methods for making native calls just as CPython can call it using its extension approach. Then IronPython, Jython (and for that matter Ruby, Javascript, etc.) can all call the C-library and leverage the code. There is some effort to do this and it's not trivial. Perhaps, there is even a way for PyPy to generate C-libraries from Python source code --- now that would be an interesting way to collaborate.

The second way forward is for PyPy to interact better with the Cython community. Support in PyPy for Cython extension modules would be a first step. There is wide agreement among NumPy developers that more of NumPy should be written at a high-level (probably using Cython). Cython already is used to implement many, many extension modules for the Sage project. William Stein's valiant efforts in that community have made Cython the de-facto standard for how most scientists and engineers are writing extensions modules for Python these days. This is a good thing for efforts like PyPy because it adds a layer of indirection that allows PyPy to make a Cython backend and avoid the Python C-API.  

I was quite astonished that Cython never came up in the panel discussion at the last PyCon when representatives from CPython, PyPy, IronPython, and Jython all talked about the Python VMs.   To me that oversite was very troublesome.   I was left doubting the PyPy community after Cython was not mentioned at all --- even when the discussion about how to manage extensions to the language came up during the panel discussion.    It shows that pure Python developers on all fronts have lost sight of what the scientific Python community is doing.   This is dangerous. I encourage Python developers to come to a SciPy conference and take a peek at what is going on.   I hope to be able to contribute more to the discussion as well. 

If you are a Python developer and want to extend an olive leaf, then put a matrix infix operator into the language.  It's way past time :-)

Speeding up Python Again

After getting a few great comments on my recent post --- especially regarding using PyPy and Fortran90 to speed up Python --- I decided my simple comparison needed an update.

The big-news is that my tests for this problem actually showed PyPy quite favorably (even a bit faster than the CPython NumPy solution). This is very interesting indeed! I knew PyPy was improving, but this shows it has really come a long way.

Also, I updated the Python-only comparison to not use NumPy arrays at all. It is well-known that NumPy arrays are not very efficient containers for doing element-by-element calculations in Python syntax. There is both more overhead for getting and setting elements than there is for simple lists, and the NumPy scalars that are returned when specific elements of NumPy arrays are selected can be a bit slow when doing scalar math computations on the Python side.

Finally, I included a Fortran 90 example based on the code and comments provided by SymPy author Ondrej Certik. Fortran 77 was part of the original comparison that Prabhu Ramanchandran put together several years ago. Fortran 90 includes some nice constructs for vectorization that make it's update code very similar to the NumPy update solution. Apparently, gfortran can optimize this kind of code very well. In fact, the Fortran 90 solution was the very best of all of the approaches I took (about 4x faster than the NumPy solution and 2x faster than the other compiled approaches).

At Prabhu's suggestion, I made the code available at github under a new GitHub repository in the SciPy project so that others could contribute and provide additional comparisons.

The new results are summarized in the following table which I updated to running on a 150x150 grid with again 8000 iterations.

Method	Time (sec)	Relative Speed
Pure Python	202	36.3
NumExpr	8.04	1.45
NumPy	5.56	1
PyPy	4.71	0.85
Weave	2.42	0.44
Cython	2.21	0.40
Looped Fortran	2.19	0.39
Vectorized Fortran	1.42	0.26

The code for both the Pure Python and the PyPy solution is laplace2.py. This code uses a list-of-lists for the storage of the values. The same code produces the Pure Python solution and the PyPy solution. The only difference is that one is run with the standard CPython and the other with the PyPy binary. Here is sys.version from the PyPy binary used to obtain these results:

'2.7.1 (b590cf6de419, Apr 30 2011, 03:30:00)\n[PyPy 1.5.0-alpha0 with GCC 4.0.1]'

This is a pretty impressive achievement for the PyPy team. Kudos!

For the other solutions, the code that was executed is located at laplace.py. The Fortran 90 module compiled and made available to Python with f2py is located at _laplace.f90. The single Cython solution is located at _laplace.pyx.

It may be of interest to some to see what the actual calculated potential field looks like. Here is an image of the 150x150 grid after 8000 iterations:

Here is a plot showing three lines from the image (at columns 30, 80, 130 respectively):

It would be interesting to add more results (from IronPython, Jython, pure C++, etc.). Feel free to check out the code from github and experiment. Alternatively, add additional problems to the speed project on SciPy and make more comparisons. It is clear that you can get squeeze that last ounce of speed out of Python by linking to machine code. It also seems clear that there is enough information in the vectorized NumPy expression to be able to produce fast machine code automatically --- even faster than is possible with an explicit loop. The PyPy project shows that generally-available JIT-technology for Python is here and the scientific computing community should grapple with how we will make use of it (and improve upon it). My prediction is that we can look forward to more of that in the coming months and years.

Speeding up Python (NumPy, Cython, and Weave)

The high-level nature of Python makes it very easy to program, read, and reason about code. Many programmers report being more productive in Python. For example, Robert Kern once told me that "Python gets out of my way" when I asked him why he likes Python. Others express it as "Python fits your brain." My experience resonates with both of these comments.

It is not rare, however, to need to do many calculations over a lot of data. No matter how fast computers get, there will always be cases where you still need the code to be as fast as you can get it. In those cases, I first reach for NumPy which provides high-level expressions of fast low-level calculations over large arrays. With NumPy's rich slicing and broadcasting capabilities, as well as its full suite of vectorized calculation routines, I can quite often do the number crunching I am trying to do with very little effort.

Even with NumPy's fast vectorized calculations, however, there are still times when either the vectorization is too complex, or it uses too much memory. It is also sometimes just easier to express the calculation with a simple loop. For those parts of the application, there are two general approaches that work really well to get you back to compiled speeds: weave or Cython.

Weave is a sub-package of SciPy and allows you to inline arbitrary C or C++ code into an extension module that is dynamically loaded into Python and executed in-line with the rest of your Python code. The code is compiled and linked at run-time the very first time the code is executed. The compiled code is then cached on-disk and made available for immediate later use if it is called again.

Cython is an extension-module generator for Python that allows you to write Python-looking code (Python syntax with type declarations) that is then pre-compiled to an extension module for later dynamic linking into the Python run-time. Cython translates Python-looking code into "not-for-human-eyes" C-code that compiles to reasonably fast C-code. Cython has been gaining a lot of momentum in recent years as people who have never learned C, can use Cython to get C-speeds exactly where they want them starting from working Python code. Even though I feel quite comfortable in C, my appreciation for Cython has been growing over the past few years, and I know am an avid supporter of the Cython community and like to help it whenever I can.

Recently I re-did the same example that Prabhu Ramachandran first created several years ago which is reported here. This example solves Laplace's equation over a 2-d rectangular grid using a simple iterative method. The code finds a two-dimensional function, u, where ∇² u = 0, given some fixed boundary conditions.

Pure Python Solution

The pure Python solution is the following:

from numpy import zeros
from scipy import weave

dx = 0.1
dy = 0.1
dx2 = dx*dx
dy2 = dy*dy

def py_update(u):
    nx, ny = u.shape
    for i in xrange(1,nx-1):
        for j in xrange(1, ny-1):
            u[i,j] = ((u[i+1, j] + u[i-1, j]) * dy2 +
                      (u[i, j+1] + u[i, j-1]) * dx2) / (2*(dx2+dy2))

def calc(N, Niter=100, func=py_update, args=()):
    u = zeros([N, N])
    u[0] = 1
    for i in range(Niter):
        func(u,*args)
    return u

This code takes a very long time to run in order to converge to the correct solution. For a 100x100 grid, visually-indistinguishable convergence occurs after about 8000 iterations. The pure Python solution took an estimated 560 seconds (9 minutes) to finish (using IPython's %timeit magic command).

NumPy Solution

Using NumPy, we can speed this code up significantly by using slicing and vectorized (automatic looping) calculations that replace the explicit loops in the Python-only solution. The NumPy update code is:

def num_update(u):
    u[1:-1,1:-1] = ((u[2:,1:-1]+u[:-2,1:-1])*dy2 + 
                    (u[1:-1,2:] + u[1:-1,:-2])*dx2) / (2*(dx2+dy2))

Using num_update as the calculation function reduced the time for 8000 iterations on a 100x100 grid to only 2.24 seconds (a 250x speed-up).   Such speed-ups are not uncommon when using NumPy to replace Python loops where the inner loop is doing simple math on basic data-types.

Quite often it is sufficient to stop there and move on to another part of the code-base.  Even though you might be able to speed up this section of code more, it may not be the critical path anymore in your over-all problem.  Programmer effort should be spent where more benefit will be obtained.  Occasionally, however, it is essential to speed-up even this kind of code.

Even though NumPy implements the calculations at compiled speeds, it is possible to get even faster code.   This is mostly because NumPy needs to create temporary arrays to hold intermediate simple calculations in expressions like the average of adjacent cells shown above.  If you were to implement such a calculation in C/C++ or Fortran, you would likely create a single loop with no intermediate temporary memory allocations and perform a more complex computation at each iteration of the loop.

In order to get an optimized version of the update function, we need a machine-code implementation  that Python can call.   Of course, we could do this manually by writing the inner call in a compilable language and using Python's extension facilities.  More simply, we can use Cython and Weave which do most of the heavy lifting for us.

Cython solution

Cython is an extension-module writing language that looks a lot like Python except for optional type declarations for variables.  These type declarations allow the Cython compiler to replace generic, highly dynamic Python code with specific and very fast compiled code that is then able to be loaded into the Python run-time dynamically.  Here is the Cython code for the update function:

cimport numpy as np

def cy_update(np.ndarray[double, ndim=2] u, double dx2, double dy2):
    cdef unsigned int i, j
    for i in xrange(1,u.shape[0]-1):
        for j in xrange(1, u.shape[1]-1):
            u[i,j] = ((u[i+1, j] + u[i-1, j]) * dy2 +
                      (u[i, j+1] + u[i, j-1]) * dx2) / (2*(dx2+dy2))

This code looks very similar to the original Python-only implementation except for the additional type-declarations.   Notice that even NumPy arrays can be declared with Cython and Cython will correctly translate Python element selection into fast memory-access macros in the generated C code.   When this function was used for each iteration in the inner calculation loop,  the 8000 iterations on a 100x100 grid took only 1.28 seconds.

For completeness, the following shows the contents of the setup.py file that was also created in order to produce a compiled-module where the cy_update function lived.

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

import numpy

ext = Extension("laplace", ["laplace.pyx"],
    include_dirs = [numpy.get_include()])
                
setup(ext_modules=[ext],
      cmdclass = {'build_ext': build_ext})

The extension module was then built using the command: python setup.py build_ext --inplace

Weave solution

An older, but still useful, approach to speeding up code is to use weave to directly embed a C or C++ implementation of the algorithm into the Python program directly.   Weave is a module that surrounds the bit of C or C++ code that you write with a template to on-the-fly create an extension module that is compiled and then dynamically loaded into the Python run-time.   Weave has a caching mechanism so that different strings or different types of inputs lead to a new extension module being created, compiled, and loaded.    The first time code using weave runs, the compilation has to take place.  Subsequent runs of the same code will load the cached extension module and run the machine code.

For this particular case, an update routine using weave looks like:

def weave_update(u):
    code = """
    int i, j;
    for (i=1; i<Nu[0]-1; i++) {
       for (j=1; j<Nu[1]-1; j++) {
           U2(i,j) = ((U2(i+1, j) + U2(i-1, j))*dy2 + \
                       (U2(i, j+1) + U2(i, j-1))*dx2) / (2*(dx2+dy2));
       }
    }
    """
    weave.inline(code, ['u', 'dx2', 'dy2'])

The inline function takes a string of C or C++ code plus a list of variable names that will be pushed from the Python namespace into the compiled code.   The inline function takes this code and the list of variables and either loads and executes a function in a previously-created extension module (if the string and types of the variables have been previously created) or else creates a new extension module before compiling, loading, and executing the code.

Notice that weave defines special macros so that U2 allows referencing the elements of the 2-d array u using simple expressions.   Weave also defines the special C-array of integers Nu to contain the shape of the u array.   There are also special macros similarly defined to access the elements of array u if it would have been a 1-, 3-, or 4-dimensional array (U1, U3, and U4).   Although not used in this snippet of code, the C-array Su containing the strides in each dimension and the integer Du defining the number of dimensions of the array are both also defined.

Using the weave_update function, 8000 iterations on a 100x100 grid took only 1.02 seconds.   This was the fastest implementation of all of the methods used.    Knowing a little C and having a compiler on hand can certainly speed up critical sections of code in a big way.

Faster Cython solution (Update)

After I originally published this post, I received some great feedback in the Comments section that encouraged me to add some parameters to the Cython solution in order to get an even faster solution. I was also reminded about pyximport and given example code to make it work more easily. Basically by adding some compiler directives to Cython to avoid some checks at each iteration of the loop, Cython generated even faster C-code. To the top of my previous Cython code, I added a few lines:

#cython: boundscheck=False
#cython: wraparound=False

I then saved this new file as _laplace.pyx, and added the following lines to the top of the Python file that was running the examples:

import pyximport
import numpy as np
pyximport.install(setup_args={'include_dirs':[np.get_include()]})
from _laplace import cy_update as cy_update2

This provided an update function cy_update2 that resulted in the very fastest implementation (943 ms) for 8000 iterations of a 100x100 grid.

Summary

The following table summarizes the results which were all obtained on a 2.66 Ghz Intel Core i7 MacBook Pro with 8GB of 1067 Mhz DDR3 Memory. The relative speed column shows the speed relative to the NumPy implementation.

Method	Time (sec)	Relative Speed
Pure Python	560	250
NumPy	2.24	1
Cython	1.28	0.57
Weave	1.02	0.45
Faster Cython	0.94	0.42

Clearly when it comes to doing a lot of heavy number crunching, Pure Python is not really an option. However, perhaps somewhat surprisingly, NumPy can get you most of the way to compiled speeds through vectorization. In situations where you still need the last ounce of speed in a critical section, or when it either requires a PhD in NumPy-ology to vectorize the solution or it results in too much memory overhead, you can reach for Cython or Weave. If you already know C/C++, then weave is a simple and speedy solution. If, however, you are not already familiar with C then you may find Cython to be exactly what you are looking for to get the speed you need out of Python.

Python Enhancement Proposals I Wish I Had Time to Champion

Today I was trying to make progress on a few different NumPy proposal enhancements and ended up frustrated knowing that come Monday morning, I will not have any time to follow-up on them. Managing a growing consulting company takes a lot of time (Enthought is over 30 people now and growing to near 50 by the end of the year). There are countless meetings devoted to new hires, program development, project reviews, customer relations, budgeting, and sales. I also take a direct role in delivering on training and select consulting projects. Someday I may get a chance to write something of use about things I've learned along the way, but that is for another day (and likely another blog). This post is to get a few ideas I've been sitting on written down in the hopes that somebody might read them and get excited about contributing. At the very least anybody that reads this post, will know (at least some of) my current opinion about a few technical proposals.

About a month ago, I had the privilege of organizing a "data-array" summit in which several people in the NumPy and SciPy community came together at the Enthought offices to discuss some ideas related to how to improve data analysis with the NumPy and SciPy stack. We spent 3-days thinking and brainstorming which led to many fruitful discussions. I expect that some of the ideas generated will result in important and interesting changes to NumPy and SciPy over the coming months and years. More information about the summit can be learned by listening to the relevant inSCIght podcast.

It's actually a very exciting time to get involved in the SciPy community as Python takes its place as one of the approaches people will be using to analyze all the data that we are generating. In that spirit, I wanted to express a few of what I consider to be important enhancements that are needed to Python and NumPy.

I will start with Python and leave NumPy to another post. Here there are really three big missing features that would really benefit those of us who use Python for technical computing. Unfortunately, I don't think there is enough representation of the Python for Science crowd in the current batch of Python developers. This is not due to any exclusion from the Python developers who have always been very accommodating. It is simply due to the scarcity of people who understand the SciPy perspective and use-cases also willing to engage with developers in the Python world. Those (like Mark Dickinson) who cross the chasm are a real gem.

If anyone has an interest in shepherding a PEP in any of the following directions, you will have my '+1' support (and any community-organizing that I can muster to help you).  Honestly, if these things were put into Python 3, there would be a serious motivation to move to Python 3 for the scientific community (which is otherwise going to lag in the great migration).

Python Enhancements I Want

Adding additional operators

We need additional operators to easily represent at least matrix multiplication, matrix power, and matrix solve). I could possibly back-off the last two if we at least had matrix multiplication. This should have been done a long time ago. If I had been able to spare the time, I would have pushed to hold off porting of NumPy to Python 3 until we got matrix multiplication operators. Yes, I know that blackmail usually backfires and thankfully Pauli Virtanen and Charles Harris acted before I even had a chance to suggest such a thing :-). But, seriously, we need this support in the language. 

The reasons are quite simple:

• Syntax matters: writing d = numpy.solve(numpy.dot(numpy.dot(a,b),c), x) is a whole lot more ugly than something like d = (a*b*c) \ x.  If the former is fine, then we should all just go back to writing LISP.  The point of having nice syntax is to minimize the line-noise and mental overhead of mapping the mental idea to working code.   For Python to be used with mental efficiency in technical computing you need to write expressions involving higher-order operations like this all the time.   
• Right now, the recommended way to do this is to convert a, b, c, and x to "matrices", perform the computation in a nice expression and then convert back to arrays.  This is clunky at best.

I've been back and forth on this for 13 years and can definitively say that we would be much better off in Python if we had a matrix multiplication operator.   Please, please, can we get one!    The relevant PEPS where this has been discussed are:  PEP 211 and PEP 225.  I think I like having more than just one operator added (ala PEP 225, but the subject would have to be re-visited by a brave soul).

Overloadable Boolean Operations

PEP 335 was a fantastic idea. I really wish we had the ability to overload and, or, and not. Among other things, this would allow the very nice syntax so that mask = 2<a<10 would generate an array of True and False values when a is an array. Currently, to generate this same mask you have to do (2<a)&(a<4). The PEP has other important use-cases as well. It would be excellent if this PEP were re-visited, championed, and put into Python 3.

Allowing slice object literals outside of []

Python's syntax allows construction of a slice object inside brackets so that one can write a[1:3] which is equivalent to a.__getitem__(slice(1,3)). Many times over the years, I have wanted to be able to specify a slice object using the syntax start:stop:step outside of the getitem. Even, if Python's parser were extended to allow the slice literal to be accepted as the input to a function it would be preferred. The biggest wart this would remove is the (ab)use of getitem to return new ranges and grids in NumPy (go use mgrid and r_ in NumPy to see what I mean). I would prefer that these were functions, but I would need mgrid(1:5, 1:5) to work.

There was a PEP for range literals (PEP 204) once upon a time.  There were some interesting aspects about that proposal, but frankly I don't want the slice syntax to produce ranges.  I would just be content for it always to produce slice objects --- just allow it outside of brackets.

While I started by lamenting my lack of time to implement NumPy enhancements, I will leave the discussion of NumPy enhancements I'm dreaming about to another post.   I would be thrilled if somebody took up the charge to push any of these Python enhancements in Python 3.   If Python 3 ends up with any of them, it would be a huge motivation to me to migrate to Python 3 entirely.