gadget-framework / gadget3

TMB-based gadget implemtation
GNU General Public License v2.0
8 stars 6 forks source link

g3a_predate_catchablity_hcr #66

Open lentinj opened 2 years ago

lentinj commented 2 years ago

Generate catchability based on a previously generated aggregate. e.g. use last 3 years' of aggregate to determine next 3 years' quota. May be some time lag in this, e.g. data aggregates for q1, but quota changed in q4.

See also https://github.com/gadget-framework/gadget2/blob/4cac91860834df254b0e228da05249ff0227d795/src/proglikelihood.cc

lentinj commented 2 years ago

Annotated copy of a prognosis likelihood component option:-

components = list(list('[component]',name = 'prognosis',
  weight = 1,
  type = "proglikelihood",
  # Fleets / stocks harvested by those fleets
  fleetnames = pre.fleet.names,
  stocknames = c(imm.stock,mat.stock),

  # Proportion of stock biomass making up fishable biomass (bio1)
  # 3 values per-stock, used to weight biomass: https://github.com/gadget-framework/gadget2/blob/4cac91860834df254b0e228da05249ff0227d795/src/stockmemberfunctions.cc#L319
  biocoeffs = c(-100,ref.cm,1,-100,ref.cm,1),

  # Proportions of stock biomass making up trigger biomass (bio2)
  # 3 values per-stock, used to weight biomass: https://github.com/gadget-framework/gadget2/blob/4cac91860834df254b0e228da05249ff0227d795/src/stockmemberfunctions.cc#L319
  triggercoeffs=c(0,0,0,0,0,1),

  # coefficents used in CalcTac
  # If trigger biomass less than first trigger value, use first harvestrate
  # If trigger biomass less than last trigger value, use last harvestrate
  # Otherwise assign harvest rate proportionally based on trigger values the trigger biomass sits within. 
  triggervalues=c(0, sprintf('#%s.btrigger',mat.stock)),
  # 2 or 3 values
  harvestrates = sprintf(c('#%s.hr.low','#%s.hr.high'),ref.fleet),

  # What steps in the year do we apply the calculated TAC to fleet's ``timemultiplier``
  quotasteps = c(3, 4, 5, 6),

  # How to distribute TAC between fished stocks
  quotaproportions = c(0.25, 0.25, 0.25, 0.25),
  # How to distribute TAC between fleets
  fleetproportions= sprintf('#fleet.prop.%s',pre.fleet.names),

  # Proportion of old TAC in new TAC
  weightoflastyearstac=0,

  # Unused
  maxchange=4,
  functionnumber=1,

  # Don't calculate the TAC before this year
  firsttacyear=min(schedule$year),
  # Recalculate TAC at on this step
  assessmentstep=2,

  # Error to apply to bio1/bio2 in TAC calculation
  asserr=gadgetfile('asserr',
    file_type = 'timevariable',
    components=list(list('asserr',
      data = data_frame(year = main[[1]]$timefile[[1]]$firstyear,
          step = main[[1]]$timefile[[1]]$firststep,
          value = "0") %>%
        bind_rows(expand.grid(year = unique(schedule$year), step = 2) %>% dplyr::mutate(value = paste0('#',mat.stock,'.asserr.',year))) %>% as.data.frame()))))))

  # Error to apply to bio1/bio2 in TAC distribution
  implerr=gadgetfile(...)

Steps are:-

lentinj commented 2 years ago

@bthe Does the above make sense? I'm a bit confused about what quotasteps is doing, but beyond that I think I've got a handle of it.

Translating into a g3a_assessment_tac and a corresponding catchability function will be reasonably easy though. Obviously we'd have to use logspace_add rather than having hard breakpoints.

lentinj commented 2 years ago

Breaking this down a bit further...

Having separate g3a_tac_assess and g3a_tac_apply actions means a small lag between assessment and application is pretty trivial, assess has run_f = ~cur_step == 3, apply has run_f = ~cur_step == 2 and we don't apply to the following spring. A whole year would be harder though, if that's needed we'd have to have some option to switch ordering of the assess/apply steps---an apply/assess order on the same step would mean you apply the results of the previous year.

History could be done using existing reporting mechanisms, which I think would be easier than updating the current and all future TACs

lentinj commented 2 years ago

I've made a stab at implementing this and showing roughly how to use in demo-ling in the above. Before going much further it'd be good to go through it and make sure it's doing what you expect.