Main Page Namespace List Class Hierarchy Compound List File List Namespace Members Compound Members File Members Related Pages

# Emd.cpp

    00001 /*
    00002
    00003     This code has been modified so that if works for our feature type. -Fatih
    00004
    00005     An implementation of the Earth Movers Distance.
    00006     Based of the solution for the Transportation problem as described in
    00007     "Introduction to Mathematical Programming" by F. S. Hillier and
    00008     G. J. Lieberman, McGraw-Hill, 1990.
    00009
    00010     Copyright (C) 1998 Yossi Rubner
    00011     Computer Science Department, Stanford University
    00012     E-Mail: rubner@cs.stanford.edu   URL: http://vision.stanford.edu/~rubner
    00013 */
    00014
    00015
    00016 #include <stdio.h>
    00017 #include <stdlib.h>
    00018 #include <math.h>
    00019
    00020 #include "emd.h"
    00021
    00022 #define DEBUG_LEVEL 0
    00023 /*
    00024  DEBUG_LEVEL:
    00025    0 = NO MESSAGES
    00026    1 = PRINT THE NUMBER OF ITERATIONS AND THE FINAL RESULT
    00027    2 = PRINT THE RESULT AFTER EVERY ITERATION
    00028    3 = PRINT ALSO THE FLOW AFTER EVERY ITERATION
    00029    4 = PRINT A LOT OF INFORMATION (PROBABLY USEFUL ONLY FOR THE AUTHOR)
    00030 */
    00031
    00032
    00033 #define MAX_SIG_SIZE1  (MAX_SIG_SIZE+1)  /* FOR THE POSIBLE DUMMY FEATURE */
    00034
    00035 /* NEW TYPES DEFINITION */
    00036
    00037 /* node1_t IS USED FOR SINGLE-LINKED LISTS */
    00038 typedef struct node1_t {
    00039   int i;
    00040   double val;
    00041   struct node1_t *Next;
    00042 } node1_t;
    00043
    00044 /* node1_t IS USED FOR DOUBLE-LINKED LISTS */
    00045 typedef struct node2_t {
    00046   int i, j;
    00047   double val;
    00048   struct node2_t *NextC;               /* NEXT COLUMN */
    00049   struct node2_t *NextR;               /* NEXT ROW */
    00050 } node2_t;
    00051
    00052
    00053
    00054 /* GLOBAL VARIABLE DECLARATION */
    00055 static int _n1, _n2;                          /* SIGNATURES SIZES */
    00056 static float _C[MAX_SIG_SIZE1][MAX_SIG_SIZE1];/* THE COST MATRIX */
    00057 static node2_t _X[MAX_SIG_SIZE1*2];            /* THE BASIC VARIABLES VECTOR */
    00058 /* VARIABLES TO HANDLE _X EFFICIENTLY */
    00059 static node2_t *_EndX, *_EnterX;
    00060 static char _IsX[MAX_SIG_SIZE1][MAX_SIG_SIZE1];
    00061 static node2_t *_RowsX[MAX_SIG_SIZE1], *_ColsX[MAX_SIG_SIZE1];
    00062 static double _maxW;
    00063 static float _maxC;
    00064
    00065 /* DECLARATION OF FUNCTIONS */
    00066 static float init(signature_t *Signature1, signature_t *Signature2,
    00067                   float (*Dist)(feature_t *, feature_t *));
    00068 static void findBasicVariables(node1_t *U, node1_t *V);
    00069 static int isOptimal(node1_t *U, node1_t *V);
    00070 static int findLoop(node2_t **Loop);
    00071 static void newSol();
    00072 static void russel(double *S, double *D);
    00073 static void addBasicVariable(int minI, int minJ, double *S, double *D,
    00074                              node1_t *PrevUMinI, node1_t *PrevVMinJ,
    00075                              node1_t *UHead);
    00076 #if DEBUG_LEVEL > 0
    00077 static void printSolution();
    00078 #endif
    00079
    00080
    00081 /******************************************************************************
    00082 double emd(signature_t *Signature1, signature_t *Signature2,
    00083           float (*Dist)(feature_t *, feature_t *), flow_t *Flow, int *FlowSize)
    00084
    00085 where
    00086
    00087    Signature1, Signature2  Pointers to signatures that their distance we want
    00088               to compute.
    00089    Dist       Pointer to the ground distance. i.e. the function that computes
    00090               the distance between two features.
    00091    Flow       (Optional) Pointer to a vector of flow_t (defined in emd.h)
    00092               where the resulting flow will be stored. Flow must have n1+n2-1
    00093               elements, where n1 and n2 are the sizes of the two signatures
    00094               respectively.
    00095               If NULL, the flow is not returned.
    00096    FlowSize   (Optional) Pointer to an integer where the number of elements in
    00097               Flow will be stored
    00098
    00099 ******************************************************************************/
    00100
    00101 double emd(signature_t *Signature1, signature_t *Signature2,
    00102           float (*Dist)(feature_t *, feature_t *),
    00103           flow_t *Flow, int *FlowSize)
    00104 {
    00105   int itr;
    00106   double totalCost;
    00107   float w;
    00108   node2_t *XP;
    00109   flow_t *FlowP;
    00110   node1_t U[MAX_SIG_SIZE1], V[MAX_SIG_SIZE1];
    00111
    00112   w = init(Signature1, Signature2, Dist);
    00113
    00114 #if DEBUG_LEVEL > 1
    00115   printf("\nINITIAL SOLUTION:\n");
    00116   printSolution();
    00117 #endif
    00118
    00119   if (_n1 > 1 && _n2 > 1)  /* IF _n1 = 1 OR _n2 = 1 THEN WE ARE DONE */
    00120     {
    00121       for (itr = 1; itr < MAX_ITERATIONS; itr++)
    00122         {
    00123           /* FIND BASIC VARIABLES */
    00124           findBasicVariables(U, V);
    00125
    00126           /* CHECK FOR OPTIMALITY */
    00127           if (isOptimal(U, V))
    00128             break;
    00129
    00130           /* IMPROVE SOLUTION */
    00131           newSol();
    00132
    00133 #if DEBUG_LEVEL > 1
    00134           printf("\nITERATION # %d \n", itr);
    00135           printSolution();
    00136 #endif
    00137         }
    00138
    00139       if (itr == MAX_ITERATIONS)
    00140         fprintf(stderr, "emd: Maximum number of iterations has been reached (%d)\n",
    00141                 MAX_ITERATIONS);
    00142     }
    00143
    00144   /* COMPUTE THE TOTAL FLOW */
    00145   totalCost = 0;
    00146   if (Flow != NULL)
    00147     FlowP = Flow;
    00148   for(XP=_X; XP < _EndX; XP++)
    00149     {
    00150
    00151         /***************************/
    00152       if (XP == _EnterX)  /* _EnterX IS THE EMPTY SLOT */
    00153         continue;
    00154       if (XP->i == Signature1->n || XP->j == Signature2->n)  /* DUMMY FEATURE */
    00155         continue;
    00156
    00157       if (XP->val == 0)  /* ZERO FLOW */
    00158         continue;
    00159
    00160       totalCost += (double)XP->val * _C[XP->i][XP->j];
    00161       if (Flow != NULL)
    00162         {
    00163           FlowP->from = XP->i;
    00164           FlowP->to = XP->j;
    00165           FlowP->amount = XP->val;
    00166           FlowP++;
    00167         }
    00168     }
    00169   if (Flow != NULL)
    00170     *FlowSize = FlowP-Flow;
    00171
    00172 #if DEBUG_LEVEL > 0
    00173   printf("\n*** OPTIMAL SOLUTION (%d ITERATIONS): %f ***\n", itr, totalCost);
    00174 #endif
    00175
    00176   /* RETURN THE NORMALIZED COST == EMD */
    00177   return (float)(totalCost / w);
    00178 }
    00179
    00180
    00181
    00182
    00183
    00184 /**********************
    00185    init
    00186 **********************/
    00187 static float init(signature_t *Signature1, signature_t *Signature2,
    00188                   float (*Dist)(feature_t *, feature_t *))
    00189 {
    00190   int i, j;
    00191   double sSum, dSum, diff;
    00192   feature_t *P1, *P2;
    00193   double S[MAX_SIG_SIZE1], D[MAX_SIG_SIZE1];
    00194
    00195   _n1 = Signature1->n;
    00196   _n2 = Signature2->n;
    00197
    00198   if (_n1 > MAX_SIG_SIZE || _n2 > MAX_SIG_SIZE)
    00199     {
    00200       fprintf(stderr, "emd: Signature size is limited to %d\n", MAX_SIG_SIZE);
    00201       exit(1);
    00202     }
    00203
    00204   /* COMPUTE THE DISTANCE MATRIX */
    00205   _maxC = 0;
    00206   for(i=0, P1=Signature1->Features; i < _n1; i++, P1++)
    00207     for(j=0, P2=Signature2->Features; j < _n2; j++, P2++)
    00208       {
    00209         _C[i][j] = Dist(P1, P2);
    00210         if (_C[i][j] > _maxC)
    00211           _maxC = _C[i][j];
    00212       }
    00213
    00214   /* SUM UP THE SUPPLY AND DEMAND */
    00215   sSum = 0.0;
    00216   for(i=0; i < _n1; i++)
    00217     {
    00218       S[i] = Signature1->Weights[i];
    00219       sSum += Signature1->Weights[i];
    00220       _RowsX[i] = NULL;
    00221     }
    00222   dSum = 0.0;
    00223   for(j=0; j < _n2; j++)
    00224     {
    00225       D[j] = Signature2->Weights[j];
    00226       dSum += Signature2->Weights[j];
    00227       _ColsX[j] = NULL;
    00228     }
    00229
    00230   /* IF SUPPLY DIFFERENT THAN THE DEMAND, ADD A ZERO-COST DUMMY CLUSTER */
    00231   diff = sSum - dSum;
    00232   if (fabs(diff) >= EPSILON * sSum)
    00233     {
    00234       if (diff < 0.0)
    00235         {
    00236           for (j=0; j < _n2; j++)
    00237             _C[_n1][j] = 0;
    00238           S[_n1] = -diff;
    00239           _RowsX[_n1] = NULL;
    00240           _n1++;
    00241         }
    00242       else
    00243         {
    00244           for (i=0; i < _n1; i++)
    00245             _C[i][_n2] = 0;
    00246           D[_n2] = diff;
    00247           _ColsX[_n2] = NULL;
    00248           _n2++;
    00249         }
    00250     }
    00251
    00252   /* INITIALIZE THE BASIC VARIABLE STRUCTURES */
    00253   for (i=0; i < _n1; i++)
    00254     for (j=0; j < _n2; j++)
    00255         _IsX[i][j] = 0;
    00256   _EndX = _X;
    00257
    00258   _maxW = sSum > dSum ? sSum : dSum;
    00259
    00260   /* FIND INITIAL SOLUTION */
    00261   russel(S, D);
    00262
    00263   _EnterX = _EndX++;  /* AN EMPTY SLOT (ONLY _n1+_n2-1 BASIC VARIABLES) */
    00264
    00265   return sSum > dSum ? dSum : sSum;
    00266 }
    00267
    00268
    00269 /**********************
    00270     findBasicVariables
    00271  **********************/
    00272 static void findBasicVariables(node1_t *U, node1_t *V)
    00273 {
    00274   int i, j, found;
    00275   int UfoundNum, VfoundNum;
    00276   node1_t u0Head, u1Head, *CurU, *PrevU;
    00277   node1_t v0Head, v1Head, *CurV, *PrevV;
    00278
    00279   /* INITIALIZE THE ROWS LIST (U) AND THE COLUMNS LIST (V) */
    00280   u0Head.Next = CurU = U;
    00281   for (i=0; i < _n1; i++)
    00282     {
    00283       CurU->i = i;
    00284       CurU->Next = CurU+1;
    00285       CurU++;
    00286     }
    00287   (--CurU)->Next = NULL;
    00288   u1Head.Next = NULL;
    00289
    00290   CurV = V+1;
    00291   v0Head.Next = _n2 > 1 ? V+1 : NULL;
    00292   for (j=1; j < _n2; j++)
    00293     {
    00294       CurV->i = j;
    00295       CurV->Next = CurV+1;
    00296       CurV++;
    00297     }
    00298   (--CurV)->Next = NULL;
    00299   v1Head.Next = NULL;
    00300
    00301   /* THERE ARE _n1+_n2 VARIABLES BUT ONLY _n1+_n2-1 INDEPENDENT EQUATIONS,
    00302      SO SET V[0]=0 */
    00303   V[0].i = 0;
    00304   V[0].val = 0;
    00305   v1Head.Next = V;
    00306   v1Head.Next->Next = NULL;
    00307
    00308   /* LOOP UNTIL ALL VARIABLES ARE FOUND */
    00309   UfoundNum=VfoundNum=0;
    00310   while (UfoundNum < _n1 || VfoundNum < _n2)
    00311     {
    00312
    00313 #if DEBUG_LEVEL > 3
    00314       printf("UfoundNum=%d/%d,VfoundNum=%d/%d\n",UfoundNum,_n1,VfoundNum,_n2);
    00315       printf("U0=");
    00316       for(CurU = u0Head.Next; CurU != NULL; CurU = CurU->Next)
    00317         printf("[%d]",CurU-U);
    00318       printf("\n");
    00319       printf("U1=");
    00320       for(CurU = u1Head.Next; CurU != NULL; CurU = CurU->Next)
    00321         printf("[%d]",CurU-U);
    00322       printf("\n");
    00323       printf("V0=");
    00324       for(CurV = v0Head.Next; CurV != NULL; CurV = CurV->Next)
    00325         printf("[%d]",CurV-V);
    00326       printf("\n");
    00327       printf("V1=");
    00328       for(CurV = v1Head.Next; CurV != NULL; CurV = CurV->Next)
    00329         printf("[%d]",CurV-V);
    00330       printf("\n\n");
    00331 #endif
    00332
    00333       found = 0;
    00334       if (VfoundNum < _n2)
    00335         {
    00336           /* LOOP OVER ALL MARKED COLUMNS */
    00337           PrevV = &v1Head;
    00338           for (CurV=v1Head.Next; CurV != NULL; CurV=CurV->Next)
    00339             {
    00340               j = CurV->i;
    00341               /* FIND THE VARIABLES IN COLUMN j */
    00342               PrevU = &u0Head;
    00343               for (CurU=u0Head.Next; CurU != NULL; CurU=CurU->Next)
    00344                 {
    00345                   i = CurU->i;
    00346                   if (_IsX[i][j])
    00347                     {
    00348                       /* COMPUTE U[i] */
    00349                       CurU->val = _C[i][j] - CurV->val;
    00350                       /* ...AND ADD IT TO THE MARKED LIST */
    00351                       PrevU->Next = CurU->Next;
    00352                       CurU->Next = u1Head.Next != NULL ? u1Head.Next : NULL;
    00353                       u1Head.Next = CurU;
    00354                       CurU = PrevU;
    00355                     }
    00356                   else
    00357                     PrevU = CurU;
    00358                 }
    00359               PrevV->Next = CurV->Next;
    00360               VfoundNum++;
    00361               found = 1;
    00362             }
    00363         }
    00364      if (UfoundNum < _n1)
    00365         {
    00366           /* LOOP OVER ALL MARKED ROWS */
    00367           PrevU = &u1Head;
    00368           for (CurU=u1Head.Next; CurU != NULL; CurU=CurU->Next)
    00369             {
    00370               i = CurU->i;
    00371               /* FIND THE VARIABLES IN ROWS i */
    00372               PrevV = &v0Head;
    00373               for (CurV=v0Head.Next; CurV != NULL; CurV=CurV->Next)
    00374                 {
    00375                   j = CurV->i;
    00376                   if (_IsX[i][j])
    00377                     {
    00378                       /* COMPUTE V[j] */
    00379                       CurV->val = _C[i][j] - CurU->val;
    00380                       /* ...AND ADD IT TO THE MARKED LIST */
    00381                       PrevV->Next = CurV->Next;
    00382                       CurV->Next = v1Head.Next != NULL ? v1Head.Next: NULL;
    00383                       v1Head.Next = CurV;
    00384                       CurV = PrevV;
    00385                     }
    00386                   else
    00387                     PrevV = CurV;
    00388                 }
    00389               PrevU->Next = CurU->Next;
    00390               UfoundNum++;
    00391               found = 1;
    00392             }
    00393         }
    00394      if (! found)
    00395        {
    00396          fprintf(stderr, "emd: Unexpected error in findBasicVariables!\n");
    00397          fprintf(stderr, "This typically happens when the EPSILON defined in\n");
    00398          fprintf(stderr, "emd.h is not right for the scale of the problem.\n");
    00399          exit(1);
    00400        }
    00401     }
    00402 }
    00403
    00404
    00405
    00406
    00407 /**********************
    00408     isOptimal
    00409  **********************/
    00410 static int isOptimal(node1_t *U, node1_t *V)
    00411 {
    00412   double delta, deltaMin;
    00413   int i, j, minI, minJ;
    00414
    00415   /* FIND THE MINIMAL Cij-Ui-Vj OVER ALL i,j */
    00416   deltaMin = INFINITY;
    00417   for(i=0; i < _n1; i++)
    00418     for(j=0; j < _n2; j++)
    00419       if (! _IsX[i][j])
    00420         {
    00421           delta = _C[i][j] - U[i].val - V[j].val;
    00422           if (deltaMin > delta)
    00423             {
    00424               deltaMin = delta;
    00425               minI = i;
    00426               minJ = j;
    00427             }
    00428         }
    00429
    00430 #if DEBUG_LEVEL > 3
    00431   printf("deltaMin=%f\n", deltaMin);
    00432 #endif
    00433
    00434    if (deltaMin == INFINITY)
    00435      {
    00436        fprintf(stderr, "emd: Unexpected error in isOptimal.\n");
    00437        exit(0);
    00438      }
    00439
    00440    _EnterX->i = minI;
    00441    _EnterX->j = minJ;
    00442
    00443    /* IF NO NEGATIVE deltaMin, WE FOUND THE OPTIMAL SOLUTION */
    00444    return deltaMin >= -EPSILON * _maxC;
    00445
    00446 /*
    00447    return deltaMin >= -EPSILON;
    00448  */
    00449 }
    00450
    00451
    00452 /**********************
    00453     newSol
    00454 **********************/
    00455 static void newSol()
    00456 {
    00457     int i, j, k;
    00458     double xMin;
    00459     int steps;
    00460     node2_t *Loop[2*MAX_SIG_SIZE1], *CurX, *LeaveX;
    00461
    00462 #if DEBUG_LEVEL > 3
    00463     printf("EnterX = (%d,%d)\n", _EnterX->i, _EnterX->j);
    00464 #endif
    00465
    00466     /* ENTER THE NEW BASIC VARIABLE */
    00467     i = _EnterX->i;
    00468     j = _EnterX->j;
    00469     _IsX[i][j] = 1;
    00470     _EnterX->NextC = _RowsX[i];
    00471     _EnterX->NextR = _ColsX[j];
    00472     _EnterX->val = 0;
    00473     _RowsX[i] = _EnterX;
    00474     _ColsX[j] = _EnterX;
    00475
    00476     /* FIND A CHAIN REACTION */
    00477     steps = findLoop(Loop);
    00478
    00479     /* FIND THE LARGEST VALUE IN THE LOOP */
    00480     xMin = INFINITY;
    00481     for (k=1; k < steps; k+=2)
    00482       {
    00483         if (Loop[k]->val < xMin)
    00484           {
    00485             LeaveX = Loop[k];
    00486             xMin = Loop[k]->val;
    00487           }
    00488       }
    00489
    00490     /* UPDATE THE LOOP */
    00491     for (k=0; k < steps; k+=2)
    00492       {
    00493         Loop[k]->val += xMin;
    00494         Loop[k+1]->val -= xMin;
    00495       }
    00496
    00497 #if DEBUG_LEVEL > 3
    00498     printf("LeaveX = (%d,%d)\n", LeaveX->i, LeaveX->j);
    00499 #endif
    00500
    00501     /* REMOVE THE LEAVING BASIC VARIABLE */
    00502     i = LeaveX->i;
    00503     j = LeaveX->j;
    00504     _IsX[i][j] = 0;
    00505     if (_RowsX[i] == LeaveX)
    00506       _RowsX[i] = LeaveX->NextC;
    00507     else
    00508       for (CurX=_RowsX[i]; CurX != NULL; CurX = CurX->NextC)
    00509         if (CurX->NextC == LeaveX)
    00510           {
    00511             CurX->NextC = CurX->NextC->NextC;
    00512             break;
    00513           }
    00514     if (_ColsX[j] == LeaveX)
    00515       _ColsX[j] = LeaveX->NextR;
    00516     else
    00517       for (CurX=_ColsX[j]; CurX != NULL; CurX = CurX->NextR)
    00518         if (CurX->NextR == LeaveX)
    00519           {
    00520             CurX->NextR = CurX->NextR->NextR;
    00521             break;
    00522           }
    00523
    00524     /* SET _EnterX TO BE THE NEW EMPTY SLOT */
    00525     _EnterX = LeaveX;
    00526 }
    00527
    00528
    00529
    00530 /**********************
    00531     findLoop
    00532 **********************/
    00533 static int findLoop(node2_t **Loop)
    00534 {
    00535   int i, steps;
    00536   node2_t **CurX, *NewX;
    00537   char IsUsed[2*MAX_SIG_SIZE1];
    00538
    00539   for (i=0; i < _n1+_n2; i++)
    00540     IsUsed[i] = 0;
    00541
    00542   CurX = Loop;
    00543   NewX = *CurX = _EnterX;
    00544   IsUsed[_EnterX-_X] = 1;
    00545   steps = 1;
    00546
    00547   do
    00548     {
    00549       if (steps%2 == 1)
    00550         {
    00551           /* FIND AN UNUSED X IN THE ROW */
    00552           NewX = _RowsX[NewX->i];
    00553           while (NewX != NULL && IsUsed[NewX-_X])
    00554             NewX = NewX->NextC;
    00555         }
    00556       else
    00557         {
    00558           /* FIND AN UNUSED X IN THE COLUMN, OR THE ENTERING X */
    00559           NewX = _ColsX[NewX->j];
    00560           while (NewX != NULL && IsUsed[NewX-_X] && NewX != _EnterX)
    00561             NewX = NewX->NextR;
    00562           if (NewX == _EnterX)
    00563             break;
    00564         }
    00565
    00566      if (NewX != NULL)  /* FOUND THE NEXT X */
    00567        {
    00568          /* ADD X TO THE LOOP */
    00569          *++CurX = NewX;
    00570          IsUsed[NewX-_X] = 1;
    00571          steps++;
    00572 #if DEBUG_LEVEL > 3
    00573          printf("steps=%d, NewX=(%d,%d)\n", steps, NewX->i, NewX->j);
    00574 #endif
    00575        }
    00576      else  /* DIDN'T FIND THE NEXT X */
    00577        {
    00578          /* BACKTRACK */
    00579          do
    00580            {
    00581              NewX = *CurX;
    00582              do
    00583                {
    00584                  if (steps%2 == 1)
    00585                    NewX = NewX->NextR;
    00586                  else
    00587                    NewX = NewX->NextC;
    00588                } while (NewX != NULL && IsUsed[NewX-_X]);
    00589
    00590              if (NewX == NULL)
    00591                {
    00592                  IsUsed[*CurX-_X] = 0;
    00593                  CurX--;
    00594                  steps--;
    00595                }
    00596          } while (NewX == NULL && CurX >= Loop);
    00597
    00598 #if DEBUG_LEVEL > 3
    00599          printf("BACKTRACKING TO: steps=%d, NewX=(%d,%d)\n",
    00600                 steps, NewX->i, NewX->j);
    00601 #endif
    00602            IsUsed[*CurX-_X] = 0;
    00603            *CurX = NewX;
    00604            IsUsed[NewX-_X] = 1;
    00605        }
    00606     } while(CurX >= Loop);
    00607
    00608   if (CurX == Loop)
    00609     {
    00610       fprintf(stderr, "emd: Unexpected error in findLoop!\n");
    00611       exit(1);
    00612     }
    00613 #if DEBUG_LEVEL > 3
    00614   printf("FOUND LOOP:\n");
    00615   for (i=0; i < steps; i++)
    00616     printf("%d: (%d,%d)\n", i, Loop[i]->i, Loop[i]->j);
    00617 #endif
    00618
    00619   return steps;
    00620 }
    00621
    00622
    00623
    00624 /**********************
    00625     russel
    00626 **********************/
    00627 static void russel(double *S, double *D)
    00628 {
    00629   int i, j, found, minI, minJ;
    00630   double deltaMin, oldVal, diff;
    00631   double Delta[MAX_SIG_SIZE1][MAX_SIG_SIZE1];
    00632   node1_t Ur[MAX_SIG_SIZE1], Vr[MAX_SIG_SIZE1];
    00633   node1_t uHead, *CurU, *PrevU;
    00634   node1_t vHead, *CurV, *PrevV;
    00635   node1_t *PrevUMinI, *PrevVMinJ, *Remember;
    00636
    00637   /* INITIALIZE THE ROWS LIST (Ur), AND THE COLUMNS LIST (Vr) */
    00638   uHead.Next = CurU = Ur;
    00639   for (i=0; i < _n1; i++)
    00640     {
    00641       CurU->i = i;
    00642       CurU->val = -INFINITY;
    00643       CurU->Next = CurU+1;
    00644       CurU++;
    00645     }
    00646   (--CurU)->Next = NULL;
    00647
    00648   vHead.Next = CurV = Vr;
    00649   for (j=0; j < _n2; j++)
    00650     {
    00651       CurV->i = j;
    00652       CurV->val = -INFINITY;
    00653       CurV->Next = CurV+1;
    00654       CurV++;
    00655     }
    00656   (--CurV)->Next = NULL;
    00657
    00658   /* FIND THE MAXIMUM ROW AND COLUMN VALUES (Ur[i] AND Vr[j]) */
    00659   for(i=0; i < _n1 ; i++)
    00660     for(j=0; j < _n2 ; j++)
    00661       {
    00662         float v;
    00663         v = _C[i][j];
    00664         if (Ur[i].val <= v)
    00665           Ur[i].val = v;
    00666         if (Vr[j].val <= v)
    00667           Vr[j].val = v;
    00668       }
    00669
    00670   /* COMPUTE THE Delta MATRIX */
    00671   for(i=0; i < _n1 ; i++)
    00672     for(j=0; j < _n2 ; j++)
    00673       Delta[i][j] = _C[i][j] - Ur[i].val - Vr[j].val;
    00674
    00675   /* FIND THE BASIC VARIABLES */
    00676   do
    00677     {
    00678 #if DEBUG_LEVEL > 3
    00679       printf("Ur=");
    00680       for(CurU = uHead.Next; CurU != NULL; CurU = CurU->Next)
    00681         printf("[%d]",CurU-Ur);
    00682       printf("\n");
    00683       printf("Vr=");
    00684       for(CurV = vHead.Next; CurV != NULL; CurV = CurV->Next)
    00685         printf("[%d]",CurV-Vr);
    00686       printf("\n");
    00687       printf("\n\n");
    00688 #endif
    00689
    00690       /* FIND THE SMALLEST Delta[i][j] */
    00691       found = 0;
    00692       deltaMin = INFINITY;
    00693       PrevU = &uHead;
    00694       for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
    00695         {
    00696           int i;
    00697           i = CurU->i;
    00698           PrevV = &vHead;
    00699           for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
    00700             {
    00701               int j;
    00702               j = CurV->i;
    00703               if (deltaMin > Delta[i][j])
    00704                 {
    00705                   deltaMin = Delta[i][j];
    00706                   minI = i;
    00707                   minJ = j;
    00708                   PrevUMinI = PrevU;
    00709                   PrevVMinJ = PrevV;
    00710                   found = 1;
    00711                 }
    00712               PrevV = CurV;
    00713             }
    00714           PrevU = CurU;
    00715         }
    00716
    00717       if (! found)
    00718         break;
    00719
    00720       /* ADD X[minI][minJ] TO THE BASIS, AND ADJUST SUPPLIES AND COST */
    00721       Remember = PrevUMinI->Next;
    00722       addBasicVariable(minI, minJ, S, D, PrevUMinI, PrevVMinJ, &uHead);
    00723
    00724       /* UPDATE THE NECESSARY Delta[][] */
    00725       if (Remember == PrevUMinI->Next)  /* LINE minI WAS DELETED */
    00726         {
    00727           for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
    00728             {
    00729               int j;
    00730               j = CurV->i;
    00731               if (CurV->val == _C[minI][j])  /* COLUMN j NEEDS UPDATING */
    00732                 {
    00733                   /* FIND THE NEW MAXIMUM VALUE IN THE COLUMN */
    00734                   oldVal = CurV->val;
    00735                   CurV->val = -INFINITY;
    00736                   for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
    00737                     {
    00738                       int i;
    00739                       i = CurU->i;
    00740                       if (CurV->val <= _C[i][j])
    00741                         CurV->val = _C[i][j];
    00742                     }
    00743
    00744                   /* IF NEEDED, ADJUST THE RELEVANT Delta[*][j] */
    00745                   diff = oldVal - CurV->val;
    00746                   if (fabs(diff) < EPSILON * _maxC)
    00747                     for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
    00748                       Delta[CurU->i][j] += diff;
    00749                 }
    00750             }
    00751         }
    00752       else  /* COLUMN minJ WAS DELETED */
    00753         {
    00754           for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
    00755             {
    00756               int i;
    00757               i = CurU->i;
    00758               if (CurU->val == _C[i][minJ])  /* ROW i NEEDS UPDATING */
    00759                 {
    00760                   /* FIND THE NEW MAXIMUM VALUE IN THE ROW */
    00761                   oldVal = CurU->val;
    00762                   CurU->val = -INFINITY;
    00763                   for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
    00764                     {
    00765                       int j;
    00766                       j = CurV->i;
    00767                       if(CurU->val <= _C[i][j])
    00768                         CurU->val = _C[i][j];
    00769                     }
    00770
    00771                   /* If NEEDED, ADJUST THE RELEVANT Delta[i][*] */
    00772                   diff = oldVal - CurU->val;
    00773                   if (fabs(diff) < EPSILON * _maxC)
    00774                     for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
    00775                       Delta[i][CurV->i] += diff;
    00776                 }
    00777             }
    00778         }
    00779     } while (uHead.Next != NULL || vHead.Next != NULL);
    00780 }
    00781
    00782
    00783
    00784
    00785 /**********************
    00786     addBasicVariable
    00787 **********************/
    00788 static void addBasicVariable(int minI, int minJ, double *S, double *D,
    00789                              node1_t *PrevUMinI, node1_t *PrevVMinJ,
    00790                              node1_t *UHead)
    00791 {
    00792   double T;
    00793
    00794   if (fabs(S[minI]-D[minJ]) <= EPSILON * _maxW)  /* DEGENERATE CASE */
    00795     {
    00796       T = S[minI];
    00797       S[minI] = 0;
    00798       D[minJ] -= T;
    00799     }
    00800   else if (S[minI] < D[minJ])  /* SUPPLY EXHAUSTED */
    00801     {
    00802       T = S[minI];
    00803       S[minI] = 0;
    00804       D[minJ] -= T;
    00805     }
    00806   else  /* DEMAND EXHAUSTED */
    00807     {
    00808       T = D[minJ];
    00809       D[minJ] = 0;
    00810       S[minI] -= T;
    00811     }
    00812
    00813   /* X(minI,minJ) IS A BASIC VARIABLE */
    00814   _IsX[minI][minJ] = 1;
    00815
    00816   _EndX->val = T;
    00817   _EndX->i = minI;
    00818   _EndX->j = minJ;
    00819   _EndX->NextC = _RowsX[minI];
    00820   _EndX->NextR = _ColsX[minJ];
    00821   _RowsX[minI] = _EndX;
    00822   _ColsX[minJ] = _EndX;
    00823   _EndX++;
    00824
    00825   /* DELETE SUPPLY ROW ONLY IF THE EMPTY, AND IF NOT LAST ROW */
    00826   if (S[minI] == 0 && UHead->Next->Next != NULL)
    00827     PrevUMinI->Next = PrevUMinI->Next->Next;  /* REMOVE ROW FROM LIST */
    00828   else
    00829     PrevVMinJ->Next = PrevVMinJ->Next->Next;  /* REMOVE COLUMN FROM LIST */
    00830 }
    00831
    00832
    00833
    00834
    00835
    00836 /**********************
    00837     printSolution
    00838 **********************/
    00839 static void printSolution()
    00840 {
    00841   node2_t *P;
    00842   double totalCost;
    00843
    00844   totalCost = 0;
    00845
    00846 #if DEBUG_LEVEL > 2
    00847   printf("SIG1\tSIG2\tFLOW\tCOST\n");
    00848 #endif
    00849   for(P=_X; P < _EndX; P++)
    00850     if (P != _EnterX && _IsX[P->i][P->j])
    00851       {
    00852 #if DEBUG_LEVEL > 2
    00853         printf("%d\t%d\t%f\t%f\n", P->i, P->j, P->val, _C[P->i][P->j]);
    00854 #endif
    00855         totalCost += (double)P->val * _C[P->i][P->j];
    00856 //      printf("MY modifications\n");
    00857 //      printf("i = %d\tj = %d\t Flow :%f\n",P->i, P->j,(double)P->val * _C[P->i][P->j]);
    00858       }
    00859
    00860   printf("COST = %f\n", totalCost);
    00861 }
    00862
    00863