This week focused on implementing support for the == and != boolean comparisons for the BSR, CSC, CSR sparse matrix types. The culmination of my work this week is this pull request.

!=

For sparse matrices, the != (not equal to) generally has an output that can be efficiently represented by sparse matrices. With the exception of comparing the whole matrix with a nonzero scalar.

==

The boolean equal to comparison in general has an output that is not efficiently represented by sparse matrices, since all the zero entries return True. The exception being comparison with a nonzero scalar.

Adding a c routine for this operation with the existing code, would be problematic. The binop routines only apply the given binop to element pairs in which one of them is nonzero. So all the co-occurring zero elements would not return True like they should. Because of this, and == being in general inefficient, I did not bother to implement the == operation in sparsetools. Instead, the == operation is computed using the != operation and inverting it.

Implementation

C++ & SWIG level

Implementing the routines for != and == in the sparsetools/ was pretty easy, using the existing code. There are handy csr_binop_csr and bsr_binop_bsr functions which take c++ a functional as an argument for the binop. The only problem is that the output type is not the same as the input type. So the output is not bool by default.

Python level

In Python, != and == operations are can be overridden by defining the methods __ne__ and __eq__. I did this to handle the case wher the other object is a scalar, dense, or a sparse matrix.

I defined these methods in compressed.py which defines the base class for BSR, CSC, and CSR.

In defined these methods in base.py, here they just convert the sparse matrix to CSR and preform the desired operation.

__bool__

There was however some complication with how these worked when a dense ndarray is the first argument, ideally it would ask the sparse matrix what to do, however this was not happening whenever the sparse matrix's __bool__ (for Python3, __nonzero__ for python2.x) method returned True or False, I still don't know why. It did however work when __bool__ raised an error or returned something other than T/F. Since NumPy's ndarrays raise an error when __bool__ is called, I thought this was a reasonable change to make to sparse matrices __bool__. Although ideally, I'd like to better understand what is going on so I don't have to make an incompatible change to SciPy's API.

This change broke one test in scipy/io/matlab/tests/test_mio.py there could be more serious issues relating to this change, currently they are being discussed on the pull request

IPCluster provides an easy way to do parallel computing in Python, allowing you to connect to and control Python processes running on multiple compute cores, either locally or across a cluster. This introduction should take about 20 minutes and will show you how to start a cluster, access engines, and execute commands. We will use our web-based Wakari analytics environment that will allow you to get started right away, without needing to install or configure anything.

I'll be in Austin, TX for SciPy next week, and I'll be giving an updated talk about tools for reproducible research.

What follows is the abstract of a talk that I really wanted to give and submitted to the general track which did not get in. I care deeply about this topic, and I'm hoping to have conversations with others about it. The talk was deferred to a poster, so if you don't bump into me at other times, seek me out during the poster session (mine's poster #15, 10:35 AM - 11:35 AM on June 27th). This post is intended as a pre-conference warm up to that.

Navigating the Scientific Python Communities - the missing guide.

The awful truth: our newcomers have a deluge of options to wade through as they begin their journey. What tools are available, how do I install them, how do I make them work together -- all hard questions facing a budding scipythonista.

On developer-friendly platforms, the popular approach is to just install numpy, scipy, matplotlib, and ipython using package management facilities provided by the operating system. On Mac OS X and Windows, the least painful, bootstrapping approach is to use a Python distribution like Python(X,Y), EPD, Anaconda, or Sage.

Both of these paths obscure a reality which must be stated explicitly: the development of packages is fundamentally decentralized. The scientific python ecosystem consists of a loose but thriving confederation of projects and communities.

The chaos of installation options and lack of centralization around a canonical solution, which on the surface appears to be a point of weakness that is detrimental to community growth, is a symptom of one of its greatest strengths. Namely, what we have is a direct democracy for user-developers. They vote with their feet: filing bug reports, testing pre-release versions, participating in mailing lists, writing and reviewing pull requests, and advertising the tools they use in talks, papers, and conversations at conferences.

I will discuss why this is the case, why this is a good thing, and how we can embrace it. The talk will present the ethos and expectations of an effective newcomer, as well as resources and strategies for incremental progress toward becoming a master scipythonista.

I had the opportunity from Packt Publishing to review the second edition of Numpy Beginner’s Guide. Many thanks to the publisher for this and let’s go to the review.

Content and opinions

The first chapter introduces the reader to the scientific ecosystem. The main modules are covered, but there are some small mistakes (ipython –pylab does not import Matplotlib, Numpy ans Scipy, just a subset of the first that gives a Matlab-like interface aith some renamed Numpy functions) and what I think is a bad habit (importing everything from Pylab and using it as is). Nothing is lost because in the second chapter, Numpy is properly introduced and imported explicitly. There is a link missing between the two chapters because I didn’t understand why Numpy was used this way and not with pylab.

So the second chapter is about the core Numpy data structure: the multidimensional array. The author browses through different ways of creating them, by stacking them, flattening them… The next chapter deals with universal functions, or put it in an other way, functions that run on array element by elements. There are many different way to do so, with specifics that are tackled properly.

Before the chapter of matrices, there is a useful chapter on correlation, convolution and polynomials. Although one may want to go up to Scipy for signal processing, more often than not, Numpy is enough. There enough examples to understand how everything work. Then, the much hated matrix class is introduced. There are many discussions online on whether this class is actually useful, and I won’t delve on this. Suffice to say that the power of this class can be seen in the examples.

The following chapter is about the different submodules inside Numpy: linear algebra, fft, random number. The proper pointers to Scipy are provided, as well as explicit samples. I only can regret the time shifting example is not perfect, as a filter is applied on the amplitude and the inverse FFT is applied on the amplitudes only. So the transformed signal loses all the phase information, which may be why the image is similar but not that similar to the original one.

The seventh chapter tackles different extra functions, mainly finance-related, and I have to say that I don’t know their use enough to comment. Of course, this is why they are in the middle of the book and not in the first pages as the other Numpy functions. Still, there are some useful functions here (some I didn’t know about), as sorting, searching…

When one code scientific code, one often forget about testing. And Numpy has a nice module for scientific testing. It is nice to know that this aspect of science is not forgotten here and has a proper introduction. (also don’t forget about version control!)

The last three chapters introduces additional modules. The first one was partial addressed before, Matplotlib, and if you need something more advanced, there is always the Matplotlib book also published by Packt. Then there are some examples with Scipy, Scikits (soon a new book on Machine Learning with scikit-learn will be available, also by Packt, for which I was a technical reviewer, and it is really great) and other tools. The final chapter is about Pygame, but I don’t code games 100% in Python, so I didn’t really read it!

Conclusion

It’s hard to be mad at the author for the import issue. But I find it also difficult not to, as the philosophy changes depending on the chapter without saying why. Still, for an introductory book to Numpy, it is great if not excellent. A lot of simple examples, a lot of checks, the good pointers to write efficient code…

Major Changes

• A dozen new plugins were added - mahotas, cffi, lxml, paramiko, Babel, PycURL and others. 
• Eliminated most duplicate dependencies between packages - HDF5, zlib, bzip2, libpng, jpeg, tiff etc. These are part of their respective python wrappers or grouped under base_libraries.
• base_python plugin created to host numerous infrastructure, python3 compatibility and utility packages (six, decorator, curses, pep8 etc).
• Image manipulation was greatly improved with the addition of FreeImage python wrappers and replacing PIL with pillow.
• WindowsXP is no longer officially supported. HDF5 cannot be compiled to be thread safe on it.

Notable Plugin Changes:

• Sypder - updated to 2.2.0
• ETS updated to 4.3.0.
• VTK - updated to 5.10.1, Enabled: TextAnalysis, QVTK, WrapSIP, Oggrtheora,
• ITK - updated to v4.4.0, Enabled: VTKGlue.
• NetCDF4 - Enabled support for DAP and HDF4.
• lxml - libxml & libxslt rebuilt with iconv, ICU, zlib and LZMA.
• Numpy, Scipy and IPython updated.

From weather readings to financial tickers, many interesting data come in the form of streams. An issue with streaming data analysis in IPython is that one often has to either be collecting the data from the stream, or analyzing it, but not both. In this tutorial, I demonstrate a single-notebook IPython workflow for simultaneous collection and analysis of realtime streaming data:

There are a few things that have taken some getting used to in working with SciPy. These things are aggregated here:

Keeping my SciPy repo up-to-date

I need to upadate my copy of SciPy regularly from upstream, to avoid merge conflicts. I have done this plenty of times, but for some reason I forget how every time. There is a good howto on stackoverflow on this. From said article:

# Add the remote, call it "upstream":

git remote add upstream git://github.com/whoever/whatever.git

# Fetch all the branches of that remote into remote-tracking branches,
# such as upstream/master:

git fetch upstream

# Make sure that you're on your master branch:

git checkout master

# Rewrite your master branch so that any commits of yours that
# aren't already in upstream/master are replayed on top of that
# other branch:

git rebase upstream/master

Pushing new local branches to a remote Git repo

Every time a make a new local branch, I forget how to push it to github as new branch. Google always brings me to this article. Basically if I have a new local branch call it new-feature and want to push it to my github repo and create a new new-feature branch there too. I do:

git push -u origin new-feature

and if I wanted to delete this remote branch, I could do:

git push origin :new-feature

and delete it locally with:

git branch -D new-feature

SciPy's commit messages

SciPy require commit messages start with some standard acronyms.

API: an (incompatible) API change
BLD: change related to building numpy
BUG: bug fix
DEP: deprecate something, or remove a deprecated object
DEV: development tool or utility
DOC: documentation
ENH: enhancement
MAINT: maintenance commit (refactoring, typos, etc.)
REV: revert an earlier commit
STY: style fix (whitespace, PEP8)
TST: addition or modification of tests
REL: related to releasing numpy

The acronyms are easy to remember, but adding the acronyms is easy to forget. Which brings me to my next topic.

Changing commit messages, after you have pushed

This is not the simplest thing to do in git, fortunately my mentor Pauli showed me a great way to do it. Assuming no one has pulled from your repo yet.

• First do git log and look for the last COMMIT befor the commit messages you want to change.

• Do git rebase -i COMMIT

• Replace pick with r

• Then Git will prompt you to change the commit messages one at a time.

• Push your changes with git push -f

I got accepted to Google Summer of Code! Thank you everyone who has helped me get this far, and Amber for making me a congratulatory google cake!

On my proposal timeline I gave myself 2 weeks to add support to sparse matrices for dtype=bool. So far, dtype=bool works implementation was not too hard, and most of my time was spent figuring out how SWIG works. However there is still a lot of work to do in testing bool support. This will be more time consuming than it should be because the sparse test suite is messy.

bool dtype implementation

sparsetools

The sparsetools/ directory contains several *.cxx, *.h, *.i, and *.py files. The header files (*.h) contain C++ routines these are wrapped by SWIG according to instructions in the *.i interface files. SWIG uses the header and interface files to generate python *.py files C++ code in *.cxx files. Once these *.cxx files are compiled and linked the functions they define can be called by the generated *.py files.

Where do I add bool support?

Contained in sparsetools/ is a file called complex_ops.h which translates numpy's complex types to C++ classes. Well, that is what I want to do, except for numpy's bool types instead! My new file bool_ops.h is much simpler though. It basically does

typedef npy_int8 npy_bool_wrapper;

Yep that's it, (for now anyway). This file basically says to treat the npy_bool_wrapper type like the npy_int8. Because boolean algebra is almost integer algebra.

npy_bool_wrapper is used in numpy.i to define a typemap from numpy's npy_bool type to our npy_bool_wrapper.

In sparsetools.i we have to declare the new data type and add it to the INSTANTIATE_ALL macro which is used in all the *.i files to call the functions we want from the *.h files.

What else?

Adding bool to supported_dtypes in sputils.py prevents it from upcast. Then tests need to be added, but this isn't too easy. The test suite does not have a way to take some general data or dtype. So I am going through adding tests one-by-one where it seems appropriate. I'm also trying to add tests for other dtypes because there don't seem to be any tests for anything other than float64 and int64.

While I was doing this I discovered a few bugs relating to how sparse handles uint dtypes.

Once I'm done adding tests and fixing what ever does not pass, I'll move on to the next stage in my proposal. Adding support for boolean operations.

One of IOPro’s more advanced features is the ability to work with csv data directly from Amazon’s S3 service. Combining this with Wakari’s notebook sharing feature gives us a powerful way to collaborate on data analysis using data stored in the cloud.

I followed a Google+ link to an article that has just come out in Science magazine. The article is "Troubling Trends in Scientific Software Use" by Lucas Joppa and others ¹.

I did not find Joppa's article very satisfying, but there were some interesting references. For example (from Joppa et al):

From fields such as high-performance computing, we learn key insights and best practices for how to develop, standardize, and implement software (11)

Reference 11 is an article by Victor Basili and others on "Understanding the High-Performance-Computing Community" ². It is a thoughtful reflection on the experience of academic software engineers working with scientists using massively parallel computing. The message I took from the article was that scientists may have good reasons to reject suggestions from academic software engineers. This is often because the solutions are too general to be useful for a particular problem. The last sentence from the article is "We've much to learn from each other".

Another couple of sentences in the Joppa article led me to interesting places:

Reliance on personal recommendations and trust is a strategy with risks to science and scientist. "End-user developers" commonly create scientific software (17, 21, 22) ¹

Reference (22) is "A Software Chasm: Software Engineering and Scientific Computing" by Diane Kelley ³. Kelley also has some interesting thoughts on the separation of science and software engineering:

Creating the chasm

In a 1997 paper, Iris Vessey described a shift in computing philosophy that occurred in the late 1960s. She quoted George E. Forsyth, the ACM's president at that time, as stating that computer science was a field of its own, separate from the application domains. The result, as Vessey points out, was the insistence that anyone who wanted to be taken seriously in the field of computing, and software engineering, must develop domain-independent techniques, methods, and paradigms. The pressure was such that claims of broad applicability became commonplace. So, the bridge between software engineering and any application domain became a single massive structure that everyone had to use.

This is the relevant quote from Vessey (1997)⁴:

Then followed a debate, which endured for over a decade, that pitted academia against the computing profession. The notions of traditional computer scientists steeped in the notion of theory and the stigma of all things "applied" led to the following statement by George E. Forsyth, President of ACM, in 1965: "... the core of computer science has become and will remain a field of its own, ahead of, and separate from the application domain specialists." This philosophy led, not only to the creation of generic "solutions" to problems, but also to claims that any "creation" was applicable in all situations, because to state otherwise would have branded the creator as lacking in academic respectability.

I wonder if this separation has had a larger effect on our thinking than we recognize.

For example, our original NIPY grant application asked for salaries for two full-time programmers to help us apply software engineering methods to writing neuroimaging software. We did hire two programmers but we had given them an impossible job. To do a good job, they had to now become experts in a scientific field they did not know. Useful code and documentation had to be written so that other scientists in the field would find it easy to follow. To do this, you need to know the field. We scientists spend more time than we think working out how to talk to our colleagues. That experience is very hard to teach.

These articles gave me a new way to think about the false separation of "software" and "science". It is mysterious that we make this separation, because my whole experience tells me they are closely linked. I wonder whether we make this separation because we have come under the unconscious influence of the movement that Kelley and Vessey describe. It is easy to see how tempting it is. How easy life would be if we really could pay a programmer to write up our ideas as we sketch them airily on a white-board and retire to the garden, to think great thoughts.

We various of the Berkeley NIPY crowd were in a Mexican restaurant yesterday lunchtime. We discussed the IPython notebook and scientific programming and reproducible research.

One question struck me as fundamental: do scientists have to learn to be software engineers?

I think the common answer to this is "No - scientific coding is different". The "No" camp points out that scientists use code in a different way to software engineers. Science is exploration, and so is scientific coding. Scientific code is provisional, often changing. There is no specfication, there is no production system that must just work. We are trying stuff out, seeing what works, adjusting to our better understanding of the data and the ideas.

Maybe the scientist's idea of the canonical software engineer is someone who writes and hosts a web application with some database behind it.

Accepting the "No", we look at things that software engineers should learn:

1. Version control. Essential for a software engineer, apparently, because they need to be able to rollback their changes, and keep track of released and development versions. Desirable maybe for a scientist, but not essential, because the code never goes into production, and we always use the latest version of the code.
2. Testing. Essential for a software engineer, otherwise they may release code that corrupts a database of financial transactions or delivers the wrong item to your address. Desirable maybe for a scientist, but not essential because the nature of the code changes so often that it would only slow things down to write exhaustive tests, only for the code to change track and make the tests obsolete.
3. Code review. Essential for a software engineer, because they work in teams with other software engineers. This is their job, to write code, and they can learn from each other. Desirable maybe for a scientist, but not essential because each scientist owns their own problem and so only they can understand their code. Explaining the code is time-consuming and slows progress in developing new ideas. Others in the lab are not software engineers, so they are not trained to review code, they will not enjoy it and they will not do a good job.
4. Releases. Essential for a software engineer because they may be paid to provide code that others can use. The release labels code that can be used. The software engineer can and should make sure the code works as advertised. Desirable maybe for a scientist, but not essential, because the code changes often, and the code may be written to solve a very specific problem that could not be of much interest to another researcher. If the other researcher is interested, they should make the effort to work out how to use the code, that is not the job of the author-scientist, they are not paid to support software, but to produce science.
5. Documentation. Essential for a software engineer because they want their code to be used correctly by other people. Desirable but not essential for the scientist because scientific code is not for distribution except to other scientists who must take responsibility for the code themselves.

I think that these set of beliefs lie at the heart of the problem for reproducible science.

The beliefs rest on the following model of writing scientific code:

A folk model of scientific code

A scientist named A writes some code. They check the code. They confirm it is working. Most likely this means the code is free of major errors.

Another scientist named B is reading some results published by A. They can assume that the results are free of major errors.

Living with the folk model

Now everything makes sense. Version control, testing, code review, releases, documentation are great, of course, but if we can assume there are no major errors, then - why would B want to see my code? If B gets my code, why would she need to understand it? Why would she need to know what version it was, or how that related to my paper? She should assume that there are no major errors.

Rejecting the folk model

The folk model is not just too simple, it is flat-out wrong. We make mistakes all the time. All the time. If you know that, then you know that B needs to check your stuff. They can't do science without checking your stuff. If they need to check your stuff, you need to write your code for production. Then, nothing is optional. We all need; version control, testing, code review, releases, documentation.

I was talking to an esteemed colleague about six months ago about how easy it was to believe that we do not make errors, unless we check.

He said that he had noticed this change when his students and research assistants started to use their own computer programs to analyze data. Before this, when he and others added up a column of numbers, they would check several times that they were correct. After, the student or RA would analyze some data with a program they had written themselves, and my friend would say "are you sure the program is right" and they would say "yes I'm sure". My friend would then go through the results and check; they would often find mistakes. He thought that there was something about using a computer that made it easy to believe that the result was right, even if the process of getting the result had the same chance of error as a simpler task like adding numbers.

To quote David Donoho:

In my own experience, error is ubiquitous in scientific computing, and one needs to work very diligently and energetically to eliminate it. One needs a very clear idea of what has been done in order to know where to look for likely sources of error. I often cannot really be sure what a student or colleague has done from his/her own presentation, and in fact often his/her description does not agree with my own understanding of what has been done, once I look carefully at the scripts. Actually, I find that researchers quite generally forget what they have done and misrepresent their computations.

Computing results are now being presented in a very loose, “breezy” way—in journal articles, in conferences, and in books. All too often one simply takes computations at face value. This is spectacularly against the evidence of my own experience. I would much rather that at talks and in referee reports, the possibility of such error were seriously examined.

-- David L. Donoho (2010). An invitation to reproducible computational research. Biostatistics Volume 11, Issue 3 Pp. 385-388

There is something about computational science that seems to make the idea of error less worrying or important:

In stark contrast to the sciences relying on deduction or empiricism, computational science is far less visibly concerned with the ubiquity of error. At conferences and in publications, it’s now completely acceptable for a researcher to simply say, “here is what I did, and here are my results.” Presenters devote almost no time to explaining why the audience should believe that they found and corrected errors in their computations. The presentation’s core isn’t about the struggle to root out error — as it would be in mature fields — but is instead a sales pitch: an enthusiastic presentation of ideas and a breezy demo of an implementation. Computational science has nothing like the elaborate mechanisms of formal proof in mathematics or meta-analysis in empirical science. Many users of scientific computing aren’t even trying to follow a systematic, rigorous discipline that would in principle allow others to verify the claims they make. How dare we imagine that computational science, as routinely practiced, is reliable!

-- David L. Donoho, Arian Maleki, Inam Ur Rahman, Morteza Shahram and Victoria Stodden (2009) Reproducible Research in Computational Harmonic Analysis. Computing in Science and Engineering 11(1) pp 8-18

So, following up on my transition to pelican, here's a little help for navigating a pelican site. If you've ever wondered why sometimes there's a highlighted menu item at the top, and other times there isn't, read on.

Categories

With Pelican, every post belongs to exactly one category - that category is 'misc' if one is not provided (and you can override that named by setting DEFAULT_CATEGORY = 'somethingelse' in your pelicanconf.py file). Alternatively, if the post is in a subdirectory, it's category becomes the subdirectory name. So if you have a situation like this:

.../content/unicorns/lady-rainicorn.md

The lady-rainicorn post will be filed under the category of unicorns. You can explicitly override this category-from-subdirectory behavior, by specifying the category in the header of the file, like so:

title: All about Lady Rainicorn
slug: lady-rainicorn
category: adventure-time

Categories are where a lot of the pelican themes get the top navigation bar, which is referred to as the menu in Pelican. You'll notice that this post is filed under technology. If you're viewing this on the main page, that won't be highlighted, but if you click on the title of this post, you'll see that technology will become highlighted. Same deal if you just click on technology, which will get you a listing of the latest post in that category, and a pagination of all previous ones.

Pages

The menu items need not be limited to categories, however. With the popular notmyidea template (which is what I've customized for my journal here), if you create any pages (by putting the file in a content/pages/ directory), then the title of that page will appear in front of all of the categories in the menu. These pages don't show up in RSS feeds, though.

Menuitems

Additionally, you can specify other MENUITEMS in pelicanconf.py, which will go in front of all pages and categories. MENUITEMS is a list of tuples, each tuple should be composed of a link name, and a url. This is how I link to the listing of all of my blog posts:

MENUITEMS = [('all', '/blog/archives.html')]

Tags

Finally, you can have zero or more tags attached to a post. Tags, by default, do not get a feed generated for them. To enable a tag feed, you'll have to set a line like this in your pelicanconf.py files:

TAG_FEED_ATOM = "feeds/tag_%s.atom.xml"

Menu highlighting

The notmyidea template only contains logic for highlighting the active category. This is why, if you click on 'All' at the top of mine, it doesn't stay highlighted when you're on that page.

On the other hand, if you click at cycling at the menu bar, it will stay highlighted, and you will see a listing of the articles in that category. But if you click on the cycling tag in one of those posts - the menu item will not be highlighted, and you will see some more articles - in particular, right now my review of Just Ride will also be listed, but it lives in the category of books, because that's where I'm keeping book reviews.

Summary

So there you have it - now you know the differences between categories, tags, pages, and menu items, as far as Pelican is concerned. For me, the bottom line is that I'll probably stick to just using categories, and will try to use tags sparingly. In my transition from WordPress, I cleaned up and removed a whole bunch of tags (since they were either redundant (the same tags kept appearing together), or were too specific (most tags were only used for one post).

Well, it's happened again, I've jumped (back) on the static blog engine bandwagon. Early versions of my site were generated literally using #define ,#include, gcc, and a Makefile...). Back then I was transitioning away from using livejournal, and decided to use WordPress so I wouldn't have to roll my own RSS feed generator.

I tolerated WP for a while - and the frequency of my posts was so low, that it wasn't much of an issue.

Except for the security upgrades. I logged into the wp-admin console way more times than I cared to just to press the little "upgrade" box. The reason is that wordpress keeps everything in a database that gets queried every time someone hits the site. I was never comfortable with the fact, because content can be lost in case of database corruption during either an upgrade or a security breech. Also, my content just isn't that dynamic. The WP-cache stuff just seemed overkill, since I don't get that many visitors.

But lately, I've found myself wanting to write more, to post more, but also shying away from it because I hate dealing with the WordPress editor. And I also hate being uncertain about whether any of it will survive the next upgrade, or the next security hole, whichever I happen to stub my toe on first.

And the thing is, I really like to use version control for everything I do. I liked my blog posts to be just text files I can check into version control. I also like typing "make" to generate the blog, and now I get to!

For added fun, I'm hoping that writing my posts in markdown will make it easier to coordinate my gopher presence, since it's pretty close.

For posterity, I'm capturing what the first version of my Pelican-based blog looked like. I did the same thing when I moved to WordPress.

I ran into some confusing things about transitioning to Pelican, so I thought I'd note them here, for the benefit of others.

Unadulturated code blocks

I like to use indentation as a proxy for the venerable <tt> tag - which uses a monospace font.

If you just want to use an indentation, but do not want the indented text parsed as a programming language, put a :::text at the top of that block. Here's what I mean. Take this Oscar Wilde quote, where I've inserted a line break for drammatic effect, for example:

A gentleman is one who never hurts anyone's feelings
unintentionally.

Now, you see that ugly red box around the apostrophe? Well, That's because all I did was indent the two lines. If I just put a :::text above the quote, indented to the same level,

    :::text
    A gentleman is one who never hurts anyone's feelings
    unintentionally.

the result will render like this.

A gentleman is one who never hurts anyone's feelings
unintentionally.

As the pelican documentation specifies, this is also the way you can also specify the specific programming language you want, so :::python would be one way to not make pygments guess. You can get a list of all supported languages here, just use one of the short names for your language of choice.

And you should really do this, even if you aren't bothered by the red marks, because the code highlighting plugin goes in and tokenizes all of those words and surrounds them in <span> tags. Here's the HTML generated for the first version:

<div class="codehilite"><pre><span class="n">A</span> <span class="n">gentleman</span> <span class="n">is</span> <span class="n">one</span> <span class="n">who</span> <span class="n">never</span> <span class="n">hurts</span> <span class="n">anyone</span><span class="err">&#39;</span><span class="n">s</span> <span class="n">feelings</span>
<span class="n">unintentionally</span><span class="p">.</span>
</pre></div>

and here's the version generated when you add the :::text line at the top:

<div class="codehilite"><pre>A gentleman is one who never hurts anyone&#39;s feelings
unintentionally.
</pre></div>

nikola?

At some point, having a dialogue with myself, I wrote in here "or should I not use pelican and use nicola instead?"

Ok, tried it - nikola takes too long to startup -

nikola --help
real    0m1.202s
user    0m0.876s
sys     0m0.300s

pelican --help
real    0m0.639s
user    0m0.496s
sys 0m0.132s

I'm sure it's a fine static blogging engine - and Damian Avila's already written IPython Notebook converters for it, but it just feels like it tries to do too many things. Constraints are good. I'll stick with Pelican for now. (Though I did use the nikola wordpress import tool to grab the wp-upload images from my WordPress blog)

another set of instructions I consulted: Kevin Deldycke's WordPress to Pelican, which is how I did get my articles out using exitwp which I patched slightly, so files got saved as .md, and preserve other format properties.

Redirects

also known as: not breaking the web

I wanted to preserve rss feeds, and also not break old WordPress style /YYYY/MM/DD urls - the Nikola wp-import script had created a url remapping scheme in a file called url_map.csv.

I don't have that many posts, so I just added them in by hand:

Options +FollowSymLinks
RewriteEngine on
RedirectMatch 301 /blog/feed/ /blog/feeds/all.atom.xml
RedirectMatch 301 /blog/2006/10/18/todd-chritien-greens-choice-voting/ /blog/todd-chritien-greens-choice-voting.html
RedirectMatch 301 /blog/2007/01/04/changelogs-with-dates-gui-goodness/ /blog/changelogs-with-dates-gui-goodness.html
...

Enabling table of contents for posts

If you want to include a table of contents within a post using [TOC], you must enabled the markdown toc processor with a line like this is your pelicanconf.py:

MD_EXTENSIONS =  [ 'toc', 'codehilite','extra']

Categories and tags

Ok, so this was never clear to me in wordpress, either - but what's the difference between a tag and a category? is it the case that a post can only belong to one category, whereas it can have any number of tags?

I think I used categories as tags on wordpress. Looks like all posts on Pelican can have at most one category. Turns out this little aside was long enough to turn into its own post, so if you're interested, pelican tags-vs-categories has got you covered.

That's it for now

Thanks to Preston Holmes (@ptone) for encouraging me to transition away from WordPress, and pointing me to this post by Gabe Weatherhead (@MacDrifter) for how to do that. It should be said that the pelican documentation itself is very good for getting you going. Additionally, I consulted this post by Steve George which has a good description to get you started, and also covers a bunch of little gotchas, and lots of pointers. Also, thanks to Jake Vanderplas (@jakevdp) for his writeup on transitioning to Pelican, which I will consult later for incorporating IPython notebooks into my markdown posts, in the future. This is good enough for now. LTS.

Welcome. This post is part of a series of Continuum Analytics Open Notebooks showcasing our projects, products, and services.

In this Continuum Open Notebook, you’ll learn more about how Numba works and how it reduces your programming effort, and see that it achieves comparable performance to C and Cython over a range of benchmarks.

A page about cycling in the East Bay (and the San Francisco Bay Area in general).

This started off as the sort of resource that I wish I had when visiting a new town: so hi there, out-of-towner, welcome to our corner of the world. Hope you enjoy your stay and find these resources helpful.

You need a bike, a water bottle, a jersey or a saddle bag (to store your snack, phone, patch kit). Everything else is optional.

Routes

Though I consider myself an avid cyclist, the reality is that I am only mildly so. I have not been on a ton of rides in the area, and don't know them well at all outside of the East Bay - I have yet to cycle across the Golden Gate Bridge - and the Marin Headlands (an omission I plan to rectify soon). I did mountain bike a bit in Marin when I was a kid, and rode a fair bit around the Peninsula while I was in high school. I got my first roadie the summer before heading off to UC Davis. It was a wonderfully light steel blue Peugeot from the 70s (you can see a picture of it here). I convinced (my then future roommate) Philip to get a matching white one. But it was taken from me my first year in grad school (2006). I hadn't been cycling much since, until I started up again recently. Below, you will find a concise guide to some of my favorite rides in the area. As far as my actual experiences, you can get a flavor of the rides I go on by reading my cycling log.

20-40 km

Here is the overview: If you want a shorter ride with a great view - I recommend this 24 km (15 mi) loop from UC Berkeley to Grizzly Peak. The ascent to Grizzly Peak Rd is steady and manageable. The road does get pretty rough when you turn off of College until you get to Broadway, but you avoid the dangerous Ashby route to Tunnel Road. Recently, I've been doing a similar ride that avoids that rough patch by going down College Ave all the way to Broadway. As an alternative option, instead of going up to Grizzly Peak, you can turn off for a slightly shorter (but really fun) descent that ends at the historic Claremont Hotel.

What's nice about these loops is that you can do them back-to-back for a total of ~50 km. But If you want to do a ride that long with more variety, read on for the longer routes.

50 km and beyond

If you have some time (about twice the time it takes to do the shorter routes above), the Redwood-Pinehurst route will take you through an incredible canyon that's filled with redwoods! You can take the loop in the other direction - Pinehurst-Redwood.

Another local favorite is the "Three Bears" ride, which has some truly scenic views all throughout (62 km, ~1 km gain).

The longest ride I've done so far was pretty pleasant, and had a wide variety of terrain, see the map (93 km, ~1.5 km gain).

Centuries

I haven't yet completed a century. But there are some local ones that I have my eye on.

Grizzly Peak Century is one that goes through some of the local East Bay routes I described above.

Best of The Bay Century likewise starts off in the East Bay routes I'm most familiar with, and heads south all the way to San Jose, before returning a bit to Fremont BART station.

Foxy's Fall Century is an organized ride based out of Davis, CA, where I was an undergrad.

Davis Double Century is in my sights as well. Maybe next year I'll get on the board of the California Triple Crown.

There are a lot of long rides listed on the Different Spokes San Francisco website. San Francisco Randonneurs run what looks like a very regular brevet series.

Routes to avoid

Ohlone Greenway: it's not fun - there are too many stops and unprotected crossings. This is a nice route to take a stroll through, on a bike it's just frustrating.

Local Bike Shops

I want to keep track of reasonable places that rent bicycles in the area (in particularly road bikes), as well as routes you should be sure to checkout (with variations, caveats, and general overview).

Missing Link Bicycle Cooperative - One can't say enough good things about this place, super friendly atmosphere. They're located right downtown Berkeley, they both sell gear and fix bikes, and even have two complete sets of tools and stands you can use to fix up your bike for free!

Blue Heron Bikes - a relatively new shop (opened in 2012), in "Westbrae" right of the Ohlone Greenway, I cycle past it every day on my commute. Soon after they first opened, I chatted with the owner, Rob, and he's a great guy.

Wheels of Justice Cyclery - The wheels of justice grind slowly, but exceedingly fine.

Resources

East Bay Bike Coalition - a very active local group.

Albany Strollers and Rollers - Albany, California, is a little 1.5 square mile city just north of Berkeley which I call home.

Berkeley Critical Mass - 2nd Friday of each month - departs at 6pm from Downtown Berkeley BART station.

East Bay Bike Party - a monthly evening ride.

Just Ride - my review of Grant Petersen's book.

Here's some ASCII art I made that I now use as my email signature:

                   _
                  / \
                A*   \^   -
             ,./   _.`\\ / \
            / ,--.S    \/   \
           /  `"~,_     \    \
     __o           ?
   _ \<,_         /:\
--(_)/-(_)----.../ | \
--------------.......J

Continuum Analytics, the premier provider of Python-based data analytics solutions and services, today announced the release of Wakari version 1.0, an easy-to-use, cloud-based, collaborative Python environment for analyzing, exploring and visualizing large data sets

Forgot to mention i am in San Francisco at the moment with Work just took this picture for the lulz:

So this is kind of a rant/political type post. When i started my project GCCPY as soon as it was announced as part of Google Summer of code 2010. A person who i should probably leave nameless nearly immediately jumped on #gcc irc.oftc.net, and started super questioning me on why the hell i would approach python in this way. To the point where he said that i should not bother and just contribute to PyPy. And the moderators in the channel had to kind of step in because it was just way too much.

So ok fair enough, then everything was fairly quiet i get the odd email from interested people then i noticed my project was on reddit http://www.reddit.com/r/Python/comments/1dc0tl/python_frontend_to_gcc_xpost_from_rgcc/

Same guy on there generally a negative slate on my project again. Partly i don’t blame them since my wiki is not that detailed. But i tend to feel people don’t think for themselves and just start believing people from established projects are rock stars, to the point were new projects are pointless to them.

Everyone is entiled to their opinion but what i will say is, think for yourself. And _dont_ compare a project like GCCPY to PyPy. Event just look at PyPy’s main page they have over 25k in donations!! What do i have, my spare time at work and doing everything from scratch in C where as pypy is basically python and jit reusing libpython.so. So before you start slating something as not worth its and cant do x,y,z. I have my own opinion that i believe this is a valid way of implementing python. Yes its _no_ where near python 2 complete! Make it a full time job it would be done in a year. I will say from my benchmarks which i want to post at the end of the year at python con IE; gccpy is looking truly amazing.

Think for yourself reddit!

Following a challenge proposed by Gael to my group I compared several implementations of Logistic Regression. The task was to implement a Logistic Regression model using standard optimization tools from scipy.optimize and compare them against state of the art implementations such as LIBLINEAR.

In this blog post I'll write down all the implementation details of this model, in the hope that not only the conclusions but also the process would be useful for future comparisons and benchmarks.

Function evaluation

The loss function for the $\ell_2$-regularized logistic regression, i.e. the function to be minimized is

$$ \mathcal{L}(w, \lambda, X, y) = - \sum_{i=1}^n \log(\phi(y_i w^T X_i)) + \lambda w^T w $$

where $\phi(t) = 1. / (1 + \exp(-t))$ is the logistic function, $\lambda w^T w$ is the regularization term and $X, y$ is the input data, with $X \in \mathbb{R}^{n \times p}$ and $y \in \{-1, 1\}^n$. However, this formulation is not great from a practical standpoint. Even for not unlikely values of $t$ such as $t = -100$, $\exp(100)$ will overflow, assigning the loss an (erroneous) value of $+\infty$. For this reason ¹, we evaluate $\log(\phi(t))$ as

$$ \log(\phi(t)) = \begin{cases} - \log(1 + \exp(-t)) \text{ if } t > 0 \\ t - \log(1 + \exp(t)) \text{ if } t \leq 0\\ \end{cases} $$

The gradient of the loss function is given by

$$ \nabla_w \mathcal{L} = \sum_{i=1}^n y_i X_i (\phi(y_i w^T X_i) - 1) + \lambda w $$

Similarly, the logistic function $\phi$ used here can be computed in a more stable way using the formula

$$ \phi(t) = \begin{cases} 1 / (1 + \exp(-t)) \text{ if } t > 0 \\ \exp(t) / (1 + \exp(t)) \text{ if } t \leq 0\\ \end{cases} $$

Finally, we will also need the Hessian for some second-order methods, which is given by

$$ \nabla_w ^2 \mathcal{L} = X^T D X + \lambda I $$

where $I$ is the identity matrix and $D$ is a diagonal matrix given by $D_{ii} = \phi(y_i w^T X_i)(1 - \phi(y_i w^T X_i))$.

In Python, these function can be written as

def phi(t):
    # logistic function, returns 1 / (1 + exp(-t))
    idx = t > 0
    out = np.empty(t.size, dtype=np.float)
    out[idx] = 1. / (1 + np.exp(-t[idx]))
    exp_t = np.exp(t[~idx])
    out[~idx] = exp_t / (1. + exp_t)
    return out

def loss(w, X, y, alpha):
    # logistic loss function, returns Sum{-log(phi(t))}
    z = X.dot(w)
    yz = y * z
    idx = yz > 0
    out = np.zeros_like(yz)
    out[idx] = np.log(1 + np.exp(-yz[idx]))
    out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx])))
    out = out.sum() + .5 * alpha * w.dot(w)
    return out

