diff --git a/pyproject.toml b/pyproject.toml index c155cf77..a1e4bfed 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,6 +49,7 @@ dependencies = [ "pydoe>=0.9.7", "scipy>=1.16.3", "seaborn>=0.13.2", + "simpy>=4.1.2", "structlog>=25.5.0", ] diff --git a/simopt/input_models.py b/simopt/input_models.py index 066f2adc..151abf48 100644 --- a/simopt/input_models.py +++ b/simopt/input_models.py @@ -10,6 +10,7 @@ P = ParamSpec("P") R = TypeVar("R", covariant=True) +T = TypeVar("T") class InputModel(Protocol[P, R]): @@ -75,7 +76,7 @@ def random(self, alpha: float, beta: float) -> float: class WeightedChoice(InputModel): """Discrete weighted choice wrapper.""" - def random(self, population: Sequence[object], weights: Sequence[float]) -> object: + def random(self, population: Sequence[T], weights: Sequence[float]) -> T: """Sample an element from ``population`` according to ``weights``. Args: @@ -83,7 +84,7 @@ def random(self, population: Sequence[object], weights: Sequence[float]) -> obje weights (Sequence[float]): Nonnegative weights for each item. Returns: - Any: A randomly selected element from ``population``. + T: A randomly selected element from ``population``. """ # Calculate cumulative weights cum_weights = list(itertools.accumulate(weights)) diff --git a/simopt/models/ambulance.py b/simopt/models/ambulance.py index 3cfb5d72..9801590b 100644 --- a/simopt/models/ambulance.py +++ b/simopt/models/ambulance.py @@ -2,9 +2,11 @@ from __future__ import annotations -from typing import Annotated, ClassVar, Self +from collections.abc import Generator +from typing import Annotated, ClassVar, Self, cast import numpy as np +import simpy from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a @@ -19,9 +21,6 @@ from simopt.input_models import Beta, Exp from simopt.utils import override -AVAILABLE = 0 -BUSY = 1 - class AmbulanceConfig(BaseModel): """Configuration for the Ambulance simulation model.""" @@ -189,46 +188,6 @@ def replicate(self) -> tuple[dict, dict]: beta_x_model = self.beta_x_model beta_y_model = self.beta_y_model - def next_arrival(curr: float) -> list: - """Generate the next arrival event. - - - Arrival inter-time: exponential. - - Call location (x, y): drawn from scaled Beta distributions. - """ - # Draw Beta-based coordinates in [0, SQUARE_WIDTH] - x_coord = beta_x_model.random(alpha_x, beta_x) * sqaure_width - y_coord = beta_y_model.random(alpha_y, beta_y) * sqaure_width - - # Generate times using Exp input model - interarrival_time = arrival_time_model.random(1.0 / mean_interval) - service_time = scene_time_model.random(1.0 / mean_scene_time) - - return [ - curr + interarrival_time, # interarrival time - 1, # event type: arrival - x_coord, # x from Beta - y_coord, # y from Beta - service_time, # scene time - ] - - # ------------------------------ - # Event list: - # For type 0: end: [time, 0, 0, 0, 0] - # For type 1: arrival: [time, 1, x_coord, y_coord, service_time] - # For type 2: service completion: [time, 2, assigned_amb_index, 0, 0] - # ------------------------------ - event_list = [] - current_time = 0.0 - - # Schedule termination and first arrival - event_list.append([sim_length, 0, 0, 0, 0]) - event_list.append(next_arrival(0)) - - # Ambulance state: [x, y, status] - ambs = np.array([[bx, by, AVAILABLE] for bx, by in bases]) - queued_calls = [] - active_calls = 0 - total_response_time = 0.0 num_calls = 0 grad_total = np.zeros((variable_base_count, 2)) @@ -237,101 +196,71 @@ def next_arrival(curr: float) -> list: # carry is used if the next call for this ambulance has to wait carry_next = np.zeros((variable_base_count, 2)) - # ------------------------------ - # Main event loop - # ------------------------------ - while event_list: - event_list.sort(key=lambda e: e[0]) - event = event_list.pop(0) - current_time = event[0] - etype = event[1] - - if etype == 0: - # End - break - - if etype == 1: - # Arrival - active_calls += 1 - if active_calls > len(bases): - queued_calls.append(event) + env = simpy.Environment() + available_ambulances = simpy.FilterStore(env, capacity=n_ambulances) + for i in range(n_ambulances): + available_ambulances.put(i) + + def call( + arrival_time: float, + x_coord: float, + y_coord: float, + service_time: float, + ) -> Generator[simpy.Event, object, None]: + nonlocal num_calls, total_response_time + queued = not available_ambulances.items + + if queued: + i = cast(int, (yield available_ambulances.get())) + else: + times = [ + ( + np.sum(np.abs(np.array(bases[i]) - [x_coord, y_coord])) / amb_speed + if i in available_ambulances.items + else float("inf") + ) + for i in range(n_ambulances) + ] + i = int(np.argmin(times)) + yield available_ambulances.get(lambda ambulance: ambulance == i) + + travel = np.sum(np.abs(np.array(bases[i]) - [x_coord, y_coord])) / amb_speed + queue_delay = env.now - arrival_time + total_response_time += travel + queue_delay + num_calls += 1 + + if ( + i >= variable_base_start_index + and i - variable_base_start_index < variable_base_count + ): + j = i - variable_base_start_index + dx = np.sign(bases[i][0] - x_coord) / amb_speed + dy = np.sign(bases[i][1] - y_coord) / amb_speed + dd = np.array([dx, dy]) + if queued: + grad_total[j] += carry_next[j] + dd + carry_next[j] = carry_next[j] + 2.0 * dd else: - # Find nearest available ambulance - times = [ - ( - np.sum(np.abs(amb[:2] - event[2:4])) / amb_speed - if amb[2] == AVAILABLE - else float("inf") - ) - for amb in ambs - ] - i = int(np.argmin(times)) - ambs[i, 2] = BUSY - response_time = times[i] - total_response_time += response_time - num_calls += 1 - - # Gradient update - # no waiting: Response time (R) = Driving time (D) - if ( - i >= variable_base_start_index - and i - variable_base_start_index < variable_base_count - ): - j = i - variable_base_start_index - # Compute travel gradient dD wrt base position - dx = np.sign(ambs[i, 0] - event[2]) / amb_speed - dy = np.sign(ambs[i, 1] - event[3]) / amb_speed - dd = np.array([dx, dy]) - # No waiting time - grad_total[j] += dd - # Set carry for next potential queued call for this ambulance - # dW(i)+2dD(i) = 0 + 2dD - carry_next[j] = 2.0 * dd - - done_time = current_time + 2 * response_time + event[4] - event_list.append([done_time, 2, i, 0, 0]) - - event_list.append(next_arrival(current_time)) - - elif etype == 2: - # service completion - i = int(event[2]) - ambs[i, 2] = AVAILABLE - active_calls -= 1 - - if queued_calls: - # dispatch first queued call - qevent = queued_calls.pop(0) - ambs[i, 2] = BUSY # mark ambulance as busy again - - travel = np.sum(np.abs(ambs[i, 0:2] - qevent[2:4])) / amb_speed - queue_delay = current_time - qevent[0] - total_response_time += travel + queue_delay - num_calls += 1 - - # Gradient update (queued: R = W + D) - # If ambulance is from a variable base: - if ( - i >= variable_base_start_index - and i - variable_base_start_index < variable_base_count - ): - j = i - variable_base_start_index - # Compute travel gradient dD wrt base position - dx = np.sign(ambs[i, 0] - qevent[2]) / amb_speed - dy = np.sign(ambs[i, 1] - qevent[3]) / amb_speed - dd = np.array([dx, dy]) - # For queued calls: - # Response time (R) = Waiting time (W) + Driving time (D) - grad_total[j] += carry_next[j] + dd - # Update carry for the next call assigned to this ambulance - # new carry = old carry + 2 * dD - carry_next[j] = carry_next[j] + 2.0 * dd - - done_time = current_time + 2 * travel + qevent[4] - event_list.append([done_time, 2, i, 0, 0]) - - # if no queue, we do not need to change carry_next here - # the next non-queued arrival will overwrite it with 2*dD + grad_total[j] += dd + carry_next[j] = 2.0 * dd + + yield env.timeout(2 * travel) + yield env.timeout(service_time) + yield available_ambulances.put(i) + + def call_arrivals() -> Generator[simpy.Event, object, None]: + while True: + # Draw Beta-based coordinates in [0, SQUARE_WIDTH] + x_coord = beta_x_model.random(alpha_x, beta_x) * sqaure_width + y_coord = beta_y_model.random(alpha_y, beta_y) * sqaure_width + interarrival_time = arrival_time_model.random(1.0 / mean_interval) + service_time = scene_time_model.random(1.0 / mean_scene_time) + + yield env.timeout(interarrival_time) + env.process(call(env.now, x_coord, y_coord, service_time)) + + env.process(call_arrivals()) + env.run(until=sim_length) if num_calls: avg_time = total_response_time / num_calls diff --git a/simopt/models/amusementpark.py b/simopt/models/amusementpark.py index 0d7c7cca..23f7c747 100644 --- a/simopt/models/amusementpark.py +++ b/simopt/models/amusementpark.py @@ -2,9 +2,10 @@ from __future__ import annotations -import math as math -from typing import Annotated, ClassVar, Final, Self, cast +from collections.abc import Generator +from typing import Annotated, ClassVar, Final, Self +import simpy from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a @@ -19,8 +20,6 @@ from simopt.input_models import Exp, Gamma, WeightedChoice from simopt.models._ext import patch_model -INF = float("inf") - # Default values for the model PARK_CAPACITY: Final[int] = 350 NUM_ATTRACTIONS: Final[int] = 7 @@ -304,28 +303,6 @@ def replicate(self) -> tuple[dict[str, float | list[float]], dict]: - dict: Gradients of the performance measures with respect to model factors. """ # noqa: E501 - def set_completion(i: int, new_time: float) -> None: - """Set the completion time for an attraction. - - Updates the minimum completion time and index if necessary. - This function doesn't offer much (if any) performance gain with small - numbers of attractions, but with larger numbers it is significantly faster. - - Args: - i (int): The index of the attraction. - new_time (float): The new completion time for the attraction. - """ - nonlocal min_completion_time, min_completion_index - completion_times[i] = new_time - - if new_time < min_completion_time: - min_completion_time = new_time - min_completion_index = i - elif i == min_completion_index: - # Grab the min index and time with one scanning pass - min_completion_time = min(completion_times) - min_completion_index = completion_times.index(min_completion_time) - # Keep local copies of factors to prevent excessive lookups num_attractions: int = self.factors["number_attractions"] arrival_gammas: list[int] = self.factors["arrival_gammas"] @@ -340,12 +317,6 @@ def set_completion(i: int, new_time: float) -> None: attraction_range = range(num_attractions) destination_range = range(num_attractions + 1) depart_idx = destination_range[-1] - # initialize lists of each attraction's next completion time - completion_times: list[float] = [INF] * num_attractions - min_completion_time = INF - min_completion_index = -1 - # initialize actual queues. - queues: list[int] = [0] * num_attractions # create external arrival probabilities for each attraction. arrival_prob_sum: float = float(sum(arrival_gammas)) @@ -353,113 +324,90 @@ def set_completion(i: int, new_time: float) -> None: arrival_gammas[i] / arrival_prob_sum for i in attraction_range ] - # Initiate clock variables for statistics tracking and event handling. - clock: float = 0.0 - previous_clock: float = 0.0 - next_arrival = self.arrival_model.random(arrival_prob_sum) - # Initialize quantities to track: total_visitors: int = 0 total_departed: int = 0 # initialize time average and utilization quantities. - in_system: int = 0 time_average: float = 0.0 cumulative_util: list[float] = [0.0] * num_attractions + previous_clock: float = 0.0 + + env = simpy.Environment() + queues = [simpy.Store(env, capacity=max(1, queue_capacities[i])) for i in attraction_range] + busy = [False] * num_attractions - # Run simulation over time horizon. - while clock < time_open: - # Count number of tourists on attractions and in queues. - riders: int = 0 - delta_time: float = clock - previous_clock + def update_statistics() -> None: + nonlocal previous_clock, time_average + delta_time = env.now - previous_clock for i in attraction_range: - if not math.isinf(completion_times[i]): - riders += 1 + if busy[i]: cumulative_util[i] += delta_time - in_system = sum(queues) + riders - time_average += in_system * (delta_time) - - previous_clock = clock - if next_arrival < min_completion_time: - # Next event is external tourist arrival. - total_visitors += 1 - # Select attraction. - attraction_selection = cast( - int, - self.attraction_model.random( - attraction_range, - arrival_probabilities, - ), - ) - # Check if attraction is currently available. - # If available, arrive at that attraction. Otherwise check queue. - if math.isinf(completion_times[attraction_selection]): - # Generate completion time if attraction available. - completion_time = next_arrival + self.service_models[ - attraction_selection - ].random( - alpha=erlang_shape[attraction_selection], - beta=erlang_scale[attraction_selection], - ) - set_completion(attraction_selection, completion_time) - # If unavailable, check if current queue is less than capacity. - # If queue is not full, join queue. - elif queues[attraction_selection] < queue_capacities[attraction_selection]: - queues[attraction_selection] += 1 - # If queue is full, leave park + 1. - else: - total_departed += 1 - # Use superposition of Poisson processes to generate next arrival time. - next_arrival += self.arrival_model.random(arrival_prob_sum) + in_system = sum(len(queue.items) for queue in queues) + sum(busy) + time_average += in_system * delta_time + previous_clock = env.now + + def admit(attraction: int) -> None: + nonlocal total_departed + if not busy[attraction]: + busy[attraction] = True + queues[attraction].put(None) + elif len(queues[attraction].items) < queue_capacities[attraction]: + queues[attraction].put(None) else: - # Next event is the completion of an attraction. - # Identify finished attraction. - finished_attraction = completion_times.index(min_completion_time) - # Check if there is a queue for that attraction. - # If so then start new completion time and subtract 1 from queue. - alpha = erlang_shape[finished_attraction] - beta = erlang_scale[finished_attraction] - if queues[finished_attraction] > 0: - completion_time = min_completion_time + self.service_models[ - finished_attraction - ].random( - alpha=alpha, - beta=beta, - ) - set_completion(finished_attraction, completion_time) - queues[finished_attraction] -= 1 - else: # If attraction queue is empty, set next completion to infinity. - set_completion(finished_attraction, INF) - # Check if that person will leave the park. - next_destination = cast( - int, - self.destination_model.random( + total_departed += 1 + + def attraction_server( + finished_attraction: int, + ) -> Generator[simpy.Event, object, None]: + while True: + yield queues[finished_attraction].get() + service_time = self.service_models[finished_attraction].random( + alpha=erlang_shape[finished_attraction], + beta=erlang_scale[finished_attraction], + ) + while True: + yield env.timeout(service_time) + update_statistics() + + has_queued_visitor = bool(queues[finished_attraction].items) + if has_queued_visitor: + yield queues[finished_attraction].get() + service_time = self.service_models[finished_attraction].random( + alpha=erlang_shape[finished_attraction], + beta=erlang_scale[finished_attraction], + ) + else: + busy[finished_attraction] = False + + next_destination = self.destination_model.random( destination_range, transition_probabilities[finished_attraction] + [depart_probabilities[finished_attraction]], - ), + ) + if next_destination != depart_idx: + admit(next_destination) + + if not has_queued_visitor: + break + + def external_arrivals() -> Generator[simpy.Event, object, None]: + nonlocal total_visitors + while True: + interarrival_time = self.arrival_model.random(arrival_prob_sum) + yield env.timeout(interarrival_time) + update_statistics() + total_visitors += 1 + + attraction_selection = self.attraction_model.random( + attraction_range, arrival_probabilities ) + admit(attraction_selection) - # Check if tourist leaves park. - if next_destination != depart_idx: - # Check if attraction is currently available. - # If available, arrive at that attraction. Otherwise check queue. - if math.isinf(completion_times[next_destination]): - # Generate completion time if attraction available. - completion_time = min_completion_time + self.service_models[ - next_destination - ].random(alpha, beta) - set_completion(next_destination, completion_time) - # If unavailable, check if current queue is less than capacity. - # If queue is not full, join queue. - elif queues[next_destination] < queue_capacities[next_destination]: - queues[next_destination] += 1 - # If queue is full, leave park + 1. - else: - total_departed += 1 - # End of while loop. - # Check if any attractions are available. - clock = min(next_arrival, min_completion_time) - # End of simulation. + for i in attraction_range: + env.process(attraction_server(i)) + env.process(external_arrivals()) + env.run(until=time_open) + update_statistics() # Calculate overall percent utilization calculation for each attraction. cumulative_util = [cumulative_util[i] / time_open for i in attraction_range] diff --git a/simopt/models/chessmm.py b/simopt/models/chessmm.py index 34f0e965..bdc79332 100644 --- a/simopt/models/chessmm.py +++ b/simopt/models/chessmm.py @@ -2,10 +2,12 @@ from __future__ import annotations +from collections.abc import Generator from random import Random -from typing import Annotated, ClassVar, Final +from typing import Annotated, ClassVar, Final, cast import numpy as np +import simpy from pydantic import BaseModel, Field from scipy import special @@ -19,7 +21,7 @@ StochasticConstraint, VariableType, ) -from simopt.input_models import InputModel, Poisson +from simopt.input_models import Exp, InputModel MEAN_ELO: Final[int] = 1200 MAX_ALLOWABLE_DIFF: Final[int] = 150 @@ -147,7 +149,7 @@ def __init__(self, fixed_factors: dict | None = None) -> None: super().__init__(fixed_factors) self.elo_model = EloInputModel() - self.arrival_model = Poisson() + self.arrival_model = Exp() def before_replicate(self, rng_list: list[MRG32k3a]) -> None: self.elo_model.set_rng(rng_list[0]) @@ -176,37 +178,51 @@ def replicate(self) -> tuple[dict, dict]: allowable_diff = self.factors["allowable_diff"] poisson_rate = self.factors["poisson_rate"] - # Generate Elo ratings (normal distribution). - player_ratings = [ - self.elo_model.random(elo_mean, elo_sd, elo_min, elo_max) for _ in num_players_range - ] - - # Generate interarrival times (Poisson distribution). - interarrival_times = [self.arrival_model.random(poisson_rate) for _ in num_players_range] - # Initialize statistics. # Incoming players are initialized with a wait time of 0. wait_times = np.zeros(num_players) - waiting_players = [] + env = simpy.Environment() + waiting_players = simpy.FilterStore(env, capacity=num_players) total_diff = 0 # TODO: make this do something elo_diffs = [] - # Simulate arrival and matching and players. - for interarrival_time, player_rating in zip( - interarrival_times, player_ratings, strict=False - ): - # Try to match the player - for i, waiting_rating in enumerate(waiting_players): - diff = abs(player_rating - waiting_rating) - if diff <= allowable_diff: - total_diff += diff - elo_diffs.append(diff) - waiting_players.pop(i) - break - wait_times[i] += interarrival_time - # If break did not execute, then the player was not matched. - else: - waiting_players.append(player_rating) + def player_arrivals() -> Generator[simpy.Event, object, None]: + nonlocal total_diff + for player_idx in num_players_range: + # Generate the player's Elo rating and interarrival time. + player_rating = self.elo_model.random(elo_mean, elo_sd, elo_min, elo_max) + interarrival_time = self.arrival_model.random(poisson_rate) + yield env.timeout(interarrival_time) + + # Try to match the player + for waiting_player in waiting_players.items: + waiting_rating = waiting_player[0] + diff = abs(player_rating - waiting_rating) + if diff <= allowable_diff: + total_diff += diff + elo_diffs.append(diff) + matched_player = cast( + tuple[float, int, float], + ( + yield waiting_players.get( + lambda player, incoming_rating=player_rating: ( + abs(incoming_rating - player[0]) <= allowable_diff + ) + ) + ), + ) + wait_times[matched_player[1]] = env.now - matched_player[2] + break + # If break did not execute, then the player was not matched. + else: + yield waiting_players.put((player_rating, player_idx, env.now)) + + env.process(player_arrivals()) + env.run() + + # Players still in the pool have waited through the end of the replication. + for _, player_idx, arrival_time in waiting_players.items: + wait_times[player_idx] = env.now - arrival_time # If there weren't any matches, the elo_diffs list will be empty. avg_diff = np.mean(elo_diffs) if elo_diffs else np.nan diff --git a/simopt/models/hotel.py b/simopt/models/hotel.py index 3619e5c4..9bd0b170 100644 --- a/simopt/models/hotel.py +++ b/simopt/models/hotel.py @@ -2,10 +2,11 @@ from __future__ import annotations -import heapq +from collections.abc import Generator from typing import Annotated, ClassVar, Self import numpy as np +import simpy from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a @@ -304,35 +305,33 @@ def replicate(self) -> tuple[dict, dict]: ] ) - # Initialize arrival times - arrival = [-time_before + arr_time[i, 0] for i in range(num_products)] - a_idx = np.ones(num_products, dtype=int) # Next interarrival index - # Precompute resource conflict matrix (bool) conflicts = (product_incidence.T @ product_incidence) >= 1 - # Min-heap for tracking next arrival events (arrival_time, product_idx) - heap = [(arrival[i], i) for i in range(num_products) if arrival[i] <= time_limit[i]] - heapq.heapify(heap) - - while heap: - current_time, product_idx = heapq.heappop(heap) - if current_time > run_length: - break - if booking_limits[product_idx] > 0: - rate = rack_rate if product_idx % 2 == 0 else discount_rate - total_revenue += rate * np.sum(product_incidence[:, product_idx]) - for i in range(num_products): - if conflicts[product_idx, i] and booking_limits[i] > 0: - booking_limits[i] -= 1 - - # Schedule next arrival for this product - next_idx = a_idx[product_idx] - if next_idx < arr_bound: - next_time = current_time + arr_time[product_idx, next_idx] - a_idx[product_idx] += 1 - if next_time <= time_limit[product_idx]: - heapq.heappush(heap, (next_time, product_idx)) + env = simpy.Environment(initial_time=-time_before) + booking_inventory = [ + simpy.Container(env, capacity=max(1, limit), init=limit) for limit in booking_limits + ] + + def booking_stream(product_idx: int) -> Generator[simpy.Event, object, None]: + nonlocal total_revenue + for interarrival_time in arr_time[product_idx]: + next_time = env.now + interarrival_time + if next_time > time_limit[product_idx] or next_time > run_length: + return + yield env.timeout(interarrival_time) + + if booking_inventory[product_idx].level > 0: + rate = rack_rate if product_idx % 2 == 0 else discount_rate + total_revenue += rate * np.sum(product_incidence[:, product_idx]) + for i in range(num_products): + inventory = booking_inventory[i] + if conflicts[product_idx, i] and inventory.level > 0: + inventory.get(1) + + for product_idx in range(num_products): + env.process(booking_stream(product_idx)) + env.run() # Compose responses and gradients. responses = {"revenue": total_revenue} diff --git a/simopt/models/mm1queue.py b/simopt/models/mm1queue.py index 77e64344..1cac8368 100644 --- a/simopt/models/mm1queue.py +++ b/simopt/models/mm1queue.py @@ -2,10 +2,12 @@ from __future__ import annotations +from collections.abc import Generator from enum import IntEnum from typing import Annotated, ClassVar import numpy as np +import simpy from pydantic import BaseModel, Field from mrg32k3a.mrg32k3a import MRG32k3a @@ -198,46 +200,53 @@ class Col(IntEnum): cust_mat = np.zeros((total, 10)) cust_mat[:, Col.ARR] = np.cumsum(arrival_times) cust_mat[:, Col.SVC] = service_times - # Input entries for first customer's queueing experience. - first_cust = cust_mat[0] - first_cust[Col.DONE] = first_cust[Col.ARR] + first_cust[Col.SVC] - first_cust[Col.SOJ] = first_cust[Col.SVC] - # first_cust[Col.WAIT] = 0 - # cfirst_cust[Col.IN_SYS] = 0 - first_cust[Col.G_SOJ_MU] = -first_cust[Col.SVC] / mu_floor - # first_cust[Col.G_WAIT_MU] = 0 - # first_cust[Col.G_SOJ_LAM] = 0 - # first_cust[Col.G_WAIT_LAM] = 0 - # Fill in entries for remaining customers' experiences. - for i in range(1, total): - # Views into the customer matrix. - # NOT copies, so be careful! - curr_cust = cust_mat[i] - prev_cust = cust_mat[i - 1] - arrival = curr_cust[Col.ARR] - prev_departure = prev_cust[Col.DONE] + env = simpy.Environment() + server = simpy.Resource(env, capacity=1) - # Completion time - curr_cust[Col.DONE] = max(arrival, prev_departure) + curr_cust[Col.SVC] - # Sojourn and waiting times - curr_cust[Col.SOJ] = curr_cust[Col.DONE] - arrival - curr_cust[Col.WAIT] = curr_cust[Col.SOJ] - curr_cust[Col.SVC] + def customer(i: int) -> Generator[simpy.Event, object, None]: + curr_cust = cust_mat[i] + arrival = curr_cust[Col.ARR] + yield env.timeout(arrival) # Number in system at arrival - lookback = int(prev_cust[Col.IN_SYS]) + 1 - arrivals_in_window = cust_mat[i - lookback : i, Col.DONE] - curr_cust[Col.IN_SYS] = np.count_nonzero(arrivals_in_window > arrival) + curr_cust[Col.IN_SYS] = server.count + len(server.queue) - # Gradients w.r.t mu - n_in_sys = int(curr_cust[Col.IN_SYS]) - grad_range = cust_mat[i - n_in_sys : i + 1, Col.SVC] - curr_cust[Col.G_SOJ_MU] = -np.sum(grad_range) / mu_floor - curr_cust[Col.G_WAIT_MU] = -np.sum(grad_range[:-1]) / mu_floor + with server.request() as request: + yield request + yield env.timeout(curr_cust[Col.SVC]) + + curr_cust[Col.DONE] = env.now + curr_cust[Col.SOJ] = curr_cust[Col.DONE] - arrival + curr_cust[Col.WAIT] = curr_cust[Col.SOJ] - curr_cust[Col.SVC] # Gradients w.r.t lambda # cust_mat[i, 8] = 0.0 # cust_mat[i, 9] = 0.0 + + for i in range(total): + env.process(customer(i)) + env.run() + + # Calculate IPA gradients with respect to mu. A customer's completion-time + # gradient carries through the entire busy period, including customers who + # have already departed. Below the service-rate floor, the simulated service + # times do not depend on mu; at the floor, use the zero-valued left derivative. + if mu > epsilon: + prev_done_grad_mu = 0.0 + for i in range(total): + arrival = cust_mat[i, Col.ARR] + service = cust_mat[i, Col.SVC] + busy = i > 0 and cust_mat[i - 1, Col.DONE] > arrival + + grad_wait_mu = prev_done_grad_mu if busy else 0.0 + grad_service_mu = -service / mu + grad_sojourn_mu = grad_wait_mu + grad_service_mu + + cust_mat[i, Col.G_WAIT_MU] = grad_wait_mu + cust_mat[i, Col.G_SOJ_MU] = grad_sojourn_mu + prev_done_grad_mu = grad_sojourn_mu + cust_mat_warmup = cust_mat[warmup:] # Compute average sojourn time and its gradient. mean_sojourn_time = np.mean(cust_mat_warmup[:, Col.SOJ]) diff --git a/simopt/models/network.py b/simopt/models/network.py index 8fb6e09d..c5e4f78f 100644 --- a/simopt/models/network.py +++ b/simopt/models/network.py @@ -2,11 +2,13 @@ from __future__ import annotations +from collections.abc import Generator from enum import IntEnum from random import Random from typing import Annotated, ClassVar, Final, Self import numpy as np +import simpy from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a @@ -288,31 +290,32 @@ class Col(IntEnum): message_mat[:, Col.ARR] = np.cumsum(arrival_times) message_mat[:, Col.ROUTE] = network_routes message_mat[:, Col.SVC] = service_times - # Fill in entries for remaining messages' metrics. - # Create a list recording the index of the last customer sent to each network. - # Starting with -1, indicating no one is in line. + # Fill in entries for messages' metrics. routes = message_mat[:, Col.ROUTE].astype(int) arrival = message_mat[:, Col.ARR] service = message_mat[:, Col.SVC] - # Initialize completion time tracking per network - last_in_line = [-1] * n_networks + env = simpy.Environment() + networks = [simpy.Resource(env, capacity=1) for _ in range(n_networks)] - for i in range(total_arrivals): + def message(i: int) -> Generator[simpy.Event, object, None]: net = routes[i] arr_i = arrival[i] svc_i = service[i] + curr_message = message_mat[i] - if last_in_line[net] == -1: - done = arr_i + svc_i - else: - done = max(arr_i, message_mat[last_in_line[net], Col.DONE]) + svc_i + yield env.timeout(arr_i) + with networks[net].request() as request: + yield request + curr_message[Col.WAIT] = env.now - arr_i + yield env.timeout(svc_i) - curr_message = message_mat[i] - curr_message[Col.DONE] = done - curr_message[Col.SOJ] = done - arr_i - curr_message[Col.WAIT] = curr_message[Col.SOJ] - svc_i - last_in_line[net] = i + curr_message[Col.DONE] = env.now + curr_message[Col.SOJ] = curr_message[Col.DONE] - arr_i + + for i in range(total_arrivals): + env.process(message(i)) + env.run() # Vectorized cost computations after SOJ is known message_mat[:, Col.PROC_COST] = np.array(cost_process)[routes] @@ -321,12 +324,6 @@ class Col(IntEnum): message_mat[:, Col.PROC_COST] + message_mat[:, Col.TIME_COST] ) - routes = message_mat[:, Col.ROUTE].astype(int) - message_mat[:, Col.PROC_COST] = np.array(cost_process)[routes] - message_mat[:, Col.TIME_COST] = np.array(cost_time)[routes] * message_mat[:, Col.SOJ] - message_mat[:, Col.TOTAL_COST] = ( - message_mat[:, Col.PROC_COST] + message_mat[:, Col.TIME_COST] - ) # Compute total costs for the simulation run. total_cost = np.sum(message_mat[:, Col.TOTAL_COST]) responses = {"total_cost": total_cost} diff --git a/simopt/models/tableallocation.py b/simopt/models/tableallocation.py index 8371e6e0..f51e2cc1 100644 --- a/simopt/models/tableallocation.py +++ b/simopt/models/tableallocation.py @@ -2,12 +2,11 @@ from __future__ import annotations -import bisect -import itertools -from collections.abc import Sequence -from typing import Annotated, ClassVar, Self, cast +from collections.abc import Generator +from typing import Annotated, ClassVar, Self import numpy as np +import simpy from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a @@ -196,31 +195,6 @@ def replicate(self) -> tuple[dict, dict]: each response. """ - def fast_weighted_choice( - population: Sequence[int], weights: Sequence[float], rng: MRG32k3a - ) -> int: - """Select a single element from a population based on weights. - - Designed to be faster than `random.choices()` when only one element - is needed. - - Args: - population (Sequence[int]): The population to select from. - weights (Sequence[float]): The weights for each element in the - population. - rng (MRG32k3a): The random number generator to use for selection. - - Returns: - int: The selected element from the population. - """ - # Calculate cumulative weights - cum_weights = list(itertools.accumulate(weights)) - # Generate a value somewhere between 0 and the sum of weights - x = rng.random() * cum_weights[-1] - # Find the index of the first cumulative weight that is >= x - # Return the corresponding element from the population - return population[bisect.bisect(cum_weights, x)] - num_tables = self.factors["num_tables"] # TODO: figure out how floats are getting into the num_tables list num_tables = [int(n) for n in num_tables] @@ -232,9 +206,6 @@ def fast_weighted_choice( table_revenue = self.factors["table_revenue"] # Track total revenue. total_rev = 0 - # Track table availability. - # (i,j) is the time that jth table of size i becomes available. - table_avail = np.zeros((4, max(num_tables))) # Generate total number of arrivals in the period n_arrivals = self.arrival_number_model.random(round(n_hours * sum(f_lambda))) # Generate arrival times in minutes @@ -245,15 +216,19 @@ def fast_weighted_choice( found = np.zeros(n_arrivals) # Precompute options for group sizes. group_size_options = list(range(1, max_table_cap + 1)) - # Pass through all arrivals of groups to the restaurants. - for n in range(n_arrivals): + env = simpy.Environment() + tables = [ + [simpy.Resource(env, capacity=1) for _ in range(num_tables[k])] + for k in range(len(num_tables)) + ] + + def group(n: int) -> Generator[simpy.Event, object, None]: + nonlocal total_rev + yield env.timeout(arrival_times[n]) + # Determine group size. - group_size = cast( - int, - self.group_size_model.random( - population=group_size_options, - weights=f_lambda, - ), + group_size = self.group_size_model.random( + population=group_size_options, weights=f_lambda ) # Find smallest table size to start search. @@ -262,28 +237,36 @@ def fast_weighted_choice( table_size_idx += 1 # Find smallest available table. - def find_table(table_size_idx: int, n: int) -> tuple[int, int] | None: - for k in range(table_size_idx, len(num_tables)): - for j in range(num_tables[k]): - # Check if table is currently available. - if table_avail[k, j] < arrival_times[n]: - return k, j - # Return None if no table is available. - return None - - result = find_table(table_size_idx, n) + result = next( + ( + (k, j) + for k in range(table_size_idx, len(num_tables)) + for j, table in enumerate(tables[k]) + if table.count == 0 + ), + None, + ) # If no table is available, move on to next group. if result is None: - continue + return k, j = result - # Mark group as seated. - found[n] = 1 - # Sample service time. - service_time = self.service_time_model.random(1 / service_time_means[group_size - 1]) - # Update table availability. - table_avail[k, j] += service_time - # Update revenue. - total_rev += table_revenue[group_size - 1] + with tables[k][j].request() as request: + yield request + # Mark group as seated. + found[n] = 1 + # Sample service time. + service_time = self.service_time_model.random( + 1 / service_time_means[group_size - 1] + ) + # Update revenue. + total_rev += table_revenue[group_size - 1] + # Hold the table for the full sampled service duration. + yield env.timeout(service_time) + + # Pass through all arrivals of groups to the restaurants. + for n in range(n_arrivals): + env.process(group(n)) + env.run() # Calculate responses from simulation data. responses = { "total_revenue": total_rev, diff --git a/test/expected_results/AMUSEMENTPARK1_RNDSRCH.pickle.zst b/test/expected_results/AMUSEMENTPARK1_RNDSRCH.pickle.zst index 22bc5342..8e53e10b 100644 Binary files a/test/expected_results/AMUSEMENTPARK1_RNDSRCH.pickle.zst and b/test/expected_results/AMUSEMENTPARK1_RNDSRCH.pickle.zst differ diff --git a/test/expected_results/CHESS1_FCSA.pickle.zst b/test/expected_results/CHESS1_FCSA.pickle.zst index 759f88a0..c22e94cb 100644 Binary files a/test/expected_results/CHESS1_FCSA.pickle.zst and b/test/expected_results/CHESS1_FCSA.pickle.zst differ diff --git a/test/expected_results/CHESS1_RNDSRCH.pickle.zst b/test/expected_results/CHESS1_RNDSRCH.pickle.zst index a9a9ffcb..b8053db0 100644 Binary files a/test/expected_results/CHESS1_RNDSRCH.pickle.zst and b/test/expected_results/CHESS1_RNDSRCH.pickle.zst differ diff --git a/test/expected_results/MM11_ADAM.pickle.zst b/test/expected_results/MM11_ADAM.pickle.zst index e75328b9..5ac1d115 100644 Binary files a/test/expected_results/MM11_ADAM.pickle.zst and b/test/expected_results/MM11_ADAM.pickle.zst differ diff --git a/test/expected_results/MM11_ALOE.pickle.zst b/test/expected_results/MM11_ALOE.pickle.zst index 8eeb1d76..40ea01e6 100644 Binary files a/test/expected_results/MM11_ALOE.pickle.zst and b/test/expected_results/MM11_ALOE.pickle.zst differ diff --git a/test/expected_results/MM11_ASTRODF_darwin_arm64.pickle.zst b/test/expected_results/MM11_ASTRODF_darwin_arm64.pickle.zst index eedaa8eb..51d2e639 100644 Binary files a/test/expected_results/MM11_ASTRODF_darwin_arm64.pickle.zst and b/test/expected_results/MM11_ASTRODF_darwin_arm64.pickle.zst differ diff --git a/test/expected_results/MM11_NELDMD.pickle.zst b/test/expected_results/MM11_NELDMD.pickle.zst index a23466ff..95fdecbd 100644 Binary files a/test/expected_results/MM11_NELDMD.pickle.zst and b/test/expected_results/MM11_NELDMD.pickle.zst differ diff --git a/test/expected_results/MM11_RNDSRCH.pickle.zst b/test/expected_results/MM11_RNDSRCH.pickle.zst index c025bc68..efbc942a 100644 Binary files a/test/expected_results/MM11_RNDSRCH.pickle.zst and b/test/expected_results/MM11_RNDSRCH.pickle.zst differ diff --git a/test/expected_results/MM11_SPSA.pickle.zst b/test/expected_results/MM11_SPSA.pickle.zst index dc6eed9c..4ba120f1 100644 Binary files a/test/expected_results/MM11_SPSA.pickle.zst and b/test/expected_results/MM11_SPSA.pickle.zst differ diff --git a/test/expected_results/MM11_STRONG.pickle.zst b/test/expected_results/MM11_STRONG.pickle.zst index f5574442..79f42d87 100644 Binary files a/test/expected_results/MM11_STRONG.pickle.zst and b/test/expected_results/MM11_STRONG.pickle.zst differ diff --git a/test/expected_results/TABLEALLOCATION1_RNDSRCH.pickle.zst b/test/expected_results/TABLEALLOCATION1_RNDSRCH.pickle.zst index 9f0b4da0..f8c70ad3 100644 Binary files a/test/expected_results/TABLEALLOCATION1_RNDSRCH.pickle.zst and b/test/expected_results/TABLEALLOCATION1_RNDSRCH.pickle.zst differ