MarkovModelModule

StateTransitionGraph: A directed graph class for modeling continuous-time Markov chains (CTMCs).

This module defines the StateTransitionGraph class, an extension of networkx.DiGraph with specialized methods and attributes for analyzing, simulating, and visualizing continuous-time Markov chains (CTMCs).

Key Features

  • Symbolic and numerical parameter support using SymPy.
  • Automatic generation and storage of the rate matrix (Q-matrix).
  • Simulation of state transitions and computation of steady-state probabilities.
  • Enhanced graph visualization with transition rates and labels.
  • Support for dynamic import of constants via a custom module (default: "Constants").

Typical Use Case

  1. Define a set of model parameters (e.g., lambda, mu, n).
  2. Create a StateTransitionGraph instance with those parameters.
  3. Add states and transitions with symbolic rates.
  4. Compute the rate matrix and steady-state probabilities.
  5. Run simulations and visualize the CTMC behavior.

Example

>>> parameters = {'lambda': lam, 'mu': 1.0, 'n': 5}
>>> G = StateTransitionGraph(parameters)
>>> # Define transitions
>>> for i in range(parameters['n']):
>>>     G.addTransition(i, i + 1, G.sym("lambda"), tid=const.ARRIVAL)
>>>     G.addTransition(i + 1, i, (i + 1) * G.sym("mu"), tid=const.DEPARTURE)
>>> # Set node properties
>>> G.setAllStateDefaultProperties()
>>> G.setStateColor(G.n, const.COLOR_NODE_BLOCKING)
>>> # Compute steady-state probabilities
>>> probs = G.calcSteadyStateProbs()
>>> # Plot bar chart of state probabilities
>>> import matplotlib.pyplot as plt
>>> plt.figure(1, clear=True)
>>> plt.bar(probs.keys(), probs.values())
>>> plt.xlabel('states')
>>> plt.ylabel('probability')
>>> # Layout for transition graph
>>> pos = {s: (s, 0) for s in G.states()}
>>> G.drawTransitionGraph(pos=pos, bended=True, label_pos=0.5, num=2)
>>> # Simulate and animate system dynamics
>>> states, times = G.simulateMarkovModel(startNode=0, numEvents=2000)

See Also

networkx.DiGraph : The base class that is extended. sympy.symbols : Used to define symbolic transition rates. matplotlib.pyplot : Used internally for graph visualization.

Author

Tobias Hossfeld tobias.hossfeld@uni-wuerzburg.de

Created

