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

// Execute a single coordinate ascent update to maximize the variational 
// lower bound for Bayesian variable selection in linear regression.
//
// [[Rcpp::export]]
void varbvsnormupdate_rcpp (const NumericMatrix& X, double sigma, double sa,
			    const NumericVector& logodds,
			    const NumericVector& xy, const NumericVector& d,
			    NumericVector& alpha, NumericVector& mu,
			    NumericVector& Xr, const IntegerVector& i) {

  // Cycle through the coordinate ascent updates.
  for (int iter = 0; iter < i.size(); iter++) {

    // j is the coordinate we update in this iteration.
    int j = i[iter];

    // Compute the variational estimate of the posterior variance.
    double s = sa*sigma/(sa*d[j] + 1);

    // Update the variational estimate of the posterior mean.
    double r = alpha[j] * mu[j];
    mu[j] = s/sigma*(xy[j] + d[j]*r - sum(X.column(j) * Xr));

    // Update the variational estimate of the posterior inclusion probability.
    alpha[j] = sigmoid_rcpp(logodds[j] + (log(s/(sa*sigma)) +
					  mu[j]*mu[j]/s)/2);

    // Update Xr = X*r.
    Xr = Xr + (alpha[j] * mu[j] - r) * X.column(j);
  }
}

// [[Rcpp::export]]
List runit_Row_Column_sugar( NumericMatrix x){
    NumericVector r0 = x.row(0) ;
    NumericVector c0 = x.column(0) ;
    return List::create(
        r0,
        c0,
        x.row(1),
        x.column(1),
        x.row(1) + x.column(1)
        ) ;
}

// [[Rcpp::export]]
List softmaxUnitCpp(NumericMatrix input)
{
  int nrows = input.nrow();
  int ncols = input.ncol();
  NumericMatrix activations = transpose(input);
  NumericMatrix derivatives(Dimension(ncols, nrows));
  
  for (int i = 0; i < nrows; i++)
  {
    activations(_, i) = exp(activations.column(i) - max(activations.column(i)));
    activations(_, i) = activations.column(i) / sum(activations.column(i));
    derivatives(_, i) = activations.column(i) * (1 - activations.column(i));
  }
  
  return List::create(transpose(activations), transpose(derivatives));
}

// for multiple markers, but only for one traits input, for multiple traits using version 2 or
// cc <- lapply(all[c(2,3)], function(x) mylm(x,marker3))
// [[Rcpp::export]]
List mylm(NumericVector y, NumericMatrix x){
	//int n = x.size(); //for a list
    int n = x.ncol();
    List out(n);
	for(int i = 0; i < n; ++i) {
        //NumericMatrix m1 = add1(x[i], y);
        NumericMatrix m1 = add1(x(_,i), y);
        NumericMatrix mm = wrap(cleanmat(m1)); //wrap for arma::mat to Rcpp::NumericMatrix
        int m = mm.nrow();
        // seems x.column(i) and x(_,i) can both extract columns
        NumericVector y2 = mm.column(0);
        NumericMatrix x2 (m, 2);
        x2(_,0) = mm.column(1); x2(_,1) = mm.column(2);
        out[i] = fastLm1(y2, x2);
        //out[i] = x2;
	}
	return out;
}

// [[Rcpp::export]]
NumericMatrix runit_NumericMatrix_cumsum( NumericMatrix input ){
    int nr = input.nrow(), nc = input.ncol() ;
    NumericMatrix output(nr, nc) ;
    NumericVector tmp( nr );
    for( int i=0; i<nc; i++){
        tmp = tmp + input.column(i) ;
        NumericMatrix::Column target( output, i ) ;
        std::copy( tmp.begin(), tmp.end(), target.begin() ) ;
    }
    return output ;
}

// [[Rcpp::export]]
List sigmoidUnitCpp(NumericMatrix input)
{
  int nrows = input.nrow();
  int ncols = input.ncol();
  NumericMatrix activations(Dimension(nrows, ncols));
  NumericMatrix derivatives(Dimension(nrows, ncols));
  
  for (int i = 0; i < ncols; i++)
  {
    activations(_, i) = 1/(1 + exp(-input.column(i)));
    derivatives(_, i) = activations.column(i) * (1 - activations.column(i));
  }
  
  return List::create(activations, derivatives);
}

