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

//Function to generate a single list of two matrices (train or test), given an index
//  [[Rcpp::export]]
Rcpp::List MatSplit( const arma::mat & A, const arma::mat & B, const arma::uvec & splitind){
  arma::mat tA = A.rows(splitind);
  arma::mat tB = B.rows(splitind);
  return Rcpp::List::create(Rcpp::Named("A")=tA,
                            Rcpp::Named("B")=tB);

}

// this function sample the loading parameters considering constant 
//  observational variance.
// [[Rcpp::export]]
arma::mat SampleDfmLoads(arma::mat y, arma::mat factors, arma::vec psi, 
  int s, double c0)
{
  /**
   * The data matrix 'y' is (T \times q).
   * Note that the factors matrix is (T+h \times k). In the models considered
   *  in the dissertation, h > s. This means that we always have the 
   *  VAR order, 'h', greater than the number of dynamic factors, 's', in the 
   *  observational equation. 
   * The argument 's' is the order of the dynamic factors present in the
   *  observational equation.
   * The prior mean is set to zero.
   */
   
  int T = y.n_rows;
  int q = y.n_cols;;
  int k = factors.n_cols;
  int TpH = factors.n_rows;
  int h = TpH - T;
  if (h < s+1){
    throw std::range_error(
      "Argument 's' greater than or equal to factors' VAR order in SampleDfmLoads.cpp"
      );
  }
  arma::mat I_q = eye(q, q);
//  std::cout << "\nh = " << h << "\n";
  arma::mat Xlambda;
  Xlambda = factors.rows(h, TpH-1);
//  std::cout << "\nNumber of rows in Xlambda is " << Xlambda.n_rows << "\n";
  if (s > 0){
    for (int i = 1; i <= s; i++){
    Xlambda = arma::join_rows(Xlambda, factors.rows(h-i, TpH-i-1));
    }
  }
  
  // posterior variance
  arma::mat Eigvec;
  arma::vec eigval;
  
  arma::eig_sym(eigval, Eigvec, arma::symmatu(Xlambda.t() * Xlambda));
  
  arma::mat kronEig = kron(I_q, Eigvec);  
  arma::vec L1eigval = 1.0/(arma::kron((1.0/psi), eigval) + 1/c0);
  arma::mat L1 = kronEig * diagmat(L1eigval) * kronEig.t();
  
  // posterior mean
  arma::vec yStar(y.begin(),  T*q);
  arma::vec l1 = L1 * (arma::kron(arma::diagmat(1.0/psi), Xlambda.t()) * yStar);
  arma::mat l1Mat(l1.begin(), k*(s+1), q); // mean in matrix form
  
  // random generation
  arma::vec vecnorm = sqrt(L1eigval) % randn(q*k*(s+1), 1);  
  arma::mat LambdaBar(vecnorm.begin(), k*(s+1), q, false);
  LambdaBar = l1Mat + Eigvec * LambdaBar;
  
  return LambdaBar.t();
}

// [[Rcpp::export]]
arma::mat BeQTL2(const arma::mat & A, const arma::mat & B, const arma::umat & Bootmat){
  int bsi= Bootmat.n_rows;
  arma::mat C(A.n_cols*B.n_cols,Bootmat.n_rows);
  arma::mat tC(A.n_rows,B.n_rows);
  for(int i=0; i<bsi; i++){
    tC = cor(A.rows(Bootmat.row(i)),B.rows(Bootmat.row(i)));
    C.col(i) = vectorise(tC,0);
  }
  C.elem(find_nonfinite(C)).zeros();

  return reshape(median(C,1),A.n_cols,B.n_cols);
}

/******************************************************************
 * 函数功能：几何校正需要调用的函数
 * 
 *
 */
