%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is an example template Matlab program.                           %
%                                                                       %
% It shows the standard steps involved in solving the generator matrix  %
% of the Markov process underlying a PEPA component.  It is assumed     %
% that the matrix is in the file "transitions.m".                       %
%                                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% First, how many states were there?  

n = 6841;

% Now assign values to the rate variables
% which appear in the PEPA program
LT=1;
PT=1;
DT=1;
IT=1;
RT=1;
ITR=1;
IR=1;
RR=1;
CT=1;
disp('Assigned values of variables');

% Now make Q a sparse all-zero matrix of the right size.

Q = sparse(n,n);
disp('Made the sparse matrix');

% Load in the transitions

transitions;
disp('Loaded the transitions');

% In order to incorporate the normalisation condition and make the 
% equations non-singular the final column of the matrix is replaced
% by a column of 1's

for i = 1:n,
    Q(i,n) = 1.0;
end
disp('Final column all 1s');

% we want to solve the matrix equation PQ = 0 where P is the steady
% state probability distribution and therefore subject to the 
% normalisation condition. We set up the vector b = 0 as follows.

b = zeros(n,1);

% the corresponding adjustment to b is to introduce a 1 in the row
% corresponding to the column of 1s in Q. 

b (n) = 1;
disp('Set up b');

% We will be solving the equation QT P = b where QT is the transpose
% of the modified Q:

QT = Q';
disp('Transposed Q');

% x = A\b is the MATLAB command to find x such that A*x = b 
% Thus in this case P will be the steady state probability distribution:

P = QT\b

% This is the end of the file