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 study_sdss_density(hemi='south'):
    grid = grid3d(hemi=hemi)
    n_data = num_sdss_data_both_catalogs(hemi, grid)
    n_rand, weight = num_sdss_rand_both_catalogs(hemi, grid)
    n_rand *= ((n_data*weight).sum() / (n_rand*weight).sum())
    delta = (n_data - n_rand) / n_rand
    delta[weight==0]=0.
    fdelta = np.fft.fftn(delta*weight)
    power = np.abs(fdelta)**2.
    ks = get_wavenumbers(delta.shape, grid.reso_mpc)
    kmag = ks[3]
    kbin = np.arange(0,0.06,0.002)
    ind = np.digitize(kmag.ravel(), kbin)
    power_ravel = power.ravel()
    power_bin = np.zeros_like(kbin)
    for i in range(len(kbin)):
        print i
        wh = np.where(ind==i)[0]
        power_bin[i] = power_ravel[wh].mean()
    #pl.clf()
    #pl.plot(kbin, power_bin)
    from cosmolopy import perturbation
    pk = perturbation.power_spectrum(kbin, 0.4, **cosmo)
    pl.clf(); pl.plot(kbin, power_bin/pk, 'b')
    pl.plot(kbin, power_bin/pk, 'bo')    
    pl.xlabel('k (1/Mpc)',fontsize=16)
    pl.ylabel('P(k) ratio, DATA/THEORY [arb. norm.]',fontsize=16)
    ipdb.set_trace()

	def transition_related_averaging_run(self, simulation_data, smoothing_kernel_width = 200, sampling_interval = [-50, 150], plot = True ):
		"""docstring for transition_related_averaging"""
		transition_occurrence_times = self.transition_occurrence_times(simulation_data = simulation_data, smoothing_kernel_width = smoothing_kernel_width)
		# make sure only valid transition_occurrence_times survive
		transition_occurrence_times = transition_occurrence_times[(transition_occurrence_times > -sampling_interval[0]) * (transition_occurrence_times < (simulation_data.shape[0] - sampling_interval[1]))]
		
		# separated into on-and off periods:
		transition_occurrence_times_separated = [transition_occurrence_times[::2], transition_occurrence_times[1::2]]
		
		mean_time_course, std_time_course = np.zeros((2, sampling_interval[1] - sampling_interval[0], 5)), np.zeros((2, sampling_interval[1] - sampling_interval[0], 5))
		if transition_occurrence_times_separated[0].shape[0] > 2:
			
			for k in [0,1]:
				averaging_interval_times = np.array([transition_occurrence_times_separated[k] + sampling_interval[0],transition_occurrence_times_separated[k] + sampling_interval[1]]).T
				interval_data = np.array([simulation_data[avit[0]:avit[1]] for avit in averaging_interval_times])
				mean_time_course[k] = interval_data.mean(axis = 0)
				std_time_course[k] = (interval_data.std(axis = 0) / np.sqrt(interval_data.shape[0]))
			
			if plot:
				f = pl.figure(figsize = (10,8))
				for i in range(simulation_data.shape[1]):
					s = f.add_subplot(simulation_data.shape[1], 1, 1 + i)
					for j in [0,1]:
						pl.plot(np.arange(mean_time_course[j].T[i].shape[0]), mean_time_course[j].T[i], ['r--','b--'][j], linewidth = 2.0 )
						pl.fill_between(np.arange(mean_time_course[j].shape[0]), mean_time_course[j].T[i] + std_time_course[j].T[i], mean_time_course[j].T[i] - std_time_course[j].T[i], ['r','b'][j], alpha = 0.2)
					s.set_title(self.variable_names[i])
				pl.draw()
			
		return (mean_time_course, std_time_course)

    def check_models(self):
        '''
        Displays a plot of the models against that taken from a
        respected website (https://www.pvlighthouse.com.au/)
        '''
        plt.figure('Intrinsic bandgap')
        t = np.linspace(1, 500)

        for author in self.available_models():

            Eg = self.update(temp=t, author=author, multiplier=1.0)
            plt.plot(t, Eg, label=author)

        test_file = os.path.join(
            os.path.dirname(os.path.realpath(__file__)),
            'Si', 'check data', 'iBg.csv')

        data = np.genfromtxt(test_file, delimiter=',', names=True)

        for temp, name in zip(data.dtype.names[0::2], data.dtype.names[1::2]):
            plt.plot(
                data[temp], data[name], '--', label=name)

        plt.xlabel('Temperature (K)')
        plt.ylabel('Intrinsic Bandgap (eV)')

        plt.legend(loc=0)
        self.update(temp=0, author=author, multiplier=1.01)

def plot_mock(mock):
    plt.clf()
    plt.plot(mock['dates'], mock['y'], marker='+', color='blue',
             label='data', markersize=9)
    plt.plot(mock['dates'], mock['y_without_seasonal'],
             color='green', alpha=0.6, linewidth=1,
             label='model without seasonal')

def compareFrequencies():
	times = generateTimes(sampleFreq, numSamples)
	signal = (80.0, 0.1)
	coherent = (60.0, 1.0)
	incoherent = (60.1, 1.0)
	highFNoise = (500.0, 0.01)
	timeData = generateTimeDomain(times, [signal, coherent, highFNoise])
	timeData2 = generateTimeDomain(times, [signal, incoherent, highFNoise])
	#timeData3 = generateTimeDomain(times, [signal, highFNoise])
	
	#timeData = generateTimeDomain(times, [(60.0, 1.0)])
	#timeData2 = generateTimeDomain(times, [(61.0, 1.0)])
	
	roi = (0, 20)
	
	freqData = map(toDb, map(dtype, map(absolute, fourier(timeData))))[roi[0]:roi[1]]
	freqData2 = map(toDb, map(dtype, map(absolute, fourier(timeData2))))[roi[0]:roi[1]]
	#freqData3 = map(toDb, map(dtype, map(absolute, fourier(timeData3))))[roi[0]:roi[1]]
	
	frequencies = generateFFTFrequencies(sampleFreq, numSamples)[roi[0]:roi[1]]
	
	#pylab.subplot(111)
	pylab.plot(frequencies, freqData)
	
	#pylab.subplot(112)
	pylab.plot(frequencies, freqData2)
	
	#pylab.plot(frequencies, freqData3)
	
	pylab.grid(True)
	pylab.show()

 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 plot_corner_posteriors(self, savefile=None, labels=["T1", "R1", "Av", "T2", "R2"]):
        '''
        Plots the corner plot of the MCMC results.
        '''
        ndim = len(self.sampler.flatchain[0,:])
        chain = self.sampler
        samples = chain.flatchain
        
        samples = samples[:,0:ndim]  
        plt.figure(figsize=(8,8))
        fig = corner.corner(samples, labels=labels[0:ndim])
        plt.title("MJD: %.2f"%self.mjd)
        name = self._get_save_path(savefile, "mcmc_posteriors")
        plt.savefig(name)
        plt.close("all")
        

        plt.figure(figsize=(8,ndim*3))
        for n in range(ndim):
            plt.subplot(ndim,1,n+1)
            chain = self.sampler.chain[:,:,n]
            nwalk, nit = chain.shape
            
            for i in np.arange(nwalk):
                plt.plot(chain[i], lw=0.1)
                plt.ylabel(labels[n])
                plt.xlabel("Iteration")
        name_walkers = self._get_save_path(savefile, "mcmc_walkers")
        plt.tight_layout()
        plt.savefig(name_walkers)
        plt.close("all")  

 def plot_fft(self,b):
     a = len(self.fullfft_dft_py_fc_0.output)
     
     for i in range(0,b):
         self.frq.append(i)
     plt.plot(self.frq,self.fullfft_dft_py_fc_0.output)
     plt.show()

def plot_locking_states(df, meta, num_joints=None):

    marker_style = dict(linestyle=':', marker='o', s=100,)
    
    def format_axes(ax):
        ax.margins(0.2)
        ax.set_axis_off()

    if num_joints is None:
        num_joints = determine_num_joints(df)

    points = np.ones(num_joints)
    
    fig, ax = plt.subplots()
    for j in range(num_joints):
        ax.text(-1.5, j, "%d" % j)
    ax.text(0, -1.5, "time")
        
    for t in df.index:
        lock_states = df.loc[t][ [ "LockingState%d" % k for k in range(num_joints) ] ].tolist()
        c = ["orange" if l else "k" for l in lock_states]
        
        ax.scatter((t+0.1) * points, range(num_joints), color=c, **marker_style)
        format_axes(ax)
        
    ax.set_title('Locking state evolution')
    ax.set_xlabel("t")
    
    plt.plot()

