Discussion of Numerical Python

headers Ralf Gommers | 1 Jun 02:19

_______________________________________________
Numpy-discussion mailing list
Numpy-discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Hi,

I have a set of n points with real coordinates between 0 and 1, given  
by two numpy arrays x and y, with a value at each point represented by  
a third array z. I am trying to then rasterize the points onto a grid  
of size npix*npix. So I can start by converting x and y to integer  
pixel coordinates ix and iy. But my question is, is there an efficient  
way to add z[i] to the pixel given by (xi[i],yi[i])? Below is what I  
am doing at the moment, but the for loop becomes very inefficient for  
large n. I would imagine that there is a way to do this without using  
a loop?

---

import numpy as np

n = 10000000
x = np.random.random(n)
y = np.random.random(n)
z = np.random.random(n)

npix = 100
ix = np.array(x*float(npix),int)
iy = np.array(y*float(npix),int)

image = np.zeros((npix,npix))
for i in range(len(ix)):
     image[ix[i],iy[i]] = image[ix[i],iy[i]] + z[i]

---

Thanks for any advice,

Thomas

Charles R Harris wrote:
>
>
> On Sun, May 31, 2009 at 11:54 AM, rob steed <rjsteed <at> talk21.com
> <mailto:rjsteed <at> talk21.com>> wrote:
>
>
>     Hi,
>     After my previous email, I have opened a ticket #1117 (correlate
>     not order dependent)
>
>     I have found that the correlate function is defined in
>     multiarraymodule.c and
>     that inputs are being swapped using the following code
>
>        n1 = ap1->dimensions[0];
>        n2 = ap2->dimensions[0];
>        if (n1 < n2) {
>            ret = ap1;
>            ap1 = ap2;
>            ap2 = ret;
>            ret = NULL;
>            i = n1;
>            n1 = n2;
>            n2 = i;
>        }
>
>     I do not know the code well enough to see whether this could just
>     be removed (I don't know c either).
>     Maybe the algorithmn requires the inputs to be length ordered? I
>     will try to work it out.
>
>
> If the correlation algorithm doesn't use an fft and is done
> explicitly, then the maximum overlap for any shift is the length of
> the shortest input.  Swapping the arrays makes that logic easier to
> implement, but it isn't necessary.

But this logic is also wrong if the swapping is not taken into account -
as the OP mentioned, correlate(a, b) is not equal to correlate(b, a) in
the general case. The output is reversed in the second case compared to
the first case.

cheers,

David

_______________________________________________
Numpy-discussion mailing list
Numpy-discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Charles R Harris wrote:
>
>
> On Sun, May 31, 2009 at 7:18 PM, David Cournapeau
> <david <at> ar.media.kyoto-u.ac.jp <mailto:david <at> ar.media.kyoto-u.ac.jp>>
> wrote:
>
>     Charles R Harris wrote:
>     >
>     >
>     > On Sun, May 31, 2009 at 11:54 AM, rob steed <rjsteed <at> talk21.com
>     <mailto:rjsteed <at> talk21.com>
>     > <mailto:rjsteed <at> talk21.com <mailto:rjsteed <at> talk21.com>>> wrote:
>     >
>     >
>     >     Hi,
>     >     After my previous email, I have opened a ticket #1117 (correlate
>     >     not order dependent)
>     >
>     >     I have found that the correlate function is defined in
>     >     multiarraymodule.c and
>     >     that inputs are being swapped using the following code
>     >
>     >        n1 = ap1->dimensions[0];
>     >        n2 = ap2->dimensions[0];
>     >        if (n1 < n2) {
>     >            ret = ap1;
>     >            ap1 = ap2;
>     >            ap2 = ret;
>     >            ret = NULL;
>     >            i = n1;
>     >            n1 = n2;
>     >            n2 = i;
>     >        }
>     >
>     >     I do not know the code well enough to see whether this could
>     just
>     >     be removed (I don't know c either).
>     >     Maybe the algorithmn requires the inputs to be length ordered? I
>     >     will try to work it out.
>     >
>     >
>     > If the correlation algorithm doesn't use an fft and is done
>     > explicitly, then the maximum overlap for any shift is the length of
>     > the shortest input.  Swapping the arrays makes that logic easier to
>     > implement, but it isn't necessary.
>
>     But this logic is also wrong if the swapping is not taken into
>     account -
>     as the OP mentioned, correlate(a, b) is not equal to correlate(b,
>     a) in
>     the general case. The output is reversed in the second case
>     compared to
>     the first case.
>
>
> I didn't say it was *correctly* implemented ;)

:) So I gave it a shot

