"""
Order-independent (OI) normalizing-constant tools.
Native Python port of the MATLAB OI functional-server routines:
pfqn_oi_nc - balanced-fairness normalizing constant for a closed network
of order-independent stations and one aggregated delay node.
pfqn_oi_fnc - OI generalization of the load-dependent functional server
(FNC) of Casale, "On Single-Class Load-Dependent Normalizing
Constant Equations", QEST 2006 (Theorem 3, Corollary 1).
References:
Original MATLAB: matlab/src/api/pfqn/pfqn_oi_nc.m, pfqn_oi_fnc.m
Bonald and Proutiere, "Insensitive bandwidth sharing in data networks"
(2003); Casale, QEST 2006.
"""
import numpy as np
from typing import Callable, List, Optional, Sequence, Tuple, Union
[docs]
def pfqn_oi_nc(Z: Sequence[float],
N: Sequence[int],
mu: Optional[Union[Callable, List[Callable]]] = None,
options=None) -> Tuple[float, float]:
"""Normalizing constant of a closed OI + single-delay product-form network.
The OI stations are analyzed by the balanced-fairness recursion of Bonald
and Proutiere (2003) combined with the multichain convolution over
stations. The exact G(N) is obtained by 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.
Parameters
----------
Z : (R,) think-time demand vector of the aggregated delay node.
N : (R,) closed population vector, finite.
mu : list of callables, one per OI station. Each ``mu[m](n)`` returns the
total service rate of station m given the per-class occupancy (count)
vector n. May be empty/None to model a pure delay network.
options : accepted for signature parity; unused.
Returns
-------
(G, lG) : normalizing constant and its natural log.
"""
N = np.round(np.asarray(N, dtype=float)).astype(int).ravel()
R = N.size
if mu is None:
mu = []
elif callable(mu):
mu = [mu]
else:
mu = list(mu)
if Z is None or len(np.asarray(Z).ravel()) == 0:
Z = np.zeros(R)
Z = np.asarray(Z, dtype=float).ravel()
if Z.size != R:
raise ValueError('pfqn_oi_nc: Z and N must have the same number of classes.')
G = _oi_nc_rec(Z, N, mu)
lG = np.log(G) if G > 0 else -np.inf
return G, lG
def _oi_nc_rec(Z: np.ndarray, N: np.ndarray, mu: List[Callable]) -> float:
"""Recursive balanced-fairness convolution over the OI stations."""
R = N.size
if int(N.sum()) == 0:
return 1.0
M = len(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.0
for r in range(R):
if N[r] > 0:
if Z[r] <= 0:
return 0.0 # population but no delay demand: infeasible
from scipy.special import gammaln
logf += N[r] * np.log(Z[r]) - gammaln(N[r] + 1)
return float(np.exp(logf))
# Step A: last OI station holds no jobs; peel it off.
mu_sub = mu[: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 - 1]
for r in range(R):
if N[r] > 0:
e_r = np.zeros(R, dtype=int)
e_r[r] = 1
mu_r = active(e_r)
if mu_r <= 0:
continue
shifted = _make_shifted(active, e_r)
mu_next = list(mu)
mu_next[M - 1] = shifted
Nm = N.copy()
Nm[r] -= 1
G += (1.0 / mu_r) * _oi_nc_rec(Z, Nm, mu_next)
return G
def _make_shifted(active: Callable, e_r: np.ndarray) -> Callable:
"""Return the rate handle ``state -> active(state + e_r)``."""
return lambda state: active(np.asarray(state) + e_r)
[docs]
def pfqn_mvaoi(Z: Sequence[float],
N: Sequence[int],
mu: Union[Callable, List[Callable]],
Dli: Optional[Sequence[Sequence[float]]] = None,
options=None) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""Mean-value analysis of a closed product-form OI network.
Mean-value counterpart of :func:`pfqn_oi_nc` and the marginal form
:func:`pfqn_mvaoi_marg`: for a closed product-form network of an aggregated
infinite-server (delay) node, any number of load-independent (LI) single-server
product-form queues, and any number of order-independent (OI) stations, it
returns the same exact per-class throughput and queue-lengths WITHOUT computing
any normalizing constant or joint marginal, using only mean quantities. It is
the composition-dependent generalization of the Conditional MVA (CMVA) of
Casale, "A Note on Stable Flow-Equivalent Aggregation in Closed Networks"
(QUESTA 2009), extended to MULTIPLE OI stations by carrying one rate-shift
vector ``s_i`` per OI station i (row ``i`` of the shift matrix ``S``).
Throughout, ``r`` and ``s`` index job classes; ``i`` indexes OI stations; ``j``
indexes LI queues. State ``(S, Nn)`` is processed by increasing ``sum(Nn)``;
each OI station keeps its own ``D^i``, ``rho^i`` and ``Q^i`` recursions driven
by the common throughput ``X^{(S)}(Nn)``, and the population conservation
aggregates every station's contribution::
Nn_r = X_r Z_r + sum_j Q^{(j)}_r + sum_i Q^{(i)}_r,
with the LI queue term ``Q^{(j)}_r = X_r D_{j,r}(1 + sum_s Q^{(j)}_s(Nn - e_r))``.
Parameters
----------
Z : (R,) think-time demand vector of the aggregated delay node.
N : (R,) closed population vector, finite.
mu : callable or list of callables ``mu_i(n)`` returning the OI total service
rate of station i for the per-class occupancy (count) vector n. A bare
callable is accepted as the single-station shorthand.
Dli : (J, R) per-class demand matrix of the LI single-server queues; None or
empty when J = 0.
options : accepted for signature parity; unused.
Returns
-------
(X, Qoi, Qli, Qdelay, Soi) : per-class throughput (R,), OI queue-lengths
(K, R), LI queue-lengths (J, R), delay queue-length (R,) = ``X * Z``, and
the per-class mean number of IN-SERVICE jobs at each OI station (K, R).
``Soi[i, r] = E[sir_r]`` counts the class-r jobs receiving a strictly
positive rank rate (see :func:`pfqn_oi_insvc`); the utilization of OI
station i is ``Soi[i, r] / c_i``. Unlike X/Qoi/Qli, which are pure
mean-value quantities, Soi is a distributional statistic and is obtained
from the OI count marginal assembled from the zero-shift throughputs
``X^{(0)}(k)`` already cached by the mean-value recursion above (no
normalizing constant is formed).
"""
if callable(mu):
mu = [mu]
mu = list(mu)
if len(mu) == 0 or not all(callable(m) for m in mu):
raise ValueError('mu must be a (nonempty) list of OI rate callables.')
K = len(mu)
Z = np.asarray(Z, dtype=float).ravel()
N = np.round(np.asarray(N, dtype=float)).astype(int).ravel()
R = N.size
if np.any(~np.isfinite(N)):
raise ValueError('pfqn_mvaoi requires finite (closed) populations.')
if Dli is None:
Dli = np.zeros((0, R))
Dli = np.asarray(Dli, dtype=float).reshape(-1, R) if np.size(Dli) else np.zeros((0, R))
J = Dli.shape[0]
ei = np.eye(R, dtype=int)
Xc: dict = {} # X^{(S)}(Nn)
Qlc: dict = {} # Qli^{(S)}(Nn) (J, R)
Dc = [dict() for _ in range(K)] # D_i^{(S)}(Nn)
Qc = [dict() for _ in range(K)] # Q_i^{(S)}(Nn)
rho_cache = [dict() for _ in range(K)]
def skey(S, Nn):
return (tuple(np.asarray(S, dtype=int).ravel()), tuple(np.asarray(Nn, dtype=int)))
def shift_plus(S, i, r):
Sp = np.array(S, dtype=int, copy=True)
Sp[i] = Sp[i] + ei[r]
return Sp
def rho(i, r, S, M):
M = np.asarray(M, dtype=int)
key = (skey(S, M), r)
v = rho_cache[i].get(key)
if v is not None:
return v
if M.sum() == 0:
rho_cache[i][key] = 1.0
return 1.0
s = next(ss for ss in range(R) if ss != r and M[ss] > 0)
xu = Xc[skey(S, M)][s]
xu2 = Xc[skey(shift_plus(S, i, r), M)][s]
val = rho(i, r, S, M - ei[s]) * (xu / xu2 if xu2 > 0 else 0.0)
rho_cache[i][key] = val
return val
zeroS = np.zeros((K, R), dtype=int)
k0 = skey(zeroS, np.zeros(R, dtype=int))
Xc[k0] = np.zeros(R)
Qlc[k0] = np.zeros((J, R))
for i in range(K):
Dc[i][k0] = np.zeros(R)
Qc[i][k0] = np.zeros(R)
# Enumerate all reachable (S, Nn) with, per class r, the K+1 buckets (K OI
# shifts, then the free bucket Nn) summing to <= N_r. Each per-class row is the
# first K+1 entries of a (K+2)-part composition of N_r whose dropped (K+2)-th
# part is the slack bucket; this mirrors MATLAB sprod(K+2, N) with the slack
# row removed. States are then processed by increasing sum(Nn).
from itertools import product
perclass = []
for r in range(R):
rows = [comp[:K + 1] for comp in _compositions(int(N[r]), K + 2)]
perclass.append(rows)
states = []
for combo in product(*perclass):
S = np.zeros((K, R), dtype=int)
Nn = np.zeros(R, dtype=int)
for r in range(R):
row = combo[r]
for i in range(K):
S[i, r] = row[i]
Nn[r] = row[K]
states.append((S, Nn))
states.sort(key=lambda st: int(st[1].sum()))
for S, Nn in states:
kk = skey(S, Nn)
if Nn.sum() == 0:
Xc[kk] = np.zeros(R)
Qlc[kk] = np.zeros((J, R))
for i in range(K):
Dc[i][kk] = np.zeros(R)
Qc[i][kk] = np.zeros(R)
continue
Dt = np.zeros((K, R))
Qsub = [np.zeros((R, R)) for _ in range(K)]
for i in range(K):
for r in range(R):
if Nn[r] == 0:
continue
if Nn[r] == 1:
mur = mu[i](S[i] + ei[r])
if mur > 0:
Dt[i, r] = (1.0 / mur) * rho(i, r, S, Nn - ei[r])
else:
Nmr = Nn - ei[r]
xs = Xc[skey(S, Nmr)][r]
xs2 = Xc[skey(shift_plus(S, i, r), Nmr)][r]
if xs2 > 0:
Dt[i, r] = (xs / xs2) * Dc[i][skey(S, Nmr)][r]
for s in range(R):
if Nn[s] > 0:
Qsub[i][s] = Qc[i][skey(shift_plus(S, i, s), Nn - ei[s])]
betaLI = np.zeros((J, R))
for r in range(R):
if Nn[r] == 0:
continue
QsubLI = Qlc[skey(S, Nn - ei[r])]
for j in range(J):
betaLI[j, r] = Dli[j, r] * (1.0 + float(np.sum(QsubLI[j, :])))
idx = [r for r in range(R) if Nn[r] > 0]
A = np.zeros((len(idx), len(idx)))
for a, r in enumerate(idx):
for b, s in enumerate(idx):
if s == r:
val = Z[r] + float(np.sum(betaLI[:, r]))
for i in range(K):
val += Dt[i, r] * (1.0 + Qsub[i][r, r])
A[a, b] = val
else:
val = 0.0
for i in range(K):
val += Dt[i, s] * Qsub[i][s, r]
A[a, b] = val
# Throughput closure via the product-form ratio identity (valid at fixed
# shift S): X_s(S,Nn)/X_r(S,Nn) = X_s(S,Nn-e_r)/X_r(S,Nn-e_s), giving the
# scalar per-class recurrence and avoiding the linear solve.
Xk = np.zeros(R)
for a, r in enumerate(idx):
denom = A[a, a]
for b, s in enumerate(idx):
if b != a:
Xner = Xc[skey(S, Nn - ei[r])] # X(S, Nn-e_r)
Xnes = Xc[skey(S, Nn - ei[s])] # X(S, Nn-e_s)
if Xnes[r] > 0:
denom += A[a, b] * (Xner[s] / Xnes[r])
if denom > 0:
Xk[r] = Nn[r] / denom
QkLI = np.zeros((J, R))
for r in range(R):
if Nn[r] == 0:
continue
QkLI[:, r] = Xk[r] * betaLI[:, r]
for i in range(K):
U = Dt[i] * Xk
Qi = np.array([U[r] + float(np.dot(U, Qsub[i][:, r])) for r in range(R)])
Qc[i][kk] = Qi
Dc[i][kk] = Dt[i].copy()
Xc[kk] = Xk
Qlc[kk] = QkLI
keyN = skey(zeroS, N)
X = Xc[keyN]
Qoi = np.array([Qc[i][keyN] for i in range(K)]) if K > 0 else np.zeros((0, R))
Qli = Qlc[keyN]
Soi = _oi_insvc_means(N, mu, Xc, skey, zeroS, K, R)
return X, Qoi, Qli, X * Z, Soi
def _oi_insvc_means(N, mu, Xc, skey, zeroS, K, R):
"""Mean number of in-service jobs per class at each OI station, E[sir_r].
Built from the OI count marginal pM_i(n|k) on the cached zero-shift
throughputs X^{(0)}(k). Exact because in product form
``pM_i(n|k) = Phi_i(n) G_{-i}(k-n)/G(k)`` and ``X_r(k) = G(k-e_r)/G(k)``, so
the recursion below reproduces the balanced-fairness identity for Phi_i::
pM_i(n|k) = (1/mu_i(n)) sum_r X_r(k) pM_i(n-e_r|k-e_r)
pM_i(0|k) = 1 - sum_{n != 0} pM_i(n|k)
"""
N = np.asarray(N, dtype=int).ravel()
shp = N + 1
stride = np.ones(R, dtype=int)
for d in range(1, R):
stride[d] = stride[d - 1] * shp[d - 1]
total = int(np.prod(shp))
subs = np.zeros((total, R), dtype=int)
for i in range(total):
li = i
for d in range(R):
subs[i, d] = li % shp[d]
li //= shp[d]
order = sorted(range(total), key=lambda i: subs[i].sum())
Xk = np.zeros((total, R))
for i in range(total):
Xk[i, :] = Xc[skey(zeroS, subs[i])]
Soi = np.zeros((K, R))
for m in range(K):
gm, _, _ = pfqn_oi_insvc(mu[m], N)
muv = np.zeros(total)
for i in range(total):
if subs[i].sum() > 0:
muv[i] = float(mu[m](subs[i].copy()))
pMv = np.zeros((total, total))
pMv[0, 0] = 1.0 # pM(0|0) = 1
for b in order:
k = subs[b]
if k.sum() == 0:
continue
acc0 = 0.0
for a in order:
n = subs[a]
if n.sum() == 0 or np.any(n > k) or muv[a] <= 0:
continue
acc = 0.0
for r in range(R):
if n[r] > 0:
acc += Xk[b, r] * pMv[a - stride[r], b - stride[r]]
pMv[a, b] = acc / muv[a]
acc0 += pMv[a, b]
pMv[0, b] = 1.0 - acc0 # empty state by complement
idxN = int(np.sum(N * stride))
for r in range(R):
Soi[m, r] = float(pMv[:, idxN] @ gm[:, r])
return Soi
def _compositions(m: int, p: int) -> List[Tuple[int, ...]]:
"""All nonnegative integer tuples of length p summing to m."""
if p == 1:
return [(m,)]
out = []
for first in range(m + 1):
for rest in _compositions(m - first, p - 1):
out.append((first,) + rest)
return out
[docs]
def pfqn_mvaoi_marg(D: Sequence[Sequence[float]],
N: Sequence[int],
isDelay: Sequence[bool],
mu: List[Optional[Callable]]
) -> Tuple[np.ndarray, np.ndarray]:
"""Exact marginal load-dependent MVA for OI networks.
Marginal-distribution counterpart of :func:`pfqn_mvaoi`. Carries, for each OI
station, its joint count-vector marginal ``pM_i(n | k)`` and closes the per-class
throughput by population conservation. Handles delay + LI product-form queues +
any number of OI stations.
Parameters
----------
D : (M, R) per-class demand at every station (OI rows ignored).
N : (R,) closed population vector, finite.
isDelay : (M,) True for infinite-server (delay) stations.
mu : length-M list; ``mu[i]`` is the OI rate callable of the count vector n, or
None for non-OI stations.
Returns
-------
(XN, QN) : per-class throughput (R,) and per-station queue-lengths (M, R).
"""
D = np.asarray(D, dtype=float)
M = D.shape[0]
N = np.round(np.asarray(N, dtype=float)).astype(int).ravel()
R = N.size
isDelay = np.asarray(isDelay, dtype=bool).ravel()
oi_list = [i for i in range(M) if mu[i] is not None]
nOI = len(oi_list)
if nOI == 0:
raise ValueError('pfqn_mvaoi_marg requires at least one OI station.')
ei = np.eye(R, dtype=int)
from itertools import product
def vecs_upto(k):
return [np.array(t, dtype=int)
for t in product(*[range(int(k[r]) + 1) for r in range(R)])]
zero = tuple(np.zeros(R, dtype=int))
Xc = {zero: np.zeros(R)}
Qc = {zero: np.zeros((M, R))}
pM = [{zero: {zero: 1.0}} for _ in range(nOI)]
def _muM(fun, n):
cls = np.repeat(np.arange(R), np.asarray(n, dtype=int))
if cls.size == 0:
return 0.0
return float(fun(cls))
def pM_given(o, k, Xk):
fun = mu[oi_list[o]]
pmk = {}
for n in vecs_upto(k):
if n.sum() == 0:
continue
rate = _muM(fun, n)
if rate <= 0:
continue
acc = 0.0
for r in range(R):
if n[r] >= 1 and k[r] >= 1:
acc += Xk[r] * pM[o][tuple(k - ei[r])].get(tuple(n - ei[r]), 0.0)
pmk[tuple(n)] = acc / rate
pmk[zero] = 1.0 - sum(pmk.values())
return pmk
pops = sorted(vecs_upto(N), key=lambda v: (v.sum(), tuple(v)))
for k in pops:
kt = tuple(k)
if k.sum() == 0:
continue
Rfix = np.zeros((M, R))
A = np.zeros(R)
for r in range(R):
if k[r] == 0:
continue
Qkr = Qc[tuple(k - ei[r])]
for i in range(M):
if mu[i] is not None:
continue
if isDelay[i]:
Rfix[i, r] = D[i, r]
else:
Rfix[i, r] = D[i, r] * (1.0 + float(np.sum(Qkr[i, :])))
A[r] += Rfix[i, r]
Xk = np.array([k[r] / (A[r] + 1.0) if k[r] > 0 else 0.0 for r in range(R)])
pmk = [None] * nOI
for _ in range(2000):
QMtot = np.zeros(R)
for o in range(nOI):
pmk[o] = pM_given(o, k, Xk)
for n_t, p in pmk[o].items():
QMtot += np.array(n_t) * p
Xnew = np.array([max((k[r] - QMtot[r]) / A[r], 0.0) if (k[r] > 0 and A[r] > 0)
else 0.0 for r in range(R)])
if np.max(np.abs(Xnew - Xk)) < 1e-13:
Xk = Xnew
break
Xk = 0.5 * Xk + 0.5 * Xnew
Qk = np.zeros((M, R))
for o in range(nOI):
pmk[o] = pM_given(o, k, Xk)
QM = np.zeros(R)
for n_t, p in pmk[o].items():
QM += np.array(n_t) * p
Qk[oi_list[o], :] = QM
for r in range(R):
for i in range(M):
if mu[i] is not None:
continue
Qk[i, r] = Xk[r] * Rfix[i, r]
Xc[kt] = Xk
Qc[kt] = Qk
for o in range(nOI):
pM[o][kt] = pmk[o]
Nt = tuple(N)
return Xc[Nt].copy(), Qc[Nt].copy()
[docs]
def pfqn_oi_insvc(oirate: Callable,
N: Sequence[int],
options=None) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Conditional mean number of in-service jobs per class at an OI station.
This is the quantity underlying the LINE utilization convention at
order-independent stations, ``U_r = E[sir_r] / c`` with ``c`` the number of
servers and ``sir_r`` the number of class-r jobs receiving a strictly
positive service rate.
In an OI station the state is the ordered list ``c = (c_1,...,c_n)`` of job
classes (position 1 = head) and the job in position p is served at the rank
rate increment ``Delta_p(c) = mu(c_1..c_p) - mu(c_1..c_{p-1})``, so the total
rate telescopes to ``mu(c)``. Position p is in service when
``Delta_p(c) > 0``, and ``sir_r(c) = #{p : c_p = r, Delta_p(c) > 0}``. Note
that ``sir_r`` counts JOBS, not servers: a single job served concurrently by
several compatible servers counts once. This matches the definition used by
the exact CTMC solver (``State.to_marginal``, PAS branch) and by LDES.
Because ``mu`` is permutation-invariant, the unnormalized weight of an
ordering c of the multiset n factorizes over its prefixes as
``w(c) = prod_p 1/mu(n(c_1..c_p))``, and ``Phi(n) = sum_c w(c)`` obeys the
balanced-fairness recursion (condition on the tail element)::
Phi(0) = 1, Phi(n) = (1/mu(n)) sum_{r: n_r>0} Phi(n - e_r)
Conditioning the same way and using
``sir_r(c) = sir_r(c_1..c_{|n|-1}) + [c_{|n|} = r] * 1{mu(n) > mu(n - e_r)}``
gives the companion recursion for ``Xi_r(n) = sum_c w(c) sir_r(c)``::
Xi_r(0) = 0
Xi_r(n) = (1/mu(n)) [ sum_{s: n_s>0} Xi_r(n - e_s)
+ 1{n_r > 0} 1{mu(n) > mu(n - e_r)} Phi(n - e_r) ]
Given n every ordering carries the same class-weight factor, so the
conditional law of the ordering is ``w(c)/Phi(n)`` and
``E[sir_r | n] = Xi_r(n)/Phi(n) =: g_r(n)``, a function of the count vector
alone. The station mean then follows from the count marginal pM as
``E[sir_r] = sum_n pM(n) g_r(n)``, or in normalizing-constant form from the
functional-server identity of :func:`pfqn_oi_fnc` applied to ``f(n) = g_r(n)``
(note ``g_r(0) = 0``, as required).
Parameters
----------
oirate : function ``mu(n)`` returning the OI total service rate for the
per-class count vector n (length R). ``mu(0)`` is taken as 0.
N : (R,) closed population vector, finite.
options : accepted for signature parity, currently unused.
Returns
-------
g : (prod(N+1), R) table, column-major over the lattice 0 <= n <= N, with
``g[1 + sum(n * stride), r] = E[sir_r | n]``.
Xi : (prod(N+1), R) table with the sir-weighted balance ``Xi_r(n)``.
Phi : (prod(N+1),) table with the OI balance function ``Phi(n)``.
See Also
--------
pfqn_oi_fnc, pfqn_oi_nc, pfqn_mvaoi, pfqn_mvaoi_marg
"""
if not callable(oirate):
raise ValueError('oirate must be a callable mu(n).')
N = np.asarray(N, dtype=int).ravel()
if np.any(N < 0):
raise ValueError('pfqn_oi_insvc requires finite nonnegative populations.')
R = N.size
shp = N + 1
stride = np.ones(R, dtype=int)
for d in range(1, R):
stride[d] = stride[d - 1] * shp[d - 1]
total = int(np.prod(shp))
subs = np.zeros((total, R), dtype=int)
for i in range(total):
li = i
for d in range(R):
subs[i, d] = li % shp[d]
li //= shp[d]
muv = np.zeros(total)
for i in range(total):
if subs[i].sum() > 0:
muv[i] = float(oirate(subs[i].copy()))
Phi = np.zeros(total)
Xi = np.zeros((total, R))
for i in range(total):
n = subs[i]
if n.sum() == 0:
Phi[i] = 1.0
continue
mun = muv[i]
if mun <= 0:
# Unreachable state (no server can serve this composition).
continue
sPhi = 0.0
sXi = np.zeros(R)
for s in range(R):
if n[s] > 0:
j = i - stride[s]
sPhi += Phi[j]
sXi += Xi[j, :]
Phi[i] = sPhi / mun
for r in range(R):
acc = sXi[r]
if n[r] > 0:
j = i - stride[r]
if mun > muv[j]:
acc += Phi[j] # the tail class-r job is in service
Xi[i, r] = acc / mun
g = np.zeros((total, R))
for i in range(total):
if Phi[i] > 0:
g[i, :] = Xi[i, :] / Phi[i]
return g, Xi, Phi
[docs]
def pfqn_oi_fnc(Phi: Sequence[float],
N: Optional[Sequence[int]] = None,
f: Optional[Callable] = None,
options=None) -> Tuple[Callable, np.ndarray, np.ndarray]:
"""OI generalization of the load-dependent functional server.
Builds an auxiliary OI station whose balance function Psi satisfies the
convolution identity ``(Psi * Phi)(n) = (1 + f(n)) Phi(n)``, then inverts
Psi to the FNC rate ``mu_f(n) = (sum_{r: n_r>0} Psi(n-e_r)) / Psi(n)``.
Parameters
----------
Phi : balance function of the existing OI station over the lattice, an
R-dimensional array of shape (N_1+1, ..., N_R+1), or a flat
column-major vector.
N : (R,) closed population vector. Optional when Phi is a full
R-dimensional array (then N = shape(Phi) - 1).
f : target queue-dependent function f(n), f(0)=0 (default f = sum(n)).
options : accepted for signature parity; unused.
Returns
-------
(muf, Psi, mu) : callable rate handle (Inf outside the lattice), FNC
balance array (lattice shape), and tabulated FNC rate array.
"""
Phi_arr = np.asarray(Phi, dtype=float)
if f is None:
f = lambda n: float(np.sum(n))
if N is None:
N = np.asarray(Phi_arr.shape) - 1
N = np.round(np.asarray(N, dtype=float)).astype(int).ravel()
R = N.size
shp = (N + 1).astype(int)
# Flatten Phi in column-major (Fortran) order to match the MATLAB Phi(:).
Phiv = Phi_arr.ravel(order='F')
total = int(np.prod(shp))
if Phiv.size != total:
raise ValueError('numel(Phi) must equal prod(N+1).')
# Column-major strides and decoded subscripts.
stride = np.ones(R, dtype=int)
for d in range(1, R):
stride[d] = stride[d - 1] * shp[d - 1]
subs = np.zeros((total, R), dtype=int)
for i in range(total):
li = i
for d in range(R):
subs[i, d] = li % shp[d]
li //= shp[d]
# Step 1: deconvolve (Psi * Phi)(n) = (1 + f(n)) Phi(n) for Psi.
Psiv = np.zeros(total)
for i in range(total):
n = subs[i, :]
acc = (1.0 + f(n)) * Phiv[i]
for j in range(i):
k = subs[j, :]
if np.all(k <= n):
idx = int(np.sum((n - k) * stride)) # linear index of n-k
acc -= Psiv[j] * Phiv[idx]
Psiv[i] = acc
# Step 2: balanced-fairness inversion of Psi to the FNC rate.
muv = np.full(total, np.inf)
for i in range(total):
n = subs[i, :]
if np.all(n == 0):
muv[i] = 0.0 # empty state
continue
denom = Psiv[i]
if denom == 0:
muv[i] = np.inf # non-physical / undefined rate
continue
num = 0.0
for r in range(R):
if n[r] > 0:
num += Psiv[i - stride[r]]
muv[i] = num / denom
if R == 1:
Psi = Psiv.copy()
mu = muv.copy()
else:
Psi = Psiv.reshape(tuple(shp), order='F')
mu = muv.reshape(tuple(shp), order='F')
muf = _make_oi_fnc_eval(muv, shp, stride)
return muf, Psi, mu
def _make_oi_fnc_eval(muv: np.ndarray, shp: np.ndarray, stride: np.ndarray) -> Callable:
"""Return the tabulated-lattice rate handle muf(n) (Inf out of lattice)."""
def oi_fnc_eval(n):
n = np.round(np.asarray(n, dtype=float)).astype(int).ravel()
if n.size != shp.size or np.any(n < 0) or np.any(n > shp - 1):
return np.inf
return float(muv[int(np.sum(n * stride))])
return oi_fnc_eval