import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from scipy.stats import gaussian_kde

plt.rcParams['svg.fonttype'] = 'none'  # <-- add this line

# Configuration
ROPE_CENTER = 0  # Center of ROPE band (in wins per 162 scale)
ROPE_HALFWIDTH = 2  # Default: ±2 wins per 162 games

# Grid layout (rows, cols)
GRID_SHAPE = (3, 3)  # 3x3 grid

# X-axis range in wins per 162 games
X_RANGE = (-10, 10)

# ========== HDI FUNCTION ==========
def calculate_hdi(draws, credible_mass=0.95):
    """
    Calculate the Highest Density Interval (HDI) - the shortest interval
    containing the specified credible mass.
    """
    sorted_draws = np.sort(draws)
    n = len(sorted_draws)
    interval_size = int(np.ceil(credible_mass * n))
    n_intervals = n - interval_size + 1
    
    # Calculate width of each possible interval
    interval_widths = sorted_draws[interval_size-1:] - sorted_draws[:n_intervals]
    
    # Find the narrowest interval
    min_idx = np.argmin(interval_widths)
    hdi_lower = sorted_draws[min_idx]
    hdi_upper = sorted_draws[min_idx + interval_size - 1]
    
    return hdi_lower, hdi_upper

# ========== LOAD DATA ==========
print("Loading posterior draws...")
data_path = r'[df_6]'
df = pd.read_csv(data_path)

# Convert to wins per 162 games, centered on average manager
df['prob_raw'] = 1 / (1 + np.exp(-df['u']))
mean_prob = df['prob_raw'].mean()  # Empirical mean across all draws
df['w162'] = (df['prob_raw'] - mean_prob) * 162

print(f"Loaded {len(df)} posterior draws for {df['mgrID'].nunique()} managers")
print()

# Get list of unique managers
managers = sorted(df['mgrID'].unique())
print(f"Available managers: {len(managers)}")
print("Sample managers:", managers[:10])
print()

# ========== INPUT MANAGER IDs ==========
print("Enter manager IDs for each row (3 managers per row)")
print("Format: mgrID1 mgrID2 mgrID3")
print("Example: mcgrj101 chanf101 clarf101")
print()

# Ask for ROPE halfwidth
rope_input = input(f"ROPE halfwidth in wins per 162 (default={ROPE_HALFWIDTH}): ").strip()
if rope_input:
    ROPE_HALFWIDTH = float(rope_input)
print(f"Using ROPE: ±{ROPE_HALFWIDTH} wins per 162 games")
print()

selected_managers = []

for row_num in range(1, 4):
    while True:
        row_input = input(f"Row {row_num} (managers {(row_num-1)*3+1}-{row_num*3}): ")
        try:
            # Parse the input - split by whitespace
            mgrs_in_row = row_input.split()
            if len(mgrs_in_row) != 3:
                print(f"  Error: Expected 3 managers, got {len(mgrs_in_row)}. Try again.")
                continue
            
            # Check if all managers exist
            invalid_mgrs = [m for m in mgrs_in_row if m not in managers]
            if invalid_mgrs:
                print(f"  Error: Manager(s) not found: {invalid_mgrs}. Try again.")
                continue
            
            selected_managers.extend(mgrs_in_row)
            print(f"  ✓ Row {row_num} added successfully")
            break
        except Exception as e:
            print(f"  Error: {e}. Try again.")
    print()

print(f"All {len(selected_managers)} managers selected!")
print("Selected:", selected_managers)
print()

# ========== CALCULATE ROPE STATISTICS ==========
print("Calculating ROPE statistics for each manager...")

rope_stats = []
for mgr in selected_managers:
    mgr_draws = df[df['mgrID'] == mgr]['w162'].values
    
    # Calculate ROPE probabilities
    p_above = np.mean(mgr_draws > ROPE_HALFWIDTH)
    p_below = np.mean(mgr_draws < -ROPE_HALFWIDTH)
    p_in = np.mean((mgr_draws >= -ROPE_HALFWIDTH) & (mgr_draws <= ROPE_HALFWIDTH))
    
    # Calculate HDI
    hdi_lower, hdi_upper = calculate_hdi(mgr_draws, credible_mass=0.95)
    
    # Calculate MAP (mode) using KDE
    kde = gaussian_kde(mgr_draws)
    x_grid = np.linspace(mgr_draws.min() - 1, mgr_draws.max() + 1, 1000)
    density_vals = kde(x_grid)
    map_estimate = x_grid[np.argmax(density_vals)]
    
    rope_stats.append({
        'mgrID': mgr,
        'p_above': p_above,
        'p_in': p_in,
        'p_below': p_below,
        'hdi_lower': hdi_lower,
        'hdi_upper': hdi_upper,
        'mean': mgr_draws.mean(),
        'median': np.median(mgr_draws),
        'map': map_estimate
    })

rope_df = pd.DataFrame(rope_stats)
print("\nROPE Statistics:")
print(rope_df[['mgrID', 'p_above', 'p_in', 'p_below', 'mean', 'median', 'map']])
print()

# ========== PLOTTING FUNCTION ==========
def plot_rope_grid_hb(df, selected_managers, grid_shape, rope_center=0, rope_hw=2, x_range=(-10, 10)):
    """
    Create a grid of posterior density plots using actual draws with continuous ROPE band
    
    Parameters:
    -----------
    df : DataFrame
        Posterior draws with columns: mgrID, w162
    selected_managers : list
        List of manager IDs to plot (in order)
    grid_shape : tuple
        (rows, cols) for subplot grid
    rope_center : float
        Center of ROPE region in wins per 162
    rope_hw : float
        Half-width of ROPE region in wins per 162
    x_range : tuple
        (min, max) for x-axis range in wins per 162
    """
    n_rows, n_cols = grid_shape
    
    # Create figure with minimal spacing between subplots
    fig, axes = plt.subplots(n_rows, n_cols, figsize=(14, 11))
    axes = axes.flatten()
    
    rope_left = rope_center - rope_hw
    rope_right = rope_center + rope_hw
    
    for idx, mgr in enumerate(selected_managers):
        if idx >= len(axes):
            break
            
        ax = axes[idx]
        
        # Get posterior draws for this manager
        mgr_draws = df[df['mgrID'] == mgr]['w162'].values
        
        # Create x range for density plot
        x = np.linspace(x_range[0], x_range[1], 1000)
        
        # Use KDE for actual posterior density
        kde = gaussian_kde(mgr_draws)
        density = kde(x)
        
        # Calculate 95% HDI
        hdi_lower, hdi_upper = calculate_hdi(mgr_draws, credible_mass=0.95)
        
        # Plot density curve
        ax.plot(x, density, 'k-', linewidth=1.5, zorder=3)
        
        # Fill under curve
        ax.fill_between(x, density, alpha=0.6, color='#00a3ff', zorder=2)
        
        # Add HDI line at bottom
        hdi_y = -0.02 * density.max()
        ax.plot([hdi_lower, hdi_upper], [hdi_y, hdi_y], 'k-', linewidth=3, zorder=5)
        
        # Set consistent x-axis range
        ax.set_xlim(x_range)
        
        # Force x-axis ticks
        ax.set_xticks([-10, -5, 0, 5, 10])
        
        # Minimal formatting
        ax.set_xlabel('')
        ax.set_ylabel('')
        ax.tick_params(axis='y', which='both', left=False, labelleft=False)
        ax.tick_params(axis='x', labelsize=18)
        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.spines['left'].set_visible(False)
        
        # Set background to white
        ax.set_facecolor('white')
    
    # Hide unused subplots
    for idx in range(len(selected_managers), len(axes)):
        axes[idx].axis('off')
    
    # Adjust spacing
    plt.tight_layout()
    
    # Add continuous vertical ROPE bands for each column
    from matplotlib.patches import Rectangle
    
    for col in range(n_cols):
        first_in_col_idx = col
        if first_in_col_idx >= len(selected_managers):
            continue
            
        ax_ref = axes[first_in_col_idx]
        
        # Get the bbox of this subplot
        bbox = ax_ref.get_position()
        
        # Convert ROPE boundaries from data coordinates to figure coordinates
        x_fig_left = ax_ref.transData.transform([(rope_left, 0)])[0][0]
        x_fig_right = ax_ref.transData.transform([(rope_right, 0)])[0][0]
        
        # Convert to figure coordinates
        x_fig_left = fig.transFigure.inverted().transform([(x_fig_left, 0)])[0][0]
        x_fig_right = fig.transFigure.inverted().transform([(x_fig_right, 0)])[0][0]
        
        # Add rectangle spanning full height
        rope_band = Rectangle((x_fig_left, 0), x_fig_right - x_fig_left, 1,
                              transform=fig.transFigure, 
                              color='gray', alpha=0.3, zorder=0)
        fig.patches.append(rope_band)
    
    return fig

# ========== CREATE PLOT ==========
print("Creating visualization...")
fig = plot_rope_grid_hb(df, selected_managers, GRID_SHAPE, ROPE_CENTER, ROPE_HALFWIDTH, X_RANGE)

# Save the figure as SVG
output_path = r'C:\Users\dmk38\Documents\x5\rope_hb_plot.svg'
fig.savefig(output_path, format='svg', bbox_inches='tight')
print(f"\n✓ Plot saved to: {output_path}")

# Save ROPE statistics
stats_path = r'C:\Users\dmk38\Documents\x5\rope_hb_stats.csv'
rope_df.to_csv(stats_path, index=False)
print(f"✓ ROPE statistics saved to: {stats_path}")

plt.show()