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

//' Parallel C++ implementation of multivariate Normal probability density function
//' 
//'Based on the implementation from Nino Hardt and Dicko Ahmadou
//'http://gallery.rcpp.org/articles/dmvnorm_arma/ 
//'(accessed in August 2014)
//'
//'@param x data matrix
//'@param mean mean vector
//'@param varcovM variance covariance matrix
//'@param logical flag for returning the log of the probability density 
//'function. Defaults is \code{TRUE}
//'@return vector of densities
//'
//'@export
//'
//'@examples
//'mvnpdf(x=matrix(1.96), mean=0, varcovM=diag(1), Log=FALSE)
//'mvnpdfC(x=matrix(1.96), mean=0, varcovM=diag(1), Log=FALSE)
//'mvnpdf(x=matrix(1.96), mean=0, varcovM=diag(1))
//'mvnpdfC(x=matrix(1.96), mean=0, varcovM=diag(1))
//'
//'library(microbenchmark)
//'microbenchmark(mvnpdf(x=matrix(1.96), mean=0, varcovM=diag(1), Log=FALSE),
//'               mvnpdfC(x=matrix(1.96), mean=0, varcovM=diag(1), Log=FALSE),
//'               times=10000L)     
//'               
// [[Rcpp::export]]
NumericVector mvnpdfC_par(NumericMatrix x, 
                      NumericVector mean, 
                      NumericMatrix varcovM,
                      bool Log=true,
                      int ncores=1){

    const mat xx = as<mat>(x);
    const mat S =  as<mat>(varcovM);
    const colvec m = as<colvec>(mean); 
    const int p = xx.n_rows;
    const int n = xx.n_cols;
    NumericVector y = NumericVector(n);
    
    const mat Rinv = inv(trimatu(chol(S)));
    const double logSqrtDetvarcovM = sum(log(Rinv.diag()));
    const double constant = - p*log2pi2;
    
    omp_set_num_threads(ncores);
    
    #pragma omp parallel for shared(y)
    for (int i=0; i < n; i++) {
        colvec x_i = xx.col(i) - m;
        rowvec xRinv = trans(x_i)*Rinv;
        double quadform = sum(xRinv%xRinv);
        y(i) = exp(-0.5*quadform + logSqrtDetvarcovM + constant);
        if (!Log) {
            y(i) = exp(-0.5*quadform + logSqrtDetvarcovM + constant);
        } else{
            y(i) = -0.5*quadform + logSqrtDetvarcovM + constant;
        }
    }
    
    return y;

}

void JacobiRotations(int n, double rho_max, mat &A, uvec diagonalA, mat &R){

    double epsilon = pow(10,-8);
    //cout << epsilon << endl;

    // Starting w/A(1,2)
    //cout << A << endl;
    int k = 0; int l = 0;
    double max_A = 0.0;
    findMaximumElementOnNonDiagonal(A, k, l, max_A, n);
    //cout << max_A << endl;

    int numberOfIterations = 0;
    int maxIterations = pow(10, 6);

    //==================================================================================
    //While-loop to change the non-diagonal elements to be approximately zero.
    //==================================================================================


    while( fabs(max_A) > epsilon && (double) numberOfIterations < maxIterations){
        double c = 0.0;
        double s = 0.0;
        findSinAndCos(A, k, l, c, s);
        // Rotate
        rotateMatrix(A, k, l, c, s, n, R);
        //Find max off-diagonal matrix element
        findMaximumElementOnNonDiagonal(A, k, l, max_A, n);
        numberOfIterations++;
    }

    //Preparing results
    diagonalA = sort_index(A.diag());

    //Printing

    //printMatlabMatrix("A", A);
    cout << "Number of iterations: " << numberOfIterations << endl;
    //cout << rho_max << endl;
    cout << "lambda0:  " << A.diag()(diagonalA(0)) << "  lambda1:  " << A.diag()(diagonalA(1)) << "  lambda2:  " << A.diag()(diagonalA(2)) << endl;

    //cout << "hallo" << endl;


}

mat nnls_solver_without_missing(mat & WtW, mat & WtA,
	const mat & A, const mat & W, const mat & W1, const mat & H2, const umat & mask, 
	const double & eta, const double & beta, int max_iter, double rel_tol, int n_threads)
{

	// A = [W, W1, W2] [H, H1, H2]^T.
	// Where A has not missing
	// WtW and WtA are auxiliary matrices, passed by referenced and can be modified 
	update_WtW(WtW, W, W1, H2);
	update_WtA(WtA, W, W1, H2, A);
	if (beta > 0) WtW += beta;
	if (eta > 0) WtW.diag() += eta;
	
	return nnls_solver(WtW, WtA, mask, max_iter, rel_tol, n_threads);
}

void fill_matrix(mat &A, vec pot, double h)
{
  A.diag() = (2/(h*h) + pot);
  A.diag(1).fill(-1/(h*h));
  A.diag(-1).fill(-1/(h*h));
}