#!/usr/local/bin/python
#
# Slotted p-persistent CSMA Medium Access
#
# Christof Meerwald <cmeerw@web.de>
# http://cmeerw.cjb.net/study/mobil-3/
#
# $RCSfile: ppers-markov.py,v $ - $Author: cmeerw $
# $Revision: 1.6 $ - $Date: 2000/06/25 15:02:44 $
#
import string, sys


# associative array to map state names to indices
index = {
    "1":   0,
    "2":   1,
    "3":   2,
    "D1":  3,
    "D1a": 4,
    "D1b": 5,
    "S1":  6,
    "D2":  7,
    "D2a": 8,
    "D2b": 9,
    "S2":  10,
    "C2":  11,
    "D3":  12,
    "D3a": 13,
    "D3b": 14,
    "S3":  15,
    "C3":  16,
    "END": 17
}


##
# generate a null n*n matrix
#
# @param n matrix width/length
# @return generated matrix
##
def generate_matrix(n):
    m = []
    row = []
    for i in range(n):
        row.append(0)
    for j in range(n):
        m.append(row[:])

    return m


##
# matrix multiplication of n*n matrices
#
# @param x matrix
# @param y matrix
# @return result matrix
##
def matrix_mul(x, y):
    z = []
    range_n = range(len(x))

    for i in range_n:
        row = []
        for j in range_n:
            sum = 0.0
            for k in range_n:
                sum = sum + x[i][k]*y[k][j]
            row.append(sum)

        z.append(row)

    return z


##
# set the probability of a transition
#
# @param matrix matrix to modifiy
# @param src source node
# @param dst destination node
# @param p transition probability
##
def transition(matrix, src, dst, p):
    matrix[index[src]][index[dst]] = p


##
# Print program usage information
##
def usage():
    print "Usage: ppers-markov.py [p [pn [n]]]"
    print "       p .. persistence value"
    print "       pn .. arrival probability"
    print "       n .. number of iterations"


##
# Main program.
#
# @param args parameter list
# @return return value
##
def main(args):
    # default parameters
    p = 0.35
    pn = 0.1
    n = 40

    if len(args) > 3:
        usage()
        return 1

    if len(args) >= 1:
        try:
            p = string.atof(args[0])
        except ValueError:
            usage()
            return 1

    if len(args) >= 2:
        try:
            pn = string.atof(args[1])
        except ValueError:
            usage()
            return 1

    if len(args) >= 3:
        try:
            n = string.atoi(args[2])
        except ValueError:
            usage()
            return 1


    # generate a null n*n matrix (n = len(index))
    matrix = generate_matrix(len(index))


    # set the transition probabilities for the matrix
    transition(matrix,  "1",   "1",  (1-p)*(1-pn)**2)
    transition(matrix,  "1",   "2",  2*(1-p) * pn*(1-pn))
    transition(matrix,  "1",   "3",  (1-p) * pn**2)
    transition(matrix,  "1",  "D1",  p)
    transition(matrix,  "D1", "D1a", 1)
    transition(matrix, "D1a", "D1b", 1)
    transition(matrix, "D1b", "S1",  1)

    transition(matrix,  "2",   "2",  (1-p)**2 * (1-pn))
    transition(matrix,  "2",   "3",  (1-p)**2 * pn)
    transition(matrix,  "2",  "D1",  p*(1-p))
    transition(matrix,  "2",  "D2",  p*(1-p))
    transition(matrix,  "2",  "C2",  p*p)
    transition(matrix, "D2",  "D2a", 1)
    transition(matrix, "D2a", "D2b", 1)
    transition(matrix, "D2b", "S2",  1)
    transition(matrix, "S2",   "1",  (1-pn)**6)
    transition(matrix, "S2",   "2",  pn*(1-pn)**5 + 1 - pn - (1-pn)**6)
    transition(matrix, "S2",   "3",  pn - pn*(1-pn)**5)
    transition(matrix, "C2",   "2",  (1-pn)**2)
    transition(matrix, "C2",   "3",  2*pn - pn*pn)

    transition(matrix,  "3",   "3",  (1-p)**3)
    transition(matrix,  "3",  "D1",  p*(1-p)**2)
    transition(matrix,  "3",  "D3",  2*p*(1-p)**2)
    transition(matrix,  "3",  "C3",  3*p*p - 2*p**3)
    transition(matrix, "D3",  "D3a", 1)
    transition(matrix, "D3a", "D3b", 1)
    transition(matrix, "D3b", "S3",  1)
    transition(matrix, "S3",   "2",  1-pn)
    transition(matrix, "S3",   "3",  pn)
    transition(matrix, "C3",   "3",  1)
    
    transition(matrix, "S1",  "END", 1)
    transition(matrix, "END", "END", 1)


    # initialize or working matrix
    m = matrix

    # Markov-model iteration
    for i in range(1, n):
        # print the probabilities for successful completion after i steps
        # for 1, 2, and 3 initial stations
        print '%3d: %.4f  %.4f  %.4f' % (i,
                                         m[index["1"]][index["END"]],
                                         m[index["2"]][index["END"]],
                                         m[index["3"]][index["END"]])

        # prepare for the next iteration
        m = matrix_mul(m, matrix)


    print '%3d: %.4f  %.4f  %.4f' % (n,
                                     m[index["1"]][index["END"]],
                                     m[index["2"]][index["END"]],
                                     m[index["3"]][index["END"]])

    return 0


main(sys.argv[1:])