1
#include <Rcpp.h>
2
// [[Rcpp::plugins(cpp11)]]
3
// [[Rcpp::depends(RcppEigen)]]
4

5
using namespace Rcpp;
6

7
#include <vector>
8
#include <RcppEigen.h>
9
#include <Rcpp.h>
10

11
// [[Rcpp::export]]
12
std::vector<double> rcpp_rrap_squared_distance(Rcpp::NumericMatrix pu_coordinates, Rcpp::NumericVector pu_probabilities, Rcpp::NumericMatrix dp_coordinates, Rcpp::NumericVector dp_weights, double failure_distance, std::size_t maximum_r_level) {
13
  //// Initialization
14
//   Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> pu_coordinates_MTX(Rcpp::as<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(pu_coordinates));
15
//   Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> dp_coordinates_MTX(Rcpp::as<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(dp_coordinates));
16

17
  Eigen::MatrixXd pu_coordinates_MTX(Rcpp::as<Eigen::MatrixXd>(pu_coordinates));
18
  Eigen::MatrixXd dp_coordinates_MTX(Rcpp::as<Eigen::MatrixXd>(dp_coordinates));
19

20
  //// Preliminary processing
21
  // declare variables
22
  std::size_t N_pu_INT = pu_coordinates_MTX.rows();
23
  std::size_t N_dp_INT = dp_coordinates_MTX.rows(); 
24
  std::size_t N_k_INT = pu_coordinates_MTX.cols(); 
25
  double currProb, curr_value;
26

27
  //// Main processing
28
  std::vector<double> d(N_dp_INT);
29
  Eigen::ArrayXd tmp_AXD(N_pu_INT);
30
  std::vector<size_t> pu_ids(N_pu_INT);
31
  
32
  for (std::size_t j=0; j < N_dp_INT; ++j) {
33
    // init
34
    std::iota(pu_ids.begin(), pu_ids.end(), 0);
35
    
36
    // calculate distances
37
    tmp_AXD = (pu_coordinates_MTX.rowwise() - dp_coordinates_MTX.row(j)).rowwise().squaredNorm().transpose();
38
    
39
    // sort pus in order of distance
40
    std::partial_sort(
41
      pu_ids.begin(), pu_ids.begin()+maximum_r_level, pu_ids.end(),
42
      [&](const std::size_t p1, const std::size_t p2) {
43
        return(tmp_AXD[p1] < tmp_AXD[p2]);
44
      }
45
    );
46
    
47
    // calculate expected weighted distances for real pus
48
    currProb=1.0;
49
    curr_value=0.0;
50
    for (auto i = pu_ids.cbegin(); i<pu_ids.cbegin()+maximum_r_level; ++i) {
51
      curr_value+=(pu_probabilities[*i] * currProb * tmp_AXD[*i] * dp_weights[j]);
52
      currProb*=(1.0 - pu_probabilities[*i]);
53
    }
54
    // calculate expected weighted distance for failure PU
55
    curr_value+=(currProb*failure_distance);
56
    
57
    // update value
58
    d[j] = curr_value;
59
  }
60
  
61
  //// Exports
62
  // return proportion held
63
  return(d);
64
}
65

Read our documentation on viewing source code .

Loading