def gradient(w, X, y, alpha):
    # gradient of the logistic loss
    z = X.dot(w)
    z = phi(y * z)
    z0 = (z - 1) * y
    grad = X.T.dot(z0) + alpha * w
    return grad

def hessian(s, w, X, y, alpha):
    # hessian of the logistic loss
    z = X.dot(w)
    z = phi(y * z)
    d = z * (1 - z)
    wa = d * X.dot(s)
    Hs = X.T.dot(wa)
    out = Hs + alpha * s
    return  out

Benchmarks

I tried several methods to estimate this $\ell_2$-regularized logistic regression. There is one first-order method (that is, it only makes use of the gradient and not of the Hessian), Conjugate Gradient whereas all the others are Quasi-Newton methods. The method I tested are:

• CG = Conjugate Gradient as implemented in scipy.optimize.fmin_cg
• TNC = Truncated Newton as implemented in scipy.optimize.fmin_tnc
• BFGS = Broyden–Fletcher–Goldfarb–Shanno method, as implemented in scipy.optimize.fmin_bfgs.
• L-BFGS = Limited-memory BFGS as implemented in scipy.optimize.fmin_l_bfgs_b. Contrary to the BFGS algorithm, which is written in Python, this one wraps a C implementation.
• Trust Region = Trust Region Newton method ¹. This is the solver used by LIBLINEAR that I've wrapped to accept any Python function in the package pytron

