def triple_plot(cccsum, cccsum_hist, trace, threshold, save=False,
                savefile=''):
    r"""Main function to make a triple plot with a day-long seismogram, \
    day-long correlation sum trace and histogram of the correlation sum to \
    show normality.

    :type cccsum: numpy.ndarray
    :param cccsum: Array of the cross-channel cross-correlation sum
    :type cccsum_hist: numpy.ndarray
    :param cccsum_hist: cccsum for histogram plotting, can be the same as \
        cccsum but included if cccsum is just an envelope.
    :type trace: obspy.Trace
    :param trace: A sample trace from the same time as cccsum
    :type threshold: float
    :param threshold: Detection threshold within cccsum
    :type save: bool, optional
    :param save: If True will svae and not plot to screen, vice-versa if False
    :type savefile: str, optional
    :param savefile: Path to save figure to, only required if save=True
    """
    if len(cccsum) != len(trace.data):
        print('cccsum is: ' +
              str(len(cccsum))+' trace is: '+str(len(trace.data)))
        msg = ' '.join(['cccsum and trace must have the',
                        'same number of data points'])
        raise ValueError(msg)
    df = trace.stats.sampling_rate
    npts = trace.stats.npts
    t = np.arange(npts, dtype=np.float32) / (df * 3600)
    # Generate the subplot for the seismic data
    ax1 = plt.subplot2grid((2, 5), (0, 0), colspan=4)
    ax1.plot(t, trace.data, 'k')
    ax1.axis('tight')
    ax1.set_ylim([-15 * np.mean(np.abs(trace.data)),
                  15 * np.mean(np.abs(trace.data))])
    # Generate the subplot for the correlation sum data
    ax2 = plt.subplot2grid((2, 5), (1, 0), colspan=4, sharex=ax1)
    # Plot the threshold values
    ax2.plot([min(t), max(t)], [threshold, threshold], color='r', lw=1,
             label="Threshold")
    ax2.plot([min(t), max(t)], [-threshold, -threshold], color='r', lw=1)
    ax2.plot(t, cccsum, 'k')
    ax2.axis('tight')
    ax2.set_ylim([-1.7 * threshold, 1.7 * threshold])
    ax2.set_xlabel("Time after %s [hr]" % trace.stats.starttime.isoformat())
    # ax2.legend()
    # Generate a small subplot for the histogram of the cccsum data
    ax3 = plt.subplot2grid((2, 5), (1, 4), sharey=ax2)
    ax3.hist(cccsum_hist, 200, normed=1, histtype='stepfilled',
             orientation='horizontal', color='black')
    ax3.set_ylim([-5, 5])
    fig = plt.gcf()
    fig.suptitle(trace.id)
    fig.canvas.draw()
    if not save:
        plt.show()
        plt.close()
    else:
        plt.savefig(savefile)
    return

def plot_feat_hist(data_name_list, filename=None):
    if len(data_name_list)>1:
        assert filename is not None

    pylab.figure(num=None, figsize=(8, 6))
    num_rows = 1 + (len(data_name_list) - 1) / 2
    num_cols = 1 if len(data_name_list) == 1 else 2
    pylab.figure(figsize=(5 * num_cols, 4 * num_rows))

    for i in range(num_rows):
        for j in range(num_cols):
            pylab.subplot(num_rows, num_cols, 1 + i * num_cols + j)
            x, name = data_name_list[i * num_cols + j]
            pylab.title(name)
            pylab.xlabel('Value')
            pylab.ylabel('Fraction')
            # the histogram of the data
            max_val = np.max(x)
            if max_val <= 1.0:
                bins = 50
            elif max_val > 50:
                bins = 50
            else:
                bins = max_val
            n, bins, patches = pylab.hist(
                x, bins=bins, normed=1, facecolor='blue', alpha=0.75)

            pylab.grid(True)

    if not filename:
        filename = "feat_hist_%s.png" % name.replace(" ", "_")

    pylab.savefig(os.path.join(CHART_DIR, filename), bbox_inches="tight")

def display_grid(grid, **kwargs):
    fig = plt.figure()
    plt.axes().set_aspect('equal')

    if kwargs.get('mark_core_cells', True):
        core_cell_coords = grid._cell_nodes[1:-1, 1:-1]
        cellx, celly = core_cell_coords[:, :, 0], core_cell_coords[:, :, 1]
        plt.plot(cellx, celly, '-o', np.transpose(cellx), np.transpose(celly), '-o', color='red')

    if kwargs.get('mark_boundary_cells', True):
        boundary_cell_coords = grid._cell_nodes[0, :], \
                               grid._cell_nodes[-1, :], \
                               grid._cell_nodes[1:-1, 0], \
                               grid._cell_nodes[1:-1, -1]

        for coords in boundary_cell_coords:
            plt.plot(coords[:, 0], coords[:, 1], '-x', color='blue')

    if kwargs.get('show', False):
        plt.show()

    f = BytesIO()
    plt.savefig(f)

    return f

 def test_simple_gen(self):
     self_con = .8
     other_con = 0.05
     g = self.gen.gen_stoch_blockmodel(min_degree=1, blocks=5, self_con=self_con, other_con=other_con,
                                       powerlaw_exp=2.1, degree_seq='powerlaw', num_nodes=1000, num_links=3000)
     deg_hist = vertex_hist(g, 'total')
     res = fit_powerlaw.Fit(g.degree_property_map('total').a, discrete=True)
     print 'powerlaw alpha:', res.power_law.alpha
     print 'powerlaw xmin:', res.power_law.xmin
     if len(deg_hist[0]) != len(deg_hist[1]):
         deg_hist[1] = deg_hist[1][:len(deg_hist[0])]
     print 'plot degree dist'
     plt.plot(deg_hist[1], deg_hist[0])
     plt.xscale('log')
     plt.xlabel('degree')
     plt.ylabel('#nodes')
     plt.yscale('log')
     plt.savefig('deg_dist_test.png')
     plt.close('all')
     print 'plot graph'
     pos = sfdp_layout(g, groups=g.vp['com'], mu=3)
     graph_draw(g, pos=pos, output='graph.png', output_size=(800, 800),
                vertex_size=prop_to_size(g.degree_property_map('total'), mi=2, ma=30), vertex_color=[0., 0., 0., 1.],
                vertex_fill_color=g.vp['com'],
                bg_color=[1., 1., 1., 1.])
     plt.close('all')
     print 'init:', self_con / (self_con + other_con), other_con / (self_con + other_con)
     print 'real:', gt_tools.get_graph_com_connectivity(g, 'com')

def make_corr1d_fig(dosave=False):
    corr = make_corr_both_hemi()
    lw=2; fs=16
    pl.figure(1)#, figsize=(8, 7))
    pl.clf()
    pl.xlim(4,300)
    pl.ylim(-400,+500)    
    lambda_titles = [r'$20 < \lambda < 30$',
                     r'$30 < \lambda < 40$',
                     r'$\lambda > 40$']
    colors = ['blue','green','red']
    for i in range(3):
        corr1d, rcen = corr_1d_from_2d(corr[i])
        ipdb.set_trace()
        pl.semilogx(rcen, corr1d*rcen**2, lw=lw, color=colors[i])
        #pl.semilogx(rcen, corr1d*rcen**2, 'o', lw=lw, color=colors[i])
    pl.xlabel(r'$s (Mpc)$',fontsize=fs)
    pl.ylabel(r'$s^2 \xi_0(s)$', fontsize=fs)    
    pl.legend(lambda_titles, 'lower left', fontsize=fs+3)
    pl.plot([.1,10000],[0,0],'k--')
    s_bao = 149.28
    pl.plot([s_bao, s_bao],[-9e9,+9e9],'k--')
    pl.text(s_bao*1.03, 420, 'BAO scale')
    pl.text(s_bao*1.03, 370, '%0.1f Mpc'%s_bao)
    if dosave: pl.savefig('xi1d_3bin.pdf')

    def test_flux(self):
        tol = 150.
        inputcat = catalog.read(os.path.join(self.args.tmp_path, 'ccd_1.cat'))
        pixradius = 3*self.target["psf"]/self.instrument["PIXEL_SCALE"]
        positions = list(zip(inputcat["X_IMAGE"]-1, inputcat["Y_IMAGE"]-1))
        fluxes = image.simple_aper_phot(self.im[1], positions, pixradius)
        sky_background = image.annulus_photometry(self.im[1], positions,
        	pixradius+5, pixradius+8)

        total_bg_pixels = np.shape(image.build_annulus_mask(pixradius+5, pixradius+8, positions[0]))[1]
        total_source_pixels = np.shape(image.build_circle_mask(pixradius,
        	positions[0]))[1]

        estimated_fluxes = fluxes - sky_background*1./total_bg_pixels*total_source_pixels

        estimated_magnitude = image.flux2mag(estimated_fluxes,
        	self.im[1].header['SIMMAGZP'], self.target["exptime"])

        expected_flux = image.mag2adu(17.5, self.target["zeropoint"][0],
        	exptime=self.target["exptime"])

        p.figure()
        p.hist(fluxes, bins=50)
        p.title('Expected flux: {:0.2f}, mean flux: {:1.2f}'.format(expected_flux, np.mean(estimated_fluxes)))
        p.savefig(os.path.join(self.figdir,'Fluxes.png'))

        assert np.all(np.abs(fluxes-expected_flux) < tol)

def Animate(g):
    for i in range(1,64,5):
        pylab.clf()
        x= g[0:i,:,:]
        y= numarray.sum(x, axis=0)
        pylab.matshow( y)
        pylab.savefig("temp/3dturb-%03i.png" % i)

def plot_part2(filename):
	"""
	Plots the result of count ones test
	"""
	fig1 = pl.figure()
	iterations, runtimes, fvals = extract(filename)
	algos = ["SA", "GA", "MIMIC"]
	iters_sa, iters_ga, iters_mimic = [np.array(iterations[a]) for a in algos]
	runtime_sa, runtime_ga, runtime_mimic = [np.array(runtimes[a]) for a in algos]
	fvals_sa, fvals_ga, fvals_mimic = [np.array(fvals[a]) for a in algos]

	plotfunc = getattr(pl, "loglog")
	plotfunc(runtime_sa, fvals_sa, "bs", mew=0)
	plotfunc(runtime_ga, fvals_ga, "gs", mew=0)
	plotfunc(runtime_mimic, fvals_mimic, "rs", mew=0)

	# plotfunc(iters_sa, fvals_sa/(runtime_sa * iters_sa), "bs", mew=0)
	# plotfunc(iters_ga, fvals_ga/(runtime_ga * iters_ga), "gs", mew=0)
	# plotfunc(iters_mimic, fvals_mimic/(runtime_mimic * iters_mimic), "rs", mew=0)

	pl.xlabel("Runtime (seconds)")
	pl.ylabel("Objective function value")
	pl.ylim([min(fvals_sa) / 2, max(fvals_mimic) * 2])
	pl.legend(["SA", "GA", "MIMIC"], loc=4)

	pl.savefig(filename.replace(".csv", ".png"), bbox_inches="tight") 

    def plot_knowledge_count(agent_network, filename):

        word_dict = dict()
        agent_list = agent_network.get_all_agents()
        for agent_item in agent_list:
            for word in agent_item.knowledge:
                if word not in word_dict:
                    word_dict[word] = 0
                word_dict[word] = word_dict[word] + 1

        word_count_tuple_list = word_dict.items()
        word_count_tuple_list = sorted(word_count_tuple_list, key=itemgetter(1))
        print word_count_tuple_list

        x = list()
        y = list()

        for item in word_count_tuple_list:
            word = item[0]
            count = item[1]
            x.append(word)
            y.append(count)

        plt.scatter(x, y, s=30, vmin = 0, vmax= 100, alpha=0.5)
        plt.savefig(filename)

def plot_size_of_c(size_of_c, path):
    xlabel('|C|')
    ylabel('Max model size |Ci|')
    grid(True)
    plot([x+1 for x in range(len(size_of_c))], size_of_c)
    savefig(os.path.join(path, 'size_of_c.png'))
    close()

def plot_q(frame,file_prefix='claw',file_format='petsc',path='./_output/',plot_pcolor=True,plot_slices=True,slices_xlimits=None):
    import sys
    sys.path.append('.')
    import gaussian_1d

    sol=Solution(frame,file_format=file_format,read_aux=False,path=path,file_prefix=file_prefix)
    x=sol.state.grid.x.centers
    mx=len(x)

    bathymetry = 0.5
    eta=sol.state.q[0,:] + bathymetry

    if frame < 10:
        str_frame = "00"+str(frame)
    elif frame < 100:
        str_frame = "0"+str(frame)
    else:
        str_frame = str(frame)

    fig = pl.figure(figsize=(40,10))
    ax = fig.add_subplot(111)
    ax.set_aspect(aspect=1)
    ax.plot(x,eta)
    #pl.title("t= "+str(sol.state.t),fontsize=20)
    #pl.xticks(size=20); pl.yticks(size=20)
    #pl.xlim([0, gaussian_1d.Lx])
    pl.ylim([0.5, 1.0])
    pl.xlim([0., 4.0])
    #pl.axis('equal')
    pl.savefig('./_plots/eta_'+str_frame+'_slices.png')
    pl.close()

def plot_running_time(running_time, path):
    xlabel('|C|')
    ylabel('MTV iteration in secs.')
    grid(True)
    plot([x for x in range(len(running_time))], running_time)
    savefig(os.path.join(path, 'running_time.png'))
    close()

