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..0b801949 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,31 +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) ){
-    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 )
   }
 
   ##############
@@ -202,15 +201,15 @@ function( object,
   ##############
 
   new_control = object$internal$control
-  new_control$run_model = TRUE
-  new_control$nlminb_loops = 0
+  new_control$run_model = FALSE
+  new_control$opt_loops = 0
   new_control$newton_loops = 0
   new_control$getsd = FALSE
   new_control$calculate_deviance_explained = FALSE
   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,
@@ -235,12 +234,15 @@ 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) ] = oldpar[ (names(oldpar) %in% tmb_fixed) ]
-  newrep = newobj$obj$report( newpar )
+  # 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,
+                          silent = TRUE,
+                          DLL = "tinyVAST" )
 
   ##############
   # Step 4: Merge ParList and ParList1
@@ -270,27 +272,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
@@ -363,105 +362,74 @@ function( object,
     return(newdelta_tc)
   }
 
-  #
-  new_parlist = newobj$tmb_inputs$tmb_par
-  #new_parlist = newobj$tmb_par
-  Q_ss = object$rep$Q_ss
-
-  # 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$tmb_inputs$tmb_par,
+                       type = "Fun",
+                       silent = TRUE,
+                       DLL = "tinyVAST" )
 
   ##############
-  # 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
+  # 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 )
 
-  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
-  )
+    # 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) )
 
-  # 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))
-  }
+    # 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
-  )
+    # 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)
 }
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
 
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" )
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..e8e987b9 100644
--- a/src/tinyVAST.cpp
+++ b/src/tinyVAST.cpp
@@ -40,7 +40,6 @@ Type objective_function<Type>::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<Type>::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
@@ -272,18 +276,18 @@ Type objective_function<Type>::operator() (){
 
   // 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 +303,9 @@ Type objective_function<Type>::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<Type> p_i( n_i );
@@ -312,9 +316,9 @@ Type objective_function<Type>::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<Type> pepsilon1_i = multiply_epsilon( Aepsilon_zz, Aepsilon_z, epsilon_stc, p_i.size() ) / exp(log_tau);
-  vector<Type> pomega1_i = multiply_omega( Aomega_zz, Aomega_z, omega_sc, p_i.size() ) / exp(log_tau);
-  vector<Type> pxi1_i = multiply_xi( A_is, xi_sl, W_il ) / exp(log_tau);
+  vector<Type> pepsilon1_i = multiply_epsilon( Aepsilon_zz, Aepsilon_z, epsilon_stc, p_i.size() );
+  vector<Type> pomega1_i = multiply_omega( Aomega_zz, Aomega_z, omega_sc, p_i.size() );
+  vector<Type> pxi1_i = multiply_xi( A_is, xi_sl, W_il );
   vector<Type> pdelta1_i = multiply_delta( delta_tc, t_i, c_i, n_i );
   p_i += palpha1_i;
   p_i += pgamma1_i;
@@ -332,9 +336,9 @@ Type objective_function<Type>::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<Type> pepsilon2_i = multiply_epsilon( Aepsilon_zz, Aepsilon_z, epsilon2_stc, p_i.size() ) / exp(log_tau);
-  vector<Type> pomega2_i = multiply_omega( Aomega_zz, Aomega_z, omega2_sc, p_i.size() ) / exp(log_tau);
-  vector<Type> pxi2_i = multiply_xi( A_is, xi2_sl, W2_il ) / exp(log_tau);
+  vector<Type> pepsilon2_i = multiply_epsilon( Aepsilon_zz, Aepsilon_z, epsilon2_stc, p_i.size() );
+  vector<Type> pomega2_i = multiply_omega( Aomega_zz, Aomega_z, omega2_sc, p_i.size() );
+  vector<Type> pxi2_i = multiply_xi( A_is, xi2_sl, W2_il );
   vector<Type> pdelta2_i = multiply_delta( delta2_tc, t_i, c_i, n_i );
   p2_i += palpha2_i;
   p2_i += pgamma2_i;
@@ -358,14 +362,18 @@ Type objective_function<Type>::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 +400,17 @@ Type objective_function<Type>::operator() (){
   if( n_g > 0 ){
     vector<Type> palpha1_g = X_gj*alpha_j;
     vector<Type> pgamma1_g = Z_gk*gamma_k;
-    vector<Type> pepsilon1_g = multiply_epsilon( AepsilonG_zz, AepsilonG_z, epsilon_stc, palpha1_g.size() ) / exp(log_tau);
-    vector<Type> pomega1_g = multiply_omega( AomegaG_zz, AomegaG_z, omega_sc, palpha1_g.size() ) / exp(log_tau);
-    vector<Type> pxi1_g = multiply_xi( A_gs, xi_sl, W_gl ) / exp(log_tau);
+    vector<Type> pepsilon1_g = multiply_epsilon( AepsilonG_zz, AepsilonG_z, epsilon_stc, palpha1_g.size() );
+    vector<Type> pomega1_g = multiply_omega( AomegaG_zz, AomegaG_z, omega_sc, palpha1_g.size() );
+    vector<Type> pxi1_g = multiply_xi( A_gs, xi_sl, W_gl );
     vector<Type> pdelta1_g = multiply_delta( delta_tc, t_g, c_g, n_g );
     vector<Type> p1_g = palpha1_g + pgamma1_g+ pepsilon1_g + pomega1_g + pdelta1_g + pxi1_g + offset_g ;
     // Second linear predictor
     vector<Type> palpha2_g = X2_gj*alpha2_j;
     vector<Type> pgamma2_g = Z2_gk*gamma2_k;
-    vector<Type> pepsilon2_g = multiply_epsilon( AepsilonG_zz, AepsilonG_z, epsilon2_stc, palpha2_g.size() ) / exp(log_tau);
-    vector<Type> pomega2_g = multiply_omega( AomegaG_zz, AomegaG_z, omega2_sc, palpha2_g.size() ) / exp(log_tau);
-    vector<Type> pxi2_g = multiply_xi( A_gs, xi2_sl, W2_gl ) / exp(log_tau);
+    vector<Type> pepsilon2_g = multiply_epsilon( AepsilonG_zz, AepsilonG_z, epsilon2_stc, palpha2_g.size() );
+    vector<Type> pomega2_g = multiply_omega( AomegaG_zz, AomegaG_z, omega2_sc, palpha2_g.size() );
+    vector<Type> pxi2_g = multiply_xi( A_gs, xi2_sl, W2_gl );
     vector<Type> pdelta2_g = multiply_delta( delta2_tc, t_g, c_g, n_g );
     vector<Type> 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<Type> omega_distribution(
     const Eigen::SparseMatrix<Type> &Gamma_cc,
     const Eigen::SparseMatrix<Type> &Gammainv_cc,
     const Eigen::SparseMatrix<Type> &Q_ss,
+    Type tau,
     Type range,
     const nngp_data_t<Type> &nngp_data,
     Type &nll ){
@@ -470,7 +471,7 @@ tmbutils::array<Type> omega_distribution(
         Eigen::SparseMatrix<Type> 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<Type> 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<Type> xi_sl,
     vector<Type> log_sigmaxi_l,
     const Eigen::SparseMatrix<Type> &Q_ss,
+    Type tau,
     Type range,
     const nngp_data_t<Type> &nngp_data ){
 
@@ -524,7 +526,7 @@ Type xi_distribution(
       }
     }else{
       for( int l=0; l<n_l; l++ ){
-        nll += SCALE( GMRF(Q_ss), exp(log_sigmaxi_l(l)) )( xi_sl.col(l) );
+        nll += SCALE( GMRF(Q_ss), exp(log_sigmaxi_l(l)) / tau )( xi_sl.col(l) );
       }
     }
   }
@@ -542,6 +544,7 @@ tmbutils::array<Type> epsilon_distribution(
     const Eigen::SparseMatrix<Type> &Gamma_hh,
     const Eigen::SparseMatrix<Type> &Gammainv_hh,
     const Eigen::SparseMatrix<Type> &Q_ss,
+    Type tau,
     Type range,
     const nngp_data_t<Type> &nngp_data,
     Type &nll ){
@@ -580,7 +583,7 @@ tmbutils::array<Type> epsilon_distribution(
       Eigen::SparseMatrix<Type> 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<Type> 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<Type>::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<Type>::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<Type>::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<Type>::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<Type>::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<Type>::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<Type>::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<Type>::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<Type>::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-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-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", {
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)