BUG: change poly() real output checking to test if all imaginary parts are zero #4073

endolith · 2013-11-24T21:45:28Z

Some arrays should produce real output but don't:

In [144]: poly([1j, -1j])
Out[144]: array([ 1.,  0.,  1.])

In [145]: poly([1j, -1j, 2j, -2j])
Out[145]: array([ 1.+0.j,  0.+0.j,  5.+0.j,  0.+0.j,  4.+0.j])

In [146]: _.imag == 0
Out[146]: array([ True,  True,  True,  True,  True], dtype=bool)

I think the check was meant to be sort_complex(conjugate(... rather than conjugate(sort_complex(..., which would work, but checking that the imaginary parts are exactly zero is equivalent and simpler anyway.

endolith · 2013-11-24T22:11:55Z

oh wait now it doesn't work for poly([+1.082j, +2.613j, -2.613j, -1.082j]). So maybe make it neg_roots = sort_complex(NX.conjugate( instead. I'll do more tests

charris · 2013-11-25T01:33:58Z

It might be better to have a real keywordarg as the user is likely to have a better sense of whether the result should be real or not. Otherwise, one relies on a certain degree of numerical precision which may not be available depending on where the roots come from.

endolith · 2013-11-25T17:42:56Z

That would be equivalent to real(poly()), right?

Regardless, if the output is exactly real, or input are exact conjugates, it should automatically not use a complex type for output.

endolith · 2013-11-27T05:36:29Z

Ok I changed it back to the conjugate then sort method, and added a bunch of tests. Does this look good?

Is there a reason it's for k in range(len(seq_of_zeros)): instead of for x in seq_of_zeros:?

…ry parts Fixed in numpy/numpy#4073, but we can't depend on a specific numpy version, so copying it here, too.

charris · 2015-01-24T20:47:33Z

numpy/lib/polynomial.py

-        if (len(pos_roots) == len(neg_roots) and
-            NX.alltrue(neg_roots == pos_roots)):
-            a = a.real.copy()
+        pos_roots = NX.compress(roots.imag > 0, roots)


I think you could just check np.all(roots.sort() == roots.conjugate().sort()).

Well, .sort() is in-place sorting and doesn't return anything, so anything.sort() == anythingelse.sort() is always true.

np.all(np.sort(roots) == np.sort(roots.conjugate()))?

np.sort_complex() just calls np.sort() internally but always outputs a complex type? Not sure if it's necessary, but just kept what was in the original code.

You're right about the sort method.

Both the polynomial and sort_complex code date back to 2002 with later refactoring, so I expect not using sort directly is just historical accident.

And since only complex types are checked, the return will be complex anyway.

I assume checking the size first is to speed up rejections for huge arrays?

Though I guess polynomials are always going to be short 1d arrays?

I'd guess that going much over 20 is going to be a problem in general. I've gone up to 60 equally spaced roots in the polynomial package, which does better than this, but the results can be sensitive to root placement.

jaimefrio · 2015-01-25T04:01:21Z

You need to rebase on master, some changes were made to the type handling a while back.

On a separate note regarding the actual object of this PR, wouldn't it make sense to detect the conjugate roots before multiplying? That way you could collapse every two conjugate roots into a strictly real second degree binomial and avoid roundoff errors creeping into the imaginary part of the result.

charris · 2015-01-25T04:26:11Z

@jaimefrio Yes, and for power series that is pretty easy to do, although the convolution here would need modification. For other polynomial basis it gets a bit more involved.

…t conjugates

Change to a deterministic test instead of using rand

endolith · 2015-01-25T17:34:20Z

Ok, I rebased it and fixed an import error. Should I definitely change to np.all(np.sort(roots) == np.sort(roots.conjugate()))?

charris · 2015-01-25T17:48:15Z

@endolith I think it would be cleaner. I don't think speed is of much importance here, but it might even be faster.

endolith · 2015-02-18T01:35:28Z

Ok I changed it.

rgommers · 2015-03-08T22:23:00Z

@charris this seems about ready to merge. ignore the test failure, that's unrelated (couldn't download Cython)

endolith · 2015-03-08T23:11:00Z

Actually instead of that long string of numbers in the test, it could just set the numpy random seed for determinism and go back to

a = rand(100) + 1j* rand(100)
assert_(np.isrealobj(np.poly(np.concatenate((a, np.conjugate(a))))))

also fixed some PEP8 issues

endolith · 2016-02-16T00:41:54Z

Was there anything else that needed to be done on this?

endolith · 2016-06-15T00:02:50Z

Is there anything wrong with this change?

charris · 2016-06-15T14:35:10Z

No, LGTM ;) Thanks @endolith .

Removed code defending against old Numpy: [1] was included in Numpy 1.12.0, and minimum Numpy now is 1.23.5. See [2] for context. [1] numpy/numpy#4073 [2] scipy#3085 (comment)

…test for complex k, real p, z (#22213) Removed code defending against old Numpy: [1] was included in Numpy 1.12.0, and minimum Numpy now is 1.23.5. See [2] for context. [1] numpy/numpy#4073 [2] #3085 (comment)

endolith mentioned this pull request Nov 24, 2013

BUG: buttap() should produce a purely real pole scipy/scipy#3081

Closed

endolith closed this Nov 27, 2013

endolith deleted the patch-2 branch November 27, 2013 05:33

endolith reopened this Nov 27, 2013

endolith added a commit to endolith/scipy that referenced this pull request Jan 19, 2014

BUG: zpkt2f sometimes produces complex output with negligible imagina…

e26799c

…ry parts Fixed in numpy/numpy#4073, but we can't depend on a specific numpy version, so copying it here, too.

charris reviewed Jan 24, 2015
View reviewed changes

endolith added 2 commits January 25, 2015 00:50

BUG: sort after conjugating so that poly() outputs real type for exac…

1fc2687

…t conjugates

TST: Add tests for np.poly()

6a88cb9

endolith force-pushed the patch-2 branch from 50b2133 to 6a88cb9 Compare January 25, 2015 05:53

TST: Fix missing imports, undeterministic test

7d9789f

Change to a deterministic test instead of using rand

Simpler checking of complex conjugates

cb000ea

rgommers added component: numpy.polynomial 03 - Maintenance labels Mar 8, 2015

TST: Set seed for deterministic random test

050f390

also fixed some PEP8 issues

endolith mentioned this pull request Jun 15, 2016

BUG: Use zpk in IIR filter design to improve accuracy scipy/scipy#3085

Merged

charris merged commit e26738a into numpy:master Jun 15, 2016

endolith changed the title ~~BUG: change real output checking to test if all imaginary parts are zero~~ BUG: change poly() real output checking to test if all imaginary parts are zero Apr 22, 2022

roryyorke mentioned this pull request Dec 30, 2024

BUG: zpk2tf works correctly with complex k, real p, z scipy/scipy#22213

Merged

Uh oh!

BUG: change poly() real output checking to test if all imaginary parts are zero #4073

BUG: change poly() real output checking to test if all imaginary parts are zero #4073

Uh oh!

Conversation

endolith commented Nov 24, 2013

Uh oh!

endolith commented Nov 24, 2013

Uh oh!

charris commented Nov 25, 2013

Uh oh!

endolith commented Nov 25, 2013

Uh oh!

endolith commented Nov 27, 2013

Uh oh!

charris Jan 24, 2015

Choose a reason for hiding this comment

Uh oh!

endolith Jan 24, 2015

Choose a reason for hiding this comment

Uh oh!

charris Jan 25, 2015

Choose a reason for hiding this comment

Uh oh!

charris Jan 25, 2015

Choose a reason for hiding this comment

Uh oh!

endolith Jan 25, 2015

Choose a reason for hiding this comment

Uh oh!

endolith Jan 25, 2015

Choose a reason for hiding this comment

Uh oh!

charris Jan 25, 2015

Choose a reason for hiding this comment

Uh oh!

jaimefrio commented Jan 25, 2015

Uh oh!

charris commented Jan 25, 2015

Uh oh!

endolith commented Jan 25, 2015

Uh oh!

charris commented Jan 25, 2015

Uh oh!

endolith commented Feb 18, 2015

Uh oh!

rgommers commented Mar 8, 2015

Uh oh!

endolith commented Mar 8, 2015

Uh oh!

endolith commented Feb 16, 2016

Uh oh!

endolith commented Jun 15, 2016

Uh oh!

charris commented Jun 15, 2016

Uh oh!

Uh oh!