amachine.am_hmm API documentation

HMM( states: Sequence[amachine.am_causal_state.CausalState], transitions: Sequence[amachine.am_transition.Transition], alphabet: Sequence[str], start_state: int = 0, composition_depth: int = 0, name: str = '', description: str = '') View Source

48    def __init__( 
49        self,
50        states      : Sequence[CausalState],
51        transitions : Sequence[Transition],
52        alphabet    : Sequence[str],
53        start_state : int = 0,
54        composition_depth : int = 0,
55        name        : str = "",
56        description : str = "" ) : 
57
58        self.unique_id = next(HMM._id_generator)
59
60        # MappingProxyType offers immutable view of dict
61        self.symbol_idx_map : MappingProxyType = MappingProxyType( {} )
62        self.state_idx_map  : MappingProxyType = MappingProxyType( {} )
63
64        self.alphabet    = tuple( alphabet    )   
65        self.states      = tuple( states      )
66        self.transitions = tuple( transitions )
67
68        self.update_symbol_idx_map()
69        self.update_state_idx_map()
70
71        self.name : str = name
72        self.description : str = description
73
74        self.composition_depth = composition_depth
75
76        self.start_state : int = start_state
77
78        # --- derived --------
79
80        self._complexity : dict[str, Any] = {}
81        self._pi_fractional = None
82        self._pi : np.ndarray | None = None
83        
84        self._T : np.ndarray | None = None
85        self._T_x  : list[np.ndarray] | None = None
86        self._msp : MSP | None = None
87
88        self._reverse_am : HMM | None = None
89        self._has_valid_rational_probabilities : bool = False
90        self._already_minimized : bool = False
91
92        # --- const ----------
93
94        self._EPS : float = 1e-12

unique_id

symbol_idx_map: mappingproxy

state_idx_map: mappingproxy

alphabet

states

transitions

name: str

description: str

composition_depth

start_state: int

def clear_cache(self): View Source

100    def clear_cache(self) :
101
102        """
103        Reset all derived properties so they will be recomputed when requested later.
104        """
105
106        self._complexity = {}
107        self._T = None
108        self._T_x = None
109        self._pi = None
110        self._pi_fractional = None
111        self._msp = None 
112        self._reverse_am = None
113
114        self._has_valid_rational_probabilities = False
115        self._already_minimized = False

Reset all derived properties so they will be recomputed when requested later.

def update_symbol_idx_map(self): View Source

117    def update_symbol_idx_map(self) :
118        new_idx_map = {}
119        for idx, symbol in enumerate( self.alphabet ) :
120            new_idx_map[ symbol ] = idx
121        self.symbol_idx_map = MappingProxyType( new_idx_map )

def update_state_idx_map(self): View Source

123    def update_state_idx_map(self) :
124        new_idx_map = {}
125        for idx, state in enumerate( self.states ) :
126            new_idx_map[ state.name ] = idx
127        self.state_idx_map = MappingProxyType( new_idx_map )

def set_states(self, states: Sequence[amachine.am_causal_state.CausalState]): View Source

129    def set_states( self, states : Sequence[CausalState] ) :
130        self.clear_cache()        
131        self.states = tuple( states )
132        self.update_state_idx_map()

def set_alphabet(self, alphabet: Sequence[str]): View Source

134    def set_alphabet( self, alphabet : Sequence[str] ) :
135
136        self.clear_cache()
137
138        old_alphabet = tuple( self.alphabet )
139
140        self.alphabet = tuple( sorted( set( alphabet ) ) )
141        self.update_symbol_idx_map()
142
143        new_transitions = []
144        for i, tr in enumerate( self.transitions ) :
145            symbol = old_alphabet[ tr.symbol_idx ]
146            new_transitions.append(  
147                tr.modified_deep_copy(
148                    symbol_idx=self.symbol_idx_map[ symbol ]
149                )
150            )
151
152        self.set_transitions( new_transitions )

def set_transitions(self, transitions: Sequence[amachine.am_transition.Transition]): View Source

154    def set_transitions( self, transitions : Sequence[Transition] ) :
155        self.clear_cache()
156        self.transitions = tuple( transitions )

def get_complexity_measure_if_exists(self, measure): View Source

158    def get_complexity_measure_if_exists(self, measure ) :
159        m = self._complexity.get( measure, None )
160        return m

def set_complexity_measure(self, measure, value): View Source

162    def set_complexity_measure(self, measure, value ) :
163        self._complexity[ measure ] = value

is_q_weighted View Source

184    @property
185    def is_q_weighted(self) :
186        return self._has_valid_rational_probabilities

def from_json_dict(self, config: dict[str, typing.Any]): View Source

192    def from_json_dict( self, config : dict[str, Any] ) :
193
194        self.name             = config.get( "name", "" )
195        self.description      = config.get( "description", "" )
196        self.start_state     = config.get( "start_state", 0 )
197        
198        json_states      = config.get( "states",      [] )
199        json_transitions = config.get( "transitions", [] )
200
201        states=[ 
202            CausalState.from_json_dict( state )
203            for state in json_states
204        ]
205
206        transitions=[ 
207            Transition.from_json_dict( tr )
208            for tr in json_transitions
209        ]
210
211        self.alphabet = config.get( "alphabet", () )
212
213        self.states      = tuple( states )
214        self.transitions = tuple( transitions )
215
216        self.update_symbol_idx_map()
217        self.update_state_idx_map()

def to_json_dict(self) -> dict[str, typing.Any]: View Source

219    def to_json_dict(self) -> dict[ str, Any ]:
220
221        """Create a dict representing the HMM configuration.
222
223        Returns:
224
225            dict[str,any]: Dictionary containing, name, description, states, transitions, alphabet.
226        """
227
228        return {
229            "name"            : self.name,
230            "description"     : self.description,
231            "start_state"     : self.start_state,
232            "states"          : [ state.to_json_dict()      for state      in self.states      ],
233            "transitions"     : [ transition.to_json_dict() for transition in self.transitions ],
234            "alphabet"        : list(self.alphabet)
235        }

Create a dict representing the HMM configuration.

Returns:

dict[str,any]: Dictionary containing, name, description, states, transitions, alphabet.

def save_config( self, output_dir: pathlib.Path | str, with_complexity: bool = False, with_non_trivial_complexity: bool = False): View Source

237    def save_config(
238        self, 
239        output_dir : Path | str, 
240        with_complexity : bool = False, 
241        with_non_trivial_complexity : bool = False ) :
242
243        output_dir = Path( output_dir )
244
245        config = self.to_json_dict()
246
247        if with_complexity :
248            
249            complexity = self.get_complexities( 
250                with_non_trivial=with_non_trivial_complexity
251            )
252
253            config[ "complexity" ] = complexity
254
255        config[ "structural_properties" ] = {
256            "unifilar"              : self.is_unifilar(),
257            "row_stochastic"        : self.is_row_stochastic(),
258            "strongly_connected"    : self.is_strongly_connected(),
259            "aperiodic"             : self.is_aperiodic() #,
260            # "minimal"               : self._is_minimal_as_dfa( topological_only=False )
261        }
262
263        with open( output_dir / "am_config.json", "w", encoding="utf-8" ) as f :
264            json.dump( config, f, ensure_ascii=False, indent=2, default=list )

def from_file(self, path: pathlib.Path): View Source

266    def from_file( self, path : Path ) :
267        with open( path / "am_config.json", "r" ) as f:
268            config = json.load(f)
269        self.from_json_dict( config )

def get_transition_list(self) -> list[list[tuple[int, float, int]]]: View Source

275    def get_transition_list(self) -> list[list[tuple[int, float, int]]] :
276        trs = [ [] for _ in range( len( self.states ) ) ]
277        for tr in self.transitions :
278            trs[ tr.origin_state_idx ].append( ( 
279                tr.symbol_idx, 
280                float( tr.prob ),
281                tr.target_state_idx ) )
282        return trs

def get_complexities(self, with_non_trivial=False): View Source

284    def get_complexities( 
285        self, 
286        with_non_trivial=False ) :
287
288        trivial = [
289            self.C_mu,
290            self.h_mu,
291            self.H_1,
292            self.rho_mu
293        ]
294
295        non_trivial = [
296            self.E, 
297            self.T_inf,
298            self.S,
299            self.chi
300        ]
301            
302        complexities = { m.__name__ : m() for m in trivial }
303
304        if with_non_trivial :
305
306            complexities |= { m.__name__ : float( m() ) for m in non_trivial }
307
308            _ = self.block_convergence()
309
310            scalar_complexity_keys = {
311                "E",    
312                "S",    
313                "T_inf"
314            }
315
316            complexities[ "block" ] = {}
317            for key in scalar_complexity_keys :
318                c = self.get_complexity_measure_if_exists( key )
319                if c is not None :
320                    complexities[ "block" ][ key ] = float( c )
321
322            block_complexity_keys = {
323                "E_L",   
324                "T_L",   
325                "S_L",   
326                "H_L",   
327                "h_mu_L",
328                "H_sync"
329            }
330
331            for key in block_complexity_keys :
332                bc = self.get_complexity_measure_if_exists( key )
333                if bc is not None :
334                    complexities[ "block" ][ key ] = [ float(x) for x in bc ]
335
336        return complexities

def get_metadata(self): View Source

340    def get_metadata(self) :
341
342        C = self.get_complexities( with_non_trivial=False )
343
344        return {
345            "name" : self.name,
346            'complexity'  : self._complexity,
347            "description" : self.description
348        }

def get_transition_matrix(self): View Source

350    def get_transition_matrix(self) :
351
352        if self._T  is not None :
353            return self._T
354
355        n_states = len( self.states )
356        T = np.zeros((n_states, n_states))
357
358        for tr in self.transitions :    
359            T[ tr.origin_state_idx, tr.target_state_idx  ] = tr.prob
360
361        self._T = T
362
363        return self._T

def get_T_X(self): View Source

367    def get_T_X(self) :
368
369        if self._T_x  is not None :
370            return self._T_x
371
372        n_states  = len( self.states )
373        n_symbols = len( self.alphabet )
374
375        T_x = [ np.zeros( ( n_states, n_states) ) for _ in range( n_symbols ) ]
376
377        for tr in self.transitions :
378            T_x[ tr.symbol_idx ][tr.origin_state_idx, tr.target_state_idx] = tr.prob
379
380        self._T_x = T_x
381        return self._T_x

def get_msp_qw( self, exact_state_cap: int = 1000, verbose: bool = True) -> amachine.am_msp.MSP: View Source

385    def get_msp_qw(
386        self,
387        exact_state_cap: int = 1000,
388        verbose: bool = True,
389    ) -> MSP :
390        if self._msp is not None:
391            return self._msp
392
393        try : 
394
395            print( "\nTrying to Compute Mixed State Presentation using Exact Fractions\n" )
396
397            self._msp = compute_msp_exact(
398                T_x=self.get_Tx_fractional(),
399                pi=self.get_fractional_stationary_distribution(),
400                n_states=len(self.states),
401                alphabet=self.alphabet,
402                exact_state_cap=exact_state_cap,
403                verbose=verbose
404            )
405
406            return self._msp 
407
408        except RuntimeError as e :
409            warnings.warn( f"Exact msp failed: {e} Falling back to msp approximation." )
410
411        return self.get_msp()

def get_msp( self, exact_state_cap: int = 1250000, verbose=True) -> amachine.am_msp.MSP: View Source

413    def get_msp(
414        self,
415        exact_state_cap: int = 1_250_000,
416        verbose = True,
417    ) -> MSP :
418
419        if self._msp is not None:
420            return self._msp
421     
422        T_x = self.get_T_X()
423        pi  = self.get_stationary_distribution()
424    
425        print( "\nComputing Mixed State Presentation..." )
426
427        self._msp = compute_msp( 
428            T_x=T_x,
429            pi=pi,
430            n_states=len(self.states),
431            alphabet=self.alphabet,
432            exact_state_cap=exact_state_cap,
433            verbose=verbose
434        )
435
436        return self._msp

def get_reverse_am(self): View Source

438    def get_reverse_am(self) :
439
440        was_q_weighted = self._has_valid_rational_probabilities
441
442        if self._reverse_am is not None:
443            return self._reverse_am
444
445        pi = self.get_stationary_distribution()
446        self._reverse_am = copy.deepcopy(self)
447        
448        new_transitions = []
449        for tr in self.transitions:
450            i = tr.target_state_idx
451            j = tr.origin_state_idx
452            
453            p_reversed = (pi[j] * tr.prob) / pi[i]
454            
455            new_transitions.append(
456                tr.modified_deep_copy(
457                    origin_state_idx=i,
458                    target_state_idx=j,
459                    prob=p_reversed,
460                    pq=None
461                )
462            )
463
464        self._reverse_am.set_transitions(new_transitions)
465
466        if self._reverse_am.is_epsilon_machine():
467            return self._reverse_am
468
469        rmsp = self._reverse_am.get_msp_qw( exact_state_cap=len(self.states)*4 )
470
471        self._reverse_am.set_states( rmsp.states )
472        self._reverse_am.set_transitions( rmsp.transitions )
473        self._reverse_am._msp = rmsp
474        self._reverse_am.start_state = 0
475
476        self._reverse_am.collapse_to_largest_strongly_connected_subgraph()
477        self._reverse_am.minimize()
478
479        if was_q_weighted :
480            self._reverse_am.to_q_weighted()
481
482        return self._reverse_am

def get_Tx_fractional(self) -> list[list[list[fractions.Fraction]]]: View Source

486    def get_Tx_fractional(self) -> list[ list[ list[ Fraction ] ] ] :
487
488        self.to_q_weighted()
489
490        n_states  = len( self.states )
491        n_symbols = len( self.alphabet )
492
493        T_x = []
494
495        for x in range( n_symbols ) :
496            T_x.append( [] )
497            for i in range( n_states ) :
498                T_x[ x ].append( [ 0 for _ in range( n_states ) ] )
499
500        for tr in self.transitions :
501            T_x[ tr.symbol_idx ][ tr.origin_state_idx ][ tr.target_state_idx ] = tr.pq
502
503        return T_x

def get_T_sympy(self): View Source

505    def get_T_sympy( self ) :
506
507        self.to_q_weighted()
508
509        n = len( self.states )
510        T = sympy.zeros( n, n )
511
512        for tr in self.transitions :
513            T[ tr.origin_state_idx, tr.target_state_idx ] = tr.pq
514
515        return T

def get_fractional_stationary_distribution(self): View Source

517    def get_fractional_stationary_distribution(self) :
518
519        T = self.get_T_sympy()
520
521        if self._pi_fractional is not None :
522            return self._pi_fractional
523
524        G = self.as_digraph()
525
526        if not nx.is_strongly_connected(G):
527            raise ValueError( "Single stationary distribution requires strongly connected HMM." )
528
529        self._pi_fractional = solve_for_pi_fractional( T )
530
531        return self._pi_fractional

def get_stationary_distribution(self): View Source

533    def get_stationary_distribution(self):
534
535        if self._pi is not None :
536            return self._pi
537
538        G = self.as_digraph()
539        
540        if not nx.is_strongly_connected(G):
541            raise ValueError( "Single stationary distribution requires strongly connected HMM." )
542
543        T = self.get_transition_matrix()
544        return solve_for_pi( T )

def C_mu(self): View Source

550    def C_mu( self ) :
551
552        """The *statistical complexity* (aka *forecasting complexity*) :
553
554        .. math::
555
556            C_{\\mu} = - \\sum_{\\sigma \\in \\mathcal{S}} \\Pr(\\sigma) \\log_2 \\Pr(\\sigma),
557
558        where :math:`\\mathcal{S}` is the set of states [^crutchfield_exact_2016], p.2.
559
560        .. note::
561
562            **Interpretations**
563
564            * The amount of historical information a process stores.
565            * The amount of structure in a process.
566
567        Returns:
568
569            float: :math:`C_{\\mu}`.
570
571        [^crutchfield_exact_2016]: Crutchfield et al., *Exact Complexity: The Spectral
572            Decomposition of Intrinsic Computation*, 2016.
573            <https://arxiv.org/abs/1309.3792>
574        """
575
576        m = self.get_complexity_measure_if_exists( "C_mu" )
577
578        if m is not None :
579            return m
580
581        pi = self.get_stationary_distribution()
582
583        h = 0
584        for i, pr in enumerate( pi ) :
585            
586            if pr < self._EPS :
587                continue
588
589            h += -pr * np.log2( pr )
590
591        self.set_complexity_measure( "C_mu", h )
592
593        return h

The statistical complexity (aka forecasting complexity) :

$$C_{\mu} = - \sum_{\sigma \in \mathcal{S}} \Pr(\sigma) \log_2 \Pr(\sigma),$$

where $ \mathcal{S} $ is the set of states ¹, p.2.

Interpretations

The amount of historical information a process stores.
The amount of structure in a process.

Returns:

float: $ C_{\mu} $.

Crutchfield et al., Exact Complexity: The Spectral Decomposition of Intrinsic Computation, 2016. https://arxiv.org/abs/1309.3792 ↩

def h_mu(self): View Source

597    def h_mu( self ) :
598
599        """The *entropy rate* :
600
601        .. math::
602
603            h_{\\mu}(\\boldsymbol{\\mathcal{S}}) = - \\sum_{\\sigma \\in \\mathcal{S}} \\Pr(\\sigma) \\sum_{x \\in \\mathcal{A}} \\Pr(x|\\sigma) \\log_2 \\Pr(x|\\sigma),
604
605        where :math:`\\mathcal{A}` is the alphabet and :math:`\\mathcal{S}` is the set of states [^crutchfield_exact_2016], p.2.
606
607        .. note::
608
609            **Interpretations**
610
611            * The lower bound on achievable loss in bits. 
612            * The irreducable randomness in the process.
613            * The intrinsic Randomness in the process.
614
615        Returns:
616            
617            float: :math:`h_{\\mu}`.
618
619        [^crutchfield_exact_2016]: Crutchfield et al., *Exact Complexity: The Spectral
620            Decomposition of Intrinsic Computation*, 2016.
621            <https://arxiv.org/abs/1309.3792>
622        """
623
624        m = self.get_complexity_measure_if_exists( "h_mu" )
625
626        if m is not None :
627            return m
628
629        T  = self.get_transition_matrix()
630        pi = self.get_stationary_distribution()
631
632        n_states = pi.size
633
634        h = 0
635        for i, pr in enumerate( pi ) :
636
637            if pr < self._EPS :
638                continue
639
640            row_entropy = 0
641            for j in range( len( pi ) ) :
642
643                if T[ i, j ]  < self._EPS :
644                    continue
645
646                row_entropy -= T[ i, j ] * np.log2( T[ i, j ] )
647
648            h += pr * row_entropy
649
650        self.set_complexity_measure( "h_mu", h )
651
652        return h

The entropy rate :

$$h_{\mu}(\boldsymbol{\mathcal{S}}) = - \sum_{\sigma \in \mathcal{S}} \Pr(\sigma) \sum_{x \in \mathcal{A}} \Pr(x|\sigma) \log_2 \Pr(x|\sigma),$$

where $ \mathcal{A} $ is the alphabet and $ \mathcal{S} $ is the set of states ¹, p.2.

Interpretations

The lower bound on achievable loss in bits.
The irreducable randomness in the process.
The intrinsic Randomness in the process.

Returns:

float: $ h_{\mu} $.

Crutchfield et al., Exact Complexity: The Spectral Decomposition of Intrinsic Computation, 2016. https://arxiv.org/abs/1309.3792 ↩

def H_1(self) -> float: View Source

656    def H_1(self) -> float :
657
658        """The *single symbol uncertainty*:
659
660        .. math::
661
662            H(1)=-\\sum_{x\\in\\mathcal{A}} \\Pr(x) \\log_2{\\Pr(x)},
663
664        where :math:`\\mathcal{A}` is the alphabet [^James_2018], p.2.
665
666        .. note::
667
668            **Interpretations**
669
670            * How uncertain you are on average about a single measurement with no context.
671
672        Returns:
673
674            float: :math:`H(1)`.
675
676        [^James_2018]: James et al., Anatomy of a Bit: Information in a Time Series Observation, 2018.
677            <https://arxiv.org/abs/1105.2988>
678        """
679
680        m = self.get_complexity_measure_if_exists("H_1")
681        if m is not None:
682            return m
683
684        pi  = self.get_stationary_distribution()
685        T_X = self.get_T_X()  # dict: symbol -> matrix
686
687        h = 0.0
688        for T_x in T_X:
689            # Pr(x) = sum_i pi[i] * sum_j T^(x)[i,j]
690            p_sym = 0.0
691            for i, pr in enumerate(pi):
692                if pr < self._EPS:
693                    continue
694                p_sym += pr * T_x[i, :].sum()
695
696            if p_sym < self._EPS:
697                continue
698            h -= p_sym * np.log2(p_sym)
699
700        self.set_complexity_measure("H_1", h)
701        return h

The single symbol uncertainty:

$$H(1)=-\sum_{x\in\mathcal{A}} \Pr(x) \log_2{\Pr(x)},$$

where $ \mathcal{A} $ is the alphabet ¹, p.2.

Interpretations

How uncertain you are on average about a single measurement with no context.

Returns:

float: $ H(1) $.

James et al., Anatomy of a Bit: Information in a Time Series Observation, 2018. https://arxiv.org/abs/1105.2988 ↩

def rho_mu(self) -> float: View Source

705    def rho_mu(self) -> float :
706        
707        """The *anticipated information* [^James_2018], p.3.:
708
709        .. math::
710
711            \\rho_{\\mu}= H(1) - h_{\\mu}
712
713        Returns:
714            
715            float: :math:`\\rho_{\\mu}`
716
717        [^James_2018]: James et al., Anatomy of a Bit: Information in a Time Series Observation, 2018.
718            <https://arxiv.org/abs/1105.2988>
719        """
720
721        m = self.get_complexity_measure_if_exists("rho_mu")
722        
723        if m is not None:
724            return m
725
726        rho = self.H_1() - self.h_mu()
727        
728        self.set_complexity_measure("rho_mu", rho)
729        
730        return rho

The anticipated information ¹, p.3.:

$$\rho_{\mu}= H(1) - h_{\mu}$$

Returns:

float: $ \rho_{\mu} $

James et al., Anatomy of a Bit: Information in a Time Series Observation, 2018. https://arxiv.org/abs/1105.2988 ↩

def block_convergence(self): View Source

734    def block_convergence( self )  :
735
736        """
737        Run [block entropy convergence](am_fast.html#block_entropy_convergence). Estimates [$\\mathbf{E}$](am_hmm.html#HMM.E), [$\\mathbf{S}$](am_hmm.html#HMM.S), [$\\mathbf{T}$](am_hmm.html#HMM.T_inf), and block measures[^crutchfield_exact_2016]:
738
739        - $\\mathbf{E}(L) = H(L) - L \\cdot h_{\\mu}$, 
740
741        - $\\mathbf{T}(L) = \\sum_{l=1}^{L} l \\left[ h_{\\mu}(l) - h_{\\mu} \\right]$, 
742
743        - $\\mathbf{S}(L) = \\sum_{l=0}^{L} \\mathcal{H}(l)$, 
744
745        - $H(L) = H[X_{0:L}]$, 
746
747        - $h_{\\mu}(L) = H(L) - H(L-1)$, and 
748
749        - $\\mathcal{H}(L) = -\\sum_{w \\in \\mathcal{A}^L} Pr(w) \\sum_{\\sigma \\in \\mathcal{S}} Pr(\\sigma|w) \\log_2 Pr(\\sigma|w)$.
750
751        You can plot these curves, and the block entropy curves using, [`amachine.HMM.draw_block_measure_curves`](am_hmm.html#HMM.draw_block_measure_curves), and [`amachine.HMM.draw_block_entropy_curve`](am_hmm.html#HMM.draw_block_entropy_curve). 
752
753        <img src="../resources/curves.png" alt="block measures plots" style="width: 100%; margin-left: 0%;">
754
755        Returns:
756        
757            ComplexityMeasures: An object containing the estimated measures with the following attributes:
758            
759            - E (float): The excess entropy.
760            - T_inf (float): The transient information ($\\mathbf{T}$).
761            - S (float): The synchronization information.
762            - E_L (numpy.ndarray): The block excess entropy ($\\mathbf{E}(L)$).
763            - T_L (numpy.ndarray): The block transient information ($\\mathbf{T}(L)$).
764            - S_L (numpy.ndarray): The block synchronization information ($\\mathbf{S}(L)$).
765            - H_L (numpy.ndarray): The block entropy ($H(L)$).
766            - h_mu_L (numpy.ndarray): The entropy rate estimates ($h_{\\mu}(L)$).
767            - H_sync (numpy.ndarray): The state-block synchronization ($\\mathcal{H}(L)$).
768            - converged (bool): True if the algorithm converged.
769
770        [^crutchfield_exact_2016]: Crutchfield et al., *Exact Complexity: The Spectral
771            Decomposition of Intrinsic Computation*, 2016.
772            <https://arxiv.org/abs/1309.3792>
773        """
774
775        trs = [ [] for _ in range( len( self.states ) ) ]
776        for tr in self.transitions :
777            trs[ tr.origin_state_idx ].append( ( 
778                tr.symbol_idx, 
779                float( tr.prob ),
780                tr.target_state_idx ) )
781
782        pi = self.get_stationary_distribution()
783
784        state_dist = [ float( pi[ i ] ) for i in range( len( self.states ) ) ]
785        branches = [(1.0, list(state_dist))]
786
787        print( "\nComputing Block Entropy\n" )
788
789        C = am_fast.block_entropy_convergence(
790            h_mu            = self.h_mu(),
791            n_states        = len( self.states ),
792            n_symbols       = len( self.alphabet ),
793            convergence_tol = 1e-8,
794            precision       = 15,
795            eps             = 1e-25,
796            branches        = branches,
797            trans           = trs,
798            max_branches    = 30_000_000
799        )
800
801        print( "Done\n" )
802
803        self.set_complexity_measure( f"E",       C.E )
804        self.set_complexity_measure( f"S",       C.S )
805        self.set_complexity_measure( f"T_inf",   C.T )
806        self.set_complexity_measure( f"E_L",     C.E_L.tolist() )
807        self.set_complexity_measure( f"T_L",     C.T_L.tolist() )
808        self.set_complexity_measure( f"S_L",     C.S_L.tolist() )
809        self.set_complexity_measure( f"H_L",     C.H_L.tolist() )
810        self.set_complexity_measure( f"h_mu_L",  C.h_mu_L.tolist() )
811        self.set_complexity_measure( f"H_sync",  C.H_sync.tolist() )
812
813        return C

Run block entropy convergence. Estimates $\mathbf{E}$, $\mathbf{S}$, $\mathbf{T}$, and block measures¹:

$\mathbf{E}(L) = H(L) - L \cdot h_{\mu}$,
$\mathbf{T}(L) = \sum_{l=1}^{L} l \left[ h_{\mu}(l) - h_{\mu} \right]$,
$\mathbf{S}(L) = \sum_{l=0}^{L} \mathcal{H}(l)$,
$H(L) = H[X_{0:L}]$,
$h_{\mu}(L) = H(L) - H(L-1)$, and
$\mathcal{H}(L) = -\sum_{w \in \mathcal{A}^L} Pr(w) \sum_{\sigma \in \mathcal{S}} Pr(\sigma|w) \log_2 Pr(\sigma|w)$.

You can plot these curves, and the block entropy curves using, amachine.HMM.draw_block_measure_curves, and amachine.HMM.draw_block_entropy_curve.

block measures plots

Returns:

ComplexityMeasures: An object containing the estimated measures with the following attributes:

E (float): The excess entropy.

T_inf (float): The transient information ($\mathbf{T}$).

S (float): The synchronization information.

E_L (numpy.ndarray): The block excess entropy ($\mathbf{E}(L)$).

T_L (numpy.ndarray): The block transient information ($\mathbf{T}(L)$).

S_L (numpy.ndarray): The block synchronization information ($\mathbf{S}(L)$).

H_L (numpy.ndarray): The block entropy ($H(L)$).

h_mu_L (numpy.ndarray): The entropy rate estimates ($h_{\mu}(L)$).

H_sync (numpy.ndarray): The state-block synchronization ($\mathcal{H}(L)$).

converged (bool): True if the algorithm converged.

Crutchfield et al., Exact Complexity: The Spectral Decomposition of Intrinsic Computation, 2016. https://arxiv.org/abs/1309.3792 ↩

def E(self) -> float: View Source

817    def E( self ) -> float :
818
819        """The *excess entropy* [^crutchfield_exact_2016], p.4:
820
821        .. math::
822
823            \\mathbf{E} \\equiv \\sum_{L=1}^{\\infty} I[X_{-\\infty:0}; X_{0:\\infty}]
824        
825        Computed via :meth:`get_msp` and :meth:`amachine.am_msp.MSP.get_E_S_T`, or :meth:`amachine.am_fast.block_entropy_convergence`
826
827        .. note::
828
829            **Interpretations**
830
831            * The information from the past that reduces uncertainty in the future [^crutchfield_exact_2016].
832            * How much information an observer must extract to synchronize to the process.
833            * Measures how long the process appears more complex than it asymptotically is.
834            * Vanishes for immediately synchronizable processes.
835
836        Returns:
837        
838            float: :math:`\\mathbf{E}`
839
840        [^crutchfield_exact_2016]: Crutchfield et al., *Exact Complexity: The Spectral
841            Decomposition of Intrinsic Computation*, 2016.
842            <https://arxiv.org/abs/1309.3792>
843        """
844
845        m = self.get_complexity_measure_if_exists( "E" )
846
847        if m is not None :
848            return m
849
850        try : 
851            msp = self.get_msp()
852            E, S, T = msp.get_E_S_T()
853            self.set_complexity_measure( "E", E )
854            self.set_complexity_measure( "S", S )
855            self.set_complexity_measure( "T_inf", T )
856            
857        except Exception as e :
858
859            print( f"MSP failed {e}" )
860
861            C = self.block_convergence()	
862            E = C.E
863            self.set_complexity_measure( "E", E )
864
865        return E

The excess entropy ¹, p.4:

$$\mathbf{E} \equiv \sum_{L=1}^{\infty} I[X_{-\infty:0}; X_{0:\infty}]$$

Computed via get_msp() and amachine.am_msp.MSP.get_E_S_T(), or amachine.am_fast.block_entropy_convergence()

Interpretations

The information from the past that reduces uncertainty in the future ¹.
How much information an observer must extract to synchronize to the process.
Measures how long the process appears more complex than it asymptotically is.
Vanishes for immediately synchronizable processes.

Returns:

float: $ \mathbf{E} $

Crutchfield et al., Exact Complexity: The Spectral Decomposition of Intrinsic Computation, 2016. https://arxiv.org/abs/1309.3792 ↩

def S(self) -> float: View Source

869    def S( self ) -> float :
870
871        """The *synchronization* information:
872
873        .. math::
874
875            \\mathbf{S} \\equiv \\sum_{L=1}^{\\infty} \\mathcal{H}(L),
876
877        where :math:`\\mathcal{H}(L)` is the average state uncertainty having seen all length-L words [^crutchfield_exact_2016], p.4.
878
879        .. note::
880
881            **Interpretations**
882
883            * The total amount of state information that an observer must extract to become synchronized [^crutchfield_exact_2016].
884
885        Computed via :meth:`get_msp` and :meth:`amachine.am_msp.MSP.get_E_S_T`, or :meth:`amachine.am_fast.block_entropy_convergence`
886
887        Returns:
888        
889            float: :math:`\\mathbf{S}`
890
891        [^crutchfield_exact_2016]: Crutchfield et al., *Exact Complexity: The Spectral
892            Decomposition of Intrinsic Computation*, 2016.
893            <https://arxiv.org/abs/1309.3792>
894        """
895
896        m = self.get_complexity_measure_if_exists( "S" )
897
898        if m is not None :
899            return m
900
901        try : 
902            msp = self.get_msp()
903            E, S, T = msp.get_E_S_T()
904            self.set_complexity_measure( "E", E )
905            self.set_complexity_measure( "S", S )
906            self.set_complexity_measure( "T_inf", T )
907
908        except Exception as e :
909            print( f"{e} \nFalling back to iterative estimation.")
910            C = self.block_convergence()	
911            S = C.S
912            self.set_complexity_measure( "S", S )
913
914        return S

The synchronization information:

$$\mathbf{S} \equiv \sum_{L=1}^{\infty} \mathcal{H}(L),$$

where $ \mathcal{H}(L) $ is the average state uncertainty having seen all length-L words ¹, p.4.

Interpretations

The total amount of state information that an observer must extract to become synchronized ¹.

Computed via get_msp() and amachine.am_msp.MSP.get_E_S_T(), or amachine.am_fast.block_entropy_convergence()

Returns:

float: $ \mathbf{S} $

Crutchfield et al., Exact Complexity: The Spectral Decomposition of Intrinsic Computation, 2016. https://arxiv.org/abs/1309.3792 ↩

def T_inf(self) -> float: View Source

918    def T_inf( self ) -> float :
919
920        """The *transient information*[^crutchfield_exact_2016], p.4:
921
922        .. math::
923
924            \\mathbf{T} \\equiv \\sum_{L=1}^{\\infty} L \\left[ h_{\\mu}(L) - h_{\\mu} \\right]
925
926        Computed via :meth:`get_msp` and :meth:`amachine.am_msp.MSP.get_E_S_T`, or :meth:`amachine.am_fast.block_entropy_convergence`
927
928        .. note::
929
930            **Interpretations**
931
932            * The amount of information one must extract from observations so that the block entropy converges to its linear asymptote[^crutchfield_exact_2016].
933
934        Returns:
935        
936            float: :math:`\\mathbf{T}`
937
938        [^crutchfield_exact_2016]: Crutchfield et al., *Exact Complexity: The Spectral
939            Decomposition of Intrinsic Computation*, 2016.
940            <https://arxiv.org/abs/1309.3792>
941        """
942
943        m = self.get_complexity_measure_if_exists( "T_inf" )
944
945        if m is not None :
946            return m
947
948        try : 
949            msp = self.get_msp()
950            E, S, T = msp.get_E_S_T()
951            self.set_complexity_measure( "E", E )
952            self.set_complexity_measure( "S", S )
953            self.set_complexity_measure( "T_inf", T )
954            T_inf = T
955
956        except Exception as e :
957            print( f"{e} \nFalling back to iterative estimation.")
958            C = self.block_convergence()	
959            T_inf = C.T
960
961        return T_inf

The transient information¹, p.4:

$$\mathbf{T} \equiv \sum_{L=1}^{\infty} L \left[ h_{\mu}(L) - h_{\mu} \right]$$

Computed via get_msp() and amachine.am_msp.MSP.get_E_S_T(), or amachine.am_fast.block_entropy_convergence()

Interpretations

The amount of information one must extract from observations so that the block entropy converges to its linear asymptote¹.

