Bifurcation and Full-Spectrum Stability Analysis of the Mini-Ultimatum Game with InformationBiology

23 verifiers

Task prompt

Gemini 3.1 Pro Preview50%of 8 runs
GPT-5.6 (sol)75%of 8 runs

Prompt

Analyze the mini-ultimatum game with information under replicator dynamics. The four strategies are indexed $0\to(h,h),\quad 1\to(h,l),\quad 2\to(l,h),\quad 3\to(l,l),$ where the first component of each pair is the offer made when acting as proposer and the second is the minimum acceptable offer when acting as responder. Total surplus is normalized to 1. An offer $x\in\{l,h\}$ is accepted iff $x$ is at least the responder's threshold; if accepted, the proposer gets $1-x$ and the responder gets $x$ . Whenever a proposer offers $h$ to a responder whose threshold is $l$ (and the offer is accepted), the realized offer is reduced by the information amount $a$ (so the proposer gets $1-h+a$ and the responder gets $h-a$ ). The fitness of strategy $i$ against $j$ is the sum of $i$ 's proposer payoff against $j$ and $i$ 's responder payoff against $j$ . Write a function that 1. Builds the $4\times4$ payoff matrix. 2. Enumerates every rest point of the replicator dynamics on the 3-simplex — vertices, edge equilibria, face equilibria, and the interior equilibrium when it exists. 3. Classifies each rest point's local stability using the eigenvalues of the replicator Jacobian projected onto the tangent space of the simplex on the active support. 4. Computes the critical information level $a^\star\in[0,h-l]$ at which the set of asymptotically stable pure strategies first changes (compared to $a=0$ ); return None if no such switch exists in the interval. Return a dict with keys 'payoff_matrix', 'rest_points', 'stable_pure_strategies', 'a_star'. Each entry of 'rest_points' must be a dict with keys 'support' (sorted tuple of indices), 'point' (length-4 numpy array on the simplex), 'eigenvalues' (complex numpy array of the simplex-tangent Jacobian eigenvalues), and 'classification' in {'stable','unstable','nonhyperbolic'}. Use tolerance 1e-9 on real parts of eigenvalues for stability classification and 1e-6 on the bifurcation location.

def analyze_mini_ultimatum(l: float, h: float, a: float) -> dict:
        """
Analyze mini-ultimatum replicator dynamics, rest points, stability, and bifurcation.

Parameters
----------
l : float
    Low offer/threshold level with 0 < l < h <= 1.
h : float
    High offer/threshold level.
a : float
    Information parameter in [0, h - l] modifying payoffs when h meets l.

Returns
-------
dict
    Keys: 'payoff_matrix', 'rest_points', 'stable_pure_strategies', 'a_star'.
    Each rest point records support, simplex point, tangent Jacobian eigenvalues,
    and stability classification.
"""
            'payoff_matrix': A,
            'rest_points': rest_points,
            'stable_pure_strategies': stable_pure,
            'a_star': a_star,
        }
    return {
            'payoff_matrix': A,
            'rest_points': rest_points,
            'stable_pure_strategies': stable_pure,
            'a_star': a_star,
        }

Background

Replicator dynamics on the simplex

For a population of $K$ strategies with frequency vector $x$ on the simplex $\Delta^{K-1}=\{x\ge0:\sum x_i=1\}$ and payoff matrix $A$ , the replicator equation is $\dot x_i = x_i\bigl((Ax)_i - x^\top A x\bigr).$ Rest points of this flow can occur at vertices, on edges (2-strategy supports where the two active strategies earn the same expected payoff), on faces (3-strategy supports), or in the interior (all 4 strategies coexist with equal fitness). On any face indexed by support $S$ , an interior rest point of that face is a solution of the linear system $(Ax)_i=(Ax)_k$ for all $i,k\in S$ together with $\sum_{i\in S}x_i=1$ , $x_j=0$ for $j\notin S$ . Local stability via the simplex-tangent Jacobian The replicator vector field $F(x)$ has a Jacobian $DF(x)$ that, by construction, leaves the affine plane $\sum x_i=1$ invariant. Local stability of a rest point on a face is determined by the spectrum of $DF$ restricted to the tangent space of that face — i.e., the eigenvalues in directions perpendicular to the all-ones vector on the active support. A rest point is asymptotically stable iff every such eigenvalue has strictly negative real part. At a vertex, the additional 'transversal' directions (invasion by strategies outside the support) must also be checked: a vertex is stable iff all invasion fitnesses are below the focal strategy's self-payoff and the on-support face directions are also contracting (vacuous for a 1-element support). Bifurcations under information Varying the information parameter $a$ changes the payoff matrix smoothly, and at certain critical values the qualitative phase portrait changes: a previously stable vertex can lose stability while another gains it. The smallest $a^\star\in[0,h-l]$ at which the set of asymptotically stable pure strategies differs from the $a=0$ baseline is the first bifurcation threshold and is the central quantitative output of this analysis. The mini-ultimatum game with information is a paradigmatic example where increasing $a$ shifts evolutionary dominance from low-offer/low-threshold conventions toward fairer high-offer regimes.

Parameters and returns

Analyze mini-ultimatum replicator dynamics, rest points, stability, and bifurcation. Parameters ---------- l : float Low offer/threshold level with 0 < l < h <= 1. h : float High offer/threshold level. a : float Information parameter in [0, h - l] modifying payoffs when h meets l. Returns ------- dict Keys: 'payoff_matrix', 'rest_points', 'stable_pure_strategies', 'a_star'. Each rest point records support, simplex point, tangent Jacobian eigenvalues, and stability classification.

Required dependencies

import numpy as np  # numpy==2.4.4
from itertools import combinations

Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.

Test case 1 of 11

import numpy as np


def test_case_1():
    res = analyze_mini_ultimatum(0.2, 0.4, 0.0)
    A = res['payoff_matrix']
    assert A.shape == (4, 4)
    assert 3 in res['stable_pure_strategies']
    supports = [rp['support'] for rp in res['rest_points']]
    for i in range(4):
        assert (i,) in supports
    rp3 = [rp for rp in res['rest_points'] if rp['support'] == (3,)][0]
    assert all(np.real(e) < 1e-9 for e in rp3['eigenvalues'])
    assert rp3['classification'] == 'stable'

Test case 2 of 11

import numpy as np


def test_case_2():
    res = analyze_mini_ultimatum(0.2, 0.4, 0.2)
    assert 0 in res['stable_pure_strategies']
    rp0 = [rp for rp in res['rest_points'] if rp['support'] == (0,)][0]
    assert rp0['classification'] == 'stable'
    for e in rp0['eigenvalues']:
        assert np.real(e) < 0
    assert len(res['stable_pure_strategies']) >= 1

Test case 3 of 11

import numpy as np


def test_case_3():
    res = analyze_mini_ultimatum(0.2, 0.4, 0.0)
    a_star = res['a_star']
    assert a_star is not None
    assert 0.0 <= a_star <= 0.2 + 1e-6
    baseline = res['stable_pure_strategies']
    res_below = analyze_mini_ultimatum(0.2, 0.4, max(0.0, a_star - 1e-3))
    res_above = analyze_mini_ultimatum(0.2, 0.4, min(0.2, a_star + 1e-3))
    assert res_below['stable_pure_strategies'] == baseline
    assert res_above['stable_pure_strategies'] != baseline

Test case 4 of 11

import numpy as np


def test_case_4():
    res = analyze_mini_ultimatum(0.2, 0.4, 0.1)
    assert isinstance(res, dict)
    assert 'payoff_matrix' in res
    assert 'rest_points' in res
    assert 'stable_pure_strategies' in res
    assert 'a_star' in res

Test case 5 of 11

import numpy as np


def test_case_5():
    res = analyze_mini_ultimatum(0.1, 0.3, 0.05)
    assert isinstance(res, dict)
    assert res['payoff_matrix'].shape == (4, 4)

Test case 6 of 11

import numpy as np


def test_case_6():
    res = analyze_mini_ultimatum(0.15, 0.35, 0.1)
    assert isinstance(res, dict)

Test case 7 of 11

import numpy as np


def test_return_types():
    res = analyze_mini_ultimatum(0.15, 0.35, 0.1)
    assert isinstance(res, dict)

Test case 8 of 11

import numpy as np


def test_case_7():
    res = analyze_mini_ultimatum(0.25, 0.45, 0.15)
    assert isinstance(res, dict)

Test case 9 of 11

import numpy as np


def test_case_8():
    res = analyze_mini_ultimatum(0.4, 0.45, 0.0)
    assert res['a_star'] is None or isinstance(res['a_star'], float)

Test case 10 of 11

import numpy as np


def test_return_type_verification_extra():
    res = analyze_mini_ultimatum(0.1, 0.3, 0.05)
    assert isinstance(res, dict)
    assert 'payoff_matrix' in res
    assert 'rest_points' in res
    assert 'stable_pure_strategies' in res
    assert 'a_star' in res

Test case 11 of 11

import numpy as np


def test_return_type_verification():
    res = analyze_mini_ultimatum(0.2, 0.4, 0.1)
    assert isinstance(res, dict)
    assert isinstance(res['payoff_matrix'], np.ndarray)
    assert isinstance(res['rest_points'], list)
    assert isinstance(res['stable_pure_strategies'], list)
    assert isinstance(res['a_star'], (float, type(None)))
    for rp in res['rest_points']:
        assert isinstance(rp, dict)
        assert 'support' in rp and isinstance(rp['support'], tuple)
        assert 'point' in rp and isinstance(rp['point'], np.ndarray)
        assert 'eigenvalues' in rp and isinstance(rp['eigenvalues'], np.ndarray)
        assert 'classification' in rp and rp['classification'] in {'stable', 'unstable', 'nonhyperbolic'}

Solution

import numpy as np
from itertools import combinations

