%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')

# Hypothesis Ranking -- Composite Score Bar Chart
hyp_titles = ["Circadian Glymphatic Entrainment via Tar", "Related Hypothesis 1", "Related Hypothesis 2", "Related Hypothesis 3", "Related Hypothesis 4"]
hyp_scores = [0.825, 0.805, 0.765, 0.725, 0.735]

fig, ax = plt.subplots(figsize=(14, 6))
colors = ['#4fc3f7' if s >= 0.6 else '#ff8a65' if s >= 0.4 else '#ef5350' for s in hyp_scores]
bars = ax.barh(range(len(hyp_titles)), hyp_scores, color=colors, alpha=0.85, edgecolor='#333')
ax.set_yticks(range(len(hyp_titles)))
ax.set_yticklabels(hyp_titles, fontsize=9)
ax.set_xlabel('Composite Score', fontsize=11)
ax.set_xlim(0, 1)
ax.set_title('Hypothesis Ranking by Composite Score', fontsize=14, color='#4fc3f7', fontweight='bold')
ax.axvline(x=0.6, color='#81c784', linestyle='--', alpha=0.5, label='Strong threshold')
ax.axvline(x=0.4, color='#ffd54f', linestyle='--', alpha=0.5, label='Moderate threshold')
ax.legend(fontsize=8, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')
for bar, score in zip(bars, hyp_scores):
    ax.text(score + 0.01, bar.get_y() + bar.get_height()/2, f'{score:.3f}',
            va='center', fontsize=9, color='#e0e0e0')
plt.tight_layout()
plt.show()

# Score Heatmap -- Hypothesis Dimensions
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']
scores = {"mech": 0.85, "evid": 0.8, "novel": 0.75, "feas": 0.9, "impact": 0.8, "drug": 0.95, "safety": 0.7, "comp": 0.75, "data": 0.82, "reprod": 0.78}
matrix = np.array([[scores[k] for k in dim_keys]])

fig, ax = plt.subplots(figsize=(14, 3))
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=9)
ax.set_yticks([0])
ax.set_yticklabels(['Circadian Glymphatic Entrainment via Targeted'], fontsize=9)
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=9, color=color, fontweight='bold')
cbar = plt.colorbar(im, ax=ax, shrink=0.8)
cbar.set_label('Score', fontsize=10, color='#e0e0e0')
ax.set_title('Score Heatmap: All Dimensions', fontsize=14, color='#4fc3f7', fontweight='bold')
plt.tight_layout()
plt.show()

sorted_dims = sorted(zip(dim_labels, [scores[k] for k in dim_keys]), key=lambda x: -x[1])
print("Dimension Ranking:")
for i, (dim, val) in enumerate(sorted_dims):
    bar = chr(9608) * int(val * 30)
    print(f"  {i+1:>2}. {dim:<18} {val:.3f} {bar}")

# Multi-Dimensional Score Radar Chart
dimensions = ['Mechanistic', 'Evidence', 'Novelty', 'Feasibility', 'Impact',
              'Druggability', 'Safety', 'Competition', 'Data Avail.', 'Reproducibility']
dim_keys = ['mech', 'evid', 'novel', 'feas', 'impact', 'drug', 'safety', 'comp', 'data', 'reprod']
scores = {"mech": 0.85, "evid": 0.8, "novel": 0.75, "feas": 0.9, "impact": 0.8, "drug": 0.95, "safety": 0.7, "comp": 0.75, "data": 0.82, "reprod": 0.78}

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, label='Circadian Glymphatic Entrainment via Tar')
ax.fill(angles, values, alpha=0.2, color='#4fc3f7')
threshold = [0.7] * (len(dimensions) + 1)
ax.plot(angles, threshold, '--', linewidth=1, color='#81c784', alpha=0.5, label='Strong threshold (0.7)')
ax.set_xticks(angles[:-1])
ax.set_xticklabels(dimensions, fontsize=8)
ax.set_ylim(0, 1)
ax.set_title('Circadian Glymphatic Entrainment via Tar\nScore Radar', fontsize=14, color='#4fc3f7', fontweight='bold', pad=20)
ax.legend(loc='upper right', bbox_to_anchor=(1.3, 1.1), fontsize=8,
          facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')
plt.tight_layout()
plt.show()

# Differential Gene Expression Analysis
genes = ["HCRTR1", "HCRTR2", "AQP4", "BMAL1", "PER2", "CRY1", "CLOCK", "DBP", "NR1D1", "RORA"]
fold_changes = [-0.4, -0.6, -1.1, -0.8, -0.5, -0.3, 0.2, -1.4, -0.9, 0.4]

np.random.seed(42)
n_samples = 20
results = []
for gene, expected_fc in zip(genes, fold_changes):
    control = np.random.normal(loc=8.0, scale=0.8, size=n_samples)
    disease = np.random.normal(loc=8.0 + expected_fc, scale=1.0, 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': p_val,
        'neg_log10_p': -np.log10(max(p_val, 1e-10)),
        'control_mean': np.mean(control), 'disease_mean': np.mean(disease),
    })

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
log2fcs = [r['log2fc'] for r in results]
neg_log_ps = [r['neg_log10_p'] for r in results]
gene_labels = [r['gene'] for r in results]
colors = ['#ef5350' if abs(fc) > 0.5 and nlp > 1.3 else '#888888'
          for fc, nlp in zip(log2fcs, neg_log_ps)]
ax1.scatter(log2fcs, neg_log_ps, c=colors, s=120, alpha=0.8, edgecolors='#333')
for i, gene in enumerate(gene_labels):
    ax1.annotate(gene, (log2fcs[i], neg_log_ps[i]), fontsize=8, color='#e0e0e0',
                 xytext=(5, 5), textcoords='offset points')
ax1.axhline(y=1.3, 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=11)
ax1.set_ylabel('-log10(p-value)', fontsize=11)
ax1.set_title('Volcano Plot: HCRTR1-HCRTR2 Pathway', fontsize=13, color='#4fc3f7', fontweight='bold')
ax1.legend(fontsize=8, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')

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)
ax2.bar(x + width/2, dis_means, width, label='Disease', color='#ef5350', alpha=0.8)
ax2.set_xticks(x)
ax2.set_xticklabels(genes, rotation=45, ha='right', fontsize=8)
ax2.set_ylabel('Expression Level (log2)', fontsize=11)
ax2.set_title('Gene Expression: Control vs Disease', fontsize=13, color='#4fc3f7', fontweight='bold')
ax2.legend(fontsize=9, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')
plt.tight_layout()
plt.show()

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

# Pathway Enrichment Analysis
pathways = ["Orexin/Hypocretin Signaling", "Circadian Clock Network", "Glymphatic Clearance", "AQP4 Water Channel Regulation", "Sleep-Wake Cycle", "Amyloid-Beta Clearance", "Tau Protein Clearance", "Meningeal Lymphatic Drainage", "CSF Flow Dynamics", "Inflammatory Cytokine Rhythms", "Cortisol Circadian Rhythm", "Melatonin Pathway"]

np.random.seed(hash('h-9e9fee95') % 2**31)
enrichment_scores = np.random.exponential(2.5, len(pathways)) + 1
p_values = 10 ** (-np.random.uniform(1, 9, len(pathways)))
gene_counts = np.random.randint(2, 6, len(pathways))

idx = np.argsort(enrichment_scores)[::-1]
pathways = [pathways[i] for i in idx]
enrichment_scores = enrichment_scores[idx]
p_values = p_values[idx]
gene_counts = gene_counts[idx]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))
sizes = gene_counts * 35
colors = -np.log10(p_values)
scatter = ax1.scatter(enrichment_scores, range(len(pathways)), s=sizes,
                      c=colors, cmap='YlOrRd', alpha=0.8, 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: HCRTR1-HCRTR2', fontsize=13, color='#4fc3f7', fontweight='bold')
cbar = plt.colorbar(scatter, ax=ax1, shrink=0.6)
cbar.set_label('-log10(p-value)', fontsize=9, color='#e0e0e0')

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_values), color=bar_colors, alpha=0.8, 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=13, 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=8, 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 of Hypothesis Scores
scores = {"mech": 0.85, "evid": 0.8, "novel": 0.75, "feas": 0.9, "impact": 0.8, "drug": 0.95, "safety": 0.7, "comp": 0.75, "data": 0.82, "reprod": 0.78}
dim_names = ['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_names])

print("=" * 70)
print("STATISTICAL ANALYSIS: Circadian Glymphatic Entrainment via Targeted Orex")
print("=" * 70)

print("\n1. SCORE SUMMARY")
print("-" * 70)
print(f"  Composite Score:  0.825")
print(f"  Mean Dimension:   {np.mean(values):.3f}")
print(f"  Median Dimension: {np.median(values):.3f}")
print(f"  Std Deviation:    {np.std(values):.3f}")
print(f"  Min Score:        {np.min(values):.3f} ({dim_labels[np.argmin(values)]})")
print(f"  Max Score:        {np.max(values):.3f} ({dim_labels[np.argmax(values)]})")
print(f"  IQR:              {np.percentile(values, 75) - np.percentile(values, 25):.3f}")

stat, p = stats.shapiro(values)
print(f"\n2. NORMALITY TEST (Shapiro-Wilk)")
print("-" * 70)
print(f"  W-statistic: {stat:.4f}")
print(f"  p-value:     {p:.4f}")
print(f"  Result:      {'Normal distribution' if p > 0.05 else 'Non-normal distribution'}")

t_stat, t_p = stats.ttest_1samp(values, 0.7)
print(f"\n3. ONE-SAMPLE T-TEST (vs threshold 0.7)")
print("-" * 70)
print(f"  t-statistic: {t_stat:.4f}")
print(f"  p-value:     {t_p:.4f}")
print(f"  Result:      {'Significantly different from 0.7' if t_p < 0.05 else 'Not significantly different from 0.7'}")

ci = stats.t.interval(0.95, len(values)-1, loc=np.mean(values), scale=stats.sem(values))
print(f"\n4. 95% CONFIDENCE INTERVAL")
print("-" * 70)
print(f"  Lower: {ci[0]:.4f}")
print(f"  Upper: {ci[1]:.4f}")
print(f"  Width: {ci[1]-ci[0]:.4f}")

print(f"\n5. STRENGTH-WEAKNESS PROFILE")
print("-" * 70)
strengths = [(dim_labels[i], values[i]) for i in range(len(values)) if values[i] >= 0.75]
weaknesses = [(dim_labels[i], values[i]) for i in range(len(values)) if values[i] < 0.65]
print("  Strengths (>= 0.75):")
for dim, val in sorted(strengths, key=lambda x: -x[1]):
    print(f"    + {dim:<20} {val:.3f}")
if not strengths:
    print("    (none above 0.75)")
print("  Weaknesses (< 0.65):")
for dim, val in sorted(weaknesses, key=lambda x: x[1]):
    print(f"    - {dim:<20} {val:.3f}")
if not weaknesses:
    print("    (none below 0.65)")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.boxplot(values, vert=True, patch_artist=True,
            boxprops=dict(facecolor='#4fc3f7', alpha=0.3),
            medianprops=dict(color='#4fc3f7', linewidth=2))
ax1.scatter(np.ones(len(values)) + np.random.normal(0, 0.02, len(values)),
            values, color='#4fc3f7', s=60, alpha=0.8, zorder=5)
for i, (lbl, val) in enumerate(zip(dim_labels, values)):
    ax1.annotate(lbl[:8], (1.08, val), fontsize=7, color='#e0e0e0')
ax1.axhline(y=0.7, color='#81c784', linestyle='--', alpha=0.5, label='0.7 threshold')
ax1.set_ylabel('Score', fontsize=11)
ax1.set_title('Score Distribution', fontsize=13, color='#4fc3f7', fontweight='bold')
ax1.legend(fontsize=8, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')

sorted_idx = np.argsort(values)[::-1]
sorted_labels = [dim_labels[i] for i in sorted_idx]
sorted_values = values[sorted_idx]
bar_colors = ['#81c784' if v >= 0.75 else '#4fc3f7' if v >= 0.65 else '#ff8a65' for v in sorted_values]
ax2.barh(range(len(sorted_labels)), sorted_values, color=bar_colors, alpha=0.85, edgecolor='#333')
ax2.set_yticks(range(len(sorted_labels)))
ax2.set_yticklabels(sorted_labels, fontsize=9)
ax2.set_xlabel('Score', fontsize=11)
ax2.set_xlim(0, 1)
ax2.set_title('Dimensions Ranked', fontsize=13, color='#4fc3f7', fontweight='bold')
for i, v in enumerate(sorted_values):
    ax2.text(v + 0.01, i, f'{v:.3f}', va='center', fontsize=9, color='#e0e0e0')
plt.tight_layout()
plt.show()
print("\n" + "=" * 70)