diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml index 192ca33..32b4b33 100644 --- a/.github/workflows/actions.yml +++ b/.github/workflows/actions.yml @@ -2,12 +2,12 @@ name: actions run-name: RePyability CI -# The one automated test run: the full suite on a pull request into master +# The one automated test run: the full suite on a pull request into main # (from dev). Nothing else runs it; release.yml publishes only what it # passed. on: pull_request: - branches: [master] + branches: [main] jobs: lint: diff --git a/.github/workflows/close-issues.yml b/.github/workflows/close-issues.yml index fcdfbcc..ac2d3d7 100644 --- a/.github/workflows/close-issues.yml +++ b/.github/workflows/close-issues.yml @@ -2,10 +2,10 @@ name: close issues # Close the issues a pull request into dev finishes, when it merges. # GitHub itself closes an issue only when the pull request that says -# "Closes #N" merges into the default branch (master), which here happens +# "Closes #N" merges into the default branch (main), which here happens # at a release: this closes them as the work reaches dev, with a comment # naming the pull request. It runs no tests (the one automated test run -# is actions.yml, on the pull request from dev into master). +# is actions.yml, on the pull request from dev into main). on: pull_request: diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 1ed7cf2..c0e84d1 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -1,7 +1,7 @@ name: docs # Build the MkDocs site and publish it to GitHub Pages on every push to -# master. Pull requests are validated by the strict docs build in +# main. Pull requests are validated by the strict docs build in # actions.yml; this workflow does the actual deployment. # # One-time repository setup (Settings -> Pages -> Build and deployment): @@ -9,7 +9,7 @@ name: docs on: push: - branches: [master] + branches: [main] workflow_dispatch: # Grant the GITHUB_TOKEN the permissions the Pages deployment needs. diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index d564615..41c7f17 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -2,12 +2,12 @@ name: release # Publishes the package to PyPI, in either of two ways. # -# 1. Run this workflow by hand on master (Actions -> release -> Run workflow; -# `gh workflow run release.yml --ref master -f version=X.Y`; or the REST +# 1. Run this workflow by hand on main (Actions -> release -> Run workflow; +# `gh workflow run release.yml --ref main -f version=X.Y`; or the REST # API). It checks that the version typed in matches # repyability/_version.py, that CHANGELOG.md has a dated section for it, # that the tag vX.Y does not exist yet, that the full suite passed on -# these files (on the pull request from dev into master) and +# these files (on the pull request from dev into main) and # that PyPI does not have the version. It then builds, publishes to PyPI, # and last creates the tag vX.Y and a GitHub Release. The release notes # open with the CHANGELOG section's summary, followed by GitHub's @@ -71,8 +71,8 @@ jobs: echo "::error::'$REQUESTED' isn't major.minor: versions have two parts (e.g. 0.11)." exit 1 fi - if [ "$GITHUB_REF" != "refs/heads/master" ]; then - echo "::error::Releases are cut from master, not $GITHUB_REF." + if [ "$GITHUB_REF" != "refs/heads/main" ]; then + echo "::error::Releases are cut from main, not $GITHUB_REF." exit 1 fi version=$(sed -n 's/^__version__ = "\(.*\)"$/\1/p' repyability/_version.py) @@ -88,9 +88,9 @@ jobs: echo "::error::The tag v$version already exists." exit 1 fi - # The full suite runs on the pull request from dev into master + # The full suite runs on the pull request from dev into main # (actions.yml), whose run is recorded against dev's commit. The - # merge commit on master holds the same files, so a pass on this + # merge commit on main holds the same files, so a pass on this # commit, or on a parent with the same files, counts. repo="repos/$GITHUB_REPOSITORY" tree=$(gh api "$repo/git/commits/$GITHUB_SHA" --jq '.tree.sha') @@ -104,7 +104,7 @@ jobs: fi done if [ -z "$passed" ]; then - echo "::error::The full suite (actions.yml) has not passed on these files: on $GITHUB_SHA or the dev commit it merged. Open the pull request from dev into master and wait for it." + echo "::error::The full suite (actions.yml) has not passed on these files: on $GITHUB_SHA or the dev commit it merged. Open the pull request from dev into main and wait for it." exit 1 fi echo "The full suite passed on $passed, with the same files as $GITHUB_SHA." diff --git a/CHANGELOG.md b/CHANGELOG.md index 39e74b5..abf071f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,42 @@ other release, fixes included, the minor. ## [Unreleased] +## [0.15] - 2026-10-17 + +surpyval 0.25, and a joint model in a diagram. A unit of a surpyval joint +model of readings and failures (`surpyval.WeibullJointPH`), given its +readings so far, is a block of a `NonRepairableRBD` (`JointNode`, or its +`DynamicPrediction` given as the node): its reliability is from now, and +it takes part in every analysis of the diagram. A `RepairableRBD` refuses +it, as the units that replace it are new and without its readings. What +0.14 deprecated is gone. + +Behaviour changes: surpyval 0.25 is required; `NodeState`'s `phase` is +refused (it is `since_scheduled`). + +### 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). + +### Removed + +- `NodeState`'s `phase`, deprecated in 0.14 and renamed `since_scheduled` + (#350): giving it raises `TypeError`, and a state has no `phase`. + ## [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 1b60227..8aba763 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -1,5 +1,14 @@ # RePyability — project notes for Claude +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.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 - **Distribution/lifetime fitting stays in surpyval, not RePyability.** surpyval @@ -9,8 +18,8 @@ "from data to system reliability" workflow is wanted, do the fitting in surpyval and pass the resulting models in — do not add fitting logic here. -- **Visualization stays out of RePyability.** The Reliafy app (separate, - open source) handles plotting/dashboards/reporting. RePyability is the +- **Visualization stays out of RePyability.** The Reliafy app (separate) + handles plotting/dashboards/reporting. RePyability is the computational reliability engine; keep it free of plotting dependencies. - **RePyability is about timelines; Reliafy is about calendars.** @@ -38,12 +47,12 @@ - **surpyval releases reach RePyability users at once.** surpyval is required with no upper bound, so its next release is what a fresh install - gets. No automated run tests surpyval's `develop` (the maintainer wants - one automated test run, below): before surpyval releases, run - RePyability's tests against its `develop` by hand (`pip install - "surpyval @ git+https://github.com/derrynknife/SurPyval@develop"`), and + gets. No automated run tests surpyval's `dev` (there is one automated + test run, on the release pull request): before surpyval releases, run + RePyability's tests against its `dev` by hand (`pip install + "surpyval @ git+https://github.com/derrynknife/SurPyval@dev"`), and when they fail, fix RePyability (working with both the released surpyval - and `develop`) and release that first. CI's `test (minimum dependencies)` + and `dev`) and release that first. CI's `test (minimum dependencies)` job tests the oldest surpyval `pyproject.toml` allows (and every other dependency at its floor); raise that minimum, rather than keep code for older versions, once RePyability needs what a @@ -423,38 +432,24 @@ that `levers()` lists what `parameter_sensitivity` reports, in its order, that each value lies in its range, and that `with_levers` refuses a value past it. -- **A deprecation gives one minor release's notice.** It warns in one - minor release and the next removes it, with a `FutureWarning` (always - shown) through `repyability/utils/deprecation.py`. What 0.11 deprecated - went in 0.12 (#149), and `test_removed_in_0_12.py` keeps it gone: the - old names and ignored arguments raise `TypeError`, a diagram refuses a - non-parametric node, and a standby or load-sharing model with no exact - or numerical reliability refuses one (`is_simulated`), where a fit to - simulated lifetimes stood in, and is left to the simulations. What 0.12 - deprecated went in 0.13, and `test_removed_in_0_13.py` keeps it gone: - calling `SparesDemand.mean()`/`std()` (#184) and those models' - `mc_samples`, `lower` and `seed` raise `TypeError`, and `StandbyModel`'s - options after them are given by name. What 0.13 deprecated went in - 0.14, and `test_removed_in_0_14.py` keeps it gone: - `optimal_inspection_intervals(offsets=)` (#222) raises `TypeError`, for - `offset_shares=`; `CapacityDistribution.mean` is a plain number or - array, not callable (#235); the exact `mean` of a `StandbyModel`, - `LoadSharingModel` or `DegradingNode` refuses `mc_samples` and `seed` - (#233, `standby_node.drawn_mean`); and `rbd._sensitivity`'s `levers` - and `Lever` are gone (#244). What 0.14 deprecates goes in 0.15 - (`NEXT_REMOVAL`): `NodeState`'s `phase`, renamed `since_scheduled` - (#350, through `deprecation.renamed` for the argument and - `renamed_attribute` for reading it); a new deprecation warns through - those, `called` or `ignored`. `test_what_0_14_deprecates_goes_in_0_15` - fails once the version reaches it. +- **Deprecations** follow CONTRIBUTING (one minor release's notice, a + `FutureWarning` through `repyability/utils/deprecation.py`). Each + release's removals stay gone through `test_removed_in_0_12.py`, + `_0_13.py`, `_0_14.py` and `_0_15.py` (`NodeState`'s `phase`, #350): + a removed argument raises `TypeError`. 0.15 deprecates nothing; what it + would deprecate goes in 0.16 (`NEXT_REMOVAL`). A new deprecation warns + through `deprecation.renamed`, `renamed_attribute`, `called` or + `ignored`, with a test that fails once the version reaches its removal + (as `test_what_0_14_deprecates_goes_in_0_15` did). - **No result depends on the unit of time (or of a covariate).** A step, grid or probe is sized from the models (a typical life, `failure_time_scale`), never in absolute time, and a proportional change of a lever is of its scale (`_sensitivity._proportional_scale`: 1 for a log-scale parameter, none for a covariate). `test_units.py` checks results scale with the unit, from 1e-6 to 1e6. A model's unit - is surpyval's to carry (SurPyval#845); until it does, every node is - taken to be in one unit, unchecked. + is surpyval's to carry: its next release labels a model with one + (`model.unit`, SurPyval#845, not converted or checked). Until RePyability + reads it, every node is taken to be in one unit, unchecked. - **Exact by default, simulation on request.** Where an analysis can be computed exactly or numerically, that is the default, and the Monte-Carlo estimate is a `method="simulate"` away (as for `NonRepairableRBD.mean`). @@ -471,48 +466,24 @@ them (a renewal sum, a band of a grid) then runs hundreds of times slower. Matrix products are left to BLAS. -## Testing - -- **One automated test run: the full suite on the pull request from dev - into master.** `actions.yml` runs on nothing else: not on pushes, not on - pull requests into dev, not on a schedule, and there is no other test - workflow. Work going into dev is checked locally by the tests of what it - changes (the files that exercise the code touched, and new tests for - it), and is merged without waiting for CI. Do not start, re-run or add - other automated test runs; run the full suite locally only when the - maintainer asks. - -## Releasing - -Releases are cut from master by `.github/workflows/release.yml`, which this -session can run: it cannot push tags or create GitHub Releases itself. Every -merge and release needs the maintainer's go-ahead. - -Versions have two parts, major.minor, from 0.11 (they had three until -0.10.1): from 1.0, a release that breaks compatibility raises the major -number, and any other release, fixes included, the minor. There are no -patch releases, and `release.yml` refuses a version that isn't -major.minor. - -1. On the working branch, bump `repyability/_version.py`. Roll the CHANGELOG - `[Unreleased]` section into `## [X.Y] - YYYY-MM-DD`, opening with a - summary paragraph and the behaviour changes (they open the release notes), - and update the version in `docs/guide/saving.md`. PR to dev, then dev to - master. A pull request into dev lists "Closes #N" for each issue it - finishes, and `.github/workflows/close-issues.yml` closes them when it - merges, with a comment naming it (GitHub itself closes an issue only on a - merge into master): commit messages' "(#N)" close nothing, and the - dev-to-master pull request need not list them again. -2. When the full suite has passed on the pull request from dev into - master, merge it (a merge commit) and run the workflow with - `actions_run_trigger`: `run_workflow`, workflow `release.yml`, ref - `master`, inputs `{"version": "X.Y", "dry_run": "true"}`. If that - passes, run it again with `"dry_run": "false"`. Then check the run, the - tag, the GitHub Release and `https://pypi.org/pypi/repyability/X.Y/json`. +## Testing and releasing + +The checks, the one automated test run and the release steps are in +CONTRIBUTING.md. What it leaves to the maintainer's say: + +- Do not start, re-run or add automated test runs, and run the full suite + locally only when the maintainer asks. Work going into `dev` is checked + locally by the tests of what it changes (`scripts/check.sh `) and + merged without waiting for CI. +- Every merge and every release needs the maintainer's go-ahead. Merge the + release pull request with a merge commit (`release.yml` finds CI's pass + on the merged `dev` commit), run `release.yml` with `dry_run` first, and afterwards check the run, the tag + `vX.Y`, the GitHub Release and `https://pypi.org/pypi/repyability/X.Y/json`. ## 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 @@ -521,8 +492,8 @@ limited-failure proportion `lfp_p` (#608) and takes `success_run`'s - **A mixture's quantile** (surpyval #821). surpyval 0.24's `MixtureModel` has a `qf` (#651), but it inverts `F` in the upper tail, losing the long lives' precision (about 1e-4), and takes some 200 times as long as - RePyability's. `_sampling.mixture_quantile` inverts the distribution - (with its components' `sf` summed above the median), which + RePyability's (fixed on surpyval's `dev`, for its next release). + `_sampling.mixture_quantile` inverts the distribution (with its components' `sf` summed above the median), which `stream_sampler` gives a `RepairableRBD`'s streams (a change to it changes seeded results). Take the mixture's own `qf` once the minimum surpyval's is as precise and fast. @@ -544,11 +515,19 @@ limited-failure proportion `lfp_p` (#608) and takes `success_run`'s change to it changes seeded results). Take surpyval's once the minimum surpyval has it. -- **Rescaling a model's time for a duty cycle** (surpyval #845). - A spec's `"duty"` puts a life in operating time on the calendar, - `R(d t)`, by `_duty._SCALED`, one rule per surpyval distribution (a new +- **Rescaling a model's time for a duty cycle.** surpyval has no + rescaling (#845 gave models only a unit label). A spec's `"duty"` puts a + life in operating time on the calendar, `R(d t)`, by `_duty._SCALED`, one rule per surpyval distribution (a new one needs its rule there, which `test_units` checks against `R(d t)`). Take surpyval's rescaling once the minimum surpyval has it. +- **A fit's offset under parameter uncertainty** (surpyval #830). + surpyval's `covariance()` leaves an offset fit's `gamma` out, so + `rbd/uncertainty.py` holds it where it was fitted, with a warning. Draw it + once the minimum surpyval's covariance includes it. + +Refused until surpyval has it: an ARA, ARI or G1 renewal model as a life, +whose next failure surpyval cannot draw yet (surpyval #833, +`rbd/_processes.py`). List each new workaround here with its surpyval issue and where it lives, so it can go once the minimum surpyval in `pyproject.toml` includes the diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 49fa6fd..b61c942 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -18,15 +18,15 @@ 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 Work branches off **`dev`** and goes back into `dev` by pull request. -Releases go from `dev` into `master`. +Releases go from `dev` into `main`. There is **one automated test run**: the full suite on the pull request from -`dev` into `master` (`.github/workflows/actions.yml`). Nothing runs on pushes +`dev` into `main` (`.github/workflows/actions.yml`). Nothing runs on pushes or on pull requests into `dev`, so check your change locally before it goes in: @@ -104,8 +104,10 @@ go once the minimum surpyval version has the fix. A deprecation warns (a `FutureWarning`, through `repyability/utils/deprecation.py`) for one minor release, and the next -removes it. `test_what_0_14_deprecates_goes_in_0_15` and its kind fail once -the version reaches the removal, as a reminder. +removes it. A test of the kind `test_what_0_14_deprecates_goes_in_0_15` +(in `test_api_consistency.py`) fails once the version reaches the removal, +as a reminder; the removed names then stay gone through +`test_removed_in_X_Y.py`. ## Releases (maintainers) @@ -118,10 +120,10 @@ release: `[Unreleased]` section into a dated `## [X.Y] - YYYY-MM-DD` section. Its opening paragraphs become the summary at the top of the release notes. Update the version in `docs/guide/saving.md`. Merge to `dev`, then open - the pull request from `dev` into `master`, which runs the full suite. -2. Once it has passed, merge it and run the `release` workflow on master + the pull request from `dev` into `main`, which runs the full suite. +2. Once it has passed, merge it and run the `release` workflow on main with the version: Actions → release → Run workflow, or - `gh workflow run release.yml --ref master -f version=X.Y`. Add + `gh workflow run release.yml --ref main -f version=X.Y`. Add `-f dry_run=true` to check and build without publishing. The workflow checks the version, the CHANGELOG section, that the tag is new, that CI passed and that PyPI does not have the version. Then it builds, publishes 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/docs/guide/saving.md b/docs/guide/saving.md index ad20eb8..8c50c0f 100644 --- a/docs/guide/saving.md +++ b/docs/guide/saving.md @@ -34,7 +34,7 @@ path = Path(tempfile.mkdtemp()) / "plant.json" rbd.to_json(path) # written; returns None NonRepairableRBD.from_json(path).sf(30) == rbd.sf(30) # True type(RBD.from_dict(data)).__name__ # 'NonRepairableRBD': the base class dispatches on type -data["type"], data["repyability_version"] # ('NonRepairableRBD', '0.14') +data["type"], data["repyability_version"] # ('NonRepairableRBD', '0.15') ``` What is saved: 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/_version.py b/repyability/_version.py index 1e8b0e3..2d6e9aa 100644 --- a/repyability/_version.py +++ b/repyability/_version.py @@ -5,4 +5,4 @@ dependencies) via ``[tool.setuptools.dynamic]`` in ``pyproject.toml``. """ -__version__ = "0.14" +__version__ = "0.15" 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/node_state.py b/repyability/rbd/node_state.py index fe88f81..7786b27 100644 --- a/repyability/rbd/node_state.py +++ b/repyability/rbd/node_state.py @@ -16,18 +16,8 @@ from dataclasses import dataclass from typing import Optional -from repyability.utils.deprecation import renamed, renamed_attribute -#: The release that refuses ``phase``, renamed ``since_scheduled`` in 0.14 -#: (#350). -_PHASE_REMOVAL = "0.15" -_PHASE_WHY = ( - "it is the time since the last scheduled replacement or test, and a " - "phased mission has phases of its own" -) - - -@dataclass(frozen=True, init=False) +@dataclass(frozen=True) class NodeState: """The current condition of a single RBD node. @@ -93,9 +83,6 @@ class NodeState: ``virtual_age`` now and ``since_failure``: give ``age= since_failure`` and ``virtual_age=virtual_age - since_failure``. The simulations take it. - phase : float, optional - Deprecated: ``since_scheduled``, which it is renamed (#350), as a - phased mission has phases of its own. Refused in 0.15. Raises ------ @@ -104,8 +91,6 @@ class NodeState: is negative or not finite, ``down_for`` or ``maintenance`` is given for a component that is alive, or another field is given with ``stationary``. - TypeError - If ``since_scheduled`` is given under both names. Notes ----- @@ -162,54 +147,6 @@ class NodeState: stationary: bool = False virtual_age: Optional[float] = None - def __init__( - self, - age: float = 0.0, - alive: bool = True, - down_for: float = 0.0, - maintenance: bool = False, - since_scheduled: Optional[float] = None, - stationary: bool = False, - virtual_age: Optional[float] = None, - *, - phase: Optional[float] = None, - ) -> None: - if phase is not None: - if since_scheduled is not None: - raise TypeError( - "NodeState takes since_scheduled alone: phase is its " - "old name." - ) - renamed( - "NodeState", - "phase", - "since_scheduled", - _PHASE_WHY, - removal=_PHASE_REMOVAL, - ) - since_scheduled = phase - object.__setattr__(self, "age", age) - object.__setattr__(self, "alive", alive) - object.__setattr__(self, "down_for", down_for) - object.__setattr__(self, "maintenance", maintenance) - object.__setattr__(self, "since_scheduled", since_scheduled) - object.__setattr__(self, "stationary", stationary) - object.__setattr__(self, "virtual_age", virtual_age) - self.__post_init__() - - @property - def phase(self) -> Optional[float]: - """``since_scheduled``, by its old name: deprecated, and refused in - 0.15 (#350).""" - renamed_attribute( - "NodeState", - "phase", - "since_scheduled", - _PHASE_WHY, - removal=_PHASE_REMOVAL, - ) - return self.since_scheduled - def __post_init__(self) -> None: if self.age < 0: raise ValueError( 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_api_consistency.py b/repyability/tests/test_api_consistency.py index 4df96b5..dd15a18 100644 --- a/repyability/tests/test_api_consistency.py +++ b/repyability/tests/test_api_consistency.py @@ -30,7 +30,6 @@ mtbf_pass_probability, mtbf_test_time, ) -from repyability._version import __version__ from repyability.rbd.ccf import with_parameters from repyability.rbd.uncertainty import draw_ccf_models from repyability.utils import deprecation @@ -261,13 +260,6 @@ def test_the_spares_demand_values_are_properties(): assert demand.std == pytest.approx(np.sqrt(10.0), rel=1e-4) -def test_what_0_14_deprecates_goes_in_0_15(): - # NodeState's phase, renamed since_scheduled (#350): giving phase= and - # reading .phase warn, and NEXT_REMOVAL (0.15) refuses them. - version = tuple(map(int, __version__.split("."))) - assert version < tuple(map(int, NEXT_REMOVAL.split("."))) - - def test_a_deprecation_warns_until_the_next_removal(): # The helpers a deprecation warns through. match = f"{NEXT_REMOVAL} will refuse it" @@ -284,29 +276,6 @@ def test_a_deprecation_warns_until_the_next_removal(): assert np.array_equal(array(), [1.0, 2.0]) -def test_a_node_state_takes_phase_as_since_scheduled_for_now(): - with pytest.warns( - FutureWarning, match="phase= is renamed since_scheduled" - ): - state = NodeState(age=10.0, phase=4.0) - assert state.since_scheduled == 4.0 - assert state == NodeState(age=10.0, since_scheduled=4.0) - assert ( - repr(state) == "NodeState(age=10.0, alive=True, since_scheduled=4.0)" - ) - with pytest.warns(FutureWarning, match="phase is renamed since_scheduled"): - assert state.phase == 4.0 - with pytest.raises(TypeError, match="since_scheduled alone"): - NodeState(phase=4.0, since_scheduled=4.0) - with warnings.catch_warnings(): - warnings.simplefilter("error") - assert NodeState(since_scheduled=4.0).since_scheduled == 4.0 - assert ( - NodeState(stationary=True, since_scheduled=4.0).since_scheduled - == 4.0 - ) - - def test_an_uncertainty_result_shows_a_summary(): pump = surv.Weibull.fit(np.linspace(200, 1800, 50)) rbd = NonRepairableRBD([("s", "pump"), ("pump", "t")], {"pump": pump}) 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 diff --git a/repyability/tests/test_removed_in_0_15.py b/repyability/tests/test_removed_in_0_15.py new file mode 100644 index 0000000..f8b9aad --- /dev/null +++ b/repyability/tests/test_removed_in_0_15.py @@ -0,0 +1,15 @@ +"""What 0.14 deprecated is gone in 0.15: ``NodeState``'s ``phase``, +renamed ``since_scheduled`` (#350), as a phased mission has phases of its +own.""" + +import pytest + +from repyability import NodeState + + +def test_phase_is_refused_for_since_scheduled(): + with pytest.raises(TypeError, match="phase"): + NodeState(age=10.0, phase=4.0) + state = NodeState(age=10.0, since_scheduled=4.0) + assert state.since_scheduled == 4.0 + assert not hasattr(state, "phase") diff --git a/repyability/utils/deprecation.py b/repyability/utils/deprecation.py index d975974..7df5dc7 100644 --- a/repyability/utils/deprecation.py +++ b/repyability/utils/deprecation.py @@ -1,28 +1,25 @@ """Deprecated names, arguments and calls. -A deprecation gives one minor release's notice: what 0.14 deprecates goes -in 0.15 (``NEXT_REMOVAL``). Each warns with a ``FutureWarning``, which +A deprecation gives one minor release's notice: what 0.15 deprecates goes +in 0.16 (``NEXT_REMOVAL``). Each warns with a ``FutureWarning``, which Python always shows, as the notice is short: an argument renamed (``renamed``; an attribute, ``renamed_attribute``), a result's value called that is now a property (``called``), or arguments a method ignores -(``ignored``). 0.14 deprecates: - -- ``NodeState``'s ``phase``, renamed ``since_scheduled`` (``renamed`` for - the argument, ``renamed_attribute`` for reading it, #350). - -What 0.13 deprecated went in 0.14: ``optimal_inspection_intervals``' -``offsets``, renamed ``offset_shares`` (#222); calling -``CapacityDistribution.mean()``, a property since 0.13 (#235); the -simulation options given to the exact ``mean`` of a ``StandbyModel``, -``LoadSharingModel`` or ``DegradingNode``, which it would ignore (#233); -and ``rbd._sensitivity``'s ``levers`` and ``Lever``, the names callers -used before the levers were public (#244). What 0.12 deprecated went in -0.13: calling ``SparesDemand.mean()`` and ``std()``, properties since 0.12 -(#184), and ``StandbyModel``'s and ``LoadSharingModel``'s ``mc_samples``, -``lower`` and ``seed``, which set a fit to simulated lifetimes that 0.12 -removed (#149). What 0.11 deprecated went in 0.12 (#149): its old names -are refused with the names that took their place -(``refuses_removed_names``, #232). +(``ignored``). 0.15 deprecates nothing. + +What 0.14 deprecated went in 0.15: ``NodeState``'s ``phase``, renamed +``since_scheduled`` (#350). What 0.13 deprecated went in 0.14: +``optimal_inspection_intervals``' ``offsets``, renamed ``offset_shares`` +(#222); calling ``CapacityDistribution.mean()``, a property since 0.13 (#235); +the simulation options given to the exact ``mean`` of a ``StandbyModel``, +``LoadSharingModel`` or ``DegradingNode``, which it would ignore (#233); and +``rbd._sensitivity``'s ``levers`` and ``Lever``, the names callers used before +the levers were public (#244). What 0.12 deprecated went in 0.13: calling +``SparesDemand.mean()`` and ``std()``, properties since 0.12 (#184), and +``StandbyModel``'s and ``LoadSharingModel``'s ``mc_samples``, ``lower`` and +``seed``, which set a fit to simulated lifetimes that 0.12 removed (#149). What +0.11 deprecated went in 0.12 (#149): its old names are refused with the names +that took their place (``refuses_removed_names``, #232). """ import functools @@ -34,8 +31,8 @@ from repyability.utils.wrappers import outside_level -#: The release that removes what 0.14 deprecates. -NEXT_REMOVAL = "0.15" +#: The release that removes what 0.15 deprecates. +NEXT_REMOVAL = "0.16" def _called_warning(name: str, removal: str) -> None: