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. Regards Rob
On Sun, May 31, 2009 at 11:54 AM, rob steed <rjsteed@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. Chuck
Charles R Harris wrote:
On Sun, May 31, 2009 at 11:54 AM, rob steed <rjsteed@talk21.com <mailto:rjsteed@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
On Sun, May 31, 2009 at 7:18 PM, David Cournapeau < david@ar.media.kyoto-u.ac.jp> wrote:
Charles R Harris wrote:
On Sun, May 31, 2009 at 11:54 AM, rob steed <rjsteed@talk21.com <mailto:rjsteed@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 ;) Chuck
Charles R Harris wrote:
On Sun, May 31, 2009 at 7:18 PM, David Cournapeau <david@ar.media.kyoto-u.ac.jp <mailto:david@ar.media.kyoto-u.ac.jp>> wrote:
Charles R Harris wrote: > > > On Sun, May 31, 2009 at 11:54 AM, rob steed <rjsteed@talk21.com <mailto:rjsteed@talk21.com> > <mailto:rjsteed@talk21.com <mailto:rjsteed@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 9:08 PM, David Cournapeau < david@ar.media.kyoto-u.ac.jp> wrote:
Charles R Harris wrote:
On Sun, May 31, 2009 at 7:18 PM, David Cournapeau <david@ar.media.kyoto-u.ac.jp <mailto:david@ar.media.kyoto-u.ac.jp>> wrote:
Charles R Harris wrote: > > > On Sun, May 31, 2009 at 11:54 AM, rob steed <rjsteed@talk21.com <mailto:rjsteed@talk21.com> > <mailto:rjsteed@talk21.com <mailto:rjsteed@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'm not clear on how correlate should apply to any of 'SUV' but it might be worth having it work for objects.
It still bothers me that correlate does not conjugate the second argument for complex arrays...
It bothers me also... Chuck
Charles R Harris wrote:
On Sun, May 31, 2009 at 9:08 PM, David Cournapeau <david@ar.media.kyoto-u.ac.jp <mailto:david@ar.media.kyoto-u.ac.jp>> wrote:
Charles R Harris wrote: > > > On Sun, May 31, 2009 at 7:18 PM, David Cournapeau > <david@ar.media.kyoto-u.ac.jp <mailto:david@ar.media.kyoto-u.ac.jp> <mailto:david@ar.media.kyoto-u.ac.jp <mailto:david@ar.media.kyoto-u.ac.jp>>> > wrote: > > Charles R Harris wrote: > > > > > > On Sun, May 31, 2009 at 11:54 AM, rob steed <rjsteed@talk21.com <mailto:rjsteed@talk21.com> > <mailto:rjsteed@talk21.com <mailto:rjsteed@talk21.com>> > > <mailto:rjsteed@talk21.com <mailto:rjsteed@talk21.com> <mailto:rjsteed@talk21.com <mailto:rjsteed@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
On Mon, Jun 1, 2009 at 00:05, David Cournapeau <david@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
On Mon, Jun 1, 2009 at 10:35 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Mon, Jun 1, 2009 at 00:05, David Cournapeau <david@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.
Maybe we need a new function. But what to call it? Chuck
On Mon, Jun 1, 2009 at 11:48 AM, Charles R Harris <charlesr.harris@gmail.com
wrote:
On Mon, Jun 1, 2009 at 10:35 AM, Robert Kern <robert.kern@gmail.com>wrote:
On Mon, Jun 1, 2009 at 00:05, David Cournapeau <david@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.
Maybe we need a new function. But what to call it?
How about introducing acorrelate and deprecating the old version? Chuck
On Mon, Jun 1, 2009 at 13:44, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jun 1, 2009 at 11:48 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jun 1, 2009 at 10:35 AM, Robert Kern <robert.kern@gmail.com> wrote:
On Mon, Jun 1, 2009 at 00:05, David Cournapeau <david@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.
Maybe we need a new function. But what to call it?
How about introducing acorrelate and deprecating the old version?
Sure. -- 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
Charles R Harris wrote:
On Mon, Jun 1, 2009 at 11:48 AM, Charles R Harris <charlesr.harris@gmail.com <mailto:charlesr.harris@gmail.com>> wrote:
On Mon, Jun 1, 2009 at 10:35 AM, Robert Kern <robert.kern@gmail.com <mailto:robert.kern@gmail.com>> wrote:
On Mon, Jun 1, 2009 at 00:05, David Cournapeau <david@ar.media.kyoto-u.ac.jp <mailto:david@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.
Maybe we need a new function. But what to call it?
How about introducing acorrelate and deprecating the old version?
This does not solve the C function problem (PyArray_Correlate). The easy solution would be to keep the current C version, deal with the problem in python for acorrelate for the time being, and replace the old C function with the 'correct' one once we remove the deprecated correlate ? David
On Mon, Jun 1, 2009 at 22:33, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
Charles R Harris wrote:
On Mon, Jun 1, 2009 at 11:48 AM, Charles R Harris <charlesr.harris@gmail.com <mailto:charlesr.harris@gmail.com>> wrote:
On Mon, Jun 1, 2009 at 10:35 AM, Robert Kern <robert.kern@gmail.com <mailto:robert.kern@gmail.com>> wrote:
On Mon, Jun 1, 2009 at 00:05, David Cournapeau <david@ar.media.kyoto-u.ac.jp <mailto:david@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.
Maybe we need a new function. But what to call it?
How about introducing acorrelate and deprecating the old version?
This does not solve the C function problem (PyArray_Correlate). The easy solution would be to keep the current C version, deal with the problem in python for acorrelate for the time being, and replace the old C function with the 'correct' one once we remove the deprecated correlate ?
No, you do the same thing at the C level. -- 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
Robert Kern wrote:
This does not solve the C function problem (PyArray_Correlate). The easy solution would be to keep the current C version, deal with the problem in python for acorrelate for the time being, and replace the old C function with the 'correct' one once we remove the deprecated correlate ?
No, you do the same thing at the C level.
Ok, David
On Tue, Jun 2, 2009 at 12:37 PM, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
Robert Kern wrote:
This does not solve the C function problem (PyArray_Correlate). The easy solution would be to keep the current C version, deal with the problem in python for acorrelate for the time being, and replace the old C function with the 'correct' one once we remove the deprecated correlate ?
No, you do the same thing at the C level.
Done in r7031 - correlate/PyArray_Correlate should be unchanged, and acorrelate/PyArray_Acorrelate implement the conventional definitions, David
On Tue, Jun 2, 2009 at 11:36 AM, David Cournapeau <cournape@gmail.com> wrote:
Done in r7031 - correlate/PyArray_Correlate should be unchanged, and acorrelate/PyArray_Acorrelate implement the conventional definitions,
I don't know if it's been discussed before but while people are thinking about/changing correlate I thought I'd like to request as a user a matlab style xcorr function (basically with the functionality of the matlab version). I don't know if this is a deliberate emission, but it is often one of the first things my colleagues try when I get them using Python, and as far as I know there isn't really a good answer. There is xcorr in pylab, but it isn't vectorised like xcorr from matlab... Cheers Robin
Robin wrote:
On Tue, Jun 2, 2009 at 11:36 AM, David Cournapeau <cournape@gmail.com> wrote:
Done in r7031 - correlate/PyArray_Correlate should be unchanged, and acorrelate/PyArray_Acorrelate implement the conventional definitions,
I don't know if it's been discussed before but while people are thinking about/changing correlate I thought I'd like to request as a user a matlab style xcorr function (basically with the functionality of the matlab version).
I don't know if this is a deliberate emission, but it is often one of the first things my colleagues try when I get them using Python, and as far as I know there isn't really a good answer. There is xcorr in pylab, but it isn't vectorised like xcorr from matlab...
There is one in the talkbox scikit: http://github.com/cournape/talkbox/blob/202135a9d848931ebd036b97302f1e10d748... It uses the fft, and bonus point, the file is independent of the rest of toolbox. There is another version which uses direct implementation (this is faster if you need only a few lags, and it takes less memory too). David
On Tue, Jun 2, 2009 at 5:59 AM, David Cournapeau < david@ar.media.kyoto-u.ac.jp> wrote:
Robin wrote:
On Tue, Jun 2, 2009 at 11:36 AM, David Cournapeau <cournape@gmail.com> wrote:
Done in r7031 - correlate/PyArray_Correlate should be unchanged, and acorrelate/PyArray_Acorrelate implement the conventional definitions,
I don't know if it's been discussed before but while people are thinking about/changing correlate I thought I'd like to request as a user a matlab style xcorr function (basically with the functionality of the matlab version).
I don't know if this is a deliberate emission, but it is often one of the first things my colleagues try when I get them using Python, and as far as I know there isn't really a good answer. There is xcorr in pylab, but it isn't vectorised like xcorr from matlab...
There is one in the talkbox scikit:
http://github.com/cournape/talkbox/blob/202135a9d848931ebd036b97302f1e10d748...
It uses the fft, and bonus point, the file is independent of the rest of toolbox. There is another version which uses direct implementation (this is faster if you need only a few lags, and it takes less memory too).
I'd be +1 on including something like this (provided it expanded to include complex-valued data). I think it's a real need, since everyone seems to keep rolling their own. I had to write my own just so that I can calculate a few lags in a vectorized fashion. Ryan -- Ryan May Graduate Research Assistant School of Meteorology University of Oklahoma
On Tue, Jun 2, 2009 at 10:56 PM, Ryan May <rmay31@gmail.com> wrote:
On Tue, Jun 2, 2009 at 5:59 AM, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
Robin wrote:
On Tue, Jun 2, 2009 at 11:36 AM, David Cournapeau <cournape@gmail.com> wrote:
Done in r7031 - correlate/PyArray_Correlate should be unchanged, and acorrelate/PyArray_Acorrelate implement the conventional definitions,
I don't know if it's been discussed before but while people are thinking about/changing correlate I thought I'd like to request as a user a matlab style xcorr function (basically with the functionality of the matlab version).
I don't know if this is a deliberate emission, but it is often one of the first things my colleagues try when I get them using Python, and as far as I know there isn't really a good answer. There is xcorr in pylab, but it isn't vectorised like xcorr from matlab...
There is one in the talkbox scikit:
http://github.com/cournape/talkbox/blob/202135a9d848931ebd036b97302f1e10d748...
It uses the fft, and bonus point, the file is independent of the rest of toolbox. There is another version which uses direct implementation (this is faster if you need only a few lags, and it takes less memory too).
I'd be +1 on including something like this (provided it expanded to include complex-valued data). I think it's a real need, since everyone seems to keep rolling their own. I had to write my own just so that I can calculate a few lags in a vectorized fashion.
The code in talkbox is not good enough for scipy. I made an attempt for scipy.signal here: http://github.com/cournape/scipy3/blob/b004d17d824f1c03921d9663207ee40adadc5... It is reasonably fast when only a few lags are needed, both double and complex double are supported, and it works on arbitrary axis and lags. Other precisions should be easy to add, but I think I need to extend the numpy code generators to support cython sources to avoid code duplication. Does that fill your need ? cheers, David
On Thu, Jun 4, 2009 at 5:14 AM, David Cournapeau <cournape@gmail.com> wrote:
On Tue, Jun 2, 2009 at 10:56 PM, Ryan May <rmay31@gmail.com> wrote:
On Tue, Jun 2, 2009 at 5:59 AM, David Cournapeau <david@ar.media.kyoto-u.ac.jp> wrote:
Robin wrote:
On Tue, Jun 2, 2009 at 11:36 AM, David Cournapeau <cournape@gmail.com
wrote:
Done in r7031 - correlate/PyArray_Correlate should be unchanged, and acorrelate/PyArray_Acorrelate implement the conventional definitions,
I don't know if it's been discussed before but while people are thinking about/changing correlate I thought I'd like to request as a user a matlab style xcorr function (basically with the functionality of the matlab version).
I don't know if this is a deliberate emission, but it is often one of the first things my colleagues try when I get them using Python, and as far as I know there isn't really a good answer. There is xcorr in pylab, but it isn't vectorised like xcorr from matlab...
There is one in the talkbox scikit:
http://github.com/cournape/talkbox/blob/202135a9d848931ebd036b97302f1e10d748...
It uses the fft, and bonus point, the file is independent of the rest of toolbox. There is another version which uses direct implementation (this is faster if you need only a few lags, and it takes less memory too).
I'd be +1 on including something like this (provided it expanded to include complex-valued data). I think it's a real need, since everyone seems to keep rolling their own. I had to write my own just so that I can calculate a few lags in a vectorized fashion.
The code in talkbox is not good enough for scipy. I made an attempt for scipy.signal here:
http://github.com/cournape/scipy3/blob/b004d17d824f1c03921d9663207ee40adadc5...
It is reasonably fast when only a few lags are needed, both double and complex double are supported, and it works on arbitrary axis and lags. Other precisions should be easy to add, but I think I need to extend the numpy code generators to support cython sources to avoid code duplication.
Does that fill your need ?
It would fill mine. Would it make sense to make y default to x, so that you can use xcorr to do the autocorrelation as: xcorr(x) ? Ryan -- Ryan May Graduate Research Assistant School of Meteorology University of Oklahoma
participants (7)
-
Charles R Harris -
David Cournapeau -
David Cournapeau -
rob steed -
Robert Kern -
Robin -
Ryan May