pfqn_ca.m

%{
%{
 % @file pfqn_ca.m
 % @brief Convolution Algorithm for exact normalizing constant computation.
%}
%}
 
function [Gn,lGn]=pfqn_ca(L,N,Z)
%{
%{
 % @brief Convolution Algorithm for exact normalizing constant computation.
 % @fn pfqn_ca(L, N, Z)
 % @param L Service demand matrix.
 % @param N Population vector.
 % @param Z Think time vector.
 % @return Gn Normalizing constant.
 % @return lGn Logarithm of the normalizing constant.
%}
%}
[M,R]=size(L);
if nargin<3 || isempty(Z)
    Z=zeros(1,R);
end
if M==0
    lGn = - sum(factln(N)) + sum(N.*log(sum(Z,1)));
    Gn = exp(lGn);
    return
end
 
if min(N)<0
    Gn=0;
    lGn=-Inf;
    return;
end
 
if sum(N)==0
    Gn=1;
    lGn=0;
    return;
end
 
% Demand scaling, so that lGn stays computable once G(N) leaves the
% double-precision range. The recursion below runs in linear space, so it
% overflows to Inf as soon as G(N) > realmax, i.e. log G > 709.78 -- and lGn
% was then returned as Inf even though log G is perfectly representable. This
% is the floating-point range problem Reiser and Lavenberg (1980, JACM 27(2),
% p.319) report for the convolution algorithm, and the scaling remedy is Lam
% (1982), "Dynamic scaling and growth behavior of queuing network normalization
% constants".
%
% Every state in G(N) carries the same total population sum(N), so dividing all
% demands and think times by a constant c divides G(N) by exactly c^sum(N):
%   G(N; L/c, Z/c) = G(N; L, Z) / c^sum(N)
% hence log G = log G_scaled + sum(N) log c, which is exact, not an
% approximation.
%
% c must be chosen to CENTRE log G_scaled near 0, not merely to shrink the
% demands. G(N) can leave the double range in EITHER direction: Reiser (1981,
% Perf. Eval. 1:7-18, Sec. 6.1) reports the unnormalized convolution
% UNDERFLOWING (< 1e-75, losing all significant digits) on his Fig. 5
% central-server model for K > 160, and OVERFLOWING (> 1e75) on the same model
% if the think time is raised. Scaling by max(L), the textbook choice, only
% shrinks G and so makes the underflow strictly worse -- on Reiser's model it
% turns a representable G(200) ~ 1e-87 into an underflow to 0.
%
% log G is estimated first, from the largest SINGLE-STATE term, which is a
% lower bound on G(N) and in practice within O(log #states) of log G:
%   all jobs at one queueing station i : sum_r N_r log L_ir
%   all jobs at the delay              : sum_r [N_r log Z_r - log N_r!]
% On Reiser's model at K=200 this gives max(-460.5, -264.1) = -264.1 against a
% true log G of about -263, i.e. accurate enough to place the scaled value well
% inside the exponent range.
%
% c is then the power of two nearest exp(lGest/sum(N)), so that
% log G_scaled = log G - sum(N) log c is near 0. A power of two matters: it
% shifts exponents only, so L/c and Z/c stay exactly representable and the
% recursion is bit-for-bit the unscaled one with shifted exponents.
Nt = sum(N);
lGest = -Inf;
for i = 1:M
    t = 0; ok = true;
    for r = 1:R
        if N(r) > 0
            if L(i,r) > 0
                t = t + N(r)*log(L(i,r));
            else
                ok = false; break
            end
        end
    end
    if ok, lGest = max(lGest, t); end
end
if any(Z > 0)                       % all jobs at the delay
    t = 0; ok = true;
    for r = 1:R
        if N(r) > 0
            if sum(Z(:,r)) > 0
                t = t + N(r)*log(sum(Z(:,r))) - factln(N(r));
            else
                ok = false; break
            end
        end
    end
    if ok, lGest = max(lGest, t); end
end
if ~isfinite(lGest)
    kscale = 0;
else
    kscale = round(lGest/(Nt*log(2)));
end
cscale = pow2(kscale);
L = L / cscale;
Z = Z / cscale;
 
G = ones(M+1,prod(N+1)); % stores G across recursion
n = pprod(N);
while sum(n)~=-1
    idxn = hashpop(n,N);
    G(1,idxn) = Fz(Z,n);
    for m=2:M+1
        G(m,idxn) = G(m-1,idxn); % norm constant with m-1 queues
        for r=1:R
            if n(r)>=1
                n(r) = n(r)-1;
                idxn_1r = hashpop(n,N);
                n(r) = n(r)+1;
                G(m,idxn) = G(m,idxn) + L(m-1,r)*G(m,idxn_1r);
            end
        end
    end
    n=pprod(n,N);
end
% Undo the scaling in log space: log G = log G_scaled + sum(N) log c. lGn is
% therefore finite whenever log G itself is, even though Gn may legitimately
% overflow to Inf (the true constant really is outside double range).
lGn = log(G(M+1,end)) + sum(N)*kscale*log(2);
% Gn is recovered by scaling the mantissa back with pow2 (an exact exponent
% adjustment), not as exp(lGn): on models that never overflowed this returns
% bit-for-bit the value the unscaled recursion used to return, and it still
% goes to Inf when the constant genuinely leaves double range -- in which case
% lGn above remains finite and usable.
Gn = pow2(G(M+1,end), sum(N)*kscale);
end
 
function idx=hashpop(n,N,R,prods)
% IDX=HASHPOP(N,N,R,PRODS)
 
% hash a population vector in n: 0<=n<=N
idx=1;
if nargin==2
    R=length(N);
    for r=1:R
        idx= idx + prod(N(1:r-1)+1)*n(r);
    end
    return
else
    for r=1:R
        idx= idx + prods(r)*n(r);
    end
end
end
 
function [n]=pprod(n,N)
% [N]=PPROD(N,N)
 
% sequentially generate all vectors n: 0<=n<=N
% n=pprod(N) - init
% n=pprod(n,N) - next state
if nargin==1
    N=n;
    n=zeros(size(N));
    return;
end
 
R=length(N);
if sum(n==N)==R
    n=-1;
    return
end
 
s=R;
while s>0 && n(s)==N(s)
    n(s)=0;
    s=s-1;
end
if s==0
    %n=-1*ones(1,R);
    return
end
n(s)=n(s)+1;
return;
end
 
function f=Fz(Z,n)
% F=FZ(Z,N)
 
R=length(n);
if sum(n)==0
    f=1;
    return
end
f=0;
for r=1:R
    if Z(r)>0
        f=f+log(Z(r))*n(r);
        f=f-gammaln(1+n(r));
    elseif n(r)>0
        f = 0;
        return
    end
end
f=exp(f);
end