make_gadget_model: Assemble an entire Gadget model with a single function

In inspktrgadget/gadgetSim: Create Gadget replicates for use in simulation studies

Description

This function is the almagamation of all the make_gadget_* functions. It allows a user to put together all aspects of a Gadget model for simulation tesing purposes. Note that this function does not currently implement the inclusion of a likelihood file. If a Gadget model optimized to data is more desirable, then you should look to the Rgadget package (specifically gadget.iterative).

Usage

make_gadget_model(...)

Details

This function utilizes the framework of gadgetSim functions to easily assemble the files necessary to simulate a Gadget model. It does not call the model directly, rather it sets up all files desired for inclusion into the Gadget model. You will still; however, have to set up the params files separately.

... should be a number of named elements corresponding to the following names: time, area, stock (or stocks), and fleet (or fleets). Each of these named elements should also correspond to their respective gadget_* classes as produced by the family of make_gadget_* functions (see example below).

Value

A list of class gadget_model. This can then be fed into write_gadget_model

Examples

#------------------------------
# setup time and area
st_year <- 1985
end_year <- 2015
time <- make_gadget_timefile(st_year, end_year, "quarterly")
area <-
  make_gadget_areafile(areas = 1, size = 1e6,
                       temp_data = expand.grid(year = st_year:end_year,
                                               step = 1:4, area = 1,
                                               temp = 3))

#------------------------------
# setup the stock
# setup basic stock information
minage <- 1
maxage <- 10
minlength <- 1
maxlength <- 100
dl <- 1
alpha <- 1e-04
beta <- 3
reflength <- seq(minlength, maxlength, dl)
stock_info <-
    list(stockname = "cod", livesonareas = 1,
         minage = minage, maxage = maxage,
         minlength = minlength, maxlength = maxlength, dl = dl)

# setup refweightfile
stock_refwgt <-
    data.frame(len = reflength,
               weight = alpha * reflength ^ beta)

# setup growth
linf <- 125
k <- 0.15
t0 <- -0.5
stock_growth <-
    list(growthfunction = "lengthvbsimple",
         growthparameters = c(linf, k, alpha, beta))

# setup naturalmortality
stock_m <- rep(0.2, 10)

# setup initial conditions
init_data <-
    normalparamfile(age = seq(minage, maxage, 1),
                    area = 1,
                    age.factor = 1e6 *
                                 exp(-cumsum(rep(0.2,
                                                 (maxage - minage + 1)))),
                    area.factor = 1,
                    mean = vb_formula("cod", minage:maxage,
                                      params = list(linf = linf,
                                                    k = k,
                                                    t0 = t0)),
                    sd = c(1:10, rep(10, 10)),
                    alpha = alpha,
                    beta = beta)
stock_initcond <- list(normalparamfile = init_data)

# setup spawning
stock_spawnfile <-
    make_gadget_spawnfile(
      stockname = "cod",
      start_year = st_year,
      end_year = end_year,
      proportionfunction = c("exponential", -0.25, 35),
      recruitment = bev_holt_formula("cod", params = c(4e08, 1.067e08)),
      stockparameters = c(20, 2, alpha, beta)
    )

# create gadget stockfile
cod <-
   make_gadget_stockfile(stock = stock_info,
                         refweightfile = stock_refwgt,
                         growth = stock_growth,
                         naturalmortality = stock_m,
                         iseaten = 1,
                         initialconditions = stock_initcond,
                         spawning = stock_spawnfile)

#------------------------------
# setup the fleet
lin_flt_data <- expand.grid(year = st_year:end_year,
                            steps = 1:4,
                            area = 1,
                            fleetname = "lin")
lin_flt_data <- lin_flt_data[order(lin_flt_data$year), ]
lin_flt_data$scaling <- c(seq(0.01,0.8, length.out = 62),
                          seq(0.8,0.2,length.out = 62))
lin_fleet <-
    list(type = "linearfleet",
         suitability = make_exponentiall50_suit("lin", "cod"),
         amount = lin_flt_data)
lin <- make_gadget_fleet(lin = lin_fleet)

#------------------------------
# simulate the Gadget model
gad_mod <-
  make_gadget_model(time = time,
                    area = area,
                    stock = cod,
                    fleet = lin)
write_gadget_model(gad_mod, path = "test_model")

# must set up params first, then can run
stock_std <- get_stock_std(path = "test_model")

#------------------------------
# see the results
plot(number ~ year, subset(stock_std$cod, age == 2, step == 1), type = "l")

# if using tidyverse
g <-
  ggplot(data=filter(stock_std$cod, step == 1), aes(x=year, y=number)) +
  geom_line() + facet_wrap(~age, scales = "free_y")

inspktrgadget/gadgetSim documentation built on May 10, 2019, 9:51 a.m.