arma::mat centering(arma::mat &x){
	arma::mat tmp = -mean(x.rows(0, 1), 1);
	arma::mat tmp2(2,2);
	tmp2.eye();
	arma::mat tmp1 = join_horiz(tmp2, tmp);
	arma::mat v;
	v << 0 << 0 << 1 << arma::endr;
	arma::mat T = join_vert(tmp1, v);
	//T.print("T =");
	arma::mat xm = T * x;
	//xm.print("xm =");
	
	//at least one pixel apart to avoid numerical problems
	//xm.print("xm =");
	double std11 = arma::stddev(xm.row(0));
	//cout << "std11:" << std11 << endl;
	double std22 = stddev(xm.row(1));
	//cout << "std22:" << std22 << endl;

    double std1 = std::max(std11, 1.0);
    double std2 = std::max(std22, 1.0);
	
	arma::mat S;
	S << 1/std1 << 0 << 0 << arma::endr
	  << 0 << 1/std2 << 0 << arma::endr
	  << 0 << 0 << 1 << arma::endr;
	arma::mat C = S * T ;
	//C.print("C =");
	return C;
}

// **********************************************************//
//                  Calculate Dyadic statistics              //
// **********************************************************//
// [[Rcpp::export]]
arma::cube degreestat (arma::cube cache, arma::mat intervals_d, int a, int A, arma::mat pd) {
	int nIP = pd.n_cols;
	arma::cube out = arma::zeros<arma::cube>(A, nIP, 6);
	for (unsigned int i = 0; i < 3; i++) {
		int mins = intervals_d(i, 0)-1;
		int maxs = intervals_d(i, 1)-1;
		arma::rowvec sendr = arma::zeros<arma::rowvec>(nIP);
		arma::mat receiver = arma::zeros<arma::mat>(A, nIP);
		arma::mat pdnew = pd.rows(mins, maxs);
		for (int r = 0; r < A; r++) {
			arma::colvec docsar = cache(arma::span(a), arma::span(r), arma::span(mins, maxs));
			sendr += docsar.t() * pdnew;
			for (int h = 0; h < A; h++) {
				arma::colvec docsra = cache(arma::span(h), arma::span(r), arma::span(mins, maxs));
				receiver.row(r) += docsra.t() * pdnew;	
				}	
		}
		for (int IP = 0; IP < nIP; IP++) {
			for (int r = 0; r < A; r++) {
				out(r, IP, i) = sendr[IP];
				out(r, IP, i+3) = receiver(r,IP);
			}
		}		
	}	
	return out;
}

//Script that takes two matrices, performs bootstrapped correlation, and returns the median
// [[Rcpp::export]]
arma::mat BeQTL(const arma::mat & A, const arma::mat & B, const arma::umat & Bootmat){
  int bsi= Bootmat.n_rows;
  Rcpp::Rcout<<"Starting Bootstrap!"<<std::endl;
  arma::mat C(A.n_cols*B.n_cols,Bootmat.n_rows);
  arma::mat tA(A.n_rows,A.n_cols);
  arma::mat tB(B.n_rows,B.n_cols);
  arma::mat tC(A.n_rows,B.n_rows);
  for(int i=0; i<bsi; i++){
    tA = A.rows(Bootmat.row(i));
    tB = B.rows(Bootmat.row(i));
    tC = cor(tA,tB);
    C.col(i) = vectorise(tC,0);
  }
  C.elem(find_nonfinite(C)).zeros();

 return reshape(median(C,1),A.n_cols,B.n_cols);
}

arma::colvec computeBetaJ_cpp(arma::colvec y, arma::mat X, int j, int m, int adjFinSample) {
  int bigT = X.n_rows;
  int numThrow = ceil((double)m / 2);
  if (adjFinSample == 0) numThrow = m ;
  mat XSparse; colvec ySparse;
  XSparse.zeros(bigT - numThrow, X.n_cols); ySparse.zeros(bigT - numThrow);
  if (j < 2) {
    XSparse.rows(0, bigT - numThrow - 1) = X.rows(numThrow, bigT - 1);
    ySparse.rows(0, bigT - numThrow - 1) = y.rows(numThrow, bigT - 1);
  } else {
    XSparse.rows(0, j - 2) = X.rows(0, j - 2); ySparse.rows(0, j - 2) = y.rows(0, j - 2);
    XSparse.rows(j - 1, bigT - numThrow - 1) = X.rows(j + numThrow - 1, bigT - 1);
    ySparse.rows(j - 1, bigT - numThrow - 1) = y.rows(j + numThrow - 1, bigT - 1);
  }
  colvec beta = runOLSFast_cpp(ySparse, XSparse);
  return(beta);
}

void
HMM::computeXiCached() {
  arma::mat temp = B_.rows(1,T_-1) % beta_.cols(1,T_-1).t();
  for(unsigned int i = 0; i < N_; ++i) {
    xi_.slice(i) = temp % 
      (alpha_(i,arma::span(0, T_-2)).t() * A_.row(i));
  }
}

static arma::vec s_recall(arma::mat &runMatrix, Qrels &qrels, int rank){
    // Use a greedy approach to find the max subtopics at rank k

    // Initialize working variables
    int numOfSubtopics = int(*qrels.subtopics.rbegin());
    int numOfRows = min(rank, qrels.numOfRelDocs);
    arma::vec maxSubtopics = arma::zeros(rank);
    arma::rowvec seenSubtopics = arma::zeros<arma::rowvec>(numOfSubtopics);
    arma::vec tmpVector;
    set<int> idealRankList;
    set<int>::iterator it;

    // Greedy approach searches for the best document at each rank
    for(int i = 0; i < numOfRows; i++){
        int maxNewSubtopics = 0;
        int pickedDoc = -1;
        // Iterate through the set of relevant documents to find
        //   the best document
        for(int j = 0; j < qrels.matrix.n_rows; j++){

            // Ignore the document already picked
            it = idealRankList.find(j);
            if(it != idealRankList.end() )
                continue;

            // Compute the number of new subtopics the document contains
            int numOfNewSubtopics = arma::sum((qrels.matrix.row(j) + seenSubtopics) >
                                arma::zeros<arma::vec>(numOfSubtopics));

            // Keep track of the best document/ max number of subtopics
            if(numOfNewSubtopics > maxNewSubtopics){
                maxNewSubtopics = numOfNewSubtopics;
                pickedDoc = j;
            }
        }

        // Add the best document to the ideal rank list
        //  and keep track of subtopics seen
        seenSubtopics += qrels.matrix.row(pickedDoc);
        idealRankList.insert(pickedDoc);
        maxSubtopics(i) = arma::sum(seenSubtopics >
                                        arma::zeros<arma::vec>(numOfSubtopics));
    }
    if(numOfRows < rank){
        for(int i=numOfRows; i < rank; i++)
            maxSubtopics(i) = maxSubtopics(numOfRows - 1);
    }

    // Consider only the top 'rank' number of documents
    arma::vec s_recall = arma::zeros(rank);
    for(int i=0; i < rank; i++){
        arma::mat matrix = runMatrix.rows(0,i);
        double retSubtopics = arma::sum(arma::sum(matrix) >
                                        arma::zeros<arma::vec>(numOfSubtopics));
        s_recall(i) = (retSubtopics/maxSubtopics(i));
    }
    return s_recall;
}

//[[Rcpp::export]]
double computeSStat_cpp(arma::colvec beta, arma::mat invS, arma::colvec y, arma::mat X,
                        int j, int m) {
  mat newX = X.rows(j - 1, j + m - 2); colvec newY = y.rows(j - 1, j + m - 2);
  mat A = trans(newX) * invS * (newY - newX * beta);
  mat V =  trans(newX) * invS * newX;
  mat S = trans(A) * inv(V) * A;
  double s = S(0, 0);
  return(s);
}

