diff --git a/CHANGELOG.md b/CHANGELOG.md index 58ce42f..71cf7ca 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -101,6 +101,34 @@ Equation numbers below are that manual's. PPU 5.74; ȳ alone gives 5.73). It equals ȳ in ADS 1. ### Fixed +- **Stratified charts take time in proportion to the data, however many strata there are.** + The stratified X/mR path looked each stratum up with a scan of the whole table, so its cost + grew with strata × rows. With every row its own condition, an X/mR pair took 48 s at 30,165 + rows, and an estimated 80 minutes at 307,683 (a gene-expression study, 6,033 genes × 51 + subjects). It now takes 0.2 s and 1.9 s. The same scan sat in the stratified Xbar/S loops and + the lane-boundary loop, the per-stratum limits were computed one `pd.Series` at a time, and + two lists of strata were searched as lists. Results are unchanged: every chart, statistic, + lane boundary and stratum list was compared before and after on the reference data and a + 30,165-row gene-effects file. `scripts/benchmark.py` adds a one-row-strata scenario, and + `tests/test_performance.py` checks that it scales linearly. Stratified Xbar/S still costs + about 2.5 ms per stratum (15 s at 6,033): linear, but not yet vectorised. +- **Charts with many lanes draw in a fraction of a second.** Each lane divider was added to the + figure with its own `add_shape` call, and plotly re-validates every shape already there, so the + cost grew with the square of the number of lanes: 800 dividers took 27 s, and an X chart of a + 6,033-gene study (one lane per gene) would have taken about 25 minutes. They are now added in + one assignment, which builds the same figure (0.03 s for 800). Above 1,000 + (`plotting.lane_boundaries.MAX_LANE_DIVIDERS`) none are drawn: the dividers would merge into a + solid band at any chart width. +- **`detect_signals()` works on phased charts.** It raised `TypeError` on every `phased=True` + chart: each phase has its own centre, so the statistics carry `center=None`, and the zone + calculation read that scalar. The chart already marks run rules not applicable, but only the + plot read the flag. Phased charts now get Rule 1 by default, judged against each phase's + limits, with the same points the chart flags; an explicit rule list runs as given, against each + row's own centre. The limits and the flagged points on the chart were always correct. +- **`value='response'` charts the response**, as the `execute()` signature and docstring say. + Every chart refused it ("Valid charts for RESPONSE: None"): residuals are told apart by a leading + `R`, which `RESPONSE` also has. `execute()`, `supports_calibration()` and `why_not()` now treat + it, in any case, exactly as `value=None`. - **Xbar and S charts with subgroups of unequal size follow the manual (Eqs 11-5, 11-16, 11-17).** sigma-hat is the average of each subgroup's S_r / c4(N_r) (was S-bar / c4(N_r)), and each subgroup's limits use its own N_r, so they step with subgroup size, as VAS draws them. N_r now diff --git a/processbehavior/analysis.py b/processbehavior/analysis.py index 050d687..8db5514 100644 --- a/processbehavior/analysis.py +++ b/processbehavior/analysis.py @@ -197,6 +197,65 @@ def _split_strata_by_sufficiency( return insufficient, published +def _stratum_row_positions(df: pd.DataFrame, stratify_col: str) -> dict: + """Row positions of every stratum in `df`, from one grouping pass. + + The per-stratum loops used to select with ``df[df[stratify_col] == stratum]``, which + compares the whole column for every stratum: O(strata x rows). With one row per + stratum that is quadratic in the data: 48s for an X/mR pair at 30,165 rows, an + estimated 80 minutes at 307,683. Look rows up here instead with :func:`_stratum_rows`. + + NaN strata are absent (groupby drops them), so they select no rows, exactly as + ``col == NaN`` did. + """ + return df.groupby(stratify_col, sort=False, observed=True).indices + + +def _stratum_rows(df: pd.DataFrame, positions: dict, stratum) -> pd.DataFrame: + """The rows of `df` in `stratum`, in `df`'s order (a view; copy before writing).""" + pos = positions.get(stratum) + return df.iloc[pos] if pos is not None else df.iloc[0:0] + + +def _mr_family_limits(mr_spec: _MRChartSpec, *, mean, mR) -> pd.DataFrame: + """Limits for one MR-family chart per row of a per-stratum or per-phase table. + + The array form of the scalar ``calculate_limits`` these sites used to call through + ``DataFrame.apply(axis=1)``: the same formulae, 10-12 manual Eqs 12-10/12-11 (X) and + 12-5 (mR). The scalar call built a ``pd.Series`` per stratum, about 0.1 ms each, so + it cost 30s per chart at 300,000 strata. ``mean`` is ignored for the mR chart. + """ + if mr_spec.limits_type == 'XmR': + return calculate_limits_vectorized('XmR', mean=mean, mR=mR) + return calculate_limits_vectorized(mr_spec.limits_type, mR=mR) + + +def _mr_stratum_statistics(mr_spec: _MRChartSpec, grouped: pd.DataFrame, stratify_col: str, + lpl_col: str, upl_col: str, round_to: int) -> dict: + """The per-stratum statistics dict of a stratified MR-family chart, in one pass. + + `grouped` has one row per stratum. Replaces a loop that found each stratum's row + with ``grouped[grouped[stratify_col] == stratum]``: O(strata^2). The values stay + numpy floats (``to_numpy``, not ``tolist``), as that loop's row lookup returned them. + """ + n = _mr_family_subgroup_n(mr_spec.chart_type) + return { + stratum: { + 'N': n, + 'center': round(center, round_to), + 'lpl': round(lpl, round_to), + 'upl': round(upl, round_to), + } + for stratum, center, lpl, upl in zip( + grouped[stratify_col].tolist(), + grouped['center'].to_numpy(), + grouped[lpl_col].to_numpy(), + grouped[upl_col].to_numpy(), + strict=True, + ) + } + + def _no_replication_alternative(chart: str) -> str: """What to reach for instead, when a subgroup chart has nothing to subgroup.""" return ( @@ -278,7 +337,9 @@ def _build_stratified_mr_payload( published_strata = list(strata) else: metadata['insufficient_strata'] = insufficient_strata or None - published_strata = [s for s in strata if s not in (insufficient_strata or [])] + # A set: against the list this was O(strata^2) when most strata are insufficient. + insufficient_set = set(insufficient_strata or ()) + published_strata = [s for s in strata if s not in insufficient_set] return { mr_spec.chart_type: { 'data': chart_out, @@ -951,6 +1012,28 @@ def _resolve_by_grouping(self, value_col: str) -> tuple[list[str], str | None, l groupby_cols = list(by) return groupby_cols, None, [] + def _stratum_lane_boundaries(self, out: pd.DataFrame, stratify_col: str, strata: list, + collapsed_factors: list[str]) -> dict: + """Lane boundaries for each stratum of a stratified chart, keyed by stratum. + + Strata with no boundary are left out. A stratum of fewer than two rows can't have + one, so it is skipped without the per-call cost, which matters when nearly every + row is its own stratum. + """ + lane_boundaries: dict = {} + if not collapsed_factors: + return lane_boundaries + positions = _stratum_row_positions(out, stratify_col) + for stratum in strata: + pos = positions.get(stratum) + if pos is None or len(pos) < 2: + continue + stratum_data = out.iloc[pos].reset_index(drop=True) + boundaries = self._calculate_lane_boundaries(stratum_data, collapsed_factors) + if boundaries: + lane_boundaries[stratum] = boundaries + return lane_boundaries + def _calculate_lane_boundaries(self, df: pd.DataFrame, collapsed_vars: list[str]) -> list[dict]: """ Calculate lane boundaries where collapsed factors change. @@ -1368,6 +1451,7 @@ def _calculate_xbar_stratified( _limits_col = self._resolve_limits_column(value_col, df) strata = df[stratify_col].unique().tolist() + positions = _stratum_row_positions(df, stratify_col) all_xbar_frames = [] chart_statistics = {} intermediates_per_stratum = {} @@ -1378,7 +1462,7 @@ def _calculate_xbar_stratified( insufficient_strata: list = [] for stratum in strata: - sdf = df[df[stratify_col] == stratum].copy() + sdf = _stratum_rows(df, positions, stratum).copy() # Group by time within this stratum # Note: we always compute mean directly from filtered data, NOT using @@ -1457,7 +1541,8 @@ def _calculate_xbar_stratified( chart_out = pd.concat(all_xbar_frames, ignore_index=True) # Publish only strata a caller can actually focus() into. - published_strata = [s for s in strata if s not in set(insufficient_strata)] + insufficient_set = set(insufficient_strata) + published_strata = [s for s in strata if s not in insufficient_set] # Select output columns group_col = groupby_cols[0] if groupby_cols else None @@ -1541,10 +1626,12 @@ def _calculate_s_stratified( df['_stratify_key'] = encode_rsg_series(df, stratify_by) stratify_col = '_stratify_key' strata = df[stratify_col].unique().tolist() + positions = _stratum_row_positions(df, stratify_col) per_stratum = None insufficient_strata = [] - all_strata = list(strata) + [s for s in insufficient_strata if s not in strata] + strata_set = set(strata) + all_strata = list(strata) + [s for s in insufficient_strata if s not in strata_set] all_s_frames = [] chart_statistics = {} @@ -1562,7 +1649,7 @@ def _calculate_s_stratified( n_max = inter['n_max'] n_avg = inter['n_avg'] else: - sdf = df[df[stratify_col] == stratum].copy() + sdf = _stratum_rows(df, positions, stratum).copy() _limits_col = self._resolve_limits_column(value_col, sdf) agg_dict = { 's': pd.NamedAgg(column=_limits_col, aggfunc='std'), @@ -1619,7 +1706,8 @@ def _calculate_s_stratified( _raise_no_replication('S', sds, insufficient_strata, len(all_strata)) chart_out = pd.concat(all_s_frames, ignore_index=True) - published_strata = [s for s in all_strata if s not in set(insufficient_strata)] + insufficient_set = set(insufficient_strata) + published_strata = [s for s in all_strata if s not in insufficient_set] group_col = groupby_cols[0] if groupby_cols else None cols_to_keep = [stratify_col, group_col, 's', 'center', 'lpl', 'upl', 'beyond_limits'] @@ -1983,17 +2071,7 @@ def _calculate_mr_chart_from_precomputed( # Recompute R-specific limits per phase within each stratum r_phase_stats = phase_stats[[stratify_col, '_phase_id', 'mR']].copy() - r_lims = r_phase_stats.apply( - lambda row: calculate_limits( - mean=0, - sd=0, - N=0, - mR=row['mR'], - limits_type=mr_spec.limits_type, - round_to=spec.round_to, - ), - axis=1, - ) + r_lims = _mr_family_limits(mr_spec, mean=0, mR=r_phase_stats['mR']) r_phase_stats = _join_limits(r_phase_stats, r_lims) r_phase_stats['center'] = r_phase_stats['mR'] @@ -2048,17 +2126,7 @@ def _calculate_mr_chart_from_precomputed( # grouped has columns: stratify_col, center (=mean), mR, lpl, upl (X limits) r_grouped = grouped.rename(columns={'center': 'xmr_center', 'mR': 'center'}) - r_lims = r_grouped.apply( - lambda row: calculate_limits( - mean=0, - sd=0, - N=0, - mR=row['center'], - limits_type=mr_spec.limits_type, - round_to=spec.round_to, - ), - axis=1, - ) + r_lims = _mr_family_limits(mr_spec, mean=0, mR=r_grouped['center']) r_grouped = _join_limits(r_grouped, r_lims, rsuffix='_r') @@ -2084,15 +2152,9 @@ def _calculate_mr_chart_from_precomputed( r_out = self._add_beyond_limits_flag(r_out, value_col=plot_col) # Build statistics - chart_statistics = {} - for stratum in strata: - row = r_grouped[r_grouped[stratify_col] == stratum].iloc[0] - chart_statistics[stratum] = { - 'N': _mr_family_subgroup_n(mr_spec.chart_type), - 'center': round(row['center'], spec.round_to), - 'lpl': round(row[r_lpl_col], spec.round_to), - 'upl': round(row[r_upl_col], spec.round_to), - } + chart_statistics = _mr_stratum_statistics( + mr_spec, r_grouped, stratify_col, r_lpl_col, r_upl_col, spec.round_to, + ) # Format output chart_out = self._build_output_columns( @@ -2186,30 +2248,10 @@ def _calculate_mr_chart_stratified( # noqa: C901 # 5. Compute per-phase limits if mr_spec.center_source == 'mean': phase_stats['center'] = phase_stats['_mean'] - phase_lims = phase_stats.apply( - lambda row: calculate_limits( - mean=row['_mean'], - sd=0, - N=0, - mR=row['mR'], - limits_type=mr_spec.limits_type, - round_to=spec.round_to, - ), - axis=1, - ) + phase_lims = _mr_family_limits(mr_spec, mean=phase_stats['_mean'], mR=phase_stats['mR']) else: # center_source == 'mR' (R chart) phase_stats['center'] = phase_stats['mR'] - phase_lims = phase_stats.apply( - lambda row: calculate_limits( - mean=0, - sd=0, - N=0, - mR=row['mR'], - limits_type=mr_spec.limits_type, - round_to=spec.round_to, - ), - axis=1, - ) + phase_lims = _mr_family_limits(mr_spec, mean=0, mR=phase_stats['mR']) phase_stats = _join_limits(phase_stats, phase_lims) # 6. Merge per-phase limits back to row-level data @@ -2228,13 +2270,7 @@ def _calculate_mr_chart_stratified( # noqa: C901 out = self._add_beyond_limits_flag(out, value_col=plot_col) # 9. Lane boundaries per stratum - lane_boundaries = {} - if collapsed_factors: - for stratum in strata: - stratum_data = out[out[stratify_col] == stratum].reset_index(drop=True) - boundaries = self._calculate_lane_boundaries(stratum_data, collapsed_factors) - if boundaries: - lane_boundaries[stratum] = boundaries + lane_boundaries = self._stratum_lane_boundaries(out, stratify_col, strata, collapsed_factors) # 10. Statistics per stratum: limits vary per phase chart_statistics = { @@ -2311,29 +2347,10 @@ def _calculate_mr_chart_stratified( # noqa: C901 # For X: center=mean, limits based on mean ± E2*mR # For mR: center=mR, limits based on 0 to D4*mR if mr_spec.center_source == 'mean': - lims = grouped.apply( - lambda row: calculate_limits( - mean=row['mean'], - sd=0, - N=0, - mR=(row['mR'] if pd.notna(row['mR']) else 0.0), - limits_type=mr_spec.limits_type, - round_to=spec.round_to, - ), - axis=1, - ) + # A one-row stratum has no moving range: its X limits collapse to the mean. + lims = _mr_family_limits(mr_spec, mean=grouped['mean'], mR=grouped['mR'].fillna(0.0)) else: # center_source == 'mR' - lims = grouped.apply( - lambda row: calculate_limits( - mean=0, - sd=0, - N=0, - mR=row['mR'], - limits_type=mr_spec.limits_type, - round_to=spec.round_to, - ), - axis=1, - ) + lims = _mr_family_limits(mr_spec, mean=0, mR=grouped['mR']) # Normalize/join lpl/upl grouped = _join_limits(grouped, lims) @@ -2363,24 +2380,10 @@ def _calculate_mr_chart_stratified( # noqa: C901 strata = grouped[stratify_col].tolist() # Lane boundaries per stratum - lane_boundaries = {} - if collapsed_factors: - for stratum in strata: - stratum_data = out[out[stratify_col] == stratum].reset_index(drop=True) - boundaries = self._calculate_lane_boundaries(stratum_data, collapsed_factors) - if boundaries: - lane_boundaries[stratum] = boundaries + lane_boundaries = self._stratum_lane_boundaries(out, stratify_col, strata, collapsed_factors) # Build statistics nested by stratum - chart_statistics = {} - for stratum in strata: - row = grouped[grouped[stratify_col] == stratum].iloc[0] - chart_statistics[stratum] = { - 'N': _mr_family_subgroup_n(mr_spec.chart_type), - 'center': round(row['center'], spec.round_to), - 'lpl': round(row['lpl'], spec.round_to), - 'upl': round(row['upl'], spec.round_to), - } + chart_statistics = _mr_stratum_statistics(mr_spec, grouped, stratify_col, 'lpl', 'upl', spec.round_to) # Format output extra_cols = [c for c in stratify_by if c != spec.rsg_var_name and c != spec.time_var] diff --git a/processbehavior/analysis_result.py b/processbehavior/analysis_result.py index 9bf8619..4a86b39 100644 --- a/processbehavior/analysis_result.py +++ b/processbehavior/analysis_result.py @@ -1039,6 +1039,10 @@ def detect_signals( ``signals.is_partial`` is True and ``signals.summary`` says which rules were evaluated. + On a phased chart (``phased=True``) each phase has its own centre line, so + the default rules are Rule 1 only, judged against each phase's limits. An + explicit rule list still runs as given, against each row's own centre. + Access violations: >>> signals.by_rule['rule_2'] # Rule 2 violations diff --git a/processbehavior/focused_result.py b/processbehavior/focused_result.py index 84085d2..1bc4982 100644 --- a/processbehavior/focused_result.py +++ b/processbehavior/focused_result.py @@ -40,14 +40,16 @@ def focus_on(result: AnalysisResult, stratum: str) -> AnalysisResult: ProcessBehaviorError If the chart lacks both stratify_col metadata and an rsg column. """ - if not result.strata: + # `strata` is a property that rebuilds the list on every read: read it once. + strata = result.strata + if not strata: raise ValidationError( "Cannot focus: this result is not stratified. " "Use result.strata to check available subgroups." ) - if stratum not in result.strata: + if stratum not in strata: raise ValidationError( - f"Stratum '{stratum}' not found. Available strata: {result.strata}" + f"Stratum '{stratum}' not found. Available strata: {strata}" ) focused_charts: dict[str, dict[str, Any]] = {} diff --git a/processbehavior/plotting/lane_boundaries.py b/processbehavior/plotting/lane_boundaries.py index 0d0d5d6..9edef09 100644 --- a/processbehavior/plotting/lane_boundaries.py +++ b/processbehavior/plotting/lane_boundaries.py @@ -10,6 +10,12 @@ from .themes import ChartTheme +# More dividers than this can't be told apart at any chart width: they merge into a solid band, +# and the browser still has to draw every one. Above it none are drawn. A chart with that many +# lanes (one per condition: 6,032 for a 6,033-gene study) is read from its points and limits. +MAX_LANE_DIVIDERS = 1_000 + + def add_lane_boundaries( fig: go.Figure, lane_boundaries: list[dict] | None, @@ -25,6 +31,7 @@ def add_lane_boundaries( Unified function for both single charts and faceted subplots. Lane boundaries show where collapsed factors change within the chart, helping distinguish groups of observations from different factor levels. + Nothing is drawn above ``MAX_LANE_DIVIDERS`` boundaries. Parameters ---------- @@ -43,7 +50,7 @@ def add_lane_boundaries( show_labels : bool, default True Whether to show factor labels at boundary positions """ - if not lane_boundaries: + if not lane_boundaries or len(lane_boundaries) > MAX_LANE_DIVIDERS: return y_min, y_max = y_range @@ -52,27 +59,36 @@ def add_lane_boundaries( if row is not None and col is not None: subplot_kwargs = {'row': row, 'col': col} + line = dict(color=theme.lane_boundary_color, dash=theme.lane_boundary_dash, width=theme.lane_boundary_width) + font = dict(size=theme.lane_boundary_annotation_size, color=theme.lane_boundary_color) + shapes, annotations = [], [] for boundary in lane_boundaries: x_pos = boundary['position'] label = boundary.get('label', '') - - fig.add_shape( - type='line', - x0=x_pos, - x1=x_pos, - y0=y_min, - y1=y_max, - line=dict(color=theme.lane_boundary_color, dash=theme.lane_boundary_dash, width=theme.lane_boundary_width), - **subplot_kwargs, - ) - + shapes.append(dict(type='line', x0=x_pos, x1=x_pos, y0=y_min, y1=y_max, line=line)) if show_labels and label: - fig.add_annotation( - x=x_pos, - y=y_max, - text=label, - showarrow=False, - yanchor='bottom', - font=dict(size=theme.lane_boundary_annotation_size, color=theme.lane_boundary_color), - **subplot_kwargs, + annotations.append( + dict(x=x_pos, y=y_max, text=label, showarrow=False, yanchor='bottom', font=font) ) + + _add_all(fig, 'shapes', fig.add_shape, shapes, subplot_kwargs) + _add_all(fig, 'annotations', fig.add_annotation, annotations, subplot_kwargs) + + +def _add_all(fig: go.Figure, prop: str, add_one, items: list[dict], subplot_kwargs: dict) -> None: + """Add layout shapes or annotations in one assignment instead of one call each. + + ``fig.add_shape`` re-validates every shape already on the figure, so adding N one at a + time costs O(N^2): 800 dividers took 27s, and 6,032 would take about 25 minutes. The first + item goes through ``add_one`` so plotly resolves the subplot's axis references for + ``row``/``col``. The rest copy those references and are appended in one assignment. The + figure is the same as one built call by call. + """ + if not items: + return + add_one(**items[0], **subplot_kwargs) + if len(items) == 1: + return + first = fig.layout[prop][-1] + refs = {'xref': first.xref, 'yref': first.yref} + fig.layout[prop] = fig.layout[prop] + tuple({**item, **refs} for item in items[1:]) diff --git a/processbehavior/result_signals.py b/processbehavior/result_signals.py index b561ab1..37f247d 100644 --- a/processbehavior/result_signals.py +++ b/processbehavior/result_signals.py @@ -146,4 +146,6 @@ def _detect_for_chart( value_col=value_col, chart_name=chart_name, chart_type=chart_type, + # Phased charts mark run rules not applicable; until now only the renderer read it. + run_rules_applicable=metadata.get("run_rules_applicable", True), ) diff --git a/processbehavior/signals/detector.py b/processbehavior/signals/detector.py index 414f056..8bb2753 100644 --- a/processbehavior/signals/detector.py +++ b/processbehavior/signals/detector.py @@ -91,6 +91,7 @@ def detect( value_col: str = 'mean', chart_name: str = 'Chart', chart_type: str = 'Xbar', + run_rules_applicable: bool = True, ) -> SignalResult: """ Detect signals in control chart data. @@ -109,6 +110,11 @@ def detect( Name of chart for reporting chart_type : str, default 'Xbar' Chart type ('Xbar', 'S', 'X', 'mR') - determines applicable rules + run_rules_applicable : bool, default True + False for charts whose centre moves along the series (phased limits, + ``metadata['run_rules_applicable']``). The chart-type default rules + narrow to Rule 1; an explicit rule list runs as given, against each + row's own centre. Returns ------- @@ -126,6 +132,13 @@ def detect( # Get applicable rules for this chart type applicable_rules = config.get_rules_for_chart(chart_type) + # Phased limits give each phase its own centre, so the chart marks run rules + # 2-8 not applicable: they are not part of the default set, rather than + # skipped (a skip would mark the evaluation partial). An explicit rule list is + # the analyst's choice and runs as given. + if not run_rules_applicable and config.enabled_rules == 'default': + applicable_rules = [r for r in applicable_rules if r == 'rule_1'] + # Log rule selection if logger.isEnabledFor(logging.INFO): all_rules = [f'rule_{i}' for i in range(1, 9)] @@ -291,7 +304,7 @@ def _build_result( 'rule_number': int(rule_name.split('_')[1]), 'description': self.RULE_DESCRIPTIONS[rule_name], 'value': data.loc[idx, value_col], - 'center': stats['center'], + 'center': self._row_center(data, stats, idx), 'upl': ucl, 'lpl': lcl, } @@ -309,6 +322,22 @@ def _limits_vary(self, stats: dict) -> bool: """Check if control limits vary (per-row limits).""" return bool(stats.get('limits_vary')) + @staticmethod + def _center(data: pd.DataFrame, stats: dict): + """The centre line: the scalar, or the per-row column when there is none. + + Phased charts give each phase its own centre, so their statistics carry + ``center=None`` and the per-phase value lives on each row. + """ + if stats.get('center') is None and 'center' in data.columns: + return data['center'] + return stats['center'] + + def _row_center(self, data: pd.DataFrame, stats: dict, idx): + """The centre line at one row.""" + center = self._center(data, stats) + return center.loc[idx] if isinstance(center, pd.Series) else center + def _calculate_per_row_zones(self, data: pd.DataFrame, stats: dict, config: SignalConfig) -> pd.DataFrame: """ Calculate per-row zone boundaries for varying process limits. @@ -334,7 +363,7 @@ def _calculate_per_row_zones(self, data: pd.DataFrame, stats: dict, config: Sign if 'upl' not in data.columns or 'lpl' not in data.columns: raise ValueError("Cannot calculate per-row zones: data missing 'upl' or 'lpl' columns") - center = stats['center'] + center = self._center(data, stats) zone_def = config.zone_definition # Calculate sigma per row diff --git a/processbehavior/signals/detectors.py b/processbehavior/signals/detectors.py index 38e18e0..2c2431b 100644 --- a/processbehavior/signals/detectors.py +++ b/processbehavior/signals/detectors.py @@ -173,6 +173,9 @@ def detect_run(data: pd.DataFrame, stats: dict, value_col: str, length: int = 8) """ values = data[value_col] center = stats.get('center', stats.get('Mean')) # Try 'center' first, fallback to 'Mean' + if center is None and 'center' in data.columns: + # Phased charts: each phase has its own centre, carried on the rows. + center = data['center'] # Above or below center above = (values > center).astype(int) diff --git a/processbehavior/spc_constants.py b/processbehavior/spc_constants.py index 459e03f..1b4323a 100644 --- a/processbehavior/spc_constants.py +++ b/processbehavior/spc_constants.py @@ -389,11 +389,12 @@ def calculate_limits_vectorized( Notes ----- **Conversion is deliberately partial.** The four primary Xbar/S limit computations in - ``analysis.py`` use this; the phase-segment and per-stratum sites still call the scalar - form. Those iterate over phases or strata (a handful of rows), not subgroups, so they - were not hot when measured — converting them would carry risk for no gain. If a future - profile shows otherwise, they are the next candidates; do not assume the absence of a - call here means a site was checked and rejected on correctness grounds. + ``analysis.py`` use this, and so do the stratified X/mR paths (through + ``_mr_family_limits``): those once iterated over a handful of strata, until a study + with one stratum per row (300,000 of them) made the scalar calls cost 30s per chart. + The ungrouped phase-segment sites still call the scalar form; they iterate over phases, + not subgroups or strata. Do not assume the absence of a call here means a site was + checked and rejected on correctness grounds. """ index = next((arg.index for arg in (mean, sd, N, mR) if isinstance(arg, pd.Series)), None) diff --git a/processbehavior/study.py b/processbehavior/study.py index fda6421..2f4a935 100644 --- a/processbehavior/study.py +++ b/processbehavior/study.py @@ -75,6 +75,17 @@ def _base_residual_code(value: str) -> str: return upper[2:] if upper.startswith('RCR') else upper +def _normalise_value(value: str | None) -> str | None: + """'response' (any case) -> None, the response's own spelling; residuals unchanged. + + Residuals are told apart by their leading 'R', which 'RESPONSE' also has, so the + documented ``value='response'`` must become None before anything tests the prefix. + """ + if isinstance(value, str) and value.lower() == 'response': + return None + return value + + @dataclass class DesignReport: """ @@ -1521,6 +1532,8 @@ def why_not(self, chart: str, value: str | None = None) -> str: >>> study.why_not('Xbar', value='R2') "'Xbar' with value='R2' is not valid for ADS 1.\\nValid charts for R2: S, X" """ + value = _normalise_value(value) + # Handle old fused syntax: why_not('R5_Xbar') -> why_not('Xbar', value='R5') if '_' in chart and value is None: parts = chart.split('_', 1) @@ -2018,6 +2031,8 @@ def execute( ) by = stratify_by + value = _normalise_value(value) + # 0. Resolve calibration (label string -> attached Calibration object). calibration_obj = self._resolve_calibration(calibration) @@ -2114,6 +2129,7 @@ def supports_calibration( 'parameter.' ) by = stratify_by + value = _normalise_value(value) from .analysis import Analysis from .formulation_spec import ChartRequest @@ -2396,7 +2412,10 @@ def _validate_recentered(self, value: str | None) -> None: ) recenterable = set(RESIDUAL_CODES) if value is None: - raise ValidationError('recentered=True requires a residual value (R1-R6), got value=None') + raise ValidationError( + 'recentered=True requires a residual value (R1-R6); the response ' + '(value=None or "response") is already on the measurement scale.' + ) if value.upper() not in recenterable: raise ValidationError(f"recentered=True requires a residual value (R1-R6), got '{value}'") diff --git a/scripts/benchmark.py b/scripts/benchmark.py index a0b0e50..7c138f8 100644 --- a/scripts/benchmark.py +++ b/scripts/benchmark.py @@ -2,8 +2,9 @@ """Performance benchmark + drift check for processbehavior. Runs a fixed set of scenarios — design {SDS-1 exact-R2, SDS-2 MA2 worst-case} × size -{10K, 100K, 1M} × op {init, formulate, execute, execute-stratified} — prints a report table, writes -``benchmarks/last_run.json``, and compares against the committed ``benchmarks/baseline.json``. +{10K, 100K, 1M} × op {init, formulate, execute, execute-stratified, and for SDS-2 one-row +strata} — prints a report table, writes ``benchmarks/last_run.json``, and compares against the +committed ``benchmarks/baseline.json``. The baseline is *same-machine* (absolute times are hardware-specific), so a regression flag is meaningful when you re-run on the same box. Cross-machine, rely on the shape checks in @@ -109,13 +110,24 @@ def run_scenario(sds: int, target_rows: int) -> list[dict]: peak = _peak_mb(lambda: ProcessBehavior(df).formulate(**spec)) common = {"sds": sds, "size": target_rows, "rows": rows} - return [ + results = [ {**common, "op": "init", "seconds": t_init}, {**common, "op": "formulate", "seconds": t_form}, {**common, "op": f"execute_{exec_chart.lower()}", "seconds": t_exec}, {**common, "op": "execute_x_stratified", "seconds": t_strat}, - {**common, "op": "formulate_peak", "seconds": None, "peak_mb": peak, "input_mb": input_mb}, ] + if sds == 2: + # One row per stratum: the time column taken as a factor makes every row its own + # condition. The scenario above has at most a few thousand strata; this one has as + # many as rows. The stratified X/mR path was quadratic in strata until 2026-10-12: + # about 80 minutes per chart at 307,683 one-row strata, unnoticed because nothing + # measured this shape. + one_row = [*DESIGNS[sds]["factors"], "time"] + one_row_study = pdf.formulate(response="y", factors=one_row) + t_one_row, _ = _time(lambda: one_row_study.execute(chart="X", by=one_row, companion=True), repeat) + results.append({**common, "op": "execute_x_one_row_strata", "seconds": t_one_row}) + results.append({**common, "op": "formulate_peak", "seconds": None, "peak_mb": peak, "input_mb": input_mb}) + return results def run_benchmarks(quick: bool = False) -> list[dict]: diff --git a/tests/test_execute_validation.py b/tests/test_execute_validation.py index 054d113..a54591c 100644 --- a/tests/test_execute_validation.py +++ b/tests/test_execute_validation.py @@ -218,3 +218,54 @@ def test_execute_with_accessor(self, sds1_study): """End-to-end: study.execute(chart=..., value=study.residuals.R5).""" result = sds1_study.execute(chart=sds1_study.charts.Xbar, value=sds1_study.residuals.R5) assert result is not None + + +# ============================================================================= +# value='response' is the response (#142) +# ============================================================================= + + +class TestValueResponse: + """The documented value='response' charts the response, like value=None. + + Residuals are told apart by a leading 'R', which 'RESPONSE' also has, so the + documented spelling was routed down the residual path and refused. + """ + + @pytest.mark.parametrize( + 'kwargs', + [ + {'chart': 'Xbar'}, + {'chart': 'S'}, + {'chart': 'Xbar', 'companion': True}, + {'chart': 'X', 'by': []}, + {'chart': 'mR', 'by': []}, + {'chart': 'Histogram'}, + ], + ids=['Xbar', 'S', 'Xbar+S', 'X', 'mR', 'Histogram'], + ) + @pytest.mark.parametrize('spelling', ['response', 'Response', 'RESPONSE']) + def test_same_result_as_none(self, sds1_study, kwargs, spelling): + from pandas.testing import assert_frame_equal + + expected = sds1_study.execute(**kwargs) + actual = sds1_study.execute(value=spelling, **kwargs) + assert set(actual.charts) == set(expected.charts) + for name in expected.charts: + assert_frame_equal(actual.get_chart(name), expected.get_chart(name)) + assert actual.get_statistics(name) == expected.get_statistics(name) + assert actual.charts[name]['metadata'] == expected.charts[name]['metadata'] + + def test_why_not_agrees(self, sds1_study): + for chart in ('Xbar', 'S', 'X', 'mR'): + assert sds1_study.why_not(chart, value='response') == sds1_study.why_not(chart) + + def test_supports_calibration_agrees(self, sds1_study): + for kwargs in ({'chart': 'X', 'by': []}, {'chart': 'Xbar'}, {'chart': 'X', 'by': ['factor 1']}): + assert sds1_study.supports_calibration(value='response', **kwargs) == ( + sds1_study.supports_calibration(**kwargs) + ) + + def test_recentered_response_still_refused(self, sds1_study): + with pytest.raises(ValidationError, match='recentered.*requires.*residual'): + sds1_study.execute(chart='Xbar', value='response', recentered=True) diff --git a/tests/test_performance.py b/tests/test_performance.py index c02e87f..ad113f0 100644 --- a/tests/test_performance.py +++ b/tests/test_performance.py @@ -592,6 +592,14 @@ def test_ma2_scales_near_linearly(self, benchmark_results): print(f'\nSDS-2 MA2 per-row: {small:.2f}us (10K) -> {large:.2f}us (100K), ratio {large / small:.2f}') assert large <= small * self.MAX_RATIO, f'MA2 per-row grew {large / small:.1f}x (>{self.MAX_RATIO}x)' + def test_one_row_strata_scale_linearly(self, benchmark_results): + """Stratified X/mR with every row its own stratum. A per-stratum table scan makes + this quadratic: 10x the rows took ~100x the time before 2026-10-12.""" + small = _per_row_us(benchmark_results, 2, 'execute_x_one_row_strata', 10_000) + large = _per_row_us(benchmark_results, 2, 'execute_x_one_row_strata', 100_000) + print(f'\nOne-row strata per-row: {small:.2f}us (10K) -> {large:.2f}us (100K), ratio {large / small:.2f}') + assert large <= small * self.MAX_RATIO, f'one-row strata per-row grew {large / small:.1f}x (>{self.MAX_RATIO}x)' + @pytest.mark.slow @pytest.mark.benchmark diff --git a/tests/test_phased_limits.py b/tests/test_phased_limits.py index ff451cf..d4c5cf6 100644 --- a/tests/test_phased_limits.py +++ b/tests/test_phased_limits.py @@ -412,3 +412,77 @@ def test_phased_residual_limits_vary(self, sds2_study): result = sds2_study.execute(chart='X', by=[], phased=True, value='R2') data = result.get_chart('X') assert data['center'].nunique() > 1 or data['upl'].nunique() > 1 + + +# ============================================================================= +# Phased + signal detection +# ============================================================================= + + +class TestPhasedSignals: + """detect_signals() on phased charts: Rule 1 against each phase's own limits. + + Phased statistics carry ``center=None`` (each phase has its own centre) and the + chart marks run rules not applicable. detect_signals() used to read the scalar + centre and raise TypeError on every phased chart. + """ + + def test_phased_x_rule_1_against_phase_limits(self, sds1_study): + result = sds1_study.execute(chart='X', by=[], phased=True) + signals = result.detect_signals(chart='X') + assert signals.rules_evaluated == ['rule_1'] + assert signals.rules_skipped == {} + assert not signals.is_partial + data = result.get_chart('X') + beyond = set(data.index[data['beyond_limits'] != 0]) + assert beyond, 'fixture should have a point beyond its phase limits' + assert signals.flagged_observations == beyond + + def test_violation_records_the_phase_centre(self, sds1_study): + result = sds1_study.execute(chart='X', by=[], phased=True) + data = result.get_chart('X') + violations = result.detect_signals(chart='X').to_dataframe() + assert violations['center'].notna().all() + for _, row in violations.iterrows(): + assert row['center'] == data.loc[row['obs_id'], 'center'] + assert row['upl'] == data.loc[row['obs_id'], 'upl'] + + def test_phased_companion_both_charts(self, sds1_study): + result = sds1_study.execute(chart='X', by=[], phased=True, companion=True) + signals = result.detect_signals() + assert set(signals) == {'X', 'mR'} + for name, sig in signals.items(): + data = result.get_chart(name) + assert sig.rules_evaluated == ['rule_1'] + assert sig.flagged_observations == set(data.index[data['beyond_limits'] != 0]) + + def test_phased_stratified_focus(self, sds1_study): + result = sds1_study.execute(chart='X', by=['factor 1'], phased=True, companion=True) + for stratum in result.strata: + signals = result.focus(stratum).detect_signals() + assert {sig.rules_evaluated[0] for sig in signals.values()} == {'rule_1'} + + def test_explicit_rules_run_against_each_phase_centre(self): + """An explicit rule list runs as given; run and zone rules read the row's centre.""" + from processbehavior.signals import SignalConfig, SignalDetector + + # Phase 1 sits above its centre 0, phase 2 below its centre 10: one run each. + # Phase 2's last two points fall in its lower zone A (7 to 8, sigma 1). + data = pd.DataFrame({ + 'x': [1.0] * 8 + [9.0] * 6 + [7.5, 7.5], + 'center': [0.0] * 8 + [10.0] * 8, + 'upl': [3.0] * 8 + [13.0] * 8, + 'lpl': [-3.0] * 8 + [7.0] * 8, + }) + stats = {'N': None, 'center': None, 'lpl': None, 'upl': None, 'limits_vary': True} + signals = SignalDetector().detect( + data, stats, SignalConfig(enabled_rules=['rule_2', 'rule_4'], min_observations=1), + value_col='x', chart_type='X', run_rules_applicable=False, + ) + assert signals.rules_evaluated == ['rule_2', 'rule_4'] + assert set(signals.by_rule['rule_4']['obs_id']) == {7, 15} + assert set(signals.by_rule['rule_2']['obs_id']) == {15} + + def test_unphased_x_keeps_every_rule(self, sds1_study): + signals = sds1_study.execute(chart='X', by=[]).detect_signals(chart='X') + assert signals.rules_evaluated == [f'rule_{i}' for i in range(1, 9)] diff --git a/tests/test_plotting_decorations.py b/tests/test_plotting_decorations.py index 68364b6..0583b5b 100644 --- a/tests/test_plotting_decorations.py +++ b/tests/test_plotting_decorations.py @@ -12,7 +12,7 @@ from plotly.subplots import make_subplots from processbehavior.plotting import get_theme -from processbehavior.plotting.lane_boundaries import add_lane_boundaries +from processbehavior.plotting.lane_boundaries import MAX_LANE_DIVIDERS, add_lane_boundaries from processbehavior.plotting.limits import ( add_stepped_limit_line, build_stepped_coordinates, @@ -224,6 +224,47 @@ def test_boundary_without_label(self, theme): assert len(fig.layout.shapes) == 1 assert len(fig.layout.annotations) == 0 + @staticmethod + def _one_call_each(fig, boundaries, y_range, theme, row=None, col=None, show_labels=True): + """The figure add_lane_boundaries built before batching: one add_shape per divider.""" + y_min, y_max = y_range + kw = {'row': row, 'col': col} if row is not None and col is not None else {} + for b in boundaries: + fig.add_shape(type='line', x0=b['position'], x1=b['position'], y0=y_min, y1=y_max, + line=dict(color=theme.lane_boundary_color, dash=theme.lane_boundary_dash, + width=theme.lane_boundary_width), **kw) + if show_labels and b.get('label', ''): + fig.add_annotation(x=b['position'], y=y_max, text=b['label'], showarrow=False, yanchor='bottom', + font=dict(size=theme.lane_boundary_annotation_size, + color=theme.lane_boundary_color), **kw) + + @pytest.mark.parametrize(('make', 'subplot'), [ + (go.Figure, {}), + (lambda: make_subplots(rows=2, cols=2), {'row': 2, 'col': 1}), + ]) + @pytest.mark.parametrize('show_labels', [True, False]) + def test_batched_figure_matches_one_call_each(self, theme, make, subplot, show_labels): + """Adding the dividers in one assignment builds the same figure, references included.""" + boundaries = [{'position': i * 3, 'label': f'L{i}' if i % 4 else ''} for i in range(1, 40)] + expected, actual = make(), make() + for fig in (expected, actual): + fig.add_shape(type='rect', x0=0, x1=1, y0=0, y1=1) # a shape already on the figure + self._one_call_each(expected, boundaries, (-1, 5), theme, show_labels=show_labels, **subplot) + add_lane_boundaries(actual, boundaries, (-1, 5), theme, show_labels=show_labels, **subplot) + assert actual.to_dict() == expected.to_dict() + + def test_draws_up_to_the_divider_limit(self, theme): + fig = go.Figure() + add_lane_boundaries(fig, [{'position': i} for i in range(MAX_LANE_DIVIDERS)], (0, 1), theme) + assert len(fig.layout.shapes) == MAX_LANE_DIVIDERS + + def test_draws_none_above_the_divider_limit(self, theme): + """Past the limit dividers merge into a solid band: one lane per gene of a 6,033-gene study.""" + fig = go.Figure() + add_lane_boundaries(fig, [{'position': i, 'label': 'g'} for i in range(6_032)], (0, 1), theme) + assert len(fig.layout.shapes) == 0 + assert len(fig.layout.annotations) == 0 + # ========================================================================= # Stats Box diff --git a/tests/test_stratified_insufficiency.py b/tests/test_stratified_insufficiency.py index d9691b2..4e9c33e 100644 --- a/tests/test_stratified_insufficiency.py +++ b/tests/test_stratified_insufficiency.py @@ -184,3 +184,28 @@ def test_sufficient_data_publishes_every_stratum(df, response, chart): assert _meta(result, chart)['insufficient_strata'] is None for stratum in result.strata: assert len(result.focus(stratum).get_chart(chart)) > 0 + + +# --------------------------------------------------------------------------- +# Every row its own stratum (X/mR) +# --------------------------------------------------------------------------- + + +def test_one_row_strata_publish_nothing_and_report_every_stratum(df): + """Every row its own condition, as in a 307,683-row study that took 80 minutes per chart. + + Each stratum holds one reading, so none has a moving range and none is published, but + every one is reported. The time column taken as a factor does it on the reference data. + The cost of this shape is guarded in tests/test_performance.py (one-row strata). + """ + by = [*FACTORS, TIME] + study = pb.formulate(df, response='PM SDS 2', factors=by) + result = study.execute(chart='X', by=by, companion=True) + + rows = int(df['PM SDS 2'].notna().sum()) + assert result.strata == [] + for chart in ('X', 'mR'): + meta = _meta(result, chart) + assert len(meta['insufficient_strata']) == rows + assert len(set(meta['insufficient_strata'])) == rows + assert meta['lane_boundaries'] is None