"""
SN Getter Functions for Parameter Extraction.
Native Python implementations for extracting parameters from
network structures including arrival rates, throughputs, and
product-form chain parameters.
Key functions:
sn_get_arvr_from_tput: Compute arrival rates from throughputs
sn_get_node_arvr_from_tput: Compute node arrival rates from throughputs
sn_get_node_tput_from_tput: Compute node throughputs from station throughputs
sn_get_product_form_chain_params: Extract chain-aggregated parameters
References:
Original MATLAB: matlab/src/api/sn/sn_get_*.m
"""
import numpy as np
from typing import Tuple, Optional
from dataclasses import dataclass
from .network_struct import NetworkStruct, NodeType
from .fluid_places import sn_is_fluid_place
from line_solver.api.io.errors import LineInputError
from line_solver.constants import GlobalConstants
@dataclass
class ChainParams:
"""Chain-aggregated product-form parameters."""
lambda_vec: np.ndarray # Chain arrival rates
D: np.ndarray # Chain service demands at queuing stations
N: np.ndarray # Chain populations
Z: np.ndarray # Chain think times
mu: np.ndarray # Load-dependent service capacity scaling
S: np.ndarray # Number of servers at queuing stations
V: np.ndarray # Chain visit ratios
def _sn_nodeparam(sn: NetworkStruct, ind: int):
"""Return the nodeparam entry of node ``ind``, or None."""
nodeparam = getattr(sn, 'nodeparam', None)
if nodeparam is None:
return None
if isinstance(nodeparam, dict):
return nodeparam.get(ind, None)
return nodeparam[ind] if ind < len(nodeparam) else None
def sn_pn_firing_rates(sn: NetworkStruct, TN: np.ndarray, tput_is_tokens: bool):
"""
Recover per-mode transition firing rates from the Place throughputs.
The firing rates of a Petri net are not carried by the network structure,
but they are determined by the Place throughputs together with the net
structure. Writing x for the vector of per-mode firing rates, two families
of equations hold at steady state, for every Place p and class k:
departure the sum over the modes consuming (p,k) of x, weighted by the
input arc multiplicity when tput_is_tokens is True and
unweighted when it is False, equals TN(p,k)
balance the sum over all modes of x times (produced minus consumed)
equals zero
The system is solved in least squares. That is deliberate: an exact solver
supplies throughputs that satisfy it exactly and the fit is then the exact
answer, whereas a simulator supplies estimates that satisfy it only up to
sampling error and the least-squares fit is the right estimator there. A
residual test would reject every simulated run.
Args:
sn: Network structure
TN: Average throughputs at stations (M x R)
tput_is_tokens: True when TN counts tokens, False when it counts firing events
Returns:
(x, consumed, produced, place_nodes) where x is the firing rate per
(transition, mode) pair and is None when undetermined, consumed and
produced are indexed (mode, place, class), and place_nodes holds the
node indices of the Places in the order used above.
References:
Original MATLAB: matlab/src/api/sn/sn_pn_firing_rates.m
"""
undetermined = (None, None, None, [])
R = sn.nclasses
if TN is None or np.size(TN) == 0:
return undetermined
TN = np.atleast_2d(np.asarray(TN, dtype=float))
nodetype = np.asarray(sn.nodetype)
place_nodes = [int(i) for i in np.where(nodetype == NodeType.PLACE)[0]]
trans_nodes = [int(i) for i in np.where(nodetype == NodeType.TRANSITION)[0]]
if len(place_nodes) == 0 or len(trans_nodes) == 0:
return undetermined
# see _kb/03-api-layer.md for rationale
if np.any(nodetype == NodeType.SOURCE) or np.any(nodetype == NodeType.SINK):
return undetermined
stateful_nodes = [int(i) for i in np.where(np.asarray(sn.isstateful).ravel() > 0)[0]]
rt_arr = np.asarray(sn.rt, dtype=float)
for pind in place_nodes:
if pind not in stateful_nodes:
return undetermined
sfp = stateful_nodes.index(pind)
for sfj, jnd in enumerate(stateful_nodes):
if sfj == sfp or nodetype[jnd] == NodeType.TRANSITION:
continue
block_out = rt_arr[sfp * R:(sfp + 1) * R, sfj * R:(sfj + 1) * R]
block_in = rt_arr[sfj * R:(sfj + 1) * R, sfp * R:(sfp + 1) * R]
if np.any(block_out > 0) or np.any(block_in > 0):
return undetermined
# Enumerate the (transition, mode) pairs: a mode is what carries a firing
# rate, and a transition may hold several.
mode_trans = []
mode_idx = []
mode_timed = []
for ind in trans_nodes:
param = _sn_nodeparam(sn, ind)
nmodes = getattr(param, 'nmodes', None) if param is not None else None
if nmodes is None:
return undetermined
timing = getattr(param, 'timingstrategies', None)
for m in range(int(nmodes)):
mode_trans.append(ind)
mode_idx.append(m)
# see _kb/03-api-layer.md for rationale
timed = True
if timing is not None and len(timing) > m:
entry = timing[m]
name = entry if isinstance(entry, str) else getattr(entry, 'name', str(entry))
timed = str(name).upper() != 'IMMEDIATE'
mode_timed.append(timed)
mode_timed = np.asarray(mode_timed, dtype=bool)
n_modes = len(mode_trans)
if n_modes == 0:
return undetermined
n_places = len(place_nodes)
consumed = np.zeros((n_modes, n_places, R))
produced = np.zeros((n_modes, n_places, R))
for mm in range(n_modes):
param = _sn_nodeparam(sn, mode_trans[mm])
enab = np.atleast_2d(np.asarray(param.enabling[mode_idx[mm]], dtype=float))
fire = np.atleast_2d(np.asarray(param.firing[mode_idx[mm]], dtype=float))
for pp, pind in enumerate(place_nodes):
for k in range(R):
consumed[mm, pp, k] = max(0.0, enab[pind, k])
produced[mm, pp, k] = max(0.0, fire[pind, k])
# see _kb/03-api-layer.md for rationale
n_eq = 2 * n_places * R
A = np.zeros((n_eq, n_modes))
b = np.zeros(n_eq)
row = 0
n_measured = 0
for pp, pind in enumerate(place_nodes):
ist = int(sn.nodeToStation[pind])
for k in range(R):
if tput_is_tokens:
arow = np.array(consumed[:, pp, k], dtype=float)
else:
arow = (consumed[:, pp, k] > 0).astype(float)
arow[~mode_timed] = 0.0
if np.any(arow != 0.0):
A[row, :] = arow
b[row] = TN[ist, k] if ist >= 0 else 0.0
row += 1
n_measured += 1
A[row, :] = produced[:, pp, k] - consumed[:, pp, k]
b[row] = 0.0
row += 1
# With no measured row the system is homogeneous and pinv returns the zero
# vector, which would report every Place as idle. Keep what the caller had.
if n_measured == 0:
return undetermined
A = A[:row, :]
b = b[:row]
xfit = np.linalg.pinv(A) @ b
# A negative firing rate means the net structure was not read as intended;
# reporting a rate that cannot occur would be worse than reporting nothing.
if np.any(xfit < -1e-6 * max(1.0, float(np.max(np.abs(xfit))) if xfit.size else 1.0)):
return undetermined
return xfit, consumed, produced, place_nodes
def sn_pn_avg_rates(sn: NetworkStruct, QN: np.ndarray, TN: np.ndarray,
AN: Optional[np.ndarray] = None,
RN: Optional[np.ndarray] = None):
"""
Place throughput, arrival rate and response time in tokens.
A Place is a station and a token is the job it holds, so a firing that
consumes two tokens is two departures, not one. The CTMC and SSA analyzers
count firing events instead, which for unit arc multiplicities is the same
number and for weighted arcs is not: the reported throughput is then not a
token rate, and QLen over it is not a sojourn time.
Rescales the Place rows to tokens:
TN(p,k) tokens consumed from the Place per unit time
AN(p,k) tokens produced into the Place per unit time
RN(p,k) QN(p,k) / TN(p,k), Little's law over the Place
Rows that do not belong to a Place are returned untouched, so a mixed
Queue/Place model keeps its queueing metrics, and so are the rows of a FLUID
place, which report a level rather than tokens. When the firing rates cannot
be recovered from the throughputs the inputs are returned unchanged rather
than replaced by a guess.
Args:
sn: Network structure
QN: Average queue lengths, i.e. mean token counts at the Places
TN: Average throughputs at stations, counting firing events
AN: Average arrival rates at stations, as computed by the caller
RN: Average response times at stations, as computed by the caller
Returns:
(TN, AN, RN) with the Place rows expressed in tokens.
References:
Original MATLAB: matlab/src/api/sn/sn_pn_avg_rates.m
"""
if TN is None or np.size(TN) == 0:
return TN, AN, RN
if not np.any(np.asarray(sn.nodetype) == NodeType.PLACE):
return TN, AN, RN
# The analyzers hand over event counts, hence the False.
x, consumed, produced, place_nodes = sn_pn_firing_rates(sn, TN, False)
if x is None:
return TN, AN, RN
R = sn.nclasses
TN = np.array(np.atleast_2d(np.asarray(TN, dtype=float)), copy=True)
if AN is not None and np.size(AN) > 0:
AN = np.array(np.atleast_2d(np.asarray(AN, dtype=float)), copy=True)
if RN is not None and np.size(RN) > 0:
RN = np.array(np.atleast_2d(np.asarray(RN, dtype=float)), copy=True)
QN = np.atleast_2d(np.asarray(QN, dtype=float))
for pp, pind in enumerate(place_nodes):
ist = int(sn.nodeToStation[pind])
if ist < 0:
continue
# A FLUID place holds a level that no token arc touches; its row is the level and outflow the solver
# measured, and there is no token count to convert it to.
if sn_is_fluid_place(sn, int(pind)):
continue
for k in range(R):
tput = float(consumed[:, pp, k] @ x)
TN[ist, k] = tput
if AN is not None and np.size(AN) > 0:
AN[ist, k] = float(produced[:, pp, k] @ x)
if RN is not None and np.size(RN) > 0:
RN[ist, k] = QN[ist, k] / tput if tput > 0 else 0.0
return TN, AN, RN
[docs]
def sn_get_arvr_from_tput(sn: NetworkStruct, TN: np.ndarray,
TH: Optional[np.ndarray] = None) -> np.ndarray:
"""
Compute average arrival rates at stations from throughputs.
Calculates the average arrival rate at each station in steady-state
from the station throughputs and routing matrix.
Args:
sn: Network structure
TN: Average throughputs at stations (M x R)
TH: Throughput handles (optional)
Returns:
AN: Average arrival rates at stations (M x R)
References:
Original MATLAB: matlab/src/api/sn/sn_get_arvr_from_tput.m
"""
M = sn.nstations
R = sn.nclasses
if TN is None or len(TN) == 0:
return np.array([])
TN = np.atleast_2d(np.asarray(TN, dtype=float))
AN = np.zeros((M, R))
# Build mapping from stateful nodes to their position in rt matrix
stateful_nodes = np.where(sn.isstateful)[0]
n_stateful = len(stateful_nodes)
# Build throughput vector for all stateful nodes (stations have TN, others need computation)
TN_stateful = np.zeros((n_stateful, R))
ucache = [] # rt positions sf*R+k of the Cache retrieval-class flows, solved with the Routers below
for sf, ind in enumerate(stateful_nodes):
ist = sn.nodeToStation[ind]
if ist >= 0 and ist < TN.shape[0]:
# This stateful node is a station - use station throughput
TN_stateful[sf, :] = TN[ist, :]
elif sn.nodetype[ind] == NodeType.CACHE:
# This stateful node is not a station (e.g., Cache): compute the
# hit/miss class throughputs from the reference station throughput
# and the actual hit/miss probabilities
nodeparam = None
if getattr(sn, 'nodeparam', None) is not None:
if isinstance(sn.nodeparam, dict):
nodeparam = sn.nodeparam.get(ind, None)
elif ind < len(sn.nodeparam):
nodeparam = sn.nodeparam[ind]
if nodeparam is None:
continue
hitclass = getattr(nodeparam, 'hitclass', None)
missclass = getattr(nodeparam, 'missclass', None)
if hitclass is None or missclass is None:
continue
actualhitprob = getattr(nodeparam, 'actualhitprob', None)
actualmissprob = getattr(nodeparam, 'actualmissprob', None)
if actualhitprob is None or np.size(actualhitprob) == 0:
# Actual probabilities not yet computed - skip this cache;
# arrival rates will be computed later when they are available
continue
hitclass = np.atleast_1d(hitclass).astype(int)
missclass = np.atleast_1d(missclass).astype(int)
actualhitprob = np.atleast_1d(np.asarray(actualhitprob, dtype=float))
actualmissprob = np.atleast_1d(np.asarray(actualmissprob, dtype=float))
# Delayed-hit (retrieval) cache: a delayed hit leaves in the hit class, and each retrieval class carries
# missclass only so that a COMPLETED fetch switches to the miss class. Its entry must not be read as a
# second read class: with actualmissprob 0 or NaN there it overwrote the read class's miss flow.
retr_cls = set(int(k) for k in (getattr(nodeparam, 'retrieval_class_indices', None) or set()))
dhp = getattr(nodeparam, 'actualdelayedhitprob', None)
dhp = np.atleast_1d(np.asarray(dhp, dtype=float)) if dhp is not None and np.size(dhp) > 0 else None
# Find the chain and reference station for this cache
for c in range(sn.nchains):
inchain = np.atleast_1d(sn.inchain[c]).astype(int).ravel()
refstat = int(sn.refstat[c])
totalTput = float(np.sum(TN[refstat, inchain]))
# Set throughput for hit/miss classes (0-based: -1 means none)
for origClass in range(len(hitclass)):
if origClass in retr_cls or origClass not in inchain:
continue
hc = int(hitclass[origClass])
mc = int(missclass[origClass])
if 0 <= hc < R and origClass < len(actualhitprob) and np.isfinite(actualhitprob[origClass]):
ph = actualhitprob[origClass]
if dhp is not None and origClass < len(dhp) and np.isfinite(dhp[origClass]):
ph = ph + dhp[origClass]
TN_stateful[sf, hc] = totalTput * ph
if 0 <= mc < R and origClass < len(actualmissprob) and np.isfinite(actualmissprob[origClass]):
TN_stateful[sf, mc] = totalTput * actualmissprob[origClass]
# A fetch leaves the Cache at the rate it returns to it, so each retrieval class's (Cache, class) flow is an
# UNKNOWN of the flow balance below, fed by whatever routes it back (fetch stations, or a Router between).
# The open analyzers book the fetch flow to the READ class, which rt does not route from the Cache; they
# supply the fetch stations' ArvR themselves (solver_{mva,nc}_retrieval_analyzer).
ucache.extend(sf * R + rc for rc in sorted(retr_cls) if 0 <= rc < R)
# Materialize rt as dense numpy array once for vectorized operations
rt_arr = np.asarray(sn.rt)
rt_rows, rt_cols = rt_arr.shape
max_sf_r = min(n_stateful * R, rt_rows)
# Non-station, non-Cache stateful nodes (e.g., Router) have no throughput of their own. Solve the flow balance
# x_u = x_k*P_ku + x_u*P_uu exactly: a sweep in node order left a Router fed by a later-declared Router at zero.
unk = [sf for sf in range(n_stateful)
if sn.nodeToStation[stateful_nodes[sf]] < 0 and sn.nodetype[stateful_nodes[sf]] != NodeType.CACHE]
if unk or ucache:
ucols = np.array([sf * R + k for sf in unk for k in range(R)] + ucache, dtype=int)
is_u = np.zeros(n_stateful * R, dtype=bool)
is_u[ucols] = True
kcols = np.where(~is_u)[0]
TN_flat = TN_stateful.reshape(-1).copy()
# a Cache's read-class rows may carry NaN (state-dependent hit/miss switch); they are not flows of rt here
rt_solve = np.nan_to_num(rt_arr, nan=0.0)
b = TN_flat[kcols] @ rt_solve[np.ix_(kcols, ucols)]
A = np.eye(len(ucols)) - rt_solve[np.ix_(ucols, ucols)]
xu = np.linalg.solve(A.T, b)
TN_flat[ucols] = xu
TN_stateful = TN_flat.reshape(n_stateful, R)
# Compute arrival rates using stateful node throughputs and rt matrix
# Vectorized: AN[ist, :] = TN_flat @ rt_block
TN_flat = TN_stateful.reshape(-1)[:max_sf_r]
for ist in range(M):
ind_ist = sn.stationToNode[ist]
if sn.nodetype[ind_ist] == NodeType.SOURCE:
AN[ist, :] = 0
else:
sf_ist_arr = np.where(stateful_nodes == ind_ist)[0]
if len(sf_ist_arr) == 0:
continue
sf_ist = sf_ist_arr[0]
col_start = sf_ist * R
col_end = min(col_start + R, rt_cols)
if col_start < rt_cols:
rt_block = rt_arr[:max_sf_r, col_start:col_end]
AN[ist, :col_end - col_start] += TN_flat @ rt_block
# Fork-join special handling: delegate to node-level arrival rate computation
# MATLAB's sn_get_arvr_from_tput lines 117-122
if hasattr(sn, 'fj') and sn.fj is not None and np.any(sn.fj):
ANn = sn_get_node_arvr_from_tput(sn, TN, TH, AN)
for ist in range(M):
ind = sn.stationToNode[ist]
if ind >= 0 and ind < ANn.shape[0]:
AN[ist, :] = ANn[ind, :]
# see _kb/03-api-layer.md for rationale
if AN.size > 0 and np.any(np.asarray(sn.nodetype) == NodeType.PLACE):
x, _, produced, place_nodes = sn_pn_firing_rates(sn, TN, True)
if x is not None:
for pp, pind in enumerate(place_nodes):
ist = int(sn.nodeToStation[pind])
if ist >= 0:
for k in range(R):
AN[ist, k] = float(produced[:, pp, k] @ x)
return AN
def sn_arvr_state_dependent_balance(sn: NetworkStruct, AN: np.ndarray, TN: np.ndarray,
WN: Optional[np.ndarray] = None):
"""Arrival rates at stations fed by a JSQ or SQ(d) dispatcher, from flow balance.
sn_get_arvr_from_tput reads the routing matrix, which for a state-dependent
dispatcher holds only the static RAND split: under JSQ it reported ArvR
0.5/0.5 at two stations whose throughputs were 0.57/0.42. At a station a
job leaves only by completing service (no finite buffer, no impatience, no
retrial) arrivals and departures balance in steady state, so ArvR = Tput
there whatever the routing. Residence times derived from the same static
visit ratios are rescaled by Tput over the STATIC arrival rate, recomputed
here from TN: the AN passed in may be measured (JMT, LDES, SSA), and scaling
by Tput/measured ArvR (~1) left WN at the static visit ratio. Other stations
are left as computed: their losses are not visible here.
Returns:
(AN, WN) with the corrected entries; WN is None when none was given.
"""
routing = getattr(sn, 'routing', None)
rtnodes = getattr(sn, 'rtnodes', None)
if routing is None or rtnodes is None or AN is None or TN is None:
return AN, WN
routing = np.asarray(routing)
rtnodes = np.asarray(rtnodes, dtype=float)
I, K = int(sn.nnodes), int(sn.nclasses)
dyn = {4, 6} # RoutingStrategy.JSQ, RoutingStrategy.SQ
fed = set()
for ind in range(min(I, routing.shape[0])):
for r in range(min(K, routing.shape[1])):
val = routing[ind, r]
if int(getattr(val, 'value', val)) not in dyn:
continue
for jnd in range(I):
for s in range(K):
if rtnodes[ind * K + r, jnd * K + s] > 0 and sn.isstation[jnd]:
fed.add((int(sn.nodeToStation[jnd]), s))
if not fed:
return AN, WN
AN = np.array(AN, dtype=float, copy=True)
WN = None if WN is None else np.array(WN, dtype=float, copy=True)
cap = np.asarray(getattr(sn, 'cap', np.full(int(sn.nstations), np.inf)), dtype=float).ravel()
classcap = getattr(sn, 'classcap', None)
imp = getattr(sn, 'impatienceClass', None)
rproc = getattr(sn, 'retrialProc', None)
ANstatic = None if WN is None else sn_get_arvr_from_tput(sn, TN, TN)
for (ist, s) in fed:
if ist >= AN.shape[0] or s >= AN.shape[1]:
continue
if ist < cap.size and np.isfinite(cap[ist]):
continue
if classcap is not None and np.isfinite(np.asarray(classcap, dtype=float)[ist, s]):
continue
if imp is not None and np.size(imp) > 0 and np.asarray(imp)[ist, s] != 0:
continue
if rproc is not None and ist < len(rproc) and rproc[ist][s] is not None:
continue
AN[ist, s] = TN[ist, s]
if ANstatic is not None and ANstatic[ist, s] > 0:
WN[ist, s] = WN[ist, s] * TN[ist, s] / ANstatic[ist, s]
return AN, WN
def sn_get_node_arvr_from_tput(sn: NetworkStruct, TN: np.ndarray,
TH: Optional[np.ndarray] = None,
AN: Optional[np.ndarray] = None) -> np.ndarray:
"""
Compute node arrival rates from station throughputs.
Handles:
- Station nodes: Uses station arrival rates directly
- Cache nodes: Only requesting classes arrive (not hit/miss classes)
- Non-station nodes (ClassSwitch, Sink): Uses nodevisits-based computation
Args:
sn: Network structure
TN: Station throughputs (M x R)
TH: Throughput handles (optional)
AN: Station arrival rates (optional, computed if not provided)
Returns:
ANn: Node arrival rates (I x R)
References:
Original MATLAB: matlab/src/api/sn/sn_get_node_arvr_from_tput.m
"""
I = sn.nnodes
M = sn.nstations
C = sn.nchains
R = sn.nclasses
ANn = np.zeros((I, R))
if TN is None or np.size(TN) == 0:
return ANn
TN = np.atleast_2d(np.asarray(TN, dtype=float))
if AN is None:
AN = sn_get_arvr_from_tput(sn, TN, TH)
AN = np.atleast_2d(np.asarray(AN, dtype=float))
# Station nodes take the station arrival rates
for ist in range(M):
ind = int(sn.stationToNode[ist])
if 0 <= ind < I:
ANn[ind, :] = AN[ist, :]
with np.errstate(divide='ignore', invalid='ignore'):
for ind in range(I):
node_type = sn.nodetype[ind]
if node_type == NodeType.SOURCE:
continue
for c in range(C):
if c not in sn.inchain:
continue
inchain = np.atleast_1d(sn.inchain[c]).astype(int).ravel()
refstat = int(sn.refstat[c])
visits_c = sn.visits[c] if c in sn.visits else None
nodevisits_c = sn.nodevisits[c] if c in sn.nodevisits else None
if visits_c is None or nodevisits_c is None:
continue
refsf = int(sn.stationToStateful[refstat])
denom = float(np.sum(visits_c[refsf, inchain]))
totalTput = float(np.sum(TN[refstat, inchain]))
for r in inchain:
if r >= R:
continue
if node_type == NodeType.CACHE:
# For cache nodes, only the requesting class arrives;
# hit/miss classes don't arrive - they leave
hitclass = []
missclass = []
node_param = None
if sn.nodeparam is not None:
if isinstance(sn.nodeparam, dict):
node_param = sn.nodeparam.get(ind, None)
elif ind < len(sn.nodeparam):
node_param = sn.nodeparam[ind]
if node_param is not None:
if getattr(node_param, 'hitclass', None) is not None:
hitclass = np.atleast_1d(node_param.hitclass).ravel()
if getattr(node_param, 'missclass', None) is not None:
missclass = np.atleast_1d(node_param.missclass).ravel()
if not ((r in hitclass) or (r in missclass)):
ANn[ind, r] = (nodevisits_c[ind, r] / denom) * totalTput
else:
# For station nodes, the arrival rate is already set from AN;
# only use nodevisits for non-station nodes (ClassSwitch, Sink, ...)
if not (sn.nodeToStation[ind] >= 0):
ANn[ind, r] = (nodevisits_c[ind, r] / denom) * totalTput
# Replace NaN with 0
ANn = np.nan_to_num(ANn, nan=0.0)
return ANn
def _cache_downstream_arvr(sn: NetworkStruct, ANn: np.ndarray, TNn: np.ndarray,
TN: np.ndarray, TH: Optional[np.ndarray] = None):
"""MATLAB NetworkSolver/getAvgNode.m post-pass, run after the two node getters.
A cache leaves its hit/miss routing as NaN until it is solved, so the node visits split the read
flow evenly; re-derive the arrivals at the non-station nodes a cache feeds (ClassSwitch, Sink, a
further Cache) from its solved hit/miss throughput along rtnodes, then the node throughputs from
those arrivals. A cache feeding a cache needs one pass per level, so iterate to the fixed point.
Returns (ANn, TNn)."""
I = sn.nnodes
R = sn.nclasses
rtn = np.asarray(sn.rtnodes)
for _ in range(I):
changed = False
for cacheInd in range(I):
if sn.nodetype[cacheInd] != NodeType.CACHE:
continue
cache_param = _sn_nodeparam(sn, cacheInd)
if cache_param is None:
continue
hitclass = np.atleast_1d(getattr(cache_param, 'hitclass', [])).astype(int).ravel()
missclass = np.atleast_1d(getattr(cache_param, 'missclass', [])).astype(int).ravel()
outs = np.unique(np.concatenate([hitclass[hitclass >= 0], missclass[missclass >= 0]]))
if outs.size == 0:
continue
for ind in range(I):
if ind == cacheInd or int(sn.nodeToStation[ind]) >= 0:
continue
for s in range(R):
p = rtn[cacheInd * R + outs, ind * R + s]
fed = p > 0
if np.any(fed):
arv = float(TNn[cacheInd, outs[fed]] @ p[fed])
if abs(ANn[ind, s] - arv) > GlobalConstants.FineTol * max(1.0, arv):
ANn[ind, s] = arv
changed = True
if not changed:
break
TNn = sn_get_node_tput_from_tput(sn, TN, TH, ANn)
# The station arrivals were derived from routing that still split each read evenly; with the
# node throughputs settled, re-derive them by flow balance into each station.
if any(sn.nodetype[i] == NodeType.CACHE for i in range(I)):
flow = np.nan_to_num(rtn, nan=0.0)
inflow = TNn.reshape(1, I * R) @ flow
for ist in range(sn.nstations):
ind = int(sn.stationToNode[ist])
if sn.nodetype[ind] != NodeType.SOURCE:
ANn[ind, :] = inflow[0, ind * R:(ind + 1) * R]
return ANn, TNn
def sn_get_node_tput_from_tput(sn: NetworkStruct, TN: np.ndarray,
TH: Optional[np.ndarray] = None,
ANn: Optional[np.ndarray] = None) -> np.ndarray:
"""
Compute node throughputs from station throughputs.
Handles:
- Station nodes: Uses station throughputs directly
- Cache nodes: Uses actual hit/miss probabilities if available
- Non-station nodes: Uses routing matrix (rtnodes) for computation
Args:
sn: Network structure
TN: Station throughputs (M x R)
TH: Throughput handles (optional)
ANn: Node arrival rates (optional, computed if not provided)
Returns:
TNn: Node throughputs (I x R)
References:
Original MATLAB: matlab/src/api/sn/sn_get_node_tput_from_tput.m
"""
I = sn.nnodes
M = sn.nstations
C = sn.nchains
R = sn.nclasses
TNn = np.zeros((I, R))
if TN is None or np.size(TN) == 0:
return TNn
TN = np.atleast_2d(np.asarray(TN, dtype=float))
if ANn is None:
ANn = sn_get_node_arvr_from_tput(sn, TN, TH)
has_ANn = ANn is not None and np.size(ANn) > 0
if has_ANn:
ANn = np.atleast_2d(np.asarray(ANn, dtype=float))
def _nodeparam(ind):
if sn.nodeparam is None:
return None
if isinstance(sn.nodeparam, dict):
return sn.nodeparam.get(ind, None)
return sn.nodeparam[ind] if ind < len(sn.nodeparam) else None
def _prob_vec(node_param, name):
v = getattr(node_param, name, None)
if v is None or np.size(v) == 0:
return None
return np.atleast_1d(np.asarray(v, dtype=float)).ravel()
# First pass: Cache nodes, from the actual hit/miss ratios when available
with np.errstate(divide='ignore', invalid='ignore'):
for ind in range(I):
if sn.nodetype[ind] != NodeType.CACHE:
continue
node_param = _nodeparam(ind)
if node_param is None:
continue
hitclass = np.atleast_1d(getattr(node_param, 'hitclass', [])).astype(int).ravel()
missclass = np.atleast_1d(getattr(node_param, 'missclass', [])).astype(int).ravel()
actual_hit_prob = _prob_vec(node_param, 'actualhitprob')
actual_miss_prob = _prob_vec(node_param, 'actualmissprob')
actual_delayed_hit_prob = _prob_vec(node_param, 'actualdelayedhitprob')
for c in range(C):
if c not in sn.inchain:
continue
inchain = np.atleast_1d(sn.inchain[c]).astype(int).ravel()
refstat = int(sn.refstat[c])
totalTput = float(np.sum(TN[refstat, inchain]))
for r in inchain:
if r >= R:
continue
if actual_hit_prob is not None:
# Use actual hit/miss probabilities from solver result.
# Delayed hits (retrieval system) depart as the hit class,
# so the hit-class throughput is (hit + delayed hit) * arrivals.
if actual_delayed_hit_prob is None:
actual_delayed_hit_prob = np.zeros_like(actual_hit_prob)
for origClass in range(len(hitclass)):
if has_ANn:
arvTput = ANn[ind, origClass]
elif origClass in inchain:
arvTput = TN[refstat, origClass]
else:
arvTput = 0.0 # read class outside this chain
if (hitclass[origClass] == r and origClass < len(actual_hit_prob)
and not np.isnan(actual_hit_prob[origClass])):
dh = 0.0
if (origClass < len(actual_delayed_hit_prob)
and not np.isnan(actual_delayed_hit_prob[origClass])):
dh = actual_delayed_hit_prob[origClass]
TNn[ind, r] += arvTput * (actual_hit_prob[origClass] + dh)
elif (actual_miss_prob is not None and missclass[origClass] == r
and origClass < len(actual_miss_prob)
and not np.isnan(actual_miss_prob[origClass])):
TNn[ind, r] += arvTput * actual_miss_prob[origClass]
else:
# Fallback to nodevisits-based calculation
if (r in hitclass) or (r in missclass):
refsf = int(sn.stationToStateful[refstat])
denom = float(np.sum(sn.visits[c][refsf, inchain]))
TNn[ind, r] = (sn.nodevisits[c][ind, r] / denom) * totalTput
# Copy station throughputs directly to station nodes
for ist in range(M):
ind = int(sn.stationToNode[ist])
if 0 <= ind < I:
TNn[ind, :] = TN[ist, :]
rtnodes = np.asarray(sn.rtnodes) if sn.rtnodes is not None else None
def _rt(ind, r, jnd, s):
if rtnodes is None or rtnodes.size == 0:
return 0.0
i0 = ind * R + r
j0 = jnd * R + s
if i0 < rtnodes.shape[0] and j0 < rtnodes.shape[1]:
return rtnodes[i0, j0]
return 0.0
def _cache_output_class(ind, r):
# hit, miss or retrieval-pending class of the Cache (MATLAB isReceived)
node_param = _nodeparam(ind)
if node_param is None:
return False
for name in ('hitclass', 'missclass'):
v = getattr(node_param, name, None)
if v is not None and r in np.atleast_1d(v).astype(int).ravel():
return True
return r in getattr(node_param, 'retrieval_class_indices', set())
cache_output_classes = {r for ind in range(I) if sn.nodetype[ind] == NodeType.CACHE
for r in range(R) if _cache_output_class(ind, r)}
# Precompute sum_{jnd} rtnodes(ind, r, jnd, s_) and sum_{jnd != ind} rtnodes(ind, r, jnd, s_)
rtn_arr = np.asarray(sn.rtnodes)
if rtn_arr.shape == (I * R, I * R):
rtn_4d = rtn_arr.reshape(I, R, I, R)
rt_sum_all = rtn_4d.sum(axis=2) # (I, R, R)
rt_diag = np.diagonal(rtn_4d, axis1=0, axis2=2).transpose(2, 0, 1) # (I, R, R)
rt_sum_neq = rt_sum_all - rt_diag
else:
rt_sum_all = None
rt_sum_neq = None
# Non-station nodes (ClassSwitch, Router, Cache, Join, ...) from the node routing matrix
for ind in range(I):
node_type = sn.nodetype[ind]
for c in range(C):
if c not in sn.inchain:
continue
inchain = np.atleast_1d(sn.inchain[c]).astype(int).ravel()
visits_c = sn.visits[c] if c in sn.visits else None
for r in inchain:
if r >= R:
continue
anystateful = visits_c is not None and np.any(visits_c[:, r] != 0)
# MATLAB's ind->ind hit/miss self-loop gives a cache output class visits
# at the Cache; Python's stochastic complement folds them away (rtnodes above)
if not anystateful and r in cache_output_classes:
anystateful = True
if not anystateful:
continue
if node_type != NodeType.SINK and node_type != NodeType.JOIN:
for s_ in inchain:
if s_ >= R:
continue
if node_type == NodeType.SOURCE:
ist = int(sn.nodeToStation[ind])
TNn[ind, s_] = TN[ist, s_]
elif node_type == NodeType.CACHE:
# MATLAB getRoutingMatrix leaves a read-class row of a Cache
# only the ind->ind hit/miss self-loop, which ind~=jnd drops;
# Python's rtnodes routes that row downstream instead, so the
# first pass alone accounts for it (else hits count twice).
if has_ANn and _cache_output_class(ind, r):
if rt_sum_neq is not None:
TNn[ind, s_] += ANn[ind, r] * rt_sum_neq[ind, r, s_]
else:
for jnd in range(I):
if ind != jnd:
TNn[ind, s_] += ANn[ind, r] * _rt(ind, r, jnd, s_)
else:
# For station nodes, throughput is already set from TN;
# only compute for non-station nodes (ClassSwitch, Router)
if not (sn.nodeToStation[ind] >= 0) and has_ANn:
if rt_sum_all is not None:
TNn[ind, s_] += ANn[ind, r] * rt_sum_all[ind, r, s_]
else:
for jnd in range(I):
TNn[ind, s_] += ANn[ind, r] * _rt(ind, r, jnd, s_)
elif node_type == NodeType.JOIN:
if has_ANn:
for s_ in inchain:
if s_ >= R:
continue
if rt_sum_all is not None:
TNn[ind, s_] += ANn[ind, r] * rt_sum_all[ind, r, s_]
else:
for jnd in range(I):
TNn[ind, s_] += ANn[ind, r] * _rt(ind, r, jnd, s_)
# Replace NaN with 0
TNn = np.nan_to_num(TNn, nan=0.0)
return TNn
def sn_get_product_form_chain_params(sn: NetworkStruct) -> ChainParams:
"""
Extract product-form parameters aggregated by chain.
Extracts parameters from a network structure and aggregates them
by chain for use in product-form analysis methods.
Args:
sn: Network structure
Returns:
ChainParams with lambda_vec, D, N, Z, mu, S, V
References:
Original MATLAB: matlab/src/api/sn/sn_get_product_form_chain_params.m
"""
from .transforms import sn_get_product_form_params
from .demands import sn_get_demands_chain
# Get base product-form parameters
params = sn_get_product_form_params(sn)
# Get chain-aggregated demands
demands = sn_get_demands_chain(sn)
# Find queue and delay indices
# Handle both list and numpy array nodetype
nodetype = sn.nodetype if isinstance(sn.nodetype, np.ndarray) else np.array([nt.value if hasattr(nt, 'value') else nt for nt in sn.nodetype])
queue_indices = np.where(nodetype == NodeType.QUEUE.value)[0]
delay_indices = np.where(nodetype == NodeType.DELAY.value)[0]
# Initialize chain parameters
nchains = sn.nchains
lambda_chains = np.zeros(nchains)
for c in range(nchains):
chain_classes = sn.inchain[c]
lambda_chains[c] = np.nansum(params.lam[chain_classes])
# Extract demands at queue and delay stations
# Note: Lchain in Python corresponds to Dchain in MATLAB
D_chains = demands.Lchain[sn.nodeToStation[queue_indices], :]
Z_chains = demands.Lchain[sn.nodeToStation[delay_indices], :] if len(delay_indices) > 0 else np.zeros((0, nchains))
# Number of servers at queuing stations
S = sn.nservers[sn.nodeToStation[queue_indices]]
# Visit ratios: drop the rows of the Source/Join stations (MATLAB
# Vchain(sn.nodeToStation(ignoreIndex),:) = [])
V = demands.Vchain.copy()
ignore_indices = np.where((nodetype == NodeType.SOURCE.value) | (nodetype == NodeType.JOIN.value))[0]
if len(ignore_indices) > 0:
ignore_stations = [int(sn.nodeToStation[i]) for i in ignore_indices if int(sn.nodeToStation[i]) >= 0]
keep_stations = [i for i in range(V.shape[0]) if i not in ignore_stations]
V = V[keep_stations, :]
# Z is the (Mz x C) matrix of chain demands at the delay stations; with no
# delay station it is a zero vector shaped like N (MATLAB: Z = 0*N)
Z = np.asarray(Z_chains, dtype=float)
if Z.size == 0:
with np.errstate(invalid='ignore'): # 0*Inf = NaN for open chains, as in MATLAB
Z = 0 * np.asarray(demands.Nchain, dtype=float)
return ChainParams(
lambda_vec=lambda_chains,
D=D_chains,
N=demands.Nchain,
Z=Z,
mu=params.mu,
S=S,
V=V
)
def sn_set_routing_prob(sn: NetworkStruct, from_stateful: int, from_class: int,
to_stateful: int, to_class: int, prob: float,
auto_refresh: bool = False) -> NetworkStruct:
"""
Set a routing probability between two stateful node-class pairs.
Updates a single entry in the rt matrix.
Args:
sn: Network structure
from_stateful: Source stateful node index (0-based)
from_class: Source class index (0-based)
to_stateful: Destination stateful node index (0-based)
to_class: Destination class index (0-based)
prob: Routing probability [0, 1]
auto_refresh: If True, refresh visit ratios (default False)
Returns:
Modified network structure
References:
Original MATLAB: matlab/src/api/sn/sn_set_routing_prob.m
"""
K = sn.nclasses
# Calculate indices in rt matrix
from_idx = from_stateful * K + from_class
to_idx = to_stateful * K + to_class
# Update rt matrix
sn.rt[from_idx, to_idx] = prob
# Auto-refresh visit ratios if requested
if auto_refresh:
from .transforms import sn_refresh_visits
sn_refresh_visits(sn)
return sn
def sn_region_members(sn, f, Rmat, memvec) -> np.ndarray:
"""Station membership mask of finite capacity region f, as a bool array of length M.
Membership is read from sn.regionmembers[f], which the region refresh records
directly from the region's node list. It cannot be derived from sn.region[f]:
-1 there means "unbounded", which is indistinguishable from "not a member", so a
region constrained only by regionlincon (or only by a memory budget) reads as
empty and is silently ignored.
Rmat and memvec provide the legacy derivation, used only for an sn built before
regionmembers existed (for instance one deserialised from an older model file).
That derivation carries the ambiguity above and is not equivalent.
"""
Rmat = np.atleast_2d(np.asarray(Rmat, dtype=float))
M = Rmat.shape[0]
members = getattr(sn, 'regionmembers', None)
if members is not None and len(members) > f and members[f] is not None:
mask = np.asarray(members[f]).ravel().astype(bool)
if mask.size >= M:
return mask[:M]
legacy = np.any(Rmat != -1, axis=1)
if memvec is not None:
mv = np.asarray(memvec, dtype=float).ravel()
if mv.size >= M:
legacy = legacy | (mv[:M] != -1)
return legacy
def sn_map_modulation(sn):
"""Collect the (D0,D1) modulation records of every non-renewal process.
A MAP with matrices (D0,D1) is a Poisson-like point process modulated by the
CTMC with generator Q = D0 + D1 (the phase process), whose conditional
intensity in phase k is lambda(k) = sum_j D1(k,j). This returns one record
per modulating process, so that a solver-agnostic transformation can replace
each of them by a random-environment stage set (see api.io.map2renv).
Only processes declared as MAP, MMPP2 or MMAP are reported: every other
distribution is stored in sn.proc in (D0,D1) form as well (Erlang, Coxian,
APH, ...), but those are renewal processes that carry no modulation and are
supported natively by the phase-type solvers.
Marked processes (MMAP) at a Source are reported as a single record whose
'classes' entry lists every marked class, since all marks share one phase
process; the per-class intensity comes from the mark-specific D1 matrices.
Mirrors matlab/src/api/sn/sn_map_modulation.m.
Args:
sn: NetworkStruct object
Returns:
List of dicts with keys ist, node, arrival, classes, D0, D1 (list, one
per entry of classes), order, is_mmpp.
"""
from .network_struct import NodeType
from ...constants import ProcessType
mods = []
procid = getattr(sn, 'procid', None)
if procid is None or sn.proc is None:
return mods
mod_types = (ProcessType.MAP, ProcessType.MMPP2, ProcessType.MMAP,
ProcessType.MPH)
markidx = getattr(sn, 'markidx', None)
def _map_of(ist, r):
proc_st = sn.proc[ist] if ist < len(sn.proc) else None
if proc_st is None or r >= len(proc_st) or proc_st[r] is None:
return None
map_ir = proc_st[r]
if not isinstance(map_ir, (list, tuple)) or len(map_ir) < 2 or map_ir[0] is None:
return None
if np.any(np.isnan(np.asarray(map_ir[0], dtype=np.float64))):
return None # disabled
return map_ir
def _mark_index(ist, r):
if markidx is None:
return -1
arr = np.asarray(markidx)
if ist >= arr.shape[0] or r >= arr.shape[1]:
return -1
return int(arr[ist, r])
def _is_diagonal(D1):
D1 = np.asarray(D1, dtype=np.float64)
off = D1 - np.diag(np.diag(D1))
return np.linalg.norm(off, 'fro') <= 1e-14 * max(1.0, np.linalg.norm(D1, 'fro'))
for ist in range(sn.nstations):
nd = int(sn.stationToNode[ist])
arrival = sn.nodetype[nd] == NodeType.SOURCE
done = [False] * sn.nclasses
for r in range(sn.nclasses):
if done[r]:
continue
pt = procid[ist][r] if not hasattr(procid, 'shape') else procid[ist, r]
if pt not in mod_types:
continue
map_ir = _map_of(ist, r)
if map_ir is None:
continue
if _mark_index(ist, r) > 0:
# MMAP: one phase process shared by every marked class of the
# station, one D1 block per mark
marked = [rr for rr in range(sn.nclasses) if _mark_index(ist, rr) > 0]
carrier = min(marked, key=lambda rr: _mark_index(ist, rr))
map_c = _map_of(ist, carrier)
if map_c is None or len(map_c) < 2 + len(marked):
raise LineInputError(
"The marked arrival process at station %d carries %d mark matrices for %d marked "
"classes; the (D0,D1,D1^(1),...,D1^(C)) form is required."
% (ist + 1, 0 if map_c is None else max(0, len(map_c) - 2), len(marked)), error_id='LINE:input:Model')
D1c = [np.asarray(map_c[1 + _mark_index(ist, rr)], dtype=np.float64) for rr in marked]
mods.append({
'ist': ist, 'node': nd, 'arrival': bool(arrival), 'classes': marked,
'D0': np.asarray(map_c[0], dtype=np.float64), 'D1': D1c,
'order': np.asarray(map_c[0]).shape[0],
'is_mmpp': all(_is_diagonal(D1k) for D1k in D1c),
})
for rr in marked:
done[rr] = True
else:
D1 = np.asarray(map_ir[1], dtype=np.float64)
mods.append({
'ist': ist, 'node': nd, 'arrival': bool(arrival), 'classes': [r],
'D0': np.asarray(map_ir[0], dtype=np.float64), 'D1': [D1],
'order': np.asarray(map_ir[0]).shape[0],
'is_mmpp': _is_diagonal(D1),
})
done[r] = True
return mods
__all__ = [
'ChainParams',
'sn_region_members',
'sn_pn_firing_rates',
'sn_pn_avg_rates',
'sn_get_arvr_from_tput',
'sn_get_node_arvr_from_tput',
'sn_get_node_tput_from_tput',
'sn_get_product_form_chain_params',
'sn_set_routing_prob',
'sn_map_modulation',
]