```python
import ambulance_game as abg

```


```python
import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
import itertools
import scipy as sci
import scipy.integrate


def build_states(threshold, system_capacity, buffer_capacity):
    if threshold > system_capacity:
        states_1 = [(0, v) for v in range(0, system_capacity + 1)]
        states_2 = [(1, system_capacity)]
        return states_1 + states_2

    states_1 = [(0, v) for v in range(0, threshold)]
    states_2 = [
        (u, v)
        for v in range(threshold, system_capacity + 1)
        for u in range(buffer_capacity + 1)
    ]
    all_states = states_1 + states_2

    return all_states


def visualise_ambulance_markov_chain(
    num_of_servers, threshold, system_capacity, buffer_capacity
):
    all_states = build_states(threshold, system_capacity, buffer_capacity)
    G = nx.MultiDiGraph()
    for _, origin_state in enumerate(all_states):
        for _, destination_state in enumerate(all_states):
            column_adjacent = (
                destination_state[0] - origin_state[0] == 1
                and destination_state[1] - origin_state[1] == 0
            )
            row_adjacent = (
                destination_state[1] - origin_state[1] == 1
                and destination_state[0] - origin_state[0] == 0
            )
            if row_adjacent or column_adjacent:
                G.add_edge(origin_state, destination_state, color="blue")

    plt.figure(figsize=(1.5 * (buffer_capacity + 1), 1.5 * system_capacity))
    pos = {state: list(state) for state in all_states}
    nx.draw_networkx_nodes(
        G,
        pos,
        node_size=2000,
        nodelist=[state for state in all_states if state[1] < num_of_servers],
    )
    nx.draw_networkx_nodes(
        G,
        pos,
        node_size=2000,
        nodelist=[state for state in all_states if state[1] >= num_of_servers],
        node_color="red",
    )
    nx.draw_networkx_edges(G, pos)
    nx.draw_networkx_labels(G, pos)

    plt.axis("off")

    return G


def get_transition_matrix_entry(
    origin, destination, threshold, lambda_2, lambda_1, Lambda, mu, num_of_servers
):
    row_diff = origin[0] - destination[0]
    column_diff = origin[1] - destination[1]

    if row_diff == 0 and column_diff == -1:
        if origin[1] < threshold:
            return Lambda
        return lambda_1
    elif row_diff == -1 and column_diff == 0:
        return lambda_2
    elif (row_diff == 0 and column_diff == 1) or (
        row_diff == 1 and column_diff == 0 and origin[1] == threshold
    ):
        if origin[1] <= num_of_servers:
            return origin[1] * mu
        else:
            return num_of_servers * mu
    else:
        return 0


def get_symbolic_transition_matrix(
    num_of_servers, threshold, system_capacity, buffer_capacity
):
    Lambda = sym.symbols("Lambda")
    lambda_1 = sym.symbols("lambda") ** sym.symbols("o")
    lambda_2 = sym.symbols("lambda") ** sym.symbols("A")
    mu = sym.symbols("mu")

    all_states = build_states(threshold, system_capacity, buffer_capacity)
    Q = sym.zeros(len(all_states))

    if threshold > system_capacity:
        threshold = system_capacity

    for (i, origin_state), (j, destination_state) in itertools.product(
        enumerate(all_states), repeat=2
    ):
        Q[i, j] = get_transition_matrix_entry(
            origin=origin_state,
            destination=destination_state,
            threshold=threshold,
            lambda_2=lambda_2,
            lambda_1=lambda_1,
            Lambda=Lambda,
            mu=mu,
            num_of_servers=num_of_servers,
        )

    sum_of_rates = -np.sum(Q, axis=1)
    Q = Q + sym.Matrix(np.diag(sum_of_rates))

    return Q


def get_transition_matrix(
    lambda_2, lambda_1, mu, num_of_servers, threshold, system_capacity, buffer_capacity
):
    all_states = build_states(threshold, system_capacity, buffer_capacity)
    size = len(all_states)
    Q = np.zeros((size, size))

    if threshold > system_capacity:
        threshold = system_capacity

    for (i, origin_state), (j, destination_state) in itertools.product(
        enumerate(all_states), repeat=2
    ):
        Q[i, j] = get_transition_matrix_entry(
            origin=origin_state,
            destination=destination_state,
            threshold=threshold,
            lambda_2=lambda_2,
            lambda_1=lambda_1,
            Lambda=lambda_2 + lambda_1,
            mu=mu,
            num_of_servers=num_of_servers,
        )
    sum_of_rates = np.sum(Q, axis=1)
    np.fill_diagonal(Q, -sum_of_rates)
    return Q


def convert_symbolic_transition_matrix(Q_sym, lambda_2, lambda_1, mu):
    sym_Lambda = sym.symbols("Lambda")
    sym_lambda_1 = sym.symbols("lambda") ** sym.symbols("o")
    sym_lambda_2 = sym.symbols("lambda") ** sym.symbols("A")
    sym_mu = sym.symbols("mu")

    Q = np.array(
        Q_sym.subs(
            {
                sym_Lambda: lambda_2 + lambda_1,
                sym_lambda_1: lambda_1,
                sym_lambda_2: lambda_2,
                sym_mu: mu,
            }
        )
    ).astype(np.float64)
    return Q


def is_steady_state(state, Q):
    return np.allclose(np.dot(state, Q), 0)


def get_steady_state_numerically(
    Q, max_t=100, number_of_timepoints=1000, integration_function=sci.integrate.odeint
):
    def derivative_odeint(x, t):
        return np.dot(x, Q)

    def derivative_solve_ivp(t, x):
        return np.dot(x, Q)

    dimension = Q.shape[0]
    state = np.ones(dimension) / dimension
    while not is_steady_state(state=state, Q=Q):
        t_span = np.linspace(0, max_t, number_of_timepoints)
        if integration_function == sci.integrate.odeint:
            sol = integration_function(func=derivative_odeint, y0=state, t=t_span)
            state = sol[-1]
        elif integration_function == sci.integrate.solve_ivp:
            sol = integration_function(
                fun=derivative_solve_ivp, y0=state, t_span=t_span
            )
            state = sol.y[:, -1]
    return state


def augment_Q(Q):
    dimension = Q.shape[0]
    M = np.vstack((Q.transpose()[:-1], np.ones(dimension)))
    b = np.vstack((np.zeros((dimension - 1, 1)), [1]))
    return M, b


def get_steady_state_algebraically(Q, algebraic_function=np.linalg.solve):
    M, b = augment_Q(Q)
    if algebraic_function == np.linalg.solve:
        state = algebraic_function(M, b).transpose()[0]
    elif algebraic_function == np.linalg.lstsq:
        state = algebraic_function(M, b, rcond=None)[0][:, 0]
    return state


def get_markov_state_probabilities(
    pi, all_states, output=np.ndarray, system_capacity=None, buffer_capacity=None
):
    if output == dict:
        states_probabilities_dictionary = {}
        for i in range(len(all_states)):
            states_probabilities_dictionary[all_states[i]] = pi[i]
        return states_probabilities_dictionary
    elif output == np.ndarray:
        if buffer_capacity == None:
            buffer_capacity = max([state[0] for state in all_states])
        if system_capacity == None:
            system_capacity = max([state[1] for state in all_states])
        states_probabilities_array = np.full(
            (buffer_capacity + 1, system_capacity + 1), np.NaN
        )
        for index in range(len(all_states)):
            states_probabilities_array[all_states[index]] = pi[index]
        return states_probabilities_array


def get_mean_number_of_patients_in_system(pi, states):
    states = np.array(states)
    mean_patients_in_system = np.sum((states[:, 0] + states[:, 1]) * pi)
    return mean_patients_in_system


def get_mean_number_of_patients_in_hospital(pi, states):
    states = np.array(states)
    mean_patients_in_hospital = np.sum(states[:, 1] * pi)
    return mean_patients_in_hospital


def get_mean_ambulances_blocked(pi, states):
    states = np.array(states)
    mean_ambulances_blocked = np.sum(states[:, 0] * pi)
    return mean_ambulances_blocked

```


```python
lambda_2 = 0.05
lambda_1 = 0.1
mu = 0.1

num_of_servers = 2
threshold = 4
system_capacity = 3
buffer_capacity = 2

```


```python
all_states = abg.markov.build_states(threshold, system_capacity, buffer_capacity)
all_states

```




    [(0, 0), (0, 1), (0, 2), (0, 3)]




```python
Q_1 = abg.markov.get_transition_matrix(
    lambda_2, lambda_1, mu, num_of_servers, threshold, system_capacity, buffer_capacity
)
Q_1

```




    array([[-0.15,  0.15,  0.  ,  0.  ],
           [ 0.1 , -0.25,  0.15,  0.  ],
           [ 0.  ,  0.2 , -0.35,  0.15],
           [ 0.  ,  0.  ,  0.2 , -0.2 ]])




```python
Q_sym = abg.markov.get_symbolic_transition_matrix(
    num_of_servers, threshold, system_capacity, buffer_capacity
)
Q_sym

```




$\displaystyle \left[\begin{matrix}- \Lambda & \Lambda & 0 & 0\\\mu & - \Lambda - \mu & \Lambda & 0\\0 & 2 \mu & - \Lambda - 2 \mu & \Lambda\\0 & 0 & 2 \mu & - 2 \mu\end{matrix}\right]$




```python
Q_2 = abg.markov.convert_symbolic_transition_matrix(Q_sym, lambda_2, lambda_1, mu)
Q_2

```




    array([[-0.15,  0.15,  0.  ,  0.  ],
           [ 0.1 , -0.25,  0.15,  0.  ],
           [ 0.  ,  0.2 , -0.35,  0.15],
           [ 0.  ,  0.  ,  0.2 , -0.2 ]])




```python
Q_1 == Q_2

```




    array([[ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True],
           [ True,  True,  True,  True]])




```python
pi = abg.markov.get_steady_state_algebraically(Q_1)
pi

```




    array([0.22377622, 0.33566434, 0.25174825, 0.18881119])




```python
abg.markov.get_markov_state_probabilities(pi, all_states)

```




    array([[0.22377622, 0.33566434, 0.25174825, 0.18881119]])




```python
abg.markov.get_markov_state_probabilities(
    pi, all_states, buffer_capacity=buffer_capacity
)

```




    array([[0.22377622, 0.33566434, 0.25174825, 0.18881119],
           [       nan,        nan,        nan,        nan],
           [       nan,        nan,        nan,        nan]])




```python

```


```python

```


```python

```


```python

```


```python

```

40fc813e094bd395711548224f54074ad46eef61
ipynb
Jupyter Notebook
nbs/src/Markov/other/fix_markov_formulation.ipynb
11michalis11/AmbulanceDecisionGame
45164ba51da0417297f715e41716cb91facc120f
MIT
nbs/src/Markov/other/fix_markov_formulation.ipynb
11michalis11/AmbulanceDecisionGame
45164ba51da0417297f715e41716cb91facc120f
MIT
2020-04-20T09:08:31.000Z
2021-09-23T11:09:25.000Z
nbs/src/Markov/other/fix_markov_formulation.ipynb
11michalis11/AmbulanceDecisionGame
45164ba51da0417297f715e41716cb91facc120f
MIT
Qwen/Qwen-72B
1. YES
2. YES