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

Last change on this file since 1064 was 1064, checked in by Maciej Komosinski, 3 years ago

Properly allocated and de-allocated dynamic arrays of size calculated in runtime, not relying on g++ "extension"

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