merge dev into master (#11)

* add option to use multiple tmb-models

* add basic yaps_sync functionality - seems to be working ok

* more functions added

* finalizing and cleaning up #1

* restructuring separate tmb-models into one

* finalizing and cleaning up #2

* finalizing and cleaning up #3

* added KØG suggestions to readme

* cleanup

* fix bug re merge.data.table

* allow slight out-of-bounds rbis

* add basic plotting

* fix issue with basic plotting

* fix issue with basic plotting

* Update syncPlotters.R

minor modifications to label plotting

* Update syncPlotters.R

* Update syncPlotters.R

* Update syncPlotters.R

modification of tick labs in plotSyncModelCheck

* Update syncPlotters.R

* suppressWarnings when silent=true

* re-compile example data ssu1

* replacing readme

* up version

* fix error in readme example

.Rbuildignore

@@ -1,2 +1,4 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^README\.Rmd$
+^README_sync\.Rmd$

.gitignore

@@ -31,6 +31,10 @@ vignettes/*.pdf
 # Temporary files created by R markdown
 *.utf8.md
 *.knit.md
+man/figures/
+/man/figures/
+man/figures/*.png
+/man/figures/*.png
 
 # Shiny token, see https://shiny.rstudio.com/articles/shinyapps.html
 rsconnect/

DESCRIPTION

@@ -1,12 +1,16 @@
 Package: yaps
 Title: Track estimation using YAPS (Yet Another Positioning Solver)
-Version: 1.1.1
-Authors@R: person("Henrik", "Baktoft", email = "[email protected]", role = c("aut", "cre"))
+Version: 1.2.0.9100
+Authors@R: c(	person("Henrik", "Baktoft", email = "[email protected]", role = c("cre", "aut")),
+				person("Karl", "Gjelland", role=c("aut")),
+				person("Uffe H.", "Thygesen", role=c("aut")),
+				person("Finn", "Økland", role=c("aut"))
+			)
 Description: Estimate track from acoustic data using YAPS (doi:10.1038/s41598-017-14278-z).
 Depends: R (>= 3.5.0)
 License: GPL-3 + file LICENSE
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 6.1.1
 LinkingTo: Rcpp, TMB, RcppEigen
-Imports: Rcpp, TMB, circular, RcppEigen, zoo, stats
+Imports: circular, data.table, ggplot2, plyr,Rcpp, RcppEigen, reshape2, splusTimeSeries, stats, tictoc, TMB, viridis,zoo

NAMESPACE

@@ -1,13 +1,59 @@
 # Generated by roxygen2: do not edit by hand
 
+export(applySync)
 export(getDatTmb)
 export(getInits)
 export(getInp)
+export(getInpSync)
 export(getParams)
+export(getSyncModel)
+export(getToaYaps)
+export(plotSyncModelCheck)
+export(plotSyncModelResids)
+export(plotYaps)
+export(prepDetections)
 export(runTmb)
+export(runYaps)
 export(simHydros)
 export(simTelemetryTrack)
 export(simToa)
 export(simTrueTrack)
 importFrom(Rcpp,sourceCpp)
-useDynLib(yaps, .registration = TRUE)
+importFrom(data.table,"%between%")
+importFrom(data.table,"%like%")
+importFrom(data.table,":=")
+importFrom(data.table,.I)
+importFrom(data.table,.N)
+importFrom(data.table,.SD)
+importFrom(ggplot2,aes)
+importFrom(ggplot2,facet_grid)
+importFrom(ggplot2,facet_wrap)
+importFrom(ggplot2,geom_bar)
+importFrom(ggplot2,geom_boxplot)
+importFrom(ggplot2,geom_histogram)
+importFrom(ggplot2,geom_hline)
+importFrom(ggplot2,geom_line)
+importFrom(ggplot2,geom_point)
+importFrom(ggplot2,geom_ribbon)
+importFrom(ggplot2,geom_smooth)
+importFrom(ggplot2,geom_tile)
+importFrom(ggplot2,geom_vline)
+importFrom(ggplot2,labs)
+importFrom(ggplot2,scale_fill_gradientn)
+importFrom(ggplot2,theme)
+importFrom(ggplot2,xlab)
+importFrom(ggplot2,ylab)
+importFrom(grDevices,dev.off)
+importFrom(grDevices,pdf)
+importFrom(graphics,abline)
+importFrom(graphics,lines)
+importFrom(graphics,matplot)
+importFrom(graphics,par)
+importFrom(graphics,plot)
+importFrom(graphics,points)
+importFrom(stats,median)
+importFrom(stats,power)
+importFrom(stats,quantile)
+importFrom(stats,rnorm)
+importFrom(utils,tail)
+useDynLib(yaps)

R/applySync.R

@@ -0,0 +1,83 @@
+#' Apply sync model to toa matrix to obtain synced data
+#' @param toa Object containing data to be synchronized. Typically a `data.table` as e.g. `ssu1$detections`, but can also be a matrix dim=(n_ping, n_hydo).
+#' @param hydros data.table formatted as `ssu1$hydros`
+#' @param sync_model Synchronization model obtained using `getSyncModel()`
+#' @export
+applySync <- function(toa, hydros="", sync_model){
+	if(is.matrix(toa)) {type <- "toa_matrix"}
+	else if(data.table::is.data.table(toa)) {type <- "detections_table"}
+
+	inp_synced <- sync_model$inp_synced
+
+	ks <- inp_synced$inp_params$offset_levels[, 1]
+	ks[1] <- ks[1] - inp_synced$inp_params$max_epo_diff
+
+
+	if(type=="toa_matrix"){
+		offset_idx_mat <- matrix(findInterval(toa, ks), ncol=ncol(toa))
+		offset_level_mat <- matrix(inp_synced$inp_params$offset_levels[offset_idx_mat, 1], ncol=ncol(offset_idx_mat))
+
+		toa_offset <- toa - offset_level_mat
+		offset_hydro_idx <- as.matrix(reshape2::melt(offset_idx_mat, value.name="idx")[,c(2,3)])
+
+		offset_mat <- matrix(sync_model$pl$OFFSET[offset_hydro_idx], ncol=ncol(toa))
+		slope_mat <- matrix(sync_model$pl$SLOPE1[offset_hydro_idx], ncol=ncol(toa))
+		slope2_mat <- matrix(sync_model$pl$SLOPE2[offset_hydro_idx], ncol=ncol(toa))
+
+		toa_synced <- offset_level_mat + toa_offset - offset_mat - slope_mat*toa_offset/1E6 - slope2_mat*(toa_offset/1E6)^2
+	}
+
+	if(type=="detections_table"){
+
+		if(!'epofrac' %in% colnames(toa)) {toa[, epofrac:=epo+frac]}
+		if(!'hydro_idx' %in% colnames(toa)){
+			toa[, hydro_idx := merge(toa, hydros[, c('serial','idx')], by='serial', sort=FALSE)$idx]
+		}
+
+		sync_dt <- data.table::data.table()
+		sync_dt[, epofrac := toa$epofrac]
+		sync_dt[, hydro_idx := toa$hydro_idx]
+		sync_dt[, id:=1:.N]
+		sync_dt[, offset_idx:=findInterval(toa$epofrac, ks)]
+		# NA those epofracs outside sync_period, i.e. offset_idx outside range 1:length(ks)
+		sync_dt[!offset_idx %in% 1:length(ks), 'offset_idx'] <- NA
+		sync_dt[, offset_level:= inp_synced$inp_params$offset_levels[offset_idx,1] ]
+		# sync_dt[, offset_hydro_idx:=toa$hydro_idx]
+
+		OFFSET_long <- data.table::data.table(reshape2::melt(sync_model$pl$OFFSET))
+		colnames(OFFSET_long) <- c('hydro_idx', 'offset_idx', 'OFFSET')
+		SLOPE1_long <- data.table::data.table(reshape2::melt(sync_model$pl$SLOPE1))
+		colnames(SLOPE1_long) <- c('hydro_idx', 'offset_idx', 'SLOPE1')
+		SLOPE2_long <- data.table::data.table(reshape2::melt(sync_model$pl$SLOPE2))
+		colnames(SLOPE2_long) <- c('hydro_idx', 'offset_idx', 'SLOPE2')
+
+		sync_dt <- merge(sync_dt, OFFSET_long, sort=FALSE, all.x=TRUE)
+		sync_dt <- merge(sync_dt, SLOPE1_long, sort=FALSE, all.x=TRUE)
+		sync_dt <- merge(sync_dt, SLOPE2_long, sort=FALSE, all.x=TRUE)
+
+		# sync_dt[, sync_model$pl$OFFSET[hydro_idx, offset_idx]]
+
+		# sync_dt[, SLOPE1:=sync_model$pl$SLOPE1[hydro_idx, offset_idx]]
+		# sync_dt[, SLOPE2:=sync_model$pl$SLOPE2[hydro_idx, offset_idx]]
+
+		# table(sync_dt$hydro_idx)
+		# table(sync_dt$offset_idx)
+		# table(sync_dt$OFFSET)
+		# table(sync_dt$OFFSET, sync_dt$hydro_idx)
+
+
+		sync_dt[, eposync := epofrac - OFFSET - SLOPE1*(epofrac - offset_level)/1E6 - SLOPE2*(((epofrac - offset_level)/1E6)^2)]
+		sync_dt[, slope1 := SLOPE1*(epofrac - offset_level)/1E6]
+		sync_dt[, slope2 := SLOPE2*(((epofrac - offset_level)/1E6)^2)]
+
+
+		toa[, eposync := sync_dt[, eposync]]
+		# toa[tag==5138]
+
+		toa_synced <- toa
+		toa_synced[]
+	}
+
+	return(toa_synced)
+}
+

R/data.R

@@ -0,0 +1,65 @@
+#' Test data from Florida Bay
+#'
+#' Small data set collected for positioning using acoustic telemetry and YAPS. \cr
+#' The data are part of a feasibility study using YAPS on Vemco PPM style data to track fish in shallow parts of Florida Bay. Data were collected using VR2 (Vemco) hydrophones. \cr
+#' Included in yaps with permission from J.S. Rehage, FIU Florida International University.
+#'
+#' @format A list containing 3 data.tables:
+#' \describe{
+#'   \item{hydros}{
+#'     \itemize{
+#'        \item serial Hydrophone serial number.
+#'        \item x,y,z Position of hydrophones in UTM.
+#'        \item sync_tag ID of co-located sync tag. Must be identical to entries in data.table detections$tag.
+#'        \item idx Unique values from 1:nrow(hydros).
+#'      }
+#'   }
+#'   \item{detections}{
+#'     \itemize{
+#'        \item ts Timestamp of detection in POSIXct().
+#'        \item tag ID of detected tag.
+#'        \item epo Timestamp as number of seconds since Unix epoch. Can be obtained using as.numeric(ts).
+#'        \item frac Sub-second part of detection timestamp in fractions of second [0-1].
+#'        \item serial Serial number of detecting hydrophone. Must match entry in data.table hydros.
+#'      }
+#'   }
+#'   \item{gps}{
+#'     \itemize{
+#'        \item ts Timestamp of gps position in POSIXct().
+#'        \item utm_x, utm_y Coordinates of position. Same projection and coordinate system as used in hydros.
+#'      }
+#'   }
+#' }
+"ssu1"
+
+#' Test data from Lake Hald, Denmark
+#'
+#' 
+#' 
+#' @format A list containing 3 data.tables:
+#' \describe{
+#'   \item{hydros}{
+#'     \itemize{
+#'        \item serial Hydrophone serial number.
+#'        \item x,y,z Position of hydrophones in UTM.
+#'        \item sync_tag ID of co-located sync tag. Must be identical to entries in data.table detections$tag.
+#'        \item idx Unique values from 1:nrow(hydros).
+#'      }
+#'   }
+#'   \item{detections}{
+#'     \itemize{
+#'        \item ts Timestamp of detection in POSIXct().
+#'        \item tag ID of detected tag.
+#'        \item epo Timestamp as number of seconds since Unix epoch. Can be obtained using as.numeric(ts).
+#'        \item frac Sub-second part of detection timestamp in fractions of second [0-1].
+#'        \item serial Serial number of detecting hydrophone. Must match entry in data.table hydros.
+#'      }
+#'   }
+#'   \item{gps}{
+#'     \itemize{
+#'        \item ts Timestamp of gps position in POSIXct().
+#'        \item utm_x, utm_y Coordinates of position. Same projection and coordinate system as used in hydros.
+#'      }
+#'   }
+#' }
+"hald"

R/getToaYaps.R

@@ -0,0 +1,79 @@
+#' Build TOA matrix from synced data.table - also do some pre-filtering of severe MP, pruning loose ends etc
+#' 
+#' @param synced_dat `data.table` containing synchronized data formatted as output from/or obtained using `applySync()`
+#' @inheritParams getInp
+#' @export
+getToaYaps <- function(synced_dat, hydros, rbi_min, rbi_max){
+
+	# remove multipaths...
+	data.table::setkey(synced_dat, hydro_idx, eposync)
+	diffs <- synced_dat[, c(diff(eposync),NA), by=hydro_idx]$V1
+	mps <- which(diffs < .5)+1
+	if(length(mps) > 0){
+		synced_dat <- synced_dat[-mps]
+	}
+
+	# build toa-matrix...
+	ts_focal <- splusTimeSeries::signalSeries(pos=floor(synced_dat[hydro_idx==1,eposync*10]), data=synced_dat[hydro_idx==1,eposync*10])
+	for(i in 2:nrow(hydros)){
+		# print(i)
+		xi <- splusTimeSeries::signalSeries(pos=floor(synced_dat[hydro_idx==i,eposync*10]), data=synced_dat[hydro_idx==i,eposync*10])
+		ts_focal <- splusTimeSeries::seriesMerge(ts_focal, xi, pos="union", matchtol=10)
+	}
+	ts_focal <- as.matrix(as.data.frame(ts_focal))
+	ts_focal <- ts_focal/10
+	toa <- ts_focal
+	dim(toa)
+	dimnames(toa) <- NULL
+
+	# remove rows with too short BI...
+	top1 <- rowMeans(toa, na.rm=TRUE)
+	diffs1 <- c(diff(top1),NA)
+	nobs <- apply(toa, 1, function(k) sum(!is.na(k)))
+	rem_idx <- which(diffs1 < rbi_min-1) # THIS NEEDS TO BE SET BASED ON USED SYSTEM - 1 IS TOO HIGH FOR HR-LIKE
+	if(length(rem_idx) > 0){
+		toa[rem_idx, ] <- NA
+	}
+
+	# remove empty rows... 
+	empty_rows <- which(apply(toa, 1, function(k) sum(!is.na(k))) == 0)
+	if(length(empty_rows) > 0){
+		toa <- toa[-empty_rows, ]
+	}
+
+	# trim toa to exclude rows in start and end with very few obs
+	nobs <- apply(toa, 1, function(k) sum(!is.na(k)))
+	first_ping <- 	which(nobs >= 2)[1]
+	last_ping <- 	rev(which(nobs >= 2))[1]
+	toa <- toa[first_ping:last_ping, ]
+
+	# remake toa-matrix to include pings missed by all hydros...
+	# top2 <- apply(toa, 1, function(k) median(k, na.rm=TRUE))
+	top2 <- rowMeans(toa, na.rm=TRUE) # we use rowMeans instead of apply(... median) - rowMeans is much faster...
+	diffs2 <- c(diff(top2),NA)
+
+	pings <- data.table::data.table(top=top2, diff=diffs2)
+	pings[, toa_idx:=1:.N]
+	pings[, ping2next := 1]
+	pings[, next_ping_too_late := diff > rbi_max+1]
+	pings[next_ping_too_late==TRUE, ping2next:=ping2next+round(diff/rbi_max)]
+	pings[, ping_idx:=cumsum(c(1,ping2next[-.N]))]
+	toa_all <- matrix(ncol=ncol(toa), nrow=max(pings$ping_idx))
+	toa_all[pings$ping_idx, ] <- toa
+
+
+	# pings[, ping2next:=round(diffs2 / ((rbi_max - rbi_min)/2 + rbi_min))]
+	# pings[, min_ping2next:=round(diffs2 / rbi_max)]
+	# pings[, max_ping2next:=round(diffs2 / rbi_min)]
+
+	# top3 <- rowMeans(toa_all, na.rm=TRUE) # we use rowMeans instead of apply(... median) - rowMeans is much faster...
+	# diffs3 <- c(diff(top3),NA)
+	# plot(diffs3)
+	# which(diffs3 > 31)
+
+	# toa_all[110:115,]
+
+	return(toa_all)
+
+
+}

R/plotYaps.R

@@ -0,0 +1,30 @@
+#' Basic plots of yaps output
+#' 
+#' @param inp Input object obtained using \code{getInp()}
+#' @param yaps_out Output from succesful run of \code{runYaps()}
+#' @param type Plot type. \code{type="map"} prodces a basic map of estimated track and hydrophones; \code{type="coord_X"}, \code{type="coord_Y"} produces plots of X and Y coordinated including +- 1 standard error.
+#' @export
+plotYaps <- function(inp, yaps_out, type="map"){
+	pl <- yaps_out$pl
+	plsd <- yaps_out$plsd
+
+	pl$X <- pl$X + inp$inp_params$Hx0
+	pl$Y <- pl$Y + inp$inp_params$Hy0
+	pl$top <- pl$top + inp$inp_params$T0
+
+	hydros <- data.frame(hx=inp$datTmb$H[,1] + inp$inp_params$Hx0, hy=inp$datTmb$H[,2] + inp$inp_params$Hy0)
+
+	if(type=="map"){
+		plot(hy~hx, data=hydros, col="green", pch=20, cex=2, asp=1, xlab="UTM_X", ylab="UTM_Y")
+		lines(pl$Y~pl$X, col="red")
+	} else if(type=="coord_X"){
+		plot(pl$top, pl$X, type="l", col="red", ylab="UTM_X", xlab="TimeOfPing")
+		lines(pl$top, pl$X-plsd$X, col="red", lty=3)
+		lines(pl$top, pl$X+plsd$X, col="red", lty=3)
+	} else if(type == "coord_Y"){
+		plot(pl$top, pl$Y, type="l", col="red", ylab="UTM_Y", xlab="TimeOfPing")
+		lines(pl$top, pl$Y-plsd$Y, col="red", lty=3)
+		lines(pl$top, pl$Y+plsd$Y, col="red", lty=3)
+	}
+
+}

R/prepFiles.R

@@ -0,0 +1,17 @@
+#' Prepare detections data.table from raw data - csv-files exported from vendor software
+#' @param raw_dat Data file from vendor supplied software
+#' @param type Type of the vendor file. Currently only 'vemco_vue' is supported.
+#' @export
+prepDetections <- function(raw_dat, type){
+	detections <- data.table::data.table()
+	if (type == "vemco_vue"){
+		detections[, ts:=as.POSIXct(raw_dat$'Date and Time (UTC)', tz="UTC")]
+		detections[, tag:=as.numeric(sapply(raw_dat$Transmitter, function(x) strsplit(x, "-")[[1]][3]))]
+		detections[, epo:=as.numeric(ts)]
+		detections[, frac:=as.numeric(sapply(raw_dat$'Date and Time (UTC)', function(x) strsplit(x, "\\.")[[1]][2]))]
+		detections[, serial:=as.numeric(sapply(raw_dat$Receiver, function(x) strsplit(x, "-")[[1]][2]))]
+	}
+	detections[]
+	return(detections)
+}
+