C++ mat::elem方法代码示例

//' @title Create a taper setup for an Eigen sparse matrix
//' @description Generate the storage information we'll need to create a sparse eigen matrix
//' @param d A \code{matrix} that is a symmetric distance matrix.
//' @param taprange A \code{double} that is used to specify the gamma.
//' @return A \code{list}
//' \itemize{
//' \item \code{n} - row/col of matrix
//' \item \code{good.dists} - distances that meet tapering requirements
//' \item \code{taps} - results of applying the taper function
//' \item \code{ja} - column index
//' \item \code{ia} - row index
//' }
//' @details This function is meant to be used to create SPAM matrices.
//' The function formats data so that it can be used with a compressed sparse column format.
//' Do not use output for Eigen matrices.
//' @examples
//' data(anom1962)
//' d = rdist_earth1(head(loc))
//' make_tapersetup_eigen(d, 20)
// [[Rcpp::export]]
Rcpp::List make_tapersetup_eigen(arma::mat d, double taprange){
  double n = d.n_rows;
  
  // Find elements that are less than taprange and return element ids.
  arma::uvec inrange_id = find(d < taprange);
  
  // Approximate how many zero entries exist per columnss
  int rescol = floor(double(inrange_id.n_elem) / n); 

  // Obtain the small distance elements
  arma::vec good_dists = d.elem(inrange_id);
  
  // Apply taper function
  arma::vec taps = wendland2_1(good_dists, taprange);
  
  // Create a row index matrix and extract out ranges
  arma::vec ja = row_arma(d).elem(inrange_id) - 1;
  
  // Create a col index matrix and extract the ranges
  arma::vec ia = col_arma(d).elem(inrange_id) - 1;
    
  // Export results as an R list object
  return Rcpp::List::create(
                          Rcpp::Named("n") = n,
                          Rcpp::Named("good.dists") = good_dists,
                          Rcpp::Named("taps") = taps,
                          Rcpp::Named("ja") = ja,
                          Rcpp::Named("ia") = ia,
                          Rcpp::Named("rescol") = rescol);
}

//[[Rcpp::export]]
Rcpp::List changeNA(arma::mat X) {
  NumericVector x(X.begin(), X.end());
  int n = x.size();
  IntegerVector y = seq_len(n)-1;//(x.begin(), x.end());
  IntegerVector z = y[is_na(x)];
  if(z.size()==0){
    std::cout << "Nenhum valor faltante.\n";
    SEXP xNULL = R_NilValue;
    return Rcpp::List::create(xNULL);  
  } else{
    arma::uvec Xna = as<arma::uvec>(z);
    NumericVector bx(Xna.n_rows, 100.0);
    arma::vec b = as<arma::vec>(bx);
    X.elem(Xna) = b;
    return Rcpp::List::create(X);
  }  
}

//' @title Create a taper setup for a SPAM matrix
//' @description Get the storage information we'll need to create a sparse matrix
//' @param d A \code{matrix} that is a symmetric distance matrix.
//' @param taprange A \code{double} that is used to specify the gamma.
//' @return A \code{list}
//' \itemize{
//' \item \code{n} - row/col of matrix
//' \item \code{good.dists} - distances that meet tapering requirements
//' \item \code{taps} - results of applying the taper function
//' \item \code{ja} - ordered column indices
//' \item \code{ia} -  pointer to the beginning of each row in the arrays entries and colindices
//' \item \code{index} - element position in matrix 
//' }
//' @details This function is meant to be used to create SPAM matrices.
//' Do not use output for Eigen matrices.
//' @examples
//' data(anom1962)
//' d = rdist_earth1(head(loc))
//' make_tapersetup_R(d, 20)
// [[Rcpp::export]]
Rcpp::List make_tapersetup_R(arma::mat d, double taprange){
  unsigned int n = d.n_rows;
  
  // Find elements that are less than (logical matrix e.g. 0/1's)
  arma::umat inrange = d < taprange;
  
  // Same idea, but return element ids.
  arma::uvec inrange_id = find(inrange);
  
  // Convert matrix to allow operations
  arma::mat temp = arma::conv_to<arma::mat>::from(inrange);

  // Obtain the small distance elements
  arma::vec good_dists = d.elem(inrange_id);
  
  // Apply taper function (this cannot be made dynamic)
  arma::vec taps = wendland2_1(good_dists, taprange);
  
  // Create ordered column indices of the nonzero values
  arma::vec ja = row_arma(d).elem(inrange_id);
  
  // Calculate pointer to the beginning of each row in the arrays' entries and colindices
  arma::vec ia(inrange.n_rows + 1);
  ia(0) = 1;
  ia.rows(1,inrange.n_rows) = 1 + arma::cumsum(arma::sum(temp,1));
  
  // Create a column index matrix and isolate ideal elements
  arma::vec index = (col_arma(d).elem(inrange_id) - 1) * n + ja;
  
  // Export results as an R list object
  return Rcpp::List::create(
                          Rcpp::Named("n") = n,
                          Rcpp::Named("good_dists") = good_dists,
                          Rcpp::Named("taps") = taps,
                          Rcpp::Named("ja") = ja,
                          Rcpp::Named("ia") = ia,
                          Rcpp::Named("index") = index);
}

