"""
Base distribution classes for LINE (pure Python).
This module provides the base classes for probability distributions
implemented in pure Python.
"""
from abc import ABC, abstractmethod
from typing import Any, Optional, Tuple, Union
import math
import numpy as np
from line_solver.constants import GlobalConstants
from line_solver.api.io.errors import LineUnsupportedError
from .._aliasing import alias_getattr as _dist_alias_getattr
def _late_alias(name):
"""Snake-case alias resolved on the INSTANCE. A plain ``get_skew = getSkew`` in a base class binds the base's
method at class creation, so every subclass overriding ``getSkew`` still answered the base's 0.0 or NotImplementedError."""
def alias(self, *args, **kwargs):
return getattr(self, name)(*args, **kwargs)
alias.__name__ = name
alias.__doc__ = "Alias of :meth:`%s`." % name
return alias
[docs]
class Distribution(ABC):
"""
Base class for all probability distributions.
Provides the common interface for all probability distributions
in LINE, including service time distributions and inter-arrival time distributions.
"""
[docs]
def __init__(self):
"""Initialize a new distribution."""
self._name = 'Distribution'
def __getattr__(self, name):
# API-parity fallback: resolve documented alias spellings
# (snake/camel/accessor-prefix) to the implemented methods
return _dist_alias_getattr(self, name)
@property
def name(self) -> str:
"""Get the distribution name."""
return self._name
def getName(self) -> str:
"""Get the name of this distribution."""
return self._name
[docs]
def get_feature_name(self) -> str:
"""The SolverFeatureSet entry this distribution is marked under.
Separate from name because that one also selects the ProcessType and the
JSON wire type: a subclass whose registry name is more specific than its
process type (a Trace, or a two-phase Coxian) says so here without
moving to a different process type. Defaults to name, so a distribution
needing no distinction is unaffected.
"""
return self._name
# MATLAB/JAR-compatible alias
getFeatureName = get_feature_name
[docs]
@abstractmethod
def getMean(self) -> float:
"""Get the mean (expected value) of the distribution."""
pass
[docs]
@abstractmethod
def getVar(self) -> float:
"""Get the variance of the distribution."""
pass
# Abstract methods for snake_case aliases
[docs]
def get_mean(self) -> float:
"""Get the mean (expected value) of the distribution."""
return self.getMean()
[docs]
def get_var(self) -> float:
"""Get the variance of the distribution."""
return self.getVar()
[docs]
def eval_pmf(self, x: int) -> float:
"""Evaluate the probability mass function at point x (snake_case alias
for evalPMF; dispatches to the concrete distribution's override)."""
return self.evalPMF(x)
[docs]
def getSCV(self) -> float:
"""
Get the squared coefficient of variation (SCV).
SCV = Var[X] / E[X]^2
"""
mean = self.getMean()
if mean == 0:
return float('inf')
return self.getVar() / (mean ** 2)
[docs]
def getRate(self) -> float:
"""Get the rate parameter (1/mean).
For immediate (zero-mean) distributions, returns GlobalConstants.Immediate
(a large but finite value) to avoid numerical issues with infinite rates.
This matches MATLAB's behavior.
"""
mean = self.getMean()
if mean == 0:
return GlobalConstants.Immediate
return 1.0 / mean
[docs]
def getSkew(self) -> float:
"""Get the skewness of the distribution."""
return 0.0 # Default, override in subclasses
[docs]
def getSupport(self) -> Tuple[float, float]:
"""Get the support range [min, max] of the distribution."""
return (0.0, float('inf'))
[docs]
def isContinuous(self) -> bool:
"""Check if this distribution is continuous."""
return True
[docs]
def isDiscrete(self) -> bool:
"""Check if this distribution is discrete."""
return False
[docs]
def isDisabled(self) -> bool:
"""Check if this distribution is equivalent to a Disabled distribution.
Mirrors MATLAB Distribution.isDisabled (isnan(getMean(self))): a
distribution whose mean is not a number carries no usable process
representation, so the struct records it as disabled rather than
mislabelling it. Returning a hardcoded False, as this did, made
refreshStruct treat such a distribution as an ordinary one.
"""
try:
mean = self.getMean()
except (NotImplementedError, ArithmeticError): # no usable mean
return False
try:
return bool(math.isnan(float(mean)))
except (TypeError, ValueError):
return False
[docs]
def evalCDF(self, x: float) -> float:
"""Evaluate the cumulative distribution function at point x."""
raise LineUnsupportedError("CDF not implemented for this distribution", error_id='LINE:unsupported:Analysis')
# eval_cdf cannot come from the generic snake->camel fallback, which would ask for evalCdf: the solvers'
# getCdfRespT family spells it Cdf, the distributions spell it CDF.
eval_cdf = _late_alias('evalCDF')
[docs]
def evalProbInterval(self, t0: float, t1: float) -> float:
"""Probability of falling in [t0, t1], F(t1) - F(t0) (Distribution.evalProbInterval)."""
if t1 < t0:
raise ValueError("CDF interval incorrectly specified (t1<t0)")
return self.evalCDF(t1) - self.evalCDF(t0)
eval_prob_interval = _late_alias('evalProbInterval')
[docs]
def evalPDF(self, x: float) -> float:
"""Evaluate the probability density function at point x."""
raise LineUnsupportedError("PDF not implemented for this distribution", error_id='LINE:unsupported:Analysis')
[docs]
def sample(self, n: int = 1, rng: Optional[np.random.Generator] = None) -> np.ndarray:
"""
Generate random samples from this distribution.
Args:
n: Number of samples to generate.
rng: Optional random number generator.
Returns:
Array of n random samples.
"""
raise LineUnsupportedError("Sampling not implemented for this distribution", error_id='LINE:unsupported:Analysis')
# Aliases (note: get_mean and get_var are already implemented as methods above)
get_name = _late_alias('getName')
get_scv = _late_alias('getSCV')
get_rate = _late_alias('getRate')
get_skew = _late_alias('getSkew')
get_support = _late_alias('getSupport')
[docs]
class ContinuousDistribution(Distribution):
"""Base class for continuous probability distributions."""
def __init__(self):
super().__init__()
[docs]
def isContinuous(self) -> bool:
return True
[docs]
def isDiscrete(self) -> bool:
return False
[docs]
def evalLST(self, s: float) -> float:
"""Laplace-Stieltjes transform E[e^{-sX}] evaluated at s.
Generic numerical fallback: rectangle-rule integration of evalPDF over
[0, 20*mean] with 1000 points. Subclasses with a closed form or a
MATLAB-matching numerical scheme (Pareto/Lognormal/Weibull) override
this. Used by the G/M/1 MVA solver's sigma-root for non-PH arrivals.
"""
import math
# A Markovian law has the CLOSED FORM alpha (sI - D0)^-1 (-D0 e). The
# rectangle rule below carries percent-level error on one (Exp(1) at
# s = 0.5 read 0.65672 against 2/3) and returns zero for a COMPLEX
# argument, since math.exp refuses it -- and complex arguments are what
# transform inversion and root location need. MRO puts this class ahead
# of Markovian for Exp and its siblings, so the dispatch is made here.
if hasattr(self, 'getD0') and hasattr(self, 'getInitProb'):
try:
D0 = np.atleast_2d(np.asarray(self.getD0(), dtype=float))
alpha = np.asarray(self.getInitProb(), dtype=float).ravel()
n = D0.shape[0]
if alpha.size == n and np.all(np.isfinite(D0)):
e = np.ones(n)
cplx = isinstance(s, complex) or np.iscomplexobj(s)
M = s * np.eye(n) - D0
val = alpha @ np.linalg.solve(M.astype(complex) if cplx else M, -(D0 @ e))
return complex(val) if cplx else float(np.real(val))
except (ValueError, NotImplementedError, ArithmeticError): # LinAlgError is a ValueError
pass
mean = self.getMean()
if not math.isfinite(mean) or mean <= 0.0:
return 1.0
n = 1000
dx = (20.0 * mean) / n
total = 0.0 + 0.0j if isinstance(s, complex) else 0.0
for i in range(1, n + 1):
x = i * dx
total += np.exp(-s * x) * self.evalPDF(x)
return (complex(total) if isinstance(s, complex) else float(total)) * dx
# CamelCase alias
evalLaplaceTransform = evalLST
[docs]
class DiscreteDistribution(Distribution):
"""Base class for discrete probability distributions."""
def __init__(self):
super().__init__()
[docs]
def isContinuous(self) -> bool:
return False
[docs]
def isDiscrete(self) -> bool:
return True
[docs]
def evalPMF(self, x: int) -> float:
"""Evaluate the probability mass function at point x."""
raise LineUnsupportedError("PMF not implemented for this distribution", error_id='LINE:unsupported:Analysis')
[docs]
class Markovian(Distribution):
"""
Base class for Markovian (phase-type) distributions.
Markovian distributions can be represented in terms of
initial probability vectors and transition rate matrices.
"""
def __init__(self):
super().__init__()
[docs]
def getD0(self) -> np.ndarray:
"""Get the D0 matrix (for MAP representations)."""
raise NotImplementedError
[docs]
def getD1(self) -> np.ndarray:
"""Get the D1 matrix (for MAP representations)."""
raise NotImplementedError
[docs]
def getMu(self) -> np.ndarray:
"""Get the service rates in each phase."""
raise NotImplementedError
[docs]
def getPhi(self) -> np.ndarray:
"""Get the completion probabilities from each phase."""
raise NotImplementedError
[docs]
def getInitProb(self) -> np.ndarray:
"""Get the initial probability vector."""
raise NotImplementedError
[docs]
def getNumberOfPhases(self) -> int:
"""Get the number of phases."""
raise NotImplementedError
[docs]
def getRepresentation(self) -> list:
"""Return the ``(D0, D1)`` matrix representation used in the theory of
Markovian arrival processes, as a list ``[D0, D1]`` (element ``k``
corresponds to matrix ``D_k``)."""
return [self.getD0(), self.getD1()]
[docs]
def plot(self, show: bool = True, ax=None):
"""Draw the phase structure (Markovian.plot): one node per phase, an edge per positive rate, labelled by
the rate; hidden transitions (D0) in black, observable ones (D1, a completion or arrival) in red.
Returns the matplotlib figure. Requires matplotlib and networkx.
"""
try:
import matplotlib.pyplot as plt
import networkx as nx
except ImportError:
raise ImportError("Matplotlib and NetworkX are required for plotting. "
"Install with: pip install matplotlib networkx")
rep = self.getRepresentation()
D0 = np.atleast_2d(np.asarray(rep[0], dtype=float))
D1 = np.atleast_2d(np.asarray(sum(np.asarray(m, dtype=float) for m in rep[1:]), dtype=float))
n = D0.shape[0]
G = nx.MultiDiGraph()
G.add_nodes_from(range(1, n + 1))
edges, colors, labels = [], [], {}
for i in range(n):
for j in range(n):
for M, col in ((D0, 'k'), (D1, 'r')):
if (col == 'r' or i != j) and M[i, j] > 0:
G.add_edge(i + 1, j + 1)
edges.append((i + 1, j + 1))
colors.append(col)
labels[(i + 1, j + 1)] = ('%s, %.4g' % (labels[(i + 1, j + 1)], M[i, j])
if (i + 1, j + 1) in labels else '%.4g' % M[i, j])
if ax is None:
fig, ax = plt.subplots()
else:
fig = ax.figure
pos = {k: (k - 1, 0.0) for k in range(1, n + 1)} if n > 0 else {}
nx.draw_networkx_nodes(G, pos, ax=ax, node_color='k', node_size=80)
nx.draw_networkx_labels(G, {k: (x, y + 0.08) for k, (x, y) in pos.items()}, ax=ax)
nx.draw_networkx_edges(G, pos, ax=ax, edgelist=edges, edge_color=colors, arrows=True,
connectionstyle='arc3,rad=0.25')
nx.draw_networkx_edge_labels(G, pos, edge_labels=labels, ax=ax, font_size=8)
ax.set_axis_off()
if show:
plt.show()
return fig
[docs]
def get_representation(self) -> list:
"""snake_case alias for :meth:`getRepresentation`."""
return self.getRepresentation()
[docs]
def evalLST(self, s):
"""
Laplace-Stieltjes transform, alpha (sI - D0)^-1 (-D0 e), as in MATLAB.
Distribution.evalLST is a rectangle rule over evalPDF: it carries
percent-level error on a phase-type law (Exp(1) at s = 0.5 read
0.65672 against 2/3) and returns ZERO for a COMPLEX argument, since it
evaluates math.exp. Every law reaching here is Markovian, so the closed
form applies, and being analytic it also serves the complex arguments
that transform inversion and root location need.
"""
D0 = np.atleast_2d(np.asarray(self.getD0(), dtype=float))
n = D0.shape[0]
e = np.ones(n)
t = -(D0 @ e)
alpha = np.asarray(self.getInitProb(), dtype=float).ravel()
cplx = isinstance(s, complex) or np.iscomplexobj(s)
M = (s * np.eye(n) - D0)
val = alpha @ np.linalg.solve(M.astype(complex) if cplx else M, t)
return complex(val) if cplx else float(np.real(val))
[docs]
def getPH(self) -> dict:
"""Return the phase-type representation as a dict {0: D0, 1: D1},
matching the JAR Map<Integer, Matrix> shape."""
return {0: self.getD0(), 1: self.getD1()}
get_ph = _late_alias('getPH')
# Aliases
get_d0 = getD0
get_d1 = getD1
get_mu = _late_alias('getMu')
get_phi = _late_alias('getPhi')
get_init_prob = _late_alias('getInitProb')
get_number_of_phases = _late_alias('getNumberOfPhases')