def computeShannon(fh, header, verbose):
    for line in fh:
        if header:
            header -= 1
            print "%s\tShannon_index" % line.rstrip()
            continue
        #-------------------------
        lineL = line.split()
        expr = [float(i) for i in lineL[1:]]
        if verbose:
            print >>sys.stderr, expr
        expr_sum = sum(expr)
        if verbose:
            print >>sys.stderr, expr_sum
        assert expr_sum != 0
        expr_R = [1.0 * i / expr_sum for i in expr]
        if verbose:
            print >>sys.stderr, expr_R
        expr_Log = []
        for i in expr_R:
            if i != 0:
                expr_Log.append(i*ln(i)/ln(2))
            else:
                expr_Log.append(i)
        if verbose:
            print >>sys.stderr, expr_Log
        shannon = -1 * sum(expr_Log)
        print "%s\t%s" % (line.strip(),str(shannon))

def test_convergenceRate():

        problem = trigonometricSolution()
        h = 0.1
        T = 1.0
        hValues = []
        errorValues = []
        for i in range(7):
                dt = h
                xnodes = ynodes = int(1/sqrt(dt))
                u_e, u = nld.runSolver(problem,T, dt, [xnodes,ynodes])
                e = u_e.vector().array() - u.vector().array()
                E = np.sqrt(np.sum(e**2)/u.vector().array().size)
                hValues.append(h)
                errorValues.append(E)
                h/=2
        
        
        r = zeros((len(hValues)))
        for i in range(1, len(hValues)):    
            r[i] = ln(errorValues[i-1]/errorValues[i])/ ln(hValues[i-1]/hValues[i])
            print "h = ",hValues[i-1], " E = ", errorValues[i-1], " r = ", r[i]
        
        if not nt.assert_almost_equal(1,r[-1],places=1):
            print "test_convergenceRate succeeded!"

 def propagate(self, context):
     l = self.lhs.getValue()
     b = self.base.getValue()
     e = self.exponent.getValue()
     if l:
         if b:
             if e: # l,b,e
                 if l==b**e:
                     pass # but overconstrained
                 else:
                     print("Error! " + self.name + " is overconstrained")
                     return False
             else: # l,b => e
                 self.exponent.setValue(math.ln(l)/math.ln(b))
         else:
             if e: #l,e => b
                 self.base.setValue(l**(1/e)) # TODO allow multiple solutions...
             else: # l only
                 pass
     else: # not l
         if b:
             if e: # b,e => l
                 self.lhs.setValue(b**e)
             else: # b only
                 pass
         else: # no l, no b
             pass
     return True

def AIC(records, SSR, p):
    """ln(SSR(p)/T + (p+1)2/T is the formula. The lower this is, the
        better the fit. This is the Akaike Information Criterion.
        I am using it because I don't think that any of the models
        accurately represents reality, so I am finding the one that
        fits best. """
    return math.ln(SSR(records, p) / math.ln(len(records))) + (p + 1) * 2 / len(records)

def main(fname = '/home/rerla/Downloads/Galaxy28-[HT100_K562_differential_results.xls'):
    """ two passes so we can estimate percentiles and estimate a kind of fisher's independent combined
    p value
    """
    dat = open(fname,'r').readlines()
    dhead = dat[0].strip().split('\t')
    dat = dat[1:]
    dat = [x.split() for x in dat if len(x.split()) > 0]
    conts=[]
    treats=[]
    for index,row in enumerate(dat): # header
       cont = row[5]
       treat = row[6]
       conts.append((float(cont),index))
       treats.append((float(treat),index)) # so can recover row
    conts.sort()
    conts.reverse()
    treats.sort()
    treats.reverse()
    treats = [(rank,x[0],x[1]) for rank,x in enumerate(treats)]
    conts = [(rank,x[0],x[1]) for rank,x in enumerate(conts)]
    tdict = dict(zip([x[2] for x in treats],treats)) # decorate so can lookup a data row
    cdict = dict(zip([x[2] for x in conts],conts)) # decorate so can lookup a data row
    res = []
    n = float(len(dat) - 1)
    for dindex in range(len(dat)): # dindex = record in d
       if dindex % 10000 == 0:
           print dindex
       treati,treat,tindex = tdict.get(dindex,(0,0,0))
       conti,cont,cindex = cdict.get(dindex,(0,0,0))
       crank = conti/n
       trank = treati/n
       try:
           logfold = math.log(treat/cont,2)
       except:
           print 'bad logfold treat=%f cont=%f' % (treat,cont)
           logfold = 0

       logA = math.log(treat+cont,2)
       try:
           logM = math.log(abs(treat-cont),2)
       except:
           print 'bad logM treat=%f cont=%f' % (treat,cont)
           logM = 0
       try:
           fish = -2.0*(math.ln(crank) + math.ln(trank))
       except:
           print "bad fisher's combined crank=%f trank=%f" % (crank,trank)
           fish = 0
       row = copy.copy(dat[dindex])
       row += ['%i' % conti, '%i' % treati,'%f' % logfold,'%f' % logA,
        '%f' % logM,'%f' % crank,'%f' % trank,'%f' % fish]
       res.append('\t'.join(row))
    h = copy.copy(dhead)
    h += ['conti','treati','logfold','logA','logM','crank','trank','fakefishers']
    res.insert(0,h)
    outfname = '%s_fixed.xls' % fname
    outf = open(outfname,'w')
    outf.write('\n'.join(res))
    outf.close()

    def _get_means_stdevs(cls, x, y):
        x_y_counter_lin = cls._convert_x_y_to_counter(x, y)
        x_y_counter = cls._convert_x_y_to_counter(x, [ln(y_i) for y_i in y])

        st_dev = {x: ln(stdev(y) if stdev(y) > 0 else 1 ** -10) for x, y in x_y_counter_lin.items()}
        mean_ = {x: mean(y) for x, y in x_y_counter.items()}
        return cls._get_mean_stdev_from_counter(x_y_counter, st_dev, mean_)

def logjointprob(Z, X, A, phi, prior):
    Z = minus1(Z)
    X = minus1(X)
    return sum([ln(prior[Z[0]]), ln(phi[X[0]][Z[0]])] + 
                 [ln(A[Z[e]][Z[e-1]]) + ln(phi[X[e]][Z[e]]) 
                     for e in xrange(1, len(Z))] + 
                 [ln(A[2][Z[-1]])])

def genColor (n, startpoint=0):
    assert n >= 1
    # This splits the 0 - 1 segment in the pizza way
    h = (2*n-1)/(2**ceil(ln(n)/ln(2)))-1
    h = (h + startpoint) % 1
    # We set saturation based on the amount of green, in the range 0.6 to 0.8
    rgb = colorsys.hsv_to_rgb(h, 1, 1)
    rgb = colorsys.hsv_to_rgb(h, 1, (1-rgb[1])*0.2+0.6)
    return rgb

    def get_mean_stdev_r_squared(self, x, y):
        # remove all x,y's if y is less than or equal to 0
        yy, xx = ExponentialDistributionFunction._remove_non_positive_values(y, x)

        x_, y_mean, y_std = ExponentialDistributionFunction._get_means_stdevs(x, [y_i for y_i in yy])
        # these have to be of the form y = c x + d, therefore we use ln(y) = b x + ln(a)
        # which (in theory) is equivalent to y = a e^(bx)
        mu = get_r_squared(x_, y_mean, self.mean_b, ln(self.mean_a))
        sigma = get_r_squared(x_, y_std, self.stdev_b, ln(self.stdev_a))
        return mu, sigma

def solver(mesh, deg):
    # Physical parameters
    dpdx = Constant(-1)
    mu = Constant(100)

    V = FunctionSpace(mesh, "Lagrange", deg)
    u = TrialFunction(V)
    v = TestFunction(V)

    # Mark boundary subdomians
    class Sides(SubDomain):
        def inside(self, x, on_boundry):
            return on_boundry

    side = Sides()

    mf = FacetFunction("size_t", mesh)
    mf.set_all(2)

    side.mark(mf, 1)
    noslip = DirichletBC(V, Constant(0), mf, 1)

    a = inner(grad(u), grad(v)) * dx
    L = -1.0 / mu * dpdx * v * dx

    u_ = Function(V)
    solve(a == L, u_, bcs=noslip)

    # Compute the flux
    Q = assemble(u_ * dx)

    # Flux from analytical expression
    mu = 100
    dpdx = -1
    F = (A ** 2 - B ** 2 + C ** 2) / (2 * C)
    M = sqrt(F ** 2 - A ** 2)
    alpha = 0.5 * ln((F + M) / (F - M))
    beta = 0.5 * ln((F - C + M) / (F - C - M))
    s = 0
    for n in range(1, 100):
        s += (n * exp(-n * (beta + alpha))) / sinh(n * beta - n * alpha)

    Q_analytical = (
        (pi / (8 * mu)) * (-dpdx) * (A ** 4 - B ** 4 - (4 * C * C * M * M) / (beta - alpha) - 8 * C * C * M * M * s)
    )

    Q_error = abs(Q - Q_analytical)

    print "Flux computed numerically : ", Q
    print "Flux computed using (3-52): ", Q_analytical

    return mesh.hmin(), Q_error

