Re: [Numpy-discussion] composition of the steering council (was Re: Governance model request)

[Nathaniel Smith]

On Fri, Sep 25, 2015 at 7:15 AM, Thomas Caswell  wrote:
> To respond to the devils advocate:
>
>  Creating this organizational framework is a one time boot-strapping event.
> You could use wording like "The initial council will include those who have
> made significant contributions to numpy in the past and want to be on it" or
> "The initial council will be constructed by invitation by Nathaniel and
> Chuck".  More objective criteria should be used going forward, but in terms
> of getting things spun up quickly doing things by fiat is probably ok.  I am
> not even sure that the method by which the initial group is formed needs to
> go into the governing document.

The problem is that according to the current text, not only is Travis
ineligible to join the council (it's a little weird to put people on
the seed council who wouldn't be eligible to join it normally, but
okay, sure) -- he's not even eligible to stay on the council once he
joins. So we would need to change the text no matter what.

Which we can do, if we decide that that's what we need to do to
accomplish what we want. It's our text, after all. I think it's
extremely important though that what we actually do, and what we write
down saying we will do, somehow match. Otherwise this whole exercise
has no point.

-n

-- 
Nathaniel J. Smith -- http://vorpus.org
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

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

[Juha Jeronen]

Hi all,

I recently developed a Cython-based, OpenMP-accelerated quartic (and 
cubic, quadratic) polynomial solver to address a personal research need 
for quickly solving a very large number of independent low-degree 
polynomial equations for both real and complex coefficients.


For example, on an Intel i7-4710MQ, running with 8 threads this code 
solves approximately 1.2e7 quartic polynomial equations per second. 
(With OMP_NUM_THREADS=4 the solver gives approximately 8.9e6 solutions 
per second, so it seems HyperThreading in this case helps about 30%.)


The algorithms for cubics and quadratics come from Numerical Recipes 
(3rd ed.), and the quartic problem is internally reduced to a cubic and 
two quadratics, using well-known standard tricks.



Since to my understanding, for solving polynomial equations NumPy only 
provides roots(), which works one problem at a time, I'm writing to 
inquire whether there would be interest to include this functionality to 
NumPy, if I contribute the code (and clean it up for integration)?



I have some previous open source contribution experience. I have 
authored the IVTC and Phosphor deinterlacers for VLC, and a modular 
postprocessing filter framework for the Panda3D game engine (posted at 
the forums on panda3d.org, projected for version 1.10).


Currently the polynomial solver is available in a git repository hosted 
by our university, the relevant files being:


https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/polysolve2.pyx
https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/polysolve.py
https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/polysolve2example.py
https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/setup.py
https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/compile.sh

I'm working on a research grant, which allows me to hold the copyright; 
license is BSD.


Polysolve2 is the fast Cython+OpenMP version, while polysolve is an 
earlier (slower) experimental version using only NumPy array operations. 
The example module contains a usage example, and setup (in the standard 
manner) instructs distutils and Cython as to how to build polysolve2 
(including setting the relevant math flags and enabling OpenMP).


(The setup script includes also some stuff specific to my solver for the 
Ziegler problem; basically, the "tworods" module can be ignored. I 
apologize for the inconvenience.)



Thoughts?


Best regards,

 -J

-
Juha Jeronen, Ph.D.
juha.jero...@jyu.fi
University of Jyväskylä
Department of Mathematical Information Technology
-

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Freddy Rietdijk]

I wouldn't know of any valid output when applying the sign function to NaN.
Therefore, I think it is correct to return a ValueError. Furthermore, I
would prefer such an error over just returning NaN since it helps you
locating where NaN is generated.

On Tue, Sep 29, 2015 at 5:13 PM, Charles R Harris  wrote:

> Hi All,
>
> Due to a recent commit, Numpy master now raises an error when applying the
> sign function to an object array containing NaN. Other options may be
> preferable, returning NaN for instance, so I would like to open the topic
> for discussion on the list.
>
> Thoughts?
>
> Chuck
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Allan Haldane]

On 09/29/2015 11:39 AM, josef.p...@gmail.com wrote:
> 
> 
> On Tue, Sep 29, 2015 at 11:25 AM, Anne Archibald  > wrote:
> 
> IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point
> arrays. Why should it be different for object arrays?
> 
> Anne
> 
> P.S. If you want exceptions when NaNs appear, that's what np.seterr
> is for. -A
> 
> 
> 
> I also think NaN should be treated the same way as floating point
> numbers (whatever that is). Otherwise it is difficult to remember when
> nan is essentially a float dtype or another dtype.  
> (given that float is the smallest dtype that can hold a nan)
> 

Note that I've reimplemented np.sign for object arrays along these lines
in this open PR:
https://github.com/numpy/numpy/pull/6320

That PR recursively uses the np.sign ufunc to evaluate object arrays
containing float and complex numbers. This way the behavior on object
arrays is identical to float/complex arrays.

Here is what the np.sign ufunc does (for arbitrary x):

np.sign(np.nan)  ->  nan
np.sign(complex(np.nan, x))  ->  complex(nan, 0)
np.sign(complex(x, np.nan))  ->  complex(nan, 0)

Allan
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

[Numpy-discussion] Sign of NaN

[Charles R Harris]

Hi All,

Due to a recent commit, Numpy master now raises an error when applying the
sign function to an object array containing NaN. Other options may be
preferable, returning NaN for instance, so I would like to open the topic
for discussion on the list.

Thoughts?

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] composition of the steering council (was Re: Governance model request)

[Travis Oliphant]

Thanks for the candid discussion and for expressing concerns freely.

I think Nathaniel's "parenting" characterization of NumPy from me is pretty
accurate.I do feel a responsibility for the *stuff* that's out there,
and that is what drives me.   I do see the many contributions from others
and really learn from them as well.

I have seen conversations on this list and others have characterizations of
history that I don't agree with which affects decisions that are being made
--- and so I feel compelled to try and share my view.

I'm in a situation now where at least for 6 months or so I can help with
NumPy more than I have been able to for 7 years.

Focusing on the initial governance text, my issues are that

1) 1 year of inactivity to be removed from the council is too little for a
long-running project like NumPy --- somewhere between 2 and 4 years would
be more appropriate.   I suppose 1 year of inactivity is fine if that is
defined only as "failure to vote on matters before the council"

2) The seed council should not just be recent contributors but should
include as many people as are willing to help who have a long history with
the project.

3) I think people who contribute significantly generally should be able to
re-join the steering council more easily than "going through the 1-year
vetting process" again --- they would have to be approved by the current
steering council but it should not be automatically disallowed (thus
requiring the equivalent of an amendment to change it).

I applaud the fact that the steering council will not be and should not be
used except when absolutely necessary and for limited functions.

Thanks,

-Travis




On Tue, Sep 29, 2015 at 4:06 AM, Nathaniel Smith  wrote:

> On Fri, Sep 25, 2015 at 7:15 AM, Thomas Caswell 
> wrote:
> > To respond to the devils advocate:
> >
> >  Creating this organizational framework is a one time boot-strapping
> event.
> > You could use wording like "The initial council will include those who
> have
> > made significant contributions to numpy in the past and want to be on
> it" or
> > "The initial council will be constructed by invitation by Nathaniel and
> > Chuck".  More objective criteria should be used going forward, but in
> terms
> > of getting things spun up quickly doing things by fiat is probably ok.
> I am
> > not even sure that the method by which the initial group is formed needs
> to
> > go into the governing document.
>
> The problem is that according to the current text, not only is Travis
> ineligible to join the council (it's a little weird to put people on
> the seed council who wouldn't be eligible to join it normally, but
> okay, sure) -- he's not even eligible to stay on the council once he
> joins. So we would need to change the text no matter what.
>
> Which we can do, if we decide that that's what we need to do to
> accomplish what we want. It's our text, after all. I think it's
> extremely important though that what we actually do, and what we write
> down saying we will do, somehow match. Otherwise this whole exercise
> has no point.
>
> -n
>
> --
> Nathaniel J. Smith -- http://vorpus.org
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>



-- 

*Travis Oliphant*
*Co-founder and CEO*


@teoliphant
512-222-5440
http://www.continuum.io
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[josef.pktd]

On Tue, Sep 29, 2015 at 11:25 AM, Anne Archibald 
wrote:

> IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point arrays.
> Why should it be different for object arrays?
>
> Anne
>
> P.S. If you want exceptions when NaNs appear, that's what np.seterr is
> for. -A
>


I also think NaN should be treated the same way as floating point numbers
(whatever that is). Otherwise it is difficult to remember when nan is
essentially a float dtype or another dtype.
(given that float is the smallest dtype that can hold a nan)

Josef



>
> On Tue, Sep 29, 2015 at 5:18 PM Freddy Rietdijk 
> wrote:
>
>> I wouldn't know of any valid output when applying the sign function to
>> NaN. Therefore, I think it is correct to return a ValueError. Furthermore,
>> I would prefer such an error over just returning NaN since it helps you
>> locating where NaN is generated.
>>
>> On Tue, Sep 29, 2015 at 5:13 PM, Charles R Harris <
>> charlesr.har...@gmail.com> wrote:
>>
>>> Hi All,
>>>
>>> Due to a recent commit, Numpy master now raises an error when applying
>>> the sign function to an object array containing NaN. Other options may be
>>> preferable, returning NaN for instance, so I would like to open the topic
>>> for discussion on the list.
>>>
>>> Thoughts?
>>>
>>> Chuck
>>>
>>> ___
>>> NumPy-Discussion mailing list
>>> NumPy-Discussion@scipy.org
>>> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>>>
>>>
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Anne Archibald]

IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point arrays.
Why should it be different for object arrays?

Anne

P.S. If you want exceptions when NaNs appear, that's what np.seterr is for.
-A

On Tue, Sep 29, 2015 at 5:18 PM Freddy Rietdijk 
wrote:

> I wouldn't know of any valid output when applying the sign function to
> NaN. Therefore, I think it is correct to return a ValueError. Furthermore,
> I would prefer such an error over just returning NaN since it helps you
> locating where NaN is generated.
>
> On Tue, Sep 29, 2015 at 5:13 PM, Charles R Harris <
> charlesr.har...@gmail.com> wrote:
>
>> Hi All,
>>
>> Due to a recent commit, Numpy master now raises an error when applying
>> the sign function to an object array containing NaN. Other options may be
>> preferable, returning NaN for instance, so I would like to open the topic
>> for discussion on the list.
>>
>> Thoughts?
>>
>> Chuck
>>
>> ___
>> NumPy-Discussion mailing list
>> NumPy-Discussion@scipy.org
>> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Antoine Pitrou]

On Tue, 29 Sep 2015 09:13:15 -0600
Charles R Harris  wrote:
> 
> Due to a recent commit, Numpy master now raises an error when applying the
> sign function to an object array containing NaN. Other options may be
> preferable, returning NaN for instance, so I would like to open the topic
> for discussion on the list.

None for example? float('nan') may be a bit weird amongst e.g. an array
of Decimals.

Regards

Antoine.


___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Joe Kington]

On Tue, Sep 29, 2015 at 11:14 AM, Antoine Pitrou 
wrote:

>
> None for example? float('nan') may be a bit weird amongst e.g. an array
> of Decimals


The downside to `None` is that it's one more thing to check for and makes
object arrays an even weirder edge case.  (Incidentally, Decimal does have
its own non-float NaN which throws a whole different wrench in the works. `
np.sign(Decimal('NaN'))` is going to raise an error no matter what.)

A float (or numpy) NaN makes more sense to return for mixed datatypes than
None does, in my opinion. At least then one can use `isfinite`, etc to
check while `np.isfinite(None)` will raise an error.  Furthermore, if
there's at least one floating point NaN in the object array, getting a
float NaN out makes sense.

Just my $0.02, anyway.
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Charles R Harris]

On Tue, Sep 29, 2015 at 9:25 AM, Anne Archibald  wrote:

> IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point arrays.
> Why should it be different for object arrays?
>

What about non-numeric objects in general ?



Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

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

[Matti Picus]

  
  
The algorithms for cubics and quadratics come from Numerical Recipes 
(3rd ed.), and the quartic problem is internally reduced to a cubic and 
two quadratics, using well-known standard tricks.Nice, wll documented code.
Just to be sure you are on safe ground, you took the *algorithms* but no *code* from Numerical Recipes, right? They are pretty clear about their licensing, see
http://www.nr.com/licenses/redistribute.html
This bit from the site might be relevant:
"We never give permission for Numerical Recipes source code to be posted on any public server, or distributed with any freeware or shareware package."
Matti

  

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Sebastian Berg]

On Di, 2015-09-29 at 11:16 -0700, Nathaniel Smith wrote:
> On Sep 29, 2015 8:25 AM, "Anne Archibald"  wrote:
> >
> > IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point
> arrays. Why should it be different for object arrays?
> 
> The argument for doing it this way would be that arbitrary python
> objects don't have a sign, and the natural way to implement something
> like np.sign's semantics using only the "object" interface is
> 
> if obj < 0:
> return -1
> elif obj > 0:
> return 1
> elif obj == 0:
> return 0
> else:
> raise
> 
> In general I'm not a big fan of trying to do all kinds of guessing
> about how to handle random objects in object arrays, the kind that
> ends up with a big chain of type checks and fallback behaviors. Pretty
> soon we find ourselves trying to extend the language with our own
> generic dispatch system for arbitrary python types, just for object
> arrays. (The current hack where for object arrays np.log will try
> calling obj.log() is particularly horrible. There is no rule in python
> that "log" is a reserved method name for "logarithm" on arbitrary
> objects. Ditto for the other ufuncs that implement this hack.)
> 
> Plus we hope that many use cases for object arrays will soon be
> supplanted by better dtype support, so now may not be the best time to
> invest heavily in making object arrays complicated and powerful. 
> 

I have the little dream here that what could happen is that we create a
PyFloatDtype kind of thing (it is a bit different from our float because
it would always convert back to a python float and maybe raises more
errors), which "registers" with the dtype system in that it says "I know
how to handle python floats and store them in an array and provide ufunc
implementations for it".

Then, the "object" dtype ufuncs would try to call the ufunc on each
element, including "conversion". They would find a "float", since it is
not an array-like container, they interpret it as a PyFloatDtype scalar
and call the scalars ufunc (the PyFloatDtype scalar would be a python
float).

Of course likely I am thinking down the wrong road, but if you want e.g.
an array of Decimals, you need some way to tell that numpy as a
PyDecimalDtype.
Now "object" would possibly be just a fallback to mean "figure out what
to use for each element". It would be a bit slower, but it would work
very generally, because numpy would not impose limits as such.

- Sebastian


> OTOH sometimes practicality beats purity, and at least object arrays
> are already kinda cordoned off from the rest of the system, so I don't
> feel as strongly as if we were talking about core functionality.
> 
> ...is there a compelling reason to even support np.sign on object
> arrays? This seems pretty far into the weeds, and that tends to lead
> to poor intuition and decision making.
> 
> -n
> 
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion



signature.asc
Description: This is a digitally signed message part
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Allan Haldane]

On 09/29/2015 02:16 PM, Nathaniel Smith wrote:
> On Sep 29, 2015 8:25 AM, "Anne Archibald"  > wrote:
>>
>> IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point
> arrays. Why should it be different for object arrays?
> 
> The argument for doing it this way would be that arbitrary python
> objects don't have a sign, and the natural way to implement something
> like np.sign's semantics using only the "object" interface is
> 
> if obj < 0:
> return -1
> elif obj > 0:
> return 1
> elif obj == 0:
> return 0
> else:
> raise
> 
> In general I'm not a big fan of trying to do all kinds of guessing about
> how to handle random objects in object arrays, the kind that ends up
> with a big chain of type checks and fallback behaviors. Pretty soon we
> find ourselves trying to extend the language with our own generic
> dispatch system for arbitrary python types, just for object arrays. (The
> current hack where for object arrays np.log will try calling obj.log()
> is particularly horrible. There is no rule in python that "log" is a
> reserved method name for "logarithm" on arbitrary objects. Ditto for the
> other ufuncs that implement this hack.)
> 
> Plus we hope that many use cases for object arrays will soon be
> supplanted by better dtype support, so now may not be the best time to
> invest heavily in making object arrays complicated and powerful.

Even though I submitted the PR to make object arrays more powerful, this
makes a lot of sense to me.

Let's say we finish a new dtype system, in which (I imagine) each dtype
specifies how to calculate each ufunc elementwise for that type. What
are the remaining use cases for generic object arrays? The only one I
see is having an array with elements of different types, which seems
like a dubious idea anyway. (Nested ndarrays of varying length could be
implemented as a dtype, as could the PyFloatDtype Sebastian mentioned,
without need for a generic 'object' dtype which has to figure out how
to call ufuncs on individual objects of different type).

Allan
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Chris Barker - NOAA Federal]

One of the usecases that has sneaked in during the last few numpy versions
is that object arrays contain numerical arrays where the shapes don't add
up to a rectangular array.


I think that's the wrong way to dve that problem -- we really should have a
"proper" ragged array implementation. But is is the easiest way at this
point.

For this, and other use-cases, special casing Numpy arrays stored in object
arrays does make sense:

"If this is s a Numpy array, pass the operation through."

-CHB


In those cases being able to apply numerical operations might be useful.

But I'm +0 since I don't work with object arrays.

Josef



> -n
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Charles R Harris]

On Tue, Sep 29, 2015 at 7:31 PM, Charles R Harris  wrote:

>
>
> On Tue, Sep 29, 2015 at 6:58 PM, Chris Barker - NOAA Federal <
> chris.bar...@noaa.gov> wrote:
>
>>
>> One of the usecases that has sneaked in during the last few numpy
>> versions is that object arrays contain numerical arrays where the shapes
>> don't add up to a rectangular array.
>>
>>
>> I think that's the wrong way to dve that problem -- we really should have
>> a "proper" ragged array implementation. But is is the easiest way at this
>> point.
>>
>> For this, and other use-cases, special casing Numpy arrays stored in
>> object arrays does make sense:
>>
>> "If this is s a Numpy array, pass the operation through."
>>
>
> Because we now (development) use rich compare, the result looks like
>
> In [1]: a = ones(3)
>
> In [2]: b = array([a, -a], object)
>
> In [3]: b
> Out[3]:
> array([[1.0, 1.0, 1.0],
>[-1.0, -1.0, -1.0]], dtype=object)
>
> In [4]: sign(b)
> Out[4]:
> array([[1L, 1L, 1L],
>[-1L, -1L, -1L]], dtype=object)
>
> The function returns long integers in order to not special case Python 3.
> Hmm, wonder if we might want to change that.
>

Oops, not what was intended. In fact it raises an error

In [7]: b
Out[7]: array([array([ 1.,  1.,  1.]), array([-1., -1., -1.])],
dtype=object)

In [8]: sign(b)
---
ValueErrorTraceback (most recent call last)
 in ()
> 1 sign(b)

ValueError: The truth value of an array with more than one element is
ambiguous. Use a.any() or a.all()

Chuck



>
> Chuck
>
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

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

[Chris Barker - NOAA Federal]

This sounds pretty cool -- and I've had a use case. So it would be
nice to get into Numpy.

But: I doubt we'd want OpenMP dependence in Numpy itself.

But maybe a pure Cython non-MP version?

Are we putting Cuthon in Numpy now? I lost track.

-CHB

Sent from my iPhone

> On Sep 29, 2015, at 7:35 AM, Juha Jeronen  wrote:
>
> Hi all,
>
> I recently developed a Cython-based, OpenMP-accelerated quartic (and cubic, 
> quadratic) polynomial solver to address a personal research need for quickly 
> solving a very large number of independent low-degree polynomial equations 
> for both real and complex coefficients.
>
> For example, on an Intel i7-4710MQ, running with 8 threads this code solves 
> approximately 1.2e7 quartic polynomial equations per second. (With 
> OMP_NUM_THREADS=4 the solver gives approximately 8.9e6 solutions per second, 
> so it seems HyperThreading in this case helps about 30%.)
>
> The algorithms for cubics and quadratics come from Numerical Recipes (3rd 
> ed.), and the quartic problem is internally reduced to a cubic and two 
> quadratics, using well-known standard tricks.
>
>
> Since to my understanding, for solving polynomial equations NumPy only 
> provides roots(), which works one problem at a time, I'm writing to inquire 
> whether there would be interest to include this functionality to NumPy, if I 
> contribute the code (and clean it up for integration)?
>
>
> I have some previous open source contribution experience. I have authored the 
> IVTC and Phosphor deinterlacers for VLC, and a modular postprocessing filter 
> framework for the Panda3D game engine (posted at the forums on panda3d.org, 
> projected for version 1.10).
>
> Currently the polynomial solver is available in a git repository hosted by 
> our university, the relevant files being:
>
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/polysolve2.pyx
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/polysolve.py
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/polysolve2example.py
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/setup.py
> https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworods/compile.sh
>
> I'm working on a research grant, which allows me to hold the copyright; 
> license is BSD.
>
> Polysolve2 is the fast Cython+OpenMP version, while polysolve is an earlier 
> (slower) experimental version using only NumPy array operations. The example 
> module contains a usage example, and setup (in the standard manner) instructs 
> distutils and Cython as to how to build polysolve2 (including setting the 
> relevant math flags and enabling OpenMP).
>
> (The setup script includes also some stuff specific to my solver for the 
> Ziegler problem; basically, the "tworods" module can be ignored. I apologize 
> for the inconvenience.)
>
>
> Thoughts?
>
>
> Best regards,
>
> -J
>
> -
> Juha Jeronen, Ph.D.
> juha.jero...@jyu.fi
> University of Jyväskylä
> Department of Mathematical Information Technology
> -
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

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

[Charles R Harris]

On Tue, Sep 29, 2015 at 6:48 PM, Chris Barker - NOAA Federal <
chris.bar...@noaa.gov> wrote:

> This sounds pretty cool -- and I've had a use case. So it would be
> nice to get into Numpy.
>
> But: I doubt we'd want OpenMP dependence in Numpy itself.
>
> But maybe a pure Cython non-MP version?
>
> Are we putting Cuthon in Numpy now? I lost track.
>

Yes, but we need to be careful of Numeric Recipes.



Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Charles R Harris]

On Tue, Sep 29, 2015 at 6:58 PM, Chris Barker - NOAA Federal <
chris.bar...@noaa.gov> wrote:

>
> One of the usecases that has sneaked in during the last few numpy versions
> is that object arrays contain numerical arrays where the shapes don't add
> up to a rectangular array.
>
>
> I think that's the wrong way to dve that problem -- we really should have
> a "proper" ragged array implementation. But is is the easiest way at this
> point.
>
> For this, and other use-cases, special casing Numpy arrays stored in
> object arrays does make sense:
>
> "If this is s a Numpy array, pass the operation through."
>

Because we now (development) use rich compare, the result looks like

In [1]: a = ones(3)

In [2]: b = array([a, -a], object)

In [3]: b
Out[3]:
array([[1.0, 1.0, 1.0],
   [-1.0, -1.0, -1.0]], dtype=object)

In [4]: sign(b)
Out[4]:
array([[1L, 1L, 1L],
   [-1L, -1L, -1L]], dtype=object)

The function returns long integers in order to not special case Python 3.
Hmm, wonder if we might want to change that.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Stephan Hoyer]

