"""
JMT solver integration.
This implementation calls JMT via subprocess (command line), matching how
MATLAB's SolverJMT works. No JVM integration in Python itself.
The solver:
1. Writes the model to JSIM/JMVA XML format
2. Calls JMT via command line
3. Parses the result XML file
Copyright (c) 2012-2026, Imperial College London
All rights reserved.
"""
import numpy as np
import pandas as pd
import os
import re
import sys
import tempfile
import shutil
from typing import Optional, Dict, Any, List, Tuple, Set
from ....api.sn.transforms import sn_get_residt_from_respt
from ....api.io.logging import line_debug, line_warning, line_ack
from ..._solver_base import OptionsDict # noqa: E402 (one definition, shared)
from dataclasses import dataclass, field, fields as _dc_fields, replace
from ....constants import default_verbose
from ....api.solvers.jmt.handler import (
solver_jmt,
SolverJMTOptions as _SolverJMTOptions,
SolverJMTReturn,
is_jmt_available,
_get_jmt_jar_path,
)
from ...base import NetworkSolver, method_label, method_type
from ...config_validate import validate_options
from line_solver.api.io.errors import LineInputError, LineUnsupportedError
#: The refusal for the removed 'replication' method, worded as in MATLAB Solver.m.
JMT_REPLICATION_REMOVED = ("the 'replication' method was removed: use a finite timespan with "
"options.config.replications")
#: Default size of the transient ensemble, options.config.replications (as LDES spells it).
JMT_DEFAULT_REPLICATIONS = 10
def jmt_replications(value):
"""Validate a transient-ensemble size: an integer >= 1 (options.config.replications)."""
try:
f = float(value)
except (TypeError, ValueError):
raise ValueError("SolverJMT: options.config.replications must be an integer >= 1, got %r."
% (value,))
if not np.isfinite(f) or f != int(f) or f < 1:
raise ValueError("SolverJMT: options.config.replications must be an integer >= 1, got %r."
% (value,))
return int(f)
#: The JMVA algorithms that solve a CLOSED product-form network only.
#:
#: RECAL, CoMoM, Chow, Bard-Schweitzer (both spellings), AQL, Linearizer and De
#: Souza-Muntz Linearizer. Measured against JMT 1.2.x: each answers an open or a
#: mixed model with ``jmt.common.exception.UnsupportedModelException: The
#: selected solver cannot handle open classes, please choose another.`` and a
#: load-dependent one with the same exception naming load-dependent stations,
#: while the exact MVA engine behind 'jmva' and 'jmva.mva' serves both.
_JMVA_CLOSED_ONLY = frozenset((
'jmva.amva', 'jmva.recal', 'jmva.comom', 'jmva.chow',
'jmva.bs', 'jmva.aql', 'jmva.lin', 'jmva.dmlin',
))
def jmva_is_closed_only(method):
"""True for a JMVA algorithm restricted to closed single-server networks."""
return str(method or '').lower() in _JMVA_CLOSED_ONLY
def jmt_method_refusal(sn, method, options=None, engine=None):
"""The structural half of SolverJMT's method gate; '' when admissible.
Four rules: immediate feedback, the removed 'replication' method name,
single-server stations for the eight closed-form JMVA algorithms, and a
binding finite buffer, which neither engine carries. None of them has a
registry feature name.
IMMEDIATE FEEDBACK: a job that self-loops keeps its server instead of
re-queueing, which neither JMT document can state. runAnalyzer used to WARN
and return no solution, so the gate called the pair runnable and the table
came back empty. The rule is keyed on the CALLER: without an ``engine``
argument the asker is SolverJMT's own gate or analyzer. The JMVA writer
passes ``'jmva'`` on behalf of the qns methods of SolverLQNS too, which refuses the same feature
in its own words (``qns_immfeed_refusal``), so this sentence must not be
handed to it.
ONE PREDICATE, TWO CALLERS:
``supportsModelMethod`` asks it so findSolver and SolverAUTO never offer a
pair that would die at run time, and the analyzer asks it again so a caller
naming the method by hand gets the same sentence rather than a JMT stack
trace. A second copy of either rule is how the gate and the run drift into
two different answers.
Everything a feature name CAN state lives in ``getMethodFeatureSet``
instead: the JMVA envelope is narrower than the JSIM one, and the
closed-only algorithms additionally drop OpenClass and LoadDependence.
"""
method = str(method or '')
asked_by_jmt = not engine
if asked_by_jmt and sn is not None:
from ....api.sn import sn_has_immfeed
if sn_has_immfeed(sn):
return ('SolverJMT does not support immediate feedback (sn.immfeed): neither the '
'JSIM nor the JMVA document can keep a self-looping job on its server. Use '
'SolverCTMC or SolverSSA, whose state space carries the self-loop.')
if method.lower() == 'replication':
# Removed: the transient ensemble is 'default' on a finite timespan, sized by
# options.config.replications. An error, not an alias, so a caller learns the new spelling.
return JMT_REPLICATION_REMOVED
if asked_by_jmt and sn is not None:
why = jmt_polling_refusal(sn)
if why:
return why
# An EDD or EDF station whose served classes carry no due date: JSIM aborts on such a job.
# JSIM only, since the JMVA feature set withdraws EDD/EDF outright.
if not method.lower().startswith('jmva'):
from ....api.solvers.jmt.handler import jmt_deadline_refusal, jmt_cache_promotion
why = jmt_deadline_refusal(sn)
if why:
return why
# every Cache the JSIM writer would refuse (FIFO above capacity 1, access graphs
# JMT cannot carry, policies it lacks), so the gate and the writer agree
from ....api.sn import NodeType
nodetype = getattr(sn, 'nodetype', None)
for ind, nt in enumerate(nodetype if nodetype is not None else []):
if nt == NodeType.CACHE and getattr(sn.nodeparam[ind], 'nitems', None) is not None:
_, why = jmt_cache_promotion(sn, ind)
if why:
return why
# no JMT region carries An <= b: the constraint was silently dropped (A, B both admitted, no loss)
lincon = getattr(sn, 'regionlincon', None) or []
if any(lc is not None for lc in lincon):
return ('SolverJMT does not support linear admission constraints on a finite capacity region (setLinearConstraints): no JMT region carries them, so the region would be simulated without its A*n <= b limit. Use SolverCTMC, SolverNC or SolverLDES.')
if jmva_is_closed_only(method):
# JMVA implements these seven algorithms for SINGLE-SERVER stations
# only, which is why the JMVA writer refuses the model rather than
# emitting an <ldstation> the algorithm cannot read. A server count is
# not a declared feature, so it cannot ride in the feature set the way
# the load-dependent scaling of the same restriction does.
nservers = getattr(sn, 'nservers', None)
if nservers is not None:
ns = np.asarray(nservers, dtype=float).ravel()
ns = ns[np.isfinite(ns)]
if ns.size and float(np.max(ns)) > 1.0:
return '%s does not support multi-server stations.' % method
# JMVA's multiclass CoMoM answers a different model every run: when its BTF system is singular,
# SolverMultiClosedCoMoM perturbs every demand with an unseeded MersenneTwister, re-solves, and
# reports the result as an exact CoMoM solution. Keyed on the chain count, as in MATLAB
# jmtMethodRefusal, because the singular case cannot be predicted from the demands; a
# single-chain model runs SolverSingleClosedRECAL and stays exact. An open or mixed model is
# left to the OpenClass refusal of getMethodFeatureSet: the closed-only CoMoM never runs on it.
nchains = int(getattr(sn, 'nchains', 0) or 0) if sn is not None else 0
njobs = np.asarray(getattr(sn, 'njobs', [np.inf]) if sn is not None else [np.inf], dtype=float).ravel()
if method.lower() in ('jmva.comom', 'jmt.jmva.comom') and nchains > 1 and np.all(np.isfinite(njobs)):
return ("%s is refused on a model with %d chains: JMT's multiclass CoMoM "
"engine silently re-solves a randomly perturbed model whenever its exact linear "
"system is singular, and reports that as an exact CoMoM solution, so the answer is "
"wrong and different on every run. Use 'jmva.recal', which is exact and "
"reproducible on the same document, 'jmva'/'jmva.mva' for exact MVA, or "
"SolverNC for a LINE-native normalizing constant." % (method, nchains))
why = jmt_jmva_document_refusal(sn, method, engine)
if why:
return why
return jmt_buffer_capacity_refusal(sn, method)
def jmt_jmva_document_refusal(sn, method, engine=None):
"""What the JMVA document cannot state; '' otherwise (MATLAB jmtMethodRefusal, isJmva arm).
The document carries a station type, a mean demand and a visit count per chain and nothing else.
A FCFS station whose classes have different mean service times is not BCMP, and the engine answers
for the PS station with the same demands (3+2 jobs, rates 2 and 0.5: QLen 1.99/1.76 against the
exact 2.42/1.65); that rule is keyed on the jmva* NAME, since SolverLQNS's qnsolver methods read
the same document with their own FCFS approximations. Fork-join and a non-exponential law at a
FCFS/SIRO/LCFS station are keyed on the ENGINE, as in MATLAB.
"""
if sn is None:
return ''
m = str(method or '').lower()
named_jmva = m.startswith('jmva') or m.startswith('jmt.jmva')
if not (named_jmva or str(engine or '').lower() == 'jmva'):
return ''
from ....api.sn import sn_has_fork_join, sn_has_multi_class_heter_fcfs
from ....constants import GlobalConstants, SchedStrategy
if named_jmva and sn_has_multi_class_heter_fcfs(sn):
return ("%s requires a product-form model, but a FCFS station serves its classes with "
"different mean service times, and the JMVA engine would answer for the "
"processor-sharing station with the same demands. Use SolverJMT 'jsim', SolverMVA "
"(whose FCFS correction handles it), SolverCTMC or SolverLDES." % method)
if sn_has_fork_join(sn):
return ("The JMVA document has no fork or join element, so the analytical engine would "
"solve the model without the synchronisation. Use a simulator (SolverJMT 'jsim', "
"SolverLDES) or SolverMVA, which approximates fork-join.")
rates = np.asarray(sn.rates, dtype=float)
scv = np.asarray(sn.scv, dtype=float)
sensitive = {int(SchedStrategy.FCFS), int(SchedStrategy.SIRO), int(SchedStrategy.LCFS)}
for ist in range(int(sn.nstations)):
sc = sn.sched[ist]
if int(sc) not in sensitive:
continue
for r in range(int(sn.nclasses)):
if np.isfinite(rates[ist, r]) and rates[ist, r] > 0 and np.isfinite(scv[ist, r]) \
and abs(scv[ist, r] - 1.0) > GlobalConstants.FineTol:
return ("Station %s serves class %s under %s with a non-exponential law (SCV %.4g). "
"The JMVA document carries the mean demand only, so the analytical engine "
"would answer for the exponential law. Use a simulator (SolverJMT 'jsim', "
"SolverLDES), SolverMVA, whose AMVA corrections read the SCV, or SolverCTMC "
"or SolverSSA." % (sn.nodenames[int(sn.stationToNode[ist])], sn.classnames[r],
SchedStrategy(int(sc)).name, scv[ist, r]))
return ''
def jmt_polling_refusal(sn):
"""EXHAUSTIVE and K-LIMITED polling, which JMT's engine mis-simulates; '' otherwise.
The export is JMT's own (ExhaustivePollingServer/LimitedPollingServer with
the matching polling get strategy, as JMT's XMLWriter writes them), but the
engine reports a queue length that grows without bound next to a sane
response time, and a throughput below the arrival rate: on two classes at
0.3/s, service Exp(1.5) and switchover Exp(1), QLen 2674 and Tput 0.283
against the exact 0.883 and 0.3 (MVA, LDES and SSA agree), and still QLen
704 with a 0.1 switchover. GATED polling is simulated correctly.
"""
from ....constants import PollingType
sched = getattr(sn, 'sched', None)
nodeparam = getattr(sn, 'nodeparam', None)
if sched is None or nodeparam is None:
return ''
for ist in range(int(sn.nstations)):
sc = sched.get(ist) if isinstance(sched, dict) else sched[ist]
if getattr(sc, 'name', '') != 'POLLING':
continue
npm = nodeparam.get(int(sn.stationToNode[ist])) if isinstance(nodeparam, dict) else None
if not npm:
continue
for par in npm.values():
pt = par.get('pollingType') if isinstance(par, dict) else None
pv = getattr(pt, 'value', pt)
if pv is not None and int(pv) in (int(PollingType.EXHAUSTIVE.value), int(PollingType.KLIMITED.value)):
return ('SolverJMT does not support %s polling: JMT\'s polling server for it reports '
'a queue length that grows without bound and a throughput below the arrival '
'rate on a stable station. Use SolverMVA, SolverLDES or SolverSSA; JMT '
'simulates GATED polling correctly.'
% ('EXHAUSTIVE' if int(pv) == int(PollingType.EXHAUSTIVE.value) else 'K-LIMITED'))
return ''
def jmt_buffer_capacity_refusal(sn, method):
"""A binding finite buffer, which NEITHER engine can carry; '' otherwise.
The two engines fail it for opposite reasons, so the binding TEST is shared
and the verdict is not.
What makes a buffer BIND is not that ``sn.cap`` is finite: ``refresh_capacity``
DERIVES a finite cap for every station nobody capped. It is that the cap is
strictly below the population that can REACH the station, which is the JSIM
writer's own test, and an infinite-server station has no buffer at all. Both
are the writer's (``_jmt_reachable_population``), so the gate binds exactly
where the writer binds.
JSIM exports the buffer, but only for the rules JMT can read, and
``_jmt_station_cap_assert`` -- the writer's own predicate, asked here without
letting it raise -- is what decides which. An open loss buffer and a declared
BAS one stay runnable; only the cases JMT would answer unconstrained go.
JMVA is refused OUTRIGHT: ``write_jmva`` emits a station type, a per-chain
service demand and a per-chain visit count and nothing else, so the document
has no capacity element for the buffer to ride in. Measured on a closed
Delay+FCFS model, N=4, cap 2: every jmva method reported 2.19 jobs at a
station that can hold 2, against the exact 1.33. This is the rule SolverMVA,
SolverNC and SolverLQNS already apply -- and SolverLQNS writes THIS SAME
DOCUMENT, so the jmva arm was the one hole in it.
WHICH ENGINE IS ASKED ABOUT matters, and this port derives it from the method
name because nothing else reaches the predicate: unlike MATLAB and the JAR,
SolverLQNS here writes its own JMVA document rather than borrowing
SolverJMT.writeJMVA. A name SolverJMT does not implement therefore gets no
verdict rather than JSIM's.
"""
from ....api.solvers.jmt.handler import (
_jmt_reachable_population, _jmt_station_cap_assert)
name = str(method or '').lower()
is_jmva = name.startswith('jmva')
if not is_jmva and name not in ('default', 'jsim'):
return ''
if getattr(sn, 'isbasdestination', None) is None:
# Derived, and the native struct refresh does not build it: populate it
# exactly as the JSIM writer does, so _jmt_station_cap_assert sees a BAS
# declared on the UPSTREAM station here as it does there (BUG-83).
from ....api.state.ctmc_ssg import _populate_isbasblocking
_populate_isbasblocking(sn)
cap = getattr(sn, 'cap', None)
if cap is None:
return ''
cap = np.asarray(cap, dtype=float).ravel()
nservers = np.asarray(getattr(sn, 'nservers', []), dtype=float).ravel()
from ....api.sn import NodeType
nodetype = getattr(sn, 'nodetype', None)
for ist in range(int(getattr(sn, 'nstations', 0))):
# A SOURCE AND A SINK HAVE NO BUFFER THAT CAN BIND. The Source IS the
# external world and the Sink absorbs, so neither ever holds a job a
# capacity could refuse, yet refresh_capacity writes them a row like any
# other station. Excluded on NODE TYPE, as the shared binding-capacity
# gate of SolverMVA/SolverNC excludes them, and not by name.
if nodetype is not None:
ntype = nodetype[int(sn.stationToNode[ist])]
if ntype in (NodeType.SOURCE, NodeType.SINK):
continue
# UNBOUNDED IS inf HERE. The JAR cannot carry inf on sn.cap -- its
# Station.cap is an int whose "no bound" value is Integer.MAX_VALUE, and
# refreshCapacity SUMS that sentinel across the classes served, so a
# mixed station comes out as 2147483647 + N there and needs
# SaveHandlers.jmtCapIsUnbounded. MATLAB, this port and C++ all default
# station.cap to inf, so isfinite is the whole test.
if ist >= cap.size or not np.isfinite(cap[ist]):
continue
if cap[ist] >= _jmt_reachable_population(sn, ist):
continue
if ist < nservers.size and np.isinf(nservers[ist]):
continue
if is_jmva:
return ("Station %s carries a finite capacity %d that binds. The JMVA document "
"has no capacity element at all, so the analytical engine would solve the model "
"as if the buffer were unbounded and report that as the answer. Use the "
"'jsim' method, which exports the buffer with its drop rule when JMT can "
"express it, or SolverCTMC, SolverSSA or SolverLDES."
% (sn.nodenames[int(sn.stationToNode[ist])], int(cap[ist])))
try:
_jmt_station_cap_assert(sn, ist)
except ValueError as exc:
return str(exc)
return ''
class _instance_or_static:
"""Method callable on an instance (self bound) or on the class (self None).
Lets SolverJMT.writeJMVA(sn, fname, options) work as MATLAB's static-style
writeJMVA and the JAR's static method do, while solver.writeJMVA(fname) keeps
working on an instance.
"""
def __init__(self, func):
self.__func__ = func
self.__doc__ = func.__doc__
self.__name__ = func.__name__
def __get__(self, obj, objtype=None):
import functools
return functools.partial(self.__func__, obj)
def _jmva_writer_options(options) -> Dict[str, Any]:
"""{'method', 'samples'} for write_jmva from a dict, a SolverOptions-like object or None."""
if options is None:
return {}
if isinstance(options, dict):
return {k: options[k] for k in ('method', 'samples') if k in options}
return {k: getattr(options, k) for k in ('method', 'samples') if getattr(options, k, None) is not None}
@dataclass
class SolverJMTOptions:
"""Options for the native JMT solver."""
method: str = 'default'
samples: int = 10000
seed: int = 23000
max_simulated_time: float = float('inf')
# options.config.replications: how many seeded JSIM runs 'default' averages on a finite
# timespan (replication k uses seed + k - 1). No effect on a steady-state run.
replications: int = JMT_DEFAULT_REPLICATIONS
conf_int: float = 0.99
# options.confint: report the half-widths of JMT's intervals next to the means (avg table *CI columns).
# False/0 off, True = 95%, a number in (0, 1) the level; when set it is also the level JMT runs at.
confint: Any = False
max_rel_err: float = 0.03
verbose: bool = field(default_factory=default_verbose)
keep: bool = False # Keep temp files after execution
timeout: float = float('inf')
# Backend selection, see api/solvers/jmt/runner.py. rest_url points at a
# JMT REST server (imperialqore/jmt-rest); container overrides the Docker
# image used when no local JVM exists. Both empty means the local JVM.
rest_url: Optional[str] = None
container: Optional[str] = None
# options.config, kept rather than dropped so it is validated against the schema and visible to
# the caller; 'replications' is ALSO lifted into the field above, which is what the runner reads.
config: Dict[str, Any] = field(default_factory=dict)
# THIS FIELD WAS MISSING, and its absence is why `SolverJMT(model,
# lang='cpp')` looked wired and was not: the kwarg went into **kwargs, was
# never read, and every `options.lang` test in this file saw None. The
# solve itself stays JSIM either way -- what lang='cpp' selects is WHICH
# wrapper drives it, this one or `line-cli -s jmt`, so the two can be
# compared. Same default resolution as every other solver.
lang: str = field(default_factory=lambda: os.environ.get('LINE_SOLVER_LANG', 'python'))
def jmt_structural_refusal_id(sn, reason):
"""'LINE:unsupported:Cache' when REASON is one of the JSIM writer's cache refusals
(jmt_cache_promotion), else 'LINE:unsupported:Structure' (MATLAB/JAR jmtStructuralRefusalId)."""
if not reason or sn is None:
return 'LINE:unsupported:Structure'
from ....api.solvers.jmt.handler import jmt_cache_promotion
from ....api.sn import NodeType
nodetype = getattr(sn, 'nodetype', None)
for ind, nt in enumerate(nodetype if nodetype is not None else []):
if nt == NodeType.CACHE and getattr(sn.nodeparam[ind], 'nitems', None) is not None:
why = jmt_cache_promotion(sn, ind)[1]
if why and why in reason:
return 'LINE:unsupported:Cache'
return 'LINE:unsupported:Structure'
[docs]
class SolverJMT(NetworkSolver):
"""
JMT solver integration.
This solver provides discrete-event simulation and analytical methods
via command line
is launched as an external process, exactly like MATLAB's SolverJMT.
Supported methods:
- 'default': JSIM; on a finite timespan, the mean of options.config.replications
seeded JSIM runs (the transient ensemble)
- 'jsim': one JSIM run (discrete-event simulation)
- 'jmva' / 'jmva.mva': Mean Value Analysis
- 'jmva.amva': Approximate MVA
- 'jmva.recal': RECALsimulation
- 'jmva.comom': CoMoM algorithm
- 'jmva.chow': Chow algorithm
- 'jmva.bs': Bard-Schweitzer
- 'jmva.aql': AQL algorithm
- 'jmva.lin': Linearizer
- 'jmva.dmlin': De Souza-Muntz Linearizer
Args:
model: Network model (Python wrapper or native structure)
method: Solution method (default: 'default')
**kwargs: Additional solver options (samples, seed, etc.)
Example:
>>> solver = SolverJMT(model, samples=10000, seed=42)
>>> solver.runAnalyzer()
>>> table = solver.getAvgTable()
"""
def __init__(self, model, method_or_options=None, **kwargs):
self.model = model
# see _kb/09-ldes-and-cache.md (Warm start) for initFromSolver contract
init_solver = None
if method_or_options is not None and hasattr(method_or_options, 'getAvgQLen'):
init_solver = method_or_options
method_or_options = None
# Handle options passed as second argument (MATLAB-style)
if method_or_options is None:
# honor method= keyword (consistent with SolverNC/SolverMVA)
self.method = str(kwargs.pop('method', 'default')).lower()
elif isinstance(method_or_options, str):
self.method = method_or_options.lower()
elif hasattr(method_or_options, 'get'):
# Dict-like options object
self.method = method_or_options.get('method', 'default')
for key in ('config', 'replications', 'timespan', 'max_simulated_time'):
if key in method_or_options:
kwargs.setdefault(key, method_or_options[key])
if 'samples' in method_or_options:
kwargs.setdefault('samples', method_or_options['samples'])
if 'seed' in method_or_options:
kwargs.setdefault('seed', method_or_options['seed'])
if 'keep' in method_or_options:
kwargs.setdefault('keep', method_or_options['keep'])
if hasattr(method_or_options, 'verbose'):
kwargs.setdefault('verbose', method_or_options.verbose)
elif isinstance(method_or_options, SolverJMTOptions):
# A whole native options object (e.g. the logged copy of sampleSysAggr): every field crosses.
self.method = str(method_or_options.method or 'default').lower()
for f in _dc_fields(SolverJMTOptions):
if f.name != 'method':
kwargs.setdefault(f.name, getattr(method_or_options, f.name))
elif hasattr(method_or_options, 'method'):
# SolverOptions-like object
self.method = getattr(method_or_options, 'method', 'default')
for key in ('config', 'replications', 'timespan', 'max_simulated_time'):
if getattr(method_or_options, key, None) is not None:
kwargs.setdefault(key, getattr(method_or_options, key))
if hasattr(method_or_options, 'samples'):
kwargs.setdefault('samples', method_or_options.samples)
if hasattr(method_or_options, 'seed'):
kwargs.setdefault('seed', method_or_options.seed)
if hasattr(method_or_options, 'verbose'):
kwargs.setdefault('verbose', method_or_options.verbose)
else:
self.method = 'default'
self.method = str(self.method or 'default').lower()
if self.method == 'replication':
raise ValueError(JMT_REPLICATION_REMOVED)
# An unknown keyword is an error, as for MVA/CTMC/SSA: it used to be
# dropped silently, so a misspelt option ran as if never given.
import dataclasses as _dc
_known = {'samples', 'seed', 'verbose', 'keep', 'conf_int', 'confint', 'max_rel_err',
'timespan', 'max_simulated_time', 'replications', 'timeout', 'rest_url',
'container', 'lang', 'config', 'method'}
_known |= {f.name for f in _dc.fields(SolverJMTOptions)}
from ..._solver_base import generic_option_keys
_known |= generic_option_keys()
for _key in kwargs:
if _key not in _known:
raise TypeError("SolverJMT got an unexpected keyword argument '%s'" % _key)
# Parse options
samples = kwargs.get('samples', 10000)
seed = kwargs.get('seed', 23000)
verbose = kwargs.get('verbose', default_verbose())
keep = kwargs.get('keep', False)
from ...base import parse_confint
confint = kwargs.get('confint', False)
ci_on, ci_level = parse_confint(confint)
conf_int = kwargs.get('conf_int', ci_level if ci_on else 0.99)
max_rel_err = kwargs.get('max_rel_err', 0.03)
timespan = kwargs.get('timespan')
max_simulated_time = kwargs.get('max_simulated_time',
float(timespan[1]) if isinstance(timespan, (list, tuple, np.ndarray))
and len(timespan) >= 2 else float('inf'))
replications = self._resolve_replications(kwargs)
timeout = kwargs.get('timeout', float('inf'))
rest_url = kwargs.get('rest_url', None)
container = kwargs.get('container', None)
lang = kwargs.get('lang', os.environ.get('LINE_SOLVER_LANG', 'python'))
self.options = SolverJMTOptions(
lang=lang,
method=self.method,
samples=samples,
seed=seed,
max_simulated_time=max_simulated_time,
replications=replications,
conf_int=conf_int,
confint=confint,
max_rel_err=max_rel_err,
verbose=verbose,
keep=keep,
timeout=timeout,
rest_url=rest_url,
container=container,
config=dict(kwargs['config']) if isinstance(kwargs.get('config'), dict) else (kwargs.get('config') or {})
)
validate_options(self.options, 'JMT')
# JSIM event budget written as <sim maxEvents>; -1 is unlimited (the sample* methods set it)
self.maxEvents = -1
self._result: Optional[SolverJMTReturn] = None
self._sn = None
# Extract network structure
self._extract_network_params()
if init_solver is not None:
self.initFromSolver(init_solver)
@staticmethod
def _resolve_replications(kwargs):
"""The flat ``replications`` when set, else options.config.replications, else 10.
``config`` may be a dict or an attribute object; the order is MATLAB's jmtReplications,
which reads the knob as SolverLDES does.
"""
value = kwargs.get('replications')
if value is None:
cfg = kwargs.get('config')
if isinstance(cfg, dict):
value = cfg.get('replications')
elif cfg is not None:
value = getattr(cfg, 'replications', None)
if value is None:
return JMT_DEFAULT_REPLICATIONS
return jmt_replications(value)
def getName(self) -> str:
"""Get the name of this solver."""
return "JMT"
get_name = getName
def _extract_network_params(self):
"""Extract parameters from the model."""
model = self.model
# see _kb/06-solver-catalog.md (Wrappers: "JMT python wrapper: export/import internals")
if hasattr(model, 'refresh_struct'):
model.refresh_struct()
if hasattr(model, '_sn') and model._sn is not None:
self._sn = model._sn
return
# Fallback: Use existing _sn if refresh_struct is not available
if hasattr(model, '_sn') and model._sn is not None:
self._sn = model._sn
return
# Native model (snake-case get_struct()); no JAR-wrapper bridge here.
if hasattr(model, 'get_struct'):
self._sn = model.get_struct()
if self._sn is not None:
return
# Already a native NetworkStruct
if hasattr(model, 'nclasses') and hasattr(model, 'nstations'):
self._sn = model
return
raise ValueError(
"Cannot extract a native NetworkStruct from model. Native solvers "
"accept only native Network / NetworkStruct inputs (no JAR wrapper).")
[docs]
def supportsTransientAnalysis(self):
"""Transient averages are available (simulation restricted to options.timespan)."""
return True
supports_transient_analysis = supportsTransientAnalysis
[docs]
def runAnalyzer(self) -> 'SolverJMT':
"""
Run the JMT analyzer.
Calls JMT via command line and stores the results.
Returns:
self for method chaining
"""
line_ack('JMT', self.options.verbose)
line_debug("JMT analyzer starting: lang=python, samples=%d, seed=%d",
self.options.samples, self.options.seed, options=self.options)
if self._sn is None:
raise RuntimeError("Network structure not available")
# see _kb/06-solver-catalog.md (Wrappers: "JMT python wrapper: export/import internals")
model = getattr(self, 'model', None)
if model is not None and hasattr(model, 'get_used_lang_features'):
self.runAnalyzerChecks(self.options)
# The structural half of the gate, asked again here so a caller who
# reaches the analyzer with the checks disabled still gets the gate's
# own sentence rather than a JMT stack trace.
structural = jmt_method_refusal(self._sn, self.options.method, self.options)
if structural:
raise LineUnsupportedError(structural, error_id=jmt_structural_refusal_id(self._sn, structural))
method = self.options.method
if method == 'default' and np.isfinite(float(self.options.max_simulated_time)):
# THE TRANSIENT ENSEMBLE: a single sample path is not the transient mean E[N](t), there
# being no time-ergodicity at a fixed t, so options.config.replications seeded JSIM runs
# are averaged onto a common grid. An explicit 'jsim' stays one run, which is what the
# logged copies of sampleSysAggr rely on to not re-enter this arm.
return self._runTransientReplications()
if method in ('jsim', 'default'):
line_debug("JMT: using JSIM method (discrete-event simulation), samples=%d, seed=%d",
self.options.samples, self.options.seed, options=self.options)
elif method.startswith('jmva'):
line_debug("JMT: using JMVA method: %s", method, options=self.options)
else:
line_debug("JMT: using method: %s", method, options=self.options)
samples = self.options.samples
if samples < 5000:
# JSIM's dynamic analyzer refuses a smaller budget with a raw IllegalArgumentException; it used to be
# passed through unchanged after a debug message claiming otherwise. MATLAB runAnalyzer.m raises it.
line_warning('SolverJMT', 'JMT requires at least 5000 samples for each metric, the current value is %d. '
'Starting the simulation with 5000 samples.' % samples)
samples = 5000
# Convert options to handler format
handler_options = _SolverJMTOptions(
method=self.options.method,
samples=samples,
seed=self.options.seed,
max_simulated_time=self.options.max_simulated_time,
conf_int=self.options.conf_int,
max_rel_err=self.options.max_rel_err,
verbose=self.options.verbose,
keep=self.options.keep,
timeout=getattr(self.options, 'timeout', float('inf')),
rest_url=getattr(self.options, 'rest_url', None),
container=getattr(self.options, 'container', None),
max_events=self.maxEvents
)
# Call the handler (pass model for FCR region support)
self._result = solver_jmt(self._sn, handler_options, self.model)
self._apply_cache_hit_probs()
# Print completion message (matches MATLAB verbose guard)
if self.options.verbose:
py_version = f"{sys.version_info.major}.{sys.version_info.minor}.{sys.version_info.micro}"
runtime = self._result.runtime if hasattr(self._result, 'runtime') else 0.0
method = self._result.method if hasattr(self._result, 'method') else self.options.method
from line_solver.solvers.base import print_solver_banner
print_solver_banner(f"JMT analysis [method: {method_label(self.options.method, method)}; type: {method_type('JMT', method_label(self.options.method, method))}; lang: python; env: {py_version}] completed in {runtime:.6f}s.")
return self
def _apply_cache_hit_probs(self):
"""Store JMT's simulated cache hit/miss split on each Cache and re-route at it.
MATLAB runAnalyzer.m sets setResultHitProb/setResultMissProb and the
struct's actualhitprob/actualmissprob from the Cache Hit Rate measure;
the read class is then re-routed in place at that split and the visits
refreshed, so the node and cache tables read the simulated split rather
than the nominal one the routing was linked with.
"""
hits = getattr(self._result, 'cacheHitProb', None)
if not hits:
return
from line_solver.api.sn.transforms import sn_refresh_cacheqn_visits
nodes = self.model.get_nodes() if self.model is not None else []
for ind, hp in hits.items():
hp = np.asarray(hp, dtype=float).copy()
mp = 1.0 - hp
param = self._sn.nodeparam[ind]
param.actualhitprob = hp
param.actualmissprob = mp
if ind < len(nodes) and hasattr(nodes[ind], 'set_result_hit_prob'):
nodes[ind].set_result_hit_prob(hp)
nodes[ind].set_result_miss_prob(mp)
sn_refresh_cacheqn_visits(self._sn)
[docs]
def getAvgTable(self) -> pd.DataFrame:
"""
Get average performance metrics as a DataFrame.
Returns:
DataFrame with columns: Station, Class, QLen, Util, RespT, Tput, ArvR
"""
tran = self._tran_avg_table() # a finite timespan reports the averages at its end, as MATLAB
if tran is not None:
return tran
if self._result is None:
self._ensureAvgResults()
M = self._sn.nstations
K = self._sn.nclasses
nodenames = self._sn.nodenames if self._sn.nodenames else [f'Station{i}' for i in range(M)]
classnames = self._sn.classnames if self._sn.classnames else [f'Class{r}' for r in range(K)]
# Compute ResidT using proper visit ratios from network structure
# This uses the correct formula: WN[ist,k] = RN[ist,k] * V[ist,k] / V[refstat,refclass]
if self._sn is not None and self._sn.visits and self._result.R is not None:
WN = sn_get_residt_from_respt(self._sn, self._result.R, None)
else:
# Fallback: ResidT = RespT (no visit information available)
WN = self._result.R.copy() if self._result.R is not None else np.zeros((M, K))
# Get station names and identify source stations
station_names = []
source_stations = set()
nodetype = self._sn.nodetype if hasattr(self._sn, 'nodetype') else None
for i in range(M):
node_idx = int(self._sn.stationToNode[i]) if self._sn.stationToNode is not None else i
if node_idx < len(nodenames):
station_names.append(nodenames[node_idx])
else:
station_names.append(f'Station{i}')
# Check if source station
if nodetype is not None and node_idx < len(nodetype):
if int(nodetype[node_idx]) == 0: # SOURCE = 0
source_stations.add(i)
# Get arrival rates from rates matrix for source stations
rates = np.asarray(self._sn.rates) if hasattr(self._sn, 'rates') and self._sn.rates is not None else None
rows = []
for i in range(M):
for r in range(K):
is_source = i in source_stations
# Get values with NaN handling
qlen = self._result.Q[i, r] if self._result.Q is not None else np.nan
util = self._result.U[i, r] if self._result.U is not None else np.nan
respt = self._result.R[i, r] if self._result.R is not None else np.nan
residt = WN[i, r] if i < WN.shape[0] and r < WN.shape[1] else respt
arvr = self._result.A[i, r] if self._result.A is not None else np.nan
tput = self._result.T[i, r] if self._result.T is not None else np.nan
# For source stations, replace NaN with 0 and set Tput to arrival rate
if is_source:
qlen = 0.0 if np.isnan(qlen) else qlen
util = 0.0 if np.isnan(util) else util
respt = 0.0 if np.isnan(respt) else respt
residt = 0.0 if np.isnan(residt) else residt
arvr = 0.0 # Source has no arrivals to itself
# Set Tput from arrival rate
if np.isnan(tput) and rates is not None:
stationToNode = np.asarray(self._sn.stationToNode).flatten()
node_idx = int(stationToNode[i])
if node_idx < rates.shape[0] and r < rates.shape[1]:
tput = rates[node_idx, r]
# Filter out rows where all metrics are zero or NaN (matching MATLAB behavior)
# Only include row if at least one metric is non-zero and not NaN
metrics = [qlen, util, respt, residt, arvr, tput]
has_significant_value = any(
(not np.isnan(v) and v > 0) for v in metrics
)
if not has_significant_value:
continue
rows.append({
'Station': station_names[i],
'JobClass': classnames[r],
'QLen': qlen,
'Util': util,
'RespT': respt,
'ResidT': residt,
'ArvR': arvr,
'Tput': tput,
})
df = pd.DataFrame(rows)
from ...base import parse_confint, append_ci_columns
if parse_confint(getattr(self.options, 'confint', False))[0]:
res = self._result
ci = {'Q': res.QNCI, 'U': res.UNCI, 'R': res.RNCI, 'T': res.TNCI, 'A': res.ANCI}
# ResidT is RespT scaled by the visit ratio, so its half-width scales the same way
if res.RNCI is not None and res.R is not None:
with np.errstate(divide='ignore', invalid='ignore'):
ci['W'] = np.where(np.asarray(res.R) > 0, np.asarray(res.RNCI) * np.asarray(WN) / np.asarray(res.R),
np.nan)
df = append_ci_columns(df, station_names, list(classnames), ci)
if not getattr(self, '_table_silent', False):
print(df.to_string(index=False))
return df
[docs]
def getAvgQLen(self) -> np.ndarray:
"""Get average queue lengths (M x K matrix)."""
if self._result is None:
self._ensureAvgResults()
return self._result.Q if self._result.Q is not None else np.array([])
[docs]
def getAvgUtil(self) -> np.ndarray:
"""Get average utilizations (M x K matrix)."""
if self._result is None:
self._ensureAvgResults()
return self._result.U if self._result.U is not None else np.array([])
[docs]
def getAvgRespT(self) -> np.ndarray:
"""Get average response times (M x K matrix)."""
if self._result is None:
self._ensureAvgResults()
return self._result.R if self._result.R is not None else np.array([])
[docs]
def getAvgTput(self) -> np.ndarray:
"""Get average throughputs (M x K matrix)."""
if self._result is None:
self._ensureAvgResults()
return self._result.T if self._result.T is not None else np.array([])
[docs]
def getAvgArvR(self) -> np.ndarray:
"""Get average arrival rates (M x K matrix)."""
if self._result is None:
self._ensureAvgResults()
return self._result.A if self._result.A is not None else np.array([])
[docs]
def getAvgFcr(self) -> np.ndarray:
"""The finite-capacity-region rows, ((nregions*6) x nclasses) or empty.
WHY A SEPARATE ACCESSOR. `getAvg` returns the STATION metrics alone, and
a region is not a station: its rows live past the last one, which is
where `getAvgNodeTable` reads them from to fill the FCR pseudo-node. A
host bridging through `getAvg` (MATLAB `lang='python'`) therefore saw no
region at all and dropped the FCR row from its node table
(fcr_mm1waitq[M2P], "row FCR1 missing").
The six blocks are stacked in the order Q, U, R, W, A, T, each
(nregions x nclasses), so one marshalled matrix carries all of them.
Util and ArvR are NaN: JMT reports neither for a region.
"""
if self._result is None:
self._ensureAvgResults()
F = int(getattr(self._sn, 'nregions', 0) or 0)
Qfcr = getattr(self._result, 'Qfcr', None)
if F <= 0 or Qfcr is None:
return np.array([])
K = int(self._sn.nclasses)
nan = np.full((F, K), np.nan)
zero = np.zeros((F, K))
def blk(name, default):
v = getattr(self._result, name, None)
return default if v is None else np.asarray(v, dtype=float).reshape(F, K)
return np.vstack([blk('Qfcr', zero), nan, blk('Rfcr', zero),
blk('Wfcr', zero), nan, blk('Tfcr', zero)])
[docs]
def sampleSysAggr(self, num_events: Optional[int] = None):
"""Sample system-wide aggregated state trajectories via JMT logging.
Faithful port of MATLAB SolverJMT.sampleSysAggr: all non-Source stations
are logged in a temporary model copy, simulated, and their per-class
queue-length trajectories reconstructed and interpolated (previous/step)
onto a common timeline (the union of all stations' event times).
Args:
num_events: Number of events to sample (default: uses options.samples)
Returns:
SamplePath (fields also readable as keys) with 't' (common timeline), 'state' (list of per-station
len(t) x nclasses matrices; Source stations = [[inf]]), 'handle',
'event' (chronological event list) and 'isaggregate'=True.
"""
import os
from ....api.solvers.jmt.handler import parse_tran_state
from ....lang.base import NodeType
sn = self._sn if self._sn is not None else self.model.get_struct()
if num_events is None:
num_events = getattr(self.options, 'samples', 10000)
K = sn.nclasses
M = sn.nstations
def _is_source(ind):
nt = sn.nodetype[ind]
nt_val = nt.value if hasattr(nt, 'value') else int(nt)
src_val = NodeType.SOURCE.value if hasattr(NodeType.SOURCE, 'value') else int(NodeType.SOURCE)
return nt_val == src_val
nnodes = self.model.get_number_of_nodes()
is_node_logged = [False] * nnodes
station_node = [int(sn.stationToNode[ist]) for ist in range(M)]
for ist in range(M):
ind = station_node[ist]
if not _is_source(ind):
is_node_logged[ind] = True
# the initial state counts as an event (MATLAB sampleSysAggr: numEvents - 1)
model_copy, log_path = self._run_logged_copy(is_node_logged, num_events - 1)
class_names = [sn.classnames[r] for r in range(K)] if sn.classnames is not None else None
stat_t = [None] * M
stat_q = [None] * M
all_events = []
for ist in range(M):
ind = station_node[ist]
if not is_node_logged[ind]:
continue
name = sn.nodenames[ind]
preload = self._node_preload(ind, K)
arv = os.path.join(log_path, f"{name}-Arv.csv")
dep = os.path.join(log_path, f"{name}-Dep.csv")
state, evtype, evclass, evjob = parse_tran_state(arv, dep, preload, class_names)
_, uniq = np.unique(state[:, 0], return_index=True)
uniq = np.sort(uniq)
stat_t[ist] = state[uniq, 0]
stat_q[ist] = state[uniq, 1:1 + K]
for e in range(len(evtype)):
if not np.isnan(evjob[e]):
all_events.append({
'event': evtype[e], 'node': ind,
'class': int(evclass[e]) if not np.isnan(evclass[e]) else None,
't': float(state[e, 0]),
'job': int(evjob[e]) if not np.isnan(evjob[e]) else None,
})
# Common timeline = union of all stations' timestamps, capped at the
# earliest station end so previous-interpolation never extrapolates.
t_union = np.array([])
maxes = []
for ist in range(M):
if stat_t[ist] is not None and len(stat_t[ist]) > 0:
t_union = np.union1d(t_union, stat_t[ist])
maxes.append(np.max(stat_t[ist]))
if maxes:
t_union = t_union[t_union <= min(maxes)]
def _prev_interp(ts, ys, tq):
if len(ts) == 0:
return np.full(len(tq), np.nan)
idx = np.searchsorted(ts, tq, side='right') - 1
idx = np.clip(idx, 0, len(ys) - 1)
return ys[idx]
state_list = []
for ist in range(M):
ind = station_node[ist]
if _is_source(ind) or stat_t[ist] is None:
state_list.append(np.array([[np.inf]]))
continue
cols = np.zeros((len(t_union), K))
for r in range(K):
cols[:, r] = _prev_interp(stat_t[ist], stat_q[ist][:, r], t_union)
state_list.append(cols)
all_events.sort(key=lambda ev: ev['t'])
handles = [self.model.get_stations()[ist] if hasattr(self.model, 'get_stations')
else ist for ist in range(M)]
self._cleanup_log_dir(log_path)
from line_solver.solvers.solver_ssa.solver_ssa import SamplePath
return SamplePath(handle=handles, t=t_union, state=state_list, event=all_events, isaggregate=True)
[docs]
def getProbSysAggr(self) -> float:
"""Get system state probability via simulation sampling.
Uses JMT simulation with logging to estimate the probability of the
system being in the current aggregated state (as set via setState).
The probability is computed as the fraction of time the system spends
in the target state during simulation.
Note: This method requires Logger nodes to be present in the model for
full functionality. If log files are not available, it attempts to
estimate probabilities from the simulation's average metrics.
Returns:
float: Estimated probability of the current system state.
Returns 0.0 if the state was not observed during simulation.
"""
if getattr(self.options, 'lang', 'python') == 'cpp':
import warnings
from ...cpp_dispatch import jmt_prob_aggr_via_cpp
r = jmt_prob_aggr_via_cpp(self)
if not r['sysStateSeen']:
warnings.warn("the system state was not seen in the simulation, "
"so its probability is reported as 0")
return r['probSysAggr']
if self._sn is None:
self._extract_network_params()
sn = self._sn
# Try to get state samples from simulation with logging
sample_result = self.sampleSysAggr()
if sample_result is not None and 't' in sample_result and 'state' in sample_result:
timestamps = sample_result['t']
states = sample_result['state']
if len(timestamps) >= 2:
# Get target state from the model's current state
nstations = sn.nstations
nclasses = sn.nclasses
target_state = np.zeros((nstations, nclasses))
# Get state from sn.state dict
if hasattr(sn, 'state') and sn.state is not None:
for ist in range(nstations):
node_idx = int(sn.stationToNode[ist]) if sn.stationToNode is not None else ist
stateful_idx = int(sn.nodeToStateful[node_idx]) if sn.nodeToStateful is not None and node_idx < len(sn.nodeToStateful) else -1
if stateful_idx >= 0 and isinstance(sn.state, dict):
for key, state_vec in sn.state.items():
key_idx = -1
if hasattr(key, 'statefulIndex'):
key_idx = key.statefulIndex
elif hasattr(key, 'getStatefulIndex'):
key_idx = key.getStatefulIndex()
if key_idx == stateful_idx:
if state_vec is not None:
state_arr = np.asarray(state_vec).flatten()
for r in range(min(nclasses, len(state_arr))):
target_state[ist, r] = state_arr[r]
break
# see _kb/06-solver-catalog.md (Wrappers: "JMT python wrapper:
# export/import internals") for the vectorized snapshot matching
time_diffs = np.diff(timestamps)
total_time = np.sum(time_diffs)
if total_time > 0:
n_steps = len(time_diffs)
tol = 0.5 + 1e-5 * np.abs(target_state)
match = np.ones(n_steps, dtype=bool)
for ist in range(nstations):
si = states[ist] if ist < len(states) else None
row = np.zeros((n_steps, nclasses))
if si is not None and np.all(np.isfinite(si)):
# Source station ([inf]) keeps the zero row.
rows_avail = min(n_steps, si.shape[0])
ncol = min(nclasses, si.shape[1])
row[:rows_avail, :ncol] = si[:rows_avail, :ncol]
match &= np.all(
np.abs(row - target_state[ist]) <= tol[ist], axis=1)
matching_time = float(np.sum(time_diffs[match]))
if matching_time > 0:
return matching_time / total_time
return 0.0
[docs]
def getProbAggr(self, station, state_a=None) -> float:
"""Probability of a per-class job distribution at a station, estimated by simulation.
Port of MATLAB @SolverJMT/getProbAggr.m: the fraction of simulated time the
station's logged per-class counts (sampleAggr) equal `state_a`, which defaults
to the station's declared state. Under lang='cpp' this is `-s jmt -a prob`,
the same dwell-weighted estimate from one logged JSIM run.
Args:
station: station object, or its 0-based station index
state_a: per-class job counts, e.g. [2, 1]; default the declared state
"""
if getattr(self.options, 'lang', 'python') == 'cpp':
import warnings
from ...cpp_dispatch import jmt_prob_aggr_via_cpp
r = jmt_prob_aggr_via_cpp(self)
ist = int(station)
if not (0 <= ist < r['probAggr'].size):
raise ValueError("station index %r is outside 0..%d"
% (station, r['probAggr'].size - 1))
if ist < len(r['stateSeen']) and not r['stateSeen'][ist]:
warnings.warn("station %d's state was not seen in the simulation, "
"so its probability is reported as 0" % ist)
return float(r['probAggr'][ist])
from ....api.state.marginal import toMarginal
sn = self.model.getStruct()
ist = int(station) if isinstance(station, (int, np.integer)) else station.get_station_index0()
ind = int(np.asarray(sn.stationToNode).flatten()[ist])
if state_a is None:
isf = int(np.asarray(sn.nodeToStateful).flatten()[ind])
_, nir, _, _ = toMarginal(sn, ind, np.asarray(sn.state[isf], dtype=float).flatten())
state_a = np.asarray(nir, dtype=float).reshape(-1)[:int(sn.nclasses)]
path = self.sampleAggr(self.model.get_nodes()[ind])
states = np.atleast_2d(np.asarray(path['state'], dtype=float))
t = np.asarray(path['t'], dtype=float).reshape(-1)
dt = np.append(np.diff(t), 0.0)
rows = np.all(np.isclose(states, np.asarray(state_a, dtype=float).reshape(1, -1)), axis=1)
return float(dt[rows].sum() / dt.sum()) if dt.sum() > 0 else 0.0
[docs]
def getRuntime(self) -> float:
"""Get solver runtime in seconds."""
if self._result is None:
return 0.0
return self._result.runtime
[docs]
def getMethod(self) -> str:
"""Get the method used."""
if self._result is None:
return self.method
return self._result.method
[docs]
def listValidMethods(self) -> List[str]:
"""List valid methods for this solver."""
# the names and their order are this solver's row of the method registry
# (doc/line-methods.json, line_solver.solvers.method_registry)
from line_solver.solvers import method_registry
return method_registry.methods('JMT')
[docs]
def isStochasticMethod(self, method):
"""Simulation-based methods (default, jsim) return
stochastic estimates. The analytical JMVA methods do not, except
for the sampling-based variants (e.g. jmva.ls).
"""
if not method:
return True # default resolves to jsim simulation
method_names = re.split(r'[./]', str(method).lower())
if 'jmva' in method_names:
return any(tok in ('ls', 'mci', 'imci', 'sampling') for tok in method_names)
return True
is_stochastic_method = isStochasticMethod
[docs]
@staticmethod
def isAvailable() -> bool:
"""Check if JMT solver is available."""
return is_jmt_available()
[docs]
@staticmethod
def getFeatureSet() -> Set[str]:
"""Get the set of features supported by this solver."""
return {
'Sink', 'Source', 'Router', 'ClassSwitch',
'Delay', 'DelayStation', 'Queue',
'Fork', 'Join', 'Forker', 'Joiner', 'Logger',
'JoinPartial', # quorum join, written out as a jmt PartialJoin
# A variable forking level: saveForkStrategy turns isSimplifiedFork
# off and writes the per-branch entries, so jmt reads the counts,
# the probabilities and the degree distribution rather than sending
# one job down every link.
'ForkFanoutVector', 'ForkFanoutRandom', 'ForkBranchProbability',
'Coxian', 'Cox2', 'APH', 'Erlang', 'Exp', 'HyperExp',
'Det', 'Gamma', 'Lognormal', 'MAP', 'MMPP2',
'Normal', 'PH', 'Pareto', 'Weibull', 'Replayer', 'Uniform',
'StatelessClassSwitcher', 'InfiniteServer', 'SharedServer',
'Buffer', 'Dispatcher', 'Server', 'JobSink', 'RandomSource',
'ServiceTunnel', 'LogTunnel', 'Linkage',
'Enabling', 'Inhibiting', 'Timing', 'Firing', 'Storage', 'Place', 'Transition',
'SchedStrategy_INF', 'SchedStrategy_PS', 'SchedStrategy_DPS',
'SchedStrategy_FCFS', 'SchedStrategy_GPS', 'SchedStrategy_SIRO',
'SchedStrategy_HOL', 'SchedStrategy_PSPRIO', 'SchedStrategy_DPSPRIO',
'SchedStrategy_GPSPRIO', 'SchedStrategy_LCFS', 'SchedStrategy_LCFSPR',
'SchedStrategy_LCFSPRIO', 'SchedStrategy_LCFSPRPRIO',
# LCFSPI is emitted by the writer (QueuePutStrategies.LCFSPIStrategy).
# The FCFS preemptive family and the priority variants are emitted
# since the writer gained MATLAB's FCFSPR/FCFSPI(Priority) and
# LCFSPIPriority branches and the PreemptiveServer promotion.
'SchedStrategy_LCFSPI', 'SchedStrategy_LCFSPIPRIO',
'SchedStrategy_FCFSPR', 'SchedStrategy_FCFSPI',
'SchedStrategy_FCFSPRPRIO', 'SchedStrategy_FCFSPIPRIO', 'SchedStrategy_FCFSPRIO',
'SchedStrategy_SEPT', 'SchedStrategy_SRPT', 'SchedStrategy_SRPTPRIO', 'SchedStrategy_LEPT',
'SchedStrategy_SJF', 'SchedStrategy_LJF', 'SchedStrategy_LPS',
'SchedStrategy_POLLING', 'SchedStrategy_EXT',
# EDD and EDF are REAL in JMT 1.2.x (EDDStrategy/EDFStrategy order the buffer by the due
# dates the writer emits as <classSoftDeadlines>); a served class without a deadline aborts
# JSIM, so jmt_deadline_refusal gates it (MATLAB SolverJMT.getFeatureSet).
'SchedStrategy_EDD', 'SchedStrategy_EDF',
'RoutingStrategy_PROB', 'RoutingStrategy_RAND',
'RoutingStrategy_RROBIN', 'RoutingStrategy_WRROBIN',
'RoutingStrategy_JSQ',
# RoutingStrategy_SQ not declared: JMT's PowerOfK draws d DISTINCT candidates (JSQ once d >= #outputs), LINE's SQ(d) draws with replacement.
'ClosedClass', 'SelfLoopingClass', 'OpenClass',
# Exported by _write_cache_node (port of MATLAB saveCacheStrategy.m);
# JMT's Cache Hit Rate measure gives the simulated hit/miss split.
# FIFO/SFIFO export only where JMT's in-list reorder is a no-op (the
# writer refuses the rest by name); CLIMB reaches the writer already
# rewritten as FIFO on unit lists, which is that case.
'Cache', 'CacheClassSwitcher',
'ReplacementStrategy_RR', 'ReplacementStrategy_FIFO',
'ReplacementStrategy_SFIFO', 'ReplacementStrategy_LRU',
'ReplacementStrategy_CLIMB',
'Region',
# Limited load dependence reaches JMT only as a SERVER COUNT: the
# JSIM writer exports max(nservers, max(alpha)) and the JMVA writer
# the matching <ldstation>. That is exact for alpha(n) = min(n,c)
# and for nothing else, so supportsModelMethod refuses any other
# scaling by name.
'LoadDependence',
# Exported as delayOffTime/setUpTime (_write_delayoff_strategy).
'SetupDelayOff',
# Exported as classParallelism (_write_class_parallelism).
'ServerParallelism',
# Heterogeneous server pools: the type names, the servers per type
# and the compatibility matrix are exported as serverTypesNames /
# serverTypesNumOfServers / serverTypesCompatibilities, so jsim
# simulates the pools rather than a station of the same total size.
'HeteroServers',
# Exported as Impatience/Reneging and Impatience/Balking strategies.
'Reneging', 'Balking',
# c-server stations (the writer emits numberOfServers) and finite
# buffers with their drop rule (a capacity plus the dropStrategy
# text): jsim exports both, and the structural capacity gate keeps
# refusing the buffers JMT would answer unconstrained (closed WAITQ,
# BBS, RSRD, retrial-with-limit); the JMVA set withdraws the buffer.
# 'Retrial' is NOT declared, and this port is the only one that
# withholds it: MATLAB, the JAR and C++ pick the retrial Queue
# constructor and write the per-class orbit delay, while this writer
# has no such branch, so the orbit would simply not be exported.
'MultiServer', 'FiniteCapacity',
}
[docs]
@staticmethod
def getJMVAFeatureSet() -> Set[str]:
"""What the JMVA ANALYTICAL engine accepts, much less than JSIM.
Derived from the writer rather than guessed: ``write_jmva`` emits, per
station, a ``<delaystation>``, a ``<listation>`` or an ``<ldstation>``,
a per-chain ``<servicetime>`` and a per-chain ``<visit>``, and at model
level the closed populations, the open arrival rates and the reference
station. NOTHING ELSE IN THE MODEL REACHES JMVA, so a construct whose
whole effect is not carried by (station type, demand, visits,
population) would be solved away silently.
Dropped from the JSIM set, and why:
* Fork/Join and the fan-out names -- no fork element exists, and a
visit ratio cannot express the join synchronization.
* Place/Transition and the Petri-net sections -- no counterpart.
* Region -- JMVA has no finite capacity region.
* Reneging/Balking -- no impatience element; the abandonment would
simply not happen.
* SetupDelayOff, ServerParallelism, HeteroServers -- each a
server-side attribute the JMVA document has no slot for.
* the non-BCMP disciplines -- the writer emits NO discipline at all,
so a priority, weighted, size-based or limited-sharing station
would be solved as an ordinary load-independent one. Only the four
BCMP station types survive the encoding, the same line SolverNC and
SolverMVA draw.
* the state-dependent routings (RROBIN, WRROBIN, JSQ, SQ) -- the
document carries mean visit counts, which is not what makes a
join-the-shortest-queue model behave as it does.
The DISTRIBUTIONS are deliberately kept: JMVA consumes a mean service
demand, so any renewal law with a finite mean is admissible, exactly as
it is for SolverMVA and SolverNC. Cache has no JMVA element either.
"""
return SolverJMT.getFeatureSet() - {
'Cache', 'CacheClassSwitcher',
'ReplacementStrategy_RR', 'ReplacementStrategy_FIFO',
'ReplacementStrategy_SFIFO', 'ReplacementStrategy_LRU',
'ReplacementStrategy_CLIMB',
'Fork', 'Join', 'Forker', 'Joiner', 'JoinPartial',
'ForkFanoutVector', 'ForkFanoutRandom', 'ForkBranchProbability',
'Place', 'Transition', 'Enabling', 'Inhibiting', 'Timing',
'Firing', 'Storage',
'Region',
'Reneging', 'Balking',
'SetupDelayOff', 'ServerParallelism', 'HeteroServers',
'SchedStrategy_DPS', 'SchedStrategy_GPS', 'SchedStrategy_HOL',
'SchedStrategy_PSPRIO', 'SchedStrategy_DPSPRIO', 'SchedStrategy_GPSPRIO',
'SchedStrategy_LCFSPI', 'SchedStrategy_LCFSPIPRIO',
'SchedStrategy_LCFSPRIO', 'SchedStrategy_LCFSPRPRIO',
'SchedStrategy_FCFSPR', 'SchedStrategy_FCFSPI',
'SchedStrategy_FCFSPRPRIO', 'SchedStrategy_FCFSPIPRIO',
'SchedStrategy_SEPT', 'SchedStrategy_LEPT',
'SchedStrategy_SJF', 'SchedStrategy_LJF',
'SchedStrategy_SRPT', 'SchedStrategy_SRPTPRIO',
'SchedStrategy_LPS', 'SchedStrategy_POLLING',
'SchedStrategy_EDD', 'SchedStrategy_EDF',
'RoutingStrategy_RROBIN', 'RoutingStrategy_WRROBIN',
'RoutingStrategy_JSQ', 'RoutingStrategy_SQ',
# the JMVA document has no capacity element at all
'FiniteCapacity',
}
get_jmva_feature_set = getJMVAFeatureSet
[docs]
def getMethodFeatureSet(self, method):
"""SolverJMT drives TWO ENGINES, and they accept different models.
'default' and 'jsim' run the JSIM SIMULATOR, whose
envelope is ``getFeatureSet``. The 'jmva.*' names run the JMVA
ANALYTICAL engine, which reads a document carrying only a station type,
a per-chain demand, a per-chain visit count, the populations or arrival
rates and a reference station -- so declaring the JSIM envelope for
jmva was a promise the writer could not keep.
Defining this is also what lets the base runAnalyzerChecks gate name the
offending features (mirrors MATLAB SolverJMT.getMethodFeatureSet):
without it the coarse supports(model) is used, which accepts every
model. A non-Network model (e.g. a LayeredNetwork) keeps the coarse
path and any structural checks that operate on such models."""
from ....lang.network import Network
model = getattr(self, 'model', None)
if not isinstance(model, Network):
return None
if str(method or '').lower().startswith('jmva'):
feats = SolverJMT.getJMVAFeatureSet()
if jmva_is_closed_only(method):
# RECAL, CoMoM, Chow, Bard-Schweitzer, AQL, Linearizer and De
# Souza-Muntz Linearizer are closed-network algorithms: JMT
# answers an open or a mixed model with "The selected solver
# cannot handle open classes" and a load-dependent one with the
# matching refusal. Exact MVA, which 'jmva' and 'jmva.mva'
# select, serves both.
# the eight closed-form algorithms are single-server ones
# (jmt_method_refusal words it); exact MVA carries the count
feats = feats - {'OpenClass', 'LoadDependence', 'MultiServer'}
return feats
return SolverJMT.getFeatureSet()
get_method_feature_set = getMethodFeatureSet
[docs]
def resolveMethod(self, options):
"""'default' runs as 'jsim', the simulator; every other name as itself.
MATLAB SolverJMT.resolveMethod."""
method = getattr(options, 'method', 'default') if options is not None else 'default'
return 'jsim' if str(method).lower() == 'default' else method
resolve_method = resolveMethod
[docs]
def supportsModelMethod(self, method):
"""Structural gate for what no registry name can state.
Three rules: the removed 'replication' method name,
the single-server restriction of the closed-form JMVA algorithms (a
server count is not a declared feature), and the one feature JMT admits
in a RE-ENCODED form only. Limited load dependence has no
representation of its own in either JMT document: the JSIM writer turns
it into a server count and the JMVA writer into the matching
<ldstation>, so alpha(n) = min(n,c) with an integer c is written
exactly and any other scaling would be solved at a service rate JMT
never saw. The first two come from ``jmt_method_refusal``, which the
analyzer asks as well."""
from ....constants import GlobalConstants
from ....lang.network import Network
model = getattr(self, 'model', None)
if isinstance(model, Network):
# The same predicate the JMVA arm and the analyzer ask, so
# this gate and those runs cannot answer differently.
structural = jmt_method_refusal(model.getStruct(), method, self.options)
if structural:
return False, structural
lld = model.getStruct().lldscaling if isinstance(model, Network) else None
if lld is not None:
lld = np.atleast_2d(np.asarray(lld, dtype=float))
for ist in range(lld.shape[0] if lld.ndim == 2 else 0):
alpha = lld[ist, :]
if alpha.size == 0 or np.all(alpha == 1.0):
continue
c = float(np.max(alpha))
shape = np.minimum(np.arange(1, alpha.size + 1, dtype=float), c)
if c != round(c) or c < 1 or np.any(np.abs(alpha - shape) > GlobalConstants.Zero):
return False, (
'Station %d uses a load-dependent scaling that is not the multiserver '
'encoding alpha(n) = min(n,c): JMT has no representation for it, since '
'both the JSIM and the JMVA writer carry the scaling as a server count, '
'and the model would be solved at the nominal service rate. Use '
'SolverCTMC, SolverNC, SolverMVA or SolverSSA, which read sn.lldscaling '
'directly.' % (ist + 1))
return super().supportsModelMethod(method)
supports_model_method = supportsModelMethod
[docs]
def unsupported_refusal_id(self, method, reason):
"""A cache refusal keeps the id the JSIM writer throws for it; else the base rule."""
if jmt_structural_refusal_id(self._sn, reason) == 'LINE:unsupported:Cache':
return 'LINE:unsupported:Cache'
return super().unsupported_refusal_id(method, reason)
[docs]
@staticmethod
def supports(model) -> bool:
"""Check if this solver supports the given model.
Mirrors MATLAB SolverJMT.supports. This previously returned True
unconditionally, so it accepted models built on features the JSIM writer
cannot represent (e.g. FCFSPR).
"""
from ...base import supports_via_featureset
from ....lang.network import Network
if not supports_via_featureset(SolverJMT, model):
return False
# The structural rules the analyzer refuses for the default (JSIM) method, so a
# model this answers True for is one JMT will run (MATLAB SolverJMT.supports).
if isinstance(model, Network):
return not jmt_method_refusal(model.getStruct(), 'default', SolverJMT.defaultOptions())
return True
[docs]
@staticmethod
def defaultOptions() -> OptionsDict:
"""Get default solver options."""
return OptionsDict({
'method': 'default',
'samples': 10000,
'seed': 23000,
'config': {'replications': JMT_DEFAULT_REPLICATIONS},
'verbose': default_verbose(),
'keep': False,
'conf_int': 0.99,
'max_rel_err': 0.03,
})
# =========================================================================
# File Management Methods (Gap 3a)
# =========================================================================
def _ensureTempDir(self, solvername: str) -> str:
"""Allocate, once, a private temporary directory for this solver instance.
Mirrors MATLAB getJSIMTempPath/getJMVATempPath, which set self.filePath =
lineTempName(solvername) on first use and reuse it thereafter. The
allocation must be cached: callers write a model with writeJSIM() and then
read the path back, so handing out a fresh directory per call would point
them at an empty one. The directory must also be private to the instance,
because the previous fixed fallback (/tmp/line_jmt/model.jsimg) made two
concurrent writeJSIM() calls collide on one file.
Whichever accessor runs first fixes the directory for both, exactly as in
MATLAB, where the second call finds self.filePath already set.
"""
if not getattr(self, '_temp_dir', None):
from ....api.io.sysutils import line_temp_name
self._temp_dir = line_temp_name(solvername)
return self._temp_dir
[docs]
def getFileName(self) -> str:
"""Get the model file name, WITHOUT directory and WITHOUT extension.
Matches MATLAB getFileName.m, whose callers build the name as
[fileName '.jsim'], and the JAR, which does fileName + ".jsim". This
previously returned 'model.jsimg', i.e. it included the extension, so it
did not compose the way the other two codebases' callers expect.
"""
return 'model'
[docs]
def getFilePath(self) -> str:
"""Get the directory holding the model file.
Returns the DIRECTORY, matching MATLAB getFilePath (out = self.filePath);
the file name is getFileName() and the joined path is getJSIMTempPath().
"""
return self._ensureTempDir('jsim')
[docs]
@staticmethod
def getJMTJarPath() -> str:
"""Get path to JMT.jar."""
from ....api.solvers.jmt.handler import _get_jmt_jar_path
return _get_jmt_jar_path()
[docs]
def getJMVATempPath(self) -> str:
"""Get path to the temporary JMVA model file."""
return os.path.join(self._ensureTempDir('jmva'), 'model.jmva')
[docs]
def getJSIMTempPath(self) -> str:
"""Get path to the temporary JSIM model file."""
return os.path.join(self._ensureTempDir('jsim'), 'model.jsim')
# =========================================================================
# Export Methods (Gap 3b)
# =========================================================================
@_instance_or_static
def writeJMVA(self, sn=None, outputFileName: str = None, options=None) -> str:
"""Write model to JMVA format.
Callable three ways:
solver.writeJMVA() # temp path, this solver's model
solver.writeJMVA('out.jmva') # this solver's model
SolverJMT.writeJMVA(sn, 'out.jmva', opts) # MATLAB/JAR form, no solver needed
The last is MATLAB's ``writeJMVA(sn, outputFileName, options)`` and the JAR's
``static String writeJMVA(NetworkStruct sn, String outputFileName, SolverOptions
options)``; a solver instance also accepts it.
Args:
sn: NetworkStruct to write; None means this solver's model. A string here is
taken as outputFileName, which keeps the one-argument instance form.
outputFileName: Output file path. If None, writes to a temp directory.
options: SolverOptions-like object or dict supplying ``method`` and ``samples``;
None means this solver's options, or the defaults with no solver.
Returns:
Path to the written file.
"""
from ..solver_lqns.jmva_writer import write_jmva
if isinstance(sn, (str, os.PathLike)) and outputFileName is None:
sn, outputFileName = None, sn
if sn is None:
if self is None:
raise ValueError('SolverJMT.writeJMVA needs a NetworkStruct when called without a solver.')
if self._sn is None:
self._extract_network_params()
sn = self._sn
if options is None and self is not None:
options = {'method': self.method, 'samples': self.options.samples}
if outputFileName is None:
if self is not None:
outputFileName = self.getJMVATempPath()
else:
outputFileName = os.path.join(tempfile.mkdtemp(prefix='line_jmva_'), 'model.jmva')
folder = os.path.dirname(outputFileName)
if folder:
os.makedirs(folder, exist_ok=True)
write_jmva(sn, outputFileName, _jmva_writer_options(options))
return outputFileName
[docs]
def writeJSIM(self, sn=None, outputFileName: str = None) -> str:
"""Write model to JSIM XML format.
Callable as ``solver.writeJSIM(outputFileName)`` or, as in MATLAB
(``writeJSIM(self, sn, outputFileName)``) and the JAR
(``writeJSIM(NetworkStruct sn, String outputFileName)``), as
``solver.writeJSIM(sn, outputFileName)``.
Args:
sn: NetworkStruct to write; None means this solver's model. A string here is
taken as outputFileName, which keeps the one-argument form.
outputFileName: Output file path. If None, writes to temp directory.
Returns:
Path to the written file.
"""
# handler exports no underscored SolverJMTOptions name; import the real one.
from ....api.solvers.jmt.handler import _write_jsim_file
if isinstance(sn, (str, os.PathLike)) and outputFileName is None:
sn, outputFileName = None, sn
if sn is None:
if self._sn is None:
self._extract_network_params()
sn = self._sn
if outputFileName is None:
outputFileName = self.getJSIMTempPath()
folder = os.path.dirname(outputFileName)
if folder:
os.makedirs(folder, exist_ok=True)
handler_options = _SolverJMTOptions(
method=self.method,
samples=self.options.samples,
seed=self.options.seed,
max_simulated_time=self.options.max_simulated_time,
conf_int=self.options.conf_int,
max_rel_err=self.options.max_rel_err,
max_events=self.maxEvents,
)
_write_jsim_file(sn, outputFileName, handler_options, model=self.model)
return outputFileName
[docs]
def QN2JSIMG(self, outputFileName: str = None) -> str:
"""Convert queueing network to JSIMG format. Wrapper for writeJSIM.
Args:
outputFileName: Output file path.
Returns:
Path to the written file.
"""
return self.writeJSIM(outputFileName)
# =========================================================================
# Transient Methods (Gap 3c)
# =========================================================================
[docs]
def getTranCdfRespT(self, R=None):
"""Get transient CDF of response times.
The same logged pipeline as getCdfRespT WITHOUT the steady-state seed:
the reference @SolverJMT/getTranCdfRespT.m starts the logged run from
the model's default initial state, so the collected samples cover the
transient, where getCdfRespT preloads the rounded steady-state queue
lengths to shorten the warmup.
"""
return self._cdfRespTPipeline(R, init_from_steady=False)
[docs]
def getTranCdfPassT(self, R=None):
"""Get transient CDF of passage times. Delegates to getTranCdfRespT,
its own name in the reference's sibling file."""
return self.getTranCdfRespT(R)
def _ensureAvgResults(self):
"""The steady-state accessors' gate, refusing where the analyzer would run the ensemble.
'default' on a finite timespan runs the transient ensemble, which yields no steady-state
averages (as in MATLAB); running it again per accessor and then reading an empty result
would be both slow and an AttributeError, so the accessor is told what to ask instead.
"""
if (self._resultStore() is None and self.options.method == 'default'
and np.isfinite(float(self.options.max_simulated_time))):
raise LineUnsupportedError(
"SolverJMT 'default' on a finite timespan runs the transient ensemble "
"(options.config.replications), which has no steady-state averages: use "
"getTranAvg, method 'jsim' for one run over the horizon, or an infinite timespan.", error_id='LINE:unsupported:Analysis')
return super()._ensureAvgResults()
def _require_finite_timespan(self, what):
"""Raise unless options.timespan is finite: a transient at an unstated horizon is not a quantity."""
if not np.isfinite(float(self.options.max_simulated_time)):
raise LineInputError(
"%s in SolverJMT requires to specify a finite timespan T, e.g., "
"SolverJMT(model, timespan=[0, T])." % what, error_id='LINE:input:Option')
def _sample_replications(self):
"""The transient ensemble: ``options.replications`` seeded system trajectories.
Replication k (1-based) runs at seed + k - 1, as MATLAB does. Returns the valid
(t, state) pairs and the union time grid capped at the MINIMUM of their maxima, so
previous-neighbour interpolation never runs past the shortest replication. The
requested and valid counts are kept on ``self._tran_replications``.
"""
reps = jmt_replications(getattr(self.options, 'replications', JMT_DEFAULT_REPLICATIONS))
init_seed = self.options.seed
paths = []
tumax = float('inf')
grid = set()
try:
for it in range(1, reps + 1):
self.options.seed = init_seed + it - 1
try:
path = self.sampleSysAggr()
except Exception as exc: # noqa: BLE001 - reported via line_warning, replication skipped
line_warning("SolverJMT", "Replication %d failed (%s), skipping.", it, exc)
continue
if not path or path.get('t') is None or len(path['t']) == 0:
line_warning("SolverJMT",
"Replication %d produced empty/invalid time series, skipping.", it)
continue
tv = np.asarray(path['t'], dtype=float).ravel()
paths.append((tv, path['state']))
tumax = min(tumax, float(np.max(tv)))
grid.update(tv.tolist())
finally:
self.options.seed = init_seed
self._tran_replications = {'requested': reps, 'valid': len(paths)}
if not paths:
raise RuntimeError("No valid replications produced. Cannot compute transient averages.")
tu = np.array(sorted(v for v in grid if v <= tumax), dtype=float)
return paths, tu
[docs]
def getTranProbAggr(self, node):
"""Empirical transient probability of the per-class aggregate states at a station.
Port of MATLAB @SolverJMT/getTranProbAggr.m: ``options.replications`` seeded runs
over the finite timespan, each interpolated (previous-neighbour) onto the union
grid; pi_t[k, j] is the fraction of replications whose state at t[k] is states[j].
Args:
node: Node object, or a 0-based station index.
Returns:
TranProbResult (Pi_t, SSnode_a): Pi_t is [t, prob] with prob len(t) x nstates, and
SSnode_a the unique aggregate states, one row each, one column per class. The parts
also read by name: 't', 'prob', 'states', 'Pi_t'.
"""
from ....lang.base import NodeType
self._require_finite_timespan('getTranProbAggr')
sn = self._sn if self._sn is not None else self.model.get_struct()
if isinstance(node, (int, np.integer)):
ist = int(node)
else:
nidx = getattr(node, '_node_index', None)
if nidx is None:
nidx = self.model.get_nodes().index(node)
ist = int(sn.nodeToStation[int(nidx)])
if ist < 0 or ist >= sn.nstations:
raise ValueError("getTranProbAggr in SolverJMT requires the input node to be a station.")
nt = sn.nodetype[int(sn.stationToNode[ist])]
if (nt.value if hasattr(nt, 'value') else int(nt)) == NodeType.SOURCE.value:
raise ValueError("getTranProbAggr in SolverJMT does not apply to Source nodes.")
paths, tu = self._sample_replications()
reps = self._tran_replications['requested']
K = sn.nclasses
nT = len(tu)
# (replication, time, class); NaN where a replication has no sample yet
node_states = np.full((len(paths), nT, K), np.nan)
for r, (tv, states) in enumerate(paths):
block = np.asarray(states[ist], dtype=float) if ist < len(states) else None
if block is None or block.ndim != 2 or block.size == 0 or np.isinf(block).any():
continue
idx = np.searchsorted(tv, tu, side='right') - 1
valid = idx >= 0
node_states[r, valid, :] = block[idx[valid], :K]
flat = node_states.reshape(-1, K)
flat = flat[~np.isnan(flat).any(axis=1)]
ss = np.unique(flat, axis=0) if flat.size else np.zeros((0, K))
pi_t = np.zeros((nT, ss.shape[0]))
if ss.shape[0] > 0:
for r in range(len(paths)):
for t_idx in range(nT):
row = node_states[r, t_idx, :]
if np.isnan(row).any():
continue
j = int(np.flatnonzero(np.all(ss == row, axis=1))[0])
pi_t[t_idx, j] += 1.0
# over the REQUESTED count, as MATLAB divides by iter_max: a failed replication is not dropped silently
pi_t /= reps
from ....constants import TranProbResult
return TranProbResult(np.column_stack([tu, pi_t]), ss)
def _runTransientReplications(self):
"""Transient averages over the ensemble of ``options.replications`` seeded JSIM runs.
Samples the seeded system trajectories (``_sample_replications``) and averages them
on the union of their time grids.
Utilization is read as ``min(n, c)/c`` at a finite server and as the raw
queue length at a delay; throughput follows it as ``U*c*mu`` and ``U*mu``.
Mirrors runTransientReplications of MATLAB @SolverJMT/runAnalyzer.m.
"""
import time as _time
from ....lang.base import NodeType
sn = self._sn
M, K = sn.nstations, sn.nclasses
self._require_finite_timespan('The transient ensemble')
t0 = _time.time()
paths, tu = self._sample_replications()
nvalid = len(paths)
QNt, UNt, TNt = {}, {}, {}
nservers = np.asarray(sn.nservers, dtype=float).ravel()
for i in range(M):
c = nservers[i] if i < len(nservers) else 1.0
for k in range(K):
q = np.zeros(len(tu))
u = np.zeros(len(tu))
for tv, states in paths:
st = states[i] if i < len(states) else None
if st is None:
continue
st = np.asarray(st, dtype=float)
if st.ndim != 2 or k >= st.shape[1] or not np.isfinite(st).any():
continue
col = st[:, k]
# previous-neighbour interpolation; a grid point before the
# first sample has no predecessor and reads 0, which is what
# the reference's own NaN-to-zero step leaves
idx = np.searchsorted(tv, tu, side='right') - 1
valid = idx >= 0
qv = np.zeros(len(tu))
qv[valid] = np.nan_to_num(col[idx[valid]])
q += qv / nvalid
if np.isfinite(c) and c > 0:
uv = np.zeros(len(tu))
uv[valid] = np.nan_to_num(np.minimum(col[idx[valid]], c) / c)
else:
uv = qv
u += uv / nvalid
rate = float(sn.rates[i, k]) if sn.rates is not None else 0.0
if not np.isfinite(rate):
rate = 0.0
scale = (c * rate) if np.isfinite(c) else rate
tput = u * scale
if sn.nodetype is not None and sn.nodetype[int(sn.stationToNode[i])] == NodeType.SOURCE:
# no queue is logged at a Source: it holds no jobs and emits at its arrival
# rate, as the MATLAB reference's transient arm reports it
q, u = np.zeros(len(tu)), np.zeros(len(tu))
tput = np.full(len(tu), rate if rate > 0 else 0.0)
QNt[(i, k)] = {'t': tu.copy(), 'metric': q}
UNt[(i, k)] = {'t': tu.copy(), 'metric': u}
TNt[(i, k)] = {'t': tu.copy(), 'metric': tput}
self._tran_avg = (QNt, UNt, TNt)
self._tran_runtime = _time.time() - t0
return self
[docs]
def getTranAvg(self):
"""Get transient average metrics from simulation.
Runs simulation with logging and extracts time series.
Returns:
Tuple of (QNt, UNt, TNt) time series dicts, or None if unavailable.
"""
# The transient ensemble produced the real transient mean; return it rather
# than the constant series the steady-state fallback below builds.
tran = getattr(self, '_tran_avg', None)
if tran is not None:
return tran
# Transient analysis runs the simulation with a finite timespan.
self._require_finite_timespan('getTranAvg')
if self.options.method == 'default':
# 'default' on a finite horizon IS the ensemble (runAnalyzer dispatches it)
self.runAnalyzer()
return self._tran_avg
if self._result is None:
self._ensureAvgResults()
# JMT steady-state results don't have time series natively.
# Return steady-state values as constant time series.
if self._result is None:
return None, None, None
M = self._sn.nstations
K = self._sn.nclasses
t = np.array([0.0, self.options.max_simulated_time if np.isfinite(self.options.max_simulated_time) else 1.0])
QNt = {}
UNt = {}
TNt = {}
for i in range(M):
for k in range(K):
q_val = self._result.Q[i, k] if self._result.Q is not None else 0.0
u_val = self._result.U[i, k] if self._result.U is not None else 0.0
t_val = self._result.T[i, k] if self._result.T is not None else 0.0
QNt[(i, k)] = {'t': t.copy(), 'metric': np.array([q_val, q_val])}
UNt[(i, k)] = {'t': t.copy(), 'metric': np.array([u_val, u_val])}
TNt[(i, k)] = {'t': t.copy(), 'metric': np.array([t_val, t_val])}
return QNt, UNt, TNt
# =========================================================================
# Probability Methods (Gap 3d)
# =========================================================================
[docs]
def getProb(self, node=None, state=None):
"""Get state probability from simulation trajectory.
Args:
node: Node index or node object.
state: Target state vector. If None, returns probability of current state.
Returns:
Float probability value, or dict of probabilities.
"""
if self._result is None:
self._ensureAvgResults()
# For simulation-based solver, compute from trajectory
result = self.sampleSysAggr(num_events=self.options.samples)
if result is None:
return 0.0
if node is None:
return self.getProbSys()
node_idx = node if isinstance(node, int) else getattr(node, '_station_index', 0)
if state is None:
# Use current state from sn
if self._sn is not None and hasattr(self._sn, 'state') and self._sn.state is not None:
state = self._sn.state[node_idx] if node_idx < len(self._sn.state) else None
if state is None or result is None:
return 0.0
# Time-weighted probability
if 'states' in result and node_idx < len(result.get('states', [])):
trajectory = result['states'][node_idx]
t = result.get('t', np.array([]))
if len(t) < 2:
return 0.0
target = np.asarray(state).flatten()
total_time = 0.0
match_time = 0.0
for idx in range(len(t) - 1):
dt = t[idx + 1] - t[idx]
total_time += dt
if np.allclose(trajectory[idx], target, atol=1e-10):
match_time += dt
return match_time / total_time if total_time > 0 else 0.0
return 0.0
prob = getProb
[docs]
def getProbSys(self):
"""Get joint system state probability.
Returns:
Float probability of the current system state.
"""
return self.getProbSysAggr()
[docs]
def getProbMarg(self, node=None, jobclass=None):
"""Get marginal state probability for a specific class at a node.
Args:
node: Node index or object.
jobclass: Class index.
Returns:
Dict mapping state values to probabilities.
"""
result = self.sampleSysAggr(num_events=self.options.samples)
if result is None:
return {}
node_idx = node if isinstance(node, int) else getattr(node, '_station_index', 0)
class_idx = jobclass if isinstance(jobclass, int) else getattr(jobclass, '_index', 0)
if 'states' not in result or node_idx >= len(result.get('states', [])):
return {}
trajectory = result['states'][node_idx]
t = result.get('t', np.array([]))
if len(t) < 2:
return {}
# Compute time-weighted histogram for the given class
prob_map = {}
total_time = 0.0
for idx in range(len(t) - 1):
dt = t[idx + 1] - t[idx]
total_time += dt
state_val = int(trajectory[idx][class_idx]) if class_idx < len(trajectory[idx]) else 0
prob_map[state_val] = prob_map.get(state_val, 0.0) + dt
if total_time > 0:
for k in prob_map:
prob_map[k] /= total_time
return prob_map
[docs]
def getProbNormConstAggr(self):
"""Log normalizing constant, from the JMVA engine only.
A simulation computes no normalizing constant, but the analytical JMVA
algorithms report one in the result file's <normconst logValue>, which
MATLAB stores as result.Prob.logNormConstAggr. It is returned here for
the jmva* methods and refused for the simulation ones rather than
handing back the NaN placeholder.
Raises:
NotImplementedError: on the simulation methods.
"""
from ....api.solvers.jmt.handler import _is_jmva_method
if not _is_jmva_method(self.method):
raise LineUnsupportedError(
"getProbNormConstAggr() is not supported by SolverJMT with method='%s'. "
"Use an analytical method (SolverJMT 'jmva'), SolverNC or SolverCTMC "
"for normalizing constant computation." % self.method, error_id='LINE:unsupported:Analysis')
if str(self.method).lower() in ('jmva', 'jmt.jmva'):
# JMVA's plain MVA algorithm writes logValue="NaN"; RECAL is exact and computes G, as in MATLAB.
return SolverJMT(self.model, method='jmva.recal', verbose=False).getProbNormConstAggr()
if self._result is None:
self._ensureAvgResults()
return getattr(self._result, 'logNormConstAggr', float('nan'))
# =========================================================================
# Sampling Methods (Gap 3e)
# =========================================================================
[docs]
def sample(self, node, numEvents: int = 1000):
"""Sample the aggregated state trajectory at a node (JMT logs only carry
per-class counts, so this is equivalent to :meth:`sampleAggr`).
Args:
node: Node index or node object.
numEvents: Number of events to sample.
Returns:
SamplePath (as sampleAggr): .t, .state; also readable as sp['t'].
"""
return self.sampleAggr(node, numEvents)
def _resolve_node(self, node):
"""Resolve a node argument to (node_index_0based, node_name)."""
sn = self._sn if self._sn is not None else self.model.get_struct()
if isinstance(node, int):
node_idx = node
node_name = sn.nodenames[node_idx]
else:
node_name = node.name if hasattr(node, 'name') else str(node)
node_idx = self.model.get_node_index(node_name) - 1 # 0-based
return node_idx, node_name
def _cleanup_log_dir(self, path):
"""Remove a temporary JMT sampling log directory unless options.keep is set.
Mirrors the JAR SolverJMT.cleanupDir and MATLAB @SolverJMT/runAnalyzer.m
(keep=false), so log-based sampling does not accumulate scratch folders.
"""
if getattr(self.options, 'keep', False):
return
if not path:
return
shutil.rmtree(path, ignore_errors=True)
def _run_logged_copy(self, is_node_logged, numEvents):
"""Simulate a logged copy of the model and return (modelCopy, log_path).
Mirrors MATLAB SolverJMT.sampleAggr: build a temporary copy, apply
linkAndLog to the requested nodes, and run JMT so the per-node
arrival/departure CSV logs are produced.
"""
import tempfile
sn = self._sn if self._sn is not None else self.model.get_struct()
model_copy = self.model.copy()
model_copy.reset_network()
Plinked = None
if hasattr(self.model, 'get_linked_routing_matrix'):
Plinked = self.model.get_linked_routing_matrix()
if Plinked is None and hasattr(sn, 'rtorig'):
Plinked = sn.rtorig
if Plinked is None:
raise RuntimeError("JMT log-based sampling requires the routing "
"matrix (rtorig).")
log_path = tempfile.mkdtemp(prefix='jmt_sample_logs_')
model_copy.link_and_log(Plinked, is_node_logged, log_path)
# ONE JSIM run over the caller's horizon: the whole options object crosses (seed, horizon,
# backend), and 'jsim' rather than 'default' keeps the copy out of the transient ensemble.
jmt = SolverJMT(model_copy, replace(self.options, method='jsim'))
# JMT cannot cap events per node; scale the total event budget as MATLAB
# does (numEvents * nnodes * nclasses) so each node yields ~numEvents.
if numEvents and numEvents > 0:
jmt.maxEvents = int(numEvents) * sn.nnodes * sn.nclasses
jmt.runAnalyzer()
return model_copy, log_path
def _node_preload(self, node_idx, K):
"""Initial per-class population at a node (INIT row of the trajectory)."""
from ....api.sn.transforms import sn_get_state_aggr
sn = self._sn if self._sn is not None else self.model.get_struct()
preload = np.zeros(K)
preload_map = sn_get_state_aggr(sn)
isf = int(sn.nodeToStateful[node_idx]) if sn.nodeToStateful is not None else -1
if isf in preload_map and preload_map[isf] is not None:
pl = np.ravel(np.asarray(preload_map[isf], dtype=float))
preload[:min(K, len(pl))] = pl[:min(K, len(pl))]
return preload
[docs]
def sampleAggr(self, node, numEvents: int = 1000):
"""Sample the aggregated (per-class count) state trajectory at a node.
Faithful port of MATLAB SolverJMT.sampleAggr: a temporary logged copy of
the model is simulated and its arrival/departure logs are reconstructed
into a piecewise-constant per-class queue-length trajectory.
Args:
node: Node index (0-based) or node object.
numEvents: Desired number of sampled events at the node.
Returns:
SamplePath (fields also readable as keys) with 'handle', 't' (event-boundary times), 'state'
(len x nclasses per-class counts), 'event' (chronological event
list) and 'isaggregate'=True.
"""
import os
from ....api.solvers.jmt.handler import parse_tran_state
sn = self._sn if self._sn is not None else self.model.get_struct()
K = sn.nclasses
node_idx, node_name = self._resolve_node(node)
nnodes = self.model.get_number_of_nodes()
is_node_logged = [False] * nnodes
is_node_logged[node_idx] = True
model_copy, log_path = self._run_logged_copy(is_node_logged, numEvents)
preload = self._node_preload(node_idx, K)
class_names = [sn.classnames[r] for r in range(K)] if sn.classnames is not None else None
arv_file = os.path.join(log_path, f"{node_name}-Arv.csv")
dep_file = os.path.join(log_path, f"{node_name}-Dep.csv")
state, evtype, evclass, evjob = parse_tran_state(
arv_file, dep_file, preload, class_names)
# see _kb/06-solver-catalog.md (Wrappers: "JMT python wrapper:
# export/import internals") for the event-boundary time shift
t_all = state[:, 0]
_, uniq_idx = np.unique(t_all, return_index=True)
uniq_idx = np.sort(uniq_idx)
t = t_all[uniq_idx]
qlen = state[uniq_idx, 1:1 + K]
if len(t) > 1:
t_shift = np.concatenate([t[1:], t[-1:]])
else:
t_shift = t.copy()
m = min(len(t_shift), 1 + numEvents) if (numEvents and numEvents > 0) else len(t_shift)
t_shift = t_shift[:m]
qlen = qlen[:m, :]
t_out = np.concatenate([[0.0], t_shift[:-1]]) if len(t_shift) > 0 else t_shift
tmax = float(t_out[-1]) if len(t_out) > 0 else np.inf
events = []
for e in range(len(evtype)):
te = float(state[e, 0])
if te > tmax:
continue
events.append({
'event': evtype[e],
'node': node_idx,
'class': (int(evclass[e]) if not np.isnan(evclass[e]) else None),
't': te,
'job': (int(evjob[e]) if not np.isnan(evjob[e]) else None),
})
self._cleanup_log_dir(log_path)
from line_solver.solvers.solver_ssa.solver_ssa import SamplePath
return SamplePath(handle=node, t=t_out, state=qlen, event=events, isaggregate=True)
[docs]
def sampleSys(self, numEvents: int = 1000):
"""Sample system-wide state trajectory.
Args:
numEvents: Number of events to sample.
Returns:
SamplePath (as sampleSysAggr): .t, .state; also readable as sp['t'].
"""
return self.sampleSysAggr(num_events=numEvents)
# =========================================================================
# Additional Metric Methods
# =========================================================================
[docs]
def getAvgResidT(self) -> np.ndarray:
"""Get average residence times (M x K)."""
if self._result is None:
self._ensureAvgResults()
if self._sn is not None and hasattr(self._sn, 'visits') and self._sn.visits:
return sn_get_residt_from_respt(self._sn, self._result.R, None)
return self._result.R.copy() if self._result.R is not None else np.array([])
[docs]
def getAvgWaitT(self) -> np.ndarray:
"""Get average waiting times (M x K). W = R - S."""
if self._result is None:
self._ensureAvgResults()
R = self._result.R.copy() if self._result.R is not None else np.array([])
if len(R) == 0:
return R
if hasattr(self._sn, 'rates') and self._sn.rates is not None:
rates = np.asarray(self._sn.rates)
S = np.zeros_like(rates)
nonzero = rates > 0
S[nonzero] = 1.0 / rates[nonzero]
W = R - S
W = np.maximum(W, 0.0)
return W
return R
[docs]
def getAvg(self):
"""Get all average metrics at once.
Returns:
Tuple of (Q, U, R, T, A, W)
"""
if self._result is None:
self._ensureAvgResults()
r = self._result
Q = r.Q if r.Q is not None else np.array([])
U = r.U if r.U is not None else np.array([])
R = r.R if r.R is not None else np.array([])
T = r.T if r.T is not None else np.array([])
A = r.A if r.A is not None else T.copy()
W = self.getAvgResidT()
return Q, U, R, T, A, W
[docs]
def getAvgNode(self):
"""Get average metrics per node.
Returns:
Tuple of (QNn, UNn, RNn, WNn, ANn, TNn) - node-level metrics
"""
if self._result is None:
self._ensureAvgResults()
sn = self._sn
I = sn.nnodes
M = sn.nstations
K = sn.nclasses
QN = self._result.Q if self._result.Q is not None else np.zeros((M, K))
UN = self._result.U if self._result.U is not None else np.zeros((M, K))
RN = self._result.R if self._result.R is not None else np.zeros((M, K))
TN = self._result.T if self._result.T is not None else np.zeros((M, K))
AN = self._result.A if self._result.A is not None else TN.copy()
WN = self.getAvgResidT() if len(RN) > 0 else np.zeros((M, K))
QNn = np.zeros((I, K))
UNn = np.zeros((I, K))
RNn = np.zeros((I, K))
WNn = np.zeros((I, K))
TNn = np.zeros((I, K))
ANn = np.zeros((I, K))
for ist in range(M):
ind = int(sn.stationToNode[ist])
if 0 <= ind < I:
QNn[ind, :] = QN[ist, :]
UNn[ind, :] = UN[ist, :]
RNn[ind, :] = RN[ist, :]
WNn[ind, :] = WN[ist, :]
TNn[ind, :] = TN[ist, :]
ANn, TNn = self._node_flows(sn, TN, AN)
return QNn, UNn, RNn, WNn, ANn, TNn
@staticmethod
def _node_flows(sn, TN, AN):
"""Node arrival and throughput rates, by MATLAB's generic NetworkSolver/getAvgNode.m rule.
SolverJMT has no getAvgNode of its own in MATLAB, so its node table is the generic one:
sn_get_node_arvr_from_tput, then sn_get_node_tput_from_tput, then the cache post-pass. The
latter routes each node's arrivals along rtnodes, so a ClassSwitch reports the class it
sends out (not the class that arrived, as Tput == ArvR did) and a Cache its hit/miss split.
JMT leaves a class it never measured at a station as NaN where MATLAB getResults.m leaves 0;
a closed chain then sums NaN over TN(refstat, inchain) and loses every non-station node of
the chain, so the node flows read the MATLAB convention.
"""
from line_solver.api.sn import sn_get_node_arvr_from_tput
from line_solver.api.sn.getters import sn_get_node_tput_from_tput, _cache_downstream_arvr
TN0 = np.nan_to_num(TN, nan=0.0)
ANn = sn_get_node_arvr_from_tput(sn, TN0, None, AN)
TNn = sn_get_node_tput_from_tput(sn, TN0, None, ANn)
return _cache_downstream_arvr(sn, ANn, TNn, TN0, None)
[docs]
def getCdfRespT(self, R=None):
"""Get response time CDF via transient simulation with logging.
Runs JMT twice:
1. First run: Get steady-state queue lengths to initialize state
2. Second run: Run with logging enabled to collect response time samples
Ported from MATLAB's SolverJMT.getCdfRespT.
Args:
R: Optional response time handles (uses defaults if None)
Returns:
List of lists where RD[station][class] is a 2D array [cdf, time]
"""
return self._cdfRespTPipeline(R, init_from_steady=True)
def _cdfRespTPipeline(self, R=None, init_from_steady=True):
"""The shared logged-run pipeline behind getCdfRespT (seeded from the
rounded steady-state queue lengths) and getTranCdfRespT (unseeded)."""
import os
import tempfile
from ....api.solvers.jmt.handler import parse_tran_resp_t
sn = self._sn
M = sn.nstations
K = sn.nclasses
# Initialize result structure
RD = [[None for _ in range(K)] for _ in range(M)]
n = None
if init_from_steady:
# Step 1: Get steady-state queue lengths (first JMT run)
QN = self.getAvgQLen()
n = QN.copy()
# Adjust job numbers based on network constraints
for r in range(K):
if np.isinf(sn.njobs[r]):
# Open class - use floor of queue lengths
for i in range(M):
n[i, r] = np.floor(QN[i, r])
else:
# Closed class - ensure total population equals njobs
for i in range(M):
n[i, r] = np.floor(QN[i, r])
total_jobs = np.sum(n[:, r])
if total_jobs < sn.njobs[r]:
# Put remaining jobs on bottleneck station
imax = np.argmax(n[:, r])
n[imax, r] = n[imax, r] + sn.njobs[r] - total_jobs
# Step 2: Copy model for CDF computation
cdfmodel = self.model.copy()
cdfmodel.reset_network()
cdfmodel.reset()
# Determine which nodes should be logged (all stations except Source/Sink)
nnodes = cdfmodel.get_number_of_nodes()
is_node_logged = [False] * nnodes
for i in range(cdfmodel.get_number_of_stations()):
station = cdfmodel.get_stations()[i]
node_idx = cdfmodel.get_node_index(station.name) - 1 # Convert to 0-based
# Don't log Source or Sink
if hasattr(station, 'node_type'):
from ....lang.base import NodeType
if station.node_type not in (NodeType.SOURCE, NodeType.SINK):
is_node_logged[node_idx] = True
else:
is_node_logged[node_idx] = True
# Get original routing matrix from the model
Plinked = None
if hasattr(self.model, 'get_linked_routing_matrix'):
Plinked = self.model.get_linked_routing_matrix()
if Plinked is None and hasattr(sn, 'rtorig'):
Plinked = sn.rtorig
if Plinked is None:
raise RuntimeError("getCdfRespT requires routing matrix (rtorig)")
# Step 3: Set up logging
log_path = tempfile.mkdtemp(prefix='jmt_cdf_logs_')
cdfmodel.link_and_log(Plinked, is_node_logged, log_path)
# Initialize model state from marginal distribution (seeded route only;
# the transient getter starts from the model's default initial state)
if init_from_steady and n is not None:
try:
cdfmodel.init_from_marginal(n)
except (ValueError, NotImplementedError):
pass # May not be supported for all models
# Step 4: Run JMT on logged model (second JMT run)
# one logged run, so 'default' must not become the transient ensemble on a finite horizon
cdf_opts = replace(self.options, method='jsim') if self.options.method == 'default' else self.options
cdf_solver = SolverJMT(cdfmodel, cdf_opts)
cdf_solver.runAnalyzer()
# Step 5: Parse logs to get response time samples
node_names = self.model.get_node_names() if hasattr(self.model, 'get_node_names') else []
station_names = sn.nodenames if sn.nodenames else []
# Get class names in model order for deterministic mapping
class_names = [sn.classnames[r] for r in range(K)] if sn.classnames is not None else None
for i in range(M):
# Get original station name
node_idx = int(sn.stationToNode[i]) if sn.stationToNode is not None else i
station_name = station_names[node_idx] if node_idx < len(station_names) else f'Station{i}'
# Check if this node was logged
if node_idx < len(is_node_logged) and is_node_logged[node_idx]:
arv_file = os.path.join(log_path, f"{station_name}-Arv.csv")
dep_file = os.path.join(log_path, f"{station_name}-Dep.csv")
if os.path.exists(arv_file) and os.path.exists(dep_file):
class_resp_t, _, _ = parse_tran_resp_t(arv_file, dep_file, class_names=class_names)
for r in range(min(K, len(class_resp_t))):
resp_times = class_resp_t[r]
if len(resp_times) > 0:
# Create empirical CDF (ecdf equivalent, matching MATLAB ecdf)
sorted_times = np.sort(resp_times)
n = len(sorted_times)
# Group duplicate values (matches MATLAB ecdf and JAR createEmpiricalCDF)
unique_vals, counts = np.unique(sorted_times, return_counts=True)
cumulative = np.cumsum(counts) / n
# Add point at F=0 for the first data point
X = np.concatenate([[unique_vals[0]], unique_vals])
F = np.concatenate([[0.0], cumulative])
# Store as [F, X] format
RD[i][r] = np.column_stack([F, X])
self._cleanup_log_dir(log_path)
return RD
[docs]
def getPerctRespT(self, percentiles=None):
"""Get response time percentiles.
Args:
percentiles: Array of percentile values (default: [90, 95, 99])
Returns:
Tuple of (PercRT, PercTable) where PercRT is list of dicts
and PercTable is a pandas DataFrame
"""
import pandas as pd
if percentiles is None:
percentiles = np.array([90, 95, 99])
else:
percentiles = np.asarray(percentiles)
if self._result is None:
self._ensureAvgResults()
R = self._result.R
M = self._sn.nstations
K = self._sn.nclasses
PercRT = []
rows = []
perc_col_names = [f'P{int(p)}' for p in percentiles]
percentiles_normalized = percentiles / 100.0
station_names = self._sn.nodenames if self._sn.nodenames else [f'Station{i}' for i in range(M)]
class_names = self._sn.classnames if self._sn.classnames else [f'Class{r}' for r in range(K)]
for i in range(M):
for r in range(K):
if R is not None and i < R.shape[0] and r < R.shape[1]:
mean_resp_t = R[i, r]
if mean_resp_t > 0 and not np.isnan(mean_resp_t):
lambda_rate = 1.0 / mean_resp_t
perc_values = -np.log(1 - percentiles_normalized) / lambda_rate
PercRT.append({
'station': i + 1,
'class': r + 1,
'percentiles': percentiles.tolist(),
'values': perc_values.tolist(),
})
node_idx = int(self._sn.stationToNode[i]) if self._sn.stationToNode is not None else i
station_name = station_names[node_idx] if node_idx < len(station_names) else f'Station{i}'
row_data = {
'Station': station_name,
'Class': class_names[r] if r < len(class_names) else f'Class{r}',
}
for perc_col, perc_val in zip(perc_col_names, perc_values):
row_data[perc_col] = perc_val
rows.append(row_data)
PercTable = pd.DataFrame(rows) if rows else pd.DataFrame()
return PercRT, PercTable
def __repr__(self) -> str:
return f"SolverJMT(method='{self.method}', samples={self.options.samples})"
# Snake case aliases for MATLAB compatibility
avg_qlen = getAvgQLen
prob_sys_aggr = getProbSysAggr
prob_aggr = getProbAggr
avg_util = getAvgUtil
avg_respt = getAvgRespT
get_avg_respt = getAvgRespT
avg_residt = getAvgResidT
avg_waitt = getAvgWaitT
avg_tput = getAvgTput
avg_arv_r = getAvgArvR
avg_fcr = getAvgFcr
run_analyzer = runAnalyzer
get_runtime = getRuntime
get_method = getMethod
list_valid_methods = listValidMethods
is_available = isAvailable
get_feature_set = getFeatureSet
default_options = defaultOptions
cdf_resp_t = getCdfRespT
cdf_respt = getCdfRespT
get_cdf_resp_t = getCdfRespT
get_tran_cdf_respt = getTranCdfRespT
get_tran_cdf_resp_t = getTranCdfRespT
get_tran_cdf_pass_t = getTranCdfPassT
perct_resp_t = getPerctRespT
perct_respt = getPerctRespT
[docs]
def getAvgNodeTable(self) -> pd.DataFrame:
"""
Per-node average performance metrics. Mirrors MATLAB and JAR output by
including non-station nodes (Source, Sink, Router, ClassSwitch) alongside
station rows. Per-node arrival rates are computed via
sn_get_node_arvr_from_tput and throughputs via sn_get_node_tput_from_tput
(see _node_flows); queue/delay rows reuse the station-level metrics.
"""
if self._result is None:
self._ensureAvgResults()
from line_solver.api.sn import sn_get_node_arvr_from_tput
from line_solver.lang.base import NodeType
sn = self._sn
I = sn.nnodes
K = sn.nclasses
nodenames = sn.nodenames if sn.nodenames else [f'Node{i}' for i in range(I)]
classnames = sn.classnames if sn.classnames else [f'Class{r}' for r in range(K)]
nodetype = list(sn.nodetype) if sn.nodetype is not None else []
TN = self._result.T if self._result.T is not None else np.zeros((sn.nstations, K))
# see _kb/06-solver-catalog.md (Wrappers: "JMT python wrapper: export/import internals")
AN = self._result.A if self._result.A is not None else TN.copy()
ANn, TNn = self._node_flows(sn, TN, AN)
stationToNode = np.asarray(sn.stationToNode).flatten() if sn.stationToNode is not None else np.array([])
# Build rows: stations first (with full Q/U/R/ResidT) then non-station nodes
if sn.visits and self._result.R is not None:
WN = sn_get_residt_from_respt(sn, self._result.R, None)
else:
WN = self._result.R.copy() if self._result.R is not None else np.zeros((sn.nstations, K))
rows = []
station_node_set = set()
for ist in range(sn.nstations):
ind = int(stationToNode[ist]) if ist < len(stationToNode) else ist
station_node_set.add(ind)
for r in range(K):
qlen = self._result.Q[ist, r] if self._result.Q is not None else 0.0
util = self._result.U[ist, r] if self._result.U is not None else 0.0
respt = self._result.R[ist, r] if self._result.R is not None else 0.0
residt = WN[ist, r] if ist < WN.shape[0] and r < WN.shape[1] else respt
arvr = ANn[ind, r] if 0 <= ind < I else 0.0
tput = TNn[ind, r] if 0 <= ind < I else 0.0
if ind < len(nodetype):
nt_v = int(nodetype[ind]) if not hasattr(nodetype[ind], 'value') else int(nodetype[ind].value)
src_v = NodeType.SOURCE.value if hasattr(NodeType.SOURCE, 'value') else int(NodeType.SOURCE)
if nt_v == src_v:
qlen = util = respt = residt = 0.0
arvr = 0.0
rows.append({
'Node': nodenames[ind] if ind < len(nodenames) else f'Node{ind}',
'JobClass': classnames[r],
'QLen': qlen, 'Util': util, 'RespT': respt, 'ResidT': residt,
'ArvR': arvr, 'Tput': tput,
})
# Non-station nodes: emit one row per node x class with passthrough flow
for ind in range(I):
if ind in station_node_set or ind >= len(nodetype):
continue
for r in range(K):
arvr = ANn[ind, r]
tput = TNn[ind, r]
if arvr == 0 and tput == 0:
continue
rows.append({
'Node': nodenames[ind] if ind < len(nodenames) else f'Node{ind}',
'JobClass': classnames[r],
'QLen': 0.0, 'Util': 0.0, 'RespT': 0.0, 'ResidT': 0.0,
'ArvR': arvr, 'Tput': tput,
})
# see _kb/06-solver-catalog.md (Wrappers: "JMT python wrapper: export/import internals")
F = int(sn.nregions) if getattr(sn, 'nregions', 0) else 0
Qfcr = getattr(self._result, 'Qfcr', None)
if F > 0 and Qfcr is not None:
Rfcr = getattr(self._result, 'Rfcr', None)
Wfcr = getattr(self._result, 'Wfcr', None)
Tfcr = getattr(self._result, 'Tfcr', None)
region_names = []
if self.model is not None and hasattr(self.model, 'get_regions'):
region_names = [rg.get_name() if hasattr(rg, 'get_name') else f'FCR{f+1}'
for f, rg in enumerate(self.model.get_regions())]
for f in range(F):
rname = region_names[f] if f < len(region_names) else f'FCR{f+1}'
for r in range(K):
qlen = Qfcr[f, r] if Qfcr is not None else 0.0
respt = Rfcr[f, r] if Rfcr is not None else 0.0
residt = Wfcr[f, r] if Wfcr is not None else respt
tput = Tfcr[f, r] if Tfcr is not None else 0.0
if np.nansum([qlen, respt, tput]) <= 0:
continue
rows.append({
'Node': rname,
'JobClass': classnames[r],
'QLen': qlen, 'Util': np.nan, 'RespT': respt, 'ResidT': residt,
'ArvR': np.nan, 'Tput': tput,
})
df = pd.DataFrame(rows)
if not getattr(self, '_table_silent', False):
print(df.to_string(index=False))
return df
avg_node_table = getAvgNodeTable
get_avg_node_table = getAvgNodeTable
[docs]
def getAvgCacheTable(self) -> pd.DataFrame:
"""Per-class cache hit/miss table from JMT's Cache Hit Rate measures (see cache_table).
MATLAB reaches it through the generic NetworkSolver.getAvgCacheTable; the
split is the simulated one _apply_cache_hit_probs stored on the struct.
"""
if self._result is None:
self._ensureAvgResults()
from ...cache_table import build_cache_avg_table
return build_cache_avg_table(self)
get_avg_cache_table = getAvgCacheTable
avg_cache_table = getAvgCacheTable
get_file_name = getFileName
get_file_path = getFilePath
get_jmt_jar_path = getJMTJarPath
write_jmva = writeJMVA
write_jsim = writeJSIM