8  Filtered solution spaces for outpatient appointment scheduling

8.1 Introduction

Outpatient appointment scheduling is a critical operational challenge in healthcare systems, striving to efficiently manage patient flow while balancing multiple, often conflicting, objectives. The core dilemma involves minimizing patient waiting times, reducing clinician idle time, and preventing clinic sessions from running excessively late (tardiness). Effectively addressing this requires robust scheduling policies that can adapt to inherent uncertainties, such as variable service times and patient no-shows.

Work by Kaandorp and Koole (2007) provides a comprehensive mathematical framework for tackling this problem. They proposed an objective function that considers a weighted average of expected patient waiting times, doctor idle time, and session tardiness, explicitly incorporating the possibility of patient no-shows. A key contribution of their research was the demonstration that, under certain conditions, their objective function exhibits multimodularity. This property is significant because it guarantees that a local search algorithm can converge to a globally optimal appointment schedule. Their model assumes discrete time intervals and exponentially distributed service times.

This document details a computational exploration and implementation inspired by the principles and challenges outlined in the work of Kaandorp and Koole. Our primary objective is to investigate the characteristics of optimal, or near-optimal, appointment schedules under various parameter settings. While the original paper focuses on the theoretical properties and a local search method, this study employs a direct enumeration and evaluation approach for a specific class of schedule structures.

Specifically, we implement a Python-based model to:

1. Generate a Service Time Distribution: Create a plausible distribution of patient service times with a target average, allowing for variability.
2. Construct Candidate Schedules: Systematically generate a variety of patient appointment schedules. The schedules explored here are characterized by a defined number of patients at the beginning (‘start’) and end (‘tail’) of the session, with the remaining patients distributed as evenly as possible across the central portion of the available time intervals. We will explore this fitered solution space.
3. Evaluate Schedule Performance: For each generated schedule, calculate an objective function value. This objective is a weighted sum, reflecting the trade-off between minimizing expected patient waiting time (EWT) and minimizing expected service provider underutilization or overtime (ESP). This evaluation incorporates the probability of patient no-shows (q_no_show) and a user-defined weight (w_weight) to adjust the relative importance of EWT and ESP.
4. Visualize and Analyze Results: Utilize interactive 3D visualizations to explore the solution space. This allows for the examination of how different parameters—such as the no-show probability (q), the objective function weight (w), and the total number of patients (N_patients)—influence the structure of optimal schedules and their corresponding objective values.

Through this computational study, we aim to gain deeper insights into the sensitivity of optimal schedules to key operational parameters and to visualize the complex trade-offs inherent in outpatient appointment systems. The findings from this enumerative approach can also serve as a benchmark and provide foundational understanding for the development or comparison of more sophisticated heuristic or exact optimization algorithms.

import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.linalg import hadamard
from scipy.optimize import minimize
from itertools import combinations
import random
import textwrap
from typing import List, Tuple, Dict, Optional, Union
import unittest
import time

from functions import (
    get_v_star,
    compute_convolutions,
    generate_weighted_list,
    calculate_objective_serv_time_lookup,
    local_search_w_timer
)

# ---- Global parameters ----
d_interval_len = 10  # Length of each discrete time interval
max_s_time = 30  # Maximum service time (discrete units)
l_target_avg_service_time = 15.0  # Target average service time
i_sorting_split = 10  # Index to split sorting of service time probabilities

# Parameters for fig_search2
N_patients_fig_search2 = 74
q_patients_fig_search2 = 0.1
w_patients_fig_search2 = 0.9
T_intervals_patients_fig_search2 = 80
# s_dist is assumed to be the globally computed service time distribution.
# d_interval_len is assumed to be the globally defined d_interval_len = 10.
# HEURISTIC_LIMIT = N_patients_fig_search2
ls_time_limit_fig_search2 = 60

# --- Helper: Generate evenly distributed patient schedule ---
def generate_evenly_distributed_schedule_intervals(n_patients: int, t_intervals: int) -> List[int]:
    """
    Generates a list representing a schedule, distributing 'n_patients'
    into 't_intervals' as evenly as possible with guaranteed placement.
    A '1' indicates a scheduled interval, and '0' an empty one.
    
    Key Behavior:
    - If n_patients == 0: Returns a list of all zeros of length 't_intervals'.
    - If t_intervals < n_patients (and n_patients > 0): Raises a ValueError.
    - If n_patients == t_intervals (and both > 0): Returns a list of all ones.
    - If n_patients < t_intervals (and both > 0): Places each patient at the 
      center of their allocated segment to ensure even distribution without collisions.
    
    **Guarantee**: This approach ensures exactly 'n_patients' will be scheduled
    with no duplicates or collisions.
    
    Args:
        n_patients (int): The number of patients to schedule.
                          Must be a non-negative integer.
        t_intervals (int): The total number of available time intervals.
                           Must be a non-negative integer.
    
    Returns:
        List[int]: A list of integers (0s and 1s) of length 't_intervals',
                   representing the schedule. '1' for a scheduled patient, '0' otherwise.
    
    Raises:
        ValueError: If 'n_patients' or 't_intervals' are not non-negative integers.
        ValueError: If 't_intervals' is less than 'n_patients' (and n_patients > 0),
                    as it's impossible to schedule more patients than available intervals.
    """
    # Validate inputs
    if not isinstance(n_patients, int) or n_patients < 0:
        raise ValueError("n_patients must be a non-negative integer.")
    if not isinstance(t_intervals, int) or t_intervals < 0:
        raise ValueError("t_intervals must be a non-negative integer.")
    
    # Handle n_patients == 0
    if n_patients == 0:
        return [0] * t_intervals
    
    # At this point, n_patients > 0
    if t_intervals < n_patients:
        raise ValueError(
            f"Cannot schedule {n_patients} patients in only {t_intervals} intervals. "
            "Not enough unique slots available."
        )
    elif n_patients == t_intervals:
        return [1] * t_intervals
    else:  # n_patients < t_intervals
        schedule = [0] * t_intervals
        
        # Place each patient at the center of their allocated segment
        # This guarantees no collisions and even distribution
        for i in range(n_patients):
            # Each patient gets a segment of size (t_intervals / n_patients)
            # Place patient at center: (i + 0.5) * segment_size
            position = int((i + 0.5) * t_intervals / n_patients)
            schedule[position] = 1
        
        return schedule

# --- Helper: Generate start and tail numbers ---
def generate_start_tail_distribution(N_patients: int, T_intervals: int) -> List[Tuple[int, int]]:
    """
    Generates valid combinations of 'start' and 'tail' patient counts for a schedule.

    'Start' refers to the number of patients scheduled consecutively at the
    beginning of the `T_intervals`. 'Tail' refers to those at the end.
    The function iterates through all possible numbers of patients at the start (from 0 to N_patients)
    and for each, all possible numbers of patients at the tail (from 0 up to remaining patients).
    A (start, tail) combination is considered valid if the remaining patients
    (N_patients - start - tail) can fit within the remaining intervals
    (T_intervals - 2, assuming start and tail each occupy one conceptual "block"
    or represent counts for the first and last actual interval/slot).

    Args:
        N_patients (int): The total number of patients to be scheduled.
                          Must be a positive integer for useful output.
        T_intervals (int): The total number of available time intervals.

    Returns:
        List[Tuple[int, int]]: A list of tuples, where each tuple is (start_patients, tail_patients).
                               Returns an empty list if `N_patients <= 0`.
    
    Constraint for validity:
        The number of patients to be scheduled in the "center" part of the schedule
        (N_patients - start - tail) must be less than or equal to the number of
        intervals available for the center part (T_intervals - 2). This assumes
        the 'start' and 'tail' numbers might represent single bulk appointments or
        counts for the very first and very last interval positions, leaving T-2 intervals
        for the "center" distribution.
    """
    if not isinstance(N_patients, int) or N_patients <= 0:
        return []
    if not isinstance(T_intervals, int) or T_intervals < 0: # Allow T_intervals = 0, 1 for edge cases
        # Consider raising error if T_intervals < 2 if the -2 logic is strict
        pass

    distribution = []
    for i in range(N_patients + 1):  # Number of patients at the start
        start = i
        for j in range(N_patients + 1 - i): # Number of patients at the tail
            tail = j
            # Number of patients to be scheduled in the center:
            center_patients = N_patients - start - tail
            # Number of intervals available for the center:
            # Assumes start and tail might refer to the first and last interval slot counts,
            # or conceptual blocks, leaving T_intervals - 2 for the middle part.
            center_intervals = T_intervals - 2
            
            if center_intervals < 0: # Not enough intervals for even start/tail representation
                if center_patients == 0: # Only if no center patients are needed
                     # This case (e.g., T_intervals=1, N_patients=1, start=1, tail=0) might be valid
                     # depending on interpretation. The condition below will handle it.
                     pass # Let the condition below decide
                else:
                    continue # Cannot schedule center patients if no center_intervals

            if center_patients <= center_intervals or (center_patients == 0 and center_intervals >= -1) : # Allow N-s-t=0 even if T-2 is -1 or 0
                distribution.append((start, tail))
    return distribution
  
def compute_schedules_and_objectives(
    q_no_show: float,
    w_weight: float,
    current_N_patients: int, # Added to make N_patients explicit
    current_T_intervals: int # Added to make T_intervals explicit
    ) -> Tuple[List[List[int]], List[float]]:
    """
    Computes various patient schedules and their corresponding objective function values.

    This function orchestrates the generation of schedules based on start/tail
    distributions and a central, evenly distributed portion. It then calculates
    an objective value for each schedule, which is a weighted sum of expected
    waiting time (EWT) and expected service provider overtime/idleness (ESP).

    Args:
        q_no_show (float): The probability of a patient not showing up.
        w_weight (float): The weighting factor for the objective function.
                          It balances EWT and ESP: `w_weight * EWT + (1 - w_weight) * ESP`.
        current_N_patients (int): The total number of patients for this computation.
        current_T_intervals (int): The total number of time intervals for this computation.

    Returns:
        Tuple[List[List[int]], List[float]]:
            - schedules (List[List[int]]): A list of generated schedules. Each schedule
              is a list of integers: `[start_count, c1, c2, ..., c(T-2), tail_count]`,
              where `c_i` is 0 or 1 for the center part. The sum of elements in
              a schedule should equal `current_N_patients`.
            - objective_values (List[float]): A list of objective function values
              corresponding to each schedule.

    Raises:
        ValueError: If `s_dist` (globally defined service time distribution) is None
                    at the time of calling `compute_convolutions` (indirectly, as it's a global).
        ValueError: If `generate_evenly_distributed_schedule_intervals` fails to produce
                    a center schedule (e.g., returns None, though its current implementation raises).
    
    Dependencies:
        - `s_dist` (global): Pre-computed service time distribution.
        - `d_interval_len` (global): Length of each discrete interval.
        - `compute_convolutions`: External function to get convolution dictionary.
        - `generate_start_tail_distribution`: Helper to get start/tail patient counts.
        - `generate_evenly_distributed_schedule_intervals`: Helper for the center part of schedule.
        - `calculate_objective_serv_time_lookup`: External function for EWT and ESP.
    """
    if s_dist is None:
        raise ValueError("s_dist (service time distribution) is None.")

    # `s_dist` from `generate_weighted_list` is np.ndarray.
    # `compute_convolutions` might expect a list.
    s_dist_list = s_dist.tolist() if isinstance(s_dist, np.ndarray) else s_dist

    convolutions_dict = compute_convolutions(s_dist_list, current_N_patients, q_no_show)

    # Generate possible (start, tail) combinations for patient counts
    # Note: generate_start_tail_distribution uses N_patients and T_intervals
    # implicitly. It's better to pass them as arguments if they can vary.
    # Here, we use current_N_patients and current_T_intervals.
    start_tail_options = generate_start_tail_distribution(current_N_patients, current_T_intervals)
    schedules = []
    objective_values = []

    for start_count, tail_count in start_tail_options:
        schedule_parts = [] # To build the full schedule: [start, center..., tail]

        # The 'start_count' itself represents the number of patients in the first conceptual slot/block.
        # In the context of `calculate_objective_serv_time_lookup`, this might be interpreted
        # as the count for the first element of the schedule list.
        schedule_parts.append(start_count)

        num_center_patients = current_N_patients - start_count - tail_count
        num_center_intervals = current_T_intervals - 2 # Intervals available for the center part

        if num_center_patients < 0: # Should not happen if generate_start_tail_distribution is correct
            # print(f"Warning: Negative center patients ({num_center_patients}) for start={start_count}, tail={tail_count}. Skipping.")
            continue
        
        if num_center_intervals < 0:
            if num_center_patients > 0:
                # print(f"Warning: Cannot schedule {num_center_patients} center patients in {num_center_intervals} intervals. Skipping.")
                continue
            # If num_center_patients is 0 and num_center_intervals is < 0, center_schedule should be empty.
            center_schedule_segment = []
        elif num_center_patients == 0: # No patients for the center part
             center_schedule_segment = [0] * num_center_intervals
        else:
            center_schedule_segment = generate_evenly_distributed_schedule_intervals(
                num_center_patients,
                num_center_intervals
            )
            # generate_evenly_distributed_schedule_intervals raises ValueError on failure
            # so no need to check for None if its contract is maintained.

        schedule_parts.extend(center_schedule_segment)
        schedule_parts.append(tail_count) # Add the tail count as the last element

        # Ensure the generated schedule has the correct total length current_T_intervals
        # The structure is [start_val] + [c1, ..., c_{T-2}] + [tail_val]
        # So, 1 (for start) + (T-2) (for center) + 1 (for tail) = T elements.
        if len(schedule_parts) != current_T_intervals:
            # This might happen if T_intervals < 2.
            # Example: T_intervals = 1. num_center_intervals = -1.
            # schedule_parts = [start_count] + [] + [tail_count] -> length 2, not 1.
            # The interpretation of 'start' and 'tail' and how they map to the schedule list
            # for `calculate_objective_serv_time_lookup` is critical here.
            # Assuming `calculate_objective_serv_time_lookup` expects a list of length `current_T_intervals`
            # where each element is a count of patients for that interval.
            # If T_intervals=1, schedule should be [N_patients_fig1].
            # If T_intervals=0, schedule should be [].
            # The current structure [start, ccc, tail] implies T_intervals >= 2.
            # Let's adjust if T_intervals < 2.
            if current_T_intervals == 1:
                if start_count + tail_count == current_N_patients and not center_schedule_segment: # Only start or tail contributes
                     # The logic of start/tail might imply distinct first/last intervals.
                     # For T_intervals=1, the schedule is just [N_patients_fig1].
                     # We need a single value for the schedule.
                     # This case is tricky with start/center/tail logic.
                     # A simpler schedule representation might be needed for T_intervals < 2.
                     # For now, let's assume `calculate_objective_serv_time_lookup` can handle it,
                     # or this case (T_intervals < 2) is not expected for this specific schedule structure.
                     # If `start_count` is meant to be schedule[0] and `tail_count` schedule[-1],
                     # then for T_intervals=1, start_count and tail_count would refer to the same slot.
                     # The current `generate_start_tail_distribution` allows (N,0) or (0,N) for T=1, N.
                     # Let's assume the schedule sent to objective should be of length T_intervals.
                     if current_N_patients == start_count and tail_count == 0:
                         final_schedule = [start_count]
                     elif current_N_patients == tail_count and start_count == 0:
                         final_schedule = [tail_count]
                     else: # Ambiguous or invalid state for T=1 with this start/tail split.
                         # print(f"Warning: Ambiguous schedule for T_intervals=1 with start={start_count}, tail={tail_count}. Skipping.")
                         continue
                else: # Mismatch or center patients exist where they shouldn't
                    # print(f"Warning: Patient count mismatch for T_intervals=1. Skipping.")
                    continue

            elif current_T_intervals == 0:
                if current_N_patients == 0:
                    final_schedule = []
                else: # Cannot schedule patients in 0 intervals
                    # print(f"Warning: Cannot schedule {current_N_patients} in 0 intervals. Skipping.")
                    continue
            else: # T_intervals >= 2, schedule_parts should be T_intervals long.
                  # If not, there's a logic error in how center_schedule_segment was made.
                  # This should ideally not happen if num_center_intervals is correct.
                  # print(f"Warning: Schedule length {len(schedule_parts)} != T_intervals {current_T_intervals}. Skipping.")
                  # print(f"  Details: N={current_N_patients}, T={current_T_intervals}, start={start_count}, tail={tail_count}, center_p={num_center_patients}, center_i={num_center_intervals}")
                  # print(f"  Center segment: {center_schedule_segment}")
                  continue
        else:
            final_schedule = schedule_parts

        # Sanity check: sum of patients in the schedule should match current_N_patients
        if sum(final_schedule) != current_N_patients:
            # print(f"Warning: Sum of patients in schedule {sum(final_schedule)} != N_patients {current_N_patients}. Schedule: {final_schedule}. Skipping.")
            continue

        schedules.append(final_schedule)

        # Calculate objective value for this schedule
        # `d_interval_len` is used from global scope.
        ewt, esp = calculate_objective_serv_time_lookup(final_schedule, d_interval_len, convolutions_dict)
        objective_value = w_weight * ewt + (1 - w_weight) * esp
        objective_values.append(objective_value)

    return schedules, objective_values

8.2 Experiment 1 - Exploration of filtered solution spaces with small instance

We will consider all possible combinations of start and tail patients, and then fill the center part with evenly distributed patients. The center part is defined as the patients that are not in the start or tail part. The start and tail parts are defined as the first and last elements of the schedule, respectively. The center part is filled with patients in such a way that the number of patients in the center part is as evenly distributed as possible across the available time intervals and no interval contains more than one patient at a time.

The maximum number of solutions in the filtered equals \((N + 2) * (N + 1) / 2\). This limits the search space to a manageable size, allowing for an exhaustive exploration of the solution space.

We will start by visualizing the service time distribution, which is used to compute the expected waiting time and expected service provider overtime/idleness.

# --- Parameters for fig1 ---
N_patients_fig1 = 24  # Number of patients for Figure 1
T_intervals_fig1 = 48  # Total number of time intervals
# --- Values for 3x3 subplot grid ---
q_values_fig1 = [0.1, 0.15, 0.2]
w_values_fig1 = [0.1, 0.5, 0.9]

s_dist = generate_weighted_list(max_s_time, l_target_avg_service_time, i_sorting_split)
if s_dist is None:
    raise ValueError("Failed to generate service time distribution (s_dist).")
mean_s_dist = np.dot(np.arange(len(s_dist)), s_dist)
  
# Plot s_dist to visualize the service time distribution
fig0 = go.Figure()
fig0.add_trace(
    go.Bar(
        x=list(range(len(s_dist))),
        y=s_dist,
        name='Service Time Distribution',
        marker_color='blue'
    )
)
fig0.update_layout(
    title=f'Service Time Distribution<sub></br>(Mean: {mean_s_dist:.2f})</sub>',
    xaxis_title='Service Time (Discrete Units)',
    yaxis_title='Probability',
    height=400,
    width=800,
    margin=dict(t=50, b=50, l=50, r=50)
)
fig0.show()

Then we will create a 3x3 subplot grid with 3D scenes for different values of q (no-show probability) and w (weighting factor). Each subplot will show the schedules and their corresponding objective values. The minimum objective value for each subplot will be highlighted, and the hover text will show the schedule and its corresponding objective value.

# Create subplots with 3D scenes
fig1 = make_subplots(
    rows=3, cols=3,
    specs=[[{'type': 'scatter3d'} for _ in range(3)] for _ in range(3)],
    subplot_titles=[f'q={q:.2f}, w={w:.1f}' for w in w_values_fig1 for q in q_values_fig1], # Placeholder titles
    horizontal_spacing=0.02,
    vertical_spacing=0.12
)

# Global color scale limits for consistency
global_z_min = float('inf')
global_z_max = float('-inf')

# First pass: compute all data and find global min/max for consistent coloring
all_data = {}
for i, w in enumerate(w_values_fig1):
    for j, q in enumerate(q_values_fig1):
        schedules, objective_values = compute_schedules_and_objectives(q, w, N_patients_fig1, T_intervals_fig1)
        all_data[(i, j)] = (schedules, objective_values)
        if objective_values: # Check if list is not empty
            global_z_min = min(global_z_min, min(objective_values))
            global_z_max = max(global_z_max, max(objective_values))

# Handle case where no objective values were found (e.g., all lists were empty)
if global_z_min == float('inf'): global_z_min = 0
if global_z_max == float('-inf'): global_z_max = 1


# Second pass: create plots
subplot_titles_with_min_info = [] # To store titles with min objective info
for i, w in enumerate(w_values_fig1):
    for j, q in enumerate(q_values_fig1):
        schedules, objective_values = all_data[(i, j)]

        current_title = f'q={q:.2f}, w={w:.1f}'
        if not schedules or not objective_values:
            subplot_titles_with_min_info.append(f'{current_title}<br><sub>No data</sub>')
            # Add an empty trace or skip plotting for this subplot
            fig.add_trace(go.Scatter3d(x=[], y=[], z=[]), row=i+1, col=j+1) # Add empty trace
            continue

        # Extract x, y, z coordinates for plotting
        # x: value of the first element of the schedule (start_count)
        # y: value of the last element of the schedule (tail_count)
        # z: objective value
        x_coords = [schedule[0] for schedule in schedules]
        y_coords = [schedule[-1] for schedule in schedules] # schedule can be shorter than T_intervals if T_intervals < 2
        # Ensure y_coords access is safe if schedule can be very short (e.g. length 1)
        y_coords = [(schedule[-1] if len(schedule)>0 else 0) for schedule in schedules]


        z_values = objective_values

        # Find minimum point
        min_obj_val = min(z_values)
        min_index = z_values.index(min_obj_val)
        # min_schedule_at_min_obj = schedules[min_index] # For debugging or detailed hover

        subplot_titles_with_min_info.append(f'{current_title}<br><sub>Min: {min_obj_val:.4f}</sub>')

        # Create hover text
        schedules_wrapped = [textwrap.wrap(str(schedule), width=75) for schedule in schedules]
        formatted_schedules_str = ["<br>".join(schedule_lines) for schedule_lines in schedules_wrapped]

        hover_texts = [
            f"Sch: {f_sch_str}, N_sch: {sum(sch_list)}<br>Obj: {obj_val:.4f}<br>"
            for sch_list, f_sch_str, obj_val in zip(schedules, formatted_schedules_str, z_values)
        ]

        # Separate normal and minimum points for distinct plotting
        normal_indices = [idx for idx in range(len(z_values)) if idx != min_index]
        normal_x = [x_coords[idx] for idx in normal_indices]
        normal_y = [y_coords[idx] for idx in normal_indices]
        normal_z = [z_values[idx] for idx in normal_indices]
        normal_hover = [hover_texts[idx] for idx in normal_indices]

        min_x = [x_coords[min_index]]
        min_y = [y_coords[min_index]]
        min_z_val = [z_values[min_index]] # Renamed to avoid conflict
        min_hover = [hover_texts[min_index]]

        # Add normal points scatter plot
        fig1.add_trace(
            go.Scatter3d(
                x=normal_x,
                y=normal_y,
                z=normal_z,
                mode='markers',
                marker=dict(
                    size=4,
                    color=normal_z,          # Color by Z value (objective)
                    colorscale='Viridis',
                    opacity=0.7,
                    cmin=global_z_min,       # Use global min for color scale
                    cmax=global_z_max,       # Use global max for color scale
                    colorbar=dict(title='Objective Value') if i == 0 and j == 2 else None, # Show colorbar once
                    showscale=(i == 0 and j == 2) # Control visibility of the color bar
                ),
                name=f'Data (q={q:.2f}, w={w:.1f})', # Legend name (if legend shown)
                hovertemplate='%{hovertext}<extra></extra>', # Custom hover text
                hovertext=normal_hover,
                showlegend=False
            ),
            row=i+1, col=j+1
        )

        # Add minimum point highlighted
        fig1.add_trace(
            go.Scatter3d(
                x=min_x,
                y=min_y,
                z=min_z_val,
                mode='markers',
                marker=dict(
                    size=6,
                    color='red',    # Distinct color for minimum point
                    opacity=1.0
                ),
                name=f'Minimum (q={q:.2f}, w={w:.1f})', # Legend name
                hovertemplate='%{hovertext}<extra></extra>', # Custom hover text
                hovertext=min_hover,
                showlegend=False
            ),
            row=i+1, col=j+1
        )

# Update subplot titles with the minimum objective information
for idx, title_text in enumerate(subplot_titles_with_min_info):
    if idx < len(fig1.layout.annotations): # Safety check
        fig1.layout.annotations[idx].text = title_text

# Update overall layout for the figure
fig1.update_layout(
    title=dict(
        text=f'3D Analysis: Schedule Optimization (N={N_patients_fig1})<br><sub>Rows: w (weighting factor), Columns: q (no-show probability)</sub>',
        x=0.5, # Center title
        y=0.98,
        font=dict(size=16)
    ),
    height=1200,
    width=1400,
    margin=dict(t=100, b=100, l=50, r=50),
    font=dict(size=10) # Default font size for other elements
)

# Update 3D scene properties for all subplots (axes titles, camera angle)
for r_idx_loop in range(1, 4): # 1-based for rows
    for c_idx_loop in range(1, 4): # 1-based for columns
        # Calculate scene name (e.g., 'scene', 'scene2', 'scene3', ..., 'scene9')
        scene_num = (r_idx_loop - 1) * 3 + c_idx_loop
        scene_name_key = f'scene{scene_num}' if scene_num > 1 else 'scene'

        fig1.update_layout(**{
            scene_name_key: dict(
                xaxis_title='Start Count (Schedule[0])',
                yaxis_title='Tail Count (Schedule[-1])',
                zaxis_title='Objective Value',
                camera=dict(eye=dict(x=-1.5, y=-1.5, z=1.5)), # Adjust camera view
                xaxis=dict(range=[0, None]), # Ensure x-axis starts at 0
                yaxis=dict(range=[0, None])  # Ensure y-axis starts at 0
            )
        })

fig1.show()

8.3 Experiment 2 - Exploration of filtered solution spaces with large instance

We repeat the same procedure as in Experiment 1, but with a larger instance. The number of patients is increased to 40, and the total number of time intervals is set to 48. We will explore the solution space for two different values of q (no-show probability) and two different values of w (weighting factor). The results will be visualized in a 2x2 subplot grid.

# --- Parameters for fig1 ---
N_patients_fig2 = 40  # Number of patients for Figure 1
T_intervals_fig2 = 48  # Total number of time intervals
# --- Values for 2x2 subplot grid ---
q_values_fig2 = [0.1, 0.2]
w_values_fig2 = [0.1, 0.9]

# --- Create 2x2 subplot grid for Figure 2 ---
fig2 = make_subplots(
    rows=len(w_values_fig2), cols=len(q_values_fig2),
    specs=[[{'type': 'scatter3d'} for _ in range(len(q_values_fig2))] for _ in range(len(w_values_fig2))],
    subplot_titles=[f'q={q_val:.2f}, w={w_val:.1f}' for w_val in w_values_fig2 for q_val in q_values_fig2], # Placeholder
    horizontal_spacing=0.02,
    vertical_spacing=0.12
)

# Global color scale limits for consistency in Figure 2
global_z_min_fig2 = float('inf')
global_z_max_fig2 = float('-inf')

# First pass: compute all data for Figure 2 and find global min/max for consistent coloring
all_data_fig2 = {}
for i, w in enumerate(w_values_fig2):
    for j, q in enumerate(q_values_fig2):
        schedules_fig2, objective_values_fig2 = compute_schedules_and_objectives(q, w, N_patients_fig2, T_intervals_fig2) # Use N_patients_fig2
        all_data_fig2[(i, j)] = (schedules_fig2, objective_values_fig2)
        if objective_values_fig2: # Ensure list is not empty
            global_z_min_fig2 = min(global_z_min_fig2, min(objective_values_fig2))
            global_z_max_fig2 = max(global_z_max_fig2, max(objective_values_fig2))
        else:
            print(f"Warning: No objective values for Figure 2 with q={q}, w={w}, N_patients={N_patients_fig2}")

if global_z_min_fig2 == float('inf'): global_z_min_fig2 = 0 # Default if no data
if global_z_max_fig2 == float('-inf'): global_z_max_fig2 = 1 # Default if no data


# Second pass: create plots for Figure 2
subplot_titles_with_min_fig2 = []
for i, w in enumerate(w_values_fig2):
    for j, q in enumerate(q_values_fig2):
        schedules, objective_values = all_data_fig2[(i, j)]

        current_title_fig2 = f'q={q:.2f}, w={w:.1f}'
        if not schedules or not objective_values:
            subplot_titles_with_min_fig2.append(f'{current_title_fig2}<br><sub>No data</sub>')
            fig2.add_trace(go.Scatter3d(x=[], y=[], z=[]), row=i+1, col=j+1) # Add empty trace
            continue

        x_coords_fig2 = [schedule[0] for schedule in schedules]
        y_coords_fig2 = [(schedule[-1] if len(schedule)>0 else 0) for schedule in schedules] # Safe access
        z_values_fig2 = objective_values

        min_obj_fig2 = min(z_values_fig2)
        min_index_fig2 = z_values_fig2.index(min_obj_fig2)

        subplot_titles_with_min_fig2.append(f'{current_title_fig2}<br><sub>Min: {min_obj_fig2:.4f}</sub>')

        schedules_wrapped_fig2 = [textwrap.wrap(str(schedule), width=75) for schedule in schedules]
        formatted_schedules_str_fig2 = ["<br>".join(sl) for sl in schedules_wrapped_fig2]

        hover_texts_fig2 = [
            f"Sch: {f_sch_str}, N_sch: {sum(sch_list)}<br>Obj: {obj_val:.4f}<br>"
            for sch_list, f_sch_str, obj_val in zip(schedules, formatted_schedules_str_fig2, z_values_fig2)
        ]

        normal_indices_fig2 = [idx for idx in range(len(z_values_fig2)) if idx != min_index_fig2]
        normal_x_fig2 = [x_coords_fig2[idx] for idx in normal_indices_fig2]
        normal_y_fig2 = [y_coords_fig2[idx] for idx in normal_indices_fig2]
        normal_z_fig2 = [z_values_fig2[idx] for idx in normal_indices_fig2]
        normal_hover_fig2 = [hover_texts_fig2[idx] for idx in normal_indices_fig2]

        min_x_fig2 = [x_coords_fig2[min_index_fig2]]
        min_y_fig2 = [y_coords_fig2[min_index_fig2]]
        min_z_val_fig2 = [z_values_fig2[min_index_fig2]]
        min_hover_fig2 = [hover_texts_fig2[min_index_fig2]]

        fig2.add_trace(
            go.Scatter3d(
                x=normal_x_fig2, y=normal_y_fig2, z=normal_z_fig2, mode='markers',
                marker=dict(
                    size=4, color=normal_z_fig2, colorscale='Viridis', opacity=0.7,
                    cmin=global_z_min_fig2, cmax=global_z_max_fig2,
                    colorbar=dict(title='Objective Value') if i == 0 and j == (len(q_values_fig2) - 1) else None,
                    showscale=(i == 0 and j == (len(q_values_fig2) - 1))
                ),
                name=f'Data (q={q:.2f}, w={w:.1f})',
                hovertemplate='%{hovertext}<extra></extra>', hovertext=normal_hover_fig2, showlegend=False
            ), row=i+1, col=j+1
        )

        fig2.add_trace(
            go.Scatter3d(
                x=min_x_fig2, y=min_y_fig2, z=min_z_val_fig2, mode='markers',
                marker=dict(size=6, color='red', opacity=1.0),
                name=f'Minimum (q={q:.2f}, w={w:.1f})',
                hovertemplate='%{hovertext}<extra></extra>', hovertext=min_hover_fig2, showlegend=False
            ), row=i+1, col=j+1
        )

# Update subplot titles for Figure 2
for idx, title_text in enumerate(subplot_titles_with_min_fig2):
    if idx < len(fig2.layout.annotations):
        fig2.layout.annotations[idx].text = title_text

# Update layout for Figure 2
fig2.update_layout(
    title=dict(
        text=f'<b>Figure 2:</b> N_patients={N_patients_fig2} - Schedule Optimization<br><sub>Rows: w (weight), Columns: q (no-show probability)</sub>',
        x=0.5, y=0.98, font=dict(size=16)
    ),
    height=1000, width=1200,
    margin=dict(t=100, b=100, l=50, r=50), font=dict(size=10)
)

# Update 3D scene properties for Figure 2
num_rows_fig2 = len(w_values_fig2)
num_cols_fig2 = len(q_values_fig2)
for r_loop_idx in range(num_rows_fig2):
    for c_loop_idx in range(num_cols_fig2):
        scene_idx_plotly = r_loop_idx * num_cols_fig2 + c_loop_idx + 1
        scene_name_key_fig2 = f'scene{scene_idx_plotly}' if scene_idx_plotly > 1 else 'scene'

        fig2.update_layout(**{
            scene_name_key_fig2: dict(
                xaxis_title='Start Count (Schedule[0])',
                yaxis_title='Tail Count (Schedule[-1])',
                zaxis_title='Objective Value',
                camera=dict(eye=dict(x=-1.5, y=-1.5, z=1.5)),
                xaxis=dict(range=[0, None]),
                yaxis=dict(range=[0, None])
            )
        })

fig2.show()

8.4 Experiment 3 - Projection of solution space

Examining the results from the previous experiments, we can observe a common pattern of some kind of convexity in multiple directions. The theoretical link between multimodularity and convexity is well-established; for instance, Altman, Gaujal, and Hordijk (2000) prove that multimodularity is a sufficient condition for integer convexity.

A function \(f\) is said to be integer convex if the following holds. For vectors \(x\) and \(v\) in \(\mathbb{Z}^T\), we have \[f(x + v) - f(x) \geq f(x) - f(x-v).\]

# --- Parameters for fig3 ---
N_patients_fig3 = 40
T_intervals_fig3 = 48
q_fig3 = 0.1
w_fig3 = 0.5
projection_level = 32

# Compute schedules and objectives for the new parameters
schedules, objective_values = compute_schedules_and_objectives(
    q_fig3, 
    w_fig3, 
    N_patients_fig3, 
    T_intervals_fig3 # Use the existing global T_intervals
)

# Create the 3D scatter plot
fig3 = go.Figure()

if not schedules or not objective_values:
    print(f"No data to plot for N={N_patients_fig3}, q={q_fig3}, w={w_fig3}")
    fig3.add_annotation(text="No data to plot for these parameters",
                              xref="paper", yref="paper",
                              x=0.5, y=0.5, showarrow=False,
                              font=dict(size=16))
    fig3.update_layout(
        title=f'3D Schedule Analysis: N={N_patients_fig3}, q={q_fig3:.2f}, w={w_fig3:.1f} - No Data',
        height=500, width=700
    )

else:
    x_coords_fig3 = [schedule[0] for schedule in schedules]
    y_coords_fig3 = [(schedule[-1] if len(schedule)>0 else 0) for schedule in schedules]
    z_values_fig3 = objective_values

    min_obj_val_fig3 = min(z_values_fig3)
    min_index_fig3 = z_values_fig3.index(min_obj_val_fig3)

    schedules_wrapped_fig3 = [textwrap.wrap(str(schedule), width=75) for schedule in schedules]
    formatted_schedules_str_fig3 = ["<br>".join(schedule_lines) for schedule_lines in schedules_wrapped_fig3]

    hover_texts_fig3 = [
        f"Sch: {f_sch_str}<br>N_sched: {sum(sch_list)}<br>Obj: {obj_val:.4f}<br>" +
        f"Start: {sch_list[0]}, Tail: {(sch_list[-1] if len(sch_list)>0 else 0)}"
        for sch_list, f_sch_str, obj_val in zip(schedules, formatted_schedules_str_fig3, z_values_fig3)
    ]

    # --- Add visual representation of the plane: start_count + tail_count = projection_level ---
    if z_values_fig3: 
        plane_z_min = np.min(z_values_fig3) 
        plane_z_max = np.max(z_values_fig3) 
        
        if plane_z_min == plane_z_max: # Add margin if all z-values are the same
            plane_z_min -= 0.5 
            plane_z_max += 0.5
        
        # Define vertices for the plane segment x + y = projection_level
        # Ensure plane vertices don't exceed the valid range [0, N_patients_fig3_case]
        max_coord = min(projection_level, N_patients_fig3)
        plane_mesh_x = [0, max_coord, max_coord, 0] 
        plane_mesh_y = [max_coord, 0, 0, max_coord] 
        plane_mesh_z = [plane_z_min, plane_z_min, plane_z_max, plane_z_max]

        fig3.add_trace(go.Mesh3d(
            x=plane_mesh_x,
            y=plane_mesh_y,
            z=plane_mesh_z,
            # Triangles using vertex indices: (0,1,2) and (0,2,3)
            i=[0, 0],  
            j=[1, 2],  
            k=[2, 3],  
            opacity=0.15,
            color='gray',
            name=f'Plane S+T={projection_level}',
            showlegend=True
        ))

    # --- Add normal data points ---
    normal_indices_fig3 = [idx for idx in range(len(z_values_fig3)) if idx != min_index_fig3]
    normal_x_fig3 = [x_coords_fig3[idx] for idx in normal_indices_fig3]
    normal_y_fig3 = [y_coords_fig3[idx] for idx in normal_indices_fig3]
    normal_z_fig3 = [z_values_fig3[idx] for idx in normal_indices_fig3]
    normal_hover_fig3 = [hover_texts_fig3[idx] for idx in normal_indices_fig3]

    fig3.add_trace(
        go.Scatter3d(
            x=normal_x_fig3, y=normal_y_fig3, z=normal_z_fig3,
            mode='markers',
            marker=dict(size=4, color=normal_z_fig3, colorscale='Viridis', opacity=0.7, colorbar=dict(title='Objective Value')),
            hovertemplate='%{hovertext}<extra></extra>',
            hovertext=normal_hover_fig3,
            name='Data points' 
        )
    )

    # --- Add minimum point ---
    min_x_fig3 = [x_coords_fig3[min_index_fig3]]
    min_y_fig3 = [y_coords_fig3[min_index_fig3]]
    min_z_val_fig3 = [z_values_fig3[min_index_fig3]]
    min_hover_fig3 = [hover_texts_fig3[min_index_fig3]]
    
    fig3.add_trace(
        go.Scatter3d(
            x=min_x_fig3, y=min_y_fig3, z=min_z_val_fig3,
            mode='markers',
            marker=dict(size=6, color='red', opacity=1.0),
            hovertemplate='%{hovertext}<extra></extra>',
            hovertext=min_hover_fig3,
            name=f'Minimum (Obj: {min_obj_val_fig3:.4f})' 
        )
    )

    # --- Add projection points for start_count + tail_count = projection_level ---
    proj_x = []
    proj_y = []
    proj_z = []
    proj_hover = []

    for idx, schedule in enumerate(schedules):
        start_count = schedule[0]
        tail_count = schedule[-1] if len(schedule) > 0 else 0
        
        if start_count + tail_count == projection_level:
            proj_x.append(start_count)
            proj_y.append(tail_count)
            proj_z.append(objective_values[idx])
            proj_hover.append(hover_texts_fig3[idx]) 

    if proj_x: 
        fig3.add_trace(
            go.Scatter3d(
                x=proj_x, y=proj_y, z=proj_z,
                mode='markers',
                marker=dict(
                    size=5, 
                    color='cyan', 
                    opacity=0.9,
                    symbol='diamond'
                ),
                hovertemplate='%{hovertext}<extra></extra>',
                hovertext=proj_hover,
                name=f'On Plane S+T={projection_level}' 
            )
        )

    fig3.update_layout(
        title=dict(
            text=f'3D Schedule Analysis: N={N_patients_fig3}, q={q_fig3:.2f}, w={w_fig3:.1f}<br><sub>Optimal Objective Value: {min_obj_val_fig3:.4f} | Projection: S+T={projection_level}</sub>',
            x=0.5, font=dict(size=16)
        ),
        scene=dict(
            xaxis_title='Start Count (Schedule[0])',
            yaxis_title='Tail Count (Schedule[-1])',
            zaxis_title='Objective Value',
            camera=dict(eye=dict(x=-1.5, y=-1.5, z=1.5)),
            xaxis=dict(range=[0, N_patients_fig3], autorange=False), 
            yaxis=dict(range=[0, N_patients_fig3], autorange=False)
        ),
        height=700, width=900,
        showlegend=True,
        legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01, bgcolor="rgba(255,255,255,0.5)") 
    )

fig3.show()

A projection of a solution space is a subset of the solution space that satisfies some additional constraints, such as a fixed sum of certain components. In this case, we are interested in projections where the sum of the start count and tail count is equal to a target value.

Let \(m \leq N\) be a target value for the sum of the start count and tail count. We can generate a set of projections where the start count and tail count sum to \(m\) by starting with a schedule that has \(m\) patients in the first interval, \(0\) patients in the last interval, and evenly distributed patients in the center intervals. Then we can iteratively add the vector \(v_0 = (-1, 0, \ldots, 0, 1)\) to this schedule until we’ve reached \(m\) iterations.

Because the function \(f\) is multimodular, the projections generated in this way will be integer convex.

def generate_schedules_w_equal_start_tail_sum(
    target_sum: int,
    N: int,
    T: int,
    d: int,
    convolutions: Dict[int, np.ndarray],
    w: float = 0.5
) -> List[Tuple[List[int], float]]:
    """
    Generate filtered solution space projections based on a target sum for start_count + tail_count.
    
    Parameters:
        target_sum (int): The target sum for start_count + tail_count.
        N (int): Total number of patients.
        T (int): Number of intervals
        
    Returns:
        projections: Filtered projections as tuple of (schedule, objective).
    """
    ## Validate inputs
    if not isinstance(target_sum, int) or target_sum < 0:
        raise ValueError("target_sum must be a non-negative integer.")
    if not isinstance(N, int) or N < 0:
        raise ValueError("N must be a non-negative integer.")
    if not isinstance(T, int) or T < 0:
        raise ValueError("T must be a non-negative integer.")
      
    ## Create v_star
    v_star_matrix = get_v_star(T)
    v_star_zero = v_star_matrix[0]  # Assuming we want the first row for T intervals
    
    ## Create center_schedule_segment
    num_center_patients = N - target_sum # Patients in the center part. Must be non-negative.
    num_center_intervals = T - 2 # Intervals available for the center part. Must be non-negative.
    if num_center_patients < 0 or num_center_intervals < 0:
        # print(f"Warning: Invalid center part with N={N}, T={T}, target_sum={target_sum}. Skipping.")
        return []  # No valid projections if center part is negative
    center_schedule_segment = generate_evenly_distributed_schedule_intervals(
                num_center_patients,
                num_center_intervals
            )
    # Create initial schedule
    initial_schedule = [target_sum] + center_schedule_segment + [0]  # Start with target_sum, center part, and tail 0
    ewt, esp = calculate_objective_serv_time_lookup(
        initial_schedule, d, convolutions
    )
    initial_objective_value = w * ewt + (1 - w) * esp
    obj_diff = None
    projections = [(np.array(initial_schedule), initial_objective_value, obj_diff)]
    for i in range(1, target_sum + 1):  # Iterate from 1 to target_sum
        new_schedule = initial_schedule + v_star_zero  # Add v_star_zero to each interval
        if np.any(new_schedule < 0):
            print(f"\nStopping at iteration {i}: Schedule became invalid with negative values.")
            print(f"Invalid schedule: {new_schedule[:5]}...")
            break # Exit the loop
        ewt, esp = calculate_objective_serv_time_lookup(
        new_schedule, d, convolutions
    )
        new_objective_value = w * ewt + (1 - w) * esp
        obj_diff = new_objective_value - initial_objective_value
        projections.append((new_schedule, new_objective_value, obj_diff))
        initial_schedule = new_schedule  # Update for next iteration
        initial_objective_value = new_objective_value  # Update objective value for next iteration
    
    return projections

convolutions_projections_fig3 = compute_convolutions(s_dist, N_patients_fig3, q_fig3)
projections_example = generate_schedules_w_equal_start_tail_sum(projection_level, N_patients_fig3, T_intervals_fig3, d_interval_len, convolutions_projections_fig3, w_fig3)

# Extract data for plotting
schedules = [p[0] for p in projections_example]
objectives = [p[1] for p in projections_example]
deltas = [p[2] for p in projections_example]
# Replace None in deltas with 0 for plotting
deltas[0] = 0

# Create the chart
fig4 = make_subplots(
    rows=2, cols=1,
    shared_xaxes=True,
    vertical_spacing=0.1,
    subplot_titles=("Objective Value per Schedule Iteration", "Change in Objective (Delta) per Iteration")
)

# Add Objective Value trace
fig4.add_trace(go.Scatter(
    x=list(range(len(objectives))),
    y=objectives,
    mode='lines+markers',
    name='Objective Value',
    hovertemplate='<b>Iteration %{x}</b><br>Objective: %{y:.2f}<br>Schedule: %{customdata}<extra></extra>',
    customdata=[str(s) for s in schedules],
    marker=dict(color='royalblue')
), row=1, col=1)

# Add Delta trace
fig4.add_trace(go.Bar(
    x=list(range(len(deltas))),
    y=deltas,
    name='Objective Delta',
    hovertemplate='<b>Iteration %{x}</b><br>Delta: %{y:.2f}<extra></extra>',
    marker=dict(color='lightcoral')
), row=2, col=1)

# Update layout
fig4.update_layout(
    title_text='Analysis of Scheduling Projections',
    height=600,
    showlegend=False,
    xaxis2_title='Schedule Iteration',
    yaxis_title='Objective Value',
    yaxis2_title='Objective Delta'
)

fig4.show()

8.5 Experiment 4 - Local search on a schedule with a fixed start and tail

In this experiment we will exploit the integer convexity of the filtered solution space to perform a heuristic search on a schedule with a fixed start and tail. We will start at the lowest possible sum of the start and tail counts. This initial schedule will be iteratively improved by exploring neighboring schedules. If we find a better schedule, we will update the current schedule and continue the search until no better schedules are found in the filtered solution space. The goal is to find the best schedule with a fixed start and tail that minimizes the objective function. This intermediate optimal schedule will be used as a starting point for an normal local search. The heuristic search and local search will be run will run for a fixed time limit.

As a benchmark we will run a normal local search on the same initial schedule for the same time limit. This will allow us to compare the performance of the standard local search with the heuristic search on the filtered solution space.

We will start with a small instance of 24 patients and 48 time intervals, with a no-show probability of 0.1 and a weighting factor of 0.5. The time limit for the local search will be set to 60 seconds.

# --- Parameters for fig3 ---
N_patients_fig5 = 24
T_intervals_fig5 = 48
q_fig5 = 0.1
w_fig5 = 0.5
convolutions_projections_fig5 = compute_convolutions(s_dist, N_patients_fig5, q_fig5)
ls_time_limit = 60
target_sum = max(0, N_patients_fig5 - T_intervals_fig5 + 2)  # Initialize start + tail sum. Start from max(0, N - (T - 2)) to ensure non-negative start
stop_sum = 30
lowest_objective_value = float('inf')  # Initialize to a very high value

while target_sum <= stop_sum:
  found_lower = False
  num_center_patients = N_patients_fig5 - target_sum # Patients in the center part
  num_center_intervals = T_intervals_fig5 - 2 # Intervals available for the center part

  if num_center_patients < 0: # Should not happen if generate_start_tail_distribution is correct
      # print(f"Warning: Negative center patients ({num_center_patients}) for start={start_count}, tail={tail_count}. Skipping.")
      continue
  
  if num_center_intervals < 0:
      if num_center_patients > 0:
          # print(f"Warning: Cannot schedule {num_center_patients} center patients in {num_center_intervals} intervals. Skipping.")
          continue
      # If num_center_patients is 0 and num_center_intervals is < 0, center_schedule should be empty.
      center_schedule_segment = []
  elif num_center_patients == 0: # No patients for the center part
       center_schedule_segment = [0] * num_center_intervals
  else:
      center_schedule_segment = generate_evenly_distributed_schedule_intervals(
          num_center_patients,
          num_center_intervals
      )
  # Try all possible start/tail combinations that sum to target_sum
  for s in range(target_sum + 1):
    schedule_parts = [s]  # Start with the start count
    t = target_sum - s
    schedule_parts.extend(center_schedule_segment)
    schedule_parts.append(t)
    ewt, esp = calculate_objective_serv_time_lookup(
        schedule_parts, d_interval_len, convolutions_projections_fig5
    )
    new_objective_value = w_fig5 * ewt + (1 - w_fig5) * esp
    if new_objective_value < lowest_objective_value:
        lowest_objective_value = new_objective_value
        best_schedule = schedule_parts
        print(f"New lowest objective value found: {lowest_objective_value} for target_sum={target_sum}, schedule={schedule_parts}. Sum schedule={sum(schedule_parts)}")
        found_lower = True
  if not found_lower:
    print(f"No lower objective value found for target_sum={target_sum}. Stopping search.\nBest schedule so far: {best_schedule} with objective value {lowest_objective_value}")
    break 
  target_sum += 1

v_star = get_v_star(T_intervals_fig5)
x_star, c_star, log = local_search_w_timer(best_schedule, d_interval_len, convolutions_projections_fig5, w_fig5, v_star, T_intervals_fig5, time_limit = ls_time_limit, process_start_time = time.time(), echo = False)
print(f"Final best schedule found: {x_star} with objective value {c_star}")

New lowest objective value found: 10.811224646234729 for target_sum=0, schedule=[0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]. Sum schedule=24
New lowest objective value found: 7.844480972416717 for target_sum=1, schedule=[1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]. Sum schedule=24
No lower objective value found for target_sum=2. Stopping search.
Best schedule so far: [1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0] with objective value 7.844480972416717
Final best schedule found: [1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
 0 1 0 1 0 1 0 1 0 1 0] with objective value 7.844480972416717

def find_optimal_schedule(
    N_patients: int,                    # Total number of patients to schedule
    T_intervals: int,                   # Total number of time intervals available
    d_interval_len: int,                # Duration of each interval (in minutes)
    w: float,                          # Weight parameter: w*EWT + (1-w)*ESP (0 ≤ w ≤ 1)
    convolutions_projections: Dict,     # Pre-computed service time distributions for objective calculation
    heuristic_stop_limit: int = 20,    # Maximum target sum to explore in heuristic search
    local_search_time_limit: int = 60, # Time limit (seconds) for local search phase
    start_with_heuristic: bool = True   # Whether to use heuristic initialization vs simple start
) -> Tuple[List[Dict], float]:
    """
    Finds an optimal patient schedule using a two-stage approach:
    1. Heuristic Search: Systematically explores schedules with different start/tail patterns
    2. Local Search: Refines the best heuristic solution using neighborhood exploration
    
    The function logs all intermediate best solutions for analysis and plotting.
    
    Returns:
        optimization_log: List of dict entries with 'schedule', 'cost', 'source', 'time_elapsed'
        total_time_taken: Total optimization time in seconds
    """
    # Initialize timing and logging
    total_start_time = time.time()
    optimization_log = []
    
    # === INITIALIZATION STRATEGY SELECTION ===
    if start_with_heuristic:
        print(f"🚀 Starting optimization for {N_patients} patients with Heuristic Search.")
        stage_two_text = "\n--- Stage 2: Local Search Refinement ---"
        
        # --- Stage 1: Heuristic Search ---
        print("\n--- Stage 1: Heuristic Search ---")
        lowest_objective_value = float('inf')
        
        # Calculate minimum start+tail sum needed to ensure feasible center distribution
        # Logic: If we have more patients than interior intervals, some must go to start/tail
        start_target_sum = max(0, N_patients - T_intervals + 2)
        num_center_intervals = T_intervals - 2  # Interior intervals (excluding first and last)
        
        # Systematic exploration of start+tail combinations
        for target_sum in range(start_target_sum, heuristic_stop_limit + 1):
            # Calculate how many patients go in the center intervals
            num_center_patients = N_patients - target_sum
            if num_center_patients < 0: 
                continue  # Skip if target_sum exceeds total patients
            
            # Generate evenly distributed schedule for center intervals
            center_schedule = generate_evenly_distributed_schedule_intervals(
                num_center_patients, num_center_intervals
            )
            
            # Try all possible start/tail combinations that sum to target_sum
            for s in range(target_sum + 1):  # s = patients in first interval
                t = target_sum - s           # t = patients in last interval
                
                # Construct full schedule: [start, center_intervals..., tail]
                candidate_schedule = [s, *center_schedule, t] if T_intervals > 1 else [N_patients]
                
                # Evaluate objective function: w*EWT + (1-w)*ESP
                ewt, esp = calculate_objective_serv_time_lookup(
                    candidate_schedule, d_interval_len, convolutions_projections
                )
                current_objective = w * ewt + (1 - w) * esp
                
                # Track new best solutions
                if current_objective < lowest_objective_value:
                    lowest_objective_value = current_objective
                    time_elapsed = time.time() - total_start_time
                    
                    # Log this improvement for analysis
                    log_entry = {
                        "schedule": candidate_schedule,
                        "cost": lowest_objective_value,
                        "source": "heuristic",
                        "time_elapsed": time_elapsed
                    }
                    optimization_log.append(log_entry)
        
        # Handle case where heuristic search fails
        if not optimization_log:
            print("\nHeuristic search could not find a valid initial schedule.")
            # Fallback: put all patients in first interval
            initial_schedule = [N_patients] + [0] * (T_intervals - 1)
        else:
            # Use the best schedule found by heuristic search
            initial_schedule = optimization_log[-1]['schedule']
        
        print(f"\n✅ Stage 1 Complete. Best Heuristic Schedule: {initial_schedule}")
    
    else:
        # === NON-HEURISTIC START: Use same initial candidate for fair comparison ===
        print(f"🚀 Starting optimization for {N_patients} patients with Local Search only.")
        stage_two_text = "\n--- Local Search Only ---"
        # Create the same starting schedule that heuristic would begin with
        if T_intervals > 1:
            # Start with minimum feasible start+tail sum
            start_target_sum = max(0, N_patients - T_intervals + 2)
            num_center_intervals = T_intervals - 2
            num_center_patients = N_patients - start_target_sum
            
            # Generate center distribution
            center_schedule = []
            if num_center_intervals > 0 and num_center_patients >= 0:
                center_schedule = generate_evenly_distributed_schedule_intervals(
                    num_center_patients, num_center_intervals
                )
            
            # Use conservative start: all target_sum patients in tail
            s = 0
            t = start_target_sum
            initial_schedule = [s, *center_schedule, t]
        else:
            # Single interval case
            initial_schedule = [N_patients]
        
        # Safety check for schedule validity
        if sum(initial_schedule) != N_patients:
            print(f"Warning: Calculated initial schedule {initial_schedule} sums to {sum(initial_schedule)}, not {N_patients}. Falling back to naive.")
            initial_schedule = [N_patients] + [0] * (T_intervals - 1)
            
        print(f"Using heuristic's first candidate as starting schedule: {initial_schedule}")
    
    # === Stage 2: Local Search Refinement ===
    print(stage_two_text)
    
    if local_search_time_limit > 0 and T_intervals > 1:
        # Perform neighborhood-based local search starting from initial_schedule
        _, _, local_search_log = local_search_w_timer(
            initial_schedule,              # Starting point
            d_interval_len,               # Interval duration
            convolutions_projections,     # Service time distributions
            w,                           # Objective weight
            v_star=get_v_star(T_intervals),  # Search neighborhood parameter
            size=T_intervals,            # Schedule size
            time_limit=local_search_time_limit,  # Time budget
            process_start_time=total_start_time, # For consistent timing
            echo=False                   # Suppress verbose output
        )
        
        # Label log entries based on initialization strategy
        source_label = "local_search" if start_with_heuristic else "local_search_from_start"
        for entry in local_search_log:
            entry["source"] = source_label
        
        # Combine heuristic and local search logs
        optimization_log.extend(local_search_log)
    
    # === COMPLETION AND SUMMARY ===
    total_time_taken = time.time() - total_start_time
    
    if optimization_log:
        best_schedule = optimization_log[-1]['schedule']
        best_cost = optimization_log[-1]['cost']
        print(f"\n🏁 Optimization complete in {total_time_taken:.2f} seconds.")
        print(f"Best schedule found: {best_schedule} with cost {best_cost:.4f}")
    else:
        print(f"\n⚠️ Optimization complete in {total_time_taken:.2f} seconds, but no valid schedule found.")
    
    return optimization_log, total_time_taken

# --- Plotting ---

# Run the optimization and get the log
# --- Run 1: Heuristic + Local Search ---
heuristic_log, total_time = find_optimal_schedule(
    N_patients=N_patients_fig5,
    T_intervals=T_intervals_fig5,
    d_interval_len=d_interval_len,
    w=w_fig5,
    convolutions_projections=convolutions_projections_fig5,
    heuristic_stop_limit=N_patients_fig5,
    local_search_time_limit=ls_time_limit
)

# --- Run 2: Local Search Only (for the same total duration) ---
ls_only_log, _ = find_optimal_schedule(
    N_patients=N_patients_fig5,
    T_intervals=T_intervals_fig5,
    d_interval_len=d_interval_len,
    w=w_fig5,
    convolutions_projections=convolutions_projections_fig5,
    heuristic_stop_limit=N_patients_fig5,
    local_search_time_limit=total_time, # Use same total time as the first run
    start_with_heuristic=False
)

# --- Combine logs and Plot ---
full_log = heuristic_log + ls_only_log

df = pd.DataFrame(full_log)
fig5 = go.Figure()
colors = {
    'heuristic': '#636EFA', 
    'local_search': '#00CC96', 
    'local_search_from_start': '#EF553B'
}

for source_name, group in df.groupby('source'):
    fig5.add_trace(go.Scatter(
        x=group['time_elapsed'], y=group['cost'],
        name=source_name, mode='lines+markers',
        line=dict(shape='hv', color=colors.get(source_name)),
        marker=dict(size=8, line=dict(width=1, color='DarkSlateGrey'))
    ))

fig5.update_layout(
    title=f"Optimization Strategy Comparison<br><sub>N = {N_patients_fig5}, T = {T_intervals_fig5}, q = {q_fig5} , w = {w_fig5}</sub>",
    xaxis_title="Time Elapsed (seconds)",
    yaxis_title="Objective Function Cost",
    legend_title_text="Search Method",
    template="plotly_white"
)

print("\nDisplaying convergence plot...")
fig5.show()

🚀 Starting optimization for 24 patients with Heuristic Search.

--- Stage 1: Heuristic Search ---

✅ Stage 1 Complete. Best Heuristic Schedule: [1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]

--- Stage 2: Local Search Refinement ---

🏁 Optimization complete in 60.12 seconds.
Best schedule found: [1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0] with cost 7.8445
🚀 Starting optimization for 24 patients with Local Search only.
Using heuristic's first candidate as starting schedule: [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]

--- Local Search Only ---

🏁 Optimization complete in 60.26 seconds.
Best schedule found: [1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0] with cost 7.8445

Displaying convergence plot...

Next we will try a slightly more extensive instance of 40 patients and 48 time intervals, with a no-show probability of 0.1 and a weighting factor of 0.5. The time limit for the local search will remain 60 seconds.

# --- Parameters for fig6 ---
N_patients_fig6 = 40
T_intervals_fig6 = 48
q_fig6 = 0.1
w_fig6 = 0.5
convolutions_projections_fig6 = compute_convolutions(s_dist, N_patients_fig6, q_fig6)
ls_time_limit = 60
target_sum = max(0, N_patients_fig6 - T_intervals_fig6 + 2)  # Initialize start + taol sum. Start from max(0, N - (T - 2)) to ensure non-negative start
stop_sum = 30
lowest_objective_value = float('inf')  # Initialize to a very high value

while target_sum <= stop_sum:
  found_lower = False
  num_center_patients = N_patients_fig6 - target_sum # Patients in the center part
  num_center_intervals = T_intervals_fig6 - 2 # Intervals available for the center part

  if num_center_patients < 0: # Should not happen if generate_start_tail_distribution is correct
      # print(f"Warning: Negative center patients ({num_center_patients}) for start={start_count}, tail={tail_count}. Skipping.")
      continue
  
  if num_center_intervals < 0:
      if num_center_patients > 0:
          # print(f"Warning: Cannot schedule {num_center_patients} center patients in {num_center_intervals} intervals. Skipping.")
          continue
      # If num_center_patients is 0 and num_center_intervals is < 0, center_schedule should be empty.
      center_schedule_segment = []
  elif num_center_patients == 0: # No patients for the center part
       center_schedule_segment = [0] * num_center_intervals
  else:
      center_schedule_segment = generate_evenly_distributed_schedule_intervals(
          num_center_patients,
          num_center_intervals
      )
  # Try all possible start/tail combinations that sum to target_sum
  for s in range(target_sum + 1):
    schedule_parts = [s]  # Start with the start count
    t = target_sum - s
    schedule_parts.extend(center_schedule_segment)
    schedule_parts.append(t)
    ewt, esp = calculate_objective_serv_time_lookup(
        schedule_parts, d_interval_len, convolutions_projections_fig6
    )
    new_objective_value = w_fig6 * ewt + (1 - w_fig6) * esp
    if new_objective_value < lowest_objective_value:
        lowest_objective_value = new_objective_value
        best_schedule = schedule_parts
        print(f"New lowest objective value found: {lowest_objective_value} for target_sum={target_sum}, schedule={schedule_parts}. Sum schedule={sum(schedule_parts)}")
        found_lower = True
  if not found_lower:
    print(f"No lower objective value found for target_sum={target_sum}. Stopping search.\nBest schedule so far: {best_schedule} with objective value {lowest_objective_value}")
    break 
  target_sum += 1

v_star = get_v_star(T_intervals_fig6)
x_star, c_star, log = local_search_w_timer(best_schedule, d_interval_len, convolutions_projections_fig6, w_fig6, v_star, T_intervals_fig6, time_limit = ls_time_limit, process_start_time = time.time(), echo = False)
print(f"Final best schedule found: {x_star} with objective value {c_star}")

New lowest objective value found: 558.3822593408088 for target_sum=0, schedule=[0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0]. Sum schedule=40
New lowest objective value found: 486.9037987872014 for target_sum=1, schedule=[0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1]. Sum schedule=40
New lowest objective value found: 486.34316835298864 for target_sum=1, schedule=[1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0]. Sum schedule=40
New lowest objective value found: 440.5059525302765 for target_sum=2, schedule=[0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 2]. Sum schedule=40
New lowest objective value found: 431.15407955776806 for target_sum=2, schedule=[1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1]. Sum schedule=40
New lowest objective value found: 393.07048086845356 for target_sum=3, schedule=[0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 3]. Sum schedule=40
New lowest objective value found: 378.1340038597886 for target_sum=3, schedule=[1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 2]. Sum schedule=40
New lowest objective value found: 358.0830843579592 for target_sum=4, schedule=[0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 4]. Sum schedule=40
New lowest objective value found: 337.3703314986522 for target_sum=4, schedule=[1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 3]. Sum schedule=40
New lowest objective value found: 334.7136889411036 for target_sum=5, schedule=[0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 5]. Sum schedule=40
New lowest objective value found: 308.11957789236874 for target_sum=5, schedule=[1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 4]. Sum schedule=40
New lowest objective value found: 301.247615363644 for target_sum=6, schedule=[1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 5]. Sum schedule=40
New lowest objective value found: 300.0751292172615 for target_sum=7, schedule=[1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 6]. Sum schedule=40
No lower objective value found for target_sum=8. Stopping search.
Best schedule so far: [1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 6] with objective value 300.0751292172615
Final best schedule found: [1 1 1 1 0 1 1 1 0 1 1 0 1 1 1 0 1 1 0 1 1 0 1 1 1 0 1 1 0 1 1 1 0 1 1 0 1
 1 1 0 1 1 1 0 1 1 1 5] with objective value 291.1896404050326

# --- Plotting ---

# Run the optimization and get the log
# --- Run 1: Heuristic + Local Search ---
heuristic_log, total_time = find_optimal_schedule(
    N_patients=N_patients_fig6,
    T_intervals=T_intervals_fig6,
    d_interval_len=d_interval_len,
    w=w_fig6,
    convolutions_projections=convolutions_projections_fig6,
    heuristic_stop_limit=N_patients_fig6,
    local_search_time_limit=ls_time_limit
)

# --- Run 2: Local Search Only (for the same total duration) ---

ls_only_log, _ = find_optimal_schedule(
    N_patients=N_patients_fig6,
    T_intervals=T_intervals_fig6,
    d_interval_len=d_interval_len,
    w=w_fig6,
    convolutions_projections=convolutions_projections_fig6,
    heuristic_stop_limit=N_patients_fig6,
    local_search_time_limit=total_time, # Use same total time as the first run
    start_with_heuristic=False
)

# --- Combine logs and Plot ---
full_log = heuristic_log + ls_only_log

df = pd.DataFrame(full_log)
fig6 = go.Figure()
colors = {
    'heuristic': '#636EFA', 
    'local_search': '#00CC96', 
    'local_search_from_start': '#EF553B'
}

for source_name, group in df.groupby('source'):
    fig6.add_trace(go.Scatter(
        x=group['time_elapsed'], y=group['cost'],
        name=source_name, mode='lines+markers',
        line=dict(shape='hv', color=colors.get(source_name)),
        marker=dict(size=8, line=dict(width=1, color='DarkSlateGrey'))
    ))

fig6.update_layout(
    title=f"Optimization Strategy Comparison<sub>N = {N_patients_fig6}, T = {T_intervals_fig6}, q = {q_fig6} , w = {w_fig6}</sub>",
    xaxis_title="Time Elapsed (seconds)",
    yaxis_title="Objective Function Cost",
    legend_title_text="Search Method",
    template="plotly_white"
)

