Reference

qmm.core.structure

Define model structure in graph, matrix or equation forms.

import_digraph

import_digraph(data, file_path=True)

Import a JSON model and convert to a NetworkX DiGraph with sign attributes.

Parameters:

Name	Type	Description	Default
data	Union[str, dict]	Path to JSON file or dictionary containing model structure	required
file_path	bool	If True, data is a file path. If False, data is a dictionary	True

Returns:

Type	Description
DiGraph	nx.DiGraph: Signed directed graph (signed digraph)

Examples:

from qmm import import_digraph
model_dict = {
"nodes": [{"id": "A"}, {"id": "B"}],
"edges": [{"from": "A", "to": "B", "arrows": {"to": {"type": "triangle"}}}]
}
list(import_digraph(model_dict, file_path=False).nodes())
# ['A', 'B']

G = import_digraph(model_dict, file_path=False)
list(G.edges(data='sign'))
# [('A', 'B', 1)]

Source code in qmm/core/structure.py

def import_digraph(data: Union[str, dict], file_path: bool = True) -> nx.DiGraph:
    """Import a JSON model and convert to a NetworkX DiGraph with sign attributes.

    Args:
        data: Path to JSON file or dictionary containing model structure
        file_path: If True, data is a file path. If False, data is a dictionary

    Returns:
        nx.DiGraph: Signed directed graph (signed digraph)

    Examples:
        ```python
        from qmm import import_digraph
        model_dict = {
        "nodes": [{"id": "A"}, {"id": "B"}],
        "edges": [{"from": "A", "to": "B", "arrows": {"to": {"type": "triangle"}}}]
        }
        list(import_digraph(model_dict, file_path=False).nodes())
        # ['A', 'B']

        G = import_digraph(model_dict, file_path=False)
        list(G.edges(data='sign'))
        # [('A', 'B', 1)]
        ```
    """
    if file_path:
        with open(data, "r") as file:
            data = json.load(file)
    G = nx.DiGraph()
    for node in data["nodes"]:
        att = {k: v for k, v in node.items() if k != "id"}
        if "title" not in att:
            att["title"] = None
        G.add_node(str(node["id"]), **att)
    for edge in data["edges"]:
        source, target = str(edge["from"]), str(edge["to"])
        att = {k: v for k, v in edge.items() if k not in ["from", "to", "arrows"]}
        arr = edge.get("arrows", {}).get("to", {})
        if isinstance(arr, dict):
            arr_type = arr.get("type")
            if arr_type == "triangle":
                att["sign"] = 1
            elif arr_type == "circle":
                att["sign"] = -1
        if "dashes" not in att:
            att["dashes"] = False
        if "title" not in att:
            att["title"] = None
        G.add_edge(source, target, **att)
    nx.set_node_attributes(G, "state", "category")
    nx.freeze(G)
    return G

create_matrix

create_matrix(G, form='symbolic', matrix_type='A')

Create an interaction matrix from a signed digraph in symbolic, signed, or binary form.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing a signed digraph model	required
form	Literal['symbolic', 'signed', 'binary']	Type of matrix elements ('symbolic', 'signed', or 'binary')	'symbolic'
matrix_type	Literal['A', 'B', 'C', 'D']	Type of matrix to create ('A', 'B', 'C', or 'D')	'A'

Returns:

Type	Description
Matrix	sp.Matrix: Interaction matrix

References

• Levins, R. (1968). Evolution in changing environments; some theoretical explorations.
• Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.

Examples:

from qmm import load_digraph, create_matrix
create_matrix(load_digraph("snowshoe_rp"), form='symbolic')
# Matrix([
# [-a_R,R, -a_R,C,      0],
# [ a_C,R,      0, -a_C,P],
# [ a_P,R,  a_P,C, -a_P,P]])

create_matrix(load_digraph("snowshoe_rp"), form='signed')
# Matrix([
# [-1, -1,  0],
# [ 1,  0, -1],
# [ 1,  1, -1]])

create_matrix(load_digraph("snowshoe_rp"), form='binary')
# Matrix([
# [1, 1, 0],
# [1, 0, 1],
# [1, 1, 1]])

Source code in qmm/core/structure.py

def create_matrix(
    G: nx.DiGraph,
    form: Literal["symbolic", "signed", "binary"] = "symbolic",
    matrix_type: Literal["A", "B", "C", "D"] = "A",
) -> sp.Matrix:
    """Create an interaction matrix from a signed digraph in symbolic, signed, or binary form.

    Args:
        G: NetworkX DiGraph representing a signed digraph model
        form: Type of matrix elements ('symbolic', 'signed', or 'binary')
        matrix_type: Type of matrix to create ('A', 'B', 'C', or 'D')

    Returns:
        sp.Matrix: Interaction matrix

    References:
        - Levins, R. (1968). Evolution in changing environments; some theoretical explorations.
        - Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.

    Examples:
        ```python
        from qmm import load_digraph, create_matrix
        create_matrix(load_digraph("snowshoe_rp"), form='symbolic')
        # Matrix([
        # [-a_R,R, -a_R,C,      0],
        # [ a_C,R,      0, -a_C,P],
        # [ a_P,R,  a_P,C, -a_P,P]])

        create_matrix(load_digraph("snowshoe_rp"), form='signed')
        # Matrix([
        # [-1, -1,  0],
        # [ 1,  0, -1],
        # [ 1,  1, -1]])

        create_matrix(load_digraph("snowshoe_rp"), form='binary')
        # Matrix([
        # [1, 1, 0],
        # [1, 0, 1],
        # [1, 1, 1]])
        ```
    """

    def sym(source: str, target: str, prefix: str) -> sp.Symbol:
        return sp.Symbol(f"{prefix}_{target},{source}")

    def sign(source: str, target: str, prefix: str) -> Union[sp.Symbol, int]:
        if form == "symbolic":
            return sym(source, target, prefix) * G[source][target].get("sign", 1)
        elif form == "signed":
            return G[source][target].get("sign", 1)
        else:
            return int(G.has_edge(source, target))

    def product(path: List[str]) -> Union[sp.Symbol, int]:
        effect = 1
        for i in range(len(path) - 1):
            effect *= sign(path[i], path[i + 1], prefix)
        return effect

    state_n = get_nodes(G, "state")
    input_n = get_nodes(G, "input")
    output_n = get_nodes(G, "output")
    matrix_configs: Dict[str, Tuple[List[str], List[str], str, str]] = {
        "A": (state_n, state_n, "a", "state"),
        "B": (state_n, input_n, "b", "input"),
        "C": (output_n, state_n, "c", "output"),
        "D": (output_n, input_n, "d", "input"),
    }
    rows, cols, prefix, category = matrix_configs[matrix_type]
    matrix = sp.zeros(len(rows), len(cols))
    if matrix_type == "D":
        _check_direct_io_edges(G)
        return matrix
    for i, target in enumerate(rows):
        for j, source in enumerate(cols):
            if matrix_type == "A":
                if G.has_edge(source, target):
                    matrix[i, j] = sign(source, target, prefix)
            else:
                paths = nx.all_simple_paths(G, source, target)
                valid = [p for p in paths if all(G.nodes[n]["category"] == category for n in p[1:-1])]
                matrix[i, j] = sum(product(path) for path in valid)
    return matrix

create_equations

create_equations(G, form='state')

Create linear system of differential equations from a signed digraph.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing a signed digraph model	required
form	Literal['state', 'output']	Type of equations to create ('state' or 'output')	'state'

Returns:

Type	Description
Matrix	sp.Matrix: Linear system of differential equations

References

Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.

Examples:

from qmm import load_digraph, create_equations
create_equations(load_digraph("snowshoe_rp"), form='state')
# Matrix([
# [           -a_R,C*x_C - a_R,R*x_R],
# [           -a_C,P*x_P + a_C,R*x_R],
# [a_P,C*x_C - a_P,P*x_P + a_P,R*x_R]])

Source code in qmm/core/structure.py

def create_equations(G: nx.DiGraph, form: Literal["state", "output"] = "state") -> sp.Matrix:
    """Create linear system of differential equations from a signed digraph.

    Args:
        G: NetworkX DiGraph representing a signed digraph model
        form: Type of equations to create ('state' or 'output')

    Returns:
        sp.Matrix: Linear system of differential equations

    References:
        Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.

    Examples:
        ```python
        from qmm import load_digraph, create_equations
        create_equations(load_digraph("snowshoe_rp"), form='state')
        # Matrix([
        # [           -a_R,C*x_C - a_R,R*x_R],
        # [           -a_C,P*x_P + a_C,R*x_R],
        # [a_P,C*x_C - a_P,P*x_P + a_P,R*x_R]])
        ```
    """
    A = create_matrix(G, form="symbolic", matrix_type="A")
    B = create_matrix(G, form="symbolic", matrix_type="B")
    C = create_matrix(G, form="symbolic", matrix_type="C")
    D = create_matrix(G, form="symbolic", matrix_type="D")
    state_nodes = get_nodes(G, "state")
    input_nodes = get_nodes(G, "input")
    output_nodes = get_nodes(G, "output")
    x = sp.Matrix([sp.Symbol(f"x_{i}") for i in state_nodes])
    u = sp.Matrix([sp.Symbol(f"u_{i}") for i in input_nodes]) if input_nodes else None
    if form == "state":
        equations = A * x
        if B.shape[1] > 0 and u is not None:
            equations += B * u
        return equations
    if form != "output":
        raise ValueError("form must be either 'state' or 'output'")
    if not output_nodes:
        raise ValueError("No output nodes found in graph")
    equations = C * x
    if D.shape[1] > 0 and u is not None:
        equations += D * u
    return equations

nodes_table

nodes_table(G)

Create a table of node metadata.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Node, label, category, and description

Examples:

from qmm import load_digraph, nodes_table
nodes_table(load_digraph("snowshoe_rp"))
#       Node Label Category Description
# 0  $x_{R}$     R    State        None
# 1  $x_{C}$     C    State        None
# 2  $x_{P}$     P    State        None

Source code in qmm/core/structure.py

def nodes_table(G: nx.DiGraph) -> pd.DataFrame:
    """Create a table of node metadata.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        pd.DataFrame: Node, label, category, and description

    Examples:
        ```python
        from qmm import load_digraph, nodes_table
        nodes_table(load_digraph("snowshoe_rp"))
        #       Node Label Category Description
        # 0  $x_{R}$     R    State        None
        # 1  $x_{C}$     C    State        None
        # 2  $x_{P}$     P    State        None
        ```
    """
    rows = []
    for node_id, data in G.nodes(data=True):
        category = data.get("category", "state")
        if category == "input":
            symbol = f"u_{{{node_id}}}"
        elif category == "output":
            symbol = f"y_{{{node_id}}}"
        else:
            symbol = f"x_{{{node_id}}}"
        if category == "input":
            category_label = "Input"
        elif category == "output":
            category_label = "Output"
        else:
            category_label = "State"
        rows.append(
            {
                "Node": f"${symbol}$",
                "Label": data.get("label", str(node_id)),
                "Category": category_label,
                "Description": data.get("title"),
            }
        )
    return pd.DataFrame(rows)

edges_table

edges_table(G)

Create a table of edge metadata.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Edge, from, sign, to, dashes, and description

Examples:

from qmm import load_digraph, edges_table
edges_table(load_digraph("snowshoe_rp"))
#         Edge From Sign To  Dashes Description
# 0  $a_{R,R}$    R    -  R   False        None
# 1  $a_{C,R}$    R    +  C   False        None
# 2  $a_{P,R}$    R    +  P   False        None
# 3  $a_{R,C}$    C    -  R   False        None
# 4  $a_{P,C}$    C    +  P   False        None
# 5  $a_{C,P}$    P    -  C   False        None
# 6  $a_{P,P}$    P    -  P   False        None

Source code in qmm/core/structure.py

def edges_table(G: nx.DiGraph) -> pd.DataFrame:
    """Create a table of edge metadata.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        pd.DataFrame: Edge, from, sign, to, dashes, and description

    Examples:
        ```python
        from qmm import load_digraph, edges_table
        edges_table(load_digraph("snowshoe_rp"))
        #         Edge From Sign To  Dashes Description
        # 0  $a_{R,R}$    R    -  R   False        None
        # 1  $a_{C,R}$    R    +  C   False        None
        # 2  $a_{P,R}$    R    +  P   False        None
        # 3  $a_{R,C}$    C    -  R   False        None
        # 4  $a_{P,C}$    C    +  P   False        None
        # 5  $a_{C,P}$    P    -  C   False        None
        # 6  $a_{P,P}$    P    -  P   False        None
        ```
    """
    rows = []
    for source, target, data in G.edges(data=True):
        prefix = _edge_prefix(G, source, target)
        sign_val = data.get("sign", 1)
        if sign_val == 1:
            sign_label = "+"
        elif sign_val == -1:
            sign_label = "-"
        else:
            sign_label = str(sign_val)
        rows.append(
            {
                "Edge": f"${prefix}_{{{target},{source}}}$",
                "From": G.nodes[source].get("label", str(source)),
                "Sign": sign_label,
                "To": G.nodes[target].get("label", str(target)),
                "Dashes": data.get("dashes", False),
                "Description": data.get("title"),
            }
        )
    return pd.DataFrame(rows)

qmm.core.stability

Analyse the stability properties of a system based on its structure.

sign_stability

sign_stability(G)

Evaluate necessary and sufficient conditions for sign stability including color test.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Test results for sign stability conditions

References

• May, R.M. (1973). Qualitative Stability in Model Ecosystems. Ecology 54, 638–641.
• Jeffries, C. (1974). Qualitative Stability and Digraphs in Model Ecosystems. Ecology 55, 1415–1419.

Examples:

from qmm import load_digraph, sign_stability
sign_stability(load_digraph("snowshoe_rp"))
#             Test                                                                     Definition  Result
# 0    Condition i                                                       No positive self-effects    True
# 1   Condition ii                                           At least one node is self-regulating    True
# 2  Condition iii                        The product of any pairwise interaction is non-positive    True
# 3   Condition iv                                              No cycles greater than length two   False
# 4    Condition v  Non-zero determinant (all nodes have at least one incoming and outgoing link)    True
# 5    Colour test                                                    Fails Jeffries' colour test    True
# 6    Sign stable               Satisfies necessary and sufficient conditions for sign stability   False

Source code in qmm/core/stability.py

def sign_stability(G: nx.DiGraph) -> pd.DataFrame:
    """Evaluate necessary and sufficient conditions for sign stability including color test.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        pd.DataFrame: Test results for sign stability conditions

    References:
        - May, R.M. (1973). Qualitative Stability in Model Ecosystems. Ecology 54, 638–641.
        - Jeffries, C. (1974). Qualitative Stability and Digraphs in Model Ecosystems. Ecology 55, 1415–1419.

    Examples:
        ```python
        from qmm import load_digraph, sign_stability
        sign_stability(load_digraph("snowshoe_rp"))
        #             Test                                                                     Definition  Result
        # 0    Condition i                                                       No positive self-effects    True
        # 1   Condition ii                                           At least one node is self-regulating    True
        # 2  Condition iii                        The product of any pairwise interaction is non-positive    True
        # 3   Condition iv                                              No cycles greater than length two   False
        # 4    Condition v  Non-zero determinant (all nodes have at least one incoming and outgoing link)    True
        # 5    Colour test                                                    Fails Jeffries' colour test    True
        # 6    Sign stable               Satisfies necessary and sufficient conditions for sign stability   False
        ```
    """
    A_signed = create_matrix(G, form="signed")
    A = sp.matrix2numpy(A_signed).astype(int)
    n = A.shape[0]
    conditions = [
        all(A[i, i] <= 0 for i in range(n)),
        any(A[i, i] < 0 for i in range(n)),
        all(A[i, j] * A[j, i] <= 0 for i in range(n) for j in range(n) if i != j),
        all(len(cycle) < 3 for cycle in nx.simple_cycles(nx.DiGraph(A))),
        bool(A_signed.det() != 0),
    ]
    colour_result = _colour_test(G) == "Fail"
    is_sign_stable = all(conditions) and colour_result
    return pd.DataFrame(
        {
            "Test": [
                "Condition i",
                "Condition ii",
                "Condition iii",
                "Condition iv",
                "Condition v",
                "Colour test",
                "Sign stable",
            ],
            "Definition": [
                "No positive self-effects",
                "At least one node is self-regulating",
                "The product of any pairwise interaction is non-positive",
                "No cycles greater than length two",
                "Non-zero determinant (all nodes have at least " + "one incoming and outgoing link)",
                "Fails Jeffries' colour test",
                "Satisfies necessary and sufficient conditions for sign stability",
            ],
            "Result": conditions + [colour_result] + [is_sign_stable],
        }
    )

system_feedback cached

system_feedback(G, level=None, form='symbolic')

Calculate the product of conjunct and disjunct feedback cycles for any level of the system (coefficients of the characteristic polynomial).

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Level of feedback to compute (None for all levels)	None
form	Literal['symbolic', 'signed']	Type of feedback ('symbolic', 'signed')	'symbolic'

Returns:

Type	Description
Matrix	sp.Matrix: Feedback cycle products at specified levels

References

• Lyapunov, A.M. (1892). The general problem of the stability of motion. International Journal of Control 55, 531–534.
• Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
• Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

Examples:

from qmm import load_digraph, system_feedback
system_feedback(load_digraph("snowshoe_rp"), level=2, form='symbolic')
# Matrix([[-a_C,P*a_P,C - a_C,R*a_R,C - a_P,P*a_R,R]])

system_feedback(load_digraph("snowshoe_rp"), level=3, form='symbolic')
# Matrix([[-a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C - a_C,R*a_P,P*a_R,C]])

system_feedback(load_digraph("snowshoe_rp"), form='symbolic')
# Matrix([
# [                                                        -1],
# [                                            -a_P,P - a_R,R],
# [                  -a_C,P*a_P,C - a_C,R*a_R,C - a_P,P*a_R,R],
# [-a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C - a_C,R*a_P,P*a_R,C]])

Source code in qmm/core/stability.py

@cache
def system_feedback(
    G: nx.DiGraph,
    level: Optional[int] = None,
    form: Literal["symbolic", "signed"] = "symbolic",
) -> sp.Matrix:
    """Calculate the product of conjunct and disjunct feedback cycles for any level of the system (coefficients of the characteristic polynomial).

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Level of feedback to compute (None for all levels)
        form: Type of feedback ('symbolic', 'signed')

    Returns:
        sp.Matrix: Feedback cycle products at specified levels

    References:
        - Lyapunov, A.M. (1892). The general problem of the stability of motion. International Journal of Control 55, 531–534.
        - Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
        - Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

    Examples:
        ```python
        from qmm import load_digraph, system_feedback
        system_feedback(load_digraph("snowshoe_rp"), level=2, form='symbolic')
        # Matrix([[-a_C,P*a_P,C - a_C,R*a_R,C - a_P,P*a_R,R]])

        system_feedback(load_digraph("snowshoe_rp"), level=3, form='symbolic')
        # Matrix([[-a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C - a_C,R*a_P,P*a_R,C]])

        system_feedback(load_digraph("snowshoe_rp"), form='symbolic')
        # Matrix([
        # [                                                        -1],
        # [                                            -a_P,P - a_R,R],
        # [                  -a_C,P*a_P,C - a_C,R*a_R,C - a_P,P*a_R,R],
        # [-a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C - a_C,R*a_P,P*a_R,C]])
        ```
    """
    A = create_matrix(G, form=form)
    if level == 0:
        return sp.Matrix([-1])
    n = A.shape[0]
    if form not in ("symbolic", "signed"):
        raise ValueError("form must be either 'symbolic' or 'signed'")
    if level is not None and (level < 0 or level > n):
        raise ValueError(f"Level must be between 0 and {n}")
    lam = sp.symbols("lambda")
    p = A.charpoly(lam).as_expr()
    if level is None:
        fb = [-p.coeff(lam, n - k) for k in range(n + 1)]
    else:
        fb = [-p.coeff(lam, n - level)]
    return sp.Matrix(fb)

net_feedback cached

net_feedback(G, level=None)

Calculate net feedback at a specified level of the system.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Level of feedback to compute (None for all levels)	None

Returns:

Type	Description
Matrix	sp.Matrix: Net feedback at specified levels

References

• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, net_feedback
net_feedback(load_digraph("snowshoe_rp"), level=2)
# Matrix([[-3]])

net_feedback(load_digraph("snowshoe_rp"))
# Matrix([
# [-1],
# [-2],
# [-3],
# [-1]])

Source code in qmm/core/stability.py

@cache
def net_feedback(G: nx.DiGraph, level: Optional[int] = None) -> sp.Matrix:
    """Calculate net feedback at a specified level of the system.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Level of feedback to compute (None for all levels)

    Returns:
        sp.Matrix: Net feedback at specified levels

    References:
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, net_feedback
        net_feedback(load_digraph("snowshoe_rp"), level=2)
        # Matrix([[-3]])

        net_feedback(load_digraph("snowshoe_rp"))
        # Matrix([
        # [-1],
        # [-2],
        # [-3],
        # [-1]])
        ```
    """
    return system_feedback(G, level=level, form="signed")

absolute_feedback cached

absolute_feedback(G, level=None, method='combinations')

Calculate absolute feedback at a specified level of the system.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Level of feedback to compute (None for all levels)	None
method	Literal['combinations', 'polynomial']	Method for computing feedback ('combinations' or 'polynomial')	'combinations'

Returns:

Type	Description
Matrix	sp.Matrix: Total number of feedback terms at specified levels

References

• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, absolute_feedback
absolute_feedback(load_digraph("snowshoe_rp"), level=2)
# Matrix([[3]])

absolute_feedback(load_digraph("snowshoe_rp"))
# Matrix([
# [1],
# [2],
# [3],
# [3]])

Source code in qmm/core/stability.py

@cache
def absolute_feedback(
    G: nx.DiGraph,
    level: Optional[int] = None,
    method: Literal["combinations", "polynomial"] = "combinations",
) -> sp.Matrix:
    """Calculate absolute feedback at a specified level of the system.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Level of feedback to compute (None for all levels)
        method: Method for computing feedback ('combinations' or 'polynomial')

    Returns:
        sp.Matrix: Total number of feedback terms at specified levels

    References:
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, absolute_feedback
        absolute_feedback(load_digraph("snowshoe_rp"), level=2)
        # Matrix([[3]])

        absolute_feedback(load_digraph("snowshoe_rp"))
        # Matrix([
        # [1],
        # [2],
        # [3],
        # [3]])
        ```
    """
    A = create_matrix(G, form="signed")
    if level == 0:
        return sp.Matrix([1])
    n = A.shape[0]
    if level is not None and (level < 0 or level > n):
        raise ValueError(f"Level must be between 0 and {n}")
    if method == "combinations":
        A = sp.matrix2numpy(A).astype(int)
        A = np.abs(A)
        if level is None:
            fb = []
            for k in range(n + 1):
                fb_k = sum(perm(A[np.ix_(c, c)], method="bbfg") for c in combinations(range(n), k))
                fb.append(int(fb_k))
        else:
            fb_k = sum(perm(A[np.ix_(c, c)], method="bbfg") for c in combinations(range(n), level))
            fb = [int(fb_k)]
        if any(v > 2**53 for v in fb):
            raise OverflowError("absolute_feedback exceeds float precision (2**53)")
    elif method == "polynomial":
        lam = sp.Symbol("lambda")
        A_abs = sp.Matrix(sp.Abs(A) + lam * sp.eye(n))
        P = sp.per(A_abs)
        if level is None:
            fb = [P.coeff(lam, n - k) for k in range(n + 1)]
        else:
            fb = [P.coeff(lam, n - level)]
    else:
        raise ValueError("method must be either 'combinations' or 'polynomial'")
    return sp.Matrix(fb)

weighted_feedback cached

weighted_feedback(G, level=None)

Calculate ratio of net to total feedback terms at each level of the system.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Level to compute weighted feedback (None for all levels)	None

Returns:

Type	Description
Matrix	sp.Matrix: Weighted feedback metrics for each level

References

• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, weighted_feedback
weighted_feedback(load_digraph("snowshoe_rp"), level=2)
# Matrix([[-1]])

weighted_feedback(load_digraph("snowshoe_rp"))
# Matrix([
# [  -1],
# [  -1],
# [  -1],
# [-1/3]])

Source code in qmm/core/stability.py

@cache
def weighted_feedback(G: nx.DiGraph, level: Optional[int] = None) -> sp.Matrix:
    """Calculate ratio of net to total feedback terms at each level of the system.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Level to compute weighted feedback (None for all levels)

    Returns:
        sp.Matrix: Weighted feedback metrics for each level

    References:
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, weighted_feedback
        weighted_feedback(load_digraph("snowshoe_rp"), level=2)
        # Matrix([[-1]])

        weighted_feedback(load_digraph("snowshoe_rp"))
        # Matrix([
        # [  -1],
        # [  -1],
        # [  -1],
        # [-1/3]])
        ```
    """
    net_fb = net_feedback(G, level=level)
    tot_fb = absolute_feedback(G, level=level)
    return get_weight(net_fb, tot_fb)

