C++ NumericMatrix类代码示例

// LMM scan of a single chromosome with interactive covariates
// this version should be fast but requires more memory
// (since we first expand the genotype probabilities to probs x intcovar)
// and this one allows weights for the individuals (the same for all phenotypes)
//
// genoprobs = 3d array of genotype probabilities (individuals x genotypes x positions)
// pheno     = matrix with one column of numeric phenotypes
//             (no missing data allowed)
// addcovar  = additive covariates (an intercept, at least)
// intcovar  = interactive covariates (should also be included in addcovar)
// eigenvec  = matrix of transposed eigenvectors of variance matrix
// weights   = vector of weights (really the SQUARE ROOT of the weights)
//
// output    = vector of log likelihood values
//
// [[Rcpp::export]]
NumericVector scan_pg_onechr_intcovar_highmem(const NumericVector& genoprobs,
                                              const NumericMatrix& pheno,
                                              const NumericMatrix& addcovar,
                                              const NumericMatrix& intcovar,
                                              const NumericMatrix& eigenvec,
                                              const NumericVector& weights,
                                              const double tol=1e-12)
{
    const unsigned int n_ind = pheno.rows();
    if(pheno.cols() != 1)
        throw std::range_error("ncol(pheno) != 1");
    const Dimension d = genoprobs.attr("dim");
    const unsigned int n_pos = d[2];
    if(n_ind != d[0])
        throw std::range_error("nrow(pheno) != nrow(genoprobs)");
    if(n_ind != addcovar.rows())
        throw std::range_error("nrow(pheno) != nrow(addcovar)");
    if(n_ind != intcovar.rows())
        throw std::range_error("nrow(pheno) != nrow(intcovar)");
    if(n_ind != weights.size())
        throw std::range_error("nrow(pheno) != length(weights)");
    if(n_ind != eigenvec.rows())
        throw std::range_error("ncol(pheno) != nrow(eigenvec)");
    if(n_ind != eigenvec.cols())
        throw std::range_error("ncol(pheno) != ncol(eigenvec)");

    // expand genotype probabilities to include geno x interactive covariate
    NumericVector genoprobs_rev = expand_genoprobs_intcovar(genoprobs, intcovar);

    // pre-multiply everything by eigenvectors
    genoprobs_rev = matrix_x_3darray(eigenvec, genoprobs_rev);
    NumericMatrix addcovar_rev = matrix_x_matrix(eigenvec, addcovar);
    NumericMatrix pheno_rev = matrix_x_matrix(eigenvec, pheno);

    // multiply everything by the (square root) of the weights
    // (weights should ALREADY be the square-root of the real weights)
    addcovar_rev = weighted_matrix(addcovar_rev, weights);
    pheno_rev = weighted_matrix(pheno_rev, weights);
    genoprobs_rev = weighted_3darray(genoprobs_rev, weights);

    // regress out the additive covariates
    genoprobs_rev = calc_resid_linreg_3d(addcovar_rev, genoprobs_rev, tol);
    pheno_rev = calc_resid_linreg(addcovar_rev, pheno_rev, tol);

    // now the scan, return RSS
    NumericMatrix rss = scan_hk_onechr_nocovar(genoprobs_rev, pheno_rev, tol);

    // 0.5*sum(log(weights)) [since these are sqrt(weights)]
    double sum_logweights = sum(log(weights));

    // calculate log likelihood
    NumericVector result(n_pos);
    for(unsigned int pos=0; pos<n_pos; pos++)
        result[pos] = -(double)n_ind/2.0*log(rss[pos]) + sum_logweights;

    return result;
}

// [[Rcpp::export]]
NumericVector colMeansC(NumericMatrix x) {
  
  // number of rows and columns
  int nCol = x.ncol();
  int nRow = x.nrow();
  
  // temporary variable of size nrow(x) to store column values in
  NumericVector nVal(nRow);
  
  // initialize output vector
  NumericVector out(nCol);
  
  // loop over each column
  for (int i = 0; i < nCol; i++) {
    
    // values in current column
    nVal = x(_, i);
    
    // store mean of current 'nVal' in 'out[i]'
    out[i] = mean(nVal);
  }
  
  return out;
}

// Summing log row
// [[Rcpp::export]]
NumericVector logplusRowMatC(const NumericMatrix & x) {

    int n = x.nrow();

    NumericVector out(n);

    for(int i = 0; i < n; i++) {

        out(i) = logplusvecC( ExtRowMatC(x, i) );

    }

    return out;

}