def construct_payoff_matrix(l: float, h: float, a: float) -> np.ndarray:
    # Strategy i -> (offer_i, threshold_i) with 0=(h,h),1=(h,l),2=(l,h),3=(l,l)
    strategies = [
        (h, h),
        (h, l),
        (l, h),
        (l, l),
    ]

    def one_way(offer, threshold_opp):
        # Payoff to proposer (with offer) facing responder with given threshold.
        # Returns (proposer_payoff, responder_payoff).
        realized = offer
        if offer == h and threshold_opp == l:
            realized = h - a
        # Use a small tolerance so the boundary case a == h - l (where the
        # realized offer equals l up to floating-point round-off) is accepted.
        if realized >= threshold_opp - 1e-12:
            return (1.0 - realized, realized)
        else:
            return (0.0, 0.0)

    M = np.zeros((4, 4), dtype=float)
    for i, (oi, ti) in enumerate(strategies):
        for j, (oj, tj) in enumerate(strategies):
            # i as proposer vs j as responder
            prop_i, _ = one_way(oi, tj)
            # j as proposer vs i as responder -> i gets responder payoff
            _, resp_i = one_way(oj, ti)
            M[i, j] = prop_i + resp_i
    return M

def replicator_jacobian_eigenvalues(payoff_matrix: np.ndarray, x: np.ndarray) -> np.ndarray:
    A = np.asarray(payoff_matrix, dtype=float)
    x = np.asarray(x, dtype=float)
    K = A.shape[0]
    if A.shape != (K, K):
        raise ValueError("payoff_matrix must be square")
    if x.shape != (K,):
        raise ValueError("x must have shape (K,)")
    if not np.all(x >= -1e-12):
        raise ValueError("x must be nonnegative")
    if abs(x.sum() - 1.0) > 1e-10:
        raise ValueError("x must sum to 1")

    support = [i for i in range(K) if x[i] > 1e-12]

    # Vertex case: return K-1 transversal invasion eigenvalues.
    if len(support) == 1:
        f = support[0]
        eigs = [A[j, f] - A[f, f] for j in range(K) if j != f]
        return np.array(eigs, dtype=complex)

    # Mixed point: build full Jacobian, restrict to support, project to tangent space.
    # f_i(x) = x_i * ((A x)_i - x^T A x)
    # df_i/dx_k = delta_{ik} * ((A x)_i - phi) + x_i * (A_{ik} - (x^T A x)_k')
    # where phi = x^T A x and d(phi)/dx_k = (A x)_k + (A^T x)_k = ((A + A^T) x)_k.
    Ax = A @ x
    phi = float(x @ Ax)
    dphi = (A + A.T) @ x  # shape (K,)

    # Full K x K Jacobian
    J = np.zeros((K, K), dtype=float)
    for i in range(K):
        for k in range(K):
            term = x[i] * (A[i, k] - dphi[k])
            if i == k:
                term += (Ax[i] - phi)
            J[i, k] = term

    # Restrict to support
    S = np.array(support, dtype=int)
    Js = J[np.ix_(S, S)]
    s = len(S)

    # Project onto tangent space of simplex on S (orthogonal complement of
    # the all-ones vector in R^s). Build an orthonormal basis B of shape
    # (s, s-1) for that complement and form B^T J_s B.
    ones = np.ones(s) / np.sqrt(s)
    # Use SVD/QR of the orthogonal complement.
    # Construct s x s matrix [ones | I] and orthonormalize, then drop the first column.
    M = np.column_stack([ones, np.eye(s)])
    Q, _ = np.linalg.qr(M)
    # Q[:,0] is parallel to ones; the next s-1 columns span the complement.
    B = Q[:, 1:s]
    Jt = B.T @ Js @ B

    eigs = np.linalg.eigvals(Jt)
    return np.asarray(eigs, dtype=complex)

def enumerate_replicator_rest_points(payoff_matrix: np.ndarray) -> list:
    A = np.asarray(payoff_matrix, dtype=float)
    K = A.shape[0]
    if A.shape != (K, K):
        raise ValueError("payoff_matrix must be square")

    results = []
    tol_pos = 1e-9

    indices = list(range(K))

    # iterate over all nonempty subsets
    for r in range(1, K + 1):
        for S in combinations(indices, r):
            S = tuple(sorted(S))
            if r == 1:
                i = S[0]
                pt = np.zeros(K)
                pt[i] = 1.0
                results.append({'support': S, 'point': pt})
                continue

            # Solve equal-payoff condition on support S.
            # Variables: y_i for i in S, with sum y = 1, and (A_S y)_i equal for all i in S.
            # Build linear system: (r-1) equations of (A_S y)_i - (A_S y)_{i0} = 0, plus sum=1.
            A_S = A[np.ix_(S, S)]
            r_S = len(S)
            # equations: M y = b
            # row k (k=0..r_S-2): A_S[k+1] y - A_S[0] y = 0
            # last row: sum y = 1
            M = np.zeros((r_S, r_S))
            b = np.zeros(r_S)
            for k in range(r_S - 1):
                M[k, :] = A_S[k + 1, :] - A_S[0, :]
            M[r_S - 1, :] = 1.0
            b[r_S - 1] = 1.0

            # Check rank to detect degenerate (non-isolated) faces.
            try:
                rank_M = np.linalg.matrix_rank(M, tol=1e-10)
            except np.linalg.LinAlgError:
                continue
            if rank_M < r_S:
                # degenerate / non-isolated; skip
                continue

            try:
                y = np.linalg.solve(M, b)
            except np.linalg.LinAlgError:
                continue

            # check strictly positive on support and nonnegative overall (zeros off support)
            if np.any(y <= tol_pos):
                continue
            if not np.isclose(np.sum(y), 1.0, atol=1e-9):
                continue

            pt = np.zeros(K)
            for idx, s in enumerate(S):
                pt[s] = y[idx]
            results.append({'support': S, 'point': pt})

    # Deduplicate: same point within 1e-9
    deduped = []
    for rp in results:
        is_dup = False
        for kept in deduped:
            if np.allclose(rp['point'], kept['point'], atol=1e-9):
                is_dup = True
                break
        if not is_dup:
            deduped.append(rp)

    return deduped

