# Experimental Plan Template (Pre-Registration Style)

## 1. Study Information

**Title:** [Descriptive title]

**Authors:** [Names]

**Date:** [Pre-registration date]

**Version:** [1.0]

## 2. Research Questions & Hypotheses

**Primary research question:** [What is the main question this study answers?]

**Hypotheses:**

- H₀ (null): [Statement of no effect]
- H₁ (alternative): [Statement of expected effect]

**Secondary questions:** [1-3 additional questions, if any]

**Exploratory questions:** [Questions without directional predictions]

## 3. Design

**Design type:** [ ] Randomized experiment [ ] Quasi-experiment [ ] Observational study [ ] Other: ______

**Design specification:** [ ] CRD [ ] RCBD [ ] Factorial [ ] Crossover [ ] Longitudinal [ ] Case-control [ ] Cohort

**Randomization unit:** [ ] Individual [ ] Cluster [ ] Block [ ] Other: ______

**Analysis unit:** [Same as randomization unit, or specify if different]

**Blinding:** [ ] None [ ] Single-blind [ ] Double-blind

## 4. Sample Size / Power

**Power analysis parameters:**

| Parameter | Value | Source |
|-----------|-------|--------|
| α (Type I error) | 0.05 | Convention |
| Power (1 − β) | 0.80 | Convention |
| Expected effect size | [value and metric] | [prior study, pilot, or MDE] |
| Minimum detectable effect | [value] | [computed] |
| Required sample size (total) | [N] | [computed] |
| Attrition adjustment | [additional N] | [expected dropout rate] |
| **Final target N** | **[N]** | |

**Power analysis output:** [Attach or reference the power analysis script output]

## 5. Participants / Subjects

**Target population:** [Who/what are we studying?]

**Inclusion criteria:**

1. [Criterion 1]
2. [Criterion 2]

**Exclusion criteria:**

1. [Criterion 1]
2. [Criterion 2]

**Recruitment:** [How will participants be identified and recruited?]

**Assignment:** [How will participants be assigned to conditions?]

## 6. Variables

### Primary Outcome

| Variable | Definition | Measurement | Type |
|----------|------------|-------------|------|
| [name] | [operational definition] | [instrument, units] | [continuous/binary/time-to-event] |

### Secondary Outcomes

| Variable | Definition | Measurement | Type |
|----------|------------|-------------|------|
| [name] | | | |

### Predictors / Treatments

| Variable | Levels | Assignment | Notes |
|----------|--------|------------|-------|
| [name] | [level1, level2, ...] | [randomized/observed] | |

### Covariates (pre-registered)

| Variable | Rationale |
|----------|-----------|
| [name] | [why this variable is included] |

## 7. Procedure

**Timeline:**

| Phase | Activity | Duration |
|-------|----------|----------|
| 1 | Recruitment | [time] |
| 2 | Pre-treatment measurement | [time] |
| 3 | Treatment administration | [time] |
| 4 | Post-treatment measurement | [time] |
| 5 | Follow-up | [time] |

**Detailed protocol:** [Step-by-step description of what happens to each participant]

## 8. Analysis Plan

### Primary Analysis

**Method:** [e.g., Independent t-test, linear regression, ANCOVA]

**Model specification:** [e.g., outcome ~ treatment + covariate1]

**Assumptions to be checked:**

| Assumption | Diagnostic | Remediation If Violated |
|------------|------------|------------------------|
| Normality | Shapiro-Wilk, Q-Q plot | Transform data or use nonparametric test |
| Equal variance | Levene's test | Welch's correction |
| Independence | Durbin-Watson | GLS / mixed model |
| [other] | [method] | [alternative] |

### Secondary Analyses

[Method for each secondary question]

### Exploratory Analyses

[Methods for exploratory questions]

### Subgroup Analyses (pre-specified)

| Subgroup | Rationale |
|----------|-----------|
| [group] | [why this subgroup might differ] |

### Missing Data Handling

**Expected missingness:** [Amount and pattern expected]

**Primary approach:** [ ] Complete case analysis [ ] Multiple imputation [ ] Maximum likelihood [ ] Last observation carried forward [ ] Other: ______

### Multiple Testing

**Correction method:** [ ] Bonferroni [ ] Holm [ ] Benjamini-Hochberg [ ] None (confirmatory study with single primary outcome) [ ] Other: ______

**Number of comparisons:** [ ]

## 9. Data Collection & Management

**Data collection tools:** [ ] Survey platform [ ] Lab equipment [ ] Administrative records [ ] Other: ______

**Data storage:** [How and where will data be stored?]

**Quality control:** [Checks for data quality during collection]

## 10. Ethics & Reporting

**Ethics approval:** [ ] Obtained (IRB #: ______) [ ] Pending [ ] Not required

**Consent:** [ ] Written [ ] Verbal [ ] Waived

**Conflicts of interest:** [None / describe]

**Data availability:** [Where will data and code be posted?]

## 11. Changes from Original Plan

[If this is an update to a pre-registration, document any changes here with date and rationale]

# Analysis Report Template

## 1. Question & Context

**Research question:** [One sentence: what are we trying to learn?]

**Motivation:** [Why does this question matter? What decision depends on the answer?]

**Data source:** [Where did the data come from? Collection method, timeframe, sample frame]

**Pre-registration:** [Link to preregistration if applicable. If not, flag any exploratory analyses.]

---

## 2. Data

**Sample size:** N = [n] ([n1] in group 1, [n2] in group 2)

**Inclusion/exclusion criteria:** [Who/what was included and excluded, and why]

**Missing data:** [Amount, pattern (MCAR/MAR/MNAR), handling method]

**Variables:**

| Variable | Type | Role | Description |
|----------|------|------|-------------|
| [name] | [continuous/binary/ordinal/etc.] | [outcome/predictor/covariate] | [description] |
| ... | | | |

---

## 3. Methods

**Analytic approach:** [e.g., two-sample t-test, linear regression with covariates, Bayesian hierarchical model]

**Justification:** [Why this method? What assumptions are we willing to make?]

**Pre-specified analyses:** [What was planned before seeing the data]

**Exploratory analyses:** [What was added after seeing the data]

**Software:** [Python 3.x with scipy/statsmodels/scikit-learn, R 4.x with package vX]

---

## 4. Assumption Checks

| Assumption | Method | Result | Status |
|------------|--------|--------|--------|
| Normality (Group 1) | Shapiro-Wilk | W = [value], p = [value] | ✓ / ✗ / N/A |
| Equal variance | Levene's test | F = [value], p = [value] | ✓ / ✗ / N/A |
| ... | | | |

**Summary:** [All assumptions met / violations detected and addressed]

---

## 5. Results

**Primary analysis:**

| Estimate | SE | 95% CI | Test Statistic | p-value | Effect Size [95% CI] |
|----------|-----|--------|----------------|---------|---------------------|
| [value] | [value] | [lower, upper] | [t/χ²/F/Z = value] | [value] | [d/η²/V/OR = value [CI]] |

**Secondary analyses:**

[Brief summary of secondary results]

**Visualization:**

[Figure: appropriate plot with clear axes, uncertainty shown, caption below]

*Figure 1: [Caption describing what the reader should see]*

---

## 6. Diagnostics & Robustness

**Sensitivity analyses:**

| Analysis | Result | Conclusion |
|----------|--------|------------|
| Main analysis | [original estimate] | — |
| Excluding outliers | [estimate] | Consistent / different |
| Alternative specification | [estimate] | Robust / sensitive |
| Different analysis method | [estimate] | Robust / sensitive |

**Residual diagnostics:** [Pattern in residuals? Influential points?]

---

## 7. Limitations

1. **[Assumption that may be violated]:** [How this could affect results]
2. **[Confounding not addressed]:** [Direction and magnitude of potential bias]
3. **[Generalizability concern]:** [Population or setting limits]
4. **[Measurement issue]:** [Reliability, validity of measures]

---

## 8. Conclusion

[One paragraph answering the original question, with appropriate uncertainty. Include:

- What we found (with effect size and precision)
- What we didn't find (null results with equivalence if applicable)
- What remains uncertain
- Practical implications

Example: "We found moderate evidence that the intervention increases response rate by 12 percentage points (95% CI [4, 20], p = 0.003, d = 0.45). This effect was robust to excluding outliers and controlling for baseline covariates. However, the result may not generalize to non-English-speaking populations, and the mechanism remains unclear."]

---

## Appendix

**Full model output:** [Link or table]

**Code:** [Link to repository]

**Data:** [Access information or note about availability]

{
  "schema_version": 1,
  "skill_name": "data-scientist",
  "evals": [
    {
      "id": "ab-test-design-power",
      "prompt": "We want to test a new onboarding flow that we believe will increase activation rate from 20% to 22%. How many users do we need in the experiment, how long should it run, and what analysis should we do at the end? I want to be rigorous and avoid a false-positive-driven launch.",
      "expected_output": "An experiment design that starts by stating the unit of randomization (user), the metric (activation rate), and the minimum detectable effect (2 percentage points), then computes required sample size per arm using a standard two-proportion power calculation (alpha 0.05, power 0.8), accounting for multiple metrics and multiple variants with a correction if applicable. It covers duration planning from expected daily traffic plus a buffer for novelty effects and seasonality, pre-registers the primary metric and stopping rule, and prescribes the analysis: check sample ratio mismatch, compute confidence interval on the effect, run sensitivity checks, and distinguish statistical significance from practical significance before launch.",
      "assertions": [
        "The response specifies the randomization unit, primary metric, baseline rate, and minimum detectable effect before computing sample size",
        "The response computes sample size with stated alpha, power, and a two-proportion formula",
        "The response plans duration from daily traffic including buffers for novelty and seasonality",
        "The response pre-registers the primary metric and stopping rule and checks for sample ratio mismatch",
        "The response distinguishes statistical from practical significance before recommending launch"
      ]
    },
    {
      "id": "causal-inference-vs-correlation",
      "prompt": "Sales data shows customers who attend our webinars churn 40% less than those who do not. My boss wants to make webinars the centerpiece of the retention strategy based on this. Is that justified, and what would it take to actually establish causality?",
      "expected_output": "A response that resists the correlational conclusion: webinar attendees differ systematically from non-attendees (they are more engaged, more likely to be on certain plans, earlier in lifecycle), so the naive comparison suffers from selection bias and confounding. It proposes the hierarchy of evidence for the question: randomized encouragement designs, natural experiments or instrument variables, difference-in-differences using a roll-out, or propensity-score/regression adjustments as weaker alternatives, and specifies what data would be needed to support each. It also states what analysis should be run now to quantify the selection bias (compare observables between groups) before any investment decision.",
      "assertions": [
        "The response flags selection bias and confounding as the core problem with the observed comparison",
        "The response explains why attendees differ systematically from non-attendees and how that undermines the causal claim",
        "The response proposes an identification strategy such as randomized encouragement, diff-in-diff, or instrumental variables",
        "The response includes a near-term analysis comparing observables between groups to quantify selection",
        "The response does not endorse the webinar strategy on the correlation alone"
      ]
    },
    {
      "id": "model-selection-task",
      "prompt": "We need to predict which accounts will churn in the next 30 days so our sales team can intervene. We have 40k accounts, 120 features with lots of missing values, class imbalance (about 5% churn), and the team has been tuning XGBoost for weeks. How should I frame model selection here, and what should drive the final choice?",
      "expected_output": "A model-selection framing that leads with the business decision context: churn prediction is a ranking task for intervention, so evaluation should use recall-at-k or precision-at-k at the intervention capacity, not raw accuracy on an imbalanced set. It recommends a baseline (logistic regression or simple heuristic) before complex models, a proper train/validation/test split that respects time ordering (no random split leaking future information), handling of missingness that is validated rather than assumed, and a cost-aware threshold choice based on the cost of a false positive versus a missed churner. The response compares the XGBoost candidate against baselines with the ranking metric and states that the choice is justified by validated lift, not tuning effort.",
      "assertions": [
        "The response reframes evaluation around ranking metrics (recall-at-k or precision-at-k) tied to intervention capacity",
        "The response mandates a time-respecting split rather than a random split",
        "The response requires a simple baseline before accepting the tuned model",
        "The response treats threshold choice as cost-aware, weighing false positives against missed churners",
        "The response rejects accuracy as the evaluation metric on an imbalanced set"
      ]
    },
    {
      "id": "bayesian-vs-frequentist",
      "prompt": "We ran an experiment and the frequentist analysis says the effect is not significant (p=0.09). A colleague says we should switch to a Bayesian analysis because it will let us conclude there is a high probability the change is positive. Is that a valid reason to switch analysis methods?",
      "expected_output": "A response that distinguishes the legitimate from the illegitimate uses of Bayesian analysis: switching after peeking because the frequentist result is not convenient is p-hacking by another name, and a Bayesian analysis with a flat prior run after the fact will not manufacture evidence. It explains that a Bayesian approach can add value when designed up front: an informative prior based on prior experiments, a decision rule on the posterior (P(effect > 0) and expected loss), and sequential monitoring that is principled. It notes that the two frameworks answer different questions and that the analysis choice must be pre-registered, and it shows how to compute the posterior probability of a positive effect and the posterior probability of a practically meaningful effect from the observed data.",
      "assertions": [
        "The response flags switching methods after seeing the p-value as post-hoc analysis rather than principled",
        "The response explains that a flat-prior Bayesian analysis run post hoc does not create evidence",
        "The response describes when Bayesian analysis is genuinely useful: informative priors, decision rules, principled sequential monitoring",
        "The response distinguishes the question each framework answers and requires pre-registration of the analysis plan",
        "The response computes or specifies computing P(effect > 0) and the posterior probability of a meaningful effect"
      ]
    },
    {
      "id": "analysis-report-uncertainty",
      "prompt": "I ran a regression analysis on customer spend and found a coefficient for the new pricing plan of +$12/month. I need to write a report for leadership. What should the report contain beyond the coefficient, and how should I communicate the uncertainty?",
      "expected_output": "An analysis report structured for decision-makers: the question and the decision it informs, the data and its limitations, the model and its key assumptions stated plainly, and the estimate with a confidence interval rather than a single point, expressed in decision-relevant language (range of plausible effects, probability of the effect being positive or economically meaningful if a Bayesian interpretation is used). The report discloses confounders and omitted-variable risk, checks robustness (alternative model specifications, sensitivity to outliers), and ends with what would change the conclusion. It avoids overprecision and states clearly what is measured versus assumed.",
      "assertions": [
        "The response structures the report around the decision the analysis informs",
        "The response communicates the estimate with a confidence interval rather than a single point",
        "The response discloses model assumptions, confounders, and omitted-variable risk",
        "The response includes robustness checks such as alternative specifications or sensitivity to outliers",
        "The response states what is measured versus assumed and what would change the conclusion"
      ]
    },
    {
      "id": "interpretability-target-and-scope",
      "prompt": "A regulator asks why a credit model declined an applicant. Design an interpretability analysis and report without confusing a debugging artifact with a customer explanation.",
      "expected_output": "Defines the audience, decision, local explanation target, model/data versions, reference background, method assumptions, sensitive fields, and a separate bounded customer explanation. It validates stability and states causal limits.",
      "assertions": ["Defines audience and explanation target", "Records model, data, and background versions", "Separates diagnostic from user-facing explanation", "States causal limits and sensitive-data controls"]
    },
    {
      "id": "interpretability-correlation-perturbation",
      "prompt": "A SHAP report says income is the main driver, but income and debt are highly correlated and synthetic perturbations create impossible applicants. How should this be interpreted?",
      "expected_output": "Flags attribution sharing among correlated features and invalid perturbations, tests alternate backgrounds or conditional approaches where justified, and reports the result as model association rather than causal effect.",
      "assertions": ["Identifies correlated-feature attribution ambiguity", "Rejects impossible perturbations as evidence", "Uses sensitivity or alternate background checks", "Avoids causal overclaiming"]
    },
    {
      "id": "interpretability-stability-disagreement",
      "prompt": "Two explanation methods disagree on a high-impact case, and randomizing labels leaves the saliency map almost unchanged. What is the decision?",
      "expected_output": "Runs method-appropriate sanity and stability checks, treats unchanged saliency after randomization as a warning, preserves method disagreement, and refuses a confident explanation until the discrepancy is investigated.",
      "assertions": ["Uses randomization sanity checks", "Treats unchanged saliency as a warning", "Preserves disagreement between methods", "Does not issue a confident explanation prematurely"]
    },
    {
      "id": "interpretability-fairness-slice",
      "prompt": "The model's overall fairness metric passes, but a small protected subgroup has much worse false-negative rates. How should the analysis be reported?",
      "expected_output": "Reports subgroup counts and uncertainty, investigates measurement and data coverage, keeps the subgroup failure visible despite the aggregate result, and avoids universal fairness claims.",
      "assertions": ["Reports subgroup counts and uncertainty", "Investigates coverage and measurement limits", "Does not average away the subgroup failure", "Avoids universal fairness claims"]
    }
  ]
}

# Bayesian Workflow Reference

## The Bayesian Workflow (Gelman et al. 2020)

Bayesian analysis is not "put a prior on everything and sample." It's an iterative process:

```
1. Model building
   │
   ├─ 2. Prior predictive check (do priors produce plausible data?)
   │  └─ If not, revise priors
   │
   ├─ 3. Inference (MCMC, VI, or exact)
   │
   ├─ 4. Posterior predictive check (does model reproduce observed data?)
   │  └─ If not, revise model
   │
   ├─ 5. Model comparison (WAIC, LOO, cross-validation)
   │
   └─ 6. Sensitivity analysis (try different priors, different likelihoods)
```

---

## 1. Prior Elicitation

### Types of Priors

| Prior Type | Description | When to Use |
|-----------|-------------|-------------|
| **Flat / Improper** | No information, ∝ 1 | Default reference, rarely recommended |
| **Weakly Informative** | Wide but finite: Normal(0, 10), Cauchy(0, 5) | Default for most parameters; shrink extreme values |
| **Regularizing** | Shrink toward zero: Normal(0, 1), Laplace(0, 1) | Many predictors, prevents overfitting |
| **Informative** | Narrowly centered on known value | Strong prior knowledge (previous studies, physical constraints) |
| **Hierarchical** | Prior on prior parameters (hyperpriors) | Multiple groups sharing information |

### Prior Recommendation Cheat Sheet

| Parameter Type | Default Prior | Why |
|---------------|--------------|-----|
| **Regression coefficient** (continuous predictor) | Normal(0, 1) or Normal(0, 2.5) | 95% of coefficients within ±2-5 units on logit/probit scale |
| **Intercept** | Normal(0, 10) | Wide enough for most scales |
| **Standard deviation** (positive) | Half-Cauchy(0, 2.5), Exponential(1), or Inverse-Gamma(3, 3) | Half-Cauchy preferred; avoids near-zero mass of inverse-gamma |
| **Correlation** (between -1 and 1) | LKJ(2) | Slightly regularizes toward zero correlation |
| **Variance parameter** (hierarchical) | Exponential(1) or Half-Cauchy(0, 2) | Robust to scale differences |
| **Proportion / probability** | Beta(1, 1) or Beta(2, 2) | Uniform vs slightly central |
| **Log-odds** (logistic) | Normal(0, 1.5) or Student-t(3, 0, 2.5) | Gelman et al. 2008: default for logistic regression |

### Prior Predictive Check Protocol

1. Sample from the prior distribution
2. Generate synthetic data from the prior-predictive distribution
3. Plot the distribution of synthetic data
4. Ask: "Could data generated this way actually occur in my problem?"
5. If the synthetic data implies impossible values (negative variances, 150-year-old humans, etc.): tighten the prior
6. If the synthetic data is too constrained (never allows realistic values): widen the prior

---

## 2. Computation

### MCMC Methods

| Method | Software | When to Use |
|--------|----------|-------------|
| **Hamiltonian Monte Carlo (HMC)** | Stan (via PyStan, CmdStanPy), PyMC (NUTS), NumPyro (NUTS) | Default for continuous parameters. Efficient exploration of posterior. |
| **Variational Inference (VI)** | PyMC (ADVI), NumPyro (SVI), Stan (VB) | Large datasets, fast approximate inference. Check fit against HMC on subset. |
| **Gibbs Sampling** | JAGS, BUGS | Legacy. Only when HMC is impractical. |
| **SMC / ABC** | PyMC (SMC), custom | Discrete parameters, likelihood-free inference |
| **Laplace Approximation** | INLA, statsmodels | Fast, accurate for latent Gaussian models (GLMM, spatial) |

### MCMC Diagnostics

| Diagnostic | What It Checks | Threshold | How to Fix |
|-----------|---------------|-----------|------------|
| **R-hat (Ř)** | Convergence across chains | < 1.01 (modern threshold, Vehtari et al. 2021) | Run longer, reparameterize, increase warmup |
| **Effective Sample Size (ESS)** | Number of independent samples | bulk-ESS > 400, tail-ESS > 400 | Run longer, thinner |
| **Trace plot** | Mixing, stationarity | No trends, good overlap between chains | Re-run with more iterations |
| **Autocorrelation** | Sampling efficiency | Lag-1 autocorrelation < 0.5 | Re-parameterize (centered vs non-centered) |
| **Divergent transitions** | HMC-specific: geometry problems | 0 divergences | Increase adapt_delta (0.80 → 0.95-0.99), reparameterize |
| **BFMI (Bayesian Fraction of Missing Information)** | Energy transition efficiency | BFMI > 0.2-0.3 | Re-parameterize non-centered |

### Reparameterization: Centered vs Non-Centered

| Type | Formula | When |
|------|---------|------|
| **Centered** | μ_normal = α + βX, σ ~ Half-Cauchy | Weak data, strong priors |
| **Non-centered** | μ_raw ~ Normal(0, 1), μ = α + βX + σ·μ_raw | Strong data, weak priors, hierarchical models |

Non-centered breaks the dependence between group-level mean and variance parameters, improving HMC geometry. Try non-centered when you see divergent transitions.

---

## 3. Model Comparison

### Information Criteria

| Criterion | What It Measures | When to Use | Formula |
|-----------|-----------------|-------------|---------|
| **WAIC** (Watanabe-Akaike) | Expected log pointwise predictive density + correction | Bayesian models, MCMC | WAIC = lppd − p_waic |
| **LOO-CV** (Leave-One-Out) | Same as WAIC but with Pareto-smoothed importance sampling | Bayesian models, MCMC (via `loo()` in ArviZ, loo package) | Computed via PSIS-LOO |
| **AIC** | −2 logL + 2k | Frequentist models, no prior | — |
| **BIC** | −2 logL + k log n | Frequentist models, nested comparison | — |
| **DIC** | Deviance + effective parameters | Bayesian models, MCMC (less stable than WAIC) | — |

### Model Comparison Decision Tree

```
Compare candidate models:
│
├─ Same likelihood, same data, nested?
│  ├─ Use information criteria (WAIC/LOO) for predictive comparison
│  └─ Use cross-validation for predictive accuracy
│
├─ Different likelihoods (e.g., Poisson vs Negative Binomial)?
│  └─ Use WAIC/LOO on same data, or cross-validation of predictive score
│
├─ Need to quantify evidence for H₀ vs H₁?
│  └─ Use Bayes Factor (BF₁₀) — but sensitive to prior choice. Report sensitivity.
│
└─ Which model is more useful given decision context?
    └─ Consider predictive performance (CV), interpretability, computational cost
```

### Bayes Factor Guidelines

| BF₁₀ | Evidence for H₁ |
|------|----------------|
| 1-3 | Weak / Anecdotal |
| 3-10 | Moderate |
| 10-30 | Strong |
| 30-100 | Very Strong |
| >100 | Extreme |

**Warning:** Bayes factors are very sensitive to prior choice. A diffuse prior can arbitrarily deflate the BF in favor of H₁. Always report sensitivity to prior width.

### Cross-Validation in Bayesian Models

- **Approximate LOO-CV** via PSIS (Pareto-smoothed importance sampling) — built into ArviZ
- **k-fold CV** — run inference on k-1 folds, evaluate on held-out fold. More robust but more expensive
- **Diagnostic:** Pareto k parameter. k < 0.5: reliable. 0.5 < k < 0.7: ok but watch. k > 0.7: unreliable (use k-fold instead)

---

## 4. Posterior Predictive Checking

### Core Idea
If the model is adequate, data generated from the posterior predictive distribution should resemble the observed data.

### Checks

| Check | What It Tests | Visualization |
|-------|--------------|--------------|
| **Mean / variance** | First and second moments | Overlay observed statistic on posterior of statistic from replicated datasets |
| **Extreme values** | Tail behavior | Probability of max(y_rep) ≥ max(y) |
| **Distribution shape** | Whole-distribution fit | Density overlay of y_rep and y |
| **Discrepancy measures** | Specific features you care about | Bayesian p-value (χ², skewness, proportion of zeros) |
| **Residual vs fitted** | Model structure violations | Plot standardized residuals from posterior mean vs fitted values |

### Bayesian p-values

A Bayesian p-value near 0.5 indicates good fit. Values near 0 or 1 indicate systematic discrepancy between model and data. Because it's the posterior probability that the test statistic is more extreme in the replicated data than in the observed, it doesn't have the same issues as frequentist p-values.

---

## 5. Hierarchical / Multilevel Bayesian Models

### Shrinkage

Hierarchical models borrow strength across groups. Group-level estimates are "shrunk" toward the population mean, with the amount of shrinkage determined by the group-level variance:

```
Groups with small N → more shrinkage (less reliable estimates pulled toward mean)
Groups with large N → less shrinkage (more data to estimate their own parameter)
```

### The 8 Schools Example (classic)

The canonical hierarchical Bayesian example: 8 schools each ran an experiment. Estimates vary widely due to small samples. Hierarchical model produces more realistic, shrunken estimates with proper uncertainty.

### When to Go Hierarchical:
- Multi-level data structure (students in classes in schools)
- You want to estimate group-level parameters while sharing information
- Small sample sizes in some groups (partial pooling)
- You need to model variation across groups explicitly

---

## 6. Software Quick Reference

| Task | PyMC | NumPyro | Stan |
|------|------|---------|------|
| Basic regression | `pm.Model()` with `pm.Normal()` | `numpyro.sample()` with `dist.Normal()` | `model { ... }` block |
| MCMC sampling | `pm.sample(1000)` | `mcmc.run(rng_key, ...)` | `./sample` via cmdstan |
| HMC NUTS | Default | `NUTS` | Default |
| Variational inference | `pm.fit(method='advi')` | `SVI` with `AutoDiagonalNormal` | `stan optimize` (VB deprecated) |
| Prior predictive | `pm.sample_prior_predictive()` | `Predictive()` with prior_samples | `generated quantities` block |
| Posterior predictive | `pm.sample_posterior_predictive()` | `Predictive()` with posterior_samples | `generated quantities` block |
| LOO/WAIC | `az.loo()`, `az.waic()` | `az.loo()`, `az.waic()` | `loo::loo()` in R |
| Visualization | ArviZ (`az.plot_trace`, etc.) | ArviZ | `bayesplot`, `shinystan` |

### Installation Notes

- **PyMC:** `pip install pymc arviz`
- **NumPyro:** `pip install numpyro arviz jax jaxlib`
- **Stan (Python):** `pip install pystan` (v3) or `cmdstanpy` (must install CmdStan separately: `install_cmdstan()`)
- **Stan (R):** `install.packages("rstan")`, `brms` for formula-based interface
- **R with Python:** use `rpy2` bridge or `reticulate` in R

# Causal Inference Framework

## Two Frameworks, Unified

Causal inference rests on two complementary traditions. A PhD-level data scientist is fluent in both.

| Tradition | Core Question | Key Concepts | Key Figures |
|-----------|--------------|--------------|-------------|
| **Structural / Graphical (SCM)** | What would happen under intervention? | DAGs, d-separation, do-calculus, structural equations | Pearl, Glymour, Jewell |
| **Potential Outcomes (Rubin Causal Model)** | What is the difference between observed and counterfactual? | Treatment effects, assignment mechanism, ignorability | Rubin, Hernán, Robins |

**Both frameworks should agree when applied correctly to the same problem.** Use them as complementary lenses.

---

## Part 1: DAGs & Graphical Approach

### DAG Construction Rules

A Directed Acyclic Graph (DAG) encodes causal assumptions. Every arrow represents a direct causal effect.

1. **Include all common causes** of any two variables in the graph (observed or unobserved)
2. **Arrows go from cause to effect** (temporal order)
3. **No cycles** (hence "acyclic")
4. **Mark unobserved variables** (dashed nodes, U variables) to track what's not measured
5. **Be parsimonious** — include only variables relevant to the causal question

### Causal DAG Vocabulary

| Concept | Definition | DAG Representation |
|---------|-----------|-------------------|
| **Confounder** | Variable that causes both treatment and outcome | C → T, C → Y (back door path T ← C → Y) |
| **Collider** | Variable caused by both treatment and outcome | T → C ← Y (conditioning on C opens path) |
| **Mediator** | Variable on causal pathway from T to Y | T → M → Y (don't adjust — blocks indirect effect) |
| **Instrumental Variable** | Variable that causes T but only affects Y through T | Z → T → Y (used when unobserved confounding) |
| **Selection Bias** | Conditioning on a collider or common effect | T → C ← U → Y (conditioning on C creates spurious association) |

### Back Door Criterion

A set of variables W satisfies the back door criterion for (T, Y) if:
1. No variable in W is a descendant of T
2. W blocks every back door path between T and Y (every path with an arrow into T)

If satisfied, the causal effect of T on Y is identifiable by adjusting for W:
```
P(Y | do(T = t)) = Σ_w P(Y | T = t, W = w) P(W = w)
```

### How to Identify Confounders in a DAG

```
Identify all paths between treatment T and outcome Y.
For each path:
│
├─ Directed path T → ... → Y: causal path (leave open)
│
├─ Back door path T ← ... → Y: potential confounding
│  ├─ Can I block it by adjusting for measured variables?
│  │  ├─ Yes → include them as covariates
│  │  └─ No → unmeasured confounding; consider IV, RDD, DID, or sensitivity
│  └─ Is there a collider on the path?
│     └─ If so, DON'T adjust for it (opens the path)
```

### Do-Calculus (Pearl)

Three rules for transforming expressions containing `do(T = t)` into ordinary conditional probabilities:

1. **Insert/delet observations**: P(Y | do(T), Z, W) = P(Y | do(T), W) if (Y ⟂ Z | T, W) in the graph where incoming arrows to T are removed
2. **Action/observation exchange**: P(Y | do(T), do(Z), W) = P(Y | do(T), Z, W) if (Y ⟂ Z | T, W) in the graph where incoming arrows to Z are removed
3. **Delet actions**: P(Y | do(T), do(Z), W) = P(Y | do(T), W) if (Y ⟂ Z | T, W) in the graph where incoming arrows to T and outgoing from Z are removed

---

## Part 2: Potential Outcomes Framework

### Core Notation

| Symbol | Meaning |
|--------|---------|
| Yᵢ(1) | Outcome for unit i if treated |
| Yᵢ(0) | Outcome for unit i if untreated (counterfactual) |
| Tᵢ | Treatment indicator (1 = treated, 0 = control) |
| Yᵢ = Tᵢ·Yᵢ(1) + (1−Tᵢ)·Yᵢ(0) | Observed outcome (fundamental problem: only one potential outcome observed) |
| ITE = Yᵢ(1) − Yᵢ(0) | Individual treatment effect (unobservable) |
| ATE = E[Y(1) − Y(0)] | Average Treatment Effect |
| ATT = E[Y(1) − Y(0) | T = 1] | Average Treatment Effect on the Treated |
| CATE = E[Y(1) − Y(0) | X = x] | Conditional Average Treatment Effect |

### Key Assumptions for Causal Identification

1. **Stable Unit Treatment Value Assumption (SUTVA):** No interference between units (one unit's treatment doesn't affect another's outcome) and only one version of treatment.
2. **Ignorability / Unconfoundedness:** Y(1), Y(0) ⟂ T | X (treatment assignment is independent of potential outcomes conditional on covariates). Also called "no unmeasured confounding."
3. **Positivity / Overlap:** 0 < P(T = 1 | X = x) < 1 for all x (every unit has non-zero probability of receiving either treatment).

**These assumptions are untestable from data alone.** They must be justified by design (randomization) or defended with subject matter knowledge and sensitivity analysis.

---

## Part 3: Identification Strategies

### Method Selection Decision Tree

```
You want to estimate the causal effect of T on Y.
│
├─ Was T randomly assigned?
│  ├─ YES → Analyze as randomized experiment
│  │  ├─ Simple: difference in means + t-test
│  │  ├─ Adjust for chance imbalance: ANCOVA
│  │  └─ Account for clustering: mixed model
│  └─ NO → Use observational causal method:
│
│     Is there a variable that determines T discontinuously?
│     ├─ YES → Regression Discontinuity Design
│     │
│     └─ NO → Is there a natural experiment / shock?
│        ├─ YES → Difference-in-Differences
│        │
│        └─ NO → Do you have a valid instrument?
│           ├─ YES → Instrumental Variables
│           │
│           └─ NO → Can you measure all confounders?
│              ├─ YES → G-formula / IPTW / Matching / Doubly Robust
│              └─ NO → Can you use sensitivity analysis?
│                 ├─ YES → E-value, Rosenbaum bounds
│                 └─ NO → Consider whether causal effect is identifiable at all
```

### 1. Randomized Experiments (Gold Standard)

**When:** Treatment assigned by investigator with randomization.
**Assumptions:** SUTVA, random assignment holds (no differential attrition, no non-compliance).
**Analysis:** ITT (intention-to-treat) primary; per-protocol and IV for non-compliance as sensitivity.
**Caveats:** External validity may be limited (treatment effect in experimental sample ≠ population).

### 2. Instrumental Variables (IV)

**When:** Z affects T, Z affects Y only through T, Z is independent of confounders of T and Y.

**Assumptions (for IV to be valid):**
1. **Relevance:** Z is correlated with T
2. **Exclusion restriction:** Z affects Y only through T (no direct path Z → Y)
3. **Independence:** Z is as good as randomly assigned (or conditional on X)
4. **Monotonicity:** Z affects T in one direction for all units (for LATE interpretation)

**Estimand:** Local Average Treatment Effect (LATE) — effect for compliers (units whose treatment status is changed by Z).

**Analysis:** Two-stage least squares (2SLS):
1. T̂ = α + βZ + γX (predict T from Z)
2. Y = α + τT̂ + γX (outcome on predicted T)

**First-stage diagnostics:** F-statistic > 10 (rule-of-thumb for strong instrument).

### 3. Regression Discontinuity Design (RDD)

**When:** Treatment assigned by a cutoff on a continuous variable (the "running variable").

**Assumptions:**
1. **Continuity:** The relationship between the running variable and potential outcomes is continuous at the cutoff
2. **No manipulation:** Units cannot precisely manipulate their running variable value around the cutoff

**Analysis:**
- Local linear regression around the cutoff
- Bandwidth selection (Imbens-Kalyanaraman, cross-validation)
- Polynomial specifications with robustness checks
- McCrary density test (check for manipulation of running variable)

**Estimand:** Local Average Treatment Effect at the cutoff (LATE).

**Key diagnostic:**
- Density plot of running variable (should be smooth at cutoff)
- Balance test on covariates at cutoff (should be smooth)
- Placebo tests at fake cutoffs (should show no effect)

### 4. Difference-in-Differences (DiD)

**When:** You have pre/post data for treated and untreated groups.

**Assumptions:**
1. **Parallel trends:** In the absence of treatment, the outcome would have evolved identically in treated and untreated groups
2. **No anticipation:** Treatment has no effect before it occurs
3. **No spillover:** Treatment doesn't affect control group outcomes

**Analysis:**
Y = β₀ + β₁·Treat + β₂·Post + β₃·(Treat × Post) + ε
where β₃ is the DiD estimate.

**Key diagnostics:**
- Plot pre-treatment trends (should be parallel)
- Placebo test: fake treatment date (should show no effect)
- Event study: plot treatment effect over time
- Sensitivity to parallel trends assumption? Use Rambachan-Roth (2023) approach

**Modern developments:**
- Two-way fixed effects (TWFE) can be biased with staggered adoption
- Use Callaway & Sant'Anna (2021), Sun & Abraham (2021), or Borusyak et al. (2024) estimators
- These handle heterogeneous treatment effects over time better

### 5. Matching & Weighting

**When:** You can measure all relevant confounders (unconfoundedness).

**Propensity Score:** e(X) = P(T = 1 | X)

| Method | Description |
|--------|-------------|
| **Propensity score matching (PSM)** | Match each treated unit to control unit(s) with nearest propensity score |
| **Nearest neighbor matching** | Match on Mahalanobis distance or other distance metric directly on covariates |
| **Coarsened Exact Matching (CEM)** | Coarsen continuous variables into strata, exact match within strata |
| **Inverse Probability of Treatment Weighting (IPTW)** | Weight by 1/e(X) for treated, 1/(1−e(X)) for control |
| **Doubly Robust** | Combine outcome regression with IPTW — consistent if either model is correct |
| **Targeted Maximum Likelihood (TMLE)** | Doubly robust with additional bias correction for efficiency |

**Post-matching diagnostics:**
- Standardized mean differences < 0.1 (balance achieved)
- Variance ratios between 0.5 and 2.0 (first moments balanced)
- Love plot (visualize balance improvement)

### 6. Synthetic Control

**When:** One or a few treated units, many potential control units.

**Idea:** Construct a weighted combination of control units that matches the treated unit's pre-treatment outcome trajectory. The post-treatment difference between treated and synthetic control is the treatment effect.

**Assumptions:**
- Pre-treatment fit is good
- No interference between treated and control units
- Control units are not affected by the treatment

**Analysis:** `Synth` (Abadie et al.) or augmented SC (Ben-Michael et al.).
**Diagnostic:** Placebo test — permute treatment across control units.

---

## Part 4: Sensitivity Analysis

**All observational causal analyses should include sensitivity analysis.** The core assumptions (ignorability, parallel trends, exclusion restriction) are untestable. Sensitivity analysis asks: **how far off must the assumption be for the conclusion to change?**

### Common Sensitivity Methods

| Method | When | What It Does |
|--------|------|-------------|
| **E-value** | General unmeasured confounding | Minimum strength of association between unmeasured confounder and both T and Y needed to explain away the observed effect |
| **Rosenbaum bounds** | Matched pairs | How large would hidden bias need to be to change significance? |
| **VanderWeele & Ding** | Binary confounders | Same as E-value but for binary exposures and outcomes |
| **Imbens' method** | General | How much of the treatment effect variance would need to be due to unconfoundedness? |
| **Cinelli & Hazlett** | OLS regression | Robustness of inference to confounding (partial R²) — how much residual variation must confounders explain? |
| **Placebo tests** | DiD, RDD | Use known fake treatment dates/cutoffs; if you find an effect, the real estimate is suspect |
| **Negative controls** | General | Use a control outcome that shouldn't be affected by treatment; if you find an effect, confounding is present |

### E-value Calculator (Mental)

```
E-value = RR_obs + sqrt(RR_obs × (RR_obs − 1))
        (for RR_obs > 1; symmetric for RR_obs < 1, use 1/RR_obs)

Interpretation: If an unmeasured confounder had an association of < E-value with both
treatment and outcome (after measured confounders), it could not explain away the observed
effect. Larger E-values = more robust findings.

Example: RR = 2.0 → E-value = 2.0 + √(2.0 × 1.0) = 2.0 + 1.41 = 3.41
Meaning: An unmeasured confounder would need >3.4× association with both T and Y to 
nullify the observed effect.
```

---

## Part 5: Causal Inference in A/B Testing & Product Analytics

| Scenario | Causal Method | Key Concern |
|----------|--------------|-------------|
| Randomized feature rollout | Difference in means | Network interference (social features) |
| Natural experiment (random error, scheduled downtime) | IV / DiD | Validity of natural experiment |
| Self-selected feature adoption | Matching / IPTW | Unmeasured confounding (motivated users) |
| Geographic rollout | DiD or Synthetic Control | Parallel trends assumption |
| Pre-post with no control | Interrupted time series | Secular trends, regression to mean |
| User chooses treatment | IV (instrument: random encouragement) | Weak instrument, LATE interpretation |

---

## Quick Reference: Causal Language

| Say This | Don't Say This |
|----------|---------------|
| "The estimate suggests a causal effect of X on Y under the assumption that..." | "X causes Y" |
| "We find that X is associated with a __ change in Y..." | "X impacts Y" |
| "In this randomized experiment, the treatment caused..." | Any causal claim from observational data without caveats |
| "The result is robust to unmeasured confounding of magnitude < E-value" | "We controlled for all confounders" |
| "The effect is identified under the assumption of parallel trends" | "We used DiD so it's causal" |

# Data Science Coding Workflow

**Source validated against:** Cookiecutter Data Science, MLflow documentation, DVC documentation, Kedro documentation, established DS project conventions.  
**Last reviewed:** 2026-05-23  
**When to load:** The campaign protocol has produced results and you need to structure them into a reproducible project; or the user asks "how should I set up this DS project?"

---

## Project Directory Structure

A consistent project structure makes experiments reproducible, results findable, and collaboration possible.

### Recommended Layout

```
project/
├── data/
│   ├── raw/               # Immutable original data
│   ├── processed/          # Cleaned, feature-engineered data
│   └── external/           # External reference data (lookups, metadata)
├── notebooks/              # Exploratory analysis, prototypes
│   └── 01-exploration.ipynb
├── src/                    # Reusable code
│   ├── __init__.py
│   ├── features/           # Feature engineering
│   ├── models/             # Model definitions, training logic
│   └── utils/              # Helper functions (logging, metrics)
├── models/                 # Trained model artifacts
│   └── run_001/
│       ├── model.pt
│       └── config.json
├── reports/                # Generated analysis, figures
│   └── figures/
├── config/
│   ├── config.yaml         # Experiment configuration
│   └── params.yaml         # Hyperparameters
├── experiments/
│   └── experiment_log.json  # Structured experiment record
├── requirements.txt
├── environment.yaml        # Conda env export
├── setup.py                # If src/ is a Python package
├── Makefile                # Common commands (make train, make test)
└── README.md
```

### Quick Bootstrap

```bash
# Using Cookiecutter Data Science (cookiecutter)
python -m pip install cookiecutter
cookiecutter https://github.com/drivendata/cookiecutter-data-science
```

### The Golden Rule

**Raw data is read-only.** Never modify `data/raw/`. Always create derived data in `data/processed/` with explicit scripts. This ensures reproducibility — any change to processing is captured in the script, not hidden in a manual edit.

---

## Configuration Management

### YAML Config Pattern

```yaml
# config/config.yaml
data:
  raw_path: "data/raw/dataset.csv"
  test_size: 0.2
  random_state: 42

preprocessing:
  scaling: "standard"
  handle_missing: "median"
  categorical_encoding: "onehot"

model:
  name: "random_forest"
  params:
    n_estimators: 200
    max_depth: 10
    random_state: 42

training:
  batch_size: 32
  learning_rate: 0.001
  epochs: 100
```

### Loading Config in Python

```python
import yaml
from pathlib import Path

with open("config/config.yaml") as f:
    config = yaml.safe_load(f)

# Use config throughout the code
model_class = config["model"]["name"]
model_params = config["model"]["params"]
```

### Sweep Config (for Hyperparameter Search)

```yaml
# config/sweep.yaml
parameters:
  learning_rate: [0.0001, 0.001, 0.01]
  batch_size: [16, 32, 64]
  n_layers: [2, 4, 6]
```

### OmegaConf / Hydra (Advanced)

For complex experiment configurations with hierarchical overrides:

```python
# pip install omegaconf
from omegaconf import OmegaConf

config = OmegaConf.create("""
model:
  name: resnet50
  pretrained: true
data:
  path: ./data
  augment: true
""")

# Override from command line or code
config.model.name = "efficientnet"
```

---

## Experiment Logging

### Why Log Experiments

Without logging, you lose the mapping between code, data, hyperparameters, and results. A year later, "run_004" means nothing. Logging solves:
- **What** hyperparameters produced this result?
- **Where** is the trained model artifact?
- **When** was it trained (data version, code version)?
- **How** does this compare to previous runs?

### Minimal Logging (JSON File)

```python
import json
from datetime import datetime
from pathlib import Path

def log_experiment(
    experiment_dir: str,
    model_name: str,
    params: dict,
    metrics: dict,
    model_path: str = None,
) -> dict:
    """Log a single experiment to a JSON file."""
    log_path = Path(experiment_dir) / "experiment_log.json"
    log_path.parent.mkdir(parents=True, exist_ok=True)

    entry = {
        "timestamp": datetime.now().isoformat(),
        "model_name": model_name,
        "params": params,
        "metrics": metrics,
        "model_path": model_path,
    }

    # Append to log
    if log_path.exists():
        with open(log_path) as f:
            log = json.load(f)
    else:
        log = []
    log.append(entry)

    with open(log_path, "w") as f:
        json.dump(log, f, indent=2)

    return entry
```

### MLflow Tracking

```python
# pip install mlflow
import mlflow

mlflow.set_experiment("customer-churn")

with mlflow.start_run(run_name="random_forest_v2"):
    # Log parameters
    mlflow.log_param("n_estimators", 200)
    mlflow.log_param("max_depth", 10)

    # Log metrics
    mlflow.log_metric("f1", 0.87)
    mlflow.log_metric("accuracy", 0.91)

    # Log model
    mlflow.sklearn.log_model(pipeline, "model")

    # Log artifacts (figures, configs)
    mlflow.log_artifact("config/config.yaml")
    mlflow.log_artifact("reports/confusion_matrix.png")

    # Log tags for searchability
    mlflow.set_tag("dataset_version", "v2.1")
    mlflow.set_tag("status", "candidate")
```

### TensorBoard (for Deep Learning)

```python
from torch.utils.tensorboard import SummaryWriter

writer = SummaryWriter(log_dir="runs/experiment_1")

# Log per-epoch metrics
for epoch in range(num_epochs):
    train_loss = train_one_epoch(model, dataloader)
    val_loss, val_acc = evaluate(model, val_loader)

    writer.add_scalar("Loss/train", train_loss, epoch)
    writer.add_scalar("Loss/val", val_loss, epoch)
    writer.add_scalar("Accuracy/val", val_acc, epoch)

    # Log model graph (once)
    if epoch == 0:
        writer.add_graph(model, example_input)

# Launch: tensorboard --logdir runs/
```

### WandB (Weights & Biases)

```python
# pip install wandb
import wandb

wandb.init(project="customer-churn", config={
    "learning_rate": 0.001,
    "batch_size": 32,
    "epochs": 100,
})

# Log metrics
for epoch in range(config["epochs"]):
    loss = train_step()
    wandb.log({"loss": loss, "epoch": epoch})

# Log model
wandb.save("model.pt")
```

**When to use what:**

| Tool | Best For | Hosting |
|---|---|---|
| JSON file | Single user, no infrastructure | Local |
| MLflow | Teams, experiment comparison | Self-hosted or Databricks |
| TensorBoard | Deep learning training curves | Local |
| WandB | Collaborative DL experiments | Cloud (SaaS) |

---

## Result Serialization

| Data Type | Format | Library | Notes |
|---|---|---|---|
| Tabular data | Parquet | `pandas.DataFrame.to_parquet()` | Fast, compressed, columnar. **Best choice for most data.** |
| Metrics / hyperparams | JSON | `json.dump()` | Human-readable, universally parseable |
| Model (sklearn) | `.pkl` / `.joblib` | `joblib.dump()` | Load with `joblib.load()` |
| Model (PyTorch) | `.pt` / `.pth` | `torch.save()` | Use state_dict format |
| Model (export) | `.onnx` | `torch.onnx.export()` | Framework-neutral, deployable anywhere |
| Figures | `.png` / `.pdf` | `matplotlib.savefig()` | 300 DPI minimum for publication |
| Intermediate data | Feather | `pandas.DataFrame.to_feather()` | Fast read/write, no compression |

```python
# Parquet — best for tabular data
df.to_parquet("data/processed/features.parquet")
df = pd.read_parquet("data/processed/features.parquet")

# JSON — best for metrics
with open("reports/metrics.json", "w") as f:
    json.dump(metrics, f, indent=2)

# Joblib — best for sklearn models
import joblib
joblib.dump(pipeline, "models/pipeline_v2.pkl")
```

---

## Reproducibility

### Seed Management

```python
import random
import numpy as np
import torch

def set_all_seeds(seed: int = 42):
    """Set seeds for all random number generators used in ML."""
    random.seed(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    if torch.cuda.is_available():
        torch.cuda.manual_seed_all(seed)
        torch.backends.cudnn.deterministic = True
        torch.backends.cudnn.benchmark = False
```

### Environment Pinning

```bash
# pip: freeze exact versions
pip freeze > requirements.txt

# conda: export full environment
conda env export > environment.yaml

# pip-compile (pip-tools): layered requirements
# requirements.in has loose deps, requirements.txt has pinned
```

### Docker for Full Reproducibility

```dockerfile
FROM python:3.12-slim
WORKDIR /app
COPY requirements.txt .
RUN pip install --no-cache-dir -r requirements.txt
COPY src/ src/
COPY config/ config/
ENTRYPOINT ["python", "src/train.py"]
```

**When Docker is overkill:** Single-script analyses, exploration, individual experiment debugging. Use pinned requirements + seed setting instead.

**When Docker is necessary:** Team projects, production deployment, sharing with non-technical stakeholders, running experiments on different hardware.

### Code Version Tracking

```python
# Embed git commit hash in experiment log
import subprocess

def get_git_commit_hash():
    """Get the current git commit hash."""
    try:
        return subprocess.run(
            ["git", "rev-parse", "--short", "HEAD"],
            capture_output=True, text=True, check=True
        ).stdout.strip()
    except (subprocess.CalledProcessError, FileNotFoundError):
        return "unknown"

# Include in experiment log:
log_entry["git_commit"] = get_git_commit_hash()
```

---

## Data Versioning

### DVC (Data Version Control)

```bash
# pip install dvc
dvc init
dvc add data/raw/dataset.csv       # Tracks dataset with .dvc file
git add data/raw/dataset.csv.dvc    # Commit pointer, not data
git commit -m "add dataset v1"

# Push to remote storage
dvc remote add myremote s3://mybucket/dvc
dvc push

# Later, pull a specific version
git checkout <commit_hash>
dvc checkout  # Restores the matching data version
```

### Without DVC: Simple Hash-Based Cache

```python
import hashlib
from pathlib import Path

def hash_file(path: Path) -> str:
    """SHA-256 hash of a file for integrity checking."""
    hasher = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(65536), b""):
            hasher.update(chunk)
    return hasher.hexdigest()

# Store hash alongside experiment results
dataset_hash = hash_file("data/raw/dataset.csv")
log_entry["data_hash"] = dataset_hash
```

---

## Unit Testing for Data Science

### Test Data Pattern

```python
def test_feature_engineering():
    """Test feature engineering with a tiny known dataset."""
    # Arrange: create 5-sample dataset with known properties
    X_test = pd.DataFrame({
        "age": [25, 30, 45, 60, 35],
        "income": [50000, 60000, 80000, 120000, 75000],
        "gender": ["M", "F", "F", "M", "F"],
    })

    # Act
    result = create_features(X_test)

    # Assert: known properties
    assert result.shape[0] == 5, "Should preserve row count"
    assert "age_scaled" in result.columns, "Should have age_scaled column"
    assert result["age_scaled"].std() > 0, "Scaled values should have variance"
```

### Model Invariance Test

```python
def test_model_output_shape():
    """Model should produce correct output shape on valid input."""
    X_sample = np.random.randn(32, 10)  # 32 samples, 10 features
    y_sample = (X_sample[:, 0] > 0).astype(int)

    model = RandomForestClassifier(n_estimators=10, random_state=42)
    model.fit(X_sample, y_sample)

    predictions = model.predict(X_sample)
    assert predictions.shape == (32,), "Should output one prediction per sample"
    assert set(predictions).issubset({0, 1}), "Should predict binary classes"
```

### Data Integrity Tests

```python
def test_no_missing_values_after_imputation():
    """Preprocessing should handle all missing values."""
    # Load a sample of processed data
    X = pd.read_parquet("data/processed/features.parquet")
    assert X.isnull().sum().sum() == 0, "No missing values should remain"

def test_target_distribution():
    """Target variable should have expected distribution."""
    df = pd.read_parquet("data/processed/train.parquet")
    class_counts = df["target"].value_counts()
    # Warn if any class has < 1% prevalence
    for cls, count in class_counts.items():
        assert count / len(df) >= 0.01, f"Class {cls} has < 1% prevalence"
```

---

## Common Pitfalls

| Pitfall | Symptom | Fix |
|---|---|---|
| Notebooks with unnumbered cells | Can't reproduce order | Number cells (01-load, 02-explore, 03-model). Convert to scripts before production. |
| Hardcoded file paths | Code breaks on different machines | Use `pathlib.Path`, config files, or `os.getenv` |
| No `random_state` | Results change each run | Set seeds at the top of every script |
| Data leakage in preprocessing | Overly optimistic results | Fit preprocessors on training data only, use `Pipeline` |
| Training on full data before evaluation | No held-out test set | Always split before any modeling |
| Git-ignored data/ directory | No one else can run the code | Use DVC or document how to obtain data |
| One giant `train.py` | Hard to debug, test, reuse | Split into `features.py`, `model.py`, `train.py`, `evaluate.py` |

---

## See Also

- `references/experimental-campaign-protocol.md` — the high-level workflow this supports
- `references/pytorch-integration.md` — training loops and model persistence
- `references/sklearn-integration.md` — pipelines and model selection
- `assets/experimental-plan-template.md` — pre-registration-style planning document
- `assets/report-template.md` — analysis report format

# Docker Experiment Isolation

**When to load this reference:** Running resource-intensive experiments (deep learning training, hyperparameter sweeps) that could destabilize the host if they go wrong. Use Docker containers to contain crashes, limit resource usage, and ensure reproducibility.

---

## Why Containerize Experiments

| Risk | Without Docker | With Docker |
|---|---|---|
| OOM kills the experiment AND your SSH session | Host runs out of swap, becomes unresponsive | Container hits memory limit, gets OOM-killed. Host is fine. |
| Experiment fills up disk | `/tmp` fills, other processes fail | Container disk limit, clean shutdown |
| Python dependency conflict | `pip install` breaks other projects | Isolated Python environment per container |
| "It works on my machine" | Hours of debugging environment differences | Same Docker image = same environment |
| GPU memory leak over many trials | VRAM fragments, subsequent trials crash | Containers freed after each trial |

---

## Quick Start: Single Experiment

### Dockerfile

```dockerfile
FROM pytorch/pytorch:2.5.0-cuda12.4-cudnn9-devel

# Install Python dependencies
RUN pip install --no-cache-dir \
    scikit-learn \
    pandas \
    numpy \
    scipy \
    optuna \
    tensorboard \
    mlflow

# Copy experiment code
COPY src/ /workspace/src/
COPY config/ /workspace/config/
WORKDIR /workspace

# Run training by default
ENTRYPOINT ["python", "src/train.py"]
```

### Build and Run with Resource Limits

```bash
# Build
docker build -t experiment-runner -f Dockerfile .

# Run with resource constraints
docker run --rm \
    --gpus '"device=0"' \          # Use GPU 0 only
    --memory=16g \                  # Hard memory limit
    --memory-swap=16g \             # Disable swap (prevents swapping to disk)
    --cpus=4 \                      # Limit to 4 CPU cores
    --shm-size=8g \                 # Increased /dev/shm (needed for DataLoader with many workers)
    -v /mnt/data:/workspace/data \  # Mount data directory
    -v $(pwd)/output:/workspace/output \
    -e CUDA_VISIBLE_DEVICES=0 \
    -e WANDB_API_KEY=$WANDB_API_KEY \
    experiment-runner \
    --config config/experiment.yaml
```

---

## Resource Limit Reference

| Flag | Recommended Setting | Why |
|---|---|---|
| `--memory` | 75% of host RAM | Leaves headroom for system processes |
| `--memory-swap` | Same as `--memory` | Disables swap — OOM kills the container instead of thrashing |
| `--cpus` | Host CPU count - 2 | Leaves CPUs for system, Docker daemon, monitoring |
| `--gpus` | `"device=0"` or `"all"` | Pin to specific GPU to avoid conflicts |
| `--shm-size` | `8g` or more | PyTorch DataLoader uses `/dev/shm` for shared memory. Default 64MB is too small. |
| `--pids-limit` | `1000` | Prevents fork bombs from runaway experiments |
| `--ulimit nofile` | `1024` | Prevents file descriptor exhaustion |
| `--storage-opt size` | `50GB` | Limits container disk usage |

### Full Resource-Constrained Run

```bash
docker run --rm \
    --gpus '"device=0"' \
    --memory=16g --memory-swap=16g \
    --cpus=4 \
    --shm-size=8g \
    --pids-limit=1000 \
    --storage-opt size=50GB \
    --ulimit nofile=1024:1024 \
    -v /mnt/data:/workspace/data:ro \    # Read-only data mount
    -v $(pwd)/output:/workspace/output \
    experiment-runner
```

---

## Log Collection Pattern

### Write Logs to stdout (Docker-Friendly)

```python
# In your training script — write all logs to stdout so docker logs works
import sys
import logging

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s [%(levelname)s] %(message)s",
    stream=sys.stdout,  # ← Write to stdout, not stderr or a file
    force=True,
)

# Now docker logs captures everything
logger = logging.getLogger(__name__)
logger.info(f"Starting training with config: {config}")
```

### Capturing Logs from Container

```bash
# Run container in background
CONTAINER_ID=$(docker run -d --gpus all --memory=16g experiment-runner)

# Stream logs live
docker logs -f $CONTAINER_ID

# Wait for completion and get exit code
EXIT_CODE=$(docker wait $CONTAINER_ID)
echo "Exit code: $EXIT_CODE"

# Save logs to file
docker logs $CONTAINER_ID > "experiment_log_$(date +%Y%m%d_%H%M%S).txt"
```

---

## Experiment Sweep (Multiple Containers)

```bash
for trial in {1..10}; do
    docker run -d \
        --gpus '"device=0"' \
        --memory=16g \
        --name "experiment_trial_$trial" \
        -v $(pwd)/output:/workspace/output \
        experiment-runner \
        --trial $trial
done

# Monitor
watch -n 30 "docker ps --filter name=experiment_trial --format 'table {{.Names}}\t{{.Status}}'"

# Wait for all to finish
echo "Waiting for trials to complete..."
while [ "$(docker ps -q --filter name=experiment_trial | wc -l)" -gt 0 ]; do
    sleep 30
done
echo "All trials complete."

# Collect results
for trial in {1..10}; do
    LOG=$(docker logs "experiment_trial_$trial" 2>&1 | tail -5)
    echo "Trial $trial: $LOG"
    docker rm "experiment_trial_$trial" > /dev/null 2>&1
done
```

### Sequence of GPUs

To parallelize across multiple GPUs without conflicts:

```bash
for trial in {1..8}; do
    GPU_ID=$((trial % NUM_GPUS))
    docker run -d \
        --gpus "\"device=$GPU_ID\"" \
        --name "trial_$trial" \
        experiment-runner \
        --trial $trial
done
```

---

## Cleanup Patterns

### Always Clean Up (Prevent Disk Overload)

```bash
# Remove container on exit (--rm flag handles this)
docker run --rm ...

# Manual cleanup of all stopped experiment containers
docker container prune --filter "name=experiment_trial" --force

# Clean dangling images (old builds)
docker image prune --force --filter "until=24h"

# Clean all build cache
docker builder prune --force

# Check disk usage
docker system df
```

### GPU Memory Leak Protection

```bash
# After each trial, verify GPU memory is freed
nvidia-smi --query-gpu=memory.used --format=csv,noheader

# If VRAM isn't released, restart Docker (last resort)
# sudo systemctl restart docker
```

---

## Docker Compose for Multi-Container Campaigns

```yaml
# docker-compose.yml
version: "3.9"

services:
  baseline:
    build: .
    command: --config config/baseline.yaml
    deploy:
      resources:
        limits:
          memory: 8g
          cpus: "4"
    volumes:
      - ./data:/workspace/data:ro
      - ./output:/workspace/output
    environment:
      - CUDA_VISIBLE_DEVICES=0
    shm_size: 4g

  experiment_a:
    build: .
    command: --config config/experiment_a.yaml
    deploy:
      resources:
        limits:
          memory: 16g
          cpus: "8"
    volumes:
      - ./data:/workspace/data:ro
      - ./output:/workspace/output
    environment:
      - CUDA_VISIBLE_DEVICES=1
    shm_size: 8g
    depends_on:
      - baseline
```

```bash
# Run all experiments
docker compose up --abort-on-container-exit

# Run specific experiment
docker compose run experiment_a
```

---

## When Docker Is Not Available

### Conda Environments

```bash
# Create isolated environment
conda create -n experiment_001 python=3.12
conda activate experiment_001
pip install -r requirements.txt

# Remove when done
conda remove -n experiment_001 --all
```

### Virtual Environments

```bash
python3 -m venv .venv_experiment
source .venv_experiment/bin/activate
pip install -r requirements.txt
# ... run experiment ...
deactivate
rm -rf .venv_experiment
```

---

## See Also

- `references/subagent-experiment-supervision.md` — automated failure recovery (works inside containers)
- `references/experimental-campaign-protocol.md` — the campaign workflow
- `scripts/detect-compute.py` — know your hardware before setting limits

# Experimental Campaign Protocol

**When to load this reference:** The user has asked you to find the best approach to a problem, run a research campaign, or optimize a model — not just answer a single question. The question classifier in the parent skill routes `CAMPAIGN`-type questions here.

**Core philosophy:** Start boring, end smart. Establish a heuristic baseline first, then apply bleeding-edge research, then iterate. Never skip directly to the fancy method.

---

## Protocol Overview

```
┌─────────────────────────────────────────────────────────┐
│  1. Problem Formulation       (What are we solving?)     │
├─────────────────────────────────────────────────────────┤
│  2. Baseline Heuristic        (What's the floor?)        │
├─────────────────────────────────────────────────────────┤
│  3. Bleeding-Edge Survey      (What's the SOTA?)         │
├─────────────────────────────────────────────────────────┤
│  4. Moonshot Experiments      (What works on our data?)  │
├─────────────────────────────────────────────────────────┤
│  5. Transfer Learning         (What's already learned?)  │
├─────────────────────────────────────────────────────────┤
│  6. Hyperparameter Search     (What are best params?)    │
├─────────────────────────────────────────────────────────┤
│  7. Distillation              (What can we trim?)        │
├─────────────────────────────────────────────────────────┤
│  8. Synthesis & Handoff       (What did we learn?)       │
└─────────────────────────────────────────────────────────┘
```

---

## Phase 1: Problem Formulation

**Goal:** Translate a vague request into a well-defined question with measurable success criteria.

**Entry criteria:** User has a data science ask. Level of specificity varies.

**Steps:**

1. **Identify data type** — numeric, categorical, time series, text, spatial, hierarchical, high-dimensional
2. **Identify question type** — descriptive (what happened?), predictive (what will happen?), causal (what causes X?), mechanistic (how does X work?), exploratory (what's interesting?)
3. **Identify decision context** — what will the answer be used for? A product launch? A policy decision? An academic paper?
4. **Identify constraints** — sample size, feature count, compute budget, time budget, interpretability requirements, regulatory constraints
5. **Formalize the question** — restate in a falsifiable form. "Is model B better than model A on metric X with statistical significance?"
6. **Define success criteria** — what constitutes a win? Lift over baseline? Statistical significance? Latency budget? AUC > 0.9?

**Exit criteria:** A one-paragraph problem statement with: data description, question type, decision context, constraints, and success criteria.

**Code integration:** Write the problem statement to a `campaign.md` file at the experiment root. This becomes the source of truth that all phases reference.

**Time budget:** 15-30 minutes. If you're stuck, you haven't narrowed the scope enough.

**Failure modes:**
- Problem is too broad → "What single decision would most benefit from a better answer here?"
- No clear success criteria → "What metric would convince you to change what you're doing today?"
- User says "I don't know" → start with exploratory analysis on available data and iterate

---

## Phase 2: Baseline Heuristic

**Goal:** Establish a performance floor using simple, well-understood methods before trying anything fancy.

**Entry criteria:** Problem formulated, data accessible.

**Principles:**
- The baseline is the bar that every bleeding-edge method must clear
- If a complex method can't beat logistic regression, it's not useful
- Baselines are fast to train, interpretable, and reliable

**Steps:**

1. **Select baseline model(s) based on problem type:**

| Problem Type | Baseline Models |
|---|---|
| Binary classification | Logistic Regression, Naive Bayes, Decision Tree (max_depth=3) |
| Multi-class classification | Multinomial Logistic Regression, Random Forest (n_estimators=100) |
| Regression | Linear Regression, Ridge, Decision Tree Regressor |
| Time series forecasting | Naive forecast (last value), Seasonal Naive, ARIMA(1,1,1) |
| Recommendation | Popularity baseline, User/Item mean, KNN with Jaccard |
| Clustering | K-Means (k=sqrt(n/2)), Hierarchical with Ward linkage |
| Anomaly detection | Isolation Forest, 3-sigma rule, IQR-based |
| Text classification | TF-IDF + Logistic Regression |
| Image classification | Simple CNN (2 conv layers + 2 dense), or raw pixel + RF |

2. **Implement using scikit-learn pipelines:**
```python
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

# Build pipeline with preprocessing
numeric_features = [...]  # column names
categorical_features = [...]  # column names

preprocessor = ColumnTransformer([
    ("num", StandardScaler(), numeric_features),
    ("cat", OneHotEncoder(handle_unknown="ignore"), categorical_features),
])

pipeline = Pipeline([
    ("preprocessor", preprocessor),
    ("classifier", LogisticRegression(max_iter=1000, random_state=42)),
])

# Cross-validated evaluation
scores = cross_val_score(pipeline, X, y, cv=5, scoring="f1_macro")
print(f"Baseline F1: {scores.mean():.4f} ± {scores.std():.4f}")
```

3. **Record baseline metrics** — accuracy, F1, precision/recall, RMSE, MAE, or domain-specific metric. Store in experiment log.

4. **Diagnose baseline failures** — if the baseline fails (e.g., logistic regression returns chance-level performance), investigate before proceeding:
   - Is there signal in the data at all? (Check class balance, feature-target correlations)
   - Is preprocessing correct? (Missing values? Scaling issues? Leakage?)
   - Is the metric appropriate? (Imbalanced classes? Asymmetric costs?)

**Exit criteria:** Baseline metrics recorded. A clear "this is the bar to beat" value.

**Extension (when baselines are known to be insufficient):**
- Ensemble baseline: Random Forest with 500 trees
- Tuned baseline: `GridSearchCV` over a small grid on the pipeline

**Failure modes:**
- Baseline performs at chance → check data quality, label sanity, feature engineering
- Baseline takes too long → subsample data for baseline, use faster estimator
- Multiple baselines disagree → use the best one as reference, document all

---

## Phase 3: Bleeding-Edge Survey

**Goal:** Identify the most promising recent approaches from research that could outperform the baseline.

**Entry criteria:** Baseline established. We know what "good enough" looks like.

**Steps:**

1. **Construct search queries** based on problem + data type:
   - `"<problem type>" AND "<data type>" AND SOTA OR state-of-the-art`
   - `"<dataset name>" benchmark leaderboard`
   - `"<problem type>" 2025 2026`
   - `best practices <problem type> <data type>`

2. **Search venues:**
   - **arXiv** — use the agent's arXiv search tool with filters for recent papers (last 2 years)
   - **Papers With Code** — benchmark leaderboards with code links
   - **HuggingFace Papers** — huggingface.co/papers
   - **GitHub** — search by topic for recent repos with high stars
   - **Model Zoos** — HuggingFace Models, PyTorch Hub, TensorFlow Hub

3. **Evaluate each candidate** — for each paper/repo found, assess:
   - **Code available?** (paper without code is 10x harder to reproduce)
   - **Benchmarked on similar data?** (same domain, size, characteristics)
   - **Computational cost?** (GPU hours, model size, inference latency)
   - **Reproducibility?** (clear hyperparameters, seeded runs, open dataset)
   - **Recency?** (within 2 years for cutting edge, 3-5 for well-established)

4. **Select top 2-3 candidates** — choose methods that:
   - Have public code (priority)
   - Match our compute budget
   - Show meaningful improvement over baseline on comparable benchmarks
   - Cover different approaches (don't pick 3 variants of the same method)

**Exit criteria:** Shortlist of 2-3 methods to try, each with: paper link, code link (or "no code available — reimplement from paper"), expected compute cost, expected improvement over baseline, and risk assessment.

**Code integration:**
```python
# Phases 3-4 typically use PyTorch for deep learning methods.
# See references/pytorch-integration.md for training loop templates.
```

**Failure modes:**
- No relevant papers found → expand search to adjacent domains, or use ensemble of baselines (stacking, boosting, bagging)
- Papers found but no code → assess reimplementation cost; skip if > 2 days
- Papers require 8x our compute → note as aspirational, focus on methods within budget
- Field moves fast (NLP/CV) → prioritize papers from last 6 months

---

## Phase 4: Moonshot Experiments

**Goal:** Implement the most promising approaches and compare to baseline.

**Entry criteria:** 2-3 candidate methods selected and understood. Compute available.

**Steps:**

1. **For each candidate method:**
   a. Clone/setup the codebase or implement the approach
   b. Configure for your data (input format, preprocessing, output mapping)
   c. Run with recommended hyperparameters from the paper
   d. Evaluate on the same metrics as the baseline

2. **Compare results:**
   - Create a comparison table (method | metric | vs baseline | compute cost | notes)
   - Flag methods that underperform baseline (common and informative!)
   - Flag methods that exceed baseline — these proceed to Phase 5

3. **Document failures:** For each method that didn't work, note:
   - What was tried? (hyperparameters, data splits, preprocessing variations)
   - What went wrong? (convergence issues, OOM, poor generalization, data mismatch)
   - Could it work with more tuning? More data? Different preprocessing?

**Exit criteria:** For each candidate, a clear "pass/fail/worth-more-tuning" assessment with evidence.

**Code integration:**
```python
import torch
import torch.nn as nn
from torch.utils.data import DataLoader, TensorDataset
from torch.utils.tensorboard import SummaryWriter

# Canonical PyTorch training loop (see pytorch-integration.md for details)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
model = YourModel().to(device)
optimizer = torch.optim.AdamW(model.parameters(), lr=1e-4)
scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, patience=3)
writer = SummaryWriter(log_dir="runs/experiment_1")

for epoch in range(num_epochs):
    model.train()
    for batch_x, batch_y in train_loader:
        batch_x, batch_y = batch_x.to(device), batch_y.to(device)
        optimizer.zero_grad()
        loss = criterion(model(batch_x), batch_y)
        loss.backward()
        torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
        optimizer.step()
        writer.add_scalar("train/loss", loss.item(), epoch)
```

**Failure modes:**
- All methods fail → revisit Phase 3 with broader search. Consider ensemble of baselines.
- Method works but needs 10x compute → note as future work, use with smaller data
- Code doesn't run → spend max 2 hours debugging; if not fixed, move to next candidate

---

## Phase 5: Transfer Learning

**Goal:** Leverage pre-trained models to improve performance with less data and compute.

**Entry criteria:** At least one promising method identified. Data is compatible with pre-trained models.

**Steps:**

1. **Identify pre-trained candidates:**
   - **Images:** `torchvision.models` (ResNet, EfficientNet, ViT, ConvNeXt)
   - **Text:** `transformers` (BERT, RoBERTa, DeBERTa, Llama, GPT)
   - **Audio:** `torchaudio` models, `transformers` (Wav2Vec2, Whisper)
   - **Time series:** PatchTST, TimesNet (less standardized; check papers)
   - **Tabular:** No strong TL tradition; skip or use FT-Transformer

2. **Choose adaptation strategy:**

| Strategy | When | Compute | Data Needed |
|---|---|---|---|
| **Feature extraction** | Target domain close to source | Low | Small |
| **Fine-tuning (full)** | Target domain diverges | High | Medium |
| **Fine-tuning (LoRA)** | Target diverges, limited VRAM | Medium | Medium |
| **Fine-tuning (last layers)** | Target domain similar | Low | Small |
| **Adapter modules** | Multi-task or parameter-efficient | Low | Small |

3. **Implement feature extraction:**
```python
import torchvision.models as models

# Load pre-trained backbone, freeze it
backbone = models.resnet50(weights=models.ResNet50_Weights.IMAGENET1K_V2)
for param in backbone.parameters():
    param.requires_grad = False

# Replace the head
backbone.fc = nn.Linear(backbone.fc.in_features, num_classes)
```

4. **Implement fine-tuning with LoRA** (when VRAM-constrained):
```python
from peft import LoraConfig, get_peft_model

lora_config = LoraConfig(
    r=8,
    lora_alpha=32,
    target_modules=["query", "value"],  # for transformer models
    lora_dropout=0.1,
)

model = AutoModelForSequenceClassification.from_pretrained("bert-base-uncased")
model = get_peft_model(model, lora_config)
# Only LoRA params are trainable
```

5. **Compare to non-transfer version:** Does TL actually help on your data? Sometimes it doesn't — especially if the pre-training distribution is very different from your target.

**Exit criteria:** Transfer-learned model metrics recorded alongside non-transfer results. Decision on whether TL is worth the complexity.

**Failure modes:**
- Pre-trained model too large for VRAM → use LoRA/QLoRA, smaller variant (DistilBERT vs BERT), or CPU offloading
- Pre-training domain too different → try a different foundation model, or skip TL
- No improvement from TL → the dataset may not benefit from pre-trained features; use Phase 4 results

---

## Phase 6: Hyperparameter Optimization

**Goal:** Systematically find the best hyperparameters for the most promising model(s).

**Entry criteria:** At least one model performing above baseline. Hyperparameter search space defined.

**Steps:**

1. **Define search space:**
   - Include model architecture choices (layers, hidden size, activation)
   - Training hyperparameters (learning rate, batch size, optimizer, scheduler)
   - Data preprocessing parameters (sequence length, augmentation, normalization)
   - Regularization (dropout, weight decay, label smoothing)

2. **Choose search strategy:**

| Strategy | When | Trials Needed | Notes |
|---|---|---|---|
| Grid search | <= 3 params, discrete values | Product of values | Exhaustive but expensive |
| Random search | 3-10 params, mixed types | 30-100 | Good default — 60 random trials covers most of the space |
| Bayesian (Optuna) | 5-20 params, complex interactions | 50-300 | Best for neural networks. TPESampler by default. |
| Hyperband / ASHA | Many trials, early stopping | 100-1000 | Prunes poor trials early. Great for deep learning. |

3. **Implement with Optuna:**
```python
import optuna
from optuna.pruners import MedianPruner
from optuna.samplers import TPESampler

def objective(trial):
    lr = trial.suggest_float("lr", 1e-5, 1e-2, log=True)
    dropout = trial.suggest_float("dropout", 0.1, 0.5)
    batch_size = trial.suggest_categorical("batch_size", [16, 32, 64])
    weight_decay = trial.suggest_float("weight_decay", 1e-6, 1e-3, log=True)
    n_layers = trial.suggest_int("n_layers", 1, 4)

    model = create_model(dropout=dropout, n_layers=n_layers)
    loader = create_loader(batch_size=batch_size)
    val_score = train_and_eval(model, loader, lr, weight_decay, trial=trial)
    return val_score  # e.g., validation F1

study = optuna.create_study(
    direction="maximize",
    sampler=TPESampler(seed=42),
    pruner=MedianPruner(n_startup_trials=10, n_warmup_steps=5),
)
study.optimize(objective, n_trials=100)

print(f"Best params: {study.best_params}")
print(f"Best value: {study.best_value:.4f}")
```

4. **Analyze results:** Use `optuna.visualization` to understand parameter importance.
   - Which params matter most? (Parameter importance plot)
   - Is there overfitting? (Train vs val performance across trials)
   - Should you expand the search space in a promising region?

**Exit criteria:** Best hyperparameters identified and validated. At least one configuration significantly exceeding baseline.

**Failure modes:**
- Random search finds no good region → expand search space, use wider ranges, add more trials
- Best parameters overfit → increase regularization, use early stopping, reduce model capacity
- Optuna study takes too long → reduce trials, use Hyperband pruner, subsample training data
- All trials converge to same value → the model may have plateaued; revisit architecture or data

---

## Phase 7: Distillation

**Goal:** Compress the best-performing model into a smaller, faster, cheaper version.

**Entry criteria:** A high-performing (but likely large) model exists. There's a need for smaller model (latency, cost, edge deployment).

**When to skip:** If the best model is already small (< 100M params), or there's no deployment constraint, skip distillation.

**Steps:**

1. **Knowledge distillation** (teacher → student):
```python
import torch.nn.functional as F

def distillation_loss(student_logits, teacher_logits, labels, T=4.0, alpha=0.7):
    """
    T: temperature - higher = softer probability distribution
    alpha: weight for distillation loss vs. standard loss
    """
    soft_loss = F.kl_div(
        F.log_softmax(student_logits / T, dim=-1),
        F.softmax(teacher_logits.detach() / T, dim=-1),
        reduction="batchmean"
    ) * (T * T)  # Scale factor to keep gradients in right range
    
    hard_loss = F.cross_entropy(student_logits, labels)
    return alpha * soft_loss + (1 - alpha) * hard_loss

# Training loop with distillation
teacher_model.eval()
for batch_x, batch_y in dataloader:
    with torch.no_grad():
        teacher_logits = teacher_model(batch_x)
    student_logits = student_model(batch_x)
    
    loss = distillation_loss(student_logits, teacher_logits, batch_y)
    
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
```

2. **Pruning** (remove low-magnitude weights):
```python
import torch.nn.utils.prune as prune

# Apply L1 unstructured pruning to all linear layers
for name, module in model.named_modules():
    if isinstance(module, nn.Linear):
        prune.l1_unstructured(module, name="weight", amount=0.3)  # Remove 30% of weights
        prune.remove(module, "weight")  # Make pruning permanent
```

3. **Quantization** (reduce precision):
```python
# Post-training quantization (simplest, often good enough)
import torch.quantization as quant
quantized_model = quant.quantize_dynamic(
    model, {nn.Linear, nn.LSTM, nn.GRU}, dtype=torch.qint8
)

# Or: QAT (Quantization-Aware Training)
# model.qconfig = quant.get_default_qat_qconfig('fbgemm')
# quant.prepare_qat(model, inplace=True)
# ... train further ...
# quant.convert(model, inplace=True)
```

4. **Validate distilled model:** Compare distilled vs full model on the same test set. Document the tradeoff: size, speed, and accuracy.

**Exit criteria:** Distilled model metrics, size comparison, inference speed comparison. Decision: is the tradeoff acceptable?

**Failure modes:**
- Student model can't match teacher → increase student capacity, try softer targets (higher T), more epochs
- Pruning degrades performance → lower pruning amount, use structured pruning, fine-tune after pruning
- Quantization errors > 1% → use QAT instead of PTQ, keep higher precision on sensitive layers
- Model is already optimal → skip distillation. Document that the full model is the final artifact.

---

## Phase 8: Synthesis & Handoff

**Goal:** Package everything learned into a clear, actionable summary.

**Entry criteria:** At least one model meets success criteria. All experiments documented.

**Steps:**

1. **Summarize findings:**

| Model | Metric | vs Baseline | Compute Cost | Model Size | Notes |
|---|---|---|---|---|---|
| Baseline (LogReg) | F1=0.82 | — | 30s | 2KB | Training only |
| Paper Method A | F1=0.89 | +8.5% | 4h (GPU) | 340MB | Best accuracy |
| Paper Method B | F1=0.87 | +6.1% | 1h (GPU) | 220MB | Good tradeoff |
| Distilled (Student) | F1=0.88 | +7.3% | 6h (distill) | 89MB | 4x smaller |

2. **Recommend production path:**
   - Which model to use and why
   - Deployment requirements (GPU vs CPU, RAM, latency)
   - Monitoring plan (performance drift, data drift)
   - Retraining cadence

3. **Document what didn't work:**
   - Methods tried that underperformed baseline
   - Experiments that failed and why
   - This is often as valuable as what worked

4. **Produce final artifacts:**
   - Trained model file(s) with version
   - Reproducible training script (with seeds)
   - Requirements file with pinned versions
   - Short README for the model card

**Exit criteria:** A written report covering: what was tried, what worked, what didn't, and what to do next.

**Failure modes:**
- Nobody reads the report → keep it to one page. Use the comparison table as the centerpiece.
- No clear winner → present the tradeoffs honestly. Sometimes the answer is "it depends."
- User wants more experiments → document what the next iteration should test and why.
- Results not reproducible → check seed settings, library versions, data splits. Use pinned requirements.

---

## Appendix: Quick Reference

### When to Skip Phases

| If... | Skip to... |
|---|---|
| User only needs a quick answer, not a campaign | Don't use this protocol at all. Use the parent skill's ADVICE or ANALYSIS path. |
| Baseline already beats production | Phase 4 (skip directly to moonshots) |
| Problem is well-studied with known best practice | Phase 4 (skip survey, known method) |
| Model is for edge deployment | Phase 7 (distillation is mandatory) |
| No labeled data available | Phase 1 (reformulate as unsupervised or few-shot) |

### Experiment Directory Structure

```
experiments/
├── campaign.md                   # Problem statement (Phase 1)
├── baselines/                    # Baseline models (Phase 2)
│   └── sklearn_pipeline.pkl
├── survey/                       # Research notes (Phase 3)
│   └── candidate_papers.md
├── runs/                         # Experiment logs (Phase 4-6)
│   ├── run_001_method_a/
│   ├── run_002_method_b/
│   └── optuna_study.db           # HP search results (Phase 6)
├── distillation/                 # Distilled models (Phase 7)
│   └── student_model.pt
├── final/                        # Production artifacts (Phase 8)
│   ├── model.pt
│   ├── model_card.md
│   └── requirements.txt
└── logs/                         # Experiment tracking
    └── experiment_log.json
```

### Always Set Your Seeds

```python
import random, numpy as np, torch

def set_seed(seed: int = 42):
    random.seed(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True
    torch.backends.cudnn.benchmark = False
```

### Experiment Logging Template

Record every experiment as a structured entry:

```json
{
  "experiment_id": "run_006",
  "date": "2026-05-23",
  "phase": 4,
  "method": "Method B with dropout=0.3",
  "model_class": "TransformerClassifier",
  "params": {"n_layers": 4, "hidden_dim": 256, "dropout": 0.3, "lr": 3e-4, "batch_size": 32},
  "dataset": {"train_size": 10000, "val_size": 2000, "test_size": 2000},
  "results": {"val_f1": 0.871, "test_f1": 0.865, "val_loss": 0.34, "train_time_s": 3600},
  "hardware": {"gpu": "RTX 5070 Ti", "vram_mb": 16384},
  "notes": "Performs well but training is slow. Try reducing hidden_dim."
}
```

---

## See Also

- `scripts/detect-compute.py` — know your hardware before starting Phase 2
- `references/pytorch-integration.md` — training loops, device management, mixed precision
- `references/sklearn-integration.md` — pipelines, model selection, preprocessing
- `references/data-science-coding-workflow.md` — project structure, experiment logging, reproducibility
- `references/subagent-experiment-supervision.md` — self-healing experiment pattern
- `references/docker-experiment-isolation.md` — safe containerized execution

# Experimental Design Reference

## Design Taxonomy

```
What kind of study?
│
├─ EXPERIMENTAL (randomized assignment)
│  ├─ Completely Randomized Design (CRD)
│  │  └─ Subjects assigned randomly to treatments. Simplest design.
│  │  └─ Use when: homogeneous experimental units, no blocking factor
│  │
│  ├─ Randomized Complete Block Design (RCBD)
│  │  └─ Subjects grouped into blocks, randomized within each block
│  │  └─ Use when: known source of variability can be blocked (batch, day, location)
│  │
│  ├─ Split-Plot Design
│  │  └─ Hard-to-change factor applied at whole-plot level, easy-to-change at sub-plot
│  │  └─ Use when: some factors are expensive/difficult to change (temperature, batch)
│  │  └─ Warning: two error terms, complex analysis
│  │
│  ├─ Latin Square Design
│  │  └─ Block in two directions (row and column), each treatment appears once per row/col
│  │  └─ Use when: two nuisance factors need blocking
│  │  └─ Example: 4 operators × 4 days, testing 4 treatments
│  │
│  ├─ Factorial Design
│  │  ├─ Full factorial: all combinations of all factor levels
│  │  │  └─ Use when: interactions are of interest, factors ≤ 4
│  │  └─ Fractional factorial: subset of combinations
│  │     └─ Use when: many factors, screening design, interactions assumed negligible
│  │
│  ├─ Crossover Design
│  │  └─ Each subject receives multiple treatments in sequence
│  │  └─ Use when: washout period feasible, within-subject comparison more powerful
│  │  └─ Warning: carryover effects, period effects
│  │
│  ├─ Sequential / Adaptive Design
│  │  └─ Interim analyses, sample size re-estimation, adaptive randomization
│  │  └─ Use when: ethical concerns (clinical trials), expensive experiments
│  │  └─ Warning: complex statistics, must pre-specify stopping rules
│  │
│  └─ A/B Testing (online experiments)
│     └─ Special case of CRD with 2 treatments
│     └─ Considerations: sample ratio mismatch, novelty effects, network interference
│
├─ QUASI-EXPERIMENTAL (no randomization, but design estimates causal effect)
│  ├─ Difference-in-Differences (DiD)
│  │  └─ Compare treated vs untreated groups before and after treatment
│  │  └─ Assumption: parallel trends in absence of treatment
│  │
│  ├─ Regression Discontinuity (RD)
│  │  └─ Treatment assigned by cutoff on continuous variable
│  │  └─ Assumption: smooth relationship between running variable and outcome
│  │
│  ├─ Interrupted Time Series
│  │  └─ Compare outcome trajectory before vs after intervention
│  │  └─ Assumption: no other changes coinciding with intervention
│  │
│  └─ Propensity Score Methods
│     └─ Matching, stratification, IPTW
│     └─ Assumption: no unmeasured confounding
│
└─ OBSERVATIONAL (no randomization, descriptive or predictive)
    ├─ Cross-sectional: single time point
    ├─ Cohort: follow forward in time
    ├─ Case-control: select on outcome, look back
    └─ Ecological: aggregate-level data
```

---

## Power Analysis

### What It Answers
- Given α and β, how many subjects do I need to detect a given effect?
- Given α, N, and design, what effect size can I detect?
- Given N and expected effect, what power do I have?

### Key Parameters
- **α** (Type I error rate): usually 0.05 (two-sided) or 0.025 (one-sided)
- **β** (Type II error rate): usually 0.20 (power = 0.80)
- **Effect size**: standardized measure of the expected effect (Cohen's d, f, w, OR, etc.)
- **Sample size (N)**: total number of units
- **Design features**: number of groups, number of measurements, number of covariates, ICC (clustered designs)

### Power by Design Type

| Design | Test Statistic | Software | Key Inputs |
|--------|---------------|----------|------------|
| **Two-group comparison** (continuous, unpaired) | Two-sample t-test | `tt_ind_solve_power()` power | d, α, power → n_per_group |
| **Two-group comparison** (continuous, paired) | Paired t-test | `tt_solve_power()` power | d, α, power → n_pairs |
| **Two proportions** | z-test, chi-square | `zt_ind_solve_power()` or `power.prop.test()` | p₁, p₂, α → n_per_group |
| **One-way ANOVA** (k groups) | F-test | `power.anova()` or `pwr.anova.test()` | f, α, k → n_per_group |
| **Multiple regression** (p predictors) | F-test for R² | `FTestRegPower()` or `pwr.f2.test()` | f², p, α → N |
| **Logistic regression** | Wald test | `power.logistic()` or `vary()` approaches | OR, p_base, α → N |
| **Survival (log-rank)** | Log-rank test | `power.survival.test()` or `power.zt.survival()` | hazard ratio, median, α → events |
| **Cluster RCT** (m clusters, n/cluster) | Mixed model | `power.sim.normal()` or `clusterPower` | ICC, cluster size, m → power |
| **ANOVA interaction** | F-test | Manual or simulation | f, α, design → N |

### Sample Size Heuristics (Use Power Analysis Instead When Possible)

| Scenario | Rough Rule of Thumb |
|----------|---------------------|
| Detect large effect (d = 0.8) | ~26 per group (t-test, α=0.05, 80% power) |
| Detect medium effect (d = 0.5) | ~64 per group |
| Detect small effect (d = 0.2) | ~394 per group |
| A/B test (10% relative increase from 10% base) | ~15,000 per arm |
| A/B test (10% relative increase from 50% base) | ~3,100 per arm |
| Cluster RCT (ICC = 0.05, 20 per cluster) | Multiply individual-sample N by ~2.7 |
| Interaction in factorial design | 4× the sample for main effect |

### Power Analysis Protocol

```
1. Determine primary outcome measure
2. State minimum clinically/practically important effect
3. Choose α (usually 0.05) and desired power (usually 0.80)
4. Select appropriate test and design
5. Compute required N
6. Check feasibility: do you have this N?
7. If no: reduce effect size, accept lower power, or change design (e.g., paired instead of independent)
8. If yes: account for attrition (inflate N by expected dropout rate)
9. Pre-register the analysis plan including power assumptions
```

---

## A/B Testing Framework

### Standard Protocol

1. **Define the metric.** Primary metric must be one, pre-specified, measurable, and tied to a business/investigator decision.
2. **Determine minimum detectable effect (MDE).** What's the smallest effect worth acting on?
3. **Compute sample size.** Account for multiple metrics with Bonferroni correction on α.
4. **Randomize properly.** At the unit of analysis level. Check for sample ratio mismatch (SRM).
5. **Pre-register.** Analysis plan including exclusion criteria, stopping rule, and primary analysis method.
6. **Run for pre-computed duration.** Don't peek (or use sequential testing).
7. **Analyze.** Intention-to-treat primary analysis, per-protocol sensitivity. Report effect size with CI.
8. **Check assumptions.** Balance checks, novelty effects, network interference (SUE/stable unit treatment value assumption violation).

### Common A/B Testing Mistakes

| Mistake | Why It's Wrong | Fix |
|---------|---------------|-----|
| Peeking at results | Inflation of Type I error rate | Sequential testing (always valid confidence intervals) |
| Stopping early when significant | Same as above | Pre-specify duration, or use sequential design |
| Multiple metrics without correction | Inflated false positive rate | Pre-specify primary, use Bonferroni/Holm on secondaries |
| Sample ratio mismatch (SRM) | Indicates randomization failure | Check χ² test on group assignment ratio |
| Novelty effect | Early effect decays as users adapt | Run long enough (2+ full business cycles) |
| Network interference | Treatment spills to control (social networks, marketplace) | Cluster randomization, design experiments at higher level |
| Segment hunting | Finding significance in subgroups | Pre-specify subgroups or correct for multiple comparisons |

### Minimum Detectable Effect by Sample Size (Continuous, 80% power, α=0.05)

| N per arm | MDE (Cohen's d) | MDE (proportion, base=50%) | MDE (proportion, base=10%) |
|-----------|-----------------|---------------------------|---------------------------|
| 100 | 0.40 | ±14% pp | ±12% pp |
| 500 | 0.18 | ±6.3% pp | ±5.4% pp |
| 1,000 | 0.13 | ±4.4% pp | ±3.8% pp |
| 5,000 | 0.06 | ±2.0% pp | ±1.7% pp |
| 10,000 | 0.04 | ±1.4% pp | ±1.2% pp |
| 50,000 | 0.02 | ±0.6% pp | ±0.5% pp |

---

## Blocking & Covariate Adjustment

### When to Block
- You have a pre-treatment variable known to affect the outcome
- You have a limited number of experimental units and want to reduce error variance
- You can group units into homogeneous blocks

### When to Use Covariate Adjustment (ANCOVA)
- Continuous pre-treatment variable correlated with outcome
- Increases statistical power beyond blocking alone
- Valid even in randomized experiments (does not introduce bias if pre-specified)

### When NOT to Adjust
- Post-treatment variables (they're outcomes, not covariates — introduces selection bias)
- Variables affected by treatment (collider bias)
- Multiple covariates without pre-specification (researcher degrees of freedom)

---

## Factorial Design Quick Reference

| Factors | Full Factorial Runs | ½ Fraction Runs | Resolution |
|---------|-------------------|----------------|------------|
| 2 | 4 | — | Full |
| 3 | 8 | 4 | III (½) |
| 4 | 16 | 8 | IV (½) |
| 5 | 32 | 16 | V (½) |
| 6 | 64 | 32 | VI (½) |
| 7 | 128 | 64 | VII (½) |

**Resolution guide:**
- **Resolution III:** Main effects may be confounded with two-way interactions. Screening only.
- **Resolution IV:** Main effects clear of two-way interactions; two-way interactions may be confounded with each other.
- **Resolution V:** Main effects and two-way interactions are clear. Three-way interactions may be confounded.

# Interpretability sources

Primary sources inform method choice; none establishes a universal explanation
quality threshold.

| Source | Use |
|---|---|
| Ribeiro, Singh, and Guestrin, “Why Should I Trust You?” (LIME), KDD 2016, https://doi.org/10.1145/2939672.2939778 | Local surrogate explanations and locality assumptions. |
| Lundberg and Lee, “A Unified Approach to Interpreting Model Predictions” (SHAP), NeurIPS 2017, https://proceedings.neurips.cc/paper/2017/hash/8a20a8621978632d76c43dfd28b67767-Abstract.html | Additive attribution framing; implementation and dependence assumptions must be checked. |
| Adebayo et al., “Sanity Checks for Saliency Maps,” NeurIPS 2018, https://papers.nips.cc/paper/8160-sanity-checks-for-saliency-maps | Parameter/data randomization checks for some saliency methods. |
| Molnar, Interpretable Machine Learning, https://christophm.github.io/interpretable-ml-book/ | Living reference for method assumptions and limitations; verify cited methods against original papers. |
| Mitchell et al., “Model Cards for Model Reporting,” FAT* 2019, https://doi.org/10.1145/3287560.3287596 | Documentation of intended use, performance, and subgroup limitations. |

# Interpretability workflow

Interpretability is evidence about model behavior under a method and reference
distribution. It is not automatically a causal explanation or a user-facing
justification.

1. State the decision, audience, stakes, and explanation target: global model
   behavior, local prediction, error diagnosis, fairness investigation, or a
   contrastive question. Define whether the audience needs a diagnostic view or
   an actionable explanation.
2. Choose scope and background deliberately. Record the model, data version,
   reference population, feature preprocessing, missingness, and perturbation or
   intervention semantics. For local explanations, state the neighborhood and
   baseline; for global explanations, state the population and aggregation.
3. Check correlated features, proxies, extrapolation, and distribution shift.
   A feature attribution may be shared among correlated variables or reflect a
   model shortcut. A perturbation may create impossible records and should be
   marked invalid rather than interpreted.
4. Validate the explanation: rerun with seeds, nearby backgrounds, and relevant
   perturbations; use label or feature randomization sanity checks where
   applicable; compare at least one materially different method or model. Record
   disagreement instead of averaging it into false certainty.
5. Inspect slices, especially protected groups and high-impact cases. An
   aggregate explanation or fairness result can hide a subgroup failure. Report
   sample sizes, uncertainty, missingness, and the limits of slice comparisons.
6. Write the conclusion in predictive language unless a causal design identifies
   an intervention effect. Never turn “the model relied on X” into “X caused Y.”
   Route causal claims to [the causal-inference framework](causal-inference-framework.md).

Explanation outputs belong in a versioned report with method, target, audience,
background, stability results, method disagreement, slice results, and known
failure modes. Do not expose sensitive feature values or internal reasoning just
to make an explanation look complete.

# PyTorch Integration Reference

**Source validated against:** PyTorch 2.12 documentation (pytorch.org/docs/stable)  
**Last reviewed:** 2026-05-23  
**When to load:** The campaign protocol (Phase 3-7) or any task requires implementing, training, or deploying a PyTorch model.

---

## Device Management

### Canonical Device Pattern

Always parameterize the device. Never hardcode `"cuda"`.

```python
import torch

device = torch.device(
    "cuda" if torch.cuda.is_available()
    else "mps" if torch.backends.mps.is_available()
    else "cpu"
)
# Usage: model.to(device); tensor.to(device)
```

### Checking Device Properties

```python
if torch.cuda.is_available():
    print(f"Device: {torch.cuda.get_device_name(0)}")
    print(f"VRAM: {torch.cuda.get_device_properties(0).total_memory / 1e9:.1f} GB")
    print(f"CUDA Capability: {torch.cuda.get_device_capability(0)}")
```

**MPS (Apple Silicon):** Available on macOS 12.3+. Not all operations are supported — check `torch.backends.mps.is_available()` and `torch.backends.mps.is_built()`. Some operations fall back to CPU automatically.

---

## Training Loop Patterns

### Basic Supervised Training Loop

```python
model.train()
for epoch in range(num_epochs):
    running_loss = 0.0
    for batch_x, batch_y in dataloader:
        batch_x, batch_y = batch_x.to(device), batch_y.to(device)

        optimizer.zero_grad()
        outputs = model(batch_x)
        loss = criterion(outputs, batch_y)
        loss.backward()

        # Gradient clipping (essential for stability)
        torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)

        optimizer.step()

        running_loss += loss.item()

    avg_loss = running_loss / len(dataloader)
    print(f"Epoch {epoch}: loss={avg_loss:.4f}")
```

### Gradient Accumulation (for large models on limited VRAM)

```python
accumulation_steps = 4  # Effective batch_size = physical_batch * accumulation

for i, (batch_x, batch_y) in enumerate(dataloader):
    outputs = model(batch_x)
    loss = criterion(outputs, batch_y) / accumulation_steps
    loss.backward()

    if (i + 1) % accumulation_steps == 0:
        torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
        optimizer.step()
        optimizer.zero_grad()
```

### Validation Loop

```python
model.eval()
val_loss = 0.0
correct = 0
total = 0

with torch.no_grad():
    for batch_x, batch_y in val_loader:
        batch_x, batch_y = batch_x.to(device), batch_y.to(device)
        outputs = model(batch_x)
        val_loss += criterion(outputs, batch_y).item()
        _, predicted = torch.max(outputs, 1)
        total += batch_y.size(0)
        correct += (predicted == batch_y).sum().item()

print(f"Val Loss: {val_loss/len(val_loader):.4f}, Acc: {100*correct/total:.2f}%")
```

---

## Dataset & DataLoader

### Custom Dataset Class

```python
from torch.utils.data import Dataset, DataLoader

class CustomDataset(Dataset):
    def __init__(self, features, labels, transform=None):
        self.features = torch.tensor(features, dtype=torch.float32)
        self.labels = torch.tensor(labels, dtype=torch.long)
        self.transform = transform

    def __len__(self):
        return len(self.labels)

    def __getitem__(self, idx):
        x = self.features[idx]
        y = self.labels[idx]
        if self.transform:
            x = self.transform(x)
        return x, y
```

### DataLoader Configuration

```python
dataloader = DataLoader(
    dataset,
    batch_size=32,
    shuffle=True,
    num_workers=4,          # Set to 0 on Windows if multiprocessing issues
    pin_memory=True,        # Speeds up GPU transfer (only with CUDA)
    persistent_workers=True if num_workers > 0 else False,  # PyTorch 2.0+
    collate_fn=None,        # Custom collation for variable-length data
)
```

**Note on `num_workers`:** On Linux, 4-8 workers is typical. On macOS, keep at 0-2. On Windows, 0 is safest. The optimal value depends on the data loading speed vs GPU speed.

### Collate Function for Variable-Length Data

```python
def collate_fn(batch):
    """Pad sequences in a batch to the same length."""
    inputs, labels = zip(*batch)
    # Pad to max length in this batch
    inputs_padded = torch.nn.utils.rnn.pad_sequence(inputs, batch_first=True)
    return inputs_padded, torch.tensor(labels)
```

---

## Model Saving & Loading

### Save/Load State Dict (Recommended)

```python
# Save
torch.save({
    "epoch": epoch,
    "model_state_dict": model.state_dict(),
    "optimizer_state_dict": optimizer.state_dict(),
    "loss": loss,
    "config": {"n_layers": 4, "hidden_dim": 256},  # metadata
}, "checkpoint.pt")

# Load
checkpoint = torch.load("checkpoint.pt", map_location=device)
model.load_state_dict(checkpoint["model_state_dict"])
optimizer.load_state_dict(checkpoint["optimizer_state_dict"])
start_epoch = checkpoint["epoch"] + 1
```

### Save for Inference (weights only)

```python
torch.save(model.state_dict(), "model_weights.pt")

# Load for inference
model = YourModel()
model.load_state_dict(torch.load("model_weights.pt", map_location="cpu"))
model.eval()
```

### TorchScript / torch.export (for production deployment)

```python
# torch.export (PyTorch 2.x preferred, producesExportedProgram)
exported_program = torch.export.export(model, (example_input,))
torch.export.save(exported_program, "model.pt2")

# TorchScript (legacy, still widely supported)
scripted_model = torch.jit.script(model)
scripted_model.save("model_scripted.pt")
```

---

## Mixed Precision (AMP)

Automatic Mixed Precision trains with `float16` (or `bfloat16`) where safe and `float32` where needed. Typically 1.5-2x faster with minimal accuracy loss.

```python
from torch.amp import autocast, GradScaler

scaler = GradScaler("cuda")  # "cuda" or "cpu"

for batch_x, batch_y in dataloader:
    batch_x, batch_y = batch_x.to(device), batch_y.to(device)

    optimizer.zero_grad()

    # Autocast context manager
    with autocast(device_type="cuda"):  # or "cpu"
        outputs = model(batch_x)
        loss = criterion(outputs, batch_y)

    # Scale loss, backward, unscale, step
    scaler.scale(loss).backward()
    scaler.unscale_(optimizer)
    torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)
    scaler.step(optimizer)
    scaler.update()
```

**Key points:**
- `autocast` wraps forward pass + loss computation only
- `GradScaler` prevents underflow in small gradients
- Use `bfloat16` on Ampere+ GPUs (less numerical issues than `float16`)
- AMP is most beneficial for large models (CNNs, Transformers) — small models may not see speedup

---

## torch.compile (PyTorch 2.x JIT Compilation)

`torch.compile` compiles model graphs for faster execution with minimal code changes.

```python
# Basic usage (reduces model — returns a compiled wrapper)
compiled_model = torch.compile(model)

# With options
compiled_model = torch.compile(
    model,
    mode="reduce-overhead",  # "default", "reduce-overhead", "max-autotune"
    fullgraph=True,          # Fail if graph breaks
    dynamic=True,            # Support dynamic tensor shapes
)

# Training is identical:
output = compiled_model(batch_x)  # First call compiles, subsequent calls are fast
```

**When to use:** Models with large tensor operations (CNNs, Transformers). Not beneficial for very small models or heavy data loading bottlenecks.

**When not to use:** Dynamic control flow, custom CUDA extensions, models compiled for TorchScript export.

**Mode selection:**
- `"default"` — conservative, works for most models
- `"reduce-overhead"` — good for inference, reduces Python overhead
- `"max-autotune"` — benchmarks and picks the best backend (slow first run)

---

## Loss Functions

| Problem Type | Loss Function | Import |
|---|---|---|
| Binary classification | `BCEWithLogitsLoss` | `torch.nn.BCEWithLogitsLoss` |
| Multi-class classification | `CrossEntropyLoss` | `torch.nn.CrossEntropyLoss` |
| Multi-label classification | `BCEWithLogitsLoss` | Combines sigmoid + BCELoss |
| Regression (MSE) | `MSELoss` | `torch.nn.MSELoss` |
| Regression (MAE) | `L1Loss` | `torch.nn.L1Loss` |
| Regression (Huber) | `HuberLoss` | `torch.nn.HuberLoss` (delta parameter) |
| Contrastive / Siamese | `TripletMarginLoss` / `ContrastiveLoss` | `torch.nn.TripletMarginLoss` |
| Imbalanced classes | Weighted `CrossEntropyLoss` | Pass `weight` tensor to constructor |
| Sequence (CTC) | `CTCLoss` | `torch.nn.CTCLoss` |

```python
# Weighted loss for imbalanced data
class_weights = torch.tensor([0.2, 0.8]).to(device)  # inverse frequency
criterion = torch.nn.CrossEntropyLoss(weight=class_weights)
```

---

## Optimizers & Schedulers

| Optimizer | When | Learning Rate |
|---|---|---|
| `AdamW` | Default for most models (Transformers, CNNs) | 1e-4 to 3e-4 |
| `Adam` | Legacy; prefer AdamW (proper weight decay) | 1e-4 to 3e-4 |
| `SGD` + momentum | When Adam overfits; vision models | 1e-2 to 1e-1 |
| `AdamW` (LoRA) | Fine-tuning with LoRA | 2e-4 to 5e-4 |

```python
import torch.optim as optim

# AdamW — the default
optimizer = optim.AdamW(model.parameters(), lr=1e-4, weight_decay=0.01)

# SGD with momentum (for vision fine-tuning)
optimizer = optim.SGD(model.parameters(), lr=0.01, momentum=0.9, weight_decay=1e-4)
```

### Learning Rate Schedulers

```python
from torch.optim.lr_scheduler import (
    ReduceLROnPlateau,  # Reduce when metric plateaus
    CosineAnnealingLR,  # Cosine decay
    OneCycleLR,         # Warmup + cosine ("super-convergence")
)

# Reduce on plateau (works well for most tasks)
scheduler = ReduceLROnPlateau(optimizer, mode="min", patience=5, factor=0.5)

# One-cycle (fast convergence, needs max_lr)
scheduler = OneCycleLR(
    optimizer,
    max_lr=1e-3,
    steps_per_epoch=len(train_loader),
    epochs=num_epochs,
)

# Cosine annealing (good for Transformers)
scheduler = CosineAnnealingLR(optimizer, T_max=num_epochs)

# Warmup + cosine (standard for LLM fine-tuning)
# Implement manually:
def get_cosine_schedule_with_warmup(optimizer, warmup_steps, total_steps):
    def lr_lambda(current_step):
        if current_step < warmup_steps:
            return float(current_step) / float(max(1, warmup_steps))
        progress = float(current_step - warmup_steps) / float(max(1, total_steps - warmup_steps))
        return 0.5 * (1.0 + math.cos(math.pi * progress))
    return torch.optim.lr_scheduler.LambdaLR(optimizer, lr_lambda)
```

---

## Transfer Learning

### Feature Extraction (Freeze Backbone)

```python
import torchvision.models as models

model = models.resnet50(weights=models.ResNet50_Weights.IMAGENET1K_V2)

# Freeze all layers
for param in model.parameters():
    param.requires_grad = False

# Replace classifier
num_features = model.fc.in_features
model.fc = torch.nn.Linear(num_features, num_classes)

# Only the new layer's params are trainable
optimizer = optim.AdamW(model.fc.parameters(), lr=1e-3)
```

### Full Fine-Tuning

```python
# Unfreeze all layers
for param in model.parameters():
    param.requires_grad = True

# Use lower learning rate
optimizer = optim.AdamW(model.parameters(), lr=2e-5)

# Often helps to freeze early layers, fine-tune later layers:
for name, param in model.named_parameters():
    if "layer1" in name or "layer2" in name:
        param.requires_grad = False
```

### LoRA for Transformer Models

LoRA (Low-Rank Adaptation) trains small rank-decomposition matrices while keeping the base model frozen. Requires the `peft` library.

```python
from peft import LoraConfig, get_peft_model

lora_config = LoraConfig(
    r=8,                    # Rank — higher = more expressiveness, more params
    lora_alpha=32,          # Scaling factor
    target_modules=["q_proj", "v_proj"],  # Which modules to apply LoRA to
    lora_dropout=0.1,
    bias="none",            # Don't train bias terms
    task_type="SEQ_CLS",    # Task type (SEQ_CLS, CAUSAL_LM, TOKEN_CLS, etc.)
)

model = get_peft_model(base_model, lora_config)

# Train with slightly higher LR
optimizer = optim.AdamW(model.parameters(), lr=3e-4)

# LoRA adds very few params (< 1% of base model)
model.print_trainable_parameters()  # e.g., "trainable params: 294,912 || all params: 110,080,512 || %: 0.2679"
```

---

## Knowledge Distillation

```python
import torch.nn.functional as F

def distillation_loss(student_logits, teacher_logits, labels, temperature=4.0, alpha=0.7):
    """
    Combined distillation + supervised loss.

    Args:
        temperature: Higher = softer probability distribution (more information from teacher)
        alpha: Weight for distillation loss (vs standard cross-entropy)
    """
    # Soft target loss (distillation)
    soft_student = F.log_softmax(student_logits / temperature, dim=-1)
    soft_teacher = F.softmax(teacher_logits.detach() / temperature, dim=-1)
    distill_loss = F.kl_div(soft_student, soft_teacher, reduction="batchmean")
    distill_loss *= temperature ** 2  # Scale to keep gradients in right range

    # Hard target loss
    hard_loss = F.cross_entropy(student_logits, labels)

    return alpha * distill_loss + (1 - alpha) * hard_loss


# Training loop with distillation
teacher_model.eval()
for batch_x, batch_y in dataloader:
    with torch.no_grad():
        teacher_logits = teacher_model(batch_x)

    student_logits = student_model(batch_x)
    loss = distillation_loss(student_logits, teacher_logits, batch_y)

    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
```

**Temperature tuning:** Start with T=4.0. Higher values produce softer targets (more small-class information). Lower values (T=1.0) collapse to standard cross-entropy.

**Alpha tuning:** α=0.7 (weight on distillation) is a common starting point. Increase α when the teacher is much better than the student.

---

## Model Pruning

```python
import torch.nn.utils.prune as prune

# Apply pruning to specific layers
prune.l1_unstructured(module=model.fc, name="weight", amount=0.3)  # Remove 30% of weights

# Make pruning permanent (removes the pruning mask)
prune.remove(module=model.fc, name="weight")

# Structured pruning (removes entire neurons/channels)
prune.ln_structured(module=model.conv1, name="weight", amount=0.2, n=2, dim=0)

# Global pruning (prune all layers together by importance)
parameters_to_prune = [
    (model.layer1, "weight"),
    (model.layer2, "weight"),
    (model.fc, "weight"),
]
prune.global_unstructured(
    parameters_to_prune,
    pruning_method=prune.L1Unstructured,
    amount=0.2,  # Remove 20% of weights globally
)
```

**After pruning:** Fine-tune the pruned model. Pruning then fine-tuning almost always recovers accuracy. Pruning without fine-tuning degrades performance significantly.

---

## Distributed Data Parallel (DDP)

For multi-GPU training. DDP wraps the model and handles gradient synchronization.

```python
import torch.distributed as dist
import torch.multiprocessing as mp
from torch.nn.parallel import DistributedDataParallel as DDP

def setup(rank, world_size):
    """Initialize the distributed process group."""
    dist.init_process_group(
        backend="nccl",              # "nccl" for NVIDIA, "gloo" for CPU
        init_method="env://",         # Use env vars MASTER_ADDR and MASTER_PORT
        rank=rank,
        world_size=world_size,
    )

def cleanup():
    dist.destroy_process_group()

def train(rank, world_size):
    setup(rank, world_size)

    # Model must be on the correct device BEFORE wrapping
    model = YourModel().to(rank)
    ddp_model = DDP(model, device_ids=[rank])

    # DataLoader must use DistributedSampler
    sampler = torch.utils.data.distributed.DistributedSampler(
        dataset, num_replicas=world_size, rank=rank
    )
    dataloader = DataLoader(dataset, batch_size=32, sampler=sampler)

    # Training loop (same structure, use ddp_model)
    for epoch in range(num_epochs):
        sampler.set_epoch(epoch)  # Shuffle each epoch
        for batch_x, batch_y in dataloader:
            batch_x, batch_y = batch_x.to(rank), batch_y.to(rank)
            outputs = ddp_model(batch_x)
            loss = criterion(outputs, batch_y)
            loss.backward()
            optimizer.step()

    cleanup()

# Launch
if __name__ == "__main__":
    world_size = torch.cuda.device_count()
    mp.spawn(train, args=(world_size,), nprocs=world_size)
```

**When DDP is worth it:** Models that take > 1 hour to train on a single GPU. For short experiments, single-GPU + AMP is often faster due to communication overhead.

---

## Debugging

### Common Failure Patterns

| Symptom | Likely Cause | Fix |
|---|---|---|
| `loss = nan` | Exploding gradients, bad learning rate | Lower LR, add gradient clipping, check for NaN in input data |
| `loss = nan` after AMP | Gradient underflow | Increase `GradScaler` init_scale, or use `bfloat16` |
| Loss doesn't decrease | Wrong LR, wrong loss function | Check LR range, verify loss function matches task |
| CUDA OOM | Batch size too large | Reduce batch size, enable gradient checkpointing, use AMP |
| `Expected all tensors to be on...` | Device mismatch | Always `.to(device)` tensors before model forward |
| `CUDA error: device-side assert` | Wrong label class (out of range) | Check label values are in `[0, num_classes)` |
| Model doesn't overfit 1 batch | Bug in model architecture | Try overfitting on a single batch (batch of 2-4 samples for 100 steps) |

### Overfit on One Batch (Diagnostic)

```python
# If model can't overfit a single batch, something is fundamentally wrong
single_batch = next(iter(dataloader))

for step in range(100):
    outputs = model(single_batch[0].to(device))
    loss = criterion(outputs, single_batch[1].to(device))
    loss.backward()
    optimizer.step()

print(f"Loss after 100 steps on 1 batch: {loss.item():.6f}")
# If loss is not near 0, model has a bug or LR is way off.
```

### Gradient Checking

```python
# Check gradient norms during training
total_norm = 0.0
for p in model.parameters():
    if p.grad is not None:
        param_norm = p.grad.data.norm(2)
        total_norm += param_norm.item() ** 2
total_norm = total_norm ** 0.5

if total_norm > 10.0:
    print(f"WARNING: Large gradient norm: {total_norm:.4f}")
```

---

## Reproducibility

```python
import random
import numpy as np
import torch

def set_seed(seed: int = 42):
    random.seed(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True
    torch.backends.cudnn.benchmark = False  # Trades off speed for determinism

set_seed(42)
```

**Note on `cudnn.deterministic`:** Setting this to `True` may reduce performance (5-15%). For experiments where exact reproducibility isn't critical, leave `benchmark=True` and accept minor stochastic variation.

---

## See Also

- `references/experimental-campaign-protocol.md` — where this reference fits in the campaign workflow
- `references/sklearn-integration.md` — for non-deep-learning methods
- `scripts/detect-compute.py` — check your hardware before choosing model size
- pytorch.org/docs/stable — official API documentation

# Regression Modeling Reference

## Model Selection Hierarchy

```
What's your outcome type?
│
├─ CONTINUOUS (unbounded, approximately normal)
│  └─ Linear regression → if assumptions violated, try:
│     ├─ Transform Y (log, Box-Cox, Yeo-Johnson)
│     ├─ Robust regression (M-estimation, Huber, quantile)
│     ├─ Generalized least squares (correlated errors)
│     └─ Nonparametric regression (GAM, kernel, Gaussian process)
│
├─ CONTINUOUS (bounded [0,1] or proportions)
│  └─ Beta regression (logit link) or fractional logit
│
├─ CONTINUOUS (positive, right-skewed)
│  └─ Gamma GLM (log link) or log-normal model
│
├─ BINARY (0/1)
│  ├─ Logistic regression (logit link — default)
│  ├─ Probit regression (normal CDF link — latent normal interpretation)
│  ├─ Complementary log-log (cloglog — asymmetric link, rare events)
│  └─ Robust Poisson (binary with common outcome, >10%)
│
├─ COUNT (non-negative integer)
│  ├─ Poisson regression (mean = variance)
│  │  └─ If overdispersed: Negative binomial regression
│  ├─ Zero-inflated model (excess structural zeros)
│  └─ Hurdle model (zero vs positive, then truncated count)
│
├─ ORDINAL (ordered categories)
│  └─ Proportional odds model (default)
│     └─ If proportional odds violated: partial proportional odds, adjacent category, or continuation ratio
│
├─ MULTINOMIAL (unordered categories)
│  └─ Multinomial logistic regression
│
├─ TIME-TO-EVENT
│  └─ Cox proportional hazards → if PH violated: stratified Cox, time-varying covariates, AFT
│
├─ TIME SERIES
│  └─ ARIMA / GARCH / state space / dynamic regression
│
├─ CLUSTERED / HIERARCHICAL
│  ├─ Linear/GLM + random effects (mixed models)
│  └─ GEE (population-averaged effects, robust to correlation misspecification)
│
└─ LONGITUDINAL
    ├─ Mixed effects models (subject-specific random effects)
    ├─ GEE (population-averaged)
    └─ Transition models (Markov-type, previous outcome as predictor)
```

---

## Linear Regression

### Assumptions & Diagnostics (in order of importance)

| Assumption | Violation Consequence | Diagnostic | Mitigation |
|-----------|----------------------|------------|------------|
| **Linearity** | Biased estimates, wrong functional form | Residuals vs fitted plot (check for patterns) | Add polynomials, splines, interaction terms, or GAM |
| **Independence of errors** | Inflated Type I error, wrong SEs | Durbin-Watson statistic (for serial correlation), plot residuals by order | GLS, cluster-robust SEs, mixed model |
| **Homoscedasticity** | Wrong SEs, inefficient estimates | Scale-location plot, Breusch-Pagan test | Heteroscedasticity-consistent SEs (HC1-HC3), weighted least squares |
| **Normality of residuals** | Invalid inference in small samples (large n → robust via CLT) | Q-Q plot, Shapiro-Wilk (n<50), Kolmogorov-Smirnov | Bootstrap inference, robust regression |
| **No influential points** | Estimates driven by few observations | Cook's distance, DFBETAS, DFFITS, leverage | Robust regression, drop/trim with documentation |
| **Multicollinearity** | Inflated SEs, unstable estimates | VIF > 5-10, condition index > 30 | Regularization (ridge/LASSO), remove/combine correlated predictors, PCA |

### Interpretation Guide

| Coefficient Type | Interpretation | Example |
|-----------------|---------------|---------|
| **Continuous (linear-linear)** | "A 1-unit increase in X is associated with a β-unit change in Y, holding other variables constant" | β = 2.3: "Each additional year of education is associated with a $2,300 increase in income" |
| **Binary (0/1)** | "The predicted Y is β units higher for the exposed group vs reference" | β = 5.1: "Women earn $5,100 more than men, controlling for other factors" |
| **Interaction (continuous × continuous)** | "The effect of X₁ on Y changes by β₃ for each 1-unit increase in X₂" | Simple slopes: plot at ±1 SD of moderator |
| **Interaction (binary × continuous)** | "The slope of X differs by β₃ between groups" | Plot separate regression lines for each group |
| **Log-transformed Y** | "A 1-unit change in X is associated with a 100·β % change in Y" | β = 0.03: "3% increase in Y per unit X" (approximate) |
| **Log-transformed X** | "A 1% increase in X is associated with β/100 unit change in Y" | β = 0.5: "1% more X → 0.005 unit more Y" |
| **Log-Log model** | "A 1% increase in X is associated with a β% change in Y" (elasticity) | β = 0.8: "1% more X → 0.8% more Y" |

### Common Pitfalls

| Pitfall | Why | Fix |
|---------|-----|-----|
| **Stepwise selection** | Inflated R², invalid inference, doesn't replicate | Use LASSO, domain knowledge, or AIC-based comparison of candidate models |
| **Interpreting coefficients when interactions are present** | Main effects are conditional (at zero of the moderator) | Center variables, plot marginal effects |
| **Ignoring nonlinearity** | Linear assumption hides U-shaped or threshold effects | Splines, GAMs, piecewise regression |
| **HARKing** (Hypothesizing After Results Known) | Inflated Type I error | Pre-register, split-sample (explore in half, confirm in half) |
| **p-value rounding** | "p = 0.051" is not trending | Report exact p-values and interpret continuously |
| **Not reporting uncertainty** | Overconfidence in point estimates | Always report CI/CrI with coefficients |

---

## Generalized Linear Models (GLMs)

### Canonical GLM Family Links

| Family | Default Link | Variance Function | Uses |
|--------|-------------|-------------------|------|
| Gaussian | Identity | σ² (constant) | Continuous outcomes |
| Binomial | Logit | μ(1−μ) | Binary, proportion |
| Poisson | Log | μ | Count data |
| Gamma | Inverse | μ² | Positive continuous, right-skewed |
| Inverse Gaussian | μ⁻² | μ³ | Positive continuous, very skewed |
| Negative Binomial | Log | μ + μ²/θ | Overdispersed counts |

### GLM Diagnostics (Beyond Linear)

- **Deviance residuals vs fitted** — check for pattern
- **DHARMa residuals** — simulated residuals for any GLM (recommended)
- **Overdispersion test** — for Poisson: residual deviance / df > 1.5 indicates overdispersion
- **Zero inflation** — compare observed vs predicted zeros for count models
- **Influence** — delta-betas and hat values (available via `statsmodels.graphics`)

---

## Mixed Effects / Hierarchical Models

### When to Use
- Repeated measures on same subjects (longitudinal)
- Data clustered in groups (students in schools, patients in hospitals)
- Crossed random effects (items and subjects in psycholinguistics)

### Standard Model Equation

```
Level 1: Yᵢⱼ = β₀ⱼ + β₁ⱼXᵢⱼ + εᵢⱼ
Level 2: β₀ⱼ = γ₀₀ + γ₀₁Wⱼ + u₀ⱼ
         β₁ⱼ = γ₁₀ + γ₁₁Wⱼ + u₁ⱼ
```

### Random Effect Structures

| Structure | Interpretation | N of Additional Parameters |
|-----------|---------------|---------------------------|
| Random intercept | Groups differ in baseline | 1 variance per grouping |
| Random slope + intercept | Groups differ in both baseline and covariate effect | 3 parameters: 2 variances, 1 covariance |
| Unstructured covariance | Full covariance matrix for repeated measures | k(k+1)/2 for k time points |

### Key Diagnostics
- **ICC** (intraclass correlation): proportion of variance due to between-group differences. ICC > 0.05 suggests multilevel modeling is beneficial.
- **Random effects Q-Q plot**: check normality of random effects
- **Centering**: group-mean centering for Level 1 predictors separates within- from between-group effects
- **Singular fit**: random effect variance estimated at zero → simplify random structure (Bates et al. 2015: keep maximal but drop zero-variance terms)

---

## Nonparametric & Semi-Parametric Regression

| Method | Use Case | Output |
|--------|----------|--------|
| **Smoothing splines** | Smooth nonlinear relationship, penalized | Fit with CV-chosen λ |
| **GAM (Generalized Additive Model)** | Multiple smoothed predictors, any GLM family | Partial dependence plots |
| **LOESS / LOWESS** | Local polynomial smoothing, one predictor | Smooth curve, no equation |
| **Kernel regression** | Nadaraya-Watson estimator | Smoothed conditional mean |
| **Gaussian Process** | Bayesian nonparametric regression | Full posterior over functions |
| **Regression Trees** | Decision tree for regression | Tree structure, interpretable |
| **Random Forest** | Ensemble of trees | Variable importance, partial dependence |
| **Gradient Boosting** | Sequential trees (XGBoost, LightGBM, CatBoost) | Usually best predictive accuracy |

### When to Choose Nonparametric Over Parametric
- The functional form is unknown and not theoretically specified
- The sample is large enough to estimate flexible relationships (>200 observations per smooth term)
- Prediction accuracy is more important than interpretability
- You've checked that parametric assumptions are violated and transformations don't fix it

# Scikit-Learn Integration Reference

**Source validated against:** scikit-learn 1.8.0 (scikit-learn.org/stable)  
**Last reviewed:** 2026-05-23  
**When to load:** The campaign protocol (Phase 2, 4, 6), baseline modeling, preprocessing, or any task involving sklearn estimators.

---

## Pipeline Composition

Pipelines chain preprocessing and modeling into a single estimator. This enables proper cross-validation (no data leakage from preprocessing) and simplifies deployment.

### Basic Pipeline

```python
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression

pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("classifier", LogisticRegression(max_iter=1000, random_state=42)),
])

# Use like a regular estimator
pipeline.fit(X_train, y_train)
y_pred = pipeline.predict(X_test)
```

### Shortcut: `make_pipeline`

```python
from sklearn.pipeline import make_pipeline

pipeline = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000))
# Step names are auto-generated: "standardscaler", "logisticregression"
```

### Accessing Step Attributes

```python
# After fitting
pipeline.fit(X_train, y_train)

# Access the trained scaler
scaler = pipeline.named_steps["scaler"]
print(f"Mean: {scaler.mean_}")

# Access coefficients from the classifier
coefs = pipeline.named_steps["classifier"].coef_
```

---

## ColumnTransformer (Heterogeneous Data)

When your data has both numeric and categorical columns, use `ColumnTransformer` to apply different preprocessing to different columns.

```python
from sklearn.compose import ColumnTransformer, make_column_selector
from sklearn.preprocessing import StandardScaler, OneHotEncoder

numeric_features = ["age", "income", "score"]
categorical_features = ["gender", "region", "education"]

preprocessor = ColumnTransformer([
    ("num", StandardScaler(), numeric_features),
    ("cat", OneHotEncoder(handle_unknown="ignore"), categorical_features),
])

# Or use column type selectors
preprocessor = ColumnTransformer([
    ("num", StandardScaler(), make_column_selector(dtype_include="number")),
    ("cat", OneHotEncoder(handle_unknown="ignore"),
     make_column_selector(dtype_include="object")),
])
```

### Full Pipeline with ColumnTransformer

```python
from sklearn.ensemble import RandomForestClassifier

pipeline = Pipeline([
    ("preprocessor", preprocessor),
    ("classifier", RandomForestClassifier(n_estimators=200, random_state=42)),
])

# Grid search over both preprocessing and model params
param_grid = {
    "preprocessor__num__with_mean": [True, False],
    "classifier__n_estimators": [100, 200, 500],
    "classifier__max_depth": [10, 20, None],
}

grid = GridSearchCV(pipeline, param_grid, cv=5, scoring="f1_macro")
grid.fit(X_train, y_train)
```

**Memory-Efficient ColumnTransformer:** Set `remainder="passthrough"` to keep columns not specified, or `remainder="drop"` (default) to drop them.

---

## Preprocessing

### Scaling & Normalization

| Scaler | Description | When |
|---|---|---|
| `StandardScaler` | Z-score: (x - μ) / σ | Default for most models. Assumes roughly Gaussian data. |
| `MinMaxScaler` | Scale to [0, 1] | When bounded ranges matter (neural nets, distance-based). |
| `RobustScaler` | Uses median and IQR | When data has outliers. More robust than StandardScaler. |
| `MaxAbsScaler` | Scale to [-1, 1] | For sparse data (preserves sparsity). |
| `Normalizer` | Unit norm per sample | Text classification, cosine similarity. |

```python
from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler

# StandardScaler is the default
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_train)

# Always fit on training, transform both train and test
X_test_scaled = scaler.transform(X_test)
```

### Encoding Categorical Features

```python
from sklearn.preprocessing import OneHotEncoder, OrdinalEncoder

# One-hot (nominal categories — no ordering)
encoder = OneHotEncoder(handle_unknown="ignore", sparse_output=False)
X_encoded = encoder.fit_transform(X_categorical)

# Ordinal (ordered categories)
encoder = OrdinalEncoder(categories=[["low", "medium", "high"]])
X_encoded = encoder.fit_transform(X_ordinal)
```

### Handling Missing Values

```python
from sklearn.impute import SimpleImputer, KNNImputer, IterativeImputer

# Simple imputation (fast)
imputer = SimpleImputer(strategy="median")  # "mean", "median", "most_frequent", "constant"

# KNN imputation (better for local patterns, slower)
imputer = KNNImputer(n_neighbors=5)

# Iterative imputation (MICE-style, best but slow)
imputer = IterativeImputer(max_iter=10, random_state=42)  # Experimental — requires explicit import

# In a pipeline:
pipeline = Pipeline([
    ("imputer", SimpleImputer(strategy="median")),
    ("scaler", StandardScaler()),
    ("classifier", LogisticRegression()),
])
```

---

## Model Selection

### Cross-Validation Strategies

| Splitter | Use Case |
|---|---|
| `KFold(n_splits=5, shuffle=True)` | Default for most tasks |
| `StratifiedKFold(n_splits=5)` | Classification — preserves class proportions |
| `GroupKFold(n_splits=5)` | When samples belong to groups (e.g., same patient) |
| `TimeSeriesSplit(n_splits=5)` | Temporal data — train on past, test on future |
| `RepeatedStratifiedKFold(n_repeats=3)` | More robust estimate, higher variance |
| `LeaveOneOut()` | Very small datasets (< 100 samples) |

```python
from sklearn.model_selection import (
    KFold, StratifiedKFold, GroupKFold, TimeSeriesSplit, cross_val_score, cross_validate
)

# Quick cross-validation score
scores = cross_val_score(pipeline, X, y, cv=StratifiedKFold(5), scoring="f1_macro")
print(f"F1: {scores.mean():.4f} ± {scores.std():.4f}")

# Detailed cross-validation
cv_results = cross_validate(
    pipeline, X, y,
    cv=StratifiedKFold(5),
    scoring=["f1_macro", "accuracy", "roc_auc"],
    return_estimator=True,  # Return fitted models for inspection
    return_train_score=True,  # Detect overfitting
)
```

### Grid Search

```python
from sklearn.model_selection import GridSearchCV, RandomizedSearchCV, HalvingGridSearchCV

# Grid search (exhaustive)
grid = GridSearchCV(
    pipeline,
    param_grid={
        "classifier__C": [0.01, 0.1, 1.0, 10.0],
        "classifier__penalty": ["l2"],
    },
    cv=5,
    scoring="f1_macro",
    n_jobs=-1,  # Use all CPU cores
    verbose=1,
)
grid.fit(X_train, y_train)

print(f"Best params: {grid.best_params_}")
print(f"Best score: {grid.best_score_:.4f}")

# Random search (better for high-dimensional spaces)
random_search = RandomizedSearchCV(
    pipeline,
    param_distributions={
        "classifier__C": [0.01, 0.1, 1.0, 10.0, 100.0],
        "classifier__max_iter": [500, 1000, 2000],
    },
    n_iter=20,  # Number of random combinations to try
    cv=5,
    scoring="f1_macro",
    n_jobs=-1,
    random_state=42,
)

# Halving search (successive halving — tries many candidates, prunes poor ones fast)
halving_search = HalvingGridSearchCV(
    pipeline,
    param_grid={"classifier__C": [0.01, 0.1, 1.0, 10.0]},
    factor=3,  # Reduce candidates by factor 3 each iteration
    cv=5,
    scoring="f1_macro",
    n_jobs=-1,
    verbose=1,
)
```

### Nested Cross-Validation (Unbiased Performance Estimate)

```python
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import GridSearchCV
from sklearn.tree import DecisionTreeClassifier

# Inner CV: model selection
inner_cv = StratifiedKFold(3, shuffle=True, random_state=42)
grid = GridSearchCV(DecisionTreeClassifier(),
                    {"max_depth": [3, 5, 10, None]},
                    cv=inner_cv)

# Outer CV: performance estimation
outer_cv = StratifiedKFold(5, shuffle=True, random_state=42)
nested_scores = cross_val_score(grid, X, y, cv=outer_cv, scoring="f1_macro")
# This gives an unbiased estimate of the tuned model's performance
print(f"Unbiased F1: {nested_scores.mean():.4f} ± {nested_scores.std():.4f}")
```

---

## Ensemble Methods

```python
from sklearn.ensemble import (
    RandomForestClassifier,
    GradientBoostingClassifier,
    StackingClassifier,
    VotingClassifier,
    AdaBoostClassifier,
    BaggingClassifier,
)

# Stacking (meta-model combines base models)
stack = StackingClassifier(
    estimators=[
        ("rf", RandomForestClassifier(n_estimators=100, random_state=42)),
        ("gb", GradientBoostingClassifier(n_estimators=100, random_state=42)),
        ("svc", LinearSVC(random_state=42)),
    ],
    final_estimator=LogisticRegression(),
    cv=5,
)

# Voting (simple majority or weighted average)
vote = VotingClassifier(
    estimators=[
        ("lr", LogisticRegression()),
        ("rf", RandomForestClassifier(n_estimators=100)),
        ("gnb", GaussianNB()),
    ],
    voting="soft",  # "hard" for majority vote, "soft" for probability average
)
```

### XGBoost / LightGBM Integration

```python
# sklearn-compatible API
import xgboost as xgb
import lightgbm as lgb

xgb_model = xgb.XGBClassifier(
    n_estimators=200,
    max_depth=6,
    learning_rate=0.1,
    eval_metric="logloss",
    use_label_encoder=False,
    random_state=42,
)

lgb_model = lgb.LGBMClassifier(
    n_estimators=200,
    num_leaves=31,
    learning_rate=0.1,
    random_state=42,
    verbose=-1,
)

# Both work in sklearn pipelines and GridSearchCV
pipeline = Pipeline([
    ("preprocessor", preprocessor),
    ("classifier", xgb_model),
])
```

---

## Custom Estimators

### Custom Transformer

```python
from sklearn.base import BaseEstimator, TransformerMixin

class LogTransformer(BaseEstimator, TransformerMixin):
    """Apply log(1 + x) to specified columns."""

    def __init__(self, columns=None):
        self.columns = columns  # None = all columns

    def fit(self, X, y=None):
        # LogTransform doesn't need fitting, but fit must return self
        return self

    def transform(self, X):
        X = X.copy()
        cols = self.columns if self.columns is not None else X.columns
        X[cols] = X[cols].applymap(lambda x: np.log1p(x))  # log1p = log(1+x)
        return X
```

### Custom Estimator

```python
from sklearn.base import BaseEstimator, ClassifierMixin

class SimpleThresholdClassifier(BaseEstimator, ClassifierMixin):
    """Classify based on a learned threshold on one feature."""

    def __init__(self, threshold=0.5):
        self.threshold = threshold

    def fit(self, X, y):
        # Learn optimal threshold
        # Implementation here
        self.is_fitted_ = True
        return self

    def predict(self, X):
        check_is_fitted(self)
        return (X[:, 0] > self.threshold).astype(int)

    def predict_proba(self, X):
        # Not implemented — raises error if called
        raise NotImplementedError("This estimator doesn't support probabilities")
```

### FunctionTransformer (Quick Custom Transform)

```python
import numpy as np
from sklearn.preprocessing import FunctionTransformer

# No class needed for simple transforms
log_transform = FunctionTransformer(func=np.log1p, validate=True)

# In a pipeline:
pipeline = Pipeline([
    ("log", log_transform),
    ("scaler", StandardScaler()),
])
```

---

## Persistence

```python
import joblib

# Save
joblib.dump(pipeline, "model.pkl")

# Load
loaded_pipeline = joblib.load("model.pkl")
predictions = loaded_pipeline.predict(X_new)
```

**⚠️ Security:** `joblib.load` can execute arbitrary code on deserialization. Only load models from trusted sources. Use `pickle` with the same caveat.

**Model portability:** sklearn models versioned with the sklearn version that created them. Cross-version compatibility is not guaranteed. Always save the sklearn version alongside the model.

---

## Imbalanced Data

### Built-in sklearn Support

```python
from sklearn.linear_model import LogisticRegression
from sklearn.utils.class_weight import compute_class_weight

# Option 1: Use class_weight parameter
model = LogisticRegression(class_weight="balanced", max_iter=1000)

# Option 2: Manual class weights
weights = compute_class_weight("balanced", classes=np.unique(y), y=y)
class_weight_dict = dict(zip(np.unique(y), weights))
model = LogisticRegression(class_weight=class_weight_dict, max_iter=1000)
```

### imbalanced-learn Library

```python
from imblearn.over_sampling import SMOTE, ADASYN, RandomOverSampler
from imblearn.under_sampling import RandomUnderSampler, NearMiss
from imblearn.pipeline import Pipeline as ImbPipeline  # Note: different import!

# SMOTE in pipeline (SMOTE + classifier)
pipeline = ImbPipeline([
    ("sampler", SMOTE(random_state=42)),
    ("classifier", RandomForestClassifier(n_estimators=200, random_state=42)),
])
```

---

## Calibration

```python
from sklearn.calibration import CalibratedClassifierCV

# Most sklearn classifiers output uncalibrated probabilities
# Calibrate after training for reliable probability estimates

# Method 1: Platt scaling (sigmoid) — default, good for SVMs, boosting
calibrated = CalibratedClassifierCV(model, method="sigmoid", cv=5)
calibrated.fit(X_train, y_train)
probabilities = calibrated.predict_proba(X_test)

# Method 2: Isotonic regression — non-parametric, needs more data
calibrated = CalibratedClassifierCV(model, method="isotonic", cv=5)
```

---

## Dimensionality Reduction

### PCA

```python
from sklearn.decomposition import PCA

pca = PCA(n_components=0.95)  # Keep 95% of variance
X_pca = pca.fit_transform(X_scaled)

print(f"Components: {pca.n_components_}")  # How many components retained
print(f"Explained variance: {pca.explained_variance_ratio_}")
```

**PCA assumptions:** Data should be scaled first (use `StandardScaler`). PCA assumes linear relationships. PCA is exploratory / descriptive, not inferential — it cannot confirm a hypothesis.

### t-SNE / UMAP (Visualization Only)

```python
from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, perplexity=30, random_state=42)
X_tsne = tsne.fit_transform(X_scaled)
```

**⚠️ t-SNE is for visualization only.** The embedding is stochastic and non-parametric. Different runs produce different results. Do not use t-SNE embeddings as input to other models.

---

## Feature Selection

```python
from sklearn.feature_selection import (
    SelectKBest,
    SelectFromModel,
    RFE,
    mutual_info_classif,
    chi2,
)

# Filter method (fast, univariate)
selector = SelectKBest(mutual_info_classif, k=20)
X_selected = selector.fit_transform(X, y)

# Wrapper method (RFE — Recursive Feature Elimination)
selector = RFE(estimator=RandomForestClassifier(), n_features_to_select=20)
X_selected = selector.fit_transform(X, y)

# Embedded method (from model coefficients)
selector = SelectFromModel(
    LogisticRegression(C=1.0, max_iter=1000, penalty="l1", solver="libao"),
    max_features=20,
    threshold="median",
)

# In a pipeline
pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("feature_selection", SelectKBest(mutual_info_classif, k=20)),
    ("classifier", RandomForestClassifier(n_estimators=200)),
])
```

---

## Common Pitfalls

| Pitfall | Symptom | Fix |
|---|---|---|
| Data leakage from preprocessing | Overly optimistic CV scores | Always use `Pipeline` for preprocessing |
| `OneHotEncoder` creates too many features | High-dimensional sparse matrix | Use `min_frequency=0.01` to group rare categories |
| `KNNImputer` on unscaled data | Poor imputation | Scale before imputing |
| `GridSearchCV` on entire parameter space | Search takes days | Use `RandomizedSearchCV` for > 5 params |
| Using `PCA` before train/test split | Data leakage | PCA in pipeline, fitted on training only |
| `stratify` parameter in `train_test_split` | Uneven class distribution in splits | Always `stratify=y` for classification |
| Not setting `random_state` | Non-reproducible results | Set `random_state=42` on every estimator |
| `joblib.load` from untrusted source | Code execution vulnerability | Only load models you trained |

---

## Reproducibility

```python
import numpy as np

# Set random state on every estimator
model = RandomForestClassifier(n_estimators=200, random_state=42)

# Set numpy seed for reproducibility in preprocessing
np.random.seed(42)

# Use the same seed in train_test_split and CV
from sklearn.model_selection import train_test_split, KFold
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
cv = KFold(n_splits=5, shuffle=True, random_state=42)
```

---

## See Also

- `references/experimental-campaign-protocol.md` — where this reference fits in the campaign workflow (Baseline phase)
- `references/pytorch-integration.md` — for deep learning methods
- `references/data-science-coding-workflow.md` — project structure, experiment logging
- scikit-learn.org/stable/user_guide — official user guide

# Statistical Methodology Reference

## Test Selection Decision Tree

```
What type of outcome (dependent variable)?
│
├─ CONTINUOUS (interval/ratio)
│  ├─ 1 group → One-sample t-test (normal) or Wilcoxon signed-rank (non-normal)
│  ├─ 2 independent groups → Independent t-test (normal, equal var) or Welch's t-test (normal, unequal var) or Mann-Whitney U (non-normal)
│  ├─ 2 paired groups → Paired t-test (normal) or Wilcoxon signed-rank (non-normal)
│  ├─ 3+ independent groups → One-way ANOVA (normal, equal var, independent) or Kruskal-Wallis (non-normal)
│  ├─ 3+ paired groups → Repeated measures ANOVA (normal, sphericity) or Friedman test (non-normal)
│  ├─ Controlling for covariates → ANCOVA (normal, equal slopes, independent)
│  └─ Multiple predictors → Linear regression / GLM
│
├─ BINARY (yes/no, success/failure)
│  ├─ 1 group / compare to known proportion → Binomial test / One-proportion z-test
│  ├─ 2 independent groups → Chi-square test of independence or Fisher's exact (small cells) or Two-proportion z-test
│  ├─ 2 paired groups → McNemar's test
│  ├─ 3+ groups → Chi-square test
│  ├─ Multiple predictors → Logistic regression
│  └─ Rare outcome → Logistic regression with Firth correction or exact methods
│
├─ COUNT (0, 1, 2, 3...)
│  ├─ Unbounded → Poisson regression (mean = variance) or Negative binomial (variance > mean)
│  ├─ Excess zeros → Zero-inflated Poisson/NB or Hurdle model
│  └─ Over time → Poisson/NB with offset for exposure
│
├─ ORDINAL (Likert, ranked)
│  ├─ Compare 2 groups → Mann-Whitney U (independent) or Wilcoxon signed-rank (paired)
│  ├─ Compare 3+ groups → Kruskal-Wallis
│  └─ Multiple predictors → Ordinal logistic regression (proportional odds)
│
├─ TIME-TO-EVENT (survival)
│  ├─ Compare 2 groups → Log-rank test
│  ├─ Multiple predictors → Cox proportional hazards
│  └─ Proportional hazards violated → Accelerated failure time models
│
├─ CATEGORICAL (nominal, 3+ levels)
│  ├─ Outcome is categorical → Chi-square test of independence or Multinomial logistic regression
│  └─ Agreement between raters → Cohen's kappa (2 raters) or Fleiss' kappa (3+)
│
├─ TIME SERIES (repeated over time)
│  ├─ Single series forecasting → ARIMA, Exponential smoothing, Prophet
│  ├─ Multiple series comparison → Structural time series, Dynamic regression
│  ├─ Seasonal patterns → Seasonal decomposition (STL), SARIMA
│  └─ Anomaly detection → Change point detection, Twitter/AnomalyDetection
│
└─ MULTIVARIATE (multiple outcomes)
    ├─ Dimensionality reduction → PCA, t-SNE, UMAP, Factor analysis
    ├─ Group structure → Cluster analysis (k-means, hierarchical, DBSCAN)
    └─ Multiple DVs by group → MANOVA (normal, equal cov matrices)
```

---

## Assumptions Reference Table

| Method | Assumptions | How to Check | What If Violated |
|--------|-------------|-------------|-------------------|
| **One-sample t-test** | Independence, normality | Q-Q plot, Shapiro-Wilk test (n < 50), Anderson-Darling | Wilcoxon signed-rank test |
| **Independent t-test** | Independence, normality, equal variance | Levene's test, F-test of variances, Q-Q plot | Welch's t-test (unequal var), Mann-Whitney (non-normal) |
| **Paired t-test** | Normality of differences, independence of pairs | Q-Q plot of differences | Wilcoxon signed-rank test |
| **One-way ANOVA** | Independence, normality (within groups), equal variance (homoscedasticity) | Levene's test, Shapiro-Wilk per group, Q-Q plot, residual plot | Welch's ANOVA (unequal var), Kruskal-Wallis (non-normal) |
| **Repeated measures ANOVA** | Normality, sphericity (equal variances of differences), independence between subjects | Mauchly's test for sphericity, ε correction (Greenhouse-Geisser, Huynh-Feldt) | Friedman test, or use mixed effects model with unstructured covariance |
| **Linear regression** | Linearity, independence of errors, homoscedasticity, normality of residuals, no multicollinearity | Residuals vs fitted plot, Q-Q plot, Breusch-Pagan test, VIF (< 5-10), Durbin-Watson, Cook's distance for influential points | Robust SEs (heteroscedasticity), weighted least squares, transformations (non-linearity), GLS (correlated errors) |
| **Logistic regression** | Linearity in logit (continuous predictors independent of log-odds), independence | Box-Tidwell test (linearity), Hosmer-Lemeshow goodness-of-fit, AUC-ROC, residual plots | Splines or polynomials (non-linearity), Firth regression (rare events/separation) |
| **Chi-square test** | Expected frequency ≥ 5 in each cell, independence of observations | Check expected counts | Fisher's exact test (2x2), Monte Carlo simulation (larger tables) |
| **Cox proportional hazards** | Proportional hazards, independent censoring | Schoenfeld residuals test (global + per covariate), log-log plots | Time-dependent covariates, stratified Cox, AFT models |
| **ANCOVA** | Normality, equal variance, independence, linear relationship with covariate, equal slopes assumption | Homogeneity of slopes test (covariate × group interaction) | Add interaction term, use nonparametric ANCOVA (Quade) |
| **MANOVA** | Multivariate normality, equal covariance matrices, independence | Box's M-test, Q-Q plots per variable | Separate ANOVAs with Bonferroni, PERMANOVA |

---

## Effect Size Guide

| Test / Design | Effect Size Measure | Interpretation | CI Available |
|--------------|-------------------|----------------|-------------|
| **t-test** (independent) | Cohen's d = (M₁−M₂)/s_pooled | 0.2=small, 0.5=medium, 0.8=large | Yes (non-central t) |
| **t-test** (paired) | Cohen's d_z = t/√n | Same conventions | Yes |
| **ANOVA** (one-way) | η² = SS_between/SS_total | 0.01=small, 0.06=medium, 0.14=large | Yes |
| **ANOVA** | Partial η² | Same (for multifactor designs) | Yes |
| **ANOVA** (fixed effects) | Cohen's f = √(η²/(1−η²)) | 0.10=small, 0.25=medium, 0.40=large | Yes |
| **Chi-square** | Cramér's V = √(χ²/(n·min(r−1,c−1))) | 0.1=small, 0.3=medium, 0.5=large | Yes (bootstrap) |
| **2×2 tables** | Odds ratio | OR=1 no effect, OR>1 increased odds | Yes (Woolf) |
| **2×2 tables** | Risk ratio / Relative risk | RR=1 no effect | Yes (Katz) |
| **2×2 tables** | Risk difference | Absolute difference in proportions | Yes (Newcombe) |
| **Correlation** | r (Pearson) | 0.1=small, 0.3=medium, 0.5=large | Yes (Fisher z) |
| **Correlation** | r_s (Spearman) | Same conventions | Yes |
| **Regression** | R² | Variance explained | Yes |
| **Regression** | Cohen's f² = R²/(1−R²) | 0.02=small, 0.15=medium, 0.35=large | Yes |
| **Bayesian** | Bayes factor BF₁₀ | 1-3=weak, 3-10=moderate, 10-30=strong, 30-100=very strong, >100=extreme evidence for H₁ | Yes (HDI) |

### Reporting Examples

- *"The mean difference was 3.2 points (95% CI [1.8, 4.6]), t(58) = 3.41, p = 0.001, Cohen's d = 0.87."*
- *"The treatment group had 2.3× the odds of recovery (OR 2.3, 95% CI [1.4, 3.8], χ²(1) = 12.4, p < 0.001)."*
- *"The model explained 34% of variance in outcome (R² = 0.34, F(3, 96) = 16.5, p < 0.001, Cohen's f² = 0.52)."*

---

## Multiple Testing Corrections

| Correction | Use When | How It Works | Power |
|-----------|---------|-------------|-------|
| **Bonferroni** | Small number of planned comparisons | α/m | Low — conservative |
| **Holm-Bonferroni** | Same, slightly less conservative | Sequential rejective | Better than Bonferroni |
| **Benjamini-Hochberg (BH)** | Exploratory analysis, many tests | Controls FDR | Higher — recommended for omics/exploratory |
| **Benjamini-Yekutieli (BY)** | Dependent tests, same as BH | Controls FDR under dependency | Lower than BH |
| **Tukey HSD** | All pairwise comparisons after ANOVA | Studentized range distribution | Good, designed for this case |
| **Dunnett** | Multiple comparisons vs a single control | Comparison-specific critical values | Good for this case |
| **Scheffé** | Post-hoc contrasts not planned in advance | Most flexible, most conservative | Low |
| **FDR (q-value)** | Thousands of tests (genomics, fMRI) | Estimates proportion of false discoveries | High |

**Rule of thumb:** For 2-5 planned comparisons, use Bonferroni or Holm. For dozens of exploratory tests, use BH at FDR=0.05. For pairwise post-ANOVA, use Tukey HSD. Never cherry-pick which p-values to correct.

---

## Bayesian Alternatives for Common Frequentist Tests

| Frequentist | Bayesian Alternative | Key Benefit |
|------------|---------------------|-------------|
| One-sample t-test | Bayesian one-sample t-test (BEST) | Can quantify evidence for H₀ via BF |
| Two-sample t-test | Bayesian two-sample t-test (BEST) | Robust to outliers via heavy-tailed likelihood |
| ANOVA | Bayesian ANOVA (BANOVA) | Model comparison via BFs, no sphericity assumption |
| Linear regression | Bayesian linear regression | Prior regularization, full posterior for coefficients |
| Logistic regression | Bayesian logistic regression (with priors) | Handles separation, shrinks extreme estimates |
| Chi-square test | Beta-Binomial model, contingency table Bayes factor | More intuitive: what's the posterior difference in proportions? |
| Correlation | Bayesian correlation (beta* prior on ρ) | Posterior distribution of ρ |
| t-test with non-inferiority | Bayesian region of practical equivalence (ROPE) | Direct probability of clinically meaningful difference |
| Meta-analysis | Bayesian hierarchical meta-analysis | Handles heterogeneity, small studies better |

---

## Best Practices & Red Flags

### Green Flags (good analysis practices)
- Pre-registered analysis plan (when possible)
- Effect sizes with CIs reported alongside p-values
- Assumption checks documented (and violations addressed)
- Sensitivity analyses reported (different specifications, outlier exclusion)
- Code and data available for reproduction
- Multiple testing corrections applied where appropriate
- Missing data mechanism discussed (MCAR, MAR, MNAR)
- Power analysis conducted before data collection

### Red Flags (poor analysis practices)
- p-values without effect sizes
- Stepwise variable selection (forward/backward)
- P-hacking: trying tests until significance
- Ignoring violations of assumptions
- No missing data handling (or claiming "no missing data" unrealistically)
- Over-interpreting non-significant results as "no effect" without equivalence testing
- Reporting only significant results (cherry-picking)
- Using parametric tests on clearly non-normal data without justification
- "p = 0.06 is marginally significant" (it's not — p = 0.06 is non-significant)

# Subagent Experiment Supervision

**When to load this reference:** Running a multi-experiment campaign (Phase 4-7 of the experimental protocol) and want automated failure recovery without manual intervention. Requires a subagent-capable harness (Hermes `delegate_task`, OpenCode subagents, or similar).

---

## Architecture

```
┌───────────────────────────────────────────────────┐
│                  Orchestrator                     │
│  (responsible for the overall experiment campaign) │
└────┬──────────────┬──────────────────┬────────────┘
     │              │                  │
     ▼              ▼                  ▼
┌──────────┐ ┌──────────┐     ┌──────────────┐
│ Worker 1 │ │ Worker 2 │ ... │ Supervisor   │
│ Runner   │ │ Runner   │     │ (one per     │
│          │ │          │     │  experiment) │
└──────────┘ └──────────┘     └──────┬───────┘
                                     │ watches logs
                                     │ auto-fixes
                                     │ escalates
```

**Key insight:** The supervisor is not a separate process. It's a lightweight monitoring loop that runs *alongside* the experiment, checking logs and applying known fixes. The orchestrator spawns one supervisor per experiment worker.

---

## Supervision Loop

```python
import re
import time
import subprocess
from pathlib import Path

FAILURE_PATTERNS = {
    "cuda_oom": {
        "pattern": r"CUDA out of memory",
        "fix": "reduce_batch_size",
        "priority": 1,
    },
    "cpu_oom": {
        "pattern": r"MemoryError|Cannot allocate memory",
        "fix": "reduce_memory",
        "priority": 2,
    },
    "nan_loss": {
        "pattern": r"loss.*nan|Loss.*NaN|nan.*loss",
        "fix": "gradient_clipping",
        "priority": 3,
    },
    "import_error": {
        "pattern": r"ModuleNotFoundError|ImportError.*No module named",
        "fix": "pip_install",
        "priority": 4,
    },
    "cuda_mismatch": {
        "pattern": r"CUDA error.*no kernel image|CUDA driver error",
        "fix": "fallback_cpu",
        "priority": 5,
    },
    "disk_full": {
        "pattern": r"No space left on device|Disk quota exceeded",
        "fix": "clean_disk",
        "priority": 6,
    },
    "timeout": {
        "pattern": r"TIMEOUT|killed|SIGTERM|SIGKILL",
        "fix": "reduce_scope",
        "priority": 7,
    },
}
```

### Supervisor Function

```python
def supervise_experiment(
    experiment_cmd: list,
    log_path: Path,
    max_retries: int = 3,
    check_interval: float = 5.0,
    timeout_hours: float = 24,
) -> dict:
    """
    Run an experiment with automated failure recovery.

    Args:
        experiment_cmd: Command to run (e.g., ["python", "train.py", "--config", "x.yaml"])
        log_path: Path to capture logs
        max_retries: Max consecutive auto-fix attempts before escalation
        check_interval: How often to check for errors (seconds)
        timeout_hours: Max wall time before considering the experiment hung

    Returns:
        Dict with status ("success", "fixed", "escalated", "timeout"), logs, fix_history
    """
    fix_history = []
    log_path.parent.mkdir(parents=True, exist_ok=True)
    start_time = time.time()

    for attempt in range(max_retries + 1):
        with open(log_path, "w") as log_file:
            process = subprocess.Popen(
                experiment_cmd,
                stdout=log_file,
                stderr=subprocess.STDOUT,
                text=True,
            )

            # Monitoring loop
            while True:
                # Check for timeout
                elapsed = (time.time() - start_time) / 3600
                if elapsed > timeout_hours:
                    process.kill()
                    return {
                        "status": "timeout",
                        "log_path": str(log_path),
                        "elapsed_hours": elapsed,
                    }

                # Check if process finished
                retcode = process.poll()
                if retcode is not None:
                    if retcode == 0:
                        return {
                            "status": "success" if attempt == 0 else "fixed",
                            "fix_history": fix_history,
                            "log_path": str(log_path),
                            "attempts": attempt,
                        }
                    break  # Process failed — try to diagnose

                time.sleep(check_interval)

        # Process failed — check logs for known patterns
        log_text = log_path.read_text()
        fix = diagnose_failure(log_text)

        if fix is None:
            # Unknown failure — escalate
            return {
                "status": "escalated",
                "log_path": str(log_path),
                "fix_history": fix_history,
                "error_snippet": log_text[-2000:],
            }

        # Apply fix
        fix_result = apply_fix(fix, experiment_cmd)
        fix_history.append({"attempt": attempt, "fix": fix, "result": fix_result})
        print(f"  [supervisor] Applied fix: {fix}")
```

### Failure Diagnosis

```python
def diagnose_failure(log_text: str) -> str | None:
    """Check logs against known failure patterns. Returns fix name or None."""
    for name, info in sorted(FAILURE_PATTERNS.items(), key=lambda x: x[1]["priority"]):
        if re.search(info["pattern"], log_text, re.IGNORECASE):
            return info["fix"]
    return None
```

---

## Failure Catalog with Fixes

| Signature | Detection (log pattern) | Fix |
|---|---|---|
| **CUDA OOM** | `CUDA out of memory` in stderr | Reduce `batch_size` by 50%, re-run |
| **CPU OOM** | `MemoryError` or system OOM killer | Halve data loading. Use `--data-fraction 0.5`. Re-run. |
| **NaN loss** | `loss: nan` or `Loss is NaN` or `nan in loss` | Add gradient clipping (`max_norm=1.0`). Reduce LR by 10x. Check for NaN in input data. Re-run. |
| **ImportError** | `ModuleNotFoundError: No module named 'X'` | `pip install X` and re-run |
| **CUDA version mismatch** | `CUDA error: no kernel image is available` or `CUDA driver version is insufficient` | Fall back to CPU: set `CUDA_VISIBLE_DEVICES=""` and re-run. Log the constraint. |
| **Disk full** | `No space left on device` or `Disk quota exceeded` | Clean temp files (`rm -rf /tmp/*.pt /tmp/__pycache__`). Alert user if < 1GB free. |
| **Timeout / Hung** | Process exceeds expected wall time with no log output for 30+ minutes | Kill process. Re-run with `--max-epochs 5 --max-steps 1000` (reduced scope). |
| **OOM during data loading** | `RuntimeError: DataLoader worker` + OOM | Reduce `num_workers` to 0. Set `persistent_workers=False`. Re-run. |
| **cuDNN init error** | `cuDNN error: CUDNN_STATUS_NOT_INITIALIZED` | Restart with fresh CUDA context. `torch.cuda.empty_cache()`. Re-run. |
| **Checkpoint corruption** | `RuntimeError: Error(s) in loading state_dict` | Remove corrupted checkpoint, restart from last known good epoch. |

---

## Fix Implementations

### Fix: Reduce Batch Size

```python
def fix_reduce_batch_size(cmd: list) -> list:
    """Modify command to use half the batch size."""
    new_cmd = cmd[:]
    for i, arg in enumerate(new_cmd):
        if arg == "--batch-size" or arg == "--batch_size":
            try:
                current = int(new_cmd[i + 1])
                new_cmd[i + 1] = str(max(1, current // 2))
                return new_cmd
            except (ValueError, IndexError):
                break
    # If no batch-size flag, append one
    new_cmd.extend(["--batch_size", "16"])
    return new_cmd
```

### Fix: Gradient Clipping

```python
def fix_gradient_clipping(cmd: list) -> list:
    """Add gradient clipping to the command."""
    if "--grad_clip" not in cmd and "--gradient-clip" not in cmd:
        cmd.extend(["--grad_clip", "1.0"])
    return cmd
```

### Fix: Reduce Scope (for timeouts)

```python
def fix_reduce_scope(cmd: list) -> list:
    """Limit epochs and data for a quick smoke test."""
    cmd.extend(["--max-epochs", "5", "--data-fraction", "0.1", "--quick-test"])
    return cmd
```

### Fix: Fall Back to CPU

```python
def fix_fallback_cpu(cmd: list) -> list:
    """Set environment to disable CUDA."""
    import os
    os.environ["CUDA_VISIBLE_DEVICES"] = ""
    return cmd  # Same command, but CUDA is now invisible
```

### Fix: Install Missing Package

```python
def fix_pip_install(log_text: str) -> str:
    """Extract missing module name from ImportError and install it."""
    match = re.search(
        r"ModuleNotFoundError: No module named '([^']+)'",
        log_text
    )
    if match:
        module = match.group(1)
        import subprocess
        subprocess.run(
            [sys.executable, "-m", "pip", "install", "--quiet", module],
            capture_output=True,
        )
        return f"installed {module}"
    return "unknown"
```

---

## Escalation Path

When 3 consecutive auto-fixes fail (or the supervisor encounters an unknown error):

1. **Save experiment context** — log file, fix history, snapshot of current state
2. **Send notification** to the user via Telegram (or configured alert channel)
3. **Include enough context** for the user to make a decision:

> **⚠️ Experiment supervision escalation**
>
> Experiment `run_007` failed after 3 auto-fix attempts.
>
> Last fix tried: `reduce_batch_size` (batch 64 → 32) — still OOM on a 24GB GPU.
>
> Next steps available:
> 1. Try gradient checkpointing (--gradient-checkpoint)
> 2. Use CPU offloading (--cpu-offload)
> 3. Abort this experiment, move to next candidate
>
> Log snippet: .../experiments/logs/run_007.log (last 50 lines attached)

**Telegram notification pattern** (when available):
```
python3 -c "
import urllib.request, json
urllib.request.urlopen(
    'https://api.telegram.org/bot<TOKEN>/sendMessage',
    data=json.dumps({
        'chat_id': '<CHAT_ID>',
        'text': '<message>',
        'parse_mode': 'Markdown'
    }).encode()
)
"
```

If Telegram is not configured, fall back to writing escalation to a file and continuing with the next experiment. The orchestrator should check for escalation files at the end of the campaign.

---

## Integration with the Campaign Protocol

### When to Use Supervision

| Protocol Phase | Supervision Value |
|---|---|
| Phase 2 (Baselines) | Low — baselines are fast enough to re-run manually |
| Phase 4 (Moonshots) | **High** — these are the most likely to fail |
| Phase 5 (Transfer Learning) | Medium — download failures, size mismatches |
| Phase 6 (HP Search) | **High** — 100+ trials, many will fail |
| Phase 7 (Distillation) | Medium — complex loss functions, convergence issues |

### What the Orchestrator Does

```python
# Pseudocode for the orchestrator's experiment loop
results = []
for candidate in shortlist:
    for trial in range(num_trials):
        supervisor = spawn_supervisor(
            experiment_cmd=["python", "train.py", "--config", candidate.config],
            log_path=f"logs/{candidate.name}_trial_{trial}.log",
        )
        result = supervisor.run()
        results.append(result)
        if result["status"] == "escalated":
            notify_user(result)
            # Continue with other candidates while waiting
```

---

## Harness-Specific Implementation Notes

### Hermes Agent (delegate_task)

```python
# The supervisor is spawned as a delegate_task that monitors the experiment
from hermes_tools import delegate_task  # if available

# Or: The orchestrator runs inline with subprocess monitoring
import subprocess, time

def run_supervised(experiment_cmd, log_path, max_retries=3):
    """Simple supervised experiment runner without subagent framework."""
    for attempt in range(max_retries + 1):
        with open(log_path, "w") as f:
            proc = subprocess.Popen(experiment_cmd, stdout=f, stderr=subprocess.STDOUT)

        while True:
            retcode = proc.poll()
            if retcode is not None:
                if retcode == 0:
                    return {"status": "success", "log_path": str(log_path)}
                break
            time.sleep(5)

        # Diagnose and fix (see functions above)
        log_text = Path(log_path).read_text()
        fix = diagnose_failure(log_text)
        if fix is None:
            return {"status": "escalated", "log_path": str(log_path)}
        experiment_cmd = apply_fix(fix, experiment_cmd)

    return {"status": "escalated", "log_path": str(log_path)}
```

### OpenCode / Claude Code

These harnesses don't have `delegate_task` but can use subprocess-based supervision. The same pattern applies — the agent runs the experiment as a subprocess and reads its logs periodically.

---

## Limitations

| Limitation | Mitigation |
|---|---|
| Only catches known failure patterns | The failure catalog is extensible — add patterns as you encounter them |
| Some fixes require code changes, not just CLI flags | For complex failures (architectural bugs), escalate immediately |
| Distributed training failures are more complex | DDP failures often require restarting the entire process group |
| Supervisor consumes monitoring overhead | Negligible (< 0.1% GPU) for the check_interval=5s pattern |
| Can't fix fundamental problems (bad architecture, wrong loss) | Escalate those — no auto-fix can rescue a fundamentally wrong approach |

---

## See Also

- `references/experimental-campaign-protocol.md` — the campaign workflow this supports
- `references/docker-experiment-isolation.md` — running experiments in containers
- `scripts/detect-compute.py` — know your hardware limits before scheduling

#!/usr/bin/env python3
"""
Model Assumption Diagnostics

Runs appropriate diagnostics on fitted models or raw data + method specs.
Returns structured report with warnings.

Usage:
  python assumption-diagnostics.py --method ttest --data data.csv --group-var group --value-var score
  python assumption-diagnostics.py --method regression --data data.csv --formula "y ~ x1 + x2"
  python assumption-diagnostics.py --method anova --data data.csv --group-var condition --value-var score
  python assumption-diagnostics.py --method mannwhitney --data data.csv --group-var group --value-var score
  python assumption-diagnostics.py ... --json
  python assumption-diagnostics.py ... --engine r
"""

import argparse
import json
import math
import sys

try:
    import numpy as np
    from scipy import stats as sp_stats
    HAS_NUMERIC = True
except ImportError:
    HAS_NUMERIC = False


def _check_deps():
    if not HAS_NUMERIC:
        print("Error: scipy and numpy required. pip install scipy numpy", file=sys.stderr)
        sys.exit(1)


METHODS = {
    "ttest": "Independent two-sample t-test",
    "ttest-paired": "Paired t-test",
    "onesample": "One-sample t-test",
    "anova": "One-way ANOVA",
    "regression": "Linear regression",
    "logistic": "Logistic regression",
    "correlation": "Pearson correlation",
    "mannwhitney": "Mann-Whitney U (nonparametric)",
    "kruskal": "Kruskal-Wallis (nonparametric)",
    "chisquare": "Chi-square test of independence",
}


def check_normality(data, method="shapiro"):
    """Test normality. Returns (statistic, p_value, is_violated)."""
    _check_deps()
    if method == "shapiro":
        if len(data) < 3:
            return None, None, True
        if len(data) > 5000:
            # Shapiro-Wilk is unreliable for n > 5000, use D'Agostino-Pearson
            stat, p = sp_stats.normaltest(data)
            method_used = "D'Agostino-Pearson"
        else:
            stat, p = sp_stats.shapiro(data)
            method_used = "Shapiro-Wilk"
    else:
        stat, p = sp_stats.normaltest(data)
        method_used = "D'Agostino-Pearson"

    return {
        "test": method_used,
        "statistic": round(stat, 4),
        "p_value": round(p, 4),
        "is_violated": p < 0.05
    }


def check_equal_variance(*groups):
    """Levene's test for equal variance across groups."""
    _check_deps()
    stat, p = sp_stats.levene(*groups)
    return {
        "test": "Levene's test",
        "statistic": round(stat, 4),
        "p_value": round(p, 4),
        "is_violated": p < 0.05
    }


def check_sphericity(data, groups, blocks):
    """Approximate sphericity check (Mauchly's test approximation)."""
    # Full sphericity requires repeated measures ANOVA structure
    return {
        "note": "Full Mauchly's test requires R (see --engine r). "
                "As a heuristic: check pairwise variance differences with Bartlett's test.",
        "is_violated": None
    }


def check_independence_durbin_watson(residuals):
    """Durbin-Watson test for autocorrelation of residuals."""
    _check_deps()
    n = len(residuals)
    dw = sum((residuals[i] - residuals[i-1])**2 for i in range(1, n)) / sum(r**2 for r in residuals)
    return {
        "test": "Durbin-Watson",
        "statistic": round(dw, 4),
        "is_violated": dw < 1.5 or dw > 2.5,
        "note": f"DW ≈ 2 = no autocorrelation. DW = {dw:.4f}"
    }


def check_linearity(x, y):
    """Check linearity via correlation ratio (eta) vs Pearson r."""
    _check_deps()
    r, _ = sp_stats.pearsonr(x, y)
    # Simple check: fit quadratic and see if it improves over linear
    # For now, report correlation and flag non-monotonic patterns
    return {
        "pearson_r": round(r, 4),
        "note": "For thorough linearity check, plot residuals vs fitted values.\n"
                "Significant non-linearity if residuals show clear U-shaped or curved pattern."
    }


def check_multicollinearity(X_matrix):
    """Approximate VIF for each predictor."""
    _check_deps()
    X = np.array(X_matrix)
    n_features = X.shape[1]
    vifs = []
    for i in range(n_features):
        y_i = X[:, i]
        X_i = np.delete(X, i, axis=1)
        try:
            # Regress feature i on all others, get R²
            X_i_with_intercept = np.column_stack([np.ones(X_i.shape[0]), X_i])
            beta = np.linalg.lstsq(X_i_with_intercept, y_i, rcond=None)[0]
            y_pred = X_i_with_intercept @ beta
            ss_res = np.sum((y_i - y_pred)**2)
            ss_tot = np.sum((y_i - np.mean(y_i))**2)
            r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 0
            vif = 1 / (1 - r2) if r2 < 1 else float('inf')
        except Exception:
            vif = float('inf')
        vifs.append(vif)
    return {
        "VIF_values": [round(v, 2) if v != float('inf') else "inf" for v in vifs],
        "high_collinearity": any(v > 10 for v in vifs if v != float('inf')),
        "note": "VIF > 5-10 indicates problematic multicollinearity."
    }


def check_outliers(data, method="iqr"):
    """Flag potential outliers."""
    _check_deps()
    data = np.array(data)
    q1, q3 = np.percentile(data, [25, 75])
    iqr = q3 - q1
    lower = q1 - 1.5 * iqr
    upper = q3 + 1.5 * iqr
    outliers = data[(data < lower) | (data > upper)]
    return {
        "method": "IQR (1.5×)",
        "n_outliers": len(outliers),
        "percent_outliers": round(len(outliers) / len(data) * 100, 1),
        "bounds": {"lower": round(lower, 4), "upper": round(upper, 4)},
        "is_violated": len(outliers) > 0,
        "note": f"{len(outliers)} potential outliers ({len(outliers)/len(data)*100:.1f}%)"
    }


def run_ttest_assumptions(group1, group2, is_paired=False):
    results = []

    n1, n2 = len(group1), len(group2)
    results.append({"check": "Sample size", "detail": f"n1 = {n1}, n2 = {n2}"})

    if not is_paired:
        results.append({"check": "Independence", "detail": "Design-based assumption (random assignment)",
                        "is_violated": False})

    # Normality per group
    norm1 = check_normality(group1)
    norm2 = check_normality(group2)
    if n1 < 30:
        results.append({"check": "Normality (Group 1)", "result": norm1,
                        "is_violated": norm1["is_violated"] if norm1 else True})
    if n2 < 30:
        results.append({"check": "Normality (Group 2)", "result": norm2,
                        "is_violated": norm2["is_violated"] if norm2 else True})

    if not is_paired and n1 >= 30 and n2 >= 30:
        results.append({"check": "Normality", "detail": "Both n ≥ 30 — CLT applies, normality not required"})

    if not is_paired:
        eqvar = check_equal_variance(group1, group2)
        results.append({"check": "Equal variance (Levene's)", "result": eqvar,
                        "is_violated": eqvar["is_violated"]})

    # Outliers per group
    out1 = check_outliers(group1)
    out2 = check_outliers(group2)
    results.append({"check": "Outliers (Group 1)", "result": out1})
    results.append({"check": "Outliers (Group 2)", "result": out2})

    if is_paired:
        diffs = np.array(group1) - np.array(group2)
        norm_diff = check_normality(diffs)
        results.append({"check": "Normality of differences", "result": norm_diff,
                        "is_violated": norm_diff["is_violated"] if norm_diff else True})

    passed = all(
        not r.get("is_violated", False)
        for r in results
        if "is_violated" in r and r["is_violated"] is not None
    )
    return {"method": "Independent t-test" if not is_paired else "Paired t-test",
            "overall_passed": passed,
            "checks": results,
            "recommendation": "All assumptions met" if passed else
                              "Violations detected. Consider: Welch's t-test (unequal var), "
                              "Mann-Whitney/Wilcoxon (non-normal), or check outliers."}


def run_regression_assumptions(X, y, residuals=None):
    """Run linear regression diagnostics."""
    results = []
    n = len(y)
    p = X.shape[1] if hasattr(X, 'shape') and len(X.shape) > 1 else 1
    results.append({"check": "Sample size", "detail": f"N = {n}, predictors = {p}, ratio = {n/p:.1f}:1",
                    "is_violated": n/p < 10})

    if residuals is not None:
        # Linearity: residuals vs fitted
        # Homoscedasticity: Breusch-Pagan approximation
        res = np.array(residuals)
        bp_stat = n * (sum(r**2 for r in res) / n) ** 2  # Simplified
        results.append({"check": "Residual normality", "result": check_normality(res)})

        # Durbin-Watson
        dw = check_independence_durbin_watson(res)
        results.append({"check": "Error independence (DW)", "result": dw, "is_violated": dw["is_violated"]})

        # Homoscedasticity via Breusch-Pagan simplified
        res2 = res ** 2
        bp_corr, _ = sp_stats.spearmanr(range(len(res2)), res2) if len(res2) > 3 else (0, 1)
        results.append({"check": "Homoscedasticity",
                        "detail": f"Spearman ρ between fitted values and |residuals| = {bp_corr:.4f}",
                        "is_violated": abs(bp_corr) > 0.15})

    results.append({"check": "Linearity", "detail": "Check residuals vs fitted plot for patterns"})

    passed = all(
        not r.get("is_violated", False)
        for r in results
        if "is_violated" in r and r["is_violated"] is not None
    )
    return {"method": "Linear regression",
            "overall_passed": passed,
            "checks": results,
            "recommendation": "All assumptions met" if passed else
                              "Violations detected. Consider: robust SEs (heteroscedasticity), "
                              "transformations (non-linearity), or GLS (correlated errors)."}


def run_anova_assumptions(groups):
    """Run one-way ANOVA diagnostics."""
    results = []
    n_groups = len(groups)
    sizes = [len(g) for g in groups]
    results.append({"check": "Sample sizes", "detail": str(sizes),
                    "is_violated": max(sizes) / min(sizes) > 2 if min(sizes) > 0 else True})

    # Normality per group (for small n)
    for i, g in enumerate(groups):
        if len(g) < 30:
            norm = check_normality(g)
            results.append({"check": f"Normality (Group {i+1})", "result": norm})

    # Equal variance
    eqvar = check_equal_variance(*groups)
    results.append({"check": "Equal variance (Levene's)", "result": eqvar,
                    "is_violated": eqvar["is_violated"]})

    # Independence
    results.append({"check": "Independence", "detail": "Design-based assumption (random assignment within blocks)"})

    passed = all(
        not r.get("is_violated", False)
        for r in results
        if "is_violated" in r and r["is_violated"] is not None
    )
    return {"method": "One-way ANOVA",
            "overall_passed": passed,
            "checks": results,
            "recommendation": "All assumptions met" if passed else
                              "Violations detected. Consider: Welch's ANOVA (unequal var), "
                              "Kruskal-Wallis (non-normal), or transform data."}


def run(args):
    # Parse data if provided
    if args.data:
        if not HAS_NUMERIC:
            print("Error: scipy + numpy required for data analysis. Install with: pip install scipy numpy",
                  file=sys.stderr)
            sys.exit(1)
        try:
            import pandas as pd
            df = pd.read_csv(args.data)
        except ImportError:
            print("Error: pandas required for CSV reading. pip install pandas", file=sys.stderr)
            sys.exit(1)
        except FileNotFoundError:
            print(f"Error: file not found: {args.data}", file=sys.stderr)
            sys.exit(1)
    else:
        df = None

    if not args.method:
        print("Error: --method required. Options: " + ", ".join(METHODS.keys()), file=sys.stderr)
        sys.exit(1)

    method_map = {
        "ttest": run_ttest_assumptions,
        "ttest-paired": lambda g1, g2: run_ttest_assumptions(g1, g2, is_paired=True),
    }

    result = {"method": args.method, "status": "ok"}

    if args.method in ("ttest", "ttest-paired") and df is not None:
        if not args.group_var or not args.value_var:
            print("Error: --group-var and --value-var required for t-test", file=sys.stderr)
            sys.exit(1)
        groups = [group[args.value_var].values for name, group in df.groupby(args.group_var)]
        if len(groups) != 2:
            print("Error: t-test requires exactly 2 groups", file=sys.stderr)
            sys.exit(1)
        result["diagnostics"] = run_ttest_assumptions(groups[0], groups[1], args.method == "ttest-paired")

    elif args.method == "anova" and df is not None:
        if not args.group_var or not args.value_var:
            print("Error: --group-var and --value-var required for ANOVA", file=sys.stderr)
            sys.exit(1)
        groups = [group[args.value_var].values for name, group in df.groupby(args.group_var)]
        result["diagnostics"] = run_anova_assumptions(groups)

    elif args.method == "regression":
        if df is not None and args.formula:
            try:
                import statsmodels.api as sm
                import statsmodels.formula.api as smf
                model = smf.ols(args.formula, data=df).fit()
                X = model.model.exog
                y = model.model.endog
                residuals = model.resid
                result["diagnostics"] = run_regression_assumptions(X, y, residuals)
                result["model_summary"] = {
                    "R_squared": round(model.rsquared, 4),
                    "adj_R_squared": round(model.rsquared_adj, 4),
                    "F_statistic": round(model.fvalue, 2),
                    "F_p_value": round(model.f_pvalue, 4),
                    "AIC": round(model.aic, 2),
                    "BIC": round(model.bic, 2),
                }
            except ImportError:
                print("Error: statsmodels required for regression diagnostics. pip install statsmodels",
                      file=sys.stderr)
                sys.exit(1)
        else:
            result["diagnostics"] = {
                "method": "Linear regression",
                "note": "Provide --data and --formula for full diagnostics",
                "checks": [
                    {"check": "Linearity", "detail": "Check residuals vs fitted plot"},
                    {"check": "Independence", "detail": "Check Durbin-Watson"},
                    {"check": "Homoscedasticity", "detail": "Check Breusch-Pagan test"},
                    {"check": "Normality of residuals", "detail": "Check Q-Q plot"},
                ]
            }

    elif args.method in ("mannwhitney", "kruskal", "chisquare", "onesample", "correlation", "logistic"):
        result["diagnostics"] = {
            "method": args.method,
            "note": f"Nonparametric and special methods require fewer assumptions. "
                    f"Provide data with --data, --group-var, --value-var for automated checks.",
            "checks": [
                {"check": "Independence", "detail": "Design-based assumption"},
            ]
        }

    else:
        print(f"Error: method '{args.method}' requires --data file. Provide CSV data.", file=sys.stderr)
        sys.exit(1)

    if args.json:
        print(json.dumps(result, indent=2, default=str))
    else:
        diag = result.get("diagnostics", {})
        print(f"## {diag.get('method', args.method)} Assumption Diagnostics")
        print(f"Status: {'✓ PASS' if diag.get('overall_passed', False) else '⚠ ISSUES FOUND'}")
        print()
        for check in diag.get("checks", []):
            status = "✓" if not check.get("is_violated") else "✗"
            print(f"{status} {check.get('check', 'Check')}")
            if "detail" in check:
                print(f"    {check['detail']}")
            if "result" in check and isinstance(check["result"], dict):
                for k, v in check["result"].items():
                    if k == "is_violated":
                        continue
                    print(f"    {k}: {v}")
            print()
        if "recommendation" in diag:
            print(f"**Recommendation:** {diag['recommendation']}")
        if "model_summary" in result:
            print(f"\nModel: R² = {result['model_summary']['R_squared']}, "
                  f"AIC = {result['model_summary']['AIC']}")
        if args.engine == "r":
            print("\n--- R equivalent ---")
            print(f"# In R, use: install.packages(c('car', 'lmtest', 'performance'))")
            print(f"library(car); library(lmtest); library(performance)")
            print(f"model <- lm({args.formula if args.formula else 'y ~ x'}, data = {args.data if args.data else 'df'})")
            print(f"check_model(model)  # Comprehensive assumptions plot")


def main():
    parser = argparse.ArgumentParser(description="Model Assumption Diagnostics")
    parser.add_argument("--method", choices=list(METHODS.keys()), help="Statistical method")
    parser.add_argument("--data", help="CSV file path")
    parser.add_argument("--group-var", help="Grouping variable name (for t-test, ANOVA)")
    parser.add_argument("--value-var", help="Value/outcome variable name")
    parser.add_argument("--formula", help="R-style formula for regression (e.g., 'y ~ x1 + x2')")
    parser.add_argument("--json", action="store_true", help="Output as JSON")
    parser.add_argument("--engine", choices=["python", "r"], default="python",
                        help="Output language")
    args = parser.parse_args()
    run(args)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
detect-compute.py — Probe hardware and software environment for ML feasibility.

Outputs structured recommendations so the agent can self-constrain its approach
based on available compute. Run before any experiment campaign to determine
what model sizes, batch sizes, and techniques are feasible.

Usage:
  python detect-compute.py                # Pretty-printed system overview
  python detect-compute.py --json          # Machine-readable JSON output
  python detect-compute.py --minimal       # Only the recommendations object
  python detect-compute.py --verbose       # Show every probe and its result
  python detect-compute.py --list-gpus     # Quick GPU inventory only

Exit codes:
  0 — Success
  1 — Probe completed but with warnings or degraded environment
"""

import json
import os
import platform
import shutil
import subprocess
import sys
import warnings

# ── Argument parsing ──────────────────────────────────────────────
FLAGS = {
    "json": False,
    "minimal": False,
    "verbose": False,
    "list_gpus": False,
}

for arg in sys.argv[1:]:
    if arg == "--json":
        FLAGS["json"] = True
    elif arg == "--minimal":
        FLAGS["minimal"] = True
    elif arg == "--verbose":
        FLAGS["verbose"] = True
    elif arg == "--list-gpus":
        FLAGS["list_gpus"] = True
    elif arg in ("-h", "--help"):
        print(__doc__.strip())
        sys.exit(0)


# ── Probe functions ───────────────────────────────────────────────

def _run(cmd: list[str], timeout: int = 15) -> tuple[str, str, int]:
    """Run a subprocess, return (stdout, stderr, exit_code)."""
    try:
        proc = subprocess.run(cmd, capture_output=True, text=True, timeout=timeout)
        return proc.stdout.strip(), proc.stderr.strip(), proc.returncode
    except FileNotFoundError:
        return "", f"command not found: {cmd[0]}", -1
    except subprocess.TimeoutExpired:
        return "", f"timed out after {timeout}s", -1


def _nvidia_smi() -> dict:
    """Parse nvidia-smi for GPU inventory. Returns empty dict if unavailable."""
    stdout, _, rc = _run(["nvidia-smi", "--query-gpu=index,name,memory.total,compute_cap",
                          "--format=csv,noheader,nounits"])
    if rc != 0:
        return {}

    gpus = []
    for line in stdout.strip().split("\n"):
        line = line.strip()
        if not line:
            continue
        parts = [p.strip() for p in line.split(", ")]
        if len(parts) >= 4:
            try:
                gpus.append({
                    "index": int(parts[0]),
                    "name": parts[1],
                    "vram_mb": int(float(parts[2])),
                    "compute_capability": parts[3],
                })
            except (ValueError, IndexError):
                continue

    if not gpus:
        return {}

    # Get driver version
    drv_stdout, _, _ = _run(["nvidia-smi", "--query-gpu=driver_version",
                             "--format=csv,noheader,nounits"])
    driver_version = drv_stdout.strip().split("\n")[0].strip() if drv_stdout else ""

    return {"gpu_count": len(gpus), "gpus": gpus, "driver_version": driver_version}


def _cuda_version() -> str:
    """Detect CUDA version from nvcc or nvidia-smi."""
    stdout, _, rc = _run(["nvcc", "--version"])
    if rc == 0:
        for line in stdout.split("\n"):
            if "release" in line:
                parts = line.split("release ")
                if len(parts) > 1:
                    return parts[1].split(",")[0].strip()
    # Fallback: try nvidia-smi topo
    stdout, _, _ = _run(["nvidia-smi"])
    for line in stdout.split("\n"):
        if "CUDA Version:" in line:
            return line.split("CUDA Version:")[-1].strip()
    return ""


def _torch_info() -> dict:
    """Probe PyTorch availability and capabilities via subprocess."""
    probe = r"""
import json, sys
try:
    import torch
    info = {
        "available": True,
        "version": torch.__version__,
        "cuda_available": torch.cuda.is_available(),
        "cuda_device_count": torch.cuda.device_count() if hasattr(torch.cuda, 'device_count') else 0,
        "mps_available": getattr(torch.backends, 'mps', None) is not None and torch.backends.mps.is_available(),
        "cuda_version": torch.version.cuda if hasattr(torch.version, 'cuda') else None,
    }
    if info["cuda_available"] and info["cuda_device_count"] > 0:
        info["current_device"] = torch.cuda.current_device()
        props = torch.cuda.get_device_properties(0)
        info["gpu_name"] = props.name
        info["vram_total_mb"] = props.total_memory // (1024 * 1024)
    print(json.dumps(info))
except Exception as e:
    print(json.dumps({"available": False, "error": str(e)}))
"""
    stdout, _, rc = _run([sys.executable, "-c", probe])
    if rc != 0 or not stdout:
        return {"available": False}
    try:
        return json.loads(stdout)
    except json.JSONDecodeError:
        return {"available": False}


def _sklearn_info() -> dict:
    """Probe scikit-learn availability and version."""
    probe = r"""
import json, sys
try:
    import sklearn
    print(json.dumps({"available": True, "version": sklearn.__version__}))
except Exception as e:
    print(json.dumps({"available": False, "error": str(e)}))
"""
    stdout, _, rc = _run([sys.executable, "-c", probe])
    if rc != 0 or not stdout:
        return {"available": False}
    try:
        return json.loads(stdout)
    except json.JSONDecodeError:
        return {"available": False}


def _get_ram_mb() -> int:
    """Get total physical RAM in MB."""
    try:
        import psutil
        return psutil.virtual_memory().total // (1024 * 1024)
    except ImportError:
        pass
    # Fallback: /proc/meminfo on Linux
    try:
        with open("/proc/meminfo") as f:
            for line in f:
                if line.startswith("MemTotal:"):
                    kb = int(line.split()[1])
                    return kb // 1024
    except (FileNotFoundError, ValueError, IndexError):
        pass
    # Fallback: sysctl on macOS
    stdout, _, _ = _run(["sysctl", "-n", "hw.memsize"])
    if stdout:
        try:
            return int(stdout.strip()) // (1024 * 1024)
        except ValueError:
            pass
    return 0


def _get_disk_free_mb(path: str = ".") -> int:
    """Get free disk space at path in MB."""
    try:
        import shutil
        _, _, free = shutil.disk_usage(path)
        return free // (1024 * 1024)
    except (ImportError, FileNotFoundError):
        pass
    # Fallback: df on Unix
    stdout, _, _ = _run(["df", "-P", path])
    for line in stdout.split("\n")[1:]:
        parts = line.strip().split()
        if len(parts) >= 4:
            try:
                return int(parts[3])  # Free in KB → convert
            except ValueError:
                pass
    return 0


def _has_jax() -> bool:
    """Check if JAX is available."""
    stdout, _, rc = _run([sys.executable, "-c", "import jax; print(jax.__version__)"])
    return rc == 0 and bool(stdout.strip())


def _has_optuna() -> bool:
    """Check if Optuna is available."""
    stdout, _, rc = _run([sys.executable, "-c", "import optuna; print(optuna.__version__)"])
    return rc == 0 and bool(stdout.strip())


# ── Recommendation engine ────────────────────────────────────────

def _recommendations(info: dict) -> dict:
    """Generate actionable recommendations based on detected hardware."""
    recs = {}
    vram_mb = 0

    # Get VRAM from the most reliable source
    torch_avail = info.get("torch", {}).get("available", False)
    if torch_avail and info["torch"].get("vram_total_mb"):
        vram_mb = info["torch"]["vram_total_mb"]
    elif info.get("nvidia", {}).get("gpus"):
        vram_mb = info["nvidia"]["gpus"][0].get("vram_mb", 0)

    has_cuda = info.get("nvidia", {}).get("gpu_count", 0) > 0
    has_torch = torch_avail
    has_sklearn = info.get("sklearn", {}).get("available", False)

    # Model size tier
    if vram_mb >= 24000:
        recs["model_size_tier"] = "13B-70B"
        recs["feasible_techniques"] = ["full_fine_tuning", "lora", "qlora", "distillation"]
    elif vram_mb >= 16000:
        recs["model_size_tier"] = "7B-13B"
        recs["feasible_techniques"] = ["full_fine_tuning", "lora", "qlora", "distillation"]
    elif vram_mb >= 8000:
        recs["model_size_tier"] = "3B-7B"
        recs["feasible_techniques"] = ["lora", "qlora", "distillation"]
        recs["notes"] = "Full fine-tuning may be tight for 7B. Prefer LoRA/QLoRA."
    elif vram_mb >= 4000:
        recs["model_size_tier"] = "up_to_3B"
        recs["feasible_techniques"] = ["qlora", "distillation"]
        recs["notes"] = "Full fine-tuning only for models <= 1.5B. Use QLoRA for larger."
    elif has_cuda:
        recs["model_size_tier"] = "up_to_1B"
        recs["feasible_techniques"] = ["qlora", "cpu_offloading"]
        recs["notes"] = "Limited VRAM. Consider cloud GPU or CPU-based methods."
    else:
        recs["model_size_tier"] = "cpu_only"
        recs["feasible_techniques"] = ["sklearn", "xgboost", "lightgbm"]
        recs["notes"] = "No GPU detected. Use sklearn/xgboost/lightgbm. No deep learning."

    # Batch size guidance
    if vram_mb >= 24000:
        recs["batch_size_guide"] = "LoRA: 128, Full FT: 32, Inference: 4096"
    elif vram_mb >= 16000:
        recs["batch_size_guide"] = "LoRA: 64, Full FT: 16, Inference: 2048"
    elif vram_mb >= 8000:
        recs["batch_size_guide"] = "LoRA: 32, Full FT: 8, Inference: 1024"
    elif vram_mb >= 4000:
        recs["batch_size_guide"] = "LoRA: 16, Full FT: 4, Inference: 512"
    elif has_cuda:
        recs["batch_size_guide"] = "LoRA: 8, Full FT: 2, Inference: 256"
    else:
        recs["batch_size_guide"] = "CPU-based. Batch size less relevant — use sklearn pipelines."

    # Quantization guidance
    if vram_mb >= 8000:
        recs["quantization_available"] = ["int8", "fp4", "fp8"]
    elif vram_mb >= 4000:
        recs["quantization_available"] = ["int8", "fp4"]
    elif has_cuda:
        recs["quantization_available"] = ["int8"]
    else:
        recs["quantization_available"] = []

    # Distillation
    recs["distillation_feasible"] = has_torch and vram_mb >= 4000

    # Fallback if no deep learning at all
    if not has_torch and not has_sklearn:
        recs["notes"] = "Neither PyTorch nor scikit-learn detected. Install: pip install torch scikit-learn"
    elif not has_torch and has_sklearn:
        recs["notes"] = recs.get("notes", "") + " PyTorch not found. sklearn/xgboost available."
    elif has_torch and not has_sklearn:
        recs["notes"] = recs.get("notes", "") + " scikit-learn not found. PyTorch available."

    return recs


# ── Main probe ────────────────────────────────────────────────────

def run_probes() -> dict:
    """Run all hardware and software probes, return structured results."""
    info = {
        "python_version": f"{sys.version_info.major}.{sys.version_info.minor}.{sys.version_info.micro}",
        "python_executable": sys.executable,
        "platform": sys.platform,
        "platform_detail": platform.platform(),
        "hostname": platform.node(),
    }

    # NVIDIA GPU probe
    nvidia = _nvidia_smi()
    info["nvidia"] = nvidia
    info["has_cuda"] = nvidia.get("gpu_count", 0) > 0
    info["cuda_version"] = _cuda_version() if info["has_cuda"] else None

    if FLAGS["verbose"] and info["has_cuda"]:
        info["_nvidia_smi_raw"] = _run(["nvidia-smi"])[0]

    # Torch
    info["torch"] = _torch_info()

    # sklearn
    info["sklearn"] = _sklearn_info()

    # JAX
    info["has_jax"] = _has_jax()

    # Optuna
    info["has_optuna"] = _has_optuna()

    # System resources
    ram_mb = _get_ram_mb()
    info["ram_mb"] = ram_mb
    info["ram_gb"] = round(ram_mb / 1024, 1) if ram_mb else 0

    disk_free_mb = _get_disk_free_mb()
    info["disk_free_mb"] = disk_free_mb
    info["disk_free_gb"] = round(disk_free_mb / 1024, 1) if disk_free_mb else 0

    # Recommendations
    info["recommendations"] = _recommendations(info)

    return info


def _format_verbose(info: dict) -> str:
    """Produce verbose human-readable output."""
    lines = []

    def kv(k: str, v: object) -> None:
        lines.append(f"  {k}: {v}")

    lines.append("── System ──────────────────────────────")
    kv("Python", info["python_version"])
    kv("Platform", info["platform"])
    kv("Host", info["hostname"])

    lines.append("\n── GPU ─────────────────────────────────")
    if info["has_cuda"]:
        for gpu in info["nvidia"]["gpus"]:
            kv(f"GPU {gpu['index']}", f"{gpu['name']} ({gpu['vram_mb']} MB VRAM, CC {gpu['compute_capability']})")
        kv("CUDA", info["cuda_version"] or "unknown")
        kv("Driver", info["nvidia"].get("driver_version", "unknown"))
    else:
        lines.append("  (none detected)")

    lines.append("\n── ML Frameworks ───────────────────────")
    t = info["torch"]
    if t.get("available"):
        kv("PyTorch", t["version"])
        kv("  CUDA avail", t.get("cuda_available", False))
        kv("  MPS avail", t.get("mps_available", False))
        if t.get("gpu_name"):
            kv("  Torch GPU", f"{t['gpu_name']} ({t.get('vram_total_mb', '?')} MB)")
    else:
        kv("PyTorch", "not installed")

    s = info["sklearn"]
    if s.get("available"):
        kv("scikit-learn", s["version"])
    else:
        kv("scikit-learn", "not installed")

    kv("JAX", "yes" if info["has_jax"] else "no")
    kv("Optuna", "yes" if info["has_optuna"] else "no")

    lines.append("\n── Resources ───────────────────────────")
    kv("RAM", f"{info['ram_gb']} GB" if info["ram_gb"] else "unknown")
    kv("Disk free", f"{info['disk_free_gb']} GB" if info["disk_free_gb"] else "unknown")

    lines.append("\n── Recommendations ─────────────────────")
    for k, v in info["recommendations"].items():
        lines.append(f"  {k}: {v}")

    return "\n".join(lines)


# ── Entry point ──────────────────────────────────────────────────

def main():
    info = run_probes()

    # Handle --list-gpus (fast path)
    if FLAGS["list_gpus"]:
        if info["has_cuda"]:
            for gpu in info["nvidia"]["gpus"]:
                print(f"GPU {gpu['index']}: {gpu['name']} ({gpu['vram_mb']} MB)")
        else:
            print("No NVIDIA GPUs detected")
        return

    # Handle output formats
    if FLAGS["json"]:
        if FLAGS["minimal"]:
            print(json.dumps(info["recommendations"], indent=2))
        else:
            print(json.dumps(info, indent=2, default=str))
    elif FLAGS["minimal"]:
        print(json.dumps(info["recommendations"], indent=2))
    elif FLAGS["verbose"]:
        print(_format_verbose(info))
    else:
        # Default: pretty human-readable, recommendations-focused
        print(json.dumps(info, indent=2, default=str))


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Effect Size Calculator

Computes effect sizes with confidence intervals for common designs.

Usage:
  python effect-size-calculator.py --design cohens-d --mean1 10 --mean2 8 --sd1 2.5 --sd2 2.8 --n1 30 --n2 30
  python effect-size-calculator.py --design cohens-d --t-stat 3.41 --df 58
  python effect-size-calculator.py --design eta-squared --ss-between 120 --ss-total 350
  python effect-size-calculator.py --design cramers-v --chi-sq 12.4 --n 200 --min-dim 2
  python effect-size-calculator.py --design odds-ratio --a 45 --b 55 --c 30 --d 70
  python effect-size-calculator.py --design correlation --r 0.45 --n 100
  python effect-size-calculator.py --design cohens-f2 --r-squared 0.34
  python effect-size-calculator.py --design r-squared --r 0.58
  python effect-size-calculator.py ... --json
  python effect-size-calculator.py ... --engine r
"""

import argparse
import json
import math
import sys

try:
    import numpy as np
    from scipy import stats as sp_stats
    HAS_NUMERIC = True
except ImportError:
    HAS_NUMERIC = False


DESIGNS = {
    "cohens-d": "Cohen's d (independent groups)",
    "cohens-d-paired": "Cohen's d_z (paired groups)",
    "hedges-g": "Hedges' g (corrected Cohen's d for small samples)",
    "eta-squared": "Eta-squared (ANOVA)",
    "partial-eta-squared": "Partial eta-squared (factorial ANOVA)",
    "cohens-f": "Cohen's f (ANOVA effect size)",
    "cohens-f2": "Cohen's f² (regression effect size)",
    "r-squared": "R-squared from correlation",
    "cramers-v": "Cramér's V (chi-square associations)",
    "odds-ratio": "Odds ratio (2×2 tables)",
    "risk-ratio": "Risk ratio / Relative risk (2×2 tables)",
    "risk-difference": "Risk difference (2×2 tables)",
    "pearson-r": "Pearson correlation",
    "cohens-h": "Cohen's h (arcsine transformation for proportions)",
    "glass-delta": "Glass's Δ (experimental SD as reference)",
}


def _check_deps():
    if not HAS_NUMERIC:
        print("Error: scipy required. pip install scipy numpy", file=sys.stderr)
        sys.exit(1)


def cohens_d_from_means(m1, m2, sd1, sd2, n1=None, n2=None):
    """Cohen's d with pooled SD."""
    pooled_sd = math.sqrt(((n1 - 1) * sd1**2 + (n2 - 1) * sd2**2) / (n1 + n2 - 2)) if n1 and n2 else math.sqrt((sd1**2 + sd2**2) / 2)
    d = (m1 - m2) / pooled_sd
    return d, pooled_sd


def cohens_d_from_t(t, df, n1, n2):
    """Cohen's d from t-statistic."""
    d = t * math.sqrt((n1 + n2) / (n1 * n2))
    return d


def hedges_g(d, n1, n2):
    """Convert Cohen's d to Hedges' g (small sample correction)."""
    df = n1 + n2 - 2
    correction = 1 - (3 / (4 * df - 1))
    return d * correction


def ci_cohens_d(d, n1, n2, alpha=0.05):
    """Confidence interval for Cohen's d using non-central t."""
    _check_deps()
    t_val = d * math.sqrt((n1 * n2) / (n1 + n2))
    df = n1 + n2 - 2
    try:
        # Non-central t confidence interval
        ncp_lower = sp_stats.nct.ppf(alpha / 2, df, t_val)
        ncp_upper = sp_stats.nct.ppf(1 - alpha / 2, df, t_val)
        d_lower = ncp_lower * math.sqrt((n1 + n2) / (n1 * n2))
        d_upper = ncp_upper * math.sqrt((n1 + n2) / (n1 * n2))
        return round(d_lower, 4), round(d_upper, 4)
    except Exception:
        # Fallback: delta method approximation
        se_d = math.sqrt((n1 + n2) / (n1 * n2) + d**2 / (2 * (n1 + n2)))
        z = sp_stats.norm.ppf(1 - alpha / 2)
        return round(d - z * se_d, 4), round(d + z * se_d, 4)


def eta_squared(ss_between, ss_total):
    """η² = SS_between / SS_total"""
    return ss_between / ss_total


def cohens_f_from_eta2(eta2):
    """f = sqrt(η² / (1 - η²))"""
    return math.sqrt(eta2 / (1 - eta2))


def cohens_f2_from_r2(r2):
    """f² = R² / (1 - R²)"""
    return r2 / (1 - r2)


def cramers_v(chi2, n, min_dim):
    """V = sqrt(χ² / (n * min(r-1, c-1)))"""
    return math.sqrt(chi2 / (n * min_dim))


def odds_ratio(a, b, c, d):
    """OR = (a/c) / (b/d) = ad / bc"""
    return (a * d) / (b * c)


def risk_ratio(a, b, c, d):
    """RR = (a/(a+b)) / (c/(c+d))"""
    p1 = a / (a + b) if (a + b) > 0 else 0
    p2 = c / (c + d) if (c + d) > 0 else 0
    return p1 / p2 if p2 > 0 else float('inf')


def risk_difference(a, b, c, d):
    """RD = a/(a+b) - c/(c+d)"""
    p1 = a / (a + b) if (a + b) > 0 else 0
    p2 = c / (c + d) if (c + d) > 0 else 0
    return p1 - p2


def ci_odds_ratio(a, b, c, d, alpha=0.05):
    """Woolf's CI for odds ratio."""
    _check_deps()
    or_val = odds_ratio(a, b, c, d)
    se_log_or = math.sqrt(1/a + 1/b + 1/c + 1/d)
    z = sp_stats.norm.ppf(1 - alpha / 2)
    lo = math.exp(math.log(or_val) - z * se_log_or)
    hi = math.exp(math.log(or_val) + z * se_log_or)
    return round(lo, 4), round(hi, 4)


def ci_pearson_r(r, n, alpha=0.05):
    """Fisher z-transformation CI for Pearson r."""
    _check_deps()
    z_r = 0.5 * math.log((1 + r) / (1 - r))
    se_z = 1 / math.sqrt(n - 3)
    z = sp_stats.norm.ppf(1 - alpha / 2)
    lo = math.tanh(z_r - z * se_z)
    hi = math.tanh(z_r + z * se_z)
    return round(lo, 4), round(hi, 4)


def cohens_h(p1, p2):
    """h = 2 * arcsin(sqrt(p1)) - 2 * arcsin(sqrt(p2))"""
    return 2 * math.asin(math.sqrt(p1)) - 2 * math.asin(math.sqrt(p2))


def r_squared_from_r(r):
    """R² = r²"""
    return r ** 2


def interpret_cohens_d(d):
    if abs(d) < 0.2:
        return "negligible"
    elif abs(d) < 0.5:
        return "small"
    elif abs(d) < 0.8:
        return "medium"
    else:
        return "large"


def run(args):
    result = {}

    if args.design == "cohens-d":
        if args.mean1 is not None and args.mean2 is not None and args.sd1 is not None and args.sd2 is not None and args.n1 and args.n2:
            d, pooled = cohens_d_from_means(args.mean1, args.mean2, args.sd1, args.sd2, args.n1, args.n2)
            ci = ci_cohens_d(d, args.n1, args.n2)
            result = {
                "effect_size": "Cohen's d",
                "d": round(d, 4),
                "pooled_sd": round(pooled, 4),
                "interpretation": interpret_cohens_d(d),
                "ci_95": ci,
                "parameters": {"mean1": args.mean1, "mean2": args.mean2, "n1": args.n1, "n2": args.n2}
            }
        elif args.t_stat is not None and args.df is not None:
            # Need n1 and n2 to compute d from t
            result = {"error": "t-statistic requires also --n1 and --n2 for Cohen's d"}
        else:
            result = {"error": "Provide --mean1 --mean2 --sd1 --sd2 --n1 --n2"}

    elif args.design == "hedges-g":
        if args.mean1 is not None and args.mean2 is not None and args.sd1 is not None and args.sd2 is not None and args.n1 and args.n2:
            d, pooled = cohens_d_from_means(args.mean1, args.mean2, args.sd1, args.sd2, args.n1, args.n2)
            g = hedges_g(d, args.n1, args.n2)
            result = {
                "effect_size": "Hedges' g",
                "g": round(g, 4),
                "cohens_d": round(d, 4),
                "interpretation": interpret_cohens_d(d),
                "parameters": {"mean1": args.mean1, "mean2": args.mean2, "n1": args.n1, "n2": args.n2}
            }
        else:
            result = {"error": "Provide --mean1 --mean2 --sd1 --sd2 --n1 --n2"}

    elif args.design == "eta-squared":
        if args.ss_between is not None and args.ss_total is not None:
            eta2 = eta_squared(args.ss_between, args.ss_total)
            f = cohens_f_from_eta2(eta2)
            result = {
                "effect_size": "Eta-squared",
                "eta_squared": round(eta2, 4),
                "cohens_f": round(f, 4),
                "interpretation": "small" if eta2 < 0.01 else ("medium" if eta2 < 0.06 else ("large" if eta2 < 0.14 else "very large")),
            }
        else:
            result = {"error": "Provide --ss-between and --ss-total"}

    elif args.design == "cohens-f2":
        if args.r_squared is not None:
            f2 = cohens_f2_from_r2(args.r_squared)
            result = {
                "effect_size": "Cohen's f²",
                "f_squared": round(f2, 4),
                "r_squared": args.r_squared,
                "interpretation": "small" if f2 < 0.02 else ("medium" if f2 < 0.15 else "large"),
            }
        else:
            result = {"error": "Provide --r-squared"}

    elif args.design == "cramers-v":
        if args.chi_sq is not None and args.n is not None and args.min_dim is not None:
            v = cramers_v(args.chi_sq, args.n, args.min_dim)
            result = {
                "effect_size": "Cramér's V",
                "V": round(v, 4),
                "parameters": {"chi_sq": args.chi_sq, "n": args.n, "min_dim": args.min_dim},
                "interpretation": "small" if v < 0.1 else ("medium" if v < 0.3 else "large"),
            }
        else:
            result = {"error": "Provide --chi-sq --n --min-dim"}

    elif args.design == "odds-ratio":
        if args.a is not None and args.b is not None and args.c is not None and args.d is not None:
            or_val = odds_ratio(args.a, args.b, args.c, args.d)
            ci = ci_odds_ratio(args.a, args.b, args.c, args.d)
            rr = risk_ratio(args.a, args.b, args.c, args.d)
            rd = risk_difference(args.a, args.b, args.c, args.d)
            result = {
                "effect_size": "Odds ratio",
                "OR": round(or_val, 4),
                "ci_95": ci,
                "RR": round(rr, 4),
                "risk_difference": round(rd, 4),
                "log_OR": round(math.log(or_val), 4),
                "interpretation": "OR = 1 (no effect)" if 0.95 <= or_val <= 1.05 else
                                 f"OR > 1 (increased odds)" if or_val > 1 else "OR < 1 (decreased odds)",
            }
        else:
            result = {"error": "Provide --a --b --c --d (2×2 table counts)"}

    elif args.design == "risk-ratio":
        if args.a is not None and args.b is not None and args.c is not None and args.d is not None:
            rr = risk_ratio(args.a, args.b, args.c, args.d)
            result = {
                "effect_size": "Risk ratio",
                "RR": round(rr, 4),
                "interpretation": "RR = 1 (no effect)" if 0.95 <= rr <= 1.05 else ("RR > 1 (increased risk)" if rr > 1 else "RR < 1 (decreased risk)"),
            }
        else:
            result = {"error": "Provide --a --b --c --d"}

    elif args.design == "risk-difference":
        if args.a is not None and args.b is not None and args.c is not None and args.d is not None:
            rd = risk_difference(args.a, args.b, args.c, args.d)
            result = {
                "effect_size": "Risk difference",
                "RD": round(rd, 4),
                "interpretation": f"Absolute difference of {abs(rd):.1%} in risk",
            }
        else:
            result = {"error": "Provide --a --b --c --d"}

    elif args.design in ("pearson-r", "correlation"):
        if args.r is not None and args.n is not None:
            ci = ci_pearson_r(args.r, args.n)
            result = {
                "effect_size": "Pearson r",
                "r": args.r,
                "r_squared": round(args.r ** 2, 4),
                "ci_95": ci,
                "n": args.n,
                "interpretation": "small" if abs(args.r) < 0.1 else ("medium" if abs(args.r) < 0.3 else "large"),
            }
        else:
            result = {"error": "Provide --r and --n"}

    elif args.design == "r-squared":
        if args.r is not None:
            r2 = r_squared_from_r(args.r)
            result = {
                "effect_size": "R-squared",
                "r_squared": round(r2, 4),
                "r": args.r,
            }
        else:
            result = {"error": "Provide --r"}

    elif args.design == "cohens-h":
        if args.p1 is not None and args.p2 is not None:
            h = cohens_h(args.p1, args.p2)
            result = {
                "effect_size": "Cohen's h",
                "h": round(h, 4),
                "parameters": {"p1": args.p1, "p2": args.p2},
                "interpretation": "small" if abs(h) < 0.2 else ("medium" if abs(h) < 0.5 else "large"),
            }
        else:
            result = {"error": "Provide --p1 and --p2"}

    else:
        result = {"error": f"Unknown design '{args.design}'"}

    if "error" in result:
        print(f"Error: {result['error']}", file=sys.stderr)
        sys.exit(1)

    if args.json:
        print(json.dumps(result, indent=2, default=str))
    else:
        print(f"## {result.get('effect_size', 'Effect Size')}")
        print(f"Value: {result.get(list(result.keys())[1], '?')}")
        if "ci_95" in result:
            print(f"95% CI: ({result['ci_95'][0]:.4f}, {result['ci_95'][1]:.4f})")
        if "interpretation" in result:
            print(f"Interpretation: {result['interpretation']}")
        if "r_squared" in result:
            print(f"R²: {result['r_squared']}")
        print()
        # Report template
        print("Reporting template:")
        if result.get("effect_size") == "Cohen's d":
            print(f"  d = {result['d']:.2f}, 95% CI [{result['ci_95'][0]:.2f}, {result['ci_95'][1]:.2f}]")
        elif result.get("effect_size") == "Odds ratio":
            print(f"  OR = {result['OR']:.2f}, 95% CI [{result['ci_95'][0]:.2f}, {result['ci_95'][1]:.2f}]")
        elif result.get("effect_size") == "Cramér's V":
            print(f"  V = {result['V']:.2f}")

        if args.engine == "r":
            print("\n--- R equivalent ---")
            print("library(effectsize)")
            if result.get("effect_size") == "Cohen's d":
                print(f"cohens_d({args.mean1}, {args.mean2}, pooled_sd = TRUE)" if args.mean1 else "# Provide data vectors")
            elif result.get("effect_size") == "Pearson r":
                print(f"library(psych); r.con(r = {args.r}, n = {args.n}, p = 0.95)")
            elif result.get("effect_size") == "Cramér's V":
                print(f"cramers_v(chi2 = {args.chi_sq}, n = {args.n}, nrow = {args.min_dim + 1})")


def main():
    parser = argparse.ArgumentParser(description="Effect Size Calculator")
    parser.add_argument("--design", choices=list(DESIGNS.keys()), required=True,
                        help="Effect size type")
    # Means/SDs
    parser.add_argument("--mean1", type=float)
    parser.add_argument("--mean2", type=float)
    parser.add_argument("--sd1", type=float)
    parser.add_argument("--sd2", type=float)
    parser.add_argument("--n1", type=int)
    parser.add_argument("--n2", type=int)
    parser.add_argument("--t-stat", type=float)
    parser.add_argument("--df", type=int)
    # ANOVA
    parser.add_argument("--ss-between", type=float)
    parser.add_argument("--ss-total", type=float)
    # Regression
    parser.add_argument("--r-squared", type=float)
    parser.add_argument("--r", type=float)
    # Categorical
    parser.add_argument("--chi-sq", type=float)
    parser.add_argument("--n", type=int)
    parser.add_argument("--min-dim", type=int, help="Cramér's V: min(rows-1, cols-1). For a 2×2 table pass 1; for 3×4 table pass 2 (min(2,3)=2)")
    parser.add_argument("--a", type=float)
    parser.add_argument("--b", type=float)
    parser.add_argument("--c", type=float)
    parser.add_argument("--d", type=float)
    # Proportions
    parser.add_argument("--p1", type=float)
    parser.add_argument("--p2", type=float)
    # General
    parser.add_argument("--json", action="store_true", help="Output as JSON")
    parser.add_argument("--engine", choices=["python", "r"], default="python",
                        help="Output language")
    args = parser.parse_args()
    run(args)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Experimental Design Generator

Generates experimental designs (randomization schedules) for common designs.
Supports Python default with --engine r for R output.

Usage:
  python experimental-design.py --design crd --treatments A B C --n-per-group 10
  python experimental-design.py --design rcbd --treatments Control Treatment --blocks 6 --n-per-block 1
  python experimental-design.py --design latin-square --treatments A B C D
  python experimental-design.py --design factorial --factors "temp:2:low,high" "pressure:2:100,200" --reps 3
  python experimental-design.py --design crossover --treatments A B C --sequences 3 --subjects 6
  python experimental-design.py --list-designs
  python experimental-design.py ... --json
  python experimental-design.py ... --output schedule.csv
  python experimental-design.py ... --engine r
"""

import argparse
import csv
import itertools
import json
import math
import random
import sys

try:
    import numpy as np
    HAS_NP = True
except ImportError:
    HAS_NP = False


DESIGNS = {
    "crd": "Completely Randomized Design",
    "rcbd": "Randomized Complete Block Design",
    "latin-square": "Latin Square Design",
    "factorial": "Factorial Design (full)",
    "crossover": "Crossover Design (2×2)",
    "split-plot": "Split-Plot Design",
}


def generate_crd(treatments, n_per_group):
    """Generate completely randomized design."""
    units = []
    for t in treatments:
        for i in range(n_per_group):
            units.append({"treatment": t, "unit_id": f"{t}_{i+1}"})
    random.shuffle(units)
    for i, u in enumerate(units, 1):
        u["run_order"] = i
    return units


def generate_rcbd(treatments, blocks, n_per_block=1):
    """Generate randomized complete block design."""
    units = []
    for b in range(1, blocks + 1):
        block_units = []
        for t in treatments:
            for r in range(n_per_block):
                block_units.append({"treatment": t, "block": b, "unit_id": f"B{b}_{t}_{r+1}"})
        random.shuffle(block_units)
        for j, u in enumerate(block_units, 1):
            u["run_within_block"] = j
        units.extend(block_units)
    return units


def generate_latin_square(treatments):
    """Generate Latin square design. Uses cyclic method for odd n, then randomize."""
    n = len(treatments)
    if n < 2:
        raise ValueError("Need at least 2 treatments")

    # Generate cyclic Latin square
    square = []
    for i in range(n):
        row = [treatments[(i + j) % n] for j in range(n)]
        square.append(row)

    # Randomize rows and columns
    random.shuffle(square)
    col_order = list(range(n))
    random.shuffle(col_order)

    units = []
    for i in range(n):
        for j in range(n):
            units.append({
                "row": i + 1,
                "column": j + 1,
                "treatment": square[i][col_order[j]],
                "unit_id": f"R{i+1}C{j+1}"
            })
    return units


def generate_factorial(factors, reps):
    """Generate full factorial design.
    factors: list of dicts with name, levels, level_labels
    """
    # Build factor levels
    factor_names = []
    factor_levels = []
    level_labels_list = []
    for f in factors:
        factor_names.append(f["name"])
        n_levels = int(f["levels"])
        labels = f.get("labels", "").split(",")
        # If explicit labels provided, use them; otherwise use coded levels
        if len(labels) == n_levels and labels[0]:
            level_labels_list.append(labels)
        else:
            level_labels_list.append([str(i+1) for i in range(n_levels)])
        factor_levels.append(list(range(n_levels)))

    # All combinations
    combos = list(itertools.product(*factor_levels))
    units = []
    rep_count = 1
    for rep in range(1, reps + 1):
        for combo in combos:
            entry = {"rep": rep}
            for i, name in enumerate(factor_names):
                entry[name] = level_labels_list[i][combo[i]]
            entry["unit_id"] = f"R{rep}_" + "_".join(str(level_labels_list[i][combo[i]]) for i in range(len(combo)))
            units.append(entry)

    random.shuffle(units)
    for i, u in enumerate(units, 1):
        u["run_order"] = i
    return units


def generate_crossover_2x2(treatments, subjects_per_seq):
    """Generate 2×2 crossover design."""
    if len(treatments) < 2:
        raise ValueError("Need at least 2 treatments for crossover")

    t = treatments[:2]
    sequences = [
        [t[0], t[1]],  # Sequence 1: A → B
        [t[1], t[0]],  # Sequence 2: B → A
    ]

    units = []
    subj_id = 1
    for seq_idx, seq in enumerate(sequences):
        for s in range(subjects_per_seq):
            for period, treatment in enumerate(seq, 1):
                units.append({
                    "subject": subj_id,
                    "sequence": seq_idx + 1,
                    "period": period,
                    "treatment": treatment,
                    "unit_id": f"S{subj_id}_P{period}"
                })
            subj_id += 1
    return units


def output_schedule(units, output_path=None, fmt="text"):
    """Output design schedule."""
    if not units:
        return

    if fmt == "json":
        output = json.dumps(units, indent=2)
    else:
        # Determine columns
        keys = list(units[0].keys())
        header = "\t".join(keys)
        rows = []
        for u in units:
            rows.append("\t".join(str(u.get(k, "")) for k in keys))
        output = header + "\n" + "\n".join(rows)

    if output_path:
        with open(output_path, "w") as f:
            if fmt == "json":
                json.dump(units, f, indent=2)
            else:
                # Use CSV for file output
                keys = list(units[0].keys())
                writer = csv.DictWriter(f, fieldnames=keys)
                writer.writeheader()
                writer.writerows(units)
        print(f"Schedule written to {output_path}")

    return output


def run(args):
    if args.list_designs:
        print("Available designs:\n")
        for key, desc in DESIGNS.items():
            print(f"  {key:20s} {desc}")
        return

    if not args.design:
        print("Error: --design required", file=sys.stderr)
        sys.exit(1)

    random.seed(args.seed) if args.seed else None

    units = []

    if args.design == "crd":
        if not args.treatments or not args.n_per_group:
            print("Error: --treatments and --n-per-group required for CRD", file=sys.stderr)
            sys.exit(1)
        units = generate_crd(args.treatments, args.n_per_group)

    elif args.design == "rcbd":
        if not args.treatments or not args.blocks:
            print("Error: --treatments and --blocks required for RCBD", file=sys.stderr)
            sys.exit(1)
        units = generate_rcbd(args.treatments, args.blocks, args.n_per_block or 1)

    elif args.design == "latin-square":
        if not args.treatments or len(args.treatments) < 2:
            print("Error: --treatments requires 2+ treatments for Latin square", file=sys.stderr)
            sys.exit(1)
        units = generate_latin_square(args.treatments)

    elif args.design == "factorial":
        if not args.factors:
            print("Error: --factors required for factorial (format: 'name:levels:label1,label2')", file=sys.stderr)
            sys.exit(1)
        factors = []
        for f_str in args.factors:
            parts = f_str.split(":")
            name = parts[0]
            levels = parts[1] if len(parts) > 1 else "2"
            labels = parts[2] if len(parts) > 2 else ""
            factors.append({"name": name, "levels": levels, "labels": labels})
        units = generate_factorial(factors, args.reps or 1)

    elif args.design == "crossover":
        if not args.treatments or len(args.treatments) < 2:
            print("Error: --treatments requires 2 treatments for crossover", file=sys.stderr)
            sys.exit(1)
        if not args.subjects:
            print("Error: --subjects required for crossover", file=sys.stderr)
            sys.exit(1)
        units = generate_crossover_2x2(args.treatments, args.subjects // 2)

    elif args.design == "split-plot":
        if not args.treatments:
            whole_plot = args.treatments[:2] if len(args.treatments) >= 2 else ["A", "B"]
            sub_plot = args.treatments[2:4] if len(args.treatments) >= 4 else ["C", "D"]
        else:
            whole_plot = ["WP1", "WP2"]
            sub_plot = ["SP1", "SP2"]
        n_reps = args.n_per_group or 3
        units = []
        for rep in range(1, n_reps + 1):
            for wp in whole_plot:
                for sp in sub_plot:
                    units.append({
                        "rep": rep,
                        "whole_plot": wp,
                        "sub_plot": sp,
                        "unit_id": f"R{rep}_{wp}_{sp}"
                    })
        random.shuffle(units)
        for i, u in enumerate(units, 1):
            u["run_order"] = i

    else:
        print(f"Error: unknown design '{args.design}'", file=sys.stderr)
        sys.exit(1)

    if not units:
        print("Error: no units generated", file=sys.stderr)
        sys.exit(1)

    # Summary
    n_total = len(units)
    n_treatments = len(set(u.get("treatment", u.get(units[0]["treatment"] if "treatment" in units[0] else "")) for u in units if "treatment" in u))

    summary = {
        "design": args.design,
        "n_units": n_total,
        "seed": args.seed,
    }

    if args.json:
        print(json.dumps({"summary": summary, "schedule": units}, indent=2, default=str))
    elif args.output:
        keys = list(units[0].keys())
        with open(args.output, "w", newline="") as f:
            writer = csv.DictWriter(f, fieldnames=keys)
            writer.writeheader()
            writer.writerows(units)
        print(f"Design: {DESIGNS.get(args.design, args.design)}")
        print(f"Schedule: {n_total} experimental units → {args.output}")
        print(f"Seed: {args.seed or 'random'}")
    else:
        print(f"Design: {DESIGNS.get(args.design, args.design)}")
        print(f"Experimental units: {n_total}")
        print(f"Seed: {args.seed or 'random'}")
        print()
        # Show first few rows
        keys = list(units[0].keys())
        print("\t".join(keys))
        for u in units[:min(20, len(units))]:
            print("\t".join(str(u.get(k, "")) for k in keys))
        if len(units) > 20:
            print(f"... and {len(units) - 20} more units")

    if args.engine == "r":
        print("\n--- R equivalent ---")
        if args.design == "crd":
            print("library(agricolae)")
            treatments_str = ", ".join(f'"{t}"' for t in args.treatments)
            print(f'treatments <- c({treatments_str})')
            print(f'design.crd(treatments, r = {args.n_per_group}, seed = {args.seed or 42})')
        elif args.design == "rcbd":
            treatments_str = ", ".join(f'"{t}"' for t in args.treatments)
            print(f'treatments <- c({treatments_str})')
            print(f'design.rcbd(treatments, r = {args.blocks}, seed = {args.seed or 42})')
        elif args.design == "latin-square":
            treatments_str = ", ".join(f'"{t}"' for t in args.treatments)
            print(f'treatments <- c({treatments_str})')
            print(f'design.lsd(treatments, seed = {args.seed or 42})')


def main():
    parser = argparse.ArgumentParser(description="Experimental Design Generator")
    parser.add_argument("--design", choices=list(DESIGNS.keys()), help="Experimental design type")
    parser.add_argument("--list-designs", action="store_true", help="List available designs")
    parser.add_argument("--treatments", nargs="+", help="Treatment level names")
    parser.add_argument("--n-per-group", type=int, help="Subjects per group (CRD) or reps")
    parser.add_argument("--blocks", type=int, help="Number of blocks (RCBD)")
    parser.add_argument("--n-per-block", type=int, default=1, help="Subjects per block per treatment")
    parser.add_argument("--factors", nargs="+", help='Factor specs for factorial: "name:levels:l1,l2"')
    parser.add_argument("--reps", type=int, default=1, help="Replicates per combination (factorial)")
    parser.add_argument("--subjects", type=int, help="Total subjects (crossover)")
    parser.add_argument("--seed", type=int, help="Random seed for reproducibility")
    parser.add_argument("--output", help="Output CSV file path")
    parser.add_argument("--json", action="store_true", help="Output as JSON")
    parser.add_argument("--engine", choices=["python", "r"], default="python",
                        help="Output language")
    args = parser.parse_args()
    run(args)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Model Comparison Tool

Compares multiple models using information criteria and cross-validation.
Accepts model fit statistics directly or fits models from data.

Usage:
  python model-comparison.py --models "OLS AIC=1200 BIC=1220 k=5 LL=-590" "GLM AIC=1190 BIC=1215 k=6 LL=-588"
  python model-comparison.py --from-data data.csv --formulas "y~x1" "y~x1+x2" "y~x1*x2"
  python model-comparison.py ... --json
  python model-comparison.py ... --engine r
"""

import argparse
import json
import math
import sys
import re


def compute_aic(ll, k):
    """AIC = -2*LL + 2*k"""
    return -2 * ll + 2 * k


def compute_bic(ll, k, n):
    """BIC = -2*LL + k*ln(n)"""
    return -2 * ll + k * math.log(n)


def compute_aicc(aic, k, n):
    """Corrected AIC for small samples."""
    return aic + (2 * k * (k + 1)) / (n - k - 1) if n > k + 1 else float('inf')


def compute_weights(values):
    """Compute Akaike weights or similar model probabilities from IC values."""
    if not values:
        return []
    min_val = min(values)
    rel_likelihoods = [math.exp(-0.5 * (v - min_val)) for v in values]
    total = sum(rel_likelihoods)
    return [rl / total for rl in rel_likelihoods]


def model_entry_from_string(s):
    """Parse 'Name AIC=1200 BIC=1220 k=5 LL=-590 [n=100]'"""
    parts = s.split()
    if not parts:
        return None
    entry = {"name": parts[0]}
    for p in parts[1:]:
        if "=" in p:
            key, val = p.split("=", 1)
            key = key.lower()
            if key in ("aic", "bic", "aicc", "ll"):
                entry[key] = float(val)
            elif key in ("k", "n"):
                entry[key] = int(float(val))
    return entry if len(entry) > 1 else None


def compute_from_formulas(data_file, formulas, family="gaussian"):
    """Fit models from R-style formulas and compute ICs."""
    try:
        import pandas as pd
        import statsmodels.api as sm
        import statsmodels.formula.api as smf
    except ImportError:
        print("Error: pandas and statsmodels required for --from-data. pip install pandas statsmodels",
              file=sys.stderr)
        sys.exit(1)

    try:
        df = pd.read_csv(data_file)
    except Exception as e:
        print(f"Error reading {data_file}: {e}", file=sys.stderr)
        sys.exit(1)

    n = len(df)
    models = []
    for formula in formulas:
        try:
            if family == "gaussian":
                model = smf.ols(formula, data=df).fit()
            elif family == "binomial":
                model = smf.logit(formula, data=df).fit(disp=0)
            else:
                model = smf.ols(formula, data=df).fit()

            k = model.df_model + 1  # +1 for intercept
            ll = model.llf
            aic = model.aic
            bic = model.bic
            models.append({
                "name": formula,
                "aic": round(aic, 2),
                "bic": round(bic, 2),
                "k": int(k),
                "ll": round(ll, 2),
                "n": n,
                "r_squared": round(getattr(model, "rsquared", 0), 4),
            })
        except Exception as e:
            print(f"Warning: model '{formula}' failed: {e}", file=sys.stderr)
            continue

    return models


def compute_loo_cv(models_data, data_file):
    """Approximate LOO-CV using information criteria (not full CV)."""
    # WAIC and LOO require full posterior samples
    # For frequentist models, AIC/based weights serve as approximations
    return {
        "note": "Full LOO-CV requires Bayesian models with MCMC posterior samples. "
                "For frequentist models, AICc weights or k-fold CV are recommended instead.",
        "alternative": "Use k-fold cross-validation or Bayesian model (PyMC/Stan/NumPyro) with loo package."
    }


def run(args):
    models = []

    if args.models:
        for s in args.models:
            entry = model_entry_from_string(s)
            if entry:
                models.append(entry)

    if args.from_data:
        if not args.formulas:
            print("Error: --formulas required with --from-data", file=sys.stderr)
            sys.exit(1)
        fitted = compute_from_formulas(args.from_data, args.formulas, args.family)
        models.extend(fitted)

    if not models:
        print("Error: no valid models provided. Use --models or --from-data with --formulas",
              file=sys.stderr)
        print("\nExample:")
        print('  python model-comparison.py --models "OLS AIC=1200 BIC=1220 k=5 LL=-590" \\')
        print('    "GLM AIC=1190 BIC=1215 k=6 LL=-588"')
        print('  python model-comparison.py --from-data data.csv --formulas "y~x1" "y~x1+x2" \\')
        print('    --formulas "y~I(x1^2)"')
        sys.exit(1)

    # Ensure all models have required fields
    found_n = args.n or max(m.get("n", 0) for m in models) if models else 0
    for m in models:
        if not m.get("n") and found_n:
            m["n"] = found_n

    # Compute derived metrics
    for m in models:
        if "aic" not in m and "ll" in m and "k" in m:
            m["aic"] = round(compute_aic(m["ll"], m["k"]), 2)
        if "bic" not in m and "ll" in m and "k" in m and m.get("n"):
            m["bic"] = round(compute_bic(m["ll"], m["k"], m["n"]), 2)
        if "aic" in m and m.get("k") and m.get("n"):
            m["aicc"] = round(compute_aicc(m["aic"], m["k"], m["n"]), 2)
            if m["aicc"] == float('inf'):
                m["aicc"] = None

    # Sort by AIC (lower is better)
    models.sort(key=lambda m: m.get("aic", float('inf')))

    # Compute Akaike weights and delta
    aic_vals = [m.get("aic", float('inf')) for m in models]
    if all(v != float('inf') for v in aic_vals):
        weights = compute_weights(aic_vals)
        for i, m in enumerate(models):
            m["delta_aic"] = round(aic_vals[i] - aic_vals[0], 2)
            m["akaike_weight"] = round(weights[i], 4)

    bic_vals = [m.get("bic", float('inf')) for m in models]
    if all(v != float('inf') for v in bic_vals):
        bic_weights = compute_weights(bic_vals)
        for i, m in enumerate(models):
            m["delta_bic"] = round(bic_vals[i] - bic_vals[0], 2)
            m["bic_weight"] = round(bic_weights[i], 4)

    result = {
        "n_models": len(models),
        "criterion": "AIC / BIC / AICc",
        "best_model": models[0]["name"] if models else None,
        "models": models,
    }

    # Ranking table
    if args.json:
        print(json.dumps(result, indent=2, default=str))
    else:
        print("Model Comparison")
        print("=" * 80)
        print(f"Best model: {result['best_model']}")
        print()
        header = f"{'Model':<25} {'AIC':>8} {'ΔAIC':>8} {'w(AIC)':>8} {'BIC':>8} {'ΔBIC':>8} {'k':>4} {'R²':>6}"
        if any(m.get("aicc") for m in models if m.get("aicc") is not None):
            header += f" {'AICc':>8}"
        print(header)
        print("-" * len(header))
        for m in models:
            row = f"{m['name']:<25} "
            row += f"{m.get('aic', '-'):>8.2f} " if isinstance(m.get('aic'), (int, float)) else f"{'NA':>8} "
            row += f"{m.get('delta_aic', '-'):>8.2f} " if isinstance(m.get('delta_aic'), (int, float)) else f"{'NA':>8} "
            row += f"{m.get('akaike_weight', '-'):>8.4f} " if isinstance(m.get('akaike_weight'), (int, float)) else f"{'NA':>8} "
            row += f"{m.get('bic', '-'):>8.2f} " if isinstance(m.get('bic'), (int, float)) else f"{'NA':>8} "
            row += f"{m.get('delta_bic', '-'):>8.2f} " if isinstance(m.get('delta_bic'), (int, float)) else f"{'NA':>8} "
            row += f"{m.get('k', '-'):>4} "
            row += f"{m.get('r_squared', '-'):>6.4f} " if isinstance(m.get('r_squared'), (int, float)) else f"{'NA':>6} "
            if any(mm.get("aicc") for mm in models if mm.get("aicc") is not None):
                aicc_val = m.get('aicc', '-')
                if aicc_val is not None and isinstance(aicc_val, (int, float)):
                    row += f"{aicc_val:>8.2f}"
                else:
                    row += f"{'NA':>8}"
            print(row)

        print()
        print("Interpretation:")
        print("  ΔAIC < 2: substantial support for model being best")
        print("  ΔAIC 4-7: considerably less support")
        print("  ΔAIC > 10: essentially no support")
        print("  w(AIC): probability that model is best (given set)")
        print()
        if args.from_data:
            loo_note = compute_loo_cv(models, args.from_data)
            print(f"Note: {loo_note['note']}")

        if args.engine == "r":
            print("\n--- R equivalent ---")
            print("# Compare models in R")
            print("library(AICcmodavg)")
            print("# Assuming model objects: m1, m2, m3")
            print("models <- list(m1=m1, m2=m2, m3=m3)")
            print("aictab(models)")
            print("print(bictab(models))")


def main():
    parser = argparse.ArgumentParser(description="Model Comparison Tool")
    parser.add_argument("--models", nargs="+", help='Model specs: "Name AIC=1200 BIC=1220 k=5 LL=-590"')
    parser.add_argument("--from-data", help="CSV file to fit models from")
    parser.add_argument("--formulas", nargs="+", help="R-style formulas for model fitting")
    parser.add_argument("--family", default="gaussian", choices=["gaussian", "binomial"],
                        help="Distribution family for model fitting")
    parser.add_argument("--n", type=int, help="Sample size (if not in model specs)")
    parser.add_argument("--json", action="store_true", help="Output as JSON")
    parser.add_argument("--engine", choices=["python", "r"], default="python",
                        help="Output language")
    args = parser.parse_args()
    run(args)


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""
Power Analysis Calculator

Computes sample size from effect size (and vice versa) for common designs.
Supports Python default with --engine r for R output.

Usage:
  python power-analysis.py --design ttest-ind --effect-size 0.5 --alpha 0.05 --power 0.80
  python power-analysis.py --design ttest-paired --n-per-group 30 --alpha 0.05 --power 0.80
  python power-analysis.py --design anova --k 3 --effect-size 0.25 --alpha 0.05 --power 0.80
  python power-analysis.py --design prop --p1 0.10 --p2 0.15 --alpha 0.05 --power 0.80
  python power-analysis.py --design correlation --effect-size 0.3 --alpha 0.05 --power 0.80
  python power-analysis.py --design regression --predictors 5 --effect-size 0.15 --alpha 0.05 --power 0.80
  python power-analysis.py --design chi-square --df 2 --effect-size 0.3 --alpha 0.05 --power 0.80
  python power-analysis.py --design equivalence --effect-size 0.5 --alpha 0.05 --power 0.80
  python power-analysis.py --list-designs
  python power-analysis.py ... --json
  python power-analysis.py ... --engine r
"""

import argparse
import math
import json
import sys

try:
    import numpy as np
    from scipy import stats as sp_stats
    HAS_NUMERIC = True
except ImportError:
    HAS_NUMERIC = False


def _check_deps():
    if not HAS_NUMERIC:
        print("Error: scipy and numpy are required. Install with: pip install scipy numpy",
              file=sys.stderr)
        sys.exit(1)


DESIGNS = {
    "ttest-ind": "Two-sample independent t-test (equal n per group)",
    "ttest-ind-unequal": "Two-sample independent t-test (unequal n, specify ratio)",
    "ttest-paired": "Paired t-test",
    "onesample": "One-sample t-test",
    "prop": "Two-proportion z-test",
    "onesample-prop": "One-sample proportion test",
    "anova": "One-way ANOVA (k groups, equal n per group)",
    "anova-interaction": "ANOVA interaction effect (2×2 factorial)",
    "correlation": "Pearson correlation test",
    "regression": "Multiple linear regression (F-test for R²)",
    "logistic": "Logistic regression (Wald test for single coefficient)",
    "chi-square": "Chi-square test of independence (contingency table)",
    "equivalence": "Two one-sided tests (TOST) for equivalence",
    "survival": "Survival analysis (log-rank test)",
}


def solve_power_ttest_ind(d, alpha=0.05, power=0.80, ratio=1.0, alternative="two-sided"):
    """Compute per-group sample size for independent t-test. Returns n_per_group."""
    _check_deps()
    if alternative == "two-sided":
        alpha /= 2
    z_beta = sp_stats.norm.ppf(power)
    z_alpha = sp_stats.norm.ppf(1 - alpha)
    n_per_group = ((z_alpha + z_beta) ** 2 * (1 + 1/ratio) / (d ** 2)) + 2
    return int(math.ceil(n_per_group))


def solve_power_ttest_paired(d, alpha=0.05, power=0.80):
    """Compute number of pairs for paired t-test."""
    _check_deps()
    z_beta = sp_stats.norm.ppf(power)
    z_alpha = sp_stats.norm.ppf(1 - alpha / 2)
    n = ((z_alpha + z_beta) ** 2 / (d ** 2)) + 2
    return int(math.ceil(n))


def solve_power_onesample(d, alpha=0.05, power=0.80):
    """Compute sample size for one-sample t-test."""
    return solve_power_ttest_paired(d, alpha, power)


def solve_power_prop(p1, p2, alpha=0.05, power=0.80, ratio=1.0):
    """Compute per-group sample size for two-proportion z-test."""
    _check_deps()
    p_bar = (p1 + ratio * p2) / (1 + ratio)
    z_beta = sp_stats.norm.ppf(power)
    z_alpha = sp_stats.norm.ppf(1 - alpha / 2)
    n = ((z_alpha + z_beta) ** 2 *
         (p1 * (1 - p1) / 1 + p2 * (1 - p2) / ratio)) / ((p1 - p2) ** 2)
    return int(math.ceil(n))


def solve_power_anova(f, k, alpha=0.05, power=0.80):
    """Compute per-group sample size for one-way ANOVA using non-central F distribution.

    More accurate than the normal approximation. Uses iterative search.
    f = Cohen's f = sqrt(η² / (1 - η²))
    Non-centrality parameter λ = n * k * f²
    df1 = k - 1,  df2 = k * (n - 1)
    """
    _check_deps()
    f2 = f ** 2

    def _power_at_n(n):
        df1 = k - 1
        df2 = k * (n - 1)
        if df2 < 1:
            return 0.0
        f_crit = sp_stats.f.ppf(1 - alpha, df1, df2)
        lam = n * k * f2
        return 1.0 - sp_stats.ncf.cdf(f_crit, df1, df2, lam)

    # Binary search for minimum n that achieves desired power
    lo, hi = 2, 10000
    while hi - lo > 1:
        mid = (lo + hi) // 2
        if _power_at_n(mid) >= power:
            hi = mid
        else:
            lo = mid
    return hi


def solve_power_correlation(r, alpha=0.05, power=0.80):
    """Compute sample size for Pearson correlation test."""
    _check_deps()
    z_beta = sp_stats.norm.ppf(power)
    z_alpha = sp_stats.norm.ppf(1 - alpha / 2)
    z_r = 0.5 * math.log((1 + r) / (1 - r))
    n = ((z_alpha + z_beta) / z_r) ** 2 + 3
    return int(math.ceil(n))


def solve_power_regression(f2, p, alpha=0.05, power=0.80):
    """Compute total sample size for multiple regression. f2 = Cohen's f² = R²/(1-R²)."""
    _check_deps()
    z_beta = sp_stats.norm.ppf(power)
    z_alpha = sp_stats.norm.ppf(1 - alpha / 2)
    n = ((z_alpha + z_beta) ** 2 / f2) + p + 1
    return int(math.ceil(n))


def solve_power_chisquare(w, df, alpha=0.05, power=0.80):
    """Compute total sample size for chi-square test. w = Cohen's w = Cramér's V × sqrt(min(r,c)-1)."""
    _check_deps()
    # Non-central chi-square approximation
    z_beta = sp_stats.norm.ppf(power)
    z_alpha = sp_stats.norm.ppf(1 - alpha)
    ncp = (z_alpha + z_beta) ** 2
    n = ncp / (w ** 2)
    return int(math.ceil(n))


def solve_power_equivalence(d, alpha=0.05, power=0.80):
    """TOST equivalence test sample size. d is equivalence bound in Cohen's d units."""
    _check_deps()
    # Two one-sided tests: approximate
    z_beta = sp_stats.norm.ppf(power)
    z_alpha = sp_stats.norm.ppf(1 - alpha)
    n = ((z_alpha + z_beta) ** 2) / (2 * (d ** 2))
    return int(math.ceil(n))


def solve_power_logistic(or_val, p_base, alpha=0.05, power=0.80):
    """Approximate per-group sample for logistic regression."""
    _check_deps()
    p1 = p_base * or_val / (1 - p_base + p_base * or_val)
    d = p1 - p_base
    p_bar = (p1 + p_base) / 2
    n = ((sp_stats.norm.ppf(1 - alpha/2) + sp_stats.norm.ppf(power)) ** 2 *
         (2 * p_bar * (1 - p_bar))) / (d ** 2)
    return int(math.ceil(n))


def run(args):
    if args.list_designs:
        print("Available designs:\n")
        for key, desc in DESIGNS.items():
            print(f"  {key:25s} {desc}")
        return

    if not args.design:
        print("Error: --design is required (use --list-designs to see options)", file=sys.stderr)
        sys.exit(1)

    design = args.design
    alpha = args.alpha
    power = args.power

    result = {"design": design, "alpha": alpha, "power": power, "parameters": {}}

    if design == "ttest-ind":
        if args.effect_size is not None:
            n = solve_power_ttest_ind(args.effect_size, alpha, power, args.ratio)
            result["type"] = "sample_size"
            result["n_per_group"] = n
            result["n_total"] = n * 2
            result["parameters"]["effect_size_d"] = args.effect_size
            result["note"] = f"Need {n} per group ({n * 2} total) for d = {args.effect_size}"
        elif args.n_per_group is not None:
            # Compute detectable effect size
            _check_deps()
            z_alpha = sp_stats.norm.ppf(1 - alpha / 2)
            z_beta = sp_stats.norm.ppf(power)
            d = (z_alpha + z_beta) / math.sqrt(args.n_per_group / 2)
            result["type"] = "detectable_effect"
            result["effect_size_d"] = round(d, 4)
            result["parameters"]["n_per_group"] = args.n_per_group
            result["note"] = f"With {args.n_per_group} per group, can detect d = {d:.4f}"
        else:
            print("Error: provide --effect-size or --n-per-group for ttest-ind", file=sys.stderr)
            sys.exit(1)

    elif design == "ttest-paired":
        if args.effect_size is not None:
            n = solve_power_ttest_paired(args.effect_size, alpha, power)
            result["type"] = "sample_size"
            result["n_pairs"] = n
            result["parameters"]["effect_size_dz"] = args.effect_size
            result["note"] = f"Need {n} pairs for d_z = {args.effect_size}"
        elif args.n_per_group is not None:
            _check_deps()
            z_alpha = sp_stats.norm.ppf(1 - alpha / 2)
            z_beta = sp_stats.norm.ppf(power)
            d = (z_alpha + z_beta) / math.sqrt(args.n_per_group - 2)
            result["type"] = "detectable_effect"
            result["effect_size_dz"] = round(d, 4)
            result["parameters"]["n_pairs"] = args.n_per_group
            result["note"] = f"With {args.n_per_group} pairs, can detect d_z = {d:.4f}"
        else:
            print("Error: provide --effect-size or --n-per-group for ttest-paired", file=sys.stderr)
            sys.exit(1)

    elif design == "prop":
        if args.p1 is not None and args.p2 is not None:
            n = solve_power_prop(args.p1, args.p2, alpha, power, args.ratio)
            mde = args.p2 - args.p1
            result["type"] = "sample_size"
            result["n_per_group"] = n
            result["n_total"] = n * 2
            result["parameters"]["p1"] = args.p1
            result["parameters"]["p2"] = args.p2
            result["note"] = f"Need {n} per group ({n * 2} total) to detect {mde:.1%} difference (base={args.p1:.1%})"
        else:
            print("Error: provide --p1 and --p2 for proportion test", file=sys.stderr)
            sys.exit(1)

    elif design == "anova":
        if args.k is None:
            print("Error: --k required for ANOVA", file=sys.stderr)
            sys.exit(1)
        if args.effect_size is not None:
            n = solve_power_anova(args.effect_size, args.k, alpha, power)
            eta2 = args.effect_size ** 2 / (1 + args.effect_size ** 2)
            result["type"] = "sample_size"
            result["n_per_group"] = n
            result["n_total"] = n * args.k
            result["parameters"]["k"] = args.k
            result["parameters"]["cohens_f"] = args.effect_size
            result["parameters"]["eta_squared"] = round(eta2, 4)
            result["note"] = f"Need {n} per group ({n * args.k} total) for f = {args.effect_size} (η² = {eta2:.4f})"
        elif args.n_per_group is not None:
            _check_deps()
            z_alpha = sp_stats.norm.ppf(1 - alpha / 2)
            z_beta = sp_stats.norm.ppf(power)
            f = (z_alpha + z_beta) / math.sqrt(args.n_per_group * args.k - 1)
            result["type"] = "detectable_effect"
            result["cohens_f"] = round(f, 4)
            result["parameters"]["k"] = args.k
            result["parameters"]["n_per_group"] = args.n_per_group
            result["note"] = f"With {args.n_per_group} per group ({args.k} groups), can detect f = {f:.4f}"
        else:
            print("Error: provide --effect-size or --n-per-group for ANOVA", file=sys.stderr)
            sys.exit(1)

    elif design == "correlation":
        if args.effect_size is not None:
            n = solve_power_correlation(args.effect_size, alpha, power)
            result["type"] = "sample_size"
            result["n_total"] = n
            result["parameters"]["r"] = args.effect_size
            result["note"] = f"Need N = {n} to detect r = {args.effect_size}"
        elif args.n_per_group is not None:
            _check_deps()
            z_alpha = sp_stats.norm.ppf(1 - alpha / 2)
            z_beta = sp_stats.norm.ppf(power)
            z_r_thresh = z_alpha + z_beta / math.sqrt(args.n_per_group - 3)
            r = math.tanh(z_r_thresh)
            result["type"] = "detectable_effect"
            result["r"] = round(r, 4)
            result["parameters"]["n"] = args.n_per_group
            result["note"] = f"With N = {args.n_per_group}, can detect r = {r:.4f}"
        else:
            print("Error: provide --effect-size or --n-per-group for correlation", file=sys.stderr)
            sys.exit(1)

    elif design == "regression":
        if args.predictors is None:
            print("Error: --predictors required for regression", file=sys.stderr)
            sys.exit(1)
        if args.effect_size is not None:
            f2 = args.effect_size  # Cohen's f²
            n = solve_power_regression(f2, args.predictors, alpha, power)
            r2 = f2 / (1 + f2)
            result["type"] = "sample_size"
            result["n_total"] = n
            result["parameters"]["predictors"] = args.predictors
            result["parameters"]["cohens_f2"] = f2
            result["parameters"]["r_squared"] = round(r2, 4)
            result["note"] = f"Need N = {n} for {args.predictors} predictors, f² = {f2} (R² = {r2:.4f})"
        elif args.n_per_group is not None:
            _check_deps()
            z_alpha = sp_stats.norm.ppf(1 - alpha / 2)
            z_beta = sp_stats.norm.ppf(power)
            f2 = ((z_alpha + z_beta) ** 2) / (args.n_per_group - args.predictors - 1)
            result["type"] = "detectable_effect"
            result["cohens_f2"] = round(f2, 4)
            result["parameters"]["predictors"] = args.predictors
            result["parameters"]["n"] = args.n_per_group
            result["note"] = f"With N = {args.n_per_group}, {args.predictors} predictors, can detect f² = {f2:.4f}"
        else:
            print("Error: provide --effect-size or --n-per-group for regression", file=sys.stderr)
            sys.exit(1)

    elif design == "chi-square":
        if args.df is None:
            print("Error: --df required for chi-square", file=sys.stderr)
            sys.exit(1)
        if args.effect_size is not None:
            n = solve_power_chisquare(args.effect_size, args.df, alpha, power)
            result["type"] = "sample_size"
            result["n_total"] = n
            result["parameters"]["df"] = args.df
            result["parameters"]["w"] = args.effect_size
            result["note"] = f"Need N = {n} for χ² test, df = {args.df}, w = {args.effect_size}"
        elif args.n_per_group is not None:
            _check_deps()
            z_alpha = sp_stats.norm.ppf(1 - alpha)
            z_beta = sp_stats.norm.ppf(power)
            w = (z_alpha + z_beta) / math.sqrt(args.n_per_group)
            result["type"] = "detectable_effect"
            result["w"] = round(w, 4)
            result["parameters"]["n"] = args.n_per_group
            result["parameters"]["df"] = args.df
            result["note"] = f"With N = {args.n_per_group}, df = {args.df}, can detect w = {w:.4f}"
        else:
            print("Error: provide --effect-size or --n-per-group for chi-square", file=sys.stderr)
            sys.exit(1)

    elif design == "equivalence":
        if args.effect_size is not None:
            n = solve_power_equivalence(args.effect_size, alpha, power)
            result["type"] = "sample_size"
            result["n_total"] = n if args.design == "onesample" else n
            result["n_per_group"] = n
            result["parameters"]["equivalence_bound_d"] = args.effect_size
            result["note"] = f"Need N = {n} per group for equivalence TOST, bound d = {args.effect_size}"
        else:
            print("Error: provide --effect-size for equivalence test", file=sys.stderr)
            sys.exit(1)

    elif design == "logistic":
        if args.or_val is not None and args.p_base is not None:
            n = solve_power_logistic(args.or_val, args.p_base, alpha, power)
            p1 = args.p_base * args.or_val / (1 - args.p_base + args.p_base * args.or_val)
            result["type"] = "sample_size"
            result["n_per_group"] = n
            result["n_total"] = n * 2
            result["parameters"]["or"] = args.or_val
            result["parameters"]["p_base"] = args.p_base
            result["parameters"]["p_treated"] = round(p1, 4)
            result["note"] = f"Need N = {n} per group to detect OR = {args.or_val} from base {args.p_base}"
        else:
            print("Error: provide --or-val and --p-base for logistic power", file=sys.stderr)
            sys.exit(1)

    else:
        print(f"Error: unknown design '{design}'. Use --list-designs.", file=sys.stderr)
        sys.exit(1)

    if args.engine == "r":
        print(_to_r_code(design, result))
    elif args.json:
        print(json.dumps(result, indent=2))
    else:
        print(result.get("note", ""))
        for k, v in result.items():
            if k not in ("note", "type", "parameters"):
                if isinstance(v, float):
                    print(f"  {k}: {v:.4f}")
                else:
                    print(f"  {k}: {v}")


def _to_r_code(design, result):
    lines = ["# R power analysis code", "# Run in R with: install.packages('pwr')", ""]
    lines.append("library(pwr)")
    lines.append("")
    if result.get("type") == "sample_size":
        n = result.get("n_per_group") or result.get("n_total")
        if design == "ttest-ind":
            lines.append(f"# Two-sample t-test: n = {n} per group")
            d = result.get("parameters", {}).get("effect_size_d", "?")
            lines.append(f'pwr.t.test(d = {d}, power = {result.get("power", 0.8)}, '
                         f'sig.level = {result.get("alpha", 0.05)}, type = "two.sample")')
        elif design == "ttest-paired":
            lines.append(f'pwr.t.test(d = {result.get("parameters", {}).get("effect_size_dz", "?")}, '
                         f'power = {result.get("power", 0.8)}, sig.level = {result.get("alpha", 0.05)}, '
                         f'type = "paired")')
        elif design == "prop":
            h = 2 * math.asin(math.sqrt(result.get("parameters", {}).get("p2", 0.15))) - \
                2 * math.asin(math.sqrt(result.get("parameters", {}).get("p1", 0.10)))
            lines.append(f'h = ES.h({result.get("parameters", {}).get("p1", 0.1)}, '
                         f'{result.get("parameters", {}).get("p2", 0.15)})')
            lines.append(f'pwr.2p.test(h = {h:.4f}, n = {n}, '
                         f'sig.level = {result.get("alpha", 0.05)}, power = {result.get("power", 0.8)})')
        elif design == "correlation":
            lines.append(f'pwr.r.test(r = {result.get("parameters", {}).get("r", "?")}, '
                         f'power = {result.get("power", 0.8)}, sig.level = {result.get("alpha", 0.05)})')
        elif design == "anova":
            lines.append(f'pwr.anova.test(k = {result.get("parameters", {}).get("k", "?")}, '
                         f'f = {result.get("parameters", {}).get("cohens_f", "?")}, '
                         f'power = {result.get("power", 0.8)}, sig.level = {result.get("alpha", 0.05)})')
    return "\n".join(lines)


def main():
    parser = argparse.ArgumentParser(description="Power Analysis Calculator")
    parser.add_argument("--design", choices=list(DESIGNS.keys()), help="Study design")
    parser.add_argument("--list-designs", action="store_true", help="List available designs")
    parser.add_argument("--effect-size", type=float, help="Standardized effect size")
    parser.add_argument("--n-per-group", type=int, help="Sample size per group (for computing detectable effect)")
    parser.add_argument("--alpha", type=float, default=0.05, help="Type I error rate")
    parser.add_argument("--power", type=float, default=0.80, help="Desired statistical power")
    parser.add_argument("--ratio", type=float, default=1.0, help="Control:treated ratio")
    parser.add_argument("--k", type=int, help="Number of groups (ANOVA)")
    parser.add_argument("--df", type=int, help="Degrees of freedom (chi-square)")
    parser.add_argument("--predictors", type=int, help="Number of predictors (regression)")
    parser.add_argument("--p1", type=float, help="Proportion in group 1 (proportion test)")
    parser.add_argument("--p2", type=float, help="Proportion in group 2 (proportion test)")
    parser.add_argument("--or-val", type=float, help="Odds ratio to detect (logistic)")
    parser.add_argument("--p-base", type=float, help="Baseline proportion (logistic)")
    parser.add_argument("--json", action="store_true", help="Output as JSON")
    parser.add_argument("--engine", choices=["python", "r"], default="python",
                        help="Output language (python = compute now, r = generate R code)")
    args = parser.parse_args()
    run(args)


if __name__ == "__main__":
    main()

#!/usr/bin/env bash
# test_campaign_protocol.sh — Structural validation for experimental-campaign-protocol.md

set -euo pipefail

SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
REPO_DIR="$(cd "$SCRIPT_DIR/.." && pwd)"
REF="$REPO_DIR/references/experimental-campaign-protocol.md"
PASS=0
FAIL=0

test_case() {
    local name="$1"
    shift
    echo "  TEST: $name"
    if "$@" 2>/dev/null; then
        echo "    ✓ PASS"
        PASS=$((PASS + 1))
    else
        echo "    ✗ FAIL"
        FAIL=$((FAIL + 1))
    fi
}

echo "══════════════════════════════════════════════"
echo "Experimental Campaign Protocol Test Suite"
echo "══════════════════════════════════════════════"
echo ""

# ── File existence ───────────────────────────────────────────────
echo "── File Structure ──────────────────────────"
echo ""

test_case "Reference document exists" test -f "$REF"

# ── Phase headings ───────────────────────────────────────────────
echo ""
echo "── Phase Coverage ──────────────────────────"
echo ""

for phase_num in 1 2 3 4 5 6 7 8; do
    test_case "Phase $phase_num has heading" \
        grep -q "Phase $phase_num:" "$REF"
done

test_case "All 8 phase headings present" \
    bash -c "grep -c '^## Phase' \"$REF\" | xargs test 8 -eq"

# ── Entry/exit criteria ──────────────────────────────────────────
echo ""
echo "── Structural Elements ─────────────────────"
echo ""

test_case "Each phase has 'Entry criteria'" \
    bash -c "grep -c 'Entry criteria' \"$REF\" | xargs test 8 -eq"

test_case "Each phase has 'Exit criteria'" \
    bash -c "grep -c 'Exit criteria' \"$REF\" | xargs test 8 -eq"

test_case "Each phase has 'Failure modes'" \
    bash -c "grep -c 'Failure modes' \"$REF\" | xargs test 8 -eq"

# ── Code integration ─────────────────────────────────────────────
echo ""
echo "── Code Integration ────────────────────────"
echo ""

test_case "Contains sklearn Pipeline example" grep -q "sklearn.pipeline" "$REF"
test_case "Contains PyTorch training loop" grep -q "torch.*optim.*AdamW\|DataLoader\|model.train()" "$REF"
test_case "Contains Optuna example" grep -q "optuna" "$REF"
test_case "Contains distillation code" grep -q "distillation\|KL.*div\|teacher" "$REF"
test_case "Contains pruning reference" grep -q "prune" "$REF"

# ── Cross-references ─────────────────────────────────────────────
echo ""
echo "── Cross-References ────────────────────────"
echo ""

for ref_name in "detect-compute" "pytorch-integration" "sklearn-integration" \
                "data-science-coding-workflow" "subagent-experiment-supervision" \
                "docker-experiment-isolation"; do
    test_case "References $ref_name" \
        grep -q "$ref_name" "$REF"
done

# ── Quick reference completeness ─────────────────────────────────
echo ""
echo "── Appendix ────────────────────────────────"
echo ""

test_case "Has skip-guide table" grep -q "Skip to" "$REF"
test_case "Has directory structure" grep -q "experiments/" "$REF"
test_case "Has seed-setting code" grep -q "set_seed\|random.seed.*np.*torch" "$REF"
test_case "Has experiment logging template" grep -q "experiment_id" "$REF"

# ── Summary ──────────────────────────────────────────────────────
echo ""
echo "══════════════════════════════════════════════"
echo "Results: $PASS passed, $FAIL failed"
echo "══════════════════════════════════════════════"

if [ "$FAIL" -gt 0 ]; then
    exit 1
fi

#!/usr/bin/env bash
# test_detect_compute.sh — Test suite for detect-compute.py
#
# Usage:
#   bash scripts/test_detect_compute.sh          # Full test suite
#   bash scripts/test_detect_compute.sh --local  # Local host only (no Docker)
#   bash scripts/test_detect_compute.sh --docker # Docker tests only

set -euo pipefail

SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
REPO_DIR="$(cd "$SCRIPT_DIR/.." && pwd)"
DETECT_SCRIPT="$REPO_DIR/scripts/detect-compute.py"
TMPDIR="${TMPDIR:-/tmp}/ds-test-$$"
PASS=0
FAIL=0
DOCKER_AVAILABLE=false

if command -v docker &>/dev/null; then
    DOCKER_AVAILABLE=true
fi

mkdir -p "$TMPDIR"

cleanup() {
    rm -rf "$TMPDIR" 2>/dev/null || true
}
trap cleanup EXIT

# ── Helpers ───────────────────────────────────────────────────────

run_captured() {
    # Run a command, capture its stdout and stderr to files, return its exit code
    local name="$1"
    shift
    "$@" > "$TMPDIR/$name.stdout" 2>"$TMPDIR/$name.stderr"
    return $?
}

read_stdout() {
    cat "$TMPDIR/$1.stdout" 2>/dev/null
}

read_stderr() {
    cat "$TMPDIR/$1.stderr" 2>/dev/null
}

test_case() {
    local name="$1"
    shift
    local safe_name="${name// /_}"
    safe_name="${safe_name//\//_}"
    echo "  TEST: $name"
    if run_captured "$safe_name" "$@"; then
        echo "    ✓ PASS"
        PASS=$((PASS + 1))
    else
        local ec=$?
        echo "    ✗ FAIL (exit code $ec)"
        echo "    stdout: $(head -c 500 < "$TMPDIR/$safe_name.stdout")"
        echo "    stderr: $(head -c 500 < "$TMPDIR/$safe_name.stderr")"
        FAIL=$((FAIL + 1))
    fi
}

assert_json_valid() {
    local name="$1"
    python3 -c "import json; json.load(open('$TMPDIR/$name.stdout'))" 2>/dev/null
}

echo "══════════════════════════════════════════════"
echo "detect-compute.py Test Suite"
echo "══════════════════════════════════════════════"
echo ""

# ── Phase 1: Local host tests ────────────────────────────────────
echo "── Phase 1: Local Host ──────────────────────"
echo ""

test_case "Script runs without error" python3 "$DETECT_SCRIPT"
assert_json_valid "Script_runs_without_error" && echo "    ✓ Default output is JSON" || true

test_case "Python version detected" python3 -c "
import json
d = json.loads(open('$TMPDIR/Script_runs_without_error.stdout').read())
assert 'python_version' in d, 'Missing python_version'
print(d['python_version'])
"

test_case "--json produces valid JSON" python3 "$DETECT_SCRIPT" --json
assert_json_valid "--json_produces_valid_JSON" && echo "    ✓ --json output is valid JSON" || true

test_case "--json output has all required keys" python3 -c "
import json
d = json.loads(open('$TMPDIR/--json_produces_valid_JSON.stdout').read())
required = ['python_version', 'platform', 'has_cuda', 'torch', 'sklearn', 'recommendations']
for k in required:
    assert k in d, f'Missing key: {k}'
print(f'All {len(required)} required keys present')
"

test_case "--minimal returns only recommendations" python3 "$DETECT_SCRIPT" --minimal
assert_json_valid "--minimal_returns_only_recommendations" && echo "    ✓ --minimal output is valid JSON" || true

test_case "--minimal output has recommendation keys" python3 -c "
import json
d = json.loads(open('$TMPDIR/--minimal_returns_only_recommendations.stdout').read())
assert 'model_size_tier' in d, 'Missing model_size_tier'
assert 'feasible_techniques' in d, 'Missing feasible_techniques'
assert 'batch_size_guide' in d, 'Missing batch_size_guide'
print(f'Tier: {d[\"model_size_tier\"]}, techniques: {d[\"feasible_techniques\"]}')
"

test_case "--list-gpus runs without error" python3 "$DETECT_SCRIPT" --list-gpus

# ── Phase 2: Docker tests (no GPU, no torch) ────────────────────
echo ""
echo "── Phase 2: Docker (no GPU, no torch) ───────"
echo ""

if [ "$DOCKER_AVAILABLE" = true ] && [ "${1:-}" != "--local" ]; then
    DOCKER_TAG="ds-detect-test-$$"

    cat > "$TMPDIR/Dockerfile.nogpu" << 'DOCKERFILE'
FROM python:3.12-slim
RUN pip install --quiet --no-cache-dir scikit-learn numpy psutil
COPY scripts/detect-compute.py /scripts/detect-compute.py
WORKDIR /
DOCKERFILE

    mkdir -p "$TMPDIR/context/scripts"
    cp "$DETECT_SCRIPT" "$TMPDIR/context/scripts/detect-compute.py"
    cp "$TMPDIR/Dockerfile.nogpu" "$TMPDIR/context/Dockerfile"

    echo "  Building Docker test image (no GPU, no torch)..."
    docker build -t "$DOCKER_TAG" -f "$TMPDIR/context/Dockerfile" "$TMPDIR/context" >/dev/null 2>&1

    test_case "Docker: script runs" \
        docker run --rm "$DOCKER_TAG" python3 /scripts/detect-compute.py --json

    test_case "Docker: has_cuda is false" \
        docker run --rm "$DOCKER_TAG" python3 -c "
import json, subprocess, sys
r = subprocess.run([sys.executable, '/scripts/detect-compute.py', '--json'], capture_output=True, text=True)
d = json.loads(r.stdout)
assert d['has_cuda'] == False, 'GPU should not be detected in plain container'
assert d['nvidia'] == {}, 'nvidia should be empty'
print('OK: has_cuda=false, nvidia={}')
"

    test_case "Docker: torch is not available" \
        docker run --rm "$DOCKER_TAG" python3 -c "
import json, subprocess, sys
r = subprocess.run([sys.executable, '/scripts/detect-compute.py', '--json'], capture_output=True, text=True)
d = json.loads(r.stdout)
assert d['torch']['available'] == False, 'torch should not be available'
print('OK: torch not available')
"

    test_case "Docker: sklearn is available" \
        docker run --rm "$DOCKER_TAG" python3 -c "
import json, subprocess, sys
r = subprocess.run([sys.executable, '/scripts/detect-compute.py', '--json'], capture_output=True, text=True)
d = json.loads(r.stdout)
assert d['sklearn']['available'] == True, 'sklearn should be available'
print(f'sklearn {d[\"sklearn\"][\"version\"]}')
"

    test_case "Docker: recommendations reflect CPU-only" \
        docker run --rm "$DOCKER_TAG" python3 -c "
import json, subprocess, sys
r = subprocess.run([sys.executable, '/scripts/detect-compute.py', '--json'], capture_output=True, text=True)
d = json.loads(r.stdout)
rec = d['recommendations']
assert rec['model_size_tier'] == 'cpu_only', f'Expected cpu_only, got {rec[\"model_size_tier\"]}'
assert 'sklearn' in str(rec['feasible_techniques']), 'Should recommend sklearn'
print(f'Tier: {rec[\"model_size_tier\"]}, techniques: {rec[\"feasible_techniques\"]}')
"

    docker image rm "$DOCKER_TAG" >/dev/null 2>&1 || true
else
    echo "  (skipping — Docker not available or --local flag)"
fi

# ── Summary ──────────────────────────────────────────────────────
echo ""
echo "══════════════════════════════════════════════"
echo "Results: $PASS passed, $FAIL failed"
echo "══════════════════════════════════════════════"

if [ "$FAIL" -gt 0 ]; then
    exit 1
fi

#!/usr/bin/env bash
# test_references_completeness.sh — Validate the researched code integration references
# 
# Checks:
# - All three reference files exist
# - Each covers the required topic areas
# - Each cross-references the parent skill
# - Source URLs are documented

set -euo pipefail

SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
REPO_DIR="$(cd "$SCRIPT_DIR/.." && pwd)"
PASS=0
FAIL=0

test_case() {
    local name="$1"
    shift
    echo "  TEST: $name"
    if "$@" 2>/dev/null; then
        echo "    ✓ PASS"
        PASS=$((PASS + 1))
    else
        echo "    ✗ FAIL"
        FAIL=$((FAIL + 1))
    fi
}

echo "══════════════════════════════════════════════"
echo "Code Integration References Test Suite"
echo "══════════════════════════════════════════════"
echo ""

# ── File existence ───────────────────────────────────────────────
echo "── File Existence ──────────────────────────"
echo ""

test_case "pytorch-integration.md exists" \
    test -f "$REPO_DIR/references/pytorch-integration.md"
test_case "sklearn-integration.md exists" \
    test -f "$REPO_DIR/references/sklearn-integration.md"
test_case "data-science-coding-workflow.md exists" \
    test -f "$REPO_DIR/references/data-science-coding-workflow.md"

# ── PyTorch reference coverage ─────────────────────────────────
echo ""
echo "── PyTorch Reference Coverage ─────────────"
echo ""

PT="$REPO_DIR/references/pytorch-integration.md"

for topic in "Device Management" "Training Loop" "DataLoader" "AMP" "Mixed Precision" \
             "torch.compile" "Transfer Learning" "LoRA" "Knowledge Distillation" \
             "Model Pruning" "DDP" "Distributed Data" "Reproducibility"; do
    test_case "Covers: $topic" grep -qi "$topic" "$PT"
done

test_case "Has device pattern (cuda/mps/cpu)" \
    grep -q "torch.device" "$PT"
test_case "Has gradient clipping" \
    grep -q "clip_grad_norm" "$PT"
test_case "Has model saving/loading pattern" \
    grep -q "torch.save\|model.state_dict" "$PT"
test_case "Has loss function table" \
    grep -q "CrossEntropyLoss\|BCEWithLogitsLoss\|MSELoss" "$PT"
test_case "Has learning rate schedulers" \
    grep -q "ReduceLROnPlateau\|CosineAnnealingLR\|OneCycleLR" "$PT"
test_case "Has debugging section" \
    grep -q "Debugging\|Common Failure" "$PT"
test_case "Source validated date present" \
    grep -q "Last reviewed\|Source validated" "$PT"
test_case "References pytorch.org" \
    grep -q "pytorch.org" "$PT"
test_case "References experimental-campaign-protocol" \
    grep -q "experimental-campaign-protocol" "$PT"

# ── sklearn reference coverage ─────────────────────────────────
echo ""
echo "── Scikit-Learn Reference Coverage ────────"
echo ""

SK="$REPO_DIR/references/sklearn-integration.md"

for topic in "Pipeline" "ColumnTransformer" "Preprocessing" "Model Selection" \
             "Cross-Validation" "GridSearchCV" "Ensemble" "Calibration" \
             "Imbalanced" "Feature Selection" "PCA" "Custom Estimator" \
             "Persistence" "Reproducibility"; do
    test_case "Covers: $topic" grep -qi "$topic" "$SK"
done

test_case "Has ColumnTransformer example" \
    grep -q "ColumnTransformer" "$SK"
test_case "Has OneHotEncoder + StandardScaler" \
    grep -q "OneHotEncoder\|StandardScaler" "$SK"
test_case "Has imputation (SimpleImputer/IterativeImputer)" \
    grep -q "SimpleImputer\|IterativeImputer" "$SK"
test_case "Has HalvingGridSearchCV" \
    grep -q "HalvingGridSearchCV" "$SK"
test_case "Has Random Forest example" \
    grep -q "RandomForest" "$SK"
test_case "Has XGBoost/LightGBM integration" \
    grep -q "XGBClassifier\|LGBMClassifier" "$SK"
test_case "Source validated date present" \
    grep -q "Last reviewed\|Source validated" "$SK"
test_case "References scikit-learn.org" \
    grep -q "scikit-learn.org" "$SK"
test_case "References experimental-campaign-protocol" \
    grep -q "experimental-campaign-protocol" "$SK"

# ── DS Coding Workflow coverage ────────────────────────────────
echo ""
echo "── DS Coding Workflow Coverage ────────────"
echo ""

WF="$REPO_DIR/references/data-science-coding-workflow.md"

for topic in "Project Directory" "Configuration" "Experiment Logging" \
             "MLflow" "Result Serialization" "Reproducibility" \
             "Data Versioning" "Unit Testing" "Docker" "Seed"; do
    test_case "Covers: $topic" grep -qi "$topic" "$WF"
done

test_case "Has directory structure layout" \
    grep -q "data/raw/\|data/processed/" "$WF"
test_case "Has MLflow example" \
    grep -q "mlflow" "$WF"
test_case "Has DVC reference" \
    grep -q "dvc\|DVC" "$WF"
test_case "Has JSON experiment log pattern" \
    grep -q "experiment_log\.json\|json" "$WF"
test_case "Has reproducibility section" \
    grep -q "Reproducibility\|random_state\|set_all_seeds" "$WF"
test_case "Has pitfalls table" \
    grep -q "Pitfall\|pitfall" "$WF"
test_case "Source validated date present" \
    grep -q "Last reviewed\|Source validated" "$WF"
test_case "References experimental-campaign-protocol" \
    grep -q "experimental-campaign-protocol" "$WF"

# ── Summary ──────────────────────────────────────────────────────
echo ""
echo "══════════════════════════════════════════════"
echo "Results: $PASS passed, $FAIL failed"
echo "══════════════════════════════════════════════"

if [ "$FAIL" -gt 0 ]; then
    exit 1
fi

#!/usr/bin/env bash
# test_supervision_protocol.sh — Validate subagent supervision and Docker isolation references

set -euo pipefail

SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
REPO_DIR="$(cd "$SCRIPT_DIR/.." && pwd)"
PASS=0
FAIL=0
DOCKER_AVAILABLE=false

if command -v docker &>/dev/null; then
    DOCKER_AVAILABLE=true
fi

test_case() {
    local name="$1"
    shift
    echo "  TEST: $name"
    if "$@" 2>/dev/null; then
        echo "    ✓ PASS"
        PASS=$((PASS + 1))
    else
        echo "    ✗ FAIL"
        FAIL=$((FAIL + 1))
    fi
}

echo "══════════════════════════════════════════════"
echo "Supervision & Isolation Test Suite"
echo "══════════════════════════════════════════════"
echo ""

# ── File existence ───────────────────────────────────────────────
echo "── File Existence ──────────────────────────"
echo ""

test_case "subagent-experiment-supervision.md exists" \
    test -f "$REPO_DIR/references/subagent-experiment-supervision.md"
test_case "docker-experiment-isolation.md exists" \
    test -f "$REPO_DIR/references/docker-experiment-isolation.md"
test_case "scripts/Dockerfile exists" \
    test -f "$REPO_DIR/scripts/Dockerfile"

# ── Subagent supervision content coverage ───────────────────────
echo ""
echo "── Supervision Reference Coverage ──────────"
echo ""

SUP="$REPO_DIR/references/subagent-experiment-supervision.md"

test_case "Describes architecture (orchestrator/worker/supervisor)" \
    grep -qi "orchestrator\|supervisor.*worker\|architecture" "$SUP"
test_case "Has failure catalog table" \
    grep -q "CUDA OOM\|NaN loss\|ImportError" "$SUP"
test_case "Has at least 7 failure entries" \
    bash -c "grep -c '| \*\*' \"$SUP\" || grep -c '^| \*\*' \"$SUP\" || grep -c 'OOM\|NaN\|ImportError\|Disk full\|cuDNN' \"$SUP\" | xargs test 7 -le"
test_case "Each failure has a fix" \
    grep -q "batch_size\|pip install\|gradient.*clip\|fallback.*CPU\|reduce.*scope\|clean" "$SUP"
test_case "Has escalation path" \
    grep -qi "escalat\|telegram\|notify\|alert" "$SUP"
test_case "Has fix implementation code" \
    grep -q "def fix_\|def apply_fix\|def diagnose" "$SUP"
test_case "Has Python supervision loop example" \
    grep -q "subprocess\.Popen\|supervise_experiment\|FAILURE_PATTERNS" "$SUP"
test_case "Has harness-specific notes" \
    grep -qi "Hermes\|delegate_task\|OpenCode\|subagent" "$SUP"
test_case "Has limitations section" \
    grep -qi "limitation\|Limitation" "$SUP"
test_case "References experimental-campaign-protocol" \
    grep -q "experimental-campaign-protocol" "$SUP"

# ── Docker isolation content coverage ───────────────────────────
echo ""
echo "── Docker Isolation Coverage ──────────────"
echo ""

DKR="$REPO_DIR/references/docker-experiment-isolation.md"

test_case "Has resource limit guidance" \
    grep -qi "memory.*limit\|--memory\|--cpus\|--gpus" "$DKR"
test_case "Has full docker run example" \
    grep -q "docker run" "$DKR"
test_case "Has log collection pattern" \
    grep -q "docker logs\|logging.*stdout" "$DKR"
test_case "Has cleanup pattern" \
    grep -q "prune\|clean\|cleanup\|docker rm" "$DKR"
test_case "Has multi-container sweep example" \
    grep -q "for trial\|docker-compose\|compose" "$DKR"
test_case "Has fallback when Docker unavailable" \
    grep -qi "conda\|venv\|When Docker is not\|fallback" "$DKR"
test_case "References subagent-experiment-supervision" \
    grep -q "subagent-experiment-supervision" "$DKR"
test_case "References experimental-campaign-protocol" \
    grep -q "experimental-campaign-protocol" "$DKR"

# ── Docker build test ───────────────────────────────────────────
echo ""
echo "── Docker Build ───────────────────────────"
echo ""

if [ "$DOCKER_AVAILABLE" = true ]; then
    echo "  Building test image (this may take a moment)..."
    if docker build -q -t ds-supervision-test -f "$REPO_DIR/scripts/Dockerfile" "$REPO_DIR" >/dev/null 2>&1; then
        PASS=$((PASS + 1))
        echo "    ✓ Docker image builds successfully"
        docker image rm ds-supervision-test >/dev/null 2>&1 || true
    else
        FAIL=$((FAIL + 1))
        echo "    ✗ Docker image build failed"
    fi
else
    echo "  ⚠ Docker not available — skipping build test"
fi

# ── Summary ──────────────────────────────────────────────────────
echo ""
echo "══════════════════════════════════════════════"
echo "Results: $PASS passed, $FAIL failed"
echo "══════════════════════════════════════════════"

if [ "$FAIL" -gt 0 ]; then
    exit 1
fi

# Interpretability report

- Decision, stakes, audience, and explanation target:
- Model/data/preprocessing versions:
- Scope: global / local / representation / fairness / contrastive:
- Reference population or local background:
- Feature correlations, proxies, perturbation validity, and shift:
- Method(s), assumptions, and implementation versions:
- Stability seeds/backgrounds/perturbations and sanity checks:
- Method disagreement and unresolved findings:

| Slice | N | Result | Uncertainty/limitations | Action |
|---|---:|---|---|---|
| Overall | | | | |
| High-impact/protected slice | | | | |

- Predictive interpretation:
- Causal claim (only if identified design supports it):
- Sensitive fields/access/retention controls:
- Reviewer and decision: