C++ NumericMatrix::begin方法代码示例

// [[Rcpp::export]]
SEXP FilterIndepOld  (Rcpp::NumericVector y_rcpp,
					Rcpp::NumericMatrix y_lagged_rcpp,
					Rcpp::NumericMatrix z_dependent_rcpp,
					Rcpp::NumericMatrix z_independent_rcpp,
					Rcpp::NumericVector beta_rcpp,
					Rcpp::NumericVector mu_rcpp,
					Rcpp::NumericVector sigma_rcpp,
					Rcpp::NumericMatrix gamma_dependent_rcpp,
					Rcpp::NumericVector gamma_independent_rcpp,
					Rcpp::NumericMatrix transition_probs_rcpp,
					Rcpp::NumericVector initial_dist_rcpp
					)
{
	int n = y_rcpp.size();
	int M = mu_rcpp.size();
	arma::mat xi_k_t(M, n); // make a transpose first for easier column operations.

	arma::colvec y(y_rcpp.begin(), y_rcpp.size(), false);
	arma::mat    y_lagged(y_lagged_rcpp.begin(),
								y_lagged_rcpp.nrow(), y_lagged_rcpp.ncol(), false);
	arma::mat    z_dependent(z_dependent_rcpp.begin(),
								z_dependent_rcpp.nrow(), z_dependent_rcpp.ncol(), false);
	arma::mat    z_independent(z_independent_rcpp.begin(),
								z_independent_rcpp.nrow(), z_independent_rcpp.ncol(), false);
	arma::colvec beta(beta_rcpp.begin(), beta_rcpp.size(), false);
	arma::colvec mu(mu_rcpp.begin(), mu_rcpp.size(), false);
	arma::colvec sigma(sigma_rcpp.begin(), sigma_rcpp.size(), false);
	arma::mat    gamma_dependent(gamma_dependent_rcpp.begin(),
								gamma_dependent_rcpp.nrow(), gamma_dependent_rcpp.ncol(), false);
	arma::colvec gamma_independent(gamma_independent_rcpp.begin(),
								gamma_independent_rcpp.size(), false);
	arma::mat    transition_probs(transition_probs_rcpp.begin(),
								transition_probs_rcpp.nrow(), transition_probs_rcpp.ncol(), false);
	arma::colvec initial_dist(initial_dist_rcpp.begin(), initial_dist_rcpp.size(), false);

	double likelihood = 0;

	SEXP eta_rcpp = EtaIndep(y_rcpp, y_lagged_rcpp,
						z_dependent_rcpp, z_independent_rcpp,
						beta_rcpp, mu_rcpp, sigma_rcpp,
						gamma_dependent_rcpp, gamma_independent_rcpp);
	arma::mat eta_t = (Rcpp::as<arma::mat>(eta_rcpp)).t();

	xi_k_t.col(0) = eta_t.col(0) % initial_dist;
	double total = sum(xi_k_t.col(0));
	xi_k_t.col(0) = xi_k_t.col(0) / total;
	likelihood += log(total);

	for (int k = 1; k < n; k++)
	{
	  	xi_k_t.col(k) = eta_t.col(k) % (transition_probs * xi_k_t.col(k-1));
	 	total = sum(xi_k_t.col(k));
		xi_k_t.col(k) = xi_k_t.col(k) / total;
		likelihood += log(total);
	}

	return Rcpp::List::create(Named("xi.k") = wrap(xi_k_t.t()),
														Named("likelihood") = wrap(likelihood));
}

    double sPredModule::solveCoef(double wrss) {
	Rcpp::NumericMatrix
	    cc = d_fac.solve(CHOLMOD_L, d_fac.solve(CHOLMOD_P, d_Vtr));
	double ans = sqrt(inner_product(cc.begin(), cc.end(),
					cc.begin(), double())/wrss);
	Rcpp::NumericMatrix
	    bb = d_fac.solve(CHOLMOD_Pt, d_fac.solve(CHOLMOD_Lt, cc));
	copy(bb.begin(), bb.end(), d_coef.begin());
	return ans;
    }

Rcpp::List RadiusSearch(Rcpp::NumericMatrix query_,
                        Rcpp::NumericMatrix ref_,
                        double radius,
                        int max_neighbour,
                        std::string build,
                        int cores,
                        int checks) {
  const std::size_t n_dim = query_.ncol();
  const std::size_t n_query = query_.nrow();
  const std::size_t n_ref = ref_.nrow();
  // Column major to row major
  arma::mat query(n_dim, n_query);
  {
    arma::mat temp_q(query_.begin(), n_query, n_dim, false);
    query = arma::trans(temp_q);
  }
  flann::Matrix<double> q_flann(query.memptr(), n_query, n_dim);
  arma::mat ref(n_dim, n_ref);
  {
    arma::mat temp_r(ref_.begin(), n_ref, n_dim, false);
    ref = arma::trans(temp_r);
  }
  flann::Matrix<double> ref_flann(ref.memptr(), n_ref, n_dim);
  // Setting the flann index params
  flann::IndexParams params;
  if (build == "kdtree") {
    params = flann::KDTreeSingleIndexParams(1);
  } else if (build == "kmeans") {
    params = flann::KMeansIndexParams(2, 10, flann::FLANN_CENTERS_RANDOM, 0.2);
  } else if (build == "linear") {
    params = flann::LinearIndexParams();
  }
  // Perform the radius search
  flann::Index<flann::L2<double> > index(ref_flann, params);
  index.buildIndex();
  std::vector< std::vector<int> >
      indices_flann(n_query, std::vector<int>(max_neighbour));
  std::vector< std::vector<double> >
      dists_flann(n_query, std::vector<double>(max_neighbour));
  flann::SearchParams search_params;
  search_params.cores = cores;
  search_params.checks = checks;
  search_params.max_neighbors = max_neighbour;
  index.radiusSearch(q_flann, indices_flann, dists_flann, radius,
                     search_params);
  return Rcpp::List::create(Rcpp::Named("indices") = indices_flann,
                            Rcpp::Named("distances") = dists_flann);
}

//' Check whether there are any non-finite values in a matrix
//'
//' The C++ functions will not work with NA values, and the calculation of the
//' summary profile will take a long time to run before crashing.
//'
//' @param matPtr matrix to check.
//' 
//' @return
//'  Throws an error if any \code{NA}, \code{NaN}, \code{Inf}, or \code{-Inf}
//'  values are found, otherwise returns silently.
//' 
// [[Rcpp::export]]
void CheckFinite(Rcpp::NumericMatrix matPtr) {
  arma::mat mat = arma::mat(matPtr.begin(), matPtr.nrow(), matPtr.ncol(), false, true);
  arma::uvec nonFiniteIdx = arma::find_nonfinite(mat);
  if (nonFiniteIdx.n_elem > 0) {
    throw Rcpp::exception("matrices cannot have non-finite or missing values");
  } 
}

void eigs_sym_shift_c(
    mat_op op, int n, int k, double sigma,
    const spectra_opts *opts, void *data,
    int *nconv, int *niter, int *nops,
    double *evals, double *evecs, int *info
)
{
    BEGIN_RCPP

    CRealShift cmat_op(op, n, data);
    Rcpp::List res;
    try {
        res = run_eigs_shift_sym((RealShift*) &cmat_op, n, k, opts->ncv, opts->rule,
                                 sigma, opts->maxitr, opts->tol, opts->retvec != 0);
        *info = 0;
    } catch(...) {
        *info = 1;  // indicates error
    }

    *nconv = Rcpp::as<int>(res["nconv"]);
    *niter = Rcpp::as<int>(res["niter"]);
    *nops  = Rcpp::as<int>(res["nops"]);
    Rcpp::NumericVector val = res["values"];
    std::copy(val.begin(), val.end(), evals);
    if(opts->retvec != 0)
    {
        Rcpp::NumericMatrix vec = res["vectors"];
        std::copy(vec.begin(), vec.end(), evecs);
    }

    VOID_END_RCPP
}

