Quadrature routines

[sglyon]

I think the quadrature routines from Miranda and Fackler are ready to be incorporated into the master branch.

I'm opening this PR to start the review process.

Eventually, we should probably update the examples notebook to provide a bit more exposition explaining how to use the functions as well as adding the new demqua##.m routines in Miranda's 2014 version of CompEcon.

[sglyon]

Hey John, I think I have gone through the points in your email from 07-15 and this is ready to merge.

I added some of Mario's new examples (in the CE 2014 source) to the example notebook.

There are still a couple of minor things to clean up on this code (a couple of docstrings and a reminder that we could beef up the pedagogy in the example notebook) and I have documented it in issue #31.

If you are ok merging, I think it is ready.

[jstac]

Well done guys, great work finishing this up.

[albop]

I was writing a comment, but I was slightly too late and John closed the issue before.
I have read parts of the code quickly and I have a few concerns. They can still be dealt with later in the quad branch and merged later. The two first should be fixed faster, since they can have side-effects on future code.

• the gridmake function does not enumerate points in a way that is consistent with default Python conventions. If you compute gridmake(array([0,1]), ([-1,2])) it produces:

array([[ 0, -1],
       [ 1, -1],
       [ 0, -2],
       [ 1, -2],

meaning that the first order varies faster. Now, if you want to represent values on this grid by a 2d array vals such that vals[i,j] contains values, then when you do vals.ravel(), you don't enumerate points in the same order because last index is supposed to vary faster. This is quite annoying.

• Multidimensional functions don't always return an array of the same dimensionality. That is not a problem in Matlab, but in Python, if I do quad. qnwnorm([2,2]) I get 2-dimensional arrays and with qnwnorm([2]) I get 1-dimensional vector. This is problematic, since in generic code, one will always have to distinguish dimension 1, from higher dimensions. My opinion here is that multidimensional routines, should always return multidimensional objects in a predictible fashion. (maybe the 1-d routines could be exported too)
• (minor) It looks like qnwnorm will fail if one column of the covariance matrix (and the corresponding line) is full of zeros. That was already the case in the compecon toolbox, but it is still a common usecase.
• (minor) in the future, we may want to have more engaging names than qnwnorm and co, don't we ?

[albop]

A possible idea would be to rename quad.py into ce_quad.py so that the latter can be left untouched and remain as close as possible from the original version. A quadrature.py could then contain Python compliant versions. What do you think ?

[albop]

As for a gridmake replacement, the function cartesian does the same thing. It is also faster. (cf http://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-all-combinations-of-two-arrays)