From 3a60c74395523a8763690bd2d17a77377ab94d41 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Tue, 29 Sep 2026 12:37:32 -0700 Subject: [PATCH 1/9] Speed up model fitting; see details - With help from Claude's Opus 5.5, did some extensive benchmarking trying to figure out why tinyVAST was slower than sdmTMB for relatively large SPDE models. Here's what we came up with. - Overall, ~25-30% faster - Somewhat reduced memory too Changes: - Vectorize apply() over A-matrix on R side (easy win) - Put tau in GMRF densities - Make DATA_UPDATE(weights_i) optional; it's slower - Change nll lines from -1 * d() to -d(); trivial but strangely seems to interact with DATA_UPDATE speed-up - Models also seem to often take more nlminb iterations to converge; Could be due to parameterization of marginal SD as leading parameter rather than 'tau', but that's a more major structural change that may not be worth it. --- R/fit.R | 16 +++++-- R/predict.R | 4 +- R/project.R | 3 +- man/tinyVASTcontrol.Rd | 7 ++- src/tinyVAST.cpp | 53 +++++++++++++-------- src/utils.h | 35 +++++++------- tests/testthat/test-index-standardization.R | 4 +- 7 files changed, 74 insertions(+), 48 deletions(-) diff --git a/R/fit.R b/R/fit.R index a5e21e8c..3337f4af 100644 --- a/R/fit.R +++ b/R/fit.R @@ -604,7 +604,7 @@ function( formula, # Aepsilon_zz = cbind(Atriplet$i, Atriplet$j, t_i[Atriplet$i], c_i[Atriplet$i]) - which_Arows = which(apply( Aepsilon_zz, MARGIN=1, FUN=function(x) all(!is.na(x)) & any(x>0) )) + which_Arows = which( rowSums(is.na(Aepsilon_zz))==0 & rowSums(Aepsilon_zz>0, na.rm=TRUE)>0 ) which_Arows = which_Arows[ which(Atriplet$x[which_Arows] > 0) ] if( (nrow(spacetime_term_ram$output$ram)==0) & (nrow(delta_spacetime_term_ram$output$ram)==0) ){ which_Arows = numeric(0) @@ -614,7 +614,7 @@ function( formula, # Aomega_zz = cbind(Atriplet$i, Atriplet$j, c_i[Atriplet$i]) - which_Arows = which(apply( Aomega_zz, MARGIN=1, FUN=function(x) all(!is.na(x)) )) + which_Arows = which( rowSums(is.na(Aomega_zz))==0 ) which_Arows = which_Arows[ which(Atriplet$x[which_Arows] > 0) ] if( (nrow(space_term_ram$output$ram)==0) & (nrow(delta_space_term_ram$output$ram)==0) ){ which_Arows = numeric(0) @@ -799,7 +799,8 @@ function( formula, ifelse( control$gmrf_parameterization=="separable", 0, 1), ifelse( isFALSE(control$get_rsr), 0, 1), ifelse( isFALSE(control$extra_reporting), 0, 1), - 0 + 0, + ifelse( isTRUE(control$update_weights), 1, 0) ) # make dat @@ -1307,6 +1308,9 @@ function( formula, #' in the SPDE method when using \code{add_mesh_covariates}. The #' default \code{barrier_stiffness = 0.01} is the value from Bakka et al. #' 2019. +#' @param update_weights Whether to allow \code{weights_i} to be changed via +#' \code{obj$env$data} without re-taping (e.g., for stochastic gradient +#' descent). This slows down fitting for large datasets. #' #' @references #' Bakka, H., Vanhatalo, J., Illian, J., Simpson, D., Rue, H. (2019). Non-stationary Gaussian models with physical barriers. Spatial Statistics, 29, 268-288. \doi{10.1016/j.spasta.2019.01.002} @@ -1343,7 +1347,8 @@ function( opt_loops = 1, extra_reporting = FALSE, use_anisotropy = FALSE, sar_adjacency = "queen", - barrier_stiffness = 0.01 ){ + barrier_stiffness = 0.01, + update_weights = FALSE ){ gmrf_parameterization = match.arg(gmrf_parameterization) #optimizer = match.arg(optimizer) @@ -1383,7 +1388,8 @@ function( opt_loops = 1, extra_reporting = extra_reporting, use_anisotropy = use_anisotropy, sar_adjacency = sar_adjacency, - barrier_stiffness = barrier_stiffness + barrier_stiffness = barrier_stiffness, + update_weights = update_weights ), class = "tinyVASTcontrol" ) } diff --git a/R/predict.R b/R/predict.R index defa5c5b..eddda794 100644 --- a/R/predict.R +++ b/R/predict.R @@ -629,7 +629,7 @@ function( object, # AepsilonG_zz = cbind(predAtriplet$i, predAtriplet$j, t_g[predAtriplet$i], c_g[predAtriplet$i]) - which_Arows = which(apply( AepsilonG_zz, MARGIN=1, FUN=\(x) all(!is.na(x)) & any(x>0) )) + which_Arows = which( rowSums(is.na(AepsilonG_zz))==0 & rowSums(AepsilonG_zz>0, na.rm=TRUE)>0 ) which_Arows = which_Arows[ which(predAtriplet$x[which_Arows] > 0) ] if( (nrow(object$internal$spacetime_term_ram$output$ram)==0) & (nrow(object$internal$delta_spacetime_term_ram$output$ram)==0) ){ which_Arows = numeric(0) @@ -639,7 +639,7 @@ function( object, # AomegaG_zz = cbind(predAtriplet$i, predAtriplet$j, c_g[predAtriplet$i]) - which_Arows = which(apply( AomegaG_zz, MARGIN=1, FUN=\(x) all(!is.na(x)) )) + which_Arows = which( rowSums(is.na(AomegaG_zz))==0 ) which_Arows = which_Arows[ which(predAtriplet$x[which_Arows] > 0) ] if( (nrow(object$internal$space_term_ram$output$ram)==0) & (nrow(object$internal$delta_space_term_ram$output$ram)==0) ){ which_Arows = numeric(0) diff --git a/R/project.R b/R/project.R index 8edbc1e7..b7314d3d 100644 --- a/R/project.R +++ b/R/project.R @@ -366,7 +366,8 @@ function( object, # new_parlist = newobj$tmb_inputs$tmb_par #new_parlist = newobj$tmb_par - Q_ss = object$rep$Q_ss + # random effects are scaled by tau, so their precision is Q_ss * tau^2 + Q_ss = object$rep$Q_ss * exp(2 * object$rep$log_tau) # Replace epsilon new_parlist$epsilon_stc = augment_epsilon( diff --git a/man/tinyVASTcontrol.Rd b/man/tinyVASTcontrol.Rd index d45a31b5..2e284e91 100644 --- a/man/tinyVASTcontrol.Rd +++ b/man/tinyVASTcontrol.Rd @@ -27,7 +27,8 @@ tinyVASTcontrol( extra_reporting = FALSE, use_anisotropy = FALSE, sar_adjacency = "queen", - barrier_stiffness = 0.01 + barrier_stiffness = 0.01, + update_weights = FALSE ) } \arguments{ @@ -116,6 +117,10 @@ rate and resulting decorrelation distance) for barriers relative to normal areas in the SPDE method when using \code{add_mesh_covariates}. The default \code{barrier_stiffness = 0.01} is the value from Bakka et al. 2019.} + +\item{update_weights}{Whether to allow \code{weights_i} to be changed via +\code{obj$env$data} without re-taping (e.g., for stochastic gradient +descent). This slows down fitting for large datasets.} } \value{ An object (list) of class \code{tinyVASTcontrol}, containing either default or diff --git a/src/tinyVAST.cpp b/src/tinyVAST.cpp index 8b58db30..ffb43334 100644 --- a/src/tinyVAST.cpp +++ b/src/tinyVAST.cpp @@ -40,7 +40,6 @@ Type objective_function::operator() (){ DATA_IVECTOR( c_i ); DATA_VECTOR( offset_i ); DATA_VECTOR( weights_i ); - DATA_UPDATE( weights_i ); // Warning: cannot use branching (weights_i > 0) when using DATA_UPDATE DATA_VECTOR( size_i ); // Ignored unless family = binomial DATA_SPARSE_MATRIX( S_kk ); // Sparse penalization matrix DATA_IVECTOR( Sdims ); // Dimensions of blockwise components of S_kk @@ -56,6 +55,11 @@ Type objective_function::operator() (){ // model_options(2)==0: no RSS; model_options(2)==1: yes RSR // model_options(3)==0: no extra reporting; model_options(3)==1: yes extra reporting // model_options(4)==0: no SE for _g; model_options(4)==1: SE for p_g; model_options(4)==2: SE for mu_g + // model_options(5)==0: fixed weights_i; model_options(5)==1: weights_i can be updated without re-taping + bool update_weights = (model_options.size() > 5) && (model_options(5) == 1); + if( update_weights ){ + DATA_UPDATE( weights_i ); // Warning: cannot use branching (weights_i > 0) when using DATA_UPDATE + } DATA_IMATRIX( Aepsilon_zz ); // NAs get converted to -2147483648 DATA_VECTOR( Aepsilon_z ); DATA_IMATRIX( Aomega_zz ); // NAs get converted to -2147483648 @@ -270,20 +274,23 @@ Type objective_function::operator() (){ Eigen::SparseMatrix Gammainv2_cc = make_ram( ram2_space_term, ram2_space_term_start, theta2_z, omega2_sc.dim(1), int(1) ); Eigen::SparseMatrix Gamma2_cc = make_ram( ram2_space_term, ram2_space_term_start, theta2_z, omega2_sc.dim(1), int(2) ); + // Random effects are on the scale of the linear predictor, with tau applied in their densities; + // scaling their projection by 1/tau instead adds per-observation terms to the random-effect Hessian + // space_term omega_sc = omega_distribution( omega_sc, model_options, Rho_cc, - Gamma_cc, Gammainv_cc, Q_ss, + Gamma_cc, Gammainv_cc, Q_ss, exp(log_tau), exp(log_kappa), nngp_data, nll ); // nngp_data, omega2_sc = omega_distribution( omega2_sc, model_options, Rho2_cc, - Gamma2_cc, Gammainv2_cc, Q_ss, + Gamma2_cc, Gammainv2_cc, Q_ss, exp(log_tau), exp(log_kappa), nngp_data, nll ); // nngp_data, // spacetime_term epsilon_stc = epsilon_distribution( epsilon_stc, model_options, Rho_hh, - Gamma_hh, Gammainv_hh, Q_ss, + Gamma_hh, Gammainv_hh, Q_ss, exp(log_tau), exp(log_kappa), nngp_data, nll ); epsilon2_stc = epsilon_distribution( epsilon2_stc, model_options, Rho2_hh, - Gamma2_hh, Gammainv2_hh, Q_ss, + Gamma2_hh, Gammainv2_hh, Q_ss, exp(log_tau), exp(log_kappa), nngp_data, nll ); // time_term @@ -299,9 +306,9 @@ Type objective_function::operator() (){ // Distribution for SVC components nll += xi_distribution( model_options, xi_sl, log_sigmaxi_l, - Q_ss, exp(log_kappa), nngp_data ); + Q_ss, exp(log_tau), exp(log_kappa), nngp_data ); nll += xi_distribution( model_options, xi2_sl, log_sigmaxi2_l, - Q_ss, exp(log_kappa), nngp_data ); + Q_ss, exp(log_tau), exp(log_kappa), nngp_data ); // Linear predictor .. keep partial effects for RSR adjustment below vector p_i( n_i ); @@ -312,9 +319,9 @@ Type objective_function::operator() (){ p_i = offset_i; if(alpha_j.size() > 0){ palpha1_i = X_ij*alpha_j; } if(gamma_k.size() > 0){ pgamma1_i = Z_ik*gamma_k; } - vector pepsilon1_i = multiply_epsilon( Aepsilon_zz, Aepsilon_z, epsilon_stc, p_i.size() ) / exp(log_tau); - vector pomega1_i = multiply_omega( Aomega_zz, Aomega_z, omega_sc, p_i.size() ) / exp(log_tau); - vector pxi1_i = multiply_xi( A_is, xi_sl, W_il ) / exp(log_tau); + vector pepsilon1_i = multiply_epsilon( Aepsilon_zz, Aepsilon_z, epsilon_stc, p_i.size() ); + vector pomega1_i = multiply_omega( Aomega_zz, Aomega_z, omega_sc, p_i.size() ); + vector pxi1_i = multiply_xi( A_is, xi_sl, W_il ); vector pdelta1_i = multiply_delta( delta_tc, t_i, c_i, n_i ); p_i += palpha1_i; p_i += pgamma1_i; @@ -332,9 +339,9 @@ Type objective_function::operator() (){ p2_i.setZero(); if(alpha2_j.size() > 0){ palpha2_i = X2_ij*alpha2_j; } if(gamma2_k.size() > 0){ pgamma2_i = Z2_ik*gamma2_k; } - vector pepsilon2_i = multiply_epsilon( Aepsilon_zz, Aepsilon_z, epsilon2_stc, p_i.size() ) / exp(log_tau); - vector pomega2_i = multiply_omega( Aomega_zz, Aomega_z, omega2_sc, p_i.size() ) / exp(log_tau); - vector pxi2_i = multiply_xi( A_is, xi2_sl, W2_il ) / exp(log_tau); + vector pepsilon2_i = multiply_epsilon( Aepsilon_zz, Aepsilon_z, epsilon2_stc, p_i.size() ); + vector pomega2_i = multiply_omega( Aomega_zz, Aomega_z, omega2_sc, p_i.size() ); + vector pxi2_i = multiply_xi( A_is, xi2_sl, W2_il ); vector pdelta2_i = multiply_delta( delta2_tc, t_i, c_i, n_i ); p2_i += palpha2_i; p2_i += pgamma2_i; @@ -358,14 +365,18 @@ Type objective_function::operator() (){ if( components_e(e_i(i))==1 ){ mu_i(i) = one_predictor_likelihood( y_i(i), p_i(i), size_i(i), link_ez(e_i(i),0), family_ez(e_i(i),0), log_sigma_segment, nll_tmp, devresid, this ); negloglik_i(i) = nll_tmp; - nll += weights_i(i) * nll_tmp; + // Multiplying by weights_i=1 is not free, because it adds per-observation terms to the random-effect Hessian + if( update_weights || (weights_i(i) != Type(1.0)) ) nll_tmp *= weights_i(i); + nll += nll_tmp; deviance += pow( devresid, 2.0 ); devresid_i(i) = devresid; } if( components_e(e_i(i))==2 ){ mu_i(i) = two_predictor_likelihood( y_i(i), p_i(i), p2_i(i), size_i(i), link_ez.row(e_i(i)), family_ez.row(e_i(i)), log_sigma_segment, poislink_e(e_i(i)), nll_tmp, dev, this ); negloglik_i(i) = nll_tmp; - nll += weights_i(i) * nll_tmp; + // Multiplying by weights_i=1 is not free, because it adds per-observation terms to the random-effect Hessian + if( update_weights || (weights_i(i) != Type(1.0)) ) nll_tmp *= weights_i(i); + nll += nll_tmp; deviance += dev; devresid_i(i) = NAN; } @@ -392,17 +403,17 @@ Type objective_function::operator() (){ if( n_g > 0 ){ vector palpha1_g = X_gj*alpha_j; vector pgamma1_g = Z_gk*gamma_k; - vector pepsilon1_g = multiply_epsilon( AepsilonG_zz, AepsilonG_z, epsilon_stc, palpha1_g.size() ) / exp(log_tau); - vector pomega1_g = multiply_omega( AomegaG_zz, AomegaG_z, omega_sc, palpha1_g.size() ) / exp(log_tau); - vector pxi1_g = multiply_xi( A_gs, xi_sl, W_gl ) / exp(log_tau); + vector pepsilon1_g = multiply_epsilon( AepsilonG_zz, AepsilonG_z, epsilon_stc, palpha1_g.size() ); + vector pomega1_g = multiply_omega( AomegaG_zz, AomegaG_z, omega_sc, palpha1_g.size() ); + vector pxi1_g = multiply_xi( A_gs, xi_sl, W_gl ); vector pdelta1_g = multiply_delta( delta_tc, t_g, c_g, n_g ); vector p1_g = palpha1_g + pgamma1_g+ pepsilon1_g + pomega1_g + pdelta1_g + pxi1_g + offset_g ; // Second linear predictor vector palpha2_g = X2_gj*alpha2_j; vector pgamma2_g = Z2_gk*gamma2_k; - vector pepsilon2_g = multiply_epsilon( AepsilonG_zz, AepsilonG_z, epsilon2_stc, palpha2_g.size() ) / exp(log_tau); - vector pomega2_g = multiply_omega( AomegaG_zz, AomegaG_z, omega2_sc, palpha2_g.size() ) / exp(log_tau); - vector pxi2_g = multiply_xi( A_gs, xi2_sl, W2_gl ) / exp(log_tau); + vector pepsilon2_g = multiply_epsilon( AepsilonG_zz, AepsilonG_z, epsilon2_stc, palpha2_g.size() ); + vector pomega2_g = multiply_omega( AomegaG_zz, AomegaG_z, omega2_sc, palpha2_g.size() ); + vector pxi2_g = multiply_xi( A_gs, xi2_sl, W2_gl ); vector pdelta2_g = multiply_delta( delta2_tc, t_g, c_g, n_g ); vector p2_g = palpha2_g + pgamma2_g + pepsilon2_g + pomega2_g + pdelta2_g + pxi2_g; // Combined diff --git a/src/utils.h b/src/utils.h index 48fb16f0..16cdf709 100644 --- a/src/utils.h +++ b/src/utils.h @@ -446,6 +446,7 @@ tmbutils::array omega_distribution( const Eigen::SparseMatrix &Gamma_cc, const Eigen::SparseMatrix &Gammainv_cc, const Eigen::SparseMatrix &Q_ss, + Type tau, Type range, const nngp_data_t &nngp_data, Type &nll ){ @@ -470,7 +471,7 @@ tmbutils::array omega_distribution( Eigen::SparseMatrix Q_cc = IminusRho_cc.transpose() * Vinv2_cc * IminusRho_cc; // GMRF for SEM: separable variable-space - nll += SEPARABLE( GMRF(Q_cc), GMRF(Q_ss) )( omega_sc ); + nll += SCALE( SEPARABLE( GMRF(Q_cc), GMRF(Q_ss) ), Type(1.0)/tau )( omega_sc ); // Including this line with Makevars below seems to cause a crash: // PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) // PKG_CXXFLAGS=$(SHLIB_OPENMP_CXXFLAGS) @@ -483,7 +484,7 @@ tmbutils::array omega_distribution( nll += NNGP( Type(1.0), range, omega_s, nngp_data ); } }else{ - nll += SEPARABLE( GMRF(I_cc), GMRF(Q_ss) )( omega_sc ); + nll += SCALE( SEPARABLE( GMRF(I_cc), GMRF(Q_ss) ), Type(1.0)/tau )( omega_sc ); } // Sparse inverse-product @@ -509,6 +510,7 @@ Type xi_distribution( tmbutils::array xi_sl, vector log_sigmaxi_l, const Eigen::SparseMatrix &Q_ss, + Type tau, Type range, const nngp_data_t &nngp_data ){ @@ -524,7 +526,7 @@ Type xi_distribution( } }else{ for( int l=0; l epsilon_distribution( const Eigen::SparseMatrix &Gamma_hh, const Eigen::SparseMatrix &Gammainv_hh, const Eigen::SparseMatrix &Q_ss, + Type tau, Type range, const nngp_data_t &nngp_data, Type &nll ){ @@ -580,7 +583,7 @@ tmbutils::array epsilon_distribution( Eigen::SparseMatrix Q_hh = IminusRho_hh.transpose() * Vinv2_hh * IminusRho_hh; // GMRF for DSEM: non-separable time-variable, with separable space - nll += SEPARABLE( GMRF(Q_ss), GMRF(Q_hh) )( epsilon_hs ); + nll += SCALE( SEPARABLE( GMRF(Q_ss), GMRF(Q_hh) ), Type(1.0)/tau )( epsilon_hs ); }else{ // Rank-deficient (projection) method if( model_options(0) == 7 ){ @@ -590,7 +593,7 @@ tmbutils::array epsilon_distribution( nll += NNGP( Type(1.0), range, epsilon_s, nngp_data ); } }else{ - nll += SEPARABLE( GMRF(Q_ss), GMRF(I_hh) )( epsilon_hs ); + nll += SCALE( SEPARABLE( GMRF(Q_ss), GMRF(I_hh) ), Type(1.0)/tau )( epsilon_hs ); } // Sparse inverse-product @@ -933,28 +936,28 @@ Type one_predictor_likelihood( // Distribution switch( family ){ case gaussian_family: - nll = -1 * dnorm( y, mu, exp(log_sigma_segment(0)), true ); + nll = -dnorm( y, mu, exp(log_sigma_segment(0)), true ); devresid = y - mu; if(isDouble::value && of->do_simulate){ y = rnorm( mu, exp(log_sigma_segment(0)) ); } break; case tweedie_family: - nll = -1 * dtweedie( y, mu, exp(log_sigma_segment(0)), 1.0 + invlogit(log_sigma_segment(1)), true ); + nll = -dtweedie( y, mu, exp(log_sigma_segment(0)), 1.0 + invlogit(log_sigma_segment(1)), true ); devresid = devresid_tweedie( y, mu, 1.0 + invlogit(log_sigma_segment(1)) ); if(isDouble::value && of->do_simulate){ y = rtweedie( mu, exp(log_sigma_segment(0)), 1.0 + invlogit(log_sigma_segment(1)) ); } break; case lognormal_family: - nll = -1 * dlnorm( y, logmu - 0.5*exp(2.0*log_sigma_segment(0)), exp(log_sigma_segment(0)), true ); + nll = -dlnorm( y, logmu - 0.5*exp(2.0*log_sigma_segment(0)), exp(log_sigma_segment(0)), true ); devresid = log(y) - ( logmu - 0.5*exp(2.0*log_sigma_segment(0)) ); if(isDouble::value && of->do_simulate){ y = exp(rnorm( logmu - 0.5*exp(2.0*log_sigma_segment(0)), exp(log_sigma_segment(0)) )); } break; case poisson_family: - nll = -1 * dpois( y, mu, true ); + nll = -dpois( y, mu, true ); devresid = sign(y - mu) * pow(2*(y*log((Type(1e-10) + y)/mu) - (y-mu)), 0.5); if(isDouble::value && of->do_simulate){ y = rpois( mu ); @@ -966,7 +969,7 @@ Type one_predictor_likelihood( //}else{ // nll = -1 * logmu; //} - nll = -1 * dbinom_custom( y * size, logmu, log_one_minus_mu, size, true ); + nll = -dbinom_custom( y * size, logmu, log_one_minus_mu, size, true ); if(isDouble::value && of->do_simulate){ y = rbinom( size, mu ); } @@ -975,7 +978,7 @@ Type one_predictor_likelihood( devresid = devresid_binom( y, size, mu ); break; case gamma_family: // shape = 1/CV^2; scale = mean*CV^2 - nll = -1 * dgamma( y, exp(-2.0*log_sigma_segment(0)), mu*exp(2.0*log_sigma_segment(0)), true ); + nll = -dgamma( y, exp(-2.0*log_sigma_segment(0)), mu*exp(2.0*log_sigma_segment(0)), true ); devresid = sign(y - mu) * pow(2 * ( (y-mu)/mu - log(y/mu) ), 0.5); if(isDouble::value && of->do_simulate){ y = rgamma( exp(-2.0*log_sigma_segment(0)), mu*exp(2.0*log_sigma_segment(0)) ); @@ -983,7 +986,7 @@ Type one_predictor_likelihood( break; case nbinom1_family: // dnbinom_robust( x, log(mu_i), log(var - mu) ) // var - mu = exp( log(mu) + log(theta) ) = theta * mu --> var = (theta+1) * mu - nll = -1 * dnbinom_robust( y, logmu, logmu + log_sigma_segment(0), true); + nll = -dnbinom_robust( y, logmu, logmu + log_sigma_segment(0), true); devresid = devresid_nbinom2( y, logmu, logmu - log_sigma_segment(0) ); // theta = mu / phi if(isDouble::value && of->do_simulate){ // rnbinom2( mu, var ) @@ -992,7 +995,7 @@ Type one_predictor_likelihood( break; case nbinom2_family: // dnbinom_robust( x, log(mu_i), log(var - mu) ) // var - mu = exp( 2 * log(mu) - log(theta) ) = mu^2 / theta --> var = mu + mu^2 / theta - nll = -1 * dnbinom_robust( y, logmu, Type(2.0) * logmu - log_sigma_segment(0), true); + nll = -dnbinom_robust( y, logmu, Type(2.0) * logmu - log_sigma_segment(0), true); devresid = devresid_nbinom2( y, logmu, log_sigma_segment(0) ); if(isDouble::value && of->do_simulate){ // rnbinom2( mu, var ) @@ -1001,7 +1004,7 @@ Type one_predictor_likelihood( break; case student_family: // dnbinom_robust( x, log(mu_i), log(var - mu) ) // var - mu = exp( 2 * log(mu) - log(theta) ) = mu^2 / theta --> var = mu + mu^2 / theta - nll = -1 * dstudent( y, mu, exp(log_sigma_segment(0)), 1.0 + exp(log_sigma_segment(1)), true); + nll = -dstudent( y, mu, exp(log_sigma_segment(0)), 1.0 + exp(log_sigma_segment(1)), true); devresid = devresid_student( y, mu, exp(log_sigma_segment(0)), 1.0 + exp(log_sigma_segment(1)) ); if(isDouble::value && of->do_simulate){ y = mu + exp(log_sigma_segment(0)) * rt(1.0 + exp(log_sigma_segment(1))); @@ -1063,11 +1066,11 @@ Type two_predictor_likelihood( y = rbinom( Type(1), mu1 ); } if( y == 0 ){ - nll = -1 * log_one_minus_mu1; + nll = -log_one_minus_mu1; dev = -2 * log_one_minus_mu1; } if( y>0 ){ // Not if-else so y>0 triggered when simulating y>0 - nll = -1 * logmu1; + nll = -logmu1; dev = -2 * logmu1; //deviance1_i(i) = -2 * log_mu1(i); switch( family(1) ){ diff --git a/tests/testthat/test-index-standardization.R b/tests/testthat/test-index-standardization.R index 73a931f7..474e7610 100644 --- a/tests/testthat/test-index-standardization.R +++ b/tests/testthat/test-index-standardization.R @@ -48,7 +48,7 @@ test_that("Basic index standardization works", { spatial_domain = mesh, family = delta_gamma(type="poisson-link"), control = tinyVASTcontrol(gmrf="sep") ) - expect_equal( my1$opt, my2$opt, tolerance=0.001 ) + expect_equal( my1$opt[c("par","objective")], my2$opt[c("par","objective")], tolerance=0.001 ) # Predicted sample-weighted total index = integrate_output( my1, @@ -70,7 +70,7 @@ test_that("Basic index standardization works", { spatial_domain = mesh, family = delta_gamma(type="poisson-link"), control = tinyVASTcontrol(gmrf="sep") ) - expect_equal( my3$opt, my4$opt, tolerance=0.001 ) + expect_equal( my3$opt[c("par","objective")], my4$opt[c("par","objective")], tolerance=0.001 ) }) test_that("Index standardization results are identical in VAST, tinyVAST, and sdmTMB", { From dda387c9ccbdab1a7ef6e4a32f03bb6dd931dd26 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Tue, 29 Sep 2026 12:58:58 -0700 Subject: [PATCH 2/9] Fix 2 unit tests --- tests/testthat/test-basic-fits.R | 14 ++++++++++++-- tests/testthat/test-smooths.R | 15 +++++++++------ 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/tests/testthat/test-basic-fits.R b/tests/testthat/test-basic-fits.R index 7220c267..318c4a42 100644 --- a/tests/testthat/test-basic-fits.R +++ b/tests/testthat/test-basic-fits.R @@ -42,8 +42,18 @@ test_that("cv::cv works", { skip_on_covr() # Run cv::cv - CV = cv::cv(out, seed = 123) - expect_equal( as.numeric(CV[['CV crit']]), 1.607874, tolerance = 0.0001 ) + # cv::cv() re-evaluates the call from inside cv + # where the test-local `mesh` isn't visible + # bquote() solves this without assigning to global env. + out_cv <- eval(bquote(tinyVAST( + data = dat, + formula = n ~ s(w), + spatial_domain = .(mesh), # .() looks up mesh in the test environment + space_term = "", + control = tinyVASTcontrol(newton_loops = 1) + ))) + CV = cv::cv(out_cv, seed = 123) + expect_equal( as.numeric(CV[['CV crit']]), 1.60805, tolerance = 0.0001 ) }) # test_that("data_colnames are robust", { diff --git a/tests/testthat/test-smooths.R b/tests/testthat/test-smooths.R index 6a57785f..469f0d5e 100644 --- a/tests/testthat/test-smooths.R +++ b/tests/testthat/test-smooths.R @@ -262,18 +262,21 @@ test_that("A model with s(x, bs = 'cc') works", { }) test_that("A model with s() by variables works", { - set.seed(1) - # For some reason, it doesn't work on CI but does locally skip_on_ci() skip_on_cran() - dat <- mgcv::gamSim(4) - m_mgcv <- mgcv::gam(y ~ fac + s(x2, by = fac) + s(x0), data = dat) + set.seed(1) + n <- 400 + x2 <- runif(n) + fac <- factor(sample(1:3, n, replace = TRUE)) + mu <- ifelse(fac == 1, sin(2 * pi * x2), ifelse(fac == 2, cos(3 * pi * x2), 8 * (x2 - 0.5)^2)) + dat <- data.frame(y = mu + rnorm(n, sd = 0.5), x2 = x2, fac = fac) + m_mgcv <- mgcv::gam(y ~ fac + s(x2, by = fac), data = dat) p_mgcv <- predict(m_mgcv) - m_s <- sdmTMB::sdmTMB(formula = y ~ fac + s(x2, by = fac) + s(x0), data = dat, spatial = "off") + m_s <- sdmTMB::sdmTMB(formula = y ~ fac + s(x2, by = fac), data = dat, spatial = "off") - m_v <- tinyVAST(formula = y ~ fac + s(x2, by = fac) + s(x0), data = dat) + m_v <- tinyVAST(formula = y ~ fac + s(x2, by = fac), data = dat) expect_s3_class(m_v, "tinyVAST") p_m <- predict(m_mgcv) From 71ef69519566b1f3cdbf48134d46712e9c18f249 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Tue, 29 Sep 2026 13:12:58 -0700 Subject: [PATCH 3/9] Delete extra comment --- src/tinyVAST.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/tinyVAST.cpp b/src/tinyVAST.cpp index ffb43334..e8e987b9 100644 --- a/src/tinyVAST.cpp +++ b/src/tinyVAST.cpp @@ -274,9 +274,6 @@ Type objective_function::operator() (){ Eigen::SparseMatrix Gammainv2_cc = make_ram( ram2_space_term, ram2_space_term_start, theta2_z, omega2_sc.dim(1), int(1) ); Eigen::SparseMatrix Gamma2_cc = make_ram( ram2_space_term, ram2_space_term_start, theta2_z, omega2_sc.dim(1), int(2) ); - // Random effects are on the scale of the linear predictor, with tau applied in their densities; - // scaling their projection by 1/tau instead adds per-observation terms to the random-effect Hessian - // space_term omega_sc = omega_distribution( omega_sc, model_options, Rho_cc, Gamma_cc, Gammainv_cc, Q_ss, exp(log_tau), From 9b89e91aca02383baeecadbee6b43f3b24677c26 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Thu, 1 Oct 2026 21:37:17 -0700 Subject: [PATCH 4/9] Fix #141 --- R/project.R | 2 +- R/utility.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/project.R b/R/project.R index b7314d3d..51077bb4 100644 --- a/R/project.R +++ b/R/project.R @@ -203,7 +203,7 @@ function( object, new_control = object$internal$control new_control$run_model = TRUE - new_control$nlminb_loops = 0 + new_control$opt_loops = 0 new_control$newton_loops = 0 new_control$getsd = FALSE new_control$calculate_deviance_explained = FALSE diff --git a/R/utility.R b/R/utility.R index 8076299e..9eaabf1a 100644 --- a/R/utility.R +++ b/R/utility.R @@ -403,7 +403,7 @@ function( x, # control_initial = control - control_initial$nlminb_loops = 0 + control_initial$opt_loops = 0 control_initial$newton_loops = 0 control_initial$run_model = FALSE From 3336e4d5bf50731a01884623fd8a5ac63c28ee99 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Fri, 2 Oct 2026 21:32:33 -0700 Subject: [PATCH 5/9] Use sampled params for future with project(parm_var = TRUE) Future Rho/Gamma and Q_ss were built from the fitted parameters rather than the sampled draw, so parameter uncertainty didn't propagate into projected process errors --- R/project.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/R/project.R b/R/project.R index 51077bb4..4a9166a1 100644 --- a/R/project.R +++ b/R/project.R @@ -173,6 +173,7 @@ function( object, ############## if( isFALSE(parm_var) & isFALSE(past_var) ){ + parvec = object$obj$env$last.par.best parlist = object$internal$parlist } if( isTRUE(parm_var) & isFALSE(past_var) ){ @@ -239,7 +240,7 @@ function( object, newpar = newobj$obj$env$last.par oldpar = object$obj$env$last.par.best tmb_fixed = setdiff(unique(names(newpar)), newobj$tmb_inputs$tmb_random) - newpar[ (names(newpar) %in% tmb_fixed) ] = oldpar[ (names(oldpar) %in% tmb_fixed) ] + newpar[ (names(newpar) %in% tmb_fixed) ] = parvec[ (names(oldpar) %in% tmb_fixed) ] newrep = newobj$obj$report( newpar ) ############## @@ -367,7 +368,7 @@ function( object, new_parlist = newobj$tmb_inputs$tmb_par #new_parlist = newobj$tmb_par # random effects are scaled by tau, so their precision is Q_ss * tau^2 - Q_ss = object$rep$Q_ss * exp(2 * object$rep$log_tau) + Q_ss = newrep$Q_ss * exp(2 * newrep$log_tau) # Replace epsilon new_parlist$epsilon_stc = augment_epsilon( From 3a1235e85276ffb155d23fd2c6df153833f878dd Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Fri, 2 Oct 2026 21:45:26 -0700 Subject: [PATCH 6/9] Add nsim to project() and build TMB objects once across samples *Much* faster --- R/project.R | 210 +++++++++++++++++++++---------------------------- man/project.Rd | 29 +++---- 2 files changed, 103 insertions(+), 136 deletions(-) diff --git a/R/project.R b/R/project.R index 4a9166a1..a432a761 100644 --- a/R/project.R +++ b/R/project.R @@ -87,9 +87,11 @@ function( Q, #' condition of the forecast #' @param parm_var logical indicating whether to re-sample fixed effects from their #' predictive distribution, thus changing the GMRF for future process errors +#' @param nsim number of samples #' #' @return -#' A vector of values corresponding to rows in \code{newdata} +#' A vector of values corresponding to rows in \code{newdata}, or a matrix +#' with a column for each sample when \code{nsim > 1} #' #' @examples #' # Convert to long-form @@ -139,18 +141,15 @@ function( Q, #' extra_times = length(x) + 1:100 #' n_sims = 10 #' newdata = data.frame( "time" = c(seq_along(x),extra_times), "var" = "y" ) -#' Y = NULL -#' for(i in seq_len(n_sims) ){ -#' tmp = project( -#' mytiny, -#' newdata = newdata, -#' extra_times = extra_times, -#' future_var = TRUE, -#' past_var = TRUE, -#' parm_var = TRUE -#' ) -#' Y = cbind(Y, tmp) -#' } +#' Y = project( +#' mytiny, +#' newdata = newdata, +#' extra_times = extra_times, +#' future_var = TRUE, +#' past_var = TRUE, +#' parm_var = TRUE, +#' nsim = n_sims +#' ) #' matplot( x = row(Y), #' y = Y, #' type = "l", lty = "solid", col = "black" ) @@ -164,32 +163,31 @@ function( object, what = "mu_g", future_var = TRUE, past_var = FALSE, - parm_var = FALSE ){ + parm_var = FALSE, + nsim = 1 ){ ############## # Step 1: Generate uncertainty from parm_var and past_var, - # and load into parlist + # and load into columns of parmat ############## if( isFALSE(parm_var) & isFALSE(past_var) ){ - parvec = object$obj$env$last.par.best - parlist = object$internal$parlist + parmat = object$obj$env$last.par.best %o% rep(1, nsim) } if( isTRUE(parm_var) & isFALSE(past_var) ){ stop("option not available") } if( isFALSE(parm_var) & isTRUE(past_var) ){ - parvec = object$obj$env$last.par.best - MC = object$obj$env$MC( keep=TRUE, n=1, antithetic=FALSE ) - parvec[object$obj$env$lrandom()] = attr(MC, "samples") - parlist = object$obj$env$parList( par = parvec ) + parmat = object$obj$env$last.par.best %o% rep(1, nsim) + MC = object$obj$env$MC( keep=TRUE, n=nsim, antithetic=FALSE ) + parmat[object$obj$env$lrandom(),] = attr(MC, "samples") } if( isTRUE(parm_var) & isTRUE(past_var) ){ if(is.null(object$sdrep$jointPrecision)) stop("Rerun with `getJointPrecision=TRUE`") - parvec = rmvnorm_prec( mu = object$obj$env$last.par.best, - prec = object$sdrep$jointPrecision ) - parlist = object$obj$env$parList( par = parvec ) + parmat = rmvnorm_prec( mu = object$obj$env$last.par.best, + prec = object$sdrep$jointPrecision, + n = nsim ) } ############## @@ -236,13 +234,6 @@ function( object, development = object$internal$development ) - # - newpar = newobj$obj$env$last.par - oldpar = object$obj$env$last.par.best - tmb_fixed = setdiff(unique(names(newpar)), newobj$tmb_inputs$tmb_random) - newpar[ (names(newpar) %in% tmb_fixed) ] = parvec[ (names(oldpar) %in% tmb_fixed) ] - newrep = newobj$obj$report( newpar ) - ############## # Step 4: Merge ParList and ParList1 ############## @@ -364,106 +355,81 @@ function( object, return(newdelta_tc) } - # - new_parlist = newobj$tmb_inputs$tmb_par - #new_parlist = newobj$tmb_par - # random effects are scaled by tau, so their precision is Q_ss * tau^2 - Q_ss = newrep$Q_ss * exp(2 * newrep$log_tau) - - # Replace epsilon - new_parlist$epsilon_stc = augment_epsilon( - #beta_z = parlist$beta_z, - eps_stc = parlist$epsilon_stc, - neweps_stc = new_parlist$epsilon_stc, - #model = object$internal$spacetime_term_ram$output$model, - linpred = 1 - ) - new_parlist$epsilon2_stc = augment_epsilon( - #beta_z = parlist$beta2_z, - eps_stc = parlist$epsilon2_stc, - neweps_stc = new_parlist$epsilon2_stc, - #model = object$internal$delta_spacetime_term_ram$output$model, - linpred = 2 - ) - - # Replace delta - new_parlist$delta_tc = augment_delta( - #nu_z = parlist$nu_z, - delta_tc = parlist$delta_tc, - newdelta_tc = new_parlist$delta_tc, - #model = object$internal$time_term_ram$output$model, - linpred = 1 - ) - new_parlist$delta2_tc = augment_delta( - #nu_z = parlist$nu2_z, - delta_tc = parlist$delta2_tc, - newdelta_tc = new_parlist$delta2_tc, - #model = object$internal$delta_time_term_ram$output$model, - linpred = 2 - ) + ############## + # Step 5: Re-build model + ############## - # Replace other variables that are not changed - same_vars = setdiff( names(new_parlist), c("epsilon_stc","epsilon2_stc","delta_tc","delta2_tc") ) - new_parlist[same_vars] = parlist[same_vars] + # Build prediction object once, and then REPORT for each sample + predobj = MakeADFun( data = add_predictions( object = newobj, newdata = newdata ), + parameters = newobj$internal$parlist, + map = newobj$tmb_inputs$tmb_map, + random = newobj$tmb_inputs$tmb_random, + profile = newobj$internal$control$profile, + DLL = "tinyVAST" ) + predobj$env$beSilent() ############## - # Step 5: Re-build model + # Step 6: simulate samples ############## - #new_control$run_model = TRUE - new_control$tmb_par = new_parlist - new_control$extra_reporting = FALSE - #new_control$nlminb_loops = 0 - #new_control$newton_loops = 0 - #new_control$getsd = FALSE - #new_control$calculate_deviance_explained = FALSE + pred = matrix( NA, nrow = nrow(newdata), ncol = nsim ) + for( i in seq_len(nsim) ){ + parvec = parmat[,i] + parlist = object$obj$env$parList( par = parvec ) - newobj = tinyVAST( - formula = object$formula, - data = object$data, - time_term = object$internal$time_term, - space_term = object$internal$space_term, - spacetime_term = object$internal$spacetime_term, - family = object$internal$family, - space_columns = object$internal$space_columns, - spatial_domain = object$spatial_domain, - time_column = object$internal$time_column, - times = all_times, - variable_column = object$internal$variable_column, - variables = object$internal$variables, - distribution_column = object$internal$distribution_column, - delta_options = list( formula = object$internal$delta_formula, - space_term = object$internal$delta_space_term, - time_term = object$internal$delta_time_term, - spacetime_term = object$internal$delta_spacetime_term, - spatial_varying = object$internal$delta_spatial_varying ), - spatial_varying = object$internal$spatial_varying, - weights = object$internal$weights, - control = new_control, - development = object$internal$development - ) + # + newpar = newobj$obj$env$last.par + oldpar = object$obj$env$last.par.best + tmb_fixed = setdiff(unique(names(newpar)), newobj$tmb_inputs$tmb_random) + newpar[ (names(newpar) %in% tmb_fixed) ] = parvec[ (names(oldpar) %in% tmb_fixed) ] + newrep = newobj$obj$report( newpar ) - # CHECK SIZES ... if throwing error - if( FALSE ){ - # Check length - range( - # NEW - sapply( newobj$tmb_inputs$tmb_par, length ) - - # OLD - sapply( new_parlist, length ) - ) # - isTRUE(all.equal(new_parlist, newobj$tmb_inputs$tmb_par, tolerance=Inf)) - } + new_parlist = newobj$tmb_inputs$tmb_par + #new_parlist = newobj$tmb_par + # random effects are scaled by tau, so their precision is Q_ss * tau^2 + Q_ss = newrep$Q_ss * exp(2 * newrep$log_tau) + + # Replace epsilon + new_parlist$epsilon_stc = augment_epsilon( + #beta_z = parlist$beta_z, + eps_stc = parlist$epsilon_stc, + neweps_stc = new_parlist$epsilon_stc, + #model = object$internal$spacetime_term_ram$output$model, + linpred = 1 + ) + new_parlist$epsilon2_stc = augment_epsilon( + #beta_z = parlist$beta2_z, + eps_stc = parlist$epsilon2_stc, + neweps_stc = new_parlist$epsilon2_stc, + #model = object$internal$delta_spacetime_term_ram$output$model, + linpred = 2 + ) - ############## - # Step 6: simulate samples - ############## + # Replace delta + new_parlist$delta_tc = augment_delta( + #nu_z = parlist$nu_z, + delta_tc = parlist$delta_tc, + newdelta_tc = new_parlist$delta_tc, + #model = object$internal$time_term_ram$output$model, + linpred = 1 + ) + new_parlist$delta2_tc = augment_delta( + #nu_z = parlist$nu2_z, + delta_tc = parlist$delta2_tc, + newdelta_tc = new_parlist$delta2_tc, + #model = object$internal$delta_time_term_ram$output$model, + linpred = 2 + ) - pred = predict( - object = newobj, - newdata = newdata, - what = what - ) + # Replace other variables that are not changed + same_vars = setdiff( names(new_parlist), c("epsilon_stc","epsilon2_stc","delta_tc","delta2_tc") ) + new_parlist[same_vars] = parlist[same_vars] + + # Load into newpar and simulate sample + for( var in newobj$tmb_inputs$tmb_random ) newpar[ names(newpar) == var ] = new_parlist[[var]] + pred[,i] = predobj$report( newpar )[[what]] + } + if( nsim == 1 ) pred = pred[,1] return(pred) } diff --git a/man/project.Rd b/man/project.Rd index f6f6ac22..7ec7f25e 100644 --- a/man/project.Rd +++ b/man/project.Rd @@ -11,7 +11,8 @@ project( what = "mu_g", future_var = TRUE, past_var = FALSE, - parm_var = FALSE + parm_var = FALSE, + nsim = 1 ) } \arguments{ @@ -46,9 +47,12 @@ condition of the forecast} \item{parm_var}{logical indicating whether to re-sample fixed effects from their predictive distribution, thus changing the GMRF for future process errors} + +\item{nsim}{number of samples} } \value{ -A vector of values corresponding to rows in \code{newdata} +A vector of values corresponding to rows in \code{newdata}, or a matrix +with a column for each sample when \code{nsim > 1} } \description{ Projects a fitted model forward in time. @@ -101,18 +105,15 @@ plot( x = seq_along(Y), extra_times = length(x) + 1:100 n_sims = 10 newdata = data.frame( "time" = c(seq_along(x),extra_times), "var" = "y" ) -Y = NULL -for(i in seq_len(n_sims) ){ - tmp = project( - mytiny, - newdata = newdata, - extra_times = extra_times, - future_var = TRUE, - past_var = TRUE, - parm_var = TRUE - ) - Y = cbind(Y, tmp) -} +Y = project( + mytiny, + newdata = newdata, + extra_times = extra_times, + future_var = TRUE, + past_var = TRUE, + parm_var = TRUE, + nsim = n_sims +) matplot( x = row(Y), y = Y, type = "l", lty = "solid", col = "black" ) From 8cdd288f6efd22ff888f1da3a7bf21a12ce50be8 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Fri, 2 Oct 2026 22:49:48 -0700 Subject: [PATCH 7/9] Condition future epsilon on separate time and space precisions in project() --- R/project.R | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/R/project.R b/R/project.R index a432a761..96bd0d6d 100644 --- a/R/project.R +++ b/R/project.R @@ -262,27 +262,24 @@ function( object, IminusRho_hh = Matrix::Diagonal(n=nrow(newrep$Rho_hh)) - newrep$Rho2_hh Q_kk = Matrix::t(IminusRho_hh) %*% Matrix::solve(Matrix::t(newrep$Gamma2_hh) %*% newrep$Gamma2_hh) %*% IminusRho_hh } - Q_hh = Matrix::kronecker( Q_kk, Q_ss ) # - grid = expand.grid( s = seq_len(dim(neweps_stc)[1]), - t = all_times, + grid = expand.grid( t = all_times, c = object$internal$variables ) - grid$num = seq_len(prod(dim(neweps_stc))) + grid$num = seq_len(nrow(grid)) observed_idx = subset( grid, t %in% object$internal$times )$num - - # - tmp = conditional_gmrf( - Q = Q_hh, - observed_idx = observed_idx, - x_obs = as.vector( eps_stc ), - n_sims = 1, - what = ifelse(future_var, "simulate", "predict") - ) + unobserved_idx = setdiff( grid$num, observed_idx ) + + # Precision is kronecker( Q_kk, Q_ss ) and past times include all sites, so the + # conditional mean only needs Q_kk and the conditional precision is kronecker( Q_uu, Q_ss ) + Q_uu = Q_kk[unobserved_idx, unobserved_idx, drop = FALSE] + Q_uo = Q_kk[unobserved_idx, observed_idx, drop = FALSE] + eps_so = matrix( eps_stc, nrow = dim(eps_stc)[1] ) + simeps_h = -1 * as.matrix( eps_so %*% Matrix::t(Matrix::solve(Q_uu, Q_uo)) ) if( isTRUE(future_var) ){ - simeps_h = tmp[,1] - }else{ - simeps_h = tmp$mean + # Columns have precision Q_ss, and are then correlated to have precision Q_uu among columns + z_su = rmvnorm_prec( prec = Q_ss, n = length(unobserved_idx) ) + simeps_h = simeps_h + t( backsolve(chol(as.matrix(Q_uu)), t(z_su)) ) } # Compile From 81159463b7b356dcf6cbcc6404621170c3eeb9d3 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Fri, 2 Oct 2026 22:58:19 -0700 Subject: [PATCH 8/9] Use report-only TMB objects in project() Faster again --- R/project.R | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/R/project.R b/R/project.R index 96bd0d6d..4b1caefc 100644 --- a/R/project.R +++ b/R/project.R @@ -201,7 +201,7 @@ function( object, ############## new_control = object$internal$control - new_control$run_model = TRUE + new_control$run_model = FALSE new_control$opt_loops = 0 new_control$newton_loops = 0 new_control$getsd = FALSE @@ -209,7 +209,7 @@ function( object, new_control$suppress_user_warnings = TRUE new_control$extra_reporting = TRUE - newobj = tinyVAST( + new_inputs = tinyVAST( formula = object$formula, data = object$data, time_term = object$internal$time_term, @@ -234,6 +234,16 @@ function( object, development = object$internal$development ) + # Copy object with padded times, and TMB object without random effects for REPORT + newobj = object + newobj$tmb_inputs = new_inputs + newobj$internal$times = all_times + newobj$obj = MakeADFun( data = new_inputs$tmb_data, + parameters = new_inputs$tmb_par, + map = new_inputs$tmb_map, + silent = TRUE, + DLL = "tinyVAST" ) + ############## # Step 4: Merge ParList and ParList1 ############## @@ -358,12 +368,11 @@ function( object, # Build prediction object once, and then REPORT for each sample predobj = MakeADFun( data = add_predictions( object = newobj, newdata = newdata ), - parameters = newobj$internal$parlist, + parameters = newobj$tmb_inputs$tmb_par, map = newobj$tmb_inputs$tmb_map, - random = newobj$tmb_inputs$tmb_random, - profile = newobj$internal$control$profile, + type = "Fun", + silent = TRUE, DLL = "tinyVAST" ) - predobj$env$beSilent() ############## # Step 6: simulate samples From deb82877d13a63b6f4c47b391d7d1dee62916891 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Mon, 5 Oct 2026 12:52:42 -0700 Subject: [PATCH 9/9] Pass full unmapped parameter lists to report objects in project() --- R/project.R | 30 ++++++++++++------------------ 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/R/project.R b/R/project.R index 4b1caefc..0b801949 100644 --- a/R/project.R +++ b/R/project.R @@ -234,13 +234,13 @@ function( object, development = object$internal$development ) - # Copy object with padded times, and TMB object without random effects for REPORT + # REPORT uses expanded parameter lists, so omit both maps and random effects. + # parList() on the fitted object restores fixed and shared mapped values. newobj = object newobj$tmb_inputs = new_inputs newobj$internal$times = all_times newobj$obj = MakeADFun( data = new_inputs$tmb_data, parameters = new_inputs$tmb_par, - map = new_inputs$tmb_map, silent = TRUE, DLL = "tinyVAST" ) @@ -369,7 +369,6 @@ function( object, # Build prediction object once, and then REPORT for each sample predobj = MakeADFun( data = add_predictions( object = newobj, newdata = newdata ), parameters = newobj$tmb_inputs$tmb_par, - map = newobj$tmb_inputs$tmb_map, type = "Fun", silent = TRUE, DLL = "tinyVAST" ) @@ -378,21 +377,21 @@ function( object, # Step 6: simulate samples ############## + # Use TMB's parameter order when flattening lists for REPORT. + par_template = newobj$obj$env$parList() + same_vars = setdiff( names(par_template), c("epsilon_stc","epsilon2_stc","delta_tc","delta2_tc") ) pred = matrix( NA, nrow = nrow(newdata), ncol = nsim ) for( i in seq_len(nsim) ){ parvec = parmat[,i] parlist = object$obj$env$parList( par = parvec ) - # - newpar = newobj$obj$env$last.par - oldpar = object$obj$env$last.par.best - tmb_fixed = setdiff(unique(names(newpar)), newobj$tmb_inputs$tmb_random) - newpar[ (names(newpar) %in% tmb_fixed) ] = parvec[ (names(oldpar) %in% tmb_fixed) ] - newrep = newobj$obj$report( newpar ) + # Copy parameters with unchanged dimensions before building future precisions. + # Process arrays retain their padded dimensions until augmented below. + new_parlist = par_template + new_parlist[same_vars] = parlist[same_vars] + newrep = newobj$obj$report( unlist(new_parlist, use.names = FALSE) ) # - new_parlist = newobj$tmb_inputs$tmb_par - #new_parlist = newobj$tmb_par # random effects are scaled by tau, so their precision is Q_ss * tau^2 Q_ss = newrep$Q_ss * exp(2 * newrep$log_tau) @@ -428,13 +427,8 @@ function( object, linpred = 2 ) - # Replace other variables that are not changed - same_vars = setdiff( names(new_parlist), c("epsilon_stc","epsilon2_stc","delta_tc","delta2_tc") ) - new_parlist[same_vars] = parlist[same_vars] - - # Load into newpar and simulate sample - for( var in newobj$tmb_inputs$tmb_random ) newpar[ names(newpar) == var ] = new_parlist[[var]] - pred[,i] = predobj$report( newpar )[[what]] + # Include every augmented effect, including effects fitted as fixed parameters. + pred[,i] = predobj$report( unlist(new_parlist, use.names = FALSE) )[[what]] } if( nsim == 1 ) pred = pred[,1] return(pred)