Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

LAPJV is not giving the correct result. #83

Open
mayejudi opened this issue Oct 13, 2023 · 5 comments
Open

LAPJV is not giving the correct result. #83

mayejudi opened this issue Oct 13, 2023 · 5 comments

Comments

@mayejudi
Copy link

mayejudi commented Oct 13, 2023

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
Copy link
Collaborator

vmarkovtsev commented Oct 13, 2023

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

  • There should be an environment variable or a boolean flag in lapjv() to print logs just like in your custom cpp. You can compare both.
  • The third returned value should be the total cost or something else useful. Check it out.
  • As you see the Python code assigned X to X, which is suspicious. You can debug the flow with gdb if the logs are not enough.
  • Check with float32.

@mayejudi
Copy link
Author

mayejudi commented Oct 13, 2023

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
Copy link
Author

mayejudi commented Oct 14, 2023

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
Copy link

For future reference, the results of scipy's LSA and of this repo's output have different formatting (see gatagat/lap#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
Copy link

weisha1991 commented Jan 19, 2024

`/// @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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants