util.hh 2.76 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-2020 Tiago de Paula Peixoto <tiago@skewed.de>
4
//
5 6 7 8
// This program is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 3 of the License, or (at your option) any
// later version.
9
//
10 11 12 13
// 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 Lesser General Public License for more
// details.
14
//
15
// You should have received a copy of the GNU Lesser General Public License
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#ifndef UTIL_HH
#define UTIL_HH

#include "config.h"

#include <cmath>

#include "cache.hh"

namespace graph_tool
{
using namespace boost;

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

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

51
template <class T1, class T2>
52
[[gnu::const]]
53
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
template <class T>
72
[[gnu::const]]
73
inline auto lbeta(T x, T y)
74
{
75
    return (std::lgamma(x) + std::lgamma(y)) - std::lgamma(x + y);
76 77
}

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

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

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

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

111 112 113 114

} // namespace graph_tool


115 116 117


#endif // UTIL_HH