#======================================================================
#                    L A D C P P R O C . B E S T L A G 
#                    doc: Tue Sep 28 21:58:48 2010
#                    dlm: Fri Mar  6 15:53:43 2015
#                    (c) 2010 A.M. Thurnherr
#                    uE-Info: 33 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

# TODO:
#	- first lag is always(?) nan unless CTD is turned on with LADCP in water

# HISTORY:
#	Sep 28, 2010: - created
#	Dec  9, 2010: - adapted to %CTD
#	Dec 10, 2010: - hardened bestlag failure test to require 1/3 agreeing lags
#	Jan  5, 2011: - changed first guess from 80% down to 10% down
#				  - added LADCP time lag to %PARAMs
#				  - added support of -i
#	Jul  7, 2011: - added code to remove window-mean of w before lagging to
#				    make it work in regions of crazy ocean w (IWISE 16007)
#	Jul 15, 2011: - changed screen-output of lag to take first guess lag into
#					account
#	Apr 11, 2011: - removed 1st guess lag consistency check based on large
#					elapsed offsets
#	May 18, 2012: - BUG: window start index was not always calculated correctly
#	Oct 19, 2012: - BUG: opt_i had wrong sign!
#	Jun 25, 2013: - adapted to :: %PARAM convention
#	Mar 19, 2014: - moved %PARAM to LADCPproc
#	Jul 19, 2014: - made lagging obey -z)oom
#	May 25, 2015: - added assertion to require numeric interpolated LADCP_w
#				  - BUG: interp_LADCP_w left gaps
#				  - added debug code to output bestLag input time series

sub interp_LADCP_w($$)
{
	my($elapsed,$ens) = @_;
	my($sc) = ($elapsed - $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME}) /
			  ($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} -
					$LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME});
	unless (numberp($LADCP{ENSEMBLE}[$ens]->{W})) {
		$nGaps++;
		$LADCP{ENSEMBLE}[$ens]->{W} = $LADCP{ENSEMBLE}[$ens-1]->{W};
	}
	return $LADCP{ENSEMBLE}[$ens-1]->{W} +
				$sc * ($LADCP{ENSEMBLE}[$ens]->{W} - $LADCP{ENSEMBLE}[$ens-1]->{W});
}

sub bestLag($)
{
	my($ws) = @_;														# window start index

	my($best);
	my($bestmad) = 9e99;												# mean absolute deviation
	for (my($Llag)=-int($opt_w/2); $Llag<int($opt_w/2); $Llag++) {
		my($mCw,$mLw,$nw) = (0,0,0);									# first calc means
		for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
			my($Li) = $Ci + $Llag;
			next if ($Li<0 || $Li>=$opt_w);
			next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
			$mCw += $CTD{w}[$ws+$Ci];
			$mLw += $LADCP_w[$ws+$Li];
			$nw++;
		}
		next unless ($nw > 0);
		$mCw /= $nw; $mLw /= $nw;
		
		my($sad) = my($nad) = 0;										# calc mad with means removed
		for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
			my($Li) = $Ci + $Llag;
			next if ($Li<0 || $Li>=$opt_w);
			next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
			$sad += abs($CTD{w}[$ws+$Ci]-$mCw - ($LADCP_w[$ws+$Li]-$mLw));
			$nad++;
		}
		if ($sad/$nad < $bestmad) {
			$best = $Llag;
			$bestmad = $sad/$nad;
		}
	}
	return $best;	
}

