C++ NumericMatrix::ncol方法代码示例

// [[Rcpp::export]]
NumericMatrix rectify(NumericMatrix m){
  NumericMatrix out(m.nrow(), m.ncol());
  for(int i = 0; i< m.ncol(); i++){
      out(_, i) = clamp(0, m(_, i), INFINITY);
    };
  return out;
}

// [[Rcpp::export]]
NumericMatrix sortMatrix(NumericMatrix bipartiteAdjMatrix) {
	int numRows = bipartiteAdjMatrix.nrow();
	int numCols = bipartiteAdjMatrix.ncol();

	std::vector<myPair> colSums(numCols);
	std::vector<myPair> rowSums(numRows);

	NumericMatrix sortedMatrix(bipartiteAdjMatrix.nrow(), bipartiteAdjMatrix.ncol());

	//calculate row and column sums
	for(int rowIdx = 0; rowIdx < numRows; rowIdx++) {
		rowSums[rowIdx].second = rowIdx;
		for(int colIdx = 0; colIdx < numCols; colIdx++) {
			colSums[colIdx].second = colIdx;

			rowSums[rowIdx].first += bipartiteAdjMatrix(rowIdx,colIdx);
			colSums[colIdx].first += bipartiteAdjMatrix(rowIdx,colIdx);
		}
	}

	//sort with custom comparator method
	std::sort(rowSums.begin(), rowSums.end(), comparator_r);
	std::sort(colSums.begin(), colSums.end(), comparator_r);

	//create sorted matrix by cross-referencing sorted indexes
	for(int rowIdx = 0; rowIdx < numRows; rowIdx++) {
		for(int colIdx = 0; colIdx < numCols; colIdx++) {
			sortedMatrix(rowIdx, colIdx) = bipartiteAdjMatrix(rowSums[rowIdx].second, colSums[colIdx].second);
		}
	}

	return sortedMatrix;
}

// [[Rcpp::export]]
Rcpp::List   setlogP(NumericMatrix logP,NumericVector NegLL,NumericMatrix cbars,NumericMatrix G3) {

    int n = logP.nrow(), k = logP.ncol();
    int l1 =cbars.ncol();
//    int l2=cbars.nrow();
    
    arma::mat logP2(logP.begin(), n, k, false); 
    NumericVector cbartemp=cbars(0,_);  
    NumericVector G3temp=G3(0,_);  
    Rcpp::NumericMatrix LLconst(n,1);
    
    arma::colvec cbarrow(cbartemp.begin(),l1,false);
    arma::colvec G3row(G3temp.begin(),l1,false);
    
//    double v = arma::as_scalar(cbarrow.t() * cbarrow);
//    LLconst[j]<--t(as.matrix(cbars[j,1:l1]))%*%t(as.matrix(G3[j,1:l1]))+NegLL[j]    
        
    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);
    }
    
    
//    return logP;
    return Rcpp::List::create(Rcpp::Named("logP")=logP,Rcpp::Named("LLconst")=LLconst);
    
}

// [[Rcpp::export]]
NumericVector fstatsC(NumericMatrix pairMat, NumericMatrix mod, NumericMatrix mod0, NumericVector outcome) {
	int nrow = pairMat.nrow(); int ncol = pairMat.ncol();
	double ss, ss0;
	int df0 = mod0.ncol(); int df = df0+1; // alternative model always has +1 column
	
	arma::mat modc(mod.begin(), mod.nrow(), mod.ncol(), false);
	arma::mat mod0c(mod0.begin(), mod0.nrow(), mod0.ncol(), false);
	arma::colvec outcomec(outcome.begin(), outcome.size(), false);
	arma::vec fstats(nrow);
	arma::vec res = arma::zeros<arma::vec>(outcome.size());
	arma::vec res0 = arma::zeros<arma::vec>(outcome.size());
    
	try{ 
		res0 = outcomec - mod0c*arma::solve(mod0c, outcomec); // The residual for the null model remains the same
		ss0 = arma::as_scalar(arma::trans(res0)*res0);
	} catch(std::exception &ex) {
		stop("WTF???");
	}
	for(int i=0; i < nrow; i++){ // loop through all rows in pairMat
		//modc.insert_cols(mod.ncol(), pairMat(i,_)); can try this later, it's not working
		for(int j=0; j < pairMat.ncol(); j++){ // this loop is for copying the ith row of pairMat into the last column of modc
			modc(j,modc.n_cols-1) = pairMat(i,j); // Here we add the current row to the model
		}
		try{
			res = outcomec - modc*arma::solve(modc, outcomec); // Calculate the residual
			ss = arma::as_scalar(arma::trans(res)*res);
			fstats(i) = ((ss0 - ss)/(df-df0))/(ss/(ncol-df));
		} catch(std::exception &ex) {
			fstats(i) = NA_REAL;
		}
	}
	return Rcpp::wrap(fstats);
}

