lap.cpp
Go to the documentation of this file.
00001 /************************************************************************
00002 *
00003 *  lap.cpp
00004    version 1.0 - 4 September 1996
00005    author: Roy Jonker @ MagicLogic Optimization Inc.
00006    e-mail: roy_jonker@magiclogic.com
00007 
00008    Code for Linear Assignment Problem, according to
00009 
00010    "A Shortest Augmenting Path Algorithm for Dense and Sparse Linear
00011     Assignment Problems," Computing 38, 325-340, 1987
00012 
00013    by
00014 
00015    R. Jonker and A. Volgenant, University of Amsterdam.
00016 *
00017 *************************************************************************/
00018 
00019 #include <stdio.h>
00020 // Replacing system with better functions to not shadow things
00021 #include "gnrl.h"
00022 #include "lap.h"
00023 
00024 namespace cpl_visual_features
00025 {
00026 LapCost lap(int dim,
00027              LapCost **assigncost,
00028              LapCol *rowsol,
00029              LapRow *colsol,
00030              LapCost *u,
00031              LapCost *v)
00032 
00033 // input:
00034 // dim        - problem size
00035 // assigncost - cost matrix
00036 
00037 // output:
00038 // rowsol     - column assigned to row in solution
00039 // colsol     - row assigned to column in solution
00040 // u          - dual variables, row reduction numbers
00041 // v          - dual variables, column reduction numbers
00042 
00043 {
00044   boolean unassignedfound;
00045   LapRow  i, imin, numfree = 0, prvnumfree, f, i0, k, freerow, *pred, *free;
00046   LapCol  j, j1, j2, endofpath, last, low, up, *collist, *matches;
00047   LapCost min, h, umin, usubmin, v2, *d;
00048 
00049   free = new LapRow[dim];       // list of unassigned rows.
00050   collist = new LapCol[dim];    // list of columns to be scanned in various ways.
00051   matches = new LapCol[dim];    // counts how many times a row could be assigned.
00052   d = new LapCost[dim];         // 'cost-distance' in augmenting path calculation.
00053   pred = new LapRow[dim];       // row-predecessor of column in augmenting/alternating path.
00054 
00055   // init how many times a row will be assigned in the column reduction.
00056   for (i = 0; i < dim; i++)
00057     matches[i] = 0;
00058 
00059   // COLUMN REDUCTION
00060   for (j = dim-1; j >= 0; j--)    // reverse order gives better results.
00061   {
00062     // find minimum cost over rows.
00063     min = assigncost[0][j];
00064     imin = 0;
00065     for (i = 1; i < dim; i++)
00066       if (assigncost[i][j] < min)
00067       {
00068         min = assigncost[i][j];
00069         imin = i;
00070       }
00071     v[j] = min;
00072 
00073     if (++matches[imin] == 1)
00074     {
00075       // init assignment if minimum row assigned for first time.
00076       rowsol[imin] = j;
00077       colsol[j] = imin;
00078     }
00079     else
00080       colsol[j] = -1;        // row already assigned, column not assigned.
00081   }
00082 
00083   // REDUCTION TRANSFER
00084   for (i = 0; i < dim; i++)
00085     if (matches[i] == 0)     // fill list of unassigned 'free' rows.
00086       free[numfree++] = i;
00087     else
00088       if (matches[i] == 1)   // transfer reduction from rows that are assigned once.
00089       {
00090         j1 = rowsol[i];
00091         min = BIG;
00092         for (j = 0; j < dim; j++)
00093           if (j != j1)
00094             if (assigncost[i][j] - v[j] < min)
00095               min = assigncost[i][j] - v[j];
00096         v[j1] = v[j1] - min;
00097       }
00098 
00099   // AUGMENTING ROW REDUCTION
00100   int loopcnt = 0;           // do-loop to be done twice.
00101   do
00102   {
00103     loopcnt++;
00104 
00105     // scan all free rows.
00106     // in some cases, a free row may be replaced with another one to be scanned next.
00107     k = 0;
00108     prvnumfree = numfree;
00109     numfree = 0;             // start list of rows still free after augmenting row reduction.
00110     while (k < prvnumfree)
00111     {
00112       i = free[k];
00113       k++;
00114 
00115       // find minimum and second minimum reduced cost over columns.
00116       umin = assigncost[i][0] - v[0];
00117       j1 = 0;
00118       usubmin = BIG;
00119       for (j = 1; j < dim; j++)
00120       {
00121         h = assigncost[i][j] - v[j];
00122         if (h < usubmin)
00123         {
00124           if (h >= umin)
00125           {
00126             usubmin = h;
00127             j2 = j;
00128           }
00129           else
00130           {
00131             usubmin = umin;
00132             umin = h;
00133             j2 = j1;
00134             j1 = j;
00135           }
00136         }
00137       }
00138 
00139       i0 = colsol[j1];
00140       if (umin < usubmin)
00141         // change the reduction of the minimum column to increase the minimum
00142         // reduced cost in the row to the subminimum.
00143         v[j1] = v[j1] - (usubmin - umin);
00144       else                   // minimum and subminimum equal.
00145         if (i0 >= 0)         // minimum column j1 is assigned.
00146         {
00147           // swap columns j1 and j2, as j2 may be unassigned.
00148           j1 = j2;
00149           i0 = colsol[j2];
00150         }
00151 
00152       // (re-)assign i to j1, possibly de-assigning an i0.
00153       rowsol[i] = j1;
00154       colsol[j1] = i;
00155 
00156       if (i0 >= 0)           // minimum column j1 assigned earlier.
00157       {
00158         if (umin < usubmin)
00159           // put in current k, and go back to that k.
00160           // continue augmenting path i - j1 with i0.
00161           free[--k] = i0;
00162         else
00163           // no further augmenting reduction possible.
00164           // store i0 in list of free rows for next phase.
00165           free[numfree++] = i0;
00166       }
00167     }
00168   }
00169   while (loopcnt < 2);       // repeat once.
00170 
00171   // AUGMENT SOLUTION for each free row.
00172   for (f = 0; f < numfree; f++)
00173   {
00174     freerow = free[f];       // start row of augmenting path.
00175 
00176     // Dijkstra shortest path algorithm.
00177     // runs until unassigned column added to shortest path tree.
00178     for (j = 0; j < dim; j++)
00179     {
00180       d[j] = assigncost[freerow][j] - v[j];
00181       pred[j] = freerow;
00182       collist[j] = j;        // init column list.
00183     }
00184 
00185     low = 0; // columns in 0..low-1 are ready, now none.
00186     up = 0;  // columns in low..up-1 are to be scanned for current minimum, now none.
00187              // columns in up..dim-1 are to be considered later to find new minimum,
00188              // at this stage the list simply contains all columns
00189     unassignedfound = FALSE;
00190     do
00191     {
00192       if (up == low)         // no more columns to be scanned for current minimum.
00193       {
00194         last = low - 1;
00195 
00196         // scan columns for up..dim-1 to find all indices for which new minimum occurs.
00197         // store these indices between low..up-1 (increasing up).
00198         min = d[collist[up++]];
00199         for (k = up; k < dim; k++)
00200         {
00201           j = collist[k];
00202           h = d[j];
00203           if (h <= min)
00204           {
00205             if (h < min)     // new minimum.
00206             {
00207               up = low;      // restart list at index low.
00208               min = h;
00209             }
00210             // new index with same minimum, put on undex up, and extend list.
00211             collist[k] = collist[up];
00212             collist[up++] = j;
00213           }
00214         }
00215 
00216         // check if any of the minimum columns happens to be unassigned.
00217         // if so, we have an augmenting path right away.
00218         for (k = low; k < up; k++)
00219           if (colsol[collist[k]] < 0)
00220           {
00221             endofpath = collist[k];
00222             unassignedfound = TRUE;
00223             break;
00224           }
00225       }
00226 
00227       if (!unassignedfound)
00228       {
00229         // update 'distances' between freerow and all unscanned columns, via next scanned column.
00230         j1 = collist[low];
00231         low++;
00232         i = colsol[j1];
00233         h = assigncost[i][j1] - v[j1] - min;
00234 
00235         for (k = up; k < dim; k++)
00236         {
00237           j = collist[k];
00238           v2 = assigncost[i][j] - v[j] - h;
00239           if (v2 < d[j])
00240           {
00241             pred[j] = i;
00242             if (v2 == min)   // new column found at same minimum value
00243             {
00244               if (colsol[j] < 0)
00245               {
00246                 // if unassigned, shortest augmenting path is complete.
00247                 endofpath = j;
00248                 unassignedfound = TRUE;
00249                 break;
00250               }
00251               // else add to list to be scanned right away.
00252               else
00253               {
00254                 collist[k] = collist[up];
00255                 collist[up++] = j;
00256               }
00257             d[j] = v2;
00258             }
00259           }
00260         }
00261       }
00262     }
00263     while (!unassignedfound);
00264 
00265     // update column prices.
00266     for (k = 0; k <= last; k++)
00267     {
00268       j1 = collist[k];
00269       v[j1] = v[j1] + d[j1] - min;
00270     }
00271 
00272     // reset row and column assignments along the alternating path.
00273     do
00274     {
00275       i = pred[endofpath];
00276       colsol[endofpath] = i;
00277       j1 = endofpath;
00278       endofpath = rowsol[i];
00279       rowsol[i] = j1;
00280     }
00281     while (i != freerow);
00282   }
00283 
00284   // calculate optimal cost.
00285   LapCost lapcost = 0;
00286   for (i = 0; i < dim; i++)
00287   {
00288     j = rowsol[i];
00289     u[i] = assigncost[i][j] - v[j];
00290     lapcost = lapcost + assigncost[i][j];
00291   }
00292 
00293   // free reserved memory.
00294   delete[] pred;
00295   delete[] free;
00296   delete[] collist;
00297   delete[] matches;
00298   delete[] d;
00299 
00300   return lapcost;
00301 }
00302 
00303 void checklap(int dim, LapCost **assigncost,
00304               LapCol *rowsol, LapRow *colsol, LapCost *u, LapCost *v)
00305 {
00306   LapRow  i;
00307   LapCol  j;
00308   LapCost lapcost = 0, redcost = 0;
00309   boolean *matched;
00310   char wait;
00311 
00312   matched = new boolean[dim];
00313 
00314   for (i = 0; i < dim; i++)
00315     for (j = 0; j < dim; j++)
00316       if ((redcost = assigncost[i][j] - u[i] - v[j]) < 0)
00317       {
00318         printf("\n");
00319         printf("negative reduced cost i %d j %d redcost %d\n", i, j, redcost);
00320         printf("\n\ndim %5d - press key\n", dim);
00321         scanf("%d", &wait);
00322         break;
00323       }
00324 
00325   for (i = 0; i < dim; i++)
00326     if ((redcost = assigncost[i][rowsol[i]] - u[i] - v[rowsol[i]]) != 0)
00327     {
00328       printf("\n");
00329       printf("non-null reduced cost i %d soli %d redcost %d\n", i, rowsol[i], redcost);
00330       printf("\n\ndim %5d - press key\n", dim);
00331       scanf("%d", &wait);
00332       break;
00333     }
00334 
00335   for (j = 0; j < dim; j++)
00336     matched[j] = FALSE;
00337 
00338   for (i = 0; i < dim; i++)
00339     if (matched[rowsol[i]])
00340     {
00341       printf("\n");
00342       printf("column matched more than once - i %d soli %d\n", i, rowsol[i]);
00343       printf("\n\ndim %5d - press key\n", dim);
00344       scanf("%d", &wait);
00345       break;
00346     }
00347     else
00348       matched[rowsol[i]] = TRUE;
00349 
00350 
00351   for (i = 0; i < dim; i++)
00352     if (colsol[rowsol[i]] != i)
00353     {
00354       printf("\n");
00355       printf("error in row solution i %d soli %d solsoli %d\n", i, rowsol[i], colsol[rowsol[i]]);
00356       printf("\n\ndim %5d - press key\n", dim);
00357       scanf("%d", &wait);
00358       break;
00359     }
00360 
00361   for (j = 0; j < dim; j++)
00362     if (rowsol[colsol[j]] != j)
00363     {
00364       printf("\n");
00365       printf("error in col solution j %d solj %d solsolj %d\n", j, colsol[j], rowsol[colsol[j]]);
00366       printf("\n\ndim %5d - press key\n", dim);
00367       scanf("%d", &wait);
00368       break;
00369     }
00370 
00371   delete[] matched;
00372   return;
00373 }
00374 };


cpl_visual_features
Author(s): Tucker Hermans
autogenerated on Wed Nov 27 2013 11:52:36