1 |
greg |
2.2 |
/* RCSid $Id: transportSimplex.h,v 2.1 2014/03/26 02:52:31 greg Exp $ */ |
2 |
greg |
2.1 |
/* |
3 |
|
|
transportSimplex.h |
4 |
|
|
|
5 |
|
|
A C++ implementation of the transportation simplex algorithm. |
6 |
|
|
|
7 |
|
|
Last edit Sept 3 2006 |
8 |
|
|
|
9 |
|
|
Copyright (C) 2006 Darren MacDonald |
10 |
|
|
[email protected] |
11 |
|
|
www.site.uottawa.ca/~dmacd070/emd |
12 |
|
|
|
13 |
|
|
except some data structures and interface adapted from |
14 |
|
|
http://ai.stanford.edu/~rubner/emd/emd.c, |
15 |
|
|
Copyright 1998 Yossi Rubner. |
16 |
|
|
|
17 |
|
|
This algorithm solves the problem of finding the least-cost way to transport |
18 |
|
|
goods from a set of source (supply) nodes to a set of sink (demand) nodes. |
19 |
|
|
|
20 |
|
|
To use the code, simply include this file in the user code and add 'using namespace t_simplex;'. |
21 |
|
|
Organize the nodes into signatures using the TsSignature<TF> class, which contains a feature |
22 |
|
|
array (of type TF) and an array of the respective weights of the features. |
23 |
|
|
|
24 |
|
|
Solve the system by calling |
25 |
|
|
transportSimplex(&Sig1, &Sig2, grndDist); |
26 |
|
|
where Sig1 is the source signature, containing an array of source features and an array of their respective |
27 |
|
|
supplies, Sig2 is the sink signature, containing the sink features and and their respective demands, and |
28 |
|
|
grndDist is a pointer to a function which accepts two features pointers as arguments and returns the unit |
29 |
|
|
cost of transporting goods between the two. |
30 |
|
|
|
31 |
|
|
transportSimplex either returns the optimal transportation cost or throws a TsError. |
32 |
|
|
|
33 |
|
|
For more information see the documentation at |
34 |
|
|
www.site.uottawa.ca/~dmacd070/emd/index.html |
35 |
|
|
and the example implemenation at |
36 |
|
|
www.site.uottawa.ca/~dmacd070/emd/main.cpp. |
37 |
|
|
*/ |
38 |
|
|
|
39 |
|
|
#ifndef _T_SIMPLEX_H |
40 |
|
|
#define _T_SIMPLEX_H |
41 |
|
|
#include <iostream> |
42 |
|
|
#include <stdlib.h> |
43 |
|
|
#include <math.h> |
44 |
|
|
#include <new> |
45 |
|
|
|
46 |
|
|
#define TSINFINITY 1e20 |
47 |
|
|
#define TSEPSILON 1e-6 |
48 |
|
|
#define TSPIVOTLIMIT 0.00 |
49 |
|
|
|
50 |
|
|
namespace t_simplex { |
51 |
|
|
|
52 |
|
|
/* DECLARATION OF DATA TYPES */ |
53 |
|
|
enum TsError { TsErrBadInput }; |
54 |
|
|
|
55 |
|
|
//TsSignature is used for inputting the source and sink signatures |
56 |
|
|
template <class TF> |
57 |
|
|
class TsSignature { |
58 |
|
|
public: |
59 |
|
|
int n; // Number of features in the signature |
60 |
|
|
TF *features; // Pointer to the features vector |
61 |
|
|
double *weights; // Pointer to the weights of the features |
62 |
|
|
TsSignature(int nin, TF *fin, double * win):n(nin), features(fin), weights(win){}; |
63 |
|
|
}; |
64 |
|
|
|
65 |
|
|
//TsFlow is used for outputting the final flow table |
66 |
|
|
typedef struct TsFlow { |
67 |
|
|
int from; // Feature number in signature 1 |
68 |
|
|
int to; // Feature number in signature 2 |
69 |
|
|
double amount; // Amount of flow from signature1.features[from] to signature2.features[to] |
70 |
|
|
} TsFlow; |
71 |
|
|
|
72 |
|
|
// TsBasic is used for 2D lists, allowing for easy navigation of the basic variables |
73 |
|
|
typedef struct TsBasic{ |
74 |
|
|
int i, j; |
75 |
|
|
double val; |
76 |
|
|
TsBasic *nextSnk, *prevSnk; // next/previous node in the column |
77 |
|
|
TsBasic *nextSrc, *prevSrc; // next/previous node in the row |
78 |
|
|
} TsBasic; |
79 |
|
|
|
80 |
|
|
// TsStone is used for _BFS |
81 |
|
|
typedef struct TsStone { |
82 |
|
|
struct TsStone *prev; |
83 |
|
|
struct TsBasic *node; |
84 |
|
|
} TsStone; |
85 |
|
|
|
86 |
|
|
// TsRusPen is used for 1D lists in _initRussel |
87 |
|
|
typedef struct TsRusPen { |
88 |
|
|
int i; |
89 |
|
|
double val; |
90 |
|
|
struct TsRusPen *next, *prev; |
91 |
|
|
} TsRusPen; |
92 |
|
|
|
93 |
|
|
// TsVogPen is used for 1D lists in _initVogel |
94 |
|
|
typedef struct TsVogPen { |
95 |
|
|
int i; |
96 |
|
|
struct TsVogPen *next, *prev; |
97 |
|
|
int one, two; |
98 |
|
|
double oneCost, twoCost; |
99 |
|
|
} TsVogPen; |
100 |
|
|
|
101 |
|
|
/* DECLARATION OF GLOBALS */ |
102 |
|
|
double ** _tsC = NULL; // Cost matrix |
103 |
|
|
double _tsMaxC; // Maximum of all costs |
104 |
|
|
double _tsMaxW; // Maximum of all weights |
105 |
|
|
|
106 |
|
|
|
107 |
|
|
/* INTERNAL FUNCTIONS */ |
108 |
|
|
double _pivot(TsBasic * basics, TsBasic ** srcBasics, TsBasic ** snkBasics, bool ** isBasic, int n1, int n2); |
109 |
|
|
TsStone * _BFS(TsStone *stoneTree, TsBasic ** srcBasics, TsBasic ** snkBasics, bool complete = false); |
110 |
|
|
void _initVogel(double *S, double *D, TsBasic * basicsEnd, TsBasic ** srcBasics, TsBasic ** snkBasics, bool ** isBasic, int n1, int n2); |
111 |
|
|
void _initRussel(double *S, double *D, TsBasic * basicsEnd, TsBasic ** srcBasics, TsBasic ** snkBasics, bool ** isBasic, int n1, int n2); |
112 |
|
|
|
113 |
|
|
/* |
114 |
|
|
transportSimplex() - Program entry point. |
115 |
|
|
|
116 |
|
|
signature1 and signature2 define the source and sink sets. Unit costs between two features are computed from |
117 |
|
|
grndDist. Signature weights must be positive and all costs must be positive. |
118 |
|
|
TsFlow is an output parameter which can be set to an array that will be filled with the final flow amounts. |
119 |
|
|
The array must be of size signature1->n + signature2->n - 1 . flowSize is a pointer to an integer which |
120 |
|
|
indicates the number of functional entries in Flow, because all spaces are not necessarily used. Flow and |
121 |
|
|
flowSize can be set to NULL or omitted from the argument list if this information is not important. |
122 |
|
|
|
123 |
|
|
The return value is the transportation cost. |
124 |
|
|
|
125 |
|
|
*/ |
126 |
|
|
template <class TF> |
127 |
|
|
double transportSimplex(TsSignature<TF> *signature1, TsSignature<TF> *signature2, double (*grndDist)(TF *, TF *), |
128 |
|
|
TsFlow *flowTable = NULL, int *flowSize = NULL) { |
129 |
|
|
|
130 |
|
|
int n1, n2; // number of features in signature1 and Signature 2 |
131 |
|
|
int i, j; |
132 |
|
|
|
133 |
|
|
double totalCost; |
134 |
|
|
double w; |
135 |
|
|
|
136 |
|
|
double srcSum, snkSum, diff; |
137 |
|
|
|
138 |
|
|
TF *P1, *P2; |
139 |
|
|
|
140 |
|
|
n1 = signature1->n; |
141 |
|
|
n2 = signature2->n; |
142 |
|
|
|
143 |
|
|
TsBasic * basics = NULL; ///Array of basic variables. |
144 |
|
|
bool ** isBasic = NULL; //Flag matrix. isBasic[i][j] is true there is flow between source i and sink j |
145 |
|
|
TsBasic **srcBasics = NULL; //Array of pointers to the first basic variable in each row |
146 |
|
|
TsBasic **snkBasics = NULL; //Array of pointers to the first basic variable in each column |
147 |
|
|
double * src = NULL; //Array of source supplies |
148 |
|
|
double * snk =NULL; //Array of sink demands |
149 |
|
|
|
150 |
|
|
|
151 |
|
|
// Equalize source and sink weights. A dummy source or sink may be added to equalize the total sink |
152 |
|
|
// and source weights. n1 = signature1->n + 1 if there is a dummy source, and n2 = signature2->n + 1 |
153 |
|
|
// if there is a dummy sink. |
154 |
|
|
srcSum = 0.0; |
155 |
|
|
for(i=0; i < n1; i++) |
156 |
|
|
srcSum += signature1->weights[i]; |
157 |
|
|
|
158 |
|
|
snkSum = 0.0; |
159 |
|
|
for(j=0; j < n2; j++) |
160 |
|
|
snkSum += signature2->weights[j]; |
161 |
|
|
|
162 |
|
|
diff = srcSum - snkSum; |
163 |
greg |
2.2 |
if (fabs(diff) > TSEPSILON * srcSum) { |
164 |
greg |
2.1 |
if (diff < 0.0) |
165 |
|
|
n1++; |
166 |
|
|
else |
167 |
|
|
n2++; |
168 |
greg |
2.2 |
} |
169 |
greg |
2.1 |
_tsMaxW = srcSum > snkSum ? srcSum : snkSum; |
170 |
|
|
w = srcSum < snkSum ? srcSum : snkSum; |
171 |
|
|
|
172 |
|
|
//Allocate memory for arrays |
173 |
|
|
try { |
174 |
|
|
|
175 |
|
|
basics = new TsBasic[n1 + n2]; |
176 |
|
|
|
177 |
|
|
isBasic = new bool*[n1]; |
178 |
|
|
for(i = 0; i < n1; i++) |
179 |
|
|
isBasic[i] = NULL; |
180 |
|
|
|
181 |
|
|
for(i = 0; i < n1; i++) { |
182 |
|
|
isBasic[i] = new bool[n2]; |
183 |
|
|
for (j=0; j < n2; j++) |
184 |
|
|
isBasic[i][j] = 0; |
185 |
|
|
} |
186 |
|
|
|
187 |
|
|
srcBasics = new TsBasic*[n1]; |
188 |
|
|
for (i = 0; i < n1; i++) |
189 |
|
|
srcBasics[i] = NULL; |
190 |
|
|
|
191 |
|
|
snkBasics = new TsBasic*[n2]; |
192 |
|
|
for (i = 0; i < n2; i++) |
193 |
|
|
snkBasics[i] = NULL; |
194 |
|
|
|
195 |
|
|
// Compute the cost matrix |
196 |
|
|
_tsMaxC = 0; |
197 |
|
|
_tsC = new double*[n1]; |
198 |
|
|
for(i=0; i < n1; i++) |
199 |
|
|
_tsC[i] = NULL; |
200 |
|
|
|
201 |
|
|
for(i=0, P1=signature1->features; i < n1; i++, P1++) { |
202 |
|
|
_tsC[i] = new double[n2]; //What happens if bad here? |
203 |
|
|
for(j=0, P2=signature2->features; j < n2; j++, P2++) { |
204 |
|
|
if (i == signature1->n || j == signature2->n) { |
205 |
|
|
_tsC[i][j] = 0; // cost is zero for flow to or from a dummy |
206 |
|
|
} else { |
207 |
|
|
_tsC[i][j] = grndDist(P1, P2); |
208 |
|
|
if (_tsC[i][j] < 0) throw TsErrBadInput; |
209 |
|
|
} |
210 |
|
|
|
211 |
|
|
|
212 |
|
|
if (_tsC[i][j] > _tsMaxC) |
213 |
|
|
_tsMaxC = _tsC[i][j]; |
214 |
|
|
} |
215 |
|
|
} |
216 |
|
|
|
217 |
|
|
src = new double[n1]; //init the source array |
218 |
|
|
for (i = 0; i < signature1->n; i++) { |
219 |
|
|
src[i] = signature1->weights[i]; |
220 |
|
|
if (src[i] < 0) throw TsErrBadInput; |
221 |
|
|
} |
222 |
|
|
if (n1 != signature1->n) |
223 |
|
|
src[signature1->n] = -diff; |
224 |
|
|
|
225 |
|
|
snk = new double[n2]; //init the sink array |
226 |
|
|
for (i = 0; i < signature2->n; i++) { |
227 |
|
|
snk[i] = signature2->weights[i]; |
228 |
|
|
if (snk[i] < 0) throw TsErrBadInput; |
229 |
|
|
} |
230 |
|
|
if (n2 != signature2->n) |
231 |
|
|
snk[signature2->n] = diff; |
232 |
|
|
|
233 |
|
|
// Find the initail basic feasible solution. Use either _initRussel or _initVogel |
234 |
|
|
|
235 |
|
|
|
236 |
|
|
_initRussel(src, snk, basics, srcBasics, snkBasics, isBasic, n1, n2); |
237 |
|
|
//_initVogel(src, snk, basics, srcBasics, snkBasics, isBasic, n1, n2); |
238 |
|
|
|
239 |
|
|
|
240 |
|
|
// Enter the main pivot loop |
241 |
|
|
totalCost = _pivot(basics, srcBasics, snkBasics, isBasic, n1, n2); |
242 |
|
|
|
243 |
|
|
} catch (...) { |
244 |
|
|
for(i = 0; i < n1; i++) |
245 |
|
|
delete[] isBasic[i]; |
246 |
|
|
delete[] isBasic; |
247 |
|
|
|
248 |
|
|
for (i = 0; i < n1; i++) |
249 |
|
|
delete[] _tsC[i]; |
250 |
|
|
delete[] _tsC; |
251 |
|
|
|
252 |
|
|
delete[] src; |
253 |
|
|
delete[] snk; |
254 |
|
|
delete[] srcBasics; |
255 |
|
|
delete[] snkBasics; |
256 |
|
|
delete[] basics; |
257 |
|
|
|
258 |
|
|
throw; |
259 |
|
|
} |
260 |
|
|
|
261 |
|
|
// Fill the Flow data structure |
262 |
|
|
TsBasic * basicPtr; |
263 |
|
|
TsFlow * flowPtr = flowTable; |
264 |
|
|
if (flowTable != NULL) { |
265 |
|
|
for (i = 0; i < n1+n2 -1; i++) { |
266 |
|
|
basicPtr = basics + i; |
267 |
|
|
if (isBasic[basicPtr->i][basicPtr->j] && basicPtr->i != signature1->n && basicPtr->j != signature2->n && basicPtr->val != 0.0) { |
268 |
|
|
flowPtr->to = basicPtr->j; |
269 |
|
|
flowPtr->from = basicPtr->i; |
270 |
|
|
flowPtr->amount = basicPtr->val; |
271 |
|
|
flowPtr++; |
272 |
|
|
} |
273 |
|
|
} |
274 |
|
|
} |
275 |
|
|
if (flowSize != NULL) { |
276 |
|
|
*flowSize = (int)(flowPtr - flowTable); |
277 |
|
|
} |
278 |
|
|
|
279 |
|
|
|
280 |
|
|
for(i = 0; i < n1; i++) |
281 |
|
|
delete[] isBasic[i]; |
282 |
|
|
delete[] isBasic; |
283 |
|
|
|
284 |
|
|
for (i = 0; i < n1; i++) |
285 |
|
|
delete[] _tsC[i]; |
286 |
|
|
delete[] _tsC; |
287 |
|
|
|
288 |
|
|
delete[] src; |
289 |
|
|
delete[] snk; |
290 |
|
|
delete[] srcBasics; |
291 |
|
|
delete[] snkBasics; |
292 |
|
|
delete[] basics; |
293 |
|
|
|
294 |
|
|
return totalCost; |
295 |
|
|
} |
296 |
|
|
|
297 |
|
|
/* |
298 |
|
|
Main pivot loop. |
299 |
|
|
Pivots until the system is optimal and return the optimal transportation cost. |
300 |
|
|
*/ |
301 |
|
|
double _pivot(TsBasic * basics, TsBasic ** srcBasics, TsBasic ** snkBasics, bool ** isBasic, int n1, int n2) { |
302 |
|
|
|
303 |
|
|
double * srcDuals = NULL; |
304 |
|
|
double * snkDuals = NULL; |
305 |
|
|
TsStone * stonePath = NULL; |
306 |
|
|
|
307 |
|
|
TsStone * spitra, * spitrb, * leaving; |
308 |
|
|
TsBasic * XP; |
309 |
|
|
TsBasic * basicsEnd = basics + n1 + n2; |
310 |
|
|
TsBasic * entering = basicsEnd - 1 ; |
311 |
|
|
TsBasic dummyBasic; |
312 |
|
|
dummyBasic.i = -1; |
313 |
|
|
dummyBasic.j = 0; |
314 |
|
|
|
315 |
|
|
int i,j, lowI, lowJ; |
316 |
|
|
double objectiveValue = TSINFINITY, oldObjectiveValue = 0; |
317 |
|
|
double lowVal; |
318 |
|
|
int numPivots = 0; |
319 |
|
|
|
320 |
|
|
try { |
321 |
|
|
srcDuals = new double[n1]; |
322 |
|
|
snkDuals = new double[n2]; |
323 |
|
|
stonePath = new TsStone[n1 + n2]; |
324 |
|
|
} catch (std::bad_alloc) { |
325 |
|
|
delete[] srcDuals; |
326 |
|
|
delete[] snkDuals; |
327 |
|
|
delete[] stonePath; |
328 |
|
|
throw; |
329 |
|
|
} |
330 |
|
|
|
331 |
|
|
while (1) { |
332 |
|
|
|
333 |
|
|
oldObjectiveValue = objectiveValue; |
334 |
|
|
objectiveValue = 0; |
335 |
|
|
|
336 |
|
|
for (XP = basics; XP != basicsEnd; XP++){ |
337 |
|
|
if (XP != entering) |
338 |
|
|
objectiveValue += _tsC[XP->i][XP->j] * XP->val; |
339 |
|
|
} |
340 |
|
|
|
341 |
|
|
// Compute the dual variables for each row and column. Begin by finding a spanning tree (stonepath) |
342 |
|
|
// of the basic variables using the breadth-first search routine seeded at an imaginary basic |
343 |
|
|
// variable in the first column. The dual variables can then be computed incrementally by traversing |
344 |
|
|
// the tree. |
345 |
|
|
stonePath[0].node = &dummyBasic; |
346 |
|
|
stonePath[0].prev = NULL; |
347 |
|
|
spitrb = _BFS(stonePath, srcBasics, snkBasics, true); |
348 |
|
|
|
349 |
|
|
spitra = stonePath; |
350 |
|
|
snkDuals[spitra->node->j] = 0; |
351 |
|
|
for (spitra++; spitra != spitrb; spitra++) { |
352 |
|
|
if (spitra->node->i == spitra->prev->node->i) { |
353 |
|
|
//node is in same row as parent |
354 |
|
|
snkDuals[spitra->node->j] = _tsC[spitra->node->i][spitra->node->j] - srcDuals[spitra->node->i]; |
355 |
|
|
} else if (spitra->node->j == spitra->prev->node->j) { |
356 |
|
|
srcDuals[spitra->node->i] = _tsC[spitra->node->i][spitra->node->j] - snkDuals[spitra->node->j]; |
357 |
|
|
} |
358 |
|
|
} |
359 |
|
|
|
360 |
|
|
// After computing the duals, find the non-basic variable that has the greatest negative value of |
361 |
|
|
// delta = _tsC[i][j] - srcDuals[i] - snkDuals[j]. This is the entering variable |
362 |
|
|
lowVal = 0.0; |
363 |
|
|
for (i = 0; i < n1; i++) |
364 |
|
|
for (j = 0; j < n2; j++) |
365 |
|
|
if (!isBasic[i][j] && _tsC[i][j] - srcDuals[i] - snkDuals[j] < lowVal) { |
366 |
|
|
lowVal = _tsC[i][j] - srcDuals[i] - snkDuals[j]; |
367 |
|
|
lowI = i; |
368 |
|
|
lowJ = j; |
369 |
|
|
} |
370 |
|
|
|
371 |
|
|
// If all delta values are non-negative, the table is optimal |
372 |
|
|
if (lowVal >= -TSEPSILON * _tsMaxC || (oldObjectiveValue - objectiveValue) < TSPIVOTLIMIT) { |
373 |
|
|
delete[] srcDuals; |
374 |
|
|
delete[] snkDuals; |
375 |
|
|
delete[] stonePath; |
376 |
|
|
//std::cout << numPivots << "\t"; |
377 |
|
|
return objectiveValue; |
378 |
|
|
} |
379 |
|
|
|
380 |
|
|
// Add the entering variable |
381 |
|
|
entering->i = lowI; |
382 |
|
|
entering->j = lowJ; |
383 |
|
|
isBasic[lowI][lowJ] = 1; |
384 |
|
|
entering->val = 0; |
385 |
|
|
|
386 |
|
|
entering->nextSrc = srcBasics[lowI]; |
387 |
|
|
if (srcBasics[lowI] != NULL) srcBasics[lowI]->prevSrc = entering; |
388 |
|
|
entering->nextSnk = snkBasics[lowJ]; |
389 |
|
|
if (snkBasics[lowJ] != NULL) snkBasics[lowJ]->prevSnk = entering; |
390 |
|
|
|
391 |
|
|
srcBasics[lowI] = entering; |
392 |
|
|
entering->prevSrc = srcBasics[lowI]; |
393 |
|
|
snkBasics[lowJ] = entering; |
394 |
|
|
entering->prevSnk = snkBasics[lowJ]; |
395 |
|
|
|
396 |
|
|
stonePath[0].node = entering; |
397 |
|
|
stonePath[0].prev = NULL; |
398 |
|
|
|
399 |
|
|
// Use breadth-first search to find a loop of basics. |
400 |
|
|
spitra = spitrb = _BFS(stonePath, srcBasics, snkBasics); |
401 |
|
|
lowVal = TSINFINITY; |
402 |
|
|
bool add = false; |
403 |
|
|
|
404 |
|
|
// Find the lowest flow along the loop (leaving variable) |
405 |
|
|
do { |
406 |
|
|
if (!add && spitrb->node->val < lowVal) { |
407 |
|
|
leaving = spitrb; |
408 |
|
|
lowVal = spitrb->node->val; |
409 |
|
|
} |
410 |
|
|
add = !add; |
411 |
greg |
2.2 |
} while ((spitrb = spitrb->prev)); |
412 |
greg |
2.1 |
|
413 |
|
|
add = false; |
414 |
|
|
spitrb = spitra; |
415 |
|
|
|
416 |
|
|
// Alternately increase and decrease flow along the loop |
417 |
|
|
do { |
418 |
|
|
if (add) spitrb->node->val += lowVal; |
419 |
|
|
else spitrb->node->val -= lowVal; |
420 |
|
|
add = !add; |
421 |
greg |
2.2 |
} while ((spitrb = spitrb->prev)); |
422 |
greg |
2.1 |
|
423 |
|
|
i = leaving->node->i; |
424 |
|
|
j = leaving->node->j; |
425 |
|
|
isBasic[i][j] = 0; |
426 |
|
|
|
427 |
|
|
|
428 |
|
|
if (srcBasics[i] == leaving->node) { |
429 |
|
|
srcBasics[i] = leaving->node->nextSrc; |
430 |
|
|
srcBasics[i]->prevSrc = NULL; |
431 |
|
|
} else { |
432 |
|
|
leaving->node->prevSrc->nextSrc = leaving->node->nextSrc; |
433 |
|
|
if (leaving->node->nextSrc != NULL) |
434 |
|
|
leaving->node->nextSrc->prevSrc = leaving->node->prevSrc; |
435 |
|
|
} |
436 |
|
|
|
437 |
|
|
if (snkBasics[j] == leaving->node) { |
438 |
|
|
snkBasics[j] = leaving->node->nextSnk; |
439 |
|
|
snkBasics[j]->prevSnk = NULL; |
440 |
|
|
} else { |
441 |
|
|
leaving->node->prevSnk->nextSnk = leaving->node->nextSnk; |
442 |
|
|
if (leaving->node->nextSnk != NULL) |
443 |
|
|
leaving->node->nextSnk->prevSnk = leaving->node->prevSnk; |
444 |
|
|
} |
445 |
|
|
entering = leaving->node; |
446 |
|
|
numPivots++; |
447 |
|
|
} |
448 |
|
|
}; |
449 |
|
|
|
450 |
|
|
|
451 |
|
|
/******************* |
452 |
|
|
_BFS |
453 |
|
|
Perform a breadth-first search of the basic varaibles. The search tree is stored in stoneTree, where each |
454 |
|
|
'stone' contains a pointer to a basic variable, and a pointer to that stone's parent. srcBasics and snkBasics |
455 |
|
|
are arrays of linked lists which allow for easy identification of a node's neighbours in the flow network. |
456 |
|
|
If complete == true, find all basics in the table. Otherwise, terminate when a basic variable which completes |
457 |
|
|
a loop has been found and return a pointer to the final stone in the loop. |
458 |
|
|
*********************/ |
459 |
|
|
TsStone * _BFS(TsStone * stoneTree, TsBasic ** srcBasics, TsBasic ** snkBasics, bool complete) { |
460 |
|
|
bool column = true; |
461 |
|
|
int jumpoffset = 0; |
462 |
|
|
TsBasic * bitr; |
463 |
|
|
TsStone * sitra = &stoneTree[0], * sitrb = &stoneTree[1]; |
464 |
|
|
do { |
465 |
|
|
if (column) { |
466 |
|
|
for (bitr = snkBasics[sitra->node->j]; bitr != NULL; bitr = bitr->nextSnk) { |
467 |
|
|
if (bitr != sitra->node){ |
468 |
|
|
sitrb->node = bitr; |
469 |
|
|
sitrb->prev = sitra; |
470 |
|
|
sitrb++; |
471 |
|
|
} |
472 |
|
|
} |
473 |
|
|
} else { |
474 |
|
|
for (bitr = srcBasics[sitra->node->i]; bitr != NULL; bitr = bitr->nextSrc) { |
475 |
|
|
if (bitr != sitra->node){ |
476 |
|
|
sitrb->node = bitr; |
477 |
|
|
sitrb->prev = sitra; |
478 |
|
|
sitrb++; |
479 |
|
|
} |
480 |
|
|
} |
481 |
|
|
} |
482 |
|
|
|
483 |
|
|
sitra++; |
484 |
|
|
if (sitra == sitrb) //no cycle found and no cycles in tree |
485 |
|
|
return sitra; |
486 |
|
|
|
487 |
|
|
if (sitra->node->i == sitra->prev->node->i) |
488 |
|
|
column = true; |
489 |
|
|
else |
490 |
|
|
column = false; |
491 |
|
|
|
492 |
|
|
// cycle found |
493 |
|
|
if (!complete && sitra->node->i == stoneTree[0].node->i && sitra->node->j != stoneTree[0].node->j && column == false) |
494 |
|
|
return sitra; |
495 |
|
|
} while(1); |
496 |
|
|
} |
497 |
|
|
|
498 |
|
|
|
499 |
|
|
// Helper function for _initVogel |
500 |
|
|
inline void addPenalty(TsVogPen * pitr, double cost, int i) { |
501 |
|
|
if (cost < pitr->oneCost) { |
502 |
|
|
pitr->twoCost = pitr->oneCost; |
503 |
|
|
pitr->two = pitr->one; |
504 |
|
|
pitr->oneCost = cost; |
505 |
|
|
pitr->one = i; |
506 |
|
|
} else if (cost < pitr->twoCost) { |
507 |
|
|
pitr->twoCost = cost; |
508 |
|
|
pitr->two = i; |
509 |
|
|
} |
510 |
|
|
} |
511 |
|
|
|
512 |
|
|
|
513 |
|
|
/********************** |
514 |
|
|
Vogel's initialization method |
515 |
|
|
**********************/ |
516 |
|
|
void _initVogel(double *S, double *D, TsBasic * basicsEnd, TsBasic ** srcBasics, TsBasic ** snkBasics, bool ** isBasic, int n1, int n2) { |
517 |
|
|
int i, j; |
518 |
|
|
TsVogPen *srcPens = NULL; |
519 |
|
|
TsVogPen *snkPens = NULL; |
520 |
|
|
TsVogPen *pitra, *pitrb; //iterators |
521 |
|
|
TsVogPen *maxPen; |
522 |
|
|
TsVogPen srcPenHead, snkPenHead; |
523 |
|
|
bool maxIsSrc; |
524 |
|
|
double lowVal; |
525 |
|
|
|
526 |
|
|
try { |
527 |
|
|
srcPens = new TsVogPen[n1]; |
528 |
|
|
snkPens = new TsVogPen[n2]; |
529 |
|
|
} catch (std::bad_alloc) { |
530 |
|
|
delete[] srcPens; |
531 |
|
|
delete[] snkPens; |
532 |
|
|
throw; |
533 |
|
|
} |
534 |
|
|
|
535 |
|
|
srcPenHead.next = pitra = srcPens; |
536 |
|
|
for (i=0; i < n1; i++) { |
537 |
|
|
pitra->i = i; |
538 |
|
|
pitra->next = pitra+1; |
539 |
|
|
pitra->prev = pitra-1; |
540 |
|
|
pitra->one = pitra->two = 0; |
541 |
|
|
pitra->oneCost = pitra->twoCost = TSINFINITY; |
542 |
|
|
pitra++; |
543 |
|
|
} |
544 |
|
|
(--pitra)->next = NULL; |
545 |
|
|
srcPens[0].prev = &srcPenHead; |
546 |
|
|
|
547 |
|
|
snkPenHead.next = pitra = snkPens; |
548 |
|
|
for (i=0; i < n2; i++) { |
549 |
|
|
pitra->i = i; |
550 |
|
|
pitra->next = pitra+1; |
551 |
|
|
pitra->prev = pitra-1; |
552 |
|
|
pitra->one = pitra->two = 0; |
553 |
|
|
pitra->oneCost = pitra->twoCost = TSINFINITY; |
554 |
|
|
pitra++; |
555 |
|
|
} |
556 |
|
|
(--pitra)->next = NULL; |
557 |
|
|
snkPens[0].prev = &snkPenHead; |
558 |
|
|
|
559 |
|
|
|
560 |
|
|
for (pitra = srcPenHead.next, i=0; pitra != NULL; pitra = pitra->next, i++) |
561 |
|
|
for (pitrb = snkPenHead.next, j=0; pitrb != NULL; pitrb = pitrb->next, j++) { |
562 |
|
|
//initialize Source Penalties; |
563 |
|
|
addPenalty(pitra, _tsC[i][j], j); |
564 |
|
|
addPenalty(pitrb, _tsC[i][j], i); |
565 |
|
|
} |
566 |
|
|
|
567 |
|
|
|
568 |
|
|
while (srcPenHead.next != NULL && snkPenHead.next != NULL) { |
569 |
|
|
maxIsSrc = true; |
570 |
|
|
for (maxPen = pitra = srcPenHead.next; pitra != NULL; pitra = pitra->next) |
571 |
|
|
if ((pitra->twoCost - pitra->oneCost) > (maxPen->twoCost - maxPen->oneCost)) |
572 |
|
|
maxPen = pitra; |
573 |
|
|
|
574 |
|
|
for (pitra = snkPenHead.next; pitra != NULL; pitra = pitra->next) |
575 |
|
|
if ((pitra->twoCost - pitra->oneCost) > (maxPen->twoCost - maxPen->oneCost)) { |
576 |
|
|
maxPen = pitra; |
577 |
|
|
maxIsSrc = false; |
578 |
|
|
} |
579 |
|
|
|
580 |
|
|
if (maxIsSrc) { |
581 |
|
|
i = maxPen->i; |
582 |
|
|
j = maxPen->one; |
583 |
|
|
} else { |
584 |
|
|
j = maxPen->i; |
585 |
|
|
i = maxPen->one; |
586 |
|
|
} |
587 |
|
|
|
588 |
|
|
if (D[j] - S[i] > _tsMaxW * TSEPSILON || (srcPenHead.next->next != NULL && fabs(S[i] - D[j]) < _tsMaxW * TSEPSILON)) { |
589 |
|
|
//delete source |
590 |
|
|
lowVal = S[i]; |
591 |
|
|
maxPen = srcPens + i; |
592 |
|
|
maxPen->prev->next = maxPen->next; |
593 |
|
|
if (maxPen->next != NULL) |
594 |
|
|
maxPen->next->prev = maxPen->prev; |
595 |
|
|
|
596 |
|
|
for (pitra = snkPenHead.next; pitra != NULL; pitra = pitra->next) { |
597 |
|
|
if (pitra->one == i || pitra->two == i){ |
598 |
|
|
pitra->oneCost = TSINFINITY; |
599 |
|
|
pitra->twoCost = TSINFINITY; |
600 |
|
|
for (pitrb = srcPenHead.next; pitrb != NULL; pitrb = pitrb->next) |
601 |
|
|
addPenalty(pitra, _tsC[pitrb->i][pitra->i], pitrb->i); |
602 |
|
|
} |
603 |
|
|
} |
604 |
|
|
} else { |
605 |
|
|
//delete sink |
606 |
|
|
lowVal = D[j]; |
607 |
|
|
maxPen = snkPens + j; |
608 |
|
|
maxPen->prev->next = maxPen->next; |
609 |
|
|
if (maxPen->next != NULL) |
610 |
|
|
maxPen->next->prev = maxPen->prev; |
611 |
|
|
|
612 |
|
|
for (pitra = srcPenHead.next; pitra != NULL; pitra = pitra->next) { |
613 |
|
|
if (pitra->one == j || pitra->two == j){ |
614 |
|
|
pitra->oneCost = TSINFINITY; |
615 |
|
|
pitra->twoCost = TSINFINITY; |
616 |
|
|
for (pitrb = snkPenHead.next; pitrb != NULL; pitrb = pitrb->next) |
617 |
|
|
addPenalty(pitra, _tsC[pitra->i][pitrb->i], pitrb->i); |
618 |
|
|
} |
619 |
|
|
} |
620 |
|
|
} |
621 |
|
|
|
622 |
|
|
S[i] -= lowVal; |
623 |
|
|
D[j] -= lowVal; |
624 |
|
|
|
625 |
|
|
isBasic[i][j] = 1; |
626 |
|
|
basicsEnd->val = lowVal; |
627 |
|
|
basicsEnd->i = i; |
628 |
|
|
basicsEnd->j = j; |
629 |
|
|
|
630 |
|
|
basicsEnd->nextSnk = snkBasics[j]; |
631 |
|
|
if (snkBasics[j] != NULL) snkBasics[j]->prevSnk = basicsEnd; |
632 |
|
|
basicsEnd->nextSrc = srcBasics[i]; |
633 |
|
|
if (srcBasics[i] != NULL) srcBasics[i]->prevSrc = basicsEnd; |
634 |
|
|
|
635 |
|
|
srcBasics[i] = basicsEnd; |
636 |
|
|
basicsEnd->prevSnk = NULL; |
637 |
|
|
snkBasics[j] = basicsEnd; |
638 |
|
|
basicsEnd->prevSrc = NULL; |
639 |
|
|
|
640 |
|
|
basicsEnd++; |
641 |
|
|
|
642 |
|
|
} |
643 |
|
|
delete[] srcPens; |
644 |
|
|
delete[] snkPens; |
645 |
|
|
} |
646 |
|
|
|
647 |
|
|
|
648 |
|
|
/********************** |
649 |
|
|
Russel's initialization method |
650 |
|
|
**********************/ |
651 |
|
|
void _initRussel(double *S, double *D, TsBasic * basicsEnd, TsBasic ** srcBasics, TsBasic ** snkBasics, bool ** isBasic, int n1, int n2) { |
652 |
|
|
double ** Delta = NULL; |
653 |
|
|
int i, j, lowI, lowJ; |
654 |
|
|
TsRusPen *U = NULL; |
655 |
|
|
TsRusPen *V = NULL; |
656 |
|
|
TsRusPen *Uhead, *Vhead; |
657 |
|
|
TsRusPen *Uitr, *Vitr; |
658 |
|
|
double cost, lowVal; |
659 |
|
|
|
660 |
|
|
try { |
661 |
|
|
Delta = new double*[n1]; |
662 |
|
|
for (i = 0; i < n1; i++) |
663 |
|
|
Delta[i] = new double[n2]; |
664 |
|
|
|
665 |
|
|
U = new TsRusPen[n1]; |
666 |
|
|
V = new TsRusPen[n2]; |
667 |
|
|
} catch (std::bad_alloc) { |
668 |
|
|
for (i = 0; i < n1; i++) |
669 |
|
|
delete[] Delta[i]; |
670 |
|
|
delete[] Delta; |
671 |
|
|
|
672 |
|
|
delete[] U; |
673 |
|
|
delete[] V; |
674 |
|
|
throw; |
675 |
|
|
} |
676 |
|
|
|
677 |
|
|
for (i = 0; i < n1; i++) { |
678 |
|
|
U[i].i = i; |
679 |
|
|
U[i].val = 0; |
680 |
|
|
U[i].next = &U[i+1]; |
681 |
|
|
U[i].prev = &U[i-1]; |
682 |
|
|
} |
683 |
|
|
U[n1-1].next = NULL; |
684 |
|
|
U[0].prev = NULL; |
685 |
|
|
|
686 |
|
|
for (i = 0; i < n2; i++) { |
687 |
|
|
V[i].i = i; |
688 |
|
|
V[i].val = 0; |
689 |
|
|
V[i].next = &V[i+1]; |
690 |
|
|
V[i].prev = &V[i-1]; |
691 |
|
|
} |
692 |
|
|
V[n2-1].next = NULL; |
693 |
|
|
V[0].prev = NULL; |
694 |
|
|
|
695 |
|
|
for (i = 0; i < n1; i++) |
696 |
|
|
for (j = 0; j < n2; j++) { |
697 |
|
|
cost = _tsC[i][j]; |
698 |
|
|
if (cost > U[i].val) |
699 |
|
|
U[i].val = cost; |
700 |
|
|
if (cost > V[j].val) |
701 |
|
|
V[j].val = cost; |
702 |
|
|
} |
703 |
|
|
|
704 |
|
|
for (i = 0; i < n1; i++) |
705 |
|
|
for (j = 0; j < n2; j++) |
706 |
|
|
Delta[i][j] = _tsC[i][j] - U[i].val - V[j].val; |
707 |
|
|
|
708 |
|
|
Uhead = U; |
709 |
|
|
Vhead = V; |
710 |
|
|
while (Uhead != NULL && Vhead != NULL) { |
711 |
|
|
|
712 |
|
|
//Find lowest Delta |
713 |
|
|
lowVal = TSINFINITY; |
714 |
|
|
for (Uitr = Uhead; Uitr != NULL; Uitr = Uitr->next) |
715 |
|
|
for (Vitr = Vhead; Vitr != NULL; Vitr = Vitr->next) |
716 |
|
|
if (Delta[Uitr->i][Vitr->i] < lowVal) { |
717 |
|
|
lowI = Uitr->i; |
718 |
|
|
lowJ = Vitr->i; |
719 |
|
|
lowVal = Delta[Uitr->i][Vitr->i]; |
720 |
|
|
} |
721 |
|
|
|
722 |
|
|
|
723 |
|
|
if (D[lowJ] - S[lowI] > _tsMaxW * TSEPSILON || (fabs(S[lowI] - D[lowJ]) < _tsMaxW * TSEPSILON && Uhead->next != NULL)) { |
724 |
|
|
//Delete Source |
725 |
|
|
if (&U[lowI] == Uhead) { //Entering variable is first in list |
726 |
|
|
Uhead = Uhead->next; |
727 |
|
|
if (Uhead != NULL) Uhead->prev = NULL; |
728 |
|
|
} else { |
729 |
|
|
U[lowI].prev->next = U[lowI].next; //Entering variable is in middle of list; |
730 |
|
|
if (U[lowI].next != NULL) //Entering variable is at the end of the list; |
731 |
|
|
U[lowI].next->prev = U[lowI].prev; |
732 |
|
|
} |
733 |
|
|
//See if this source was the maximum cost for any dest |
734 |
|
|
for (Vitr = Vhead; Vitr != NULL; Vitr = Vitr->next) { |
735 |
|
|
if (Vitr->val == _tsC[lowI][Vitr->i]) { |
736 |
|
|
//it is; update the dest |
737 |
|
|
//find maximum cost in the dest |
738 |
|
|
Vitr->val = 0; |
739 |
|
|
for (Uitr = Uhead; Uitr != NULL; Uitr = Uitr->next) |
740 |
|
|
if (_tsC[Uitr->i][Vitr->i] > Vitr->val) |
741 |
|
|
Vitr->val = _tsC[Uitr->i][Vitr->i]; |
742 |
|
|
//update Delta |
743 |
|
|
for (Uitr = Uhead; Uitr != NULL; Uitr = Uitr->next) |
744 |
|
|
Delta[Uitr->i][Vitr->i] = _tsC[Uitr->i][Vitr->i] - Uitr->val - Vitr->val; |
745 |
|
|
|
746 |
|
|
} |
747 |
|
|
} |
748 |
|
|
lowVal = S[lowI]; |
749 |
|
|
|
750 |
|
|
} else { |
751 |
|
|
//Delete Dest |
752 |
|
|
if (&V[lowJ] == Vhead) { //Entering variable is first in list |
753 |
|
|
Vhead = Vhead->next; |
754 |
|
|
if (Vhead != NULL) Vhead->prev = NULL; |
755 |
|
|
} else { |
756 |
|
|
V[lowJ].prev->next = V[lowJ].next; //Entering variable is in middle of list; |
757 |
|
|
if (V[lowJ].next != NULL) //Entering variable is at the end of the list; |
758 |
|
|
V[lowJ].next->prev = V[lowJ].prev; |
759 |
|
|
} |
760 |
|
|
//See if this source was the maximum cost for any dest |
761 |
|
|
for (Uitr = Uhead; Uitr != NULL; Uitr = Uitr->next) { |
762 |
|
|
if (Uitr->val == _tsC[Uitr->i][lowJ]) { |
763 |
|
|
//it is; update the dest |
764 |
|
|
//find maximum cost in the dest |
765 |
|
|
Uitr->val = 0; |
766 |
|
|
for (Vitr = Vhead; Vitr != NULL; Vitr = Vitr->next) |
767 |
|
|
if (_tsC[Uitr->i][Vitr->i] > Uitr->val) |
768 |
|
|
Uitr->val = _tsC[Uitr->i][Vitr->i]; |
769 |
|
|
//update Delta |
770 |
|
|
for (Vitr = Vhead; Vitr != NULL; Vitr = Vitr->next) |
771 |
|
|
Delta[Uitr->i][Vitr->i] = _tsC[Uitr->i][Vitr->i] - Uitr->val - Vitr->val; |
772 |
|
|
|
773 |
|
|
} |
774 |
|
|
} |
775 |
|
|
lowVal = D[lowJ]; |
776 |
|
|
} |
777 |
|
|
|
778 |
|
|
S[lowI] -= lowVal; |
779 |
|
|
D[lowJ] -= lowVal; |
780 |
|
|
|
781 |
|
|
isBasic[lowI][lowJ] = 1; |
782 |
|
|
basicsEnd->val = lowVal; |
783 |
|
|
basicsEnd->i = lowI; |
784 |
|
|
basicsEnd->j = lowJ; |
785 |
|
|
|
786 |
|
|
basicsEnd->nextSnk = snkBasics[lowJ]; |
787 |
|
|
if (snkBasics[lowJ] != NULL) snkBasics[lowJ]->prevSnk = basicsEnd; |
788 |
|
|
basicsEnd->nextSrc = srcBasics[lowI]; |
789 |
|
|
if (srcBasics[lowI] != NULL) srcBasics[lowI]->prevSrc = basicsEnd; |
790 |
|
|
|
791 |
|
|
srcBasics[lowI] = basicsEnd; |
792 |
|
|
basicsEnd->prevSnk = NULL; |
793 |
|
|
snkBasics[lowJ] = basicsEnd; |
794 |
|
|
basicsEnd->prevSrc = NULL; |
795 |
|
|
|
796 |
|
|
basicsEnd++; |
797 |
|
|
|
798 |
|
|
} |
799 |
|
|
|
800 |
|
|
delete[] U; |
801 |
|
|
delete[] V; |
802 |
|
|
for (i = 0; i < n1; i++) |
803 |
|
|
delete[] Delta[i]; |
804 |
|
|
delete[] Delta; |
805 |
|
|
} |
806 |
|
|
} |
807 |
|
|
#endif |
808 |
|
|
|