def plot_heuristic(heuristic, path):
    xlabel('|C|')
    ylabel('h')
    grid(True)
    plot(heuristic)
    savefig(os.path.join(path, 'heuristic.png'))
    close()

def plot_BIC_score(BIC_SCORE, path):
    xlabel('|C|')
    ylabel('BIC score')
    grid(True)
    plot(BIC_SCORE)
    savefig(os.path.join(path, 'BIC.png'))
    close()

def plotHist(outBaseStr,outData,motifID):
    """Plot histograms for each motif results."""
    def _chooseBins():
        """Uses one of a number of heuristics to 
        choose the bins for each hist.
        
        More discussed here:
        http://en.wikipedia.org/wiki/Histogram#Number_of_bins_and_width"""
        # sqrtChoice
        # k = int(sqrt(n))+1
        return int(math.sqrt(len(ctrls)))*3 

        
    outFile = '%s_%s_hist.pdf' % (outBaseStr,motifID)
    ctrls   = [x for x in outData[motifID][1:]]
    yLab    = 'Controls [n=%s]' % (len(ctrls))
    xLab    = 'Enrichment [p-value]'
    
    fig  = pl.figure()
    ax = fig.add_subplot(111)
    ax.set_ylabel(yLab)
    ax.set_xlabel(xLab)
    ax.set_xscale('log')
    
    bins = _chooseBins()
    hData = ax.hist(ctrls, bins=bins,histtype='stepfilled',color='grey',cumulative=1)
    ax.axvline(x=outData[motifID][0] ,linewidth=2, color='r')
    pl.savefig(outFile)

def plotFeaturePDF(ift, pft, outbase, fmin=0.0, fmax=1.0, fstep=0.01):
    """
    Plot a comparison between the input feature distribution and the 
    feature distribution of the predicted halos
    """
    plt.clf()
    nfbins = ( fmax - fmin ) / fstep
    fbins = np.logspace( fmin, fmax, nfbins )
    fcen = ( fbins[:-1] + fbins[1:] ) / 2

    plt.xscale( 'log', nonposx='clip' )
    plt.yscale( 'log', nonposy='clip' )
    
    ic, e, p = plt.hist( ift, fbins, label='Original Halos', alpha=0.5, normed=True )
    pc, e, p = plt.hist( pft, fbins, label='Added Halos', alpha=0.5, normed=True )

    plt.legend()
    plt.xlabel( r'$\delta$' )
    plt.savefig( outbase+'_fpdf.png' )

    fdtype = np.dtype( [ ('fcen', float), ('ifcounts', float), ('pfcounts', float) ] )
    fd = np.ndarray( len(fcen), dtype = fdtype )
    fd[ 'mcen' ] = fcen
    fd[ 'imcounts' ] = ic
    fd[ 'pmcounts' ] = pc

    fitsio.write( outbase+'_fpdf.fit', fd )

    def visualization2(self, sp_to_vis=None):
        if sp_to_vis:
            species_ready = list(set(sp_to_vis).intersection(self.all_sp_signatures.keys()))
        else:
            raise Exception('list of driver species must be defined')

        if not species_ready:
            raise Exception('None of the input species is a driver')

        for sp in species_ready:
            # Setting up figure
            plt.figure()
            plt.subplot(313)

            mon_val = OrderedDict()
            signature = self.all_sp_signatures[sp]
            for idx, mon in enumerate(list(set(signature))):
                if mon[0] == 'C':
                    mon_val[self.all_comb[sp][mon] + (-1,)] = idx
                else:
                    mon_val[self.all_comb[sp][mon]] = idx

            mon_rep = [0] * len(signature)
            for i, m in enumerate(signature):
                if m[0] == 'C':
                    mon_rep[i] = mon_val[self.all_comb[sp][m] + (-1,)]
                else:
                    mon_rep[i] = mon_val[self.all_comb[sp][m]]
            # mon_rep = [mon_val[self.all_comb[sp][m]] for m in signature]

            y_pos = numpy.arange(len(mon_val.keys()))
            plt.scatter(self.tspan[1:], mon_rep)
            plt.yticks(y_pos, mon_val.keys())
            plt.ylabel('Monomials', fontsize=16)
            plt.xlabel('Time(s)', fontsize=16)
            plt.xlim(0, self.tspan[-1])
            plt.ylim(0, max(y_pos))

            plt.subplot(312)

            for name in self.model.odes[sp].as_coefficients_dict():
                mon = name
                mon = mon.subs(self.param_values)
                var_to_study = [atom for atom in mon.atoms(sympy.Symbol)]
                arg_f1 = [numpy.maximum(self.mach_eps, self.y[str(va)][1:]) for va in var_to_study]
                f1 = sympy.lambdify(var_to_study, mon)
                mon_values = f1(*arg_f1)
                mon_name = str(name).partition('__')[2]
                plt.plot(self.tspan[1:], mon_values, label=mon_name)
            plt.ylabel('Rate(m/sec)', fontsize=16)
            plt.legend(bbox_to_anchor=(-0.1, 0.85), loc='upper right', ncol=1)

            plt.subplot(311)
            plt.plot(self.tspan[1:], self.y['__s%d' % sp][1:], label=parse_name(self.model.species[sp]))
            plt.ylabel('Molecules', fontsize=16)
            plt.legend(bbox_to_anchor=(-0.15, 0.85), loc='upper right', ncol=1)
            plt.suptitle('Tropicalization' + ' ' + str(self.model.species[sp]))

            # plt.show()
            plt.savefig('s%d' % sp + '.png', bbox_inches='tight', dpi=400)

def plotMassFunction(im, pm, outbase, mmin=9, mmax=13, mstep=0.05):
    """
    Make a comparison plot between the input mass function and the 
    predicted projected correlation function
    """
    plt.clf()

    nmbins = ( mmax - mmin ) / mstep
    mbins = np.logspace( mmin, mmax, nmbins )
    mcen = ( mbins[:-1] + mbins[1:] ) /2
    
    plt.xscale( 'log', nonposx = 'clip' )
    plt.yscale( 'log', nonposy = 'clip' )
    
    ic, e, p = plt.hist( im, mbins, label='Original Halos', alpha=0.5, normed = True)
    pc, e, p = plt.hist( pm, mbins, label='Added Halos', alpha=0.5, normed = True)
    
    plt.legend()
    plt.xlabel( r'$M_{vir}$' )
    plt.ylabel( r'$\frac{dN}{dM}$' )
    #plt.tight_layout()
    plt.savefig( outbase+'_mfcn.png' )
    
    mdtype = np.dtype( [ ('mcen', float), ('imcounts', float), ('pmcounts', float) ] )
    mf = np.ndarray( len(mcen), dtype = mdtype )
    mf[ 'mcen' ] = mcen
    mf[ 'imcounts' ] = ic
    mf[ 'pmcounts' ] = pc

    fitsio.write( outbase+'_mfcn.fit', mf )

def PlotTurbulenceIllustr(a):

    """
    Can generate the grid with
    g=kolmogorovutils.GenerateKolmogorov3D( 1025, 129, 129)
    a=kolmogorovutils.GridToNumarray(g)

    """

    for x in [1,10,100]:

        suba= numarray.sum(a[:,:,0:x], axis=2)

        suba.transpose()
        pylab.clf()
        pylab.matshow(suba)
        pylab.savefig("temp/turb3d-sum%03i.eps" % x)

    for x in [1,10,100]:

        for j in [0,1,2]:

            suba= numarray.sum(a[:,:200,j*x:(j+1)*x], axis=2)

            suba.transpose()
            pylab.clf()
            pylab.matshow(suba)
            pylab.savefig("temp/turb3d-sum%03i-s%i.eps" % (x,j))        

    def handle(self, *args, **options):
        try:
            from matplotlib import pylab as pl
            import numpy as np
        except ImportError:
            raise Exception('Be sure to install requirements_scipy.txt before running this.')

        all_names_and_counts = RawCommitteeTransactions.objects.all().values('attest_by_name').annotate(total=Count('attest_by_name')).order_by('-total')
        all_names_and_counts_as_tuple_and_sorted = sorted([(row['attest_by_name'], row['total']) for row in all_names_and_counts], key=lambda row: row[1])
        print "top ten attestors:  (name, number of transactions they attest for)"
        for row in all_names_and_counts_as_tuple_and_sorted[-10:]:
            print row

        n_bins = 100
        filename = 'attestor_participation_distribution.png'

        x_max = all_names_and_counts_as_tuple_and_sorted[-31][1]  # eliminate top outliers from hist
        x_min = all_names_and_counts_as_tuple_and_sorted[0][1]

        counts = [row['total'] for row in all_names_and_counts]
        pl.figure(1, figsize=(18, 6))
        pl.hist(counts, bins=np.arange(x_min, x_max, (float(x_max)-x_min)/100) )
        pl.title('Histogram of Attestor Participation in RawCommitteeTransactions')
        pl.xlabel('Number of transactions a person attested for')
        pl.ylabel('Number of people')
        pl.savefig(filename)