source: cpp/frams/model/similarity/EMD/emd.c @ 1271

/*
    emd.c

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

    An implementation of the Earth Movers Distance.
    Based of the solution for the Transportation problem as described in
    "Introduction to Mathematical Programming" by F. S. Hillier and 
    G. J. Lieberman, McGraw-Hill, 1990.

    Copyright (C) 1998 Yossi Rubner
    Computer Science Department, Stanford University
    E-Mail: rubner@cs.stanford.edu   URL: http://robotics.stanford.edu/~rubner/emd/default.htm
*/

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <algorithm>

#include "emd.h"

#define DEBUG_LEVEL 0
/*
 DEBUG_LEVEL:
   0 = NO MESSAGES
   1 = PRINT THE NUMBER OF ITERATIONS AND THE FINAL RESULT
   2 = PRINT THE RESULT AFTER EVERY ITERATION
   3 = PRINT ALSO THE FLOW AFTER EVERY ITERATION
   4 = PRINT A LOT OF INFORMATION (PROBABLY USEFUL ONLY FOR THE AUTHOR)
*/

/* NEW TYPES DEFINITION */

/* node1_t IS USED FOR SINGLE-LINKED LISTS */
typedef struct node1_t {
  int i;
  double val;
  struct node1_t *Next;
} node1_t;

/* node1_t IS USED FOR DOUBLE-LINKED LISTS */
typedef struct node2_t {
  int i, j;
  double val;
  struct node2_t *NextC;               /* NEXT COLUMN */
  struct node2_t *NextR;               /* NEXT ROW */
} node2_t;



/* GLOBAL VARIABLE DECLARATION */

/* VARIABLES TO HANDLE _X EFFICIENTLY */
static node2_t *_EndX, *_EnterX;
static double _maxW;
static float _maxC;

/* DECLARATION OF FUNCTIONS */
static float init(signature_t *Signature1, signature_t *Signature2,
                  float (*Dist)(feature_t *, feature_t *), int _n1, int _n2,
        float **_CM, node2_t *_XV, char **_IsX, node2_t **_RowsX, node2_t **_ColsX);
static void findBasicVariables(node1_t *U, node1_t *V, int _n1, int _n2, float **_CM, char **_IsX);
static int isOptimal(node1_t *U, node1_t *V, int _n1, int _n2, float **_CM, char **_IsX);
static int findLoop(node2_t **Loop, int _n1, int _n2, node2_t *_XV, node2_t **_RowsX, node2_t **_ColsX);
static void newSol(int _n1, int _n2, node2_t *_XV, char **_IsX, node2_t **_RowsX, node2_t **_ColsX);
static void russel(double *S, double *D, int _n1, int _n2, float **_CM, char **_IsX, node2_t **_RowsX, node2_t **_ColsX);
static void addBasicVariable(int minI, int minJ, double *S, double *D, 
                             node1_t *PrevUMinI, node1_t *PrevVMinJ,
                             node1_t *UHead, char **_IsX, node2_t **_RowsX, node2_t **_ColsX);
#if DEBUG_LEVEL > 0
static void printSolution();
#endif


/******************************************************************************
float emd(signature_t *Signature1, signature_t *Signature2,
          float (*Dist)(feature_t *, feature_t *), flow_t *Flow, int *FlowSize)
  
where

   Signature1, Signature2  Pointers to signatures that their distance we want
              to compute.
   Dist       Pointer to the ground distance. i.e. the function that computes
              the distance between two features.
   Flow       (Optional) Pointer to a vector of flow_t (defined in emd.h) 
              where the resulting flow will be stored. Flow must have n1+n2-1
              elements, where n1 and n2 are the sizes of the two signatures
              respectively.
              If NULL, the flow is not returned.
   FlowSize   (Optional) Pointer to an integer where the number of elements in
              Flow will be stored
              
******************************************************************************/

float emd(signature_t *Signature1, signature_t *Signature2,
          float (*Dist)(feature_t *, feature_t *),
          flow_t *Flow, int *FlowSize, int _n1, int _n2)
{
  int itr;
  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.
  double totalCost;
  float w;
  node2_t *XP;
  flow_t *FlowP;
  node1_t *U=new node1_t[max_n], *V=new node1_t[max_n];
  /* THE COST MATRIX */
  float** _CM = new float*[_n1];
  char** _IsX = new char*[_n1];
  for(int k = 0; k < _n1; ++k)
  {
    _CM[k] = new float[_n2];
    _IsX[k] = new char[_n2];
  }
  /* THE BASIC VARIABLES VECTOR */  
  node2_t *_XV = new node2_t[max_n*2];
  node2_t **_RowsX = new node2_t*[max_n*2];
  node2_t **_ColsX = new node2_t*[max_n*2];
  
  w = init(Signature1, Signature2, Dist, _n1, _n2, _CM, _XV, _IsX, _RowsX, _ColsX);

#if DEBUG_LEVEL > 1
  printf("\nINITIAL SOLUTION:\n");
  printSolution();
#endif
 
  if (_n1 > 1 && _n2 > 1)  /* IF _n1 = 1 OR _n2 = 1 THEN WE ARE DONE */
    {
      for (itr = 1; itr < MAX_ITERATIONS; itr++)
        {
          /* FIND BASIC VARIABLES */
          findBasicVariables(U, V, _n1, _n2, _CM, _IsX);
          
          /* CHECK FOR OPTIMALITY */
          if (isOptimal(U, V, _n1, _n2, _CM, _IsX))
            break;
          
          /* IMPROVE SOLUTION */
          newSol(_n1, _n2, _XV, _IsX, _RowsX, _ColsX);
          
#if DEBUG_LEVEL > 1
          printf("\nITERATION # %d \n", itr);
          printSolution();
#endif
        }

          if (itr == MAX_ITERATIONS) //print a warning message
          {
#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.
#ifdef EMD_LIMIT_WARNING_MESSAGES
                  static int warncounter = 0;
                  static int warnnextlimit = 1;
                  warncounter++;
                  if (warncounter == warnnextlimit)
                  {
                          fprintf(stderr, "[repeated %dx in total, next warning at %dx] ", warncounter, warnnextlimit * EMD_LIMIT_WARNING_MESSAGES);
#endif
                          fprintf(stderr, "emd: Maximum number of iterations (%d) has been reached\n", MAX_ITERATIONS);
#ifdef EMD_LIMIT_WARNING_MESSAGES
                          warnnextlimit *= EMD_LIMIT_WARNING_MESSAGES;
                  }
#endif
          }
  }

  /* COMPUTE THE TOTAL FLOW */
  totalCost = 0;
  if (Flow != NULL)
    FlowP = Flow;
  for(XP=_XV; XP < _EndX; XP++)
    {
      if (XP == _EnterX)  /* _EnterX IS THE EMPTY SLOT */
        continue;
      if (XP->i == Signature1->n || XP->j == Signature2->n)  /* DUMMY FEATURE */
        continue;
      
      if (XP->val == 0)  /* ZERO FLOW */
        continue;

      totalCost += (double)XP->val * _CM[XP->i][XP->j];
      if (Flow != NULL)
        {
          FlowP->from = XP->i;
          FlowP->to = XP->j;
          FlowP->amount = XP->val;
          FlowP++;
        }
    }
  if (Flow != NULL)
    *FlowSize = FlowP-Flow;

#if DEBUG_LEVEL > 0
  printf("\n*** OPTIMAL SOLUTION (%d ITERATIONS): %f ***\n", itr, totalCost);
#endif

  for(int k = 0; k < _n1; ++k)
  {
    delete[] _CM[k];
    delete[] _IsX[k];
  }
  delete[] _CM;
  delete[] _IsX;
  delete[] _XV;
  delete[] _RowsX;
  delete[] _ColsX;  

  delete[] U;
  delete[] V;
  
  /* RETURN THE NORMALIZED COST == EMD */    
  return (float)(totalCost / w);
}





/**********************
   init
**********************/
static float init(signature_t *Signature1, signature_t *Signature2, 
                  float (*Dist)(feature_t *, feature_t *), int _n1, int _n2,
        float **_CM, node2_t *_XV, char **_IsX, node2_t **_RowsX, node2_t **_ColsX)
{
  int i, j;
  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.
  double sSum, dSum, diff;
  feature_t *P1, *P2;
  double *S=new double[max_n], *D=new double[max_n];

  
  /* COMPUTE THE DISTANCE MATRIX */
  _maxC = 0;
  for(i=0, P1=Signature1->Features; i < _n1; i++, P1++)
    for(j=0, P2=Signature2->Features; j < _n2; j++, P2++) 
      {
        _CM[i][j] = Dist(P1, P2);
        if (_CM[i][j] > _maxC)
          _maxC = _CM[i][j];
      }
        
  /* SUM UP THE SUPPLY AND DEMAND */
  sSum = 0.0;
  for(i=0; i < _n1; i++)
    {
      S[i] = Signature1->Weights[i];
      sSum += Signature1->Weights[i];
      _RowsX[i] = NULL;
    }
  dSum = 0.0;
  for(j=0; j < _n2; j++)
    {
      D[j] = Signature2->Weights[j];
      dSum += Signature2->Weights[j];
      _ColsX[j] = NULL;
    }

  /* IF SUPPLY DIFFERENT THAN THE DEMAND, ADD A ZERO-COST DUMMY CLUSTER */
  diff = sSum - dSum;
  if (fabs(diff) >= EPSILON * sSum)
    {
      if (diff < 0.0)
        {
          for (j=0; j < _n2; j++)
            _CM[_n1][j] = 0;
          S[_n1] = -diff;
          _RowsX[_n1] = NULL;
          _n1++;
        }
      else
        {
          for (i=0; i < _n1; i++)
            _CM[i][_n2] = 0;
          D[_n2] = diff;
          _ColsX[_n2] = NULL;
          _n2++;
        }
    }

  /* INITIALIZE THE BASIC VARIABLE STRUCTURES */
  for (i=0; i < _n1; i++)
    for (j=0; j < _n2; j++)
        _IsX[i][j] = 0;
  _EndX = _XV;
   
  _maxW = sSum > dSum ? sSum : dSum;

  /* FIND INITIAL SOLUTION */
  russel(S, D, _n1, _n2, _CM, _IsX, _RowsX, _ColsX);

  _EnterX = _EndX++;  /* AN EMPTY SLOT (ONLY _n1+_n2-1 BASIC VARIABLES) */

  delete[] S;
  delete[] D;

  return sSum > dSum ? dSum : sSum;
}


/**********************
    findBasicVariables
 **********************/
static void findBasicVariables(node1_t *U, node1_t *V, int _n1, int _n2, float **_CM, char **_IsX)
{
  int i, j, found;
  int UfoundNum, VfoundNum;
  node1_t u0Head, u1Head, *CurU, *PrevU;
  node1_t v0Head, v1Head, *CurV, *PrevV;

  /* INITIALIZE THE ROWS LIST (U) AND THE COLUMNS LIST (V) */
  u0Head.Next = CurU = U;
  for (i=0; i < _n1; i++)
    {
      CurU->i = i;
      CurU->Next = CurU+1;
      CurU++;
    }
  (--CurU)->Next = NULL;
  u1Head.Next = NULL;

  CurV = V+1;
  v0Head.Next = _n2 > 1 ? V+1 : NULL;
  for (j=1; j < _n2; j++)
    {
      CurV->i = j;
      CurV->Next = CurV+1;
      CurV++;
    }
  (--CurV)->Next = NULL;
  v1Head.Next = NULL;

  /* THERE ARE _n1+_n2 VARIABLES BUT ONLY _n1+_n2-1 INDEPENDENT EQUATIONS,
     SO SET V[0]=0 */
  V[0].i = 0;
  V[0].val = 0;
  v1Head.Next = V;
  v1Head.Next->Next = NULL;

  /* LOOP UNTIL ALL VARIABLES ARE FOUND */
  UfoundNum=VfoundNum=0;
  while (UfoundNum < _n1 || VfoundNum < _n2)
    {

#if DEBUG_LEVEL > 3
      printf("UfoundNum=%d/%d,VfoundNum=%d/%d\n",UfoundNum,_n1,VfoundNum,_n2);
      printf("U0=");
      for(CurU = u0Head.Next; CurU != NULL; CurU = CurU->Next)
        printf("[%d]",CurU-U);
      printf("\n");
      printf("U1=");
      for(CurU = u1Head.Next; CurU != NULL; CurU = CurU->Next)
        printf("[%d]",CurU-U);
      printf("\n");
      printf("V0=");
      for(CurV = v0Head.Next; CurV != NULL; CurV = CurV->Next)
        printf("[%d]",CurV-V);
      printf("\n");
      printf("V1=");
      for(CurV = v1Head.Next; CurV != NULL; CurV = CurV->Next)
        printf("[%d]",CurV-V);
      printf("\n\n");
#endif
      
      found = 0;
      if (VfoundNum < _n2)
        {
          /* LOOP OVER ALL MARKED COLUMNS */
          PrevV = &v1Head;
          for (CurV=v1Head.Next; CurV != NULL; CurV=CurV->Next)
            {
              j = CurV->i;
              /* FIND THE VARIABLES IN COLUMN j */
              PrevU = &u0Head;
              for (CurU=u0Head.Next; CurU != NULL; CurU=CurU->Next)
                {
                  i = CurU->i;
                  if (_IsX[i][j])
                    {
                      /* COMPUTE U[i] */
                      CurU->val = _CM[i][j] - CurV->val;
                      /* ...AND ADD IT TO THE MARKED LIST */
                      PrevU->Next = CurU->Next;
                      CurU->Next = u1Head.Next != NULL ? u1Head.Next : NULL;
                      u1Head.Next = CurU;
                      CurU = PrevU;
                    }
                  else
                    PrevU = CurU;
                }
              PrevV->Next = CurV->Next;
              VfoundNum++;
              found = 1;
            }
        }
     if (UfoundNum < _n1)
        {
          /* LOOP OVER ALL MARKED ROWS */
          PrevU = &u1Head;
          for (CurU=u1Head.Next; CurU != NULL; CurU=CurU->Next)
            {
              i = CurU->i;
              /* FIND THE VARIABLES IN ROWS i */
              PrevV = &v0Head;
              for (CurV=v0Head.Next; CurV != NULL; CurV=CurV->Next)
                {
                  j = CurV->i;
                  if (_IsX[i][j])
                    {
                      /* COMPUTE V[j] */
                      CurV->val = _CM[i][j] - CurU->val;
                      /* ...AND ADD IT TO THE MARKED LIST */
                      PrevV->Next = CurV->Next;
                      CurV->Next = v1Head.Next != NULL ? v1Head.Next: NULL;
                      v1Head.Next = CurV;
                      CurV = PrevV;
                    }
                  else
                    PrevV = CurV;
                }
              PrevU->Next = CurU->Next;
              UfoundNum++;
              found = 1;
            }
        }
     if (! found)
       {
         fprintf(stderr, "emd: Unexpected error in findBasicVariables!\n");
         fprintf(stderr, "This typically happens when the EPSILON defined in\n");
         fprintf(stderr, "emd.h is not right for the scale of the problem.\n");
         exit(1);
       }
    }
}




/**********************
    isOptimal
 **********************/
static int isOptimal(node1_t *U, node1_t *V, int _n1, int _n2, float **_CM, char **_IsX)
{    
  double delta, deltaMin;
  int i, j, minI, minJ;

  /* FIND THE MINIMAL Cij-Ui-Vj OVER ALL i,j */
  deltaMin = INFINITY;
  for(i=0; i < _n1; i++)
    for(j=0; j < _n2; j++)
      if (! _IsX[i][j])
        {
          delta = _CM[i][j] - U[i].val - V[j].val;
          if (deltaMin > delta)
            {
              deltaMin = delta;
              minI = i;
              minJ = j;
            }
        }

#if DEBUG_LEVEL > 3
  printf("deltaMin=%f\n", deltaMin);
#endif

   if (deltaMin == INFINITY)
     {
       fprintf(stderr, "emd: Unexpected error in isOptimal.\n");
       exit(0);
     }
   
   _EnterX->i = minI;
   _EnterX->j = minJ;
   
   /* IF NO NEGATIVE deltaMin, WE FOUND THE OPTIMAL SOLUTION */
   return deltaMin >= -EPSILON * _maxC;

/*
   return deltaMin >= -EPSILON;
 */
}


/**********************
    newSol
**********************/
static void newSol(int _n1, int _n2, node2_t * _XV, char **_IsX, node2_t **_RowsX, node2_t **_ColsX)
{
    int i, j, k;
    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.
    double xMin;
    int steps;
    node2_t **Loop=new node2_t*[2*max_n], *CurX, *LeaveX;
 
#if DEBUG_LEVEL > 3
    printf("EnterX = (%d,%d)\n", _EnterX->i, _EnterX->j);
#endif

    /* ENTER THE NEW BASIC VARIABLE */
    i = _EnterX->i;
    j = _EnterX->j;
    _IsX[i][j] = 1;
    _EnterX->NextC = _RowsX[i];
    _EnterX->NextR = _ColsX[j];
    _EnterX->val = 0;
    _RowsX[i] = _EnterX;
    _ColsX[j] = _EnterX;

    /* FIND A CHAIN REACTION */
    steps = findLoop(Loop, _n1, _n2, _XV, _RowsX, _ColsX);

    /* FIND THE LARGEST VALUE IN THE LOOP */
    xMin = INFINITY;
    for (k=1; k < steps; k+=2)
      {
        if (Loop[k]->val < xMin)
          {
            LeaveX = Loop[k];
            xMin = Loop[k]->val;
          }
      }

    /* UPDATE THE LOOP */
    for (k=0; k < steps; k+=2)
      {
        Loop[k]->val += xMin;
        Loop[k+1]->val -= xMin;
      }

#if DEBUG_LEVEL > 3
    printf("LeaveX = (%d,%d)\n", LeaveX->i, LeaveX->j);
#endif

    /* REMOVE THE LEAVING BASIC VARIABLE */
    i = LeaveX->i;
    j = LeaveX->j;
    _IsX[i][j] = 0;
    if (_RowsX[i] == LeaveX)
      _RowsX[i] = LeaveX->NextC;
    else
      for (CurX=_RowsX[i]; CurX != NULL; CurX = CurX->NextC)
        if (CurX->NextC == LeaveX)
          {
            CurX->NextC = CurX->NextC->NextC;
            break;
          }
    if (_ColsX[j] == LeaveX)
      _ColsX[j] = LeaveX->NextR;
    else
      for (CurX=_ColsX[j]; CurX != NULL; CurX = CurX->NextR)
        if (CurX->NextR == LeaveX)
          {
            CurX->NextR = CurX->NextR->NextR;
            break;
          }

    /* SET _EnterX TO BE THE NEW EMPTY SLOT */
    _EnterX = LeaveX;

        delete[] Loop;
}



/**********************
    findLoop
**********************/
static int findLoop(node2_t **Loop, int _n1, int _n2, node2_t *_XV, node2_t **_RowsX, node2_t **_ColsX)
{
  int i, steps;
  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.
  node2_t **CurX, *NewX;
  char *IsUsed=new char[2*max_n]; 
 
  for (i=0; i < _n1+_n2; i++)
    IsUsed[i] = 0;

  CurX = Loop;
  NewX = *CurX = _EnterX;
  IsUsed[_EnterX-_XV] = 1;
  steps = 1;

  do
    {
      if (steps%2 == 1)
        {
          /* FIND AN UNUSED X IN THE ROW */
          NewX = _RowsX[NewX->i];
          while (NewX != NULL && IsUsed[NewX-_XV])
            NewX = NewX->NextC;
        }
      else
        {
          /* FIND AN UNUSED X IN THE COLUMN, OR THE ENTERING X */
          NewX = _ColsX[NewX->j];
          while (NewX != NULL && IsUsed[NewX-_XV] && NewX != _EnterX)
            NewX = NewX->NextR;
          if (NewX == _EnterX)
            break;
        }

     if (NewX != NULL)  /* FOUND THE NEXT X */
       {
         /* ADD X TO THE LOOP */
         *++CurX = NewX;
         IsUsed[NewX-_XV] = 1;
         steps++;
#if DEBUG_LEVEL > 3
         printf("steps=%d, NewX=(%d,%d)\n", steps, NewX->i, NewX->j);    
#endif
       }
     else  /* DIDN'T FIND THE NEXT X */
       {
         /* BACKTRACK */
         do
           {
             NewX = *CurX;
             do 
               {
                 if (steps%2 == 1)
                   NewX = NewX->NextR;
                 else
                   NewX = NewX->NextC;
               } while (NewX != NULL && IsUsed[NewX-_XV]);
             
             if (NewX == NULL)
               {
                 IsUsed[*CurX-_XV] = 0;
                 CurX--;
                 steps--;
               }
         } while (NewX == NULL && CurX >= Loop);
         
#if DEBUG_LEVEL > 3
         printf("BACKTRACKING TO: steps=%d, NewX=(%d,%d)\n",
                steps, NewX->i, NewX->j);    
#endif
           IsUsed[*CurX-_XV] = 0;
           *CurX = NewX;
           IsUsed[NewX-_XV] = 1;
       }     
    } while(CurX >= Loop);
  
  if (CurX == Loop)
    {
      fprintf(stderr, "emd: Unexpected error in findLoop!\n");
      exit(1);
    }
#if DEBUG_LEVEL > 3
  printf("FOUND LOOP:\n");
  for (i=0; i < steps; i++)
    printf("%d: (%d,%d)\n", i, Loop[i]->i, Loop[i]->j);
#endif

  delete[] IsUsed;

  return steps;
}



/**********************
    russel
**********************/
static void russel(double *S, double *D, int _n1, int _n2, float **_CM, char **_IsX, node2_t **_RowsX, node2_t **_ColsX)
{
  int i, j, found, minI, minJ;
  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.
  double deltaMin, oldVal, diff;
  double** Delta = new double*[_n1];
  for(int k = 0; k < _n1; ++k)
    Delta[k] = new double[_n2];
  node1_t *Ur=new node1_t[max_n], *Vr=new node1_t[max_n];
  node1_t uHead, *CurU, *PrevU;
  node1_t vHead, *CurV, *PrevV;
  node1_t *PrevUMinI, *PrevVMinJ, *Remember;

  /* INITIALIZE THE ROWS LIST (Ur), AND THE COLUMNS LIST (Vr) */
  uHead.Next = CurU = Ur;
  for (i=0; i < _n1; i++)
    {
      CurU->i = i;
      CurU->val = -INFINITY;
      CurU->Next = CurU+1;
      CurU++;
    }
  (--CurU)->Next = NULL;
  
  vHead.Next = CurV = Vr;
  for (j=0; j < _n2; j++)
    {
      CurV->i = j;
      CurV->val = -INFINITY;
      CurV->Next = CurV+1;
      CurV++;
    }
  (--CurV)->Next = NULL;
  
  /* FIND THE MAXIMUM ROW AND COLUMN VALUES (Ur[i] AND Vr[j]) */
  for(i=0; i < _n1 ; i++)
    for(j=0; j < _n2 ; j++)
      {
        float v;
        v = _CM[i][j];
        if (Ur[i].val <= v)
          Ur[i].val = v;
        if (Vr[j].val <= v)
          Vr[j].val = v;
      }
  
  /* COMPUTE THE Delta MATRIX */
  for(i=0; i < _n1 ; i++)
    for(j=0; j < _n2 ; j++)
      Delta[i][j] = _CM[i][j] - Ur[i].val - Vr[j].val;

  /* FIND THE BASIC VARIABLES */
  do
    {
#if DEBUG_LEVEL > 3
      printf("Ur=");
      for(CurU = uHead.Next; CurU != NULL; CurU = CurU->Next)
        printf("[%d]",CurU-Ur);
      printf("\n");
      printf("Vr=");
      for(CurV = vHead.Next; CurV != NULL; CurV = CurV->Next)
        printf("[%d]",CurV-Vr);
      printf("\n");
      printf("\n\n");
#endif
 
      /* FIND THE SMALLEST Delta[i][j] */
      found = 0; 
      deltaMin = INFINITY;      
      PrevU = &uHead;
      for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
        {
          int i;
          i = CurU->i;
          PrevV = &vHead;
          for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
            {
              int j;
              j = CurV->i;
              if (deltaMin > Delta[i][j])
                {
                  deltaMin = Delta[i][j];
                  minI = i;
                  minJ = j;
                  PrevUMinI = PrevU;
                  PrevVMinJ = PrevV;
                  found = 1;
                }
              PrevV = CurV;
            }
          PrevU = CurU;
        }
      
      if (! found)
        break;

      /* ADD X[minI][minJ] TO THE BASIS, AND ADJUST SUPPLIES AND COST */
      Remember = PrevUMinI->Next;
      addBasicVariable(minI, minJ, S, D, PrevUMinI, PrevVMinJ, &uHead, _IsX, _RowsX, _ColsX);

      /* UPDATE THE NECESSARY Delta[][] */
      if (Remember == PrevUMinI->Next)  /* LINE minI WAS DELETED */
        {
          for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
            {
              int j;
              j = CurV->i;
              if (CurV->val == _CM[minI][j])  /* COLUMN j NEEDS UPDATING */
                {
                  /* FIND THE NEW MAXIMUM VALUE IN THE COLUMN */
                  oldVal = CurV->val;
                  CurV->val = -INFINITY;
                  for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
                    {
                      int i;
                      i = CurU->i;
                      if (CurV->val <= _CM[i][j])
                        CurV->val = _CM[i][j];
                    }
                  
                  /* IF NEEDED, ADJUST THE RELEVANT Delta[*][j] */
                  diff = oldVal - CurV->val;
                  if (fabs(diff) < EPSILON * _maxC)
                    for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
                      Delta[CurU->i][j] += diff;
                }
            }
        }
      else  /* COLUMN minJ WAS DELETED */
        {
          for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
            {
              int i;
              i = CurU->i;
              if (CurU->val == _CM[i][minJ])  /* ROW i NEEDS UPDATING */
                {
                  /* FIND THE NEW MAXIMUM VALUE IN THE ROW */
                  oldVal = CurU->val;
                  CurU->val = -INFINITY;
                  for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
                    {
                      int j;
                      j = CurV->i;
                      if(CurU->val <= _CM[i][j])
                        CurU->val = _CM[i][j];
                    }
                  
                  /* If NEEDED, ADJUST THE RELEVANT Delta[i][*] */
                  diff = oldVal - CurU->val;
                  if (fabs(diff) < EPSILON * _maxC)
                    for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
                      Delta[i][CurV->i] += diff;
                }
            }
        }
    } while (uHead.Next != NULL || vHead.Next != NULL);

        delete[] Ur;
        delete[] Vr;
    for(int k = 0; k < _n1; ++k)
      delete[] Delta[k];
    delete[] Delta;
}




/**********************
    addBasicVariable
**********************/
static void addBasicVariable(int minI, int minJ, double *S, double *D, 
                             node1_t *PrevUMinI, node1_t *PrevVMinJ,
                             node1_t *UHead, char **_IsX, node2_t **_RowsX, node2_t **_ColsX)
{
  double T;
  
  if (fabs(S[minI]-D[minJ]) <= EPSILON * _maxW)  /* DEGENERATE CASE */
    {
      T = S[minI];
      S[minI] = 0;
      D[minJ] -= T; 
    }
  else if (S[minI] < D[minJ])  /* SUPPLY EXHAUSTED */
    {
      T = S[minI];
      S[minI] = 0;
      D[minJ] -= T; 
    }
  else  /* DEMAND EXHAUSTED */
    {
      T = D[minJ];
      D[minJ] = 0; 
      S[minI] -= T; 
    }

  /* X(minI,minJ) IS A BASIC VARIABLE */
  _IsX[minI][minJ] = 1; 

  _EndX->val = T;
  _EndX->i = minI;
  _EndX->j = minJ;
  _EndX->NextC = _RowsX[minI];
  _EndX->NextR = _ColsX[minJ];
  _RowsX[minI] = _EndX;
  _ColsX[minJ] = _EndX;
  _EndX++;

  /* DELETE SUPPLY ROW ONLY IF THE EMPTY, AND IF NOT LAST ROW */
  if (S[minI] == 0 && UHead->Next->Next != NULL)
    PrevUMinI->Next = PrevUMinI->Next->Next;  /* REMOVE ROW FROM LIST */
  else
    PrevVMinJ->Next = PrevVMinJ->Next->Next;  /* REMOVE COLUMN FROM LIST */
}