1 #!/usr/bin/env python
   2 
   3 labda, mu1, mu2 = 1., 1.01, 1.001
   4 N1, N2 = 50, 50
   5 size = N1*N2
   6 
   7 from numpy import ones, zeros, empty
   8 import scipy.sparse as sp
   9 import  pysparse 
  10 from pylab import matshow, savefig
  11 from scipy.linalg import norm
  12 import time
  13 
  14 def state(i,j):
  15     return j*N1 + i
  16 
  17 def fillOffDiagonal(Q):
  18     # labda
  19     for i in range(0,N1-1):
  20         for j in range(0,N2):
  21             Q[(state(i,j),state(i+1,j))]= labda
  22     # mu2
  23     for i in range(0,N1):
  24         for j in range(1,N2):
  25             Q[(state(i,j),state(i,j-1))]= mu2
  26     # mu1
  27     for i in range(1,N1):
  28         for j in range(0,N2-1):
  29             Q[(state(i,j),state(i-1,j+1))]= mu1
  30     print "ready filling"
  31 
  32 def computePiMethod1():
  33     e0 = time.time()
  34     Q = sp.dok_matrix((size,size)) 
  35     fillOffDiagonal(Q)
  36     # Set the diagonal of Q such that the row sums are zero
  37     Q.setdiag( -Q*ones(size) )
  38     # Compute a suitable stochastic matrix by means of uniformization
  39     l = min(Q.values())*1.001  # avoid periodicity, see trivedi's book
  40     P = sp.speye(size, size) - Q/l
  41     # compute Pi
  42     P =  P.tocsr()
  43     pi = zeros(size);  pi1 = zeros(size)
  44     pi[0] = 1;
  45     n = norm(pi - pi1,1); i = 0; 
  46     while n > 1e-3 and i < 1e5:
  47         pi1 = pi*P
  48         pi = pi1*P   # avoid copying pi1 to pi
  49         n = norm(pi - pi1,1); i += 1
  50     print "Method 1: ", time.time() - e0, i
  51     return pi
  52 
  53 def computePiMethod2():
  54     e0 = time.time()
  55     Q = pysparse.spmatrix.ll_mat(size,size)
  56     fillOffDiagonal(Q)
  57     # fill diagonal
  58     x =  empty(size)
  59     Q.matvec(ones(size),x)
  60     Q.put(-x)
  61     # uniformize 
  62     l = min(Q.values())*1.001
  63     P = pysparse.spmatrix.ll_mat(size,size)
  64     P.put(ones(size))
  65     P.shift(-1./l, Q)
  66     # Compute pi
  67     P = P.to_csr()
  68     pi = zeros(size);  pi1 = zeros(size)
  69     pi[0] = 1;
  70     n = norm(pi - pi1,1); i = 0; 
  71     while n > 1e-3 and i < 1e5:
  72         P.matvec_transp(pi,pi1)
  73         P.matvec_transp(pi1,pi) 
  74         n = norm(pi - pi1,1); i += 1
  75     print "Method 2: ", time.time() - e0, i
  76     return pi
  77 
  78 def plotPi(pi):
  79     pi = pi.reshape(N2,N1)
  80     matshow(pi)
  81     savefig("pi.png")
  82 
  83 if __name__ == "__main__":
  84     pi = computePiMethod1()
  85     pi = computePiMethod2()
  86     plotPi(pi)