Open kristinmarshall-NOAA opened 3 years ago
Hi Kristin, I'm just back and will have a look at this issue once I get the hake data together for the assessment. We need to make sure the R vector is not overwritten elsewhere. Maybe a plot showing these would help confirm.
Here is a plot of the relative SSB trajectory - this is 3 runs, which I think confirms that my approach above set recruitment deviations to zero because all three trajectories are identical.
However, these do not reach the unfished level as I would expect. I think that is because Catches are not set to zero for some reason that I haven't figured out yet. The catch time series is plotted below.
So, I need to look into why a constant catch is being removed during the projections, rather than catch=0.
I think I fixed this: 8ede49aae4345f4365839fbe913bed4cb6dbf172
It's late so I didn't have time to run it through a bunch of tests but I did 30 runs and the resulting biomass approaches unfished but doesn't quite get there. Need to figure out why that is. Use set_rec_devs_zero = TRUE
in argument list below. The doc text for it is in load_data_om()
.
run_oms(ss_model = ss_model,
n_runs = 30,
yr_future = 100,
fns = fns,
n_surveys = 2,
b_futures = 0,
c_increases = 0,
m_increases = 0,
sel_changes = 0,
catch_in = 0,
set_rec_devs_zero = TRUE,
plot_names = plotnames,
random_seed = 12345,
results_root_dir = results_root_dir,
results_dir = results_dir)
Thanks @cgrandin. Your fix ultimately had the same effect as my approach, but is obviously a better implementation.
My previous problem was actually that I hadn't set catch_in=0. I missed it because catch_in=0 was the default the last time I ran the bias adjustment scenarios, which were my template for this run, but that changed back in e33001159281389a0f56af7168425789149f2d88. I'm all squared away now.
I think your plot above doesn't reach SSB0 because it's mid-year spawning biomass. I plotted beginning of the year spawning biomass and it's much closer. At mid-year, half the natural mortality would have been applied for the year.
Ahh of course, mid-year biomass. I'm relieved that it's all working properly without needing a major investigation!
Circling back on these runs, I pushed run files and plots to the runhakemse repo. Extracted from that and pasted here is a plot of seasonal ssb by country. The black lines are the spatial ssb0 values specified in the initial conditions as 0.75 and 0.25 of the coastwide ssb0.
This is output from ssb_all in the season_loop: https://github.com/pacific-hake/pacifichakemse/blob/4ad9f370690addd1ddd6c9be29d1ed20a11f880a/R/run_season_loop_om.R#L238
In our team call just now, I was confused by the seasonal cycles because ssb_all is saved after the movement and mortality have occurred in the code. But after looking more closely at the indexing, I see now that, ssb is calculated from the numbers at age in the same season (e.g., season 1 ssb is calculated from season 1 numbers at age). Looking up in the movement code in earlier lines, this matches the indexing for numbers at age at the beginning of the season, BEFORE the movement and mortality have been applied.
So, I think this all aligns as it should, and the seasonal cycles of high biomass in US and low biomass in Canada in season 1 are at the beginning of season 1. Now, I just need to fix the indexing over in the seasonal plot to show this is ssb at the BEGINNING of each season, not at the end as I had originally interpreted.
One of the SRG's suggestions was to explore in more detail the equilibrium assumptions of the spatial movement and recruitment dynamics. I think it may be helpful to build this up from the simplest case, OM-only with no fishing and no recruitment deviations. I think this can be accomplished by starting with the bias adjustment runs that we already have, and setting the recruitment deviations to zero during the projection years. Then, with recruitment deviations set to zero we can check whether long-term equilibrium SSB in the two areas matches matches the spatial SSB0 implied by the spatial R0.
After looking at the code, it looks to me like a one-off way way to accomplish this is by doing a run where we partially overwriting om$r_y after it's read into the year loop: https://github.com/pacific-hake/pacifichakemse/blob/49cdb088cffc8ec062d0e3399e128c74df6e937b/R/run_year_loop_om.R#L29
I ran this locally, but the results don't make sense to me yet so I need to dig in more when I get back from leave.