def compute_stats(te_diffs, gene_diffs, plot_dir):
    pvals = []
    table_lines = []

    for te_or in te_diffs:
        rep, fam, orient = te_or
        
        for sample_key in te_diffs[te_or]:        
            sample1, sample2 = sample_key

            # if enough data
            if len(te_diffs[te_or][sample_key]) >= 10:
                wo_te = list((gene_diffs[sample_key] - te_diffs[te_or][sample_key]).elements())
                w_te = list(te_diffs[te_or][sample_key].elements())

                wo_mean = stats.mean(wo_te)
                w_mean = stats.mean(w_te)

                z, p = stats.mannwhitneyu(w_te, wo_te)

                cols = (rep, fam, orient, sample1, sample2, len(w_te), w_mean, wo_mean, z, p)
                table_lines.append('%-17s %-17s  %1s  %-10s %-10s %6d %9.2f %9.2f %8.2f %10.2e' % cols)

                pvals.append(p)

                # plot ...
                if rep in ['*'] and fam in ['*','LINE/L1','SINE/Alu','LTR/ERV1','LTR/ERVL-MaLR','LINE/L2','LTR/ERVL','SINE/MIR','DNA/hAT-Charlie','LTR/ERVK','DNA/TcMar-Tigger']:
                    out_pdf = '%s/%s_%s_%s_%s-%s.pdf' % (plot_dir,rep.replace('/','-'),fam.replace('/','-'),orient,sample1,sample2)
                    cdf_plot(te_or, w_te, wo_te, out_pdf)

    return table_lines, pvals

 def testMean(self):
     """
     Check that mean works as expected.
     """
     self.assertAlmostEqual(stats.mean(self.dataA), self.meanA, 5)
     self.assertAlmostEqual(stats.mean(self.dataB), self.meanB, 5)
     return

def least_squares_fit(x, y):
    """given training values for x and y,
    find the least-squares values of alpha and beta"""
    beta = stats.correlation(x, y) * \
        stats.standard_deviation(y) / stats.standard_deviation(x)
    alpha = stats.mean(y) - beta * stats.mean(x)
    return alpha, beta

def sync_check():
#    print 'Checking sync...'
    max_mcnt_difference=4
    mcnts=dict()
    mcnts_list=[]
    mcnt_tot=0

    for f,fpga in enumerate(fpgas):
        mcnts[f]=dict()
        try:
            hdr_index=bram_oob[f]['hdr'].index(1)
        except:
            print 'ERR: No headers found in BRAM. Are the F engines properly connected?'
            exit()

        pkt_64bit = struct.unpack('>Q',bram_dmp['bram_msb'][f]['data'][(4*hdr_index):(4*hdr_index)+4]+bram_dmp['bram_lsb'][f]['data'][(4*hdr_index):(4*hdr_index)+4])[0]
        mcnts[f]['mcnt'] =(pkt_64bit&((2**64)-(2**16)))>>16
        mcnts_list.append(mcnts[f]['mcnt'])
#        print '[%s] MCNT: %i'%(servers[f],mcnts[f]['mcnt'])

    mcnts['mean']=stats.mean(mcnts_list)
    mcnts['median']=stats.median(mcnts_list)
    mcnts['mode']=stats.mode(mcnts_list)
    mcnts['modalmean']=stats.mean(mcnts['mode'][1])

