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)

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)
    """
    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"}
        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
        G.add_edge(source, target, **att)
    nx.set_node_attributes(G, "state", "category")
    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	str	Type of matrix elements ('symbolic', 'signed', or 'binary')	'symbolic'
matrix_type	str	Type of matrix to create ('A', 'B', 'C', or 'D')	'A'

Returns:

Type	Description
Matrix	sp.Matrix: Interaction matrix

Source code in qmm/core/structure.py

def create_matrix(G: nx.DiGraph, form: str = "symbolic", matrix_type: str = "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
    """

    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:  # form == 'binary'
            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))
    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	str	Type of equations to create ('state' or 'output')	'state'

Returns:

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

Source code in qmm/core/structure.py

def create_equations(G: nx.DiGraph, form: str = "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
    """
    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 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

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

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
    """
    A = sp.matrix2numpy(create_matrix(G, form="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))),
        np.linalg.det(A) != 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	str	Type of feedback ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

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

Source code in qmm/core/stability.py

@cache
def system_feedback(G: nx.DiGraph, level: Optional[int] = None, form: str = "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', or 'binary')

    Returns:
        sp.Matrix: Feedback cycle products at specified levels
    """
    A = create_matrix(G, form=form)
    if level == 0:
        return sp.Matrix([-1])
    n = A.shape[0]
    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

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
    """
    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	str	Method for computing feedback ('combinations' or 'polynomial')	'combinations'

Returns:

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

Source code in qmm/core/stability.py

@cache
def absolute_feedback(G: nx.DiGraph, level: Optional[int] = None, method: str = "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
    """
    A = create_matrix(G, form="signed")
    if level == 0:
        return sp.Matrix([1])
    n = A.shape[0]
    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="glynn") for c in combinations(range(n), k))
                fb.append(int(fb_k))
        else:
            fb_k = sum(perm(A[np.ix_(c, c)], method="glynn") for c in combinations(range(n), level))
            fb = [int(fb_k)]
    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)]
    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

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
    """
    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

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
    """
    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	str	Type of computation ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

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

Source code in qmm/core/stability.py

@cache
def hurwitz_determinants(G: nx.DiGraph, level: Optional[int] = None, form: str = "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', 'signed', or 'binary')

    Returns:
        sp.Matrix: Hurwitz determinants at specified levels
    """
    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 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

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
    """
    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

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
    """
    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(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(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

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
    """
    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

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
    """
    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

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
    """
    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 cached

simulation_stability(G, n_sim=10000)

Analyse stability using randomly sampled interaction strengths from a uniform 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

Returns:

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

Source code in qmm/core/stability.py

@cache
def simulation_stability(G: nx.DiGraph, n_sim: int = 10000) -> pd.DataFrame:
    """Analyse stability using randomly sampled interaction strengths from a uniform distribution.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        n_sim: Number of simulations to perform (default 10000)

    Returns:
        pd.DataFrame: Proportion of stable matrices and proportion that fail Hurwitz criteria
    """
    A = create_matrix(G, "signed")
    A = sp.matrix2numpy(A).astype(int)
    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 _ in range(n_sim):
        M = np.random.rand(*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 i in range(1, n):
            for j in range(1, n):
                index = 2 * j - i
                if 0 <= index < n:
                    H[i - 1, j - 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	str	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

Source code in qmm/core/press.py

@cache
def adjoint_matrix(G: nx.DiGraph, form: str = "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
    """
    A = create_matrix(G, form=form)
    A = sp.Matrix(-A)
    nodes = get_nodes(G, "state")
    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

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
    """
    A = create_matrix(G, form="binary")
    A_np = np.array(sp.matrix2numpy(A), dtype=int)
    nodes = get_nodes(G, "state")
    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.astype(float)))
        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.astype(float)))
    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

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
    """
    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	str	Method for computing determinacy ('average', '95_bound', 'simulation')	'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

Source code in qmm/core/press.py

@cache
def sign_determinacy_matrix(G: nx.DiGraph, method: str = "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', 'simulation')
        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
    """
    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	str	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 of proportions	False
match_adjoint	bool	If true, counts the proportion of simulations matching the sign of the adjoint matrix predictions.	False

Returns:

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

Raises:

Type	Description
ValueError	If invalid parameter combinations are used.

Source code in qmm/core/press.py

@cache
def numerical_simulations(G: nx.DiGraph, n_sim: int = 10000, dist: str = "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 of proportions
        match_adjoint: If true, counts the proportion of simulations matching the sign of the adjoint matrix predictions.

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

    Raises:
        ValueError: If invalid parameter combinations are used.
    """
    if positive_only and not as_nan:
        raise ValueError("Invalid parameter combination: positive_only=True requires as_nan=False")
    if as_abs and not as_nan:
        raise ValueError("Invalid parameter combination: as_abs=True requires as_nan=True")
    if match_adjoint:
        if as_abs:
            raise ValueError("Invalid parameter combination: match_adjoint=True requires as_abs=False")
        if not as_nan:
            raise ValueError("Invalid parameter combination: match_adjoint=True requires as_nan=True")
        if positive_only:
            raise ValueError("Invalid parameter combination: match_adjoint=True requires positive_only=False")

    np.random.seed(seed)
    A = create_matrix(G, form="symbolic", matrix_type="A")
    state_nodes = get_nodes(G, "state")
    n = len(state_nodes)
    symbols = list(A.free_symbols)
    A_sp = sp.lambdify(symbols, A)
    dist_funcs = {
        "uniform": lambda size: np.random.uniform(0, 1, size),
        "weak": lambda size: np.random.beta(1, 3, size),
        "moderate": lambda size: np.random.beta(2, 2, size),
        "strong": lambda size: np.random.beta(3, 1, size),
        "normal_weak": lambda size: truncnorm.rvs(a=0, b=3, loc=0, scale=1/3, size=size),
        "normal_moderate": lambda size: truncnorm.rvs(a=-3, b=3, loc=0.5, scale=1/6, size=size),
        "normal_strong": lambda size: truncnorm.rvs(a=-3, b=0, loc=1, scale=1/3, size=size)
    }
    if match_adjoint:
        adj_mat = adjoint_matrix(G, form="signed")
        adj_sign_np = np.array(adj_mat.applyfunc(sp.sign).tolist(), dtype=float)
        matches = np.zeros((n, n), dtype=int)
    else:
        positive = np.zeros((n, n), dtype=int)
        negative = np.zeros((n, n), dtype=int)
    total_simulations = 0
    while total_simulations < n_sim:
        values = dist_funcs[dist](len(symbols))
        sim_A = A_sp(*values)
        if np.all(np.real(np.linalg.eigvals(sim_A)) < 0):
            try:
                inv_A = np.linalg.inv(-sim_A)
                if match_adjoint:
                    matches += (np.sign(inv_A) == adj_sign_np)
                else:
                    positive += inv_A > 0
                    negative += inv_A < 0
                total_simulations += 1
            except np.linalg.LinAlgError:
                continue
    if total_simulations == 0:
        smat = np.full((n, n), np.nan)
    elif match_adjoint:
        smat = matches / total_simulations
    elif positive_only:
        smat = positive / total_simulations
    else:
        smat = np.where(negative > positive, -negative / total_simulations, positive / total_simulations)
    smat = sp.Matrix(smat.astype(float).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 and not match_adjoint:
            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:
        smat = sp.Matrix([[0 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.

table_of_predictions

table_of_predictions(M, t1=0.8, t2=1.0, 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 likely predictions	0.8
t2	float	Upper threshold for determined predictions	1.0
index	Optional[List[str]]	Row labels	None
columns	Optional[List[str]]	Column labels	None

Returns:

Type	Description
DataFrame	pd.DataFrame: Qualitative predictions

Source code in qmm/core/prediction.py

def table_of_predictions(M: Union[sp.Matrix, np.ndarray], t1: float = 0.8, t2: float = 1.0,
                        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 likely predictions
        t2 (float): Upper threshold for determined predictions
        index (Optional[List[str]]): Row labels
        columns (Optional[List[str]]): Column labels

    Returns:
        pd.DataFrame: Qualitative predictions
    """
    if isinstance(M, sp.Matrix):
        M = sp.matrix2numpy(M, dtype=float)
    conditions = [
        (lambda x: x == sp.nan, "0"),
        (lambda x: x >= t2, "+"),
        (lambda x: x >= t1 and x < t2, "(+)"),
        (lambda x: x > -t1 and x < t1, "?"),
        (lambda x: x > -t2 and x <= -t1, "(\u2212)"),
        (lambda x: x <= -t2, "\u2212"),
    ]
    predictions = np.empty(M.shape, dtype=object)
    for i in range(M.shape[0]):
        for j in range(M.shape[1]):
            value = M[i, j]
            predictions[i, j] = next((val for cond, val in conditions if cond(value)), "0")
    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

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
    """
    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.

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.
    """
    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]))
    nx.set_node_attributes(G, "state", "category")
    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.

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.
    """
    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	str	Type of nodes to extract ('state' 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.

Source code in qmm/core/helper.py

def get_nodes(G: nx.DiGraph, node_type: str = "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' or 'all').
        labels: If True, return node labels instead of node ids.

    Returns:
        List of node identifiers or dictionaries containing node data.
    """
    if not isinstance(G, nx.DiGraph):
        raise TypeError("Input must be a networkx.DiGraph.")

    if node_type == "all":
        return list(G.nodes()) if not labels else list(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.

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.
    """
    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.

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.
    """
    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.

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.
    """
    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	str	Method to use for probability calculation ('average' or '95_bound').	'average'

Returns:

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

Source code in qmm/core/helper.py

def sign_determinacy(wmat: sp.Matrix, tmat: sp.Matrix, method: str = "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.
    """

    MAX_PROB = sp.Float('0.999999')

    def compute_prob(w, t, method):
        if w == sp.Integer(0):
            return sp.Rational(1, 2)
        elif w == sp.Integer(1):
            return sp.Integer(1)
        elif w == sp.Integer(-1):
            return sp.Integer(-1)
        elif 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:  # exp(700) is near the float64 limit
            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

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

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
    """
    A = create_matrix(G, "signed")
    n = A.shape[0]
    fcp = system_feedback(G)[1:]
    if level is None:
        level = n
    S = sp.zeros(n, n)
    nodes = get_nodes(G, "state")
    for i in range(n):
        for j in range(n):
            if A[i, j] != 0:
                sG = nx.DiGraph(G)
                sG[nodes[j]][nodes[i]]["sign"] = 0
                scp = system_feedback(sG)[1:]
                if level <= len(fcp) and level <= len(scp):
                    N = fcp[level - 1] - scp[level - 1]
                    S[i, j] = N
    return S

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

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
    """
    A = create_matrix(G, "signed")
    n = A.shape[0]
    fcp = net_feedback(G)[1:]
    if level is None:
        level = n
    S = sp.zeros(n, n)
    nodes = get_nodes(G, "state")
    for i in range(n):
        for j in range(n):
            if A[i, j] != 0:
                sG = nx.DiGraph(G)
                sG[nodes[j]][nodes[i]]["sign"] = 0
                scp = net_feedback(sG)[1:]
                if level <= len(fcp) and level <= len(scp):
                    N = fcp[level - 1] - scp[level - 1]
                    S[i, j] = N
    return S

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

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
    """
    A = create_matrix(G, "signed")
    n = A.shape[0]
    fcp = absolute_feedback(G)[1:]
    if level is None:
        level = n
    S = sp.zeros(n, n)
    nodes = get_nodes(G, "state")
    for i in range(n):
        for j in range(n):
            if A[i, j] != 0:
                sG = nx.DiGraph(G)
                sG[nodes[j]][nodes[i]]["sign"] = 0
                scp = absolute_feedback(sG)[1:]
                if level <= len(fcp) and level <= len(scp):
                    N = fcp[level - 1] - scp[level - 1]
                    S[i, j] = N
    return S

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

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
    """
    A = create_matrix(G, "signed")
    n = A.shape[0]
    if level is None:
        level = n
    net = sp.Matrix(net_structural_sensitivity(G, level))
    absolute = sp.Matrix(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	str	Type of computation ('symbolic', 'signed')	'symbolic'

Returns:

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

Source code in qmm/extensions/life.py

@cache
def birth_matrix(G: nx.DiGraph, form: str = "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
    """
    A_sgn = create_matrix(G, form="signed")
    A_sym = create_matrix(G, form="symbolic")
    nodes = get_nodes(G, "state")
    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:  # form == 'signed'
            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	str	Type of computation ('symbolic', 'signed')	'symbolic'

Returns:

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

Source code in qmm/extensions/life.py

@cache
def death_matrix(G: nx.DiGraph, form: str = "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
    """
    A_sgn = create_matrix(G, form="signed")
    A_sym = create_matrix(G, form="symbolic")
    nodes = get_nodes(G, "state")
    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:  # form == 'signed'
            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	str	Change in birth or death rate ('birth' or 'death')	'birth'
form	str	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

Source code in qmm/extensions/life.py

@cache
def life_expectancy_change(G: nx.DiGraph, form: str = "symbolic", type: str = "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
    """
    amat = adjoint_matrix(G, form=form)
    if type == "birth":
        matrix = death_matrix(G, form=form)
    else:  # type == 'death'
        matrix = birth_matrix(G, form=form)
    result = sp.expand(sp.Integer(-1) * matrix * amat)
    if perturb is not None:
        nodes = get_nodes(G, "state")
        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	str	Change in birth or death rate ('birth' or 'death')	'birth'

Returns:

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

Source code in qmm/extensions/life.py

@cache
def net_life_expectancy_change(G: nx.DiGraph, type: str = "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
    """
    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	str	Change in birth or death rate ('birth' or 'death')	'birth'

Returns:

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

Source code in qmm/extensions/life.py

@cache
def absolute_life_expectancy_change(G: nx.DiGraph, type: str = "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
    """
    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	str	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

Source code in qmm/extensions/life.py

@cache
def weighted_predictions_life_expectancy(G: nx.DiGraph, type: str = "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
    """
    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

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
    """
    A = create_matrix(G, form="symbolic")
    nodes = get_nodes(G, "state")
    node_id = {n: i for i, n in enumerate(nodes)}
    cycle_list = nx.simple_cycles(G)
    cycle_nodes = sorted([c for c in cycle_list], key=lambda x: len(x))
    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

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
    """
    cycle_nodes = sorted([path for path in nx.simple_cycles(G)], key=lambda x: (len(x), x))
    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	required
form	str	Type of path products ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

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

Source code in qmm/extensions/paths.py

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

    Args:
        G: NetworkX DiGraph representing signed digraph model
        source: Source node
        target: Target node
        form: Type of path products ('symbolic', 'signed', or 'binary')

    Returns:
        sp.Matrix: Products of interactions along each path
    """
    nodes = get_nodes(G, "all")
    A = create_matrix(G, form=form)
    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 sp.Matrix([sp.Integer(0)])
    path_nodes = list(nx.all_simple_paths(G, source, target))
    paths = [sp.prod(A[nodes.index(p[i + 1]), nodes.index(p[i])] for i in range(len(p) - 1)) for p in path_nodes]
    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

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
    """
    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))
    if not paths:
        return None
    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 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	required
form	str	Type of feedback ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

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

Source code in qmm/extensions/paths.py

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

    Args:
        G: NetworkX DiGraph representing signed digraph model
        source: Source node
        target: Target node
        form: Type of feedback ('symbolic', 'signed', or 'binary')

    Returns:
        sp.Matrix: Feedback cycles in complementary subsystem
    """
    nodes = get_nodes(G, "state")
    n = len(nodes)
    if source not in nodes or target not in nodes or source == target:
        raise ValueError("Invalid source or target node")
    path_nodes = list(nx.all_simple_paths(G, source, target))
    feedback = []
    for path in path_nodes:
        path_nodes_set = set(path)
        subsystem_nodes = [node for node in nodes if node not in path_nodes_set]
        if not subsystem_nodes:
            if form == "binary":
                feedback.append(sp.Integer(1))
            else:
                feedback.append(sp.Integer(-1))
            continue
        subsystem = G.subgraph(subsystem_nodes).copy()
        level = n - len(path)
        if form == "symbolic":
            feedback.append(system_feedback(subsystem, level=level)[0])
        elif form == "signed":
            feedback.append(net_feedback(subsystem, level=level)[0])
        elif form == "binary":
            feedback.append(absolute_feedback(subsystem, level=level)[0])
        else:
            raise ValueError("Invalid form. Choose 'symbolic', 'signed', or 'binary'.")
    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	required
form	str	Type of computation ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

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

Source code in qmm/extensions/paths.py

@cache
def system_paths(G: nx.DiGraph, source: str, target: str, form: str = "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
        form: Type of computation ('symbolic', 'signed', or 'binary')

    Returns:
        sp.Matrix: Total effects including paths and feedback
    """
    if source == target:
        raise ValueError("Source and target must be different nodes")
    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

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
    """
    nodes = get_nodes(G, "state")
    A_sgn = create_matrix(G, form="signed")
    if source not in nodes or target not in 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 = [node for node in nodes if node not in path]
        if not subsystem_nodes:
            feedback = sp.Integer(-1)
        else:
            subsystem = G.subgraph(subsystem_nodes).copy()
            feedback = weighted_feedback(subsystem, level=len(nodes) - len(path))
            if feedback[0] == sp.nan:
                feedback = sp.Integer(0)
        sign = sp.prod(A_sgn[nodes.index(path[i + 1]), nodes.index(path[i])] for i in range(len(path) - 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

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
    """
    nodes = get_nodes(G, "state")
    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 pd.DataFrame()
    paths = list(nx.all_simple_paths(G, source, target))
    complementary_nodes = [[node for node in nodes if node 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)
    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 complementary_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)],
        }
    )

    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)

Define model components as state variables, inputs and outputs.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
remove_disconnected	bool	Remove disconnected components	True

Returns:

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

Source code in qmm/extensions/effects.py

def define_input_output(G: nx.DiGraph, remove_disconnected: bool = True) -> nx.DiGraph:
    """Define model components as state variables, inputs and outputs.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        remove_disconnected: Remove disconnected components

    Returns:
        nx.DiGraph: Model with input, state and output classification
    """
    if not isinstance(G, nx.DiGraph):
        raise TypeError("Input must be a networkx.DiGraph.")
    G_def = G.copy()
    if remove_disconnected:
        G_undirected = G_def.to_undirected()
        connected = list(nx.connected_components(G_undirected))
        if len(connected) > 1:
            largest = max(connected, key=len)
            nodes_to_remove = [node for system in connected if system != largest for node in system]
            G_def.remove_nodes_from(nodes_to_remove)
    nx.set_node_attributes(G_def, "state", "category")
    while True:
        reclassified = False
        for node in list(G_def.nodes()):
            if G_def.nodes[node]["category"] == "state":
                if all(G_def.nodes[pred]["category"] == "input" for pred in G_def.predecessors(node)):
                    G_def.nodes[node]["category"] = "input"
                    reclassified = True
                elif all(G_def.nodes[succ]["category"] == "output" for succ in G_def.successors(node)):
                    G_def.nodes[node]["category"] = "output"
                    reclassified = True
        if not reclassified:
            break
    return G_def

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	str	Type of computation ('symbolic', 'signed', or 'binary')	'symbolic'

Returns:

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

Source code in qmm/extensions/effects.py

@cache
def cumulative_effects(G: nx.DiGraph, form: str = "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
    """
    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")
    if form == "symbolic":
        effects = adjoint_matrix(G, form="symbolic")
    elif form in "signed":
        effects = adjoint_matrix(G, form="signed")
    elif form in "binary":
        effects = absolute_feedback_matrix(G)
    else:
        raise ValueError("Invalid form. Choose 'symbolic', 'signed', 'binary'.")
    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)

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

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
    """
    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

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
    """
    net = cumulative_effects(G, form="signed")
    absolute = cumulative_effects(G, form="binary")
    return get_weight(net, absolute)

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	str	Method for computing determinacy ('average', '95_bound')	'average'

Returns:

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

Source code in qmm/extensions/effects.py

@cache
def sign_determinacy_effects(G: nx.DiGraph, method: str = "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
    """
    weighted = weighted_effects(G)
    absolute = cumulative_effects(G, form="binary")
    return sign_determinacy(weighted, absolute, method=method)

get_simulations cached

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

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	str	Distribution for sampling	'uniform'
seed	int	Random seed	42
perturb	Optional[Tuple[str, int]]	Optional tuple of (node, sign) to perturb	None
observe	Optional[Tuple[Tuple[str, int], ...]]	Optional tuple of observations as (node, sign) tuples	None

Returns:

Type	Description
Dict[str, Any]	Dict containing effects, valid_sims, all_nodes, and tmat

Source code in qmm/extensions/effects.py

@cache
def get_simulations(G: nx.DiGraph, n_sim: int = 10000, dist: str = "uniform", seed: int = 42, perturb: Optional[Tuple[str, int]] = None, observe: Optional[Tuple[Tuple[str, int], ...]] = None) -> 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 tuple of (node, sign) to perturb
        observe: Optional tuple of observations as (node, sign) tuples

    Returns:
        Dict containing effects, valid_sims, all_nodes, and tmat
    """
    np.random.seed(seed)
    A, B, C, D = [create_matrix(G, form="symbolic", matrix_type=m) for m in "ABCD"]
    state_nodes, input_nodes, output_nodes = [get_nodes(G, t) for t in ["state", "input", "output"]]
    all_nodes = state_nodes + input_nodes + output_nodes
    node_idx = {node: i for i, node in enumerate(all_nodes)}
    tmat = sp.matrix2numpy(absolute_effects(G)).astype(int)
    dist_funcs = {
        "uniform": lambda size: np.random.uniform(0, 1, size),
        "weak": lambda size: np.random.beta(1, 3, size),
        "moderate": lambda size: np.random.beta(2, 2, size),
        "strong": lambda size: np.random.beta(3, 1, size),
        "normal_weak": lambda size: truncnorm.rvs(a=0, b=3, loc=0, scale=1/3, size=size),
        "normal_moderate": lambda size: truncnorm.rvs(a=-3, b=3, loc=0.5, scale=1/6, size=size),
        "normal_strong": lambda size: truncnorm.rvs(a=-3, b=0, loc=1, scale=1/3, size=size)
    }
    pert_idx, perturb_sign = (node_idx[perturb[0]], perturb[1]) if perturb else (None, 1)
    n_state, n_input, n_output = len(state_nodes), len(input_nodes), len(output_nodes)
    n_rows, n_cols = n_state + n_output, n_state + n_input
    symbols = list(set(A.free_symbols) | set(B.free_symbols) | set(C.free_symbols) | set(D.free_symbols))
    A_sp = sp.lambdify(symbols, A)
    B_sp = sp.lambdify(symbols, B) if n_input > 0 else None
    C_sp = sp.lambdify(symbols, C) if n_output > 0 else None
    D_sp = sp.lambdify(symbols, D) if D.shape != (0, 0) else None
    effects, valid_sims = [], []
    for _ in range(n_sim * 100):
        values = dist_funcs[dist](len(symbols))
        sim_A = A_sp(*values)
        if np.all(np.real(np.linalg.eigvals(sim_A)) < 0):
            try:
                inv_A = np.linalg.inv(-sim_A)
                B_np = B_sp(*values) if B_sp else np.array([]).reshape(n_state, 0)
                D_np = D_sp(*values) if D_sp else np.array([]).reshape(n_output, n_input)
                C_np = C_sp(*values) if C_sp else np.array([]).reshape(0, n_state)
                effect_matrix = np.zeros((n_rows, n_cols))
                if n_state > 0:
                    effect_matrix[:n_state, :n_state] = inv_A
                    if n_input > 0:
                        effect_matrix[:n_state, n_state:] = inv_A @ B_np
                if n_output > 0:
                    if n_state > 0:
                        effect_matrix[n_state:, :n_state] = C_np @ inv_A
                        if n_input > 0:
                            effect_matrix[n_state:, n_state:] = C_np @ inv_A @ B_np + D_np
                    elif n_input > 0:
                        effect_matrix[n_state:, n_state:] = D_np
                effect = effect_matrix[:, pert_idx] * perturb_sign if pert_idx is not None else effect_matrix
                effects.append(effect)
                if observe:
                    valid = all(
                        (
                            node in state_nodes
                            and (
                                (tmat[node_idx[node], pert_idx] == 0 and obs == 0)
                                or (obs != 0 and np.sign(effect[node_idx[node]]) == obs)
                            )
                        )
                        or (
                            node in output_nodes
                            and (
                                (tmat[len(state_nodes) + output_nodes.index(node), pert_idx] == 0 and obs == 0)
                                or (obs != 0 and np.sign(effect[len(state_nodes) + output_nodes.index(node)]) == obs)
                            )
                        )
                        for node, obs in observe
                    )
                    valid_sims.append(valid)
                else:
                    valid_sims.append(True)
                if len(effects) == n_sim:
                    break
            except np.linalg.LinAlgError:
                continue
    else:
        raise RuntimeError(f"Maximum iterations reached. Stable proportion: {len(effects) / (n_sim * 100):.4f}")
    return {"effects": effects, "valid_sims": valid_sims, "all_nodes": all_nodes, "tmat": tmat}

simulation_effects

simulation_effects(
    G, n_sim=10000, dist="uniform", seed=42, positive_only=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	str	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

Returns:

Type	Description
Matrix	SymPy Matrix containing simulation results

Source code in qmm/extensions/effects.py

def simulation_effects(G: nx.DiGraph, n_sim: int = 10000, dist: str = "uniform", seed: int = 42, positive_only: 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

    Returns:
        SymPy Matrix containing simulation results
    """
    sims = get_simulations(G, n_sim, dist, seed)
    tmat = sims["tmat"]
    state_nodes = get_nodes(G, "state")
    input_nodes = get_nodes(G, "input")
    output_nodes = get_nodes(G, "output")
    n_state, n_input, n_output = len(state_nodes), len(input_nodes), len(output_nodes)
    n_rows, n_cols = n_state + n_output, n_state + n_input
    positive = np.zeros((n_rows, n_cols), dtype=int)
    negative = np.zeros((n_rows, n_cols), dtype=int)
    for effect in sims["effects"]:
        positive += effect > 0
        negative += effect < 0
    if positive_only:
        smat = positive / n_sim
    else:
        smat = np.where(negative > positive, -negative / n_sim, positive / n_sim)
    smat = sp.Matrix(smat)
    smat = sp.Matrix([[sp.nan if not tmat[i, j] else smat[i, j] for j in range(smat.cols)] for i in range(smat.rows)])
    return smat

qmm.extensions.indicators

Identify informative indicators from press perturbations.

mutual_information cached

mutual_information(models, perturb, n_sim=10000, seed=42)

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	required
n_sim	int	Number of simulations	10000
seed	int	Random seed	42

Returns:

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

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) -> 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
        n_sim: Number of simulations
        seed: Random seed

    Returns:
        pd.DataFrame: Mutual information for indicator selection
    """
    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:
        sims = get_simulations(G, n_sim=n_sim, seed=seed, perturb=_NodeSign.from_str(perturb).to_tuple())
        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))
    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]
        if len(node_effects) == 0:
            mi_vals.append(0)
            continue
        joint, _, _ = np.histogram2d(labels, node_effects > 0, bins=(len(models), 2))
        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.log2(joint_p[i, j] / (l_p[i] * e_p[j] + 1e-10))
                 for i in range(len(models))
                 for j in range(2)
                 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, distribution="uniform", seed=42
)

Calculate proportion of simulations matching qualitative observations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
perturb	str	Node and sign to perturb	required
observe	str	String of observations	required
n_sim	int	Number of simulations	10000
distribution	str	Distribution for sampling	'uniform'
seed	int	Random seed	42

Returns:

Name	Type	Description
float	float	Marginal likelihood

Source code in qmm/extensions/validation.py

@cache
def marginal_likelihood(G: nx.DiGraph, perturb: str, observe: str, n_sim: int = 10000, distribution: str = "uniform", seed: int = 42) -> float:
    """Calculate proportion of simulations matching qualitative observations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        perturb: Node and sign to perturb
        observe: String of observations
        n_sim: Number of simulations
        distribution: Distribution for sampling 
        seed: Random seed

    Returns:
        float: Marginal likelihood
    """
    sims = get_simulations(G, n_sim=n_sim, dist=distribution, seed=seed,
                          perturb=_NodeSign.from_str(perturb).to_tuple(),
                          observe=_parse_observations(observe) if observe else None)
    return sum(sims["valid_sims"]) / n_sim

model_validation cached

model_validation(
    G, perturb, observe, n_sim=10000, distribution="uniform", seed=42
)

Compare marginal likelihoods from alternative model structures.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
perturb	str	Node and sign to perturb	required
observe	str	String of observations	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Marginal likelihood comparison for model variants

Source code in qmm/extensions/validation.py

@cache
def model_validation(G: nx.DiGraph, perturb: str, observe: str, n_sim: int = 10000, distribution: str = "uniform", seed: int = 42) -> pd.DataFrame:
    """Compare marginal likelihoods from alternative model structures.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        perturb: Node and sign to perturb
        observe: String of observations

    Returns:
        pd.DataFrame: Marginal likelihood comparison for model variants
    """
    dashed_edges = [(u, v) for u, v, d in G.edges(data=True) if d.get("dashes", False)]
    variants = []
    edge_presence = []
    for i in range(2 ** len(dashed_edges)):
        G_variant = G.copy()
        presence = []
        for j, (u, v) in enumerate(dashed_edges):
            if not (i & (1 << j)):
                G_variant.remove_edge(u, v)
            presence.append(bool(i & (1 << j)))
        variants.append(G_variant)
        edge_presence.append(presence)

    likelihoods = [marginal_likelihood(G, perturb, observe, n_sim, distribution, seed) for G in variants]
    columns = ["Marginal likelihood"] + [_arrows(G, [u, v]) for u, v in dashed_edges]
    rows = []
    for i in range(len(variants)):
        row = {
            "Marginal likelihood": likelihoods[i]
        }
        for j, (u, v) in enumerate(dashed_edges):
            row[_arrows(G, [u, v])] = "\u2713" if edge_presence[i][j] else ""
        rows.append(row)
    df = pd.DataFrame(rows, columns=columns).sort_values("Marginal likelihood", ascending=False).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
)

Calculate model predictions conditioned on observations.

Parameters:

Name	Type	Description	Default
G	DiGraph	NetworkX DiGraph representing signed digraph model	required
perturb	str	Node and sign to perturb	required
observe	str	String of observations	''
n_sim	int	Number of simulations	10000
dist	str	Distribution for sampling	'uniform'
seed	int	Random seed	42

Returns:

Type	Description
Matrix	sp.Matrix: Predictions conditioned on observations

Source code in qmm/extensions/validation.py

@cache
def posterior_predictions(G: nx.DiGraph, perturb: str, observe: str = "", n_sim: int = 10000, dist: str = "uniform", seed: int = 42) -> sp.Matrix:
    """Calculate model predictions conditioned on observations.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        perturb: Node and sign to perturb
        observe: String of observations
        n_sim: Number of simulations
        dist: Distribution for sampling
        seed: Random seed

    Returns:
        sp.Matrix: Predictions conditioned on observations
    """
    sims = get_simulations(G, n_sim=n_sim, dist=dist, seed=seed,
                          perturb=_NodeSign.from_str(perturb).to_tuple(),
                          observe=_parse_observations(observe) if observe else None)
    state_nodes, output_nodes = get_nodes(G, "state"), get_nodes(G, "output")
    n, m = len(state_nodes), len(output_nodes)
    valid_count = sum(sims["valid_sims"])
    tmat = sims["tmat"]
    if valid_count == 0:
        return sp.Matrix([np.nan] * (n + m))
    effects = np.array(
        [
            e[: n + m] if len(e) >= n + m else np.pad(e, (0, n + m - len(e)))
            for e, v in zip(sims["effects"], sims["valid_sims"])
            if v
        ]
    )
    positive = np.sum(effects > 0, axis=0)
    negative = np.sum(effects < 0, axis=0)
    smat = positive / valid_count
    tmat_np = np.array(tmat.tolist(), dtype=bool)
    perturb_node = _NodeSign.from_str(perturb).node
    perturb_index = sims["all_nodes"].index(perturb_node)
    smat = [np.nan if not tmat_np[i, perturb_index] else smat[i] for i in range(n + m)]

    if observe:
        observations = _parse_observations(observe)
        for node, value in observations:
            index = (
                state_nodes.index(node)
                if node in state_nodes
                else (n + output_nodes.index(node) if node in output_nodes else None)
            )
            if index is not None:
                smat[index] = 1 if value > 0 else (0 if value < 0 else np.nan)

    smat = np.where(negative > positive, -negative / valid_count, smat)
    return sp.Matrix(smat)

diagnose_observations cached

diagnose_observations(G, observe, n_sim=10000, distribution='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	String of observations	required

Returns:

Type	Description
DataFrame	pd.DataFrame: Ranked perturbations matching observations

Source code in qmm/extensions/validation.py

@cache
def diagnose_observations(G: nx.DiGraph, observe: str, n_sim: int = 10000, distribution: str = "uniform", seed: int = 42) -> pd.DataFrame:
    """Identify possible perturbations from marginal likelihoods.

    Args:
        G: NetworkX DiGraph representing signed digraph model
        observe: String of observations

    Returns:
        pd.DataFrame: Ranked perturbations matching observations
    """
    perturb_nodes = get_nodes(G, "state") + get_nodes(G, "input")
    results = []
    for node in perturb_nodes:
        for sign in ["+", "-"]:
            try:
                likelihood = marginal_likelihood(G, f"{node}:{sign}", observe, n_sim, distribution, seed)
                results.append({"Input": node, "Sign": sign, "Marginal likelihood": likelihood})
            except Exception as e:
                print(f"Error for node {node} with sign {sign}: {str(e)}")
    return pd.DataFrame(results).sort_values("Marginal likelihood", ascending=False).reset_index(drop=True)

bayes_factors

bayes_factors(
    G_list,
    perturb,
    observe,
    n_sim=10000,
    distribution="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	Node and sign to perturb	required
observe	str	String of observations	required
n_sim	int	Number of simulations	10000
distribution	str	Distribution for sampling	'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

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, distribution: str = "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: Node and sign to perturb
        observe: String of observations
        n_sim: Number of simulations
        distribution: Distribution for sampling
        seed: Random seed
        names: Optional list of model names

    Returns:
        pd.DataFrame: DataFrame containing Bayes factors
    """
    # Convert tuple to list if needed
    G_list = list(G_list) if isinstance(G_list, tuple) else G_list
    likelihoods = [marginal_likelihood(G, perturb, observe, n_sim, distribution, seed) for G in G_list]
    model_names = names if names and len(names) == len(G_list) else [f"Model {chr(65+i)}" for i in range(len(G_list))]
    bayes_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 in range(len(G_list)) for j in range(i + 1, len(G_list))}

    return pd.DataFrame({
        "Model comparison": list(bayes_factors.keys()),
        "Likelihood 1": [likelihoods[i] for i in range(len(G_list)) for j in range(i + 1, len(G_list))],
        "Likelihood 2": [likelihoods[j] for i in range(len(G_list)) for j in range(i + 1, len(G_list))],
        "Bayes factor": list(bayes_factors.values()),
    })