MarkovModelModule

   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