This document demonstrate how to use the library to define a "density dependent population process" and to compute its mean-field approximation and refined mean-field approximation
In [1]:
# To load the library
import src.rmf_tool as rmf
import importlib
importlib.reload(rmf)
# To plot the results
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
This example illustrate two ways of constructing models:
The transitions for this model are:
with the convention : $x_{-1} = 1$ and $x_{K} = 0$ for some $K$.
In [2]:
# This code creates an object that represents a "density dependent population process"
ddpp = rmf.DDPP()
K = 9 # We use K=9 otherwise the tool is too slow for the 1/N^2 term of transient analysis.
# If one is only interested in the 1/N term, one can set K=20 or more.
rho=0.9
# The vector 'e(i)' is a vector where the $i$th coordinate is equal to $1$ (the other being equal to $0$)
def e(i):
l = np.zeros(K)
l[i] = 1
return(l)
# We then add the transitions :
for i in range(K):
if i>=1:
ddpp.add_transition(e(i),eval('lambda x: rho*(x[{}]*x[{}] - x[{}]*x[{}] )'.format(i-1,i-1,i,i) ))
if i<K-1:
ddpp.add_transition(-e(i),eval('lambda x: (x[{}] - x[{}])'.format(i,i+1) ))
ddpp.add_transition(e(0), lambda x : rho*(1-x[0]*x[0]))
ddpp.add_transition(-e(K-1), lambda x : x[K-1])
In [ ]:
In [ ]:
In [18]:
rho = 0.95
ddpp.set_initial_state(e(0)) # We first need to define an initial stater
%time pi,V,W = ddpp.meanFieldExpansionSteadyState(order=1)
T,X = ddpp.simulate(100,time=30) # We first plot a trajectory for $N=100$
plt.plot(T,X)
%time T,X = ddpp.simulate(1000,time=30) # Then for $N=1000$
plt.plot(T,X,'--')
Out[18]:
In [4]:
plt.figure()
ddpp.plot_ODE_vs_simulation(N=50)
plt.figure()
ddpp.plot_ODE_vs_simulation(N=1000)
(reference to be added)
This class also contains some functions to compute the fixed point of the mean-field approximation, to compute the "refined mean-field approximation" and to compare it with simulations.
If $\pi$ is the fixed point of the ODE, and $V$ the constant calculated by the function "meanFieldExpansionSteadyState(order=1)", then we have
$$E[X^N] = \pi + \frac1N V + o(\frac1N) $$To compute these constants :
In [17]:
rho=0.9
pi,V,W = ddpp.meanFieldExpansionSteadyState(order=1)
print(pi,V)
In [6]:
import misc.jsqD_simulate.average_valueJSQ as jsq
myN = [10,20,30,50,100]
Xs = {}
Vs = {}
Xrmf = {}
for N in myN:
# For steady-state simulation, 2 options :
# Option 1 : we could use the tool (to do so uncomment the line below)
#%time Xs[N],Vs[N] = ddpp.steady_state_simulation(N=N,time=100000/N)
# Option 2 : we use simulations from an external tool
Xs[N],Vs[N] = jsq.loadSteadyStateDistributionQueueLength(N=N,d=2,rho=rho,returnConfInterval=True)
Xrmf[N] = pi+V/N
print(N,'done')
In [7]:
# Average queue length :
for N in myN :
print(sum(Xs[N]),'+/-',sum(Vs[N]),'(Simulation, N={})'.format(N))
print(sum(Xrmf[N]),'(refined mean-field, N={})'.format(N))
In [8]:
for N in myN:
plt.figure()
plt.plot((Xs[N]),'+-')
plt.plot((Xrmf[N]),'*-')
plt.plot((pi),'--')
plt.legend(('Simulation','Refined mean-field','Mean-field'))
plt.xlabel('Queue length $i$')
plt.ylabel('Probability that a server has i or more jobs (in steady-state)')
plt.title('$N={}$'.format(N))
plt.xlim([0,10])
In [9]:
myN = [10,20,30,50,100]
print(' ',end='\t')
for N in myN : print(' N={} (simu)'.format(N),end='\t')
print()
for rho in [0.7,0.9,0.95]:
print('rho={} '.format(rho),end='\t')
pi,V,A, VWABCD = ddpp.meanFieldExpansionSteadyState(order=2)
for N in myN:
print('{0:.4f} ({1:.4f})'.format(sum(pi+V/N+A/N**2),
jsq.loadSteadyStateAverageQueueLength(rho=rho,N=N,d=2)),end='\t')
print('')
In [10]:
rho=0.9
pi,V,A, VWABCD = ddpp.meanFieldExpansionSteadyState(order=2)
for N in [10,20,30]:
Xsimu = jsq.loadSteadyStateDistributionQueueLength(N=N,d=2,rho=rho)
plt.figure()
for i in [1,2]:
plt.subplot(1,2,i)
if i==1: plt.plot((Xsimu),'+-'); plt.ylabel('Probability that a server has i or more jobs (in steady-state)')
else: plt.semilogy(Xsimu,'+-')
plt.plot((pi),'*-')
plt.plot((pi+V/N),'--')
plt.plot((pi+V/N+A/N**2),'--')
plt.legend(('Simulation','Mean field','1/N-expansion','1/N^2-expansion'))
plt.xlim([0,8])
plt.ylim([1e-5,1])
plt.xlabel('Queue length $i$') ;
plt.title('$N={}$'.format(N))
In [11]:
x0 = np.zeros(K)
x0[0]=1
x0[1]=1
x0[2]=0.8
rho=0.9
ddpp.set_initial_state(x0)
%time T,X,V,XVW=ddpp.meanFieldExpansionTransient(order=1,time=50)
In [12]:
N=10
Tsimu,Xsimu = jsq.loadTransientSimu(rho=0.9,d=2,N=10,initial_number_of_jobs=2.8)
plt.plot(Tsimu,np.sum(Xsimu,1)-1)
plt.plot(T,np.sum(XVW[:,0:K],1),'--')
plt.plot(T,np.sum(XVW[:,0:K]+XVW[:,K:2*K]/N,1),'-.')
plt.legend(('Simu','Mean field','Refined mean field ($1/N$-expansion)'))
Out[12]:
In [13]:
%time T,X,V,A,XVWABCD=ddpp.meanFieldExpansionTransient(order=2,time=50)
# Warning : can be slow
In [14]:
Tsimu,Xsimu = jsq.loadTransientSimu(rho=0.9,d=2,N=10,initial_number_of_jobs=2.8)
plt.plot(Tsimu,np.sum(Xsimu,1)-1)
X = XVWABCD[:,0:K]
V = XVWABCD[:,K:2*K]
A = XVWABCD[:,2*K+K**2:3*K+K**2]
plt.plot(T,np.sum(X,1),'--')
plt.plot(T,np.sum(X+V/N,1),'-.')
plt.plot(T,np.sum(X+V/N+A/N**2,1),'-.')
plt.legend(('Simu','Mean field','Refined mean field ($1/N$-expansion)','Refined mean field ($1/N^2$-expansion)'))
Out[14]:
In [ ]: