The CyMeAlg Software (Release 0.01)
cyclemean.h
Go to the documentation of this file.
1 /// @addtogroup cymealg
2 /// @{
3 
4 /////////////////////////////////////////////////////////////////////////////
5 ///
6 /// @file
7 ///
8 /// This header file contains an implementation of a few variants
9 /// of algorithms for the computation of the minimum cycle mean of a digraph.
10 ///
11 /// @author Pawel Pilarczyk
12 ///
13 /////////////////////////////////////////////////////////////////////////////
14 
15 // Copyright (C) 1997-2020 by Pawel Pilarczyk.
16 //
17 // This file is part of my research software. This is free software:
18 // you can redistribute it and/or modify it under the terms of the GNU
19 // General Public License as published by the Free Software Foundation,
20 // either version 3 of the License, or (at your option) any later version.
21 //
22 // This software is distributed in the hope that it will be useful,
23 // but WITHOUT ANY WARRANTY; without even the implied warranty of
24 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 // GNU General Public License for more details.
26 //
27 // You should have received a copy of the GNU General Public License
28 // along with this software; see the file "license.txt". If not,
29 // please, see <http://www.gnu.org/licenses/>.
30 
31 // Started in January 2006. Last revision: August 23, 2019.
32 
33 
34 #ifndef _CYMEALG_CYCLEMEAN_H_
35 #define _CYMEALG_CYCLEMEAN_H_
36 
37 // include selected header files from the CHomP library
38 #include "chomp/system/config.h"
39 #include "chomp/system/textfile.h"
40 
41 namespace cymealg {
42 
43 /// Runs the Karp algorithm for each strongly connected component
44 /// of the graph and returns the minimum mean cycle weight,
45 /// which can be negative.
46 /// As a byproduct, saves the transposed graph, if requested to.
47 template <class wType>
49  diGraph<wType> *transposed = 0)
50 {
51  // find the strongly connected components of the graph
52  chomp::homology::multitable<int_t> compVertices, compEnds;
53  bool copyweights = !!transposed;
54  int_t countSCC = SCC (g, compVertices, compEnds,
55  static_cast<diGraph<wType> *> (0), transposed, copyweights);
56  if (countSCC <= 0)
57  return 0;
58 
59  // compute the maximum size of each strongly connected component
60  int_t maxCompSize = compEnds [0];
61  for (int_t comp = 1; comp < countSCC; ++ comp)
62  {
63  int_t compSize = compEnds [comp] - compEnds [comp - 1];
64  if (maxCompSize < compSize)
65  maxCompSize = compSize;
66  }
67 
68  // allocate arrays for weights and bit fields
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)
73  {
74  F [i] = new wType [maxCompSize];
75  finite [i]. allocate (maxCompSize);
76  }
77 
78  // compute the numbers of vertices in each component
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)
83  components [i] = -1;
84  int_t offset = 0;
85  for (int_t comp = 0; comp < countSCC; ++ comp)
86  {
87  int_t maxOffset = compEnds [comp];
88  int_t pos = offset;
89  for (; pos < maxOffset; ++ pos)
90  {
91  numbers [compVertices [pos]] = pos - offset;
92  components [compVertices [pos]] = comp;
93  }
94  offset = pos;
95  }
96 
97  // compute the minimum mean cycle weight for each component
98  wType minWeight (0);
99  for (int_t comp = 0; comp < countSCC; ++ comp)
100  {
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);
106  F [0] [0] = 0;
107  finite [0]. set (0);
108 
109  // when to show the next progress indicator?
110  int_t nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
111  short percentage = 0;
112 
113  // compute path progressions of given length
114  for (int_t len = 1; len <= compSize; ++ len)
115  {
116  // process source vertices
117  for (int_t i = 0; i < compSize; ++ i)
118  {
119  if (!finite [len - 1]. test (i))
120  continue;
121  wType prevWeight = F [len - 1] [i];
122  int_t source = compVertices [offset + i];
123 
124  // process target vertices (and edges)
125  int_t nEdges = g. countEdges (source);
126  for (int edge = 0; edge < nEdges; ++ edge)
127  {
128  int_t target =
129  g. getEdge (source, edge);
130  if (components [target] != comp)
131  continue;
132  int_t j = numbers [target];
133  wType newWeight = prevWeight +
134  g. getWeight (source, edge);
135  if (!finite [len]. test (j))
136  {
137  finite [len]. set (j);
138  F [len] [j] = newWeight;
139  }
140  else if (newWeight < F [len] [j])
141  F [len] [j] = newWeight;
142  }
143  }
144 
145  // show progress indicator if necessary
146  if (len >= nextProgress)
147  {
148  chomp::homology::scon << std::setw (3) <<
149  percentage << "%\b\b\b\b";
150  ++ percentage;
151  nextProgress = (percentage * compSize) / 75;
152  }
153  }
154 
155  // compute the minimum mean cycle weight for this component
156  nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
157  percentage = 0;
158  wType minCompWeight (0);
159  bool firstMin = true;
160  for (int_t vert = 0; vert < compSize; ++ vert)
161  {
162  if (!finite [compSize]. test (vert))
163  continue;
164  bool firstAverage = true;
165  wType maxAverage = 0;
166  for (int_t first = 0; first < compSize; ++ first)
167  {
168  if (!finite [first]. test (vert))
169  continue;
170  wType average = (F [compSize] [vert] -
171  F [first] [vert]) /
172  (compSize - first);
173  if (firstAverage)
174  {
175  maxAverage = average;
176  firstAverage = false;
177  }
178  else if (maxAverage < average)
179  maxAverage = average;
180  }
181  if (firstMin || (maxAverage < minCompWeight))
182  {
183  if (firstAverage)
184  throw "Bug in Karp's algorithm";
185  minCompWeight = maxAverage;
186  firstMin = false;
187  }
188 
189  // show progress indicator if necessary
190  if (vert >= nextProgress)
191  {
192  chomp::homology::scon << std::setw (3) <<
193  (percentage + 75) << "%\b\b\b\b";
194  ++ percentage;
195  nextProgress = (percentage * compSize) / 25;
196  }
197  }
198 
199  // make a correction to the total minimum if necessary
200  if (!comp || (minCompWeight < minWeight))
201  minWeight = minCompWeight;
202  }
203 
204  // release allocated memory
205  delete [] components;
206  delete [] numbers;
207  for (int_t i = 0; i <= maxCompSize; ++ i)
208  {
209  finite [i]. free ();
210  delete [] (F [i]);
211  }
212  delete [] finite;
213  delete [] F;
214 
215  // return the computed minimum
216  return minWeight;
217 } /* minMeanCycleWeight */
218 
219 
220 // --------------------------------------------------
221 
222 /// A version of Karp algorithm modified for the purpose
223 /// of interval arithmetic to provide the correct lower bound
224 /// for the minimum mean cycle weight in a graph.
225 /// This specialization is necessary, because of the subtraction
226 /// operation used in Karp's algorithm. Therefore, upper and lower
227 /// bounds for the minimum path progression weights must be computed
228 /// independently.
229 /// The intervals are compared by comparing their lower bounds only.
230 template <class wType, class roundType>
232  const roundType &rounding, diGraph<wType> *transposed)
233 {
234  // find the strongly connected components of the graph
235  chomp::homology::multitable<int_t> compVertices, compEnds;
236  bool copyweights = !!transposed;
237  int_t countSCC = SCC (g, compVertices, compEnds,
238  static_cast<diGraph<wType> *> (0), transposed, copyweights);
239  if (countSCC <= 0)
240  return 0;
241 
242  // compute the maximum size of each strongly connected component
243  int_t maxCompSize = compEnds [0];
244  for (int_t comp = 1; comp < countSCC; ++ comp)
245  {
246  int_t compSize = compEnds [comp] - compEnds [comp - 1];
247  if (maxCompSize < compSize)
248  maxCompSize = compSize;
249  }
250 
251  // allocate arrays for weights and bit fields
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)
257  {
258  FUpper [i] = new wType [maxCompSize];
259  FLower [i] = new wType [maxCompSize];
260  finite [i]. allocate (maxCompSize);
261  }
262 
263  // compute the numbers of vertices in each component
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)
268  components [i] = -1;
269  int_t offset = 0;
270  for (int_t comp = 0; comp < countSCC; ++ comp)
271  {
272  int_t maxOffset = compEnds [comp];
273  int_t pos = offset;
274  for (; pos < maxOffset; ++ pos)
275  {
276  numbers [compVertices [pos]] = pos - offset;
277  components [compVertices [pos]] = comp;
278  }
279  offset = pos;
280  }
281 
282  // compute the minimum mean cycle weight for each component
283  wType minWeight (0);
284  for (int_t comp = 0; comp < countSCC; ++ comp)
285  {
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);
291  FUpper [0] [0] = 0;
292  FLower [0] [0] = 0;
293  finite [0]. set (0);
294 
295  // when to show the next progress indicator?
296  int_t nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
297  short percentage = 0;
298 
299  // compute path progressions of given length
300  for (int_t len = 1; len <= compSize; ++ len)
301  {
302  // process source vertices
303  for (int_t i = 0; i < compSize; ++ i)
304  {
305  if (!finite [len - 1]. test (i))
306  continue;
307  wType prevUpper = FUpper [len - 1] [i];
308  wType prevLower = FLower [len - 1] [i];
309  int_t source = compVertices [offset + i];
310 
311  // process target vertices (and edges)
312  int_t nEdges = g. countEdges (source);
313  for (int edge = 0; edge < nEdges; ++ edge)
314  {
315  int_t target =
316  g. getEdge (source, edge);
317  if (components [target] != comp)
318  continue;
319  int_t j = numbers [target];
320  wType newUpper = rounding. add_up
321  (prevUpper,
322  g. getWeight (source, edge));
323  wType newLower = rounding. add_down
324  (prevLower,
325  g. getWeight (source, edge));
326  if (!finite [len]. test (j))
327  {
328  finite [len]. set (j);
329  FUpper [len] [j] = newUpper;
330  FLower [len] [j] = newLower;
331  }
332  else
333  {
334  wType &curUpper =
335  FUpper [len] [j];
336  if (newUpper < curUpper)
337  curUpper = newUpper;
338  wType &curLower =
339  FLower [len] [j];
340  if (newLower < curLower)
341  curLower = newLower;
342  }
343  }
344  }
345 
346  // show progress indicator if necessary
347  if (len >= nextProgress)
348  {
349  chomp::homology::scon << std::setw (3) <<
350  percentage << "%\b\b\b\b";
351  ++ percentage;
352  nextProgress = (percentage * compSize) / 75;
353  }
354  }
355 
356  // compute the minimum mean cycle weight for this component
357  nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
358  percentage = 0;
359  wType minCompWeight (0);
360  bool firstMin = true;
361  for (int_t vert = 0; vert < compSize; ++ vert)
362  {
363  if (!finite [compSize]. test (vert))
364  continue;
365  bool firstAverage = true;
366  wType maxAverage = 0;
367  for (int_t first = 0; first < compSize; ++ first)
368  {
369  if (!finite [first]. test (vert))
370  continue;
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);
376  if (firstAverage)
377  {
378  maxAverage = average;
379  firstAverage = false;
380  }
381  else if (maxAverage < average)
382  maxAverage = average;
383  }
384  if (firstMin || (maxAverage < minCompWeight))
385  {
386  if (firstAverage)
387  throw "Bug in Karp's algorithm";
388  minCompWeight = maxAverage;
389  firstMin = false;
390  }
391 
392  // show progress indicator if necessary
393  if (vert >= nextProgress)
394  {
395  chomp::homology::scon << std::setw (3) <<
396  (percentage + 75) << "%\b\b\b\b";
397  ++ percentage;
398  nextProgress = (percentage * compSize) / 25;
399  }
400  }
401 
402  // make a correction to the total minimum if necessary
403  if (!comp || (minCompWeight < minWeight))
404  minWeight = minCompWeight;
405  }
406 
407  // release allocated memory
408  delete [] components;
409  delete [] numbers;
410  for (int_t i = 0; i <= maxCompSize; ++ i)
411  {
412  finite [i]. free ();
413  delete [] (FUpper [i]);
414  delete [] (FLower [i]);
415  }
416  delete [] finite;
417  delete [] FUpper;
418  delete [] FLower;
419 
420  // return the computed minimum
421  return minWeight;
422 } /* minMeanCycleWeight */
423 
424 
425 // --------------------------------------------------
426 
427 /// A rigorous numerics version of Karp algorithm
428 /// modified in such a way that the memory usage is minimized,
429 /// at the cost of slower execution (up to twice slower).
430 /// Note that a memory conservative version of the algorithm
431 /// that uses non-rigorous numerics is not programmed yet.
432 template <class wType, class roundType>
434  const roundType &rounding, diGraph<wType> *transposed)
435 {
436 // NOTE: Support for CHECK_ROUNDING_WIDTH is currently in preparation.
437 
438  // find the strongly connected components of the graph
439  chomp::homology::multitable<int_t> compVertices, compEnds;
440  bool copyweights = !!transposed;
441  int_t countSCC = SCC (g, compVertices, compEnds,
442  static_cast<diGraph<wType> *> (0), transposed, copyweights);
443  if (countSCC <= 0)
444  return 0;
445 
446  // compute the maximum size of each strongly connected component
447  int_t maxCompSize = compEnds [0];
448  for (int_t comp = 1; comp < countSCC; ++ comp)
449  {
450  int_t compSize = compEnds [comp] - compEnds [comp - 1];
451  if (maxCompSize < compSize)
452  maxCompSize = compSize;
453  }
454 
455  // allocate arrays for weights and bit fields;
456  // table 0 for F [0], tables 1 and 2 for F [k - 1] and F [k],
457  // table 3 for F [n], table 4 for max of the quotients
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];
464 #endif
465  chomp::homology::BitField *finite =
466  new chomp::homology::BitField [nTables];
467  for (int_t i = 0; i < nTables; ++ i)
468  {
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];
474 #endif
475  finite [i]. allocate (maxCompSize);
476  }
477 
478  // compute the numbers of vertices in each component
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)
483  components [i] = -1;
484  int_t offset = 0;
485  for (int_t comp = 0; comp < countSCC; ++ comp)
486  {
487  int_t maxOffset = compEnds [comp];
488  int_t pos = offset;
489  for (; pos < maxOffset; ++ pos)
490  {
491  numbers [compVertices [pos]] = pos - offset;
492  components [compVertices [pos]] = comp;
493  }
494  offset = pos;
495  }
496 
497  // compute the minimum mean cycle weight for each component
498  wType minWeight (0);
499 #ifdef CHECK_ROUNDING_WIDTH
500  wType minWeightE (0);
501 #endif
502  for (int_t comp = 0; comp < countSCC; ++ comp)
503  {
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);
509  FUpper [0] [0] = 0;
510  FLower [0] [0] = 0;
511 #ifdef CHECK_ROUNDING_WIDTH
512  FUpperE [0] [0] = 0;
513  FLowerE [0] [0] = 0;
514 #endif
515  finite [0]. set (0);
516 
517  // when to show the next progress indicator?
518  int_t nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
519  short percentage = 0;
520 
521  // Pass 1: compute path progressions of maximal length
522  const int_t maxLength1 (compSize + 1);
523  for (int_t len = 1; len < maxLength1; ++ len)
524  {
525  // determine tables to use and clear the current one
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);
531 
532  // process all the edges starting at vertices
533  // to which finite paths are already known
534  for (int_t i = finite [prevIndex].
535  find (0, compSize);
536  (i >= 0) && (i < compSize);
537  i = finite [prevIndex]. find
538  (i + 1, compSize))
539  // for (int_t i = 0; i < compSize; ++ i)
540  {
541  // if (!finite [prevIndex]. test (i))
542  // throw "Programming bug: infinite.";
543  int_t source = compVertices [offset + i];
544 
545  // process target vertices (and edges)
546  int_t nEdges = g. countEdges (source);
547  for (int edge = 0; edge < nEdges; ++ edge)
548  {
549  int_t target =
550  g. getEdge (source, edge);
551  if (components [target] != comp)
552  continue;
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));
561 #endif
562 
563  // if this is the first time that
564  // a path of this length
565  // to the vertex was encountered
566  if (!finite [curIndex]. test (j))
567  {
568  finite [curIndex]. set (j);
569  FLower [curIndex] [j] =
570  newLower;
571 #ifdef CHECK_ROUNDING_WIDTH
572  FLowerE [curIndex] [j] =
573  newLowerE;
574 #endif
575  }
576  // if a better path of this length
577  // to the vertex was found
578  else if (newLower <
579  FLower [curIndex] [j])
580  {
581  FLower [curIndex] [j] =
582  newLower;
583 #ifdef CHECK_ROUNDING_WIDTH
584  FLowerE [curIndex] [j] =
585  newLowerE;
586 #endif
587  }
588  }
589  }
590 
591  // show progress indicator if necessary
592  if (len >= nextProgress)
593  {
594  chomp::homology::scon << std::setw (3) <<
595  percentage << "%\b\b\b\b";
596  ++ percentage;
597  nextProgress = (percentage * compSize) / 50;
598  }
599  }
600 
601  // set up the mean average at level 0
602  finite [4]. set (0);
603  FLower [4] [0] = rounding. div_down (FLower [3] [0],
604  compSize);
605 #ifdef CHECK_ROUNDING_WIDTH
606  FLowerE [4] [0] = rounding. div_up (FLowerE [3] [0],
607  compSize);
608 #endif
609 
610  // Pass 2: compute path progressions of all the lengths and
611  // update the max values of the quotients for each vertex
612  nextProgress = (compSize < 1000) ? (compSize + 1) : 0;
613  percentage = 0;
614  const int_t maxLength2 (compSize);
615  for (int_t len = 1; len < maxLength2; ++ len)
616  {
617  // determine tables to use and clear the current one
618  int_t prevIndex ((len == 1) ?
619  0 : (((len - 1) & 1) + 1));
620  int_t curIndex ((len & 1) + 1);
621  finite [curIndex]. clearall (compSize);
622 
623  // process all the edges starting at vertices
624  // to which finite paths are already known
625  for (int_t i = finite [prevIndex].
626  find (0, compSize);
627  (i >= 0) && (i < compSize);
628  i = finite [prevIndex]. find
629  (i + 1, compSize))
630  // for (int_t i = 0; i < compSize; ++ i)
631  {
632  if (!finite [prevIndex]. test (i))
633  throw "Programming bug: infinite.";
634  int_t source = compVertices [offset + i];
635 
636  // process target vertices (and edges)
637  int_t nEdges = g. countEdges (source);
638  for (int edge = 0; edge < nEdges; ++ edge)
639  {
640  int_t target =
641  g. getEdge (source, edge);
642  if (components [target] != comp)
643  continue;
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));
652 #endif
653  bool checkAverage = false;
654 
655  // if this is the first time that
656  // a path of this length
657  // to the vertex was encountered
658  if (!finite [curIndex]. test (j))
659  {
660  finite [curIndex]. set (j);
661  FUpper [curIndex] [j] =
662  newUpper;
663 #ifdef CHECK_ROUNDING_WIDTH
664  FUpperE [curIndex] [j] =
665  newUpperE;
666 #endif
667  checkAverage = true;
668  }
669  // if a better path of this length
670  // to the vertex was found
671  else if (newUpper <
672  FUpper [curIndex] [j])
673  {
674  FUpper [curIndex] [j] =
675  newUpper;
676 #ifdef CHECK_ROUNDING_WIDTH
677  FUpperE [curIndex] [j] =
678  newUpperE;
679 #endif
680  checkAverage = true;
681  }
682 
683  // update the max average weight
684  if (!checkAverage)
685  continue;
686  const wType diff = rounding. sub_down
687  (FLower [3] [j],
688  FUpper [curIndex] [j]);
689 #ifdef CHECK_ROUNDING_WIDTH
690  const wType diffE = rounding. sub_up
691  (FLowerE [3] [j],
692  FUpperE [curIndex] [j]);
693 #endif
694  const wType average = rounding.
695  div_down (diff,
696  compSize - len);
697 #ifdef CHECK_ROUNDING_WIDTH
698  const wType averageE = rounding.
699  div_up (diffE,
700  compSize - len);
701 #endif
702  if (!finite [4]. test (j))
703  {
704  finite [4]. set (j);
705  FLower [4] [j] = average;
706 #ifdef CHECK_ROUNDING_WIDTH
707  FLowerE [4] [j] = averageE;
708 #endif
709  }
710  else if (FLower [4] [j] < average)
711  {
712  FLower [4] [j] = average;
713 #ifdef CHECK_ROUNDING_WIDTH
714  FLowerE [4] [j] = averageE;
715 #endif
716  }
717  }
718  }
719 
720  // show progress indicator if necessary
721  if (len >= nextProgress)
722  {
723  chomp::homology::scon << std::setw (3) <<
724  (percentage + 50) << "%\b\b\b\b";
725  ++ percentage;
726  nextProgress = (percentage * compSize) / 50;
727  }
728  }
729 
730  // compute the minimum mean cycle weight for this component
731  wType minCompWeight (0);
732 #ifdef CHECK_ROUNDING_WIDTH
733  wType minCompWeightE (0);
734 #endif
735  for (int_t vert = 0; vert < compSize; ++ vert)
736  {
737  // if (!finite [4]. test (vert))
738  // throw "Programming bug: Undefined average.";
739  if (!vert || (FLower [4] [vert] < minCompWeight))
740  {
741  minCompWeight = FLower [4] [vert];
742 #ifdef CHECK_ROUNDING_WIDTH
743  minCompWeightE = FLowerE [4] [vert];
744 #endif
745  }
746  }
747 
748  // make a correction to the total minimum if necessary
749  if (!comp || (minCompWeight < minWeight))
750  {
751  minWeight = minCompWeight;
752 #ifdef CHECK_ROUNDING_WIDTH
753  minWeightE = minCompWeightE;
754 #endif
755  }
756  }
757 
758  // release allocated memory
759  delete [] components;
760  delete [] numbers;
761  for (int_t i = 0; i < nTables; ++ i)
762  {
763  finite [i]. free ();
764  delete [] (FUpper [i]);
765  delete [] (FLower [i]);
766 #ifdef CHECK_ROUNDING_WIDTH
767  delete [] (FUpperE [i]);
768  delete [] (FLowerE [i]);
769 #endif
770  }
771  delete [] finite;
772  delete [] FUpper;
773  delete [] FLower;
774 #ifdef CHECK_ROUNDING_WIDTH
775  delete [] FUpperE;
776  delete [] FLowerE;
777 #endif
778 
779 #ifdef CHECK_ROUNDING_WIDTH
780  // show the loss of precision caused by rounding
781  chomp::homology::sbug << "\nLoss of precision info: [" <<
782  minWeight << ", " << minWeightE <<
783  "], width: " << (minWeightE - minWeight) << ".\n";
784 #endif
785 
786  // return the computed minimum
787  return minWeight;
788 } /* minMeanCycleWeightMem */
789 
790 
791 } // namespace cymealg
792 
793 #endif // _CYMEALG_CYCLEMEAN_H_
794 
795 /// @}
796 
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...
Definition: cyclemean.h:433
This class defines a weighted directed graph with very limited number of operations.
Definition: digraph.h:75
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 &#39;g&#39;.
Definition: scc.h:59
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...
Definition: cyclemean.h:48