5 April 2025

   1# -*- coding: utf-8 -*-
   2"""
   3StateTransitionGraph: A directed graph class for modeling continuous-time Markov chains (CTMCs).
   4
   5This module defines the `StateTransitionGraph` class, an extension of `networkx.DiGraph` 
   6with specialized methods and attributes for analyzing, simulating, and visualizing 
   7continuous-time Markov chains (CTMCs).
   8
   9Key Features
  10------------
  11- Symbolic and numerical parameter support using SymPy.
  12- Automatic generation and storage of the rate matrix (Q-matrix).
  13- Simulation of state transitions and computation of steady-state probabilities.
  14- Enhanced graph visualization with transition rates and labels.
  15- Support for dynamic import of constants via a custom module (default: "Constants").
  16
  17Typical Use Case
  18----------------
  191. Define a set of model parameters (e.g., lambda, mu, n).
  202. Create a `StateTransitionGraph` instance with those parameters.
  213. Add states and transitions with symbolic rates.
  224. Compute the rate matrix and steady-state probabilities.
  235. Run simulations and visualize the CTMC behavior.
  24
  25
  26Example
  27-------
  28>>> parameters = {'lambda': lam, 'mu': 1.0, 'n': 5}
  29>>> G = StateTransitionGraph(parameters)
  30
  31>>> # Define transitions
  32>>> for i in range(parameters['n']):
  33>>>     G.addTransition(i, i + 1, G.sym("lambda"), tid=const.ARRIVAL)
  34>>>     G.addTransition(i + 1, i, (i + 1) * G.sym("mu"), tid=const.DEPARTURE)
  35
  36>>> # Set node properties
  37>>> G.setAllStateDefaultProperties()
  38>>> G.setStateColor(G.n, const.COLOR_NODE_BLOCKING)
  39
  40>>> # Compute steady-state probabilities
  41>>> probs = G.calcSteadyStateProbs()
  42
  43>>> # Plot bar chart of state probabilities
  44>>> import matplotlib.pyplot as plt
  45>>> plt.figure(1, clear=True)
  46>>> plt.bar(probs.keys(), probs.values())
  47>>> plt.xlabel('states')
  48>>> plt.ylabel('probability')
  49
  50>>> # Layout for transition graph
  51>>> pos = {s: (s, 0) for s in G.states()}
  52>>> G.drawTransitionGraph(pos=pos, bended=True, label_pos=0.5, num=2)
  53
  54>>> # Simulate and animate system dynamics
  55>>> states, times = G.simulateMarkovModel(startNode=0, numEvents=2000)
  56
  57See Also
  58--------
  59networkx.DiGraph : The base class that is extended.
  60sympy.symbols   : Used to define symbolic transition rates.
  61matplotlib.pyplot : Used internally for graph visualization.
  62
  63
  64Author
  65------
  66Tobias Hossfeld <tobias.hossfeld@uni-wuerzburg.de>
  67
  68Created
  69-------
  705 April 2025 
  71
  72"""
  73
  74
  75import networkx as nx  # For the magic
  76import matplotlib.pyplot as plt  # For plotting
  77import numpy as np
  78from scipy import linalg
  79#import itertools
  80import sympy as sp
  81import warnings
  82import random
  83import importlib
  84from scipy.linalg import expm
  85from IPython.display import display, clear_output
  86from openpyxl import load_workbook
  87from openpyxl.styles import PatternFill
  88import matplotlib.colors as mcolors
  89import pandas as pd
  90import os
  91#import Constants as const
  92#from collections import namedtuple
  93#%%
  94class StateTransitionGraph(nx.DiGraph):
  95    """
  96    Initialize the StateTransitionGraph for continuous time Markov chains (CTMC). 
  97    Extends the networkx.Graph class to include additional methods and properties 
  98    for the analysis, visualization, and simulation of CTMCs. 
  99    The analysis allows deriving the state probabilities in the steady-state 
 100    and in the transient phase.
 101
 102    This constructor sets up symbolic and numeric parameters for the state transition graph,
 103    imports a constants module dynamically, and initializes internal data structures for
 104    rate matrix computation and steady-state probability analysis.
 105
 106    Parameters
 107    ----------
 108    - **parameters** : dict  
 109      A dictionary mapping parameter names to their numerical values.     
 110        Example: {'lambda': 1.0, 'mu': 1.0}. 
 111        For each parameter key:
 112        - A symbolic variable is created using sympy.symbols. Can be accessed as `G.sym("key")` or `G.sym_key`
 113        - A variable is created which can be accessed via `G.key` (and the numerical value is assigned)
 114        - For the key="lambda", `G.lam` and `G.sym_lambda` are created
 115    
 116        The parameters can be accessed via property `parameters`
 117    - **constants** : str, optional
 118        The name of a Python module (as string) to import dynamically and assign to `self.const`.
 119        This can be used to store global or user-defined constants.
 120        Default is "Constants".
 121    - `*args` : tuple
 122        Additional positional arguments passed to the base `networkx.DiGraph` constructor.
 123    - `**kwargs` : dict
 124        Additional keyword arguments passed to the base `networkx.DiGraph` constructor.
 125
 126    Attributes
 127    ----------
 128    - `_rate_matrix` : ndarray or None  
 129      The continuous-time rate matrix `Q` of the CTMC. Computed when calling `createRateMatrix()`.
 130    - `_state_probabilities` : dict or None  
 131      Dictionary of steady-state probabilities, computed via `calcSteadyStateProbs()`.
 132    - `_n2i` : dict or None  
 133      Mapping from node identifiers to matrix row/column indices used for rate matrix construction.
 134    - `_sym` : dict  
 135      Mapping of parameter names (strings) to their symbolic `sympy` representations.
 136    - `_subs` : dict  
 137      Mapping from symbolic parameters to their numerical values (used for substitutions).
 138    - `_parameters` : dict  
 139      Copy of the original parameter dictionary passed at initialization.
 140    - `_transitionProbabilitesComputed` : bool  
 141      Internal flag indicating whether transition probabilities for the embedded Markov chain
 142      have been computed.
 143    - `const` : module  
 144      Dynamically imported module (default `"Constants"`), used for storing user-defined constants 
 145      such as colors or transition IDs.
 146        
 147    Inherits
 148    --------
 149        networkx.Graph: The base graph class from NetworkX.
 150    """
 151
 152    def __init__(self, parameters, constants="Constants", *args, **kwargs):
 153        """
 154        Initialize the state transition graph with symbolic parameters.
 155
 156        This constructor sets up the symbolic and numerical environment for modeling
 157        a CTMC. Parameters are stored as both symbolic expressions and numerical values,
 158        and convenient attribute access is provided for both forms.
 159
 160        Parameters
 161        ----------
 162        - **parameters** : dict  
 163          A dictionary mapping parameter names to their numerical values.     
 164            Example: {'lambda': 1.0, 'mu': 1.0}. 
 165            For each parameter key:
 166            - A symbolic variable is created using sympy.symbols. Can be accessed as `G.sym("key")` or `G.sym_key`
 167            - A variable is created which can be accessed via `G.key` (and the numerical value is assigned)
 168            - For the key="lambda", `G.lam` and `G.sym_lambda` are created
 169        
 170            The parameters can be accessed via property `parameters`
 171        - **constants** : str, optional
 172            The name of a Python module (as string) to import dynamically and assign to `self.const`.
 173            This can be used to store global or user-defined constants.
 174            Default is "Constants".
 175        - `*args` : tuple
 176            Additional positional arguments passed to the base `networkx.DiGraph` constructor.
 177        - `**kwargs` : dict
 178            Additional keyword arguments passed to the base `networkx.DiGraph` constructor.
 179        """        
 180        super().__init__(*args, **kwargs)
 181        self._rate_matrix = None  # Rate Matrix Q of this transition graph
 182        self._state_probabilities = None  # steady state proabilities of this CTMC
 183        self._n2i = None # internal index for generating the rate matrix Q: dictionary (node to index)
 184        
 185        #parameters = {'lambda': 1.0, 'mu':1.0, 'n':5}
 186        self._parameters = parameters
 187        self._sym = {key: sp.symbols(key) for key in parameters}
 188        self._subs = {} # contains symbolic parameter as key and provides numerical value for computation.
 189        for k in parameters:
 190            self._subs[ self._sym[k] ] = parameters[k] 
 191
 192        # access to parameters and symbols using dot-notation, e.g. G.lam, G.sym_lambda, or G.sym("lambda")
 193        for key, value in parameters.items():
 194            if key=='lambda':
 195                setattr(self, 'lam', value)
 196            else:
 197                setattr(self, key, value)
 198            
 199            setattr(self, "sym_"+key, self._sym[key])
 200            
 201        self._transitionProbabilitesComputed = False
 202        
 203        self._const = importlib.import_module(constants)
 204
 205        
 206        
 207    def addState(self, state, color=None, label=None):
 208        """
 209        Add a state (node) to the graph with optional color and label attributes.
 210    
 211        If the node already exists, its attributes are updated with the provided values.
 212    
 213        Parameters
 214        ----------
 215        - **state** : hashable
 216            The identifier for the state. Can be a string, int, or a tuple. If no label
 217            is provided, the label will be auto-generated as:
 218            - comma-separated string for tuples (e.g., `(1, 2)` -> `'1,2'`)
 219            - stringified value for other types
 220    
 221        - **color** : str, optional
 222            The color assigned to the state (used for visualization).
 223            If not specified, defaults to `self.const.COLOR_NODE_DEFAULT`.
 224    
 225        - **label** : str, optional
 226            A label for the node (used for display or annotation). If not provided,
 227            it will be auto-generated from the state identifier.
 228    
 229        Notes
 230        -----
 231        - This method wraps `self.add_node()` from NetworkX.
 232        - If the node already exists in the graph, this method does not raise an error;
 233          it will update the node’s existing attributes.
 234        """
 235        if label is None:
 236            if isinstance(state, tuple):
 237                label = ','.join(map(str, state)) # if state is a tuple
 238            else: 
 239                label = str(state)
 240            #self.nodes[state]["label"] = label
 241        
 242        if color is None:
 243            color = self._const.COLOR_NODE_DEFAULT
 244        
 245        # no error if the node already exists.Attributes get updated if you pass them again.
 246        self.add_node(state, color=color, label=label)
 247        
 248    def setStateProperties(self, state, color=None, label=None, **kwargs):
 249        """
 250        Update visual or metadata properties for an existing state (node) in the graph.
 251    
 252        This method updates the node's attributes such as color, label, and any additional
 253        custom key-value pairs. The node must already exist in the graph.
 254    
 255        Parameters
 256        ----------
 257        - **state** : hashable
 258            Identifier of the state to update. Must be a node already present in the graph.
 259    
 260        - **color** : str, optional
 261            Color assigned to the node, typically used for visualization. If not provided,
 262            the default node color (`self._const.COLOR_NODE_DEFAULT`) is used.
 263    
 264        - **label** : str, optional
 265            Label for the node. If not provided, it is auto-generated from the state ID
 266            (especially for tuples).
 267    
 268        - `**kwargs` : dict
 269            Arbitrary additional key-value pairs to set as node attributes. These can include
 270            any custom metadata.
 271    
 272        Raises
 273        ------
 274        KeyError
 275            If the specified state is not present in the graph.
 276    
 277        Notes
 278        -----
 279        - This method ensures that required attributes (like `color` and `label`) are set for visualization. 
 280        - Existing node attributes are updated or extended with `kwargs`.
 281        """        
 282        if state not in self:
 283            raise KeyError(f"State '{state}' does not exist in the graph.")
 284        
 285        # mandatory for plotting
 286        self.addState(state, color, label)
 287        
 288        # # Add any additional attributes
 289        self.nodes[state].update(kwargs)
 290        
 291    def setAllStateDefaultProperties(self, **kwargs):
 292        """
 293        Set default visual and metadata properties for all states (nodes) in the graph.
 294    
 295        This method applies `setStateProperties()` to every node, ensuring each state
 296        has at least the required default attributes (e.g., color and label).
 297        Any additional keyword arguments are passed to `setStateProperties()`.
 298    
 299        Parameters
 300        ----------
 301        `**kwargs` : dict, optional
 302            Optional keyword arguments to apply uniformly to all nodes,
 303            such as color='gray', group='core', etc.
 304    
 305        Notes
 306        -----
 307        - This method is useful for initializing node attributes before plotting
 308          or performing other graph-based computations.
 309        - Nodes without existing attributes will receive defaults via `addState()`.
 310        """
 311        for state in self.nodes():
 312            self.setStateProperties(state, **kwargs)
 313
 314        
 315    #def addTransition(self, origin_state, list_destination_state, sym_rate, color=plt.cm.tab10(5), tid = None):        
 316    def addTransition(self, origin_state, destination_state, sym_rate, tid = None, color=None, label=None):     
 317        """
 318        Add a directed transition (edge) between two states with a symbolic transition rate.
 319    
 320        This method adds a directed edge from `origin_state` to `destination_state`,
 321        stores both the symbolic and numerical rate, and sets optional display properties
 322        such as color, label, and transition ID. If the edge already exists, a warning is issued.
 323    
 324        Parameters
 325        ----------
 326        - **origin_state** : hashable
 327            The source node (state) of the transition.
 328    
 329        - **destination_state** : hashable
 330            The target node (state) of the transition.
 331    
 332        - **sym_rate** : sympy.Expr
 333            A symbolic expression representing the transition rate.
 334            This will be numerically evaluated using internal substitutions (`self._subs`).
 335    
 336        - **tid** : str or int, optional
 337            A transition identifier (e.g. for grouping or styling).
 338            If provided, the color will default to `self._const.COLOR_TRANSITIONS[tid]` if not set explicitly.
 339    
 340        - **color** : str, optional
 341            Color for the edge (used in visualization). Defaults to the transition-specific color
 342            if `tid` is provided.
 343    
 344        - **label** : str, optional
 345            LaTeX-formatted string for labeling the edge. If not provided, a default label is
 346            generated using `sympy.latex(sym_rate)`.
 347    
 348        Notes
 349        -----
 350        - The numeric rate is computed by substituting known parameter values into `sym_rate`.
 351        - The following attributes are stored on the edge:
 352            - `rate` (float): numeric value of the transition rate
 353            - `sym_rate` (sympy.Expr): original symbolic expression
 354            - `label` (str): display label
 355            - `color` (str): edge color
 356            - `tid` (str or int): transition ID
 357        - If the edge already exists, a `UserWarning` is issued but the edge is still overwritten.
 358        """        
 359        if tid is not None:
 360            if color is None:
 361                color = self._const.COLOR_TRANSITIONS[ tid ]
 362                
 363        rate = sym_rate.subs( self._subs )
 364        
 365        if label is None: string = "$"+sp.latex(sym_rate)+"$"
 366                
 367        if self.has_edge(origin_state, destination_state):
 368            theedge = self[origin_state][destination_state]
 369            warningstr = f"Edge already exists: {origin_state} -> {destination_state} with rate {theedge['label']}"
 370            warnings.warn(warningstr, category=UserWarning)
 371        self.add_edge(origin_state, destination_state, rate=float(rate), label=string, 
 372                   color=color, sym_rate=sym_rate, tid=tid)
 373    
 374    def sym(self, key):
 375        """
 376        Retrieve the symbolic representation of a model parameter.
 377    
 378        Parameters
 379        ----------
 380        **key** : str
 381            The name of the parameter (e.g., 'lambda', 'mu', 'n').
 382    
 383        Returns
 384        -------
 385        **sympy.Symbol**
 386            The symbolic variable corresponding to the given key.
 387    
 388        Raises
 389        ------
 390        KeyError
 391            If the key does not exist in the internal symbolic mapping (`_sym`).
 392    
 393        Examples
 394        --------
 395        >>> G.sym('lambda')
 396        lambda
 397    
 398        Notes
 399        -----
 400        - The symbolic mapping is initialized during object construction based on the `parameters` dictionary.
 401        - To access the symbol using dot-notation, use `.sym_key` attributes created during initialization.
 402        """
 403        return self._sym[key]
 404
 405        
 406    @property
 407    def parameters(self):
 408        """
 409        Return the model parameters as a dictionary.
 410    
 411        This includes the original numerical values provided during initialization.
 412    
 413        Returns
 414        -------
 415        **dict**
 416            Dictionary of parameter names and their corresponding numeric values.
 417    
 418        Examples
 419        --------
 420        >>> G.parameters
 421        {'lambda': 1.0, 'mu': 0.5, 'n': 3}
 422    
 423        Notes
 424        -----
 425        - This property does not return symbolic values. Use `self.sym(key)` for symbolic access.
 426        - Parameter symbols can also be accessed using attributes like `sym_lambda` or the `sym()` method.
 427        """
 428        return self._parameters
 429
 430    
 431
 432    @property
 433    def rate_matrix(self):
 434        """
 435        Return the continuous-time rate matrix (Q matrix) of the state transition graph.
 436    
 437        This matrix defines the transition rates between all pairs of states in the
 438        continuous-time Markov chain (CTMC) represented by this graph.
 439    
 440        Returns
 441        -------
 442        **numpy.ndarray** or None
 443            The rate matrix `Q`, where `Q[i, j]` is the transition rate from state `i` to `j`.
 444            Returns `None` if the rate matrix has not been computed yet.
 445    
 446        Notes
 447        -----
 448        - The matrix is typically generated via a method like `computeRateMatrix()` or similar.
 449        - Internal state indexing is handled by `_n2i` (node-to-index mapping).
 450        - The rate matrix is stored internally as `_rate_matrix`.
 451        """
 452        return self._rate_matrix
 453
 454
 455
 456    @property
 457    def state_probabilities(self):
 458        """
 459        Return the steady-state probabilities of the continuous-time Markov chain (CTMC).
 460    
 461        These probabilities represent the long-run proportion of time the system spends
 462        in each state, assuming the CTMC is irreducible and positive recurrent.
 463    
 464        Returns
 465        -------
 466        **dict** or None
 467            A dictionary mapping each state to its steady-state probability.
 468            Returns `None` if the probabilities have not yet been computed.
 469    
 470        Notes
 471        -----
 472        - The probabilities are stored internally after being computed via a dedicated method
 473          `calcSteadyStateProbs()`.
 474        - The result is stored in the `_state_probabilities` attribute.
 475        """
 476        return self._state_probabilities
 477
 478     
 479    
 480    def printEdges(self, state):
 481        """
 482        Print all outgoing edges from a given state in the graph.
 483    
 484        For each edge from the specified `state`, this method prints the destination
 485        state and the corresponding transition label (typically the symbolic rate).
 486    
 487        Parameters
 488        ----------
 489        **state** : hashable
 490            The source node (state) from which outgoing edges will be printed.            
 491    
 492        Returns
 493        -------
 494        None
 495            This method prints output directly to the console and does not return a value.
 496    
 497        Notes
 498        -----
 499        - The transition label is retrieved from the edge attribute `'label'`.
 500        - Assumes `state` exists in the graph; otherwise, a KeyError will be raised.
 501        """
 502        print(f"from {state}")
 503        node = self[state]
 504        for e in node:
 505            rate = node[e]['label']
 506            print(f" -> {e} with {rate}")
 507
 508            
 509    def createRateMatrix(self):
 510        """
 511        Construct the continuous-time Markov chain (CTMC) rate matrix Q from the graph.
 512    
 513        This method generates the generator matrix `Q` corresponding to the structure and
 514        rates defined in the directed graph. Each node represents a system state, and each
 515        directed edge must have a `'rate'` attribute representing the transition rate from
 516        one state to another.
 517    
 518        The matrix `Q` is constructed such that:
 519        - `Q[i, j]` is the transition rate from state `i` to state `j`
 520        - Diagonal entries satisfy `Q[i, i] = -sum(Q[i, j]) for j ≠ i`, so that each row sums to 0
 521    
 522        The resulting matrix and the internal node-to-index mapping (`_n2i`) are stored
 523        as attributes of the graph.
 524    
 525        Returns
 526        -------
 527        None
 528            The method updates internal attributes:
 529            - `_rate_matrix` : numpy.ndarray
 530                The CTMC rate matrix Q
 531            - `_n2i` : dict
 532                Mapping from node identifiers to matrix row/column indices
 533    
 534        Notes
 535        -----
 536        - Assumes all edges in the graph have a `'rate'` attribute.
 537        - The node order used in the matrix is consistent with the order from `list(self.nodes)`.
 538        - Use `self.rate_matrix` to access the generated matrix after this method is called.
 539        """
 540        n2i = {}
 541        nodes = list(self.nodes)
 542        for i,node in enumerate(nodes):
 543            n2i[node] = i         
 544        Q = np.zeros((len(nodes),len(nodes))) # rate matrix
 545        
 546        for edge in self.edges:
 547            i0 = n2i[edge[0]]
 548            i1 = n2i[edge[1]]
 549            Q[i0,i1] = self[edge[0]][edge[1]]["rate"] 
 550            
 551        np.fill_diagonal(Q, -Q.sum(axis=1))
 552        self._rate_matrix = Q
 553        self._n2i = n2i
 554        #return Q, n2i
 555    
 556    def calcSteadyStateProbs(self, verbose=False):        
 557        """
 558        Compute the steady-state probabilities of the continuous-time Markov chain (CTMC).
 559    
 560        This method calculates the stationary distribution of the CTMC defined by the rate
 561        matrix of the state transition graph. The computation solves the linear system
 562        `XQ = 0` with the constraint that the probabilities X sum to 1, using a modified
 563        version of the rate matrix `Q`.
 564    
 565        Parameters
 566        ----------
 567        `verbose` : bool, optional
 568            If True, prints the rate matrix, the modified matrix, the right-hand side vector,
 569            and the resulting steady-state probabilities.
 570    
 571        Returns
 572        -------
 573        **probs** : dict
 574            A dictionary mapping each state (node in the graph) to its steady-state probability.
 575    
 576        Notes
 577        -----
 578        - The graph must already contain valid transition rates on its edges.
 579        - If the rate matrix has not yet been created, `createRateMatrix()` is called automatically.
 580        - The method modifies the last column of the rate matrix to impose the normalization constraint
 581          `sum(X) = 1`.
 582        - The inverse of the modified matrix is computed directly, so the graph should be small enough
 583          to avoid numerical instability or performance issues.
 584        - The resulting probabilities are stored internally in `_state_probabilities` and can be accessed
 585          via the `state_probabilities` property.
 586        """
 587        
 588        # compute transition matrix if not already done
 589        if self._rate_matrix is None:  self.createRateMatrix()
 590        Q2 =  self._rate_matrix.copy()
 591        if verbose:
 592            print("Rate matrix Q")
 593            print(f'Q=\n{Q2}\n')
 594        
 595        Q2[:, -1] = 1
 596        if verbose:        
 597            print(f'Matrix is changed to\nQ2=\n{Q2}\n')
 598
 599        b = np.zeros(len(Q2))
 600        b[-1] = 1
 601        if verbose:
 602            print("Solve X = b @ Q2")
 603            print(f'b=\n{b}\n')
 604        
 605        # state probabilities
 606        X = b @ linalg.inv(Q2) # compute the matrix inverse
 607            
 608        # Generate a matrix with P(X,Y)
 609        matrix = {}
 610        for n in self.nodes:
 611            matrix[n] = X[ self._n2i[n] ]
 612        
 613        if verbose:
 614            print("Steady-state probabilities X")
 615            print(f'X=\n{matrix}\n')
 616        self._state_probabilities = matrix
 617        return matrix
 618        
 619    def calcTransitionProbabilities(self):
 620        """
 621        Compute and store transition probabilities for each node in the graph of the embedded discrete Markov chain.
 622    
 623        This method processes each node in the graph, treating outgoing edges as
 624        transitions in an embedded discrete-time Markov chain derived from a continuous-time
 625        Markov chain (CTMC). Transition *rates* on outgoing edges are normalized into
 626        *probabilities*, and relevant data is stored in the node attributes.
 627    
 628        For each node, the following attributes are set:
 629        - `"transitionProbs"` : np.ndarray
 630            Array of transition probabilities to successor nodes.
 631        - `"transitionNodes"` : list
 632            List of successor node identifiers (ordered as in the probability array).
 633        - `"sumOutgoingRates"` : float
 634            Sum of all outgoing transition rates. Used to model the exponential waiting time
 635            in continuous-time simulation (`T ~ Exp(sumOutgoingRates)`).
 636    
 637        Notes
 638        -----
 639        - The graph must be a directed graph where each edge has a `'rate'` attribute.
 640        - Nodes with no outgoing edges are considered absorbing states.
 641        - Useful for discrete-event simulation and Markov chain Monte Carlo (MCMC) modeling.
 642        - Sets an internal flag `_transitionProbabilitesComputed = True` upon completion.
 643    
 644        Raises
 645        ------
 646        KeyError
 647            If any outgoing edge lacks a `'rate'` attribute.
 648        """
 649        for node in self:
 650            successors = list(self.successors(node))
 651            rates = np.array([self[node][nbr]['rate'] for nbr in successors])
 652            probs = rates / rates.sum()  # Normalize to probabilities
 653            self.nodes[node]["transitionProbs"] = probs # probabilities
 654            self.nodes[node]["transitionNodes"] = successors # list of nodes
 655            self.nodes[node]["sumOutgoingRates"] = rates.sum() # time in state ~ EXP(sumOutgoingRates)
 656            numel = len(successors)
 657            if numel==0: 
 658                print(f"Absorbing state: {node}")
 659        self._transitionProbabilitesComputed = True
 660                
 661    def simulateMarkovModel(self, startNode=None, numEvents=10):
 662        """
 663        Simulate a trajectory through the state transition graph as a continuous-time Markov chain (CTMC).
 664        
 665        This method performs a discrete-event simulation of the CTMC by generating a random walk
 666        through the graph, starting at `startNode` (or a default node if not specified). Transitions
 667        are chosen based on the normalized transition probabilities computed from edge rates, and
 668        dwell times are sampled from exponential distributions based on the total outgoing rate
 669        from each state.
 670        
 671        Parameters
 672        ----------
 673        - **startNode** : hashable, optional
 674            The node at which the simulation starts. If not specified, the first node in the graph
 675            (based on insertion order) is used.
 676        
 677        - **numEvents** : int, optional
 678            The number of transition events (steps) to simulate. Default is 10.
 679        
 680        Returns
 681        -------
 682        **states** : list
 683            A list of visited states, representing the sequence of nodes traversed in the simulation.
 684        
 685        **times** : numpy.ndarray
 686            An array of dwell times in each state, sampled from exponential distributions
 687            with rate equal to the sum of outgoing transition rates from each state.
 688        
 689        Notes
 690        -----
 691        - This method automatically calls `calcTransitionProbabilities()` if transition probabilities
 692          have not yet been computed.
 693        - If an absorbing state (i.e., a node with no outgoing edges) is reached, the simulation stops early.
 694        - The cumulative time can be obtained by `np.cumsum(times)`.
 695        - State durations follow `Exp(sumOutgoingRates)` distributions per CTMC theory.
 696        """
 697        if not self._transitionProbabilitesComputed: self.calcTransitionProbabilities()
 698        
 699        # simulate random walk through Graph = Markov Chain
 700        if startNode is None:
 701            startNode = next(iter(self))
 702        
 703
 704        states = [] # list of states
 705        times = np.zeros(numEvents) # time per state
 706
 707        node = startNode
 708        for i in range(numEvents):
 709            outgoingRate = self.nodes[node]["sumOutgoingRates"]  # time in state ~ EXP(sumOutgoingRates)
 710            times[i] = random.expovariate(outgoingRate) # exponentially distributed time in state
 711            states.append(node)
 712            
 713            # get next node    
 714            probs = self.nodes[node]["transitionProbs"] # probabilities
 715            successors = self.nodes[node]["transitionNodes"] # list of nodes        
 716            numel = len(successors)
 717            if numel==0: 
 718                print(f"Absorbing state: {node}")
 719                break
 720            elif numel==1:
 721                nextNode = successors[0]
 722            else:
 723                nextNode = successors[ np.random.choice(numel, p=probs) ]
 724            # edge = G[node][nextNode]
 725            node = nextNode
 726            
 727        return states, times
 728    
 729    def calcProbabilitiesSimulationRun(self, states, times):
 730        """
 731        Estimate empirical state probabilities from a simulation run of the CTMC.
 732    
 733        This method computes the fraction of total time spent in each state during a simulation.
 734        It uses the sequence of visited states and corresponding sojourn times to estimate the
 735        empirical (time-weighted) probability distribution over states.
 736    
 737        Parameters
 738        ----------
 739        - **states** : list of tuple
 740            A list of state identifiers visited during the simulation.    
 741        - **times** : numpy.ndarray
 742            An array of sojourn times corresponding to each visited state, same length as `states`.
 743    
 744        Returns
 745        -------
 746        **simProbs** : dict
 747            A dictionary mapping each unique state to its estimated empirical probability,
 748            based on total time spent in that state relative to the simulation duration.
 749    
 750        Notes
 751        -----
 752        - Internally, states are mapped to flat indices for efficient counting using NumPy.
 753        - This method assumes that `states` contains only valid node identifiers from the graph.
 754        - If `_n2i` (node-to-index mapping) has not been initialized, it will be computed here.
 755        - This is especially useful for validating analytical steady-state probabilities
 756          against sampled simulation results. Or in case of state-space explosions.
 757        """        
 758        if self._n2i is None:
 759            n2i = {}
 760            nodes = list(self.nodes)
 761            for i,node in enumerate(nodes):
 762                n2i[node] = i     
 763            self._n2i = n2i
 764                
 765        
 766        # Map tuples to flat indices
 767        data = np.array(states)
 768        max_vals = data.max(axis=0) + 1
 769        indices = np.ravel_multi_index(data.T, dims=max_vals)
 770
 771        # Weighted bincount
 772        bincount = np.bincount(indices, weights=times)
 773
 774        # Decode back to tuples
 775        tuple_indices = np.array(np.unravel_index(np.nonzero(bincount)[0], shape=max_vals)).T
 776        timePerState = bincount[bincount > 0]
 777
 778        simProbs = {}
 779        totalTime = times.sum()
 780        for t, c in zip(tuple_indices, timePerState):
 781            simProbs[tuple(t)] = c / totalTime
 782            #print(f"Tuple {tuple(t)} → weighted count: {c}")
 783        return simProbs
 784    
 785    
 786    def _draw_networkx_edge_labels(
 787        self,
 788        pos,
 789        edge_labels=None, label_pos=0.75,
 790        font_size=10, font_color="k", font_family="sans-serif", font_weight="normal",
 791        alpha=None, bbox=None,
 792        horizontalalignment="center", verticalalignment="center",
 793        ax=None, rotate=True, clip_on=True, rad=0):
 794        """
 795        Draw edge labels for a directed graph with optional curved edge support.
 796    
 797        This method extends NetworkX's default edge label drawing by supporting
 798        curved (bended) edges through quadratic Bézier control points. It places
 799        labels dynamically at a specified position along the curve and optionally
 800        rotates them to follow the edge orientation.
 801    
 802        Parameters
 803        ----------
 804        pos : dict
 805            Dictionary mapping nodes to positions (2D coordinates). Required.
 806    
 807        edge_labels : dict, optional
 808            Dictionary mapping edge tuples `(u, v)` to label strings. If None,
 809            edge data is used and converted to strings automatically.
 810    
 811        label_pos : float, optional
 812            Position along the edge to place the label (0=head, 1=tail).
 813            Default is 0.75 (closer to the target node).
 814    
 815        font_size : int, optional
 816            Font size of the edge labels. Default is 10.
 817    
 818        font_color : str, optional
 819            Font color for the edge labels. Default is 'k' (black).
 820    
 821        font_family : str, optional
 822            Font family for the edge labels. Default is 'sans-serif'.
 823    
 824        font_weight : str, optional
 825            Font weight for the edge labels (e.g., 'normal', 'bold'). Default is 'normal'.
 826    
 827        alpha : float, optional
 828            Opacity for the edge labels. If None, labels are fully opaque.
 829    
 830        bbox : dict, optional
 831            Matplotlib-style bbox dictionary to style the label background.
 832            Default is a white rounded box.
 833    
 834        horizontalalignment : str, optional
 835            Horizontal alignment of the text. Default is 'center'.
 836    
 837        verticalalignment : str, optional
 838            Vertical alignment of the text. Default is 'center'.
 839    
 840        ax : matplotlib.axes.Axes, optional
 841            Axes on which to draw the labels. If None, uses the current axes.
 842    
 843        rotate : bool, optional
 844            Whether to rotate the label to match the edge direction. Default is True.
 845    
 846        clip_on : bool, optional
 847            Whether to clip labels at the plot boundary. Default is True.
 848    
 849        rad : float, optional
 850            Curvature of the edge (used to determine Bézier control point).
 851            Positive values bend counterclockwise; negative values clockwise. Default is 0 (straight).
 852    
 853        Returns
 854        -------
 855        dict
 856            Dictionary mapping edge tuples `(u, v)` to the corresponding matplotlib Text objects.
 857    
 858        Notes
 859        -----
 860        - If `rotate=True`, label text is automatically aligned with the direction of the edge.
 861        - Curvature (`rad`) enables visualization of bidirectional transitions using bent edges.
 862        - This method is intended for internal use and supports bended edge label placement to match custom edge rendering.
 863    
 864        See Also
 865        --------
 866        draw
 867        draw_networkx
 868        draw_networkx_nodes
 869        draw_networkx_edges
 870        draw_networkx_labels
 871        """                
 872        if ax is None:
 873            ax = plt.gca()
 874        if edge_labels is None:
 875            labels = {(u, v): d for u, v, d in self.edges(data=True)}
 876        else:
 877            labels = edge_labels
 878        text_items = {}
 879        for (n1, n2), label in labels.items():
 880            (x1, y1) = pos[n1]
 881            (x2, y2) = pos[n2]
 882            (x, y) = (
 883                x1 * label_pos + x2 * (1.0 - label_pos),
 884                y1 * label_pos + y2 * (1.0 - label_pos),
 885            )
 886            pos_1 = ax.transData.transform(np.array(pos[n1]))
 887            pos_2 = ax.transData.transform(np.array(pos[n2]))
 888            #linear_mid = 0.5*pos_1 + 0.5*pos_2
 889            linear_mid = (1-label_pos)*pos_1 + label_pos*pos_2
 890            d_pos = pos_2 - pos_1
 891            rotation_matrix = np.array([(0,1), (-1,0)])
 892            ctrl_1 = linear_mid + rad*rotation_matrix@d_pos
 893            ctrl_mid_1 = 0.5*pos_1 + 0.5*ctrl_1
 894            ctrl_mid_2 = 0.5*pos_2 + 0.5*ctrl_1
 895            bezier_mid = 0.5*ctrl_mid_1 + 0.5*ctrl_mid_2
 896            (x, y) = ax.transData.inverted().transform(bezier_mid)
 897
 898            if rotate:
 899                # in degrees
 900                angle = np.arctan2(y2 - y1, x2 - x1) / (2.0 * np.pi) * 360
 901                # make label orientation "right-side-up"
 902                if angle > 90:
 903                    angle -= 180
 904                if angle < -90:
 905                    angle += 180
 906                # transform data coordinate angle to screen coordinate angle
 907                xy = np.array((x, y))
 908                trans_angle = ax.transData.transform_angles(
 909                    np.array((angle,)), xy.reshape((1, 2))
 910                )[0]
 911            else:
 912                trans_angle = 0.0
 913            # use default box of white with white border
 914            if bbox is None:
 915                bbox = dict(boxstyle="round", ec=(1.0, 1.0, 1.0), fc=(1.0, 1.0, 1.0))
 916            if not isinstance(label, str):
 917                label = str(label)  # this makes "1" and 1 labeled the same
 918
 919            t = ax.text(
 920                x, y, label,
 921                size=font_size, color=font_color, family=font_family, weight=font_weight,
 922                alpha=alpha,
 923                horizontalalignment=horizontalalignment, verticalalignment=verticalalignment,
 924                rotation=trans_angle, transform=ax.transData,
 925                bbox=bbox, zorder=1,clip_on=clip_on)
 926            text_items[(n1, n2)] = t
 927
 928        ax.tick_params(axis="both", which="both",bottom=False, left=False, labelbottom=False, labelleft=False)
 929
 930        return text_items    
 931
 932    
 933    def drawTransitionGraph(self, pos=None, 
 934                        bended=False, node_size=1000, num=1, rad=-0.2,
 935                        edge_labels=None, node_shapes='o', figsize=(14, 7),
 936                        clear=True, fontsize=10,
 937                        fontcolor='black', label_pos=0.75):        
 938        """
 939        Visualize the state transition graph with nodes and directed edges.
 940    
 941        This method draws the graph using `matplotlib` and `networkx`, including labels,
 942        node colors, and optional curved (bended) edges for better visualization of
 943        bidirectional transitions. It supports layout customization and is useful for
 944        understanding the structure of the Markov model.
 945    
 946        Parameters
 947        ----------
 948        - pos : dict, optional
 949            A dictionary mapping nodes to positions. If None, a Kamada-Kawai layout is used.
 950    
 951        - bended : bool, optional
 952            If True, edges are drawn with curvature using arcs. Useful for distinguishing
 953            bidirectional transitions. Default is False.
 954    
 955        - node_size : int, optional
 956            Size of the nodes in the plot. Default is 1000.
 957    
 958        - num : int, optional
 959            Figure number for matplotlib (useful for managing multiple figures). Default is 1.
 960    
 961        - rad : float, optional
 962            The curvature radius for bended edges. Only used if `bended=True`. Default is -0.2.
 963    
 964        - edge_labels : dict, optional
 965            Optional dictionary mapping edges `(u, v)` to labels. If None, the `'label'`
 966            attribute from each edge is used.
 967    
 968        - node_shapes : str, optional
 969            Shape of the nodes (e.g., `'o'` for circle, `'s'` for square). Default is `'o'`.
 970    
 971        - figsize : tuple, optional
 972            Size of the matplotlib figure. Default is (14, 7).
 973    
 974        - clear : bool, optional
 975            If True, clears the figure before plotting. Default is True.
 976    
 977        - fontsize : int, optional
 978            Font size used for node labels. Default is 10.
 979    
 980        - fontcolor : str, optional
 981            Font color for node labels. Default is `'black'`.
 982    
 983        - label_pos : float, optional
 984            Position of edge labels along the edge (0=start, 1=end). Default is 0.75.
 985    
 986        Returns
 987        -------
 988        - fig : matplotlib.figure.Figure
 989            The created matplotlib figure object.
 990    
 991        - ax : matplotlib.axes.Axes
 992            The axes object where the graph is drawn.
 993    
 994        - pos : dict
 995            The dictionary of node positions used for drawing.
 996    
 997        Notes
 998        -----
 999        - Node labels and colors must be set beforehand using `addState()` or `setStateProperties()`.
1000        - Edge attributes must include `'label'` and `'color'`.
1001        - This method modifies the current matplotlib figure and is intended for interactive or inline use.
1002        """        
1003        node_labels = {node: data["label"] for node, data in self.nodes(data=True)}        
1004        node_colors = {node: data["color"] for node, data in self.nodes(data=True)}        
1005        node_colors = [data["color"] for node, data in self.nodes(data=True)]
1006        #node_shapes = {(hw,sw): "o" if hw+sw<G.N else "s" for hw,sw in G.nodes()}        
1007            
1008        edge_cols = {(u,v): self[u][v]["color"] for u,v in self.edges}
1009            
1010        if edge_labels is None:
1011            edge_labels = {(u,v): self[u][v]["label"] for u,v in self.edges}
1012            
1013        if pos is None:
1014            pos = nx.kamada_kawai_layout(self)    
1015            
1016        plt.figure(figsize=figsize, num=num, clear=clear)                        
1017        nx.draw_networkx_nodes(self, pos, node_color=node_colors, node_shape = node_shapes, node_size=node_size)                    
1018        
1019        nx.draw_networkx_labels(self, pos, labels=node_labels, font_size=fontsize, font_color=fontcolor)  # Draw node labels
1020        if bended: 
1021            nx.draw_networkx_edges(self, pos,  width=1, edge_color=[edge_cols[edge] for edge in self.edges], 
1022                                   node_size = node_size, 
1023                                   arrows = True, arrowstyle = '-|>',
1024                                   connectionstyle=f"arc3,rad={rad}")        
1025            self._draw_networkx_edge_labels(pos, ax=plt.gca(), edge_labels=edge_labels,rotate=False, rad = rad, label_pos=label_pos)
1026        else:
1027            nx.draw_networkx_edges(self, pos,  width=1, edge_color=[edge_cols[edge] for edge in self.edges], 
1028                               node_size = node_size, 
1029                               #min_target_margin=17, arrowsize=15,
1030                               arrows = True, arrowstyle = '-|>')                           
1031            nx.draw_networkx_edge_labels(self, pos, edge_labels, label_pos=label_pos) #, verticalalignment='center',)        
1032        plt.axis('off');    
1033        
1034        return plt.gcf(), plt.gca(), pos
1035    
1036    def animateSimulation(self, expectedTimePerState=1, inNotebook=False, **kwargs):
1037        """
1038        Animate a simulation trajectory through the state transition graph.
1039        The animation of the CTMC simulation run either in a Jupyter notebook or as a regular script.
1040    
1041        This method visualizes the path of a simulated or predefined sequence of states by
1042        temporarily highlighting each visited state over time. The time spent in each state
1043        is scaled relative to the average sojourn time to produce a visually smooth animation.
1044    
1045        Parameters
1046        ----------
1047        - **expectedTimePerState** : float, optional
1048            Approximate duration (in seconds) for an average state visit in the animation.
1049            The actual pause duration for each state is scaled proportionally to its sojourn time.
1050            Default is 1 second.
1051            
1052        - **inNotebook** : bool, optional
1053            If True, uses Jupyter-compatible animation (with display + clear_output).
1054            If False, uses standard matplotlib interactive animation. Default is True.
1055    
1056        - `**kwargs` : dict
1057            Additional keyword arguments including:
1058                - states : list
1059                    A list of visited states (nodes) to animate.
1060                - times : list or np.ndarray
1061                    Sojourn times in each state (same length as `states`).
1062                - startNode : hashable
1063                    Optional start node for automatic simulation if `states` and `times` are not provided.
1064                - numEvents : int
1065                    Number of events (state transitions) to simulate.
1066                - Additional drawing-related kwargs are passed to `drawTransitionGraph()`.
1067    
1068        Raises
1069        ------
1070        ValueError
1071            If neither a `(states, times)` pair nor `(startNode, numEvents)` are provided.
1072    
1073        Returns
1074        -------
1075        None
1076            The animation is shown interactively using matplotlib but no data is returned.
1077    
1078        Notes
1079        -----
1080        - If `states` and `times` are not supplied, a simulation is run via `simulateMarkovModel()`.
1081        - The function uses matplotlib's interactive mode (`plt.ion()`) to animate transitions.
1082        - The average sojourn time across all states is used to normalize animation speed.
1083        - Each visited state is highlighted in red for its corresponding (scaled) dwell time.
1084        - Drawing arguments like layout, font size, or color can be customized via kwargs.
1085        """        
1086        # expTimePerState: in seconds
1087        if "states" in kwargs and "times" in kwargs:
1088            states = kwargs["states"]
1089            times = kwargs["times"]
1090        else:
1091            # Run default simulation if data not provided
1092            if "startNode" in kwargs and "numEvents" in kwargs:
1093                states, times = self.simulateMarkovModel(startNode=None, numEvents=10)
1094            else:
1095                raise ValueError("Must provide either ('states' and 'times') or ('startNode' and 'numEvents').")
1096        
1097        #avg_OutGoingRate = np.mean([self.nodes[n].get("sumOutgoingRates", 0) for n in self.nodes])
1098        avg_Time = np.mean([1.0/self.nodes[n].get("sumOutgoingRates", 0) for n in self.nodes])
1099        
1100            
1101        # Remove simulation-related keys before drawing
1102        draw_kwargs = {k: v for k, v in kwargs.items() if k not in {"states", "times", "startNode", "numEvents"}}
1103                
1104        fig, ax, pos = self.drawTransitionGraph(**draw_kwargs)
1105        plt.tight_layout()
1106        
1107        if not inNotebook:
1108            plt.ion()
1109                                
1110        for node, time in zip(states,times):
1111            if inNotebook:
1112                clear_output(wait=True)
1113                artist = nx.draw_networkx_nodes(self, pos, nodelist=[node], ax=ax,
1114                                                node_color=self._const.COLOR_NODE_ANIMATION_HIGHLIGHT, node_size=1000)
1115                display(fig)
1116            else:
1117                artist = nx.draw_networkx_nodes(self, pos, nodelist=[node],ax=ax,
1118                                       node_color=self._const.COLOR_NODE_ANIMATION_HIGHLIGHT, node_size=1000)
1119                plt.draw()                        
1120            plt.pause(time/avg_Time*expectedTimePerState)
1121            
1122            # Remove highlighted node    
1123            artist.remove()
1124        
1125        pass
1126    
1127        
1128    
1129    def states(self, data=False):
1130        """
1131        Return a view of the graph's states (nodes), optionally with attributes.
1132    
1133        This method is functionally identical to `self.nodes()` in NetworkX, but is
1134        renamed to `states()` to reflect the semantics of a state transition graph or
1135        Markov model, where nodes represent system states.
1136    
1137        Parameters
1138        ----------
1139        **data** : bool, optional
1140            If True, returns a view of `(node, attribute_dict)` pairs.
1141            If False (default), returns a view of node identifiers only.
1142    
1143        Returns
1144        -------
1145        *networkx.classes.reportviews.NodeView* or NodeDataView
1146            A view over the graph’s nodes or `(node, data)` pairs, depending on the `data` flag.
1147    
1148        Examples
1149        --------
1150        >>> G.states()
1151        ['A', 'B', 'C']
1152    
1153        >>> list(G.states(data=True))
1154        [('A', {'color': 'red', 'label': 'Start'}), ('B', {...}), ...]
1155    
1156        Notes
1157        -----
1158        - This is a convenience wrapper for semantic clarity in models where nodes represent states.
1159        - The view is dynamic: changes to the graph are reflected in the returned view.
1160        """
1161        return self.nodes(data=data)
1162
1163    
1164    def setStateColor(self, state, color):
1165        """
1166        Set the color attribute of a specific state (node).
1167    
1168        Parameters
1169        ----------
1170        **state** : hashable
1171            The identifier of the state whose color is to be updated.
1172    
1173        **color** : str
1174            The new color to assign to the state (used for visualization).
1175    
1176        Returns
1177        -------
1178        None
1179        """
1180        self.nodes[state]["color"] = color
1181        
1182        
1183    def setStateLabel(self, state, label):
1184        """
1185        Set the label attribute of a specific state (node).
1186    
1187        Parameters
1188        ----------
1189        **state** : hashable
1190            The identifier of the state whose label is to be updated.
1191    
1192        **label** : str
1193            The new label to assign to the state (used for visualization or annotation).
1194    
1195        Returns
1196        -------
1197        None
1198        """
1199        self.nodes[state]["label"] = label
1200
1201    def prob(self, state):
1202        """
1203        Return the steady-state probability of a specific state.
1204    
1205        Parameters
1206        ----------
1207        **state** : hashable
1208            The identifier of the state whose probability is to be retrieved.
1209    
1210        Returns
1211        -------
1212        **float**
1213            The steady-state probability of the specified state.
1214    
1215        Raises
1216        ------
1217        KeyError
1218            If the state is not found in the steady-state probability dictionary.
1219        """
1220        return self.state_probabilities[state]
1221        
1222    
1223    def probs(self):
1224        """
1225        Return the steady-state probabilities for all states.
1226    
1227        Returns
1228        -------
1229        **dict**
1230            A dictionary mapping each state to its steady-state probability.
1231    
1232        Notes
1233        -----
1234        - The probabilities are computed via `calcSteadyStateProbs()` and stored internally.
1235        - Use this method to access the full steady-state distribution.
1236        """
1237        return self.state_probabilities
1238
1239    def getEmptySystemProbs(self, state=None):
1240        """
1241        Create an initial state distribution where all probability mass is in a single state.
1242    
1243        By default, this method returns a distribution where the entire probability mass
1244        is placed on the first node in the graph. Alternatively, a specific starting state
1245        can be provided.
1246    
1247        Parameters
1248        ----------
1249        **state** : hashable, optional
1250            The state to initialize with probability 1. If None, the first node
1251            in the graph (based on insertion order) is used.
1252    
1253        Returns
1254        -------
1255        **dict**
1256            A dictionary representing the initial distribution over states.
1257            The selected state has probability 1, and all others have 0.
1258    
1259        Examples
1260        --------
1261        >>> X0 = G.getEmptySystem()
1262        >>> sum(X0.values())  # 1.0
1263    
1264        Notes
1265        -----
1266        - This method is useful for setting up deterministic initial conditions
1267          in transient simulations of a CTMC.
1268        - The return format is compatible with methods like `transientProbs()`.
1269        """
1270        if state is None:
1271            first_node = next(iter(self.nodes))
1272        X0 = {s: 1 if s == first_node else 0 for s in self}
1273        return X0
1274
1275        
1276
1277    def transientProbs(self, initialDistribution, t):
1278        """
1279        Compute the transient state probabilities at time `t` for the CTMC.
1280    
1281        This method solves the forward equation of a continuous-time Markov chain (CTMC):
1282        X(t) = X0 · expm(Q·t), where `Q` is the generator (rate) matrix and `X0` is the
1283        initial state distribution. It returns a dictionary mapping each state to its
1284        probability at time `t`.
1285    
1286        Parameters
1287        ----------
1288        - **initialDistribution** : dict 
1289            A dictionary mapping states (nodes) to their initial probabilities at time `t=0`.
1290            The keys must match the nodes in the graph, and the values should sum to 1.    
1291            
1292        - **t** : float
1293            The time at which the transient distribution is to be evaluated.
1294    
1295        Returns
1296        -------
1297        **dict**
1298            A dictionary mapping each state (node) to its probability at time `t`.
1299    
1300        Notes
1301        -----
1302        - The rate matrix `Q` is generated automatically if not yet computed.
1303        - The state order in the rate matrix corresponds to the internal `_n2i` mapping.
1304        - Internally, the dictionary `initialDistribution` is converted into a vector
1305          aligned with the state indexing.
1306        - The transient solution is computed using `scipy.linalg.expm` for matrix exponentiation.
1307    
1308        Raises
1309        ------
1310        ValueError
1311            If the input dictionary has missing states or probabilities that do not sum to 1.
1312            (This is not enforced but recommended for correctness.)
1313        """        
1314        if self._rate_matrix is None:  self.createRateMatrix()
1315        Q = self._rate_matrix
1316        X0 = np.zeros(len(initialDistribution))
1317        
1318        for i,n in enumerate(self.nodes):
1319            X0[ self._n2i[n] ] = initialDistribution[n]
1320                
1321        
1322        X = X0 @ expm(Q*t)    
1323        
1324        matrix = {}
1325        for n in self.nodes:
1326            matrix[n] = X[ self._n2i[n] ]
1327                
1328        return matrix
1329        
1330    def symSolveMarkovModel(self):
1331        """
1332        Symbolically solve the global balance equations of the CTMC.
1333    
1334        This method constructs and solves the system of symbolic equations for the
1335        steady-state probabilities of each state in the CTMC. It uses symbolic transition
1336        rates stored on the edges to form balance equations for each node, along with
1337        the normalization constraint that all probabilities sum to 1.
1338    
1339        Returns
1340        -------
1341        - **solution** : dict
1342            A symbolic solution dictionary mapping SymPy symbols (e.g., x_{A}) to
1343            expressions in terms of model parameters.
1344    
1345        - **num_solution** : dict
1346            A numerical dictionary mapping each state (node) to its steady-state probability,
1347            computed by substituting the numeric values from `self._subs` into the symbolic solution.
1348    
1349        Notes
1350        -----
1351        - For each node `s`, a symbolic variable `x_{label}` is created using the node's label attribute.
1352        - One balance equation is created per state: total inflow = total outflow.
1353        - An additional constraint `sum(x_i) = 1` ensures proper normalization.
1354        - The symbolic system is solved with `sympy.solve`, and the results are simplified.
1355        - Numeric values are computed by substituting known parameter values (`self._subs`) into the symbolic solution.
1356    
1357        Examples
1358        --------
1359        >>> sym_sol, num_sol = G.symSolveMarkovModel()
1360        >>> sym_sol  # symbolic expressions for each state
1361        >>> num_sol  # evaluated numerical values for each state
1362        """
1363        X={}
1364        to_solve = []
1365        for s in self.nodes():
1366            X[s] = sp.Symbol(f'x_{{{self.nodes[s]["label"]}}}')
1367            to_solve.append(X[s])
1368        #%
1369        eqs = []
1370        for s in self.nodes():
1371             succ = self.successors(s) # raus
1372             out_rate = 0
1373             for z in list(succ):
1374                 #out_rate += X[s]*self.edges[s, z]["sym_rate"]                 
1375                 out_rate += X[s]*self[s][z]["sym_rate"]
1376             
1377             in_rate = 0
1378             for r in list(self.predecessors(s)):
1379                 #in_rate += X[r]*self.edges[r,s]["sym_rate"] 
1380                 in_rate += X[r]*self[r][s]["sym_rate"] 
1381             eqs.append( sp.Eq(out_rate, in_rate) )
1382        #% 
1383        Xsum = 0    
1384        for s in X:
1385            Xsum += X[s]
1386        
1387        eqs.append( sp.Eq(Xsum, 1) )    
1388        #%
1389        solution = sp.solve(eqs, to_solve)
1390        simplified_solution = sp.simplify(solution)
1391        #print(simplified_solution)
1392        num_solution = simplified_solution.subs(self._subs)
1393        #num_dict =  {str(var): str(sol) for var, sol in simplified_solution.items()}
1394        num_dict = {s: num_solution[X[s]] for s in X}                
1395        
1396        return simplified_solution, num_dict 
1397    
1398    def exportTransitionsToExcel(self, excel_path = "output.xlsx", num_colors=9, open_excel=False, verbose=True):
1399        """
1400        Export transition data from the state transition graph to an Excel file.
1401        
1402        This method collects all transitions (edges) from the graph and writes them
1403        to an Excel spreadsheet, optionally applying background colors to distinguish
1404        outgoing transitions from different source states. This export is helpful for
1405        manually verifying the transition structure of the model.
1406
1407        The transition type is defined by the `tid` field of each edge, and mapped to
1408        a human-readable string via the `EXPORT_NAME_TRANSITIONS` dictionary in the
1409        constants module (e.g., `Constants.EXPORT_NAME_TRANSITIONS[tid]`).
1410
1411        See Also
1412        --------
1413        StateTransitionGraph.addTransition : Method to add edges with symbolic rates and `tid`.
1414        StateTransitionGraph.__init__ : Initializes the graph and loads the constants module.
1415
1416        Parameters
1417        ----------
1418        - **excel_path** : str, optional
1419            Path to the output Excel file. Default is "output.xlsx".
1420        
1421        - num_colors : int, optional
1422            Number of unique background colors to cycle through for different states.
1423            Default is 9.
1424        
1425        - open_excel : bool, optional
1426            If True, automatically opens the generated Excel file after export.
1427            Default is False.
1428        
1429        - verbose : bool, optional
1430            If True, prints summary information about the export process.
1431            Default is True.
1432        
1433        Returns
1434        -------
1435        None
1436            The function writes data to disk and optionally opens Excel,
1437            but does not return a value.
1438
1439        Notes
1440        -----
1441        - Each edge must contain a `tid` attribute for transition type and a `label` for the rate.
1442        - Background coloring groups all outgoing transitions from the same source state.
1443        - Requires the constants module to define `EXPORT_NAME_TRANSITIONS`.
1444        """
1445
1446        rows = []
1447        for node in self.nodes():
1448            edges = list(self.edges(node))
1449            if verbose: print(f"Edges connected to node '{node}':")
1450            
1451            for i, (u, v) in enumerate(edges):
1452                latex = self[u][v]["label"]
1453                rate = sp.pretty(self[u][v]["sym_rate"])
1454                tid = self[u][v]["tid"]
1455                
1456                if tid: # 
1457                    if hasattr(self._const, 'EXPORT_NAME_TRANSITIONS'):
1458                        tid_name = self._const.EXPORT_NAME_TRANSITIONS[tid]
1459                    else:
1460                        tid_name  = ''                               
1461                
1462                rows.append({"from": u, "type": tid_name, "to": v, "rate": rate, "latex": latex})
1463            
1464        df = pd.DataFrame(rows, columns=["from", "type", "to", "rate", "latex"])
1465        df["checked"] = ""
1466        df.to_excel(excel_path, index=False)
1467        
1468        # Load workbook with openpyxl
1469        wb = load_workbook(excel_path)
1470        ws = wb.active
1471        
1472        # Define colors for different "from" values
1473        cmap = map(plt.cm.Pastel1, range(num_colors))
1474        # Convert RGBA to hex
1475        color_list = [mcolors.to_hex(color, keep_alpha=False).upper().replace("#", "") for color in cmap]
1476        
1477        # Assign color to each unique "from" value using modulo
1478        unique_from = (df["from"].unique())  # Sorted for consistent order
1479        from_color_map = {
1480            value: color_list[i % len(color_list)] for i, value in enumerate(unique_from)
1481        }
1482        
1483        # Apply coloring row by row
1484        for i,row in enumerate(ws.iter_rows(min_row=2, max_row=ws.max_row)):  # Skip header
1485            from_value = row[0].value  # "from" is first column
1486            #fill_color = color_list[i % len(color_list)]
1487            if isinstance(from_value, int):
1488                fill_color = from_color_map[from_value]
1489            else:
1490                fill_color = from_color_map.get(eval(from_value))
1491            
1492            if fill_color:
1493                fill = PatternFill(start_color=fill_color, end_color=fill_color, fill_type="solid")
1494                for cell in row:
1495                    cell.fill = fill
1496        
1497        # Save styled Excel file
1498        wb.save(excel_path)
1499        if verbose: print(f"File {excel_path} created.")
1500        if open_excel: 
1501            print(f"Opening the file {excel_path}.")
1502            os.startfile(excel_path)
1503    
class StateTransitionGraph(networkx.classes.digraph.DiGraph):
  95class StateTransitionGraph(nx.DiGraph):
  96    """
  97    Initialize the StateTransitionGraph for continuous time Markov chains (CTMC). 
  98    Extends the networkx.Graph class to include additional methods and properties 
  99    for the analysis, visualization, and simulation of CTMCs. 
 100    The analysis allows deriving the state probabilities in the steady-state 
 101    and in the transient phase.
 102
 103    This constructor sets up symbolic and numeric parameters for the state transition graph,
 104    imports a constants module dynamically, and initializes internal data structures for
 105    rate matrix computation and steady-state probability analysis.
 106
 107    Parameters
 108    ----------
 109    - **parameters** : dict  
 110      A dictionary mapping parameter names to their numerical values.     
 111        Example: {'lambda': 1.0, 'mu': 1.0}. 
 112        For each parameter key:
 113        - A symbolic variable is created using sympy.symbols. Can be accessed as `G.sym("key")` or `G.sym_key`
 114        - A variable is created which can be accessed via `G.key` (and the numerical value is assigned)
 115        - For the key="lambda", `G.lam` and `G.sym_lambda` are created
 116    
 117        The parameters can be accessed via property `parameters`
 118    - **constants** : str, optional
 119        The name of a Python module (as string) to import dynamically and assign to `self.const`.
 120        This can be used to store global or user-defined constants.
 121        Default is "Constants".
 122    - `*args` : tuple
 123        Additional positional arguments passed to the base `networkx.DiGraph` constructor.
 124    - `**kwargs` : dict
 125        Additional keyword arguments passed to the base `networkx.DiGraph` constructor.
 126
 127    Attributes
 128    ----------
 129    - `_rate_matrix` : ndarray or None  
 130      The continuous-time rate matrix `Q` of the CTMC. Computed when calling `createRateMatrix()`.
 131    - `_state_probabilities` : dict or None  
 132      Dictionary of steady-state probabilities, computed via `calcSteadyStateProbs()`.
 133    - `_n2i` : dict or None  
 134      Mapping from node identifiers to matrix row/column indices used for rate matrix construction.
 135    - `_sym` : dict  
 136      Mapping of parameter names (strings) to their symbolic `sympy` representations.
 137    - `_subs` : dict  
 138      Mapping from symbolic parameters to their numerical values (used for substitutions).
 139    - `_parameters` : dict  
 140      Copy of the original parameter dictionary passed at initialization.
 141    - `_transitionProbabilitesComputed` : bool  
 142      Internal flag indicating whether transition probabilities for the embedded Markov chain
 143      have been computed.
 144    - `const` : module  
 145      Dynamically imported module (default `"Constants"`), used for storing user-defined constants 
 146      such as colors or transition IDs.
 147        
 148    Inherits
 149    --------
 150        networkx.Graph: The base graph class from NetworkX.
 151    """
 152
 153    def __init__(self, parameters, constants="Constants", *args, **kwargs):
 154        """
 155        Initialize the state transition graph with symbolic parameters.
 156
 157        This constructor sets up the symbolic and numerical environment for modeling
 158        a CTMC. Parameters are stored as both symbolic expressions and numerical values,
 159        and convenient attribute access is provided for both forms.
 160
 161        Parameters
 162        ----------
 163        - **parameters** : dict  
 164          A dictionary mapping parameter names to their numerical values.     
 165            Example: {'lambda': 1.0, 'mu': 1.0}. 
 166            For each parameter key:
 167            - A symbolic variable is created using sympy.symbols. Can be accessed as `G.sym("key")` or `G.sym_key`
 168            - A variable is created which can be accessed via `G.key` (and the numerical value is assigned)
 169            - For the key="lambda", `G.lam` and `G.sym_lambda` are created
 170        
 171            The parameters can be accessed via property `parameters`
 172        - **constants** : str, optional
 173            The name of a Python module (as string) to import dynamically and assign to `self.const`.
 174            This can be used to store global or user-defined constants.
 175            Default is "Constants".
 176        - `*args` : tuple
 177            Additional positional arguments passed to the base `networkx.DiGraph` constructor.
 178        - `**kwargs` : dict
 179            Additional keyword arguments passed to the base `networkx.DiGraph` constructor.
 180        """        
 181        super().__init__(*args, **kwargs)
 182        self._rate_matrix = None  # Rate Matrix Q of this transition graph
 183        self._state_probabilities = None  # steady state proabilities of this CTMC
 184        self._n2i = None # internal index for generating the rate matrix Q: dictionary (node to index)
 185        
 186        #parameters = {'lambda': 1.0, 'mu':1.0, 'n':5}
 187        self._parameters = parameters
 188        self._sym = {key: sp.symbols(key) for key in parameters}
 189        self._subs = {} # contains symbolic parameter as key and provides numerical value for computation.
 190        for k in parameters:
 191            self._subs[ self._sym[k] ] = parameters[k] 
 192
 193        # access to parameters and symbols using dot-notation, e.g. G.lam, G.sym_lambda, or G.sym("lambda")
 194        for key, value in parameters.items():
 195            if key=='lambda':
 196                setattr(self, 'lam', value)
 197            else:
 198                setattr(self, key, value)
 199            
 200            setattr(self, "sym_"+key, self._sym[key])
 201            
 202        self._transitionProbabilitesComputed = False
 203        
 204        self._const = importlib.import_module(constants)
 205
 206        
 207        
 208    def addState(self, state, color=None, label=None):
 209        """
 210        Add a state (node) to the graph with optional color and label attributes.
 211    
 212        If the node already exists, its attributes are updated with the provided values.
 213    
 214        Parameters
 215        ----------
 216        - **state** : hashable
 217            The identifier for the state. Can be a string, int, or a tuple. If no label
 218            is provided, the label will be auto-generated as:
 219            - comma-separated string for tuples (e.g., `(1, 2)` -> `'1,2'`)
 220            - stringified value for other types
 221    
 222        - **color** : str, optional
 223            The color assigned to the state (used for visualization).
 224            If not specified, defaults to `self.const.COLOR_NODE_DEFAULT`.
 225    
 226        - **label** : str, optional
 227            A label for the node (used for display or annotation). If not provided,
 228            it will be auto-generated from the state identifier.
 229    
 230        Notes
 231        -----
 232        - This method wraps `self.add_node()` from NetworkX.
 233        - If the node already exists in the graph, this method does not raise an error;
 234          it will update the node’s existing attributes.
 235        """
 236        if label is None:
 237            if isinstance(state, tuple):
 238                label = ','.join(map(str, state)) # if state is a tuple
 239            else: 
 240                label = str(state)
 241            #self.nodes[state]["label"] = label
 242        
 243        if color is None:
 244            color = self._const.COLOR_NODE_DEFAULT
 245        
 246        # no error if the node already exists.Attributes get updated if you pass them again.
 247        self.add_node(state, color=color, label=label)
 248        
 249    def setStateProperties(self, state, color=None, label=None, **kwargs):
 250        """
 251        Update visual or metadata properties for an existing state (node) in the graph.
 252    
 253        This method updates the node's attributes such as color, label, and any additional
 254        custom key-value pairs. The node must already exist in the graph.
 255    
 256        Parameters
 257        ----------
 258        - **state** : hashable
 259            Identifier of the state to update. Must be a node already present in the graph.
 260    
 261        - **color** : str, optional
 262            Color assigned to the node, typically used for visualization. If not provided,
 263            the default node color (`self._const.COLOR_NODE_DEFAULT`) is used.
 264    
 265        - **label** : str, optional
 266            Label for the node. If not provided, it is auto-generated from the state ID
 267            (especially for tuples).
 268    
 269        - `**kwargs` : dict
 270            Arbitrary additional key-value pairs to set as node attributes. These can include
 271            any custom metadata.
 272    
 273        Raises
 274        ------
 275        KeyError
 276            If the specified state is not present in the graph.
 277    
 278        Notes
 279        -----
 280        - This method ensures that required attributes (like `color` and `label`) are set for visualization. 
 281        - Existing node attributes are updated or extended with `kwargs`.
 282        """        
 283        if state not in self:
 284            raise KeyError(f"State '{state}' does not exist in the graph.")
 285        
 286        # mandatory for plotting
 287        self.addState(state, color, label)
 288        
 289        # # Add any additional attributes
 290        self.nodes[state].update(kwargs)
 291        
 292    def setAllStateDefaultProperties(self, **kwargs):
 293        """
 294        Set default visual and metadata properties for all states (nodes) in the graph.
 295    
 296        This method applies `setStateProperties()` to every node, ensuring each state
 297        has at least the required default attributes (e.g., color and label).
 298        Any additional keyword arguments are passed to `setStateProperties()`.
 299    
 300        Parameters
 301        ----------
 302        `**kwargs` : dict, optional
 303            Optional keyword arguments to apply uniformly to all nodes,
 304            such as color='gray', group='core', etc.
 305    
 306        Notes
 307        -----
 308        - This method is useful for initializing node attributes before plotting
 309          or performing other graph-based computations.
 310        - Nodes without existing attributes will receive defaults via `addState()`.
 311        """
 312        for state in self.nodes():
 313            self.setStateProperties(state, **kwargs)
 314
 315        
 316    #def addTransition(self, origin_state, list_destination_state, sym_rate, color=plt.cm.tab10(5), tid = None):        
 317    def addTransition(self, origin_state, destination_state, sym_rate, tid = None, color=None, label=None):     
 318        """
 319        Add a directed transition (edge) between two states with a symbolic transition rate.
 320    
 321        This method adds a directed edge from `origin_state` to `destination_state`,
 322        stores both the symbolic and numerical rate, and sets optional display properties
 323        such as color, label, and transition ID. If the edge already exists, a warning is issued.
 324    
 325        Parameters
 326        ----------
 327        - **origin_state** : hashable
 328            The source node (state) of the transition.
 329    
 330        - **destination_state** : hashable
 331            The target node (state) of the transition.
 332    
 333        - **sym_rate** : sympy.Expr
 334            A symbolic expression representing the transition rate.
 335            This will be numerically evaluated using internal substitutions (`self._subs`).
 336    
 337        - **tid** : str or int, optional
 338            A transition identifier (e.g. for grouping or styling).
 339            If provided, the color will default to `self._const.COLOR_TRANSITIONS[tid]` if not set explicitly.
 340    
 341        - **color** : str, optional
 342            Color for the edge (used in visualization). Defaults to the transition-specific color
 343            if `tid` is provided.
 344    
 345        - **label** : str, optional
 346            LaTeX-formatted string for labeling the edge. If not provided, a default label is
 347            generated using `sympy.latex(sym_rate)`.
 348    
 349        Notes
 350        -----
 351        - The numeric rate is computed by substituting known parameter values into `sym_rate`.
 352        - The following attributes are stored on the edge:
 353            - `rate` (float): numeric value of the transition rate
 354            - `sym_rate` (sympy.Expr): original symbolic expression
 355            - `label` (str): display label
 356            - `color` (str): edge color
 357            - `tid` (str or int): transition ID
 358        - If the edge already exists, a `UserWarning` is issued but the edge is still overwritten.
 359        """        
 360        if tid is not None:
 361            if color is None:
 362                color = self._const.COLOR_TRANSITIONS[ tid ]
 363                
 364        rate = sym_rate.subs( self._subs )
 365        
 366        if label is None: string = "$"+sp.latex(sym_rate)+"$"
 367                
 368        if self.has_edge(origin_state, destination_state):
 369            theedge = self[origin_state][destination_state]
 370            warningstr = f"Edge already exists: {origin_state} -> {destination_state} with rate {theedge['label']}"
 371            warnings.warn(warningstr, category=UserWarning)
 372        self.add_edge(origin_state, destination_state, rate=float(rate), label=string, 
 373                   color=color, sym_rate=sym_rate, tid=tid)
 374    
 375    def sym(self, key):
 376        """
 377        Retrieve the symbolic representation of a model parameter.
 378    
 379        Parameters
 380        ----------
 381        **key** : str
 382            The name of the parameter (e.g., 'lambda', 'mu', 'n').
 383    
 384        Returns
 385        -------
 386        **sympy.Symbol**
 387            The symbolic variable corresponding to the given key.
 388    
 389        Raises
 390        ------
 391        KeyError
 392            If the key does not exist in the internal symbolic mapping (`_sym`).
 393    
 394        Examples
 395        --------
 396        >>> G.sym('lambda')
 397        lambda
 398    
 399        Notes
 400        -----
 401        - The symbolic mapping is initialized during object construction based on the `parameters` dictionary.
 402        - To access the symbol using dot-notation, use `.sym_key` attributes created during initialization.
 403        """
 404        return self._sym[key]
 405
 406        
 407    @property
 408    def parameters(self):
 409        """
 410        Return the model parameters as a dictionary.
 411    
 412        This includes the original numerical values provided during initialization.
 413    
 414        Returns
 415        -------
 416        **dict**
 417            Dictionary of parameter names and their corresponding numeric values.
 418    
 419        Examples
 420        --------
 421        >>> G.parameters
 422        {'lambda': 1.0, 'mu': 0.5, 'n': 3}
 423    
 424        Notes
 425        -----
 426        - This property does not return symbolic values. Use `self.sym(key)` for symbolic access.
 427        - Parameter symbols can also be accessed using attributes like `sym_lambda` or the `sym()` method.
 428        """
 429        return self._parameters
 430
 431    
 432
 433    @property
 434    def rate_matrix(self):
 435        """
 436        Return the continuous-time rate matrix (Q matrix) of the state transition graph.
 437    
 438        This matrix defines the transition rates between all pairs of states in the
 439        continuous-time Markov chain (CTMC) represented by this graph.
 440    
 441        Returns
 442        -------
 443        **numpy.ndarray** or None
 444            The rate matrix `Q`, where `Q[i, j]` is the transition rate from state `i` to `j`.
 445            Returns `None` if the rate matrix has not been computed yet.
 446    
 447        Notes
 448        -----
 449        - The matrix is typically generated via a method like `computeRateMatrix()` or similar.
 450        - Internal state indexing is handled by `_n2i` (node-to-index mapping).
 451        - The rate matrix is stored internally as `_rate_matrix`.
 452        """
 453        return self._rate_matrix
 454
 455
 456
 457    @property
 458    def state_probabilities(self):
 459        """
 460        Return the steady-state probabilities of the continuous-time Markov chain (CTMC).
 461    
 462        These probabilities represent the long-run proportion of time the system spends
 463        in each state, assuming the CTMC is irreducible and positive recurrent.
 464    
 465        Returns
 466        -------
 467        **dict** or None
 468            A dictionary mapping each state to its steady-state probability.
 469            Returns `None` if the probabilities have not yet been computed.
 470    
 471        Notes
 472        -----
 473        - The probabilities are stored internally after being computed via a dedicated method
 474          `calcSteadyStateProbs()`.
 475        - The result is stored in the `_state_probabilities` attribute.
 476        """
 477        return self._state_probabilities
 478
 479     
 480    
 481    def printEdges(self, state):
 482        """
 483        Print all outgoing edges from a given state in the graph.
 484    
 485        For each edge from the specified `state`, this method prints the destination
 486        state and the corresponding transition label (typically the symbolic rate).
 487    
 488        Parameters
 489        ----------
 490        **state** : hashable
 491            The source node (state) from which outgoing edges will be printed.            
 492    
 493        Returns
 494        -------
 495        None
 496            This method prints output directly to the console and does not return a value.
 497    
 498        Notes
 499        -----
 500        - The transition label is retrieved from the edge attribute `'label'`.
 501        - Assumes `state` exists in the graph; otherwise, a KeyError will be raised.
 502        """
 503        print(f"from {state}")
 504        node = self[state]
 505        for e in node:
 506            rate = node[e]['label']
 507            print(f" -> {e} with {rate}")
 508
 509            
 510    def createRateMatrix(self):
 511        """
 512        Construct the continuous-time Markov chain (CTMC) rate matrix Q from the graph.
 513    
 514        This method generates the generator matrix `Q` corresponding to the structure and
 515        rates defined in the directed graph. Each node represents a system state, and each
 516        directed edge must have a `'rate'` attribute representing the transition rate from
 517        one state to another.
 518    
 519        The matrix `Q` is constructed such that:
 520        - `Q[i, j]` is the transition rate from state `i` to state `j`
 521        - Diagonal entries satisfy `Q[i, i] = -sum(Q[i, j]) for j ≠ i`, so that each row sums to 0
 522    
 523        The resulting matrix and the internal node-to-index mapping (`_n2i`) are stored
 524        as attributes of the graph.
 525    
 526        Returns
 527        -------
 528        None
 529            The method updates internal attributes:
 530            - `_rate_matrix` : numpy.ndarray
 531                The CTMC rate matrix Q
 532            - `_n2i` : dict
 533                Mapping from node identifiers to matrix row/column indices
 534    
 535        Notes
 536        -----
 537        - Assumes all edges in the graph have a `'rate'` attribute.
 538        - The node order used in the matrix is consistent with the order from `list(self.nodes)`.
 539        - Use `self.rate_matrix` to access the generated matrix after this method is called.
 540        """
 541        n2i = {}
 542        nodes = list(self.nodes)
 543        for i,node in enumerate(nodes):
 544            n2i[node] = i         
 545        Q = np.zeros((len(nodes),len(nodes))) # rate matrix
 546        
 547        for edge in self.edges:
 548            i0 = n2i[edge[0]]
 549            i1 = n2i[edge[1]]
 550            Q[i0,i1] = self[edge[0]][edge[1]]["rate"] 
 551            
 552        np.fill_diagonal(Q, -Q.sum(axis=1))
 553        self._rate_matrix = Q
 554        self._n2i = n2i
 555        #return Q, n2i
 556    
 557    def calcSteadyStateProbs(self, verbose=False):        
 558        """
 559        Compute the steady-state probabilities of the continuous-time Markov chain (CTMC).
 560    
 561        This method calculates the stationary distribution of the CTMC defined by the rate
 562        matrix of the state transition graph. The computation solves the linear system
 563        `XQ = 0` with the constraint that the probabilities X sum to 1, using a modified
 564        version of the rate matrix `Q`.
 565    
 566        Parameters
 567        ----------
 568        `verbose` : bool, optional
 569            If True, prints the rate matrix, the modified matrix, the right-hand side vector,
 570            and the resulting steady-state probabilities.
 571    
 572        Returns
 573        -------
 574        **probs** : dict
 575            A dictionary mapping each state (node in the graph) to its steady-state probability.
 576    
 577        Notes
 578        -----
 579        - The graph must already contain valid transition rates on its edges.
 580        - If the rate matrix has not yet been created, `createRateMatrix()` is called automatically.
 581        - The method modifies the last column of the rate matrix to impose the normalization constraint
 582          `sum(X) = 1`.
 583        - The inverse of the modified matrix is computed directly, so the graph should be small enough
 584          to avoid numerical instability or performance issues.
 585        - The resulting probabilities are stored internally in `_state_probabilities` and can be accessed
 586          via the `state_probabilities` property.
 587        """
 588        
 589        # compute transition matrix if not already done
 590        if self._rate_matrix is None:  self.createRateMatrix()
 591        Q2 =  self._rate_matrix.copy()
 592        if verbose:
 593            print("Rate matrix Q")
 594            print(f'Q=\n{Q2}\n')
 595        
 596        Q2[:, -1] = 1
 597        if verbose:        
 598            print(f'Matrix is changed to\nQ2=\n{Q2}\n')
 599
 600        b = np.zeros(len(Q2))
 601        b[-1] = 1
 602        if verbose:
 603            print("Solve X = b @ Q2")
 604            print(f'b=\n{b}\n')
 605        
 606        # state probabilities
 607        X = b @ linalg.inv(Q2) # compute the matrix inverse
 608            
 609        # Generate a matrix with P(X,Y)
 610        matrix = {}
 611        for n in self.nodes:
 612            matrix[n] = X[ self._n2i[n] ]
 613        
 614        if verbose:
 615            print("Steady-state probabilities X")
 616            print(f'X=\n{matrix}\n')
 617        self._state_probabilities = matrix
 618        return matrix
 619        
 620    def calcTransitionProbabilities(self):
 621        """
 622        Compute and store transition probabilities for each node in the graph of the embedded discrete Markov chain.
 623    
 624        This method processes each node in the graph, treating outgoing edges as
 625        transitions in an embedded discrete-time Markov chain derived from a continuous-time
 626        Markov chain (CTMC). Transition *rates* on outgoing edges are normalized into
 627        *probabilities*, and relevant data is stored in the node attributes.
 628    
 629        For each node, the following attributes are set:
 630        - `"transitionProbs"` : np.ndarray
 631            Array of transition probabilities to successor nodes.
 632        - `"transitionNodes"` : list
 633            List of successor node identifiers (ordered as in the probability array).
 634        - `"sumOutgoingRates"` : float
 635            Sum of all outgoing transition rates. Used to model the exponential waiting time
 636            in continuous-time simulation (`T ~ Exp(sumOutgoingRates)`).
 637    
 638        Notes
 639        -----
 640        - The graph must be a directed graph where each edge has a `'rate'` attribute.
 641        - Nodes with no outgoing edges are considered absorbing states.
 642        - Useful for discrete-event simulation and Markov chain Monte Carlo (MCMC) modeling.
 643        - Sets an internal flag `_transitionProbabilitesComputed = True` upon completion.
 644    
 645        Raises
 646        ------
 647        KeyError
 648            If any outgoing edge lacks a `'rate'` attribute.
 649        """
 650        for node in self:
 651            successors = list(self.successors(node))
 652            rates = np.array([self[node][nbr]['rate'] for nbr in successors])
 653            probs = rates / rates.sum()  # Normalize to probabilities
 654            self.nodes[node]["transitionProbs"] = probs # probabilities
 655            self.nodes[node]["transitionNodes"] = successors # list of nodes
 656            self.nodes[node]["sumOutgoingRates"] = rates.sum() # time in state ~ EXP(sumOutgoingRates)
 657            numel = len(successors)
 658            if numel==0: 
 659                print(f"Absorbing state: {node}")
 660        self._transitionProbabilitesComputed = True
 661                
 662    def simulateMarkovModel(self, startNode=None, numEvents=10):
 663        """
 664        Simulate a trajectory through the state transition graph as a continuous-time Markov chain (CTMC).
 665        
 666        This method performs a discrete-event simulation of the CTMC by generating a random walk
 667        through the graph, starting at `startNode` (or a default node if not specified). Transitions
 668        are chosen based on the normalized transition probabilities computed from edge rates, and
 669        dwell times are sampled from exponential distributions based on the total outgoing rate
 670        from each state.
 671        
 672        Parameters
 673        ----------
 674        - **startNode** : hashable, optional
 675            The node at which the simulation starts. If not specified, the first node in the graph
 676            (based on insertion order) is used.
 677        
 678        - **numEvents** : int, optional
 679            The number of transition events (steps) to simulate. Default is 10.
 680        
 681        Returns
 682        -------
 683        **states** : list
 684            A list of visited states, representing the sequence of nodes traversed in the simulation.
 685        
 686        **times** : numpy.ndarray
 687            An array of dwell times in each state, sampled from exponential distributions
 688            with rate equal to the sum of outgoing transition rates from each state.
 689        
 690        Notes
 691        -----
 692        - This method automatically calls `calcTransitionProbabilities()` if transition probabilities
 693          have not yet been computed.
 694        - If an absorbing state (i.e., a node with no outgoing edges) is reached, the simulation stops early.
 695        - The cumulative time can be obtained by `np.cumsum(times)`.
 696        - State durations follow `Exp(sumOutgoingRates)` distributions per CTMC theory.
 697        """
 698        if not self._transitionProbabilitesComputed: self.calcTransitionProbabilities()
 699        
 700        # simulate random walk through Graph = Markov Chain
 701        if startNode is None:
 702            startNode = next(iter(self))
 703        
 704
 705        states = [] # list of states
 706        times = np.zeros(numEvents) # time per state
 707
 708        node = startNode
 709        for i in range(numEvents):
 710            outgoingRate = self.nodes[node]["sumOutgoingRates"]  # time in state ~ EXP(sumOutgoingRates)
 711            times[i] = random.expovariate(outgoingRate) # exponentially distributed time in state
 712            states.append(node)
 713            
 714            # get next node    
 715            probs = self.nodes[node]["transitionProbs"] # probabilities
 716            successors = self.nodes[node]["transitionNodes"] # list of nodes        
 717            numel = len(successors)
 718            if numel==0: 
 719                print(f"Absorbing state: {node}")
 720                break
 721            elif numel==1:
 722                nextNode = successors[0]
 723            else:
 724                nextNode = successors[ np.random.choice(numel, p=probs) ]
 725            # edge = G[node][nextNode]
 726            node = nextNode
 727            
 728        return states, times
 729    
 730    def calcProbabilitiesSimulationRun(self, states, times):
 731        """
 732        Estimate empirical state probabilities from a simulation run of the CTMC.
 733    
 734        This method computes the fraction of total time spent in each state during a simulation.
 735        It uses the sequence of visited states and corresponding sojourn times to estimate the
 736        empirical (time-weighted) probability distribution over states.
 737    
 738        Parameters
 739        ----------
 740        - **states** : list of tuple
 741            A list of state identifiers visited during the simulation.    
 742        - **times** : numpy.ndarray
 743            An array of sojourn times corresponding to each visited state, same length as `states`.
 744    
 745        Returns
 746        -------
 747        **simProbs** : dict
 748            A dictionary mapping each unique state to its estimated empirical probability,
 749            based on total time spent in that state relative to the simulation duration.
 750    
 751        Notes
 752        -----
 753        - Internally, states are mapped to flat indices for efficient counting using NumPy.
 754        - This method assumes that `states` contains only valid node identifiers from the graph.
 755        - If `_n2i` (node-to-index mapping) has not been initialized, it will be computed here.
 756        - This is especially useful for validating analytical steady-state probabilities
 757          against sampled simulation results. Or in case of state-space explosions.
 758        """        
 759        if self._n2i is None:
 760            n2i = {}
 761            nodes = list(self.nodes)
 762            for i,node in enumerate(nodes):
 763                n2i[node] = i     
 764            self._n2i = n2i
 765                
 766        
 767        # Map tuples to flat indices
 768        data = np.array(states)
 769        max_vals = data.max(axis=0) + 1
 770        indices = np.ravel_multi_index(data.T, dims=max_vals)
 771
 772        # Weighted bincount
 773        bincount = np.bincount(indices, weights=times)
 774
 775        # Decode back to tuples
 776        tuple_indices = np.array(np.unravel_index(np.nonzero(bincount)[0], shape=max_vals)).T
 777        timePerState = bincount[bincount > 0]
 778
 779        simProbs = {}
 780        totalTime = times.sum()
 781        for t, c in zip(tuple_indices, timePerState):
 782            simProbs[tuple(t)] = c / totalTime
 783            #print(f"Tuple {tuple(t)} → weighted count: {c}")
 784        return simProbs
 785    
 786    
 787    def _draw_networkx_edge_labels(
 788        self,
 789        pos,
 790        edge_labels=None, label_pos=0.75,
 791        font_size=10, font_color="k", font_family="sans-serif", font_weight="normal",
 792        alpha=None, bbox=None,
 793        horizontalalignment="center", verticalalignment="center",
 794        ax=None, rotate=True, clip_on=True, rad=0):
 795        """
 796        Draw edge labels for a directed graph with optional curved edge support.
 797    
 798        This method extends NetworkX's default edge label drawing by supporting
 799        curved (bended) edges through quadratic Bézier control points. It places
 800        labels dynamically at a specified position along the curve and optionally
 801        rotates them to follow the edge orientation.
 802    
 803        Parameters
 804        ----------
 805        pos : dict
 806            Dictionary mapping nodes to positions (2D coordinates). Required.
 807    
 808        edge_labels : dict, optional
 809            Dictionary mapping edge tuples `(u, v)` to label strings. If None,
 810            edge data is used and converted to strings automatically.
 811    
 812        label_pos : float, optional
 813            Position along the edge to place the label (0=head, 1=tail).
 814            Default is 0.75 (closer to the target node).
 815    
 816        font_size : int, optional
 817            Font size of the edge labels. Default is 10.
 818    
 819        font_color : str, optional
 820            Font color for the edge labels. Default is 'k' (black).
 821    
 822        font_family : str, optional
 823            Font family for the edge labels. Default is 'sans-serif'.
 824    
 825        font_weight : str, optional
 826            Font weight for the edge labels (e.g., 'normal', 'bold'). Default is 'normal'.
 827    
 828        alpha : float, optional
 829            Opacity for the edge labels. If None, labels are fully opaque.
 830    
 831        bbox : dict, optional
 832            Matplotlib-style bbox dictionary to style the label background.
 833            Default is a white rounded box.
 834    
 835        horizontalalignment : str, optional
 836            Horizontal alignment of the text. Default is 'center'.
 837    
 838        verticalalignment : str, optional
 839            Vertical alignment of the text. Default is 'center'.
 840    
 841        ax : matplotlib.axes.Axes, optional
 842            Axes on which to draw the labels. If None, uses the current axes.
 843    
 844        rotate : bool, optional
 845            Whether to rotate the label to match the edge direction. Default is True.
 846    
 847        clip_on : bool, optional
 848            Whether to clip labels at the plot boundary. Default is True.
 849    
 850        rad : float, optional
 851            Curvature of the edge (used to determine Bézier control point).
 852            Positive values bend counterclockwise; negative values clockwise. Default is 0 (straight).
 853    
 854        Returns
 855        -------
 856        dict
 857            Dictionary mapping edge tuples `(u, v)` to the corresponding matplotlib Text objects.
 858    
 859        Notes
 860        -----
 861        - If `rotate=True`, label text is automatically aligned with the direction of the edge.
 862        - Curvature (`rad`) enables visualization of bidirectional transitions using bent edges.
 863        - This method is intended for internal use and supports bended edge label placement to match custom edge rendering.
 864    
 865        See Also
 866        --------
 867        draw
 868        draw_networkx
 869        draw_networkx_nodes
 870        draw_networkx_edges
 871        draw_networkx_labels
 872        """                
 873        if ax is None:
 874            ax = plt.gca()
 875        if edge_labels is None:
 876            labels = {(u, v): d for u, v, d in self.edges(data=True)}
 877        else:
 878            labels = edge_labels
 879        text_items = {}
 880        for (n1, n2), label in labels.items():
 881            (x1, y1) = pos[n1]
 882            (x2, y2) = pos[n2]
 883            (x, y) = (
 884                x1 * label_pos + x2 * (1.0 - label_pos),
 885                y1 * label_pos + y2 * (1.0 - label_pos),
 886            )
 887            pos_1 = ax.transData.transform(np.array(pos[n1]))
 888            pos_2 = ax.transData.transform(np.array(pos[n2]))
 889            #linear_mid = 0.5*pos_1 + 0.5*pos_2
 890            linear_mid = (1-label_pos)*pos_1 + label_pos*pos_2
 891            d_pos = pos_2 - pos_1
 892            rotation_matrix = np.array([(0,1), (-1,0)])
 893            ctrl_1 = linear_mid + rad*rotation_matrix@d_pos
 894            ctrl_mid_1 = 0.5*pos_1 + 0.5*ctrl_1
 895            ctrl_mid_2 = 0.5*pos_2 + 0.5*ctrl_1
 896            bezier_mid = 0.5*ctrl_mid_1 + 0.5*ctrl_mid_2
 897            (x, y) = ax.transData.inverted().transform(bezier_mid)
 898
 899            if rotate:
 900                # in degrees
 901                angle = np.arctan2(y2 - y1, x2 - x1) / (2.0 * np.pi) * 360
 902                # make label orientation "right-side-up"
 903                if angle > 90:
 904                    angle -= 180
 905                if angle < -90:
 906                    angle += 180
 907                # transform data coordinate angle to screen coordinate angle
 908                xy = np.array((x, y))
 909                trans_angle = ax.transData.transform_angles(
 910                    np.array((angle,)), xy.reshape((1, 2))
 911                )[0]
 912            else:
 913                trans_angle = 0.0
 914            # use default box of white with white border
 915            if bbox is None:
 916                bbox = dict(boxstyle="round", ec=(1.0, 1.0, 1.0), fc=(1.0, 1.0, 1.0))
 917            if not isinstance(label, str):
 918                label = str(label)  # this makes "1" and 1 labeled the same
 919
 920            t = ax.text(
 921                x, y, label,
 922                size=font_size, color=font_color, family=font_family, weight=font_weight,
 923                alpha=alpha,
 924                horizontalalignment=horizontalalignment, verticalalignment=verticalalignment,
 925                rotation=trans_angle, transform=ax.transData,
 926                bbox=bbox, zorder=1,clip_on=clip_on)
 927            text_items[(n1, n2)] = t
 928
 929        ax.tick_params(axis="both", which="both",bottom=False, left=False, labelbottom=False, labelleft=False)
 930
 931        return text_items    
 932
 933    
 934    def drawTransitionGraph(self, pos=None, 
 935                        bended=False, node_size=1000, num=1, rad=-0.2,
 936                        edge_labels=None, node_shapes='o', figsize=(14, 7),
 937                        clear=True, fontsize=10,
 938                        fontcolor='black', label_pos=0.75):        
 939        """
 940        Visualize the state transition graph with nodes and directed edges.
 941    
 942        This method draws the graph using `matplotlib` and `networkx`, including labels,
 943        node colors, and optional curved (bended) edges for better visualization of
 944        bidirectional transitions. It supports layout customization and is useful for
 945        understanding the structure of the Markov model.
 946    
 947        Parameters
 948        ----------
 949        - pos : dict, optional
 950            A dictionary mapping nodes to positions. If None, a Kamada-Kawai layout is used.
 951    
 952        - bended : bool, optional
 953            If True, edges are drawn with curvature using arcs. Useful for distinguishing
 954            bidirectional transitions. Default is False.
 955    
 956        - node_size : int, optional
 957            Size of the nodes in the plot. Default is 1000.
 958    
 959        - num : int, optional
 960            Figure number for matplotlib (useful for managing multiple figures). Default is 1.
 961    
 962        - rad : float, optional
 963            The curvature radius for bended edges. Only used if `bended=True`. Default is -0.2.
 964    
 965        - edge_labels : dict, optional
 966            Optional dictionary mapping edges `(u, v)` to labels. If None, the `'label'`
 967            attribute from each edge is used.
 968    
 969        - node_shapes : str, optional
 970            Shape of the nodes (e.g., `'o'` for circle, `'s'` for square). Default is `'o'`.
 971    
 972        - figsize : tuple, optional
 973            Size of the matplotlib figure. Default is (14, 7).
 974    
 975        - clear : bool, optional
 976            If True, clears the figure before plotting. Default is True.
 977    
 978        - fontsize : int, optional
 979            Font size used for node labels. Default is 10.
 980    
 981        - fontcolor : str, optional
 982            Font color for node labels. Default is `'black'`.
 983    
 984        - label_pos : float, optional
 985            Position of edge labels along the edge (0=start, 1=end). Default is 0.75.
 986    
 987        Returns
 988        -------
 989        - fig : matplotlib.figure.Figure
 990            The created matplotlib figure object.
 991    
 992        - ax : matplotlib.axes.Axes
 993            The axes object where the graph is drawn.
 994    
 995        - pos : dict
 996            The dictionary of node positions used for drawing.
 997    
 998        Notes
 999        -----
1000        - Node labels and colors must be set beforehand using `addState()` or `setStateProperties()`.
1001        - Edge attributes must include `'label'` and `'color'`.
1002        - This method modifies the current matplotlib figure and is intended for interactive or inline use.
1003        """        
1004        node_labels = {node: data["label"] for node, data in self.nodes(data=True)}        
1005        node_colors = {node: data["color"] for node, data in self.nodes(data=True)}        
1006        node_colors = [data["color"] for node, data in self.nodes(data=True)]
1007        #node_shapes = {(hw,sw): "o" if hw+sw<G.N else "s" for hw,sw in G.nodes()}        
1008            
1009        edge_cols = {(u,v): self[u][v]["color"] for u,v in self.edges}
1010            
1011        if edge_labels is None:
1012            edge_labels = {(u,v): self[u][v]["label"] for u,v in self.edges}
1013            
1014        if pos is None:
1015            pos = nx.kamada_kawai_layout(self)    
1016            
1017        plt.figure(figsize=figsize, num=num, clear=clear)                        
1018        nx.draw_networkx_nodes(self, pos, node_color=node_colors, node_shape = node_shapes, node_size=node_size)                    
1019        
1020        nx.draw_networkx_labels(self, pos, labels=node_labels, font_size=fontsize, font_color=fontcolor)  # Draw node labels
1021        if bended: 
1022            nx.draw_networkx_edges(self, pos,  width=1, edge_color=[edge_cols[edge] for edge in self.edges], 
1023                                   node_size = node_size, 
1024                                   arrows = True, arrowstyle = '-|>',
1025                                   connectionstyle=f"arc3,rad={rad}")        
1026            self._draw_networkx_edge_labels(pos, ax=plt.gca(), edge_labels=edge_labels,rotate=False, rad = rad, label_pos=label_pos)
1027        else:
1028            nx.draw_networkx_edges(self, pos,  width=1, edge_color=[edge_cols[edge] for edge in self.edges], 
1029                               node_size = node_size, 
1030                               #min_target_margin=17, arrowsize=15,
1031                               arrows = True, arrowstyle = '-|>')                           
1032            nx.draw_networkx_edge_labels(self, pos, edge_labels, label_pos=label_pos) #, verticalalignment='center',)        
1033        plt.axis('off');    
1034        
1035        return plt.gcf(), plt.gca(), pos
1036    
1037    def animateSimulation(self, expectedTimePerState=1, inNotebook=False, **kwargs):
1038        """
1039        Animate a simulation trajectory through the state transition graph.
1040        The animation of the CTMC simulation run either in a Jupyter notebook or as a regular script.
1041    
1042        This method visualizes the path of a simulated or predefined sequence of states by
1043        temporarily highlighting each visited state over time. The time spent in each state
1044        is scaled relative to the average sojourn time to produce a visually smooth animation.
1045    
1046        Parameters
1047        ----------
1048        - **expectedTimePerState** : float, optional
1049            Approximate duration (in seconds) for an average state visit in the animation.
1050            The actual pause duration for each state is scaled proportionally to its sojourn time.
1051            Default is 1 second.
1052            
1053        - **inNotebook** : bool, optional
1054            If True, uses Jupyter-compatible animation (with display + clear_output).
1055            If False, uses standard matplotlib interactive animation. Default is True.
1056    
1057        - `**kwargs` : dict
1058            Additional keyword arguments including:
1059                - states : list
1060                    A list of visited states (nodes) to animate.
1061                - times : list or np.ndarray
1062                    Sojourn times in each state (same length as `states`).
1063                - startNode : hashable
1064                    Optional start node for automatic simulation if `states` and `times` are not provided.
1065                - numEvents : int
1066                    Number of events (state transitions) to simulate.
1067                - Additional drawing-related kwargs are passed to `drawTransitionGraph()`.
1068    
1069        Raises
1070        ------
1071        ValueError
1072            If neither a `(states, times)` pair nor `(startNode, numEvents)` are provided.
1073    
1074        Returns
1075        -------
1076        None
1077            The animation is shown interactively using matplotlib but no data is returned.
1078    
1079        Notes
1080        -----
1081        - If `states` and `times` are not supplied, a simulation is run via `simulateMarkovModel()`.
1082        - The function uses matplotlib's interactive mode (`plt.ion()`) to animate transitions.
1083        - The average sojourn time across all states is used to normalize animation speed.
1084        - Each visited state is highlighted in red for its corresponding (scaled) dwell time.
1085        - Drawing arguments like layout, font size, or color can be customized via kwargs.
1086        """        
1087        # expTimePerState: in seconds
1088        if "states" in kwargs and "times" in kwargs:
1089            states = kwargs["states"]
1090            times = kwargs["times"]
1091        else:
1092            # Run default simulation if data not provided
1093            if "startNode" in kwargs and "numEvents" in kwargs:
1094                states, times = self.simulateMarkovModel(startNode=None, numEvents=10)
1095            else:
1096                raise ValueError("Must provide either ('states' and 'times') or ('startNode' and 'numEvents').")
1097        
1098        #avg_OutGoingRate = np.mean([self.nodes[n].get("sumOutgoingRates", 0) for n in self.nodes])
1099        avg_Time = np.mean([1.0/self.nodes[n].get("sumOutgoingRates", 0) for n in self.nodes])
1100        
1101            
1102        # Remove simulation-related keys before drawing
1103        draw_kwargs = {k: v for k, v in kwargs.items() if k not in {"states", "times", "startNode", "numEvents"}}
1104                
1105        fig, ax, pos = self.drawTransitionGraph(**draw_kwargs)
1106        plt.tight_layout()
1107        
1108        if not inNotebook:
1109            plt.ion()
1110                                
1111        for node, time in zip(states,times):
1112            if inNotebook:
1113                clear_output(wait=True)
1114                artist = nx.draw_networkx_nodes(self, pos, nodelist=[node], ax=ax,
1115                                                node_color=self._const.COLOR_NODE_ANIMATION_HIGHLIGHT, node_size=1000)
1116                display(fig)
1117            else:
1118                artist = nx.draw_networkx_nodes(self, pos, nodelist=[node],ax=ax,
1119                                       node_color=self._const.COLOR_NODE_ANIMATION_HIGHLIGHT, node_size=1000)
1120                plt.draw()                        
1121            plt.pause(time/avg_Time*expectedTimePerState)
1122            
1123            # Remove highlighted node    
1124            artist.remove()
1125        
1126        pass
1127    
1128        
1129    
1130    def states(self, data=False):
1131        """
1132        Return a view of the graph's states (nodes), optionally with attributes.
1133    
1134        This method is functionally identical to `self.nodes()` in NetworkX, but is
1135        renamed to `states()` to reflect the semantics of a state transition graph or
1136        Markov model, where nodes represent system states.
1137    
1138        Parameters
1139        ----------
1140        **data** : bool, optional
1141            If True, returns a view of `(node, attribute_dict)` pairs.
1142            If False (default), returns a view of node identifiers only.
1143    
1144        Returns
1145        -------
1146        *networkx.classes.reportviews.NodeView* or NodeDataView
1147            A view over the graph’s nodes or `(node, data)` pairs, depending on the `data` flag.
1148    
1149        Examples
1150        --------
1151        >>> G.states()
1152        ['A', 'B', 'C']
1153    
1154        >>> list(G.states(data=True))
1155        [('A', {'color': 'red', 'label': 'Start'}), ('B', {...}), ...]
1156    
1157        Notes
1158        -----
1159        - This is a convenience wrapper for semantic clarity in models where nodes represent states.
1160        - The view is dynamic: changes to the graph are reflected in the returned view.
1161        """
1162        return self.nodes(data=data)
1163
1164    
1165    def setStateColor(self, state, color):
1166        """
1167        Set the color attribute of a specific state (node).
1168    
1169        Parameters
1170        ----------
1171        **state** : hashable
1172            The identifier of the state whose color is to be updated.
1173    
1174        **color** : str
1175            The new color to assign to the state (used for visualization).
1176    
1177        Returns
1178        -------
1179        None
1180        """
1181        self.nodes[state]["color"] = color
1182        
1183        
1184    def setStateLabel(self, state, label):
1185        """
1186        Set the label attribute of a specific state (node).
1187    
1188        Parameters
1189        ----------
1190        **state** : hashable
1191            The identifier of the state whose label is to be updated.
1192    
1193        **label** : str
1194            The new label to assign to the state (used for visualization or annotation).
1195    
1196        Returns
1197        -------
1198        None
1199        """
1200        self.nodes[state]["label"] = label
1201
1202    def prob(self, state):
1203        """
1204        Return the steady-state probability of a specific state.
1205    
1206        Parameters
1207        ----------
1208        **state** : hashable
1209            The identifier of the state whose probability is to be retrieved.
1210    
1211        Returns
1212        -------
1213        **float**
1214            The steady-state probability of the specified state.
1215    
1216        Raises
1217        ------
1218        KeyError
1219            If the state is not found in the steady-state probability dictionary.
1220        """
1221        return self.state_probabilities[state]
1222        
1223    
1224    def probs(self):
1225        """
1226        Return the steady-state probabilities for all states.
1227    
1228        Returns
1229        -------
1230        **dict**
1231            A dictionary mapping each state to its steady-state probability.
1232    
1233        Notes
1234        -----
1235        - The probabilities are computed via `calcSteadyStateProbs()` and stored internally.
1236        - Use this method to access the full steady-state distribution.
1237        """
1238        return self.state_probabilities
1239
1240    def getEmptySystemProbs(self, state=None):
1241        """
1242        Create an initial state distribution where all probability mass is in a single state.
1243    
1244        By default, this method returns a distribution where the entire probability mass
1245        is placed on the first node in the graph. Alternatively, a specific starting state
1246        can be provided.
1247    
1248        Parameters
1249        ----------
1250        **state** : hashable, optional
1251            The state to initialize with probability 1. If None, the first node
1252            in the graph (based on insertion order) is used.
1253    
1254        Returns
1255        -------
1256        **dict**
1257            A dictionary representing the initial distribution over states.
1258            The selected state has probability 1, and all others have 0.
1259    
1260        Examples
1261        --------
1262        >>> X0 = G.getEmptySystem()
1263        >>> sum(X0.values())  # 1.0
1264    
1265        Notes
1266        -----
1267        - This method is useful for setting up deterministic initial conditions
1268          in transient simulations of a CTMC.
1269        - The return format is compatible with methods like `transientProbs()`.
1270        """
1271        if state is None:
1272            first_node = next(iter(self.nodes))
1273        X0 = {s: 1 if s == first_node else 0 for s in self}
1274        return X0
1275
1276        
1277
1278    def transientProbs(self, initialDistribution, t):
1279        """
1280        Compute the transient state probabilities at time `t` for the CTMC.
1281    
1282        This method solves the forward equation of a continuous-time Markov chain (CTMC):
1283        X(t) = X0 · expm(Q·t), where `Q` is the generator (rate) matrix and `X0` is the
1284        initial state distribution. It returns a dictionary mapping each state to its
1285        probability at time `t`.
1286    
1287        Parameters
1288        ----------
1289        - **initialDistribution** : dict 
1290            A dictionary mapping states (nodes) to their initial probabilities at time `t=0`.
1291            The keys must match the nodes in the graph, and the values should sum to 1.    
1292            
1293        - **t** : float
1294            The time at which the transient distribution is to be evaluated.
1295    
1296        Returns
1297        -------
1298        **dict**
1299            A dictionary mapping each state (node) to its probability at time `t`.
1300    
1301        Notes
1302        -----
1303        - The rate matrix `Q` is generated automatically if not yet computed.
1304        - The state order in the rate matrix corresponds to the internal `_n2i` mapping.
1305        - Internally, the dictionary `initialDistribution` is converted into a vector
1306          aligned with the state indexing.
1307        - The transient solution is computed using `scipy.linalg.expm` for matrix exponentiation.
1308    
1309        Raises
1310        ------
1311        ValueError
1312            If the input dictionary has missing states or probabilities that do not sum to 1.
1313            (This is not enforced but recommended for correctness.)
1314        """        
1315        if self._rate_matrix is None:  self.createRateMatrix()
1316        Q = self._rate_matrix
1317        X0 = np.zeros(len(initialDistribution))
1318        
1319        for i,n in enumerate(self.nodes):
1320            X0[ self._n2i[n] ] = initialDistribution[n]
1321                
1322        
1323        X = X0 @ expm(Q*t)    
1324        
1325        matrix = {}
1326        for n in self.nodes:
1327            matrix[n] = X[ self._n2i[n] ]
1328                
1329        return matrix
1330        
1331    def symSolveMarkovModel(self):
1332        """
1333        Symbolically solve the global balance equations of the CTMC.
1334    
1335        This method constructs and solves the system of symbolic equations for the
1336        steady-state probabilities of each state in the CTMC. It uses symbolic transition
1337        rates stored on the edges to form balance equations for each node, along with
1338        the normalization constraint that all probabilities sum to 1.
1339    
1340        Returns
1341        -------
1342        - **solution** : dict
1343            A symbolic solution dictionary mapping SymPy symbols (e.g., x_{A}) to
1344            expressions in terms of model parameters.
1345    
1346        - **num_solution** : dict
1347            A numerical dictionary mapping each state (node) to its steady-state probability,
1348            computed by substituting the numeric values from `self._subs` into the symbolic solution.
1349    
1350        Notes
1351        -----
1352        - For each node `s`, a symbolic variable `x_{label}` is created using the node's label attribute.
1353        - One balance equation is created per state: total inflow = total outflow.
1354        - An additional constraint `sum(x_i) = 1` ensures proper normalization.
1355        - The symbolic system is solved with `sympy.solve`, and the results are simplified.
1356        - Numeric values are computed by substituting known parameter values (`self._subs`) into the symbolic solution.
1357    
1358        Examples
1359        --------
1360        >>> sym_sol, num_sol = G.symSolveMarkovModel()
1361        >>> sym_sol  # symbolic expressions for each state
1362        >>> num_sol  # evaluated numerical values for each state
1363        """
1364        X={}
1365        to_solve = []
1366        for s in self.nodes():
1367            X[s] = sp.Symbol(f'x_{{{self.nodes[s]["label"]}}}')
1368            to_solve.append(X[s])
1369        #%
1370        eqs = []
1371        for s in self.nodes():
1372             succ = self.successors(s) # raus
1373             out_rate = 0
1374             for z in list(succ):
1375                 #out_rate += X[s]*self.edges[s, z]["sym_rate"]                 
1376                 out_rate += X[s]*self[s][z]["sym_rate"]
1377             
1378             in_rate = 0
1379             for r in list(self.predecessors(s)):
1380                 #in_rate += X[r]*self.edges[r,s]["sym_rate"] 
1381                 in_rate += X[r]*self[r][s]["sym_rate"] 
1382             eqs.append( sp.Eq(out_rate, in_rate) )
1383        #% 
1384        Xsum = 0    
1385        for s in X:
1386            Xsum += X[s]
1387        
1388        eqs.append( sp.Eq(Xsum, 1) )    
1389        #%
1390        solution = sp.solve(eqs, to_solve)
1391        simplified_solution = sp.simplify(solution)
1392        #print(simplified_solution)
1393        num_solution = simplified_solution.subs(self._subs)
1394        #num_dict =  {str(var): str(sol) for var, sol in simplified_solution.items()}
1395        num_dict = {s: num_solution[X[s]] for s in X}                
1396        
1397        return simplified_solution, num_dict 
1398    
1399    def exportTransitionsToExcel(self, excel_path = "output.xlsx", num_colors=9, open_excel=False, verbose=True):
1400        """
1401        Export transition data from the state transition graph to an Excel file.
1402        
1403        This method collects all transitions (edges) from the graph and writes them
1404        to an Excel spreadsheet, optionally applying background colors to distinguish
1405        outgoing transitions from different source states. This export is helpful for
1406        manually verifying the transition structure of the model.
1407
1408        The transition type is defined by the `tid` field of each edge, and mapped to
1409        a human-readable string via the `EXPORT_NAME_TRANSITIONS` dictionary in the
1410        constants module (e.g., `Constants.EXPORT_NAME_TRANSITIONS[tid]`).
1411
1412        See Also
1413        --------
1414        StateTransitionGraph.addTransition : Method to add edges with symbolic rates and `tid`.
1415        StateTransitionGraph.__init__ : Initializes the graph and loads the constants module.
1416
1417        Parameters
1418        ----------
1419        - **excel_path** : str, optional
1420            Path to the output Excel file. Default is "output.xlsx".
1421        
1422        - num_colors : int, optional
1423            Number of unique background colors to cycle through for different states.
1424            Default is 9.
1425        
1426        - open_excel : bool, optional
1427            If True, automatically opens the generated Excel file after export.
1428            Default is False.
1429        
1430        - verbose : bool, optional
1431            If True, prints summary information about the export process.
1432            Default is True.
1433        
1434        Returns
1435        -------
1436        None
1437            The function writes data to disk and optionally opens Excel,
1438            but does not return a value.
1439
1440        Notes
1441        -----
1442        - Each edge must contain a `tid` attribute for transition type and a `label` for the rate.
1443        - Background coloring groups all outgoing transitions from the same source state.
1444        - Requires the constants module to define `EXPORT_NAME_TRANSITIONS`.
1445        """
1446
1447        rows = []
1448        for node in self.nodes():
1449            edges = list(self.edges(node))
1450            if verbose: print(f"Edges connected to node '{node}':")
1451            
1452            for i, (u, v) in enumerate(edges):
1453                latex = self[u][v]["label"]
1454                rate = sp.pretty(self[u][v]["sym_rate"])
1455                tid = self[u][v]["tid"]
1456                
1457                if tid: # 
1458                    if hasattr(self._const, 'EXPORT_NAME_TRANSITIONS'):
1459                        tid_name = self._const.EXPORT_NAME_TRANSITIONS[tid]
1460                    else:
1461                        tid_name  = ''                               
1462                
1463                rows.append({"from": u, "type": tid_name, "to": v, "rate": rate, "latex": latex})
1464            
1465        df = pd.DataFrame(rows, columns=["from", "type", "to", "rate", "latex"])
1466        df["checked"] = ""
1467        df.to_excel(excel_path, index=False)
1468        
1469        # Load workbook with openpyxl
1470        wb = load_workbook(excel_path)
1471        ws = wb.active
1472        
1473        # Define colors for different "from" values
1474        cmap = map(plt.cm.Pastel1, range(num_colors))
1475        # Convert RGBA to hex
1476        color_list = [mcolors.to_hex(color, keep_alpha=False).upper().replace("#", "") for color in cmap]
1477        
1478        # Assign color to each unique "from" value using modulo
1479        unique_from = (df["from"].unique())  # Sorted for consistent order
1480        from_color_map = {
1481            value: color_list[i % len(color_list)] for i, value in enumerate(unique_from)
1482        }
1483        
1484        # Apply coloring row by row
1485        for i,row in enumerate(ws.iter_rows(min_row=2, max_row=ws.max_row)):  # Skip header
1486            from_value = row[0].value  # "from" is first column
1487            #fill_color = color_list[i % len(color_list)]
1488            if isinstance(from_value, int):
1489                fill_color = from_color_map[from_value]
1490            else:
1491                fill_color = from_color_map.get(eval(from_value))
1492            
1493            if fill_color:
1494                fill = PatternFill(start_color=fill_color, end_color=fill_color, fill_type="solid")
1495                for cell in row:
1496                    cell.fill = fill
1497        
1498        # Save styled Excel file
1499        wb.save(excel_path)
1500        if verbose: print(f"File {excel_path} created.")
1501        if open_excel: 
1502            print(f"Opening the file {excel_path}.")
1503            os.startfile(excel_path)

Initialize the StateTransitionGraph for continuous time Markov chains (CTMC). Extends the networkx.Graph class to include additional methods and properties for the analysis, visualization, and simulation of CTMCs. The analysis allows deriving the state probabilities in the steady-state and in the transient phase.

This constructor sets up symbolic and numeric parameters for the state transition graph, imports a constants module dynamically, and initializes internal data structures for rate matrix computation and steady-state probability analysis.

Parameters

  • parameters : dict
    A dictionary mapping parameter names to their numerical values.
    Example: {'lambda': 1.0, 'mu': 1.0}. For each parameter key:

    • A symbolic variable is created using sympy.symbols. Can be accessed as G.sym("key") or G.sym_key
    • A variable is created which can be accessed via G.key (and the numerical value is assigned)
    • For the key="lambda", G.lam and G.sym_lambda are created

      The parameters can be accessed via property parameters

  • constants : str, optional The name of a Python module (as string) to import dynamically and assign to self.const. This can be used to store global or user-defined constants. Default is "Constants".
  • *args : tuple Additional positional arguments passed to the base networkx.DiGraph constructor.
  • **kwargs : dict Additional keyword arguments passed to the base networkx.DiGraph constructor.

Attributes

  • _rate_matrix : ndarray or None
    The continuous-time rate matrix Q of the CTMC. Computed when calling createRateMatrix().
  • _state_probabilities : dict or None
    Dictionary of steady-state probabilities, computed via calcSteadyStateProbs().
  • _n2i : dict or None
    Mapping from node identifiers to matrix row/column indices used for rate matrix construction.
  • _sym : dict
    Mapping of parameter names (strings) to their symbolic sympy representations.
  • _subs : dict
    Mapping from symbolic parameters to their numerical values (used for substitutions).
  • _parameters : dict
    Copy of the original parameter dictionary passed at initialization.
  • _transitionProbabilitesComputed : bool
    Internal flag indicating whether transition probabilities for the embedded Markov chain have been computed.
  • const : module
    Dynamically imported module (default "Constants"), used for storing user-defined constants such as colors or transition IDs.

