[Numpy-discussion] failure to register ufunc loops for user defined types

[Geoffrey Irving]

On Mon, Dec 5, 2011 at 6:59 AM, Charles R Harris
<charlesr.harris at gmail.com> wrote:
> Hi Geoffrey,
>
> On Mon, Dec 5, 2011 at 12:37 AM, Geoffrey Irving <irving at naml.us> wrote:
>>
>> On Sun, Dec 4, 2011 at 6:45 PM, Charles R Harris
>> <charlesr.harris at gmail.com> wrote:
>> >
>> >
>> > On Sun, Dec 4, 2011 at 6:59 PM, Geoffrey Irving <irving at naml.us> wrote:
>> >>
>> >> On Sun, Dec 4, 2011 at 5:18 PM, Charles R Harris
>> >> <charlesr.harris at gmail.com> wrote:
>> >> >
>> >> >
>> >> > On Sun, Dec 4, 2011 at 5:41 PM, Geoffrey Irving <irving at naml.us>
>> >> > wrote:
>> >> >>
>> >> >> This may be the problem.  Simple diffs are pleasant.  I'm guessing
>> >> >> this code doesn't get a lot of testing.  Glad it's there, though!
>> >> >>
>> >> >> Geoffrey
>> >> >>
>> >> >> diff --git a/numpy/core/src/umath/ufunc_type_resolution.c
>> >> >> b/numpy/core/src/umath/ufunc_type_resolution.c
>> >> >> index 0d6cf19..a93eda1 100644
>> >> >> --- a/numpy/core/src/umath/ufunc_type_resolution.c
>> >> >> +++ b/numpy/core/src/umath/ufunc_type_resolution.c
>> >> >> @@ -1866,7 +1866,7 @@ linear_search_type_resolver(PyUFuncObject
>> >> >> *self,
>> >> >>             case -1:
>> >> >>                 return -1;
>> >> >>             /* A loop was found */
>> >> >> -            case 1:
>> >> >> +            case 0:
>> >> >>                 return 0;
>> >> >>         }
>> >> >>     }
>> >> >>
>> >> >
>> >> > Heh. Can you verify that this fixes the problem? That function is
>> >> > only
>> >> > called once  and its return value is passed up the chain, but the
>> >> > documented
>> >> > return values of that calling function are -1, 0. So the
>> >> > documentation
>> >> > needs
>> >> > to be changed if this is the right thing to do.
>> >>
>> >> Actually, that patch was wrong, since
>> >> linear_search_userloop_type_resolver needs to return three values
>> >> (error, not-found, success).  A better patch follows.  I can confirm
>> >> that this gets me further, but I get other failures down the line, so
>> >> more fixes may follow.  I'll push the branch with all my fixes for
>> >> convenience once I have everything working.
>> >>
>> >> > Speaking of tests... I was wondering if you could be talked into
>> >> > putting
>> >> > together a simple user type for including in the tests?
>> >>
>> >> Yep, though likely not for a couple weeks.  If there's interest, I
>> >> could also be convinced to sanitize my entire rational class so you
>> >> could include that directly.  Currently it's both C++ and uses some
>> >> gcc specific features like __int128_t.  Basically it's
>> >> numerator/denominator, where both are 64 bit integers, and an
>> >> OverflowError is thrown if anything can't be represented as such
>> >> (possibly a different exception would be better in cases like
>> >> (1<<64)/((1<<64)+1)).  It would be easy to generalize it to rational32
>> >> vs. rational64 as well.
>> >>
>> >> If you want tests but not rational, it would be straightforward to
>> >> strip what I have down to a bare bones test case.
>> >>
>> >
>> > We'll see how much interest there is. If it becomes official you may get
>> > more feedback on features. There are some advantages to having some user
>> > types in numpy. One is that otherwise they tend to get lost, another is
>> > that
>> > having a working example or two provides a templates for others to work
>> > from, and finally they provide test material. Because official user
>> > types
>> > aren't assigned anywhere there might also be some conflicts. Maybe
>> > something
>> > like an extension types module would be a way around that. In any case,
>> > I
>> > think both rational numbers and quaternions would be useful to have and
>> > I
>> > hope there is some discussion of how to do that. Rationals may be a bit
>> > trickier than quaternions though, as usually they are used to provide
>> > exact
>> > arithmetic without concern for precision. I don't know how restrictive
>> > the
>> > 64 bit limitation will be in practice. What are you using them for?
>>
>> I'm using them for frivolous analysis of poker Nash equilibria.  I'll
>> let others decide if it has any non-toy uses.  64 bits seems to be
>> enough for me, though it's possible that I'll run in trouble with
>> other examples.  It still exact, though, in the sense that it throws
>> an exception rather than doing anything weird if it overflows.  And it
>> has the key advantage of being orders of magnitude faster than object
>> arrays of Fractions.
>>
>> Back to the bugs: here's a branch with all the changes I needed to get
>> rational arithmetic to work:
>>
>>    https://github.com/girving/numpy
>>
>> I discovered two more after the last email.  One is another simple 0
>> vs. 1 bug, and another is somewhat optional:
>>
>> commit 730b05a892371d6f18d9317e5ae6dc306c0211b0
>> Author: Geoffrey Irving <irving at naml.us>
>> Date:   Sun Dec 4 20:03:46 2011 -0800
>>
>>    After loops, check for PyErr_Occurred() even if needs_api is 0
>>
>>    For certain types of user defined classes, casting and ufunc loops
>>    normally run without the Python API, but occasionally need to throw
>>    an error.  Currently we assume that !needs_api means no error occur.
>>    However, the fastest way to implement such loops is to run without
>>    the GIL normally and use PyGILState_Ensure/Release if an error occurs.
>>
>>    In order to support this usage pattern, change all post-loop checks
>> from
>>
>>        needs_api && PyErr_Occurred()
>>
>>    to simply
>>
>>        PyErr_Occurred()
>>
>
> Thanks. Could you put this work into a separate branch, say fixuserloops,
> and enter a pull request? It's best not to work in master and a pull request
> make things easier to discuss and merge.

Done:

    https://github.com/numpy/numpy/pull/175

Thanks,
Geoffrey