Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Views: 60

Measuring the Price of Anarchy in CCU Interactions.

(Knight, Komenda, Griffiths)

This worksheet contains code used to obtain the graphics for the paper entitled: 'Measuring the Price of Anarchy in CCU Interaction'.

The abstract of the paper is:

"Hospital throughput is often studied and optimised in isolation, ignoring the interactions between hospitals. In this paper Critical Care Unit (CCU) interaction is placed within a game theoretic framework. The methodology involves the use of a normal form game underpinned by a two dimensional continuous Markov chain. A theoretical result is given that ensures the existence of a Nash equilibrium.

The effect of target policies is investigated justifying the use of these to align the interests of individual hospitals and social welfare. In particular, we identify the lowest value of a utilisation target that aligns these."

Best response for t=.8

NHbr = [[0, 7],
 [1, 7],
 [2, 7],
 [3, 7],
 [4, 7],
 [5, 8],
 [6, 8],
 [7, 8],
 [8, 8],
 [9, 8],
 [10, 8],
 [11, 8],
 [12, 8],
 [13, 8],
 [14, 8],
 [15, 8],
 [16, 8]]
RGbr = [[14, 0],
 [15, 1],
 [15, 2],
 [16, 3],
 [16, 4],
 [16, 5],
 [16, 6],
 [16, 7],
 [16, 8]]
p = list_plot(NHbr, plotjoined=True, linestyle='', marker='D', legend_label='Best response to RG', color='grey', ticks=[range(1,17),range(1,9)])
p += list_plot(RGbr, plotjoined=True, linestyle='', marker='x', color='black', legend_label='Best response to NH')
p += circle([16,8],.25)
p.axes_labels(['$K_{RG}$', '$K_{NH}$'])
p.set_legend_options(loc=[.1,.1])
p.axes_labels_size(3)
p.save('best_responses_model1.pdf')
p

Best response graph for t=.6

NHbr = [[0, 5],
 [1, 5],
 [2, 5],
 [3, 5],
 [4, 5],
 [5, 6],
 [6, 6],
 [7, 6],
 [8, 6],
 [9, 6],
 [10, 6],
 [11, 6],
 [12, 6],
 [13, 6],
 [14, 6],
 [15, 7],
 [16, 7]]
RGbr = [[10, 0],
 [11, 1],
 [11, 2],
 [11, 3],
 [11, 4],
 [11, 5],
 [11, 6],
 [12, 7],
 [12, 8]]
p = list_plot(NHbr, plotjoined=True, linestyle='', marker='D', legend_label='Best response to RG', color='grey', ticks=[range(1,17),range(1,9)])
p += list_plot(RGbr, plotjoined=True, linestyle='', marker='x', color='black', legend_label='Best response to NH')
p += circle([11,6],.25)
p.axes_labels(['$K_{RG}$', '$K_{NH}$'])
p.set_legend_options(loc=[.1,.1])
p.axes_labels_size(3)
p.save('best_responses_model1_t=0_6.pdf')
p

Function to plot heatmap

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
from numpy.random import randn
from pylab import cm, get_cmap

def cmap_map(function,cmap):
    """ Applies function (which should operate on vectors of shape 3:
    [r, g, b], on colormap cmap. This routine will break any discontinuous     points in a colormap.
    """
    cdict = cmap._segmentdata
    step_dict = {}
    # Firt get the list of points where the segments start or end
    for key in ('red','green','blue'):         step_dict[key] = map(lambda x: x[0], cdict[key])
    step_list = sum(step_dict.values(), [])
    step_list = array(list(set(step_list)))
    # Then compute the LUT, and apply the function to the LUT
    reduced_cmap = lambda step : array(cmap(step)[0:3])
    old_LUT = array(map( reduced_cmap, step_list))
    new_LUT = array(map( function, old_LUT))
    # Now try to make a minimal segment definition of the new LUT
    cdict = {}
    for i,key in enumerate(('red','green','blue')):
        this_cdict = {}
        for j,step in enumerate(step_list):
            if step in step_dict[key]:
                this_cdict[step] = new_LUT[j,i]
            elif new_LUT[j,i]!=old_LUT[j,i]:
                this_cdict[step] = new_LUT[j,i]
        colorvector=  map(lambda x: x + (x[1], ), this_cdict.items())
        colorvector.sort()
        cdict[key] = colorvector

    return matplotlib.colors.LinearSegmentedColormap('colormap',cdict,1024)

def makedict(lst):
    d = {}
    for row in lst:
        if row[0] in d:
            d[row[0]][row[1]] = row[2]
        else:
            d[row[0]] = {row[1] : row[2]}
    return d

def makematrix(lst):
    d = makedict(lst)
    xset = sorted(list(set(([k[0] for k in lst]))))
    yset = sorted(list(set(([k[1] for k in lst]))))
    m = len(xset)
    n = len(yset)
    M = matrix(RR, m, n)
    indexdict = {}
    for i in range(m):
        for j in range(n):
            if xset[i] in d:
                if yset[j] in d[xset[i]]:
                    M[i, j] = d[xset[i]][yset[j]]
    return M

def heatplot(lst, title='heatplot.pdf', xlabel='$x$', ylabel='$y$', titlelabel=False, redoticks=True, incrementticks=False, clr='gist_yarg',cbticks=False,cblabels=False):
    fig, ax = plt.subplots()

    xvalues = [row[1] for row in lst]
    yvalues = [row[0] for row in lst]
    xmin = min(xvalues)
    xmax = max(xvalues)
    ymin = min(yvalues)
    ymax = max(yvalues)

    data = makematrix(lst)

    cax = ax.imshow(data, interpolation='nearest', cmap=get_cmap(clr), aspect='auto')
    if titlelabel:
        ax.set_title(titlelabel, fontsize="xx-large")

    cbar = fig.colorbar(cax)
    if cbticks:
        cbar.set_ticks(cbticks)
    if cblabels:
        cbar.set_ticklabels(cblabels)
    if redoticks:
        nbrxlabels = len(plt.xticks()[0])
        xlabs = srange(xmin, xmax, (xmax - xmin) / nbrxlabels)

        nbrylabels =len(plt.yticks()[0])
        ylabs = srange(ymin, ymax, (ymax - ymin) / (nbrylabels))

        ax.set_xticklabels([round(x,2) for x in xlabs])
        ax.set_yticklabels([round(y,2) for y in ylabs])


        ax.set_xticks(srange(0, max(plt.xticks()[0]), max(plt.xticks()[0])/nbrxlabels))
        ax.set_yticks(srange(0, max(plt.yticks()[0]), max(plt.yticks()[0])/nbrylabels))

    if incrementticks:
        xlabs = plt.xticks()[0]

        ylabs = plt.yticks()[0]

        ax.set_xticklabels([int(x + incrementticks) for x in xlabs])
        ax.set_yticklabels([int(y + incrementticks) for y in ylabs])

    plt.ylabel(xlabel, fontsize="x-large")
    plt.xlabel(ylabel, fontsize="x-large")
    plt.savefig(title)

Model 1

PoA

import csv
DATA = 'foo.data/'
f = open(DATA+'M1_target_demand', 'r')
data = csv.reader(f)
data = [row for row in data]
f.close()
data = [[round(eval(k),2) for k in row] for row in data[1:]]

heatplot([[row[0], row[1], ln(row[2])] for row in data if row[0] >= .1 and row[2] != +infinity], title='lnPoAmodel1targetvdemand.pdf', xlabel='$t$', ylabel='$x$', titlelabel='PoA',cbticks=[0,.25,.5,.75,1,1.25,1.5,1.75,2], cblabels=["%.01f" % e^x for x in [0,.25,.5,.75,1,1.25,1.5,1.75,2]])

heatplot([row for row in data if row[0] < .4], title='model1targetvdemandfort<4.pdf', xlabel='$t$', ylabel='$x$', titlelabel='PoA')

heatplot([row for row in data if row[0] > .25], title='model1targetvdemandfort>25.pdf', xlabel='$t$', ylabel='$x$', titlelabel='PoA')

heatplot([row for row in data if row[0] > .5], title='model1targetvdemandfort>5.pdf', xlabel='$t$', ylabel='$x$', titlelabel='PoA')

Lowest target for PoA = 1

def PoAvxdict(lst):
    d = {}
    for row in lst:
        if row[1] in d:
            d[row[1]][row[0]] = row[2]
        else:
            d[row[1]] = {row[0] : row[2]}
    return d