/*******************************************************************
 * 函数功能：几何校正需要调用的函数
 * 
 *
 */
arma::mat toAffinity(arma::mat &f){
	arma::mat A;
	arma::mat v;
	v << 0 << 0 << 1 << arma::endr;
	int flag = f.n_rows;
	switch(flag){
		case 6:{ // oriented ellipses
			arma::mat T = f.rows(0, 1);
			arma::mat tmp = join_horiz(f.rows(2, 3), f.rows(4, 5));
			arma::mat tmp1 = join_horiz(tmp, T);
			A = join_vert(tmp1, v);
			break;}
		case 4:{   // oriented discs
			arma::mat T = f.rows(0, 1);
			double s = f.at(2,0);
			double th = f.at(3,0);
			arma::mat S = arma::randu<arma::mat>(2,2);
			/*S.at(0, 0) = s*cos(th);
			S.at(0, 1) = -s*sin(th);
			S.at(1, 0) = s*sin(th);
			S.at(1, 1) = s*cos(th);*/
			S << s*cos(th) << -s*sin(th) << arma::endr
			  << s*sin(th) << s*cos(th)  << arma::endr;
			arma::mat tmp1 = join_horiz(S, T);
			A = join_vert(tmp1, v);
			//A.print("A =");
			break;}
		/*case 3:{    // discs
			mat T = f.rows(0, 1);
			mat s = f.row(2);
			int th = 0 ;
			A = [s*[cos(th) -sin(th) ; sin(th) cos(th)], T ; 0 0 1] ;
			   }
		case 5:{ // ellipses
			mat T = f.rows(0, 1);
			A = [mapFromS(f(3:5)), T ; 0 0 1] ;
			   }*/
		default:
            std::cout << "出错啦！" << std::endl;
			break;
	}
	return A;
}

// **********************************************************//
//                  Calculate Triadic statistics              //
// **********************************************************//
// [[Rcpp::export]]
arma::cube triadicstat (arma::cube cache, arma::mat intervals_d, int a, int A, arma::mat pd) {
	int nIP = pd.n_cols;
	arma::cube out = arma::zeros<arma::cube>(A, nIP, 36);
	int it = 0;
	for (unsigned int i = 0; i < 3; i++) {
		for (unsigned int j = 0; j < 3; j++) {
		int mins1 = intervals_d(i, 0)-1;
		int maxs1 = intervals_d(i, 1)-1;
		int mins2 = intervals_d(j, 0)-1;
		int maxs2 = intervals_d(j, 1)-1;
		arma::mat tsend = arma::zeros<arma::mat>(A, nIP);
		arma::mat treceive = arma::zeros<arma::mat>(A, nIP);
		arma::mat sibling = arma::zeros<arma::mat>(A, nIP);
		arma::mat cosibling = arma::zeros<arma::mat>(A, nIP);
		arma::mat pdnew1 = pd.rows(mins1, maxs1);
		arma::mat pdnew2 = pd.rows(mins2, maxs2);
		for (int r = 0; r < A; r++) {
			for (int h = 0; h < A; h++) {
				arma::colvec docsah = cache(arma::span(a), arma::span(h), arma::span(mins1, maxs1));
				arma::colvec docsrh = cache(arma::span(r), arma::span(h), arma::span(mins2, maxs2));
				arma::colvec docsha = cache(arma::span(h), arma::span(a), arma::span(mins1, maxs1));
				arma::colvec docshr = cache(arma::span(h), arma::span(r), arma::span(mins2, maxs2));
				tsend.row(r) += (docsah.t() * pdnew1) % (docshr.t() * pdnew2);				
				treceive.row(r) += (docsha.t() * pdnew1) % (docsrh.t() * pdnew2);	
				sibling.row(r) += (docsha.t() * pdnew1) % (docshr.t() * pdnew2);				
				cosibling.row(r) += (docsah.t() * pdnew1) % (docsrh.t() * pdnew2);	
				}	
		}
		for (int IP = 0; IP < nIP; IP++) {
			for (int r = 0; r < A; r++) {
				out(r, IP, it) = tsend(r, IP);
				out(r, IP, it+9) = treceive(r,IP);
				out(r, IP, it+18) = sibling(r, IP);
				out(r, IP, it+27) = cosibling(r,IP);
			}
		}
		it +=1;		
	}
	}
	arma::cube out2 = arma::zeros<arma::cube>(A, nIP, 12);
	for (int t = 0; t < 4; t++) {
		int add = 9*t;
	for (int IP = 0; IP < nIP; IP++) {
		for (int r = 0; r < A; r++) {
			out2(r, IP, 3*t+0) = out(r, IP, add+0);
			out2(r, IP, 3*t+1) = out(r, IP, add+1) + out(r, IP, add+3)+ out(r, IP, add+4);
			out2(r, IP, 3*t+2) = out(r, IP, add+2) + out(r, IP, add+5)+ out(r, IP, add+6)+ out(r, IP, add+7)+ out(r, IP, add+8);
			}
		}
	}		
	return out2 / 10;
}

// This piece of code computes the flowlength before taking into account 
// the dimensions of pixels/slope. 
//[[Rcpp::export]]
double fl_internal(arma::mat m) { 
  
  uword ny = m.n_cols; 
  uword nx = m.n_rows; 
  
  rowvec flcol = zeros(1, ny); 
  for ( uword r=0; r<nx; r++ ) { 
    rowvec a = sum(cumprod(m.rows(r, nx - 1), 0), 0); 
    flcol += a; 
  } 
  
  return( accu(flcol)/((double)nx * (double)ny) ); 
}

void RNN<
LayerTypes, OutputLayerType, InitializationRuleType, PerformanceFunction
>::Gradient(const arma::mat& /* unused */,
            const size_t i,
            arma::mat& gradient)
{
  if (gradient.is_empty())
  {
    gradient = arma::zeros<arma::mat>(parameter.n_rows, parameter.n_cols);
  }
  else
  {
    gradient.zeros();
  }

  Evaluate(parameter, i, false);

  arma::mat currentGradient = arma::mat(gradient.n_rows, gradient.n_cols);
  NetworkGradients(currentGradient, network);

  const arma::mat input = arma::mat(predictors.colptr(i), predictors.n_rows,
      1, false, true);

  // Iterate through the input sequence and perform the feed backward pass.
  for (seqNum = seqLen - 1; seqNum >= 0; seqNum--)
  {
    // Load the network activation for the upcoming backward pass.
    LoadActivations(input.rows(seqNum * inputSize, (seqNum + 1) *
        inputSize - 1), network);

    // Perform the backward pass.
    if (seqOutput)
    {
      arma::mat seqError = error.unsafe_col(seqNum);
      Backward(seqError, network);
    }
    else
    {
      Backward(error, network);
    }

    // Link the parameters and update the gradients.
    LinkParameter(network);
    UpdateGradients<>(network);

    // Update the overall gradient.
    gradient += currentGradient;

    if (seqNum == 0) break;
  }
}

arma::colvec computeTestStatDistribution_cpp(arma::colvec y, arma::mat X, int breakDate,
                                             int m, int adjFinSample, arma::mat invS) {
  int bigT = X.n_rows;
  colvec sJ = zeros(breakDate - m + 1);
  for (int j = 0; j < (breakDate - m + 1); j++) {
    colvec betaJ = computeBetaJ_cpp(y.rows(0, bigT - 1 - m), X.rows(0, bigT - 1 - m),
                                    j + 1, m, adjFinSample);
    sJ(j) = computeSStat_cpp(betaJ, invS, y, X, j + 1, m);
//     if (j == 44) {
//       Rcout << "betaJ = "<< betaJ << "sJ(j) = " << sJ(j) << endl;
//     }
  }
  return(sJ);
}