C++ NumericVector::attr方法代码示例

// Bellman recursion using row rearrangement
//[[Rcpp::export]]
Rcpp::List Bellman(const arma::mat& grid,
                   Rcpp::NumericVector reward_,
                   const arma::cube& scrap,
                   Rcpp::NumericVector control_,
                   const arma::cube& disturb,
                   const arma::vec& weight) {
  // Passing R objects to C++
  const std::size_t n_grid = grid.n_rows;
  const std::size_t n_dim = grid.n_cols;
  const arma::ivec r_dims = reward_.attr("dim");
  const std::size_t n_pos = r_dims(3);
  const std::size_t n_action = r_dims(2);
  const std::size_t n_dec = r_dims(4) + 1;
  const arma::cube
      reward(reward_.begin(), n_grid, n_dim * n_action * n_pos, n_dec - 1, false);
  const arma::ivec c_dims = control_.attr("dim");
  arma::cube control2;
  arma::imat control;
  bool full_control;
  if (c_dims.n_elem == 3) {
    full_control = false;
    arma::cube temp_control2(control_.begin(), n_pos, n_action, n_pos, false);
    control2 = temp_control2;
  } else {
    full_control = true;
    arma::mat temp_control(control_.begin(), n_pos, n_action, false);
    control = arma::conv_to<arma::imat>::from(temp_control);
  }
  const std::size_t n_disturb = disturb.n_slices;
  // Bellman recursion
  arma::cube value(n_grid, n_dim * n_pos, n_dec);
  arma::cube cont(n_grid, n_dim * n_pos, n_dec - 1, arma::fill::zeros);
  arma::mat d_value(n_grid, n_dim);
  Rcpp::Rcout << "At dec: " << n_dec - 1 << "...";
  for (std::size_t pp = 0; pp < n_pos; pp++) {
    value.slice(n_dec - 1).cols(n_dim * pp, n_dim * (pp + 1) - 1) =
        scrap.slice(pp);
  }
  for (int tt = (n_dec - 2); tt >= 0; tt--) {
    Rcpp::Rcout << tt;
    // Approximating the continuation value
    for (std::size_t pp = 0; pp < n_pos; pp++) {
      cont.slice(tt).cols(n_dim * pp, n_dim * (pp + 1) - 1) =
          Expected(grid,
                   value.slice(tt + 1).cols(pp * n_dim, n_dim * (pp + 1) - 1),
                   disturb, weight);
    }
    Rcpp::Rcout << "..";
    // Optimise value function
    if (full_control) {
      BellmanOptimal(grid, control, value, reward, cont, tt);
    } else {
      BellmanOptimal2(grid, control2, value, reward, cont, tt);
    }
    Rcpp::Rcout << ".";
  }
  return Rcpp::List::create(Rcpp::Named("value") = value,
                            Rcpp::Named("expected") = cont);
}

//' Update eigen values, vectors, and inverse matrices for irms
//'
//' @param tpm_prods array of transition probability matrix products
//' @param tpms array of transition probability matrices
//' @param pop_mat population level bookkeeping matrix
//' @param obstimes vector of observation times
//'
//' @return Updated eigenvalues, eigenvectors, and inverse matrices
// [[Rcpp::export]]
void tpmProdSeqs(Rcpp::NumericVector& tpm_prods, Rcpp::NumericVector& tpms, const Rcpp::IntegerVector obs_time_inds) {

        // Get indices
        int n_obs = obs_time_inds.size();
        Rcpp::IntegerVector tpmDims = tpms.attr("dim");

        // Ensure obs_time_inds starts at 0
        Rcpp::IntegerVector obs_inds = obs_time_inds - 1;

        // Set array pointers
        arma::cube prod_arr(tpm_prods.begin(), tpmDims[0], tpmDims[1], tpmDims[2], false);
        arma::cube tpm_arr(tpms.begin(), tpmDims[0], tpmDims[1], tpmDims[2], false);

        // Generate tpm products and store them in the appropriate spots in the tpm product array
        for(int j = 0; j < (n_obs - 1); ++j) {

                prod_arr.slice(obs_inds[j+1] - 1) = tpm_arr.slice(obs_inds[j+1] - 1);

                for(int k = (obs_inds[j+1] - 2); k > (obs_inds[j] - 1); --k) {

                        prod_arr.slice(k) = tpm_arr.slice(k) * prod_arr.slice(k + 1);
                }

        }
}

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

//[[Rcpp::export]]
Rcpp::List map_data_set(Rcpp::NumericMatrix data_, Rcpp::NumericVector inpar,bool lc_) {

  svec parnames = Rcpp::as<svec>(inpar.attr("names"));

  dataobject* data = new dataobject(data_,parnames);

  data->map_uid();
  data->locate_tran(lc_);
  data -> idata_row();
  Rcpp::List ans = data->ex_port();

  delete data;

  return(ans);
}

//[[Rcpp::export]]
arma::cube PathDisturb(const arma::vec& start,
                       Rcpp::NumericVector disturb_) {
  // R objects to C++
  const arma::ivec d_dims = disturb_.attr("dim");
  const std::size_t n_dim = d_dims(0);
  const std::size_t n_dec = d_dims(3) + 1;
  const std::size_t n_path = d_dims(2);
  const arma::cube disturb(disturb_.begin(), n_dim, n_dim * n_path, n_dec - 1,
			  false);
  // Simulating the sample paths
  arma::cube path(n_path, n_dim, n_dec);
  // Assign starting values
  for (std::size_t ii = 0; ii < n_dim; ii++) {
    path.slice(0).col(ii).fill(start(ii));
  }
  // Disturb the paths
  for (std::size_t pp = 0; pp < n_path; pp++) {
    for (std::size_t tt = 1; tt < n_dec; tt++) {
      path.slice(tt).row(pp) = path.slice(tt - 1).row(pp) *
          arma::trans(disturb.slice(tt - 1).cols(n_dim * pp, n_dim * (pp + 1) - 1));
    }
  }
  return path;
}

// Calculate the martingale increments using the row rearrangement
//[[Rcpp::export]]
arma::cube AddDual(const arma::cube& path,
                   Rcpp::NumericVector subsim_,
                   const arma::vec& weight,
                   Rcpp::NumericVector value_,
                   Rcpp::Function Scrap_) {
  // R objects to C++
  const std::size_t n_path = path.n_rows;
  const arma::ivec v_dims = value_.attr("dim");
  const std::size_t n_grid = v_dims(0);
  const std::size_t n_dim = v_dims(1);
  const std::size_t n_pos = v_dims(2);
  const std::size_t n_dec = v_dims(3);
  const arma::cube value(value_.begin(), n_grid, n_dim * n_pos, n_dec, false);
  const arma::ivec s_dims = subsim_.attr("dim");
  const std::size_t n_subsim = s_dims(2);
  const arma::cube subsim(subsim_.begin(), n_dim, n_dim * n_subsim * n_path, n_dec - 1, false);
  // Duals
  arma::cube mart(n_path, n_pos, n_dec - 1);
  arma::mat temp_state(n_subsim * n_path, n_dim);
  arma::mat fitted(n_grid, n_dim);
  std::size_t ll;
  Rcpp::Rcout << "Additive duals at dec: ";
  // Find averaged value
  for (std::size_t tt = 0; tt < (n_dec - 2); tt++) {
    Rcpp::Rcout << tt << "...";
    // 1 step subsimulation
#pragma omp parallel for private(ll)
    for (std::size_t ii = 0; ii < n_path; ii++) {
      for (std::size_t ss = 0; ss < n_subsim; ss++) {
        ll = n_subsim * ii + ss;
        temp_state.row(ll) = weight(ss) * path.slice(tt).row(ii) *
            arma::trans(subsim.slice(tt).cols(n_dim * ll, n_dim * (ll + 1) - 1));       
      }
    }
    // Averaging
    for (std::size_t pp = 0; pp < n_pos; pp++) {
      fitted = value.slice(tt + 1).cols(n_dim * pp, n_dim * (pp + 1) - 1);
      mart.slice(tt).col(pp) = arma::conv_to<arma::vec>::from(arma::sum(arma::reshape(
          OptimalValue(temp_state, fitted), n_subsim, n_path)));
      // Subtract the path realisation
      mart.slice(tt).col(pp) -= OptimalValue(path.slice(tt + 1), fitted);
    }
  }
  // Scrap value
  Rcpp::Rcout << n_dec - 1 << "...";
  // 1 step subsimulation
#pragma omp parallel for private(ll)
  for (std::size_t ii = 0; ii < n_path; ii++) {
    for (std::size_t ss = 0; ss < n_subsim; ss++) {
      ll = n_subsim * ii + ss;
      temp_state.row(n_path * ss + ii) = path.slice(n_dec - 2).row(ii) *
          arma::trans(subsim.slice(n_dec - 2).cols(n_dim * ll, n_dim * (ll + 1) - 1));
    }
  }
  // Averaging
  arma::mat subsim_scrap(n_subsim * n_path, n_pos);
  subsim_scrap = Rcpp::as<arma::mat>(Scrap_(
      Rcpp::as<Rcpp::NumericMatrix>(Rcpp::wrap(temp_state))));
  arma::mat scrap(n_path, n_pos);
  scrap = Rcpp::as<arma::mat>(Scrap_(
      Rcpp::as<Rcpp::NumericMatrix>(Rcpp::wrap(path.slice(n_dec - 1)))));
  for (std::size_t pp = 0; pp < n_pos; pp++) {
    mart.slice(n_dec - 2).col(pp) =
        arma::reshape(subsim_scrap.col(pp), n_path, n_subsim) * weight;
    // Subtract the path realisation
    mart.slice(n_dec - 2).col(pp) -= scrap.col(pp);
  }
  Rcpp::Rcout << "done\n";
  return mart;
}

//[[Rcpp::export]]
Rcpp::List AddDualBounds(const arma::cube& path,
                         Rcpp::NumericVector control_,
                         Rcpp::Function Reward_,
                         Rcpp::Function Scrap_,
                         const arma::cube& mart,
                         const arma::ucube& path_action) {
  // Extract parameters
  const std::size_t n_dec = path.n_slices;
  const std::size_t n_path = path.n_rows;
  const std::size_t n_dim = path.n_cols;
  const arma::ivec c_dims = control_.attr("dim");
  const std::size_t n_pos = c_dims(0);
  const std::size_t n_action = c_dims(1);
  arma::imat control;  // full control
  arma::cube control2;  // partial control
  bool full_control;
  if (c_dims.n_elem == 3) {
    full_control = false;
    arma::cube temp_control2(control_.begin(), n_pos, n_action, n_pos, false);
    control2 = temp_control2;
  } else {
    full_control = true;
    arma::mat temp_control(control_.begin(), n_pos, n_action, false);
    control = arma::conv_to<arma::imat>::from(temp_control);
  }
  // Initialise with scrap value
  arma::cube primals(n_path, n_pos, n_dec);
  primals.slice(n_dec - 1) = Rcpp::as<arma::mat>(
      Scrap_(Rcpp::as<Rcpp::NumericMatrix>(Rcpp::wrap(path.slice(n_dec - 1)))));
  arma::cube duals = primals;
  // Perform the backward induction
  arma::uword policy;
  arma::cube reward(n_path, n_action, n_pos);
  if (full_control) {  // For the full control case
    arma::uword next;
    for (int tt = (n_dec - 2); tt >= 0; tt--) {
      reward = Rcpp::as<arma::cube>(Reward_(
          Rcpp::as<Rcpp::NumericMatrix>(Rcpp::wrap(path.slice(tt))), tt + 1));
#pragma omp parallel for private(policy, next)
      for (std::size_t ii = 0; ii < n_path; ii++) {
        for (std::size_t pp = 0; pp < n_pos; pp++) {
          // Primal values
          policy = path_action(ii, pp, tt) - 1;  // R to C indexing
          next = control(pp, policy) - 1;
          primals(ii, pp, tt) = reward(ii, policy, pp) + mart(ii, next, tt)
              + primals(ii, next, tt + 1);
          // Dual values
          next = control(pp, 0) - 1;
          duals(ii, pp, tt) = reward(ii, 0, pp) + mart(ii, next, tt)
              + duals(ii, next, tt + 1);
          for (std::size_t aa = 1; aa < n_action; aa++) {
            next = control(pp, aa) - 1;
            duals(ii, pp, tt) = std::max(reward(ii, aa, pp) + mart(ii, next, tt)
                         + duals(ii, next, tt + 1), duals(ii, pp, tt));
          }
        }
      }
    }
  } else {  // Positions evolve randomly
    arma::rowvec mod(n_pos);
    arma::rowvec prob_weight(n_pos);
    for (int tt = (n_dec - 2); tt >= 0; tt--) {
      reward = Rcpp::as<arma::cube>(Reward_(
          Rcpp::as<Rcpp::NumericMatrix>(Rcpp::wrap(path.slice(tt))), tt + 1));
#pragma omp parallel for private(policy, prob_weight, mod)
      for (std::size_t ii = 0; ii < n_path; ii++) {
        for (std::size_t pp = 0; pp < n_pos; pp++) {
          //  Primal values
          mod = primals.slice(tt + 1).row(ii) + mart.slice(tt).row(ii);
          policy = path_action(ii, pp, tt) - 1;
          prob_weight = control2.tube(pp, policy);
          primals(ii, pp, tt) =
              reward(ii, policy, pp) + arma::sum(prob_weight % mod);
          // Dual values
          mod = duals.slice(tt + 1).row(ii) + mart.slice(tt).row(ii);
          prob_weight = control2.tube(pp, 0);
          duals(ii, pp, tt) = reward(ii, 0, pp) + arma::sum(prob_weight % mod);
          for (std::size_t aa = 1; aa < n_action; aa++) {
            prob_weight = control2.tube(pp, aa);
            duals(ii, pp, tt) =
                std::max(reward(ii, aa, pp) + arma::sum(prob_weight % mod),
                         duals(ii, pp, tt));
          }
        }
      }
    }
  }
  return Rcpp::List::create(Rcpp::Named("primal") = primals,
                            Rcpp::Named("dual") = duals);
}

// [[Rcpp::export]]
XPtrImage magick_image_readbitmap_double(Rcpp::NumericVector x){
  Rcpp::IntegerVector dims(x.attr("dim"));
  return magick_image_bitmap(x.begin(), Magick::DoublePixel, dims[0], dims[1], dims[2]);
}