In [ ]:

    

from __future__ import division

from scipy import optimize
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 [ ]:

    

renormalize = lambda x: x / x.sum()


    

In [ ]:

    

num_decks = 5
review_rate_budget = 0.1902
arrival_rate = 0.015
difficulty = 0.0077
capacity = 100
service_rates = np.append(renormalize(1 / np.sqrt(np.arange(1, num_decks + 1, 1))) * review_rate_budget, arrival_rate)


    

In [ ]:

    

def mva_throughput(init_throughput, service_rates=service_rates, capacity=capacity):
    P = np.zeros((num_decks + 1, num_decks + 1))
    for i in xrange(num_decks):
        P[i, i + 1] = (service_rates[i] - init_throughput) / (service_rates[i] - init_throughput + difficulty / (i + 1))
    P[num_decks, 0] = 1
    P[0, 0] = 1 - (service_rates[0] - init_throughput) / (service_rates[0] - init_throughput + difficulty)
    for i in xrange(1, num_decks):
        P[i, i - 1] = 1 - (service_rates[i] - init_throughput) / (service_rates[i] - init_throughput + difficulty / (i + 1))
                
    v = renormalize(np.linalg.eig(P.T)[1][:, 0]).astype(float)
        
    L = np.zeros(num_decks + 1)
    for m in xrange(1, capacity + 1):
        W = (L + 1) / service_rates
        actual_throughput = m / W.dot(v)
        L = actual_throughput * W * v

    return actual_throughput


    

In [ ]:

    

max_throughput = service_rates.min()
init_throughputs = np.arange(0, max_throughput, max_throughput / 100)


    

In [ ]:

    

actual_throughputs = [mva_throughput(x) for x in init_throughputs]


    

In [ ]:

    

plt.xlabel(r'$\lambda_{in}$')
plt.ylabel(r'$\lambda_{out}$')
mx = max(max(init_throughputs), max(actual_throughputs))
mx = 1.1 * mx
plt.xlim([0, mx])
plt.ylim([0, mx])
plt.plot([0, mx], [0, mx], linestyle='--')
plt.scatter(init_throughputs, actual_throughputs, linewidth=0)
plt.show()


    

In [ ]:

    

def mva_fixpoint(capacity, service_rates=service_rates):
    try:
        return optimize.fixed_point(lambda x: mva_throughput(x, service_rates=service_rates, capacity=capacity), 0)
    except:
        return np.nan


    

In [ ]:

    

capacities = np.arange(1, 1000, 10)


    

In [ ]:

    

throughputs = [mva_fixpoint(x) for x in capacities]


    

In [ ]:

    

plt.xlabel('Capacity')
plt.ylabel('Throughput')
plt.scatter(capacities, throughputs)
plt.show()


    

In [ ]: