init commit of v.2.0.0.9000

baktoft · Sep 13, 2023 · 5a60d11 · 5a60d11
1 parent 42c8c0a
commit 5a60d11
Show file tree

Hide file tree

Showing 24 changed files with 562 additions and 390 deletions.
diff --git 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 = "[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")),

diff --git 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
@@ -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
diff --git 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,22 +62,18 @@ 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))
 	sds <- summary(jointrep)[,2]
 	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
@@ -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
@@ -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
@@ -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)
+}
+