too many changes...

baktoft · Sep 13, 2023 · 42c8c0a · 42c8c0a
1 parent 17d079e
commit 42c8c0a
Show file tree

Hide file tree

Showing 19 changed files with 71 additions and 44 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -12,7 +12,7 @@ License: GPL-3
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.1.1
+RoxygenNote: 7.2.1
 LinkingTo: Rcpp, TMB, RcppEigen
 Imports: circular, cowplot, data.table, ggplot2, ggrepel, nloptr, plyr, Rcpp, reshape2, splusTimeSeries, stats, TMB, viridis, zoo
 Suggests: 

diff --git a/R/prepTmb.R b/R/prepTmb.R
@@ -96,7 +96,7 @@ getParams <- function(datTmb){
 		  # X = 0 + stats::rnorm(ncol(datTmb$toa), sd=10)
 		# , Y = 0 + stats::rnorm(ncol(datTmb$toa), sd=10)
 		, top = zoo::na.approx(apply(datTmb$toa, 2, function(k) {stats::median(k, na.rm=TRUE)}), rule=2)	#time of ping
-		# , ss=stats::rnorm(datTmb$n_ss, 1450, 5) 	#speed of sound
+		# , SS=stats::rnorm(datTmb$n_ss, 1450, 5) 	#speed of sound
 		, logD_xy = 0				#diffusivity of transmitter movement (D_xy in ms)
 		# , logD_v = 0				#diffusivity of speed of sound (D_v in ms)
 		# , logSigma_toa = 0			#sigma for Gaussian
@@ -114,7 +114,7 @@ getParams <- function(datTmb){
 	# # # ss related
 	if(datTmb$ss_data_what == 'est'){
 		out$logD_v <- 0				#diffusivity of speed of sound (D_v in ms)
-		out$ss <- stats::rnorm(datTmb$n_ss, 1450, 5) 	#speed of sound
+		out$SS <- stats::rnorm(datTmb$n_ss, 1450, 5) 	#speed of sound
 	}
 
 	# # # Edist related

diff --git a/R/runYaps.R b/R/runYaps.R
@@ -38,7 +38,7 @@ runYaps <- function(inp, maxIter=1000, getPlsd=TRUE, getRep=TRUE, silent=TRUE, o
 		random <- c(random, "Z")
 	}	
 	if(inp$datTmb$ss_data_what == 'est'){
-		random <- c(random, 'ss')
+		random <- c(random, 'SS')
 	}
 
 	if(inp$datTmb$pingType == 'pbi'){

diff --git a/man/applySync.Rd b/man/applySync.Rd
diff --git a/man/checkInp.Rd b/man/checkInp.Rd
diff --git a/man/checkInpSync.Rd b/man/checkInpSync.Rd
diff --git a/man/fineTuneSyncModel.Rd b/man/fineTuneSyncModel.Rd
diff --git a/man/getInp.Rd b/man/getInp.Rd
diff --git a/man/getInpSync.Rd b/man/getInpSync.Rd
diff --git a/man/getSyncCoverage.Rd b/man/getSyncCoverage.Rd
diff --git a/man/getSyncModel.Rd b/man/getSyncModel.Rd
diff --git a/man/getToaYaps.Rd b/man/getToaYaps.Rd
diff --git a/man/runYaps.Rd b/man/runYaps.Rd
diff --git a/man/simTelemetryTrack.Rd b/man/simTelemetryTrack.Rd
diff --git a/src/nll_mix-dist.h b/src/nll_mix-dist.h
@@ -7,19 +7,28 @@
 
 	PARAMETER(log_t_part);		// t-part of mixture model 
 	// Type t_part = exp(log_t_part);
+
+	nll -= dnorm(logSigma_toa, Type(-5), Type(5), true);
+	nll -= dnorm(logScale, Type(-5), Type(5), true);
+	nll -= dnorm(log_t_part, Type(-5), Type(5), true);
+
+
 	Type t_part = invlogit(log_t_part);
 	Type G_part = Type(1.0) - t_part; //Gaussian part of mixture model
 
 
 	for(int i=0; i<np; ++i){ //iterate pings
 		for(int h=0; h<nh; ++h){ //iterate hydros
 			if(!isNA(toa(h,i))){ //ignore NA's...
-				nll -= Edist(1) * log( G_part * dnorm(eps(h,i), Type(0),sigma_toa,false) + 		//Gaussian part
-					t_part * dt(eps(h,i)/scale, Type(3.0), false) );					//t part
-					// t_part * dt(eps(h,i)/scale, Type(3.0), false) / scale );					//t part
-					// t_part * dt(eps(h,i), Type(3.0), false));					//t part
+				nll -= Edist(1) * 
+							// ( 1 / ( 1 + 1E3 * (softplus(Type(0.0) - delta_t(h,i), Type(1E-6)))) ) *	// Attempt to keep toas before top out of estimation
+//							( 1 / ( 1 + 1E3 * (	)) ) * 
+							log( 
+								G_part * dnorm(eps(h,i), Type(0),sigma_toa,false) + 		//Gaussian part
+								t_part * dt(eps(h,i)/scale, Type(3.0), false) / scale				//t part
+							);
+						// t_part * dt(eps(h,i)/scale, Type(3.0), false) / scale );					//t part
+						// t_part * dt(eps(h,i), Type(3.0), false));					//t part
 			}
 		}
 	}
-
-
diff --git a/src/nll_pingtype_sbi.h b/src/nll_pingtype_sbi.h
@@ -2,9 +2,17 @@
 	PARAMETER(logSigma_bi);		// Sigma for burst interval
 	Type sigma_bi = exp(logSigma_bi);
 
-	nll -= dnorm(top(0),Type(0.0),Type(4.0),true);
-	nll -= dnorm(top(1),Type(approxBI),Type(4.0),true);
-	for(int i = 2; i < np; ++i)	{
-		nll -= dnorm(top(i)-2*top(i-1)+top(i-2), Type(0), sigma_bi, true);
-		// nll += bi_penalty * (softplus(0 - (top(i) - top(i-1)), bi_epsilon));
+	// nll -= dnorm(top(0),Type(0.0),Type(4.0),true);
+	// nll -= dnorm(top(1),Type(approxBI),Type(4.0),true);
+	// for(int i = 2; i < np; ++i)	{
+		// nll -= dnorm(top(i)-2*top(i-1)+top(i-2), Type(0), sigma_bi, true);
+		// // nll += bi_penalty * (softplus(0 - (top(i) - top(i-1)), bi_epsilon));
+	// }
+
+
+	nll -= dnorm(top(0), Type(0), Type(0.1), true);
+	nll -= dnorm(top(0) - (-1 * top(2) + 2*top(1)), Type(0), exp(logSigma_bi), true);
+	nll -= dnorm(top(1) - (-1 * top(3) + 2*top(2)), Type(0), exp(logSigma_bi), true);
+	for(int i=2; i < np; i++){
+	    nll -= dnorm(top(i)-2*top(i-1)+top(i-2), Type(0), exp(logSigma_bi), true);
 	}
diff --git a/src/nll_ss_est.h b/src/nll_ss_est.h
@@ -1,17 +1,20 @@
 	PARAMETER(logD_v);    		// Diffusivity of sound speed
 	Type D_v = exp(logD_v);
+	nll -= dnorm(logD_v, Type(0), Type(5), true);
+
+	PARAMETER_VECTOR(SS);		// Estimated speed of sound
 
-	PARAMETER_VECTOR(ss);		// Estimated speed of sound
-
-	nll -= dnorm(ss(0),Type(1475.0),Type(100.0),true);		
+	// nll -= dnorm(SS(0),Type(1475.0),Type(100.0),true);		
+	for(int i = 0; i < n_ss; ++i){
+		nll -= dnorm(SS(i),Type(1475.0),Type(100.0),true);	
+	}
 	for(int i = 1; i < n_ss; ++i){
-		nll -= dnorm(ss(i), ss(i-1),sqrt(2*D_v), true);
-		// nll -= dnorm(ss(i),Type(1475.0),Type(100.0),true);	
+		nll -= dnorm(SS(i), SS(i-1),sqrt(2*D_v), true);
 	}
 	for(int i = 0; i < n_ss; ++i){
-		nll += Type(1E6) * (softplus(ss(i) - 1600, Type(1E-6)) + softplus(1400 - ss(i), Type(1E-6)));
+		nll += Type(1E6) * (softplus(SS(i) - 1600, Type(1E-6)) + softplus(1400 - SS(i), Type(1E-6)));
 	}
 
 	for(int i = 0; i < np; i++){
-		ss_i(i) = ss(ss_idx(i));
+		ss_i(i) = SS(ss_idx(i));
 	}
diff --git a/src/nll_sync_ss_est.h b/src/nll_sync_ss_est.h
@@ -1,20 +1,21 @@
 	// 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 < np; i++){
-		ss_i(i) = SS(ss_idx(i));
+		nll -= dnorm(ss(i), Type(1475.0), Type(100), true);
 	}
 
 	// nll -= dnorm(ss(0),Type(1475.0),Type(100.0),true);		
 	// for(int i = 1; i < n_ss; ++i){
 		// nll -= dnorm(ss(i), ss(i-1),sqrt(2*D_v), true);
 		// nll -= dnorm(ss(i),Type(1475.0),Type(100.0),true);		
 	// }
+
+	for(int i = 0; i < np; i++){
+		ss_i(i) = ss(ss_idx(i));
+	}
+
 
diff --git a/src/yaps_track.h b/src/yaps_track.h
@@ -30,6 +30,7 @@
 
 
 	array<Type> mu_toa(nh,np);  // mu-matrix
+	array<Type> delta_t(nh,np);  // matrix for delta toa - top
 	array<Type> dist(nh,np);	// dist-matrix
 	array<Type> eps(nh,np);		// eps-matrix
 	vector<Type> ss_i(np);
@@ -53,6 +54,11 @@
 					}
 					mu_toa(h,i) = top(i) +  dist(h,i)/ss_i(i);
 					eps(h,i) = toa(h,i) - mu_toa(h,i);
+
+					// // Making sure that all toas are later than top
+					// delta_t(h,i) = (toa(h,i) - top(i));
+					// nll += bi_penalty * (softplus(0 - delta_t(h,i), bi_epsilon));
+					// nll += Type(1E6) * (softplus(delta_t(h,i) - 2, Type(1E-6)) + softplus(-1 - delta_t(h,i), Type(1E-6)));
 				}
 			}
 		}
@@ -75,6 +81,10 @@
 	for(int i=1; i<np; ++i)	{
 		nll -= dnorm(X(i), X(i-1),sqrt(2*D_xy*(top(i) - top(i-1))),true);	
 		nll -= dnorm(Y(i), Y(i-1),sqrt(2*D_xy*(top(i) - top(i-1))),true);
+
+		// nll -= dnorm(X(i),Type(0),Type(10000),true);
+		// nll -= dnorm(Y(i),Type(0),Type(10000),true);
+
 	}
 
 	// spatial constraints