v1.2.4 (#34)

* bug fix re BIs in getToaYaps() & insert mandatory checkInp() in runYaps()

* redo to make getToaYaps() backwards comp

* robustify tmb using bounds - EXPERIMENTAL

* relax pre-flight checks for ping_type='sbi'

* v1.2.3.9005 docs and examples

* reduce track length in example to run faster

* reduce track length in example to run faster

* bump version

* v1.2.4 - CRAN

.Rbuildignore

@@ -15,3 +15,4 @@
 ^codecov\.yml$
 ^\.travis\.yml$
 ^cran-comments\.md$
+^CRAN-RELEASE$

.gitignore

@@ -43,3 +43,5 @@ inst/doc
 debug.log
 
 README.html
+cran-comments.md
+CRAN-RELEASE

.travis.yml

@@ -6,4 +6,3 @@ addons:
   apt:
     packages:
       - libgit2-dev
-

DESCRIPTION

@@ -1,21 +1,23 @@
 Package: yaps
 Title: Track Estimation using YAPS (Yet Another Positioning Solver)
-Version: 1.2.3
+Version: 1.2.4
 Authors@R: c(	person("Henrik", "Baktoft", email = "[email protected]", role = c("cre", "aut"), comment=c(ORCID = "0000-0002-3644-4960")),
-				person("Karl", "Gjelland", role=c("aut")),
-				person("Uffe H.", "Thygesen", role=c("aut")),
-				person("Finn", "Økland", role=c("aut"))
+				person("Karl", "Gjelland", role=c("aut"), comment=c(ORCID = "0000-0003-4036-4207")),
+				person("Uffe H.", "Thygesen", role=c("aut"), comment=c(ORCID = "0000-0002-4311-6324")),
+				person("Finn", "Økland", role=c("aut"), comment=c(ORCID = "0000-0002-1938-5460"))
 			)
 Description: Estimate tracks of animals tagged with acoustic transmitters. 'yaps' was introduced in 2017 as a transparent open-source tool to estimate positions of fish (and other aquatic animals) tagged with acoustic transmitters. Based on registrations of acoustic transmitters on hydrophones positioned in a fixed array, 'yaps' enables users to synchronize the collected data (i.e. correcting for drift in the internal clocks of the hydrophones/receivers) and subsequently to estimate tracks of the tagged animals. The paper introducing 'yaps' is available in open access at Baktoft, Gjelland, Økland & Thygesen (2017) <doi:10.1038/s41598-017-14278-z>. Also check out our cookbook with a completely worked through example at Baktoft, Gjelland, Økland, Rehage, Rodemann, Corujo, Viadero & Thygesen (2019) <DOI:10.1101/2019.12.16.877688>. Additional tutorials will eventually make their way onto the project website at <https://baktoft.github.io/yaps/>.
 Depends: R (>= 3.5.0)
 License: GPL-3
 Encoding: UTF-8
 LazyData: true
+Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.1.1
 LinkingTo: Rcpp, TMB, RcppEigen
 Imports: circular, cowplot, data.table, ggplot2, ggrepel, nloptr, plyr, Rcpp, reshape2, splusTimeSeries, stats, tictoc, TMB, viridis, zoo
 Suggests: 
-    covr,
+    caTools,
+	covr,
     knitr,
     rmarkdown,
     testthat (>= 2.1.0),

NEWS.md

@@ -1,8 +1,21 @@
+# yaps v1.2.4
+
+## New stuff
+* Now on CRAN
+* More checks in checkInp() to catch typical errors in format of inp.
+* EXPERIMENTAL Attempt to robustify runYaps() - use with care.
+
+## Bug fixes
+* Fix bug in getToaYaps() re number of empty pings.
+* Docs and examples fixed to meet requirements.
+* Make getToaYaps() aware of pingType
+
+
 # yaps v1.2.3
 
 ## New stuff
-* Moved example data `hald` to an external package with yaps example data `yapsdata`. Available from github using devtools::install_github('baktoft/yapsdata')
-* Lots of examples and tests added
+* Moved example data `hald` to an external package with yaps example data `yapsdata`. Available from github using devtools::install_github('baktoft/yapsdata').
+* Lots of examples and tests added.
 
 
 # yaps v1.2.2

R/alignBurstSeq.R

@@ -8,8 +8,8 @@
 #' @param seq_lng_min Minimum length of sequence of consecutive pings to use for the alignment. Finds first occurence of sequence of this length in the data and compare to the known burst sequence
 #' @param rbi_min,rbi_max Minimum and maximum burst interval of the transmitter. Used to identify sequence of consecutive pings in the data
 #' @param plot_diag Logical indicating if visual diagnosis plots should be created.
-#' return data.table like input, but with extra columns seq_ping_idx and seq_epo
 #' @export
+#' @return `data.table` like the input `synced_dat`, but with extra columns seq_ping_idx and seq_epo
 #' @example man/examples/example-alignBurstSeq.R
 alignBurstSeq <- function(synced_dat, burst_seq, seq_lng_min=10, rbi_min, rbi_max, plot_diag=TRUE){
 	burst_seq_dt <- data.table::data.table(bi=burst_seq)
@@ -52,12 +52,14 @@ alignBurstSeq <- function(synced_dat, burst_seq, seq_lng_min=10, rbi_min, rbi_ma
 
 	# plot if plot_diag == TRUE
 	if(plot_diag){
+		oldpar <- par(no.readonly = TRUE) 
+		on.exit(par(oldpar))
+
 		par(mfrow=c(1,2))
 		plot(log(seq_diffs))
 		points(log(seq_diffs[seq_fix_idx]) ~ seq_fix_idx, col="red", pch=20, cex=2)
 
 		plot(synced_dat[, eposync - seq_epo] ~ synced_dat$ts, pch=".")
-		par(mfrow=c(1,1))
 	}
 
 	return(synced_dat)

R/applySync.R

@@ -1,10 +1,12 @@
 #' Apply sync model to toa matrix to obtain synced data
+#'
 #' @param toa Object containing data to be synchronized. Typically a `data.table` as e.g. `ssu1$detections`, but can also be a matrix dim=(n_ping, n_hydo).
 #' @param hydros data.table formatted as `ssu1$hydros`
 #' @param sync_model Synchronization model obtained using `getSyncModel()`
-#' @example man/examples/example-yaps_ssu1.R
-
+#'
 #' @export
+#' @return A `data.table` with the now synchronized time-of-arrivals in column `eposync`. 
+#' @example man/examples/example-yaps_ssu1.R
 applySync <- function(toa, hydros="", sync_model){
 	if(is.matrix(toa)) {type <- "toa_matrix"
 	} else if(data.table::is.data.table(toa)) {type <- "detections_table"}

R/checkInp.R

@@ -2,18 +2,27 @@
 #' 
 #' @param inp Object obtained using `getInp()`
 #' @export
+#' @return No return value, but prints errors/warnings if issues with `inp` is detected.
 #' @example man/examples/example-yaps_ssu1.R
 checkInp <- function(inp){
+
+	# check that all BIs are in range of values in the model
+	# only relevant for ping_types 'rbi' and 'pbi'?
+	if(inp$datTmb$pingType != 'sbi'){
+		stopifnot(inp$datTmb$rbi_min <= min(diff(inp$params$top)))
+		stopifnot(inp$datTmb$rbi_max >= max(diff(inp$params$top)))
+	}
 
 	stopifnot(ncol(inp$datTmb$toa) == inp$datTmb$np)
 	stopifnot(nrow(inp$datTmb$toa) == inp$datTmb$nh)
 
+	stopifnot(dim(inp$datTmb$H)[2] == 3)
+
 	# if z_vec != NULL 
 	if(inp$datTmb$how_3d != 'none'){
-		stopifnot(dim(inp$datTmb$H)[2] == 3)
 		stopifnot(length(inp$datTmb$z_vec) == inp$datTmb$np)
 	}
 
-	print("checkInp passed!")
+	print("Pre-flight checkInp() passed!")
 
 }

R/checkInpSync.R

@@ -3,6 +3,7 @@
 #' @inheritParams getInpSync
 #' @param inp_sync Object obtained using `getInpSync()`
 #' @export
+#' @return No return value, but prints errors/warnings if issues with `inp_sync` is detected.
 #' @example man/examples/example-yaps_ssu1.R
 checkInpSync <- function(inp_sync, silent_check){
 	# speed of sound stuff
@@ -34,14 +35,15 @@ checkInpSync <- function(inp_sync, silent_check){
 	} else if(!silent_check & min(sync_coverage$N) < 50) {
 		cat("NOTE: At least one hydro has less than 50 pings in an offset_idx - try getSyncCoverage(inp_sync, plot=TRUE) for visual\n and rerun getInpSync() with increased keep_rate\n")
 	}
-	return(sync_coverage)
 }
 
 #' Quick overview to check if all hydros have enough data within each offset period.
 #'
 #' @inheritParams checkInpSync
 #' @param plot Logical indicating whether to plot a visual or not.
 #' @export
+#' @return A data.table containing number of pings included in each hydro x offset combination.
+#' @example man/examples/example-yaps_ssu1.R
 getSyncCoverage <- function(inp_sync, plot=FALSE){
 	toa <- inp_sync$dat_tmb_sync$toa
 	nh <- ncol(toa)

R/data.R

@@ -19,7 +19,7 @@
 #'        \item ts Timestamp of detection in POSIXct().
 #'        \item tag ID of detected tag.
 #'        \item epo Timestamp as number of seconds since Unix epoch. Can be obtained using as.numeric(ts).
-#'        \item frac Sub-second part of detection timestamp in fractions of second [0-1].
+#'        \item frac Sub-second part of detection timestamp in fractions of second (0-1).
 #'        \item serial Serial number of detecting hydrophone. Must match entry in data.table hydros.
 #'      }
 #'   }

R/fineTuneSyncModel.R

@@ -3,8 +3,8 @@
 #' @param sync_model sync_model obtained using getSyncModel()
 #' @param eps_threshold Maximum value of residual measured in meter assuming speed of sound = 1450 m/s
 #' @param silent logical whether to make getSyncModel() silent
-#' @example man/examples/example-yaps_ssu1.R
 #' @export
+#' @return Fine tuned `sync_model`. See `?getSyncModel` for more info.
 #' @example man/examples/example-yaps_ssu1.R
 fineTuneSyncModel <- function(sync_model, eps_threshold, silent=TRUE){
 	# original inp_sync

R/getBbox.R

@@ -1,12 +1,12 @@
 #' Get a standard bounding box to impose spatial constraints
 #'
 #' Standard is a rectangle based on coordinates of outer hydros +- the buffer in meters
-#' Returns a vector of lenght 6: c(x_min, x_max, y_min, y_max, eps, pen). Limits are given in UTM coordinates.
 #' @param buffer Number of meters the spatial domain extends beyound the outer hydros.
 #' @param eps Specifies how well-defined the borders are (eps=1E-2 is very sharp, eps=100 is very soft). 
 #' @param pen Specifies the penalty multiplier.
 #' @inheritParams getInp
 #' @export
+#' @return Vector of lenght 6: c(x_min, x_max, y_min, y_max, eps, pen). Limits are given in UTM coordinates.
 #' @example man/examples/example-bbox.R
 getBbox <- function(hydros, buffer=5, eps=1E-3, pen=1){
 	x_min <- hydros[which.min(hydros$hx), hx] - buffer

R/getInp.R

@@ -0,0 +1,34 @@
+#' Get prepared inp-object for use in TMB-call
+#'
+#' Wrapper-function to compile a list of input needed to run TMB
+#' @param hydros Dataframe from simHydros() or Dataframe with columns hx and hy containing positions of the receivers. Translate the coordinates to get the grid centre close to (0;0).
+#' @param toa TOA-matrix: matrix with receivers in rows and detections in columns. Make sure that the receivers are in the same order as in hydros, and that the matrix is very regular: one ping per column (inlude empty columns if a ping is not detected).
+#' @param E_dist Which distribution to use in the model - "Gaus" = Gaussian, "Mixture" = mixture of Gaussian and t or "t" = pure t-distribution
+#' @param n_ss Number of soundspeed estimates: one estimate per hour is usually enough
+#' @param pingType Type of transmitter to simulate - either stable burst interval ('sbi'), random burst interval ('rbi') or random burst interval but where the random sequence is known a priori
+#' @param rbi_min,rbi_max Minimum and maximum BI for random burst interval transmitters
+#' @param sdInits If >0 initial values will be randomized around the normally fixed value using rnorm(length(inits), mean=inits, sd=sdInits)
+#' @param ss_data_what What speed of sound (ss) data to be used. Default ss_data_what='est': ss is estimated by the model. Alternatively, if ss_data_what='data': ss_data must be provided and length(ss_data) == ncol(toa)
+#' @param ss_data Vector of ss-data to be used if ss_data_what = 'est'. Otherwise ss_data <- 0 (default)
+#' @param biTable Table of known burst intervals. Only used when pingType == "pbi". Default=NULL
+#' @param z_vec Vector of known depth values (positive real). Default=NULL is which case no 3D is assumed. If z_vec = "est" depth will be estimated.
+#' @param bbox Spatial constraints in the form of a bounding box. See ?getBbox for details.
+
+#' @return List of input data ready for use in `runYaps()`
+#' @export
+#' @example man/examples/example-yaps_ssu1.R
+getInp <- function(hydros, toa, E_dist, n_ss, pingType, sdInits=1, rbi_min=0, rbi_max=0, ss_data_what='est', ss_data=0, biTable=NULL, z_vec=NULL, bbox=NULL){
+	stopifnot(pingType %in% c('sbi', 'pbi', 'rbi'))
+	inp_params 	<- getInpParams(hydros, toa, pingType)
+	datTmb 		<- getDatTmb(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max, ss_data_what, ss_data, biTable, inp_params, z_vec, bbox)
+	params 		<- getParams(datTmb)
+	inits 		<- getInits(datTmb, sdInits)
+	return(list(
+		datTmb = datTmb,
+		params= params,
+		inits = inits,
+		inp_params = inp_params
+		)
+	)
+}
+

R/getInpSync.R

@@ -0,0 +1,75 @@
+#' Get object inp for synchronization
+#'
+#' @param sync_dat List containing data.tables with hydrophone information and detections. See e.g. `?ssu1` for example
+#' @param max_epo_diff Sets the upper threshold for differences in TOA of sync tags. Best parameter value depends on burst rate of sync tags and how far apart the internal clocks of the hydros are prior to synchronization. A bit less than half of minimum sync tag burst rate is a good starting choice.
+#' @param min_hydros Sets the lower threshold of how many hydrophones need to detect each sync tag ping in order to be included in the sync process. Should be as high as possible while observing that all hydrosphones are contributing. If too low, isolated hydrophones risk falling out completely. Future versions will work towards automising this.
+#' @param time_keeper_idx Index of the hydrophone to use as time keeper. Could e.g. be the one with smallest overall clock-drift.
+#' @param fixed_hydros_idx Vector of hydro idx's for all hydrophones where the position is assumed to be known with adequate accuracy and precission. Include as many as possible as fixed hydros to reduce overall computation time and reduce overall variability. As a bare minimum two hydros need to be fixed, but we strongly advice to use more than two.
+#' @param n_offset_day Specifies the number of hydrophone specific quadratic polynomials to use per day. For PPM based systems, 1 or 2 is often adeqaute.
+#' @param n_ss_day Specifies number of speed of sound to estimate per day if no ss data is supplied. It is recommended to use logged water temperature instead. However, estimating SS gives an extra option for sanity-checking the final sync-model.
+#' @param ss_data_what Indicates whether to estimate ("est") speed of sound or to use data based on logged water temperature ("data").
+#' @param ss_data data.table containing timestamp and speed of sound for the entire period to by synchronised. Must contain columns 'ts' (POSIXct timestamp) and 'ss' speed of sound in m/s (typical values range 1400 - 1550).
+#' @param keep_rate Syncing large data sets can take a really long time. However, there is typically an excess number of sync tag detections 
+#'			and a sub-sample is typically enough for good synchronization. 
+#'			This parameter EITHER specifies a proportion (0-1) of data to keep when sub-sampling 
+#'			OR (if keep_rate > 10) number of pings (approximate) to keep in each hydro X offset_idx combination if enough exists.
+#' @param excl_self_detect Logical whether to excluded detections of sync tags on the hydros they are co-located with. Sometimes self detections can introduce excessive residuals in the sync model in which case they should be excluded.
+#' @param lin_corr_coeffs Matrix of coefficients used for pre-sync linear correction. `dim(lin_corr_coeffs)=(#hydros, 2)`. 
+#' @param silent_check Logical whether to get output from `checkInpSync()`. Default is FALSE
+#'
+#' @export
+#' @return List of input data ready for use in `getSyncModel()`
+#' @example man/examples/example-yaps_ssu1.R
+getInpSync <- function(sync_dat, max_epo_diff, min_hydros, time_keeper_idx, fixed_hydros_idx, n_offset_day, n_ss_day, keep_rate=1, excl_self_detect=TRUE, lin_corr_coeffs=NA, ss_data_what="est", ss_data=c(0), silent_check=FALSE){
+	if(length(unique(sync_dat$hydros$serial)) != nrow(sync_dat$hydros)){
+		print(sync_dat$hydros[, .N, by=serial][N>=2])
+		stop("ERROR: At least one hydrophone serial number is used more than once in sync_dat$hydros!\n")
+	}
+
+	if(keep_rate <=0 | (keep_rate > 1 & keep_rate < 10) | (keep_rate >= 10 & keep_rate %% 1 != 0)){
+		stop("ERROR: Invalid keep_rate! Must be either ]0;1] or integer >= 10\n")
+	}
+
+	sync_dat <- appendDetections(sync_dat)
+
+	if(is.na(lin_corr_coeffs[1])){
+		lin_corr_coeffs <- matrix(0, nrow=nrow(sync_dat$hydros), ncol=2, byrow=TRUE)
+	}
+
+	sync_dat <- applyLinCorCoeffsInpSync(sync_dat, lin_corr_coeffs)
+
+	T0 <- min(sync_dat$detections$epo)
+
+	inp_H_info <- getInpSyncHInfo(sync_dat)
+
+	inp_toa_list_all		<- getInpSyncToaList(sync_dat, max_epo_diff, min_hydros, excl_self_detect)
+	fixed_hydros_vec 		<- getFixedHydrosVec(sync_dat, fixed_hydros_idx)
+	offset_vals_all			<- getOffsetVals(inp_toa_list_all, n_offset_day)
+	inp_toa_list 			<- getDownsampledToaList(inp_toa_list_all, offset_vals_all, keep_rate)
+	offset_vals				<- getOffsetVals(inp_toa_list, n_offset_day)
+	ss_vals 				<- getSsVals(inp_toa_list, n_ss_day)
+	if(ss_data_what == "data"){
+		ss_data_vec		<- getSsDataVec(inp_toa_list, ss_data)
+	} else {
+		ss_data_vec <- c(0)
+	}
+
+	dat_tmb_sync 		<- getDatTmbSync(sync_dat, time_keeper_idx, inp_toa_list, fixed_hydros_vec, offset_vals, ss_vals, inp_H_info, T0, ss_data_what, ss_data_vec)
+	params_tmb_sync 	<- getParamsTmbSync(dat_tmb_sync, ss_data_what)
+	if(ss_data_what == "est"){
+		random_tmb_sync <- c("TOP", "OFFSET", "SLOPE1", "SLOPE2", "SS", "TRUE_H")
+	} else {
+		random_tmb_sync <- c("TOP", "OFFSET", "SLOPE1", "SLOPE2", "TRUE_H")
+	}
+	# inits_tmb_sync <- c(3, rep(-3,dat_tmb_sync$nh))
+	inits_tmb_sync <- c(-3)
+	inp_params <- list(toa=inp_toa_list$toa, T0=T0, Hx0=inp_H_info$Hx0, Hy0=inp_H_info$Hy0, offset_levels=offset_vals$offset_levels, 
+		ss_levels=ss_vals$ss_levels, max_epo_diff=max_epo_diff, hydros=sync_dat$hydros,
+		lin_corr_coeffs=lin_corr_coeffs, min_hydros=min_hydros, ss_data=ss_data
+	)
+
+	inp_sync <- list(dat_tmb_sync=dat_tmb_sync, params_tmb_sync=params_tmb_sync, random_tmb_sync=random_tmb_sync, inits_tmb_sync=inits_tmb_sync, inp_params=inp_params)
+	inp_sync$inp_params$sync_coverage <- checkInpSync(inp_sync, silent_check)
+	return(inp_sync)
+
+}