print("\nDisplaying convergence plot...")
fig6.show()

🚀 Starting optimization for 40 patients with Heuristic Search.

--- Stage 1: Heuristic Search ---

✅ Stage 1 Complete. Best Heuristic Schedule: [1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 6]

--- Stage 2: Local Search Refinement ---

🏁 Optimization complete in 60.18 seconds.
Best schedule found: [1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 5] with cost 291.1896
🚀 Starting optimization for 40 patients with Local Search only.
Using heuristic's first candidate as starting schedule: [0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0]

--- Local Search Only ---

🏁 Optimization complete in 60.37 seconds.
Best schedule found: [1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 4] with cost 297.2457

Displaying convergence plot...

Finally we will examine a large instance of 70 patients and 80 time intervals, with a no-show probability of 0.1 and a weighting factor of 0.5. The time limit for the local search will remain 60 seconds. The search space of this instance is too large to be explored exhaustively \(3.55 * 10^{43}\) possible solutions).

# --- Parameters for fig7 ---
N_patients_fig7 = 70
T_intervals_fig7 = 80
q_fig7 = 0.1
w_fig7 = 0.5
convolutions_projections_fig7 = compute_convolutions(s_dist, N_patients_fig7, q_fig7)
ls_time_limit = 60
target_sum = max(0, N_patients_fig7 - T_intervals_fig7 + 2)  # Initialize start + taol sum. Start from max(0, N - (T - 2)) to ensure non-negative start
stop_sum = 30
lowest_objective_value = float('inf')  # Initialize to a very high value

while target_sum <= stop_sum:
  found_lower = False
  num_center_patients = N_patients_fig7 - target_sum # Patients in the center part
  num_center_intervals = T_intervals_fig7 - 2 # Intervals available for the center part

  if num_center_patients < 0: # Should not happen if generate_start_tail_distribution is correct
      # print(f"Warning: Negative center patients ({num_center_patients}) for start={start_count}, tail={tail_count}. Skipping.")
      continue
  
  if num_center_intervals < 0:
      if num_center_patients > 0:
          # print(f"Warning: Cannot schedule {num_center_patients} center patients in {num_center_intervals} intervals. Skipping.")
          continue
      # If num_center_patients is 0 and num_center_intervals is < 0, center_schedule should be empty.
      center_schedule_segment = []
  elif num_center_patients == 0: # No patients for the center part
       center_schedule_segment = [0] * num_center_intervals
  else:
      center_schedule_segment = generate_evenly_distributed_schedule_intervals(
          num_center_patients,
          num_center_intervals
      )
  # Try all possible start/tail combinations that sum to target_sum
  for s in range(target_sum + 1):
    schedule_parts = [s]  # Start with the start count
    t = target_sum - s
    schedule_parts.extend(center_schedule_segment)
    schedule_parts.append(t)
    ewt, esp = calculate_objective_serv_time_lookup(
        schedule_parts, d_interval_len, convolutions_projections_fig7
    )
    new_objective_value = w_fig7 * ewt + (1 - w_fig7) * esp
    if new_objective_value < lowest_objective_value:
        lowest_objective_value = new_objective_value
        best_schedule = schedule_parts
        print(f"New lowest objective value found: {lowest_objective_value} for target_sum={target_sum}, schedule={schedule_parts}. Sum schedule={sum(schedule_parts)}")
        found_lower = True
  if not found_lower:
    print(f"No lower objective value found for target_sum={target_sum}. Stopping search.\nBest schedule so far: {best_schedule} with objective value {lowest_objective_value}")
    break 
  target_sum += 1

v_star = get_v_star(T_intervals_fig7)
x_star, c_star, log = local_search_w_timer(best_schedule, d_interval_len, convolutions_projections_fig7, w_fig7, v_star, T_intervals_fig7, time_limit = ls_time_limit, process_start_time = time.time(), echo = False)
print(f"Final best schedule found: {x_star} with objective value {c_star}")

New lowest objective value found: 1771.827236447373 for target_sum=0, schedule=[0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0]. Sum schedule=70
New lowest objective value found: 1620.4063636957135 for target_sum=1, schedule=[0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1]. Sum schedule=70
New lowest objective value found: 1507.063469941892 for target_sum=2, schedule=[0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 2]. Sum schedule=70
New lowest objective value found: 1500.9548912628877 for target_sum=2, schedule=[1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1]. Sum schedule=70
New lowest objective value found: 1380.2943800415296 for target_sum=3, schedule=[0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 3]. Sum schedule=70
New lowest objective value found: 1369.3656794891144 for target_sum=3, schedule=[1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 2]. Sum schedule=70
New lowest objective value found: 1264.0796111656696 for target_sum=4, schedule=[0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 4]. Sum schedule=70
New lowest objective value found: 1248.1616064148534 for target_sum=4, schedule=[1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 3]. Sum schedule=70
New lowest objective value found: 1173.4808033199697 for target_sum=5, schedule=[0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 5]. Sum schedule=70
New lowest objective value found: 1145.5298976670124 for target_sum=5, schedule=[1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 4]. Sum schedule=70
New lowest objective value found: 1094.309137300939 for target_sum=6, schedule=[0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 6]. Sum schedule=70
New lowest objective value found: 1062.2782285622934 for target_sum=6, schedule=[1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 5]. Sum schedule=70
New lowest objective value found: 1022.1716506665317 for target_sum=7, schedule=[0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 7]. Sum schedule=70
New lowest objective value found: 983.8478580230372 for target_sum=7, schedule=[1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 6]. Sum schedule=70
New lowest objective value found: 962.1106585777768 for target_sum=8, schedule=[0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 8]. Sum schedule=70
New lowest objective value found: 918.5589667881738 for target_sum=8, schedule=[1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 7]. Sum schedule=70
New lowest objective value found: 876.4091655813788 for target_sum=9, schedule=[1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 8]. Sum schedule=70
New lowest objective value found: 847.2180859119298 for target_sum=10, schedule=[1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 9]. Sum schedule=70
New lowest objective value found: 831.997487781327 for target_sum=11, schedule=[1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 10]. Sum schedule=70
New lowest objective value found: 824.2552651460912 for target_sum=12, schedule=[2, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 10]. Sum schedule=70
No lower objective value found for target_sum=13. Stopping search.
Best schedule so far: [2, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 10] with objective value 824.2552651460912
Final best schedule found: [2 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1
 1 0 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 1
 0 1 1 1 1 9] with objective value 820.9355371642945

# --- Plotting ---

# Run the optimization and get the log
# --- Run 1: Heuristic + Local Search ---
heuristic_log, total_time = find_optimal_schedule(
    N_patients=N_patients_fig7,
    T_intervals=T_intervals_fig7,
    d_interval_len=d_interval_len,
    w=w_fig7,
    convolutions_projections=convolutions_projections_fig7,
    heuristic_stop_limit=N_patients_fig7,
    local_search_time_limit=ls_time_limit
)

# --- Run 2: Local Search Only (for the same total duration) ---

ls_only_log, _ = find_optimal_schedule(
    N_patients=N_patients_fig7,
    T_intervals=T_intervals_fig7,
    d_interval_len=d_interval_len,
    w=w_fig7,
    convolutions_projections=convolutions_projections_fig7,
    heuristic_stop_limit=N_patients_fig7,
    local_search_time_limit=total_time, # Use same total time as the first run
    start_with_heuristic=False
)

# --- Combine logs and Plot ---
full_log = heuristic_log + ls_only_log

df = pd.DataFrame(full_log)
fig7 = go.Figure()
colors = {
    'heuristic': '#636EFA', 
    'local_search': '#00CC96', 
    'local_search_from_start': '#EF553B'
}

for source_name, group in df.groupby('source'):
    fig7.add_trace(go.Scatter(
        x=group['time_elapsed'], y=group['cost'],
        name=source_name, mode='lines+markers',
        line=dict(shape='hv', color=colors.get(source_name)),
        marker=dict(size=8, line=dict(width=1, color='DarkSlateGrey'))
    ))

fig7.update_layout(
    title=f"Optimization Strategy Comparison<sub>N = {N_patients_fig7}, T = {T_intervals_fig7}, q = {q_fig7} , w = {w_fig7}</sub>",
    xaxis_title="Time Elapsed (seconds)",
    yaxis_title="Objective Function Cost",
    legend_title_text="Search Method",
    template="plotly_white"
)

print("\nDisplaying convergence plot...")
fig7.show()

🚀 Starting optimization for 70 patients with Heuristic Search.

--- Stage 1: Heuristic Search ---

✅ Stage 1 Complete. Best Heuristic Schedule: [2, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 10]

--- Stage 2: Local Search Refinement ---

🏁 Optimization complete in 60.18 seconds.
Best schedule found: [2, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 9] with cost 820.9355
🚀 Starting optimization for 70 patients with Local Search only.
Using heuristic's first candidate as starting schedule: [0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0]

--- Local Search Only ---

🏁 Optimization complete in 60.77 seconds.
Best schedule found: [1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 7] with cost 859.4951

Displaying convergence plot...

Altman, Eitan, Bruno Gaujal, and Arie Hordijk. 2000. “Multimodularity, Convexity, and Optimization Properties.” Mathematics of Operations Research 25 (2): 324–47. https://www.jstor.org/stable/3690584.

Kaandorp, Guido C., and Ger Koole. 2007. “Optimal Outpatient Appointment Scheduling.” Health Care Management Science 10 (3): 217–29. https://doi.org/10.1007/s10729-007-9015-x.