Skip to content

Commit

Permalink
adapting getInp() - round 1
Browse files Browse the repository at this point in the history
  • Loading branch information
@baktoft
baktoft committed Oct 20, 2023
1 parent a88bb1e commit ac45026
Show file tree
Hide file tree
Showing 10 changed files with 161 additions and 91 deletions.
8 changes: 4 additions & 4 deletions R/getBbox.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,10 @@
#' @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=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
y_max <- hydros[which.max(hydros$hy), hy] + buffer
x_min <- hydros[which.min(hydros$h_x), h_x] - buffer
x_max <- hydros[which.max(hydros$h_x), h_x] + buffer
y_min <- hydros[which.min(hydros$h_y), h_y] - buffer
y_max <- hydros[which.max(hydros$h_y), h_y] + buffer
bbox <- c(x_min, x_max, y_min, y_max, eps, pen)
return(bbox)
}
24 changes: 12 additions & 12 deletions R/getBounds.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,18 +4,18 @@
#' @param datTmb Object obtained using getDatTmb()
#' @return Matrix of bounds restricting the optimizer when running runYaps().
#' @noRd
getBounds <- function(datTmb) {
getBounds <- function(dat_tmb) {
lu_logD_xy <- c(-50, 2)
lu_logD_z <- c(-50, 2)

lu_logSigma_toa <- c(-12, -2)
if(datTmb$Edist[2] == 1){ # mixture
if(dat_tmb$E_dist_vec[2] == 1){ # mixture
lu_logScale <- c(-30, 10)
} else if (datTmb$Edist[3] == 1) { # t
} else if (dat_tmb$E_dist_vec[3] == 1) { # t
lu_logScale <- c(-10,2)
}
lu_log_t_part <- c(-100, 100)
if(datTmb$pingType == 'rbi'){
if(dat_tmb$ping_type == 'rbi'){
lu_logSigma_bi <- c(-20, 20)
} else {
lu_logSigma_bi <- c(-20, -2)
Expand All @@ -25,27 +25,27 @@ getBounds <- function(datTmb) {
bounds <- c()
bounds <- rbind(bounds, lu_logD_xy)

if(datTmb$how_3d == 'est'){
if(dat_tmb$how_3d == 'est'){
bounds <- rbind(bounds, lu_logD_z)
}

if(datTmb$ss_data_what == 'est'){
if(dat_tmb$ss_data == 'none'){
bounds <- rbind(bounds, lu_logD_v)
}

if(datTmb$Edist[1] == 1){
if(dat_tmb$E_dist_vec[1] == 1){
bounds <- rbind(bounds, lu_logSigma_toa)
} else if(datTmb$Edist[2] == 1){
} else if(dat_tmb$E_dist_vec[2] == 1){
bounds <- rbind(bounds, lu_logSigma_toa, lu_logScale, lu_log_t_part)
} else if(datTmb$Edist[3] == 1){
} else if(dat_tmb$E_dist_vec[3] == 1){
bounds <- rbind(bounds, lu_logScale)
}

if(datTmb$pingType == 'sbi'){
if(dat_tmb$ping_type == 'sbi'){
bounds <- rbind(bounds, lu_logSigma_bi)
} else if (datTmb$pingType == 'rbi'){
} else if (dat_tmb$ping_type == 'rbi'){
bounds <- rbind(bounds, lu_logSigma_bi)
} else if (datTmb$pingType == 'pbi'){
} else if (dat_tmb$ping_type == 'pbi'){
bounds <- rbind(bounds, lu_logSigma_bi)
}

Expand Down
65 changes: 33 additions & 32 deletions R/getDatTmb.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,49 +6,51 @@
#'
#' @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){
# getDatTmb <- function(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max, ss_data_what, ss_data, biTable, inp_params, z_vec, bbox){
getDatTmb <- function(hydros, toa, ss_data, yaps_params, inp_params){

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
rbi_min <- yaps_params$bi[1] - yaps_params$bi[1] * 0.05
rbi_max <- yaps_params$bi[2] + yaps_params$bi[2] * 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

if(yaps_params$n_ss > 1){
ss_idx <- cut(1:nrow(toa), yaps_params$n_ssn_ss, labels=FALSE) - 1 #-1 because zero-indexing in TMB
} else {
ss_idx <- rep(0, ncol(toa))
ss_idx <- rep(0, nrow(toa))
}
approxBI <- mean(diff(colMeans(toa, na.rm=TRUE), na.rm=TRUE), na.rm=TRUE)

approx_bi <- mean(diff(rowMeans(toa, na.rm=TRUE), na.rm=TRUE), na.rm=TRUE)

if(ss_data_what == 'data') { stopifnot(length(ss_data) == ncol(toa))}
if(ss_data != 'none' & length(ss_data) != nrow(toa) ){
cat("ERROR: Seems like ss_data is provided, but length(ss_data) != nrow(toa) \n")
cat("...getInp() found ss_data != 'none' and ",length(ss_data)," != ",nrow(toa)," \n")
stopSilent()
}

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}
E_dist_vec <- rep(0,3)
if(yaps_params$E_dist == "Gaus") {E_dist_vec[1] <- 1}
if(yaps_params$E_dist == "Mixture") {E_dist_vec[2] <- 1}
if(yaps_params$E_dist == "t") {E_dist_vec[3] <- 1}

if(is.null(z_vec)){
if(z_data[1] == 'none'){
how_3d <- 'none'
z_vec <- c(1)
} else if(z_vec[1] == "est") {
} else if(z_data[1] == 'est'){
how_3d <- 'est'
z_vec <- c(1)
} else {
how_3d <- 'data'
}

if(is.null(bbox)){
if(yaps_params$bbox != TRUE){
bbox <- NA
} else {
bbox <- getBbox(hydros)
bbox[1] <- bbox[1] - inp_params$Hx0
bbox[2] <- bbox[2] - inp_params$Hx0
bbox[3] <- bbox[3] - inp_params$Hy0
Expand All @@ -57,27 +59,26 @@ getDatTmb <- function(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max, ss_

datTmb <- list(
model = "yaps_track",
H = matrix(c(hydros$h_x-Hx0, hydros$h_y-Hy0, hydros$h_z), ncol=3),
H = matrix(c(hydros$h_x, hydros$h_y, hydros$h_z), ncol=3),
nh = nrow(hydros),
np = ncol(toa),
Edist = Edist,
np = nrow(toa),
E_dist_vec = E_dist_vec,
toa = toa,
bi_epsilon = 1E-6,
bi_penalty = 1E9,
rbi_min = rbi_min,
rbi_max = rbi_max,
pingType = pingType,
n_ss = n_ss,
ping_type = yaps_params$ping_type,
n_ss = yaps_params$n_ss,
ss_idx = ss_idx,
ss_data_what = ss_data_what,
ss_data = ss_data,
approxBI = approxBI,
biTable = c(1),
approx_bi = approx_bi,
# biTable = c(1), # NOT IMPLEMENTED YET
how_3d = how_3d,
z_vec = z_vec,
z_data = z_data,
bbox = bbox
)
if(pingType == 'pbi') {datTmb$biTable = biTable}
# if(pingType == 'pbi') {datTmb$biTable = biTable} # NOT IMPLEMENTED YET

return(datTmb)
}
28 changes: 14 additions & 14 deletions R/getInits.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,18 +6,18 @@
#' @inheritParams getInp
#' @return Vector of initial values to use in TMB
#' @noRd
getInits <- function(datTmb, sdInits=1) {
getInits <- function(dat_tmb, yaps_params) {
init_logD_xy <- -1

if(datTmb$pingType == 'sbi') {
if(dat_tmb$ping_type == 'sbi') {
init_logSigma_bi <- -6
} else if(datTmb$pingType == 'rbi'){
if(datTmb$rbi_max >= 10){
} else if(dat_tmb$ping_type == 'rbi'){
if(dat_tmb$rbi_max >= 10){
init_logSigma_bi <- 4
} else {
init_logSigma_bi <- -3
}
} else if(datTmb$pingType == 'pbi'){
} else if(dat_tmb$ping_type == 'pbi'){
init_logSigma_bi <- -5
}

Expand All @@ -28,32 +28,32 @@ getInits <- function(datTmb, sdInits=1) {

inits <- c(init_logD_xy)

if(datTmb$how_3d == 'est'){
if(dat_tmb$how_3d == 'est'){
init_logD_z <- 0
inits <- c(inits, init_logD_z)
}

if(datTmb$ss_data_what == 'est'){
if(dat_tmb$ss_data == 'none'){
inits <- c(inits, init_logD_v)
}


if(datTmb$Edist[1] == 1){
if(dat_tmb$E_dist_vec[1] == 1){
inits <- c(inits, init_logSigma_toa)
} else if(datTmb$Edist[2] == 1){
} else if(dat_tmb$E_dist_vec[2] == 1){
inits <- c(inits, init_logSigma_toa, init_logScale, init_log_t_part)
} else if(datTmb$Edist[3] == 1){
} else if(dat_tmb$E_dist_vec[3] == 1){
inits <- c(inits, init_logScale)
}

if(datTmb$pingType == 'sbi'){
if(dat_tmb$ping_type == 'sbi'){
inits <- c(inits, init_logSigma_bi)#, init_logD_v)#, init_logSigma_toa, init_logScale, init_log_t_part)
} else if (datTmb$pingType == 'rbi'){
} else if (dat_tmb$ping_type == 'rbi'){
inits <- c(inits, init_logSigma_bi)#, init_logD_v)#, init_logSigma_toa, init_logScale, init_log_t_part)
} else if (datTmb$pingType == 'pbi'){
} else if (dat_tmb$ping_type == '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)
inits <- stats::rnorm(length(inits), mean=inits, sd=yaps_params$sd_inits)
return(inits)
}
28 changes: 21 additions & 7 deletions R/getInp.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,16 +17,30 @@
#' @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'))
# 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){
getInp <- function(hydros, toa, ss_data='none', z_data='none', yaps_params){

yaps_params <- prepYapsParams(yaps_params, ss_data, z_data)
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)
inits <- getInits(datTmb, sdInits)
bounds <- getBounds(datTmb)
if(ncol(toa) != nrow(hydros)){
cat("ERROR: ncol(toa) != nrow(hydros) \n")
cat("...getInp found ",ncol(toa)," != ",nrow(hydros),"\n")
stopSilent()
}

if(z_data[1] != 'none' & length(z_data) != nrow(toa)){
cat("ERROR: nrow(toa) != length(z_data) \n")
cat("...getInp found ",nrow(toa)," != ",length(z_data),"\n")
stopSilent()
}


inp_params <- getInpParams(hydros, toa)
dat_tmb <- getDatTmb(hydros, toa, ss_data, yaps_params, inp_params)
params <- getParams(dat_tmb)
inits <- getInits(dat_tmb, yaps_params)
bounds <- getBounds(dat_tmb)
return(list(
datTmb = datTmb,
params= params,
Expand Down
2 changes: 1 addition & 1 deletion R/getInpParams.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
#' Compile a list of relevant parameters (e.g. T0) to use later on
#' @inheritParams getInp
#' @noRd
getInpParams <- function(hydros, toa, pingType){
getInpParams <- function(hydros, toa){
T0 <- min(toa, na.rm=TRUE)

Hx0 <- hydros[1,h_x]
Expand Down
30 changes: 15 additions & 15 deletions R/getParams.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,14 @@
#' @param datTmb Object obtained using getDatTmb()
#' @return List of params for use in TMB
#' @noRd
getParams <- function(datTmb){
params_XY <- yaps:::getParamsXYFromCOA(datTmb)
getParams <- function(dat_tmb){
params_XY <- getParamsXYFromCOA(dat_tmb)
out <- list(
X = params_XY$X + stats::rnorm(ncol(datTmb$toa), sd=10)
, Y = params_XY$Y + stats::rnorm(ncol(datTmb$toa), sd=10)
X = params_XY$X + stats::rnorm(nrow(dat_tmb$toa), sd=10)
, Y = params_XY$Y + stats::rnorm(nrow(dat_tmb$toa), sd=10)
# X = 0 + stats::rnorm(ncol(datTmb$toa), sd=10)
# , Y = 0 + stats::rnorm(ncol(datTmb$toa), sd=10)
, top = zoo::na.approx(apply(datTmb$toa, 2, function(k) {stats::median(k, na.rm=TRUE)}), rule=2) #time of ping
, TOP = zoo::na.approx(apply(dat_tmb$toa, 1, function(k) {stats::median(k, na.rm=TRUE)}), rule=2) #time of ping
# , SS=stats::rnorm(datTmb$n_ss, 1450, 5) #speed of sound
, logD_xy = 0 #diffusivity of transmitter movement (D_xy in ms)
# , logD_v = 0 #diffusivity of speed of sound (D_v in ms)
Expand All @@ -21,37 +21,37 @@ getParams <- function(datTmb){
)

# # # If estimating 3D
if(datTmb$how_3d == "est"){
out$Z <- stats::runif(ncol(datTmb$toa), -10, 0)
if(dat_tmb$how_3d == "est"){
out$Z <- stats::runif(nrow(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'){
if(dat_tmb$ss_data == 'none'){
out$logD_v <- 0 #diffusivity of speed of sound (D_v in ms)
out$SS <- stats::rnorm(datTmb$n_ss, 1450, 5) #speed of sound
out$SS <- stats::rnorm(dat_tmb$n_ss, 1450, 5) #speed of sound
}

# # # Edist related
if(datTmb$Edist[1] == 1){
if(dat_tmb$E_dist_vec[1] == 1){
out$logSigma_toa = 0 #sigma for Gaussian
} else if(datTmb$Edist[2] == 1){
} else if(dat_tmb$E_dist_vec[2] == 1){
out$logSigma_toa = 0 #sigma for Gaussian
out$logScale = 0 #scale parameter for t-distribution
out$log_t_part = 0 #Mixture ratio between Gaussian and t
} else if(datTmb$Edist[3] == 1){
} else if(dat_tmb$E_dist_vec[3] == 1){
out$logScale = 0 #scale parameter for t-distribution
}

# # # Ping type related
if(datTmb$pingType %in% c('sbi', 'sbi_double', 'rbi')){
if(dat_tmb$ping_type %in% c('sbi', 'sbi_double', 'rbi')){
out$logSigma_bi <- 0 #sigma burst interval (sigma_bi in ms)
}

if(datTmb$pingType == 'pbi'){
if(dat_tmb$ping_type == 'pbi'){
out$logSigma_bi <- 0 #sigma burst interval (sigma_bi in ms)
out$tag_drift <- stats::rnorm(datTmb$np, 0, 1e-2)
out$tag_drift <- stats::rnorm(dat_tmb$np, 0, 1e-2)
}


Expand Down
10 changes: 5 additions & 5 deletions R/getParamsXYFromCOA.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,15 +3,15 @@
#' Attempts to give meaningful initial values for X and Y based on which hydros detected each ping
#' @inheritParams getInp
#' @noRd
getParamsXYFromCOA <- function(datTmb){
toa <- datTmb$toa
hydros <- datTmb$H
getParamsXYFromCOA <- function(dat_tmb){
toa <- dat_tmb$toa
hydros <- dat_tmb$H

toa_detect <- toa
toa_detect[!is.na(toa_detect)] <- 1

X <- zoo::na.approx(colMeans((toa_detect) * hydros[,1], na.rm=TRUE))
Y <- zoo::na.approx(colMeans((toa_detect) * hydros[,2], na.rm=TRUE))
X <- zoo::na.approx(rowMeans((toa_detect) * hydros[,1], na.rm=TRUE))
Y <- zoo::na.approx(rowMeans((toa_detect) * hydros[,2], na.rm=TRUE))

return(list(X=X, Y=Y))

Expand Down
Loading

0 comments on commit ac45026

Please sign in to comment.