Returns:

float: $ \mathbf{T} $

Crutchfield et al., Exact Complexity: The Spectral Decomposition of Intrinsic Computation, 2016. https://arxiv.org/abs/1309.3792 ↩

def chi(self) -> float: View Source

 965    def chi( self ) -> float :
 966
 967        """The foward crypticity[^crutchfield_crypticity_2009][^Mahoney_crypticity_2021], p.2:
 968
 969        .. math::
 970
 971            \\chi = C_{\\mu} - \\mathbf{E}
 972
 973        :math:`C_{\\mu}` is trivially computed from the stationary distribution in :meth:`C_mu` and :math:`\\mathbf{E}` in :meth:`E`.
 974
 975        .. note::
 976
 977            **Interpretations**
 978
 979            * Difference between internal stored information and apparent information to an observer.
 980            * How muching information is hiding in the system.
 981
 982        Returns:
 983        
 984            float: :math:`\\chi`
 985
 986        [^crutchfield_crypticity_2009]: Crutchfield et al., Time’s barbed arrow: Irreversibility, crypticity, and stored information, 2009.
 987            <https://arxiv.org/abs/0902.1209>
 988
 989        [^Mahoney_crypticity_2021]: Mahoney et al., Information Accessibility and Cryptic Processes, 2021.
 990            <https://arxiv.org/abs/0905.4787>
 991        """
 992
 993        m = self.get_complexity_measure_if_exists( "chi" )
 994
 995        if m is not None :
 996            return m
 997
 998        chi = self.C_mu() - self.E()
 999
1000        if chi < 0 :
1001            
1002            # if chi is 0, accumulated floating point error can result in small negative values
1003            if chi < -1e-5:
1004                warnings.warn(f"Crypticity is negative ({chi:.6e}).")
1005            
1006            chi = np.clip( chi, 0 )
1007
1008        self.set_complexity_measure( "chi", chi )
1009
1010        return chi

The foward crypticity¹², p.2:

$$\chi = C_{\mu} - \mathbf{E}$$

$ C_{\mu} $ is trivially computed from the stationary distribution in C_mu() and $ \mathbf{E} $ in E().

Interpretations

Difference between internal stored information and apparent information to an observer.
How muching information is hiding in the system.

Returns:

float: $ \chi $

Crutchfield et al., Time’s barbed arrow: Irreversibility, crypticity, and stored information, 2009. https://arxiv.org/abs/0902.1209 ↩
Mahoney et al., Information Accessibility and Cryptic Processes, 2021. https://arxiv.org/abs/0905.4787 ↩

def is_row_stochastic(self): View Source

1016    def is_row_stochastic(self) :
1017
1018        """
1019        Check that all states have outgoing transition probabilities that sum to 1.
1020        """
1021
1022        sums = np.zeros( len( self.states ) )
1023        for tr in self.transitions :
1024            sums[ tr.origin_state_idx ] += tr.prob
1025
1026        return np.allclose( sums, 1.0 )

Check that all states have outgoing transition probabilities that sum to 1.

def is_unifilar(self): View Source

1030    def is_unifilar(self):
1031        """
1032        Check that no state emits the same symbol on transitions to different states.
1033        """
1034
1035        seen_state_symbols = set()
1036        for tr in self.transitions:
1037
1038            sym_key = (tr.origin_state_idx, tr.symbol_idx)
1039
1040            # If we've seen this origin+symbol or origin+target before
1041            if sym_key in seen_state_symbols :
1042                return False
1043
1044            seen_state_symbols.add(sym_key)
1045
1046        return True

Check that no state emits the same symbol on transitions to different states.

def is_strongly_connected(self): View Source

1050    def is_strongly_connected(self) :
1051
1052        """
1053        Check if every state is reachable from every other state. Relies on [nx.is_strongly_connected](https://networkx.org/documentation/latest/reference/algorithms/generated/networkx.algorithms.components.is_strongly_connected.html).
1054        """
1055
1056        return nx.is_strongly_connected( self.as_digraph() )

Check if every state is reachable from every other state. Relies on nx.is_strongly_connected.

def is_aperiodic(self): View Source

1060    def is_aperiodic(self) :
1061
1062        """
1063        Checks if machine is periodic. Relies on [nx.is_aperiodic](https://networkx.org/documentation/latest/reference/algorithms/generated/networkx.algorithms.dag.is_aperiodic.html), "A strongly connected directed graph is aperiodic if there is no integer k > 1 that divides the length of every cycle in the graph."
1064        """
1065
1066        return nx.is_aperiodic( self.as_digraph() )

Checks if machine is periodic. Relies on nx.is_aperiodic, "A strongly connected directed graph is aperiodic if there is no integer k > 1 that divides the length of every cycle in the graph."

def is_topological_epsilon_machine(self, verbose=True): View Source

1089    def is_topological_epsilon_machine( self, verbose=True ) :
1090
1091        """
1092        Checks if the HMM is a topological $\\epsilon$-machine [^1].
1093
1094        [^1]: Johnson et al, Enumerating Finitary Processes, 2024.
1095            <https://arxiv.org/abs/1011.0036>
1096        """
1097
1098        if not ( self.is_unifilar() and self.is_strongly_connected() ) :
1099            if verbose : 
1100                print( f"Either non unifilar or not strongly connected" )
1101            return False
1102        else :
1103            return self._is_minimal_as_dfa( topological_only=True, verbose=verbose )

Checks if the HMM is a topological $\epsilon$-machine ¹.

Johnson et al, Enumerating Finitary Processes, 2024. https://arxiv.org/abs/1011.0036 ↩

def is_epsilon_machine(self, verbose=True): View Source

1105    def is_epsilon_machine( self, verbose=True ) :
1106
1107        if not ( self.is_unifilar() and self.is_strongly_connected() ) :
1108            if verbose : 
1109                print( f"Either non unifilar or not strongly connected" )
1110            return False
1111        else :
1112            return self._is_minimal_as_dfa( topological_only=False, verbose=verbose )

def minimize(self, retain_names: bool = True, verbose=False): View Source

1118    def minimize(self, retain_names: bool = True, verbose=False):
1119
1120        """
1121        Minimizes the HMM, resulting in an :math:`\\epsilon-`machine if the HMM
1122        is unifilar and strongly connected. Converts the HMM to a DFA with symbols
1123        labeled jointly with symbols and probabilities, and uses Myhill-Nerode 
1124        equivalence for minimization. Relies on `automata_lib` and uses
1125         `automata.fa.dfa.DFA.minify` with `allow_partial=True`, and all states
1126         final.
1127
1128        Args:
1129            retain_names (bool): If `True`, the merged states will be named by their union, e.g. `{s_0, s_1}`, and other states will retain their origion names. Otherwise, they will be relabled `{ '0', '1', ..., 'n-1' }`.
1130
1131        Returns:
1132        
1133            automata.fa.dfa.DFA : the resulting DFA.
1134        """
1135
1136        if self._already_minimized :
1137            return
1138
1139        start = time.perf_counter()
1140
1141        if not self.is_unifilar():
1142            raise ValueError(
1143                "DFA minimization is not valid for non-unifilar HMMs"
1144            )
1145
1146        was_strongly_connected = self.is_strongly_connected()
1147        was_row_stochastic     = self.is_row_stochastic()
1148        was_unifilar           = self.is_unifilar()
1149
1150        n_states_before = len(self.states)
1151
1152        dfa = self.as_dfa(with_probs=True)
1153
1154        min_dfa = self.as_dfa(with_probs=True).minify(retain_names=True)
1155        #min_dfa = am_fast.minify_cpp( dfa, retain_names=True )
1156
1157        # Build lookup from original state index -> CausalState object
1158        orig_state   = {i: s for i, s in enumerate(self.states)}
1159        eq_list      = list(min_dfa.states)
1160
1161        start_eq = min_dfa.initial_state
1162        
1163        # Separate the start state, then sort the rest by the 
1164        # smallest original state index inside each equivalence class.
1165        other_eqs = [eq for eq in eq_list if eq != start_eq]
1166        other_eqs.sort(key=lambda eq: min(eq))
1167
1168        # Recombine so start eq comes first, followed by the sorted remaining classes
1169        eq_list = [start_eq] + other_eqs
1170        # ----------------------------------------------------------
1171
1172        # Recompute eq_to_idx with the new ordering
1173        eq_to_idx = {eq: i for i, eq in enumerate(eq_list)}
1174
1175        # new_start is now guaranteed to be 0
1176        new_start = 0
1177
1178        # Map each original state index -> its equivalence class
1179        # Guard: minify() silently drops unreachable states
1180        orig_to_eq = {s: eq for eq in min_dfa.states for s in eq}
1181
1182        # Build lookup from original state index -> its transitions
1183        orig_trs = defaultdict(list)
1184        for t in self.transitions:
1185            orig_trs[t.origin_state_idx].append(t)
1186
1187        new_trs = []
1188        for eq in min_dfa.states:
1189            rep        = next(iter(eq))
1190            origin_idx = eq_to_idx[eq]
1191            for t in orig_trs[rep]:
1192
1193                target_eq  = orig_to_eq[t.target_state_idx]
1194                target_idx = eq_to_idx[target_eq]
1195                
1196                new_trs.append(
1197                    t.modified_deep_copy(
1198                        origin_state_idx = origin_idx,
1199                        target_state_idx = target_idx
1200                    )
1201                )
1202
1203        members_list = [[orig_state[i] for i in sorted(eq)] for eq in eq_list]  # sorted for determinism
1204
1205        # Compute new names
1206        if retain_names:
1207            new_names = [
1208                "{" + ",".join(str(m.name) for m in members) + "}" if len(members) > 1
1209                else members[0].name
1210                for members in members_list
1211            ]
1212        else:
1213            new_names = [str(j) for j in range(len(eq_list))]
1214
1215        old_name_to_new_name = {
1216            m.name: new_names[j]
1217            for j, members in enumerate(members_list)
1218            for m in members
1219        }
1220
1221        # Build the new states, preserving classes and isomorphs regardless of naming
1222        new_states = []
1223        for j, (eq, members, name) in enumerate(zip(eq_list, members_list, new_names)):
1224            
1225            classes : defaultdict[int,set[str]]= {}
1226            for m in members:
1227                classes |= m.classes
1228            
1229            isomorphs = {
1230                old_name_to_new_name.get(iso, iso)
1231                for m in members
1232                for iso in m.isomorphs
1233                if old_name_to_new_name.get(iso, iso) != name
1234            }
1235            
1236            pseudo_isomorphs = {
1237                old_name_to_new_name.get(iso, iso)
1238                for m in members
1239                for iso in m.pseudo_isomorphs
1240                if old_name_to_new_name.get(iso, iso) != name
1241            }
1242
1243            new_states.append(CausalState(
1244                name      = name,
1245                classes   = classes,
1246                isomorphs = isomorphs,
1247                pseudo_isomorphs = pseudo_isomorphs
1248            ))
1249
1250        self.set_states(new_states)
1251        self.set_transitions(new_trs)
1252        self.start_state = new_start
1253
1254        if n_states_before == len(new_states) and verbose :
1255            print( f"{n_states_before} state HMM was already minimal.\n" )
1256        elif verbose :
1257            print( f"Minimized from {n_states_before} to {len(new_states)}\n" )
1258
1259        if not ( was_strongly_connected ==  self.is_strongly_connected() ) :
1260            raise RuntimeError(
1261                f"Minimization broke strongly connected"
1262            )
1263
1264        if not ( was_row_stochastic ==  self.is_row_stochastic() ) :
1265            raise RuntimeError(
1266                f"Minimization broke row stochasticity"
1267            )
1268
1269        if not( was_unifilar == self.is_unifilar() ) :
1270            raise RuntimeError(
1271                f"Minimization broke unifilarity"
1272            )
1273
1274        self._already_minimized = True

Minimizes the HMM, resulting in an $ \epsilon- $machine if the HMM is unifilar and strongly connected. Converts the HMM to a DFA with symbols labeled jointly with symbols and probabilities, and uses Myhill-Nerode equivalence for minimization. Relies on automata_lib and uses automata.fa.dfa.DFA.minify with allow_partial=True, and all states final.

Arguments:

retain_names (bool): If True, the merged states will be named by their union, e.g. {s_0, s_1}, and other states will retain their origion names. Otherwise, they will be relabled { '0', '1', ..., 'n-1' }.

Returns:

automata.fa.dfa.DFA : the resulting DFA.

def collapse_to_largest_strongly_connected_subgraph(self, rename_states=True): View Source

1281    def collapse_to_largest_strongly_connected_subgraph( self, rename_states=True ) :
1282
1283        was_q_weighted = self._has_valid_rational_probabilities
1284
1285        # get equivalent networkx graph
1286        G = self.as_digraph()
1287
1288        # if already strongly connected, nothing to do
1289        if not nx.is_strongly_connected( G ) :
1290
1291            start = time.perf_counter()
1292            subgraph_nodes = list( nx.strongly_connected_components( G ) )
1293
1294            # decompose into strongly connected components and sort by length
1295            # subgraph_nodes = list(nx.strongly_connected_components( G ))
1296            subgraph_nodes.sort(key=len)
1297            component_state_set = subgraph_nodes[-1]
1298
1299            # Take the largest strongly connected component (as list of state names)
1300            component_states = sorted( list( component_state_set ) )
1301
1302            # make temporary copies of the old transitions and states
1303            old_transitions = [
1304                tr.deepcopy()
1305                for tr in self.transitions
1306            ]
1307
1308            old_states = [
1309                s.deepcopy()
1310                for s in self.states
1311            ]
1312
1313            self.set_states(
1314                states=[ 
1315                    state
1316                    for i, state in enumerate( old_states ) if i in component_state_set
1317                ]
1318            )
1319
1320            # we will build new transition list based on those belonging to the component
1321            self.set_transitions( transitions= [] )
1322
1323            # for tracking which new transitions leave each state
1324            transitions_from_state = { state : set() for state in component_states }
1325            new_transitions = []
1326
1327            for tr in old_transitions :
1328
1329                origin_state_name = old_states[ tr.origin_state_idx ].name
1330                target_state_name = old_states[ tr.target_state_idx ].name
1331
1332                # skip transitions that connect separate strongly connected components
1333                if not ( tr.origin_state_idx in component_state_set and tr.target_state_idx in component_state_set ) :
1334                    continue
1335
1336                # track transitions (by index in new transitions list) that leave this state
1337                transitions_from_state[ tr.origin_state_idx ].add( len( new_transitions ) )
1338
1339                my_origin_state_idx = self.state_idx_map[ origin_state_name ]
1340                my_target_state_idx = self.state_idx_map[ target_state_name ]
1341
1342                new_transitions.append( 
1343                    tr.modified_deep_copy(  
1344                        origin_state_idx=my_origin_state_idx,
1345                        target_state_idx=my_target_state_idx
1346                    )
1347                )
1348
1349            self.set_transitions( transitions=new_transitions )
1350
1351            # if we removed an outgoing transition from a state, we need to distribute its probability 
1352            # among the remaining outgoing transitions from the state
1353            
1354            transition_list = list( self.transitions )
1355            
1356            for state in component_states :
1357                
1358                # get the set of transitions leaving this state
1359                state_trs = transitions_from_state[ state ]
1360
1361                # sum the probabilities of the outgoing transitions from the state
1362                p_sum = np.sum( [ self.transitions[ i ].prob for i in state_trs ] )
1363
1364                # how much probability is missing
1365                diff = 1.0 - p_sum
1366
1367                # if significant difference
1368                if abs( diff ) > self._EPS : 
1369
1370                    # calculate how much of the difference each transition gets
1371                    adjustment = diff / len( state_trs )
1372                    
1373                    new_transitions = []
1374
1375                    # update the transitions
1376                    for i in state_trs :
1377
1378                        # adjusted probability
1379                        transition_list[ i ] = self.transitions[ i ].modified_deep_copy(
1380                            prob=self.transitions[ i ].prob + adjustment,
1381                            pq=None
1382                        )
1383
1384            self.set_transitions( transition_list )
1385
1386            if rename_states :
1387                self.set_states( [
1388                    s.modified_deep_copy( name=f"{i}" )
1389                    for s in self.states
1390                ] )
1391
1392        if was_q_weighted :
1393            self.to_q_weighted()

def to_q_weighted(self, denominator_limit=1000): View Source

1395    def to_q_weighted( self, denominator_limit=1000 ) :
1396
1397        """
1398        Approximates the existing transition probabilities with exact fractions, stores 
1399        the fractional probabilities as Fraction in Transition.pq, and sets the floating
1400        point probabilty to `float(pq)`. If `denominator_limit` is too small for a sane 
1401        conversion, the function recurses with `denominator_limit=denominator_limit*10`.
1402
1403        Args:
1404            denominator_limit (int): The initial input to :meth:`Fraction.limit_denominator` in 
1405            the conversion.
1406        """
1407
1408        if self._has_valid_rational_probabilities :
1409            return
1410
1411        if not self.is_row_stochastic() :
1412            raise ValueError( "Cannot convert to q-weighted because not row stochastic" )
1413
1414        t_from = [[] for _ in range(len(self.states))]
1415
1416        for i, tr in enumerate( self.transitions ) :
1417            t_from[ tr.origin_state_idx ].append( i )
1418
1419        new_transitions = []
1420        for t_list in t_from :
1421            
1422            if not t_list : 
1423                continue
1424
1425            p_q_sum = Fraction(0,1)
1426            p_qs = []
1427
1428            for t_idx in t_list :
1429            
1430                p_q = Fraction( self.transitions[ t_idx ].prob ).limit_denominator( denominator_limit )
1431                p_q_sum += p_q
1432                p_qs.append( p_q )
1433
1434            if p_q_sum != Fraction(1,1) :
1435                
1436                max_pq_i = np.argmax( p_qs )
1437                max_oq = p_qs[ max_pq_i ]
1438
1439                diff = p_q_sum - Fraction(1,1)
1440
1441                # If recurse with higher resolution
1442                if diff > max_oq :
1443                    return self.to_q_weighted( denominator_limit*10 )
1444                else :
1445                    p_qs[ max_pq_i ] -= diff
1446
1447            for i, t_idx in enumerate( t_list ) :	
1448                new_transitions.append( 
1449                    self.transitions[  t_idx ].modified_deep_copy(
1450                        prob=float(p_qs[ i ]),
1451                        pq=p_qs[ i ]
1452                    )
1453                )
1454
1455        self.set_transitions( new_transitions )
1456        self._has_valid_rational_probabilities = True

Approximates the existing transition probabilities with exact fractions, stores the fractional probabilities as Fraction in Transition.pq, and sets the floating point probabilty to float(pq). If denominator_limit is too small for a sane conversion, the function recurses with denominator_limit=denominator_limit*10.

Arguments:

denominator_limit (int): The initial input to Fraction.limit_denominator() in
the conversion.

def isomorphic_shift( self, input_symbol_indices: numpy.ndarray, input_state_indices: numpy.ndarray, shift: int = 1) -> dict[str, numpy.ndarray]: View Source

1462    def isomorphic_shift(
1463        self,
1464        input_symbol_indices: np.ndarray,
1465        input_state_indices:  np.ndarray,
1466        shift : int = 1
1467    ) -> dict[str, np.ndarray]:
1468
1469        """
1470        Generates a new sequence of symbols that are permuted with the symbols emitted by
1471        isomorphic states, if they exists.
1472
1473        :math:`\\sigma_o = \\mathcal{S}\\left[\\texttt{input\\_state\\_indices}[i]\\right]`<br>
1474        :math:`\\sigma_t = \\mathcal{S}\\left[\\texttt{input\\_state\\_indices}[i+1]\\right]`
1475 
1476        :math:`\\mathcal{I}(\\sigma_o) = \\\\{\\sigma^0_o,\\, \\sigma^1_o,\\, \\dots,\\, \\sigma^{n-1}_o \\\\}`<br>
1477        :math:`\\mathcal{I}(\\sigma_t) = \\\\{\\sigma^0_t,\\, \\sigma^1_t,\\, \\dots,\\, \\sigma^{n-1}_t \\\\}`
1478 
1479        :math:`k = \\bigl(\\mathcal{I}(\\sigma_o).\\texttt{index}(\\sigma_o) + \\texttt{shift}\\bigr) \\bmod n`
1480 
1481        :math:`\\texttt{output\\_symbol\\_indices}[i]   := T(\\sigma_o^k,\\, \\sigma_t^k).\\text{symbol\\_index}`<br>
1482        :math:`\\texttt{output\\_state\\_indices}[i]    := \\mathcal{S}.\\texttt{index}(\\sigma_o^k)`<br>
1483        :math:`\\texttt{output\\_state\\_indices}[i+1]  := \\mathcal{S}.\\texttt{index}(\\sigma_t^k)`
1484
1485        Where :math:`\\mathcal{I}(\\sigma)`` is the ordered set of states isomorphic to :math:`\\sigma` including :math:`\\sigma` itself.
1486
1487        Args:
1488            input_symbol_indices (np.ndarray): The sequence of generated symbols.
1489            input_state_indices (np.ndarray): The sequence of states that generated symbols with the final state at the end.
1490            shift : int: How much to shift the symbols across the isomorphic states.
1491        """
1492        
1493        if not any( state.isomorphs for state in self.states ):
1494            raise ValueError("HMM has no states with isomorphs")
1495
1496        inputs = np.asarray(input_symbol_indices)
1497        states = np.asarray(input_state_indices)
1498
1499        n_states = len(self.states)
1500
1501        tr_sym_table    = np.full((n_states, n_states), -1, dtype=np.int32)
1502        is_pseudo_table = np.zeros((n_states, n_states), dtype=bool)
1503        tr_cross_table  = np.zeros((n_states, n_states), dtype=bool)
1504
1505        for tr in self.transitions:
1506            tr_sym_table[ tr.origin_state_idx, tr.target_state_idx ] = tr.symbol_idx
1507            tr_cross_table[tr.origin_state_idx, tr.target_state_idx] = tr.composition_depth > 0
1508
1509        for i, state in enumerate(self.states):
1510            for p_iso in state.pseudo_isomorphs:
1511                j = self.state_idx_map[p_iso]
1512                is_pseudo_table[i, j] = True
1513
1514        # Build isomorph remapping: identity by default, overridden where isomorphs exist
1515        iso_table = np.arange( n_states, dtype=np.int32 )
1516
1517        for i, state in enumerate(self.states):
1518            effective_isos = state.isomorphs | state.pseudo_isomorphs
1519            if len(effective_isos) > 0:
1520                isormorphs_with_identity = sorted([i] + [self.state_idx_map[iso] for iso in effective_isos])
1521                pos = isormorphs_with_identity.index(i)
1522                iso_table[i] = isormorphs_with_identity[(pos + shift) % len(isormorphs_with_identity)]
1523
1524        origins = states[:-1]
1525        targets = states[1:]
1526
1527        out_origins = iso_table[origins]
1528        out_targets = iso_table[targets]
1529
1530        rotated_symbols = tr_sym_table[ out_origins, out_targets ]
1531
1532        ###########################################################################
1533        # handle attempts to rotate symbols emitted by pseudo-isomorphic states
1534        
1535        # are the isomorphs used for the origin and target shift pseudo-isomorphism?
1536        origin_is_pseudo = is_pseudo_table[origins, out_origins]
1537        target_is_pseudo = is_pseudo_table[targets, out_targets]
1538
1539        # does the edge cross a component
1540        is_cross = tr_cross_table[ origins, targets ]
1541
1542        # is_cross is sufficient condition for single-level compositions, 
1543        # but condition the combined is required for multi-level composition
1544        invalid_shift = is_cross & origin_is_pseudo & target_is_pseudo
1545
1546        final_symbols = np.where( invalid_shift, inputs, rotated_symbols)
1547
1548        if np.any(final_symbols == -1):
1549            raise RuntimeError(
1550                "Invalid isomorphic shift: topology invariant violated."
1551            )
1552
1553        ############################################################################
1554
1555        invalid_state_shift = np.append( invalid_shift, invalid_shift[-1] )
1556
1557        rotated_states = np.empty(states.size, dtype=states.dtype)
1558        rotated_states[:-1] = out_origins
1559        rotated_states[-1]  = out_targets[-1]
1560
1561        rotated_states = np.where( invalid_state_shift, states, rotated_states )
1562
1563        return {
1564            "symbol_index": final_symbols.astype(inputs.dtype),
1565            "state_index":  rotated_states,
1566        }

Generates a new sequence of symbols that are permuted with the symbols emitted by isomorphic states, if they exists.

$ \sigma_o = \mathcal{S}\left[\texttt{input_state_indices}[i]\right] $
$ \sigma_t = \mathcal{S}\left[\texttt{input_state_indices}[i+1]\right] $

$ \mathcal{I}(\sigma_o) = \{\sigma^0_o,\, \sigma^1_o,\, \dots,\, \sigma^{n-1}_o \} $
$ \mathcal{I}(\sigma_t) = \{\sigma^0_t,\, \sigma^1_t,\, \dots,\, \sigma^{n-1}_t \} $

$ k = \bigl(\mathcal{I}(\sigma_o).\texttt{index}(\sigma_o) + \texttt{shift}\bigr) \bmod n $

$ \texttt{output_symbol_indices}[i] := T(\sigma_o^k,\, \sigma_t^k).\text{symbol_index} $
$ \texttt{output_state_indices}[i] := \mathcal{S}.\texttt{index}(\sigma_o^k) $
$ \texttt{output_state_indices}[i+1] := \mathcal{S}.\texttt{index}(\sigma_t^k) $

Where $ \mathcal{I}(\sigma) $` is the ordered set of states isomorphic to $ \sigma $ including $ \sigma $ itself.

Arguments:

input_symbol_indices (np.ndarray): The sequence of generated symbols.
input_state_indices (np.ndarray): The sequence of states that generated symbols with the final state at the end.
shift : int: How much to shift the symbols across the isomorphic states.

def generate_belief_trajectory_from(self, symbols: numpy.ndarray): View Source

1568    def generate_belief_trajectory_from( 
1569        self, 
1570        symbols : np.ndarray ) :
1571
1572        """
1573        Tracks belief states as the symbols are observed (Bayesian updates on state probability distribution).
1574
1575        Returns:
1576            np.ndarray: The belief state sequence. 
1577        """
1578
1579        T_x = self.get_T_X()
1580        pi  = self.get_stationary_distribution()
1581
1582        mu        = pi.copy()
1583        states    = np.zeros(( len(symbols) + 1, len(mu)))
1584        states[0] = mu
1585        
1586        for i, x in enumerate( symbols ) :
1587            mu = mu @ T_x[x]
1588            mu = mu / mu.sum()
1589            states[i+1] = mu
1590        return states

Tracks belief states as the symbols are observed (Bayesian updates on state probability distribution).

Returns:

np.ndarray: The belief state sequence.

def generate_belief_trajectory(self, n_steps: int, random_seed: int = 42) -> numpy.ndarray: View Source

1592    def generate_belief_trajectory( self, 
1593        n_steps : int, 
1594        random_seed : int=42 ) -> np.ndarray:
1595        """
1596        Generates `n_steps` symbols from the HMM, then tracks belief states as the symbols are observed (Bayesian updates on state probability distribution).
1597
1598        Returns:
1599            np.ndarray: The belief state sequence. 
1600        """
1601
1602        trs = self.get_transition_list()
1603
1604        data = am_fast.generate_data(
1605            n_gen=n_steps,
1606            start_state=self.start_state,
1607            transitions=trs,
1608            alphabet=sorted(list(self.alphabet)),
1609            include_states=False,
1610            random_seed=random_seed
1611        )
1612
1613        return self.generate_belief_trajectory_from( data["symbol_index"] )

Generates n_steps symbols from the HMM, then tracks belief states as the symbols are observed (Bayesian updates on state probability distribution).

Returns:

np.ndarray: The belief state sequence.

def generate_data( self, file_prefix: str, n_gen: int, include_states: bool, row_size: int | None = None, include_belief_states: bool = False, isomorphic_shifts: set[int] | None = None, with_component_map: bool = True, random_seed: int = 42) -> dict[str, typing.Any]: View Source

1615    def generate_data(
1616        self,
1617        file_prefix: str,
1618        n_gen: int,
1619        include_states: bool,
1620        row_size : int | None = None,
1621        include_belief_states : bool=False,
1622        isomorphic_shifts : set[int] | None = None,
1623        with_component_map : bool = True, 
1624        random_seed : int=42 ) -> dict[str,Any] : 
1625
1626        if isomorphic_shifts is not None and not include_states :
1627            raise ValueError( "Isomorphic inversion requires include_states=True" )
1628
1629        trs = self.get_transition_list()
1630
1631        data = am_fast.generate_data(
1632            n_gen=n_gen,
1633            start_state=self.start_state,
1634            transitions=trs,
1635            alphabet=sorted(list(self.alphabet)),
1636            include_states=include_states,
1637            random_seed=random_seed
1638        )
1639
1640        if isomorphic_shifts is not None :
1641
1642            data[ "isomorphic_shifts" ] = {}
1643
1644            for shift in isomorphic_shifts :
1645
1646                try : 
1647
1648                    shifted = self.isomorphic_shift(
1649                        input_symbol_indices=data[ "symbol_index" ], 
1650                        input_state_indices=data[ "state_index" ], 
1651                        shift=shift
1652                    )
1653
1654                    data[ "isomorphic_shifts" ][ f"{shift}" ] = {
1655                        "symbol_index" : shifted[ "symbol_index" ],
1656                        "state_index"  : shifted[ "state_index" ]
1657                    }
1658
1659                except Exception as e :
1660                    print( f"Exception {e}" )
1661
1662        if include_belief_states :
1663            belief_states = self.generate_belief_trajectory_from( data[ "symbol_index" ] )
1664            data[ "belief_states" ] = belief_states
1665
1666        metadata = self.get_metadata()
1667
1668        # Maps global state index to the component instance id the state belongs to at 
1669        # a given depth in the composition. 
1670        # somewhat hackish, namely the string based encoding, parsing, and fact that 
1671        # class cmi is only added in compositionally constructed machines
1672        if with_component_map :
1673            
1674            cmp_map = defaultdict(dict)
1675            all_cmis : set[tuple[str,str]] = set()
1676            
1677            for state_index, state in enumerate( self.states ) : 
1678                for level, l_classes in state.classes.items() :
1679                    
1680                    cmi = None
1681
1682                    level_str = str(level)
1683                    state_idx_str = str(state_index)
1684
1685                    for cls in l_classes : 
1686                        cls_tpy = cls.split( "_" )[ 0 ] 
1687                        if cls_tpy == "cmi" :
1688                            cmi = cls
1689                            break
1690                    
1691                    if cmi is not None : 
1692                        cmp_map[ state_idx_str ][ level_str ] = cmi
1693
1694                    if state_idx_str in cmp_map :
1695                        all_cmis.add( ( level_str, cmp_map[ state_idx_str ][ level_str ] ) )
1696            
1697            cmi_indices : dict[ str, dict[str, int ] ] = defaultdict(dict)
1698            for idx, cmi in enumerate( sorted( all_cmis ) ) :
1699                cmi_indices[ cmi[0] ][ cmi[1] ] = idx
1700
1701            metadata["state_component_map" ] = cmp_map
1702            metadata["all_components"      ] = sorted(all_cmis)
1703            metadata["component_indices"   ] = cmi_indices
1704
1705        if "state_index" in data :
1706            index_histogram(
1707                data=data["state_index"], 
1708                output_path=file_prefix + "_state_frequency",
1709                title="State Frequency",
1710                x_label="state index",
1711                show=False
1712            )
1713
1714        am_fast.save_data(
1715            data=data,
1716            file_prefix=file_prefix,
1717            alphabet=sorted(list(self.alphabet)),
1718            n_states=len( self.states ),
1719            start_state=self.start_state,
1720            random_seed=random_seed,
1721            row_size=row_size,
1722            machine_metadata=metadata )
1723
1724        return data

def as_digraph(self) -> networkx.classes.digraph.DiGraph: View Source

1730    def as_digraph( self ) -> nx.DiGraph :
1731
1732        """
1733        Builds a [networkx.DiGraph](https://networkx.org/documentation/stable/reference/classes/digraph.html) constructed from the machine's transitions with no edge symbols or weights.
1734
1735        Returns:
1736        
1737            networkx.DiGraph : the resulting graph.
1738        """
1739
1740        G = nx.DiGraph()
1741        G.add_nodes_from( [ i for i, s in enumerate( self.states ) ] )
1742
1743        for tr in self.transitions :
1744            G.add_edge( tr.origin_state_idx, tr.target_state_idx )
1745
1746        return G

Builds a networkx.DiGraph constructed from the machine's transitions with no edge symbols or weights.

Returns:

networkx.DiGraph : the resulting graph.

def as_dfa(self, with_probs: bool): View Source

1748    def as_dfa( self, with_probs : bool ) :
1749
1750        """
1751        Builds an [automata.fa.dfa.DFA](https://caleb531.github.io/automata/api/fa/class-dfa/) constructed from the machine's transitions.
1752
1753        Args:
1754            with_probs (bool): If true the DFA transitions are labeled based on
1755                the symbol of the machines transition concatenated with its
1756                probability, othwise, the only the symbols.
1757
1758        Returns:
1759        
1760            automata.fa.dfa.DFA : the resulting DFA.
1761        """
1762
1763        precision=8
1764
1765        def edge_label( symb, prob ) :
1766            return f"({symb},{round(prob, precision)})"
1767
1768        # Build states, symbols, and transitions 
1769        dfa_states  = { i for i, _ in enumerate( self.states ) }
1770            
1771        if not with_probs :
1772            dfa_symbols = set( { str(t.symbol_idx) for t in self.transitions } )
1773        else : 
1774            dfa_symbols = set( { edge_label( t.symbol_idx, t.prob ) for t in self.transitions } )
1775
1776        dfa_transitions = defaultdict(dict)
1777
1778        if not with_probs :
1779            for t in self.transitions :
1780                dfa_transitions[ t.origin_state_idx ][ t.symbol_idx ] = t.target_state_idx
1781        else :
1782            for t in self.transitions :
1783                dfa_transitions[ t.origin_state_idx ][ edge_label( t.symbol_idx, t.prob ) ] = t.target_state_idx
1784
1785        # Construct the DFA
1786        return DFA(
1787            states=dfa_states,
1788            input_symbols=dfa_symbols,
1789            transitions=dfa_transitions,
1790            initial_state=self.start_state,
1791            allow_partial=True,
1792            final_states={ 
1793                s for s in dfa_states
1794            }
1795        )

Builds an automata.fa.dfa.DFA constructed from the machine's transitions.

Arguments:

with_probs (bool): If true the DFA transitions are labeled based on the symbol of the machines transition concatenated with its probability, othwise, the only the symbols.

Returns:

automata.fa.dfa.DFA : the resulting DFA.