Assignment Algorithms

Algorithms for data association and assignment problems in multi-target tracking.

Assignment algorithms for data association in target tracking.

This module provides: - 2D assignment algorithms (Hungarian, Auction) - K-best 2D assignment (Murty’s algorithm) - 3D assignment algorithms (Lagrangian relaxation, Auction) - Gating methods (ellipsoidal, rectangular) - Data association algorithms (GNN, JPDA)

2D Assignment

Optimal assignment algorithms for bipartite matching problems.

Hungarian Algorithm

Gating

Validation region (gating) functions for measurement association.

Gating functions for data association in target tracking.

This module provides gating methods to determine which measurements fall within a validation region around predicted track states.

pytcl.assignment_algorithms.gating.mahalanobis_distance(innovation, innovation_covariance)[source]

Compute the squared Mahalanobis distance.

The Mahalanobis distance measures how many standard deviations a point is from the center of a distribution.

Parameters:
  • innovation (array_like) – Innovation (measurement residual) vector of shape (m,).

  • innovation_covariance (array_like) – Innovation covariance matrix of shape (m, m).

Returns:

Squared Mahalanobis distance.

Return type:

float

Examples

>>> innovation = np.array([1.0, 0.5])
>>> S = np.array([[2.0, 0.0], [0.0, 1.0]])
>>> d2 = mahalanobis_distance(innovation, S)
>>> d2
0.75

Notes

The squared Mahalanobis distance is defined as:

d^2 = (z - z_pred)^T @ S^{-1} @ (z - z_pred)

where S is the innovation covariance matrix.

Dispatch (measured; see task-C2-report.md in the v2.5.0-region-lcd-perf campaign for the full before/after numbers and the behavior-equality bounds):

  • n == 2 or n == 3: closed-form njit kernels (_mahalanobis_distance_2d/ _3d, fed by _invert_2x2/_invert_3x3) that avoid np.linalg.solve’s generic LAPACK dispatch entirely (measured ~4x faster than the generic path at these sizes).

  • 1 <= n <= 10 otherwise: np.linalg.inv + the njit _mahalanobis_distance_general quadratic-form kernel, measured ~10-25% faster than np.linalg.solve in this range (a single inv call plus a tight njit loop beats solve’s per-call LAPACK dispatch overhead here; the constant-overhead advantage disappears and reverses for n > 12, where inv’s O(n^3) cost with a larger constant than solving one right-hand side starts to dominate – hence the n <= 10 cutoff, comfortably inside the measured-winning range).

  • n > 10, or an exactly-singular 2D/3D covariance: the original generic np.linalg.solve path. Every branch raises the same numpy.linalg.LinAlgError on an exactly singular covariance.

pytcl.assignment_algorithms.gating.mahalanobis_batch(innovations, S_inv, output)[source]

Compute Mahalanobis distances for a batch of innovations.

JIT-compiled for performance. Computes squared Mahalanobis distances for multiple innovations against a single covariance matrix.

Parameters:
  • innovations (ndarray) – Innovations of shape (n_measurements, dim).

  • S_inv (ndarray) – Inverse of innovation covariance matrix of shape (dim, dim).

  • output (ndarray) – Output array of shape (n_measurements,) to store distances.

pytcl.assignment_algorithms.gating.ellipsoidal_gate(innovation, innovation_covariance, gate_threshold)[source]

Test if a measurement passes an ellipsoidal gate.

The ellipsoidal gate defines a validation region based on the chi-squared distribution of the squared Mahalanobis distance.

Parameters:
  • innovation (array_like) – Innovation vector of shape (m,).

  • innovation_covariance (array_like) – Innovation covariance matrix of shape (m, m).

  • gate_threshold (float) – Gate threshold (chi-squared value). Common values: - 9.21 for 99% probability with 2 measurements - 11.34 for 99% probability with 3 measurements - 16.27 for 99% probability with 4 measurements

Returns:

True if measurement passes the gate (is inside the ellipsoid).

Return type:

bool

Examples

>>> innovation = np.array([1.0, 0.5])
>>> S = np.array([[2.0, 0.0], [0.0, 1.0]])
>>> ellipsoidal_gate(innovation, S, gate_threshold=9.21)
True

See also

chi2_gate_threshold

Compute threshold from probability.

pytcl.assignment_algorithms.gating.rectangular_gate(innovation, innovation_covariance, num_sigmas=3.0)[source]

Test if a measurement passes a rectangular gate.

The rectangular gate defines a validation region as a hypercube based on the marginal standard deviations.

Parameters:
  • innovation (array_like) – Innovation vector of shape (m,).

  • innovation_covariance (array_like) – Innovation covariance matrix of shape (m, m).

  • num_sigmas (float, optional) – Number of standard deviations for gate bounds (default: 3.0).

Returns:

True if measurement passes the gate.

Return type:

bool

Examples

>>> innovation = np.array([1.0, 0.5])
>>> S = np.array([[4.0, 0.0], [0.0, 1.0]])
>>> rectangular_gate(innovation, S, num_sigmas=3.0)
True

Notes

Rectangular gating is computationally cheaper but less tight than ellipsoidal gating. It may pass more false measurements.

pytcl.assignment_algorithms.gating.gate_measurements(predicted_measurement, innovation_covariance, measurements, gate_threshold, gate_type='ellipsoidal')[source]

Gate multiple measurements against a predicted track state.

Parameters:
  • predicted_measurement (array_like) – Predicted measurement of shape (m,).

  • innovation_covariance (array_like) – Innovation covariance matrix of shape (m, m).

  • measurements (array_like) – Array of measurements of shape (n_meas, m).

  • gate_threshold (float) – Gate threshold. For ellipsoidal gates, this is the chi-squared value. For rectangular gates, this is the number of sigmas.

  • gate_type (str, optional) – Type of gate: “ellipsoidal” or “rectangular” (default: “ellipsoidal”).

Returns:

  • valid_indices (ndarray) – Indices of measurements that pass the gate.

  • distances (ndarray) – Squared Mahalanobis distances for valid measurements.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[int64]], ndarray[tuple[Any, …], dtype[float64]]]

Examples

>>> z_pred = np.array([0.0, 0.0])
>>> S = np.eye(2)
>>> measurements = np.array([[0.5, 0.5], [5.0, 5.0], [1.0, -1.0]])
>>> valid_idx, dists = gate_measurements(z_pred, S, measurements, 9.21)
>>> valid_idx
array([0, 2])

Notes

This function efficiently gates multiple measurements against a single track prediction, which is common in multi-target tracking.

pytcl.assignment_algorithms.gating.chi2_gate_threshold(probability, num_dimensions)[source]

Compute chi-squared gate threshold for a given probability.

Parameters:
  • probability (float) – Gate probability (e.g., 0.99 for 99% of true measurements to pass).

  • num_dimensions (int) – Measurement dimension (degrees of freedom).

Returns:

Chi-squared threshold value.

Return type:

float

Examples

>>> round(chi2_gate_threshold(0.99, 2), 6)  # 2D measurement, 99% probability
9.21034
>>> round(chi2_gate_threshold(0.99, 3), 6)  # 3D measurement, 99% probability
11.344867
pytcl.assignment_algorithms.gating.compute_gate_volume(innovation_covariance, gate_threshold)[source]

Compute the volume of an ellipsoidal gate.

Parameters:
  • innovation_covariance (array_like) – Innovation covariance matrix of shape (m, m).

  • gate_threshold (float) – Chi-squared gate threshold.

Returns:

Volume of the ellipsoidal gate region.

Return type:

float

Notes

The gate volume is used in probabilistic data association methods to compute the clutter density.

For an m-dimensional ellipsoid with threshold gamma:

V = c_m * sqrt(det(S)) * gamma^(m/2)

where c_m is the volume of the unit hypersphere in m dimensions.

Examples

Compute gate volume for a 2D measurement with 99% gate probability:

>>> import numpy as np
>>> from scipy.stats import chi2
>>> S = np.array([[4.0, 0.0], [0.0, 1.0]])  # innovation covariance
>>> gate_prob = 0.99
>>> threshold = chi2.ppf(gate_prob, df=2)
>>> volume = compute_gate_volume(S, threshold)
>>> volume > 0
True

See also

ellipsoidal_gate

Test if measurement passes gate.

mahalanobis_distance

Compute distance used in gating.

Data Association

Track-to-measurement association algorithms.

Global Nearest Neighbor (GNN)

Data association algorithms for multi-target tracking.

This module provides algorithms for associating measurements to tracks, including Global Nearest Neighbor (GNN) and related methods.

class pytcl.assignment_algorithms.data_association.AssociationResult(track_to_measurement, measurement_to_track, costs, total_cost)[source]

Bases: NamedTuple

Result of data association.

Variables:
  • track_to_measurement (ndarray) – Array of shape (n_tracks,). track_to_measurement[i] gives the measurement index assigned to track i, or -1 if unassigned.

  • measurement_to_track (ndarray) – Array of shape (n_measurements,). measurement_to_track[j] gives the track index assigned to measurement j, or -1 if unassigned.

  • costs (ndarray) – Array of association costs (squared Mahalanobis distances) for each assigned pair.

  • total_cost (float) – Total assignment cost.

track_to_measurement: ndarray[tuple[Any, ...], dtype[int64]]

Alias for field number 0

measurement_to_track: ndarray[tuple[Any, ...], dtype[int64]]

Alias for field number 1

costs: ndarray[tuple[Any, ...], dtype[float64]]

Alias for field number 2

total_cost: float

Alias for field number 3

pytcl.assignment_algorithms.data_association.compute_association_cost(track_predictions, track_covariances, measurements, measurement_models=None)[source]

Compute cost matrix for track-to-measurement association.

Parameters:
  • track_predictions (array_like) – Predicted track states of shape (n_tracks, n_state).

  • track_covariances (array_like) – Track covariance matrices of shape (n_tracks, n_state, n_state).

  • measurements (array_like) – Measurements of shape (n_measurements, n_meas).

  • measurement_models (array_like, optional) – Measurement matrices of shape (n_tracks, n_meas, n_state) or (n_meas, n_state) if same for all tracks. If None, assumes direct measurement of first n_meas states.

Returns:

cost_matrix – Cost matrix of shape (n_tracks, n_measurements). cost_matrix[i, j] is the squared Mahalanobis distance from track i to measurement j.

Return type:

ndarray

Examples

>>> # 2 tracks, 3 measurements, 2D state [x, vx]
>>> predictions = np.array([[0.0, 1.0], [5.0, -1.0]])
>>> covariances = np.array([np.eye(2), np.eye(2)])
>>> measurements = np.array([[0.1], [4.9], [10.0]])
>>> H = np.array([[1.0, 0.0]])  # Measure position only
>>> costs = compute_association_cost(predictions, covariances, measurements, H)
pytcl.assignment_algorithms.data_association.nearest_neighbor(cost_matrix, gate_threshold=inf)[source]

Simple nearest neighbor data association.

Each track is assigned to its closest measurement, but a measurement can only be assigned to one track. Greedy assignment, not globally optimal.

Parameters:
  • cost_matrix (array_like) – Cost matrix of shape (n_tracks, n_measurements).

  • gate_threshold (float, optional) – Maximum cost for valid assignment. Assignments with cost above this threshold are rejected.

Returns:

result – Association result with track and measurement assignments.

Return type:

AssociationResult

Notes

This is a greedy algorithm that processes tracks in order of their minimum cost. It does not guarantee globally optimal assignment. Use gnn_association for globally optimal assignment.

Examples

>>> cost = np.array([[1.0, 5.0], [4.0, 2.0]])
>>> result = nearest_neighbor(cost, gate_threshold=10.0)
>>> result.track_to_measurement
array([0, 1])
pytcl.assignment_algorithms.data_association.gnn_association(cost_matrix, gate_threshold=inf, cost_of_non_assignment=None)[source]

Global Nearest Neighbor (GNN) data association.

Finds the globally optimal assignment of tracks to measurements that minimizes the total association cost.

Parameters:
  • cost_matrix (array_like) – Cost matrix of shape (n_tracks, n_measurements). cost_matrix[i, j] is the cost of assigning track i to measurement j.

  • gate_threshold (float, optional) – Maximum cost for valid assignment. Entries above this threshold are set to infinity before optimization.

  • cost_of_non_assignment (float, optional) – Cost for not assigning a track or measurement. If None, all tracks/measurements that can be assigned will be assigned.

Returns:

result – Association result with track and measurement assignments.

Return type:

AssociationResult

Examples

>>> cost = np.array([[1.0, 5.0, 2.0],
...                  [4.0, 2.0, 3.0]])
>>> result = gnn_association(cost, gate_threshold=10.0)
>>> result.track_to_measurement
array([0, 1])
>>> result.total_cost
3.0

Notes

GNN uses the Hungarian algorithm to find the globally optimal assignment. It is more computationally expensive than nearest neighbor but guarantees optimality.

The algorithm handles rectangular cost matrices (different numbers of tracks and measurements) and allows tracks/measurements to remain unassigned if cost_of_non_assignment is specified.

References

  • Blackman, S.S. and Popoli, R., “Design and Analysis of Modern Tracking Systems”, Artech House, 1999.

pytcl.assignment_algorithms.data_association.gated_gnn_association(track_predictions, track_covariances, measurements, measurement_models=None, gate_probability=0.99, cost_of_non_assignment=None)[source]

GNN association with automatic gating.

Combines gating and GNN association in a single function for convenience.

Parameters:
  • track_predictions (array_like) – Predicted track states of shape (n_tracks, n_state).

  • track_covariances (array_like) – Track covariance matrices of shape (n_tracks, n_state, n_state).

  • measurements (array_like) – Measurements of shape (n_measurements, n_meas).

  • measurement_models (array_like, optional) – Measurement matrices. See compute_association_cost for details.

  • gate_probability (float, optional) – Probability for chi-squared gate threshold (default: 0.99).

  • cost_of_non_assignment (float, optional) – Cost for not assigning a track or measurement.

Returns:

result – Association result.

Return type:

AssociationResult

Examples

>>> predictions = np.array([[0.0, 1.0], [5.0, -1.0]])
>>> covariances = np.array([0.1 * np.eye(2), 0.1 * np.eye(2)])
>>> measurements = np.array([[0.1], [4.9]])
>>> H = np.array([[1.0, 0.0]])
>>> result = gated_gnn_association(
...     predictions, covariances, measurements, H,
...     gate_probability=0.99
... )

Joint Probabilistic Data Association (JPDA)

The JPDA algorithm computes association probabilities for all feasible track-measurement pairings and updates each track with a weighted combination of innovations.

Joint Probabilistic Data Association (JPDA) algorithm.

JPDA computes association probabilities between tracks and measurements by considering all possible joint association hypotheses, then performs a probabilistically weighted update for each track.

This is more sophisticated than GNN which makes hard assignment decisions, as JPDA can handle measurement origin uncertainty in cluttered environments.

class pytcl.assignment_algorithms.jpda.JPDAResult(association_probs, marginal_probs, likelihood_matrix, gated)[source]

Bases: NamedTuple

Result of JPDA algorithm.

Variables:
  • association_probs (ndarray[Any]) – Association probability matrix of shape (n_tracks, n_measurements + 1). association_probs[i, j] is the probability that track i is associated with measurement j. The last column (j = n_measurements) represents the probability that track i has no measurement.

  • marginal_probs (list of ndarray) – List of marginal association probabilities for each track. marginal_probs[i][j] = P(measurement j originated from track i).

  • likelihood_matrix (ndarray[Any]) – Measurement likelihood matrix of shape (n_tracks, n_measurements).

  • gated (ndarray[Any]) – Boolean matrix indicating which track-measurement pairs passed gating.

association_probs: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

marginal_probs: List[ndarray[tuple[Any, ...], dtype[floating]]]

Alias for field number 1

likelihood_matrix: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

gated: ndarray[tuple[Any, ...], dtype[bool]]

Alias for field number 3

class pytcl.assignment_algorithms.jpda.JPDAUpdate(states, covariances, association_probs, innovations)[source]

Bases: NamedTuple

Result of JPDA-based track update.

Variables:
  • states (list of ndarray) – Updated state estimates for each track.

  • covariances (list of ndarray) – Updated covariances for each track (includes spread of means).

  • association_probs (ndarray[Any]) – Association probability matrix.

  • innovations (list of ndarray) – Combined weighted innovations for each track.

states: List[ndarray[tuple[Any, ...], dtype[floating]]]

Alias for field number 0

covariances: List[ndarray[tuple[Any, ...], dtype[floating]]]

Alias for field number 1

association_probs: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 2

innovations: List[ndarray[tuple[Any, ...], dtype[floating]]]

Alias for field number 3

pytcl.assignment_algorithms.jpda.compute_measurement_likelihood(innovation, innovation_cov, detection_prob=1.0)[source]

Compute measurement likelihood for a track-measurement pair.

Parameters:
  • innovation (ndarray[Any]) – Measurement innovation (residual), shape (m,).

  • innovation_cov (ndarray[Any]) – Innovation covariance, shape (m, m).

  • detection_prob (float) – Probability of detection (Pd).

Returns:

likelihood – Measurement likelihood.

Return type:

float

pytcl.assignment_algorithms.jpda.compute_likelihood_matrix(track_states, track_covariances, measurements, H, R, detection_prob=1.0, gate_threshold=None)[source]

Compute likelihood matrix for all track-measurement pairs.

Parameters:
  • track_states (list of ndarray) – State estimates for each track.

  • track_covariances (list of ndarray) – Covariances for each track.

  • measurements (ndarray[Any]) – Measurements, shape (n_meas, m).

  • H (ndarray[Any]) – Measurement matrix, shape (m, n).

  • R (ndarray[Any]) – Measurement noise covariance, shape (m, m).

  • detection_prob (float) – Probability of detection.

  • gate_threshold (float, optional) – Mahalanobis distance threshold for gating.

Returns:

  • likelihood_matrix (ndarray[Any]) – Likelihood values, shape (n_tracks, n_meas).

  • gated (ndarray[Any]) – Boolean gating matrix, shape (n_tracks, n_meas).

Return type:

tuple[ndarray[tuple[Any, …], dtype[Any]], ndarray[tuple[Any, …], dtype[Any]]]

Examples

>>> import numpy as np
>>> # Two tracks, two measurements
>>> states = [np.array([0.0, 1.0]), np.array([5.0, 0.0])]
>>> covs = [np.eye(2) * 0.5, np.eye(2) * 0.5]
>>> measurements = np.array([[0.1], [5.2]])
>>> H = np.array([[1, 0]])
>>> R = np.array([[0.1]])
>>> L, gated = compute_likelihood_matrix(states, covs, measurements, H, R)
>>> L.shape
(2, 2)
>>> # Track 0 should have high likelihood for measurement 0
>>> L[0, 0] > L[0, 1]
True
pytcl.assignment_algorithms.jpda.jpda_probabilities(likelihood_matrix, gated, detection_prob=1.0, clutter_density=1e-06)[source]

Compute JPDA association probabilities.

Uses an efficient algorithm that avoids explicit enumeration of all joint hypotheses when possible.

Parameters:
  • likelihood_matrix (ndarray[Any]) – Likelihood values, shape (n_tracks, n_meas).

  • gated (ndarray[Any]) – Boolean gating matrix, shape (n_tracks, n_meas).

  • detection_prob (float) – Probability of detection (Pd).

  • clutter_density (float) – Spatial density of clutter (lambda_c).

Returns:

beta – Association probability matrix, shape (n_tracks, n_meas + 1). beta[i, j] = P(measurement j is from track i) for j < n_meas. beta[i, n_meas] = P(track i has no measurement).

Return type:

ndarray[Any]

Examples

>>> import numpy as np
>>> # Likelihood matrix: 2 tracks, 2 measurements
>>> # Track 0 has high likelihood for meas 0
>>> # Track 1 has high likelihood for meas 1
>>> likelihood = np.array([[0.9, 0.1],
...                        [0.1, 0.8]])
>>> gated = np.array([[True, True],
...                   [True, True]])
>>> beta = jpda_probabilities(likelihood, gated, detection_prob=0.9)
>>> beta.shape  # 2 tracks, 3 columns (2 meas + 1 miss)
(2, 3)
>>> # Track 0 most likely associated with measurement 0
>>> np.argmax(beta[0, :2])
0
pytcl.assignment_algorithms.jpda.jpda_update(track_states, track_covariances, measurements, H, R, detection_prob=0.9, clutter_density=1e-06, gate_probability=0.99)[source]

Perform JPDA-based track update.

Parameters:
  • track_states (list of array_like) – Predicted state estimates for each track.

  • track_covariances (list of array_like) – Predicted covariances for each track.

  • measurements (array_like) – Measurements, shape (n_meas, m).

  • H (array_like) – Measurement matrix, shape (m, n).

  • R (array_like) – Measurement noise covariance, shape (m, m).

  • detection_prob (float) – Probability of detection (Pd). Default 0.9.

  • clutter_density (float) – Spatial density of clutter. Default 1e-6.

  • gate_probability (float) – Probability for chi-squared gate. Default 0.99.

Returns:

result – Updated states and covariances with association probabilities.

Return type:

JPDAUpdate

Examples

>>> import numpy as np
>>> # Two tracks, three measurements
>>> x1 = np.array([0., 1.])  # [position, velocity]
>>> x2 = np.array([5., -1.])
>>> P = np.eye(2) * 0.1
>>> measurements = np.array([[0.1], [5.2], [10.0]])
>>> H = np.array([[1., 0.]])  # Measure position
>>> R = np.array([[0.1]])
>>> result = jpda_update([x1, x2], [P, P], measurements, H, R)
>>> len(result.states)
2

Notes

The JPDA update consists of: 1. Compute measurement likelihoods and gating 2. Compute association probabilities 3. Compute combined innovations for each track 4. Update each track with weighted innovation 5. Compute covariance with spread of means term

References

  • Bar-Shalom, Y. and Fortmann, T.E., “Tracking and Data Association”, Academic Press, 1988.

pytcl.assignment_algorithms.jpda.jpda(track_states, track_covariances, measurements, H, R, detection_prob=0.9, clutter_density=1e-06, gate_probability=0.99)[source]

Compute JPDA association probabilities.

This is a convenience function that computes association probabilities without performing the state update.

Parameters:
  • track_states (list of array_like) – Predicted state estimates for each track.

  • track_covariances (list of array_like) – Predicted covariances for each track.

  • measurements (array_like) – Measurements, shape (n_meas, m).

  • H (array_like) – Measurement matrix, shape (m, n).

  • R (array_like) – Measurement noise covariance, shape (m, m).

  • detection_prob (float) – Probability of detection. Default 0.9.

  • clutter_density (float) – Spatial density of clutter. Default 1e-6.

  • gate_probability (float) – Probability for chi-squared gate. Default 0.99.

Returns:

result – Association probabilities and related information.

Return type:

JPDAResult

Examples

Compute association probabilities for 2 tracks and 3 measurements:

>>> import numpy as np
>>> # Two tracks with [x, vx] state
>>> states = [np.array([0.0, 1.0]), np.array([10.0, -0.5])]
>>> covariances = [np.eye(2) * 0.5, np.eye(2) * 0.5]
>>> # Three position measurements
>>> measurements = np.array([[0.1], [9.8], [5.0]])
>>> H = np.array([[1, 0]])  # measure position only
>>> R = np.array([[0.1]])
>>> result = jpda(states, covariances, measurements, H, R)
>>> result.association_probs.shape  # (2 tracks, 4 columns: 3 meas + miss)
(2, 4)
>>> # Track 0 should have high prob for measurement 0
>>> result.association_probs[0, 0] > 0.5
True

See also

jpda_update

JPDA with state update.

Nd Assignment

N-dimensional assignment algorithms (4D and higher).

This module extends the 3D assignment solver to arbitrary dimensions, enabling more complex assignment scenarios such as: - 4D: Measurements × Tracks × Hypotheses × Sensors - 5D+: Additional dimensions for time frames, maneuver classes, etc.

The module provides a unified interface for solving high-dimensional assignment problems using generalized relaxation methods.

Performance Notes

For sparse cost tensors (mostly invalid assignments), use SparseCostTensor to reduce memory usage by up to 50% and improve performance on large problems.

