varargout=compute_blocking_kovacevic(varargin)

function varargout=compute_blocking_kovacevic(varargin)
%COMPUTE_BLOCKING_KOVACEVIC Compute blocking propability in a circuit switched
%network using a generalized reduced load approximation simplified by
%Kovacevi and Acampora [5].
%   [BLOCKING]=COMPUTE_BLOCKING_BIRMAN(NAME,NUMBER,EPSILON,METHOD,D,C) is a
%   function which computes blocking probability using generalized reduced
%   load approximation method (RLA) first introduced by Kelly [4], extended
%   by Chung et al., implemented for WDM network by Birman [1] (full and no
%   wavelength conversion) and Tripathi [2] (sparse wavelength conversion).
%   BLOCKING is given in % and simplified by Kovacevi and Acampora [5].
%
%   Because normal RLA requires O(C) operations for circuit-switched case,
%   and O(C^(n+1)) operations for wavelength routed cases where length of
%   path is 2 it is worth to use simple models which may give not so prcise
%   solutions but can be much faster than normal RLA. 
%
%   NAME is a string containing name of analized network saved in
%   MATLABROOT/networks directory, NUMBER is the number of dimension in
%   TRAFFIC structure. This means, that for analisys TRAFFIC{:,:,NUMBER}
%   will be choosed.
%
%   METHOD is the type of conversion and D is the degree of sparse
%   conversion (for details see help P_M).
%
%   EPSILON is some treshold value given by user. For details see [1], [2].
%   C is the channel number which we want to use to evaluate blocking.
%
%   See also COMPUTE_BLOCKING_BIRMAN, P_M, CREATE_NETWORK,
%   CREATE_LOAD_MANUAL, LINK_SIM_CONV.
%   
%   References:
%       [1] Alexander Birman, "Computing Approximate Blocking Probabilities
%       for a Class of All-Optical Network", IEEE Journal on
%       Selected Areas in Communications, vol. 14, June 1996 
%       [2] Tushar Tripathi, Kumar N. Sivarajan, "Computing Approximate
%       Blocking Probabilities in Wavelength Routed All_optical Networks
%       with Limited-Range Wavelength Conversion", IEEE Journal on
%       Selected Areas in Communications, vol. 18, October 2000 
%       [3] S. Chung, A. Kashper, K.W. Ross, "Computing Approximate
%       Blocking Propabilities for large loss networks", IEEE/ACM
%       Transactions on Networking, vol. 1, 1993
%       [4] F.P. Kelly, "Routing and capcity allocation in networks with
%       trunk reservation", Matematics of Operation Research, vol. 15, no.
%       4, 1990 
%       [5] Milan Kovacevi, Anthony Acampora, "Benefits of Wavelength
%       translation in All-Optical Clear Channel Networks", IEEE Journal on 
%       Selected Areas in Communications, vol. 14, June 1996 

%   Copyright 2003-2004 Przemyslaw Pawelczak
%   E-mail: [e-mail]
%   Wroclaw University of Technology, Poland
%   $Revision: 1.0 $  $Date: 2004/05/17 17:21:32 $

%Exeption handling
message=nargchk(6,6,nargin);
if ~isempty(message)
    error('MATLAB:COMPUTE_BLOCKING_KOVACEVIC:NumberOfInputArguments',...
        message);
end
testing_int=[varargin{2},varargin{3},varargin{4},varargin{5},varargin{6}];
if find(isnan(testing_int))
    error('MATLAB:COMPUTE_BLOCKING_KOVACEVIC:ArgumentType',...
        'Argument must be numbers.');
end
%Check if VARARGIN are positive integers
if sum([testing_int<0,...
            fix([testing_int(1),testing_int(3:5)])~=...
            [testing_int(1),testing_int(3:5)]])~=0
    error('MATLAB:COMPUTE_BLOCKING_KOVACEVIC:ArgumentType',...
        'Wrong type of argument.');
end
if ~isstr(varargin{1})
    error('MATLAB:COMPUTE_BLOCKING_KOVACEVIC:ArgumentType',...
        'Argument must be string.');
end

%Change string into number (for all number arguments)
for k=2:5
    if isstr(varargin{k})==1
        varargin{2}=str2num(varargin{2});
    end
end

%List all '*.mat' files in current working directory
%Change directory to MATLABROOT\networks - catalog must exist
cd([matlabroot,['\work\networks']]);
directory=what;
directory=struct2cell(directory);
directory=cellstr(directory{3});

%Check if specified '*.mat' file already exist
if find(strcmp(directory,[varargin{1},'.mat'])==1)
    %Load *.MAT file
    load(varargin{1});
    fields_vector=[exist('description'),exist('structure'),exist('traffic')];
    if length(find(fields_vector))~=3
        error('MATLAB:COMPUTE_BLOCKING_KOVACEVIC:BadFieldNames',...
            'Bad field names.'); 
    end
    %Check if BLOCKING definition already exists in this file
    %Initialize flag
    blocking_exists=0;
    varias=who;
    varias_dimension=cellfun('ndims',varias);
    varias_dimension=length(varias_dimension');
    %Check
    for k=1:varias_dimension
        if isequal(cell2mat(varias(k)),'blocking')
            blocking_exists=1;
        end
    end
    %If definition does exist - add another
    if blocking_exists==1
        %Check how many BLOCKING definitions are there
        blocking_dimension=cellfun('ndims',blocking);
        blocking_dimension=length(blocking_dimension)+1;
        %Load BLOCKING procedure
        blocking=kovacevic(traffic,varargin{2},varargin{3},varargin{4},...
            varargin{5},varargin{6},blocking,blocking_dimension);
        varargout{1}=save_network(description,structure,traffic,blocking,...
            varargin{1});
        %If definition does not exist - make first one
    elseif blocking_exists==0;
        blocking{1,1}=[];
        blocking_dimension=1;
        %Load MAKE_LOAD procedure
        blocking=kovacevic(traffic,varargin{2},varargin{3},varargin{4},...
            varargin{5},varargin{6},blocking,blocking_dimension);
        varargout{1}=save_network(description,structure,traffic,blocking,...
            varargin{1});
    end 
else
    error('MATLAB:COMPUTE_BLOCKING_KOVACEVIC:BadFiledName',...
        'Bad file name.'); 
end


function [blocking]=kovacevic(traffic,number,epsilon,method,d,C,blocking,blocking_dimension);

%Check how many routes does TRAFFIC{:,:,NUMBER} have
%Check how many nodes are in the TRAFFIC{:,:,NUMBER}
traffic_dim=cellfun('ndims',traffic);
traffic_length=length(traffic_dim(:,1,1));
route_number=0;
for k=1:traffic_length
    if ~isempty(traffic{k,1,number})
        route_number=route_number+1;
        max_node_number(k)=[max(traffic{k,1,number})];
    end
end
node_number=max(max_node_number);

%Name links in the network with numbers
nodes_matrix=zeros(node_number,node_number);
node_names(1:(node_number*(node_number-1)/2))=1:(node_number*(node_number-1)/2);
for k=1:(node_number-1)
    places=(1+(k-1)*node_number+(k-1)+1):(node_number*k);
    nodes_matrix(places)=node_names(1:(node_number-k));
    node_names(1:(node_number-k))=[];
end
nodes_matrix=nodes_matrix+nodes_matrix';

%Check which links exist in TRAFFIC definition and change routes from edge
%adjacency to link adjacency
%Memory initialization
link_routes{1}=[];
for k=1:route_number
    route=traffic{k,1,number};
    %Initialization (vector where links of the route will be stored)
    link_route=[];
    for w=1:(length(route)-1)
        %Find name of the link in NODES_MATRIX
        link_route=[link_route,nodes_matrix(route(w),route(w+1))];
        %One link can traverse path for only one time
        link_route=unique(link_route);
        %Check how long current route is
        route_length(k)=length(link_route);
    end
    link_routes{k}=link_route;
end
%Find which path is the longest (if >1 longest paths - choose one randomly)
%PATH_LENGTH_POSITION will be lately used in main loop
[path_length,path_length_position]=min(route_length);

%Check which routes use specified links
%Memory initialization
for k=1:(node_number*(node_number-1)/2)
    which_routes_use_link{k}=[];
end
for k=1:route_number
    path=link_routes{k};
    for w=1:length(path)
        which_routes_use_link{path(w)}=[which_routes_use_link{path(w)},k];
    end
end

%Lr and Lr_prim initialization
Pl(1:route_number)=0;
Pl_prim(1:route_number)=Inf;
%Ls initialization (offered load on link s)
Ls(1:(node_number*(node_number-1)/2))=0;
%Bs initialization (blocking on link s)
Bs(1:(node_number*(node_number-1)/2))=0;

%Main loop
while abs(Pl(path_length_position)-Pl_prim(path_length_position))>epsilon
    %Assign Pl_prim to Pl
    Pl_prim=Pl;
    
    %Compute Ls
    for s=1:(node_number*(node_number-1)/2)
        if ~isempty(which_routes_use_link{s})
            Al_vector=[];
            Pl_vector=[];
            for g=1:length(which_routes_use_link{s})
                Al_vector=[Al_vector,traffic{which_routes_use_link{s}(g),2,number}];
                Pl_vector=[Pl_vector,1-Pl(which_routes_use_link{s}(g))];
            end
            %Equation (7) in [5]
            Ls(s)=sum(Al_vector.*Pl_vector)/(1-Bs(s));
        end
    end
    
    %Compute Bs
    for s=1:(node_number*(node_number-1)/2)
        if ~isempty(which_routes_use_link{s})
            Bs(s)=erlang_b(Ls(s),C,C);
        end
    end

    %Compute Pl
    for k=1:route_number
        L_vector=[];
        for w=1:length(link_routes{k})
            L_vector=[L_vector,Ls(link_routes{k}(w))];
        end
        Pl(k)=p_l(L_vector,C,C,d,method);
    end
    
end

%Blocking in %
Pl=Pl*100;
%Assign blocking to VARARGOUT of function
blocking{1,blocking_dimension}=Pl;


%Saving to structure
function [net_out]=save_network(description,structure,traffic,blocking,name);
save(name,'description','structure','traffic','blocking');
%Cannot make structure with cells of different dimension (sic!) - it needs
%checking
net_out=struct('blocking',blocking);