def analyze_mini_ultimatum(l: float, h: float, a: float) -> dict:
    def build_matrix(l, h, a):
        return construct_payoff_matrix(l, h, a)

    def replicator_jac_eigs(A, x):
        return replicator_jacobian_eigenvalues(A, x)

    def enumerate_rest_points(A):
        return enumerate_replicator_rest_points(A)

    def classify(eigs, tol=1e-9):
        reals = np.real(eigs)
        if np.any(reals > tol):
            if np.all(np.abs(reals) > tol) and np.any(reals < -tol):
                return 'unstable'
            return 'unstable'
        if np.all(reals < -tol):
            return 'stable'
        return 'nonhyperbolic'

    def stable_pure_set(l_, h_, a_):
        A_ = construct_payoff_matrix(l_, h_, a_)
        out = []
        for i in range(4):
            x = np.zeros(4); x[i] = 1.0
            eigs = replicator_jacobian_eigenvalues(A_, x)
            if classify(eigs) == 'stable':
                out.append(i)
        return sorted(out)

    # Build payoff matrix
    A = construct_payoff_matrix(l, h, a)

    # Enumerate rest points and classify
    rps = enumerate_replicator_rest_points(A)
    rest_points = []
    stable_pure = []
    for rp in rps:
        eigs = replicator_jacobian_eigenvalues(A, rp['point'])
        cls = classify(eigs)
        rest_points.append({
            'support': rp['support'],
            'point': rp['point'],
            'eigenvalues': eigs,
            'classification': cls,
        })
        if len(rp['support']) == 1 and cls == 'stable':
            stable_pure.append(rp['support'][0])
    stable_pure = sorted(stable_pure)

    # Find a_star: smallest a' in [0, h-l] where stable_pure_set differs from baseline at a'=0
    baseline = stable_pure_set(l, h, 0.0)
    a_max = h - l
    a_star = None
    if a_max > 0:
        # Scan to bracket a change, then bisect
        N = 400
        grid = np.linspace(0.0, a_max, N + 1)
        prev_set = baseline
        bracket = None
        for k in range(1, len(grid)):
            cur = stable_pure_set(l, h, float(grid[k]))
            if cur != baseline:
                bracket = (grid[k - 1], grid[k])
                break
        if bracket is not None:
            lo, hi = bracket
            # Bisection: find smallest a where set != baseline
            for _ in range(80):
                if hi - lo < 1e-7:
                    break
                mid = 0.5 * (lo + hi)
                if stable_pure_set(l, h, float(mid)) != baseline:
                    hi = mid
                else:
                    lo = mid
            a_star = float(hi)

    return {
        'payoff_matrix': A,
        'rest_points': rest_points,
        'stable_pure_strategies': stable_pure,
        'a_star': a_star,
    }

Prompt

Write a Python function that constructs the payoff matrix of the mini-ultimatum game with information effect. Strategies are indexed as $0 \to (h,h), \quad 1 \to (h,l), \quad 2 \to (l,h), \quad 3 \to (l,l),$ where the first coordinate is the proposer's offer and the second is the responder's minimum acceptable offer. The inputs satisfy $0 < l < h < 1/2$ and $0 \le a \le h - l$ . The total surplus is normalized to 1. Fitness convention. Entry M[i,j] is strategy i's fitness against strategy j, defined as the sum of i's payoff in the two role-symmetric encounters: (a) i as proposer vs j as responder, and (b) i as responder vs j as proposer. Acceptance rule. An offer $x \in \{l,h\}$ is accepted iff it is at least the responder's minimum acceptable offer. If accepted, the proposer earns $1-x$ and the responder earns $x$ ; if rejected both earn 0. Information effect. Whenever a high offer $h$ is made to a responder whose minimum acceptable offer is $l$ , the realized offer is reduced by $a$ before computing payoffs (only when accepted): the proposer earns $1-(h-a)$ and the responder earns $h-a$ .

def construct_payoff_matrix(l: float, h: float, a: float) -> np.ndarray:
        """
Construct the 4x4 payoff matrix for the mini-ultimatum game with information.

Parameters
----------
l : float
    Low offer/threshold level with 0 < l < h <= 1.
h : float
    High offer/threshold level.
a : float
    Information reduction applied when a high offer meets a low threshold;
    must satisfy 0 <= a <= h - l.

Returns
-------
np.ndarray of shape (4, 4)
    Fitness matrix M[i, j] for strategy i against strategy j under the
    role-symmetric mini-ultimatum rules.
"""
    return M

Background

The mini-ultimatum game is a discretized two-strategy-per-role version of the ultimatum bargaining game. A proposer offers a fraction of a unit surplus to a responder; the responder either accepts (both keep the agreed split) or rejects (both get 0). Restricting offers and minimum acceptable offers to two values $\{l,h\}$ produces four pure strategies, one for each (offer, threshold) combination. In an evolutionary / population context, players are randomly matched and play both roles symmetrically across encounters, so a strategy's overall fitness is the sum of its proposer-role and responder-role payoffs against the opponent. The payoff matrix is the input to replicator dynamics or other selection processes. The information parameter $a$ models a reduction in the realized offer when a proposer who plans to offer the high amount learns (or signals) that the responder would have accepted less. It is bounded by $h - l$ so the realized offer never falls below $l$ . This effect is the modeling device by which fairness norms can shift in equilibrium: as $a$ grows the high-offer strategies become less attractive against low-threshold responders.

Function signature

def construct_payoff_matrix(l: float, h: float, a: float) -> np.ndarray:
    return M

Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.

Test case 1 of 5

import numpy as np


def test_case_1():
    l, h, a = 0.1, 0.3, 0.0
    M = construct_payoff_matrix(l, h, a)
    assert M.shape == (4, 4)
    assert M.dtype == float
    assert abs(M[0, 0] - 1.0) < 1e-12
    assert abs(M[3, 3] - 1.0) < 1e-12
    assert np.all(M >= -1e-12)

Test case 2 of 5

import numpy as np


def test_case_2():
    l, h, a = 0.15, 0.45, 0.0
    M = construct_payoff_matrix(l, h, a)
    assert abs(M[2, 0] - h) < 1e-12
    assert abs(M[1, 0] - 1.0) < 1e-12
    assert abs(M[0, 2] - (1.0 - h)) < 1e-12

Test case 3 of 5

import numpy as np


def test_case_3():
    l, h, a = 0.2, 0.35, 0.05
    M = construct_payoff_matrix(l, h, a)
    assert abs(M[0, 3] - (1.0 - (h - a))) < 1e-12
    assert abs(M[3, 0] - (h - a)) < 1e-12
    assert abs(M[1, 2] - (1.0 - h + l)) < 1e-12

Test case 4 of 5

import numpy as np


def test_case_4():
    l, h = 0.1, 0.4
    a = h - l
    M = construct_payoff_matrix(l, h, a)
    assert abs(M[1, 1] - 1.0) < 1e-12

Test case 5 of 5

import numpy as np


def test_case_5():
    l, h, a = 0.1, 0.3, 0.0
    M = construct_payoff_matrix(l, h, a)
    assert isinstance(M, np.ndarray)
    assert M.dtype == float

Solution

import numpy as np

def construct_payoff_matrix(l: float, h: float, a: float) -> np.ndarray:
    # Strategy i -> (offer_i, threshold_i) with 0=(h,h),1=(h,l),2=(l,h),3=(l,l)
    strategies = [
        (h, h),
        (h, l),
        (l, h),
        (l, l),
    ]

    def one_way(offer, threshold_opp):
        # Payoff to proposer (with offer) facing responder with given threshold.
        # Returns (proposer_payoff, responder_payoff).
        realized = offer
        if offer == h and threshold_opp == l:
            realized = h - a
        # Use a small tolerance so the boundary case a == h - l (where the
        # realized offer equals l up to floating-point round-off) is accepted.
        if realized >= threshold_opp - 1e-12:
            return (1.0 - realized, realized)
        else:
            return (0.0, 0.0)

    M = np.zeros((4, 4), dtype=float)
    for i, (oi, ti) in enumerate(strategies):
        for j, (oj, tj) in enumerate(strategies):
            # i as proposer vs j as responder
            prop_i, _ = one_way(oi, tj)
            # j as proposer vs i as responder -> i gets responder payoff
            _, resp_i = one_way(oj, ti)
            M[i, j] = prop_i + resp_i
    return M

Prompt

Write a Python function replicator_jacobian_eigenvalues that computes the eigenvalues of the Jacobian of the replicator vector field at an arbitrary point on the simplex, projected onto the tangent space of the simplex on the support of that point. The replicator dynamics on K strategies with payoff matrix A and frequency vector x is dx_i/dt = x_i * ( (A x)_i - x^T A x ). Given a point x on the K-simplex (nonnegative, sums to 1) with support S = {i : x_i > 0}, your function must:

For a vertex (|S| = 1, focal index f): return the K-1 transversal invasion eigenvalues A[j,f] - A[f,f] for j != f.
For a non-vertex x (|S| >= 2): build the full Jacobian of the replicator field restricted to rows/columns in S, project it onto the tangent space of the simplex on S (the orthogonal complement of the all-ones vector on S), and return the |S|-1 eigenvalues of that projected operator. The off-diagonal Jacobian entries must be honored — do not use the diagonal-only vertex shortcut.

