diff --git a/R/MizerParams-class.R b/R/MizerParams-class.R index 2bb60be1..74445ccf 100644 --- a/R/MizerParams-class.R +++ b/R/MizerParams-class.R @@ -829,3 +829,11 @@ get_w_min_idx <- function(species_params, w) { names(w_min_idx) <- as.character(species_params$species) w_min_idx } + +# helper function to calculate w_max_idx +get_w_max_idx <- function(species_params, w) { + unlist(lapply(species_params$w_max, + function(w_max, wx) min(max(which(wx <= w_max)) + 1, + length(w)), + wx = w)) +} diff --git a/R/RcppExports.R b/R/RcppExports.R index 4c88e810..79841c09 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,7 +1,7 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -inner_project_loop <- function(no_sp, no_w, n, A, B, S, w_min_idx) { - .Call('_mizer_inner_project_loop', PACKAGE = 'mizer', no_sp, no_w, n, A, B, S, w_min_idx) +inner_project_loop <- function(no_sp, n, A, B, S, w_min_idx, w_max_idx) { + .Call('_mizer_inner_project_loop', PACKAGE = 'mizer', no_sp, n, A, B, S, w_min_idx, w_max_idx) } diff --git a/R/project.R b/R/project.R index eb235e74..fa25ef21 100644 --- a/R/project.R +++ b/R/project.R @@ -317,6 +317,7 @@ project_simple <- no_sp <- nrow(params@species_params) # number of species no_w <- length(params@w) # number of fish size bins idx <- 2:no_w + w_max_idx <- params@other_params$w_max_idx # Hacky shortcut to access the correct element of a 2D array using 1D # notation # This references the egg size bracket for all species, so for example @@ -377,9 +378,10 @@ project_simple <- # for (j in (params@w_min_idx[i]+1):no_w) # n[i,j] <- (S[i,j] - A[i,j]*n[i,j-1]) / B[i,j] # This is implemented via Rcpp - n <- inner_project_loop(no_sp = no_sp, no_w = no_w, n = n, + n <- inner_project_loop(no_sp = no_sp, n = n, A = a, B = b, S = S, - w_min_idx = params@w_min_idx) + w_min_idx = params@w_min_idx, + w_max_idx = w_max_idx) # * Update time ---- t <- t + dt diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 2a79d037..b4c40441 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11,19 +11,19 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif // inner_project_loop -NumericMatrix inner_project_loop(int no_sp, int no_w, NumericMatrix n, NumericMatrix A, NumericMatrix B, NumericMatrix S, NumericVector w_min_idx); -RcppExport SEXP _mizer_inner_project_loop(SEXP no_spSEXP, SEXP no_wSEXP, SEXP nSEXP, SEXP ASEXP, SEXP BSEXP, SEXP SSEXP, SEXP w_min_idxSEXP) { +NumericMatrix inner_project_loop(int no_sp, NumericMatrix n, NumericMatrix A, NumericMatrix B, NumericMatrix S, NumericVector w_min_idx, NumericVector w_max_idx); +RcppExport SEXP _mizer_inner_project_loop(SEXP no_spSEXP, SEXP nSEXP, SEXP ASEXP, SEXP BSEXP, SEXP SSEXP, SEXP w_min_idxSEXP, SEXP w_max_idxSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< int >::type no_sp(no_spSEXP); - Rcpp::traits::input_parameter< int >::type no_w(no_wSEXP); Rcpp::traits::input_parameter< NumericMatrix >::type n(nSEXP); Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP); Rcpp::traits::input_parameter< NumericMatrix >::type B(BSEXP); Rcpp::traits::input_parameter< NumericMatrix >::type S(SSEXP); Rcpp::traits::input_parameter< NumericVector >::type w_min_idx(w_min_idxSEXP); - rcpp_result_gen = Rcpp::wrap(inner_project_loop(no_sp, no_w, n, A, B, S, w_min_idx)); + Rcpp::traits::input_parameter< NumericVector >::type w_max_idx(w_max_idxSEXP); + rcpp_result_gen = Rcpp::wrap(inner_project_loop(no_sp, n, A, B, S, w_min_idx, w_max_idx)); return rcpp_result_gen; END_RCPP } diff --git a/src/inner_project_loop.cpp b/src/inner_project_loop.cpp index 442f2ef4..788da6aa 100644 --- a/src/inner_project_loop.cpp +++ b/src/inner_project_loop.cpp @@ -3,12 +3,16 @@ using namespace Rcpp; // [[Rcpp::export]] -NumericMatrix inner_project_loop(int no_sp, int no_w, - NumericMatrix n, NumericMatrix A, NumericMatrix B, - NumericMatrix S, NumericVector w_min_idx) { +NumericMatrix inner_project_loop(int no_sp, + NumericMatrix n, + NumericMatrix A, + NumericMatrix B, + NumericMatrix S, + NumericVector w_min_idx, + NumericVector w_max_idx) { for (int i = 0; i < no_sp; i++) { - for (int j = w_min_idx[i]; j < no_w; j++) { + for (int j = w_min_idx[i]; j < w_max_idx[i]; j++) { n(i,j) = (S(i,j) - A(i,j)*n(i,j-1)) / B(i,j); } }