To assure the most accurate results across implementations, all timings were collected by callback functions that were called from the algorithm on each iteration. Finally, I plot the maximum absolute value of the gradient (=the infinity norm of the gradient) with respect to time.

The synthetic data used in the benchmarks was generated as described in ² and consists primarily of the design matrix $X$ being Gaussian noise, the vector of coefficients is drawn also from a Gaussian distribution and the explained variable $y$ is generated as $y = \text{sign}(X w)$. We then perturb matrix $X$ by adding Gaussian noise with covariance 0.8. The number of samples and features was fixed to $10^4$ and $10^3$ respectively. The penalization parameter $\lambda$ was fixed to 1.

In this setting variables are typically uncorrelated and most solvers perform decently:

Here, the Trust Region and L-BFGS solver perform almost equally good, with Conjugate Gradient and Truncated Newton falling shortly behind. I was surprised by the difference between BFGS and L-BFGS, I would have thought that when memory was not an issue both algorithms should perform similarly.

To make things more interesting, we now make the design to be slightly more correlated. We do so by adding a constant term of 1 to the matrix $X$ and adding also a column vector of ones this matrix to account for the intercept. These are the results:

Here, we already see that second-order methods dominate over first-order methods (well, except for BFGS), with Trust Region clearly dominating the picture but with TNC not far behind.

