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

//---------------------------------------------------------------------
//
//---------------------------------------------------------------------
// [[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;
  
}

// **********************************************************//
//              Calculate mu matrix for document d           //
// **********************************************************//
// [[Rcpp::export]]
NumericVector mu_vec(arma::vec p_d, NumericMatrix xi) {
    int nIP = xi.ncol();
    NumericVector muvec(xi.nrow());
    for (int IP = 0; IP < nIP; IP++) {
        double pdIP = p_d[IP];
        for (int i = 0; i < xi.nrow(); i++) {
        	if (pdIP > 0) {
           		muvec[i] += p_d[IP]*xi(i,IP);
        	}
        }
    }
    return muvec;
}

// [[Rcpp::export]]
NumericMatrix rcpp_parallel_js_distance(NumericMatrix mat) {
  
   // allocate the matrix we will return
   NumericMatrix rmat(mat.nrow(), mat.nrow());

   // create the worker
   JsDistance jsDistance(mat, rmat);
     
   // call it with parallelFor
   parallelFor(0, mat.nrow(), jsDistance);

   return rmat;
}

// [[Rcpp::export]]
NumericMatrix residLm(NumericMatrix Yr, NumericMatrix Xr) {
   
   int nX = Xr.nrow(), nY = Yr.nrow();
   arma::mat Y(Yr.begin(), nY, nX, false); // reuses memory and avoids extra copy
   arma::mat X(Xr.begin(), nX, 2, false);
   arma::mat resid(nX,nY);
   arma::colvec y;
   for(int i = 0; i < nY; i++){     
     y = Y.row(i).t();
     resid.col(i) = y - X*arma::solve(X, y);    
  }
     
   return wrap(resid.t());
}

// [[Rcpp::export]]
NumericMatrix parallelMatrixSqrt(NumericMatrix x) {
  
  // allocate the output matrix
  NumericMatrix output(x.nrow(), x.ncol());
  
  // SquareRoot functor (pass input and output matrixes)
  SquareRoot squareRoot(x, output);
  
  // call parallelFor to do the work
  parallelFor(0, x.nrow() * x.ncol(), squareRoot);
  
  // return the output matrix
  return output;
}

// [[Rcpp::export]]
List gemmFitRcppI(int n, NumericMatrix betas, NumericMatrix data, int p, NumericVector kCor, CharacterVector correction, bool isTauB) {

  NumericVector fit;
  NumericVector fitR(betas.nrow()), fitTauA(betas.nrow()), fitTauB(betas.nrow()), fitBIC(betas.nrow()), fitBICr(betas.nrow()), fitAIC(betas.nrow()), fitAICr(betas.nrow());
  NumericVector fitTauPairs(betas.nrow()), fitTauTies1(betas.nrow()), fitTauTies2(betas.nrow());
  NumericVector fitTauTiesBoth(betas.nrow()), fitTauDis(betas.nrow()), fitTauCon(betas.nrow());

  NumericMatrix data2(clone(data));
  
  bool corType; 
  
  /*std::string s_targetname = as<std::string>(targetname)
  */
  
  (as<std::string>(correction)=="knp") ? (corType=true) : (corType=false);


for (int i=0; i < betas.nrow(); i++) {
  fit = gemmFitRcpp(n, betas(i,_), data2, p, kCor(i), corType, isTauB);

    fitBIC(i) = fit[1];
    fitBICr(i) = fit[2];
    fitR(i) = fit[0];
    fitAIC(i) = fit[4];
    fitAICr(i) = fit[5];
    fitTauA(i) = fit[6];
    fitTauB(i) = fit[7];
    fitTauPairs(i) = fit[8];
    fitTauTies1(i) = fit[9];
    fitTauTies2(i) = fit[10];
    fitTauTiesBoth(i) = fit[11];
    fitTauDis(i) = fit[12];
    fitTauCon(i) = fit[13];
  }


  return Rcpp::List::create(Rcpp::Named("r") = fitR,
                            Rcpp::Named("bic") = fitBIC,
                            Rcpp::Named("tau.a") = fitTauA,
                            Rcpp::Named("tau.b") = fitTauB,
                            Rcpp::Named("tau.n.pairs") = fitTauPairs,
                            Rcpp::Named("tau.n.ties.1") = fitTauTies1,
                            Rcpp::Named("tau.n.ties.2") = fitTauTies2,
                            Rcpp::Named("tau.n.ties.both") = fitTauTiesBoth,
                            Rcpp::Named("tau.n.dis") = fitTauDis,
                            Rcpp::Named("tau.n.con") = fitTauCon,
                            Rcpp::Named("bic.r") = fitBICr,
                            Rcpp::Named("aic") = fitAIC,
                            Rcpp::Named("aic.r") = fitAICr);
}

// Sends employed Bees
void SendEmployedBees(
    double & GlobalMin,
    NumericVector & GlobalParams, 
    NumericVector & fitness,
    NumericVector & f,
    IntegerVector & trial,
    Function & fun,
    NumericMatrix & Foods,
    const NumericVector & lb,
    const NumericVector & ub
  ) {
  
  int param2change, neighbour;
  double ObjValSol, FitnessSol;
  NumericVector solution(Foods.ncol());
  for (int i=0;i<Foods.nrow();i++) {
    // Random parameter to change
    param2change = (int)(unif_rand()*Foods.ncol());
    
    // Random neighbour to select
    neighbour    = i;
    while (neighbour == i) 
      neighbour = (int)(unif_rand()*Foods.nrow());
    
    // Suggesting new solution
    solution = Foods(i,_);
    solution[param2change] = Foods(i,param2change) + 
      (Foods(i, param2change) - Foods(neighbour, param2change))*(unif_rand()-.5)*2;
    
    // Truncating
    if (solution[param2change] < lb.at(param2change))
      solution[param2change] = lb.at(param2change);
    if (solution[param2change] > ub.at(param2change))
      solution[param2change] = ub.at(param2change);
    
    // Comparing current solution with new one
    ObjValSol  = as<double>(fun(solution));
    FitnessSol = CalculateFitness(ObjValSol);
    if (FitnessSol > fitness[i]) {
      Foods(i,_) = solution;
      fitness[i] = FitnessSol;
      f[i]       = ObjValSol;
      trial[i]   = 0;
    } else {
      trial[i]+=1;
    }
  } 
  
  return;
}

//' Apply crossprod and rowSums
//'
//' @description Fast computation of crossprod(rowSums(X),Y)
//' @param X A matrix with dimensions n*k. Hence the result of \code{rowSums(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:10,nrow=5)
//' y <- matrix(1:20,ncol=4)
//' rowSumsCrossprod(x,y,0)
//'
//' x <- matrix(1:10,nrow=5)
//' y <- matrix(1:20,ncol=5)
//' rowSumsCrossprod(x,y,1)
//' @export
// [[Rcpp::export]]
NumericMatrix rowSumsCrossprod(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
  // since sum(A,1) is a matrix with one column
  // we transpose before multiplication
  if (transposeY)
    result = arma::sum(A,1).t()*B.t();
  else
    result = arma::sum(A,1).t()*B;
  return wrap(result); 
}

void CalculateProbabilities(
    NumericMatrix & Foods,
    NumericVector & fitness,
    NumericVector & prob
  ) {
  double maxfit = fitness[0];
  for (int i=0;i<Foods.nrow();i++) 
    if (fitness[i] > maxfit) maxfit = fitness[i];
  
  for (int i=0;i<Foods.nrow();i++) 
    prob[i] = (0.9*((fitness[i] + 1e-40)/(maxfit + 1e-40))) + 0.1;
  
  return;
  
}

// **********************************************************//
//      Multiply list of matrix X and list of vector B       //
// **********************************************************//
// [[Rcpp::export]]
List MultiplyXB(List X, NumericMatrix B){
	List XB(B.nrow());
	for (int IP = 0; IP < B.nrow(); IP++) {
		arma::mat XB_IP(X.size(), X.size());
		arma::vec B_IP = B(IP, _);
		for (unsigned int n = 0; n < X.size(); n++) {
			List X_n = X[n];
			arma::mat X_n_IP = X_n[IP];
			arma::vec rows = X_n_IP*B_IP;
			XB_IP.row(n) = rows.t();
		}
		XB[IP] = XB_IP;
	}	
  return XB;
}

// [[Rcpp::export]]
NumericVector hazr(double t, double m, double bm, double bt, NumericMatrix dat){
  // MAKE SURE dat IS SORTED BY TIME OR THIS FUNCTION WILL NOT WORK!!!
  // bm should be on marker scale
  // bt should be on time scale
  // ---------------------------------------------------------------------------
  NumericVector out(2);
  
  // Collect subset of observations in marker trip in mstrip
  // ---------------------------------------------------------------------------
  int j = 0; // number of values in marker strip
  NumericMatrix mstrip(dat.nrow(), 2); // put obsevations in marker strip here
  for (int i = 0; i < dat.nrow(); i++){
    // populate mstrip with observations
    if ((dat(i, 2) < m+bm) && (dat(i,2) > m-bm)){
      mstrip(j, 0) = dat(i, 0);
      mstrip(j, 1) = dat(i, 1);
      j += 1;
    }
  }
  if (j == 0) return(out);
  
  // - Do the local-in-marker hazard and N-A calculations
  // ---------------------------------------------------------------------------
  NumericVector na(j); // create vector for local N-A estimates
  double nac = 0;      // holder for N-A contributions
  na[0] = mstrip(0, 1)/j; // N-A estimate at earliest time in strip
  if (fabs(t - mstrip(0, 0)) <= bt){
    // if first value is close enough to t, include in hazard calculation
    out[0] = na[0];
    out[1] = na[0]; // set N-A estimator
  }
  for (int i = 1; i < j; i++){
    // include rest of qualified observations into hazard calculation
    // calculate remaining N-A estimates
    nac = (mstrip(i, 1)/(j - i)); // censoring indicat. divided by risk set size
    na[i] = na[i-1] + nac;
    if (fabs(t - mstrip(i, 0)) <= bt){
      // if close enough to t, include in hazard calculation
      out[0] += nac;
    }
    if (t >= mstrip(i, 0)){
      // if we have reached t, return this for N-A estimate
      out[1] = na[i];
    }
  }
  out[0] = out[0]/(2*bt);
  return(out);
}

// [[Rcpp::export]]
NumericMatrix imp_neighbour_avg(NumericMatrix x, double k) {
  // input matrix is expected to have >= 3 columns
  NumericMatrix ans = clone(x);
  int nr = ans.nrow(), nc = ans.ncol();

  for(int i = 0; i < nr; i++) {
    // first and last values are set to 0 if NA
    if (R_IsNA(ans(i, 0))) ans(i, 0) = k;
    if (R_IsNA(ans(i, nc-1))) ans(i, nc-1) = k;
    
    for(int j = 1; j < (nc-1); j++) {
      if (R_IsNA(ans(i,j))) {
	// if the next value is NA and all previous values are 0
	// then we set to 0	
	if (R_IsNA(ans(i,j+1))) {
	  NumericVector v = subset(ans.row(i), j);
	  if (allZero(v, k)) ans(i,j) = k;
	} else { // next is not NA, set to mean of neighbours
	  ans(i,j) = (ans(i,j-1) + ans(i,j+1))/2;
	}
      }
    }
  }    
  return(ans);
}

// @export
// [[Rcpp::export]]
NumericMatrix cpp_pMatrix(const NumericMatrix& ExpressionSet,const NumericVector& AgeVector)
{
        
        int nRows = ExpressionSet.nrow();
        int nCols = ExpressionSet.ncol(); 
        NumericMatrix results(nRows,nCols);
        NumericVector DivisorVector(nCols);
        
        for(int stage = 0; stage < nCols; stage++) {
                double divisor = 0;
                for(int gene = 0; gene < nRows; gene++) {
                        divisor  += ExpressionSet(gene, stage);
                }
                
                DivisorVector[stage] = divisor;
        }
        
        for(int stage = 0; stage < nCols; stage++){
                for(int gene = 0; gene < nRows; gene++){
                        results(gene,stage) = (double) AgeVector[gene] * (ExpressionSet(gene, stage)/DivisorVector[stage]);
                }
        }
        return results;
        
}

// [[Rcpp::export]]
void appendRcpp(  List fillVecs, NumericVector newLengths, NumericMatrix retmat, NumericVector retmatLengths ) {
  /* appendRcpp
	Fills numeric matrix
	Loop through rows, filling retmat in with the vectors in list
	then update return matrix size to index the next free
	*/
	
	// Declare vars
	NumericVector fillTmp;
	int sizeOld, sizeNew;
    
	// Pull out dimensions of return matrix to fill
	int nrow = retmat.nrow();
  int ncol = retmat.ncol();
    
	// Check that dimensions match
	if ( nrow != retmatLengths.size() || nrow != fillVecs.size() ) { 
        throw std::range_error("In appendC(): dimension mismatch");
    }
    
	// Traverse ddimensions
	for (int ii=0; ii<ncol; ii++) {
        throw std::range_error("In appendC(): exceeded max cols");
		
	// Iterator for row to fill
        NumericMatrix::Row retRow = retmat(ii, _);
	
  // Fill row of return matrix, starting at first non-zero element
        std::copy( fillTmp.begin(), fillTmp.end(), retRow.begin() + sizeOld );
  
  // Update size of return matrix
        retmatLengths[ii] = sizeNew;
		}
	}

void Chain_Factorial::update_mu(NumericMatrix Y){
  double alpha0 = 1.0/K;
  
  //fit_linear_model(X, y, n, K, mu);
  //fit_Bayesian_linear_model(X, y, n, K, mu, sigma);
  NumericVector w_i, u_i, y_i;
  for(int i=0; i<Y.nrow(); i++){
    w_i = w(i, _);
    u_i = u(i, _);
    y_i = Y(i, _);
    lambdas = calculate_mean_for_all_t(X, w_i, h, K, n);
    double logprob_current = calculate_posterior_prob(y_i, lambdas, w_i, alpha0, K, n);
    
    // propose a new alpha
    // double alpha_proposed = random_walk_log_scale(alpha0, 0.3);
    // propose a new w
    NumericVector w_unnorm_proposed = RWMH(u_i, K, 0.5);
    NumericVector w_proposed = w_unnorm_proposed / sum(w_unnorm_proposed);
    NumericVector lambdas_proposed = calculate_mean_for_all_t(X, w_proposed, h, K, n);
    double logprob_proposed = calculate_posterior_prob(y_i, lambdas_proposed, w_proposed, alpha0, K, n);
    //printf("current %f, proposed %f, accept prob: %f\n\n", logprob_current, logprob_proposed, exp(logprob_proposed - logprob_current));
    
    
    // accept or reject
    if(R::runif(0, 1) < exp(logprob_proposed - logprob_current)){
      u(i, _) = clone(w_unnorm_proposed);
      w(i, _) = clone(w_proposed);
    }
  }
  
}