Use ss data (#6)

* Add option to use measured ss_data

* Update example in readme.md

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,14 +33,16 @@ 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 {
 		ss_idx <- rep(0, ncol(toa))
 	}
 	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)

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")
-
 ```

src/yaps.cpp

@@ -29,6 +29,8 @@ Type objective_function<Type>::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<Type>::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<Type>::operator() ()
 		for(int h=0; h<nh; ++h){ //iterate hydros
 			if(!isNA(toa(h,i))){ //ignore NA's...
 				dist(h,i) = sqrt((H(h,0)-X(i))*(H(h,0)-X(i)) + (H(h,1)-Y(i))*(H(h,1)-Y(i)));
-				mu_toa(h,i) = top(i) +  dist(h,i)/ss(ss_idx(i));
+				if(ss_data_what == "est"){
+					mu_toa(h,i) = top(i) +  dist(h,i)/ss(ss_idx(i));
+				} else {
+					mu_toa(h,i) = top(i) +  dist(h,i)/ss_data(i);
+				}
 				Type eps = toa(h,i) - mu_toa(h,i);
 
 				nll -= Edist(0) * dnorm(eps, Type(0), sigma_toa, true); 					//Gaussian part					
@@ -85,6 +91,7 @@ Type objective_function<Type>::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);