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

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
}

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
}

SEXP assignRawSymmetricMatrixFromEstimateRF(SEXP destination_, SEXP rowIndices_, SEXP columnIndices_, SEXP source_)
{
BEGIN_RCPP
	Rcpp::S4 destination = destination_;
	Rcpp::RawVector source = source_;
	Rcpp::RawVector destinationData = destination.slot("data");
	Rcpp::IntegerVector rowIndices = rowIndices_;
	Rcpp::IntegerVector columnIndices = columnIndices_;

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

	std::vector<int> markerRows, markerColumns;
	markerRows = Rcpp::as<std::vector<int> >(rowIndices);
	markerColumns = Rcpp::as<std::vector<int> >(columnIndices);
	if(countValuesToEstimate(markerRows, markerColumns) != (unsigned long long)source.size())
	{
		throw std::runtime_error("Mismatch between index length and source object size");
	}

	triangularIterator iterator(markerRows, markerColumns);
	R_xlen_t counter = 0;
	for(; !iterator.isDone(); iterator.next())
	{
		std::pair<int, int> markerPair = iterator.get();
		R_xlen_t markerRow = markerPair.first, markerColumn = markerPair.second;
		destinationData((markerColumn*(markerColumn-(R_xlen_t)1))/(R_xlen_t)2 + (markerRow - (R_xlen_t)1)) = source(counter);
		counter++;
	}
	return R_NilValue;
END_RCPP
}

    glmerResp::glmerResp(Rcpp::S4 xp) throw (std::runtime_error)
	: merResp(xp),
	  d_fam(SEXP(xp.slot("family"))),
	  d_n(       xp.slot("n")),
	  d_eta(     xp.slot("eta")) {
	updateWts();
    }

// construct cell type from equivalent class in R
CellType::CellType(unsigned id, const Rcpp::S4& type)
{
    mID = id;
    mName = Rcpp::as<std::string>(type.slot("name"));
    mSize = Rcpp::as<double>(type.slot("size"));
    mMinCycle = Rcpp::as<double>(type.slot("minCycle"));
    mCellTypeClass = type;
}

// [[Rcpp::export]]
Rcpp::NumericVector marginal_theta(Rcpp::S4 xmod) {
  RNGScope scope ;
  Rcpp::S4 model_(xmod) ;
  Rcpp::S4 model = clone(model_) ;
  Rcpp::S4 params=model.slot("mcmc.params") ;
  Rcpp::S4 chains(model.slot("mcmc.chains")) ;  
  int S = params.slot("iter") ;
  // Assume current values are the modes (in R, useModes(object) ensures this)
  // List modes = model.slot("modes") ;
  // NumericVector thetastar = as<NumericVector>(modes["theta"]) ;
  List modes = model.slot("modes") ;
  NumericVector theta_ = as<NumericVector>(modes["theta"]) ;
  NumericVector thetastar=clone(theta_) ;
  int K = thetastar.size() ;
  NumericVector p_theta(S) ;
  NumericVector muc = chains.slot("mu") ;
  NumericVector tau2c = chains.slot("tau2") ;
  NumericMatrix sigma2 = chains.slot("sigma2") ;
  NumericVector tauc = sqrt(tau2c) ;
  NumericVector tmp(K) ;

  IntegerMatrix Z = chains.slot("z") ;
  IntegerVector zz ;

  double tau2_tilde ;
  NumericVector sigma2_tilde(K) ;

  // this should be updated for each iteration
  NumericVector data_mean(K) ;
  IntegerVector nn(K) ;
  double post_prec;
  double tau_n;
  double mu_n;
  double w1;
  double w2;

  for(int s=0; s < S; ++s){
    zz = Z(s, _) ;
    model.slot("z") = zz ;
    nn = tableZ(K, zz) ;
    data_mean = compute_means(model) ;
    tau2_tilde = 1/tau2c[s] ;
    sigma2_tilde = 1.0/sigma2(s, _) ;
    //tmp = dnorm(thetastar, muc[s], tauc[s]) ;
    double prod = 1.0;
    for(int k = 0; k < K; ++k) {
      post_prec = tau2_tilde + sigma2_tilde[k] * nn[k];
      tau_n = sqrt(1/post_prec);
      w1 = tau2_tilde/post_prec;
      w2 = nn[k]*sigma2_tilde[k]/post_prec;
      mu_n = w1*muc[s] + w2*data_mean[k];
      tmp = dnorm(thetastar, mu_n, tau_n) ;
      prod = prod * tmp[k] ;
    }
    p_theta[s] = prod ;
  }
  return p_theta ;
}

void parse_gb_feature_table(
    Rcpp::S4& gb_feature,
    const std::vector<std::string>& feature_string,
    std::string& accession )
{   
    std::smatch m;
    std::string this_line("");
    std::vector<std::string> merged_lines;
    // merge qualifier values that occupy more than one line into "merged_lines"
    // Define qualifier position: Optional FT + blanks followed by forward slash
    static const std::regex QUAL_POS("^(FT)?\\s+/");
    static const std::regex TRANS("^translation=");
    std::vector<std::string>::const_iterator s_it = std::begin(feature_string);
    for ( ; s_it != std::end(feature_string); s_it++ ) {
        std::regex_search( std::begin(*s_it), std::end(*s_it), m, QUAL_POS ); 
        if ( m[0].matched ) {
            merged_lines.push_back(this_line);
            // If we have an embl file trim FT, whitespace and /
            // If we have a gbk file trim whitespace and /
            this_line = *s_it;
            this_line = ltrim(this_line, "FT");
            this_line = ltrim(this_line, " /");
            this_line = rtrim(this_line, " \n\r\t");
        } else {
            std::string tmp = *s_it;
            tmp = ltrim(tmp, "FT");
            tmp = ltrim(tmp, " /");
            tmp = rtrim(tmp, " \n\r\t");
            if ( std::regex_search( this_line, m, TRANS ) ) {
                // don't add whitespace when joining translations
                this_line += tmp;
            } else {
                this_line += " " + tmp;
            }
        }
        // Push the last qualifier
        if ( s_it == std::end(feature_string) - 1 ) {
            merged_lines.push_back(this_line);
        }
    }
    
    // Get key
    int begin_key = merged_lines[0].find_first_not_of(" ");
    int end_key = merged_lines[0].find_first_of( " ", begin_key );
    std::string key = merged_lines[0].substr( begin_key, end_key - 1 );
    gb_feature.slot("key") = key;
    
    // Get location
    Rcpp::S4 gb_location = Rcpp::S4("gbLocation");
    std::string gb_base_span = merged_lines[0].substr( end_key, merged_lines[0].length() - end_key );
    parse_gb_location( gb_location, gb_base_span, accession );
    gb_feature.slot("location") = gb_location;
    
    merged_lines.erase( std::begin(merged_lines) );
    gb_feature.slot("qualifiers") = parse_gb_qualifiers(merged_lines);
}

SEXP rawSymmetricMatrixSubsetByMatrix(SEXP object_, SEXP index_)
{
BEGIN_RCPP
	Rcpp::S4 object;
	try
	{
		object = object_;
	}
	catch(...)
	{
		throw std::runtime_error("Input object must be an S4 object");
	}

	Rcpp::RawVector data;
	try
	{
		data = Rcpp::as<Rcpp::RawVector>(object.slot("data"));
	}
	catch(...)
	{
		throw std::runtime_error("Slot [email protected] must be a raw vector");
	}

	Rcpp::NumericVector levels;
	try
	{
		levels = Rcpp::as<Rcpp::NumericVector>(object.slot("levels"));
	}
	catch(...)
	{
		throw std::runtime_error("Slot [email protected] must be a numeric vector");
	}

	Rcpp::IntegerMatrix index;
	try
	{
		index = index_;
	}
	catch(...)
	{
		throw std::runtime_error("Input index must be an integer matrix");
	}

	int nIndices = index.nrow();
	Rcpp::NumericVector output(nIndices);
	for(int row = 0; row < nIndices; row++)
	{
		R_xlen_t i = index(row, 0);
		R_xlen_t j = index(row, 1);
		if(i > j) std::swap(i, j);
		output(row) = levels[data[(j*(j-(R_xlen_t)1))/(R_xlen_t)2 + i-(R_xlen_t)1]];
	}
	return output;
END_RCPP
}

    reModule::reModule(Rcpp::S4 xp)
	: d_L(     S4(clone(SEXP(xp.slot("L"))))),
	  d_Lambda(S4(clone(SEXP(xp.slot("Lambda"))))),
	  d_Zt(          S4(xp.slot("Zt"))),
	  d_Lind(           xp.slot("Lind")),
	  d_lower(          xp.slot("lower")),
	  d_theta(          xp.slot("theta")),
	  d_u0(             d_L.n),
	  d_incr(           d_L.n),
	  d_u(              d_L.n),
	  d_cu(             d_L.n) {
	d_Ut = (CHM_SP)NULL;
    }

double cpp_luT(
	Rcpp::S4 xR,
	Rcpp::S4 dR) {

	const bool BisVCL=1;
	const int ctx_id = INTEGER(xR.slot(".context_index"))[0]-1;

	std::shared_ptr<viennacl::matrix<T> > vclX = getVCLptr<T>(xR.slot("address"), BisVCL, ctx_id);
	std::shared_ptr<viennacl::vector_base<T> > vclD = getVCLVecptr<T>(dR.slot("address"), BisVCL, ctx_id);


	return(luT(*vclX, *vclD));
}

//Returns the value of any(([email protected] >= length([email protected])) & [email protected] != as.raw(255))
SEXP checkRawSymmetricMatrix(SEXP rawSymmetric_)
{
BEGIN_RCPP
	Rcpp::S4 rawSymmetric = rawSymmetric_;
	Rcpp::NumericVector levels = Rcpp::as<Rcpp::NumericVector>(rawSymmetric.slot("levels"));
	Rcpp::RawVector data = Rcpp::as<Rcpp::RawVector>(rawSymmetric.slot("data"));
	R_xlen_t size = data.size(), levelsSize = levels.size();
	for(R_xlen_t i = 0; i < size; i++)
	{
		if(data[i] >= levelsSize && data[i] != 0xff) return Rcpp::wrap(true);
	}
	return Rcpp::wrap(false);
END_RCPP
}

// [[Rcpp::export]]
Rcpp::NumericVector p_mu_reduced(Rcpp::S4 xmod) {
  RNGScope scope ;
  Rcpp::S4 model(xmod) ;
  Rcpp::S4 mcmcp = model.slot("mcmc.params") ;
  Rcpp::S4 chains = model.slot("mcmc.chains") ;
  Rcpp::S4 hypp = model.slot("hyperparams") ;
  List modes = model.slot("modes") ;
  //
  //
  NumericVector x = model.slot("data") ;    
  int K = hypp.slot("k") ;
  int S = mcmcp.slot("iter") ;  
  int N = x.size() ;
  //
  NumericVector p_=as<NumericVector>(modes["mixprob"]) ;
  NumericVector theta_=as<NumericVector>(modes["theta"]) ;
  NumericVector mu_=as<NumericVector>(modes["mu"]) ;
  NumericVector pstar = clone(p_) ;
  NumericVector mustar = clone(mu_) ;
  NumericVector thetastar = clone(theta_) ;
  IntegerMatrix Z = chains.slot("z") ;
  NumericVector tau2chain = chains.slot("tau2") ;
  IntegerVector zz(N) ;

  IntegerVector nn(K) ;
  NumericVector mu0=hypp.slot("mu.0") ;
  double mu_0 = mu0[0] ;
  NumericVector tau2_0 = hypp.slot("tau2.0") ;
  double tau20_tilde = 1.0/tau2_0[0] ;

  double mu_k ;
  NumericVector tau2_tilde = 1.0/tau2chain ;
  double w1 ;
  double w2 ;
  double tau_k ;
  NumericVector p_mu(S) ;
  NumericVector tmp(1) ;
  for(int s = 0; s < S; ++s){
    zz = Z(s, _) ;
    nn = tableZ(K, zz) ;
    double total = 0.0 ;
    double thetabar = 0.0 ;    
    for(int k = 0; k < K; k++) total += nn[k] ;    
    for(int k = 0; k < K; k++) thetabar += nn[k] * thetastar[k] / total ;
    double post_prec = tau20_tilde + K*tau2_tilde[s] ;
    w1 = tau20_tilde/post_prec ;
    w2 = K*tau2_tilde[s]/post_prec ;
    mu_k =  w1*mu_0 +  w2*thetabar ;
    tau_k = sqrt(1.0/post_prec) ;
    tmp  = dnorm(mustar, mu_k, tau_k) ;
    p_mu[s] = tmp[0] ;
  }
  return p_mu ;
}

    merResp::merResp(Rcpp::S4 xp)  throw (std::runtime_error)
	: d_y(                    xp.slot("y")),
	  d_weights(        xp.slot("weights")),
	  d_offset(          xp.slot("offset")),
	  d_mu(                 d_y.size(), 0.),
	  d_sqrtrwt(                d_y.size()),
	  d_wtres(                  d_y.size()),
	  d_sqrtXwt(d_y.size(), d_offset.size()/d_y.size()) {
	int n = d_y.size(), os = d_offset.size();
	if (d_mu.size() != n || d_weights.size() != n || d_sqrtrwt.size() != n)
	    throw std::runtime_error("y, mu, sqrtrwt and weights slots must have equal lengths");
	if (os < 1 || os % n)
	    throw std::runtime_error("length(offset) must be a positive multiple of length(y)");
	init();
    }

SEXP constructDissimilarityMatrix(SEXP object, SEXP clusters_)
{
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"));
	int nMarkers = markers.size();
	std::vector<double> levelsCopied = Rcpp::as<std::vector<double> >(levels);
	
	std::vector<int> permutation(nMarkers);
	for(int i = 0; i < nMarkers; i++) permutation[i] = i;
	return constructDissimilarityMatrixInternal(&(data(0)), levelsCopied, nMarkers, clusters_, 0, permutation);
END_RCPP
}

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
}