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

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

NumericVector logLikMixHMM(NumericVector transitionMatrix, NumericVector emissionArray, 
  NumericVector initialProbs, IntegerVector obsArray, NumericMatrix coefs, 
  NumericMatrix X_, IntegerVector numberOfStates) {  
  
  
  IntegerVector eDims = emissionArray.attr("dim"); //m,p,r
  IntegerVector oDims = obsArray.attr("dim"); //k,n,r
  
  int q = coefs.nrow();
  arma::mat coef(coefs.begin(),q,coefs.ncol());
  coef.col(0).zeros();
  arma::mat X(X_.begin(),oDims[0],q);
  arma::mat lweights = exp(X*coef).t();
  if(!lweights.is_finite()){
    return wrap(-std::numeric_limits<double>::max());
    
  }
  lweights.each_row() /= sum(lweights,0);
  arma::colvec init(initialProbs.begin(),eDims[0], true);
  arma::mat transition(transitionMatrix.begin(),eDims[0],eDims[0], true);
  arma::cube emission(emissionArray.begin(), eDims[0], eDims[1], eDims[2], true);
  arma::icube obs(obsArray.begin(), oDims[0], oDims[1], oDims[2], false);
  
  arma::vec alpha(eDims[0]);
  NumericVector ll(oDims[0]);  
  double tmp;
  arma::vec initk(eDims[0]);
  
  for(int k = 0; k < oDims[0]; k++){    
    initk = init % reparma(lweights.col(k),numberOfStates);
    
    for(int i=0; i < eDims[0]; i++){      
      alpha(i) = initk(i);
      for(int r = 0; r < oDims[2]; r++){
        alpha(i) *= emission(i,obs(k,0,r),r);
      }
    }    
    
    tmp = sum(alpha);
    ll(k) = log(tmp);
    alpha /= tmp;
    
    arma::vec alphatmp(eDims[0]);
    
    for(int t = 1; t < oDims[1]; t++){  
      for(int i = 0; i < eDims[0]; i++){
        alphatmp(i) = arma::dot(transition.col(i), alpha);
        for(int r = 0; r < oDims[2]; r++){
          alphatmp(i) *= emission(i,obs(k,t,r),r);
        }
      }
      tmp = sum(alphatmp);
      ll(k) += log(tmp);
      alpha = alphatmp/tmp;
    }
  } 
  
  return ll;
}

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

//' 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]]
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;
  
}

//' 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); 
}

// [[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("multiKernel_noCon_cpp")]]
SEXP multiKernel_noCon(NumericMatrix Yr, NumericMatrix Kr, double tau) {
    int n = Yr.nrow(), p = Yr.ncol();

    arma::mat K(Kr.begin(), n, n, false);       // reuses memory and avoids extra copy
    arma::mat Y(Yr.begin(), n, p, 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 alpha_mat = inv(K + tuning_mat) * Y;

    return Rcpp::List::create(Rcpp::Named("alpha") = alpha_mat);
}

// [[Rcpp::export]]
NumericMatrix CPP_dsm_score_dense(NumericMatrix f, NumericVector f1, NumericVector f2, double N, int am_code, int sparse, int transform_code) {
  if (am_code < 0 || am_code >= am_table_entries)
    stop("internal error -- invalid AM code");
  am_func AM = am_table[am_code]; /* selected association measure */

  int nr = f.nrow(), nc = f.ncol();
  if (am_code != 0 && (nr != f1.size() || nc != f2.size()))
    stop("internal error -- marginal vectors f1 and f2 not conformable with matrix f");

  NumericMatrix scores(nr, nc);

  NumericMatrix::iterator _f = f.begin();
  NumericVector::iterator _f1 = f1.begin();
  NumericVector::iterator _f2 = f2.begin();
  NumericMatrix::iterator _scores = scores.begin();

  int i = 0;
  for (int col = 0; col < nc; col++) {
    for (int row = 0; row < nr; row++) {
      /* frequeny measure (am_code == 0) is a special case, since marginals may not be available (in "reweight" mode) */
      double score = (am_code == 0) ? _f[i] : AM(_f[i], _f1[row], _f2[col], N, sparse);
      _scores[i] = (transform_code) ? transform(score, transform_code) : score;
      i++;
    }
  }

  return scores;
}

 void sapply_master(Function infun, bool realloc, bool cppfun) {
 // apply user-supplied function infun to each vector, update in place
 // realloc=true: allow function input and return dim to differ
 // cppfun=false: infun is regular R function 
 // cppfun=true: infun returns XPtr to C++ function,
 //      infun takes arma::vec&, returns void
     std::size_t icol, thisLen;
     funcPtr fun;
     if (cppfun) {
         SEXP funPtrfun = infun();
         XPtr<funcPtr> xpfun(funPtrfun);
         fun = *xpfun;
     }
     for ( icol = 0; icol < nvec; icol++){
         thisLen = lengths[icol];
         NumericMatrix::iterator colStart = data.begin() + (icol * this->allocLen);
         // copy_aux_mem=false: reuse existing memory, modify data in-place
         // strict=true, arma enforces dim match, no size change allowed
         arma::vec dataVec(colStart, thisLen, false, !realloc);
         if(cppfun) {
             fun(dataVec);
         } else {
             FunFromR(infun, dataVec);
         }
         if ( dataVec.size() != thisLen) {
             // won't make it here unless realloc=true
             // change dimension of data matrix, reset colStart
             grow_assign_sapply(icol, thisLen, dataVec, colStart);
         }
     }
 }

// [[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::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]]
NumericMatrix samplehouseholds_format2(NumericMatrix phi, NumericMatrix w, NumericVector pi,
                               NumericVector d, List lambda,
                               int currrentbatch, int nHouseholds,  int householdsize) {

  int K = w.nrow();
  int L = w.ncol();
  int p = d.length();
  int n_lambdas = lambda.length();
  int *lambda_columns = new int[n_lambdas];
  double **lambdas = new double*[n_lambdas];
  int maxDDtp = phi.nrow();
  int maxdd = maxDDtp / p;

  //int ncol = householdsize * DIM + 1 + householdsize;
  //output data: zero-based
  //column 0: household unique index
  //column 1: member index within the household (pernum:person number?)
  //column 2 to 2+p-1: individual data
  //column 2+p: 2+p+n_lambdas-2: household level data
  //column 2+p+n_lambdas-1: household group indicator
  //column last hh_size: individual group indicator
  int DIM = 2 + p + n_lambdas - 1; //not output the household size
  int ncol = DIM * householdsize + householdsize  + 1;
  NumericMatrix data(nHouseholds, ncol);

  //copy data from list of matrices to C++
  for (int i = 0; i < n_lambdas; i++) {
    NumericMatrix l = lambda[i];
    lambda_columns[i] = l.ncol();
    lambdas[i] = new double[l.length()];
    std::copy(l.begin(), l.end(), lambdas[i]);
  }
  //printf("in samplehouseholds\n");
  NumericVector rand = runif(nHouseholds * ncol); //at most this many
  sampleHouseholds_imp_format2(data.begin(), rand.begin(), lambdas, lambda_columns, w.begin(),
                       phi.begin(), pi.begin(),d.begin(),
                       nHouseholds, householdsize, K, L,maxdd,p, currrentbatch,n_lambdas);

  //clean up
  delete [] lambda_columns;
  for (int i = 0; i < n_lambdas; i++) {
    delete [] lambdas[i];
  }
  delete [] lambdas;
  //printf("done samplehouseholds\n");
  return data;
}

//---------------------------------------------------------------------
//R -> openCV  Like an as. Try to minimize copying?? NOTE: Does not convert 
//NumericMatrix from 64bit. Do that yourself with convertTo on the Mat object
//---------------------------------------------------------------------
void NumericMatrix2openCVMat(NumericMatrix dmat, Mat &odmat) {
  
  //NumericMatrix stored row-wise, but openCV Mat stored column wise?? This seems to work combind with the transpose below.
  Mat otmp(Size(dmat.nrow(), dmat.ncol()), CV_64F, dmat.begin(), Mat::AUTO_STEP);
  
  odmat = otmp.t(); //COULD THERE BE A PROBLEM HERE????? WHY DO WE HAVE TO DO THIS????
  
}