def check_isometry(G, chart, nseeds=100, verbose = 0):
    """
    A simple check of the Isometry:
    look whether the output distance match the intput distances
    for nseeds points
    
    Returns
    -------
    a scaling factor between the proposed and the true metrics
    """
    nseeds = np.minimum(nseeds, G.V)
    aux = np.argsort(nr.rand(nseeds))
    seeds =  aux[:nseeds]
    dY = Euclidian_distance(chart[seeds],chart)
    dx = G.floyd(seeds)

    dY = np.reshape(dY,np.size(dY))
    dx = np.reshape(dx,np.size(dx))

    if verbose:
        import matplotlib.pylab as mp
        mp.figure()
        mp.plot(dx,dY,'.')
        mp.show()

    scale = np.dot(dx,dY)/np.dot(dx,dx)
    return scale

def fdr(p_values=None, verbose=0):
    """Returns the FDR associated with each p value

    Parameters
    -----------
    p_values : ndarray of shape (n)
        The samples p-value

    Returns
    -------
    q : array of shape(n)
        The corresponding fdr values
    """
    p_values = check_p_values(p_values)
    n_samples = p_values.size
    order = p_values.argsort()
    sp_values = p_values[order]

    # compute q while in ascending order
    q = np.minimum(1, n_samples * sp_values / np.arange(1, n_samples + 1))
    for i in range(n_samples - 1, 0, - 1):
        q[i - 1] = min(q[i], q[i - 1])

    # reorder the results
    inverse_order = np.arange(n_samples)
    inverse_order[order] = np.arange(n_samples)
    q = q[inverse_order]

    if verbose:
        import matplotlib.pylab as mp
        mp.figure()
        mp.xlabel('Input p-value')
        mp.plot(p_values, q, '.')
        mp.ylabel('Associated fdr')
    return q

def plot_q(model='cem', r_min=0.0, r_max=6371.0, dr=1.0):
    """
    Plot a radiallysymmetric Q model.

    plot_q(model='cem', r_min=0.0, r_max=6371.0, dr=1.0):

    r_min=minimum radius [km], r_max=maximum radius [km], dr=radius
    increment [km]

    Currently available models (model): cem, prem, ql6
    """
    import matplotlib.pylab as plt

    r = np.arange(r_min, r_max + dr, dr)
    q = np.zeros(len(r))

    for k in range(len(r)):

        if model == 'cem':
            q[k] = q_cem(r[k])
        elif model == 'ql6':
            q[k] = q_ql6(r[k])
        elif model == 'prem':
            q[k] = q_prem(r[k])

    plt.plot(r, q, 'k')
    plt.xlim((0.0, r_max))
    plt.xlabel('radius [km]')
    plt.ylabel('Q')
    plt.show()

    def check_models(self):
        plt.figure('Bandgap narrowing')
        Na = np.logspace(12, 20)
        Nd = 0.
        dn = 1e14
        temp = 300.

        for author in self.available_models():
            BGN = self.update(Na=Na, Nd=Nd, nxc=dn,
                              author=author,
                              temp=temp)

            if not np.all(BGN == 0):
                plt.plot(Na, BGN, label=author)

        test_file = os.path.join(
            os.path.dirname(os.path.realpath(__file__)),
            'Si', 'check data', 'Bgn.csv')

        data = np.genfromtxt(test_file, delimiter=',', names=True)

        for name in data.dtype.names[1:]:
            plt.plot(
                data['N'], data[name], 'r--',
                label='PV-lighthouse\'s: ' + name)

        plt.semilogx()
        plt.xlabel('Doping (cm$^{-3}$)')
        plt.ylabel('Bandgap narrowing (K)')

        plt.legend(loc=0)

def existe_croche_blanche_mesure(img,img2,liste,ecart):
	for elt in liste:
		#Si on a une noire en haut ou (exclusif) en bas
		if (not(elt[3]) and elt[4]) or (elt[3] and not(elt[4])):
			if elt[3]:
				elt.append(existe_croche_haut(img,ecart,elt[0],elt[2]))
			else:
				elt.append(existe_croche_bas(img,ecart,elt[1],elt[2]))
			#on regarde s'il y a d'autres croches
			liste = existe_autre_croche(img,liste,ecart)
			elt.extend([False,False])
			
		#s'il n'y a pas de noire
		elif (not(elt[3]) and not(elt[4])):
			#on met le nombre de croches à zéro
			elt.append(0)
			elt.append(existe_note(img2,ecart,elt[1],elt[2],pc_blan,'magenta'))
			elt.append(existe_note(img2,ecart,elt[0],elt[2],pc_blan,'magenta'))
			
			#c'est une barre de mesure (ni noire, ni blanche)
			if (not(elt[6]) and not(elt[7])):
				#elt.extend('m')
				x = [elt[2],elt[2]]
				y = [elt[0],elt[1]]
				plt.plot(x,y,'b')
	return liste

 def sanity_RotBroadProfileExample(self):
   """
     Example of rotational broadening.
   """
   import numpy as np
   import matplotlib.pylab as plt
   from PyAstronomy import modelSuite as ms
   
   # Get an instance of the model ...
   x = ms.RotBroadProfile()
   # ... and define some starting value
   x["xmax"] = 60.0
   x["A"] = 1.0
   x["eps"] = 0.8
   x["off"] = 0.0
   
   # Define a radial velocity axis
   vv = np.linspace(-90.,90.,200)
   
   # Construct some "data" and ...
   data = x.evaluate(vv)
   # ... add noise
   data += np.random.normal(0.0, 1e-3, data.size)
   
   # Fit the model using A, xmax, and eps as free
   # parameters ...
   x.thaw(["A", "xmax", "eps"])
   x.fit(vv, data)
   # ... and show the resulting parameter values.
   x.parameterSummary()
   
   # Plot the data and the model
   plt.plot(vv, data, 'bp')
   plt.plot(vv, x.model, 'r--')

def histograma(hist):
    
    hist=hist.histogram(255)
##    hist.save("hola4Hist.txt")
    pylab.plot(hist)
    pylab.draw()
    pylab.pause(0.0001)

	def plot_cell(self,cell_number=0,label='insert_label'):

		current_cell = self.cell_list[cell_number]
		temp = current_cell.temp
		cd_signal = current_cell.cd_signal
		cd_calc = current_cell.cd_calc()
		
		ax = pylab.gca()

		pylab.plot(temp,cd_signal,'o',color='black')
                pylab.plot(temp,cd_calc,color='black')
		pylab.xlabel(r'Temperature ($^{\circ}$C)')
		pylab.ylabel('mdeg')
		pylab.ylim([-25,-4])
		dH = numpy.round(current_cell.dH, decimals=1)
		Tm = numpy.round(current_cell.Tm-273.15, decimals=1)
		nf = current_cell.nf
		nu = current_cell.nu
		textstr_dH = '${\Delta}H_{m}$ = %.1f kcal/mol' %dH
		textstr_Tm ='$T_{m}$ = %.1f $^{\circ}$C' %Tm
		textstr_nf ='$N_{folded}$ = %d' %nf
		textstr_nu ='$N_{unfolded}$ = %d'%nu
		ax.text(8,-6,textstr_dH, fontsize=16,ha='left',va='top')
		ax.text(8,-7.5,textstr_Tm, fontsize=16,ha='left',va='top')
		ax.text(8,-9,textstr_nf, fontsize=16,ha='left',va='top')
		ax.text(8,-10.5,textstr_nu, fontsize=16,ha='left',va='top')
		pylab.title(label)		
		pylab.show()

		return

 def plot_values(self, TITLE, SAVE):
     plot(self.list_of_densities, self.list_of_pressures)
     title(TITLE)
     xlabel("Densities")
     ylabel("Pressure")
     savefig(SAVE)
     show()

def fancy_dendrogram(*args, **kwargs):
    '''
    Source: https://joernhees.de/blog/2015/08/26/scipy-hierarchical-clustering-and-dendrogram-tutorial/
    '''
    from scipy.cluster import hierarchy
    import matplotlib.pylab as plt
    
    max_d = kwargs.pop('max_d', None)
    if max_d and 'color_threshold' not in kwargs:
        kwargs['color_threshold'] = max_d
    annotate_above = kwargs.pop('annotate_above', 0)

    ddata = hierarchy.dendrogram(*args, **kwargs)

    if not kwargs.get('no_plot', False):
        plt.title('Hierarchical Clustering Dendrogram (truncated)')
        plt.xlabel('sample index or (cluster size)')
        plt.ylabel('distance')
        for i, d, c in zip(ddata['icoord'], ddata['dcoord'], ddata['color_list']):
            x = 0.5 * sum(i[1:3])
            y = d[1]
            if y > annotate_above:
                plt.plot(x, y, 'o', c=c)
                plt.annotate("%.3g" % y, (x, y), xytext=(0, -5),
                             textcoords='offset points',
                             va='top', ha='center')
        if max_d:
            plt.axhline(y=max_d, c='k')
    return ddata