def cnt_map(region, table = 'sample', draw = True):
    """Draw a region map of tweets"""
    twt_lst = dataset.loadrows(GEOTWEET, ('lat', 'lng'),
            ('MBRContains({0}, geo)'.format(dataset.geo_rect(*region)),), table)
    lat = list();
    lng = list();
    for twt in twt_lst:
        lat.append(twt['lat'])
        lng.append(twt['lng'])
    if draw:
        x = np.array(lng)
        y = np.array(lat)
        xmin = x.min()
        xmax = x.max()
        ymin = y.min()
        ymax = y.max()

        plt.hexbin(x,y, gridsize=200, cmap=cm.jet)
        plt.axis([xmin, xmax, ymin, ymax])
        plt.title("Hexagon binning")
        cb = plt.colorbar()
        cb.set_label('counts')

        plt.show()
    return lat, lng

def MakeFigures(pool, firsts, others):
    """Creates several figures for the book."""


    # plot CDFs of birth ages for first babies and others
    line_options = [
                    dict(linewidth=0.5),
                    dict(linewidth=0.5)
                    ]

    myplot.Cdfs([firsts.age_cdf, others.age_cdf], 
                root='nsfg_age_cdf',
                line_options=line_options, 
                title="Mother's age CDF",
                xlabel='age (years)',
                ylabel='probability')

    # make a scatterplot of ages and weights
    ages, weights = GetAgeWeight(pool)
    pyplot.clf()
    #pyplot.scatter(ages, weights, alpha=0.2)
    pyplot.hexbin(ages, weights, cmap=matplotlib.cm.gray_r)
    myplot.Save(root='age_scatter',
                xlabel='Age (years)',
                ylabel='Birth weight (oz)',
                legend=False)

def psych_chart(T, W=None, RH=None, heatplot=False, lims = (0, 90, 0, 0.02), **kwargs):
    if (W==None) & (RH==None):
        return None # TODO exception?
    elif W==None:
        W = humidity_ratio(RH, T)
    fig = plt.figure()
    if heatplot:
        plt.hexbin(T, W, extent=lims, **kwargs)
    else:
        plt.scatter(T, W, **kwargs)
    ts = np.linspace(lims[0],lims[1],21)
    if heatplot:
        linecolor = 'w'
    else:
        linecolor = 'k'
    for rh in np.linspace(0.1,1,10):
        plt.plot(ts, humidity_ratio(rh, ts), linecolor)
    plt.xlabel('Temperature [degF]')
    plt.ylabel('Humidity Ratio')
    plt.ylim(lims[2], lims[3])
    plt.xlim(lims[0], lims[1])
    # RH labels on right hand side
    ax1 = plt.axes()
    ax2 = ax1.twinx()
    ax2.get_yaxis().set_major_locator(ticker.FixedLocator(humidity_ratio(np.linspace(0,1,11), lims[1])/(lims[3]-lims[2])))
    ax2.get_yaxis().set_major_formatter(ticker.FixedFormatter(np.linspace(0,100,11)))
    plt.ylabel('Relative Humidity [%]')
    return ax1

def SummarizeWeight(rows, input_limit=None):
    years = [1981, 1982, 1985, 1986, 1988, 1989, 1990, 1992, 
             1993, 1994, 1996, 1998, 2000, 2002, 2004, 2006, 2008]
    
    all_diffs = []
    
    for i, row in enumerate(rows):
        if i == input_limit:
            break
        
        id, race, sex = row[:3]
        weights = row[3:]
        print id
        diffs = Differences(years, weights, jitter=3)
        all_diffs.extend(diffs)
        
    weights, changes = zip(*all_diffs)
    
    print 'Mean weight', thinkstats.Mean(weights)
    print 'Mean change', thinkstats.Mean(changes)
    print numpy.corrcoef(weights, changes)

    pyplot.hexbin(weights, changes, cmap=matplotlib.cm.gray_r)
    myplot.Save('nlsy_scatter',
                title = 'Weight change vs. weight',
                xlabel = 'Current weight (pounds)',
                ylabel = 'Weight change (pounds)',
                axis = [70, 270, -25, 25],
                legend=False,
                show=True,
                )

 def _hex_bin(root, heights, weights, cmap=matplotlib.cm.Blues):
     pyplot.hexbin(heights, weights, cmap=cmap)
     _05_myplot._save(root=root,
                      xlabel='Height (cm)',
                      ylabel='Weight (kg)',
                      axis=[140, 210, 20, 200],
                      legend=False)

def draw2():

    import numpy as np
    import matplotlib.pyplot as plt

    np.random.seed(0)
    n = 100000
    x = np.random.standard_normal(n)
    y = 2.0 + 3.0 * x + 4.0 * np.random.standard_normal(n)
    xmin = x.min()
    xmax = x.max()
    ymin = y.min()
    ymax = y.max()

    plt.subplots_adjust(hspace=0.5)
    plt.subplot(121)
    plt.hexbin(x, y, cmap=plt.cm.YlOrRd_r)
    plt.axis([xmin, xmax, ymin, ymax])
    plt.title("Hexagon binning")
    cb = plt.colorbar()
    cb.set_label('counts')
    plt.subplot(122)
    plt.hexbin(x, y, bins='log', cmap=plt.cm.YlOrRd_r)
    plt.axis([xmin, xmax, ymin, ymax])
    plt.title("With a log color scale")
    cb = plt.colorbar()
    cb.set_label('log10(N)')
    plt.show()

def make_figure(standard_id, scores, collection):
	print_fig_info()
	fig_file = os.path.join(args.output, '{0}_{1}_{2}.pdf'.format(args.db, collection, standard_id))
	x = [100.0 - s[0] for s in scores]
	y = [s[1] for s in scores]
	xmin = min(x)
	xmax = max(x)
	ymin = min(y)
	# ymax = max(y)
	# plot params
	plt.subplots_adjust(hspace=0.95)
	plt.subplot(111)
	plt.hexbin(x, y, bins='log', cmap=mpl.cm.jet, mincnt=2, gridsize=100)
	plt.title(standard_id, fontsize=18)
	# set and label axes
	plt.axis([xmin-2, xmax+2, ymin-2, 102])
	# plt.gca().invert_xaxis()
	plt.xlabel('Germline divergence')
	plt.ylabel('{0} identity'.format(standard_id))
	# make and label the colorbar
	cb = plt.colorbar()
	cb.set_label('Sequence count (log10)', labelpad=10)
	# save figure and close
	plt.savefig(fig_file)
	plt.close()
	print_done()

def rama_plot(selection='all', from_state=1, to_state=1, step=1, scatter=True):
    """ 
    Makes a scatter plot with the phi and psi angle pairs
    """
    first, last = pose_from_pdb(selection)
    if first or last:
        bonds = get_glyco_bonds(first, last)
    
        con_matrix = writer(bonds)
    
        phi = []
        psi = []
        for state in range(from_state, to_state+1, step):
            for element in con_matrix:
                phi.append(get_phi(selection, element, state))
                psi.append(get_psi(selection, element, state))

        if scatter:
            plt.scatter(phi, psi)
        else:
            gridsize=100
            #gridsize = int(2*len(phi)**(1/3))
            #if gridsize < 36:
            #    gridsize = 36
            plt.hexbin(phi, psi, gridsize=gridsize, cmap=plt.cm.summer, mincnt=1)
    
        plt.xlabel('$\phi$', fontsize=16)
        plt.ylabel('$\psi$', fontsize=16, rotation=0)
        plt.xlim(-180, 180) 
        plt.ylim(-180, 180)
        plt.show()

def chunks():
	kincr = []
	pincr = []
	for doc in collection.find(parms):
		ktr = trajectory(doc['_id'],False,False,'karma')
		tr = trajectory(doc['_id'],False,False)
		for i in range(1,len(tr)):
			kincr.append(ktr[i] - ktr[i-1])
			pincr.append(tr[i] - tr[i-1])
	
	#log scaler
	nkincr = []
	for k in kincr:
		if k < 0: nkincr.append(-log(-k+1))
		elif k == 0: nkincr.append(0)
		elif k > 0 : nkincr.append(log(k+1))
	npincr = []
	for p in pincr:
		if p < 0: npincr.append(-log(-p+1))
		elif p == 0: npincr.append(0)
		elif p > 0 : npincr.append(log(p+1))
	plt.hexbin(npincr,nkincr,bins='log',mincnt=1)
	plt.colorbar()
	plt.title('Trajectory movement')
	plt.xlabel('log scaled change of position')
	plt.ylabel('log scaled change of karma')
	plt.show()