Finally, if we force the matrix to be even more correlated (we add 10. to the design matrix $X$), then we have:

Here, the Trust-Region method has the same timing as before, but all other methods have got substantially worse.The Trust Region method, unlike the other methods is surprisingly robust to correlated designs.

To sum up, the Trust Region method performs extremely well for optimizing the Logistic Regression model under different conditionings of the design matrix. The LIBLINEAR software uses this solver and thus has similar performance, with the sole exception that the evaluation of the logistic function and its derivatives is done in C++ instead of Python. In practice, however, due to the small number of iterations of this solver I haven't seen any significant difference.

This is a pretty good book - though it's more of a set of short little blog posts combined into a book. It's nice to read the very opinionated thoughts of someone who's been dealing with bicycles for a very long time.

The whole time I was reading this, I felt like this is the type of stuff that my pal Jon G would have said, and in a similar style to the way he would have said it. I look forward to reading Jon's treatises on cycling in the future.

Reading this book made me proud to and convinced me to continue to never wear clipless pedals and cycling shoes (I have toe clips), and having a kickstand on my roadie. I'm an unracer - and to quote the title of one of the little chapters: "Racing ruins the breed".

I also learned how to corner and turn properly - which has been invaluable in my windy descents lately. Finally, I was reminded to not just count miles, because every ride counts ("No ride too short"). Also it makes as much sense to count time spent on the bike, and the amount of elevation gain, and the number of days biked, period. Thanks to my commute, this means my days biked per week stays above 5 during the most of the year in California.

I also learned a new word: Beausage (byoo-sidj). My current bike doesn't have any, but my last laptop certainly does.

It also felt kind of cool that the author is local. Grant Petersen founded and still runs Rivendell Bicycle Works in Walnut Creek, CA.

Accrual Ratios are an important index for assessing the performance and continued viability of every publicly traded company. Using Wakari, TR CONNECT, and Thomson Reuters Financial Data I calculated accrual ratios for a variety of companies and explored news events, which may correlate with strong signals in the data. You can easily download and edit the notebook used in this post into your Wakari account.

The Nipy World blog used to be on Blogspot at http://nipyworld.blogspot.com. I'm restarting the blog here because it's easier and more fun to write the blog using Pelican. I hope that it will be easier for my fellow NIPistas to add their own posts using Github.

The Nipy World blog used to be on Blogspot at http://nipyworld.blogspot.com. I'm restarting the blog here because it's easier and more fun to write the blog using Pelican. I hope that it will be easier for my fellow NIPistas to add their own posts using Github.

I’m thrilled to announce that my paper “Block Coordinate Descent Algorithms for Large-scale Sparse Multiclass Classification” (published in the Machine Learning journal) is now online: PDF, BibTeX [*].

Abstract

Over the past decade, l1 regularization has emerged as a powerful way to learn classifiers with implicit feature selection. More recently, mixed-norm (e.g., l1/l2) regularization has been utilized as a way to select entire groups of features. In this paper, we propose a novel direct multiclass formulation specifically designed for large-scale and high-dimensional problems such as document classification. Based on a multiclass extension of the squared hinge loss, our formulation employs l1/l2 regularization so as to force weights corresponding to the same features to be zero across all classes, resulting in compact and fast-to-evaluate multiclass models. For optimization, we employ two globally-convergent variants of block coordinate descent, one with line search (Tseng and Yun, 2009) and the other without (Richtárik and Takáč, 2012). We present the two variants in a unified manner and develop the core components needed to efficiently solve our formulation. The end result is a couple of block coordinate descent algorithms specifically tailored to our multiclass formulation. Experimentally, we show that block coordinate descent performs favorably to other solvers such as FOBOS, FISTA and SpaRSA. Furthermore, we show that our formulation obtains very compact multiclass models and outperforms l1/l2- regularized multiclass logistic regression in terms of training speed, while achieving comparable test accuracy.

Code

The code of the proposed multiclass method is available in my Python/Cython machine learning library, lightning. Below is an example of how to use it on the News20 dataset.

from sklearn.datasets import fetch_20newsgroups_vectorized
from lightning.primal_cd import CDClassifier
 
bunch = fetch_20newsgroups_vectorized(subset="all")
X = bunch.data
y = bunch.target
 
clf = CDClassifier(penalty="l1/l2",
                   loss="squared_hinge",
                   multiclass=True,
                   max_iter=20,
                   alpha=1e-4,
                   C=1.0 / X.shape[0],
                   tol=1e-3)
clf.fit(X, y)
# accuracy
print clf.score(X, y) 
# percentage of selected features
print clf.n_nonzero(percentage=True)

To use the variant without line search (as presented in the paper), add the max_steps=0 option to CDClassifier.

Data

I also released the Amazon7 dataset used in the paper. It contains 1,362,109 reviews of Amazon products. Each review may belong to one of 7 categories (apparel, book, dvd, electronics, kitchen & housewares, music, video) and is represented as a 262,144-dimensional vector. It is, to my knowledge, one of the largest publically available multiclass classification dataset.

[*] The final publication is available here.

I’m please to announce a new version for scikits.optimization. The main focus of this iteration was to finish usual unconstrained optimization algorithms.

Changelog

• Fixes on the Simplex state implementation
• Added several Quasi-Newton steps (BFGS, rank 1 update…)

The scikit can be installed with pip/easy_install or downloaded from PyPI

Old announces:

• 0.2
• 0.1

Continuum Analytics, the premier provider of Python-based data analytics solutions and services, announced today the release of the latest version of Anaconda, its premium collection of libraries for Python. Anaconda enables big data management, analysis, and cross-platform visualization for business intelligence, scientific analysis, engineering, machine learning, and more. The latest release, version 1.5, includes the addition of netCDF4, as well as updates to Python, NumPy, SciPy, IPython, MatPlotLib, Pandas, and Cython.

After nine months on Octopress, I've decided to move on.

I should start by saying that Octopress is a great platform for static blogging: it's powerful, flexible, well-supported, well-integrated with GitHub pages, and has tools and plugins to do just about anything you might imagine. There's only one problem:

It's written in Ruby.

Now I don't have anything against Ruby per se. However, it was starting to seem a bit awkward that a blog called Pythonic Perambulations was built with Ruby, especially given the availability of so many excellent Python-based static site generators (Hyde, Nikola, and Pelican in particular).

Additionally, a few things with Octopress were starting to become difficult: first, I wanted a way to easily insert IPython notebooks into posts. Sure, I developed a hackish solution to notebooks in Octopress which had worked well enough for a while, but a cleaner method would have involved digging into the Ruby source code and writing a full-fledged Octopress extension, based on nbconvert. This would have involved a fair bit of effort to learn Ruby and figure out how to best interface it with the Python nbconvert code. Second, Ruby has so many strange and difficult pieces: GemFiles, RVM, rake... and I never took the time to really understand the real purpose of all of them (self-reflection: what parts of Python would seem strange and difficult if I hadn't been using them for so many years?). The black-box nature of the process, at least in my own case, was starting to bother me.

But the kicker was this: In January I got a new computer, and after a reasonable amount of effort was unable to successfully build my blog. I've been writing my posts exclusively on my old laptop which I somehow managed to successfully set up last August. But that laptop now has a sorely outdated Ubuntu distro that I couldn't upgrade for fear of losing the ability to update my blog. Needless to say, this was not the most effective setup.

It was time to switch my blog engine to Python.

Choosing a Python Static Generator

I started asking around, and found that there were three solid contenders for a Python-based platform to replace Octopress: Hyde, Nikola, and Pelican. I gave Hyde a test-run a few weeks ago in re-making my website, and I really like it: it's clean, straightforward, powerful, and easy to use. The documentation is a bit lacking, though, and I think it would take a fair bit more effort at this point to build a more complicated site with Hyde.

Nikola and Pelican both seem to be well-loved by their users, but I had to choose one. I went with Pelican for one simple reason: it has more GitHub forks. I'm sure this is entirely unfair to Nikola and all the contributors who have poured their energy into the project, but I had to choose one way or another. I'm pleased to say that Pelican has not been a disappointment: I've found it to be flexible and powerful. It has an active developer-base, and makes available a wide array of themes and plugins. For the few extra pieces I needed, I found the plugin and theming API to be well-documented and straightforward to use.

Migrating to Pelican from Octopress

I won't attempt to write a one-size-fits-all guide to migrating to Pelican from Octopress: there are too many possibile combinations of formats, plugins, themes, etc. But I will walk through my own process in some detail, in hopes that it might help others who find themselves in a similar predicament.

I had several goals when doing this migration:

• I wanted, as much as possible, to maintain the look and feel of the blog. I like the default Octopress theme: it's simple, clean, compact, and includes all the aspects I need for a good blog.
• I wanted, as much as possible, to leave the source of my posts unmodified: luckily Pelican supports writing posts in markdown and allows easy insertion of custom plugins, so this was relatively easy to accommodate.
• I wanted to maintain the history of Disqus comments for each page, as well as the Twitter and Google Pages tools.
• I wanted, as much as possible, to maintain the same URLs for all content, including posts, notebooks, images, and videos.
• I wanted a clean way to insert html-formatted IPython notebooks into blog posts. Nearly half my posts are written as notebooks, and the old way of including them was becoming much too cumbersome.

I was able to suitably address all these goals with Pelican in a few evenings' effort. Some of it was already built-in to the Pelican settings architecture, some required customization of themes and extensions, and some required writing some brand new plugins. I'll summarize these aspects below:

Blog theme

As I mentioned, I wanted to keep the look and feel of the blog consistent. Luckily, someone had gone before me and created an octopress Pelican theme which did most of the heavy lifting. I contributed a few additional features, including the ability to specify Disqus tags and maintain comment history, to add Twitter, Google Plus, and Facebook links in the sidebar and footer, to add a custom site search box which appears in the upper right of the navigation panel, as well as a few other tweaks. The result is what you see here: nearly identical to the old layout, with all the bells and whistles included.

Octopress Markdown to Pelican Markdown

Octopress has a few plugins which add some syntactic sugar to the markdown language. These are tags specified in Liquid-style syntax:

{% tag arg1 arg2 ... %}

I have made extensive use of these in my octopress posts, primarily to insert videos, images, and code blocks from file. In order to accommodate this in Pelican, I wrote a Pelican plugin which wraps a custom Markdown preprocessor written via the Markdown extension API which can correctly interpret these types of tags. The tags ported from octopress thus far are:

The Image Tag

The image tag allows insertion of an image into the post with a specified size and position:

{% img [position] /url/to/img.png [width] [height] [title] [alt] %}

Here is an example of the result of the image tag:

The Video Tag

The video tag allows embedding of an HTML5/Flash-compatible video into the post:

{% video /url/to/video.mp4 [width] [height] [/url/to/poster.png] %}

Here is an example of the output of the video tag:

(see this post for a description of this video).

The Code Include Tag

The include_code tag allows the insertion of formatted lines from a file into the post, with a title and a link to the source file:

{% include_code filename.py [lang:python] [title] %}

Here is an example of the output of the code include tag:

import sys
import os
print("hello_world")

For more information on using these tags, refer to the module doc-strings.

Maintaining Disqus Comment Threads

Static blogs are fast, lightweight, and easy to deploy. A disadvantage, though, is the inability to natively include dynamic elements such as comment threads. Disqus is a third-party service that skirts this disadvantage very seamlessly. All it takes is to add a small javascript snippet with some identifiers in the appropriate place on your page, and Disqus takes care of the rest. To keep the comment history on each page required assuring that the site identifier and page identifiers remained the same between blog versions. This part happens within the theme, and my Disqus PR to the Pelican Octopress theme made this work correctly.

Maintaining the URL structure

By default, Octopress stores posts with a structure looking like blog/YYYY/MM/DD/article-slug/. The Pelican default is different, but easy enough to change. In the pelicanconf.py settings file, this corresponds to the following:

ARTICLE_URL = 'blog/{date:%Y}/{date:%m}/{date:%d}/{slug}/'
ARTICLE_SAVE_AS = 'blog/{date:%Y}/{date:%m}/{date:%d}/{slug}/index.html'

Next, at the top of the markdown file for each article, the metadata needs to be slightly modified from the form used by Octopress -- here is the actual metadata used in the document that generates this page:

Title: Migrating from Octopress to Pelican
date: 2013-05-07 17:00
slug: migrating-from-octopress-to-pelican

Additionally, the static elements of the blog (images, videos, IPython notebooks, code snippets, etc.) must be put within the correct directory structure. These static files should be put in paths which are specified via the STATIC_PATHS setting:

STATIC_PATHS = ['images', 'figures', 'downloads']

Pelican presented a challenge here: as of the time of this writing, Pelican has a hard-coded 'static' subdirectory where these static paths are saved. I submitted a pull request to Pelican that replaces this hard-coded setting with a configurable path: because the change conflicts with a bigger refactoring of the code which is ongoing, the PR will not be merged. But until this new refactoring is finished, I'll be using my own branch of Pelican to make this blog, and specify the correct static paths.

Inserting Notebooks

The ability to seamlessly insert IPython notebooks into posts was one of the biggest drivers of my switch to Pelican. Pelican has an ipython notebook plugin available, but I wasn't completely happy with it. The plugin implements a reader which treats the notebooks themselves as the source of a post, leading to the requirement to insert blog metadata into the notebook itself. This is a suitable solution, but for my own purposes I much prefer a solution in which the content of a notebook could be inserted into a stand-alone post, such that the notebook and the blog metadata are completely separate.

To accomplish this, I added a submodule to my liquid_tags Pelican plugin which allows the insertion of notebooks using the following syntax:

{% notebook path/to/notebook.ipynb [cells[i:j]] %}

This inserts the notebook at the given location in the post, optionally selecting a given range of notebook cells with standard Python list slicing syntax.

The formatting of notebooks requires some special CSS styles which must be inserted into the header of each page where notebooks are shown. Rather than requiring the theme to be customized to support notebooks, I decided on a solution where an EXTRA_HEADER setting is used to specify html and CSS which should be added to the header of the main page. The notebook plugin saves the required header to a file called _nb_header.html within the main source directory. To insert the appropriate formatting, we add the following lines to the settings file, pelicanconf.py:

EXTRA_HEADER = open('_nb_header.html').read().decode('utf-8')

In the theme file, within the header tag, we add the following:

 {% if EXTRA_HEADER %}
 {{ EXTRA_HEADER }}
 {% endif %}

Here is the result: a short notebook inserted directly into the post:

import numpy
print numpy.random.random(10)

[ 0.25463203  0.55637185  0.02743774  0.57380221  0.52378531  0.95099357
  0.70975568  0.19575853  0.589227    0.06959599]

With all those things in place, the blog was ready to be built. The result is right in front of you: you're reading it. If you'd like to see the source from which this blog built, it's available at http://www.github.com/jakevdp/PythonicPerambulations. Feel free to adapt the configurations and theme to suit your own needs.

I'm glad to be working purely in Python from now on!

Welcome. This post is part of a series of Continuum Analytics Open Notebooks showcasing our projects, products, and services.

In this Continuum Open Notebook, you’ll see how Numba accelerates the performance of the GrowCut image segmentation window function by three orders of magnitude in two lines of Python.

In the next version of scikits.optimization, I’ve added some Quasi-Newton steps. Before this version is released, I thought I would compare several methods of optimizing the Rosenbrock function.

Optimizers

What is great with the Rosenbrock cost function can be summed up in a few points:

1. It is hard to optimize
2. Gradient can be easily computed

I’ve decided to compare the number of function and gradient calls as well as the cost behavior for several usual optimization algorithms. So the contestants will be:

• A Nelder-Mead/Polytope/simplex optimizer
• SSA, a simplex with simulated annealing (think of amotsa from Numerical Recipes)
• GA, a genetic algorithm
• a gradient-based optimizer
• CG, a conjugate-gradient optimizer
• BFGS, a quasi-Newton optimizer

The first 4 are from scikits.optimization, the last 2 are based on proprietary code that cannot be published, but it’s interesting to compare with other tools that are used to compare gradient-free complex cost functions.

Results

I’ve made a small slideshow with the derivative-free algorithms. First you have for each of the three algorithms the number of function calls versus iteration, the cost versus iteration and finally the location of testing parameters.
Click to view slideshow.

This slideshow is for the derivative-based algorithms.
Click to view slideshow.

Conclusion

I was quite surprised by some algorithm behaviors. Clearly, the conjugate gradient algorithm behaves far better than the simple gradient, but the BFGS followed the Rosenbrock valley far better. A good quasi-Newton can be really efficient (not a brent because it needs to solve a linear equation), although a conjugate gradient can be enough in some cases.

For the gradient-free algorithms, SSA really behaved badly. This is mainly due because the hyperparameters that must be adequately tuned. This function is quite simple, but my first trial at setting these parameters was far more efficient for GA or the simplex than for SSA. So I would go for GA for gradient-free optimization: few and easy hyper parameters and a good browse of the search space.

The code for the 4 first tests and display plots

TL;DR: I've implemented a logistic ordinal regression or proportional odds model. Here is the Python code

The logistic ordinal regression model, also known as the proportional odds was introduced in the early 80s by McCullagh [¹, ²] and is a generalized linear model specially tailored for the case of predicting ordinal variables, that is, variables that are discrete (as in classification) but which can be ordered (as in regression). It can be seen as an extension of the logistic regression model to the ordinal setting.

Given $X \in \mathbb{R}^{n \times p}$ input data and $y \in \mathbb{N}^n$ target values. For simplicity we assume $y$ is a non-decreasing vector, that is, $y_1 \leq y_2 \leq ...$. Just as the logistic regression models posterior probability $P(y=j|X_i)$ as the logistic function, in the logistic ordinal regression we model the cummulative probability as the logistic function. That is,

$$ P(y \leq j|X_i) = \phi(\theta_j - w^T X_i) = \frac{1}{1 + \exp(w^T X_i - \theta_j)} $$

where $w, \theta$ are vectors to be estimated from the data and $\phi$ is the logistic function defined as $\phi(t) = 1 / (1 + \exp(-t))$.

Compared to multiclass logistic regression, we have added the constrain that the hyperplanes that separate the different classes are parallel for all classes, that is, the vector $w$ is common across classes. To decide to which class will $X_i$ be predicted we make use of the vector of thresholds $\theta$. If there are $K$ different classes, $\theta$ is a non-decreasing vector (that is, $\theta_1 \leq \theta_2 \leq ... \leq \theta_{K-1}$) of size $K-1$. We will then assign the class $j$ if the prediction $w^T X$ (recall that it's a linear model) lies in the interval $[\theta_{j-1}, \theta_{j}[$. In order to keep the same definition for extremal classes, we define $\theta_{0} = - \infty$ and $\theta_K = + \infty$.

The intuition is that we are seeking a vector $w$ such that $X w$ produces a set of values that are well separated into the different classes by the different thresholds $\theta$. We choose a logistic function to model the probability $P(y \leq j|X_i)$ but other choices are possible. In the proportional hazards model ¹ the probability is modeled as $-\log(1 - P(y \leq j | X_i)) = \exp(\theta_j - w^T X_i)$. Other link functions are possible, where the link function satisfies $\text{link}(P(y \leq j | X_i)) = \theta_j - w^T X_i$. Under this framework, the logistic ordinal regression model has a logistic link function and the proportional hazards model has a log-log link function.

The logistic ordinal regression model is also known as the proportional odds model, because the ratio of corresponding odds for two different samples $X_1$ and $X_2$ is $\exp(w^T(X_1 - X_2))$ and so does not depend on the class $j$ but only on the difference between the samples $X_1$ and $X_2$.

Optimization

Model estimation can be posed as an optimization problem. Here, we minimize the loss function for the model, defined as minus the log-likelihood:

$\begin{align} \mathcal{L}(w, \theta) &= - \sum_{i=1}^n \log(\phi(\theta_{y_i} - w^T X_i) - \phi(\theta_{y_i -1} - w^T X_i)) \\ &= \sum_{i=1}^n w^T X_i - \log(\exp(\theta_{y_i}) - \exp(\theta_{y_i-1})) + \log(\phi(\theta_{y_i -1} - w^T X_i)) + \log(\phi(\theta_{y_i} - w^T X_i)) \\ \end{align}$

In this sum all terms are convex on $w$, thus the loss function is convex over $w$. It might be also jointly convex over $w$ and $\theta$, although I haven't checked. I use the function fmin_slsqp in scipy.optimize to optimize $\mathcal{L}$ under the constraint that $\theta$ is a non-decreasing vector. There might be better options, I don't know. If you do know, please leave a comment!.

Using the formula $\log(\phi(t))^\prime = (1 - \phi(t))$, we can compute the gradient of the loss function as

$\begin{align} \nabla_w \mathcal{L}(w, \theta) &= \sum_{i=1}^n X_i (1 - \phi(\theta_{y_i} - w^T X_i)) - \phi(\theta_{y_i-1} - w^T X_i)) \\ % \nabla_\theta \mathcal{L}(w, \theta) &= \sum_{i=1}^n - \frac{e_{y_i} \exp(\theta_{y_i}) - e_{y_i -1} \exp(\theta_{y_i -1})}{\exp(\theta_{y_i}) - \exp(\theta_{y_i-1})} \\ \nabla_\theta \mathcal{L}(w, \theta) &= \sum_{i=1}^n e_{y_i} \left(1 - \phi(\theta_{y_i} - w^T X_i) - \frac{1}{1 - \exp(\theta_{y_i -1} - \theta_{y_i})}\right) \\ & \qquad + e_{y_i -1}\left(1 - \phi(\theta_{y_i -1} - w^T X_i) - \frac{1}{1 - \exp(- (\theta_{y_i-1} - \theta_{y_i}))}\right) \end{align}$

where $e_i$ is the $i$th canonical vector.

Code

I've implemented a Python version of this algorithm using Scipy's optimize.fmin_slsqp function. This takes as arguments the loss function, the gradient denoted before and a function that is > 0 when the inequalities on $\theta$ are satisfied.

Code can be found here as part of the minirank package, which is my sandbox for code related to ranking and ordinal regression. At some point I would like to submit it to scikit-learn but right now the I don't know how the code will scale to medium-scale problems, but I suspect not great. On top of that I'm not sure if there is a real demand of these models for scikit-learn and I don't want to bloat the package with unused features.

Performance

I compared the prediction accuracy of this model in the sense of mean absolute error (IPython notebook) on the boston house-prices dataset. To have an ordinal variable, I rounded the values to the closest integer, which gave me a problem of size 506 $\times$ 13 with 46 different target values. Although not a huge increase in accuracy, this model did give me better results on this particular dataset:

Here, ordinal logistic regression is the best-performing model, followed by a Linear Regression model and a One-versus-All Logistic regression model as implemented in scikit-learn.

We are happy to announce a new Anaconda add-on product called “MKL Optimizations”, which allows packages in Anaconda to take advantage of the Math Kernel Library (MKL) by Intel.