C++ IntegerVector::size方法代码示例

SEXP assignRawSymmetricMatrixDiagonal(SEXP destination_, SEXP indices_, SEXP source_)
{
BEGIN_RCPP
	Rcpp::S4 destination = destination_;
	Rcpp::RawVector source = source_;
	Rcpp::RawVector destinationData = destination.slot("data");
	Rcpp::IntegerVector indices = indices_;

	if(&(source(0)) == &(destinationData(0)))
	{
		throw std::runtime_error("Source and destination cannot be the same in assignRawSymmetricMatrixDiagonal");
	}

	if((indices.size()*(indices.size()+(R_xlen_t)1))/(R_xlen_t)2 != source.size())
	{
		throw std::runtime_error("Mismatch between index length and source object size");
	}
	for(R_xlen_t column = 0; column < indices.size(); column++)
	{
		for(R_xlen_t row = 0; row <= column; row++)
		{
			R_xlen_t rowIndex = indices[row];
			R_xlen_t columnIndex = indices[column];
			if(rowIndex > columnIndex)
			{
				std::swap(rowIndex, columnIndex);
			}
			destinationData((columnIndex*(columnIndex-(R_xlen_t)1))/(R_xlen_t)2+rowIndex-(R_xlen_t)1) = source((column*(column+(R_xlen_t)1))/(R_xlen_t)2 + row);
		}
	}
END_RCPP
}

// [[Rcpp::export]]
Rcpp::RawVector magick_image_write( XPtrImage input, Rcpp::CharacterVector format, Rcpp::IntegerVector quality,
                                    Rcpp::IntegerVector depth, Rcpp::CharacterVector density, Rcpp::CharacterVector comment){
  if(!input->size())
    return Rcpp::RawVector(0);
  XPtrImage image = copy(input);
#if MagickLibVersion >= 0x691
  //suppress write warnings see #74 and #116
  image->front().quiet(true);
#endif
  if(format.size())
    for_each ( image->begin(), image->end(), Magick::magickImage(std::string(format[0])));
  if(quality.size())
    for_each ( image->begin(), image->end(), Magick::qualityImage(quality[0]));
  if(depth.size())
    for_each ( image->begin(), image->end(), Magick::depthImage(depth[0]));
  if(density.size()){
    for_each ( image->begin(), image->end(), Magick::resolutionUnitsImage(Magick::PixelsPerInchResolution));
    for_each ( image->begin(), image->end(), Magick::densityImage(Point(density[0])));
  }
  if(comment.size())
    for_each ( image->begin(), image->end(), Magick::commentImage(std::string(comment.at(0))));
  Magick::Blob output;
  writeImages( image->begin(), image->end(),  &output );
  Rcpp::RawVector res(output.length());
  std::memcpy(res.begin(), output.data(), output.length());
  return res;
}

SEXP hclustThetaMatrix(SEXP mpcrossRF_, SEXP preClusterResults_)
{
BEGIN_RCPP
	Rcpp::List preClusterResults = preClusterResults_;
	bool noDuplicates;
	R_xlen_t preClusterMarkers = countPreClusterMarkers(preClusterResults_, noDuplicates);
	if(!noDuplicates)
	{
		throw std::runtime_error("Duplicate marker indices in call to hclustThetaMatrix");
	}

	Rcpp::S4 mpcrossRF = mpcrossRF_;
	Rcpp::S4 rf = mpcrossRF.slot("rf");

	Rcpp::S4 theta = rf.slot("theta");
	Rcpp::RawVector data = theta.slot("data");
	Rcpp::NumericVector levels = theta.slot("levels");
	Rcpp::CharacterVector markers = theta.slot("markers");
	if(markers.size() != preClusterMarkers)
	{
		throw std::runtime_error("Number of markers in precluster object was inconsistent with number of markers in mpcrossRF object");
	}
	R_xlen_t resultDimension = preClusterResults.size();
	//Allocate enough storage. This symmetric matrix stores the *LOWER* triangular part, in column-major storage. Excluding the diagonal. 
	Rcpp::NumericVector result(((resultDimension-(R_xlen_t)1)*resultDimension)/(R_xlen_t)2);
	for(R_xlen_t column = 0; column < resultDimension; column++)
	{
		Rcpp::IntegerVector columnMarkers = preClusterResults(column);
		for(R_xlen_t row = column + 1; row < resultDimension; row++)
		{
			Rcpp::IntegerVector rowMarkers = preClusterResults(row);
			double total = 0;
			R_xlen_t counter = 0;
			for(R_xlen_t columnMarkerCounter = 0; columnMarkerCounter < columnMarkers.size(); columnMarkerCounter++)
			{
				R_xlen_t marker1 = columnMarkers[columnMarkerCounter]-(R_xlen_t)1;
				for(R_xlen_t rowMarkerCounter = 0; rowMarkerCounter < rowMarkers.size(); rowMarkerCounter++)
				{
					R_xlen_t marker2 = rowMarkers[rowMarkerCounter]-(R_xlen_t)1;
					R_xlen_t column = std::max(marker1, marker2);
					R_xlen_t row = std::min(marker1, marker2);
					Rbyte thetaDataValue = data((column*(column+(R_xlen_t)1))/(R_xlen_t)2 + row);
					if(thetaDataValue != 0xFF)
					{
						total += levels(thetaDataValue);
						counter++;
					}
				}
			}
			if(counter == 0) total = 0.5;
			else total /= counter;
			result(((resultDimension-(R_xlen_t)1)*resultDimension)/(R_xlen_t)2 - ((resultDimension - column)*(resultDimension-column-(R_xlen_t)1))/(R_xlen_t)2 + row-column-(R_xlen_t)1) = total;
		}
	}
	return result;
END_RCPP
}

static Rcpp::NumericMatrix x_tilde(Rcpp::NumericMatrix X, 
				   Rcpp::IntegerVector nk,
				   Rcpp::IntegerMatrix groups, 
				   Rcpp::NumericVector d_cur, 
				   Rcpp::NumericMatrix eta_cur)
{
  int K = nk.size();
  int n_tot = X.nrow();
  int p = X.ncol();
  int L = groups.ncol();
  Rcpp::NumericMatrix result(n_tot, p * K);

  int idx = 0;
  for (int k = 0; k < K; k++) {
    int n = nk[k];
    for (int j = 0; j < p; j++) {
      //calculate sum for column j
      double sum = 0.0;
      for (int l = 0; l < L; l++) {
	if (elem(groups, j, l)) {
	  sum += d_cur[l] * elem(eta_cur, l, k);
	} 
      }
      
      //multiply column j in submatrix k of X with sum
      for (int i = 0; i < n; i++) {
	elem(result, idx + i, p * k + j) = elem(X, idx + i, j) * sum;
      }
    }
    idx += n;
  }
  return result;
}

static void print(Rcpp::IntegerVector a)
{
  for (int i = 0; i < a.size(); i++) {
    printf("%d ", a[i]);
  }
  putchar('\n');
}

static double bic_logistic(Rcpp::NumericMatrix X, 
			   Rcpp::NumericVector y, 
			   Rcpp::NumericMatrix beta_new, 
			   double eps, 
			   Rcpp::IntegerVector nk)
{
  int n_tot = X.nrow();
  int p = X.ncol();
  int K = nk.size();
    
  int idx = 0;
  double ll = 0.0;
  for (int k = 0; k < K; k++) {
    int n = nk[k];
    for (int i = 0; i < n; i++) {
      double lp = 0.0;
      for (int j = 0; j < p; j++) {
	lp += elem(X, idx+i, j) * elem(beta_new, j, k);
      }
      ll += y[idx+i] * lp - log(1.0 + exp(lp));
    }
    idx += n;
  }
  
  double bic = -2.0 * ll + df(beta_new, eps) * log(n_tot);
  return bic;
}

static double bic_linear(Rcpp::NumericMatrix X, 
			 Rcpp::NumericVector y, 
			 Rcpp::NumericMatrix beta_new, 
			 double eps, 
			 Rcpp::IntegerVector nk)
{
  int n_tot = X.nrow();
  int p = X.ncol();
  int K = nk.size();
 

  /*calculate SSe*/
  double SSe = 0.0;  
  int idx = 0;
  
  for (int k = 0; k < K; k++) {
    int n = nk[k];
    for (int i = 0; i < n; i++) {
      double Xrow_betacol = 0.0;
      for (int j = 0; j < p; j++) {
	Xrow_betacol += elem(X, idx+i, j) * elem(beta_new, j, k);
      }
      SSe += pow(y[idx+i] - Xrow_betacol, 2);
    }
    idx += n;
  }
  
  double ll = -n_tot / 2.0 * (log(SSe) - log(n_tot) + log(2.0 * M_PI) + 1);
  double bic = -2 * ll + df(beta_new, eps) * log(n_tot);

  return bic;
}

// Calculate CPO for spatial Copula Cox PH
arma::vec LinvSpCopulaCox(Rcpp::NumericVector tobs, Rcpp::IntegerVector delta, Rcpp::NumericVector Xbeta, 
                          Rcpp::NumericVector h, Rcpp::NumericVector d, arma::mat Cinv, arma::vec z){
  int n = delta.size();
  arma::vec res(n);
  for(int i=0; i<n; ++i){
    double Cii = Cinv(i,i);
    double s2i = 1.0/Cii;
    double mui = -s2i*( arma::dot(Cinv.col(i), z) - Cii*z[i] );
    double si = std::sqrt(s2i);
    double Fi = Foft(tobs[i], h, d, Xbeta[i]);
    double PinvFyi = Rf_qnorm5(Fi, 0, 1, true, false);
    double newzi = (PinvFyi-mui)/si;
    if(delta[i]==0){
      double St = 1.0-Rf_pnorm5(newzi, 0, 1, true, false);
      //Rprintf( "St=%f\n", St );
      res(i) = 1.0/St;
    } else {
      double fi = foft(tobs[i], h, d, Xbeta[i]);
      double diff = -0.5*std::pow(newzi, 2) + 0.5*std::pow(PinvFyi, 2);
      double ft = (1.0/si*std::exp(diff)*fi);
      res(i) = 1.0/ft;
      //Rprintf( "ft=%f\n", ft );
    }
  }
  return(res);
}

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

        }
}

void check_topset(const Rcpp::IntegerVector& top) {
    for (size_t t=1; t<top.size(); ++t) {
        if (top[t] < top[t-1]) { 
            throw std::runtime_error("numbers of top genes must be sorted"); 
        }
    }
    return;
}

//[[Rcpp::export]]
Rcpp::CharacterVector edgeIdCpp (Rcpp::IntegerMatrix edge, std::string type) {
    Rcpp::IntegerVector ances = getAnces(edge);
    Rcpp::IntegerVector desc = getDesc(edge);
    int nedge;
        
    if (type == "tip" || type == "internal") {
	Rcpp::IntegerVector tips = tipsFast(ances);
	nedge = tips.size();
	Rcpp::IntegerVector ans = match(tips, desc);
	if (type == "tip") {
            Rcpp::IntegerVector tmpAnces(nedge);
            Rcpp::IntegerVector tmpDesc(nedge);
	    for (int j = 0; j < nedge; j++) {
                tmpAnces[j] = ances[ans[j]-1];
                tmpDesc[j] = desc[ans[j]-1];
	    }
            Rcpp::CharacterVector c1(nedge);
            c1 = edgeIdCppInternal(tmpAnces, tmpDesc);
            return c1;
	}
	else if (type == "internal") {
	    int allEdges = ances.size();
	    Rcpp::IntegerVector idEdge = Rcpp::seq_len(allEdges);
	    Rcpp::IntegerVector intnd = Rcpp::setdiff(idEdge, ans);
	    nedge = intnd.size();
            Rcpp::IntegerVector tmpAnces(nedge);
            Rcpp::IntegerVector tmpDesc(nedge);
	    for (int j = 0; j < nedge; j++) {
                tmpAnces[j] = ances[intnd[j]-1];
                tmpDesc[j] = desc[intnd[j]-1];
            }
            Rcpp::CharacterVector c1(nedge);
            c1 = edgeIdCppInternal(tmpAnces, tmpDesc);
            return c1;
	}
    }
    else {
        nedge = ances.size();
        Rcpp::IntegerVector tmpAnces = ances;
        Rcpp::IntegerVector tmpDesc = desc;
        Rcpp::CharacterVector c1(nedge);
        c1 = edgeIdCppInternal(tmpAnces, tmpDesc);
        return c1;
    }
    return "";
}

// Calculate mk = sum_i I(M(ti)=k), 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 Getmk(Rcpp::IntegerVector& mk, const Rcpp::IntegerVector& Mt){
  int n = Mt.size();
  std::fill(mk.begin(), mk.end(), 0);
  for (int i=0; i<n; ++i){
    int k = Mt[i];
    mk[k] +=1;
  }
}

Rcpp::IntegerVector match_and_substract(Rcpp::IntegerVector x, Rcpp::IntegerVector yInt) {
    int y(1);
    y = yInt[0];
    for (unsigned k=0; k < x.size(); ++k) {
	if (x[k] > y)
	    x[k] = x[k] - 1;
    }
    return x;
}

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

        void fill(Rcpp::IntegerVector idx, SEXP A){
            viennacl::vector_range<viennacl::vector_base<T> > v_sub(*shptr.get(), r);
#ifdef UNDEF            
            Eigen::Matrix<T, Eigen::Dynamic, 1> Am;
            Am = Rcpp::as<Eigen::Matrix<T, Eigen::Dynamic, 1> >(A);
            
            for(int i = 0; i < idx.size(); i++) {
                v_sub(idx[i]) = Am(i);
            }
#endif            
        }