//----------------------------------------------------------------
// [[Rcpp::export]]
LogicalMatrix matchspatial(NumericMatrix locs, NumericMatrix x, int ncells, int state) {
  /*
  
    x is vector of simulated locations
    locs are cell locations (row, column numbers) to match
    ncells is the number of cells (distance) used to match location
    state is the value to check for
    return matrix z contains TRUE if match otherwise FALSE 
    -1 taken off locs row and col numbers to convert to C++ indexing
  */
  int n = x.nrow();
  int m = x.ncol();
  LogicalMatrix z(n,m); 
  int nlocs = locs.nrow();
    for(int k = 0; k < nlocs; k++) {
      for(int ii = imax(-ncells,-(locs(k,0) - 1)); ii <= imin(n-(locs(k,0) - 1),ncells); ii++) {
          for(int jj= imax(-ncells,-(locs(k,1) - 1)); jj <= imin(m-(locs(k,1) - 1),ncells); jj++) {
            if(x((locs(k,0) + ii - 1),(locs(k,1) + jj - 1)) == state) 
               z((locs(k,0) +ii  - 1),(locs(k,1) + jj - 1)) = true;
          }
        }
      }
  return(z);
}

void setlogP_C2(NumericMatrix logP,NumericVector NegLL,NumericMatrix cbars,NumericMatrix G3,NumericMatrix LLconst){
 
    int n = logP.nrow(), k = logP.ncol();
    int l1 =cbars.ncol();
    
      arma::mat logP2(logP.begin(), n, k, false); 
      NumericVector cbartemp=cbars(0,_);  
      NumericVector G3temp=G3(0,_);  
    
      arma::colvec cbarrow(cbartemp.begin(),l1,false);
      arma::colvec G3row(G3temp.begin(),l1,false);
    
        
      for(int i=0;i<n;i++){
        cbartemp=cbars(i,_);  
        G3temp=G3(i,_);  

        logP(i,1)=logP(i,0)-NegLL(i)+0.5*arma::as_scalar(cbarrow.t() * cbarrow)+arma::as_scalar(G3row.t() * cbarrow);

      LLconst(i,0)=NegLL(i)-arma::as_scalar(G3row.t() * cbarrow);
      }

 
}

void update_marginal_distr(ListOf<NumericMatrix> Q, NumericMatrix res, int k, int n){
  arma::mat out(res.begin(), k, n, false);
  for(int t=1; t<=n-1; t++){
    // calculate rowsums of Q[t]
    arma::mat B(Q[t].begin(), k, k, false);
    arma::colvec rowsums = sum(B, 1);
    // assign rowsums to res(_, t-1)
    out.col(t-1) += rowsums;
  }
  // calculate colsums of Q[n-1]
  arma::mat B(Q[n-1].begin(), k, k, false);
  arma::rowvec colsums = sum(B, 0);
  arma::colvec temp = arma::vec(colsums.begin(), k, false, false);
  out.col(n-1) += temp;
}

// [[Rcpp::export]]
double lnL_scalar(NumericVector y, List x, NumericMatrix groups,
           NumericMatrix beta, List gamma, NumericVector delta,
           NumericMatrix z, double alpha, double sigma2,
           List bs_beta) {
  double out;
  
  //Calculate mean value for each observation
  double sum_sq = 0;
  for (int i = 0; i < y.size(); i++) {
    double mu = alpha;
    for (int k = 0; k < beta.nrow(); k++) {
      NumericMatrix bs_beta_k = bs_beta[k];
      NumericMatrix x_k = x[k];
      for (int j = 0; j < x_k.ncol(); j++) {
        double mu_tmp = 0;
        for (int jj = 0; jj < beta.ncol(); jj++) {
          mu_tmp += bs_beta_k(j, jj) * beta(k, jj);
        }
        mu += x_k(i, j) * mu_tmp;
      }
    }
    for (int k = 1; k < z.ncol(); k++) {
      mu += delta[k] * z(i, k);
    }
    for (int q = 0; q < gamma.size(); q++) {
      NumericVector gamma_q = gamma[q];
      int group_set = groups(i, q);
      mu += gamma_q(group_set);
    }
    sum_sq += (y(i) - mu) * (y(i) - mu);
  }

  //Calculate loglik
  out = (- (double) y.size() / 2.0) * log(2.0 * PI * sigma2) - (1.0 / (2.0 * sigma2)) * sum_sq;
  return out;
}

    double nlmerResp::updateMu(Rcpp::NumericVector const &gamma) throw(std::runtime_error) {
	int n = d_y.size();
#ifdef USE_RCPP_SUGAR
	Rcpp::NumericVector gam = gamma + d_offset;
#else
	NumericVector gam(d_offset.size());
	std::transform(gamma.begin(), gamma.end(), d_offset.begin(),
		       gam.begin(), std::plus<double>());
#endif
	double *gg = gam.begin();

	for (int p = 0; p < d_pnames.size(); p++) {
	    std::string pn(d_pnames[p]);
	    Rcpp::NumericVector pp = d_nlenv.get(pn);
	    std::copy(gg + n * p, gg + n * (p + 1), pp.begin());
	}
	NumericVector rr = d_nlmod.eval(SEXP(d_nlenv));
	if (rr.size() != n)
	    throw std::runtime_error("dimension mismatch");
	std::copy(rr.begin(), rr.end(), d_mu.begin());
	NumericMatrix rrg = rr.attr("gradient");
	std::copy(rrg.begin(), rrg.end(), d_sqrtXwt.begin());
	return updateWrss();
    }

// **********************************************************//
//              Likelihood evaluation of Timepart            //
// **********************************************************//
// [[Rcpp::export]]
double Timepartsum (NumericMatrix mumat, double sigma_tau,
                   IntegerVector senders, NumericVector timestamps){
    int D = senders.size();
    double timesum = 0;
    for (unsigned int d = 0; d < D; d++) {
        int a_d = senders[d] - 1;
        timesum += R::dlnorm(timestamps[d], mumat(d, a_d), sigma_tau, TRUE);
        for (unsigned int i = 0; i < mumat.ncol(); i++) {
        if (i != a_d) {
            timesum += R::plnorm(timestamps[d], mumat(d, i), sigma_tau, FALSE, TRUE);
    	  }
        }
    }
    return timesum;
}

// optimal proposal
void LinearGaussian::transition_and_weight(Tree* tree, NumericMatrix & xparticles, NumericVector & logweights, int step){
  RNGScope scope;
  int nparticles = xparticles.nrow();
  NumericVector trans_noise = rnorm(nparticles, 0, 1);
  double var = sigma2*tau2 / (sigma2 + tau2);
  double sd = sqrt(var);
  double Minus1Div2Sd = -1.0 / (2*(sigma2 + tau2));
  for (int iparticle = 0; iparticle < nparticles; iparticle++){
    int j = tree->l_star(iparticle);
    xparticles(iparticle, 0) = var * ((rho * tree->x_star(j, 0)) / sigma2 + observations(step, 0) / tau2);
    xparticles(iparticle, 0) = xparticles(iparticle, 0) + sd *  trans_noise(iparticle);
    logweights(iparticle) = Minus1Div2Sd * (rho * tree->x_star(j, 0) - observations(step, 0)) *
      (rho * tree->x_star(j, 0) - observations(step, 0));
  }
}

NumericMatrix updatePPLCpp(NumericMatrix x, NumericMatrix dm, int idx) {
  int ncolx = x.ncol(), nrowx = x.nrow(), ncoldm = dm.ncol(), nrowdm = dm.nrow();
  NumericVector d(nrowx, 0.0000);
  NumericMatrix x2(1, ncolx);
  NumericMatrix res(nrowdm, ncoldm);
  
  // Get the data of the distance matrix
  // This is needed so that the object passed to 'dm' is not replaced in the 
  // global environment.
  for (int i = 0; i < nrowdm; i++) {
    for (int j = 0; j < ncoldm; j++) {
      res(i, j) = dm(i, j);
    }
  }
  
  // Get the coordinates of the new point
  idx -= 1;
  for (int i = 0; i < ncolx; i++) {
    x2(0, i) = x(idx, i);
  }
  
  // Calculate distances
  for (int i = 0; i < nrowx; i++) {
    for (int j = 0; j < ncolx; j++) {
      d[i] += (x[nrowx * j + i] - x2[j]) * (x[nrowx * j + i] - x2[j]);
    }
    d[i] = pow(d[i], 0.5);
  }
  
  // replace the values in the distance matrix
  for (int i = 0; i < nrowx; i++) {
    res(i, idx) = d[i];
    res(idx, i) = d[i];
  }
  return (res);
}

