pfqn_ncld.m

%{
%{
 % @file pfqn_ncld.m
 % @brief Normalizing constant for load-dependent closed networks.
%}
%}
 
%{
%{
 % @brief Normalizing constant for load-dependent closed networks.
 % @fn pfqn_ncld(L, N, Z, mu, varargin)
 % @param L Service demand matrix.
 % @param N Population vector.
 % @param Z Think time vector.
 % @param mu Load-dependent rate matrix.
 % @param varargin Optional solver parameters.
 % @return lG Logarithm of normalizing constant.
 % @return G Normalizing constant.
 % @return method Method used for computation.
%}
%}
function [lG,G,method] = pfqn_ncld(L,N,Z,mu,varargin)
% [LGN,G,METHOD] = PFQN_NCLD(L,N,Z,VARARGIN)
 
options = Solver.parseOptions(varargin, SolverNC.defaultOptions);
lG = NaN;
G = NaN;
method = options.method;
 
% backup initial parameters
 
mu = mu(:,1:sum(N));
% first remove empty classes
nnzClasses = find(N);
L = L(:,nnzClasses);
N = N(:,nnzClasses);
Z = Z(:,nnzClasses);
 
% then scale demands in [0,1], importat that stays before the other
% simplications in case both D and Z are all very small or very large in a
% given class, in which case the may look to filter but not if all of them
% are at the same scale
R = length(N);
scalevec = ones(1,R);
for r=1:R
    scalevec(r) = max([L(:,r);Z(:,r)]);
end
L = L ./ repmat(scalevec,size(L,1),1);
Z = Z ./ scalevec;
 
% remove stations with no demand
Lsum = sum(L,2);
Lmax = max(L,[],2);
demStations = find((Lmax./Lsum)>GlobalConstants.FineTol);
L = L(demStations,:);
mu = mu(demStations,:);
 
% if there is a class with jobs but with L and Z all zero
if any(N((sum(L,1) + sum(Z,1)) == 0)>0)
    line_warning(mfilename,'The model has no positive demands in any class.\n');
    if isempty(Z) || sum(Z(:))<options.tol
        lG = 0;
    else
        lG = - sum(factln(N)) + sum(N.*log(sum(Z,1))) + N*log(scalevec)';
    end
    G = NaN;
 
    return
end
 
% update M and R
[M,R]=size(L);
 
% return immediately if the model is a degenerate case
if isempty(L) || sum(L(:))<options.tol % all demands are zero
    if isempty(Z) || sum(Z(:))<options.tol
        lG = 0;
    else
        lG = - sum(factln(N)) + sum(N.*log(sum(Z,1))) + N*log(scalevec)';
    end
    return
elseif M==1 && (isempty(Z) || sum(Z(:))<options.tol) % single node and no think time
    lG = factln(sum(N)) - sum(factln(N)) + sum(N.*log(sum(L,1))) + N*log(scalevec)' - sum(log(mu(1:sum(N))));
    return
end
 
% determine contribution from jobs that permanently loop at delay
zeroDemandClasses = find(sum(L,1)<options.tol); % all jobs in delay
nonzeroDemandClasses = setdiff(1:R, zeroDemandClasses);
 
if isempty(sum(Z,1)) || all(sum(Z(:,zeroDemandClasses),1)<options.tol)
    lGzdem = 0;
    Nz = 0;
else
    if isempty(zeroDemandClasses) % for old MATLAB release compatibility
        lGzdem = 0;
        Nz = 0;
    else
        Nz = N(zeroDemandClasses);
        lGzdem = - sum(factln(Nz)) + sum(Nz.*log(sum(Z(:,zeroDemandClasses),1))) + Nz*log(scalevec(zeroDemandClasses))';
    end
end
L = L(:,nonzeroDemandClasses);
N = N(nonzeroDemandClasses);
Z = Z(:,nonzeroDemandClasses);
scalevecz = scalevec(nonzeroDemandClasses);
% compute G for classes No with non-zero demand
if any(N<0)
    lGnnzdem = 0;
else
    [lGnnzdem,method] = compute_norm_const_ld(L, N, Z, mu, options);
end
% scale back to original demands
lG = lGnnzdem + lGzdem + N*log(scalevecz)';
G = exp(lG);
end
 