def main():

    print "Loading file ..."
    with open("inverse_kinematics.txt", "rb") as handle:
        iv = pickle.loads(handle.read())
    print "Done."

    coord_x = []
    coord_y = []

    for coord in iv.keys():
        for _ in range(len(iv[coord])):
            coord_x.append(coord[0])
            coord_y.append(coord[1])

    x = np.array(coord_x)
    y = np.array(coord_y)

    xmin = x.min()
    xmax = x.max()
    ymin = y.min()
    ymax = y.max()

    gridsize = 60

    plt.hexbin(x, y, gridsize=gridsize, cmap=plt.cm.jet, bins=None)
    plt.axis([xmin, xmax, ymin, ymax])
    plt.subplot(111)
    cb = plt.colorbar()

    plt.show()

def plot_column_pair(i, num_columns, save_dir, titles, data, gmm_means, refcoords):
  for j in range(i+1, num_columns):
    plt.hexbin(data[:,i],  data[:,j], bins = 'log', mincnt=1)
    print(gmm_means)
    print(gmm_means[i])
    print(gmm_means[j])
    for mean in gmm_means[i]:
      plt.axvline(x=mean,color='k',ls='dashed')
    for mean in gmm_means[j]:
      plt.axhline(y=mean,color='k',ls='dashed')
    if refcoords is not None:
      plt.scatter([refcoords[0,i]], [refcoords[0,j]], marker = 's', c='w',s=15)
      plt.scatter([refcoords[1,i]], [refcoords[1,j]], marker = 'v', c='k',s=15)
    if titles is not None: 
      pp = PdfPages("%s/%s_%s.pdf" %(save_dir, titles[i], titles[j]))
      plt.xlabel(titles[i])
      plt.ylabel(titles[j])
      pp.savefig()
      pp.close()
      plt.clf()
    else:
      pp = PdfPages("%s/tIC.%d_tIC.%d.pdf" %(save_dir, i+1, j+1))
      plt.xlabel("tIC.%d" %(i+1))
      plt.ylabel("tIC.%d" %(j+1))
      pp.savefig()
      pp.close()
      plt.clf()

def bivar_hexbin_latlon(lat, lon):
    fig, ax = plt.subplots(figsize=(13, 7))
    # set major ticks every 20
    x_major_ticks = np.arange(0, 361, 20)
    ax.set_xticks(x_major_ticks)
    y_major_ticks = np.arange(-90, 91, 10)
    ax.set_yticks(y_major_ticks)

    ax.set_xlim([0, 360])
    ax.set_ylim([-90, 90])

    plt.hexbin(lon,
               lat,
               gridsize=140,
               cmap=mpl.cm.plasma,
               bins=None)

    cb = plt.colorbar()
    cb.set_label('Frequency')

    ax.set_title("Crater Distribution on Mars Surface (System: Planetocentric)",fontsize=14)
    ax.set_xlabel("LONGITUDE East (degrees)",fontsize=12)
    ax.set_ylabel("LATITUDE (degrees)",fontsize=12)

    # show plot figure
    #plt.show()
    # save plot figure
    filename = 'graph_6.png'
    plt.savefig(filename, dpi = 150)
    plt.close()

def  my_plot_histogram_hexbin(x, y, 
                              bins='log', 
                              gridsize=None,
                              ttitle="", txlabel='x', tylabel='y', ymin=None,
                              file_name=None):
    """Plot a 2d-histogram using hexbins
    """
    fig = plt.figure(figsize = (9,8))
    ax = fig.add_subplot(111)
    xdata = x.ravel()
    ydata = y.ravel()
    if gridsize is not None:
        plt.hexbin(xdata,ydata,cmap=plt.cm.gray_r,bins='log',gridsize=gridsize)
    else:
        plt.hexbin(xdata,ydata,cmap=plt.cm.gray_r,bins='log')    
    plt.colorbar()
    plt.title(ttitle, color='k')
    if ymin is not None:
        ax.set_ylim(ymin,np.max(ydata))
    plt.xlabel(txlabel)
    plt.ylabel(tylabel)
    if file_name is not None:
        plt.savefig(file_name)
    plt.show()
    return 1

def main():
    trajs = load_trajs()
    pca = fit_pca(trajs)


    cmap2 = matplotlib.cm.hot_r
    #cmap1 = brewer2mpl.get_map('OrRd', 'Sequential', 9).mpl_colormap

    for i, mfn in enumerate(model_fns):
        msm = MarkovStateModel.load(mfn)
        try:
            X = np.vstack([trajs[str(t)] for t in msm.traj_filenames])
        except KeyError as e:
            logging.exception('Not found? round %d', i)
            continue

        X_reduced = pca.transform(X)
        pp.hexbin(X_reduced[:,0], X_reduced[:,1], cmap=cmap2, bins='log',
                  vmin=0, vmax=VMAX)

        pp.xlabel('PC 1')
        pp.ylabel('PC 2')
        pp.title('Round %d PCA: %d frames' % (i, X_reduced.shape[0]))
        pp.xlim(*XLIM)
        pp.ylim(*YLIM)
        cb = pp.colorbar()
        cb.set_label('log10(N)')
        
        fn = 'plot_%04d.png' % i
        print fn
        pp.savefig(fn)
        pp.clf()

    subprocess.check_output('echo y | avconv -i plot_%04d.png -c:v libx264 -preset slow -crf 1 output.avi', shell=True)

def plot_tica_and_clusters(component_j, transformed_data, clusterer, lag_time, component_i, label = "dot", active_cluster_ids = [], intermediate_cluster_ids = [], inactive_cluster_ids = [], tica_dir = ""):

	trajs = np.concatenate(transformed_data)
	plt.hexbin(trajs[:,component_i], trajs[:,component_j], bins='log', mincnt=1)
	plt.xlabel("tIC %d" %(component_i + 1))
	plt.ylabel('tIC %d' %(component_j+1))
	centers = clusterer.cluster_centers_
	indices = [j for j in range(0,len(active_cluster_ids),1)]

	for i in [active_cluster_ids[j] for j in indices]:
		center = centers[i,:]
		if label == "dot":
			plt.scatter([center[component_i]],[center[component_j]],  marker='v', c='k', s=10)
		else:
			plt.annotate('%d' %i, xy=(center[component_i],center[component_j]), xytext=(center[component_i], center[component_j]),size=6)
	indices = [j for j in range(0,len(intermediate_cluster_ids),5)]
	for i in [intermediate_cluster_ids[j] for j in indices]:
		center = centers[i,:]
		if label == "dot":
			plt.scatter([center[component_i]],[center[component_j]],  marker='8', c='m', s=10)
		else:
			plt.annotate('%d' %i, xy=(center[component_i],center[component_j]), xytext=(center[component_i], center[component_j]),size=6)
	indices = [j for j in range(0,len(inactive_cluster_ids),5)]
	for i in [inactive_cluster_ids[j] for j in indices]:
		center = centers[i,:]
		if label == "dot":
			plt.scatter([center[component_i]],[center[component_j]],  marker='s', c='w', s=10)
		else:
			plt.annotate('%d' %i, xy=(center[component_i],center[component_j]), xytext=(center[component_i], center[component_j]),size=6)


	pp = PdfPages("%s/c%d_c%d_clusters%d.pdf" %(tica_dir, component_i, component_j, np.shape(centers)[0]))
	pp.savefig()
	pp.close()
	plt.clf()

def plot_pnas_vs_tics(pnas_dir, tic_dir, pnas_names, directory, scale = 7.14, refcoords_file = None):
	pnas = np.concatenate(load_file(pnas_dir))
	pnas[:,0] *= scale
	print(np.shape(pnas))
	print(len(pnas_names))
	if("ktICA" in tic_dir):
		tics = load_dataset(tic_dir)
	else:
		tics = verboseload(tic_dir)
	print(np.shape(tics))
	tics = np.concatenate(tics)
	print(np.shape(tics))
	if len(pnas_names) != np.shape(pnas)[1]:
		print("Invalid pnas names")
		return

	for i in range(0,np.shape(pnas)[1]):
		for j in range(0,np.shape(tics)[1]):
			tic = tics[:,j]
			pnas_coord = pnas[:,i]
			plt.hexbin(tic, pnas_coord, bins = 'log', mincnt=1)
			coord_name = pnas_names[i]
			tic_name = "tIC.%d" %(j+1)
			plt.xlabel(tic_name)
			plt.ylabel(coord_name)
			pp = PdfPages("%s/%s_%s_hexbin.pdf" %(directory, tic_name, coord_name))
			pp.savefig()
			pp.close()
			plt.clf()

	return

