From 01450a613c68706f87ea927addeb63da040e91a1 Mon Sep 17 00:00:00 2001 From: Henrik Baktoft Date: Fri, 1 Feb 2019 14:29:24 +0100 Subject: [PATCH] Use ss data (#6) * Add option to use measured ss_data * Update example in readme.md --- R/prepTmb.R | 12 +++++++++--- README.md | 14 +++++++++----- man/getDatTmb.Rd | 7 ++++++- man/getInp.Rd | 6 +++++- src/yaps.cpp | 11 +++++++++-- 5 files changed, 38 insertions(+), 12 deletions(-) diff --git a/R/prepTmb.R b/R/prepTmb.R index 59a04f2..4c4af4d 100644 --- a/R/prepTmb.R +++ b/R/prepTmb.R @@ -8,12 +8,14 @@ #' @param pingType Type of transmitter to simulate - either stable burst interval ('sbi') or random burst interval ('rbi') #' @param rbi_min,rbi_max Minimum and maximum BI for random burst interval transmitters #' @param sdInits If >0 initial values will be randomized around the normally fixed value using rnorm(length(inits), mean=inits, sd=sdInits) +#' @param ss_data_what What speed of sound (ss) data to be used. Default ss_data_what='est': ss is estimated by the model. Alternatively, if ss_data_what='data': ss_data must be provided and length(ss_data) == ncol(toa) +#' @param ss_data Vector of ss-data to be used if ss_data_what = 'est'. Otherwise ss_data <- 0 (default) #' @return List of input data ready for use in TMB-call #' @export -getInp <- function(hydros, toa, E_dist, n_ss, pingType, sdInits=1, rbi_min=0, rbi_max=0){ - datTmb <- getDatTmb(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max) +getInp <- function(hydros, toa, E_dist, n_ss, pingType, sdInits=1, rbi_min=0, rbi_max=0, ss_data_what='est', ss_data=0){ + datTmb <- getDatTmb(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max, ss_data_what, ss_data) params <- getParams(datTmb) inits <- getInits(pingType, sdInits) return(list( @@ -31,7 +33,7 @@ getInp <- function(hydros, toa, E_dist, n_ss, pingType, sdInits=1, rbi_min=0, rb #' #' @return List for use in TMB. #' @export -getDatTmb <- function(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max){ +getDatTmb <- function(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max, ss_data_what, ss_data){ if(n_ss > 1){ ss_idx <- cut(1:ncol(toa), n_ss, labels=FALSE) - 1 #-1 because zero-indexing in TMB } else { @@ -39,6 +41,8 @@ getDatTmb <- function(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max){ } approxBI <- mean(diff(toa[1,]), na.rm=TRUE) + if(ss_data_what == 'data') { stopifnot(length(ss_data) == ncol(toa))} + Edist <- rep(0,2) if(E_dist == "Gaus") {Edist[1] <- 1} if(E_dist == "Mixture") {Edist[2] <- 1} @@ -56,6 +60,8 @@ getDatTmb <- function(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max){ pingType = pingType, n_ss = n_ss, ss_idx = ss_idx, + ss_data_what = ss_data_what, + ss_data = ss_data, approxBI = approxBI ) return(datTmb) diff --git a/README.md b/README.md index 8456161..05526a0 100644 --- a/README.md +++ b/README.md @@ -26,10 +26,10 @@ trueTrack <- simTrueTrack(model='crw', n = 15000, deltaTime=1, shape=1, scale=0. # Simulate telemetry observations from true track. # Format and parameters depend on type of transmitter burst interval (BI) - stable (sbi) or random (rbi). -pingType <- 'rbi' +pingType <- 'sbi' if(pingType == 'sbi') { # stable BI - sbi_mean <- 5; sbi_sd <- 1e-4; + sbi_mean <- 30; sbi_sd <- 1e-4; teleTrack <- simTelemetryTrack(trueTrack, ss='rw', pingType=pingType, sbi_mean=sbi_mean, sbi_sd=sbi_sd) } else if(pingType == 'rbi'){ # random BI pingType <- 'rbi'; rbi_min <- 20; rbi_max <- 40; @@ -41,10 +41,15 @@ hydros <- simHydros(auto=TRUE, trueTrack=trueTrack) toa_list <- simToa(teleTrack, hydros, pingType, sigmaToa=1e-4, pNA=0.25, pMP=0.01) toa <- toa_list$toa +# Specify whether to use ss_data from measured water temperature (ss_data_what <- 'data') or to estimate ss in the model (ss_data_what <- 'est') +ss_data_what <- 'data' +if(ss_data_what == 'data') {ss_data <- teleTrack$ss} else {ss_data <- 0} + + if(pingType == 'sbi'){ - inp <- getInp(hydros, toa, E_dist="Mixture", n_ss=10, pingType=pingType, sdInits=0) + inp <- getInp(hydros, toa, E_dist="Mixture", n_ss=10, pingType=pingType, sdInits=0, ss_data_what=ss_data_what, ss_data=ss_data) } else if(pingType == 'rbi'){ - inp <- getInp(hydros, toa, E_dist="Mixture", n_ss=10, pingType=pingType, sdInits=0, rbi_min=rbi_min, rbi_max=rbi_max) + inp <- getInp(hydros, toa, E_dist="Mixture", n_ss=10, pingType=pingType, sdInits=0, rbi_min=rbi_min, rbi_max=rbi_max, ss_data_what=ss_data_what, ss_data=ss_data) } str(inp) @@ -64,5 +69,4 @@ plot(y~x, data=trueTrack, type="l", xlim=range(hydros$hx), ylim=range(hydros$hy) lines(y~x, data=teleTrack) points(hy~hx, data=hydros, col="green", pch=20, cex=3) lines(pl$Y~pl$X, col="red") - ``` diff --git a/man/getDatTmb.Rd b/man/getDatTmb.Rd index a351dca..51c7818 100644 --- a/man/getDatTmb.Rd +++ b/man/getDatTmb.Rd @@ -4,7 +4,8 @@ \alias{getDatTmb} \title{Get data for input to TMB} \usage{ -getDatTmb(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max) +getDatTmb(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max, + ss_data_what, ss_data) } \arguments{ \item{hydros}{Dataframe from simHydros() or Dataframe with columns hx and hy containing positions of the receivers. Translate the coordinates to get the grid centre close to (0;0).} @@ -20,6 +21,10 @@ getDatTmb(hydros, toa, E_dist, n_ss, pingType, rbi_min, rbi_max) \item{rbi_min}{Minimum and maximum BI for random burst interval transmitters} \item{rbi_max}{Minimum and maximum BI for random burst interval transmitters} + +\item{ss_data_what}{What speed of sound (ss) data to be used. Default ss_data_what='est': ss is estimated by the model. Alternatively, if ss_data_what='data': ss_data must be provided and length(ss_data) == ncol(toa)} + +\item{ss_data}{Vector of ss-data to be used if ss_data_what = 'est'. Otherwise ss_data <- 0 (default)} } \value{ List for use in TMB. diff --git a/man/getInp.Rd b/man/getInp.Rd index e6a9135..8bbb586 100644 --- a/man/getInp.Rd +++ b/man/getInp.Rd @@ -5,7 +5,7 @@ \title{Get prepared inp-object for use in TMB-call} \usage{ getInp(hydros, toa, E_dist, n_ss, pingType, sdInits = 1, rbi_min = 0, - rbi_max = 0) + rbi_max = 0, ss_data_what = "est", ss_data = 0) } \arguments{ \item{hydros}{Dataframe from simHydros() or Dataframe with columns hx and hy containing positions of the receivers. Translate the coordinates to get the grid centre close to (0;0).} @@ -21,6 +21,10 @@ getInp(hydros, toa, E_dist, n_ss, pingType, sdInits = 1, rbi_min = 0, \item{sdInits}{If >0 initial values will be randomized around the normally fixed value using rnorm(length(inits), mean=inits, sd=sdInits)} \item{rbi_min, rbi_max}{Minimum and maximum BI for random burst interval transmitters} + +\item{ss_data_what}{What speed of sound (ss) data to be used. Default ss_data_what='est': ss is estimated by the model. Alternatively, if ss_data_what='data': ss_data must be provided and length(ss_data) == ncol(toa)} + +\item{ss_data}{Vector of ss-data to be used if ss_data_what = 'est'. Otherwise ss_data <- 0 (default)} } \value{ List of input data ready for use in TMB-call diff --git a/src/yaps.cpp b/src/yaps.cpp index 8254a8a..9c5cf5d 100644 --- a/src/yaps.cpp +++ b/src/yaps.cpp @@ -29,6 +29,8 @@ Type objective_function::operator() () DATA_SCALAR(bi_penalty); DATA_SCALAR(rbi_min); DATA_SCALAR(rbi_max); + DATA_STRING(ss_data_what); // 'est' = estimate SS-data; 'data' = Use SS-data + DATA_VECTOR(ss_data); // Vector of SS-data if used - length(ss_data) = np DATA_IVECTOR(ss_idx); DATA_INTEGER(n_ss); DATA_SCALAR(approxBI); @@ -38,7 +40,7 @@ Type objective_function::operator() () PARAMETER_VECTOR(X); //Position at time of ping PARAMETER_VECTOR(Y); //Position at time of ping PARAMETER_VECTOR(top); // Estimated time of pings - PARAMETER_VECTOR(ss); + PARAMETER_VECTOR(ss); // Estimated speed of sound PARAMETER(logD_xy); // Diffusivity of fish Type D_xy = exp(logD_xy); @@ -69,7 +71,11 @@ Type objective_function::operator() () for(int h=0; h::operator() () nll -= dnorm(logSigma_toa, Type(0), Type(25), true); nll -= dnorm(logScale, Type(0), Type(25), true); nll -= dnorm(logSigma_bi, Type(0), Type(25), true); + nll -= dnorm(logD_v, Type(0), Type(25), true); //position component nll -= dnorm(X(0),Type(0),Type(100),true);