Merge pull request #36 from baktoft/v1.2.4.9000

V1.2.5

.github/workflows/check-standard.yaml

@@ -17,8 +17,8 @@ jobs:
           - {os: windows-latest, r: '3.6'}
           - {os: windows-latest, r: 'release'}
           - {os: windows-latest, r: 'devel'}
-          - {os: macOS-latest, r: '3.6'}
           - {os: macOS-latest, r: 'release'}
+          - {os: macOS-latest, r: 'devel'}
           - {os: ubuntu-20.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
           - {os: ubuntu-20.04, r: 'devel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
 

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: yaps
 Title: Track Estimation using YAPS (Yet Another Positioning Solver)
-Version: 1.2.4
+Version: 1.2.5
 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"), comment=c(ORCID = "0000-0003-4036-4207")),
 				person("Uffe H.", "Thygesen", role=c("aut"), comment=c(ORCID = "0000-0002-4311-6324")),

NEWS.md

@@ -1,3 +1,15 @@
+# yaps v1.2.5
+
+## New stuff
+* Add support for ping_type = 'sbi_doulbe'. Special case - needs carefull contrsuction of TOA-matrix.
+* More robustification of the optimizer.
+* More pre-flight checks to catch issues with inp.
+* Use separate diffusivity for Z when estimating 3D
+
+## Bug fixes
+* Bug fix - constraint hitting very fast transmitters with BI < 1 relaxed
+* Implement better start of top sequence when ping_type = sbi_double
+
 # yaps v1.2.4
 
 ## New stuff

R/checkInp.R

@@ -13,13 +13,22 @@ checkInp <- function(inp){
 		stopifnot(inp$datTmb$rbi_max >= max(diff(inp$params$top)))
 	}
 
+	# check correct dimensions...
 	stopifnot(ncol(inp$datTmb$toa) == inp$datTmb$np)
 	stopifnot(nrow(inp$datTmb$toa) == inp$datTmb$nh)
+	stopifnot(dim(inp$datTmb$H)[2] == 3)
+
+	# check for correct number of lower and upper bounds
+	stopifnot(nrow(inp$bounds) == length(inp$inits))
+
+	# check that neither first nor last row in toa is all NAs
+	stopifnot(sum(!is.na(inp$datTmb$toa[,1])) > 0)
+	stopifnot(sum(!is.na(inp$datTmb$toa[,ncol(inp$datTmb$toa)])) > 0)
 
 	stopifnot(dim(inp$datTmb$H)[2] == 3)
 
 	# if z_vec != NULL 
-	if(inp$datTmb$how_3d != 'none'){
+	if(!inp$datTmb$how_3d %in% c('none', 'est')){
 		stopifnot(length(inp$datTmb$z_vec) == inp$datTmb$np)
 	}
 

R/getBbox.R

@@ -8,7 +8,7 @@
 #' @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){
+getBbox <- function(hydros, buffer=100, eps=1E-3, pen=1E6){
 	x_min <- hydros[which.min(hydros$hx), hx] - buffer
 	x_max <- hydros[which.max(hydros$hx), hx] + buffer
 	y_min <- hydros[which.min(hydros$hy), hy] - buffer

R/getInp.R

@@ -23,11 +23,13 @@ getInp <- function(hydros, toa, E_dist, n_ss, pingType, sdInits=1, rbi_min=0, rb
 	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)
+	bounds 		<- getBounds(datTmb)
 	return(list(
 		datTmb = datTmb,
 		params= params,
 		inits = inits,
-		inp_params = inp_params
+		inp_params = inp_params,
+		bounds = bounds
 		)
 	)
 }

R/plotBbox.R

@@ -7,10 +7,10 @@ plotBbox <- function(hydros, bbox){
 	Var1 <- Var2 <- NULL
 	softplus <- function(x, eps){return (0.5*(x+sqrt(x*x+eps*eps)))	}
 
-	x_space_min <- floor(hydros[which.min(hydros$hx), hx] - 20)
-	x_space_max <- ceiling(hydros[which.max(hydros$hx), hx] + 20)
-	y_space_min <- floor(hydros[which.min(hydros$hy), hy] - 20)
-	y_space_max <- ceiling(hydros[which.max(hydros$hy), hy] + 20)
+	x_space_min <- floor(hydros[which.min(hydros$hx), hx] - 100)
+	x_space_max <- ceiling(hydros[which.max(hydros$hx), hx] + 100)
+	y_space_min <- floor(hydros[which.min(hydros$hy), hy] - 100)
+	y_space_max <- ceiling(hydros[which.max(hydros$hy), hy] + 100)
 
 	xs <- x_space_min:x_space_max
 	ys <- y_space_min:y_space_max
@@ -33,7 +33,10 @@ plotBbox <- function(hydros, bbox){
 	pen_long[, y := ys[Var2]]
 	pen_long[, pen := value]
 
-	p <- ggplot2::ggplot(pen_long) + geom_tile(aes(x, y, fill=(pen))) + coord_fixed(ratio=1) + viridis::scale_fill_viridis() + geom_point(data=hydros, aes(x=hx, y=hy), col="red", size=2)
+	p <- ggplot2::ggplot(pen_long) + geom_tile(aes(x, y, fill=log(pen))) + 
+		coord_fixed(ratio=1, xlim=range(c(xs, x_min, x_max)), ylim=range(c(ys, y_min, y_max))) + 
+		viridis::scale_fill_viridis() + 
+		geom_point(data=hydros, aes(x=hx, y=hy), col="red", size=2)
 	print(p)
 }
 

R/prepTmb.R

@@ -62,8 +62,8 @@ getDatTmb <- function(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max, ss_
 		np = ncol(toa),
 		Edist = Edist,
 		toa = toa,
-		bi_epsilon = 1e-6,
-		bi_penalty = 1*1e9,
+		bi_epsilon = 1E-6,
+		bi_penalty = 1E9,
 		rbi_min = rbi_min,
 		rbi_max = rbi_max,
 		pingType = pingType,
@@ -104,10 +104,13 @@ getParams <- function(datTmb){
 		# , log_t_part = 0				#Mixture ratio between Gaussian and t
 	)
 
+	# # # If estimating 3D	
 	if(datTmb$how_3d == "est"){
 		out$Z <- stats::runif(ncol(datTmb$toa), -10, 0)
+		out$logD_z <- 0				#diffusivity of transmitter vertical movement (D_z in ms)
 	}
 
+
 	# # # ss related
 	if(datTmb$ss_data_what == 'est'){
 		out$logD_v <- 0				#diffusivity of speed of sound (D_v in ms)
@@ -126,16 +129,15 @@ getParams <- function(datTmb){
 	}
 
 	# # # Ping type related
-	if(datTmb$pingType == 'sbi'){
+	if(datTmb$pingType %in% c('sbi', 'sbi_double', 'rbi')){
 		out$logSigma_bi <- 0			#sigma  burst interval (sigma_bi in ms)
 	}
 
 	if(datTmb$pingType == 'pbi'){
 		out$logSigma_bi <- 0			#sigma  burst interval (sigma_bi in ms)
 		out$tag_drift <- stats::rnorm(datTmb$np, 0, 1e-2)
 	}
-
-
+
 
 	return(out)
 }
@@ -184,6 +186,11 @@ getInits <- function(datTmb, sdInits=1) {
 	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)
@@ -201,16 +208,69 @@ getInits <- function(datTmb, sdInits=1) {
 	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_logD_v)#, init_logSigma_toa, init_logScale, init_log_t_part)
+		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)
 }
 
+#' 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)
+}
+
 #' Get parameters for this specific data set
 #'
 #' Compile a list of relevant parameters (e.g. T0) to use later on

R/runYaps.R

@@ -8,8 +8,7 @@
 #' @param opt_controls List of controls passed to optimization function. For instances, tolerances such as `x.tol=1E-8`. \cr 
 #'   If `opt_fun = 'nloptr'`, `opt_controls` must be a list formatted appropriately. For instance: \cr
 #'  `opt_controls <- list( algorithm="NLOPT_LD_AUGLAG", xtol_abs=1e-12, maxeval=2E+4, print_level = 1, local_opts= list(algorithm="NLOPT_LD_AUGLAG_EQ", xtol_rel=1e-4) )`. \cr
-#'  See `?nloptr` and the NLopt site https://nlopt.readthedocs.io/en/latest/ for more info. Some algorithms in `nloptr` require bounded parameters - see `bounds`.
-#' @param bounds List of two vectors specifying lower and upper bounds of fixed parameters. Length of each vector must be equal to number of fixed parameters. For instance, `bounds = list(lb = c(-3, -1, -2), ub = c(2,0,1) )`.
+#'  See `?nloptr` and the NLopt site https://nlopt.readthedocs.io/en/latest/ for more info. Some algorithms in `nloptr` require bounded parameters - this is not currently implemented.
 #' @param tmb_smartsearch Logical whether to use the TMB smartsearch in the inner optimizer (see `?TMB::MakeADFun` for info). Default and original implementation is TRUE. However, there seems to be an issue with recent versions of `Matrix` that requires `tmb_smartsearch=FALSE`. 
 #'
 #' @export
@@ -27,7 +26,7 @@
 #'  }
 #'
 #' @example man/examples/example-yaps_ssu1.R
-runYaps <- function(inp, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE, opt_fun='nlminb', opt_controls=list(), bounds=list(), tmb_smartsearch=TRUE){
+runYaps <- function(inp, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE, opt_fun='nlminb', opt_controls=list(), tmb_smartsearch=TRUE){
 
 	# making sure inp is correct...
 	checkInp(inp)
@@ -51,39 +50,55 @@ runYaps <- function(inp, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE, o
 			parameters = inp$params,
 			random = random,
 			DLL = "yaps",
-			inner.control = list(maxit = maxIter), 
+			inner.control = list(maxit = maxIter),
 			silent=silent
 		)
+
+	# Attempt to robustify the inner optim problem.
+	# Refuse to optimize if gradient is too steep. Default is 1E60
+	# TMB::newtonOption(obj, tol10=1)
+
+	# # # EXPERIMENTAL
+	if(inp$datTmb$pingType == 'rbi'){
+		TMB::newtonOption(obj, mgcmax=1E6)
+	}
+
+	if(!silent){
+		obj$env$tracepar = TRUE
+		obj$env$tracemgc = TRUE
+	}
+
 	if(!tmb_smartsearch){
 		obj$fn(obj$par) 
 		TMB::newtonOption(obj, smartsearch=FALSE)
 	}
 
-	if(	!is.null(opt_controls[['use_bounds']])){
-		lower <- opt_controls[['lower']]
-		upper <- opt_controls[['upper']]
-		opt_controls <- list()
-	} else {
-		lower <- -Inf
-		upper <- Inf
-	}
+	# if(	!is.null(opt_controls[['use_bounds']])){
+		# lower <- opt_controls[['lower']]
+		# upper <- opt_controls[['upper']]
+		# opt_controls <- list()
+	# } else {
+		# lower <- -Inf
+		# upper <- Inf
+	# }
 
 	if(opt_fun == 'nloptr'){
 		opts <- opt_controls
-		if(length(bounds) == 0){
-			opt <- nloptr::nloptr(inp$inits,obj$fn,obj$gr,opts=opts,...=NA)
-		} else {
-			opt <- nloptr::nloptr(inp$inits,obj$fn,obj$gr,opts=opts,...=NA, lb=bounds$lb, ub=bounds$ub)
-		}
+		# if(length(bounds) == 0){
+			# opt <- nloptr::nloptr(inp$inits,obj$fn,obj$gr,opts=opts,...=NA)
+		# } else {
+			# opt <- nloptr::nloptr(inp$inits,obj$fn,obj$gr,opts=opts,...=NA, lb=inp$bounds[,1], ub=inp$bounds[,2])
+		# }
+		opt <- nloptr::nloptr(inp$inits,obj$fn,obj$gr,opts=opts,...=NA)#, lb=inp$bounds[,1], ub=inp$bounds[,2])
 
 	} else if(opt_fun == 'nlminb'){
 		control_list <- opt_controls
 		if(!silent){
 				# opt <- stats::nlminb(inp$inits,obj$fn,obj$gr, control = control_list)
 				# opt <- stats::nlminb(inp$inits,obj$fn,obj$gr, control = control_list, lower=c(-50,-15, -100, -50, -20), upper= c(2, 2, 100, 2, -2))
-				opt <- stats::nlminb(inp$inits,obj$fn,obj$gr, control = control_list, lower=lower, upper=upper)
+				opt <- stats::nlminb(inp$inits,obj$fn,obj$gr, control = control_list, lower=inp$bounds[,1], upper=inp$bounds[,2])
 		} else {
-			suppressWarnings(opt <- stats::nlminb(inp$inits,obj$fn,obj$gr, control = control_list, lower=lower, upper=upper))
+			suppressWarnings(opt <- stats::nlminb(inp$inits,obj$fn,obj$gr, control = control_list, lower=inp$bounds[,1], upper=inp$bounds[,2]))
 		}
 	}
 

R/syncPlotters.R

@@ -14,7 +14,7 @@ plotSyncModelResids <- function(sync_model, by='overall'){
 		p <- p + xlab("Residual (m)")
 	} else if(by == 'sync_tag'){
 		p <- ggplot2::ggplot(data=eps_long) 
-		p <- p + geom_violin(aes(x=factor(hydro_idx), y=E_m), scale="count", draw_quantiles=c(.25,.5, .75)) 
+		p <- p + geom_violin(aes(x=factor(hydro_idx), y=E_m), scale="count") 
 		p <- p + geom_vline(xintercept=seq(0, length(unique(eps_long$hydro_idx)), by=5), lty=3, col="red")
 		p <- p + geom_hline(yintercept=0, col="red") 
 		p <- p + facet_wrap(~sync_tag_idx)
@@ -23,7 +23,7 @@ plotSyncModelResids <- function(sync_model, by='overall'){
 		p <- p + xlab("hydro_idx")
 		p <- p + ggplot2::scale_x_discrete(breaks = pretty(unique(eps_long$hydro_idx)))
 	} else if(by == 'hydro'){
-		p <- ggplot2::ggplot(data=eps_long) + geom_violin(aes(x=factor(sync_tag_idx), y=E_m), scale="count", draw_quantiles=c(.25,.5, .75)) 
+		p <- ggplot2::ggplot(data=eps_long) + geom_violin(aes(x=factor(sync_tag_idx), y=E_m), scale="count") 
 		p <- p + geom_hline(yintercept=0, col="red") + facet_wrap(~hydro_idx) 
 		# p <- p + labs(title="Sync model residuals - facet by hydro_idx") + xlab("sync_tag_idx")
 		p <- p + xlab("sync_tag_idx")
@@ -109,7 +109,7 @@ plotSyncModelCheck <- function(sync_model, by=""){
 		p <- p + xlab("Sync period")
 	} else if(by == "sync_bin_hydro"){
 		plot_dat <- sync_check_dat[, .(med_delta=median((delta))), by=c('sync_tag_idx', 'focal_hydro_idx', 'hydro_idx', 'offset_idx')]
-		p <- ggplot2::ggplot(data=plot_dat) + geom_violin(aes(x=factor(offset_idx), y=(med_delta)), scale="count", draw_quantiles=c(.25,.5, .75))  + facet_wrap(~focal_hydro_idx)
+		p <- ggplot2::ggplot(data=plot_dat) + geom_violin(aes(x=factor(offset_idx), y=(med_delta)), scale="count")  + facet_wrap(~focal_hydro_idx)
 		p <- p + xlab("Sync period")
 	} else if(by == "sync_bin_sync_smooth"){
 		plot_dat <- sync_check_dat[, .(med_delta=median((delta))), by=c('sync_tag_idx', 'focal_hydro_idx', 'hydro_idx', 'offset_idx')]

R/testYaps.R

@@ -11,15 +11,17 @@
 #' @return If `return_yaps == TRUE`, the fitted `yaps` object. See `?runYaps` for further info.
 
 #' @examples
+#' \donttest{
 #' #' # To test  basic functionality of yaps using simulated data
 #' testYaps()
 #' # # # Three pingTypes are availabe: 
 #' # # #	fixed burst interval (testYaps(pingType='sbi')), 
 #' # # # 	random burst interval with UNKNOWN burst interval sequence('testYaps(pingType='rbi')), 
 #' # # # 	random burst interval with KNOWN burst interval sequence (testYaps(pingType='pbi'))
-testYaps <- function(silent=TRUE, pingType='sbi', est_ss=TRUE, opt_fun='nlminb', opt_controls=list(), bounds=list(), return_yaps=FALSE, tmb_smartsearch=TRUE){
+#' }
+testYaps <- function(silent=TRUE, pingType='sbi', est_ss=TRUE, opt_fun='nlminb', opt_controls=list(), return_yaps=FALSE, tmb_smartsearch=TRUE){
 	set.seed(42)
-	trueTrack <- simTrueTrack(model='crw', n = 2000, deltaTime=1, shape=1, scale=0.5, addDielPattern=TRUE, ss='rw')
+	trueTrack <- simTrueTrack(model='crw', n = 2500, deltaTime=1, shape=1, scale=0.5, addDielPattern=FALSE, ss='rw')
 	if(pingType == 'sbi'){
 		sbi_mean <- 20; sbi_sd <- 1e-3;
 		rbi_min <- sbi_mean; 
@@ -49,7 +51,7 @@ testYaps <- function(silent=TRUE, pingType='sbi', est_ss=TRUE, opt_fun='nlminb',
 	}
 	inp <- getInp(hydros, toa, E_dist="Mixture", n_ss=5, pingType=pingType, sdInits=0, ss_data_what=ss_data_what, ss_data=ss_data, rbi_min=rbi_min, rbi_max=rbi_max, biTable=biTable)
 	maxIter <- 500
-	yaps <- runYaps(inp, maxIter=maxIter, getPlsd=TRUE, getRep=TRUE, silent=silent, opt_fun=opt_fun, opt_controls, bounds, tmb_smartsearch)
+	yaps <- runYaps(inp, maxIter=maxIter, getPlsd=TRUE, getRep=TRUE, silent=silent, opt_fun=opt_fun, opt_controls, tmb_smartsearch)
 	pl <- yaps$pl
 	yaps_out <- data.table::data.table(X=pl$X + inp$inp_params$Hx0, Y=pl$Y + inp$inp_params$Hy0)
 	plsd <- yaps$plsd