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

[Anne Archibald]

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?

[Anne Archibald]

On Mon, Mar 13, 2017 at 12:57 PM Eric Wieser 
wrote:

> > `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?

[Anne Archibald]

On Thu, Mar 9, 2017, 11:27 Nico Schlömer  wrote:

> 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

[Anne Archibald]

On Mon, Jan 23, 2017 at 3:34 PM Robert Kern  wrote:

> 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

[Anne Archibald]

On Wed, Jan 18, 2017 at 4:13 PM Nadav Har'El  wrote:

> 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

[Anne Archibald]

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

[Anne Archibald]

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 Harris 
wrote:

> 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

[Anne Archibald]

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 Rietdijk 
wrote:

> 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

[Anne Archibald]

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 Harris 
wrote:

> 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

[Anne Archibald]

On Fri, Sep 11, 2015 at 3:20 PM Sebastian Berg 
wrote:

> 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

[Anne Archibald]

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)

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?)

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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!?

[Anne Archibald]

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?

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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)?

[Anne Archibald]

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)?

[Anne Archibald]

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)?

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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?

[Anne Archibald]

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?

[Anne Archibald]

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()

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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

[Anne Archibald]

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?)

[Anne Archibald]

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?)

[Anne Archibald]

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?

[Anne Archibald]

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