Mailman 3 October 2010 - NumPy-Discussion

Hi Travis
by Charles R Harris 30 Oct '10

30 Oct '10

Hi Travis, Could you look over the thread about dtype comparisons<http://thread.gmane.org/gmane.comp.python.numeric.general/41074/focus=41082>and offer an opinion? TIA, Chuck

3 3

ValueError with np.polynomial.Polynomial([1]).roots()
by josef.pktd＠gmail.com 29 Oct '10

29 Oct '10

For the lag polynomials I have cases where the order is zero, but roots raises a ValueError. Is this intended? Then, I need to catch it in my code. I haven't checked yet if my other parts will go through with empty roots. >>> np.polynomial.Polynomial([1]) Polynomial([ 1.], [-1., 1.]) >>> np.polynomial.Polynomial([1]).roots() Traceback (most recent call last): File "<pyshell#140>", line 1, in <module> np.polynomial.Polynomial([1]).roots() File "<string>", line 485, in roots File "C:\Programs\Python25\lib\site-packages\numpy\polynomial\polyutils.py", line 280, in mapdomain [x] = as_series([x], trim=False) File "C:\Programs\Python25\lib\site-packages\numpy\polynomial\polyutils.py", line 137, in as_series raise ValueError("Coefficient array is empty") ValueError: Coefficient array is empty Josef

2 2

ndarray __getattr__ to perform __getitem__
by Ian Stokes-Rees 29 Oct '10

29 Oct '10

I have an ndarray with named dimensions. I find myself writing some fairly laborious code with lots of square brackets and quotes. It seems like it wouldn't be such a big deal to overload __getattribute__ so instead of doing: r = genfromtxt('results.dat',dtype=[('a','int'), ('b', 'f8'), ('c','int'), ('d', 'a20')]) scatter(r[r['d'] == 'OK']['a'], r[r['d'] == 'OK']['b']) I could do: scatter(r[r.d == 'OK'].a, r[r.d == 'OK'].b) which is really a lot clearer. Is something like this already possible somehow? Is there some reason not to map __getattr__ to __getitem__? class ndarray(): ... __getattr__ = __getitem__ Or something somewhat more sophisticated like: def __getattr__(self, attr): if self.hasitem(attr): return self[attr] else: raise AttributeError(attr) TIA for any comments or thoughts on this. Ian

5 9

C extension compiling question
by Henry Gomersall 29 Oct '10

29 Oct '10

I'm trying to get a really simple toy example for a numpy extension working (you may notice its based on the example in the numpy docs and the python extension docs). The code is given below. The problem I am having is running the module segfaults at any attempt to access &PyArray_Type (so, as presented, the segfault occurs at the call to PyArg_ParseTuple). I've been assuming that importing numpy into the calling python code implicitly provides the numpy libraries that are necessary for linking. I can only assume I'm missing something. Do i need some additional includes etc in my compile flags? I'm currently using the distutils.core module and not the distutils module that comes with numpy if that makes a difference. I'm actually pretty keen to document what I've learnt as a 5 minute how-to for writing basic numpy extensions, if there is demand for it. It doesn't seem to exist as a coherent whole at the moment. Thanks, Henry #include <Python.h> #include <numpy/arrayobject.h> PyMODINIT_FUNC initspam(void); static PyObject * spam_system(PyObject *self, PyObject *args); void worker_function(double *input, double *output, npy_intp size); static PyMethodDef SpamMethods[] = { {"system", spam_system, METH_VARARGS, "Execute a shell command."}, {NULL, NULL, 0, NULL} /* Sentinel */ }; PyMODINIT_FUNC initspam(void) { (void) Py_InitModule("spam", SpamMethods); } static PyObject * spam_system(PyObject *self, PyObject *args) { int n; int ndims; PyObject *inarg=NULL, *outarg=NULL; PyObject *inarr=NULL, *outarr=NULL; /* Both input and output need to be PyArray_Type type*/ if (!PyArg_ParseTuple(args, "O!O!", &PyArray_Type, &inarg, &PyArray_Type, &outarg)); return NULL; inarr = PyArray_FROM_OTF(inarg, NPY_DOUBLE, NPY_IN_ARRAY); if (inarr == NULL) goto fail; outarr = PyArray_FROM_OTF(outarg, NPY_DOUBLE, NPY_OUT_ARRAY); if (outarr == NULL) goto fail; /* Check the shape of the two arrays is the same */ ndims = ((PyArrayObject *)outarr)->nd; if (((PyArrayObject *)inarr)->nd != ndims) goto fail; /* ...still checking...*/ for (n=1; n<ndims; n++) { if (PyArray_DIM(inarr,n) != PyArray_DIM(outarr,n)) goto fail; } /* Now actually do something with the arrays */ worker_function((double *)PyArray_DATA(inarr), (double *)PyArray_DATA(outarr), PyArray_Size(inarr)); Py_DECREF(inarr); Py_DECREF(outarr); return Py_None; fail: Py_XDECREF(inarr); Py_XDECREF(outarr); return NULL; } void worker_function(double *input, double *output, npy_intp size) { npy_intp n; for (n=0;n<size;n++) { output[n] = input[n] + 10; } return; }

4 4

genfromtxt behaviour
by Matt Studley 29 Oct '10

29 Oct '10

<snip> >> How can I do my nice 2d slicing on the latter? >> >> array([('a', 2, 3), ('b', 5, 6), ('c', 8, 9)], >> dtype=[('f0', '|S1'), ('f1', '<i4'), ('f2', '<i4')]) >Select a column by its name: >yourarray['f0'] Super! So I would need to get the dtype object... myData[ myData.dtype.names[0] ] in order to index by column. Matt

2 1

genfromtxt behaviour
by Matt Studley 29 Oct '10

29 Oct '10

Hi all first, please forgive me for my ignorance - I am taking my first stumbling steps with numpy and scipy. I am having some difficulty with the behaviour of genfromtxt. s = SIO.StringIO("""1, 2, 3 4, 5, 6 7, 8, 9""") g= genfromtxt(s, delimiter=', ', dtype=None) print g[:,0] This produces the output I expected, a slice for column 0... array([1, 4, 7]) BUT.... s = SIO.StringIO("""a, 2, 3 b, 5, 6 c, 8, 9""") g= genfromtxt(s, delimiter=', ', dtype=None) g[:,0] Produces the following error: IndexError: invalid index .... In the first case, genfromtxt returns me a 2d array (list of lists), in the second it returns a list of tuples with an associated dtype list of tuples. How can I do my nice 2d slicing on the latter? array([('a', 2, 3), ('b', 5, 6), ('c', 8, 9)], dtype=[('f0', '|S1'), ('f1', '<i4'), ('f2', '<i4')]) thanks in advance for any assistance Matt

2 1

Inconsistency with __index__() for rank-1 arrays?
by Francesc Alted 29 Oct '10

29 Oct '10

Hi, I find this a bit misleading: >>> a = np.arange(10) >>> a[np.array(0)] 0 >>> a[np.array([0])] array([0]) >>> a[[0]] array([0]) But, for regular python lists we have: >>> l = a.tolist() >>> l[np.array(0)] 0 >>> l[np.array([0])] 0 i.e. indexing with a rank-0 array and a rank-1 array with one single element return the same result, which I find inconsistent with the expected behaviour for this case, i.e.: >>> l[[0]] --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /tmp/tables-2.2/<ipython console> in <module>() TypeError: list indices must be integers, not list The ultimate reason for this behaviour is this: >>> np.array(0).__index__() 0 >>> np.array([0]).__index__() 0 But I wonder why NumPy needs the latter behaviour, instead of the more logical: >>> np.array([0]).__index__() --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /tmp/tables-2.2/<ipython console> in <module>() TypeError: only rank-0 integer arrays can be converted to an index This inconsistency has indeed introduced a bug in my application and for solving this I'd need something like: """ def is_idx(index): """Check if an object can work as an index or not.""" if hasattr(index, "__index__"): # Only works on Python 2.5 on if (hasattr(index, "shape") and index.shape == (1,)): return False try: # (as per PEP 357) idx = index.__index__() return True except TypeError: return False return False """ i.e. for determining if an object can be an index or not, I need to explicitly check for a shape different from (1,), which is unnecessarily complicated. So I find the current behaviour prone to introduce errors in apps and I'm wondering why exactly np.array([1]) should work as an index at all. It would not be better if that would raise a ``TypeError``? Thanks, -- Francesc Alted

