Registration is now open for the 5th BrainScaleS/FACETS CodeJam workshop, which will take place March 14th-16th, 2011 in Edinburgh.

The goal of the CodeJam workshops is to catalyze open-source, collaborative software development in computational and systems neuroscience and neuroinformatics (especially, but not exclusively, using Python), by bringing together researchers, students and engineers to share ideas, present their work, and write code together. The general format of the workshops is to dedicate the mornings to invited and contributed talks, leaving the afternoons free for discussions and code sprints.

For the 5th BrainScaleS/FACETS CodeJam, the main theme of the meeting will be convergence in computational neuroscience software: recent developments to promote interoperability of modelling, simulation and data analysis tools and future efforts to develop common tools and libraries. We are planning sessions on:

• Neuronal network modelling
• Code generation for neuron and synapse models
• Multicompartmental neuron modelling in Python
• Data analysis tools for computational and systems neuroscience
• Best practices for running an open-source software project

More details on the program and invited speakers will follow soon.

We invite contributions on any topic related to software in neuroscience, but especially on topics related to the main theme and planned sessions. If you have ideas for organising code sprints, whether a feature that you would like to see added to an existing tool or an idea for new software, please also let us know.

The meeting is being organised by Andrew Davison, Mike Hull, Abigail Morrison, Eilif Muller, Miha Pelko and Laurent Perrinet.

Registration & Further Information

The registration deadline is 19th February 2012, and is limited to 50 participants.

Please consult the meeting website at http://neuralensemble.org/meetings/CodeJam5 for registration and further information.

• The FE quadrature is not accurate enough
• The condition number of the matrix is high, thus LAPACK doesn't return very accurate eigenvalues
• Bug in the assembly/solver (like single/double corruption in Fortran, or some other subtle bug)

It’s been a while since I’ve started considering buying a real photo camera. And I’ve decided on a Canon 7D. As usual, the user guide covers the camera, but I thought it would be better to have a more complete guide to the 7D.

Content and opinions

You don’t have to read the book to know how to take a picture, but each chapter adds a new angle o the way one can use the 7D. The first chapter is a ten configuration tips list, and I think each one of them is relevant. They have already helped me enhance my pictures. The second chapter is a list of things to check and understand before shooting. I have to say I didn’t apply them all (updating the firmware for instance), and there is also an overview of the basic technical elements in a camera (shutter, iso…).

Then the serious work begins. The third chapter browses through the different capture modes. Each mode is explained with the associated available changes (exposure, ISO…). Then two chapters tackle what I would say are the most common styles: portrait and landscape. I learnt about a lot of things offered by the 7D in those chapters, like exposure metering or the electronic level. There are a lot of tips that are available in those two chapters. The sixth chapter is about moving targets, a subject that is hard, but it seems that the camera has some useful functions for that.

Then, two chapters tackle subjects that are common to all modes: light and composition. There are more advanced books on these topics, but the author manages to cover the most useful ground, from hints to warnings. The last chapter but one is on video, it is interesting, but I didn’t get my 7D for video. Still, there is a warning that I did find very interesting on this subject. I won’t tell it here, go and grab the book instead.

The ending chapter is about some “advanced” techniques, but it’s more about explaining some things you can do with the 7D, but there is 90% chance that you won’t ever use any of them.

Finally, it is worth mentioning that each chapter starts with one or two pictures and some associated comments. This is a great help when you don’t have much experience. After each chapter, a small list sums up all new information.

Conclusion

This book is a perfect complement to the user guide. The latter explains one button or setting after the other, whereas the former links everything together. When you don’t know about the different modes of the camera (aperture or shutter for instance), it is a must have to take great pictures.

Just after the 0.1 release, I’ve worked to add some few tricks and fix a few bugs (see QtMosaic on Github). The most important change is a better image search.

The Antipole Tree

An antipole tree is a binary search tree that is constructed on antipole pairs. An antipole pair is composed of the furthest points in a set. When constructing the tree, one starts with the full set of images/thumbnails, then recursively, an antipole pair is chosen from the (sub)set (figure 1) and two new clusters are created, images are aggregated to the nearest image of the pair (figure 2). When the subset is small enough, the recursion stops.

Figure 1: Original set with the first antipole pair

Figure 2: Two subsets based on the antipole pair in the original set

The search uses a distance-sorted list. When a leaf node is visited, one looks for the nearest image and stores the result. When an internal node is visited, its children will be added to the sorted list, and the distance will be the distance to the center of the cluster minus the radius (figure 3). This way, if the current best is nearer to the reference than the nearest cluster, no nearest image will be found. It’s a usual distance search in a binary tree.

Figure 3: distance from a new image to the different subsets

Result

Of course, building the tree is costly. Although it can be parallelized, it is still slow. The upper side is that using a tree for the search is O(log(n)) instead of the usual O(n). It is actually more efficient for mosaic databases with several thousands of images or more. The result is of course exactly the same as with the linear search.

There are some differences between the original publication (see the reference) and the current code. In the reference, the antipole pair has a special status, as the center of the clusters and they are not tested with the rest of the subset in a leaf node. So in the original code, each tested cluster tests at least one image with the reference, whereas only leaf nodes give the closest images search in my case. The difference is that the center of each cluster will be the real center of the cluster, meaning that the minimum distance will be more close to the reality, and I hope that this will translate into speed-up. The other addition is that the visiting code is also simpler, which is better for maintenance.

Conclusion

This tool has now a logarithmic image search. It should be faster than linear search for huge mosaic databases. Also, the final photomosaic was enhanced by changing the result mosaic mean to the original image mosaic.

Reference : Antipole Tree Indexing to Support Range Search and K-Nearest Neighbor Search in Metric Spaces

Someone recently asked me about creating cross-tabulations and contingency tables using pandas. I wrote a bit about this in October after implementing the pivot_table function for DataFrame. So I thought I would give a few more examples and show R code vs. the equivalent pandas code which will be helpful for those already used to the R way of doing things.

When calling help(xtabs) (or help(table)) in R, you get a number of examples for using table, xtabs, and ftable:

> head(esoph)
  agegp     alcgp    tobgp ncases ncontrols
1 25-34 0-39g/day 0-9g/day      0        40
2 25-34 0-39g/day    10-19      0        10
3 25-34 0-39g/day    20-29      0         6
4 25-34 0-39g/day      30+      0         5
5 25-34     40-79 0-9g/day      0        27
6 25-34     40-79    10-19      0         7

> ftable(xtabs(cbind(ncases, ncontrols) ~ agegp + tobgp, data = esoph))
                ncases ncontrols
agegp tobgp                     
25-34 0-9g/day       0        70
      10-19          1        19
      20-29          0        11
      30+            0        16
35-44 0-9g/day       2       109
      10-19          4        46
      20-29          3        27
      30+            0        17
45-54 0-9g/day      14       104
      10-19         13        57
      20-29          8        33
      30+           11        19
55-64 0-9g/day      25       117
      10-19         23        65
      20-29         12        38
      30+           16        22
65-74 0-9g/day      31        99
      10-19         12        38
      20-29         10        20
      30+            2         4
75+   0-9g/day       6        26
      10-19          5        11
      20-29          0         3
      30+            2         4

> ftable(xtabs(cbind(ncases, ncontrols) ~ agegp, data = esoph))
       ncases ncontrols
agegp                  
25-34       1       116
35-44       9       199
45-54      46       213
55-64      76       242
65-74      55       161
75+        13        44

Here are the same examples using pandas:

In [19]: table = esoph.pivot_table(['ncases', 'ncontrols'], 
                                    rows=['agegp', 'tobgp'], 
                                   aggfunc=np.sum)

In [20]: table
Out[20]: 
                ncases  ncontrols
agegp tobgp                      
25-34 0-9g/day  0.      70.      
      10-19     1.      19.      
      20-29     0.      11.      
      30+       0.      16.      
35-44 0-9g/day  2.      109      
      10-19     4.      46.      
      20-29     3.      27.      
      30+       0.      17.      
45-54 0-9g/day  14      104      
      10-19     13      57.      
      20-29     8.      33.      
      30+       11      19.      
55-64 0-9g/day  25      117      
      10-19     23      65.      
      20-29     12      38.      
      30+       16      22.      
65-74 0-9g/day  31      99.      
      10-19     12      38.      
      20-29     10      20.      
      30+       2.      4.0      
75+   0-9g/day  6.      26.      
      10-19     5.      11.      
      20-29     0.      3.0      
      30+       2.      4.0     

In [22]: table2 = esoph.pivot_table(['ncases', 'ncontrols'], 
                                    rows='agegp', aggfunc=np.sum)
In [23]: table2
Out[23]: 
       ncases  ncontrols
agegp                   
25-34  1.      116      
35-44  9.      199      
45-54  46      213      
55-64  76      242      
65-74  55      161      
75+    13      44.

One of the really great things about pandas is that the object produced is still a DataFrame (the R object is of a special class)! So we could do anything with it that we do with normal DataFrame objects:

In [26]: table['ncases'].unstack('agegp')
Out[26]: 
agegp     25-34  35-44  45-54  55-64  65-74  75+
tobgp                                           
0-9g/day  0      2      14     25     31     6  
10-19     1      4      13     23     12     5  
20-29     0      3      8.     12     10     0  
30+       0      0      11     16     2.     2

Here is another fun example:

In [31]: import pandas.rpy.common as com
In [32]: wp = com.load_data('warpbreaks')
In [33]: wp['replicate'] = np.tile(range(1, 10), 6)

In [34]: wp.pivot_table('breaks', rows=['wool', 'tension'], 
                        cols='replicate', aggfunc=np.sum)
Out[34]: 
replicate     1   2   3   4   5   6   7   8   9 
wool tension                                    
A    H        36  21  24  18  10  43  28  15  26
     L        26  30  54  25  70  52  51  26  67
     M        18  21  29  17  12  18  35  30  36
B    H        20  21  24  17  13  15  15  16  28
     L        27  14  29  19  29  31  41  20  44
     M        42  26  19  16  39  28  21  39  29

Here is the equivalent R code:

> warpbreaks$replicate <- rep(1:9, len = 54)
> ftable(xtabs(breaks ~ wool + tension + replicate, data = warpbreaks))
             replicate  1  2  3  4  5  6  7  8  9
wool tension                                     
A    L                 26 30 54 25 70 52 51 26 67
     M                 18 21 29 17 12 18 35 30 36
     H                 36 21 24 18 10 43 28 15 26
B    L                 27 14 29 19 29 31 41 20 44
     M                 42 26 19 16 39 28 21 39 29
     H                 20 21 24 17 13 15 15 16 28

As soon as we add a formula framework to statsmodels I could see incorporating that into this kind of function in Python.

I had the privilege of speaking last night at the NYCPython meetup group. I’ve given tons of “use pandas!” talks so I thought I would take a slightly different angle and talk about some of the design and implementation work that I’ve done for getting good performance in critical data manipulations. I’ll turn some of this material into some blog articles in the near future.

Wes McKinney: pandas design and development from Adam Klein on Vimeo.

Here’s some more video footable shot by my awesome friend Emily Paup with her HDSLR.

I did a little interactive demo (using the ever-amazing IPython HTML Notebook) on Ashley Williams’s Food Nutrient JSON Database:

Link to PDF output of Demo

Link to IPython Notebook file

If you want to run the code in the IPython notebook, you’ll have to download the food database file above.

The audience and I learned from this demo that if you’re after Tryptophan, “Sea lion, Steller, meat with fat (Alaska Native)” is where to get it (in the highest density).

Packt publishing sent me a copy of NumPy 1.5 Beginner’s guide by Ivan Idris.

The book actually covers more than only numpy: it is a full introduction to numerical computing with Python. The table of contents is the following:

• NumPy Quick Start
• Beginning with NumPy Fundamentals
• Get into Terms with Commonly Used Functions
• Convenience Functions for Your Convenience
• Working with Matrices and ufuncs
• Move Further with NumPy Modules
• Peeking Into Special Routines
• Assure Quality with Testing
• Plotting with Matplotlib
• When NumPy is Not Enough: SciPy and Beyond

The book is easy to read, as it requires no specific expertise other than knowing basic Python programming. It is full of examples and exercises, which is really great for learning. I find the style of the author, Ivan Idris, particularly amusing and relaxing, engaging the reader with questions, challenges, or even jokes (“Have a go hero”).

With regards to the formatting and the print, the book is written in large fonts, with sectioning information, tips and exercises clearly standing out.

It is full of practical information, such as how to install the software, or where to get help. Finally, One thing that I appreciated, is that the examples are typed in IPython. Each time I teach, I like to use IPython, because it is full of features to help plotting, debugging and profiling numerical code. The book even has a little introduction to some useful IPython features.

After an introduction to the work flow, the book explores array manipulation such as creation or reshaping, followed by some simple numerics and the battery of array-based operations on functions and polynomials. Then it presents linear algebra and signal processing basics (FFT). It also covers the financial functions that are present in numpy and mentions testing, which is very important to achieve quality code. The book finishes with matplotlib and scipy, two modules that are important to know to go further.

The examples are mostly drawn from statistics or financial applications, such as computing running averages on stock quotes. Basic math explanations, such as the definition of the Moore-Penrose pseudo-inverse, are given when needed.

To conclude, I enjoyed this book and I think that it is a nice addition to my library. It answers exactly it’s title: it is well-suited for beginners wanting to learn numpy. On the other hand, I would not recommend it as a reference material, or as a book to learn more general scientific or numerical computing with Python.

Joblib 0.6: better I/O and Python 3 support

Happy new year, every one. I have just released Joblib 0.6.0 beta. The highlights of the 0.6 release are a reworked enhanced pickler, and Python 3 support.

Many thanks go to the contributors to the 0.5.X series (Fabian Pedregosa, Yaroslav Halchenko, Kenneth C. Arnold, Alexandre Gramfort, Lars Buitinck, Bala Subrahmanyam Varanasi, Olivier Grisel, Ralf Gommers, Juan Manuel Caicedo Carvajal, and myself). In particular Fabian made sure that Joblib worked under Python 3.

In this blog post, I’d like to discuss a bit more the compressed persistence engine, as it illustrates well key factors in implementing and using compressed serialization.

Fast compressed persistence

One of the key components of joblib is it’s ability to persist arbitrary Python objects, and read them back very quickly. It is particularly efficient for containers that do their heavy lifting with numpy arrays. The trick to achieving great speed has been to save in separate files the numpy arrays, and load them via memmapping.

However, one drawback of joblib, is that the caching mechanism may end up using a lot of disk space. As a result, there is strong interest in having compressed storage, provided it doesn’t slow down the library too much. Another use case that I have in mind for fast compressed persistence, is implementing out of core computation.

There are some great compressed I/O libraries for Python, for instance Pytables. You may wonder why the need to code yet another one. The answer is that joblib is pure Python, depending only on the standard library (numpy is optional), but also that the goal here is black-box persistence of arbitrary objects.

Comparing I/O speed and compression to other libraries

Implementing efficient compressed storage was a bit of a struggle and I learned a lot. Rather than going into the details straight away, let me first discuss a few benchmarks of the resulting code. Benching such feature is very hard, first because you are fighting with the disk cache, second because they performances depends very much on the data at hand (some data compress better than others), last because they are three interesting metrics: disk space used, write speed, and read speed.

Dataset used - I chose to compare the different strategies on some datasets that I work with, namely the probabilistic brain atlases MNI 1mm (62Mb uncompressed) and Juelich 2mm (105Mb uncompressed). Whether the data is represented as a Fortran-ordered array, or a C-ordered array is important for the I/O performance. This data is normally stored to disk compressed using the domain-specific Nifti format (.nii files), accessed in Python with the Nibabel library.

Libraries used - I benched different compression strategies in joblib against Nibabel’s Nifti I/O, compressed or not, and against using Pytables to store the data buffer (without the meta-informations). Pytables exposed a variety of compression strategies, with different speed compromises. In addition, I benched numpy’s builtin save_compressed.

I would like to stress that I am comparing a general purpose persistence engine (joblib) to specific I/O libraries either optimized for the data (Nifti), or requiring some massaging to enable persistence (pytables).

Comparing to other libraries

Actual numbers can be found here.

Take home messages - The graphs are not crystal-clear, but a few tendencies appear:

• Pytables with LZO or blosc compression is the king of the hill for read and write speed.
• I/O of compressed data is often faster than with uncompressed data for a good compression algorithm.
• Joblib with Zlib compression level 1 performs honorably in terms of speed with only the Python standard library and no compiled code.
• Read time of memmapping (with nibabel or joblib) is negligeable (it is tiny on the graphs), however the loading time appears when you start accessing the data.
• Passing in arrays with a memory layout (Fortran versus C order) that the I/O library doesn’t expect can really slow down writing.
• Compressing with Zlib compression-level 1 gets you most of the disk space gains for a reasonable cost in write/read speed.
• Compressing with Zlib compression-level 9 (not shown on the figures) doesn’t buy you much in disk space, but costs a lot in writing time.

Benching datasets richer than pure arrays

The datasets used so far are pretty much composed of one big array, a 4D smooth spatial map. I wanted to test on more datasets, to see how the performances varied with data type and richness. For this, I used the datasets of the scikit-learn, real life data of various nature, described here:

• 20 news - 20 usenet news group: this data mainly consists of text, and not numpy arrays.
• LFW people - Labeled faces in the wild, many pictures of different people’s face.
• LFW pairs - Labeled faces in the wild, pairs of pictures for each individual. This is a high entropy dataset, it does not have much redundant information.
• Olivetti - Olivetti dataset: centered pictures of faces.
• Juelich(F) - Our previous Juelich atlas
• Big people - The LFW people dataset, but repeated 4 times, to put a strain on memory resources.
• MNI(F) - Our previous MNI atlas
• Species - Occurence of species measured in latin America, with a lot of missing data.

Actual numbers can be found here.

