import numpy as np
from scipy.stats.stats import pearsonr
from scipy.stats import norm
def thresholdcompare(countfs='/home/scottn/sixmercounts.csv', phenofs='/home/scottn/runwithgemma/mouse.pheno', headerfs='/home/scottn/sixmerlist.csv', threshold=None, pctile=.1, includethresholds=True):
	## Syntax: thresholdcompare(countfs, phenofs, headerfs, threshold, pctile)
	## Output: List, with each row containing the header, correlation coefficient, and p-value.
	## All fields optional.  Default values:
	##     countfs = '/home/scottn/sixmercounts.csv' (str)
	##         This file provides the x values for the analysis.
	##     phenofs = '/home/scottn/runwithgemma/mouse.pheno' (str)
	##         This file provides the y values for the analysis.
	##     headerfs = '/home/scottn/sixmerlist.csv' (str)
	##         This file provides the headers for the output.
	##     threshold = None (float)
	##         When a value is specified, this program will automatically filter anything less than that value.  When no value is specified, this program will compute the threshold for each row based on the value of pctile.
	##     pctile = 0.1 (float, strictly between 0 and 1)
	##         The program will use this percentile to generate the threshold for each row assuming the counts follow a normal distribution.
	##     includethresholds = True (bol)
	##         If True, this program will return two outputs: the correlations and the thresholds.  If False, the program will only return the standard output.
	if not(0 < pctile < 1):
		raise ValueError('Percentile must be strictly between 0 and 1.')
	print('Reading files...')
	countfile = open(countfs)
	phenofile = open(phenofs)
	counts = np.genfromtxt(countfile, delimiter=None, dtype=np.float64, names=None)
	phenofail = np.genfromtxt(phenofile, delimiter=',', dtype=None, names=None)
	phenos = []
	headerfile = open(headerfs)
	headers = np.genfromtxt(headerfile, delimiter='\n', dtype='|S6', names=None)
	for row in phenofail:
		phenos.append(row[1])
	print('File read complete!')
	output = []
	n = len(headers)
#	n = 50
	threshes = []
	for j in range(0,n):
		trunk = []
		skunk = []
		if threshold == None:
			mu = np.mean(counts[:,j], dtype=np.float64)
			sx = np.std(counts[:,j], dtype=np.float64)
			threshold = norm.ppf(pctile, mu, sx)
		threshes.append(threshold)
		for i in range(0, len(counts)):
			if counts[i,j] > threshold:
				trunk.append(counts[i,j])
				skunk.append(phenos[i])
		#if len(trunk) > 2:
		(c, p) = pearsonr(trunk, skunk)
		#else:
	#	(c, p) = (np.nan, np.nan)
		output.append([headers[j],c,p])
		try:
			mu + 0
		except:
			pass
		else:
			threshold = None
	if includethresholds:
		return output, threshes
	else:
		return output
def binarycompare(binfs='/home/scottn/sixmersbin.csv', phenofs='/home/scottn/runwithgemma/mouse.pheno', headerfs='/home/scottn/sixmerlist.csv'):
	## Syntax: thresholdcompare(countfs, phenofs, headerfs)
	## Output: List, with each row containing the header, mean of motif absent, mean of motif present, difference, and p-value.
	## All fields optional.  Default values:
	##     countfs = '/home/scottn/sixmercounts.csv' (str)
	##         This file provides the x values for the analysis.
	##     phenofs = '/home/scottn/runwithgemma/mouse.pheno' (str)
	##         This file provides the y values for the analysis.
	##     headerfs = '/home/scottn/sixmerlist.csv' (str)
	##         This file provides the headers for the output.
	print('Reading files...')
	binfile = open(binfs)
	phenofile = open(phenofs)
	bins = np.genfromtxt(binfile, delimiter=',', dtype=np.float64, names=None)
	phenofail = np.genfromtxt(phenofile, delimiter=',', dtype=None, names=None)
	phenos = []
	headerfile = open(headerfs)
	headers = np.genfromtxt(headerfile, delimiter='\n', dtype='|S6', names=None)
	for row in phenofail:
		phenos.append(row[1])
	output = []
	print('File read complete!')
	n = len(headers)
#	n = 50
	for j in range(0,n):
		trunk = []
		skunk = []
		for i in range(0, len(bins)):
			if bins[i,j] == 0:
				trunk.append(phenos[i])
			else:
				skunk.append(phenos[i])
		if trunk == []:
			mu1 = 0
		else:
			mu1 = np.mean(trunk, dtype=np.float64)
		if skunk == []:
			mu2 = 0
		else:
			mu2 = np.mean(skunk, dtype=np.float64)
		output.append([headers[j], mu1, mu2, mu1-mu2])
	diffs = []
	for row in output:
		diffs.append(row[3])
	# mu = np.mean(diffs, dtype=np.float64)
	sx = np.std(diffs, dtype=np.float64)
	for row in output:
		if row[3] > 0:
			row.append(norm.cdf(-row[3], 0, sx))
		else:
			row.append(norm.cdf(row[3], 0, sx))
	return output