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

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

Similarity measures code refactored. Distribution-based similarity measure added.

File size: 18.9 KB
RevLine 
[1044]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 _C[MAX_SIG_SIZE1][MAX_SIG_SIZE1];/* THE COST MATRIX */
57static node2_t _X[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,
75                             node1_t *UHead);
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=_X; 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 * _C[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        _C[i][j] = Dist(P1, P2);
208        if (_C[i][j] > _maxC)
209          _maxC = _C[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            _C[_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            _C[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 = _X;
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;
274  node1_t u0Head, u1Head, *CurU, *PrevU;
275  node1_t v0Head, v1Head, *CurV, *PrevV;
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;
286  u1Head.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;
297  v1Head.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;
303  v1Head.Next = V;
304  v1Head.Next->Next = NULL;
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 */
335          PrevV = &v1Head;
336          for (CurV=v1Head.Next; CurV != NULL; CurV=CurV->Next)
337            {
338              j = CurV->i;
339              /* FIND THE VARIABLES IN COLUMN j */
340              PrevU = &u0Head;
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 = _C[i][j] - CurV->val;
348                      /* ...AND ADD IT TO THE MARKED LIST */
349                      PrevU->Next = CurU->Next;
350                      CurU->Next = u1Head.Next != NULL ? u1Head.Next : NULL;
351                      u1Head.Next = CurU;
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 */
365          PrevU = &u1Head;
366          for (CurU=u1Head.Next; CurU != NULL; CurU=CurU->Next)
367            {
368              i = CurU->i;
369              /* FIND THE VARIABLES IN ROWS i */
370              PrevV = &v0Head;
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 = _C[i][j] - CurU->val;
378                      /* ...AND ADD IT TO THE MARKED LIST */
379                      PrevV->Next = CurV->Next;
380                      CurV->Next = v1Head.Next != NULL ? v1Head.Next: NULL;
381                      v1Head.Next = CurV;
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 = _C[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-_X] = 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-_X])
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-_X] && 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-_X] = 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-_X]);
587             
588             if (NewX == NULL)
589               {
590                 IsUsed[*CurX-_X] = 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-_X] = 0;
601           *CurX = NewX;
602           IsUsed[NewX-_X] = 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[MAX_SIG_SIZE1][MAX_SIG_SIZE1];
630  node1_t Ur[MAX_SIG_SIZE1], Vr[MAX_SIG_SIZE1];
631  node1_t uHead, *CurU, *PrevU;
632  node1_t vHead, *CurV, *PrevV;
633  node1_t *PrevUMinI, *PrevVMinJ, *Remember;
634
635  /* INITIALIZE THE ROWS LIST (Ur), AND THE COLUMNS LIST (Vr) */
636  uHead.Next = CurU = Ur;
637  for (i=0; i < _n1; i++)
638    {
639      CurU->i = i;
640      CurU->val = -INFINITY;
641      CurU->Next = CurU+1;
642      CurU++;
643    }
644  (--CurU)->Next = NULL;
645 
646  vHead.Next = CurV = Vr;
647  for (j=0; j < _n2; j++)
648    {
649      CurV->i = j;
650      CurV->val = -INFINITY;
651      CurV->Next = CurV+1;
652      CurV++;
653    }
654  (--CurV)->Next = NULL;
655 
656  /* FIND THE MAXIMUM ROW AND COLUMN VALUES (Ur[i] AND Vr[j]) */
657  for(i=0; i < _n1 ; i++)
658    for(j=0; j < _n2 ; j++)
659      {
660        float v;
661        v = _C[i][j];
662        if (Ur[i].val <= v)
663          Ur[i].val = v;
664        if (Vr[j].val <= v)
665          Vr[j].val = v;
666      }
667 
668  /* COMPUTE THE Delta MATRIX */
669  for(i=0; i < _n1 ; i++)
670    for(j=0; j < _n2 ; j++)
671      Delta[i][j] = _C[i][j] - Ur[i].val - Vr[j].val;
672
673  /* FIND THE BASIC VARIABLES */
674  do
675    {
676#if DEBUG_LEVEL > 3
677      printf("Ur=");
678      for(CurU = uHead.Next; CurU != NULL; CurU = CurU->Next)
679        printf("[%d]",CurU-Ur);
680      printf("\n");
681      printf("Vr=");
682      for(CurV = vHead.Next; CurV != NULL; CurV = CurV->Next)
683        printf("[%d]",CurV-Vr);
684      printf("\n");
685      printf("\n\n");
686#endif
687 
688      /* FIND THE SMALLEST Delta[i][j] */
689      found = 0;
690      deltaMin = INFINITY;     
691      PrevU = &uHead;
692      for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
693        {
694          int i;
695          i = CurU->i;
696          PrevV = &vHead;
697          for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
698            {
699              int j;
700              j = CurV->i;
701              if (deltaMin > Delta[i][j])
702                {
703                  deltaMin = Delta[i][j];
704                  minI = i;
705                  minJ = j;
706                  PrevUMinI = PrevU;
707                  PrevVMinJ = PrevV;
708                  found = 1;
709                }
710              PrevV = CurV;
711            }
712          PrevU = CurU;
713        }
714     
715      if (! found)
716        break;
717
718      /* ADD X[minI][minJ] TO THE BASIS, AND ADJUST SUPPLIES AND COST */
719      Remember = PrevUMinI->Next;
720      addBasicVariable(minI, minJ, S, D, PrevUMinI, PrevVMinJ, &uHead);
721
722      /* UPDATE THE NECESSARY Delta[][] */
723      if (Remember == PrevUMinI->Next)  /* LINE minI WAS DELETED */
724        {
725          for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
726            {
727              int j;
728              j = CurV->i;
729              if (CurV->val == _C[minI][j])  /* COLUMN j NEEDS UPDATING */
730                {
731                  /* FIND THE NEW MAXIMUM VALUE IN THE COLUMN */
732                  oldVal = CurV->val;
733                  CurV->val = -INFINITY;
734                  for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
735                    {
736                      int i;
737                      i = CurU->i;
738                      if (CurV->val <= _C[i][j])
739                        CurV->val = _C[i][j];
740                    }
741                 
742                  /* IF NEEDED, ADJUST THE RELEVANT Delta[*][j] */
743                  diff = oldVal - CurV->val;
744                  if (fabs(diff) < EPSILON * _maxC)
745                    for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
746                      Delta[CurU->i][j] += diff;
747                }
748            }
749        }
750      else  /* COLUMN minJ WAS DELETED */
751        {
752          for (CurU=uHead.Next; CurU != NULL; CurU=CurU->Next)
753            {
754              int i;
755              i = CurU->i;
756              if (CurU->val == _C[i][minJ])  /* ROW i NEEDS UPDATING */
757                {
758                  /* FIND THE NEW MAXIMUM VALUE IN THE ROW */
759                  oldVal = CurU->val;
760                  CurU->val = -INFINITY;
761                  for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
762                    {
763                      int j;
764                      j = CurV->i;
765                      if(CurU->val <= _C[i][j])
766                        CurU->val = _C[i][j];
767                    }
768                 
769                  /* If NEEDED, ADJUST THE RELEVANT Delta[i][*] */
770                  diff = oldVal - CurU->val;
771                  if (fabs(diff) < EPSILON * _maxC)
772                    for (CurV=vHead.Next; CurV != NULL; CurV=CurV->Next)
773                      Delta[i][CurV->i] += diff;
774                }
775            }
776        }
777    } while (uHead.Next != NULL || vHead.Next != NULL);
778}
779
780
781
782
783/**********************
784    addBasicVariable
785**********************/
786static void addBasicVariable(int minI, int minJ, double *S, double *D,
787                             node1_t *PrevUMinI, node1_t *PrevVMinJ,
788                             node1_t *UHead)
789{
790  double T;
791 
792  if (fabs(S[minI]-D[minJ]) <= EPSILON * _maxW)  /* DEGENERATE CASE */
793    {
794      T = S[minI];
795      S[minI] = 0;
796      D[minJ] -= T;
797    }
798  else if (S[minI] < D[minJ])  /* SUPPLY EXHAUSTED */
799    {
800      T = S[minI];
801      S[minI] = 0;
802      D[minJ] -= T;
803    }
804  else  /* DEMAND EXHAUSTED */
805    {
806      T = D[minJ];
807      D[minJ] = 0;
808      S[minI] -= T;
809    }
810
811  /* X(minI,minJ) IS A BASIC VARIABLE */
812  _IsX[minI][minJ] = 1;
813
814  _EndX->val = T;
815  _EndX->i = minI;
816  _EndX->j = minJ;
817  _EndX->NextC = _RowsX[minI];
818  _EndX->NextR = _ColsX[minJ];
819  _RowsX[minI] = _EndX;
820  _ColsX[minJ] = _EndX;
821  _EndX++;
822
823  /* DELETE SUPPLY ROW ONLY IF THE EMPTY, AND IF NOT LAST ROW */
824  if (S[minI] == 0 && UHead->Next->Next != NULL)
825    PrevUMinI->Next = PrevUMinI->Next->Next;  /* REMOVE ROW FROM LIST */
826  else
827    PrevVMinJ->Next = PrevVMinJ->Next->Next;  /* REMOVE COLUMN FROM LIST */
828}
Note: See TracBrowser for help on using the repository browser.