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

void NetSim::prop_spike(size_t id, double t) {
    size_t c_id = id-1;
    NumericVector sp = net[c_id];
    sp.push_back(t);
    net[c_id] = sp;

    for(size_t con_i=0; con_i<cons[c_id].size(); con_i++) {
//        cout << "id: " << id <<  " c_id: " << c_id << " con_i: " << con_i << " i: " << cons[c_id][con_i].nid << "\n";
        queue_of_spikes[ cons[c_id][con_i].nid ].push_back(TSynSpike(t, cons[c_id][con_i].syn_id));
    }
}

List homogeneousPoolCoxCD1(double r, double c, double tauc, double dt, int ndt, int n) {
	// Clock-driven generation of homogeneously correlated spike trains
	// (Cox processes)
	// r = rate (Hz)
	// c = total correlation strength (in [0,1])
	// tauc = correlation time constant (ms)
	// dt = timestep (ms)
	// ndt = number of timesteps
	// n = number of neurons
	// spike = array of 0/1 integers, each row is a timestep, each column is a neuron
    List net(n);
    for(size_t ni=0; ni<n; ni++) {
        net[ni] = NumericVector();
    }
	double sigmar,sigma,s,lambda,x,mu;
	int i,j;
	
	// Correction of mu and sigma
	sigmar=sqrt(c*r/(2.*tauc*0.001));
	invRectifiedGaussian(r,sigmar,&mu,&sigma);
	
	// Simulation
	x=0.;
	lambda=exp(-dt/tauc);
	s=sigma*sqrt(1-exp(-2.*dt/tauc));
	mu=mu*dt*0.001;
	s=s*dt*0.001;
    NumericVector nrm = rnorm(ndt, 0, s);
	for(j=0;j<ndt;j++) {
		x=x*lambda+nrm[j];
        NumericVector coins = runif(n);
		for(i=0;i<n;i++)
			if (coins[i]<x+mu) {
                NumericVector sp = as<NumericVector>(net[i]);
                sp.push_back(j);
                net[i] = sp;
			}
	}
    return net;
}

/* Function to run breadth-first search, returning only sets of connected 
   components that reside on the boundary of the districts */
