R/projected.exploitation.rate.r

In jae0/snowcrab: Snow crab assessment suite using aegis* and associated projects and data streams.

  projected.exploitation.rate = function( target.landings=NULL, fm.reg, numbers.initial, sizes, yr, region, scenario="mature.only") {
 
    eps = 1e-5
    objective.threshold = 0.1  # tons
    CV.abundance = 0.1 # assume 10% coefficient of variation in numbers and biomass (minimal assumption)  ... this is really a dummy argument to allow biomass estimation, the error is recalculated at a later step
  
    fm.mean = apply(fm.reg, c(1), function(b) { mean( b[(b > 0) & (b< 1)], na.rm=T, trim=0.1 ) } )

    if (scenario=="same.as.last.year") {
        Constraint.ER = fm.mean
    } else if (scenario=="historic" ) {
        fm.binary = ifelse( fm.mean > eps, 1, 0 )
        n = fm.binary * fm.mean
        Constraint.ER = n / sum (n, na.rm=T)
    } else if (scenario=="mature.only" ) {
        fm.binary = fm.mean * 0
        fm.binary[ which( fm.mean >= eps ) ] = 1
        fm.binary[ grep("imm", names(fm.mean)) ] = 0 
        fm.binary[ grep("CC1to2",  names(fm.mean)) ] = 0 
        n = fm.binary * fm.mean
        Constraint.ER = n / sum (n, na.rm=T)
    } else {
        Constraint.ER = fm.mean 
    }
 
    if (target.landings==0) Constraint.ER = fm.mean * 0
  
    Constraint.ER [ !is.finite(Constraint.ER)] = 0
  
    # NOTE for future reference: the optimization functions must live internal to the calling environment !!!
    optim.model.core = function(inits,  ...) {
      exploitation.rate.optim = inits * Constraint.ER 
      exploitation.rate.optim [ which( exploitation.rate.optim > 1 ) ] = 1  # make sure solution does not exceed 100% exploitation
      exploitation.rate.optim [ which( exploitation.rate.optim < eps ) ] = eps  # make sure solution does not exceed 100% exploitation
      numbers.optim = numbers.initial * exploitation.rate.optim   
      biomass.optim = estimate.biomass( numbers.optim, numbers.optim*CV.abundance, sizes, y=yr, r=region )   
      biomass.optim.fishable = subset.biomass (biomass.optim$x, biomass.optim$x*CV.abundance, names(biomass.optim$x), type=scenario)[1]
      error = abs( target.landings - biomass.optim.fishable )
      V = list( error=error, exploitation.rate.optim=exploitation.rate.optim, numbers.optim=numbers.optim)
      return (V)
    }

    optim.model = function(inits, ... ) {
      res = optim.model.core(inits)
     return (res$error)
    } 

    
   # optim.methods = c( "Nelder-Mead", "CG", "BFGS", "L-BFGS-B", "SANN", "PORT", "NLM-Newton" )
   optim.methods = c( "BFGS", "L-BFGS-B", "PORT" ) # PORT seems most robust

    inits = 0.5
    lower = 0
    upper = 10

    params = NULL
    
    for (om in optim.methods) {
      res = NULL
      if (om == "PORT" ) {
        res = nlminb( start=inits, objective=optim.model, lower=lower, upper=upper, 
          control=list(eval.max=500, iter.max=500 ))
        res$value = res$objective
      } else if (om == "NLM-Newton" ) {
        res = nlm( p=inits, f=optim.model)
        res$convergence = ifelse (res$code %in% c(1,2), 0, 1)
        res$value = res$minimum
      } else {
        res = optim( par=inits, fn=optim.model, gr=NULL, method=om, Constraint.ER=Constraint.ER )
      }
      print(om)
      print(res)
      params = rbind ( params, list( 
        optim.method=om, convergence=res$convergence, error=res$value, par=res$par  # convergence if == 0
      ) ) 
    }
    V = extract.solution (params, optim.model.core)
    return ( V )
  }