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)

With feedback

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 $$

Without feedback

$$ \lambda_{ext} = \lambda < \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 [ ]: