archive.txt

######################################################################################################
######################################################################################################

# big loop to find optimal values for stuff (starting with temperature increases)
def get_stats(data, consec_points):
	
	# determining the depth of actual HB flag:
	for i in range(0,len(flags.flag)):
			if (flags.flag[i] == "HB"):
				HB_depth = flags.depth[i]
			else:
 				continue
	
	try:
		# calculate peak accuracy (spike function returns int 0 if no peak found)
		temp_peak = temp_increase(data, consec_points) 
		temp_depth = temp_peak[0][0]
		difference = HB_depth - temp_depth
	
		# categorising (1=good, 0=bad)
		if (abs(difference) < 5):
			temp_stats = 1
		else:
			temp_stats = 0
	except:
		temp_stats = 0
		pass

	return(temp_stats)

# big loop to find optimal values for stuff (starting with temperature increases)
def get_stats_all(data, consec_points):
	
	# determining the depth of actual HB flag:
	for i in range(0,len(flags.flag)):
			if (flags.flag[i] == "HB"):
				HB_depth = flags.depth[i]
			else:
 				continue
	
	try:
		# calculate peak accuracy (spike function returns int 0 if no peak found)
		peaks = temp_increase(data, consec_points) 
		count = 0
	
		# categorising (1=good, 0=bad)
		for j in range(0,len(peaks)):
			difference = peaks[j][0] - HB_depth
			if (abs(difference) < 5):
				count = count + 1
	except:
		count = 0
		pass

	return(count)

# for loop to go over all data
consec_points = np.arange(10,160,10)
t_consec = []
tall_consec = []

for j in range(0,len(consec_points)):
	print(j)
	total_count = 0 
	total_all_counts = 0
	for i in range(0,len(name_array)):
		print(i)
	
		# reading in file here
		filename = name_array[i]
		read_data(filename)
		total_count = total_count + get_stats(data, consec_points[j])
		total_all_counts = total_all_counts + get_stats_all(data, consec_points[j])
		
	# store in bigger array to know which is best
	t_consec.append(total_count)
	tall_consec.append(total_all_counts)

# printing overall counts per number of consecutive points
print(t_consec)
print(tall_consec)

######################################################################################################
######################################################################################################

def get_stats(data, gradient, consec_points, detection_threshold):
	
	# determining the depth of actual HB flag:
	for i in range(0,len(flags.flag)):
			if (flags.flag[i] == "HB"):
				HB_depth = flags.depth[i]
			else:
 				continue
	
	try:
		# calculate peak accuracy (spike function returns int 0 if no peak found)
		peak = const_temp(data, gradient, consec_points, detection_threshold)
		temp_stats = 0	
	
		for j in range(0,len(peak)):
			difference = peak[j][0] - HB_depth
			if (abs(difference) < 5):
				temp_stats = temp_stats + 1
	except:
		temp_stats = 0
		pass

	return(temp_stats)

# for loop to go over all data
detection_threshold = np.logspace(-8,-1,8)
all_count = []

for j in range(0,len(detection_threshold)):
	print(j)
	total_count = 0 
	for i in range(0,len(name_array)):
		print(i)
	
		# reading in file here
		filename = name_array[i]
		read_data(filename)
		total_count = total_count + get_stats(data, gradient, 100, detection_threshold[j])
		
	# store in bigger array to know which is best
	all_count.append(total_count)

# printing overall counts per number of consecutive points
print(all_count)



######################################################################################################
######################################################################################################

# plotting additonal features
		consecpts = const_temp(data,gradient,100,0.001)
		plt.plot(consecpts[:,1],consecpts[:,0],'go')
		grow = temp_increase(data,50)
		plt.plot(grow[:,1],grow[:,0],'bo')	
		bath_z = bath_depth(latitude, longitude, bath_lon, bath_lat, bath_height)
		plt.axhline(y=bath_z, hold=None, color='g')	
		spikes = grad_spike(data,gradient, 3)
		if (type(spikes) != int):
			plt.plot(spikes[:,1], spikes[:,0],'ro')
		small_spikes = T_spike(data, 0.05)
		plt.plot(small_spikes[:,1],small_spikes[:,0],'yo')		


######################################################################################################
######################################################################################################

	# looking at the second derivative
	d2Tdz2 = []
	depth_secgrad = []
	
	for jj in range(0,n-2):
		depth_secgrad.append(gradient[jj][0])
		der = float((gradient[jj+1][1]-gradient[jj][1])/(gradient[jj+1][1]-gradient[jj][1]))
		d2Tdz2.append(der)
	secDer = np.column_stack((depth_secgrad,d2Tdz2))

	# taking a moving average for the temperature values
	depth_9pt = []
	temp9pt = []
	n = len(gradient[:,0])
	for jj in range(4,n-4):
		depth_9pt.append(gradient[jj][0])
		Tav = (gradient[jj-4][1]+gradient[jj-3][1]+gradient[jj-2][1]+gradient[jj-1][1]+gradient[jj][1]+gradient[jj+1][1]+gradient[jj+2][1]+gradient[jj+3][1]+gradient[jj+4][1])/float(9.0)
		temp9pt.append(Tav)
	dT9pt = np.column_stack((depth_9pt,temp9pt))


######################################################################################################
######################################################################################################

		# plotting temperature 9 point moving average
		plt.subplot(1,3,3)
		plt.plot(dT9pt[:,1], dT9pt[:,0])
		plt.ylabel("Depth [m]")
		plt.xlabel("T - 9pt moving av [degrees C]")
		plt.gca().invert_yaxis()
		plt.title("T 9pt MA")
		for i in range(0,len(flags.flag)):
			if (flags.flag[i] == "HB"):
				ref = flags.depth[i]
				plt.axhline(y=ref, hold=None, color='r')
			else:
 				continue

######################################################################################################
######################################################################################################
"""
This is a section of code that is removed from the hitbottom.py file after being used
"""

# statistics for the spikes (independent)
"""
collect information about what is the best standard deviation threshold (function input) 
to use to get the fewest false detections
The match precision for a "good match" will be set to +- 5m in depth (chosen arbitrarily)
Bad is if no match is detected or if difference is more than 5m
"""
"""
def get_stats(data, gradient, threshold):
	
	# determining the depth of actual HB flag:
	for i in range(0,len(flags.flag)):
			if (flags.flag[i] == "HB"):
				HB_depth = flags.depth[i]
			else:
 				continue
	
	# calculate peak accuracy (spike function returns int 0 if no peak found)
	spikes = spike(data, gradient, threshold) 
	if (type(spikes) != int):
		spike_depth = spikes[0][0]
		difference = HB_depth - spike_depth
	else:
		difference = 999
	
	# categorising (1=good, 0=bad)
	if (abs(difference) < 5):
		spike_stats = 1
	else:
		spike_stats = 0

	return(spike_stats)
"""

# importing libraries
import numpy as np

""" 
Opening file for reading the fraction of good detections found with each threshold
"""

threshold = []
good = []
bad = []

# reading file
with open("stats.txt") as f:
	next(f) #skipping header
	for line in f:
		line = line.split(',')
		threshold.append(int(line[0]))
		good.append(float(line[1]))
		bad.append(float(line[2].rstrip('\n')))

f.close()

threshold = np.array(threshold)
good = np.array(good)
bad = np.array(bad)

# calculation of the fraction of the hit bottoms that are being identified accurately with the
# temperature spikes method
frac = good/(good+bad)

# writing to terminal
for i in range(0,len(threshold)):
	print("for "+str(threshold[i])+"*sigma threshold, HB identified = "+str(frac[i])+" accuracy.")

"""
Find that the 3 sigma threshold for detections is the most ideal for determining HB
"""


######################################################################################################
######################################################################################################

	'''
	# true positive detection rate
	trueHB = 0
	detectHB = 0
	for i in range(0,m):
		# checking if the result is truly a "good" HB point
		if (ytest[i] == 1):
			trueHB += 1
			# counting the number of points identified out of these
			if (predictions[i][0] > predictions[i][1]):
				detectHB += 1
				identIndex[i] = 1
			else:
				continue
		else:
			continue
	true_pos = float(detectHB)/float(trueHB)
	if (vb == 1):
		print("True positive detection rate: "+str(true_pos))

	# true negative detection rate
	trueNeg = 0
	detectNeg = 0
	for i in range(0,m):
		# checking all points that are not hit bottoms
		if (ytest[i] == 0):
			trueNeg += 1
			# counting the number of these points that are identified as not HB points
			if (predictions[i][1] > predictions[i][0]):
				detectNeg += 1
				identIndex[i] = 1
			else:
				continue
		else:
			continue
	true_neg = float(detectNeg)/float(trueNeg)
	if (vb == 1):
		print("True negative detection rate: "+str(true_neg))

	# false positive rate
	falsePos = 0
	detectFalseNotHB = 0 
	for i in range(0,m):
		# checking all points that are not HB points
		if (ytest[i] == 0):
			falsePos += 1
			# taking fraction that are actually detected as HB
			if (predictions[i][0] > predictions[i][1]):
				detectFalseNotHB += 1
			else:
				continue
		else:
			continue
	false_pos = float(detectFalseNotHB)/float(falsePos)
	if (vb == 1):
		print("False positive detection rate: "+str(false_pos))

	# false negative rate
	falseNeg = 0
	detectFalseHB = 0
	for i in range(0,m):
		# checking all of the points that are HB
		if (ytest[i] == 1):
			falseNeg += 1
			# taking fraction of these that are not detected as HB
			if (predictions[i][1] > predictions[i][0]):
				detectFalseHB += 1
			else:
				continue
		else:
			continue
	false_neg = float(detectFalseHB)/float(falseNeg)
	if (vb == 1):
		print("False negative detection rate: "+str(false_neg))

	# pulling out other statistics: precision, recall and F-Score
	precision = float(detectHB)/float(detectHB+detectFalseNotHB)
	recall = float(detectHB)/float(detectFalseHB+detectHB)
	Fscore = float(2*precision*recall)/float(precision+recall)

	# printing to screen the results of this statistics call
	if (vb == 1):
		print("Precison: "+str(precision))
		print("Recall: "+str(recall))
	print("F-score: "+str(Fscore))	
	print("\n")
	'''

######################################################################################################
######################################################################################################

					"""
					# try all error points to find location	of hit bottom
					else:
						print("one entry is no good - going to look through error_pts")
						for j in range(0,len(error_pts)):
							print(str(j))
							below = 0
							# find correct index
							eq_index = 0
							for i in range(0,n):
								if (error_pts[j][0]==data[i][0]):
									eq_index = i
								else:
									continue
					
							for i in range(eq_index,n):
								for k in range(0,len(error_pts)):
									if (error_pts[k][0] == data[i][0]):
										below = below + 1
									else:
										continue
							frac_above = below/float(len(error_pts))
							if (frac_above > fracAbove):
								predHB = error_pts[j][0]
							else: 
								continue		
						"""
			"""
			# if there is no predicted hit bottom location (try each point of the error points)
			else:	
				print("no potential hb points found - looking through error points")
				for j in range(0,len(error_pts)):
					print(str(j))
					below = 0
					# find correct index
					eq_index = 0
					for i in range(0,n):
						if (error_pts[j][0]==data[i][0]):
							eq_index = i
						else:
							continue

					for i in range(eq_index,n):
						for k in range(0,len(error_pts)):
							if (error_pts[k][0] == data[i][0]):
								below = below + 1
							else:
								continue
					frac_above = below/float(len(error_pts))
					if (frac_above > fracAbove):
						predHB = error_pts[j][0]
					else: 
						continue
				"""

######################################################################################################
######################################################################################################

	# this code finds where the variation in the gradient is low
	# need to change it to find the regions where the deviation in the gradient is near the 
	# lowest values in the profile
	'''
	# calculation of the mean and standard deviation for that profile
	mean = sum(gradient[:,1])/float(len(gradient[:,1]))
	stddev = 0	
	n = len(gradient[:,1])
	for i in range(0,n):
		stddev = stddev + (float(gradient[i][1])-float(mean))**2
	stddev = stddev/float(len(gradient[:,1]))

	# finding the average gradient value over a specified range of points
	pointz = []
	pointT = []
	n = len(gradient[:,1])
	for i in range(0,n-consec):
		# mean and std.dev initialisation
		mean_sample = sum(gradient[i:i+consec,1])/float(consec)
		stddev_sample = 0
		# computing the standard dev for the next consec points
		for j in range(0,consec):
			stddev_sample = stddev_sample + (float(gradient[i+j][1])-float(mean_sample))**2
		stddev_sample = stddev_sample/float(consec)
		
		# saving the points that are within the threshold
		if (float(stddev_sample)/float(stddev) < threshold):
			pointz.append(float(data[i][0]))
			pointT.append(float(data[i][1]))
	
	# returning the points that fall within the threshold
	grad_drop = np.column_stack((pointz,pointT))
	return(grad_drop)
	'''

######################################################################################################
######################################################################################################



######################################################################################################
######################################################################################################



######################################################################################################
######################################################################################################



######################################################################################################
######################################################################################################



######################################################################################################
######################################################################################################



######################################################################################################
######################################################################################################



######################################################################################################
######################################################################################################



######################################################################################################
######################################################################################################



######################################################################################################
######################################################################################################


######################################################################################################
######################################################################################################



######################################################################################################
######################################################################################################