"""
Loss Network Analysis via Erlang Fixed Point.
Native Python implementations of fixed-point algorithms for analyzing
loss networks using Erlang formulas. Provides efficient computation
of blocking probabilities for circuit-switched networks.
References:
Kelly, F.P., "Blocking probabilities in large circuit-switched networks."
Advances in Applied Probability, 1986.
Kaufman, J.S., "Blocking in a shared resource environment." IEEE Trans.
Communications, 1981; Roberts, J.W., "A service system with heterogeneous
user requirements", 1981 (the per-link occupancy recursion).
"""
import numpy as np
from scipy import special
from typing import Tuple
def _erlang_b(nu: float, C: int) -> float:
"""
Calculate the Erlang B formula for blocking probability.
The Erlang B formula computes the blocking probability in a loss system
where calls arrive according to a Poisson process.
Args:
nu: Traffic intensity (offered load = arrival rate * service time).
C: Number of servers (capacity).
Returns:
Blocking probability.
"""
if C < 0:
return 1.0
if nu <= 0:
return 0.0
# Use log-space computation for numerical stability
# log(nu^C / C!) - log(sum_{i=0}^C nu^i / i!)
log_terms = np.array([i * np.log(nu) - special.gammaln(i + 1) for i in range(C + 1)])
log_max = np.max(log_terms)
log_sum = log_max + np.log(np.sum(np.exp(log_terms - log_max)))
log_block = C * np.log(nu) - special.gammaln(C + 1) - log_sum
return np.exp(log_block)
def _rat_den(x: float, tol: float = 1e-12) -> int:
"""Denominator of the continued-fraction approximation of x within absolute tol, as MATLAB rat(x, tol)."""
if not np.isfinite(x) or x == np.floor(x):
return 1
p0, q0, p1, q1 = 0, 1, 1, 0
y = float(x)
for _ in range(64):
a = np.floor(y)
p0, p1 = p1, int(a) * p1 + p0
q0, q1 = q1, int(a) * q1 + q0
if abs(x - p1 / q1) <= tol or y == a:
break
y = 1.0 / (y - a)
return abs(q1)
def _kaufman_roberts_blocking(rho: np.ndarray, a: np.ndarray, c: int) -> np.ndarray:
"""Per-route blocking on one link of capacity c, routes asking a[r] units at offered load rho[r]:
n q(n) = sum_r a[r] rho[r] q(n - a[r]); route r is blocked when fewer than a[r] units are free."""
R = len(rho)
q = np.zeros(c + 1)
q[0] = 1.0
for n in range(1, c + 1):
s = 0.0
for r in range(R):
if a[r] > 0 and n - a[r] >= 0:
s += a[r] * rho[r] * q[n - a[r]]
q[n] = s / n
if q[n] > 1e250:
q = q / q[n] # rescale; only the normalized distribution is used
q = q / np.sum(q)
b = np.zeros(R)
for r in range(R):
if a[r] > 0:
b[r] = float(np.sum(q[max(0, c - a[r] + 1):]))
return b
[docs]
def lossn_erlangfp(nu: np.ndarray, A: np.ndarray, c: np.ndarray,
tol: float = 1e-10, max_iter: int = 10000
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, int]:
"""
Erlang fixed point (reduced-load) approximation for loss networks.
Calls (jobs) on route (class) r arrive according to Poisson rate nu_r.
Call service times on route r have unit mean. The link capacity
requirements are::
sum_r A[j,r] * n[j,r] <= c[j]
for all links j, where n[j,r] counts calls on route r on link j.
Link j sees route r offered at nu_r thinned by the blocking of r on every
OTHER link, and its occupancy distribution is the Kaufman-Roberts one
(Kaufman 1981, Roberts 1981), so B[j,r] = P(more than c[j]-A[j,r] units
busy). Unit requirements reduce this to Kelly's Erlang fixed point. The
earlier form fed the multi-unit load to Erlang B, which counted a call of
size a as a single unit. Mirrors MATLAB lossn_erlangfp.m.
Args:
nu: Arrival rates vector (R,) for each route.
A: Capacity requirement matrix (J, R).
c: Capacity vector (J,).
tol: Convergence tolerance (max abs change of the blocking matrix).
max_iter: Maximum iterations.
Returns:
Tuple of (qlen, loss, eblock, niter) where:
- qlen: Mean queue-length (carried load) for each route (R,)
- loss: Loss probability for each route (R,)
- eblock: Blocking probability of route r on link j (J, R)
- niter: Number of iterations
Example:
>>> nu = np.array([0.3, 0.1])
>>> A = np.array([[1, 1], [1, 4]]) # 2 links, 2 routes
>>> c = np.array([1, 3])
>>> qlen, loss, eblock, niter = lossn_erlangfp(nu, A, c)
"""
nu = np.asarray(nu, dtype=np.float64).ravel()
A = np.asarray(A, dtype=np.float64)
c = np.asarray(c, dtype=np.float64).ravel()
R = len(nu) # Number of routes
J = len(c) # Number of links
if A.shape != (J, R):
raise ValueError(f"A must have shape ({J}, {R}), got {A.shape}")
# A fractional row is the same constraint scaled to integers: multiply row j by the least common
# denominator of A[j,:] and c[j], which leaves the admissible set unchanged and gives the
# recursion whole units to count
from math import gcd
Aint = np.zeros((J, R), dtype=int)
Cint = np.zeros(J, dtype=int)
for j in range(J):
m = 1
for x in list(A[j, :]) + [c[j]]:
d = _rat_den(float(x))
m = m * d // gcd(m, d)
Aint[j, :] = np.round(m * A[j, :]).astype(int)
Cint[j] = int(np.floor(m * c[j] + 1e-9))
from ..da import da_fpi
def rl_sweep(B, itnum):
Bnew = B.copy()
for j in range(J):
rho = nu.copy()
for r in range(R):
for i in range(J):
if i != j and Aint[i, r] > 0:
rho[r] = rho[r] * (1 - B[i, r])
Bnew[j, :] = _kaufman_roberts_blocking(rho, Aint[j, :], int(Cint[j]))
return Bnew, B
E, niter, _ = da_fpi(rl_sweep, np.zeros((J, R)), max_iter, tol,
norm=lambda xn, xr: float(np.max(np.abs(xn - xr))),
nanstop=True)
carried = np.prod(1 - E, axis=0)
qlen = nu * carried
loss = 1 - carried # blocking probability = 1 - carried/offered
return qlen, loss, E, niter
[docs]
def erlang_b(offered_load: float, servers: int) -> float:
"""
Compute Erlang B blocking probability.
The Erlang B formula gives the probability that an arriving call
is blocked in an M/M/c/c loss system.
Args:
offered_load: Traffic intensity (arrival rate * service time).
servers: Number of servers.
Returns:
Blocking probability.
Example:
>>> erlang_b(10.0, 12) # Offered load 10 Erlang, 12 channels
0.1054...
"""
return _erlang_b(offered_load, servers)
[docs]
def erlang_c(offered_load: float, servers: int) -> float:
"""
Compute Erlang C delay probability.
The Erlang C formula gives the probability that an arriving call
must wait in an M/M/c queue.
Args:
offered_load: Traffic intensity (arrival rate * service time).
servers: Number of servers.
Returns:
Delay probability (probability of waiting).
Example:
>>> erlang_c(10.0, 12) # Offered load 10 Erlang, 12 agents
"""
if servers < 1:
return 1.0
if offered_load <= 0:
return 0.0
rho = offered_load / servers
if rho >= 1:
return 1.0
# Erlang C = c * B(c, A) / (c - A * (1 - B(c, A)))
B = _erlang_b(offered_load, servers)
C = servers * B / (servers - offered_load * (1 - B))
return min(1.0, max(0.0, C))
__all__ = [
'lossn_erlangfp',
'erlang_b',
'erlang_c',
]