Line data Source code
1 : """
2 : trypush!(spline, t, y, dy)
3 :
4 : Push new step `(t, y, dy)` onto `spline`, but check if time `t` is
5 : consistent: Do **not** push if either the same `t` is currently the
6 : last knot of `spline` (ignore silently) or if `t < last(spline)`. The
7 : latter may yield a warning, because there may be a flaw in the
8 :
9 : output (cascade).
10 :
11 : """
12 142830 : function trypush!(s::HermiteSpline.AbstractHermiteSpline, t, y, dy)
13 283657 : if HermiteSpline.ismonotone(s, t)
14 142826 : push!(s, t, y, dy)
15 4 : elseif t != last(s)
16 0 : @warn "Cannot push t=$t after at $(last(s))"
17 : end
18 142830 : nothing
19 : end
20 :
21 : #-----------------------------------------------------------------------------
22 :
23 : """
24 : sink(t, y, dy) # collect data from solver
25 :
26 : Data structure that collects results from [`Solver`](@ref)/
27 :
28 : A *sink* can be used as "end-point" of an [`RK43.AbstractOutput`](@ref).
29 :
30 : See also [`RK43.AbstractOutput`](@ref), [`RK43.TryPush`](@ref)
31 : """
32 : abstract type AbstractSink end
33 :
34 0 : (::AbstractSink)(t, y, dy) = @show((t=t, y=y, dy=dy))
35 :
36 : """
37 : [`RK43.AbstractSink`](@ref) that calls [`RK43.trypush!`](@ref) on provided data
38 : structure (e.g., a `HermiteSpline`).
39 :
40 : See also [`RK43.AbstractSink`](@ref), [`pushinto`](@ref)
41 : """
42 : struct TryPush{S} <: AbstractSink
43 : sink::S
44 2005 : TryPush{S}(sink::S) where {S} = new(sink)
45 : end
46 :
47 142830 : (s::TryPush{S})(t, y, dy) where {S} = RK43.trypush!(s.sink, t, y, dy)
48 :
49 : """
50 : sink = indexedfill(array, Val(2))
51 :
52 : [`RK43.AbstractSink`](@ref) that calls fills an array or array view for
53 : indexed times/knots with `index` starting from `1` for every event.
54 :
55 : The type parameter `P::Int` defines a projection of the event tuple
56 : `(t, y, dy)`, e.g., `P=2` stores positions `y`.
57 :
58 : !!! note
59 : Use with dense [`sample`](@ref)s to fill (part of) an array, e.g.,
60 : for sampling a *flow map*.
61 :
62 : See also [`RK43.AbstractSink`](@ref), [`RK43.indexedfill`](@ref), [`sample`](@ref)
63 : """
64 : mutable struct IndexedFill{A, P} <: AbstractSink
65 : const ary::A
66 : idx::Int
67 2862 : IndexedFill{A, P}(ary::A) where {A<:AbstractArray, P} = new(ary, 1)
68 : end
69 :
70 10000 : _projection(::Val{1}, t, _, _) = t
71 99444 : _projection(::Val{2}, _, y, _) = y
72 10000 : _projection(::Val{3}, _, _, dy) = dy
73 :
74 119444 : function (s::IndexedFill{A, P})(t, y, dy) where {A, P}
75 119444 : s.ary[s.idx] = _projection(Val(P), t, y, dy)
76 119444 : s.idx += 1
77 : end
78 :
79 : """
80 : sink = indexedfill(array, Val(2))
81 :
82 : Create [`IndexedFill`](@ref) sink.
83 :
84 : See also [`RK43.AbstractSink`](@ref), [`RK43.IndexedFill`](@ref), [`sample`](@ref)
85 : """
86 2863 : indexedfill(a::A, proj::Val{N} = Val(2)) where {A<:AbstractArray, N} =
87 : IndexedFill{A, N}(a)
88 :
89 : #-----------------------------------------------------------------------------
90 :
91 : """
92 : out(solver) # process current state of solver
93 :
94 : Process events/results from [`Solver`](@ref).
95 :
96 : 1. Steer `solver`: `out(solver)` may return
97 : - [`RK43.State`](@ref): `Ok` "commits" step, `RK43.Stopped` stops integration,
98 : `OutOfDomain` repeats step with adaptation or bisection.
99 : - `Bool`, where `true` and `false are mapped to `Ok` and `Stopped`,
100 : respectively.
101 : - `nothing`, which is equivalent to `Ok`
102 : 2. Store results if `out` is associated with an [`RK43.AbstractSink`](@ref).
103 :
104 : !!! note
105 : Multiple outputs may be cascaded to a [`Pipeline`](@ref). A
106 : [`RK43.Pipeline`](@ref) is an [`RK43.AbstractOutput`](@ref).
107 :
108 : See also [`RK43.Pipeline`](@ref), [`RK43.AbstractSink`](@ref), [`RK43.Output`](@ref),
109 : [`Callback`](@ref), [`NoOutput`](@ref)
110 : """
111 : abstract type AbstractOutput end
112 :
113 0 : (::AbstractOutput)(::Solver{N, T, S}) where {N, T, S} = Ok
114 :
115 : """
116 : No output: discard all events.
117 :
118 : See also [`RK43.nooutput`](@ref), [`RK43.AbstractOutput`](@ref)
119 : """
120 920 : struct NoOutput <: AbstractOutput
121 : end
122 :
123 : """
124 : No output: discard all events
125 :
126 : See also [`NoOutput`](@ref)
127 : """
128 26 : nooutput() = NoOutput()
129 :
130 : """
131 : Output that is associated with `Snk`, an [`RK43.AbstractSink`](@ref).
132 :
133 : See also [`RK43.AbstractOutput`](@ref), [`RK43.AbstractSink`](@ref),
134 : [`RK43.output`](@ref), [`RK43.Callback`](@ref)
135 : """
136 : struct Output{Snk} <:AbstractOutput
137 : sink::Snk
138 1066 : Output{Snk}(sink::Snk) where {Snk<:AbstractSink} = new(sink)
139 : end
140 :
141 48930 : function (out::Output{Snk})(solver::RK43.Solver{T, N, S}) where {T, N, S, Snk}
142 48930 : out.sink(solver.tnew, solver.ynew, solver.dynew)
143 48930 : RK43.Ok
144 : end
145 :
146 : """
147 : out = output(sink)
148 :
149 : Create an [`RK43.Output`](@ref) associated with `sink`.
150 :
151 : See also [`RK43.Output`](@ref), [`pushinto`](@ref)
152 : """
153 1066 : output(sink::S) where {S<:AbstractSink} = Output{S}(sink)
154 :
155 : """
156 : out = pushinto(spline)
157 :
158 : Create an [`RK43.Output`](@ref) and a [`RK43.TryPush`](@Ref) sink associated with
159 : `spline`.
160 :
161 : See also [`RK43.Output`](@ref), [`RK43.TryPush`](@ref)
162 : """
163 1066 : pushinto(s::S) where {S} = output(TryPush{S}(s))
164 :
165 : # TODO: [alg:upstream] filteredpushinto -- push or replace last if Δh < thresh (don't record excessive bisection)
166 :
167 : """
168 : cb = callback() do solver
169 : # ...
170 : Ok # return state
171 : end
172 :
173 : In contrast to [`RK43.Output`](@ref) `Callback` has no [`RK43.AbstractSink`](@ref).
174 : The called function must return [`RK43.State`](@ref), `Bool` or `nothing` (see
175 : [`RK43.AbstractOutput`](@ref)).
176 :
177 : See also See also [`RK43.AbstractOutput`](@ref), [`RK43.AbstractSink`](@ref),
178 : [`callback`](@ref)
179 : """
180 : struct Callback{F} <: AbstractOutput
181 : callback::F
182 37 : Callback{F}(f::F) where {F} = new(f)
183 : end
184 :
185 473 : function (out::Callback{F})(solver::RK43.Solver{T, N, S}) where {T, N, S, F}
186 726 : out.callback(solver)
187 : end
188 :
189 : """
190 : Create [`Callback`](@ref)
191 : """
192 37 : callback(f::F) where {F} = Callback{F}(f)
193 :
194 229 : _characteristic(s::Bool) = s
195 330 : _characteristic(s::T) where {T<:Real} = (s >= 0)::Bool
196 :
197 6 : function isstepinside(f::F) where {F}
198 6 : callback() do solver
199 133 : (f(solver.ynew) |> _characteristic) ? AcceptStep : OutOfDomain
200 : end
201 : end
202 :
203 5 : function isstepinside(f::F, n::Int) where {F}
204 5 : callback() do solver
205 194 : for t in range(solver.t, solver.tnew, length=n)
206 157 : y, _ = _interpolate(solver, t)
207 157 : (f(y) |> _characteristic) || return OutOfDomain
208 75 : end
209 0 : AcceptStep
210 : end
211 : end
212 :
213 10 : function isstepinside(f::F, ::Val{N}) where {F, N}
214 10 : callback() do solver
215 388 : for t in range(solver.t, solver.tnew, length=N+1)[2:end]
216 269 : y, _ = _interpolate(solver, t)
217 269 : (f(y) |> _characteristic) || return OutOfDomain
218 105 : end
219 0 : AcceptStep
220 : end
221 : end
222 :
223 : # NOTE:
224 : # - Chose not to implement isstepinside via samplestep: avoid passing
225 : # generator even though it's probably zero cost.
226 : # - Special case ::Val{2} makes no sense (need only testing ynew)
227 :
228 : """
229 : chi(y) = norm(y) <= 1 # characteristic function of unit circle
230 : chi(y) = 1 - norm(y) # alternative: >= 0 inside unit circle
231 :
232 : cb = isstepinside(chi)
233 : cb = isstepinside(chi, n)
234 : cb = isstepinside(chi, Val(N))
235 :
236 : Cretae [`callback`](@ref) that checks if the current, not yet
237 : [`commit!`](@ref)ed step is inside domain by *sampling* the
238 : characteristic function of the domain `chi`.
239 :
240 : The characteristic function either returns `Bool` (`true` ⇔ inside) or
241 : a `sign::Real` (`s >= 0` ⇔ inside).
242 :
243 : !!! warning "precondition"
244 : Assumes that the current starting position for the step is inside.
245 : (This is *not* tested.)
246 :
247 : !!! warning
248 : Testing only the end point of the current, tentative step is
249 : sufficient (but *without guarantees*) for **convex** domains
250 : **only**!
251 :
252 : !!! warning
253 : Testing `n` points by sampling the new, tentative trajectory piece
254 : gives **no guarantee** that a boundary intersection is detected!
255 :
256 : !!! note
257 : In order to have a guarantee for non-convex domains, one must
258 : *intersect* the new, tentative trajectory piece with the boundary
259 : curve! This is **not** done by `isstepinside`!
260 :
261 : !!! note
262 : Use [`samplestep`](@ref) or [`cpstep`](@ref) if all samples or
263 : the control polygon are required, e.g., for intersection tests.
264 :
265 :
266 : See also [`callback`](@ref), [`samplestep`](@ref)
267 : """ isstepinside
268 :
269 1 : function samplestep(f::F, ::Val{N}) where {F, N}
270 1 : callback() do solver
271 5 : ts = range(solver.t, solver.tnew, length=N)
272 5 : ys = (first(_interpolate(solver, t)) for t in ts)
273 5 : f(ys)
274 : end
275 : end
276 :
277 2 : function samplestep(f::F, ::Val{2}) where {F}
278 2 : callback() do solver
279 10 : f((solver.y, solver.ynew))
280 : end
281 : end
282 :
283 1 : samplestep(f::F) where {F} = samplestep(f, Val(2))
284 :
285 : """
286 : cb = samplestep(Val(N)) do ys
287 : # check positions ys sampled from tentative piece
288 : return AccepStep # or OutOfDomain
289 : end
290 :
291 : cb = samplestep(f) # default: Val(2)
292 :
293 : Sample trajectory piece from current, tentative step
294 : (with pending [`commit!`](@ref)) at `N` regularly spaced samples
295 : passed as tuple `ys`.
296 :
297 : If `Val(N)` is omitted, `samplestep` passes the two endpoints `(y,
298 : ynew)` of the polynomial piece as a tuple.
299 :
300 : !!! note
301 : Samples are passed as generator or tuple.
302 :
303 : See also [`callback`](@ref), [`isstepinside`](@ref), [`cpstep`](@ref)
304 : """ samplestep
305 :
306 : """
307 : cb = cpstep() do bs
308 : bs :: Tuple
309 : # check control polygon of tentative step
310 : return AccepStep # or OutOfDomain
311 : end
312 :
313 : Provide control polygon defined by Bernstein-Bezier form of tentative
314 : step ("Bezier points") parametrized over `[0, 1]` for testing
315 : `OutOfDomain`: The *convex hull property* can provide a guarantee that
316 : the step is inside (but cannot decide for any step whether it leaves
317 : the domain).
318 :
319 : See also [`callback`](@ref), [`isstepinside`](@ref), [`samplestep`](@ref)
320 : """
321 1 : function cpstep(f::F) where {F}
322 1 : callback() do solver
323 5 : @massert order0(solver.stepper) == 4 "assume cubic piece"
324 5 : (solver.t == solver.tnew) && return RK43.Ok
325 :
326 4 : h = solver.tnew - solver.t
327 : # @show h
328 :
329 4 : b0 = solver.y
330 4 : b1 = b0 + solver.dy*(h/3)
331 4 : b3 = solver.ynew
332 4 : b2 = b3 - solver.dynew*(h/3)
333 :
334 4 : f((b0, b1, b2, b3))
335 : end
336 : end
337 :
338 : """
339 : out = sample(sink, range(0, 1, length=4))
340 : out = sample(sink, [0, 0.25, 0.5, 0.75, 1])
341 : out = sample(array, times[, Val(P)])
342 :
343 : Sample solution at given times during integration. Samples can be specified,
344 : e.g., as range or `Vector`.
345 :
346 : !!! warning "precondition"
347 : The sequence of samples must be *increasing* for forward integration
348 : **or** *decreasing* for backward integration.
349 :
350 : See also [`RK43.AbstractOutput`](@ref), [`RK43.AbstractSink`](@ref),
351 : [`sample`](@ref)
352 : """
353 : mutable struct Sample{Snk, Ts} <: AbstractOutput
354 : const sink::Snk
355 : const times::Ts
356 : idx::Int
357 3800 : Sample{Snk, Ts}(sink::Snk, times::Ts) where {Snk, Ts} =
358 : new(sink, times, 1)
359 : end
360 :
361 147653 : function (out::Sample{Snk, Ts})(solver::RK43.Solver{T, N, S}) where {T, N, S, Snk, Ts}
362 147653 : n = length(out.times)
363 147653 : i = out.idx
364 :
365 147653 : (i <= n) || return Ok
366 :
367 147653 : if isforward(solver)
368 141313 : (solver.tnew < out.times[i]) && return Ok
369 :
370 94275 : while (i <= n) && (out.times[i] < solver.t)
371 0 : i+= 1
372 0 : end
373 : else
374 6340 : (solver.tnew > out.times[i]) && return Ok
375 :
376 3700 : while (i <= n) && (out.times[i] > solver.t)
377 0 : i+= 1
378 0 : end
379 : end
380 :
381 335843 : while (i <= n) && isinside(solver, out.times[i])
382 213344 : out.sink(out.times[i], _interpolate(solver, out.times[i])...)
383 213344 : i += 1
384 213344 : end
385 :
386 97975 : out.idx = i
387 0 : Ok
388 : end
389 :
390 : # TODO: [alg] make sample(...) emit/repeat last point at boundaries
391 :
392 2861 : sample(sink::Snk, times::Ts) where {Snk<:AbstractSink, Ts} =
393 : Sample{Snk, Ts}(sink, times)
394 :
395 939 : sample(s::S, times::Ts) where {S, Ts} =
396 : Sample{TryPush{S}, Ts}(TryPush{S}(s), times)
397 :
398 5522 : function sample(ary::AbstractArray, times::Ts, ::Val{N} = Val(2)) where {Ts, N}
399 5522 : sample(indexedfill(ary, Val(N)), times)
400 : end
401 :
402 100 : sample(ary::AbstractArray, times::Ts, n::Int) where {Ts} =
403 : sample(ary, times, Val(n))
404 :
405 : """
406 : out = sample(sink, range(0, 1, length=4))
407 : out = sample(sink, [0, 0.25, 0.5, 0.75, 1])
408 :
409 : Sample solution at given times during integration. Samples can be specified,
410 : e.g., as range or `Vector`.
411 :
412 : out = sample(array, times[, Val(P)])
413 :
414 : Same as above but create an [`indexedfill`](@ref) sink, i.e., store
415 : events in `array` (`P` defines a projection, i.e., `(t, y, dy)[P]`
416 : refers to the posotion `y` for `P=2`.
417 :
418 : !!! warning "precondition"
419 : The sequence of samples must be *increasing* for forward integration
420 : **or** *decreasing* for backward integration.
421 :
422 : See also [`RK43.Sample`](@ref), [`RK43.AbstractOutput`](@ref), [`RK43.AbstractSink`](@ref),
423 : [`indexedfill`](@ref)
424 : """ sample
425 :
426 : #-----------------------------------------------------------------------------
427 :
428 : """
429 : p = pipeline(out1, out2, ...)
430 : p(solver)
431 : # calls and checks out1(solver), out2(solver), ...
432 :
433 : Cascade outputs to a new "pipeline output".
434 :
435 : `Pipeline` stages are processed in order. An action onto a `Pipeline`
436 : is successful only if all stages yield `Ok`:
437 :
438 : - The next stage is evaluated only if the current stage returned `Ok`.
439 : - If the current stage does not return `OK`, its return value is
440 : returned by the pipeline.
441 :
442 : See also [`RK43.AbstractOutput`](@ref), [`pipeline`](@ref)
443 : """
444 : struct Pipeline{Args} <: AbstractOutput
445 : args::Args
446 4898 : function Pipeline{Args}(args) where {Args}
447 0 : @massert all(a -> isa(a, AbstractOutput), args)
448 4898 : new(args)
449 : end
450 : end
451 :
452 13 : _check_output(state::Bool) = state ? RK43.Ok : RK43.Stopped
453 1275827 : _check_output(state::RK43.State) = state
454 21 : _check_output(_) = RK43.Ok
455 :
456 196820 : function (ch::Pipeline{Args})(solver::RK43.Solver{T, N, S}) where {Args, T, N, S}
457 196820 : for a in ch.args
458 590468 : state = _check_output(a(solver))
459 393983 : RK43.issuccess(state) || return state
460 589783 : end
461 :
462 196480 : RK43.Ok
463 : end
464 :
465 4898 : Base.pipeline(out::Out, args...) where {Out<:AbstractOutput} =
466 : Pipeline{typeof(tuple(out, args...))}(tuple(out, args...))
467 :
468 1 : Base.pipeline(ch::Pipeline{Args}, out::Out) where {Args, Out<:AbstractOutput} =
469 : pipeline(ch.args..., out)
470 :
471 3 : Base.pipeline(out::Out, ch::Pipeline{Args}) where {Args, Out<:AbstractOutput} =
472 : pipeline(out, ch.args...)
473 :
474 1 : Base.pipeline(ch1::Pipeline{Args1}, ch2::Pipeline{Args2}) where {Args1, Args2} =
475 : pipeline(ch1.args..., ch2.args...)
476 :
477 : #-----------------------------------------------------------------------------
478 :
479 1 : Base.show(io::IO, ::Snk) where {Snk<:AbstractSink} = print(io, "sink")
480 :
481 1 : Base.show(io::IO, ::TryPush) = print(io, "pushinto")
482 :
483 :
484 1 : Base.show(io::IO, ::IndexedFill) = print(io, "indexedfill")
485 :
486 1 : Base.show(io::IO, ::Out) where {Out<:AbstractOutput} =
487 : print(io, Out)
488 :
489 3 : Base.show(io::IO, ::NoOutput) = print(io, "nooutput")
490 :
491 1 : Base.show(io::IO, out::Output{Snk}) where {Snk<:AbstractSink} =
492 : show(io, out.sink)
493 :
494 1 : function Base.show(io::IO, ch::Pipeline{Args}) where {Args}
495 1 : print(io, "Pipeline{")
496 :
497 1 : out= map(ch.args) do a
498 2 : buf = IOBuffer()
499 2 : show(buf, a)
500 2 : String(take!(buf))
501 : end
502 1 : print(io, join(out, ", "))
503 :
504 1 : print(io, "}")
505 : end
506 :
507 : #-----------------------------------------------------------------------------
508 :
509 : # TODO: [test] check performance (JET & BenchmarkTools, mind @assert and string interpolation)
510 :
511 : # TODO: [test] test dense sampling (forward, backward, edge cases)
|
