diff --git a/processbehavior/analysis.py b/processbehavior/analysis.py index b04585e..f229664 100644 --- a/processbehavior/analysis.py +++ b/processbehavior/analysis.py @@ -211,6 +211,38 @@ def _stratum_row_positions(df: pd.DataFrame, stratify_col: str) -> dict: return df.groupby(stratify_col, sort=False, observed=True).indices +@dataclass(frozen=True) +class _GroupFit: + """What one group of Xbar/S subgroups (a chart, a stratum or a phase) needs for its limits. + + Built by :meth:`Analysis._fit_xbar_group` from an aggregated subgroup frame; consumed by + :meth:`Analysis._apply_xbar_limits` and, for the companion S chart, by + :meth:`Analysis._apply_s_limits`. + """ + + out: pd.DataFrame # the charted subgroups (n >= 1); N overwritten under n_mode='average' + fit: pd.Series # mask over ``out``: the subgroups with n >= 2 that sigma-hat is estimated from + s_bar: float # mean S over the fitted subgroups + n_to_use: str # 'N' when the Xbar chart's subgroup sizes are constant, else 'n' + s_n_to_use: str # the same decision for the S chart, which has no one-reading subgroups + n_max: int + n_avg: float | None # mean fitted n under n_mode='average', else None + + +def _lane_labels(df: pd.DataFrame, collapsed_vars: list[str]) -> np.ndarray: + """One display label per row of ``df``: the collapsed columns as ``astype(str)``, joined with ``'_'``. + + Column-wise rather than ``.agg('_'.join, axis=1)``, which ran once per row (1.05s at 1M). + Deliberately NOT encode_rsg_series: that renders a datetime as str(Timestamp) + ('2024-01-01 00:00:00') where DataFrame.astype(str) gives '2024-01-01'. These are display + labels for a chart divider or a phase, not stratum identity, so they keep the astype + rendering — do not "unify" the two. + """ + parts = [df[col].astype(str) for col in collapsed_vars] + combined = parts[0] if len(parts) == 1 else reduce(lambda left, right: left + '_' + right, parts) + return combined.to_numpy() + + 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) @@ -1063,19 +1095,11 @@ def _calculate_lane_boundaries(self, df: pd.DataFrame, collapsed_vars: list[str] if not collapsed_vars or df.empty: return [] - # Combine column-wise rather than with `.agg('_'.join, axis=1)`, which ran once per - # row (1.05s at 1M). Deliberately NOT encode_rsg_series: that renders a datetime as - # str(Timestamp) ('2024-01-01 00:00:00') where DataFrame.astype(str) gives - # '2024-01-01'. These are display labels for a chart divider, not stratum identity, - # so they keep the astype rendering — do not "unify" the two. - parts = [df[col].astype(str) for col in collapsed_vars] - combined = parts[0] if len(parts) == 1 else reduce(lambda left, right: left + '_' + right, parts) - # Positional comparison against the previous row. The old form built a boolean mask, # took df.index[mask], dropped the first entry, then called df.index.get_loc() per # boundary — and get_loc is not O(1) on these frames (~1.3ms a call, ~1.8s total). # Comparing the array to itself shifted gives the 0-based positions directly. - labels = combined.to_numpy() + labels = _lane_labels(df, collapsed_vars) changed = np.empty(len(labels), dtype=bool) changed[0] = False # the first row is a start, not a boundary changed[1:] = labels[1:] != labels[:-1] @@ -1187,6 +1211,117 @@ def _xbar_limits(self, out: pd.DataFrame, s_bar: float, n_to_use: str) -> pd.Dat sigma_multiplier=self.request.n_sigma, ) + def _aggregate_subgroups( + self, + df: pd.DataFrame, + groupby_cols: list[str], + value_col: str, + limits_col: str, + ybar_col: str | None = None, + ) -> pd.DataFrame: + """Aggregate ``df`` into Xbar/S subgroups, one row per ``groupby_cols`` combination. + + Columns: the grouping columns; ``s`` and ``n`` from ``limits_col`` (for an effect + residual that is R2, so a reading without one leaves the subgroup one smaller, as VAS + treats it); ``xbar`` from ``value_col``; ``s_value`` when the two columns differ; and + ``N``, the largest ``n``. With ``ybar_col`` the residual calculator's cached cell mean + stands in for the subgroup mean. + """ + agg_dict = { + 's': pd.NamedAgg(column=limits_col, aggfunc='std'), + 'mean': pd.NamedAgg(column=value_col, aggfunc='mean'), + # N_r counts the values S_r is computed from: a residual missing on a row (R2's + # first value) leaves that subgroup one smaller, as VAS treats it. + 'n': pd.NamedAgg(column=limits_col, aggfunc='count'), + } + if limits_col != value_col: + agg_dict['s_value'] = pd.NamedAgg(column=value_col, aggfunc='std') + out = df.groupby(groupby_cols, as_index=False, observed=True).agg(**agg_dict) + + if ybar_col and ybar_col in df.columns: + # Get unique Ybar values per group + ybar_by_group = df.groupby(groupby_cols, observed=True)[ybar_col].first().reset_index() + out = out.merge(ybar_by_group, on=groupby_cols, how='left') + out['xbar'] = out[ybar_col] + out = out.drop(columns=[ybar_col, 'mean'], errors='ignore') + else: + out = out.rename(columns={'mean': 'xbar'}) + + out['N'] = out['n'].max() + return out + + def _fit_xbar_group(self, out: pd.DataFrame) -> _GroupFit | None: + """Decide which subgroups of ``out`` are charted and which estimate sigma-hat. + + Every subgroup with a value is charted. A one-reading subgroup has no S, so sigma-hat + comes from the subgroups with n >= 2, and its limits are center ± 3·sigma-hat (10-12 + manual Eqs 11-16, 11-17 with N_r = 1), as VAS draws them; the companion S chart keeps + only n >= 2. n counts the limits column, so a reading without one (R2's first value) is + not charted, as in VAS. Returns None when no subgroup has n >= 2; the caller decides + whether that raises (one chart) or is recorded (a stratum, a phase). Under + ``n_mode='average'`` the ``N`` column becomes the mean fitted n and both charts take + the constant-N path (Eq 11-20). + """ + out = out[out['n'].ge(1)].copy() # n=0: every charted value missing + fit = out['n'].ge(2) + if not fit.any(): + return None + + s_bar = out.loc[fit, 's'].mean() + # Constant or variable subgroup sizes: on the Xbar chart, and on the S chart, which has + # no one-reading subgroups + n_to_use, n_max = self._determine_n_to_use(out) + s_n_to_use, _ = self._determine_n_to_use(out[fit]) + + n_avg = None + if self.request.n_mode == 'average': + n_avg = out.loc[fit, 'n'].mean() + out['N'] = n_avg # overwrite N column with average + n_to_use = s_n_to_use = 'N' # force constant-N path + return _GroupFit( + out=out, fit=fit, s_bar=s_bar, n_to_use=n_to_use, s_n_to_use=s_n_to_use, n_max=n_max, n_avg=n_avg + ) + + def _apply_xbar_limits(self, out: pd.DataFrame, center: float, group: _GroupFit) -> pd.DataFrame: + """Centre, natural process limits, beyond-limits flag and rounding for one Xbar group.""" + out['center'] = center + out[['lpl', 'upl']] = self._xbar_limits(out, group.s_bar, group.n_to_use) + out = self._add_beyond_limits_flag(out, value_col='xbar') + return out.round(self.spec.round_to) + + def _apply_s_limits(self, out: pd.DataFrame, s_bar: float, n_to_use: str) -> pd.DataFrame: + """Centre (S-bar), B3/B4 limits per row, beyond-limits flag and rounding for one S group.""" + out['center'] = s_bar + out[['lpl', 'upl']] = calculate_limits_vectorized( + 'S', + sd=out['center'], + N=out[n_to_use], + sigma_multiplier=self.request.n_sigma, + ) + out = self._add_beyond_limits_flag(out, value_col='s') + return out.round(self.spec.round_to) + + def _subgroup_label_column(self, out: pd.DataFrame, groupby_cols: list[str]) -> tuple[pd.DataFrame, str | None]: + """Name the subgroup label column of an Xbar/S frame, building it when it spans several columns. + + An explicit single ``by`` keeps the factor's own column (cast to str so Plotly treats + coded factors as categories); several ``by`` columns join into ``subgroup``; the + default grouping joins into ``group``; a single default column is used as is. + """ + by = self.request.by + if by is not None and len(by) > 0: + if len(by) == 1: + out[groupby_cols[0]] = out[groupby_cols[0]].astype(str) + return out, groupby_cols[0] + out['subgroup'] = out[groupby_cols].astype(str).agg('_'.join, axis=1) + return out, 'subgroup' + if len(groupby_cols) > 1: + out['group'] = out[groupby_cols].astype(str).agg('_'.join, axis=1) + return out, 'group' + if groupby_cols: + return out, groupby_cols[0] + return out, None + def _calculate_xbar( # noqa: C901 self, value_col: str | None = None, _return_intermediates: bool = False ) -> dict: @@ -1269,81 +1404,37 @@ def _calculate_xbar( # noqa: C901 out_dict['s_value'] = [df[value_col].std()] out = pd.DataFrame(out_dict) else: - # Group by specified columns - agg_dict = { - 's': pd.NamedAgg(column=_limits_col, aggfunc='std'), - 'mean': pd.NamedAgg(column=value_col, aggfunc='mean'), - # N_r counts the values S_r is computed from: a residual missing on a row (R2's - # first value) leaves that subgroup one smaller, as VAS treats it. - 'n': pd.NamedAgg(column=_limits_col, aggfunc='count'), - } - if _limits_col != value_col: - agg_dict['s_value'] = pd.NamedAgg(column=value_col, aggfunc='std') - out = df.groupby(groupby_cols, as_index=False, observed=True).agg(**agg_dict) - - # If we have pre-calculated Ybar, use it for xbar values - if ybar_col and ybar_col in df.columns: - # Get unique Ybar values per group - ybar_by_group = df.groupby(groupby_cols, observed=True)[ybar_col].first().reset_index() - out = out.merge(ybar_by_group, on=groupby_cols, how='left') - out['xbar'] = out[ybar_col] - out = out.drop(columns=[ybar_col, 'mean'], errors='ignore') - else: - # Rename columns for consistency: mean→xbar (values) - out = out.rename(columns={'mean': 'xbar'}) - - _N = out['n'].max() - out['N'] = _N - - # Every subgroup with a value is charted. A one-reading subgroup has no S, so sigma-hat - # comes from the subgroups with n >= 2, and its limits are center ± 3·sigma-hat (10-1 - # manual Eqs 11-16/17 with N_r = 1), as VAS draws them; the companion S chart keeps only - # n >= 2. n counts the limits column, so a reading without one (R2's first value) is not - # charted, as in VAS. Xbar stays in valid_charts on ADS 2 — the chart family is valid - # there (Xbar of a residual, or of a pooled subgroup, computes fine); it is only charting - # *the response* that has nothing to subgroup. See Study._response_pair_problem, which - # reports that in advance. - out = out[out['n'].ge(1)].copy() # n=0: every charted value missing - fit = out['n'].ge(2) - if not fit.any(): + out = self._aggregate_subgroups(df, groupby_cols, value_col, _limits_col, ybar_col) + + # Which subgroups are charted and which estimate sigma-hat: _fit_xbar_group. Xbar stays + # in valid_charts on ADS 2 — the chart family is valid there (Xbar of a residual, or of + # a pooled subgroup, computes fine); it is only charting *the response* that has nothing + # to subgroup. See Study._response_pair_problem, which reports that in advance. + group = self._fit_xbar_group(out) + if group is None: sds = self.ads._ads_result.sds if self.ads._ads_result else '?' _raise_no_replicated_subgroups('Xbar', sds) + out, fit = group.out, group.fit + _S, n_to_use, s_n_to_use = group.s_bar, group.n_to_use, group.s_n_to_use + n_max, n_avg = group.n_max, group.n_avg + n_mode = self.request.n_mode # Use Bishop VAS grand mean (mean of cell means on value_col) as center _Xbar = _Ybar - _S = out.loc[fit, 's'].mean() _N = out.loc[fit, 'n'].max() - if 'N' not in out.columns: - out['N'] = _N - - # Determine if subgroup sizes are constant or variable: on the Xbar chart, and on the - # S chart, which has no one-reading subgroups - n_to_use, n_max = self._determine_n_to_use(out) - s_n_to_use, _ = self._determine_n_to_use(out[fit]) - - # Override n_to_use if n_mode="average" - n_mode = self.request.n_mode - n_avg = None - if n_mode == 'average': - n_avg = out.loc[fit, 'n'].mean() - out['N'] = n_avg # overwrite N column with average - n_to_use = s_n_to_use = 'N' # force constant-N path # CALCULATE XBAR - xbar = out.copy() if self.calibration is not None: # Standards-given Xbar: center=mean, mean ± n_sigma·sigma/√N (per-row N). + xbar = out.copy() xbar[['center', 'lpl', 'upl']] = xbar.apply( lambda row: self._calibrated_row('Xbar', N=row[n_to_use]), axis=1, ) + xbar = self._add_beyond_limits_flag(xbar, value_col='xbar') + xbar = xbar.round(spec.round_to) else: - xbar['center'] = _Xbar # Add center column for Xbar chart - xbar[['lpl', 'upl']] = self._xbar_limits(xbar, _S, n_to_use) - - # Detect beyond limits signals - xbar = self._add_beyond_limits_flag(xbar, value_col='xbar') - xbar = xbar.round(spec.round_to) + xbar = self._apply_xbar_limits(out.copy(), _Xbar, group) _center_stat = xbar['center'].iloc[0] if self.calibration is not None else _Xbar statistics['center'] = round(_center_stat, spec.round_to) @@ -1357,27 +1448,7 @@ def _calculate_xbar( # noqa: C901 statistics['upl'] = None statistics['limits_vary'] = True - # Determine the grouping column for output - by = self.request.by - if by is not None and len(by) > 0: - # Explicit by list - if len(by) == 1: - # Single factor: preserve original column name for axis labeling. - # Cast to string so Plotly treats coded factors (e.g., 1,2,3) as categories. - xbar[groupby_cols[0]] = xbar[groupby_cols[0]].astype(str) - group_col = groupby_cols[0] - else: - xbar['subgroup'] = xbar[groupby_cols].astype(str).agg('_'.join, axis=1) - group_col = 'subgroup' - elif len(groupby_cols) > 1: - # by=None - create 'group' from multiple columns (original behavior) - xbar['group'] = xbar[groupby_cols].astype(str).agg('_'.join, axis=1) - group_col = 'group' - elif groupby_cols: - # Single groupby column - group_col = groupby_cols[0] - else: - group_col = None + xbar, group_col = self._subgroup_label_column(xbar, groupby_cols) # n is each subgroup's own count from the aggregation above: the size its limits # used. A pooled subgroup spans several cells, so no per-cell size could stand in (#143). @@ -1470,46 +1541,21 @@ def _calculate_xbar_stratified( # Group by time within this stratum # Note: we always compute mean directly from filtered data, NOT using # Ybar_t (which is the marginal time mean across ALL factors). - agg_dict = { - 's': pd.NamedAgg(column=_limits_col, aggfunc='std'), - 'xbar': pd.NamedAgg(column=value_col, aggfunc='mean'), - 'n': pd.NamedAgg(column=_limits_col, aggfunc='count'), - } - if _limits_col != value_col: - agg_dict['s_value'] = pd.NamedAgg(column=value_col, aggfunc='std') - out = sdf.groupby(groupby_cols, as_index=False, observed=True).agg(**agg_dict) - - out['N'] = out['n'].max() + out = self._aggregate_subgroups(sdf, groupby_cols, value_col, _limits_col) # Centre first: Bishop's unweighted mean of the subgroup means, every # subgroup counted once — a one-observation subgroup is still a cell. _Xbar = out['xbar'].mean() # Every subgroup with a value is charted; sigma-hat comes from those with n >= 2 - # (see _calculate_xbar) - out = out[out['n'].ge(1)].copy() - fit = out['n'].ge(2) - if not fit.any(): + # (see _fit_xbar_group) + group = self._fit_xbar_group(out) + if group is None: insufficient_strata.append(stratum) continue - - # Per-stratum statistics - _S = out.loc[fit, 's'].mean() - n_to_use, n_max = self._determine_n_to_use(out) - s_n_to_use, _ = self._determine_n_to_use(out[fit]) - - n_mode = self.request.n_mode - n_avg = None - if n_mode == 'average': - n_avg = out.loc[fit, 'n'].mean() - out['N'] = n_avg - n_to_use = s_n_to_use = 'N' - - out['center'] = _Xbar - out[['lpl', 'upl']] = self._xbar_limits(out, _S, n_to_use) - - out = self._add_beyond_limits_flag(out, value_col='xbar') - out = out.round(spec.round_to) + _S, n_to_use, s_n_to_use = group.s_bar, group.n_to_use, group.s_n_to_use + n_max, n_avg = group.n_max, group.n_avg + out = self._apply_xbar_limits(group.out, _Xbar, group) # Add stratum identifier out[stratify_col] = stratum @@ -1678,16 +1724,7 @@ def _calculate_s_stratified( if groupby_cols: out[groupby_cols[0]] = out[groupby_cols[0]].astype(str) - out['center'] = _S - out[['lpl', 'upl']] = calculate_limits_vectorized( - 'S', - sd=out['center'], - N=out[n_to_use], - sigma_multiplier=self.request.n_sigma, - ) - - out = self._add_beyond_limits_flag(out, value_col='s') - out = out.round(spec.round_to) + out = self._apply_s_limits(out, _S, n_to_use) out[stratify_col] = stratum all_s_frames.append(out) @@ -1833,22 +1870,9 @@ def _calculate_s( # noqa: C901 out['N'] = n_avg n_to_use = 'N' - # Determine group_col - by = self.request.by - if by is not None and len(by) > 0: - if len(by) == 1: - out[groupby_cols[0]] = out[groupby_cols[0]].astype(str) - group_col = groupby_cols[0] - else: - out['subgroup'] = out[groupby_cols].astype(str).agg('_'.join, axis=1) - group_col = 'subgroup' - elif len(groupby_cols) > 1: - out['group'] = out[groupby_cols].astype(str).agg('_'.join, axis=1) - group_col = 'group' - elif groupby_cols: - group_col = groupby_cols[0] - else: - group_col = None + # Label the subgroups. Xbar labelled its own copy; the shared intermediate frame is + # still unlabelled, so the companion path labels it here too. + out, group_col = self._subgroup_label_column(out, groupby_cols) # Calculate S chart (common path for both precomputed and independent) out['groups'] = out['n'].count() @@ -1860,41 +1884,15 @@ def _calculate_s( # noqa: C901 lambda row: self._calibrated_row('S', N=row[n_to_use]), axis=1, ) + out = self._add_beyond_limits_flag(out, value_col='s') + out = out.round(spec.round_to) else: - out['center'] = _S - # Add limits columns - out[['lpl', 'upl']] = calculate_limits_vectorized( - 'S', - sd=out['center'], - N=out[n_to_use], - sigma_multiplier=self.request.n_sigma, - ) - - # Detect beyond limits signals - out = self._add_beyond_limits_flag(out, value_col='s') - - # Apply grouping column transformation if needed and not already done - if _precomputed is None: - # Already handled in independent path above - pass - else: - # For precomputed path, apply group_col transformation - by = self.request.by - if by is not None and len(by) > 0 and 'subgroup' not in out.columns: - if len(by) == 1: - out[groupby_cols[0]] = out[groupby_cols[0]].astype(str) - group_col = groupby_cols[0] - else: - out['subgroup'] = out[groupby_cols].astype(str).agg('_'.join, axis=1) - group_col = 'subgroup' - elif len(groupby_cols) > 1 and group_col == 'group' and 'group' not in out.columns: - out['group'] = out[groupby_cols].astype(str).agg('_'.join, axis=1) + out = self._apply_s_limits(out, _S, n_to_use) # n: each subgroup's own count, as on Xbar (#143). cols_to_keep = [group_col, 'n', 's', 'center', 'lpl', 'upl', 'beyond_limits'] cols_to_keep = [c for c in cols_to_keep if c in out.columns] out = out[cols_to_keep] - out = out.round(spec.round_to) # Build statistics _S_stat = out['center'].iloc[0] if len(out) > 0 else None diff --git a/processbehavior/spc_constants.py b/processbehavior/spc_constants.py index adfe7af..3eb50ec 100644 --- a/processbehavior/spc_constants.py +++ b/processbehavior/spc_constants.py @@ -374,7 +374,8 @@ def calculate_limits_vectorized( ``DataFrame.apply(axis=1)``, each call constructing a ``pd.Series`` to carry two numbers. At 1M rows (~50,000 subgroups) that was ~2.5s per chart; this is ~0.001s. - ``mean``, ``sd``, ``N`` and ``mR`` accept scalars or array-likes and broadcast. + ``mean``, ``sd``, ``N``, ``mR`` and ``sigma`` accept scalars or array-likes and broadcast, + so one chart can carry a sigma-hat per row (a phased chart: one per phase). Subgroup-size constants (c4/b3/b4) are looked up once per *distinct* N and mapped, so an N of 1 still raises ``ValueError`` from :func:`c4` exactly as the scalar path does. @@ -396,7 +397,7 @@ def calculate_limits_vectorized( 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) + index = next((arg.index for arg in (mean, sd, N, mR, sigma) if isinstance(arg, pd.Series)), None) def _per_n(func, sizes, *args): """Evaluate a subgroup-size constant once per distinct N, then broadcast.""" @@ -413,7 +414,7 @@ def _per_n(func, sizes, *args): sizes = np.asarray(N) if sigma is not None: # Eqs 11-16/11-17: half-width = multiplier * sigma_hat / sqrt(N_r) - half = (sigma_multiplier * float(sigma)) / np.sqrt(sizes) + half = (sigma_multiplier * np.asarray(sigma, dtype=float)) / np.sqrt(sizes) else: # Wd = S / c4(n), then half-width = (multiplier * Wd) / sqrt(n) half = (sigma_multiplier * (np.asarray(sd) / _per_n(c4, sizes))) / np.sqrt(sizes) diff --git a/tests/test_xbar_s_helpers.py b/tests/test_xbar_s_helpers.py new file mode 100644 index 0000000..a47f13e --- /dev/null +++ b/tests/test_xbar_s_helpers.py @@ -0,0 +1,153 @@ +"""The Xbar/S helpers extracted for issue #144 (its first PR), tested on their own. + +The behaviour they carry is pinned end to end by the golden masters and by +``validation/xbar_s_sweep.py`` (identical fingerprints before and after the extraction). These +tests pin the pieces the phased Xbar/S chart will reuse: one label rule for lane dividers and +phases, a per-row sigma-hat in ``calculate_limits_vectorized``, and the per-group fit. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest +from conftest import detect_sds_for_test, make_spec + +from processbehavior.analysis import Analysis, _lane_labels +from processbehavior.datasets.synthetic import make_design +from processbehavior.formulation_spec import ChartRequest +from processbehavior.spc_constants import b3, b4, calculate_limits_vectorized + +SPEC = {'rsg_vars': ['factor 1', 'factor 2'], 'response_var': 'y', 'time_var': 'time', 'round_to': 3} + + +def _analysis(**request) -> Analysis: + df = make_design(1, K1=2, K2=2, T=3, n_min=2, n_max=3, seed=1) + sds = detect_sds_for_test(df, SPEC) + return Analysis(spec=make_spec(SPEC), request=ChartRequest(chart='Xbar', **request), sds=sds, df=df) + + +class TestLaneLabels: + def test_single_column_keeps_astype_rendering(self): + df = pd.DataFrame({'d': pd.to_datetime(['2024-01-01', '2024-01-02', '2024-01-02'])}) + # astype(str) renders a date without the time part; str(Timestamp) would not. + assert list(_lane_labels(df, ['d'])) == ['2024-01-01', '2024-01-02', '2024-01-02'] + + def test_several_columns_join_with_underscore(self): + df = pd.DataFrame({'a': [1, 1, 2], 'b': ['x', 'y', 'y']}) + labels = _lane_labels(df, ['a', 'b']) + assert isinstance(labels, np.ndarray) + assert list(labels) == ['1_x', '1_y', '2_y'] + + +class TestVectorizedSigmaPerRow: + def test_series_sigma_matches_the_formula_row_by_row(self): + mean = pd.Series([10.0, 10.0, 10.0], index=[5, 6, 7]) + sigma = pd.Series([1.0, 2.0, 0.5], index=[5, 6, 7]) + n = pd.Series([4, 4, 1], index=[5, 6, 7]) + limits = calculate_limits_vectorized('Xbar', mean=mean, sigma=sigma, N=n) + expected_half = 3 * sigma / np.sqrt(n) + pd.testing.assert_series_equal(limits['lpl'], mean - expected_half, check_names=False) + pd.testing.assert_series_equal(limits['upl'], mean + expected_half, check_names=False) + assert list(limits.index) == [5, 6, 7] + + def test_scalar_sigma_equals_a_constant_series(self): + n = pd.Series([2, 3, 5, 1]) + scalar = calculate_limits_vectorized('Xbar', mean=4.5, sigma=0.7, N=n) + series = calculate_limits_vectorized('Xbar', mean=4.5, sigma=pd.Series([0.7] * 4), N=n) + assert np.array_equal(scalar.to_numpy(), series.to_numpy()) + + def test_series_sigma_supplies_the_index_when_nothing_else_does(self): + sigma = pd.Series([1.0, 2.0], index=[30, 40]) + limits = calculate_limits_vectorized('Xbar', mean=0.0, sigma=sigma, N=np.array([4, 4])) + assert list(limits.index) == [30, 40] + + +class TestFitXbarGroup: + def test_none_when_no_subgroup_is_replicated(self): + out = pd.DataFrame({'s': [np.nan, np.nan], 'xbar': [1.0, 2.0], 'n': [1, 1], 'N': [1, 1]}) + assert _analysis()._fit_xbar_group(out) is None + + def test_one_reading_subgroups_are_charted_but_not_fitted(self): + out = pd.DataFrame( + {'s': [np.nan, 1.0, 3.0, 2.0], 'xbar': [9.0, 10.0, 11.0, 12.0], 'n': [1, 3, 3, 0], 'N': [3] * 4} + ) + group = _analysis()._fit_xbar_group(out) + assert group is not None + assert list(group.out['n']) == [1, 3, 3] # n = 0 dropped, n = 1 kept + assert list(group.fit) == [False, True, True] + assert group.s_bar == pytest.approx(2.0) + assert group.n_to_use == 'n' # sizes vary on the Xbar chart (1, 3, 3) ... + assert group.s_n_to_use == 'N' # ... but not on the S chart (3, 3) + assert group.n_max == 3 + assert group.n_avg is None + assert list(group.out['N']) == [3, 3, 3] + + def test_average_n_mode_forces_the_constant_n_path(self): + out = pd.DataFrame({'s': [1.0, 3.0, 2.0], 'xbar': [10.0, 11.0, 12.0], 'n': [2, 4, 3], 'N': [4] * 3}) + group = _analysis(n_mode='average')._fit_xbar_group(out) + assert group is not None + assert group.n_avg == pytest.approx(3.0) + assert group.n_to_use == group.s_n_to_use == 'N' + assert list(group.out['N']) == [3.0, 3.0, 3.0] + assert list(group.out['n']) == [2, 4, 3] # actual sizes survive on n + + +class TestApplySLimits: + def test_b3_b4_per_row_around_s_bar(self): + out = pd.DataFrame({'s': [0.5, 4.0, 1.2], 'n': [3, 5, 5], 'N': [5] * 3}) + result = _analysis()._apply_s_limits(out.copy(), s_bar=1.5, n_to_use='n') + assert list(result['center']) == [1.5, 1.5, 1.5] + assert list(result['lpl']) == [round(b3(3) * 1.5, 3), round(b3(5) * 1.5, 3), round(b3(5) * 1.5, 3)] + assert list(result['upl']) == [round(b4(3) * 1.5, 3), round(b4(5) * 1.5, 3), round(b4(5) * 1.5, 3)] + assert list(result['beyond_limits']) == [0, 1, 0] + + +class TestAggregateSubgroups: + def test_cached_cell_mean_replaces_the_subgroup_mean(self): + df = pd.DataFrame( + { + 'k': ['A', 'A', 'A', 'B', 'B'], + 't': [1, 1, 2, 1, 1], + 'y': [1.0, 3.0, 5.0, 2.0, 4.0], + 'r2': [np.nan, 0.5, -0.5, 0.2, -0.2], + 'Ybar_kt': [2.0, 2.0, 5.0, 3.0, 3.0], + } + ) + analysis = _analysis() + plain = analysis._aggregate_subgroups(df, ['k', 't'], 'y', 'y') + assert list(plain.columns) == ['k', 't', 's', 'xbar', 'n', 'N'] + assert list(plain['xbar']) == [2.0, 5.0, 3.0] + assert list(plain['n']) == [2, 1, 2] + assert list(plain['N']) == [2, 2, 2] + + cached = analysis._aggregate_subgroups(df, ['k', 't'], 'y', 'y', ybar_col='Ybar_kt') + assert list(cached['xbar']) == [2.0, 5.0, 3.0] + assert 'Ybar_kt' not in cached.columns and 'mean' not in cached.columns + + residual = analysis._aggregate_subgroups(df, ['k', 't'], 'y', 'r2') + # n and s count the limits column: the subgroup missing its R2 is one smaller + assert list(residual['n']) == [1, 1, 2] + assert 's_value' in residual.columns + + +class TestSubgroupLabelColumn: + def test_default_grouping_joins_into_group(self): + out = pd.DataFrame({'factor 1': [1, 2], 'time': [10, 10], 'xbar': [0.0, 0.0]}) + labelled, col = _analysis(by=None)._subgroup_label_column(out, ['factor 1', 'time']) + assert col == 'group' and list(labelled['group']) == ['1_10', '2_10'] + + def test_single_by_keeps_the_column_as_strings(self): + out = pd.DataFrame({'factor 1': [1, 2], 'xbar': [0.0, 0.0]}) + labelled, col = _analysis(by=('factor 1',))._subgroup_label_column(out, ['factor 1']) + assert col == 'factor 1' and list(labelled['factor 1']) == ['1', '2'] + + def test_several_by_columns_join_into_subgroup(self): + out = pd.DataFrame({'factor 1': [1, 2], 'time': [10, 10], 'xbar': [0.0, 0.0]}) + _, col = _analysis(by=('factor 1', 'time'))._subgroup_label_column(out, ['factor 1', 'time']) + assert col == 'subgroup' + + def test_no_grouping_has_no_label_column(self): + out = pd.DataFrame({'group': ['All'], 'xbar': [0.0]}) + _, col = _analysis(by=None)._subgroup_label_column(out, []) + assert col is None diff --git a/validation/README.txt b/validation/README.txt index 5496519..0013389 100644 --- a/validation/README.txt +++ b/validation/README.txt @@ -51,6 +51,13 @@ What lives here and why mR signals) and which do not (the X chart's "not one process" conclusion). Run from the repo root: python validation/mr_permutation_invariance.py validation/aco_per_capita_expenditure.csv +- xbar_s_sweep.py - fingerprints every Xbar/S request on PM SDS 1 and PM SDS 3 + of the T100 database (chart x companion x seven by= shapes x value x n_mode, + 504 requests): one SHA-256 per request over the chart frames, statistics and + metadata, or the exception a refused request raises. Run it on two revisions + and diff the outputs; a refactor that changes nothing prints identical files + (issue #144, PR 1). Run from the repo root: + python validation/xbar_s_sweep.py > sweep.txt Removed in the 2026-08 cleanup (recoverable from git history): sds1-6_analysis.xlsx (Excel exports generated by the package itself - not diff --git a/validation/xbar_s_sweep.py b/validation/xbar_s_sweep.py new file mode 100644 index 0000000..cf1559e --- /dev/null +++ b/validation/xbar_s_sweep.py @@ -0,0 +1,92 @@ +#!/usr/bin/env python3 +""" +Fingerprint every Xbar/S request on Tom Bishop's reference data. + +One line per request: the request, then a SHA-256 over the chart frames (every value at 17 +significant digits), the statistics and the metadata, or the exception a refused request +raises. Run it on two revisions and diff the outputs: a refactor that changes nothing prints +identical files. Written for issue #144, whose first PR extracts the per-group Xbar/S +estimator and must leave every existing chart bit-identical. + +Grid: response PM SDS 1 (ADS 1) and PM SDS 3 (ADS 3) x chart {Xbar, S} x companion {off, on} +x by {None, [], one factor, both factors, time, one factor + time, both factors + time} +x value {response, R5, R2} x n_mode {actual, average}. + +Usage, from the repo root: + python validation/xbar_s_sweep.py > sweep.txt +""" + +from __future__ import annotations + +import hashlib +import json +import sys +from pathlib import Path + +import numpy as np +import pandas as pd + +import processbehavior as pb +from processbehavior.exceptions import ProcessBehaviorError + +ROOT = Path(__file__).resolve().parent +REFERENCE = ROOT / 'PBTESTDATABASE_T100.csv' +FACTORS = ['FACTOR 1', 'FACTOR 2'] +TIME = 'PRODUCTION TIME' +RESPONSES = ('PM SDS 1', 'PM SDS 3') +BY_SHAPES = (None, [], ['FACTOR 1'], FACTORS, [TIME], ['FACTOR 1', TIME], FACTORS + [TIME]) +VALUES = (None, 'R5', 'R2') +N_MODES = ('actual', 'average') + + +def _normalise(obj): + """JSON-ready copy: numpy scalars to Python, NaN to the string 'nan', dict keys to str.""" + if isinstance(obj, dict): + return {str(k): _normalise(v) for k, v in obj.items()} + if isinstance(obj, (list, tuple)): + return [_normalise(v) for v in obj] + if isinstance(obj, np.generic): + obj = obj.item() + if isinstance(obj, float) and np.isnan(obj): + return 'nan' + return obj + + +def fingerprint(result) -> str: + digest = hashlib.sha256() + for name in sorted(result.charts): + payload = result.charts[name] + digest.update(name.encode()) + digest.update(payload['data'].to_csv(index=True, float_format='%.17g').encode()) + digest.update(json.dumps(_normalise(payload['statistics']), sort_keys=True).encode()) + digest.update(json.dumps(_normalise(payload['metadata']), sort_keys=True).encode()) + digest.update(json.dumps(_normalise(payload.get('strata')), sort_keys=True).encode()) + return digest.hexdigest() + + +def main() -> int: + df = pd.read_csv(REFERENCE, na_values=['*']) + for response in RESPONSES: + study = pb.formulate(df, response=response, factors=FACTORS, time=TIME) + for chart in ('Xbar', 'S'): + for companion in (False, True): + for by in BY_SHAPES: + for value in VALUES: + for n_mode in N_MODES: + request = ( + f'{response} | chart={chart} companion={companion} by={by} ' + f'value={value} n_mode={n_mode}' + ) + try: + result = study.execute( + chart=chart, by=by, value=value, companion=companion, n_mode=n_mode + ) + except ProcessBehaviorError as exc: + print(f'{request} -> refused {type(exc).__name__}: {exc}') + continue + print(f'{request} -> {fingerprint(result)}') + return 0 + + +if __name__ == '__main__': + sys.exit(main())