diff --git a/DESCRIPTION b/DESCRIPTION index 9e1787e5..8ed33a9a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BAS Version: 1.7.2 -Date: 2024-9-14 +Date: 2024-9-16 Title: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling Authors@R: c(person("Merlise", "Clyde", email="clyde@duke.edu", role=c("aut","cre", "cph"), diff --git a/src/glm_mcmcbas.c b/src/glm_mcmcbas.c index c405b9ba..cdd9cd28 100644 --- a/src/glm_mcmcbas.c +++ b/src/glm_mcmcbas.c @@ -32,6 +32,38 @@ SEXP glm_mcmcbas(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights, SEXP NumUnique = PROTECT(allocVector(INTSXP, 1)); ++nProtected; + + PROTECT_INDEX counts_idx; + PROTECT_WITH_INDEX(counts, &counts_idx); + PROTECT_INDEX R2_idx; + PROTECT_WITH_INDEX(R2, &R2_idx); + PROTECT_INDEX shrinkage_idx; + PROTECT_WITH_INDEX(shrinkage, &shrinkage_idx); + PROTECT_INDEX modelspace_idx; + PROTECT_WITH_INDEX(modelspace, &modelspace_idx); + PROTECT_INDEX modeldim_idx; + PROTECT_WITH_INDEX(modeldim, &modeldim_idx); + // PROTECT_INDEX rank_idx; + // PROTECT_WITH_INDEX(rank, &rank_idx); + PROTECT_INDEX beta_idx; + PROTECT_WITH_INDEX(beta, &beta_idx); + PROTECT_INDEX se_idx; + PROTECT_WITH_INDEX(se, &se_idx); + PROTECT_INDEX deviance_idx; + PROTECT_WITH_INDEX(deviance, &deviance_idx); + PROTECT_INDEX modelprobs_idx; + PROTECT_WITH_INDEX(modelprobs, &modelprobs_idx); + PROTECT_INDEX priorprobs_idx; + PROTECT_WITH_INDEX(priorprobs, &priorprobs_idx); + PROTECT_INDEX logmarg_idx; + PROTECT_WITH_INDEX(logmarg, &logmarg_idx); + PROTECT_INDEX sampleprobs_idx; + PROTECT_WITH_INDEX(sampleprobs, &sampleprobs_idx); + PROTECT_INDEX Q_idx; + PROTECT_WITH_INDEX(R2, &Q_idx); + PROTECT_INDEX Rintercept_idx; + PROTECT_WITH_INDEX(Rintercept, &Rintercept_idx); + double *probs, MH=0.0, prior_m=1.0,shrinkage_m, logmargy, postold, postnew; int i, m, n, pmodel_old, *bestmodel; int mcurrent, n_sure; @@ -260,23 +292,29 @@ SEXP glm_mcmcbas(SEXP Y, SEXP X, SEXP Roffset, SEXP Rweights, else {mcurrent = k;} } -// # nocov start + // truncate vectors; legacy code from MCMC should not get to following but + // keep in case other prior choices create models with zero probabilities that + // need to be dropped and mcurrent < k; + // # nocov start if (mcurrent < k) { // truncate vectors; legacy code from MCMC should not get here - SETLENGTH(modelspace, mcurrent); - SETLENGTH(logmarg, mcurrent); - SETLENGTH(modelprobs, mcurrent); - SETLENGTH(priorprobs, mcurrent); - SETLENGTH(sampleprobs, mcurrent); - SETLENGTH(counts, mcurrent); - SETLENGTH(MCMCprobs, mcurrent); - SETLENGTH(beta, mcurrent); - SETLENGTH(se, mcurrent); - SETLENGTH(deviance, mcurrent); - SETLENGTH(Q, mcurrent); - SETLENGTH(shrinkage, mcurrent); - SETLENGTH(modeldim, mcurrent); - SETLENGTH(R2, mcurrent); - SETLENGTH(Rintercept, mcurrent); + k = mcurrent; + REPROTECT(logmarg= Rf_lengthgets(logmarg, mcurrent), logmarg_idx); + REPROTECT(modelprobs= Rf_lengthgets(modelprobs, mcurrent), modelprobs_idx); + REPROTECT(priorprobs= Rf_lengthgets(priorprobs, mcurrent), priorprobs_idx); + REPROTECT(sampleprobs= Rf_lengthgets(sampleprobs, mcurrent), sampleprobs_idx); + REPROTECT(deviance = Rf_lengthgets(deviance, mcurrent), deviance_idx); + REPROTECT(shrinkage = Rf_lengthgets(shrinkage, mcurrent), shrinkage_idx); + REPROTECT(modeldim= Rf_lengthgets(modeldim, mcurrent), modeldim_idx); + REPROTECT(R2= Rf_lengthgets(R2, mcurrent), R2_idx); + REPROTECT(se= Rf_lengthgets(se, mcurrent), se_idx); + // REPROTECT(rank = Rf_lengthgets(rank, mcurrent), rank_idx); + REPROTECT(modelspace = Rf_lengthgets(modelspace, mcurrent), modelspace_idx); + REPROTECT(beta = Rf_lengthgets(beta, mcurrent), beta_idx); + REPROTECT(se= Rf_lengthgets(se, mcurrent), se_idx); + REPROTECT(Q= Rf_lengthgets(Q, mcurrent), Q_idx); + REPROTECT(Rintercept= Rf_lengthgets(Rintercept, mcurrent), Rintercept_idx); + REPROTECT(counts= Rf_lengthgets(counts, mcurrent), counts_idx); + } // # nocov end diff --git a/src/lm_amcmc.c b/src/lm_amcmc.c index 126f4639..99107804 100644 --- a/src/lm_amcmc.c +++ b/src/lm_amcmc.c @@ -36,6 +36,35 @@ SEXP amcmc(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeldim, SEXP sampleprobs = PROTECT(allocVector(REALSXP, nModels)); ++nProtected; SEXP NumUnique = PROTECT(allocVector(INTSXP, 1)); ++nProtected; + + PROTECT_INDEX counts_idx; + PROTECT_WITH_INDEX(counts, &counts_idx); + PROTECT_INDEX R2_idx; + PROTECT_WITH_INDEX(R2, &R2_idx); + PROTECT_INDEX shrinkage_idx; + PROTECT_WITH_INDEX(shrinkage, &shrinkage_idx); + PROTECT_INDEX modelspace_idx; + PROTECT_WITH_INDEX(modelspace, &modelspace_idx); + PROTECT_INDEX modeldim_idx; + PROTECT_WITH_INDEX(modeldim, &modeldim_idx); + PROTECT_INDEX rank_idx; + PROTECT_WITH_INDEX(rank, &rank_idx); + PROTECT_INDEX beta_idx; + PROTECT_WITH_INDEX(beta, &beta_idx); + PROTECT_INDEX se_idx; + PROTECT_WITH_INDEX(se, &se_idx); + PROTECT_INDEX mse_idx; + PROTECT_WITH_INDEX(mse, &mse_idx); + PROTECT_INDEX modelprobs_idx; + PROTECT_WITH_INDEX(modelprobs, &modelprobs_idx); + PROTECT_INDEX priorprobs_idx; + PROTECT_WITH_INDEX(priorprobs, &priorprobs_idx); + PROTECT_INDEX logmarg_idx; + PROTECT_WITH_INDEX(logmarg, &logmarg_idx); + PROTECT_INDEX sampleprobs_idx; + PROTECT_WITH_INDEX(sampleprobs, &sampleprobs_idx); + + Rprintf("Allocating Space for %d Models AMCMC\n", nModels) ; double *Xwork, *Ywork,*wts, *probs, shrinkage_m, mse_m, MH=0.0, prior_m=1.0, @@ -412,19 +441,22 @@ SEXP amcmc(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeldim, SET_STRING_ELT(ANS_names, 0, mkChar("probne0")); if (nUnique < nModels) { - SETLENGTH(modelspace, nUnique); - SETLENGTH(logmarg, nUnique); - SETLENGTH(modelprobs, nUnique); - SETLENGTH(priorprobs, nUnique); - SETLENGTH(sampleprobs, nUnique); - SETLENGTH(counts, nUnique); - SETLENGTH(beta, nUnique); - SETLENGTH(se, nUnique); - SETLENGTH(mse, nUnique); - SETLENGTH(shrinkage, nUnique); - SETLENGTH(modeldim, nUnique); - SETLENGTH(R2, nUnique); - SETLENGTH(rank, nUnique); + nModels = nUnique; + REPROTECT(logmarg= Rf_lengthgets(logmarg, nUnique), logmarg_idx); + REPROTECT(modelprobs= Rf_lengthgets(modelprobs, nUnique), modelprobs_idx); + REPROTECT(priorprobs= Rf_lengthgets(priorprobs, nUnique), priorprobs_idx); + REPROTECT(sampleprobs= Rf_lengthgets(sampleprobs, nUnique), sampleprobs_idx); + REPROTECT(shrinkage = Rf_lengthgets(shrinkage, nUnique), shrinkage_idx); + REPROTECT(modeldim= Rf_lengthgets(modeldim, nUnique), modeldim_idx); + REPROTECT(R2= Rf_lengthgets(R2, nUnique), R2_idx); + REPROTECT(se= Rf_lengthgets(se, nUnique), se_idx); + REPROTECT(mse= Rf_lengthgets(se, nUnique), mse_idx); + REPROTECT(rank = Rf_lengthgets(rank, nUnique), rank_idx); + REPROTECT(modelspace = Rf_lengthgets(modelspace, nUnique), modelspace_idx); + REPROTECT(beta = Rf_lengthgets(beta, nUnique), beta_idx); + REPROTECT(se= Rf_lengthgets(se, nUnique), se_idx); + REPROTECT(counts= Rf_lengthgets(counts, nUnique), counts_idx); + } SET_VECTOR_ELT(ANS, 1, modelspace); SET_STRING_ELT(ANS_names, 1, mkChar("which")); diff --git a/src/lm_mcmc.c b/src/lm_mcmc.c index e9715fa5..0890276b 100644 --- a/src/lm_mcmc.c +++ b/src/lm_mcmc.c @@ -276,7 +276,6 @@ SEXP mcmc_new(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeldim, if (nUnique < nModels) { nModels = nUnique; - SETLENGTH(counts, nUnique); REPROTECT(counts= Rf_lengthgets(counts, nUnique), counts_idx); REPROTECT(logmarg= Rf_lengthgets(logmarg, nUnique), logmarg_idx); REPROTECT(modelprobs= Rf_lengthgets(modelprobs, nUnique), modelprobs_idx); diff --git a/src/lm_mcmcbas.c b/src/lm_mcmcbas.c index 93cf14b8..738d352c 100644 --- a/src/lm_mcmcbas.c +++ b/src/lm_mcmcbas.c @@ -45,6 +45,34 @@ SEXP mcmcbas(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeldim, SEXP logmarg = PROTECT(allocVector(REALSXP, nModels)); ++nProtected; SEXP sampleprobs = PROTECT(allocVector(REALSXP, nModels)); ++nProtected; SEXP NumUnique = PROTECT(allocVector(INTSXP, 1)); ++nProtected; + + PROTECT_INDEX counts_idx; + PROTECT_WITH_INDEX(counts, &counts_idx); + PROTECT_INDEX R2_idx; + PROTECT_WITH_INDEX(R2, &R2_idx); + PROTECT_INDEX shrinkage_idx; + PROTECT_WITH_INDEX(shrinkage, &shrinkage_idx); + PROTECT_INDEX modelspace_idx; + PROTECT_WITH_INDEX(modelspace, &modelspace_idx); + PROTECT_INDEX modeldim_idx; + PROTECT_WITH_INDEX(modeldim, &modeldim_idx); + PROTECT_INDEX rank_idx; + PROTECT_WITH_INDEX(rank, &rank_idx); + PROTECT_INDEX beta_idx; + PROTECT_WITH_INDEX(beta, &beta_idx); + PROTECT_INDEX se_idx; + PROTECT_WITH_INDEX(se, &se_idx); + PROTECT_INDEX mse_idx; + PROTECT_WITH_INDEX(mse, &mse_idx); + PROTECT_INDEX modelprobs_idx; + PROTECT_WITH_INDEX(modelprobs, &modelprobs_idx); + PROTECT_INDEX priorprobs_idx; + PROTECT_WITH_INDEX(priorprobs, &priorprobs_idx); + PROTECT_INDEX logmarg_idx; + PROTECT_WITH_INDEX(logmarg, &logmarg_idx); + PROTECT_INDEX sampleprobs_idx; + PROTECT_WITH_INDEX(sampleprobs, &sampleprobs_idx); + double *Xwork, *Ywork, *wts, *probs, shrinkage_m, mse_m, R2_m, RSquareFull, Rbestmarg, logmargy, MH=0.0, postold, postnew; @@ -385,18 +413,24 @@ SEXP mcmcbas(SEXP Y, SEXP X, SEXP Rweights, SEXP Rprobinit, SEXP Rmodeldim, /*. for when add heridity if (m < k) { // resize k = m; - SETLENGTH(modelspace, m); - SETLENGTH(logmarg, m); - SETLENGTH(modelprobs, m); - SETLENGTH(priorprobs, m); - SETLENGTH(sampleprobs, m); - SETLENGTH(beta, m); - SETLENGTH(se, m); - SETLENGTH(mse, m); - SETLENGTH(shrinkage, m); - SETLENGTH(modeldim, m); - SETLENGTH(R2, m); - SETLENGTH(rank, m); + nModels = nUnique; + REPROTECT(logmarg= Rf_lengthgets(logmarg, m), logmarg_idx); + REPROTECT(modelprobs= Rf_lengthgets(modelprobs, m), modelprobs_idx); + REPROTECT(priorprobs= Rf_lengthgets(priorprobs, m), priorprobs_idx); + REPROTECT(sampleprobs= Rf_lengthgets(sampleprobs, m), sampleprobs_idx); + REPROTECT(deviance = Rf_lengthgets(deviance, m), deviance_idx); + REPROTECT(shrinkage = Rf_lengthgets(shrinkage, m), shrinkage_idx); + REPROTECT(modeldim= Rf_lengthgets(modeldim, m), modeldim_idx); + REPROTECT(R2= Rf_lengthgets(R2, m), R2_idx); + REPROTECT(se= Rf_lengthgets(se, m), se_idx); + // REPROTECT(rank = Rf_lengthgets(rank, m), rank_idx); + REPROTECT(modelspace = Rf_lengthgets(modelspace, m), modelspace_idx); + REPROTECT(beta = Rf_lengthgets(beta, m), beta_idx); + REPROTECT(se= Rf_lengthgets(se, m), se_idx); + REPROTECT(Q= Rf_lengthgets(Q, m), Q_idx); + REPROTECT(Rintercept= Rf_lengthgets(Rintercept, m), Rintercept_idx); + REPROTECT(counts= Rf_lengthgets(counts, m), counts_idx); + // Rprintf("m %d k %d", m, LENGTH(modelprobs)); } */