What this tells us - The main message from these benchmarks is that datasets with redundant information, i.e. that compress well, give fast I/O. This is not surprising. In particular, good compression can give good I/O on text (20 news). Another result, more of a sanity check, is that compressed I/O on big data (Big people, ) works as well as on smaller data. Earlier code would start to swap. Finally, I conclude from these graphs, that compression levels from 1 to 3 buy you most of the gains for reasonable costs, and that going up to 9 is not recommended, unless you know that your data can be compressed a lot (species).

Lessons learned

I’ll keep this paragraph short, because the information is really in joblib’s code and comments. Don’t hesitate to have a look, it’s BSD-licenced, so you are free to borrow what you please.

1. Memory copies, of arrays, but also of strings and byte streams can really slow you down with big data.
2. To avoid copies with numpy arrays, fully embrace numpy’s strided memory model. For instance, you do not need to save arrays in C order, if they are given to you in a different order. Accessing the memory in the wrong striding direction explains the poor write performance of pytables on Fortran-ordered Juelich.
3. When dealing with the file system, the OS makes so much magic (e.g. prefetching) that clever hacks tend not to work: always benchmark.
4. Depending on the size of the data, it may be more efficient to store subsets in different files: it introduces ‘chunk’ that avoid filling in the memory too much (parameter cache_size in joblib’s code). In addition, data of a same nature tends to compress better.
5. The I/O stream or file object interfaces are abstractions that can hide the data movement and the creation of large temporaries. After experiments with GZipFile and StringIO/BytesIO I found it more efficient to fall back to passing around big buffer object, numpy arrays, or strings.
6. For reasons 4 and 5, I ended up avoiding the gzip module: raw access to the zlib with buffers gives more control. This explains a good part of the differences in read speed for pure arrays with numpy’s save_compressed.

One of my conclusions for joblib, is that I’ll probably use Pytables as an optional backend for persistence in a future release.

Details on the benchmarks

These benchmarks where run on a Dell Lattitude D630 laptop. That’s a dual-core Intel Core2 Duo box, with 2M of CPU cache.

The code for the benchmarks below can be found on a gist.

Thanks

I’d like to that Francesc Alted for very useful feedback he gave on this topics. In particular, the following thread on the pytables mailing-list may be of interest to the reader.

Happy new year, friends!

I’ve made a New Year’s resolution to build a better web presence and make better use of the domain that I previously only used for mail.

This has prompted me to move my blog over to http://blog.vene.ro which hopefully is shorter, better, faster and stronger!

Over the last week I have completely retooled pandas’s “database” join infrastructure / algorithms in order to support the full gamut of SQL-style many-to-many merges (pandas has had one-to-one and many-to-one joins for a long time). I was curious about the performance with reasonably large data sets as compared with base::merge.data.frame which I’ve used many times in R. So I just ran a little benchmark for joining a 100000-row DataFrame with a 10000-row DataFrame on two keys with about 10000 unique key combinations overall. Simulating a somewhat large SQL join.

Note this new functionality will be shipping in the upcoming 0.7.0 release (!).

There is a major factor affecting performance of the algorithms in pandas, namely whether the result data needs to be sorted lexicographically (which is the default behavior) by the join columns. R also offers the same option so it’s completely apples to apples. These are mean runtimes in seconds:

As you can see, the sorting time in pandas completely dominates the runtime of the join (sorting 10000 strings is a lot of string comparisons). This isn’t really an indictment of R’s merge algorithm, but rather speaks to the strength of the merge strategy I devised in pandas. After spending a few days working on this problem I definitely think R could do a lot better. I haven’t run benchmarks against SQLite or another SQL database yet; that will happen eventually.

I’m going to write a blog article in the next few days going into algorithmic detail about how I got merging pretty big data sets to be so fast–it was not easy!

Here is the R code.

And the Python code

Do you know how fast your code is? Is it faster than it was last week? Or a month ago? How do you know if you accidentally made a function slower by changes elsewhere? Unintentional performance regressions are extremely common in my experience: it’s hard to unit test the performance of your code. Over time I have gotten tired of playing the game of “performance whack-a-mole”. Thus, I started hacking together a little weekend project that I’m calling vbench. If someone thinks up a cleverer name, I’m all ears.

Link to pandas benchmarks page produced using vbench

What is vbench?

vbench is a super-lightweight Python library for running a collection of performance benchmarks over the course of your source repository’s history. Since I’m a GitHub user, it only does git for now, but it could be generalized to support other VCSs. Basically, you define a benchmark:

Then you write down the information about your repository and how to build any relevant DLLs, etc., that vary from revision to revision:

Then you pass this info, plus a list of your benchmark objects, to the BenchmarkRunner class:

Now, the BenchmarkRunner makes a clone of your repo, then runs all of the benchmarks once for each revision in the repository (or some other rule, e.g. I’ve set run_option='eod' to only take the last snapshot on each day). It persists the results in a SQLite database so that you can rerun the process and it will skip benchmarks it’s already run (this is key when you add new benchmarks, only the new ones will be updated). Benchmarks are uniquely identified by the MD5 hash of their source code.

This is the resulting plot over time for the above GroupBy benchmark related to some Cython code that I worked on late last week (where I made a major performance improvement in this case):

Here is a fully-formed vbench suite in the pandas git repository.

Kind of like codespeed and speed.pypy.org?

Before starting to write a new project I looked briefly at codespeed and the excellent work that the PyPy guys have done with speed.pypy.org. But then I started thinking, you know, Django web application, JSON requests to upload benchmark results? Seemed like far too much work to do something relatively simple. The dealbreaker is that codespeed is just a web application. It doesn’t actually (to my knowledge, someone correct me if I’m wrong?) have any kind of a framework for orchestrating the running of benchmarks throughout your code history. That is what this new project is for. I actually see a natural connection between vbench and codespeed, all you need to do is write a script to upload your vbench results to a codespeed web app!

At some point I’d like to build a simple web front end or wx/Qt viewer for the generated vbench database. I’ve never done any JavaScript, but it would be a good opportunity to learn. Knowing me, I might break down and hack out a stupid little wxPython app with an embedded matplotlib widget anyway.

Anyway, I’m really excited about this project. It’s very prototype-y at the moment but I tried to design it in a nice and extensible way. I also plan to put all my git repo analysis tools in there (like code churn graphs etc.) so it should become a nice little collection of tools.

Other ideas for extending

1. "sudo apt-get install sloccount".
2. "cp -rv SAGE_ROOT/devel/sage-main /tmp/x"
3. Use a script [1] to rename all .pyx and .pxi files to .py.
4. Ran "sloccount *" in the /tmp/x directory, which ignores autogenerated .c/.cpp files coming from Cython.

Totals grouped by language (dominant language first):
python:      530370 (96.41%)
ansic:        14538 (2.64%)
cpp:           5188 (0.94%)

SLOC    Directory       SLOC-by-Language (Sorted)
88903   rings           python=87720,cpp=1183
72913   combinat        python=71629,cpp=1284
47747   schemes         python=46255,cpp=1492
39815   graphs          python=28377,ansic=11438
31540   matrix          python=31540
31019   modular         python=31012,ansic=7
24475   libs            python=21171,ansic=2845,cpp=459
20517   misc            python=20383,ansic=134
18006   interfaces      python=18006
17577   geometry        python=16936,cpp=641
12775   categories      python=12775
12093   server          python=12093
11971   groups          python=11971
11961   plot            python=11961
10686   crypto          python=10686
9920    modules         python=9920
8389    symbolic        python=8260,cpp=129
8150    algebras        python=8150
7260    ext             python=7198,ansic=62
7093    structure       python=7093
6364    coding          python=6364
5670    functions       python=5670
5249    homology        python=5249
4798    numerical       python=4798
4323    quadratic_forms python=4323
3919    gsl             python=3919
3911    calculus        python=3911
3879    sandpiles       python=3879
3003    sets            python=3003
2647    databases       python=2647
2074    logic           python=2074
1736    finance         python=1736
1608    games           python=1608
1465    monoids         python=1465
1435    tests           python=1383,ansic=52
1370    stats           python=1370
971     interacts       python=971
959     tensor          python=959
906     lfunctions      python=906
308     parallel        python=308
275     probability     python=275
219     media           python=219
197     top_dir         python=197

#!/usr/bin/env python

import os, shutil

for dirpath, dirnames, filenames in os.walk('.'):
    for f in filenames:
        if f.endswith('.pyx') or f.endswith('.pxi'):
            print f
            shutil.move(os.path.join(dirpath, f),
                        os.path.join(dirpath, os.path.splitext(f)[0] + '.py'))

During SuperComputing 2011, I stumbled across a paper on parallel random number generation. And it brought back some memories of a blog post I wanted to do here on oversampling in my Interactive Raytracer.

Random numbers are important in a raytracer because one needs oversampling to have anti-aliasing, and to have a more realistic rendering, we need to add some jitter inside the oversampling. I won’t talk about this here, it is not my point (yet). So I’ve decided to resurrect IRT from the dead, at least long enough to test parallel random number generators. It is now on Github, and I’ve also changed the matrix library to Eigen. Some optimizations later, it is slightly faster than before, although it is still very far from Bikker’s work.

The Python wrappers are still working BTW.