def plotmint(lst):
    l = []
    d = PoAvxdict(lst)
    for x in d:
        minPoA = min([d[x][k] for k in d[x]])
        l.append([x, min([k for k in d[x] if d[x][k] == minPoA])])
    return list_plot(l, size=30, color='black')

p = plotmint(data)
p.axes_labels(['$x$', 'argmin$_t($PoA$(x))$'])
p.ymin(0)
p.ymax(1)
p.axes_labels_size(2)
p.save('argminPoAmodel1.pdf')
p

Model 2

import csv
f = open(DATA+'M2_target_demand', 'r')
data = csv.reader(f)
data = [row for row in data]
f.close()
data = [[round(eval(k),2) for k in row] for row in data[1:]]

heatplot([[row[0], row[1], row[2]] for row in data if row[0] >= .1 and row[2] != +infinity], title='model2targetvdemand.pdf', xlabel='$t$', ylabel='$x$', titlelabel='PoA')

p = plotmint(data)
p.axes_labels(['$x$', 'argmin$_t($PoA$(x))$'])
p.ymin(0)
p.ymax(1)
p.axes_labels_size(2)
p.save('argminPoAmodel2.pdf')
p

Capacity Plots

import csv
f = open(DATA+'M1_capacity', 'r')
data = csv.reader(f)
data = [row for row in data]
f.close()
FALSE = False
data = [[round(eval(k),2) for k in row] for row in data[1:]]
data = [row for row in data if row[0]>1 and row[1]>1]

heatplot([row for row in data if row[1]<22], title='model1capacity.pdf', xlabel='$c_{NH}$', ylabel='$c_{RG}$', titlelabel='PoA', redoticks=False, incrementticks = 2)

import csv
f = open(DATA+'M2_capacity', 'r')
data = csv.reader(f)
data = [row for row in data]
f.close()
data = [[round(eval(k),2) for k in row] for row in data[1:]]
data = [row for row in data if row[0]>1 and row[1]>1]

heatplot([row for row in data if row[1]<22], title='model2capacity.pdf', xlabel='$c_{NH}$', ylabel='$c_{RG}$', titlelabel='PoA', redoticks=False, incrementticks = 2)

Steady state probabilities

NH = [0.0024678535857431323,
 0.015726290629619416,
 0.05091770938993938,
 0.11226713892826898,
 0.19131907048179278,
 0.27338253952862523,
 0.3539193974560111,
 0.0,
 0.0]
RG =  [6.429788681200632e-06,
 7.867929699401477e-05,
 0.0004834366461007307,
 0.0019895883556981793,
 0.006173269294301033,
 0.015414004918794444,
 0.03229069468890144,
 0.0584485711739349,
 0.09349055513077142,
 0.13464969759037496,
 0.1777308502600037,
 0.21930666447309408,
 0.2599375583823499,
 0.0,
 0.0,
 0.0,
 0.0]

p = list_plot(NH, size=40, axes_labels=['$n$', '$P(n)$'], color='black',ticks=[range(1,9),[0,.2,.4,.6,.8,1]])
p += line([[6,0], [6,1]], color='red', linestyle='dashed')
p += text('$K_{NH}$',(6.2,1.02), color='black', fontsize=20)
p.ymax(1)
p.axes_labels_size(2)
p.save('NH.pdf')
p

p = list_plot(RG, size=40, axes_labels=['$n$', '$P(n)$'], color='black', ticks=[range(1,17),[0,.2,.4,.6,.8,1]])
p += line([[12,0], [12,1]], color='red', linestyle='dashed')
p += text('$K_{RG}$',(12.5,1.02), color='black', fontsize=20)
p.ymax(1)
p.axes_labels_size(2)
p.save('RG.pdf')
p

NH = [0.11446786319342364, 0.8855321368065764, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
RG = [2.846728005972963e-07,
 4.723444806693106e-06,
 3.943074788857625e-05,
 0.00022076875718588,
 0.0009325224253835039,
 0.003169530627486097,
 0.009029855901531656,
 0.02218313581230128,
 0.047989623116679016,
 0.09294363028933814,
 0.16338141261963893,
 0.26387851375924976,
 0.39622656782570986,
 0.0,
 0.0,
 0.0,
 0.0]

p = list_plot(NH, size=40, axes_labels=['$n$', '$P(n)$'], color='black',ticks=[range(1,9),[0,.2,.4,.6,.8,1]])
p += line([[1,0], [1,1]], color='red', linestyle='dashed')
p += text('$K_{NH}$',(1.2,1.02), color='black', fontsize=20)
p.ymax(1)
p.axes_labels_size(2)
p.save('NH_1_12.pdf')
p

p = list_plot(RG, size=40, axes_labels=['$n$', '$P(n)$'], color='black', ticks=[range(1,17),[0,.2,.4,.6,.8,1]])
p += line([[12,0], [12,1]], color='red', linestyle='dashed')
p += text('$K_{RG}$',(12.5,1.02), color='black', fontsize=20)
p.ymax(1)
p.axes_labels_size(2)
p.save('RG_1_12.pdf')
p

Example for Proof

NHbr = [[0, 2], [1, 4], [2, 4], [3, 4], [4, 4], [5, 5], [6, 5], [7, 5], [8, 5], [9, 5], [10, 6], [11, 6], [12, 6], [13, 7], [14, 7], [15, 8], [16, 8]]
RGbr = [[6, 0], [6, 1], [7, 2], [7, 3], [7, 4], [8, 5], [8, 6], [8, 7], [8, 8]]
p = list_plot(NHbr, plotjoined=True, linestyle='', marker='D', color='grey', ticks=[range(1,17),range(1,9)])
p.ymin(0)
p.ymax(8)
p.axes_labels(['$K_{RG}$', '$K_{NH}$'])
p.axes_labels_size(2)
p.save('best_responses_ex_for_proof_NH.pdf')
p

q = list_plot(RGbr, plotjoined=True, linestyle='',marker='x', color='black', ticks=[range(1,17),range(1,9)])
q.xmax(16)
q.xmin(0)
q.axes_labels(['$K_{RG}$', '$K_{NH}$'])
p.axes_labels_size(3)
q.save('best_responses_ex_for_proof_RG.pdf')
q

p = p+q
p += circle([8,5],.25)
p.axes_labels(['$K_{RG}$', '$K_{NH}$'])
p.set_legend_options(loc=[.1,.1])
p.axes_labels_size(2)
p.save('best_responses_ex_for_proof.pdf')
p

Picture for proof

m = [[1, .9], [2, .75], [3, .65], [4, .58], [5, .55], [6, .52], [7, .49], [8, .47], [9, .45], [10, .44]]
mminusone = [[row[0], row[1] -  .1*row[1]] for row in m]
mplusone = [[row[0], min(row[1] + .25*row[1],1)] for row in m]
mplustwo = [[row[0], min(row[1] + .15*row[1],1)] for row in mplusone]

p = list_plot(mplustwo, plotjoined=True, linestyle='',marker='^', color='black', legend_label='$K_{NH}=m+2$')
p += list_plot(mplusone, plotjoined=True, linestyle='',marker='x', color='black', legend_label='$K_{NH}=m+1$', ticks=[range(11),[.1,.2,.3,.4,.5,.6,.7,.8,.9,1]])
p += list_plot(m, plotjoined=True, linestyle='',marker='D', color='black', legend_label='$K_{NH}=m$')
p += list_plot(mminusone, plotjoined=True, linestyle='',marker='o', color='black', legend_label='$K_{NH}=m-1$')

p += line([[0,.7],[10,.7]], linestyle='dashed', color='red')
p += text('$t$', [10,.71], color='black')

locn = [4, 1]
p += line([[3,0],[3,.65]], linestyle='dashed', color='blue')
p += arrow(vector(locn) - vector([0,.025]),[3.05,.66], color='black')
p += text('$f_{NH}(3)=m$',locn, color='black')

locn = [4.65,.95]
p += line([[4,0],[4,.73]], linestyle='dashed', color='blue')
p += arrow(vector(locn) - vector([0,.025]),[4.05,.74], color='black')
p += text('$f_{NH}(4)=m+1$', locn, color='black')

locn = [5.55, .9]
p += line([[5,0],[5,.55*1.25]], linestyle='dashed', color='blue')
p += arrow(vector(locn) - vector([0,.025]),[5.05,.55*1.25+.01], color='black')
p += text('$f_{NH}(5)=m+1$', locn, color='black')

p.axes_labels(["$K_{RG}$", "$U_{NH}$"])
p.ymin(0)
p.axes_labels_size(2)
p.save('proof.pdf')
p