%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Copyright (c) 2004.  Takayuki Osogami.  All Rights Reserved.   %
%                                                                   %
% Permission to use, copy, modify, and distribute this source code  %
% and its documentation for any purpose, without fee, and without   %
% a written agreement, is hereby granted, provided that the above   %
% copyright notice, this paragraph and the following two paragraphs %
% appear in all copies, modifications, and distributions.           %
% Created by Takayuki Osogami, Department of Computer Science,      %
% Carnegie Mellon University.                                       %
%                                                                   %
% IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR ANY       %
% DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS   %
% SOURCE CODE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN    %
% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                        %
%                                                                   %
% THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING,      %
% BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY     %
% AND FITNESS FOR A PARTICULAR PURPOSE.  THE SOURCE CODE AND        %
% ACCOMPANYING DOCUMENTATION, IF ANY, PROVIDED HEREUNDER IS         %
% PROVIDED "AS IS."  THE AUTHOR HAS NO OBLIGATION TO PROVIDE        %
% MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Usage: [MRI] = MRImatrix(M,r,i,G)
%
% Theorem 2 in the following paper:
%
% Marchel F. Neuts, "Moment formulas for the Markov renewal branching
% processes," Advances in Applied Probabilities, Vol. 8, pp. 690-711 (1978)
%
% MRI is the r-th moment of the busy period started by i jobs.

function [MRI] = MRImatrix(M,r,i,G)

S = size(M{1});
dim = S(1);
MRI = zeros(dim,dim);

N = zeros(dim,dim,4);
N(:,:,1) = G;
for j = 1:r
 N(:,:,j+1) = M{j};
end

if( i == 1)

 MRI = M{r};

else

 if( r == 1 )

  for j = 0:i-1
   MRI = MRI + G^j * M{1} * G^(i-1-j);
  end

 elseif( r == 2 )

  if( i == 2 )
 
   for k1 = 0:r
    k2 = r - k1;
    MRI = MRI + factorial(r) / (factorial(k1) * factorial(k2)) * N(:,:,k1+1) * N(:,:,k2+1);
   end

  elseif( i == 3 )

   for k1 = 0:r
    for k2 = 0:r-k1
     k3 = r - k1 - k2;
     MRI = MRI + factorial(r) / (factorial(k1) * factorial(k2) * factorial(k3)) * N(:,:,k1+1) * N(:,:,k2+1) * N(:,:,k3+1);
    end
   end

  end

 end 
end