NumericMatrix col_erase (NumericMatrix& x, IntegerVector& colID) {
  colID = colID.sort();

  NumericMatrix x2(Dimension(x.nrow(), x.ncol()- colID.size()));
  int iter = 0;
  int del = 1;
  for (int i = 0; i < x.ncol(); i++) {
    if (i != colID[del - 1]) {
      x2.column(iter) = x.column(i);
      iter++;
    } else {
      del++;
    }
  }
  return x2;
}

// [[Rcpp::export]]
double runit_NumericMatrix_column( NumericMatrix m ){
    NumericMatrix::Column col = m.column(0) ;
    return std::accumulate( col.begin(), col.end(), 0.0 ) ;
}

// [[Rcpp::export]]
List cxxMixEM(NumericMatrix matrix_lik, NumericVector prior, NumericVector pi_init, double tol=0.0001, int maxiter=5000){//note: no default pi_init=NULL  
  int n=matrix_lik.nrow(), k=matrix_lik.ncol(), j=0;
  bool converged=NA_LOGICAL;
  NumericVector pi(k);

  if(Rf_isNull(pi_init))
    std::fill(pi.begin(), pi.end(), 1./(double)k);
  else{
    pi=clone(pi_init);
    for (int i=0;i<k;i++)//set any estimates that are very small to be very small    
      pi[i]=std::max(1e-5, pi[i]); 
    pi=pi/sum(pi); //normalize pi
  }
  NumericMatrix m(n,k);
  NumericVector m_rowsum(n);
  NumericMatrix classprob(m);
  std::vector<double> loglik, lpriordens, penloglik;
  loglik.reserve(maxiter);
  lpriordens.reserve(maxiter); 
  penloglik.reserve(maxiter);
  
  for (int i=0;i<k;i++){
    m.column(i)=pi[i]*matrix_lik.column(i);
    m_rowsum=m_rowsum+m.column(i);
  }
  for (int i=0;i<k;i++)//can vectorize this with sugar?
    classprob.column(i)=classprob.column(i)/m_rowsum;
  loglik.push_back(sum(log(m_rowsum)));
  lpriordens.push_back(sum((prior-1.)*log(pi)));
  penloglik.push_back(sum(log(m_rowsum)) + sum((prior-1.)*log(pi)));
  
  for(j=1;j<maxiter;j++){
    //update pi
    //to do: can vectorize this with sugar?
    for (int i=0;i<k;i++)//set any estimates that are less than zero, which can happen with prior<1, to 0
      pi[i]=std::max(1e-5, sum(classprob.column(i))+prior[i]-1.);
    pi=pi/sum(pi); //normalize pi
    
    //Now re-estimate pi
    std::fill(m_rowsum.begin(), m_rowsum.end(), 0);
    for (int i=0;i<k;i++){
      m.column(i)=pi[i]*matrix_lik.column(i);
      m_rowsum=m_rowsum+m.column(i);
    }
    for (int i=0;i<k;i++)
      classprob.column(i)=classprob.column(i)/m_rowsum;
    loglik.push_back(sum(log(m_rowsum)));
    lpriordens.push_back(sum((prior-1.)*log(pi))); 
    penloglik.push_back(sum(log(m_rowsum)) + sum((prior-1.)*log(pi)));
 
    converged=(bool) (std::abs(loglik[j]+lpriordens[j]-loglik[j-1]-lpriordens[j-1])<tol);
    if(converged)
      break;
  }
  if (j==maxiter) 
    j-=1;
  
  return(List::create(Named("pihat")=pi,
                      Named("B")=loglik,
                      Named("penloglik")=penloglik,
                      Named("niter")=wrap(j+1),
                      Named("converged")=wrap(converged)));
}

