def generateToy():

  print 'loading values'
  if not os.path.isfile('values2.p'):
    z_data = np.loadtxt('values2.dat')
    pkl.dump( z_data, open( 'values2.p', "wb" ),pkl.HIGHEST_PROTOCOL )
  else:
    z_data = pkl.load(open('values2.p',"rb"))
  print 'loaded'

  #x = np.random.normal(size=1000)
  z_data_subset = z_data[0:20000]
  plot_range = [50,400]
  print 'max',max(z_data_subset),'min',min(z_data_subset)
  plt.yscale('log', nonposy='clip')
  plt.axes().set_ylim(0.0000001,0.17)
  hist(z_data_subset,range=plot_range,bins=100,normed=1,histtype='stepfilled',
      color=['lightgrey'], label=['100 bins'])
  #hist(z_data_subset,range=plot_range,bins='knuth',normed=1,histtype='step',linewidth=1.5,
  #    color=['navy'], label=['knuth'])
  hist(z_data_subset,range=plot_range,bins='blocks',normed=1,histtype='step',linewidth=2.0,
      color=['crimson'], label=['b blocks'])
  plt.legend()
  #plt.yscale('log', nonposy='clip')
  #plt.axes().set_ylim(0.0000001,0.17)
  plt.xlabel(r'$m_{\ell\ell}$ (GeV)')
  plt.ylabel('A.U.')
  plt.title(r'Z$\to\mu\mu$ Data')
  plt.savefig('z_data_hist_comp.png')
  plt.show()

    def triangle(self):
        assert self.sample_invoked == True, \
            'Must sample first! Use sample(iter, burn)'

        fig = plt.figure(figsize=(4, 4))
        ax = fig.add_subplot(111)

        hist(self.trace[:].flatten(), bins='knuth', normed=True,
             histtype='step', color='k', ax=ax)
        plt.xlabel('$b$')

def histogram(totalDF, column, axis, binN=10, incl_untrunc=True):
	trunc = totalDF[(totalDF.untruncated == False) & (totalDF.h2 > 0) & (totalDF.M1 < 50)][column].dropna()
	if not incl_untrunc:
		untrunc = trunc
	else:
		untrunc = totalDF[(totalDF.untruncated == True) & (totalDF.h2 > 0) & (totalDF.M1 < 50)][column].dropna()
	bins = get_bins(binN, untrunc, trunc)
	wuntrunc, wtrunc = [1./len(untrunc)]*len(untrunc), [1./len(trunc)]*len(trunc)
	hist(trunc, bins, ax=axis, alpha=0.5, label='Truncated', weights=wtrunc)
	if incl_untrunc:
		hist(untrunc, bins, ax=axis, alpha=0.5, label='Untruncated', weights=wuntrunc)

def bin_edges_f(bin_method, mag_col_cl):
    '''
    Obtain bin edges for each photometric dimension using the cluster region
    diagram.
    '''
    bin_edges = []
    if bin_method in ['sturges', 'sqrt']:
        if bin_method == 'sturges':
            b_num = 1 + np.log2(len(mag_col_cl[0]))
        else:
            b_num = np.sqrt(len(mag_col_cl[0]))

        for mag_col in mag_col_cl:
            bin_edges.append(np.histogram(mag_col, bins=b_num)[1])

    elif bin_method == 'bb':
        # Based on Bonatto & Bica (2007) 377, 3, 1301-1323. Fixed bin width
        # of 0.25 for colors and 0.5 for magnitudes.
        b_num = [(max(mag_col_cl[0]) - min(mag_col_cl[0])) / 0.25,
                 (max(mag_col_cl[1]) - min(mag_col_cl[1])) / 0.5]

        for i, mag_col in enumerate(mag_col_cl):
            bin_edges.append(np.histogram(mag_col, bins=b_num[i])[1])

    else:
        for mag_col in mag_col_cl:
            bin_edges.append(hist(mag_col, bins=bin_method)[1])

    return bin_edges

def get_index(df):
    """
    A bit of a wasteful, hackish way to get an index to use for the various
    lags.
    """
    _, idx, _ = hist(df.diff(1).loc[30], bins='scott', alpha=.35)
    return idx

	def histo(self):
		#------------------------------------------------------------
		# First figure: show normal histogram binning
		fig = plt.figure(figsize=(10, 4))
		fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
	
		ax1 = fig.add_subplot(121)
		ax1.hist(self.entropy, bins=15, histtype='stepfilled', alpha=0.2, normed=True)
		ax1.set_xlabel('entropy bins=15')
		ax1.set_ylabel('Count(t)')
	
		ax2 = fig.add_subplot(122)
		ax2.hist(self.entropy, bins=200, histtype='stepfilled', alpha=0.2, normed=True)
		ax2.set_xlabel('entropy bins=200')
		ax2.set_ylabel('Count(t)')
	
		#------------------------------------------------------------
		# Second & Third figure: Knuth bins & Bayesian Blocks
		fig = plt.figure(figsize=(10, 4))
		fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)
		
		for bins, title, subplot in zip(['knuth', 'blocks'],
		                                ["Knuth's rule-fixed bin-width", 'Bayesian blocks variable width'],
		                                [121, 122]):
		    ax = fig.add_subplot(subplot)
		
		    # plot a standard histogram in the background, with alpha transparency
		    hist(self.entropy, bins=200, histtype='stepfilled',
		         alpha=0.2, normed=True, label='standard histogram')
		
		    # plot an adaptive-width histogram on top
		    hist(self.entropy, bins='blocks', ax=ax, color='black',
		         histtype='step', normed=True, label=title)
		
		    ax.legend(prop=dict(size=12))
		    ax.set_xlabel('entropy bins')
		    ax.set_ylabel('C(t)')
		
		plt.show()

    def auto_discretize(self, num_data, method, range_min_max):

        """
        Perform automatic discretization of a selected feature; a method
        (bayesian blocks, scott method or fixed bin number) along the desired data range is passed to a special version of hist which gives cutpoints for discretization and returns the "categorized" version of the original data
        """
        hist_data = hist(num_data, bins=method, range=range_min_max)
        plt.close("all")
        leng = len(hist_data[1])
        # fix cutoff to make sure outliers are properly categorized as well if necessary
        hist_data[1][leng - 1] = num_data.max()
        # hist_data[1][0] = num_data.min()
        # automatically assign category labels of '1','2',etc
        cat_data = pandas.cut(num_data, hist_data[1], labels=range(1, leng), include_lowest="TRUE")
        return pandas.Series(cat_data).astype(str)

def noise(fname, x0 = 100, y0 = 100, maxrad = 30):
    from astroML.plotting import hist
    hdulist = pf.open(fname)
    im = hdulist[0].data
    #print np.mean(im), np.min(im), np.max(im)
    #print im[95:105, 95:105]
    # x0, y0 = 100, 100
    xi, yi = np.indices(im.shape)
    R = np.sqrt( (yi - int(y0))**2. + (xi - int(x0))**2. )
    phot_a = np.zeros(maxrad + 1)
    phot_a[0] = 0
    
    bmasked = im * ((R > maxrad) * (R < maxrad + 20.))
    bdata = bmasked.flatten()
    #print bdata[bdata != 0.]
    #print len(bdata[bdata != 0.])
    #print len(bdata)
    
    plt.subplot(3, 1, 1)
    hist(bdata[bdata != 0.], bins = 'blocks')
    plt.xlabel('Flux')
    plt.ylabel('(Bayesian Blocks)')
    plt.title('Noise')
    #plt.show()
    
    plt.subplot(3, 1, 2)
    hist(bdata[bdata != 0.], bins = 50)
    plt.xlabel('Flux')
    plt.ylabel('(50 bins)')
    #plt.title('Noise (50 bins)')
    #plt.show()
    
    plt.subplot(3, 1, 3)
    hist(bdata[bdata != 0.], bins = 'knuth')
    plt.xlabel('Flux')
    plt.ylabel('(Knuth\'s Rule)')
    #plt.title('Noise (Knuth\'s Rule)')
    plt.show()
    
    A2, crit, sig = anderson(bdata[bdata != 0.], dist = 'norm')
    print 'A-D Statistic:', A2
    print ' CVs \t  Sig.'
    print np.vstack((crit, sig)).T

    normality = normaltest(bdata[bdata != 0.])
    print 'Normality:', normality

    skewness = skewtest(bdata[bdata != 0.])
    print 'Skewness:', skewness

    kurtosis = kurtosistest(bdata[bdata != 0.])
    print 'Kurtosis:', kurtosis

    print 'Mean:', np.mean(bdata[bdata != 0.])
    print 'Median:', np.median(bdata[bdata != 0.])

