Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Gadget projections #99

Open
bthe opened this issue Feb 4, 2020 · 8 comments
Open

Gadget projections #99

bthe opened this issue Feb 4, 2020 · 8 comments
Assignees

Comments

@bthe
Copy link
Collaborator

bthe commented Feb 4, 2020

I've started the redevelopment of the forward projections functionality with the more flexible gadget_project_* functions. The logic behind these functions is that the user defines the projections, unlike gadget.forward, in a stepwise manner thus giving greater control over individual processes:

  • gadget_project_time defines the projection horizon, and behind the scene extends the area and time files
  • gadget_project_stock links together two stocks in a SSB - Rec relationships by defining a (TimeVariable) stochastic spawning process. Currently only hockey-stick implemented but should accommodate all spawn function quite easily.
  • gadget_project_fleet defines a projection fleet by lifting the suitability from a already defined fleet in the model. The type of the model can be any of the fleet types in Gadget.
  • gadget_project_recruitment projects stochastically the recruitment based on a AR(1) regression of the estimated recruitment series
  • gadget_project_advice projects the advice error as a log-normally distributed AR(1) process with a given mean, CV and rho
  • gadget_project_ref_point defines the reference points for the simulation (currently only Blim) but can in theory update other parameter values for the spawn function. Fleet based reference points (Btrigger) cannot be defined yet due to a shortcoming of Gadget

All these functions could be applied repeatedly if the number of stocks in the model is greater than 2. @vbartolino and @pamelajwoods would you mind giving these functions a go and report on any shortcomings.

@bthe bthe self-assigned this Feb 4, 2020
@vbartolino
Copy link
Contributor

when recruitment is scheduled in a timestep different from step 1 the following error is returned:

Error: Can't subset columns that don't exist.
✖ The columns `venimm.rec.pre.2019.1`, `venimm.rec.pre.2020.1`, `venimm.rec.pre.2021.1`, `venimm.rec.pre.2022.1`, `venimm.rec.pre.2023.1`, etc. don't exist.
Run `rlang::last_error()` to see where the error occurred.

I think the issue is that gadget_project_recruitment still looks for projecting recruitment in the first step rather than the correct one

dplyr::mutate(year = sprintf('%s.rec.pre.%s.1',stock,.data$year))

@bthe
Copy link
Collaborator Author

bthe commented Feb 6, 2020

Actually the recruitment step is not changed, meaning that if the immature stock was recruited at step 3 it will continue to recruit at step 3. These functions, however, only assume one recruitment step.

The error you have looks like you haven't yet called gadget on the new model variant, so directly after gadget_project_fleets you need to call gadget with gadget_evaluate to get the update parameter file:

gadget_project_time() %>% 
  gadget_project_stocks(imm.file = 'gssimm',mat.file = 'gssmat') %>% 
  gadget_project_fleets() %>% 
  gadget_evaluate(params.out = paste(attr(.,'variant_dir'),'params.pre',sep='/'),
                  params.in = 'WGTS/params.final')

ps: stock venimm sounds pretty ominous.

@vbartolino
Copy link
Contributor

this is the code I'm using (borrowed from the example) and that returns the error

res <- 
      gadget_project_time() %>% 
      gadget_project_stocks(imm.file='venimm', mat.file='venmat') %>% 
      gadget_project_fleets(pre_fleet='comven2', fleet_type="linearfleet") %>% 
      gadget_evaluate(params.out = paste(attr(.,'variant_dir'),'params.pre',sep='/'),
                      params.in = 'WGTS/params.final') %>% 
      gadget_project_recruitment(stock = 'venimm', 
                                 recruitment = fit$stock.recruitment,
                                 params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/'))

If I check the file params.pre which is created along the way I can see that the recruitment parameters

...
venimm.rec.pre.2019.2	           1	   -9999     9999        0
venimm.rec.pre.2020.2	           1	   -9999     9999        0
venimm.rec.pre.2021.2	           1	   -9999     9999        0
...

recruitment happens in step 2 in this model like correctly reported by these parameters but judging from the error the function is looking for venimm.rec.pre.2019.1`, `venimm.rec.pre.2020.1, ....
If I hack the function gadget_projection_recruitment and change the line
dplyr::mutate(year = sprintf('%s.rec.pre.%s.1',stock,.data$year))
with
dplyr::mutate(year = sprintf('%s.rec.pre.%s.2',stock,.data$year))
at least this error disappear. Although I notice that the format of the file params.pre is completely changed...

bthe pushed a commit that referenced this issue Feb 6, 2020
@bthe
Copy link
Collaborator Author

bthe commented Feb 6, 2020

The most recent version I've fixed this issue in gadget_project_recruitment among a few other brainfarts.

@vbartolino
Copy link
Contributor

The fix does the job, great!
Few more to report. Based on the code below (adapted from your example):

res <- 
      gadget_project_time(num_years=15) %>% 
      gadget_project_stocks(imm.file='venimm', mat.file='venmat') %>% 
      gadget_project_fleets(pre_fleet='comven2', fleet_type="linearfleet") %>% 
      gadget_evaluate(params.out = paste(attr(.,'variant_dir'),'params.pre',sep='/'),
                      params.in = 'WGTS/params.final') %>% 
      gadget_project_recruitment(stock = 'venimm', n_replicates=100,
                                 recruitment = fit$stock.recruitment,
                                 params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/')) %>%
      gadget_project_advice(harvest_rate = 1:10/100, n_replicates=100,
                       pre_fleet = "comven2",
                       params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/')) %>%
      gadget_project_output(imm.file='venimm',mat.file='venmat', pre_fleets="comven2") %>% 
      gadget_evaluate(params.in = paste(attr(.,'variant_dir'),'params.pre',sep='/'))  %>%
      read.printfiles(paste(attr(.,'variant_dir'),'out',sep='/'))

I need to specify the pre_fleets argument in gadget_project_output otherwise the gadget_evaluate that follows doesn't work and return the following message

Warning message:
In system(run.string, intern = TRUE, ignore.stderr = ignore.stderr) :
  running command '/usr/local/R350/lib/R/library/gadget2/bin/gadget  -s    -i prj/params.pre  -main prj/main          2>/dev/null' had status 1

The code above run until the second last line. It prints the files under 'out' and at a quick check they look fine, but when I try to read them with read.printfiles I get the following error:

Error in !suppress : invalid argument type
In addition: Warning message:
In .f(.x[[i]], ...) :
  Old style printfile detected, update Gadget to the most recent version

which is surprising given that I've just updated gadget using the installation via the R package gadget2 - very handy by the way ;)

@bthe
Copy link
Collaborator Author

bthe commented Feb 7, 2020

Yes, this is my bad. Note that the notation x %>% f(y=g(.)) is equivalent to x %>% f(a=., y=g(.)). If you don't want x to be considered as the first argument of the function call you need to wrap the call in curly braces i.e like this x %>% {f(y=g(.))}. In our case we need to wrap the call to read.printfiles in curly braces:

... %>%
{read.printfiles(paste(attr(.,'variant_dir'),'out',sep='/'))} %>% 
... 

it will otherwise think that the output from the preceding steps is the output directory and the output directory is supplied into the suppress argument.
Hope this fixes the issue.

@vbartolino
Copy link
Contributor

yes, it fix it, thanks!
Interesting use of {...}, I didn't think about it.
I need to do some more tests, but so far so good and it really improves forecasts functionalities.

@bthe
Copy link
Collaborator Author

bthe commented Apr 22, 2020

I've added a functionality to gadget_project_fleets that allows the user to define a common multiplier to the fleet data:

... %>% 
gadget_project_fleets(pre_fleet = 'fleet1',
                        common_mult = 'effort.pre', ## defining a common multiplier for (all) fleets
                        pre_propotion = 0.5, ## fleet specific scalar
                        ...) %>%
... %>%

You can then just use the gadget_project_advice like the group of fleets is a single fleet:

... %>%
  gadget_project_advice(pre_fleet = 'effort',
                        harvest_rate = 0.2,
                        params.file = paste(attr(.,'variant_dir'),'params.pre',sep='/'),
                        advice_cv = 0.2) %>%
...                   

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants