Symbolic Lowering & Calculus (mkpp.lowering)

mkpp.lowering

SparsityOptimizer

Analyzes Jacobian sparsity for fill-in prediction, reordering, and block detection.

Source code in src/mkpp/lowering.py

class SparsityOptimizer:
    """Analyzes Jacobian sparsity for fill-in prediction, reordering, and block detection."""

    def __init__(self, jacobian_structure: "set[tuple[int, int]]", n: int):
        self.structure = jacobian_structure
        self.n = n

    def predict_fill_in(self) -> "set[tuple[int, int]]":
        """
        Graph-reachability fill-in prediction (symbolic Gaussian elimination).
        For Doolittle LU, position (i,j) fills in if there exists k < min(i,j)
        such that both (i,k) and (k,j) are structurally non-zero (transitively).
        """
        active = set(self.structure)
        active |= {(i, i) for i in range(self.n)}  # diagonal always present
        fill = set()

        for k in range(self.n):
            rows_with_k = [i for i in range(k + 1, self.n) if (i, k) in active]
            cols_with_k = [j for j in range(k + 1, self.n) if (k, j) in active]

            for i in rows_with_k:
                for j in cols_with_k:
                    if (i, j) not in active:
                        fill.add((i, j))
                        active.add((i, j))

        return fill

    def compute_rcm_ordering(self) -> "list[int]":
        """
        Reverse Cuthill-McKee on the symmetrized structure graph.
        Returns permutation vector p where new_index = p[old_index].
        """
        G = nx.Graph()
        G.add_nodes_from(range(self.n))
        for i, j in self.structure:
            if i != j:
                G.add_edge(i, j)

        # Handle disconnected graphs (no edges means no bandwidth to reduce)
        if G.number_of_edges() == 0:
            return list(range(self.n))

        # NetworkX provides Cuthill-McKee ordering directly
        # cuthill_mckee_ordering returns nodes in CM order; reverse for RCM
        cm_order = list(nx.utils.rcm.cuthill_mckee_ordering(G))
        # Reverse for RCM (Reverse Cuthill-McKee)
        rcm_order = list(reversed(cm_order))

        def _bw(order):
            inv = [0] * self.n
            for new_idx, old_idx in enumerate(order):
                inv[old_idx] = new_idx
            return max(abs(inv[i] - inv[j]) for i, j in self.structure) if self.structure else 0

        orig_order = list(range(self.n))
        if _bw(rcm_order) > _bw(orig_order):
            return orig_order

        return rcm_order

    def detect_blocks(self) -> "list[list[int]]":
        """
        Run Tarjan SCC on the directed Jacobian structure graph.
        Each SCC with no cross-block edges becomes an independent block.
        """
        G = nx.DiGraph()
        G.add_nodes_from(range(self.n))
        for i, j in self.structure:
            if i != j:
                G.add_edge(i, j)

        sccs = list(nx.strongly_connected_components(G))
        # Sort blocks deterministically by minimum index
        blocks = sorted([sorted(list(scc)) for scc in sccs], key=lambda b: b[0])
        return blocks

    def _check_block_independence(self, blocks: "list[list[int]]") -> bool:
        """
        Check if blocks are independent (no cross-block non-zeros after fill-in).
        Returns True if the system is truly block-diagonal.
        """
        if len(blocks) <= 1:
            return False  # Single block is not "block-diagonal"

        # Build a mapping from species index to block index
        species_to_block: dict[int, int] = {}
        for block_idx, block in enumerate(blocks):
            for species_idx in block:
                species_to_block[species_idx] = block_idx

        # Include fill-in positions in the check
        fill = self.predict_fill_in()
        all_positions = self.structure | fill | {(i, i) for i in range(self.n)}

        # Check for cross-block non-zeros
        for i, j in all_positions:
            if i == j:
                continue
            if species_to_block.get(i) != species_to_block.get(j):
                return False  # Cross-block coupling exists

        return True

    def analyze(self) -> "SparsityAnalysis":
        """Full sparsity analysis pipeline."""
        from .model import SparsityAnalysis

        fill = self.predict_fill_in()
        perm = self.compute_rcm_ordering()
        inv_perm = [0] * self.n
        for new_idx, old_idx in enumerate(perm):
            inv_perm[old_idx] = new_idx

        blocks = self.detect_blocks()
        is_block_diag = self._check_block_independence(blocks)

        return SparsityAnalysis(
            original_nnz=len(self.structure),
            fill_in_positions=fill,
            total_nnz_after_fill=len(self.structure) + len(fill),
            permutation=perm,
            inverse_permutation=inv_perm,
            blocks=blocks,
            is_block_diagonal=is_block_diag,
        )

