%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    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: [B0,B1,B2] =  infinite2finite(A0,A1,A2,n,k,epsilon,PH)
%
% Given a QBD process with infinite levels,
% it returns an approximate QBD process with finite levels.
% Here, the generator matrix of the input 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.
%
% The generator matrix of the output QBD process is
%
% [ B1{1} B0{1}                                          ]
% [ B2{2} B1{2} B0{2}                                    ]
% [       B2{3} B1{3} B0{3}                              ]
% [              :     :                                 ]
% [                        B2{k} B1{k}   B0{k}           ]
% [                              B2{k+1} B1{k+1} B0{k+1} ]
%
%
% 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}
%        n: level that the QBD process starts repeating
%        k: the number of levels (k+1) in the output QBD process
%        epsilon: the error bound (epsilon = 1e-6: recommended).
%        PH: number of phases used in approximate PH distributions (1 or 2)
%
% OUTPUT: B0,B1,B2: cells for B0{i}, B1{i}, B2{i}
%
% author: Takayuki Osogami
%         Department of Computer Science
%         Carnegie Mellon University
%         osogami@cs.cmu.edu
% date: August 15, 2004

function [B0,B1,B2] =  infinite2finite(A0,A1,A2,n,k,epsilon,PH)

%%%%%%%%%%%%%%%%%%%%%
%% "non-busy" part %%
%%%%%%%%%%%%%%%%%%%%%

for i = 1:k
  if( i <= n )
    B1{i} = A1{i};
    B2{i} = A2{i};
  else
    B1{i} = A1{n};
    B2{i} = A2{n};
  end
end
for i = 1:k-1
  if( i <= n )
    B0{i} = A0{i};
  else
    B0{i} = A0{n};
  end
end


%%%%%%%%%%%%%%%%%
%% "busy" part %%
%%%%%%%%%%%%%%%%%

%% calculating moments of busy periods %%

[R,U,G] = RUGmatrices(A0,A1,A2,n,epsilon);

if( k + 1 >= n )
  GG = G{n};
else
  GG = G{k+1};
end

FF = A0{k};
dim = size(FF);
N = dim(1);
for i = 1:N
  total = sum(FF(i,:));
  if( total > 0 )
    FF(i,:) = FF(i,:) / total;
  end
end

PP = FF * GG;

DD = A0{k} * GG;

HH = Neuts3(A0,A1,A2,G,n,k+1,epsilon);

for h = 1:3

 EE{h} = A0{k} * HH{h};

 for i = 1:N
  for j = 1:N
   if( DD(i,j) > 0 )
    MM{h}(i,j) = EE{h}(i,j) / DD(i,j);
   else
    MM{h}(i,j) = 0;
   end
  end
 end
end


%% approximating busy periods by PH distributions %%

if( PH == 2 ) % matching 3 moments
  for i = 1:N
   for j = 1:N
    if( MM{1}(i,j) > 0 )
     [tau{i,j},T{i,j}] = coxian2(MM{1}(i,j),MM{2}(i,j),MM{3}(i,j));
     tau{i,j} = tau{i,j} * PP(i,j);
     dim = size(T{i,j});
     t{i,j} = - T{i,j} * ones(dim(1),1);
    end
   end
  end
elseif( PH == 1) % matching mean (1 moment)
  for i = 1:N
   for j = 1:N
    if( MM{1}(i,j) > 0 )
     T{i,j} = - 1 / MM{1}(i,j);
     tau{i,j} = 1 * PP(i,j);
     t{i,j} = - T{i,j};
    end
   end
  end
else
 not_supported = 1
end


%% specifying B0{k}, B1{k+1}, B2{k+1} %%
 
base = 0;
for i = 1:N
  for j = 1:N
   if( MM{1}(i,j) > 0 )
    B1{k+1}(base+1:base+PH,base+1:base+PH) = T{i,j};
    base = base + PH;
   end
  end
end
    
dim0 = size(B1{k});
dim1 = size(B1{k+1});
B0{k} = zeros(dim0(1),dim1(2));
B2{k+1} = zeros(dim1(1),dim0(2));

base = 0;
for i = 1:N
  for j = 1:N
   if( MM{1}(i,j) > 0 )
    B2{k+1}(base(1)+1:base(1)+PH,j) = t{i,j};
    base = base + PH;
   end
  end
end

dim = size(A0{k});
lambda = A0{k} * ones(dim(2),1);

base = 0;
for i = 1:N
  for j = 1:N
   if( MM{1}(i,j) > 0 )
    B0{k}(i,base+1:base+PH) = lambda(i) * tau{i,j};
    base = base + PH;
   end
  end
end