#!/usr/bin/env python3 """ Univariate analysis of bicorder readings: per-protocol and per-gradient averages. Reproduces the ad-hoc averages workflow described in README.md for any readings CSV, version-agnostic: gradient columns are discovered from the file itself (canonicalized via bicorder_common.COLUMN_RENAMES), so every run directory can be analyzed with the same command. Outputs (default /analysis/): plots/protocol_averages.png — protocol averages, ascending plots/gradient_averages.png — gradient averages with gaps between the three sets plots/averages_histogram.png — distribution of protocol averages data/protocol_averages.csv — per-protocol mean over all gradients data/gradient_averages.csv — per-gradient mean/median/coverage reports/univariate_summary.txt — printed summary With --img, also publishes the three summary PNGs to an image directory (analysis/img/ by default — the charts the README links from img/). Usage: python3 scripts/univariate_analysis.py data/synthetic_1.4.0/readings.csv python3 scripts/univariate_analysis.py data/synthetic_1.4.0/readings.csv --output data/synthetic_1.4.0/analysis python3 scripts/univariate_analysis.py data/synthetic_1.4.0/readings.csv --img # also refresh img/ """ import argparse import shutil import sys from pathlib import Path import warnings warnings.filterwarnings('ignore') import pandas as pd import numpy as np import matplotlib.pyplot as plt from bicorder_common import apply_renames, csv_version, dimension_columns MIDPOINT = 5 # all gradients run 1 (hard) .. 9 (soft) def main(): parser = argparse.ArgumentParser( description='Univariate analysis of Protocol Bicorder readings', formatter_class=argparse.RawDescriptionHelpFormatter, epilog=""" Examples: python3 scripts/univariate_analysis.py data/synthetic_1.4.0/readings.csv python3 scripts/univariate_analysis.py data/synthetic_1.2.6/readings.csv --output data/synthetic_1.2.6/analysis """, ) parser.add_argument('csv_file', help='Readings CSV (e.g. data/synthetic_1.4.0/readings.csv)') parser.add_argument('--output', '-o', default=None, help='Output directory (default: /analysis)') parser.add_argument('--img', nargs='?', const='img', default=None, help='Publish the three summary PNGs to an image directory in addition ' 'to the run analysis outputs (default with --img: img/ at the ' "analysis root — the charts the README links from img/)") args = parser.parse_args() if not Path(args.csv_file).exists(): print(f"Error: File not found: {args.csv_file}") sys.exit(1) dataset_dir = Path(args.csv_file).parent output_dir = Path(args.output) if args.output else dataset_dir / 'analysis' (output_dir / 'plots').mkdir(parents=True, exist_ok=True) (output_dir / 'data').mkdir(parents=True, exist_ok=True) (output_dir / 'reports').mkdir(parents=True, exist_ok=True) version = csv_version(args.csv_file) print("=" * 80) print("PROTOCOL BICORDER - UNIVARIATE ANALYSIS" + (f" (bicorder v{version})" if version else "")) print(f"Source: {args.csv_file}") print("=" * 80) df = pd.read_csv(args.csv_file) df = apply_renames(df) # Identify gradient columns, grouped by set in bicorder.json order all_dims = dimension_columns(df.columns.tolist()) groups = [[c for c in all_dims if c.startswith(prefix)] for prefix in ('Design_', 'Entanglement_', 'Experience_')] dimension_cols = [c for g in groups for c in g] if not dimension_cols: print("Error: no gradient columns found") sys.exit(1) values = df[dimension_cols].apply(pd.to_numeric, errors='coerce') # --- Per-protocol averages --- protocol = pd.DataFrame({ 'Descriptor': df.get('Descriptor', pd.Series(range(len(df)))), 'average': values.mean(axis=1), 'n_gradients_scored': values.notna().sum(axis=1), }).dropna(subset=['average']).sort_values('average').reset_index(drop=True) protocol.to_csv(output_dir / 'data' / 'protocol_averages.csv', index=False) print(f"\nSaved: {output_dir / 'data' / 'protocol_averages.csv'}") # --- Per-gradient averages --- gradient = pd.DataFrame({ 'gradient': dimension_cols, 'set': [c.split('_', 1)[0] for c in dimension_cols], 'mean': [values[c].mean() for c in dimension_cols], 'median': [values[c].median() for c in dimension_cols], 'coverage': [values[c].notna().mean() for c in dimension_cols], }).sort_values('mean', ascending=False) gradient.to_csv(output_dir / 'data' / 'gradient_averages.csv', index=False) print(f"Saved: {output_dir / 'data' / 'gradient_averages.csv'}") # --- Summary statistics on protocol averages --- avg = protocol['average'] mean, median, std = avg.mean(), avg.median(), avg.std() pearson_skew = 3 * (mean - median) / std if std else float('nan') moment_skew = avg.skew() # bias-corrected Fisher moment skewness lines = [] lines.append(f"Readings analyzed: {len(protocol)} protocols x {len(dimension_cols)} gradients") if version: lines.append(f"Bicorder version: {version}") lines.append("") lines.append("Protocol averages:") for label, value in [('mean', mean), ('median', median), ('std', std), ('min', avg.min()), ('max', avg.max())]: lines.append(f" {label:8s} = {value:.3f}") lines.append(f" midpoint = {MIDPOINT}") lines.append(f" deviation from midpoint = {mean - MIDPOINT:+.3f} " f"(normalized over half-range 4: {(mean - MIDPOINT) / 4:+.3f})") lines.append(f" Pearson skew (3*(mean-median)/std) = {pearson_skew:+.3f}") lines.append(f" Fisher moment skew = {moment_skew:+.3f}") lines.append("") lines.append("Gradient averages (all sets):") for _, row in gradient.iterrows(): lines.append(f" {row['gradient'][:60]:60s} mean={row['mean']:.2f} coverage={row['coverage']:.0%}") lines.append("") lines.append(f"Extremes: highest = {gradient.iloc[0]['gradient']} ({gradient.iloc[0]['mean']:.2f}); " f"lowest = {gradient.iloc[-1]['gradient']} ({gradient.iloc[-1]['mean']:.2f})") summary_text = "\n".join(lines) report_path = output_dir / 'reports' / 'univariate_summary.txt' report_path.write_text(summary_text + "\n") print(f"\nSaved: {report_path}\n") print(summary_text) # --- Plots --- # Protocol averages, ascending (cf. img/protocol_averages.png) fig, ax = plt.subplots(figsize=(12, 8)) ax.plot(range(len(protocol)), protocol['average'], linewidth=1.2) ax.axhline(MIDPOINT, color='gray', linestyle='--', linewidth=1, label=f'midpoint ({MIDPOINT})') ax.set_title('Protocol averages (ascending order)' + (f' — bicorder v{version}' if version else '')) ax.set_xlabel('Protocol (ranked by average)') ax.set_ylabel('Average gradient value') ax.set_ylim(0, 10) ax.legend() plt.tight_layout() path = output_dir / 'plots' / 'protocol_averages.png' plt.savefig(path, dpi=300, bbox_inches='tight') plt.close() print(f"\nSaved: {path}") # Histogram (cf. img/averages_histogram.png) fig, ax = plt.subplots(figsize=(10, 6)) ax.hist(avg, bins=40, color='#4C72B0', edgecolor='white') ax.axvline(MIDPOINT, color='gray', linestyle='--', linewidth=1, label=f'midpoint ({MIDPOINT})') ax.axvline(mean, color='#C44E52', linewidth=1, label=f'mean ({mean:.2f})') ax.set_title('Distribution of protocol averages') ax.set_xlabel('Average gradient value') ax.set_ylabel('Number of protocols') ax.legend() plt.tight_layout() path = output_dir / 'plots' / 'averages_histogram.png' plt.savefig(path, dpi=300, bbox_inches='tight') plt.close() print(f"Saved: {path}") # Gradient averages, sorted within each set, with visual gaps between sets # (cf. img/gradient_averages.png) palette = {'Design': '#4C72B0', 'Entanglement': '#55A868', 'Experience': '#C44E52'} ordered_cols = [c for g in groups for c in sorted(g, key=lambda cv: values[cv].mean(), reverse=True)] means = [values[c].mean() for c in ordered_cols] colors = [palette[c.split('_', 1)[0]] for c in ordered_cols] x_pos, pos = [], 0.0 for gi, group in enumerate(groups): group_sorted = sorted(group, key=lambda cv: values[cv].mean(), reverse=True) if gi > 0: pos += 1.5 # visual gap between gradient sets x_pos.extend(pos + i for i in range(len(group_sorted))) pos += len(group_sorted) fig, ax = plt.subplots(figsize=(14, 6)) ax.bar(x_pos, means, color=colors) ax.axhline(MIDPOINT, color='gray', linestyle='--', linewidth=1, label=f'midpoint ({MIDPOINT})') ax.set_title('Gradient averages (gaps separate the three gradient sets)' + (f' — bicorder v{version}' if version else '')) ax.set_xticks(x_pos) ax.set_xticklabels([c.replace('_vs_', '\nvs\n') for c in ordered_cols], fontsize=7) ax.set_ylabel('Average value') ax.set_ylim(0, 10) ax.legend() plt.tight_layout() path = output_dir / 'plots' / 'gradient_averages.png' plt.savefig(path, dpi=300, bbox_inches='tight') plt.close() print(f"Saved: {path}") # Optionally publish the summary charts for docs (the README links img/…) if args.img: img_dir = Path(args.img) if not img_dir.is_absolute(): img_dir = Path(__file__).resolve().parent.parent / img_dir img_dir.mkdir(parents=True, exist_ok=True) published = [] for name in ('protocol_averages.png', 'averages_histogram.png', 'gradient_averages.png'): shutil.copy2(output_dir / 'plots' / name, img_dir / name) published.append(str(img_dir / name)) print(f"Published charts → {', '.join(published)}") print("\nDone.") if __name__ == '__main__': main()