Inherits

networkx.Graph: The base graph class from NetworkX.
StateTransitionGraph(parameters, constants='Constants', *args, **kwargs)
153    def __init__(self, parameters, constants="Constants", *args, **kwargs):
154        """
155        Initialize the state transition graph with symbolic parameters.
156
157        This constructor sets up the symbolic and numerical environment for modeling
158        a CTMC. Parameters are stored as both symbolic expressions and numerical values,
159        and convenient attribute access is provided for both forms.
160
161        Parameters
162        ----------
163        - **parameters** : dict  
164          A dictionary mapping parameter names to their numerical values.     
165            Example: {'lambda': 1.0, 'mu': 1.0}. 
166            For each parameter key:
167            - A symbolic variable is created using sympy.symbols. Can be accessed as `G.sym("key")` or `G.sym_key`
168            - A variable is created which can be accessed via `G.key` (and the numerical value is assigned)
169            - For the key="lambda", `G.lam` and `G.sym_lambda` are created
170        
171            The parameters can be accessed via property `parameters`
172        - **constants** : str, optional
173            The name of a Python module (as string) to import dynamically and assign to `self.const`.
174            This can be used to store global or user-defined constants.
175            Default is "Constants".
176        - `*args` : tuple
177            Additional positional arguments passed to the base `networkx.DiGraph` constructor.
178        - `**kwargs` : dict
179            Additional keyword arguments passed to the base `networkx.DiGraph` constructor.
180        """        
181        super().__init__(*args, **kwargs)
182        self._rate_matrix = None  # Rate Matrix Q of this transition graph
183        self._state_probabilities = None  # steady state proabilities of this CTMC
184        self._n2i = None # internal index for generating the rate matrix Q: dictionary (node to index)
185        
186        #parameters = {'lambda': 1.0, 'mu':1.0, 'n':5}
187        self._parameters = parameters
188        self._sym = {key: sp.symbols(key) for key in parameters}
189        self._subs = {} # contains symbolic parameter as key and provides numerical value for computation.
190        for k in parameters:
191            self._subs[ self._sym[k] ] = parameters[k] 
192
193        # access to parameters and symbols using dot-notation, e.g. G.lam, G.sym_lambda, or G.sym("lambda")
194        for key, value in parameters.items():
195            if key=='lambda':
196                setattr(self, 'lam', value)
197            else:
198                setattr(self, key, value)
199            
200            setattr(self, "sym_"+key, self._sym[key])
201            
202        self._transitionProbabilitesComputed = False
203        
204        self._const = importlib.import_module(constants)

Initialize the state transition graph with symbolic parameters.

This constructor sets up the symbolic and numerical environment for modeling a CTMC. Parameters are stored as both symbolic expressions and numerical values, and convenient attribute access is provided for both forms.

Parameters

  • parameters : dict
    A dictionary mapping parameter names to their numerical values.
    Example: {'lambda': 1.0, 'mu': 1.0}. For each parameter key:

    • A symbolic variable is created using sympy.symbols. Can be accessed as G.sym("key") or G.sym_key
    • A variable is created which can be accessed via G.key (and the numerical value is assigned)
    • For the key="lambda", G.lam and G.sym_lambda are created

      The parameters can be accessed via property parameters

  • constants : str, optional The name of a Python module (as string) to import dynamically and assign to self.const. This can be used to store global or user-defined constants. Default is "Constants".
  • *args : tuple Additional positional arguments passed to the base networkx.DiGraph constructor.
  • **kwargs : dict Additional keyword arguments passed to the base networkx.DiGraph constructor.
def addState(self, state, color=None, label=None):
208    def addState(self, state, color=None, label=None):
209        """
210        Add a state (node) to the graph with optional color and label attributes.
211    
212        If the node already exists, its attributes are updated with the provided values.
213    
214        Parameters
215        ----------
216        - **state** : hashable
217            The identifier for the state. Can be a string, int, or a tuple. If no label
218            is provided, the label will be auto-generated as:
219            - comma-separated string for tuples (e.g., `(1, 2)` -> `'1,2'`)
220            - stringified value for other types
221    
222        - **color** : str, optional
223            The color assigned to the state (used for visualization).
224            If not specified, defaults to `self.const.COLOR_NODE_DEFAULT`.
225    
226        - **label** : str, optional
227            A label for the node (used for display or annotation). If not provided,
228            it will be auto-generated from the state identifier.
229    
230        Notes
231        -----
232        - This method wraps `self.add_node()` from NetworkX.
233        - If the node already exists in the graph, this method does not raise an error;
234          it will update the node’s existing attributes.
235        """
236        if label is None:
237            if isinstance(state, tuple):
238                label = ','.join(map(str, state)) # if state is a tuple
239            else: 
240                label = str(state)
241            #self.nodes[state]["label"] = label
242        
243        if color is None:
244            color = self._const.COLOR_NODE_DEFAULT
245        
246        # no error if the node already exists.Attributes get updated if you pass them again.
247        self.add_node(state, color=color, label=label)

