import numpy as np
import pandas as pd
import hashlib
import json
import datetime
from pathlib import Path

# 1. Set random seeds for reproducibility
np.random.seed(42)

# 2. Create a reproducible analysis pipeline
class ReproducibleAnalysis:
    def __init__(self, seed=42):
        self.seed = seed
        self.artifacts = {}
        self.log = []

    def log_step(self, step_name, description):
        entry = {
            "step": step_name,
            "description": description,
            "timestamp": datetime.datetime.now().isoformat()
        }
        self.log.append(entry)
        return entry

    def compute_hash(self, data, algorithm='sha256'):
        """Compute hash of data for integrity verification."""
        if isinstance(data, np.ndarray):
            data_bytes = data.tobytes()
        elif isinstance(data, pd.DataFrame):
            data_bytes = pd.util.hash_pandas_object(data).values.tobytes()
        else:
            data_bytes = str(data).encode()
        return hashlib.new(algorithm, data_bytes).hexdigest()

    def save_artifact(self, name, data, directory='artifacts'):
        """Save and hash an analysis artifact."""
        Path(directory).mkdir(exist_ok=True)
        filepath = Path(directory) / f"{name}.csv"

        if isinstance(data, pd.DataFrame):
            data.to_csv(filepath, index=False)
        else:
            pd.DataFrame(data).to_csv(filepath, index=False)

        file_hash = self.compute_hash(data)
        self.artifacts[name] = {
            "path": str(filepath),
            "hash": file_hash,
            "rows": len(data) if hasattr(data, '__len__') else None
        }
        self.log_step(f"save_{name}", f"Saved {filepath} (hash: {file_hash[:16]}...)")
        return filepath

    def generate_report(self):
        """Generate a reproducibility report."""
        report = {
            "seed": self.seed,
            "numpy_version": np.__version__,
            "pandas_version": pd.__version__,
            "artifacts": self.artifacts,
            "execution_log": self.log,
            "generated_at": datetime.datetime.now().isoformat()
        }
        return report

# 3. Example analysis with full traceability
analysis = ReproducibleAnalysis(seed=42)
analysis.log_step("init", "Initialize reproducible analysis with seed=42")

# Simulate data
n = 500
data = pd.DataFrame({
    'x': np.random.normal(0, 1, n),
    'group': np.random.choice(['A', 'B', 'C'], n)
})
data['y'] = 2.5 * data['x'] + np.random.normal(0, 1, n)

# Save raw data
analysis.save_artifact("raw_data", data)

# Analysis step
analysis.log_step("regression", "Fit OLS: y ~ x")
from numpy.linalg import lstsq
X = np.column_stack([np.ones(n), data['x'].values])
beta, residuals, rank, sv = lstsq(X, data['y'].values, rcond=None)

# Save results
results = pd.DataFrame({'coefficient': ['intercept', 'slope'], 'estimate': beta})
analysis.save_artifact("regression_results", results)

# 4. Generate reproducibility report
report = analysis.generate_report()
print(json.dumps(report, indent=2))

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

def p_curve(p_values):
    """
    Analyze the distribution of significant p-values.
    A right-skewed p-curve suggests evidential value.
    A flat or left-skewed curve suggests p-hacking.
    """
    sig_pvals = p_values[p_values < 0.05]

    # Binomial test: under H1, p-values should be right-skewed
    bins = [0, 0.01, 0.02, 0.03, 0.04, 0.05]
    observed_counts, _ = np.histogram(sig_pvals, bins=bins)

    # Expected under uniform (no effect)
    expected = np.full_like(observed_counts, len(sig_pvals) / 5, dtype=float)

    # Chi-square test for uniformity
    chi2, p_uniform = stats.chisquare(observed_counts, f_exp=expected)

    # Skewness of significant p-values
    if len(sig_pvals) > 2:
        skewness = stats.skew(sig_pvals)
    else:
        skewness = np.nan

    return {
        'n_significant': len(sig_pvals),
        'n_total': len(p_values),
        'chi2_uniformity': chi2,
        'p_uniformity': p_uniform,
        'skewness': skewness,
        'counts': observed_counts,
        'bins': bins
    }

# Simulate p-values with and without p-hacking
np.random.seed(42)
n_studies = 200

# Scenario 1: Genuine effects (right-skewed p-curve)
genuine_pvals = np.concatenate([
    np.random.beta(0.5, 20, 80),   # Strong effects
    np.random.uniform(0, 1, 120)    # Null effects
])

# Scenario 2: p-hacked (uniform or left-skewed)
p_hacked_pvals = np.random.uniform(0.01, 0.05, 100)  # Only p < .05 reported

# Analyze both
result_genuine = p_curve(genuine_pvals)
result_hacked = p_curve(p_hacked_pvals)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

axes[0].bar(range(5), result_genuine['counts'], color='steelblue', edgecolor='black')
axes[0].set_xticks(range(5))
axes[0].set_xticklabels(['0-0.01', '0.01-0.02', '0.02-0.03', '0.03-0.04', '0.04-0.05'])
axes[0].set_title(f'Genuine Effects\nSkewness: {result_genuine["skewness"]:.2f}')
axes[0].set_ylabel('Frequency')

axes[1].bar(range(5), result_hacked['counts'], color='coral', edgecolor='black')
axes[1].set_xticks(range(5))
axes[1].set_xticklabels(['0-0.01', '0.01-0.02', '0.02-0.03', '0.03-0.04', '0.04-0.05'])
axes[1].set_title(f'p-Hacked Data\nSkewness: {result_hacked["skewness"]:.2f}')

plt.tight_layout()
plt.savefig('p_curve_analysis.png', dpi=150)
plt.show()