On Tue, Sep 29, 2015 at 8:13 AM, Charles R Harris  wrote:

> Due to a recent commit, Numpy master now raises an error when applying the
> sign function to an object array containing NaN. Other options may be
> preferable, returning NaN for instance, so I would like to open the topic
> for discussion on the list.
>

We discussed this last month on the list and on GitHub:
https://mail.scipy.org/pipermail/numpy-discussion/2015-August/073503.html
https://github.com/numpy/numpy/issues/6265
https://github.com/numpy/numpy/pull/6269/files

The discussion was focused on what to do in the generic fallback case. Now
that I think about this more, I think it makes sense to explicitly check
for NaN in the unorderable case, and return NaN is the input is NaN. I
would not return NaN in general from unorderable objects, though -- in
general we should raise an error.

It sounds like Allan has already fixed this in his PR, but it also would
not be hard to add that logic to the existing code. Is this code in the
NumPy 1.10?

Stephan
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Charles R Harris]

On Tue, Sep 29, 2015 at 11:59 AM, Stephan Hoyer  wrote:

> On Tue, Sep 29, 2015 at 8:13 AM, Charles R Harris <
> charlesr.har...@gmail.com> wrote:
>
>> Due to a recent commit, Numpy master now raises an error when applying
>> the sign function to an object array containing NaN. Other options may be
>> preferable, returning NaN for instance, so I would like to open the topic
>> for discussion on the list.
>>
>
> We discussed this last month on the list and on GitHub:
> https://mail.scipy.org/pipermail/numpy-discussion/2015-August/073503.html
> https://github.com/numpy/numpy/issues/6265
> https://github.com/numpy/numpy/pull/6269/files
>
> The discussion was focused on what to do in the generic fallback case. Now
> that I think about this more, I think it makes sense to explicitly check
> for NaN in the unorderable case, and return NaN is the input is NaN. I
> would not return NaN in general from unorderable objects, though -- in
> general we should raise an error.
>
> It sounds like Allan has already fixed this in his PR, but it also would
> not be hard to add that logic to the existing code. Is this code in the
> NumPy 1.10?
>

No. NumPy 1.10 also has differing behavior between python 2 and python 3.
The reason I raise the question now is that current master has replace use
of PyObject_Compare by PyObject_RichCompare for both python 2 and 3. I
would be easy to extend it. I'm less sure of Allan's work, on a quick look
it seems more complicated.

charris@fc [~]$ python3
Python 3.4.2 (default, Jul  9 2015, 17:24:30)
[GCC 5.1.1 20150618 (Red Hat 5.1.1-4)] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> np.sign(np.array([float('nan')]*3, np.object))
array([None, None, None], dtype=object)
>>>
charris@fc [~]$ python2
Python 2.7.10 (default, Jul  5 2015, 14:15:43)
[GCC 5.1.1 20150618 (Red Hat 5.1.1-4)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> np.sign(np.array([float('nan')]*3, np.object))
array([-1, -1, -1], dtype=object)

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

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

[Juha Jeronen]

Hi,

On 29.09.2015 19:16, Matti Picus wrote:

The algorithms for cubics and quadratics come from Numerical Recipes
(3rd ed.), and the quartic problem is internally reduced to a cubic and
two quadratics, using well-known standard tricks.

Nice, wll documented code.
Just to be sure you are on safe ground, you took the *algorithms* but no *code* 
from Numerical Recipes, right? They are pretty clear about their licensing, see
http://www.nr.com/licenses/redistribute.html
This bit from the site might be relevant:
"We never give permission for Numerical Recipes source code to be posted on any 
public server, or distributed with any freeware or shareware package."
Matti


Yes, the code is original (written from scratch), only the algorithms 
come from the book.



 -J

-
Juha Jeronen, Ph.D.
juha.jero...@jyu.fi
University of Jyväskylä
Department of Mathematical Information Technology
-

___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Nathaniel Smith]

On Sep 29, 2015 8:25 AM, "Anne Archibald"  wrote:
>
> IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point arrays.
Why should it be different for object arrays?

The argument for doing it this way would be that arbitrary python objects
don't have a sign, and the natural way to implement something like
np.sign's semantics using only the "object" interface is

if obj < 0:
return -1
elif obj > 0:
return 1
elif obj == 0:
return 0
else:
raise

In general I'm not a big fan of trying to do all kinds of guessing about
how to handle random objects in object arrays, the kind that ends up with a
big chain of type checks and fallback behaviors. Pretty soon we find
ourselves trying to extend the language with our own generic dispatch
system for arbitrary python types, just for object arrays. (The current
hack where for object arrays np.log will try calling obj.log() is
particularly horrible. There is no rule in python that "log" is a reserved
method name for "logarithm" on arbitrary objects. Ditto for the other
ufuncs that implement this hack.)

Plus we hope that many use cases for object arrays will soon be supplanted
by better dtype support, so now may not be the best time to invest heavily
in making object arrays complicated and powerful.

OTOH sometimes practicality beats purity, and at least object arrays are
already kinda cordoned off from the rest of the system, so I don't feel as
strongly as if we were talking about core functionality.

...is there a compelling reason to even support np.sign on object arrays?
This seems pretty far into the weeds, and that tends to lead to poor
intuition and decision making.

-n
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[Charles R Harris]

On Tue, Sep 29, 2015 at 12:16 PM, Nathaniel Smith  wrote:

> On Sep 29, 2015 8:25 AM, "Anne Archibald"  wrote:
> >
> > IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point arrays.
> Why should it be different for object arrays?
>
> The argument for doing it this way would be that arbitrary python objects
> don't have a sign, and the natural way to implement something like
> np.sign's semantics using only the "object" interface is
>
> if obj < 0:
> return -1
> elif obj > 0:
> return 1
> elif obj == 0:
> return 0
> else:
> raise
>

That is what current master does,  using PyObject_RichCompareBool for the
comparisons.

Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion

Re: [Numpy-discussion] Sign of NaN

[josef.pktd]

On Tue, Sep 29, 2015 at 2:16 PM, Nathaniel Smith  wrote:

> On Sep 29, 2015 8:25 AM, "Anne Archibald"  wrote:
> >
> > IEEE 754 has signum(NaN)->NaN. So does np.sign on floating-point arrays.
> Why should it be different for object arrays?
>
> The argument for doing it this way would be that arbitrary python objects
> don't have a sign, and the natural way to implement something like
> np.sign's semantics using only the "object" interface is
>
> if obj < 0:
> return -1
> elif obj > 0:
> return 1
> elif obj == 0:
> return 0
> else:
> raise
>
> In general I'm not a big fan of trying to do all kinds of guessing about
> how to handle random objects in object arrays, the kind that ends up with a
> big chain of type checks and fallback behaviors. Pretty soon we find
> ourselves trying to extend the language with our own generic dispatch
> system for arbitrary python types, just for object arrays. (The current
> hack where for object arrays np.log will try calling obj.log() is
> particularly horrible. There is no rule in python that "log" is a reserved
> method name for "logarithm" on arbitrary objects. Ditto for the other
> ufuncs that implement this hack.)
>
> Plus we hope that many use cases for object arrays will soon be supplanted
> by better dtype support, so now may not be the best time to invest heavily
> in making object arrays complicated and powerful.
>
> OTOH sometimes practicality beats purity, and at least object arrays are
> already kinda cordoned off from the rest of the system, so I don't feel as
> strongly as if we were talking about core functionality.
>
> ...is there a compelling reason to even support np.sign on object arrays?
> This seems pretty far into the weeds, and that tends to lead to poor
> intuition and decision making.
>

One of the usecases that has sneaked in during the last few numpy versions
is that object arrays contain numerical arrays where the shapes don't add
up to a rectangular array.
In those cases being able to apply numerical operations might be useful.

But I'm +0 since I don't work with object arrays.

Josef



> -n
>
> ___
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion