def badness(self):
     """
     How bad is the error?
     
     Returns the max of the absolute logarithmic errors of kf and Kc
     """
     return max(abs(math.log10(self.k_ratio)), abs(math.log10(self.Keq_ratio)))

def metal_con(filename, distances, real_dist, bins=35, limits=(-3.1, 0.2),
              avgs=1, detection=1, tag="out"):
    """ main bit """
    if filename[-4:] == '.csv':  delim = ','
    else:  delim = None
    data = shift_data(fi.read_data(filename, delim), real_dist, distances[0])
    mod_actual = 5.*(ma.log10(real_dist*1000) - 1.)
    mod_new = 5.*(ma.log10(distances[0]*1000) - 1.)
    mod = mod_actual - mod_new
    print "Effective Magnitude Shift = {0}, average g={1}".format(mod, sc.mean(data[:,2]))
    new_data = cut_data(data, 4,2,3,5, deff=0, modulus=mod, full=1)
    FeH = get_photo_metal(new_data[:,4],new_data[:,2],new_data[:,3])
    ref_hist = np.histogram(FeH, bins, limits)
    hist = []
    #Also iterate over several runs and average
    for i in range(len(distances)):
        print "#- Convolving to distance {0} kpc".format(distances[i])
        if i==0:  deff=0
        else: deff=detection
        temp_hist = []
        for j in range(avgs):
            #holds dist constant, applies appropriate errors for new distance
            new_data = con.convolve(data, real_dist, distances[i])
            #shift data so detection efficiency works correctly;  has no noticable effect if deff=0
            new_data = shift_data(new_data, distances[0], distances[i])
            # apply color cuts and detection efficiency to shifted and convolved data
            new_data = cut_data(new_data, 4,2,3,5, deff=deff, modulus=None, full=0)
            print "Average g = {0}, total stars = {1}".format(sc.mean(new_data[:,2]), len(new_data[:,0]))
            FeH = get_photo_metal(new_data[:,4],new_data[:,2],new_data[:,3])
            temp_hist.append(np.histogram(FeH, bins, limits))
        new_hist = avg_hists(temp_hist)
        hist.append(new_hist)
    plot_hists(hist, ref_hist, distances, tag)
    return hist

 def auto_precision(self):
     """Automatically find the precision based on the format of the tick
     label.
     """
     
     if self.ticklabel.format.lower()=="general":
         self.ticklabel.prec = 0 # doesn't matter
     elif self.ticklabel.format.lower()=="decimal":
         x = self.tick.major
         p = int(math.floor(math.log10(x)))
         if p>=0:
             self.ticklabel.prec = 0
         else:
             z = math.floor(x/(10**p))
             y = x-z*10**p
             if y==0.0:
                 self.ticklabel.prec = -p
             else:
                 prec = -int(math.floor(math.log10(y)+0.0000001))
                 self.ticklabel.prec = prec
     elif self.ticklabel.format.lower()=="power":
         self.ticklabel.prec = 0
     elif self.ticklabel.format.lower() in ["exponential","scientific"]:
         x = self.tick.major
         p = int(math.floor(math.log10(x)))
         z = math.floor(x/(10**p))
         y = x-z*10**p
         if y==0.0:
             self.ticklabel.prec = 0
         else:
             prec = -int(math.floor(math.log10(y/10**p)+0.0000001))
             self.ticklabel.prec = prec
     else:
         self.ticklabel.prec = 1

def lporder (freq1, freq2, delta_p, delta_s):
    '''
    FIR lowpass filter length estimator.  freq1 and freq2 are
    normalized to the sampling frequency.  delta_p is the passband
    deviation (ripple), delta_s is the stopband deviation (ripple).

    Note, this works for high pass filters too (freq1 > freq2), but
    doesnt work well if the transition is near f == 0 or f == fs/2

    From Herrmann et al (1973), Practical design rules for optimum
    finite impulse response filters.  Bell System Technical J., 52, 769-99
    '''
    df = abs (freq2 - freq1)
    ddp = math.log10 (delta_p)
    dds = math.log10 (delta_s)

    a1 = 5.309e-3
    a2 = 7.114e-2
    a3 = -4.761e-1
    a4 = -2.66e-3
    a5 = -5.941e-1
    a6 = -4.278e-1

    b1 = 11.01217
    b2 = 0.5124401

    t1 = a1 * ddp * ddp
    t2 = a2 * ddp
    t3 = a4 * ddp * ddp
    t4 = a5 * ddp

    dinf=((t1 + t2 + a3) * dds) + (t3 + t4 + a6)
    ff = b1 + b2 * (ddp - dds)
    n = dinf / df - ff * df + 1
    return n

def add_stat_text(axhi, weights, bins) :
    #mean, rms, err_mean, err_rms, neff = proc_stat(weights,bins)
    mean, rms, err_mean, err_rms, neff, skew, kurt, err_err = proc_stat(weights,bins)
    pm = r'$\pm$' 
    txt  = 'Mean=%.2f%s%.2f\nRMS=%.2f%s%.2f\n' % (mean, pm, err_mean, rms, pm, err_rms)
    txt += r'$\gamma1$=%.3f  $\gamma2$=%.3f' % (skew, kurt)
    #txt += '\nErr of err=%8.2f' % (err_err)
    xb,xe = axhi.get_xlim()     
    yb,ye = axhi.get_ylim()     
    #x = xb + (xe-xb)*0.84
    #y = yb + (ye-yb)*0.66
    #axhi.text(x, y, txt, fontsize=10, color='k', ha='center', rotation=0)
    x = xb + (xe-xb)*0.98
    y = yb + (ye-yb)*0.95

    if axhi.get_yscale() is 'log' :
        #print 'axhi.get_yscale():', axhi.get_yscale()
        log_yb, log_ye = log10(yb), log10(ye)
        log_y = log_yb + (log_ye-log_yb)*0.95
        y = 10**log_y

    axhi.text(x, y, txt, fontsize=10, color='k',
              horizontalalignment='right',
              verticalalignment='top',
              rotation=0)

    def __init__(self, image, normalize=True, title=None, parent=None):
        qt.QFrame.__init__(self, parent)
        self.viewer = ImageViewer(image, normalize, title, parent=self)
        self.setWindowTitle(self.viewer.windowTitle())

        self._captionCoords = 0, 0
        self._xplaces = int(math.log10(self.viewer.image.width) + 1.0)
        self._yplaces = int(math.log10(self.viewer.image.height) + 1.0)
        self._valueplaces = self.viewer.image.channels * 5

        self.label = qt.QLabel(self)
        font = qt.QFont()
        font.setPointSize(10)
        font.setStyleHint(qt.QFont.TypeWriter)
        self.label.setFont(font)

        self._layout = qt.QVBoxLayout(self)
        self._layout.setSpacing(5)
        self._layout.addWidget(self.viewer, 1)
        self._layout.addWidget(self.label)

        self.connect(self.viewer, SIGNAL('mouseOver(int, int)'), self.updateCaption)
        self.connect(self.viewer.cursorAction, SIGNAL('triggered()'), self._toggleCaptionSignals)

        self.updateCaption()

    def get_median_mag(self, area, rake):
        """
        Return magnitude (Mw) given the area and rake.

        Setting the rake to ``None`` causes their "All" rupture-types
        to be applied.

        :param area:
            Area in square km.
        :param rake:
            Rake angle (the rupture propagation direction) in degrees,
            from -180 to 180.
        """
        assert rake is None or -180 <= rake <= 180
        if rake is None:
            # their "All" case
            return 4.07 + 0.98 * log10(area)
        elif (-45 <= rake <= 45) or (rake > 135) or (rake < -135):
            # strike slip
            return 3.98 + 1.02 * log10(area)
        elif rake > 0:
            # thrust/reverse
            return 4.33 + 0.90 * log10(area)
        else:
            # normal
            return 3.93 + 1.02 * log10(area)

 def score_response(self, comment, response):
     """
     This function can be modified to give a good internal scoring for a
     response. If negative, we won't post. This is useful when there is more
     than one possible response in our database.
     """
     # Discard the obviously bad responses
     if response.body.strip() == "[deleted]":
         return -1
     simple_body = rewriter.simplify_body(comment.body)
     if response.score < config.good_comment_threshold:
         return -1
     # Derive our base score. We use a logarithm, because reddit scores are
     # roughly logrithmic <http://amix.dk/blog/post/19588>
     base_score = math.log10(response.score)
     # A raw pentalty to subtract for the comment being in a different
     # context from it's parent
     response_parent = self.__get_parent(response)
     if response_parent is not None:
         similarity = self._get_parent_similarity_ratio(comment, response_parent)
         difference_penalty = math.log10(10000) * (1 - similarity) ** 10
     else:
         difference_penalty = math.log10(10000)
     # give it some points for length
     length_reward = math.log10(len(simple_body))
     # throw in some randomness for good luck
     fuzz_multiplier = random.gauss(mu=1, sigma=0.05)
     # put it all together
     final_score = (base_score - difference_penalty + length_reward) * fuzz_multiplier
     return final_score

def _rescale(lo,hi,step,pt=None,bal=None,scale='linear'):
    """
    Rescale (lo,hi) by step, returning the new (lo,hi)
    The scaling is centered on pt, with positive values of step
    driving lo/hi away from pt and negative values pulling them in.
    If bal is given instead of point, it is already in [0,1] coordinates.

    This is a helper function for step-based zooming.
    """
    # Convert values into the correct scale for a linear transformation
    # TODO: use proper scale transformers
    if scale=='log':
        lo,hi = log10(lo),log10(hi)
        if pt is not None: pt = log10(pt)

    # Compute delta from axis range * %, or 1-% if persent is negative
    if step > 0:
        delta = float(hi-lo)*step/100
    else:
        delta = float(hi-lo)*step/(100-step)

    # Add scale factor proportionally to the lo and hi values, preserving the
    # point under the mouse
    if bal is None:
        bal = float(pt-lo)/(hi-lo)
    lo = lo - bal*delta
    hi = hi + (1-bal)*delta

    # Convert transformed values back to the original scale
    if scale=='log':
        lo,hi = pow(10.,lo),pow(10.,hi)

    return (lo,hi)

 def p_powerset(self, p):
     '''powerset : POWER EQ default
                 | POWER EQ value
                 | POWER EQ value dbm
                 | POWER EQ value w
                 | POWER EQ value mw
                 | POWER EQ value uw
                 | POWER EQ value nw'''
     from math import log10
     if len(p) == 4: # Default or unitless value
         p[0] = {'power' : p[3]} # Return power, even if just 'default', since user explicitly set it
     else:   # Units specified
         if p[4] == 'dbm':
             p[0] = {'power' : p[3]}
         elif p[4] == 'w':
             # Convert value in W to dBm
             p[0] = {'power' : 10*log10(p[3]*1e+3)}
         elif p[4] == 'mw':
             # Convert value in mW to dBm
             p[0] = {'power' : 10*log10(p[3])}
         elif p[4] == 'uw':
             # Convert value in uW to dBm
             p[0] = {'power' : 10*log10(p[3]*1e-3)}
         elif p[4] == 'nw':
             # Convert value in nW to dBm
             p[0] = {'power' : 10*log10(p[3]*1e-6)}
         else: # How'd we get here?! Assume units of dBm
             print("ERROR: Undefined case in grammar of powerset. Assuming units of dBm.")   # TODO: throw exception here
             p[0] = {'power' : p[3]}

    def interact(self):
        # counterclockwise
        if self.left_pressed:
            self.angle = (self.angle + 10) % 360

        if self.right_pressed:
            self.angle = self.angle - 10
            if self.angle < 0:
                self.angle += 360

        if self.up_pressed:
            self.speed = math.log10(self.speed + 2) * 5
        elif self.speed > 1:
            self.speed = math.log10(self.speed) * 5
        else:
            self.speed = 0
        
        if self.space_pressed:
            self.bomb(0, 2)
            self.bomb(30, 5)
            self.bomb(-30 , 5)
            self.bomb(15, 15)
            self.bomb(-15 , 15)
            self.bomb(50, 30)
            self.bomb(-50 , 30)

        self.bombs += 0.2

        self.rotate()
        return True

def main():
    ht = 50
    hr = 2
    f = 900 * 10**6
    c = 3 * 10**8
    gr = [1, 0.316, 0.1, 0.01]
    gl = 1
    r = -1

    distance = np.arange(1, 100001, 1, dtype=float)
    lambd = c / float(f)
    dc = 4 * ht * hr / lambd
    reflect = (distance**2 + (ht + hr)**2)**0.5
    los = (distance**2 + (ht - hr)**2)**0.5
    phi = 2 * pi * (reflect - los) / lambd
    flat = distance[:ht]
    decline = distance[ht:(dc + 1)]
    steep = distance[dc:]

    for i in range(len(gr)):
        temp = gl**0.5 / los + r * (gr[i]**0.5) * np.exp(phi * -1J) / reflect
        pr = (lambd / 4 / pi)**2 * (abs(temp)**2)
        plt.subplot(220 + i + 1)
        plt.plot(10*log(distance),10*log(pr)-10*log10(pr[0]),'b', \
                10*log(flat), np.zeros(len(flat)), 'y', \
                10*log(decline),-20*log(decline),'g', 10*log(steep),-40*log(steep),'r')
        plt.axvline(x=10 * log10(ht), linestyle='-.')
        plt.axvline(x=10 * log10(dc), linestyle='-.')
        plt.title("Gr = %s" % gr[i])

    plt.show()

 def Apparent_Magnitude_Absolute_Magnitude_Distance_Relation(self, var, val1, val2):
     if var == 'm': #val1 = M, val2 = d
         return(val1 + 5 * math.log10(val2 - 5))
     if var == 'M': #val1 = m, val2 = d
         return(val1 - 5 * math.log10(val2 - 5))
     if var == 'd': #val1 = m, val2 = M
         return(math.pow(10, (val1 - val2) / 5) + 5)

def fm_heat(gene2heat, fm_threshold, cis_threshold=0.01, CIS=False):
    print "* Creating oncodrive heat map..."
    if CIS:
        print "\tIncluding CIS scores at threshold", cis_threshold, "..."
    heat = dict()
    src_fm, src_cis_amp, src_cis_del = 0, 0, 0
    for g, scores in gene2heat.items():
        if CIS:
            del_score = scores["del"] if scores["del"] < cis_threshold else NULL
            amp_score = scores["amp"] if scores["amp"] < cis_threshold else NULL
            fm_score = scores["fm"] if scores["fm"] < fm_threshold else NULL
            if fm_score == NULL and amp_score == NULL and del_score == NULL:
                continue
            min_val = min(del_score, amp_score, fm_score)
            heat[g] = -log10(min_val)
            if min_val == scores["fm"]:
                src_fm += 1
            elif min_val == scores["amp"]:
                src_cis_amp += 1
            elif min_val == scores["del"]:
                src_cis_del += 1
        else:
            if scores["fm"] >= fm_threshold:
                continue
            heat[g] = -log10(scores["fm"])
            src_fm += 1
    print "\t- Genes using FM score:", src_fm
    print "\t- Genes using CIS AMP score:", src_cis_amp
    print "\t- Genes using CIS DEL score:", src_cis_del

    return heat

