6 #include <gsl/gsl_randist.h> 41 template <
typename RandomIterator>
42 RandomIterator
choice(RandomIterator cum_prob_begin,
43 RandomIterator cum_prob_end,
46 if (cum_prob_begin == cum_prob_end) {
51 auto max_val = *(
std::prev(cum_prob_end));
52 double p = gsl_ran_flat(gsl_rng.
get(), 0, max_val);
87 template <
typename Real,
typename Integer>
90 BNMF_ASSERT(prob.cols() == count.cols(),
91 "Number of event probabilities and counts differ in " 92 "util::multinomial_mode");
95 "Number of trials must be nonnegative in util::multinomial_mode");
97 if (prob.cols() == 0) {
102 const auto num_events = prob.cols();
107 normalized_probs = normalized_probs.array() / normalized_probs.sum();
110 (num_trials + 0.5 * num_events) * normalized_probs;
112 count_real = freq.array().floor();
113 diff = freq - count_real;
114 count = count_real.template cast<Integer>();
117 const Integer total_count = count.sum();
118 if (total_count == num_trials) {
127 if (total_count < num_trials) {
128 fractions = (1 - diff.array()) / (count_real.array() + 1);
129 }
else if (total_count > num_trials) {
130 fractions = diff.array() / (count_real.array() + eps);
134 for (
int i = 0; i < num_events; ++i) {
143 return elem_left.second > elem_right.second;
150 for (Integer i = total_count; i < num_trials; ++i) {
153 int min_idx = fraction_heap.
back().first;
158 fraction_heap.
back().second =
159 (1 - diff(min_idx)) / (count(min_idx) + 1);
166 for (Integer i = total_count; i > num_trials; --i) {
169 int min_idx = fraction_heap.
back().first;
174 fraction_heap.
back().second = diff(min_idx) / (count(min_idx) + eps);
204 template <
typename T>
207 BNMF_ASSERT(norm_axis == 0 || norm_axis == 1,
208 "Axis must be 0 or 1 in util::dirichlet_mat");
210 matrix_t<T> result(gamma_shape_mat.rows(), gamma_shape_mat.cols());
214 for (
long i = 0; i < result.rows(); ++i) {
215 for (
long j = 0; j < result.cols(); ++j) {
217 gsl_ran_gamma(rnd_gen.get(), gamma_shape_mat(i, j), 1);
222 if (norm_axis == 0) {
223 result = result.array().rowwise() / result.colwise().sum().array();
225 result = result.array().colwise() / result.rowwise().sum().array();
296 no_repl_comp([](const HeapElem& left, const HeapElem& right) {
298 return left.timestamp >= right.timestamp;
300 rnd_gen(gsl_rng_alloc(gsl_rng_taus), gsl_rng_free) {
302 for (
int i = 0; i < X.rows(); ++i) {
303 for (
int j = 0; j < X.cols(); ++j) {
304 Scalar entry = X(i, j);
307 gsl_ran_beta(rnd_gen.get(), entry + 1, 1);
341 if (m_heap.empty()) {
346 auto& heap_entry = m_heap.back();
351 --heap_entry.matrix_entry;
352 if (heap_entry.matrix_entry <= 0) {
355 heap_entry.timestamp +=
356 gsl_ran_beta(rnd_gen.get(), heap_entry.matrix_entry + 1, 1);
362 prev_val.second = jj;
371 : timestamp(timestamp), matrix_entry(matrix_entry),
441 : cum_prob(), X_cols(X.cols()), X_sum(X.array().sum()),
446 auto X_arr = X.array();
447 cum_prob(0) = X_arr(0);
448 for (
int i = 1; i < X_arr.size(); ++i) {
449 cum_prob(i) = cum_prob(i - 1) + X_arr(i);
476 Scalar* begin = cum_prob.data();
477 Scalar* end = cum_prob.data() + cum_prob.cols();
481 prev_val.first =
static_cast<int>(m / X_cols);
482 prev_val.second =
static_cast<int>(m % X_cols);
494 BNMF_ASSERT(X.array().size() != 0,
495 "Matrix X must have at least one element");
496 BNMF_ASSERT((X.array() >= 0).all(),
"Matrix X must be nonnegative");
498 "Matrix X must have at least one nonzero element");
524 template <
typename T>
529 auto num_samples =
static_cast<size_t>(X.array().sum());
532 util::Generator<std::pair<int, int>,
534 gen(init_val, num_samples + 1,
565 template <
typename T>
573 gen(init_val, num_samples + 1,
SampleOnesNoReplaceComputer(const matrix_t< Scalar > &X)
Construct a new SampleOnesNoReplaceComputer.
Definition: sampling.hpp:294
A template generator that generates values from an initial value by applying a computation to the pre...
Definition: generator.hpp:279
void operator()(size_t curr_step, std::pair< int, int > &prev_val)
Function call operator that will compute the next sample in-place.
Definition: sampling.hpp:474
util::Generator< std::pair< int, int >, details::SampleOnesNoReplaceComputer< T > > sample_ones_noreplace(const matrix_t< T > &X)
Return a bnmf_algs::util::Generator that will generate a sequence of samples from the given matrix wi...
Definition: sampling.hpp:526
SampleOnesReplaceComputer(const matrix_t< Scalar > &X)
Construct a new SampleOnesComputer.
Definition: sampling.hpp:440
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > matrix_t
Matrix type used in the computations.
Definition: defs.hpp:41
RandomIterator choice(RandomIterator cum_prob_begin, RandomIterator cum_prob_end, const bnmf_algs::util::gsl_rng_wrapper &gsl_rng)
Choose a random sample according to the given cumulative probability distribution.
Definition: sampling.hpp:42
util::Generator< std::pair< int, int >, details::SampleOnesReplaceComputer< T > > sample_ones_replace(const matrix_t< T > &X, size_t num_samples)
Return a bnmf_algs::util::Generator that will generate a sequence of samples from the given matrix wi...
Definition: sampling.hpp:567
iter_type begin()
Return the iterator pointing to the previously computed value.
Definition: generator.hpp:421
void check_sample_ones_params(const matrix_t< T > &X)
Definition: sampling.hpp:493
void operator()(size_t curr_step, std::pair< int, int > &prev_val)
Function call operator that will compute the next sample without replacement in-place.
Definition: sampling.hpp:338
void multinomial_mode(Integer num_trials, const vector_t< Real > &prob, vector_t< Integer > &count, double eps=1e-50)
Find the mode of a multinomial distribution using Finucan's algorithm.
Definition: sampling.hpp:88
std::unique_ptr< gsl_rng, decltype(&gsl_rng_free)> gsl_rng_wrapper
RAII wrapper around gsl_rng types.
Definition: defs.hpp:73
Computer type that will compute the next sample from a given matrix with replacement when its operato...
Definition: sampling.hpp:426
Eigen::Matrix< Scalar, 1, Eigen::Dynamic, Eigen::RowMajor > vector_t
Vector type used in the computations.
Definition: defs.hpp:27
Main namespace for bnmf-algs library.
Definition: alloc_model_funcs.hpp:12
Computer type that will compute the next sample from a given matrix without replacement when its oper...
Definition: sampling.hpp:279
matrix_t< T > dirichlet_mat(const matrix_t< T > &gamma_shape_mat, size_t norm_axis)
Construct a matrix containing Dirichlet distribution samples on its rows or columns.
Definition: sampling.hpp:205