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

with(linalg):
with(plots):

PhPhInfNetSojourn  := proc(A, lambda, B, mu, K, Pr, arrivalT)
#
# this procedure returns the mean, variance and a plot structure of
# the cdf of the sojourn time of a virtual arrival at time arrivalT
# to the [Ph_t/Ph_t/Infinity]^K network as given in Nelson and 
# Taaffe (2000)
#
# last update 7/28/00 08/20/01 10/15/03
# 
# 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
# arrivalT = arrival time of virtual arrival
#
#
# Output:
# [1] mean of sojourn time
# [2] variance of sojourn time
# [3] plot structure for the cdf of sojourn time
#
           local  mA, mB, N, k, i, j, g1, a, c, d, b, x,
                  InitP, InitN, InitNA, InitNN, mdeq, cdf, first,
                  second, Tinf, BB, Bmu, offi, offj, v, g, h, r;
#
# create index mapping function
           g1 := proc(i,j)
              RETURN(N||i||N||j);
           end;

# check if there is only a single node
	if K = 1 then
	   Return(PhPhInfSojourn(A, lambda, B, mu, arrivalT));

	else

# else it is a network, so use the network code
#
# determine number of phases in phase representations
# and offsets for nodes
          mA := rowdim(A) -1;


# first extract key dimensions
	   v := vector(K);
	   v[1] := rowdim(B[1]) - 1;
	   mB := v[1];
	   for r from 2 to K do
	      v[r] := rowdim(B[r]) - 1;
	      mB := mB + v[r];
	   od;

	   BB := matrix(mB+1, mB+1, 0);

# build the composite service-rate vector
	   Bmu := vector(mB);
           offi := 0;
           for i from 1 to K do
	      for g from 1 to v[i] do
	         Bmu[offi + g] := mu[i][g];
 	      od;
	      offi := offi + v[i]
	   od; 

# build the composite phase transition matrix
	   offi := 0; offj := 0;

#      build rows and columns 1 through mB
	   for i from 1 to K do
	      for j from 1 to K do
	         for g from 1 to v[i] do
	            for h from 1 to v[j] do
	               BB[offi + g, offj + h] := Pr[i,j]*B[i][g, v[i]+1]*B[j][v[j]+1,h];
	            od;
	         od;
	         offj := offj + v[j];
              od;
	      offi := offi + v[i];
	      offj := 0;
           od;

#      build the last row
	   offj := 0;
	   for j from 1 to K do
	      for h from 1 to v[j] do
	         BB[mB+1, offj + h] := Pr[K+1,j]*B[j][v[j]+1,h];
              od;
              offj := offj + v[j];
           od;

#      build the last column
	   offi := 0;
	   for i from 1 to K do
	      for g from 1 to v[i] do
	         BB[offi + g, mB+1] := Pr[i, K+1]*B[i][g, v[i]+1];
              od;
	      offi := offi + v[i];
           od;

	   BB[mB+1, mB+1] := 0;

#      add the terms representing immediate feedback
	   BB := matadd(BB, MatDiag(B, K));


# extract beta
          b := transpose(submatrix(BB, mB+1..mB+1, 1..mB));
# create the initial conditions
          InitP := vector(mA, 0);
          InitN := vector(mB);
          for i from 1 to mB do
             InitN[i] := b[i,1];
          od;
          InitNA := matrix(mA, mB, 0);
          InitNN := vector(mA);
          for i from 1 to mA do
             InitNN[i] := matrix(mB, mB, 0);
          od;

# create the numerical procedure
          mdeq := PhPhInfNumericInit(A, lambda, BB, Bmu, arrivalT, InitP, InitN, InitNA, InitNN); 
#
# create the N vector for plotting
          N := matrix(mB, 1, [N||(1..mB)]);
          N := N(t);
#
# compute first and second moment
         cdf := subs(mdeq, 1 - sum(N[a,1], a=1..mB)):
         first := subs(mdeq, M1(t));
         second := subs(mdeq, M2(t));
         Tinf := arrivalT;
         while cdf(Tinf) < 0.99999 do
          Tinf := Tinf + 0.1
         od;
         RETURN(first(Tinf), second(Tinf)- first(Tinf)^2, 
         odeplot(mdeq, [t, 1 - sum(N[a,1], a=1..mB)], arrivalT..arrivalT + Tinf));

	fi;

         end;



MatDiag := proc(x, k)

# auxiliary proc to build a matrix with other matrices on the diagonal

	local m, lower, i, j, l, d, M;

	m := rowdim(x[1]) - 1;
  	for l from 2 to k do
	  m := m + rowdim(x[l]) - 1;
	od;
	M := matrix(m+1, m+1, 0);

	lower := 0;

	for l from 1 to k do
	   d := rowdim(x[l]) - 1;
	   for i from 1 to d do
	      for j from 1 to d do
	         M[lower + i, lower + j] := x[l][i,j];
              od;
	   od;
	   lower := lower + d;
	od;
	RETURN(M);
	end;