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

// [[Rcpp::export]]
Rcpp::List subsetCounts(Rcpp::IntegerVector counts, Rcpp::IntegerVector start, Rcpp::IntegerVector width, Rcpp::LogicalVector strand){
	if (start.length() != width.length() || start.length() != strand.length()) Rcpp::stop("provided vectors have different lengths...");
	int nr = start.length();
	int len = counts.length();
	int tot = 0;
	int* S = start.begin(); int* W = width.begin();
	for (int i = 0; i < nr; ++i){
		int s = S[i] - 1;
		int w = W[i]; 
		if (s < 0) Rcpp::stop("negative start positions are invalid");
		if (s + w > len) Rcpp::stop("range exceeds the lengths of the counts vector");
		tot += w;
	}
	
	Rcpp::IntegerVector res(tot); 
	Rcpp::IntegerVector nstart(nr);
	Rcpp::IntegerVector nend(nr);
	int* R = res.begin(); int* C = counts.begin(); int* ST = strand.begin();
	int* NS = nstart.begin(); int* NE = nend.begin();
	int currpos = 0;
	for (int i = 0; i < nr; ++i){
		NS[i] = currpos + 1;
		int w = W[i];
		if (ST[i]) std::copy(C + S[i]-1, C + S[i]-1 + w, R + currpos);
		else std::reverse_copy(C + S[i]-1, C + S[i]-1 + w, R + currpos);
		currpos += w;
		NE[i] = currpos;
	}
	return List::create(_("counts")=res, _("starts")=nstart, _("ends")=nend);
}

// [[Rcpp::export]]
Rcpp::NumericMatrix testColPost(Rcpp::NumericMatrix post, Rcpp::List m2u, int nthreads){
    Rcpp::IntegerVector values = Rcpp::as<Rcpp::IntegerVector>(m2u["values"]);
    Rcpp::IntegerVector map = Rcpp::as<Rcpp::IntegerVector>(m2u["map"]);
    if (post.ncol() != map.length()) Rcpp::stop("posteriors doesn't match with m2u");
    
    Rcpp::NumericMatrix smallerPost(post.nrow(), values.length());
    Vec<double> foo; NMPreproc preproc(asVec(values), asVec(map), foo);
    collapsePosteriors_core(asMat(smallerPost), asMat(post), preproc);
    return smallerPost;
}

SEXP getAllFunnels(SEXP Rmpcross)
{
	char* stackmem;
	{
		std::string error;
		{
			int nFounders;
			Rcpp::RObject mpcross_ = Rmpcross;
			bool valid = validateMPCross(mpcross_, nFounders, error, true, false, false);
			if(!valid)
			{
				goto signal_error;
			}
			Rcpp::List mpcross = Rmpcross;
			Rcpp::DataFrame pedigree(mpcross["pedigree"]);
			Rcpp::IntegerVector id = mpcross["id"];
			int nFinals = id.length();
			std::vector<int> fid = Rcpp::as<std::vector<int> >(mpcross["fid"]);
			Rcpp::IntegerMatrix output(id.length(), nFounders);
			std::vector<int> nIntercrossingGenerations;
			nIntercrossingGenerations.resize(nFinals, 0);
			//get number of intercrossing generations
			bool ok = intercrossingGenerations(pedigree, nFounders, id, nIntercrossingGenerations);
			if(!ok)
			{
				error = "Problem determining number of intercrossing generations";
				goto signal_error;
			}
			//now get the actual funnels from the pedigree
			int funnel[16];
			for(int i = 0; i < id.length(); i++)
			{
				ok = getFunnel(id[i], pedigree, fid, nIntercrossingGenerations[i], funnel, pedigree.nrows(), nFounders);
				if(!ok)
				{
					std::stringstream ss;
					ss << "Problem with pedigree, for individual number " << (i+1) << ", having id " << id[i];
					error = ss.str();
					goto signal_error;
				}
				for(int j = 0; j < nFounders; j++) output(i, j) = funnel[j];
			}
			return output;
		}
	signal_error:
		stackmem = (char*)alloca(error.size() + 4);
		memset(stackmem, 0, error.size() + 4);
		memcpy(stackmem, error.c_str(), error.size());
	}
	Rf_error(stackmem);
	return R_NilValue;
}

// [[Rcpp::export]]
Rcpp::IntegerVector testSortCounts(Rcpp::IntegerVector v){
    Rcpp::IntegerVector res(v.length());
    Vec<int> v2 = asVec(v);
    Vec<int> v3 = asVec(res);
    sortCounts(v2, v3);
    return res;
}

// [[Rcpp::export]]
Rcpp::IntegerMatrix quantileNorm(Rcpp::IntegerMatrix mat, Rcpp::IntegerVector ref, int nthreads=1, int seed=13){
    if (mat.nrow() != ref.length()) Rcpp::stop("incompatible arrays...");
    if (!std::is_sorted(ref.begin(), ref.end())) Rcpp::stop("ref must be sorted");
    int ncol = mat.ncol();
    int nrow = mat.nrow();
    //allocate new matrix
    Rcpp::IntegerMatrix res(nrow, ncol);
    Mat<int> oldmat = asMat(mat); 
    Mat<int> newmat = asMat(res);
    Vec<int> ref2 = asVec(ref);
    //allocate a seed for each column
    std::seed_seq sseq{seed};
    std::vector<std::uint32_t> seeds(ncol);
    sseq.generate(seeds.begin(), seeds.end());
    
    #pragma omp parallel num_threads(nthreads)
    {
        std::vector<std::pair<int, int> > storage(nrow);//pairs <value, index>
        #pragma omp for 
        for (int col = 0; col < ncol; ++col){
            std::mt19937 gen(seeds[col]);
            qtlnorm(oldmat.getCol(col), ref2, newmat.getCol(col), storage, gen);
        }
    }
    
    res.attr("dimnames") = mat.attr("dimnames");
    return res;
}

// [[Rcpp::export]]
Rcpp::IntegerVector countInSubset(Rcpp::IntegerVector counts, Rcpp::IntegerVector start, Rcpp::IntegerVector width){
	if (start.length() != width.length()) Rcpp::stop("provided vectors have different lengths...");
	int nr = start.length();
	int len = counts.length();
	Rcpp::IntegerVector res(nr); 
	int* R = res.begin(); int* C = counts.begin();
	int* S = start.begin(); int* W = width.begin();
	for (int i = 0; i < nr; ++i){
		int s = S[i] - 1;
		int w = W[i]; 
		if (s < 0) Rcpp::stop("negative start positions are invalid");
		if (s + w > len) Rcpp::stop("range exceeds the lengths of the counts vector");
		R[i] = sum(C + s, w);
	}
	
	return res;
}

//Function to summarize results from BeQTL and return a dataframe
// [[Rcpp::export]]
Rcpp::DataFrame SumRes(const arma::mat  & cormat, const arma::mat & errmat, const Rcpp::DataFrame SnpDF, const Rcpp::DataFrame Genedf, const int samplesize, const double tcutoff){
  Rcpp::Rcout<<"Generating Tmat"<<std::endl;
  arma::mat tmat = sqrt(samplesize-2)*(cormat/(sqrt(1-square(cormat))));
  Rcpp::Rcout<<"Finding strong t"<<std::endl;
  arma::uvec goods = find(abs(tmat)>tcutoff);
  Rcpp::Rcout<<"Generating matrix index"<<std::endl;
  arma::umat goodmat = Ind(tmat.n_rows,goods);
  Rcpp::Rcout<<"Subsetting tmat"<<std::endl;
  Rcpp::Rcout<<"This many good results"<<goodmat.n_rows<<std::endl;
  arma::vec tvec = tmat(goods);
  Rcpp::Rcout<<"Subsetting Errmat"<<std::endl;
  arma::vec errvec = errmat(goods);
  Rcpp::Rcout<<"Generating SNP and Gene lists"<<std::endl;
  Rcpp::IntegerVector GoodGenes = Rcpp::wrap(arma::conv_to<arma::ivec>::from(goodmat.col(0)));
  Rcpp::IntegerVector GoodSNPs = Rcpp::wrap(arma::conv_to<arma::ivec>::from(goodmat.col(1)));

//Subset SNP anno

  Rcpp::CharacterVector SNPnames = SnpDF["rsid"];
  SNPnames = SNPnames[GoodSNPs];
  arma::ivec SNPchr = Rcpp::as<arma::ivec>(SnpDF["Chrom"]);
  SNPchr = SNPchr(goodmat.col(1));
  arma::ivec SNPpos = Rcpp::as<arma::ivec>(SnpDF["Pos"]);
  SNPpos = SNPpos(goodmat.col(1));
//Subset Geneanno
  Rcpp::CharacterVector GeneNames = Genedf["Symbol"];
  GeneNames = GeneNames[GoodGenes];
  arma::ivec Genechr = Rcpp::as<arma::ivec>(Genedf["Chrom"]);
  Genechr = Genechr(goodmat.col(0));
  arma::ivec Genestart = Rcpp::as<arma::ivec>(Genedf["Start"]);
  Genestart = Genestart(goodmat.col(0));
  arma::ivec Genestop = Rcpp::as<arma::ivec>(Genedf["Stop"]);
  Genestop = Genestop(goodmat.col(0));

  arma::ivec CisDist(GoodGenes.length());
  Rcpp::Rcout<<"Calculating Cisdist"<<std::endl;
  CisDist = arma::min(arma::join_cols(abs(Genestop-SNPpos),abs(Genestart-SNPpos)),1);
  Rcpp::Rcout<<"CisDist Calculated"<<std::endl;
  CisDist.elem(find(Genechr!=SNPchr)).fill(-1);
  return  Rcpp::DataFrame::create(Rcpp::Named("SNP")=SNPnames,
                                  Rcpp::Named("Gene")=GeneNames,
                                  Rcpp::Named("t-stat")=Rcpp::wrap(tvec),
                                  Rcpp::Named("err")=Rcpp::wrap(errvec),
                                  Rcpp::Named("CisDist")=Rcpp::wrap(CisDist));

}

//' The Gillespie algorithm for simulating continuous time Markov chains.
//'
//' This function is called by \code{\link{Simulate.EventSim}} 
//'
//' @param start_states an integer vector giving the start states for each patient
//' @param start_times a numeric vector containing the calendar time at which each patient is recruited 
//' @param n_state an integer, the number of states on the underlying DAG 
//' @param rates A list of \code{n_state * n_state * length(patientEndTime)} transition rates
//' @param patientEndTimes The ith \code{n_state*n_state} block of rates represents the transition
//' rate matrix until patient time patientEndTimes[i]
//' @param calendarEndTimes The ith \code{n_state*n_state} block of rates represents the transition
//' rate matrix until calendar time calendarEndTimes[i]
//' @param shape The shape parameter for the edges in the DAG. A matrix where shape[i,j] is the Weibull
//' shape parameter for the edge from node i to node j. If there is no edge between i and j
//' then shape[i,j] = 0.
//' @param resetEdges A vector (from1,to1,from2,to2,...) of edges for which if a subject traverses 
//' one of these edges then patient time switches are reset. See \code{SetIsResetEdge.ProgressionGraph} for further 
//' details
//' @param duration The total (patient) time each subject is to be simulated 
//' @return A data frame with columns "id", "state" and "patient_transition_time" listing the (patient) time at
//' which patients transition into new states. For example
//' 
//' id state patient_transition_time
//' 1   1      0.0
//' 1   2      1.4
//' 
//' patient 1 transitions to state 1 at time 0 (w.r.t. patient time) and then at time 1.4 transitions
//' to state 2
//[[Rcpp::export]]
Rcpp::DataFrame gillespie(Rcpp::IntegerVector start_states,Rcpp::NumericVector start_times,
                          const int n_state,Rcpp::NumericVector rates, Rcpp::NumericVector patientEndTimes,
                          Rcpp::NumericVector calendarEndTimes,Rcpp::NumericMatrix shape, Rcpp::NumericVector
                          resetEdges, const double duration){
    //set up output vectors
    std::vector<int> Id;
    std::vector<int> transition_to_states;
    std::vector<double> transition_times;
    
    //Random number generator
    Rcpp::RNGScope scope;
    
    //lots of validation and set up here
    Transitions transitions(n_state,rates,patientEndTimes,calendarEndTimes,shape,resetEdges);
    
    //for each subject
    for(R_len_t subject_id=0; subject_id != start_states.length(); ++subject_id){
        
        //Initialize patient-specific variables
        int current_Id = subject_id+1;
        int current_state = start_states[subject_id];
                
        SubjectTime subjectTime(start_times[subject_id]);
              
        //store subject recruitment in output vectors
        Id.push_back(current_Id);
        transition_to_states.push_back(current_state);
        transition_times.push_back(0);
    
        //get intial transition from current state         
        double out_time=0; 
        double time_left_before_switch = transitions.getNextSwitch(subjectTime);
        int proposed_new_state = transitions.ProposeNewState(out_time,subjectTime, current_state);
        
            
        //keep going until, no more rate switching and the rate leaving current state = 0
        while(!(out_time == INFINITY && time_left_before_switch == INFINITY ) &&
                subjectTime.getCurrentPatientTime()  <=  duration){
            
            //std::cout << time_left_before_switch << " " << out_time << " " << subjectTime.getCurrentPatientTime() 
            //          << " " << proposed_new_state << std::endl;
                        
            if(out_time >= time_left_before_switch ||
               subjectTime.getCurrentPatientTime() + out_time > duration){ //subject switch rate matrices before transition
                 subjectTime.IncreasePatientTime(time_left_before_switch);   
                            
            }
            else{ //subject transitions at time current_patient_time + out_time
                subjectTime.IncreasePatientTime(out_time);
                
                int old_state = current_state;
               
                current_state = proposed_new_state;
               
                if(current_state != old_state){
                    Id.push_back(current_Id);
                    transition_to_states.push_back(current_state);
                    transition_times.push_back(subjectTime.getCurrentPatientTime());
                }
                //are we reseting the patient time switches (i.e. have we crossed 
                //an edge with isResetEdge = TRUE?)
                if(transitions.resetPatientTimeForSwitches(old_state,current_state)){
                    subjectTime.ResetPatientSwitches();
                }
                          
                
            }
            
             proposed_new_state = transitions.ProposeNewState(out_time,subjectTime, current_state);
             time_left_before_switch = transitions.getNextSwitch(subjectTime );
             
//.........这里部分代码省略.........