Skip to content

Markov source

A Markov source, a sequence of symbols in which each symbol depends on the one before it.

  • Name: markov_source
  • Field: mathematics/stochastic_processes
  • Formalism: discrete map
  • References: Lindgren 2024

Factory

@register_model(
    field='mathematics/stochastic_processes',
    tags=('markov chain', 'information theory', 'entropy rate', 'discrete map'),
    references=('Lindgren 2024',),
    defaults={'until': 300, 'dt': 1},
)
def markov_source(transition: Matrix = CHAIN, *, symbol: int = 0) -> System:
    """Return a Markov source, a sequence of symbols in which each symbol depends on the one before it.

    The symbol after `i` is `j` with probability `P_ij`. The entropy rate of the source is
    `h = -Σ_i π_i Σ_j P_ij log2 P_ij` bits per symbol, where `π` is the stationary distribution of the chain, and the
    block entropies of the sequence approach it from above.

    Args:
        transition:
            The transition matrix `P`, one row and one column per symbol, with rows that sum to 1.

        symbol:
            The initial symbol.

    Returns:
        system:
            The model, with the symbols `source.symbol`.

    Examples:
        >>> from skcomplex.models.mathematics.stochastic_processes import markov_source
        >>> markov_source().interactions[0].state
        'symbol'
    """
    chain = TransitionInteraction('chain', target='source', state='symbol', matrix=lambda s: transition)
    return System(elements=[MarkovSourceElement('source', symbol=symbol)], interactions=[chain])

Components

class MarkovSourceElement(Element):
    """Source of symbols `0` to `n - 1`, whose symbol a transition matrix steps in discrete time.

    Args:
        label:
            The label of the source.

        visual:
            How the element is drawn, or `None` for `default_visual_`.

        symbol:
            The initial symbol.

    Examples:
        >>> from skcomplex.models.mathematics.stochastic_processes import MarkovSourceElement
        >>> MarkovSourceElement('source', symbol=1).symbol
        1
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        symbol: Annotated[float, StateSpec(domain=NN)] = 0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.symbol = symbol

class TransitionInteraction(Interaction):
    """Move of a discrete state of the target from one value to another, or between all its values by a matrix.

    A member whose state holds `source` moves to `destination` at a `rate`, as a jump in continuous time, with a
    `probability` in every step of discrete time, or in the step in which a `condition` holds. With `over='neighbors'`
    the rate and the probability are those of one neighbour, written with `namespace.target.other` for the neighbour:
    the rates of the neighbours add up, and their probabilities combine as `1 - prod(1 - p)`, the probability that at
    least one of them moves the member. With a matrix instead, the state takes the values `0` to `n - 1` and steps by
    the rows of a transition matrix in discrete time, the row of its current value. An `emission` matrix then draws the
    `observed` state of the target from the row that the next value of the state selects, as a hidden Markov source.
    With a residence time `after`, a member moves by a probability or a condition only once it has held the source
    value for `after` steps, at the earliest at the end of the step that completes them.

    In the `'symbolic'` form, every function reads `namespace.target`, the states and parameters of the target, and the
    time `namespace.t`. In the `'agent'` form, the target is an agent group and every function takes the agent and
    the model view, `f(agent, model)`, and returns a number, a truth value, or the rows of a matrix, and a rate `r`
    moves an agent with probability `1 - exp(-r)` in a step. With `over='cell'`, every member acts on the other members
    of its cell: it moves each one that holds the source value with the probability, the rate, or the condition
    `f(agent, other, model)` of that member `agent` and the acting member `other`, each by its own draw.

    The `schedule` of the agent form orders the members. Under `'synchronous'` every member decides from the states at
    the start of the step, under `'sequential'` the members act one after the other in their order, and under
    `'random-sequential'` in a new random order every step, each reading the values the members before it set. Under
    `'member'` every member acts in its own activation, together with the rules of its element and the other member
    interactions, in the order the system lists them, so one rule of a member runs its movement, transfers,
    transitions, and life cycle.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        target:
            The path of the element or group whose state moves.

        state:
            The name of the discrete state of the target.

        source:
            The value the state moves from, a collection of such values, or a function of the value of the state that
            returns whether the state moves from it, such as `lambda value: value >= 0`, a condition in the symbolic
            form. A state that moves from several values moves to the destination from each.

        destination:
            The value the state moves to.

        rate:
            The function of the rate of the move, or `None`.

        probability:
            The function of the probability of the move in a step, or `None`.

        condition:
            The function of the condition under which the state moves in a step, or `None`.

        matrix:
            The function of the transition matrix, one row and one column per value of the state, or `None`.

        emission:
            With a matrix, the function of the emission matrix, one row per value of the state and one column per value
            of the observed state, or `None`.

        observed:
            The name of the state of the target that the emission matrix draws, or `None`.

        over:
            `'member'` for a rate or a probability of the member, `'neighbors'` for that of one neighbour, or `'cell'`
            for the agent form whose members act on the members of their cell.

        neighborhood:
            The neighbourhood of `over='neighbors'`, or `None` for the neighbours of the network or the lattice.

        after:
            The number of steps a member holds the source value before it can move, or `None`.

        schedule:
            The order of the members of the agent form, `'synchronous'`, `'sequential'`, `'random-sequential'`, or
            `'member'`.

        form:
            `'symbolic'` for laws of the states of the target, or `'agent'` for a rule over an agent group.

        phase:
            The phase of the step the agent form runs in, or `None` for the default phase.

    Raises:
        ValidationError:
            The transition has not exactly one kind, the source and the destination are equal, an emission has no
            matrix or no observed state, `over` is unknown or goes with a form or a kind it cannot take, the
            residence time is negative or goes with a rate, a matrix, or members acting on their cell, the form or the
            schedule is unknown, the symbolic form has a schedule, or a function is not callable.

    Examples:
        >>> from skcomplex.core.interactions import TransitionInteraction
        >>> recovery = TransitionInteraction(
        ...     'recovery', target='hosts', state='state', source=1, destination=0, rate=lambda s: 0.5
        ... )
        >>> recovery.check().over
        'member'
        >>> agents = TransitionInteraction(
        ...     'recovery', target='hosts', state='state', source=1, destination=2, probability=lambda a, m: 0.1,
        ...     after=3, form='agent'
        ... )
        >>> type(agents.check().prepare().run(SimpleNamespace()).laws_[0]).__name__
        'RuleLaw'
        >>> separation = TransitionInteraction(
        ...     'separation', target='workers', state='employer', source=lambda value: value >= 0, destination=-1,
        ...     probability=lambda a, m: 0.05, form='agent'
        ... )
        >>> separation.check().destination
        -1
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        state: str = 'state',
        source: SourceValues = 0.0,
        destination: float = 1.0,
        rate: Callable[..., SupportsFloat] | None = None,
        probability: Callable[..., SupportsFloat] | None = None,
        condition: Callable[..., object] | None = None,
        matrix: Callable[..., Sequence[Sequence[object]]] | None = None,
        emission: Callable[..., Sequence[Sequence[object]]] | None = None,
        observed: str | None = None,
        over: str = 'member',
        neighborhood: BaseNeighborhood | None = None,
        after: float | None = None,
        schedule: str = 'synchronous',
        form: str = 'symbolic',
        phase: str | None = None,
    ) -> None:
        super().__init__(label)
        self.target = target
        self.state = state
        self.source = source
        self.destination = destination
        self.rate = rate
        self.probability = probability
        self.condition = condition
        self.matrix = matrix
        self.emission = emission
        self.observed = observed
        self.over = over
        self.neighborhood = neighborhood
        self.after = after
        self.schedule = schedule
        self.form = form
        self.phase = phase

    def _check_kind(self) -> str:
        """Raise an error for a transition without exactly one callable kind, or with an emission it cannot take."""
        kinds = [kind for kind in KINDS if getattr(self, kind) is not None]
        if len(kinds) != 1 or not callable(getattr(self, kinds[0])):
            error_msg = (
                f'The transition {self.label!r} has the forms {kinds}. Give exactly one of {", ".join(KINDS)}, as a '
                'function of the namespace.'
            )
            raise ValidationError(error_msg)
        unpaired = (self.emission is None) != (self.observed is None)
        if unpaired or (self.emission is not None and self.matrix is None):
            error_msg = (
                f'The transition {self.label!r} has the emission {self.emission!r} of the state {self.observed!r}. '
                'Give both with a matrix, or neither.'
            )
            raise ValidationError(error_msg)
        return kinds[0]

    def _check_form(self) -> None:
        """Raise an error for a transition without exactly one callable kind, or with arguments its kind ignores."""
        kind = self._check_kind()
        if self.over not in OVERS:
            error_msg = f'The transition {self.label!r} runs over {self.over!r}. Use one of {", ".join(OVERS)}.'
            raise ValidationError(error_msg)
        if self.form not in FORMS:
            error_msg = f'The form of transition {self.label!r} is {self.form!r}. Use one of {", ".join(FORMS)}.'
            raise ValidationError(error_msg)
        if self.schedule not in SCHEDULES or (self.form == 'symbolic' and self.schedule != 'synchronous'):
            error_msg = (
                f'The schedule of transition {self.label!r} is {self.schedule!r}. Use one of {", ".join(SCHEDULES)}, '
                'and the agent form for another than `synchronous`.'
            )
            raise ValidationError(error_msg)
        if self.over == 'cell' and (self.form != 'agent' or kind == 'matrix' or self.after is not None):
            error_msg = (
                f'The members of transition {self.label!r} act on their cell, which only the agent form does, by a '
                'rate, a probability, or a condition and without a residence time.'
            )
            raise ValidationError(error_msg)
        if self.over == 'neighbors' and (kind not in {'rate', 'probability'} or self.form == 'agent'):
            error_msg = (
                f'The transition {self.label!r} runs over neighbours, which only a rate or a probability of the '
                'symbolic form does. In the agent form, read the neighbours of the agent in the function.'
            )
            raise ValidationError(error_msg)
        single = not (callable(self.source) or isinstance(self.source, Collection))
        if self.matrix is None and single and self.source == self.destination:
            error_msg = f'The transition {self.label!r} moves from {self.source} to itself. Give distinct values.'
            raise ValidationError(error_msg)
        if self.after is not None and (kind not in {'probability', 'condition'} or not self.after >= 0):
            error_msg = (
                f'The transition {self.label!r} has the residence time {self.after!r}. Give a number of steps at least '
                '0 with a probability or a condition.'
            )
            raise ValidationError(error_msg)

    def _check_transition(self) -> None:
        """Raise an error for arguments the transition cannot take."""
        self._check_form()

    def _check(self) -> None:
        """Raise an error for a transition without exactly one kind, or with arguments its kind ignores.

        Raises:
            ValidationError:
                The transition has not exactly one kind, the source and the destination are equal, an emission has
                no matrix or no observed state, `over` is unknown or goes with a condition, a matrix, or the agent
                form, the residence time is negative or goes with a rate or a matrix, the form is unknown, or a
                function is not callable.
        """
        super()._check()
        self._check_transition()

    def _get_moving(self, value: Expression, source: SourceValues) -> Expression:
        """Return 1 for a member whose state holds a source value, and 0 otherwise."""
        if callable(source):
            return where(SR(source(value)), 1, 0)
        if isinstance(source, Collection):
            return sum((where(value == item, 1, 0) for item in set(source)), SR(0))
        return where(value == source, 1, 0)

    def _get_jump(self, value: Expression, source: SourceValues, destination: float, rate: Expression) -> JumpLaw:
        """Return the jump of the state from a source value to the destination at a rate."""
        if callable(source) or isinstance(source, Collection):
            return JumpLaw(rate * self._get_moving(value, source), {value: destination - value})
        return JumpLaw(rate * self._get_moving(value, source), {value: destination - source})

    def _get_step(self, value: Expression, source: SourceValues, destination: float, moves: Expression) -> UpdateLaw:
        """Return the update of the state that moves a member from a source value where a condition holds."""
        if callable(source) or isinstance(source, Collection):
            return UpdateLaw(value, where(self._get_moving(value, source) > 0, where(moves, destination, value), value))
        return UpdateLaw(value, where(value == source, where(moves, destination, source), value), source=source)

    def _holds_source(self, value: float) -> bool:
        """Return whether the value of the state of an agent is a source value."""
        if callable(self.source):
            return bool(self.source(value))
        if isinstance(self.source, Collection):
            return value in set(self.source)
        return value == self.source

    def _combine_rates(self, rate: Expression) -> Expression:
        """Return the rate of the member, or the sum of the rates of its neighbours."""
        return neighbors(rate, over=self.neighborhood) if self.over == 'neighbors' else rate

    def _combine_probabilities(self, probability: Expression) -> Expression:
        """Return the probability of the member, or the probability that at least one neighbour moves it."""
        if self.over == 'member':
            return probability
        return 1 - exp(neighbors(log(1 - probability), over=self.neighborhood))

    def _get_draw(self) -> Expression:
        """Return the uniform draw of the transition in a step."""
        return uniform(re.sub(r'\W', '_', f'transition_{self.label}'))

    def _get_chain(self, namespace: SimpleNamespace, value: Expression) -> list[BasePart]:
        """Return the transition law of the state by its matrix, with the emission law of the observed state."""
        laws: list[BasePart] = []
        if self.matrix is not None:
            laws.append(TransitionLaw(value, self.matrix(namespace), label=self.label))
        if self.emission is not None and self.observed is not None:
            observed = getattr(self._get_member(namespace, 'target'), self.observed)
            laws.append(EmissionLaw(observed, value, self.emission(namespace), label=self.label))
        return laws

    def _get_held(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> int:
        """Return the number of whole steps an agent has held its value, counted in the storage of the model."""
        value = float(getattr(agent, self.state))
        record = model.storage.get(f'{self.label}.{agent.id}')
        held = record[1] + 1 if isinstance(record, tuple) and record[0] == value else 0
        model.storage[f'{self.label}.{agent.id}'] = (value, held)
        return held

    def _draw_row(self, row: Sequence[object], model: ModelViewProtocol) -> float:
        """Return a value drawn by the probabilities of a row of a matrix."""
        weights = np.asarray(row, dtype=float)
        return float(model.rng.choice(weights.size, p=weights / weights.sum()))

    def _evaluate[T](
        self,
        function: Callable[..., T],
        agent: AgentViewProtocol,
        other: AgentViewProtocol | None,
        model: ModelViewProtocol,
    ) -> T:
        """Return a function of an agent, or of an agent and the member of its cell that acts on it."""
        return function(agent, model) if other is None else function(agent, other, model)

    def _moves_agent(
        self,
        agent: AgentViewProtocol,
        model: ModelViewProtocol,
        other: AgentViewProtocol | None = None,
    ) -> bool:
        """Return whether an agent that holds the source value moves in this step, moved by another one or alone."""
        if self.condition is not None:
            return bool(self._evaluate(self.condition, agent, other, model))
        if self.rate is not None:
            probability = 1.0 - math.exp(-float(self._evaluate(self.rate, agent, other, model)))
        elif self.probability is not None:
            probability = float(self._evaluate(self.probability, agent, other, model))
        else:
            probability = 0.0
        return bool(model.rng.random() < probability)

    def _get_next(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> float | None:
        """Return the next value of the state of an agent, or `None` when it keeps its value."""
        value = float(getattr(agent, self.state))
        if self.matrix is not None:
            return self._draw_row(self.matrix(agent, model)[int(value)], model)
        held = self._get_held(agent, model) if self.after is not None else 0
        if not self._holds_source(value) or (self.after is not None and held + 1 < self.after):
            return None
        return float(self.destination) if self._moves_agent(agent, model) else None

    def _get_moves(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> list[tuple[AgentViewProtocol, float]]:
        """Return the members an agent moves in this step, itself or the members of its cell, with their values."""
        if self.over != 'cell':
            value = self._get_next(agent, model)
            return [] if value is None else [(agent, value)]
        return [
            (other, float(self.destination))
            for other in agent.agents_here(agent.group)
            if other is not agent
            and self._holds_source(float(getattr(other, self.state)))
            and self._moves_agent(other, model, agent)
        ]

    def _apply(self, moves: list[tuple[AgentViewProtocol, float]], model: ModelViewProtocol) -> None:
        """Set the moved states, and draw the observed state of an emission."""
        for agent, value in moves:
            setattr(agent, self.state, value)
            if self.emission is not None and self.observed is not None:
                setattr(agent, self.observed, self._draw_row(self.emission(agent, model)[int(value)], model))

    def _move_member(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
        """Move the state of one member, or of the members of its cell, from the states the members before it set."""
        if agent.alive:
            self._apply(self._get_moves(agent, model), model)

    def move(self, model: ModelViewProtocol) -> None:
        """Move the state of every agent of the target that moves in this step, by the schedule.

        Args:
            model:
                The model view of the interaction.
        """
        agents = list(model.agents(self._get_role('target')))
        if self.schedule == 'synchronous':
            moves = [move for agent in agents for move in self._get_moves(agent, model)]
            self._apply(moves, model)
            return
        order = model.rng.permutation(len(agents)) if self.schedule == 'random-sequential' else range(len(agents))
        for index in order:
            self._move_member(agents[int(index)], model)

    def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
        """Return the jump, the update, the transition or emission law of the state, or the rule of the agent form.

        Args:
            namespace:
                The namespace of the symbols the laws can use, with `namespace.target`.

        Returns:
            laws:
                The law of the move, or the rule that moves the agents.
        """
        if self.form == 'agent' and self.schedule == 'member':
            return [RuleLaw(self._move_member, phase=self.phase, members=self._get_role('target'))]
        if self.form == 'agent':
            return [RuleLaw(self.move, phase=self.phase)]
        value = getattr(self._get_member(namespace, 'target'), self.state)
        if self.matrix is not None:
            return self._get_chain(namespace, value)
        if self.rate is not None:
            rate = self._combine_rates(SR(self.rate(namespace)))
            return [self._get_jump(value, self.source, self.destination, rate)]
        if self.probability is not None:
            probability = self._combine_probabilities(SR(self.probability(namespace)))
            moves = self._get_draw() < probability
        else:
            moves = SR(self.condition(namespace)) if self.condition is not None else SR(0)
        if self.after is not None:
            moves = both(residence(value) + 1 >= self.after, moves)
        return [self._get_step(value, self.source, self.destination, moves)]

Equations

source.symbol: symbol_source(t + 1) == where(categorical__transition_symbol_source_chain*(where(symbol_source(t) == 0, 0.9, 0.4) + where(symbol_source(t) == 0, 0.1, 0.6)) >= where(symbol_source(t) == 0, 0.9, 0.4), 1, 0)