inputFileName='4Client2Freqc.aut';
fid = fopen( inputFileName,'rt');
if (fid < 0)
    error('could not open file' + inputFileName);
end;
line = strrep( fgetl(fid), 'des (', '');
line = strrep( line, ',', ' ');
systemSize = sscanf( line, '%d'); % extracting the transition system size
display( systemSize);
stateNum = systemSize(3); % number of states in the model
display(stateNum);
Q = sparse( stateNum, stateNum);
b = sparse( stateNum, 1);
reward = ones( 1, stateNum);
% variables for the rates
x = 0;
y = 0;
while ~feof(fid) 
    remain = fgetl(fid);
    % display(remain);
    % extracting the transition information from the file
    [start,remain] = strtok( remain, ',)( "');
    [mcrl_action,remain] = strtok( remain, ',)( "');
    [pepa_action,remain] = strtok( remain, ',)( "');
    [type,remain] = strtok( remain, ',)( "');
    [rate,remain] = strtok( remain, ',)( "');
    [finish,remain] = strtok( remain, ',)( "');    
    dim1 =  str2num(start) + 1;
    dim2 =  str2num(finish) + 1;
    if( strcmp( pepa_action, 'p_col1')) 
        reward(dim1)=0;    %rewarding zero for collision
    end
    if (strcmp( rate, 'x'))
        Q( dim1, dim2) = -10;
    elseif (strcmp( rate, 'y'))
        Q( dim1, dim2) = -20;
    else
        Q( dim1, dim2) = str2double( rate);
    end
end;
Q_orig = Q;
for x = 1:20
    for y = 1:20
        Q = Q_orig;
        [m,n] = size(Q);
        for i = 1:m
            for j = 1:n
                if( Q(i,j) == -10) Q(i,j) = x;
                    elseif( Q(i,j) == -20) Q(i,j) = y;
                end
            end
        end
        %display(Q);
                    
        
        
        rowsSum = sum(Q,2);
        %display(rowsSum);
        for i=1:stateNum,
            Q(i,i)=-rowsSum(i);
        end;
        Qt = Q';                        % transpose the generator matrix
        for i=1:stateNum                % replace the last equation by the normalisation constant
          Qt(stateNum,i) = 1.0;                       
        end;					  
        b(stateNum) = 1.0;    
        p = Qt\b;                       % solve to find the steady state distribution
        %display(det(Qt));
        %display( reward);
        display('*******************************');
        display( x);
        display( y);
        display(reward * p);
        Result( x, y) = reward * p;
        display('*******************************');        
    end
end
display( Result);
fclose(fid);