### Re: [Numpy-discussion] caching large allocations on gnu/linux

On Mon, Mar 13, 2017 at 12:21 PM Julian Taylor < jtaylor.deb...@googlemail.com> wrote: Should it be agreed that caching is worthwhile I would propose a very > simple implementation. We only really need to cache a small handful of > array data pointers for the fast allocate deallocate cycle that appear > in common numpy usage. > For example a small list of maybe 4 pointers storing the 4 largest > recent deallocations. New allocations just pick the first memory block > of sufficient size. > The cache would only be active on systems that support MADV_FREE (which > is linux 4.5 and probably BSD too). > > So what do you think of this idea? > This is an interesting thought, and potentially a nontrivial speedup with zero user effort. But coming up with an appropriate caching policy is going to be tricky. The thing is, for each array, numpy grabs a block "the right size", and that size can easily vary by orders of magnitude, even within the temporaries of a single expression as a result of broadcasting. So simply giving each new array the smallest cached block that will fit could easily result in small arrays in giant allocated blocks, wasting non-reclaimable memory. So really you want to recycle blocks of the same size, or nearly, which argues for a fairly large cache, with smart indexing of some kind. How much difference is this likely to make? Note that numpy is now in some cases able to eliminate allocation of temporary arrays. I think the only way to answer these questions is to set up a trial implementation, with user-switchable behaviour (which should include the ability for users to switch it on even when MADV_FREE is not available) and sensible statistics reporting. Then volunteers can run various numpy workloads past it. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] float16/32: wrong number of digits?

On Mon, Mar 13, 2017 at 12:57 PM Eric Wieserwrote: > > `float(repr(a)) == a` is guaranteed for Python `float` > > And `np.float16(repr(a)) == a` is guaranteed for `np.float16`(and the same > is true up to `float128`, which can be platform-dependent). Your code > doesn't work because you're deserializing to a higher precision format than > you serialized to. > I would hesitate to make this guarantee - certainly for old versions of numpy, np.float128(repr(x))!=x in many cases. I submitted a patch, now accepted, that probably accomplishes this on most systems (in fact this is now in the test suite) but if you are using a version of numpy that is a couple of years old, there is no way to convert long doubles to human-readable or back that doesn't lose precision. To repeat: only in recent versions of numpy can long doubles be converted to human-readable and back without passing through doubles. It is still not possible to use % or format() on them without discarding all precision beyond doubles. If you actually need long doubles (and if you don't, why use them?) make sure your application includes a test for this ability. I recommend checking repr(1+np.finfo(np.longdouble).eps). Anne P.S. You can write (I have) a short piece of cython code that will reliably repr and back long doubles, but on old versions of numpy it's just not possible from within python. -A ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] float16/32: wrong number of digits?

On Thu, Mar 9, 2017, 11:27 Nico Schlömerwrote: > Hi everyone, > > I wondered how to express a numpy float exactly in terms of format, and > found `%r` quite useful: `float(repr(a)) == a` is guaranteed for Python > `float`s. When trying the same thing with lower-precision Python floats, I > found this identity not quite fulfilled: > ``` > import numpy > b = numpy.array([1.0 / 3.0], dtype=np.float16) > float(repr(b[0])) - b[0] > Out[12]: -1.953125093259e-06 > ``` > Indeed, > ``` > b > Out[6]: array([ 0.33325195], dtype=float16) > ``` > ``` > repr(b[0]) > Out[7]: '0.33325' > ``` > When counting the bits, a float16 should hold 4.8 decimal digits, so > `repr()` seems right. Where does the garbage tail -1.953125093259e-06 > come from though? > Even more troubling, the high precision numpy types - np.longdouble and its complex version - lose intimation when used with repr. The basic problem is (roughly) that all floating-point numbers are converted to python floats before printing. I put some effort into cleaning this up, but the code is messy (actually there are several independent code paths for converting numbers to strings) and the algorithms python uses to make repr work out nicely are nontrivial. Anne > Cheers, > Nico > ___ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > https://mail.scipy.org/mailman/listinfo/numpy-discussion > ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Question about numpy.random.choice with probabilties

On Mon, Jan 23, 2017 at 3:34 PM Robert Kernwrote: > I don't object to some Notes, but I would probably phrase it more like we > are providing the standard definition of the jargon term "sampling without > replacement" in the case of non-uniform probabilities. To my mind (or more > accurately, with my background), "replace=False" obviously picks out the > implemented procedure, and I would have been incredibly surprised if it did > anything else. If the option were named "unique=True", then I would have > needed some more documentation to let me know exactly how it was > implemented. > It is what I would have expected too, but we have a concrete example of a user who expected otherwise; where one user speaks up, there are probably more who didn't (some of whom probably have code that's not doing what they think it does). So for the cost of adding a Note, why not help some of them? As for the standardness of the definition: I don't know, have you a reference where it is defined? More natural to me would be to have a list of items with integer multiplicities (as in: "cat" 3 times, "dog" 1 time). I'm hesitant to claim ours is a standard definition unless it's in a textbook somewhere. But I don't insist on my phrasing. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Question about numpy.random.choice with probabilties

On Wed, Jan 18, 2017 at 4:13 PM Nadav Har'Elwrote: > On Wed, Jan 18, 2017 at 4:30 PM, wrote: > > > > Having more sampling schemes would be useful, but it's not possible to > implement sampling schemes with impossible properties. > > > > BTW: sampling 3 out of 3 without replacement is even worse > > No matter what sampling scheme and what selection probabilities we use, we > always have every element with probability 1 in the sample. > > > I agree. The random-sample function of the type I envisioned will be able > to reproduce the desired probabilities in some cases (like the example I > gave) but not in others. Because doing this correctly involves a set of n > linear equations in comb(n,k) variables, it can have no solution, or many > solutions, depending on the n and k, and the desired probabilities. A > function of this sort could return an error if it can't achieve the desired > probabilities. > It seems to me that the basic problem here is that the numpy.random.choice docstring fails to explain what the function actually does when called with weights and without replacement. Clearly there are different expectations; I think numpy.random.choice chose one that is easy to explain and implement but not necessarily what everyone expects. So the docstring should be clarified. Perhaps a Notes section: When numpy.random.choice is called with replace=False and non-uniform probabilities, the resulting distribution of samples is not obvious. numpy.random.choice effectively follows the procedure: when choosing the kth element in a set, the probability of element i occurring is p[i] divided by the total probability of all not-yet-chosen (and therefore eligible) elements. This approach is always possible as long as the sample size is no larger than the population, but it means that the probability that element i occurs in the sample is not exactly p[i]. Anne > ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Fast vectorized arithmetic with ~32 significant digits under Numpy

On Fri, Dec 11, 2015, 18:04 David Cournapeau <courn...@gmail.com> wrote: On Fri, Dec 11, 2015 at 4:22 PM, Anne Archibald <archib...@astron.nl> wrote: Actually, GCC implements 128-bit floats in software and provides them as __float128; there are also quad-precision versions of the usual functions. The Intel compiler provides this as well, I think, but I don't think Microsoft compilers do. A portable quad-precision library might be less painful. The cleanest way to add extended precision to numpy is by adding a C-implemented dtype. This can be done in an extension module; see the quaternion and half-precision modules online. We actually used __float128 dtype as an example of how to create a custom dtype for a numpy C tutorial we did w/ Stefan Van der Walt a few years ago at SciPy. IIRC, one of the issue to make it more than a PoC was that numpy hardcoded things like long double being the higest precision, etc... But that may has been fixed since then. I did some work on numpy's long-double support, partly to better understand what would be needed to make quads work. The main obstacle is, I think, the same: python floats are only 64-bit, and many functions are stuck passing through them. It takes a lot of fiddling to make string conversions work without passing through python floats, for example, and it takes some care to produce scalars of the appropriate type. There are a few places where you'd want to modify the guts of numpy if you had a higher precision available than long doubles. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Fast vectorized arithmetic with ~32 significant digits under Numpy

Actually, GCC implements 128-bit floats in software and provides them as __float128; there are also quad-precision versions of the usual functions. The Intel compiler provides this as well, I think, but I don't think Microsoft compilers do. A portable quad-precision library might be less painful. The cleanest way to add extended precision to numpy is by adding a C-implemented dtype. This can be done in an extension module; see the quaternion and half-precision modules online. Anne On Fri, Dec 11, 2015, 16:46 Charles R Harriswrote: > On Fri, Dec 11, 2015 at 6:25 AM, Thomas Baruchel wrote: > >> From time to time it is asked on forums how to extend precision of >> computation on Numpy array. The most common answer >> given to this question is: use the dtype=object with some arbitrary >> precision module like mpmath or gmpy. >> See >> http://stackoverflow.com/questions/6876377/numpy-arbitrary-precision-linear-algebra >> or >> http://stackoverflow.com/questions/21165745/precision-loss-numpy-mpmath >> or >> http://stackoverflow.com/questions/15307589/numpy-array-with-mpz-mpfr-values >> >> While this is obviously the most relevant answer for many users because >> it will allow them to use Numpy arrays exactly >> as they would have used them with native types, the wrong thing is that >> from some point of view "true" vectorization >> will be lost. >> >> With years I got very familiar with the extended double-double type which >> has (for usual architectures) about 32 accurate >> digits with faster arithmetic than "arbitrary precision types". I even >> used it for research purpose in number theory and >> I got convinced that it is a very wonderful type as long as such >> precision is suitable. >> >> I often implemented it partially under Numpy, most of the time by trying >> to vectorize at a low-level the libqd library. >> >> But I recently thought that a very nice and portable way of implementing >> it under Numpy would be to use the existing layer >> of vectorization on floats for computing the arithmetic operations by >> "columns containing half of the numbers" rather than >> by "full numbers". As a proof of concept I wrote the following file: >> https://gist.github.com/baruchel/c86ed748939534d8910d >> >> I converted and vectorized the Algol 60 codes from >> http://szmoore.net/ipdf/documents/references/dekker1971afloating.pdf >> (Dekker, 1971). >> >> A test is provided at the end; for inverting 100,000 numbers, my type is >> about 3 or 4 times faster than GMPY and almost >> 50 times faster than MPmath. It should be even faster for some other >> operations since I had to create another np.ones >> array for testing this type because inversion isn't implemented here >> (which could of course be done). You can run this file by yourself >> (maybe you will have to discard mpmath or gmpy if you don't have it). >> >> I would like to discuss about the way to make available something related >> to that. >> >> a) Would it be relevant to include that in Numpy ? (I would think to some >> "contribution"-tool rather than including it in >> the core of Numpy because it would be painful to code all ufuncs; on the >> other hand I am pretty sure that many would be happy >> to perform several arithmetic operations by knowing that they can't use >> cos/sin/etc. on this type; in other words, I am not >> sure it would be a good idea to embed it as an every-day type but I think >> it would be nice to have it quickly available >> in some way). If you agree with that, in which way should I code it (the >> current link only is a "proof of concept"; I would >> be very happy to code it in some cleaner way)? >> >> b) Do you think such attempt should remain something external to Numpy >> itself and be released on my Github account without being >> integrated to Numpy? >> > > I think astropy does something similar for time and dates. There has also > been some talk of adding a user type for ieee 128 bit doubles. I've looked > once for relevant code for the latter and, IIRC, the available packages > were GPL :(. > > Chuck > ___ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > https://mail.scipy.org/mailman/listinfo/numpy-discussion > ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Sign of NaN

IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point arrays. Why should it be different for object arrays? Anne P.S. If you want exceptions when NaNs appear, that's what np.seterr is for. -A On Tue, Sep 29, 2015 at 5:18 PM Freddy Rietdijkwrote: > I wouldn't know of any valid output when applying the sign function to > NaN. Therefore, I think it is correct to return a ValueError. Furthermore, > I would prefer such an error over just returning NaN since it helps you > locating where NaN is generated. > > On Tue, Sep 29, 2015 at 5:13 PM, Charles R Harris < > charlesr.har...@gmail.com> wrote: > >> Hi All, >> >> Due to a recent commit, Numpy master now raises an error when applying >> the sign function to an object array containing NaN. Other options may be >> preferable, returning NaN for instance, so I would like to open the topic >> for discussion on the list. >> >> Thoughts? >> >> Chuck >> >> ___ >> NumPy-Discussion mailing list >> NumPy-Discussion@scipy.org >> https://mail.scipy.org/mailman/listinfo/numpy-discussion >> >> > ___ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > https://mail.scipy.org/mailman/listinfo/numpy-discussion > ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Python needs goto

goto! and comefrom! Together with exceptions, threads, lambda, super, generators, and coroutines, all we're lacking is call-with-current-continuation for the full list of impenetrable control-flow constructs. Oh, and lisp-style resumable exception handling. (Suggested syntax: drop(exception, value) to return control to where the exception was raised and make the raise statement return value.) On Thu, Sep 24, 2015 at 8:42 PM Charles R Harriswrote: > On Thu, Sep 24, 2015 at 12:13 PM, Yarko Tymciurak > wrote: > >> >> >> I think there are more valid uses - I've read that "goto" basically is >> what a state machine does. >> Have a read of the brief implementation notes for "goto" in golang, for >> example. Goto may not be unreasonable to use, just most people would >> abuse. Sort of like "everyone shouldn't write assembly, but if you >> understand the machine, you can make good things happen". Without >> compiler/interpreter checks, more responsibility rests on the coder to keep >> out of trouble. >> > > I would agree about state machines. When implemented using the standard > control flow constructs they always look a bit artificial. > > That depends what your "standard" control flow constructs are. Has anyone tried implementing a state machine using coroutines? They seem like a rather natural setup: each state is a coroutine that loops, doing the appropriate actions for the state and then handing control over to the coroutine for the next state. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] feature request - increment counter on write check

On Fri, Sep 11, 2015 at 3:20 PM Sebastian Bergwrote: > On Fr, 2015-09-11 at 13:10 +, Daniel Manson wrote: > > Originally posted as issue 6301 on github. > > > > > > Presumably any block of code that modifies an ndarray's buffer is > > wrapped in a (thread safe?) check of the writable flag. Would it be > > possible to hold a counter rather than a simple bool flag and then > > increment the counter whenever you test the flag? Hopefully this would > > introduce only a tiny additional overhead, but would permit > > pseudo-hashing to test whether a mutable array has changed since you > > last encountered it. > > > > Just a quick note. This is a bit more complex then it might appear. The > reason being that when a view of the array is changed, you would have to > "notify" the array itself that it has changed. So propagation from top > to bottom does not seem straight forward to me. (the other way is fine, > since on check you could check all parents, but you cannot check all > children). > Actually not so much. Like the writable flag, you'd make the counter be a per-buffer piece of information. Each array already has a pointer to the array object that "owns" the buffer, so you'd just go there in one hop. This does mean that modifying one view would affect the modified flag on all views sharing the same buffer, whether there's data overlap or not, but for caching purposes that's not so bad. I think a more serious concern is that it may be customary to simply check the writable flag by hand rather than calling an is_writable function, so to make this idea work you'd have to change all code that checks the writable flag, including user code. You'd also have to make sure that all code that tries to write to an array really checks the writable flag. Rather than making this happen for all arrays, does it make sense to use an array subclass with a "dirty flag", maybe even if this requires manual setting in some cases? Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Defining a white noise process using numpy

On Thu, Aug 27, 2015 at 12:51 AM Daniel Bliss daniel.p.bl...@gmail.com wrote: Can anyone give me some advice for translating this equation into code using numpy? eta(t) = lim(dt - 0) N(0, 1/sqrt(dt)), where N(a, b) is a Gaussian random variable of mean a and variance b**2. This is a heuristic definition of a white noise process. This is an abstract definition. How to express it in numpy will depend on what you want to do with it. The easiest and most likely thing you could want would be a time series, with N time steps dt, in which sample i is the average value of the white noise process from i*dt to (i+1)*dt. This is very easy to write in numpy: 1/np.sqrt(dt) * np.random.randn(N) Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### [Numpy-discussion] Development workflow (not git tutorial)

Hi, What is a sensible way to work on (modify, compile, and test) numpy? There is documentation about contributing to numpy at: http://docs.scipy.org/doc/numpy-dev/dev/index.html and: http://docs.scipy.org/doc/numpy-dev/dev/gitwash/development_workflow.html but these are entirely focused on using git. I have no problem with that aspect. It is building and testing that I am looking for the Right Way to do. My current approach is to build an empty virtualenv, pip install nose, and from the numpy root directory do python setup.py build_ext --inplace and python -c 'import numpy; numpy.test()'. This works, for my stock system python, though I get a lot of weird messages suggesting distutils problems (for example python setup.py develop, although suggested by setup.py itself, claims that develop is not a command). But I don't know how (for example) to test with python3 without starting from a separate clean source tree. What do you recommend: use virtualenvs? Is building inplace the way to go? Is there a better way to run all tests? Are there other packages that should go into the virtualenv? What is the best way to test on multiple python versions? Switch cleanly between feature branches? Surely I can't be the only person wishing for advice on a sensible way to work with an in-development version of numpy? Perhaps this would be a good addition to CONTRIBUTING.md or the website? Thanks, Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] DEP: Deprecate boolean array indices with non-matching shape #4353

On Fri, Jun 5, 2015 at 5:45 PM Sebastian Berg sebast...@sipsolutions.net wrote: On Fr, 2015-06-05 at 08:36 -0400, josef.p...@gmail.com wrote: snip What is actually being deprecated? It looks like there are different examples. wrong length: Nathaniels first example above, where the mask is not broadcastable to original array because mask is longer or shorter than shape[axis]. I also wouldn't have expected this to work, although I use np.nozero and boolean mask indexing interchangeably, I would assume we need the correct length for the mask. For the moment we are only talking about wrong length (along a given dimension). Not about wrong number of dimensions or multiple boolean indices. I am pro-deprecation then, definitely. I don't see a use case for padding a wrong-shaped boolean array with Falses, and the padding has burned me in the past. It's not orthogonal to the wrong-number-of-dimensions issue, though, because if your Boolean array has a dimension of length 1, broadcasting says duplicate it along that axis to match the indexee, and wrong-length says pad it with Falses. This ambiguity/pitfall disappears if the padding never happens, and that kind of broadcasting is very useful. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Backwards-incompatible improvements to numpy.random.RandomState

Do we want a deprecation-like approach, so that eventually people who want replicability will specify versions, and everyone else gets bug fixes and improvements? This would presumably take several major versions, but it might avoid people getting unintentionally trapped on this version. Incidentally, bug fixes are complicated: if a bug fix uses more or fewer raw random numbers, it breaks repeatability not just for the call that got fixed but for all successive random number generations. Anne On Sun, May 24, 2015 at 5:04 PM josef.p...@gmail.com wrote: On Sun, May 24, 2015 at 9:08 AM, Alan G Isaac alan.is...@gmail.com wrote: On 5/24/2015 8:47 AM, Ralf Gommers wrote: Values only change if you leave out the call to seed() OK, but this claim seems to conflict with the following language: the global RandomState object should use the latest implementation of the methods. I take it that this is what Nathan meant by I think this is just a bug in the description of the proposal here, not in the proposal itself. So, is the correct phrasing the global RandomState object should use the latest implementation of the methods, unless explicitly seeded? that's how I understand it. I don't see any problems with the clarified proposal for the use cases that I know of. Can we choose the version also for the global random state, for example to fix both version and seed in unit tests, with version 0? BTW: I would expect that bug fixes are still exempt from backwards compatibility. fixing #5851 should be independent of the version, (without having looked at the issue) (If you need to replicate bugs, then use an old version of a package.) Josef Thanks, Alan ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] supporting quad precision

Looking at the rational module, I think you're right: it really shouldn't be too hard to get quads working as a user type using gcc's __float128 type, which will provide hardware arithmetic in the unlikely case that the user has hardware quads. Alternatively, probably more work, one could use a package like qd to provide portable quad precision (and quad-double). I'll take a look. Anne On Jun 5, 2013 7:10 PM, David Cournapeau courn...@gmail.com wrote: On Wed, Jun 5, 2013 at 5:21 PM, Charles R Harris charlesr.har...@gmail.com wrote: Hi Anne, Long time no see ;) On Wed, Jun 5, 2013 at 10:07 AM, Anne Archibald archib...@astron.nl wrote: Hi folks, I recently came across an application I needed quad precision for (high-accuracy solution of a differential equation). I found a C++ library (odeint) that worked for the integration itself, but unfortunately it appears numpy is not able to work with quad precision arrays. For my application the quad precision is only really needed for integrator state, so I can manually convert my data to long doubles as I go to and from numpy, but it's a pain. So quad precision support would be nice. There's a thread discussing quad support: http://mail.scipy.org/pipermail/numpy-discussion/2012-February/061080.html Essentially, there isn't any, but since gcc = 4.6 supports them on Intel hardware (in software), it should be possible. (Then the thread got bogged down in bike-shedding about what to call them.) What would it take to support quads in numpy? I looked into the numpy base dtype definitions, and it's a complex arrangement involving detection of platform support and templatized C code; in the end I couldn't really see where to start. But it seems to me all the basics are available: native C syntax for basic arithmetic, qabs-style versions of all the basic functions, NaNs and Infs. So how would one go about adding quad support? There are some improvements for user types committed in 1.8-dev. Perhaps quad support could be added as a user type as it is still platform/compiler dependent. The rational type added to numpy could supply a template for adding the new type. I would be in support of that direction as well: let it live separately until CPU/compiler support is coming up. Anne, will you be at scipy conference ? Improving user data type internal API is something I'd like to work on as well David ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### [Numpy-discussion] supporting quad precision

Hi folks, I recently came across an application I needed quad precision for (high-accuracy solution of a differential equation). I found a C++ library (odeint) that worked for the integration itself, but unfortunately it appears numpy is not able to work with quad precision arrays. For my application the quad precision is only really needed for integrator state, so I can manually convert my data to long doubles as I go to and from numpy, but it's a pain. So quad precision support would be nice. There's a thread discussing quad support: http://mail.scipy.org/pipermail/numpy-discussion/2012-February/061080.html Essentially, there isn't any, but since gcc = 4.6 supports them on Intel hardware (in software), it should be possible. (Then the thread got bogged down in bike-shedding about what to call them.) What would it take to support quads in numpy? I looked into the numpy base dtype definitions, and it's a complex arrangement involving detection of platform support and templatized C code; in the end I couldn't really see where to start. But it seems to me all the basics are available: native C syntax for basic arithmetic, qabs-style versions of all the basic functions, NaNs and Infs. So how would one go about adding quad support? Thanks, Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Efficient way to load a 1Gb file?

There was also some work on a semi-mutable array type that allowed appending along one axis, then 'freezing' to yield a normal numpy array (unfortunately I'm not sure how to find it in the mailing list archives). One could write such a setup by hand, using mmap() or realloc(), but I'd be inclined to simply write a filter that converted the text file to some sort of binary file on the fly, value by value. Then the file can be loaded in or mmap()ed. A 1 Gb text file is a miserable object anyway, so it might be desirable to convert to (say) HDF5 and then throw away the text file. Anne On 10 August 2011 15:43, Derek Homeier de...@astro.physik.uni-goettingen.de wrote: On 10 Aug 2011, at 19:22, Russell E. Owen wrote: A coworker is trying to load a 1Gb text data file into a numpy array using numpy.loadtxt, but he says it is using up all of his machine's 6Gb of RAM. Is there a more efficient way to read such text data files? The npyio routines (loadtxt as well as genfromtxt) first read in the entire data as lists, which creates of course significant overhead, but is not easy to circumvent, since numpy arrays are immutable - so you have to first store the numbers in some kind of mutable object. One could write a custom parser that tries to be somewhat more efficient, e.g. first reading in sub-arrays from a smaller buffer. Concatenating those sub-arrays would still require about twice the memory of the final array. I don't know if using the array.array type (which is mutable) is much more efficient than a list... To really avoid any excess memory usage you'd have to know the total data size in advance - either by reading in the file in a first pass to count the rows, or explicitly specifying it to a custom reader. Basically, assuming a completely regular file without missing values etc., you could then read in the data like X = np.zeros((n_lines, n_columns), dtype=float) delimiter = ' ' for n, line in enumerate(file(fname, 'r')): X[n] = np.array(line.split(delimiter), dtype=float) (adjust delimiter and dtype as needed...) HTH, Derek ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] [SciPy-User] recommendation for saving data

In astronomy we tend to use FITS, which is well-supported by pyfits, but a little limited. Some new instruments are beginning to use HDF5. All these generic formats allow very general data storage, so you will need to come up with a standrdized way to represent your own data. Used well, these formats can be self-describing enough that generic tools can be very useful (e.g. display images, build histograms) but it takes some thought when designing files. Anne On 8/1/11, Christopher Barker chris.bar...@noaa.gov wrote: On 7/31/11 5:48 AM, Brian Blais wrote: I was wondering if there are any recommendations for formats for saving scientific data. every field has it's own standards -- I'd try to find one that is likely to be used by folks that may care about your results. For Oceanographic and Atmospheric modeling data, netcdf is a good option. I like the NetCDF4 python lib: http://code.google.com/p/netcdf4-python/ (there are others) For broader use, and a bit more flexibility, HDF is a good option. There are at least two ways to use it with numpy: PyTables: http://www.pytables.org (Nice higher-level interface) hf5py: http://alfven.org/wp/hdf5-for-python/ (a more raw HDF5 wrapper) There is also the npz format, built in to numpy, if you are happy with requiring python to read the data. -Chris I am running a simulation, which has many somewhat-indepedent parts which have their own internal state and parameters. I've been using pickle (gzipped) to save the entire object (which contains subobjects, etc...), but it is getting too unwieldy and I think it is time to look for a more robust solution. Ideally I'd like to have something where I can call a save method on the simulation object, and it will call the save methods on all the children, on down the line all saving into one file. It'd also be nice if it were cross-platform, and I could depend on the files being readable into the future for a while. Are there any good standards for this? What do you use for saving scientific data? thank you, Brian Blais -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/ORR(206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception chris.bar...@noaa.gov ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- Sent from my mobile device ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Can I index array starting with 1?

Don't forget the everything-looks-like-a-nail approach: make all your arrays one bigger than you need and ignore element zero. Anne On 7/28/11, Stéfan van der Walt ste...@sun.ac.za wrote: Hi Jeremy On Thu, Jul 28, 2011 at 3:19 PM, Jeremy Conlin jlcon...@gmail.com wrote: I have a need to index my array(s) starting with a 1 instead of a 0. The reason for this is to be consistent with the documentation of a format I'm accessing. I know I can just '-1' from what the documentation says, but that can get cumbersome. Is there a magic flag I can pass to a numpy array (maybe when it is created) so I can index starting at 1 instead of the Python default? You may want to have a look at some of the labeled array packages out there, such as larry, datarray, pandas, etc. I'm sure at least one of them allows integer re-labelling, although I'm not certain whether it can be done in a programmatic fashion. An alternative may be to create an indexing function that remaps the input space, e.g.: def ix(n): return n - 1 x[ix(3), ix(5)] But that looks pretty nasty, and you'll have to expand ix quite a bit to handle slices, etc. :/ Stéfan ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- Sent from my mobile device ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Can I index array starting with 1?

The can is open and the worms are everywhere, so: The big problem with one-based indexing for numpy is interpretation. In python indexing, -1 is the last element of the array, and ranges have a specific meaning. In a hypothetical one-based indexing scheme, would the last element be element 0? if not, what does looking up zero do? What about ranges - do ranges still include the first endpoint and not the second? I suppose one could choose the most pythonic of the 1-based conventions, but do any of them provide from-the-end indexing without special syntax? Once one had decided what to do, implementation would be pretty easy - just make a subclass of ndarray that replaces the indexing function. Anne On 28 July 2011 19:26, Derek Homeier de...@astro.physik.uni-goettingen.de wrote: On 29.07.2011, at 1:19AM, Stéfan van der Walt wrote: On Thu, Jul 28, 2011 at 4:10 PM, Anne Archibald aarch...@physics.mcgill.ca wrote: Don't forget the everything-looks-like-a-nail approach: make all your arrays one bigger than you need and ignore element zero. Hehe, why didn't I think of that :) I guess the kind of problem I struggle with more frequently is books written with summations over -m to +n. In those cases, it's often convenient to use the mapping function, so that I can enter the formulas as they occur. I don't want to open any cans of worms at this point, but given that Fortran90 supports such indexing (arbitrary limits, including negative ones), there definitely are use cases for it (or rather, instances where it is very convenient at least, like in Stéfan's books). So I am wondering how much it would take to implement such an enhancement for the standard ndarray... Cheers, Derek ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Quaternion dtype for NumPy - initial implementation available

What a useful package! Apart from helping all the people who know they need quaternions, this package removes one major family of use cases for vectorized small-matrix operations, namely, 3D rotations. Quaternions are the canonical way to represent orientation and rotation in three dimensions, and their multiplication gives (with some fiddling) composition of rotations. The next interesting question is, how well does scipy.interpolate deal with them? For really good rotational paths I seem to recall you want specialized splines, but simply interpolating in the quaternion domain is not a bad quick and dirty approach. Anne (now awaiting octonions, though I've never heard of a practical use for them) On 16 July 2011 10:50, Martin Ling martin-nu...@earth.li wrote: Hi all, I have just pushed a package to GitHub which adds a quaternion dtype to NumPy: https://github.com/martinling/numpy_quaternion Some backstory: on Wednesday I gave a talk at SciPy 2011 about an inertial sensing simulation package I have been working on (http://www.imusim.org/). One component I suggested might be reusable from that code was the quaternion math implementation, written in Cython. One of its features is a wrapper class for Nx4 NumPy arrays that supports efficient operations using arrays of quaternion values. Travis Oliphant suggested that a quaternion dtype would be a better solution, and got me talking to Mark Weibe about this. With Mark's help I completed this initial version at yesterday's sprint session. Incidentally, how to do something like this isn't well documented and I would have had little hope without both Mark's in-person help and his previous code (for adding a half-precision float dtype) to refer to. I don't know what the consensus is about whether people writing custom dtypes is a desirable thing, but if it is then the process needs to be made a lot easier. That said, the fact this is doable without patching the numpy core at all is really, really nice. Example usage: import numpy as np import quaternion np.quaternion(1,0,0,0) quaternion(1, 0, 0, 0) q1 = np.quaternion(1,2,3,4) q2 = np.quaternion(5,6,7,8) q1 * q2 quaternion(-60, 12, 30, 24) a = np.array([q1, q2]) a array([quaternion(1, 2, 3, 4), quaternion(5, 6, 7, 8)], dtype=quaternion) exp(a) array([quaternion(1.69392, -0.78956, -1.18434, -1.57912), quaternion(138.909, -25.6861, -29.9671, -34.2481)], dtype=quaternion) The following ufuncs are implemented: add, subtract, multiply, divide, log, exp, power, negative, conjugate, copysign, equal, not_equal, less, less_equal, isnan, isinf, isfinite, absolute Quaternion components are stored as doubles. The package could be extended to support e.g. qfloat, qdouble, qlongdouble Comparison operations follow the same lexicographic ordering as tuples. The unary tests isnan, isinf and isfinite return true if they would return true for any individual component. Real types may be cast to quaternions, giving quaternions with zero for all three imaginary components. Complex types may also be cast to quaternions, with their single imaginary component becoming the first imaginary component of the quaternion. Quaternions may not be cast to real or complex types. Comments very welcome. This is my first attempt at NumPy hacking :-) Martin ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Numeric integration of higher order integrals

When the dimensionality gets high, grid methods like you're describing start to be a problem (the curse of dimensionality). The standard approaches are simple Monte Carlo integration or its refinements (Metropolis-Hasings, for example). These converge somewhat slowly, but are not much affected by the dimensionality. Anne On 1 June 2011 05:44, Mario Bettenbuehl bette...@uni-potsdam.de wrote: Hello everyone, I am currently tackling the issue to numerically solve an integral of higher dimensions numerically. I am comparing models and their dimension increase with 2^n order. Taking a closer look to its projections along the axes, down to a two dimensions picture, the projections are of Gaussian nature, thus they show a Gaussian bump. I already used to approaches: 1. brute force: Process the values at discrete grid points and calculate the area of the obtained rectangle, cube, ... with a grid of 5x5x5x5 for a 4th order equation. 2. Gaussian quad: Cascading Gaussian quadrature given from numpy/ scipy with a grid size of 100x100x... The problem I have: For 1: How reliable are the results and does anyone have experience with equations whose projections are Gaussian like and solved these with the straight-forward-method? But how large should the grid be. For 2: How large do I need to choose the grid to still obtain reliable results? Is a grid of 10x10 sufficiently large? Thanks in advance for any reply. If needed, I'll directly provide further informations about the problem. Greetings, Mario ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] RFC: Detecting array changes (NumPy 2.0?)

On 11 March 2011 15:34, Charles R Harris charlesr.har...@gmail.com wrote: On Fri, Mar 11, 2011 at 1:06 PM, Dag Sverre Seljebotn d.s.seljeb...@astro.uio.no wrote: On Fri, 11 Mar 2011 19:37:42 + (UTC), Pauli Virtanen p...@iki.fi wrote: On Fri, 11 Mar 2011 11:47:58 -0700, Charles R Harris wrote: [clip] What about views? Wouldn't it be easier to write another object wrapping an ndarray? I think the buffer interfaces and all other various ways Numpy provides exports for arrays make keeping tabs on modification impossible to do completely reliably. Not to mention all the pain of making sure the arrays are wrapped and stay wrapped in the first place. In particular in combination with other array wrappers. I wasn't saying this is absolutely needed, just that it'd be a really convenient feature helpful for caching. Sometimes, introducing fast caching this way can remove a lot of logic from the code. Introducing a Python-space visible wrapper object kind of defeats the purpose for me. Well, starting with a wrapped object would allow you to experiment and discover what it is you really need. A smallish specialized object is probably a better starting point for development than a big solution. Operating systems do this sort of thing with the VM, but they have hardware assistance down at the lowest level and rather extensive structures to track status. Furthermore, the memory is organized into blocks and that makes it a lot easier to monitor than strided memory. In fact, I think you might want to set up your own memory subsystem and have the arrays sit on top of that. In fact, on many systems, using malloc on large contiguous blocks of memory returns a freshly-mmaped region. It's possible that with a little deviousness (and, sadly, some system-specific code) one could arrange to allocate some arrays in a way that would trigger modification-count updating by the VM system. If you're serious about detecting modifications, this sort of thing may be the only way to go - a modification-detection system that misses some modifications might be worse than none at all. An internal numpy setup is going to be a nightmare even if all you have to worry about is views and you're willing to allow non-overlapping views to count as modifying each other - you'd have to add a modification count to the ultimate base array (the one whose deletion triggers disposal of the memory arena), and then every modification to a view would have to walk the linked list of views all the way up to the top to increment the modification counter. You'll also be triggering increments of the modification counter on all sorts of non-modifications that occur in C code. Doable, but a huge job for dubious benefit. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] ragged array implimentation

On 7 March 2011 15:29, Christopher Barker chris.bar...@noaa.gov wrote: On 3/7/11 11:18 AM, Francesc Alted wrote: but, instead of returning a numpy array of 'object' elements, plain python lists are returned instead. which gives you the append option -- I can see how that would be usefull. Though I'd kind of like to have numpy ufunc/broadcasting performance. i.e.: vlarray * some_constant Be fast, and work out of the box! How about this? Represent them as object arrays (or lists) of variable-length 1d arrays. But create the 1d arrays by building one large 1d array with them all concatenated, then taking slices. Keep the single array around, perhaps in an attribute of the object array; after all, most of the numpy operations are elementwise and don't care about array structure. You'd want to be able to convert between a flat array plus length vector and object array in some reasonably efficient manner, and I think this allows that. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] sample without replacement

I know this question came up on the mailing list some time ago (19/09/2008), and the conclusion was that yes, you can do it more or less efficiently in pure python; the trick is to use two different methods. If your sample is more than, say, a quarter the size of the set you're drawing from, you permute the set and take the first few. If your sample is smaller, you draw with replacement, then redraw the duplicates, and repeat until there aren't any more duplicates. Since you only do this when your sample is much smaller than the population you don't need to repeat many times. Here's the code I posted to the previous discussion (not tested this time around) with comments: ''' def choose_without_replacement(m,n,repeats=None): Choose n nonnegative integers less than m without replacement Returns an array of shape n, or (n,repeats). if repeats is None: r = 1 else: r = repeats if nm: raise ValueError, Cannot find %d nonnegative integers less than %d % (n,m) if nm/2: res = np.sort(np.random.rand(m,r).argsort(axis=0)[:n,:],axis=0) else: res = np.random.random_integers(m,size=(n,r)) while True: res = np.sort(res,axis=0) w = np.nonzero(np.diff(res,axis=0)==0) nr = len(w[0]) if nr==0: break res[w] = np.random.random_integers(m,size=nr) if repeats is None: return res[:,0] else: return res For really large values of repeats it does too much sorting; I didn't have the energy to make it pull all the ones with repeats to the beginning so that only they need to be re-sorted the next time through. Still, the expected number of trips through the while loop grows only logarithmically with repeats, so it shouldn't be too bad. ''' Anne On 20 December 2010 12:13, John Salvatier jsalv...@u.washington.edu wrote: I think this is not possible to do efficiently with just numpy. If you want to do this efficiently, I wrote a no-replacement sampler in Cython some time ago (below). I hearby release it to the public domain. ''' Created on Oct 24, 2009 http://stackoverflow.com/questions/311703/algorithm-for-sampling-without-replacement @author: johnsalvatier ''' from __future__ import division import numpy def random_no_replace(sampleSize, populationSize, numSamples): samples = numpy.zeros((numSamples, sampleSize),dtype=int) # Use Knuth's variable names cdef int n = sampleSize cdef int N = populationSize cdef i = 0 cdef int t = 0 # total input records dealt with cdef int m = 0 # number of items selected so far cdef double u while i numSamples: t = 0 m = 0 while m n : u = numpy.random.uniform() # call a uniform(0,1) random number generator if (N - t)*u = n - m : t += 1 else: samples[i,m] = t t += 1 m += 1 i += 1 return samples On Mon, Dec 20, 2010 at 8:28 AM, Alan G Isaac alan.is...@gmail.com wrote: I want to sample *without* replacement from a vector (as with Python's random.sample). I don't see a direct replacement for this, and I don't want to carry two PRNG's around. Is the best way something like this? permutation(myvector)[:samplesize] Thanks, Alan Isaac ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Solving Ax = b: inverse vs cholesky factorization

On 8 November 2010 14:38, Joon groups.and.li...@gmail.com wrote: Oh I see. So I guess in invA = solve(Ax, I) and then x = dot(invA, b) case, there are more places where numerical errors occur, than just x = solve(Ax, b) case. That's the heart of the matter, but one can be more specific. You can think of a matrix by how it acts on vectors. Taking the inverse amounts to solving Ax=b for all the standard basis vectors (0,...,0,1,0,...,0); multiplying by the inverse amounts to expressing your vector in terms of these, finding where they go, and adding them together. But it can happen that when you break your vector up like that, the images of the components are large but almost cancel. This sort of near-cancellation amplifies numerical errors tremendously. In comparison, solving directly, if you're using a stable algorithm, is able to avoid ever constructing these nearly-cancelling combinations explicitly. The standard reason for trying to construct an inverse is that you want to solve equations for many vectors with the same matrix. But most solution methods are implemented as a matrix factorization followed by a single cheap operation, so if this is your goal, it's better to simply keep the matrix factorization around. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Stacking a 2d array onto a 3d array

On 26 October 2010 21:02, Dewald Pieterse dewald.piete...@gmail.com wrote: I see my slicing was the problem, np.vstack((test[:1], test)) works perfectly. Yes and no. np.newaxis (or None for short) is a very useful tool; you just stick it in an index expression and it adds an axis of length one there. If what you really wanted to do was pull out one plane of the array, then indexing with a number was the right thing to do. If you want to stack that plane back on the array, just add an axis of length one to it. S = A[1,...] A = np.vstack((S[np.newaxis,...],A)) As a side note, np.newaxis is actually just None, but I find the longer name much clearer, so I try to use it in my own code, and I always use it in examples I'm showing other people. I suppose a third option would be import numpy.newaxis as na. Anne On Wed, Oct 27, 2010 at 12:55 AM, josef.p...@gmail.com wrote: On Tue, Oct 26, 2010 at 8:15 PM, Dewald Pieterse dewald.piete...@gmail.com wrote: 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? keep the first dimension or add it back in test = np.array([[[1,1,1],[1,1,1]],[[2,2,2],[2,2,2]],[[3,3,3],[3,3,3]]]) np.vstack((test[:1], test)) array([[[1, 1, 1], [1, 1, 1]], [[1, 1, 1], [1, 1, 1]], [[2, 2, 2], [2, 2, 2]], [[3, 3, 3], [3, 3, 3]]]) np.vstack((test[0][None,...], test)) array([[[1, 1, 1], [1, 1, 1]], [[1, 1, 1], [1, 1, 1]], [[2, 2, 2], [2, 2, 2]], [[3, 3, 3], [3, 3, 3]]]) np.vstack((test[0][None,:,:], test)) array([[[1, 1, 1], [1, 1, 1]], [[1, 1, 1], [1, 1, 1]], [[2, 2, 2], [2, 2, 2]], [[3, 3, 3], [3, 3, 3]]]) I like expand_dims for arbitrary axis, e.g. ax=1 np.concatenate((np.expand_dims(test[:,0,:],ax), test), ax) array([[[1, 1, 1], [1, 1, 1], [1, 1, 1]], [[2, 2, 2], [2, 2, 2], [2, 2, 2]], [[3, 3, 3], [3, 3, 3], [3, 3, 3]]]) Josef -- Dewald Pieterse ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- Dewald Pieterse A democracy is nothing more than mob rule, where fifty-one percent of the people take away the rights of the other forty-nine. ~ Thomas Jefferson ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Assigning complex value to real array

On 7 October 2010 13:01, Pauli Virtanen p...@iki.fi wrote: to, 2010-10-07 kello 12:08 -0400, Andrew P. Mullhaupt kirjoitti: [clip] No. You can define the arrays as backed by mapped files with real and imaginary parts separated. Then the imaginary part, being initially zero, is a sparse part of the file, takes only a fraction of the space (and, on decent machine doesn't incur memory bandwidth costs either). You can then slipstream the cost of testing for whether the imaginary part has been subsequently assigned to zero (so you can re-sparsify the representation of a page) with any operation that examines all the values on that page. Consistency would be provided by the OS, so there wouldn't really be much numpy-specific code involved. So there is at least one efficient way to implement my suggestion. Interesting idea. Most OSes offer also page-allocated memory not backed in files. In fact, Glibc's malloc works just like this on Linux for large memory blocks. It would work automatically like this with complex arrays, if the imaginary part was stored after the real part, and additional branches were added to not write zeros to memory. But to implement this, you'd have to rewrite large parts of Numpy since the separated storage of re/im conflicts with its memory model. I believe this will simply not be done, since there seems to be little need for such a feature. Years ago MATLAB did just this - store real and complex parts of arrays separately (maybe it still does, I haven't used it in a long time). It caused us terrible performance headaches, since it meant that individual complex values were spread over two different memory areas (parts of the disk in fact, since we were using gigabyte arrays in 1996), so that writing operations in the natural way tended to cause disk thrashing. As for what's right for numpy, well, I think it makes a lot more sense to simply raise an exception when assigning a complex value to a real array (or store a NaN). Usually when this happens it means you either didn't know how numpy worked or you were feeding bogus values to some special function so it went complex on you. If you actually wanted potentially-complex values, you'd create complex arrays; as the OP says, there's little extra cost. Fortunately, this latter strategy is the way that numpy is already headed; currently I believe it emits a warning, and ISTR this is intended to be strengthened to an exception or NaN soon. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Assigning complex value to real array

On 7 October 2010 19:46, Andrew P. Mullhaupt d...@zen-pharaohs.com wrote: It wouldn't be the first time I suggested rewriting the select and choose operations. I spent months trying to get Guido to let anything more than slice indexing in arrays. And now, in the technologically advanced future, we can index a numpy array with a list, not just slices. I'm not claiming any credit for that though - I don't know who actually got that ball over the finish line. Someone who wrote code instead of complaining. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] A proposed change to rollaxis() behavior for negative 'start' values

On 23 September 2010 02:20, Ralf Gommers ralf.gomm...@googlemail.com wrote: On Wed, Sep 22, 2010 at 4:14 AM, Anne Archibald aarch...@physics.mcgill.ca wrote: Hi Ken, This is a tricky one. The current behaviour of rollaxis is to remove the requested axis from the list of axes and then insert it before the axis specified. This is exactly how python's list insertion works: In [1]: a = range(10) In [3]: a.insert(-1,'a') In [4]: a Out[4]: [0, 1, 2, 3, 4, 5, 6, 7, 8, 'a', 9] And indeed, there's no clean way to add something to the end of a list using insert (apart from the obvious a.insert(len(a),'b') ). For this you have .append(). Unfortunately numpy's rollaxis, while it agrees with insert in its behaviour, doesn't have a move_axis_to_end. The situation is also somewhat muddied by the fact that rollaxis also removes the axis from the original list of axes, so that the interpretation of index numbers is a little more subtle. But I think your suggested behaviour would be confusing because of the conflict with python's insert. How about allowing the string end as an argument to rollaxis to specify that the axis should go at the end? Allowing end is an easy solution, but note that moving an axis to the end is already possible: a = np.ones((3,4,5,6)) np.rollaxis(a, 2, len(a)+1).shape # roll axis to to last position (3, 4, 6, 5) Not consistent with insert though, there you would use len(a) instead of len(a)+1. It's a little ugly, but perhaps just documenting this is no worse than allowing a string or adding yet another function. Just a quick correction: len(a) gives a.shape[0], while what you want is actually len(a.shape). So: In [1]: a = np.zeros((2,3,4,5,6)) In [2]: len(a) Out[2]: 2 In [8]: np.rollaxis(a,0,len(a.shape)).shape Out[8]: (3, 4, 5, 6, 2) So it behaves just like insert. But len(a.shape) is rather cumbersome, especially if you haven't given a a name yet: d = (a+b*c).rollaxis(2,'end') Anne Ralf Anne On 21 September 2010 15:48, Ken Basye kbas...@jhu.edu wrote: Hi Numpy Folks, A while back, I filed this ticket: http://projects.scipy.org/numpy/ticket/1441 suggesting a change to rollaxis() and some fixes to the doc and error reporting. Ralf Gommers suggested I float the behavior change here, so that's what I'm doing. The motivation for the change comes because it seems like there should be a simpler way to get some axis into the last position than to do this: a = np.ones((3,4,5,6)) b = np.rollaxis(a, axis=0, start=len(a.shape)) b.shape (4, 5, 6, 3) But currently it seems there isn't because when you specify -1 as the 'start' argument, the axis is moved into the second-to-last position. My proposed change, which you can see on the ticket, would change that so that using -1 referred to the end position. Note that the use of negative 'start' arguments isn't currently documented and, in its current form, doesn't seem very useful. My proposal wouldn't change the behavior for positive 'start' values at all, and the interpretation of 'axis' arguments is also unaffected. If that's going to break too much code, here's a pathway that might be acceptable: Add a new function moveaxis() which works the way rollaxis() does for positive arguments but in the new way for negative arguments. Eventually, rollaxis could be deprecated to keep things tidy. This has the added advantage of using a name that seems to fit what the function does better - 'rollaxis' suggests a behavior like the roll() function which affects other axes, which isn't what happens. Thanks for listening; I'm a big fan of Numpy. Best, Ken Basye ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] slicing / indexing question

On 21 September 2010 19:20, Timothy W. Hilton hil...@meteo.psu.edu wrote: I have an 80x1200x1200 nd.array of floats this_par. I have a 1200x1200 boolean array idx, and an 80-element float array pars. For each element of idx that is True, I wish to replace the corresponding 80x1x1 slice of this_par with the elements of pars. I've tried lots of variations on the theme of this_par[idx[np.newaxis, ...]] = pars[:, np.newaxis, np.newaxis] but so far, no dice. How about this? In [1]: A = np.zeros((2,3,5)) In [2]: B = np.array([1,2]) In [3]: C = np.zeros((3,5), dtype=np.bool) In [4]: C[1,1] = True In [5]: C[2,3] = True In [6]: C Out[6]: array([[False, False, False, False, False], [False, True, False, False, False], [False, False, False, True, False]], dtype=bool) In [7]: A[:,C] = B[:,np.newaxis] In [8]: A Out[8]: array([[[ 0., 0., 0., 0., 0.], [ 0., 1., 0., 0., 0.], [ 0., 0., 0., 1., 0.]], [[ 0., 0., 0., 0., 0.], [ 0., 2., 0., 0., 0.], [ 0., 0., 0., 2., 0.]]]) The key is that indexing with C replaces the two axes C is indexing with only one; boolean indexing necessarily flattens the relevant axes. You can check this with (e.g.) A[:,C].shape. Be careful with these mixed indexing modes (partly fancy indexing, partly slicing) as they can sometimes seem to reorder your axes for you: In [16]: np.zeros((2,3,7))[:,np.ones(5,dtype=int),np.ones(5,dtype=int)].shape Out[16]: (2, 5) In [17]: np.zeros((2,3,7))[np.ones(5,dtype=int),:,np.ones(5,dtype=int)].shape Out[17]: (5, 3) In [18]: np.zeros((2,3,7))[np.ones(5,dtype=int),np.ones(5,dtype=int),:].shape Out[18]: (5, 7) Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Indexing and lookup issues

On 22 September 2010 16:38, Ross Williamson rosswilliamson@gmail.com wrote: Hi everyone I suspect this is easy but I'm stuck say I have a 1D array: t = [10,11,12] and a 2D array: id = [[0,1,0] [0,2,0] [2,0,2]] In could in IDL do y = t[id] which would produce: y = [[10,11,10] [10,12,10] [12,10,12]] i.e. use the indexes in id on the lookup array t. Is there an easy way to do this in numpy? In [1]: t = np.array([10,11,12]) In [2]: id = np.array([[0,1,0], [0,2,0], [2,0,2]]) In [3]: t[id] Out[3]: array([[10, 11, 10], [10, 12, 10], [12, 10, 12]]) Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] A proposed change to rollaxis() behavior for negative 'start' values

Hi Ken, This is a tricky one. The current behaviour of rollaxis is to remove the requested axis from the list of axes and then insert it before the axis specified. This is exactly how python's list insertion works: In [1]: a = range(10) In [3]: a.insert(-1,'a') In [4]: a Out[4]: [0, 1, 2, 3, 4, 5, 6, 7, 8, 'a', 9] And indeed, there's no clean way to add something to the end of a list using insert (apart from the obvious a.insert(len(a),'b') ). For this you have .append(). Unfortunately numpy's rollaxis, while it agrees with insert in its behaviour, doesn't have a move_axis_to_end. The situation is also somewhat muddied by the fact that rollaxis also removes the axis from the original list of axes, so that the interpretation of index numbers is a little more subtle. But I think your suggested behaviour would be confusing because of the conflict with python's insert. How about allowing the string end as an argument to rollaxis to specify that the axis should go at the end? Anne On 21 September 2010 15:48, Ken Basye kbas...@jhu.edu wrote: Hi Numpy Folks, A while back, I filed this ticket: http://projects.scipy.org/numpy/ticket/1441 suggesting a change to rollaxis() and some fixes to the doc and error reporting. Ralf Gommers suggested I float the behavior change here, so that's what I'm doing. The motivation for the change comes because it seems like there should be a simpler way to get some axis into the last position than to do this: a = np.ones((3,4,5,6)) b = np.rollaxis(a, axis=0, start=len(a.shape)) b.shape (4, 5, 6, 3) But currently it seems there isn't because when you specify -1 as the 'start' argument, the axis is moved into the second-to-last position. My proposed change, which you can see on the ticket, would change that so that using -1 referred to the end position. Note that the use of negative 'start' arguments isn't currently documented and, in its current form, doesn't seem very useful. My proposal wouldn't change the behavior for positive 'start' values at all, and the interpretation of 'axis' arguments is also unaffected. If that's going to break too much code, here's a pathway that might be acceptable: Add a new function moveaxis() which works the way rollaxis() does for positive arguments but in the new way for negative arguments. Eventually, rollaxis could be deprecated to keep things tidy. This has the added advantage of using a name that seems to fit what the function does better - 'rollaxis' suggests a behavior like the roll() function which affects other axes, which isn't what happens. Thanks for listening; I'm a big fan of Numpy. Best, Ken Basye ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] restrictions on fancy indexing

On 17 September 2010 13:47, Neal Becker ndbeck...@gmail.com wrote: It's nice I can do: f = np.linspace (0, 1, 100) u[f.1] = 0 cool, this seems to work also: u[np.abs(f).1] = 0 cool! But exactly what kind of expressions are possible here? Certainly not arbitrary code. The short answer is, anything that yields a boolean or integer array. There's no syntactical magic here. It might be clearer to write it as: c = np.abs(f).1 u[c] = 0 As for what generates boolean arrays, well, they're just numpy arrays, you can mangle them any way you want. But in particular, == != are operators that take two arrays and yield a boolean array. Also useful are ~ | and , which are the logical operators on boolean arrays. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] indexing with booleans without making a copy?

2010/9/8 Ernest Adrogué eadro...@gmx.net: I have a sorted, flat array: In [139]: a =np.array([0,1,2,2,2,3]) Basically, I want views of the areas where there are repeated numbers (since the array is sorted, they will be contiguous). But, of course, to find the numbers that are repeated I have to use comparison operations that return boolean arrays, so I suppose the problem is converting the boolean array into a slice. Well, you're going to have to do some allocation, but how's this? Use unique1d to get an array of unique values. Then use searchsorted twice, once to find the leftmost end of each hunk, and once to find the rightmost end of each hunk. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] IEEE 754-2008 decimal floating point support

On 8 September 2010 16:33, Robert Kern robert.k...@gmail.com wrote: On Wed, Sep 8, 2010 at 15:10, Michael Gilbert michael.s.gilb...@gmail.com wrote: On Wed, 8 Sep 2010 15:04:17 -0500, Robert Kern wrote: On Wed, Sep 8, 2010 at 14:44, Michael Gilbert michael.s.gilb...@gmail.com wrote: Just wanted to say that numpy object arrays + decimal solved all of my problems, which were all caused by the disconnect between decimal and binary representation of floating point numbers. Are you sure? Unless if I'm failing to think through this properly, catastrophic cancellation for large numbers is an intrinsic property of fixed-precision floating point regardless of the base. decimal and mpmath both help with that problem because they have arbitrary precision. Here is an example: 0.3/3.0 - 0.1 -1.3877787807814457e-17 mpmath.mpf( '0.3' )/mpmath.mpf( '3.0' ) - mpmath.mpf( '0.1' ) mpf('-1.3877787807814457e-17') decimal.Decimal( '0.3' )/decimal.Decimal( '3.0' ) - decimal.Decimal ( '0.1' ) Decimal(0.0) Decimal solves the problem; whereas mpmath doesn't. Okay, that's not an example of catastrophic cancellation, just a representation issue. Indeed - and as I understand it the motivation for decimal numbers is not that they suffer less from roundoff than binary, but because the round-off they suffer from is better suited to (for example) financial applications, where representing exact decimal numbers can be important. If your problem is the fact of roundoff error, try using as your test case, taking the square root of two and squaring it again. This will suffer from the same sort of roundoff problems in decimal as binary. Anne -- Robert Kern I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth. -- Umberto Eco ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] test if two arrays share the same data

2010/9/5 Ernest Adrogué eadro...@gmx.net: 5/09/10 @ 15:59 (-0500), thus spake Robert Kern: 2010/9/5 Ernest Adrogué eadro...@gmx.net: 5/09/10 @ 21:25 (+0200), thus spake Gael Varoquaux: On Sun, Sep 05, 2010 at 09:12:34PM +0200, Ernest Adrogué wrote: Hi, How can it be done? np.may_share_memory Thanks Gael and Puneeth. I think the .base attribute is enough for what I want. No, you really want may_share_memory(). Ok, I trust you :) Just to elaborate, one of the problems is that .base gives only the direct ancestor of the current array: In [2]: a = np.arange(10) In [3]: b = a[:] In [4]: c = b[:] In [5]: c.base Out[5]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) In [6]: c.base is b Out[6]: True In [7]: c.base is a Out[7]: False To check whether the arrays use the same memory arena, you have to walk the chain all the way back to the original base array. may_share_memory does this, and checks whether the arrays use overlapping ranges of memory as well. (It's may because you can have things like one being the even elements and the other being the odd elements, but checking for this is highly nontrivial, so may_share_memory just returns True.) Anne P.S. if you're thinking that this definition of base can cause memory leaks, then yes, you're right, but they are leaks only of the array descriptor objects, not the underlying memory. Still, you can exhaust memory doing things like: a = np.arange(10) while True: a = a[:] So don't do that. -A Ernest ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Array slices and number of dimensions

On 1 September 2010 17:54, Thomas Robitaille thomas.robitai...@gmail.com wrote: Hi, I'm trying to extract sub-sections of a multidimensional array while keeping the number of dimensions the same. If I just select a specific element along a given direction, then the number of dimensions goes down by one: import numpy as np a = np.zeros((10,10,10)) a.shape (10, 10, 10) a[0,:,:].shape (10, 10) This makes sense to me. If I want to retain the initial number of dimensions, I can do a[[0],:,:].shape (1, 10, 10) However, if I try and do this along two directions, I do get a reduction in the number of dimensions: a[[0],:,[5]].shape (1, 10) I'm wondering if this is normal, or is a bug? In fact, I can get what I want by doing: a[[0],:,:][:,:,[5]].shape (1, 10, 1) so I can get around the issue, but just wanted to check whether the issue with a[[0],:,[5]] is a bug? No, it's not a bug. The key problem is that supplying lists does not extract a slice - it uses fancy indexing. This implies, among other things, that the data must be copied. When you supply two lists, that means something very different in fancy indexing. When you are supplying arrays in all index slots, what you get back has the same shape as the arrays you put in; so if you supply one-dimensional lists, like A[[1,2,3],[1,4,5],[7,6,2]] what you get is [A[1,1,7], A[2,4,6], A[3,5,2]] When you supply slices in some slots, what you get is complicated, and maybe not well-defined. In particular, I think the fancy-indexing dimensions always wind up at the front, and any slice dimensions are left at the end. In short, fancy indexing is not the way to go with your problem. I generally use np.newaxis: a[7,np.newaxis,:,8,np.newaxis] but you can also use slices of length one: a[7:8, :, 8:9] Anne Thanks, Tom ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] inversion of large matrices

Hi Melissa, On 30 August 2010 17:42, Melissa Mendonça meliss...@gmail.com wrote: I've been lurking for a while here but never really introduced myself. I'm a mathematician in Brazil working with optimization and numerical analysis and I'm looking into scipy/numpy basically because I want to ditch matlab. Welcome to the list! I hope you will find the switch to numpy and scipy rewarding above and beyond not-being-matlab. Please feel free to ask questions on the list; as you've probably seen, we get lots of questions with simple but non-obvious answers, and a few really tough questions. I'm just curious as to why you say scipy.linalg.solve(), NOT numpy.linalg.solve(). Can you explain the reason for this? I find myself looking for information such as this on the internet but I rarely find real documentation for these things, and I seem to have so many performance issues with python... I'm curious to see what I'm missing here. I agree that the documentation is a little hard to find one's way around sometimes. The numpy documentation project has done a wonderful job of providing detailed documentation for all the functions in numpy, but there's not nearly as much documentation giving a general orientation (which functions are better, what the relationship is with scipy and matplotlib). The scipy documentation project is unfortunately still getting started. What sorts of performance issues have you had? python has some important limitations, but there are often ways to work around them. Sometimes there isn't an efficient way to do things, though, and in those cases we'd like to think about whether numpy/scipy can be improved. (Bulk linear algebra - solving large numbers of small problems - stands out as an example.) Numerical optimization is another - we know the optimizers built into scipy have serious limitations, and welcome ideas on improving them. Thanks, and sorry if I hijacked the thread, No problem, and welcome again to the list. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Boolean arrays

On 27 August 2010 16:17, Robert Kern robert.k...@gmail.com wrote: On Fri, Aug 27, 2010 at 15:10, Ken Watford kwatford+sc...@gmail.com wrote: On Fri, Aug 27, 2010 at 3:58 PM, Brett Olsen brett.ol...@gmail.com wrote: Hello, I have an array of non-numeric data, and I want to create a boolean array denoting whether each element in this array is a valid value or not. This is straightforward if there's only one possible valid value: import numpy as N ar = N.array((a, b, c, b, b, a, d, c, a)) ar == a array([ True, False, False, False, False, True, False, False, True], dtype=bool) If there's multiple possible valid values, I've come up with a couple possible methods, but they all seem to be inefficient or kludges: valid = N.array((a, c)) (ar == valid[0]) | (ar == valid[1]) array([ True, False, True, False, False, True, False, True, True], dtype=bool) N.array(map(lambda x: x in valid, ar)) array([ True, False, True, False, False, True, False, True, True], dtype=bool) Is there a numpy-appropriate way to do this? Thanks, Brett Olsen amap: Like Map, but for arrays. ar = numpy.array((a, b, c, b, b, a, d, c, a)) valid = ('a', 'c') numpy.amap(lambda x: x in valid, ar) array([ True, False, True, False, False, True, False, True, True], dtype=bool) I'm not sure what version of numpy this would be in; I've never seen it. But in any case, that would be very slow for large arrays since it would invoke a Python function call for every value in ar. Instead, iterate over the valid array, which is much shorter: mask = np.zeros(ar.shape, dtype=bool) for good in valid: mask |= (ar == good) Wrap that up into a function and you're good to go. That's about as efficient as it gets unless if the valid array gets large. The problem here is really one of how you specify which values are valid. If your only specification is with a python function, then you're stuck calling that python function once for each possible value, no way around it. But it could happen that you have an array of possible values and a corresponding boolean array that says whether they're valid or not. Then there's a shortcut that's probably faster than oring as Robert suggests: In [3]: A = np.array([1,2,6,4,4,2,1,7,8,2,2,1]) In [4]: B = np.unique1d(A) In [5]: B Out[5]: array([1, 2, 4, 6, 7, 8]) Here C specifies which ones are valid. C could be computed using some sort of validity function (which it may be possible to vectorize). In any case it's only the distinct values, and they're sorted (so you can use ranges). In [6]: C = np.array([True,True,True,False,False,True]) Now to compute validity of A: In [10]: C[np.searchsorted(B,A)] Out[10]: array([ True, True, False, True, True, True, True, False, True, True, True, True], dtype=bool) Anne -- Robert Kern I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth. -- Umberto Eco ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] np.asfortranarray: unnecessary copying?

This seems to me to be a bug, or rather, two bugs. 1D arrays are automatically Fortran-ordered, so isfortran should return True for them (incidentally, the documentation should be edited to indicate that the data must also be contiguous in memory). Whether or not this change is made, there's no point in asfortranarray making a copy of a 1D array, since the copy isn't any more Fortran-ordered than the input array. Another kind of iffy case is axes of length one. These should not affect C/Fortran order, since the length of their strides doesn't matter, but they do; if you use newaxis to add an axis to an array, it's still just as C/Fortran ordered as it was, but np.isfortran reports False. (Is there a np.isc or equivalent function?) Incidentally, there is a subtle misconception in your example code: when reshaping an array, the order='F' has a different meaning. It has nothing direct to do with the memory layout; what it does is define the logical arrangement of elements used while reshaping the array. The array returned will be in C order if a copy must be made, or in whatever arbitrarily-strided order is necessary if the reshape can be done without a copy. As it happens, in your example, the latter case occurs and works out to Fortran order. Anne On 30 July 2010 13:50, Kurt Smith kwmsm...@gmail.com wrote: What are the rules for when 'np.asarray' and 'np.asfortranarray' make a copy? This makes sense to me: In [3]: carr = np.arange(3) In [6]: carr2 = np.asarray(carr) In [8]: carr2[0] = 1 In [9]: carr Out[9]: array([1, 1, 2]) No copy is made. But doing the same with a fortran array makes a copy: In [10]: farr = np.arange(3).copy('F') In [12]: farr2 = np.asfortranarray(farr) In [13]: farr2[0] = 1 In [14]: farr Out[14]: array([0, 1, 2]) Could it be a 1D thing, since it's both C contiguous F contiguous? Here's a 2D example: In [15]: f2D = np.arange(10).reshape((2,5), order='F') In [17]: f2D2 = np.asfortranarray(f2D) In [19]: f2D2[0,0] = 10 In [20]: f2D Out[20]: array([[10, 2, 4, 6, 8], [ 1, 3, 5, 7, 9]]) So it looks like np.asfortranarray makes an unnecessary copy if the array is simultaneously 1D, C contiguous and F contiguous. Coercing the array with np.atleast_2d() makes asfortranarry behave. Looking further, np.isfortran always returns false if the array is 1D, even if it's Fortran contiguous (and np.isfortran is documented as such). What is the rationale here? Is it a 'column' vs. 'row' thing? Kurt ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] 3 dim array unpacking

On 12 July 2010 13:24, K.-Michael Aye kmichael@gmail.com wrote: Dear numpy hackers, I can't find the syntax for unpacking the 3 dimensions of a rgb array. so i have a MxNx3 image array 'img' and would like to do: red, green, blue = img[magical_slicing] Which slicing magic do I need to apply? Not slicing exactly; unpacking happens along the first dimension, so you need to reorder your axes so your array is 3xMxN. np.rollaxis is handy for this, as it pulls the axis you specify to the front: red, green, blue = np.rollaxis(img,2) Anne Thanks for your help! BR, Michael ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] numpy.load raising IOError but EOFError expected

On 28 June 2010 10:52, Ruben Salvador rsalvador...@gmail.com wrote: Sorry I had no access during these days. Thanks for the answer Friedrich, I had already checked numpy.savez, but unfortunately I cannot make use of it. I don't have all the data needed to be saved at the same time...it is produced each time I run a test. I think people are uncomfortable because .npy files are not designed to contain more than one array. It's a bit like concatenating a whole lot of .PNG files together - while a good decoder could pick them apart again, it's a highly misleading file since the headers do not contain any information about all the other files. npy files are similarly self-describing, and so concatenating them is a peculiar sort of thing to do. Why not simply save a separate file each time, so that you have a directory full of files? Or, if you must have just one file, use np.savez (loading the old one each time then saving the expanded object). Come to think of it, it's possible to append files to an existing zip file without rewriting the whole thing. Does numpy.savez allow this mode? That said, good exception hygiene argues that np.load should throw EOFErrors rather than the more generic IOErrors, but I don't know how difficult this would be to achieve. Anne Thanks anyway! Any other idea why this is happening? Is it expected behavior? On Thu, Jun 24, 2010 at 7:30 PM, Friedrich Romstedt friedrichromst...@gmail.com wrote: 2010/6/23 Ruben Salvador rsalvador...@gmail.com: Therefore, is this a bug? Shouldn't EOFError be raised instead of IOError? Or am I missunderstanding something? If this is not a bug, how can I detect the EOF to stop reading (I expect a way for this to work without tweaking the code with saving first in the file the number of dumps done)? Maybe you can make use of numpy.savez, http://docs.scipy.org/doc/numpy/reference/generated/numpy.savez.html#numpy-savez . Friedrich ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- Rubén Salvador PhD student @ Centro de Electrónica Industrial (CEI) http://www.cei.upm.es Blog: http://aesatcei.wordpress.com ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Solving a NLLSQ Problem by Pieces?

The basic problem with nonlinear least squares fitting, as with other nonlinear minimization problems, is that the standard algorithms find only a local minimum. It's easy to miss the global minimum and instead settle on a local minimum that is in fact a horrible fit. To deal with this, there are global optimization techniques - simulated annealing, genetic algorithms, and what have you (including the simplest, explore the whole space with a dense enough grid then fine-tune the best one with a local optimizer) - but they're computationally very expensive and not always reliable. So when it's possible to use domain-specific knowledge to make sure what you're settling on is the real minimum, this is a better solution. In this specific case, as in many optimization problems, the problem can be viewed as a series of nested approximations. The crudest approximation might even be linear in some cases; but the point is, you solve the first approximation first because it has fewer parameters and a solution space you understand better (e.g. maybe you can be sure it only has one local minimum). Then because your approximations are by assumption nested, adding more parameters complicates the space you're solving over, but you can be reasonably confident that you're close to the right solution already. (If your approximations are orthogonal in the right sense, you can even fix the parameters from previous stages and optimize only the new ones; be careful with this, though.) This approach is a very good way to incorporate domain-specific knowledge into your code, but you need to parameterize your problem as a series of nested approximations, and if it turns out the corrections are not small you can still get yourself into trouble. (Or, for that matter, if the initial solution space is complex enough you can still get yourself in trouble. Or if you're not careful your solver can take your sensible initial guess at some stage and zip off into never-never land instead of exploring nearby.) If you're interested in how other people solve this particular problem, you could take a look at the open-source panorama stitcher hugin, which fits for a potentially very large number of parameters, including a camera model. To bring this back nearer to on-topic, you will naturally not find domain-specific knowledge built into scipy or numpy, but you will find various local and global optimizers, some of which are specialized for the case of least-squares. So if you wanted to implement this sort of thing with scipy, you could write the domain-specific code yourself and simply call into one of scipy's optimizers. You could also look at OpenOpt, a scikit containing a number of global optimizers. Anne P.S. This question would be better suited to scipy-user or astropy than numpy-discussion. -A On 26 June 2010 13:12, Wayne Watson sierra_mtnv...@sbcglobal.net wrote: The context here is astronomy and optics. The real point is in the last paragraph. I'm looking at a paper that deals with 5 NL (nonlinear) equations and 8 unknown parameters. A. a=a0+arctan((y-y0)/(x-x0) B. z=V*r+S*e**(D*r) r=sqrt((x-x0)**2+(y-y0)**2) and C. cos(z)=cos(u)*cos(z)-sin(u)*sin(ep)*cos(b) sin(a-E) = sin(b)*sin(u)/sin(z) He's trying to estimate parameters of a fisheye lens which has taken star images on a photo plate. For example, x0,y0 is the offset of the center of projection from the zenith (camera not pointing straight up in the sky.) Eq. 2 expresses some nonlinearity in the lens. a0, xo, y0, V, S, D, ep, and E are the parameters. It looks like he uses gradient descent (NLLSQ is nonlinear least squares in Subject.), and takes each equation in turn using the parameter values from the preceding one in the next, B. He provides reasonable initial estimates. A final step uses all eight parameters. He re-examines ep and E, and assigns new estimates. For all (star positions) on the photo plate, he minimizes SUM (Fi**2*Gi) using values from the step for A and B, except for x0,y0. He then does some more dithering, which I'll skip. What I've presented is probably a bit difficult to understand without a good optics understanding, but my question is something like is this commonly done to solve a system of NLLSQ? It looks a bit wild. I guess if one knows his subject well, then bringing some extra knowledge to the process helps. As I understand it, he solves parameters in A, then uses them in B, and so on. I guess that's a reasonable way to do it. -- Wayne Watson (Watson Adventures, Prop., Nevada City, CA) (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time) Obz Site: 39° 15' 7 N, 121° 2' 32 W, 2700 feet Air travel safety Plus Three/Minus 8 rule. Eighty % of crashes take place in the first 3 minutes and last 8 minutes. Survival chances are good in you're paying attention. No hard drinks prior to those periods.

### Re: [Numpy-discussion] Ufunc memory access optimization

Correct me if I'm wrong, but this code still doesn't seem to make the optimization of flattening arrays as much as possible. The array you get out of np.zeros((100,100)) can be iterated over as an array of shape (1,), which should yield very substantial speedups. Since most arrays one operates on are like this, there's potentially a large speedup here. (On the other hand, if this optimization is being done, then these tests are somewhat deceptive.) On the other hand, it seems to me there's still some question about how to optimize execution order when the ufunc is dealing with two or more arrays with different memory layouts. In such a case, which one do you reorder in favour of? Is it acceptable to return freshly-allocated arrays that are not C-contiguous? Anne On 15 June 2010 07:37, Pauli Virtanen p...@iki.fi wrote: pe, 2010-06-11 kello 10:52 +0200, Hans Meine kirjoitti: On Friday 11 June 2010 10:38:28 Pauli Virtanen wrote: [clip] I think I there was some code ready to implement this shuffling. I'll try to dig it out and implement the shuffling. That would be great! Ullrich Köthe has implemented this for our VIGRA/numpy bindings: http://tinyurl.com/fast-ufunc At the bottom you can see that he basically wraps all numpy.ufuncs he can find in the numpy top-level namespace automatically. Ok, here's the branch: http://github.com/pv/numpy-work/compare/master...feature;ufunc-memory-access-speedup Some samples: (the reported times in braces are times without the optimization) x = np.zeros([100,100]) %timeit x + x 1 loops, best of 3: 106 us (99.1 us) per loop %timeit x.T + x.T 1 loops, best of 3: 114 us (164 us) per loop %timeit x.T + x 1 loops, best of 3: 176 us (171 us) per loop x = np.zeros([100,5,5]) %timeit x.T + x.T 1 loops, best of 3: 47.7 us (61 us) per loop x = np.zeros([100,5,100]).transpose(2,0,1) %timeit np.cos(x) 100 loops, best of 3: 3.77 ms (9 ms) per loop As expected, some improvement can be seen. There's also appears to be an additional 5 us (~ 700 inner loop operations it seems) overhead coming from somewhere; perhaps this can still be reduced. -- Pauli Virtanen ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Ufunc memory access optimization

On 15 June 2010 11:16, Pauli Virtanen p...@iki.fi wrote: ti, 2010-06-15 kello 10:10 -0400, Anne Archibald kirjoitti: Correct me if I'm wrong, but this code still doesn't seem to make the optimization of flattening arrays as much as possible. The array you get out of np.zeros((100,100)) can be iterated over as an array of shape (1,), which should yield very substantial speedups. Since most arrays one operates on are like this, there's potentially a large speedup here. (On the other hand, if this optimization is being done, then these tests are somewhat deceptive.) It does perform this optimization, and unravels the loop as much as possible. If all arrays are wholly contiguous, iterators are not even used in the ufunc loop. Check the part after /* Determine how many of the trailing dimensions are contiguous */ However, in practice it seems that this typically is not a significant win -- I don't get speedups over the unoptimized numpy code even for shapes (2,)*20 where you'd think that the iterator overhead could be important: x = np.zeros((2,)*20) %timeit np.cos(x) 10 loops, best of 3: 89.9 ms (90.2 ms) per loop This is actually consistent with the result x = np.zeros((2,)*20) y = x.flatten() %timeit x+x 10 loops, best of 3: 20.9 ms per loop %timeit y+y 10 loops, best of 3: 21 ms per loop that you can probably verify with any Numpy installation. This result may actually be because the Numpy ufunc inner loops are not especially optimized -- they don't use fixed stride sizes etc., and are so not much more efficient than using array iterators. This is a bit strange. I think I'd still vote for including this optimization, since one hopes the inner loop will get faster at some point. (If nothing else, the code generator can probably be made to generate specialized loops for common cases). On the other hand, it seems to me there's still some question about how to optimize execution order when the ufunc is dealing with two or more arrays with different memory layouts. In such a case, which one do you reorder in favour of? The axis permutation is chosen so that the sum (over different arrays) of strides (for each dimension) is decreasing towards the inner loops. I'm not sure, however, that this is the optimal choice. Things may depend quite a bit on the cache architecture here. I suspect that it isn't optimal. As I understand it, the key restriction imposed by the cache architecture is that an entire cache line - 64 bytes or so - must be loaded at once. For strides that are larger than 64 bytes I suspect that it simply doesn't matter how big they are. (There may be some subtle issues with cache associativity but this would be extremely architecture-dependent.) So I would say, first do any dimensions in which some or all strides are less than 64 bytes, starting from the smallest. After that, any order you like is fine. Is it acceptable to return freshly-allocated arrays that are not C-contiguous? Yes, it returns non-C-contiguous arrays. The contiguity is determined so that the ufunc loop itself can access the memory with a single stride. This may cause some speed loss in some expressions. Perhaps there should be a subtle bias towards C-order? But I'm not sure this is worth the bother. I'm more worried this may violate some users' assumptions. If a user knows they need an array to be in C order, really they should use ascontiguousarray. But as it stands it's enough to make sure that it's newly-allocated as the result of an arithmetic expression. Incautious users could suddenly start seeing copies happening, or even transposed arrays being passed to, C and Fortran code. It's also worth exploring common usage patterns to make sure that numpy still gravitates towards using C contiguous arrays for everything. I'm imagining a user who at some point adds a transposed array to a normal array, then does a long series of computations on the result. We want as many of those operations as possible to operate on contiguous arrays, but it's possible that an out-of-order array could propagate indefinitely, forcing all loops to be done with one array having large strides, and resulting in output that is stil out-of-order. Some preference for C contiguous output is worth adding. It would also be valuable to build some kind of test harness to track the memory layout of the arrays generated in some typical calculations as well as the ordering of the loop used. More generally, this problem is exactly what ATLAS is for - finding cache-optimal orders for linear algebra. So we shouldn't expect it to be simple. Anne -- Pauli Virtanen ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy

### Re: [Numpy-discussion] NumPy re-factoring project

On 11 June 2010 11:12, Benjamin Root ben.r...@ou.edu wrote: On Fri, Jun 11, 2010 at 8:31 AM, Sturla Molden stu...@molden.no wrote: It would also make sence to evaluate expressions like y = b*x + a without a temporary array for b*x. I know roughly how to do it, but don't have time to look at it before next year. (Yes I know about numexpr, I am talking about plain Python code.) Sturla If I may chime in here with my own experience with NumPy code... I typically use older, weaker computers for my work. I am not doing real-time modeling or some other really advanced, difficult computations. For me, NumPy works fast enough, even on an EeePC. My main issue is the one given above by Sturla. I find that NumPy's memory usage can go out-of-control very easily in long mathematical expressions. With a mix of constants and large matricies, each step in the order of operations seems to take up more memory. Often, I would run into a major slow-down from trashing the swap. This is fairly trivial to get around by operating over slices of the matrices at a time, but -- to me -- all of this talk about optimizing the speed of the operations without addressing the temporaries issue is like trying to tune-up the engine of a car without bothering to take the lead weights out of the trunk. I should say, though, that I've gone through the process of removing all temporary allocation using ufunc output arguments (np.add(a,b,c)) only to discover that it didn't actually save any memory and it was slower. The nice thing about temporaries is that they're, well, temporary; they go away. On the other hand, since memory reads are very slow, optimizations that do more calculation per load/store could make a very big difference, eliminating temporaries as a side effect. Anne Just my 2 cents. Ben Root ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] numpy.savez does /not/ compress!?

On 8 June 2010 06:11, Pauli Virtanen p...@iki.fi wrote: ti, 2010-06-08 kello 12:03 +0200, Hans Meine kirjoitti: On Tuesday 08 June 2010 11:40:59 Scott Sinclair wrote: The savez docstring should probably be clarified to provide this information. I would prefer to actually offer compression to the user. Unfortunately, adding another argument to this function will never be 100% secure, since currently, all kwargs will be saved into the zip, so it could count as behaviour change. Yep, that's the only question to be resolved. I suppose compression is not so usual as a variable name, so it probably wouldn't break anyone's code. This sounds like trouble, not just now but for any future additions to the interface. Perhaps it would be better to provide a function with a different, more extensible interface? (For example, one that accepts an explicit dictionary?) I'm also a little dubious about making compression the default. np.savez provides a feature - storing multiple arrays - that is not otherwise available. I suspect many users care more about speed than size. Anne -- Pauli Virtanen ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] C vs. Fortran order -- misleading documentation?

On 8 June 2010 14:16, Eric Firing efir...@hawaii.edu wrote: On 06/08/2010 05:50 AM, Charles R Harris wrote: On Tue, Jun 8, 2010 at 9:39 AM, David Goldsmith d.l.goldsm...@gmail.com mailto:d.l.goldsm...@gmail.com wrote: On Tue, Jun 8, 2010 at 8:27 AM, Pavel Bazant maxpla...@seznam.cz mailto:maxpla...@seznam.cz wrote: Correct me if I am wrong, but the paragraph Note to those used to IDL or Fortran memory order as it relates to indexing. Numpy uses C-order indexing. That means that the last index usually (see xxx for exceptions) represents the most rapidly changing memory location, unlike Fortran or IDL, where the first index represents the most rapidly changing location in memory. This difference represents a great potential for confusion. in http://docs.scipy.org/doc/numpy/user/basics.indexing.html is quite misleading, as C-order means that the last index changes rapidly, not the memory location. Any index can change rapidly, depending on whether is in an inner loop or not. The important distinction between C and Fortran order is how indices translate to memory locations. The documentation seems correct to me, although it might make more sense to say the last index addresses a contiguous range of memory. Of course, with modern processors, actual physical memory can be mapped all over the place. Chuck To me, saying that the last index represents the most rapidly changing memory location means that if I change the last index, the memory location changes a lot, which is not true for C-order. So for C-order, supposed one scans the memory linearly (the desired scenario), it is the last *index* that changes most rapidly. The inverted picture looks like this: For C-order, changing the first index leads to the most rapid jump in *memory*. Still have the feeling the doc is very misleading at this important issue. Pavel The distinction between your two perspectives is that one is using for-loop traversal of indices, the other is using pointer-increment traversal of memory; from each of your perspectives, your conclusions are correct, but my inclination is that the pointer-increment traversal of memory perspective is closer to the spirit of the docstring, no? I think the confusion is in most rapidly changing memory location, which is kind of ambiguous because a change in the indices is always a change in memory location if one hasn't used index tricks and such. So from a time perspective it means nothing, while from a memory perspective the largest address changes come from the leftmost indices. Exactly. Rate of change with respect to what, or as you do what? I suggest something like the following wording, if you don't mind the verbosity as a means of conjuring up an image (although putting in diagrams would make it even clearer--undoubtedly there are already good illustrations somewhere on the web): Note to those used to Matlab, IDL, or Fortran memory order as it relates to indexing. Numpy uses C-order indexing by default, although a numpy array can be designated as using Fortran order. [With C-order, sequential memory locations are accessed by incrementing the last index.] For a two-dimensional array, think if it as a table. With C-order indexing the table is stored as a series of rows, so that one is reading from left to right, incrementing the column (last) index, and jumping ahead in memory to the next row by incrementing the row (first) index. With Fortran order, the table is stored as a series of columns, so one reads memory sequentially from top to bottom, incrementing the first index, and jumps ahead in memory to the next column by incrementing the last index. One more difference to be aware of: numpy, like python and C, uses zero-based indexing; Matlab, [IDL???], and Fortran start from one. - If you want to keep it short, the key wording is in the sentence in brackets, and you can chop out the table illustration. I'd just like to point out a few warnings to keep in mind while rewriting this section: Numpy arrays can have any configuration of memory strides, including some that are zero; C and Fortran contiguous arrays are simply those that have special arrangements of the strides. The actual stride values is normally almost irrelevant to python code. There is a second meaning of C and Fortran order: when you are reshaping an array, you can specify one order or the

### Re: [Numpy-discussion] C vs. Fortran order -- misleading documentation?

On 8 June 2010 17:17, David Goldsmith d.l.goldsm...@gmail.com wrote: On Tue, Jun 8, 2010 at 1:56 PM, Benjamin Root ben.r...@ou.edu wrote: On Tue, Jun 8, 2010 at 1:36 PM, Eric Firing efir...@hawaii.edu wrote: On 06/08/2010 08:16 AM, Eric Firing wrote: On 06/08/2010 05:50 AM, Charles R Harris wrote: On Tue, Jun 8, 2010 at 9:39 AM, David Goldsmithd.l.goldsm...@gmail.com mailto:d.l.goldsm...@gmail.com wrote: On Tue, Jun 8, 2010 at 8:27 AM, Pavel Bazantmaxpla...@seznam.cz mailto:maxpla...@seznam.cz wrote: Correct me if I am wrong, but the paragraph Note to those used to IDL or Fortran memory order as it relates to indexing. Numpy uses C-order indexing. That means that the last index usually (see xxx for exceptions) represents the most rapidly changing memory location, unlike Fortran or IDL, where the first index represents the most rapidly changing location in memory. This difference represents a great potential for confusion. in http://docs.scipy.org/doc/numpy/user/basics.indexing.html is quite misleading, as C-order means that the last index changes rapidly, not the memory location. Any index can change rapidly, depending on whether is in an inner loop or not. The important distinction between C and Fortran order is how indices translate to memory locations. The documentation seems correct to me, although it might make more sense to say the last index addresses a contiguous range of memory. Of course, with modern processors, actual physical memory can be mapped all over the place. Chuck To me, saying that the last index represents the most rapidly changing memory location means that if I change the last index, the memory location changes a lot, which is not true for C-order. So for C-order, supposed one scans the memory linearly (the desired scenario), it is the last *index* that changes most rapidly. The inverted picture looks like this: For C-order, changing the first index leads to the most rapid jump in *memory*. Still have the feeling the doc is very misleading at this important issue. Pavel The distinction between your two perspectives is that one is using for-loop traversal of indices, the other is using pointer-increment traversal of memory; from each of your perspectives, your conclusions are correct, but my inclination is that the pointer-increment traversal of memory perspective is closer to the spirit of the docstring, no? I think the confusion is in most rapidly changing memory location, which is kind of ambiguous because a change in the indices is always a change in memory location if one hasn't used index tricks and such. So from a time perspective it means nothing, while from a memory perspective the largest address changes come from the leftmost indices. Exactly. Rate of change with respect to what, or as you do what? I suggest something like the following wording, if you don't mind the verbosity as a means of conjuring up an image (although putting in diagrams would make it even clearer--undoubtedly there are already good illustrations somewhere on the web): Note to those used to Matlab, IDL, or Fortran memory order as it relates to indexing. Numpy uses C-order indexing by default, although a numpy array can be designated as using Fortran order. [With C-order, sequential memory locations are accessed by incrementing the last Maybe change sequential to contiguous. I was thinking maybe subsequent might be a better word. IMV, contiguous has more of a physical connotation. (That just isn't valid in Numpy, correct?) So I'd prefer subsequent as an alternative to sequential. In the end, we need to communicate this clearly. No matter which language, I have always found it difficult to get new programmers to understand the importance of knowing the difference between row-major and column-major. A thick paragraph isn't going to help to get the idea across to a person who doesn't even know that a problem exists. Maybe a car analogy would be good here... Maybe if one imagine city streets (where many of the streets are one-way), and need to drop off mail at each address. Would it be more efficient to go up and back a street or to drop off mail at the

### Re: [Numpy-discussion] Dynamic convolution in Numpy

On 6 June 2010 04:44, David Cournapeau courn...@gmail.com wrote: On Thu, Jun 3, 2010 at 7:49 PM, arthur de conihout arthurdeconih...@gmail.com wrote: I don't know if i made myself very clear. if anyone has suggestions or has already operated a dynamic filtering i would be well interested. Does fade-in/fade-out actually works ? I would have thought that it had killed the properties of your filters ? There are two issues: - how to do convolution fast - how to go from one filter to the other I think the kicker is here: what is the right way to interpolate between filters? If you have, or can generate, a lot of filters, then at least you can evaluate the quality of the interpolation. The right way to understand what kind of interpolation to is to have some idea of the physics you're dealing with. In this case, as I understand it, you're dealing with the auditory effects of the head, seen from different angles. I would say that the ear senses something integrated power in logarithmically-spaced frequency bins over roughly twentieth-of-a-second time intervals. So the relevant effects you care about are amplitude absorption and time delays, if they are as long as a twentieth of a second. Doing simple linear interpolation, unfortunately, will probably get you in trouble - imagine, for example, that two impulse responses have the same amplitude at 440 Hz but different phases. A linear interpolation will change the amplitude (for example if they're 180 degrees apart it'll pass through zero). You might do better with interpolating in polar coordinates, though you might have phase wrapping issues. A really thorough approach might be to take a granular-synthesis approach to the impulse responses, breaking them up into orthogonal time-domain channels within which the response is defined by an amplitude and a time delay, which you'd interpolate in the natural way. I'd try polar interpolation on the FFTs of the amplitudes first, though (since in fact it's the same thing with the minimal possible frequency channels). I suspect that some interpolation, however unrealistic (even linear interpolation) is necessary or listeners may perceive sounds snapping from place to place in the aural field. The main issue with changing filters is that your system is not LTI anymore. If your filters have finite impulse answer, I guess it should not be too much of an issue. To do convolution quickly, you need to use FFT, which is a bit tricky if you want to do things in real-time, as you need to partition the impulse response. Using partitioned impulse answer as keywords should give you plenty of references on how to do it, As far as convolution, as David says, take a look at existing algorithms and maybe even music software - there's a trade-off between the n^2 computation of a brute-force FIR filter and the delay introduced by an FFT approach, but on-the-fly convolution is a well-studied problem. Anne David ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Numerical Recipes (for Python)?

On 4 June 2010 00:24, Wayne Watson sierra_mtnv...@sbcglobal.net wrote: The link below leads me to http://numpy.scipy.org/, with or without the whatever. IRAF is not mentioned on the home page. Um. I was not being specific. For a concrete example of what I mean, suppose you wanted to solve an ordinary differential equation. I would recommend you read the chapter on ODEs in Numerical Recipes in (say) C. This would talk about adaptive versus fixed step sizes, how to convert higher-order ODEs into first-order ODEs, how to formulate and solve boundary value problems, and so on. It would also describe in detail one particular adaptive integrator, a Runge-Kutta 4/5 integrator. My recommendation would be to take that understanding of what integrators can and can't do and how they should be treated, and then use scipy.integrate.odeint or scipy.integrate.ode to solve your actual problem. These two packages contain careful thoroughly-tested implementations of adaptive integrators of the sort described in NR. They will correctly handle all sorts of awkward special cases, and are fairly hard to fool. If these are not sufficient (and I know their interface in scipy is not ideal) I'd recommend going to pydstool, which has a much more flexible interface, better performance, and more modern algorithms under the hood. Only in extremis would I consider implementing my own ODE solver: perhaps if I needed one with special features (a symplectic integrator, perhaps) and I couldn't find public well-tested code to do that. So: read Numerical Recipes, by all means, in any programming language you like; but use, if at all possible, existing libraries rather than implementing anything described in NR. Getting numerical code right is really hard. Let someone else do it for you. In the case of python, scipy itself is pretty much a library providing what's in NR. Anne On 6/1/2010 9:04 PM, Anne Archibald wrote: On 2 June 2010 00:33, Wayne Watsonsierra_mtnv...@sbcglobal.net wrote: Subject is a book title from some many years ago, I wonder if it ever got to Python? I know there were C and Fortran versions. There is no Numerical Recipes for python. The main reason there isn't a NR for python is that practically everything they discuss is already implemented as python libraries, and most of it is in numpy and/or scipy. (Their algorithms are also not suitable for pure-python implementation, but that's a whole other discussion.) I should also say that while NR is justifiably famous for its explanations of numerical issues, its code is not under a free license (so you may not use it without the authors' permission) and many people feel it has many bugs. The algorithms they discuss are also not always the best available. I generally recommend that people doing scientific programming read all or part of NR to understand the algorithms' limitations but then use the implementations available in numpy/scipy/scikits/IRAF/whatever. Anne -- Wayne Watson (Watson Adventures, Prop., Nevada City, CA) (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time) Obz Site: 39° 15' 7 N, 121° 2' 32 W, 2700 feet Science and democracy are based on the rejection of dogma. -- Dick Taverne, The March of Unreason Web Page:www.speckledwithstars.net/ ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- Wayne Watson (Watson Adventures, Prop., Nevada City, CA) (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time) Obz Site: 39° 15' 7 N, 121° 2' 32 W, 2700 feet Science and democracy are based on the rejection of dogma. -- Dick Taverne, The March of Unreason Web Page:www.speckledwithstars.net/ ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Numerical Recipes (for Python)?

On 4 June 2010 14:32, Wayne Watson sierra_mtnv...@sbcglobal.net wrote: At one point in my career I was very familiar, and that's an understatement :-), with many of these methods (NR and beyond). I have zero interest in implementing them.I do not need explanations of the theory behind them. What I need to know is where some of these methods exist in libraries? Optimization (linear, nonlinear), regression (multiple, stepwise, and others), matrix inverse, eigenvalues, Fourier transforms, ..., on and onI would expect to find a site that lists all of them, and I can pick the ones I need. Python of course. Read the scipy documentation. Anne On 6/3/2010 11:09 PM, Anne Archibald wrote: On 4 June 2010 00:24, Wayne Watsonsierra_mtnv...@sbcglobal.net wrote: The link below leads me to http://numpy.scipy.org/, with or without the whatever. IRAF is not mentioned on the home page. Um. I was not being specific. For a concrete example of what I mean, suppose you wanted to solve an ordinary differential equation. I would recommend you read the chapter on ODEs in Numerical Recipes in (say) C. This would talk about adaptive versus fixed step sizes, how to convert higher-order ODEs into first-order ODEs, how to formulate and solve boundary value problems, and so on. It would also describe in detail one particular adaptive integrator, a Runge-Kutta 4/5 integrator. My recommendation would be to take that understanding of what integrators can and can't do and how they should be treated, and then use scipy.integrate.odeint or scipy.integrate.ode to solve your actual problem. These two packages contain careful thoroughly-tested implementations of adaptive integrators of the sort described in NR. They will correctly handle all sorts of awkward special cases, and are fairly hard to fool. If these are not sufficient (and I know their interface in scipy is not ideal) I'd recommend going to pydstool, which has a much more flexible interface, better performance, and more modern algorithms under the hood. Only in extremis would I consider implementing my own ODE solver: perhaps if I needed one with special features (a symplectic integrator, perhaps) and I couldn't find public well-tested code to do that. So: read Numerical Recipes, by all means, in any programming language you like; but use, if at all possible, existing libraries rather than implementing anything described in NR. Getting numerical code right is really hard. Let someone else do it for you. In the case of python, scipy itself is pretty much a library providing what's in NR. Anne On 6/1/2010 9:04 PM, Anne Archibald wrote: On 2 June 2010 00:33, Wayne Watsonsierra_mtnv...@sbcglobal.netwrote: Subject is a book title from some many years ago, I wonder if it ever got to Python? I know there were C and Fortran versions. There is no Numerical Recipes for python. The main reason there isn't a NR for python is that practically everything they discuss is already implemented as python libraries, and most of it is in numpy and/or scipy. (Their algorithms are also not suitable for pure-python implementation, but that's a whole other discussion.) I should also say that while NR is justifiably famous for its explanations of numerical issues, its code is not under a free license (so you may not use it without the authors' permission) and many people feel it has many bugs. The algorithms they discuss are also not always the best available. I generally recommend that people doing scientific programming read all or part of NR to understand the algorithms' limitations but then use the implementations available in numpy/scipy/scikits/IRAF/whatever. Anne -- Wayne Watson (Watson Adventures, Prop., Nevada City, CA) (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time) Obz Site: 39° 15' 7 N, 121° 2' 32 W, 2700 feet Science and democracy are based on the rejection of dogma. -- Dick Taverne, The March of Unreason Web Page:www.speckledwithstars.net/ ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion -- Wayne Watson (Watson Adventures, Prop., Nevada City, CA) (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time) Obz Site: 39° 15' 7 N, 121° 2' 32 W, 2700 feet Science and democracy are based on the rejection of dogma. -- Dick Taverne, The March of Unreason Web Page:www.speckledwithstars.net/ ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org

### Re: [Numpy-discussion] Numerical Recipes (for Python)?

On 2 June 2010 00:33, Wayne Watson sierra_mtnv...@sbcglobal.net wrote: Subject is a book title from some many years ago, I wonder if it ever got to Python? I know there were C and Fortran versions. There is no Numerical Recipes for python. The main reason there isn't a NR for python is that practically everything they discuss is already implemented as python libraries, and most of it is in numpy and/or scipy. (Their algorithms are also not suitable for pure-python implementation, but that's a whole other discussion.) I should also say that while NR is justifiably famous for its explanations of numerical issues, its code is not under a free license (so you may not use it without the authors' permission) and many people feel it has many bugs. The algorithms they discuss are also not always the best available. I generally recommend that people doing scientific programming read all or part of NR to understand the algorithms' limitations but then use the implementations available in numpy/scipy/scikits/IRAF/whatever. Anne -- Wayne Watson (Watson Adventures, Prop., Nevada City, CA) (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time) Obz Site: 39° 15' 7 N, 121° 2' 32 W, 2700 feet Science and democracy are based on the rejection of dogma. -- Dick Taverne, The March of Unreason Web Page:www.speckledwithstars.net/ ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] ix_ and copies

On 30 May 2010 11:25, Keith Goodman kwgood...@gmail.com wrote: On Sat, May 29, 2010 at 2:49 PM, Anne Archibald aarch...@physics.mcgill.ca wrote: On 29 May 2010 15:09, Robert Kern robert.k...@gmail.com wrote: On Sat, May 29, 2010 at 12:27, Keith Goodman kwgood...@gmail.com wrote: Will making changes to arr2 never change arr1 if arr2 = arr1[np.ix_(*lists)] where lists is a list of (index) lists? np.ix_ returns a tuple of arrays so I'm guessing (and hoping) the answer is yes. Correct. To expand: any time you do fancy indexing - that is, index by anything but a tuple of integers or slice objects - you get back a copy. I have never seen such a simple and clear definition of the line between regular and fancy indexing. To make sure I understand I'll try to expand. Is the following right? Regular indexing (no copy made): int float bool slice tuple of any combination of the above I think you should not include bool on this list. Strictly speaking I believe you can use bools as if they were integers, but that's a little limited. Normally when one indexes with bools one is indexing with an array of bools, as a sort of condition index; that is fancy indexing. The underlying reason for the copy/no copy distinction is that numpy arrays must be evenly strided, that is, as you move along any axis, the space between data elements must not vary. So slicing is no problem, and supplying an integer is no problem. Supplying a float is kind of bogus but might work anyway. Supplying None or np.newaxis also works, since this just adds an axis of length one. Fancy indexing (copy made): Any indexing that is not regular indexing The only two options here are (essentially) indexing with (tuples of) arrays of indices or indexing with boolean (condition) arrays. Mixed modes, where you supply a tuple containing some arrays of indices and/or some booleans but also some slice objects or integers, may work but may do something unexpected or may simply fail to work. There was last time I looked no systematic testing of such constructions, and the implementation was erratic. (This is largely a definitional issue; given the way numpy's arrays of indices and boolean indexing work it's not clear how one should interpret such a mixed indexing operation.) There have been occasional calls for a pure-python implementation of numpy indexing for reference purposes. I think such a thing would be fun to write, but I haven't had time. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] ix_ and copies

On 29 May 2010 15:09, Robert Kern robert.k...@gmail.com wrote: On Sat, May 29, 2010 at 12:27, Keith Goodman kwgood...@gmail.com wrote: Will making changes to arr2 never change arr1 if arr2 = arr1[np.ix_(*lists)] where lists is a list of (index) lists? np.ix_ returns a tuple of arrays so I'm guessing (and hoping) the answer is yes. Correct. To expand: any time you do fancy indexing - that is, index by anything but a tuple of integers or slice objects - you get back a copy. Anne -- Robert Kern I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth. -- Umberto Eco ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Finding Star Images on a Photo (Video chip) Plate?

On 28 May 2010 21:09, Charles R Harris charlesr.har...@gmail.com wrote: On Fri, May 28, 2010 at 5:45 PM, Wayne Watson sierra_mtnv...@sbcglobal.net wrote: Suppose I have a 640x480 pixel video chip and would like to find star images on it, possible planets and the moon. A possibility of noise exits, or bright pixels. Is there a known method for finding the centroids of these astro objects? You can threshold the image and then cluster the pixels in objects. I've done this on occasion using my own software, but I think there might be something in scipy/ndimage that does the same. Someone here will know. There are sort of two passes here - the first is to find all the stars, and the second is to fine down their positions, ideally to less than a pixel. For the former, thresholding and clumping is probably the way to go. For the latter I think a standard approach is PSF fitting - that is, you fit (say) a two-dimensional Gaussian to the pixels near your star. You'll fit for at least central (subpixel) position, probably radius, and maybe eccentricity and orientation. You might even fit for a more sophisticated PSF (doughnuts are natural for Schmidt-Cassegrain telescopes, or the diffraction pattern of your spider). Any spot whose best-fit PSF is just one pixel wide is noise or a cosmic ray hit or a hotpixel; any spot whose best-fit PSF is huge is a detector splodge or a planet or galaxy. All this assumes that your CCD has more resolution than your optics; if this is not the case you're more or less stuck, since a star is then just a bright pixel. In this case your problem is one of combining multiple offset images, dark skies, and dome flats to try to distinguish detector crud and cosmic ray hits from actual stars. It can be done, but it will be a colossal pain if your pointing accuracy is not subpixel (which it probably won't be). In any case, my optical friends tell me that the Right Way to do all this is to use all the code built into IRAF (or its python wrapper, pyraf) that does all this difficult work for you. Anne P.S. if your images have been fed through JPEG or some other lossy compression the process will become an utter nightmare. -A Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] curious about how people would feel about moving to github

On 27 May 2010 04:43, Charles R Harris charlesr.har...@gmail.com wrote: Maybe most importantly, distributed revision control places any possible contributor on equal footing with those with commit access; this is one important step in making contributors feel valued. Well, not quite. They can't commit to the main repository. I think the main thing is to be responsive: fast review, quick commit. And quick to offer commit rights to anyone who sends in more that a couple of decent patches. Maybe we should take a vow to review one patch a week. Okay. Suppose we wanted to replicate the current permissions arrangement as closely as possible with git. It seems to me it would look something like this: * Set up a git repository somewhere on scipy.org. * Give everyone who currently has permission to commit to SVN permission to write to this repository. * git submissions would become possible: a user would make some changes but instead of posting a patch would link to a particular git state. The changes could be reviewed and incorporated like a patch, but with easier merging and better history. If the changes became out of date the user could easily merge from the central repository and resolve the conflict themselves. * Patch submissions would be reviewed as now and committed to git by one of the people who do this now. Alternatively they could be integrated to the mainline by someone without write access and published as a git change, to be incorporated (easily) as above by someone with write access. * if review and inclusion were slow it would nevertheless be easy for users to pull from each other and build on each other's changes without making the eventual merge a nightmare. So, no major change to who controls what. The nipy/ipython model takes this a step further, reasoning that git makes branching and merging so easy there's no need for such a large group of people with write access to the central repository, but if that doesn't work for numpy/scipy we don't need to do it. And we can change in either direction at any time with no major changes to infrastructure or workflow. To get back to the original point of the thread: nobody has yet objected to git, and all we have are some debates about the ultimate workflow that don't make much difference to whether or how git should be adopted. Is this a fair description? Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] curious about how people would feel about moving to github

Hi Jarrod, I'm in favour of the switch, though I don't use Windows. I find git far more convenient to use than SVN; I've been using git-svn, and in spite of the headaches it's caused me I still prefer it to raw SVN. It seems to me that git's flexibility in how people collaborate means we can do a certain amount of figuring out after the switch. My experience with a small project has been that anyone who wants to make major changes just clones the repository on github and makes the changes; then we email the main author to ask him to pull particular branches into the main repo. It works well enough. Anne On 26 May 2010 19:47, Jarrod Millman mill...@berkeley.edu wrote: Hello, I changed the subject line for this thread, since I didn't want to hijack another thread. Anyway, I am not proposing that we actually decide whether to move to git and github now, but I am just curious how people would feel. We had a conversation about this a few years ago and it was quite contentious at the time. Since then, I believe a number of us have started using git and github for most of our work. And there are a number of developers using git-svn to develop numpy now. So I was curious to get a feeling for what people would think about it, if we moved to git. (I don't want to rehash the arguments for the move.) Anyway, Chuck listed the main concerns we had previously when we discussed moving from svn to git. See the discussion below. Are there any other concerns? Am I right in thinking that most of the developers would prefer git at this point? Or are there still a number of developers who would prefer using svn still? On Wed, May 26, 2010 at 12:54 PM, Charles R Harris charlesr.har...@gmail.com wrote: I think the main problem has been windows compatibility. Git is best from the command line whereas the windows command line is an afterthought. Another box that needs a check-mark is the buildbot. If svn clients are supported then it may be that neither of those are going to be a problem. I was under the impression that there were a number of decent git clients for Windows now, but I don't know anyone who develops on Windows. Are there any NumPy developers who use Windows who could check out the current situation? Pulling from github with an svn client works very well, so buildbot could continue working as is: http://github.com/blog/626-announcing-svn-support And if it turns out the Windows clients are still not good enough, we could look into the recently add svn write support to github: http://github.com/blog/644-subversion-write-support No need for us to make any changes immediately. I am just curious how people would feel about it at this point. Jarrod ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] curious about how people would feel about moving to github

On 27 May 2010 01:22, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, May 26, 2010 at 9:49 PM, Jarrod Millman mill...@berkeley.edu wrote: On Wed, May 26, 2010 at 8:08 PM, Matthew Brett matthew.br...@gmail.com wrote: That's the model we've gone for in nipy and ipython too. We wrote it up in a workflow doc project. Here are the example docs giving the git workflow for ipython: https://cirl.berkeley.edu/mb312/gitwash/ and in particular: https://cirl.berkeley.edu/mb312/gitwash/development_workflow.html I would highly recommend using this workflow. Ideally, we should use the same git workflow for all the scipy-related projects. That way developers can switch between projects without having to switch workflows. The model that Matthew and Fernando developed for nipy and ipython seem like a very reasonable place to start. __ I wouldn't. Who is going to be the gate keeper and pull the stuff? No vacations for him/her, on 24 hour call, yes? They might as well run a dairy. And do we really want all pull requests cross-posted to the list? Linus works full time as gatekeeper for Linux and gets paid for the effort. I think a central repository model would work better for us. I don't think this is as big a problem as it sounds. If the gatekeeper takes a week-long vacation, so what? People keep working on their changes independently and they can get merged when the gatekeeper gets around to it. If they want to accelerate the ultimate merging they can pull the central repository into their own and resolve all conflicts, so that the pull into the central repository goes smoothly. If the gatekeeper's away and the users want to swap patches, well, they just pull from each other's public git repositories. Anne Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] curious about how people would feel about moving to github

On 27 May 2010 01:55, Matthew Brett matthew.br...@gmail.com wrote: Hi, Linux has Linus, ipython has Fernando, nipy has... well, I'm sure it is somebody. Numpy and Scipy no longer have a central figure and I like it that way. There is no reason that DVCS has to inevitably lead to a central authority. I think I was trying to say that the way it looks as if it will be - before you try it - is very different from the way it actually is when you get there. Anne put the idea very well - but I still think it is very hard to understand, without trying it, just how liberating the workflow is from anxieties about central authorities and so on. You can just get on with what you want to do, talk with or merge from whoever you want, and the whole development process becomes much more fluid and productive. And I know that sounds chaotic but - it just works. Really really well. One way to think of it is that there is no main line of development. The only time the central repository needs to pull from the others is when a release is being prepared. As it stands we do have a single release manager, though it's not necessarily the same for each version. So if we wanted, they could just go and pull and merge the repositories of everyone who's made a useful change, then release the results. Of course, this will be vastly easier if all those other people have already merged each other's results (into different branches if appropriate). But just like now, it's the release manager's decision which changes end up in the next version. This is not the only way to do git development; it's the only one I have experience with, so I can't speak for the effectiveness of others. But I have no doubt that we can find some way that works, and I don't think we necessarily need to decide what that is any time soon. Anne See you, Matthew ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] efficient way to manage a set of floats?

On 12 May 2010 20:09, Dr. Phillip M. Feldman pfeld...@verizon.net wrote: Warren Weckesser-3 wrote: A couple questions: How many floats will you be storing? When you test for membership, will you want to allow for a numerical tolerance, so that if the value 1 - 0.7 is added to the set, a test for the value 0.3 returns True? (0.3 is actually 0.2, while 1-0.7 is 0.30004) Warren Anne- Thanks for that absolutely beautiful explanation!! Warren- I had not initially thought about numerical tolerance, but this could potentially be an issue, in which case the management of the data would have to be completely different. Thanks for pointing this out! I might have as many as 50,000 values. If you want one-dimensional sets with numerical tolerances, then either a sorted-array implementation looks more appealing. A sorted-tree implementation will be a little awkward, since you will often need to explore two branches to find out the nearest neighbour of a query point. In fact what you have is a one-dimensional kd-tree, which is helpfully provided by scipy.spatial, albeit without insertion or deletion operators. I should also point out that when you start wanting approximate matches, which you will as soon as you do any sort of arithmetic on your floats, your idea of a set becomes extremely messy. For example, suppose you try to insert a float that's one part in a million different from one that's in the table. Does it get inserted too or is it equal to what's there? When it comes time to remove it, your query will probably have a value slightly different from either previous value - which one, or both, do you remove? Or do you raise an exception? Resolving these questions satisfactorily will probably require you to know the scales that are relevant in your problem and implement sensible handling of scales larger or smaller than this (but beware of the teapot in a stadium problem, of wildly different scales in the same data set). Even so, you will want to write algorithms that are robust to imprecision, duplication, and disappearance of points in your sets. (If this sounds like the voice of bitter experience, well, I discovered while writing a commercial ray-tracer that when you shoot billions of rays into millions of triangles, all sorts of astonishing limitations of floating-point turn into graphical artifacts. Which are always *highly* visible. It was during this period that the interval-arithmetic camp nearly gained a convert.) Anne Phillip -- View this message in context: http://old.nabble.com/efficient-way-to-manage-a-set-of-floats--tp28518014p28542439.html Sent from the Numpy-discussion mailing list archive at Nabble.com. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] efficient way to manage a set of floats?

On 10 May 2010 18:53, Dr. Phillip M. Feldman pfeld...@verizon.net wrote: I have an application that involves managing sets of floats. I can use Python's built-in set type, but a data structure that is optimized for fixed-size objects that can be compared without hashing should be more efficient than a more general set construct. Is something like this available? You might not find this as useful as you think - on a 32-bit machine, the space overhead is roughly a 32-bit object pointer or two for each float, plus about twice the number of floats times 32-bit pointers for the table. And hashing might be worthwhile anyway - you could easily have a series of floats with related bit patterns you'd want to scatter all over hash space. Plus python's set object has seen a fair amount of performance tweaing. That said, there's support in numpy for many operations which use a sorted 1D array to represent a set of floats. There's searchsorted for lookup, plus IIRC union and intersection operators; I'm not sure about set difference. The big thing missing is updates, though if you can batch them, concatenate followed by sorting should be reasonable. Removal can be done with fancy indexing, though again batching is recommended. Maybe these should be regarded as analogous to python's frozensets. In terms of speed, the numpy functions are obviously not as asymptotically efficient as hash tables, though I suspect memory coherence is more of an issue than O(1) versus O(log(n)). The numpy functions allow vectorized lookups and vector operations on sets, which could be handy. Anne P.S. if you want sets with fuzzy queries, it occurs to me that scipy's kdtrees will actually do an okay job, and in compiled code. No updates there either, though. -A View this message in context: http://old.nabble.com/efficient-way-to-manage-a-set-of-floats--tp28518014p28518014.html Sent from the Numpy-discussion mailing list archive at Nabble.com. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] efficient way to manage a set of floats?

On 10 May 2010 21:56, Dr. Phillip M. Feldman pfeld...@verizon.net wrote: Anne Archibald-2 wrote: on a 32-bit machine, the space overhead is roughly a 32-bit object pointer or two for each float, plus about twice the number of floats times 32-bit pointers for the table. Hello Anne, I'm a bit confused by the above. It sounds as though the hash table approach might occupy 4 times as much storage as a single array. Is that right? Probably. Hash tables usually operate at about half-full for efficiency, so you'd have twice as many entries as you do objects. If you're using a python hash table, you also have type tagging and malloc information for each object. If you had a custom implementation, you'd have to have some way to mark hash cells as empty. In any case, expect at least a doubling of the space needed for a simple array. Also, I don't understand why hashing would be useful for the set application. The reason hash tables rely on hashing is not just to obtain a number for a potentially complex object; a good hash function should produce effectively random numbers for the user's objects. These numbers are reduced modulo the size of the hash table to determine where the object should go. If the user supplies a whole bunch of objects that all happen to hash to the same value, they'll all try to go into the same bin, and the hash table degenerates to an O(n) object as it has to search through the whole list each time. If you are using floating-point objects, well, for example the exponent may well be all the same, or take only a few values. If, after reduction modulo the table size, it's what gets used to determine where your numbers go, they'll all go in the same bin. Or you could be using all integers, which usually end with lots of zeros in floating-point representation, so that all your numbers go in exactly the same bin. You could try using non-power-of-two table sizes, but that just sort of hides the problem: someone someday will provide a collection of numbers, let's say a, a+b, a+2b, ... that reduce to the same value modulo your table size, and suddenly your hash table is agonizingly slow. There's kind of an art to designing good general-purpose hash functions; it's very handy that python provides one. It seems as though a red-black tree might be a good implementation for a set of floats if all that one wants to do is add, delete, and test membership. (I will also need to do unions). If you're implementing it from scratch, you could go with a red-black tree, but a hash table is probably faster. I'd go with a simple hash table with linked-list buckets. Managing insertions and deletions should be only a minor pain compared to implementing a whole tree structure. You can probably find a nice hash function for floats with a bit of googling the literature. I should say, try just using python sets first, and only go into all this if they prove to be the slowest part of your program. Thanks for the help! Good luck, Anne Phillip -- View this message in context: http://old.nabble.com/efficient-way-to-manage-a-set-of-floats--tp28518014p28519085.html Sent from the Numpy-discussion mailing list archive at Nabble.com. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Question about numpy.arange()

On 1 May 2010 16:36, Gökhan Sever gokhanse...@gmail.com wrote: Hello, Is b an expected value? I am suspecting another floating point arithmetic issue. I[1]: a = np.arange(1.6, 1.8, 0.1, dtype='float32') I[2]: a O[2]: array([ 1.6002, 1.7005], dtype=float32) I[3]: b = np.arange(1.7, 1.8, 0.1, dtype='float32') I[4]: b O[4]: array([ 1.7005, 1.7995], dtype=float32) A bit conflicting with the np.arange docstring: Values are generated within the half-open interval ``[start, stop)`` (in other words, the interval including `start` but excluding `stop`). This is a floating-point issue; since 1.7995 does not actually equal 1.8, it is included. This arises because 0.1, 1.7, and 1.8 cannot be exactly represented in floating-point. A good rule to avoid being annoyed by this is: only use arange for integers. Use linspace if you want floating-point. Anne Thanks. -- Gökhan ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] proposing a beware of [as]matrix() warning

On 28 April 2010 14:30, Alan G Isaac ais...@american.edu wrote: On 4/28/2010 12:08 PM, Dag Sverre Seljebotn wrote: it would be good to deprecate the matrix class from NumPy Please let us not have this discussion all over again. I think you may be too late on this, but it's worth a try. The matrix class is very useful for teaching. In economics for example, the use of matrix algebra is widespread, while algebra with arrays that are not matrices is very rare. I can (and do) use NumPy matrices even in undergraduate courses. If you do not like them, do not use them. This is the problem: lots of people start using numpy and think hmm, I want to store two-dimensional data so I'll use a matrix, and have no idea that matrix means anything different from two-dimensional array. It was this that inspired David's original post, and it's this that we're trying to find a solution for. If you want `matrix` replaced with a better matrix object, offer a replacement for community consideration. Thank you, Alan Isaac PS There is one change I would not mind: let A * M be undefined if A is an ndarray and M is a NumPy matrix. I can definitely vote for this, in the interest of catching as many inadvertent matrix users as possible. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Aggregate memmap

On 21 April 2010 15:41, Matthew Turk matthewt...@gmail.com wrote: Hi there, I've quite a bit of unformatted fortran data that I'd like to use as input to a memmap, as sort of a staging area for selection of subregions to be loaded into RAM. Unfortunately, what I'm running into is that the data was output as a set of slices through a 3D cube, instead of a single 3D cube -- the end result being that each slice also contains a record delimiter. I was wondering if there's a way to either specify that every traversal through the least-frequent dimension requires an additional offset or to calculate individual offsets for each slice myself and then aggregate these into a super memmap. If you had a delimiter for each entry in a record array you could just make a dummy record field for the delimiter, but it sounds like you have one delimiter per slice instead. I think the way to go is to memmap it in as a flat uint8 array and then construct a view into that array with appropriate strides (which includes skipping the delimiter). This may well be what you meant by a super memmap. Thanks for any suggestions you might have! Good luck, Anne -Matt ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] How to combine a pair of 1D arrays?

On 14 April 2010 11:34, Robert Kern robert.k...@gmail.com wrote: On Wed, Apr 14, 2010 at 10:25, Peter Shinners p...@shinners.org wrote: Is there a way to combine two 1D arrays with the same size into a 2D array? It seems like the internal pointers and strides could be combined. My primary goal is to not make any copies of the data. There is absolutely no way to get around that, I am afraid. Well, I'm not quite sure I agree with this. The best way, of course, is to allocate the original two arrays as subarrays of one larger array, that is, start with the fused array and select your two 1D arrays as subarrays. Of course, this depends on how you're generating the 1D arrays; if they're simply being returned to you from a black-box function, you're stuck with them. But if it's a ufunc-like operation, it may have an output argument that lets you write to a supplied array rather than allocating a new one. If they're coming from a file, you may be able to map the whole file (or a large chunk) as an array and select them as subarrays (even if alignment and type are issues). You should also keep in mind that allocating arrays and copying data really isn't very expensive - malloc() is extremely fast, especially for large arrays which are just grabbed as blocks from the OS - and copying arrays is also cheap, and can be used to reorder data into a more cache-friendly order. If the problem is that your arrays almost fill available memory, you will already have noticed that using numpy is kind of a pain, because many operations involve copies. But if you really have to do this, it may be possible. numpy arrays are specified by giving a data area, and strides into that data area. The steps along each axis must be uniform, but if you really have two arrays with the same stride, you may be able to use a gruesome hack to make it work. Each of your two arrays has a data pointer, which essentially points to the first element. So if you make up your two-row array using the same data pointer as your first array, and give it a stride along the second dimension equal to the difference between pointers, this should work. Of course, you have to make sure python doesn't deallocate the second array out from under you, and you may have to defeat some error checking, but in principle it should be possible. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Find indices of largest elements

On 14 April 2010 16:56, Keith Goodman kwgood...@gmail.com wrote: On Wed, Apr 14, 2010 at 12:39 PM, Nikolaus Rath nikol...@rath.org wrote: Keith Goodman kwgood...@gmail.com writes: On Wed, Apr 14, 2010 at 8:49 AM, Keith Goodman kwgood...@gmail.com wrote: On Wed, Apr 14, 2010 at 8:16 AM, Nikolaus Rath nikol...@rath.org wrote: Hello, How do I best find out the indices of the largest x elements in an array? Example: a = [ [1,8,2], [2,1,3] ] magic_function(a, 2) == [ (0,1), (1,2) ] Since the largest 2 elements are at positions (0,1) and (1,2). Here's a quick way to rank the data if there are no ties and no NaNs: ...or if you need the indices in order: shape = (3,2) x = np.random.rand(*shape) x array([[ 0.52420123, 0.43231286], [ 0.97995333, 0.87416228], [ 0.71604075, 0.66018382]]) r = x.reshape(-1).argsort().argsort() I don't understand why this works. Why do you call argsort() twice? Doesn't that give you the indices of the sorted indices? It is confusing. Let's look at an example: x = np.random.rand(4) x array([ 0.37412289, 0.68248559, 0.12935131, 0.42510212]) If we call argsort once we get the index that will sort x: idx = x.argsort() idx array([2, 0, 3, 1]) x[idx] array([ 0.12935131, 0.37412289, 0.42510212, 0.68248559]) Notice that the first element of idx is 2. That's because element x[2] is the min of x. But that's not what we want. We want the first element to be the rank of the first element of x. So we need to shuffle idx around so that the order aligns with x. How do we do that? We sort it! Unless I am mistaken, what you are doing here is inverting the permutation returned by the first argsort. The second argsort is an n log n method, though, and permutations can be inverted in linear time: ix = np.argsort(X) ixinv = np.zeros_like(ix) ixinv[ix] = np.arange(len(ix)) This works because if ix is a permutation and ixinv is its inverse, A = B[ix] is the same as A[ixinv] = B This also means that you can often do without the actual inverse by moving the indexing operation to the other side of the equal sign. (Not in the OP's case, though.) I'll also add that if the OP needs the top m for 1mn, sorting the whole input array is not the most efficient algorithm; there are priority-queue-based schemes that are asymptotically more efficient, but none exists in numpy. Since numpy's sorting is quite fast, I personally would just use the sorting. Anne idx.argsort() array([1, 3, 0, 2]) The min value of x is x[2], that's why 2 is the first element of idx which means that we want ranked(x) to contain a 0 at position 2 which it does. Bah, it's all magic. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] binomial coefficient, factorial

On 14 April 2010 17:37, Christopher Barker chris.bar...@noaa.gov wrote: jah wrote: Is there any chance that a binomial coefficent and factorial function can make their way into NumPy? probably not -- numpy is over-populated already I know these exist in Scipy, but I don't want to have to install SciPy just to have something so basic. The problem is that everyone has a different basic. What we really need is an easier way for folks to use sub-packages of scipy. I've found myself hand-extracting just what I need, too. Maybe I've been spoiled by using Linux, but my answer to this sort of thing is always just install scipy already. If it's not already packaged for your operating system (e.g. Linux, OSX, or Windows), is it really difficult to compile it yourself? (If so, perhaps post a report to scipy-user and see if we can fix it so it isn't.) Anne -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/ORR(206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception chris.bar...@noaa.gov ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Best dtype for Boolean values

On 12 April 2010 11:59, John Jack itsmilesda...@gmail.com wrote: Hello all. I am (relatively) new to python, and 100% new to numpy. I need a way to store arrays of booleans and compare the arrays for equality. I assume I want arrays of dtype Boolean, and I should compare the arrays with array_equal tmp.all_states array([False, True, False], dtype=bool) tmp1.all_states array([False, False, False], dtype=bool) tmp1.all_states[1]=True tmp1.all_states array([False, True, False], dtype=bool) array_equal(tmp.all_states,tmp1.all_states) True any(tmp.all_states) True Would this be (a) the cheapest way (w.r.t. memory) to store Booleans and (b) the most efficient way to compare two lists of Booleans? The short answer is yes and yes. The longer answer is, that uses one byte per Boolean, which is a tradeoff. In some sense, modern machines are happier working with 32- or 64-bit quantities, so loading a one-byte Boolean requires a small amount of byte-shuffling. On the other hand, if you're really short of memory, 8 bits for a Boolean is really wasteful. In fact, since modern machines are almost always limited by memory bandwidth, a packed Boolean data structure would probably be much faster for almost all operations in spite of the bit-fiddling required. But such a representation is incompatible with the whole infrastructure of numpy and so would require a great deal of messy code to support. So yes, it's the best representation of Booleans available, unless you're dealing with mind-bogglingly large arrays of them, in which case some sort of packed-Boolean representation would be better. This can even be partially supported by numpy, using uint8s, bitwise operations, and manually-specified bitmasks. There are probably not many applications for which this is worth the pain. Anne P.S. There's actually at least one python package for bit vectors, outside numpy; I can't speak for how good it is, though. -A Thanks for your advice. -JJ. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Proposal for new ufunc functionality

On 12 April 2010 18:26, Travis Oliphant oliph...@enthought.com wrote: We should collect all of these proposals into a NEP. Or several NEPs, since I think they are quasi-orthogonal. To clarify what I mean by group-by behavior. Suppose I have an array of floats and an array of integers. Each element in the array of integers represents a region in the float array of a certain kind. The reduction should take place over like-kind values: Example: add.reduceby(array=[1,2,3,4,5,6,7,8,9], by=[0,1,0,1,2,0,0,2,2]) results in the calculations: 1 + 3 + 6 + 7 2 + 4 5 + 8 + 9 and therefore the output (notice the two arrays --- perhaps a structured array should be returned instead...) [0,1,2], [17, 6, 22] The real value is when you have tabular data and you want to do reductions in one field based on values in another field. This happens all the time in relational algebra and would be a relatively straightforward thing to support in ufuncs. As an example, if I understand correctly, this would allow the histogram functions to be replaced by a one-liner, e.g.: add.reduceby(array=1, by=((A-min)*n/(max-min)).astype(int)) It would also be valuable to support output arguments of some sort, so that, for example, reduceby could be used to accumulate values into an output array at supplied indices. (I leave the value of using this on matrix multiplication or arctan2 to the reader's imagination.) Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Proposal for new ufunc functionality

2010/4/10 Stéfan van der Walt ste...@sun.ac.za: On 10 April 2010 19:45, Pauli Virtanen p...@iki.fi wrote: Another addition to ufuncs that should be though about is specifying the Python-side interface to generalized ufuncs. This is an interesting idea; what do you have in mind? I can see two different kinds of answer to this question: one is a tool like vectorize/frompyfunc that allows construction of generalized ufuncs from python functions, and the other is thinking out what methods and support functions generalized ufuncs need. The former would be very handy for prototyping gufunc-based libraries before delving into the templated C required to make them actually efficient. The latter is more essential in the long run: it'd be nice to have a reduce-like function, but obviously only when the arity and dimensions work out right (which I think means (shape1,shape2)-(shape2) ). This could be applied along an axis or over a whole array. reduceat and the other, more sophisticated, schemes might also be worth supporting. At a more elementary level, gufunc objects should have good introspection - docstrings, shape specification accessible from python, named formal arguments, et cetera. (So should ufuncs, for that matter.) Anne Regards Stéfan ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Do ufuncs returned by frompyfunc() have the out arg?

On 6 April 2010 15:42, Ken Basye kbas...@jhu.edu wrote: Folks, I hope this is a simple question. When I created a ufunc with np.frompyfunc(), I got an error when I called the result with an 'out' argument: In fact, ordinary ufuncs do not accept names for their arguments. This is annoying, but fixing it involves rooting around in the bowels of the ufunc machinery, which are not hacker-friendly. Anne def foo(x): return x * x + 1 ufoo = np.frompyfunc(foo, 1, 1) arr = np.arange(9).reshape(3,3) ufoo(arr, out=arr) Traceback (most recent call last): File stdin, line 1, in module TypeError: 'out' is an invalid keyword to foo (vectorized) But I notice that if I just put the array there as a second argument, it seems to work: ufoo(arr, arr) array([[2, 5, 26], [101, 290, 677], [1370, 2501, 4226]], dtype=object) # and now arr is the same as the return value Is it reasonable to conclude that there is an out-arg in the resulting ufunc and I just don't know the right name for it? I also tried putting some other right-shaped array as a second argument and it did indeed get filled in. Thanks as always, Ken ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Bug in logaddexp2.reduce

On 1 April 2010 01:59, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:46 PM, Anne Archibald peridot.face...@gmail.com wrote: On 1 April 2010 01:40, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:25 PM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:22 AM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:17 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 6:08 PM, josef.p...@gmail.com wrote: On Wed, Mar 31, 2010 at 7:37 PM, Warren Weckesser warren.weckes...@enthought.com wrote: T J wrote: On Wed, Mar 31, 2010 at 1:21 PM, Charles R Harris charlesr.har...@gmail.com wrote: Looks like roundoff error. So this is expected behavior? In [1]: np.logaddexp2(-1.5849625007211563, -53.584962500721154) Out[1]: -1.5849625007211561 In [2]: np.logaddexp2(-0.5849625007211563, -53.584962500721154) Out[2]: nan Is any able to reproduce this? I don't get 'nan' in either 1.4.0 or 2.0.0.dev8313 (32 bit Mac OSX). In an earlier email T J reported using 1.5.0.dev8106. np.logaddexp2(-0.5849625007211563, -53.584962500721154) nan np.logaddexp2(-1.5849625007211563, -53.584962500721154) -1.5849625007211561 np.version.version '1.4.0' WindowsXP 32 What compiler? Mingw? yes, mingw 3.4.5. , official binaries release 1.4.0 by David sse2 Pentium M Can you try the exp2/log2 functions with the problem data and see if something goes wrong? Works fine for me. If it helps clarify things, the difference between the two problem input values is exactly 53 (and that's what logaddexp2 does an exp2 of); so I can provide a simpler example: In [23]: np.logaddexp2(0, -53) Out[23]: nan Of course, for me it fails in both orders. Ah, that's progress then ;) The effective number of bits in a double is 53 (52 + implicit bit). That shouldn't cause problems but it sure looks suspicious. Indeed, that's what led me to the totally wrong suspicion that denormals have something to do with the problem. More data points: In [38]: np.logaddexp2(63.999, 0) Out[38]: nan In [39]: np.logaddexp2(64, 0) Out[39]: 64.0 In [42]: np.logaddexp2(52.999, 0) Out[42]: 52.9990002 In [43]: np.logaddexp2(53, 0) Out[43]: nan It looks to me like perhaps the NaNs are appearing when the smaller term affects only the extra bits provided by the FPU's internal larger-than-double representation. Some such issue would explain why the problem seems to be hardware- and compiler-dependent. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Bug in logaddexp2.reduce

On 1 April 2010 02:24, Charles R Harris charlesr.har...@gmail.com wrote: On Thu, Apr 1, 2010 at 12:04 AM, Anne Archibald peridot.face...@gmail.com wrote: On 1 April 2010 01:59, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:46 PM, Anne Archibald peridot.face...@gmail.com wrote: On 1 April 2010 01:40, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:25 PM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:22 AM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:17 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 6:08 PM, josef.p...@gmail.com wrote: On Wed, Mar 31, 2010 at 7:37 PM, Warren Weckesser warren.weckes...@enthought.com wrote: T J wrote: On Wed, Mar 31, 2010 at 1:21 PM, Charles R Harris charlesr.har...@gmail.com wrote: Looks like roundoff error. So this is expected behavior? In [1]: np.logaddexp2(-1.5849625007211563, -53.584962500721154) Out[1]: -1.5849625007211561 In [2]: np.logaddexp2(-0.5849625007211563, -53.584962500721154) Out[2]: nan Is any able to reproduce this? I don't get 'nan' in either 1.4.0 or 2.0.0.dev8313 (32 bit Mac OSX). In an earlier email T J reported using 1.5.0.dev8106. np.logaddexp2(-0.5849625007211563, -53.584962500721154) nan np.logaddexp2(-1.5849625007211563, -53.584962500721154) -1.5849625007211561 np.version.version '1.4.0' WindowsXP 32 What compiler? Mingw? yes, mingw 3.4.5. , official binaries release 1.4.0 by David sse2 Pentium M Can you try the exp2/log2 functions with the problem data and see if something goes wrong? Works fine for me. If it helps clarify things, the difference between the two problem input values is exactly 53 (and that's what logaddexp2 does an exp2 of); so I can provide a simpler example: In [23]: np.logaddexp2(0, -53) Out[23]: nan Of course, for me it fails in both orders. Ah, that's progress then ;) The effective number of bits in a double is 53 (52 + implicit bit). That shouldn't cause problems but it sure looks suspicious. Indeed, that's what led me to the totally wrong suspicion that denormals have something to do with the problem. More data points: In [38]: np.logaddexp2(63.999, 0) Out[38]: nan In [39]: np.logaddexp2(64, 0) Out[39]: 64.0 In [42]: np.logaddexp2(52.999, 0) Out[42]: 52.9990002 In [43]: np.logaddexp2(53, 0) Out[43]: nan It looks to me like perhaps the NaNs are appearing when the smaller term affects only the extra bits provided by the FPU's internal larger-than-double representation. Some such issue would explain why the problem seems to be hardware- and compiler-dependent. Hmm, that 63.999 is kinda strange. Here is a bit of code that might confuse a compiler working with different size mantissas: @type@ npy_log2...@c@(@type@ x) { @type@ u = 1 + x; if (u == 1) { return LOG2E*x; } else { return npy_l...@c@(u) * x / (u - 1); } } It might be that u != 1 does not imply u-1 != 0. That does indeed look highly suspicious. I'm not entirely sure how to work around it. GSL uses a volatile declaration: http://www.google.ca/codesearch/p?hl=en#p9nGS4eQGUI/gnu/gsl/gsl-1.8.tar.gz%7C8VCQSLJ5jR8/gsl-1.8/sys/log1p.cq=log1p On the other hand boost declares itself defeated by optimizing compilers and uses a Taylor series: http://www.google.ca/codesearch/p?hl=en#sdP2GRSfgKo/dcplusplus/trunk/boost/boost/math/special_functions/log1p.hppq=log1psa=Ncd=7ct=rc While R makes no mention of the corrected formula or optimizing compilers but takes the same approach, only with Chebyshev series: http://www.google.ca/codesearch/p?hl=en#gBBSWbwZmuk/src/base/R-2/R-2.3.1.tar.gz%7CVuh8XhRbUi8/R-2.3.1/src/nmath/log1p.cq=log1p Since, at least on my machine, ordinary log1p appears to work fine, is there any reason not to have log2_1p call it and scale the result? Or does the compiler make a hash of our log1p too? Anne Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Bug in logaddexp2.reduce

On 1 April 2010 02:49, David Cournapeau da...@silveregg.co.jp wrote: Charles R Harris wrote: On Thu, Apr 1, 2010 at 12:04 AM, Anne Archibald peridot.face...@gmail.comwrote: On 1 April 2010 01:59, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:46 PM, Anne Archibald peridot.face...@gmail.com wrote: On 1 April 2010 01:40, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:25 PM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:22 AM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:17 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 6:08 PM, josef.p...@gmail.com wrote: On Wed, Mar 31, 2010 at 7:37 PM, Warren Weckesser warren.weckes...@enthought.com wrote: T J wrote: On Wed, Mar 31, 2010 at 1:21 PM, Charles R Harris charlesr.har...@gmail.com wrote: Looks like roundoff error. So this is expected behavior? In [1]: np.logaddexp2(-1.5849625007211563, -53.584962500721154) Out[1]: -1.5849625007211561 In [2]: np.logaddexp2(-0.5849625007211563, -53.584962500721154) Out[2]: nan Is any able to reproduce this? I don't get 'nan' in either 1.4.0 or 2.0.0.dev8313 (32 bit Mac OSX). In an earlier email T J reported using 1.5.0.dev8106. np.logaddexp2(-0.5849625007211563, -53.584962500721154) nan np.logaddexp2(-1.5849625007211563, -53.584962500721154) -1.5849625007211561 np.version.version '1.4.0' WindowsXP 32 What compiler? Mingw? yes, mingw 3.4.5. , official binaries release 1.4.0 by David sse2 Pentium M Can you try the exp2/log2 functions with the problem data and see if something goes wrong? Works fine for me. If it helps clarify things, the difference between the two problem input values is exactly 53 (and that's what logaddexp2 does an exp2 of); so I can provide a simpler example: In [23]: np.logaddexp2(0, -53) Out[23]: nan Of course, for me it fails in both orders. Ah, that's progress then ;) The effective number of bits in a double is 53 (52 + implicit bit). That shouldn't cause problems but it sure looks suspicious. Indeed, that's what led me to the totally wrong suspicion that denormals have something to do with the problem. More data points: In [38]: np.logaddexp2(63.999, 0) Out[38]: nan In [39]: np.logaddexp2(64, 0) Out[39]: 64.0 In [42]: np.logaddexp2(52.999, 0) Out[42]: 52.9990002 In [43]: np.logaddexp2(53, 0) Out[43]: nan It looks to me like perhaps the NaNs are appearing when the smaller term affects only the extra bits provided by the FPU's internal larger-than-double representation. Some such issue would explain why the problem seems to be hardware- and compiler-dependent. Hmm, that 63.999 is kinda strange. Here is a bit of code that might confuse a compiler working with different size mantissas: @type@ npy_log2...@c@(@type@ x) { @type@ u = 1 + x; if (u == 1) { return LOG2E*x; } else { return npy_l...@c@(u) * x / (u - 1); } } It might be that u != 1 does not imply u-1 != 0. I don't think that's true for floating point, since the spacing at 1 and at 0 are quite different, and that's what defines whether two floating point numbers are close or not. But I think you're right about the problem. I played a bit on my netbook where I have the same problem (Ubuntu 10.04 beta, gc 4.4.3, 32 bits). The problem is in npy_log2_1p (for example, try npy_log2_1p(1e-18) vs npy_log2_p1(1e-15)). If I compile with -ffloat-store, the problem also disappears. I think the fix is to do : @type@ npy_log2...@c@(@type@ x) { �...@type@ u = 1 + x; if ((u-1) == 0) { return LOG2E*x; } else { return npy_l...@c@(u) * x / (u - 1); } } Particularly given the comments in the boost source code, I'm leery of this fix; who knows what an optimizing compiler will do with it? What if, for example, u-1 and u get treated differently - one gets truncated to double, but the other doesn't, for example? Now the correction is completely bogus. Or, of course, the compiler might rewrite the expressions. At the least an explicit volatile variable or two is necessary, and even so making it compiler-proof seems like a real challenge. At the least it's letting us in for all sorts of subtle trouble down the road (though at least unit tests can catch it). Since this is a subtle issue, I vote for delegating it (and log10_1p) to log1p. If *that* has trouble, well, at least it'll only be on non-C99 machines, and we can try compiler gymnastics. this works for a separate C program, but for some reason, once I try this in numpy, it does not work, and I have no battery anymore, Clearly the optimizing compiler is inserting the DRDB (drain David's battery) opcode. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Bug in logaddexp2.reduce

On 1 April 2010 03:15, David Cournapeau da...@silveregg.co.jp wrote: Anne Archibald wrote: Particularly given the comments in the boost source code, I'm leery of this fix; who knows what an optimizing compiler will do with it? But the current code *is* wrong: it is not true that u == 1 implies u - 1 == 0 (and that (u-1) != 0 - u != 1), because the spacing between two consecutive floats is much bigger at 1 than at 0. And the current code relies on this wrong assumption: at least with the correction, we test for what we care about. I don't think this is true for IEEE floats, at least in the case we're interested in where u is approximately 1. After all, the issue is representation: if you take an x very close to zero and add 1 to it, you get exactly 1 (which is what you're pointing out); but if you then subtract 1 again, you get exactly zero. After all, u==1 means that the bit patterns for u and 1 are identical, so if you subtract 1 from u you get exactly zero. If u is not equal to 1 (and we're talking about IEEE floats) and you subtract 1, you will get something that is not zero. The problem is that what we're dealing with is not IEEE floating-point: it resembles IEEE floating-point, but arithmetic is often done with extra digits, which are discarded or not at the whim of the compiler. If the extra digits were kept we'd still have (u-1)==0 iff u==1, but if the compiler discards the extra digits in one expression but not the other - which is a question about FPU register scheduling - funny things can happen. So I would say the code correction I suggest is at least better in that respect. Now the correction is completely bogus. I am not sure to understand why ? It is at least not worse than the current code. Well, it would be hard to be worse than inappropriate NaNs. But the point of the correction is that, through some deviousness I don't fully understand, the roundoff error involved in the log is compensated for by the roundoff error in the calculation of (1+x)-1. If the compiler uses a different number of digits to calculate the log than it uses to caiculate (1+x)-1, the correction will make a hash of things. Since this is a subtle issue, I vote for delegating it (and log10_1p) to log1p. If *that* has trouble, well, at least it'll only be on non-C99 machines, and we can try compiler gymnastics. Yes, we could do that. I note that on glibc, the function called is an intrinsic for log1p (FYL2XP1) if x is sufficiently small. Even if we end up having to implement this function by hand, we should at least implement it only once - log1p - and not three times (log1p, log10_1p, log2_1p). Clearly the optimizing compiler is inserting the DRDB (drain David's battery) opcode. :) That was tongue-in-cheek, but there was a real point: optimizing compilers do all sorts of peculiar things, particularly with floating-point. For one thing, many Intel CPUs have two (logically-)separate units for doing double-precision calculations - the FPU, which keeps extra digits, and SSE2, which doesn't. Who knows which one gets used for which operations? -ffloat-store is supposed to force the constant discarding of extra digits, but we certainly don't want to apply that to all of numpy. So if we have to get log1p working ourselves it'll be an exercise in suppressing specific optimizations in all optimizing compilers that might compile numpy. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] ufuncs on funny strides; also isnan, isfinite, etc on a variety of dtypes

On 1 April 2010 14:30, M Trumpis mtrum...@berkeley.edu wrote: Hi all, disclaimer: pardon my vast ignorance on the subject of ufuncs, that will explain the naivety of the following questions This morning I was looking at this line of code, which was running quite slow for me and making me think data_has_nan = numpy.isnan(data_array).any() I knew that the array was of type uint8. This raised the question of how much slower would a function like isnan() run on integer typed data, and furthermore if it's integer typed data then why bother looking at the values at all? Yes, in fact there is no NaN for integer types, so this will never come out True. Actually I realized later that the main slow-down comes from the fact that my array was strided in fortran order (ascending strides). But from the point of view of a ufunc that is operating independently at each value, why would it need to respect striding? This is a subject that has come up a number of times. In an ideal world, numpy would be free to run ufuncs in a cache-optimal order. As it stands, though, numpy does some limited and possibly ill-advised reordering, and in fact the result of certain operations can depend on the order in which numpy carries out an operation, so it's kind of a pain to fix this. But we'd like to. The conclusion seems to be to do it in two steps: first, declare by fiat that all operations reading from and writing to the same memory return the same result as if the result were saved to a new array and then copied into the result. This changes the result of certain currently ill-defined operations, and requires real copies in some cases, but would remove undefined behaviour. The second step, once we have the freedom to reorder operations, would be to write some ufunc reordering code that flattens arrays as much as possible and loops over cache-coherent strides first if possible. Somebody needs to write the code, though. Anne And a last mini question, it doesn't appear that any() is doing short circuit evaluation. It runs in appx the same time whether an array is sparsely nonzero, fully zero, or fully nonzero. Kind regards, Mike ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Bug in logaddexp2.reduce

On 1 April 2010 13:38, Charles R Harris charlesr.har...@gmail.com wrote: On Thu, Apr 1, 2010 at 8:37 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Thu, Apr 1, 2010 at 12:46 AM, Anne Archibald peridot.face...@gmail.com wrote: On 1 April 2010 02:24, Charles R Harris charlesr.har...@gmail.com wrote: On Thu, Apr 1, 2010 at 12:04 AM, Anne Archibald peridot.face...@gmail.com wrote: On 1 April 2010 01:59, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:46 PM, Anne Archibald peridot.face...@gmail.com wrote: On 1 April 2010 01:40, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:25 PM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:22 AM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:17 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 6:08 PM, josef.p...@gmail.com wrote: On Wed, Mar 31, 2010 at 7:37 PM, Warren Weckesser warren.weckes...@enthought.com wrote: T J wrote: On Wed, Mar 31, 2010 at 1:21 PM, Charles R Harris charlesr.har...@gmail.com wrote: Looks like roundoff error. So this is expected behavior? In [1]: np.logaddexp2(-1.5849625007211563, -53.584962500721154) Out[1]: -1.5849625007211561 In [2]: np.logaddexp2(-0.5849625007211563, -53.584962500721154) Out[2]: nan Is any able to reproduce this? I don't get 'nan' in either 1.4.0 or 2.0.0.dev8313 (32 bit Mac OSX). In an earlier email T J reported using 1.5.0.dev8106. np.logaddexp2(-0.5849625007211563, -53.584962500721154) nan np.logaddexp2(-1.5849625007211563, -53.584962500721154) -1.5849625007211561 np.version.version '1.4.0' WindowsXP 32 What compiler? Mingw? yes, mingw 3.4.5. , official binaries release 1.4.0 by David sse2 Pentium M Can you try the exp2/log2 functions with the problem data and see if something goes wrong? Works fine for me. If it helps clarify things, the difference between the two problem input values is exactly 53 (and that's what logaddexp2 does an exp2 of); so I can provide a simpler example: In [23]: np.logaddexp2(0, -53) Out[23]: nan Of course, for me it fails in both orders. Ah, that's progress then ;) The effective number of bits in a double is 53 (52 + implicit bit). That shouldn't cause problems but it sure looks suspicious. Indeed, that's what led me to the totally wrong suspicion that denormals have something to do with the problem. More data points: In [38]: np.logaddexp2(63.999, 0) Out[38]: nan In [39]: np.logaddexp2(64, 0) Out[39]: 64.0 In [42]: np.logaddexp2(52.999, 0) Out[42]: 52.9990002 In [43]: np.logaddexp2(53, 0) Out[43]: nan It looks to me like perhaps the NaNs are appearing when the smaller term affects only the extra bits provided by the FPU's internal larger-than-double representation. Some such issue would explain why the problem seems to be hardware- and compiler-dependent. Hmm, that 63.999 is kinda strange. Here is a bit of code that might confuse a compiler working with different size mantissas: @type@ npy_log2...@c@(@type@ x) { @type@ u = 1 + x; if (u == 1) { return LOG2E*x; } else { return npy_l...@c@(u) * x / (u - 1); } } It might be that u != 1 does not imply u-1 != 0. That does indeed look highly suspicious. I'm not entirely sure how to work around it. GSL uses a volatile declaration: http://www.google.ca/codesearch/p?hl=en#p9nGS4eQGUI/gnu/gsl/gsl-1.8.tar.gz%7C8VCQSLJ5jR8/gsl-1.8/sys/log1p.cq=log1p On the other hand boost declares itself defeated by optimizing compilers and uses a Taylor series: http://www.google.ca/codesearch/p?hl=en#sdP2GRSfgKo/dcplusplus/trunk/boost/boost/math/special_functions/log1p.hppq=log1psa=Ncd=7ct=rc While R makes no mention of the corrected formula or optimizing compilers but takes the same approach, only with Chebyshev series: http://www.google.ca/codesearch/p?hl=en#gBBSWbwZmuk/src/base/R-2/R-2.3.1.tar.gz%7CVuh8XhRbUi8/R-2.3.1/src/nmath/log1p.cq=log1p Since, at least on my machine, ordinary log1p appears to work fine, is there any reason not to have log2_1p call it and scale the result? Or does the compiler make a hash of our log1p too? Calling log1p and scaling looks like the right thing to do here. And our log1p needs improvement. Tinkering a bit, I think we should implement the auxiliary function f(p) = log((1+p)/(1 - p)), which is antisymmetric and has the expansion 2p*(1 + p^2/3 + p^4/5 + ...). The series in the parens is increasing, so it is easy to terminate. Note that for p = +/- 1

### Re: [Numpy-discussion] Bug in logaddexp2.reduce

On 31 March 2010 16:21, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:38 AM, T J tjhn...@gmail.com wrote: On Wed, Mar 31, 2010 at 10:30 AM, T J tjhn...@gmail.com wrote: Hi, I'm getting some strange behavior with logaddexp2.reduce: from itertools import permutations import numpy as np x = np.array([-53.584962500721154, -1.5849625007211563, -0.5849625007211563]) for p in permutations([0,1,2]): print p, np.logaddexp2.reduce(x[list(p)]) Essentially, the result depends on the order of the array...and we get nans in the bad orders. Likely, this also affects logaddexp. Sorry, forgot version information: $ python -c import numpy;print numpy.__version__ 1.5.0.dev8106 __ Looks like roundoff error. No. The whole point of logaddexp and logaddexp2 is that they function correctly - in particular, avoid unnecessary underflow and overflow - in this sort of situation. This is a genuine problem. I see it using the stock numpy in Ubuntu karmic: In [3]: np.logaddexp2(-0.5849625007211563, -53.584962500721154) Out[3]: nan In [4]: np.logaddexp2(-53.584962500721154, -0.5849625007211563) Out[4]: nan In [5]: np.log2(2**(-53.584962500721154) + 2**(-0.5849625007211563)) Out[5]: -0.58496250072115608 In [6]: numpy.__version__ Out[6]: '1.3.0' In [8]: !uname -a Linux octopus 2.6.31-20-generic #58-Ubuntu SMP Fri Mar 12 05:23:09 UTC 2010 i686 GNU/Linux (the machine is a 32-bit Pentium M laptop.) I do not see a similar problem with logaddexp. Anne Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Bug in logaddexp2.reduce

On 1 April 2010 01:11, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:03 PM, Anne Archibald peridot.face...@gmail.com wrote: On 31 March 2010 16:21, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:38 AM, T J tjhn...@gmail.com wrote: On Wed, Mar 31, 2010 at 10:30 AM, T J tjhn...@gmail.com wrote: Hi, I'm getting some strange behavior with logaddexp2.reduce: from itertools import permutations import numpy as np x = np.array([-53.584962500721154, -1.5849625007211563, -0.5849625007211563]) for p in permutations([0,1,2]): print p, np.logaddexp2.reduce(x[list(p)]) Essentially, the result depends on the order of the array...and we get nans in the bad orders. Likely, this also affects logaddexp. Sorry, forgot version information: $ python -c import numpy;print numpy.__version__ 1.5.0.dev8106 __ Looks like roundoff error. No. The whole point of logaddexp and logaddexp2 is that they function correctly - in particular, avoid unnecessary underflow and overflow - in this sort of situation. This is a genuine problem. Well, it did function correctly for me ;) Just some roundoff error. I see it using the stock numpy in Ubuntu karmic: In [3]: np.logaddexp2(-0.5849625007211563, -53.584962500721154) Out[3]: nan In [4]: np.logaddexp2(-53.584962500721154, -0.5849625007211563) Out[4]: nan In [5]: np.log2(2**(-53.584962500721154) + 2**(-0.5849625007211563)) Out[5]: -0.58496250072115608 In [6]: numpy.__version__ Out[6]: '1.3.0' In [8]: !uname -a Linux octopus 2.6.31-20-generic #58-Ubuntu SMP Fri Mar 12 05:23:09 UTC 2010 i686 GNU/Linux (the machine is a 32-bit Pentium M laptop.) I think this is related to 32 bits somehow. The code for the logaddexp and logaddexp2 functions is almost identical. Maybe the compiler? Possibly, or possibly I just didn't find values that trigger the bug for logaddexp. I wonder if the problem is the handling of denormalized floats? Anne Chuck ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Bug in logaddexp2.reduce

On 1 April 2010 01:40, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 11:25 PM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:22 AM, josef.p...@gmail.com wrote: On Thu, Apr 1, 2010 at 1:17 AM, Charles R Harris charlesr.har...@gmail.com wrote: On Wed, Mar 31, 2010 at 6:08 PM, josef.p...@gmail.com wrote: On Wed, Mar 31, 2010 at 7:37 PM, Warren Weckesser warren.weckes...@enthought.com wrote: T J wrote: On Wed, Mar 31, 2010 at 1:21 PM, Charles R Harris charlesr.har...@gmail.com wrote: Looks like roundoff error. So this is expected behavior? In [1]: np.logaddexp2(-1.5849625007211563, -53.584962500721154) Out[1]: -1.5849625007211561 In [2]: np.logaddexp2(-0.5849625007211563, -53.584962500721154) Out[2]: nan Is any able to reproduce this? I don't get 'nan' in either 1.4.0 or 2.0.0.dev8313 (32 bit Mac OSX). In an earlier email T J reported using 1.5.0.dev8106. np.logaddexp2(-0.5849625007211563, -53.584962500721154) nan np.logaddexp2(-1.5849625007211563, -53.584962500721154) -1.5849625007211561 np.version.version '1.4.0' WindowsXP 32 What compiler? Mingw? yes, mingw 3.4.5. , official binaries release 1.4.0 by David sse2 Pentium M Can you try the exp2/log2 functions with the problem data and see if something goes wrong? Works fine for me. If it helps clarify things, the difference between the two problem input values is exactly 53 (and that's what logaddexp2 does an exp2 of); so I can provide a simpler example: In [23]: np.logaddexp2(0, -53) Out[23]: nan Of course, for me it fails in both orders. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Dealing with roundoff error

On 27 March 2010 19:38, Mike Sarahan msara...@gmail.com wrote: Hi all, I have run into some roundoff problems trying to line up some experimental spectra. The x coordinates are given in intervals of 0.1 units. I read the data in from a text file using np.loadtxt(). I think Robert's post here explains why the problem exists: http://mail.scipy.org/pipermail/numpy-discussion/2007-June/028133.html However, even linspace shows roundoff error: a=np.linspace(0.0,10.0,endpoint=False) b=np.linspace(0.1,10.1,endpoint=False) np.sum(a[1:]==b[:-1]) # Gives me 72, no 100 What is the best way to deal with it? Multiply the intervals by 10, then convert them to ints? It is almost never a good idea to compare floats for equality. (Exceptions include mostly situations where the float is not being operated on at all.) If your problem is that your spectra are really sampled at the same points but the floats coming out are slightly different, it's probably enough to test for abs(x-y)abs(x+y)*1e-13 ; as long as you don't to too many operations on the doubles, this should be enough elbow room to cover roundoff error. If your spectra are sampled at genuinely different points, you may want to look into some sort of interpolation to resample them to the same points. Anne Thanks, Mike ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Interpolation question

On 27 March 2010 20:24, Andrea Gavana andrea.gav...@gmail.com wrote: Hi All, I have an interpolation problem and I am having some difficulties in tackling it. I hope I can explain myself clearly enough. Basically, I have a whole bunch of 3D fluid flow simulations (close to 1000), and they are a result of different combinations of parameters. I was planning to use the Radial Basis Functions in scipy, but for the moment let's assume, to simplify things, that I am dealing only with one parameter (x). In 1000 simulations, this parameter x has 1000 values, obviously. The problem is, the outcome of every single simulation is a vector of oil production over time (let's say 40 values per simulation, one per year), and I would like to be able to interpolate my x parameter (1000 values) against all the simulations (1000x40) and get an approximating function that, given another x parameter (of size 1x1) will give me back an interpolated production profile (of size 1x40). If I understand your problem correctly, you have a function taking one value as input (or one 3D vector) and returning a vector of length 40. You want to know whether there are tools in scipy to support this. I'll say first that it's not strictly necessary for there to be: you could always just build 40 different interpolators, one for each component of the output. After all, there's no interaction in the calculations between the output coordinates. This is of course awkward, in that you'd like to just call F(x) and get back a vector of length 40, but that can be remedied by writing a short wrapper function that simply calls all 40 interpolators. A problem that may be more serious is that the python loop over the 40 interpolators can be slow, while a C implementation would give vector-valued results rapidly. To make this work, you're going to have to find a compiled-code interpolator that returns vector values. This is not in principle complicated, since they just need to run the same interpolator code on 40 sets of coefficients. But I don't think many of scipy's interpolators do this. The only ones I'm sure are able to do this are the polynomial interpolators I wrote, which work only on univariate inputs (and provide no kind of smoothing). If you're going to use these I recommend using scipy's spline functions to construct smoothing splines, which you'd then convert to a piecewise cubic. But I'd say, give the 40 interpolators a try. If they're slow, try interpolating many points at once: rather than giving just one x value, call the interpolators with a thousand (or however many you need) at a time. This will speed up scipy's interpolators, and it will make the overhead of a forty-element loop negligible. Anne Something along these lines: import numpy as np from scipy.interpolate import Rbf # x.shape = (1000, 1) # y.shape = (1000, 40) rbf = Rbf(x, y) # New result with xi.shape = (1, 1) -- fi.shape = (1, 40) fi = rbf(xi) Does anyone have a suggestion on how I could implement this? Sorry if it sounds confused... Please feel free to correct any wrong assumptions I have made, or to propose other approaches if you think RBFs are not suitable for this kind of problems. Thank you in advance for your suggestions. Andrea. Imagination Is The Only Weapon In The War Against Reality. http://xoomer.alice.it/infinity77/ == Never *EVER* use RemovalGroup for your house removal. You'll regret it forever. http://thedoomedcity.blogspot.com/2010/03/removal-group-nightmare.html == ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] [OT] Starving CPUs article featured in IEEE's ComputingNow portal

On 22 March 2010 14:42, Pauli Virtanen p...@iki.fi wrote: la, 2010-03-20 kello 17:36 -0400, Anne Archibald kirjoitti: I was in on that discussion. My recollection of the conclusion was that on the one hand they're useful, carefully applied, while on the other hand they're very difficult to reliably detect (since you don't want to forbid operations on non-overlapping slices of the same array). I think one alternative brought up was copy if unsure whether the slices overlap which would make A[whatever] = A[whatever2] be always identical in functionality to A[whatever] = A[whatever2].copy() which is how things should work. This would permit optimizing simple cases (at least 1D), and avoids running into NP-completeness (for numpy, the exponential growth is however limited by NPY_MAXDIMS which is 64, IIRC). It can produce surprise copies, but I could certainly live with this. Or maybe a slight modification: always produces values equivalent to using a copy, to allow handling the common A[:-1]=A[1:] without a copy. Of course, we'd have to wait for someone to implement it... This would be a change in semantics, but in a very obscure corner that hopefully nobody relies on. It would certainly be nice to replace unspecified behaviour by specified behaviour if it can be done with minimal cost. And I think it could be, with some careful programming. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] [OT] Starving CPUs article featured in IEEE's ComputingNow portal

On 20 March 2010 06:32, Francesc Alted fal...@pytables.org wrote: A Friday 19 March 2010 18:13:33 Anne Archibald escrigué: [clip] What I didn't go into in detail in the article was that there's a trade-off of processing versus memory access available: we could reduce the memory load by a factor of eight by doing interpolation on the fly instead of all at once in a giant FFT. But that would cost cache space and flops, and we're not memory-dominated. One thing I didn't try, and should: running four of these jobs at once on a four-core machine. If I correctly understand the architecture, that won't affect the cache issues, but it will effectively quadruple the memory bandwidth needed, without increasing the memory bandwidth available. (Which, honestly, makes me wonder what the point is of building multicore machines.) Maybe I should look into that interpolation stuff. Please do. Although you may be increasing the data rate by 4x, your program is already very efficient in how it handles data, so chances are that you still get a good speed-up. I'd glad to hear you back on your experience. The thing is, it reduces the data rate from memory, but at the cost of additional FFTs (to implement convolutions). If my program is already spending all its time doing FFTs, and the loads from memory are happening while the CPU does FFTs, then there's no win in runtime from reducing the memory load, and there's a runtime cost of doing those convolutions - not just more flops but also more cache pressure (to store the interpolated array and the convolution kernels). One could go a further step and do interpolation directly, without convolution, but that adds really a lot of flops, which translates directly to runtime. On the other hand, if it doesn't completely blow out the cache, we do have non-interpolated FFTs already on disk (with red noise adjustments already applied), so we might save on the relatively minor cost of the giant FFT. I'll have to do some time trials. Anne -- Francesc Alted ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] [OT] Starving CPUs article featured in IEEE's ComputingNow portal

On 20 March 2010 14:56, Dag Sverre Seljebotn da...@student.matnat.uio.no wrote: Pauli Virtanen wrote: Anne Archibald wrote: I'm not knocking numpy; it does (almost) the best it can. (I'm not sure of the optimality of the order in which ufuncs are executed; I think some optimizations there are possible.) Ufuncs and reductions are not performed in a cache-optimal fashion, IIRC dimensions are always traversed in order from left to right. Large speedups are possible in some cases, but in a quick try I didn't manage to come up with an algorithm that would always improve the speed (there was a thread about this last year or so, and there's a ticket). Things varied between computers, so this probably depends a lot on the actual cache arrangement. But perhaps numexpr has such heuristics, and we could steal them? At least in MultiIter (and I always assumed ufuncs too, but perhaps not) there's functionality to remove the largest dimension so that it can be put innermost in a loop. In many situations, removing the dimension with the smallest stride from the iterator would probably work much better. It's all about balancing iterator overhead and memory overhead. Something simple like select the dimension with length 200 which has smallest stride, or the dimension with largest length if none are above 200 would perhaps work well? There's more to it than that: there's no point (I think) going with the smallest stride if that stride is more than a cache line (usually 64 bytes). There are further optimizations - often stride[i]*shape[i]==shape[j] for some (i,j), and in that case these two dimensions can be treated as one with stride stride[i] and length shape[i]*shape[j]. Applying this would unroll all C- or Fortran-contiguous arrays to a single non-nested loop. Of course, reordering ufuncs is potentially hazardous in terms of breaking user code: the effect of A[1:]*=A[:-1] depends on whether you run through A in index order or reverse index order (which might be the order it's stored in memory, potentially faster from a cache point of view). Anne Dag Sverre ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] [OT] Starving CPUs article featured in IEEE's ComputingNow portal

On 20 March 2010 16:18, Sebastian Haase seb.ha...@gmail.com wrote: On Sat, Mar 20, 2010 at 8:22 PM, Anne Archibald peridot.face...@gmail.com wrote: On 20 March 2010 14:56, Dag Sverre Seljebotn da...@student.matnat.uio.no wrote: Pauli Virtanen wrote: Anne Archibald wrote: I'm not knocking numpy; it does (almost) the best it can. (I'm not sure of the optimality of the order in which ufuncs are executed; I think some optimizations there are possible.) Ufuncs and reductions are not performed in a cache-optimal fashion, IIRC dimensions are always traversed in order from left to right. Large speedups are possible in some cases, but in a quick try I didn't manage to come up with an algorithm that would always improve the speed (there was a thread about this last year or so, and there's a ticket). Things varied between computers, so this probably depends a lot on the actual cache arrangement. But perhaps numexpr has such heuristics, and we could steal them? At least in MultiIter (and I always assumed ufuncs too, but perhaps not) there's functionality to remove the largest dimension so that it can be put innermost in a loop. In many situations, removing the dimension with the smallest stride from the iterator would probably work much better. It's all about balancing iterator overhead and memory overhead. Something simple like select the dimension with length 200 which has smallest stride, or the dimension with largest length if none are above 200 would perhaps work well? There's more to it than that: there's no point (I think) going with the smallest stride if that stride is more than a cache line (usually 64 bytes). There are further optimizations - often stride[i]*shape[i]==shape[j] for some (i,j), and in that case these two dimensions can be treated as one with stride stride[i] and length shape[i]*shape[j]. Applying this would unroll all C- or Fortran-contiguous arrays to a single non-nested loop. Of course, reordering ufuncs is potentially hazardous in terms of breaking user code: the effect of A[1:]*=A[:-1] depends on whether you run through A in index order or reverse index order (which might be the order it's stored in memory, potentially faster from a cache point of view). I think Travis O suggested at some point to make this kind of overlapping inplace operations invalid. (maybe it was someone else ...) The idea was to protect new (uninitiated) users from unexpected results. Anyway, maybe one could use (something like) np.seterror to determine how to handle these cases... I was in on that discussion. My recollection of the conclusion was that on the one hand they're useful, carefully applied, while on the other hand they're very difficult to reliably detect (since you don't want to forbid operations on non-overlapping slices of the same array). Anne -Sebastian ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] [OT] Starving CPUs article featured in IEEE's ComputingNow portal

On 18 March 2010 13:53, Francesc Alted fal...@pytables.org wrote: A Thursday 18 March 2010 16:26:09 Anne Archibald escrigué: Speak for your own CPUs :). But seriously, congratulations on the wide publication of the article; it's an important issue we often don't think enough about. I'm just a little snarky because this exact issue came up for us recently - a visiting astro speaker put it as flops are free - and so I did some tests and found that even without optimizing for memory access, our tasks are already CPU-bound: http://lighthouseinthesky.blogspot.com/2010/03/flops.html Well, I thought that my introduction was enough to convince anybody about the problem, but forgot that you, the scientists, always try to demonstrate things experimentally :-/ Snrk. Well, technically, that is our job description... Seriously, your example is a clear example of what I'm recommending in the article, i.e. always try to use libraries that are already leverage the blocking technique (that is, taking advantage of both temporal and spatial locality). Don't know about FFTW (never used it, sorry), but after having a look at its home page, I'm pretty convinced that its authors are very conscious about these techniques. Being said this, it seems that, in addition, you are applying the blocking technique yourself also: get the data in bunches (256 floating point elements, which fits perfectly well on modern L1 caches), apply your computation (in this case, FFTW) and put the result back in memory. A perfect example of what I wanted to show to the readers so, congratulations! you made it without the need to read my article (so perhaps the article was not so necessary after all :-) What I didn't go into in detail in the article was that there's a trade-off of processing versus memory access available: we could reduce the memory load by a factor of eight by doing interpolation on the fly instead of all at once in a giant FFT. But that would cost cache space and flops, and we're not memory-dominated. One thing I didn't try, and should: running four of these jobs at once on a four-core machine. If I correctly understand the architecture, that won't affect the cache issues, but it will effectively quadruple the memory bandwidth needed, without increasing the memory bandwidth available. (Which, honestly, makes me wonder what the point is of building multicore machines.) Maybe I should look into that interpolation stuff. Heh. Indeed numexpr is a good tool for this sort of thing; it's an unfortunate fact that simple use of numpy tends to do operations in the pessimal order... Well, to honor the truth, NumPy does not have control in the order of the operations in expressions and how temporaries are managed: it is Python who decides that. NumPy only can do what Python wants it to do, and do it as good as possible. And NumPy plays its role reasonably well here, but of course, this is not enough for providing performance. In fact, this problem probably affects to all interpreted languages out there, unless they implement a JIT compiler optimised for evaluating expressions --and this is basically what numexpr is. I'm not knocking numpy; it does (almost) the best it can. (I'm not sure of the optimality of the order in which ufuncs are executed; I think some optimizations there are possible.) But a language designed from scratch for vector calculations could certainly compile expressions into a form that would save a lot of memory accesses, particularly if an optimizer combined many lines of code. I've actually thought about whether such a thing could be done in python; I think the way to do it would be to build expression objects from variable objects, then have a single apply function that fed values in to all the variables. The same framework would support automatic differentiation and other cool things, but I'm not sure it would be useful enough to be worth the implementation complexity. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] [OT] Starving CPUs article featured in IEEE's ComputingNow portal

On 18 March 2010 09:57, Francesc Alted fal...@pytables.org wrote: Hi, Konrad Hinsen has just told me that my article Why Modern CPUs Are Starving and What Can Be Done About It, which has just released on the March/April issue of Computing in Science and Engineering, also made into this month's free-access selection on IEEE's ComputingNow portal: Speak for your own CPUs :). But seriously, congratulations on the wide publication of the article; it's an important issue we often don't think enough about. I'm just a little snarky because this exact issue came up for us recently - a visiting astro speaker put it as flops are free - and so I did some tests and found that even without optimizing for memory access, our tasks are already CPU-bound: http://lighthouseinthesky.blogspot.com/2010/03/flops.html In terms of specifics, I was a little surprised you didn't mention FFTW among your software tools that optimize memory access. FFTW's planning scheme seems ideal for ensuring memory locality, as much as possible, during large FFTs. (And in fact I also found that for really large FFTs, reducing padding - memory size - at the cost of a non-power-of-two size was also worth it.) http://www.computer.org/portal/web/computingnow http://www.computer.org/portal/web/computingnow/0310/whatsnew/cise On it, I discuss one of my favourite subjects, memory access, and why it is important for nowadays computing. There are also recommendations for people wanting to squeeze the maximum of performance out of their computers. And, although I tried to be as language-agnostic as I could, there can be seen some Python references here and there :-). Heh. Indeed numexpr is a good tool for this sort of thing; it's an unfortunate fact that simple use of numpy tends to do operations in the pessimal order... Ane Well, sorry about this semi-OT but I could not resist :-) -- Francesc Alted ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] matrix division without a loop

On 12 March 2010 13:54, gerardo.berbeglia gberbeg...@gmail.com wrote: Hello, I want to divide an n x n (2-dimension) numpy array matrix A by a n (1-dimension) array d as follows: Look up broadcasting in the numpy docs. The short version is this: operations like division act elementwise on arrays of the same shape. If two arrays' shapes differ, but only in that one has a 1 where the other has an n, the one whose dimension is 1 will be effectively repeated n times along that axis so they match. Finally, if two arrays have shapes of different lengths, the shorter is extended by prepending as many dimensions of size one as necessary, then the previous rules are applied. Put this together with the ability to add extra dimensions of length 1 to an array by indexing with np.newaxis, and your problem becomes simple: In [13]: A = np.array([[1,2],[3,4]]); B = np.array([1.,2.]) In [14]: A.shape Out[14]: (2, 2) In [15]: B.shape Out[15]: (2,) In [16]: A/B Out[16]: array([[ 1., 1.], [ 3., 2.]]) In [17]: A/B[:,np.newaxis] Out[17]: array([[ 1. , 2. ], [ 1.5, 2. ]]) If you're using matrices, don't. Stick to arrays. Anne Take n = 2. Let A= 2 3 1 10 and let d = [ 3 2 ] Then i would like to have A/d = 2/3 3/3 1/2 10/2 This is to avoid loops to improve the speed. Thank you in advance! -- View this message in context: http://old.nabble.com/matrix-division-without-a-loop-tp27881350p27881350.html Sent from the Numpy-discussion mailing list archive at Nabble.com. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] arange including stop value?

On 11 March 2010 19:30, Tom K. t...@kraussfamily.org wrote: davefallest wrote: ... In [3]: np.arange(1.01, 1.1, 0.01) Out[3]: array([ 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.1 ]) Why does the ... np.arange command end up including my stop value? Don't use arange for floating-point values. Use linspace instead. Anne From the help for arange: For floating point arguments, the length of the result is ``ceil((stop - start)/step)``. Because of floating point overflow, this rule may result in the last element of `out` being greater than `stop`. -- View this message in context: http://old.nabble.com/arange-including-stop-value--tp27866607p27872069.html Sent from the Numpy-discussion mailing list archive at Nabble.com. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] fast duplicate of array

2010/1/23 Alan G Isaac ais...@american.edu: Suppose x and y are conformable 2d arrays. I now want x to become a duplicate of y. I could create a new array: x = y.copy() or I could assign the values of y to x: x[:,:] = y As expected the latter is faster (no array creation). Are there better ways? If both arrays are C contiguous, or more generally contiguous blocks of memory with the same strided structure, you might get faster copying by flattening them first, so that it can go in a single memcpy(). For really large arrays that use complete pages, some low-level hackery involving memmap() might be able to make a shared copy-on-write copy at almost no cost until you start modifying one array or the other. But both of these tricks are intended for the regime where copying the data is the expensive part, not fabricating the array object; for that, I'm not sure you can accelerate things much. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] fast duplicate of array

2010/1/23 Alan G Isaac ais...@american.edu: On 1/23/2010 5:01 PM, Anne Archibald wrote: If both arrays are C contiguous, or more generally contiguous blocks of memory with the same strided structure, you might get faster copying by flattening them first, so that it can go in a single memcpy(). I may misuderstand this. Did you just mean x.flat = y.flat ? No, .flat constructs an iterator that traverses the object as if it were flat. I had in mind accessing the underlying data through views that were flat: In [3]: x = np.random.random((1000,1000)) In [4]: y = np.random.random((1000,1000)) In [5]: xf = x.view() In [6]: xf.shape = (-1,) In [7]: yf = y.view() In [8]: yf.shape = (-1,) In [9]: yf[:] = xf[:] This may still use a loop instead of a memcpy(), in which case you'd want to look for an explicit memcpy()-based implementation, but when manipulating multidimensional arrays you have (in principle, anyway) nested loops which may not be executed in the cache-optimal order. Ideally numpy would automatically notice when operations can be done on flattened versions of arrays and get rid of some of the looping and indexing, but I wouldn't count on it. At one point I remember finding that the loops were reordered not for cache coherence but to make the inner loop over the biggest dimension (to minimize looping overhead). Anne If so, I find that to be *much* slower. Thanks, Alan x = np.random.random((1000,1000)) y = x.copy() t0 = time.clock() for t in range(1000): x = y.copy() print(time.clock() - t0) t0 = time.clock() for t in range(1000): x[:,:] = y print(time.clock() - t0) t0 = time.clock() for t in range(1000): x.flat = y.flat print(time.clock() - t0) ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Efficient orthogonalisation with scipy/numpy

2010/1/19 Gael Varoquaux gael.varoqu...@normalesup.org: For the google-completness of this thread, to get a speed gain, one needs to use the 'econ=True' flag to qr. Be warned that in some installations (in particular some using ATLAS), supplying econ=True can cause a segfault under certain conditions (I think only when the arrays are misaligned, e.g. coming from unpickling, and even then only with certain data); the bugfix for this may or may not be complete, and involves copying misaligned arrays. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] Efficient orthogonalisation with scipy/numpy

2010/1/19 Charles R Harris charlesr.har...@gmail.com: Note that if you apply the QR algorithm to a Vandermonde matrix with the columns properly ordered you can get a collection of graded orthogonal polynomials over a given set of points. Or, if you want the polynomials in some other representation - by values, or in terms of some basis of orthogonal polynomials - you can construct an appropriate Vandermonde-style matrix and use QR. (When I tried this, switching from the power basis to the Chebyshev basis let me go from tens to hundreds of polynomials, and now Chebyshev polynomials are first-class objects.) Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] LA improvements (was: dot function or dot notation, matrices, arrays?)

2009/12/23 David Goldsmith d.l.goldsm...@gmail.com: Starting a new thread for this. On Tue, Dec 22, 2009 at 7:13 PM, Anne Archibald peridot.face...@gmail.com wrote: I think we have one major lacuna: vectorized linear algebra. If I have to solve a whole whack of four-dimensional linear systems, right now I need to either write a python loop and use linear algebra on them one by one, or implement my own linear algebra. It's a frustrating lacuna, because all the machinery is there: generalized ufuncs and LAPACK wrappers. Somebody just needs to glue them together. I've even tried making a start on it, but numpy's ufunc machinery and generic type system is just too much of a pain for me to make any progress as is. Please be more specific: what (which aspects) have been too much of a pain? (I ask out of ignorance, not out of challenging your opinion/experience.) It's been a little while since I took a really close look at it, but I'll try to describe the problems I had. Chiefly I had problems with documentation - the only way I could figure out how to build additional gufuncs was monkey-see-monkey-do, just copying an existing one in an existing file and hoping the build system figured it out. It was also not at all clear how to, say, link to LAPACK, let alone decide based on input types which arguments to promote and how to call out to LAPACK. I'm not saying this is impossible, just that it was enough frustrating no-progress to defeat my initial hey, I could do that impulse. I think if someone wanted to start building a low-level Again, please be more specific: what do you mean by this? (I know generally what is meant by low level, but I'd like you to spell out a little more fully what you mean by this in this context.) Sure. Let me first say that all this is kind of beside the point - the hard part is not designing an API, so it's a bit silly to dream up an API without implementing anything. I had pictured two interfaces to the vectorized linear algebra code. The first would simply provide more-or-less direct access to vectorized versions of the linear algebra functions we have now, with no dimension inference. Thus inv, pinv, svd, lu factor, lu solve, et cetera - but not dot. Dot would have to be split up into vector-vector, vector-matrix, matrix-vector, and matrix-matrix products, since one can no longer use the dimensionality of the inputs to figure out what is wanted. The key idea would be that the linear algebra dimensions would always be the last one(s); this is fairly easy to arrange with rollaxis when it isn't already true, would tend to reduce copying on input to LAPACK, and is what the gufunc API wants. This is mostly what I meant by low-level. (A second generation would do things like combine many vector-vector products into a single LAPACK matrix-vector product.) The higher-level API I was imagining - remember, vaporware here - had a Matrix and a Vector class, each holding an arbitrarily-dimensioned array of the relevant object. The point of this is to avoid having to constantly specify whether you want a matrix-vector or matrix-matrix product; it also tidily avoids the always-two-dimensional nuisance of the current matrix API. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] LA improvements (was: dot function or dot notation, matrices, arrays?)

2009/12/23 David Warde-Farley d...@cs.toronto.edu: On 23-Dec-09, at 10:34 AM, Anne Archibald wrote: The key idea would be that the linear algebra dimensions would always be the last one(s); this is fairly easy to arrange with rollaxis when it isn't already true, would tend to reduce copying on input to LAPACK, and is what the gufunc API wants. Would it actually reduce copying if you were using default C-ordered arrays? Maybe I'm mistaken but I thought one almost always had to copy in order to translate C to Fortran order except for a few functions that can take row-ordered stuff. That's a good point. But even if you need to do a transpose, it'll be faster to transpose data in a contiguous block than data scattered all over memory. Maybe more to the point, broadcasting adds axes to the beginning, so that (say) two-dimensional arrays can act as matrix scalars. Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

### Re: [Numpy-discussion] dot function or dot notation, matrices, arrays?

2009/12/21 David Goldsmith d.l.goldsm...@gmail.com: On Mon, Dec 21, 2009 at 9:57 AM, Christopher Barker chris.bar...@noaa.gov wrote: Dag Sverre Seljebotn wrote: I recently got motivated to get better linear algebra for Python; wonderful! To me that seems like the ideal way to split up code -- let NumPy/SciPy deal with the array-oriented world and Sage the closer-to-mathematics notation. well, maybe -- but there is a lot of call for pure-computational linear algebra. I do hope you'll consider building the computational portion of it in a way that might be included in numpy or scipy by itself in the future. My personal opinion is that the LA status quo is acceptably good: there's maybe a bit of an adjustment to make for newbies, but I don't see it as a very big one, and this list strikes me as very efficient at getting people over little bumps (e.g., someone emails in: how do you matrix-multiply two arrays? within minutes (:-)) Robert or Charles replies with np.dot: np.dot([[1,2],[3,4]],[[1,2],[3,4]]) = array([[7,10],[15,22]])). Certainly any significant changes to the base should need to run the gauntlet of an NEP process. I think we have one major lacuna: vectorized linear algebra. If I have to solve a whole whack of four-dimensional linear systems, right now I need to either write a python loop and use linear algebra on them one by one, or implement my own linear algebra. It's a frustrating lacuna, because all the machinery is there: generalized ufuncs and LAPACK wrappers. Somebody just needs to glue them together. I've even tried making a start on it, but numpy's ufunc machinery and generic type system is just too much of a pain for me to make any progress as is. I think if someone wanted to start building a low-level generalized ufunc library interface to LAPACK, that would be a big improvement in numpy/scipy's linear algebra. Pretty much everything else strikes me as a question of notation. (Not to trivialize it: good notation makes a tremendous difference.) Anne ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion