Opened 4 years ago
Last modified 3 weeks ago
#19033 needs_info enhancement
RIF/CIF: search and subintervals by a bisectionalgorithm
Reported by:  dkrenn  Owned by:  

Priority:  major  Milestone:  sage9.1 
Component:  numerical  Keywords:  
Cc:  cheuberg  Merged in:  
Authors:  Daniel Krenn  Reviewers:  
Report Upstream:  N/A  Work issues:  
Branch:  u/dkrenn/rif/bisect (Commits)  Commit:  4f4b2c9fbf56c032e5d3c5cfc31826deb11377f8 
Dependencies:  Stopgaps: 
Description
This ticket implements a bisection algorithm used with interval field elements.
Change History (53)
comment:1 Changed 4 years ago by
 Branch set to u/dkrenn/rif/bisect
comment:2 Changed 4 years ago by
 Commit set to 788fb330983d2ab71f23b264753d2c810aa6d142
 Status changed from new to needs_review
comment:3 followup: ↓ 4 Changed 4 years ago by
Hello Daniel,
This is cool!
I guess that the line
verbose('iteration %s with results in %s of %s cells' % (iteration, len(result), len(open)), level=2)
is quite bad for performances. The string formatting is done whatever the verbose function is doing. Did you do some profiling?
Just for curiosity, do you know how good is fast_callable
with interval fields? It looks pretty bad on the following example
sage: f(x) = log(exp(x*sin(x)) + exp(x*cos(x))) sage: F = fast_callable(f, domain=RIF) sage: r = RIF(0.1,0.2) sage: %timeit F(r) 10000 loops, best of 3: 54.6 µs per loop sage: %timeit ((r*r.sin()).exp() + (r*r.cos()).exp()).log() 10000 loops, best of 3: 32.8 µs per loop
compared to reals
sage: F = fast_callable(f, domain=RR) sage: r = 0.1 sage: %timeit F(r) 100000 loops, best of 3: 12.4 µs per loop sage: %timeit ((r*r.sin()).exp() + (r*r.cos()).exp()).log() 10000 loops, best of 3: 16.5 µs per loop
Vincent
comment:4 in reply to: ↑ 3 ; followup: ↓ 5 Changed 4 years ago by
Hi Vincent,
Replying to vdelecroix:
I guess that the line
verbose('iteration %s with results in %s of %s cells' % (iteration, len(result), len(open)), level=2)is quite bad for performances. The string formatting is done whatever the verbose function is doing. > Did you do some profiling?
No, I didn't do any with this line included (I did on my old code, which used if and print für debugging).
Just for curiosity, do you know how good is
fast_callable
with interval fields?
Don't know. I used it with some larger symbolic expressions, and what I remember, I had a significant speedup there.
It looks pretty bad on the following example
sage: f(x) = log(exp(x*sin(x)) + exp(x*cos(x))) sage: F = fast_callable(f, domain=RIF) sage: r = RIF(0.1,0.2) sage: %timeit F(r) 10000 loops, best of 3: 54.6 µs per loop sage: %timeit ((r*r.sin()).exp() + (r*r.cos()).exp()).log() 10000 loops, best of 3: 32.8 µs per loopcompared to reals
sage: F = fast_callable(f, domain=RR) sage: r = 0.1 sage: %timeit F(r) 100000 loops, best of 3: 12.4 µs per loop sage: %timeit ((r*r.sin()).exp() + (r*r.cos()).exp()).log() 10000 loops, best of 3: 16.5 µs per loop
Not so good...
Any suggestions on the bisectcode using fast_callable?
Daniel
comment:5 in reply to: ↑ 4 ; followups: ↓ 7 ↓ 10 ↓ 11 ↓ 16 Changed 4 years ago by
Replying to dkrenn:
Hi Vincent,
Replying to vdelecroix:
I guess that the line
verbose('iteration %s with results in %s of %s cells' % (iteration, len(result), len(open)), level=2)is quite bad for performances. The string formatting is done whatever the verbose function is doing. > Did you do some profiling?
No, I didn't do any with this line included (I did on my old code, which used if and print für debugging).
If you really want these verbose you should use lazy formatting. I do not remember what is the state of the art, but you can have a look at sage.misc.lazy_string
and sage.misc.lazy_format
. But you are really creating a string within the critical loop. I do find it not very reasonable (the other is fine though).
Just for curiosity, do you know how good is
fast_callable
with interval fields?Don't know. I used it with some larger symbolic expressions, and what I remember, I had a significant speedup there.
I am pretty sure that fast_callable
is better than symbolic and you seem to confirm that. However, you can see in my small snippet above that fast_callable
looks worse than using Python. I should try on a larger expression but then you have to translate into an object call oriented expression and I am lazy...
Any suggestions on the bisectcode using fast_callable?
Check what fast_callable
is doing with domain=RIF
and domain=CIF
(quickly looking at ext/fast_callable.pyx
I would say nothing). It can be done independently of this ticket.
Other comments:
 What is your variable
result
for?
 did you try to more typing. In other words
open = [] > cdef list open = []
anditeration = 0 > cdef size_t iteration = 0
, etc. Cython does good job with Python list when it knows that these are lists.
comment:6 Changed 4 years ago by
 Commit changed from 788fb330983d2ab71f23b264753d2c810aa6d142 to a97ddceab8ea9f0c93aa320edefc39d9d176e6ec
Branch pushed to git repo; I updated commit sha1. New commits:
a97ddce  reintroduced flag bisect_on_success

comment:7 in reply to: ↑ 5 Changed 4 years ago by
Replying to vdelecroix:
 What is your variable
result
for?
It collects the already finished cells. I've rewritten the code and introduced a new option (this was already there in a preversion of this code).
comment:8 followup: ↓ 12 Changed 4 years ago by
Bad syntax for INPUT::
, it should be INPUT:
comment:9 Changed 4 years ago by
 Commit changed from a97ddceab8ea9f0c93aa320edefc39d9d176e6ec to ac3428ca3c6cb3b658bf4d5899df3d81be77d13c
comment:10 in reply to: ↑ 5 Changed 4 years ago by
Replying to vdelecroix:
If you really want these verbose you should use lazy formatting. I do not remember what is the state of the art, but you can have a look at
sage.misc.lazy_string
andsage.misc.lazy_format
. But you are really creating a string within the critical loop. I do find it not very reasonable (the other is fine though).
LazyFormat did not bring something on the performance. Now checking verbosity level directly; this is faster now.
comment:11 in reply to: ↑ 5 Changed 4 years ago by
Replying to vdelecroix:
 did you try to more typing. In other words
open = [] > cdef list open = []
anditeration = 0 > cdef size_t iteration = 0
, etc. Cython does good job with Python list when it knows that these are lists.
Done. Is faster now.
comment:12 in reply to: ↑ 8 Changed 4 years ago by
Replying to chapoton:
Bad syntax for
INPUT::
, it should beINPUT:
Corrected. Also fixed some broken links in the files (in existing code).
Again, needs_review :)
comment:13 followup: ↓ 18 Changed 4 years ago by
 Milestone changed from sage6.9 to sage6.10
 Status changed from needs_review to needs_work
Unless I'm missing something, the case f is None
is completely undocumented and untested.
comment:14 followup: ↓ 19 Changed 4 years ago by
Also, I guess the test is meant to satisfy the assumption that if an interval A
passes, any interval containing A
should also pass. This should be documented.
comment:15 followup: ↓ 20 Changed 4 years ago by
Another detail: you can replace
from sage.rings.real_mpfi import bisect return bisect(f, self, test, **kwds)
by
return real_mpfi.bisect(f, self, test, **kwds)
comment:16 in reply to: ↑ 5 ; followup: ↓ 21 Changed 4 years ago by
Replying to vdelecroix:
 did you try to more typing. In other words
open = [] > cdef list open = []
anditeration = 0 > cdef size_t iteration = 0
, etc. Cython does good job with Python list when it knows that these are lists.
I agree, there should be a lot more typing. All variables which are boolean should be typed bint
and all integers a suitable integer type (Py_ssize_t
for lengths).
comment:17 Changed 4 years ago by
 Commit changed from ac3428ca3c6cb3b658bf4d5899df3d81be77d13c to dd7cb139c411ebda5993f1a8cea320dbbf9b49ed
Branch pushed to git repo; I updated commit sha1. New commits:
565c482  Merge tag '6.9' into t/19033/rif/bisect

dfbbccb  Trac #19033, comment 15: simplify import

1381c29  Trac #19033, comment 16: types specified

6ccfcf0  Trac #19033, comment 13: remove case f=None

dd7cb13  Trac #19033. comment 14: add note on monotonicity of parameter test

comment:18 in reply to: ↑ 13 Changed 4 years ago by
Replying to jdemeyer:
Unless I'm missing something, the case
f is None
is completely undocumented and untested.
Is now removed (since not needed).
comment:19 in reply to: ↑ 14 Changed 4 years ago by
Replying to jdemeyer:
Also, I guess the test is meant to satisfy the assumption that if an interval
A
passes, any interval containingA
should also pass. This should be documented.
I've added a note block explaining this.
comment:20 in reply to: ↑ 15 Changed 4 years ago by
Replying to jdemeyer:
Another detail: you can replace
from sage.rings.real_mpfi import bisect return bisect(f, self, test, **kwds)by
return real_mpfi.bisect(f, self, test, **kwds)
Done.
comment:21 in reply to: ↑ 16 Changed 4 years ago by
Replying to jdemeyer:
Replying to vdelecroix:
 did you try to more typing. In other words
open = [] > cdef list open = []
anditeration = 0 > cdef size_t iteration = 0
, etc. Cython does good job with Python list when it knows that these are lists.I agree, there should be a lot more typing. All variables which are boolean should be typed
bint
and all integers a suitable integer type (Py_ssize_t
for lengths).
I've added cdef bint success
; this saves about 5 percent in time in my tested examples. I've tried to add Py_ssize_t
and other bint
...it does not seem that this brought something. Maybe I am using it wrong...
comment:22 Changed 4 years ago by
 Status changed from needs_work to needs_review
I've also merged in 6.9 (I am sorry for this, since it was not needed, but I didn't have the previous version available and did not want to want until it was compiled).
Set it back to needs_review.
comment:23 followup: ↓ 26 Changed 4 years ago by
 Status changed from needs_review to needs_work
one failing doctest
comment:24 Changed 4 years ago by
 Commit changed from dd7cb139c411ebda5993f1a8cea320dbbf9b49ed to 99c3e9dd55fa1a99e9d726b1ad4aea253b0115f8
Branch pushed to git repo; I updated commit sha1. New commits:
99c3e9d  Revert "Trac #19033, comment 15: simplify import" since imports changed on 6.10.beta0

comment:25 Changed 4 years ago by
 Commit changed from 99c3e9dd55fa1a99e9d726b1ad4aea253b0115f8 to faacb5fec0a838d4ab3aa1ce95c0198b491d59da
Branch pushed to git repo; I updated commit sha1. New commits:
faacb5f  Merge tag '6.10.beta0' into t/19033/rif/bisect

comment:26 in reply to: ↑ 23 Changed 4 years ago by
 Status changed from needs_work to needs_review
Replying to chapoton:
one failing doctest
Was not there on 6.9.
Merged 6.10.beta0 and reverted the changes from comment 15 of this ticket (change of imports). Back to needs_review.
comment:27 Changed 4 years ago by
 Commit changed from faacb5fec0a838d4ab3aa1ce95c0198b491d59da to 0dfca39ab7ffe8693e9757d09c782fd938b0d3f4
Branch pushed to git repo; I updated commit sha1. New commits:
0dfca39  Merge tag '7.1.beta3' into t/19033/rif/bisect

comment:28 Changed 4 years ago by
Merged 7.1.beta3.
comment:29 followup: ↓ 31 Changed 4 years ago by
 Status changed from needs_review to needs_work
does not apply
comment:30 Changed 4 years ago by
 Commit changed from 0dfca39ab7ffe8693e9757d09c782fd938b0d3f4 to d378b1f1b591dde7af4d73b2ea3506efbdc1421b
Branch pushed to git repo; I updated commit sha1. New commits:
d378b1f  Merge tag '7.1.beta4' into t/19033/rif/bisect

comment:31 in reply to: ↑ 29 Changed 4 years ago by
 Status changed from needs_work to needs_review
comment:32 Changed 2 years ago by
 Milestone changed from sage6.10 to sage8.2
 Status changed from needs_review to needs_work
does not apply
comment:33 followup: ↓ 35 Changed 2 years ago by
Note that arb already offers a bissection/newton scheme for root finding of analytic functions, see arb_calc.
comment:34 Changed 10 months ago by
 Commit changed from d378b1f1b591dde7af4d73b2ea3506efbdc1421b to 4f4b2c9fbf56c032e5d3c5cfc31826deb11377f8
Branch pushed to git repo; I updated commit sha1. New commits:
4f4b2c9  Merge tag '8.6' into u/dkrenn/rif/bisect

comment:35 in reply to: ↑ 33 Changed 10 months ago by
Replying to vdelecroix:
Note that arb already offers a bissection/newton scheme for root finding of analytic functions, see arb_calc.
Good to know, thanks.
comment:36 Changed 10 months ago by
 Status changed from needs_work to needs_review
Merged in 8.6, so back to needs review again.
When I see this correctly, all comments given on this ticket so far have been incorporated or answered, but noone gave a positive review three years ago.....I guess, however, we are close to a positive review after all. Can someone have a look again, please.
comment:37 Changed 10 months ago by
 Milestone changed from sage8.2 to sage8.8
comment:38 followup: ↓ 43 Changed 10 months ago by
15 is not addressed. real_mpfi
is already globally cimported in complex_interval.pyx
.
comment:39 followup: ↓ 44 Changed 10 months ago by
I don't think that "bisection on success"/"bisection on failure" is enough. In many instances you have better than a boolean test. Let me consider the situation where you want to isolate all roots of the function f
(e.g. a squarefree polynomial). Interval arithmetic allows you to:
 prove that a cell contains a single zero (if
f(left)
andf(right)
have different signs andf'(cell)
does not vanish)  prove that a cell does not contain a zero (if
f(cell)
does not contain zero)
In this situation, you want to perform bisection if the two above certificate failed. It is a True/False/Unknown
situation.
I think that the bisection should be general enough to handle this search because it does provide a proof of what you are looking for. The return value would be a pairs of lists: the one you are interested in and the pieces for which bisection failed to determinate their status.
Also, isolating roots in an example such as
sage: def contains_zero(fct, cell): ....: return fct(cell).contains_zero() sage: result = bisect(lambda x: x^34*x+2, RIF(10, 10), contains_zero)
would be faster since you would not refine the already valid cells.
comment:40 Changed 10 months ago by
(for all this, you should have a look at arb_calc_isolate_roots
that was already mentioned)
comment:41 followup: ↓ 45 Changed 10 months ago by
The following is a big waste.
if join_neighboring_cells: verbose('joining neighboring cells...', level=1) cells = result result = [] while len(cells) > 0: joined = cells.pop(0) k = 0 while k < len(cells): try: joined.intersection(cells[k]) except ValueError: k += 1 else: joined = joined.union(cells.pop(k)) result.append(joined)
The intervals should be stored in order (a binary tree seems the most appropriate).
comment:42 Changed 10 months ago by
 Status changed from needs_review to needs_work
comment:43 in reply to: ↑ 38 Changed 10 months ago by
Replying to vdelecroix:
15 is not addressed.
real_mpfi
is already globally cimported incomplex_interval.pyx
.
Doesn't work:
return real_mpfi.bisect(f, self, test, **kwds) ^  sage/rings/complex_interval.pyx:2201:24: cimported module has no attribute 'bisect'
comment:44 in reply to: ↑ 39 ; followup: ↓ 48 Changed 10 months ago by
Replying to vdelecroix:
I don't think that "bisection on success"/"bisection on failure" is enough. In many instances you have better than a boolean test. Let me consider the situation where you want to isolate all roots of the function
f
(e.g. a squarefree polynomial). Interval arithmetic allows you to:
 prove that a cell contains a single zero (if
f(left)
andf(right)
have different signs andf'(cell)
does not vanish) prove that a cell does not contain a zero (if
f(cell)
does not contain zero)In this situation, you want to perform bisection if the two above certificate failed. It is a
True/False/Unknown
situation.I think that the bisection should be general enough to handle this search because it does provide a proof of what you are looking for. The return value would be a pairs of lists: the one you are interested in and the pieces for which bisection failed to determinate their status.
I agree that at some point this would be an extension to do. However, I think that it should not be now and that the current implementation has its advantages:
 The test at the moment is compatible and in the sense as comparisons etc. work in RIF/CIF:
sage: RIF(1,3) < RIF(2,4) False
returnsTrue/False
and not unknown. I think there were a lot of discussions going on about this during Py3 comparison discussions. Apparently this is still so in SageMath, so this should be done consistently in all of the module.
 It is easy to use in the sense what the input testfunction should be/return and what the result of
bisect
will be. Therefore this function is suitable for the end user; one main goal of this function.
 The proposed change has a different result; see below.
Also, isolating roots in an example such as
sage: def contains_zero(fct, cell): ....: return fct(cell).contains_zero() sage: result = bisect(lambda x: x^34*x+2, RIF(10, 10), contains_zero)would be faster since you would not refine the already valid cells.
Faster yes, but gives a different result. The current implementation results in small intervals that satisfy the condition, whereas your proposed solution would find some possible large intervals. One could, of course, do this by the proposed functionality, but rather complicated by searching for the inverse and then inverting the result again, which speaks against the simplicity in usage.
comment:45 in reply to: ↑ 41 ; followup: ↓ 47 Changed 10 months ago by
Replying to vdelecroix:
The following is a big waste.
if join_neighboring_cells: verbose('joining neighboring cells...', level=1) cells = result result = [] while len(cells) > 0: joined = cells.pop(0) k = 0 while k < len(cells): try: joined.intersection(cells[k]) except ValueError: k += 1 else: joined = joined.union(cells.pop(k)) result.append(joined)The intervals should be stored in order (a binary tree seems the most appropriate).
This can only be simplified if we talk about real intervals, but not if we talk about complex intervals, or do I miss something here?
comment:46 Changed 10 months ago by
 Status changed from needs_work to needs_info
comment:47 in reply to: ↑ 45 ; followup: ↓ 50 Changed 10 months ago by
Replying to dkrenn:
Replying to vdelecroix:
The following is a big waste.
if join_neighboring_cells: verbose('joining neighboring cells...', level=1) cells = result result = [] while len(cells) > 0: joined = cells.pop(0) k = 0 while k < len(cells): try: joined.intersection(cells[k]) except ValueError: k += 1 else: joined = joined.union(cells.pop(k)) result.append(joined)The intervals should be stored in order (a binary tree seems the most appropriate).
This can only be simplified if we talk about real intervals, but not if we talk about complex intervals, or do I miss something here?
Indeed. For complex intervals you want a quadtree not a binary tree.
comment:48 in reply to: ↑ 44 ; followups: ↓ 49 ↓ 51 Changed 10 months ago by
Replying to dkrenn:
Replying to vdelecroix:
I don't think that "bisection on success"/"bisection on failure" is enough. In many instances you have better than a boolean test. Let me consider the situation where you want to isolate all roots of the function
f
(e.g. a squarefree polynomial). Interval arithmetic allows you to:
 prove that a cell contains a single zero (if
f(left)
andf(right)
have different signs andf'(cell)
does not vanish) prove that a cell does not contain a zero (if
f(cell)
does not contain zero)In this situation, you want to perform bisection if the two above certificate failed. It is a
True/False/Unknown
situation.I think that the bisection should be general enough to handle this search because it does provide a proof of what you are looking for. The return value would be a pairs of lists: the one you are interested in and the pieces for which bisection failed to determinate their status.
I agree that at some point this would be an extension to do.
Good. Then it shoud be clear that all the functionalities will be forward compatible.
However, I think that it should not be now and that the current implementation has its advantages:
 The test at the moment is compatible and in the sense as comparisons etc. work in RIF/CIF:
sage: RIF(1,3) < RIF(2,4) FalsereturnsTrue/False
and not unknown. I think there were a lot of discussions going on about this during Py3 comparison discussions. Apparently this is still so in SageMath, so this should be done consistently in all of the module.
The workaround is very simple
sage: def is_really_less_than(a,b): ....: if a < b: return True ....: elif b <= a: return False ....: else: return Unknown sage: is_really_less_than(RIF(1,2), RIF(3,4)) True sage: is_really_less_than(RIF(1,3), RIF(2,4)) Unknown
 It is easy to use in the sense what the input testfunction should be/return and what the result of
bisect
will be. Therefore this function is suitable for the end user; one main goal of this function.
 The proposed change has a different result; see below.
Also, isolating roots in an example such as
sage: def contains_zero(fct, cell): ....: return fct(cell).contains_zero() sage: result = bisect(lambda x: x^34*x+2, RIF(10, 10), contains_zero)would be faster since you would not refine the already valid cells.
Faster yes, but gives a different result. The current implementation results in small intervals that satisfy the condition, whereas your proposed solution would find some possible large intervals. One could, of course, do this by the proposed functionality, but rather complicated by searching for the inverse and then inverting the result again, which speaks against the simplicity in usage.
When roots are isolated you typically switch to something else. If you know that f has a unique root in [0,1] you would not implement the bisection the way this function proceed. Moreover, if you want really small intervals, Newton method is better suited. So the fact that the code in the current branch does more than the minimum amount of work is a downside to me.
comment:49 in reply to: ↑ 48 Changed 10 months ago by
Replying to vdelecroix:
Replying to dkrenn:
Faster yes, but gives a different result. The current implementation results in small intervals that satisfy the condition, whereas your proposed solution would find some possible large intervals. One could, of course, do this by the proposed functionality, but rather complicated by searching for the inverse and then inverting the result again, which speaks against the simplicity in usage.
When roots are isolated you typically switch to something else. If you know that f has a unique root in [0,1] you would not implement the bisection the way this function proceed. Moreover, if you want really small intervals, Newton method is better suited. So the fact that the code in the current branch does more than the minimum amount of work is a downside to me.
Newton can only be done if the function is differentiable, which is not assumed by the current implementation.
comment:50 in reply to: ↑ 47 Changed 10 months ago by
Replying to vdelecroix:
Replying to dkrenn:
Replying to vdelecroix:
The following is a big waste.
if join_neighboring_cells: verbose('joining neighboring cells...', level=1) cells = result result = [] while len(cells) > 0: joined = cells.pop(0) k = 0 while k < len(cells): try: joined.intersection(cells[k]) except ValueError: k += 1 else: joined = joined.union(cells.pop(k)) result.append(joined)The intervals should be stored in order (a binary tree seems the most appropriate).
This can only be simplified if we talk about real intervals, but not if we talk about complex intervals, or do I miss something here?
Indeed. For complex intervals you want a quadtree not a binary tree.
Followup ticket #27595.
comment:51 in reply to: ↑ 48 Changed 10 months ago by
Replying to vdelecroix:
Replying to dkrenn:
Replying to vdelecroix:
I don't think that "bisection on success"/"bisection on failure" is enough. In many instances you have better than a boolean test. Let me consider the situation where you want to isolate all roots of the function
f
(e.g. a squarefree polynomial). Interval arithmetic allows you to:
 prove that a cell contains a single zero (if
f(left)
andf(right)
have different signs andf'(cell)
does not vanish) prove that a cell does not contain a zero (if
f(cell)
does not contain zero)In this situation, you want to perform bisection if the two above certificate failed. It is a
True/False/Unknown
situation.I think that the bisection should be general enough to handle this search because it does provide a proof of what you are looking for. The return value would be a pairs of lists: the one you are interested in and the pieces for which bisection failed to determinate their status.
I agree that at some point this would be an extension to do.
Good. Then it shoud be clear that all the functionalities will be forward compatible.
After all, this might be difficult, as already the output and the desired functionality seems to be different in our two cases. Do you have something particular in mind here?
comment:52 Changed 7 months ago by
 Milestone changed from sage8.8 to sage8.9
Tickets still needing working or clarification should be moved to the next release milestone at the soonest (please feel free to revert if you think the ticket is close to being resolved).
comment:53 Changed 3 weeks ago by
 Milestone changed from sage8.9 to sage9.1
Ticket retargeted after milestone closed
New commits:
write bisect function
write RIF/CIF.find_subintervals
small improvement in docstring