## Script for stockassessment.org
.libPaths(".")

## **stockassessment.org**
## =======================
## - Avoid installing local versions of at least 'TMB', 'Rcpp', 'RcppEigen'
##   If you make that mistake simply remove your local versions by out-commenting:
##     remove.packages("TMB")
##     remove.packages("Rcpp")
##     remove.packages("RcppEigen")

devtools:::install_github("kaskr/gridConstruct/gridConstruct", upgrade=FALSE) ## upgrade=FALSE just in case...
devtools:::install_github("DTUAqua/limfjord/mussel", upgrade=FALSE) ## upgrade=FALSE => DON'T upgrade dependencies (TMB, Rcpp,...)

## Load packages and data
library(DATRAS)
library(rgdal)
library(sp)
library(mussel)
d <- readExchange("data/data.zip")

## Remove invalid hauls (I/V)
d <- subset(d, HaulVal == "V")

## Select Gear codes
d <- subset(d, Gear %in% c("MSK", "MSKOES"))

## Select Species = "Mytilus edulis" (not yet in DATRAS Worms table)
d <- subset(d, SpecCode==140480)

## Space time subset:
d <- subset(d,lon<9.5 & lat>56)
latestYear <- max(levels(d$Year))
d <- subset(d, Year %in% 2004:as.numeric(latestYear))

## Add swept area in m^2 to dataset:
Width <- c("MSK" = 1.00, "MSKOES" = 0.96) ## FIXME: New gear 1.00 width ?
d$sweptArea <- d$Distance * Width[as.character(d$Gear)]

## Add response variable: Total weight (g)
w <- tapply(d[["HL"]]$CatCatchWgt,d[["HL"]]$haul.id,function(x)x[1])
w[is.na(w)] <- 0
d[["HH"]]$Weight <- w[as.character(d$haul.id)]

## Standardize response to kg/m^2
d[["HH"]]$Weight <- d[["HH"]]$Weight / d[["HH"]]$sweptArea
d$sweptArea[] <- 1

## Print DATRASraw
print(d)

#############################################################
## 3. Add spatial grid
#############################################################

data(grid, package="mussel") ## Get default grid

#############################################################
## 4. Fit model
#############################################################

env <- fitModel(d)

#############################################################
## 5. Save results
#############################################################

save.image("res.RData.exe") ## .exe == don't check in :)