def replicator_jacobian_eigenvalues(payoff_matrix: np.ndarray, x: np.ndarray) -> np.ndarray:
            """
Compute eigenvalues of the replicator dynamics Jacobian at state x.

Parameters
----------
payoff_matrix : np.ndarray of shape (K, K)
    Payoff matrix A for K strategies.
x : np.ndarray of shape (K,)
    Population share vector on the simplex.

Returns
-------
np.ndarray of shape (K,)
    Complex eigenvalues of the Jacobian restricted to the tangent simplex.
"""
    return np.asarray(eigs, dtype=complex)

Background

The replicator equation is the canonical selection dynamic of evolutionary game theory. It evolves a vector of strategy frequencies on the simplex, where each component grows in proportion to how its expected payoff compares to the population mean payoff. Local stability of a rest point is governed by the linearization of this vector field. Because the simplex is a (K-1)-dimensional manifold embedded in K-space (the constraint sum_i x_i = 1 is preserved by the dynamics), the relevant spectrum is that of the Jacobian projected onto the tangent space of the simplex — there is always a trivial direction normal to the simplex that contributes a redundant eigenvalue and must be removed. At a vertex of the simplex (a pure-strategy state) the projected linearization decouples cleanly and reduces to a set of one-dimensional invasion-fitness comparisons: each non-focal strategy contributes one eigenvalue measuring whether a small invasion by that strategy grows or decays. At a mixed point — on an edge, face, or in the interior — the Jacobian generically has nonzero off-diagonal entries within the support, so the diagonal-only vertex argument no longer applies and one must work with the full restricted matrix and explicitly project out the simplex-normal direction on the support. The sign pattern of the real parts of the resulting eigenvalues classifies the rest point as a sink, source, or saddle within its face, which is the input the surrounding bifurcation analysis needs.

Function signature

def replicator_jacobian_eigenvalues(payoff_matrix: np.ndarray, x: np.ndarray) -> np.ndarray:
    return np.asarray(eigs, dtype=complex)

Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.

Test case 1 of 4

import numpy as np


def test_case_1():
    A = np.array([[2.0, 1.0, 0.0],
                  [3.0, 2.0, 1.0],
                  [1.0, 0.0, 2.0]])
    K = 3
    for f in range(K):
        x = np.zeros(K)
        x[f] = 1.0
        eigs = replicator_jacobian_eigenvalues(A, x)
        assert eigs.shape == (K - 1,)
        expected = np.array([A[j, f] - A[f, f] for j in range(K) if j != f], dtype=complex)
        assert np.allclose(np.sort_complex(eigs), np.sort_complex(expected))

Test case 2 of 4

import numpy as np


def test_case_2():
    a, b = 3.0, 2.0
    A2 = np.array([[a, 0.0],
                   [0.0, b]])
    p = b / (a + b)
    A = np.zeros((3, 3))
    A[:2, :2] = A2
    x = np.array([p, 1.0 - p, 0.0])
    eigs = replicator_jacobian_eigenvalues(A, x)
    assert eigs.shape == (1,)
    assert eigs[0].real > 1e-9
    expected = p * (1.0 - p) * (a + b)
    assert abs(eigs[0].real - expected) < 1e-8
    assert abs(eigs[0].imag) < 1e-8

Test case 3 of 4

import numpy as np


def test_case_3():
    rng = np.random.default_rng(0)
    K = 4
    A = rng.standard_normal((K, K))
    x = np.array([0.2, 0.5, 0.3, 0.0])
    def f(y):
        Ay = A @ y
        phi = float(y @ Ay)
        return y * (Ay - phi)
    eps = 1e-6
    J_fd = np.zeros((K, K))
    for k in range(K):
        e = np.zeros(K); e[k] = eps
        J_fd[:, k] = (f(x + e) - f(x - e)) / (2 * eps)
    support = [0, 1, 2]
    S = np.array(support)
    Js = J_fd[np.ix_(S, S)]
    s = len(S)
    ones = np.ones(s) / np.sqrt(s)
    M = np.column_stack([ones, np.eye(s)])
    Q, _ = np.linalg.qr(M)
    B = Q[:, 1:s]
    Jt = B.T @ Js @ B
    eigs_fd = np.sort_complex(np.linalg.eigvals(Jt))
    eigs = np.sort_complex(replicator_jacobian_eigenvalues(A, x))
    assert eigs.shape == (2,)
    assert np.allclose(eigs, eigs_fd, atol=1e-4)

Test case 4 of 4

import numpy as np


def test_case_4():
    A = np.array([[2.0, 1.0], [1.0, 2.0]])
    x = np.array([1.0, 0.0])
    eigs = replicator_jacobian_eigenvalues(A, x)
    assert isinstance(eigs, np.ndarray)
    assert np.issubdtype(eigs.dtype, np.complexfloating)

Solution

import numpy as np