sub lagLADCP2CTD()
{
	#------------------------------------------------------------------------
	# find 1st rec & ensemble >=10% down to max depth & make 1st guess at lag
	#------------------------------------------------------------------------
	
	my($first_guess_lag);											# in units of CTD records
	
	if (defined($opt_i)) {
		$first_guess_lag = $opt_i / $CTD{sampint};
	} else {
		my($CTD_10pct_down) = 0;
		$CTD_10pct_down++											# "until" formulation allows for missing pressures
			until ($CTD{press}[$CTD_10pct_down]-$CTD{press}[0] >= 0.1*($CTD{maxpress}-$CTD{press}[0]));
	    
		my($LADCP_10pct_down) = 0;
		$LADCP_10pct_down++
			while ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH} < 0.1*$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
	    
		$first_guess_lag = ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED_TIME} -
							   $CTD_10pct_down*$CTD{sampint}) / $CTD{sampint};
	    
		printf(STDERR "\n\t1st guess offset [CTD pressure, LADCP estimated depth] = %ds [%ddbar, %dm]\n",
				$first_guess_lag*$CTD{sampint},$CTD{press}[$CTD_10pct_down],$LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH})
					if ($opt_d);
	}
	
	#------------------------------------------------------------------------------------
	# Linearly interpolate LADCP time series onto a new grid with $CTD{sampint} resolution
	#	ALSO: apply first_guess_lag to make lags small, which keeps the bestlag data
	#		  chunks large
	#------------------------------------------------------------------------------------
	
	$nGaps = 0;
	
	for (my($ens)=$LADCP_start,my($r)=0; $ens<=$LADCP_end; $ens++) {
		while ($r*$CTD{sampint} < $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}) {
			unless ($first_guess_lag > $r) {
				my($w) = interp_LADCP_w($r*$CTD{sampint},$ens);
				next if (!defined($firstValid) && !defined($w));
				$firstValid = $r-$first_guess_lag unless defined($firstValid);
				die("assertion failed") unless defined($w);
				$LADCP_w[$r-$first_guess_lag] = $w;
				$nValid++;
			}
			$r++;
		}
	}
	
	print(STDERR "\t$nGaps gaps in w timeseries")
		if ($opt_d);
	
	print(STDERR "\n");

	#----------------------------------------------------------------------
	# Output w Time Series
	#----------------------------------------------------------------------

#	open(F,'>bestLag.out');
#	print(F "#ANTS#FIELDS# {rec} {LADCP_w} {CTD_w}\n");
#	for (my($r)=$firstValid; $r<$firstValid+$nValid; $r++) {
#		print(F "$r $LADCP_w[$r] $CTD{w}[$r]\n");
#	}
#	close(F);

	#----------------------------------------------------------------------
	# Calculate lags
	#----------------------------------------------------------------------

	printf(STDERR "\tcalculating $opt_n lags from %ds-long windows [s]:",$opt_w);
	$opt_w = int($opt_w / $CTD{sampint});

	#---------------------------------------------------------------
	# carry out opt_n lag correlations and keep tally of the results
	#---------------------------------------------------------------
	my(%nBest);
	my($nLags) = 0;
	my($lags) = '';
	my($lastLag) = 9e99; my($nSame) = 1;
	for (my($window)=0; $window<$opt_n; $window++) {
		my($ws) = $firstValid + $window * int($nValid/$opt_n);	# window start
		$ws = @LADCP_w-$opt_w if ($ws+$opt_w >= @LADCP_w);
		$bestLag = bestLag($ws);
		if (defined($bestLag)) {
			if (defined($lastLag) && $bestLag == $lastLag) {
				$nSame++;
			} else {
				printf(STDERR "(x%d)",$nSame)
					if ($nSame > 1);
				printf(STDERR " %d",$bestLag*$CTD{sampint});
				$nSame = 1;
				$lastLag = $bestLag;
			}
			$lags .= sprintf(" %s",$bestLag*$CTD{sampint});
			$nBest{$bestLag}++;
			$nLags++;
		} else {
			if (!defined($lastLag)) {
				$nSame++;
			} else {
				printf(STDERR "(x%d)",$nSame)
					if ($nSame > 1);
				printf(STDERR " nan");
				$nSame = 1;
				$lastLag = $bestLag;
			}
		}
	}
	printf(STDERR "(x%d)",$nSame) if ($nSame > 1);
    &antsAddParams('LADCPproc::time_lags',$lags);
	
	#----------------------
	# find most popular lag
    #----------------------
	my($best_lag);
	foreach my $i (keys(%nBest)) {
		$best_lag = $i if ($nBest{$i} > $nBest{$best_lag});
	}
	croak("\n$0: cannot determine a valid lag\n")
		unless ($opt_z || $nBest{$best_lag}>$opt_n/3);
	print(STDERR "\n\n\t\tWARNING: only $nBest{$best_lag} of the lag estimates agree!\n")
		if ($nBest{$best_lag} < $opt_n/2);

	if ($nBest{$best_lag} == $nLags) {
		printf(STDERR "\n\t\tunanimous lag = %ds\n",($first_guess_lag+$best_lag)*$CTD{sampint});
	} else {
		printf(STDERR "\n\t\tmost popular lag = %ds\n",($first_guess_lag+$best_lag)*$CTD{sampint});
	}

	return ($first_guess_lag + $best_lag) * $CTD{sampint};
}

1;