// [[Rcpp::export]]
NumericMatrix viterbi_ths(NumericVector &theta, NumericMatrix &data,
	                  NumericVector &integrControl) {
  // theta lambda0, lambda1, lambda2, sigma, p
  // data diff of t and x
  int n = data.nrow(); int dim = data.ncol() - 1;
  double lambda0 = theta[0], lambda1 = theta[1], lambda2 = theta[2];
  double p = theta[4];
  if (lambda1 < lambda2) return NA_REAL;
  double ps0 = 1. / lambda0 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps1 = p / lambda1 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps2 = (1 - p) / lambda2 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  NumericVector tt = data.column(0);
  NumericMatrix x = data(Range(0, n - 1), Range(1, dim));

  // result matrix: three cols stand for Viterbi prob of state 0,1,2
  //                at current time points. For numerical reason,
  //                the log-prob is returned.
  NumericMatrix result(n + 1, 3);
  result(0, 0) = log(ps0); result(0, 1) = log(ps1); result(0, 2) = log(ps2);
  NumericVector cartV = result.row(0);
  NumericVector cartW(3);

  // calculate all h functions
  NumericVector
    hresult00 = ths_h00(x, tt, theta, integrControl),
    hresult01 = ths_h01(x, tt, theta, integrControl),
    hresult02 = ths_h02(x, tt, theta, integrControl),
    hresult10 = ths_h10(x, tt, theta, integrControl),
    hresult11 = ths_h11(x, tt, theta, integrControl),
    hresult12 = ths_h12(x, tt, theta, integrControl),
    hresult20 = ths_h20(x, tt, theta, integrControl),
    hresult21 = ths_h21(x, tt, theta, integrControl),
    hresult22 = ths_h22(x, tt, theta, integrControl);
  
  for (int i = 0; i < n; i++) {
    NumericVector crow = x.row(i);
    if (is_true(all(crow == 0.))) {
      hresult00[i] = 0.;
      hresult01[i] = 0.;
      hresult02[i] = 0.;
      hresult10[i] = 0.;
      hresult11[i] = exp(-lambda1 * tt[i]);
      hresult12[i] = 0.;
      hresult20[i] = 0.;
      hresult21[i] = 0.;
      hresult22[i] = exp(-lambda2 * tt[i]);
    }
  }

  // calculate Viterbi path
  for (int i = 1; i <= n; i++) {
    cartW[0] = cartV[0] + log(hresult00[i - 1]);
    cartW[1] = cartV[1] + log(hresult10[i - 1]);
    cartW[2] = cartV[2] + log(hresult20[i - 1]);
    result(i, 0) = max(cartW);

    cartW[0] = cartV[0] + log(hresult01[i - 1]);
    cartW[1] = cartV[1] + log(hresult11[i - 1]);
    cartW[2] = cartV[2] + log(hresult21[i - 1]);
    result(i, 1) = max(cartW);

    cartW[0] = cartV[0] + log(hresult02[i - 1]);
    cartW[1] = cartV[1] + log(hresult12[i - 1]);
    cartW[2] = cartV[2] + log(hresult22[i - 1]);
    result(i, 2) = max(cartW);

    cartV = result.row(i);
  }

  return result;
}

// [[Rcpp::export]]
NumericMatrix partial_viterbi_ths(NumericVector &theta, NumericMatrix &data,
				  NumericVector &integrControl,
				  int &startpoint, int &pathlength){
  // theta lambda0, lambda1, lambda2, sigma, p
  // data diff of t and x
  // startpoint the start time point, note that
  //            the first time point in data is t0
  // pathlength the length of partial viterbi path
  int n = data.nrow(); int dim = data.ncol() - 1;
  double lambda0 = theta[0], lambda1 = theta[1], lambda2 = theta[2];
  double p = theta[4];
  double ps0 = 1. / lambda0 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps1 = p / lambda1 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps2 = (1 - p) / lambda2 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  NumericVector tt = data.column(0);
  NumericMatrix x = data(Range(0, n - 1), Range(1, dim));

  // bf_result matrix: frist three col for forward
  //                   last  three col for backward
  NumericMatrix bf_result(n + 1, 6);
  bf_result(0, 0) = ps0; bf_result(0, 1) = ps1; bf_result(0, 2) = ps2;
  bf_result(n, 3) =   1; bf_result(n, 4) =   1; bf_result(n, 5) =   1;
  NumericVector dx(n);
  
  // result matrix: three cols stand for Viterbi prob of state 0,1,2
  //                at current time points. For numerical reason,
  //                the log-prob is returned.
  NumericMatrix result(pathlength, 3);
  NumericVector cartV(3);
  NumericVector cartW(3);
  
  

  // calculate all h functions
  NumericVector
    hresult00 = ths_h00(x, tt, theta, integrControl),
    hresult01 = ths_h01(x, tt, theta, integrControl),
    hresult02 = ths_h02(x, tt, theta, integrControl),
    hresult10 = ths_h10(x, tt, theta, integrControl),
    hresult11 = ths_h11(x, tt, theta, integrControl),
    hresult12 = ths_h12(x, tt, theta, integrControl),
    hresult20 = ths_h20(x, tt, theta, integrControl),
    hresult21 = ths_h21(x, tt, theta, integrControl),
    hresult22 = ths_h22(x, tt, theta, integrControl);
  
  for (int i = 0; i < n; i++) {
    NumericVector crow = x.row(i);
    if (is_true(all(crow == 0.))) {
      hresult00[i] = 0.;
      hresult01[i] = 0.;
      hresult02[i] = 0.;
      hresult10[i] = 0.;
      hresult11[i] = exp(-lambda1 * tt[i]);
      hresult12[i] = 0.;
      hresult20[i] = 0.;
      hresult21[i] = 0.;
      hresult22[i] = exp(-lambda2 * tt[i]);
    }
  }

  // forward algorithm
  for (int i = 0; i < n; i++) {
    double sumf0 = bf_result(i, 0) * hresult00[i] + bf_result(i, 1) * hresult10[i] + bf_result(i, 2) * hresult20[i];
    double sumf1 = bf_result(i, 0) * hresult01[i] + bf_result(i, 1) * hresult11[i] + bf_result(i, 2) * hresult21[i];
    double sumf2 = bf_result(i, 0) * hresult02[i] + bf_result(i, 1) * hresult12[i] + bf_result(i, 2) * hresult22[i];
    dx[i] = sumf0 + sumf1 + sumf2;
    bf_result(i + 1, 0) = sumf0 / dx[i];
    bf_result(i + 1, 1) = sumf1 / dx[i];
    bf_result(i + 1, 2) = sumf2 / dx[i];
  }

  // backward algorithm
  for (int i = 0; i < n; i++) {
    double sumb0 = bf_result(n-i, 3) * hresult00[n-i-1] + bf_result(n-i, 4) * hresult01[n-i-1] + bf_result(n-i, 5) * hresult02[n-i-1];
    double sumb1 = bf_result(n-i, 3) * hresult10[n-i-1] + bf_result(n-i, 4) * hresult11[n-i-1] + bf_result(n-i, 5) * hresult12[n-i-1];
    double sumb2 = bf_result(n-i, 3) * hresult20[n-i-1] + bf_result(n-i, 4) * hresult21[n-i-1] + bf_result(n-i, 5) * hresult22[n-i-1];
    bf_result(n-i-1, 3) = sumb0 / dx[n-i-1];
    bf_result(n-i-1, 4) = sumb1 / dx[n-i-1];
    bf_result(n-i-1, 5) = sumb2 / dx[n-i-1];
  }

  // prepare for viterbi path
  result(0, 0) = log(bf_result(startpoint, 0));
  result(0, 1) = log(bf_result(startpoint, 1));
  result(0, 2) = log(bf_result(startpoint, 2));
  cartV = result.row(0);
  int ite_stop = startpoint + pathlength - 2;

  // viterbi algorithm
  for (int i = startpoint; i < ite_stop; i++) {
    cartW[0] = cartV[0] + log(hresult00[i]);
    cartW[1] = cartV[1] + log(hresult10[i]);
    cartW[2] = cartV[2] + log(hresult20[i]);
    result(i - startpoint + 1, 0) = max(cartW);

    cartW[0] = cartV[0] + log(hresult01[i]);
    cartW[1] = cartV[1] + log(hresult11[i]);
    cartW[2] = cartV[2] + log(hresult21[i]);
    result(i - startpoint + 1, 1) = max(cartW);
//.........这里部分代码省略.........

// [[Rcpp::export]]
NumericMatrix fwd_bwd_ths(NumericVector &theta, NumericMatrix &data,
	                  NumericVector &integrControl) {
  // theta lambda0, lambda1, lambda2, sigma, p
  // data diff of t and x
  int n = data.nrow(); int dim = data.ncol() - 1;
  double lambda0 = theta[0], lambda1 = theta[1], lambda2 = theta[2];
  double p = theta[4];
  if (lambda1 < lambda2) return NA_REAL;
  double ps0 = 1. / lambda0 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps1 = p / lambda1 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  double ps2 = (1 - p) / lambda2 / (1. / lambda0 + p / lambda1 + (1 - p) / lambda2);
  NumericVector tt = data.column(0);
  NumericMatrix x = data(Range(0, n - 1), Range(1, dim));

  // result matrix: frist three col for forward
  //                last  three col for backward
  NumericMatrix result(n + 1, 6);
  result(0, 0) = ps0; result(0, 1) = ps1; result(0, 2) = ps2;
  result(n, 3) =   1; result(n, 4) =   1; result(n, 5) =   1;
  NumericVector dx(n);

  // calculate all h functions
  NumericVector
    hresult00 = ths_h00(x, tt, theta, integrControl),
    hresult01 = ths_h01(x, tt, theta, integrControl),
    hresult02 = ths_h02(x, tt, theta, integrControl),
    hresult10 = ths_h10(x, tt, theta, integrControl),
    hresult11 = ths_h11(x, tt, theta, integrControl),
    hresult12 = ths_h12(x, tt, theta, integrControl),
    hresult20 = ths_h20(x, tt, theta, integrControl),
    hresult21 = ths_h21(x, tt, theta, integrControl),
    hresult22 = ths_h22(x, tt, theta, integrControl);
  
  for (int i = 0; i < n; i++) {
    NumericVector crow = x.row(i);
    if (is_true(all(crow == 0.))) {
      hresult00[i] = 0.;
      hresult01[i] = 0.;
      hresult02[i] = 0.;
      hresult10[i] = 0.;
      hresult11[i] = exp(-lambda1 * tt[i]);
      hresult12[i] = 0.;
      hresult20[i] = 0.;
      hresult21[i] = 0.;
      hresult22[i] = exp(-lambda2 * tt[i]);
    }
  }

  // forward algorithm
  for (int i = 0; i < n; i++) {
    double sumf0 = result(i, 0) * hresult00[i] + result(i, 1) * hresult10[i] + result(i, 2) * hresult20[i];
    double sumf1 = result(i, 0) * hresult01[i] + result(i, 1) * hresult11[i] + result(i, 2) * hresult21[i];
    double sumf2 = result(i, 0) * hresult02[i] + result(i, 1) * hresult12[i] + result(i, 2) * hresult22[i];
    dx[i] = sumf0 + sumf1 + sumf2;
    result(i + 1, 0) = sumf0 / dx[i];
    result(i + 1, 1) = sumf1 / dx[i];
    result(i + 1, 2) = sumf2 / dx[i];
  }

  //backward algorithm
  for (int i = 0; i < n; i++) {
    double sumb0 = result(n-i, 3) * hresult00[n-i-1] + result(n-i, 4) * hresult01[n-i-1] + result(n-i, 5) * hresult02[n-i-1];
    double sumb1 = result(n-i, 3) * hresult10[n-i-1] + result(n-i, 4) * hresult11[n-i-1] + result(n-i, 5) * hresult12[n-i-1];
    double sumb2 = result(n-i, 3) * hresult20[n-i-1] + result(n-i, 4) * hresult21[n-i-1] + result(n-i, 5) * hresult22[n-i-1];
    result(n-i-1, 3) = sumb0 / dx[n-i-1];
    result(n-i-1, 4) = sumb1 / dx[n-i-1];
    result(n-i-1, 5) = sumb2 / dx[n-i-1];
  }

  return result;
}

//' Viterbi algorithm
//'
//' Standard viterbi algorithm in the log space
//' @param initP matrix of initial probabilities: each column corresponds to a sequence
//' @param trans transition matrix (rows are previous state, columns are next state)
//' @param lliks matrix with emission probabilities for each datapoint and each state.
//' Columns are datapoints and rows are states.
//' @param seqlens length of each subsequence of datapoints (set this to ncol(lliks)
//' if there is only one sequence).
//' @return a list with the following arguments:
//'    \item{vpath}{viterbi path}
//'    \item{vllik}{log-likelihood of the viterbi path}
//' @export
// [[Rcpp::export]]
List viterbi(NumericMatrix initP, NumericMatrix trans, NumericMatrix lliks, NumericVector seqlens){
    int nmod = initP.nrow();
    double totlen = Rcpp::sum(seqlens);
    if (nmod != trans.nrow() || nmod != trans.ncol() || nmod != lliks.nrow()) Rcpp::stop("Unable to figure out the number of models");
    if (((double) lliks.ncol()) != totlen) Rcpp::stop("Sequence lengths don't match with the provided matrix");
    
    int ncol = lliks.ncol();
    IntegerVector vpath(ncol);
    IntegerMatrix backtrack(nmod, max(seqlens));
    std::vector<long double> scores(nmod);
    std::vector<long double> new_scores(nmod);
    
    /* log-transform the transition probabilities */
    NumericMatrix ltrans(nmod,nmod);
    for (diter curr = ltrans.begin(), currt = trans.begin(); curr < ltrans.end(); ++curr, ++currt){
        *curr = log(*currt);
    }
    
    /* Viterbi independently on each chunk */
    double tot_maxscore = 0;
    for (int o = 0, chunk_start = 0; o < seqlens.length(); chunk_start += seqlens[o], ++o){
        int chunk_end = chunk_start + seqlens[o];
        /* dynamic programming */
        {
            MatrixColumn<REALSXP> llikcol = lliks.column(chunk_start);
            MatrixColumn<REALSXP> curr_initP = initP.column(o);
            for (int t = 0; t < nmod; ++t){
                scores[t] = llikcol[t] + log(curr_initP[t]);
            }
        }
        for (int i = chunk_start + 1; i < chunk_end; ++i){
            
            MatrixColumn<REALSXP> llikcol = lliks.column(i);
            MatrixColumn<INTSXP> backtrackcol = backtrack.column(i-chunk_start);
            
            for (int t = 0; t < nmod; ++t){
                int maxs = 0;
                long double maxscore = scores[0] + ltrans(0, t);
                for (int s = 1; s < nmod; ++s){
                    long double currscore = scores[s] + ltrans(s,t);
                    if (currscore > maxscore){
                        maxscore = currscore;
                        maxs = s;
                    }
                }
                backtrackcol[t] = maxs;
                new_scores[t] = llikcol[t] + maxscore;
            }
            
            memcpy(scores.data(), new_scores.data(), sizeof(long double)*nmod);
        }
        
        /* backtracking */
        int maxp = 0;
        double maxscore = scores[0];
        for (int p = 1; p < nmod; ++p){
            if (scores[p] > maxscore){
                maxscore = scores[p];
                maxp = p;
            }
        }
        tot_maxscore += maxscore;
        vpath[chunk_end - 1] = maxp + 1;
        for (int i = chunk_end - 2; i >= chunk_start; --i){
            maxp = backtrack(maxp, i - chunk_start + 1);
            vpath[i] = maxp + 1; //in R indices are 1-based
        }
    }
    return List::create(_("vpath")=vpath, _("vllik")=tot_maxscore);
}