Showing posts with label numpy. Show all posts
Showing posts with label numpy. Show all posts

Tuesday, April 15, 2014

NumPy on PyPy - Status Update

Work on NumPy on PyPy continued in March, though at a lighter pace than the previous few months. Progress was made on both compatibility and speed fronts. Several behavioral issues reported to the bug tracker were resolved. The most significant of these was probably the correction of casting to built-in Python types. Previously, int/long conversions of numpy scalars such as inf/nan/1e100 would return bogus results. Now, they raise or return values, as appropriate.

On the speed front, enhancements to the PyPy JIT were made to support virtualizing the raw_store/raw_load memory operations used in numpy arrays. Further work remains here in virtualizing the alloc_raw_storage when possible. This will allow scalars to have storages but still be virtualized when possible in loops.

Aside from continued work on compatibility/speed of existing code, we also hope to begin implementing the C-level components of other numpy modules such as mtrand, nditer, linalg, and so on. Several approaches could be taken to get C-level code in these modules working, ranging from reimplementing in RPython to interfacing with existing code with CFFI, if possible. The appropriate approach depends on many factors and will probably vary from module to module.

To try out PyPy + NumPy, grab a nightly PyPy and install our NumPy fork. Feel free to report comments/issues to IRC, our mailing list, or bug tracker. Thanks to the contributors to the NumPy on PyPy proposal for supporting this work.

Friday, March 7, 2014

NumPy on PyPy - Progress in February

More progress was made on the NumPy front in the past month. On the compatibility front, we now pass ~130 more tests from NumPy's suite since the end of January. Currently, we pass 2336 tests out of 3265 tests run, with many of the failures representing portions of NumPy that we don't plan to implement in the near future (object dtypes, unicode, etc). There are still some failures that do represent issues, such as special indexing cases and failures to respect subclassed ndarrays in return values, which we do plan to resolve. There are also some unimplemented components and ufuncs remaining which we hope to implement, such as nditer and mtrand. Overall, the most common array functionality should be working.

Additionally, I began to take a look at some of the loops generated by our code. One widely used loop is dot, and we were running about 5x slower than NumPy's C version. I was able to optimize the dot loop and also the general array iterator to get us to ~1.5x NumPy C time on dot operations of various sizes. Further progress in this area could be made by using CFFI to tie into BLAS libraries, when available. Also, work remains in examining traces generated for our other loops and checking for potential optimizations.

To try out PyPy + NumPy, grab a nightly PyPy and install our NumPy fork. Feel free to report comments/issues to IRC, our mailing list, or bug tracker. Thanks to the contributors to the NumPy on PyPy proposal for supporting this work.

Cheers,
Brian

Thursday, February 6, 2014

NumPy Status Update - December/January

Work continued on the NumPy + PyPy front steadily in December and more lightly in January. The continued focus was compatibility, targeting incorrect or unimplemented features that appeared in multiple NumPy test suite failures. We now pass ~2/3 of the NumPy test suite. The biggest improvements were made in these areas:

- Bugs in conversions of arrays/scalars to/from native types
- Fix cases where we would choose incorrect dtypes when initializing or computing results
- Improve handling of subclasses of ndarray through computations
- Support some optional arguments for array methods that are used in the pure-python part of NumPy
- Support additional attributes in arrays, array.flags, and dtypes
- Fix some indexing corner cases that arise in NumPy testing
- Implemented part of numpy.fft (cffti and cfftf)

Looking forward, we plan to continue improving the correctness of the existing implemented NumPy functionality, while also beginning to look at performance. The initial focus for performance will be to look at areas where we are significantly worse than CPython+NumPy. Those interested in trying these improvements out will need a PyPy nightly, and an install of the PyPy NumPy fork. Thanks again to the NumPy on PyPy donors for funding this work.

Tuesday, December 10, 2013

NumPy Status Update - November

Since the PyPy 2.2 release last month, more progress has been made on the NumPy compatibility front. Initial work has been directed by running the NumPy test suite and targeting failures that appear most frequently, along with fixing the few bugs reported on the bug tracker.

Improvements were made in these areas:
- Many missing/broken scalar functionalities were added/fixed. The scalar API should match up more closely with arrays now.
- Some missing dtype functionality was added (newbyteorder, hasobject, descr, etc)
- Support for optional arguments (axis, order) was added to some ndarray functions
- Fixed some corner cases for string/record types

Most of these improvements went onto trunk after 2.2 was split, so if you're interested in trying them out or running into problems on 2.2, try the nightly.

Thanks again to the NumPy on PyPy donors who make this continued progress possible.

Cheers,
Brian

Saturday, November 30, 2013

PyPy Leysin Winter Sprint (11-19st January 2014)

The next PyPy sprint will be in Leysin, Switzerland, for the ninth time. This is a fully public sprint: newcomers and topics other than those proposed below are welcome.

Goals and topics of the sprint

  • Py3k: work towards supporting Python 3 in PyPy
  • NumPyPy: work towards supporting the numpy module in PyPy
  • STM: work towards supporting Software Transactional Memory
  • And as usual, the main side goal is to have fun in winter sports :-) We can take a day off for ski.

Exact times

For a change, and as an attempt to simplify things, I specified the dates as 11-19 January 2014, where 11 and 19 are travel days. We will work full days between the 12 and the 18. You are of course allowed to show up for a part of that time only, too.

Location & Accomodation

Leysin, Switzerland, "same place as before". Let me refresh your memory: both the sprint venue and the lodging will be in a very spacious pair of chalets built specifically for bed & breakfast: http://www.ermina.ch/. The place has a good ADSL Internet connexion with wireless installed. You can of course arrange your own lodging anywhere (as long as you are in Leysin, you cannot be more than a 15 minutes walk away from the sprint venue), but I definitely recommend lodging there too -- you won't find a better view anywhere else (though you probably won't get much worse ones easily, either :-)

Please confirm that you are coming so that we can adjust the reservations as appropriate. The rate so far has been around 60 CHF a night all included in 2-person rooms, with breakfast. There are larger rooms too (less expensive per person) and maybe the possibility to get a single room if you really want to.

Please register by Mercurial:

https://bitbucket.org/pypy/extradoc/
https://bitbucket.org/pypy/extradoc/raw/extradoc/sprintinfo/leysin-winter-2014

or on the pypy-dev mailing list if you do not yet have check-in rights:

http://mail.python.org/mailman/listinfo/pypy-dev

You need a Swiss-to-(insert country here) power adapter. There will be some Swiss-to-EU adapters around -- bring a EU-format power strip if you have one.

Friday, November 15, 2013

NumPy status update

Here is what has been happening with NumPy in PyPy in October thanks to the people who donated to the NumPyPy proposal:

The biggest change is that we shifted to using an external fork of numpy rather than a minimal numpypy module. The idea is that we will be able to reuse most of the upstream pure-python numpy components, replacing the C modules with appropriate RPython micronumpy pieces at the correct places in the module namespace.

The numpy fork should work just as well as the old numpypy for functionality that existed previously, and also include much new functionality from the pure-python numpy pieces that simply hadn't been imported yet in numpypy. However, this new functionality will not have been "hand picked" to only include pieces that work, so you may run into functionality that relies on unimplemented components (which should fail with user-level exceptions).

This setup also allows us to run the entire numpy test suite, which will help in directing future compatibility development. The recent PyPy release includes these changes, so download it and let us know how it works! And if you want to live on the edge, the nightly includes even more numpy progress made in November.

To install the fork, download the latest release, and then install numpy either separately with a virtualenv: pip install git+https://bitbucket.org/pypy/numpy.git; or directly: git clone https://bitbucket.org/pypy/numpy.git; cd numpy; pypy setup.py install.

EDIT: if you install numpy as root, you may need to also import it once as root before it works: sudo pypy -c 'import numpy'

Along with this change, progress was made in fixing internal micronumpy bugs and increasing compatibility:
  • Fixed a bug with strings in record dtypes
  • Fixed a bug where the multiplication of an ndarray with a Python int or float resulted in loss of the array's dtype
  • Fixed several segfaults encountered in the numpy test suite (suite should run now without segfaulting)

We also began working on __array_prepare__ and __array_wrap__, which are necessary pieces for a working matplotlib module.

Cheers,
Romain and Brian

Wednesday, September 25, 2013

Numpy Status Update

Hi everyone

Thanks to the people who donated money to the numpy proposal, here is what I've been working on recently :

- Fixed conversion from a numpy complex number to a python complex number
- Implement the rint ufunc
- Make numpy.character usable as a dtype
- Fix ndarray(dtype=str).fill()
- Various fixes on boolean and fancy indexing

Cheers
Romain

Tuesday, August 27, 2013

