Using projectToSteady()

The grid calculation in Part A of “toyexample1.Rmd” took a looooong time on my machine. The problem is that for each iteration you run project() for 200 years. A lot of time can be saved by running the projection only until there is not much change any more. This time will be shortened if each time we start at the steady state of the previous set of parameters. The following code does that. It uses projectToSteady() and also uses an environment in which it updates the initial state in the params object at every iteration.

First we need to do some initialisation.

library(mizerExperimental)
if (packageVersion("mizerExperimental") < "0.1.1.9005") {
  remotes::install_github("sizespectrum/mizerExperimental")
}
library(tidyverse)
library(ggplot2)
library(plotly)

#read in fish community size spectrum data for 2000
css<-read.csv("data/commsizespectrum.csv",header=T)[,-1]

We define a function that runs to steady state and returns the distance to the data. It uses an environment to keep its state.

fn <- function (par, dat = css$logden, env = state, tol = 0.1) {
  kappa <- 10^(par[[2]])
  effort <- par[[1]]
  # if kappa has changed a lot, start from scratch
  if (abs(log10(resource_params(env$params)$kappa) - log10(kappa)) > 0.4) {
    state$params <- set_community_model(knife_edge_size = 10, kappa = kappa)
  } else {
    # otherwise keep the previous initial state and just update kappa
    resource_params(env$params)$kappa <- kappa
  }
  env$params@initial_effort <- effort
  
  # run to steady state and update params
  env$params <- 
    projectToSteady(env$params, distance_func = distanceSSLogN,
                    tol = tol, t_max = 200)
  
  # Calculate distance to data
  sse <- sum((log(initialN(env$params)[w_sel]) - css$logden)^2)
  message("log_kappa = ", log10(kappa), " and effort = ", effort,
          " gives sse = ", sse)
  sse
}

We set up the environment in which we’ll keep the state

state <- new.env(parent = emptyenv())
state$params <- set_community_model(knife_edge_size = 10, kappa = 0.05)

Note: No m column in species data frame so using m = 1.

w_sel <- which(state$params@w  >= 16 & state$params@w <= 20000)

We create the grid and evaluate the function at every grid point. We choose to do the equally spaced grid in log10(kappa) instead of kappa.

log_kappa <- seq(from = -4, to = 0, by = 0.5)
effort <- seq(from = 0, to = 2, by = 0.1)
grid <- expand.grid(effort = effort, log_kappa = log_kappa)
grid$SSE <- apply(grid, 1, fn)

A plot of the result:

ggplot(grid, aes(x=effort, y=SSE,col=log_kappa) ) +
  geom_point() 

A better plot of the result using plotly

library(plotly)
z <- reshape2::acast(grid, effort ~ log_kappa)

Using SSE as value column: use value.var to override.

plot_ly(x = ~effort, y = ~log_kappa, z = ~z) %>% 
  add_surface()

The parameters that lead to the least distance to the data:

grid[which.min(grid$SSE), ]

Note that this approach unfortunately does not work well for optim().

Note: No m column in species data frame so using m = 1.
Convergence was achieved in 100.5 years.
log_kappa = -1 and effort = 0.1 gives sse = 923.760027188315
Convergence was achieved in 51 years.
log_kappa = -1 and effort = 0.101 gives sse = 922.648437753245
Convergence was achieved in 64.5 years.
log_kappa = -1 and effort = 0.099 gives sse = 924.864233640294
Convergence was achieved in 49.5 years.
log_kappa = -0.999 and effort = 0.1 gives sse = 924.095325984494
Convergence was achieved in 13.5 years.
log_kappa = -1.001 and effort = 0.1 gives sse = 923.452813663231
Note: No m column in species data frame so using m = 1.
Convergence was achieved in 159 years.
log_kappa = -3 and effort = 3 gives sse = 10776.8821197806
Convergence was achieved in 1.5 years.
log_kappa = -3 and effort = 3 gives sse = 10776.8796731337
Convergence was achieved in 13.5 years.
log_kappa = -3 and effort = 2.999 gives sse = 10771.2521320624
Convergence was achieved in 27 years.
log_kappa = -2.999 and effort = 3 gives sse = 10775.8079040793
Convergence was achieved in 16.5 years.
log_kappa = -3 and effort = 3 gives sse = 10776.8829757653
Note: No m column in species data frame so using m = 1.
Convergence was achieved in 67.5 years.
log_kappa = -1.21441127933837 and effort = 0.410896355040632 gives sse = 489.466966106295
Convergence was achieved in 39 years.
log_kappa = -1.21441127933837 and effort = 0.411896355040632 gives sse = 488.439102529185
Convergence was achieved in 46.5 years.
log_kappa = -1.21441127933837 and effort = 0.409896355040632 gives sse = 490.498480950453
Convergence was achieved in 37.5 years.
log_kappa = -1.21341127933837 and effort = 0.410896355040632 gives sse = 489.675802076089
Convergence was achieved in 9 years.
log_kappa = -1.21541127933837 and effort = 0.410896355040632 gives sse = 489.3290929637
Simulation run did not converge after 199.5 years. Value returned by the distance function was: 0.000538502371229113
log_kappa = -1.54824317805333 and effort = 3 gives sse = 7596.37611278224
Convergence was achieved in 84 years.
log_kappa = -1.54824317805333 and effort = 3 gives sse = 7596.39317190674
Convergence was achieved in 15 years.
log_kappa = -1.54824317805333 and effort = 2.999 gives sse = 7591.79460826313
Convergence was achieved in 37.5 years.
log_kappa = -1.54724317805333 and effort = 3 gives sse = 7595.76604297761
Convergence was achieved in 51 years.
log_kappa = -1.54924317805333 and effort = 3 gives sse = 7597.01680768313
Convergence was achieved in 55.5 years.
log_kappa = -1.24700970942616 and effort = 0.663720349647343 gives sse = 354.824233047493
Convergence was achieved in 31.5 years.
log_kappa = -1.24700970942616 and effort = 0.664720349647343 gives sse = 354.85563496204
Convergence was achieved in 37.5 years.
log_kappa = -1.24700970942616 and effort = 0.662720349647343 gives sse = 354.797086552004
Convergence was achieved in 31.5 years.
log_kappa = -1.24600970942616 and effort = 0.663720349647343 gives sse = 354.936050592302
Convergence was achieved in 9 years.
log_kappa = -1.24800970942616 and effort = 0.663720349647343 gives sse = 354.666211710749
Convergence was achieved in 48 years.
log_kappa = -1.33359423731248 and effort = 0.657741912882538 gives sse = 345.352226047175
Convergence was achieved in 31.5 years.
log_kappa = -1.33359423731248 and effort = 0.658741912882538 gives sse = 345.390611311402
Convergence was achieved in 37.5 years.
log_kappa = -1.33359423731248 and effort = 0.656741912882538 gives sse = 345.31811640736
Convergence was achieved in 31.5 years.
log_kappa = -1.33259423731248 and effort = 0.657741912882538 gives sse = 345.454742093095
Convergence was achieved in 9 years.
log_kappa = -1.33459423731248 and effort = 0.657741912882538 gives sse = 345.20620852828
Convergence was achieved in 54 years.
log_kappa = -1.67993234885775 and effort = 0.633828165823316 gives sse = 315.574858750225
Convergence was achieved in 33 years.
log_kappa = -1.67993234885775 and effort = 0.634828165823316 gives sse = 315.615990740151
Convergence was achieved in 39 years.
log_kappa = -1.67993234885775 and effort = 0.632828165823316 gives sse = 315.538074861688
Convergence was achieved in 31.5 years.
log_kappa = -1.67893234885775 and effort = 0.633828165823316 gives sse = 315.638483689336
Convergence was achieved in 7.5 years.
log_kappa = -1.68093234885775 and effort = 0.633828165823316 gives sse = 315.493530213342
Note: No m column in species data frame so using m = 1.
Convergence was achieved in 58.5 years.
log_kappa = -2.20866553851669 and effort = 0.584862466321276 gives sse = 154.942150119733
Convergence was achieved in 34.5 years.
log_kappa = -2.20866553851669 and effort = 0.585862466321276 gives sse = 155.59883380954
Convergence was achieved in 40.5 years.
log_kappa = -2.20866553851669 and effort = 0.583862466321276 gives sse = 154.289936221254
Convergence was achieved in 33 years.
log_kappa = -2.20766553851669 and effort = 0.584862466321276 gives sse = 154.951029310567
Convergence was achieved in 9 years.
log_kappa = -2.20966553851669 and effort = 0.584862466321276 gives sse = 154.838086262208
Note: No m column in species data frame so using m = 1.
Convergence was achieved in 61.5 years.
log_kappa = -3 and effort = 0.51157740269288 gives sse = 145.829045029557
Convergence was achieved in 36 years.
log_kappa = -3 and effort = 0.51257740269288 gives sse = 146.916974557407
Convergence was achieved in 42 years.
log_kappa = -3 and effort = 0.51057740269288 gives sse = 144.745892627699
Convergence was achieved in 34.5 years.
log_kappa = -2.999 and effort = 0.51157740269288 gives sse = 145.757942312481
Convergence was achieved in 7.5 years.
log_kappa = -3 and effort = 0.51157740269288 gives sse = 145.765077256925
Convergence was achieved in 79.5 years.
log_kappa = -3 and effort = 0.245401621107392 gives sse = 23.8900404297287
Convergence was achieved in 45 years.
log_kappa = -3 and effort = 0.246401621107392 gives sse = 23.7708840506153
Convergence was achieved in 54 years.
log_kappa = -3 and effort = 0.244401621107392 gives sse = 24.0127423953676
Convergence was achieved in 43.5 years.
log_kappa = -2.999 and effort = 0.245401621107392 gives sse = 23.9111517658352
Convergence was achieved in 7.5 years.
log_kappa = -3 and effort = 0.245401621107392 gives sse = 23.8941597218156
Convergence was achieved in 70.5 years.
log_kappa = -3 and effort = 0.27208145012297 gives sse = 21.9737348795937
Convergence was achieved in 43.5 years.
log_kappa = -3 and effort = 0.27308145012297 gives sse = 21.9546785879314
Convergence was achieved in 52.5 years.
log_kappa = -3 and effort = 0.27108145012297 gives sse = 21.9967195701925
Convergence was achieved in 42 years.
log_kappa = -2.999 and effort = 0.27208145012297 gives sse = 21.9872071936021
Convergence was achieved in 7.5 years.
log_kappa = -3 and effort = 0.27208145012297 gives sse = 21.9728004174819
Convergence was achieved in 61.5 years.
log_kappa = -3 and effort = 0.277694807262647 gives sse = 21.917970362921
Convergence was achieved in 43.5 years.
log_kappa = -3 and effort = 0.278694807262647 gives sse = 21.9211895968837
Convergence was achieved in 52.5 years.
log_kappa = -3 and effort = 0.276694807262647 gives sse = 21.9187479365356
Convergence was achieved in 42 years.
log_kappa = -2.999 and effort = 0.277694807262647 gives sse = 21.9297756511815
Convergence was achieved in 7.5 years.
log_kappa = -3 and effort = 0.277694807262647 gives sse = 21.915898133522

$par
[1]  0.2776948 -3.0000000

$value
[1] 21.91797

$counts
function gradient 
      12       12 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

I think the reason this is so slow to converge and needs such a high precision in the steady state is that the calculation of the gradients is too noisy. In the numerical gradient calculations, optim makes small changes of the parameters to approximate the gradient though finite differences. These small parameter changes lead to small changes in the steady state that are smaller than we are able to detect with our way of running to steady state.