[Numpy-discussion] Cython-based OpenMP-accelerated quartic polynomial solver

[Juha Jeronen]

Hi all,

Thanks Jon and Ryan for the suggestions. Both asanyarray() or 
atleast_1d() sound good.

There's the technical detail that Cython needs to know the object type 
(e.g. in the parameters to quartic_z_array()), so I think atleast_1d() 
may be safer. I've taken this option for now.

This simplified the code somewhat. The *_scalar() functions have been 
removed, as they are no longer needed. The *_array() versions have been 
renamed, removing the _array suffix.

The return values have stayed as they were - if there is only one 
problem to solve, the singleton dimension is dropped, and otherwise a 2D 
array is returned.

(The exception is linear(), which does not need the second dimension, 
since the solution for each problem is unique. It will return a scalar 
in the case of a single problem, or a 1D array in the case of multiple 
problems.)


I've pushed the new version. It's available from the same repository:

https://yousource.it.jyu.fi/jjrandom2/miniprojects/trees/master/misc/polysolve_for_numpy


Other comments? The next step?

  -J


P.S.:

I'm not sure how exact the object type must be - i.e. whether Cython 
wants to know that the object stores its data somewhat like an ndarray, 
or that its C API exactly matches that of an ndarray.

In Cython there are some surprising details regarding this kind of 
things, such as the ctypedef'ing of scalar types. For example, see the 
comments about complex support near the beginning of polysolve2.pyx, and 
the commented-out SSE2 intrinsics experiment in 
https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/tworods.pyx 
.

In the latter, it was somewhat tricky to get Cython to recognize __m128d 
- turns out it's close enough for Cython to know that it behaves like a 
double, although it actually contains two doubles.

(Note that these ctypedefs never end up in the generated C code; Cython 
uses them as extra context information when mapping the Python code into C.)

(And no need to worry, I'm not planning to put any SSE into polysolve2 :) )



On 02.10.2015 17:18, Ryan May wrote:
> numpy.asanyarray() would be my preferred goto, as it will leave 
> subclasses of ndarray untouched; asarray() and atleast_1d() force 
> ndarray. It's nice to do the whenever possible.
>
> Ryan
>
> On Fri, Oct 2, 2015 at 6:52 AM, Slavin, Jonathan 
> <jslavin at cfa.harvard.edu <mailto:jslavin at cfa.harvard.edu>> wrote:
>
>     ​ Personally I like atleast_1d, which will convert a scalar into a
>     1d array but will leave arrays untouched (i.e. won't change the
>     dimensions.  Not sure what the advantages/disadvantages are
>     relative to asarray.
>
>     Jon​
>
>
>     On Fri, Oct 2, 2015 at 7:05 AM,
>     <numpy-discussion-request at scipy.org
>     <mailto:numpy-discussion-request at scipy.org>> wrote:
>
>         From: Juha Jeronen <juha.jeronen at jyu.fi
>         <mailto:juha.jeronen at jyu.fi>>
>         To: Discussion of Numerical Python <numpy-discussion at scipy.org
>         <mailto:numpy-discussion at scipy.org>>
>         Cc:
>         Date: Fri, 2 Oct 2015 13:31:47 +0300
>         Subject: Re: [Numpy-discussion] Cython-based
>         OpenMP-accelerated quartic polynomial solver
>
>         On 02.10.2015 13:07, Daπid wrote:
>>
>>         On 2 October 2015 at 11:58, Juha Jeronen <juha.jeronen at jyu.fi
>>         <mailto:juha.jeronen at jyu.fi>> wrote:
>>
>>
>>
>>             First version done and uploaded:
>>
>>             https://yousource.it.jyu.fi/jjrandom2/miniprojects/trees/master/misc/polysolve_for_numpy
>>
>>
>>         Small comment: now you are checking if the input is a scalar
>>         or a ndarray, but it should also accept any array-like. If I
>>         pass a list, I expect it to work, internally converting it
>>         into an array.
>
>         Good catch.
>
>         Is there an official way to test for array-likes? Or should I
>         always convert with asarray()? Or something else?
>
>
>          -J
>
>
>
>
>
>     -- 
>     ________________________________________________________
>     Jonathan D. Slavin Harvard-Smithsonian CfA
>     jslavin at cfa.harvard.edu <mailto:jslavin at cfa.harvard.edu>    60
>     Garden Street, MS 83
>     phone: (617) 496-7981 <tel:%28617%29%20496-7981>       Cambridge,
>     MA 02138-1516
>     cell: (781) 363-0035 <tel:%28781%29%20363-0035>             USA
>     ________________________________________________________
>
>
>     _______________________________________________
>     NumPy-Discussion mailing list
>     NumPy-Discussion at scipy.org <mailto:NumPy-Discussion at scipy.org>
>     https://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
>
>
> -- 
> Ryan May
>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20151006/a056133d/attachment.html>