diff --git a/README.md b/README.md index 0a9065c..552a2c7 100755 --- a/README.md +++ b/README.md @@ -173,6 +173,7 @@ Which analyses each model type has (the method; — where it has none): | **Costs and spares** | | | | | | | Long-run cost rate | — | `expected_cost_rate` | — | — | — | | Expected cost in a window | — | `expected_cost` | — | — | — | +| Cost in each of a run of windows | — | `cost_by_window` | — | — | — | | Simulated cost (distribution) | — | `cost` | — | — | — | | Total cost of ownership | — | `total_cost` | — | — | — | | Spares demand and stock | — | `spares_demand` | — | — | — | diff --git a/docs/api/results.md b/docs/api/results.md index 9b915da..1e4bb3c 100644 --- a/docs/api/results.md +++ b/docs/api/results.md @@ -22,6 +22,10 @@ ::: repyability.ExpectedEvents +::: repyability.EventsByWindow + +::: repyability.CostByWindow + ::: repyability.RateBreakdown ::: repyability.ExpectedCost diff --git a/repyability/__init__.py b/repyability/__init__.py index 44d1b6d..1a61e74 100644 --- a/repyability/__init__.py +++ b/repyability/__init__.py @@ -99,9 +99,11 @@ ConditionalRun, ConfidenceInterval, ControlVariate, + CostByWindow, CostResult, Criticalities, DemonstrationPlan, + EventsByWindow, ExpectedCost, ExpectedEvents, FailureCriticalityIndex, @@ -186,6 +188,8 @@ "CostResult", "ExpectedCost", "ExpectedEvents", + "EventsByWindow", + "CostByWindow", "StructureSummary", "Lever", "RedundancyAllocation", diff --git a/repyability/rbd/_repairable_routes.py b/repyability/rbd/_repairable_routes.py index 61ccf02..1a36dc2 100644 --- a/repyability/rbd/_repairable_routes.py +++ b/repyability/rbd/_repairable_routes.py @@ -691,7 +691,23 @@ def from_long_run(route, reason): "shards", ): out[name] = r.refused(meshed) - out = rbd._meshed_routes(out, free) + # A run of windows: the exact values over time (the window's own + # route), and for the cost its spread simulated from one seeded run. + out["expected_events_by_window"] = out["expected_events"] + if out["expected_cost"].route not in (r.EXACT, r.NUMERICAL): + out["cost_by_window"] = out["expected_cost"] + elif out["cost"].route == r.REFUSED: + out["cost_by_window"] = out["cost"] + else: + out["cost_by_window"] = dataclasses.replace( + out["cost"], + reason="Each window's mean is the difference of expected_cost " + "(exact); its spread comes from the simulations cost runs, " + "one seeded set of histories cut at every window's end.", + twin="", + ) + drawn = {"cost_by_window"} if rbd.has_costs else set() + out = rbd._meshed_routes(out, free, drawn) # Parameter uncertainty (#200): each draw's diagram rebuilt and its # value worked out as the diagram's own, so refused where that is # (a structure too meshed included). diff --git a/repyability/rbd/_window_report.py b/repyability/rbd/_window_report.py new file mode 100644 index 0000000..566c93d --- /dev/null +++ b/repyability/rbd/_window_report.py @@ -0,0 +1,157 @@ +"""A repairable diagram's work and cost in each of a run of consecutive +windows from its state now. + +The exact counts and means over ``[0, t)`` (``expected_events``, +``expected_cost``) add over windows, so a window's own value is the +difference of the values at its two ends. The spread of a cost cannot be +differenced: it comes from simulated histories, each cut at every +window's end. The simulations of a seeded run draw the same lives +whatever the window, so the runs to each end are the same histories cut +off there, and a simulation's cost in a window is its cost to the +window's end less its cost to the window's start. +""" + +import math +from typing import Sequence, Tuple + +import numpy as np + +from repyability.rbd.results import CostByWindow, EventsByWindow +from repyability.utils.checks import real_array + +#: Numerical noise below this (relative to the largest value) is shown as 0. +_NOISE = 1e-9 + + +def window_ends(ends, width, count) -> np.ndarray: + """The windows' end times from ``ends`` (increasing, more than 0), or + ``count`` windows of length ``width``.""" + if ends is None: + if width is None or count is None: + raise ValueError( + "Give the windows' ends, or width and count (windows of " + "that length)." + ) + if width is True or count is True: + raise TypeError("width and count must be numbers.") + width = float(real_array(width, "width")) + if not (math.isfinite(width) and width > 0.0): + raise ValueError(f"width must be more than 0; got {width!r}.") + if int(count) != count or count < 1: + raise ValueError("count must be a whole number, 1 or more.") + return width * np.arange(1, int(count) + 1) + if width is not None or count is not None: + raise ValueError( + "Give the windows' ends, or width and count, not both." + ) + out = np.atleast_1d(real_array(ends, "ends")) + if out.ndim != 1 or out.size == 0 or not np.all(np.isfinite(out)): + raise ValueError(f"ends must be a list of times; got {ends!r}.") + if out[0] <= 0.0 or np.any(np.diff(out) <= 0.0): + raise ValueError( + "ends must increase, the first more than 0 (the windows start " + f"from now); got {ends!r}." + ) + return out.astype(float) + + +def _between(cumulative) -> np.ndarray: + """Each window's own value from the values at its ends, the tiny + negatives numerics leave in a window with nothing in it set to 0.""" + values = np.diff(np.asarray(cumulative, dtype=float), prepend=0.0) + scale = max(1.0, float(np.max(np.abs(cumulative)))) + values[(values < 0.0) & (values > -_NOISE * scale)] = 0.0 + return values + + +def expected_events_by_window( + rbd, ends, width, count, state, working_nodes, broken_nodes, method +) -> EventsByWindow: + """See ``RepairableRBD.expected_events_by_window``.""" + stops = window_ends(ends, width, count) + total = rbd.expected_events( + stops, + working_nodes=working_nodes, + broken_nodes=broken_nodes, + method=method, + state=state, + ) + + def per_node(table: dict) -> dict: + return {node: _between(values) for node, values in table.items()} + + return EventsByWindow( + starts=np.concatenate(([0.0], stops[:-1])), + ends=stops, + system_failures=_between(total.system_failures), + system_planned_outages=_between(total.system_planned_outages), + system_downtime=_between(total.system_downtime), + node_failures=per_node(total.node_failures), + node_corrective=per_node(total.node_corrective), + node_preventive=per_node(total.node_preventive), + node_inspections=per_node(total.node_inspections), + node_downtime=per_node(total.node_downtime), + ) + + +def cost_by_window( + rbd, + ends, + width, + count, + state, + mc_samples, + seed, + engine, + levels: Sequence[float], + working_nodes, + broken_nodes, + method, +) -> CostByWindow: + """See ``RepairableRBD.cost_by_window``.""" + stops = window_ends(ends, width, count) + levels_t: Tuple[float, ...] = tuple( + float(q) for q in np.atleast_1d(real_array(levels, "levels")) + ) + if not all(0.0 <= q <= 100.0 for q in levels_t): + raise ValueError( + f"levels must be percentiles, 0 to 100; got {levels!r}." + ) + exact = rbd.expected_cost( + stops, + working_nodes=working_nodes, + broken_nodes=broken_nodes, + method=method, + state=state, + ).mean + cumulative = [] + for stop in stops: + run = rbd.cost( + float(stop), + working_nodes=working_nodes, + broken_nodes=broken_nodes, + method=method, + mc_samples=mc_samples, + seed=seed, + engine=engine, + state=state, + control_variate=False, + conditional=False, + ) + cumulative.append(np.asarray(run.samples, dtype=float)) + to_end = np.vstack(cumulative) + samples = np.diff(to_end, axis=0, prepend=0.0) + scale = max(1.0, float(np.max(np.abs(to_end)))) + samples[(samples < 0.0) & (samples > -_NOISE * scale)] = 0.0 + n = samples.shape[1] + return CostByWindow( + starts=np.concatenate(([0.0], stops[:-1])), + ends=stops, + mean=_between(exact), + simulated_mean=samples.mean(axis=1), + standard_error=samples.std(axis=1, ddof=1) / math.sqrt(n), + samples=samples, + n_simulations=n, + levels=levels_t, + percentiles=np.percentile(samples, levels_t, axis=1).T, + ) diff --git a/repyability/rbd/repairable_rbd.py b/repyability/rbd/repairable_rbd.py index 1a1c042..71731a1 100644 --- a/repyability/rbd/repairable_rbd.py +++ b/repyability/rbd/repairable_rbd.py @@ -46,7 +46,9 @@ AvailabilityResult, CapacityDistribution, ConfidenceInterval, + CostByWindow, CostResult, + EventsByWindow, ExpectedCost, ExpectedEvents, Lever, @@ -79,6 +81,7 @@ _repairable_uncertainty, _runs, _spec, + _window_report, _windows, ) from repyability.rbd._common import ( @@ -3523,6 +3526,165 @@ def expected_events( state=state, ) + def expected_events_by_window( + self, + ends=None, + *, + width: Optional[float] = None, + count: Optional[int] = None, + state=None, + working_nodes: Optional[Collection[Hashable]] = None, + broken_nodes: Optional[Collection[Hashable]] = None, + method: str = "p", + ) -> EventsByWindow: + """The work expected in each of a run of consecutive windows from + the system's state now: how many replacements, repairs and tests + each component has in each window, and how many failures and + planned outages the system has. + + A planner building next year's work and budget asks exactly this: + "in each of the next twelve windows of 730 hours, what is due?" + Each window's number is its own (not a running total), exact, with + no simulation: the difference of + [`expected_events`][repyability.RepairableRBD.expected_events] at + the window's two ends, which add up over windows. Give the + windows' ``ends`` (uneven windows are fine), or a ``width`` and a + ``count``. The first window starts now. + + A component's task and the system's outage are different counts: + a task done while the system is already down (a hidden failure + found earlier, say) causes no new outage. + + Parameters + ---------- + ends : array-like, optional + The windows' end times from now, increasing. + width : float, optional + Instead of ``ends``: the length of each window. + count : int, optional + Instead of ``ends``: how many windows of ``width``. + state : dict or str, optional + Where the components are now, as for + [`expected_events`][repyability.RepairableRBD.expected_events]. + working_nodes : list, optional + As for `expected_events`. + broken_nodes : list, optional + As for `expected_events`. + method : str, optional + As for `expected_events`. + + Returns + ------- + EventsByWindow + The windows and each count in each (arrays); a hair below 0 + from numerics, in a window with nothing in it, is shown as 0. + + Raises + ------ + ValueError, NotImplementedError + As for `expected_events`, or if the windows are not given + properly. + + Examples + -------- + See [`EventsByWindow`][repyability.EventsByWindow]. + """ + return _window_report.expected_events_by_window( + self, + ends, + width, + count, + state, + working_nodes, + broken_nodes, + method, + ) + + def cost_by_window( + self, + ends=None, + *, + width: Optional[float] = None, + count: Optional[int] = None, + state=None, + mc_samples: Optional[int] = None, + seed: Optional[int] = None, + engine: str = "auto", + levels=(10.0, 50.0, 90.0), + working_nodes: Optional[Collection[Hashable]] = None, + broken_nodes: Optional[Collection[Hashable]] = None, + method: str = "p", + ) -> CostByWindow: + """What the system costs in each of a run of consecutive windows + from its state now, and how much that can vary. + + A budget starts from the plant as it is, and needs each window's + own spend, not the running total, and what could happen in it. The + mean of each window is exact (the difference of + [`expected_cost`][repyability.RepairableRBD.expected_cost] at its + ends). The spread is simulated, as a spread must be, from one + seeded set of histories cut at every window's end, so the windows + are of the same histories: a simulation's cost in a window is its + cost to the window's end less its cost to the window's start. Give + a ``seed`` for a repeatable answer. + + Costs usually come in steps (a replacement at a fixed price), so + percentiles can hide that most windows cost exactly one amount: + read the outcomes (``outcomes(window)``) and the chance of + exceeding a budget (``probability_above(budget)``) too. + + Parameters + ---------- + ends : array-like, optional + The windows' end times from now, increasing. + width : float, optional + Instead of ``ends``: the length of each window. + count : int, optional + Instead of ``ends``: how many windows of ``width``. + state : dict or str, optional + Where the components are now, as for + [`expected_cost`][repyability.RepairableRBD.expected_cost]. + mc_samples : int, optional + The number of simulations, as for + [`cost`][repyability.RepairableRBD.cost]. + seed : int, optional + The random seed, as for `cost`. + engine : str, optional + The simulation engine, as for `cost`. + levels : array-like, optional + The percentiles to report, by default 10, 50 and 90. + working_nodes : list, optional + As for `expected_cost`. + broken_nodes : list, optional + As for `expected_cost`. + method : str, optional + As for `expected_cost`. + + Returns + ------- + CostByWindow + The exact mean, the simulations' own mean and standard error, + every simulation's cost in each window, and the percentiles. + + Examples + -------- + See [`CostByWindow`][repyability.CostByWindow]. + """ + return _window_report.cost_by_window( + self, + ends, + width, + count, + state, + mc_samples, + seed, + engine, + levels, + working_nodes, + broken_nodes, + method, + ) + def expected_cost( self, t, diff --git a/repyability/rbd/results.py b/repyability/rbd/results.py index 229af2b..5f42b0e 100644 --- a/repyability/rbd/results.py +++ b/repyability/rbd/results.py @@ -217,6 +217,182 @@ class StructureSummary(_ResultMapping): too_meshed: Optional[str] +@dataclass +class EventsByWindow(_ResultMapping): + """What a system is expected to do in each of a run of consecutive + windows from now: returned by + [`RepairableRBD.expected_events_by_window`][repyability.RepairableRBD.expected_events_by_window]. + + The windows are ``[starts[i], ends[i])``, the first starting at 0 + (now, or the state given). Every count is exact (see + [`ExpectedEvents`][repyability.ExpectedEvents]) and is the window's + own, not a running total: the difference of the expected counts at its + two ends, which add up over windows. A value that numerics leave a + hair below 0 for a window with no events is shown as 0. + + Attributes + ---------- + starts : numpy.ndarray + Each window's start, in the models' time unit. + ends : numpy.ndarray + Each window's end. + system_failures : numpy.ndarray + The expected system failures in each window. + system_planned_outages : numpy.ndarray + The expected planned outages of the system in each window. + system_downtime : numpy.ndarray + The expected time the system is down in each window. + node_failures : dict + Each component's expected failures in each window (an array per + component). + node_corrective : dict + Each component's expected corrective actions in each window. + node_preventive : dict + Each component's expected preventive replacements in each window. + node_inspections : dict + Each component's expected tests in each window. + node_downtime : dict + Each component's expected down time in each window. + + Examples + -------- + One component with a block replacement every 100 hours, over three + windows of 100 hours from 40 hours into its interval: its replacement + is due in 60 hours, and then every 100: + + >>> import surpyval as surv + >>> from repyability import NodeState, RepairableRBD + >>> rbd = RepairableRBD( + ... [("s", "c"), ("c", "t")], + ... {"c": { + ... "reliability": surv.Exponential.from_params([1e-9]), + ... "repairability": "instant", + ... "preventive": {"interval": 100.0, "policy": "block"}, + ... }}, + ... ) + >>> state = {"c": NodeState(since_scheduled=40.0)} + >>> windows = rbd.expected_events_by_window( + ... width=100.0, count=3, state=state + ... ) + >>> windows.node_preventive["c"].round(3).tolist() + [1.0, 1.0, 1.0] + """ + + starts: np.ndarray + ends: np.ndarray + system_failures: np.ndarray + system_planned_outages: np.ndarray + system_downtime: np.ndarray + node_failures: Dict[Hashable, np.ndarray] + node_corrective: Dict[Hashable, np.ndarray] + node_preventive: Dict[Hashable, np.ndarray] + node_inspections: Dict[Hashable, np.ndarray] + node_downtime: Dict[Hashable, np.ndarray] + + +@dataclass +class CostByWindow(_ResultMapping): + """What a system costs in each of a run of consecutive windows from + now, and how much that can vary: returned by + [`RepairableRBD.cost_by_window`][repyability.RepairableRBD.cost_by_window]. + + The windows are ``[starts[i], ends[i])``. ``mean`` is exact. The spread + comes from one set of simulated histories cut at each window's end, so + every window's costs are of the same histories (a window's cost is + the cost to its end less the cost to its start). A cost that comes + in steps (a replacement at a fixed price, say) has few distinct + values, and percentiles hide that: ``outcomes`` gives each cost and + its share, and ``probability_above`` the chance of overspending a + budget. Ranges do not add over windows: use the ``total`` row for a + run's total. + + Attributes + ---------- + starts : numpy.ndarray + Each window's start. + ends : numpy.ndarray + Each window's end. + mean : numpy.ndarray + The exact expected cost of each window. + simulated_mean : numpy.ndarray + The simulations' own mean of each window's cost (it agrees with + ``mean`` to within ``standard_error``). + standard_error : numpy.ndarray + The standard error of ``simulated_mean``. + samples : numpy.ndarray + Each simulation's cost in each window, shape ``(windows, + simulations)``; none is below 0. + n_simulations : int + The number of simulations. + levels : tuple of float + The percentiles of ``percentiles``. + percentiles : numpy.ndarray + Each window's cost percentiles at ``levels``, shape ``(windows, + len(levels))``. + """ + + starts: np.ndarray + ends: np.ndarray + mean: np.ndarray + simulated_mean: np.ndarray + standard_error: np.ndarray + samples: np.ndarray + n_simulations: int + levels: Tuple[float, ...] + percentiles: np.ndarray + + def percentile(self, q: float) -> np.ndarray: + """Each window's ``q``-th percentile of cost. + + Parameters + ---------- + q : float + The percentile, from 0 to 100. + """ + return np.percentile(self.samples, q, axis=1) + + def total_percentile(self, q: float) -> float: + """The ``q``-th percentile of the cost of all the windows + together: not the sum of the windows' (ranges do not add). + + Parameters + ---------- + q : float + The percentile, from 0 to 100. + """ + return float(np.percentile(self.samples.sum(axis=0), q)) + + def probability_above(self, budget: float) -> np.ndarray: + """Each window's chance of costing more than ``budget``. + + Parameters + ---------- + budget : float + The cost to exceed. + """ + return np.mean(self.samples > budget, axis=1) + + def outcomes(self, window: int, decimals: int = 6) -> list: + """The distinct costs of a window (counted from 0) and how often + each happened, cheapest first: ``[{"cost": ..., "share": ...}]`` + (costs rounded to ``decimals`` places). + + Parameters + ---------- + window : int + The window, counted from 0. + decimals : int, optional + The places costs are rounded to, by default 6. + """ + costs, counts = np.unique( + np.round(self.samples[window], decimals), return_counts=True + ) + return [ + {"cost": float(c), "share": float(n) / self.n_simulations} + for c, n in zip(costs, counts) + ] + + @dataclass class ConfidenceInterval(_ResultMapping): """A Monte-Carlo estimate with its sampling uncertainty. diff --git a/repyability/tests/capabilities.py b/repyability/tests/capabilities.py index d968211..6f13ead 100644 --- a/repyability/tests/capabilities.py +++ b/repyability/tests/capabilities.py @@ -96,7 +96,11 @@ ("Mean up and down times", ("mean_up_time", "mean_down_time")), ( "Expected failures and events in a window", - ("expected_events", "expected_failures"), + ( + "expected_events", + "expected_failures", + "expected_events_by_window", + ), ), ( "Simulated availability (distribution)", @@ -116,6 +120,7 @@ [ ("Long-run cost rate", ("expected_cost_rate",)), ("Expected cost in a window", ("expected_cost",)), + ("Cost in each of a run of windows", ("cost_by_window",)), ("Simulated cost (distribution)", ("cost",)), ("Total cost of ownership", ("total_cost",)), ("Spares demand and stock", ("spares_demand", "spares_stock")), diff --git a/repyability/tests/test_analysis_routes.py b/repyability/tests/test_analysis_routes.py index cd45e42..87d0887 100644 --- a/repyability/tests/test_analysis_routes.py +++ b/repyability/tests/test_analysis_routes.py @@ -255,6 +255,12 @@ def stepped(rbd): "expected_failures": lambda rbd: rbd.expected_failures([10.0, 200.0]), "expected_events": lambda rbd: rbd.expected_events(200.0), "expected_cost": lambda rbd: rbd.expected_cost(200.0), + "expected_events_by_window": lambda rbd: rbd.expected_events_by_window( + width=100.0, count=2 + ), + "cost_by_window": lambda rbd: rbd.cost_by_window( + width=100.0, count=2, mc_samples=20, seed=1 + ), "point_capacity": lambda rbd: rbd.point_capacity([10.0, 200.0]), "mission_capacity": lambda rbd: rbd.mission_capacity(200.0), "availability": lambda rbd: rbd.availability(200.0, mc_samples=20, seed=1), diff --git a/repyability/tests/test_window_report.py b/repyability/tests/test_window_report.py new file mode 100644 index 0000000..5b1759e --- /dev/null +++ b/repyability/tests/test_window_report.py @@ -0,0 +1,246 @@ +"""The work and cost in each of a run of windows from a state +(``expected_events_by_window``, ``cost_by_window``): checked against the +exact counts over the whole span (windows add), hand-worked schedules, +a closed form, and the simulation.""" + +import numpy as np +import pytest +import surpyval as surv + +from repyability import CostByWindow, EventsByWindow, NodeState, RepairableRBD + +W = surv.Weibull.from_params +E = surv.Exponential.from_params +FIXED = surv.ExactEventTime.from_params + + +def plant(): + return RepairableRBD( + [ + ("s", "pump"), + ("pump", "valve"), + ("valve", "sensor"), + ("sensor", "t"), + ], + { + "pump": { + "reliability": W([9000.0, 2.5]), + "repairability": surv.LogNormal.from_params( + [np.log(8.0), 0.5] + ), + "preventive": {"interval": 6000.0, "duration": FIXED([4.0])}, + }, + "valve": { + "reliability": W([20000.0, 1.5]), + "repairability": E([1 / 12.0]), + "preventive": { + "interval": 8760.0, + "policy": "block", + "duration": FIXED([6.0]), + }, + }, + "sensor": { + "reliability": E([1 / 50000.0]), + "repairability": "instant", + "inspection": {"interval": 4380.0}, + }, + }, + ) + + +STATE = { + "pump": NodeState(age=5200.0), + "valve": NodeState(age=7000.0, since_scheduled=7000.0), + "sensor": NodeState(age=3000.0, since_scheduled=3000.0), +} + + +def bearing(): + return RepairableRBD( + [("s", "b"), ("b", "t")], + { + "b": { + "reliability": W([365.0, 2.0]), + "repairability": "instant", + "replace_cost": 400.0, + "preventive": { + "interval": 180.0, + "policy": "block", + "cost": 100.0, + }, + } + }, + ) + + +BEARING_STATE = {"b": NodeState(age=224.0, since_scheduled=159.0)} + + +def test_the_windows_add_up_to_the_whole_span(): + rbd = plant() + windows = rbd.expected_events_by_window(width=730.0, count=12, state=STATE) + whole = rbd.expected_events(730.0 * 12, state=STATE) + assert isinstance(windows, EventsByWindow) + for table, total in ( + (windows.node_preventive, whole.node_preventive), + (windows.node_failures, whole.node_failures), + (windows.node_inspections, whole.node_inspections), + ): + for node in total: + assert table[node].sum() == pytest.approx(total[node], abs=1e-9) + assert windows.system_planned_outages.sum() == pytest.approx( + whole.system_planned_outages, abs=1e-9 + ) + assert min(v.min() for v in windows.node_failures.values()) >= 0.0 + + +def test_the_planners_year_matches_the_hand_worked_table(): + windows = plant().expected_events_by_window( + width=730.0, count=12, state=STATE + ) + pump = windows.node_preventive["pump"].round(2) + assert pump.tolist() == [0, 0.9, 0, 0, 0, 0, 0, 0, 0.05, 0.65, 0, 0.01] + # The pump's first chance is R(6000) / R(5200), by hand. + life = W([9000.0, 2.5]) + assert windows.node_preventive["pump"][1] == pytest.approx( + life.sf(6000.0) / life.sf(5200.0), abs=1e-6 + ) + # Block replacement of the valve at 8760: certain, 1760 h away, in + # window 3 (1460 to 2190); the sensor's tests at 4380: 1380 h away, + # in window 2, then 4380 h later, in window 8. + # (skipped while the valve is down in a repair, so a hair below 1) + assert windows.node_preventive["valve"].round(2).tolist() == ( + [0, 0, 1] + [0] * 9 + ) + assert windows.node_inspections["sensor"].round(6).tolist() == ( + [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0] + ) + # A task and a system outage are different counts. + assert ( + windows.system_planned_outages[1] < windows.node_preventive["pump"][1] + ) + + +def test_uneven_windows_and_width_and_count_agree(): + rbd = plant() + even = rbd.expected_events_by_window(width=500.0, count=4, state=STATE) + listed = rbd.expected_events_by_window( + [500.0, 1000.0, 1500.0, 2000.0], state=STATE + ) + np.testing.assert_allclose( + even.node_failures["pump"], listed.node_failures["pump"] + ) + uneven = rbd.expected_events_by_window([300.0, 2000.0], state=STATE) + assert uneven.starts.tolist() == [0.0, 300.0] + assert uneven.node_failures["pump"].sum() == pytest.approx( + even.node_failures["pump"].sum(), abs=1e-9 + ) + + +def test_a_unit_in_its_long_run_state_fails_at_a_steady_rate(): + lam, mu = 0.01, 0.2 + rbd = RepairableRBD( + [("s", "c"), ("c", "t")], + {"c": {"reliability": E([lam]), "repairability": E([mu])}}, + ) + windows = rbd.expected_events_by_window( + width=50.0, count=4, state="stationary" + ) + # lam * availability * time, in every window. + expected = lam * mu / (lam + mu) * 50.0 + np.testing.assert_allclose(windows.system_failures, expected, rtol=1e-5) + + +def test_the_windows_must_be_given_properly(): + rbd = plant() + for bad in ( + dict(), + dict(width=10.0), + dict(ends=[10.0], width=10.0, count=1), + dict(ends=[0.0, 5.0]), + dict(ends=[5.0, 5.0]), + dict(width=-1.0, count=2), + dict(width=10.0, count=0), + ): + with pytest.raises(ValueError): + rbd.expected_events_by_window(**bad) + + +def test_costs_by_hand_on_a_nearly_unfailing_bearing(): + # Shuts every 180 days, the next in 21: at 21, 201, 381, ... A bearing + # that never fails costs exactly $100 a shut. + rbd = RepairableRBD( + [("s", "b"), ("b", "t")], + { + "b": { + "reliability": E([1e-12]), + "repairability": "instant", + "replace_cost": 400.0, + "preventive": { + "interval": 180.0, + "policy": "block", + "cost": 100.0, + }, + } + }, + ) + out = rbd.cost_by_window( + [82.0, 447.0], state=BEARING_STATE, mc_samples=50, seed=1 + ) + assert out.mean.round(6).tolist() == [100.0, 200.0] + assert out.samples.tolist() == [[100.0] * 50, [200.0] * 50] + assert out.outcomes(1) == [{"cost": 200.0, "share": 1.0}] + assert out.probability_above(150.0).tolist() == [0.0, 1.0] + + +def test_cost_means_agree_with_the_simulation_and_add_up(): + rbd = bearing() + ends = [82.0, 447.0, 813.0, 1086.0] + out = rbd.cost_by_window( + ends, state=BEARING_STATE, mc_samples=20_000, seed=1 + ) + assert isinstance(out, CostByWindow) + whole = rbd.expected_cost(1086.0, state=BEARING_STATE).mean + assert out.mean.sum() == pytest.approx(whole, rel=1e-12) + # README's every-shut policy: 1,155 over 1,086 days. + assert whole == pytest.approx(1155.0, abs=0.5) + assert np.all( + np.abs(out.simulated_mean - out.mean) < 4 * out.standard_error + ) + assert out.samples.min() >= 0.0 + # The same histories cut at each end: the windows add to one run. + run = rbd.cost( + 1086.0, + state=BEARING_STATE, + mc_samples=20_000, + seed=1, + control_variate=False, + conditional=False, + ) + np.testing.assert_allclose(out.samples.sum(axis=0), run.samples) + + +def test_outcomes_percentiles_and_budgets(): + rbd = bearing() + out = rbd.cost_by_window( + [263.0, 629.0], state=BEARING_STATE, mc_samples=20_000, seed=1 + ) + first = out.outcomes(0) + assert sum(o["share"] for o in first) == pytest.approx(1.0) + assert [o["cost"] for o in first] == sorted(o["cost"] for o in first) + # Costs come in steps of $100 shut and $400 failure. + assert all(o["cost"] % 100.0 == 0.0 for o in first) + # README: more than $200 (any failure) in about 29% of the first window. + assert out.probability_above(200.0)[0] == pytest.approx(0.29, abs=0.02) + assert out.percentiles.shape == (2, 3) + assert out.levels == (10.0, 50.0, 90.0) + np.testing.assert_allclose(out.percentiles[:, 1], out.percentile(50.0)) + # Ranges do not add. + assert out.total_percentile(90.0) <= out.percentile(90.0).sum() + assert out.to_dict()["n_simulations"] == 20_000 + + +def test_the_same_seed_gives_the_same_answer(): + a = bearing().cost_by_window(width=100.0, count=3, mc_samples=200, seed=5) + b = bearing().cost_by_window(width=100.0, count=3, mc_samples=200, seed=5) + np.testing.assert_array_equal(a.samples, b.samples)