Merged develop into master

.gitignore

@@ -0,0 +1,14 @@
+html/
+src/build/*
+src/equilibrium.obj
+src/gmacs.cpp
+src/gmacs.exe
+src/gmacs.htp
+src/gmacs.obj
+src/moltIncrement.cpp
+src/multinomial.obj
+src/nloglike.obj
+src/robust_multi.obj
+src/spr.obj
+src/tailcompression.obj
+.Rproj.user

CHANGES.md

@@ -0,0 +1,54 @@
+# Gmacs Changes
+## List of changes to Gmacs since pilot version (V1.0):
+
+1. 	Generic functions at end of code, now exported to Cstar.
+2. 	Selectivity functions exported to Cstar.
+3. 	Cstar selectivity functions used in code (more options available).
+4. 	When number of classes in data and model are same, class link matrix not required. Automatically generated instead.
+5.	When number of classes in the data and the model differ by an integer factor, class link matrix is automatically generated. 
+6.	When number of classes in the data is not a multiple of number of classes in the model, read in class_link matrix.
+7.	Names of fleets and surveys now printed correctly to echoinput file: (bug fix).
+8.	Functionality improved for writeR: Fleet and survey names now exported to R and used for plots.
+9.	All references in model are now to size-classes rather than length classes.
+10. Gmacs R Functions now wrapped in simple 'gmplot_all' function, can be used with any model.
+11.	Fleet control section of data file is now extended to include surveys, catch units and multipliers, all entered in the one matrix. This is to prepare for future options where a fishing fleet might also have a 'survey' such as a CPUE index, and where a survey might have some catch.
+12.	Data file now expects extra dimensions to be specified, such as number of shell and maturity types.
+13. Initial numbers can now be specified by estimating early recruitments. Initial numbers options input via control file. NOTE: Currently the lognin_parms have to be entered and have phase set to -ve so as not to be estimated. Later: Make reading these lines conditional.
+14. Retention function can now be selected from among multiple options. This includes a logistic function. Currently borrowed from cstar::selex functions.
+15.	Growth functions can now be selected from among multiple options. This includes a linear growth relationship with a gamma distribution about each size class.
+16. Internal calculations have been modified so that multiple copies of selectivity, retention, or size-transition matrix patterns are not created nor stored as before. This make the code mode efficient. NOTE: Will do the same for selectivity, but waiting until further selectivity updates are made [to selex_fleet and selex_survey]
+
+* Cstar functions for selectivity (and some others) have now been written in Template style code.
+* Many updates to Gmacs R functions: Can now be used for any Gmacs model. 
+
+---
+
+## Priority changes for Gmacs: Development of a complete example for BBRKC.
+
+1. A 20 size class model requires several changes from the basic structure presented in January:
+  * Change to available selectivity functions, beyond 'parameter-per-class'. DONE
+  * Change to initial numbers estimation from 'parameter-per-class' option to early-recruitment build-up option. DONE
+  * Change to available retention functions, beyond 'parameter-per-class'. MODIFY FROM SELECTIVITY
+  * Change to growth estimation from 'parameter-per-class' to parametric approach. 
+2. Change population dynamics calculations to include more dimensions:
+  * Set-up dimensions of N matrix via input numbers of sex, shell, maturity. DONE
+  * Read in data for each dimension as necessary: DONE
+  * Make sure predicted and observed vales have the same generalized structure and complete calculations. 
+  * Change pop dynamics equations to loop over different dimensions:
+3. Remove penalty E3 in the final estimation phase.
+
+4. Include molting probability (inc. time-varying) into the calculation of the growth transition matrix.
+  * Might be able to cheat at start by using growth transition matrix from Jie.
+
+5. Create output file with estimated parameters, starting values, final estimates, bounds, phases, as well as asymptotic standard errors. 
+  * Automatically highlight cases where parameters are on or close to the bounds in the input.
+
+6. Update R plots for Gmacs Assessment Report 
+  * Add confidence intervals to plots of data points which are considered uncertain. 
+    - DONE for Recruitment estimates
+    - Others with SD functions?
+
+7. Check weightings and likelihoods are working correctly with more generalised model.
+
+8. See NPRB model for 'biological parameters' set up, which uses a counter to allow for differing numbers of parms. 
+  This could replace the current 'theta' block of parameters.

README.md

@@ -1,36 +1,25 @@
-# Gmacs Version 1.0 #
+# Gmacs
 
-This is the pilot release of Gmacs. Currently posted source files are compilable using ADMB 11.1 and have been tested using the BBRKC model available in the examples folder. This release will remain active until the current 'under development' version is released. **Updated February 2014, by Athol Whitten**
+Gmacs is currently under development. A simple working release version of Gmacs is available via `Tag V1.0` and has been tested using the BBRKC model available in the examples folder. The next major release of Gmacs is planned for September 2014.
 
-## Generalized Modeling for Alaskan Crab Stocks ##
+## Table of contents
+- [About Gmacs](#about-gmacs)
+- [Gmacs R Package](#r-package-for-gmacs)
+- [Development](#development)
 
-This repository holds source code, instructions, examples, and associated scripts for **Gmacs** (Generalized Modeling for Alaskan Crab Stocks), a generic size-based stock assessment model. 
+## About Gmacs
+**Gmacs** is a generalized size-structured stock assessment modelling framework. The framework is designed with similar flexibility to that provided by age-structured stock assessment modelling frameworks like Stock Synthesis and CASAL. Gmacs can fit to a wide-variety of data for single sex or sex-specific population dynamics and fishery models: data can include survey and fishery indices of abundance and fishery- and survey-based size-composition data.
 
-## Modeling structure and format ##
+### Data Requirements
+Data must be supplied via the `model.dat` file in a *flat format* to enable easy indexing and simple preparation. Each record for catch, abundance, length-structure etc. should be held in an individual row, with information relating to year, fleet, sex and more.
 
-Gmacs implements a size-structured modelling framework with flexibility similar to that provided by other general stock assessment modelling frameworks. Some effort has been made to maintain consistency with data and control file formats familiar to users of Stock Synthesis.
+Model specifcations are controlled through the `model.ctl` file. Information read from files is printed to a separate file called `echoinput.rep` allowing users to check and debug their input. For more information see the [Gmacs Wiki](https://github.com/seacode/gmacs/wiki).
 
-### Input file structure
+## R Package for Gmacs
+An R package, called `gmr` is under development in support of Gmacs. The package provides functions for creating plots from Gmacs output files. A full pilot version is intended for release in September 2014, timed to coincide with the next stable release of Gmacs. Current development versions of the package can be downloaded from Github directly through R, see https://github.com/seacode/gmr for more details.
 
-Data are supplied via the `model.dat` file in a *flat format* to enable easy indexing and simple preparation using spreadsheet software. Each record for catch, abundance, length-structure etc. should be held in an individual row, with information relating to year, fleet, sex and more:
+## Simulation Mode
+A simulation-estimation procedure can be performed with Gmacs, by using the `gmacs -sim` flag. For example, try `gmacs -sim 123`, where 123 is a random number seed.
 
-####  Catch data structure
-
-  * Year, Season, Fleet, Sex, Observation    
-
-####  Survey data structure
-
-  * Year, Season, Survey, Observation, Error
-
-####  Length frequency data structure  
-
-  * Year, Season, Fleet/Survey, Sex, Maturity, Shell Condition, No. Samples, Data
-
-Gmacs allows for the inclusion of an optional growth data file `growth.dat` to specify a fixed growth transtion matrix or year-specific growth transtion matrices. The program also reads a `starter.gm` file for specifying the overall model run conditions, and a control file `model.ctl` for specifications relating to parameter estimation. Finally, a `forecast.gm` file is read to specify the calculation of relevant reference points. This file will allow users to specify model projection options in later versions of Gmacs.
-
-During the read-in procedure, helpful messages are printed to screen and the information read in is printed to a separate file called `echoinput.gm` allowing users to check and debug their data and control files. 
-
-A general user-guide to the program is under development and will be made available with future releases.
-
-## Development ##
-This software is under development and is not yet intended for general use. If you would like to contribute to the project, please contact [Athol Whitten](mailto:[email protected]). 
+## Development
+This software is under development and is not yet intended for general use. If you would like to contribute to the project, please see the [Developers Guide](https://github.com/seacode/gmacs/wiki/5.-Developers).

Rsrc/DESCRIPTION

@@ -0,0 +1,19 @@
+Package: gmr
+Title: Plotting and analysis tools for the Gmacs stock assessment framework
+Description: gmr is a set of tools for analysing data and outputs from Gmacs
+    stock assessment models. The package streamlines the process of taking
+    text-based ADMB output files and creating visual plots of both input data
+    and fitted model results.
+Authors@R: c(person("Whitten", "Athol", email = "[email protected]", role
+    = c("aut","cre")), person("Ianelli", "Jim", email = "[email protected]",
+    role = c("aut","cre")), person("Martell", "Steve", role = "aut"))
+URL: https://github.com/seacode/gmr
+Version: 0.2
+Depends:
+    R (>= 3.0.0),
+    ggplot2 (>= 0.9.3.1)
+Imports:
+    gdata (>= 2.13.3),
+    reshape2 (>= 1.2.2),
+    shiny (>= 0.10.2.2),
+License: MIT

Rsrc/NAMESPACE

@@ -0,0 +1,31 @@
+# Generated by roxygen2 (4.0.2): do not edit by hand
+
+export(get_cpue)
+export(get_recruitment)
+export(get_selectivity)
+export(get_sizecomp)
+export(get_ssb)
+export(plot_catch)
+export(plot_cpue)
+export(plot_cpue_res)
+export(plot_datarange)
+export(plot_growth)
+export(plot_growth_inc)
+export(plot_multiple)
+export(plot_naturalmortality)
+export(plot_recruitment)
+export(plot_selectivity)
+export(plot_sizecomp)
+export(plot_sizecomp_res)
+export(plot_sizetransition)
+export(plot_ssb)
+export(read_admb)
+export(read_fit)
+export(read_psv)
+export(read_rep)
+export(set_ggtheme)
+export(shiny_gmacs)
+import(gdata)
+import(ggplot2)
+import(reshape2)
+import(shiny)

Rsrc/R/get-cpue.R

@@ -0,0 +1,16 @@
+#' Plot cpue or other indices
+#'
+#' @param replist List object created by read_admb function
+#' @return dataframe of observed and predicted indices and residuals
+#' @export
+get_cpue <- function(replist){
+	A <- replist
+	df <- as.data.frame(A$dSurveyData)
+	colnames(df) <- c("year","seas","fleet","sex","cpue","cv","units")
+	sd <- sqrt(log(1+df$cv^2))
+	df$lb <- exp(log(df$cpue)-1.96*sd)
+	df$ub <- exp(log(df$cpue)+1.96*sd)
+	df$pred <- na.exclude(as.vector(t(A$pre_cpue)))
+	df$resd <- na.exclude(as.vector(t(A$res_cpue)))
+  return(df)
+}

Rsrc/R/get-recruitment.R

@@ -0,0 +1,20 @@
+#' Plot predicted recruitment and approximate asymptotic error-bars
+#'
+#'
+#' @param replist List object created by read_admb function
+#' @return Dataframe of recruitment
+#' @export
+get_recruitment <- function(replist){
+  A <- replist
+  if(is.null(A$fit$logDetHess)) {
+      stop("Appears that the Hessian was not positive definite\n
+            thus estimates of recruitment do not exist.\n
+            See this in replist$fit.")
+  }
+  dfpar   <- data.frame(par=A$fit$names,log_rec=A$fit$est,log_sd=A$fit$std)
+  df      <- subset(dfpar,par=="sd_log_recruits")[,-1]
+  df$year <- A$mod_yrs
+  df$lb   <- exp(df$log_rec - 1.96*df$log_sd)
+  df$ub   <- exp(df$log_rec + 1.96*df$log_sd)
+  return(df)
+}

Rsrc/R/get-ssb.R

@@ -0,0 +1,16 @@
+#' Plot predicted spawning stock biomass (ssb)
+#'
+#' Spawning biomass may be defined as all males or some combination of males and females
+#'
+#' @param replist List object created by read_admb function
+#' @return Dataframe of spawning biomass
+#' @export
+get_ssb <- function(replist){
+  A <- replist
+  dfpar   <- data.frame(par=A$fit$names,log_mmb=A$fit$est,log_sd=A$fit$std)
+  df      <- subset(dfpar,par=="sd_log_mmb")[,-1]
+  df$year <- A$mod_yrs
+  df$lb   <- exp(df$log_mmb - 1.96*df$log_sd)
+  df$ub   <- exp(df$log_mmb + 1.96*df$log_sd)
+  return(df)
+}

Rsrc/R/gmr.R

@@ -0,0 +1,13 @@
+#' gmr: R code for Gmacs
+#'
+#' gmr is a set of tools for analysing data and outputs
+#' related to Gmacs stock assessment models.
+#'
+#'
+#' @name gmr
+#' @docType package
+#' @import ggplot2
+#' @import shiny
+#' @import reshape2
+#' @import gdata
+NULL

Rsrc/R/plot-catch.R

@@ -0,0 +1,42 @@
+#' Plot observed and predicted catch values
+#'
+#' @param replist List object created by read_admb function
+#' @param plot_res plot residuals only (default=F)
+#' @return Plot of catch history (observed) and predicted values
+#' @export
+plot_catch <- function(replist, plot_res=FALSE)
+{
+  A  <- replist
+  df <- as.data.frame(A$dCatchData)
+  colnames(df)<- c("year","seas","fleet","sex","obs","cv","type","units","mult","effort")
+  df$residuals <- na.omit(as.vector(t(A$res_catch)))
+
+  #Loop over retained and discarded catch.
+  type = unique(df$type)
+  ldf  = list()
+  for(i in type)
+  {
+    ldf[[i]] <- subset(df,type %in% i)
+  }
+  if (plot_res)
+  {
+    # Residuals
+    p <- ggplot(df,aes(x=factor(year),y=residuals,fill=factor(sex)))
+    p <- p + geom_bar(stat = "identity", position="dodge")
+    p <- p + scale_x_discrete(breaks=pretty(df$year))
+    p <- p + labs(x="Year",y="Residuals ln(kt)",fill="Sex")
+    p <- p + facet_wrap(~fleet~type,scales="free")
+  }
+  else
+  {
+    p <- ggplot(df,aes(x=factor(year),y=obs,fill=factor(sex)))
+    p <- p + geom_bar(stat = "identity")
+    p <- p + scale_x_discrete(breaks=pretty(df$year))
+    p <- p + labs(x="Year",y="Catch (kt)",fill="Sex")
+    p <- p + facet_wrap(~fleet,scales="free")
+  }
+  # This line applies the plotting over all unique types...
+  pCatch <- lapply(ldf,FUN = function(x,p){p %+% x},p=p)
+  pCatch <- p + ggtheme
+  return(pCatch)
+}

Rsrc/R/plot-cpue.R

@@ -0,0 +1,32 @@
+#' Plot cpue or other indices
+#'
+#' @param replist List object created by read_admb function
+#' @return Plot of all observed and predicted incices
+#' @export
+plot_cpue <- function(replist){
+	df <- get_cpue(replist)
+	p  <- ggplot(df,aes(year,cpue))
+# p  <- p + geom_point(aes(col=sex))
+	p  <- p + geom_pointrange(aes(year,cpue,ymax=ub,ymin=lb,col=sex))
+	p  <- p + labs(x="Year",y="CPUE",col="Sex")
+	pCPUE  <- p + facet_wrap(~fleet+sex,scales="free")
+# Fitted CPUE
+	pCPUEfit <- pCPUE + geom_line(data=df,aes(year,pred))
+  return(pCPUEfit)
+}
+
+#' Plot residuals of cpue or other indices
+#'
+#' @param replist List object created by read_admb function
+#' @return Plot of fit indices residuals
+#' @export
+plot_cpue_res <- function(replist){
+# CPUE residuals
+	df <- get_cpue(replist)
+	p  <- ggplot(df,aes(factor(year),resd))
+	p  <- p + geom_bar(aes(fill=factor(sex)),stat = "identity", position="dodge")
+	p  <- p + scale_x_discrete(breaks=pretty(df$year))
+	p  <- p + labs(x="Year",y="CPUE Residuals",fill="Sex")
+	pCPUEres  <- p + facet_wrap(~fleet,scales="free_x")
+  return(pCPUEres)
+}