Let’s hope I find some time to finish this before getting back at VSTs

A while back, Bob L. Sturm blogged about a similar implementation of OMP to the one in scikit-learn. Instead of using the Cholesky decomposition like we did, his Matlab code uses the QR decomposition, to a similar (or maybe even identical) outcome, in theory. So lucky that Alejandro pointed out to him the existence of the scikit-learn implementation, and that Bob’s code exposed a bug that all the test coverage didn’t catch! This plot should increase, certainly not decrease! Something is clearly wrong here.

Luckily we were able to find it and fix it very quickly. I have updated the old entries I wrote on the OMP optimizations, so they no longer include the bug. But I take this opportunity to explain what exactly went wrong.

A key part of the optimization was that slicing out arbitrary columns out of an array is slow when they are passed to BLAS functions like matrix multiplication. In order to make the most out of your code, the data should have a contiguous layout. We achieved this by swapping active dictionary atoms (columns) to the beginning of the array.

Something that can happen, but won’t happen very often, is that after an atom is selected as active, the atom that takes its place after swapping needs to be selected. This is rare because dictionaries have many columns, out of which only very very few will be active. But when it happens, because the code didn’t keep track of swapped indices, the corresponding coefficient of the solution would get updated twice, leading to more zero entries than we should have. A keen eye could have noticed that the first `n_nonzero_coefs` entries in OMP solution vectors were never non-zero. But alas, my eye was not a keen one at all.

In other words, the following test (that was written after the bug was found, unfortunately) was failing:

def test_swapped_regressors():
    gamma = np.zeros(n_features)
    # X[:, 21] should be selected first, then X[:, 0] selected second,
    # which will take X[:, 21]'s place in case the algorithm does
    # column swapping for optimization (which is the case at the moment)
    gamma[21] = 1.0
    gamma[0] = 0.5
    new_y = np.dot(X, gamma)
    new_Xy = np.dot(X.T, new_y)
    gamma_hat = orthogonal_mp(X, new_y, 2)
    gamma_hat_gram = orthogonal_mp_gram(G, new_Xy, 2)
    # active indices should be [0, 21], but prior to the bugfix
    # the algorithm would update only [21] but twice
    assert_equal(np.flatnonzero(gamma_hat), [0, 21])
    assert_equal(np.flatnonzero(gamma_hat_gram), [0, 21])

Note that this bug has been fixed for a while, but I didn’t get the free time to write this post until now. Good news is: we fixed it, and did so very quickly after the report. So you can still trust me, I guess!

The NIPS conference: time for a sprint. The NIPS conference, one of the major conferences in machine learning, is hosted in Granada this year. I believe that it is the first time that it is hosted in Europe. As many of the scikit-learn developers are part of the wider NIPS community, but also many live in Europe, we jumped on the occasion to organize a truly international sprint: the NIPS 2011 scikit-learn sprint.

Finding money. As often with open source development, a lot of our contributors are young people, investing their free time outside of any request from their hierarchy. In such a situation, it can be hard to find travel money. So we started looking for sponsors. We needed to find a decent sum of money, as we were flying people in from places such as the West coast of the US, or even Japan. The good news is that we found money, and between supervisors pitching in, universities giving travel grants, and our generous sponsors, there will be an impressive list of contributors from all over the world at the sprint.

Thanks to our sponsors. The first people that we need to thank are Google, who gave us a sizable sponsorship, and the PSF, who made Google’s sponsorship possible through their accounting and sprints programs. We also need to thanks our other sponsors, namely Tinyclues. Thanks to these sponsors, and additional investment from many universities and research group, we have been able to gather a total of 12 contributors in Granada, a handful coming from overseas. Also, we are indebted to the University of Granada, and the Gnu/Linux Granada Group (GGG), who are providing hosting for the sprint, as well as Régine Bricquet, from INRIA, who did a lot of the trip planing for the sponsored people.

I am very much looking forward to the sprint. It will be the first time that meet in real life many of the contributors, and judging by the warmness of the on-line exchanges, it will be a great moment. Besides, Granada is known to be a lively and historical city.

If you are around and want to join us, to work on Python in machine learning, send us a mail on the mailing list.