Add a state (node) to the graph with optional color and label attributes.

If the node already exists, its attributes are updated with the provided values.

Parameters

  • state : hashable The identifier for the state. Can be a string, int, or a tuple. If no label is provided, the label will be auto-generated as:

    • comma-separated string for tuples (e.g., (1, 2) -> '1,2')
    • stringified value for other types
  • color : str, optional The color assigned to the state (used for visualization). If not specified, defaults to self.const.COLOR_NODE_DEFAULT.

  • label : str, optional A label for the node (used for display or annotation). If not provided, it will be auto-generated from the state identifier.

Notes

  • This method wraps self.add_node() from NetworkX.
  • If the node already exists in the graph, this method does not raise an error; it will update the node’s existing attributes.
def setStateProperties(self, state, color=None, label=None, **kwargs):
249    def setStateProperties(self, state, color=None, label=None, **kwargs):
250        """
251        Update visual or metadata properties for an existing state (node) in the graph.
252    
253        This method updates the node's attributes such as color, label, and any additional
254        custom key-value pairs. The node must already exist in the graph.
255    
256        Parameters
257        ----------
258        - **state** : hashable
259            Identifier of the state to update. Must be a node already present in the graph.
260    
261        - **color** : str, optional
262            Color assigned to the node, typically used for visualization. If not provided,
263            the default node color (`self._const.COLOR_NODE_DEFAULT`) is used.
264    
265        - **label** : str, optional
266            Label for the node. If not provided, it is auto-generated from the state ID
267            (especially for tuples).
268    
269        - `**kwargs` : dict
270            Arbitrary additional key-value pairs to set as node attributes. These can include
271            any custom metadata.
272    
273        Raises
274        ------
275        KeyError
276            If the specified state is not present in the graph.
277    
278        Notes
279        -----
280        - This method ensures that required attributes (like `color` and `label`) are set for visualization. 
281        - Existing node attributes are updated or extended with `kwargs`.
282        """        
283        if state not in self:
284            raise KeyError(f"State '{state}' does not exist in the graph.")
285        
286        # mandatory for plotting
287        self.addState(state, color, label)
288        
289        # # Add any additional attributes
290        self.nodes[state].update(kwargs)

Update visual or metadata properties for an existing state (node) in the graph.

This method updates the node's attributes such as color, label, and any additional custom key-value pairs. The node must already exist in the graph.

Parameters

  • state : hashable Identifier of the state to update. Must be a node already present in the graph.

  • color : str, optional Color assigned to the node, typically used for visualization. If not provided, the default node color (self._const.COLOR_NODE_DEFAULT) is used.

  • label : str, optional Label for the node. If not provided, it is auto-generated from the state ID (especially for tuples).

  • **kwargs : dict Arbitrary additional key-value pairs to set as node attributes. These can include any custom metadata.

Raises

KeyError If the specified state is not present in the graph.

Notes

  • This method ensures that required attributes (like color and label) are set for visualization.
  • Existing node attributes are updated or extended with kwargs.
def setAllStateDefaultProperties(self, **kwargs):
292    def setAllStateDefaultProperties(self, **kwargs):
293        """
294        Set default visual and metadata properties for all states (nodes) in the graph.
295    
296        This method applies `setStateProperties()` to every node, ensuring each state
297        has at least the required default attributes (e.g., color and label).
298        Any additional keyword arguments are passed to `setStateProperties()`.
299    
300        Parameters
301        ----------
302        `**kwargs` : dict, optional
303            Optional keyword arguments to apply uniformly to all nodes,
304            such as color='gray', group='core', etc.
305    
306        Notes
307        -----
308        - This method is useful for initializing node attributes before plotting
309          or performing other graph-based computations.
310        - Nodes without existing attributes will receive defaults via `addState()`.
311        """
312        for state in self.nodes():
313            self.setStateProperties(state, **kwargs)

Set default visual and metadata properties for all states (nodes) in the graph.

This method applies setStateProperties() to every node, ensuring each state has at least the required default attributes (e.g., color and label). Any additional keyword arguments are passed to setStateProperties().

Parameters

**kwargs : dict, optional Optional keyword arguments to apply uniformly to all nodes, such as color='gray', group='core', etc.

Notes

  • This method is useful for initializing node attributes before plotting or performing other graph-based computations.
  • Nodes without existing attributes will receive defaults via addState().
def addTransition( self, origin_state, destination_state, sym_rate, tid=None, color=None, label=None):
317    def addTransition(self, origin_state, destination_state, sym_rate, tid = None, color=None, label=None):     
318        """
319        Add a directed transition (edge) between two states with a symbolic transition rate.
320    
321        This method adds a directed edge from `origin_state` to `destination_state`,
322        stores both the symbolic and numerical rate, and sets optional display properties
323        such as color, label, and transition ID. If the edge already exists, a warning is issued.
324    
325        Parameters
326        ----------
327        - **origin_state** : hashable
328            The source node (state) of the transition.
329    
330        - **destination_state** : hashable
331            The target node (state) of the transition.
332    
333        - **sym_rate** : sympy.Expr
334            A symbolic expression representing the transition rate.
335            This will be numerically evaluated using internal substitutions (`self._subs`).
336    
337        - **tid** : str or int, optional
338            A transition identifier (e.g. for grouping or styling).
339            If provided, the color will default to `self._const.COLOR_TRANSITIONS[tid]` if not set explicitly.
340    
341        - **color** : str, optional
342            Color for the edge (used in visualization). Defaults to the transition-specific color
343            if `tid` is provided.
344    
345        - **label** : str, optional
346            LaTeX-formatted string for labeling the edge. If not provided, a default label is
347            generated using `sympy.latex(sym_rate)`.
348    
349        Notes
350        -----
351        - The numeric rate is computed by substituting known parameter values into `sym_rate`.
352        - The following attributes are stored on the edge:
353            - `rate` (float): numeric value of the transition rate
354            - `sym_rate` (sympy.Expr): original symbolic expression
355            - `label` (str): display label
356            - `color` (str): edge color
357            - `tid` (str or int): transition ID
358        - If the edge already exists, a `UserWarning` is issued but the edge is still overwritten.
359        """        
360        if tid is not None:
361            if color is None:
362                color = self._const.COLOR_TRANSITIONS[ tid ]
363                
364        rate = sym_rate.subs( self._subs )
365        
366        if label is None: string = "$"+sp.latex(sym_rate)+"$"
367                
368        if self.has_edge(origin_state, destination_state):
369            theedge = self[origin_state][destination_state]
370            warningstr = f"Edge already exists: {origin_state} -> {destination_state} with rate {theedge['label']}"
371            warnings.warn(warningstr, category=UserWarning)
372        self.add_edge(origin_state, destination_state, rate=float(rate), label=string, 
373                   color=color, sym_rate=sym_rate, tid=tid)

Add a directed transition (edge) between two states with a symbolic transition rate.

This method adds a directed edge from origin_state to destination_state, stores both the symbolic and numerical rate, and sets optional display properties such as color, label, and transition ID. If the edge already exists, a warning is issued.

Parameters

  • origin_state : hashable The source node (state) of the transition.

  • destination_state : hashable The target node (state) of the transition.

  • sym_rate : sympy.Expr A symbolic expression representing the transition rate. This will be numerically evaluated using internal substitutions (self._subs).

  • tid : str or int, optional A transition identifier (e.g. for grouping or styling). If provided, the color will default to self._const.COLOR_TRANSITIONS[tid] if not set explicitly.

  • color : str, optional Color for the edge (used in visualization). Defaults to the transition-specific color if tid is provided.

  • label : str, optional LaTeX-formatted string for labeling the edge. If not provided, a default label is generated using sympy.latex(sym_rate).

Notes

  • The numeric rate is computed by substituting known parameter values into sym_rate.
  • The following attributes are stored on the edge:
    • rate (float): numeric value of the transition rate
    • sym_rate (sympy.Expr): original symbolic expression
    • label (str): display label
    • color (str): edge color
    • tid (str or int): transition ID
  • If the edge already exists, a UserWarning is issued but the edge is still overwritten.
def sym(self, key):
375    def sym(self, key):
376        """
377        Retrieve the symbolic representation of a model parameter.
378    
379        Parameters
380        ----------
381        **key** : str
382            The name of the parameter (e.g., 'lambda', 'mu', 'n').
383    
384        Returns
385        -------
386        **sympy.Symbol**
387            The symbolic variable corresponding to the given key.
388    
389        Raises
390        ------
391        KeyError
392            If the key does not exist in the internal symbolic mapping (`_sym`).
393    
394        Examples
395        --------
396        >>> G.sym('lambda')
397        lambda
398    
399        Notes
400        -----
401        - The symbolic mapping is initialized during object construction based on the `parameters` dictionary.
402        - To access the symbol using dot-notation, use `.sym_key` attributes created during initialization.
403        """
404        return self._sym[key]

Retrieve the symbolic representation of a model parameter.

Parameters

key : str The name of the parameter (e.g., 'lambda', 'mu', 'n').

Returns

sympy.Symbol The symbolic variable corresponding to the given key.

Raises

KeyError If the key does not exist in the internal symbolic mapping (_sym).

Examples

>>> G.sym('lambda')
lambda

Notes

  • The symbolic mapping is initialized during object construction based on the parameters dictionary.
  • To access the symbol using dot-notation, use .sym_key attributes created during initialization.
parameters
407    @property
408    def parameters(self):
409        """
410        Return the model parameters as a dictionary.
411    
412        This includes the original numerical values provided during initialization.
413    
414        Returns
415        -------
416        **dict**
417            Dictionary of parameter names and their corresponding numeric values.
418    
419        Examples
420        --------
421        >>> G.parameters
422        {'lambda': 1.0, 'mu': 0.5, 'n': 3}
423    
424        Notes
425        -----
426        - This property does not return symbolic values. Use `self.sym(key)` for symbolic access.
427        - Parameter symbols can also be accessed using attributes like `sym_lambda` or the `sym()` method.
428        """
429        return self._parameters

Return the model parameters as a dictionary.

This includes the original numerical values provided during initialization.

Returns

dict Dictionary of parameter names and their corresponding numeric values.

Examples

>>> G.parameters
{'lambda': 1.0, 'mu': 0.5, 'n': 3}

Notes

  • This property does not return symbolic values. Use self.sym(key) for symbolic access.
  • Parameter symbols can also be accessed using attributes like sym_lambda or the sym() method.
rate_matrix
433    @property
434    def rate_matrix(self):
435        """
436        Return the continuous-time rate matrix (Q matrix) of the state transition graph.
437    
438        This matrix defines the transition rates between all pairs of states in the
439        continuous-time Markov chain (CTMC) represented by this graph.
440    
441        Returns
442        -------
443        **numpy.ndarray** or None
444            The rate matrix `Q`, where `Q[i, j]` is the transition rate from state `i` to `j`.
445            Returns `None` if the rate matrix has not been computed yet.
446    
447        Notes
448        -----
449        - The matrix is typically generated via a method like `computeRateMatrix()` or similar.
450        - Internal state indexing is handled by `_n2i` (node-to-index mapping).
451        - The rate matrix is stored internally as `_rate_matrix`.
452        """
453        return self._rate_matrix

Return the continuous-time rate matrix (Q matrix) of the state transition graph.

This matrix defines the transition rates between all pairs of states in the continuous-time Markov chain (CTMC) represented by this graph.

Returns

numpy.ndarray or None The rate matrix Q, where Q[i, j] is the transition rate from state i to j. Returns None if the rate matrix has not been computed yet.

Notes

  • The matrix is typically generated via a method like computeRateMatrix() or similar.
  • Internal state indexing is handled by _n2i (node-to-index mapping).
  • The rate matrix is stored internally as _rate_matrix.
state_probabilities
457    @property
458    def state_probabilities(self):
459        """
460        Return the steady-state probabilities of the continuous-time Markov chain (CTMC).
461    
462        These probabilities represent the long-run proportion of time the system spends
463        in each state, assuming the CTMC is irreducible and positive recurrent.
464    
465        Returns
466        -------
467        **dict** or None
468            A dictionary mapping each state to its steady-state probability.
469            Returns `None` if the probabilities have not yet been computed.
470    
471        Notes
472        -----
473        - The probabilities are stored internally after being computed via a dedicated method
474          `calcSteadyStateProbs()`.
475        - The result is stored in the `_state_probabilities` attribute.
476        """
477        return self._state_probabilities

Return the steady-state probabilities of the continuous-time Markov chain (CTMC).

These probabilities represent the long-run proportion of time the system spends in each state, assuming the CTMC is irreducible and positive recurrent.

Returns

dict or None A dictionary mapping each state to its steady-state probability. Returns None if the probabilities have not yet been computed.

Notes

  • The probabilities are stored internally after being computed via a dedicated method calcSteadyStateProbs().
  • The result is stored in the _state_probabilities attribute.
def printEdges(self, state):
481    def printEdges(self, state):
482        """
483        Print all outgoing edges from a given state in the graph.
484    
485        For each edge from the specified `state`, this method prints the destination
486        state and the corresponding transition label (typically the symbolic rate).
487    
488        Parameters
489        ----------
490        **state** : hashable
491            The source node (state) from which outgoing edges will be printed.            
492    
493        Returns
494        -------
495        None
496            This method prints output directly to the console and does not return a value.
497    
498        Notes
499        -----
500        - The transition label is retrieved from the edge attribute `'label'`.
501        - Assumes `state` exists in the graph; otherwise, a KeyError will be raised.
502        """
503        print(f"from {state}")
504        node = self[state]
505        for e in node:
506            rate = node[e]['label']
507            print(f" -> {e} with {rate}")

Print all outgoing edges from a given state in the graph.

For each edge from the specified state, this method prints the destination state and the corresponding transition label (typically the symbolic rate).

Parameters

state : hashable The source node (state) from which outgoing edges will be printed.

Returns

None This method prints output directly to the console and does not return a value.

Notes

  • The transition label is retrieved from the edge attribute 'label'.
  • Assumes state exists in the graph; otherwise, a KeyError will be raised.
def createRateMatrix(self):
510    def createRateMatrix(self):
511        """
512        Construct the continuous-time Markov chain (CTMC) rate matrix Q from the graph.
513    
514        This method generates the generator matrix `Q` corresponding to the structure and
515        rates defined in the directed graph. Each node represents a system state, and each
516        directed edge must have a `'rate'` attribute representing the transition rate from
517        one state to another.
518    
519        The matrix `Q` is constructed such that:
520        - `Q[i, j]` is the transition rate from state `i` to state `j`
521        - Diagonal entries satisfy `Q[i, i] = -sum(Q[i, j]) for j ≠ i`, so that each row sums to 0
522    
523        The resulting matrix and the internal node-to-index mapping (`_n2i`) are stored
524        as attributes of the graph.
525    
526        Returns
527        -------
528        None
529            The method updates internal attributes:
530            - `_rate_matrix` : numpy.ndarray
531                The CTMC rate matrix Q
532            - `_n2i` : dict
533                Mapping from node identifiers to matrix row/column indices
534    
535        Notes
536        -----
537        - Assumes all edges in the graph have a `'rate'` attribute.
538        - The node order used in the matrix is consistent with the order from `list(self.nodes)`.
539        - Use `self.rate_matrix` to access the generated matrix after this method is called.
540        """
541        n2i = {}
542        nodes = list(self.nodes)
543        for i,node in enumerate(nodes):
544            n2i[node] = i         
545        Q = np.zeros((len(nodes),len(nodes))) # rate matrix
546        
547        for edge in self.edges:
548            i0 = n2i[edge[0]]
549            i1 = n2i[edge[1]]
550            Q[i0,i1] = self[edge[0]][edge[1]]["rate"] 
551            
552        np.fill_diagonal(Q, -Q.sum(axis=1))
553        self._rate_matrix = Q
554        self._n2i = n2i
555        #return Q, n2i

