'''
'''

from reading_lib import *
from commands import getoutput

###########################################################
#                                                         #
#                    Basic Math Stuff                     #
#                                                         #
###########################################################

def pythagorean(lat1, lon1, lat2, lon2):
	# Finds the distance between two points
	
	a = lon1-lon2
	b = lat1-lat2
	
	return sqrt(a**2 + b**2)


def moving_avg(list, points):
	
	if points%2 == 0:
		padding = (points)/2
	else:
		padding = (points - 1)/2
	
	moving_avgs = []
	for i in range(padding):
		moving_avgs.append(list[i])
	
	for i in range(padding, len(list)-padding):
		
		sample = list[i-padding:i+padding]
		moving_avgs.append(average(sample))
		
	for i in range(padding):
		moving_avgs.append(float('nan'))
			
	return moving_avgs
	

def list_avgs(lst):
	
	for i in range(len(lst)):
		lst[i] = asarray(lst[i])
	
	sums = lst[0]
	for i in range(len(lst)-1):
		sums = sums + lst[i+1]
		
	return sums / float(len(lst))
	
	
def sum_squared(set):
	
	sum = 0
	
	for x in set:
		sum+=x**2
		
	return sum	
	
	
def sum_product(set1, set2):
	
	sum = 0
	
	for i in range(len(set1)):
		sum+=(set1[i]*set2[i])
		
	return sum
	

def t_test(sample1, sample2):
	
	N1 = float(len(sample1))
	N2 = float(len(sample2))
	
	M1 = mean(sample1)
	S1 = std(sample1)
	M2 = mean(sample2)
	S2 = std(sample2)
	
	return (M1-M2)/sqrt((S1**2/N1)+(S2**2/N2))
	

def R_squared(xs, ys):

	N = len(xs)
	prods = []
	sqrd_x, sqrd_y = [], []
	for i in range(N):
		prods.append(xs[i]*ys[i])
		sqrd_x.append(xs[i]**2)
		sqrd_y.append(ys[i]**2)

	r = (N*sum(prods)-sum(ys)*sum(xs)) /\
	sqrt((N*sum(sqrd_x)-sum(xs)**2)*(N*sum(sqrd_y)-sum(ys)**2))
	
	return r**2
	

def best_fit_line(Xs, Ys):
	
	denominator = sum_squared(Xs) - mean(Xs) * sum(Xs)
	
	m = (sum_product(Xs,Ys) - mean(Ys) * sum(Xs)) / denominator
	b = (mean(Ys) * sum_squared(Xs) -  mean(Xs) * sum_product(Xs,Ys)) / denominator
	
	return m, b


###########################################################
#                                                         #
#                        Booleans                         #
#                                                         #
###########################################################

def is_DPZ(loc, SL1, SL2):
	# Is this location the the depopulated zone?
	
	x = loc[0]
	y = loc[1]
	
	if SL1[int(x)] > -1 and SL2[int(x)] > -1:
		if y >= SL1[int(x)] and y < SL2[int(x)]:
			DPZ = True
		else:
			DPZ = False
	else:
		DPZ = False
		
	return DPZ
	
def is_overSnow(loc, SL2):
	# Is this location over the adjusted snow extent?
	
	x = loc[0]
	y = loc[1]
	
	if SL2[int(x)] > -1:
		if y > SL2[int(x)]:
			over = True
		else:
			over = False
	else:
		over = False
		
	return over
	
def is_clipper(case_no):
	
	clipper_list = [True, False, False, True, False, True, True, False,
	True, False, False, False, True, True, False, False, False, True, 
	True, False]
	
	return clipper_list[case_no]


###########################################################
#                                                         #
#                   Cyclone Functions                     #
#                                                         #
###########################################################

def MTD(path1, path2):
	# Mean trajectory deviation (per interval)
	
	if path1 < path2:
		path2 = path2[:len(path1)+1]
	elif path1 > path2:
		path1 = path1[:len(path2)+1]
	
	dists = []
	for i in range(len(path1)):
		loc1 = path1[i]
		loc2 = path2[i]
		
		dist = pythagorean(loc1[0], loc1[1], loc2[0], loc2[1])
		dists.append(dist*30)
	
	return sum(dists)/(len(dists))

def mean_poleward(path1, path2):
	# Mean poleward deviation of cyclone trajectory
	
	y1s, y2s = [], []
	
	for i in range(len(path1)):
		
		y1s.append(path1[i][1])
		y2s.append(path2[i][1])
		
	return (average(y2s) - average(y1s))*30
	
	
def mean_snowward(path1, path2, SL):
	# Mean snowward deviation of cyclone trajectory
	
	ys = []
	
	for i in range(len(path1)):
		
		x = path2[i][0]
		if x > 200:
			break
		snowline = SL[int(x)]
		if snowline > -1:
			pass
		else:
			break
		
		ol_dist = abs(snowline - path1[i][1])
		nw_dist = abs(snowline - path2[i][1])
		ys.append(ol_dist - nw_dist)
		
	return average(ys)*30