feedback_metrics cached

feedback_metrics(G)

Calculate net, absolute and weighted feedback metrics at each level of the system.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Feedback metrics for each system level

References

• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, feedback_metrics
feedback_metrics(load_digraph("snowshoe_rp"))
#   Feedback level Net Absolute Positive Negative Weighted
# 0              0  -1        1        0        1       -1
# 1              1  -2        2        0        2       -1
# 2              2  -3        3        0        3       -1
# 3              3  -1        3        1        2     -1/3

Source code in qmm/core/stability.py

@cache
def feedback_metrics(G: nx.DiGraph) -> pd.DataFrame:
    """Calculate net, absolute and weighted feedback metrics at each level of the system.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        pd.DataFrame: Feedback metrics for each system level

    References:
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, feedback_metrics
        feedback_metrics(load_digraph("snowshoe_rp"))
        #   Feedback level Net Absolute Positive Negative Weighted
        # 0              0  -1        1        0        1       -1
        # 1              1  -2        2        0        2       -1
        # 2              2  -3        3        0        3       -1
        # 3              3  -1        3        1        2     -1/3
        ```
    """
    net = net_feedback(G)
    absolute = absolute_feedback(G)
    positive = get_positive(net, absolute)
    negative = get_negative(net, absolute)
    weighted = weighted_feedback(G)
    n = len(positive)
    levels = [str(i) for i in range(n)]

    df = {
        "Feedback level": levels,
        "Net": [net[i, 0] for i in range(n)],
        "Absolute": [absolute[i, 0] for i in range(n)],
        "Positive": [positive[i, 0] for i in range(n)],
        "Negative": [negative[i, 0] for i in range(n)],
        "Weighted": [weighted[i, 0] for i in range(n)],
    }

    return pd.DataFrame(df)

hurwitz_determinants cached

hurwitz_determinants(G, level=None, form='symbolic')

Calculate Hurwitz determinants for analysing system stability.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Level to compute determinants (None for all Hurwitz determinants)	None
form	Literal['symbolic', 'signed']	Type of computation ('symbolic' or 'signed')	'symbolic'

Returns:

Type	Description
Matrix	sp.Matrix: Hurwitz determinants at specified levels

References

• Hurwitz, A. (1895). On the conditions under which an equation has only roots with negative real parts. Mathematische Annelen 65, 273–284.
• Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

Examples:

from qmm import load_digraph, hurwitz_determinants
hurwitz_determinants(load_digraph("snowshoe_rp"), level=2, form='symbolic')
# Matrix([[a_C,P*a_P,C*a_P,P + a_C,P*a_P,R*a_R,C + a_C,R*a_R,C*a_R,R + a_P,P**2*a_R,R + a_P,P*a_R,R**2]])

Source code in qmm/core/stability.py

@cache
def hurwitz_determinants(
    G: nx.DiGraph,
    level: Optional[int] = None,
    form: Literal["symbolic", "signed"] = "symbolic",
) -> sp.Matrix:
    """Calculate Hurwitz determinants for analysing system stability.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Level to compute determinants (None for all Hurwitz determinants)
        form: Type of computation ('symbolic' or 'signed')

    Returns:
        sp.Matrix: Hurwitz determinants at specified levels

    References:
        - Hurwitz, A. (1895). On the conditions under which an equation has only roots with negative real parts. Mathematische Annelen 65, 273–284.
        - Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

    Examples:
        ```python
        from qmm import load_digraph, hurwitz_determinants
        hurwitz_determinants(load_digraph("snowshoe_rp"), level=2, form='symbolic')
        # Matrix([[a_C,P*a_P,C*a_P,P + a_C,P*a_P,R*a_R,C + a_C,R*a_R,C*a_R,R + a_P,P**2*a_R,R + a_P,P*a_R,R**2]])
        ```
    """
    fb = system_feedback(G, level=None, form=form)
    n = len(fb) - 1
    if n > 5 and form == "symbolic":
        raise ValueError("Limited to systems with five or fewer variables.")
    if level is not None and (level < 0 or level > n):
        raise ValueError(f"Level must be between 0 and {n}")
    if level is None:
        h = _hurwitz_matrix(fb, n)
        hd = sp.Matrix([sp.det(h[:k, :k]) for k in range(0, n + 1)])
    else:
        h = _hurwitz_matrix(fb, level)
        hd = sp.Matrix([sp.det(h[:level, :level])])
    return sp.Matrix(hd)

net_determinants cached

net_determinants(G, level=None)

Calculate net terms in Hurwitz determinants.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Level to compute determinants (None for all Hurwitz determinants)	None

Returns:

Type	Description
Matrix	sp.Matrix: Net terms in Hurwitz determinants

References

• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, net_determinants
net_determinants(load_digraph("snowshoe_rp"), level=2)
# Matrix([[5]])

net_determinants(load_digraph("snowshoe_rp"))
# Matrix([
# [1],
# [2],
# [5],
# [5]])

Source code in qmm/core/stability.py

@cache
def net_determinants(G: nx.DiGraph, level: Optional[int] = None) -> sp.Matrix:
    """Calculate net terms in Hurwitz determinants.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Level to compute determinants (None for all Hurwitz determinants)

    Returns:
        sp.Matrix: Net terms in Hurwitz determinants

    References:
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, net_determinants
        net_determinants(load_digraph("snowshoe_rp"), level=2)
        # Matrix([[5]])

        net_determinants(load_digraph("snowshoe_rp"))
        # Matrix([
        # [1],
        # [2],
        # [5],
        # [5]])
        ```
    """
    return hurwitz_determinants(G, level=level, form="signed")

absolute_determinants cached

absolute_determinants(G, level=None)

Calculate absolute terms in Hurwitz determinants.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Level to compute determinants (None for all Hurwitz determinants)	None

Returns:

Type	Description
Matrix	sp.Matrix: Absolute terms in Hurwitz determinants

References

• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, absolute_determinants
absolute_determinants(load_digraph("snowshoe_rp"), level=2)
# Matrix([[9]])

absolute_determinants(load_digraph("snowshoe_rp"))
# Matrix([
# [ 1],
# [ 2],
# [ 9],
# [27]])

Source code in qmm/core/stability.py

@cache
def absolute_determinants(G: nx.DiGraph, level: Optional[int] = None) -> sp.Matrix:
    """Calculate absolute terms in Hurwitz determinants.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Level to compute determinants (None for all Hurwitz determinants)

    Returns:
        sp.Matrix: Absolute terms in Hurwitz determinants

    References:
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, absolute_determinants
        absolute_determinants(load_digraph("snowshoe_rp"), level=2)
        # Matrix([[9]])

        absolute_determinants(load_digraph("snowshoe_rp"))
        # Matrix([
        # [ 1],
        # [ 2],
        # [ 9],
        # [27]])
        ```
    """
    tot_fb = absolute_feedback(G)
    n = tot_fb.shape[0] - 1
    h = _hurwitz_matrix(tot_fb, n)
    if level is None:
        td = [sp.Integer(1)]
        for k in range(1, n + 1):
            h_k = np.array(h[:k, :k].tolist(), dtype=float)
            td.append(sp.Abs(sp.Integer(int(round(perm(h_k))))))
    else:
        if level < 0 or level > n:
            raise ValueError(f"Level must be between 0 and {n}")
        if level == 0:
            td = [sp.Integer(1)]
        else:
            H_k = np.array(h[:level, :level].tolist(), dtype=float)
            td = [sp.Abs(sp.Integer(int(round(perm(H_k)))))]
    return sp.Matrix(td)

weighted_determinants cached

weighted_determinants(G, level=None)

Calculate ratio of net to total terms for Hurwitz determinants.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Level to compute determinants (None for all Hurwitz determinants)	None

Returns:

Type	Description
Matrix	sp.Matrix: Ratio of net to total terms for Hurwitz determinants

References

• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, weighted_determinants
weighted_determinants(load_digraph("snowshoe_rp"), level=2)
# Matrix([[5/9]])

weighted_determinants(load_digraph("snowshoe_rp"))
# Matrix([
# [   1],
# [   1],
# [ 5/9],
# [5/27]])

Source code in qmm/core/stability.py

@cache
def weighted_determinants(G: nx.DiGraph, level: Optional[int] = None) -> sp.Matrix:
    """Calculate ratio of net to total terms for Hurwitz determinants.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Level to compute determinants (None for all Hurwitz determinants)

    Returns:
        sp.Matrix: Ratio of net to total terms for Hurwitz determinants

    References:
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, weighted_determinants
        weighted_determinants(load_digraph("snowshoe_rp"), level=2)
        # Matrix([[5/9]])

        weighted_determinants(load_digraph("snowshoe_rp"))
        # Matrix([
        # [   1],
        # [   1],
        # [ 5/9],
        # [5/27]])
        ```
    """
    net_det = net_determinants(G, level=level)
    tot_det = absolute_determinants(G, level=level)
    wgt_det = get_weight(net_det, tot_det)
    return wgt_det

determinants_metrics cached

determinants_metrics(G)

Calculate net, absolute and weighted Hurwitz determinant metrics.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Hurwitz determinant metrics

References

• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, determinants_metrics
determinants_metrics(load_digraph("snowshoe_rp"))
#   Hurwitz determinant Net Absolute Weighted
# 0                   0   1        1        1
# 1                   1   2        2        1
# 2                   2   5        9      5/9
# 3                   3   5       27     5/27

Source code in qmm/core/stability.py

@cache
def determinants_metrics(G: nx.DiGraph) -> pd.DataFrame:
    """Calculate net, absolute and weighted Hurwitz determinant metrics.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        pd.DataFrame: Hurwitz determinant metrics

    References:
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, determinants_metrics
        determinants_metrics(load_digraph("snowshoe_rp"))
        #   Hurwitz determinant Net Absolute Weighted
        # 0                   0   1        1        1
        # 1                   1   2        2        1
        # 2                   2   5        9      5/9
        # 3                   3   5       27     5/27
        ```
    """
    net = net_determinants(G)
    absolute = absolute_determinants(G)
    weighted = weighted_determinants(G)
    n = len(net)
    levels = [str(i) for i in range(n)]
    df = {
        "Hurwitz determinant": levels,
        "Net": [net[i, 0] for i in range(n)],
        "Absolute": [absolute[i, 0] for i in range(n)],
        "Weighted": [weighted[i, 0] for i in range(n)],
    }
    return pd.DataFrame(df)

conditional_stability cached

conditional_stability(G)

Analyse conditional stability metrics and model stability class.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Conditional stability metrics and model class

References

• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, conditional_stability
conditional_stability(load_digraph("snowshoe_rp"))
#                       Test                                                 Definition   Result
# 0        Weighted feedback                        Maximum weighted feedback (level 3)    -0.33
# 1     Weighted determinant                          n-1 weighted determinant at level     0.56
# 2  Ratio to model-c system                           Ratio to a 'model-c' type system      1.7
# 3              Model class  Class of the model based on conditional stability metrics  Class I

Source code in qmm/core/stability.py

@cache
def conditional_stability(G: nx.DiGraph) -> pd.DataFrame:
    """Analyse conditional stability metrics and model stability class.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        pd.DataFrame: Conditional stability metrics and model class

    References:
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, conditional_stability
        conditional_stability(load_digraph("snowshoe_rp"))
        #                       Test                                                 Definition   Result
        # 0        Weighted feedback                        Maximum weighted feedback (level 3)    -0.33
        # 1     Weighted determinant                          n-1 weighted determinant at level     0.56
        # 2  Ratio to model-c system                           Ratio to a 'model-c' type system      1.7
        # 3              Model class  Class of the model based on conditional stability metrics  Class I
        ```
    """
    A = create_matrix(G, form="signed")
    n = A.shape[0]
    w_fb = weighted_feedback(G)
    w_det = weighted_determinants(G, level=n - 1)[0]
    C = _create_model_c(n)
    w_det_c = weighted_determinants(C, level=n - 1)[0]
    ratio_C = w_det / w_det_c
    max_fb_n = np.max(w_fb) == w_fb[-1]
    kmax = len(w_fb) - 1 - np.argmax(w_fb[::-1])
    is_sign_stable = sign_stability(G)["Result"].iloc[-1]
    if is_sign_stable:
        model_class = "Sign stable"
    elif max_fb_n and ratio_C >= 1:
        model_class = "Class I"
    else:
        model_class = "Class II"
    stability_metrics = pd.DataFrame(
        {
            "Test": [
                "Weighted feedback",
                "Weighted determinant",
                "Ratio to model-c system",
                "Model class",
            ],
            "Definition": [
                f"Maximum weighted feedback (level {kmax})",
                "n-1 weighted determinant at level",
                "Ratio to a 'model-c' type system",
                "Class of the model based on conditional stability metrics",
            ],
            "Result": [
                np.max(w_fb).evalf(2),
                w_det.evalf(2),
                ratio_C.evalf(2),
                model_class,
            ],
        }
    )
    return stability_metrics

simulation_stability

simulation_stability(G, n_sim=10000, dist='uniform', seed=42, presample=None)

Analyse stability using randomly sampled interaction strengths from a specified distribution.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
n_sim	int	Number of simulations to perform (default 10000)	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution to sample from (default 'uniform'): - "uniform": Uniform(0, 1) - no assumption about interaction strength - "weak": Beta(1, 3) - weak interactions predominate - "moderate": Beta(2, 2) - moderate interactions predominate - "strong": Beta(3, 1) - strong interactions predominate - "uniform_two_oom": Uniform(0.01, 1)	'uniform'
seed	int	Random seed	42
presample	Optional[ndarray]	Optional pre-sampled values of shape (n_sim, n, n) to use instead of random sampling	None

Returns:

Type	Description
DataFrame	pd.DataFrame: Proportion of stable matrices and proportion that fail Hurwitz criteria

References

• May, R.M. (1973). Qualitative Stability in Model Ecosystems. Ecology 54, 638–641.
• Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

Examples:

from qmm import load_digraph, simulation_stability
simulation_stability(load_digraph("snowshoe_rp"), n_sim=1000)
#                         Test                                                             Definition  Result
# 0            Stable matrices              Proportion where all eigenvalues have negative real parts  76.70%
# 1          Unstable matrices      Proportion where one or more eigenvalues have positive real parts  23.30%
# 2        Hurwitz criterion i  Proportion where polynomial coefficients are not all of the same sign  23.30%
# 3       Hurwitz criterion ii             Proportion where Hurwitz determinants are not all positive   0.00%
# 4   Hurwitz criterion i only                        Proportion where only Hurwitz criterion i fails  23.30%
# 5  Hurwitz criterion ii only                       Proportion where only Hurwitz criterion ii fails   0.00%

Source code in qmm/core/stability.py

def simulation_stability(
    G: nx.DiGraph,
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
    presample: Optional[np.ndarray] = None,
) -> pd.DataFrame:
    """Analyse stability using randomly sampled interaction strengths from a specified distribution.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        n_sim: Number of simulations to perform (default 10000)
        dist: Distribution to sample from (default 'uniform'):
            - "uniform": Uniform(0, 1) - no assumption about interaction strength
            - "weak": Beta(1, 3) - weak interactions predominate
            - "moderate": Beta(2, 2) - moderate interactions predominate
            - "strong": Beta(3, 1) - strong interactions predominate
            - "uniform_two_oom": Uniform(0.01, 1)
        seed: Random seed
        presample: Optional pre-sampled values of shape (n_sim, n, n) to use instead of random sampling

    Returns:
        pd.DataFrame: Proportion of stable matrices and proportion that fail Hurwitz criteria

    References:
        - May, R.M. (1973). Qualitative Stability in Model Ecosystems. Ecology 54, 638–641.
        - Dambacher, J.M., Luh, H.-K., Li, H.W., Rossignol, P.A. (2003). Qualitative stability and ambiguity in model ecosystems. The American Naturalist 161, 876–888.

    Examples:
        ```python
        from qmm import load_digraph, simulation_stability
        simulation_stability(load_digraph("snowshoe_rp"), n_sim=1000)
        #                         Test                                                             Definition  Result
        # 0            Stable matrices              Proportion where all eigenvalues have negative real parts  76.70%
        # 1          Unstable matrices      Proportion where one or more eigenvalues have positive real parts  23.30%
        # 2        Hurwitz criterion i  Proportion where polynomial coefficients are not all of the same sign  23.30%
        # 3       Hurwitz criterion ii             Proportion where Hurwitz determinants are not all positive   0.00%
        # 4   Hurwitz criterion i only                        Proportion where only Hurwitz criterion i fails  23.30%
        # 5  Hurwitz criterion ii only                       Proportion where only Hurwitz criterion ii fails   0.00%
        ```
    """
    A = create_matrix(G, "signed")
    A = sp.matrix2numpy(A).astype(int)

    if presample is not None:
        if presample.shape != (n_sim, *A.shape):
            raise ValueError(f"presample must have shape ({n_sim}, {A.shape[0]}, {A.shape[1]})")

    rng = np.random.RandomState(seed)
    n_stable = 0
    n_unstable = 0
    n_hurwitz_i_fail = 0
    n_hurwitz_ii_fail = 0
    n_hurwitz_i_only_fail = 0
    n_hurwitz_ii_only_fail = 0

    for i in range(n_sim):
        if presample is not None:
            M = presample[i]
        else:
            M = _random_sampler(dist, A.size, rng).reshape(A.shape)
        S = A * M
        if np.all(np.real(np.linalg.eigvals(S)) < 0):
            n_stable += 1
        else:
            n_unstable += 1
        pc = np.poly(S)
        hurwitz_i = np.all(pc[1:] > 0) or np.all(pc[1:] < 0)
        n = len(pc)
        H = np.zeros((n - 1, n - 1))
        for r in range(1, n):
            for c in range(1, n):
                index = 2 * c - r
                if 0 <= index < n:
                    H[r - 1, c - 1] = pc[index]
        hd = [np.linalg.det(H[: k + 1, : k + 1]) for k in range(n - 1)]
        hurwitz_ii = np.all(np.array(hd[1:-1]) > 0)
        if not hurwitz_i:
            n_hurwitz_i_fail += 1
            if hurwitz_ii:
                n_hurwitz_i_only_fail += 1
        if not hurwitz_ii:
            n_hurwitz_ii_fail += 1
            if hurwitz_i:
                n_hurwitz_ii_only_fail += 1

    prop_stable = n_stable / n_sim
    prop_unstable = n_unstable / n_sim
    prop_hurwitz_i_fail = n_hurwitz_i_fail / n_sim
    prop_hurwitz_ii_fail = n_hurwitz_ii_fail / n_sim
    prop_hurwitz_i_only_fail = n_hurwitz_i_only_fail / n_sim
    prop_hurwitz_ii_only_fail = n_hurwitz_ii_only_fail / n_sim

    sim_df = pd.DataFrame(
        {
            "Test": [
                "Stable matrices",
                "Unstable matrices",
                "Hurwitz criterion i",
                "Hurwitz criterion ii",
                "Hurwitz criterion i only",
                "Hurwitz criterion ii only",
            ],
            "Definition": [
                "Proportion where all eigenvalues have negative real parts",
                "Proportion where one or more eigenvalues have positive real parts",
                "Proportion where polynomial coefficients are not " + "all of the same sign",
                "Proportion where Hurwitz determinants are not all positive",
                "Proportion where only Hurwitz criterion i fails",
                "Proportion where only Hurwitz criterion ii fails",
            ],
            "Result": [
                f"{prop_stable:.2%}",
                f"{prop_unstable:.2%}",
                f"{prop_hurwitz_i_fail:.2%}",
                f"{prop_hurwitz_ii_fail:.2%}",
                f"{prop_hurwitz_i_only_fail:.2%}",
                f"{prop_hurwitz_ii_only_fail:.2%}",
            ],
        }
    )
    return sim_df

qmm.core.press

Analyse direct and indirect effects of press perturbations.

adjoint_matrix cached

adjoint_matrix(G, form='symbolic', perturb=None)

Calculate elements of classical adjoint matrix for press perturbation response.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
form	Literal['symbolic', 'signed']	Type of computation ('symbolic', 'signed')	'symbolic'
perturb	Optional[str]	Node to perturb (None for full matrix)	None

Returns:

Type	Description
Matrix	sp.Matrix: Classical adjoint matrix elements

References

• Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.

Examples:

from qmm import load_digraph, adjoint_matrix
adjoint_matrix(load_digraph("snowshoe_rp"), form='symbolic')
# Matrix([
# [               a_C,P*a_P,C,              -a_P,P*a_R,C,  a_C,P*a_R,C],
# [-a_C,P*a_P,R + a_C,R*a_P,P,               a_P,P*a_R,R, -a_C,P*a_R,R],
# [               a_C,R*a_P,C, a_P,C*a_R,R - a_P,R*a_R,C,  a_C,R*a_R,C]])

adjoint_matrix(load_digraph("snowshoe_rp"), form='symbolic', perturb='R')
# Matrix([
# [               a_C,P*a_P,C],
# [-a_C,P*a_P,R + a_C,R*a_P,P],
# [               a_C,R*a_P,C]])

adjoint_matrix(load_digraph("snowshoe_rp"), form='signed')
# Matrix([
# [1, -1,  1],
# [0,  1, -1],
# [1,  0,  1]])

adjoint_matrix(load_digraph("snowshoe_rp"), form='signed', perturb='R')
# Matrix([
# [1],
# [0],
# [1]])

Source code in qmm/core/press.py

@cache
def adjoint_matrix(
    G: nx.DiGraph,
    form: Literal["symbolic", "signed"] = "symbolic",
    perturb: Optional[str] = None,
) -> sp.Matrix:
    """Calculate elements of classical adjoint matrix for press perturbation response.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        form: Type of computation ('symbolic', 'signed')
        perturb: Node to perturb (None for full matrix)

    Returns:
        sp.Matrix: Classical adjoint matrix elements

    References:
        - Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.

    Examples:
        ```python
        from qmm import load_digraph, adjoint_matrix
        adjoint_matrix(load_digraph("snowshoe_rp"), form='symbolic')
        # Matrix([
        # [               a_C,P*a_P,C,              -a_P,P*a_R,C,  a_C,P*a_R,C],
        # [-a_C,P*a_P,R + a_C,R*a_P,P,               a_P,P*a_R,R, -a_C,P*a_R,R],
        # [               a_C,R*a_P,C, a_P,C*a_R,R - a_P,R*a_R,C,  a_C,R*a_R,C]])

        adjoint_matrix(load_digraph("snowshoe_rp"), form='symbolic', perturb='R')
        # Matrix([
        # [               a_C,P*a_P,C],
        # [-a_C,P*a_P,R + a_C,R*a_P,P],
        # [               a_C,R*a_P,C]])

        adjoint_matrix(load_digraph("snowshoe_rp"), form='signed')
        # Matrix([
        # [1, -1,  1],
        # [0,  1, -1],
        # [1,  0,  1]])

        adjoint_matrix(load_digraph("snowshoe_rp"), form='signed', perturb='R')
        # Matrix([
        # [1],
        # [0],
        # [1]])
        ```
    """
    A = create_matrix(G, form=form)
    A = sp.Matrix(-A)
    nodes = get_nodes(G, "state")
    if perturb is not None and perturb not in nodes:
        raise ValueError(f"Perturbation node must be one of: {nodes}")
    n = len(nodes)
    if perturb is not None:
        src_id = nodes.index(perturb)
        return sp.Matrix([sp.Integer(-1) ** (src_id + j) * A.minor(src_id, j) for j in range(n)])
    adjoint_matrix = sp.expand(A.adjugate())
    return sp.Matrix(adjoint_matrix)

absolute_feedback_matrix cached

absolute_feedback_matrix(G, perturb=None)

Calculate total number of both positive and negative terms for press perturbation response.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
perturb	Optional[str]	Node to perturb (None for full matrix)	None

Returns:

Type	Description
Matrix	sp.Matrix: Absolute feedback matrix elements

References

• Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.

Examples:

from qmm import load_digraph, absolute_feedback_matrix
absolute_feedback_matrix(load_digraph("snowshoe_rp"))
# Matrix([
# [1, 1, 1],
# [2, 1, 1],
# [1, 2, 1]])

absolute_feedback_matrix(load_digraph("snowshoe_rp"), perturb='R')
# Matrix([
# [1],
# [2],
# [1]])

Source code in qmm/core/press.py

@cache
def absolute_feedback_matrix(G: nx.DiGraph, perturb: Optional[str] = None) -> sp.Matrix:
    """Calculate total number of both positive and negative terms for press perturbation response.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        perturb: Node to perturb (None for full matrix)

    Returns:
        sp.Matrix: Absolute feedback matrix elements

    References:
        - Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.

    Examples:
        ```python
        from qmm import load_digraph, absolute_feedback_matrix
        absolute_feedback_matrix(load_digraph("snowshoe_rp"))
        # Matrix([
        # [1, 1, 1],
        # [2, 1, 1],
        # [1, 2, 1]])

        absolute_feedback_matrix(load_digraph("snowshoe_rp"), perturb='R')
        # Matrix([
        # [1],
        # [2],
        # [1]])
        ```
    """
    A = create_matrix(G, form="binary")
    A_np = np.array(sp.matrix2numpy(A), dtype=float)
    nodes = get_nodes(G, "state")
    if perturb is not None and perturb not in nodes:
        raise ValueError(f"Perturbation node must be one of: {nodes}")
    n = A_np.shape[0]
    if perturb is not None:
        perturb_index = nodes.index(perturb)
        result = np.zeros(n, dtype=int)
        for j in range(n):
            minor = np.delete(np.delete(A_np, perturb_index, 0), j, 1)
            result[j] = int(perm(minor))
        return sp.Matrix(result)
    tmat = np.zeros((n, n), dtype=int)
    for i in range(n):
        for j in range(n):
            minor = np.delete(np.delete(A_np, j, 0), i, 1)
            tmat[i, j] = int(perm(minor))
    return sp.Matrix(tmat)

weighted_predictions_matrix cached

weighted_predictions_matrix(G, as_nan=True, as_abs=False, perturb=None)

Calculate ratio of net to total terms for a press perturbation response.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
as_nan	bool	Return NaN for undefined ratios	True
as_abs	bool	Return absolute values	False
perturb	Optional[str]	Node to perturb (None for full matrix)	None

Returns:

Type	Description
Matrix	sp.Matrix: Prediction weights

References

• Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.

Examples:

from qmm import load_digraph, weighted_predictions_matrix
weighted_predictions_matrix(load_digraph("snowshoe_rp"))
# Matrix([
# [1, -1,  1],
# [0,  1, -1],
# [1,  0,  1]])

weighted_predictions_matrix(load_digraph("snowshoe_rp"), perturb='R')
# Matrix([
# [1],
# [0],
# [1]])

weighted_predictions_matrix(load_digraph("snowshoe_rp"), as_abs=True)
# Matrix([
# [1, 1, 1],
# [0, 1, 1],
# [1, 0, 1]])

weighted_predictions_matrix(load_digraph("snowshoe_rp"), as_abs=False)
# Matrix([
# [1, -1,  1],
# [0,  1, -1],
# [1,  0,  1]])

Source code in qmm/core/press.py

@cache
def weighted_predictions_matrix(G: nx.DiGraph, as_nan: bool = True, as_abs: bool = False, perturb: Optional[str] = None) -> sp.Matrix:
    """Calculate ratio of net to total terms for a press perturbation response.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        as_nan: Return NaN for undefined ratios
        as_abs: Return absolute values
        perturb: Node to perturb (None for full matrix)

    Returns:
        sp.Matrix: Prediction weights

    References:
        - Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.

    Examples:
        ```python
        from qmm import load_digraph, weighted_predictions_matrix
        weighted_predictions_matrix(load_digraph("snowshoe_rp"))
        # Matrix([
        # [1, -1,  1],
        # [0,  1, -1],
        # [1,  0,  1]])

        weighted_predictions_matrix(load_digraph("snowshoe_rp"), perturb='R')
        # Matrix([
        # [1],
        # [0],
        # [1]])

        weighted_predictions_matrix(load_digraph("snowshoe_rp"), as_abs=True)
        # Matrix([
        # [1, 1, 1],
        # [0, 1, 1],
        # [1, 0, 1]])

        weighted_predictions_matrix(load_digraph("snowshoe_rp"), as_abs=False)
        # Matrix([
        # [1, -1,  1],
        # [0,  1, -1],
        # [1,  0,  1]])
        ```
    """
    amat = adjoint_matrix(G, perturb=perturb, form="signed")
    if as_abs:
        amat = sp.Abs(amat)
    tmat = absolute_feedback_matrix(G, perturb=perturb)
    if as_nan:
        wmat = get_weight(amat, tmat)
    else:
        wmat = get_weight(amat, tmat, sp.Integer(1))
    return sp.Matrix(wmat)

sign_determinacy_matrix cached

sign_determinacy_matrix(
    G, method="average", as_nan=True, as_abs=False, perturb=None
)

Calculate probability of a correct sign prediction (matches adjoint).

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
method	Literal['average', '95_bound']	Method for computing determinacy ('average', '95_bound')	'average'
as_nan	bool	Return NaN for undefined ratios	True
as_abs	bool	Return absolute values	False
perturb	Optional[str]	Node to perturb (None for full matrix)	None

Returns:

Type	Description
Matrix	sp.Matrix: Probability of sign determinacy

References

• Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.

Examples:

from qmm import load_digraph, sign_determinacy_matrix
sign_determinacy_matrix(load_digraph("snowshoe_rp"), method='average')
# Matrix([
# [  1,  -1,  1],
# [1/2,   1, -1],
# [  1, 1/2,  1]])

sign_determinacy_matrix(load_digraph("snowshoe_rp"), method='average', perturb='R')
# Matrix([
# [  1],
# [1/2],
# [  1]])

sign_determinacy_matrix(load_digraph("snowshoe_rp"), method='average', as_abs=True)
# Matrix([
# [  1,   1, 1],
# [1/2,   1, 1],
# [  1, 1/2, 1]])

sign_determinacy_matrix(load_digraph("snowshoe_rp"), method='average', as_abs=False)
# Matrix([
# [  1,  -1,  1],
# [1/2,   1, -1],
# [  1, 1/2,  1]])

Source code in qmm/core/press.py

@cache
def sign_determinacy_matrix(
    G: nx.DiGraph,
    method: Literal["average", "95_bound"] = "average",
    as_nan: bool = True,
    as_abs: bool = False,
    perturb: Optional[str] = None,
) -> sp.Matrix:
    """Calculate probability of a correct sign prediction (matches adjoint).

    Args:
        G: NetworkX DiGraph representing signed digraph model
        method: Method for computing determinacy ('average', '95_bound')
        as_nan: Return NaN for undefined ratios
        as_abs: Return absolute values
        perturb: Node to perturb (None for full matrix)

    Returns:
        sp.Matrix: Probability of sign determinacy

    References:
        - Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.

    Examples:
        ```python
        from qmm import load_digraph, sign_determinacy_matrix
        sign_determinacy_matrix(load_digraph("snowshoe_rp"), method='average')
        # Matrix([
        # [  1,  -1,  1],
        # [1/2,   1, -1],
        # [  1, 1/2,  1]])

        sign_determinacy_matrix(load_digraph("snowshoe_rp"), method='average', perturb='R')
        # Matrix([
        # [  1],
        # [1/2],
        # [  1]])

        sign_determinacy_matrix(load_digraph("snowshoe_rp"), method='average', as_abs=True)
        # Matrix([
        # [  1,   1, 1],
        # [1/2,   1, 1],
        # [  1, 1/2, 1]])

        sign_determinacy_matrix(load_digraph("snowshoe_rp"), method='average', as_abs=False)
        # Matrix([
        # [  1,  -1,  1],
        # [1/2,   1, -1],
        # [  1, 1/2,  1]])
        ```
    """
    wmat = weighted_predictions_matrix(G, perturb=perturb, as_nan=as_nan, as_abs=as_abs)
    tmat = sp.Matrix(absolute_feedback_matrix(G, perturb=perturb))
    pmat = sign_determinacy(wmat, tmat, method)
    return sp.Matrix(pmat)

numerical_simulations cached

numerical_simulations(
    G,
    n_sim=10000,
    dist="uniform",
    seed=42,
    as_nan=True,
    as_abs=False,
    positive_only=False,
    match_adjoint=False,
)

Calculate proportion of positive and negative responses from stable simulations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
n_sim	int	Number of simulations	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution for sampling ('uniform', 'weak', 'moderate', 'strong')	'uniform'
seed	int	Random seed	42
as_nan	bool	Return NaN for undefined ratios	True
positive_only	bool	Return just the proportion of positive responses instead of sign-dominant proportions.	False
as_abs	bool	Return absolute values	False
match_adjoint	bool	Return proportion of simulations matching the sign of the adjoint. Values are always between 0 and 1. Entries where the adjoint sign is ambiguous (0) return 0.5. Incompatible with positive_only and as_abs.	False

References

• Dambacher, J.M., Li, H.W., Rossignol, P.A. (2003). Qualitative predictions in model ecosystems. Ecological Modelling 161, 79–93.

Returns:

Type	Description
Matrix	sp.Matrix: Average proportion of positive and negative responses

Raises:

Type	Description
ValueError	If invalid parameter combinations are used.

Examples:

from qmm import load_digraph, numerical_simulations
numerical_simulations(load_digraph("snowshoe_rp"), n_sim=1000, seed=42)
# Matrix([
# [  1.0,  -1.0,  1.0],
# [0.621,   1.0, -1.0],
# [  1.0, 0.639,  1.0]])

numerical_simulations(load_digraph("snowshoe_rp"), n_sim=1000, seed=42, as_abs=True)
# Matrix([
# [  1.0,   1.0, 1.0],
# [0.621,   1.0, 1.0],
# [  1.0, 0.639, 1.0]])

numerical_simulations(load_digraph("snowshoe_rp"), n_sim=1000, seed=42, as_abs=False)
# Matrix([
# [  1.0,  -1.0,  1.0],
# [0.621,   1.0, -1.0],
# [  1.0, 0.639,  1.0]])

numerical_simulations(load_digraph("snowshoe_rp"), n_sim=1000, seed=42, positive_only=True)
# Matrix([
# [  1.0,   0.0, 1.0],
# [0.621,   1.0, 0.0],
# [  1.0, 0.639, 1.0]])

Source code in qmm/core/press.py

@cache
def numerical_simulations(
    G: nx.DiGraph,
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
    as_nan: bool = True,
    as_abs: bool = False,
    positive_only: bool = False,
    match_adjoint: bool = False,
) -> sp.Matrix:
    """Calculate proportion of positive and negative responses from stable simulations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        n_sim: Number of simulations
        dist: Distribution for sampling ('uniform', 'weak', 'moderate', 'strong')
        seed: Random seed
        as_nan: Return NaN for undefined ratios
        positive_only: Return just the proportion of positive responses instead of sign-dominant proportions.
        as_abs: Return absolute values
        match_adjoint: Return proportion of simulations matching the sign of the adjoint.
            Values are always between 0 and 1. Entries where the adjoint sign is
            ambiguous (0) return 0.5. Incompatible with positive_only and as_abs.

    References:
        - Dambacher, J.M., Li, H.W., Rossignol, P.A. (2003). Qualitative predictions in model ecosystems. Ecological Modelling 161, 79–93.

    Returns:
        sp.Matrix: Average proportion of positive and negative responses

    Raises:
        ValueError: If invalid parameter combinations are used.

    Examples:
        ```python
        from qmm import load_digraph, numerical_simulations
        numerical_simulations(load_digraph("snowshoe_rp"), n_sim=1000, seed=42)
        # Matrix([
        # [  1.0,  -1.0,  1.0],
        # [0.621,   1.0, -1.0],
        # [  1.0, 0.639,  1.0]])

        numerical_simulations(load_digraph("snowshoe_rp"), n_sim=1000, seed=42, as_abs=True)
        # Matrix([
        # [  1.0,   1.0, 1.0],
        # [0.621,   1.0, 1.0],
        # [  1.0, 0.639, 1.0]])

        numerical_simulations(load_digraph("snowshoe_rp"), n_sim=1000, seed=42, as_abs=False)
        # Matrix([
        # [  1.0,  -1.0,  1.0],
        # [0.621,   1.0, -1.0],
        # [  1.0, 0.639,  1.0]])

        numerical_simulations(load_digraph("snowshoe_rp"), n_sim=1000, seed=42, positive_only=True)
        # Matrix([
        # [  1.0,   0.0, 1.0],
        # [0.621,   1.0, 0.0],
        # [  1.0, 0.639, 1.0]])
        ```
    """
    if positive_only and not as_nan:
        raise ValueError("Invalid parameter combination: positive_only=True requires as_nan=True")
    if as_abs and not as_nan:
        raise ValueError("Invalid parameter combination: as_abs=True requires as_nan=True")
    if match_adjoint and positive_only:
        raise ValueError("Invalid parameter combination: match_adjoint=True is incompatible with positive_only=True")
    if match_adjoint and as_abs:
        raise ValueError("Invalid parameter combination: match_adjoint=True is incompatible with as_abs=True")

    rng = np.random.RandomState(seed)
    A = create_matrix(G, form="symbolic", matrix_type="A")
    state_nodes = get_nodes(G, "state")
    n = len(state_nodes)
    symbols = sorted(list(A.free_symbols), key=str)
    A_sp = sp.lambdify(symbols, A)
    positive = np.zeros((n, n), dtype=int)
    negative = np.zeros((n, n), dtype=int)
    total_simulations = 0
    if match_adjoint:
        adjoint_signs_np = np.array(
            adjoint_matrix(G, form="signed").tolist(), dtype=float
        )
    attempts, max_attempts = 0, n_sim * 100
    while total_simulations < n_sim and attempts < max_attempts:
        attempts += 1
        values = _random_sampler(dist, len(symbols), rng)
        sim_A = A_sp(*values)
        if np.all(np.real(np.linalg.eigvals(sim_A)) < 0):
            try:
                inv_A = np.linalg.inv(-sim_A)
                positive += inv_A > 0
                negative += inv_A < 0
                total_simulations += 1
            except np.linalg.LinAlgError:
                continue
    if total_simulations < n_sim and n_sim > 0:
        raise RuntimeError(f"Maximum iterations reached. Stable proportion: {total_simulations / max_attempts:.4f}")
    if total_simulations == 0:
        smat = np.full((n, n), np.nan)
    elif positive_only:
        smat = positive / total_simulations
    elif match_adjoint:
        matches = np.where(adjoint_signs_np > 0, positive,
                  np.where(adjoint_signs_np < 0, negative, 0))
        smat = np.where(adjoint_signs_np == 0, 0.5,
                        matches / total_simulations)
    else:
        smat = np.where(negative > positive, -negative / total_simulations, positive / total_simulations)
    smat = sp.Matrix(smat.tolist())
    if total_simulations > 0:
        tmat = absolute_feedback_matrix(G)
        tmat_np = np.array(tmat.tolist(), dtype=bool)
        smat = sp.Matrix([[sp.nan if not tmat_np[i, j] else smat[i, j] for j in range(n)] for i in range(n)])
        if as_abs:
            smat = sp.Matrix([[sp.Abs(x) if x != sp.nan else sp.nan for x in row] for row in smat.tolist()])

    if not as_nan:
        fill = sp.Rational(1, 2) if match_adjoint else 0
        smat = sp.Matrix([[fill if sp.nan == x else x for x in row] for row in smat.tolist()])
    return sp.Matrix(smat)

qmm.core.prediction

Generate qualitative predictions of system response to press perturbations with thresholds for ambiguity.

qualitative_predictions

qualitative_predictions(G, generator=numerical_simulations, t1=0.8, t2=0.95)

Create a sympy matrix of qualitative predictions with thresholds for ambiguity.

Parameters:

Name	Type	Description	Default
G	DiGraph	Graph input for the matrix generator	required
generator	Callable	Matrix generator function from core.press module	numerical_simulations
t1	float	Lower threshold for predictions	0.8
t2	float	Higher threshold for predictions	0.95

Returns:

Type	Description
Matrix	sp.Matrix: Qualitative predictions

Raises:

Type	Description
ValueError	If generator is not callable or thresholds are invalid

References

• Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.
• Dambacher, J.M., Li, H.W., Rossignol, P.A. (2003). Qualitative predictions in model ecosystems. Ecological Modelling 161, 79–93.
• Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.

Examples:

from qmm import load_digraph, qualitative_predictions, weighted_predictions_matrix
qualitative_predictions(load_digraph("snowshoe_rp"), generator=weighted_predictions_matrix, t1=0.5, t2=1.0)
# Matrix([
# [+, −, +],
# [?, +, −],
# [+, ?, +]])

Source code in qmm/core/prediction.py

def qualitative_predictions(
    G: nx.DiGraph,
    generator: Union[
        Callable[..., Union[sp.Matrix, np.ndarray, pd.DataFrame]],
        Literal[
            "weighted_predictions_matrix",
            "sign_determinacy_matrix",
            "numerical_simulations",
        ],
    ] = numerical_simulations,
    t1: float = 0.8,
    t2: float = 0.95,
) -> sp.Matrix:
    """Create a sympy matrix of qualitative predictions with thresholds for ambiguity.

    Args:
        G (nx.DiGraph): Graph input for the matrix generator
        generator (Callable): Matrix generator function from core.press module
        t1 (float): Lower threshold for predictions
        t2 (float): Higher threshold for predictions

    Returns:
        sp.Matrix: Qualitative predictions

    Raises:
        ValueError: If generator is not callable or thresholds are invalid

    References:
        - Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.
        - Dambacher, J.M., Li, H.W., Rossignol, P.A. (2003). Qualitative predictions in model ecosystems. Ecological Modelling 161, 79–93.
        - Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.

    Examples:
        ```python
        from qmm import load_digraph, qualitative_predictions, weighted_predictions_matrix
        qualitative_predictions(load_digraph("snowshoe_rp"), generator=weighted_predictions_matrix, t1=0.5, t2=1.0)
        # Matrix([
        # [+, −, +],
        # [?, +, −],
        # [+, ?, +]])
        ```
    """
    if isinstance(generator, str):
        generator = {
            "weighted_predictions_matrix": weighted_predictions_matrix,
            "sign_determinacy_matrix": sign_determinacy_matrix,
            "numerical_simulations": numerical_simulations,
        }.get(generator)
    if not callable(generator):
        raise ValueError(f"Generator must be callable, got: {type(generator)}")

    matrix = generator(G)
    predictions = _apply_thresholds(matrix, t1, t2)
    rows = [[sp.Integer(0) if vals == "0" else sp.Symbol(vals) for vals in row] for row in predictions]
    return sp.Matrix(rows)

table_of_predictions

table_of_predictions(M, t1=0.8, t2=0.95, index=None, columns=None)

Create a table of qualitative predictions with thresholds for ambiguity.

Parameters:

Name	Type	Description	Default
M	Union[Matrix, ndarray]	Matrix of predictions from press perturbation analysis	required
t1	float	Lower threshold for predictions	0.8
t2	float	Higher threshold for predictions	0.95
index	Optional[List[str]]	Row labels	None
columns	Optional[List[str]]	Column labels	None

Returns:

Type	Description
DataFrame	pd.DataFrame: Qualitative predictions

References

• Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.
• Dambacher, J.M., Li, H.W., Rossignol, P.A. (2003). Qualitative predictions in model ecosystems. Ecological Modelling 161, 79–93.
• Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.

Examples:

from qmm import load_digraph, weighted_predictions_matrix, table_of_predictions
W = weighted_predictions_matrix(load_digraph("snowshoe_rp"))
nodes = ['R', 'C', 'P']
table_of_predictions(W, t1=0.5, t2=1.0, index=nodes, columns=nodes)
#    R  C  P
# R  +  −  +
# C  ?  +  −
# P  +  ?  +

Source code in qmm/core/prediction.py

def table_of_predictions(
    M: Union[sp.Matrix, np.ndarray],
    t1: float = 0.8,
    t2: float = 0.95,
    index: Optional[List[str]] = None,
    columns: Optional[List[str]] = None,
) -> pd.DataFrame:
    """Create a table of qualitative predictions with thresholds for ambiguity.

    Args:
        M (Union[sp.Matrix, np.ndarray]): Matrix of predictions from press perturbation analysis
        t1 (float): Lower threshold for predictions
        t2 (float): Higher threshold for predictions
        index (Optional[List[str]]): Row labels
        columns (Optional[List[str]]): Column labels

    Returns:
        pd.DataFrame: Qualitative predictions

    References:
        - Dambacher, J.M., Li, H.W., Rossignol, P.A. (2002). Relevance of Community Structure in Assessing Indeterminacy of Ecological Predictions. Ecology 83, 1372–1385.
        - Dambacher, J.M., Li, H.W., Rossignol, P.A. (2003). Qualitative predictions in model ecosystems. Ecological Modelling 161, 79–93.
        - Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.

    Examples:
        ```python
        from qmm import load_digraph, weighted_predictions_matrix, table_of_predictions
        W = weighted_predictions_matrix(load_digraph("snowshoe_rp"))
        nodes = ['R', 'C', 'P']
        table_of_predictions(W, t1=0.5, t2=1.0, index=nodes, columns=nodes)
        #    R  C  P
        # R  +  −  +
        # C  ?  +  −
        # P  +  ?  +
        ```
    """
    if isinstance(M, sp.Matrix) and any(isinstance(v, (sp.Symbol, str)) for v in M):
        return pd.DataFrame(M.tolist(), index=index, columns=columns)
    predictions = _apply_thresholds(M, t1, t2)
    return pd.DataFrame(predictions, index=index, columns=columns)

compare_predictions

compare_predictions(M1, M2)

Compare predictions between alternative models or prediction methods.

Parameters:

Name	Type	Description	Default
M1	DataFrame	First matrix of predictions	required
M2	DataFrame	Second matrix of predictions	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Combined predictions showing differences and agreements

Examples:

from qmm import load_digraph, weighted_predictions_matrix, table_of_predictions, compare_predictions
G1 = load_digraph("snowshoe_rp")
G2 = G1.copy()
G2.remove_edge('C', 'P')
nodes = ['R', 'C', 'P']
W1 = weighted_predictions_matrix(G1)
W2 = weighted_predictions_matrix(G2)
pred1 = table_of_predictions(W1, t1=0.5, t2=1.0, index=nodes, columns=nodes)
pred2 = table_of_predictions(W2, t1=0.5, t2=1.0, index=nodes, columns=nodes)
compare_predictions(pred1, pred2)
#       R     C  P
# R  +, 0     −  +
# C     ?     +  −
# P  +, 0  ?, −  +

Source code in qmm/core/prediction.py

def compare_predictions(M1: pd.DataFrame, M2: pd.DataFrame) -> pd.DataFrame:
    """Compare predictions between alternative models or prediction methods.

    Args:
        M1 (pd.DataFrame): First matrix of predictions
        M2 (pd.DataFrame): Second matrix of predictions

    Returns:
        pd.DataFrame: Combined predictions showing differences and agreements

    Examples:
        ```python
        from qmm import load_digraph, weighted_predictions_matrix, table_of_predictions, compare_predictions
        G1 = load_digraph("snowshoe_rp")
        G2 = G1.copy()
        G2.remove_edge('C', 'P')
        nodes = ['R', 'C', 'P']
        W1 = weighted_predictions_matrix(G1)
        W2 = weighted_predictions_matrix(G2)
        pred1 = table_of_predictions(W1, t1=0.5, t2=1.0, index=nodes, columns=nodes)
        pred2 = table_of_predictions(W2, t1=0.5, t2=1.0, index=nodes, columns=nodes)
        compare_predictions(pred1, pred2)
        #       R     C  P
        # R  +, 0     −  +
        # C     ?     +  −
        # P  +, 0  ?, −  +
        ```
    """
    if not M1.index.equals(M2.index) or not M1.columns.equals(M2.columns):
        raise ValueError("M1 and M2 must have the same index and columns")
    M1_str = M1.astype(str)
    M2_str = M2.astype(str)
    combined = pd.DataFrame(
        index=M1.index,
        columns=M1.columns,
        data=np.where(M1_str.values == M2_str.values, M1_str.values, M1_str.values + ", " + M2_str.values),
    )
    return combined

qmm.core.helper

Utility functions for model development and analysis.

list_to_digraph

list_to_digraph(matrix, ids=None)

Convert an adjacency matrix to a directed graph.

Parameters:

Name	Type	Description	Default
matrix	Union[List[List[int]], ndarray]	A square matrix (list of lists or numpy array) representing the adjacency matrix. Non-zero values indicate edges, where the value represents the sign of the edge.	required
ids	Optional[List[str]]	Optional list of node identifiers. If None, nodes will be labeled 1 to n.	None

Returns:

Type	Description
DiGraph	nx.DiGraph: A NetworkX directed graph with signed edges.

Examples:

from qmm import list_to_digraph
G = list_to_digraph([[-1, -1, 0], [1, 0, -1], [1, 1, -1]])
list(G.nodes())
# ['1', '2', '3']