List bsearch_boundary(List aList,
		      arma::vec boundary)
{

  /* Inputs to function:
     aList: adjacency list

     boundary: vector of boundary element indicators (as arma)
   */

  // Get indices of boundary units
  arma::uvec boundary_indices = find(boundary == 1);

  // Container - outputted of breadth search, a list
  List bsearch;

  // Container - partition vector, gets added to bsearch when queue is empty
  NumericVector partition;

  // Set mark vector - ledger of which indices have been reached
  NumericVector mark(aList.size());

  // Set queue vector
  NumericVector q;

  // Initialize breadth search with first element in boundary_indices
  mark(boundary_indices(0)) = boundary_indices(0);
  partition.push_back(boundary_indices(0));
  q = aList(boundary_indices(0));

  // Initialize objects inside loop
  int u; bool in_part; NumericVector adj_u; int i; int v; 

  // Begin do{} loop - run until number of elements in boundary_indices is 0
  do{

    // Begin while{} loop - run until q is empty
    while(q.size() > 0){
      
      // Dequeue first element in queue
      u = q(0);

      // Mark that element in ledger
      mark(u) = u;

      // Check if element is in the partition - add to partition if false
      in_part = is_true(any(partition == u));
      if(in_part == false){
	partition.push_back(u);
      }
      
      // Get adjacency vector for unit u
      adj_u = aList(u);

      // Loop through elements of adj_u, add to queue and mark if not reached
      if(adj_u.size() > 0){
	
	// Start loop
	for(i = 0; i < adj_u.size(); i++){
	  
	  // Reach element v
	  v = adj_u(i);

	  /* Check if already reached - if false, mark, add to partition, and
	     add to queue */
	  if(is_true(any(mark == v)) == FALSE){
	    mark(v) = v;
	    partition.push_back(v);
	    q.push_back(v);
	  }

	}

      }

      // Erase dequeued element from queue when done searching
      q.erase(q.begin());

    }

    // Handling an empty queue
    if(q.size() == 0){

      /* First, find boundary units that are in the reached partition and
	 remove them from boundary_units vector */
      for(i = boundary_indices.n_elem - 1; i >= 0; i--){
	if(is_true(any(partition == boundary_indices(i))) == TRUE){
	  boundary_indices.shed_row(i);
	}
      }
      
      // Store the partition, clear partition vector
      bsearch.push_back(partition);
      partition.erase(partition.begin(), partition.end());

      // Re-initialize breadth search from new starting value if nonempty
      if(boundary_indices.n_elem > 0){
	q = aList(boundary_indices(0));
//.........这里部分代码省略.........

// Function to cut edges of adjacency list probabilistically - Step 2 of swMH
List cut_edges(List aList_con,
	       double eprob)
{

  /* Inputs to function:
     aList_con: adjacency list within cong district

     eprob: edgecut probability (transformed into 1-eprob in function)
   */

  // Create threshold
  double threshold_prob = 1 - eprob;

  // Define lists to store cut-edge and uncut-edge vectors
  List aList_uncut(aList_con.size());
  List aList_cut(aList_con.size());

  // Initialize inside loop
  int i; NumericVector cc_vec_i_all; NumericVector cc_vec_i;
  arma::vec draws;

  // Define list to store output of both lists

  // Loop through elements of aList_con
  for(i = 0; i < aList_con.size(); i++){

    // Extract i'th vector in list
    cc_vec_i_all = aList_con(i);

    // Subset cc_vec_i to elements > i
    cc_vec_i = cc_vec_i_all[cc_vec_i_all > i];

    // For each element in vector, take random draw from [0,1] uniform
    draws = runif(cc_vec_i.size());

    // Create container vectors of cut and uncut edges
    NumericVector cut;
    NumericVector uncut;

    // Loop through elements of cc_vec_i and compare to entry in draws
    for(int j = 0; j < cc_vec_i.size(); j++){
      
      // Compare to threshold_prob - if draws < thresh, cut edge, else uncut
      if(draws(j) < threshold_prob){
	cut.push_back(cc_vec_i(j));
      } else{
	uncut.push_back(cc_vec_i(j));
      }

    }

    /* Here - look at lines 1201-1212 in original code. Modifying original
       alConnected to remove edges that are cut, but isn't this just the 
       uncut list (which will be aList_postcut? Skipping this bit for now */
    
    // Store vectors in container lists
    aList_uncut(i) = uncut;
    aList_cut(i) = cut;

  }
  
  // Add ties to aList_uncut, aList_cut
  List aList_uncut_bd = add_ties(aList_uncut);
  List aList_cut_bd = add_ties(aList_cut);

  // Return contents
  List out;
  out["connectedlist"] = aList_uncut_bd;
  out["cutedgelist"] = aList_cut_bd;
  
  return out;

}

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
List inseq(mat M, bool adjust=true)
{
  int i, m, trun;
  mat mu=mean(M);
  //center the rows
  M.each_row() -= mu;
  int n=M.n_rows, p=M.n_cols;
  //Dtm is the vector of det(Sig)'s
  NumericVector Dtm;
  //gam_0 and gam_1 are for gam_2m and gam_2m+1, resp.
  mat gam0(p,p), gam1(p,p), Gam(p,p), Sig(p,p), Sig1(p,p), Gamadj(p,p), eigvec(p,p);
  //for adjustment, set initial increment in Gam=0
  Gamadj.zeros();
  //store the eigenvalues and eigenvectors of each Gam
  vec eigval(p),eigvalneg(p);
  double dtm;
  int sn= n/2; 
  for (m=0; m<n/2; m++)
  {
    gam0.zeros(); gam1.zeros();
    //calculate gam_2m (gam0) and gam_2m+1 (gam1)
    for(i=0; i<(n-2*m);i++) gam0+=trans(M.row(i))*M.row(i+2*m);
    gam0=gam0/n;
    for(i=0; i<(n-2*m-1);i++) gam1+=trans(M.row(i))*M.row(i+2*m+1);
    gam1=gam1/n;
    //Gam_m=gam_2m+gam_2m+1, then symmetrize
    Gam=gam0+gam1; Gam=(Gam+Gam.t())/2;
    
    if (m==0) Sig=-gam0+2*Gam;
    else Sig=Sig+2*Gam;
    
    if (eig_sym(Sig)(0)>0)
    {
      sn=m;
      break;
    }
  }
  if (sn>n/2-1) 
  {
    stop("Not enough samples.");
  }
  Dtm=det(Sig);
  for (m=sn+1; m<n/2; m++)
  {
    gam0.zeros(); gam1.zeros();
    //calculate gam_2m (gam0) and gam_2m+1 (gam1)
    for(i=0; i<(n-2*m);i++) gam0+=trans(M.row(i))*M.row(i+2*m);
    gam0=gam0/n;
    for(i=0; i<(n-2*m-1);i++) gam1+=trans(M.row(i))*M.row(i+2*m+1);
    gam1=gam1/n;
    //Gam_m=gam_2m+gam_2m+1, then symmetrize
    Gam=gam0+gam1; Gam=(Gam+Gam.t())/2;
    
    //Sig_m=Sig_m-1+2Gam_m
    Sig1=Sig+2*Gam;
    dtm=det(Sig1);
    //if dtm1>dtm, continue
    if (dtm<=Dtm(m-sn-1)) break;
    //update Sig
    Sig=Sig1;
    //record dtm
    Dtm.push_back(dtm);

    //to adjust the original Sig, subtract the negative part of Gam
    if (adjust) 
    {
      //calculate eigenvalues and eigenvectors of Gam
      eig_sym(eigval,eigvec,Gam);
      eigvalneg=eigval;
      eigvalneg.elem(find(eigvalneg>0)).zeros();
      Gamadj-=eigvec*diagmat(eigvalneg)*eigvec.t();
    }
  }
  trun = Dtm.size()-1+sn;
  List res; res["Sig"]=Sig; res["Dtm"]=Dtm; res["trunc"]= trun; res["sn"]=sn;
  if (adjust) res["Sigadj"]=Sig+2*Gamadj;
  return res;
}