def plot_tica(transformed_data_dir, lag_time):
	transformed_data = verboseload(transformed_data_dir)
	trajs = np.concatenate(transformed_data)
	plt.hexbin(trajs[:,0], trajs[:,1], bins='log', mincnt=1)
	pp = PdfPages("/scratch/users/enf/b2ar_analysis/tica_phi_psi_chi2_t%d.pdf" %lag_time)
	pp.savefig()
	pp.close()

def main():
    fnm = 'prob3.data'
    data = md.read_data(fnm)
    D1 = data[0:8,].T
    D2 = data[8:,].T

    u1 = np.matrix((np.mean(D1[0,:]), np.mean(D1[1,:]))).T
    u2 = np.matrix((np.mean(D2[0,:]), np.mean(D2[1,:]))).T

    sigma1 = np.asmatrix(np.cov(D1, bias=1))
    sigma2 = np.asmatrix(np.cov(D1, bias=1))

    g1 = discrim_func(u1, sigma1)
    g2 = discrim_func(u2, sigma2)

    steps = 100
    x = np.linspace(-2,2,steps)
    y = np.linspace(-6,6,steps)

    X,Y = np.meshgrid(x,y)
    z = [g1(X[r,c], Y[r,c]) - g2(X[r,c], Y[r,c])
         for r in range(0,steps) for c in range(0,steps)]
    Z = np.array(z)
    px = X.ravel()
    py = Y.ravel()
    pz = Z.ravel()
    gridsize = 50
    plot = plt.subplot(111)
    plt.hexbin(px,py,C=pz, gridsize=gridsize, cmap=cm.jet, bins=None)
    cb = plt.colorbar()
    cb.set_label('g1 minus g2')
    return plot

def plot(x, y, c, filename):
    xmin=-180
    xmax = 180
    ymin = -300
    ymax = 300

    val = np.clip(c, -0.005, 1.05)
    #print val
    fig = plt.figure(1, (12,6))
    # plt.subplots_adjust(hspace=0.5)
    plt.subplot(111, axisbg='black')
    plt.hexbin(x, y, val, cmap=plt.cm.gist_heat_r)
    # plt.scatter(x, y, cmap=plt.cm.YlOrRd_r)
    plt.axis([xmax, xmin, ymin, ymax])
    plt.title("Longitude-Velocity")
    plt.xlabel('Galactic longitude (l)')
    plt.ylabel('LSR Velocity (km/s)')
    cb = plt.colorbar()
    cb.set_label(r'$e^{(-\tau)}$')

    #plt.subplot(122)
    #plt.hexbin(x, y, c, bins='log', cmap=plt.cm.YlOrRd_r)
    #plt.axis([xmin, xmax, ymin, ymax])
    #plt.title("With a log color scale")
    #cb = plt.colorbar()
    #cb.set_label('log10(N)')

    plt.savefig(filename)

    # plt.show()
    return

def execute(model, data, savepath,  lwr_data, *args, **kwargs):


    parameters = {'alpha' : np.logspace(-6,0,20), 'gamma': np.logspace(-1.5,1.5,20)}
    grid_scores = []

    for a in parameters['alpha']:
        for y in parameters['gamma']:
            model = KernelRidge(alpha= a, gamma= y, kernel='rbf')
            rms = leaveout_cv(model, data, num_runs = 200)
            #rms = kfold_cv(model, data, num_folds = 5, num_runs = 50)
            #rms = alloy_cv(model, data)
            #rms = atr2_extrap(model, data)
            #rms = lwr_extrap(model, data, lwr_data)
            grid_scores.append((a, y, rms))

    grid_scores = np.asarray(grid_scores)

    with open(savepath.replace(".png","").format("grid_scores.csv"),'w') as f:
        writer = csv.writer(f, lineterminator = '\n')
        x = ["alpha", "gamma", "rms"]
        writer.writerow(x)
        for i in grid_scores:
            writer.writerow(i)

    #Heatmap of RMS scores vs the alpha and gamma
    plt.figure(1)
    plt.hexbin(np.log10(grid_scores[:,0]), np.log10(grid_scores[:,1]), C = grid_scores[:,2], gridsize=15, cmap=cm.plasma, bins=None, vmax = 60)
    plt.xlabel('log alpha')
    plt.ylabel('log gamma')
    cb = plt.colorbar()
    cb.set_label('rms')
    plt.savefig(savepath.format("alphagammahex"), dpi=200, bbox_inches='tight')
    plt.close()