made rank proposal class

src/augmentation_class.cpp

@@ -2,6 +2,7 @@
 #include "missing_data.h"
 #include "pairwise_comparisons.h"
 #include "distances.h"
+#include "rank_proposal.h"
 using namespace arma;
 
 Augmentation::Augmentation(
@@ -41,8 +42,11 @@ void Augmentation::augment_pairwise(
     vec proposal;
     int g_diff{};
     if(pars.error_model == "none"){
-      proposal = propose_pairwise_augmentation(
+      auto prop = std::make_unique<LeapAndShift>(1);
+      RankProposal rp = prop->propose(
         dat.rankings.col(i), dat.items_above[i], dat.items_below[i]);
+
+      proposal = rp.rankings;
     } else if(pars.error_model == "bernoulli"){
       proposal = propose_swap(dat.rankings.col(i), dat.items_above[i],
                               dat.items_below[i], g_diff, swap_leap);

src/missing_data.h

@@ -19,16 +19,18 @@ arma::mat initialize_missing_ranks(
 
 struct RankProposal{
   RankProposal() {};
+  RankProposal(const arma::vec& rankings) : rankings { rankings } {}
   RankProposal(
     const arma::vec& rankings,
-    double probability,
+    double prob_forward, double prob_backward,
     const arma::uvec& mutated_items) :
-  rankings { rankings }, probability { probability },
-  mutated_items { mutated_items } {}
+  rankings { rankings }, prob_forward { prob_forward },
+  prob_backward { prob_backward }, mutated_items { mutated_items } {}
   ~RankProposal() = default;
 
   arma::vec rankings{};
-  double probability{};
+  double prob_forward{};
+  double prob_backward{};
   arma::uvec mutated_items{};
 };
 

src/missing_data_functions.cpp

@@ -56,7 +56,7 @@ vec make_new_augmentation(const vec& rankings, const uvec& missing_indicator,
   } else {
     pprop = make_pseudo_proposal(
       rankings, missing_indicator, alpha, rho, pseudo_aug_distance);
-    log_hastings_correction = -std::log(pprop.probability) + log_aug_prob;
+    log_hastings_correction = -std::log(pprop.prob_forward) + log_aug_prob;
   }
 
   double u = std::log(R::runif(0, 1));
@@ -68,7 +68,7 @@ vec make_new_augmentation(const vec& rankings, const uvec& missing_indicator,
     log_hastings_correction;
 
   if(ratio > u){
-    log_aug_prob = std::log(pprop.probability);
+    log_aug_prob = std::log(pprop.prob_forward);
     return pprop.rankings;
   } else {
     return rankings;
@@ -81,7 +81,7 @@ RankProposal make_uniform_proposal(const vec& ranks, const uvec& indicator) noex
   vec mutable_ranks = ranks(missing_inds);
   ivec inds = Rcpp::sample(mutable_ranks.size(), mutable_ranks.size()) - 1;
   proposal(missing_inds) = mutable_ranks(conv_to<uvec>::from(inds));
-  return RankProposal(proposal, 1, missing_inds);
+  return RankProposal(proposal, 1, 1, missing_inds);
 }
 
 RankProposal make_pseudo_proposal(
@@ -117,5 +117,5 @@ RankProposal make_pseudo_proposal(
       available_rankings, available_rankings(span(ranking_chosen)));
   }
 
-  return RankProposal(ranks, prob, missing_inds);
+  return RankProposal(ranks, prob, prob, missing_inds);
 }

src/pairwise_comparison_functions.cpp

@@ -1,50 +1,12 @@
 #include <RcppArmadillo.h>
 #include "classes.h"
-#include "leapandshift.h"
 #include "distances.h"
+#include "rank_proposal.h"
 
 using namespace arma;
 
 // [[Rcpp::depends(RcppArmadillo)]]
 
-int find_lower_limit(int item, const uvec& items_above_item, const vec& current_ranking) {
-  if(items_above_item.size() > 0) {
-    return max(current_ranking.elem(items_above_item - 1)) + 1;
-  } else {
-    return 1;
-  }
-}
-
-int find_upper_limit(int item, const uvec& items_below_item, const vec& current_ranking) {
-  if(items_below_item.size() > 0) {
-    return min(current_ranking.elem(items_below_item - 1)) - 1;
-  } else {
-    return current_ranking.size();
-  }
-}
-
-vec propose_pairwise_augmentation(
-    const vec& ranking, const doubly_nested& items_above,
-    const doubly_nested& items_below) {
-  int n_items = ranking.n_elem;
-
-  ivec a = Rcpp::sample(n_items, 1) - 1;
-  int item = a(0);
-
-  int lower_limit = find_lower_limit(item, items_above[item], ranking);
-  int upper_limit = find_upper_limit(item, items_below[item], ranking);
-
-  Rcpp::IntegerVector b = Rcpp::seq(lower_limit, upper_limit);
-  ivec d = Rcpp::sample(b, 1);
-  int proposed_rank = d(0);
-
-  LeapShiftObject ls{ranking};
-  ls.rho_proposal(item) = proposed_rank;
-
-  // Do the shift step
-  ls.shift_step(ranking, item);
-  return ls.rho_proposal;
-}
 
 vec propose_swap(
     const vec& ranking,

src/pairwise_comparisons.h

@@ -2,11 +2,6 @@
 
 #include "classes.h"
 
-arma::vec propose_pairwise_augmentation(
-    const arma::vec& ranking,
-    const doubly_nested& items_above,
-    const doubly_nested & items_below);
-
 arma::vec propose_swap(
     const arma::vec& ranking,
     const doubly_nested& items_above,

src/proposal_functions.cpp

@@ -1,5 +1,5 @@
 #include "classes.h"
-#include "leapandshift.h"
+#include "rank_proposal.h"
 #include "proposal_functions.h"
 using namespace arma;
 
@@ -38,20 +38,23 @@ vec make_new_rho(
 
   int n_items = current_rho.n_elem;
 
-  LeapShiftObject ls = leap_and_shift(current_rho, leap_size, distfun);
+  auto prop = std::make_unique<LeapAndShift>(n_items);
+  RankProposal rp = prop->propose(current_rho, distfun);
 
   double dist_new = arma::sum(
-    distfun->matdist(rankings.rows(ls.indices), ls.rho_proposal(ls.indices)) % observation_frequency
+    distfun->matdist(rankings.rows(rp.mutated_items), rp.rankings(rp.mutated_items)) %
+      observation_frequency
   );
   double dist_old = arma::sum(
-    distfun->matdist(rankings.rows(ls.indices), current_rho(ls.indices)) % observation_frequency
+    distfun->matdist(rankings.rows(rp.mutated_items), current_rho(rp.mutated_items)) %
+      observation_frequency
   );
 
   double ratio = - alpha_old / n_items * (dist_new - dist_old) +
-    std::log(ls.prob_backward) - std::log(ls.prob_forward);
+    std::log(rp.prob_backward) - std::log(rp.prob_forward);
 
   if(ratio > std::log(R::unif_rand())){
-    return(ls.rho_proposal);
+    return(rp.rankings);
   } else {
     return(current_rho);
   }

src/rank_proposal.cpp

@@ -0,0 +1,104 @@
+#include "rank_proposal.h"
+#include "distances.h"
+
+using namespace arma;
+
+LeapAndShift::LeapAndShift(int leap_size) : leap_size { leap_size } {};
+
+int LeapAndShift::find_lower_limit(int item, const uvec& items_above_item,
+                     const vec& current_ranking) {
+  if(items_above_item.size() > 0) {
+    return max(current_ranking.elem(items_above_item - 1)) + 1;
+  } else {
+    return 1;
+  }
+}
+
+int LeapAndShift::find_upper_limit(int item, const uvec& items_below_item, const vec& current_ranking) {
+  if(items_below_item.size() > 0) {
+    return min(current_ranking.elem(items_below_item - 1)) - 1;
+  } else {
+    return current_ranking.size();
+  }
+}
+
+RankProposal LeapAndShift::shift(
+    const RankProposal& rp_in, const vec& current_rank, int u) {
+  RankProposal rp{rp_in};
+  double delta_r = rp.rankings(u) - current_rank(u);
+  rp.mutated_items = zeros<uvec>(std::abs(delta_r) + 1);
+  rp.mutated_items[0] = u;
+
+  if(delta_r > 0){
+    for(int k = 1; k <= delta_r; ++k){
+      int index = as_scalar(find(current_rank == current_rank(u) + k));
+      rp.rankings(index) -= 1;
+      rp.mutated_items[k] = index;
+    }
+  } else if(delta_r < 0) {
+    for(int k = (-1); k >= delta_r; --k){
+      int index = as_scalar(find(current_rank == current_rank(u) + k));
+      rp.rankings(index) += 1;
+      rp.mutated_items[-(k)] = index;
+    }
+  }
+
+  return rp;
+}
+
+RankProposal LeapAndShift::propose(
+    const vec& current_rank,
+    const std::unique_ptr<Distance>& distfun) {
+  RankProposal rp{current_rank};
+  int n_items = current_rank.n_elem;
+
+  ivec a = Rcpp::sample(n_items, 1) - 1;
+  int u = a(0);
+
+  vec support = join_cols(
+    regspace(std::max(1.0, current_rank(u) - leap_size), 1, current_rank(u) - 1),
+    regspace(current_rank(u) + 1, 1, std::min(n_items * 1.0, current_rank(u) + leap_size)));
+
+  ivec b = Rcpp::sample(support.n_elem, 1) - 1;
+  int index = b(0);
+  rp.rankings(u) = support(index);
+
+  double support_new = std::min(rp.rankings(u) - 1, leap_size * 1.0) +
+    std::min(n_items - rp.rankings(u), leap_size * 1.0);
+
+  rp = shift(rp, current_rank, u);
+
+  distfun->update_leap_and_shift_indices(rp.mutated_items, n_items);
+
+  if(std::abs(rp.rankings(u) - current_rank(u)) == 1){
+    rp.prob_forward = 1.0 / (n_items * support.n_elem) + 1.0 / (n_items * support_new);
+    rp.prob_backward = rp.prob_forward;
+  } else {
+    rp.prob_forward = 1.0 / (n_items * support.n_elem);
+    rp.prob_backward = 1.0 / (n_items * support_new);
+  }
+
+  return rp;
+}
+
+RankProposal LeapAndShift::propose(
+    const vec& current_rank, const doubly_nested& items_above,
+    const doubly_nested& items_below
+) {
+  int n_items = current_rank.n_elem;
+
+  ivec a = Rcpp::sample(n_items, 1) - 1;
+  int item = a(0);
+
+  int lower_limit = find_lower_limit(item, items_above[item], current_rank);
+  int upper_limit = find_upper_limit(item, items_below[item], current_rank);
+
+  Rcpp::IntegerVector b = Rcpp::seq(lower_limit, upper_limit);
+  ivec d = Rcpp::sample(b, 1);
+  int proposed_rank = d(0);
+
+  RankProposal rp{current_rank};
+  rp.rankings(item) = proposed_rank;
+  rp = shift(rp, current_rank, item);
+  return rp;
+}