C++ NumericVector::attr方法代码示例

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

// [[Rcpp::export]]
NumericVector attribs() {
  NumericVector out = NumericVector::create(1, 2, 3);

  out.names() = CharacterVector::create("a", "b", "c");
  out.attr("my-attr") = "my-value";
  out.attr("class") = "my-class";

  return out;
}

  void defineVariable(NumericVector x, std::string name) {
    readstat_label_set_t* labelSet = NULL;
    if (rClass(x) == "labelled") {
      labelSet = readstat_add_label_set(writer_, READSTAT_TYPE_DOUBLE, name.c_str());

      NumericVector values = as<NumericVector>(x.attr("labels"));
      CharacterVector labels = as<CharacterVector>(values.attr("names"));

      for (int i = 0; i < values.size(); ++i)
        readstat_label_double_value(labelSet, values[i], std::string(labels[i]).c_str());
    }

    readstat_add_variable(writer_, READSTAT_TYPE_DOUBLE, 0, name.c_str(),
      var_label(x), NULL, labelSet);
  }

// [[Rcpp::export]]
arma::cube gibbsCPP(NumericVector p2, int T, int M, double rho, double x0, double y0) 
{
  IntegerVector dim_p2=p2.attr("dim");
  arma::cube p(p2.begin(), dim_p2[0], dim_p2[1], dim_p2[2]);
  double x, y;
  arma::vec n1 = rnormC(76543,M*T);
  //arma::vec n2 = rnormC(76543,M*T);
  for(int m=0;m<M;m++)
  {
    x = x0;
    y = y0;
    p(0,0,m) = x;
    p(1,0,m) = y;
    for(int t=1;t<T;t++)
    {
      //x = rnormC(1,M*T)[m*t];
      //x = x * sqrt(1-pow(rho,2)) + rho * y;
      //y = rnormC(76543,M*T)[(m*t)+1];
      //y = y * sqrt(1-pow(rho,2)) + rho * x;
      x = n1[m*T];
      x = x * sqrt(1-pow(rho,2)) + rho * y;
      y = n1[m*T+1];
      y = y * sqrt(1-pow(rho,2)) + rho * x;
      p(0,t,m) = x;
      p(1,t,m) = y; 
    }
  }
  return(p);
}

// [[Rcpp::export]]
NumericVector row_weights(NumericMatrix x, NumericVector weight) {
  int nX = x.ncol();
  int nY = x.nrow();
  NumericVector v = no_init(nY);
  
  #pragma omp parallel for schedule(static)
  for (int i=0; i < nY; i++) {
    NumericMatrix::Row row = x(i, _);

    double w = 0;
    for (int j=0; j < nX; j++) {
      if(row[j]!=0 && !R_IsNA(row[j])) {
        w += weight[j];
      }
    }

    double o = 0;
    if (w!=0) {
      for(int j=0; j < nX; j++) {
        if(row[j]!=0 && !R_IsNA(row[j])) {
          o += row[j]*weight[j] / w;
        }
      }
    }
    
    v[i] = o;
  }
  
  v.attr("names") = rownames(x);
  return v;
}

// [[Rcpp::export]]
List pwTabMerge(List hsum, NumericMatrix pw) {
    
    int n = pw.ncol() ;
    
    List out(n) ;
    
    for(int k = 0; k < n; k++){
      
      RCPP_UNORDERED_MAP<std::string,double> preout ;
      
      List ip(2) ;
      ip[0] = hsum[pw(0,k)] ;
      ip[1] = hsum[pw(1,k)] ;
      
      for(int i = 0; i < 2; i++){
        NumericVector x = ip[i] ;
        CharacterVector names = x.attr("names") ;
        
        int m = x.size() ;
        
        for(int j = 0; j < m; j++){
          String name = names[j] ;
          preout[ name ] += x[j] ;
        }
      }
      out[k] = preout ;
    }
  return wrap(out) ;
}

//' Do values in a numeric vector fall in specified range?
//'
//' This is a shortcut for `x >= left & x <= right`, implemented
//' efficiently in C++ for local values, and translated to the
//' appropriate SQL for remote tables.
//'
//' @param x A numeric vector of values
//' @param left,right Boundary values
//' @export
//' @examples
//' between(1:12, 7, 9)
//'
//' x <- rnorm(1e2)
//' x[between(x, -1, 1)]
// [[Rcpp::export]]
LogicalVector between(NumericVector x, double left, double right) {
  int n = x.size();
  LogicalVector out(no_init(n));

  // Assume users know what they're doing with date/times. In the future
  // should ensure that left and right are the correct class too.
  if (x.attr("class") != R_NilValue && !Rf_inherits(x, "Date") && !Rf_inherits(x, "POSIXct")) {
    warningcall(R_NilValue, "between() called on numeric vector with S3 class");
  }

  if (NumericVector::is_na(left) || NumericVector::is_na(right)) {
    for (int i = 0; i < n; ++i)
      out[i] = NA_LOGICAL;
    return out;
  }

  for (int i = 0; i < n; ++i) {
    if (NumericVector::is_na(x[i])) {
      out[i] = NA_LOGICAL;
    } else if ((x[i] >= left) && (x[i] <= right)) {
      out[i] = true;
    } else {
      out[i] = false;
    }
  }

  return out;
}

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

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

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

// 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 int nrowx = X.rows();
    if(Rf_isNull(P.attr("dim")))
        throw std::invalid_argument("P should be a 3d array but has no dim attribute");
    const Dimension d = P.attr("dim");
    if(d.size() != 3)
        throw std::invalid_argument("P should be a 3d array");
    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;
}

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

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

//' FITS image writer
//' 
//' Writes a vector, matrix or 3D array to a FITS file as an image.
//' The data is written to the primary HDU.
//' 
// [[Rcpp::export]]
int gv_writefits_img(NumericVector img, CharacterVector fits_name, CharacterVector hdu_name = "")
{
  IntegerVector dim;
  if (!img.hasAttribute("dim"))
  {
    REprintf("ERROR: image has not been dimensioned.\n");
    return 1;
  }
  dim = img.attr("dim");
  if (dim.length() > 3)
	{
		REprintf("ERROR: dimension of more than 3 unsupported.\n");
		return 1;
	}

  fitsfile *pfits=NULL;
	int err=0;

	std::string fname = as<std::string>(fits_name[0]);
	fits_create_file(&pfits, (char *) fname.c_str(), &err);
	if (err)
	{
		gv_print_fits_err(err);
		return err;
	}

#ifdef GV_DEBUG
  Rcout << "Number of dim: " << dim.length() << std::endl;
	for (int i=0; i<dim.length(); i++)
	{
		Rcout << "Dim[" << i << "]: " << dim[i] << std::endl;
	}
	Rcout << "Number of elements: " << img.length() << std::endl;
  double *p = &(*img.begin());
  for (int i=0; i<img.length(); i++)
  {
    Rcout << "*(p+" << i << ") = " << *(p+i) << std::endl;
  }
#endif
  
  long longdim[3], startpix[3] = {1,1,1}; // default start
  for (int i=0; i<dim.length(); i++) longdim[i] = (long) dim[i];
  
	// start writing to file
	fits_create_img(pfits, DOUBLE_IMG, dim.length(), longdim, &err);
	fits_write_pix(pfits, TDOUBLE, startpix, img.length(), &(*img.begin()), &err);

	fits_close_file(pfits, &err);
  return err;
}