1. Commercial server hosting:  Some Amazon EC2 instances
2. An employee (or later, employees) to maintain the servers: at the beginning, this would be me in my free time, since I have a lot of experience with this.
3. Advertising that sagenb.com exists and we want users!:  Unlike sagenb.org, we would very strongly encourage as many people as possible to use sagenb.com.   The advertising and landing page would explain that though sagenb.com generates money, all that money is all given back to the Sage development community (see below).
4. Hire employees to improve the notebook: Fix bugs, implement features, etc. There would be a public commitment that all such work would be open sourced and get included with Sage. This would include adding support for a for-pay "pro" subscription version, adding nicer "offline support" (via a Sage install on the user's computer), integrated spreadsheets and better data visualization and manipulation tools for Science and business applications, and enabling development of Sage's core library and patches submission entirely through the notebook, etc.

A little experiment to see what low rank approximation looks like. These are the best rank-k approximations (in the Frobenius norm) to the a natural image for increasing values of k and an original image of rank 512.

Python code can be found here. GIF animation made using ImageMagic’s convert script.

• Large performance and stability improvements
• PySide support (PyQt is no longer exclusively required)
• PyQt API v1 and v2 support
• New profiler plugin (thanks to Santiago Jaramillo, a new contributor)
• Experimental support for IPython v0.11+
• And many other changes

I’ll be presenting this paper in a couple of weeks at SC11. It’s the most up-to-date white paper about pandas.

PDF Download Link

pandas: a Foundational Python Library for Data Analysis and Statistics

A colleague of mine give me this book, as I use “third generation” VCSs. I decided to check on the author approach on Version Control and his opinion on the matter. The book explores the different approaches of the latest VCS tools, with their advantages and drawbacks. Also, it delves into some algorithmic designs of Distributed VCSs. I’ve already discussed some of these tools, but the book is not a flame war against one DVCS, but more an explanation of all of them.

Content and opinions

The focus on the book is how to use a Version Control system, so several chapters are dedicated to how to use the most well-known one – and the one the author develops.

The book is split in three parts: “old” VCSs, new VCSs and finally how they work. The first chapter is a glossary of the special VCS words and their meaning. Once this is done, the first example chapter starts with Subversion. A small situation is depicted with two coders that are locally closed.

The third chapter starts the second part and it introduces new concepts for distributed VCSs. A complete chapter is dedicated to the advantages of DVCSs, and it doesn’t fall in the common pitfalls (those are in fact a reason to some humor from the author). The next chapter does exactly the opposite, i.e. explains the drawbacks of a DVCS. Although I disagree on GUI part (there are wonderful GUIs for both Mercurial and GIT on Windows, and on Linux they have a great integration in Eclipse for instance), I do agree with the author on the main topics here. The next four chapters are the practical application seen in the second chapter applied to Mercurial, Git and Veracity (two chapter for this one, one to explain the differences with Mercurial and Git, and one for the scenario), the tool the author develops.

The last part can be called “how to use them”. Even for the second generation tools like SVN, I was flabbergasted to see that a lot of coders didn’t even know how to use it, and even worse, that experienced managers didn’t know how to implement an appropriate workflow. So I think that the first chapter in this part should be something that everyone HAS to read if one wants to be called a coder. Seriously. Then, as I’ve said, the author explains how the distributed tools work and some of their internals. It’s not for the average coder, but it’s really interesting – for the culture. The last chapter is a set of golden rules, and like the first chapter, it’s something mandatory to read and use.

Conclusion

One of the most easy-to-read introductory books and a book on a sensible subject, this book should be read by everyone searching for a robust tool to support one’s development. It is also a very good argumentation on why you should use one.

Over the last several months, I’ve invested a great deal in the GroupBy and indexing infrastructure of pandas. GroupBy in particular could still use more work to be made even higher performance, but even for quite large data sets, most users will find it more than satisfactory compared with other Python alternatives (which are mainly DIY approaches these days) or even in other data analysis packages, e.g. those in R, which I’ll talk a bit about in this article.

An easy application of GroupBy is creating pivot tables from a tabular data set. This involves aggregating a data set on some number of keys, then reshaping it to form a 2D table with the stratified aggregated values. I added a new function “pivot_table“ which provides a high level interface for creating pivot tables. Let’s consider the tips data set used in some examples of the Hadley Wickham’s excellent R reshape2 package:

Since most readers won’t have rpy2, you can get the data file here and read it with:

All of this new functionality is available now on GitHub but will be a part of the upcoming pandas 0.5.0 release. Note that I am only using rpy2 here to pull the R data set into my Python sesson– all of the below code is implemented in pure Python (with a little Cython magic to make things fast).

Pivot table basics

To use pivot_table, pass the pandas DataFrame and specify the column you wish to aggregate and the columns to group on the rows and columns of the table. Note the default aggregation is mean, though this can be specified:

Note that the returned object is still a DataFrame, so you can use all of the standard indexing technology to select portions of the resulting table:

If you don’t specify a particular values column, it will try to aggregate all the other columns:

Other aggregations can be performed, like using len to get group counts:

Once you’ve done the aggregation, you’re free to use the built-in reshaping capability to reshape the data in the table:

If you know a bit more about groupby, you can get clever and pass a dict of functions to perform multiple aggregations at once:

I thought that was pretty cool.

Comparison with reshape2

R’s reshape2 package relies on two functions, melt and cast, to do the reshaping. Its implementation is, at least in my opinion, a slightly more ad hoc approach than in pandas (the actual pivot_table function is only about 10 lines and basically falls out of the groupby and hierarchical-indexing based reshaping). However, it can take advantage of R’s built-in formulas which makes the syntax very intuitive. Here are some of the same examples using reshape2:

R lacks hierarchical indexing as far as I know, so reshape2 munges the group labels together into a row label.

In R there is also the xtabs function for doing cross-tabulation:

ftable just creates a flat table which can be more easily viewed. This table could be created with pandas by:

Performance

pandas is very fast as I’ve invested a great deal in optimizing the indexing infrastructure and other core algorithms related to things such as this. I don’t have a lot of points of comparison, but here is a simple benchmark of reshape2 versus pandas.pivot_table on a data set with 100000 entries and 25 groups. Here is the R code for the benchmark:

Here is what the actual table looks like:

Here’s the equivalent Python code (sans timings, which I’ll do in IPython):

Similarly:

And the R timing

And pandas:

So it’s about 3-4x faster than reshape2. This obviously depends on the data set and how you structure the aggregation. Here’s a couple slightly different benchmarks on the same data:

and pandas:

I suspect one source of slowdown reshape/reshape2 approach is that the melt operation is very expensive, as is variable subsetting. However, it enables easier computation of group margins, something which I have not yet addressed (but will, eventually) in pandas.

I’m pleased to announce the 1.0 version of QtSimpleEQ, a plugin with one low-pass, two peak and one high pass second-order filters. Nothing fancy in the algorithms, it’s mainly another show case for Qt VST plugins.

The code is available under the GPL2 on github and on Sourceforge.

The plugin can be download on the Sourceforge project page.

The plugin was tested with Tracktion 3 (Windows XP).

Snapshot:

QtSimpleEQ UI

In scipy’s development version there’s a new function closely related to the QR-decomposition of a matrix and to the least-squares solution of a linear system.

What this function does is to compute the QR-decomposition of a matrix and then multiply the resulting orthogonal factor by another arbitrary matrix. In pseudocode:

but unlike this naive implementation, qr_multiply is able to do all this without explicitly computing the orthogonal Q matrix, resulting both in memory and time saving. In the following picture I measured the memory consumption as a function of time of running this computation on a 1.000 x 1.000 matrix X and a vector Y (full code can be found here):

It can be seen that not only qr_multiply is almost twice as fast as the naive approach, but also that the memory consumption is significantly reduced, since the orthogonal factor is never explicitly computed.

Credit for implementing the qr_multiply function goes to Martin Teichmann.

Just for fun, I’m pleased to announce a working version of QtMosaic.
It allows to create a mosaics database and then generate a photomosaic from this database.

For instance, here is an original image:

Original photo

Rendered photomosaic

The code is available on GitHub, as usual, and it has been tested on Windows and Linux.

After a brief (!) absence, we’re back with a new and shiny version of scikits.image, the image processing toolbox for SciPy.

This release runs under all major operating systems where Python (>=2.6 or 3.x), NumPy and SciPy can be installed.

For more information, visit our website

http://scikits-image.org

or the examples gallery at

http://scikits-image.org/docs/0.3/auto_examples/

New Features

• Shortest paths
• Total variation denoising
• Hough and probabilistic Hough transforms
• Radon transform with reconstruction
• Histogram of gradients
• Morphology, including watershed, connected components
• Faster homography transformations (rotations, zoom, etc.)
• Image dtype conversion routines
• Line drawing
• Better image collection handling
• Constant time median filter
• Edge detection (canny, sobel, etc.)
• IO: Freeimage, FITS, Qt and other image loaders; video support.
• SIFT feature loader
• Example data-sets

… as well as many bug fixes and minor updates.

Contributors for this release

• Martin Bergtholdt
• Luis Pedro Coelho
• Chris Colbert
• Damian Eads
• Dan Farmer
• Emmanuelle Gouillart
• Brian Holt
• Pieter Holtzhausen
• Thouis (Ray) Jones
• Lee Kamentsky
• Almar Klein
• Kyle Mandli
• Andreas Mueller
• Neil Muller
• Zachary Pincus
• James Turner
• Stefan van der Walt
• Gael Varoquaux
• Tony Yu

Permalink | Leave a comment  »

One year ago I had the chance to take a class on Monte Carlo simulation with prof. Ion Văduva, and my assignment for the class was to implement exactly what it says in the title of the blog post. I am going to walk you through the idea behind this.

General formulation

The ratio-of-uniforms is a method that can be applied to many density functions. Essentially, given a density function over \( \mathbb{R}^m\), \( f(x) = \frac{h(x)}{H}\) where \(H\) is a normalization constant (ie. \( h(x) \geq 0\), \( H = \int h(x)dX\)). Given a parameter \( c > 0 \) and a parametrization \(\phi\) from \( \mathbb{R}^{m+1}\) to \( \mathbb{R}^{m}\) expressed as: $$ \phi(v_0, v_1, …, v_m) = \left ( \frac{v_1}{v_0^c}, \frac{v_2}{v_0^c}, …, \frac{v_m}{v_0^c} \right )$$
Define the set \( \mathcal{C} = \{\mathbf{v} \big | \gamma(\mathbf{v}) \leq 0, v_0 > 0\} \in \mathbb{R}^{m+1}\) where
$$\gamma(\mathbf{v}) = \log v_0 – \frac{\log h(\phi(v_0, v_1, …, v_m))}{mc + 1}$$ If \( \mathcal{C}\) is bounded and we sample a uniform vector \( \mathbf{V} \sim \text{Uniform}(\mathcal{C})\) then \( \phi(\mathbf{V}) \sim f(x)\). Also note that the measure (volume) of the set \( \mathcal{C}\) is \( \frac{H}{mc + 1}\). I do not have any references for the proof, except for a book in Romanian, but if you are interested, just leave me a comment and I’ll do a follow-up post with the proofs.

Univariate scenario

For the univariate case, all the above simplifies to $$ \mathcal{C} = \left \{(u, v) \Big | 0 u \sqrt
{h\left (\frac{v}{u^c}\right )} \right \} $$. We generate \( (U, V) \sim \text{Uniform}(\mathcal{C})\) and take \( \frac{V}{U^c} \sim f(x)\).
Since we are looking at the (univariate) Gamma distribution, described by: $$ f(x; \nu, \theta) = \frac{x^{\mu - 1} \exp(\frac{-x}{\theta})}{\theta^k\Gamma(k)}$$ \( \nu\) is the shape parameter and \( \theta\) is the scale parameter.
But because of the property that if \( X \sim \text{Gamma}(\nu, \theta)\), then for any \( k > 0\), \( kX \sim \text{Gamma}(\nu, k\theta)\), we conclude that we can fix \( \theta\) to 1 without loss of generality. Replacing in the style of the definition in the previous section, we have \( h(x; \nu) = x^{\nu-1}e^{-x}\) and \( H_\nu = \Gamma(\nu)\).
This allows us to compute the equation of the boundary of the set \( \mathcal{C}\) which ends up being described by \(\gamma(u, v) = \log{u} – \frac{\nu – 1}{c + 1} \log{\left(\frac{v}{u^c}\right)} + \frac{1}{c+1} \frac{v}{u^c}\). For visualisation purposes, here is how it would look like for \( \nu=6, c=1\) (plotted using Wolfram Alpha):

Sampling algorithm

In order to uniformly sample from this set, we can apply basic rejection sampling: just uniformly sample from a rectangular region surrounding the set, and reject the points that do not satisfy the condition. In order to do this as efficiently as possible, we need to compute the minimal bounding box, which can be done by solving a couple of optimization problems using Lagrange multipliers and the KKT conditions. Also by looking closely at the image, you can see that the lower left corner is exactly the origin: this turns out not to be a coincidence. I won’t go into detail here, but here are the bounds I derived:
$$ 0 u (\nu - 1)^\frac{\nu - 1}{c + 1} e ^ {-\frac{\nu - 1}{c + 1}} \text{ and } 0 v \left(\frac{c\nu + 1}{c}\right)^{\frac{c\nu + 1}{c + 1}} e ^ {- \frac {c\nu + 1}{c+1}}$$

The probability of acceptance (which can be seen as the efficiency) of the rejection sampling method is given by the ratio of the areas of the set \( \mathcal{C}\) and the bounding box. The larger this probability, the less points we throw away and the more efficient the algorithm is. Using the values derived above, this probability is: $$ p(\nu, c) = \frac{\Gamma(\nu)e^{\nu}}{(c+1) (\nu - 1)^{\frac{\nu - 1}{c + 1}} \left(\frac{c\nu + 1}{c}\right)^{\frac{c\nu + 1}{c + 1}}}$$

Personally I got stumped here. The idea would be to determine the ideal \( c\) for a given \( \nu\) in order to maximize the probability, but I didn't manage to do it (I leave it as an exercise for the reader ). Anyway, this is enough to proceed with an implementation, so I'm gonna give the Python code for it. Note that I used the name k for the shape parameter instead of \( \nu\). Also note that the case when \( 0 \nu 1\) needed to be treated separately, which I did using the following property: Let \( \nu \in (0, 1)\). If \( X' \sim \text{Gamma}(1+\nu, 1), U \sim \text{Uniform}(0, 1)\) then $$ X = X' \cdot \sqrt[\nu]{U} \sim \text{Gamma}(\nu, 1)$$ For a proof of this fact, see [1], which is a great article on generating Gamma variates.

Implementation

from import numpy as np

def _cond(u, v, k, c):
    """Identity function describing the acceptance region"""
    x = v / u ** c
    return (c + 1) * np.log(u) <= (k - 1) * np.log(x) - x

def vn_standard_gamma(k, c=1.0, rng=np.random):
    """Generates a single standard gamma random variate"""
    if k <= 0:
        raise ValueError("Gamma shape should be positive")
    elif k < 1:
        return vn_standard_gamma(1 + k, c, rng) * rng.uniform() ** (1 / k)
    elif k == 1:
        return rng.standard_exponential()
    else:
        a, b = get_bounds(k, c)
        while True:
            u, v = rng.uniform(0, a), rng.uniform(0, b)
            if _cond(u, v, k, c):
                break;
        return v / u ** c

def vn_gamma(k, t, shape=1, c=1.0, rng=np.random):
    """Vectorized function to generate multiple gamma variates"""
    generator = lambda x: t * vn_standard_gamma(k, c, rng)
    generator = np.vectorize(generator)
    return generator(np.empty(shape))

def get_bounds(k, c=1.0):
    """Computes the minimal upper bounds surrounding the acceptance region"""
    a = ((k - 1) / np.e) ** ((k - 1) / (c + 1))
    b = ((c * k + 1) / (c * np.e)) ** ((c * k + 1) / (c + 1))
    return a, b

def prob_acc(k, c=1.0):
    """Calculates the probability of acceptance for the given parameters"""
    from scipy.special import gamma
    a, b = get_bounds(k, c)
    return gamma(k) / ((c + 1) * a * b)

Results

And of course I should show you that it works. Here are some histograms for various values of \( \nu\), with the theoretical density plotted in dotted red, after sampling \( 10^5\) values. The y-axis is the frequency (sorry for labeling in Romanian), and for the red dotted line it can be interpreted as the theoretical probability. You can clearly see the goodness of fit is excellent.

[1]: George Marsaglia and Wai Wan Tsang. 1998. The Monty Python method for generating random variables. ACM Trans. Math. Softw. 24, 3 (September 1998), 341-350. DOI=10.1145/292395.292453

In case you didn’t know it – today is Ada Lovelace Day!

Now, as any self-respecting Computer Science degree-wielding person should, I, too, think it’s important to celebrate the day named after the world’s very first programmer.

For me, the first math teacher I remember making a big difference was Shirley Theis – who taught me Algebra in 8th grade at McKinley Middle School in Redwood City, CA. Mrs Theis, an energetic dynamo in her mid fifties, was a deeply motivated and caring teacher, who expected a lot out of her students, but never in a disciplinary manner. She was full of enthusiasm, which projected out and infected even the most timid or disaffected student: in her class, you couldn’t be just a sack of potatoes planted in your seat.

She often lead class in a nearly theatrical manner – pacing back and forth, egging students on by eagerly repeating their partial responses, getting exponentially more excited if the student was on the right track, barely containing herself from jumping up and down in anticipation of that lightbulb going off — and yet just as quickly waning in her enthusiasm,becoming a personified caricature of hopelessness and despair to let you know the instant a response was starting to go astray.

It may have been the only math class I’ve ever taken where there were group assignments – we would work with a partner or a few classmates in trying to figure out an assignment, first trying it solo, and then putting our heads together to figure out why our answers disagree and which is the right one. I believe it was Mrs. Theis who succinctly captured a value I hold in high regard: “it’s not about how far you go – it’s about how many people you bring with you.”

There was one other mathematics teacher I had in my life who clearly stands out: it was Professor Evelyn Silvia who had a comparable level of enthusiasm and energy, and from whom I had the pleasure of taking the first upper-division math course (Math 108 – Intro to Abstract Math) during my second quarter at UC Davis. Dr. Silvia was the real deal – she cared, gesticulated, encouraged us to question why something was true, and had an approach to the demanded we each take ownership of our education. The book for the course, Introduction to Abstract Mathematics: A Working Excursion by D.O. Cutler and E.M. Silvia was a blue workbook – each of us had our own copy, and there were blanks left out for us to write our own answers to the exercises. The fact that the book had blanks for me to fill in was so inviting, there was a kind of “working mathematician” approach that came with it with that it made me really enjoy and look forward to working through the material. I still have mine.

Dr. Silvia was incredibly sharp, not just intellectually but also interpersonally. Not only could she could gauge when the class was lost, but she also had a knack for spotting if something was affecting you outside of class. She was really committed to helping you not just as a student, but as a person. I remember spending hours at Mishka’s, or Cafe Roma, or the CoHo, reading and writing, wanting to do well and not let Silvia down, because she invested so much energy in placed a great deal of trust in us.

So thank you both, Shirley Theis and Evelyn Silvia – you both encouraged me to grow a lot as a person, challenged my concept of what it means to be a student, and by your example provided a template of what it means to be an effective teacher, which I’ve imitated and embraced with pleasure in my own teaching.

(tagged scipy to spread word of Ada Lovelace day to Planet SciPy)

Today, I wanted to announce my new plugin, a 4-bands EQ, but when I started a test with pyVST, I encountered strange things:

• The first is my fault, as the code of the EQ disappeared from my Git repository, so I have to code it again. Mainly it is just plugin the correct actions between the filters and the GUI.
• The second is an error at the end of the test. I’ve updated my Qt version from 4.7.1 to 4.7.4, and since this update (or perhaps since I updated to Python 2.7 for pyVST), I found that even a recompiled QtSimpleOverdrive has the same behavior. It did not when I released it. It seems that Qt is complaining about events that are bounced between different threads, but the actual error message is more cryptic, and impossible to debug the application at this point.

I hope to fix these mistakes this month, I really hope I can get QtVST to work again.

Last week we released a new version of scikit-learn. The Changelog is particularly impressive, yet personally this release is important for other reasons.

This will probably be my last release as a paid engineer. I’m starting a PhD next month, and although I plan to continue contributing to the project and make a few more releases, I will certainly have less time to devote to it. Luckily, I received a lot of help from the community while preparing the release, from Changelog writing to build of Windows binaries, thus I expect the transition to go smoothly.

Almost two years have elapsed since the first 0.1 release. During this time, we did a lot of refactoring and broke the API several times. However, I’ve seen some concerns about API stability both at the EuroScipy conference and in the mailing list where I’ve realized we need to provide an API that does not break in every release, and do this in a way that the project remains fun for developers.

That’s why I’m extremely glad to see that although this release is big in changes, these have been made in a more organized manner. Yes, we’ve broken the API once again, but now there’s a compatibility layer that ensures that code written for 0.8 will continue working with the new release.

R is much more magical than Python. What do I mean by this? In R, things like this are a part of everyday life:

If you’re a seasoned Python programmer, you might have the sort of visceral negative reaction that I do to this. Seriously, just where in the hell did those variable names come from? So when I say magic here I’m talking about abusing the language’s parser. There is nothing special about R that makes the above behavior possible, but rather taking a fundamentally different design philosophy to, say, Python. As any Python programmer knows: Explicit is better than implicit. I happen to agree. There is also a bit of a semantic difference in R versus Python in that assignment in R typically copies data, whereas variables in Python are simply references (labels) for a particular object. So you could make the argument that the names a and b above are more strongly linked to the underlying data.

While building pandas over the last several years, I occasionally grapple with issues like the above. Maybe I should just break from Python ethos and embrace magic? I mean, how hard would it be to get the above behavior in Python? Python gives you stack frames and the ast module after all. So I went down the rabbit hole and wrote this little code snippet:

While this is woefully unpythonic, it’s also kind of cool:

This can even parse and format more complicated expressions (harder than it looks, because you have to walk the whole AST):

Now, I am *not* suggesting we do this any time soon. I’m going to prefer the explicit approach (cf. the Zen of Python) any day of the week:

Colleagues who are exposing a numerical C code in Python asked me for some advice on the best way to pass arrays from C to Python avoiding copies. They had Cython in mind, and I must agree with them that I have found the Cython code to be more maintainable than hand-written Python C-API code.

When writing my answer, I found out that there was no self-contained example of creating numpy arrays from existing data in Cython. Thus I created my own. The full code with readme build and demo scripts is available on a gist. Here I only give an executive summary.

The core functionality is implemented by the PyArray_SimpleNewFromData function of the C API of numpy that can create an ndarray from a pointer to the data, a simple data type, and the shape of the data. The Cython file just builds around that function:

Many scientific Python users are surprised when I tell them that ndarray.take is faster than __getitem__-based (a.k.a. “fancy” as I call it) indexing.

It’s actually kind of unbelievable when you think about it. What’s going on here that take is almost 10x faster? I really should take a closer at the internals of what __getitem__ does because this has always struck me as pretty bad. Maybe I shouldn’t be complaining? I mean, R 2.13′s indexing falls somewhere in the middle:

So 656 microseconds per iteration. (In an earlier version of this post I used rpy2 to do the benchmark and got 1.05 ms, but there was apparently some overhead from rpy2)

Another peculiarity that I noticed with take is that performance gets worse when you use the out argument, which tells the function to use an array you pass in to write out the result:

EDIT: I’ve been informed that using mode='clip' or mode='wrap' makes this run as fast as without the out argument.

Weird! I was dissatisfied by this, so I got curious how fast a hand-coded little Cython function can do this:

Don’t worry about the -1 thing– that’s a specialization that I’m using inside pandas. Curiously, this function is a lot faster than take using out but faster than the regular take by a handful of microseconds.

Very interesting.

TL;DR

Last week was marked by the international RANLP (Recent Advances in Natural Language Processing) conference, taking place in a nice spa in Hissar, Bulgaria. The excellent folks from the computational linguistics group at the University of Wolverhampton were behind it, together with the Institute of Information and Communication Technologies from the Bulgarian Academy of Sciences.

A fountain with warm mineral spring water in Hissarya.
Picture by Miranda

The schedule was busy, with three tracks going in parallel, in order to cover a wide range of topics in computational linguistics. The student workshop should also be noted for the excellent quality of the works there.

Of course it would be infeasible to write about all the great people I met and their high quality work. And if I were to write about all the fun we had, it would probably make this post look unprofessional . This doesn’t mean I forgot about any of you, and as soon as I get the chance to work on something related, I will most certainly write about it, and you.

So, if I would have to summarize the trends and the ideas stated during the conference and especially during the keynotes, I would say:

• When talking about word sense disambiguation, it’s wrong to speak about the different meanings of a word, but rather about the potential a word has for bringing a certain meaning in a certain context. See Patrick Hanks‘ Corpus Pattern Analysis. Without something like this, to have good WSD you need to heavily adjust the overlapping meanings from a Wordnet-style ontology.
• Certain relations, such as temporal and spacial ones, can naturally be modeled by complex domain-specific logics (see Inderjeet Mani‘s new book, Interpreting Motion: Grounded Representations for Spatial Language, that is due for publishing). But these only appear in a small subset of human communication. The attempt to map human language to a complete logic in which to do general-purpose inference seems futile: Ido Dagan suggests textual entailment: reasoning directly in natural language, and only abstracting away to a formal logic system when need arises.
• If you have a large enough sample of n-gram frequency data, you can eventually beat the performance you can get with a limited amount of labeled data, and most importantly: it generalizes much better when going out of the domain you trained on. Apparently the best tool for this at the moment is the Google n-gram data, which has some limitations. In time, we can easily extend this data by huge amounts by mining n-grams from Wikipedia (which allegedly has a higher count of distinct n-grams than the Google dataset), and more importantly, by aligning multi-language data, making use of transliterations and cognate identification.

Please note that I may be (or most probably, am) ignorant of older instances of similar ideas, and I may have misunderstood certain claims. Please feel free to discuss in the comments whether you think I forgot about something important, or whether I am plain wrong about something. In particular, I seem to have been completely ignorant of the existence of the Google n-gram data, which has been around for quite a while, so I must have missed other important things as well.

Take care, kind readers, and express your opinion!
V

It has been a while since my last post here, but I’m back! I had access to the French version of this book, thanks to the publisher.

CUDA is now in the trend, and there are several books, one of them I’ve also reviewed.

Content and opinions

The first two chapters are dedicated to introducing CUDA in computing science: why CUDA, why it can help someone and what nVidia offers. Of course, now nVidia offers more than what the book tells us (even nVidia offered more when the book was published if my memory serves me right).

Actual CUDA code is introduced next. The author adds small pieces in each chapter: first the specific syntax, then parallelism and threads, blocks and their cooperation with shared memory, constant memory and streams, texture memory.

The last step is the optimization of the code. The author starts with the interaction with OpenGL (and DirectX is left as an exercise), then atomicity, streams and finally multi GPUs. Some of these are only available on new GPUs, but it is not mentioned explicitly (you have to test your card with the provided code).

The last chapter browses through the different supports that are available: CUDA tools, documentation…

There are a lot of stuff that are “missing”. The most important thing is an explanation of how memory should be accessed. Coalescing memory is not even mentioned, and although nVidia made huge improvements since the first CUDA GPU, this is still one of the biggest ways to optimize your code. Also, the debuggers are mentioned, but there are sadly no real tutorial on this. You will have to fill some voids in your new knowledge (for instance the texture memory has more functionalities than what is described).

Conclusion

So CUDA By Example is a good introductory book with flaws. Contrary to the other book I read, this one (unfortunately?) does not dive into architectural details. You will still have to learn from the nVidia’s tutorials, but you will have an almost complete overview of CUDA.

Thanks to Olivier, Gaël and Alex, who reviewed the code heavily the last couple of days, and with apologies for my lack of activity during a sequence of conferences, Dictionary learning has officially been merged into scikit-learn master, and just in time for the new scikit-learn 0.9 release. Here are some glimpses of the examples you can run for yourself:

The stars of this new release are: the manifold learning module by Jake Vanderplas and Fabian Pedregosa, the Dirichlet process gaussian mixture model by Alexandre Passos, and many others, as you can see from the development changelog (as soon as the release is made, I will update this post with permanent links).

The release is due tomorrow. I will also be in charge with building the Windows installers for this release, let’s hope I do a good job and you can think of me and smile when installing!

Top notch scientific conferences are starting to add Python tracks to their program. This is good news. Indeed, it scientific Python conferences (namely Scipy, EuroSciPy and Scipy India) are doing great to get together people who have already heard about Python for science, but we need to reach out to specific Python communities to maximize impact.

ESCO 2012 - European Seminar on Coupled Problems

ESCO 2012 is the 3rd event in a series of interdisciplineary meetings dedicated to computational science challenges in multi-physics and PDEs.

I was invited as ESCO last year. It was an aboslute pleasure, because it is a small conference that is very focused on discussions. I learned a lot and could sit down with people who code top notch PDE libraries such as FEniCS and have technical discussions. Besides, it is hosted in the historical brewery where the Pilsner was invented. Plenty of great beer.

Application areas Theoretical results as well as applications are welcome. Application areas include, but are not limited to: Computational electromagnetics, Civil engineering, Nuclear engineering, Mechanical engineering, Computational fluid dynamics, Computational geophysics, Geomechanics and rock mechanics, Computational hydrology, Subsurface modeling, Biomechanics, Computational chemistry, Climate and weather modeling, Wave propagation, Acoustics, Stochastic differential equations, and Uncertainty quantification.

Minisymposia

• Multiphysics and Multiscale Problems in Civil Engineering
• Modern Numerical Methods for ODE
• Porous Media Hydrodynamics
• Nuclear Fuel Recycling Simulations
• Adaptive Methods for Eigenproblems
• Discontinuous Galerkin Methods for Electromagnetics
• Undergraduate Projects in Technical Computing

Software afternoon Important part of each ESCO conference is a software afternoon featuring software projects by participants. Presented can be any computational software that has reached certain level of maturity, i.e., it is used outside of the author’s institution, and it has a web page and a user documentation. If you would like to present your software project, let us know soon.

Proceedings For each ESCO we strive to reserve a special issue of an international journal with impact factor. Proceedings of ESCO 2008 appeared in Math. Comput. Simul., proceedings of ESCO 2010 in CiCP and Appl. Math. Comput. Proceedings of ESCO 2012 will appear in Computing.

Important Dates

• December 15, 2011: Abstract submission deadline.
• December 15, 2011: Minisymposia proposals.
• January 15, 2012: Notification of acceptance.

PyHPC: Python for High performance computing

If you are doing super computing, SC11, the Super Computing conference is the reference conference. This year there will a workshop on high performance computing with Python: PyHPC.

At the scipy conference, I was having a discussion with some of the attendees on how people often still do process management and I/O with Fortran in the big computing environment. This is counter productive. However, has success stories of supercomputing folks using high-level languages are not advertized, this is bound to stay. Come and tell us how you use Python for high performance computing!

Topics

• Python-based scientific applications and libraries
• High performance computing
• Parallel Python-based programming languages
• Scientific visualization
• Scientific computing education
• Python performance and language issues
• Problem solving environments with Python
• Performance analysis tools for Python application

Papers We invite you to submit a paper of up to 10 pages via the submission site. Authors are encouraged to use IEEE two column format.

Important Dates

• Full paper submission: September 19, 2011
• Notification of acceptance: October 7, 2011
• Camera-ready papers: October 31, 2011

1. I have fixed the PyQt input hook issue on Windows platforms. Before this revision, the console was crashing (non-responsive actually) on Windows platforms after trying to manipulate Qt objects interactively, like interactive plotting with Matplotlib for example. Now, when enabling the "Replace PyQt input hook by Spyder's" option (which is enabled by default on Windows platforms), the standard Python interpreter may be used to do interactive plotting, even on Windows.
2. Plus, I have implemented some basic special commands in the standard Python interpreter (see here more details), adding support for %pwd, %ls, %clear or !dir (or !ls on UNIX platforms), etc. (I haven't implemented the %edit, %run and other similar commands which have really no use within Spyder).
3. And I have fixed a bug with PYTHONSTARTUP substitution option (the packages imported in the startup script were not available in the Python interpreter).

• Matplotlib's startup script:

For all my Python zealotry, I never claimed it was perfect. In fact, there are a number of places you can shoot yourself in the foot if you aren’t careful. But I was kind of bummed about this one:

Try it.

EDIT: I’ve been informed by an astute reader that the GC is actually capable of dealing with the above case. Yang Zhang says “The collector actually works with cycles, including this case. (Actually, most objects are managed via reference-counting; the gc is specifically for cyclic references.) The problem is that it’s triggered by # allocations, not # bytes. In this case, each cats array is large enough that you need several GB of allocations before the threshold is triggered. Try adding gc.collect() or gc.get_count() in the loop (and inspect gc.get_threshold()).”

All is not lost however. This pattern can be implemented using weak references:

At the request of a friend, I am putting up some of the posters that I recently presented at conferences.

Large-scale functional-connectivity graphical models for individual subjects using population prior. This is a poster for our NIPS work

Multi-subject dictionary learning to segment an atlas of brain spontaneous activity. This is a poster for our IPMI work

Mayavi for 3D visualization of neuroimaging data: powerful scripting and reusable components in Python.

Machine learning for fMRI in Python: inverse inference with scikit-learn.

Anybody reading my blog should have expected me to blog about the end of my GSoC. Sorry to disappoint, but I simply did not experience anything similar to an ending. On the contrary, I feel like things have barely started. Also, I apologize for one of the few posts here without pretty pictures!

For the last two weeks, I’ve been traveling. I attended the EuroScipy conference thanks to Fabian, who offered me a place to sleep during the week. We sprinted hard, we discussed tricky APIs, we drank a lot of coffee, beer, and ate well in lovely Paris. It was great to meet all of the celebrities, the people who keep the scientific Python globe turning.

Many thanks to Gael and Emmanuelle, who worked very, very hard on organizing everything, so they weren’t around and I didn’t get to say goodbye when I ran to catch my plane last Sunday.

I was in a hurry, heading to Tarragona, a beautiful city on the Catalan coast, where the public university organized the 2011 summer school in linguistics and speech technologies (SSLST). This was a great opportunity to meet many fellow young researchers working in computational linguistics. I will not go into details now, because I plan expand on this, but I would like to state a couple of things. Firsty, even though NLP seems to be mostly a Java-dominated affair (note for example Stanford’s NLP toolkit and Sheffield’s GATE), the computational linguistics and psycholinguistics research center (CLiPS) at the University of Antwerp actually briefly manifested its devotion to Python and NLTK via its research director, Walter Daelemans.

It was good to see a little love for Python in this field. NLTK is very underrepresented in the SciPy community, I couldn’t find anybody at the EuroScipy conference knowing too much about it or about the people behind it.

Another lab that has done a lot of cool work is Cental at the Catholic University of Louvain, and they also use Python for natural language processing. Maybe in the coming years, we will see a Python for Computational Linguistics sattelite, along with Physics and Neuroscience. Doesn’t it sound more fun?

Secondly, I wish SSLST were organized by someone like Gael! As the discussion at dinner regarding who will organize next year’s EuroScipy went, it is imperative that the organizers be actively involved in the community, and generally passionate about it. Even though I’m comparing apples and oranges, Carlos Martin-Vide behaved in this context like a old, tired, emotionless academic, not taking into account even lunch breaks for the whole group, not to mention any sort of getting together or even a group photo (which, alas, we were not able to take, apart from small groups.) They said it couldn’t be done. Of course it could, they just didn’t want it hard enough.

Finally, before signing off, I would like to announce that because the Romanian Ministry of Education failed to specify the allocated number of public positions for masters’ programmes, the admission exam at the University of Bucharest will be delayed by a couple of weeks. Luckily, this will allow me to attend RANLP 2011 in Hissar, Bulgaria a week from now, where I will present my poster entitled:
“Can alternations be learned? A machine learning approach to Romanian verb conjugation” by Liviu P. Dinu, Emil Ionescu, Vlad Niculae and Octavia-Maria Sulea. See you in Hissar!

I’ve been working lately in improving the scikit-learn example gallery to show also a small thumbnail of the plotted result. Here is what the gallery looks like now

And the real thing should be already displayed in the development documentation. The next thing is to add a static image to those that don’t generate any result, examples such as the SVM GUI should have an image to display.

Once again, we are looking for a junior developer to work on the scikit-learn. Below is the official job posting. As a personal remark, I would like to stress that this is a unique opportunity to be payed for two years to work on learning and improving the scientific Python toolstack.

Job Description

INRIA is looking to hire a young graduate on a 2-year position to help with the community-driven development of the open source machine learning in Python library, scikit-learn. The scikit-learn is one of the majormajor machine-learning libraries in Python. It aims to be state-of-the-art on mid-size to large datasets by harnessing the power of the scientific Python toolstack.

Speaking French is not a requirement, as it is an international team.

Requirements

• Programming skills in Python and C/C++
• Understanding of quality assurance in software development: test-driven programming, version control, technical documentation.
• Some knowledge of Linux/Unix
• Software design skills
• Knowledge of open-source development and community-driven environments
• Good technical English level
• An experience in statistical learning or a mathematical-oriented mindset is a plus
• We can only hire a young-graduate that has received a masters or equivalent degree at most a year ago.

About INRIA

INRIA is the French computer science research institute. It recognized word-wide as one of the leading research institutions and has a strong expertise in machine learning. You will be working in the Parietal team that makes a heavy use of Python for brain imaging analysis.

Parietal is a small research team (around 10 people) with an excellent
technical knowledge of scientific and numerical computing in Python asas
well as a fine understanding of algorithmic issues in machine learning and statistics. Parietal is committed to investing in scikit-learn.

Working at Parietal is a unique opportunity to improve your skills in machine learning and numerical computing in Python. In addition, working full time on the scikit-learn, a very active open-source project, will give you premium experience of open source community management and collaborative project development.

Contact Info:

• Technical Contact: Bertand Thirion
• E-mail contact: bertrand dotnospam thirion atnospam inria dotnospam fr
• HR Contact: Marie Domingues
• E-mail Contact: marie dotnospam domingues atnospam inria dotnospam fr
• No telecommuting

• http://nicolashoening.de/?blog&nr=114
• http://khinsen.wordpress.com/2011/08/30/euroscipy-2011/
• http://www.visixion.com/?p=608

Motivation for Starting Sage

sage: crt(2, 1, 3, 5)  # Chinese Remainder Theorem

    11 

    sage: crt?        # ? = documentation and examples

    Returns a solution to a Chinese Remainder Theorem...

    ...

    sage: crt??       # ?? = source code

    def crt(...):

    ...

        g, alpha, beta = XGCD(m, n)

        q, r = (b - a).quo_rem(g)

        if r != 0:

            raise ValueError("No solution ...")

        return (a + q*alpha*m) % lcm(m, n)

def python_sum2(n):

    s = int(0)

    for i in xrange(1, n+1):

        s += i*i

    return s

%cython

def cython_sum2(long n):

    cdef long i, s = 0

    for i in range(1, n+1):

        s += i*i

    return s

sage: timeit('python_sum2(2*10^6)')

5 loops, best of 3: 154 ms per loop

sage: timeit('cython_sum2(2*10^6)')

125 loops, best of 3: 1.76 ms per loop

sage: 154/1.76

87.5

sage: var('k, n')

sage: factor(sum(k^2, k, 1, n))

1/6*(n + 1)*(2*n + 1)*n

def sum2(n):

    return n*(2*n+1)*(n+1)/6

%cython

def c_sum2(long n):

    return n*(2*n+1)*(n+1)/6

sage: n = 2*10^6

sage: timeit('sum2(n)')

625 loops, best of 3: 1.41 microseconds per loop

sage: timeit('c_sum2(n)')

625 loops, best of 3: 0.145 microseconds per loop

sage: 1.41/.145

9.72413793103448

sage: c_sum2(2*10^6)   # WARNING: overflow

-407788678951258603

What is Sage?

History

sage: R. = PolynomialRing(ZZ)

    sage: f = x + 1/2; f

    x + 1/2

    sage: parent(f)

    Univariate Polynomial Ring in x over Rational Field

> R := PolynomialRing(IntegerRing());

    > x + 1/2

        ^

    Runtime error in '+': Bad argument types

    Argument types given: RngUPolElt[RngInt], FldRatElt

sage: f(x,y) = sin(x - y) * y * cos(x)

    sage: plot3d(f, (x,-3,3), (y,-3,3), color='red')