References

  • Poore, A. B., “Multidimensional Assignment Problem and Data Association,” IEEE Transactions on Aerospace and Electronic Systems, 2013.

  • Cramer, R. D., et al., “The Emerging Role of Chemical Similarity in Drug Discovery,” Perspectives in Drug Discovery and Design, 2003.

class pytcl.assignment_algorithms.nd_assignment.AssignmentNDResult(assignments, cost, converged, n_iterations, gap)[source]

Bases: NamedTuple

Result of an N-dimensional assignment problem.

Variables:
  • assignments (ndarray) – Array of shape (n_assignments, n_dimensions) containing assigned index tuples. Each row is an n-tuple of indices.

  • cost (float) – Total assignment cost.

  • converged (bool) – Whether the algorithm converged (for iterative methods).

  • n_iterations (int) – Number of iterations used (for iterative methods).

  • gap (float) – Optimality gap (upper_bound - lower_bound) for relaxation methods.

assignments: ndarray[tuple[Any, ...], dtype[int64]]

Alias for field number 0

cost: float

Alias for field number 1

converged: bool

Alias for field number 2

n_iterations: int

Alias for field number 3

gap: float

Alias for field number 4

class pytcl.assignment_algorithms.nd_assignment.SparseCostTensor(dims, indices, costs, default_cost=inf)[source]

Bases: object

Sparse representation of N-dimensional cost tensor.

For assignment problems where most entries represent invalid assignments (infinite cost), storing only valid entries reduces memory by 50% or more and speeds up greedy algorithms.

Variables:
  • dims (tuple) – Shape of the full tensor (n1, n2, …, nk).

  • indices (ndarray) – Array of shape (n_valid, n_dims) with valid entry indices.

  • costs (ndarray) – Array of shape (n_valid,) with costs for valid entries.

  • default_cost (float) – Cost for entries not explicitly stored (default: inf).

Examples

>>> import numpy as np
>>> # Create sparse tensor for 10x10x10 problem with 50 valid entries
>>> dims = (10, 10, 10)
>>> valid_indices = np.random.randint(0, 10, size=(50, 3))
>>> valid_costs = np.random.rand(50)
>>> sparse = SparseCostTensor(dims, valid_indices, valid_costs)
>>> sparse.n_valid
50
>>> sparse.sparsity  # Fraction of valid entries
0.05
>>> # Convert from dense tensor with inf for invalid
>>> dense = np.full((5, 5, 5), np.inf)
>>> dense[0, 0, 0] = 1.0
>>> dense[1, 1, 1] = 2.0
>>> sparse = SparseCostTensor.from_dense(dense)
>>> sparse.n_valid
2
__init__(dims, indices, costs, default_cost=inf)[source]

Initialize sparse cost tensor.

Parameters:
  • dims (tuple) – Shape of the full tensor.

  • indices (ndarray) – Valid entry indices, shape (n_valid, n_dims).

  • costs (ndarray) – Costs for valid entries, shape (n_valid,).

  • default_cost (float) – Cost for invalid (unstored) entries.

property n_dims: int

Number of dimensions.

property n_valid: int

Number of valid (finite cost) entries.

property sparsity: float

Fraction of tensor that is valid (0 to 1).

property memory_savings: float

Estimated memory savings vs dense representation (0 to 1).

get_cost(index)[source]

Get cost for a specific index tuple.

to_dense()[source]

Convert to dense tensor representation.

Returns:

dense – Full tensor with default_cost for unstored entries.

Return type:

ndarray

Notes

May use significant memory for large tensors.

classmethod from_dense(dense, threshold=10000000000.0)[source]

Create sparse tensor from dense array.

Parameters:
  • dense (ndarray) – Dense cost tensor.

  • threshold (float) – Entries above this value are considered invalid. Default 1e10 (catches np.inf and large values).

Returns:

Sparse representation.

Return type:

SparseCostTensor

Examples

>>> import numpy as np
>>> dense = np.array([[[1, np.inf], [np.inf, 2]],
...                   [[np.inf, 3], [4, np.inf]]])
>>> sparse = SparseCostTensor.from_dense(dense)
>>> sparse.n_valid
4
pytcl.assignment_algorithms.nd_assignment.validate_cost_tensor(cost_tensor)[source]

Validate cost tensor and return dimensions.

Parameters:

cost_tensor (ndarray) – Cost tensor of arbitrary dimension.

Returns:

dims – Dimensions of the cost tensor.

Return type:

tuple

Raises:

ValueError – If tensor has fewer than 2 dimensions.

Examples

>>> import numpy as np
>>> cost = np.random.rand(3, 4, 5)
>>> dims = validate_cost_tensor(cost)
>>> dims
(3, 4, 5)
>>> # 1D tensor should raise error
>>> try:
...     validate_cost_tensor(np.array([1, 2, 3]))
... except ValueError:
...     print("Caught expected error")
Caught expected error
pytcl.assignment_algorithms.nd_assignment.greedy_assignment_nd(cost_tensor, max_assignments=None)[source]

Greedy solver for N-dimensional assignment.

Selects minimum-cost tuples in order until no more valid assignments exist (no dimension index is repeated).

Parameters:
  • cost_tensor (ndarray) – Cost tensor of shape (n1, n2, …, nk).

  • max_assignments (int, optional) – Maximum number of assignments to find (default: min(dimensions)).

Returns:

Assignments, total cost, and algorithm info.

Return type:

AssignmentNDResult

Examples

>>> import numpy as np
>>> # 3D cost tensor: 3 measurements x 2 tracks x 2 hypotheses
>>> cost = np.array([
...     [[1.0, 5.0], [3.0, 2.0]],   # meas 0
...     [[4.0, 1.0], [2.0, 6.0]],   # meas 1
...     [[2.0, 3.0], [5.0, 1.0]],   # meas 2
... ])
>>> result = greedy_assignment_nd(cost)
>>> result.cost  # Total cost of greedy solution
2.0
>>> len(result.assignments)  # Number of assignments made
2

Notes

Greedy assignment is fast O(n log n) but not optimal. Used as heuristic or starting solution for optimization methods.

pytcl.assignment_algorithms.nd_assignment.greedy_assignment_nd_sparse(sparse_cost, max_assignments=None)[source]

Greedy solver for sparse N-dimensional assignment.

Selects minimum-cost tuples from valid entries only, which is much faster than dense greedy when sparsity < 0.5.

Parameters:
  • sparse_cost (SparseCostTensor) – Sparse cost tensor with valid entries only.

  • max_assignments (int, optional) – Maximum number of assignments (default: min(dimensions)).

Returns:

Assignments, total cost, and algorithm info.

Return type:

AssignmentNDResult

Examples

>>> import numpy as np
>>> # Create sparse problem
>>> dims = (10, 10, 10)
>>> # Only 20 valid assignments out of 1000
>>> indices = np.array([[i, i, i] for i in range(10)] +
...                    [[i, (i+1)%10, (i+2)%10] for i in range(10)])
>>> costs = np.random.rand(20)
>>> sparse = SparseCostTensor(dims, indices, costs)
>>> result = greedy_assignment_nd_sparse(sparse)
>>> result.converged
True

Notes

Time complexity is O(n_valid * log(n_valid)) vs O(total_size * log(total_size)) for dense greedy. For a 10x10x10 tensor with 50 valid entries, this is 50*log(50) vs 1000*log(1000), about 20x faster.

pytcl.assignment_algorithms.nd_assignment.relaxation_assignment_nd(cost_tensor, max_iterations=100, tolerance=1e-06, verbose=False)[source]

Lagrangian relaxation solver for N-dimensional assignment.

Relaxes the constraints on dimensions 3..N and solves the resulting two-dimensional problem exactly, which yields a valid lower bound on the optimal cost. Subgradient ascent on the multipliers tightens the bound; a feasible solution is recovered at each iteration to provide an upper bound.

Parameters:
  • cost_tensor (ndarray) – Cost tensor of shape (n1, n2, …, nk).

  • max_iterations (int, optional) – Maximum iterations (default 100).

  • tolerance (float, optional) – Convergence tolerance for the optimality gap (default 1e-6).

  • verbose (bool, optional) – Print iteration info (default False).

Returns:

Assignments, total cost, convergence info, and optimality gap. gap is best_upper_bound - best_lower_bound; when converged is True the returned assignment is provably within tolerance of optimal.

Return type:

AssignmentNDResult

Examples

>>> import numpy as np
>>> # 3x3x3 assignment problem
>>> rng = np.random.default_rng(42)
>>> cost = rng.random((3, 3, 3))
>>> result = relaxation_assignment_nd(cost, max_iterations=50)
>>> bool(result.gap >= -1e-9)  # gap is a genuine bound
True
>>> result.assignments.shape[1]  # 3D assignments
3

Notes

The relaxation follows Poore’s formulation (see the module references):

  1. Relax the constraints on dimensions 3..N with multipliers lambda.

  2. The inner problem separates: minimizing over the relaxed dimensions for each (i, j) pair leaves a 2-D assignment problem, solved exactly. Solving it exactly is what makes L(lambda) a valid lower bound – a greedy inner solve does not bound the optimum and can certify optimality for suboptimal answers.

  3. Recover a feasible solution (upper bound) by assigning each relaxed dimension with an exact 2-D assignment.

  4. Ascend the multipliers along the constraint-violation subgradient.

The reported gap uses the best bounds seen across all iterations, so it is a genuine optimality certificate rather than a heuristic score.

pytcl.assignment_algorithms.nd_assignment.auction_assignment_nd(cost_tensor, max_iterations=100, epsilon=0.01, verbose=False)[source]

Auction algorithm for N-dimensional assignment.

Inspired by the classical auction algorithm for 2D assignment, adapted to higher dimensions. Objects bid for assignments based on relative costs.

Parameters:
  • cost_tensor (ndarray) – Cost tensor of shape (n1, n2, …, nk).

  • max_iterations (int, optional) – Maximum iterations (default 100).

  • epsilon (float, optional) – Bid increment (default 0.01). Larger epsilon → fewer iterations, worse solution; smaller epsilon → more iterations, better solution.

  • verbose (bool, optional) – Print iteration info (default False).

Returns:

Assignments, total cost, convergence info, gap estimate.

Return type:

AssignmentNDResult

Examples

>>> import numpy as np
>>> # 4D assignment: sensors x measurements x tracks x hypotheses
>>> np.random.seed(123)
>>> cost = np.random.rand(2, 3, 3, 2) * 10
>>> result = auction_assignment_nd(cost, max_iterations=50, epsilon=0.1)
>>> len(result.assignments) > 0
True
>>> result.n_iterations <= 50
True

Notes

The algorithm maintains a “price” for each index and allows bidding (price adjustment) to maximize value. Converges to epsilon-optimal solution in finite iterations.

pytcl.assignment_algorithms.nd_assignment.detect_dimension_conflicts(assignments, dims)[source]

Check if assignments violate dimension uniqueness.

For valid assignment, each index should appear at most once per dimension.

Parameters:
  • assignments (ndarray) – Array of shape (n_assignments, n_dimensions) with assignments.

  • dims (tuple) – Dimensions of the cost tensor.

Returns:

has_conflicts – True if any index appears more than once in any dimension.

Return type:

bool

Examples

>>> import numpy as np
>>> # Valid assignment: no index repeated in any dimension
>>> assignments = np.array([[0, 0], [1, 1]])
>>> detect_dimension_conflicts(assignments, (3, 3))
False
>>> # Invalid: index 0 used twice in first dimension
>>> assignments = np.array([[0, 0], [0, 1]])
>>> detect_dimension_conflicts(assignments, (3, 3))
True
pytcl.assignment_algorithms.nd_assignment.assignment_nd(cost, method='auto', max_assignments=None, max_iterations=100, tolerance=1e-06, epsilon=0.01, verbose=False)[source]

Unified interface for N-dimensional assignment.

Automatically selects between dense and sparse algorithms based on input type and sparsity.

Parameters:
  • cost (ndarray or SparseCostTensor) – Cost tensor (dense) or sparse cost representation.

  • method (str) – Algorithm to use: ‘auto’, ‘greedy’, ‘relaxation’, ‘auction’. ‘auto’ selects greedy for sparse, relaxation for dense.

  • max_assignments (int, optional) – Maximum number of assignments for greedy methods.

  • max_iterations (int) – Maximum iterations for iterative methods.

  • tolerance (float) – Convergence tolerance for relaxation.

  • epsilon (float) – Price increment for auction algorithm.

  • verbose (bool) – Print progress information.

Returns:

Assignment solution.

Return type:

AssignmentNDResult

Examples

>>> import numpy as np
>>> # Dense usage
>>> cost = np.random.rand(4, 4, 4)
>>> result = assignment_nd(cost, method='greedy')
>>> result.converged
True
>>> # Sparse usage (more efficient for large sparse problems)
>>> dense = np.full((20, 20, 20), np.inf)
>>> for i in range(20):
...     dense[i, i, i] = np.random.rand()
>>> sparse = SparseCostTensor.from_dense(dense)
>>> result = assignment_nd(sparse, method='auto')
>>> result.converged
True

See also

greedy_assignment_nd

Dense greedy algorithm.

greedy_assignment_nd_sparse

Sparse greedy algorithm.

relaxation_assignment_nd

Lagrangian relaxation.

auction_assignment_nd

Auction algorithm.

Network Flow

Network flow solutions for assignment problems.

This module provides min-cost flow formulations for assignment problems, offering an alternative to Hungarian algorithm and relaxation methods.

A min-cost flow approach: 1. Models assignment as flow network 2. Uses cost edges for penalties 3. Enforces supply/demand constraints 4. Finds minimum-cost flow solution 5. Extracts assignment from flow

References

  • Ahuja, R. K., Magnanti, T. L., & Orlin, J. B. (1993). Network Flows: Theory, Algorithms, and Applications. Prentice-Hall.

  • Costain, G., & Liang, H. (2012). An Auction Algorithm for the Minimum Cost Flow Problem. CoRR, abs/1208.4859.

class pytcl.assignment_algorithms.network_flow.FlowStatus(value)[source]

Bases: Enum

Status of min-cost flow computation.

OPTIMAL = 0
UNBOUNDED = 1
INFEASIBLE = 2
TIMEOUT = 3
class pytcl.assignment_algorithms.network_flow.MinCostFlowResult(flow, cost, status, iterations)[source]

Bases: NamedTuple

Result of min-cost flow computation.

Variables:
  • flow (ndarray) – Flow values on each edge, shape (n_edges,).

  • cost (float) – Total flow cost.

  • status (FlowStatus) – Optimization status.

  • iterations (int) – Number of iterations used.

flow: ndarray[tuple[Any, ...], dtype[float64]]

Alias for field number 0

cost: float

Alias for field number 1

status: FlowStatus

Alias for field number 2

iterations: int

Alias for field number 3

class pytcl.assignment_algorithms.network_flow.FlowEdge(from_node, to_node, capacity, cost)[source]

Bases: NamedTuple

Edge in a flow network.

Variables:
  • from_node (int) – Source node index.

  • to_node (int) – Destination node index.

  • capacity (float) – Maximum flow on edge (default 1.0 for assignment).

  • cost (float) – Cost per unit flow.

from_node: int

Alias for field number 0

to_node: int

Alias for field number 1

capacity: float

Alias for field number 2

cost: float

Alias for field number 3

pytcl.assignment_algorithms.network_flow.assignment_to_flow_network(cost_matrix)[source]

Convert 2D assignment problem to min-cost flow network.

Network structure: - Source node (0) supplies all workers - Worker nodes (1 to m) demand 1 unit each - Task nodes (m+1 to m+n) supply 1 unit each - Sink node (m+n+1) collects all completed tasks

Parameters:

cost_matrix (ndarray) – Cost matrix of shape (m, n) where cost[i,j] is cost of assigning worker i to task j.

Returns:

  • edges (list[FlowEdge]) – List of edges in the flow network.

  • supplies (ndarray) – Supply/demand at each node (shape n_nodes,). Positive = supply, negative = demand.

  • node_names (ndarray) – Names of nodes for reference.

Return type:

Tuple[list[FlowEdge], ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[Any]]]

Examples

