econiac.routing.attribution

econiac.routing.attribution

Thermal Shapley values: differentiable attribution with latent bottleneck detection.

φ_i(β) = Gibbs-weighted average marginal contribution of player i. Λ_i(β) = |∂φ_i/∂β| — latent bottleneck index. As β→∞: Laplace concentration on the bottleneck permutation. Pacioli constraint: Σ φ_i(β) = 0 for inter-sectoral SFC attribution.

Reference: Buckley (2026) Thermal Attribution, doi:10.5281/zenodo.20236870

ShapleyResult dataclass

Container for thermal Shapley analysis outputs.

Source code in src/econiac/routing/attribution.py

@dataclass
class ShapleyResult:
    """Container for thermal Shapley analysis outputs."""
    phi:          jax.Array    # shape (n_players,)  — Shapley values
    beta:         float
    n_players:    int

    @property
    def bottleneck_player(self) -> int:
        """Player with the largest absolute Shapley value."""
        return int(jnp.argmax(jnp.abs(self.phi)))

    @property
    def total_value(self) -> float:
        """Sum of Shapley values = v(grand coalition)."""
        return float(self.phi.sum())

    def __repr__(self) -> str:
        phi_str = ", ".join(f"{float(v):.4f}" for v in self.phi)
        return (
            f"ShapleyResult(β={self.beta:.3f}, "
            f"φ=[{phi_str}], "
            f"bottleneck=player_{self.bottleneck_player})"
        )

bottleneck_player property

Player with the largest absolute Shapley value.

total_value property

Sum of Shapley values = v(grand coalition).

thermal_shapley(value_fn, n_players, beta=1.0)

Thermal Shapley values φ_i(β) for each player.

The standard Shapley value is the uniform average of marginal contributions across all n! permutations. The thermal generalisation weights permutations by Gibbs weights proportional to exp(β · v(σ)), where v(σ) is the total value generated by permutation σ.

At β=0: uniform weights → classical Shapley values. At β→∞: all weight on the max-value permutation (bottleneck).

Parameters:

Name	Type	Description	Default
value_fn	Callable[[frozenset], float]	characteristic function v: frozenset → float. v(∅) = 0 is assumed.	required
n_players	int	number of players (0, 1, …, n_players-1)	required
beta	float	inverse temperature	1.0

Returns:

Type	Description
Array	shape (n_players,) — thermal Shapley values, summing to v(grand coalition).

Source code in src/econiac/routing/attribution.py

def thermal_shapley(
    value_fn:  Callable[[frozenset], float],
    n_players: int,
    beta:      float = 1.0,
) -> jax.Array:
    """
    Thermal Shapley values φ_i(β) for each player.

    The standard Shapley value is the uniform average of marginal contributions
    across all n! permutations. The thermal generalisation weights permutations
    by Gibbs weights proportional to exp(β · v(σ)), where v(σ) is the
    total value generated by permutation σ.

    At β=0: uniform weights → classical Shapley values.
    At β→∞: all weight on the max-value permutation (bottleneck).

    Args:
        value_fn:  characteristic function v: frozenset → float.
                   v(∅) = 0 is assumed.
        n_players: number of players (0, 1, …, n_players-1)
        beta:      inverse temperature

    Returns:
        shape (n_players,) — thermal Shapley values, summing to v(grand coalition).
    """
    players = list(range(n_players))
    all_perms = list(permutations(players))

    perm_values = jnp.array(
        [_permutation_utility(value_fn, perm) for perm in all_perms],
        dtype=jnp.float32,
    )
    w = gibbs_weights(perm_values, beta=float(beta))

    phi = jnp.zeros(n_players)
    for k, perm in enumerate(all_perms):
        mc = jnp.array(
            [_marginal_contribution(value_fn, player, perm) for player in players],
            dtype=jnp.float32,
        )
        phi = phi + w[k] * mc

    return phi

bottleneck_index(value_fn, n_players, beta_range)

Latent bottleneck index Λ_i(β) = |∂φ_i/∂β| for each player over beta_range.

A sharp peak in Λ_i near some β* indicates player i is the bottleneck that dominates value creation as the system cools (β increases).

Uses finite differences over the provided beta_range.

Parameters:

Name	Type	Description	Default
value_fn	Callable[[frozenset], float]	characteristic function	required
n_players	int	number of players	required
beta_range	Array	shape (m,) — beta values at which to evaluate Λ	required

Returns:

Type	Description
Array	shape (m-1, n_players) — |∂φ_i/∂β| at interior beta points.

Source code in src/econiac/routing/attribution.py

def bottleneck_index(
    value_fn:   Callable[[frozenset], float],
    n_players:  int,
    beta_range: jax.Array,
) -> jax.Array:
    """
    Latent bottleneck index Λ_i(β) = |∂φ_i/∂β| for each player over beta_range.

    A sharp peak in Λ_i near some β* indicates player i is the bottleneck
    that dominates value creation as the system cools (β increases).

    Uses finite differences over the provided beta_range.

    Args:
        value_fn:   characteristic function
        n_players:  number of players
        beta_range: shape (m,) — beta values at which to evaluate Λ

    Returns:
        shape (m-1, n_players) — |∂φ_i/∂β| at interior beta points.
    """
    betas = np.array(beta_range)
    phis  = np.array([
        np.array(thermal_shapley(value_fn, n_players, float(b)))
        for b in betas
    ])
    dbeta  = np.diff(betas)[:, None]        # (m-1, 1)
    dphi   = np.diff(phis, axis=0)          # (m-1, n_players)
    Lambda = np.abs(dphi / (dbeta + 1e-10))
    return jnp.array(Lambda)

tropical_limit(value_fn, n_players)

Bottleneck player index in the tropical (β→∞) limit.

As β→∞, Gibbs weights concentrate on the permutation σ with maximum total value. The bottleneck player is the one with the largest marginal contribution in σ (Laplace method: the contribution that most determines whether σ* beats all other permutations).

Returns:

Type	Description
int	Index of the bottleneck player (0-indexed).

Source code in src/econiac/routing/attribution.py

def tropical_limit(
    value_fn:  Callable[[frozenset], float],
    n_players: int,
) -> int:
    """
    Bottleneck player index in the tropical (β→∞) limit.

    As β→∞, Gibbs weights concentrate on the permutation σ* with maximum
    total value. The bottleneck player is the one with the largest marginal
    contribution in σ* (Laplace method: the contribution that most determines
    whether σ* beats all other permutations).

    Returns:
        Index of the bottleneck player (0-indexed).
    """
    players  = list(range(n_players))
    all_perm = list(permutations(players))

    perm_vals = [_permutation_utility(value_fn, perm) for perm in all_perm]
    best_perm = all_perm[int(np.argmax(perm_vals))]

    mc_best = [
        _marginal_contribution(value_fn, player, best_perm)
        for player in players
    ]
    return int(np.argmax(mc_best))

nonassociative_shapley(value_fn, n_players, beta=1.0)

Non-associative Shapley values: average thermal_shapley over all Catalan trees.

The C_{n-1} rooted binary trees on n players encode different coalition-formation orders. For each tree, only coalitions that appear as subtrees are admissible. The non-associative Shapley value is the average over all tree topologies.

At n=2: C_1=1 tree → coincides with classical Shapley. At n=3: C_2=2 trees → average of left- and right-associative brackets.

Parameters:

Name	Type	Description	Default
value_fn	Callable[[frozenset], float]	characteristic function (called only on admissible coalitions)	required
n_players	int	number of players	required
beta	float	inverse temperature for Gibbs weighting within each tree	1.0

Returns:

Type	Description
Array	shape (n_players,) — non-associative Shapley values.

Source code in src/econiac/routing/attribution.py

def nonassociative_shapley(
    value_fn:  Callable[[frozenset], float],
    n_players: int,
    beta:      float = 1.0,
) -> jax.Array:
    """
    Non-associative Shapley values: average thermal_shapley over all Catalan trees.

    The C_{n-1} rooted binary trees on n players encode different coalition-formation
    orders. For each tree, only coalitions that appear as subtrees are admissible.
    The non-associative Shapley value is the average over all tree topologies.

    At n=2: C_1=1 tree → coincides with classical Shapley.
    At n=3: C_2=2 trees → average of left- and right-associative brackets.

    Args:
        value_fn:  characteristic function (called only on admissible coalitions)
        n_players: number of players
        beta:      inverse temperature for Gibbs weighting within each tree

    Returns:
        shape (n_players,) — non-associative Shapley values.
    """
    trees = _all_binary_trees(list(range(n_players)))

    def tree_shapley(tree) -> jax.Array:
        admissible_leaves = sorted(_tree_leaves(tree))
        n_adm = len(admissible_leaves)
        all_perms = list(permutations(admissible_leaves))
        perm_values = jnp.array(
            [_permutation_utility(value_fn, perm) for perm in all_perms],
            dtype=jnp.float32,
        )
        w = gibbs_weights(perm_values, beta=float(beta))
        phi = jnp.zeros(n_players)
        for k, perm in enumerate(all_perms):
            mc = jnp.zeros(n_players)
            for player in admissible_leaves:
                mc_val = _marginal_contribution(value_fn, player, perm)
                mc = mc.at[player].set(mc_val)
            phi = phi + w[k] * mc
        return phi

    tree_phis = jnp.stack([tree_shapley(t) for t in trees])
    return tree_phis.mean(axis=0)

pacioli_attribution(value_fn, n_sectors, beta=1.0)

Pacioli-constrained thermal Shapley values.

In an SFC model the grand coalition has v(N) = 0 (double-entry: every asset is someone else's liability). Standard Shapley values then sum to 0 by efficiency: Σ φ_i = v(N) = 0.

This function computes thermal_shapley and enforces the constraint explicitly by subtracting the mean — a small correction for numerical drift at finite β.

The Pacioli constraint Σ φ_i = 0 corresponds to ∂²=0 (Godley's conservation law) in the SFC balance sheet.

Parameters:

Name	Type	Description	Default
value_fn	Callable[[frozenset], float]	characteristic function with v(grand_coalition) ≈ 0	required
n_sectors	int	number of sectors (players)	required
beta	float	inverse temperature	1.0

Returns:

Type	Description
Array	shape (n_sectors,) — sector attributions with zero sum (Pacioli-normalised).

Source code in src/econiac/routing/attribution.py

def pacioli_attribution(
    value_fn:  Callable[[frozenset], float],
    n_sectors: int,
    beta:      float = 1.0,
) -> jax.Array:
    """
    Pacioli-constrained thermal Shapley values.

    In an SFC model the grand coalition has v(N) = 0 (double-entry:
    every asset is someone else's liability). Standard Shapley values
    then sum to 0 by efficiency: Σ φ_i = v(N) = 0.

    This function computes thermal_shapley and enforces the constraint
    explicitly by subtracting the mean — a small correction for
    numerical drift at finite β.

    The Pacioli constraint Σ φ_i = 0 corresponds to ∂²=0 (Godley's
    conservation law) in the SFC balance sheet.

    Args:
        value_fn:  characteristic function with v(grand_coalition) ≈ 0
        n_sectors: number of sectors (players)
        beta:      inverse temperature

    Returns:
        shape (n_sectors,) — sector attributions with zero sum (Pacioli-normalised).
    """
    phi = thermal_shapley(value_fn, n_sectors, beta)
    return phi - phi.mean()

full_shapley_analysis(value_fn, n_players, beta=1.0, beta_range=None)

Run full thermal Shapley analysis.

Returns:

Type	Description
ShapleyResult	(ShapleyResult, bottleneck_index_array or None)
Optional[Array]	bottleneck_index_array has shape (m-1, n_players) if beta_range given.

Source code in src/econiac/routing/attribution.py

def full_shapley_analysis(
    value_fn:   Callable[[frozenset], float],
    n_players:  int,
    beta:       float = 1.0,
    beta_range: Optional[jax.Array] = None,
) -> tuple[ShapleyResult, Optional[jax.Array]]:
    """
    Run full thermal Shapley analysis.

    Returns:
        (ShapleyResult, bottleneck_index_array or None)
        bottleneck_index_array has shape (m-1, n_players) if beta_range given.
    """
    phi = thermal_shapley(value_fn, n_players, beta)
    result = ShapleyResult(phi=phi, beta=beta, n_players=n_players)

    Lambda = None
    if beta_range is not None:
        Lambda = bottleneck_index(value_fn, n_players, beta_range)

    return result, Lambda