import math
from scipy import optimize
from scipy.stats.stats import pearsonr

def covar(xlist, ylist):
	try:
		sumlist = [xlist[i] + ylist[i] for i in range(0, max(len(xlist), len(ylist)))]
	except IndexError:
		print('Dimension mismatch')
		raise
	else:
		n = len(xlist)
	prodlist = []
	for i in range(0,n):
		prodlist.append(ylist[i] * xlist[i])
	return sum(prodlist)/n - (sum(xlist) * sum(ylist))/(n**2)

def var(xlist):
	return covar(xlist, xlist)

def correl(xlist, ylist):
	return covar(xlist, ylist)/math.sqrt(var(xlist)*var(ylist))

def mean(list):
	return sum(list)/len(list)
	
def linregbx(xlist, ylist, xvar):
	b = covar(xlist, ylist)/var(xlist)
	a = mean(ylist) - b*mean(xlist)
	return b*xvar + a
	
def residsq(xlist, ylist, f):
	m = []
	for item in resid(xlist, ylist, f):
		m.append(item**2)
	return sum(m)

def partdiff(f, xvect, h = 1e-8):
	try:
		blow = f(xvect)
	except IndexError:
		print('Not enough arguments')
		raise
	except TypeError:
		print('Argument must be a vector')
		raise
	diff = []
	for i in range(0, len(xvect)):
		holding = xvect[0:i]
		holding.append(xvect[i] + h)
		for k in range(i+1, len(xvect)):
			holding.append(xvect[k])
		diff.append((f(holding) - f(xvect))/h)
	return diff
	
def multireg(xmat, ylist, verbose = False):
	if verbose:
		print('Testing dimensions...')
	for i in range(1, len(xmat)):
		try:
			blow = xmat[0] + xmat[i]
		except TypeError:
			print('Not a valid matrix')
			raise

	if verbose:
		print('Dimension test okay!')
	currylist = ylist
	coeff = [0,0]
	for index in range(0, len(xmat[0])):
		def line(brack):
			temp = 0
			for i in range(1, len(brack)):
				temp = temp + sum([brack[i]*xmat[k][i-1] + brack[0] - currylist[i-1] for k in range(0, len(ylist))])**2
			return math.sqrt(temp)
		least = optimize.linearmixing(line, coeff).x.tolist()
		coeff = least
		if verbose:
			print('Current coefficient list is ', coeff)
		def f(x):
			temp = coeff[0]
			for i in range(1, len(coeff)):
				temp = temp + coeff[i]*x[i-1]
			return temp
		reglist = [f(xmat[k]) for k in range(0, len(ylist))]
		if verbose:
			print('Current x list is ', currxlist)
		currylist = [currylist[i] - reglist[i] for i in range(0, len(ylist))]
		if verbose:
			print('Current residuals list is ', currylist)
		coeff.append(0)
	print('The furst entry in the output is the constant term.')
	return coeff
	
def multilinreg(xmat, ylist, verbose = False):
	if verbose:
		print('Testing dimensions...')
	for i in range(1, len(xmat)):
		try:
			blow = xmat[0] + xmat[i]
		except TypeError:
			print('Not a valid matrix')
			raise
	for k in range(0, max(len(xmat), len(ylist))):
		try:
			blow = xmat[k][0] + ylist[k]
		except IndexError:
			print('Dimension mismatch')
			raise
	if verbose:
		print('Dimension test okay!')
	def resid(xmat, ylist):
		return lambda lincoeffs: [sum([lincoeffs[i]*xmat[k][i-1] for i in range(1, len(xmat[0])+1)]) + lincoeffs[0] - ylist[k] for k in range(0, len(ylist))] # Residuals
	f = resid(xmat, ylist)
	def sumresid(xmat, ylist):
		return lambda lincoeffs: sum(f(lincoeffs))
	g = sumresid(xmat, ylist)
	def residsq(xmat, ylist):
		return lambda lincoeffs: sum([f(lincoeffs)[i]**2 for i in range(0, len(lincoeffs) - 1)])
	h = residsq(xmat, ylist)
	def mean(x):
		return sum(x)/len(x)
	coeff = [mean(ylist)]
	for i in range(0, len(xmat[0])):
		coeff.append(0)
	if verbose:
		print('Initial residuals: ', f(coeff))
		print('Initial residual square sum: ', h(coeff))
		print('Initial sum of residuals: ', g(coeff))
	least = optimize.linearmixing(h, coeff)
	# Choose the coefficients that minimize the sum of the residuals squared subject to keeping the sum of the residuals equal to 0.
	if verbose:
		print(least)
	coeff = least.x.tolist()
	if verbose:
		print('Final residuals: ', f(coeff))
		print('Final residual square sum: ', h(coeff))
		print('Final sum of residuals: ', g(coeff))
	if verbose:
		print('Coefficients: ', coeff)
	def out(x):
		return coeff[0] + sum([coeff[i]*x[i-1] for i in range(1, len(x) + 1)])
	return out
	
def factorial(x):
	try:
		brony = range(0,x)
	except:
		print('Arg must be type int')
		raise
	if x < 0:
		raise ValueError('Arg must be nonnegative')
	return int(math.gamma(x+1))

def factorial2(x):
	try:
		brony = range(0,x)
	except:
		print('Arg must be type int')
		raise
	if x < 0:
		raise ValueError('Arg must be nonnegative')
	elif x == 0:
		return 1
	else:
		output = 1
		for i in range(1,x+1):
			output = output * i
		return output

def corrmat(xmat, verbose = False):
	try:
		verbose and not verbose
	except:
		print('Second positional argument must be Boolean.')
		raise
	print('Correlation matrix of row vectors')
	result = []
	n = len(xmat)
	for i in range(0,n):
		if verbose:
			print('Iteration', i)
		templist = []
		for j in range(0,n):
			templist.append(correl(xmat[i], xmat[j]))
		result.append(templist)
		if verbose:
			print('Correlation row stored.')
	return result