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

Last change on this file since 1061 was 1061, checked in by oriona, 3 years ago

Static array removed.

File size: 19.1 KB
Line
1/*
2    emd.c
3
4    Last update: 3/14/98
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, McGraw-Hill, 1990.
10
11    Copyright (C) 1998 Yossi Rubner
12    Computer Science Department, Stanford University
13    E-Mail: rubner@cs.stanford.edu   URL: http://vision.stanford.edu/~rubner
14*/
15
16#include <stdio.h>
17#include <stdlib.h>
18#include <math.h>
19
20#include "emd.h"
21
22#define DEBUG_LEVEL 0
23/*
24 DEBUG_LEVEL:
25   0 = NO MESSAGES
26   1 = PRINT THE NUMBER OF ITERATIONS AND THE FINAL RESULT
27   2 = PRINT THE RESULT AFTER EVERY ITERATION
28   3 = PRINT ALSO THE FLOW AFTER EVERY ITERATION
29   4 = PRINT A LOT OF INFORMATION (PROBABLY USEFUL ONLY FOR THE AUTHOR)
30*/
31
32
33#define MAX_SIG_SIZE1  (MAX_SIG_SIZE+1)  /* FOR THE POSIBLE DUMMY FEATURE */
34
35/* NEW TYPES DEFINITION */
36
37/* node1_t IS USED FOR SINGLE-LINKED LISTS */
38typedef struct node1_t {
39  int i;
40  double val;
41  struct node1_t *Next;
42} node1_t;
43
44/* node1_t IS USED FOR DOUBLE-LINKED LISTS */
45typedef 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 */
55static int _n1, _n2;                          /* SIGNATURES SIZES */
56static float _CM[MAX_SIG_SIZE1][MAX_SIG_SIZE1];/* THE COST MATRIX */
57static node2_t _XV[MAX_SIG_SIZE1*2];            /* THE BASIC VARIABLES VECTOR */
58/* VARIABLES TO HANDLE _X EFFICIENTLY */
59static node2_t *_EndX, *_EnterX;
60static char _IsX[MAX_SIG_SIZE1][MAX_SIG_SIZE1];
61static node2_t *_RowsX[MAX_SIG_SIZE1], *_ColsX[MAX_SIG_SIZE1];
62static double _maxW;
63static float _maxC;
64
65/* DECLARATION OF FUNCTIONS */
66static float init(signature_t *Signature1, signature_t *Signature2,
67                  float (*Dist)(feature_t *, feature_t *));
68static void findBasicVariables(node1_t *U, node1_t *V);
69static int isOptimal(node1_t *U, node1_t *V);
70static int findLoop(node2_t **Loop);
71static void newSol();
72static void russel(double *S, double *D);
73static void addBasicVariable(int minI, int minJ, double *S, double *D,
74                             node1_t *PrevUMinI, node1_t *PrevVMinJ,
76#if DEBUG_LEVEL > 0
77static void printSolution();
78#endif
79
80
81/******************************************************************************
82float emd(signature_t *Signature1, signature_t *Signature2,
83          float (*Dist)(feature_t *, feature_t *), flow_t *Flow, int *FlowSize)
84
85where
86
87   Signature1, Signature2  Pointers to signatures that their distance we want
88              to compute.
89   Dist       Pointer to the ground distance. i.e. the function that computes
90              the distance between two features.
91   Flow       (Optional) Pointer to a vector of flow_t (defined in emd.h)
92              where the resulting flow will be stored. Flow must have n1+n2-1
93              elements, where n1 and n2 are the sizes of the two signatures
94              respectively.
95              If NULL, the flow is not returned.
96   FlowSize   (Optional) Pointer to an integer where the number of elements in
97              Flow will be stored
98
99******************************************************************************/
100
101float emd(signature_t *Signature1, signature_t *Signature2,
102          float (*Dist)(feature_t *, feature_t *),
103          flow_t *Flow, int *FlowSize)
104{
105  int itr;
106  double totalCost;
107  float w;
108  node2_t *XP;
109  flow_t *FlowP;
110  node1_t U[MAX_SIG_SIZE1], V[MAX_SIG_SIZE1];
111
112  w = init(Signature1, Signature2, Dist);
113
114#if DEBUG_LEVEL > 1
115  printf("\nINITIAL SOLUTION:\n");
116  printSolution();
117#endif
118
119  if (_n1 > 1 && _n2 > 1)  /* IF _n1 = 1 OR _n2 = 1 THEN WE ARE DONE */
120    {
121      for (itr = 1; itr < MAX_ITERATIONS; itr++)
122        {
123          /* FIND BASIC VARIABLES */
124          findBasicVariables(U, V);
125
126          /* CHECK FOR OPTIMALITY */
127          if (isOptimal(U, V))
128            break;
129
130          /* IMPROVE SOLUTION */
131          newSol();
132
133#if DEBUG_LEVEL > 1
134          printf("\nITERATION # %d \n", itr);
135          printSolution();
136#endif
137        }
138
139      if (itr == MAX_ITERATIONS)
140        fprintf(stderr, "emd: Maximum number of iterations has been reached (%d)\n",
141                MAX_ITERATIONS);
142    }
143
144  /* COMPUTE THE TOTAL FLOW */
145  totalCost = 0;
146  if (Flow != NULL)
147    FlowP = Flow;
148  for(XP=_XV; XP < _EndX; XP++)
149    {
150      if (XP == _EnterX)  /* _EnterX IS THE EMPTY SLOT */
151        continue;
152      if (XP->i == Signature1->n || XP->j == Signature2->n)  /* DUMMY FEATURE */
153        continue;
154
155      if (XP->val == 0)  /* ZERO FLOW */
156        continue;
157
158      totalCost += (double)XP->val * _CM[XP->i][XP->j];
159      if (Flow != NULL)
160        {
161          FlowP->from = XP->i;
162          FlowP->to = XP->j;
163          FlowP->amount = XP->val;
164          FlowP++;
165        }
166    }
167  if (Flow != NULL)
168    *FlowSize = FlowP-Flow;
169
170#if DEBUG_LEVEL > 0
171  printf("\n*** OPTIMAL SOLUTION (%d ITERATIONS): %f ***\n", itr, totalCost);
172#endif
173
174  /* RETURN THE NORMALIZED COST == EMD */
175  return (float)(totalCost / w);
176}
177
178
179
180
181
182/**********************
183   init
184**********************/
185static float init(signature_t *Signature1, signature_t *Signature2,
186                  float (*Dist)(feature_t *, feature_t *))
187{
188  int i, j;
189  double sSum, dSum, diff;
190  feature_t *P1, *P2;
191  double S[MAX_SIG_SIZE1], D[MAX_SIG_SIZE1];
192
193  _n1 = Signature1->n;
194  _n2 = Signature2->n;
195
196  if (_n1 > MAX_SIG_SIZE || _n2 > MAX_SIG_SIZE)
197    {
198      fprintf(stderr, "emd: Signature size is limited to %d\n", MAX_SIG_SIZE);
199      exit(1);
200    }
201
202  /* COMPUTE THE DISTANCE MATRIX */
203  _maxC = 0;
204  for(i=0, P1=Signature1->Features; i < _n1; i++, P1++)
205    for(j=0, P2=Signature2->Features; j < _n2; j++, P2++)
206      {
207        _CM[i][j] = Dist(P1, P2);
208        if (_CM[i][j] > _maxC)
209          _maxC = _CM[i][j];
210      }
211
212  /* SUM UP THE SUPPLY AND DEMAND */
213  sSum = 0.0;
214  for(i=0; i < _n1; i++)
215    {
216      S[i] = Signature1->Weights[i];
217      sSum += Signature1->Weights[i];
218      _RowsX[i] = NULL;
219    }
220  dSum = 0.0;
221  for(j=0; j < _n2; j++)
222    {
223      D[j] = Signature2->Weights[j];
224      dSum += Signature2->Weights[j];
225      _ColsX[j] = NULL;
226    }
227
228  /* IF SUPPLY DIFFERENT THAN THE DEMAND, ADD A ZERO-COST DUMMY CLUSTER */
229  diff = sSum - dSum;
230  if (fabs(diff) >= EPSILON * sSum)
231    {
232      if (diff < 0.0)
233        {
234          for (j=0; j < _n2; j++)
235            _CM[_n1][j] = 0;
236          S[_n1] = -diff;
237          _RowsX[_n1] = NULL;
238          _n1++;
239        }
240      else
241        {
242          for (i=0; i < _n1; i++)
243            _CM[i][_n2] = 0;
244          D[_n2] = diff;
245          _ColsX[_n2] = NULL;
246          _n2++;
247        }
248    }
249
250  /* INITIALIZE THE BASIC VARIABLE STRUCTURES */
251  for (i=0; i < _n1; i++)
252    for (j=0; j < _n2; j++)
253        _IsX[i][j] = 0;
254  _EndX = _XV;
255
256  _maxW = sSum > dSum ? sSum : dSum;
257
258  /* FIND INITIAL SOLUTION */
259  russel(S, D);
260
261  _EnterX = _EndX++;  /* AN EMPTY SLOT (ONLY _n1+_n2-1 BASIC VARIABLES) */
262
263  return sSum > dSum ? dSum : sSum;
264}
265
266
267/**********************
268    findBasicVariables
269 **********************/
270static void findBasicVariables(node1_t *U, node1_t *V)
271{
272  int i, j, found;
273  int UfoundNum, VfoundNum;
276
277  /* INITIALIZE THE ROWS LIST (U) AND THE COLUMNS LIST (V) */
278  u0Head.Next = CurU = U;
279  for (i=0; i < _n1; i++)
280    {
281      CurU->i = i;
282      CurU->Next = CurU+1;
283      CurU++;
284    }
285  (--CurU)->Next = NULL;
287
288  CurV = V+1;
289  v0Head.Next = _n2 > 1 ? V+1 : NULL;
290  for (j=1; j < _n2; j++)
291    {
292      CurV->i = j;
293      CurV->Next = CurV+1;
294      CurV++;
295    }
296  (--CurV)->Next = NULL;
298
299  /* THERE ARE _n1+_n2 VARIABLES BUT ONLY _n1+_n2-1 INDEPENDENT EQUATIONS,
300     SO SET V[0]=0 */
301  V[0].i = 0;
302  V[0].val = 0;
305
306  /* LOOP UNTIL ALL VARIABLES ARE FOUND */
307  UfoundNum=VfoundNum=0;
308  while (UfoundNum < _n1 || VfoundNum < _n2)
309    {
310
311#if DEBUG_LEVEL > 3
312      printf("UfoundNum=%d/%d,VfoundNum=%d/%d\n",UfoundNum,_n1,VfoundNum,_n2);
313      printf("U0=");
314      for(CurU = u0Head.Next; CurU != NULL; CurU = CurU->Next)
315        printf("[%d]",CurU-U);
316      printf("\n");
317      printf("U1=");
318      for(CurU = u1Head.Next; CurU != NULL; CurU = CurU->Next)
319        printf("[%d]",CurU-U);
320      printf("\n");
321      printf("V0=");
322      for(CurV = v0Head.Next; CurV != NULL; CurV = CurV->Next)
323        printf("[%d]",CurV-V);
324      printf("\n");
325      printf("V1=");
326      for(CurV = v1Head.Next; CurV != NULL; CurV = CurV->Next)
327        printf("[%d]",CurV-V);
328      printf("\n\n");
329#endif
330
331      found = 0;
332      if (VfoundNum < _n2)
333        {
334          /* LOOP OVER ALL MARKED COLUMNS */
336          for (CurV=v1Head.Next; CurV != NULL; CurV=CurV->Next)
337            {
338              j = CurV->i;
339              /* FIND THE VARIABLES IN COLUMN j */
341              for (CurU=u0Head.Next; CurU != NULL; CurU=CurU->Next)
342                {
343                  i = CurU->i;
344                  if (_IsX[i][j])
345                    {
346                      /* COMPUTE U[i] */
347                      CurU->val = _CM[i][j] - CurV->val;
348                      /* ...AND ADD IT TO THE MARKED LIST */
349                      PrevU->Next = CurU->Next;
352                      CurU = PrevU;
353                    }
354                  else
355                    PrevU = CurU;
356                }
357              PrevV->Next = CurV->Next;
358              VfoundNum++;
359              found = 1;
360            }
361        }
362     if (UfoundNum < _n1)
363        {
364          /* LOOP OVER ALL MARKED ROWS */
366          for (CurU=u1Head.Next; CurU != NULL; CurU=CurU->Next)
367            {
368              i = CurU->i;
369              /* FIND THE VARIABLES IN ROWS i */
371              for (CurV=v0Head.Next; CurV != NULL; CurV=CurV->Next)
372                {
373                  j = CurV->i;
374                  if (_IsX[i][j])
375                    {
376                      /* COMPUTE V[j] */
377                      CurV->val = _CM[i][j] - CurU->val;
378                      /* ...AND ADD IT TO THE MARKED LIST */
379                      PrevV->Next = CurV->Next;
382                      CurV = PrevV;
383                    }
384                  else
385                    PrevV = CurV;
386                }
387              PrevU->Next = CurU->Next;
388              UfoundNum++;
389              found = 1;
390            }
391        }
392     if (! found)
393       {
394         fprintf(stderr, "emd: Unexpected error in findBasicVariables!\n");
395         fprintf(stderr, "This typically happens when the EPSILON defined in\n");
396         fprintf(stderr, "emd.h is not right for the scale of the problem.\n");
397         exit(1);
398       }
399    }
400}
401
402
403
404
405/**********************
406    isOptimal
407 **********************/
408static int isOptimal(node1_t *U, node1_t *V)
409{
410  double delta, deltaMin;
411  int i, j, minI, minJ;
412
413  /* FIND THE MINIMAL Cij-Ui-Vj OVER ALL i,j */
414  deltaMin = INFINITY;
415  for(i=0; i < _n1; i++)
416    for(j=0; j < _n2; j++)
417      if (! _IsX[i][j])
418        {
419          delta = _CM[i][j] - U[i].val - V[j].val;
420          if (deltaMin > delta)
421            {
422              deltaMin = delta;
423              minI = i;
424              minJ = j;
425            }
426        }
427
428#if DEBUG_LEVEL > 3
429  printf("deltaMin=%f\n", deltaMin);
430#endif
431
432   if (deltaMin == INFINITY)
433     {
434       fprintf(stderr, "emd: Unexpected error in isOptimal.\n");
435       exit(0);
436     }
437
438   _EnterX->i = minI;
439   _EnterX->j = minJ;
440
441   /* IF NO NEGATIVE deltaMin, WE FOUND THE OPTIMAL SOLUTION */
442   return deltaMin >= -EPSILON * _maxC;
443
444/*
445   return deltaMin >= -EPSILON;
446 */
447}
448
449
450/**********************
451    newSol
452**********************/
453static void newSol()
454{
455    int i, j, k;
456    double xMin;
457    int steps;
458    node2_t *Loop[2*MAX_SIG_SIZE1], *CurX, *LeaveX;
459
460#if DEBUG_LEVEL > 3
461    printf("EnterX = (%d,%d)\n", _EnterX->i, _EnterX->j);
462#endif
463
464    /* ENTER THE NEW BASIC VARIABLE */
465    i = _EnterX->i;
466    j = _EnterX->j;
467    _IsX[i][j] = 1;
468    _EnterX->NextC = _RowsX[i];
469    _EnterX->NextR = _ColsX[j];
470    _EnterX->val = 0;
471    _RowsX[i] = _EnterX;
472    _ColsX[j] = _EnterX;
473
474    /* FIND A CHAIN REACTION */
475    steps = findLoop(Loop);
476
477    /* FIND THE LARGEST VALUE IN THE LOOP */
478    xMin = INFINITY;
479    for (k=1; k < steps; k+=2)
480      {
481        if (Loop[k]->val < xMin)
482          {
483            LeaveX = Loop[k];
484            xMin = Loop[k]->val;
485          }
486      }
487
488    /* UPDATE THE LOOP */
489    for (k=0; k < steps; k+=2)
490      {
491        Loop[k]->val += xMin;
492        Loop[k+1]->val -= xMin;
493      }
494
495#if DEBUG_LEVEL > 3
496    printf("LeaveX = (%d,%d)\n", LeaveX->i, LeaveX->j);
497#endif
498
499    /* REMOVE THE LEAVING BASIC VARIABLE */
500    i = LeaveX->i;
501    j = LeaveX->j;
502    _IsX[i][j] = 0;
503    if (_RowsX[i] == LeaveX)
504      _RowsX[i] = LeaveX->NextC;
505    else
506      for (CurX=_RowsX[i]; CurX != NULL; CurX = CurX->NextC)
507        if (CurX->NextC == LeaveX)
508          {
509            CurX->NextC = CurX->NextC->NextC;
510            break;
511          }
512    if (_ColsX[j] == LeaveX)
513      _ColsX[j] = LeaveX->NextR;
514    else
515      for (CurX=_ColsX[j]; CurX != NULL; CurX = CurX->NextR)
516        if (CurX->NextR == LeaveX)
517          {
518            CurX->NextR = CurX->NextR->NextR;
519            break;
520          }
521
522    /* SET _EnterX TO BE THE NEW EMPTY SLOT */
523    _EnterX = LeaveX;
524}
525
526
527
528/**********************
529    findLoop
530**********************/
531static int findLoop(node2_t **Loop)
532{
533  int i, steps;
534  node2_t **CurX, *NewX;
535  char IsUsed[2*MAX_SIG_SIZE1];
536
537  for (i=0; i < _n1+_n2; i++)
538    IsUsed[i] = 0;
539
540  CurX = Loop;
541  NewX = *CurX = _EnterX;
542  IsUsed[_EnterX-_XV] = 1;
543  steps = 1;
544
545  do
546    {
547      if (steps%2 == 1)
548        {
549          /* FIND AN UNUSED X IN THE ROW */
550          NewX = _RowsX[NewX->i];
551          while (NewX != NULL && IsUsed[NewX-_XV])
552            NewX = NewX->NextC;
553        }
554      else
555        {
556          /* FIND AN UNUSED X IN THE COLUMN, OR THE ENTERING X */
557          NewX = _ColsX[NewX->j];
558          while (NewX != NULL && IsUsed[NewX-_XV] && NewX != _EnterX)
559            NewX = NewX->NextR;
560          if (NewX == _EnterX)
561            break;
562        }
563
564     if (NewX != NULL)  /* FOUND THE NEXT X */
565       {
566         /* ADD X TO THE LOOP */
567         *++CurX = NewX;
568         IsUsed[NewX-_XV] = 1;
569         steps++;
570#if DEBUG_LEVEL > 3
571         printf("steps=%d, NewX=(%d,%d)\n", steps, NewX->i, NewX->j);
572#endif
573       }
574     else  /* DIDN'T FIND THE NEXT X */
575       {
576         /* BACKTRACK */
577         do
578           {
579             NewX = *CurX;
580             do
581               {
582                 if (steps%2 == 1)
583                   NewX = NewX->NextR;
584                 else
585                   NewX = NewX->NextC;
586               } while (NewX != NULL && IsUsed[NewX-_XV]);
587
588             if (NewX == NULL)
589               {
590                 IsUsed[*CurX-_XV] = 0;
591                 CurX--;
592                 steps--;
593               }
594         } while (NewX == NULL && CurX >= Loop);
595
596#if DEBUG_LEVEL > 3
597         printf("BACKTRACKING TO: steps=%d, NewX=(%d,%d)\n",
598                steps, NewX->i, NewX->j);
599#endif
600           IsUsed[*CurX-_XV] = 0;
601           *CurX = NewX;
602           IsUsed[NewX-_XV] = 1;
603       }
604    } while(CurX >= Loop);
605
606  if (CurX == Loop)
607    {
608      fprintf(stderr, "emd: Unexpected error in findLoop!\n");
609      exit(1);
610    }
611#if DEBUG_LEVEL > 3
612  printf("FOUND LOOP:\n");
613  for (i=0; i < steps; i++)
614    printf("%d: (%d,%d)\n", i, Loop[i]->i, Loop[i]->j);
615#endif
616
617  return steps;
618}
619
620
621
622/**********************
623    russel
624**********************/
625static void russel(double *S, double *D)
626{
627  int i, j, found, minI, minJ;
628  double deltaMin, oldVal, diff;
629  double** Delta = new double*[_n1];
630  for(int k = 0; k < _n1; ++k)
631    Delta[k] = new double[_n2];
632  node1_t Ur[MAX_SIG_SIZE1], Vr[MAX_SIG_SIZE1];
635  node1_t *PrevUMinI, *PrevVMinJ, *Remember;
636
637  /* INITIALIZE THE ROWS LIST (Ur), AND THE COLUMNS LIST (Vr) */
638  uHead.Next = CurU = Ur;
639  for (i=0; i < _n1; i++)
640    {
641      CurU->i = i;
642      CurU->val = -INFINITY;
643      CurU->Next = CurU+1;
644      CurU++;
645    }
646  (--CurU)->Next = NULL;
647
648  vHead.Next = CurV = Vr;
649  for (j=0; j < _n2; j++)
650    {
651      CurV->i = j;
652      CurV->val = -INFINITY;
653      CurV->Next = CurV+1;
654      CurV++;
655    }
656  (--CurV)->Next = NULL;
657
658  /* FIND THE MAXIMUM ROW AND COLUMN VALUES (Ur[i] AND Vr[j]) */
659  for(i=0; i < _n1 ; i++)
660    for(j=0; j < _n2 ; j++)
661      {
662        float v;
663        v = _CM[i][j];
664        if (Ur[i].val <= v)
665          Ur[i].val = v;
666        if (Vr[j].val <= v)
667          Vr[j].val = v;
668      }
669
670  /* COMPUTE THE Delta MATRIX */
671  for(i=0; i < _n1 ; i++)
672    for(j=0; j < _n2 ; j++)
673      Delta[i][j] = _CM[i][j] - Ur[i].val - Vr[j].val;
674
675  /* FIND THE BASIC VARIABLES */
676  do
677    {
678#if DEBUG_LEVEL > 3
679      printf("Ur=");
680      for(CurU = uHead.Next; CurU != NULL; CurU = CurU->Next)
681        printf("[%d]",CurU-Ur);
682      printf("\n");
683      printf("Vr=");
684      for(CurV = vHead.Next; CurV != NULL; CurV = CurV->Next)
685        printf("[%d]",CurV-Vr);
686      printf("\n");
687      printf("\n\n");
688#endif
689
690      /* FIND THE SMALLEST Delta[i][j] */
691      found = 0;
692      deltaMin = INFINITY;
694      for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
695        {
696          int i;
697          i = CurU->i;
699          for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
700            {
701              int j;
702              j = CurV->i;
703              if (deltaMin > Delta[i][j])
704                {
705                  deltaMin = Delta[i][j];
706                  minI = i;
707                  minJ = j;
708                  PrevUMinI = PrevU;
709                  PrevVMinJ = PrevV;
710                  found = 1;
711                }
712              PrevV = CurV;
713            }
714          PrevU = CurU;
715        }
716
717      if (! found)
718        break;
719
720      /* ADD X[minI][minJ] TO THE BASIS, AND ADJUST SUPPLIES AND COST */
721      Remember = PrevUMinI->Next;
723
724      /* UPDATE THE NECESSARY Delta[][] */
725      if (Remember == PrevUMinI->Next)  /* LINE minI WAS DELETED */
726        {
727          for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
728            {
729              int j;
730              j = CurV->i;
731              if (CurV->val == _CM[minI][j])  /* COLUMN j NEEDS UPDATING */
732                {
733                  /* FIND THE NEW MAXIMUM VALUE IN THE COLUMN */
734                  oldVal = CurV->val;
735                  CurV->val = -INFINITY;
736                  for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
737                    {
738                      int i;
739                      i = CurU->i;
740                      if (CurV->val <= _CM[i][j])
741                        CurV->val = _CM[i][j];
742                    }
743
744                  /* IF NEEDED, ADJUST THE RELEVANT Delta[*][j] */
745                  diff = oldVal - CurV->val;
746                  if (fabs(diff) < EPSILON * _maxC)
747                    for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
748                      Delta[CurU->i][j] += diff;
749                }
750            }
751        }
752      else  /* COLUMN minJ WAS DELETED */
753        {
754          for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
755            {
756              int i;
757              i = CurU->i;
758              if (CurU->val == _CM[i][minJ])  /* ROW i NEEDS UPDATING */
759                {
760                  /* FIND THE NEW MAXIMUM VALUE IN THE ROW */
761                  oldVal = CurU->val;
762                  CurU->val = -INFINITY;
763                  for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
764                    {
765                      int j;
766                      j = CurV->i;
767                      if(CurU->val <= _CM[i][j])
768                        CurU->val = _CM[i][j];
769                    }
770
771                  /* If NEEDED, ADJUST THE RELEVANT Delta[i][*] */
772                  diff = oldVal - CurU->val;
773                  if (fabs(diff) < EPSILON * _maxC)
774                    for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
775                      Delta[i][CurV->i] += diff;
776                }
777            }
778        }
780    for(int k = 0; k < _n1; ++k)
781      delete [] Delta[k];
782    delete [] Delta;
783}
784
785
786
787
788/**********************
790**********************/
791static void addBasicVariable(int minI, int minJ, double *S, double *D,
792                             node1_t *PrevUMinI, node1_t *PrevVMinJ,
794{
795  double T;
796
797  if (fabs(S[minI]-D[minJ]) <= EPSILON * _maxW)  /* DEGENERATE CASE */
798    {
799      T = S[minI];
800      S[minI] = 0;
801      D[minJ] -= T;
802    }
803  else if (S[minI] < D[minJ])  /* SUPPLY EXHAUSTED */
804    {
805      T = S[minI];
806      S[minI] = 0;
807      D[minJ] -= T;
808    }
809  else  /* DEMAND EXHAUSTED */
810    {
811      T = D[minJ];
812      D[minJ] = 0;
813      S[minI] -= T;
814    }
815
816  /* X(minI,minJ) IS A BASIC VARIABLE */
817  _IsX[minI][minJ] = 1;
818
819  _EndX->val = T;
820  _EndX->i = minI;
821  _EndX->j = minJ;
822  _EndX->NextC = _RowsX[minI];
823  _EndX->NextR = _ColsX[minJ];
824  _RowsX[minI] = _EndX;
825  _ColsX[minJ] = _EndX;
826  _EndX++;
827
828  /* DELETE SUPPLY ROW ONLY IF THE EMPTY, AND IF NOT LAST ROW */
829  if (S[minI] == 0 && UHead->Next->Next != NULL)
830    PrevUMinI->Next = PrevUMinI->Next->Next;  /* REMOVE ROW FROM LIST */
831  else
832    PrevVMinJ->Next = PrevVMinJ->Next->Next;  /* REMOVE COLUMN FROM LIST */
833}
Note: See TracBrowser for help on using the repository browser.