too many adjustments...

DESCRIPTION

@@ -12,9 +12,9 @@ License: GPL-3
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 LinkingTo: Rcpp, TMB, RcppEigen
-Imports: circular, cowplot, data.table, ggplot2, ggrepel, nloptr, plyr, Rcpp, reshape2, splusTimeSeries, stats, TMB, viridis, zoo
+Imports: circular, cowplot, data.table, ggplot2, ggrepel, nloptr, plyr, Rcpp, reshape2, stats, TMB, viridis, zoo
 Suggests: 
     caTools,
 	covr,

NAMESPACE

@@ -4,6 +4,7 @@ export(alignBurstSeq)
 export(applyLinCor)
 export(applySync)
 export(applySync_old)
+export(checkHydros)
 export(checkInp)
 export(checkInpSync)
 export(checkInpSyncData)
@@ -29,6 +30,7 @@ export(getSyncCheckDat)
 export(getSyncCoverage)
 export(getSyncModel)
 export(getSyncToa)
+export(getToaRbi)
 export(getToaYaps)
 export(plotBbox)
 export(plotSyncCheck)
@@ -41,6 +43,7 @@ export(plotSyncNetwork)
 export(plotYaps)
 export(prepDetections)
 export(prepSyncParams)
+export(removeMultipath)
 export(runTmb)
 export(runYaps)
 export(simHydros)

R/checkHydros.R

@@ -0,0 +1,21 @@
+#' Internal function to check hydros data.tables
+#' @export
+checkHydros <- function(hydros){
+	if(is.null(attr(hydros, 'yapsified'))){
+		cat("ERROR: The hydros data.table is not yapsified. Run e.g. hydros <- yapsify(hydros) to do so. \n")
+		stopSilent()
+	}
+
+	if(!is.null(attr(hydros, 'yapsified')) & attr(hydros, 'yapsified') == FALSE){
+		cat("ERROR: The hydros data.table is not yapsified. Run e.g. hydros <- yapsify(hydros) to do so. \n")
+		stopSilent()
+	}
+
+	if(length(unique(hydros$h_sn)) != nrow(hydros)){
+		cat("ERROR: At least one hydrophone serial number is used more than once in sync_dat$hydros!\n")
+		stopSilent()
+	}
+
+
+
+}

R/checkInp.R

@@ -10,13 +10,17 @@ checkInp <- function(inp){
 	# only relevant for ping_types 'rbi' and 'pbi'?
 	if(inp$datTmb$pingType != 'sbi'){
 		if(inp$datTmb$rbi_min > min(diff(inp$params$top)))	{
-			print(paste0("WARNING: inp$datTmb$rbi_min > min(diff(inp$params$top)) | ",inp$datTmb$rbi_min," > ",min(diff(inp$params$top)),""))
-			inp$datTmb$rbi_min <- min(diff(inp$params$top)) * 0.90
+			cat(paste0("ERROR: inp$datTmb$rbi_min > min(diff(inp$params$top)) | ",inp$datTmb$rbi_min," > ",min(diff(inp$params$top)),"\n"))
+			# inp$datTmb$rbi_min <- min(diff(inp$params$top)) * 0.90
+			# cat("...inp$datTmb$rbi_min adjusted to ", inp$datTmb$rbi_min, "\n")
+			stopSilent()
 		}
 
 		if(inp$datTmb$rbi_max < max(diff(inp$params$top))){
-			print(paste0("WARNING: inp$datTmb$rbi_max < max(diff(inp$params$top)) | ",inp$datTmb$rbi_max," < ",max(diff(inp$params$top)),""))
-			inp$datTmb$rbi_max <- max(diff(inp$params$top)) * 1.10
+			cat(paste0("ERROR: inp$datTmb$rbi_max < max(diff(inp$params$top)) | ",inp$datTmb$rbi_max," < ",max(diff(inp$params$top)),"\n"))
+			# inp$datTmb$rbi_max <- max(diff(inp$params$top)) * 1.10
+			# cat("...inp$datTmb$rbi_max adjusted to ", inp$datTmb$rbi_max, "\n")
+			stopSilent()
 		}
 	}
 
@@ -40,7 +44,4 @@ checkInp <- function(inp){
 	}
 
 	print("Pre-flight checkInp() passed!")
-
-	return(inp)
-
 }

R/checkInpSync.R

@@ -56,33 +56,3 @@ checkInpSync <- function(inp_sync, silent_check){
 
 }
 
-#' 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)
-	offset_idx <- inp_sync$dat_tmb_sync$offset_idx
-
-	toa_long <- data.table::data.table(reshape2::melt(toa))
-	colnames(toa_long) <- c('ping', 'h','toa')
-	toa_long[, offset_idx := rep(offset_idx, nh)]
-	sync_coverage <- data.table::data.table(reshape2::melt(with(toa_long[!is.na(toa)], table(h, offset_idx))))
-	colnames(sync_coverage) <- c('h', 'offset_idx' ,'N')
-
-	if(plot){
-		p <- ggplot2::ggplot(sync_coverage) 
-		p <- p + geom_point(aes(offset_idx, N), col="steelblue") 
-		p <- p + geom_point(data=sync_coverage[N < 50], aes(offset_idx, N), col="blue", size=2)
-		p <- p + geom_point(data=sync_coverage[N < 10], aes(offset_idx, N), col="orange", size=2)
-		p <- p + geom_point(data=sync_coverage[N <= 5], aes(offset_idx, N), col="red", size=2)
-		p <- p + facet_wrap(~h) + ylim(0, max(sync_coverage$N))
-		print(p)
-	}
-
-	return(sync_coverage)
-}

R/getBounds.R

@@ -0,0 +1,54 @@
+#' Get bounds restricting the optimizer
+#* 
+#' Compile a matrix of lower (bounds[,1]) and upper (bounds[,2]) bounds for the parameters to be estimated.
+#' @param datTmb Object obtained using getDatTmb()
+#' @return Matrix of bounds restricting the optimizer when running runYaps().
+#' @noRd
+getBounds <- function(datTmb) {
+	lu_logD_xy 			<- c(-50,  2)
+	lu_logD_z 			<- c(-50,  2)
+
+	lu_logSigma_toa 	<- c(-12, -2)
+	if(datTmb$Edist[2] == 1){ # mixture
+		lu_logScale 	<- c(-30, 10)
+	} else if (datTmb$Edist[3] == 1) { # t
+		lu_logScale 	<- c(-10,2)
+	}
+	lu_log_t_part 		<- c(-100, 100)
+	if(datTmb$pingType == 'rbi'){
+		lu_logSigma_bi 		<- c(-20, 20)
+	} else {
+		lu_logSigma_bi 		<- c(-20, -2)
+	}
+	lu_logD_v 			<- c(-20,  2)
+
+	bounds <- c()
+	bounds <- rbind(bounds, lu_logD_xy)
+
+	if(datTmb$how_3d == 'est'){
+		bounds <- rbind(bounds, lu_logD_z)
+	}
+
+	if(datTmb$ss_data_what == 'est'){
+		bounds <- rbind(bounds, lu_logD_v)
+	}
+
+	if(datTmb$Edist[1] == 1){
+		bounds <- rbind(bounds, lu_logSigma_toa)
+	} else if(datTmb$Edist[2] == 1){
+		bounds <- rbind(bounds, lu_logSigma_toa, lu_logScale, lu_log_t_part)
+	} else if(datTmb$Edist[3] == 1){
+		bounds <- rbind(bounds, lu_logScale)
+	}
+
+	if(datTmb$pingType == 'sbi'){
+		bounds <- rbind(bounds, lu_logSigma_bi)
+	} else if (datTmb$pingType == 'rbi'){
+		bounds <- rbind(bounds, lu_logSigma_bi)
+	} else if (datTmb$pingType == 'pbi'){
+		bounds <- rbind(bounds, lu_logSigma_bi)
+	}
+
+	return(bounds)
+}
+

R/getDatTmb.R

@@ -0,0 +1,83 @@
+#' Internal function - get data for input to TMB
+#'
+#' Compile data for input to TMB.
+#' @param inp_params Selection of parameters used to setup and run YAPS.
+#' @inheritParams getInp
+#'
+#' @return List for use in TMB.
+#' @noRd
+getDatTmb <- function(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max, ss_data_what, ss_data, biTable, inp_params, z_vec, bbox){
+	T0 <- inp_params$T0
+	Hx0 <- inp_params$Hx0
+	Hy0 <- inp_params$Hy0
+
+	toa <- toa - T0
+
+	# allowing slight out-of-bounds BIs
+	rbi_min <- rbi_min - rbi_min * 0.05
+	rbi_max <- rbi_max + rbi_max * 0.05
+
+	# attempting to make sure toa is oriented correct
+	if(!nrow(toa) == nrow(hydros)){
+		toa <- t(toa)
+	}
+
+	if(n_ss > 1){
+		ss_idx <- cut(1:ncol(toa), n_ss, labels=FALSE) - 1 #-1 because zero-indexing in TMB
+	} else {
+		ss_idx <- rep(0, ncol(toa))
+	}
+	approxBI <- mean(diff(colMeans(toa, na.rm=TRUE), na.rm=TRUE), na.rm=TRUE)
+
+	if(ss_data_what == 'data') { stopifnot(length(ss_data) == ncol(toa))}
+
+	Edist <- rep(0,3)
+	if(E_dist == "Gaus") {Edist[1] <- 1}
+	if(E_dist == "Mixture") {Edist[2] <- 1}
+	if(E_dist == "t") {Edist[3] <- 1}
+
+	if(is.null(z_vec)){
+		how_3d <- 'none'
+		z_vec <- c(1)
+	} else if(z_vec[1] == "est") {
+		how_3d <- 'est'
+		z_vec <- c(1)
+	} else {
+		how_3d <- 'data'
+	}
+
+	if(is.null(bbox)){
+		bbox <- NA
+	} else {
+		bbox[1] <- bbox[1] - inp_params$Hx0
+		bbox[2] <- bbox[2] - inp_params$Hx0
+		bbox[3] <- bbox[3] - inp_params$Hy0
+		bbox[4] <- bbox[4] - inp_params$Hy0
+	}
+
+	datTmb <- list(
+		model = "yaps_track",
+		H = matrix(c(hydros$h_x-Hx0, hydros$h_y-Hy0, hydros$h_z), ncol=3),
+		nh = nrow(hydros),
+		np = ncol(toa),
+		Edist = Edist,
+		toa = toa,
+		bi_epsilon = 1E-6,
+		bi_penalty = 1E9,
+		rbi_min = rbi_min,
+		rbi_max = rbi_max,
+		pingType = pingType,
+		n_ss = n_ss,
+		ss_idx = ss_idx,
+		ss_data_what = ss_data_what,
+		ss_data = ss_data,
+		approxBI = approxBI,
+		biTable = c(1),
+		how_3d = how_3d,
+		z_vec = z_vec,
+		bbox = bbox
+	)
+	if(pingType == 'pbi') {datTmb$biTable = biTable}
+
+	return(datTmb)
+}

R/getDistMat.R

@@ -6,7 +6,7 @@ getDistMat <- function(hydros){
 	colnames(dist_mat) <- hydros$h_sn
 	for(hx in 1:nrow(hydros)){
 		for(hy in 1:nrow(hydros)){
-			dist_mat[hx, hy] <- sqrt((hydros[hx, x] - hydros[hy, x] )^2 + (hydros[hx, y] - hydros[hy, y] )^2 + (hydros[hx, z] - hydros[hy, z] )^2)
+			dist_mat[hx, hy] <- sqrt((hydros[hx, h_x] - hydros[hy, h_x] )^2 + (hydros[hx, h_y] - hydros[hy, h_y] )^2 + (hydros[hx, h_z] - hydros[hy, h_z] )^2)
 		}
 	}
 	return(dist_mat)

R/getInits.R

@@ -0,0 +1,59 @@
+#' Get inits for use in TMB
+#'
+#' Compile a vector of initial values to use in TMB. One value for each estimated parameter (not random effects).
+#' Should all be in a credible range.
+#' @param datTmb Object obtained using getDatTmb()
+#' @inheritParams getInp
+#' @return Vector of initial values to use in TMB
+#' @noRd
+getInits <- function(datTmb, sdInits=1) {
+	init_logD_xy <- -1
+
+	if(datTmb$pingType == 'sbi') {
+		init_logSigma_bi <- -6
+	} else if(datTmb$pingType == 'rbi'){
+		if(datTmb$rbi_max >= 10){
+			init_logSigma_bi <- 4 
+		} else {
+			init_logSigma_bi <- -3
+		}
+	} else if(datTmb$pingType == 'pbi'){
+		init_logSigma_bi <- -5
+	}
+
+	init_logD_v <- 0
+	init_logSigma_toa <- -3 # used in Gaussian and mixture
+	init_logScale <- 1		# used in mixture and pure t
+	init_log_t_part <- -4	# only used in mixture
+
+	inits <- c(init_logD_xy)
+
+	if(datTmb$how_3d == 'est'){
+		init_logD_z <- 0
+		inits <- c(inits, init_logD_z)
+	}
+
+	if(datTmb$ss_data_what == 'est'){
+		inits <- c(inits, init_logD_v)
+	}
+
+
+	if(datTmb$Edist[1] == 1){
+		inits <- c(inits, init_logSigma_toa)
+	} else if(datTmb$Edist[2] == 1){
+		inits <- c(inits, init_logSigma_toa, init_logScale, init_log_t_part)
+	} else if(datTmb$Edist[3] == 1){
+		inits <- c(inits, init_logScale)
+	}
+
+	if(datTmb$pingType == 'sbi'){
+		inits <- c(inits, init_logSigma_bi)#,  init_logD_v)#, init_logSigma_toa, init_logScale, init_log_t_part)
+	} else if (datTmb$pingType == 'rbi'){
+		inits <- c(inits, init_logSigma_bi)#, 					init_logD_v)#, init_logSigma_toa, init_logScale, init_log_t_part)
+	} else if (datTmb$pingType == 'pbi'){
+		inits <- c(inits, init_logSigma_bi)#,  init_logD_v)#, init_logSigma_toa, init_logScale, init_log_t_part)
+	}
+
+	inits <- stats::rnorm(length(inits), mean=inits, sd=sdInits)
+	return(inits)
+}

R/getInp.R

@@ -19,6 +19,9 @@
 #' @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'))
+
+	checkHydros(hydros)
+
 	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)

R/getInpParams.R

@@ -0,0 +1,14 @@
+#' Get parameters for this specific data set
+#'
+#' Compile a list of relevant parameters (e.g. T0) to use later on
+#' @inheritParams getInp
+#' @noRd
+getInpParams <- function(hydros, toa, pingType){
+	T0 <- min(toa, na.rm=TRUE)
+
+	Hx0 <- hydros[1,h_x]
+	Hy0 <- hydros[1,h_y]
+
+	return(list(T0=T0, Hx0=Hx0, Hy0=Hy0))
+
+}