Repository: dstl/Stone-Soup

Issue:
Low efficiency in GOSPA metric computations
I'm using [the GOSPA metric](https://stonesoup.readthedocs.io/en/latest/stonesoup.metricgenerator.ospametric.html#stonesoup.metricgenerator.ospametric.GOSPAMetric) to evaluate the tracking results in a scenario that has generated 40 tracks over 30,000 time steps, and I've been waiting for results for more than a couple of hours. I understand that this isn't a trivial scenario, but it made me think about how efficient is the metric's implementation. 

The description indicates that the metric relies on the auction algorithm, which may not be the fastest option here. Angel Garcia-Fernandez [provides](https://github.com/Agarciafernandez/MTT/blob/master/GOSPA%20code/GOSPA.m) an improved (MATLAB) implementation that uses the [`assign2D`](https://github.com/schwehr/Tracker-Component-Library/blob/master/Assignment%20Algorithms/2D%20Assignment/assign2D.m) function instead. It would be great to have it in Stone Soup, provided that `assign2D` has already been [implemented](https://stonesoup.readthedocs.io/en/latest/stonesoup.dataassociator.html#stonesoup.dataassociator.neighbour.GNNWith2DAssignment) to speed up the GNN algorithm. 

Discussion:

Patch:
diff --git a/stonesoup/metricgenerator/ospametric.py b/stonesoup/metricgenerator/ospametric.py
index 26560b9a..c27195af 100644
--- a/stonesoup/metricgenerator/ospametric.py
+++ b/stonesoup/metricgenerator/ospametric.py
@@ -56,7 +56,7 @@ class GOSPAMetric(MetricGenerator):
     """
     Computes the Generalized Optimal SubPattern Assignment (GOSPA) metric
     for two sets of :class:`~.Track` objects. This implementation of GOSPA
-    is based on the auction algorithm.
+    is based on the modified Jonker-Volgenant algorithm.
 
     The GOSPA metric is calculated at each time step in which a
     :class:`~.Track` object is present
@@ -205,14 +205,13 @@ class GOSPAMetric(MetricGenerator):
                 time_range=TimeRange(min(timestamps), max(timestamps)),
                 generator=self)
 
-    def compute_assignments(self, cost_matrix, max_iter):
-        """Compute assignments using Auction Algorithm.
+    def compute_assignments(self, cost_matrix):
+        """Compute assignments using modified Jonker-Volgenant algorithm
 
         Parameters
         ----------
         cost_matrix: Matrix (size mxn) that denotes the cost of assigning
                             mth truth state to each of the n measured states.
-        max_iter: Maximum number of iterations to perform
 
         Returns
         -------
@@ -228,78 +227,13 @@ class GOSPAMetric(MetricGenerator):
         # Index for objects that will be left un-assigned.
         unassigned_idx = -1
 
-        opt_cost = 0.0
         measured_to_truth = np.full((n_measured, ), unassigned_idx)
         truth_to_measured = np.full((m_truth, ), unassigned_idx)
 
-        if m_truth == 1:
-            # Corner case 1: if there is only one truth state.
-            max_cost_idx = np.argmax(cost_matrix, axis=1).item()
-            opt_cost = cost_matrix[0, max_cost_idx]
-            truth_to_measured[0] = max_cost_idx
-            measured_to_truth[max_cost_idx] = 0
-
-            return truth_to_measured, measured_to_truth, opt_cost
-
-        if n_measured == 1:
-            # Corner case 1: if there is only one measured state.
-            max_cost_idx = np.argmax(cost_matrix, axis=0).item()
-            opt_cost = cost_matrix[max_cost_idx, 0]
-            measured_to_truth[0] = max_cost_idx
-            truth_to_measured[max_cost_idx] = 0
-
-            return truth_to_measured, measured_to_truth, opt_cost
-
-        swap_dim_flag = False
-        epsil = 1. / np.max([m_truth, n_measured])
-
-        if n_measured < m_truth:
-            # The implementation only works when
-            # m_truth <= n_measured
-            # So swap cost matrix
-            cost_matrix = cost_matrix.transpose()
-            m_truth, n_measured = cost_matrix.shape
-            measured_to_truth, truth_to_measured = truth_to_measured, measured_to_truth
-            swap_dim_flag = True
-
-        # Initial cost for each measured state
-        c_measured = np.zeros((n_measured, ))
-        k_iter = 0
-
-        while not np.all(truth_to_measured != unassigned_idx) and k_iter <= max_iter:
-            for i in range(m_truth):
-                if truth_to_measured[i] == unassigned_idx:
-                    # Unassigned truth object 'i' bids for the best
-                    # measured object j_star
-
-                    # Value for each measured object for truth 'i'
-                    tmp_mat = cost_matrix[i, :] - c_measured
-                    j = np.argsort(tmp_mat)[::-1]
-                    # Best measurement for truth 'i'
-                    j_star = j[0]
-                    # 1st and 2nd best value for truth 'i'
-                    v_i_j_star, w_i_j_star = tmp_mat[j[:2]]
-
-                    # Bid for measured j_star
-                    if w_i_j_star != -np.inf:
-                        c_measured[j_star] += v_i_j_star - w_i_j_star + epsil
-                    else:
-                        c_measured[j_star] += v_i_j_star + epsil
-
-                    # If j_star is unassigned
-                    if measured_to_truth[j_star] != unassigned_idx:
-                        opt_cost -= cost_matrix[measured_to_truth[j_star], j_star]
-                        truth_to_measured[measured_to_truth[j_star]] = unassigned_idx
-
-                    measured_to_truth[j_star] = i
-                    truth_to_measured[i] = j_star
-
-                    # update the cost of new assignment
-                    opt_cost += cost_matrix[i, j_star]
-            k_iter += 1
-
-        if swap_dim_flag:
-            measured_to_truth, truth_to_measured = truth_to_measured, measured_to_truth
+        row_ind, col_ind = linear_sum_assignment(-cost_matrix)
+        opt_cost = cost_matrix[row_ind, col_ind].sum()
+        truth_to_measured[row_ind] = col_ind
+        measured_to_truth[col_ind] = row_ind
 
         return truth_to_measured, measured_to_truth, opt_cost
 
@@ -391,12 +325,11 @@ class GOSPAMetric(MetricGenerator):
             if self.alpha == 2:
                 gospa_metric['missed'] = opt_cost
         else:
-            # Use auction algorithm when both truth_states
+            # Use assignment algorithm when both truth_states
             # and measured_states are non-empty
             cost_matrix = -1. * np.power(cost_matrix, self.p)
             truth_to_measured_assignment, measured_to_truth_assignment, _ =\
-                self.compute_assignments(cost_matrix,
-                                         10 * num_truth_states * num_measured_states)
+                self.compute_assignments(cost_matrix)
 
             opt_cost -= np.sum(measured_to_truth_assignment == unassigned_index) * dummy_cost
             if self.alpha == 2:


Test patch:
diff --git a/stonesoup/metricgenerator/tests/test_ospametric.py b/stonesoup/metricgenerator/tests/test_ospametric.py
index 0dadce50..814b85c7 100644
--- a/stonesoup/metricgenerator/tests/test_ospametric.py
+++ b/stonesoup/metricgenerator/tests/test_ospametric.py
@@ -51,8 +51,7 @@ def test_gospametric_compute_assignments(num_states):
                                                 truth_obj.states)
     neg_cost_matrix = -1.*cost_matrix
     meas_to_truth, truth_to_meas, opt_cost =\
-        generator.compute_assignments(neg_cost_matrix,
-                                      10 * num_states * num_states)
+        generator.compute_assignments(neg_cost_matrix)
 
     assert opt_cost == 0.0
     assert np.array_equal(meas_to_truth, truth_to_meas)
@@ -64,8 +63,7 @@ def test_gospametric_compute_assignments(num_states):
                                                 truth_obj.states)
     neg_cost_matrix = -1.*cost_matrix
     meas_to_truth, truth_to_meas, opt_cost = \
-        generator.compute_assignments(neg_cost_matrix,
-                                      10 * num_states * num_states)
+        generator.compute_assignments(neg_cost_matrix)
 
     assert opt_cost == 0.0
     assert np.array_equal(meas_to_truth, truth_to_meas[:-1])
@@ -78,8 +76,7 @@ def test_gospametric_compute_assignments(num_states):
                                                 truth_obj.states[:-1])
     neg_cost_matrix = -1.*cost_matrix
     meas_to_truth, truth_to_meas, opt_cost = \
-        generator.compute_assignments(neg_cost_matrix,
-                                      10 * num_states * num_states)
+        generator.compute_assignments(neg_cost_matrix)
 
     assert opt_cost == 0.0
     assert np.array_equal(meas_to_truth[:-1], truth_to_meas)