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 0
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 0
  Eigen::MatrixXd pu_coordinates_MTX(Rcpp::as<Eigen::MatrixXd>(pu_coordinates));
18 0
  Eigen::MatrixXd dp_coordinates_MTX(Rcpp::as<Eigen::MatrixXd>(dp_coordinates));
19

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

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

Read our documentation on viewing source code .

Loading