New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix sorting of complex roots #633

Merged
merged 1 commit into from Jun 3, 2018

Conversation

2 participants
@skirpichev
Copy link
Collaborator

skirpichev commented Mar 31, 2018

No description provided.

@skirpichev skirpichev added this to the 0.10 milestone Mar 31, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-rootof-sort branch from b4db1bc to f7d4043 Apr 5, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-rootof-sort branch from a5d10ba to c1ddc28 May 6, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-rootof-sort branch 3 times, most recently from 3ae2239 to 8cee005 May 15, 2018

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented May 20, 2018

@jksuom, this is a proof of concept implementation of the robust mathematical sorting for polynomial roots, asked in sympy/sympy#14293. I hope, it follows your suggestions in the issue thread.

So far, it doesn't break anything in the test suite, except test_RootOf_evalf() (which now randomly fails without commited workarround). test_RootOf_algebraic_domain() (abcent in sympy) - also much more slow now, but caching resultants helps a bit.

Please, review this work if it does make sense for you as a variant for sympy's bugfix.

@skirpichev skirpichev changed the title [poc] Fix sorting of complex roots [wip] Fix sorting of complex roots May 20, 2018

assert u.domain == v.domain
dom = u.domain
if tuple(u.f) not in res_cache:
res_cache[tuple(u.f)] = dmp_resultant(dmp_permute(u.f1, [1, 0], 1, dom),

This comment has been minimized.

@jksuom

jksuom May 20, 2018

Contributor

Where are u.f, u.f1, u.f2 defined?

This comment has been minimized.

@skirpichev

skirpichev May 20, 2018

Collaborator

@jksuom, they are used (or not, as for f) - in the ComplexInterval class for bookeeping. f - original polynomial (AFAIK, not saved in Sympy). f1 and f2 - its real and imaginary parts.

return False
else:
a, b = roots[0][0]
if a == b:

This comment has been minimized.

@skirpichev

skirpichev May 20, 2018

Collaborator

Beware, that part seems to be redundant. But dup_count_complex_roots sometimes fails in case, where for x-component inf==sup.

For example:

In [4]: dup_count_complex_roots([1, 0, 1], ZZ_python, inf=0, sup=0)
Out[4]: 2

In [5]: dup_count_complex_roots([1, 0, 1], ZZ_python, inf=(0, -1), sup=0)
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-5-d0b6fcdb6ef6> in <module>()
----> 1 dup_count_complex_roots([1, 0, 1], ZZ_python, inf=(0, -1), sup=0)

~/src/diofant/diofant/polys/rootisolation.py in dup_count_complex_roots(f, K, inf, sup, exclude)
   1303     Q_L4 = _intervals_to_quadrants(I_L4, f1L4F, f2L4F, t, v, F)
   1304 
-> 1305     T = _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4, exclude=exclude)
   1306 
   1307     return _winding_number(T, F)

~/src/diofant/diofant/polys/rootisolation.py in _traverse_quadrants(Q_L1, Q_L2, Q_L3, Q_L4, exclude)
   1181                 rules.append((_rules_ambiguous[qq], corners[(j, i)]))
   1182             else:  # pragma: no cover
-> 1183                 raise NotImplementedError("3 element rule (corner): " + str(qq))
   1184 
   1185         q1, k = Q[0], 1

NotImplementedError: 3 element rule (corner): ('OO', 'OO', 'A1')

This comment has been minimized.

@skirpichev

skirpichev May 23, 2018

Collaborator

BTW, the answer here should be 2, I believe. We can consider this degenerate case as limit:

In [11]: dup_count_complex_roots([1, 0, 1], ZZ_python, inf=(-1/100000000,-1), sup=1/100000)
Out[11]: 2
In [12]: dup_count_complex_roots([1, 0, 1], ZZ_python, inf=(0,-1), sup=1/100000)
Out[12]: 2
In [13]: dup_count_complex_roots([1, 0, 1], ZZ_python, inf=(-1/100000000,-1), sup=0)
Out[13]: 2

See sympy/sympy#14738.

This comment has been minimized.

@jksuom

jksuom May 23, 2018

Contributor

The corners are particularly tricky. They may lead to degenerate intervals. Maybe they should be extracted and treated separately.

This comment has been minimized.

@skirpichev

skirpichev May 23, 2018

Collaborator

I've added a rule, that removes need in this workaround. Not sure if it's correct, yet.

BTW, it seems that only line segments, which are parallel to imaginary axis - poses problems. Other cases of degenerate isolation rectangles works ok.

res_cache[tuple(v.f)] = dmp_resultant(dmp_permute(v.f1, [1, 0], 1, dom),
dmp_permute(v.f2, [1, 0], 1, dom),
1, dom)
g = dmp_gcd(res_cache[tuple(u.f)], res_cache[tuple(v.f)], 0, dom)

This comment has been minimized.

@jksuom

jksuom May 22, 2018

Contributor

I would have expected u.f and v.f be the same. Otherwise the real parts of u and/or v might not necessarily be roots of g.

This comment has been minimized.

@skirpichev

skirpichev May 22, 2018

Collaborator

Why do you expect so? In general, u.f and v.f are different factors of polynomial. For example, x**2 + 1 and x**2 + 2:

In [11]: p1 = x**2 + 1; p2 = x**2 + 2;

In [12]: var('u v', real=True)
Out[12]: (u, v)

In [13]: r1, i1 = p1.subs({x: u + I*v}).as_real_imag(); r2, i2 = p2.subs({x: u + I*v}).as_real_imag();

In [14]: resultant(r1, i1, v)
Out[14]: 
   4      2
4⋅u  + 4⋅u 

In [15]: resultant(r2, i2, v)
Out[15]: 
   4      2
4⋅u  + 8⋅u 

In [16]: gcd(_14, _15)
Out[16]: 
   2
4⋅u 

This comment has been minimized.

@jksuom

jksuom May 22, 2018

Contributor

I hadn't been thinking about what should be done that case. It could be possible that the two resultants have a common real root but that is not the real part of one (or any) of the roots. Perhaps the intervals should be refined separately for each root until there is only one real root left. Then I think that gcd would be ok.

This comment has been minimized.

@skirpichev

skirpichev May 22, 2018

Collaborator

It could be possible that the two resultants have a common real root but that is not the real part of one (or any) of the roots.

I hope, this case should be covered by return False on L344. (And L341, for a special case, which dup_count_complex_roots() currently can't handle.) In other words, currently we test that gcd has single real root in the given interval (intersection of real bounds from original isolation rectangles) and then (if so) - that same is true for complex roots of every factor (imaginary bounds are taken from original isolation rectangles), i.e. u.f and v.f.

If there is more than one gcd root - match_same_real_part() will return False and refinement process continue.

Perhaps the intervals should be refined separately for each root until there is only one real root left.

Maybe it could be made more efficiently (at present, both intervals will be refined if match_same_real_part() return False). But I still don't understand why the current implementation is wrong here. If gcd has single real roots in the given interval and we tested successfully that both polynomial have single complex roots with real part in that interval - what did I miss?

This comment has been minimized.

@skirpichev

skirpichev May 26, 2018

Collaborator

@jksuom, I did some cleanup for helper is_same_real_part(), which (hopefuly) make things more transparent and clear. Let me know if you are still interested in review.

I think, that your objection "It could be possible that the two resultants have a common real root but that is not the real part of one (or any) of the roots." is not valid. If it's the case - is_same_real_part() return False. If you disagree, could you explain what did I miss?

This comment has been minimized.

@skirpichev

skirpichev May 27, 2018

Collaborator

Could it be possible that the real parts are the same and, in addition, each polynomial has another root with same but slightly different real part that happens to be in the same real interval

How this could be? If there are other pair of roots with slightly different real part - we should continue refinement of u and v (by returning False) and eventually intervals for different multiplets of roots with same real part will be disjoint (as refinement process will shrink intersection of u and v). That's why I do think that we will have len(gcd_roots) == 1 at some point (if u and v belongs to the multiplet of roots with same real part).

to be in the same real interval though the imaginary parts may not be in the (ay, by) intervals of u and v

Why we ever should care about roots, outside of intervals for u and v? is_same_real_part() is a decision procedure: for given isolation rectangles u and v we want to decide if these particular roots have same real part or not (in which case we must continue refinement).

This comment has been minimized.

@jksuom

jksuom May 27, 2018

Contributor

if these particular roots have same real part or not (in which case we must continue refinement).

I see. You do not expect is_same_real_part returning False to mean that the real parts cannot possibly be same. The loop on line 329 will continue refinement. In view of that, I think that this should work correctly.

This comment has been minimized.

@skirpichev

skirpichev May 27, 2018

Collaborator

You do not expect is_same_real_part returning False to mean that the real parts cannot possibly be same.

Yes, False just mean: continue refinement. This process will continue until roots will be either disjoint by real part or gcd will have exactly one root. In the later case we will check if both polys have complex root with real part in the interval of gcd's root. If yes - we have found a pair of roots with same real part and should stop refinement. If not, we will return False and continue refinement again: eventually real parts also will be different.

I think that this should work correctly.

So, no issues with math you have found so far?

PS: Workaround for sympy/sympy#14738 simplifies code a lot, so I would like to include the real fix for that issue as well. But original code for a==b case from 2c706b8 is sufficiently straightforward (substitute x -> a + I*x and search for real solutions), so review of that part is less important.

This comment has been minimized.

@jksuom

jksuom May 27, 2018

Contributor

No math issues. This looks good. The name is_same_real_part was misleading me. Maybe it could be something like confirm_same_real_part.

This comment has been minimized.

@skirpichev

skirpichev May 27, 2018

Collaborator

Maybe it could be something like confirm_same_real_part.

Yes, that sounds better. But I'm thinking about moving this functional from that helper function to the ComplexInterval.is_disjoint() method. Some option could enable all this.

dmp_permute(v.f2, [1, 0], 1, dom),
1, dom)
g = dmp_gcd(res_cache[tuple(u.f)], res_cache[tuple(v.f)], 0, dom)
a, b = max(u.ax, v.ax), min(u.bx, v.bx)

This comment has been minimized.

@jksuom

jksuom May 22, 2018

Contributor

I think that max and min should be interchanged. There could be a single root in the intersection of the intervals but that could be the wrong real part if there are several roots in the union of the intervals.

This comment has been minimized.

@skirpichev

skirpichev May 22, 2018

Collaborator

But comon roots should be in the intersection of isolation rectangles, isn't? If there is any, of course. So, this code refines initial bounds for real part, nothing more. If there are some roots of g outside of this interval (namely, in the union, not in intersection) - following code expected to return False.

This comment has been minimized.

@jksuom

jksuom May 22, 2018

Contributor

If it is known at this point that there is a single root in a rectangle, then I think this is ok.

This comment has been minimized.

@skirpichev

skirpichev May 22, 2018

Collaborator

Yes, that should be known after L327. If we got something different (not a single root) - procedure will return False and refinement will be continued on L351.

@skirpichev skirpichev force-pushed the skirpichev:new-rootof-sort branch from 768ab0e to 988b843 May 23, 2018

@@ -1860,15 +1862,14 @@ def conjugate(self):
self.F1, self.F2, self.f1, self.f2,
self.f, self.domain, conj=not self.conj)

def is_disjoint(self, other):
def is_disjoint(self, other, strict=False):

This comment has been minimized.

@skirpichev

skirpichev May 27, 2018

Collaborator

strict, probably, a wrong name. Maybe re_distinct?

@skirpichev skirpichev force-pushed the skirpichev:new-rootof-sort branch from ae9a1f2 to 8815ba4 May 28, 2018

@skirpichev skirpichev changed the title [wip] Fix sorting of complex roots Fix sorting of complex roots May 28, 2018

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented May 28, 2018

Now all logic for testing real parts of roots moved to the ComplexInterval.is_disjoint() method. It could be enabled by re_distinct boolean option (see docstring). I'm a bit unsure about naming of this flag.

@jksuom, let me know if you would like to review part, related to the sympy/sympy#14738.

@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented May 29, 2018

part, related to the sympy/sympy#14738.

It seems that this is in another PR. Yes, I think I would like to see that. (Though I'm not sure if I can be helpful within reasonable time. The subject is very low level and rather tedious to work with.)

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented May 29, 2018

It seems that this is in another PR.

You mean, this belongs to another PR? Yeah, but fix should be very simple one-liner and I prefer to solve this issue, not introduce ugly workaround. Ok, I'll try to handle this myself.

Any final remarks? re_distinct flag really misnamed, through I don't see better option yet.

@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented May 29, 2018

this belongs to another PR

Yes, though I understood that you might already have that prepared.

is_disjoint seems to return True when either the rectangles really are disjoint (except possibly at boundary) or re_distinct is True and the roots are seen to have the same real part (regardless of whether the rectangles overlap or not). That might be confusing to a user as this is a public method.

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented May 29, 2018

is_disjoint seems to return True when either the rectangles really are disjoint (except possibly at boundary) or re_distinct is True and the roots are seen to have the same real part

is_disjoint() works as before (test that isolation rectangles are disjoint on complex plane), if re_distinct is False (default).

If re_distinct is True - is_disjoint() will return True either if we can confirm that roots have different real parts OR that they have common real part. Otherwise, it will return False. That's documented in the docstring, but option indeed seems to have confusing name. split_real_part maybe?

@skirpichev skirpichev changed the title Fix sorting of complex roots [wip] Fix sorting of complex roots May 30, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-rootof-sort branch from 96f6757 to 13f9f4b Jun 2, 2018

@skirpichev skirpichev force-pushed the skirpichev:new-rootof-sort branch from ca78f92 to 5138161 Jun 3, 2018

@skirpichev skirpichev changed the title [wip] Fix sorting of complex roots Fix sorting of complex roots Jun 3, 2018

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Jun 3, 2018

@jksuom, thank you for your helpful review. I think that finally may be merged.

@jksuom

This comment has been minimized.

Copy link
Contributor

jksuom commented Jun 3, 2018

Thank you, though I am not sure that I was able to help.

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Jun 3, 2018

I am not sure that I was able to help.

If I not misunderstood, you said no math issues were found you so far. That view helps, for sure.

Naming of options (which now is slightly better, I hope) and fixing sympy/sympy#14738 - less important issues.

@skirpichev skirpichev merged commit f7bf793 into diofant:master Jun 3, 2018

3 checks passed

codecov/patch 100% of diff hit (target 97%)
Details
codecov/project Absolute coverage decreased by -<1% but relative coverage increased by +3% compared to f4692df
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

@skirpichev skirpichev deleted the skirpichev:new-rootof-sort branch Jun 3, 2018

@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Jul 20, 2018

FYI: sorting with conjugated pairs of roots corrected in #658.

skirpichev added a commit to skirpichev/diofant that referenced this pull request Dec 7, 2018

polys: optimize RootOf._eval_conjugate()
After introduction of stable sorting of polynomial roots (diofant#633) and
proper detection of imaginary roots (diofant#625) - it's possible to
apply more efficient method of taking conjugated root for
polynomials with integer coefficients.

skirpichev added a commit to skirpichev/diofant that referenced this pull request Dec 7, 2018

polys: optimize RootOf._eval_conjugate()
After introduction of stable sorting of polynomial roots (diofant#633) and
proper detection of imaginary roots (diofant#625) - it's possible to
apply more efficient method of taking conjugated root for
polynomials with integer coefficients.

This commit also add coverage tests for ComplexInterval.

skirpichev added a commit to skirpichev/diofant that referenced this pull request Dec 7, 2018

polys: optimize RootOf._eval_conjugate()
After introduction of stable sorting of polynomial roots (diofant#633) and
proper detection of imaginary roots (diofant#625) - it's possible to
apply more efficient method of taking conjugated root for
polynomials with integer coefficients.

This commit also add coverage tests for ComplexInterval.

skirpichev added a commit to skirpichev/diofant that referenced this pull request Dec 7, 2018

polys: optimize RootOf._eval_conjugate()
After introduction of stable sorting of polynomial roots (diofant#633) and
proper detection of imaginary roots (diofant#625) - it's possible to
apply more efficient method of taking conjugated root for
polynomials with integer coefficients.

This commit also add coverage tests for ComplexInterval.
@skirpichev

This comment has been minimized.

Copy link
Collaborator

skirpichev commented Dec 28, 2018

FYI: #741 fix some performance problems with caching, introduced by correct sorting.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment