mri-schedule.py

from collections import namedtuple
from pprint import pprint

import numpy as np
np.random.seed(12345)

import matplotlib.pyplot as plt


Appointment = namedtuple("Appointment", "scheduled late_by emergency")
DayStats = namedtuple("DayStats", "occupancy patient_wait_times emergency_wait_times")

def next_appointment(appointments, time_now, no_show):
    if not appointments:
        return None
        
    while True:
        if appointments and appointments[0].scheduled + appointments[0].late_by <= time_now:
            app = appointments.pop(0)
            if np.random.uniform() <= no_show:
                continue
            else:
                return app

        else:
            return None

def simulate_day(day=8*60, no_show=0.01, app_duration=30):
    # before the day begins appointments are schedule
    # each one is planned to last 30minutes, including
    # time to switch patients etc. The scheduled
    # duration will be taken as the basis to generate
    # how long it really took
    appointments = [Appointment(n,
                                int(round(np.random.normal(scale=2))),
                                False) for n in np.arange(0,
                                                          day,
                                                          step=app_duration)]
    
    # Time left before this appointment is done
    time_left = 0
    in_use = False
    # minutes the MRI was occupied for during the day
    occupied = 0
    
    # Patient waiting times
    patient_wait_times = []
    emergency_wait_times = []
    
    for minute in xrange(day):
        if np.random.uniform() <= 0.10 / float(day):
            appointments.insert(0, Appointment(minute, 0, True))
        
        if time_left:
            time_left -= 1

        else:
            app = next_appointment(appointments, minute, no_show)
            if app is None:
                in_use = False
            
            else:
                # minutes the patient had to actually wait
                waited = minute - app.scheduled - app.late_by
                if app.emergency:
                    emergency_wait_times.append(waited)
                else:
                    patient_wait_times.append(waited)
                    
                real_duration = app_duration + int(round(np.random.normal(scale=2)))
                in_use = True
            
                # subtract one as the first minute is already over
                time_left = real_duration-1

        if in_use:
            occupied += 1

    occupancy = occupied/float(day)
    
    return DayStats(occupancy, patient_wait_times, emergency_wait_times)

daily_wait = []
occupancies = []
for n in xrange(1000):
    day = simulate_day()
    daily_wait.append(np.mean(day.patient_wait_times))
    occupancies.append(day.occupancy*100)

hist = lambda ax,x: ax.hist(x, histtype="stepfilled",
                            normed=True, bins=20)

with plt.xkcd():
    f, (ax1, ax2) = plt.subplots(1, 2)
    hist(ax1, daily_wait)
    m_wait = np.mean(daily_wait)
    ax1.vlines(m_wait, 0,1)
    ax1.text(1.1*m_wait, 0.8, "%.2f"%(m_wait))
    ax1.set_xlabel("average daily wait [m]")
    
    hist(ax2, occupancies)
    m_occ = np.mean(occupancies)
    ax2.vlines(m_occ, 0,1)
    ax2.text(0.95*m_occ, 0.8, "%.2f %%"%(m_occ))
    ax2.set_xlabel("MRI occupancy [%]")
    plt.show()