function [lG,method] = compute_norm_const_ld(L,N,Z,mu,options)
% LG = COMPUTE_NORM_CONST_LD(L,N,Z,OPTIONS)
[M,R] = size(L);
method = options.method;
switch options.method
    case {'default','exact'}
        % In 'default' mode prefer the Choudhury-Leung-Whitt generating-function
        % inversion (pfqn_clw_lld) for genuinely multi-station, low-class-count
        % models. Two independent limits, both calibrated by profiling the
        % canonical JAR, gate its use (else fall back to exact/gld):
        %  1) Speed. Cost is exactly prod_j 2*l_j*N_j contour points (l defaults
        %     [1,2,2,3,3,...]) at a steady ~1e7 points/s, so runtime ~=
        %     clwPredCost/1e7 s. It grows as the product of the per-class
        %     populations (degree R in the total population for balanced classes:
        %     ~2*Ntot^2 at R=2, ~Ntot^3 at R=3) and exponentially in R. Cap at a
        %     ~2s budget (clwMaxCost).
        %  2) Numerical validity. The inversion sums 2*N_j alternating-sign
        %     contour points per class, so it loses accuracy to catastrophic
        %     cancellation as the population grows: clw_lld returns NaN
        %     intermittently above per-class ~450 (JAR; native ~700). Restrict
        %     to total population sum(N) <= clwMaxPop = 200, well inside the
        %     reliable regime; larger models revert to the exact default.
        % The explicit 'exact' method always takes the convolution path below;
        % 'clw' forces the inversion irrespective of these caps.
        clwMaxClasses = 5;    % class-count gate (clw cost is exponential in R)
        clwMaxPop = 200;      % total-population cap (numerical validity; NaN onset ~450)
        clwMaxCost = 2e7;     % contour-point budget (~2s at ~1e7 pts/s; profiler)
        clwPredCost = prod(2 * clw_lattice(R) .* N(:).');
        if strcmp(options.method,'default') && M > 1 && R >= 2 ...
                && R <= clwMaxClasses && sum(N) <= clwMaxPop ...
                && clwPredCost <= clwMaxCost
            [~,lG] = pfqn_clw_lld(L, N, sum(Z,1), mu, options);
            method = 'clw';
        else
            if sum(Z(:))<GlobalConstants.FineTol
                Lz = L;
                muz = mu;
            else
                D = size(Z,1); % number of delays
                Lz = [L;Z];
                muz = [mu; repmat(1:size(mu,2),D,1)];
            end
            if R==1
                [lG] = pfqn_gldsingle(Lz, N, muz, options);
                method = 'exact/gld';
            elseif M==1 && any(Z>0)
                [~,lG]= pfqn_comomrm_ld(L, N, Z, mu, options);
                method = 'exact/comomld';
            elseif M==1 && max(Z(:)) < GlobalConstants.FineTol
                % M counts QUEUEING stations; pfqn_comomrm_ld is the finite
                % repairman routine (one queueing station, plus an optional
                % INF-server row it folds into the think time). This guard read
                % M==2 before; it is now M==1 so a genuine two-queueing-station
                % model (no delay) falls through to pfqn_gld below, which matches
                % pfqn_ca to machine precision on those models.
                [~,lG] = pfqn_comomrm_ld(L, N, 0*N, mu, options);
                method = 'exact/comomld';
            else
                [~,lG] = pfqn_gld(Lz, N, muz, options);
                method = 'exact/gld';
            end
        end
    case {'is'}
        % Importance sampling for a load-dependent closed network: the
        % sample-an-ordering estimator of pfqn_ld_is (the LD counterpart of
        % pfqn_is / pfqn_oi_is / pfqn_pas_is).
        [~,lG] = pfqn_ld_is(L,N,sum(Z,1),mu,options);
        method = 'is';
    case 'clw'
        % Choudhury-Leung-Whitt generating-function inversion extended to
        % limited load-dependent stations (Bertozzi-McKenna transforms); the
        % delay term is passed as the aggregate IS demand sum(Z,1).
        [~,lG] = pfqn_clw_lld(L, N, sum(Z,1), mu, options);
        method = 'clw';
    case 'rd'
        [lG] = pfqn_rd(L, N, Z, mu, options);
    case 'nrp'
        [lG] = pfqn_nrp(L, N, Z, mu, options);
    case 'nrl'
        [lG] = pfqn_nrl(L, N, Z, mu, options);
    case 'comomld'
        if M<=1 || sum(Z) < GlobalConstants.Zero
            [~,lG]= pfqn_comomrm_ld(L, N, Z, mu, options);
        else
            line_warning(mfilename,'Load-dependent CoMoM is available only in models with a delay and m identical stations, running the ''rd'' algorithm instead.\n');
            [lG] = pfqn_rd(L, N, Z, mu, options);
            method = 'rd';
        end
    otherwise
        line_error(mfilename,sprintf('Unrecognized method for solving load-dependent models: %s',options.method));
end
return
end
 
function l = clw_lattice(p)
% Default CLW inner lattice parameters l_j (roundoff control): l_1=1,
% l_2=l_3=2, l_j>=4 = 3. Used to predict the clw_lld contour-point cost.
l = 3 * ones(1, p);
l(1) = 1;
if p >= 2, l(2) = 2; end
if p >= 3, l(3) = 2; end
end