#    print 'mean: %i, median: %i, modal mean: %i mode:'%(mcnts['mean'],mcnts['median'],mcnts['modalmean']),mcnts['mode']
    
    for f,fpga in enumerate(fpgas):
        if mcnts[f]['mcnt']>(mcnts['modalmean']+max_mcnt_difference) or mcnts[f]['mcnt'] < (mcnts['modalmean']-max_mcnt_difference):
            print '%s OUT OF SYNC!!'%servers[f]
            mcnts[f]['sync_status']='FAIL with error of %i'%(mcnts[f]['mcnt']-mcnts['modalmean'])
        else:
            mcnts[f]['sync_status']='PASS'

    return mcnts

  def statsex(self, objects):

    """
	Do some statistics on a source list
	Return dictionary
    """

    import stats, pstat
    
    # Return if we have no objects
    if len(objects) == 0:
      return 0	 

    # Define dictionary to hold statistics	
    stat = {}

    # Get number of objects
    stat['N'] = str(len(objects))

    # Define list (float) of FWHM values
    fwhm = [ float(obj[7]) for obj in objects ]
 
    # Define list (float) of ELLIPTICITY values
    el = [ float(obj[6]) for obj in objects ]

    # Define list (float) of THETA_IMAGE values
    pa = [ float(obj[5]) for obj in objects ]

    # Define list (float) of 'Stella-like' values
    stella = [ float(obj[9]) for obj in objects ]	

    # Create a histogram of FWHM values of binsize 1 pixel
    hfwhm = stats.histogram(fwhm,40,[0,40])[0]
    
    stat['medianFWHM'] = "%.2f" % stats.median(fwhm)
    stat['meanFWHM']   = "%.2f" % stats.mean(fwhm)
    stat['modeFWHM']   = "%.2f" % float(hfwhm.index(max(hfwhm))+0.5)

    try:	
       stat['stdevFWHM']  = "%.2f" % stats.stdev(fwhm)
    except ZeroDivisionError:
       stat['stdevFWHM'] = '0.00';

    stat['medianEL'] = "%.2f" % stats.median(el)
    stat['meanEL']   = "%.2f" % stats.mean(el)

    try:
      stat['stdevEL']  = "%.2f" % stats.stdev(el)
    except ZeroDivisionError:
      stat['stdevEL']  = '0.00' 

    # Histogram of Ellipticity PA (-180 to 180, bins of 45 deg)
    #stat['histoTHETA'] = stats.histogram(pa,8,[-180,180])[0]

    # Histogram of Stellarity (0 to 1, bins of 0.05)
    #stat['histoStella']  = stats.histogram(stella,20,[0,1.01])[0]   

    return stat

 def testOnTuples(self):
     """
     Checks that methods also work on tuples.
     """
     self.assertAlmostEqual(stats.mean(tuple(self.dataA)), self.meanA, 5)
     self.assertAlmostEqual(stats.mean(tuple(self.dataB)), self.meanB, 5)
     self.assertAlmostEqual(stats.stddev(tuple(self.dataA)), self.stddevA, 5)
     self.assertAlmostEqual(stats.stddev(tuple(self.dataB)), self.stddevB, 5)
     return

def corr(xdata, ydata):
    """corr(xydata) -> float
    corr(xdata, ydata) -> float

    Return the sample Pearson's Correlation Coefficient of (x,y) data.

    If ydata is None or not given, then xdata must be an iterable of (x, y)
    pairs. Otherwise, both xdata and ydata must be iterables of values, which
    will be truncated to the shorter of the two.

    >>> corr([(0.1, 2.3), (0.5, 2.7), (1.2, 3.1), (1.7, 2.9)])
    ... #doctest: +ELLIPSIS
    0.827429009335...

    The Pearson correlation is +1 in the case of a perfect positive
    correlation (i.e. an increasing linear relationship), -1 in the case of
    a perfect anti-correlation (i.e. a decreasing linear relationship), and
    some value between -1 and 1 in all other cases, indicating the degree
    of linear dependence between the variables.

    >>> xdata = [1, 2, 3, 4, 5, 6]
    >>> ydata = [2*x for x in xdata]  # Perfect correlation.
    >>> corr(xdata, ydata)
    1.0
    >>> corr(xdata, [5-y for y in ydata])  # Perfect anti-correlation.
    -1.0

    If there are not at least two data points, or if either all the x values
    or all the y values are equal, StatsError is raised.
    """
    n = len(xdata)
    assert n == len(ydata)
    if n < 2:
        raise StatsError(
            'correlation requires at least two data points, got %d' % n)
    # First pass is to determine the means.
    mx = stats.mean(xdata)
    my = stats.mean(ydata)
    # Second pass to determine the standard deviations.
    sx = stats.stdev(xdata, mx)
    sy = stats.stdev(ydata, my)
    if sx == 0:
        raise StatsError('all x values are equal')
    if sy == 0:
        raise StatsError('all y values are equal')
    # Third pass to calculate the correlation coefficient.
    ap = add_partial
    total = []
    for x, y in zip(xdata, ydata):
        term = ((x-mx)/sx) * ((y-my)/sy)
        ap(term, total)
    r = math.fsum(total)/(n-1)
    assert -1 <= r <= r
    return r

 def get_modules(self, cutoff=.05):
     modules = []
     for e in self:
         if e.val < min(e.lo_min, e.hi_min, cutoff):
             if self.datatype=="continuous":
                 e.desc = "lo" if mean(e.a) < mean(e.b) else "hi"
             else:
                 e.desc = "enriched"
             modules.append(e)
         else:
             modules += e.get_modules(cutoff=cutoff)
     return modules

 def check_basic(self):
     a = [3,4,5,10,-3,-5,6]
     af = [3.,4,5,10,-3,-5,-6]
     Na = len(a)
     Naf = len(af)
     mn1 = 0.0
     for el in a:
         mn1 += el / float(Na)
     assert_almost_equal(stats.mean(a),mn1,11)
     mn2 = 0.0
     for el in af:
         mn2 += el / float(Naf)
     assert_almost_equal(stats.mean(af),mn2,11)

def _SP(xdata, mx, ydata, my):
    """SP = sum of product of deviations.
    Helper function for calculating covariance directly.
    """
    if mx is None:
        # Two pass algorithm.
        xdata = as_sequence(xdata)
        mx = stats.mean(xdata)
    if my is None:
        # Two pass algorithm.
        ydata = as_sequence(ydata)
        my = stats.mean(ydata)
    return _generalised_sum(zip(xdata, ydata), lambda t: (t[0]-mx)*(t[1]-my))

 def check_2d(self):
     a = [[1.0, 2.0, 3.0],
          [2.0, 4.0, 6.0],
          [8.0, 12.0, 7.0]]
     A = array(a,'d')
     N1,N2 = (3,3)
     mn1 = zeros(N2,'d')
     for k in range(N1):
         mn1 += A[k,:] / N1
     allclose(stats.mean(a),mn1,rtol=1e-13,atol=1e-13)
     mn2 = zeros(N1,'d')
     for k in range(N2):
         mn2 += A[:,k] / N2
     allclose(stats.mean(a,axis=0),mn2,rtol=1e-13,atol=1e-13)

def calculateDividingLine(gestures, maybeGestures, nonGestures):
	numFolds = min(TESTING_FOLDS, len(gestures))
	
	allGestureDistances = []
	allNonGestureDistances = []
	
	for foldNum in range(numFolds):	
		trainingGestures = [gesture for i, gesture in enumerate(gestures) if i % numFolds != foldNum]
		testingGestures = [localTimeGestures for i, localTimeGestures in enumerate(maybeGestures) if i % numFolds == foldNum]
		
		#print 'train, test #s: ', len(trainingGestures), len(testingGestures)
		
		#make a distance calculator based on the subset of hte training data
		distanceCalculator = gestureDistanceCalculator.GestureDistanceCalculator(trainingGestures)
		
		#each localTimeGestures is a list of the closest times to when a gesture was identified in training
		#since the output can be triggered at slightly different times, we should look for a minimum near where
		#the gesture is known to have happened, compared to the training gestures
		gestureDistances = []
		#print testingGestures
		for localTimeGestureSet in testingGestures:
			
			closestDistance = min(map(distanceCalculator.getDistance, localTimeGestureSet))
			gestureDistances.append(closestDistance)
		
		#gestureDistances = map(distanceCalculator.getDistance, testingGestures)
		#print gestureDistances
		nonGestureDistances = map(distanceCalculator.getDistance, nonGestures)
		#print gestureDistances
		
		allGestureDistances += gestureDistances
		allNonGestureDistances += nonGestureDistances
		#break
		
	#print len(allGestureDistances), len(allNonGestureDistances)
	print 'means: ', stats.mean(allGestureDistances), stats.mean(allNonGestureDistances)
	
	print 'std devs: ', stats.stdDev(allGestureDistances), stats.stdDev(allNonGestureDistances)
	
	meanGesture = stats.mean(allGestureDistances)
	meanNon = stats.mean(allNonGestureDistances)
	
	devGesture = stats.stdDev(allGestureDistances)
	devNon = stats.stdDev(allNonGestureDistances)
	
	line = (meanGesture * devNon + meanNon * devGesture) / ( devGesture + devNon)
	
	#print line
	return line

def diff_fpkm(diff_file, pseudocount):
    gene_fpkms = {}

    diff_in = open(diff_file)
    diff_in.readline()
    for line in diff_in:
        a = line.split('\t')

        gene_id = a[0]
        sample1 = a[4]
        sample2 = a[5]
        status = a[6]
        fpkm1 = float(a[7])
        fpkm2 = float(a[8])

        if status == 'OK':
            if gene_id in gene_fpkms:
                gene_fpkms[gene_id] += [fpkm1, fpkm2]
            else:
                gene_fpkms[gene_id] = [fpkm1, fpkm2]

    diff_in.close()

    gene_fpkm = {}
    for gene_id in gene_fpkms:
        log_fpkms = [math.log(fpkm+pseudocount,2) for fpkm in gene_fpkms[gene_id]]
        gene_fpkm[gene_id] = stats.mean(log_fpkms)

    return gene_fpkm

def wordSummary(db, table):
	f = open("wordSummary_%s.txt" % table, 'w')
	d = {}
	header = "word, length, rtAVG, rtSTD, total, percCorrect\n"
	f.write(header)
	wordList = []
	sql = "SELECT DISTINCT(word) FROM %s" % table
	for w in db.query(sql):
		wordList.append(w[0])
		
	for word in wordList:
		sql = "SELECT RT FROM %s WHERE word = '%s' AND incorrect = 0" % (table, word)
		wordLen = len(word)
		rtList = []
		zList = []

		for rt in db.query(sql):
			rtList.append(rt[0])

		rtAVG = stats.mean(rtList)
		rtSTD = stats.samplestdev(rtList)


		total = db.query("SELECT COUNT(*) FROM %s WHERE word = '%s'" % (table, word))[0][0]
		percCorrect = float(len(rtList)) / float(total) * 100.0

		print len(rtList), total

		myString = "%s, %i, %f, %f, %i, %f\n" % (word, wordLen, rtAVG, rtSTD, total, percCorrect)
		print myString
		f.write(myString)


	f.close()

def main():
    [(stat, first), (stat, second)] = load_stats(sys.argv[1:])

    # Attempt to increase robustness by dropping the outlying 10% of values.
    first = trim(first, 0.1)
    second = trim(second, 0.1)

    fmean = stats.mean(first)
    smean = stats.mean(second)
    p = 1 - ttest_1samp(second, fmean)[1][0]
    if p >= 0.95:
        # rejected the null hypothesis
        print sys.argv[1], 'mean of', fmean, 'differs from', sys.argv[2], 'mean of', smean, '(%2.0f%%)' % (p * 100,)
    else:
        # failed to reject the null hypothesis
        print 'cannot prove means (%s, %s) differ (%2.0f%%)' % (fmean, smean, p * 100,)

def scaleTestMinFinding():
	xs = range(10)
	distances = []
	noise = 3.5
	n = 1000000
	for i in range(n):
		a = random()
		b = random()
		c = random()
		ys = [x*x*a + x*b + c + random() * noise for x in xs]
		
		#print a, b, c, polynomialFit(xs, ys)[::-1]
		minExp, unc = polynomialFindMinimum(xs, ys, returnErrors = True)
		minCalc = -b/(2.0*a)
		dist = (minCalc - minExp) / unc
		#print minCalc, minExp, unc, dist
		distances.append(dist)
		
	print 'mean: %f' % stats.mean(distances)
	print 'stdDev: %f' % stats.stdDev(distances)
	for sigma in [1, 2, 3]:
		print 'With %d sigma: %f%%' % (sigma, 100.0 * sum([int(abs(d) < sigma) for d in distances]) / n)
	
	pylab.hist(distances, bins = 50, range = (-5, 5))
	pylab.show()

def scale(data_matrix):
    num_rows, num_cols = shape(data_matrix)
    means = [mean(get_column(data_matrix,j))
             for j in range(num_cols)]
    stdevs = [standard_deviation(get_column(data_matrix,j))
              for j in range(num_cols)]
    return means, stdevs

def linear_regression(x, y):
  """
  See: https://www.khanacademy.org/math/probability/regression/regression-correlation/v/regression-line-example
  """
  xy_mean = s.xy_mean(x, y)
  print xy_mean
  x_mean = s.mean(x)
  y_mean = s.mean(y)
  x_squared_mean = s.mean([xi ** 2 for xi in x])

  # Slope.
  a = (x_mean * y_mean - xy_mean) / (x_mean ** 2 - x_squared_mean)
  # Intercept.
  b = y_mean - a * x_mean

  return (a, b)

 def testVariance(self):
     data = [1, 2, 3]
     assert stats.mean(data) == 2
     self.assertEqual(stats.pvariance(data), 2/3)
     self.assertEqual(stats.variance(data), 1.0)
     self.assertEqual(stats.pstdev(data), math.sqrt(2/3))
     self.assertEqual(stats.stdev(data), 1.0)

def test_mean1():
    obs = mean([0, 0, 0, 0])
    exp = 0
    assert_equal(obs, exp)

    obs = mean([0, 200])
    exp = 100
    assert_equal(obs, exp)

    obs = mean([0, -200])
    exp = -100
    assert_equal(obs, exp)

    obs = mean([0]) 
    exp = 0
    assert_equal(obs, exp)