# function to convert age-based to length-based FLQuants # using male and female growth curve and a reference age # for the switch # :: # extra feature is that this randomises the age for a given # Monte Carlo sample: age a is actually U[a-0.5,a+0.5] # Richard Hillary r.hillary@imperial.ac.uk

a2l ← function(a.q,lbins,g.m,g.f,ref.age='missing',ref.len='missing',cv.lvb.m='missing',cv.lvb.f='missing') {

# some simple checks

if(quant(a.q) != 'age')
  stop("Error in a2l: input FLQuant is not age-based")
if(missing(ref.age))
  ref.age <- floor(mean(as.numeric(dimnames(a.q$age))))
  
if(missing(ref.len))
  ref.len <- 40
  
if(missing(cv.lvb.m))
  cv.lvb.m <- 0
  
if(missing(cv.lvb.f))
  cv.lvb.f <- 0
  
# OK create the length-based FLQuant

dm <- dim(a.q)
dms <- dimnames(a.q)[2:6]
dms$length <- as.character(lbins)
l.q <- FLQuant(quant='length',dimnames=dms)

# next is to use the growth curve to find which length bin
# the relevant ages sit in

lvb <- function(a,mod) {
  l <- mod[['Linf']] * (1 - exp(- mod[['k']] * (a - mod[['t0']])))
  return(l)
}

# so we have ages (upper and lower) and associated lengths (upper and lower)
# which assumes that the growth curve is monotonically increasing

ages <- as.numeric(dimnames(a.q)$age)
ages.u <- ages + 0.5
ages.l <- ages - 0.5
# OK so now we assign each age to a given length bin

a.flg.m.u <- ages
a.flg.m.l <- ages
a.flg.f.u <- ages
a.flg.f.l <- ages

# now we need to find out the length bins covered by the upper and 
# lower lengths given the upper and lower ages

dms.l <- dimnames(l.q)
dms.a <- dimnames(a.q)
l.q <- array(dim=dim(l.q))
l.q[] <- 0
aa.q <- array(dim=dim(a.q))
aa.q[] <- a.q@.Data[]

for(y in 1:dm[2]) { 
  
  l.m.u <- lvb(ages.u,g.m)*rlnorm(length(ages.u),0,sqrt(log(1+cv.lvb.m^2)))
  l.m.l <- lvb(ages.l,g.m)*rlnorm(length(ages.l),0,sqrt(log(1+cv.lvb.m^2)))

  # for female ages switch to female growth increments at reference age

  l.f.u <- l.m.u
  l.f.l <- l.m.l
  for(a in 1:length(ages.u)) {
    if(ages.u[a] < ref.age)
      l.f.u[a] <- lvb(ages.u[a],g.m)*rlnorm(1,0,log(sqrt(1+cv.lvb.m^2)))
    if(ages.u[a] >= ref.age && l.f.u[a] < ref.len)
      l.f.u[a] <- (l.f.u[a-1] + g.f[['Linf']] * exp(-g.f[['k']]*(ages.u[a]-g.f[['t0']]))*(1 - exp(-g.f[['k']]))) * rlnorm(1,0,log(sqrt(1+cv.lvb.m^2)))      
    if(ages.u[a] >= ref.age && l.f.u[a] >= ref.len)
      l.f.u[a] <- (l.f.u[a-1] + g.f[['Linf']] * exp(-g.f[['k']]*(ages.u[a]-g.f[['t0']]))*(1 - exp(-g.f[['k']]))) * rlnorm(1,0,log(sqrt(1+cv.lvb.f^2)))      
    
  }

  for(a in 1:length(ages.l)) {
    if(ages.l[a] < ref.age)
      l.f.l[a] <- lvb(ages.l[a],g.m)*rlnorm(1,0,log(sqrt(1+cv.lvb.m^2)))
    if(ages.l[a] >= ref.age && l.f.l[a] < ref.len)
      l.f.l[a] <- (l.f.l[a-1] + g.f[['Linf']] * exp(-g.f[['k']]*(ages.l[a]-g.f[['t0']]))*(1 - exp(-g.f[['k']]))) * rlnorm(1,0,log(sqrt(1+cv.lvb.m^2)))
    if(ages.l[a] >= ref.age && l.f.l[a] >= ref.len)
      l.f.l[a] <- (l.f.l[a-1] + g.f[['Linf']] * exp(-g.f[['k']]*(ages.l[a]-g.f[['t0']]))*(1 - exp(-g.f[['k']]))) * rlnorm(1,0,log(sqrt(1+cv.lvb.f^2)))
  }
  
  for(l in 1:length(l.m.u)) {
    if(l.m.u[l] <= l.m.l[l])
      l.m.u[l] <- l.m.l[l]
    if(l.f.u[l] <= l.f.l[l])
      l.f.u[l] <- l.f.l[l]
  }
    
  for(i in 1:length(ages)) {

    for(j in 2:length(lbins)) {
    
      # upper bins
  
      if(l.m.u[i] >= lbins[j-1] && l.m.u[i] < lbins[j])
       a.flg.m.u[i] <- j-1
      if(l.f.u[i] >= lbins[j-1] && l.f.u[i] < lbins[j])
        a.flg.f.u[i] <- j-1
      
      # lower bins
    
      if(l.m.l[i] >= lbins[j-1] && l.m.l[i] < lbins[j])
        a.flg.m.l[i] <- j-1
      if(l.f.l[i] >= lbins[j-1] && l.f.l[i] < lbins[j])
        a.flg.f.l[i] <- j-1
      # plus and minus group stuff
    
      # upper
    
      if(l.m.u[i] < lbins[1])
        a.flg.m.u[i] <- 1
      if(l.f.u[i] < lbins[1])
        a.flg.f.u[i] <- 1
      if(l.m.u[i] >= lbins[length(lbins)])
        a.flg.m.u[i] <- length(lbins)
      if(l.f.u[i] >= lbins[length(lbins)])
        a.flg.f.u[i] <- length(lbins)
      
      # lower
    
      if(l.m.l[i] < lbins[1])
        a.flg.m.l[i] <- 1
      if(l.f.l[i] < lbins[1])
        a.flg.f.l[i] <- 1
      if(l.m.l[i] >= lbins[length(lbins)])
        a.flg.m.l[i] <- length(lbins)
      if(l.f.l[i] >= lbins[length(lbins)])
        a.flg.f.l[i] <- length(lbins)
    }
  }
           
  for(a in 1:length(ages)) {

    ll.m <- a.flg.m.l[a]
    lu.m <- a.flg.m.u[a]
  
    ll.f <- a.flg.f.l[a]
    lu.f <- a.flg.f.u[a]
  
    # randomise the proportion of a given age 'smear' that goes into the selection
    # of length bins with which it coincides with
  
    p.m <- vector(length=(lu.m-ll.m+1))
    p.m[] <- 1/(lu.m-ll.m+1)
    p.m <- p.m * rlnorm(length(p.m),0,sqrt(log(1+0.1^2)))
    p.m <- p.m/sum(p.m)
    p.f <- vector(length=(lu.f-ll.f+1))
    p.f[] <- 1/(lu.f-ll.f+1)
    p.f <- p.f * rlnorm(length(p.f),0,sqrt(log(1+0.1^2)))
    p.f <- p.f/sum(p.f)
  
    for(l in ll.m:lu.m)
      l.q[l,y,1,,,] <- l.q[l,y,1,,,] + aa.q[a,y,1,,,] * p.m[l-ll.m+1]
    for(l in ll.f:lu.f) 
      l.q[l,y,2,,,] <- l.q[l,y,2,,,] + aa.q[a,y,2,,,] * p.f[l-ll.f+1]
  }
}

l.q <- FLQuant(quant='length',l.q,dimnames=dms.l)

return(l.q)

}

 
efimas1/wp4/cs4/main/nepha2l.txt · Last modified: 2008/08/20 17:01 by hd
 
Except where otherwise noted, content on this wiki is licensed under the following license:CC Attribution-Noncommercial-Share Alike 3.0 Unported
Recent changes RSS feed Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki