C++ rcpp::NumericVector类代码示例

// [[Rcpp::export]]
NumericVector computeAllLogLkhd(
	NumericVector observedCases, NumericVector expectedCases, 
	List nearestNeighborsList, int nZones, String logLkhdType) {

	NumericVector allLogLkhd(nZones);
	int nAreas = expectedCases.size(), C = sum(observedCases), 
	N = sum(expectedCases), index = 0;
	int nNeighbors = 0;
	double cz = 0.0, nz = 0.0;

	for (int i = 0; i < nAreas; ++i) {
		cz = 0;
		nz = 0;

		Rcpp::NumericVector nearestNeighbors = nearestNeighborsList[i];
		nNeighbors =  nearestNeighbors.size();

		// For each area's nearest neighbors
		for(int j = 0; j < nNeighbors; ++j) { 
			// Watch off by 1 vector indexing in C as opposed to R
			cz += observedCases[nearestNeighbors[j]-1];
			nz += expectedCases[nearestNeighbors[j]-1];

			if (logLkhdType=="poisson") {
				allLogLkhd[index] = poissonLogLkhd(cz, nz, N, C);
			} else if (logLkhdType == "binomial" ) {
				allLogLkhd[index] = binomialLogLkhd(cz, nz, N, C);
			}
			index++;
		}
	}
	return allLogLkhd;
}

// [[Rcpp::export]]
Rcpp::NumericVector LearnJ(Rcpp::NumericVector x, 
                           Rcpp::NumericVector y,
                           Rcpp::NumericVector extent,
                           Rcpp::NumericVector image_labels){
  
  assert(x.size() == y.size());
  assert(y.size() == image_labels.size());
  const int size_x = extent[1] - extent[0] + 1;
  const int size_y = extent[3] - extent[2] + 1;
  int x0 = extent[0];
  int y0 = extent[2];
  std::cout << "Size of grid: "<< size_x << " " << size_y << std::endl;
  // Construct the grid of image_labels from the 1D data
  std::vector<std::vector<int> > grid(size_y, std::vector<int>(size_x,0));
  
  for (int i = 0 ; i < x.size() ; i++){
		int posy = y[i] - y0;
		int posx = x[i] - x0;
    //std::cout << "i: " << i << " x, y: ";
		//std::cout << posy << ", ";
    //std::cout << posx << std::endl;
    grid[posy][posx] = image_labels[i];
  }
  return 0;
}

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
}

SEXP rawSymmetricMatrixToDist(SEXP object)
{
BEGIN_RCPP
	Rcpp::S4 rawSymmetric = object;
	Rcpp::NumericVector levels = Rcpp::as<Rcpp::NumericVector>(rawSymmetric.slot("levels"));
	Rcpp::CharacterVector markers = Rcpp::as<Rcpp::CharacterVector>(rawSymmetric.slot("markers"));
	Rcpp::RawVector data = Rcpp::as<Rcpp::RawVector>(rawSymmetric.slot("data"));
	R_xlen_t size = markers.size(), levelsSize = levels.size();

	Rcpp::NumericVector result(size*(size - 1)/2, 0);
	int counter = 0;
	for(R_xlen_t row = 0; row < size; row++)
	{
		for(R_xlen_t column = row+1; column < size; column++)
		{
			int byte = data(column * (column + (R_xlen_t)1)/(R_xlen_t)2 + row);
			if(byte == 255) result(counter) = std::numeric_limits<double>::quiet_NaN();
			else result(counter) = levels(byte);
			counter++;
		}
	}
	result.attr("Size") = (int)size;
	result.attr("Labels") = markers;
	result.attr("Diag") = false;
	result.attr("Upper") = false;
	result.attr("class") = "dist";
	return result;
END_RCPP
}

