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

// Scan 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 of numeric phenotypes (individuals x phenotypes)
//             (no missing data allowed)
// addcovar  = additive covariates (an intercept, at least)
// intcovar  = interactive covariates (should also be included in addcovar)
// weights   = vector of weights
//
// output    = matrix of residual sums of squares (RSS) (phenotypes x positions)
//
// [[Rcpp::export]]
NumericMatrix scan_binary_onechr_intcovar_weighted_highmem(const NumericVector& genoprobs,
                                                           const NumericMatrix& pheno,
                                                           const NumericMatrix& addcovar,
                                                           const NumericMatrix& intcovar,
                                                           const NumericVector& weights,
                                                           const int maxit=100,
                                                           const double tol=1e-6,
                                                           const double qr_tol=1e-12)
{
    const int n_ind = pheno.rows();
    if(Rf_isNull(genoprobs.attr("dim")))
        throw std::invalid_argument("genoprobs should be a 3d array but has no dim attribute");
    const Dimension d = genoprobs.attr("dim");
    if(d.size() != 3)
        throw std::invalid_argument("genoprobs should be a 3d array");
    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)");

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

    // genotype can
    return scan_binary_onechr_weighted(genoprobs_rev, pheno, addcovar, weights, maxit, tol, qr_tol);
}

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

// LMM scan of a single chromosome with interactive covariates
// this version uses less memory but will be slower
// (since we need to work with each position, one at a time)
// 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_lowmem(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("ncol(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)");

    NumericVector result(n_pos);

    // multiply pheno by eigenvectors and then by weights
    NumericMatrix pheno_rev = matrix_x_matrix(eigenvec, pheno);
    pheno_rev = weighted_matrix(pheno_rev, weights);

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

    for(unsigned int pos=0; pos<n_pos; pos++) {
        Rcpp::checkUserInterrupt();  // check for ^C from user

        // form X matrix
        NumericMatrix X = formX_intcovar(genoprobs, addcovar, intcovar, pos, true);

        // multiply by eigenvectors
        X = matrix_x_matrix(eigenvec, X);

        // multiply by weights
        X = weighted_matrix(X, weights);

        // do regression
        NumericVector rss = calc_rss_linreg(X, pheno_rev, tol);
        result[pos] = -(double)n_ind/2.0*log(rss[0]) + sum_logweights;
    }

    return result;
}

NumericVector WithDoubleFun(NumericMatrix x, double (*CombineFun)(double, double)) {
  NumericVector result(1);
  for (int i = 0; i < x.rows(); i++) {
    result = result + CombineFun(x(i,0), x(i,1));
  }
  return result;
}

NumericVector WithMatrixFun(NumericMatrix x, double (*CombineFun)(NumericMatrix, int)) {
  NumericVector result(1);
  for (int i = 0; i < x.rows(); i++) {
    result = result + CombineFun(x, i);
  }
  return result;
}

// Scan a single chromosome with additive covariates and weights
//
// genoprobs = 3d array of genotype probabilities (individuals x genotypes x positions)
// pheno     = matrix of numeric phenotypes (individuals x phenotypes)
//             (no missing data allowed, values should be in [0,1])
// addcovar  = additive covariates (an intercept, at least)
// weights   = vector of weights
//
// output    = matrix of (weighted) residual sums of squares (RSS) (phenotypes x positions)
//
// [[Rcpp::export]]
NumericMatrix scan_binary_onechr_weighted(const NumericVector& genoprobs,
                                          const NumericMatrix& pheno,
                                          const NumericMatrix& addcovar,
                                          const NumericVector& weights,
                                          const int maxit=100,
                                          const double tol=1e-6,
                                          const double qr_tol=1e-12,
                                          const double eta_max=30.0)
{
    const int n_ind = pheno.rows();
    if(Rf_isNull(genoprobs.attr("dim")))
        throw std::invalid_argument("genoprobs should be a 3d array but has no dim attribute");
    const Dimension d = genoprobs.attr("dim");
    if(d.size() != 3)
        throw std::invalid_argument("genoprobs should be a 3d array");
    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 != weights.size())
        throw std::range_error("nrow(pheno) != length(weights)");
    const int n_pos = d[2];
    const int n_gen = d[1];
    const int n_add = addcovar.cols();
    const int g_size = n_ind * n_gen;
    const int n_phe = pheno.cols();

    NumericMatrix result(n_phe, n_pos);
    NumericMatrix X(n_ind, n_gen+n_add);

    if(n_add > 0) // paste in covariates, if present
        std::copy(addcovar.begin(), addcovar.end(), X.begin() + g_size);

    for(int pos=0, offset=0; pos<n_pos; pos++, offset += g_size) {
        Rcpp::checkUserInterrupt();  // check for ^C from user

        // copy genoprobs for this pos into a matrix
        std::copy(genoprobs.begin() + offset, genoprobs.begin() + offset + g_size, X.begin());

        for(int phe=0; phe<n_phe; phe++) {
            // calc rss and paste into ith column of result
            result(phe,pos) = calc_ll_binreg_weighted(X, pheno(_,phe), weights, maxit, tol, qr_tol, eta_max);
        }
    }

    return result;
}

// calculate distance between two lists of points
//
// pt1 and pt2 are each matrices with two columns
// [[Rcpp::export]]
NumericMatrix calc_dist_pt2pt(NumericMatrix pt1, NumericMatrix pt2)
{
    int nr1 = pt1.rows(), nr2 = pt2.rows();
    int nc1 = pt1.cols(), nc2 = pt2.cols();
    if(nc1 != 2) throw std::invalid_argument("pt1 should have two columns");
    if(nc2 != 2) throw std::invalid_argument("pt2 should have two columns");

    NumericMatrix result(nr1,nr2);

    for(int i=0; i<nr1; i++) {
        for(int j=0; j<nr2; j++) {
            result(i,j) = calc_dist_pt2pt_one(pt1(i,_), pt2(j,_));
        }
    }

    return result;
}

// [[Rcpp::export]]
NumericVector Raw(NumericMatrix x) {
  NumericVector result(1);
  for (int i = 0; i < x.rows(); i++) {
    // It's the allocation of a new numeric vector, not the function
    // pointer, that slows down the other one.
    //NumericVector y = x(i, _);
    result = result + (x(i, 0) * x(i, 1));
  }
  return result;
}

// logistic regression by Qr decomposition with column pivoting
// return just the log likelihood
// [[Rcpp::export]]
double calc_ll_binreg_eigenqr(const NumericMatrix& X, const NumericVector& y,
                              const int maxit=100, const double tol=1e-6,
                              const double qr_tol=1e-12)
{
    const int n_ind = y.size();
    #ifndef RQTL2_NODEBUG
    if(n_ind != X.rows())
        throw std::invalid_argument("nrow(X) != length(y)");
    #endif

    double curllik = 0.0;
    NumericVector pi(n_ind), wt(n_ind), nu(n_ind), z(n_ind);

    for(int ind=0; ind<n_ind; ind++) {
        pi[ind] = (y[ind] + 0.5)/2;
        wt[ind] = sqrt(pi[ind] * (1-pi[ind]));
        nu[ind] = log(pi[ind]) - log(1-pi[ind]);
        z[ind] = nu[ind]*wt[ind] + (y[ind] - pi[ind])/wt[ind];
        curllik += y[ind] * log10(pi[ind]) + (1.0-y[ind])*log10(1.0-pi[ind]);
    }

    NumericMatrix XX = weighted_matrix(X, wt);

    bool converged=false;
    double llik=0.0;

    for(int it=0; it<maxit; it++) {
        Rcpp::checkUserInterrupt();  // check for ^C from user

        // fitted values using weighted XX; will need to divide by previous weights
        nu = calc_fitted_linreg_eigenqr(XX, z, qr_tol);

        llik = 0.0;
        for(int ind=0; ind<n_ind; ind++) {
            nu[ind] /= wt[ind]; // need to divide by previous weights
            pi[ind] = exp(nu[ind])/(1.0 + exp(nu[ind]));
            wt[ind] = sqrt(pi[ind] * (1.0 - pi[ind]));
            z[ind] = nu[ind]*wt[ind] + (y[ind] - pi[ind])/wt[ind];
            llik += y[ind] * log10(pi[ind]) + (1.0-y[ind])*log10(1.0-pi[ind]);
        }

        if(fabs(llik - curllik) < tol) { // converged
            converged = true;
            break;
        }

        XX = weighted_matrix(X, wt);
        curllik = llik;

    } // end iterations

    if(!converged) r_warning("binreg didn't converge");

    return llik;
}

//[[Rcpp::export]]
SEXP do_getcq_dense( NumericMatrix X, const IntegerVector& mcs0idx){

  List vnl = clone(List(X.attr("dimnames")));
  CharacterVector vn=vnl[0];

  int nrX = X.rows(), i, ii, j, k, l, past;
  
  IntegerVector pas( nrX ), vec_s( nrX ), vec2_s( nrX );
  IntegerVector ggg( nrX );
  
  // pas.setZero();
  for (i=0; i<nrX; ++i){
    j = mcs0idx[i];  // Rprintf("i=%d, j=%d, past=%d\n", i, j, past);
    vec_s = X(_, j );
    vec2_s = vec_s * pas ;
    past = sum( vec2_s );
    pas[i] = 1;
    ggg[ mcs0idx[i] ] = past;
  }
  
  // // cout << vec_s.transpose() << endl; cout << pas.transpose() << endl;

  IntegerVector ladder( nrX );
  for (i=0; i<nrX-1; ++i){
    if( ggg[i]+1>ggg[i+1]) ladder[i] = 1;
  }
  ladder[nrX-1]=1; //Rprintf("ladder: "); Rf_PrintValue( ladder );
  int ncq = sum( ladder );
  List cqlist(ncq);
  for (i=0; i<nrX; ++i) pas[i]=0;
  // pas.setZero();
  l=0;
  for (i=0; i<nrX; ++i){
    if (ladder[i]>0){
      j = mcs0idx[i];
      vec_s  = X(_, j );
      vec2_s = vec_s * pas;
      past = sum( vec2_s ) ;   //Rprintf("i=%d, j=%d, past=%d\n", i, j, past);
      IntegerVector cq(past+1);
      //cout << "vec2_s " << vec2_s.transpose() << endl; Rf_PrintValue( cq );
      k=0;
      for (ii=0; ii<nrX; ++ii){
	if (vec2_s[ii] != 0)
	  cq[k++] = ii;
      }
      cq[past] = j;
      CharacterVector cq2(past+1);
      for (k=0; k<past+1;++k) cq2[k]=vn[cq[k]];
      cqlist[l++] = cq2;     //Rf_PrintValue( cq );
    }
    pas[i] = 1;
  }
  return cqlist; //List::create( cqlist );
  //return List::create(1);
}

// [[Rcpp::export]]
NumericVector row_kth(NumericMatrix toSort, int k) {
  int n = toSort.rows();
  NumericVector meds = NumericVector(n);
  for (int i = 0; i < n; i++) {
    NumericVector curRow = toSort.row(i);
    std::nth_element(curRow.begin(), curRow.begin() + k, curRow.end());
    meds[i] = curRow[k];
  }

  return meds;
}

// Scan a single chromosome with interactive covariates
// this version uses less memory but will be slower
// (since we need to work with each position, one at a time)
// 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 of numeric phenotypes (individuals x phenotypes)
//             (no missing data allowed)
// addcovar  = additive covariates (an intercept, at least)
// intcovar  = interactive covariates (should also be included in addcovar)
// weights   = vector of weights
//
// output    = matrix of residual sums of squares (RSS) (phenotypes x positions)
//
// [[Rcpp::export]]
NumericMatrix scan_binary_onechr_intcovar_weighted_lowmem(const NumericVector& genoprobs,
                                                          const NumericMatrix& pheno,
                                                          const NumericMatrix& addcovar,
                                                          const NumericMatrix& intcovar,
                                                          const NumericVector& weights,
                                                          const int maxit=100,
                                                          const double tol=1e-6,
                                                          const double qr_tol=1e-12,
                                                          const double eta_max=30.0)
{
    const int n_ind = pheno.rows();
    if(Rf_isNull(genoprobs.attr("dim")))
        throw std::invalid_argument("genoprobs should be a 3d array but has no dim attribute");
    const Dimension d = genoprobs.attr("dim");
    if(d.size() != 3)
        throw std::invalid_argument("genoprobs should be a 3d array");
    const int n_pos = d[2];
    const int n_phe = pheno.cols();
    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)");

    NumericMatrix result(n_phe, n_pos);

    for(int pos=0; pos<n_pos; pos++) {
        Rcpp::checkUserInterrupt();  // check for ^C from user

        // form X matrix
        NumericMatrix X = formX_intcovar(genoprobs, addcovar, intcovar, pos, true);

        for(int phe=0; phe<n_phe; phe++) {
            // do regression
            result(phe,pos) = calc_ll_binreg_weighted(X, pheno(_,phe), weights, maxit, tol, qr_tol, eta_max);
        }
    }

    return result;
}

// use calc_resid_linreg for a 3-dim array
// [[Rcpp::export]]
NumericVector calc_resid_linreg_3d(const NumericMatrix& X, const NumericVector& P)
{
    int nrowx = X.rows();
    int sizep = P.size();

    NumericMatrix pr(nrowx, sizep/nrowx);
    std::copy(P.begin(), P.end(), pr.begin()); // FIXME I shouldn't need to copy

    NumericMatrix result = calc_resid_linreg(X, pr);
    result.attr("dim") = P.attr("dim");

    return result;
}

// use calc_resid_linreg for a 3-dim array
// [[Rcpp::export]]
NumericVector calc_resid_linreg_3d(const NumericMatrix& X, const NumericVector& P,
                                   const double tol=1e-12)
{
    const unsigned int nrowx = X.rows();
    const Dimension d = P.attr("dim");
    if(d[0] != nrowx)
        throw std::range_error("nrow(X) != nrow(P)");

    NumericMatrix pr(nrowx, d[1]*d[2]);
    std::copy(P.begin(), P.end(), pr.begin()); // FIXME I shouldn't need to copy

    NumericMatrix result = calc_resid_eigenqr(X, pr, tol);
    result.attr("dim") = d;

    return result;
}

// [[Rcpp::export]]
NumericMatrix term_score(NumericMatrix beta) {
  // calculate term score (Blei and Lafferty 2009)
  int W = beta.rows();
  int K = beta.cols();
  NumericMatrix term_score(W, K);
  
  for (int w = 0; w < W; w++) {
    double product = 0.0;
    for (int k = 0; k < K; k++) {
      product *= beta(w, k);
    }
    for (int k = 0; k < K; k++) {
      term_score(w, k) = beta(w, k) * log(beta(w, k) / pow(product, (1 / K)));
    }
  }
  
  return term_score;
}