C++ mat::t方法代码示例

void CosineTreeBuilder::CTNode(arma::mat A, CosineTree& root)
{
  A = A.t();
  Log::Info<<"CTNode"<<std::endl;
  //Calculating Centroid
  arma::rowvec centroid = CalculateCentroid(A);
  //Calculating sampling probabilities
  arma::vec probabilities = arma::zeros<arma::vec>(A.n_rows,1);
  LSSampling(A,probabilities);
  //Setting Values
  root.Probabilities(probabilities);
  root.Data(A);
  root.Centroid(centroid);
}

// [[Rcpp::export]]
arma::vec innerLoop(arma::vec resp,
                    arma::vec beta, double intercept,
                    double tol, int max_iter,
                    arma::mat x_mat,
                    int n, arma::vec weights) {
  bool converge = false;
  int counter = 0;

  arma::vec temp_resp;
  arma::vec beta_new = beta;
  double intercept_new;

  while(!converge && counter < max_iter) {

    temp_resp = x_mat.t() * ((resp - intercept) / n);

    // The argmin is hence equivalent to that of
    // (1/2) * ||v - beta||_2^2 + 4 * lam * Pen(beta).
    // In main func weights will be weights.col(i)

    beta_new = GetProxOne(temp_resp, 4 * weights);
    intercept_new =  mean(resp - x_mat * beta_new);


    double change1 = pow(intercept_new - intercept, 2) + sum(square(beta_new - beta));
    double change2 = pow(intercept_new, 2) + sum(square(beta_new));;


    if( pow(change1, 0.5) / pow(change2, 0.5) < tol ) {
      beta = beta_new;
      intercept = intercept_new;
      converge = true;
    } else {
      beta = beta_new;
      intercept  = intercept_new;
      counter = counter + 1;
      if(counter == max_iter) {
        beta = beta_new;
        intercept = intercept_new;
        Function warning("warning");
        warning("Function did not converge for inner loop for some lambda.");
      }
    }

  }
  arma::vec inter_vec(1);
  inter_vec(0) = intercept;
  return join_vert(inter_vec, beta);

}

// [[Rcpp::export]]
arma::mat expskewC(arma::mat M){
  /*This function takes a 3-by-3 skew symmetric matrix (in so(3)) and
  returns the exponential, a 3-by-3 rotations (in SO(3))*/
  
  double MMt = sum(sum(M-M.t()));
  
  if(fabs(MMt)>0.01){
    throw Rcpp::exception("The exp.skew function is expecting a 3-by-3 skew symmetric matrix");
  }
  
  arma::mat expM(3,3);
  expM.eye();
  
  double a = pow(0.5*trace(M.t()*M),0.5);
  
  if(a < 0.000001 && a > -0.000001){
    return expM;
  }
   
  expM = expM + (sin(a)/a) * M + (1-cos(a))*pow(a,-2)*M*M;
  
  return expM;
  
}

//evaluate Hessian operator on Z
//eucH is the ordinary Hessian matrix on Z in the ambient space
//return <Z, Hessian*Z>_Y
double specialLinear::evalHessian(const arma::mat & eucH,const arma::mat & Z){
  arma::mat YU,eucH_proj,Weingarten,Weingarten_part,temp;
  //project euclidian Hessian onto tangent space
  YU=eucH*Y.i();
  YU=1.0/n*(arma::trace(YU))*Y;
  eucH_proj=eucH-YU;
  //Weingarten Map
  Weingarten=-1.0/n*(arma::dot(eucH.t(),Y.i()))*Z;
  Weingarten_part=1.0/n*eucH*Y.i()*Z;
  temp=1.0/n*(arma::dot(Weingarten_part.t(),Y.i()))*Y;
  Weingarten_part-=temp;
  Weingarten+=Weingarten_part;
  hessian_Z=eucH_proj+Weingarten;
  Z_hessian_Z=arma::dot(hessian_Z,Z);
  return Z_hessian_Z;
}

arma::vec cpp_UpdateGamma (arma::mat x_big, arma::vec weights_big,
                           arma::vec y_big, 
                           arma::vec gamma, arma::vec beta, 
                           double lambda) {
  // Another helper function used by the main algorithm.
  // In this algorithm we loop through and do a coordinate wise update of
  // gamma.
  // 
  // Agrs: 
  //    beta/gamma: The 'current' parameter vectors.
  //    x/weights/y_big: Design, weights and response for the big population.
  //    lambda: The value of the penalty parameter.
  // Returns:
  //    gamma: The updated gamma vector. 
  

  // Obtain value of dimension.
  int p = gamma.size();
  
  // Initialize matrix X * beta. 
  arma::mat x_beta = x_big * beta;
  
  // We need the diagonal entries of X^T * W * X.
  // This is easy when W is a diagonal matrix.
  arma::mat W_times_X(x_big.begin(), x_big.n_rows, x_big.n_cols, true);
  W_times_X.each_col() %= weights_big;
  // We obtain the diagonal entries of X^T * W * X for the denominators.
  arma::mat Xt_W_X = x_big.t() * W_times_X;
  arma::vec denoms = Xt_W_X.diag();
  
  // Some temporary vectors which we use later.
  arma::vec temp_gamma(gamma.begin(), gamma.size(), true);
  double temp1;
  
  // Loop through the different gamma values which need to be updated.
  for(int i = 0; i < p; i++) {
    // Instead of excluding a column we simply set gamma_j = 0 to get the fitted
    // value without the effect of gamma_j.
    temp_gamma(i) = 0;
    temp1 = as_scalar(x_big.col(i).t() * (weights_big % (y_big - x_beta - x_big * temp_gamma)));
    gamma(i) = (cpp_sign(temp1) * cpp_max(abs(temp1) - lambda, 0.0))/denoms(i);
    
    temp_gamma(i) = gamma(i);
    
  }
  return gamma;
}

// [[Rcpp::export]]
arma::mat logSO3C(arma::mat R){
  
  arma::mat I(3,3), logR(3,3);
  I.eye();
  
  double theta = acos(0.5*trace(R)-0.5);
  
  if(theta < 0.0001){
    logR.zeros();
    return logR;
  }
  
  logR = (R-R.t())*theta/(2*sin(theta));
  
  return logR;
  
}

// Uses singular value decomposition function svd() from Armadillo.
// Thus, be sure to define the following in /path/to/armadillo/include/armadillo_bits/config.hpp:
// ARMA_USE_LAPACK
// ARMA_USE_BLAS
// ARMA_BLAS_UNDERSCORE
int PCA::Run(const DataUnlabeledNormalized &data_unlabeled_normalized) {
  const int kNumTrainEx = data_unlabeled_normalized.num_train_ex();
  assert(kNumTrainEx >= 1);
  const arma::mat kTrainingFeatures = \
    data_unlabeled_normalized.training_features_normalized();
  const arma::mat kCovMat = \
    (1.0/(float)kNumTrainEx)*kTrainingFeatures.t()*kTrainingFeatures;
  const int kNumRows = kCovMat.n_rows;
  arma::mat left_sing_vec = arma::zeros<arma::mat>(kNumRows,kNumRows);
  arma::vec sing_val = arma::zeros<arma::vec>(kNumRows,1);
  arma::mat right_sing_vec = arma::zeros<arma::mat>(kNumRows,kNumRows);
  arma::svd(left_sing_vec,sing_val,right_sing_vec,kCovMat);
  this->set_left_sing_vec(left_sing_vec);
  this->set_sing_val(sing_val);

  return 0;
}

// // [[Rcpp::export()]]
arma::mat getEx2x2_ordIRT(const arma::mat &Ex,
                   const arma::mat &Vx,
                   const int N) {

    arma::mat Ex2x2(2, 2, arma::fill::zeros) ;
    arma::mat tmp(1, 1) ;
    tmp.fill(N) ;

    arma::mat S = tmp % Vx + Ex.t() * Ex ;

    Ex2x2(0, 0) = N ;
    Ex2x2(1, 1) = S(0,0) ;

    arma::mat sums = sum(Ex, 0) ;
    Ex2x2(1, 0) = sums(0,0) ;
    Ex2x2(0, 1) = sums(0,0) ;

    return(Ex2x2) ;

}

double
HMM::forwardProcedureCached() {

  //initialisation
  alpha_.col(0) = arma::trans(pi_ % B_.row(0));

  c_(0) = arma::accu(alpha_.col(0));
  alpha_.col(0) /= arma::as_scalar(c_(0));

  //alpha_.print("alpha");
  //c_.print("scale");
  //iteration
  for(unsigned int t = 1; t < T_; ++t) {
    alpha_.col(t) = (A_.t() * alpha_.col(t-1)) % arma::trans(B_.row(t));
    c_(t) = arma::accu(alpha_.col(t));
    alpha_.col(t) /= arma::as_scalar(c_(t)); 
  }

  pprob_ = arma::accu(arma::log(c_));
  return pprob_;
}

// [[Rcpp::export]]   
arma::rowvec meanQ4C(arma::mat Q) { 
	//Compute the projected mean of the sample Q
	
	NumericMatrix Qss = as<NumericMatrix>(wrap(Q));
	int cq4 = checkQ4(Qss);
	if(cq4){
		throw Rcpp::exception("The data are not in Q4.");
	}
	
	arma::mat Qsq=Q.t()*Q;
	arma::mat eigvec;
	arma::vec eigval;
  arma::eig_sym(eigval,eigvec,Qsq);   
  arma::vec qhat=eigvec.col(3);
  
  if(qhat[0]<0){
  	qhat = -qhat;
  }
  
  return qhat.t(); //Want to return it in a row vector so transpose it
}

/**
 * Compute
 *
 *     alpha = min(1, tau * alphahat(A, dA))
 *
 * where
 *
 *     alphahat = sup{ alphahat : A + dA is psd }
 *
 * See (2.18) of [AHO98] for more details.
 */
static inline double
Alpha(const arma::mat& A, const arma::mat& dA, double tau)
{
  // On Armadillo < 4.500, the "lower" option isn't available.
#if (ARMA_VERSION_MAJOR < 4) || \
    ((ARMA_VERSION_MAJOR == 4) && (ARMA_VERSION_MINOR < 500))
  const arma::mat L = arma::chol(A).t(); // This is less efficient.
#else
  const arma::mat L = arma::chol(A, "lower");
#endif
  const arma::mat Linv = arma::inv(arma::trimatl(L));
  // TODO(stephentu): We only want the top eigenvalue, we should
  // be able to do better than full eigen-decomposition.
  const arma::vec evals = arma::eig_sym(-Linv * dA * Linv.t());
  const double alphahatinv = evals(evals.n_elem - 1);
  double alphahat = 1. / alphahatinv;
  if (alphahat < 0.)
    // dA is PSD already
    alphahat = 1.;
  return std::min(1., tau * alphahat);
}

// [[Rcpp::export]]
arma::rowvec HnCpp(arma::mat Qs){
  //Compute the Hn tests statistics
  
  int n = Qs.n_rows, i=0;
  arma::mat T = Qs.t()*Qs;
  arma::mat eigvec, eigvecJ;
  arma::vec eigval, eigvalJ;
  arma::eig_sym(eigval,eigvec,T);
  arma::rowvec Hn(n);
  arma::rowvec Qj;
  arma::mat Tj;

  for(i = 0;i<n; i++){
    Qj = Qs.row(i);
    
    Tj = T-Qj.t()*Qj;
    arma::eig_sym(eigvalJ,eigvecJ,Tj);
    Hn(i)=(n-2)*(1+eigvalJ(3)-eigval(3))/(n-1-eigvalJ(3));
    
  }
  return Hn;
}

NumericVector logLikMixHMM(const arma::mat& transition, const arma::cube& emission,
  const arma::vec& init, const arma::ucube& obs, const arma::mat& coef, const arma::mat& X,
  const arma::uvec& numberOfStates, unsigned int threads) {
  
  arma::mat weights = exp(X * coef).t();
  if (!weights.is_finite()) {
    return wrap(-arma::datum::inf);
  }
  weights.each_row() /= sum(weights, 0);
  
  arma::vec ll(obs.n_slices);
  arma::sp_mat transition_t(transition.t());
#pragma omp parallel for if(obs.n_slices >= threads) schedule(static) num_threads(threads) \
  default(none) shared(ll, obs, weights, init, emission, transition_t, numberOfStates)
    for (unsigned int k = 0; k < obs.n_slices; k++) {
      arma::vec alpha = init % reparma(weights.col(k), numberOfStates);
      
      for (unsigned int r = 0; r < obs.n_rows; r++) {
        alpha %= emission.slice(r).col(obs(r, 0, k));
      }
      
      double tmp = sum(alpha);
      ll(k) = log(tmp);
      alpha /= tmp;
      
      for (unsigned int t = 1; t < obs.n_cols; t++) {
        alpha = transition_t * alpha;
        for (unsigned int r = 0; r < obs.n_rows; r++) {
          alpha %= emission.slice(r).col(obs(r, t, k));
        }
        
        tmp = sum(alpha);
        ll(k) += log(tmp);
        alpha /= tmp;
      }
    }
    return wrap(ll);
}

QUIC_SVD::QUIC_SVD(const arma::mat& dataset,
                   arma::mat& u,
                   arma::mat& v,
                   arma::mat& sigma,
                   const double epsilon,
                   const double delta) :
    dataset(dataset)
{
  // Since columns are sample in the implementation, the matrix is transposed if
  // necessary for maximum speedup.
  CosineTree* ctree;
  if (dataset.n_cols > dataset.n_rows)
    ctree = new CosineTree(dataset, epsilon, delta);
  else
    ctree = new CosineTree(dataset.t(), epsilon, delta);

  // Get subspace basis by creating the cosine tree.
  ctree->GetFinalBasis(basis);

  // Use the ExtractSVD algorithm mentioned in the paper to extract the SVD of
  // the original dataset in the obtained subspace.
  ExtractSVD(u, v, sigma);
}

// // [[Rcpp::export()]]
arma::mat getEx2x2(const arma::mat &Ex,
                   const arma::mat &Vx,
                   const int N,
                   const int D
                   ) {
    arma::mat Ex2x2(D + 1, D + 1, arma::fill::zeros) ;
    arma::mat tmp(D, D) ;
    tmp.fill(N) ;
    arma::mat S = tmp % Vx + Ex.t() * Ex ;

    // S.print("S") ;
    // (tmp % Vx).print("tmp") ;

    Ex2x2(0, 0) = N ;
    Ex2x2.submat(1, 1, D, D) = S ;

    arma::mat sums = sum(Ex, 0) ;
    for (int d = 0; d < D ; d++) {
        Ex2x2(d + 1, 0) = sums(0,d) ;
        Ex2x2(0, d + 1) = sums(0,d) ;
    }
    return(Ex2x2) ;
}