void getEystar_endorseIRT (const arma::mat &alpha,
                           const arma::mat &beta,
                           const arma::mat &w,
                           const arma::mat &theta,
                           const arma::mat &gamma,
                           const arma::mat &y,
                           const int N,
                           const int J,
                           arma::mat &ystars,
                           const arma::mat &theta2,
                           const arma::mat &w2
                           ) {
    // arma::mat ystars(N, J, arma::fill::zeros) ;

    // Main Calculation

#pragma omp parallel for
    for (int n = 0; n < N; n++) {

        arma::mat theseystars = arma::mat(1, J, arma::fill::zeros) ;

        for (int j = 0; j < J; j++) {
            // defaults to untruncated
            double low = y(n,j) == 1 ? 0.0 : R_NegInf ;
            double high = y(n,j) == 0 ? 0.0 : R_PosInf ;

            // mu

            double mu =
                alpha(j, 0) +
                beta(n, 0) -
                gamma(0, 0) * (
                               pow(theta(n, 0), 2) +
                               pow(w(j, 0), 2) -
                               2 * theta(n, 0) * w(j, 0)
                               ) ;

            double exp = etn1(mu, 1.0, low, high) ;

            double exp2 = exp ;

            // if (isinf(exp2) || isnan(exp2)) {
            //     if (low == 0.0) exp2 =   .0001 ;
            //     if (high == 0.0) exp2 = -.0001 ;
            // }


            // if (isnan(exp2)) {
            //     exp2 = 0 ;
            // }

            // if (isinf(exp)) {
            //     Rcout << n << " " << j << " " << low << " " << high << " " << mu << " " << exp << " " << exp2 << std::endl ;
            // }

            // if (isnan(exp)) {
            //     Rcout << n << " " << j << " " << low << " " << high << " " << mu << " " << exp << " " << exp2 << std::endl ;
            // }

            theseystars(0, j) = exp2 ;
        }
        ystars.row(n) = theseystars ;
    }

    arma::uvec check = find_nonfinite(ystars) ;

    if (check.n_elem > 0) {
        ystars.elem(check).print("ystar check") ;
    }

    // return(ystars) ;
}

//' Map N(0,1) stochastic perturbations to an LNA path.
//'
//' @param pathmat matrix where the LNA path should be stored
//' @param draws matrix of N(0,1) draws to be mapped to an LNA path
//' @param lna_times vector of interval endpoint times
//' @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 svd_d vector in which to store SVD singular values
//' @param svd_U matrix in which to store the U matrix of the SVD
//' @param svd_V matrix in which to store the V matrix of the SVD
//' @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 LNA integration function.
//' @param set_pars_pointer external pointer to the function for setting the LNA
//'   parameters.
//'
//' @return fill out pathmat with the LNA path corresponding to the stochastic
//'   perturbations.
//'
//' @export
// [[Rcpp::export]]
void map_draws_2_lna(arma::mat& pathmat,
                     const arma::mat& draws,
                     const arma::rowvec& lna_times,
                     const Rcpp::NumericMatrix& lna_pars,
                     Rcpp::NumericVector& lna_param_vec,
                     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,
                     arma::vec& svd_d,
                     arma::mat& svd_U,
                     arma::mat& svd_V,
                     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;

        // vector of parameters, initial compartment columes, constants, and time-varying covariates
        std::copy(lna_pars.row(0).begin(), lna_pars.row(0).end(), lna_param_vec.begin());

        CALL_SET_ODE_PARAMS(lna_param_vec, set_pars_pointer); // set the parameters in the odeintr namespace

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

        // initialize the LNA objects
        bool good_svd = true;

        Rcpp::NumericVector lna_state_vec(n_odes); // the vector for storing the current state of the LNA ODEs

        arma::vec lna_drift(n_events, arma::fill::zeros);               // incidence mean vector (log 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

        // 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;
              }
        }

        // iterate over the time sequence, solving the LNA over each interval
        for(int j=0; j < (n_times-1); ++j) {
//.........这里部分代码省略.........