pfqn_oi_nc.m

function [G, lG] = pfqn_oi_nc(Z, N, mu, options)
% [G, LG] = PFQN_OI_NC(Z, N, MU, OPTIONS)
%
% Normalizing constant for a closed product-form queueing network that
% comprises a single aggregated infinite-server (delay) node and an
% arbitrary number of order-independent (OI) / pass-and-swap stations with
% empty swap graph.
%
% The OI stations are analyzed by the balanced-fairness recursion of
% Bonald & Proutiere (2003), "Insensitive bandwidth sharing in data
% networks", combined with the multichain convolution over stations. For a
% single OI station with rank rate mu(supp(n)) the balance function is
%   Phi(0) = 1,  Phi(n) = (1/mu(supp n)) * sum_{r: n_r>0} Phi(n - e_r),
% and G(N) is obtained by convolving the per-station balance functions with
% the multinomial delay factor F_Z(n) = prod_r Z_r^{n_r} / n_r!.
%
% The routine returns the exact G(N) via the equivalent recursive
% peeling: peel the last OI station (empty) then, per class r with N_r>0,
% place one class-r job at that station and recurse on N - e_r with the
% station rate function shifted by e_r. This reproduces the balanced-
% fairness convolution without materializing the full Phi tables.
%
% Parameters:
%   Z  - (1 x R) think-time demand vector of the aggregated delay node.
%        Z(r) = 1/sigma_r for a delay with per-class rate sigma_r.
%   N  - (1 x R) closed population vector, finite.
%   mu - cell array {1 x M} of function handles, one per OI station. Each
%        mu{m}(n) returns the total service rate of station m given the
%        per-class occupancy (count) vector n (1 x R). For an order-
%        independent station this rate depends only on the support of n
%        (which classes are present), i.e. mu{m}(n) = sum of the capacities
%        of the servers compatible with the classes present in n. May be
%        empty to model a pure delay network.
%   options - solver options (optional, currently unused; accepted for
%             signature parity with the other pfqn_* routines).
%
% Returns:
%   G  - Normalizing constant G(N).
%   lG - log(G(N)).
%
% Example (IS + two OI stations, R classes):
%   oirate  = @(n) sum(mu1(any(compat1(:, find(n>0)) ~= 0, 2)));
%   oirate2 = @(n) sum(mu2(any(compat2(:, find(n>0)) ~= 0, 2)));
%   G = pfqn_oi_nc(1./sigma, N, {oirate, oirate2});
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
if nargin < 4
    options = struct();
end
if nargin < 3 || isempty(mu)
    mu = {};
end
if ~iscell(mu)
    mu = {mu};
end
 
R = numel(N);
if isempty(Z)
    Z = zeros(1, R);
end
if numel(Z) ~= R
    line_error(mfilename, 'Z and N must have the same number of classes.');
end
if any(~isfinite(N))
    line_error(mfilename, 'pfqn_oi_nc requires finite (closed) populations.');
end
 
N = round(N(:)');
Z = Z(:)';
 
G = oi_nc_rec(Z, N, mu);
lG = log(G);
end
 
function G = oi_nc_rec(Z, N, mu)
% Recursive balanced-fairness convolution over the OI stations.
R = numel(N);
 
if sum(N) == 0
    G = 1;
    return
end
 
M = numel(mu);
if M == 0
    % Base case: only the aggregated delay node remains. Its unnormalized
    % weight is the multinomial factor prod_r Z_r^{N_r} / N_r!.
    logf = 0;
    for r = 1:R
        if N(r) > 0
            if Z(r) <= 0
                G = 0; % class r has population but no delay demand: infeasible weight
                return
            end
            logf = logf + N(r) * log(Z(r)) - gammaln(N(r) + 1);
        end
    end
    G = exp(logf);
    return
end
 
% Step A: last OI station holds no jobs; peel it off.
mu_sub = mu(1:M-1);
G = oi_nc_rec(Z, N, mu_sub);
 
% Step B: place one job of each present class at the active OI station and
% recurse with the station rate function shifted to account for that job.
active = mu{M};
for r = 1:R
    if N(r) > 0
        e_r = zeros(1, R);
        e_r(r) = 1;
        mu_r = active(e_r);
        if mu_r <= 0
            continue
        end
        shifted = @(state) active(state + e_r);
        mu_next = mu;
        mu_next{M} = shifted;
        Nm = N;
        Nm(r) = Nm(r) - 1;
        G = G + (1 / mu_r) * oi_nc_rec(Z, Nm, mu_next);
    end
end
end