list(G.edges(data='sign'))
# [('1', '1', -1), ('1', '2', 1), ('1', '3', 1), ('2', '1', -1), ('2', '3', 1), ('3', '2', -1), ('3', '3', -1)]

Source code in qmm/core/helper.py

def list_to_digraph(matrix: Union[List[List[int]], np.ndarray], ids: Optional[List[str]] = None) -> nx.DiGraph:
    """Convert an adjacency matrix to a directed graph.

    Args:
        matrix: A square matrix (list of lists or numpy array) representing the adjacency matrix.
            Non-zero values indicate edges, where the value represents the sign of the edge.
        ids: Optional list of node identifiers. If None, nodes will be labeled 1 to n.

    Returns:
        nx.DiGraph: A NetworkX directed graph with signed edges.

    Examples:
        ```python
        from qmm import list_to_digraph
        G = list_to_digraph([[-1, -1, 0], [1, 0, -1], [1, 1, -1]])
        list(G.nodes())
        # ['1', '2', '3']

        list(G.edges(data='sign'))
        # [('1', '1', -1), ('1', '2', 1), ('1', '3', 1), ('2', '1', -1), ('2', '3', 1), ('3', '2', -1), ('3', '3', -1)]
        ```
    """
    if not isinstance(matrix, (list, np.ndarray)):
        raise ValueError("Input must be a list of lists or a numpy array")
    if isinstance(matrix, list):
        matrix = np.array(matrix)
    if matrix.ndim != 2 or matrix.shape[0] != matrix.shape[1]:
        raise ValueError("Input must be a square matrix")
    G = nx.DiGraph()
    n = matrix.shape[0]
    if ids is None:
        node_ids = [str(i) for i in range(1, n + 1)]
    else:
        if len(ids) != n:
            raise ValueError("Number of ids must match matrix dimensions")
        node_ids = ids
    G.add_nodes_from(node_ids)
    for i in range(n):
        for j in range(n):
            if matrix[i][j] != 0:
                G.add_edge(node_ids[j], node_ids[i], sign=int(matrix[i][j]))
    _check_signs(G)
    nx.set_node_attributes(G, "state", "category")
    nx.freeze(G)
    return G

load_digraph

load_digraph(model)

Load a built-in example model as a signed directed graph.

Parameters:

Name	Type	Description	Default
model	str	Name of the built-in model to load. Available models: - "snowshoe": Simple 3-node predator-prey model (R, C, P) - "snowshoe_rp": Snowshoe model with an added positive R->P link - "snowshoe_io": Snowshoe_rp model with input/output nodes - "chain": 5-node linear chain with self-effects - "mesocosm": 8-node complex ecosystem model	required

Returns:

Type	Description
DiGraph	nx.DiGraph: A NetworkX directed graph with signed edges.

Raises:

Type	Description
ValueError	If model name is not recognized.

Examples:

from qmm import load_digraph
G = load_digraph("snowshoe_rp")
list(G.nodes())
# ['R', 'C', 'P']

list(G.edges(data='sign'))
# [('R', 'R', -1), ('R', 'C', 1), ('R', 'P', 1), ('C', 'R', -1), ('C', 'P', 1), ('P', 'C', -1), ('P', 'P', -1)]

Source code in qmm/core/helper.py

def load_digraph(model: str) -> nx.DiGraph:
    """Load a built-in example model as a signed directed graph.

    Args:
        model: Name of the built-in model to load. Available models:
            - "snowshoe": Simple 3-node predator-prey model (R, C, P)
            - "snowshoe_rp": Snowshoe model with an added positive R->P link
            - "snowshoe_io": Snowshoe_rp model with input/output nodes
            - "chain": 5-node linear chain with self-effects
            - "mesocosm": 8-node complex ecosystem model

    Returns:
        nx.DiGraph: A NetworkX directed graph with signed edges.

    Raises:
        ValueError: If model name is not recognized.

    Examples:
        ```python
        from qmm import load_digraph
        G = load_digraph("snowshoe_rp")
        list(G.nodes())
        # ['R', 'C', 'P']

        list(G.edges(data='sign'))
        # [('R', 'R', -1), ('R', 'C', 1), ('R', 'P', 1), ('C', 'R', -1), ('C', 'P', 1), ('P', 'C', -1), ('P', 'P', -1)]
        ```
    """
    models = {
        "snowshoe": {
            "matrix": [[-1, -1, 0], [1, 0, -1], [0, 1, -1]],
            "labels": ['R', 'C', 'P']
        },
        "snowshoe_rp": {
            "matrix": [[-1, -1, 0], [1, 0, -1], [1, 1, -1]],
            "labels": ['R', 'C', 'P']
        },
        "chain": {
            "matrix": [[-1, -1, 0, 0, 0], [1, -1, -1, 0, 0], [0, 1, -1, -1, 0], [0, 0, 1, -1, -1], [0, 0, 0, 1, -1]],
            "labels": ['1', '2', '3', '4', '5']
        },
        "mesocosm": {
            "matrix": [
                [-1, -1, -1, -1, 0, 0, 0, 0],
                [1, 0, 0, 0, -1, -1, 0, 0],
                [1, 0, 0, 0, 0, -1, 0, 0],
                [1, 0, 0, -1, 0, 0, 0, 0],
                [0, 1, 0, 0, 0, 0, -1, -1],
                [0, 1, 1, 0, 0, 0, 0, -1],
                [0, 0, 0, 0, 1, 0, 0, -1],
                [0, 0, 0, 0, 1, 1, 1, -1],
            ],
            "labels": ['P', 'A1', 'A2', 'AP', 'H1', 'H2', 'C1', 'C2']
        }
    }

    if model == "snowshoe_io":
        G = nx.DiGraph()
        for node in ['R', 'C', 'P']:
            G.add_node(node, category='state')
        for node in ['Inp1', 'Inp2']:
            G.add_node(node, category='input')
        for node in ['Out1', 'Out2']:
            G.add_node(node, category='output')
        edges = [
            ('R', 'R', -1),
            ('R', 'C', 1),
            ('C', 'R', -1),
            ('C', 'P', 1),
            ('P', 'C', -1),
            ('P', 'P', -1),
            ('Inp1', 'R', 1),
            ('Inp1', 'C', -1),
            ('Inp2', 'P', -1),
            ('C', 'Out1', -1),
            ('C', 'Out2', 1),
            ('P', 'Out1', 1),
        ]
        for source, target, sign in edges:
            G.add_edge(source, target, sign=sign)
        nx.freeze(G)
        return G

    if model not in models:
        available_models = list(models.keys()) + ["snowshoe_io"]
        available = ', '.join(f'"{m}"' for m in sorted(available_models))
        raise ValueError(f"Model '{model}' not found. Available models: {available}")

    m = models[model]
    G = list_to_digraph(m["matrix"], m["labels"])
    return G

digraph_to_list

digraph_to_list(G)

Convert a directed graph to an adjacency matrix string representation.

Parameters:

Name	Type	Description	Default
G	DiGraph	A NetworkX directed graph with signed edges.	required

Returns:

Name	Type	Description
str	str	String representation of the adjacency matrix.

Examples:

from qmm import load_digraph, digraph_to_list
digraph_to_list(load_digraph("snowshoe_rp"))
# '[[0, -1, 1], [1, -1, 1], [-1, 0, -1]]'

Source code in qmm/core/helper.py

def digraph_to_list(G: nx.DiGraph) -> str:
    """Convert a directed graph to an adjacency matrix string representation.

    Args:
        G: A NetworkX directed graph with signed edges.

    Returns:
        str: String representation of the adjacency matrix.

    Examples:
        ```python
        from qmm import load_digraph, digraph_to_list
        digraph_to_list(load_digraph("snowshoe_rp"))
        # '[[0, -1, 1], [1, -1, 1], [-1, 0, -1]]'
        ```
    """
    if not isinstance(G, nx.DiGraph):
        raise TypeError("Input must be a networkx.DiGraph.")
    n = G.number_of_nodes()
    nodes = sorted(G.nodes())
    node_to_index = {node: i for i, node in enumerate(nodes)}
    matrix = [[0 for _ in range(n)] for _ in range(n)]
    for source, target, data in G.edges(data=True):
        i, j = node_to_index[source], node_to_index[target]
        sign = data.get("sign", 1)
        matrix[j][i] = sign
    return str(matrix)

get_nodes

get_nodes(G, node_type='state', labels=False)

Get nodes of a specific type from a directed graph.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX directed graph to extract nodes from.	required
node_type	Literal['state', 'input', 'output', 'all']	Type of nodes to extract ('state', 'input', 'output', or 'all').	'state'
labels	bool	If True, return node labels instead of node ids.	False

Returns:

Type	Description
List[Union[str, Dict[str, Any]]]	List of node identifiers or dictionaries containing node data.

Examples:

from qmm import load_digraph, get_nodes
get_nodes(load_digraph("snowshoe_rp"), "state")
# ['R', 'C', 'P']

Source code in qmm/core/helper.py

def get_nodes(
    G: nx.DiGraph,
    node_type: Literal["state", "input", "output", "all"] = "state",
    labels: bool = False,
) -> List[Union[str, Dict[str, Any]]]:
    """Get nodes of a specific type from a directed graph.

    Args:
        G: NetworkX directed graph to extract nodes from.
        node_type: Type of nodes to extract ('state', 'input', 'output', or 'all').
        labels: If True, return node labels instead of node ids.

    Returns:
        List of node identifiers or dictionaries containing node data.

    Examples:
        ```python
        from qmm import load_digraph, get_nodes
        get_nodes(load_digraph("snowshoe_rp"), "state")
        # ['R', 'C', 'P']
        ```
    """
    if not isinstance(G, nx.DiGraph):
        raise TypeError("Input must be a networkx.DiGraph.")

    if node_type == "all":
        return [n if not labels else d.get("label", n) for n, d in G.nodes(data=True)]
    else:
        return [n if not labels else d.get("label", n) for n, d in G.nodes(data=True) if d.get("category") == node_type]

get_weight

get_weight(net, absolute, no_effect=sp.nan)

Calculate weight matrix by dividing net effect by absolute effect.

Parameters:

Name	Type	Description	Default
net	Matrix	Matrix of net terms.	required
absolute	Matrix	Matrix of absolute terms.	required
no_effect	Union[Basic, float]	Value to use when absolute terms is 0 (default: sympy.nan).	nan

Returns:

Type	Description
Matrix	sympy.Matrix: Matrix of weights.

Examples:

import sympy as sp
from qmm import get_weight
net = sp.Matrix([[2, -2], [1, 0]])
absolute = sp.Matrix([[4, 2], [1, 0]])
get_weight(net, absolute)
# Matrix([
# [1/2,  -1],
# [  1, nan]])

Source code in qmm/core/helper.py

def get_weight(net: sp.Matrix, absolute: sp.Matrix, no_effect: Union[sp.Basic, float] = sp.nan) -> sp.Matrix:
    """Calculate weight matrix by dividing net effect by absolute effect.

    Args:
        net: Matrix of net terms.
        absolute: Matrix of absolute terms.
        no_effect: Value to use when absolute terms is 0 (default: sympy.nan).

    Returns:
        sympy.Matrix: Matrix of weights.

    Examples:
        ```python
        import sympy as sp
        from qmm import get_weight
        net = sp.Matrix([[2, -2], [1, 0]])
        absolute = sp.Matrix([[4, 2], [1, 0]])
        get_weight(net, absolute)
        # Matrix([
        # [1/2,  -1],
        # [  1, nan]])
        ```
    """
    if net.shape != absolute.shape:
        raise ValueError("Matrices must have the same shape")
    result = sp.zeros(*net.shape)
    for i in range(net.shape[0]):
        for j in range(net.shape[1]):
            if absolute[i, j] == 0:
                result[i, j] = no_effect
            else:
                result[i, j] = net[i, j] / absolute[i, j]
    return result

get_positive

get_positive(net, absolute)

Calculate matrix of positive terms.

Parameters:

Name	Type	Description	Default
net	Matrix	Matrix of net terms.	required
absolute	Matrix	Matrix of absolute terms.	required

Returns:

Type	Description
Matrix	sympy.Matrix: Matrix of positive terms.

Examples:

import sympy as sp
from qmm import get_positive
net = sp.Matrix([[3, -2], [1, 0]])
absolute = sp.Matrix([[4, 2], [1, 0]])
get_positive(net, absolute)
# Matrix([
# [3, 0],
# [1, 0]])

Source code in qmm/core/helper.py

def get_positive(net: sp.Matrix, absolute: sp.Matrix) -> sp.Matrix:
    """Calculate matrix of positive terms.

    Args:
        net: Matrix of net terms.
        absolute: Matrix of absolute terms.

    Returns:
        sympy.Matrix: Matrix of positive terms.

    Examples:
        ```python
        import sympy as sp
        from qmm import get_positive
        net = sp.Matrix([[3, -2], [1, 0]])
        absolute = sp.Matrix([[4, 2], [1, 0]])
        get_positive(net, absolute)
        # Matrix([
        # [3, 0],
        # [1, 0]])
        ```
    """
    if net.shape != absolute.shape:
        raise ValueError("Matrices must have the same shape")
    result = sp.zeros(*net.shape)
    for i in range(net.shape[0]):
        for j in range(net.shape[1]):
            result[i, j] = (net[i, j] + absolute[i, j]) // 2
    return result

get_negative

get_negative(net, absolute)

Calculate matrix of negative terms.

Parameters:

Name	Type	Description	Default
net	Matrix	Matrix of net terms.	required
absolute	Matrix	Matrix of absolute terms.	required

Returns:

Type	Description
Matrix	sympy.Matrix: Matrix of negative terms.

Examples:

import sympy as sp
from qmm import get_negative
net = sp.Matrix([[3, -2], [1, 0]])
absolute = sp.Matrix([[4, 2], [1, 0]])
get_negative(net, absolute)
# Matrix([
# [0, 2],
# [0, 0]])

Source code in qmm/core/helper.py

def get_negative(net: sp.Matrix, absolute: sp.Matrix) -> sp.Matrix:
    """Calculate matrix of negative terms.

    Args:
        net: Matrix of net terms.
        absolute: Matrix of absolute terms.

    Returns:
        sympy.Matrix: Matrix of negative terms.

    Examples:
        ```python
        import sympy as sp
        from qmm import get_negative
        net = sp.Matrix([[3, -2], [1, 0]])
        absolute = sp.Matrix([[4, 2], [1, 0]])
        get_negative(net, absolute)
        # Matrix([
        # [0, 2],
        # [0, 0]])
        ```
    """
    if net.shape != absolute.shape:
        raise ValueError("Matrices must have the same shape")
    result = sp.zeros(*net.shape)
    for i in range(net.shape[0]):
        for j in range(net.shape[1]):
            result[i, j] = (absolute[i, j] - net[i, j]) // 2
    return result

sign_determinacy

sign_determinacy(wmat, tmat, method='average')

Calculate sign determinacy matrix from prediction weights.

Parameters:

Name	Type	Description	Default
wmat	Matrix	Matrix of prediction weights.	required
tmat	Matrix	Matrix of absolute feedback.	required
method	Literal['average', '95_bound']	Method to use for probability calculation ('average' or '95_bound').	'average'

Returns:

Type	Description
Matrix	sympy.Matrix: Probability of sign determinacy.

References

• Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.

Examples:

from qmm import load_digraph, weighted_predictions_matrix, absolute_feedback_matrix, sign_determinacy
G = load_digraph("snowshoe_rp")
wmat = weighted_predictions_matrix(G)
tmat = absolute_feedback_matrix(G)
sign_determinacy(wmat, tmat, method='average')
# Matrix([
# [  1,  -1,  1],
# [1/2,   1, -1],
# [  1, 1/2,  1]])

Source code in qmm/core/helper.py

def sign_determinacy(
    wmat: sp.Matrix,
    tmat: sp.Matrix,
    method: Literal["average", "95_bound"] = "average",
) -> sp.Matrix:
    """Calculate sign determinacy matrix from prediction weights.

    Args:
        wmat: Matrix of prediction weights.
        tmat: Matrix of absolute feedback.
        method: Method to use for probability calculation ('average' or '95_bound').

    Returns:
        sympy.Matrix: Probability of sign determinacy.

    References:
        - Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.

    Examples:
        ```python
        from qmm import load_digraph, weighted_predictions_matrix, absolute_feedback_matrix, sign_determinacy
        G = load_digraph("snowshoe_rp")
        wmat = weighted_predictions_matrix(G)
        tmat = absolute_feedback_matrix(G)
        sign_determinacy(wmat, tmat, method='average')
        # Matrix([
        # [  1,  -1,  1],
        # [1/2,   1, -1],
        # [  1, 1/2,  1]])
        ```
    """

    MAX_PROB = sp.Float('0.999999')

    def compute_prob(w, t, method):
        if t == sp.Integer(0):
            return sp.nan
        return compute_prob_average(w, t) if method == "average" else compute_prob_95_bound(w, t)

    def compute_prob_average(w, t):
        bw = 3.45962
        bwt = 0.03417
        w_float = float(w)
        t_float = float(t)
        exponent = bw * w_float + bwt * w_float * t_float

        if exponent > 700:
            return MAX_PROB

        prob_float = np.exp(exponent) / (1 + np.exp(exponent))
        prob = sp.Float(prob_float)

        prob = max(sp.Rational(1, 2), prob)

        if prob >= MAX_PROB:
            prob = MAX_PROB
        return prob

    def compute_prob_95_bound(w, t):
        bw = 9.766
        bwt = 0.139
        w_float = float(w)
        t_float = float(t)
        exponent = bw * w_float + bwt * w_float * t_float

        if exponent > 700:
            return MAX_PROB

        prob_float = np.exp(exponent) / (1253.992 + np.exp(exponent))
        prob = sp.Float(prob_float)

        prob = max(sp.Rational(1, 2), prob)
        if prob >= MAX_PROB:
            prob = MAX_PROB
        return prob

    if method not in ["average", "95_bound"]:
        raise ValueError("Invalid method. Choose 'average' or '95_bound'.")
    rows, cols = wmat.shape
    def calc_prob(i, j):
        w, t = wmat[i, j], tmat[i, j]
        if w.is_zero:
            return sp.Rational(1, 2)
        if sp.Abs(w) == sp.Integer(1):
            return sp.sign(w) * sp.Integer(1)
        prob = compute_prob(sp.Abs(w), t, method)
        return sp.sign(w) * prob if prob is not None else sp.nan

    pmat = sp.Matrix(rows, cols, lambda i, j: calc_prob(i, j))
    return pmat

perm

perm(A, method='bbfg', decompose=True)

Compute the permanent of a square matrix.

The permanent is similar to the determinant but uses only addition (no sign alternation). This implementation is based on the algorithms from thewalrus library (https://github.com/XanaduAI/thewalrus).

Parameters:

Name	Type	Description	Default
A	ndarray	A square numpy array (float or complex).	required
method	Literal['bbfg', 'ryser']	Algorithm to use - "bbfg" for BBFG formula (default, faster) or "ryser" for Ryser formula. Any other value uses Ryser.	'bbfg'
decompose	bool	Dulmage-Mendelsohn decomposition (default True) for sparse matrices.	True

Returns:

Type	Description
float	The permanent of matrix A.

Raises:

Type	Description
TypeError	If input is not a numpy array.
ValueError	If matrix is not square or contains NaNs.

References

• Ryser, H.J. (1963). Combinatorial Mathematics.
• Glynn, D.G. (2010). The permanent of a square matrix. European Journal of Combinatorics 31, 1887–1891.

Examples:

import numpy as np
from qmm.core.helper import perm
perm(np.array([[1, 2], [3, 4]]), method='bbfg')
# 10

Source code in qmm/core/helper.py

def perm(
    A: np.ndarray, method: Literal["bbfg", "ryser"] = "bbfg", decompose: bool = True
) -> float:
    """Compute the permanent of a square matrix.

    The permanent is similar to the determinant but uses only addition
    (no sign alternation). This implementation is based on the algorithms
    from thewalrus library (https://github.com/XanaduAI/thewalrus).

    Args:
        A: A square numpy array (float or complex).
        method: Algorithm to use - "bbfg" for BBFG formula (default, faster)
                or "ryser" for Ryser formula. Any other value uses Ryser.
        decompose: Dulmage-Mendelsohn decomposition (default True) for sparse matrices.

    Returns:
        The permanent of matrix A.

    Raises:
        TypeError: If input is not a numpy array.
        ValueError: If matrix is not square or contains NaNs.

    References:
        - Ryser, H.J. (1963). Combinatorial Mathematics.
        - Glynn, D.G. (2010). The permanent of a square matrix. European Journal of Combinatorics 31, 1887–1891.

    Examples:
        ```python
        import numpy as np
        from qmm.core.helper import perm
        perm(np.array([[1, 2], [3, 4]]), method='bbfg')
        # 10
        ```
    """
    if not isinstance(A, np.ndarray):
        raise TypeError("Input matrix must be a NumPy array.")

    matshape = A.shape
    if matshape[0] != matshape[1]:
        raise ValueError("Input matrix must be square.")
    if np.isnan(A).any():
        raise ValueError("Input matrix must not contain NaNs.")

    if matshape[0] == 0:
        return A.dtype.type(1.0)
    if matshape[0] == 1:
        return A[0, 0]
    if matshape[0] == 2:
        return A[0, 0] * A[1, 1] + A[0, 1] * A[1, 0]
    if matshape[0] == 3:
        return (
            A[0, 2] * A[1, 1] * A[2, 0]
            + A[0, 1] * A[1, 2] * A[2, 0]
            + A[0, 2] * A[1, 0] * A[2, 1]
            + A[0, 0] * A[1, 2] * A[2, 1]
            + A[0, 1] * A[1, 0] * A[2, 2]
            + A[0, 0] * A[1, 1] * A[2, 2]
        )

    overflow = bool(np.prod(np.abs(A).sum(axis=0, dtype=float)) > 2.0**53)
    if overflow and int((A != 0).sum()) <= 8 * matshape[0] and not np.mod(A, 1).any():
        return _perm_int(A)

    if decompose:
        S = csr_matrix(A != 0)
        col_of = maximum_bipartite_matching(S, perm_type="column")
        if (col_of < 0).any():
            return 0
        nb, labels = connected_components(S[:, col_of], connection="strong")
        if nb > 1:
            result = 1
            for b in range(nb):
                rk = np.flatnonzero(labels == b)
                blk = np.ascontiguousarray(A[np.ix_(rk, col_of[rk])])
                result *= perm(blk, method, decompose=False)
            return result
    if overflow:
        raise OverflowError("perm exceeds float precision (2**53)")
    return _perm_bbfg(A) if method == "bbfg" else _perm_ryser(A)

get_dashed_alternatives

get_dashed_alternatives(G, combinations=True)

Generate all alternative model structures based on dashed edges.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph with potentially dashed edges (edges with dashes=True attribute)	required
combinations	bool	If True, return all 2^n combinations of dashed edges. If False, return base graph (no dashed edges) plus variants with each single dashed edge added.	True

Returns:

Type	Description
List[DiGraph]	List[nx.DiGraph]: List of graph variants with different dashed edge configurations. If no dashed edges exist, returns a list containing only the original graph.

References

• Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.

Examples:

from qmm import load_digraph
from qmm.core.helper import get_dashed_alternatives
G_mod = load_digraph("snowshoe_rp").copy()
G_mod.remove_edge('C', 'P')
G_mod.add_edge('C', 'P', sign=1, dashes=True)
variants = get_dashed_alternatives(G_mod, combinations=True)

len(variants)
# 2

Source code in qmm/core/helper.py

def get_dashed_alternatives(G: nx.DiGraph, combinations: bool = True) -> List[nx.DiGraph]:
    """Generate all alternative model structures based on dashed edges.

    Args:
        G: NetworkX DiGraph with potentially dashed edges (edges with dashes=True attribute)
        combinations: If True, return all 2^n combinations of dashed edges.
                     If False, return base graph (no dashed edges) plus variants with each single dashed edge added.

    Returns:
        List[nx.DiGraph]: List of graph variants with different dashed edge configurations.
                         If no dashed edges exist, returns a list containing only the original graph.

    References:
        - Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.

    Examples:
        ```python
        from qmm import load_digraph
        from qmm.core.helper import get_dashed_alternatives
        G_mod = load_digraph("snowshoe_rp").copy()
        G_mod.remove_edge('C', 'P')
        G_mod.add_edge('C', 'P', sign=1, dashes=True)
        variants = get_dashed_alternatives(G_mod, combinations=True)

        len(variants)
        # 2
        ```
    """
    dashed_edges = [(j, i, d) for j, i, d in G.edges(data=True) if d.get("dashes", False)]

    if not dashed_edges:
        return [G]

    if combinations:
        mask_values = range(2 ** len(dashed_edges))
        variants = []
        for mask in mask_values:
            G_variant = G.copy()
            for idx, (j, i, _) in enumerate(dashed_edges):
                include_edge = bool(mask & (1 << idx))
                if not include_edge:
                    G_variant.remove_edge(j, i)
            variants.append(G_variant)
    else:
        G_base = G.copy()
        for j, i, _ in dashed_edges:
            G_base.remove_edge(j, i)

        variants = [G_base]

        for j, i, edge_data in dashed_edges:
            G_variant = G_base.copy()
            G_variant.add_edge(j, i, **edge_data)
            variants.append(G_variant)

    return variants

Extension module

qmm.extensions.senstability

Analyse the sensitivity of system stability to direct effects within feedback cycles.

structural_sensitivity cached

structural_sensitivity(G, level=None)

Calculate contribution of direct effects to stabilising and destabilising feedback.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Feedback level (None for highest level)	None

Returns:

Type	Description
Matrix	sp.Matrix: Ratio of net to total feedback terms for each direct effect

References

• Hosack, G.R., Li, H.W., Rossignol, P.A. (2009). Sensitivity of system stability to model structure. Ecological Modelling 220, 1054–1062.

Examples:

from qmm import structural_sensitivity, load_digraph
structural_sensitivity(load_digraph("snowshoe_rp"))
# Matrix([
# [-a_C,P*a_P,C*a_R,R, a_C,P*a_P,R*a_R,C - a_C,R*a_P,P*a_R,C,                                      0],
# [-a_C,R*a_P,P*a_R,C,                                     0, -a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C],
# [ a_C,P*a_P,R*a_R,C,                    -a_C,P*a_P,C*a_R,R,                     -a_C,R*a_P,P*a_R,C]])

Source code in qmm/extensions/senstability.py

@cache
def structural_sensitivity(G: nx.DiGraph, level: Optional[int] = None) -> sp.Matrix:
    """Calculate contribution of direct effects to stabilising and destabilising feedback.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Feedback level (None for highest level)

    Returns:
        sp.Matrix: Ratio of net to total feedback terms for each direct effect

    References:
        - Hosack, G.R., Li, H.W., Rossignol, P.A. (2009). Sensitivity of system stability to model structure. Ecological Modelling 220, 1054–1062.

    Examples:
        ```python
        from qmm import structural_sensitivity, load_digraph
        structural_sensitivity(load_digraph("snowshoe_rp"))
        # Matrix([
        # [-a_C,P*a_P,C*a_R,R, a_C,P*a_P,R*a_R,C - a_C,R*a_P,P*a_R,C,                                      0],
        # [-a_C,R*a_P,P*a_R,C,                                     0, -a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C],
        # [ a_C,P*a_P,R*a_R,C,                    -a_C,P*a_P,C*a_R,R,                     -a_C,R*a_P,P*a_R,C]])
        ```
    """
    return _structural_sensitivity(G, level, system_feedback)

net_structural_sensitivity cached

net_structural_sensitivity(G, level=None)

Calculate net contribution of direct effects to system feedback.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Feedback level (None for highest level)	None

Returns:

Type	Description
Matrix	sp.Matrix: Net feedback terms containing each direct effect

References

• Hosack, G.R., Li, H.W., Rossignol, P.A. (2009). Sensitivity of system stability to model structure. Ecological Modelling 220, 1054–1062.

Examples:

from qmm import net_structural_sensitivity, load_digraph
net_structural_sensitivity(load_digraph("snowshoe_rp"))
# Matrix([
# [-1,  0,  0],
# [-1,  0,  0],
# [ 1, -1, -1]])

Source code in qmm/extensions/senstability.py

@cache
def net_structural_sensitivity(G: nx.DiGraph, level: Optional[int] = None) -> sp.Matrix:
    """Calculate net contribution of direct effects to system feedback.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Feedback level (None for highest level)

    Returns:
        sp.Matrix: Net feedback terms containing each direct effect

    References:
        - Hosack, G.R., Li, H.W., Rossignol, P.A. (2009). Sensitivity of system stability to model structure. Ecological Modelling 220, 1054–1062.

    Examples:
        ```python
        from qmm import net_structural_sensitivity, load_digraph
        net_structural_sensitivity(load_digraph("snowshoe_rp"))
        # Matrix([
        # [-1,  0,  0],
        # [-1,  0,  0],
        # [ 1, -1, -1]])
        ```
    """
    return _structural_sensitivity(G, level, net_feedback)

absolute_structural_sensitivity cached

absolute_structural_sensitivity(G, level=None)

Calculate total contribution of direct effects to system feedback.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Feedback level (None for highest level)	None

Returns:

Type	Description
Matrix	sp.Matrix: Total feedback terms containing each direct effect

References

• Hosack, G.R., Li, H.W., Rossignol, P.A. (2009). Sensitivity of system stability to model structure. Ecological Modelling 220, 1054–1062.

Examples:

from qmm import absolute_structural_sensitivity, load_digraph
absolute_structural_sensitivity(load_digraph("snowshoe_rp"))
# Matrix([
# [1, 2, 0],
# [1, 0, 2],
# [1, 1, 1]])

Source code in qmm/extensions/senstability.py

@cache
def absolute_structural_sensitivity(G: nx.DiGraph, level: Optional[int] = None) -> sp.Matrix:
    """Calculate total contribution of direct effects to system feedback.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Feedback level (None for highest level)

    Returns:
        sp.Matrix: Total feedback terms containing each direct effect

    References:
        - Hosack, G.R., Li, H.W., Rossignol, P.A. (2009). Sensitivity of system stability to model structure. Ecological Modelling 220, 1054–1062.

    Examples:
        ```python
        from qmm import absolute_structural_sensitivity, load_digraph
        absolute_structural_sensitivity(load_digraph("snowshoe_rp"))
        # Matrix([
        # [1, 2, 0],
        # [1, 0, 2],
        # [1, 1, 1]])
        ```
    """
    return _structural_sensitivity(G, level, absolute_feedback)

weighted_structural_sensitivity cached

weighted_structural_sensitivity(G, level=None)

Calculate weighted structual sensitvity for each direct effect.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
level	Optional[int]	Feedback level (None for highest level)	None

Returns:

Type	Description
Matrix	sp.Matrix: Weighted structural sensitivity of each direct effect

References

• Hosack, G.R., Li, H.W., Rossignol, P.A. (2009). Sensitivity of system stability to model structure. Ecological Modelling 220, 1054–1062.

Examples:

from qmm import weighted_structural_sensitivity, load_digraph
weighted_structural_sensitivity(load_digraph("snowshoe_rp"))
# Matrix([
# [-1,   0, nan],
# [-1, nan,   0],
# [ 1,  -1,  -1]])

Source code in qmm/extensions/senstability.py

@cache
def weighted_structural_sensitivity(G: nx.DiGraph, level: Optional[int] = None) -> sp.Matrix:
    """Calculate weighted structual sensitvity for each direct effect.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        level: Feedback level (None for highest level)

    Returns:
        sp.Matrix: Weighted structural sensitivity of each direct effect

    References:
        - Hosack, G.R., Li, H.W., Rossignol, P.A. (2009). Sensitivity of system stability to model structure. Ecological Modelling 220, 1054–1062.

    Examples:
        ```python
        from qmm import weighted_structural_sensitivity, load_digraph
        weighted_structural_sensitivity(load_digraph("snowshoe_rp"))
        # Matrix([
        # [-1,   0, nan],
        # [-1, nan,   0],
        # [ 1,  -1,  -1]])
        ```
    """
    net = net_structural_sensitivity(G, level)
    absolute = absolute_structural_sensitivity(G, level)
    return get_weight(net, absolute)

qmm.extensions.life

Analyse change in life expectancy from press perturbations.

birth_matrix cached

birth_matrix(G, form='symbolic', perturb=None)

Create matrix of direct effects on birth rate from press perturbations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
form	Literal['symbolic', 'signed']	Type of computation ('symbolic', 'signed')	'symbolic'

Returns:

Type	Description
Matrix	sp.Matrix: Positive direct effects on birth rate

References

• Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

Examples:

from qmm import birth_matrix, load_digraph
birth_matrix(load_digraph("snowshoe_rp"), form='symbolic')
# Matrix([
# [    0,     0, 0],
# [a_C,R,     0, 0],
# [a_P,R, a_P,C, 0]])

birth_matrix(load_digraph("snowshoe_rp"), form='signed')
# Matrix([
# [0, 0, 0],
# [1, 0, 0],
# [1, 1, 0]])

Source code in qmm/extensions/life.py

@cache
def birth_matrix(
    G: nx.DiGraph,
    form: Literal["symbolic", "signed"] = "symbolic",
    perturb: Optional[str] = None,
) -> sp.Matrix:
    """Create matrix of direct effects on birth rate from press perturbations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        form: Type of computation ('symbolic', 'signed')

    Returns:
        sp.Matrix: Positive direct effects on birth rate

    References:
        - Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

    Examples:
        ```python
        from qmm import birth_matrix, load_digraph
        birth_matrix(load_digraph("snowshoe_rp"), form='symbolic')
        # Matrix([
        # [    0,     0, 0],
        # [a_C,R,     0, 0],
        # [a_P,R, a_P,C, 0]])

        birth_matrix(load_digraph("snowshoe_rp"), form='signed')
        # Matrix([
        # [0, 0, 0],
        # [1, 0, 0],
        # [1, 1, 0]])
        ```
    """
    A_sgn = create_matrix(G, form="signed")
    A_sym = create_matrix(G, form="symbolic")
    nodes = get_nodes(G, "state")
    if perturb is not None and perturb not in nodes:
        raise ValueError(f"Perturbation node must be one of: {nodes}")
    n = len(nodes)
    def birth_element(i, j):
        if form == "symbolic":
            return A_sym[i, j] if A_sgn[i, j] > 0 else 0
        else:
            return sp.Integer(1) if A_sgn[i, j] > 0 else 0
    if perturb is not None:
        src_id = nodes.index(perturb)
        return sp.Matrix(n, 1, lambda i, j: birth_element(i, src_id))
    else:
        return sp.Matrix(n, n, lambda i, j: birth_element(i, j))

death_matrix cached

death_matrix(G, form='symbolic', perturb=None)

Create matrix of direct effects on death rate from press perturbations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
form	Literal['symbolic', 'signed']	Type of computation ('symbolic', 'signed')	'symbolic'

Returns:

Type	Description
Matrix	sp.Matrix: Positive direct effects on death rate

References

• Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

Examples:

from qmm import death_matrix, load_digraph
death_matrix(load_digraph("snowshoe_rp"), form='symbolic')
# Matrix([
# [a_R,R, a_R,C,     0],
# [    0,     0, a_C,P],
# [    0,     0, a_P,P]])

death_matrix(load_digraph("snowshoe_rp"), form='signed')
# Matrix([
# [1, 1, 0],
# [0, 0, 1],
# [0, 0, 1]])

Source code in qmm/extensions/life.py

@cache
def death_matrix(
    G: nx.DiGraph,
    form: Literal["symbolic", "signed"] = "symbolic",
    perturb: Optional[str] = None,
) -> sp.Matrix:
    """Create matrix of direct effects on death rate from press perturbations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        form: Type of computation ('symbolic', 'signed')

    Returns:
        sp.Matrix: Positive direct effects on death rate

    References:
        - Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

    Examples:
        ```python
        from qmm import death_matrix, load_digraph
        death_matrix(load_digraph("snowshoe_rp"), form='symbolic')
        # Matrix([
        # [a_R,R, a_R,C,     0],
        # [    0,     0, a_C,P],
        # [    0,     0, a_P,P]])

        death_matrix(load_digraph("snowshoe_rp"), form='signed')
        # Matrix([
        # [1, 1, 0],
        # [0, 0, 1],
        # [0, 0, 1]])
        ```
    """
    A_sgn = create_matrix(G, form="signed")
    A_sym = create_matrix(G, form="symbolic")
    nodes = get_nodes(G, "state")
    if perturb is not None and perturb not in nodes:
        raise ValueError(f"Perturbation node must be one of: {nodes}")
    n = len(nodes)
    def death_element(i, j):
        if form == "symbolic":
            return A_sym[i, j] * sp.Integer(-1) if A_sgn[i, j] < 0 else 0
        else:
            return sp.Integer(1) if A_sgn[i, j] < 0 else 0
    if perturb is not None:
        src_id = nodes.index(perturb)
        return sp.Matrix(n, 1, lambda i, j: death_element(i, src_id))
    else:
        return sp.Matrix(n, n, lambda i, j: death_element(i, j))

life_expectancy_change cached

life_expectancy_change(G, form='symbolic', type='birth', perturb=None)

Calculate change in life expectancy from press perturbations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
type	Literal['birth', 'death']	Change in birth or death rate ('birth' or 'death')	'birth'
form	Literal['symbolic', 'signed']	Type of computation ('symbolic', 'signed')	'symbolic'
perturb	Optional[str]	Node to perturb (None for full matrix)	None

Returns:

Type	Description
Matrix	sp.Matrix: Change in life expectancy for each component

References

• Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

Examples:

from qmm import life_expectancy_change, load_digraph
life_expectancy_change(load_digraph("snowshoe_rp"), form='symbolic', type='birth')
# Matrix([
# [-a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C - a_C,R*a_P,P*a_R,C,                                      0,                  0],
# [                                        -a_C,P*a_C,R*a_P,C, -a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C, -a_C,P*a_C,R*a_R,C],
# [                                        -a_C,R*a_P,C*a_P,P, -a_P,C*a_P,P*a_R,R + a_P,P*a_P,R*a_R,C, -a_C,R*a_P,P*a_R,C]])

life_expectancy_change(load_digraph("snowshoe_rp"), form='symbolic', type='death')
# Matrix([
# [                 0,                                      0,                                     0],
# [-a_C,P*a_C,R*a_P,C,                      a_C,R*a_P,P*a_R,C,                    -a_C,P*a_C,R*a_R,C],
# [-a_C,R*a_P,C*a_P,P, -a_P,C*a_P,P*a_R,R + a_P,P*a_P,R*a_R,C, a_C,P*a_P,C*a_R,R - a_C,P*a_P,R*a_R,C]])

Source code in qmm/extensions/life.py

@cache
def life_expectancy_change(
    G: nx.DiGraph,
    form: Literal["symbolic", "signed"] = "symbolic",
    type: Literal["birth", "death"] = "birth",
    perturb: Optional[str] = None,
) -> sp.Matrix:
    """Calculate change in life expectancy from press perturbations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        type: Change in birth or death rate ('birth' or 'death')
        form: Type of computation ('symbolic', 'signed')
        perturb: Node to perturb (None for full matrix)

    Returns:
        sp.Matrix: Change in life expectancy for each component

    References:
        - Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

    Examples:
        ```python
        from qmm import life_expectancy_change, load_digraph
        life_expectancy_change(load_digraph("snowshoe_rp"), form='symbolic', type='birth')
        # Matrix([
        # [-a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C - a_C,R*a_P,P*a_R,C,                                      0,                  0],
        # [                                        -a_C,P*a_C,R*a_P,C, -a_C,P*a_P,C*a_R,R + a_C,P*a_P,R*a_R,C, -a_C,P*a_C,R*a_R,C],
        # [                                        -a_C,R*a_P,C*a_P,P, -a_P,C*a_P,P*a_R,R + a_P,P*a_P,R*a_R,C, -a_C,R*a_P,P*a_R,C]])

        life_expectancy_change(load_digraph("snowshoe_rp"), form='symbolic', type='death')
        # Matrix([
        # [                 0,                                      0,                                     0],
        # [-a_C,P*a_C,R*a_P,C,                      a_C,R*a_P,P*a_R,C,                    -a_C,P*a_C,R*a_R,C],
        # [-a_C,R*a_P,C*a_P,P, -a_P,C*a_P,P*a_R,R + a_P,P*a_P,R*a_R,C, a_C,P*a_P,C*a_R,R - a_C,P*a_P,R*a_R,C]])
        ```
    """
    amat = adjoint_matrix(G, form=form)
    if type == "birth":
        matrix = death_matrix(G, form=form)
    elif type == "death":
        matrix = birth_matrix(G, form=form)
    else:
        raise ValueError("type must be either 'birth' or 'death'")
    result = sp.expand(sp.Integer(-1) * matrix * amat)
    if perturb is not None:
        nodes = get_nodes(G, "state")
        if perturb not in nodes:
            raise ValueError(f"Perturbation node must be one of: {nodes}")
        perturb_index = nodes.index(perturb)
        return result.col(perturb_index)
    return result

net_life_expectancy_change cached

net_life_expectancy_change(G, type='birth')

Calculate net terms in life expectancy change from press perturbations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
type	Literal['birth', 'death']	Change in birth or death rate ('birth' or 'death')	'birth'

Returns:

Type	Description
Matrix	sp.Matrix: Net life expectancy change for each component

References

• Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

Examples:

from qmm import net_life_expectancy_change, load_digraph
net_life_expectancy_change(load_digraph("snowshoe_rp"), type='birth')
# Matrix([
# [-1, 0,  0],
# [-1, 0, -1],
# [-1, 0, -1]])

net_life_expectancy_change(load_digraph("snowshoe_rp"), type='death')
# Matrix([
# [ 0, 0,  0],
# [-1, 1, -1],
# [-1, 0,  0]])

Source code in qmm/extensions/life.py

@cache
def net_life_expectancy_change(
    G: nx.DiGraph,
    type: Literal["birth", "death"] = "birth",
) -> sp.Matrix:
    """Calculate net terms in life expectancy change from press perturbations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        type: Change in birth or death rate ('birth' or 'death')

    Returns:
        sp.Matrix: Net life expectancy change for each component

    References:
        - Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

    Examples:
        ```python
        from qmm import net_life_expectancy_change, load_digraph
        net_life_expectancy_change(load_digraph("snowshoe_rp"), type='birth')
        # Matrix([
        # [-1, 0,  0],
        # [-1, 0, -1],
        # [-1, 0, -1]])

        net_life_expectancy_change(load_digraph("snowshoe_rp"), type='death')
        # Matrix([
        # [ 0, 0,  0],
        # [-1, 1, -1],
        # [-1, 0,  0]])
        ```
    """
    amat = adjoint_matrix(G, form="signed")
    birth = birth_matrix(G, form="signed")
    death = death_matrix(G, form="signed")
    delta_birth = death * amat * sp.Integer(-1)
    delta_death = birth * amat * sp.Integer(-1)
    if type == "birth":
        return delta_birth
    else:
        return delta_death

absolute_life_expectancy_change cached

absolute_life_expectancy_change(G, type='birth')

Calculate absolute terms in life expectancy change from press perturbations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
type	Literal['birth', 'death']	Change in birth or death rate ('birth' or 'death')	'birth'

Returns:

Type	Description
Matrix	sp.Matrix: Absolute life expectancy change for each component

References

• Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

Examples:

from qmm import absolute_life_expectancy_change, load_digraph
absolute_life_expectancy_change(load_digraph("snowshoe_rp"), type='birth')
# Matrix([
# [3, 0, 0],
# [1, 2, 1],
# [1, 2, 1]])

absolute_life_expectancy_change(load_digraph("snowshoe_rp"), type='death')
# Matrix([
# [0, 0, 0],
# [1, 1, 1],
# [1, 2, 2]])

Source code in qmm/extensions/life.py

@cache
def absolute_life_expectancy_change(
    G: nx.DiGraph,
    type: Literal["birth", "death"] = "birth",
) -> sp.Matrix:
    """Calculate absolute terms in life expectancy change from press perturbations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        type: Change in birth or death rate ('birth' or 'death')

    Returns:
        sp.Matrix: Absolute life expectancy change for each component

    References:
        - Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

    Examples:
        ```python
        from qmm import absolute_life_expectancy_change, load_digraph
        absolute_life_expectancy_change(load_digraph("snowshoe_rp"), type='birth')
        # Matrix([
        # [3, 0, 0],
        # [1, 2, 1],
        # [1, 2, 1]])

        absolute_life_expectancy_change(load_digraph("snowshoe_rp"), type='death')
        # Matrix([
        # [0, 0, 0],
        # [1, 1, 1],
        # [1, 2, 2]])
        ```
    """
    sym_amat = adjoint_matrix(G, form="symbolic")
    n = sym_amat.shape[0]
    sym_birth = birth_matrix(G, form="symbolic")
    sym_death = death_matrix(G, form="symbolic")
    sym_delta_birth = sp.expand(sp.Integer(-1) * sym_death * sym_amat)
    sym_delta_death = sp.expand(sp.Integer(-1) * sym_birth * sym_amat)

    def count_symbols(matrix_element):
        return sum(matrix_element.count(sym) for sym in matrix_element.free_symbols)

    def create_abs_matrix(sym_delta_matrix, n):
        return sp.Matrix(n, n, lambda i, j: count_symbols(sym_delta_matrix[i, j]) // n)

    abs_birth = create_abs_matrix(sym_delta_birth, n)
    abs_death = create_abs_matrix(sym_delta_death, n)
    if type == "birth":
        return abs_birth
    else:
        return abs_death

weighted_predictions_life_expectancy cached

weighted_predictions_life_expectancy(
    G, type="birth", as_nan=True, as_abs=False
)

Calculate ratio of net to total change in life expectancy.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
type	Literal['birth', 'death']	Change in birth or death rate ('birth' or 'death')	'birth'
as_nan	bool	Return NaN for undefined ratios	True
as_abs	bool	Return absolute values	False

Returns:

Type	Description
Matrix	sp.Matrix: Net-to-total ratios for life expectancy predictions

References

• Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

Examples:

from qmm import weighted_predictions_life_expectancy, load_digraph
weighted_predictions_life_expectancy(load_digraph("snowshoe_rp"), type='birth', as_abs=False, as_nan=True)
# Matrix([
# [-1/3, nan, nan],
# [  -1,   0,  -1],
# [  -1,   0,  -1]])

weighted_predictions_life_expectancy(load_digraph("snowshoe_rp"), type='birth', as_abs=True, as_nan=False)
# Matrix([
# [1/3, 1, 1],
# [  1, 0, 1],
# [  1, 0, 1]])

weighted_predictions_life_expectancy(load_digraph("snowshoe_rp"), type='death', as_abs=False, as_nan=True)
# Matrix([
# [nan, nan, nan],
# [ -1,   1,  -1],
# [ -1,   0,   0]])

weighted_predictions_life_expectancy(load_digraph("snowshoe_rp"), type='death', as_abs=True, as_nan=False)
# Matrix([
# [1, 1, 1],
# [1, 1, 1],
# [1, 0, 0]])

Source code in qmm/extensions/life.py

@cache
def weighted_predictions_life_expectancy(
    G: nx.DiGraph,
    type: Literal["birth", "death"] = "birth",
    as_nan: bool = True,
    as_abs: bool = False,
) -> sp.Matrix:
    """Calculate ratio of net to total change in life expectancy.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        type: Change in birth or death rate ('birth' or 'death')
        as_nan: Return NaN for undefined ratios
        as_abs: Return absolute values

    Returns:
        sp.Matrix: Net-to-total ratios for life expectancy predictions

    References:
        - Dambacher, J.M., Levins, R., Rossignol, P.A. (2005). Life expectancy change in perturbed communities: Derivation and qualitative analysis. Mathematical Biosciences 197, 1–14.

    Examples:
        ```python
        from qmm import weighted_predictions_life_expectancy, load_digraph
        weighted_predictions_life_expectancy(load_digraph("snowshoe_rp"), type='birth', as_abs=False, as_nan=True)
        # Matrix([
        # [-1/3, nan, nan],
        # [  -1,   0,  -1],
        # [  -1,   0,  -1]])

        weighted_predictions_life_expectancy(load_digraph("snowshoe_rp"), type='birth', as_abs=True, as_nan=False)
        # Matrix([
        # [1/3, 1, 1],
        # [  1, 0, 1],
        # [  1, 0, 1]])

        weighted_predictions_life_expectancy(load_digraph("snowshoe_rp"), type='death', as_abs=False, as_nan=True)
        # Matrix([
        # [nan, nan, nan],
        # [ -1,   1,  -1],
        # [ -1,   0,   0]])

        weighted_predictions_life_expectancy(load_digraph("snowshoe_rp"), type='death', as_abs=True, as_nan=False)
        # Matrix([
        # [1, 1, 1],
        # [1, 1, 1],
        # [1, 0, 0]])
        ```
    """
    if type == "birth":
        net = net_life_expectancy_change(G, type="birth")
        absolute = absolute_life_expectancy_change(G, type="birth")
    elif type == "death":
        net = net_life_expectancy_change(G, type="death")
        absolute = absolute_life_expectancy_change(G, type="death")
    else:
        raise ValueError("type must be either 'birth' or 'death'")
    if as_nan:
        weighted = get_weight(net, absolute)
    else:
        weighted = get_weight(net, absolute, sp.Integer(1))
    if as_abs:
        weighted = sp.Abs(weighted)
    return weighted

qmm.extensions.paths

Analyse causal pathways, cycles and complementary feedback.

get_cycles cached

get_cycles(G)

Find all feedback cycles in the signed digraph.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
Matrix	sp.Matrix: Products of interactions along each cycle

References

• Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
• Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
• Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

Examples:

from qmm import get_cycles, load_digraph
get_cycles(load_digraph("snowshoe_rp"))
# Matrix([
# [           -a_P,P],
# [           -a_R,R],
# [     -a_C,P*a_P,C],
# [     -a_C,R*a_R,C],
# [a_C,P*a_P,R*a_R,C]])

Source code in qmm/extensions/paths.py

@cache
def get_cycles(G: nx.DiGraph) -> sp.Matrix:
    """Find all feedback cycles in the signed digraph.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        sp.Matrix: Products of interactions along each cycle

    References:
        - Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
        - Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
        - Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

    Examples:
        ```python
        from qmm import get_cycles, load_digraph
        get_cycles(load_digraph("snowshoe_rp"))
        # Matrix([
        # [           -a_P,P],
        # [           -a_R,R],
        # [     -a_C,P*a_P,C],
        # [     -a_C,R*a_R,C],
        # [a_C,P*a_P,R*a_R,C]])
        ```
    """
    A = create_matrix(G, form="symbolic")
    nodes = get_nodes(G, "state")
    node_id = {n: i for i, n in enumerate(nodes)}
    cycle_nodes = _sorted_cycles(G)
    C = [c + [c[0]] for c in cycle_nodes]
    cycles = sp.Matrix([sp.prod([A[node_id[c[i + 1]], node_id[c[i]]] for i in range(len(c) - 1)]) for c in C])
    return cycles

cycles_table cached

cycles_table(G)

Find all feedback cycles in the signed digraph.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Table with cycle length, path representation, and sign

References

• Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
• Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
• Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

Examples:

from qmm import cycles_table, load_digraph
cycles_table(load_digraph("snowshoe_rp"))
#    Length                                          Cycle Sign
# 0       1                                P $\multimap$ P    −
# 1       1                                R $\multimap$ R    −
# 2       2                C $\rightarrow$ P $\multimap$ C    −
# 3       2                C $\multimap$ R $\rightarrow$ C    −
# 4       3  C $\multimap$ R $\rightarrow$ P $\multimap$ C    +

Source code in qmm/extensions/paths.py

@cache
def cycles_table(G: nx.DiGraph) -> pd.DataFrame:
    """Find all feedback cycles in the signed digraph.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        pd.DataFrame: Table with cycle length, path representation, and sign

    References:
        - Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
        - Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
        - Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

    Examples:
        ```python
        from qmm import cycles_table, load_digraph
        cycles_table(load_digraph("snowshoe_rp"))
        #    Length                                          Cycle Sign
        # 0       1                                P $\\multimap$ P    −
        # 1       1                                R $\\multimap$ R    −
        # 2       2                C $\\rightarrow$ P $\\multimap$ C    −
        # 3       2                C $\\multimap$ R $\\rightarrow$ C    −
        # 4       3  C $\\multimap$ R $\\rightarrow$ P $\\multimap$ C    +
        ```
    """
    cycle_nodes = _sorted_cycles(G)
    all_cycles = [cycle + [cycle[0]] for cycle in cycle_nodes]
    cycle_signs = [_sign_string(G, path) for path in all_cycles]
    cycles_df = pd.DataFrame(
        {
            "Length": [len(nodes) for nodes in cycle_nodes],
            "Cycle": [_arrows(G, path) for path in all_cycles],
            "Sign": cycle_signs,
        }
    )
    return cycles_df

get_paths cached

get_paths(G, source, target, form='symbolic')

Find all causal pathways between two nodes.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
source	str	Source node	required
target	str	Target node (can be same as source for self-effect)	required
form	Literal['symbolic', 'signed', 'binary']	Type of path products ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

Type	Description
Matrix	sp.Matrix: Products of interactions along each path

References

• Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
• Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
• Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

Examples:

from qmm import get_paths, load_digraph
get_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='symbolic')
# Matrix([
# [a_C,R*a_P,C*b_R,Inp1*c_Out1,P],
# [     -a_C,R*b_R,Inp1*c_Out1,C],
# [     -a_P,C*b_C,Inp1*c_Out1,P],
# [            b_C,Inp1*c_Out1,C]])

get_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='signed')
# Matrix([
# [ 1],
# [-1],
# [-1],
# [ 1]])

get_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='binary')
# Matrix([
# [1],
# [1],
# [1],
# [1]])

Source code in qmm/extensions/paths.py

@cache
def get_paths(
    G: nx.DiGraph,
    source: str,
    target: str,
    form: Literal["symbolic", "signed", "binary"] = "symbolic",
) -> sp.Matrix:
    """Find all causal pathways between two nodes.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        source: Source node
        target: Target node (can be same as source for self-effect)
        form: Type of path products ('symbolic', 'signed', or 'binary')

    Returns:
        sp.Matrix: Products of interactions along each path

    References:
        - Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
        - Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
        - Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

    Examples:
        ```python
        from qmm import get_paths, load_digraph
        get_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='symbolic')
        # Matrix([
        # [a_C,R*a_P,C*b_R,Inp1*c_Out1,P],
        # [     -a_C,R*b_R,Inp1*c_Out1,C],
        # [     -a_P,C*b_C,Inp1*c_Out1,P],
        # [            b_C,Inp1*c_Out1,C]])

        get_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='signed')
        # Matrix([
        # [ 1],
        # [-1],
        # [-1],
        # [ 1]])

        get_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='binary')
        # Matrix([
        # [1],
        # [1],
        # [1],
        # [1]])
        ```
    """
    all_nodes = get_nodes(G, "all")
    if source not in all_nodes or target not in all_nodes:
        raise ValueError("Invalid source or target node")
    if source == target:
        return sp.Matrix([sp.Integer(1)])
    if not nx.has_path(G, source, target):
        return sp.Matrix([sp.Integer(0)])
    path_nodes = list(nx.all_simple_paths(G, source, target))
    if form == "binary":
        return sp.ones(len(path_nodes), 1)
    paths = []
    for p in path_nodes:
        effect = sp.Integer(1)
        for i in range(len(p) - 1):
            u, v = p[i], p[i + 1]
            s = G[u][v].get("sign", 1)
            if form == "signed":
                effect *= sp.Integer(s)
            else:
                effect *= sp.Symbol(f"{_edge_prefix(G, u, v)}_{v},{u}") * s
        paths.append(effect)
    return sp.Matrix(paths)

paths_table cached

paths_table(G, source, target)

Create table of paths between nodes.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
source	str	Source node	required
target	str	Target node	required

Returns:

Type	Description
Optional[DataFrame]	Optional[pd.DataFrame]: DataFrame containing path information or None if no paths exist

References

• Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
• Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
• Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

Examples:

from qmm import paths_table, load_digraph
paths_table(load_digraph("snowshoe_io"), 'Inp1', 'Out1')
#    Length                                                                     Path Sign
# 0       4  Inp1 $\rightarrow$ R $\rightarrow$ C $\rightarrow$ P $\rightarrow$ Out1    +
# 1       3                    Inp1 $\rightarrow$ R $\rightarrow$ C $\multimap$ Out1    −
# 2       3                    Inp1 $\multimap$ C $\rightarrow$ P $\rightarrow$ Out1    −
# 3       2                                      Inp1 $\multimap$ C $\multimap$ Out1    +

Source code in qmm/extensions/paths.py

@cache
def paths_table(G: nx.DiGraph, source: str, target: str) -> Optional[pd.DataFrame]:
    """Create table of paths between nodes.

    Args:
        G (nx.DiGraph): NetworkX DiGraph representing signed digraph model
        source (str): Source node
        target (str): Target node

    Returns:
        Optional[pd.DataFrame]: DataFrame containing path information or None if no paths exist

    References:
        - Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
        - Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
        - Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

    Examples:
        ```python
        from qmm import paths_table, load_digraph
        paths_table(load_digraph("snowshoe_io"), 'Inp1', 'Out1')
        #    Length                                                                     Path Sign
        # 0       4  Inp1 $\\rightarrow$ R $\\rightarrow$ C $\\rightarrow$ P $\\rightarrow$ Out1    +
        # 1       3                    Inp1 $\\rightarrow$ R $\\rightarrow$ C $\\multimap$ Out1    −
        # 2       3                    Inp1 $\\multimap$ C $\\rightarrow$ P $\\rightarrow$ Out1    −
        # 3       2                                      Inp1 $\\multimap$ C $\\multimap$ Out1    +
        ```
    """
    nodes = get_nodes(G, "all")
    if source not in nodes or target not in nodes or source == target:
        raise ValueError("Invalid source or target node")
    if not nx.has_path(G, source, target):
        return None
    paths = list(nx.all_simple_paths(G, source, target))
    paths_df = pd.DataFrame(
        {
            "Length": [len(path) - 1 for path in paths],
            "Path": [_arrows(G, path) for path in paths],
            "Sign": [_sign_string(G, path) for path in paths],
        }
    )
    return paths_df

complementary_feedback cached

complementary_feedback(G, source, target, form='symbolic')

Calculate feedback from state nodes not on paths between source and target.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
source	str	Source node	required
target	str	Target node (can be same as source for self-effect)	required
form	Literal['symbolic', 'signed', 'binary']	Type of feedback ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

Type	Description
Matrix	sp.Matrix: Feedback cycles in complementary subsystem

References

• Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
• Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
• Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

Examples:

from qmm import complementary_feedback, load_digraph
complementary_feedback(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='symbolic')
# Matrix([
# [          -1],
# [      -a_P,P],
# [      -a_R,R],
# [-a_P,P*a_R,R]])

complementary_feedback(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='signed')
# Matrix([
# [-1],
# [-1],
# [-1],
# [-1]])

complementary_feedback(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='binary')
# Matrix([
# [1],
# [1],
# [1],
# [1]])

Source code in qmm/extensions/paths.py

@cache
def complementary_feedback(
    G: nx.DiGraph,
    source: str,
    target: str,
    form: Literal["symbolic", "signed", "binary"] = "symbolic",
) -> sp.Matrix:
    """Calculate feedback from state nodes not on paths between source and target.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        source: Source node
        target: Target node (can be same as source for self-effect)
        form: Type of feedback ('symbolic', 'signed', or 'binary')

    Returns:
        sp.Matrix: Feedback cycles in complementary subsystem

    References:
        - Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
        - Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
        - Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

    Examples:
        ```python
        from qmm import complementary_feedback, load_digraph
        complementary_feedback(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='symbolic')
        # Matrix([
        # [          -1],
        # [      -a_P,P],
        # [      -a_R,R],
        # [-a_P,P*a_R,R]])

        complementary_feedback(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='signed')
        # Matrix([
        # [-1],
        # [-1],
        # [-1],
        # [-1]])

        complementary_feedback(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='binary')
        # Matrix([
        # [1],
        # [1],
        # [1],
        # [1]])
        ```
    """
    state_nodes = get_nodes(G, "state")
    all_nodes = get_nodes(G, "all")
    if source not in all_nodes or target not in all_nodes:
        raise ValueError("Invalid source or target node")
    feedback_funcs = {"symbolic": system_feedback, "signed": net_feedback, "binary": absolute_feedback}
    if form not in feedback_funcs:
        raise ValueError("Invalid form. Choose 'symbolic', 'signed', or 'binary'.")
    if source == target:
        paths = [[source]]
    elif not nx.has_path(G, source, target):
        return sp.Matrix([sp.Integer(0)])
    else:
        paths = list(nx.all_simple_paths(G, source, target))

    feedback = []
    for path in paths:
        subsystem_nodes = [n for n in state_nodes if n not in path]
        if not subsystem_nodes:
            feedback.append(sp.Integer(1) if form == "binary" else sp.Integer(-1))
        else:
            subsystem = G.subgraph(subsystem_nodes).copy()
            feedback.append(feedback_funcs[form](subsystem, level=len(subsystem_nodes))[0])
    return sp.Matrix([sp.expand_mul(f) for f in feedback])

system_paths cached

system_paths(G, source, target, form='symbolic')

Calculate combined effect of paths and complementary feedback.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
source	str	Source node	required
target	str	Target node (can be same as source for self-effect)	required
form	Literal['symbolic', 'signed', 'binary']	Type of computation ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

Type	Description
Matrix	sp.Matrix: Total effects including paths and feedback

References

• Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
• Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
• Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

Examples:

from qmm import system_paths, load_digraph
system_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='symbolic')
# Matrix([
# [ a_C,R*a_P,C*b_R,Inp1*c_Out1,P],
# [-a_C,R*a_P,P*b_R,Inp1*c_Out1,C],
# [-a_P,C*a_R,R*b_C,Inp1*c_Out1,P],
# [ a_P,P*a_R,R*b_C,Inp1*c_Out1,C]])

system_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='signed')
# Matrix([
# [ 1],
# [-1],
# [-1],
# [ 1]])

system_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='binary')
# Matrix([
# [1],
# [1],
# [1],
# [1]])

Source code in qmm/extensions/paths.py

@cache
def system_paths(
    G: nx.DiGraph,
    source: str,
    target: str,
    form: Literal["symbolic", "signed", "binary"] = "symbolic",
) -> sp.Matrix:
    """Calculate combined effect of paths and complementary feedback.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        source: Source node
        target: Target node (can be same as source for self-effect)
        form: Type of computation ('symbolic', 'signed', or 'binary')

    Returns:
        sp.Matrix: Total effects including paths and feedback

    References:
        - Mason, S.J. (1953). Feedback Theory-Some Properties of Signal Flow Graphs. Proceedings of the IRE 41, 1144–1156.
        - Levins, R. (1974). The qualitative analysis of partially specified systems. Annals of the New York Academy of Sciences 231, 123–138.
        - Puccia, C.J., Levins, R. (1985). Qualitative Modeling of Complex Systems: An Introduction to Loop Analysis and Time Averaging. Harvard University Press.

    Examples:
        ```python
        from qmm import system_paths, load_digraph
        system_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='symbolic')
        # Matrix([
        # [ a_C,R*a_P,C*b_R,Inp1*c_Out1,P],
        # [-a_C,R*a_P,P*b_R,Inp1*c_Out1,C],
        # [-a_P,C*a_R,R*b_C,Inp1*c_Out1,P],
        # [ a_P,P*a_R,R*b_C,Inp1*c_Out1,C]])

        system_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='signed')
        # Matrix([
        # [ 1],
        # [-1],
        # [-1],
        # [ 1]])

        system_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1', form='binary')
        # Matrix([
        # [1],
        # [1],
        # [1],
        # [1]])
        ```
    """
    _check_direct_io_edges(G)
    path = get_paths(G, source, target, form=form)
    feedback = complementary_feedback(G, source, target, form=form)
    if form == "binary":
        effect = path.multiply_elementwise(feedback)
    else:
        effect = path.multiply_elementwise(feedback) / sp.Integer(-1)
    return sp.Matrix([sp.expand_mul(e) for e in effect])

weighted_paths cached

weighted_paths(G, source, target)

Calculate ratio of net to total path effects.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
source	str	Source node	required
target	str	Target node	required

Returns:

Type	Description
Matrix	sp.Matrix: Net-to-total ratios for path predictions

Examples:

from qmm import weighted_paths, load_digraph
weighted_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1')
# Matrix([
# [ 1],
# [-1],
# [-1],
# [ 1]])

Source code in qmm/extensions/paths.py

@cache
def weighted_paths(G: nx.DiGraph, source: str, target: str) -> sp.Matrix:
    """Calculate ratio of net to total path effects.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        source: Source node
        target: Target node

    Returns:
        sp.Matrix: Net-to-total ratios for path predictions

    Examples:
        ```python
        from qmm import weighted_paths, load_digraph
        weighted_paths(load_digraph("snowshoe_io"), 'Inp1', 'Out1')
        # Matrix([
        # [ 1],
        # [-1],
        # [-1],
        # [ 1]])
        ```
    """
    _check_direct_io_edges(G)
    state_nodes = get_nodes(G, "state")
    all_nodes = get_nodes(G, "all")
    if source not in all_nodes or target not in all_nodes or source == target:
        raise ValueError("Invalid source or target node")
    path_nodes = list(nx.all_simple_paths(G, source, target))
    wgt_effects = []
    for path in path_nodes:
        subsystem_nodes = [n for n in state_nodes if n not in path]
        if not subsystem_nodes:
            feedback = sp.Integer(-1)
        else:
            subsystem = G.subgraph(subsystem_nodes).copy()
            feedback = weighted_feedback(subsystem, level=len(subsystem_nodes))
            if feedback[0] == sp.nan:
                feedback = sp.Integer(0)
        sign = 1
        for i in range(len(path) - 1):
            sign *= G[path[i]][path[i + 1]].get('sign', 1)
        wgt_effect = sp.Integer(-1) * sign * feedback
        wgt_effects.append(wgt_effect)
    return sp.Matrix(wgt_effects)

path_metrics cached

path_metrics(G, source, target)

Calculate comprehensive metrics for paths between nodes.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
source	str	Source node	required
target	str	Target node	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Metrics including path length, sign, and feedback

Examples:

from qmm import path_metrics, load_digraph
path_metrics(load_digraph("snowshoe_io"), 'Inp1', 'Out1')
#    Length                 Path Path sign Complementary subsystem Net feedback Absolute feedback Positive feedback Negative feedback Weighted feedback Weighted path System path
# 0       4  Inp1, R, C, P, Out1         +                    None           -1                 1                 0                 1                -1             1           1
# 1       3     Inp1, R, C, Out1         −                       P           -1                 1                 0                 1                -1            -1          -1
# 2       3     Inp1, C, P, Out1         −                       R           -1                 1                 0                 1                -1            -1          -1
# 3       2        Inp1, C, Out1         +                    R, P           -1                 1                 0                 1                -1             1           1

Source code in qmm/extensions/paths.py

@cache
def path_metrics(G: nx.DiGraph, source: str, target: str) -> pd.DataFrame:
    """Calculate comprehensive metrics for paths between nodes.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        source: Source node
        target: Target node

    Returns:
        pd.DataFrame: Metrics including path length, sign, and feedback

    Examples:
        ```python
        from qmm import path_metrics, load_digraph
        path_metrics(load_digraph("snowshoe_io"), 'Inp1', 'Out1')
        #    Length                 Path Path sign Complementary subsystem Net feedback Absolute feedback Positive feedback Negative feedback Weighted feedback Weighted path System path
        # 0       4  Inp1, R, C, P, Out1         +                    None           -1                 1                 0                 1                -1             1           1
        # 1       3     Inp1, R, C, Out1         −                       P           -1                 1                 0                 1                -1            -1          -1
        # 2       3     Inp1, C, P, Out1         −                       R           -1                 1                 0                 1                -1            -1          -1
        # 3       2        Inp1, C, Out1         +                    R, P           -1                 1                 0                 1                -1             1           1
        ```
    """
    _check_direct_io_edges(G)
    state_nodes = get_nodes(G, "state")
    all_nodes = get_nodes(G, "all")
    if source not in all_nodes or target not in all_nodes or source == target:
        raise ValueError("Invalid source or target node")
    if not nx.has_path(G, source, target):
        return pd.DataFrame()
    paths = list(nx.all_simple_paths(G, source, target))
    subsystem_nodes = [[n for n in state_nodes if n not in set(path)] for path in paths]
    net_fb = complementary_feedback(G, source=source, target=target, form="signed")
    absolute_fb = complementary_feedback(G, source=source, target=target, form="binary")
    path_signs = get_paths(G, source=source, target=target, form="signed")
    weighted_fb = get_weight(net_fb, absolute_fb, sp.Integer(0))
    positive_fb = get_positive(net_fb, absolute_fb)
    negative_fb = get_negative(net_fb, absolute_fb)
    weighted_path = weighted_paths(G, source, target)
    system_path = system_paths(G, source, target, form="signed")
    n = len(paths)
    paths_df = pd.DataFrame(
        {
            "Length": [len(path) - 1 for path in paths],
            "Path": [", ".join(str(x) for x in path) for path in paths],
            "Path sign": ["+" if sign == 1 else "−" for sign in path_signs[:n]],
            "Complementary subsystem": [
                ", ".join(str(x) for x in nodes) if nodes else None for nodes in subsystem_nodes
            ],
            "Net feedback": [net_fb[i] for i in range(n)],
            "Absolute feedback": [absolute_fb[i] for i in range(n)],
            "Positive feedback": [positive_fb[i] for i in range(n)],
            "Negative feedback": [negative_fb[i] for i in range(n)],
            "Weighted feedback": [weighted_fb[i] for i in range(n)],
            "Weighted path": [weighted_path[i] for i in range(n)],
            "System path": [system_path[i] for i in range(n)],
        }
    )
    return paths_df

qmm.extensions.effects

Analyse cumulative effects from perturbation scenarios with multiple-inputs and multiple-outputs.

define_input_output

define_input_output(G, remove_disconnected=True)

Classify nodes as state, input or output from topology (any pre-set category is overwritten).

Sources (and source-chains) become inputs, sinks (and sink-chains) outputs; a self-loop or feedback cycle keeps a node as state.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
remove_disconnected	bool	Remove all but the largest weakly-connected component (warns about dropped nodes)	True

Returns:

Type	Description
DiGraph	nx.DiGraph: Model with input, state and output classification

Examples:

from qmm import define_input_output, load_digraph
G = load_digraph("snowshoe_io")
G_io = define_input_output(G)
(G_io.nodes['Inp1']['category'], G_io.nodes['Out1']['category'])
# ('input', 'output')

Source code in qmm/extensions/effects.py

def define_input_output(G: nx.DiGraph, remove_disconnected: bool = True) -> nx.DiGraph:
    """Classify nodes as state, input or output from topology (any pre-set category is overwritten).

    Sources (and source-chains) become inputs, sinks (and sink-chains) outputs; a
    self-loop or feedback cycle keeps a node as state.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        remove_disconnected: Remove all but the largest weakly-connected
            component (warns about dropped nodes)

    Returns:
        nx.DiGraph: Model with input, state and output classification

    Examples:
        ```python
        from qmm import define_input_output, load_digraph
        G = load_digraph("snowshoe_io")
        G_io = define_input_output(G)
        (G_io.nodes['Inp1']['category'], G_io.nodes['Out1']['category'])
        # ('input', 'output')
        ```
    """
    if not isinstance(G, nx.DiGraph):
        raise TypeError("Input must be a networkx.DiGraph.")
    _check_signs(G)
    G_def = G.copy()
    if remove_disconnected:
        components = list(nx.connected_components(G_def.to_undirected()))
        if len(components) > 1:
            largest = max(components, key=lambda c: (len(c), sorted(c)))
            dropped = sorted(n for c in components if c != largest for n in c)
            warnings.warn(
                f"define_input_output: dropping {len(dropped)} node(s) in "
                f"{len(components) - 1} smaller disconnected component(s): {dropped}"
            )
            G_def.remove_nodes_from(dropped)
    nx.set_node_attributes(G_def, "state", "category")

    # Inputs then outputs, each a fixpoint (order-independent); self-loop/feedback nodes stay state.
    def classify(role, here, there):
        changed = True
        while changed:
            changed = False
            for node in G_def.nodes():
                if G_def.nodes[node]["category"] != "state" or G_def.has_edge(node, node) or not list(there(node)):
                    continue
                anchor = list(here(node))
                if not anchor or all(G_def.nodes[n]["category"] == role for n in anchor):
                    G_def.nodes[node]["category"] = role
                    changed = True

    classify("input", G_def.predecessors, G_def.successors)
    classify("output", G_def.successors, G_def.predecessors)

    _check_direct_io_edges(G_def)
    nx.freeze(G_def)
    return G_def

direct_effects cached

direct_effects(G, form='net')

Calculate direct effects from the signed digraph structure.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
form	Literal['net', 'absolute', 'positive', 'negative']	Type of direct effects ('net', 'absolute', 'positive', 'negative')	'net'

Returns:

Type	Description
Matrix	sp.Matrix: Direct effects for state/input columns and state/output rows

Source code in qmm/extensions/effects.py

@cache
def direct_effects(
    G: nx.DiGraph,
    form: Literal["net", "absolute", "positive", "negative"] = "net",
) -> sp.Matrix:
    """Calculate direct effects from the signed digraph structure.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        form: Type of direct effects ('net', 'absolute', 'positive', 'negative')

    Returns:
        sp.Matrix: Direct effects for state/input columns and state/output rows
    """
    if form not in ("net", "absolute", "positive", "negative"):
        raise ValueError("Invalid form. Choose 'net', 'absolute', 'positive', 'negative'.")

    def block(f):
        m = {t: create_matrix(G, form=f, matrix_type=t) for t in "ABCD"}
        return sp.BlockMatrix([[m["A"], m["B"]], [m["C"], m["D"]]]).as_explicit()

    if form == "net":
        return block("signed")
    if form == "absolute":
        return block("binary")
    signed, binary = block("signed"), block("binary")
    return get_positive(signed, binary) if form == "positive" else get_negative(signed, binary)

cumulative_effects cached

cumulative_effects(G, form='symbolic')

Calculate cumulative effects to multiple inputs using state-space representation.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
form	Literal['symbolic', 'signed', 'binary']	Type of computation ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

Type	Description
Matrix	sp.Matrix: Cumulative effects on state variables and outputs

Examples:

from qmm import load_digraph, cumulative_effects
cumulative_effects(load_digraph("snowshoe_io"), form='symbolic')[3, 3]
# a_C,R*a_P,C*b_R,Inp1*c_Out1,P - a_C,R*a_P,P*b_R,Inp1*c_Out1,C - a_P,C*a_R,R*b_C,Inp1*c_Out1,P + a_P,P*a_R,R*b_C,Inp1*c_Out1,C

cumulative_effects(load_digraph("snowshoe_io"), form='signed')
# Matrix([
# [1, -1,  1, 2, -1],
# [1,  1, -1, 0,  1],
# [1,  1,  1, 0, -1],
# [0,  0,  2, 0, -2],
# [1,  1, -1, 0,  1]])

cumulative_effects(load_digraph("snowshoe_io"), form='binary')
# Matrix([
# [1, 1, 1, 2, 1],
# [1, 1, 1, 2, 1],
# [1, 1, 1, 2, 1],
# [2, 2, 2, 4, 2],
# [1, 1, 1, 2, 1]])

Source code in qmm/extensions/effects.py

@cache
def cumulative_effects(
    G: nx.DiGraph,
    form: Literal["symbolic", "signed", "binary"] = "symbolic",
) -> sp.Matrix:
    """Calculate cumulative effects to multiple inputs using state-space representation.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        form: Type of computation ('symbolic', 'signed', or 'binary')

    Returns:
        sp.Matrix: Cumulative effects on state variables and outputs

    Examples:
        ```python
        from qmm import load_digraph, cumulative_effects
        cumulative_effects(load_digraph("snowshoe_io"), form='symbolic')[3, 3]
        # a_C,R*a_P,C*b_R,Inp1*c_Out1,P - a_C,R*a_P,P*b_R,Inp1*c_Out1,C - a_P,C*a_R,R*b_C,Inp1*c_Out1,P + a_P,P*a_R,R*b_C,Inp1*c_Out1,C

        cumulative_effects(load_digraph("snowshoe_io"), form='signed')
        # Matrix([
        # [1, -1,  1, 2, -1],
        # [1,  1, -1, 0,  1],
        # [1,  1,  1, 0, -1],
        # [0,  0,  2, 0, -2],
        # [1,  1, -1, 0,  1]])

        cumulative_effects(load_digraph("snowshoe_io"), form='binary')
        # Matrix([
        # [1, 1, 1, 2, 1],
        # [1, 1, 1, 2, 1],
        # [1, 1, 1, 2, 1],
        # [2, 2, 2, 4, 2],
        # [1, 1, 1, 2, 1]])
        ```
    """
    if form not in ("symbolic", "signed", "binary"):
        raise ValueError("Invalid form. Choose 'symbolic', 'signed', 'binary'.")
    B = create_matrix(G, form=form, matrix_type="B")
    C = create_matrix(G, form=form, matrix_type="C")
    D = create_matrix(G, form=form, matrix_type="D")
    effects = absolute_feedback_matrix(G) if form == "binary" else adjoint_matrix(G, form=form)
    cemat = sp.BlockMatrix([[effects, effects * B], [C * effects, C * effects * B + D]]).as_explicit()
    if form != "symbolic":
        cemat = cemat.subs({sym: 1 for sym in cemat.free_symbols})
    return sp.expand(cemat)

net_effects cached

net_effects(G)

Calculate net cumulative effects from multiple inputs.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
Matrix	sp.Matrix: Net effects on state variables and outputs

Examples:

from qmm import load_digraph, net_effects
net_effects(load_digraph("snowshoe_io"))
# Matrix([
# [1, -1,  1, 2, -1],
# [1,  1, -1, 0,  1],
# [1,  1,  1, 0, -1],
# [0,  0,  2, 0, -2],
# [1,  1, -1, 0,  1]])

Source code in qmm/extensions/effects.py

@cache
def net_effects(G: nx.DiGraph) -> sp.Matrix:
    """Calculate net cumulative effects from multiple inputs.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        sp.Matrix: Net effects on state variables and outputs

    Examples:
        ```python
        from qmm import load_digraph, net_effects
        net_effects(load_digraph("snowshoe_io"))
        # Matrix([
        # [1, -1,  1, 2, -1],
        # [1,  1, -1, 0,  1],
        # [1,  1,  1, 0, -1],
        # [0,  0,  2, 0, -2],
        # [1,  1, -1, 0,  1]])
        ```
    """
    return cumulative_effects(G, form="signed")

absolute_effects cached

absolute_effects(G)

Calculate absolute effects from multiple inputs.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
Matrix	sp.Matrix: Total effects on state variables and outputs

Examples:

from qmm import load_digraph, absolute_effects
absolute_effects(load_digraph("snowshoe_io"))
# Matrix([
# [1, 1, 1, 2, 1],
# [1, 1, 1, 2, 1],
# [1, 1, 1, 2, 1],
# [2, 2, 2, 4, 2],
# [1, 1, 1, 2, 1]])

Source code in qmm/extensions/effects.py

@cache
def absolute_effects(G: nx.DiGraph) -> sp.Matrix:
    """Calculate absolute effects from multiple inputs.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        sp.Matrix: Total effects on state variables and outputs

    Examples:
        ```python
        from qmm import load_digraph, absolute_effects
        absolute_effects(load_digraph("snowshoe_io"))
        # Matrix([
        # [1, 1, 1, 2, 1],
        # [1, 1, 1, 2, 1],
        # [1, 1, 1, 2, 1],
        # [2, 2, 2, 4, 2],
        # [1, 1, 1, 2, 1]])
        ```
    """
    return cumulative_effects(G, form="binary")

weighted_effects cached

weighted_effects(G)

Calculate ratio of net to total terms for predicting cumulative effects.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required

Returns:

Type	Description
Matrix	sp.Matrix: Ratio of net to total effects

Examples:

from qmm import load_digraph, weighted_effects
weighted_effects(load_digraph("snowshoe_io"))
# Matrix([
# [1, -1,  1, 1, -1],
# [1,  1, -1, 0,  1],
# [1,  1,  1, 0, -1],
# [0,  0,  1, 0, -1],
# [1,  1, -1, 0,  1]])

Source code in qmm/extensions/effects.py

@cache
def weighted_effects(G: nx.DiGraph) -> sp.Matrix:
    """Calculate ratio of net to total terms for predicting cumulative effects.

    Args:
        G: NetworkX DiGraph representing signed digraph model

    Returns:
        sp.Matrix: Ratio of net to total effects

    Examples:
        ```python
        from qmm import load_digraph, weighted_effects
        weighted_effects(load_digraph("snowshoe_io"))
        # Matrix([
        # [1, -1,  1, 1, -1],
        # [1,  1, -1, 0,  1],
        # [1,  1,  1, 0, -1],
        # [0,  0,  1, 0, -1],
        # [1,  1, -1, 0,  1]])
        ```
    """
    return get_weight(net_effects(G), absolute_effects(G))

sign_determinacy_effects cached

sign_determinacy_effects(G, method='average')

Calculate probability of correct sign prediction for cumulative effects.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
method	Literal['average', '95_bound']	Method for computing determinacy ('average', '95_bound')	'average'

Returns:

Type	Description
Matrix	sp.Matrix: Sign determinacy probabilities for effects

Examples:

from qmm import load_digraph, sign_determinacy_effects
sign_determinacy_effects(load_digraph("snowshoe_io"), method='average')
# Matrix([
# [  1,  -1,  1,   1, -1],
# [  1,   1, -1, 1/2,  1],
# [  1,   1,  1, 1/2, -1],
# [1/2, 1/2,  1, 1/2, -1],
# [  1,   1, -1, 1/2,  1]])

sign_determinacy_effects(load_digraph("snowshoe_io"), method='95_bound')
# Matrix([
# [  1,  -1,  1,   1, -1],
# [  1,   1, -1, 1/2,  1],
# [  1,   1,  1, 1/2, -1],
# [1/2, 1/2,  1, 1/2, -1],
# [  1,   1, -1, 1/2,  1]])

Source code in qmm/extensions/effects.py

@cache
def sign_determinacy_effects(
    G: nx.DiGraph,
    method: Literal["average", "95_bound"] = "average",
) -> sp.Matrix:
    """Calculate probability of correct sign prediction for cumulative effects.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        method: Method for computing determinacy ('average', '95_bound')

    Returns:
        sp.Matrix: Sign determinacy probabilities for effects

    Examples:
        ```python
        from qmm import load_digraph, sign_determinacy_effects
        sign_determinacy_effects(load_digraph("snowshoe_io"), method='average')
        # Matrix([
        # [  1,  -1,  1,   1, -1],
        # [  1,   1, -1, 1/2,  1],
        # [  1,   1,  1, 1/2, -1],
        # [1/2, 1/2,  1, 1/2, -1],
        # [  1,   1, -1, 1/2,  1]])

        sign_determinacy_effects(load_digraph("snowshoe_io"), method='95_bound')
        # Matrix([
        # [  1,  -1,  1,   1, -1],
        # [  1,   1, -1, 1/2,  1],
        # [  1,   1,  1, 1/2, -1],
        # [1/2, 1/2,  1, 1/2, -1],
        # [  1,   1, -1, 1/2,  1]])
        ```
    """
    absolute = absolute_effects(G)
    return sign_determinacy(weighted_effects(G), absolute, method=method)

get_simulations cached

get_simulations(
    G,
    n_sim=10000,
    dist="uniform",
    seed=42,
    perturb=None,
    observe=None,
    presample=None,
    return_samples=False,
    average_uncertain=False,
)

Calculate average proportion of positive and negative effects from stable numerical simulations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
n_sim	int	Number of simulations	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution for sampling	'uniform'
seed	int	Random seed	42
perturb	Optional[Tuple[str, int]]	Optional perturbations (node, sign)	None
observe	Optional[Tuple[Tuple[str, int], ...]]	Optional observations (node, sign)	None
presample	Optional[Callable[[Tuple[Symbol, ...]], Dict[Symbol, Any]]]	Optional callable that receives the tuple of free symbols and returns a mapping of symbol substitutions to apply before sampling.	None
return_samples	bool	If True, include dict mapping symbol names to arrays of sampled values	False
average_uncertain	bool	If True, sample edges marked dashes=True in/out each draw (structure averaging)	False

Returns:

Type	Description
Dict[str, Any]	Dict containing effects, valid_sims, all_nodes, tmat, prop_stable, attempts, and optionally samples.

Examples:

from qmm import get_simulations, load_digraph
result = get_simulations(load_digraph("snowshoe_io"), n_sim=1000, perturb=('Inp1', 1))
result['effects'][0]
# array([ 1.29385476,  2.55918625,  3.12917778, -1.74767706,  2.48217995])

Source code in qmm/extensions/effects.py

@cache
def get_simulations(
    G: nx.DiGraph,
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
    perturb: Optional[Tuple[str, int]] = None,
    observe: Optional[Tuple[Tuple[str, int], ...]] = None,
    presample: Optional[Callable[[Tuple[sp.Symbol, ...]], Dict[sp.Symbol, Any]]] = None,
    return_samples: bool = False,
    average_uncertain: bool = False,
) -> Dict[str, Any]:
    """Calculate average proportion of positive and negative effects from stable numerical simulations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        n_sim: Number of simulations
        dist: Distribution for sampling
        seed: Random seed
        perturb: Optional perturbations (node, sign)
        observe: Optional observations (node, sign)
        presample: Optional callable that receives the tuple of free symbols and
            returns a mapping of symbol substitutions to apply before sampling.
        return_samples: If True, include dict mapping symbol names to arrays of sampled values
        average_uncertain: If True, sample edges marked dashes=True in/out each draw (structure averaging)

    Returns:
        Dict containing effects, valid_sims, all_nodes, tmat, prop_stable, attempts, and optionally samples.

    Examples:
        ```python
        from qmm import get_simulations, load_digraph
        result = get_simulations(load_digraph("snowshoe_io"), n_sim=1000, perturb=('Inp1', 1))
        result['effects'][0]
        # array([ 1.29385476,  2.55918625,  3.12917778, -1.74767706,  2.48217995])
        ```
    """

    rng = np.random.RandomState(seed)

    A_sym, B_sym, C_sym, D_sym = (create_matrix(G, form="symbolic", matrix_type=m) for m in "ABCD")
    symbols_all = {s for m in (A_sym, B_sym, C_sym, D_sym) if m for s in m.free_symbols}
    symbols_all = tuple(sorted(symbols_all, key=str))
    fixed_subs = {}

    if presample and symbols_all and (subs := presample(symbols_all)):
        fixed_subs = {sym: subs[sym] for sym in symbols_all if sym in subs}
        A_sym, B_sym, C_sym, D_sym = (m.subs(subs) if m else None for m in (A_sym, B_sym, C_sym, D_sym))
        symbols_sampled = tuple(sorted({s for m in (A_sym, B_sym, C_sym, D_sym) if m for s in m.free_symbols}, key=str))
    else:
        symbols_sampled = symbols_all

    state, inputs, outputs = (get_nodes(G, t) for t in ("state", "input", "output"))
    all_nodes = state + inputs + outputs
    n_x, n_u, n_y = len(state), len(inputs), len(outputs)

    response_idx = {node: i for i, node in enumerate(state + outputs)}
    if observe:
        unknown = [n for n, _ in observe if n not in response_idx]
        if unknown:
            raise ValueError(f"Unknown observation node(s): {unknown}. Valid response nodes: {list(response_idx)}")
    perturb_nodes = state + inputs
    if perturb and perturb[0] not in perturb_nodes:
        raise ValueError(f"Perturbation node '{perturb[0]}' not found.")
    p_idx, p_sign = (perturb_nodes.index(perturb[0]), perturb[1]) if perturb else (None, 1)

    tmat = sp.matrix2numpy(absolute_effects(G)).astype(int)

    A_fn = sp.lambdify(symbols_sampled, A_sym)
    B_fn = sp.lambdify(symbols_sampled, B_sym) if n_u and B_sym else None
    C_fn = sp.lambdify(symbols_sampled, C_sym) if n_y and C_sym else None
    D_fn = sp.lambdify(symbols_sampled, D_sym) if D_sym and D_sym.shape != (0, 0) else None

    def compute_sample(values):
        A = A_fn(*values)
        if not np.all(np.real(np.linalg.eigvals(A)) < 0):
            return None
        try:
            inv_A = np.linalg.inv(-A)
        except np.linalg.LinAlgError:
            return None
        B = B_fn(*values) if B_fn else np.zeros((n_x, 0))
        C = C_fn(*values) if C_fn else np.zeros((0, n_x))
        D = D_fn(*values) if D_fn else np.zeros((n_y, n_u))
        E = np.block([[inv_A, inv_A @ B], [C @ inv_A, C @ inv_A @ B + D]]) if n_x else np.zeros((n_y, n_u))
        effect = E[:, p_idx] * p_sign if p_idx is not None else E
        return effect

    def is_valid(effect, tmat_ref):
        if not observe or tmat_ref is None:
            return True
        for node, obs in observe:
            idx = response_idx[node]
            expected = tmat_ref[idx, p_idx] != 0
            if (expected and (obs == 0 or np.sign(effect[idx]) != obs)) or (not expected and obs != 0):
                return False
        return True

    uncertain = [(u, v, symbols_sampled.index(sp.Symbol(f"{_edge_prefix(G, u, v)}_{v},{u}")))
                 for u, v, d in G.edges(data=True) if d.get("dashes")] if average_uncertain else []
    base_cls, checked = {n: d.get("category", "state") for n, d in G.nodes(data=True)}, set()
    effects, valid_sims, samples = [], [], []
    attempts, max_attempts = 0, n_sim * 100
    while len(effects) < n_sim and attempts < max_attempts:
        attempts += 1
        values = _random_sampler(dist, len(symbols_sampled), rng)
        if uncertain:
            keep = rng.uniform(size=len(uncertain)) < rng.uniform()
            if (k := tuple(keep)) not in checked:
                variant = nx.DiGraph(G)
                variant.remove_edges_from([(u, v) for ke, (u, v, _) in zip(keep, uncertain) if not ke])
                if {n: d.get("category", "state") for n, d in define_input_output(variant, remove_disconnected=False).nodes(data=True)} != base_cls:
                    raise ValueError("Excluding an uncertain edge re-classifies a node; structure averaging halted.")
                checked.add(k)
            values[[i for ke, (_, _, i) in zip(keep, uncertain) if not ke]] = 0.0
        if (effect := compute_sample(values)) is None:
            continue
        effects.append(effect)
        valid_sims.append(is_valid(effect, tmat))
        if return_samples:
            samples.append(values)

    if len(effects) < n_sim:
        raise RuntimeError(f"Maximum iterations reached. Stable proportion: {len(effects) / max_attempts:.4f}")

    result = {
        "effects": effects,
        "valid_sims": valid_sims,
        "all_nodes": all_nodes,
        "tmat": tmat,
        "prop_stable": len(effects) / attempts,
        "attempts": attempts,
        "n_stable": len(effects),
    }
    if return_samples:
        n_samples = len(samples)
        sampled_index = {sym: i for i, sym in enumerate(symbols_sampled)}
        result_samples = {}
        for sym in symbols_all:
            if sym in sampled_index:
                idx = sampled_index[sym]
                result_samples[str(sym)] = np.array([s[idx] for s in samples])
            elif sym in fixed_subs:
                result_samples[str(sym)] = np.full(n_samples, fixed_subs[sym])
        result["samples"] = result_samples
    return result

simulation_effects

simulation_effects(
    G,
    n_sim=10000,
    dist="uniform",
    seed=42,
    positive_only=False,
    presample=None,
    average_uncertain=False,
)

Performs numerical simulations of cumulative effects using random interaction strengths.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
n_sim	int	Number of simulations	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution for sampling ("uniform", "weak", "moderate", "strong")	'uniform'
seed	int	Random seed	42
positive_only	bool	Return just the proportion of positive responses instead of sign-dominant proportions	False
presample	Optional[Callable[[Tuple[Symbol, ...]], Dict[Symbol, Any]]]	Optional callable passed through to get_simulations	None
average_uncertain	bool	Passed through to get_simulations (structure averaging over uncertain links)	False

Returns:

Type	Description
Matrix	SymPy Matrix containing simulation results

Examples:

from qmm import load_digraph, simulation_effects
simulation_effects(load_digraph("snowshoe_io"), n_sim=1000)
# Matrix([
# [   1.0,   -1.0,  1.0,    1.0, -1.0],
# [   1.0,    1.0, -1.0, -0.513,  1.0],
# [   1.0,    1.0,  1.0, -0.513, -1.0],
# [-0.517, -0.517,  1.0,  0.526, -1.0],
# [   1.0,    1.0, -1.0, -0.513,  1.0]])

simulation_effects(load_digraph("snowshoe_io"), n_sim=1000, positive_only=True)
# Matrix([
# [  1.0,   0.0, 1.0,   1.0, 0.0],
# [  1.0,   1.0, 0.0, 0.487, 1.0],
# [  1.0,   1.0, 1.0, 0.487, 0.0],
# [0.483, 0.483, 1.0, 0.526, 0.0],
# [  1.0,   1.0, 0.0, 0.487, 1.0]])

Source code in qmm/extensions/effects.py

def simulation_effects(
    G: nx.DiGraph,
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
    positive_only: bool = False,
    presample: Optional[Callable[[Tuple[sp.Symbol, ...]], Dict[sp.Symbol, Any]]] = None,
    average_uncertain: bool = False,
) -> sp.Matrix:
    """Performs numerical simulations of cumulative effects using random interaction strengths.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        n_sim: Number of simulations
        dist: Distribution for sampling ("uniform", "weak", "moderate", "strong")
        seed: Random seed
        positive_only: Return just the proportion of positive responses instead of sign-dominant proportions
        presample: Optional callable passed through to get_simulations
        average_uncertain: Passed through to get_simulations (structure averaging over uncertain links)

    Returns:
        SymPy Matrix containing simulation results

    Examples:
        ```python
        from qmm import load_digraph, simulation_effects
        simulation_effects(load_digraph("snowshoe_io"), n_sim=1000)
        # Matrix([
        # [   1.0,   -1.0,  1.0,    1.0, -1.0],
        # [   1.0,    1.0, -1.0, -0.513,  1.0],
        # [   1.0,    1.0,  1.0, -0.513, -1.0],
        # [-0.517, -0.517,  1.0,  0.526, -1.0],
        # [   1.0,    1.0, -1.0, -0.513,  1.0]])

        simulation_effects(load_digraph("snowshoe_io"), n_sim=1000, positive_only=True)
        # Matrix([
        # [  1.0,   0.0, 1.0,   1.0, 0.0],
        # [  1.0,   1.0, 0.0, 0.487, 1.0],
        # [  1.0,   1.0, 1.0, 0.487, 0.0],
        # [0.483, 0.483, 1.0, 0.526, 0.0],
        # [  1.0,   1.0, 0.0, 0.487, 1.0]])
        ```
    """
    sims = get_simulations(G, n_sim, dist, seed, presample=presample, average_uncertain=average_uncertain)
    tmat = sims["tmat"]
    n_rows, n_cols = tmat.shape

    positive, negative = _sign_counts(sims["effects"])

    smat = positive / n_sim if positive_only else np.where(
        negative > positive, -negative / n_sim, positive / n_sim
    )
    smat = [[sp.nan if not tmat[i, j] else smat[i, j] for j in range(n_cols)] for i in range(n_rows)]
    return sp.Matrix(smat)

simulations_table

simulations_table(
    G,
    perturb,
    observe="",
    n_sim=10000,
    dist="uniform",
    seed=42,
    combinations=True,
    presample=None,
)

Summarise simulation effects across model variants for each response node.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
perturb	str	Node and sign to perturb (perturbation string)	required
observe	str	Observation string (node:sign, comma-separated allowed) to filter simulations	''
n_sim	int	Number of simulations	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution for sampling	'uniform'
seed	int	Random seed	42
combinations	bool	If True, evaluate every combination of dashed edges	True
presample	Optional[Callable[[Tuple[Symbol, ...]], Dict[Symbol, Any]]]	Optional callable passed through to get_simulations	None

Returns:

Type	Description
DataFrame	pd.DataFrame: Table of counts for negative, no effect, and positive responses

Examples:

from qmm import load_digraph, simulations_table
simulations_table(load_digraph("snowshoe_io"), perturb='Inp1:+', n_sim=1000)
#    model effect_on  negative  no_effect  positive  valid_sims  stable_sims  attempts
# 0      1         R         0          0      1000        1000         1000      1000
# 1      1         C       513          0       487        1000         1000      1000
# 2      1         P       513          0       487        1000         1000      1000
# 3      1      Out1       474          0       526        1000         1000      1000
# 4      1      Out2       513          0       487        1000         1000      1000

Source code in qmm/extensions/effects.py

def simulations_table(
    G: nx.DiGraph,
    perturb: str,
    observe: str = "",
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
    combinations: bool = True,
    presample: Optional[Callable[[Tuple[sp.Symbol, ...]], Dict[sp.Symbol, Any]]] = None,
) -> pd.DataFrame:
    """Summarise simulation effects across model variants for each response node.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        perturb: Node and sign to perturb (perturbation string)
        observe: Observation string (node:sign, comma-separated allowed) to filter simulations
        n_sim: Number of simulations
        dist: Distribution for sampling
        seed: Random seed
        combinations: If True, evaluate every combination of dashed edges
        presample: Optional callable passed through to get_simulations

    Returns:
        pd.DataFrame: Table of counts for negative, no effect, and positive responses

    Examples:
        ```python
        from qmm import load_digraph, simulations_table
        simulations_table(load_digraph("snowshoe_io"), perturb='Inp1:+', n_sim=1000)
        #    model effect_on  negative  no_effect  positive  valid_sims  stable_sims  attempts
        # 0      1         R         0          0      1000        1000         1000      1000
        # 1      1         C       513          0       487        1000         1000      1000
        # 2      1         P       513          0       487        1000         1000      1000
        # 3      1      Out1       474          0       526        1000         1000      1000
        # 4      1      Out2       513          0       487        1000         1000      1000
        ```
    """
    variants = get_dashed_alternatives(G, combinations=combinations)
    observations = _parse_observations(observe) if observe else None
    rows = []

    for model_idx, g in enumerate(variants, start=1):
        graph, pert = _parse_perturbations(g, perturb)
        sims = get_simulations(
            graph,
            n_sim=n_sim,
            dist=dist,
            seed=seed,
            perturb=pert,
            observe=observations,
            presample=presample,
        )

        response_nodes = get_nodes(g, "state") + get_nodes(g, "output")
        if not response_nodes:
            continue
        node_count = len(response_nodes)
        p_idx = sims["all_nodes"].index(pert[0])
        tmat = sims["tmat"][:node_count, :]

        valid_effects = [effect[:node_count] for effect, valid in zip(sims["effects"], sims["valid_sims"]) if valid]
        valid_count = len(valid_effects)
        if valid_count:
            positive, negative = (c.astype(int) for c in _sign_counts(valid_effects))
        else:
            negative = np.zeros(node_count, dtype=int)
            positive = np.zeros(node_count, dtype=int)
        has_effect = tmat[:, p_idx] != 0
        no_effect = np.where(has_effect, 0, valid_count).astype(int)
        negative = np.where(has_effect, negative, 0).astype(int)
        positive = np.where(has_effect, positive, 0).astype(int)

        for i, node in enumerate(response_nodes):
            row = {
                "model": model_idx,
                "effect_on": node,
                "negative": int(negative[i]),
                "no_effect": int(no_effect[i]),
                "positive": int(positive[i]),
                "valid_sims": int(valid_count),
                "stable_sims": int(sims["n_stable"]),
                "attempts": int(sims["attempts"]),
            }
            rows.append(row)

    cols = ["model", "effect_on", "negative", "no_effect", "positive", "valid_sims", "stable_sims", "attempts"]
    return pd.DataFrame(rows, columns=cols)

table_of_direct_effects

table_of_direct_effects(G, form='net')

Create a table of direct effects with state/input columns and state/output rows.

Source code in qmm/extensions/effects.py

def table_of_direct_effects(
    G: nx.DiGraph,
    form: Literal["net", "absolute", "positive", "negative"] = "net",
) -> pd.DataFrame:
    """Create a table of direct effects with state/input columns and state/output rows."""
    return _tabulate_effects(G, direct_effects(G, form=form))

table_of_effects

table_of_effects(G, generator=net_effects, decimals=None, **kwargs)

Source code in qmm/extensions/effects.py

def table_of_effects(
    G: nx.DiGraph,
    generator: Union[
        Callable[..., Union[sp.Matrix, np.ndarray, pd.DataFrame]],
        Literal[
            "net_effects",
            "absolute_effects",
            "weighted_effects",
            "sign_determinacy_effects",
            "simulation_effects",
        ],
    ] = net_effects,
    decimals: Optional[int] = None,
    **kwargs: Any,
) -> pd.DataFrame:
    if isinstance(generator, str):
        generator = _EFFECT_GENERATORS.get(generator)
    if not callable(generator):
        raise ValueError(f"Generator must be callable, got: {type(generator)}")
    effects = generator(G, **kwargs)
    if decimals is not None and isinstance(effects, sp.MatrixBase):
        effects = effects.evalf(decimals)
    return _tabulate_effects(G, effects)

qmm.extensions.indicators

Identify informative indicators from press perturbations.

mutual_information cached

mutual_information(
    models,
    perturb,
    n_sim=10000,
    seed=42,
    include_null=False,
    average_uncertain=False,
)

Calculate mutual information of variables for alternative models.

Parameters:

Name	Type	Description	Default
models	Union[DiGraph, List[DiGraph]]	One or more NetworkX DiGraphs representing alternative models	required
perturb	str	Node and sign to perturb (can be comma-separated for multiple perturbations)	required
n_sim	int	Number of simulations	10000
seed	int	Random seed	42
include_null	bool	If True, include a null model with equal probability (1/3) of positive, negative, or NaN response across simulations	False

Returns:

Type	Description
DataFrame	pd.DataFrame: Mutual information for indicator selection

References

• Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.
• Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

Examples:

from qmm import mutual_information, load_digraph
G1 = load_digraph("snowshoe_rp")
G2 = G1.copy()
G2.remove_edge('C', 'P')
mutual_information((G1, G2), perturb='R:+', n_sim=1000)
#   Node  Mutual Information
# 0    R            0.513827
# 1    P            0.456516
# 2    C            0.153711

Source code in qmm/extensions/indicators.py

@cache
def mutual_information(models: Union[nx.DiGraph, List[nx.DiGraph]], perturb: str, n_sim: int = 10000, seed: int = 42, include_null: bool = False, average_uncertain: bool = False) -> pd.DataFrame:
    """Calculate mutual information of variables for alternative models.

    Args:
        models: One or more NetworkX DiGraphs representing alternative models
        perturb: Node and sign to perturb (can be comma-separated for multiple perturbations)
        n_sim: Number of simulations
        seed: Random seed
        include_null: If True, include a null model with equal probability (1/3)
            of positive, negative, or NaN response across simulations

    Returns:
        pd.DataFrame: Mutual information for indicator selection

    References:
        - Hosack, G.R., Hayes, K.R., Dambacher, J.M. (2008). Assessing Model Structure Uncertainty Through an Analysis of System Feedback and Bayesian Networks. Ecological Applications 18, 1070–1082.
        - Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

    Examples:
        ```python
        from qmm import mutual_information, load_digraph
        G1 = load_digraph("snowshoe_rp")
        G2 = G1.copy()
        G2.remove_edge('C', 'P')
        mutual_information((G1, G2), perturb='R:+', n_sim=1000)
        #   Node  Mutual Information
        # 0    R            0.513827
        # 1    P            0.456516
        # 2    C            0.153711
        ```
    """
    models = [models] if not isinstance(models, (list, tuple)) else list(models)
    models = [nx.DiGraph(G) if not isinstance(G, nx.DiGraph) else G for G in models]
    nodes = sorted(set(node for G in models for node in get_nodes(G, "state") + get_nodes(G, "output")))
    all_effects = []
    for G in models:
        G_modified, perturb_tuple = _parse_perturbations(G, perturb)
        sims = get_simulations(G_modified, n_sim=n_sim, seed=seed, perturb=perturb_tuple, average_uncertain=average_uncertain)
        response_nodes = get_nodes(G, "state") + get_nodes(G, "output")
        node_map = {node: i for i, node in enumerate(response_nodes)}
        sim_effects = []
        for effect in sims["effects"]:
            sim_effects.append([effect[node_map[n]] if n in node_map and node_map[n] < len(effect) else np.nan for n in nodes])
        all_effects.append(np.array(sim_effects))
    if include_null:
        rng = np.random.RandomState(seed)
        choices = rng.choice([1.0, -1.0, 0.0], size=(n_sim, len(nodes)))
        all_effects.append(choices)
    n_models = len(all_effects)
    mi_vals = []
    for i, node in enumerate(nodes):
        node_effects = [effects[:, i] for effects in all_effects]
        if any(np.all(np.isnan(effects)) for effects in node_effects):
            mi_vals.append(0)
            continue
        node_effects = np.concatenate(node_effects)
        labels = np.concatenate([np.full(effects.shape[0], i) for i, effects in enumerate(all_effects)])
        valid = ~np.isnan(node_effects)
        node_effects, labels = node_effects[valid], labels[valid]
        effect_signs = np.sign(node_effects) + 1
        joint, _, _ = np.histogram2d(labels, effect_signs, bins=(n_models, 3))
        joint_p = joint / joint.sum()
        l_p, e_p = joint_p.sum(axis=1), joint_p.sum(axis=0)
        mi = sum(joint_p[i, j] * np.log(joint_p[i, j] / (l_p[i] * e_p[j]))
                 for i in range(n_models)
                 for j in range(3)
                 if joint_p[i, j] > 0)
        mi_vals.append(max(0, mi))
    return pd.DataFrame({"Node": nodes, "Mutual Information": mi_vals}).sort_values("Mutual Information", ascending=False).reset_index(drop=True)

qmm.extensions.validation

Validate qualitative predictions of system response to press perturbations from observations.

marginal_likelihood cached

marginal_likelihood(
    G,
    perturb,
    observe,
    n_sim=10000,
    dist="uniform",
    seed=42,
    average_uncertain=False,
)

Calculate proportion of simulations matching qualitative observations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
perturb	str	Perturbation string (node:sign, comma-separated allowed)	required
observe	str	Observation string (node:sign, comma-separated allowed)	required
n_sim	int	Number of simulations	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution for sampling ('uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom')	'uniform'
seed	int	Random seed	42
average_uncertain	bool	Passed through to get_simulations (structure averaging over uncertain links)	False

Returns:

Name	Type	Description
float	float	Marginal likelihood

References

• Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.
• Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

Examples:

from qmm import marginal_likelihood, load_digraph
marginal_likelihood(load_digraph("snowshoe_io"), perturb='Inp1:+', observe='Out1:+', n_sim=1000)
# 0.526

Source code in qmm/extensions/validation.py

@cache
def marginal_likelihood(
    G: nx.DiGraph,
    perturb: str,
    observe: str,
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
    average_uncertain: bool = False,
) -> float:
    """Calculate proportion of simulations matching qualitative observations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        perturb: Perturbation string (node:sign, comma-separated allowed)
        observe: Observation string (node:sign, comma-separated allowed)
        n_sim: Number of simulations
        dist: Distribution for sampling ('uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom')
        seed: Random seed
        average_uncertain: Passed through to get_simulations (structure averaging over uncertain links)

    Returns:
        float: Marginal likelihood

    References:
        - Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.
        - Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

    Examples:
        ```python
        from qmm import marginal_likelihood, load_digraph
        marginal_likelihood(load_digraph("snowshoe_io"), perturb='Inp1:+', observe='Out1:+', n_sim=1000)
        # 0.526
        ```
    """
    graph, pert = _parse_perturbations(G, perturb)
    sims = get_simulations(graph, n_sim=n_sim, dist=dist, seed=seed,
                          perturb=pert,
                          observe=_parse_observations(observe) if observe else None,
                          average_uncertain=average_uncertain)
    return sum(sims["valid_sims"]) / n_sim

model_validation cached

model_validation(
    G, perturb, observe, n_sim=10000, dist="uniform", seed=42, combinations=True
)

Compare marginal likelihoods from alternative model structures.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
perturb	str	Perturbation string (node:sign, comma-separated allowed)	required
observe	str	Observation string (node:sign, comma-separated allowed)	required
n_sim	int	Number of simulations	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution for sampling	'uniform'
seed	int	Random seed	42
combinations	bool	If True, evaluate every combination of dashed edges. If False, only compare the full model vs. all dashed edges removed.	True

Returns:

Type	Description
DataFrame	pd.DataFrame: Marginal likelihood comparison for requested dashed-edge configurations.

References

• Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.
• Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

Examples:

from qmm import model_validation, load_digraph
import networkx as nx
G = nx.DiGraph(load_digraph("snowshoe_io"))
G.add_edge('R', 'P', sign=1, dashes=True)
model_validation(G, perturb='Inp1:+', observe='Out1:+', n_sim=1000, combinations=False)
#   Marginal likelihood R $\rightarrow$ P
# 0               0.815                 ✓
# 1               0.526

Source code in qmm/extensions/validation.py

@cache
def model_validation(
    G: nx.DiGraph,
    perturb: str,
    observe: str,
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
    combinations: bool = True,
) -> pd.DataFrame:
    """Compare marginal likelihoods from alternative model structures.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        perturb: Perturbation string (node:sign, comma-separated allowed)
        observe: Observation string (node:sign, comma-separated allowed)
        n_sim: Number of simulations
        dist: Distribution for sampling
        seed: Random seed
        combinations: If True, evaluate every combination of dashed edges. If False, only compare the full model vs. all dashed edges removed.

    Returns:
        pd.DataFrame: Marginal likelihood comparison for requested dashed-edge configurations.

    References:
        - Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.
        - Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

    Examples:
        ```python
        from qmm import model_validation, load_digraph
        import networkx as nx
        G = nx.DiGraph(load_digraph("snowshoe_io"))
        G.add_edge('R', 'P', sign=1, dashes=True)
        model_validation(G, perturb='Inp1:+', observe='Out1:+', n_sim=1000, combinations=False)
        #   Marginal likelihood R $\\rightarrow$ P
        # 0               0.815                 ✓
        # 1               0.526
        ```
    """
    dashed_edges = [(u, v) for u, v, d in G.edges(data=True) if d.get("dashes", False)]
    if not dashed_edges:
        mask_values = [0]
    elif combinations:
        mask_values = range(2 ** len(dashed_edges))
    else:
        mask_values = [(1 << len(dashed_edges)) - 1, 0]

    variants, edge_presence = [], []
    for mask in mask_values:
        G_variant = G.copy()
        presence = [bool(mask & (1 << j)) for j in range(len(dashed_edges))]
        for j, (u, v) in enumerate(dashed_edges):
            if not presence[j]:
                G_variant.remove_edge(u, v)
        variants.append(G_variant)
        edge_presence.append(presence)

    likelihoods = [marginal_likelihood(g, perturb, observe, n_sim, dist, seed) for g in variants]
    edge_cols = [_arrows(G, [u, v]) for u, v in dashed_edges]
    rows = [
        {"Marginal likelihood": likelihoods[i], **{edge_cols[j]: "\u2713" if edge_presence[i][j] else "" for j in range(len(dashed_edges))}}
        for i in range(len(variants))
    ]
    df = pd.DataFrame(rows, columns=["Marginal likelihood"] + edge_cols)
    df = df.sort_values("Marginal likelihood", ascending=False, kind="mergesort").reset_index(drop=True)
    df["Marginal likelihood"] = df["Marginal likelihood"].apply(lambda x: f"{x:.3f}")
    return df

posterior_predictions cached

posterior_predictions(
    G,
    perturb,
    observe="",
    n_sim=10000,
    dist="uniform",
    seed=42,
    positive_only=False,
    presample=None,
    average_uncertain=False,
)

Calculate model predictions conditioned on observations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
perturb	str	Perturbation string (node:sign, comma-separated allowed)	required
observe	str	Observation string (node:sign, comma-separated allowed)	''
n_sim	int	Number of simulations	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution for sampling	'uniform'
seed	int	Random seed	42
positive_only	bool	Return just the proportion of positive responses instead of sign-dominant proportions	False
presample	Optional[Callable[[Tuple[Symbol, ...]], Dict[Symbol, Any]]]	Optional callable passed through to get_simulations	None
average_uncertain	bool	Passed through to get_simulations (structure averaging over uncertain links)	False

Returns:

Type	Description
Matrix	sp.Matrix: Predictions conditioned on observations

References

• Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.
• Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

Examples:

from qmm import posterior_predictions, load_digraph
posterior_predictions(load_digraph("snowshoe_io"), perturb='Inp1:+', observe='Out1:+', n_sim=1000)
# Matrix([
# [              1.0],
# [-0.52851711026616],
# [-0.52851711026616],
# [              1.0],
# [-0.52851711026616]])

Source code in qmm/extensions/validation.py

@cache
def posterior_predictions(
    G: nx.DiGraph,
    perturb: str,
    observe: str = "",
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
    positive_only: bool = False,
    presample: Optional[Callable[[Tuple[sp.Symbol, ...]], Dict[sp.Symbol, Any]]] = None,
    average_uncertain: bool = False,
) -> sp.Matrix:
    """Calculate model predictions conditioned on observations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        perturb: Perturbation string (node:sign, comma-separated allowed)
        observe: Observation string (node:sign, comma-separated allowed)
        n_sim: Number of simulations
        dist: Distribution for sampling
        seed: Random seed
        positive_only: Return just the proportion of positive responses instead of sign-dominant proportions
        presample: Optional callable passed through to get_simulations
        average_uncertain: Passed through to get_simulations (structure averaging over uncertain links)

    Returns:
        sp.Matrix: Predictions conditioned on observations

    References:
        - Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.
        - Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

    Examples:
        ```python
        from qmm import posterior_predictions, load_digraph
        posterior_predictions(load_digraph("snowshoe_io"), perturb='Inp1:+', observe='Out1:+', n_sim=1000)
        # Matrix([
        # [              1.0],
        # [-0.52851711026616],
        # [-0.52851711026616],
        # [              1.0],
        # [-0.52851711026616]])
        ```
    """
    graph, pert = _parse_perturbations(G, perturb)
    observations = _parse_observations(observe) if observe else None
    sims = get_simulations(graph, n_sim=n_sim, dist=dist, seed=seed,
                          perturb=pert, observe=observations, presample=presample,
                          average_uncertain=average_uncertain)

    state, outputs = get_nodes(G, "state"), get_nodes(G, "output")
    n_total = len(state) + len(outputs)
    valid_indices = [i for i, v in enumerate(sims["valid_sims"]) if v]
    valid_count = len(valid_indices)

    if valid_count == 0:
        return sp.Matrix([np.nan] * n_total)

    effects = np.array([sims["effects"][i][:n_total] for i in valid_indices])
    positive = np.sum(effects > 0, axis=0)
    negative = np.sum(effects < 0, axis=0)

    smat = positive / valid_count if positive_only else np.where(
        negative > positive, -negative / valid_count, positive / valid_count
    )

    p_idx = sims["all_nodes"].index(pert[0])
    tmat = sims["tmat"]
    smat = [sp.nan if not tmat[i, p_idx] else smat[i] for i in range(n_total)]

    return sp.Matrix(smat)

diagnose_observations

diagnose_observations(
    G, observe, perturb_nodes=None, n_sim=10000, dist="uniform", seed=42
)

Identify possible perturbations from marginal likelihoods.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
observe	str	Observation string (node:sign, comma-separated allowed)	required
perturb_nodes	Union[str, List[str]]	Node subset to test - comma-separated string, 'state', 'input', or list of nodes	None
n_sim	int	Number of simulations	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution for sampling	'uniform'
seed	int	Random seed	42

Returns:

Type	Description
DataFrame	pd.DataFrame: Ranked perturbations matching observations

Examples:

from qmm import diagnose_observations, load_digraph
diagnose_observations(load_digraph("snowshoe_io"), observe='Out1:+', perturb_nodes='input', n_sim=1000)
#   Input Sign  Marginal likelihood
# 0  Inp2    -                1.000
# 1  Inp1    +                0.526
# 2  Inp1    -                0.474
# 3  Inp2    +                0.000

Source code in qmm/extensions/validation.py

def diagnose_observations(
    G: nx.DiGraph,
    observe: str,
    perturb_nodes: Union[str, List[str]] = None,
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
) -> pd.DataFrame:
    """Identify possible perturbations from marginal likelihoods.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        observe: Observation string (node:sign, comma-separated allowed)
        perturb_nodes: Node subset to test - comma-separated string, 'state', 'input', or list of nodes
        n_sim: Number of simulations
        dist: Distribution for sampling
        seed: Random seed

    Returns:
        pd.DataFrame: Ranked perturbations matching observations

    Examples:
        ```python
        from qmm import diagnose_observations, load_digraph
        diagnose_observations(load_digraph("snowshoe_io"), observe='Out1:+', perturb_nodes='input', n_sim=1000)
        #   Input Sign  Marginal likelihood
        # 0  Inp2    -                1.000
        # 1  Inp1    +                0.526
        # 2  Inp1    -                0.474
        # 3  Inp2    +                0.000
        ```
    """
    if perturb_nodes is None:
        perturb_nodes = get_nodes(G, "state") + get_nodes(G, "input")
    elif isinstance(perturb_nodes, str):
        if perturb_nodes == "state":
            perturb_nodes = get_nodes(G, "state")
        elif perturb_nodes == "input":
            perturb_nodes = get_nodes(G, "input")
        else:
            perturb_nodes = [node.strip() for node in perturb_nodes.split(",")]

    results = []
    for node in perturb_nodes:
        for sign in ["+", "-"]:
            try:
                likelihood = marginal_likelihood(G, f"{node}:{sign}", observe, n_sim, dist, seed)
            except RuntimeError:
                likelihood = np.nan
            results.append({"Input": node, "Sign": sign, "Marginal likelihood": likelihood})

    if not results:
        return pd.DataFrame(columns=["Input", "Sign", "Marginal likelihood"])
    return pd.DataFrame(results).sort_values(
        "Marginal likelihood", ascending=False, na_position="last"
    ).reset_index(drop=True)

bayes_factors

bayes_factors(
    G_list, perturb, observe, n_sim=10000, dist="uniform", seed=42, names=None
)

Calculate Bayes factors from the ratio of marginal likelihoods of alternative models.

Parameters:

Name	Type	Description	Default
G_list	Union[List[DiGraph], Tuple[DiGraph, ...]]	List or tuple of NetworkX DiGraphs representing alternative models	required
perturb	str	Perturbation string (node:sign, comma-separated allowed)	required
observe	str	Observation string (node:sign, comma-separated allowed)	required
n_sim	int	Number of simulations	10000
dist	Literal['uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom']	Distribution for sampling ('uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom')	'uniform'
seed	int	Random seed	42
names	Optional[List[str]]	Optional list of model names	None

Returns:

Type	Description
DataFrame	pd.DataFrame: DataFrame containing Bayes factors

References

• Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.
• Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

Examples:

from qmm import bayes_factors, load_digraph
G1 = load_digraph("snowshoe_io")
G2 = G1.copy()
G2.remove_edge('C', 'P')
bayes_factors([G1, G2], perturb='Inp1:+', observe='Out1:+', n_sim=1000)
#   Model comparison  Likelihood 1  Likelihood 2  Bayes factor
# 0  Model A/Model B         0.526         0.474      1.109705

Source code in qmm/extensions/validation.py

def bayes_factors(
    G_list: Union[List[nx.DiGraph], Tuple[nx.DiGraph, ...]],
    perturb: str,
    observe: str,
    n_sim: int = 10000,
    dist: Literal["uniform", "weak", "moderate", "strong", "uniform_two_oom"] = "uniform",
    seed: int = 42,
    names: Optional[List[str]] = None,
) -> pd.DataFrame:
    """Calculate Bayes factors from the ratio of marginal likelihoods of alternative models.

    Args:
        G_list: List or tuple of NetworkX DiGraphs representing alternative models
        perturb: Perturbation string (node:sign, comma-separated allowed)
        observe: Observation string (node:sign, comma-separated allowed)
        n_sim: Number of simulations
        dist: Distribution for sampling ('uniform', 'weak', 'moderate', 'strong', 'uniform_two_oom')
        seed: Random seed
        names: Optional list of model names

    Returns:
        pd.DataFrame: DataFrame containing Bayes factors

    References:
        - Raymond, B., McInnes, J., Dambacher, J.M., Way, S., Bergstrom, D.M. (2011). Qualitative modelling of invasive species eradication on subantarctic Macquarie Island. Journal of Applied Ecology 48, 181–191.
        - Melbourne-Thomas, J., Wotherspoon, S., Raymond, B., Constable, A. (2012). Comprehensive evaluation of model uncertainty in qualitative network analyses. Ecological Monographs 82, 505–519.

    Examples:
        ```python
        from qmm import bayes_factors, load_digraph
        G1 = load_digraph("snowshoe_io")
        G2 = G1.copy()
        G2.remove_edge('C', 'P')
        bayes_factors([G1, G2], perturb='Inp1:+', observe='Out1:+', n_sim=1000)
        #   Model comparison  Likelihood 1  Likelihood 2  Bayes factor
        # 0  Model A/Model B         0.526         0.474      1.109705
        ```
    """
    graphs = list(G_list) if isinstance(G_list, tuple) else G_list
    likelihoods = [marginal_likelihood(g, perturb, observe, n_sim, dist, seed) for g in graphs]
    model_names = names if names and len(names) == len(graphs) else [f"Model {chr(65+i)}" for i in range(len(graphs))]

    comparisons = [(i, j) for i in range(len(graphs)) for j in range(i + 1, len(graphs))]
    factors = {
        f"{model_names[i]}/{model_names[j]}": (
            float("inf") if likelihoods[j] == 0 and likelihoods[i] > 0 else
            0 if likelihoods[j] == 0 else likelihoods[i] / likelihoods[j]
        ) for i, j in comparisons
    }

    return pd.DataFrame({
        "Model comparison": list(factors.keys()),
        "Likelihood 1": [likelihoods[i] for i, _ in comparisons],
        "Likelihood 2": [likelihoods[j] for _, j in comparisons],
        "Bayes factor": list(factors.values()),
    })