def psnr(sp_img_1, sp_img_2):
    """
    Peak Signal To Noise Ratio - measure of image quality

    Parameters
    ----------
    :param *SpatialImage* sp_img_1: *SpatialImage* image

    :param *SpatialImage* sp_img_2: *SpatialImage* image

    Returns
    ----------
    :return: float psnr -- psnr value (dB)
    """
    if isinstance(sp_img_1, SpatialImage) and isinstance(sp_img_2, SpatialImage):
        if sp_img_1.itemsize==sp_img_2.itemsize:
            maxi =  2**(sp_img_1.itemsize*8) - 1
            mse = mean_squared_error(sp_img_1, sp_img_2)
            if mse!=0:
                psnr = 20.0*log10(maxi) - 10*log10(mse)
            elif mse==0:
                psnr = np.inf
            return psnr
        else:
            print('sp_img_1 and sp_img_2 does not have the same type')
            return
    else:
        print('sp_img_1 and sp_img_2 must be SpatialImage instances')
        return

def find_tickmarks_log(interval, n=11):
    '''
    Find logarithmicly spaced tickmarks.

    Arguments:

        * `interval` -- interval to produce tickmarks for (two element
          iterable)
        * `n` -- Number of tickmarks to produce (default 11) - Is ignored,
          present only for compatibility with the non logarithmic tickmark
          finder.

    Note:
        Tickmarks returned are a list [(value_1, size_1), ...,
        (value_n, size_n)].
    '''
    l, h = interval
    tickmarks = []
    for i in range(int(log10(l) - 1), int(log10(h)) + 1):
        for ii in range(1, 10):
            tickmark_position = ii * 10 ** i
            if l <= tickmark_position <= h:
                if ii == 1:
                    size = 10
                else:
                    size = ii
                tickmarks.append((tickmark_position, size))
    return tickmarks

    def config_axes(self, xlog, ylog):
        if hasattr(self, "_rng"):
            (i1, j1, i2, j2) = self.visible_area()
            zoomed = 1
        else:
            zoomed = 0

        self._xlog = xlog
        self._ylog = ylog
        if xlog:
            self._rng = [log10(x) for x in self._original_rng]
        else:
            self._rng = self._original_rng
        if ylog:
            self._vals = [log10(x) for x in self._original_vals]
        else:
            self._vals = self._original_vals

        self._imin = min(self._rng)
        self._imax = max(self._rng)
        if self._imax == self._imin:
            self._imin -= 1
            self._imax += 1
        self._jmin = min(self._vals)
        self._jmax = max(self._vals)
        if self._jmax == self._jmin:
            self._jmin -= 1
            self._jmax += 1

        if zoomed:
            self.zoom(i1, j1, i2, j2)
        else:
            self.zoom(self._imin, self._jmin, self._imax, self._jmax)

def slice_frequencies_into_log_pockets(bin_key, bins):

	'''
	Within each bin, separates the frequencies again into logarithmically built pockets, 
	using the global NUM_POCKETS variable as determined above. Note that frequency is defined
	by location, hence the use of enumerate.

	'''

	bin_location = bins[bin_key]
	max_frequencies_in_bin = []
	amplitudes_in_pocket = []
	frequencies_in_bin = []

	max_log_idx = math.log10(len(bin_location))
	pocket_size = float(max_log_idx)/NUM_POCKETS

	#pockets is a list of lists--that is each of the pockets in a bin. 
	pockets = [ [] for x in range(NUM_POCKETS) ]
	
	for frequency, amplitude in enumerate(bin_location):
		if frequency == 0:
			continue
		log_index = math.log10(frequency)
		pocket_idx = int(log_index/pocket_size)
		pockets[min(pocket_idx, NUM_POCKETS-1)].append((abs(amplitude), frequency))
	return pockets

def plot_contours(f, x1, x2, y1, y2, z1, z2,
                  z_logscale, new_figure=True):
	"""Plots contour lines of the given objective function. The 
plotting range is [x1,x2]x[y1,y2]x[z1,z2], and logarithmic 
scaling of z-axis is specified with the boolean argument 
z_logscale. If new_figure is set to true, this function 
opens a new window for the plot.
"""
	X,Y,Z = tabulate_function(f, 300, (x1, x2), (y1, y2))
  
	if new_figure == True:
		fig = figure()
  
	if z_logscale == True:
		V = logspace(math.log10(z1), math.log10(z2), 15)
	else:
		V = arange(z1, z2, (z2 - z1) / 20)
	
	contour(X, Y, Z, V, colors='k', linewidths=0.25)
  
	xlim(x1, x2)
	ylim(y1, y2)
  
	#if isinstance(f, TestFunction):
		#title(f.name)
	#elif isinstance(f, native.Function) and f.has_symbolic_expression():
		#title(f.get_symbolic_expression())
	
	if new_figure == True:
		show()

def compute_cluster_score(args):
    """Computes the cluster score for a given set type"""
    cluster, cutoff = args
    set_type = SET_SET_TYPE
    matrix = SET_MATRIX
    ref_matrix = SET_REF_MATRIX
    cluster_rows = SET_MEMBERSHIP.rows_for_cluster(cluster)
    cluster_genes = [gene for gene in cluster_rows if gene in set_type.genes()]
    overlap_sizes = []
    set_sizes = []

    for set_name, eset in set_type.sets.items():
        set_genes = eset.genes_above_cutoff()
        intersect = set(cluster_genes).intersection(set_genes)
        overlap_sizes.append(len(intersect))
        set_sizes.append(len(set_genes))

    num_sets = len(set_type.sets)
    phyper_n = np.array([len(set_type.genes()) for _ in xrange(num_sets)]) - np.array(set_sizes)
    phyper_n = [value for value in phyper_n]
    phyper_k = [len(cluster_genes) for _ in xrange(num_sets)]
    enrichment_pvalues = list(util.phyper(overlap_sizes, set_sizes, phyper_n, phyper_k))
    min_pvalue = min(enrichment_pvalues)
    min_index = enrichment_pvalues.index(min_pvalue)
    min_set = set_type.sets.keys()[min_index]
    min_set_overlap = overlap_sizes[min_index]
    if min_set_overlap > 0:
        scores = [0.0 for _ in xrange(matrix.num_rows())]
        min_genes = set_type.sets[min_set].genes
        min_genes = [gene for gene in min_genes if gene in matrix.row_names]
        min_indexes = matrix.row_indexes(min_genes)

        if set_type.sets[min_set].cutoff == "discrete":
            overlap_genes = set(cluster_genes).intersection(set(min_genes))
            overlap_indexes = matrix.row_indexes(overlap_genes)
            for index in min_indexes:
                scores[index] = 0.5
            for index in overlap_indexes:
                scores[index] = 1.0
        else:
            min_set_weights = []
            for index in min_indexes:
                min_set_weights.append(set_type.sets[min_set].weights[index])
            min_weight = min(min_set_weights)
            max_weight = max(min_set_weights)
            for index in min_indexes:
                scores[index] = min_set_weights[index] - min_weight
                scores[index] = min_set_weights[index] / max_weight

        dampened_pvalue = enrichment_pvalues[min_index]
        if dampened_pvalue <= cutoff:
            dampened_pvalue = 1
        else:
            dampened_pvalue = math.log10(dampened_pvalue) / math.log10(cutoff)
        scores = [dampened_pvalue / score if score != 0.0 else score for score in scores]
        min_ref_score = ref_matrix.min()
        scores = [score * min_ref_score for score in scores]
    else:
        scores = [0.0 for _ in xrange(matrix.num_rows())]
    return scores, min_set, min_pvalue