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