Improvements for quadrature routines

Since the quad branch is closed, I'm opening a new issue with the comments I made about it, so that we can start a fresh discussion.
There are currently a few issues with the quadrature routines. The two items can have side-effects on future code and should be fixed quickly.

• 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 one is quite annoying.

• Multidimensional functions don't always return an array with the same number of dimensions. 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 vectors. This is problematic, since in generic code, it will force one to always 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 ?

A possible way to deal with these issues 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 (including Fortran order). A quadrature.py could then contain Python compliant versions, possibly with more explicit names.

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

Actually all these issues concern only the multidimensional functions.

Start by reviewing quad.py and the existing gridmake and qnwnorm behavior described in the issue. Compare the proposed cartesian replacement and the possible ce_quad.py and quadrature.py split, then clarify the API for multidimensional outputs and zero covariance columns. Done means the project agrees on a Python-consistent design for all listed concerns.