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

//' Fast addition of vector to each row of matrix
//'
//' @description Fast addition of vector to each row of a matrix. This corresponds to t(t(x) + v)
//' @param x A matrix with dimensions n*k.
//' @param v A vector of length k.
//' @return A matrix of dimension n*k where v is added to each row of x
//' @author Claus Ekstrom <[email protected]@rprimer.dk>
//' @examples
//'
//' A <- matrix(1:12, ncol=3)
//' B <- c(1, 2, 3)
//'
//' add_torows(A, B)
//'
//' @export
// [[Rcpp::export]]
arma::mat add_torows(const arma::mat &x, const arma::rowvec &v) {

  if (x.n_cols != v.n_elem) {
     stop("Non-conformable objects. Length of v does not match columns in x");
  }

  arma::mat m = x.each_row() + v;
  
  return (m);
}

arma::cx_mat generate_visibilities_from_local_skymodel(arma::mat& skymodel, arma::mat& uvw_baselines)
{
    arma::cx_mat model_vis(uvw_baselines.n_rows, 1);
    model_vis.zeros();

    skymodel.each_row([&](arma::rowvec& src_entry) {
        model_vis += visibilities_for_point_source(uvw_baselines, src_entry[0], src_entry[1], src_entry[2]);
    });

    return std::move(model_vis);
}

unsigned int optCoef(arma::mat& weights, const arma::icube& obs, const arma::cube& emission,
    const arma::mat& initk, const arma::cube& beta, const arma::mat& scales, arma::mat& coef,
    const arma::mat& X, const arma::ivec& cumsumstate, const arma::ivec& numberOfStates,
    int trace) {

  int iter = 0;
  double change = 1.0;
  while ((change > 1e-10) & (iter < 100)) {
    arma::vec tmpvec(X.n_cols * (weights.n_rows - 1));
    bool solve_ok = arma::solve(tmpvec, hCoef(weights, X),
        gCoef(obs, beta, scales, emission, initk, weights, X, cumsumstate, numberOfStates));
    if (solve_ok == false) {
      return (4);
    }

    arma::mat coefnew(coef.n_rows, coef.n_cols - 1);
    for (unsigned int i = 0; i < (weights.n_rows - 1); i++) {
      coefnew.col(i) = coef.col(i + 1) - tmpvec.subvec(i * X.n_cols, (i + 1) * X.n_cols - 1);
    }
    change = arma::accu(arma::abs(coef.submat(0, 1, coef.n_rows - 1, coef.n_cols - 1) - coefnew))
        / coefnew.n_elem;
    coef.submat(0, 1, coef.n_rows - 1, coef.n_cols - 1) = coefnew;
    iter++;
    if (trace == 3) {
      Rcout << "coefficient optimization iter: " << iter;
      Rcout << " new coefficients: " << std::endl << coefnew << std::endl;
      Rcout << " relative change: " << change << std::endl;
    }
    weights = exp(X * coef).t();
    if (!weights.is_finite()) {
      return (5);
    }
    weights.each_row() /= sum(weights, 0);

  }
  return (0);
}

///
/// \brief Vespucci::MathDimensionReduction::plsregress PLS Regression Using SIMPLS algorithm.
/// This is essentially a line-for-line rewrite of plsregress from the Octave
/// statistics package.
/// Copyright (C) 2012 Fernando Damian Nieuwveldt <[email protected]>
/// This is an implementation of the SIMPLS algorithm:
/// Reference: SIMPLS: An alternative approach to partial least squares regression.
/// Chemometrics and Intelligent Laboratory Systems (1993)
/// \param x
/// \param y
/// \param components Number of PLS components
/// \param x_loadings "Output" values
/// \param y_loadings
/// \param x_scores
/// \param y_scores
/// \param coefficients
/// \param fitted_data
/// \return
///
bool Vespucci::Math::DimensionReduction::plsregress(arma::mat X, arma::mat Y,
                                                    int components, arma::mat &X_loadings,
                                                    arma::mat &Y_loadings,
                                                    arma::mat &X_scores,
                                                    arma::mat &Y_scores,
                                                    arma::mat &coefficients,
                                                    arma::mat &percent_variance,
                                                    arma::mat &fitted)
{
    using namespace arma;
    uword observations = X.n_rows;
    uword predictors = X.n_cols;
    uword responses = Y.n_cols;
    //Mean centering
    mat Xmeans = arma::mean(X);
    mat Ymeans = arma::mean(Y);
    //Octave does below with bsxfun. After optimization, this should hopefully not
    // be slower.
    X.each_row() -= Xmeans; //element-wise subtraction of mean
    Y.each_row() -= Ymeans; //same

    mat S = trans(X) * Y;
    mat R = zeros(predictors, components);
    mat P = R;
    mat V = R;
    mat T = zeros(observations, components);
    mat U = T;
    mat Q = zeros(responses, components);
    mat eigvec;
    vec eigval;

    mat q;
    mat r;
    mat t;
    mat nt;
    mat p;
    mat u;
    mat v;
    double max_eigval;
    for (int i = 0; i < components; ++i){
        eig_sym(eigval, eigvec, (trans(S) * S));
        max_eigval = eigval.max();
        uvec dom_index = find(eigval >= max_eigval);
        uword dominant_index = dom_index(0);

        q = eigvec.col(dominant_index);

        r = S*q; //X block factor weights
        t = X*r; //X block factor weights
        t.each_row() -= mean(t); //center t
        nt = arma::sqrt(t.t()*t); //compute norm (is arma::norm() the same?)
        t.each_row() /= nt;
        r.each_row() /= nt; //normalize

        p = X.t()*t; //X block factor loadings
        q = Y.t()*t; //Y block factor loadings
        u = Y*q; //Y block factor scores
        v = p;

        //Ensure orthogonality
        if (i > 0){
            v = v - V*(V.t()*p);
            u = u - T*(T.t()*u);
        }
        v.each_row() /= arma::sqrt(trans(v) * v); //normalize orthogonal loadings
        S = S - v * (trans(v)*S); //deflate S wrt loadings
        R.col(i) = r;
        T.col(i) = t;
        P.col(i) = p;
        Q.col(i) = q;
        U.col(i) = u;
        V.col(i) = v;
    }

    //Regression coefficients
    mat B = R*trans(Q);
    fitted = T*trans(Q);
    fitted.each_row() += Ymeans;

    //Octave creates copies from inputs before sending to output. Doing same
    //here just to be safe.
//.........这里部分代码省略.........

//[[Rcpp::export(.emplik_intern)]]
List emplik_intern(arma::mat z, arma::colvec mu, arma::vec lam, double eps,
	double M = 1e30, double thresh = 1e-12, int itermax = 1000){
//(arma::mat z, arma::vec mu  = vec(z.n_cols,fill::zeros), double eps = 1/z.nrows, double M = datum::inf);
// # Backtracking line search parameters [Tweak only with extreme caution.]
// # See Boyd and Vandenberghe, pp 464-466.
double ALPHA = 0.3; // seems better than 0.01 on some 2d test data (sometimes fewer iters)
double BETA  = 0.8;
// # We need  0 < ALPHA < 1/2   and 0 < BETA < 1
// # Backtrack threshold: you can miss by this much.
double BACKEPS = 0;
// # Consider replacing 0 by 1e-10 if backtracking seems to be
// # failing due to round off.
int n = z.n_rows;
int d = z.n_cols;
 z.each_row() -= trans(mu);
//  Use lam = 0 or initial lam, whichever is best
arma::vec onen = arma::vec(n,arma::fill::ones);
arma::mat init0 = mllog(onen, eps, M, 2);
arma::mat init(init0.n_rows, init0.n_cols);
if(!any(lam)){
  init = init0;
} else {
  init = mllog(onen+z*lam, eps, M, 2);
	if(sum(init0.col(0) < sum(init.col(0)))){
    lam = arma::rowvec(z.n_cols,arma::fill::zeros);
    init = init0;
  }
}
// # Initial f, g
double fold = sum(init.col(0));
arma::rowvec gold = sum(z.each_col() % init.col(1),0);
bool converged = false;
int iter = 0;
arma::mat oldvals = init;
arma::mat newvals(oldvals.n_rows, oldvals.n_cols);
arma::vec rootllpp;
arma::mat zt(z.n_rows, z.n_cols);
arma::vec yt(z.n_rows);
arma::vec wts;
double fnew; double targ;
double logelr;
arma::vec step(z.n_cols);
int s; double ndec; double gradnorm;
bool backtrack = false;
while(!converged){
  iter += 1;
  //   # Get Newton Step
  rootllpp = sqrt(oldvals.col(2));  //# sqrt 2nd deriv of -llog lik
  zt = z;
     for(int j=0; j<d; j++){
      zt.col(j) %= rootllpp;
    }
    yt   = oldvals.col(1) / rootllpp;
    step = -svdlm(zt,yt);
    backtrack = false;
    s = 1;
    while(!backtrack){
         newvals = mllog(onen + z * (lam+s*step), eps, M, 2);
        fnew = sum(newvals.col(0));
        targ = fold + ALPHA * s * sum(trans(gold) % step) + BACKEPS; // (BACKEPS for roundoff, should not be needed)
         if(fnew <= targ){ // backtracking has converged
          backtrack = true;
          oldvals = newvals;
          fold = fnew;
     			gold = sum(z.each_col() % oldvals.col(1),0);
          lam = lam + s*step;
         } else{
          s = s * BETA;
         }
       }
    //   # Newton decrement and gradient norm
      ndec     = sqrt(sum(square(step % trans(gold))));
      gradnorm = sqrt(sum(square(gold)));
      converged = (ndec * ndec <= thresh);
       if( iter > itermax )
       	break;
     }
     wts = pow(1 + z * lam, -1)/n;
		logelr = sum(mllog(onen + z * lam, eps, M, 0).col(0));
    return Rcpp::List::create(Named("logelr")=logelr, Named("lam") = lam, Named("wts") = wts,
        Named("conv") = converged, Named("niter") = iter,Named("ndec") = ndec, Named("gradnorm") = gradnorm);
     }

//.........这里部分代码省略.........

        // iterate over the time sequence, solving the LNA over each interval
        for(int j=0; j < (n_times-1); ++j) {

                // set the times of the interval endpoints
                t_L = lna_times[j];
                t_R = lna_times[j+1];

                // Reset the LNA state vector and integrate the LNA ODEs over the next interval to 0
                std::fill(lna_state_vec.begin(), lna_state_vec.end(), 0.0);
                CALL_INTEGRATE_STEM_ODE(lna_state_vec, t_L, t_R, step_size, lna_pointer);

                // transfer the elements of the lna_state_vec to the process objects
                std::copy(lna_state_vec.begin(), lna_state_vec.begin() + n_events, lna_drift.begin());
                std::copy(lna_state_vec.begin() + n_events, lna_state_vec.end(), lna_diffusion.begin());

                // ensure symmetry of the diffusion matrix
                lna_diffusion = arma::symmatu(lna_diffusion);

                // map the stochastic perturbation to the LNA path on its natural scale
                try{
                        if(lna_drift.has_nan() || lna_diffusion.has_nan()) {
                                good_svd = false;
                                throw std::runtime_error("Integration failed.");
                        } else {
                                good_svd = arma::svd(svd_U, svd_d, svd_V, lna_diffusion); // compute the SVD
                        }

                        if(!good_svd) {
                                throw std::runtime_error("SVD failed.");

                        } else {
                                svd_d.elem(arma::find(svd_d < 0)).zeros();          // zero out negative sing. vals
                                svd_V.each_row() %= arma::sqrt(svd_d).t();          // multiply rows of V by sqrt(d)
                                svd_U *= svd_V.t();                                 // complete svd_sqrt
                                svd_U.elem(arma::find(lna_diffusion == 0)).zeros(); // zero out numerical errors

                                log_lna = lna_drift + svd_U * draws.col(j);         // map the LNA draws
                        }

                } catch(std::exception & err) {

                        // reinstatiate the SVD objects
                        arma::vec svd_d(n_events, arma::fill::zeros);
                        arma::mat svd_U(n_events, n_events, arma::fill::zeros);
                        arma::mat svd_V(n_events, n_events, arma::fill::zeros);

                        // forward the exception
                        forward_exception_to_r(err);

                } catch(...) {
                        ::Rf_error("c++ exception (unknown reason)");
                }

                // compute the LNA increment
                nat_lna = arma::vec(expm1(Rcpp::NumericVector(log_lna.begin(), log_lna.end())));

                // save the LNA increment
                pathmat(j+1, arma::span(1, n_events)) = nat_lna.t();

                // update the initial volumes
                init_volumes += stoich_matrix * nat_lna;

                // if any increments or volumes are negative, throw an error
                try{
                        if(any(nat_lna < 0)) {

// [[Rcpp::export]]
List solveHierBasis(arma::mat design_mat,
                  arma::vec y,
                  arma::vec ak,arma::mat weights,
                  int n, double lam_min_ratio, int nlam,
                  double max_lambda) {

  // This function returns the argmin of the following function:
  // (0.5) * ||y - X %*% beta||_2^2 + P(beta), where
  // P(\beta) = \sum_{j=1}^{J_n} weights[j, lam_indx] * norm(beta[j:J_n])
  // where j indexes the basis function and lam_indx indexes the different
  // lambda values.
  //
  // Args:
  //    design_mat: A design matrix of size n * J_n, where J_n is number of
  //                basis functions. This matrix should be centered.
  //    y: The centered response vector.
  //    ak: A J_n-vector of weights where ak[j] = j^m - (j-1)^m for a smoothness
  //        level m.
  //    weights: A J_n * nlam matrix which is simply the concatenation [ak,ak,...,ak].
  //    n: The number of observations, equal to length(y).
  //    lam_min_ratio: Ratio of min_lambda to max_lambda.
  //    nlam: Number of lambda values.
  //    max_lambda: A double specifying the maximum lambda, if NA then function
  //                selects a max_lambda.
  // Returns:
  //    beta_hat2: The solution matrix, a J_n * nlam  matrix.
  //    lambdas: Sequence of lambda values used for fitting.

  // Generate the x_mat, this is our orthonormal design matrix.
  arma::mat x_mat, R_mat;
  arma::qr_econ(x_mat, R_mat, design_mat);

  x_mat = x_mat * sqrt(n);
  R_mat  = R_mat / sqrt(n);
  // arma::sp_mat R_mat2 = sp_mat(R_mat / sqrt(n));


  // Note that for a orthogonal design with X^T %*% X = nI the
  // two problems are equivalent:
  // 1. (1/2n)*|Y - X %*% beta|^2 + lambda * Pen(beta),
  // 2. (1/2)*|t(X) %*% Y/n - beta|^2 + lambda * Pen(beta).
  //
  // Thus all we need, is to solve the proximal problem.
  // Define v_temp = t(X) %*% Y/n for the prox problem.

  arma::vec v_temp = x_mat.t() * (y /n);

  // We evaluate a max_lambda value if not specified by user.
  if(R_IsNA(max_lambda)) {
    // Followed by the maximum lambda value.
     max_lambda = max(abs(v_temp) / ak);
    // max_lambda = norm(v_temp, 2);
  }

  // Generate the full lambda sequence.
  arma::vec lambdas = linspace<vec>(log10(max_lambda),
                                    log10(max_lambda * lam_min_ratio),
                                    nlam);
  lambdas = exp10(lambdas);

  // Generate matrix of weights.
  weights.each_row() %= lambdas.t();

  // Obtain beta.hat, solve prox problem.
  arma::sp_mat beta_hat =  GetProx(v_temp, weights);

  // Take the inverse of the R_mat to obtain solution on the original scale
  arma::mat beta_hat2 = solve(trimatu(R_mat), mat(beta_hat));

  ///////////// Removed after version 0.7.1, moved to a different R function.//////////

  // We also obtain the DOF for the different methods.
  //arma::vec dof = getDof(x_mat, weights, beta_hat, nlam, n);

  //////////// Removed after version 0.7.1, moved to a different R function.//////////


  // Return fitted values.
  arma::mat yhat = x_mat * mat(beta_hat);

  return List::create(Named("beta") = beta_hat2,
                      Named("lambdas") = lambdas,
                      //Named("dof") = dof,
                      Named("yhat") = yhat,
                      Named(".beta_ortho") = beta_hat);
}

// [[Rcpp::export]]
List solveHierLogistic(arma::mat design_mat,
                       arma::vec y,
                       arma::vec ak,arma::mat weights,
                       int n, int nlam, int J,
                       double max_lambda, double lam_min_ratio,
                       double tol, int max_iter,
                       double tol_inner, int max_iter_inner) {

  // This function returns the argmin of the following function:
  // -L(beta) + P(beta), where
  // P(\beta) = \sum_{j=1}^{J_n} weights[j, lam_indx] * norm(beta[j:J_n])
  // where j indexes the basis function and lam_indx indexes the different
  // lambda values. 'L' in this case is the log-likelihood for a binomial model.
  //
  // Args:
  //    design_mat: A design matrix of size n * J_n, where J_n is number of
  //                basis functions.
  //    y: The response vector.
  //    ak: A J_n-vector of weights where ak[j] = j^m - (j-1)^m for a smoothness
  //        level m.
  //    weights: A J_n * nlam matrix which is simply the concatenation [ak,ak,...,ak].
  //    n: The number of observations, equal to length(y).
  //    lam_min_ratio: Ratio of min_lambda to max_lambda.
  //    nlam: Number of lambda values.
  //    max_lambda: A double specifying the maximum lambda. A max lambda needs to
  //                be specified for logistic regression.
  //    tol, tol_inner: The tolerance for the outer and inner loop respectively.
  //    max_iter, max_iter_inner: The maximum number of iterations for outer and
  //                              inner loops respectively.
  // Returns:
  //    beta_final: The sparse solution matrix, a J_n * nlam  matrix.
  //    intercept: The vector of fitted intercepts of size nlam.
  //    lambdas: Sequence of lambda values used for fitting.
  //    fitted: The fitted probabilities.
  //

  // Generate the x_mat, this is our orthonormal design matrix.
  arma::mat x_mat, R_mat;
  arma::qr_econ(x_mat, R_mat, design_mat);

  x_mat = x_mat * sqrt(n);
  R_mat  = R_mat / sqrt(n);

  // We evaluate a max_lambda value if not specified by user.
  if(R_IsNA(max_lambda)) {

    // A simple calculation to evelaute the max lambda,
    // If we wish to find a max lambda, then beta = 0 and hence the
    // temp fitted probabilities are 0.5. The current quadratic approximation
    // is given by
    // y_tilde =  intercept + X * beta + (y - p_hat)/(p_jhat * (1 - p_hat)).
    // In this case this is simply 4 * (y - 0.5).
    // Finally max_lambda is give  first finding
    // v_temp = t(x_mat) * (4 * (y - 0.5))/n followed by
    // max(abs(v_temp)) / (4*ak).
    arma::vec v_temp = x_mat.t() * ((y - 0.5) / n);

    // Followed by the maximum lambda value.
    max_lambda = max(abs(v_temp) / ak);
  }


  // Generate the full lambda sequence.
  arma::vec lambdas = linspace<vec>(log10(max_lambda),
                                    log10(max_lambda * lam_min_ratio),
                                    nlam);
  lambdas = exp10(lambdas);

  // Generate matrix of weights.
  weights.each_row() %= lambdas.t();

  // The final matrix we will use to store beta values.
  arma::mat beta_ans(J, nlam);
  // The final vector to store values for the intercept.
  arma::vec intercept_ans(nlam);

  // The vectors we will use to check convergence.
  // We also use this for the sake of warm starts.
  arma::vec beta(J, fill::zeros);
  //  arma::vec beta_new(J, fill::zeros);


  double intercept = 0;
  //  double intercept_new = 0;

  // The full parameter vector which containts the intercept + covariates.
  arma::vec pars(J + 1, fill::zeros);
  arma::vec pars_new(J + 1, fill::zeros);

  // For each lambda and each middle iteration number,
  // We design a new vector of probabilities and consequently a new
  // response vector.
  arma::vec temp_prob;
  arma::vec temp_resp;

  // We also store a temporary linear part, given by X * beta.
  arma::vec temp_xbeta;

  // Begin outer loop to decrement lambdas.
//.........这里部分代码省略.........