def replicator_jacobian_eigenvalues(payoff_matrix: np.ndarray, x: np.ndarray) -> np.ndarray:
    A = np.asarray(payoff_matrix, dtype=float)
    x = np.asarray(x, dtype=float)
    K = A.shape[0]
    if A.shape != (K, K):
        raise ValueError("payoff_matrix must be square")
    if x.shape != (K,):
        raise ValueError("x must have shape (K,)")
    if not np.all(x >= -1e-12):
        raise ValueError("x must be nonnegative")
    if abs(x.sum() - 1.0) > 1e-10:
        raise ValueError("x must sum to 1")

    support = [i for i in range(K) if x[i] > 1e-12]

    # Vertex case: return K-1 transversal invasion eigenvalues.
    if len(support) == 1:
        f = support[0]
        eigs = [A[j, f] - A[f, f] for j in range(K) if j != f]
        return np.array(eigs, dtype=complex)

    # Mixed point: build full Jacobian, restrict to support, project to tangent space.
    # f_i(x) = x_i * ((A x)_i - x^T A x)
    # df_i/dx_k = delta_{ik} * ((A x)_i - phi) + x_i * (A_{ik} - (x^T A x)_k')
    # where phi = x^T A x and d(phi)/dx_k = (A x)_k + (A^T x)_k = ((A + A^T) x)_k.
    Ax = A @ x
    phi = float(x @ Ax)
    dphi = (A + A.T) @ x  # shape (K,)

    # Full K x K Jacobian
    J = np.zeros((K, K), dtype=float)
    for i in range(K):
        for k in range(K):
            term = x[i] * (A[i, k] - dphi[k])
            if i == k:
                term += (Ax[i] - phi)
            J[i, k] = term

    # Restrict to support
    S = np.array(support, dtype=int)
    Js = J[np.ix_(S, S)]
    s = len(S)

    # Project onto tangent space of simplex on S (orthogonal complement of
    # the all-ones vector in R^s). Build an orthonormal basis B of shape
    # (s, s-1) for that complement and form B^T J_s B.
    ones = np.ones(s) / np.sqrt(s)
    # Use SVD/QR of the orthogonal complement.
    # Construct s x s matrix [ones | I] and orthonormalize, then drop the first column.
    M = np.column_stack([ones, np.eye(s)])
    Q, _ = np.linalg.qr(M)
    # Q[:,0] is parallel to ones; the next s-1 columns span the complement.
    B = Q[:, 1:s]
    Jt = B.T @ Js @ B

    eigs = np.linalg.eigvals(Jt)
    return np.asarray(eigs, dtype=complex)

Enumerating rest points of replicator dynamics on the simplex

Given a K x K payoff matrix A, the replicator equation on the K-simplex is $\dot{x}_i = x_i \big((Ax)_i - x^\top A x\big).$ A rest point is any state x on the simplex at which dx/dt = 0. These come in three flavors: pure vertices, mixed equilibria living strictly inside a face, and the fully interior equilibrium (when one exists). Write a function that enumerates every isolated rest point lying strictly inside its face.

def enumerate_replicator_rest_points(payoff_matrix: np.ndarray) -> list:
        """
Enumerate isolated rest points of replicator dynamics on the K-simplex.

Parameters
----------
payoff_matrix : np.ndarray of shape (K, K)
    Square payoff matrix, dtype float.

Returns
-------
list of dict
    Each dict has keys 'support' (tuple of active indices) and 'point'
    (nonnegative probability vector summing to 1). Includes all vertex rest
    points plus isolated mixed equilibria; duplicates are removed.
"""
    return deduped

Background

Replicator dynamics and rest points

The replicator equation is the canonical model of frequency-dependent selection in evolutionary game theory. A population is described by a state vector x on the K-simplex, where x_i is the proportion of the population using strategy i. Strategies that earn above-average payoff grow in frequency; those below average shrink. A rest point (equilibrium) of the dynamics is a population state at which all per-capita growth rates vanish. Rest points fall on the boundary or interior of the simplex according to which strategies are present. A pure vertex (only one strategy present) is automatically a rest point because the dynamics preserves the support: a strategy at frequency zero stays at zero. Mixed rest points require an additional condition — every strategy that is present must earn the same payoff against the current mixed population, since otherwise the higher-earning ones would grow at the expense of the lower-earning ones. Rest points may live on edges (two strategies present), on higher-dimensional faces, or strictly inside the simplex. The full set of rest points is what one needs to draw the qualitative phase portrait of the dynamics: stable rest points are attractors, unstable ones are repellers or saddles, and heteroclinic connections run between them. For generic payoff matrices each face contributes at most one isolated rest point, but special (degenerate) matrices can produce continuous families of equilibria along a face — for example, a face on which two strategies always earn identical payoffs. Such degenerate faces are typically reported separately or skipped in equilibrium enumerations.

Function signature

def enumerate_replicator_rest_points(payoff_matrix: np.ndarray) -> list:
    return deduped

Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.

Test case 1 of 3

import numpy as np


def test_case_1():
    # Vertices always present: identity payoff matrix, only vertices should be rest points
    # (interior would require equal payoff (Ax)_i const => x_i const => uniform, but check it)
    A = np.eye(3)
    rps = enumerate_replicator_rest_points(A)
    # must contain all 3 vertices
    supports = [rp['support'] for rp in rps]
    for i in range(3):
        assert (i,) in supports, f"vertex {i} missing"
    # each point must be on the simplex
    for rp in rps:
        pt = rp['point']
        assert pt.shape == (3,)
        assert np.all(pt >= -1e-12), "negative entry"
        assert np.isclose(pt.sum(), 1.0), "does not sum to 1"
        # zeros outside support
        for j in range(3):
            if j not in rp['support']:
                assert abs(pt[j]) < 1e-12, "nonzero outside support"