3 6

Re: [Numpy-discussion] problems with numdifftools
by Nicolai Heitz 28 Oct '10

28 Oct '10

m 27.10.2010 02:02, schrieb Sebastian Walter: > On Wed, Oct 27, 2010 at 12:59 AM, Pauli Virtanen<pav(a)iki.fi> wrote: >> Tue, 26 Oct 2010 14:24:39 -0700, Nicolai Heitz wrote: >>>> http://mail.scipy.org/mailman/listinfo/scipy-user >>> I contacted them already but they didn't responded so far and I was >>> forwarded to that list which was supposed to be more appropriated. >> I think you are thinking here about some other list -- scipy-user >> is the correct place for this discussion (and I don't remember seeing >> your mail there). I was pretty sure that I put it there. Unfortunately it had a different name there: Errors in calculation of the Hessian using numdifftools. But I can't find the post by myself at the moment so maybe something went wrong. >> [clip] >>> 1) Can I make it run/fix it, so that it is also going to work for the SI >>> scaling? >> Based on a brief look, it seems that uniform scaling will not help you, >> as you have two very different length scales in the problem, >> >> 1/sqrt(m w^2)>> C >> >> If you go to CM+relative coordinates you might be able to scale them >> separately, but that's fiddly and might not work for larger N. >> >> In your problem, things go wrong when the ratio between the >> length scales approaches 1e-15 which happens to be the machine epsilon. >> This implies that the algorithm runs into some problems caused by the >> finite precision of floating-point numbers. >> >> What exactly goes wrong and how to fix it, no idea --- I didn't look into >> how Numdifftools is implemented. Probably you are right. I converted it to natural units and it worked out in the 2 ion case. Increasing the number of ions leads to problems again and I have no idea where those problems come from. >>> 2) How can I be sure that increasing the number of ions or adding a >>> somehow more complicated term to the potential energy is not causing the >>> same problems even in natural units? >>> >>> 3) In which range is numdifftools working properly. >> That depends on the algorithm and the problem. Personally, I wouldn't >> trust numerical differentiation if the problem has significantly >> different length scales, it is important to capture all of them >> accurately, and it is not clear how to scale them to the same size. >> Writing ND software that works as expected all the time is probably >> not possible even in theory. >> >> Numerical differentiation is not the only game in the town. I'd look >> into automatic differentiation (AD) -- there are libraries available >> for Python also for that, and it is numerically stable. >> >> E.g. >> >> http://en.wikipedia.org/wiki/Automatic_differentiation#Software >> >> has a list of Python libraries. I don't know which of them would be >> the best ones, though. >> > they all have their pro's and con's. > Being (co-)author of some of these tools, my personal and very biased advice is: > if you are on Linux, I would go for PYADOLC. it provides bindings to a > feature-rich and well-tested C++ library. > However, the installation is a little tricker than a "setup.py build" > since you will need to compile ADOL-C and get Boost::Python to work. > PYADOLC can also differentiate much more complicated code than your > example in a relatively efficient manner. Is there by chance any possibility to make PYADOLC run on a (lame) windows xp engine. If not what else would u recommend (besides switching to Linux, what I am going to do soon). > For your example the code looks like: > > --------------------- code --------------------------- > > .... > > c=classicalHamiltonian() > xopt = optimize.fmin(c.potential, c.initialposition(), xtol = 1e-10) > > import adolc; import numpy > > # trace the computation > adolc.trace_on(0) > x = adolc.adouble(c.initialposition()) > adolc.independent(x) > y = c.potential(x) > adolc.dependent(y) > adolc.trace_off() > > hessian = adolc.hessian(0, xopt) > eigenvalues = numpy.linalg.eigh(hessian)[0] > normal_modes = c.normal_modes(eigenvalues) > print 'hessian=\n',hessian > print 'eigenvalues=\n',eigenvalues > print 'normal_modes=\n',normal_modes > --------------------- code --------------------------- > > and you get as output > Optimization terminated successfully. > Current function value: 0.000000 > Iterations: 81 > Function evaluations: 153 > hessian= > [[ 5.23748399e-12 -2.61873843e-12] > [ -2.61873843e-12 5.23748399e-12]] > eigenvalues= > [ 2.61874556e-12 7.85622242e-12] > normal_modes= > [ 6283185.30717959 10882786.30440101] > > > Also, you should use an eigenvalue solver for symmetric matrices, e.g. > numpy.linalg.eigh. > Your code example looks awesome and leads to the correct results. Thank you very much. I try to make it work on my pc as well. > regards, > Sebastian > >> -- >> Pauli Virtanen >> >> _______________________________________________ >> NumPy-Discussion mailing list >> NumPy-Discussion(a)scipy.org >> http://mail.scipy.org/mailman/listinfo/numpy-discussion >> > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion(a)scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion

2 2

quadratic function
by Brennan Williams 28 Oct '10

28 Oct '10

I have used both linear least squares and radial basis functions as a proxy equation, calculated from the results of computer simulations which are calculating some objective function value based on a number of varied input parameters. As an alternative option I want to add a quadratic function so if there are parameters/variables x,y,z then rather than just having a linear function f=a+bx+cy+dz I'll have f=a+bx+cx**2 + dxy + .... I'd like to have the option not to include all the different second order terms. Where should I start looking? Thanks Brennan

2 5

Stacking a 2d array onto a 3d array
by Dewald Pieterse 28 Oct '10

28 Oct '10

Starting with: In [93]: test = > numpy.array([[[1,1,1],[1,1,1]],[[2,2,2],[2,2,2]],[[3,3,3],[3,3,3]]]) > > In [94]: test > Out[94]: > array([[[1, 1, 1], > [1, 1, 1]], > > [[2, 2, 2], > [2, 2, 2]], > > [[3, 3, 3], > [3, 3, 3]]]) > > Slicing the complete first row: > > In [95]: firstrow = test[0,:,:] > > In [96]: firstrow > Out[96]: > array([[1, 1, 1], > [1, 1, 1]]) > I want to stack firstrow onto test to end up with: ([[[1, 1, 1], > [1, 1, 1]], > > [[1, 1, 1], > [1, 1, 1]], > > [[2, 2, 2], > [2, 2, 2]], > > [[3, 3, 3], > [3, 3, 3]]] > vstack wants the array dimensions to be the same, is this possible without doing 1 dimensional reshape, the actual data I want to do this on is some what larger. numpy.vstack((firstrow,test)) > --------------------------------------------------------------------------- > ValueError Traceback (most recent call last) > > /mnt/home/home/bmeagle/M/programme/analiseerverwerkteprent.py in <module>() > ----> 1 > 2 > 3 > 4 > 5 > > /usr/lib64/python2.6/site-packages/numpy/core/shape_base.py in vstack(tup) > 212 > 213 """ > --> 214 return _nx.concatenate(map(atleast_2d,tup),0) > 215 > 216 def hstack(tup): > > ValueError: arrays must have same number of dimensions > What is the correct python way to do this? -- Dewald Pieterse

3 3