void ScoreGaussL0PenScatter::setData(Rcpp::List& data)
{
	std::vector<int>::iterator vi;
	//uint i;

	// Cast preprocessed data from R list
	dout.level(2) << "Casting preprocessed data...\n";
	_dataCount = Rcpp::as<std::vector<int> >(data["data.count"]);
	dout.level(3) << "# samples per vertex: " << _dataCount << "\n";
	_totalDataCount = Rcpp::as<uint>(data["total.data.count"]);
	dout.level(3) << "Total # samples: " << _totalDataCount << "\n";
	Rcpp::List scatter = data["scatter"];
	Rcpp::NumericMatrix scatterMat;
	_disjointScatterMatrices.resize(scatter.size());
	dout.level(3) << "# disjoint scatter matrices: " << scatter.size() << "\n";
	for (R_len_t i = 0; i < scatter.size(); ++i) {
		scatterMat = Rcpp::NumericMatrix((SEXP)(scatter[i]));
		_disjointScatterMatrices[i] = arma::mat(scatterMat.begin(), scatterMat.nrow(), scatterMat.ncol(), false);
	}

	// Cast index of scatter matrices, adjust R indexing convention to C++
	std::vector<int> scatterIndex = Rcpp::as<std::vector<int> >(data["scatter.index"]);
	for (std::size_t i = 0; i < scatterIndex.size(); ++i)
		_scatterMatrices[i] = &(_disjointScatterMatrices[scatterIndex[i] - 1]);

	// Cast lambda: penalty constant
	_lambda = Rcpp::as<double>(data["lambda"]);
	dout.level(3) << "Penalty parameter lambda: " << _lambda << "\n";

	// Check whether an intercept should be calculated
	_allowIntercept = Rcpp::as<bool>(data["intercept"]);
	dout.level(3) << "Include intercept: " << _allowIntercept << "\n";
}

// [[Rcpp::export]]
SEXP EtaIndep  (Rcpp::NumericVector y_rcpp,
				Rcpp::NumericMatrix y_lagged_rcpp,
				Rcpp::NumericMatrix z_dependent_rcpp,
				Rcpp::NumericMatrix z_independent_rcpp,
				Rcpp::NumericVector beta_rcpp,
				Rcpp::NumericVector mu_rcpp,
				Rcpp::NumericVector sigma_rcpp,
				Rcpp::NumericMatrix gamma_dependent_rcpp,
				Rcpp::NumericVector gamma_independent_rcpp)
{
	int M = mu_rcpp.size();
	int n = y_rcpp.size();
	arma::mat eta(n, M);

	arma::colvec y(y_rcpp.begin(), y_rcpp.size(), false);
	arma::mat    y_lagged(y_lagged_rcpp.begin(),
												y_lagged_rcpp.nrow(), y_lagged_rcpp.ncol(), false);
	arma::mat    z_dependent(z_dependent_rcpp.begin(),
													z_dependent_rcpp.nrow(), z_dependent_rcpp.ncol(), false);
	arma::mat    z_independent(z_independent_rcpp.begin(),
														z_independent_rcpp.nrow(), z_independent_rcpp.ncol(), false);
	arma::colvec beta(beta_rcpp.begin(), beta_rcpp.size(), false);
	arma::colvec mu(mu_rcpp.begin(), mu_rcpp.size(), false);
	arma::colvec sigma(sigma_rcpp.begin(), sigma_rcpp.size(), false);
	arma::mat    gamma_dependent(gamma_dependent_rcpp.begin(),
																gamma_dependent_rcpp.nrow(), gamma_dependent_rcpp.ncol(), false);
	arma::colvec gamma_independent(gamma_independent_rcpp.begin(),
																	gamma_independent_rcpp.size(), false);


	for (int j = 0; j < M; j++)
	{
  	eta.col(j) = y - y_lagged * beta -
								z_dependent * gamma_dependent.col(j) -
								z_independent * gamma_independent - mu(j);
		eta.col(j) = eta.col(j) % eta.col(j); // element-wise multiplication
		eta.col(j) = exp(-eta.col(j) / (2 * (sigma(j) * sigma(j))));
		eta.col(j) = eta.col(j) / (SQRT2PI * sigma(j));
	}

	return (wrap(eta));
}

// [[Rcpp::export]]
arma::mat MVGAUSS(Rcpp::NumericMatrix OMEGA_, int n, int seed) {

  //  std::srand(12523);

  arma::mat OMEGA( OMEGA_.begin(), OMEGA_.nrow(), OMEGA_.ncol(), false );

  arma::vec eigval;
  arma::mat eigvec;
  arma::eig_sym(eigval,eigvec, OMEGA);

 int ncol = OMEGA.n_cols;

  arma::mat X = arma::randn<arma::mat>(n,ncol);

  eigval = arma::sqrt(eigval);
  arma::mat Z = arma::diagmat(eigval);
  X = eigvec * Z * X.t();
  return(X.t());
}

RcppExport SEXP rthxpos(SEXP m) 
{
   Rcpp::NumericMatrix tmpm = Rcpp::NumericMatrix(m);
   int nr = tmpm.nrow();
   int nc = tmpm.ncol();
   thrust::device_vector<double> dmat(tmpm.begin(),tmpm.end());
   // make space for the transpose
   thrust::device_vector<double> dxp(nr*nc);
   // iterator to march through the matrix elements
   thrust::counting_iterator<int> seqb(0);
   thrust::counting_iterator<int> seqe = seqb + nr*nc;
   // for each i in seq, copy the matrix elt to its spot in the
   // transpose
   thrust::for_each(seqb,seqe,
      copyelt2xp(dmat.begin(),dxp.begin(),nr,nc));
   // prepare the R output, and return it
   Rcpp::NumericVector routmat(nc*nr);
   thrust::copy(dxp.begin(),dxp.end(),routmat.begin());
   return routmat;
}

    /** 
     * Update L, Ut and cu for new weights.
     *
     * Update Ut from Zt and sqrtXwt, then L from Lambda and Ut
     * Update cu from wtres, Lambda and Ut.
     * 
     * @param Xwt Matrix of weights for the model matrix
     * @param wtres weighted residuals
     */
    void reModule::reweight(Rcpp::NumericMatrix const&   Xwt,
			    Rcpp::NumericVector const& wtres) {
	double mone = -1., one = 1.; 
	int Wnc = Xwt.ncol();
	if (d_Ut) M_cholmod_free_sparse(&d_Ut, &c);
	if (Wnc == 1) {
	    d_Ut = M_cholmod_copy_sparse(&d_Zt, &c);
	    chmDn csqrtX(Xwt);
	    M_cholmod_scale(&csqrtX, CHOLMOD_COL, d_Ut, &c);
	} else {
	    int n = Xwt.nrow();
	    CHM_TR tr = M_cholmod_sparse_to_triplet(&d_Zt, &c);
	    int *j = (int*)tr->j, nnz = tr->nnz;
	    double *x = (double*)tr->x, *W = Xwt.begin();
	    for (int k = 0; k < nnz; k++) {
		x[k] *= W[j[k]];
		j[k] = j[k] % n;
	    }
	    tr->ncol = (size_t)n;

	    d_Ut = M_cholmod_triplet_to_sparse(tr, nnz, &c);
	    M_cholmod_free_triplet(&tr, &c);
	}
				// update the factor L
	CHM_SP LambdatUt = d_Lambda.crossprod(d_Ut);
	d_L.update(*LambdatUt, 1.);
	d_ldL2 = d_L.logDet2();
				// update cu
	chmDn ccu(d_cu), cwtres(wtres);
	copy(d_u0.begin(), d_u0.end(), d_cu.begin());
	M_cholmod_sdmult(LambdatUt, 0/*trans*/, &one, &mone, &cwtres, &ccu, &c);
	M_cholmod_free_sparse(&LambdatUt, &c);
	NumericMatrix
	    ans = d_L.solve(CHOLMOD_L, d_L.solve(CHOLMOD_P, d_cu));
	copy(ans.begin(), ans.end(), d_cu.begin());
	d_CcNumer = inner_product(d_cu.begin(), d_cu.end(), d_cu.begin(), double());
    }

///' Calculate the network properties, data matrix not provided
///' 
///' @details
///' \subsection{Input expectations:}{
///'   Note that this function expects all inputs to be sensible, as checked by
///'   the R function 'checkUserInput' and processed by 'networkProperties'. 
///'   
///'   These requirements are:
///'   \itemize{
///'   \item{'net' is a square matrix, and its rownames are identical to its 
///'         column names.}
///'   \item{'moduleAssigments' is a named character vector, where the names
///'         represent node labels found in the discovery dataset. Unlike 
///'         'PermutationProcedure', these may include nodes that are not 
///'         present in 'data' and 'net'.}
///'   \item{The module labels specified in 'modules' must occur in 
///'         'moduleAssignments'.}
///'   }
///' }
///' 
///' @param net adjacency matrix of network edge weights between all pairs of 
///'   nodes in the dataset in which to calculate the network properties.
///' @param moduleAssignments a named character vector containing the module 
///'   each node belongs to in the discovery dataset. 
///' @param modules a character vector of modules for which to calculate the 
///'   network properties for.
///' 
///' @return a list containing the summary profile, node contribution, module
///'   coherence, weighted degree, and average edge weight for each 'module'.
///'   
///' @keywords internal
// [[Rcpp::export]]
Rcpp::List NetPropsNoData (
    Rcpp::NumericMatrix net, 
    Rcpp::CharacterVector moduleAssignments,
    Rcpp::CharacterVector modules
) {
  // convert the colnames / rownames to C++ equivalents
  const std::vector<std::string> nodeNames (Rcpp::as<std::vector<std::string>>(colnames(net)));
  unsigned int nNodes = net.ncol();
  
  R_CheckUserInterrupt(); 
  
  /* Next, we need to create two mappings:
  *  - From node IDs to indices in the dataset of interest
  *  - From modules to node IDs
  *  - From modules to only node IDs present in the dataset of interest
  */
  const namemap nodeIdxMap = MakeIdxMap(nodeNames);
  const stringmap modNodeMap = MakeModMap(moduleAssignments);
  const stringmap modNodePresentMap = MakeModMap(moduleAssignments, nodeIdxMap);
  
  // What modules do we actually want to analyse?
  const std::vector<std::string> mods (Rcpp::as<std::vector<std::string>>(modules));
  
  R_CheckUserInterrupt(); 
  
  // Calculate the network properties for each module
  std::string mod; // iterators
  unsigned int mNodesPresent, mNodes;
  arma::uvec nodeIdx, propIdx, nodeRank;
  namemap propIdxMap;
  std::vector<std::string> modNodeNames; 
  arma::vec WD; // results containers
  double avgWeight; 
  Rcpp::NumericVector degree; // for casting to R equivalents
  Rcpp::List results; // final storage container
  for (auto mi = mods.begin(); mi != mods.end(); ++mi) {
    // What nodes are in this module?
    // modNodeNames = names(moduleAssignments[moduleAssignments == mod])
    mod = *mi;
    modNodeNames = GetModNodeNames(mod, modNodeMap);
    
    // initialise results containers with NA values for nodes not present in
    // the dataset we're calculating the network properties in.
    degree = Rcpp::NumericVector(modNodeNames.size(), NA_REAL);
    avgWeight = NA_REAL;
    degree.names() = modNodeNames;
    
    // Create a mapping between node names and the result vectors
    propIdxMap = MakeIdxMap(modNodeNames);
    
    // Get just the indices of nodes that are present in the requested dataset
    nodeIdx = GetNodeIdx(mod, modNodePresentMap, nodeIdxMap);
    mNodesPresent = nodeIdx.n_elem;
    
    // And a mapping of those nodes to the initialised vectors
    propIdx = GetNodeIdx(mod, modNodePresentMap, propIdxMap);
    mNodes = propIdx.n_elem;

    // Calculate the properties if the module has nodes in the test dataset
    if (nodeIdx.n_elem > 0) {
      // sort the node indices for sequential memory access
      nodeRank = SortNodes(nodeIdx.memptr(), mNodesPresent);
      
      WD = WeightedDegree(net.begin(), nNodes, nodeIdx.memptr(), mNodesPresent);
      WD = WD(nodeRank); // reorder results
      
      avgWeight = AverageEdgeWeight(WD.memptr(), WD.n_elem);
      R_CheckUserInterrupt(); 
//.........这里部分代码省略.........

// [[Rcpp::export]]
Rcpp::List rangerCpp(uint treetype, std::string dependent_variable_name,
    Rcpp::NumericMatrix input_data, std::vector<std::string> variable_names, uint mtry, uint num_trees, bool verbose,
    uint seed, uint num_threads, bool write_forest, uint importance_mode_r, uint min_node_size,
    std::vector<std::vector<double>>& split_select_weights, bool use_split_select_weights,
    std::vector<std::string>& always_split_variable_names, bool use_always_split_variable_names,
    std::string status_variable_name, bool prediction_mode, Rcpp::List loaded_forest, Rcpp::RawMatrix sparse_data,
    bool sample_with_replacement, bool probability, std::vector<std::string>& unordered_variable_names,
    bool use_unordered_variable_names, bool save_memory, uint splitrule_r, 
    std::vector<double>& case_weights, bool use_case_weights, bool predict_all, 
    bool keep_inbag, double sample_fraction, double alpha, double minprop, bool holdout, uint prediction_type_r) {

  Rcpp::List result;
  Forest* forest = 0;
  Data* data = 0;
  try {

    // Empty split select weights and always split variables if not used
    if (!use_split_select_weights) {
      split_select_weights.clear();
    }
    if (!use_always_split_variable_names) {
      always_split_variable_names.clear();
    }
    if (!use_unordered_variable_names) {
      unordered_variable_names.clear();
    }
    if (!use_case_weights) {
      case_weights.clear();
    }

    std::ostream* verbose_out;
    if (verbose) {
      verbose_out = &Rcpp::Rcout;
    } else {
      verbose_out = new std::stringstream;
    }

    size_t num_rows = input_data.nrow();
    size_t num_cols = input_data.ncol();

    // Initialize data with double memmode
    data = new DataDouble(input_data.begin(), variable_names, num_rows, num_cols);

    // If there is sparse data, add it
    if (sparse_data.nrow() > 1) {
      data->addSparseData(sparse_data.begin(), sparse_data.ncol());
    }

    switch (treetype) {
    case TREE_CLASSIFICATION:
      if (probability) {
        forest = new ForestProbability;
      } else {
        forest = new ForestClassification;
      }
      break;
    case TREE_REGRESSION:
      forest = new ForestRegression;
      break;
    case TREE_SURVIVAL:
      forest = new ForestSurvival;
      break;
    case TREE_PROBABILITY:
      forest = new ForestProbability;
      break;
    }

    ImportanceMode importance_mode = (ImportanceMode) importance_mode_r;
    SplitRule splitrule = (SplitRule) splitrule_r;
    PredictionType prediction_type = (PredictionType) prediction_type_r;

    // Init Ranger
    forest->initR(dependent_variable_name, data, mtry, num_trees, verbose_out, seed, num_threads,
        importance_mode, min_node_size, split_select_weights, always_split_variable_names, status_variable_name,
        prediction_mode, sample_with_replacement, unordered_variable_names, save_memory, splitrule, case_weights, 
        predict_all, keep_inbag, sample_fraction, alpha, minprop, holdout, prediction_type);

    // Load forest object if in prediction mode
    if (prediction_mode) {
      size_t dependent_varID = loaded_forest["dependent.varID"];
      //size_t num_trees = loaded_forest["num.trees"];
      std::vector<std::vector<std::vector<size_t>> > child_nodeIDs = loaded_forest["child.nodeIDs"];
      std::vector<std::vector<size_t>> split_varIDs = loaded_forest["split.varIDs"];
      std::vector<std::vector<double>> split_values = loaded_forest["split.values"];
      std::vector<bool> is_ordered = loaded_forest["is.ordered"];

      if (treetype == TREE_CLASSIFICATION) {
        std::vector<double> class_values = loaded_forest["class.values"];
        ((ForestClassification*) forest)->loadForest(dependent_varID, num_trees, child_nodeIDs, split_varIDs,
            split_values, class_values, is_ordered);
      } else if (treetype == TREE_REGRESSION) {
        ((ForestRegression*) forest)->loadForest(dependent_varID, num_trees, child_nodeIDs, split_varIDs, split_values,
            is_ordered);
      } else if (treetype == TREE_SURVIVAL) {
        size_t status_varID = loaded_forest["status.varID"];
        std::vector<std::vector<std::vector<double>> > chf = loaded_forest["chf"];
        std::vector<double> unique_timepoints = loaded_forest["unique.death.times"];
        ((ForestSurvival*) forest)->loadForest(dependent_varID, num_trees, child_nodeIDs, split_varIDs, split_values,
            status_varID, chf, unique_timepoints, is_ordered);
//.........这里部分代码省略.........

	IndepTestGauss(uint sampleSize, Rcpp::NumericMatrix& cor) :
		_sampleSize(sampleSize),
		_correlation(cor.begin(), cor.nrow(), cor.ncol(), false) {}

//' rcpp_get_polygons
//'
//' Extracts all polygons from an overpass API query
//'
//' @param st Text contents of an overpass API query
//' @return A \code{SpatialLinesDataFrame} contains all polygons and associated data
// [[Rcpp::export]]
Rcpp::S4 rcpp_get_polygons (const std::string& st)
{
#ifdef DUMP_INPUT
    {
        std::ofstream dump ("./get-polygons.xml");
        if (dump.is_open())
        {
            dump.write (st.c_str(), st.size());
        }
    }
#endif

    XmlPolys xml (st);

    const std::map <osmid_t, Node>& nodes = xml.nodes ();
    const std::map <osmid_t, OneWay>& ways = xml.ways ();
    const std::vector <Relation>& rels = xml.relations ();

    int count = 0;
    float xmin = FLOAT_MAX, xmax = -FLOAT_MAX,
          ymin = FLOAT_MAX, ymax = -FLOAT_MAX;
    std::vector <float> lons, lats;
    std::unordered_set <std::string> idset; // see TODO below
    std::vector <std::string> colnames, rownames, polynames;
    std::set <std::string> varnames;
    Rcpp::List dimnames (0), dummy_list (0);
    Rcpp::NumericMatrix nmat (Rcpp::Dimension (0, 0));

    idset.clear ();

    colnames.push_back ("lon");
    colnames.push_back ("lat");
    varnames.insert ("name");
    // other varnames added below

    /*
     * NOTE: Nodes are first loaded into the 2 vectors of (lon, lat), and these
     * are then copied into nmat. This intermediate can be avoided by loading
     * directly into nmat using direct indexing rather than iterators, however
     * this does *NOT* make the routine any faster, and so the current version
     * which more safely uses iterators is kept instead.
     */

    Rcpp::Environment sp_env = Rcpp::Environment::namespace_env ("sp");
    Rcpp::Function Polygon = sp_env ["Polygon"];
    Rcpp::Language polygons_call ("new", "Polygons");
    Rcpp::S4 polygons;

    /*
     * Polygons are extracted from the XmlPolys class in three setps:
     *  1. Get the names of all polygons that are part of multipolygon relations
     *  2. Get the names of any remaining ways that are polygonal (start == end)
     *  3. From the resultant list, extract the actual polygonal ways
     *
     * NOTE: OSM polygons are stored as ways, and thus all objects in the class
     * xmlPolys are rightly referred to as ways. Here within this Rcpp function,
     * these are referred to as Polygons, but the iteration is over the actual
     * polygonal ways.
     */

    // Step#1
    std::set <osmid_t> the_ways; // Set will only insert unique values
    for (auto it = rels.begin (); it != rels.end (); ++it)
        for (auto itw = (*it).ways.begin (); itw != (*it).ways.end (); ++itw)
        {
            assert (ways.find (itw->first) != ways.end ());
            the_ways.insert (itw->first);
        }

    // Step#2
    //const std::map <osmid_t, OneWay>& ways = xml.ways ();
    for (auto it = ways.begin (); it != ways.end (); ++it)
    {
        if (the_ways.find ((*it).first) == the_ways.end ())
            if ((*it).second.nodes.begin () == (*it).second.nodes.end ())
                the_ways.insert ((*it).first);
    }
    // Step#2b - Erase any ways that contain no data (should not happen).
    for (auto it = the_ways.begin (); it != the_ways.end (); )
    {
        auto itw = ways.find (*it);
        if (itw->second.nodes.size () == 0)
            it = the_ways.erase (it);
        else
            ++it;
    }
    Rcpp::List polyList (the_ways.size ());

    // Step#3 - Extract and store the_ways
    for (auto it = the_ways.begin (); it != the_ways.end (); ++it)
    {
        auto itw = ways.find (*it);
        // Collect all unique keys
//.........这里部分代码省略.........

///' Calculate the network properties 
///' 
///' @details
///' \subsection{Input expectations:}{
///'   Note that this function expects all inputs to be sensible, as checked by
///'   the R function 'checkUserInput' and processed by 'networkProperties'. 
///'   
///'   These requirements are:
///'   \itemize{
///'   \item{The ordering of node names across 'data' and 'net' is consistent.}
///'   \item{The columns of 'data' are the nodes.}
///'   \item{'net' is a square matrix, and its rownames are identical to its 
///'         column names.}
///'   \item{'moduleAssigments' is a named character vector, where the names
///'         represent node labels found in the discovery dataset. Unlike 
///'         'PermutationProcedure', these may include nodes that are not 
///'         present in 'data' and 'net'.}
///'   \item{The module labels specified in 'modules' must occur in 
///'         'moduleAssignments'.}
///'   }
///' }
///' 
///' @param data data matrix from the dataset in which to calculate the network
///'   properties.
///' @param net adjacency matrix of network edge weights between all pairs of 
///'   nodes in the dataset in which to calculate the network properties.
///' @param moduleAssignments a named character vector containing the module 
///'   each node belongs to in the discovery dataset. 
///' @param modules a character vector of modules for which to calculate the 
///'   network properties for.
///' 
///' @return a list containing the summary profile, node contribution, module
///'   coherence, weighted degree, and average edge weight for each 'module'.
///'   
///' @keywords internal
// [[Rcpp::export]]
Rcpp::List NetProps (
    Rcpp::NumericMatrix data, Rcpp::NumericMatrix net, 
    Rcpp::CharacterVector moduleAssignments,
    Rcpp::CharacterVector modules
) {
  // First, scale the matrix data
  unsigned int nSamples = data.nrow();
  unsigned int nNodes = data.ncol();
  arma::mat scaledData = Scale(data.begin(), nSamples, nNodes);
  
  R_CheckUserInterrupt(); 
  
  // convert the colnames / rownames to C++ equivalents
  const std::vector<std::string> nodeNames (Rcpp::as<std::vector<std::string>>(colnames(net)));
  const std::vector<std::string> sampleNames (Rcpp::as<std::vector<std::string>>(rownames(data)));
  
  /* Next, we need to create two mappings:
  *  - From node IDs to indices in the dataset of interest
  *  - From modules to node IDs
  *  - From modules to only node IDs present in the dataset of interest
  */
  const namemap nodeIdxMap = MakeIdxMap(nodeNames);
  const stringmap modNodeMap = MakeModMap(moduleAssignments);
  const stringmap modNodePresentMap = MakeModMap(moduleAssignments, nodeIdxMap);
  
  // What modules do we actually want to analyse?
  const std::vector<std::string> mods (Rcpp::as<std::vector<std::string>>(modules));
  
  R_CheckUserInterrupt(); 
  
  // Calculate the network properties for each module
  std::string mod; // iterators
  unsigned int mNodesPresent, mNodes;
  arma::uvec nodeIdx, propIdx, nodeRank;
  namemap propIdxMap;
  std::vector<std::string> modNodeNames; 
  arma::vec WD, SP, NC; // results containers
  double avgWeight, coherence; 
  Rcpp::NumericVector degree, summary, contribution; // for casting to R equivalents
  Rcpp::List results; // final storage container
  for (auto mi = mods.begin(); mi != mods.end(); ++mi) {
    // What nodes are in this module?
    mod = *mi;
    modNodeNames = GetModNodeNames(mod, modNodeMap);
    
    // initialise results containers with NA values for nodes not present in
    // the dataset we're calculating the network properties in.
    degree = Rcpp::NumericVector(modNodeNames.size(), NA_REAL);
    contribution = Rcpp::NumericVector(modNodeNames.size(), NA_REAL);
    summary = Rcpp::NumericVector(nSamples, NA_REAL);
    avgWeight = NA_REAL;
    coherence = NA_REAL;
    degree.names() = modNodeNames;
    contribution.names() = modNodeNames;
    
    // Create a mapping between node names and the result vectors
    propIdxMap = MakeIdxMap(modNodeNames);
    
    // Get just the indices of nodes that are present in the requested dataset
    nodeIdx = GetNodeIdx(mod, modNodePresentMap, nodeIdxMap);
    mNodesPresent = nodeIdx.n_elem;
    
    // And a mapping of those nodes to the initialised vectors
    propIdx = GetNodeIdx(mod, modNodePresentMap, propIdxMap);
//.........这里部分代码省略.........