Changes in / [1196:959d78f3fe0e:1195:73e29215aaa4] in lemon-main
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
CMakeLists.txt
r1185 r1138 87 87 IF(LEMON_ENABLE_GLPK) 88 88 FIND_PACKAGE(GLPK 4.33) 89 IF(GLPK_FOUND)90 SET(LEMON_HAVE_LP TRUE)91 SET(LEMON_HAVE_MIP TRUE)92 SET(LEMON_HAVE_GLPK TRUE)93 ENDIF(GLPK_FOUND)94 89 ENDIF(LEMON_ENABLE_GLPK) 95 90 IF(LEMON_ENABLE_ILOG) 96 91 FIND_PACKAGE(ILOG) 97 IF(ILOG_FOUND)98 SET(LEMON_HAVE_LP TRUE)99 SET(LEMON_HAVE_MIP TRUE)100 SET(LEMON_HAVE_CPLEX TRUE)101 ENDIF(ILOG_FOUND)102 92 ENDIF(LEMON_ENABLE_ILOG) 103 93 IF(LEMON_ENABLE_COIN) 104 94 FIND_PACKAGE(COIN) 105 IF(COIN_FOUND)106 SET(LEMON_HAVE_LP TRUE)107 SET(LEMON_HAVE_MIP TRUE)108 SET(LEMON_HAVE_CLP TRUE)109 SET(LEMON_HAVE_CBC TRUE)110 ENDIF(COIN_FOUND)111 95 ENDIF(LEMON_ENABLE_COIN) 112 96 IF(LEMON_ENABLE_SOPLEX) 113 97 FIND_PACKAGE(SOPLEX) 114 IF(SOPLEX_FOUND)115 SET(LEMON_HAVE_LP TRUE)116 SET(LEMON_HAVE_SOPLEX TRUE)117 ENDIF(SOPLEX_FOUND)118 98 ENDIF(LEMON_ENABLE_SOPLEX) 99 100 IF(GLPK_FOUND) 101 SET(LEMON_HAVE_LP TRUE) 102 SET(LEMON_HAVE_MIP TRUE) 103 SET(LEMON_HAVE_GLPK TRUE) 104 ENDIF(GLPK_FOUND) 105 IF(ILOG_FOUND) 106 SET(LEMON_HAVE_LP TRUE) 107 SET(LEMON_HAVE_MIP TRUE) 108 SET(LEMON_HAVE_CPLEX TRUE) 109 ENDIF(ILOG_FOUND) 110 IF(COIN_FOUND) 111 SET(LEMON_HAVE_LP TRUE) 112 SET(LEMON_HAVE_MIP TRUE) 113 SET(LEMON_HAVE_CLP TRUE) 114 SET(LEMON_HAVE_CBC TRUE) 115 ENDIF(COIN_FOUND) 116 IF(SOPLEX_FOUND) 117 SET(LEMON_HAVE_LP TRUE) 118 SET(LEMON_HAVE_SOPLEX TRUE) 119 ENDIF(SOPLEX_FOUND) 119 120 120 121 IF(ILOG_FOUND) -
cmake/FindCOIN.cmake
r1180 r1064 66 66 ) 67 67 68 FIND_LIBRARY(COIN_PTHREADS_LIBRARY69 NAMES pthreads libpthreads70 HINTS ${COIN_ROOT_DIR}/lib/coin71 HINTS ${COIN_ROOT_DIR}/lib72 )73 74 68 INCLUDE(FindPackageHandleStandardArgs) 75 69 FIND_PACKAGE_HANDLE_STANDARD_ARGS(COIN DEFAULT_MSG … … 89 83 IF(COIN_FOUND) 90 84 SET(COIN_INCLUDE_DIRS ${COIN_INCLUDE_DIR}) 91 SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARY};${COIN_COIN_UTILS_LIBRARY} ")85 SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARY};${COIN_COIN_UTILS_LIBRARY};${COIN_ZLIB_LIBRARY};${COIN_BZ2_LIBRARY}") 92 86 IF(COIN_ZLIB_LIBRARY) 93 87 SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARIES};${COIN_ZLIB_LIBRARY}") … … 96 90 SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARIES};${COIN_BZ2_LIBRARY}") 97 91 ENDIF(COIN_BZ2_LIBRARY) 98 IF(COIN_PTHREADS_LIBRARY) 99 SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARIES};${COIN_PTHREADS_LIBRARY}") 100 ENDIF(COIN_PTHREADS_LIBRARY) 101 SET(COIN_CBC_LIBRARIES "${COIN_CBC_LIBRARY};${COIN_CBC_SOLVER_LIBRARY};${COIN_CGL_LIBRARY};${COIN_OSI_LIBRARY};${COIN_OSI_CBC_LIBRARY};${COIN_OSI_CLP_LIBRARY};${COIN_CLP_LIBRARIES}") 92 SET(COIN_CBC_LIBRARIES "${COIN_CBC_LIBRARY};${COIN_CBC_SOLVER_LIBRARY};${COIN_CGL_LIBRARY};${COIN_OSI_LIBRARY};${COIN_OSI_CBC_LIBRARY};${COIN_OSI_CLP_LIBRARY};${COIN_ZLIB_LIBRARY};${COIN_BZ2_LIBRARY};${COIN_CLP_LIBRARIES}") 102 93 SET(COIN_LIBRARIES ${COIN_CBC_LIBRARIES}) 103 94 ENDIF(COIN_FOUND) -
lemon/adaptors.h
r1184 r1092 3447 3447 /// This map adaptor class adapts two node maps of the original digraph 3448 3448 /// to get a node map of the split digraph. 3449 /// Its value type is inherited from the first node map type (\c I n).3450 /// \tparam I nThe type of the node map for the in-nodes.3451 /// \tparam O utThe type of the node map for the out-nodes.3452 template <typename I n, typename Out>3449 /// Its value type is inherited from the first node map type (\c IN). 3450 /// \tparam IN The type of the node map for the in-nodes. 3451 /// \tparam OUT The type of the node map for the out-nodes. 3452 template <typename IN, typename OUT> 3453 3453 class CombinedNodeMap { 3454 3454 public: … … 3457 3457 typedef Node Key; 3458 3458 /// The value type of the map 3459 typedef typename I n::Value Value;3460 3461 typedef typename MapTraits<I n>::ReferenceMapTag ReferenceMapTag;3462 typedef typename MapTraits<I n>::ReturnValue ReturnValue;3463 typedef typename MapTraits<I n>::ConstReturnValue ConstReturnValue;3464 typedef typename MapTraits<I n>::ReturnValue Reference;3465 typedef typename MapTraits<I n>::ConstReturnValue ConstReference;3459 typedef typename IN::Value Value; 3460 3461 typedef typename MapTraits<IN>::ReferenceMapTag ReferenceMapTag; 3462 typedef typename MapTraits<IN>::ReturnValue ReturnValue; 3463 typedef typename MapTraits<IN>::ConstReturnValue ConstReturnValue; 3464 typedef typename MapTraits<IN>::ReturnValue Reference; 3465 typedef typename MapTraits<IN>::ConstReturnValue ConstReference; 3466 3466 3467 3467 /// Constructor 3468 CombinedNodeMap(I n& in_map, Out& out_map)3468 CombinedNodeMap(IN& in_map, OUT& out_map) 3469 3469 : _in_map(in_map), _out_map(out_map) {} 3470 3470 … … 3498 3498 private: 3499 3499 3500 I n& _in_map;3501 O ut& _out_map;3500 IN& _in_map; 3501 OUT& _out_map; 3502 3502 3503 3503 }; … … 3507 3507 /// 3508 3508 /// This function just returns a combined node map. 3509 template <typename I n, typename Out>3510 static CombinedNodeMap<I n, Out>3511 combinedNodeMap(I n& in_map, Out& out_map) {3512 return CombinedNodeMap<I n, Out>(in_map, out_map);3513 } 3514 3515 template <typename I n, typename Out>3516 static CombinedNodeMap<const I n, Out>3517 combinedNodeMap(const I n& in_map, Out& out_map) {3518 return CombinedNodeMap<const I n, Out>(in_map, out_map);3519 } 3520 3521 template <typename I n, typename Out>3522 static CombinedNodeMap<I n, const Out>3523 combinedNodeMap(I n& in_map, const Out& out_map) {3524 return CombinedNodeMap<I n, const Out>(in_map, out_map);3525 } 3526 3527 template <typename I n, typename Out>3528 static CombinedNodeMap<const I n, const Out>3529 combinedNodeMap(const I n& in_map, const Out& out_map) {3530 return CombinedNodeMap<const I n, const Out>(in_map, out_map);3509 template <typename IN, typename OUT> 3510 static CombinedNodeMap<IN, OUT> 3511 combinedNodeMap(IN& in_map, OUT& out_map) { 3512 return CombinedNodeMap<IN, OUT>(in_map, out_map); 3513 } 3514 3515 template <typename IN, typename OUT> 3516 static CombinedNodeMap<const IN, OUT> 3517 combinedNodeMap(const IN& in_map, OUT& out_map) { 3518 return CombinedNodeMap<const IN, OUT>(in_map, out_map); 3519 } 3520 3521 template <typename IN, typename OUT> 3522 static CombinedNodeMap<IN, const OUT> 3523 combinedNodeMap(IN& in_map, const OUT& out_map) { 3524 return CombinedNodeMap<IN, const OUT>(in_map, out_map); 3525 } 3526 3527 template <typename IN, typename OUT> 3528 static CombinedNodeMap<const IN, const OUT> 3529 combinedNodeMap(const IN& in_map, const OUT& out_map) { 3530 return CombinedNodeMap<const IN, const OUT>(in_map, out_map); 3531 3531 } 3532 3532 -
lemon/arg_parser.cc
r1179 r877 222 222 { 223 223 Opts::iterator o = _opts.find(opt); 224 Opts::iterator s = _opts.find(syn); 224 225 LEMON_ASSERT(o!=_opts.end(), "Unknown option: '"+opt+"'"); 225 LEMON_ASSERT(_opts.find(syn)==_opts.end(), 226 "Option already used: '"+syn+"'"); 226 LEMON_ASSERT(s==_opts.end(), "Option already used: '"+syn+"'"); 227 227 ParData p; 228 228 p.help=opt; -
lemon/lgf_reader.h
r1161 r1150 452 452 /// 453 453 ///\code 454 /// DigraphReader<DGR>(digraph, std::cin) 455 /// .nodeMap("coordinates", coord_map)456 /// .arcMap("capacity", cap_map)457 /// .node("source", src)458 /// .node("target", trg)459 /// .attribute("caption", caption)460 /// .run();454 /// DigraphReader<DGR>(digraph, std::cin). 455 /// nodeMap("coordinates", coord_map). 456 /// arcMap("capacity", cap_map). 457 /// node("source", src). 458 /// node("target", trg). 459 /// attribute("caption", caption). 460 /// run(); 461 461 ///\endcode 462 462 /// … … 1247 1247 ///ListDigraph::ArcMap<int> cap(digraph); 1248 1248 ///ListDigraph::Node src, trg; 1249 ///digraphReader(digraph, std::cin) 1250 /// .arcMap("capacity", cap)1251 /// .node("source", src)1252 /// .node("target", trg)1253 /// .run();1249 ///digraphReader(digraph, std::cin). 1250 /// arcMap("capacity", cap). 1251 /// node("source", src). 1252 /// node("target", trg). 1253 /// run(); 1254 1254 ///\endcode 1255 1255 /// … … 2124 2124 ///ListGraph graph; 2125 2125 ///ListGraph::EdgeMap<int> weight(graph); 2126 ///graphReader(graph, std::cin) 2127 /// .edgeMap("weight", weight)2128 /// .run();2126 ///graphReader(graph, std::cin). 2127 /// edgeMap("weight", weight). 2128 /// run(); 2129 2129 ///\endcode 2130 2130 /// … … 3192 3192 ///ListBpGraph graph; 3193 3193 ///ListBpGraph::EdgeMap<int> weight(graph); 3194 ///bpGraphReader(graph, std::cin) 3195 /// .edgeMap("weight", weight)3196 /// .run();3194 ///bpGraphReader(graph, std::cin). 3195 /// edgeMap("weight", weight). 3196 /// run(); 3197 3197 ///\endcode 3198 3198 /// -
lemon/lgf_writer.h
r1161 r1092 409 409 /// 410 410 ///\code 411 /// DigraphWriter<DGR>(digraph, std::cout) 412 /// .nodeMap("coordinates", coord_map)413 /// .nodeMap("size", size)414 /// .nodeMap("title", title)415 /// .arcMap("capacity", cap_map)416 /// .node("source", src)417 /// .node("target", trg)418 /// .attribute("caption", caption)419 /// .run();411 /// DigraphWriter<DGR>(digraph, std::cout). 412 /// nodeMap("coordinates", coord_map). 413 /// nodeMap("size", size). 414 /// nodeMap("title", title). 415 /// arcMap("capacity", cap_map). 416 /// node("source", src). 417 /// node("target", trg). 418 /// attribute("caption", caption). 419 /// run(); 420 420 ///\endcode 421 421 /// … … 962 962 ///ListDigraph::Node src, trg; 963 963 /// // Setting the capacity map and source and target nodes 964 ///digraphWriter(digraph, std::cout) 965 /// .arcMap("capacity", cap)966 /// .node("source", src)967 /// .node("target", trg)968 /// .run();964 ///digraphWriter(digraph, std::cout). 965 /// arcMap("capacity", cap). 966 /// node("source", src). 967 /// node("target", trg). 968 /// run(); 969 969 ///\endcode 970 970 /// … … 1600 1600 ///ListGraph::EdgeMap<int> weight(graph); 1601 1601 /// // Setting the weight map 1602 ///graphWriter(graph, std::cout) 1603 /// .edgeMap("weight", weight)1604 /// .run();1602 ///graphWriter(graph, std::cout). 1603 /// edgeMap("weight", weight). 1604 /// run(); 1605 1605 ///\endcode 1606 1606 /// … … 2420 2420 ///ListBpGraph::EdgeMap<int> weight(graph); 2421 2421 /// // Setting the weight map 2422 ///bpGraphWriter(graph, std::cout) 2423 /// .edgeMap("weight", weight)2424 /// .run();2422 ///bpGraphWriter(graph, std::cout). 2423 /// edgeMap("weight", weight). 2424 /// run(); 2425 2425 ///\endcode 2426 2426 /// -
lemon/planarity.h
r1182 r1092 2384 2384 if (!pe.run()) return false; 2385 2385 2386 run(pe .embeddingMap());2386 run(pe); 2387 2387 return true; 2388 2388 } … … 2398 2398 void run(const EmbeddingMap& embedding) { 2399 2399 typedef SmartEdgeSet<Graph> AuxGraph; 2400 2401 if (countNodes(_graph) < 3) {2402 int y = 0;2403 for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {2404 _point_map[n].x = 0;2405 _point_map[n].y = y++;2406 }2407 return;2408 }2409 2400 2410 2401 if (3 * countNodes(_graph) - 6 == countEdges(_graph)) { -
lemon/preflow.h
r1169 r1092 477 477 /// flow to the given \c flowMap. The \c flowMap should contain a 478 478 /// flow or at least a preflow, i.e. at each node excluding the 479 /// source node the incoming flow should begreater or equal to the479 /// source node the incoming flow should greater or equal to the 480 480 /// outgoing flow. 481 481 /// \return \c false if the given \c flowMap is not a preflow. … … 496 496 excess -= (*_flow)[e]; 497 497 } 498 if ( _tolerance.negative(excess)&& n != _source) return false;498 if (excess < 0 && n != _source) return false; 499 499 (*_excess)[n] = excess; 500 500 } … … 640 640 (*_excess)[n] = excess; 641 641 642 if ( _tolerance.nonZero(excess)) {642 if (excess != 0) { 643 643 if (new_level + 1 < _level->maxLevel()) { 644 644 _level->liftHighestActive(new_level + 1); … … 721 721 (*_excess)[n] = excess; 722 722 723 if ( _tolerance.nonZero(excess)) {723 if (excess != 0) { 724 724 if (new_level + 1 < _level->maxLevel()) { 725 725 _level->liftActiveOn(level, new_level + 1); … … 792 792 if (!reached[n]) { 793 793 _level->dirtyTopButOne(n); 794 } else if ( _tolerance.positive((*_excess)[n])&& _target != n) {794 } else if ((*_excess)[n] > 0 && _target != n) { 795 795 _level->activate(n); 796 796 } … … 853 853 (*_excess)[n] = excess; 854 854 855 if ( _tolerance.nonZero(excess)) {855 if (excess != 0) { 856 856 if (new_level + 1 < _level->maxLevel()) { 857 857 _level->liftHighestActive(new_level + 1); -
lemon/random.h
r1185 r1136 111 111 static const Word loMask = (1u << 31) - 1; 112 112 static const Word hiMask = ~loMask; 113 113 114 114 115 static Word tempering(Word rnd) { … … 242 243 243 244 private: 245 244 246 245 247 void fillState() { … … 270 272 } 271 273 274 272 275 Word *current; 273 276 Word state[length]; … … 294 297 template <typename Result, typename Word, 295 298 int rest = std::numeric_limits<Result>::digits, int shift = 0, 296 bool last = (rest <= std::numeric_limits<Word>::digits)>299 bool last = rest <= std::numeric_limits<Word>::digits> 297 300 struct IntConversion { 298 301 static const int bits = std::numeric_limits<Word>::digits; … … 340 343 num = rnd() & mask; 341 344 } while (num > max); 342 return static_cast<Result>(num);345 return num; 343 346 } 344 347 }; … … 463 466 }; 464 467 465 /// \ingroup misc 466 /// 467 /// \brief Mersenne Twister random number generator 468 /// 469 /// The Mersenne Twister is a twisted generalized feedback 470 /// shift-register generator of Matsumoto and Nishimura. The period 471 /// of this generator is \f$ 2^{19937} - 1\f$ and it is 472 /// equi-distributed in 623 dimensions for 32-bit numbers. The time 473 /// performance of this generator is comparable to the commonly used 474 /// generators. 475 /// 476 /// This is a template implementation of both 32-bit and 477 /// 64-bit architecture optimized versions. The generators differ 478 /// sligthly in the initialization and generation phase so they 479 /// produce two completly different sequences. 480 /// 481 /// \alert Do not use this class directly, but instead one of \ref 482 /// Random, \ref Random32 or \ref Random64. 483 /// 484 /// The generator gives back random numbers of serveral types. To 485 /// get a random number from a range of a floating point type, you 486 /// can use one form of the \c operator() or the \c real() member 487 /// function. If you want to get random number from the {0, 1, ..., 488 /// n-1} integer range, use the \c operator[] or the \c integer() 489 /// method. And to get random number from the whole range of an 490 /// integer type, you can use the argumentless \c integer() or 491 /// \c uinteger() functions. Finally, you can get random bool with 492 /// equal chance of true and false or with given probability of true 493 /// result using the \c boolean() member functions. 494 /// 495 /// Various non-uniform distributions are also supported: normal (Gauss), 496 /// exponential, gamma, Poisson, etc.; and a few two-dimensional 497 /// distributions, too. 498 /// 499 ///\code 500 /// // The commented code is identical to the other 501 /// double a = rnd(); // [0.0, 1.0) 502 /// // double a = rnd.real(); // [0.0, 1.0) 503 /// double b = rnd(100.0); // [0.0, 100.0) 504 /// // double b = rnd.real(100.0); // [0.0, 100.0) 505 /// double c = rnd(1.0, 2.0); // [1.0, 2.0) 506 /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) 507 /// int d = rnd[100000]; // 0..99999 508 /// // int d = rnd.integer(100000); // 0..99999 509 /// int e = rnd[6] + 1; // 1..6 510 /// // int e = rnd.integer(1, 1 + 6); // 1..6 511 /// int b = rnd.uinteger<int>(); // 0 .. 2^31 - 1 512 /// int c = rnd.integer<int>(); // - 2^31 .. 2^31 - 1 513 /// bool g = rnd.boolean(); // P(g = true) = 0.5 514 /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 515 ///\endcode 516 /// 517 /// LEMON provides a global instance of the random number generator: 518 /// \ref lemon::rnd "rnd". In most cases, it is a good practice 519 /// to use this global generator to get random numbers. 520 /// 521 /// \sa \ref Random, \ref Random32 or \ref Random64. 522 template<class Word> 523 class Random { 524 private: 525 526 _random_bits::RandomCore<Word> core; 527 _random_bits::BoolProducer<Word> bool_producer; 528 529 530 public: 531 532 ///\name Initialization 533 /// 534 /// @{ 535 536 /// \brief Default constructor 537 /// 538 /// Constructor with constant seeding. 539 Random() { core.initState(); } 540 541 /// \brief Constructor with seed 542 /// 543 /// Constructor with seed. The current number type will be converted 544 /// to the architecture word type. 545 template <typename Number> 546 Random(Number seed) { 547 _random_bits::Initializer<Number, Word>::init(core, seed); 548 } 549 550 /// \brief Constructor with array seeding 551 /// 552 /// Constructor with array seeding. The given range should contain 553 /// any number type and the numbers will be converted to the 554 /// architecture word type. 555 template <typename Iterator> 556 Random(Iterator begin, Iterator end) { 557 typedef typename std::iterator_traits<Iterator>::value_type Number; 558 _random_bits::Initializer<Number, Word>::init(core, begin, end); 559 } 560 561 /// \brief Copy constructor 562 /// 563 /// Copy constructor. The generated sequence will be identical to 564 /// the other sequence. It can be used to save the current state 565 /// of the generator and later use it to generate the same 566 /// sequence. 567 Random(const Random& other) { 568 core.copyState(other.core); 569 } 570 571 /// \brief Assign operator 572 /// 573 /// Assign operator. The generated sequence will be identical to 574 /// the other sequence. It can be used to save the current state 575 /// of the generator and later use it to generate the same 576 /// sequence. 577 Random& operator=(const Random& other) { 578 if (&other != this) { 579 core.copyState(other.core); 580 } 581 return *this; 582 } 583 584 /// \brief Seeding random sequence 585 /// 586 /// Seeding the random sequence. The current number type will be 587 /// converted to the architecture word type. 588 template <typename Number> 589 void seed(Number seed) { 590 _random_bits::Initializer<Number, Word>::init(core, seed); 591 } 592 593 /// \brief Seeding random sequence 594 /// 595 /// Seeding the random sequence. The given range should contain 596 /// any number type and the numbers will be converted to the 597 /// architecture word type. 598 template <typename Iterator> 599 void seed(Iterator begin, Iterator end) { 600 typedef typename std::iterator_traits<Iterator>::value_type Number; 601 _random_bits::Initializer<Number, Word>::init(core, begin, end); 602 } 603 604 /// \brief Seeding from file or from process id and time 605 /// 606 /// By default, this function calls the \c seedFromFile() member 607 /// function with the <tt>/dev/urandom</tt> file. If it does not success, 608 /// it uses the \c seedFromTime(). 609 /// \return Currently always \c true. 610 bool seed() { 611 #ifndef LEMON_WIN32 612 if (seedFromFile("/dev/urandom", 0)) return true; 613 #endif 614 if (seedFromTime()) return true; 615 return false; 616 } 617 618 /// \brief Seeding from file 619 /// 620 /// Seeding the random sequence from file. The linux kernel has two 621 /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which 622 /// could give good seed values for pseudo random generators (The 623 /// difference between two devices is that the <tt>random</tt> may 624 /// block the reading operation while the kernel can give good 625 /// source of randomness, while the <tt>urandom</tt> does not 626 /// block the input, but it could give back bytes with worse 627 /// entropy). 628 /// \param file The source file 629 /// \param offset The offset, from the file read. 630 /// \return \c true when the seeding successes. 631 #ifndef LEMON_WIN32 632 bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) 633 #else 634 bool seedFromFile(const std::string& file = "", int offset = 0) 635 #endif 636 { 637 std::ifstream rs(file.c_str()); 638 const int size = 4; 639 Word buf[size]; 640 if (offset != 0 && !rs.seekg(offset)) return false; 641 if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false; 642 seed(buf, buf + size); 643 return true; 644 } 645 646 /// \brief Seeding from process id and time 647 /// 648 /// Seeding from process id and time. This function uses the 649 /// current process id and the current time for initialize the 650 /// random sequence. 651 /// \return Currently always \c true. 652 bool seedFromTime() { 653 #ifndef LEMON_WIN32 654 timeval tv; 655 gettimeofday(&tv, 0); 656 seed(getpid() + tv.tv_sec + tv.tv_usec); 657 #else 658 seed(bits::getWinRndSeed()); 659 #endif 660 return true; 661 } 662 663 /// @} 664 665 ///\name Uniform Distributions 666 /// 667 /// @{ 668 669 /// \brief Returns a random real number from the range [0, 1) 670 /// 671 /// It returns a random real number from the range [0, 1). The 672 /// default Number type is \c double. 673 template <typename Number> 674 Number real() { 675 return _random_bits::RealConversion<Number, Word>::convert(core); 676 } 677 678 double real() { 679 return real<double>(); 680 } 681 682 /// \brief Returns a random real number from the range [0, 1) 683 /// 684 /// It returns a random double from the range [0, 1). 685 double operator()() { 686 return real<double>(); 687 } 688 689 /// \brief Returns a random real number from the range [0, b) 690 /// 691 /// It returns a random real number from the range [0, b). 692 double operator()(double b) { 693 return real<double>() * b; 694 } 695 696 /// \brief Returns a random real number from the range [a, b) 697 /// 698 /// It returns a random real number from the range [a, b). 699 double operator()(double a, double b) { 700 return real<double>() * (b - a) + a; 701 } 702 703 /// \brief Returns a random integer from a range 704 /// 705 /// It returns a random integer from the range {0, 1, ..., b - 1}. 706 template <typename Number> 707 Number integer(Number b) { 708 return _random_bits::Mapping<Number, Word>::map(core, b); 709 } 710 711 /// \brief Returns a random integer from a range 712 /// 713 /// It returns a random integer from the range {a, a + 1, ..., b - 1}. 714 template <typename Number> 715 Number integer(Number a, Number b) { 716 return _random_bits::Mapping<Number, Word>::map(core, b - a) + a; 717 } 718 719 /// \brief Returns a random integer from a range 720 /// 721 /// It returns a random integer from the range {0, 1, ..., b - 1}. 722 template <typename Number> 723 Number operator[](Number b) { 724 return _random_bits::Mapping<Number, Word>::map(core, b); 725 } 726 727 /// \brief Returns a random non-negative integer 728 /// 729 /// It returns a random non-negative integer uniformly from the 730 /// whole range of the current \c Number type. The default result 731 /// type of this function is <tt>unsigned int</tt>. 732 template <typename Number> 733 Number uinteger() { 734 return _random_bits::IntConversion<Number, Word>::convert(core); 735 } 736 737 unsigned int uinteger() { 738 return uinteger<unsigned int>(); 739 } 740 741 /// \brief Returns a random integer 742 /// 743 /// It returns a random integer uniformly from the whole range of 744 /// the current \c Number type. The default result type of this 745 /// function is \c int. 746 template <typename Number> 747 Number integer() { 748 static const int nb = std::numeric_limits<Number>::digits + 749 (std::numeric_limits<Number>::is_signed ? 1 : 0); 750 return _random_bits::IntConversion<Number, Word, nb>::convert(core); 751 } 752 753 int integer() { 754 return integer<int>(); 755 } 756 757 /// \brief Returns a random bool 758 /// 759 /// It returns a random bool. The generator holds a buffer for 760 /// random bits. Every time when it become empty the generator makes 761 /// a new random word and fill the buffer up. 762 bool boolean() { 763 return bool_producer.convert(core); 764 } 765 766 /// @} 767 768 ///\name Non-uniform Distributions 769 /// 770 ///@{ 771 772 /// \brief Returns a random bool with given probability of true result. 773 /// 774 /// It returns a random bool with given probability of true result. 775 bool boolean(double p) { 776 return operator()() < p; 777 } 778 779 /// Standard normal (Gauss) distribution 780 781 /// Standard normal (Gauss) distribution. 782 /// \note The Cartesian form of the Box-Muller 783 /// transformation is used to generate a random normal distribution. 784 double gauss() 785 { 786 double V1,V2,S; 787 do { 788 V1=2*real<double>()-1; 789 V2=2*real<double>()-1; 790 S=V1*V1+V2*V2; 791 } while(S>=1); 792 return std::sqrt(-2*std::log(S)/S)*V1; 793 } 794 /// Normal (Gauss) distribution with given mean and standard deviation 795 796 /// Normal (Gauss) distribution with given mean and standard deviation. 797 /// \sa gauss() 798 double gauss(double mean,double std_dev) 799 { 800 return gauss()*std_dev+mean; 801 } 802 803 /// Lognormal distribution 804 805 /// Lognormal distribution. The parameters are the mean and the standard 806 /// deviation of <tt>exp(X)</tt>. 807 /// 808 double lognormal(double n_mean,double n_std_dev) 809 { 810 return std::exp(gauss(n_mean,n_std_dev)); 811 } 812 /// Lognormal distribution 813 814 /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of 815 /// the mean and the standard deviation of <tt>exp(X)</tt>. 816 /// 817 double lognormal(const std::pair<double,double> ¶ms) 818 { 819 return std::exp(gauss(params.first,params.second)); 820 } 821 /// Compute the lognormal parameters from mean and standard deviation 822 823 /// This function computes the lognormal parameters from mean and 824 /// standard deviation. The return value can direcly be passed to 825 /// lognormal(). 826 std::pair<double,double> lognormalParamsFromMD(double mean, 827 double std_dev) 828 { 829 double fr=std_dev/mean; 830 fr*=fr; 831 double lg=std::log(1+fr); 832 return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg)); 833 } 834 /// Lognormal distribution with given mean and standard deviation 835 836 /// Lognormal distribution with given mean and standard deviation. 837 /// 838 double lognormalMD(double mean,double std_dev) 839 { 840 return lognormal(lognormalParamsFromMD(mean,std_dev)); 841 } 842 843 /// Exponential distribution with given mean 844 845 /// This function generates an exponential distribution random number 846 /// with mean <tt>1/lambda</tt>. 847 /// 848 double exponential(double lambda=1.0) 849 { 850 return -std::log(1.0-real<double>())/lambda; 851 } 852 853 /// Gamma distribution with given integer shape 854 855 /// This function generates a gamma distribution random number. 856 /// 857 ///\param k shape parameter (<tt>k>0</tt> integer) 858 double gamma(int k) 859 { 860 double s = 0; 861 for(int i=0;i<k;i++) s-=std::log(1.0-real<double>()); 862 return s; 863 } 864 865 /// Gamma distribution with given shape and scale parameter 866 867 /// This function generates a gamma distribution random number. 868 /// 869 ///\param k shape parameter (<tt>k>0</tt>) 870 ///\param theta scale parameter 871 /// 872 double gamma(double k,double theta=1.0) 873 { 874 double xi,nu; 875 const double delta = k-std::floor(k); 876 const double v0=E/(E-delta); 877 do { 878 double V0=1.0-real<double>(); 879 double V1=1.0-real<double>(); 880 double V2=1.0-real<double>(); 881 if(V2<=v0) 882 { 883 xi=std::pow(V1,1.0/delta); 884 nu=V0*std::pow(xi,delta-1.0); 885 } 886 else 887 { 888 xi=1.0-std::log(V1); 889 nu=V0*std::exp(-xi); 890 } 891 } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); 892 return theta*(xi+gamma(int(std::floor(k)))); 893 } 894 895 /// Weibull distribution 896 897 /// This function generates a Weibull distribution random number. 898 /// 899 ///\param k shape parameter (<tt>k>0</tt>) 900 ///\param lambda scale parameter (<tt>lambda>0</tt>) 901 /// 902 double weibull(double k,double lambda) 903 { 904 return lambda*pow(-std::log(1.0-real<double>()),1.0/k); 905 } 906 907 /// Pareto distribution 908 909 /// This function generates a Pareto distribution random number. 910 /// 911 ///\param k shape parameter (<tt>k>0</tt>) 912 ///\param x_min location parameter (<tt>x_min>0</tt>) 913 /// 914 double pareto(double k,double x_min) 915 { 916 return exponential(gamma(k,1.0/x_min))+x_min; 917 } 918 919 /// Poisson distribution 920 921 /// This function generates a Poisson distribution random number with 922 /// parameter \c lambda. 923 /// 924 /// The probability mass function of this distribusion is 925 /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] 926 /// \note The algorithm is taken from the book of Donald E. Knuth titled 927 /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the 928 /// return value. 929 930 int poisson(double lambda) 931 { 932 const double l = std::exp(-lambda); 933 int k=0; 934 double p = 1.0; 935 do { 936 k++; 937 p*=real<double>(); 938 } while (p>=l); 939 return k-1; 940 } 941 942 ///@} 943 944 ///\name Two-Dimensional Distributions 945 /// 946 ///@{ 947 948 /// Uniform distribution on the full unit circle 949 950 /// Uniform distribution on the full unit circle. 951 /// 952 dim2::Point<double> disc() 953 { 954 double V1,V2; 955 do { 956 V1=2*real<double>()-1; 957 V2=2*real<double>()-1; 958 959 } while(V1*V1+V2*V2>=1); 960 return dim2::Point<double>(V1,V2); 961 } 962 /// A kind of two-dimensional normal (Gauss) distribution 963 964 /// This function provides a turning symmetric two-dimensional distribution. 965 /// Both coordinates are of standard normal distribution, but they are not 966 /// independent. 967 /// 968 /// \note The coordinates are the two random variables provided by 969 /// the Box-Muller method. 970 dim2::Point<double> gauss2() 971 { 972 double V1,V2,S; 973 do { 974 V1=2*real<double>()-1; 975 V2=2*real<double>()-1; 976 S=V1*V1+V2*V2; 977 } while(S>=1); 978 double W=std::sqrt(-2*std::log(S)/S); 979 return dim2::Point<double>(W*V1,W*V2); 980 } 981 /// A kind of two-dimensional exponential distribution 982 983 /// This function provides a turning symmetric two-dimensional distribution. 984 /// The x-coordinate is of conditionally exponential distribution 985 /// with the condition that x is positive and y=0. If x is negative and 986 /// y=0 then, -x is of exponential distribution. The same is true for the 987 /// y-coordinate. 988 dim2::Point<double> exponential2() 989 { 990 double V1,V2,S; 991 do { 992 V1=2*real<double>()-1; 993 V2=2*real<double>()-1; 994 S=V1*V1+V2*V2; 995 } while(S>=1); 996 double W=-std::log(S)/S; 997 return dim2::Point<double>(W*V1,W*V2); 998 } 999 1000 ///@} 1001 }; 1002 1003 1004 }; 468 } 1005 469 1006 470 /// \ingroup misc … … 1008 472 /// \brief Mersenne Twister random number generator 1009 473 /// 1010 /// Th is class implements either the 32-bit or the 64-bit version of1011 /// the Mersenne Twister random number generator algorithm1012 /// depending on the system architecture.1013 /// 1014 /// For the API description, see its base class1015 /// \ref _random_bits::Random.474 /// The Mersenne Twister is a twisted generalized feedback 475 /// shift-register generator of Matsumoto and Nishimura. The period 476 /// of this generator is \f$ 2^{19937} - 1 \f$ and it is 477 /// equi-distributed in 623 dimensions for 32-bit numbers. The time 478 /// performance of this generator is comparable to the commonly used 479 /// generators. 1016 480 /// 1017 /// \sa \ref _random_bits::Random1018 typedef _random_bits::Random<unsigned long> Random;1019 1020 /// \ingroup misc481 /// This implementation is specialized for both 32-bit and 64-bit 482 /// architectures. The generators differ sligthly in the 483 /// initialization and generation phase so they produce two 484 /// completly different sequences. 1021 485 /// 1022 /// \brief Mersenne Twister random number generator (32-bit version) 486 /// The generator gives back random numbers of serveral types. To 487 /// get a random number from a range of a floating point type you 488 /// can use one form of the \c operator() or the \c real() member 489 /// function. If you want to get random number from the {0, 1, ..., 490 /// n-1} integer range use the \c operator[] or the \c integer() 491 /// method. And to get random number from the whole range of an 492 /// integer type you can use the argumentless \c integer() or \c 493 /// uinteger() functions. After all you can get random bool with 494 /// equal chance of true and false or given probability of true 495 /// result with the \c boolean() member functions. 1023 496 /// 1024 /// This class implements the 32-bit version of the Mersenne Twister 1025 /// random number generator algorithm. It is recommended to be used 1026 /// when someone wants to make sure that the \e same pseudo random 1027 /// sequence will be generated on every platfrom. 497 ///\code 498 /// // The commented code is identical to the other 499 /// double a = rnd(); // [0.0, 1.0) 500 /// // double a = rnd.real(); // [0.0, 1.0) 501 /// double b = rnd(100.0); // [0.0, 100.0) 502 /// // double b = rnd.real(100.0); // [0.0, 100.0) 503 /// double c = rnd(1.0, 2.0); // [1.0, 2.0) 504 /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) 505 /// int d = rnd[100000]; // 0..99999 506 /// // int d = rnd.integer(100000); // 0..99999 507 /// int e = rnd[6] + 1; // 1..6 508 /// // int e = rnd.integer(1, 1 + 6); // 1..6 509 /// int b = rnd.uinteger<int>(); // 0 .. 2^31 - 1 510 /// int c = rnd.integer<int>(); // - 2^31 .. 2^31 - 1 511 /// bool g = rnd.boolean(); // P(g = true) = 0.5 512 /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 513 ///\endcode 1028 514 /// 1029 /// For the API description, see its base class 1030 /// \ref _random_bits::Random. 1031 /// 1032 /// \sa \ref _random_bits::Random 1033 typedef _random_bits::Random<unsigned int> Random32; 1034 1035 /// \ingroup misc 1036 /// 1037 /// \brief Mersenne Twister random number generator (64-bit version) 1038 /// 1039 /// This class implements the 64-bit version of the Mersenne Twister 1040 /// random number generator algorithm. (Even though it runs 1041 /// on 32-bit architectures, too.) It is recommended to be used when 1042 /// someone wants to make sure that the \e same pseudo random sequence 1043 /// will be generated on every platfrom. 1044 /// 1045 /// For the API description, see its base class 1046 /// \ref _random_bits::Random. 1047 /// 1048 /// \sa \ref _random_bits::Random 1049 typedef _random_bits::Random<unsigned long long> Random64; 515 /// LEMON provides a global instance of the random number 516 /// generator which name is \ref lemon::rnd "rnd". Usually it is a 517 /// good programming convenience to use this global generator to get 518 /// random numbers. 519 class Random { 520 private: 521 522 // Architecture word 523 typedef unsigned long Word; 524 525 _random_bits::RandomCore<Word> core; 526 _random_bits::BoolProducer<Word> bool_producer; 527 528 529 public: 530 531 ///\name Initialization 532 /// 533 /// @{ 534 535 /// \brief Default constructor 536 /// 537 /// Constructor with constant seeding. 538 Random() { core.initState(); } 539 540 /// \brief Constructor with seed 541 /// 542 /// Constructor with seed. The current number type will be converted 543 /// to the architecture word type. 544 template <typename Number> 545 Random(Number seed) { 546 _random_bits::Initializer<Number, Word>::init(core, seed); 547 } 548 549 /// \brief Constructor with array seeding 550 /// 551 /// Constructor with array seeding. The given range should contain 552 /// any number type and the numbers will be converted to the 553 /// architecture word type. 554 template <typename Iterator> 555 Random(Iterator begin, Iterator end) { 556 typedef typename std::iterator_traits<Iterator>::value_type Number; 557 _random_bits::Initializer<Number, Word>::init(core, begin, end); 558 } 559 560 /// \brief Copy constructor 561 /// 562 /// Copy constructor. The generated sequence will be identical to 563 /// the other sequence. It can be used to save the current state 564 /// of the generator and later use it to generate the same 565 /// sequence. 566 Random(const Random& other) { 567 core.copyState(other.core); 568 } 569 570 /// \brief Assign operator 571 /// 572 /// Assign operator. The generated sequence will be identical to 573 /// the other sequence. It can be used to save the current state 574 /// of the generator and later use it to generate the same 575 /// sequence. 576 Random& operator=(const Random& other) { 577 if (&other != this) { 578 core.copyState(other.core); 579 } 580 return *this; 581 } 582 583 /// \brief Seeding random sequence 584 /// 585 /// Seeding the random sequence. The current number type will be 586 /// converted to the architecture word type. 587 template <typename Number> 588 void seed(Number seed) { 589 _random_bits::Initializer<Number, Word>::init(core, seed); 590 } 591 592 /// \brief Seeding random sequence 593 /// 594 /// Seeding the random sequence. The given range should contain 595 /// any number type and the numbers will be converted to the 596 /// architecture word type. 597 template <typename Iterator> 598 void seed(Iterator begin, Iterator end) { 599 typedef typename std::iterator_traits<Iterator>::value_type Number; 600 _random_bits::Initializer<Number, Word>::init(core, begin, end); 601 } 602 603 /// \brief Seeding from file or from process id and time 604 /// 605 /// By default, this function calls the \c seedFromFile() member 606 /// function with the <tt>/dev/urandom</tt> file. If it does not success, 607 /// it uses the \c seedFromTime(). 608 /// \return Currently always \c true. 609 bool seed() { 610 #ifndef LEMON_WIN32 611 if (seedFromFile("/dev/urandom", 0)) return true; 612 #endif 613 if (seedFromTime()) return true; 614 return false; 615 } 616 617 /// \brief Seeding from file 618 /// 619 /// Seeding the random sequence from file. The linux kernel has two 620 /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which 621 /// could give good seed values for pseudo random generators (The 622 /// difference between two devices is that the <tt>random</tt> may 623 /// block the reading operation while the kernel can give good 624 /// source of randomness, while the <tt>urandom</tt> does not 625 /// block the input, but it could give back bytes with worse 626 /// entropy). 627 /// \param file The source file 628 /// \param offset The offset, from the file read. 629 /// \return \c true when the seeding successes. 630 #ifndef LEMON_WIN32 631 bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) 632 #else 633 bool seedFromFile(const std::string& file = "", int offset = 0) 634 #endif 635 { 636 std::ifstream rs(file.c_str()); 637 const int size = 4; 638 Word buf[size]; 639 if (offset != 0 && !rs.seekg(offset)) return false; 640 if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false; 641 seed(buf, buf + size); 642 return true; 643 } 644 645 /// \brief Seding from process id and time 646 /// 647 /// Seding from process id and time. This function uses the 648 /// current process id and the current time for initialize the 649 /// random sequence. 650 /// \return Currently always \c true. 651 bool seedFromTime() { 652 #ifndef LEMON_WIN32 653 timeval tv; 654 gettimeofday(&tv, 0); 655 seed(getpid() + tv.tv_sec + tv.tv_usec); 656 #else 657 seed(bits::getWinRndSeed()); 658 #endif 659 return true; 660 } 661 662 /// @} 663 664 ///\name Uniform Distributions 665 /// 666 /// @{ 667 668 /// \brief Returns a random real number from the range [0, 1) 669 /// 670 /// It returns a random real number from the range [0, 1). The 671 /// default Number type is \c double. 672 template <typename Number> 673 Number real() { 674 return _random_bits::RealConversion<Number, Word>::convert(core); 675 } 676 677 double real() { 678 return real<double>(); 679 } 680 681 /// \brief Returns a random real number from the range [0, 1) 682 /// 683 /// It returns a random double from the range [0, 1). 684 double operator()() { 685 return real<double>(); 686 } 687 688 /// \brief Returns a random real number from the range [0, b) 689 /// 690 /// It returns a random real number from the range [0, b). 691 double operator()(double b) { 692 return real<double>() * b; 693 } 694 695 /// \brief Returns a random real number from the range [a, b) 696 /// 697 /// It returns a random real number from the range [a, b). 698 double operator()(double a, double b) { 699 return real<double>() * (b - a) + a; 700 } 701 702 /// \brief Returns a random integer from a range 703 /// 704 /// It returns a random integer from the range {0, 1, ..., b - 1}. 705 template <typename Number> 706 Number integer(Number b) { 707 return _random_bits::Mapping<Number, Word>::map(core, b); 708 } 709 710 /// \brief Returns a random integer from a range 711 /// 712 /// It returns a random integer from the range {a, a + 1, ..., b - 1}. 713 template <typename Number> 714 Number integer(Number a, Number b) { 715 return _random_bits::Mapping<Number, Word>::map(core, b - a) + a; 716 } 717 718 /// \brief Returns a random integer from a range 719 /// 720 /// It returns a random integer from the range {0, 1, ..., b - 1}. 721 template <typename Number> 722 Number operator[](Number b) { 723 return _random_bits::Mapping<Number, Word>::map(core, b); 724 } 725 726 /// \brief Returns a random non-negative integer 727 /// 728 /// It returns a random non-negative integer uniformly from the 729 /// whole range of the current \c Number type. The default result 730 /// type of this function is <tt>unsigned int</tt>. 731 template <typename Number> 732 Number uinteger() { 733 return _random_bits::IntConversion<Number, Word>::convert(core); 734 } 735 736 unsigned int uinteger() { 737 return uinteger<unsigned int>(); 738 } 739 740 /// \brief Returns a random integer 741 /// 742 /// It returns a random integer uniformly from the whole range of 743 /// the current \c Number type. The default result type of this 744 /// function is \c int. 745 template <typename Number> 746 Number integer() { 747 static const int nb = std::numeric_limits<Number>::digits + 748 (std::numeric_limits<Number>::is_signed ? 1 : 0); 749 return _random_bits::IntConversion<Number, Word, nb>::convert(core); 750 } 751 752 int integer() { 753 return integer<int>(); 754 } 755 756 /// \brief Returns a random bool 757 /// 758 /// It returns a random bool. The generator holds a buffer for 759 /// random bits. Every time when it become empty the generator makes 760 /// a new random word and fill the buffer up. 761 bool boolean() { 762 return bool_producer.convert(core); 763 } 764 765 /// @} 766 767 ///\name Non-uniform Distributions 768 /// 769 ///@{ 770 771 /// \brief Returns a random bool with given probability of true result. 772 /// 773 /// It returns a random bool with given probability of true result. 774 bool boolean(double p) { 775 return operator()() < p; 776 } 777 778 /// Standard normal (Gauss) distribution 779 780 /// Standard normal (Gauss) distribution. 781 /// \note The Cartesian form of the Box-Muller 782 /// transformation is used to generate a random normal distribution. 783 double gauss() 784 { 785 double V1,V2,S; 786 do { 787 V1=2*real<double>()-1; 788 V2=2*real<double>()-1; 789 S=V1*V1+V2*V2; 790 } while(S>=1); 791 return std::sqrt(-2*std::log(S)/S)*V1; 792 } 793 /// Normal (Gauss) distribution with given mean and standard deviation 794 795 /// Normal (Gauss) distribution with given mean and standard deviation. 796 /// \sa gauss() 797 double gauss(double mean,double std_dev) 798 { 799 return gauss()*std_dev+mean; 800 } 801 802 /// Lognormal distribution 803 804 /// Lognormal distribution. The parameters are the mean and the standard 805 /// deviation of <tt>exp(X)</tt>. 806 /// 807 double lognormal(double n_mean,double n_std_dev) 808 { 809 return std::exp(gauss(n_mean,n_std_dev)); 810 } 811 /// Lognormal distribution 812 813 /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of 814 /// the mean and the standard deviation of <tt>exp(X)</tt>. 815 /// 816 double lognormal(const std::pair<double,double> ¶ms) 817 { 818 return std::exp(gauss(params.first,params.second)); 819 } 820 /// Compute the lognormal parameters from mean and standard deviation 821 822 /// This function computes the lognormal parameters from mean and 823 /// standard deviation. The return value can direcly be passed to 824 /// lognormal(). 825 std::pair<double,double> lognormalParamsFromMD(double mean, 826 double std_dev) 827 { 828 double fr=std_dev/mean; 829 fr*=fr; 830 double lg=std::log(1+fr); 831 return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg)); 832 } 833 /// Lognormal distribution with given mean and standard deviation 834 835 /// Lognormal distribution with given mean and standard deviation. 836 /// 837 double lognormalMD(double mean,double std_dev) 838 { 839 return lognormal(lognormalParamsFromMD(mean,std_dev)); 840 } 841 842 /// Exponential distribution with given mean 843 844 /// This function generates an exponential distribution random number 845 /// with mean <tt>1/lambda</tt>. 846 /// 847 double exponential(double lambda=1.0) 848 { 849 return -std::log(1.0-real<double>())/lambda; 850 } 851 852 /// Gamma distribution with given integer shape 853 854 /// This function generates a gamma distribution random number. 855 /// 856 ///\param k shape parameter (<tt>k>0</tt> integer) 857 double gamma(int k) 858 { 859 double s = 0; 860 for(int i=0;i<k;i++) s-=std::log(1.0-real<double>()); 861 return s; 862 } 863 864 /// Gamma distribution with given shape and scale parameter 865 866 /// This function generates a gamma distribution random number. 867 /// 868 ///\param k shape parameter (<tt>k>0</tt>) 869 ///\param theta scale parameter 870 /// 871 double gamma(double k,double theta=1.0) 872 { 873 double xi,nu; 874 const double delta = k-std::floor(k); 875 const double v0=E/(E-delta); 876 do { 877 double V0=1.0-real<double>(); 878 double V1=1.0-real<double>(); 879 double V2=1.0-real<double>(); 880 if(V2<=v0) 881 { 882 xi=std::pow(V1,1.0/delta); 883 nu=V0*std::pow(xi,delta-1.0); 884 } 885 else 886 { 887 xi=1.0-std::log(V1); 888 nu=V0*std::exp(-xi); 889 } 890 } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); 891 return theta*(xi+gamma(int(std::floor(k)))); 892 } 893 894 /// Weibull distribution 895 896 /// This function generates a Weibull distribution random number. 897 /// 898 ///\param k shape parameter (<tt>k>0</tt>) 899 ///\param lambda scale parameter (<tt>lambda>0</tt>) 900 /// 901 double weibull(double k,double lambda) 902 { 903 return lambda*pow(-std::log(1.0-real<double>()),1.0/k); 904 } 905 906 /// Pareto distribution 907 908 /// This function generates a Pareto distribution random number. 909 /// 910 ///\param k shape parameter (<tt>k>0</tt>) 911 ///\param x_min location parameter (<tt>x_min>0</tt>) 912 /// 913 double pareto(double k,double x_min) 914 { 915 return exponential(gamma(k,1.0/x_min))+x_min; 916 } 917 918 /// Poisson distribution 919 920 /// This function generates a Poisson distribution random number with 921 /// parameter \c lambda. 922 /// 923 /// The probability mass function of this distribusion is 924 /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] 925 /// \note The algorithm is taken from the book of Donald E. Knuth titled 926 /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the 927 /// return value. 928 929 int poisson(double lambda) 930 { 931 const double l = std::exp(-lambda); 932 int k=0; 933 double p = 1.0; 934 do { 935 k++; 936 p*=real<double>(); 937 } while (p>=l); 938 return k-1; 939 } 940 941 ///@} 942 943 ///\name Two Dimensional Distributions 944 /// 945 ///@{ 946 947 /// Uniform distribution on the full unit circle 948 949 /// Uniform distribution on the full unit circle. 950 /// 951 dim2::Point<double> disc() 952 { 953 double V1,V2; 954 do { 955 V1=2*real<double>()-1; 956 V2=2*real<double>()-1; 957 958 } while(V1*V1+V2*V2>=1); 959 return dim2::Point<double>(V1,V2); 960 } 961 /// A kind of two dimensional normal (Gauss) distribution 962 963 /// This function provides a turning symmetric two-dimensional distribution. 964 /// Both coordinates are of standard normal distribution, but they are not 965 /// independent. 966 /// 967 /// \note The coordinates are the two random variables provided by 968 /// the Box-Muller method. 969 dim2::Point<double> gauss2() 970 { 971 double V1,V2,S; 972 do { 973 V1=2*real<double>()-1; 974 V2=2*real<double>()-1; 975 S=V1*V1+V2*V2; 976 } while(S>=1); 977 double W=std::sqrt(-2*std::log(S)/S); 978 return dim2::Point<double>(W*V1,W*V2); 979 } 980 /// A kind of two dimensional exponential distribution 981 982 /// This function provides a turning symmetric two-dimensional distribution. 983 /// The x-coordinate is of conditionally exponential distribution 984 /// with the condition that x is positive and y=0. If x is negative and 985 /// y=0 then, -x is of exponential distribution. The same is true for the 986 /// y-coordinate. 987 dim2::Point<double> exponential2() 988 { 989 double V1,V2,S; 990 do { 991 V1=2*real<double>()-1; 992 V2=2*real<double>()-1; 993 S=V1*V1+V2*V2; 994 } while(S>=1); 995 double W=-std::log(S)/S; 996 return dim2::Point<double>(W*V1,W*V2); 997 } 998 999 ///@} 1000 }; 1001 1050 1002 1051 1003 extern Random rnd; 1052 1004 1053 1005 } 1054 1006 -
test/bellman_ford_test.cc
r1162 r1131 101 101 pp = const_bf_test.negativeCycle(); 102 102 103 #ifdef LEMON_CXX11104 103 for (BF::ActiveIt it(const_bf_test); it != INVALID; ++it) {} 105 104 for (auto n: const_bf_test.activeNodes()) { ::lemon::ignore_unused_variable_warning(n); } … … 107 106 ::lemon::ignore_unused_variable_warning(n); 108 107 } 109 #endif110 108 } 111 109 { -
test/max_flow_test.cc
r1178 r1092 27 27 #include <lemon/lgf_reader.h> 28 28 #include <lemon/elevator.h> 29 #include <lemon/tolerance.h>30 29 31 30 using namespace lemon; … … 67 66 "target 8\n"; 68 67 69 char test_lgf_float[] =70 "@nodes\n"71 "label\n"72 "0\n"73 "1\n"74 "2\n"75 "3\n"76 "4\n"77 "5\n"78 "6\n"79 "7\n"80 "8\n"81 "9\n"82 "@arcs\n"83 " capacity\n"84 "0 1 0.1\n"85 "0 2 0.1\n"86 "0 3 0.1\n"87 "1 4 0.1\n"88 "2 4 0.1\n"89 "3 4 0.1\n"90 "4 5 0.3\n"91 "5 6 0.1\n"92 "5 7 0.1\n"93 "5 8 0.1\n"94 "6 9 0.1\n"95 "7 9 0.1\n"96 "8 9 0.1\n"97 "@attributes\n"98 "source 0\n"99 "target 9\n";100 68 101 69 // Checks the general interface of a max flow algorithm … … 198 166 typedef concepts::Digraph Digraph; 199 167 typedef concepts::ReadMap<Digraph::Arc, Value> CapMap; 168 typedef Elevator<Digraph, Digraph::Node> Elev; 169 typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev; 200 170 201 171 Digraph g; … … 215 185 216 186 template <typename T> 217 T cutValue (const SmartDigraph& g,218 const SmartDigraph::NodeMap<bool>& cut,219 const SmartDigraph::ArcMap<T>& cap) {220 221 T c =0;222 for (SmartDigraph::ArcIt e(g); e !=INVALID; ++e) {223 if (cut[g.source(e)] && !cut[g.target(e)]) c +=cap[e];187 T cutValue (const SmartDigraph& g, 188 const SmartDigraph::NodeMap<bool>& cut, 189 const SmartDigraph::ArcMap<T>& cap) { 190 191 T c=0; 192 for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) { 193 if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e]; 224 194 } 225 195 return c; … … 230 200 const SmartDigraph::ArcMap<T>& flow, 231 201 const SmartDigraph::ArcMap<T>& cap, 232 SmartDigraph::Node s, SmartDigraph::Node t, 233 const Tolerance<T>& tol) { 202 SmartDigraph::Node s, SmartDigraph::Node t) { 234 203 235 204 for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) { 236 if ( tol.negative(flow[e]) || tol.less(cap[e], flow[e])) return false;205 if (flow[e] < 0 || flow[e] > cap[e]) return false; 237 206 } 238 207 … … 246 215 sum -= flow[e]; 247 216 } 248 if ( tol.nonZero(sum)) return false;217 if (sum != 0) return false; 249 218 } 250 219 return true; 251 220 } 252 221 253 void checkInitPreflow()222 void initFlowTest() 254 223 { 255 224 DIGRAPH_TYPEDEFS(SmartDigraph); 256 225 257 226 SmartDigraph g; 258 SmartDigraph::ArcMap<int> cap(g), 259 Node s = g.addNode(); Node t =g.addNode();260 Node n1 = g.addNode(); Node n2 =g.addNode();227 SmartDigraph::ArcMap<int> cap(g),iflow(g); 228 Node s=g.addNode(); Node t=g.addNode(); 229 Node n1=g.addNode(); Node n2=g.addNode(); 261 230 Arc a; 262 a = g.addArc(s, n1); cap[a] = 20; iflow[a] =20;263 a = g.addArc(n1, n2); cap[a] = 10; iflow[a] =0;264 a = g.addArc(n2, t); cap[a] = 20; iflow[a] =0;265 266 Preflow<SmartDigraph> pre(g, cap, s,t);231 a=g.addArc(s,n1); cap[a]=20; iflow[a]=20; 232 a=g.addArc(n1,n2); cap[a]=10; iflow[a]=0; 233 a=g.addArc(n2,t); cap[a]=20; iflow[a]=0; 234 235 Preflow<SmartDigraph> pre(g,cap,s,t); 267 236 pre.init(iflow); 268 237 pre.startFirstPhase(); 269 270 check(pre.flowValue() == 10, "Incorrect max flow value."); 238 check(pre.flowValue() == 10, "The incorrect max flow value."); 271 239 check(pre.minCut(s), "Wrong min cut (Node s)."); 272 240 check(pre.minCut(n1), "Wrong min cut (Node n1)."); … … 276 244 277 245 template <typename MF, typename SF> 278 void checkMaxFlowAlg( const char *input_lgf, typename MF::Value expected) {246 void checkMaxFlowAlg() { 279 247 typedef SmartDigraph Digraph; 280 248 DIGRAPH_TYPEDEFS(Digraph); … … 285 253 typedef BoolNodeMap CutMap; 286 254 287 Tolerance<Value> tol;288 289 255 Digraph g; 290 256 Node s, t; 291 257 CapMap cap(g); 292 std::istringstream input( input_lgf);293 DigraphReader<Digraph>(g, 258 std::istringstream input(test_lgf); 259 DigraphReader<Digraph>(g,input) 294 260 .arcMap("capacity", cap) 295 .node("source", 296 .node("target", 261 .node("source",s) 262 .node("target",t) 297 263 .run(); 298 264 … … 300 266 max_flow.run(); 301 267 302 check(!tol.different(expected, max_flow.flowValue()), 303 "Incorrect max flow value."); 304 check(checkFlow(g, max_flow.flowMap(), cap, s, t, tol), 268 check(checkFlow(g, max_flow.flowMap(), cap, s, t), 305 269 "The flow is not feasible."); 306 270 … … 309 273 Value min_cut_value = cutValue(g, min_cut, cap); 310 274 311 check( !tol.different(expected, min_cut_value),312 " Incorrectmin cut value.");275 check(max_flow.flowValue() == min_cut_value, 276 "The max flow value is not equal to the min cut value."); 313 277 314 278 FlowMap flow(g); 315 for (ArcIt e(g); e != INVALID; ++e) flow[e] = 13 * max_flow.flowMap()[e]; 316 for (ArcIt e(g); e != INVALID; ++e) cap[e] = 17 * cap[e]; 279 for (ArcIt e(g); e != INVALID; ++e) flow[e] = max_flow.flowMap()[e]; 280 281 Value flow_value = max_flow.flowValue(); 282 283 for (ArcIt e(g); e != INVALID; ++e) cap[e] = 2 * cap[e]; 317 284 max_flow.init(flow); 318 285 … … 323 290 min_cut_value = cutValue(g, min_cut1, cap); 324 291 325 check(!tol.different(17 * expected, max_flow.flowValue()), 326 "Incorrect max flow value."); 327 check(!tol.different(17 * expected, min_cut_value), 328 "Incorrect min cut value."); 292 check(max_flow.flowValue() == min_cut_value && 293 min_cut_value == 2 * flow_value, 294 "The max flow value or the min cut value is wrong."); 329 295 330 296 SF::startSecondPhase(max_flow); // start second phase of the algorithm 331 297 332 check(checkFlow(g, max_flow.flowMap(), cap, s, t , tol),298 check(checkFlow(g, max_flow.flowMap(), cap, s, t), 333 299 "The flow is not feasible."); 334 300 … … 337 303 min_cut_value = cutValue(g, min_cut2, cap); 338 304 339 check( !tol.different(17 * expected, max_flow.flowValue()),340 "Incorrect max flow value.");341 check(!tol.different(17 * expected, min_cut_value),342 "Incorrect min cut value."); 305 check(max_flow.flowValue() == min_cut_value && 306 min_cut_value == 2 * flow_value, 307 "The max flow value or the min cut value was not doubled"); 308 343 309 344 310 max_flow.flowMap(flow); … … 359 325 CutMap min_cut3(g); 360 326 max_flow.minCutMap(min_cut3); 361 min_cut_value =cutValue(g, min_cut3, cap);362 363 check( !tol.different(max_flow.flowValue(), min_cut_value),327 min_cut_value=cutValue(g, min_cut3, cap); 328 329 check(max_flow.flowValue() == min_cut_value, 364 330 "The max flow value or the min cut value is wrong."); 365 331 } … … 414 380 typedef Preflow<SmartDigraph, SmartDigraph::ArcMap<int> > PType1; 415 381 typedef Preflow<SmartDigraph, SmartDigraph::ArcMap<float> > PType2; 416 typedef Preflow<SmartDigraph, SmartDigraph::ArcMap<double> > PType3; 417 418 checkMaxFlowAlg<PType1, PreflowStartFunctions<PType1> >(test_lgf, 13); 419 checkMaxFlowAlg<PType2, PreflowStartFunctions<PType2> >(test_lgf, 13); 420 checkMaxFlowAlg<PType3, PreflowStartFunctions<PType3> >(test_lgf, 13); 421 422 checkMaxFlowAlg<PType2, PreflowStartFunctions<PType2> >(test_lgf_float, 0.3f); 423 checkMaxFlowAlg<PType3, PreflowStartFunctions<PType3> >(test_lgf_float, 0.3); 424 425 checkInitPreflow(); 382 checkMaxFlowAlg<PType1, PreflowStartFunctions<PType1> >(); 383 checkMaxFlowAlg<PType2, PreflowStartFunctions<PType2> >(); 384 initFlowTest(); 426 385 427 386 // Check EdmondsKarp 428 387 typedef EdmondsKarp<SmartDigraph, SmartDigraph::ArcMap<int> > EKType1; 429 388 typedef EdmondsKarp<SmartDigraph, SmartDigraph::ArcMap<float> > EKType2; 430 typedef EdmondsKarp<SmartDigraph, SmartDigraph::ArcMap<double> > EKType3; 431 432 checkMaxFlowAlg<EKType1, GeneralStartFunctions<EKType1> >(test_lgf, 13); 433 checkMaxFlowAlg<EKType2, GeneralStartFunctions<EKType2> >(test_lgf, 13); 434 checkMaxFlowAlg<EKType3, GeneralStartFunctions<EKType3> >(test_lgf, 13); 435 436 checkMaxFlowAlg<EKType2, GeneralStartFunctions<EKType2> >(test_lgf_float, 0.3f); 437 checkMaxFlowAlg<EKType3, GeneralStartFunctions<EKType3> >(test_lgf_float, 0.3); 389 checkMaxFlowAlg<EKType1, GeneralStartFunctions<EKType1> >(); 390 checkMaxFlowAlg<EKType2, GeneralStartFunctions<EKType2> >(); 391 392 initFlowTest(); 438 393 439 394 return 0; -
test/planarity_test.cc
r1182 r798 31 31 using namespace lemon::dim2; 32 32 33 const int lgfn = 8;33 const int lgfn = 4; 34 34 const std::string lgf[lgfn] = { 35 "@nodes\n"36 "label\n"37 "@edges\n"38 " label\n",39 40 "@nodes\n"41 "label\n"42 "0\n"43 "@edges\n"44 " label\n",45 46 "@nodes\n"47 "label\n"48 "0\n"49 "1\n"50 "@edges\n"51 " label\n"52 "0 1 0\n",53 54 "@nodes\n"55 "label\n"56 "0\n"57 "1\n"58 "2\n"59 "@edges\n"60 " label\n"61 "0 1 0\n"62 "1 2 1\n"63 "2 0 2\n",64 65 35 "@nodes\n" 66 36 "label\n" … … 167 137 } 168 138 } 169 170 if (face_num != 0) { 171 check(face_num + countNodes(graph) - countConnectedComponents(graph) == 172 countEdges(graph) + 1, "Euler test does not passed"); 173 } 139 check(face_num + countNodes(graph) - countConnectedComponents(graph) == 140 countEdges(graph) + 1, "Euler test does not passed"); 174 141 } 175 142 … … 279 246 checkEmbedding(graph, pe); 280 247 281 { 282 PlanarDrawing<Graph> pd(graph); 283 pd.run(pe.embeddingMap()); 284 checkDrawing(graph, pd); 285 } 286 287 { 288 PlanarDrawing<Graph> pd(graph); 289 pd.run(); 290 checkDrawing(graph, pd); 291 } 292 293 { 294 PlanarColoring<Graph> pc(graph); 295 pc.runFiveColoring(pe.embeddingMap()); 296 checkColoring(graph, pc, 5); 297 } 298 299 { 300 PlanarColoring<Graph> pc(graph); 301 pc.runFiveColoring(); 302 checkColoring(graph, pc, 5); 303 } 248 PlanarDrawing<Graph> pd(graph); 249 pd.run(pe.embeddingMap()); 250 checkDrawing(graph, pd); 251 252 PlanarColoring<Graph> pc(graph); 253 pc.runFiveColoring(pe.embeddingMap()); 254 checkColoring(graph, pc, 5); 304 255 305 256 } else { -
test/random_test.cc
r1163 r440 22 22 int seed_array[] = {1, 2}; 23 23 24 int rnd_seq32[] = {25 2732, 43567, 42613, 52416, 45891, 21243, 30403, 32103,26 62501, 33003, 12172, 5192, 32511, 50057, 43723, 7813,27 23720, 35343, 6637, 30280, 44566, 31019, 18898, 33867,28 5994, 1688, 11513, 59011, 48056, 25544, 39168, 25365,29 17530, 8366, 27063, 49861, 55169, 63848, 11863, 4960830 };31 int rnd_seq64[] = {32 56382, 63883, 59577, 64750, 9644, 59886, 57647, 18152,33 28520, 64078, 17818, 49294, 26424, 26697, 53684, 19209,34 35404, 12121, 12837, 11827, 32156, 58333, 62553, 7907,35 64427, 39399, 21971, 48789, 46981, 15716, 53335, 65256,36 12999, 15308, 10906, 42162, 47587, 43006, 53921, 1871637 };38 39 void seq_test() {40 for(int i=0;i<5;i++) {41 lemon::Random32 r32(i);42 lemon::Random64 r64(i);43 for(int j=0;j<8;j++) {44 check(r32[65536]==rnd_seq32[i*8+j], "Wrong random sequence");45 check(r64[65536]==rnd_seq64[i*8+j], "Wrong random sequence");46 }47 }48 }49 50 51 24 int main() 52 25 { … … 64 37 (sizeof(seed_array) / sizeof(seed_array[0]))); 65 38 66 seq_test();67 39 return 0; 68 40 } -
tools/dimacs-solver.cc
r1170 r1093 223 223 throw IoError("Cannot open the file for writing", ap.files()[1]); 224 224 } 225 // fall through226 225 case 1: 227 226 input.open(ap.files()[0].c_str()); … … 229 228 throw IoError("File cannot be found", ap.files()[0]); 230 229 } 231 // fall through232 230 case 0: 233 231 break; … … 254 252 case DimacsDescriptor::SP: 255 253 std::cout << "sp"; 256 break;257 254 case DimacsDescriptor::MAT: 258 255 std::cout << "mat"; -
tools/dimacs-to-lgf.cc
r1179 r584 74 74 throw IoError("Cannot open the file for writing", ap.files()[1]); 75 75 } 76 // fall through77 76 case 1: 78 77 input.open(ap.files()[0].c_str()); … … 80 79 throw IoError("File cannot be found", ap.files()[0]); 81 80 } 82 // fall through83 81 case 0: 84 82 break;
Note: See TracChangeset
for help on using the changeset viewer.