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:
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:
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_thresholdCompute 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:
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:
- Returns:
Chi-squared threshold value.
- Return type:
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:
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_gateTest if measurement passes gate.
mahalanobis_distanceCompute 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:
NamedTupleResult 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.
- 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:
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:
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:
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:
NamedTupleResult 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.
- class pytcl.assignment_algorithms.jpda.JPDAUpdate(states, covariances, association_probs, innovations)[source]
Bases:
NamedTupleResult 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.
- pytcl.assignment_algorithms.jpda.compute_measurement_likelihood(innovation, innovation_cov, detection_prob=1.0)[source]
Compute measurement likelihood for a track-measurement pair.
- 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:
- 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:
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:
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_updateJPDA 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:
NamedTupleResult 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.
- class pytcl.assignment_algorithms.nd_assignment.SparseCostTensor(dims, indices, costs, default_cost=inf)[source]
Bases:
objectSparse 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:
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
- 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:
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:
- 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:
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:
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:
- Returns:
Assignments, total cost, convergence info, and optimality gap.
gapisbest_upper_bound - best_lower_bound; whenconvergedis True the returned assignment is provably withintoleranceof optimal.- Return type:
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):
Relax the constraints on dimensions 3..N with multipliers
lambda.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.Recover a feasible solution (upper bound) by assigning each relaxed dimension with an exact 2-D assignment.
Ascend the multipliers along the constraint-violation subgradient.
The reported
gapuses 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:
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:
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:
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_ndDense greedy algorithm.
greedy_assignment_nd_sparseSparse greedy algorithm.
relaxation_assignment_ndLagrangian relaxation.
auction_assignment_ndAuction 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:
EnumStatus 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:
NamedTupleResult 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.
- status: FlowStatus
Alias for field number 2
- class pytcl.assignment_algorithms.network_flow.FlowEdge(from_node, to_node, capacity, cost)[source]
Bases:
NamedTupleEdge in a flow network.
- Variables:
- 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:
- Returns:
Solution with flow values, cost, status, and iterations.
- Return type:
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:
- Returns:
Solution with flow values, cost, status, and iterations.
- Return type:
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:
- Returns:
assignment (ndarray) – Assignment array of shape (n_assignments, 2) with [worker, task].
cost (float) – Total assignment cost.
- Return type:
- 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:
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:
NamedTupleResult 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.
- 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:
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:
- Returns:
result – Assignment result.
- Return type:
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.
gapisbest_upper_bound - best_lower_bound; whenconvergedis True the returned assignment is provably withintolof optimal.- Return type:
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: minimizingcost[i, j, k] - mu[k]overkfor each(i, j)leaves a 2-D assignment problem, which is solved exactly. Solving that inner problem exactly is what makesL(mu)a valid lower bound – an approximate inner solve does not bound the optimum and can certify optimality for suboptimal answers.The reported
gapuses the best bounds seen across all iterations, so it is a genuine optimality certificate:gap == 0means 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:
- Returns:
result – Assignment result.
- Return type:
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:
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_lagrangianLagrangian relaxation method
assign3d_auctionAuction-based method
greedy_3dGreedy heuristic
decompose_to_2dSequential 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:
NamedTupleResult 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.
- 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:
- 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 ton * epsilonfor ann-row problem.The bound is achieved exactly: for integer costs with
epsilon < 1/nthe 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 usehungarianwhen exactness matters more than speed.See also
hungarianExact 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:
- 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:
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:
NamedTupleResult 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
- 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:
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:
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:
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.