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

with(linalg):
with(plots):

PhPhInfNetMoments  := proc(A, lambda, B, mu, K, Pr, nodelist, lower, upper)
#
# this procedure returns a plot of the mean and variance
# of number in at each node for the [Ph_t/Ph_t/Infinity]^K network as given
# in Nelson and Taaffe (2000) 
# 
# initial conditions for the network are empty and idle
#
# this version is designed to work with Maple V release 5
#              
# last update 8/1/00
#
# Inputs:
# A = Markov chain transition matrix for the arrival process
# lambda = arrival rate vector
# B = Markov chain transition matrix for the service process
#     [vector of matrices of dimension K]
# mu = service rate vector
#     [vector of vectors of dimension K]
# K = number of nodes in the network
# Pr = Markov routing matrix
# nodelist = list of nodes whose performance is to be plotted
# 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
#
#    [1] = mean plot
#    [2] = variance plot


           local  mA, mB, N, NN, k, i, j, g1, a, b, c, d, mdeqs, v, off, meanplots,
                  varianceplots, r, M, V;
#
# create index mapping function
           g1 := proc(i,j)
              RETURN(N.i.N.j);
           end;

# first check if there is only a single node
	if K = 1 then
	     RETURN(PhPhInfMoments(A, lambda, B, mu, lower, upper));
	else

# else it is a network, so use the network code
#
# obtain procedure to evaluate moment differential equations
          mdeqs := PhPhInfNetNumeric(A, lambda, B, mu, K, Pr):
#
# determine number of phases in phase representations
# and offsets for nodes
          mA := rowdim(A) -1;
	  v := vector(K);
	  off := vector(K+1);
	  v[1] := rowdim(B[1]) - 1;
	  off[1] := 0;
	  mB := v[1];
	  for r from 2 to K do
	     v[r] := rowdim(B[r]) - 1;
	     off[r] := off[r-1] + v[r-1];
	  od;
	  mB := off[K] + v[K];
	  off[K+1] := mB;

# 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 list of desired plots
	  meanplots := [t, sum(N[a,1], a=1..mB)];
	  varianceplots := [t, sum( sum( sum(NN[b][c,d], b=1..mA) -
                               N[c,1]*N[d,1], c=1..mB), d=1..mB) ];
	  for i in nodelist do
	     meanplots := meanplots, [ t, sum(N[a,1], a=off[i]+1..off[i+1]) ];

	     varianceplots := varianceplots, 
                       [t, sum( sum( sum(NN[b][c,d], b=1..mA) -
                          N[c,1]*N[d,1], c=off[i]+1..off[i+1]), d=off[i]+1..off[i+1]) ];
	  od;

# create the mean and variance plot
         M := odeplot(mdeqs, [ meanplots ], lower..upper, labels=[t, EN], color=black):
         V := odeplot(mdeqs, [ varianceplots ], lower..upper, labels=[t, VarN], color=black):
         RETURN(M, V);


	fi;
        end;