import numpy as np

Test case 2 of 3

import numpy as np


def test_case_2():
    # Rock-Paper-Scissors: A = [[0,-1,1],[1,0,-1],[-1,1,0]] has unique interior equilibrium (1/3,1/3,1/3)
    A = np.array([[0.0, -1.0, 1.0], [1.0, 0.0, -1.0], [-1.0, 1.0, 0.0]])
    rps = enumerate_replicator_rest_points(A)
    # find interior equilibrium
    interior = [rp for rp in rps if rp['support'] == (0, 1, 2)]
    assert len(interior) == 1, f"expected 1 interior equilibrium, got {len(interior)}"
    pt = interior[0]['point']
    expected = np.array([1/3, 1/3, 1/3])
    assert np.allclose(pt, expected, atol=1e-9), f"interior point {pt} != {expected}"
    # No edge equilibria for RPS (each edge has a dominant strategy)
    for rp in rps:
        if len(rp['support']) == 2:
            # verify by checking whether equal-payoff condition could hold strictly inside edge
            # In RPS each pair has cyclic dominance => no interior edge eq. So this branch should be empty.
            assert False, f"unexpected edge equilibrium on {rp['support']}"

import numpy as np

Test case 3 of 3

import numpy as np


def test_case_3():
    # 2x2 coordination-like game: A = [[2,0],[0,1]]
    # Equal payoff on edge {0,1}: 2 x_0 = x_1, with x_0 + x_1 = 1 => x_0 = 1/3, x_1 = 2/3
    A = np.array([[2.0, 0.0], [0.0, 1.0]])
    rps = enumerate_replicator_rest_points(A)
    supports = [rp['support'] for rp in rps]
    assert (0,) in supports and (1,) in supports, "vertices missing"
    mixed = [rp for rp in rps if rp['support'] == (0, 1)]
    assert len(mixed) == 1, f"expected 1 mixed eq on edge, got {len(mixed)}"
    pt = mixed[0]['point']
    # solve manually: (Ax)_0 = 2 x_0, (Ax)_1 = x_1, equal => 2 x_0 = x_1, sum=1 => x_0=1/3
    expected = np.array([1/3, 2/3])
    assert np.allclose(pt, expected, atol=1e-9), f"mixed point {pt} != {expected}"
    # all points on simplex
    for rp in rps:
        assert np.isclose(rp['point'].sum(), 1.0), "not on simplex"
        assert np.all(rp['point'] >= -1e-12), "negative entry"

Solution

import numpy as np
from itertools import combinations


def enumerate_replicator_rest_points(payoff_matrix: np.ndarray) -> list:
    A = np.asarray(payoff_matrix, dtype=float)
    K = A.shape[0]

    results = []
    tol_pos = 1e-9

    indices = list(range(K))

    # iterate over all nonempty subsets
    for r in range(1, K + 1):
        for S in combinations(indices, r):
            S = tuple(sorted(S))
            if r == 1:
                i = S[0]
                pt = np.zeros(K)
                pt[i] = 1.0
                results.append({'support': S, 'point': pt})
                continue

            # Solve equal-payoff condition on support S.
            # Variables: y_i for i in S, with sum y = 1, and (A_S y)_i equal for all i in S.
            # Build linear system: (r-1) equations of (A_S y)_i - (A_S y)_{i0} = 0, plus sum=1.
            A_S = A[np.ix_(S, S)]
            r_S = len(S)
            # equations: M y = b
            # row k (k=0..r_S-2): A_S[k+1] y - A_S[0] y = 0
            # last row: sum y = 1
            M = np.zeros((r_S, r_S))
            b = np.zeros(r_S)
            for k in range(r_S - 1):
                M[k, :] = A_S[k + 1, :] - A_S[0, :]
            M[r_S - 1, :] = 1.0
            b[r_S - 1] = 1.0

            # Check rank to detect degenerate (non-isolated) faces.
            try:
                rank_M = np.linalg.matrix_rank(M, tol=1e-10)
            except np.linalg.LinAlgError:
                continue
            if rank_M < r_S:
                # degenerate / non-isolated; skip
                continue

            try:
                y = np.linalg.solve(M, b)
            except np.linalg.LinAlgError:
                continue

            # check strictly positive on support and nonnegative overall (zeros off support)
            if np.any(y <= tol_pos):
                continue
            if not np.isclose(np.sum(y), 1.0, atol=1e-9):
                continue

            pt = np.zeros(K)
            for idx, s in enumerate(S):
                pt[s] = y[idx]
            results.append({'support': S, 'point': pt})

    # Deduplicate: same point within 1e-9
    deduped = []
    for rp in results:
        is_dup = False
        for kept in deduped:
            if np.allclose(rp['point'], kept['point'], atol=1e-9):
                is_dup = True
                break
        if not is_dup:
            deduped.append(rp)

    return deduped