analyze()

Full sparsity analysis pipeline.

Source code in src/mkpp/lowering.py

def analyze(self) -> "SparsityAnalysis":
    """Full sparsity analysis pipeline."""
    from .model import SparsityAnalysis

    fill = self.predict_fill_in()
    perm = self.compute_rcm_ordering()
    inv_perm = [0] * self.n
    for new_idx, old_idx in enumerate(perm):
        inv_perm[old_idx] = new_idx

    blocks = self.detect_blocks()
    is_block_diag = self._check_block_independence(blocks)

    return SparsityAnalysis(
        original_nnz=len(self.structure),
        fill_in_positions=fill,
        total_nnz_after_fill=len(self.structure) + len(fill),
        permutation=perm,
        inverse_permutation=inv_perm,
        blocks=blocks,
        is_block_diagonal=is_block_diag,
    )

compute_rcm_ordering()

Reverse Cuthill-McKee on the symmetrized structure graph. Returns permutation vector p where new_index = p[old_index].

Source code in src/mkpp/lowering.py

def compute_rcm_ordering(self) -> "list[int]":
    """
    Reverse Cuthill-McKee on the symmetrized structure graph.
    Returns permutation vector p where new_index = p[old_index].
    """
    G = nx.Graph()
    G.add_nodes_from(range(self.n))
    for i, j in self.structure:
        if i != j:
            G.add_edge(i, j)

    # Handle disconnected graphs (no edges means no bandwidth to reduce)
    if G.number_of_edges() == 0:
        return list(range(self.n))

    # NetworkX provides Cuthill-McKee ordering directly
    # cuthill_mckee_ordering returns nodes in CM order; reverse for RCM
    cm_order = list(nx.utils.rcm.cuthill_mckee_ordering(G))
    # Reverse for RCM (Reverse Cuthill-McKee)
    rcm_order = list(reversed(cm_order))

    def _bw(order):
        inv = [0] * self.n
        for new_idx, old_idx in enumerate(order):
            inv[old_idx] = new_idx
        return max(abs(inv[i] - inv[j]) for i, j in self.structure) if self.structure else 0

    orig_order = list(range(self.n))
    if _bw(rcm_order) > _bw(orig_order):
        return orig_order

    return rcm_order

detect_blocks()

Run Tarjan SCC on the directed Jacobian structure graph. Each SCC with no cross-block edges becomes an independent block.

Source code in src/mkpp/lowering.py

def detect_blocks(self) -> "list[list[int]]":
    """
    Run Tarjan SCC on the directed Jacobian structure graph.
    Each SCC with no cross-block edges becomes an independent block.
    """
    G = nx.DiGraph()
    G.add_nodes_from(range(self.n))
    for i, j in self.structure:
        if i != j:
            G.add_edge(i, j)

    sccs = list(nx.strongly_connected_components(G))
    # Sort blocks deterministically by minimum index
    blocks = sorted([sorted(list(scc)) for scc in sccs], key=lambda b: b[0])
    return blocks

predict_fill_in()

Graph-reachability fill-in prediction (symbolic Gaussian elimination). For Doolittle LU, position (i,j) fills in if there exists k < min(i,j) such that both (i,k) and (k,j) are structurally non-zero (transitively).

Source code in src/mkpp/lowering.py

def predict_fill_in(self) -> "set[tuple[int, int]]":
    """
    Graph-reachability fill-in prediction (symbolic Gaussian elimination).
    For Doolittle LU, position (i,j) fills in if there exists k < min(i,j)
    such that both (i,k) and (k,j) are structurally non-zero (transitively).
    """
    active = set(self.structure)
    active |= {(i, i) for i in range(self.n)}  # diagonal always present
    fill = set()

    for k in range(self.n):
        rows_with_k = [i for i in range(k + 1, self.n) if (i, k) in active]
        cols_with_k = [j for j in range(k + 1, self.n) if (k, j) in active]

        for i in rows_with_k:
            for j in cols_with_k:
                if (i, j) not in active:
                    fill.add((i, j))
                    active.add((i, j))

    return fill

annotate_lu_expressions(lu_plan)

Annotate each LU expression with the set of species indices it depends on.

For each expression in lu_expressions_ordered, determines which species indices affect that entry. An expression at (row, col) directly depends on species row and col. Additionally, any W_i_j, L_i_j, or U_i_j references in the expression add species i and j to the dependency set.

Parameters:

Name	Type	Description	Default
lu_plan	SymbolicLUPlan	The symbolic LU plan with lu_expressions_ordered populated.	required

Returns:

Type	Description
list[AnnotatedLUExpression]	List of AnnotatedLUExpression with depends_on sets populated.

Source code in src/mkpp/lowering.py

def annotate_lu_expressions(lu_plan: SymbolicLUPlan) -> "list[AnnotatedLUExpression]":
    """Annotate each LU expression with the set of species indices it depends on.

    For each expression in lu_expressions_ordered, determines which species indices
    affect that entry. An expression at (row, col) directly depends on species row
    and col. Additionally, any W_i_j, L_i_j, or U_i_j references in the expression
    add species i and j to the dependency set.

    Args:
        lu_plan: The symbolic LU plan with lu_expressions_ordered populated.

    Returns:
        List of AnnotatedLUExpression with depends_on sets populated.
    """
    import re

    from .model import AnnotatedLUExpression

    annotated = []
    for kind, row, col, expr_str in lu_plan.lu_expressions_ordered:
        # Direct dependencies: the species at row and column positions
        depends_on: set[int] = {row, col}

        # Also check for references to other W/L/U entries
        for match in re.finditer(r"[WLU]_(\d+)_(\d+)", expr_str):
            depends_on.add(int(match.group(1)))
            depends_on.add(int(match.group(2)))

        annotated.append(
            AnnotatedLUExpression(
                kind=kind,
                row=row,
                col=col,
                expr=expr_str,
                depends_on=depends_on,
            )
        )

    return annotated

apply_cse_to_plan(lu_plan, f_vector)

Apply sympy.cse() to all expressions in the LU plan + rate vector.

Returns:

Name	Type	Description
replacements	list	List of (Symbol, expression) tuples for CSE temporaries
reduced	list	List of simplified expressions with CSE symbols substituted

Source code in src/mkpp/lowering.py

def apply_cse_to_plan(lu_plan: SymbolicLUPlan, f_vector: sp.Matrix) -> tuple[list, list]:
    """Apply sympy.cse() to all expressions in the LU plan + rate vector.

    Returns:
        replacements: List of (Symbol, expression) tuples for CSE temporaries
        reduced: List of simplified expressions with CSE symbols substituted
    """
    all_exprs = []
    for i, j, expr_str in lu_plan.non_zero_jacobian:
        all_exprs.append(sp.sympify(expr_str))
    for expr in f_vector:
        all_exprs.append(expr)

    replacements, reduced = sp.cse(all_exprs, optimizations="basic")
    return replacements, reduced

build_sympy_matrices(mech)

Lowering function to compute unified Jacobian and symbolic sparse LU plan. Attaches results to mechanism metadata.

Source code in src/mkpp/lowering.py

def build_sympy_matrices(mech: MechanismDefinition) -> dict[str, Any]:
    """
    Lowering function to compute unified Jacobian and symbolic sparse LU plan.
    Attaches results to mechanism metadata.
    """
    res = prepare_unified_jacobian(mech)
    plan = res.get("symbolic_lu_plan")
    if plan is None:
        plan = compute_symbolic_lu_decomposition(res["jacobian_matrix"], res["species_map"])
        res["symbolic_lu_plan"] = plan

    if getattr(mech, "metadata", None) is None:
        mech.metadata = {}
    mech.metadata["sympy_metadata"] = res
    mech.metadata["symbolic_lu_plan"] = plan
    return res

compute_symbolic_lu_decomposition(J_matrix, species_map, permutation=None, blocks=None, is_block_diagonal=False)

Computes build-time symbolic sparse LU factorization schedule (Doolittle method) on W = inv_g_dt * I - J. Extracts flat scalar L, U matrix entry expressions and forward/backward substitution steps referencing previously computed scalar variables.

When is_block_diagonal=True and blocks is provided, compute per-block LU plans independently and combine results into a single SymbolicLUPlan with block metadata.

When permutation is provided (but not block-diagonal), apply the permutation to J before computing LU and store the permutation in the resulting plan.

Source code in src/mkpp/lowering.py

def compute_symbolic_lu_decomposition(
    J_matrix: sp.Matrix,
    species_map: list[str],
    permutation: "list[int] | None" = None,
    blocks: "list[list[int]] | None" = None,
    is_block_diagonal: bool = False,
) -> SymbolicLUPlan:
    """
    Computes build-time symbolic sparse LU factorization schedule (Doolittle method) on W = inv_g_dt * I - J.
    Extracts flat scalar L, U matrix entry expressions and forward/backward substitution steps referencing previously
    computed scalar variables.

    When is_block_diagonal=True and blocks is provided, compute per-block LU plans independently
    and combine results into a single SymbolicLUPlan with block metadata.

    When permutation is provided (but not block-diagonal), apply the permutation to J before
    computing LU and store the permutation in the resulting plan.
    """
    N = J_matrix.shape[0]
    if N == 0:
        raise ValueError("Cannot perform symbolic LU decomposition on empty matrix.")

    # Handle block-diagonal case: compute per-block LU independently
    if is_block_diagonal and blocks is not None and len(blocks) > 1:
        return _compute_block_diagonal_lu(J_matrix, species_map, blocks)

    # Handle permutation case: reorder J before LU
    if permutation is not None:
        J_perm = sp.zeros(N, N)
        for i in range(N):
            for j in range(N):
                J_perm[i, j] = J_matrix[permutation[i], permutation[j]]
        perm_species_map = [species_map[permutation[i]] for i in range(N)]
        plan = _compute_lu_core(J_perm, perm_species_map)
        plan.permutation = permutation
        plan.blocks = blocks
        return plan

    return _compute_lu_core(J_matrix, species_map)

compute_transposed_lu_plan(lu_plan)

Given an existing SymbolicLUPlan, compute the transposed substitution steps.

W * x = b => L * U * x = b

• Forward sub: L * y = b (lower-triangular)
• Backward sub: U * x = y (upper-triangular)

W^T * x = b => U^T * L^T * x = b

• Forward sub with U^T: U^T * y = b (U^T is lower-triangular)
• Backward sub with L^T: L^T * x = y (L^T is upper-triangular)

The transposed steps are stored in

• lu_plan.transpose_forward_sub_steps
• lu_plan.transpose_backward_sub_steps

Requirements: 5.1, 5.2

Source code in src/mkpp/lowering.py

def compute_transposed_lu_plan(lu_plan: SymbolicLUPlan) -> SymbolicLUPlan:
    """
    Given an existing SymbolicLUPlan, compute the transposed substitution steps.

    For the forward LU solve:  W * x = b  =>  L * U * x = b
      - Forward sub:  L * y = b  (lower-triangular)
      - Backward sub: U * x = y  (upper-triangular)

    For the transposed solve: W^T * x = b  =>  U^T * L^T * x = b
      - Forward sub with U^T:  U^T * y = b  (U^T is lower-triangular)
      - Backward sub with L^T: L^T * x = y  (L^T is upper-triangular)

    The transposed steps are stored in:
      - lu_plan.transpose_forward_sub_steps
      - lu_plan.transpose_backward_sub_steps

    Requirements: 5.1, 5.2
    """
    N = lu_plan.num_species

    # Build non-zero structure sets for L and U from the existing plan
    nz_L = set()  # (row, col) pairs where L is non-zero
    nz_U = set()  # (row, col) pairs where U is non-zero

    for i, j, _expr in lu_plan.l_expressions:
        nz_L.add((i, j))
    # L diagonal is always 1 (unit lower triangular)
    for i in range(N):
        nz_L.add((i, i))

    for i, j, _expr in lu_plan.u_expressions:
        nz_U.add((i, j))

    # Transposed forward substitution: solve U^T * y = b
    # U^T is lower-triangular. U^T[i,j] = U[j,i].
    # For i = 0, 1, ..., N-1:
    #   U^T[i,i] * y_i + Σ_{k<i} U^T[i,k] * y_k = b_i
    #   => U[i,i] * y_i + Σ_{k<i} U[k,i] * y_k = b_i
    #   => y_i = (b_i - Σ_{k<i} U[k,i] * y_k) / U[i,i]
    transpose_forward_steps = []
    for i in range(N):
        sub_terms = []
        for k in range(i):
            # U^T[i,k] = U[k,i] — check if (k, i) is in nz_U
            if (k, i) in nz_U:
                sub_terms.append(f"U_{k}_{i} * y_{k}")

        if sub_terms:
            num_str = f"b_{i} - " + " - ".join(sub_terms)
            expr_str = f"({num_str}) / U_{i}_{i}"
        else:
            expr_str = f"b_{i} / U_{i}_{i}"
        transpose_forward_steps.append((i, expr_str))

    # Transposed backward substitution: solve L^T * x = y
    # L^T is upper-triangular. L^T[i,j] = L[j,i].
    # L is unit lower-triangular so L^T[i,i] = 1.
    # For i = N-1, N-2, ..., 0:
    #   L^T[i,i] * x_i + Σ_{k>i} L^T[i,k] * x_k = y_i
    #   => x_i + Σ_{k>i} L[k,i] * x_k = y_i
    #   => x_i = y_i - Σ_{k>i} L[k,i] * x_k
    transpose_backward_steps = []
    for i in range(N - 1, -1, -1):
        sub_terms = []
        for k in range(i + 1, N):
            # L^T[i,k] = L[k,i] — check if (k, i) is in nz_L
            if (k, i) in nz_L:
                sub_terms.append(f"L_{k}_{i} * x_{k}")

        if sub_terms:
            expr_str = f"y_{i} - " + " - ".join(sub_terms)
        else:
            expr_str = f"y_{i}"
        transpose_backward_steps.append((i, expr_str))

    # Store the transposed steps in the plan
    lu_plan.transpose_forward_sub_steps = transpose_forward_steps
    lu_plan.transpose_backward_sub_steps = transpose_backward_steps

    return lu_plan

partition_reactions(mech)

Partition reactions into implicit (stiff) and explicit (non-stiff) deterministic blocks using Tarjan's Strongly Connected Components (SCC) algorithm.

Source code in src/mkpp/lowering.py

def partition_reactions(mech: MechanismDefinition) -> dict[str, list[ReactionDefinition]]:
    """
    Partition reactions into implicit (stiff) and explicit (non-stiff) deterministic blocks
    using Tarjan's Strongly Connected Components (SCC) algorithm.
    """
    blocks = {"implicit": [], "explicit": []}

    # 1. Build the directed species dependency graph
    G = nx.DiGraph()
    for r in mech.reactions:
        for reactant in r.reactants:
            for product in r.products:
                G.add_edge(reactant, product)

    # 2. Find cycles (SCCs with more than 1 node, or self-loops)
    sccs = list(nx.strongly_connected_components(G))
    stiff_species = set()
    for scc in sccs:
        if len(scc) > 1:
            stiff_species.update(scc)
        elif len(scc) == 1:
            # Check for self-loop
            node = list(scc)[0]
            if G.has_edge(node, node):
                stiff_species.add(node)

    # 3. Partition reactions based on topology
    for r in mech.reactions:
        # A reaction belongs to the stiff manifold if it connects species within the stiff network
        is_stiff_topology = any(reactant in stiff_species for reactant in r.reactants) and any(
            product in stiff_species for product in r.products
        )

        # We also respect manual overrides (r.stiff) if the user forces it
        if r.stiff or is_stiff_topology:
            blocks["implicit"].append(r)
        else:
            blocks["explicit"].append(r)

    # Sort blocks deterministically by reaction type then expression
    blocks["implicit"].sort(key=lambda x: (x.reaction_type, x.rate_expression))
    blocks["explicit"].sort(key=lambda x: (x.reaction_type, x.rate_expression))

    # T026: Inject deterministic solver partition metadata
    blocks["metadata"] = {
        "sza_sorted": True,
        "micro_blocks": {"implicit": len(blocks["implicit"]), "explicit": len(blocks["explicit"])},
        "scc_count": len([s for s in sccs if len(s) > 1]),
    }

    return blocks

prepare_adjoint_and_tlm(mech)

T015: Symbolic lowering hooks for analytical Jacobian, Adjoint, and Tangent-Linear models. For the MVP, this validates that the mechanism is differentiable.

Source code in src/mkpp/lowering.py

def prepare_adjoint_and_tlm(mech: MechanismDefinition) -> dict[str, bool]:
    """
    T015: Symbolic lowering hooks for analytical Jacobian, Adjoint, and Tangent-Linear models.
    For the MVP, this validates that the mechanism is differentiable.
    """
    # Verify no discontinuous thermodynamic operators are present
    for r in mech.reactions:
        if not r.continuous_transition and r.reaction_type.lower() in (
            "condensation",
            "phase_change",
        ):
            raise ValueError(f"Reaction {r.rate_expression} lacks continuous transition for analytical differentiation.")

    # Check that equilibrium reactions have continuous_transition = True
    if hasattr(mech, "equilibrium_reactions") and mech.equilibrium_reactions:
        for eq_def in mech.equilibrium_reactions:
            if not eq_def.continuous_transition:
                raise ValueError(
                    f"Equilibrium system '{eq_def.system}' requires continuous_transition=True " f"for analytical differentiation."
                )

    return {"adjoint_ready": True, "tlm_ready": True}

prepare_unified_jacobian_parallel(mech)

Same as prepare_unified_jacobian but with parallel Jacobian column computation.

Uses multiprocessing.Pool to compute each column df/dC_j independently, then assembles the full Jacobian in deterministic column order. Falls back to sequential computation if multiprocessing raises an error.

Source code in src/mkpp/lowering.py

def prepare_unified_jacobian_parallel(mech: MechanismDefinition) -> dict[str, Any]:
    """Same as prepare_unified_jacobian but with parallel Jacobian column computation.

    Uses multiprocessing.Pool to compute each column df/dC_j independently, then
    assembles the full Jacobian in deterministic column order.  Falls back to
    sequential computation if multiprocessing raises an error.
    """
    built = _build_f_total(mech)
    ordered_species = built["ordered_species"]
    f_implicit = built["f_implicit"]
    f_explicit = built["f_explicit"]
    f_total = built["f_total"]
    c_vector = built["c_vector"]
    built["species_symbols"]

    N = len(ordered_species)

    # For mechanisms with < 200 species, the serialization overhead of
    # multiprocessing (pickling SymPy srepr strings, spawning workers, IPC)
    # exceeds the differentiation speedup. Fall back to sequential path.
    if N < 200:
        return prepare_unified_jacobian(mech)

    # Serialize f_total expressions for multiprocessing (SymPy objects aren't picklable)
    f_total_serialized = [sp.srepr(expr) for expr in f_total]
    c_sreprs = [sp.srepr(sym) for sym in c_vector]

    try:
        with multiprocessing.Pool() as pool:
            args = [(j, f_total_serialized, c_sreprs[j]) for j in range(N)]
            results = pool.map(_compute_jacobian_column, args)

        # Assemble in deterministic column order
        results.sort(key=lambda x: x[0])
        columns = []
        for _j, col_strs in results:
            columns.append(sp.Matrix([sp.sympify(s) for s in col_strs]))
        jacobian_matrix = sp.Matrix.hstack(*[col.reshape(N, 1) for col in columns])
    except Exception as e:
        warnings.warn(f"Parallel Jacobian computation failed, falling back to sequential: {e}")
        jacobian_matrix = f_total.jacobian(c_vector)

    adjoint_matrix = jacobian_matrix.transpose()

    # Mass conservation projector
    unique_elements = sorted(list(set(elem for s in mech.species for elem in s.elements.keys())))
    if unique_elements:
        E_matrix = sp.zeros(len(unique_elements), len(ordered_species))
        for j, sp_name in enumerate(ordered_species):
            species_def = next(s for s in mech.species if s.name == sp_name)
            for i, elem in enumerate(unique_elements):
                E_matrix[i, j] = species_def.elements.get(elem, 0)
        try:
            E_E_T = E_matrix * E_matrix.transpose()
            mass_projector = E_matrix.transpose() * E_E_T.pinv()
        except Exception:
            mass_projector = sp.zeros(len(ordered_species), len(unique_elements))
    else:
        E_matrix = sp.zeros(1, len(ordered_species))
        mass_projector = sp.zeros(len(ordered_species), 1)

    result = {
        "species_map": ordered_species,
        "f_implicit": f_implicit,
        "f_explicit": f_explicit,
        "jacobian_matrix": jacobian_matrix,
        "adjoint_matrix": adjoint_matrix,
        "mass_projector": mass_projector,
        "element_map": unique_elements,
    }

    try:
        N = jacobian_matrix.shape[0]
        jacobian_structure = set()
        for i in range(N):
            for j in range(N):
                if jacobian_matrix[i, j] != 0:
                    jacobian_structure.add((i, j))

        sparsity = SparsityOptimizer(jacobian_structure, N)
        analysis = sparsity.analyze()

        lu_plan = compute_symbolic_lu_decomposition(
            jacobian_matrix,
            ordered_species,
            permutation=analysis.permutation,
            blocks=analysis.blocks,
            is_block_diagonal=analysis.is_block_diagonal,
        )
        result["symbolic_lu_plan"] = lu_plan
        result["sparsity_analysis"] = analysis
    except Exception:
        try:
            lu_plan = compute_symbolic_lu_decomposition(jacobian_matrix, ordered_species)
            result["symbolic_lu_plan"] = lu_plan
        except Exception:
            pass

    return result