# copyright Barry L. Nelson and Michael R. Taaffe 1999

with(linalg):
with(plots):

MapPhInfMoments  := proc(A, vA, lambda, B, mu, lower, upper)
#
# this procedure returns a plot of the mean and variance
# of number in queue for the Map_t/Ph_t/Infinity queue as given
# in Nelson and Taaffe (2000c) 
# 
# initial conditions for the queue are empty and idle
#
# this version is designed to work with Maple V release 5
#              
# last update 11/1/00  7/29/04
#
# Inputs:
# A = Markov chain transition matrix for the arrival process
# vA = number of absorbing states in Map representation
# lambda = arrival rate vector
# B = Markov chain transition matrix for the service process
# mu = service rate vector
# mdeqs = numerical procedure for moment differential equations
# lower = lower time limit for plot
# upper = upper time limit for plot
#
#
# Output:
# odeplot of mean and variance of number in queue for
#    lower < t < upper

           local  mA, mB, N, NN, k, i, j, g1, a, b, c, d, mdeqs, M, V;
#
# create index mapping function
           g1 := proc(i,j)
              RETURN(N||i||N||j);
           end;
#
# obtain procedure to evaluate moment differential equations
          mdeqs := MapPhInfNumericEmpty(A, vA, lambda, B, mu):
#
# determine number of phases in phase representations
          mA := rowdim(A) -vA;
          mB := rowdim(B) -1;
# create the N vector
          N := matrix(mB, 1, [N||(1..mB)]);
          N := N(t);
# create the NN  matrix
          NN := vector(mA);
          for k from 1 to mA do
             NN[k] := map(cat, matrix(mB,mB,g1), "A"||k)(t); 
          od;
# create the mean and variance plot
          M := odeplot(mdeqs, [t, sum(N[a,1], a=1..mB)], lower..upper, labels=[t, EN], color=black);
          V := odeplot(mdeqs,[t, sum( sum( sum(NN[b][c,d], b=1..mA) - N[c,1]*N[d,1], c=1..mB), d=1..mB)], lower..upper, labels=[t, VarN], color=black);
          RETURN(M, V);
         end;