NumPy road forward

Hello everyone.

This is the roadmap for numpy effort in PyPy as discussed on the London sprint. First, the highest on our priority list is to finish the low-level part of the numpy module. What we'll do is to finish the RPython part of numpy and provide a pip installable numpypy repository that includes the pure python part of Numpy. This would contain the original Numpy with a few minor changes.

Second, we need to work on the JIT support that will make NumPy on PyPy faster. In detail:

  • reenable the lazy loop evaluation
  • optimize bridges, which is depending on optimizer refactorings
  • SSE support

On the compatibility front, there were some independent attempts into making the following stuff working:

  • f2py
  • C API (in fact, PyArray_* API is partly present in the nightly builds of PyPy)
  • matplotlib (both using PyArray_* API and embedding CPython runtime in PyPy)
  • scipy

In order to make all of the above happen faster, it would be helpful to raise more funds. You can donate to PyPy's NumPy project on our website. Note that PyPy is a member of SFC which is a 501(c)(3) US non-profit, so donations from US companies can be tax-deducted.

Cheers,
fijal, arigo, ronan, rguillebert, anto and others


Thursday, August 8, 2013

NumPyPy Status Update

Hello everyone

As expected, nditer is a lot of work. I'm going to pause my work on it for now and focus on simpler and more important things, here is a list of what I implemented :
  • Fixed a bug on 32 bit that made int32(123).dtype == dtype("int32") fail
  • Fixed a bug on the pickling of array slices
  • The external loop flag is implemented on the nditer class
  • The c_index, f_index and multi_index flags are also implemented
  • Add dtype("double") and dtype("str")
  • C-style iteration is available for nditer
Cheers
Romain Guillebert

Monday, June 3, 2013

NumPyPy status update

Hello everyone,

May was the first month I was paid to work on NumPyPy (thanks to all who donated!), here is what I worked on during this period :

  • It is now possible to use subarrays.
  • It is now possible to pickle ndarrays (including those using subarrays), dtypes and scalars, the pickling protocol is the same as numpy's.

For June, I plan to work on the nditer class, it seems that there's enough work for an entire month.

Cheers
Romain Guillebert

Saturday, May 11, 2013

Numpy Status Update

Hello Everyone,

I've started to work on NumPyPy since the end of April and here is a short update :

  • I implemented pickling support on ndarrays and dtypes, it will be compatible with numpy's pickling protocol when the "numpypy" module will be renamed to "numpy".
  • I am now working on subarrays.

I would also like to thank everyone who donated and allowed me to work on this.

Cheers,
Romain Guillebert

Monday, March 18, 2013

Numpy status update and developer announcement

Hello, some good news!

First the update:

  • dtype support - NumPy on PyPy now supports non-native storage formats. Due to a lack of true support for longdoubles in rpython, we decided to back out the support of longdouble-as-double which was misleading.
  • missing ndarray attributes - work has been made toward supporting the complete set of attributes on ndarrays. We are progressing alphabetically, and have made it to d. Unsupported attributes, and unsupported arguments to attribute calls will raise a NotImplementedError.
  • pickling support for numarray - hasn't started yet, but next on the list
  • There has been some work on exposing FFI routines in numpypy.
  • Brian Kearns has made progress in improving the numpypy namespace. The python numpypy submodules now more closely resemble their numpy counterparts. Also, translated _numpypy submodules are now more properly mapped to the numpy core c-based submodules, furthering the goal of being able to install numpy as a pure-python module with few modifications.

And now the good news:

While our funding drive over 2012 did not reach our goal, we still managed to raise a fair amount of money in donations. So far we only managed to spend around $10 000 of it. We issued a call for additional developers, and are glad to welcome Romain Guillebert and Ronan Lamy to the numpypy team. Hopefully we will be able to report on speedier progress soon.

Cheers,
Matti Picus, Maciej Fijalkowski


Saturday, January 26, 2013

NumPyPy 2013 Developer Position

Introduction

Proposed herein is a part-time fellowship for developing NumPy in PyPy. The work will initially consist of 100 hours with the possibility of extension, until the funds run out. Development and improvement of PyPy's NumPyPy (as with most Open Source and Free Software) is done as a collaborative process between volunteer, paid, and academic contributors. Due to a successful funding drive but a lack of contributors willing to work directly for PyPy, we find ourselves in the enviable situation of being able to offer this position.

Background

PyPy's developers make all PyPy software available to the public without charge, under PyPy's Open Source copyright license, the permissive MIT License. PyPy's license assures that PyPy is equally available to everyone freely on terms that allow both non-commercial and commercial activity. This license allows for academics, for-profit software developers, volunteers and enthusiasts alike to collaborate together to make a better Python implementation for everyone.

NumPy support for PyPy is licensed similarly, and therefore NumPy in PyPy support can directly help researchers and developers who seek to do numeric computing but want an easier programming language to use than Fortan or C, which is typically used for these applications. Being licensed freely to the general public means that opportunities to use, improve and learn about how NumPy in PyPy works itself will be generally available to everyone.

The Need for a Part-Time Developer

NumPy project in PyPy has seen some slow, but steady progress since we started working about a year ago. On one hand, it's actually impressive what we could deliver with the effort undertaken, on the other hand, we would like to see the development accelerated.

PyPy has strict coding, testing, documentation, and review standards, which ensures excellent code quality, continually improving documentation and code test coverage, and minimal regressions. A part-time developer will be able to bring us closer to the goal of full numpy-api implementation and speed improvements.

Work Plan

The current proposal is split into two parts:

  • Compatibility:

    This part covers the core NumPy Python API. We'll implement most NumPy APIs that are officially documented and we'll pass most of NumPy's tests that cover documented APIs and are not implementation details. Specifically, we don't plan to:

    • implement NumPy's C API
    • implement other scientific libraries, like SciPy, matplotlib or biopython
    • implement details that are otherwise agreed by consensus to not have a place in PyPy's implementation of NumPy or agreed with NumPy community to be implementation details
  • Speed:

    This part will cover significant speed improvements in the JIT that would make numeric computations faster. This includes, but is not necesarilly limited to:

    • write a set of benchmarks covering various use cases
    • teaching the JIT backend (or multiple backends) how to deal with vector operations, like SSE
    • experiments with automatic parallelization using multiple threads, akin to numexpr
    • improving the JIT register allocator that will make a difference, especially for tight loops

    As with all speed improvements, it's relatively hard to predict exactly how it'll cope, however we expect the results to be withing an order of magnitude of handwritten C equivalent.

Position Candidate

We would like people who are proficient in NumPy and PyPy (but don't have to be core developers of either) to step up. The developer selection will be done by consensus of PyPy core developers and consulted with the Software Freedom Conservancy for lack of conflict of interest. The main criterium will be past contributions to the PyPy project, but they don't have to be significant in size.

A candidate for the Developer position will demonstrate the following:

  • The ability to write clear, stable, suitable and tested code
  • The ability to understand and extend the JIT capabilities used in NumPyPy.
  • A positive presence in PyPy's online community on IRC and the mailing list.

Ideally the Developer will also:

  • Have familiarity with the infrastructure of the PyPy project (including bug tracker and buildbot).
  • Have Worked to provide education or outreach on PyPy in other forums such as workshops, conferences, and user groups.

Conservancy and PyPy are excited to announce the Developer Position. Renumeration for the position will be at the rate of 60 USD per hour, through the Software Freedom Conservancy.

PyPy community is promising to provide necessary guidance and help into the current codebase, however we expect a successful candidate to be able to review code and incorporate external patches within two months of the starting date of the contract.

Candidates should submit their proposal (including their CV) to:

pypy-z@python.org

The deadline for this initial round of proposals is February 1, 2013.


Thursday, November 1, 2012

NumPy status update #5

Hello.

I'm quite excited to inform that work on NumPy in PyPy has been restarted and there has been quite a bit of progress on the NumPy front in PyPy in the past two months. Things that happened:

  • complex dtype support - thanks to matti picus, NumPy on PyPy now supports complex dtype (only complex128 so far, there is work on the other part)
  • big refactoring - probably the biggest issue we did was finishing a big refactoring that disabled some speedups (notably lazy computation of arrays), but lowered the barrier of implementing cool new features.
  • fancy indexing support - all fancy indexing tricks should now work, including a[b] where b is an array of integers.
  • newaxis support - now you can use newaxis features
  • improvements to ``intp``, ``uintp``, ``void``, ``string`` and record dtypes

Features that have active branches, but hasn't been merged:

  • float16 dtype support
  • missing ndarray attributes - this is a branch to finish all attributes on ndarray, hence ending one chapter.
  • pickling support for numarray - hasn't started yet, but next on the list

More importantly, we're getting very close to able to import the python part of the original numpy with only import modifications and running it's tests. Most tests will fail at this point, however it'll be a good start for another chapter :-)

Cheers,
fijal


Tuesday, September 4, 2012

NumPy on PyPy status update

Hello everyone.

It's been a while since we posted a numpy work update, but I'm pleased to inform you that work on it has been restarted. A lot of the work has been done by Matti Picus, who is one of the newest contributors to the PyPy project. None of the work below has been merged so far, it's work in progress:

  • Complex dtype support.
  • Fixing incompatibilities between numpy and pypy's version.
  • Refactoring numpypy to simplify the code and make it easier for new contributors.
  • Reuse most of the numpy's pure python code without modifications.

Finishing this is also the plan for the next month.

Cheers,
fijal


Tuesday, August 7, 2012

NumPyPy non-progress report

Hello everyone.

Not much has happened in the past few months with numpypy development. A part of the reason was doing other stuff for me, a part of the reason was various unexpected visa-related admin, a part of the reason was EuroPython and a part was long-awaited holiday.

The thing that's maybe worth mentioning is that it does not mean the donations disappeared in the mist. PyPy developers are being paid to work on NumPyPy on an hourly basis - that means if I decide to take holidays or work on something else, the money is simply staying in the account until later.

Thanks again for all the donations, I hope to get back to this topic soon!

Cheers,
fijal


Tuesday, April 17, 2012

NumPy on PyPy progress report

Hello.

A lot of things happened in March, like pycon. I was also busy doing other things (pictured), so apologies for the late numpy status update.

However, a lot of things have happened and numpy continues to be one of the main points of entry for hacking on PyPy. Apologies to all the people whose patches I don't review in timely manner, but seriously, you do a lot of work.

This list of changes is definitely not exhaustive, and I might be forgetting important contributions. In a loose order:

  • Matti Picus made out parameter work for a lot of (but not all) functions.

  • We merged record dtypes support. The only missing dtypes left are complex (important), datetime (less important) and object (which will probably never be implemented because it makes very little sense and is a mess with moving GCs).

  • Taavi Burns and others implemented lots of details, including lots of ufuncs. On the completely unscientific measure of "implemented functions" on numpypy status page, we're close to 50% of numpy working. In reality it might be more or less, but after complex dtypes we're getting very close to running real programs.

  • Bool indexing of arrays of the same size should work, leaving only arrays-of-ints indexing as the last missing element of fancy indexing.

  • I did some very early experiments on SSE. This work is seriously preliminary - in fact the only implemented operation is addition of float single-dimension numpy arrays. However, results are encouraging, given that our assembler generator is far from ideal:

     

    Numpy

    PyPy SSE

    PyPy

    GCC non-looped

    GCC looped

    a+b

    0.6s

    0.3s

    0.4s

    0.3s

    0.25s

    a+b+c

    1.9s

    0.35s

    0.5s

    0.7s

    0.32s

    a+b+c+d+e

    3.2s

    0.36s

    0.8s

    1.7s

    0.51s

    The benchmark repo is available. GCC was run with -O3, no further options specified. PyPy was run with default options, the SSE branch is under backend-vector-ops, but it's not working completely yet.

    One might argue that C and Python is not the same code - indeed it is not. It just shows some possible approach to writing numeric code.

Next step would be to just continue implementing missing features such as

  • specialised arrays i.e. masked arrays and matrixes
  • core modules such as fft, linalg, random.
  • numpy's testing framework

The future is hard to predict, but we're not far off!

Cheers,
fijal

UPDATE:Indeed, string and unicode dtypes are not supported yet. They're as important as complex dtype

Saturday, January 28, 2012

NumPyPy status update

Hello.

This is just a quick status update on the NumPy in PyPy project that very recently became my day job. I should give my thanks once again to Getco, Nate Lawson and other contributors who donated above $40000 towards the goal.

Recently we (Alex Gaynor, Matti Picus and me) implemented a few interesting things that a lot of people use:

  • more ufuncs
  • most ufuncs now accept the axis parameter (except all and any)
  • fixed string representation of arrays, now it's identical to numpy (uses pretty much the same code)
  • ndarray.flat should be working correctly
  • ndarray.flatten, ndarray.ravel, ndarray.take
  • indexing arrays by boolean arrays of the same size
  • and various bugfixes.

We would also like to introduce the nightly report of numpy status. This is an automated tool that does package introspection. While it gives some sort of idea how much of numpy is implemented, it's not by far the authority. Your tests should be the authority. It won't report whether functions support all kinds of parameters (for example masked arrays and out parameter are completely unsupported) or that functions work at all. We also reserve the right to incorporate jokes in that website, so don't treat it that seriously overall :-)

Thanks, and stay tuned. We hope to post here regular updates on the progress.

Cheers,
fijal & the PyPy team

Friday, January 27, 2012

Py3k and Numpy First Stage: Thanks to all who Gave

Last year was quite successful for PyPy fundraising through the Software Freedom Conservancy, and Conservancy and PyPy are very excited to announce that enough was raised to begin the first stages on the Py3k and Numpy grant proposals.

As of the end of 2011, 135 different individuals gave to the Py3k campaign, and 114 to the Numpy campaign. We thank each of you who donated to help make this work possible. Meanwhile, if you haven't given to support these projects, we do hope you'll give generously now to help fund their second stages later this year!

We're also particularly excited that a few donors gave particularly large donations to support this work; those big donations really filled in the gap to help us get started!

Specifically, we're pleased to announce that Google donated $35000 towards implementing Python 3 in PyPy. Google's general support of the Python community is well known, and their specific support of our grant proposal is much appreciated.

Meanwhile, Numpy was supported in part by contributions from Nate Lawson, Cantab Capital Partners, and Getco, as well as more than a hundred other contributors.

With these donations combined with many others, we're now starting work on both projects. This week, the Conservancy signed contracts with Antonio Cuni and Benjamin Peterson to work towards the Stage 1.1 goals in Py3k proposal (and is negotiating for another contractor as well), and with Maciej Fijałkowski to work towards the Stage 1 goals in the Numpy proposal.

In 2012, PyPy will continue regular sprint meetings, at which Py3K and Numpy efforts will certainly have a place. We have some limited funds to fund travels of contributors to those meetings.

We're very thankful for all who donated so far to support these efforts, and we hope that now that work has begun, even more donors will come forward to help us finish the job. In the meantime, watch for the commits showing up from these developers and other contributors in the PyPy repositories!

Cheers, The PyPy Team

Tuesday, January 10, 2012

NumPyPy progress report - running benchmarks

Hello.

We're excited to let you know about some of the great progress we've made on NumPyPy: both completeness and performance. In this blog entry we mostly will talk about performance and how much progress we have made so far.

Word of warning: this work is in progress -- we're maybe half way to where we want to be and there are many trivial and not so trivial optimizations to be written. (For example, we haven't even started to implement important optimizations, like vectorization.)

Benchmark

We chose a laplace equation solver, based on SciPy's PerformancePython wiki. Unfortunately, the different implementations on the wiki page accidentally use two different algorithms, which have different convergences, and very different performance characteristics on modern computers. As a result, we implemented our own versions in both C and Python (with and without NumPy). The full source can be found in fijal's hack repo, all these benchmarks were performed at revision 18502dbbcdb3.

First, let me describe various algorithms used. Note that some of them contain PyPy-specific hacks to work around limitations in the current implementation. These hacks will go away eventually and the performance will improve. Numerically the algorithms used are identical, however exact data layout in memory differs between them.

A note about all the benchmarks: they each were run once, but the performance is very stable across runs.

Starting with the C version, it implements a trivial laplace transform using two loops and double-reference memory (array of int*). The double reference does not matter for performance and the two algorithms are implemented in inline-laplace.c and laplace.c. They were both compiled with gcc 4.4.5 at -O3. The inline version modifies array in-place while the non-inline version stores results in a copy. That makes them converge at different rate, hence different number of iterations

A straightforward version of those in Python is implemented in laplace.py using, respectively, inline_slow_time_step and slow_time_step. slow_2_time_step does the same thing, except it copies arrays in-place instead of creating new copies. Table below compares running PyPy against C:

bench number of iterations time per iteration
laplace C 219 6.3ms
inline-laplace C 278 20ms
slow python 219 17ms
slow 2 python 219 14ms
inline_slow python 278 23.7ms

An important thing to notice is the data dependency of the inline version causes a huge slowdown for the C versions. This is not a severe disadvantage for us though -- the brain-dead Python version takes longer and PyPy is not able to take advantage of the knowledge that the data is independent. The results are in the same ballpark as the C versions -- 15% - 170% slower, but the algorithm one chooses matters more than the language. By comparison, the slow versions take about 5.75s each on CPython 2.6 per iteration and, by estimation, are about 200x slower than the PyPy equivalent, if I had the patience to measure the full run.

The next step is to use NumPy expressions. The first problem we run into is that computing the error requires walking the entire array a second time. This is fairly inefficient in terms of cache access, so I took the liberty of computing the errors every 15 steps. This results in the convergence being rounded to the nearest 15 iterations, but speeds things up considerably. numeric_time_step takes the most braindead approach of replacing the array with itself, like this:

u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 +
                       (u[1:-1,0:-2] + u[1:-1, 2:])*dx2)*dnr_inv

We need 3 arrays here -- one is an intermediate (PyPy only needs one, for all of those subexpressions), one is a copy for computing the error, and one is the result. This works automatically because in NumPy + or * creates an intermediate, while NumPyPy avoids allocating the intermediate if possible.

numeric_2_time_step works in pretty much the same way:

src = self.u
self.u = src.copy()
self.u[1:-1, 1:-1] = ((src[0:-2, 1:-1] + src[2:, 1:-1])*dy2 +
                      (src[1:-1,0:-2] + src[1:-1, 2:])*dx2)*dnr_inv

except the copy is now explicit rather than implicit.

numeric_3_time_step does the same thing, but notice one doesn't have to copy the entire array, it's enough to copy the border pieces and fill rest with zeros:

src = self.u
self.u = numpy.zeros((self.nx, self.ny), 'd')
self.u[0] = src[0]
self.u[-1] = src[-1]
self.u[:, 0] = src[:, 0]
self.u[:, -1] = src[:, -1]
self.u[1:-1, 1:-1] = ((src[0:-2, 1:-1] + src[2:, 1:-1])*dy2 +
                      (src[1:-1,0:-2] + src[1:-1, 2:])*dx2)*dnr_inv

numeric_4_time_step is the one that tries hardest to resemble the C version. Instead of doing an array copy, it actually notices that one can alternate between two arrays. This is exactly what the C version does. The remove_invalidates call is a PyPy specific hack - we hope to remove this call in the near future, but, in short, it promises "I don't have any unbuilt intermediates that depend on the value of the argument", which means one doesn't have to compute sub-expressions one is not actually using:

remove_invalidates(self.old_u)
remove_invalidates(self.u)
self.old_u[:,:] = self.u
src = self.old_u
self.u[1:-1, 1:-1] = ((src[0:-2, 1:-1] + src[2:, 1:-1])*dy2 +
                      (src[1:-1,0:-2] + src[1:-1, 2:])*dx2)*dnr_inv

This one is the most comparable to the C version.

numeric_5_time_step does the same thing, but notices one doesn't have to copy the entire array, it's enough to just copy the edges. This is an optimization that was not done in the C version:

remove_invalidates(self.old_u)
remove_invalidates(self.u)
src = self.u
self.old_u, self.u = self.u, self.old_u
self.u[0] = src[0]
self.u[-1] = src[-1]
self.u[:, 0] = src[:, 0]
self.u[:, -1] = src[:, -1]
self.u[1:-1, 1:-1] = ((src[0:-2, 1:-1] + src[2:, 1:-1])*dy2 +
                      (src[1:-1,0:-2] + src[1:-1, 2:])*dx2)*dnr_inv

Let's look at the table of runs. As before, gcc 4.4.5, compiled at -O3, and PyPy nightly 7bb8b38d8563, on an x86-64 machine. All of the numeric methods run for 226 steps, slightly more than the 219, rounding to the next 15 when the error is computed.

benchmark PyPy CPython
numeric 21ms 35ms
numeric 2 14ms 37ms
numeric 3 13ms 29ms
numeric 4 11ms 31ms
numeric 5 9.3ms 21ms

We think that these preliminary results are pretty good. They're not as fast as the C version (or as fast as we'd like them to be), but we're already much faster than NumPy on CPython -- almost always by more than 2x on this relatively real-world example. This is not the end, though. In fact, it's hardly the beginning! As we continue work, we hope to make even more use of the high level information that we have. Looking at the assembler generated by gcc for this example, it's pretty clear we can outperform it thanks to better aliasing information and hence better possibilities for vectorization. Stay tuned.

EDIT: fixed the benchmark name

EDIT2: added info that first table is about PyPy

Cheers, fijal