def plot_singleCurrStepResp(filename, condition, Voltages, meanVolt, fitData, Rin, tauFit_ms, plateau, DC, Cap, pre, step, volt_dt, ax0, ay0, dx, dy, xwidth, xswidth, ywidth, yswidth, plotPos):

    #
    # plot voltage traces and fit
    # 

    n,m=plotPos

    xAxis = np.linspace(1., len(meanVolt), len(meanVolt))
    axesPos = [ax0+(m*(dx+xwidth)), 1-ay0-n*(dy+ywidth), xwidth, ywidth]   
    ax = plt.axes(axesPos)

    for k in np.arange(np.shape(Voltages)[1]):
        plt.plot(xAxis,Voltages[:,k], color = '0.8')

    plt.plot(meanVolt, color = 'k', linewidth = 3)
    plt.plot(xAxis[pre+step:pre+step+fitData.size], fitData, color = 'g', linewidth = 3)

    ax.set_autoscale_on(False)
    plotSpace = (np.max(meanVolt) - np.min(meanVolt))*0.4
    ax.set_ylim(np.min(meanVolt)-plotSpace, np.max(meanVolt)+plotSpace)
    # ax.set_ylim((wcCurr[pre + (0.1*pre):pre + (1.8*pre),:].min())*1.3, (wcCurr[pre + (0.1*pre):pre + (1.8*pre),:].max()))
    # ax.set_ylim((wcCurr[2*pre:3*pre,:].min())*1.5, (wcCurr[2*pre:3*pre,:].max())*1.3)
    # ax.set_xlim(pre*0.5, wcCurr.shape[0])
    ax.set_frame_on(False)
    ax.set_xticks([])
    ax.set_yticks([])

    #
    # plot scale bars
    #

    # specify size of voltage and time scalebar:
    voltSclb = 5 # in mV
    timeSclb = 0.1 # in seconds

    xRange = ax.get_xlim()[1]
    yRange = ax.get_ylim()[1] - ax.get_ylim()[0]

    # scalebar coordinates (Lx: left x, Rx: right x, Ly: lower y, Uy: upper y):
    sclbLx = ax.get_xlim()[1]-(pre)
    sclbRx = sclbLx + timeSclb/volt_dt
    sclbLy = ax.get_ylim()[1]-yRange*0.5
    sclbUy = sclbLy + voltSclb

    # plt.hlines(ax.get_ylim()[0], pre*0.52, pre*0.52 + 0.1/volt_dt, colors = 'k', linestyles = 'solid', linewidth = 3)
    plt.hlines(sclbLy, sclbLx, sclbRx, colors = 'k', linestyles = 'solid', linewidth = 3)
    plt.text(sclbLx + 0.02/volt_dt, sclbLy - 0.1*yRange,'{0} ms'.format(timeSclb*1000), fontsize = 8,)
    plt.vlines(sclbLx, sclbLy, sclbUy, colors = 'k', linestyles = 'solid', linewidth = 3)
    # minYax = ax.get_ylim()[0]
    # maxYax = ax.get_ylim()[1]
    plt.text(sclbLx - 0.06*xRange, sclbLy + 0.3,'{0} mV'.format(voltSclb), fontsize = 8, rotation = 'vertical')
    # plt.text(pre+1.2*step, maxYax-(0.3*(maxYax-minYax)), 'Rin = {0}  MOhm'.format(Rin), fontsize = 12)
    # plt.text(pre+1.2*step, maxYax-(0.4*(maxYax-minYax)), 'Tau = {0} ms'.format(tauFit_ms), fontsize = 12)
    # plt.text(pre+1.2*step, maxYax-(0.5*(maxYax-minYax)), 'Cap = {0} pF'.format(Cap), fontsize = 12)
    # plt.text(pre+1.2*step, maxYax-(0.6*(maxYax-minYax)), 'DC = {0} pA'.format(DC), fontsize = 12)
    plt.text(pre+step*0.5, ax.get_ylim()[1] - (ax.get_ylim()[1]-ax.get_ylim()[0])*0.0, 'file {0},{1} baseline mp: {2:.1f}mV {3} plateuiness: {4}'.format(filename, '\n', condition, '\n', plateau), fontsize = 10)
    # plt.text(pre*0.5, ax.get_ylim()[1] + (ax.get_ylim()[1]-ax.get_ylim()[0])*0.03, ('passive responses ' + conditions[n]), fontsize = 10)

    return

def show_fixed_lag_numberline():
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_xlim(0,10)
    ax.set_ylim(0,10)

    # draw lines
    xmin = 1
    xmax = 9
    y = 5
    height = 1

    plt.hlines(y, xmin, xmax)
    plt.vlines(xmin, y - height / 2., y + height / 2.)
    plt.vlines(4.5, y - height / 2., y + height / 2.)
    plt.vlines(6, y - height / 2., y + height / 2.)
    plt.vlines(xmax, y - height / 2., y + height / 2.)
    plt.vlines(xmax-1, y - height / 2., y + height / 2.)

    # add numbers
    plt.text(xmin, y-1.1, '$x_0$', fontsize=20, horizontalalignment='center')
    plt.text(xmax, y-1.1, '$x_k$', fontsize=20, horizontalalignment='center')
    plt.text(xmax-1, y-1.1, '$x_{k-1}$', fontsize=20, horizontalalignment='center')
    plt.text(4.5, y-1.1, '$x_{k-N+1}$', fontsize=20, horizontalalignment='center')
    plt.text(6, y-1.1, '$x_{k-N+2}$', fontsize=20, horizontalalignment='center')
    plt.text(2.7, y-1.1, '.....', fontsize=20, horizontalalignment='center')
    plt.text(7.2, y-1.1, '.....', fontsize=20, horizontalalignment='center')

    plt.axis('off')
    plt.show()

def plot_percentiles(data, numbins, xlim, ylim, vert = True, color = 'k', linestyle = 'solid', linew = 2):
    perc = 1. / numbins 
    for i in range(1, numbins):
        if vert:
            plt.vlines(sts.scoreatpercentile(data, i * perc * 100.), ylim[0], ylim[1], color, linestyle, linewidth = linew)
        else:
            plt.hlines(sts.scoreatpercentile(data, i * perc * 100.), xlim[0], xlim[1], color, linestyle, linewidth = linew)

def width_feature(im, f2, plot=False):
    '''
    f7 is median of the gap lengths
    f8 is f2 / f7
    '''
    transitions = []
    for i in xrange(len(im)):
        groups = itertools.groupby(list(im[i]))
        transitions.append((i, len([x for x in groups])))

    max_line, max_trans = max(transitions, key=operator.itemgetter(1))

    gap_lengths = []
    zeros_start = 0
    zeros = True
    for i in xrange(len(im[max_line])):
        pixel = im[max_line][i]
        if pixel == 1 and zeros is True:
            gap_lengths.append(i - zeros_start)
            zeros = False
        if pixel == 0 and zeros is not True:
            zeros_start = i
            zeros = True

    if plot:
        print gap_lengths
        plt.imshow(im, cmap=cm.Greys_r)
        plt.hlines(numpy.array([max_line]), 0, im.shape[1], linewidth=3, colors='y')
        plt.show()

    f7 = numpy.median(gap_lengths)
    return f7, float(f2) / f7

def _expand_dendrogram(cNode,swc_tree,off_x,off_y,radius,transform='plain') :
    global max_width,max_height # middle name d.i.r.t.y.
    '''
    Gold old fashioned recursion... sys.setrecursionlimit()!
    '''
    place_holder_h = H_SPACE
    max_degree = swc_tree.degree_of_node(cNode)
    required_h_space = max_degree * place_holder_h
    start_x = off_x-(required_h_space/2.0)
    if(required_h_space > max_width) :
        max_width = required_h_space
    
    if swc_tree.is_root(cNode) :
        print 'i am expanding the root'
        cNode.children.remove(swc_tree.get_node_with_index(2))
        cNode.children.remove(swc_tree.get_node_with_index(3))
    
    for cChild in cNode.children :
        l = _path_between(swc_tree,cChild,cNode,transform=transform)
        r = cChild.content['p3d'].radius

        cChild_degree = swc_tree.degree_of_node(cChild)
        new_off_x = start_x + ( (cChild_degree/2.0)*place_holder_h )
        new_off_y = off_y+(V_SPACE*2)+l
        r = r if radius  else 1
        plt.vlines(new_off_x,off_y+V_SPACE,new_off_y,linewidth=r,colors=C)
        if((off_y+(V_SPACE*2)+l) > max_height) :
            max_height = off_y+(V_SPACE*2)+l

        _expand_dendrogram(cChild,swc_tree,new_off_x,new_off_y,radius=radius,transform=transform)

        start_x = start_x + (cChild_degree*place_holder_h)
        plt.hlines(off_y+V_SPACE,off_x,new_off_x,colors=C)

def OScrossPower(angSep, crossCorr, crossCorrErr, tex=True):
    if tex == True:
        plt.rcParams['text.usetex'] = True

    fig, ax = plt.subplots()

    seps = []
    rows,cols = angSep.shape
    for ii in range(rows):
        for jj in range(ii+1,cols):
           seps.append(np.arccos(angSep[ii,jj]))
            
    #ax.plot(seps, crossCorr, 'k', color='#1B2ACC')
    #ax.fill_between(seps, crossCorr-crossCorrErr, crossCorr-crossCorrErr,
           #alpha=0.2, edgecolor='#1B2ACC', facecolor='#089FFF', linewidth=4,
           #linestyle='dashdot', antialiased=True)
    ax.errorbar(np.array(seps)*180.0/np.pi, crossCorr, yerr=crossCorrErr,
                fmt='o', color='#089FFF', ecolor='#089FFF', capsize=8, linewidth=3)
    plt.hlines(y=0.0, xmin=0.0, xmax=180.0, linewidth=3.0, linestyle='solid', color='black')
    ax.set_xlabel('Pulsar angular separation [degrees]', fontsize=20)
    ax.set_ylabel('Cross power', fontsize=20)
    ax.minorticks_on()
    plt.tick_params(labelsize=18)
    plt.grid(which='major')
    plt.grid(which='minor')
    plt.xlim(0.0,180.0)
    plt.title('Cross-power measurements')
    plt.show()

def MissingMassVolume(path):
	mvd_pt = open(path+'mvd_pt.pkl','rb')
	data = pickle.load(mvd_pt)
	mvd_pt.close() 
	plt.figure()
	plt.ylabel('sum voxels')
	x = np.arange(0,len(data['amountRm']))
	plt.scatter(x,data['amountRm'],c = 'r', label='Vol voxels removed')
	plt.hold(True)
	plt.scatter(x,data['sum_before'], c='b',label='original seg Vol') 
	plt.scatter(x,data['vol_vox'], c = 'k',label='Final volume')

	m = np.mean(data['amountRm'])
	print m 
	plt.hlines(m,0,len(data['amountRm']))
	plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.05),
      fancybox=True, shadow=True, ncol=5) 
	plt.savefig(path+'MissingVol.png')
	plt.figure()
	plt.ylabel('sum Linear Attenuation')
	plt.scatter(x,data['mass'], c = 'k', label= 'Final mass')  
	plt.scatter(x,data['mass_incompressible'], c= 'g',label='Mass "incompressible" removed') 
	plt.hlines(np.mean(data['mass_incompressible']),0,len(data['amountRm']))
	plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.05),
      fancybox=True, shadow=True, ncol=5) 
	plt.savefig(path+'Missing.png')

 def draw_clade(clade, x_start, color='k', lw=1):
     """Recursively draw a tree, down from the given clade."""
     x_here = x_posns[clade]
     y_here = y_posns[clade]
     # phyloXML-only graphics annotations
     if hasattr(clade, 'color') and clade.color is not None:
         color = clade.color.to_hex()
     if hasattr(clade, 'width') and clade.width is not None:
         lw = clade.width
     # Draw a horizontal line from start to here
     plt.hlines(y_here, x_start, x_here, color=color, lw=lw)
     # Add node/taxon labels
     label = label_func(clade)
     if label not in (None, clade.__class__.__name__):
         plt.text(x_here, y_here, ' ' + label,
                 fontsize=10, verticalalignment='center')
     # Add confidence
     if hasattr(clade, 'confidences'):
         # phyloXML supports multiple confidences
         conf_label = '/'.join(map(str, map(float, clade.confidences)))
     elif clade.confidence is not None:
         conf_label = str(clade.confidence)
     else:
         conf_label = None
     if conf_label:
         plt.text(x_start, y_here, str(float(clade.confidence)), fontsize=9)
     if clade.clades:
         # Draw a vertical line connecting all children
         y_top = y_posns[clade.clades[0]]
         y_bot = y_posns[clade.clades[-1]]
         # Only apply widths to horizontal lines, like Archaeopteryx
         plt.vlines(x_here, y_bot, y_top, color=color)
         # Draw descendents
         for child in clade:
             draw_clade(child, x_here, color=color, lw=lw)

def plot_turb_lines(linestyles="solid", linewidth=2, color="gray"):
    plt.hlines(0.5, -1, 1, linestyles=linestyles, colors=color,
               linewidth=linewidth)
    plt.vlines(-1, -0.2, 0.5, linestyles=linestyles, colors=color,
               linewidth=linewidth)
    plt.vlines(1, -0.2, 0.5, linestyles=linestyles, colors=color,
               linewidth=linewidth)

 def show(self):
     import numpy as np
     import matplotlib.pyplot as plt
     from scipy.signal import freqz
     
     # Sample rate and desired cutoff frequencies (in Hz).
     fs = 5000.0
     lowcut = 500.0
     highcut = 1250.0
     order = 5
     
     bpf = BandPassFilter(lowcut, highcut, fs, order)
     
     # Plot the frequency response for a few different orders.
     plt.figure(1)
     plt.clf()
     for order in [3, 6, 9]:
         b, a = bpf.butter_bandpass(order)
         w, h = freqz(b, a, worN=2000)
         plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)
     
     plt.plot([0, 0.5 * fs], [np.sqrt(0.5), np.sqrt(0.5)], '--', label='sqrt(0.5)')
     plt.xlabel('Frequency (Hz)')
     plt.ylabel('Gain')
     plt.grid(True)
     plt.legend(loc='best')
     
     
     # Filter a noisy signal.
     T = 0.05
     nsamples = T * fs
     t = np.linspace(0, T, nsamples, endpoint=False)
     a = 0.02
     f0 = 600.0
     x = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
     x += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
     x += a * np.cos(2 * np.pi * f0 * t + .11)
     x += 0.03 * np.cos(2 * np.pi * 2000 * t)
     
     """
     wf = wave.open('../Alarm01.wav', 'rb')
     signal = wf.readframes(1024)
     signal = np.fromstring(signal, dtype=np.int16)
     x = signal[0::2]
     right = signal[1::2]
     """
     
     y = bpf.butter_bandpass_filter(x, order)
     
     plt.figure(2)
     plt.clf()
     plt.plot(t, x, label='Noisy signal')
     plt.plot(t, y, label='Filtered signal (%g Hz)' % f0)
     plt.xlabel('time (seconds)')
     plt.hlines([-a, a], 0, T, linestyles='--')
     plt.grid(True)
     plt.axis('tight')
     plt.legend(loc='upper left')
     
     plt.show()

def PlotProfile(label1, label2, data1, data2, filename, smooth, binsz):
    plt.clf()
    plt.figure(figsize=(5,4))
    plt.rc('text', usetex=True)
    plt.rc('font', family='serif')
    
    counts1 = np.append(data1['FLUX'].data/(EXPOSURE*PIXEL_SA), 0)
    glats1 = np.append(data1['GLAT_MIN'][0], data1['GLAT_MAX'].data)

    profile1 = np.histogram(glats1, bins=np.sort(glats1), weights=counts1)
    xstep=binsz
    y_smooth_1 = gaussian_filter(profile1[0], smooth / xstep)
    print y_smooth_1.mean()
    
    counts2 = np.append(data2['FLUX'].data/(EXPOSURE*PIXEL_SA*(data2['FLUX'].data.size())), 0)
    glats2 = np.append(data2['GLAT_MIN'][0], data2['GLAT_MAX'].data)

    profile2 = np.histogram(glats2, bins=np.sort(glats2), weights=counts2)
    xstep=binsz
    y_smooth_2 = gaussian_filter(profile2[0], smooth / xstep)
    print y_smooth_2.mean()
    
    x1 = 0.5 * (glats1[1:] + glats1[:-1])
    plt.plot(x1, y_smooth_1, label='{0}'.format(label1))
    plt.plot(x1, y_smooth_2, label='{0}'.format(label2))
    plt.hlines(0, LatLow+binsz, LatHigh-binsz)
    plt.xlabel(r'Galactic Latitude/$deg$', fontsize=10)
    plt.ylabel(r'Surface Brightness/ph cm$^{-2}$ s$^{-1} sr^{-1}$', fontsize=10)
    plt.xlim([LatLow+binsz, LatHigh-binsz])
    plt.tick_params(axis='x', labelsize=10)
    plt.grid(b=True, which='major', color='0.75', linewidth=0.5)
    plt.legend(prop={'size':8})
    plt.savefig(filename)

def plot_mdotdata_mdotsteadyline(simul,tnorm,mdot,args):
    '''1 plot per sink, so nsink plots
    this is a check that steady mdot calculation looks reasonable for each sink's mdot data'''
    time= simul.time/tBHA
    mdotBHA= simul.mdot.copy()/mdot.norm/args.answer
    for sid in range(8): #range(args.nsinks):
        lab=""
        if args.vrms_vs_time: lab= " vrms=f(t) corrected"
        plt.plot(time,mdotBHA,'k-',label='data'+lab)
        steady=simul.tavg_mdot[sid]/mdot.norm/args.answer
        plt.hlines(simul.tavg_mdot[sid]/mdot.norm/args.answer,time.min(),time.max(),\
            linestyles='dashed',colors='b',label="steady= %.2g" % steady)
        plt.xlabel("$\mathbf{ t/t_{bh} }$")
        ylab= r"$\mathbf{ \dot{M} / \dot{M}_0 / \dot{M}_{ans}  }$"
        ylab += r" (<$\mathbf{ \dot{M} }$> after t= %.2f tBH)" % args.tReachSteady
        plt.ylabel(ylab)
        plt.yscale('log')
    #     ax.legend(loc=4)
        if args.ylim: plt.ylim(args.ylim[0],args.ylim[1])
        if args.xlim: plt.xlim(args.xlim[0],args.xlim[1])
        if args.LinearY: plt.yscale('linear')
        plt.legend(loc=(1.01,0.01),ncol=2,fontsize='xx-small')
        fsave="mdot_%d.png" % sid
        plt.savefig(fsave) 
        plt.close()

def graphMembershipIncomeOverTime(membertimedata, memberinfos, membership):
    ydates, yvalues = zip(*membertimedata)
    ydates = list(ydates)
    yvalues = map(lambda x:x[1],yvalues)
    plotlabel = u"Membership Income over time"
    plt.plot(ydates, yvalues, '^g',linewidth=2, markevery=1)
    plt.ylabel("Euro")
    plt.xlabel("Month")
    plt.grid(True)
    plt.title(plotlabel)
    ## label with +x-y members per month
    membersinmonth = dict(membertimedata)
    #print "\n".join([ "%s:%s" % (x[0],str(x[1])) for x in extractPlusMinusFromMemberships(membership).items()])
    pm_dates, pm_fee_dates = extractPlusMinusFromMemberships(membership)
    for astrdate, tpl in pm_fee_dates.items():
        adate = datetime.datetime.strptime(astrdate, dateformat_monthonly_).date()- dateutil.relativedelta.relativedelta(months=1)
        assert(adate.day==1)
        if adate in membersinmonth:
            plt.vlines(adate + dateutil.relativedelta.relativedelta(days=11),tpl.minus+membersinmonth[adate][1],tpl.plus+membersinmonth[adate][1])
            plt.hlines(membersinmonth[adate][1], adate + dateutil.relativedelta.relativedelta(days=10),adate + dateutil.relativedelta.relativedelta(days=12))
    xstart,xend = plt.xlim()
    locs, labels = plt.xticks(np.arange(xstart,xend,61))
    plt.gca().xaxis.set_major_formatter(plt.matplotlib.dates.DateFormatter(dateformat_monthonly_, tz=None))
    plt.setp(labels, rotation=80)
    plt.subplots_adjust(left=0.06, bottom=0.08, right=0.99, top=0.95)

 def plot(self, slice=None):
     if slice == None:
         slice = len(self.cutoffs)
     plt.hlines(1,0,1)
     plt.eventplot(self.cutoffs[:slice], orientation='horizontal', colors='b')
     plt.show()
     plt.pause(0.0001)

    def plot(self):
        global TESTING, UNITS, HEIGHTS
        '''Prepares the data using self.prepare_data and then
        graphs the data on a plot.'''
        self.prepare_data()

        plt.plot(HEIGHTS, self.data)
        plt.hlines(self.significant_shear, 0, HEIGHTS[-1])
        plt.vlines(self.significant_shear_height, -1, 2)
        print 'Significant shear at image {0}'.format(self.x_significant_shear)
        if not TESTING:
            print 'Theoretical significant shear at height {0} {1}'.format(self.significant_shear_height, UNITS)

        plt.ylim([-1, 2])
        plt.xlim([HEIGHTS[0], HEIGHTS[-1]])
        
        plt.xlabel('Height ({0})'.format(UNITS))
        plt.ylabel('Coverage')
        plt.title(self.dp_path.split('/')[-1])
        
        try:
            os.mkdir('{0}/res'.format(self.dp_path))
        except:
            pass
        plt.savefig('{0}/res/results.png'.format(self.dp_path))
        with open('{0}/res/results.txt'.format(self.dp_path), 'w') as f:
            global MODE
            f.write('{0}\nMODE {1}\n'.format(str(self.significant_shear_height), MODE))

def singlet_return_with_dc_noise(**params):
	plt.figure()

	operators = ['evolution']
	#operators.append('J_hf') # UNCOMMENT THIS LINE TO INCLUDE HF NOISE AS WELL AS DC

	times = dqd.p.range('t',t=(0,'3*T_2s',1000),**params)
	times2 = dqd.p.convert(times,output='ns')
	results = dqd.measure.amplitude.iterate(ranges={'de':('-3*e_deviation','3*e_deviation',30)},int_times=times,int_initial=['ground'],params=params,int_operators=operators)
	ranges,ranges_eval,r = results.ranges, results.ranges_eval, results.results['amplitude']
	des = ranges_eval['de']
	weights = stats.norm.pdf(des,loc=0,scale=dqd.p('e_deviation',**params))
	amps = np.average(r['amplitude'],axis=0,weights=weights)[0]

	for i in xrange(r['amplitude'].shape[3]):
		plt.plot(times2, amps[:,i], label='amplitude_%d'%i)

	plt.plot(times2,[dqd.p(p_DC,t=t,**params) for t in times])
	plt.vlines([dqd.p('_T_2s',**params).value],0,1,linestyles='dashed')
	plt.grid()
	xext = plt.xlim()
	plt.xlabel("ns")
	plt.ylim(-0.1,1.1)
	plt.hlines([0.5+0.5/math.e,0.5-0.5/math.e],*xext)

	plt.savefig('noise_DC.pdf')

 def plot(self):
     plt.hlines(self.baseCr, self.crsAll.index.min(), self.crsAll.index.max(), linestyles='dotted')
     plt.plot(self.crsAll.index, self.crsAll.values)        
     plt.plot(self.crsAll.index, self.crsAll.values, 'g.')
     plt.title(str(self))
     plt.xlabel("Date")
     plt.ylabel("Creatinine")

def animal_scatterplot(json_filename, lw=3, s=60, hline_width=0.16):
   """
   Reproduces the animal data from Kaschube et al. 2010 as a scatter
   plot in a similar style. One key difference is that the area of the
   markers no longer signifies anything.
   """
   fig = plt.figure(figsize=(8,6))
   ax = plt.subplot(111)
   plt.axhline(y=math.pi, linestyle='dotted')
   d = json.load(open(json_filename,'r'))['Kaschube 2010']

   animals = ('Ferret','Tree shrew','Galago')
   colors = ('#27833a','#7fcbf1', '#f9ab27')
   for (animal, color) in zip(animals,colors):
      point_kws = dict(edgecolor=color, marker='D', facecolor=(0,0,0,0))
      plt.scatter(d[animal]['x'], d[animal]['y'], s=s, lw=lw, **point_kws)
      medx = np.median(d[animal]['x'])
      medy = np.median(d[animal]['y'])
      hline_kws = dict(colors=color, linestyles='solid', lw=lw)
      plt.hlines(medy, medx-hline_width, medx+hline_width, **hline_kws)

   ax.get_xaxis().set_visible(False)
   ax.get_yaxis().set_visible(False)
   [spine.set_visible(False) for spine in ax.spines.values()]
   plt.ylim((0, 12)); plt.xlim((0, 1.1))
   return fig

def plot_simp(n):
    xi, wi = qnwsimp(n+1, xmin, xmax)
    
    fig = plt.figure()
    plt.plot(x, f(x), linewidth=3, label=r'$f(x)$')
    
    for k in range(n//2):
        xii = xi[(2*k):(2*k+3)]
        xiii = linspace(xii[0], xii[2], 125)
        p = fitquad(xii)
        plt.fill_between(xiii, p(xiii), color='yellow')    
        if k==0:
            plt.plot(xiii, p(xiii),'r--', label=r'$\tilde{f}_{%d}(x)$' % (n+1))
        else:
            plt.plot(xiii, p(xiii),'r--')
    
    plt.vlines(xi, 0, f(xi),'k', linestyle=':')
    plt.hlines(0,xmin-0.1, xmax+0.1,'k',linewidth=2)
    plt.xlim(xmin-0.1, xmax+0.1)
    xtl = ['$x_{%d}$' % i for i in range(n+1)]
    xtl[0] += '=a'
    xtl[n] += '=b'
    plt.xticks(xi, xtl)
    plt.yticks([0],['0'])
    plt.legend()
    return fig

def model_scatterplot(iterables, colors=['r','b'],
                      mean_mm2_ferret=0.75,
                      lw=3, s=60, hline_width=0.16):
   """
   Scatterplot of simulation data for the GCAL and L models to overlay
   over the animal data. The two iterables are to contain zipped kmax
   and pinwheel densities.

   mean_mm2_ferret - mean hypercolumn size in millimeters squared
                     Sets a value for the models (uncalibrated).
   """
   fig = plt.figure(figsize=(8,6))
   ax = plt.subplot(111)
   plt.axhline(y=math.pi, linestyle='dotted')

   for iterable, color in zip(iterables, colors):
      units_per_hc, pw_densities = zip(*iterable)
      units_per_hc2 = np.array(units_per_hc)**2

      # Center both clusters on the given mean hypercolumn size for ferret.
      HCmm2 = units_per_hc2 * (mean_mm2_ferret / units_per_hc2.mean())
      point_kws = dict(edgecolor=color, marker='o', facecolor=(0,0,0,0))
      plt.scatter(HCmm2, pw_densities, s=s,lw=lw, **point_kws)

      medx = np.median(pw_densities); medy = np.median(HCmm2)
      hline_kws = dict(colors=color, linestyles='solid', lw=lw)
      plt.hlines(medy, medx-hline_width, medx+hline_width, **hline_kws)

   ax.get_xaxis().set_visible(False)
   ax.get_yaxis().set_visible(False)
   [spine.set_visible(False) for spine in ax.spines.values()]
   plt.ylim((0, 12)); plt.xlim((0, 1.1))
   return fig