Export calc_eqdist

NAMESPACE

@@ -21,6 +21,7 @@ export(calc_HSP)
 export(calc_LAK)
 export(calc_NPR)
 export(calc_POP)
+export(calc_eqdist)
 export(calc_fsel_age)
 export(calc_growth)
 export(calc_index)

R/model_int.R

@@ -522,6 +522,26 @@ conv_mov <- function(x, g, v, na = dim(x)[1], nr = dim(x)[2], aref = ceiling(0.5
   return(aperm(mov_rra, 3:1))
 }
 
+#' Equilibrium distribution from movement matrix
+#'
+#' Applies the movement matrix several times in order to obtain the equilibrium spatial distribution of a movement matrix.
+#' Not used in the model but useful for reporting.
+#' @return Numeric vector of length `nr`
+#' @param x Movement matrix, a square matrix with rows corresponding to origin (sum to 1), and columns corresponding to destination
+#' @param nr Number of regions
+#' @param start The initial distribution. Vector of length `nr`
+#' @param nit Integer, the number of times the movement matrix will be applied
+#' @export
+calc_eqdist <- function(x, nr = dim(x)[2], start = rep(1/nr, nr), nit = 20) {
+  if (inherits(start, "advector") || inherits(x, "advector")) {
+    `[<-` <- RTMB::ADoverload("[<-")
+  }
+
+  N <- array(NA_real_, c(nit, nr))
+  N[1, ] <- start
+  for (i in 2:nit - 1) N[i+1, ] <- colSums(N[i, ] * x)
+  return(N[nit-1, ])
+}
 
 #' Predict the probability of CKMR kinship pairs
 #'

R/report-int-state.R

@@ -670,10 +670,6 @@ plot_mov <- function(fit, s = 1, y, a) {
   mname <- dat@Dlabel@season
 
   mov <- array(fit@report$mov_ymarrs[y, , a, , , s], c(nm, nr, nr))
-  dist_eq <- mov_proj(mov, start = fit@report$recdist_rs[, s])
-
-  #dist_start <- ifelse(dat@Dstock@presence_rs[, s], 1/sum(dat@Dstock@presence_rs[, s]), 0)
-  #dist_eq <- mov_proj(mov, start = dist_start)
 
   if (nm > 1) {
     old_mar <- par()$mar
@@ -685,7 +681,12 @@ plot_mov <- function(fit, s = 1, y, a) {
   }
 
   for(m in 1:nm) {
-    .plot_mov(m = mov[m, , ], p = dist_eq[m, ], rname = rname, xlab = "", ylab = "")
+    dist_eq <- calc_eqdist(mov[m, , ], start = fit@report$recdist_rs[, s])
+
+    #dist_start <- ifelse(dat@Dstock@presence_rs[, s], 1/sum(dat@Dstock@presence_rs[, s]), 0)
+    #dist_eq <- calc_eqdist(mov[m, , ], start = dist_start)
+
+    .plot_mov(m = mov[m, , ], p = dist_eq, rname = rname, xlab = "", ylab = "")
     if (nm > 1) title(mname[m], font.main = 1)
   }
 
@@ -762,19 +763,3 @@ plot_recdist <- function(fit) {
 
   invisible()
 }
-
-
-mov_proj <- function(x, nm = dim(x)[1], nr = dim(x)[2], start = rep(1/nr, nr), nit = 20) {
-  N <- array(NA, c(nit, nm, nr))
-  N[1, 1, ] <- start
-  for (i in 2:nit - 1) {
-    for(m in 1:nm) {
-      if (m == nm) {
-        N[i+1, 1, ] <- N[i, m, ] %*% x[m, , ]
-      } else {
-        N[i, m+1, ] <- N[i, m, ] %*% x[m, , ]
-      }
-    }
-  }
-  return(array(N[nit-1, , ], c(nm, nr)))
-}

_pkgdown.yml

@@ -63,6 +63,7 @@ reference:
   - calc_recruitment
   - conv_mov
   - conv_selpar
+  - calc_eqdist
 
 - title: Likelihood and prior functions
   desc: Calculate the log-likelihood or log-prior density functions from data and estimated parameters.