http://github.com/cournape/numpy/commits/fix_correlate

(It took me a while to realize that PyArray_ISFLEXIBLE returns false for
array object. Is this expected ? The documentation concerning copyswap
says that it is necessary for flexible arrays, but I think it is
necessary  for object arrays as well).

It still bothers me that correlate does not conjugate the second
argument for complex arrays...

cheers,

David

On Sun, May 31, 2009 at 8:39 PM, Thomas Robitaille
<thomas.robitaille <at> gmail.com> wrote:
> Hi,
>
> I have a set of n points with real coordinates between 0 and 1, given
> by two numpy arrays x and y, with a value at each point represented by
> a third array z. I am trying to then rasterize the points onto a grid
> of size npix*npix. So I can start by converting x and y to integer
> pixel coordinates ix and iy. But my question is, is there an efficient
> way to add z[i] to the pixel given by (xi[i],yi[i])? Below is what I
> am doing at the moment, but the for loop becomes very inefficient for
> large n. I would imagine that there is a way to do this without using
> a loop?
>
> ---
>
> import numpy as np
>
> n = 10000000
> x = np.random.random(n)
> y = np.random.random(n)
> z = np.random.random(n)
>
> npix = 100
> ix = np.array(x*float(npix),int)
> iy = np.array(y*float(npix),int)
>
> image = np.zeros((npix,npix))
> for i in range(len(ix)):
>     image[ix[i],iy[i]] = image[ix[i],iy[i]] + z[i]
>

There might be a better way, but histogram2d() seems like a good fit:

import numpy as np

n = 1000000
x = np.random.random(n)
y = np.random.random(n)
z = np.random.random(n)

npix = 100

bins = np.linspace(0, 1.0, npix + 1)
image = np.histogram2d(x, y, bins=bins, weights=z)[0]

--

-- 
Nathan Bell wnbell <at> gmail.com
http://www.wnbell.com/

_______________________________________________
Numpy-discussion mailing list
Numpy-discussion <at> scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion

Charles R Harris wrote:
>
>
> On Sun, May 31, 2009 at 9:08 PM, David Cournapeau
> <david <at> ar.media.kyoto-u.ac.jp <mailto:david <at> ar.media.kyoto-u.ac.jp>>
> wrote:
>
>     Charles R Harris wrote:
>     >
>     >
>     > On Sun, May 31, 2009 at 7:18 PM, David Cournapeau
>     > <david <at> ar.media.kyoto-u.ac.jp
>     <mailto:david <at> ar.media.kyoto-u.ac.jp>
>     <mailto:david <at> ar.media.kyoto-u.ac.jp
>     <mailto:david <at> ar.media.kyoto-u.ac.jp>>>
>     > wrote:
>     >
>     >     Charles R Harris wrote:
>     >     >
>     >     >
>     >     > On Sun, May 31, 2009 at 11:54 AM, rob steed
>     <rjsteed <at> talk21.com <mailto:rjsteed <at> talk21.com>
>     >     <mailto:rjsteed <at> talk21.com <mailto:rjsteed <at> talk21.com>>
>     >     > <mailto:rjsteed <at> talk21.com <mailto:rjsteed <at> talk21.com>
>     <mailto:rjsteed <at> talk21.com <mailto:rjsteed <at> talk21.com>>>> wrote:
>     >     >
>     >     >
>     >     >     Hi,
>     >     >     After my previous email, I have opened a ticket #1117
>     (correlate
>     >     >     not order dependent)
>     >     >
>     >     >     I have found that the correlate function is defined in
>     >     >     multiarraymodule.c and
>     >     >     that inputs are being swapped using the following code
>     >     >
>     >     >        n1 = ap1->dimensions[0];
>     >     >        n2 = ap2->dimensions[0];
>     >     >        if (n1 < n2) {
>     >     >            ret = ap1;
>     >     >            ap1 = ap2;
>     >     >            ap2 = ret;
>     >     >            ret = NULL;
>     >     >            i = n1;
>     >     >            n1 = n2;
>     >     >            n2 = i;
>     >     >        }
>     >     >
>     >     >     I do not know the code well enough to see whether this
>     could
>     >     just
>     >     >     be removed (I don't know c either).
>     >     >     Maybe the algorithmn requires the inputs to be length
>     ordered? I
>     >     >     will try to work it out.
>     >     >
>     >     >
>     >     > If the correlation algorithm doesn't use an fft and is done
>     >     > explicitly, then the maximum overlap for any shift is the
>     length of
>     >     > the shortest input.  Swapping the arrays makes that logic
>     easier to
>     >     > implement, but it isn't necessary.
>     >
>     >     But this logic is also wrong if the swapping is not taken into
>     >     account -
>     >     as the OP mentioned, correlate(a, b) is not equal to
>     correlate(b,
>     >     a) in
>     >     the general case. The output is reversed in the second case
>     >     compared to
>     >     the first case.
>     >
>     >
>     > I didn't say it was *correctly* implemented ;)
>
>     :) So I gave it a shot
>
>     http://github.com/cournape/numpy/commits/fix_correlate
>
>     (It took me a while to realize that PyArray_ISFLEXIBLE returns
>     false for
>     array object. Is this expected ? The documentation concerning copyswap
>     says that it is necessary for flexible arrays, but I think it is
>     necessary  for object arrays as well).
>
>
> Don't know. PyArray_ISFLEXIBLE looks like a macro... 
>
> #define PyArray_ISFLEXIBLE(obj) PyTypeNum_ISFLEXIBLE(PyArray_TYPE(obj))
>
> #define PyTypeNum_ISFLEXIBLE(type) (((type) >=NPY_STRING) &&  \
>                                     ((type) <=NPY_VOID))
>
> And the typecodes are '?bhilqpBHILQPfdgFDGSUVO'. So 'SUV' are flexible
> and O is not.

I re-read the copyswap documentation, and realized I did not read it
correctly. Now, I am not sure when to use copyswap vs memcpy (memcpy
should be much faster on basic types, as memcpy should be inlined
generally, whereas I doubt copyswap can).

> I'm not clear on how correlate should apply to any of 'SUV' but it
> might be worth having it work for objects.

It already does (I added a couple of unit tests in the branch, since
there were no test for correlate, and one is for Decimal object arrays).

>
>
>     It still bothers me that correlate does not conjugate the second
>     argument for complex arrays...
>
>
> It bothers me also...

I think we should just fix it to use conjugate - I will do this in the
branch, and I will integrate it in the trunk later unless someone stands
up vehemently against the change. I opened up a ticket to track this,
though,

cheers,

David

Hi,

In the process of adding single precision support to Numexpr, I'm 
experimenting a divergence between Numexpr and NumPy computations.  It all 
boils down to the fact that my implementation defined single precision 
functions completely.  As for one, consider my version of expm1f:

inline static float expm1f(float x)
{
    float u = expf(x);
    if (u == 1.0) {
        return x;
    } else if (u-1.0 == -1.0) {
        return -1;
    } else {
        return (u-1.0) * x/logf(u);
    }
}

while NumPy seems to declare expm1f as:

static float expm1f(float x)
{
    return (float) expm1((double)x);
}

This leads to different results on Windows when computing expm1(x) for large 
values of x (like 99.), where my approach returns a 'nan', while NumPy returns 
an 'inf'.  Curiously, on Linux both approaches returns 'inf'.

I suppose that the NumPy crew already experimented this divergence and finally 
used the cast approach for computing the single precision functions.  However, 
this is effectively preventing the use of optimized functions for single 
precision (i.e. double precision 'exp' and 'log' are used instead of single 
precision specific 'expf' and 'logf'), which could perform potentially better.  
So, I'm wondering if it would not be better to use a native implementation 
instead.  Thoughts?

Thanks,

--

-- 
Francesc Alted

On Mon, Jun 1, 2009 at 00:05, David Cournapeau
<david <at> ar.media.kyoto-u.ac.jp> wrote:

> I think we should just fix it to use conjugate - I will do this in the
> branch, and I will integrate it in the trunk later unless someone stands
> up vehemently against the change. I opened up a ticket to track this,
> though,

It breaks everyone's code that works around the current behavior.

--

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

Mon	Tue	Wed	Thu	Fri	Sat	Sun

1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30

239	February 2012
366	January 2012
397	December 2011
259	November 2011
475	October 2011
280	September 2011
413	August 2011
536	July 2011
833	June 2011
369	May 2011
340	April 2011
729	March 2011
328	February 2011
359	January 2011
259	December 2010
463	November 2010
623	October 2010
503	September 2010
488	August 2010
787	July 2010
576	June 2010
440	May 2010
502	April 2010
641	March 2010
848	February 2010
441	January 2010
655	December 2009
735	November 2009
497	October 2009
911	September 2009
608	August 2009
594	July 2009
842	June 2009
590	May 2009
665	April 2009
897	March 2009
681	February 2009
587	January 2009
504	December 2008
532	November 2008
557	October 2008
778	September 2008
758	August 2008
1019	July 2008
620	June 2008
1244	May 2008
1110	April 2008
620	March 2008
551	February 2008
672	January 2008
317	December 2007
398	November 2007
391	October 2007
293	September 2007
319	August 2007
351	July 2007
282	June 2007
504	May 2007
617	April 2007
720	March 2007
435	February 2007
594	January 2007
361	December 2006
500	November 2006
873	October 2006
530	September 2006
636	August 2006
803	July 2006
664	June 2006
436	May 2006
543	April 2006
584	March 2006
673	February 2006
439	January 2006
61	December 2005
118	November 2005
145	October 2005
65	September 2005
56	August 2005
120	July 2005
58	June 2005
78	May 2005
205	April 2005
295	March 2005
318	February 2005
185	January 2005
60	December 2004
55	November 2004
131	October 2004
75	September 2004
84	August 2004
150	July 2004
149	June 2004
72	May 2004
49	April 2004
70	March 2004
75	February 2004
110	January 2004
71	December 2003
83	November 2003
61	October 2003
107	September 2003
24	August 2003
48	July 2003
12	June 2003
76	May 2003
6	April 2003
42	March 2003
72	February 2003
110	January 2003
104	December 2002
23	November 2002
36	October 2002
59	September 2002
34	August 2002
36	July 2002
123	June 2002
34	May 2002
62	April 2002
126	March 2002
37	February 2002
72	January 2002
54	December 2001
81	November 2001
33	October 2001
56	September 2001
92	August 2001
52	July 2001
33	June 2001
50	May 2001
67	April 2001
42	March 2001
54	February 2001
56	January 2001
24	December 2000
32	November 2000
60	October 2000
44	September 2000
16	August 2000
16	July 2000
28	June 2000
37	May 2000
28	April 2000
48	March 2000
49	February 2000
8	January 2000

resetting set_string_function

Rasterizing points onto an array

Re: Problem with correlate

Re: Problem with correlate

Re: Problem with correlate

Re: Rasterizing points onto an array

Re: Problem with correlate

Re: Problem with correlate

Single precision equivalents of missing C99 functions

Re: Problem with correlate