From ddee760a5f4f4cf3794cb5df4bf576fd4728d10e Mon Sep 17 00:00:00 2001 From: derrynknife Date: Sat, 17 Oct 2026 21:57:28 +1000 Subject: [PATCH] A joint model's unit as a block of a non-repairable diagram A unit of a surpyval joint model of readings and failures, given its readings, is a node of a NonRepairableRBD: a JointNode, or the unit's DynamicPrediction given as the node. Its reliability is from now, pred.sf(age + t); its lifetimes pred.qf(u) - age; its mean the unit's mean life left. A NodeState with another age than the prediction's is refused, as it would condition the unit again. A RepairableRBD and a NonRepairable refuse a joint model: the units after the first are new, without its readings. It saves with the diagram (the model in surpyval's format and the unit's inputs). surpyval 0.25 is required, for its joint models. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_019Km2eWBiPJ5JNkme3vRg4B --- CHANGELOG.md | 18 ++ CLAUDE.md | 9 +- CONTRIBUTING.md | 2 +- README.md | 6 +- docs/api/nodes.md | 2 + docs/guide/condition-based.md | 89 ++++++ pyproject.toml | 14 +- repyability/__init__.py | 2 + repyability/non_repairable.py | 3 + repyability/rbd/_spec.py | 3 + repyability/rbd/joint_node.py | 407 ++++++++++++++++++++++++++ repyability/rbd/non_repairable_rbd.py | 38 ++- repyability/rbd/routes.py | 13 + repyability/rbd/serialisation.py | 10 +- repyability/tests/capabilities.py | 1 + repyability/tests/catalogue.py | 41 +++ repyability/tests/test_catalogue.py | 4 + repyability/tests/test_joint_node.py | 215 ++++++++++++++ 18 files changed, 862 insertions(+), 15 deletions(-) create mode 100644 repyability/rbd/joint_node.py create mode 100644 repyability/tests/test_joint_node.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 39e74b5..7bca3af 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,24 @@ other release, fixes included, the minor. ## [Unreleased] +### Added + +- A unit of a surpyval joint model of readings and failures, given its + readings, is a block of a `NonRepairableRBD`: a `JointNode`, or the + unit's `DynamicPrediction` (`joint_model.predict(Z, markers, given)`) + given as the node. Its reliability is from now, `pred.sf(age + t)`; its + lifetimes `pred.random() - age`; its mean the unit's mean life left. It + takes part in the reliability, MTTF, B-lives, importance measures, + states, simulations and saving. Its age is its prediction's: a + `NodeState` with another age is refused, as it would condition the unit + again. A `RepairableRBD` and a `NonRepairable` refuse a joint model, + whose units after the first are new and without its readings. See the + condition-based evaluation guide. + +### Changed + +- surpyval 0.25 is required (for its joint models). + ## [0.14] - 2026-10-16 Fitted models in, every kind of model brought up to the diagram, and diff --git a/CLAUDE.md b/CLAUDE.md index 2e859f2..6a05810 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -4,9 +4,9 @@ Setup, branches, the checks a change needs, where tests go, pull requests, deprecations and releases are in [CONTRIBUTING.md](CONTRIBUTING.md); read it first. This file adds the internals and conventions it does not cover. -**Versions (17 Oct 2026):** 0.13 is the latest release (PyPI, `v0.13`). 0.14 -is prepared on `dev` (`_version.py`, CHANGELOG dated 2026-10-16) and waits on -the release pull request from `dev` into `main`. `dev` needs surpyval 0.24 +**Versions (17 Oct 2026):** 0.14 is the latest release (PyPI, `v0.14`). 0.15 +is prepared on `dev` (`_version.py`, CHANGELOG dated 2026-10-17) and waits on +the release pull request from `dev` into `main`. `dev` needs surpyval 0.25 or later. ## Architecture & scope @@ -483,7 +483,8 @@ CONTRIBUTING.md. What it leaves to the maintainer's say: ## surpyval workarounds to remove -surpyval 0.24 (the minimum) gives a mixture a precise tail (`sf`, #671) and +surpyval 0.25 (the minimum) has the joint models a `JointNode` takes; 0.24 +gives a mixture a precise tail (`sf`, #671) and has `conditional_gaps` read a limited failure population's `lfp_p`, not a mixture's `p` (#626); 0.23 pickles its fits (#573), names the limited-failure proportion `lfp_p` (#608) and takes `success_run`'s diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 133dca4..853af90 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -18,7 +18,7 @@ pip install -e .[dev] # add ,fast for numba, ,docs for mkdocs pre-commit install # run the formatters and linters on commit ``` -RePyability needs Python 3.11 or later and surpyval 0.24 or later. +RePyability needs Python 3.11 or later and surpyval 0.25 or later. ## Branches and the one automated test run diff --git a/README.md b/README.md index 52718a8..0052993 100755 --- a/README.md +++ b/README.md @@ -226,6 +226,7 @@ How each family of analyses is worked out on a diagram with each kind of node or | A probability per demand | E | – | – | E | E | N | S | E N | – | S | | A mixture (`MixtureModel`) | E | N | N | E N | E N | N | S | E N | – | S | | A regression node (covariates) | E | N | N | E N | E N | N | S | E N | – | S | +| A joint model's unit, given its readings | N | N | N | N | N | N | S | N | – | S | | A degrading node | N | N | N | N | N | N | S | N | N | S | | A fitted degradation process | E | N | N | E N | E N | N | S | E N | – | S | | Cold standby (convolved) | N | N | N | N | N | N | S | N | – | S | @@ -303,8 +304,9 @@ places, nested diagrams): - **components:** surpyval lifetime distributions, degradation processes (Wiener and gamma, as the time to their threshold) and fixed probabilities (a diagram refuses a non-parametric fit, such as - Kaplan–Meier: fit a parametric distribution in surpyval); repeated - nodes; cold standby with one or + Kaplan–Meier: fit a parametric distribution in surpyval); a unit of a + surpyval joint model of readings and failures, given its readings + (`JointNode`, from now); repeated nodes; cold standby with one or two units operating (any units) or more (identical units); warm standby with one unit operating, and hot standby (any units); load sharing of identical units; common-cause groups, for the system's reliability, diff --git a/docs/api/nodes.md b/docs/api/nodes.md index c353c52..cd91479 100644 --- a/docs/api/nodes.md +++ b/docs/api/nodes.md @@ -18,6 +18,8 @@ models provided here (see [Building an RBD](../guide/building.md) and ::: repyability.RegressionNode +::: repyability.JointNode + ::: repyability.PerfectReliability ::: repyability.PerfectUnreliability diff --git a/docs/guide/condition-based.md b/docs/guide/condition-based.md index 0b8fb65..e56f057 100644 --- a/docs/guide/condition-based.md +++ b/docs/guide/condition-based.md @@ -111,6 +111,9 @@ and the valve works. - Parametric distributions and `RegressionNode`s age. For them, `NodeState(age=0)` is the same as leaving the node out. +- A [`JointNode`][repyability.JointNode] (a unit of a joint model, given + its readings) takes its age from its prediction: see + [below](#a-unit-with-readings-joint-models). - A fixed-probability component given any state in which it is alive contributes 1 from now on: its chance of failing has been resolved by observing it working. Left out of the state, it keeps its usual @@ -228,3 +231,89 @@ ramped.sf_given_state(1000, {"m": NodeState(age=2500)}) # -> 0.7036 semi-parametric Cox model has no defined mean and says so. - The node saves with the RBD; the fitted model round-trips through `surpyval.from_dict`. + +## A unit with readings: joint models + +Some components are watched: a vibration level, a wear measurement or an +oil analysis is read every so often, and a rising reading warns of a +failure to come. A **joint model** of the readings and the failures, +fitted in surpyval (`surpyval.WeibullJointPH`, or any +`surpyval.JointPH(distribution, path_model)`), says how. Its `predict` +gives one unit's life given its covariates, its readings so far and that +it is still running at its age now, and that prediction can be a node: + +```python +from surpyval import ItemLifetime, SurpyvalData, WeibullJointPH +from repyability import JointNode + +# A fleet of 60 pumps, each with its vibration read every 10 hours. +rng = np.random.default_rng(3) +pumps = [] +for k in range(60): + start, rate = rng.normal(2.0, 0.4), rng.normal(0.02, 0.005) + life = 100 * rng.weibull(2.0) * np.exp(-0.9 * (start - 2)) + end = min(life, 120.0) + ages = np.arange(0.0, end, 10.0) + vibration = start + rate * ages + rng.normal(0, 0.2, ages.size) + outcome = "observed" if life < 120 else "right_censored" + pumps.append( + ItemLifetime(k, **{outcome: end}, markers=dict(zip(ages, vibration))) + ) +joint = WeibullJointPH.fit(SurpyvalData.from_items(pumps)) + +# Our pump: 40 hours old, its vibration rising faster than most. +pump = joint.predict(markers={0: 2.1, 20: 2.5, 40: 3.0}, given=40) +valve = surv.Weibull.from_params([80, 1.5]) +line = NonRepairableRBD( + [("in", "pump"), ("pump", "valve"), ("valve", "out")], + {"pump": pump, "valve": valve}, +) +line.sf(20) # -> 0.6953 pump.sf(60) * valve.sf(20) +``` + +The diagram's time is **time from now**, while a prediction's functions +take ages since new. The node is the prediction moved to now: at mission +time `t` the pump's reliability is `pump.sf(40 + t)`, which surpyval has +already conditioned on the pump reaching 40 hours, its lifetimes are +`pump.random() - 40` and its mean is its mean life left. A prediction given +as a node becomes a [`JointNode`][repyability.JointNode]; +`JointNode(joint, Z=..., markers=..., age=...)` makes the same node from the +fitted model and the unit's covariates, readings and age. + +It takes part in everything an ordinary node does: the system's +reliability, MTTF and B-lives, the importance measures and the +simulations. + +```python +line.mean() # -> 35.50 +line.birnbaum_importance(20)["valve"] # -> 0.7879 the pump's sf(60) +line.sf_given_state(20, {"valve": NodeState(age=30)}) # -> 0.6048 +``` + +- **Its age is the prediction's.** In a state, the pump may be given + `alive=False`, or its prediction's age (`NodeState(age=40)`, the same as + leaving it out); any other age raises `ValueError`, as it would condition + the pump a second time. For the pump later on, predict again with + `given=` its new age and its readings since. +- **Its mean can be infinite.** With a straight-line marker path whose + slope varies from unit to unit, a few units' paths fall away and they + never fail; the model leaves even this pump a chance of that (about + 5e-9), so `JointNode(pump).mean()` is `inf`, as the mean of such a + model is. The reliability, the B-lives and a system whose other parts + wear out (`line` above) are not affected. +- **Parameter uncertainty** (`sf_uncertainty` and the others) holds the + joint model at its fit, unless the pump is given a list of predictions + to draw from (from joint models refitted to bootstrap resamples, say); + surpyval's `pump.cb(...)` bounds the pump's own reliability for its + model's. +- **A new unit** is `JointNode(joint, Z=...)`, with no readings: the + model's life over the fleet's marker paths. +- **Not in a repairable diagram.** A [`RepairableRBD`][repyability.RepairableRBD] + (and a `NonRepairable`) refuses a joint model: the prediction is of the + unit in service now, given its readings, but the units that replace it + are new and have none, while such a component draws every one of its + lives from one surpyval distribution. Use a `NonRepairableRBD` for the + mission from now. +- The node saves with the RBD: the joint model in surpyval's format + (`JointModel.to_dict()`, a built-in path model) and the unit's + covariates, readings and age. diff --git a/pyproject.toml b/pyproject.toml index e561c58..d07f2f4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,12 +43,14 @@ dependencies = [ # (sf/ff/qf/random/mean, extras/with_params, regression sf/sf_tvc, # to_dict/from_dict). A surpyval release reaches every fresh install at # once, so run the tests against surpyval's development branch before it - # is released, to show a break first. CI also tests this minimum. 0.24 - # gives a mixture a precise tail (sf) and has conditional_gaps read a - # limited failure population's lfp_p, not a mixture's p; 0.23 named the - # proportion lfp_p, takes success_run's level as alpha_ci, and pickles - # its fits (for n_jobs' worker processes), as RePyability relies on. - "surpyval>=0.24", + # is released, to show a break first. CI also tests this minimum. 0.25 + # has the joint models (surpyval.univariate.joint) a JointNode takes; + # 0.24 gives a mixture a precise tail (sf) and has conditional_gaps read + # a limited failure population's lfp_p, not a mixture's p; 0.23 named + # the proportion lfp_p, takes success_run's level as alpha_ci, and + # pickles its fits (for n_jobs' worker processes), as RePyability + # relies on. + "surpyval>=0.25", "tqdm>=4.64", ] diff --git a/repyability/__init__.py b/repyability/__init__.py index 5d3348d..4ece488 100644 --- a/repyability/__init__.py +++ b/repyability/__init__.py @@ -61,6 +61,7 @@ PerfectReliability, PerfectUnreliability, ) +from repyability.rbd.joint_node import JointNode from repyability.rbd.load_sharing_node import LoadSharingModel from repyability.rbd.node_state import NodeState from repyability.rbd.non_repairable_rbd import NonRepairableRBD @@ -130,6 +131,7 @@ "PerfectUnreliability", "NodeState", "RegressionNode", + "JointNode", "BetaFactor", "MGL", "AlphaFactor", diff --git a/repyability/non_repairable.py b/repyability/non_repairable.py index 0c18767..3886ac7 100644 --- a/repyability/non_repairable.py +++ b/repyability/non_repairable.py @@ -17,6 +17,7 @@ never_fails, ) from repyability.rbd.helper_classes import PerfectReliability, perfect_class +from repyability.rbd.joint_node import NOT_REPAIRABLE, is_joint from repyability.rbd.results import FailureForecast from repyability.rbd.standby_node import StandbyModel from repyability.utils.checks import is_number, no_distribution, real_array @@ -235,6 +236,8 @@ def __init__( # is only simulated, and its sf refuses. self.model_parameterization = "standby" self.reliability_function = _scalar_sf(reliability) + elif is_joint(reliability): + raise ValueError(NOT_REPAIRABLE) elif reliability is PerfectReliability: raise ValueError( "A NonRepairable's life must end, and PerfectReliability " diff --git a/repyability/rbd/_spec.py b/repyability/rbd/_spec.py index 0d913a8..072ee77 100644 --- a/repyability/rbd/_spec.py +++ b/repyability/rbd/_spec.py @@ -361,6 +361,7 @@ def _validate_groups( def _unknown_component(node, component) -> str: """Why ``component`` cannot be node ``node``: it is none of the kinds a node can be.""" + from repyability.rbd.joint_node import NOT_REPAIRABLE, is_joint from repyability.repairable import Repairable kinds = ( @@ -369,6 +370,8 @@ def _unknown_component(node, component) -> str: "that never fails, such as a k-out-of-n vote point, " "PerfectReliability)." ) + if is_joint(component): + return f"Component {node!r}: {NOT_REPAIRABLE}" if isinstance(component, Repairable): return ( f"Component {node!r} is a Repairable, which models the " diff --git a/repyability/rbd/joint_node.py b/repyability/rbd/joint_node.py new file mode 100644 index 0000000..ff9122c --- /dev/null +++ b/repyability/rbd/joint_node.py @@ -0,0 +1,407 @@ +"""Joint node: one unit of a fitted joint model, given its readings, as a +block of a non-repairable diagram. + +A surpyval joint model (``surpyval.WeibullJointPH``, or any +``surpyval.JointPH(dist, path_model)`` fitter) models a marker read +repeatedly on each unit (a vibration level, a wear measurement) together +with the unit's life. Its ``predict`` gives one unit's +``DynamicPrediction``: the unit's failure-age distribution given its +covariates, its readings so far and that it is still running at its age +now. A ``JointNode`` puts that unit into a diagram. + +Time in a diagram is time from now, and a prediction's functions take ages +since new. The node is the prediction moved to now: at mission time ``t`` +its reliability is ``pred.sf(age + t)``, which surpyval has already +conditioned on the unit surviving to ``age``; its lifetimes are +``pred.random() - age`` and its mean is the unit's mean residual life. The +age is the prediction's, so the node is not conditioned again on a +``NodeState`` age (see ``NonRepairableRBD.sf_given_state``). +""" + +from typing import Any, Optional + +import numpy as np +from numpy.typing import ArrayLike + +from ._sampling import RowSampler + +#: Why a repairable diagram takes no joint node, for its messages. +NOT_REPAIRABLE = ( + "A component that is replaced (a NonRepairable, or a component of a " + "RepairableRBD) takes no joint model as its life: the prediction is " + "of the unit in service now, given its readings, while each unit that " + "replaces it is new, without them, and every life of such a " + "component is drawn from one surpyval distribution. Use a " + "NonRepairableRBD for the mission from now (a JointNode), and give a " + "replaced component a distribution fitted to the population's lives." +) + + +def _joint_classes() -> tuple: + """surpyval's ``JointModel`` and ``DynamicPrediction``.""" + from surpyval.univariate.joint import DynamicPrediction, JointModel + + return JointModel, DynamicPrediction + + +def is_joint(model) -> bool: + """Whether ``model`` is a joint model's: a ``JointNode``, or surpyval's + ``JointModel`` or ``DynamicPrediction``.""" + return isinstance(model, (JointNode, *_joint_classes())) + + +def as_node(model): + """``model``, but a surpyval ``DynamicPrediction`` as the + ``JointNode`` of its unit (a diagram takes a prediction as a node).""" + _, prediction = _joint_classes() + if isinstance(model, prediction): + return JointNode(model) + return model + + +class JointNode: + """One unit of a fitted joint model, given its readings, as a node. + + surpyval's joint models (``surpyval.WeibullJointPH``) model a marker + read repeatedly on each unit together with the unit's life; a fitted + model's ``predict(Z, markers, given)`` is one unit's failure-age + distribution given its covariates, its readings so far and that it + is still running at age ``given``. This node is that unit in a + [`NonRepairableRBD`][repyability.NonRepairableRBD], from now: at time + ``t`` from now its reliability is ``pred.sf(age + t)`` (already + conditioned on the unit reaching ``age``), its lifetimes are + ``pred.random() - age`` and its ``mean`` is the unit's mean residual + life. A ``DynamicPrediction`` given as a diagram's node is taken as + the ``JointNode`` of its unit. + + The node's age is the prediction's: a ``NodeState`` for it in + ``sf_given_state`` and the other condition-based methods may say + that it has failed, or give that same age, but not another, which + would condition it again (predict again, with ``given`` the new age + and the readings since). A + [`RepairableRBD`][repyability.RepairableRBD] refuses it: the units + that replace this one are new, without its readings. See the guide's + joint-model page. + + Parameters + ---------- + model : surpyval DynamicPrediction or JointModel + The unit's prediction (``joint_model.predict(...)``), or the + fitted joint model, with the unit's ``Z``, ``markers`` and + ``age``. + Z : dict, array_like or surpyval ItemLifetime, optional + With a ``JointModel``: the unit's covariates (``{name: value}``, + or a row in the order of the model's ``feature_names``), needed + where the model has covariates; or an ``ItemLifetime`` still + running, which holds the unit's covariates, readings and age. + markers : dict, optional + With a ``JointModel``: the unit's readings so far, ``{age: + value}``. + age : float, optional + With a ``JointModel``: the unit's age now (surpyval's ``given``), + at or after its last reading; 0 (new) by default, when it has no + readings. + + Raises + ------ + TypeError + If ``model`` is neither, or ``Z``, ``markers`` or ``age`` is given + with a prediction, which holds them. + ValueError + If the model's checks refuse the unit's inputs, or a covariate is + missing (``nan``). + + Examples + -------- + Fit a joint model to a small fleet in surpyval, and put a pump whose + marker has risen into a diagram, 40 hours old: + + >>> import numpy as np + >>> from surpyval import ItemLifetime, SurpyvalData, WeibullJointPH + >>> from repyability import JointNode, NonRepairableRBD + >>> rng = np.random.default_rng(3) + >>> items = [] + >>> for k in range(60): + ... start, rate = rng.normal(2.0, 0.4), rng.normal(0.02, 0.005) + ... life = 100 * rng.weibull(2.0) * np.exp(-0.9 * (start - 2)) + ... end = min(life, 120.0) + ... ages = np.arange(0.0, end, 10.0) + ... y = start + rate * ages + rng.normal(0, 0.2, ages.size) + ... outcome = "observed" if life < 120 else "right_censored" + ... items.append( + ... ItemLifetime(k, **{outcome: end}, + ... markers=dict(zip(ages, y))) + ... ) + >>> model = WeibullJointPH.fit(SurpyvalData.from_items(items)) + >>> pred = model.predict(markers={0: 2.1, 20: 2.5, 40: 3.0}, given=40) + >>> pump = JointNode(pred) + >>> bool(pump.sf(20.0) == pred.sf(60.0)) + True + >>> rbd = NonRepairableRBD([("in", "pump"), ("pump", "out")], + ... {"pump": pump}) + >>> bool(rbd.sf(20.0) == pred.sf(60.0)) + True + """ + + def __init__( + self, + model: Any, + Z: Any = None, + markers: Optional[dict] = None, + age: Optional[float] = None, + ): + joint_model, prediction = _joint_classes() + if isinstance(model, prediction): + given = [ + name + for name, value in ( + ("Z", Z), + ("markers", markers), + ("age", age), + ) + if value is not None + ] + if given: + raise TypeError( + "A JointNode of a DynamicPrediction takes the unit's " + "covariates, readings and age from the prediction: " + f"leave out {', '.join(given)} (or give the " + "JointModel with them)." + ) + pred = model + elif isinstance(model, joint_model): + if markers and age is None: + raise ValueError( + "Give the unit's age now, age=..., at or after its " + "last reading: the readings alone do not say it." + ) + pred = model.predict(Z, markers, given=age) + else: + raise TypeError( + "A JointNode takes a surpyval joint model's prediction for " + "one unit (joint_model.predict(Z, markers, given)), or the " + "fitted JointModel with the unit's Z, markers and age; got " + f"a {type(model).__name__}." + ) + if not np.all(np.isfinite(np.asarray(pred.covariates, dtype=float))): + raise ValueError( + "The unit's covariates are missing (nan), so its life has " + "no prediction: give each of them." + ) + self.prediction = pred + self.model = pred.model + self.age = float(pred.age) + + # -- Node reliability interface --------------------------------------- + + def _at(self, name: str, x: ArrayLike) -> Any: + """The prediction's function ``name`` at the ages ``age + x`` (the + times ``x`` from now), with the shape of ``x``.""" + t = np.asarray(x, dtype=float) + out = np.asarray(getattr(self.prediction, name)(self.age + t), float) + return out.reshape(t.shape)[()] + + def sf(self, x: ArrayLike) -> Any: + """The probability that the unit survives a further ``x``: + ``pred.sf(age + x)``. + + Parameters + ---------- + x : array_like + The time(s) from now. + + Returns + ------- + float or numpy.ndarray + With the shape of ``x``. + """ + return self._at("sf", x) + + def ff(self, x: ArrayLike) -> Any: + """The probability that the unit fails within a further ``x``: + ``pred.ff(age + x)``, worked out in its own right. + + Parameters + ---------- + x : array_like + The time(s) from now. + """ + return self._at("ff", x) + + def Hf(self, x: ArrayLike) -> Any: + """The cumulative hazard over the next ``x``, ``-log sf(x)``. + + Parameters + ---------- + x : array_like + The time(s) from now. + """ + return self._at("Hf", x) + + def hf(self, x: ArrayLike) -> Any: + """The unit's hazard ``x`` from now, given its readings. + + Parameters + ---------- + x : array_like + The time(s) from now. + """ + return self._at("hf", x) + + def df(self, x: ArrayLike) -> Any: + """The density of the unit's life left at ``x`` from now. + + Parameters + ---------- + x : array_like + The time(s) from now. + """ + return self._at("df", x) + + def qf(self, p: ArrayLike) -> Any: + """The life left that the unit outlives with probability ``1 - + p``: ``pred.qf(p) - age`` (``inf`` where the survival never falls + to ``1 - p``). + + Parameters + ---------- + p : array_like + The probabilities of failing. + """ + q = np.asarray(p, dtype=float) + out = np.asarray(self.prediction.qf(q), dtype=float) - self.age + return out.reshape(q.shape)[()] + + def _knots(self) -> np.ndarray: + """Times marking the curve's scale, for its integral: its + quantiles.""" + from ._point_availability import knots + + return knots(self) + + def mean(self) -> float: + """The unit's mean residual life: the integral of ``sf`` from now, + to a relative 1e-10 on pieces split at its quantiles; ``inf`` where + its survival stays above 0 for ever (a marker path that falls + away can leave a chance of never failing). + + Returns + ------- + float + The mean life left. + """ + from ._mean_lifetime import mean_lifetime + + return float(mean_lifetime(self.sf, self._knots())) + + def _draw(self, u: np.ndarray) -> np.ndarray: + """The lives left ``qf(u)`` of the uniforms ``u``.""" + return np.asarray(self.qf(np.asarray(u, dtype=float)), dtype=float) + + def random(self, size: int) -> np.ndarray: + """Draw ``size`` lives left, from now: ``pred.qf(u) - age`` of + uniforms ``u`` from numpy's global RNG, as ``pred.random`` draws + failure ages. Seed it (``np.random.seed``) to reproduce the draws; + an RBD's seeded methods do. + + Parameters + ---------- + size : int + The number of draws. + + Returns + ------- + numpy.ndarray + The ``size`` lives left, shape ``(size,)``. + """ + return self._draw(np.random.uniform(size=size)) + + def _row_sampler(self) -> RowSampler: + """``random(1)`` as a :class:`~._sampling.RowSampler` (one uniform + per draw), so an RBD with this node batches its draws.""" + return RowSampler( + 1, lambda u: self._draw(np.ascontiguousarray(u[:, 0])) + ) + + # -- Serialisation ---------------------------------------------------- + + def to_dict(self) -> dict: + """The node as a JSON-friendly dict: the joint model's own + ``to_dict()`` and the unit's inputs, its covariates (in the order + of the model's ``feature_names``), its readings as ``[age, + value]`` pairs and its age. ``from_dict`` rebuilds it, and an + RBD's ``to_dict``/``to_json`` save a joint node so. A joint model + on a path model of one's own does not save (surpyval raises + ``NotImplementedError``). + + Returns + ------- + dict + ``{"model": ..., "covariates": [...], "markers": [[age, + value], ...], "age": ...}``. + """ + pred = self.prediction + return { + "model": self.model.to_dict(), + "covariates": [float(z) for z in np.ravel(pred.covariates)], + "markers": [[float(a), float(v)] for a, v in pred.markers.items()], + "age": self.age, + } + + @classmethod + def from_dict(cls, d: dict) -> "JointNode": + """Rebuild a node from :meth:`to_dict`'s dict. + + Parameters + ---------- + d : dict + A dict ``to_dict`` made, e.g. after a JSON round trip. + + Returns + ------- + JointNode + The node, its prediction made again from the saved model and + the unit's inputs. + + Examples + -------- + >>> import json + >>> import numpy as np + >>> from surpyval import ItemLifetime, SurpyvalData, WeibullJointPH + >>> from repyability import JointNode + >>> rng = np.random.default_rng(3) + >>> items = [] + >>> for k in range(60): + ... start, rate = rng.normal(2.0, 0.4), rng.normal(0.02, 0.005) + ... life = 100 * rng.weibull(2.0) * np.exp(-0.9 * (start - 2)) + ... end = min(life, 120.0) + ... ages = np.arange(0.0, end, 10.0) + ... y = start + rate * ages + rng.normal(0, 0.2, ages.size) + ... outcome = "observed" if life < 120 else "right_censored" + ... items.append( + ... ItemLifetime(k, **{outcome: end}, + ... markers=dict(zip(ages, y))) + ... ) + >>> model = WeibullJointPH.fit(SurpyvalData.from_items(items)) + >>> node = JointNode(model, markers={0: 2.1, 40: 3.0}, age=40) + >>> back = JointNode.from_dict(json.loads(json.dumps(node.to_dict()))) + >>> bool(back.sf(20.0) == node.sf(20.0)) + True + """ + import surpyval + + model = surpyval.from_dict(d["model"]) + covariates = list(d.get("covariates") or []) + markers = {float(a): float(v) for a, v in d.get("markers") or []} + return cls( + model, + Z=covariates if covariates else None, + markers=markers or None, + age=float(d["age"]), + ) + + def __repr__(self) -> str: + return ( + f"JointNode(age={self.age:g}, readings=" + f"{len(self.prediction.markers)})" + ) diff --git a/repyability/rbd/non_repairable_rbd.py b/repyability/rbd/non_repairable_rbd.py index 09fecd4..4edd971 100644 --- a/repyability/rbd/non_repairable_rbd.py +++ b/repyability/rbd/non_repairable_rbd.py @@ -83,6 +83,7 @@ PerfectUnreliability, perfect_class, ) +from .joint_node import JointNode, as_node, is_joint from .load_sharing_node import LoadSharingModel from .modular import GraphStructure from .node_state import NodeState @@ -361,12 +362,37 @@ def at_letter(value, name=name): return out +def _check_joint_state(node, model, state: NodeState) -> None: + """Refuse a state that gives a joint node an age other than its + prediction's: its reliability is from that age already, and a state's + age would condition it again.""" + if not isinstance(model, JointNode) or not state.alive: + return + if not math.isclose(state.age, model.age, rel_tol=1e-12, abs_tol=0.0): + raise ValueError( + f"Node {node!r} is a joint model's unit, at age {model.age:g} " + f"now (its prediction's age), and its reliability is from " + f"then: a state of age {state.age:g} would condition it again. " + "Leave it out of the state (or give it its prediction's age, " + "or alive=False); for the unit at another age, predict again " + "with given= that age and its readings since." + ) + + def _check_model(node, model) -> None: """Raise if a node's model is not one (it has no ``sf``, or it is a surpyval distribution itself, whose ``sf`` needs parameters, #233), saying what to give instead; a number most likely means a probability.""" no_distribution(model, f"The model of node {node!r}") + if is_joint(model) and not isinstance(model, JointNode): + raise TypeError( + f"The model of node {node!r} is a surpyval JointModel, whose " + "survival is of a unit given its covariates, readings and age: " + "give the unit's prediction (joint_model.predict(Z, markers, " + "given)), or JointNode(joint_model, Z=..., markers=..., " + "age=...)." + ) if is_degradation(model): check_life(model, f"The model of node {node!r}") if callable(getattr(model, "sf", None)): @@ -433,6 +459,9 @@ class NonRepairableRBD(RBD): [`RepeatedStandbyNode`][repyability.RepeatedStandbyNode], [`RegressionNode`][repyability.RegressionNode] or another ``NonRepairableRBD`` nested as a single node; + - one unit of a fitted surpyval joint model (readings and failures), + given its readings: its ``DynamicPrediction`` (``predict``), or a + [`JointNode`][repyability.JointNode], whose reliability is from now; - [`PerfectReliability`][repyability.PerfectReliability] or [`PerfectUnreliability`][repyability.PerfectUnreliability]. @@ -624,8 +653,10 @@ def __init__( edges = list(edges) ccf_groups = as_groups(ccf_groups) # PerfectReliability() stands for the class (#232). + # A joint model's prediction is the JointNode of its unit. reliabilities = { - node: perfect_class(model) for node, model in reliabilities.items() + node: perfect_class(as_node(model)) + for node, model in reliabilities.items() } self._init_args = { "edges": [tuple(e) for e in edges], @@ -7012,6 +7043,7 @@ def _validate_state(self, state) -> dict: f"distribution components in this release; node {node!r} " f"is a {type(self.reliabilities[node]).__name__}." ) + _check_joint_state(node, self.reliabilities[node], node_state) return state def _state_node_probabilities(self, x, state) -> Dict[Any, ArrayLike]: @@ -7033,6 +7065,10 @@ def _state_node_probabilities(self, x, state) -> Dict[Any, ArrayLike]: node_probabilities[node_name] = np.atleast_1d( PerfectUnreliability.sf(x) ) + elif isinstance(model, JointNode): + # Its reliability is from now already, at its age (which + # the state repeats, _check_joint_state). + node_probabilities[node_name] = np.atleast_1d(model.sf(x)) else: node_probabilities[node_name] = np.atleast_1d( conditional_survival(model, x, node_state.age) diff --git a/repyability/rbd/routes.py b/repyability/rbd/routes.py index 8aa28f2..898d84e 100644 --- a/repyability/rbd/routes.py +++ b/repyability/rbd/routes.py @@ -157,12 +157,20 @@ def refusal(check: Callable[[], object]) -> Optional[str]: return None +#: How a joint node's reliability is worked out. +JOINT_HOW = ( + "its joint model's prediction, by quadrature over the unit's marker " + "path" +) + + def model_route(model) -> Tuple[str, str]: """How a node model's reliability (its ``sf``) is obtained: the route, and a phrase saying how.""" from surpyval import Parametric from .helper_classes import PerfectReliability, PerfectUnreliability + from .joint_node import JointNode from .load_sharing_node import LoadSharingModel from .non_repairable_rbd import NonRepairableRBD from .numerical_convolution import ConvolvedSurvival @@ -200,6 +208,8 @@ def model_route(model) -> Tuple[str, str]: return worst, f"a nested RBD with {worst} nodes {_names(which)}" if isinstance(model, Parametric): return EXACT, "a closed form" + if isinstance(model, JointNode): + return NUMERICAL, JOINT_HOW return EXACT, "its own sf" @@ -207,6 +217,7 @@ def mean_route(model) -> Tuple[str, str]: """How a node model's mean lifetime (its ``mean()``) is found: the route, and a phrase saying how, or the message its ``mean()`` would raise.""" from .degrading_node import DegradingNode + from .joint_node import JointNode from .load_sharing_node import LoadSharingModel from .non_repairable_rbd import NonRepairableRBD from .regression_node import RegressionNode @@ -239,6 +250,8 @@ def mean_route(model) -> Tuple[str, str]: return model_route(model) if isinstance(model, RegressionNode): return NUMERICAL, "the mean of its survival function on a grid" + if isinstance(model, JointNode): + return NUMERICAL, "the area under its prediction's survival" return EXACT, "its model's mean" diff --git a/repyability/rbd/serialisation.py b/repyability/rbd/serialisation.py index a2781bf..e612fd0 100644 --- a/repyability/rbd/serialisation.py +++ b/repyability/rbd/serialisation.py @@ -38,6 +38,7 @@ PerfectReliability, PerfectUnreliability, ) +from repyability.rbd.joint_node import JointNode, as_node from repyability.rbd.load_sharing_node import LoadSharingModel from repyability.rbd.rbd import RBD from repyability.rbd.regression_node import RegressionNode @@ -57,6 +58,7 @@ "repeated_node": RepeatedNode, "non_repairable": NonRepairable, "regression_node": RegressionNode, + "joint_node": JointNode, "load_sharing": LoadSharingModel, } @@ -158,6 +160,10 @@ def _serialise_model(model: Any) -> dict: } if isinstance(model, RegressionNode): return {"kind": "regression_node", **model.to_dict()} + if isinstance(as_node(model), JointNode): + # A joint model's unit (or its prediction, which loads as one): + # the model in surpyval's format and the unit's inputs. + return {"kind": "joint_node", **as_node(model).to_dict()} if isinstance(model, LoadSharingModel): return { "kind": "load_sharing", @@ -178,7 +184,7 @@ def _serialise_model(model: Any) -> dict: f"Cannot serialise a node model of type {type(model).__name__}. " "Only surpyval's parametric models, the " "RePyability node wrappers (standby, degrading, repeated, " - "NonRepairable, load-sharing, regression), perfect " + "NonRepairable, load-sharing, regression, joint), perfect " "reliability/unreliability and nested RBDs are supported." ) @@ -238,6 +244,8 @@ def deserialise_model(d: dict) -> Any: ) if kind == "regression_node": return RegressionNode.from_dict(d) + if kind == "joint_node": + return JointNode.from_dict(d) if kind == "load_sharing": return LoadSharingModel( [surpyval.from_dict(md) for md in d["models"]], diff --git a/repyability/tests/capabilities.py b/repyability/tests/capabilities.py index b1c33f0..5cf0695 100644 --- a/repyability/tests/capabilities.py +++ b/repyability/tests/capabilities.py @@ -344,6 +344,7 @@ ("A probability per demand", "fixed"), ("A mixture (`MixtureModel`)", "mixture"), ("A regression node (covariates)", "regression"), + ("A joint model's unit, given its readings", "joint"), ("A degrading node", "degrading"), ("A fitted degradation process", "degradation process"), ("Cold standby (convolved)", "convolved standby"), diff --git a/repyability/tests/catalogue.py b/repyability/tests/catalogue.py index c2f8f7f..83bf3dd 100644 --- a/repyability/tests/catalogue.py +++ b/repyability/tests/catalogue.py @@ -22,6 +22,7 @@ CCFGroup, DegradingNode, FaultTree, + JointNode, LoadSharingModel, NonRepairable, NonRepairableRBD, @@ -103,6 +104,45 @@ def aft(): return surv.ExponentialAFT.fit(x, Z=load) +@_fitted_once +def joint_model(): + """A Weibull joint model (surpyval's WeibullJointPH) of a fleet of 60 + pumps, each with a marker read every 10 hours and a load: a pump's + marker starts near 2 and rises by about 0.02 an hour, and a higher + marker or load shortens its life.""" + from surpyval import ItemLifetime, SurpyvalData, WeibullJointPH + + rng = np.random.default_rng(3) + items = [] + for k in range(60): + start, rate = rng.normal(2.0, 0.4), rng.normal(0.02, 0.005) + load = float(k % 2) + life = 100 * rng.weibull(2.0) * np.exp(-0.9 * (start - 2) - 0.3 * load) + end = min(life, 120.0) + ages = np.arange(0.0, end, 10.0) + y = start + rate * ages + rng.normal(0, 0.2, ages.size) + outcome = "observed" if life < 120 else "right_censored" + items.append( + ItemLifetime( + k, + **{outcome: end}, + markers=dict(zip(ages, y)), + covariates={"load": load}, + ) + ) + return WeibullJointPH.fit(SurpyvalData.from_items(items)) + + +#: A pump of the joint model's fleet, 40 hours old, at load 1, whose +#: marker has risen faster than most. +PUMP = {"Z": {"load": 1.0}, "markers": {0: 2.1, 20: 2.5, 40: 3.0}, "age": 40} + + +def joint_unit(): + """The pump ``PUMP`` of ``joint_model``, as a JointNode.""" + return JointNode(joint_model(), **PUMP) + + @_fitted_once def fitted_life(): """A Weibull fitted to 30 failures: a fit with a parameter covariance.""" @@ -268,6 +308,7 @@ def nonrepairable_kinds(): "regression": NonRepairableRBD( EDGES, {"a": RegressionNode(aft(), covariates=[1.0]), **rest} ), + "joint": NonRepairableRBD(EDGES, {"a": joint_unit(), **rest}), "fitted": NonRepairableRBD(EDGES, {"a": fitted_life(), **rest}), "fitted shares": NonRepairableRBD(EDGES, {"a": shares_life(), **rest}), "degradation process": NonRepairableRBD(EDGES, {"a": wear(), **rest}), diff --git a/repyability/tests/test_catalogue.py b/repyability/tests/test_catalogue.py index d60871b..ea04502 100644 --- a/repyability/tests/test_catalogue.py +++ b/repyability/tests/test_catalogue.py @@ -24,6 +24,7 @@ import repyability from repyability import ( DegradingNode, + JointNode, LoadSharingModel, NonRepairable, NonRepairableRBD, @@ -83,6 +84,7 @@ RepeatedNode, RepeatedStandbyNode, StandbyModel, + JointNode, } #: The node classes a RepairableRBD takes, as a component or its life. REPAIRABLE_NODES = { @@ -159,6 +161,8 @@ def model_kinds(model) -> set: return {cls.__name__} if isinstance(model, NonRepairable): return {"NonRepairable"} | model_kinds(model.reliability) + if isinstance(model, JointNode): + return {"JointNode"} if is_mixture(model): return {"surpyval MixtureModel"} if is_degradation(model): diff --git a/repyability/tests/test_joint_node.py b/repyability/tests/test_joint_node.py new file mode 100644 index 0000000..19d23d9 --- /dev/null +++ b/repyability/tests/test_joint_node.py @@ -0,0 +1,215 @@ +"""A unit of a surpyval joint model (readings and failures) as a block of +a non-repairable diagram: ``JointNode``. + +The node is the unit's dynamic prediction moved to now: its reliability at +mission time ``t`` is ``pred.sf(age + t)`` (surpyval has conditioned it on +the unit reaching ``age``), its lifetimes ``pred.qf(u) - age``. The checks +are against those values worked out by hand, and the simulations against +them within their sampling error.""" + +import json + +import numpy as np +import pytest +import surpyval as surv + +from repyability import ( + JointNode, + NodeState, + NonRepairable, + NonRepairableRBD, + RepairableRBD, +) +from repyability.rbd import routes +from repyability.rbd.serialisation import exactly, serialise_model +from repyability.tests.catalogue import PUMP, joint_model, joint_unit + +SERIES = [("in", "pump"), ("pump", "valve"), ("valve", "out")] +PARALLEL = [("in", "pump"), ("in", "valve"), ("pump", "out"), ("valve", "out")] +VALVE = surv.Weibull.from_params([80.0, 1.5]) +T = np.array([5.0, 20.0, 45.0, 90.0]) + + +@pytest.fixture(scope="module") +def pred(): + """The pump's prediction: 40 hours old, at load 1, readings rising.""" + model = joint_model() + return model.predict(PUMP["Z"], PUMP["markers"], given=PUMP["age"]) + + +def test_the_node_is_the_prediction_from_now(pred): + node = JointNode(pred) + assert node.age == 40.0 + np.testing.assert_array_equal(node.sf(T), pred.sf(40.0 + T)) + np.testing.assert_array_equal(node.ff(T), pred.ff(40.0 + T)) + np.testing.assert_array_equal(node.Hf(T), pred.Hf(40.0 + T)) + np.testing.assert_array_equal(node.df(T), pred.df(40.0 + T)) + np.testing.assert_array_equal(node.hf(T), pred.hf(40.0 + T)) + # Now it is surely running: the prediction is conditioned on its age. + assert node.sf(0.0) == 1.0 + assert np.ndim(node.sf(20.0)) == 0 + p = np.array([0.1, 0.5, 0.9]) + np.testing.assert_array_equal(node.qf(p), pred.qf(p) - 40.0) + np.testing.assert_allclose(node.sf(node.qf(p)), 1.0 - p, rtol=1e-8) + + +def test_a_model_with_the_inputs_is_its_prediction(pred): + node = joint_unit() + np.testing.assert_array_equal(node.sf(T), JointNode(pred).sf(T)) + assert node.prediction.markers == pred.markers + + +def test_series_and_parallel_are_the_products(pred): + pump = pred.sf(40.0 + T) + valve = VALVE.sf(T) + series = NonRepairableRBD(SERIES, {"pump": pred, "valve": VALVE}) + parallel = NonRepairableRBD(PARALLEL, {"pump": pred, "valve": VALVE}) + # A prediction given as a node is its unit's JointNode. + assert isinstance(series.reliabilities["pump"], JointNode) + np.testing.assert_allclose(series.sf(T), pump * valve, rtol=1e-12) + np.testing.assert_allclose( + parallel.ff(T), (1.0 - pump) * (1.0 - valve), rtol=1e-12 + ) + # Birnbaum: the system's reliability with the pump working less with + # it failed, in series the valve's reliability. + np.testing.assert_allclose( + series.birnbaum_importance(20.0)["pump"], VALVE.sf(20.0), rtol=1e-12 + ) + + +def test_its_draws_are_its_life_left(pred): + node = JointNode(pred) + np.random.seed(11) + drawn = node.random(500) + np.random.seed(11) + u = np.random.uniform(size=500) + np.testing.assert_array_equal(drawn, pred.qf(u) - 40.0) + assert np.all(drawn >= 0.0) + + +def test_the_simulations_agree_with_the_exact_values(pred): + rbd = NonRepairableRBD(SERIES, {"pump": pred, "valve": VALVE}) + n = 2000 + lives = rbd.random(n, seed=3) + exact = rbd.sf(T) + simulated = (lives[:, None] > T[None, :]).mean(axis=0) + error = np.sqrt(exact * (1.0 - exact) / n) + assert np.all(np.abs(simulated - exact) <= 4.0 * error + 1e-12) + # The MTTF from now: exact against simulated. + mean = rbd.mean() + assert np.isfinite(mean) + assert abs(lives.mean() - mean) <= 4.0 * lives.std() / np.sqrt(n) + + +def test_its_mean_is_the_mean_life_left(pred): + # A straight marker path whose slope varies between units leaves the + # unit a chance of never failing (a path that falls away), however + # small: its mean life left is then infinite, as the model says. + node = JointNode(pred) + assert pred.sf(1e6) > 1e-14 + assert node.mean() == np.inf + # In series with a valve, whose life ends, the system's is finite: the + # area under the products. + rbd = NonRepairableRBD(SERIES, {"pump": node, "valve": VALVE}) + from scipy.integrate import simpson + + t = np.linspace(0.0, 800.0, 8001) + area = simpson(rbd.sf(t), x=t) + assert rbd.mean() == pytest.approx(area, rel=1e-7) + + +def test_a_state_takes_its_age_from_the_prediction(pred): + rbd = NonRepairableRBD(SERIES, {"pump": pred, "valve": VALVE}) + valve = VALVE.sf(30.0 + T) / VALVE.sf(30.0) + expected = pred.sf(40.0 + T) * valve + left_out = rbd.sf_given_state(T, {"valve": NodeState(age=30.0)}) + np.testing.assert_allclose(left_out, expected, rtol=1e-12) + # Its own age again is the same, not conditioned twice. + same = rbd.sf_given_state( + T, {"pump": NodeState(age=40.0), "valve": NodeState(age=30.0)} + ) + np.testing.assert_allclose(same, expected, rtol=1e-12) + failed = rbd.sf_given_state(T, {"pump": NodeState(alive=False)}) + np.testing.assert_array_equal(failed, 0.0) + for age in (0.0, 10.0, 60.0): + with pytest.raises(ValueError, match="condition it again"): + rbd.sf_given_state(T, {"pump": NodeState(age=age)}) + with pytest.raises(ValueError, match="condition it again"): + rbd.importances_given_state(20.0, {"pump": NodeState(age=5.0)}) + + +def test_a_repairable_component_takes_no_joint_model(pred): + repairs = surv.Exponential.from_params([1.0]) + for life in (pred, JointNode(pred)): + with pytest.raises(ValueError, match="takes no joint model"): + NonRepairable(life, repairs) + with pytest.raises(ValueError, match="takes no joint model"): + RepairableRBD( + SERIES, + { + "pump": {"reliability": life, "repairability": repairs}, + "valve": {"reliability": VALVE, "repairability": repairs}, + }, + ) + with pytest.raises(TypeError, match="takes no joint model"): + RepairableRBD( + SERIES, + { + "pump": life, + "valve": {"reliability": VALVE, "repairability": repairs}, + }, + ) + + +def test_a_diagram_with_a_joint_node_round_trips(pred): + for given in (pred, JointNode(pred)): + rbd = NonRepairableRBD(SERIES, {"pump": given, "valve": VALVE}) + text = rbd.to_json() + back = NonRepairableRBD.from_json(text) + node = back.reliabilities["pump"] + assert isinstance(node, JointNode) + assert node.age == 40.0 + assert node.prediction.markers == pred.markers + np.testing.assert_allclose(back.sf(T), rbd.sf(T), rtol=1e-12) + # The saved form is plain JSON, and loads as the class it is. + with exactly(): + saved = serialise_model(rbd.reliabilities["pump"]) + assert saved["kind"] == "joint_node" + assert json.loads(json.dumps(saved)) == saved + assert saved["covariates"] == [1.0] + + +def test_its_inputs_are_checked(pred): + model = joint_model() + with pytest.raises(TypeError, match="leave out Z"): + JointNode(pred, Z={"load": 1.0}) + with pytest.raises(ValueError, match="age now"): + JointNode(model, Z={"load": 1.0}, markers={0: 2.0}) + with pytest.raises(ValueError, match="missing"): + JointNode(model, Z={"load": np.nan}) + with pytest.raises(TypeError, match="joint model's prediction"): + JointNode(VALVE) + with pytest.raises(TypeError, match="JointNode"): + NonRepairableRBD(SERIES, {"pump": model, "valve": VALVE}) + + +def test_a_new_unit_is_the_population_life(): + # No readings and age 0: a new unit at its load, the model's own + # survival over the population's marker paths. + model = joint_model() + node = JointNode(model, Z={"load": 1.0}) + assert node.age == 0.0 + np.testing.assert_allclose( + node.sf(T), model.sf(T, Z={"load": 1.0}), rtol=1e-12 + ) + + +def test_the_routes_say_it_is_numerical(pred): + route, how = routes.model_route(JointNode(pred)) + assert route == routes.NUMERICAL + assert "quadrature" in how + report = NonRepairableRBD( + SERIES, {"pump": pred, "valve": VALVE} + ).analysis_routes() + assert report["sf"].route == routes.NUMERICAL + assert "pump" in report["sf"].reason