protocol-bicorder/analysis/scripts/univariate_analysis.py

#!/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 <dataset_dir>/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: <dataset_dir>/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()