statistical_power.py

#!/usr/bin/env python3
"""
GENUS-2 EXPERIMENT: STATISTICAL POWER ANALYSIS (CORRECTED)

CRITICAL: Balance noise is ABSOLUTE (μg), not RELATIVE (ppm).

A 0.01 mg balance has ~10 μg noise regardless of sample mass.
This means:
  - 100 mg sample: 10 μg noise = 100 ppm effective noise
  - 1000 mg sample: 10 μg noise = 10 ppm effective noise

The signal scales with mass:
  - Genus-1: 86 ppm × mass
  - Genus-2: 365 ppm × mass

Therefore: LARGER SAMPLES = BETTER SNR

Author: Timothy McGirl
Email: grapheneaffiliates@gmail.com
"""

import math
import random

# =============================================================================
# CONSTANTS
# =============================================================================

PHI = (1 + math.sqrt(5)) / 2  # Golden ratio = 1.618...
PHI_CUBED = PHI ** 3          # 4.2360679...

# Predicted fractional effects (ppm) if framework is correct
ZETA_G1_PPM = 86.0    # Genus-1 (torus) mass anomaly
ZETA_G2_PPM = 365.0   # Genus-2 (figure-8) mass anomaly
TRUE_RATIO = PHI_CUBED  # 4.236...

# =============================================================================
# REALISTIC BALANCE NOISE MODEL
# =============================================================================
#
# Key insight: Balance noise is in ABSOLUTE micrograms, not relative ppm.
#
# Typical analytical balance (0.01 mg resolution):
#   - Quantization: 10 μg
#   - Repeatability: ~10-20 μg (1σ)
#   - Total noise per reading: ~10-30 μg
#
# Microbalance (0.001 mg resolution):
#   - Quantization: 1 μg
#   - Repeatability: ~1-5 μg (1σ)
#   - Total noise per reading: ~2-5 μg
#
# =============================================================================

def simulate_experiment(sample_mass_mg, balance_noise_ug, n_readings, n_cycles,
                        true_g1_ppm=ZETA_G1_PPM, true_g2_ppm=ZETA_G2_PPM):
    """
    Simulate a complete genus-2 ratio experiment.
    
    Args:
        sample_mass_mg: Sample mass in milligrams
        balance_noise_ug: 1σ balance noise in micrograms (ABSOLUTE)
        n_readings: Number of readings per measurement (for averaging)
        n_cycles: Number of warm/cold cycles
        true_g1_ppm: True genus-1 effect in ppm
        true_g2_ppm: True genus-2 effect in ppm
    
    Returns:
        mean_R, std_R, list of R values per cycle
    """
    # Convert ppm to absolute mass change in μg
    # Δm = (ppm / 1e6) × mass_mg × 1000 μg/mg = ppm × mass_mg / 1000
    signal_g1_ug = true_g1_ppm * sample_mass_mg / 1000  # μg
    signal_g2_ug = true_g2_ppm * sample_mass_mg / 1000  # μg
    
    # Effective noise after averaging n_readings
    effective_noise_ug = balance_noise_ug / math.sqrt(n_readings)
    
    ratios = []
    
    for _ in range(n_cycles):
        # Simulate measurements (in μg)
        measured_g1_ug = random.gauss(signal_g1_ug, effective_noise_ug)
        measured_g2_ug = random.gauss(signal_g2_ug, effective_noise_ug)
        
        # Convert back to ppm for ratio calculation
        measured_g1_ppm = measured_g1_ug * 1000 / sample_mass_mg
        measured_g2_ppm = measured_g2_ug * 1000 / sample_mass_mg
        
        # Compute ratio (skip if g1 is zero or negative)
        if measured_g1_ppm > 0:
            R = measured_g2_ppm / measured_g1_ppm
            if R > 0:
                ratios.append(R)
    
    if len(ratios) < 2:
        return None, None, []
    
    mean_R = sum(ratios) / len(ratios)
    var_R = sum((r - mean_R)**2 for r in ratios) / (len(ratios) - 1)
    std_R = math.sqrt(var_R)
    
    return mean_R, std_R, ratios


def run_monte_carlo(n_experiments, sample_mass_mg, balance_noise_ug, 
                    n_readings, n_cycles, scenario="correct"):
    """
    Run Monte Carlo simulation for many experiments.
    
    scenario: "correct" (true R = φ³) or "wrong" (true R = 1)
    """
    results = []
    
    for _ in range(n_experiments):
        if scenario == "correct":
            mean_R, std_R, _ = simulate_experiment(
                sample_mass_mg, balance_noise_ug, n_readings, n_cycles,
                ZETA_G1_PPM, ZETA_G2_PPM
            )
        else:  # scenario == "wrong"
            # No topological scaling: both have same effect
            mean_R, std_R, _ = simulate_experiment(
                sample_mass_mg, balance_noise_ug, n_readings, n_cycles,
                ZETA_G1_PPM, ZETA_G1_PPM  # Both genus-1
            )
        
        if mean_R is not None:
            results.append({
                'mean_R': mean_R,
                'std_R': std_R,
                'passes': mean_R > 3.5 and mean_R < 5.0
            })
    
    return results


def main():
    random.seed(42)
    
    print("=" * 72)
    print("GENUS-2 EXPERIMENT: STATISTICAL POWER ANALYSIS (CORRECTED)")
    print("=" * 72)
    print()
    print("CRITICAL INSIGHT: Balance noise is ABSOLUTE (μg), not RELATIVE (ppm)")
    print()
    print("Signal scales with mass:")
    print(f"  Genus-1: {ZETA_G1_PPM} ppm × mass")
    print(f"  Genus-2: {ZETA_G2_PPM} ppm × mass")
    print(f"  Predicted ratio: φ³ = {PHI_CUBED:.4f}")
    print()
    
    # Standard analytical balance parameters
    balance_noise_ug = 10.0  # μg per reading
    n_readings = 10
    n_cycles = 10
    n_experiments = 10000
    
    print("BALANCE PARAMETERS")
    print("-" * 40)
    print(f"  Balance noise (1σ):    {balance_noise_ug} μg per reading")
    print(f"  Readings/measurement:  {n_readings}")
    print(f"  Effective noise:       {balance_noise_ug/math.sqrt(n_readings):.1f} μg")
    print(f"  Measurement cycles:    {n_cycles}")
    print()
    
    print("=" * 72)
    print("SCENARIO COMPARISON: 100 mg vs 1000 mg SAMPLES")
    print("=" * 72)
    print()
    
    for mass_mg in [100, 500, 1000, 2000]:
        signal_g1 = ZETA_G1_PPM * mass_mg / 1000
        signal_g2 = ZETA_G2_PPM * mass_mg / 1000
        eff_noise = balance_noise_ug / math.sqrt(n_readings)
        
        snr_g1 = signal_g1 / eff_noise
        snr_g2 = signal_g2 / eff_noise
        
        print(f"SAMPLE MASS: {mass_mg} mg")
        print("-" * 40)
        print(f"  Genus-1 signal: {signal_g1:.1f} μg  (SNR = {snr_g1:.1f})")
        print(f"  Genus-2 signal: {signal_g2:.1f} μg  (SNR = {snr_g2:.1f})")
        print(f"  Effective noise: {eff_noise:.1f} μg")
        print()
        
        # Run Monte Carlo if framework is correct
        results_correct = run_monte_carlo(
            n_experiments, mass_mg, balance_noise_ug, n_readings, n_cycles, "correct"
        )
        
        # Run Monte Carlo if framework is wrong
        results_wrong = run_monte_carlo(
            n_experiments, mass_mg, balance_noise_ug, n_readings, n_cycles, "wrong"
        )
        
        if results_correct:
            mean_Rs = [r['mean_R'] for r in results_correct]
            std_Rs = [r['std_R'] for r in results_correct]
            pass_rate = sum(1 for r in results_correct if r['passes']) / len(results_correct) * 100
            
            avg_R = sum(mean_Rs) / len(mean_Rs)
            avg_std = sum(std_Rs) / len(std_Rs)
            
            print(f"  IF FRAMEWORK CORRECT:")
            print(f"    Mean R: {avg_R:.2f} ± {avg_std:.2f}")
            print(f"    Pass rate (3.5 < R < 5.0): {pass_rate:.1f}%")
        
        if results_wrong:
            mean_Rs = [r['mean_R'] for r in results_wrong]
            false_pos = sum(1 for r in results_wrong if r['passes']) / len(results_wrong) * 100
            avg_R = sum(mean_Rs) / len(mean_Rs)
            
            print(f"  IF FRAMEWORK WRONG:")
            print(f"    Mean R: {avg_R:.2f}")
            print(f"    False positive rate: {false_pos:.2f}%")
        
        # Verdict
        if mass_mg == 100:
            print(f"  VERDICT: ❌ FAILS - Signal buried in noise")
        elif mass_mg >= 1000:
            print(f"  VERDICT: ✓ WORKS - Clear signal separation")
        else:
            print(f"  VERDICT: ⚠ MARGINAL - May work with careful averaging")
        
        print()
    
    print("=" * 72)
    print("RECOMMENDED CONFIGURATION")
    print("=" * 72)
    print()
    print("  Sample mass:     1000 mg (1 gram) minimum")
    print("  Balance:         0.01 mg analytical balance (~10 μg noise)")
    print("  Readings:        10 per measurement")
    print("  Cycles:          10 warm/cold")
    print()
    print("  Expected results if framework correct:")
    
    results = run_monte_carlo(10000, 1000, 10.0, 10, 10, "correct")
    mean_Rs = [r['mean_R'] for r in results]
    std_Rs = [r['std_R'] for r in results]
    avg_R = sum(mean_Rs) / len(mean_Rs)
    avg_std = sum(std_Rs) / len(std_Rs)
    pass_rate = sum(1 for r in results if r['passes']) / len(results) * 100
    
    print(f"    R = {avg_R:.2f} ± {avg_std:.2f}")
    print(f"    Pass rate: {pass_rate:.1f}%")
    print(f"    σ from null (R=1): {(avg_R - 1) / avg_std:.1f}σ")
    print(f"    σ from prediction (R=φ³): {abs(avg_R - PHI_CUBED) / avg_std:.1f}σ")
    print()
    
    print("=" * 72)
    print("BOTTOM LINE")
    print("=" * 72)
    print()
    print("  ❌ 100 mg samples: Signal buried in noise. WILL NOT WORK.")
    print("  ✓  1000 mg samples: Clear separation. WILL WORK.")
    print()
    print("  Use at least 1 gram of niobium per sample.")
    print("  A standard university analytical balance is sufficient.")
    print()
    print("=" * 72)


if __name__ == "__main__":
    main()