search

flr / doc

FLR documentation: tutorials, manuals, examples and others
10 stars 10 forks source link

fwdcontrol using "catch" as quantity #25

Open piaschuchert opened 2 years ago

piaschuchert commented 2 years ago

I am using fwdcontrol and fwd for short term forecasts. It works quite well, except when I use the "catch" quantity. I use the fsq for the intermediate year, and then want to have a catch restrict on the following year (such as a TAC rollover). However, the output value for catches is different than the one that I use as a value.

I am attaching my stock object and the code I am running. It does include examples how the SSB and F quantities work perfectly, however not the "catch" one.

Any ideas? FLR issues.zip

yvermard commented 2 years ago

Hi, looks like it's coming from your catch.wt that do not have wt for older ages (after age 7 and already quite different friom landing.wt at age 6). catch(stf1)[,ac(2023)] returns 206 which is the TAC you're asking for? Seems like you should solve catch.wt prb first? Best

iagomosqueira commented 2 years ago

Are you referring to the runs in, say, lines 163-164?

ctrl <- fwdControl(data.frame(year=max(years)+1:3,val=c(fsq,TAC,0),quantity=c('f','catch','f')))
stf1  <- fwd(stf0, ctrl=ctrl, sr=srr)

As Youen pointed out there are ages with zero wt. If you run

verify(stock)

you can see that both catch.wt and stock.wt have weights not greater than 0.

piaschuchert commented 2 years ago

Hello Iago and Youen,

I will look into this tomorrow morning. I cannot do this from my work laptop.

Thanks already,

Pia

From: Iago Mosqueira @.> Sent: 07 July 2022 12:59 To: flr/doc @.> Cc: Schuchert, Pia @.>; Author @.> Subject: Re: [flr/doc] fwdcontrol using "catch" as quantity (Issue #25)

Caution – This email has been received from outside the NICS network. Please ensure you can verify the sender’s name and email address. Treat all attachments and links with caution. FOR INTERNAL NICS STAFF ONLY - If you have any concerns regarding the email please forward to @.**@.>.

Are you referring to the runs in, say, lines 163-164?

ctrl <- fwdControl(data.frame(year=max(years)+1:3,val=c(fsq,TAC,0),quantity=c('f','catch','f')))

stf1 <- fwd(stf0, ctrl=ctrl, sr=srr)

As Youen pointed out there are ages with zero wt. If you run

verify(stock)

you can see that both catch.wt and stock.wt have weights not greater than 0.

— Reply to this email directly, view it on GitHubhttps://github.com/flr/doc/issues/25#issuecomment-1177495421, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AF2B4GEJGCLU2UTGOG7WBY3VS3BBPANCNFSM524UMGJA. You are receiving this because you authored the thread.Message ID: @.***> If you are not the intended recipient of this email (and any attachment), please inform the sender by return email and destroy all copies. Unauthorised access, use, disclosure, storage or copying is not permitted. Any email including its content may be monitored by AFBI. For more information on our Privacy Notice, please see the following link: AFBI Privacy Notice. Communication by internet email is not secure as messages may be intercepted and read by someone else. AFBI cannot guarantee that this message or any attachment is virus free or has not been intercepted and amended. You should perform your own virus checks

piaschuchert commented 2 years ago

Dear Iago and Youen, I did just update the stock object and re-ran the script. However, the error persists. Indeed, it seems that the fields

,Catch=round(c(landings(stf1)[,nyears+2]+discards(stf1)[,nyears+2])) ,Land=round(c(landings(stf1)[,nyears+2])) ,Dis=round(c(discards(stf1)[,nyears+2]))

Catch + Dis seem to add up to the desired catch control, i.e. I want to set the catch control to the TAC of 206 tons, I get Discards 10 tons, total catch 196 tons, and Landings 185 tons (i.e. total catch- discards).

I try to attach the new stock object.... cod7a_forecast.zip

pdolder commented 2 years ago

I'm not familiar with the internals of fwd using a catch constraint, but the stock object has some internal consistency.

e.g.

computeLandings(stock)[,ac(2020)] An object of class "FLQuant" , , unit = unique, season = all, area = unique

 year

age 2020 all 203.58

units: NA

computeDiscards(stock)[,ac(2020)] An object of class "FLQuant" , , unit = unique, season = all, area = unique

 year

age 2020 all 4.7186

units: NA

computeCatch(stock)[,ac(2020)] An object of class "FLQuant" , , unit = unique, season = all, area = unique

 year

age 2020 all 276.42

units: NA

computeLandings(stock)[,ac(2020)]+computeDiscards(stock)[,ac(2020)] An object of class "FLQuant" , , unit = unique, season = all, area = unique

 year

age 2020 all 208.29

units: NA

On Fri, 8 Jul 2022 at 08:09, Pia Schuchert @.***> wrote:

Dear Iago and Youen, I did just update the stock object and re-ran the script. However, the error persists. Indeed, it seems that the fields

,Catch=round(c(landings(stf1)[,nyears+2]+discards(stf1)[,nyears+2])) ,Land=round(c(landings(stf1)[,nyears+2])) ,Dis=round(c(discards(stf1)[,nyears+2]))

Catch + Dis seem to add up to the desired catch control, i.e. I want to set the catch control to the TAC of 206 tons, I get Discards 10 tons, total catch 196 tons, and Landings 185 tons (i.e. total catch- discards).

I try to attach the new stock object.... cod7a_forecast.zip https://github.com/flr/doc/files/9069680/cod7a_forecast.zip

— Reply to this email directly, view it on GitHub https://github.com/flr/doc/issues/25#issuecomment-1178635114, or unsubscribe https://github.com/notifications/unsubscribe-auth/AC4F4SKVAG3G4BISTAFKVCDVS7HZVANCNFSM524UMGJA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

iagomosqueira commented 2 years ago

The catch target in FLash simply solves for the Fbar that will give that catch. Then the catch.n is split between landings.n and discards.n according to the status quo ratios set by stf().

There is a difference in landings.wt, discards.wt and catch.wt in the older ages and projection years that I think is making it fail in the splitting. catch.wt is 10.2470, while discards.wt and landings.wt are 10.3780. catch.wt is assumed to be a weighted mean of the other two, but that does not hold here. If I set landings.wt and discards.wt in 2023 to be the same as catch.wt, then the result is correct.

> catch(stf1)[,'2023']
An object of class "FLQuant"
, , unit = unique, season = all, area = unique

     year
age   2023
  all 206

units:  NA
> landings(stf1)[,'2023'] + discards(stf1)[,'2023']
An object of class "FLQuant"
, , unit = unique, season = all, area = unique

     year
age   2023
  all 206

units:  NA

Running the projection with FLasher returns the correct total catch but slightly different values for landings and discards. This is because the correction done in a slightly different way, but I will have to look carefully and compare,

If you get the weights set correctly, the results should be as expected.

iagomosqueira commented 2 years ago

Forgot to add that this was pointed by the calls to compute() that @pdolder mentioned. Those functions recompute the overal catch, landings and discards slots from the .,n and *.wt ones.

> catch(stf1)[,'2022']
An object of class "FLQuant"
, , unit = unique, season = all, area = unique

     year
age   2022
  all 172.31

units:  NA
> computeCatch(stf1)[,'2022']
An object of class "FLQuant"
, , unit = unique, season = all, area = unique

     year
age   2022
  all 172.31

units:  NA
> computeCatch(stf1, 'all')$catch[,'2022']
An object of class "FLQuant"
, , unit = unique, season = all, area = unique

     year
age   2022
  all 165.61

units:  NA

> computeCatch(stf1, 'all')$catch.wt[,'2022']
An object of class "FLQuant"
, , unit = unique, season = all, area = unique

    year
age  2022
  0   0.10000
  1   0.28333
  2   1.75133
  3   4.02767
  4   5.95333
  5   7.76067
  6   8.59700
  7  10.37800
  8  10.37800
  9  10.37800
  10 10.37800

units:  NA

The last two commands, with the 'all' option, recalculate the three catch slots from landings and discards. As you can see, the weighted average catch.wt does not match what the stock has.

piaschuchert commented 2 years ago

Hello,

thanks, I finally got it sorted... There seems to be something tricky with the first stock object getting results and fitted estimates from the assessment model and then having to re-add the raw values for number at age and weights... The assessment model does not discriminate between landings and discards, the forecast does...

But now I got it finally sorted, thanks a lot!

Pia