[1172]  1  /*


 2  emd.c


 3 


 4  Last update: 3/14/98 (but see emd.h for a list newer changes)


 5 


 6  An implementation of the Earth Movers Distance.


 7  Based of the solution for the Transportation problem as described in


 8  "Introduction to Mathematical Programming" by F. S. Hillier and


 9  G. J. Lieberman, McGrawHill, 1990.


 10 


 11  Copyright (C) 1998 Yossi Rubner


 12  Computer Science Department, Stanford University


 13  EMail: rubner@cs.stanford.edu URL: http://robotics.stanford.edu/~rubner/emd/default.htm


 14  */


 15 


 16  // For a list of changes since 2020, see emd.h


 17 


 18  #include <stdio.h>


 19  #include <stdlib.h>


 20  #include <math.h>


 21  #include <algorithm>


 22 


 23  #include "emd.h"


 24 


 25  #define DEBUG_LEVEL 0


 26  /*


 27  DEBUG_LEVEL:


 28  0 = NO MESSAGES


 29  1 = PRINT THE NUMBER OF ITERATIONS AND THE FINAL RESULT


 30  2 = PRINT THE RESULT AFTER EVERY ITERATION


 31  3 = PRINT ALSO THE FLOW AFTER EVERY ITERATION


 32  4 = PRINT A LOT OF INFORMATION (PROBABLY USEFUL ONLY FOR THE AUTHOR)


 33  */


 34 


 35  /* NEW TYPES DEFINITION */


 36 


 37  /* node1_t IS USED FOR SINGLELINKED LISTS */


 38  typedef struct node1_t {


 39  int i;


 40  double val;


 41  struct node1_t *Next;


 42  } node1_t;


 43 


 44  /* node1_t IS USED FOR DOUBLELINKED LISTS */


 45  typedef struct node2_t {


 46  int i, j;


 47  double val;


 48  struct node2_t *NextC; /* NEXT COLUMN */


 49  struct node2_t *NextR; /* NEXT ROW */


 50  } node2_t;


 51 


 52 


 53 


 54  /* GLOBAL VARIABLE DECLARATION */


 55 


 56  /* VARIABLES TO HANDLE _X EFFICIENTLY */


 57  static node2_t *_EndX, *_EnterX;


 58  static double _maxW;


 59  static float _maxC;


 60 


 61  /* DECLARATION OF FUNCTIONS */


 62  static float init(signature_t *Signature1, signature_t *Signature2,


 63  float (*Dist)(feature_t *, feature_t *), int _n1, int _n2,


 64  float **_CM, node2_t *_XV, char **_IsX, node2_t **_RowsX, node2_t **_ColsX);


 65  static void findBasicVariables(node1_t *U, node1_t *V, int _n1, int _n2, float **_CM, char **_IsX);


 66  static int isOptimal(node1_t *U, node1_t *V, int _n1, int _n2, float **_CM, char **_IsX);


 67  static int findLoop(node2_t **Loop, int _n1, int _n2, node2_t *_XV, node2_t **_RowsX, node2_t **_ColsX);


 68  static void newSol(int _n1, int _n2, node2_t *_XV, char **_IsX, node2_t **_RowsX, node2_t **_ColsX);


 69  static void russel(double *S, double *D, int _n1, int _n2, float **_CM, char **_IsX, node2_t **_RowsX, node2_t **_ColsX);


 70  static void addBasicVariable(int minI, int minJ, double *S, double *D,


 71  node1_t *PrevUMinI, node1_t *PrevVMinJ,


 72  node1_t *UHead, char **_IsX, node2_t **_RowsX, node2_t **_ColsX);


 73  #if DEBUG_LEVEL > 0


 74  static void printSolution();


 75  #endif


 76 


 77 


 78  /******************************************************************************


 79  float emd(signature_t *Signature1, signature_t *Signature2,


 80  float (*Dist)(feature_t *, feature_t *), flow_t *Flow, int *FlowSize)


 81 


 82  where


 83 


 84  Signature1, Signature2 Pointers to signatures that their distance we want


 85  to compute.


 86  Dist Pointer to the ground distance. i.e. the function that computes


 87  the distance between two features.


 88  Flow (Optional) Pointer to a vector of flow_t (defined in emd.h)


 89  where the resulting flow will be stored. Flow must have n1+n21


 90  elements, where n1 and n2 are the sizes of the two signatures


 91  respectively.


 92  If NULL, the flow is not returned.


 93  FlowSize (Optional) Pointer to an integer where the number of elements in


 94  Flow will be stored


 95 


 96  ******************************************************************************/


 97 


 98  float emd(signature_t *Signature1, signature_t *Signature2,


 99  float (*Dist)(feature_t *, feature_t *),


 100  flow_t *Flow, int *FlowSize, int _n1, int _n2)


 101  {


 102  int itr;


 103  int max_n = std::max(_n1, _n2); //max_n was introduced in r1062 instead of the #defined constant MAX_SIG_SIZE1=1000 in the original implementation. max_n is better than the constant, but it would be even better to use either _n1 or _n2, if we only knew what size each individual array should precisely have.


 104  double totalCost;


 105  float w;


 106  node2_t *XP;


 107  flow_t *FlowP;


 108  node1_t *U=new node1_t[max_n], *V=new node1_t[max_n];


 109  /* THE COST MATRIX */


 110  float** _CM = new float*[_n1];


 111  char** _IsX = new char*[_n1];


 112  for(int k = 0; k < _n1; ++k)


 113  {


 114  _CM[k] = new float[_n2];


 115  _IsX[k] = new char[_n2];


 116  }


 117  /* THE BASIC VARIABLES VECTOR */


 118  node2_t *_XV = new node2_t[max_n*2];


 119  node2_t **_RowsX = new node2_t*[max_n*2];


 120  node2_t **_ColsX = new node2_t*[max_n*2];


 121 


 122  w = init(Signature1, Signature2, Dist, _n1, _n2, _CM, _XV, _IsX, _RowsX, _ColsX);


 123 


 124  #if DEBUG_LEVEL > 1


 125  printf("\nINITIAL SOLUTION:\n");


 126  printSolution();


 127  #endif


 128 


 129  if (_n1 > 1 && _n2 > 1) /* IF _n1 = 1 OR _n2 = 1 THEN WE ARE DONE */


 130  {


 131  for (itr = 1; itr < MAX_ITERATIONS; itr++)


 132  {


 133  /* FIND BASIC VARIABLES */


 134  findBasicVariables(U, V, _n1, _n2, _CM, _IsX);


 135 


 136  /* CHECK FOR OPTIMALITY */


 137  if (isOptimal(U, V, _n1, _n2, _CM, _IsX))


 138  break;


 139 


 140  /* IMPROVE SOLUTION */


 141  newSol(_n1, _n2, _XV, _IsX, _RowsX, _ColsX);


 142 


 143  #if DEBUG_LEVEL > 1


 144  printf("\nITERATION # %d \n", itr);


 145  printSolution();


 146  #endif


 147  }


 148 


 149  if (itr == MAX_ITERATIONS) //print a warning message


 150  {


 151  #define EMD_LIMIT_WARNING_MESSAGES 10 //if thousands of such messages are printed, we group them to avoid spamming stderr. Undefine (comment out) to disable.


 152  #ifdef EMD_LIMIT_WARNING_MESSAGES


 153  static int warncounter = 0;


 154  static int warnnextlimit = 1;


 155  warncounter++;


 156  if (warncounter == warnnextlimit)


 157  {


 158  fprintf(stderr, "[repeated %dx in total, next warning at %dx] ", warncounter, warnnextlimit * EMD_LIMIT_WARNING_MESSAGES);


 159  #endif


 160  fprintf(stderr, "emd: Maximum number of iterations (%d) has been reached\n", MAX_ITERATIONS);


 161  #ifdef EMD_LIMIT_WARNING_MESSAGES


 162  warnnextlimit *= EMD_LIMIT_WARNING_MESSAGES;


 163  }


 164  #endif


 165  }


 166  }


 167 


 168  /* COMPUTE THE TOTAL FLOW */


 169  totalCost = 0;


 170  if (Flow != NULL)


 171  FlowP = Flow;


 172  for(XP=_XV; XP < _EndX; XP++)


 173  {


 174  if (XP == _EnterX) /* _EnterX IS THE EMPTY SLOT */


 175  continue;


 176  if (XP>i == Signature1>n  XP>j == Signature2>n) /* DUMMY FEATURE */


 177  continue;


 178 


 179  if (XP>val == 0) /* ZERO FLOW */


 180  continue;


 181 


 182  totalCost += (double)XP>val * _CM[XP>i][XP>j];


 183  if (Flow != NULL)


 184  {


 185  FlowP>from = XP>i;


 186  FlowP>to = XP>j;


 187  FlowP>amount = XP>val;


 188  FlowP++;


 189  }


 190  }


 191  if (Flow != NULL)


 192  *FlowSize = FlowPFlow;


 193 


 194  #if DEBUG_LEVEL > 0


 195  printf("\n*** OPTIMAL SOLUTION (%d ITERATIONS): %f ***\n", itr, totalCost);


 196  #endif


 197 


 198  for(int k = 0; k < _n1; ++k)


 199  {


 200  delete[] _CM[k];


 201  delete[] _IsX[k];


 202  }


 203  delete[] _CM;


 204  delete[] _IsX;


 205  delete[] _XV;


 206  delete[] _RowsX;


 207  delete[] _ColsX;


 208 


 209  delete[] U;


 210  delete[] V;


 211 


 212  /* RETURN THE NORMALIZED COST == EMD */


 213  return (float)(totalCost / w);


 214  }


 215 


 216 


 217 


 218 


 219 


 220  /**********************


 221  init


 222  **********************/


 223  static float init(signature_t *Signature1, signature_t *Signature2,


 224  float (*Dist)(feature_t *, feature_t *), int _n1, int _n2,


 225  float **_CM, node2_t *_XV, char **_IsX, node2_t **_RowsX, node2_t **_ColsX)


 226  {


 227  int i, j;


 228  int max_n = std::max(_n1, _n2); //max_n was introduced in r1062 instead of the #defined constant MAX_SIG_SIZE1=1000 in the original implementation. max_n is better than the constant, but it would be even better to use either _n1 or _n2, if we only knew what size each individual array should precisely have.


 229  double sSum, dSum, diff;


 230  feature_t *P1, *P2;


 231  double *S=new double[max_n], *D=new double[max_n];


 232 


 233 


 234  /* COMPUTE THE DISTANCE MATRIX */


 235  _maxC = 0;


 236  for(i=0, P1=Signature1>Features; i < _n1; i++, P1++)


 237  for(j=0, P2=Signature2>Features; j < _n2; j++, P2++)


 238  {


 239  _CM[i][j] = Dist(P1, P2);


 240  if (_CM[i][j] > _maxC)


 241  _maxC = _CM[i][j];


 242  }


 243 


 244  /* SUM UP THE SUPPLY AND DEMAND */


 245  sSum = 0.0;


 246  for(i=0; i < _n1; i++)


 247  {


 248  S[i] = Signature1>Weights[i];


 249  sSum += Signature1>Weights[i];


 250  _RowsX[i] = NULL;


 251  }


 252  dSum = 0.0;


 253  for(j=0; j < _n2; j++)


 254  {


 255  D[j] = Signature2>Weights[j];


 256  dSum += Signature2>Weights[j];


 257  _ColsX[j] = NULL;


 258  }


 259 


 260  /* IF SUPPLY DIFFERENT THAN THE DEMAND, ADD A ZEROCOST DUMMY CLUSTER */


 261  diff = sSum  dSum;


 262  if (fabs(diff) >= EPSILON * sSum)


 263  {


 264  if (diff < 0.0)


 265  {


 266  for (j=0; j < _n2; j++)


 267  _CM[_n1][j] = 0;


 268  S[_n1] = diff;


 269  _RowsX[_n1] = NULL;


 270  _n1++;


 271  }


 272  else


 273  {


 274  for (i=0; i < _n1; i++)


 275  _CM[i][_n2] = 0;


 276  D[_n2] = diff;


 277  _ColsX[_n2] = NULL;


 278  _n2++;


 279  }


 280  }


 281 


 282  /* INITIALIZE THE BASIC VARIABLE STRUCTURES */


 283  for (i=0; i < _n1; i++)


 284  for (j=0; j < _n2; j++)


 285  _IsX[i][j] = 0;


 286  _EndX = _XV;


 287 


 288  _maxW = sSum > dSum ? sSum : dSum;


 289 


 290  /* FIND INITIAL SOLUTION */


 291  russel(S, D, _n1, _n2, _CM, _IsX, _RowsX, _ColsX);


 292 


 293  _EnterX = _EndX++; /* AN EMPTY SLOT (ONLY _n1+_n21 BASIC VARIABLES) */


 294 


 295  delete[] S;


 296  delete[] D;


 297 


 298  return sSum > dSum ? dSum : sSum;


 299  }


 300 


 301 


 302  /**********************


 303  findBasicVariables


 304  **********************/


 305  static void findBasicVariables(node1_t *U, node1_t *V, int _n1, int _n2, float **_CM, char **_IsX)


 306  {


 307  int i, j, found;


 308  int UfoundNum, VfoundNum;


 309  node1_t u0Head, u1Head, *CurU, *PrevU;


 310  node1_t v0Head, v1Head, *CurV, *PrevV;


 311 


 312  /* INITIALIZE THE ROWS LIST (U) AND THE COLUMNS LIST (V) */


 313  u0Head.Next = CurU = U;


 314  for (i=0; i < _n1; i++)


 315  {


 316  CurU>i = i;


 317  CurU>Next = CurU+1;


 318  CurU++;


 319  }


 320  (CurU)>Next = NULL;


 321  u1Head.Next = NULL;


 322 


 323  CurV = V+1;


 324  v0Head.Next = _n2 > 1 ? V+1 : NULL;


 325  for (j=1; j < _n2; j++)


 326  {


 327  CurV>i = j;


 328  CurV>Next = CurV+1;


 329  CurV++;


 330  }


 331  (CurV)>Next = NULL;


 332  v1Head.Next = NULL;


 333 


 334  /* THERE ARE _n1+_n2 VARIABLES BUT ONLY _n1+_n21 INDEPENDENT EQUATIONS,


 335  SO SET V[0]=0 */


 336  V[0].i = 0;


 337  V[0].val = 0;


 338  v1Head.Next = V;


 339  v1Head.Next>Next = NULL;


 340 


 341  /* LOOP UNTIL ALL VARIABLES ARE FOUND */


 342  UfoundNum=VfoundNum=0;


 343  while (UfoundNum < _n1  VfoundNum < _n2)


 344  {


 345 


 346  #if DEBUG_LEVEL > 3


 347  printf("UfoundNum=%d/%d,VfoundNum=%d/%d\n",UfoundNum,_n1,VfoundNum,_n2);


 348  printf("U0=");


 349  for(CurU = u0Head.Next; CurU != NULL; CurU = CurU>Next)


 350  printf("[%d]",CurUU);


 351  printf("\n");


 352  printf("U1=");


 353  for(CurU = u1Head.Next; CurU != NULL; CurU = CurU>Next)


 354  printf("[%d]",CurUU);


 355  printf("\n");


 356  printf("V0=");


 357  for(CurV = v0Head.Next; CurV != NULL; CurV = CurV>Next)


 358  printf("[%d]",CurVV);


 359  printf("\n");


 360  printf("V1=");


 361  for(CurV = v1Head.Next; CurV != NULL; CurV = CurV>Next)


 362  printf("[%d]",CurVV);


 363  printf("\n\n");


 364  #endif


 365 


 366  found = 0;


 367  if (VfoundNum < _n2)


 368  {


 369  /* LOOP OVER ALL MARKED COLUMNS */


 370  PrevV = &v1Head;


 371  for (CurV=v1Head.Next; CurV != NULL; CurV=CurV>Next)


 372  {


 373  j = CurV>i;


 374  /* FIND THE VARIABLES IN COLUMN j */


 375  PrevU = &u0Head;


 376  for (CurU=u0Head.Next; CurU != NULL; CurU=CurU>Next)


 377  {


 378  i = CurU>i;


 379  if (_IsX[i][j])


 380  {


 381  /* COMPUTE U[i] */


 382  CurU>val = _CM[i][j]  CurV>val;


 383  /* ...AND ADD IT TO THE MARKED LIST */


 384  PrevU>Next = CurU>Next;


 385  CurU>Next = u1Head.Next != NULL ? u1Head.Next : NULL;


 386  u1Head.Next = CurU;


 387  CurU = PrevU;


 388  }


 389  else


 390  PrevU = CurU;


 391  }


 392  PrevV>Next = CurV>Next;


 393  VfoundNum++;


 394  found = 1;


 395  }


 396  }


 397  if (UfoundNum < _n1)


 398  {


 399  /* LOOP OVER ALL MARKED ROWS */


 400  PrevU = &u1Head;


 401  for (CurU=u1Head.Next; CurU != NULL; CurU=CurU>Next)


 402  {


 403  i = CurU>i;


 404  /* FIND THE VARIABLES IN ROWS i */


 405  PrevV = &v0Head;


 406  for (CurV=v0Head.Next; CurV != NULL; CurV=CurV>Next)


 407  {


 408  j = CurV>i;


 409  if (_IsX[i][j])


 410  {


 411  /* COMPUTE V[j] */


 412  CurV>val = _CM[i][j]  CurU>val;


 413  /* ...AND ADD IT TO THE MARKED LIST */


 414  PrevV>Next = CurV>Next;


 415  CurV>Next = v1Head.Next != NULL ? v1Head.Next: NULL;


 416  v1Head.Next = CurV;


 417  CurV = PrevV;


 418  }


 419  else


 420  PrevV = CurV;


 421  }


 422  PrevU>Next = CurU>Next;


 423  UfoundNum++;


 424  found = 1;


 425  }


 426  }


 427  if (! found)


 428  {


 429  fprintf(stderr, "emd: Unexpected error in findBasicVariables!\n");


 430  fprintf(stderr, "This typically happens when the EPSILON defined in\n");


 431  fprintf(stderr, "emd.h is not right for the scale of the problem.\n");


 432  exit(1);


 433  }


 434  }


 435  }


 436 


 437 


 438 


 439 


 440  /**********************


 441  isOptimal


 442  **********************/


 443  static int isOptimal(node1_t *U, node1_t *V, int _n1, int _n2, float **_CM, char **_IsX)


 444  {


 445  double delta, deltaMin;


 446  int i, j, minI, minJ;


 447 


 448  /* FIND THE MINIMAL CijUiVj OVER ALL i,j */


 449  deltaMin = INFINITY;


 450  for(i=0; i < _n1; i++)


 451  for(j=0; j < _n2; j++)


 452  if (! _IsX[i][j])


 453  {


 454  delta = _CM[i][j]  U[i].val  V[j].val;


 455  if (deltaMin > delta)


 456  {


 457  deltaMin = delta;


 458  minI = i;


 459  minJ = j;


 460  }


 461  }


 462 


 463  #if DEBUG_LEVEL > 3


 464  printf("deltaMin=%f\n", deltaMin);


 465  #endif


 466 


 467  if (deltaMin == INFINITY)


 468  {


 469  fprintf(stderr, "emd: Unexpected error in isOptimal.\n");


 470  exit(0);


 471  }


 472 


 473  _EnterX>i = minI;


 474  _EnterX>j = minJ;


 475 


 476  /* IF NO NEGATIVE deltaMin, WE FOUND THE OPTIMAL SOLUTION */


 477  return deltaMin >= EPSILON * _maxC;


 478 


 479  /*


 480  return deltaMin >= EPSILON;


 481  */


 482  }


 483 


 484 


 485  /**********************


 486  newSol


 487  **********************/


 488  static void newSol(int _n1, int _n2, node2_t * _XV, char **_IsX, node2_t **_RowsX, node2_t **_ColsX)


 489  {


 490  int i, j, k;


 491  int max_n = std::max(_n1, _n2); //max_n was introduced in r1062 instead of the #defined constant MAX_SIG_SIZE1=1000 in the original implementation. max_n is better than the constant, but it would be even better to use either _n1 or _n2, if we only knew what size each individual array should precisely have.


 492  double xMin;


 493  int steps;


 494  node2_t **Loop=new node2_t*[2*max_n], *CurX, *LeaveX;


 495 


 496  #if DEBUG_LEVEL > 3


 497  printf("EnterX = (%d,%d)\n", _EnterX>i, _EnterX>j);


 498  #endif


 499 


 500  /* ENTER THE NEW BASIC VARIABLE */


 501  i = _EnterX>i;


 502  j = _EnterX>j;


 503  _IsX[i][j] = 1;


 504  _EnterX>NextC = _RowsX[i];


 505  _EnterX>NextR = _ColsX[j];


 506  _EnterX>val = 0;


 507  _RowsX[i] = _EnterX;


 508  _ColsX[j] = _EnterX;


 509 


 510  /* FIND A CHAIN REACTION */


 511  steps = findLoop(Loop, _n1, _n2, _XV, _RowsX, _ColsX);


 512 


 513  /* FIND THE LARGEST VALUE IN THE LOOP */


 514  xMin = INFINITY;


 515  for (k=1; k < steps; k+=2)


 516  {


 517  if (Loop[k]>val < xMin)


 518  {


 519  LeaveX = Loop[k];


 520  xMin = Loop[k]>val;


 521  }


 522  }


 523 


 524  /* UPDATE THE LOOP */


 525  for (k=0; k < steps; k+=2)


 526  {


 527  Loop[k]>val += xMin;


 528  Loop[k+1]>val = xMin;


 529  }


 530 


 531  #if DEBUG_LEVEL > 3


 532  printf("LeaveX = (%d,%d)\n", LeaveX>i, LeaveX>j);


 533  #endif


 534 


 535  /* REMOVE THE LEAVING BASIC VARIABLE */


 536  i = LeaveX>i;


 537  j = LeaveX>j;


 538  _IsX[i][j] = 0;


 539  if (_RowsX[i] == LeaveX)


 540  _RowsX[i] = LeaveX>NextC;


 541  else


 542  for (CurX=_RowsX[i]; CurX != NULL; CurX = CurX>NextC)


 543  if (CurX>NextC == LeaveX)


 544  {


 545  CurX>NextC = CurX>NextC>NextC;


 546  break;


 547  }


 548  if (_ColsX[j] == LeaveX)


 549  _ColsX[j] = LeaveX>NextR;


 550  else


 551  for (CurX=_ColsX[j]; CurX != NULL; CurX = CurX>NextR)


 552  if (CurX>NextR == LeaveX)


 553  {


 554  CurX>NextR = CurX>NextR>NextR;


 555  break;


 556  }


 557 


 558  /* SET _EnterX TO BE THE NEW EMPTY SLOT */


 559  _EnterX = LeaveX;


 560 


 561  delete[] Loop;


 562  }


 563 


 564 


 565 


 566  /**********************


 567  findLoop


 568  **********************/


 569  static int findLoop(node2_t **Loop, int _n1, int _n2, node2_t *_XV, node2_t **_RowsX, node2_t **_ColsX)


 570  {


 571  int i, steps;


 572  int max_n = std::max(_n1, _n2); //max_n was introduced in r1062 instead of the #defined constant MAX_SIG_SIZE1=1000 in the original implementation. max_n is better than the constant, but it would be even better to use either _n1 or _n2, if we only knew what size each individual array should precisely have.


 573  node2_t **CurX, *NewX;


 574  char *IsUsed=new char[2*max_n];


 575 


 576  for (i=0; i < _n1+_n2; i++)


 577  IsUsed[i] = 0;


 578 


 579  CurX = Loop;


 580  NewX = *CurX = _EnterX;


 581  IsUsed[_EnterX_XV] = 1;


 582  steps = 1;


 583 


 584  do


 585  {


 586  if (steps%2 == 1)


 587  {


 588  /* FIND AN UNUSED X IN THE ROW */


 589  NewX = _RowsX[NewX>i];


 590  while (NewX != NULL && IsUsed[NewX_XV])


 591  NewX = NewX>NextC;


 592  }


 593  else


 594  {


 595  /* FIND AN UNUSED X IN THE COLUMN, OR THE ENTERING X */


 596  NewX = _ColsX[NewX>j];


 597  while (NewX != NULL && IsUsed[NewX_XV] && NewX != _EnterX)


 598  NewX = NewX>NextR;


 599  if (NewX == _EnterX)


 600  break;


 601  }


 602 


 603  if (NewX != NULL) /* FOUND THE NEXT X */


 604  {


 605  /* ADD X TO THE LOOP */


 606  *++CurX = NewX;


 607  IsUsed[NewX_XV] = 1;


 608  steps++;


 609  #if DEBUG_LEVEL > 3


 610  printf("steps=%d, NewX=(%d,%d)\n", steps, NewX>i, NewX>j);


 611  #endif


 612  }


 613  else /* DIDN'T FIND THE NEXT X */


 614  {


 615  /* BACKTRACK */


 616  do


 617  {


 618  NewX = *CurX;


 619  do


 620  {


 621  if (steps%2 == 1)


 622  NewX = NewX>NextR;


 623  else


 624  NewX = NewX>NextC;


 625  } while (NewX != NULL && IsUsed[NewX_XV]);


 626 


 627  if (NewX == NULL)


 628  {


 629  IsUsed[*CurX_XV] = 0;


 630  CurX;


 631  steps;


 632  }


 633  } while (NewX == NULL && CurX >= Loop);


 634 


 635  #if DEBUG_LEVEL > 3


 636  printf("BACKTRACKING TO: steps=%d, NewX=(%d,%d)\n",


 637  steps, NewX>i, NewX>j);


 638  #endif


 639  IsUsed[*CurX_XV] = 0;


 640  *CurX = NewX;


 641  IsUsed[NewX_XV] = 1;


 642  }


 643  } while(CurX >= Loop);


 644 


 645  if (CurX == Loop)


 646  {


 647  fprintf(stderr, "emd: Unexpected error in findLoop!\n");


 648  exit(1);


 649  }


 650  #if DEBUG_LEVEL > 3


 651  printf("FOUND LOOP:\n");


 652  for (i=0; i < steps; i++)


 653  printf("%d: (%d,%d)\n", i, Loop[i]>i, Loop[i]>j);


 654  #endif


 655 


 656  delete[] IsUsed;


 657 


 658  return steps;


 659  }


 660 


 661 


 662 


 663  /**********************


 664  russel


 665  **********************/


 666  static void russel(double *S, double *D, int _n1, int _n2, float **_CM, char **_IsX, node2_t **_RowsX, node2_t **_ColsX)


 667  {


 668  int i, j, found, minI, minJ;


 669  int max_n = std::max(_n1, _n2); //max_n was introduced in r1062 instead of the #defined constant MAX_SIG_SIZE1=1000 in the original implementation. max_n is better than the constant, but it would be even better to use either _n1 or _n2, if we only knew what size each individual array should precisely have.


 670  double deltaMin, oldVal, diff;


 671  double** Delta = new double*[_n1];


 672  for(int k = 0; k < _n1; ++k)


 673  Delta[k] = new double[_n2];


 674  node1_t *Ur=new node1_t[max_n], *Vr=new node1_t[max_n];


 675  node1_t uHead, *CurU, *PrevU;


 676  node1_t vHead, *CurV, *PrevV;


 677  node1_t *PrevUMinI, *PrevVMinJ, *Remember;


 678 


 679  /* INITIALIZE THE ROWS LIST (Ur), AND THE COLUMNS LIST (Vr) */


 680  uHead.Next = CurU = Ur;


 681  for (i=0; i < _n1; i++)


 682  {


 683  CurU>i = i;


 684  CurU>val = INFINITY;


 685  CurU>Next = CurU+1;


 686  CurU++;


 687  }


 688  (CurU)>Next = NULL;


 689 


 690  vHead.Next = CurV = Vr;


 691  for (j=0; j < _n2; j++)


 692  {


 693  CurV>i = j;


 694  CurV>val = INFINITY;


 695  CurV>Next = CurV+1;


 696  CurV++;


 697  }


 698  (CurV)>Next = NULL;


 699 


 700  /* FIND THE MAXIMUM ROW AND COLUMN VALUES (Ur[i] AND Vr[j]) */


 701  for(i=0; i < _n1 ; i++)


 702  for(j=0; j < _n2 ; j++)


 703  {


 704  float v;


 705  v = _CM[i][j];


 706  if (Ur[i].val <= v)


 707  Ur[i].val = v;


 708  if (Vr[j].val <= v)


 709  Vr[j].val = v;


 710  }


 711 


 712  /* COMPUTE THE Delta MATRIX */


 713  for(i=0; i < _n1 ; i++)


 714  for(j=0; j < _n2 ; j++)


 715  Delta[i][j] = _CM[i][j]  Ur[i].val  Vr[j].val;


 716 


 717  /* FIND THE BASIC VARIABLES */


 718  do


 719  {


 720  #if DEBUG_LEVEL > 3


 721  printf("Ur=");


 722  for(CurU = uHead.Next; CurU != NULL; CurU = CurU>Next)


 723  printf("[%d]",CurUUr);


 724  printf("\n");


 725  printf("Vr=");


 726  for(CurV = vHead.Next; CurV != NULL; CurV = CurV>Next)


 727  printf("[%d]",CurVVr);


 728  printf("\n");


 729  printf("\n\n");


 730  #endif


 731 


 732  /* FIND THE SMALLEST Delta[i][j] */


 733  found = 0;


 734  deltaMin = INFINITY;


 735  PrevU = &uHead;


 736  for (CurU=uHead.Next; CurU != NULL; CurU=CurU>Next)


 737  {


 738  int i;


 739  i = CurU>i;


 740  PrevV = &vHead;


 741  for (CurV=vHead.Next; CurV != NULL; CurV=CurV>Next)


 742  {


 743  int j;


 744  j = CurV>i;


 745  if (deltaMin > Delta[i][j])


 746  {


 747  deltaMin = Delta[i][j];


 748  minI = i;


 749  minJ = j;


 750  PrevUMinI = PrevU;


 751  PrevVMinJ = PrevV;


 752  found = 1;


 753  }


 754  PrevV = CurV;


 755  }


 756  PrevU = CurU;


 757  }


 758 


 759  if (! found)


 760  break;


 761 


 762  /* ADD X[minI][minJ] TO THE BASIS, AND ADJUST SUPPLIES AND COST */


 763  Remember = PrevUMinI>Next;


 764  addBasicVariable(minI, minJ, S, D, PrevUMinI, PrevVMinJ, &uHead, _IsX, _RowsX, _ColsX);


 765 


 766  /* UPDATE THE NECESSARY Delta[][] */


 767  if (Remember == PrevUMinI>Next) /* LINE minI WAS DELETED */


 768  {


 769  for (CurV=vHead.Next; CurV != NULL; CurV=CurV>Next)


 770  {


 771  int j;


 772  j = CurV>i;


 773  if (CurV>val == _CM[minI][j]) /* COLUMN j NEEDS UPDATING */


 774  {


 775  /* FIND THE NEW MAXIMUM VALUE IN THE COLUMN */


 776  oldVal = CurV>val;


 777  CurV>val = INFINITY;


 778  for (CurU=uHead.Next; CurU != NULL; CurU=CurU>Next)


 779  {


 780  int i;


 781  i = CurU>i;


 782  if (CurV>val <= _CM[i][j])


 783  CurV>val = _CM[i][j];


 784  }


 785 


 786  /* IF NEEDED, ADJUST THE RELEVANT Delta[*][j] */


 787  diff = oldVal  CurV>val;


 788  if (fabs(diff) < EPSILON * _maxC)


 789  for (CurU=uHead.Next; CurU != NULL; CurU=CurU>Next)


 790  Delta[CurU>i][j] += diff;


 791  }


 792  }


 793  }


 794  else /* COLUMN minJ WAS DELETED */


 795  {


 796  for (CurU=uHead.Next; CurU != NULL; CurU=CurU>Next)


 797  {


 798  int i;


 799  i = CurU>i;


 800  if (CurU>val == _CM[i][minJ]) /* ROW i NEEDS UPDATING */


 801  {


 802  /* FIND THE NEW MAXIMUM VALUE IN THE ROW */


 803  oldVal = CurU>val;


 804  CurU>val = INFINITY;


 805  for (CurV=vHead.Next; CurV != NULL; CurV=CurV>Next)


 806  {


 807  int j;


 808  j = CurV>i;


 809  if(CurU>val <= _CM[i][j])


 810  CurU>val = _CM[i][j];


 811  }


 812 


 813  /* If NEEDED, ADJUST THE RELEVANT Delta[i][*] */


 814  diff = oldVal  CurU>val;


 815  if (fabs(diff) < EPSILON * _maxC)


 816  for (CurV=vHead.Next; CurV != NULL; CurV=CurV>Next)


 817  Delta[i][CurV>i] += diff;


 818  }


 819  }


 820  }


 821  } while (uHead.Next != NULL  vHead.Next != NULL);


 822 


 823  delete[] Ur;


 824  delete[] Vr;


 825  for(int k = 0; k < _n1; ++k)


 826  delete[] Delta[k];


 827  delete[] Delta;


 828  }


 829 


 830 


 831 


 832 


 833  /**********************


 834  addBasicVariable


 835  **********************/


 836  static void addBasicVariable(int minI, int minJ, double *S, double *D,


 837  node1_t *PrevUMinI, node1_t *PrevVMinJ,


 838  node1_t *UHead, char **_IsX, node2_t **_RowsX, node2_t **_ColsX)


 839  {


 840  double T;


 841 


 842  if (fabs(S[minI]D[minJ]) <= EPSILON * _maxW) /* DEGENERATE CASE */


 843  {


 844  T = S[minI];


 845  S[minI] = 0;


 846  D[minJ] = T;


 847  }


 848  else if (S[minI] < D[minJ]) /* SUPPLY EXHAUSTED */


 849  {


 850  T = S[minI];


 851  S[minI] = 0;


 852  D[minJ] = T;


 853  }


 854  else /* DEMAND EXHAUSTED */


 855  {


 856  T = D[minJ];


 857  D[minJ] = 0;


 858  S[minI] = T;


 859  }


 860 


 861  /* X(minI,minJ) IS A BASIC VARIABLE */


 862  _IsX[minI][minJ] = 1;


 863 


 864  _EndX>val = T;


 865  _EndX>i = minI;


 866  _EndX>j = minJ;


 867  _EndX>NextC = _RowsX[minI];


 868  _EndX>NextR = _ColsX[minJ];


 869  _RowsX[minI] = _EndX;


 870  _ColsX[minJ] = _EndX;


 871  _EndX++;


 872 


 873  /* DELETE SUPPLY ROW ONLY IF THE EMPTY, AND IF NOT LAST ROW */


 874  if (S[minI] == 0 && UHead>Next>Next != NULL)


 875  PrevUMinI>Next = PrevUMinI>Next>Next; /* REMOVE ROW FROM LIST */


 876  else


 877  PrevVMinJ>Next = PrevVMinJ>Next>Next; /* REMOVE COLUMN FROM LIST */


 878  }