Construct the continuous-time Markov chain (CTMC) rate matrix Q from the graph.

This method generates the generator matrix Q corresponding to the structure and rates defined in the directed graph. Each node represents a system state, and each directed edge must have a 'rate' attribute representing the transition rate from one state to another.

The matrix Q is constructed such that:

  • Q[i, j] is the transition rate from state i to state j
  • Diagonal entries satisfy Q[i, i] = -sum(Q[i, j]) for j ≠ i, so that each row sums to 0

The resulting matrix and the internal node-to-index mapping (_n2i) are stored as attributes of the graph.

Returns

None The method updates internal attributes: - _rate_matrix : numpy.ndarray The CTMC rate matrix Q - _n2i : dict Mapping from node identifiers to matrix row/column indices

Notes

  • Assumes all edges in the graph have a 'rate' attribute.
  • The node order used in the matrix is consistent with the order from list(self.nodes).
  • Use self.rate_matrix to access the generated matrix after this method is called.
def calcSteadyStateProbs(self, verbose=False):
557    def calcSteadyStateProbs(self, verbose=False):        
558        """
559        Compute the steady-state probabilities of the continuous-time Markov chain (CTMC).
560    
561        This method calculates the stationary distribution of the CTMC defined by the rate
562        matrix of the state transition graph. The computation solves the linear system
563        `XQ = 0` with the constraint that the probabilities X sum to 1, using a modified
564        version of the rate matrix `Q`.
565    
566        Parameters
567        ----------
568        `verbose` : bool, optional
569            If True, prints the rate matrix, the modified matrix, the right-hand side vector,
570            and the resulting steady-state probabilities.
571    
572        Returns
573        -------
574        **probs** : dict
575            A dictionary mapping each state (node in the graph) to its steady-state probability.
576    
577        Notes
578        -----
579        - The graph must already contain valid transition rates on its edges.
580        - If the rate matrix has not yet been created, `createRateMatrix()` is called automatically.
581        - The method modifies the last column of the rate matrix to impose the normalization constraint
582          `sum(X) = 1`.
583        - The inverse of the modified matrix is computed directly, so the graph should be small enough
584          to avoid numerical instability or performance issues.
585        - The resulting probabilities are stored internally in `_state_probabilities` and can be accessed
586          via the `state_probabilities` property.
587        """
588        
589        # compute transition matrix if not already done
590        if self._rate_matrix is None:  self.createRateMatrix()
591        Q2 =  self._rate_matrix.copy()
592        if verbose:
593            print("Rate matrix Q")
594            print(f'Q=\n{Q2}\n')
595        
596        Q2[:, -1] = 1
597        if verbose:        
598            print(f'Matrix is changed to\nQ2=\n{Q2}\n')
599
600        b = np.zeros(len(Q2))
601        b[-1] = 1
602        if verbose:
603            print("Solve X = b @ Q2")
604            print(f'b=\n{b}\n')
605        
606        # state probabilities
607        X = b @ linalg.inv(Q2) # compute the matrix inverse
608            
609        # Generate a matrix with P(X,Y)
610        matrix = {}
611        for n in self.nodes:
612            matrix[n] = X[ self._n2i[n] ]
613        
614        if verbose:
615            print("Steady-state probabilities X")
616            print(f'X=\n{matrix}\n')
617        self._state_probabilities = matrix
618        return matrix

Compute the steady-state probabilities of the continuous-time Markov chain (CTMC).

This method calculates the stationary distribution of the CTMC defined by the rate matrix of the state transition graph. The computation solves the linear system XQ = 0 with the constraint that the probabilities X sum to 1, using a modified version of the rate matrix Q.

Parameters

verbose : bool, optional If True, prints the rate matrix, the modified matrix, the right-hand side vector, and the resulting steady-state probabilities.

Returns

probs : dict A dictionary mapping each state (node in the graph) to its steady-state probability.

Notes

  • The graph must already contain valid transition rates on its edges.
  • If the rate matrix has not yet been created, createRateMatrix() is called automatically.
  • The method modifies the last column of the rate matrix to impose the normalization constraint sum(X) = 1.
  • The inverse of the modified matrix is computed directly, so the graph should be small enough to avoid numerical instability or performance issues.
  • The resulting probabilities are stored internally in _state_probabilities and can be accessed via the state_probabilities property.
def calcTransitionProbabilities(self):
620    def calcTransitionProbabilities(self):
621        """
622        Compute and store transition probabilities for each node in the graph of the embedded discrete Markov chain.
623    
624        This method processes each node in the graph, treating outgoing edges as
625        transitions in an embedded discrete-time Markov chain derived from a continuous-time
626        Markov chain (CTMC). Transition *rates* on outgoing edges are normalized into
627        *probabilities*, and relevant data is stored in the node attributes.
628    
629        For each node, the following attributes are set:
630        - `"transitionProbs"` : np.ndarray
631            Array of transition probabilities to successor nodes.
632        - `"transitionNodes"` : list
633            List of successor node identifiers (ordered as in the probability array).
634        - `"sumOutgoingRates"` : float
635            Sum of all outgoing transition rates. Used to model the exponential waiting time
636            in continuous-time simulation (`T ~ Exp(sumOutgoingRates)`).
637    
638        Notes
639        -----
640        - The graph must be a directed graph where each edge has a `'rate'` attribute.
641        - Nodes with no outgoing edges are considered absorbing states.
642        - Useful for discrete-event simulation and Markov chain Monte Carlo (MCMC) modeling.
643        - Sets an internal flag `_transitionProbabilitesComputed = True` upon completion.
644    
645        Raises
646        ------
647        KeyError
648            If any outgoing edge lacks a `'rate'` attribute.
649        """
650        for node in self:
651            successors = list(self.successors(node))
652            rates = np.array([self[node][nbr]['rate'] for nbr in successors])
653            probs = rates / rates.sum()  # Normalize to probabilities
654            self.nodes[node]["transitionProbs"] = probs # probabilities
655            self.nodes[node]["transitionNodes"] = successors # list of nodes
656            self.nodes[node]["sumOutgoingRates"] = rates.sum() # time in state ~ EXP(sumOutgoingRates)
657            numel = len(successors)
658            if numel==0: 
659                print(f"Absorbing state: {node}")
660        self._transitionProbabilitesComputed = True

Compute and store transition probabilities for each node in the graph of the embedded discrete Markov chain.

This method processes each node in the graph, treating outgoing edges as transitions in an embedded discrete-time Markov chain derived from a continuous-time Markov chain (CTMC). Transition rates on outgoing edges are normalized into probabilities, and relevant data is stored in the node attributes.

For each node, the following attributes are set:

  • "transitionProbs" : np.ndarray Array of transition probabilities to successor nodes.
  • "transitionNodes" : list List of successor node identifiers (ordered as in the probability array).
  • "sumOutgoingRates" : float Sum of all outgoing transition rates. Used to model the exponential waiting time in continuous-time simulation (T ~ Exp(sumOutgoingRates)).

Notes

  • The graph must be a directed graph where each edge has a 'rate' attribute.
  • Nodes with no outgoing edges are considered absorbing states.
  • Useful for discrete-event simulation and Markov chain Monte Carlo (MCMC) modeling.
  • Sets an internal flag _transitionProbabilitesComputed = True upon completion.

Raises

KeyError If any outgoing edge lacks a 'rate' attribute.

def simulateMarkovModel(self, startNode=None, numEvents=10):
662    def simulateMarkovModel(self, startNode=None, numEvents=10):
663        """
664        Simulate a trajectory through the state transition graph as a continuous-time Markov chain (CTMC).
665        
666        This method performs a discrete-event simulation of the CTMC by generating a random walk
667        through the graph, starting at `startNode` (or a default node if not specified). Transitions
668        are chosen based on the normalized transition probabilities computed from edge rates, and
669        dwell times are sampled from exponential distributions based on the total outgoing rate
670        from each state.
671        
672        Parameters
673        ----------
674        - **startNode** : hashable, optional
675            The node at which the simulation starts. If not specified, the first node in the graph
676            (based on insertion order) is used.
677        
678        - **numEvents** : int, optional
679            The number of transition events (steps) to simulate. Default is 10.
680        
681        Returns
682        -------
683        **states** : list
684            A list of visited states, representing the sequence of nodes traversed in the simulation.
685        
686        **times** : numpy.ndarray
687            An array of dwell times in each state, sampled from exponential distributions
688            with rate equal to the sum of outgoing transition rates from each state.
689        
690        Notes
691        -----
692        - This method automatically calls `calcTransitionProbabilities()` if transition probabilities
693          have not yet been computed.
694        - If an absorbing state (i.e., a node with no outgoing edges) is reached, the simulation stops early.
695        - The cumulative time can be obtained by `np.cumsum(times)`.
696        - State durations follow `Exp(sumOutgoingRates)` distributions per CTMC theory.
697        """
698        if not self._transitionProbabilitesComputed: self.calcTransitionProbabilities()
699        
700        # simulate random walk through Graph = Markov Chain
701        if startNode is None:
702            startNode = next(iter(self))
703        
704
705        states = [] # list of states
706        times = np.zeros(numEvents) # time per state
707
708        node = startNode
709        for i in range(numEvents):
710            outgoingRate = self.nodes[node]["sumOutgoingRates"]  # time in state ~ EXP(sumOutgoingRates)
711            times[i] = random.expovariate(outgoingRate) # exponentially distributed time in state
712            states.append(node)
713            
714            # get next node    
715            probs = self.nodes[node]["transitionProbs"] # probabilities
716            successors = self.nodes[node]["transitionNodes"] # list of nodes        
717            numel = len(successors)
718            if numel==0: 
719                print(f"Absorbing state: {node}")
720                break
721            elif numel==1:
722                nextNode = successors[0]
723            else:
724                nextNode = successors[ np.random.choice(numel, p=probs) ]
725            # edge = G[node][nextNode]
726            node = nextNode
727            
728        return states, times

Simulate a trajectory through the state transition graph as a continuous-time Markov chain (CTMC).

This method performs a discrete-event simulation of the CTMC by generating a random walk through the graph, starting at startNode (or a default node if not specified). Transitions are chosen based on the normalized transition probabilities computed from edge rates, and dwell times are sampled from exponential distributions based on the total outgoing rate from each state.

Parameters

  • startNode : hashable, optional The node at which the simulation starts. If not specified, the first node in the graph (based on insertion order) is used.

  • numEvents : int, optional The number of transition events (steps) to simulate. Default is 10.

Returns

states : list A list of visited states, representing the sequence of nodes traversed in the simulation.

times : numpy.ndarray An array of dwell times in each state, sampled from exponential distributions with rate equal to the sum of outgoing transition rates from each state.

Notes

  • This method automatically calls calcTransitionProbabilities() if transition probabilities have not yet been computed.
  • If an absorbing state (i.e., a node with no outgoing edges) is reached, the simulation stops early.
  • The cumulative time can be obtained by np.cumsum(times).
  • State durations follow Exp(sumOutgoingRates) distributions per CTMC theory.
def calcProbabilitiesSimulationRun(self, states, times):
730    def calcProbabilitiesSimulationRun(self, states, times):
731        """
732        Estimate empirical state probabilities from a simulation run of the CTMC.
733    
734        This method computes the fraction of total time spent in each state during a simulation.
735        It uses the sequence of visited states and corresponding sojourn times to estimate the
736        empirical (time-weighted) probability distribution over states.
737    
738        Parameters
739        ----------
740        - **states** : list of tuple
741            A list of state identifiers visited during the simulation.    
742        - **times** : numpy.ndarray
743            An array of sojourn times corresponding to each visited state, same length as `states`.
744    
745        Returns
746        -------
747        **simProbs** : dict
748            A dictionary mapping each unique state to its estimated empirical probability,
749            based on total time spent in that state relative to the simulation duration.
750    
751        Notes
752        -----
753        - Internally, states are mapped to flat indices for efficient counting using NumPy.
754        - This method assumes that `states` contains only valid node identifiers from the graph.
755        - If `_n2i` (node-to-index mapping) has not been initialized, it will be computed here.
756        - This is especially useful for validating analytical steady-state probabilities
757          against sampled simulation results. Or in case of state-space explosions.
758        """        
759        if self._n2i is None:
760            n2i = {}
761            nodes = list(self.nodes)
762            for i,node in enumerate(nodes):
763                n2i[node] = i     
764            self._n2i = n2i
765                
766        
767        # Map tuples to flat indices
768        data = np.array(states)
769        max_vals = data.max(axis=0) + 1
770        indices = np.ravel_multi_index(data.T, dims=max_vals)
771
772        # Weighted bincount
773        bincount = np.bincount(indices, weights=times)
774
775        # Decode back to tuples
776        tuple_indices = np.array(np.unravel_index(np.nonzero(bincount)[0], shape=max_vals)).T
777        timePerState = bincount[bincount > 0]
778
779        simProbs = {}
780        totalTime = times.sum()
781        for t, c in zip(tuple_indices, timePerState):
782            simProbs[tuple(t)] = c / totalTime
783            #print(f"Tuple {tuple(t)} → weighted count: {c}")
784        return simProbs

Estimate empirical state probabilities from a simulation run of the CTMC.

This method computes the fraction of total time spent in each state during a simulation. It uses the sequence of visited states and corresponding sojourn times to estimate the empirical (time-weighted) probability distribution over states.

Parameters

  • states : list of tuple A list of state identifiers visited during the simulation.
  • times : numpy.ndarray An array of sojourn times corresponding to each visited state, same length as states.

Returns

simProbs : dict A dictionary mapping each unique state to its estimated empirical probability, based on total time spent in that state relative to the simulation duration.

Notes

  • Internally, states are mapped to flat indices for efficient counting using NumPy.
  • This method assumes that states contains only valid node identifiers from the graph.
  • If _n2i (node-to-index mapping) has not been initialized, it will be computed here.
  • This is especially useful for validating analytical steady-state probabilities against sampled simulation results. Or in case of state-space explosions.
def drawTransitionGraph( self, pos=None, bended=False, node_size=1000, num=1, rad=-0.2, edge_labels=None, node_shapes='o', figsize=(14, 7), clear=True, fontsize=10, fontcolor='black', label_pos=0.75):
 934    def drawTransitionGraph(self, pos=None, 
 935                        bended=False, node_size=1000, num=1, rad=-0.2,
 936                        edge_labels=None, node_shapes='o', figsize=(14, 7),
 937                        clear=True, fontsize=10,
 938                        fontcolor='black', label_pos=0.75):        
 939        """
 940        Visualize the state transition graph with nodes and directed edges.
 941    
 942        This method draws the graph using `matplotlib` and `networkx`, including labels,
 943        node colors, and optional curved (bended) edges for better visualization of
 944        bidirectional transitions. It supports layout customization and is useful for
 945        understanding the structure of the Markov model.
 946    
 947        Parameters
 948        ----------
 949        - pos : dict, optional
 950            A dictionary mapping nodes to positions. If None, a Kamada-Kawai layout is used.
 951    
 952        - bended : bool, optional
 953            If True, edges are drawn with curvature using arcs. Useful for distinguishing
 954            bidirectional transitions. Default is False.
 955    
 956        - node_size : int, optional
 957            Size of the nodes in the plot. Default is 1000.
 958    
 959        - num : int, optional
 960            Figure number for matplotlib (useful for managing multiple figures). Default is 1.
 961    
 962        - rad : float, optional
 963            The curvature radius for bended edges. Only used if `bended=True`. Default is -0.2.
 964    
 965        - edge_labels : dict, optional
 966            Optional dictionary mapping edges `(u, v)` to labels. If None, the `'label'`
 967            attribute from each edge is used.
 968    
 969        - node_shapes : str, optional
 970            Shape of the nodes (e.g., `'o'` for circle, `'s'` for square). Default is `'o'`.
 971    
 972        - figsize : tuple, optional
 973            Size of the matplotlib figure. Default is (14, 7).
 974    
 975        - clear : bool, optional
 976            If True, clears the figure before plotting. Default is True.
 977    
 978        - fontsize : int, optional
 979            Font size used for node labels. Default is 10.
 980    
 981        - fontcolor : str, optional
 982            Font color for node labels. Default is `'black'`.
 983    
 984        - label_pos : float, optional
 985            Position of edge labels along the edge (0=start, 1=end). Default is 0.75.
 986    
 987        Returns
 988        -------
 989        - fig : matplotlib.figure.Figure
 990            The created matplotlib figure object.
 991    
 992        - ax : matplotlib.axes.Axes
 993            The axes object where the graph is drawn.
 994    
 995        - pos : dict
 996            The dictionary of node positions used for drawing.
 997    
 998        Notes
 999        -----
1000        - Node labels and colors must be set beforehand using `addState()` or `setStateProperties()`.
1001        - Edge attributes must include `'label'` and `'color'`.
1002        - This method modifies the current matplotlib figure and is intended for interactive or inline use.
1003        """        
1004        node_labels = {node: data["label"] for node, data in self.nodes(data=True)}        
1005        node_colors = {node: data["color"] for node, data in self.nodes(data=True)}        
1006        node_colors = [data["color"] for node, data in self.nodes(data=True)]
1007        #node_shapes = {(hw,sw): "o" if hw+sw<G.N else "s" for hw,sw in G.nodes()}        
1008            
1009        edge_cols = {(u,v): self[u][v]["color"] for u,v in self.edges}
1010            
1011        if edge_labels is None:
1012            edge_labels = {(u,v): self[u][v]["label"] for u,v in self.edges}
1013            
1014        if pos is None:
1015            pos = nx.kamada_kawai_layout(self)    
1016            
1017        plt.figure(figsize=figsize, num=num, clear=clear)                        
1018        nx.draw_networkx_nodes(self, pos, node_color=node_colors, node_shape = node_shapes, node_size=node_size)                    
1019        
1020        nx.draw_networkx_labels(self, pos, labels=node_labels, font_size=fontsize, font_color=fontcolor)  # Draw node labels
1021        if bended: 
1022            nx.draw_networkx_edges(self, pos,  width=1, edge_color=[edge_cols[edge] for edge in self.edges], 
1023                                   node_size = node_size, 
1024                                   arrows = True, arrowstyle = '-|>',
1025                                   connectionstyle=f"arc3,rad={rad}")        
1026            self._draw_networkx_edge_labels(pos, ax=plt.gca(), edge_labels=edge_labels,rotate=False, rad = rad, label_pos=label_pos)
1027        else:
1028            nx.draw_networkx_edges(self, pos,  width=1, edge_color=[edge_cols[edge] for edge in self.edges], 
1029                               node_size = node_size, 
1030                               #min_target_margin=17, arrowsize=15,
1031                               arrows = True, arrowstyle = '-|>')                           
1032            nx.draw_networkx_edge_labels(self, pos, edge_labels, label_pos=label_pos) #, verticalalignment='center',)        
1033        plt.axis('off');    
1034        
1035        return plt.gcf(), plt.gca(), pos

Visualize the state transition graph with nodes and directed edges.

This method draws the graph using matplotlib and networkx, including labels, node colors, and optional curved (bended) edges for better visualization of bidirectional transitions. It supports layout customization and is useful for understanding the structure of the Markov model.

Parameters

  • pos : dict, optional A dictionary mapping nodes to positions. If None, a Kamada-Kawai layout is used.

  • bended : bool, optional If True, edges are drawn with curvature using arcs. Useful for distinguishing bidirectional transitions. Default is False.

  • node_size : int, optional Size of the nodes in the plot. Default is 1000.

  • num : int, optional Figure number for matplotlib (useful for managing multiple figures). Default is 1.

  • rad : float, optional The curvature radius for bended edges. Only used if bended=True. Default is -0.2.

  • edge_labels : dict, optional Optional dictionary mapping edges (u, v) to labels. If None, the 'label' attribute from each edge is used.

  • node_shapes : str, optional Shape of the nodes (e.g., 'o' for circle, 's' for square). Default is 'o'.

  • figsize : tuple, optional Size of the matplotlib figure. Default is (14, 7).

  • clear : bool, optional If True, clears the figure before plotting. Default is True.

  • fontsize : int, optional Font size used for node labels. Default is 10.

  • fontcolor : str, optional Font color for node labels. Default is 'black'.

  • label_pos : float, optional Position of edge labels along the edge (0=start, 1=end). Default is 0.75.

Returns

  • fig : matplotlib.figure.Figure The created matplotlib figure object.

  • ax : matplotlib.axes.Axes The axes object where the graph is drawn.

  • pos : dict The dictionary of node positions used for drawing.

Notes

  • Node labels and colors must be set beforehand using addState() or setStateProperties().
  • Edge attributes must include 'label' and 'color'.
  • This method modifies the current matplotlib figure and is intended for interactive or inline use.
