Introducing w_max_idx which can speed up the inner loop, see #278.

The savings are minimal however.
sizespectrum · Dec 1, 2023 · a769cb3 · a769cb3
1 parent 0aa6d22
commit a769cb3
Show file tree

Hide file tree

Showing 5 changed files with 26 additions and 12 deletions.
diff --git 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
@@ -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
@@ -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
@@ -11,19 +11,19 @@ Rcpp::Rostream<false>& 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
@@ -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);
         }
     }