def plot_labeled_histogram(style, data, name,
                           x, pdf_true, ax=None,
                           hide_x=False,
                           hide_y=False):
    if ax is not None:
        ax = plt.axes(ax)

    counts, bins, patches = hist(data, bins=style, ax=ax,
                                 color='k', histtype='step', normed=True)
    ax.text(0.99, 0.95, '%s:\n%i bins' % (name, len(counts)),
            transform=ax.transAxes,
            ha='right', va='top', fontsize=12)

    ax.fill(x, pdf_true, '-', color='#CCCCCC', zorder=0)

    if hide_x:
        ax.xaxis.set_major_formatter(plt.NullFormatter())
    if hide_y:
        ax.yaxis.set_major_formatter(plt.NullFormatter())

    ax.set_xlim(-5, 5)

    return ax

def plot_wage_change_dist(df, pi, lambda_, nperiods=4, log=True, figkwargs=None,
                          axkwargs=None):
    """
    Make and save the figure for the distribution of wage changes.

    figkwargs is a dict passed to the fig constructor.
    axkwargs is a dict passed to the axes constructor.
    """
    idx = get_index(df)
    SOME_SS = 30  # just some period in the steady state.
    diffs = range(1, nperiods + 1)

    if log:
        df = np.log(df)
        strlog = ', (log scale),'  # see title formatting
    else:
        strlog = ''

    t = pd.concat([df.diff(x).iloc[SOME_SS] for x in diffs],
                  axis=1, keys=diffs)
    _figkwargs = {'figsize': (13, 8)}  # Leading _ is internal.
    if figkwargs is not None:
        _figkwargs.update(figkwargs)

    _axkwargs = {}
    if axkwargs is not None:
        _axkwargs.update(axkwargs)

    fig, ax = plt.subplots(**_figkwargs)
    cts, idx, other = hist(t.values, histtype='bar', bins=idx,
                           label=['lag={}'.format(i) for i in diffs],
                           ax=ax, normed=True, **_axkwargs)

    ax.set_title('Across Periods{0} $\pi={1:.3f}$, $\lambda={2:.3f}$'.format(
        strlog, pi, lambda_))
    ax.legend()
    return fig, ax

def plot_redshifts():
    
    gz2main_fitsfile = '/Users/willettk/Astronomy/Research/GalaxyZoo/fits/gz2main_table_sample.fits'
    hdulist = pyfits.open(gz2main_fitsfile)
    gz2main = hdulist[1].data
    hdulist.close()

    redshift_all = gz2main['redshift']
    redshift_finite = redshift_all[np.isfinite(redshift_all)]

    fig = plt.figure()
    ax = fig.add_subplot(111)
#   hist(redshift_finite, bins='blocks', ax=ax, histtype='stepfilled', color='r', ec='r', normed=True)
    hist(redshift_finite, bins='scotts', ax=ax, histtype='step', color='b',normed=True)
    hist(redshift_finite, bins='freedman', ax=ax, histtype='step', color='y',normed=True)
    hist(redshift_finite, bins='knuth', ax=ax, histtype='step', color='y',normed=True)
    ax.set_xlabel('Redshift')
    ax.set_ylabel('Frequency')
    
    plt.show()

    return None

import mpl_toolkits

pl.rcParams['font.size'] = 20


nh2 = r'$\log(n(H_2))$ [cm$^{-3}$]'

dens = fits.getdata('W51_H2CO11_to_22_logdensity_supersampled.fits')
cube = pyspeckit.Cube('W51_H2CO11_to_22_logdensity_supersampled.fits')
densOK = dens==dens

pl.figure(2)
pl.clf()
ax = pl.subplot(1,2,1)
densp = dens[densOK]
counts,bins,patches = ampl.hist(densp, bins=100, log=True, histtype='step', linewidth=2, alpha=0.8, color='k')
ylim = ax.get_ylim()

sp = pyspeckit.Spectrum(xarr=(bins[1:]+bins[:-1])/2.,data=counts)
sp.specfit(guesses= [660.23122694399035,
                     3.1516848752486522,
                     0.33836811902343894,
                     396.62714060001434,
                     2.5539176548294318,
                     0.32129608858734149,
                     199.13259679527025,
                     3.730112763513838,
                     0.4073913996012487])

def ntuples(lst, n):
    return zip(*[lst[i::n]+lst[:i:n] for i in range(n)])

    t = np.linspace(-10, 30, 1000)

    # Compute density with KDE
    kde = KDE('gaussian', h=0.1).fit(xN[:, None])
    dens_kde = kde.eval(t[:, None]) / N

    # Compute density with Bayesian nearest neighbors
    nbrs = KNeighborsDensity('bayesian', n_neighbors=k).fit(xN[:, None])
    dens_nbrs = nbrs.eval(t[:, None]) / N

    # plot the results
    #ax.plot(t, true_pdf(t), ':', color='black', zorder=3,
    #        label="Generating Distribution")
    ax.plot(xN, -0.005 * np.ones(len(xN)), '|k', lw=1.5)
    hist(xN, bins='blocks', ax=ax, normed=True, zorder=1,
         histtype='stepfilled', lw=1.5, color='k', alpha=0.2,
         label="Bayesian Blocks")
    ax.plot(t, dens_nbrs, '-', lw=2, color='gray', zorder=2,
            label="Nearest Neighbors (k=%i)" % k)
    ax.plot(t, dens_kde, '-', color='black', zorder=3,
            label="Kernel Density (h=0.1)")

    # label the plot
    ax.text(0.02, 0.95, "%i points" % N, ha='left', va='top',
            transform=ax.transAxes)
    ax.set_ylabel('$p(x)$')
    ax.legend(loc='upper right', prop=dict(size=12))

    if subplot == 212:
        ax.set_xlabel('$x$')

from paths import dpath,fpath

pl.rcParams['font.size'] = 20

h2co11 = fits.getdata(dpath('W51_H2CO11_taucube_supersampled.fits'))
h2co22 = fits.getdata(dpath('W51_H2CO22_pyproc_taucube_lores_supersampled.fits'))

ratio = fits.getdata(dpath('W51_H2CO11_to_22_tau_ratio_supersampled_neighbors.fits'))

pl.close(1)
pl.figure(1, figsize=(10,10))
pl.clf()

ax1 = pl.subplot(3,1,1)
oneone = h2co11[h2co11==h2co11]
counts, bins, patches = ampl.hist(oneone, bins=100, log=True, histtype='step',
                                  linewidth=2, alpha=0.8, color='k')
ylim = ax1.get_ylim()
med, mad = np.median(oneone),MAD(oneone)
pl.plot(bins,counts.max()*np.exp(-(bins-med)**2/(2*mad**2)),'r--')
ax1.set_ylim(*ylim)
ax1.set_xlabel("$\\tau_{obs}($H$_2$CO 1-1$)$", labelpad=10)
ax1.set_ylabel("$N($voxels$)$")

ax2 = pl.subplot(3,1,2)
twotwo = h2co22[h2co22==h2co22]
counts, bins, patches = ampl.hist(twotwo, bins=100, log=True, histtype='step', linewidth=2, alpha=0.8, color='k')
ylim = ax2.get_ylim()
med, mad = np.median(twotwo),MAD(twotwo)
pl.plot(bins,counts.max()*np.exp(-(bins-med)**2/(2*mad**2)),'r--')
ax2.set_ylim(*ylim)
ax2.set_xlabel("$\\tau_{obs}($H$_2$CO 2-2$)$", labelpad=10)

imX = np.empty((len(image_data), 2), dtype=np.float64)
imX[:, 0] = image_data['ra']
imX[:, 1] = image_data['dec']

# get standard stars
standards_data = fetch_sdss_S82standards()
stX = np.empty((len(standards_data), 2), dtype=np.float64)
stX[:, 0] = standards_data['RA']
stX[:, 1] = standards_data['DEC']

# crossmatch catalogs
max_radius = 1. / 3600  # 1 arcsec
dist, ind = crossmatch_angular(imX, stX, max_radius)
match = ~np.isinf(dist)

dist_match = dist[match]
dist_match *= 3600

ax = plt.axes()
hist(dist_match, bins='knuth', ax=ax,
     histtype='stepfilled', ec='k', fc='#AAAAAA')
ax.set_xlabel('radius of match (arcsec)')
ax.set_ylabel('N(r, r+dr)')
ax.text(0.95, 0.95,
        "Total objects: %i\nNumber with match: %i" % (imX.shape[0],
                                                      np.sum(match)),
        ha='right', va='top', transform=ax.transAxes)
ax.set_xlim(0, 0.2)

plt.show()

def bayes_block(x_data, filename, format, x_label = '', title = '',\
log_x = False, log_y = False):
	'''
	Description
	    This function takes the given data and produces a Bayesian Block 
	    histogram of it. The given axis label and title are applied, and then
	    the histogram is saved using the given filename and format.
	
	Required Input
	    x_data: The data array to be graphed. Numpy array or list of floats.
	            This array is flattened to one dimension before creating 
	            the histogram.
	    filename: The filename (including extension) to use when saving the
		      image. Provide as a string.
	    format: The format (e.g. png, jpeg) in which to save the image. This
	            is a string.
	    x_label: String specifying the x-axis label.
	    title: String specifying the title of the graph.
	    log_x: If this is True, then logarithmic binning is used for the
	           histogram, and the x-axis of the saved image is logarithmic.
	           If this is False (default) then linear binning is used.
	    log_y: If this is True, then a logarithmic scale is used for the
	           y-axis of the histogram. If this is False (default), then
	           a linear scale is used.
	
	Output
	    A histogram is automatically saved using the given data and labels,
	    in the specified format. 1 is returned if the code performs to
	    completion.
	'''
	
	# First make a figure object with matplotlib (default size)
	fig = plt.figure()
	# Create an axis object to go with this figure
	ax = fig.add_subplot(111)
	
	# Check to see if the x-axis of the histogram needs to be logarithmic
	if log_x == True:
		# Set the x-axis scale of the histogram to be logarithmic
		ax.set_xscale('log')

    # Make a histogram of the given data, with the specified number of 
    # bins. Note that the data array is flattened to one dimension.
    # Do we need to normalise to account for the bin sizes being different?
	aML.hist(x_data.flatten(), bins = 'blocks', normed = True, log = log_y)    
	
	# Add the specified x-axis label to the plot
	plt.xlabel(x_label)
	# Add a y-axis label to the plot
	plt.ylabel('Counts')
	# Add the specified title to the plot
	plt.title(title)
	
	# Save the figure using the title given by the user
	plt.savefig(filename, format = format)
	
	# Print a message to the screen saying that the image was created
	print filename + ' created successfully.\n'
	
	# Close the figure so that it does not take up memory
	plt.close()
	
	# Now that the graph has been produced, return 1
	return 1

ax1.set_ylabel('P(entropy)')

ax2 = fig.add_subplot(122)
ax2.hist(ent, bins=1000, histtype='stepfilled', alpha=0.2, normed=True)
ax2.set_xlabel('entropy')
ax2.set_ylabel('P(entropy)')

#------------------------------------------------------------
# Second & Third figure: Knuth bins & Bayesian Blocks
fig = plt.figure(figsize=(10, 4))
fig.subplots_adjust(left=0.1, right=0.95, bottom=0.15)

for bins, title, subplot in zip(['knuth', 'blocks'],
                                ["Knuth's rule", 'Bayesian blocks'],
                                [121, 122]):
    ax = fig.add_subplot(subplot)

    # plot a standard histogram in the background, with alpha transparency
    hist(ent, bins=200, histtype='stepfilled',
         alpha=0.2, normed=True, label='standard histogram')

    # plot an adaptive-width histogram on top
    hist(ent, bins=bins, ax=ax, color='black',
         histtype='step', normed=True, label=title)

    ax.legend(prop=dict(size=12))
    ax.set_xlabel('WTS')
    ax.set_ylabel('P(WTS)')

plt.show()

# print j


# peakabsmagvalue = [peakabsmagvalueb,peakabsmagvaluec]
# print peakabsmagvalueb
# hist(peakabsmagvalueb, bins = 'knuth', label = str(ib) + ' Ib datapoints', color = 'blue', histtype='stepfilled', alpha=0.2)#, stacked=True)


# hist(peakabsmagvaluec, bins = 'knuth', label = str(ic) + ' Ic datapoints', color ='green', histtype='stepfilled', alpha=0.2, des)#, stacked=True)


# plotting best fit gaussian
plt.subplot(221)

result = hist(
    peakabsmagvalueb, bins="knuth", label=str(ib) + " Ib datapoints", color="blue", histtype="stepfilled", alpha=0.2
)  # , stacked=True)

mean = np.mean(peakabsmagvalueb)
variance = np.var(peakabsmagvalueb)
sigma = np.sqrt(variance)
x = np.linspace(-23, -13, 100)
dx = result[1][1] - result[1][0]
scale = len(peakabsmagvalueb) * dx

plt.plot(
    x,
    mlab.normpdf(x, mean, sigma) * scale,
    label="Best fit, mean: " + str(round(mean, 3)) + " sigma: " + str(round(sigma, 3)),
)

def plot_hist(x, filename):
  fig = plt.figure(figsize=(10, 10))
  hist(x, bins='scott')
  plt.savefig(filename)
  plt.close()	

'GRPCZ', 'FC', 'LOGMH', 'DEN1MPC'.
"""

data = np.genfromtxt("ECO_dr1_subset.csv", delimiter=",", dtype=None, names=True)
name = data['NAME']
logmstar = data ['LOGMSTAR']
urcolor = data['MODELU_RCORR']
cz = data['CZ']
goodur = (urcolor > -99) & (logmstar > 10.)

colors=urcolor[goodur]

# First plot histograms of u-r color with different bin width "rules"
plt.figure(1)
plt.clf()
hist(colors,bins='freedman',label='freedman',normed=1,histtype='stepfilled',color='green',alpha=0.5)
hist(colors,bins='scott',label='scott',normed=1,histtype='step',color='purple',alpha=0.5,hatch='///')
# note the different format used below so as to save the bin info for Knuth's rule
n0, bins0, patches0 = hist(colors,bins='knuth',label='knuth',normed=1,histtype='stepfilled',color='blue',alpha=0.25)
plt.xlim(0,3)
plt.xlabel("u-r color (mag)")
plt.title("Galaxy Color Distribution")
plt.legend(loc="best")
# As in Fig. 5.20 (p. 227), Scott's rule makes broader bins.

# Now give Kernel Density Estimation a try. KDE is shown in Ivezic+ Fig. 6.1 
# but we'll use this newer version: sklearn.neighbors.KernelDensity -- see
# http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html

bw = 0.5*(bins0[2]-bins0[1]) 
# initially using 0.5*Knuth binsize from above as bandwidth; should test other values