%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    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: [pi] = pivectors(A0,A1,A2,R,U,n)
%
% This program computes the pi vectors for a QBD process.
% Here, the generator matrix of the QBD process is
%
% [ A1{1} A0{1}                                    ]
% [ A2{2} A1{2} A0{2}                              ]
% [       A2{3} A1{3} A0{3}                        ]
% [              :     :                           ]
% [                        A2{n} A1{n} A0{n}       ]
% [                              A2{n} A1{n} A0{n} ]
% [                                     :     :    ]
%
% Note that the QBD process repeats after level n.
%
% We define cells, A0, A1, and A2:
% A0 = {A0{1} A0{2} ... A0{n}}
% A1 = {A1{1} A1{2} ... A1{n}}
% A2 = {A2{1} A2{2} ... A2{n}},
%
% INPUT: A0,A1,A2: cells for A0{i}, A1{i}, A2{i}
%        R: cells for the R matrices (use RUGmatrices.m to get these)
%        U: cells for the U matrices (use RUGmatrices.m to get these)
%        n: level that the QBD process starts repeating
%
% OUTPUT: pi: cells for pi{i}
%             pi{i} = pi vector for level i for i <= n
%             For i > n, pi{i} can be obtained by pi{i} = pi{i-1} * R{n}
%
%
% author: Takayuki Osogami
%         Department of Computer Science
%         Carnegie Mellon University
%         osogami@cs.cmu.edu
% date: April 28, 2004

function [pi,meannum] = pivectors(A0,A1,A2,R,U0,n)

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % get pi vectors without normalization %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  dim = size(U0);
  M = - U0;
  M(:,dim(1)) = ones(dim(1),1);
  vec = zeros(1,dim(1));
  vec(1,dim(1)) = 1;
  pi{1} = vec * inv(M);

  for i = 2:n;
    pi{i} = pi{i-1} * R{i};
  end


 %%%%%%%%%%%%%%%%%
 % normalization %
 %%%%%%%%%%%%%%%%%

  total = 0;
  meannum = 0;
  for i = 1:n-1
    p = pi{i} * ones(size(pi{i}))';
    total = total + p;
    meannum = meannum + p * (i-1);
  end

  Tmp = inv(eye(size(R{n}))-R{n});
  E = ones(size(pi{n}))';
  total = total + pi{n} * Tmp * E;
  meannum = meannum + (n-1) * pi{n} * Tmp * E + pi{n} * R{n} * Tmp^2 * E;

  for i = 1:n
    pi{i} = pi{i} / total;
  end
  meannum = meannum / total;