diff --git a/README.md b/README.md index bb5ac22..bf0b45b 100644 --- a/README.md +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 291fa73..a20bb73 100644 --- a/test/runtests.jl +++ b/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 diff --git a/vcc.jl b/vcc.jl index 944a565..b729b8c 100644 --- a/vcc.jl +++ b/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) diff --git a/vccrun.jl b/vccrun.jl index 5969725..44b6cad 100644 --- a/vccrun.jl +++ b/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