util.hh 2.74 KB
Newer Older
1 2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2006-2018 Tiago de Paula Peixoto <tiago@skewed.de>
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#ifndef UTIL_HH
#define UTIL_HH

#include "config.h"

#include <cmath>
#include <iostream>

#include "cache.hh"

28 29
#include <boost/range/counting_range.hpp>

30 31 32 33
namespace graph_tool
{
using namespace boost;

34 35
template <class T1, class T2>
inline double lbinom(T1 N, T2 k)
36
{
37
    if (N == 0 || k == 0 || k >= N)
38
        return 0;
39 40
    assert(N > 0);
    assert(k > 0);
41
    return ((std::lgamma(N + 1) - std::lgamma(k + 1)) - std::lgamma(N - k + 1));
42 43
}

44 45
template <bool Init=true, class T1, class T2>
inline double lbinom_fast(T1 N, T2 k)
46 47 48
{
    if (N == 0 || k == 0 || k > N)
        return 0;
49
    return ((lgamma_fast<Init>(N + 1) - lgamma_fast<Init>(k + 1)) - lgamma_fast<Init>(N - k + 1));
50 51
}

52 53
template <class T1, class T2>
inline double lbinom_careful(T1 N, T2 k)
54 55 56
{
    if (N == 0 || k == 0 || k >= N)
        return 0;
57 58
    double lgN = std::lgamma(N + 1);
    double lgk = std::lgamma(k + 1);
59 60 61 62 63 64 65 66
    if (lgN - lgk > 1e8)
    {
        // We have N >> k. Use Stirling's approximation: ln N! ~ N ln N - N
        // and reorder
        return - N * log1p(-k / N) - k * log1p(-k / N) - k - lgk + k * log(N);
    }
    else
    {
67
        return lgN - std::lgamma(N - k + 1) - lgk;
68 69 70
    }
}

71 72
template <class T>
inline auto lbeta(T x, T y)
73
{
74
    return (std::lgamma(x) + std::lgamma(y)) - std::lgamma(x + y);
75 76
}

77 78 79 80 81 82 83 84
template <class T>
T log_sum(T a, T b)
{
    if (a < b)
        std::swap(a, b);
    return a + std::log1p(exp(b-a));
}

85
template <class Vec, class PosMap, class Val>
86
void remove_element(Vec& vec, PosMap& pos, Val val)
87
{
88 89 90 91 92 93
    auto& back = vec.back();
    auto& j = pos[back];
    auto i = pos[val];
    vec[i] = back;
    j = i;
    vec.pop_back();
94 95
}

96
template <class Vec, class PosMap, class Val>
97
void add_element(Vec& vec, PosMap& pos, Val val)
98
{
99 100
    pos[val] = vec.size();
    vec.push_back(val);
101 102
}

103
template <class Vec, class PosMap, class Val>
104
bool has_element(Vec& vec, PosMap& pos, Val val)
105
{
106
    size_t i = pos[val];
107
    return (i < vec.size() && vec[i] == val);
108 109
}

110 111 112 113

} // namespace graph_tool


114 115 116


#endif // UTIL_HH