src-d / lapjv

Linear Assignmment Problem solver using Jonker-Volgenant algorithm - Python 3 native module.
MIT License
255 stars 29 forks source link

LAPJV is not giving the correct result. #83

Open mayejudi opened 1 year ago

mayejudi commented 1 year ago

I followed the instructions and installed the package using pip3 install lapjv. However, when I test it with a simple matrix, the results are incorrect. I compared the results with scipy linear_sum_assignment.

from lapjv import lapjv
from scipy.optimize import linear_sum_assignment as lsa
import numpy as np
# cost_matrix = np.array([[9.0, 7.6, 7.5,100.0],[3.5, 8.5, 5.5,100.0],[12.5, 9.5, 9.0,100.0],[4.5, 11.0, 9.5,100.0]])
cost_matrix = np.array([[9.0, 7.6, 7.5,7.0],[3.5, 8.5, 5.5,6.5],[12.5, 9.5, 9.0,10.5],[4.5, 11.0, 9.5,11.5]])
row_ind, col_ind, _ = lapjv(cost_matrix)
cost_tot = 0
print("===LAPJV sol: ")
for id in range (len(row_ind)):
    i = row_ind[id]
    j = col_ind[id]
    print (f"track {i} match det {j} cost {cost_matrix[i,j]}")
    cost_tot += cost_matrix[i,j]
print (f"cost total = {cost_tot}")

row_ind1, col_ind1= lsa(cost_matrix)
print("===LSA sol: ")
cost_tot = 0
for id in range (len(row_ind)):
    i = row_ind1[id]
    j = col_ind1[id]
    print (f"track {i} match det {j} cost {cost_matrix[i,j]}")

    cost_tot += cost_matrix[i,j]
print (f"cost total = {cost_tot}")

I get

===LAPJV sol: 
track 3 match det 3 cost 11.5
track 2 match det 2 cost 9.0
track 1 match det 1 cost 8.5
track 0 match det 0 cost 9.0
cost total = 38.0
===LSA sol: 
track 0 match det 3 cost 7.0
track 1 match det 2 cost 5.5
track 2 match det 1 cost 9.5
track 3 match det 0 cost 4.5
cost total = 26.5

I also took the cpp code in the repo created a main. The project is here lapjv-cpp.zip and run the same squared matrix and I get this

9.000000 7.600000 7.500000 7.000000 
3.500000 8.500000 5.500000 6.500000 
12.500000 9.500000 9.000000 10.500000 
4.500000 11.000000 9.500000 11.500000 
start
Beginning lapjv method.
AVX2: enabled
lapjv: COLUMN REDUCTION finished
lapjv: REDUCTION TRANSFER finished
lapjv: AUGMENTING ROW REDUCTION 1 / 2
lapjv: AUGMENTING ROW REDUCTION 2 / 2
lapjv: AUGMENT SOLUTION finished
lapjv: optimal cost calculated
Track 2 assigned to detect 3.  Cost: 10.5
Track 1 assigned to detect 1.  Cost: 8.5
Track 3 assigned to detect 0.  Cost: 4.5
Track 0 assigned to detect 2.  Cost: 7.5
Total cost = 31

dim     4 - lap cost  3.500 - runtime  0.017 ms

Something is up with this code. I would appreciate any help with it. Am I missing something?

vmarkovtsev commented 1 year ago

Hi @mayejudi It's been many years since I last touched this so I don't remember much.

mayejudi commented 1 year ago

Thank you for your prompt answer. I think the problem is in the initial code released by Jonker. This repo was very helpful. I ran the original Jonker code and produced the same results as this repo cpp results. The difference is that Jonker original is integer squared cost matrices, while this repo is for float/ double square matrices. And I might have to read his paper because this might be very important. There are other algorithms like Linear Sum assignment that theoretically they guarantee results for matrices of integers but not floats, because with floats/doubles you can always make it an epsilon smaller (if you are minimizing). If I make it work, I post the code back and you might merge it. Also, the assignments (which is what I am interested) use :

for (int i = 0; i < dim; i++){
        j = rowsol[i];

instead of

for (int k = 0; k < dim; k++){
        i = rowsol[k];
        j= colsol[k];

That is why the assignment comes X to X

mayejudi commented 1 year ago

I did the same modifications suggested in this repo to the C code, I did not touch the functions related to AVX2. Also I did not touched the python code. The results did not change value. So, it might be that the compiler is saying is doing AVX2 and things do not go well with the results on those calls. The modification were more related in the order of things are calculated. I leave my modifications here in case someone else want to pay with this.lapjv-corrected.zip

CSautier commented 11 months ago

For future reference, the results of scipy's LSA and of this repo's output have different formatting (see https://github.com/gatagat/lap/issues/10 for reference). You can correct your code as follows, and you'll see it gives correct results with this library.

from lapjv import lapjv
from scipy.optimize import linear_sum_assignment as lsa
import numpy as np
# cost_matrix = np.array([[9.0, 7.6, 7.5,100.0],[3.5, 8.5, 5.5,100.0],[12.5, 9.5, 9.0,100.0],[4.5, 11.0, 9.5,100.0]])
cost_matrix = np.array([[9.0, 7.6, 7.5,7.0],[3.5, 8.5, 5.5,6.5],[12.5, 9.5, 9.0,10.5],[4.5, 11.0, 9.5,11.5]])
row_ind, col_ind, _ = lapjv(cost_matrix)
cost_tot = 0
print("===LAPJV sol: ")
for id in range (len(row_ind)):
    i = id
    j = row_ind[id]
    print (f"track {i} match det {j} cost {cost_matrix[i,j]}")
    cost_tot += cost_matrix[i,j]
print (f"cost total = {cost_tot}")

row_ind1, col_ind1= lsa(cost_matrix)
print("===LSA sol: ")
cost_tot = 0
for id in range (len(row_ind)):
    i = row_ind1[id]
    j = col_ind1[id]
    print (f"track {i} match det {j} cost {cost_matrix[i,j]}")

    cost_tot += cost_matrix[i,j]
print (f"cost total = {cost_tot}")
weisha1991 commented 10 months ago
`/// @brief Exact Jonker-Volgenant algorithm.
/// @param dim in problem size
/// @param assign_cost in cost matrix
/// @param verbose in indicates whether to report the progress to stdout
/// @param rowsol out column assigned to row in solution / size dim
/// @param colsol out row assigned to column in solution / size dim
/// @param u out dual variables, row reduction numbers / size dim
/// @param v out dual variables, column reduction numbers / size dim
/// @return achieved minimum assignment cost
template <bool avx2, typename idx, typename cost>
cost lap(int dim, **const cost *restrict assign_cost**, bool verbose,
         idx *restrict rowsol, idx *restrict colsol,
         cost *restrict u, cost *restrict v) ;

@mayejudi @CSautier Please be careful when call the funciton lap, the second args assign_cost should be 1- dimensional array