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

static Rcpp::IntegerVector nz_vec(Rcpp::NumericMatrix alpha_new,
				  Rcpp::NumericMatrix eta_new,
				  Rcpp::NumericVector d_new,
				  double eps)
{
  int K = alpha_new.ncol();
  int p = alpha_new.nrow();
  int L = eta_new.nrow();
  
  Rcpp::IntegerVector result(p*K + L*K + L);
 	
  for (int i = 0; i < p*K; i++) {
    result[i] = nz(alpha_new[i],eps);
  }

  for (int i = 0; i < L*K; i++) {
    result[p*K + i] = nz(eta_new[i],eps);
  }

  for (int i = 0; i < L; i++) {
    result[p*K + L*K + i] = nz(d_new[i],eps);
  }

  return result;
}

// [[Rcpp::export]]
SEXP FilterIndepOld  (Rcpp::NumericVector y_rcpp,
					Rcpp::NumericMatrix y_lagged_rcpp,
					Rcpp::NumericMatrix z_dependent_rcpp,
					Rcpp::NumericMatrix z_independent_rcpp,
					Rcpp::NumericVector beta_rcpp,
					Rcpp::NumericVector mu_rcpp,
					Rcpp::NumericVector sigma_rcpp,
					Rcpp::NumericMatrix gamma_dependent_rcpp,
					Rcpp::NumericVector gamma_independent_rcpp,
					Rcpp::NumericMatrix transition_probs_rcpp,
					Rcpp::NumericVector initial_dist_rcpp
					)
{
	int n = y_rcpp.size();
	int M = mu_rcpp.size();
	arma::mat xi_k_t(M, n); // make a transpose first for easier column operations.

	arma::colvec y(y_rcpp.begin(), y_rcpp.size(), false);
	arma::mat    y_lagged(y_lagged_rcpp.begin(),
								y_lagged_rcpp.nrow(), y_lagged_rcpp.ncol(), false);
	arma::mat    z_dependent(z_dependent_rcpp.begin(),
								z_dependent_rcpp.nrow(), z_dependent_rcpp.ncol(), false);
	arma::mat    z_independent(z_independent_rcpp.begin(),
								z_independent_rcpp.nrow(), z_independent_rcpp.ncol(), false);
	arma::colvec beta(beta_rcpp.begin(), beta_rcpp.size(), false);
	arma::colvec mu(mu_rcpp.begin(), mu_rcpp.size(), false);
	arma::colvec sigma(sigma_rcpp.begin(), sigma_rcpp.size(), false);
	arma::mat    gamma_dependent(gamma_dependent_rcpp.begin(),
								gamma_dependent_rcpp.nrow(), gamma_dependent_rcpp.ncol(), false);
	arma::colvec gamma_independent(gamma_independent_rcpp.begin(),
								gamma_independent_rcpp.size(), false);
	arma::mat    transition_probs(transition_probs_rcpp.begin(),
								transition_probs_rcpp.nrow(), transition_probs_rcpp.ncol(), false);
	arma::colvec initial_dist(initial_dist_rcpp.begin(), initial_dist_rcpp.size(), false);

	double likelihood = 0;

	SEXP eta_rcpp = EtaIndep(y_rcpp, y_lagged_rcpp,
						z_dependent_rcpp, z_independent_rcpp,
						beta_rcpp, mu_rcpp, sigma_rcpp,
						gamma_dependent_rcpp, gamma_independent_rcpp);
	arma::mat eta_t = (Rcpp::as<arma::mat>(eta_rcpp)).t();

	xi_k_t.col(0) = eta_t.col(0) % initial_dist;
	double total = sum(xi_k_t.col(0));
	xi_k_t.col(0) = xi_k_t.col(0) / total;
	likelihood += log(total);

	for (int k = 1; k < n; k++)
	{
	  	xi_k_t.col(k) = eta_t.col(k) % (transition_probs * xi_k_t.col(k-1));
	 	total = sum(xi_k_t.col(k));
		xi_k_t.col(k) = xi_k_t.col(k) / total;
		likelihood += log(total);
	}

	return Rcpp::List::create(Named("xi.k") = wrap(xi_k_t.t()),
														Named("likelihood") = wrap(likelihood));
}

// kernel Dist function on a Grid
// [[Rcpp::export]]
Rcpp::NumericVector
KdeDist(const Rcpp::NumericMatrix & X
      , const Rcpp::NumericMatrix & Grid
      , const double                h
      , const Rcpp::NumericVector & weight
      , const bool printProgress
	) {
	const unsigned sampleNum = X.nrow();
	const unsigned dimension = Grid.ncol();
	const unsigned gridNum = Grid.nrow();
	// first = sum K_h(X_i, X_j), second = K_h(x, x), third = sum K_h(x, X_i)
	std::vector< double > firstValue;
	const double second = 1.0;
	std::vector< double > thirdValue;
	double firstmean;
	Rcpp::NumericVector kdeDistValue(gridNum);
	int counter = 0, percentageFloor = 0;
	int totalCount = sampleNum + gridNum;

	if (printProgress) {
		printProgressFrame(Rprintf);
	}

	firstValue = computeKernel< std::vector< double > >(
			X, X, h, weight, printProgress, Rprintf, counter, totalCount,
			percentageFloor);

	if (dimension <= 1) {
		thirdValue = computeKernel< std::vector< double > >(
				X, Grid, h, weight, printProgress, Rprintf, counter, totalCount,
				percentageFloor);
	}
	else {
		thirdValue = computeGaussOuter< std::vector< double > >(
				X, Grid, h, weight, printProgress, Rprintf, counter, totalCount,
				percentageFloor);
	}

	if (weight.size() == 1) {
		firstmean = std::accumulate(firstValue.begin(), firstValue.end(), 0.0) / sampleNum;
	}
	else {
		firstmean = std::inner_product(
				firstValue.begin(), firstValue.end(), weight.begin(), 0.0) / 
				std::accumulate(weight.begin(), weight.end(), 0.0);
	}

	for (unsigned gridIdx = 0; gridIdx < gridNum; ++gridIdx) {
		kdeDistValue[gridIdx] = std::sqrt(firstmean + second - 2 * thirdValue[gridIdx]);
	}

	if (printProgress) {
		Rprintf("\n");
	}

	return kdeDistValue;
}

// [[Rcpp::export]]
Rcpp::List SIMRE(int n1, Rcpp::NumericMatrix OMEGA, int n2, Rcpp::NumericMatrix SIGMA, int seed) {

  arma::mat eta;
  arma::mat eps;
  if(OMEGA.nrow() > 0) eta = MVGAUSS(OMEGA,n1,-1);
  if(SIGMA.nrow() > 0) eps = MVGAUSS(SIGMA,n2,-1);

  Rcpp::List ans;
  ans["eta"] = eta;
  ans["eps"] = eps;

  return(ans);
}

// [[Rcpp::export]]
Rcpp::NumericMatrix rcpp_sweep_(Rcpp::NumericMatrix x, Rcpp::NumericVector vec)
{
  Rcpp::NumericMatrix ret(x.nrow(), x.ncol());
  
  #pragma omp parallel for default(shared)
  for (int j=0; j<x.ncol(); j++)
  {
    #pragma omp simd
    for (int i=0; i<x.nrow(); i++)
      ret(i, j) = x(i, j) - vec(i);
  }
  
  return ret;
}

std::vector<std::vector<float> > matrix_to_array_v(Rcpp::NumericMatrix& mat)
{
  std::vector<std::vector<float> > v;
  if(!mat.ncol() || !mat.nrow())
    return(v);
  v.resize(mat.nrow());
  std::vector<float> v_row(mat.ncol()); // = new float[ mat.ncol() * mat.nrow() ];
  for(unsigned int i=0; i < (unsigned int)mat.nrow(); ++i){
    for(unsigned int j=0; j < (unsigned int)mat.ncol(); ++j){
      v_row[j] = (float)mat(i, j);
    }
    v[i] = v_row;
  }
  return(v);
}

Rcpp::List RadiusSearch(Rcpp::NumericMatrix query_,
                        Rcpp::NumericMatrix ref_,
                        double radius,
                        int max_neighbour,
                        std::string build,
                        int cores,
                        int checks) {
  const std::size_t n_dim = query_.ncol();
  const std::size_t n_query = query_.nrow();
  const std::size_t n_ref = ref_.nrow();
  // Column major to row major
  arma::mat query(n_dim, n_query);
  {
    arma::mat temp_q(query_.begin(), n_query, n_dim, false);
    query = arma::trans(temp_q);
  }
  flann::Matrix<double> q_flann(query.memptr(), n_query, n_dim);
  arma::mat ref(n_dim, n_ref);
  {
    arma::mat temp_r(ref_.begin(), n_ref, n_dim, false);
    ref = arma::trans(temp_r);
  }
  flann::Matrix<double> ref_flann(ref.memptr(), n_ref, n_dim);
  // Setting the flann index params
  flann::IndexParams params;
  if (build == "kdtree") {
    params = flann::KDTreeSingleIndexParams(1);
  } else if (build == "kmeans") {
    params = flann::KMeansIndexParams(2, 10, flann::FLANN_CENTERS_RANDOM, 0.2);
  } else if (build == "linear") {
    params = flann::LinearIndexParams();
  }
  // Perform the radius search
  flann::Index<flann::L2<double> > index(ref_flann, params);
  index.buildIndex();
  std::vector< std::vector<int> >
      indices_flann(n_query, std::vector<int>(max_neighbour));
  std::vector< std::vector<double> >
      dists_flann(n_query, std::vector<double>(max_neighbour));
  flann::SearchParams search_params;
  search_params.cores = cores;
  search_params.checks = checks;
  search_params.max_neighbors = max_neighbour;
  index.radiusSearch(q_flann, indices_flann, dists_flann, radius,
                     search_params);
  return Rcpp::List::create(Rcpp::Named("indices") = indices_flann,
                            Rcpp::Named("distances") = dists_flann);
}

void ScoreGaussL0PenScatter::setData(Rcpp::List& data)
{
	std::vector<int>::iterator vi;
	//uint i;

	// Cast preprocessed data from R list
	dout.level(2) << "Casting preprocessed data...\n";
	_dataCount = Rcpp::as<std::vector<int> >(data["data.count"]);
	dout.level(3) << "# samples per vertex: " << _dataCount << "\n";
	_totalDataCount = Rcpp::as<uint>(data["total.data.count"]);
	dout.level(3) << "Total # samples: " << _totalDataCount << "\n";
	Rcpp::List scatter = data["scatter"];
	Rcpp::NumericMatrix scatterMat;
	_disjointScatterMatrices.resize(scatter.size());
	dout.level(3) << "# disjoint scatter matrices: " << scatter.size() << "\n";
	for (R_len_t i = 0; i < scatter.size(); ++i) {
		scatterMat = Rcpp::NumericMatrix((SEXP)(scatter[i]));
		_disjointScatterMatrices[i] = arma::mat(scatterMat.begin(), scatterMat.nrow(), scatterMat.ncol(), false);
	}

	// Cast index of scatter matrices, adjust R indexing convention to C++
	std::vector<int> scatterIndex = Rcpp::as<std::vector<int> >(data["scatter.index"]);
	for (std::size_t i = 0; i < scatterIndex.size(); ++i)
		_scatterMatrices[i] = &(_disjointScatterMatrices[scatterIndex[i] - 1]);

	// Cast lambda: penalty constant
	_lambda = Rcpp::as<double>(data["lambda"]);
	dout.level(3) << "Penalty parameter lambda: " << _lambda << "\n";

	// Check whether an intercept should be calculated
	_allowIntercept = Rcpp::as<bool>(data["intercept"]);
	dout.level(3) << "Include intercept: " << _allowIntercept << "\n";
}

//' Marginal correlation matrix
//' 
//' Various workhorse functions to compute the marginal (or unconditional) 
//' correlations (and cross-correlation) estimates efficiently. 
//' They are (almost) 
//' equivalent implementations of \code{\link[stats]{cor}} in Rcpp, 
//' RcppArmadillo, and RcppEigen.
//' 
//' @rdname corFamily
//' @aliases corFamily
//'   corRcpp xcorRcpp corArma xcorArma corEigen xcorEigen
//' @param X A numeric matrix.
//' @param Y A numeric matrix of compatible dimension with the \code{X}, i.e. 
//'   \code{nrow(X)} equals \code{nrow(Y)}.
//' @return
//'   The \code{corXX} family returns a numeric correlation matrix of size 
//'   \code{ncol(X)} times \code{ncol(X)}.
//'   
//'   The \code{xcorXX} family returns a numeric cross-correlation matrix 
//'   of size \code{ncol(X)} times \code{ncol(Y)}.
//' @details
//'   Functions almost like \code{\link{cor}}.
//'   For the \code{xcorXX} functions, the \code{i}'th and \code{j}'th 
//'   entry of the output matrix is the correlation between \code{X[i, ]} and 
//'   \code{X[j, ]}.
//'   Likewise, for the \code{xcorXX} functions, the \code{i}'th and
//'   \code{j}'th entry of the output is the correlation between \code{X[i, ]} 
//'   and \code{Y[j, ]}.
//' @note 
//'   \code{NA}s in \code{X} or \code{Y} will yield \code{NA}s in the correlation matrix.
//'   This also includes the diagonal unlike the behavior of 
//'   \code{\link[stats]{cor}}.
//' @author Anders Ellern Bilgrau <anders.ellern.bilgrau (at) gmail.com>
//' @export
// [[Rcpp::export]]
Rcpp::NumericMatrix corRcpp(Rcpp::NumericMatrix & X) {
  
  const int m = X.ncol();
  const int n = X.nrow();
  
  // Centering the matrix
  X = centerNumericMatrix(X);
  
  Rcpp::NumericMatrix cor(m, m);
  
  // Degenerate case
  if (n == 0) {
    std::fill(cor.begin(), cor.end(), Rcpp::NumericVector::get_na());
    return cor; 
  }
  
  // Compute 1 over the sample standard deviation
  Rcpp::NumericVector inv_sqrt_ss(m);
  for (int i = 0; i < m; ++i) {
    inv_sqrt_ss(i) = 1/sqrt(Rcpp::sum(X(Rcpp::_, i)*X(Rcpp::_, i)));
  }
  
  // Computing the correlation matrix
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j <= i; ++j) {
      cor(i, j) = Rcpp::sum(X(Rcpp::_,i)*X(Rcpp::_,j)) *
        inv_sqrt_ss(i) * inv_sqrt_ss(j);
      cor(j, i) = cor(i, j);
    }
  }

  return cor;
}

static double bic_logistic(Rcpp::NumericMatrix X, 
			   Rcpp::NumericVector y, 
			   Rcpp::NumericMatrix beta_new, 
			   double eps, 
			   Rcpp::IntegerVector nk)
{
  int n_tot = X.nrow();
  int p = X.ncol();
  int K = nk.size();
    
  int idx = 0;
  double ll = 0.0;
  for (int k = 0; k < K; k++) {
    int n = nk[k];
    for (int i = 0; i < n; i++) {
      double lp = 0.0;
      for (int j = 0; j < p; j++) {
	lp += elem(X, idx+i, j) * elem(beta_new, j, k);
      }
      ll += y[idx+i] * lp - log(1.0 + exp(lp));
    }
    idx += n;
  }
  
  double bic = -2.0 * ll + df(beta_new, eps) * log(n_tot);
  return bic;
}

///********************************************************************
///**  cdm_rcpp_irt_classify_individuals
// [[Rcpp::export]]
Rcpp::List cdm_rcpp_irt_classify_individuals( Rcpp::NumericMatrix like )
{
    int N = like.nrow();
    int TP = like.ncol();
    Rcpp::IntegerVector class_index(N);
    Rcpp::NumericVector class_maxval(N);
    double val=0;
    int ind=0;
    for (int nn=0; nn<N; nn++){
        val=0;
        ind=0;
        for (int tt=0; tt<TP; tt++){
            if ( like(nn,tt) > val ){
                val = like(nn,tt);
                ind = tt;
            }
        }
        class_index[nn] = ind + 1;
        class_maxval[nn] = val;
    }
    //---- OUTPUT:
    return Rcpp::List::create(
                Rcpp::Named("class_index") = class_index,
                Rcpp::Named("class_maxval") = class_maxval
        );
}

// [[Rcpp::export]]
Rcpp::List europeanOptionArraysEngine(std::string type, Rcpp::NumericMatrix par) {

    QuantLib::Option::Type optionType = getOptionType(type);
    int n = par.nrow();
    Rcpp::NumericVector value(n), delta(n), gamma(n), vega(n), theta(n), rho(n), divrho(n);

    QuantLib::Date today = QuantLib::Date::todaysDate();
    QuantLib::Settings::instance().evaluationDate() = today;

    QuantLib::DayCounter dc = QuantLib::Actual360();

    for (int i=0; i<n; i++) {

        double underlying    = par(i, 0);    // first column
        double strike        = par(i, 1);    // second column
        QuantLib::Spread dividendYield = par(i, 2);    // third column
        QuantLib::Rate riskFreeRate    = par(i, 3);    // fourth column
        QuantLib::Time maturity        = par(i, 4);    // fifth column
#ifdef QL_HIGH_RESOLUTION_DATE    
        // in minutes
        boost::posix_time::time_duration length = boost::posix_time::minutes(boost::uint64_t(maturity * 360 * 24 * 60)); 
#else
        int length           = int(maturity*360 + 0.5); // FIXME: this could be better
#endif
        double volatility    = par(i, 5);    // sixth column
    
        boost::shared_ptr<QuantLib::SimpleQuote> spot(new QuantLib::SimpleQuote( underlying ));
        boost::shared_ptr<QuantLib::SimpleQuote> vol(new QuantLib::SimpleQuote( volatility ));
        boost::shared_ptr<QuantLib::BlackVolTermStructure> volTS = flatVol(today, vol, dc);
        boost::shared_ptr<QuantLib::SimpleQuote> qRate(new QuantLib::SimpleQuote( dividendYield ));
        boost::shared_ptr<QuantLib::YieldTermStructure> qTS = flatRate(today, qRate, dc);
        boost::shared_ptr<QuantLib::SimpleQuote> rRate(new QuantLib::SimpleQuote( riskFreeRate ));
        boost::shared_ptr<QuantLib::YieldTermStructure> rTS = flatRate(today, rRate, dc);
        
#ifdef QL_HIGH_RESOLUTION_DATE
    QuantLib::Date exDate(today.dateTime() + length);
#else
    QuantLib::Date exDate = today + length;
#endif    
        boost::shared_ptr<QuantLib::Exercise> exercise(new QuantLib::EuropeanExercise(exDate));
        
        boost::shared_ptr<QuantLib::StrikedTypePayoff> payoff(new QuantLib::PlainVanillaPayoff(optionType, strike));
        boost::shared_ptr<QuantLib::VanillaOption> option = makeOption(payoff, exercise, spot, qTS, rTS, volTS);
        
        value[i]  = option->NPV();
        delta[i]  = option->delta();
        gamma[i]  = option->gamma();
        vega[i]   = option->vega();
        theta[i]  = option->theta();
        rho[i]    = option->rho();
        divrho[i] = option->dividendRho();
    }
    return Rcpp::List::create(Rcpp::Named("value")  = value,
                              Rcpp::Named("delta")  = delta,
                              Rcpp::Named("gamma")  = gamma,
                              Rcpp::Named("vega")   = vega,
                              Rcpp::Named("theta")  = theta,
                              Rcpp::Named("rho")    = rho,
                              Rcpp::Named("divRho") = divrho);
}

static double bic_linear(Rcpp::NumericMatrix X, 
			 Rcpp::NumericVector y, 
			 Rcpp::NumericMatrix beta_new, 
			 double eps, 
			 Rcpp::IntegerVector nk)
{
  int n_tot = X.nrow();
  int p = X.ncol();
  int K = nk.size();
 

  /*calculate SSe*/
  double SSe = 0.0;  
  int idx = 0;
  
  for (int k = 0; k < K; k++) {
    int n = nk[k];
    for (int i = 0; i < n; i++) {
      double Xrow_betacol = 0.0;
      for (int j = 0; j < p; j++) {
	Xrow_betacol += elem(X, idx+i, j) * elem(beta_new, j, k);
      }
      SSe += pow(y[idx+i] - Xrow_betacol, 2);
    }
    idx += n;
  }
  
  double ll = -n_tot / 2.0 * (log(SSe) - log(n_tot) + log(2.0 * M_PI) + 1);
  double bic = -2 * ll + df(beta_new, eps) * log(n_tot);

  return bic;
}

//' Evaluate an h-function corresponding to a copula density estimate
//' 
//' @param uev mx2 matrix of evaluation points
//' @param cond_var either 1 or 2; the variable to condition on.
//' @param vals matrix of density estimate evaluated on a kxk grid.
//' @param grid the grid points (1-dim) on which vals has been computed.
//' 
//' @return H-function estimate evaluated at uev.
//' 
//' @noRd
// [[Rcpp::export]]
Rcpp::NumericVector eval_hfunc_2d(const Rcpp::NumericMatrix& uev,
                                  const int& cond_var, 
                                  const Rcpp::NumericMatrix& vals,
                                  const Rcpp::NumericVector& grid) 
{
    int N = uev.nrow();
    int m = grid.size();
    NumericVector tmpvals(m), out(N), tmpa(4), tmpb(4);
    NumericMatrix tmpgrid(m, 2);
    double upr = 0.0;
    double tmpint, int1; 
    tmpint = 0.0;
    
    for (int n = 0; n < N; ++n) {
        if (cond_var == 1) {
            upr = uev(n, 1);
            tmpgrid(_, 0) = rep(uev(n, 0), m);
            tmpgrid(_, 1) = grid;
        } else if (cond_var == 2) {
            upr = uev(n, 0);
            tmpgrid(_, 0) = grid;
            tmpgrid(_, 1) = rep(uev(n, 1), m);
        }
        tmpvals = interp_2d(tmpgrid, vals, grid, tmpa, tmpb);
        tmpint = int_on_grid(upr, tmpvals, grid);
        int1 =  int_on_grid(1.0, tmpvals, grid);
        out[n] = tmpint/int1;
        out[n] = fmax(out[n], 1e-10);
        out[n] = fmin(out[n], 1-1e-10);
    }
    return out;
}

///********************************************************************
///** scale2_NA_C
// [[Rcpp::export]]
Rcpp::NumericMatrix scale2_NA_C( Rcpp::NumericMatrix x )
{
    int n = x.nrow();
    int p = x.ncol();
    Rcpp::NumericMatrix y(n,p);
    double mvv=0;
    double sdvv=0;
    double nvv=0;
    double eps_add = 1e-10;
    for (int vv=0;vv<p;vv++){
        mvv=0;
        sdvv=0;
        nvv=0;
        for (int ii=0;ii<n;ii++){
            if (! R_IsNA(x(ii,vv)) ) {
                mvv += x(ii,vv);
                sdvv += std::pow( x(ii,vv), 2.0 );
                nvv ++;
            }
        }
        mvv = mvv / nvv;
        sdvv = std::sqrt( ( sdvv - nvv * mvv*mvv )/(nvv - 1.0 ) );
        // define standardization
        y(_,vv) = ( x(_,vv) - mvv ) / ( sdvv + eps_add );
    }
    //--- output
    return y;
}