Rcpp::List   setlogP_C(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);
      }
    
    
    return Rcpp::List::create(Rcpp::Named("logP")=logP,Rcpp::Named("LLconst")=LLconst);
    
}

// [[Rcpp::export("multiKernel_cpp")]]
SEXP multiKernel(NumericMatrix Yr, NumericMatrix Zr, NumericMatrix Kr, double tau) {
    int n = Yr.nrow(), p = Yr.ncol(), q = Zr.ncol();

    arma::mat K(Kr.begin(), n, n, false);       // reuses memory and avoids extra copy
    arma::mat Y(Yr.begin(), n, p, false);
    arma::mat Z(Zr.begin(), n, q, false);

    // Initialize matrices
    arma::vec tuning_vect(n);
    tuning_vect.fill(tau);
    arma::mat tuning_mat(n, n, fill::zeros);
    tuning_mat.diag() = tuning_vect;
    arma::mat weight_mat = inv(K + tuning_mat);

    arma::mat B = inv(trans(Z) * weight_mat * Z) * trans(Z) * weight_mat * Y;
    arma::mat alpha_mat = weight_mat * (Y - Z * B);

    double newLS = computeLeastSq(Y, K, alpha_mat, Z, B);
    double BIC = 2 * newLS + log(n) * as_scalar(accu(B > 0) + accu(alpha_mat > 0));

    return Rcpp::List::create(
               Rcpp::Named("alpha") = alpha_mat,
               Rcpp::Named("B") = B,
               Rcpp::Named("LS") = newLS,
               Rcpp::Named("BIC") = BIC);
}

// [[Rcpp::export]]
NumericVector smooth_nd(const NumericMatrix& grid_in, 
                        const NumericVector& z_in, 
                        const NumericVector& w_in_,
                        const NumericMatrix& grid_out, 
                        const NumericVector h) {

  if (grid_in.nrow() != z_in.size()) stop("Incompatible input lengths");
  if (grid_in.ncol() != grid_out.ncol()) stop("Incompatible grid sizes");
  if (h.size() != grid_in.ncol()) stop("Incorrect h length");

  int n_in = grid_in.nrow(), n_out = grid_out.nrow(), d = grid_in.ncol();
  NumericVector w_in = (w_in_.size() > 0) ? w_in_ : 
    rep(NumericVector::create(1), n_in);
  NumericVector z_out(n_out), w_out(n_out);

  for (int i = 0; i < n_in; ++i) {
    for(int j = 0; j < n_out; ++j) {
      double w = 1;
      for (int k = 0; k < d; ++k) {
        double dist = (grid_in(i, k) - grid_out(j, k)) / h[k];
        w *= tricube(dist);
      }
      w *= w_in[i];

      w_out[j] += w;
      z_out[j] += z_in[i] * w;
    }
  }

  for(int j = 0; j < n_out; ++j) {
    z_out[j] /= w_out[j];
  }

  return z_out;
}

// Function to calculate deviation from parity from cd matrix
NumericVector distParity(NumericMatrix mat,
			 NumericVector popvec)
{

  /* Inputs to function:
     mat: Matrix of congressional district assignments

     popvec: vector of geographic unit populations

   */

  // Get the unique cd labels
  NumericVector labs = unique(mat(_,0));

  // Calculate parity
  double parity = (double)sum(popvec) / labs.size();

  // Container vector of plan deviations
  NumericVector plandevs(mat.ncol());

  // Loop through plans
  for(int i = 0; i < mat.ncol(); i++){

    // Get plan
    arma::vec plan = mat(_,i);

    // Loop through assignments
    double maxdev = 0.0;
    for(int j = 0; j < labs.size(); j++){

      arma::uvec assignments = find(plan == labs(j));

      // Loop over precincts in plan
      int distpop = 0;
      for(int k = 0; k < assignments.size(); k++){

	distpop += popvec(assignments(k));

      }

      // Get deviation from parity
      double plandev = std::abs((double)((double)distpop - parity) / parity);
      if(plandev > maxdev){
	maxdev = plandev;
      }

    }

    // Store maxdev in plandevs
    plandevs(i) = maxdev;

  }

  return plandevs;

}

// [[Rcpp::export]]
NumericVector getLengths(NumericMatrix xs, NumericMatrix ys) {
  NumericVector x(xs.ncol()), y(ys.ncol());
  NumericVector Lengths(xs.nrow());
  for(int i = 0; i < xs.nrow(); i++){
    x = xs(i,_);
    y = ys(i,_);
    Lengths[i] = getLength(x,y);
    }
  return Lengths;
  }

// [[Rcpp::export]]
double sample_gamma_scalar(NumericVector y, List x, NumericMatrix groups,
                    NumericMatrix beta, List gamma, NumericVector delta,
                    NumericMatrix z, double alpha, double sigma2,
                    List bs_beta, NumericVector sigma2_gamma,
                    int q_id, int G_id) {
  double out;
  NumericVector params(2);
  
  double sum_term = 0.0;
  double sum_term2 = 0.0;
  double sum_term2a = 0.0;
  double sum_term2b = 0.0;
  double sum_term3 = 0.0;
  double sum_term4 = 0.0;
  for (int i = 0; i < y.size(); i++) {
    if (groups(i, q_id) == G_id) {
      double term = 0.0;
      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++) {
          for (int p = 0; p < beta.ncol(); p++) {
            term += beta(k, p) * bs_beta_k(j, p) * x_k(i, j);
          }
        }
      }
      sum_term += 1;
      sum_term2 += y[i];
      sum_term2a += alpha;
      sum_term2b += term;
      for (int k = 0; k < z.ncol(); k++) {
        sum_term4 += delta[k] * z(i, k);
      }
      for (int q = 0; q < gamma.size(); q++) {
        if (q != q_id) {
          NumericVector gamma_q = gamma[q];
          int group_set = groups(i, q);
          sum_term3 += gamma_q[group_set];
        }
      }
    }
  }
  
  //Calculate params[1]
  params[1] = (sigma2 * sigma2_gamma[q_id]) / (sigma2_gamma[q_id] * sum_term + sigma2);
  
  //Calculate params[0]
  params[0] = params[1] * (-1.0 / (2.0 * sigma2)) * (-2.0 * sum_term2 + 2.0 * sum_term2a +
                                                     2.0 * sum_term2b + 2.0 * sum_term3 +
                                                     2.0 * sum_term4);
  
  //Calculate loglik
  out = rnorm(1, params[0], sqrt(params[1]))[0];
  return out;
}

//' Fast computation of trace of matrix product
//'
//' @description Fast computation of the trace of the matrix product trace(t(A) %*% B)
//' @param A A matrix with dimensions n*k.
//' @param B A matrix with dimenions n*k.
//' @return The trace of the matrix product
//' @author Claus Ekstrom <[email protected]@rprimer.dk>
//' @examples
//'
//' A <- matrix(1:12, ncol=3)
//' tracemp(A, A)
//'
//' @export
// [[Rcpp::export]]
double tracemp(NumericMatrix A, NumericMatrix B) {

  if ((A.nrow() != B.nrow()) || (A.ncol() != B.ncol()))
    Rcpp::stop("the two matrices must have the same dimensions");

  arma::mat X(A.begin(), A.nrow(), A.ncol(), false);
  arma::mat Y(B.begin(), B.nrow(), B.ncol(), false);
  double res = arma::as_scalar(accu(X % Y));
    
  return res;
}

// [[Rcpp::export]]
NumericMatrix cpp_rdirichlet(
    const int& n,
    const NumericMatrix& alpha
  ) {
  
  if (std::min({alpha.nrow(), alpha.ncol()}) < 1) {
    Rcpp::warning("NAs produced");
    NumericMatrix out(n, alpha.ncol());
    std::fill(out.begin(), out.end(), NA_REAL);
    return out;
  }

  int k = alpha.ncol();
  NumericMatrix x(n, k);
  
  bool throw_warning = false;
  
  if (k < 2)
    Rcpp::stop("number of columns in alpha should be >= 2");
  
  double row_sum, sum_alpha;
  bool wrong_values;

  for (int i = 0; i < n; i++) {
    sum_alpha = 0.0;
    row_sum = 0.0;
    wrong_values = false;

    for (int j = 0; j < k; j++) {
      sum_alpha += GETM(alpha, i, j);
      if (GETM(alpha, i, j) <= 0.0) {
        wrong_values = true;
        break;
      }
      
      x(i, j) = R::rgamma(GETM(alpha, i, j), 1.0);
      row_sum += x(i, j);
    }

    if (ISNAN(sum_alpha) || wrong_values) {
      throw_warning = true;
      for (int j = 0; j < k; j++)
        x(i, j) = NA_REAL;
    } else {
      for (int j = 0; j < k; j++)
        x(i, j) /= row_sum;
    }
  }
  
  if (throw_warning)
    Rcpp::warning("NAs produced");

  return x;
}

//' Apply crossprod and colSums 
//'
//' @description Fast computation of crossprod(colSums(X),Y) 
//' @param X A matrix with dimensions k*n. Hence the result of \code{colSums(X)} has length n.
//' @param Y A matrix with dimenions n*m. Can be a matrix with dimension m*n but then \code{transposeY} should be \code{TRUE}.
//' @param transposeY Logical. If \code{TRUE} transpose Y before matrix multiplication.
//' @return A vector of length m.
//' @author Thomas Alexander Gerds <[email protected]@biostat.ku.dk>
//' @examples
//' x <- matrix(1:8,ncol=2)
//' y <- matrix(1:16,ncol=8)
//' colSumsCrossprod(x,y,0)
//' 
//' x <- matrix(1:8,ncol=2)
//' y <- matrix(1:16,ncol=2)
//' colSumsCrossprod(x,y,1)
//' @export
// [[Rcpp::export]]
NumericMatrix colSumsCrossprod(NumericMatrix X, NumericMatrix Y, bool transposeY){
  arma::mat A(X.begin(), X.nrow(), X.ncol(), false);
  arma::mat B(Y.begin(), Y.nrow(), Y.ncol(), false);
  arma::rowvec result;
  // result of colSums(A) has to be a matrix
  // with one row and as many columns as B has rows
  if (transposeY)
    result = arma::sum(A,0)*B.t();
  else
    result = arma::sum(A,0)*B;
  return wrap(result); 
}

// [[Rcpp::export]]
NumericVector smooth_nd_1(const NumericMatrix& grid_in, 
                          const NumericVector& z_in, 
                          const NumericVector& w_in_,
                          const NumericMatrix& grid_out, 
                          const int var, const double h,
                          const std::string type = "mean") {

  if (var < 0) stop("var < 0");
  if (var >= grid_in.ncol()) stop("var too large");
  if (h <= 0) stop("h <= 0");
  if (grid_in.ncol() != grid_out.ncol()) stop("Incompatible grid sizes");

  int n_in = grid_in.nrow(), n_out = grid_out.nrow(), d = grid_in.ncol();
  NumericVector w_in = (w_in_.size() > 0) ? w_in_ : 
    rep(NumericVector::create(1), n_in);
  NumericVector z_out(n_out), w_out(n_out);

  // Will be much more efficient to slice up by input dimension:
  // and most efficient way of doing that will be to bin with / bw
  // My data structure: sparse grids
  // 
  // And once we're smoothing in one direction, with guaranteed e2venly spaced
  // grid can avoid many kernel evaluations and can also compute more
  // efficient running sum

  for(int j = 0; j < n_out; ++j) {
    boost::shared_ptr<Summary2d> summary = createSummary(type);
    for (int i = 0; i < n_in; ++i) {
      // Check that all variables (apart from var) are equal
      bool equiv = true;
      for (int k = 0; k < d; ++k) {
        if (k == var) continue;

        double in = grid_in(i, k), out = grid_out(j, k);
        if (in != out && !both_na(in, out)) {
          equiv = false;
          break;
        }
      };
      if (!equiv) continue;

      double in = grid_in(i, var), out = grid_out(j, var);
      double dist = both_na(in, out) ? 0 : in - out;
      double w = tricube(dist / h) * w_in[i];
      if (w == 0) continue;

      summary->push(dist, z_in[i], w);
    }
    z_out[j] = summary->compute();
  }

  return z_out;
}

//---------------------------------------------------------------------
//
//---------------------------------------------------------------------
// [[Rcpp::export]]
void numerator_trick(NumericMatrix A, NumericMatrix B) {
    
  arma::mat aA(A.begin(), A.nrow(), A.ncol(), false);
  arma::mat aB(B.begin(), B.nrow(), B.ncol(), false);
  
  arma::mat numerator_mat = arma::conv2(aA,arma::fliplr(arma::flipud((aB))));
  arma::cx_mat res = arma::ifft2(arma::fft2(aA) % arma::fft2(arma::fliplr(arma::flipud(aB))));
  res.print();
  numerator_mat.print();
    
  //return numerator_mat;
  
}