Coverage for pySDC/implementations/sweeper_classes/delta_form.py: 99%
206 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1r"""
2Deferred-correction ("delta-form") SDC sweepers.
4A standard SDC sweep
6.. math::
7 u^{k+1}_m = u_0 + \tau_m + \Delta t (Q f^k)_m
8 + \Delta t \sum_j Q^\Delta_{mj}\,(f^{k+1}_j - f^k_j)
10is algebraically identical to, with :math:`\delta_m = u^{k+1}_m - u^k_m` and the collocation
11residual :math:`r_m = u_0 + \tau_m + \Delta t (Q f^k)_m - u^k_m`,
13.. math::
14 \delta_m = r_m + \Delta t \sum_j Q^\Delta_{mj}\,\Delta f_j,
15 \qquad \Delta f_j = f(u^k_j + \delta_j) - f(u^k_j).
17Written this way, every sweep is iterative refinement: a high-precision residual, a correction
18solve, and a high-precision update ``u <- u + delta``. No Jacobian appears, so an IMEX splitting
19survives unchanged.
21The point of the reformulation is that the quantity handed to the node-local solver is a
22*correction*. Its magnitude tends to zero as the sweeps converge, so a reduced-precision solve
23introduces an error proportional to :math:`|\delta|` rather than to :math:`|u|` and therefore does
24not cap the attainable accuracy. What the sweeper and the solver guarantee each other for that to
25hold is the correction-solve contract, documented on :class:`~pySDC.core.problem.Problem`.
27Three node-local strategies are supported, selected automatically:
29``solve_system_delta``
30 Used when the problem provides it. Solves
31 :math:`\delta - \alpha[f(w+\delta) - f(w)] = r` for the correction directly. This is the only
32 option for a nonlinear implicit operator, and the only one that hands a reduced-precision
33 solver a small unknown.
35``linear_implicit=True``
36 For a linear or affine implicit operator, :math:`f(w+\delta) - f(w) = A\delta`, so the stock
37 ``solve_system`` already solves the correction equation once the affine part
38 :math:`\alpha f(0, t)` is removed from the right-hand side. No problem class needs changing.
39 The right-hand side is scaled to unit size around the call, so the solver may work in any
40 precision without the shrinking correction falling out of its format's range.
42fallback
43 Otherwise the substitution :math:`y = u^k_m + \delta_m` reduces the correction equation to the
44 ordinary implicit solve. Always correct and identical to
45 :class:`~pySDC.implementations.sweeper_classes.generic_implicit.generic_implicit`, but the solver sees an
46 :math:`\mathcal{O}(1)` unknown, so there is no precision benefit.
48``correction_precision`` additionally stores the small quantities
49(:math:`r`, :math:`\delta`, :math:`\Delta f`) in a reduced-precision datatype built from
50the problem's own ``init`` tuple, scaled by the residual's magnitude so the format's mantissa is
51used however small the correction gets.
53The increment :math:`\Delta f` is formed by subtracting two stored right-hand sides unless the
54problem provides ``eval_f_increment(base, delta, t)``, which expands it analytically. That
55subtraction is a cancellation carrying the operator norm, so it binds as soon as a level runs below
56backend precision.
58The same sweeper serves any number of levels, and which it is doing is decided by the transfer, not
59by the choice of sweeper. On the finest level it computes its own residual, which is where the
60accuracy of the whole iteration is set. Given a transfer that hands one down -- by setting
61``eps_in`` on the level's sweeper -- it uses that instead and advances it in place, so nothing is
62ever rebuilt out of :math:`\mathcal{O}(1)` coarse state. :class:`~pySDC.core.base_transfer.BaseTransfer` hands
63nothing down, so on a stock hierarchy this sweeper reproduces MLSDC and PFASST exactly.
64"""
66import numpy as np
68from pySDC.core.errors import ParameterError
69from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
70from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
73class DeltaFormMixin:
74 """Shared machinery for the delta-form sweepers."""
76 def _delta_setup(self):
77 """Read the optional sweeper parameters, and arm the per-sweep recorders."""
78 token = getattr(self.params, 'correction_precision', None)
79 self._work_dtype = None if token is None else np.dtype(token)
80 init = self.level.prob.init
81 if (
82 self._work_dtype is not None
83 and self._work_dtype.kind != 'c'
84 and isinstance(init, tuple)
85 and np.dtype(init[-1]).kind == 'c'
86 ):
87 # storing a complex correction in a real array keeps only its real part
88 raise ParameterError(
89 f'correction_precision={self._work_dtype} is real, but the problem is complex ({np.dtype(init[-1])}). '
90 'Use a complex correction_precision such as complex64.'
91 )
92 self._linear_implicit = bool(getattr(self.params, 'linear_implicit', False))
93 self._deltas, self._dfs = [], []
95 eps_in = None
96 """Residual handed down by the transfer, or ``None`` on the finest level."""
98 delta_acc = None
99 """Corrections this level has accumulated since the last restriction."""
101 def sync_initial_value(self):
102 r"""
103 Follow a change of :math:`u_0` made after the residual was handed down.
105 PFASST receives the initial value from the predecessor *after* the restriction, directly
106 into ``u[0]``. The residual depends on it additively, so following it is one addition:
107 :math:`r_m \leftarrow r_m + (u_0 - u_0^{\mathrm{ref}})`. Exactly zero
108 when nothing arrived, which is every serial run.
110 This is the one place the hierarchy still differences two :math:`\mathcal{O}(1)` values, and
111 it is the reason a reduced-precision coarse level buys less under PFASST than under MLSDC:
112 the difference is small -- it is the coarse-versus-fine discrepancy at the step interface,
113 and it converges to zero -- but it is *formed* by cancelling two values of size
114 :math:`|u|`. Putting the step-to-step exchange itself in delta form is what would remove it.
116 Returns
117 -------
118 None
119 """
120 lvl = self.level
121 if self.eps_in is None or lvl.u0_reference is None:
122 return
123 shift = lvl.u[0] - lvl.u0_reference
124 self.eps_in = [eps + shift for eps in self.eps_in]
125 lvl.u0_reference = lvl.prob.dtype_u(lvl.u[0])
126 return None
128 def compute_residual(self, stage=''):
129 r"""
130 Report the residual this level is already tracking, instead of rebuilding one.
132 A level with an inherited residual carries it forward through every sweep and prolongation,
133 so recomputing it from :math:`\mathcal{O}(1)` state would be both redundant and less
134 accurate. It is also what makes the :math:`\tau` term unnecessary: the only other reader is
135 :meth:`compute_end_point`, and only when the end point comes from the quadrature update.
137 Parameters
138 ----------
139 stage : str
140 The stage of the step this level belongs to.
142 Returns
143 -------
144 None
145 """
146 if self.eps_in is None:
147 return super().compute_residual(stage=stage)
149 lvl = self.level
150 if stage in self.params.skip_residual_computation:
151 lvl.status.residual = 0.0 if lvl.status.residual is None else lvl.status.residual
152 return None
154 lvl.residual = [lvl.prob.dtype_u(eps) for eps in self.eps_in]
155 norms = [abs(eps) for eps in self.eps_in]
156 kind = lvl.params.residual_type
157 if kind not in ('full_abs', 'last_abs', 'full_rel', 'last_rel'):
158 raise NotImplementedError(f'residual type "{kind}" not implemented!')
159 value = norms[-1] if kind.startswith('last') else max(norms)
160 lvl.status.residual = value / abs(lvl.u[0]) if kind.endswith('rel') else value
161 lvl.status.updated = False
162 return None
164 def advance_residual(self, eps, deltas, dfs):
165 r"""
166 Advance a residual by an update: :math:`r \leftarrow r - \delta
167 + \Delta t (Q \Delta f)`.
169 Every term is small, so this never cancels. It is exact for any update to the nodal values,
170 which is why it serves both a sweep and a prolongation.
172 Parameters
173 ----------
174 eps : list
175 The residual, one entry per node.
176 deltas : list
177 The update applied to the nodal values, one entry per node.
178 dfs : list
179 The resulting right-hand side increments, one entry per node.
181 Returns
182 -------
183 list
184 The advanced residual.
185 """
186 dt, Q = self.level.dt, self.coll.Qmat
187 out = []
188 for m in range(len(eps)):
189 acc = eps[m] - deltas[m]
190 for j in range(len(dfs)):
191 if Q[m + 1, j + 1] != 0.0:
192 acc += self._coeff(dt * Q[m + 1, j + 1]) * dfs[j]
193 out.append(acc)
194 return out
196 def accumulate(self, deltas):
197 """Add one round of corrections to what this level owes upwards."""
198 self.delta_acc = (
199 deltas if self.delta_acc is None else [a + d for a, d in zip(self.delta_acc, deltas, strict=True)]
200 )
202 def update_nodes(self):
203 r"""
204 Sweep, and keep the level's bookkeeping straight if it is part of a hierarchy.
206 The sweep itself is :meth:`_sweep_nodes`, which each concrete sweeper supplies. Around it:
207 follow any change of :math:`u_0` that arrived after the residual was handed down, advance
208 that residual by the update just applied, and bank the corrections for a transfer to prolong.
210 A level that computes its own residual has no bookkeeping to do, so on a single-level run
211 every line below the sweep is a no-op. That is why there is one sweeper rather than two.
213 Returns
214 -------
215 None
216 """
217 self.sync_initial_value()
218 self._sweep_nodes()
219 if self.eps_in is None:
220 return None
222 deltas = [self._to_work(self.level.prob, d) for d in self._deltas]
223 self.accumulate(deltas)
224 self.eps_in = self.advance_residual(self.eps_in, deltas, self._dfs)
225 return None
227 _work_scale = 1.0
228 r"""Shared divisor applied to the correction quantities before they are stored."""
230 def _work_init(self, prob):
231 """Build the problem's ``init`` tuple with the correction dtype substituted."""
232 return (prob.init[0], prob.init[1], self._work_dtype)
234 def _scales_corrections(self):
235 """
236 Whether this sweeper may choose a scale for its correction quantities.
238 A level that computes its own residual may: everything it stores is derived from that
239 residual within the same sweep, so one divisor keeps them all commensurate and ordinary
240 arithmetic on them stays correct. A level that *inherits* a residual may not, because the
241 inherited value was scaled by whoever produced it and is carried across sweeps.
242 """
243 return self.eps_in is None
245 def _set_work_scale(self, values):
246 r"""
247 Choose the divisor for this sweep, from the residual the corrections will be built out of.
249 This is what lets a correction be stored below ``float16``'s smallest normal, 6.1e-5. The
250 delta form drives :math:`r` and :math:`\delta` towards zero on purpose, and half
251 precision has almost no mantissa left down there -- 1.3e-2 relative at 1e-6, 1.9e-1 at 1e-7 --
252 so an unscaled correction turns to noise exactly when it starts to matter. Dividing by the
253 residual's own magnitude keeps the stored values at :math:`\mathcal{O}(1)`, which is block
254 floating point, and is what half-precision hardware does anyway.
256 Parameters
257 ----------
258 values : list
259 The residual at the collocation nodes, in backend units.
260 """
261 if self._work_dtype is None or not self._scales_corrections():
262 self._work_scale = 1.0
263 return
264 biggest = max((abs(value) for value in values), default=0.0)
265 self._work_scale = float(biggest) if biggest > 0.0 else 1.0
267 def _to_work(self, prob, value):
268 """
269 Store a small correction quantity in a reduced-precision datatype.
271 Returns the value unchanged when no reduced precision was requested, so the default path
272 makes no assumption about the datatype and works with any pySDC backend.
274 Raises
275 ------
276 NotImplementedError
277 If ``correction_precision`` was requested but the datatype cannot be built at another
278 precision, as is the case for datatypes not backed by a numpy array.
279 """
280 if self._work_dtype is None:
281 return value
282 try:
283 me = prob.dtype_u(self._work_init(prob))
284 me[:] = value if self._work_scale == 1.0 else value / self._work_scale
285 except (TypeError, NotImplementedError) as error:
286 raise NotImplementedError(
287 f'correction_precision is not supported for {prob.dtype_u.__name__}: it cannot be '
288 f'built at a different precision from the problem init tuple ({error})'
289 ) from error
290 return me
292 def _to_backend(self, prob, value):
293 """
294 Lift a possibly reduced-precision quantity back to backend precision.
296 A no-op when no reduced precision is in play, which keeps the default path datatype-agnostic.
297 """
298 if self._work_dtype is None:
299 return value
300 me = prob.dtype_u(prob.init)
301 me[:] = value
302 if self._work_scale != 1.0:
303 # widen first, scale second. The other order multiplies at the reduced precision, and a
304 # scale of 1e-8 then lands the result below float16's smallest subnormal on the way out
305 # -- the value is destroyed before it ever reaches the backend-precision array.
306 me *= self._work_scale
307 return me
309 def _coeff(self, value):
310 """
311 Cast a scalar coefficient so an accumulation stays at the precision it should.
313 The ``float`` matters. A coefficient taken out of a numpy array is an ``np.float64``, and
314 multiplying a reduced-precision level quantity by one of those upcasts the result to
315 ``float64`` under NEP 50 -- NumPy 2's rule -- which would quietly put the whole level back at
316 backend precision. A plain Python float is weak under both the old and the new rule, so the
317 array's own dtype wins.
318 """
319 if self._work_dtype is None:
320 return float(value)
321 return self._work_dtype.type(value)
323 def _residual_nodes(self):
324 r"""
325 Compute :math:`r_m = u_0 + \tau_m + \Delta t (Q f^k)_m - u^k_m`.
327 This is the high-precision residual of iterative refinement. It is a difference of
328 :math:`\mathcal{O}(1)` quantities and is therefore always formed in backend precision --
329 unless a transfer handed one down, in which case that one is already exact and small, and
330 rebuilding it here is what the delta-form hierarchy exists to avoid.
332 Returns
333 -------
334 list
335 One residual per collocation node.
336 """
337 if self.eps_in is not None:
338 return self.eps_in
340 lvl = self.level
341 eps = self.integrate()
342 for m in range(self.coll.num_nodes):
343 eps[m] += lvl.u[0]
344 eps[m] -= lvl.u[m + 1]
345 if lvl.tau[m] is not None:
346 eps[m] += lvl.tau[m]
347 return eps
349 def _f_increment(self, prob, f_new, f_old, u_old, delta, t_node):
350 r"""
351 Form the right-hand side increment :math:`\Delta f = f(w+\delta) - f(w)`.
353 Formed by subtraction unless the problem can expand it analytically. The subtraction is a
354 cancellation of two :math:`\mathcal{O}(|f|)` quantities, so its absolute error is
355 :math:`\varepsilon |f|` in whatever precision the level stores ``f`` at, and :math:`|f|`
356 carries the operator norm. That is harmless while the level is at backend precision and
357 becomes the binding term as soon as it is not, which is why a problem living on a
358 reduced-precision level should provide ``eval_f_increment``.
360 Parameters
361 ----------
362 prob : pySDC.core.problem.Problem
363 The problem on this level.
364 f_new, f_old : dtype_f
365 ``f`` at :math:`w+\delta` and at :math:`w`, used by the subtraction fallback.
366 u_old : dtype_u
367 The base state :math:`w`.
368 delta : dtype_u
369 The correction :math:`\delta`.
370 t_node : float
371 Physical time of the collocation node.
373 Returns
374 -------
375 dtype_f
376 The increment, with the same splitting as ``eval_f``.
377 """
378 if hasattr(prob, 'eval_f_increment'):
379 increment = prob.eval_f_increment(u_old, delta, t_node)
380 else:
381 increment = prob.dtype_f(f_new)
382 increment -= f_old
383 # recorded for the residual recursion, so a transfer never recovers it by subtraction
384 self._dfs.append(total_increment(prob, increment))
385 return increment
387 def _solve_correction(self, rhs_corr, alpha, u_old, f_old, t_node, implicit_part=None):
388 r"""
389 Solve the node-local correction equation.
391 Parameters
392 ----------
393 rhs_corr : dtype_u
394 Right-hand side :math:`r` of the correction equation.
395 alpha : float
396 Implicit prefactor :math:`\alpha = \Delta t Q^\Delta_{mm}`.
397 u_old : dtype_u
398 Current nodal value :math:`u^k_m`, the base state of the correction.
399 f_old : dtype_f
400 ``f`` evaluated at ``u_old``; already stored on the level, so it costs nothing.
401 t_node : float
402 Physical time of the collocation node.
403 implicit_part : dtype_u, optional
404 The implicit component of ``f_old`` for IMEX problems. Defaults to ``f_old``.
406 Returns
407 -------
408 dtype_u
409 The correction :math:`\delta_m`.
410 """
411 prob = self.level.prob
412 f_impl_old = f_old if implicit_part is None else implicit_part
414 rhs_phys = self._to_backend(prob, rhs_corr)
416 if alpha == 0:
417 # explicit node: the correction is the residual itself
418 delta = rhs_phys
419 elif hasattr(prob, 'solve_system_delta'):
420 delta = prob.solve_system_delta(rhs_phys, alpha, u_old, f_old, t_node)
421 elif self._linear_implicit:
422 # f(w+d) - f(w) = A d, so solve_system already solves the correction equation once the
423 # affine part f(0, t) has been removed. f(0, t) vanishes for a homogeneous operator.
424 affine = prob.eval_f(prob.dtype_u(prob.init, val=0.0), t_node)
425 affine = affine if implicit_part is None else affine.impl
426 delta = self._solve_scaled(prob, rhs_phys, alpha, affine, t_node)
427 else:
428 # Fallback: substitute y = u_old + delta. Always correct, but the solver sees an O(1)
429 # unknown, so there is no precision benefit.
430 rhs_phys += u_old
431 rhs_phys -= alpha * f_impl_old
432 solution = prob.solve_system(rhs_phys, alpha, u_old, t_node)
433 delta = prob.dtype_u(solution)
434 delta -= u_old
435 # recorded so a transfer never has to recover the correction by subtraction
436 self._deltas.append(delta)
437 return delta
439 @staticmethod
440 def _solve_scaled(prob, rhs, alpha, affine, t_node):
441 r"""
442 Hand a linear correction solve its right-hand side at unit size, and take back any precision.
444 This is the sweeper's side of the correction-solve contract documented on
445 :class:`~pySDC.core.problem.Problem`. The right-hand side shrinks with the iteration, which
446 is the point of the delta form and also what pushes it out of a reduced-precision format's
447 range: a correction of 1e-10 is below ``float16``'s smallest subnormal. Dividing by its
448 maximum before the solve and multiplying the result by it afterwards is exact for a linear
449 solve, and makes an absolute solver tolerance a relative one, which the contract asks for
450 anyway. The affine part :math:`\alpha f(0, t)` is removed *after* scaling: the stock solve
451 adds it back in unscaled, and the two only cancel at the same scale. The result may come back in any precision; it is widened before it is scaled back,
452 since scaling first would multiply at the solver's precision.
454 Parameters
455 ----------
456 prob : pySDC.core.problem.Problem
457 The level's problem.
458 rhs : dtype_u
459 Right-hand side of the correction equation, modified in place.
460 alpha : float
461 Implicit prefactor.
462 affine : dtype_u
463 The implicit right-hand side at zero, :math:`f(0, t)`, which a homogeneous operator
464 makes zero.
465 t_node : float
466 Physical time of the collocation node.
468 Returns
469 -------
470 dtype_u
471 The correction, at the level's precision.
472 """
473 zero = prob.dtype_u(prob.init, val=0.0)
474 scale = float(abs(rhs))
475 if scale == 0.0:
476 return zero
477 rhs /= scale
478 rhs -= alpha * affine
479 solution = prob.solve_system(rhs, alpha, zero, t_node)
480 if isinstance(prob.init, tuple) and getattr(solution, 'dtype', None) not in (None, np.dtype(prob.init[-1])):
481 delta = prob.dtype_u(prob.init)
482 delta[:] = solution
483 else:
484 delta = prob.dtype_u(solution)
485 delta *= scale
486 return delta
489class delta_implicit(DeltaFormMixin, generic_implicit):
490 """
491 Delta-form counterpart of :class:`~pySDC.implementations.sweeper_classes.generic_implicit.generic_implicit`.
493 Mathematically identical to the standard sweep; see the module docstring for the sweeper
494 parameters ``correction_precision`` and ``linear_implicit``.
495 """
497 def _sweep_nodes(self):
498 """
499 Perform one delta-form sweep over all collocation nodes.
501 Returns
502 -------
503 None
504 """
505 lvl = self.level
506 prob = lvl.prob
507 assert lvl.status.unlocked
508 num_nodes = self.coll.num_nodes
509 self._delta_setup()
511 residual = self._residual_nodes()
512 self._set_work_scale(residual)
513 eps = [self._to_work(prob, value) for value in residual]
514 df = [None] * (num_nodes + 1)
516 for m in range(num_nodes):
517 t_node = lvl.time + lvl.dt * self.coll.nodes[m]
519 rhs_corr = type(eps[m])(eps[m])
520 for j in range(1, m + 1):
521 if self.QI[m + 1, j] != 0.0:
522 rhs_corr += self._coeff(lvl.dt * self.QI[m + 1, j]) * df[j]
524 alpha = lvl.dt * self.QI[m + 1, m + 1]
525 u_old = prob.dtype_u(lvl.u[m + 1])
526 f_old = prob.dtype_f(lvl.f[m + 1])
528 delta = self._solve_correction(rhs_corr, alpha, u_old, f_old, t_node)
530 lvl.u[m + 1] = u_old + self._to_backend(prob, self._to_work(prob, delta))
531 lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], t_node)
533 increment = self._f_increment(prob, lvl.f[m + 1], f_old, u_old, delta, t_node)
534 df[m + 1] = self._to_work(prob, increment)
536 lvl.status.updated = True
537 return None
540class delta_imex_1st_order(DeltaFormMixin, imex_1st_order):
541 """
542 Delta-form counterpart of :class:`~pySDC.implementations.sweeper_classes.imex_1st_order.imex_1st_order`.
544 The correction equation contains only differences of ``f``, never a Jacobian, so the
545 explicit/implicit splitting is untouched.
546 """
548 def _sweep_nodes(self):
549 """
550 Perform one delta-form IMEX sweep over all collocation nodes.
552 ``QE`` is strictly lower triangular, which
553 :class:`~pySDC.implementations.sweeper_classes.imex_1st_order.imex_1st_order` already enforces, so the
554 explicit part never contributes to the node-local solve.
556 Returns
557 -------
558 None
559 """
560 lvl = self.level
561 prob = lvl.prob
562 assert lvl.status.unlocked
563 num_nodes = self.coll.num_nodes
564 self._delta_setup()
566 residual = self._residual_nodes()
567 self._set_work_scale(residual)
568 eps = [self._to_work(prob, value) for value in residual]
569 df_impl = [None] * (num_nodes + 1)
570 df_expl = [None] * (num_nodes + 1)
572 for m in range(num_nodes):
573 t_node = lvl.time + lvl.dt * self.coll.nodes[m]
575 rhs_corr = type(eps[m])(eps[m])
576 for j in range(1, m + 1):
577 if self.QI[m + 1, j] != 0.0:
578 rhs_corr += self._coeff(lvl.dt * self.QI[m + 1, j]) * df_impl[j]
579 if self.QE[m + 1, j] != 0.0:
580 rhs_corr += self._coeff(lvl.dt * self.QE[m + 1, j]) * df_expl[j]
582 alpha = lvl.dt * self.QI[m + 1, m + 1]
583 u_old = prob.dtype_u(lvl.u[m + 1])
584 f_old = prob.dtype_f(lvl.f[m + 1])
586 delta = self._solve_correction(rhs_corr, alpha, u_old, f_old, t_node, implicit_part=f_old.impl)
588 lvl.u[m + 1] = u_old + self._to_backend(prob, self._to_work(prob, delta))
589 lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], t_node)
591 increment = self._f_increment(prob, lvl.f[m + 1], f_old, u_old, delta, t_node)
592 df_impl[m + 1] = self._to_work(prob, prob.dtype_u(increment.impl))
593 df_expl[m + 1] = self._to_work(prob, prob.dtype_u(increment.expl))
595 lvl.status.updated = True
596 return None
599def total_increment(prob, increment):
600 """
601 The full right-hand side increment, recombining an IMEX splitting.
603 Parameters
604 ----------
605 prob : pySDC.core.problem.Problem
606 The problem the increment belongs to.
607 increment : dtype_f
608 The increment, possibly split into ``impl`` and ``expl``.
610 Returns
611 -------
612 dtype_u
613 The sum of the parts.
614 """
615 if not hasattr(increment, 'impl'):
616 return increment
617 return prob.dtype_u(increment.impl) + increment.expl