>>> import numpy as np
>>> # 2 workers, 3 tasks
>>> cost = np.array([[1.0, 2.0, 3.0],
...                  [4.0, 5.0, 6.0]])
>>> edges, supplies, names = assignment_to_flow_network(cost)
>>> len(edges)  # source->workers + workers->tasks + tasks->sink
11
>>> supplies[0]  # source supplies 2 (num workers)
2.0
>>> names[0]
'source'
pytcl.assignment_algorithms.network_flow.min_cost_flow_successive_shortest_paths(edges, supplies, max_iterations=1000)[source]

Solve min-cost flow using successive shortest paths with cost scaling.

Algorithm:

1. Initialize potentials using Bellman-Ford
2. While there is excess supply:
   - Find shortest path using reduced costs (Dijkstra with potentials)
   - Push unit flow along path
   - Update node potentials
   - Recompute shortest paths to maintain optimality

This is the standard min-cost flow algorithm that guarantees optimality and convergence. It uses Dijkstra’s algorithm with potentials, which maintains the dual feasibility (reduced cost property).

Parameters:
  • edges (list[FlowEdge]) – List of edges with capacities and costs.

  • supplies (ndarray) – Supply/demand at each node.

  • max_iterations (int, optional) – Maximum iterations (default 1000).

Returns:

Solution with flow values, cost, status, and iterations.

Return type:

MinCostFlowResult

Examples

>>> import numpy as np
>>> from pytcl.assignment_algorithms.network_flow import (
...     FlowEdge, assignment_to_flow_network
... )
>>> cost = np.array([[1.0, 5.0], [4.0, 2.0]])
>>> edges, supplies, _ = assignment_to_flow_network(cost)
>>> result = min_cost_flow_successive_shortest_paths(edges, supplies)
>>> result.status == FlowStatus.OPTIMAL
True
>>> result.cost  # Optimal is 1+2=3 (diagonal)
3.0

Notes

This implementation uses successive shortest paths with potentials. The algorithm is guaranteed to find the optimal solution for any feasible min-cost flow problem.

For rectangular assignment problems (m < n or m > n), all m units of flow must be satisfied. The algorithm ensures this by finding augmenting paths until all supply is routed.

pytcl.assignment_algorithms.network_flow.min_cost_flow_simplex(edges, supplies, max_iterations=10000)[source]

Solve min-cost flow using Dijkstra-based successive shortest paths.

This optimized version uses:

  • Dijkstra’s algorithm (O(E log V)) per shortest-path iteration, after one initial Bellman-Ford relaxation (O(VE)) to establish valid node potentials

  • Node potentials to maintain non-negative edge costs

  • Johnson’s technique for cost adjustment

This is significantly faster than Bellman-Ford while maintaining guaranteed correctness and optimality.

Time complexity: O(V*E + K*E log V), where the first term is the potential-initialising Bellman-Ford pass and K is the number of shortest paths. Space complexity: O(V + E).

Parameters:
  • edges (list[FlowEdge]) – List of edges with capacities and costs.

  • supplies (ndarray) – Supply/demand at each node.

  • max_iterations (int, optional) – Maximum iterations (default 10000).

Returns:

Solution with flow values, cost, status, and iterations.

Return type:

MinCostFlowResult

References

  • Ahuja, R. K., Magnanti, T. L., & Orlin, J. B. (1993). Network Flows: Theory, Algorithms, and Applications. (Chapter on successive shortest paths with potentials)

  • Johnson, D. B. (1977). Efficient All-Pairs Shortest Paths in Weighted Graphs.

pytcl.assignment_algorithms.network_flow.assignment_from_flow_solution(flow, edges, cost_matrix_shape)[source]

Extract assignment from flow network solution.

A valid flow solution for assignment should have: - Exactly 1 unit of flow from each worker to some task - Exactly 1 unit of flow to each task from some worker - No negative flows on worker->task edges (those are cancellations)

This function extracts the actual assignment by identifying which worker->task edges carry the net positive flow.

Parameters:
  • flow (ndarray) – Flow values on each edge.

  • edges (list[FlowEdge]) – List of edges used in network.

  • cost_matrix_shape (tuple) – Shape of original cost matrix (m, n).

Returns:

  • assignment (ndarray) – Assignment array of shape (n_assignments, 2) with [worker, task].

  • cost (float) – Total assignment cost.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[int64]], float]

pytcl.assignment_algorithms.network_flow.min_cost_assignment_via_flow(cost_matrix, use_simplex=True)[source]

Solve 2D assignment problem via min-cost flow network.

Uses Dijkstra-optimized successive shortest paths (Phase 1B) by default. Falls back to Bellman-Ford if needed.

Parameters:
  • cost_matrix (ndarray) – Cost matrix of shape (m, n).

  • use_simplex (bool, optional) – Use Dijkstra-optimized algorithm (default True) or Bellman-Ford based successive shortest paths (False).

Returns:

  • assignment (ndarray) – Assignment array of shape (n_assignments, 2).

  • total_cost (float) – Total assignment cost.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[int64]], float]

Examples

>>> import numpy as np
>>> cost = np.array([[1.0, 5.0, 9.0],
...                  [3.0, 2.0, 8.0],
...                  [7.0, 6.0, 4.0]])
>>> assignment, total_cost = min_cost_assignment_via_flow(cost)
>>> total_cost  # Optimal assignment: (0,0), (1,1), (2,2) = 1+2+4 = 7
7.0
>>> len(assignment)
3

Notes

Phase 1B: Dijkstra-based optimization provides O(K*E log V) vs Bellman-Ford O(K*V*E), where K is number of shortest paths needed.

Assignment

Three-dimensional assignment algorithms.

This module provides algorithms for solving the 3D assignment (axial 3-index assignment) problem, which arises in multi-sensor data fusion and multi-scan target tracking.

The 3D assignment problem extends the 2D problem to three dimensions: given a cost tensor C[i,j,k], find an assignment that minimizes the total cost subject to the constraint that each index appears in at most one selected tuple.

class pytcl.assignment_algorithms.three_dimensional.assignment.Assignment3DResult(tuples, cost, converged, n_iterations, gap)[source]

Bases: NamedTuple

Result of a 3D assignment problem.

Variables:
  • tuples (ndarray) – Array of shape (n_assignments, 3) containing assigned index tuples. Each row is (i, j, k) representing an assignment.

  • cost (float) – Total assignment cost.

  • converged (bool) – Whether the algorithm converged (for iterative methods).

  • n_iterations (int) – Number of iterations used (for iterative methods).

  • gap (float) – Optimality gap (upper_bound - lower_bound) for relaxation methods.

tuples: ndarray[tuple[Any, ...], dtype[int64]]

Alias for field number 0

cost: float

Alias for field number 1

converged: bool

Alias for field number 2

n_iterations: int

Alias for field number 3

gap: float

Alias for field number 4

pytcl.assignment_algorithms.three_dimensional.assignment.greedy_3d(cost_tensor, maximize=False)[source]

Solve 3D assignment using a greedy algorithm.

This algorithm iteratively selects the lowest-cost unassigned tuple until no more valid assignments can be made. It provides a fast but suboptimal solution.

Parameters:
  • cost_tensor (array_like) – Cost tensor of shape (n1, n2, n3).

  • maximize (bool, optional) – If True, solve maximization problem (default: False).

Returns:

result – Assignment result with tuples, cost, and metadata.

Return type:

Assignment3DResult

Examples

>>> import numpy as np
>>> cost = np.random.rand(3, 3, 3)
>>> result = greedy_3d(cost)
>>> result.tuples.shape
(3, 3)

Notes

Time complexity is O(n1 * n2 * n3 * min(n1, n2, n3)) in the worst case. The greedy solution provides a starting point for more sophisticated algorithms but is generally not optimal.

pytcl.assignment_algorithms.three_dimensional.assignment.decompose_to_2d(cost_tensor, fixed_dimension=0, maximize=False)[source]

Solve 3D assignment by decomposing into sequential 2D problems.

This heuristic fixes one dimension and solves a sequence of 2D assignment problems. While not optimal, it provides a polynomial-time approximation.

Parameters:
  • cost_tensor (array_like) – Cost tensor of shape (n1, n2, n3).

  • fixed_dimension (int, optional) – Which dimension to iterate over (0, 1, or 2). Default: 0.

  • maximize (bool, optional) – If True, solve maximization problem (default: False).

Returns:

result – Assignment result.

Return type:

Assignment3DResult

Examples

>>> import numpy as np
>>> cost = np.random.rand(4, 4, 4)
>>> result = decompose_to_2d(cost, fixed_dimension=0)
>>> result.tuples.shape[0] <= 4
True

Notes

The algorithm iterates over the fixed dimension and solves a 2D assignment for each slice. Indices that have been assigned in previous slices are excluded from subsequent problems.

Time complexity is O(n * m^3) where n is the size of the fixed dimension and m is the maximum of the other two dimensions.

pytcl.assignment_algorithms.three_dimensional.assignment.assign3d_lagrangian(cost_tensor, max_iter=100, tol=1e-06, step_size=1.0, maximize=False)[source]

Solve 3D assignment using Lagrangian relaxation.

Relaxes the constraint on the third dimension and solves the resulting two-dimensional problem exactly, yielding a valid lower bound on the optimal cost. Subgradient ascent tightens the bound while feasible solutions supply the upper bound.

Parameters:
  • cost_tensor (array_like) – Cost tensor of shape (n1, n2, n3).

  • max_iter (int, optional) – Maximum number of iterations (default: 100).

  • tol (float, optional) – Convergence tolerance for gap (default: 1e-6).

  • step_size (float, optional) – Initial step size for the subgradient update (default: 1.0). Used as a fallback when the Polyak step is not usable.

  • maximize (bool, optional) – If True, solve maximization problem (default: False).

Returns:

result – Assignment result with optimality gap information. gap is best_upper_bound - best_lower_bound; when converged is True the returned assignment is provably within tol of optimal.

Return type:

Assignment3DResult

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> cost = rng.random((5, 5, 5))
>>> result = assign3d_lagrangian(cost, max_iter=50)
>>> result.tuples.shape[1]
3
>>> bool(result.gap >= -1e-9)  # gap is a genuine bound
True

Notes

The Lagrangian dualizes the third-dimension constraint. For multipliers mu, the inner problem separates: minimizing cost[i, j, k] - mu[k] over k for each (i, j) leaves a 2-D assignment problem, which is solved exactly. Solving that inner problem exactly is what makes L(mu) a valid lower bound – an approximate inner solve does not bound the optimum and can certify optimality for suboptimal answers.

The reported gap uses the best bounds seen across all iterations, so it is a genuine optimality certificate: gap == 0 means the returned assignment is optimal.

References

  • Poore, A.B., “Multidimensional assignment formulation of data association problems arising from multitarget and multisensor tracking”, Computational Optimization and Applications, 1994.

pytcl.assignment_algorithms.three_dimensional.assignment.assign3d_auction(cost_tensor, epsilon=None, max_iter=1000, maximize=False)[source]

Solve 3D assignment using an auction-based algorithm.

This extends the 2D auction algorithm to 3D by treating the problem as a tripartite matching with iterative bidding.

Parameters:
  • cost_tensor (array_like) – Cost tensor of shape (n1, n2, n3).

  • epsilon (float, optional) – Price increment. If None, uses adaptive epsilon.

  • max_iter (int, optional) – Maximum iterations (default: 1000).

  • maximize (bool, optional) – If True, solve maximization problem (default: False).

Returns:

result – Assignment result.

Return type:

Assignment3DResult

Examples

>>> import numpy as np
>>> cost = np.random.rand(4, 4, 4)
>>> result = assign3d_auction(cost)
>>> len(result.tuples) <= 4
True

Notes

The auction algorithm for 3D assignment alternates between phases: 1. Bidding phase: unassigned elements bid for their preferred matches 2. Assignment phase: matches are updated based on bids

This is a heuristic extension of the 2D auction algorithm and may not find the global optimum.

References

  • Bertsekas, D.P., “The auction algorithm for assignment and other network flow problems”, Interfaces, 1990.

pytcl.assignment_algorithms.three_dimensional.assignment.assign3d(cost_tensor, method='lagrangian', maximize=False, **kwargs)[source]

Solve 3D assignment problem.

This is the main entry point for 3D assignment, providing access to multiple algorithms through a unified interface.

Parameters:
  • cost_tensor (array_like) – Cost tensor of shape (n1, n2, n3).

  • method (str, optional) – Algorithm to use: - “lagrangian”: Lagrangian relaxation (default) - “auction”: Auction-based algorithm - “greedy”: Fast greedy heuristic - “decompose”: Sequential 2D decomposition

  • maximize (bool, optional) – If True, solve maximization problem (default: False).

  • **kwargs (Any) – Additional arguments passed to the specific method.

Returns:

result – Assignment result.

Return type:

Assignment3DResult

Examples

>>> import numpy as np
>>> np.random.seed(42)
>>> cost = np.random.rand(5, 5, 5)
>>> result = assign3d(cost, method="lagrangian")
>>> result.tuples.shape
(5, 3)

See also

assign3d_lagrangian

Lagrangian relaxation method

assign3d_auction

Auction-based method

greedy_3d

Greedy heuristic

decompose_to_2d

Sequential decomposition

Notes

The 3D assignment problem is NP-hard, so all methods provide approximate solutions. The Lagrangian relaxation method generally provides the best quality with reasonable computation time.

Assignment

Two-dimensional assignment algorithms.

This module provides algorithms for solving the 2D assignment (bipartite matching) problem, which is fundamental to data association in target tracking.

class pytcl.assignment_algorithms.two_dimensional.assignment.AssignmentResult(row_indices, col_indices, cost, unassigned_rows, unassigned_cols)[source]

Bases: NamedTuple

Result of a 2D assignment problem.

Variables:
  • row_indices (ndarray) – Indices of assigned rows.

  • col_indices (ndarray) – Indices of assigned columns (col_indices[i] is assigned to row_indices[i]).

  • cost (float) – Total assignment cost.

  • unassigned_rows (ndarray) – Indices of unassigned rows.

  • unassigned_cols (ndarray) – Indices of unassigned columns.

row_indices: ndarray[tuple[Any, ...], dtype[int64]]

Alias for field number 0

col_indices: ndarray[tuple[Any, ...], dtype[int64]]

Alias for field number 1

cost: float

Alias for field number 2

unassigned_rows: ndarray[tuple[Any, ...], dtype[int64]]

Alias for field number 3

unassigned_cols: ndarray[tuple[Any, ...], dtype[int64]]

Alias for field number 4

pytcl.assignment_algorithms.two_dimensional.assignment.linear_sum_assignment(cost_matrix, maximize=False)[source]

Solve the linear sum assignment problem (wrapper around scipy).

Parameters:
  • cost_matrix (array_like) – Cost matrix of shape (n, m).

  • maximize (bool, optional) – If True, solve maximization problem instead of minimization.

Returns:

  • row_ind (ndarray) – Row indices of optimal assignment.

  • col_ind (ndarray) – Column indices of optimal assignment.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[int64]], ndarray[tuple[Any, …], dtype[int64]]]

Examples

>>> cost = np.array([[4, 1, 3], [2, 0, 5], [3, 2, 2]])
>>> row_ind, col_ind = linear_sum_assignment(cost)
>>> row_ind
array([0, 1, 2])
>>> col_ind
array([1, 0, 2])
>>> cost[row_ind, col_ind].sum()
5

Notes

This is a thin wrapper around scipy.optimize.linear_sum_assignment for API consistency.

pytcl.assignment_algorithms.two_dimensional.assignment.hungarian(cost_matrix, maximize=False)[source]

Solve 2D assignment using Hungarian (Kuhn-Munkres) algorithm.

The Hungarian algorithm finds the optimal assignment that minimizes (or maximizes) the total cost.

Parameters:
  • cost_matrix (array_like) – Cost matrix of shape (n, m). Can be rectangular.

  • maximize (bool, optional) – If True, solve maximization problem (default: False).

Returns:

  • row_ind (ndarray) – Row indices of optimal assignment.

  • col_ind (ndarray) – Column indices of optimal assignment.

  • total_cost (float) – Total cost of the assignment.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[int64]], ndarray[tuple[Any, …], dtype[int64]], float]

Examples

>>> cost = np.array([[4, 1, 3], [2, 0, 5], [3, 2, 2]])
>>> row_ind, col_ind, total_cost = hungarian(cost)
>>> total_cost
5.0

Notes

Time complexity is O(n^3) for an n x n matrix.

References

  • Kuhn, H.W., “The Hungarian Method for the assignment problem”, Naval Research Logistics Quarterly, 1955.

pytcl.assignment_algorithms.two_dimensional.assignment.auction(cost_matrix, epsilon=None, max_iter=1000, maximize=False)[source]

Solve 2D assignment using the Auction algorithm.

The auction algorithm is an iterative method that can be faster than Hungarian for sparse problems and is naturally parallelizable.

Parameters:
  • cost_matrix (array_like) – Cost matrix of shape (n, m).

  • epsilon (float, optional) – Price increment. If None, uses 1/(n+1) where n is matrix size.

  • max_iter (int, optional) – Maximum iterations (default: 1000).

  • maximize (bool, optional) – If True, solve maximization problem (default: False).

Returns:

  • row_ind (ndarray) – Row indices of the assignment.

  • col_ind (ndarray) – Column indices of the assignment.

  • total_cost (float) – Total cost of the assignment.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[int64]], ndarray[tuple[Any, …], dtype[int64]], float]

Examples

>>> cost = np.array([[4, 1, 3], [2, 0, 5], [3, 2, 2]])
>>> row_ind, col_ind, total_cost = auction(cost)

Notes

The auction algorithm treats rows as “bidders” and columns as “objects”. Each iteration, unassigned bidders bid for their most desirable objects, and objects are assigned to the highest bidder.

The result is epsilon-optimal, not optimal. This docstring used to call it an optimal assignment (gh-20). The auction algorithm terminates when no bidder can improve by more than epsilon, so the returned cost may exceed the true minimum by up to n * epsilon for an n-row problem.

The bound is achieved exactly: for integer costs with epsilon < 1/n the gap is smaller than one unit, so the assignment is genuinely optimal. Both the bound and that integer-cost guarantee were verified during the v2 audit. For real-valued costs, scale epsilon down or use hungarian when exactness matters more than speed.

See also

hungarian

Exact O(n^3) assignment.

References

  • Bertsekas, D.P., “The auction algorithm: A distributed relaxation method for the assignment problem”, Annals of Operations Research, 1988.

pytcl.assignment_algorithms.two_dimensional.assignment.assign2d(cost_matrix, cost_of_non_assignment=inf, maximize=False)[source]

Solve 2D assignment with cost of non-assignment.

This extends the basic assignment problem to allow tracks and measurements to remain unassigned at a specified cost.

Parameters:
  • cost_matrix (array_like) – Cost matrix of shape (n_tracks, n_measurements).

  • cost_of_non_assignment (float, optional) – Cost for leaving a track or measurement unassigned. If inf, all must be assigned (default: inf).

  • maximize (bool, optional) – If True, solve maximization problem (default: False).

Returns:

result – Named tuple containing: - row_indices: Indices of assigned rows (tracks) - col_indices: Indices of assigned columns (measurements) - cost: Total assignment cost - unassigned_rows: Indices of unassigned rows - unassigned_cols: Indices of unassigned columns

Return type:

AssignmentResult

Examples

>>> cost = np.array([[1, 10], [10, 1], [5, 5]])
>>> result = assign2d(cost, cost_of_non_assignment=3)
>>> result.row_indices
array([0, 1])
>>> result.col_indices
array([0, 1])
>>> result.unassigned_rows
array([2])

Notes

The algorithm augments the cost matrix with dummy rows and columns representing non-assignment options, then solves the augmented problem.

Kbest

K-best 2D assignment algorithms.

This module provides algorithms for finding the k best solutions to the 2D assignment problem, which is essential for hypothesis generation in Multiple Hypothesis Tracking (MHT).

class pytcl.assignment_algorithms.two_dimensional.kbest.KBestResult(assignments, costs, n_found)[source]

Bases: NamedTuple

Result of k-best assignment problem.

Variables:
  • assignments (List[AssignmentResult]) – List of k best assignments, sorted by cost (ascending for minimization).

  • costs (ndarray) – Array of costs for each assignment.

  • n_found (int) – Number of assignments found (may be less than k if fewer exist).

assignments: List[AssignmentResult]

Alias for field number 0

costs: ndarray[tuple[Any, ...], dtype[float64]]

Alias for field number 1

n_found: int

Alias for field number 2

pytcl.assignment_algorithms.two_dimensional.kbest.murty(cost_matrix, k, cost_of_non_assignment=inf, maximize=False)[source]

Find k-best assignments using Murty’s algorithm.

Murty’s algorithm systematically partitions the solution space to enumerate assignments in order of increasing cost. It is widely used in Multiple Hypothesis Tracking (MHT) for hypothesis generation.

Parameters:
  • cost_matrix (array_like) – Cost matrix of shape (n_tracks, n_measurements).

  • k (int) – Number of best assignments to find.

  • cost_of_non_assignment (float, optional) – Cost for leaving a track or measurement unassigned. If inf, all must be assigned (default: inf).

  • maximize (bool, optional) – If True, find k-best in descending order (default: False).

Returns:

result – Named tuple containing: - assignments: List of k best AssignmentResult objects - costs: Array of costs for each assignment - n_found: Number of assignments actually found

Return type:

KBestResult

Examples

>>> import numpy as np
>>> cost = np.array([[10, 5, 13], [3, 15, 8], [12, 7, 9]])
>>> result = murty(cost, k=3)
>>> result.n_found
3
>>> result.costs  # Three lowest-cost assignments
array([17., 23., 25.])

Notes

The algorithm has time complexity O(k * n^3) for an n x n matrix, as each partition requires solving a constrained assignment problem.

For large k, consider using lazy evaluation techniques or early termination based on cost thresholds.

References

  • Murty, K.G., “An algorithm for ranking all the assignments in order of increasing cost”, Operations Research, 1968.

  • Miller, M.L., Stone, H.S., and Cox, I.J., “Optimizing Murty’s ranked assignment method”, IEEE Trans. Aerospace and Electronic Systems, 1997.

pytcl.assignment_algorithms.two_dimensional.kbest.kbest_assign2d(cost_matrix, k, cost_of_non_assignment=inf, maximize=False, cost_threshold=None)[source]

Find k-best assignments with optional cost threshold.

This is an enhanced version of Murty’s algorithm that supports early termination when costs exceed a threshold.

Parameters:
  • cost_matrix (array_like) – Cost matrix of shape (n_tracks, n_measurements).

  • k (int) – Maximum number of assignments to find.

  • cost_of_non_assignment (float, optional) – Cost for leaving a track or measurement unassigned.

  • maximize (bool, optional) – If True, find k-best in descending order.

  • cost_threshold (float, optional) – Stop when assignment cost exceeds this threshold (for minimization) or falls below (for maximization). If None, no threshold.

Returns:

result – Named tuple containing assignments, costs, and count.

Return type:

KBestResult

Examples

>>> cost = np.array([[10, 5, 13], [3, 15, 8], [12, 7, 9]])
>>> result = kbest_assign2d(cost, k=10, cost_threshold=20)
>>> result.n_found  # Only assignments with cost <= 20
1

Notes

The cost threshold enables efficient pruning for MHT where hypotheses with very low probability (high cost) can be discarded.

pytcl.assignment_algorithms.two_dimensional.kbest.ranked_assignments(cost_matrix, max_assignments=100, cost_threshold=None, maximize=False)[source]

Enumerate assignments in ranked order (best to worst).

This is a convenience function for MHT hypothesis generation that provides sensible defaults.

Parameters:
  • cost_matrix (array_like) – Cost matrix of shape (n_tracks, n_measurements).

  • max_assignments (int, optional) – Maximum number of assignments to enumerate (default: 100).

  • cost_threshold (float, optional) – Stop when cost exceeds threshold (minimization) or falls below (maximization).

  • maximize (bool, optional) – If True, rank in descending order of cost.

Returns:

result – Ranked assignments.

Return type:

KBestResult

Examples

>>> cost = np.array([[10, 5], [3, 15]])
>>> result = ranked_assignments(cost, max_assignments=5)
>>> len(result.assignments)
2
>>> result.costs
array([ 8., 25.])

Notes

This function assumes all tracks must be assigned to measurements (no non-assignment option). For track-oriented MHT with missed detections, use kbest_assign2d with an appropriate cost_of_non_assignment.