pfqn_oi_is.m

function [G, lG, Q] = pfqn_oi_is(N, mu, options)
% [G, LG, Q] = PFQN_OI_IS(N, MU, OPTIONS)
%
% Importance-sampling (IS) estimate of the normalizing constant of a closed
% two-station order-independent (OI) tandem. This is PFQN_PAS_IS with an EMPTY
% swap graph: it does not restrict the sampled orderings to a pass-and-swap
% communicating class but samples over ALL microstates (every ordering of the
% present jobs is feasible), thereby recovering the plain OI normalizing
% constant of PFQN_OI_NC (which it estimates by Monte Carlo rather than exact
% balanced-fairness convolution).
%
% Model. Two OI stations (1 = upstream, 2 = downstream) hold all N jobs (no
% delay). The ordered-state chain is irreducible, so the communicating class is
% the full set of orderings D = all permutations of the job multiset. Writing
% Phi_m for the balanced-fairness balance function of station m,
%   G = sum_{c in D} sum_{k=0}^{ell} Phi_1(c_{1..k}) Phi_2(c_{ell..k+1}),
% with the ordered product Phi_m(q) = prod_{p=1}^{|q|} 1/mu_m(supp(q_{1..p})).
%
% Auto-normalized IS. Orderings c are drawn by placing, at each step, a
% uniformly random present class (no placement constraint); the draw probability
% p(c) is the product of the reciprocal branching factors. Then
%   G[xi] = E_{C~p}[ (sum_k xi(C,k) Phi_1(C_{1..k}) Phi_2(C_{ell..k+1})) / p(C) ],
% and E[xi] = G[xi]/G[1] reuses the same samples (auto-normalized IS). Taking
% xi = number of class-r jobs in the prefix yields the class-r mean queue length
% at station 1.
%
% Parameters:
%   N       - (1 x R) closed population vector (the macrostate), finite.
%   mu      - cell {1 x 2} of function handles. mu{m}(n) returns the total OI
%             rank rate of station m for the per-class occupancy (count) vector
%             n (1 x R); OI, so it depends only on supp(n). This is the
%             svcRateFun stored on an OI station.
%   options - solver options (optional). Fields used:
%               .samples  number of IS samples (default 1e4);
%               .seed     RNG seed for reproducibility (optional);
%               .verbose  print progress (default false).
%
% Returns:
%   G  - IS estimate of the OI normalizing constant (== PFQN_OI_NC with Z=0).
%   lG - log(G).
%   Q  - (2 x R) IS estimate of the mean per-class queue length; Q(1,:) at
%        station 1, Q(2,:) = N - Q(1,:) at station 2.
%
% Example (two OI queues, R=5, no swap graph):
%   mu1 = [1 2 0.5]; nb1 = {1,2,3,[1 3],[2 3]};
%   mu2 = [1 2];     nb2 = {1,2,[1 2],1,2};
%   rate = @(nb,muv,n) sum(muv(unique([nb{n>0}])));
%   mu = {@(n) rate(nb1,mu1,n), @(n) rate(nb2,mu2,n)};
%   [G,lG,Q] = pfqn_oi_is([1 1 1 3 3], mu, struct('samples',1e5));
%
% See also PFQN_PAS_IS, PFQN_OI_NC, PFQN_OI_FNC.
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
if nargin < 3, options = struct(); end
if ~iscell(mu), mu = {mu}; end
if numel(mu) ~= 2
    line_error(mfilename, 'pfqn_oi_is models a two-station OI tandem: mu must have exactly two rate functions.');
end
 
R = numel(N);
N = round(N(:)');
if any(~isfinite(N))
    line_error(mfilename, 'pfqn_oi_is requires finite (closed) populations.');
end
 
nsamples = 1e4;
if isfield(options, 'samples') && ~isempty(options.samples)
    nsamples = round(options.samples);
end
if isfield(options, 'seed') && ~isempty(options.seed)
    rng(options.seed);
end
verbose = isfield(options, 'verbose') && ~isempty(options.verbose) && options.verbose;
 
ell = sum(N);
if ell == 0
    G = 1; lG = 0; Q = zeros(2, R);
    return
end
 
nCoef = R + 1;                 % xi = [1, n_{1,1}, ..., n_{1,R}]
accum = zeros(1, nCoef);
 
for s = 1:nsamples
    % ---- draw an ordering c by placing a uniformly random present class -----
    x = N;                     % remaining per-class jobs to place
    c = zeros(1, ell);
    logp = 0;
    for ppos = 1:ell
        avail = find(x > 0);   % OI: every present class is placeable
        na = numel(avail);
        pick = avail(randi(na));
        c(ppos) = pick;
        logp = logp - log(na);
        x(pick) = x(pick) - 1;
    end
    p_c = exp(logp);
 
    % ---- prefix balance Phi_1 and per-class counts at each cut -------------
    Phi1 = ones(1, ell + 1);
    cnt1 = zeros(ell + 1, R);
    supp = zeros(1, R); phi = 1; occ = zeros(1, R);
    for k = 1:ell
        cls = c(k);
        supp(cls) = 1; occ(cls) = occ(cls) + 1;
        phi = phi / mu{1}(supp);
        Phi1(k + 1) = phi;
        cnt1(k + 1, :) = occ;
    end
    % ---- reversed-suffix balance Phi_2 (station 2 = reversed suffix) -------
    Phi2cut = ones(1, ell + 1);
    supp = zeros(1, R); phi = 1;
    for k = ell:-1:1
        cls = c(k);
        supp(cls) = 1;
        phi = phi / mu{2}(supp);
        Phi2cut(k) = phi;
    end
 
    % ---- split-convolution sample value for every coefficient -------------
    sv = zeros(1, nCoef);
    for k = 0:ell
        if k >= ell
            w2 = 1;
        else
            w2 = Phi2cut(k + 1);
        end
        w = Phi1(k + 1) * w2;
        sv(1) = sv(1) + w;
        if k > 0
            sv(2:end) = sv(2:end) + w * cnt1(k + 1, :);
        end
    end
    accum = accum + sv / p_c;
 
    if verbose && mod(s, max(1, floor(nsamples / 10))) == 0
        line_printf('\npfqn_oi_is: %d/%d samples', s, nsamples);
    end
end
 
est = accum / nsamples;
G = est(1);
lG = log(G);
Q = zeros(2, R);
if G > 0
    Q(1, :) = est(2:end) / G;
end
Q(2, :) = N - Q(1, :);
end