%matplotlib inline
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

# SciDEX dark theme for all plots
plt.rcParams.update({
    'figure.facecolor': '#0a0a14',
    'axes.facecolor': '#151525',
    'axes.edgecolor': '#333',
    'axes.labelcolor': '#e0e0e0',
    'text.color': '#e0e0e0',
    'xtick.color': '#888',
    'ytick.color': '#888',
    'legend.facecolor': '#151525',
    'legend.edgecolor': '#333',
    'figure.dpi': 120,
    'savefig.dpi': 120,
})
print('Environment ready: numpy, matplotlib, scipy')

# Score Heatmap — Hypothesis Dimensions
scores = {"composite": 0.767541, "mech": 0.85, "evid": 0.82, "novel": 0.78, "feas": 0.88, "impact": 0.8, "drug": 0.75, "safety": 0.9, "comp": 0.7, "data": 0.85, "reprod": 0.82}

dim_labels = ['Mechanistic', 'Evidence', 'Novelty', 'Feasibility', 'Impact',
              'Druggability', 'Safety', 'Competition', 'Data Avail.', 'Reproducibility']

matrix = np.array([[scores[k] for k in ['mech', 'evid', 'novel', 'feas', 'impact',
                                          'drug', 'safety', 'comp', 'data', 'reprod']]])

fig, ax = plt.subplots(figsize=(14, 2.5))
im = ax.imshow(matrix, cmap='RdYlGn', aspect='auto', vmin=0, vmax=1)

ax.set_xticks(range(len(dim_labels)))
ax.set_xticklabels(dim_labels, rotation=45, ha='right', fontsize=10)
ax.set_yticks([0])
ax.set_yticklabels(['Score'], fontsize=10)

for j in range(len(dim_labels)):
    val = matrix[0, j]
    color = '#000' if val > 0.5 else '#fff'
    ax.text(j, 0, f'{val:.2f}', ha='center', va='center', fontsize=11,
            fontweight='bold', color=color)

cbar = plt.colorbar(im, ax=ax, shrink=0.8, pad=0.02)
cbar.set_label('Score (0-1)', fontsize=10, color='#e0e0e0')
cbar.ax.yaxis.set_tick_params(color='#888')

ax.set_title('Hypothesis Score Profile — 10 Dimensions', fontsize=14,
             color='#4fc3f7', fontweight='bold')
plt.tight_layout()
plt.show()

# Multi-Dimensional Score Radar Chart
scores = {"composite": 0.767541, "mech": 0.85, "evid": 0.82, "novel": 0.78, "feas": 0.88, "impact": 0.8, "drug": 0.75, "safety": 0.9, "comp": 0.7, "data": 0.85, "reprod": 0.82}
title = "Gamma entrainment therapy to restore hippocampal-cortical synchrony"

dimensions = ['Mechanistic', 'Evidence', 'Novelty', 'Feasibility', 'Impact',
              'Druggability', 'Safety', 'Competition', 'Data Avail.', 'Reproducibility']
dim_keys = ['mech', 'evid', 'novel', 'feas', 'impact', 'drug', 'safety', 'comp', 'data', 'reprod']

fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(polar=True))
angles = np.linspace(0, 2 * np.pi, len(dimensions), endpoint=False).tolist()
angles += angles[:1]

values = [scores[k] for k in dim_keys]
values += values[:1]

ax.plot(angles, values, 'o-', linewidth=2.5, color='#4fc3f7', alpha=0.9, markersize=8)
ax.fill(angles, values, alpha=0.2, color='#4fc3f7')

# Add threshold ring
threshold = [0.7] * (len(dimensions) + 1)
ax.plot(angles, threshold, '--', linewidth=1, color='#81c784', alpha=0.5, label='Strong (0.7)')

ax.set_xticks(angles[:-1])
ax.set_xticklabels(dimensions, fontsize=9)
ax.set_ylim(0, 1)
ax.set_title(f'Score Radar: {title[:50]}', fontsize=14, color='#4fc3f7',
             fontweight='bold', pad=20)
ax.legend(loc='upper right', bbox_to_anchor=(1.3, 1.1), fontsize=9,
          facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')
plt.tight_layout()
plt.show()

print(f"\nComposite Score: {scores.get('composite', 0):.3f}")
print(f"Strongest dimension: {dimensions[np.argmax([scores[k] for k in dim_keys])]}")
print(f"Weakest dimension: {dimensions[np.argmin([scores[k] for k in dim_keys])]}")

# Differential Gene Expression — Volcano Plot
np.random.seed(42)
genes = ["SST", "PVALB", "GAD1", "GRIN2A", "GRIN2B", "BDNF", "ARC", "TREM2"]
fold_changes = {"SST": -2.3, "PVALB": -2.8, "GAD1": -1.6, "GRIN2A": -1.2, "GRIN2B": -0.9, "BDNF": -1.8, "ARC": -2.0, "TREM2": 1.4}

n_samples = 25
results = []
for gene in genes:
    fc = fold_changes.get(gene, np.random.uniform(-1.5, 1.5))
    control = np.random.normal(loc=8.0, scale=0.6, size=n_samples)
    disease = np.random.normal(loc=8.0 + fc, scale=0.8, size=n_samples)
    t_stat, p_val = stats.ttest_ind(control, disease)
    log2fc = np.mean(disease) - np.mean(control)
    results.append({
        'gene': gene, 'log2fc': log2fc, 'p_value': max(p_val, 1e-12),
        'neg_log10_p': -np.log10(max(p_val, 1e-12)),
        'control_mean': np.mean(control), 'disease_mean': np.mean(disease),
    })

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 7))

# Volcano plot
for r in results:
    is_sig = abs(r['log2fc']) > 0.5 and r['p_value'] < 0.05
    color = '#ef5350' if r['log2fc'] > 0.5 and is_sig else '#4fc3f7' if r['log2fc'] < -0.5 and is_sig else '#555'
    ax1.scatter(r['log2fc'], r['neg_log10_p'], c=color, s=120, alpha=0.85, edgecolors='#333', zorder=3)
    ax1.annotate(r['gene'], (r['log2fc'], r['neg_log10_p']), fontsize=8, color='#e0e0e0',
                 xytext=(5, 5), textcoords='offset points')

ax1.axhline(y=-np.log10(0.05), color='#ffd54f', linestyle='--', alpha=0.5, label='p=0.05')
ax1.axvline(x=-0.5, color='#888', linestyle='--', alpha=0.3)
ax1.axvline(x=0.5, color='#888', linestyle='--', alpha=0.3)
ax1.set_xlabel('log2(Fold Change)', fontsize=12)
ax1.set_ylabel('-log10(p-value)', fontsize=12)
ax1.set_title('Volcano Plot: Disease vs Control', fontsize=14, color='#4fc3f7', fontweight='bold')
ax1.legend(fontsize=9, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')

# Expression barplot
x = np.arange(len(genes))
width = 0.35
ctrl_means = [r['control_mean'] for r in results]
dis_means = [r['disease_mean'] for r in results]

ax2.bar(x - width/2, ctrl_means, width, label='Control', color='#4fc3f7', alpha=0.8, edgecolor='#333')
ax2.bar(x + width/2, dis_means, width, label='Disease', color='#ef5350', alpha=0.8, edgecolor='#333')
ax2.set_xticks(x)
ax2.set_xticklabels(genes, rotation=45, ha='right', fontsize=9)
ax2.set_ylabel('Expression Level (log2 CPM)', fontsize=11)
ax2.set_title('Gene Expression: Control vs Disease', fontsize=14, color='#4fc3f7', fontweight='bold')
ax2.legend(fontsize=10, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')

plt.tight_layout()
plt.show()

# Summary table
print("\nDifferential Expression Summary")
print("=" * 70)
print(f"{'Gene':<12} {'log2FC':>10} {'p-value':>12} {'Direction':>12} {'Significant':>12}")
print("-" * 70)
for r in sorted(results, key=lambda x: x['p_value']):
    sig = 'YES ***' if abs(r['log2fc']) > 1.0 and r['p_value'] < 0.001 else 'YES *' if abs(r['log2fc']) > 0.5 and r['p_value'] < 0.05 else 'no'
    direction = 'UP' if r['log2fc'] > 0 else 'DOWN'
    print(f"{r['gene']:<12} {r['log2fc']:>10.3f} {r['p_value']:>12.2e} {direction:>12} {sig:>12}")

# Pathway Enrichment Analysis
pathways_data = [["GABAergic Synaptic Transmission", 8.8, 5e-08, 6], ["Gamma Oscillation Generation", 8.1, 2e-07, 5], ["Microglial Phagocytosis", 7.4, 8e-07, 4], ["Synaptic Plasticity (LTP)", 6.9, 3e-06, 5], ["NMDA Receptor Signaling", 6.3, 1e-05, 4], ["BDNF-TrkB Pathway", 5.8, 4e-05, 3], ["Neurovascular Coupling", 5.3, 0.0001, 3], ["Theta-Gamma Coupling", 4.9, 0.0003, 4], ["Glymphatic Clearance", 4.4, 0.0008, 3], ["Calcium Homeostasis", 4.0, 0.002, 4], ["PV+ Interneuron Function", 3.5, 0.005, 2], ["Tau Phosphorylation", 3.1, 0.01, 3]]

pathways = [p[0] for p in pathways_data]
enrichment_scores = [p[1] for p in pathways_data]
p_values = [p[2] for p in pathways_data]
gene_counts = [p[3] for p in pathways_data]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))

# Dot plot
sizes = [gc * 40 for gc in gene_counts]
colors = [-np.log10(p) for p in p_values]
scatter = ax1.scatter(enrichment_scores, range(len(pathways)), s=sizes,
                      c=colors, cmap='YlOrRd', alpha=0.85, edgecolors='#333')
ax1.set_yticks(range(len(pathways)))
ax1.set_yticklabels(pathways, fontsize=9)
ax1.set_xlabel('Enrichment Score', fontsize=11)
ax1.set_title('Pathway Enrichment (Dot Plot)', fontsize=14, color='#4fc3f7', fontweight='bold')
cbar = plt.colorbar(scatter, ax=ax1, shrink=0.6)
cbar.set_label('-log10(p-value)', fontsize=9, color='#e0e0e0')
cbar.ax.yaxis.set_tick_params(color='#888')

# Significance bars
bar_colors = ['#ef5350' if p < 0.001 else '#ff8a65' if p < 0.01 else '#ffd54f' if p < 0.05 else '#888'
              for p in p_values]
ax2.barh(range(len(pathways)), [-np.log10(p) for p in p_values],
         color=bar_colors, alpha=0.85, edgecolor='#333')
ax2.set_yticks(range(len(pathways)))
ax2.set_yticklabels(pathways, fontsize=9)
ax2.set_xlabel('-log10(p-value)', fontsize=11)
ax2.set_title('Statistical Significance', fontsize=14, color='#4fc3f7', fontweight='bold')
ax2.axvline(x=-np.log10(0.05), color='#ffd54f', linestyle='--', alpha=0.7, label='p=0.05')
ax2.axvline(x=-np.log10(0.001), color='#ef5350', linestyle='--', alpha=0.7, label='p=0.001')
ax2.legend(fontsize=9, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')

plt.tight_layout()
plt.show()

print("\nPathway Enrichment Summary")
print("=" * 80)
print(f"{'Pathway':<35} {'Enrichment':>12} {'p-value':>12} {'Genes':>8}")
print("-" * 80)
for pw, es, pv, gc in zip(pathways, enrichment_scores, p_values, gene_counts):
    print(f"{pw:<35} {es:>12.2f} {pv:>12.2e} {gc:>8}")

# Statistical Analysis — Hypothesis Score Profile
scores = {"composite": 0.767541, "mech": 0.85, "evid": 0.82, "novel": 0.78, "feas": 0.88, "impact": 0.8, "drug": 0.75, "safety": 0.9, "comp": 0.7, "data": 0.85, "reprod": 0.82}

dim_keys = ['mech', 'evid', 'novel', 'feas', 'impact', 'drug', 'safety', 'comp', 'data', 'reprod']
dim_labels = ['Mechanistic', 'Evidence', 'Novelty', 'Feasibility', 'Impact',
              'Druggability', 'Safety', 'Competition', 'Data Avail.', 'Reproducibility']

values = np.array([scores[k] for k in dim_keys])

print("=" * 60)
print("STATISTICAL PROFILE ANALYSIS")
print("=" * 60)

print(f"\nComposite Score: {scores.get('composite', 0):.4f}")
print(f"Mean dimension score: {np.mean(values):.4f}")
print(f"Median: {np.median(values):.4f}")
print(f"Std Dev: {np.std(values):.4f}")
print(f"CV (coefficient of variation): {np.std(values)/np.mean(values)*100:.1f}%")
print(f"Range: {np.min(values):.3f} — {np.max(values):.3f}")

print(f"\nSTRENGTH ANALYSIS")
print("-" * 60)
for i, (dim, val) in enumerate(zip(dim_labels, values)):
    bar = '█' * int(val * 30) + '░' * (30 - int(val * 30))
    label = '★★★' if val >= 0.9 else '★★' if val >= 0.8 else '★' if val >= 0.7 else ''
    print(f"{dim:<20} {bar} {val:.2f} {label}")

# Percentile analysis (compared to typical hypothesis scores)
print(f"\nPERCENTILE ANALYSIS (vs typical SciDEX hypotheses)")
print("-" * 60)
# Typical score distribution centers around 0.6 with std 0.15
for dim, val in zip(dim_labels, values):
    percentile = stats.norm.cdf(val, loc=0.6, scale=0.15) * 100
    print(f"{dim:<20} Score: {val:.2f}  Percentile: {percentile:>5.1f}th")

# Correlation with composite
print(f"\nDIMENSION CONTRIBUTION ANALYSIS")
print("-" * 60)
# Estimate contribution as deviation from mean * weight
weights = np.array([0.15, 0.12, 0.10, 0.10, 0.12, 0.10, 0.08, 0.08, 0.08, 0.07])
contributions = values * weights
norm_contrib = contributions / contributions.sum() * 100
for dim, val, contrib in zip(dim_labels, values, norm_contrib):
    print(f"{dim:<20} {val:.2f} → {contrib:>5.1f}% of composite")