Skip to content

Commit

Permalink
Fixed warnings related to updateRP_genomic_cpp.h
Browse files Browse the repository at this point in the history
  • Loading branch information
@wleoncio
wleoncio committed Oct 22, 2024
1 parent 2ce0cd7 commit 83698c9
Show file tree
Hide file tree
Showing 2 changed files with 130 additions and 59 deletions.
92 changes: 71 additions & 21 deletions src/updateRP_genomic_cpp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,10 +3,60 @@

#include "updateRP_genomic_cpp.h"

double be_prop_me_ini = 0.;
double be_prop_sd_ini = 0.;
double D1 = 0.;
double D2 = 0.;
double D1_prop = 0.;
double D2_prop = 0.;
double loglh_ini = 0.;
double loglh_prop = 0.;
double logprior_prop = 0.;
double logprior_ini = 0.;
double logprop_prop = 0.;
double logprop_ini = 0.;
double logR = 0.;

double be_prop_me = 1.;
double be_prop_sd = 1.;

arma::vec be_prop;
arma::vec exp_xbeta;
arma::mat h_exp_xbeta_mat;
arma::mat h_exp_xbeta_prop_mat;
arma::mat exp_h_exp_xbeta_mat;
arma::mat exp_h_exp_xbeta_prop_mat;
arma::vec first_sum;
arma::vec second_sum;
arma::vec first_sum_prop;
arma::vec second_sum_prop;
arma::vec x_exp_xbeta;
arma::vec xbeta_prop;
arma::vec exp_xbeta_prop;
arma::vec x_exp_xbeta_prop;
arma::vec x_sq_exp_xbeta;
arma::vec x_sq_exp_xbeta_prop;
arma::vec D1_1st;
arma::vec D1_2nd;
arma::vec D1_1st_prop;
arma::vec D1_2nd_prop;
arma::vec D2_1st;
arma::vec D2_2nd;
arma::vec D2_1st_prop;
arma::vec D2_2nd_prop;
arma::mat D1_2nd_den;
arma::mat D1_2nd_num;
arma::mat D1_2nd_den_prop;
arma::mat D1_2nd_num_prop;
arma::mat D2_2nd_num;
arma::mat D2_2nd_den;
arma::mat D2_2nd_den_prop;
arma::mat D2_2nd_num_prop;

arma::mat matProdVec(const arma::mat x, const arma::vec y)
{
// multiply (element-wise) a matrix to a expanded vector

arma::mat mat_y = arma::zeros<arma::mat>(y.n_elem, x.n_cols);
mat_y.each_col() = y;
arma::mat spanMat = x % mat_y; // elementwise product
Expand All @@ -16,7 +66,7 @@ arma::mat matProdVec(const arma::mat x, const arma::vec y)
arma::vec sumMatProdVec(const arma::mat x, const arma::vec y)
{
// compute "arma::sum( matProdVec( ind_r_d_, exp_xbeta ).t(), 1 );"

arma::vec spanVec = arma::zeros(x.n_cols);
for (unsigned int i = 0; i < x.n_cols; ++i)
spanVec(i) = arma::dot(x.col(i), y);
Expand All @@ -34,7 +84,7 @@ arma::vec updateBH_cpp(const arma::mat x_,
{
// update cumulative baseline harzard
// update the increment h_j in the cumulative baseline hazard in each interval

arma::vec xbeta_ = x_ * beta_;
xbeta_.elem(arma::find(xbeta_ > 700)).fill(700.);
arma::vec h_rate = c0_ + sumMatProdVec(ind_r_d_, arma::exp( xbeta_ ));
Expand All @@ -60,30 +110,30 @@ Rcpp::List updateBH_list_cpp(const Rcpp::List x_,
{
// update cumulative baseline harzard
// update the increment h_j in the cumulative baseline hazard in each interval
int S = J_.size();

unsigned int S = J_.size();
Rcpp::List h_(S);

for (unsigned int g = 0; g < S; ++g)
{
arma::mat x_tmp = x_[g];
arma::vec beta_tmp = beta_[g];
int J_tmp = J_[g];
unsigned int J_tmp = J_[g];
arma::mat ind_r_d_tmp = ind_r_d_[g];
arma::vec hPriorSh_tmp = hPriorSh_[g];
arma::vec d_tmp = d_[g];

arma::vec xbeta_ = x_tmp * beta_tmp;
xbeta_.elem(arma::find(xbeta_ > 700)).fill(700.);
arma::vec h_rate = c0_ + sumMatProdVec(ind_r_d_tmp, arma::exp( xbeta_ ));

arma::vec h_tmp = arma::zeros<arma::vec>(J_tmp);
for (unsigned int j = 0; j < J_tmp; ++j)
h_tmp(j) = arma::randg(arma::distr_param(hPriorSh_tmp(j) + d_tmp(j), 1. / h_rate(j)));

h_[g] = h_tmp;
}

return h_;
}

Expand All @@ -99,32 +149,32 @@ Rcpp::List calJpost_helper_cpp(const arma::vec cbtau,
const arma::mat ind_d_)
{
// subfunction to update joint posterior distribution

arma::vec xbeta_ = x_ * beta_;
xbeta_.elem(arma::find(xbeta_ > 700)).fill(700.);
arma::vec exp_xbeta = arma::exp(xbeta_);

double first_sum_ini = arma::accu(-h_ % sumMatProdVec(ind_r_d_, exp_xbeta));

arma::mat h_exp_xbeta_mat = -arma::kron(exp_xbeta, h_.t());
h_exp_xbeta_mat.elem(arma::find(h_exp_xbeta_mat > -1.0e-7)).fill(-1.0e-7);
h_exp_xbeta_mat.elem(arma::find(h_exp_xbeta_mat > -1.0e-7)).fill(be_prop_me_ini-1.0e-7);
h_exp_xbeta_mat = arma::log(1.0 - arma::exp(h_exp_xbeta_mat));
// double second_sum_ini = arma::accu(arma::sum((h_exp_xbeta_mat % ind_d_).t(), 1));
double second_sum_ini = arma::accu(h_exp_xbeta_mat % ind_d_);
double loglike1 = first_sum_ini + second_sum_ini;

double logpriorBeta1 = 0.;
for (unsigned int j = 0; j < beta_.size(); ++j)
{
logpriorBeta1 += arma::log_normpdf( beta_(j), 0.0, cbtau(j) );
}

double logpriorH1 = 0.;
for (unsigned int j = 0; j < h_.size(); ++j)
{
logpriorH1 += R::dgamma( h_(j), hPriorSh_(j), 1. / c0_, true );
}

return Rcpp::List::create(
Rcpp::Named("loglike1") = loglike1,
Rcpp::Named("logpriorBeta1") = logpriorBeta1,
Expand All @@ -148,15 +198,15 @@ Rcpp::List updateRP_genomic_cpp(const unsigned int p,
const double cb)
{
// update coefficients of genomic variables via a MH sampler

arma::uvec updatej = arma::randperm(p);

arma::vec xbeta_ = x_ * be_;
arma::vec sd_be_ = arma::ones<arma::vec>(ga_.n_elem);
sd_be_.elem(arma::find(ga_ == 1.)).fill(cb);
sd_be_ = sd_be_ * tau;
arma::uvec sampleRPg_accept_ = arma::zeros<arma::uvec>(p);

unsigned int j = 0;
for (unsigned int j_id = 0; j_id < p; ++j_id)
{
Expand Down Expand Up @@ -243,7 +293,7 @@ Rcpp::List updateRP_genomic_cpp(const unsigned int p,
sampleRPg_accept_(j) = sampleRPg_accept_(j) + 1;
}
}

return Rcpp::List::create(
Rcpp::Named("be.ini") = be_,
Rcpp::Named("acceptl") = sampleRPg_accept_
Expand Down
97 changes: 59 additions & 38 deletions src/updateRP_genomic_cpp.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,46 +13,67 @@

arma::mat matProdVec(const arma::mat, const arma::vec);
arma::vec sumMatProdVec(const arma::mat, const arma::vec);
Rcpp::List updateRP_genomic_cpp(const unsigned int p,
const arma::mat,
const unsigned int,
arma::mat,
arma::mat,
arma::mat,
arma::vec,
const double,
arma::vec,
arma::vec,
const double,
const double);
Rcpp::List updateRP_genomic_cpp(const unsigned int p,
const arma::mat,
const unsigned int,
arma::mat,
arma::mat,
arma::mat,
arma::vec,
const double,
arma::vec,
arma::vec,
const double,
const double);

double be_prop_me_ini = 0.;
double be_prop_sd_ini = 0.;
double D1 = 0.;
double D2 = 0.;
double D1_prop = 0.;
double D2_prop = 0.;
double loglh_ini = 0.;
double loglh_prop = 0.;
double logprior_prop = 0.;
double logprior_ini = 0.;
double logprop_prop = 0.;
double logprop_ini = 0.;
double logR = 0.;
extern double be_prop_me_ini;
extern double be_prop_sd_ini;
extern double D1;
extern double D2;
extern double D1_prop;
extern double D2_prop;
extern double loglh_ini;
extern double loglh_prop;
extern double logprior_prop;
extern double logprior_ini;
extern double logprop_prop;
extern double logprop_ini;
extern double logR;

double be_prop_me = 1.;
double be_prop_sd = 1.;
extern double be_prop_me;
extern double be_prop_sd;

arma::vec be_prop;
arma::vec exp_xbeta;
arma::mat h_exp_xbeta_mat, h_exp_xbeta_prop_mat;
arma::mat exp_h_exp_xbeta_mat, exp_h_exp_xbeta_prop_mat; //, exp_xbeta_mat;
arma::vec first_sum, second_sum;
arma::vec first_sum_prop, second_sum_prop;
arma::vec x_exp_xbeta, xbeta_prop, exp_xbeta_prop, x_exp_xbeta_prop, x_sq_exp_xbeta, x_sq_exp_xbeta_prop;
arma::vec D1_1st, D1_2nd, D1_1st_prop, D1_2nd_prop;
arma::vec D2_1st, D2_2nd, D2_1st_prop, D2_2nd_prop;
arma::mat D1_2nd_den, D1_2nd_num, D1_2nd_den_prop, D1_2nd_num_prop;
arma::mat D2_2nd_num, D2_2nd_den, D2_2nd_den_prop, D2_2nd_num_prop;
extern arma::vec be_prop;
extern arma::vec exp_xbeta;
extern arma::mat h_exp_xbeta_mat;
extern arma::mat h_exp_xbeta_prop_mat;
extern arma::mat exp_h_exp_xbeta_mat;
extern arma::mat exp_h_exp_xbeta_prop_mat;
extern arma::vec first_sum;
extern arma::vec second_sum;
extern arma::vec first_sum_prop;
extern arma::vec second_sum_prop;
extern arma::vec x_exp_xbeta;
extern arma::vec xbeta_prop;
extern arma::vec exp_xbeta_prop;
extern arma::vec x_exp_xbeta_prop;
extern arma::vec x_sq_exp_xbeta;
extern arma::vec x_sq_exp_xbeta_prop;
extern arma::vec D1_1st;
extern arma::vec D1_2nd;
extern arma::vec D1_1st_prop;
extern arma::vec D1_2nd_prop;
extern arma::vec D2_1st;
extern arma::vec D2_2nd;
extern arma::vec D2_1st_prop;
extern arma::vec D2_2nd_prop;
extern arma::mat D1_2nd_den;
extern arma::mat D1_2nd_num;
extern arma::mat D1_2nd_den_prop;
extern arma::mat D1_2nd_num_prop;
extern arma::mat D2_2nd_num;
extern arma::mat D2_2nd_den;
extern arma::mat D2_2nd_den_prop;
extern arma::mat D2_2nd_num_prop;

#endif

0 comments on commit 83698c9

Please sign in to comment.