// [[Rcpp::export]]
NumericVector splitpointC(NumericMatrix Dm, NumericVector x, Function f) {
  int J = Dm.ncol() - 1;
  NumericVector ux = unique(x).sort(), chisq(J+2), chisqtemp(J+3, -1.0);
  int nux = ux.size(), id = -1;
  for (int i=0; i <= nux-2; i++) {
    chisq = f(Dm, x >= (ux[i]+ux[i+1])/2);
    if (chisq[0] >= chisqtemp[1]) {
       chisqtemp[0] = (ux[i]+ux[i+1])/2;
       for (int j=1; j<J+3; j++) {
         chisqtemp[j] = chisq[j-1];
       }
    }  
  }
  return chisqtemp; 
}

// Function for pair indexing
// [[Rcpp::export]]
std::vector<int> indexPairs(NumericMatrix & X,
                            const String op = "==",
                            const double ref = 0){

  // Index pairs by operator and reference
  std::vector<int> index;
  int nfeats = X.nrow();
  for(int i = 1; i < nfeats; i++){
    for(int j = 0; j < i; j++){

      if(op == "==" || op == "="){
        if(X(i, j) == ref){
          index.push_back(j * nfeats + i + 1);
        }
      }else if(op == ">"){
        if(X(i, j) > ref){
          index.push_back(j * nfeats + i + 1);
        }
      }else if(op == ">="){
        if(X(i, j) >= ref){
          index.push_back(j * nfeats + i + 1);
        }
      }else if(op == "<"){
        if(X(i, j) < ref){
          index.push_back(j * nfeats + i + 1);
        }
      }else if(op == "<="){
        if(X(i, j) <= ref){
          index.push_back(j * nfeats + i + 1);
        }
      }else if(op == "!="){
        if(X(i, j) != ref){
          index.push_back(j * nfeats + i + 1);
        }

      }else if(op == "all"){

        index.push_back(j * nfeats + i + 1);

      }else{

        stop("Operator not found.");
      }
    }
  }

  return index;
}

// **********************************************************//
//              Construct the history of interaction         //
// **********************************************************//
// [[Rcpp::export]]
List History(List edge, NumericVector timestamps, NumericMatrix p_d, IntegerVector node, int d, double timeunit) {
  int nIP = p_d.ncol();
  int A = node.size();
  List IPmat(nIP);
  for (int IP = 0; IP < nIP; IP++) {
  	List IPlist_IP(3);
  	for (unsigned int l = 0; l < 3; l++){
  		NumericMatrix IP_l(A, A);
  		IPlist_IP[l] = IP_l;
  	}
  		IPmat[IP] = IPlist_IP;
  }
    double time1 = timestamps[d-2]-384*timeunit;
	double time2 = timestamps[d-2]-96*timeunit;
    double time3 = timestamps[d-2]-24*timeunit;
    int iter = which_num(time1, timestamps);

	  for (int i = iter-1; i < (d-1); i++) {
	    List document2 = edge[i];
	    int sender = document2[0];
	    IntegerVector receiver = document2[1];
	    double time = timestamps[i];
	    for (unsigned int r = 0; r < receiver.size(); r++){
	       for (int IP = 0; IP < nIP; IP++) {
  			  List IPlist_IP = IPmat[IP];
               if (time < time2) {
                   NumericMatrix IP_l = IPlist_IP[2];
                   IP_l(sender-1, receiver[r]-1) += p_d(i, IP);
                   IPlist_IP[2] = IP_l;
               } else {
                   if (time >= time2 && time < time3) {
                       NumericMatrix IP_l = IPlist_IP[1];
                       IP_l(sender-1, receiver[r]-1) += p_d(i, IP);
                       IPlist_IP[1] = IP_l;
                   } else {
                       if (time >= time3) {
                       NumericMatrix IP_l = IPlist_IP[0];
                       IP_l(sender-1, receiver[r]-1) += p_d(i, IP);
                       IPlist_IP[0] = IP_l;
                       }
                   }
               }
			    IPmat[IP] = IPlist_IP;
  			}
         }
	  }
	return IPmat;
}

// [[Rcpp::export]]
NumericMatrix dist_by_closeness(NumericMatrix mat) {
	int n = mat.nrow();

	NumericMatrix dist(n, n);

	for(int i = 0; i < n - 1; i ++) {
		for(int j = i+1; j < n; j ++) {
			dist(i, j) = closeness(mat(i, _), mat(j, _));
			dist(j, i) = dist(i, j);
		}
	}
	for(int i = 0; i < n; i ++) {
		dist(i, i) = 0;
	}
	return dist;
}