In [ ]:
from __future__ import division
from collections import defaultdict, Counter
from Queue import Queue, LifoQueue
import numpy as np
In [ ]:
from matplotlib import pyplot as plt
import seaborn as sns
%matplotlib inline
In [ ]:
import matplotlib as mpl
mpl.rc('savefig', dpi=300)
mpl.rc('text', usetex=True)
In the steady state, we have the following traffic equation. $$ \lambda_{ext} + (1 - \mathbb{E}[e^{-\theta D}])\lambda = \lambda < \mu $$
Assuming a FIFO service discipline, we can approximate delay as $D \sim Exponential(\mu - \lambda)$. Noticing that the expected recall rate is the moment-generating function for $D$,
$$ \mathbb{E}[e^{-\theta D}] = \frac{\mu - \lambda}{\mu - \lambda + \theta} $$Plugging this expression into the traffic equation,
$$ \lambda = \frac{\mu - \lambda + \theta}{\mu - \lambda} \cdot \lambda_{ext} < \mu $$
In [ ]:
T = 100000
In [ ]:
theta = 0.01
arrival_rate = 0.2
service_rate = 0.4
In [ ]:
using_feedback = True
In [ ]:
if using_feedback:
flow_rate = (service_rate + arrival_rate - np.sqrt((service_rate + arrival_rate)**2 - 4 * arrival_rate * (service_rate + theta))) / 2
else:
flow_rate = arrival_rate
print flow_rate
In [ ]:
if using_feedback:
expected_recall_rate = (service_rate - flow_rate) / (service_rate - flow_rate + theta)
else:
expected_recall_rate = 1
print expected_recall_rate
In [ ]:
assert service_rate > flow_rate / expected_recall_rate
In [ ]:
using_clocked_delays = True
using_fifo = True
using_balking = False
In [ ]:
feedback_threshold = np.log(service_rate / arrival_rate) / np.log(1 + theta / service_rate)
print feedback_threshold
In [ ]:
q = Queue() if using_fifo else LifoQueue()
exits = []
delays = []
In [ ]:
# initialize queue with some items
size = int(np.log(service_rate / arrival_rate) / np.log(1 + theta / service_rate)) + 10 # initial queue size
iits = np.random.exponential(1 / (arrival_rate + service_rate), size) # inter-arrival times
t = -sum(iits)
for iit in iits:
t += iit
q.put(t)
In [ ]:
# initialize empty queue
size = 0
In [ ]:
sizes = [size]
In [ ]:
t = 0
outcomes_at_queue_size = defaultdict(list)
for _ in xrange(T):
t += np.random.exponential(1 / (arrival_rate + service_rate))
if np.random.random() < arrival_rate / (arrival_rate + service_rate) and (not using_balking or size <= feedback_threshold):
q.put(t)
size += 1
elif not q.empty():
if using_clocked_delays:
delay = t - q.get()
else:
q.get()
delay = np.random.exponential(1 / (service_rate - flow_rate))
delays.append(delay)
if np.random.random() < np.exp(-theta*delay):
outcomes_at_queue_size[size].append(1)
exits.append(t)
size -= 1
elif using_feedback:
outcomes_at_queue_size[size].append(0)
q.put(t)
sizes.append(size)
In [ ]:
timesteps_at_queue_size = Counter(sizes)
In [ ]:
plt.xlabel(r'Queue Size $i$ (Number of Items)')
plt.ylabel(r'Fraction of Time Spent in State $i$')
queue_sizes, num_timesteps = zip(*timesteps_at_queue_size.items())
z = sum(num_timesteps)
plt.scatter(queue_sizes, [x / z for x in num_timesteps], label='Simulated')
queue_size_states = np.arange(0, max(queue_sizes)+1, 1)
pis = ((arrival_rate / service_rate)**(queue_size_states))*((1+theta/service_rate)**(queue_size_states*(queue_size_states+1)/2))
pis /= np.sum(pis)
plt.plot(queue_size_states, pis, label=r'$\pi_i$')
plt.legend(loc='best')
plt.show()
In [ ]:
plt.xlabel(r'Queue Size $i$ (Number of Items)')
plt.ylabel('Recall Likelihood at Head of Queue')
queue_sizes, outcomes = zip(*outcomes_at_queue_size.items())
plt.scatter(queue_sizes, [np.mean(x) for x in outcomes], label='Simulated')
queue_size_states = np.arange(0, max(queue_sizes)+1, 1)
plt.plot(queue_size_states, (1+theta/service_rate)**(-queue_size_states), label=r'$\left(1 + \frac{\theta}{\mu}\right)^{-i}$')
plt.legend(loc='lower left')
plt.show()
In [ ]:
plt.xlabel('Time')
plt.ylabel('Queue Size (Number of Items)')
plt.plot(sizes)
plt.show()
In [ ]:
num_reps = 5
#arrival_rates = np.arange(2*service_rate/20, 2*service_rate, 2*service_rate/20)
#thresholds = [10, 50, 100]
arrival_rates = [0.1, 0.3]#[0.25, 0.3, 0.35]
thresholds = np.arange(1, 251, 1)
In [ ]:
throughputs = [[[] for _ in arrival_rates] for _ in thresholds]
In [ ]:
for threshold_idx, threshold in enumerate(thresholds):
for arrival_rate_idx, arrival_rate in enumerate(arrival_rates):
print "%d %d %d %d" % (threshold_idx, len(thresholds), arrival_rate_idx, len(arrival_rates))
for _ in xrange(num_reps):
q = Queue()
num_exits = 0
t = 0
size = 0
for _ in xrange(T):
t += np.random.exponential(1 / (arrival_rate + service_rate))
if np.random.random() < arrival_rate / (arrival_rate + service_rate):
if size < threshold:
size += 1
q.put(t)
elif not q.empty():
if np.random.random() < np.exp(-theta*(t - q.get())):
size -= 1
num_exits += 1
else:
q.put(t)
throughputs[threshold_idx][arrival_rate_idx].append(num_exits / t)
In [ ]:
plt.xlabel(r'Arrival Rate $\lambda_{ext}$')
plt.ylabel(r'Throughput $\lambda_{out}$')
for threshold, throughputs_for_threshold in zip(thresholds, throughputs):
plt.errorbar(arrival_rates, [np.mean(x) for x in throughputs_for_threshold], yerr=[np.std(x)/np.sqrt(len(x)) for x in throughputs_for_threshold], label=r'$k = %d$' % threshold)
plt.legend(loc='best')
plt.show()
In [ ]:
num_reps = 5
thresholds = np.arange(1, 251, 5)
arrival_rate = 0.3
In [ ]:
throughputs = [[] for _ in thresholds]
In [ ]:
for threshold_idx, threshold in enumerate(thresholds):
print "%d %d" % (threshold_idx, len(thresholds))
for _ in xrange(num_reps):
q = Queue()
num_exits = 0
t = 0
size = 0
for _ in xrange(T):
t += np.random.exponential(1 / (arrival_rate + service_rate))
if np.random.random() < arrival_rate / (arrival_rate + service_rate):
if size < threshold:
size += 1
q.put(t)
elif not q.empty():
if np.random.random() < np.exp(-theta*(t - q.get())):
size -= 1
num_exits += 1
else:
q.put(t)
throughputs[threshold_idx].append(num_exits / t)
In [ ]:
pis_approx = np.zeros(len(thresholds))
pis_ubound = np.zeros(len(thresholds))
pis_lbound = np.zeros(len(thresholds))
for i, threshold in enumerate(thresholds):
queue_size_states = np.arange(0, threshold+1, 1)
pis = ((arrival_rate / service_rate)**(queue_size_states))*((1+theta/service_rate)**(queue_size_states*(queue_size_states+1)/2))
pis /= np.sum(pis)
pis_approx[i] = pis[-1]
pis = ((arrival_rate / service_rate)**(queue_size_states))*((1+theta/(arrival_rate+service_rate))**(queue_size_states*(queue_size_states+1)/2))
pis /= np.sum(pis)
pis_ubound[i] = pis[-1]
pis = ((arrival_rate / service_rate)**(queue_size_states))*((1+theta/(arrival_rate))**(queue_size_states*(queue_size_states+1)/2))
pis /= np.sum(pis)
pis_lbound[i] = pis[-1]
In [ ]:
plt.xlabel(r'Upper Bound on Queue Size $k$')
plt.ylabel(r'Throughput $\lambda_{out}$')
plt.plot(thresholds, arrival_rate*(1-pis_approx), label=r'$\lambda_{ext}(1-\pi^{\mu}_k)$')
plt.plot(thresholds, arrival_rate*(1-pis_ubound), label=r'$\lambda_{ext}(1-\pi^{\lambda+\mu}_k)$')
plt.plot(thresholds, arrival_rate*(1-pis_lbound), label=r'$\lambda_{ext}(1-\pi^{\lambda}_k)$')
plt.errorbar(thresholds, [np.mean(x) for x in throughputs], yerr=[np.std(x)/np.sqrt(len(x)) for x in throughputs], label='Simulated')
plt.legend(loc='best')
plt.show()
In [ ]:
num_reps = 50
init_sizes = range(int(2 * feedback_threshold))
In [ ]:
summ_sizes = [[] for _ in init_sizes]
In [ ]:
for init_size_idx, init_size in enumerate(init_sizes):
print "%d %d" % (init_size_idx, len(init_sizes))
for _ in xrange(num_reps):
q = Queue()
size = init_size
iits = np.random.exponential(1 / (arrival_rate + service_rate), size)
t = -sum(iits)
for iit in iits:
t += iit
q.put(t)
sizes = [size]
t = 0
for _ in xrange(T):
t += np.random.exponential(1 / (arrival_rate + service_rate))
if np.random.random() < arrival_rate / (arrival_rate + service_rate):
q.put(t)
size += 1
elif not q.empty():
if np.random.random() < np.exp(-theta*(t - q.get())):
size -= 1
else:
q.put(t)
sizes.append(size)
summ_sizes[init_size_idx].append((np.max(sizes), np.min(sizes)))
In [ ]:
plt.xlabel('Initial Queue Size (Number of Items)')
plt.ylabel('Extinction Probability')
plt.errorbar(init_sizes, [np.mean([1 if z[1] == 0 else 0 for z in x]) for x in summ_sizes], yerr=[np.std([1 if z[1] == 0 else 0 for z in x])/np.sqrt(len(x)) for x in summ_sizes])
plt.axvline(x=feedback_threshold, linestyle='--', label='Feedback Threshold')
plt.legend(loc='best')
plt.show()
In [ ]:
plt.xlabel('Simulated Delay')
plt.ylabel('Frequency (Number of Interactions)')
plt.hist(delays, bins=20, linewidth=0, alpha=0.5)
plt.show()
In [ ]:
w = np.zeros(int(np.ceil(exits[-1])))
for x in exits:
w[int(x)] += 1
In [ ]:
s = np.random.poisson(flow_rate * expected_recall_rate, 100000)
In [ ]:
plt.xlabel('Number of Exits in 1-Unit Window')
plt.ylabel('Frequency (Number of Windows)')
plt.hist(s, bins=20, alpha=0.5, linewidth=0, normed=True, label=r'Poisson($\lambda = %0.3f$)' % (flow_rate * expected_recall_rate))
plt.hist(w, bins=20, alpha=0.5, linewidth=0, normed=True, label='Simulated')
plt.legend(loc='best')
plt.show()
In [ ]:
ex = np.array(exits)
interarrival_times = ex[1:] - ex[:-1]
In [ ]:
s = np.random.exponential(1 / (flow_rate * expected_recall_rate), 100000)
In [ ]:
plt.xlabel('Inter-Arrival Time')
plt.ylabel('Frequency (Number of Exits)')
plt.hist(interarrival_times, bins=20, alpha=0.5, linewidth=0, normed=True, label=r'Exponential($\lambda = %0.3f$)' % (flow_rate * expected_recall_rate))
plt.hist(interarrival_times, bins=20, alpha=0.5, linewidth=0, normed=True, label='Simulated')
plt.legend(loc='best')
plt.show()
In [ ]:
scaling_factors = np.arange(1, 100, 2)
In [ ]:
simulated_exit_rates = []
for scaling_factor in scaling_factors:
q = Queue()
exits = []
t = 0
for _ in xrange(T):
t += np.random.exponential(1 / (scaling_factor * (arrival_rate + service_rate)))
if np.random.random() < arrival_rate / (arrival_rate + service_rate):
q.put(t)
elif not q.empty():
if np.random.random() < np.exp(-theta*scaling_factor*(t - q.get())):
exits.append(t)
elif using_feedback:
q.put(t)
simulated_exit_rates.append(len(exits) / (scaling_factor * t))
In [ ]:
plt.xlabel('Scaling Factor')
plt.ylabel('Exit Rate')
plt.plot(scaling_factors, simulated_exit_rates, label='Simulated')
plt.plot(scaling_factors, [flow_rate * expected_recall_rate] * len(scaling_factors), '--', label='Predicted')
plt.legend(loc='best')
plt.show()
In [ ]: