diff --git a/DESCRIPTION b/DESCRIPTION index 39e4cf8..9d1093a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: yaps Title: Track Estimation using YAPS (Yet Another Positioning Solver) -Version: 1.2.5.9100 +Version: 2.0.0.9000 Authors@R: c( person("Henrik", "Baktoft", email = "hba@aqua.dtu.dk", 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")), diff --git a/NAMESPACE b/NAMESPACE index c80130f..9d7045e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,15 +1,24 @@ # Generated by roxygen2: do not edit by hand export(alignBurstSeq) +export(applyLinCor) export(applySync) export(checkInp) export(checkInpSync) +export(checkInpSyncData) +export(downsampleToaList_random) +export(downsampleToaList_selective) export(fineTuneSyncModel) export(getBbox) +export(getDownsampledToaList) export(getInp) export(getInpSync) +export(getOffsetVals) +export(getParamsTmbSync) +export(getRandomTmbSync) export(getSyncCoverage) export(getSyncModel) +export(getSyncToa) export(getToaYaps) export(plotBbox) export(plotSyncModelCheck) @@ -25,6 +34,7 @@ export(simToa) export(simTrueTrack) export(tempToSs) export(testYaps) +export(yapsifyHydros) importFrom(Rcpp,sourceCpp) importFrom(data.table,"%between%") importFrom(data.table,"%like%") diff --git a/NEWS.md b/NEWS.md index 766f666..3f1dfe1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,14 @@ +# yaps v.2.0.0.9000 + +* Work in progress - still in early phase +* Complete rework of all/most/many functions +* Will not be backward compatible +* Will break current code +* Improve functions to obtain sync models +* Attempt to improve workflow +* Remove dependence of deprecated package splusTimeSeries + + # yaps v.1.2.5.9000 ## New stuff diff --git a/R/getDatTmbSync.R b/R/getDatTmbSync.R deleted file mode 100644 index 18bfe6b..0000000 --- a/R/getDatTmbSync.R +++ /dev/null @@ -1,39 +0,0 @@ -#' Internal function to get dat for TMB sync -#' @inheritParams getInpSync -#' @noRd -getDatTmbSync <- function(sync_type, sync_dat, keep_rate, time_keeper_idx, inp_toa_list, fixed_hydros_vec, offset_vals, ss_vals, inp_H_info, T0, ss_data_what, ss_data_vec){ - H <- as.matrix(inp_H_info$inp_H) - dimnames(H) <- NULL - - toa <- inp_toa_list$toa - T0 - toa_offset <- inp_toa_list$toa - offset_vals$offset_levels[offset_vals$offset_idx] - - if(sync_type == 'delta'){ - toa_delta <- getToaDelta(toa, inp_toa_list, keep_rate, time_keeper_idx, offset_vals, ss_vals) - } else { - toa_delta <- matrix(0) - } - - dat_tmb_sync <- list( - model = 'yaps_sync', - sync_type = sync_type, - H=H, - toa = toa, - toa_offset=toa_offset, - toa_delta = toa_delta, - sync_tag_idx_vec = inp_toa_list$sync_tag_idx_vec, - np = nrow(inp_toa_list$toa), - nh = ncol(inp_toa_list$toa), - ndelta = nrow(toa_delta), - tk = time_keeper_idx, - fixed_hydros_vec = fixed_hydros_vec, - offset_idx = offset_vals$offset_idx, - n_offset_idx = offset_vals$n_offset_idx, - ss_idx = ss_vals$ss_idx, - n_ss_idx = ss_vals$n_ss_idx, - ss_data_what = ss_data_what, - ss_data_vec = ss_data_vec - ) - return(dat_tmb_sync) -} - diff --git a/R/getSyncModel.R b/R/getSyncModel.R index f6e2c26..b8789ce 100644 --- a/R/getSyncModel.R +++ b/R/getSyncModel.R @@ -9,8 +9,7 @@ #' @export #' @return List containing relevant data constituting the `sync_model` ready for use in `fineTuneSyncModel()` if needed or in `applySync()` #' @example man/examples/example-yaps_ssu1.R -getSyncModel <- function(inp_sync, silent=TRUE, fine_tune=FALSE, max_iter=100, tmb_smartsearch=TRUE){ - inp_sync$inp_params$tmb_smartsearch <- tmb_smartsearch +getSyncModel <- function(inp_sync, silent=TRUE, max_iter=500){ inp_sync$inp_params$max_iter <- max_iter dat_tmb <- inp_sync$dat_tmb_sync @@ -21,7 +20,6 @@ getSyncModel <- function(inp_sync, silent=TRUE, fine_tune=FALSE, max_iter=100, t cat(paste0(Sys.time(), " \n")) cat(". Running optimization of the sync model. Please be patient - this can take a long time. \n") - if(fine_tune){cat(".... fine tuning is enabled, but is getting deprecated in future version. Consider to use the function fineTuneSyncModel() instead. See ?fineTuneSyncModel for info. \n")} opt <- c() pl <- c() @@ -44,11 +42,6 @@ getSyncModel <- function(inp_sync, silent=TRUE, fine_tune=FALSE, max_iter=100, t obj$env$tracemgc = TRUE } - - if(!tmb_smartsearch){ - TMB::newtonOption(obj, smartsearch=FALSE) - } - lower <- c(-10) upper <- c(-2) @@ -69,13 +62,9 @@ getSyncModel <- function(inp_sync, silent=TRUE, fine_tune=FALSE, max_iter=100, t crazy_outliers <- which(abs(report$eps)*1450 > 10000) fine_outliers <- which(abs(report$eps)*1450 > 1000) if(length(crazy_outliers > 0)){ - cat(".... some extreme outliers potentially affecting the model where identified \n Consider to run fineTuneSyncModel(sync_model, eps_threshold=10000). See ?fineTuneSyncModel for more info. \n") + cat(".... some extreme outliers potentially affecting the model where identified \n Consider running fineTuneSyncModel(sync_model, eps_threshold=10000). See ?fineTuneSyncModel for more info. \n") # dat_tmb$toa_offset[crazy_outliers] <- NA - } else if(fine_tune){ - cat(".... fine tuning is enabled, but is deprecated. Use the function fineTuneSyncModel() instead. See ?fineTuneSyncModel for info. \n") - # dat_tmb$toa_offset[fine_outliers] <- NA - } - + } jointrep <- try(TMB::sdreport(obj, getJointPrecision=TRUE), silent=silent) param_names <- rownames(summary(jointrep)) @@ -83,8 +72,8 @@ getSyncModel <- function(inp_sync, silent=TRUE, fine_tune=FALSE, max_iter=100, t summ <- data.frame(param=param_names, sd=sds) plsd <- split(summ[,2], f=summ$param) - pl$TRUE_H[,1] <- pl$TRUE_H[,1] + inp_params$Hx0 - pl$TRUE_H[,2] <- pl$TRUE_H[,2] + inp_params$Hy0 + pl$TRUE_H[,1] <- pl$TRUE_H[,1] + attr(inp_params$hydros, 'Hx0') + pl$TRUE_H[,2] <- pl$TRUE_H[,2] + attr(inp_params$hydros, 'Hy0') eps_long <- getEpsLong(report, pl, inp_sync) offset_nas <- which(pl$OFFSET == 0) diff --git a/R/syncGetters.R b/R/syncGetters.R index ab7e961..62577c1 100644 --- a/R/syncGetters.R +++ b/R/syncGetters.R @@ -9,141 +9,15 @@ getSsDataVec <- function(inp_toa_list, ss_data){ return(ss_data_vec) } -#' Internal function. Apply linear correction matrix to epofrac before sync -#' @inheritParams getInpSync -#' @noRd -applyLinCorCoeffsInpSync <- function(sync_dat, lin_corr_coeffs){ - h_idxs <- sync_dat$hydros$idx - - for(h in 1:length(h_idxs)){ - h_idx <- h_idxs[h] - # sync_dat$detections[] - sync_dat$detections[hydro_idx == h_idx, epofrac := epofrac - lin_corr_coeffs[h_idx, 1] - epofrac * lin_corr_coeffs[h_idx, 2]] - } - return(sync_dat) -} - - -#' Internal function. Get toa for sync from sync_dat -#' @inheritParams getInpSync -#' @noRd -getInpSyncToaList <- function(sync_dat, max_epo_diff, min_hydros, excl_self_detect, keep_rate, lin_corr_coeffs){ - toa_list_gross <- buildToaListGross(sync_dat, excl_self_detect) - toa_list_pruned <- pruneToaListGross(toa_list_gross, max_epo_diff, min_hydros) - - return(toa_list_pruned) -} - -#' Internal warpper function to do downsampling of the toa -#' @inheritParams getInpSync -#' @noRd -getDownsampledToaList <- function(inp_toa_list_all, offset_vals_all, keep_rate){ - if(keep_rate > 0 & keep_rate <= 1){ - toa_list_downsampled <- downsampleToaList_random(inp_toa_list_all, keep_rate) - } else if(keep_rate >= 10){ - toa_list_downsampled <- downsampleToaList_selective(inp_toa_list_all, offset_vals_all, keep_rate) - } - return(toa_list_downsampled) -} - -# Internal function to selectively downsample the toa-matrix for inp_sync -#' @param inp_toa_list_all Output from `getInpSyncToaList` -#' @inheritParams getInpSync -#' @noRd -downsampleToaList_selective <- function(inp_toa_list_all, offset_vals_all, keep_rate){ - toa <- inp_toa_list_all$toa - offset_idx <- offset_vals_all$offset_idx - - toa_long <- data.table::data.table(reshape2::melt(toa), rep(offset_idx, times=ncol(toa))) - colnames(toa_long) <- c('ping', 'h_idx','toa','offset_idx') - # nobs_per_offset <- toa_long[!is.na(toa), .N, by=c('h_idx' ,'offset_idx')] - - nobs_per_offset <- data.table::data.table(reshape2::melt(with(toa_long[!is.na(toa)], table(h_idx, offset_idx)), value.name="N")) - - - keep_pings <- c() - for(i in 1:length(unique(offset_idx))){ - toa_long_i <- toa_long[offset_idx == i] - keep_pings_i <- c() - # h_order <- order(nobs_per_offset[offset_idx == i, N]) - - h_dets <- nobs_per_offset[offset_idx == i, N, by=h_idx] - data.table::setorder(h_dets, N) - h_order <- h_dets$h_idx - - for(h in 1:length(h_order)){ - already_in_keeps <- nrow(toa_long_i[!is.na(toa) & ping %in% keep_pings_i & offset_idx == i & h_idx==h_order[h]]) - if(already_in_keeps >= keep_rate){ - # already enough with this hydro... - next - } else { - need <- keep_rate - already_in_keeps + 1 - pings_h <- toa_long_i[!is.na(toa) & offset_idx == i & h_idx==h_order[h], ping] - if(length(pings_h) < need){ - keep_pings_h <- pings_h - } else { - keep_pings_h <- sample(pings_h, size=need) - } - keep_pings_i <- c(keep_pings_i, keep_pings_h) - } - } - keep_pings <- c(keep_pings, keep_pings_i) - } - - toa_list_downsampled <- list( toa = inp_toa_list_all$toa[keep_pings, ], - sync_tag_idx_vec = inp_toa_list_all$sync_tag_idx[keep_pings], - epo_self_vec = inp_toa_list_all$epo_self_vec[keep_pings]) - return(toa_list_downsampled) -} - - -# Internal function to randomly downsample the toa-matrix for inp_sync -#' @param inp_toa_list_all Output from `getInpSyncToaList` -#' @inheritParams getInpSync -#' @noRd -downsampleToaList_random <- function(inp_toa_list_all, keep_rate){ - toa_list_downsampled <- list() - keeps_idx <- which(stats::rbinom(nrow(inp_toa_list_all$toa), 1, keep_rate) == 1) - toa_list_downsampled$toa <- inp_toa_list_all$toa[keeps_idx,] - toa_list_downsampled$epo_self_vec <- inp_toa_list_all$epo_self_vec[keeps_idx] - toa_list_downsampled$sync_tag_idx_vec <- inp_toa_list_all$sync_tag_idx_vec[keeps_idx] - - return(toa_list_downsampled) -} - - -#' Internal function to get vector of which hydros are fixed -#' @inheritParams getInpSync -#' @noRd -getFixedHydrosVec <- function(sync_dat, fixed_hydros_idx){ - fixed_hydros_vec <- rep(0, times=nrow(sync_dat$hydros)) - fixed_hydros_vec[fixed_hydros_idx] <- 1 - - return(fixed_hydros_vec) -} - - -#' Internal function to get info relating to sync offsets per day -#' @inheritParams getInpSync -#' @noRd -getOffsetVals <- function(inp_toa_list, n_offset_day){ - epo_self_vec <- inp_toa_list$epo_self_vec - epo_start <- min(epo_self_vec)-10 - epo_end <- max(epo_self_vec)+10 - - n_offset_idx <- ceiling((epo_end - epo_start)/(24*60*60)) * n_offset_day - offset_cuts <- cut(epo_self_vec, breaks=n_offset_idx, dig.lab=10) - offset_idx <- as.numeric(offset_cuts) - offset_labs <- levels(offset_cuts) - offset_levels <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", offset_labs) ), upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", offset_labs) )) - - offset_levels[1,1] <- epo_start - offset_levels[n_offset_idx,2] <- epo_end +# #' Internal function. Get toa for sync from sync_dat +# #' @inheritParams getInpSync +# #' @noRd +# getInpSyncToaList <- function(sync_dat, max_epo_diff, min_hydros, excl_self_detect, keep_rate, lin_corr_coeffs){ + # toa_list_gross <- buildToaListGross(sync_dat, excl_self_detect) + # toa_list_pruned <- pruneToaListGross(toa_list_gross, max_epo_diff, min_hydros) - dimnames(offset_levels) <- NULL - return(list(n_offset_idx=n_offset_idx, offset_idx=offset_idx, offset_levels=offset_levels)) -} - + # return(toa_list_pruned) +# } #' Internal function to get info relating to sync ss per day #' @inheritParams getInpSync @@ -167,48 +41,6 @@ getSsVals <- function(inp_toa_list, n_ss_day){ return(list(n_ss_idx=n_ss_idx, ss_idx=ss_idx, ss_levels=ss_levels)) } -#' Internal function to get params for TMB sync -#' @inheritParams getInpSync -#' @noRd -getParamsTmbSync <- function(dat_tmb_sync, ss_data_what){ - params_tmb_sync <- list() - - params_tmb_sync$OFFSET <- matrix(rnorm(dat_tmb_sync$nh*dat_tmb_sync$n_offset_idx, 0, 3), nrow=dat_tmb_sync$nh, ncol=dat_tmb_sync$n_offset_idx) - params_tmb_sync$SLOPE1 <- matrix(rnorm(dat_tmb_sync$nh*dat_tmb_sync$n_offset_idx, 0, 3), nrow=dat_tmb_sync$nh, ncol=dat_tmb_sync$n_offset_idx) - params_tmb_sync$SLOPE2 <- matrix(rnorm(dat_tmb_sync$nh*dat_tmb_sync$n_offset_idx, 0, 3), nrow=dat_tmb_sync$nh, ncol=dat_tmb_sync$n_offset_idx) - params_tmb_sync$TRUE_H <- as.matrix(cbind(dat_tmb_sync$H[,1], dat_tmb_sync$H[,2], dat_tmb_sync$H[,3])) - params_tmb_sync$LOG_SIGMA_TOA <- 0 - # LOG_SIGMA_HYDROS_XY = rnorm(dat_tmb_sync$nh,-3,1) - - if(ss_data_what == "est"){ - params_tmb_sync$SS <- rnorm(dat_tmb_sync$n_ss_idx, 1420, 1) - } - - if(dat_tmb_sync$sync_type == 'top'){ - params_tmb_sync$TOP <- rowMeans(dat_tmb_sync$toa_offset, na.rm=TRUE) - } - - return(params_tmb_sync) -} - -#' Internal function to get random params for TMB sync -#' @inheritParams getInpSync -#' @noRd -getRandomTmbSync <- function(dat_tmb_sync, ss_data_what){ - if(dat_tmb_sync$sync_type == "top"){ - random_tmb_sync <- c("TOP") - } else { - random_tmb_sync <- c() - } - - random_tmb_sync <- c(random_tmb_sync, "OFFSET", "SLOPE1", "SLOPE2", "TRUE_H") - - if(ss_data_what == "est"){ - random_tmb_sync <- c(random_tmb_sync, "SS") - } - return(random_tmb_sync) -} - #' Internal function to get residuals from sync_model in long format #' @inheritParams getInpSync #' @noRd diff --git a/R/syncHelperFuncs.R b/R/syncHelperFuncs.R index 436851f..2a57182 100644 --- a/R/syncHelperFuncs.R +++ b/R/syncHelperFuncs.R @@ -1,19 +1,3 @@ -#' Internal function. Get hydros for sync from sync_dat -#' @inheritParams getInpSync -#' @noRd -getInpSyncHInfo <- function(sync_dat){ - Hx0 <- sync_dat$hydros[1,x] - Hy0 <- sync_dat$hydros[1,y] - - inp_H <- sync_dat$hydros[, c('x','y','z')] - inp_H[, x:=x-Hx0] - inp_H[, y:=y-Hy0] - inp_H[] - - return(list(inp_H=inp_H, Hx0=Hx0, Hy0=Hy0)) -} - - #' Internal function to append needed columns to table detections diff --git a/R/sync_applyLinCor.R b/R/sync_applyLinCor.R new file mode 100644 index 0000000..1de5d95 --- /dev/null +++ b/R/sync_applyLinCor.R @@ -0,0 +1,27 @@ + +#' Internal function. Apply linear correction matrix to epofrac before sync +#' @inheritParams getInpSync +#' @noRd +#' @export +applyLinCor <- function(hydros, dat_sync, sync_params){ + if(is.na(sync_params$lin_corr[1])){ + lin_corr <- matrix(0, nrow=nrow(hydros), ncol=2, byrow=TRUE) + rownames(lin_corr) <- hydros$h_sn + attr(dat_sync, 'lin_corrected') <- FALSE + } else { + lin_corr <- sync_params$lin_corr + attr(dat_sync, 'lin_corrected') <- TRUE + } + + h_sns <- hydros$h_sn + + for(h in 1:length(h_sns)){ + sn_h <- h_sns[h] + dat_sync[h_sn == sn_h, epofrac := epofrac - lin_corr[rownames(lin_corr) == sn_h, 1] - epofrac * lin_corr[rownames(lin_corr) == sn_h, 2]] + } + + + dat_sync[] + return(dat_sync) +} + diff --git a/R/sync_checkInpSyncData.R b/R/sync_checkInpSyncData.R new file mode 100644 index 0000000..23a2881 --- /dev/null +++ b/R/sync_checkInpSyncData.R @@ -0,0 +1,56 @@ +#' Checks input data for sync models... +#' @export +checkInpSyncData <- function(hydros, dat_sync, dat_ss, sync_params){ + # # # data checks... + Es <- 0 + Ws <- 0 + + if(is.null(attr(hydros, 'yapsified'))){ + Es <- Es + 1 + stop("ERROR: The hydros data.table is not yapsified. Run e.g. hydros <- yapsifyHydros(hydros) to do so. \n") + } + + if(!is.null(attr(hydros, 'yapsified')) & attr(hydros, 'yapsified') == FALSE){ + Es <- Es + 1 + stop("ERROR: The hydros data.table is not yapsified. Run e.g. hydros <- yapsifyHydros(hydros) to do so. \n") + } + + if(!"epo" %in% colnames(dat_sync)){ + cat("Note: Adding column epo to dat_sync \n") + dat_sync[, epo := as.numeric(ts)] + } + + if(!"epofrac" %in% colnames(dat_sync)){ + cat("Note: Adding column epofrac to dat_sync \n") + dat_sync[, epofrac := epo+frac] + } + + + if(length(unique(hydros$h_sn)) != nrow(hydros)){ + Es <- Es+1 + stop("ERROR: At least one hydrophone serial number is used more than once in sync_dat$hydros!\n") + } + + if(sync_params$keep_rate <=0 | (sync_params$keep_rate > 1 & sync_params$keep_rate < 10) | (sync_params$keep_rate >= 10 & sync_params$keep_rate %% 1 != 0)){ + Es <- Es+1 + stop("ERROR: Invalid sync_params$keep_rate! Must be either ]0;1] or integer >= 10\n") + } + + if(sum(!hydros$h_sn %in% unique(dat_sync$h_sn)) != 0){ + det_hydros <- unique(dat_sync$h_sn) + not_dets <- hydros$h_sn[!hydros$h_sn %in% det_hydros] + cat(paste0("WARNING: These hydro serials were not found in the detection tabel. They cannot be synced. Please fix or remove. \n",paste(not_dets, collapse=", "),"\n")) + Ws <- Ws + 1 + } + + if(sum(!unique(dat_sync$h_sn) %in% hydros$h_sn) != 0){ + det_hydros <- unique(dat_sync$h_sn) + not_dets <- det_hydros[!hydros$h_sn %in% det_hydros] + cat(paste0("WARNING: These hydro serials are present in the detection tabel, but not in hydros. Please fix or remove. \n",paste(not_dets, collapse=", "),"\n")) + Ws <- Ws + 1 + } + + if(Es == 0 & Ws == 0){ + # cat("No errors or warnings found in InpSyncData\n") + } +} \ No newline at end of file diff --git a/R/sync_downsampleToaList.R b/R/sync_downsampleToaList.R new file mode 100644 index 0000000..b3c9b7f --- /dev/null +++ b/R/sync_downsampleToaList.R @@ -0,0 +1,81 @@ + +#' Internal warpper function to do downsampling of the toa +#' @inheritParams getInpSync +#' @noRd +#' @export +getDownsampledToaList <- function(inp_toa_list_all, offset_vals_all, keep_rate){ + if(keep_rate > 0 & keep_rate <= 1){ + toa_list_downsampled <- downsampleToaList_random(inp_toa_list_all, keep_rate) + } else if(keep_rate >= 10){ + toa_list_downsampled <- downsampleToaList_selective(inp_toa_list_all, offset_vals_all, keep_rate) + } + return(toa_list_downsampled) +} + +# Internal function to selectively downsample the toa-matrix for inp_sync +#' @param inp_toa_list_all Output from `getInpSyncToaList` +#' @inheritParams getInpSync +#' @noRd +#' @export +downsampleToaList_selective <- function(inp_toa_list_all, offset_vals_all, keep_rate){ + toa <- inp_toa_list_all$toa + offset_idx <- offset_vals_all$offset_idx + + toa_long <- data.table::data.table(reshape2::melt(toa), rep(offset_idx, times=ncol(toa))) + colnames(toa_long) <- c('ping', 'h_idx','toa','offset_idx') + # nobs_per_offset <- toa_long[!is.na(toa), .N, by=c('h_idx' ,'offset_idx')] + + nobs_per_offset <- data.table::data.table(reshape2::melt(with(toa_long[!is.na(toa)], table(h_idx, offset_idx)), value.name="N")) + + + keep_pings <- c() + for(i in 1:length(unique(offset_idx))){ + toa_long_i <- toa_long[offset_idx == i] + keep_pings_i <- c() + # h_order <- order(nobs_per_offset[offset_idx == i, N]) + + h_dets <- nobs_per_offset[offset_idx == i, N, by=h_idx] + data.table::setorder(h_dets, N) + h_order <- h_dets$h_idx + + for(h in 1:length(h_order)){ + already_in_keeps <- nrow(toa_long_i[!is.na(toa) & ping %in% keep_pings_i & offset_idx == i & h_idx==h_order[h]]) + if(already_in_keeps >= keep_rate){ + # already enough with this hydro... + next + } else { + need <- keep_rate - already_in_keeps + 1 + pings_h <- toa_long_i[!is.na(toa) & offset_idx == i & h_idx==h_order[h], ping] + if(length(pings_h) < need){ + keep_pings_h <- pings_h + } else { + keep_pings_h <- sample(pings_h, size=need) + } + keep_pings_i <- c(keep_pings_i, keep_pings_h) + } + } + keep_pings <- c(keep_pings, keep_pings_i) + } + + toa_list_downsampled <- list( toa = inp_toa_list_all$toa[keep_pings, ], + sync_tag_idx_vec = inp_toa_list_all$sync_tag_idx[keep_pings], + epo_self_vec = inp_toa_list_all$epo_self_vec[keep_pings]) + return(toa_list_downsampled) +} + + +# Internal function to randomly downsample the toa-matrix for inp_sync +#' @param inp_toa_list_all Output from `getInpSyncToaList` +#' @inheritParams getInpSync +#' @noRd +#' @export +downsampleToaList_random <- function(inp_toa_list_all, keep_rate){ + toa_list_downsampled <- list() + keeps_idx <- which(stats::rbinom(nrow(inp_toa_list_all$toa), 1, keep_rate) == 1) + toa_list_downsampled$toa <- inp_toa_list_all$toa[keeps_idx,] + toa_list_downsampled$epo_self_vec <- inp_toa_list_all$epo_self_vec[keeps_idx] + toa_list_downsampled$sync_tag_idx_vec <- inp_toa_list_all$sync_tag_idx_vec[keeps_idx] + + return(toa_list_downsampled) +} + diff --git a/R/sync_getDatTmbSync.R b/R/sync_getDatTmbSync.R new file mode 100644 index 0000000..8294485 --- /dev/null +++ b/R/sync_getDatTmbSync.R @@ -0,0 +1,47 @@ +#' Internal function to get dat for TMB sync +#' @inheritParams getInpSync +#' @noRd +#' @export +getDatTmbSync <- function(hydros, dat_sync, sync_params, inp_toa_list, offset_vals, T0, ss_data_vec){ + # H <- as.matrix(inp_H_info$inp_H) + H <- as.matrix(hydros[, .(x,y,z)]) + dimnames(H) <- NULL + + toa <- inp_toa_list$toa - T0 + toa_offset <- inp_toa_list$toa - offset_vals$offset_levels[offset_vals$offset_idx] + + if(sync_params$E_dist == 'Gaus' | sync_params$E_dist == 'gaus'){ + E_model <- 1 + } else if(sync_params$E_dist == 't'){ + E_model <- 2 + } + + time_keeper_idx <- hydros[h_sn == sync_params$time_keeper, h_idx] + + fixed_hydros_vec <- getFixedHydrosVec(hydros, sync_params) + + if(ss_data_vec[1] == 0){ + ss_data_how <- "est" + } else{ + ss_data_how <- 'data' + } + + + dat_tmb_sync <- list( + model = 'yaps_sync', + H=H, + toa = toa, + toa_offset = toa_offset, + sync_tag_idx_vec = inp_toa_list$sync_tag_idx_vec, + np = nrow(inp_toa_list$toa), + nh = ncol(inp_toa_list$toa), + tk = time_keeper_idx, + fixed_hydros_vec = fixed_hydros_vec, + offset_idx = offset_vals$offset_idx, + n_offset_idx = offset_vals$n_offset_idx, + ss_data_vec = ss_data_vec, + E_model = E_model + ) + return(dat_tmb_sync) +} + diff --git a/R/sync_getFixedHydrosVec.R b/R/sync_getFixedHydrosVec.R new file mode 100644 index 0000000..e0b03af --- /dev/null +++ b/R/sync_getFixedHydrosVec.R @@ -0,0 +1,17 @@ + +#' Internal function to get vector of which hydros are fixed +#' @inheritParams getInpSync +#' @noRd +#' @export +getFixedHydrosVec <- function(hydros, sync_params){ + + if(sync_params$fixed_hydros[1] == "all"){ + fixed_hydros_vec <- rep(1, times=nrow(hydros)) + } else { + fixed_hydros_vec <- rep(0, times=nrow(hydros)) + fixed_hydros_vec[hydros$h_sn %in% sync_params$fixed_hydros] <- 1 + } + + return(fixed_hydros_vec) +} + diff --git a/R/getInpSync.R b/R/sync_getInpSync.R similarity index 52% rename from R/getInpSync.R rename to R/sync_getInpSync.R index 9e43ffa..afe4273 100644 --- a/R/getInpSync.R +++ b/R/sync_getInpSync.R @@ -21,64 +21,51 @@ #' @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, sync_type='top'){ - 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") - } +# 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, sync_type='top', Edist_sync){ + +getInpSync <- function(hydros, dat_sync, dat_ss=NA, sync_params){ + # run some checks on the data going in... + checkInpSyncData(hydros, dat_sync, dat_ss, sync_params) + + # apply linear corrections if corrections values are provided - otherwise ignore + dat_sync <- applyLinCor(hydros, dat_sync, sync_params) - 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") - } + # get time 0 + t0 <- min(dat_sync$epo) - if(sum(!sync_dat$hydros$serial %in% unique(sync_dat$detections$serial)) != 0){ - det_hydros <- unique(sync_dat$detections$serial) - not_dets <- sync_dat$hydros$serial[!sync_dat$hydros$serial %in% det_hydros] - cat(paste0("WARNING: These hydro serials were not found in the detection tabel. They cannot be synced. Please fix or remove. \n",paste(not_dets, collapse=", "),"\n")) - } - if(sum(!unique(sync_dat$detections$serial) %in% sync_dat$hydros$serial) != 0){ - det_hydros <- unique(sync_dat$detections$serial) - not_dets <- det_hydros[!sync_dat$hydros$serial %in% det_hydros] - cat(paste0("WARNING: These hydro serials are present in the detection tabel, but not in hydros. Please fix or remove. \n",paste(not_dets, collapse=", "),"\n")) - } + inp_toa_list_all <- getSyncToa(hydros, dat_sync, sync_params) + offset_vals_all <- getOffsetVals(inp_toa_list = inp_toa_list_all, sync_params) - 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) + inp_toa_list <- getDownsampledToaList(inp_toa_list_all, offset_vals_all, keep_rate=sync_params$keep_rate) + offset_vals <- getOffsetVals(inp_toa_list, sync_params) - 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 { + + # tk_idx <- hydros[h_sn == sync_params$time_keeper, h_idx] + # init_offsets <- colMeans(inp_toa_list$toa - inp_toa_list$toa[, tk_idx], na.rm=TRUE) + + # gnu <- t(t(inp_toa_list$toa) - init_offsets) + # inp_toa_list$toa <- gnu + + if(is.na(dat_ss)){ ss_data_vec <- c(0) + } else { + ss_data_vec <- getSsDataVec(inp_toa_list, ss_data) # not upgraded to ver2 yet!!! } - dat_tmb_sync <- getDatTmbSync(sync_type, sync_dat, keep_rate, 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) - random_tmb_sync <- getRandomTmbSync(dat_tmb_sync, ss_data_what) - # 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 - ) + dat_tmb_sync <- getDatTmbSync(hydros, dat_sync, sync_params, inp_toa_list, offset_vals, t0, ss_data_vec) + params_tmb_sync <- getParamsTmbSync(dat_tmb_sync, ss_data_vec) + random_tmb_sync <- getRandomTmbSync(dat_tmb_sync, ss_data_vec) + + inits_tmb_sync <- c(-3, -3) # inits for LOG_SIGMA_TOA and LOG_SIGMA_OFFSET + + + + inp_params <- list(hydros=hydros, toa=inp_toa_list$toa, offset_levels=offset_vals$offset_levels, dat_ss=dat_ss, t0 = t0) - inp_sync <- list(sync_type=sync_type, 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) + 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, sync_params=sync_params) + # inp_sync$inp_params$sync_coverage <- checkInpSync(inp_sync, silent_check) return(inp_sync) } diff --git a/R/sync_getOffsetVals.R b/R/sync_getOffsetVals.R new file mode 100644 index 0000000..272136f --- /dev/null +++ b/R/sync_getOffsetVals.R @@ -0,0 +1,23 @@ + +#' Internal function to get info relating to sync offsets per day +#' @inheritParams getInpSync +#' @noRd +#' @export +getOffsetVals <- function(inp_toa_list, sync_params){ + epo_self_vec <- inp_toa_list$epo_self_vec + epo_start <- min(epo_self_vec)-10 + epo_end <- max(epo_self_vec)+10 + + n_offset_idx <- ceiling((epo_end - epo_start)/(24*60*60)) * sync_params$n_per_day + offset_cuts <- cut(epo_self_vec, breaks=n_offset_idx, dig.lab=10) + offset_idx <- as.numeric(offset_cuts) + offset_labs <- levels(offset_cuts) + offset_levels <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", offset_labs) ), upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", offset_labs) )) + + offset_levels[1,1] <- epo_start + offset_levels[n_offset_idx,2] <- epo_end + + dimnames(offset_levels) <- NULL + return(list(n_offset_idx=n_offset_idx, offset_idx=offset_idx, offset_levels=offset_levels)) +} + diff --git a/R/sync_getParamsTmbSync.R b/R/sync_getParamsTmbSync.R new file mode 100644 index 0000000..30d0d6a --- /dev/null +++ b/R/sync_getParamsTmbSync.R @@ -0,0 +1,22 @@ + +#' Internal function to get params for TMB sync +#' @inheritParams getInpSync +#' @noRd +#' @export +getParamsTmbSync <- function(dat_tmb_sync, ss_data_vec){ + params_tmb_sync <- list() + + params_tmb_sync$OFFSET <- matrix(rnorm(dat_tmb_sync$nh*dat_tmb_sync$n_offset_idx, 0, 3), nrow=dat_tmb_sync$nh, ncol=dat_tmb_sync$n_offset_idx) + params_tmb_sync$TRUE_H <- as.matrix(cbind(dat_tmb_sync$H[,1], dat_tmb_sync$H[,2], dat_tmb_sync$H[,3])) + params_tmb_sync$LOG_SIGMA_TOA <- 0 + params_tmb_sync$LOG_SIGMA_OFFSET <- 0 + + if(ss_data_vec[1] == 0){ + params_tmb_sync$SS <- rnorm(dat_tmb_sync$n_offset_idx, 1420, 1) + } + + params_tmb_sync$TOP <- rowMeans(dat_tmb_sync$toa_offset, na.rm=TRUE) + + return(params_tmb_sync) +} + diff --git a/R/sync_getRandomTmbSync.R b/R/sync_getRandomTmbSync.R new file mode 100644 index 0000000..2c63ff6 --- /dev/null +++ b/R/sync_getRandomTmbSync.R @@ -0,0 +1,15 @@ +#' Internal function to get random params for TMB sync +#' @inheritParams getInpSync +#' @noRd +#' @export +getRandomTmbSync <- function(dat_tmb_sync, ss_data_vec){ + random_tmb_sync <- c("TOP") + + # random_tmb_sync <- c(random_tmb_sync, "OFFSET", "SLOPE1", "SLOPE2", "TRUE_H") + random_tmb_sync <- c(random_tmb_sync, "OFFSET", "TRUE_H") + + if(ss_data_vec[1] == 0){ + random_tmb_sync <- c(random_tmb_sync, "SS") + } + return(random_tmb_sync) +} diff --git a/R/sync_getSyncToa.R b/R/sync_getSyncToa.R new file mode 100644 index 0000000..d4c3378 --- /dev/null +++ b/R/sync_getSyncToa.R @@ -0,0 +1,90 @@ + +#' Internal function to build sync TOA matrix, sync_tag_vec and epo_self_vec - the matrix needs pruning before use in sync function +#' @inheritParams getInpSync +#' @noRd +#' @export +# getSyncToa <- function(hydros, dat_sync, excl_self_detect, max_epo_diff, min_hydros){ +getSyncToa <- function(hydros, dat_sync, sync_params){ + + hs <- yapsifyHydros(hydros) + dets <- dat_sync + + dets <- merge(dets, hs[, .(h_sn, h_idx)]) + + sync_tags <- hs[!is.na(sync_tag), sync_tag] + + toa <- matrix(nrow=0, ncol=nrow(hs)) + colnames(toa) <- hs$h_sn + sync_tag_idx_vec <- c() + epo_self_vec <-c() + for(i in 1:length(sync_tags)){ + hydro_i_idx <- hs[sync_tag == sync_tags[i], h_idx] + sync_tag_i_idx <- hydro_i_idx + # hydro_st_serial <- hs[sync_tag == sync_tags[st], serial] + self_detections <- dets[tag==sync_tags[i] & h_idx==hydro_i_idx] + # detects_on_others <- detections[tag==sync_tags[i] & hydro_idx!=hydro_i_idx] + + # hack added to support ref tags to serve as beacons without a hydro + if(nrow(self_detections) == 0){ + num_detects_on_hydro <- dets[tag==sync_tags[i], .N, by=h_idx] + best_hydro_idx <- num_detects_on_hydro[which.max(num_detects_on_hydro$N), h_idx] + if(length(best_hydro_idx) == 0){ + self_detections <- NA #ugly hack added to support situations where hydros are completely ignored - e.g. due to very crappy performance + } else { + best_other_detections <- dets[tag==sync_tags[i] & h_idx==best_hydro_idx] + self_detections <- best_other_detections + } + } + + other_detections <- dets[tag==sync_tags[i] & h_idx!=hydro_i_idx] + + # if no other hydro detected the sync_tag it is useless... + if(nrow(other_detections) == 0){ + next + } + + self_detections[, epo_roll:=epofrac] + other_detections[, epo_roll:=epofrac] + toa_i <- t(plyr::daply(.data=other_detections, .variables="h_idx", .fun=function(k){ + k <- data.table::data.table(k) + return(as.numeric(k[self_detections, roll="nearest", on=.(epo_roll)]$epofrac)) + })) + + if(nrow(toa_i) == 1 ){ + toa_i <- t(toa_i) + } + + # idx <- unique(other_detections$hydro_idx) + # toa_i <- matrix(nrow=nrow(self_detections), ncol=length(idx)) + # colnames(toa_i) <- idx + # for(j in 1:length(idx)){ + # k <- other_detections[hydro_idx == idx[j]] + # toa_i[, j] <- + # as.numeric(k[self_detections, roll="nearest", on=.(epo_roll)]$epofrac) + # } + + toa_i_mat <- matrix(nrow=nrow(toa_i), ncol=nrow(hs)) + toa_i_mat[, as.numeric(colnames(toa_i))] <- toa_i + toa_i_mat[, hydro_i_idx] <- self_detections$epofrac + + if(sync_params$excl_self_detect){ + toa_i_mat[, hydro_i_idx] <- NA + } + + + toa <- rbind(toa, toa_i_mat) + sync_tag_idx_vec <- c(sync_tag_idx_vec, rep(sync_tag_i_idx, nrow(toa_i_mat))) + epo_self_vec <- c(epo_self_vec, self_detections$epofrac) + } + + # pruning the gross list... + toa[which(abs(toa - epo_self_vec) > sync_params$max_epo_diff)] <- NA + + nobs <- apply(toa, 1, function(k) sum(!is.na(k))) + keepers <- which(nobs >= sync_params$min_hydros) + toa <- toa[keepers ,] + sync_tag_idx_vec <- sync_tag_idx_vec[keepers] + epo_self_vec <- epo_self_vec[keepers] + + return(list(toa=toa, sync_tag_idx_vec=sync_tag_idx_vec, epo_self_vec=epo_self_vec)) +} diff --git a/R/yapsifyHydros.R b/R/yapsifyHydros.R new file mode 100644 index 0000000..f66ebdf --- /dev/null +++ b/R/yapsifyHydros.R @@ -0,0 +1,23 @@ +#' Internal function. yapsify / standarize hydros for sync and yaps +#' @noRd +#' @export +yapsifyHydros <- function(hydros){ + Hx0 <- hydros[1,x] + Hy0 <- hydros[1,y] + + std_h <- hydros[, c('h_sn', 'x','y','z','sync_tag')] + setorder(std_h, h_sn) + std_h[, h_idx := 1:.N] + + std_h[, x:=x-Hx0] + std_h[, y:=y-Hy0] + std_h[] + + attr(std_h, 'yapsified') <- TRUE + attr(std_h, 'Hx0') <- Hx0 + attr(std_h, 'Hy0') <- Hy0 + + return(std_h=std_h) +} + + diff --git a/README.Rmd b/README.Rmd index 4b7e80e..c973e10 100644 --- a/README.Rmd +++ b/README.Rmd @@ -26,9 +26,10 @@ knitr::opts_chunk$set( [![CRAN RStudio mirror downloads](https://cranlogs.r-pkg.org/badges/grand-total/yaps)](https://cran.r-project.org/package=yaps) - # YAPS - (Yet Another Positioning Solver) +YAPS was recently archived on CRAN due to dependence of a deprecated package (splusTimeSeries). We are currently working on an updated and improved version. + Welcome to the `yaps` repository. The `yaps` package is based on the original YAPS presented in Baktoft, Gjelland, Økland & Thygesen (2017): [Positioning of aquatic animals based on time-of-arrival and random walk models using YAPS (Yet Another Positioning Solver)](https://www.nature.com/articles/s41598-017-14278-z.pdf) To use `yaps` on own data, you need to compile a TOA-matrix based on synchronized hydrophone data and replace the hydros dataframe with actual hydrophone positions. A complete step-by-step guide on how to do this, can be found in our pre-print paper [Opening the black box of fish tracking using acoustic telemetry](https://www.biorxiv.org/content/10.1101/2019.12.16.877688v1). The example in this guide is based on data collected using a 69 kHz PPM-based system (Vemco VR2). We are working towards adding examples based on data collected using other manufacturers. diff --git a/man/checkInpSyncData.Rd b/man/checkInpSyncData.Rd new file mode 100644 index 0000000..ed835bd --- /dev/null +++ b/man/checkInpSyncData.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sync_checkInpSyncData.R +\name{checkInpSyncData} +\alias{checkInpSyncData} +\title{Checks input data for sync models...} +\usage{ +checkInpSyncData(hydros, dat_sync, dat_ss, sync_params) +} +\description{ +Checks input data for sync models... +} diff --git a/man/getInpSync.Rd b/man/getInpSync.Rd index 9e28507..ccd5351 100644 --- a/man/getInpSync.Rd +++ b/man/getInpSync.Rd @@ -1,25 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/getInpSync.R +% Please edit documentation in R/sync_getInpSync.R \name{getInpSync} \alias{getInpSync} \title{Get object inp for synchronization} \usage{ -getInpSync( - 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, - sync_type = "top" -) +getInpSync(hydros, dat_sync, dat_ss = NA, sync_params) } \arguments{ \item{sync_dat}{List containing data.tables with hydrophone information and detections. See e.g. \code{?ssu1} for example} @@ -36,6 +21,10 @@ getInpSync( \item{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.} +\item{ss_data_what}{Indicates whether to estimate ("est") speed of sound or to use data based on logged water temperature ("data").} + +\item{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).} + \item{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 @@ -45,10 +34,6 @@ OR (if keep_rate > 10) number of pings (approximate) to keep in each hydro X off \item{lin_corr_coeffs}{Matrix of coefficients used for pre-sync linear correction. \verb{dim(lin_corr_coeffs)=(#hydros, 2)}.} -\item{ss_data_what}{Indicates whether to estimate ("est") speed of sound or to use data based on logged water temperature ("data").} - -\item{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).} - \item{silent_check}{Logical whether to get output from \code{checkInpSync()}. Default is FALSE} \item{sync_type}{String specifying which type og sync model to use. Default sync_type='top' is the original included in yaps. Experimental sync_type='delta' is the ealternative - seems to works just as good, but might be much faster.} diff --git a/man/getSyncModel.Rd b/man/getSyncModel.Rd index 2c5aa0e..0c57117 100644 --- a/man/getSyncModel.Rd +++ b/man/getSyncModel.Rd @@ -4,23 +4,17 @@ \alias{getSyncModel} \title{Get sync model from inp_sync object obtained by \code{getInpSync()}} \usage{ -getSyncModel( - inp_sync, - silent = TRUE, - fine_tune = FALSE, - max_iter = 100, - tmb_smartsearch = TRUE -) +getSyncModel(inp_sync, silent = TRUE, max_iter = 500) } \arguments{ \item{inp_sync}{Input data prepared for the sync model using \code{getInpSync()}} \item{silent}{Keep TMB quiet} -\item{fine_tune}{Logical. Whether to re-run the sync model excluding residual outliers. \strong{Deprecated} use fineTuneSyncModel() instead.} - \item{max_iter}{Max number of iterations to run TMB. Default=100 seems to work in most cases.} +\item{fine_tune}{Logical. Whether to re-run the sync model excluding residual outliers. \strong{Deprecated} use fineTuneSyncModel() instead.} + \item{tmb_smartsearch}{Logical whether to use the TMB smartsearch in the inner optimizer (see \code{?TMB::MakeADFun} for info). Default and original implementation is TRUE. However, there seems to be an issue with some versions of \code{Matrix} that requires \code{tmb_smartsearch=FALSE}.} } \value{ diff --git a/src/nll_sync_ss_est.h b/src/nll_sync_ss_est.h index 283c4db..70b752d 100644 --- a/src/nll_sync_ss_est.h +++ b/src/nll_sync_ss_est.h @@ -1,11 +1,11 @@ // PARAMETER(logD_v); // Diffusivity of sound speed // Type D_v = exp(logD_v); - PARAMETER_VECTOR(ss); // Estimated speed of sound + PARAMETER_VECTOR(SS); // Estimated speed of sound - for(int i = 0; i < n_ss_idx; ++i){ - nll -= dnorm(ss(i), Type(1475.0), Type(100), true); + for(int i = 0; i < n_offset_idx; ++i){ + nll -= dnorm(SS(i), Type(1475.0), Type(100), true); } // nll -= dnorm(ss(0),Type(1475.0),Type(100.0),true); @@ -15,7 +15,7 @@ // } for(int i = 0; i < np; i++){ - ss_i(i) = ss(ss_idx(i)); + ss_i(i) = SS(offset_idx(i)); } diff --git a/src/yaps_sync.h b/src/yaps_sync.h index fda9824..f5269d1 100644 --- a/src/yaps_sync.h +++ b/src/yaps_sync.h @@ -1,52 +1,51 @@ - DATA_STRING(sync_type); DATA_ARRAY(H); + DATA_ARRAY(toa); DATA_ARRAY(toa_offset); - DATA_ARRAY(toa_delta); DATA_IVECTOR(sync_tag_idx_vec); DATA_INTEGER(np); DATA_INTEGER(nh); - DATA_INTEGER(ndelta); DATA_INTEGER(tk); DATA_VECTOR(fixed_hydros_vec); DATA_IVECTOR(offset_idx); DATA_INTEGER(n_offset_idx); - - DATA_STRING(ss_data_what); // 'est' = estimate SS-data; 'data' = Use SS-data + DATA_VECTOR(ss_data_vec); // Vector of SS-data if used - length(ss_data) = np - DATA_IVECTOR(ss_idx); - DATA_INTEGER(n_ss_idx); + + DATA_INTEGER(E_model); //1 => Gaus; 2 => t PARAMETER_ARRAY(OFFSET); - PARAMETER_ARRAY(SLOPE1); - PARAMETER_ARRAY(SLOPE2); - // PARAMETER_VECTOR(SS); + // // PARAMETER_ARRAY(SLOPE1); + // // PARAMETER_ARRAY(SLOPE2); + // // PARAMETER_VECTOR(SS); PARAMETER_ARRAY(TRUE_H); PARAMETER(LOG_SIGMA_TOA); Type SIGMA_TOA = exp(LOG_SIGMA_TOA); - // PARAMETER_VECTOR(LOG_SIGMA_HYDROS_XY); - // vector SIGMA_HYDROS_XY = exp(LOG_SIGMA_HYDROS_XY); + PARAMETER(LOG_SIGMA_OFFSET); + Type SIGMA_OFFSET = exp(LOG_SIGMA_OFFSET); + + // // PARAMETER_VECTOR(LOG_SIGMA_HYDROS_XY); + // // vector SIGMA_HYDROS_XY = exp(LOG_SIGMA_HYDROS_XY); - // if(sync_type == 'top'){ - // array mu(np,nh); - // array eps(np,nh); - // } else { - // vector mu(ndelta); - // vector eps(ndelta); - // } + // // if(sync_type == 'top'){ + // // array mu(np,nh); + // // array eps(np,nh); + // // } else { + // // vector mu(ndelta); + // // vector eps(ndelta); + // // } array dist_mat(nh,nh); vector ss_i(np); - // Adapting to zero-based index in cpp + // // Adapting to zero-based index in cpp sync_tag_idx_vec = sync_tag_idx_vec - 1; offset_idx = offset_idx - 1; - ss_idx = ss_idx - 1; tk = tk - 1; - // Type nll = 0.0; + // // Type nll = 0.0; parallel_accumulator nll(this); for(int h1=0; h1 < nh; ++h1){ @@ -59,10 +58,10 @@ } } - if(ss_data_what != "est"){ - #include "nll_sync_ss_data.h" - } else { + if(ss_data_vec(0) == 0){ #include "nll_sync_ss_est.h" + } else { + #include "nll_sync_ss_data.h" } @@ -71,18 +70,22 @@ for(int h=0; h mu(np,nh); @@ -92,39 +95,42 @@ for(int p = 0; p < np; ++p){ // iterate pings for(int h=0; h < nh; ++h){ // iterate hydros in ping p if(!isNA(toa_offset(p,h))){ - mu(p,h) = TOP(p) + dist_mat(sync_tag_idx_vec(p), h)/ss_i(p) + OFFSET(h, offset_idx(p)) + SLOPE1(h, offset_idx(p))*(toa_offset(p,h)/1E6) + SLOPE2(h, offset_idx(p))*pow((toa_offset(p,h)/1E6),2); + mu(p,h) = TOP(p) + dist_mat(sync_tag_idx_vec(p), h)/ss_i(p) + OFFSET(h, offset_idx(p)); eps(p,h) = toa_offset(p,h) - mu(p,h); - // nll -= dnorm(eps(p,h), Type(0.0), SIGMA_TOA, true); - // nll -= log(dt(eps(p,h)/SCALE, Type(3.0), false)/SCALE); - nll -= log(dt(eps(p,h)/SIGMA_TOA, Type(3.0), false)/SIGMA_TOA); + + if(E_model == 1){ //gaus sync model + nll -= dnorm(eps(p,h), Type(0.0), SIGMA_TOA, true); + } else if(E_model == 2){ //t sync model + nll -= log(dt(eps(p,h)/SIGMA_TOA, Type(3.0), false)/SIGMA_TOA); + } } } } REPORT(eps); - } else { - vector mu(ndelta); - vector eps(ndelta); - - - for(int d = 0; d < ndelta; ++d){ // iterate pings - int H1 = CppAD::Integer(toa_delta(d,0)) - 1; - int H2 = CppAD::Integer(toa_delta(d,1)) - 1; - int ping_idx = CppAD::Integer(toa_delta(d,3)) - 1; - int sync_tag_idx = CppAD::Integer(toa_delta(d,4)) - 1; - int offset_idx = CppAD::Integer(toa_delta(d,5)) -1; - int ss_idx = CppAD::Integer(toa_delta(d,6)) - 1; - - // mu(d) = (dist_mat(sync_tag_idx, H1)/ss_i(ping_idx) + OFFSET(H1, 0) + SLOPE1(H1, 0)*(toa(ping_idx, H1)/1E6) + SLOPE2(H1, 0)*pow(toa(ping_idx, H1)/1E6, 2) + SLOPE3(H1, 0)*pow(toa(ping_idx, H1)/1E6, 3)) - - // (dist_mat(sync_tag_idx, H2)/ss_i(ping_idx) + OFFSET(H2, 0) + SLOPE1(H2, 0)*(toa(ping_idx, H2)/1E6) + SLOPE2(H2, 0)*pow(toa(ping_idx, H2)/1E6, 2) + SLOPE3(H2, 0)*pow(toa(ping_idx, H2)/1E6, 3)); - mu(d) = (dist_mat(sync_tag_idx, H1)/ss_i(ping_idx) + OFFSET(H1, offset_idx) + SLOPE1(H1, offset_idx)*(toa_offset(ping_idx, H1)/1E6) + SLOPE2(H1, offset_idx)*pow(toa_offset(ping_idx, H1)/1E6, 2)) - - (dist_mat(sync_tag_idx, H2)/ss_i(ping_idx) + OFFSET(H2, offset_idx) + SLOPE1(H2, offset_idx)*(toa_offset(ping_idx, H2)/1E6) + SLOPE2(H2, offset_idx)*pow(toa_offset(ping_idx, H2)/1E6, 2)); - - eps(d) = toa_delta(d, 2) - mu(d); - nll -= log(dt(eps(d)/SIGMA_TOA, Type(3.0), false)/SIGMA_TOA); - - } - REPORT(eps); - } + // // } else { + // // vector mu(ndelta); + // // vector eps(ndelta); + + + // // for(int d = 0; d < ndelta; ++d){ // iterate pings + // // int H1 = CppAD::Integer(toa_delta(d,0)) - 1; + // // int H2 = CppAD::Integer(toa_delta(d,1)) - 1; + // // int ping_idx = CppAD::Integer(toa_delta(d,3)) - 1; + // // int sync_tag_idx = CppAD::Integer(toa_delta(d,4)) - 1; + // // int offset_idx = CppAD::Integer(toa_delta(d,5)) -1; + // // int ss_idx = CppAD::Integer(toa_delta(d,6)) - 1; + + // // // mu(d) = (dist_mat(sync_tag_idx, H1)/ss_i(ping_idx) + OFFSET(H1, offset_idx) + SLOPE1(H1, offset_idx)*(toa_offset(ping_idx, H1)/1E6) + SLOPE2(H1, offset_idx)*pow(toa_offset(ping_idx, H1)/1E6, 2)) - + // // // (dist_mat(sync_tag_idx, H2)/ss_i(ping_idx) + OFFSET(H2, offset_idx) + SLOPE1(H2, offset_idx)*(toa_offset(ping_idx, H2)/1E6) + SLOPE2(H2, offset_idx)*pow(toa_offset(ping_idx, H2)/1E6, 2)); + // // mu(d) = (dist_mat(sync_tag_idx, H1)/ss_i(ping_idx) + OFFSET(H1, offset_idx) ) - + // // (dist_mat(sync_tag_idx, H2)/ss_i(ping_idx) + OFFSET(H2, offset_idx) ); + + // // eps(d) = toa_delta(d, 2) - mu(d); + // // nll -= log(dt(eps(d)/SIGMA_TOA, Type(3.0), false)/SIGMA_TOA); + + // // } + // // REPORT(eps); + // // } for(int h=0; h