Visualising SimPy Results: Making Data Speak

import simpy
import random
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# Style
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['figure.figsize'] = (10, 6)

def plot_queue_length(time_series_df):
    plt.figure()
    plt.step(time_series_df['time'], time_series_df['queue_length'],
             where='post', linewidth=1)
    plt.fill_between(time_series_df['time'], time_series_df['queue_length'],
                     step='post', alpha=0.3)
    plt.xlabel('Time')
    plt.ylabel('Queue Length')
    plt.title('Queue Length Over Time')
    plt.tight_layout()
    plt.savefig('queue_length.png', dpi=150)
    plt.show()

def plot_wait_distribution(wait_times):
    plt.figure()
    plt.hist(wait_times, bins=30, edgecolor='white', alpha=0.7)
    plt.axvline(np.mean(wait_times), color='red', linestyle='--',
                label=f'Mean: {np.mean(wait_times):.2f}')
    plt.axvline(np.percentile(wait_times, 95), color='orange', linestyle='--',
                label=f'95th percentile: {np.percentile(wait_times, 95):.2f}')
    plt.xlabel('Wait Time')
    plt.ylabel('Frequency')
    plt.title('Distribution of Wait Times')
    plt.legend()
    plt.tight_layout()
    plt.savefig('wait_distribution.png', dpi=150)
    plt.show()

def plot_utilisation(time_series_df, resource_name='server'):
    plt.figure()
    plt.plot(time_series_df['time'],
             time_series_df[f'{resource_name}_utilisation'],
             linewidth=1)
    plt.fill_between(time_series_df['time'],
                     time_series_df[f'{resource_name}_utilisation'],
                     alpha=0.3)
    avg_util = time_series_df[f'{resource_name}_utilisation'].mean()
    plt.axhline(avg_util, color='red', linestyle='--',
                label=f'Average: {avg_util:.1%}')
    plt.xlabel('Time')
    plt.ylabel('Utilisation')
    plt.ylim(0, 1.1)
    plt.title(f'{resource_name.title()} Utilisation Over Time')
    plt.legend()
    plt.tight_layout()
    plt.savefig('utilisation.png', dpi=150)
    plt.show()

def plot_scenario_comparison(results_dict):
    """results_dict = {'Scenario A': wait_times_a, 'Scenario B': wait_times_b}"""
    plt.figure()
    data = list(results_dict.values())
    labels = list(results_dict.keys())
    plt.boxplot(data, labels=labels)
    plt.ylabel('Wait Time')
    plt.title('Wait Time Comparison Across Scenarios')
    plt.tight_layout()
    plt.savefig('comparison.png', dpi=150)
    plt.show()

def plot_throughput(entity_df, window=50):
    """Cumulative completions and rolling throughput."""
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

    # Cumulative completions
    entity_df = entity_df.sort_values('departure')
    entity_df['cumulative'] = range(1, len(entity_df) + 1)
    ax1.plot(entity_df['departure'], entity_df['cumulative'])
    ax1.set_xlabel('Time')
    ax1.set_ylabel('Cumulative Completions')
    ax1.set_title('Cumulative Throughput')

    # Rolling throughput
    entity_df['throughput'] = entity_df['cumulative'] / entity_df['departure']
    ax2.plot(entity_df['departure'], entity_df['throughput'])
    ax2.set_xlabel('Time')
    ax2.set_ylabel('Throughput (per time unit)')
    ax2.set_title('Throughput Rate')

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

def plot_multiple_resources(time_series_df, resources):
    fig, axes = plt.subplots(len(resources), 1, figsize=(10, 4*len(resources)))

    for ax, resource in zip(axes, resources):
        ax.step(time_series_df['time'],
                time_series_df[f'{resource}_queue'],
                where='post', label='Queue')
        ax.step(time_series_df['time'],
                time_series_df[f'{resource}_busy'],
                where='post', label='In Service')
        ax.set_xlabel('Time')
        ax.set_ylabel('Count')
        ax.set_title(f'{resource.title()}')
        ax.legend()

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

def plot_arrival_heatmap(entity_df, bins_x=24, bins_y=7):
    """Show arrival patterns (e.g., hour of day vs day of week)."""
    # Assuming arrival times can be converted to hour/day
    entity_df['hour'] = (entity_df['arrival'] % 24).astype(int)
    entity_df['day'] = (entity_df['arrival'] // 24 % 7).astype(int)

    heatmap_data = entity_df.groupby(['day', 'hour']).size().unstack(fill_value=0)

    plt.figure(figsize=(12, 5))
    plt.imshow(heatmap_data, aspect='auto', cmap='YlOrRd')
    plt.colorbar(label='Arrivals')
    plt.xlabel('Hour of Day')
    plt.ylabel('Day of Week')
    plt.title('Arrival Pattern Heatmap')
    plt.tight_layout()
    plt.savefig('heatmap.png', dpi=150)
    plt.show()

def plot_cdf(wait_times, label='Wait Time'):
    sorted_data = np.sort(wait_times)
    cdf = np.arange(1, len(sorted_data) + 1) / len(sorted_data)

    plt.figure()
    plt.plot(sorted_data, cdf, linewidth=2)
    plt.xlabel(label)
    plt.ylabel('Cumulative Probability')
    plt.title(f'Cumulative Distribution of {label}')

    # Mark key percentiles
    for p in [50, 90, 95, 99]:
        val = np.percentile(wait_times, p)
        plt.axhline(p/100, color='gray', linestyle=':', alpha=0.5)
        plt.axvline(val, color='gray', linestyle=':', alpha=0.5)
        plt.annotate(f'{p}%: {val:.1f}', xy=(val, p/100),
                    xytext=(val + 1, p/100 - 0.05))

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

def plot_confidence_intervals(replication_results):
    """replication_results = list of mean wait times from each replication."""
    n = len(replication_results)
    mean = np.mean(replication_results)
    std = np.std(replication_results, ddof=1)
    ci = 1.96 * std / np.sqrt(n)

    plt.figure()
    plt.bar(['Mean Wait Time'], [mean], yerr=[ci], capsize=10, color='steelblue')
    plt.ylabel('Time')
    plt.title(f'Mean Wait Time with 95% CI\n(n={n} replications)')
    plt.tight_layout()
    plt.savefig('confidence.png', dpi=150)
    plt.show()

import plotly.express as px
import plotly.graph_objects as go

def interactive_queue(time_series_df):
    fig = px.line(time_series_df, x='time', y='queue_length',
                  title='Queue Length Over Time')
    fig.update_traces(line_shape='hv')  # Step plot
    fig.write_html('queue_interactive.html')
    fig.show()

def generate_all_plots(entity_df, time_series_df, resources):
    """Generate all standard visualisations."""

    # 1. Queue length
    for resource in resources:
        plot_queue_length(time_series_df.rename(
            columns={f'{resource}_queue': 'queue_length'}
        ))

    # 2. Wait distribution
    plot_wait_distribution(entity_df['wait'].values)

    # 3. Utilisation
    for resource in resources:
        plot_utilisation(time_series_df, resource)

    # 4. CDF
    plot_cdf(entity_df['wait'].values)

    # 5. Throughput
    plot_throughput(entity_df)

    print("All plots generated!")