//[[Rcpp::export]]
extern "C" SEXP notbayesEstimation(SEXP Xs, SEXP Ys, SEXP Zs, SEXP Vs) {
Rcpp::NumericVector Xr(Xs);
Rcpp::NumericMatrix Yr(Ys);
Rcpp::NumericMatrix Zr(Zs);
Rcpp::NumericMatrix Vr(Vs);
int n = Yr.nrow(), k = Yr.ncol();
int l = Vr.nrow(), m = Vr.ncol();
arma::mat x(Xr.begin(),n,k,false); //fit2$sigma
arma::mat y(Yr.begin(),n,k,false); //fit2b$stdev.unscaled
arma::mat z(Zr.begin(),1,1,false); //min.variance.factor
arma::mat v(Vr.begin(),l,m,false); //fit2$df.residual


arma::mat sda(n,1);
arma::vec dof(n);
arma::vec sd(n);

sd  = sqrt(  (y%x)%(y%x) + as_scalar(z));
sda = sd/(x%y);
dof = v;
Rcpp::NumericVector Sd = Rcpp::wrap(sd);
Rcpp::NumericVector Sda = Rcpp::wrap(sda);
Rcpp::NumericVector DOF = Rcpp::wrap(dof);
   Sd.names() = Rcpp::List(Yr.attr("dimnames"))[0];
   Sda.names() = Rcpp::List(Yr.attr("dimnames"))[0];
 
return Rcpp::List::create( Rcpp::Named("SD") = Sd,
                           Rcpp::Named("DOF") =  DOF,
                           Rcpp::Named("sd.alpha") = Sda);
}

// [[Rcpp::export]]
Rcpp::NumericVector randomSample(Rcpp::NumericVector a, int n) {
    // clone a into b to leave a alone
    Rcpp::NumericVector b(n);
    __gnu_cxx::random_sample(a.begin(), a.end(), 
                             b.begin(), b.end(), randWrapper);
    return b;
}

void MapTematikQuantile::createMapTematikQuantile()
{
    Rcpp::NumericVector quantile ;
    QString command;
    try {
        command = QString("q <- quantile(dframe[[\"%1\"]], prob = seq(0, 1, length = (%2 + 1)), type = 7);").arg(var).arg(typeMap);
        rconn.parseEvalQ(command.toStdString());

        quantile =  rconn["q"];
    } catch (...) {

    }

    QList<int> temp[quantile.size()-1];

    for(int i=0; i<numvar.size(); i++){
        if(numvar[i] <= quantile[1]){
            temp[0].append(table->verticalHeaderItem(i)->text().toInt());
        }else{
            for(int j=2; j<quantile.size(); j++){
                if(numvar[i] > quantile[j-1] && numvar[i] <= quantile[j]){
                    temp[j-1].append(table->verticalHeaderItem(i)->text().toInt());
                }
            }
        }
    }

    QList<QList<int> > temp2;
    for(int i=0; i<quantile.size()-1; i++){
        temp2.append(temp[i]);
    }

    MapTematikConfig* configWidget = new MapTematikConfig(mviewResult,vv,rconn,temp2,var,typeMap.toInt());
    setupResultViewVariableTypeChooser("Quantile Map",var, temp2,quantile,configWidget);
}

//' rcpp_neutral_hotspots_ntests
//'
//' Performs repeated neutral tests to yield average distributions of both
//' hotspot values and spatial autocorrelation statistics.
//'
//' @param nbs An \code{spdep} \code{nb} object listing all neighbours of each
//' point. 
//' @param wts Weighting factors for each neighbour; must have same length as
//' nbs. 
//' @param nbsi List of matrices as returned from \code{get_nbsi}. each element
//' of which contains the i-th nearest neighbour to each point.
//' @param alpha Strength of spatial autocorrelation
//' @param sd0 Standard deviation of truncated normal distribution used to model
//' environmental variation (with mean of 1)
//' @param nt Number of successive layers of temporal and spatial autocorrelation
//' used to generate final modelled values
//' @param ntests Number of tests used to obtain average values
//' @param ac_type Character string specifying type of aucorrelation
//' (\code{moran}, \code{geary}, or code{getis-ord}).
//'
//' @return A matrix of dimension (size, 2), with first column containing
//' sorted and re-scaled hotspot values, and second column containing sorted and
//' re-scaled spatial autocorrelation statistics.
//'
// [[Rcpp::export]]
Rcpp::NumericMatrix rcpp_neutral_hotspots_ntests (Rcpp::List nbs, 
        Rcpp::List wts, Rcpp::List nbsi, double alpha, double sd0, int niters,
        std::string ac_type, bool log_scale, int ntests)
{
    const int size = nbs.size ();

    Rcpp::NumericMatrix hs1;
    Rcpp::NumericVector z (size), z1 (size), ac (size), ac1 (size);
    std::fill (ac.begin (), ac.end (), 0.0);
    std::fill (z.begin (), z.end (), 0.0);

    for (int n=0; n<ntests; n++)
    {
        hs1 = rcpp_neutral_hotspots (nbs, wts, nbsi, alpha, sd0, log_scale,
                niters, ac_type);

        z += hs1 (Rcpp::_, 0);
        ac += hs1 (Rcpp::_, 1);
    }

    Rcpp::NumericMatrix result (size, 2);
    result (Rcpp::_, 0) = z / (double) ntests;
    result (Rcpp::_, 1) = ac / (double) ntests;
    Rcpp::colnames (result) = Rcpp::CharacterVector::create ("z", "ac");

    return result;
}

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

        }
}

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

// kernel Dist function on a Grid
// [[Rcpp::export]]
Rcpp::NumericVector
KdeDist(const Rcpp::NumericMatrix & X
      , const Rcpp::NumericMatrix & Grid
      , const double                h
      , const Rcpp::NumericVector & weight
      , const bool printProgress
	) {
	const unsigned sampleNum = X.nrow();
	const unsigned dimension = Grid.ncol();
	const unsigned gridNum = Grid.nrow();
	// first = sum K_h(X_i, X_j), second = K_h(x, x), third = sum K_h(x, X_i)
	std::vector< double > firstValue;
	const double second = 1.0;
	std::vector< double > thirdValue;
	double firstmean;
	Rcpp::NumericVector kdeDistValue(gridNum);
	int counter = 0, percentageFloor = 0;
	int totalCount = sampleNum + gridNum;

	if (printProgress) {
		printProgressFrame(Rprintf);
	}

	firstValue = computeKernel< std::vector< double > >(
			X, X, h, weight, printProgress, Rprintf, counter, totalCount,
			percentageFloor);

	if (dimension <= 1) {
		thirdValue = computeKernel< std::vector< double > >(
				X, Grid, h, weight, printProgress, Rprintf, counter, totalCount,
				percentageFloor);
	}
	else {
		thirdValue = computeGaussOuter< std::vector< double > >(
				X, Grid, h, weight, printProgress, Rprintf, counter, totalCount,
				percentageFloor);
	}

	if (weight.size() == 1) {
		firstmean = std::accumulate(firstValue.begin(), firstValue.end(), 0.0) / sampleNum;
	}
	else {
		firstmean = std::inner_product(
				firstValue.begin(), firstValue.end(), weight.begin(), 0.0) / 
				std::accumulate(weight.begin(), weight.end(), 0.0);
	}

	for (unsigned gridIdx = 0; gridIdx < gridNum; ++gridIdx) {
		kdeDistValue[gridIdx] = std::sqrt(firstmean + second - 2 * thirdValue[gridIdx]);
	}

	if (printProgress) {
		Rprintf("\n");
	}

	return kdeDistValue;
}

    // redis "set a vector" -- without R serialization, without attributes, ...
    // this is somewhat experimental
    std::string setVector(std::string key, Rcpp::NumericVector x) {

        redisReply *reply = 
            static_cast<redisReply*>(redisCommand(prc_, "SET %s %b", 
                                                  key.c_str(), x.begin(), x.size()*szdb));
        std::string res(reply->str);                                                
        freeReplyObject(reply);
        return(res);
    }

    double lmerResp::updateMu(const Rcpp::NumericVector& gamma) {
#ifdef USE_RCPP_SUGAR
	d_mu = d_offset + gamma;
#else
	std::transform(gamma.begin(), gamma.end(), d_offset.begin(),
		       d_mu.begin(), std::plus<double>());
#endif
	return updateWrss();
    }

// [[Rcpp::export]]
arma::mat rmvnorm_arma(int n,
                  const arma::vec& mu,
                  const arma::mat& Sigma) {
    unsigned int p = Sigma.n_cols;
    Rcpp::NumericVector draw = Rcpp::rnorm(n*p);
    arma::mat Z = arma::mat(draw.begin(), n, p, false, true);
    arma::mat Y = arma::repmat(mu, 1, n).t()+Z * arma::chol(Sigma);
    return Y;
}

// Calculate lk, k=1, ..., M with m0=0;
// where h=(h0, h1, ..., hM) with h0=0 and d=(d0, d1, ..., dM) with d0=0, dM=R_PosInf
void Getlk(Rcpp::NumericVector& lk, const Rcpp::IntegerVector& Mt, int M1, Rcpp::NumericVector d, 
           const Rcpp::NumericVector& t, const Rcpp::NumericVector& Xbeta){
  int n = Mt.size();
  std::fill(lk.begin(), lk.end(), 0);
  for (int k=1; k<M1; ++k){
    for (int i=0; i<n; ++i){
      if(Mt[i]>=k) lk[k] += (std::min(d[k],t[i])-d[k-1])*std::exp(Xbeta[i]);
    }
  }
}