34 #ifndef _CYMEALG_CYCLEMEAN_H_ 35 #define _CYMEALG_CYCLEMEAN_H_ 38 #include "chomp/system/config.h" 39 #include "chomp/system/textfile.h" 47 template <
class wType>
52 chomp::homology::multitable<int_t> compVertices, compEnds;
53 bool copyweights = !!transposed;
54 int_t countSCC =
SCC (g, compVertices, compEnds,
60 int_t maxCompSize = compEnds [0];
61 for (int_t comp = 1; comp < countSCC; ++ comp)
63 int_t compSize = compEnds [comp] - compEnds [comp - 1];
64 if (maxCompSize < compSize)
65 maxCompSize = compSize;
69 wType **F =
new wType * [maxCompSize + 1];
70 chomp::homology::BitField *finite =
71 new chomp::homology::BitField [maxCompSize + 1];
72 for (int_t i = 0; i <= maxCompSize; ++ i)
74 F [i] =
new wType [maxCompSize];
75 finite [i]. allocate (maxCompSize);
79 int_t nVertices (g. countVertices ());
80 int_t *numbers =
new int_t [nVertices];
81 int_t *components =
new int_t [nVertices];
82 for (int_t i = 0; i < nVertices; ++ i)
85 for (int_t comp = 0; comp < countSCC; ++ comp)
87 int_t maxOffset = compEnds [comp];
89 for (; pos < maxOffset; ++ pos)
91 numbers [compVertices [pos]] = pos - offset;
92 components [compVertices [pos]] = comp;
99 for (int_t comp = 0; comp < countSCC; ++ comp)
101 int_t offset = comp ? compEnds [comp - 1] :
102 static_cast<int_t
> (0);
103 int_t compSize = compEnds [comp] - offset;
104 for (int_t i = 0; i <= compSize; ++ i)
105 finite [i]. clearall (compSize);
110 int_t nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
111 short percentage = 0;
114 for (int_t len = 1; len <= compSize; ++ len)
117 for (int_t i = 0; i < compSize; ++ i)
119 if (!finite [len - 1]. test (i))
121 wType prevWeight = F [len - 1] [i];
122 int_t source = compVertices [offset + i];
125 int_t nEdges = g. countEdges (source);
126 for (
int edge = 0; edge < nEdges; ++ edge)
129 g. getEdge (source, edge);
130 if (components [target] != comp)
132 int_t j = numbers [target];
133 wType newWeight = prevWeight +
134 g. getWeight (source, edge);
135 if (!finite [len]. test (j))
137 finite [len].
set (j);
138 F [len] [j] = newWeight;
140 else if (newWeight < F [len] [j])
141 F [len] [j] = newWeight;
146 if (len >= nextProgress)
148 chomp::homology::scon << std::setw (3) <<
149 percentage <<
"%\b\b\b\b";
151 nextProgress = (percentage * compSize) / 75;
156 nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
158 wType minCompWeight (0);
159 bool firstMin =
true;
160 for (int_t vert = 0; vert < compSize; ++ vert)
162 if (!finite [compSize]. test (vert))
164 bool firstAverage =
true;
165 wType maxAverage = 0;
166 for (int_t first = 0; first < compSize; ++ first)
168 if (!finite [first]. test (vert))
170 wType average = (F [compSize] [vert] -
175 maxAverage = average;
176 firstAverage =
false;
178 else if (maxAverage < average)
179 maxAverage = average;
181 if (firstMin || (maxAverage < minCompWeight))
184 throw "Bug in Karp's algorithm";
185 minCompWeight = maxAverage;
190 if (vert >= nextProgress)
192 chomp::homology::scon << std::setw (3) <<
193 (percentage + 75) <<
"%\b\b\b\b";
195 nextProgress = (percentage * compSize) / 25;
200 if (!comp || (minCompWeight < minWeight))
201 minWeight = minCompWeight;
205 delete [] components;
207 for (int_t i = 0; i <= maxCompSize; ++ i)
230 template <
class wType,
class roundType>
235 chomp::homology::multitable<int_t> compVertices, compEnds;
236 bool copyweights = !!transposed;
237 int_t countSCC =
SCC (g, compVertices, compEnds,
243 int_t maxCompSize = compEnds [0];
244 for (int_t comp = 1; comp < countSCC; ++ comp)
246 int_t compSize = compEnds [comp] - compEnds [comp - 1];
247 if (maxCompSize < compSize)
248 maxCompSize = compSize;
252 wType **FUpper =
new wType * [maxCompSize + 1];
253 wType **FLower =
new wType * [maxCompSize + 1];
254 chomp::homology::BitField *finite =
255 new chomp::homology::BitField [maxCompSize + 1];
256 for (int_t i = 0; i <= maxCompSize; ++ i)
258 FUpper [i] =
new wType [maxCompSize];
259 FLower [i] =
new wType [maxCompSize];
260 finite [i]. allocate (maxCompSize);
264 int_t nVertices (g. countVertices ());
265 int_t *numbers =
new int_t [nVertices];
266 int_t *components =
new int_t [nVertices];
267 for (int_t i = 0; i < nVertices; ++ i)
270 for (int_t comp = 0; comp < countSCC; ++ comp)
272 int_t maxOffset = compEnds [comp];
274 for (; pos < maxOffset; ++ pos)
276 numbers [compVertices [pos]] = pos - offset;
277 components [compVertices [pos]] = comp;
284 for (int_t comp = 0; comp < countSCC; ++ comp)
286 int_t offset = comp ? compEnds [comp - 1] :
287 static_cast<int_t
> (0);
288 int_t compSize = compEnds [comp] - offset;
289 for (int_t i = 0; i <= compSize; ++ i)
290 finite [i]. clearall (compSize);
296 int_t nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
297 short percentage = 0;
300 for (int_t len = 1; len <= compSize; ++ len)
303 for (int_t i = 0; i < compSize; ++ i)
305 if (!finite [len - 1]. test (i))
307 wType prevUpper = FUpper [len - 1] [i];
308 wType prevLower = FLower [len - 1] [i];
309 int_t source = compVertices [offset + i];
312 int_t nEdges = g. countEdges (source);
313 for (
int edge = 0; edge < nEdges; ++ edge)
316 g. getEdge (source, edge);
317 if (components [target] != comp)
319 int_t j = numbers [target];
320 wType newUpper = rounding. add_up
322 g. getWeight (source, edge));
323 wType newLower = rounding. add_down
325 g. getWeight (source, edge));
326 if (!finite [len]. test (j))
328 finite [len].
set (j);
329 FUpper [len] [j] = newUpper;
330 FLower [len] [j] = newLower;
336 if (newUpper < curUpper)
340 if (newLower < curLower)
347 if (len >= nextProgress)
349 chomp::homology::scon << std::setw (3) <<
350 percentage <<
"%\b\b\b\b";
352 nextProgress = (percentage * compSize) / 75;
357 nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
359 wType minCompWeight (0);
360 bool firstMin =
true;
361 for (int_t vert = 0; vert < compSize; ++ vert)
363 if (!finite [compSize]. test (vert))
365 bool firstAverage =
true;
366 wType maxAverage = 0;
367 for (int_t first = 0; first < compSize; ++ first)
369 if (!finite [first]. test (vert))
371 const wType diff = rounding. sub_down
372 (FLower [compSize] [vert],
373 FUpper [first] [vert]);
374 wType average = rounding. div_down
375 (diff, compSize - first);
378 maxAverage = average;
379 firstAverage =
false;
381 else if (maxAverage < average)
382 maxAverage = average;
384 if (firstMin || (maxAverage < minCompWeight))
387 throw "Bug in Karp's algorithm";
388 minCompWeight = maxAverage;
393 if (vert >= nextProgress)
395 chomp::homology::scon << std::setw (3) <<
396 (percentage + 75) <<
"%\b\b\b\b";
398 nextProgress = (percentage * compSize) / 25;
403 if (!comp || (minCompWeight < minWeight))
404 minWeight = minCompWeight;
408 delete [] components;
410 for (int_t i = 0; i <= maxCompSize; ++ i)
413 delete [] (FUpper [i]);
414 delete [] (FLower [i]);
432 template <
class wType,
class roundType>
439 chomp::homology::multitable<int_t> compVertices, compEnds;
440 bool copyweights = !!transposed;
441 int_t countSCC =
SCC (g, compVertices, compEnds,
447 int_t maxCompSize = compEnds [0];
448 for (int_t comp = 1; comp < countSCC; ++ comp)
450 int_t compSize = compEnds [comp] - compEnds [comp - 1];
451 if (maxCompSize < compSize)
452 maxCompSize = compSize;
458 const int nTables (5);
459 wType **FUpper =
new wType * [nTables];
460 wType **FLower =
new wType * [nTables];
461 #ifdef CHECK_ROUNDING_WIDTH 462 wType **FUpperE =
new wType * [nTables];
463 wType **FLowerE =
new wType * [nTables];
465 chomp::homology::BitField *finite =
466 new chomp::homology::BitField [nTables];
467 for (int_t i = 0; i < nTables; ++ i)
469 FUpper [i] =
new wType [maxCompSize];
470 FLower [i] =
new wType [maxCompSize];
471 #ifdef CHECK_ROUNDING_WIDTH 472 FUpperE [i] =
new wType [maxCompSize];
473 FLowerE [i] =
new wType [maxCompSize];
475 finite [i]. allocate (maxCompSize);
479 int_t nVertices (g. countVertices ());
480 int_t *numbers =
new int_t [nVertices];
481 int_t *components =
new int_t [nVertices];
482 for (int_t i = 0; i < nVertices; ++ i)
485 for (int_t comp = 0; comp < countSCC; ++ comp)
487 int_t maxOffset = compEnds [comp];
489 for (; pos < maxOffset; ++ pos)
491 numbers [compVertices [pos]] = pos - offset;
492 components [compVertices [pos]] = comp;
499 #ifdef CHECK_ROUNDING_WIDTH 500 wType minWeightE (0);
502 for (int_t comp = 0; comp < countSCC; ++ comp)
504 int_t offset = comp ? compEnds [comp - 1] :
505 static_cast<int_t
> (0);
506 int_t compSize = compEnds [comp] - offset;
507 finite [0]. clearall (compSize);
508 finite [4]. clearall (compSize);
511 #ifdef CHECK_ROUNDING_WIDTH 518 int_t nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
519 short percentage = 0;
522 const int_t maxLength1 (compSize + 1);
523 for (int_t len = 1; len < maxLength1; ++ len)
526 int_t prevIndex ((len == 1) ?
527 0 : (((len - 1) & 1) + 1));
528 int_t curIndex ((len == compSize) ?
529 3 : ((len & 1) + 1));
530 finite [curIndex]. clearall (compSize);
534 for (int_t i = finite [prevIndex].
536 (i >= 0) && (i < compSize);
537 i = finite [prevIndex]. find
543 int_t source = compVertices [offset + i];
546 int_t nEdges = g. countEdges (source);
547 for (
int edge = 0; edge < nEdges; ++ edge)
550 g. getEdge (source, edge);
551 if (components [target] != comp)
553 int_t j = numbers [target];
554 wType newLower = rounding. add_down
555 (FLower [prevIndex] [i],
556 g. getWeight (source, edge));
557 #ifdef CHECK_ROUNDING_WIDTH 558 wType newLowerE = rounding. add_up
559 (FLowerE [prevIndex] [i],
560 g. getWeight (source, edge));
566 if (!finite [curIndex]. test (j))
568 finite [curIndex].
set (j);
569 FLower [curIndex] [j] =
571 #ifdef CHECK_ROUNDING_WIDTH 572 FLowerE [curIndex] [j] =
579 FLower [curIndex] [j])
581 FLower [curIndex] [j] =
583 #ifdef CHECK_ROUNDING_WIDTH 584 FLowerE [curIndex] [j] =
592 if (len >= nextProgress)
594 chomp::homology::scon << std::setw (3) <<
595 percentage <<
"%\b\b\b\b";
597 nextProgress = (percentage * compSize) / 50;
603 FLower [4] [0] = rounding. div_down (FLower [3] [0],
605 #ifdef CHECK_ROUNDING_WIDTH 606 FLowerE [4] [0] = rounding. div_up (FLowerE [3] [0],
612 nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
614 const int_t maxLength2 (compSize);
615 for (int_t len = 1; len < maxLength2; ++ len)
618 int_t prevIndex ((len == 1) ?
619 0 : (((len - 1) & 1) + 1));
620 int_t curIndex ((len & 1) + 1);
621 finite [curIndex]. clearall (compSize);
625 for (int_t i = finite [prevIndex].
627 (i >= 0) && (i < compSize);
628 i = finite [prevIndex]. find
632 if (!finite [prevIndex]. test (i))
633 throw "Programming bug: infinite.";
634 int_t source = compVertices [offset + i];
637 int_t nEdges = g. countEdges (source);
638 for (
int edge = 0; edge < nEdges; ++ edge)
641 g. getEdge (source, edge);
642 if (components [target] != comp)
644 int_t j = numbers [target];
645 wType newUpper = rounding. add_up
646 (FUpper [prevIndex] [i],
647 g. getWeight (source, edge));
648 #ifdef CHECK_ROUNDING_WIDTH 649 wType newUpperE = rounding. add_down
650 (FUpperE [prevIndex] [i],
651 g. getWeight (source, edge));
653 bool checkAverage =
false;
658 if (!finite [curIndex]. test (j))
660 finite [curIndex].
set (j);
661 FUpper [curIndex] [j] =
663 #ifdef CHECK_ROUNDING_WIDTH 664 FUpperE [curIndex] [j] =
672 FUpper [curIndex] [j])
674 FUpper [curIndex] [j] =
676 #ifdef CHECK_ROUNDING_WIDTH 677 FUpperE [curIndex] [j] =
686 const wType diff = rounding. sub_down
688 FUpper [curIndex] [j]);
689 #ifdef CHECK_ROUNDING_WIDTH 690 const wType diffE = rounding. sub_up
692 FUpperE [curIndex] [j]);
694 const wType average = rounding.
697 #ifdef CHECK_ROUNDING_WIDTH 698 const wType averageE = rounding.
702 if (!finite [4]. test (j))
705 FLower [4] [j] = average;
706 #ifdef CHECK_ROUNDING_WIDTH 707 FLowerE [4] [j] = averageE;
710 else if (FLower [4] [j] < average)
712 FLower [4] [j] = average;
713 #ifdef CHECK_ROUNDING_WIDTH 714 FLowerE [4] [j] = averageE;
721 if (len >= nextProgress)
723 chomp::homology::scon << std::setw (3) <<
724 (percentage + 50) <<
"%\b\b\b\b";
726 nextProgress = (percentage * compSize) / 50;
731 wType minCompWeight (0);
732 #ifdef CHECK_ROUNDING_WIDTH 733 wType minCompWeightE (0);
735 for (int_t vert = 0; vert < compSize; ++ vert)
739 if (!vert || (FLower [4] [vert] < minCompWeight))
741 minCompWeight = FLower [4] [vert];
742 #ifdef CHECK_ROUNDING_WIDTH 743 minCompWeightE = FLowerE [4] [vert];
749 if (!comp || (minCompWeight < minWeight))
751 minWeight = minCompWeight;
752 #ifdef CHECK_ROUNDING_WIDTH 753 minWeightE = minCompWeightE;
759 delete [] components;
761 for (int_t i = 0; i < nTables; ++ i)
764 delete [] (FUpper [i]);
765 delete [] (FLower [i]);
766 #ifdef CHECK_ROUNDING_WIDTH 767 delete [] (FUpperE [i]);
768 delete [] (FLowerE [i]);
774 #ifdef CHECK_ROUNDING_WIDTH 779 #ifdef CHECK_ROUNDING_WIDTH 781 chomp::homology::sbug <<
"\nLoss of precision info: [" <<
782 minWeight <<
", " << minWeightE <<
783 "], width: " << (minWeightE - minWeight) <<
".\n";
793 #endif // _CYMEALG_CYCLEMEAN_H_ wType minMeanCycleWeightMem(const diGraph< wType > &g, const roundType &rounding, diGraph< wType > *transposed)
A rigorous numerics version of Karp algorithm modified in such a way that the memory usage is minimiz...
This class defines a weighted directed graph with very limited number of operations.
int_t SCC(const diGraph< wType > &g, Table1 &compVertices, Table2 &compEnds, diGraph< wType > *scc=0, diGraph< wType > *transposed=0, bool copyweights=false)
Computes strongly connected path components of the graph 'g'.
wType minMeanCycleWeight(const diGraph< wType > &g, diGraph< wType > *transposed=0)
Runs the Karp algorithm for each strongly connected component of the graph and returns the minimum me...