7 #include <gsl/gsl_sf_psi.h> 68 template <
typename T,
typename Scalar>
71 size_t max_iter = 1000,
double eps = 1e-50) {
74 long x = X.rows(), y = X.cols();
82 for (
int i = 0; i < x; ++i) {
83 for (
int j = 0; j < y; ++j) {
84 gsl_ran_dirichlet(rand_gen.
get(), z, dirichlet_params.
data(),
85 dirichlet_variates.
data());
86 for (
size_t k = 0; k < z; ++k) {
87 S(i, j, k) = X(i, j) * dirichlet_variates[k];
101 eps_tensor.setConstant(eps);
103 const auto eta = [](
const size_t step) ->
double {
104 return 0.1 /
std::pow(step + 1, 0.55);
109 Eigen::ThreadPoolDevice thread_dev(&tp,
113 for (
size_t eph = 0; eph < max_iter; ++eph) {
116 S_pjk.device(thread_dev) = S.sum(
shape<1>({0}));
117 S_ipk.device(thread_dev) = S.sum(
shape<1>({1}));
123 #pragma omp parallel for schedule(static) 124 for (
size_t k = 0; k < z; ++k) {
125 for (
int i = 0; i < x; ++i) {
126 alpha_eph(i, k) = model_params.
alpha[i] + S_ipk(i, k);
128 for (
int j = 0; j < y; ++j) {
129 beta_eph(j, k) = model_params.
beta[k] + S_pjk(j, k);
134 vector_t<T> alpha_eph_sum = alpha_eph.colwise().sum();
135 #pragma omp parallel for schedule(static) 136 for (
int i = 0; i < x; ++i) {
137 for (
int j = 0; j < y; ++j) {
138 for (
size_t k = 0; k < z; ++k) {
148 Z.device(thread_dev) = S.log() + eps_tensor;
149 S_mult.device(thread_dev) = (S * grad_S).sum(
shape<1>({2}));
151 Z = S.log() + eps_tensor;
152 S_mult = (S * grad_S).sum(
shape<1>({2}));
155 #pragma omp parallel for schedule(static) 156 for (
int i = 0; i < x; ++i) {
157 for (
int j = 0; j < y; ++j) {
158 for (
size_t k = 0; k < z; ++k) {
159 Z(i, j, k) += eta(eph) * S(i, j, k) *
160 (X(i, j) * grad_S(i, j, k) - S_mult(i, j));
165 Z.device(thread_dev) = Z.exp();
172 #pragma omp parallel for schedule(static) 173 for (
int i = 0; i < x; ++i) {
174 for (
int j = 0; j < y; ++j) {
175 for (
size_t k = 0; k < z; ++k) {
176 S(i, j, k) = X(i, j) * Z(i, j, k);
Structure to hold the parameters for the Allocation Model .
Definition: alloc_model_params.hpp:25
Eigen::Tensor< Scalar, N, Eigen::RowMajor > tensor_t
Tensor type used in the computations.
Definition: defs.hpp:52
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > matrix_t
Matrix type used in the computations.
Definition: defs.hpp:41
T hardware_concurrency(T...args)
Eigen::array< size_t, N > shape
Shape of vectors, matrices, tensors, etc.
Definition: defs.hpp:66
void normalize(tensor_t< Scalar, N > &input, size_t axis, NormType type=NormType::L1, double eps=1e-50)
Normalize a tensor in-place along the given axis.
Definition: util.hpp:321
std::vector< Scalar > alpha
Parameter vector of Dirichlet prior for matrix of size .
Definition: alloc_model_params.hpp:37
Real psi_appr(Real x)
Calculate digamma function approximately using the first 8 terms of the asymptotic expansion of digam...
Definition: util.hpp:440
tensor_t< T, 3 > bld_add(const matrix_t< T > &X, size_t z, const alloc_model::Params< Scalar > &model_params, size_t max_iter=1000, double eps=1e-50)
Compute tensor , the solution of BLD problem , from matrix using additive gradient ascent updates...
Definition: bld_add.hpp:69
void check_bld_params(const matrix_t< T > &X, size_t z, const alloc_model::Params< Scalar > &model_params)
Do parameter checks on BLD computing function parameters and throw an assertion error if the paramete...
Definition: util_details.hpp:20
std::vector< Scalar > beta
Parameter vector of Dirichlet prior for matrix of size .
Definition: alloc_model_params.hpp:42
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