def ext(n):
  n = int(n * ln(n) + n * (ln(ln(n))))

  candidates = list(range(n+1))
  fin = int(n**0.5)

  # Loop over the candidates, marking out each multiple.
  for i in xrange(2, fin+1):
      if candidates[i]:
          candidates[2*i::i] = [None] * (n//i - 1)

  # Filter out non-primes and return the list.
  return [i for i in candidates[2:] if i]

		def K(feed_size,interest):
			scale_point = float(self.personalization_scale_point)
			how_many = float(self.personalization_scale_how_many)
			desired_ratio = float(1)#float(overall_interest)
			#print feed_size, 1- (		(ln(desired_ratio*std_size/feed_size) / float(self.interest_ratio_x)) 	)	
			
			return 		-(ln(desired_ratio*how_many/feed_size) / float(scale_point))

 def read(self, addr):
     debug('debugHIH6130', self.name, "read", addr)
     try:
         data = self.interface.readBlock((self.addr, 0), 4)
         status = (data[0] & 0xc0) >> 6
         humidity = ((((data[0] & 0x3f) << 8) + data[1]) * 100.0) / 16383.0
         temp = (((data[2] & 0xff) << 8) + (data[3] & 0xfc)) / 4
         tempC = (temp / 16384.0) * 165.0 - 40.0
         tempF = tempC * 1.8 + 32
         # https://en.wikipedia.org/wiki/Dew_point
         # 0C <= T <= +50C
         b = 17.368
         c = 238.88
         gamma = ln(humidity / 100) + ((b * tempC) / (c + tempC))
         dewpointC = c * (gamma) / (b - gamma)
         dewpointF = dewpointC * 1.8 + 32
         debug("debugHIH6130", "humidity:", humidity, "tempC:", tempC, "tempF:", tempF, "dewpointC:", dewpointC, "dewpointF:", dewpointF)
         if addr == "humidity":
             return humidity
         elif addr == "temp":
             return tempF
         elif addr == "dewpoint":
             return dewpointF
         else:
             return 0
     except:
         return 0

 def create_from_x_y_coordinates(cls, x, y, distribution_type: type(Distribution)=NormalDistribution):
     xx, yy = LogLinearDistributionFunction._remove_non_positive_values(x, y)
     x_, y_mean, y_std = LogLinearDistributionFunction._get_means_stdevs([ln(x_i) for x_i in xx],
                                                                         [y_i / xx[i] for i, y_i in enumerate(yy)])
     mean_a, mean_b = linear_regression(x_, y_mean)
     stdev_a, stdev_b = linear_regression(x_, y_std)
     return cls(mean_a, mean_b, stdev_a, stdev_b)

def _probabilistic_sum(number_of_dice, sides):
    """
    For inordinately large numbers of dice, we can approximate their
    sum by picking a random point under the bell-curve that represents
    the probabilistic sums.

    We accomplish this by picking a random y value on the curve, then
    picking a random x value from the bounds of that y intercept.
    """
    n = number_of_dice
    s = sides
    u = ((s + 1.) / 2) * n  # mean
    B = (1.7 * (n ** .5) * ((2 * pi) ** .5))
    max_y = 1. / B
    min_y = (e ** ((-(n - u) ** 2) / (2 * 1.7 * 1.7 * n))) / B
    Y = random.uniform(min_y, max_y)

    try:
        T = ln(Y * B) * (2 * (1.7 * 1.7) * n)
    except ValueError:
        # Too close to 0, rounding off
        T = 0
        min_x, max_x = n, n * s
    else:
        min_x, max_x = _quadratic(1, -2 * u, T + u ** 2)
    return int(round(random.uniform(min_x, max_x)))

def classify(X, prior, mu, sigma, covdiag=True):
    h = [0] * len(X)
    Ls = [np.linalg.cholesky(sigma[:, :, k]) for k in range(len(mu))]
    for i, x_star in enumerate(X):
        maxval = float("-inf")
        selectedClass = -1
        for k, muk in enumerate(mu):
            sigma_k = sigma[:, :, k]
            if covdiag == True:
                diff = x_star - muk
                y = np.linalg.solve(sigma_k, diff.T)
            else:
                L = Ls[k]
                diff = x_star - muk
                v = np.linalg.solve(L, diff.T)
                y = np.linalg.solve(L.T, v)
            (sign, logdet) = numpy.linalg.slogdet(sigma_k)
            a = -0.5 * sign * logdet
            b = -0.5 * numpy.dot(diff, y)
            c = ln(prior[k])
            total = a + b + c
            if total > maxval:
                selectedClass = k
                maxval = total
        h[i] = selectedClass
    return numpy.array(h)

 def setValue(self, value, context):
     base = self.argA.getValue(context)
     exp = self.argB.getValue(context)
     if base:
         if exp: # a,b
             if value==base**exp:
                 pass # but overconstrained
             else:
                 print("Error! " + self.name + " is overconstrained")
         else: # base => exp
             self.argB.setValue(math.ln(value)/math.ln(base), context)
     else:
         if exp: #exp => base
             self.argA.setValue(value**(1/exp), context) # TODO allow multiple solutions...
         else: # neither => can't do anything
             pass

def get_log_lbf_from_regression_coefficients(a: float, b: float, max_x: float=130000)->tuple:
    x = list()
    y = list()
    for i in float_range(0.1, max_x, 0.5):
        x.append(i)
        y.append(a * ln(i) + b)
    return x, y

def mine(n):
  # Upper bound on the nth prime value
  # http://en.wikipedia.org/wiki/Prime_number_theorem#Approximations_for_the_nth_prime_number
  n = int(n * ln(n) + n * (ln(ln(n))))
  primes = [2]
  
  def is_prime(primes, x):
    for y in primes:
      if x % y == 0:
        return False
    return True
  
  for x in xrange(3, n, 2):
    if is_prime(primes, x):
      primes.append(x)
  return primes

def computeShannon(aList, plus=1):
    #print aList
    if len(aList) < 2:
        print >>sys.stderr, "You may need to specify \
-I parameter if the program stops."
    expr = [float(i)+1 for i in aList]
    expr_sum = sum(expr)
    assert expr_sum != 0
    expr_R = [1.0 * i / expr_sum for i in expr]
    expr_Log = []
    for i in expr_R:
        if i != 0:
            expr_Log.append(i*ln(i)/ln(2))
        else:
            expr_Log.append(i)
    shannon = -1 * sum(expr_Log)
    return shannon