[Numpy-discussion] Problem with correlate

[David Cournapeau]

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