def animateSimulation(self, expectedTimePerState=1, inNotebook=False, **kwargs):
1037    def animateSimulation(self, expectedTimePerState=1, inNotebook=False, **kwargs):
1038        """
1039        Animate a simulation trajectory through the state transition graph.
1040        The animation of the CTMC simulation run either in a Jupyter notebook or as a regular script.
1041    
1042        This method visualizes the path of a simulated or predefined sequence of states by
1043        temporarily highlighting each visited state over time. The time spent in each state
1044        is scaled relative to the average sojourn time to produce a visually smooth animation.
1045    
1046        Parameters
1047        ----------
1048        - **expectedTimePerState** : float, optional
1049            Approximate duration (in seconds) for an average state visit in the animation.
1050            The actual pause duration for each state is scaled proportionally to its sojourn time.
1051            Default is 1 second.
1052            
1053        - **inNotebook** : bool, optional
1054            If True, uses Jupyter-compatible animation (with display + clear_output).
1055            If False, uses standard matplotlib interactive animation. Default is True.
1056    
1057        - `**kwargs` : dict
1058            Additional keyword arguments including:
1059                - states : list
1060                    A list of visited states (nodes) to animate.
1061                - times : list or np.ndarray
1062                    Sojourn times in each state (same length as `states`).
1063                - startNode : hashable
1064                    Optional start node for automatic simulation if `states` and `times` are not provided.
1065                - numEvents : int
1066                    Number of events (state transitions) to simulate.
1067                - Additional drawing-related kwargs are passed to `drawTransitionGraph()`.
1068    
1069        Raises
1070        ------
1071        ValueError
1072            If neither a `(states, times)` pair nor `(startNode, numEvents)` are provided.
1073    
1074        Returns
1075        -------
1076        None
1077            The animation is shown interactively using matplotlib but no data is returned.
1078    
1079        Notes
1080        -----
1081        - If `states` and `times` are not supplied, a simulation is run via `simulateMarkovModel()`.
1082        - The function uses matplotlib's interactive mode (`plt.ion()`) to animate transitions.
1083        - The average sojourn time across all states is used to normalize animation speed.
1084        - Each visited state is highlighted in red for its corresponding (scaled) dwell time.
1085        - Drawing arguments like layout, font size, or color can be customized via kwargs.
1086        """        
1087        # expTimePerState: in seconds
1088        if "states" in kwargs and "times" in kwargs:
1089            states = kwargs["states"]
1090            times = kwargs["times"]
1091        else:
1092            # Run default simulation if data not provided
1093            if "startNode" in kwargs and "numEvents" in kwargs:
1094                states, times = self.simulateMarkovModel(startNode=None, numEvents=10)
1095            else:
1096                raise ValueError("Must provide either ('states' and 'times') or ('startNode' and 'numEvents').")
1097        
1098        #avg_OutGoingRate = np.mean([self.nodes[n].get("sumOutgoingRates", 0) for n in self.nodes])
1099        avg_Time = np.mean([1.0/self.nodes[n].get("sumOutgoingRates", 0) for n in self.nodes])
1100        
1101            
1102        # Remove simulation-related keys before drawing
1103        draw_kwargs = {k: v for k, v in kwargs.items() if k not in {"states", "times", "startNode", "numEvents"}}
1104                
1105        fig, ax, pos = self.drawTransitionGraph(**draw_kwargs)
1106        plt.tight_layout()
1107        
1108        if not inNotebook:
1109            plt.ion()
1110                                
1111        for node, time in zip(states,times):
1112            if inNotebook:
1113                clear_output(wait=True)
1114                artist = nx.draw_networkx_nodes(self, pos, nodelist=[node], ax=ax,
1115                                                node_color=self._const.COLOR_NODE_ANIMATION_HIGHLIGHT, node_size=1000)
1116                display(fig)
1117            else:
1118                artist = nx.draw_networkx_nodes(self, pos, nodelist=[node],ax=ax,
1119                                       node_color=self._const.COLOR_NODE_ANIMATION_HIGHLIGHT, node_size=1000)
1120                plt.draw()                        
1121            plt.pause(time/avg_Time*expectedTimePerState)
1122            
1123            # Remove highlighted node    
1124            artist.remove()
1125        
1126        pass

Animate a simulation trajectory through the state transition graph. The animation of the CTMC simulation run either in a Jupyter notebook or as a regular script.

This method visualizes the path of a simulated or predefined sequence of states by temporarily highlighting each visited state over time. The time spent in each state is scaled relative to the average sojourn time to produce a visually smooth animation.

Parameters

  • expectedTimePerState : float, optional Approximate duration (in seconds) for an average state visit in the animation. The actual pause duration for each state is scaled proportionally to its sojourn time. Default is 1 second.

  • inNotebook : bool, optional If True, uses Jupyter-compatible animation (with display + clear_output). If False, uses standard matplotlib interactive animation. Default is True.

  • **kwargs : dict Additional keyword arguments including: - states : list A list of visited states (nodes) to animate. - times : list or np.ndarray Sojourn times in each state (same length as states). - startNode : hashable Optional start node for automatic simulation if states and times are not provided. - numEvents : int Number of events (state transitions) to simulate. - Additional drawing-related kwargs are passed to drawTransitionGraph().

Raises

ValueError If neither a (states, times) pair nor (startNode, numEvents) are provided.

Returns

None The animation is shown interactively using matplotlib but no data is returned.

Notes

  • If states and times are not supplied, a simulation is run via simulateMarkovModel().
  • The function uses matplotlib's interactive mode (plt.ion()) to animate transitions.
  • The average sojourn time across all states is used to normalize animation speed.
  • Each visited state is highlighted in red for its corresponding (scaled) dwell time.
  • Drawing arguments like layout, font size, or color can be customized via kwargs.
def states(self, data=False):
1130    def states(self, data=False):
1131        """
1132        Return a view of the graph's states (nodes), optionally with attributes.
1133    
1134        This method is functionally identical to `self.nodes()` in NetworkX, but is
1135        renamed to `states()` to reflect the semantics of a state transition graph or
1136        Markov model, where nodes represent system states.
1137    
1138        Parameters
1139        ----------
1140        **data** : bool, optional
1141            If True, returns a view of `(node, attribute_dict)` pairs.
1142            If False (default), returns a view of node identifiers only.
1143    
1144        Returns
1145        -------
1146        *networkx.classes.reportviews.NodeView* or NodeDataView
1147            A view over the graph’s nodes or `(node, data)` pairs, depending on the `data` flag.
1148    
1149        Examples
1150        --------
1151        >>> G.states()
1152        ['A', 'B', 'C']
1153    
1154        >>> list(G.states(data=True))
1155        [('A', {'color': 'red', 'label': 'Start'}), ('B', {...}), ...]
1156    
1157        Notes
1158        -----
1159        - This is a convenience wrapper for semantic clarity in models where nodes represent states.
1160        - The view is dynamic: changes to the graph are reflected in the returned view.
1161        """
1162        return self.nodes(data=data)

Return a view of the graph's states (nodes), optionally with attributes.

This method is functionally identical to self.nodes() in NetworkX, but is renamed to states() to reflect the semantics of a state transition graph or Markov model, where nodes represent system states.

Parameters

data : bool, optional If True, returns a view of (node, attribute_dict) pairs. If False (default), returns a view of node identifiers only.

Returns

networkx.classes.reportviews.NodeView or NodeDataView A view over the graph’s nodes or (node, data) pairs, depending on the data flag.

Examples

>>> G.states()
['A', 'B', 'C']
>>> list(G.states(data=True))
[('A', {'color': 'red', 'label': 'Start'}), ('B', {...}), ...]

Notes

  • This is a convenience wrapper for semantic clarity in models where nodes represent states.
  • The view is dynamic: changes to the graph are reflected in the returned view.
def setStateColor(self, state, color):
1165    def setStateColor(self, state, color):
1166        """
1167        Set the color attribute of a specific state (node).
1168    
1169        Parameters
1170        ----------
1171        **state** : hashable
1172            The identifier of the state whose color is to be updated.
1173    
1174        **color** : str
1175            The new color to assign to the state (used for visualization).
1176    
1177        Returns
1178        -------
1179        None
1180        """
1181        self.nodes[state]["color"] = color

Set the color attribute of a specific state (node).

Parameters

state : hashable The identifier of the state whose color is to be updated.

color : str The new color to assign to the state (used for visualization).

Returns

None

def setStateLabel(self, state, label):
1184    def setStateLabel(self, state, label):
1185        """
1186        Set the label attribute of a specific state (node).
1187    
1188        Parameters
1189        ----------
1190        **state** : hashable
1191            The identifier of the state whose label is to be updated.
1192    
1193        **label** : str
1194            The new label to assign to the state (used for visualization or annotation).
1195    
1196        Returns
1197        -------
1198        None
1199        """
1200        self.nodes[state]["label"] = label

Set the label attribute of a specific state (node).

Parameters

state : hashable The identifier of the state whose label is to be updated.

label : str The new label to assign to the state (used for visualization or annotation).

Returns

None

def prob(self, state):
1202    def prob(self, state):
1203        """
1204        Return the steady-state probability of a specific state.
1205    
1206        Parameters
1207        ----------
1208        **state** : hashable
1209            The identifier of the state whose probability is to be retrieved.
1210    
1211        Returns
1212        -------
1213        **float**
1214            The steady-state probability of the specified state.
1215    
1216        Raises
1217        ------
1218        KeyError
1219            If the state is not found in the steady-state probability dictionary.
1220        """
1221        return self.state_probabilities[state]

Return the steady-state probability of a specific state.

Parameters

state : hashable The identifier of the state whose probability is to be retrieved.

Returns

float The steady-state probability of the specified state.

Raises

KeyError If the state is not found in the steady-state probability dictionary.

def probs(self):
1224    def probs(self):
1225        """
1226        Return the steady-state probabilities for all states.
1227    
1228        Returns
1229        -------
1230        **dict**
1231            A dictionary mapping each state to its steady-state probability.
1232    
1233        Notes
1234        -----
1235        - The probabilities are computed via `calcSteadyStateProbs()` and stored internally.
1236        - Use this method to access the full steady-state distribution.
1237        """
1238        return self.state_probabilities

Return the steady-state probabilities for all states.

Returns

dict A dictionary mapping each state to its steady-state probability.

Notes

  • The probabilities are computed via calcSteadyStateProbs() and stored internally.
  • Use this method to access the full steady-state distribution.
def getEmptySystemProbs(self, state=None):
1240    def getEmptySystemProbs(self, state=None):
1241        """
1242        Create an initial state distribution where all probability mass is in a single state.
1243    
1244        By default, this method returns a distribution where the entire probability mass
1245        is placed on the first node in the graph. Alternatively, a specific starting state
1246        can be provided.
1247    
1248        Parameters
1249        ----------
1250        **state** : hashable, optional
1251            The state to initialize with probability 1. If None, the first node
1252            in the graph (based on insertion order) is used.
1253    
1254        Returns
1255        -------
1256        **dict**
1257            A dictionary representing the initial distribution over states.
1258            The selected state has probability 1, and all others have 0.
1259    
1260        Examples
1261        --------
1262        >>> X0 = G.getEmptySystem()
1263        >>> sum(X0.values())  # 1.0
1264    
1265        Notes
1266        -----
1267        - This method is useful for setting up deterministic initial conditions
1268          in transient simulations of a CTMC.
1269        - The return format is compatible with methods like `transientProbs()`.
1270        """
1271        if state is None:
1272            first_node = next(iter(self.nodes))
1273        X0 = {s: 1 if s == first_node else 0 for s in self}
1274        return X0

Create an initial state distribution where all probability mass is in a single state.

By default, this method returns a distribution where the entire probability mass is placed on the first node in the graph. Alternatively, a specific starting state can be provided.

Parameters

state : hashable, optional The state to initialize with probability 1. If None, the first node in the graph (based on insertion order) is used.

Returns

dict A dictionary representing the initial distribution over states. The selected state has probability 1, and all others have 0.

Examples

>>> X0 = G.getEmptySystem()
>>> sum(X0.values())  # 1.0

Notes

  • This method is useful for setting up deterministic initial conditions in transient simulations of a CTMC.
  • The return format is compatible with methods like transientProbs().
def transientProbs(self, initialDistribution, t):
1278    def transientProbs(self, initialDistribution, t):
1279        """
1280        Compute the transient state probabilities at time `t` for the CTMC.
1281    
1282        This method solves the forward equation of a continuous-time Markov chain (CTMC):
1283        X(t) = X0 · expm(Q·t), where `Q` is the generator (rate) matrix and `X0` is the
1284        initial state distribution. It returns a dictionary mapping each state to its
1285        probability at time `t`.
1286    
1287        Parameters
1288        ----------
1289        - **initialDistribution** : dict 
1290            A dictionary mapping states (nodes) to their initial probabilities at time `t=0`.
1291            The keys must match the nodes in the graph, and the values should sum to 1.    
1292            
1293        - **t** : float
1294            The time at which the transient distribution is to be evaluated.
1295    
1296        Returns
1297        -------
1298        **dict**
1299            A dictionary mapping each state (node) to its probability at time `t`.
1300    
1301        Notes
1302        -----
1303        - The rate matrix `Q` is generated automatically if not yet computed.
1304        - The state order in the rate matrix corresponds to the internal `_n2i` mapping.
1305        - Internally, the dictionary `initialDistribution` is converted into a vector
1306          aligned with the state indexing.
1307        - The transient solution is computed using `scipy.linalg.expm` for matrix exponentiation.
1308    
1309        Raises
1310        ------
1311        ValueError
1312            If the input dictionary has missing states or probabilities that do not sum to 1.
1313            (This is not enforced but recommended for correctness.)
1314        """        
1315        if self._rate_matrix is None:  self.createRateMatrix()
1316        Q = self._rate_matrix
1317        X0 = np.zeros(len(initialDistribution))
1318        
1319        for i,n in enumerate(self.nodes):
1320            X0[ self._n2i[n] ] = initialDistribution[n]
1321                
1322        
1323        X = X0 @ expm(Q*t)    
1324        
1325        matrix = {}
1326        for n in self.nodes:
1327            matrix[n] = X[ self._n2i[n] ]
1328                
1329        return matrix

Compute the transient state probabilities at time t for the CTMC.

This method solves the forward equation of a continuous-time Markov chain (CTMC): X(t) = X0 · expm(Q·t), where Q is the generator (rate) matrix and X0 is the initial state distribution. It returns a dictionary mapping each state to its probability at time t.

Parameters

  • initialDistribution : dict A dictionary mapping states (nodes) to their initial probabilities at time t=0. The keys must match the nodes in the graph, and the values should sum to 1.

  • t : float The time at which the transient distribution is to be evaluated.

Returns

dict A dictionary mapping each state (node) to its probability at time t.

Notes

  • The rate matrix Q is generated automatically if not yet computed.
  • The state order in the rate matrix corresponds to the internal _n2i mapping.
  • Internally, the dictionary initialDistribution is converted into a vector aligned with the state indexing.
  • The transient solution is computed using scipy.linalg.expm for matrix exponentiation.

Raises

ValueError If the input dictionary has missing states or probabilities that do not sum to 1. (This is not enforced but recommended for correctness.)

def symSolveMarkovModel(self):
1331    def symSolveMarkovModel(self):
1332        """
1333        Symbolically solve the global balance equations of the CTMC.
1334    
1335        This method constructs and solves the system of symbolic equations for the
1336        steady-state probabilities of each state in the CTMC. It uses symbolic transition
1337        rates stored on the edges to form balance equations for each node, along with
1338        the normalization constraint that all probabilities sum to 1.
1339    
1340        Returns
1341        -------
1342        - **solution** : dict
1343            A symbolic solution dictionary mapping SymPy symbols (e.g., x_{A}) to
1344            expressions in terms of model parameters.
1345    
1346        - **num_solution** : dict
1347            A numerical dictionary mapping each state (node) to its steady-state probability,
1348            computed by substituting the numeric values from `self._subs` into the symbolic solution.
1349    
1350        Notes
1351        -----
1352        - For each node `s`, a symbolic variable `x_{label}` is created using the node's label attribute.
1353        - One balance equation is created per state: total inflow = total outflow.
1354        - An additional constraint `sum(x_i) = 1` ensures proper normalization.
1355        - The symbolic system is solved with `sympy.solve`, and the results are simplified.
1356        - Numeric values are computed by substituting known parameter values (`self._subs`) into the symbolic solution.
1357    
1358        Examples
1359        --------
1360        >>> sym_sol, num_sol = G.symSolveMarkovModel()
1361        >>> sym_sol  # symbolic expressions for each state
1362        >>> num_sol  # evaluated numerical values for each state
1363        """
1364        X={}
1365        to_solve = []
1366        for s in self.nodes():
1367            X[s] = sp.Symbol(f'x_{{{self.nodes[s]["label"]}}}')
1368            to_solve.append(X[s])
1369        #%
1370        eqs = []
1371        for s in self.nodes():
1372             succ = self.successors(s) # raus
1373             out_rate = 0
1374             for z in list(succ):
1375                 #out_rate += X[s]*self.edges[s, z]["sym_rate"]                 
1376                 out_rate += X[s]*self[s][z]["sym_rate"]
1377             
1378             in_rate = 0
1379             for r in list(self.predecessors(s)):
1380                 #in_rate += X[r]*self.edges[r,s]["sym_rate"] 
1381                 in_rate += X[r]*self[r][s]["sym_rate"] 
1382             eqs.append( sp.Eq(out_rate, in_rate) )
1383        #% 
1384        Xsum = 0    
1385        for s in X:
1386            Xsum += X[s]
1387        
1388        eqs.append( sp.Eq(Xsum, 1) )    
1389        #%
1390        solution = sp.solve(eqs, to_solve)
1391        simplified_solution = sp.simplify(solution)
1392        #print(simplified_solution)
1393        num_solution = simplified_solution.subs(self._subs)
1394        #num_dict =  {str(var): str(sol) for var, sol in simplified_solution.items()}
1395        num_dict = {s: num_solution[X[s]] for s in X}                
1396        
1397        return simplified_solution, num_dict 

Symbolically solve the global balance equations of the CTMC.

This method constructs and solves the system of symbolic equations for the steady-state probabilities of each state in the CTMC. It uses symbolic transition rates stored on the edges to form balance equations for each node, along with the normalization constraint that all probabilities sum to 1.

Returns

  • solution : dict A symbolic solution dictionary mapping SymPy symbols (e.g., x_{A}) to expressions in terms of model parameters.

  • num_solution : dict A numerical dictionary mapping each state (node) to its steady-state probability, computed by substituting the numeric values from self._subs into the symbolic solution.

Notes

  • For each node s, a symbolic variable x_{label} is created using the node's label attribute.
  • One balance equation is created per state: total inflow = total outflow.
  • An additional constraint sum(x_i) = 1 ensures proper normalization.
  • The symbolic system is solved with sympy.solve, and the results are simplified.
  • Numeric values are computed by substituting known parameter values (self._subs) into the symbolic solution.

Examples

>>> sym_sol, num_sol = G.symSolveMarkovModel()
>>> sym_sol  # symbolic expressions for each state
>>> num_sol  # evaluated numerical values for each state
def exportTransitionsToExcel( self, excel_path='output.xlsx', num_colors=9, open_excel=False, verbose=True):
1399    def exportTransitionsToExcel(self, excel_path = "output.xlsx", num_colors=9, open_excel=False, verbose=True):
1400        """
1401        Export transition data from the state transition graph to an Excel file.
1402        
1403        This method collects all transitions (edges) from the graph and writes them
1404        to an Excel spreadsheet, optionally applying background colors to distinguish
1405        outgoing transitions from different source states. This export is helpful for
1406        manually verifying the transition structure of the model.
1407
1408        The transition type is defined by the `tid` field of each edge, and mapped to
1409        a human-readable string via the `EXPORT_NAME_TRANSITIONS` dictionary in the
1410        constants module (e.g., `Constants.EXPORT_NAME_TRANSITIONS[tid]`).
1411
1412        See Also
1413        --------
1414        StateTransitionGraph.addTransition : Method to add edges with symbolic rates and `tid`.
1415        StateTransitionGraph.__init__ : Initializes the graph and loads the constants module.
1416
1417        Parameters
1418        ----------
1419        - **excel_path** : str, optional
1420            Path to the output Excel file. Default is "output.xlsx".
1421        
1422        - num_colors : int, optional
1423            Number of unique background colors to cycle through for different states.
1424            Default is 9.
1425        
1426        - open_excel : bool, optional
1427            If True, automatically opens the generated Excel file after export.
1428            Default is False.
1429        
1430        - verbose : bool, optional
1431            If True, prints summary information about the export process.
1432            Default is True.
1433        
1434        Returns
1435        -------
1436        None
1437            The function writes data to disk and optionally opens Excel,
1438            but does not return a value.
1439
1440        Notes
1441        -----
1442        - Each edge must contain a `tid` attribute for transition type and a `label` for the rate.
1443        - Background coloring groups all outgoing transitions from the same source state.
1444        - Requires the constants module to define `EXPORT_NAME_TRANSITIONS`.
1445        """
1446
1447        rows = []
1448        for node in self.nodes():
1449            edges = list(self.edges(node))
1450            if verbose: print(f"Edges connected to node '{node}':")
1451            
1452            for i, (u, v) in enumerate(edges):
1453                latex = self[u][v]["label"]
1454                rate = sp.pretty(self[u][v]["sym_rate"])
1455                tid = self[u][v]["tid"]
1456                
1457                if tid: # 
1458                    if hasattr(self._const, 'EXPORT_NAME_TRANSITIONS'):
1459                        tid_name = self._const.EXPORT_NAME_TRANSITIONS[tid]
1460                    else:
1461                        tid_name  = ''                               
1462                
1463                rows.append({"from": u, "type": tid_name, "to": v, "rate": rate, "latex": latex})
1464            
1465        df = pd.DataFrame(rows, columns=["from", "type", "to", "rate", "latex"])
1466        df["checked"] = ""
1467        df.to_excel(excel_path, index=False)
1468        
1469        # Load workbook with openpyxl
1470        wb = load_workbook(excel_path)
1471        ws = wb.active
1472        
1473        # Define colors for different "from" values
1474        cmap = map(plt.cm.Pastel1, range(num_colors))
1475        # Convert RGBA to hex
1476        color_list = [mcolors.to_hex(color, keep_alpha=False).upper().replace("#", "") for color in cmap]
1477        
1478        # Assign color to each unique "from" value using modulo
1479        unique_from = (df["from"].unique())  # Sorted for consistent order
1480        from_color_map = {
1481            value: color_list[i % len(color_list)] for i, value in enumerate(unique_from)
1482        }
1483        
1484        # Apply coloring row by row
1485        for i,row in enumerate(ws.iter_rows(min_row=2, max_row=ws.max_row)):  # Skip header
1486            from_value = row[0].value  # "from" is first column
1487            #fill_color = color_list[i % len(color_list)]
1488            if isinstance(from_value, int):
1489                fill_color = from_color_map[from_value]
1490            else:
1491                fill_color = from_color_map.get(eval(from_value))
1492            
1493            if fill_color:
1494                fill = PatternFill(start_color=fill_color, end_color=fill_color, fill_type="solid")
1495                for cell in row:
1496                    cell.fill = fill
1497        
1498        # Save styled Excel file
1499        wb.save(excel_path)
1500        if verbose: print(f"File {excel_path} created.")
1501        if open_excel: 
1502            print(f"Opening the file {excel_path}.")
1503            os.startfile(excel_path)

Export transition data from the state transition graph to an Excel file.

This method collects all transitions (edges) from the graph and writes them to an Excel spreadsheet, optionally applying background colors to distinguish outgoing transitions from different source states. This export is helpful for manually verifying the transition structure of the model.

The transition type is defined by the tid field of each edge, and mapped to a human-readable string via the EXPORT_NAME_TRANSITIONS dictionary in the constants module (e.g., Constants.EXPORT_NAME_TRANSITIONS[tid]).

See Also

StateTransitionGraph.addTransition : Method to add edges with symbolic rates and tid. StateTransitionGraph.__init__ : Initializes the graph and loads the constants module.

Parameters

  • excel_path : str, optional Path to the output Excel file. Default is "output.xlsx".

  • num_colors : int, optional Number of unique background colors to cycle through for different states. Default is 9.

  • open_excel : bool, optional If True, automatically opens the generated Excel file after export. Default is False.

  • verbose : bool, optional If True, prints summary information about the export process. Default is True.

Returns

None The function writes data to disk and optionally opens Excel, but does not return a value.

Notes

  • Each edge must contain a tid attribute for transition type and a label for the rate.
  • Background coloring groups all outgoing transitions from the same source state.
  • Requires the constants module to define EXPORT_NAME_TRANSITIONS.