fix: harden simulation invariants

README.md

@@ -22,18 +22,21 @@ A fixed seed controls patient attributes, service selection, arrivals and servic
 
 ## Model guarantees
 
-- Exactly `n` patients arrive and are served.
+- Exactly `n` patients arrive and are served; `n = 0` returns an empty completed state.
+- Timing means and standard deviations must be finite and strictly positive; invalid parameters fail with `ArgumentError` before random distributions are built.
 - Every arrival schedules its successor independently of nurse availability.
 - Waiting patients enter service once a nurse becomes free.
 - Service times are positive and patient IDs are unique.
+- Same-time events retain their true timestamp and use deterministic insertion order; simulation time is never adjusted to break a tie.
+- Assigned service intervals do not overlap for an individual nurse.
 
 ## Analysis script
 
-`vccrun.jl` is an optional exploratory sweep and requires additional CSV, DataFrames and Makie packages. It is intentionally separate from the minimal tested core environment.
+`vccrun.jl` is an optional exploratory sweep and requires additional CSV, DataFrames and Makie packages. It is intentionally separate from the minimal tested core environment. Its 100 runs per cell use the deterministic, distinct seed series `1234:1333`; they are reproducible replicates rather than repeated copies of the default trajectory.
 
 ## Scope
 
-The distributions and urgency rules are illustrative modelling assumptions, not a validated clinical staffing model. Validate parameters and outcomes before operational use.
+The distributions, positive-time clamping, and urgency rules are illustrative modelling assumptions, not a validated clinical staffing model. The invariant tests establish program behavior only; they do not validate inputs, outcomes, staffing requirements, or operational decisions. Validate parameters and outcomes independently before any operational use.
 
 ## Licence
 

test/runtests.jl

@@ -17,3 +17,59 @@ const TEST_PARAMETERS = Parameters(20, 2.0, 0.2, 14.7, 3.2, 37.5, 10.6, 16.72, 3
     @test [patient.procedure_code for patient in state.served_patients] ==
           [patient.procedure_code for patient in repeated.served_patients]
 end
+
+@testset "Event scheduling" begin
+    state = initialise(TEST_PARAMETERS)
+    patient = make_patient(1, 3.0, RandomNGs(7, TEST_PARAMETERS))
+    nurse = state.nurses[1]
+    first_scheduled = StartService(3.0, patient, nurse)
+    second_scheduled = StartService(3.0, patient, nurse)
+    schedule!(state, first_scheduled, 3.0)
+    schedule!(state, second_scheduled, 3.0)
+
+    first_time, first_event = take_next_event!(state)
+    second_time, second_event = take_next_event!(state)
+
+    @test first_time == 3.0
+    @test second_time == 3.0
+    @test first_event === first_scheduled
+    @test second_event === second_scheduled
+    @test first_event.time == 3.0
+    @test second_event.time == 3.0
+end
+
+@testset "Replicate seeds" begin
+    seeds = replicate_seeds(900, 3)
+    @test seeds == [900, 901, 902]
+    @test length(unique(seeds)) == 3
+    @test_throws ArgumentError replicate_seeds(900, 0)
+end
+
+@testset "Parameter validation" begin
+    nonfinite = Parameters(1, NaN, 0.2, 14.7, 3.2, 37.5, 10.6, 16.72, 3.2, 22.27, 5.0, 1)
+    nonpositive_stddev = Parameters(1, 2.0, 0.0, 14.7, 3.2, 37.5, 10.6, 16.72, 3.2, 22.27, 5.0, 1)
+
+    @test_throws ArgumentError run_simulation(nonfinite, 1)
+    @test_throws ArgumentError run_simulation(nonpositive_stddev, 1)
+end
+
+@testset "Analytical invariants" begin
+    empty_parameters = Parameters(0, 2.0, 0.2, 14.7, 3.2, 37.5, 10.6, 16.72, 3.2, 22.27, 5.0, 2)
+    empty_state = run_simulation(empty_parameters, 42)
+    @test empty_state.n_arrived == 0
+    @test empty_state.n_served == 0
+    @test isempty(empty_state.served_patients)
+
+    state = run_simulation(TEST_PARAMETERS, 42)
+    @test state.n_arrived - state.n_served == 0
+    for nurse_id in unique(patient.served_by_nurse for patient in state.served_patients)
+        assignments = sort(
+            [patient for patient in state.served_patients if patient.served_by_nurse == nurse_id],
+            by = patient -> patient.service_start,
+        )
+        @test all(
+            assignments[index].service_start >= assignments[index - 1].departure_time
+            for index in 2:length(assignments)
+        )
+    end
+end

vcc.jl

@@ -6,6 +6,16 @@ using Dates
 using StatsBase
 abstract type Event end
 abstract type PatientEvent <: Event end
+
+struct ScheduledEvent
+    event::Event
+    sequence::Int
+end
+
+function Base.isless(a::ScheduledEvent, b::ScheduledEvent)
+    a.event.time == b.event.time ? a.sequence < b.sequence : a.event.time < b.event.time
+end
+
 # Entities
 mutable struct Patient
     ID::Int
@@ -40,7 +50,8 @@ mutable struct State
     n_arrived::Int
     n_served::Int
     t::Float64
-    event_queue::PriorityQueue{Float64, Event, Base.Order.ForwardOrdering}
+    event_queue::PriorityQueue{Int, ScheduledEvent, Base.Order.ForwardOrdering}
+    next_event_sequence::Int
     patients::Vector{Patient}
     nurses::Vector{Nurse}
     served_patients::Vector{Patient}  
@@ -115,6 +126,33 @@ mutable struct Parameters
     std_service_time_carecoordination::Float64
     n_nurses::Int
 end
+
+function replicate_seeds(base_seed::Int, n_replicates::Int)
+    n_replicates > 0 || throw(ArgumentError("n_replicates must be positive"))
+    return [Base.checked_add(base_seed, offset) for offset in 0:(n_replicates - 1)]
+end
+
+function validate_parameters(P::Parameters)
+    P.n >= 0 || throw(ArgumentError("n must be non-negative"))
+    P.n_nurses > 0 || throw(ArgumentError("n_nurses must be positive"))
+    for name in (
+        :mean_interarrival_time,
+        :std_interarrival_time,
+        :mean_service_time_review,
+        :std_service_time_review,
+        :mean_service_time_assessment,
+        :std_service_time_assessment,
+        :mean_service_time_nursing,
+        :std_service_time_nursing,
+        :mean_service_time_carecoordination,
+        :std_service_time_carecoordination,
+    )
+        value = getfield(P, name)
+        isfinite(value) || throw(ArgumentError("$name must be finite"))
+        value > 0 || throw(ArgumentError("$name must be positive"))
+    end
+end
+
 function Base.isless(a::Event, b::Event)
     return a.time < b.time
 end
@@ -191,20 +229,22 @@ function initialise(P::Parameters)
     n_arrived = 0
     n_served = 0
     t = 0.0
-    event_queue = PriorityQueue{Float64, Event}()
+    event_queue = PriorityQueue{Int, ScheduledEvent}()
     patients = Vector{Patient}()
     nurses = [Nurse(i, false, 0.0) for i in 1:P.n_nurses]
     served_patients = Vector{Patient}()
-    return State(n_arrived, n_served, t, event_queue, patients, nurses, served_patients)
+    return State(n_arrived, n_served, t, event_queue, 0, patients, nurses, served_patients)
 end
 
 function schedule!(state::State, event::Event, time::Float64)
-    # Ensure unique timestamps by adding a small epsilon if a time collision occurs
-    adjusted_time = time
-    while haskey(state.event_queue, adjusted_time)
-        adjusted_time += 1e-6  
-    end
-    enqueue!(state.event_queue, adjusted_time, event)
+    event.time == time || throw(ArgumentError("event time must match scheduled time"))
+    state.next_event_sequence += 1
+    enqueue!(state.event_queue, state.next_event_sequence, ScheduledEvent(event, state.next_event_sequence))
+end
+
+function take_next_event!(state::State)
+    _, scheduled = dequeue_pair!(state.event_queue)
+    return scheduled.event.time, scheduled.event
 end
 
 function make_patient(id::Int, arrival_time::Float64, rngs::RandomNGs)
@@ -289,8 +329,7 @@ function handle_nurse_busy(state::State, event::NurseBusy, rngs::RandomNGs, P::P
     end
 end
 function run_simulation(P::Parameters, rng_seed::Int = 1234)
-    P.n >= 0 || throw(ArgumentError("n must be non-negative"))
-    P.n_nurses > 0 || throw(ArgumentError("n_nurses must be positive"))
+    validate_parameters(P)
     state = initialise(P)
     P.n == 0 && return state
     rngs = RandomNGs(rng_seed, P)
@@ -299,7 +338,7 @@ function run_simulation(P::Parameters, rng_seed::Int = 1234)
     schedule!(state, Arrival(0.0, first_patient), 0.0)
 
     while state.n_served < P.n && !isempty(state.event_queue)
-        current_time, event = dequeue_pair!(state.event_queue)
+        current_time, event = take_next_event!(state)
         state.t = current_time
         if event isa Arrival
             handle_arrival(state, event, rngs, P)

vccrun.jl

@@ -9,6 +9,7 @@ NURSE_NUMBERS = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
 
 # number of patients list
 PATIENT_NUMBERS = collect(1:100)
+BASE_SEED = 1234
 
 
 
@@ -21,7 +22,7 @@ for nurse_number in NURSE_NUMBERS
         total_times = Float64[]  # Store total times for each simulation run
 
         # Run 100 simulations for each nurse-patient combination
-        for i in 1:100
+        for seed in replicate_seeds(BASE_SEED, 100)
             # Define parameters for the simulation
             P = Parameters(
                 patient_number,  # Number of patients
@@ -39,7 +40,7 @@ for nurse_number in NURSE_NUMBERS
             )
 
             # Run the simulation
-            final_state = run_simulation(P)
+            final_state = run_simulation(P, seed)
             push!(total_times, final_state.t)  # Store total simulation time
         end