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

//' Evaluate the log-density of the measurement process by calling measurement
//' process density functions via external Xptr.
//'
//' @param emitmat matrix of emission probabilities
//' @param obsmat matrix containing the data
//' @param statemat matrix containing the compartment counts at the observation
//'   times
//' @param measproc_indmat logical matrix indicating which compartments are
//'   observed at every observation time
//' @param parameters numeric vector of parameter values
//' @param constants numeric vector of constants
//' @param tcovar_censusmat numeric vector of time-varying covariate values
//' @param d_meas_ptr external pointer to measurement process density function
//'
//' @export
// [[Rcpp::export]]
void evaluate_d_measure(Rcpp::NumericMatrix& emitmat, const Rcpp::NumericMatrix& obsmat,
                        const Rcpp::NumericMatrix& statemat, const Rcpp::LogicalMatrix& measproc_indmat,
                        const Rcpp::NumericVector& parameters, const Rcpp::NumericVector& constants,
                        const Rcpp::NumericMatrix& tcovar_censusmat, SEXP d_meas_ptr) {

        // get dimensions
        Rcpp::IntegerVector emit_dims = emitmat.attr("dim");

        // evaluate the densities
        for(int j=0; j < emit_dims[0]; ++j) {
                // args: emitmat, emit_inds, record_ind, record, state, parameters, constants, tcovar, pointer
                CALL_D_MEASURE(emitmat, measproc_indmat.row(j), j, obsmat.row(j), statemat.row(j), parameters, constants, tcovar_censusmat.row(j), d_meas_ptr);
        }
}

//' @title Calculating validation scores between two adjacency matrices
//' 
//' @description
//' This function calculates the validation scores between two adjacency matrices.
//' 
//' @param inf_mat matrix. It should be adjacency matrix of inferred network.
//' @param true_mat matrix. It should be adjacency matrix of true network.
// [[Rcpp::export]]
Rcpp::NumericVector rcpp_validate(Rcpp::NumericMatrix inf_mat, Rcpp::NumericMatrix true_mat) {
  if(inf_mat.ncol() != true_mat.ncol()) {
    throw std::invalid_argument( "Two input matrices should have the same number of columns." );
  }
  if(inf_mat.nrow() != true_mat.nrow()) {
    throw std::invalid_argument( "Two input matrices should have the same number of rows." );
  }
  
  int tp=0;
  int tn=0;
  int fp=0;
  int fn=0;
  for(signed int i=0; i<inf_mat.nrow(); i++) {
    //Convert R objects into C++ objects.
    Rcpp::NumericVector xr = inf_mat.row(i);
    Rcpp::NumericVector yr = true_mat.row(i);
    std::vector<int> x = Rcpp::as<std::vector <int> >(xr);
    std::vector<int> y = Rcpp::as<std::vector <int> >(yr);
    std::vector<int> z;
    
    //Calculate the frequency of numbers.
    //tp=true positive [1,1], tn=true negative [0,0], fp=false positive [1,0], fn=false negative [0,1].
    for(unsigned int k=0; k<x.size(); k++) {
      z.push_back(x[k] + y[k]); //Calculate the summation of x and y between each element.
                
      if(z[k] == 2) {
        tp += 1;
      } else if(z[k] == 0) {
        tn += 1;
      } else if(z[k] == 1) {
        if(x[k] == 0) {
          fp += 1;
        } else {
          fn += 1;
        }
      } else {
        throw std::invalid_argument("Error in calculating the contigency table.");
      }
    }
  }
  
  //std::vector<int> output{tp, tn, fp, fn}; c++11 only
  int tmp_arr[4] = {tp, tn, fp, fn};
  std::vector<int> output(&tmp_arr[0], &tmp_arr[0]+4);

  return Rcpp::wrap(output);
}

// gradient-descent ------------------------------------------------------------
//' @title Regular gradient descent for latent factor model with covariates
//' @param X The ratings matrix. Unobserved entries must be marked NA. Users
//' must be along rows, and tracks must be along columns.
//' @param Z The covariates associated with each pair. This must be shaped as an
//' array with users along rows, tracks along columns, and features along
//' slices.
//' @param k_factors The number of latent factors for the problem.
//' @param lambdas The regularization parameters for P, Q, and beta, respectively.
//' @param n_iter The number of gradient descent iterations to run.
//' @param batch_samples The proportion of the observed entries to sample in the
//' gradient for each factor in the gradient descent.
//' @param batch_factors The proportion of latent factors to update at each
//' iteration.
//' @param gamma The step-size in the gradient descent.
//' @param verbose Print the objective at each iteration of the descent?
//' @return A list with the following elements, \cr
//'   $P The learned user latent factors. \cr
//'   $Q The learned track latent factors. \cr
//'   $beta The learned regression coefficients.
//' @export
// [[Rcpp::export]]
Rcpp::List svd_cov(Rcpp::NumericMatrix X, Rcpp::NumericVector Z_vec,
		   int k_factors, Rcpp::NumericVector lambdas, int n_iter,
		   double batch_samples, double batch_factors,
		   double gamma_pq, double gamma_beta, bool verbose) {

  // convert to arma
  Rcpp::IntegerVector Z_dim = Z_vec.attr("dim");
  arma::cube Z(Z_vec.begin(), Z_dim[0], Z_dim[1], Z_dim[2], false);

  // initialize results
  int m = X.nrow();
  int n = X.ncol();
  int l = Z_dim[2];

  arma::mat P = 1.0 / sqrt(m) * arma::randn(m, k_factors);
  arma::mat Q = 1.0 / sqrt(n) * arma::randn(n, k_factors);
  arma::vec beta = 1.0 / sqrt(l) * arma::randn(l);
  arma::vec objs = arma::zeros(n_iter);

  // perform gradient descent steps
  for(int cur_iter = 0; cur_iter < n_iter; cur_iter++) {
    // update user factors
    arma::uvec cur_factors = get_batch_ix(m, batch_factors);
    for(int i = 0; i < cur_factors.n_elem; i++) {
      int cur_ix = cur_factors(i);
      P.row(cur_ix) = update_factor(X.row(cur_ix),
				    Z(arma::span(cur_ix), arma::span(), arma::span()),
				    arma::conv_to<arma::vec>::from(P.row(cur_ix)), Q, beta,
				    batch_samples, lambdas[0], gamma_pq).t();
    }

    // update track factors
    cur_factors = get_batch_ix(n, batch_factors);
    for(int j = 0; j < cur_factors.n_elem; j++) {
      int cur_ix = cur_factors(j);
      Q.row(cur_ix) = update_factor(X.column(cur_ix),
				    Z(arma::span(), arma::span(cur_ix), arma::span()),
				    arma::conv_to<arma::vec>::from(Q.row(cur_ix)), P, beta,
				    batch_samples, lambdas[1], gamma_pq).t();
    }

    // update regression coefficients
    beta = update_beta(Rcpp::as<arma::mat>(X), Z, P, Q, beta, lambdas[2], gamma_beta);
    objs[cur_iter] = objective_fun(Rcpp::as<arma::mat>(X), Z, P, Q, beta, lambdas);
    if(verbose) {
      printf("iteration %d \n", cur_iter);
      printf("Objective: %g \n", objs[cur_iter]);
    }
  }

  return Rcpp::List::create(Rcpp::Named("P") = P,
			    Rcpp::Named("Q") = Q,
			    Rcpp::Named("beta") = beta,
			    Rcpp::Named("objs") = objs);
}

//' Simulate an LNA path using a non-centered parameterization for the
//' log-transformed counting process LNA.
//'
//' @param lna_times vector of interval endpoint times
//' @param lna_draws vector of N(0,1) draws to be mapped to the path
//' @param lna_pars numeric matrix of parameters, constants, and time-varying
//'   covariates at each of the lna_times
//' @param init_start index in the parameter vector where the initial compartment
//'   volumes start
//' @param param_update_inds logical vector indicating at which of the times the
//'   LNA parameters need to be updated.
//' @param stoich_matrix stoichiometry matrix giving the changes to compartments
//'   from each reaction
//' @param forcing_inds logical vector of indicating at which times in the
//'   time-varying covariance matrix a forcing is applied.
//' @param forcing_matrix matrix containing the forcings.
//' @param max_attempts maximum number of tries if the first increment is rejected
//' @param step_size initial step size for the ODE solver (adapted internally,
//' but too large of an initial step can lead to failure in stiff systems).
//' @param lna_pointer external pointer to the compiled LNA integration function.
//' @param set_pars_pointer external pointer to the function for setting LNA pars.
//' @return list containing the stochastic perturbations (i.i.d. N(0,1) draws) and
//' the LNA path on its natural scale which is determined by the perturbations.
//'
//' @export
// [[Rcpp::export]]
Rcpp::List propose_lna(const arma::rowvec& lna_times,
                       const Rcpp::NumericVector& lna_draws,
                       const Rcpp::NumericMatrix& lna_pars,
                       const Rcpp::IntegerVector& lna_param_inds,
                       const Rcpp::IntegerVector& lna_tcovar_inds,
                       const int init_start,
                       const Rcpp::LogicalVector& param_update_inds,
                       const arma::mat& stoich_matrix,
                       const Rcpp::LogicalVector& forcing_inds,
                       const arma::uvec& forcing_tcov_inds,
                       const arma::mat& forcings_out,
                       const arma::cube& forcing_transfers,
                       int max_attempts,
                       double step_size,
                       SEXP lna_pointer,
                       SEXP set_pars_pointer) {
      
        // get the dimensions of various objects
        int n_events = stoich_matrix.n_cols;         // number of transition events, e.g., S2I, I2R
        int n_comps  = stoich_matrix.n_rows;         // number of model compartments (all strata)
        int n_odes   = n_events + n_events*n_events; // number of ODEs
        int n_times  = lna_times.n_elem;             // number of times at which the LNA must be evaluated
        int n_tcovar = lna_tcovar_inds.size();       // number of time-varying covariates or parameters
        int n_forcings = forcing_tcov_inds.n_elem;   // number of forcings
        
        // for use with forcings
        double forcing_flow = 0;
        arma::vec forcing_distvec(n_comps, arma::fill::zeros);

        // initialize the objects used in each time interval
        double t_L = 0;
        double t_R = 0;
        Rcpp::NumericVector current_params = lna_pars.row(0);   // vector for storing the current parameter values
        CALL_SET_ODE_PARAMS(current_params, set_pars_pointer);  // set the parameters in the odeintr namespace

        // initial state vector - copy elements from the current parameter vector
        arma::vec init_volumes(current_params.begin() + init_start, n_comps);

        // initialize the LNA objects - the vector for storing the current state
        Rcpp::NumericVector lna_state_vec(n_odes);   // vector to store the results of the ODEs
        
        arma::vec lna_drift(n_events, arma::fill::zeros); // incidence mean vector (natural scale)
        arma::mat lna_diffusion(n_events, n_events, arma::fill::zeros); // diffusion matrix
        
        arma::vec log_lna(n_events, arma::fill::zeros);  // LNA increment, log scale
        arma::vec nat_lna(n_events, arma::fill::zeros);  // LNA increment, natural scale
        
        // Objects for computing the eigen decomposition of the LNA covariance matrices
        arma::vec svd_d(n_events, arma::fill::zeros);
        arma::mat svd_U(n_events, n_events, arma::fill::zeros);
        arma::mat svd_V(n_events, n_events, arma::fill::zeros);
        bool good_svd = true;
        
        // matrix in which to store the LNA path
        arma::mat lna_path(n_events+1, n_times, arma::fill::zeros); // incidence path
        arma::mat prev_path(n_comps+1, n_times, arma::fill::zeros); // prevalence path (compartment volumes)
        lna_path.row(0) = lna_times;
        prev_path.row(0) = lna_times;
        prev_path(arma::span(1,n_comps), 0) = init_volumes;
        
        // apply forcings if called for - applied after censusing at the first time
        if(forcing_inds[0]) {
              
              // distribute the forcings proportionally to the compartment counts in the applicable states
              for(int j=0; j < n_forcings; ++j) {
                    
                    forcing_flow       = lna_pars(0, forcing_tcov_inds[j]);
                    forcing_distvec    = forcing_flow * normalise(forcings_out.col(j) % init_volumes, 1);
                    init_volumes      += forcing_transfers.slice(j) * forcing_distvec;
              }
        }
        
        // sample the stochastic perturbations - use Rcpp RNG for safety
        Rcpp::NumericVector draws_rcpp(Rcpp::clone(lna_draws));
//.........这里部分代码省略.........