Skip to content

Commit

Permalink
Removal of arma:: qualifier from functions (#205)
Browse files Browse the repository at this point in the history
* replaced arma::*class* in mixtures.cpp

* removed arma:: in distances.cpp

* replaced arma:: in importance_sampling.cpp

* leapandshift

* removed more arma:: in misc.cpp

* removed arma::*class* in missing_data.cpp

* removing arma:: in pairwise_comparisons.cpp

* removing arma:: in parameterupdates.cpp

* removing arma:: in partitionfuns.cpp

* removing arma:: in rmallows.cpp

* removing arma:: in run_mcmc.cpp

* removed arma::*class* from smc functions

* replaced arma::ones with ones

* replaced arma::zeros and arma::rand* with zeros and rand*

* replaced arma::fill* with fill*

* replaced arma::span with span

* replaced arma::exp with exp

* replaced arma::log with log

* replaced arma::accu with accu

* replaced arma::conv_to with conv_to

* replaced arma::max and arma::min with max and min

* replaced arma::sum with sum

* replaced arma::normalise with normalise

* replaced arma::find with find

* replaced arma::regspace and arma::linspace with regspace and linspace

* replaced arma::shuffle with shuffle

* replaced arma::as_scalar with as_scalar

* removed arma:: from various not-so-frequently-used functions

* removal of arma:: qualifier from more functions
  • Loading branch information
osorensen authored Apr 27, 2022
1 parent d1f4654 commit 3ec2432
Show file tree
Hide file tree
Showing 26 changed files with 490 additions and 464 deletions.
36 changes: 19 additions & 17 deletions src/distances.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2,36 +2,38 @@
#include "subset.h"
#include "distances.h"

using namespace arma;

// [[Rcpp::depends(RcppArmadillo)]]

double cayley_distance(const arma::vec& r1, const arma::vec& r2){
double cayley_distance(const vec& r1, const vec& r2){
double distance = 0;
int n = r1.n_elem;
double tmp1;
arma::vec tmp2 = r1;
vec tmp2 = r1;

// This is a C++ translation of Rankcluster::distCayley
for(int i = 0; i < n; ++i){
if(tmp2(i) != r2(i)) {
distance += 1;
tmp1 = tmp2(i);
tmp2(i) = r2(i);
arma::uvec inds = arma::find(tmp2 == r2(i));
uvec inds = find(tmp2 == r2(i));
tmp2.elem(inds).fill(tmp1);
}
}
return distance;
}

double footrule_distance(const arma::vec& r1, const arma::vec& r2){
return arma::norm(r1 - r2, 1);
double footrule_distance(const vec& r1, const vec& r2){
return norm(r1 - r2, 1);
}

double hamming_distance(const arma::vec& r1, const arma::vec& r2){
return arma::sum(r1 != r2);
double hamming_distance(const vec& r1, const vec& r2){
return sum(r1 != r2);
}

double kendall_distance(const arma::vec& r1, const arma::vec& r2){
double kendall_distance(const vec& r1, const vec& r2){
double distance = 0;
int n = r1.n_elem;

Expand All @@ -46,25 +48,25 @@ double kendall_distance(const arma::vec& r1, const arma::vec& r2){
return distance;
}

double spearman_distance(const arma::vec& r1, const arma::vec& r2){
return std::pow(arma::norm(r1 - r2, 2), 2.0);
double spearman_distance(const vec& r1, const vec& r2){
return std::pow(norm(r1 - r2, 2), 2.0);
}

double ulam_distance(const arma::vec& r1, const arma::vec& r2){
double ulam_distance(const vec& r1, const vec& r2){

int N = r1.n_elem;

arma::ivec a = arma::conv_to<arma::ivec>::from(r1);
arma::ivec b = arma::conv_to<arma::ivec>::from(r2);
ivec a = conv_to<ivec>::from(r1);
ivec b = conv_to<ivec>::from(r2);

int *p1 = (int*) calloc(N, sizeof (int));
int *p2 = (int*) calloc(N, sizeof (int));

int distance;

for(int i = 0; i < N; ++i){
p1[i] = static_cast<int>(arma::as_scalar(a(i)) - 1);
p2[i] = static_cast<int>(arma::as_scalar(b(i)) - 1);
p1[i] = static_cast<int>(as_scalar(a(i)) - 1);
p2[i] = static_cast<int>(as_scalar(b(i)) - 1);
}

distance = perm0_distance ( N, p1, p2 );
Expand Down Expand Up @@ -122,7 +124,7 @@ double get_rank_distance(arma::vec r1, arma::vec r2, std::string metric){
// [[Rcpp::export]]
double rank_dist_sum(const arma::mat& rankings, const arma::vec& rho,
const std::string& metric, const arma::vec& obs_freq){
return arma::sum(rank_dist_vec(rankings, rho, metric, obs_freq));
return sum(rank_dist_vec(rankings, rho, metric, obs_freq));
}


Expand All @@ -134,7 +136,7 @@ arma::vec rank_dist_vec(const arma::mat& rankings,
const arma::vec& obs_freq){

int n = rankings.n_cols;
arma::vec result = arma::zeros(n);
vec result = zeros(n);

for(int i = 0; i < n; ++i){
result(i) = get_rank_distance(rankings.col(i), rho, metric) * obs_freq(i);
Expand Down
40 changes: 21 additions & 19 deletions src/importance_sampling.cpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
#include "RcppArmadillo.h"
#include "distances.h"

using namespace arma;

// via the depends attribute we tell Rcpp to create hooks for
// RcppArmadillo so that the build process will know what to do
//
Expand All @@ -26,57 +28,57 @@ arma::vec compute_importance_sampling_estimate(arma::vec alpha_vector, int n_ite
int n_alphas = alpha_vector.n_elem;

// The reference ranking
arma::vec rho = arma::regspace<arma::vec>(1, n_items);
vec rho = regspace<vec>(1, n_items);

// Vector which holds the result for all alphas
arma::vec logZ = arma::zeros(n_alphas);
vec logZ = zeros(n_alphas);

// Loop over the values of alpha
for(int t = 0; t < n_alphas; ++t){
// The current value of alpha
double alpha = alpha_vector(t);
// The current value of the partition function
arma::vec partfun(nmc);
vec partfun(nmc);

// Loop over the Monte Carlo samples
for(int i = 0; i < nmc; ++i){
// Support set of the proposal distribution
arma::vec support = arma::regspace<arma::vec>(1, n_items);
vec support = regspace<vec>(1, n_items);

// Vector which holds the proposed ranks
arma::vec ranks = arma::zeros(n_items);
arma::vec ranks2 = arma::zeros(n_items);
vec ranks = zeros(n_items);
vec ranks2 = zeros(n_items);

// Probability of the ranks we get
double log_q = 0;

// n_items random uniform numbers
arma::vec u = arma::log(arma::randu(n_items));
vec u = log(randu(n_items));

// Loop over possible values given to item j in random order
arma::vec myind = arma::shuffle(arma::regspace(0, n_items - 1));
vec myind = shuffle(regspace(0, n_items - 1));

for(int j = 0; j < n_items; ++j){
int jj = myind(j);
// Find the elements that have not been taken yet
arma::uvec inds = arma::find(support != 0);
arma::vec log_prob(inds.size());
uvec inds = find(support != 0);
vec log_prob(inds.size());

// Number of elements
int k_max = inds.n_elem;

// Reference vector
arma::vec r1 = inds + arma::ones(k_max);
vec r1 = inds + ones(k_max);
// Sampled vector
arma::vec r2 = rho(jj) * arma::ones(k_max);
vec r2 = rho(jj) * ones(k_max);
// Probability of sample. Note that this is a vector quantity.
log_prob = - alpha / n_items * arma::pow(arma::abs(r1 - r2), (metric == "footrule") ? 1. : 2.);
log_prob = log_prob - std::log(arma::accu(arma::exp(log_prob)));
arma::vec log_cpd = arma::log(arma::cumsum(arma::exp(log_prob)));
log_prob = - alpha / n_items * pow(abs(r1 - r2), (metric == "footrule") ? 1. : 2.);
log_prob = log_prob - std::log(accu(exp(log_prob)));
vec log_cpd = log(cumsum(exp(log_prob)));

// Draw a random sample
int item_index = arma::as_scalar(arma::find(log_cpd > u(jj), 1));
ranks(jj) = arma::as_scalar(inds(item_index)) + 1;
int item_index = as_scalar(find(log_cpd > u(jj), 1));
ranks(jj) = as_scalar(inds(item_index)) + 1;

log_q += log_prob(item_index);

Expand All @@ -89,8 +91,8 @@ arma::vec compute_importance_sampling_estimate(arma::vec alpha_vector, int n_ite
}
// Average over the Monte Carlo samples
// Using this trick: https://www.xarg.org/2016/06/the-log-sum-exp-trick-in-machine-learning/
double maxval = arma::max(partfun);
logZ(t) = maxval + std::log(arma::accu(arma::exp(partfun - maxval))) - std::log(static_cast<double>(nmc));
double maxval = max(partfun);
logZ(t) = maxval + std::log(accu(exp(partfun - maxval))) - std::log(static_cast<double>(nmc));
}
return logZ;
}
31 changes: 16 additions & 15 deletions src/leapandshift.cpp
Original file line number Diff line number Diff line change
@@ -1,40 +1,41 @@
#include "RcppArmadillo.h"
using namespace arma;

// [[Rcpp::depends(RcppArmadillo)]]

void shift_step(arma::vec& rho_proposal, const arma::vec& rho,
const int& u, double& delta_r, arma::uvec& indices){
void shift_step(vec& rho_proposal, const vec& rho,
const int& u, double& delta_r, uvec& indices){
// Shift step:
delta_r = rho_proposal(u - 1) - rho(u - 1);
indices = arma::zeros<arma::uvec>(std::abs(delta_r) + 1);
indices = zeros<uvec>(std::abs(delta_r) + 1);
indices[0] = u-1;
int index;

if(delta_r > 0){
for(int k = 1; k <= delta_r; ++k){
index = arma::as_scalar(arma::find(rho == rho(u-1) + k));
index = as_scalar(find(rho == rho(u-1) + k));
rho_proposal(index) -= 1;
indices[k] = index;
}
} else if(delta_r < 0) {
for(int k = (-1); k >= delta_r; --k){
index = arma::as_scalar(arma::find(rho == rho(u-1) + k));
index = as_scalar(find(rho == rho(u-1) + k));
rho_proposal(index) += 1;
indices[-(k)] = index;
}
}
}


void leap_and_shift(arma::vec& rho_proposal, arma::uvec& indices,
void leap_and_shift(vec& rho_proposal, uvec& indices,
double& prob_backward, double& prob_forward,
const arma::vec& rho, int leap_size, bool reduce_indices){
const vec& rho, int leap_size, bool reduce_indices){

// Set proposal equal to current
rho_proposal = rho;

// Help vectors
arma::vec support;
vec support;

// Number of items
int n = rho.n_elem;
Expand All @@ -45,7 +46,7 @@ void leap_and_shift(arma::vec& rho_proposal, arma::uvec& indices,

// Leap 1
// 1, sample u randomly between 1 and n
u = arma::as_scalar(arma::randi(1, arma::distr_param(1, n)));
u = as_scalar(randi(1, distr_param(1, n)));

// 2, compute the set S for sampling the new rank
double dobL = static_cast<double>(leap_size);
Expand All @@ -56,17 +57,17 @@ void leap_and_shift(arma::vec& rho_proposal, arma::uvec& indices,
double length2 = std::min(n - rho(u - 1), dobL);

if ((rho(u - 1) > 1) && (rho(u - 1) < n)) {
support = arma::join_cols(arma::linspace(
support = join_cols(linspace(
std::max(1.0, rho(u - 1) - leap_size), rho(u - 1) - 1, length1),
arma::linspace(rho(u - 1) + 1, std::min(dobn, rho(u - 1) + leap_size), length2));
linspace(rho(u - 1) + 1, std::min(dobn, rho(u - 1) + leap_size), length2));
} else if(rho(u - 1) == 1){
support = arma::linspace(rho(u - 1) + 1, std::min(dobn, rho(u - 1) + leap_size), length2);
support = linspace(rho(u - 1) + 1, std::min(dobn, rho(u - 1) + leap_size), length2);
} else if(rho(u - 1) == n){
support = arma::linspace(std::max(1.0, rho(u - 1) - leap_size), rho(u - 1) - 1, length1);
support = linspace(std::max(1.0, rho(u - 1) - leap_size), rho(u - 1) - 1, length1);
}

// 3. assign a random element of the support set, this completes the leap step
index = arma::as_scalar(arma::randi(1, arma::distr_param(0, support.n_elem-1)));
index = as_scalar(randi(1, distr_param(0, support.n_elem-1)));
// Picked element index-1 from the support set
rho_proposal(u-1) = support(index);

Expand All @@ -87,6 +88,6 @@ void leap_and_shift(arma::vec& rho_proposal, arma::uvec& indices,
shift_step(rho_proposal, rho, u, delta_r, indices);

if(!reduce_indices){
indices = arma::regspace<arma::uvec>(0, n - 1);
indices = regspace<uvec>(0, n - 1);
}
}
Loading

0 comments on commit 3ec2432

Please sign in to comment.