milestones for natural continuation #6
Conversation
Codecov Report
@@ Coverage Diff @@
## master #6 +/- ##
==========================================
- Coverage 87.93% 84.32% -3.61%
==========================================
Files 5 5
Lines 174 185 +11
==========================================
+ Hits 153 156 +3
- Misses 21 29 +8
Continue to review full report at Codecov.
|
|
The idea here is that one could check
|
|
I'll check this out next week. |
|
Thoughts? I've been using this a bit myself and am happy with the interface and behaviour. |
pacopy/natural.py
Outdated
| @@ -76,7 +82,10 @@ def natural( | |||
| ) | |||
|
|
|||
| # Predictor | |||
| lmbda += lambda_stepsize | |||
| if milestones is None: | |||
There was a problem hiding this comment.
Bit shorter:
lmbda += lambda_stepsize
if milestones:
lmbda = min(lmbda, milestone)
| """ | ||
| lmbda = lambda0 | ||
| if milestones is not None: |
There was a problem hiding this comment.
milestones = [] if milestones is None else milestones
There was a problem hiding this comment.
Ah. The motivation for wrapping milestones in iter is that a user will probably pass either a list or a numpy.ndarray, as in
and one can't call next on either of those.
Lists and NumPy ndarrarys are Iterable
isinstance(np.array([8e2]), typing.Iterable)but not iterators.
There was a problem hiding this comment.
There's a catch here in dispensing with the new name lambdas for the processed argument; see below.
There was a problem hiding this comment.
Note that iter is idempotent and therefore harmless in the case where a user does actually pass an iterator.
A reasonable use-case there is to generate an infinite sequence of milestones, e.g. passing milestones=itertools.count(), presumably terminating with max_steps #8 or by raising an Exception in the callback.
| @@ -109,5 +118,10 @@ def natural( | |||
|
|
|||
| callback(k, lmbda, u) | |||
| k += 1 | |||
| if milestones is not None and lmbda == milestone: | |||
There was a problem hiding this comment.
The float comparison is a bit iffy.
There was a problem hiding this comment.
Yes, interesting. Usually one avoids testing two floats for equality; here however I think it's O. K.
There are two dangers:
- false negative, missing the milestone; and
- false positive, matching something that isn't the milestone.
I don't think the false positive matters. If lmbda is close enough to the milestone to pass an equality test, it's good enough to be considered as it. I think this is fairly unlikely anyway; I have not been able to conceive a working example.
The false negative one is trickier. I think it's O. K., this won't happen because it's not really two floats, but, in the case that the milestone is reached, two (shallow) copies of the same float, being one of the specified milestones. The left term lmbda is assigned here (in the case that the milestone is reached) from the second argument of min.
So the question is whether min changes its second argument before assigning it to the left-hand side. It doesn't appear to; i.e.
milestone = 8e2
lmbda = min(1e3, milestone)
assert lmbda == milestoneIn fact, one even has object identity
assert lmbda is milestonewhile
assert lmbda is not 8e2Perhaps using is here would be clearer?
I don't know. What do you think? Is it safe enough? (Could the behaviour of min differ?) Add a reassuring comment?
Or is there another way to implement this that avoids the issue?
I don't think the usual expedient of a tolerance as in numpy.allclose is appropriate as I can see that producing false negatives which would be messy and annoying. Especially as the other thing is that in the envisaged use-case, the callback will always have another such test; e.g.
There was a problem hiding this comment.
Thinking about this some more. (Thank you for raising it. I intend to make much use of this function and feature, so don't mind spending time now getting it right. I appreciate the rigour.)
Possible 'worst case scenario' (?). Say in
that lmbda is numerically close to milestone but is judged by min to be the lesser. Then, say we successfully solve at this λ and get to
but then fail the second test, the floating-point equality. What will happen then is that we'll go through the next iteration of the while-loop with
and this time lmbda will go well past milestone and so the latter will win the min and we'll be trying to solve for lmbda = milestone that was very close to the last one. This will be solved easily by pacopy.newton because the previous solution at the close lmbda will be an excellent initial guess and we'll be on our way again.
Side-effect: lambda_stepsize might be increased because newton_steps was small. Hmm. Is it worth catching or correcting for this? I don't think that we'd want to reset lambda_stepsize when we truncate a step to land on a milestone as the old lambda_stepsize is a better estimate of the right one thereafter.
There was a problem hiding this comment.
The adjustment of lambda_stepsize is moved into the else clause in ec1b7fc ; i.e. it is not adjusted when the actual step was shortened to meet a milestone. That solves that last side-effect.
|
Oops, hang on, the last commit 457c5ec is wrong. In conflating the argument |
|
Default no-milestones case fixed in d11f9f2. |
|
This is ready for re-review. |
|
Alright, let's merge this. We might tune this or that in the future, but the fact that we have a test gives me some confidence that we won't mess up anyone's applications. |
Implemented using
collections.deque, though the parameter can be passed in as anyIterable, e.g.numpy.ndarray(withndim==1), as intest/test_bratu1d.py, modified to demonstrate.As discussed in #5, if milestones have been specified, the natural parameter continuation terminates after the last has been reached.
The overall effect can be seen in the (λ, || u ||²) plot, with an × on each vertical grid line (Δ λ = 0.5) during the first, natural, pass.
Some other minor changes were made to the example to make it easier to run; e.g.
newton_tolinnaturalas ineuler_newton(for same reason)raise RangeExceptionto terminate subsequenteuler_newtoncontinuation after the region of interest has been leftRangeException, so that the plot can be inspected (previously there was plenty of time for this)