def setDiffsPlot(CF,d0,ySym = True):
    CFtext = [str(j)+',' for j in CF]
    bText = [CFtext[0],r'...,$a_i$+c,...',CFtext[-1]]
    CFtext = ''.join(CFtext)
    bText= ''.join(bText)
    CFtext = '[' + CFtext[:-1] + ']'
    bText = '[' + bText[:-1] + ']'
    print(CFtext)

    plt.ylabel(r'$d^{crit}_b - d^{crit}_a$',fontsize=20)
    plt.xlabel(r'Element changed',fontsize=20)
    xmin, xmax, ymin, ymax = plt.axis()
    if ySym:
        plt.yscale('symlog',linthreshy=1e-15)
        yLoc = [y*ymax for y in [.1, .01, .001]]
    else:
        yLoc = [y*(ymax-ymin)+ymin for y in [0.95, 0.85, 0.75]]
    plt.plot([0,xmax],[0,0],'k--',label='_')
    plt.text((xmax-xmin)*0.15,yLoc[0],r'$a = [a_i] =$'+CFtext,fontsize=15)
    plt.text((xmax-xmin)*0.15,yLoc[1],r'$b_i = $'+bText,fontsize=15)
    plt.text((xmax-xmin)*0.15,yLoc[2],r'$d_a^{crit} = $'+str(float(d0)),fontsize=15)
    plt.legend(loc='best')
    plt.xscale('symlog',linthreshx=1e-14)
    # plt.yscale('log')
    plt.show()

 def plot_probability_calibration_curves(self):
 
     """ Compute true and predicted probabilities for a calibration plot 
         fraction_of_positives - The true probability in each bin (fraction of positives).
         mean_predicted_value - The mean predicted probability in each bin.
     """
     
     fig = plt.figure()
     ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
     ax2 = plt.subplot2grid((3, 1), (2, 0), rowspan=2)
     
     ax1.set_ylabel("Fraction of positives")
     ax1.set_ylim([-0.05, 1.05])
     ax1.legend(loc="lower right")
     ax1.set_title('Calibration plots  (reliability curve) ' + self.description)
 
     ax2.set_xlabel("Mean predicted value")
     ax2.set_ylabel("Count")
     ax2.legend(loc="upper center", ncol=2)
     
     clf_score = brier_score_loss(self.y_true, self.y_pred, pos_label=1)
     
     
     fraction_of_positives, mean_predicted_value = calibration_curve(self.y_true, self.y_pred, n_bins=50)
     
     ax1.plot(mean_predicted_value, fraction_of_positives, "s-", color="#660066",  alpha = 0.6, label="%s (%1.3f)" % (self.description, clf_score))
     ax2.hist(self.y_pred, range=(0, 1), bins=50, color="#660066", linewidth=2.0 , alpha = 0.6, label="%s (%1.3f)" % (self.description, clf_score), histtype="step", lw=2)
     plt.yscale('log')
     return

def plot_max_avg_Rho_lev012(lev0,lev1,lev2, ic, spath):
    '''lev -- TimeProfQs() object, whose lev.convert() attribute has been called'''
    #avgRho
    Tratio = lev0.t / ic.tCr
    fig, ax = plt.subplots()
    #initial
    plt.hlines(ic.rho0, Tratio.min(),Tratio.max(), colors="magenta", linestyles='dashed',label="initial")
    #lev0
    plt.plot(Tratio,lev0.maxRho,ls="-",c="black",lw=2.,label="max Lev0")
    plt.plot(Tratio,lev0.minRho,ls="-",c="green",lw=2.,label="min Lev0")
    plt.plot(Tratio,lev0.avgRho,ls="-",c="blue",lw=2.,label="avg Lev0")
    plt.plot(Tratio,lev0.avgRho_HiPa,ls="-",c="red",lw=2.,label=r"avg ($\rho > \rho_0$)")
    #lev1
    plt.plot(Tratio,lev1.maxRho,ls="--",c="black",lw=2.,label="max Lev1")
    plt.plot(Tratio,lev1.minRho,ls="--",c="green",lw=2.,label="min Lev1")
    plt.plot(Tratio,lev1.avgRho,ls="--",c="blue",lw=2.,label="avg Lev1")
    plt.plot(Tratio,lev1.avgRho_HiPa,ls="--",c="red",lw=2.,label=r"avg Lev1 ($\rho > \rho_0$)")
    #lev2
    plt.plot(Tratio,lev2.maxRho,ls=":",c="black",lw=2.,label="max Lev2")
    plt.plot(Tratio,lev2.minRho,ls=":",c="green",lw=2.,label="min Lev2")
    plt.plot(Tratio,lev2.avgRho,ls=":",c="blue",lw=2.,label="avg Lev2")
    plt.plot(Tratio,lev0.avgRho_HiPa,ls=":",c="red",lw=2.,label=r"avg Lev1 ($\rho > \rho_0$)")
    #finish
    plt.yscale("log")
    plt.xlabel("t / t_cross")
    plt.ylabel("Densities [g/cm^3]")
    plt.title(r"Max & Avg Densities, Lev0")
    py.legend(loc=4, fontsize="small")
    name = spath+"max_avg_Rho_Lev012.pdf"
    plt.savefig(name,format="pdf")
    plt.close()

def bar_graph_dict(dict_to_plot, plot_title="", xlab="", ylab="", log_scale=False, col="#71cce6", sort_key_list=None,
                   min_count=1):
    """
    Plots a bar graph of the provided dictionary.
    Params:
    dict_to_plot (dict): should have the format {'label': count}
    plot_title (str), xlab (str), ylab (str), log_scale (bool): fairly self-explanatory plot customization
    col (str): colour for the bars
    sort_key_list (list): the keys of the dictionary are assumed to match based its first item and reordered as such
    min_count (int): do not plot items with less than this in the count
    """
    # Sort dictionary & convert to list using custom keys if needed
    if not sort_key_list:
        list_to_plot = sorted(dict_to_plot.items())
    else:
        list_to_plot = sorted(dict_to_plot.items(), key=lambda x: sort_key_list.index(x[0]))

    # Remove list items with less than min_count
    if min_count != 1:
        list_to_plot = [dd for dd in list_to_plot if dd[1] >= min_count]

    # Bar plot of secondary structure regions containing mutants in each mutant Cas9
    bar_width = 0.45
    plt.bar(np.arange(len(list_to_plot)), [dd[1] for dd in list_to_plot], width=bar_width, align='center', color=col)
    plt.xticks(range(len(list_to_plot)), [dd[0] for dd in list_to_plot], rotation=45, ha='right')
    plt.title(plot_title)
    plt.xlabel(xlab)
    plt.ylabel(ylab)
    if log_scale:
        plt.yscale('log')
        plt.ylim(min_count-0.1)  # Show values with just the minimum count
    plt.show()

def main():

    sample='q'
    sm_bin='10.0_10.5'
    catalogue = 'sm_9.5_s0.2_sfr_c-0.75_250'

    #load in fiducial mock
    filepath = './'
    filename = 'sm_9.5_s0.2_sfr_c-0.8_Chinchilla_250_wp_fiducial_'+sample+'_'+sm_bin+'_cov.npy'
    cov = np.matrix(np.load(filepath+filename))
    diag = np.diagonal(cov)
    filepath = cu.get_output_path() + 'analysis/central_quenching/observables/'
    filename = 'sm_9.5_s0.2_sfr_c-0.8_Chinchilla_250_wp_fiducial_'+sample+'_'+sm_bin+'.dat'
    data = ascii.read(filepath+filename)
    rbins = np.array(data['r'])
    mu = np.array(data['wp'])
    
    #load in comparison mock
    
    
    
    
    plt.figure()
    plt.errorbar(rbins, mu, yerr=np.sqrt(np.diagonal(cov)), color='black')
    plt.plot(rbins, wp,  color='red')
    plt.xscale('log')
    plt.yscale('log')
    plt.show()
    
    inv_cov = cov.I
    Y = np.matrix((wp-mu))
    
    X = Y*inv_cov*Y.T
    
    print(X)

 def plot_samples(self):
     """ Plot the samples requested for each arm """
     plt.clf()
     plt.scatter(range(0, self.K), self.sample_size)
     plt.yscale('log')
     plt.ioff()
     plt.show()

def plottamelo3():
    a = np.arange(0,30)
    nw = NWsc(matr1,True)
    db = DBsc(matr1,True)
    dbPR = productDBsc(matr1,True)
    INsc = INsc1(matr1,True)
    plt.yscale('log')
    plt.plot(range(0,nw.iter),nw.res, color='blue', lw=2, label="Newton Sclaled")
    #plt.plot(range(0,db.iter),db.res, color='red', lw=2, label="DB Sclaled")
    #plt.plot(range(0,dbPR.iter),dbPR.res, color='green', lw=2, label="DB Product Sclaled")
    #plt.plot(range(0,INsc.iter),INsc.res, color='orange', lw=2, label="IN Sclaled")
    nw = newton(matr1,True)
    db = DB(matr1,True)
    dbPR = productDB(matr1,True)
    IN = INiteration(matr1,True)
    CR = CRiteration(matr1,True)
    plt.yscale('log')
    print len(range(0,nw.iter))
    print len(nw.res)
    plt.plot(range(0,nw.iter+1),nw.res, color='blue', lw=2, label="Newton")
    #plt.plot(range(0,db.iter+1),db.res, color='red', lw=2, label="DB")
    #plt.plot(range(0,dbPR.iter+1),dbPR.res, color='green', lw=2, label="DB Product")
    #plt.plot(range(0,IN.iter+1),IN.res, color='orange', lw=2, label="IN")
    #plt.plot(range(0,CR.iter+1),CR.res, color='brown', lw=2, label="CR")

    plt.legend(loc='upper right')
    plt.show()

def plot_gradient_over_time(points, get_grad_over_time):
    fig = plt.figure(figsize=(6.5, 4))
    # Remove the plot frame lines. They are unnecessary chartjunk.    
    ax = plt.subplot(111)    
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    # ax.xaxis.set_major_locator(plt.MultipleLocator(1.0))
    # ax.xaxis.set_minor_locator(plt.MultipleLocator(0.1))
    # ax.yaxis.set_major_locator(plt.MultipleLocator(1.0))
    # ax.yaxis.set_minor_locator(plt.MultipleLocator(0.1))
    # ax.grid(which='major', axis='x', linewidth=0.75, linestyle='-', color='0.75')
    # ax.grid(which='minor', axis='x', linewidth=0.25, linestyle='-', color='0.75')
    # ax.grid(which='major', axis='y', linewidth=0.75, linestyle='-', color='0.75')
    # ax.grid(which='minor', axis='y', linewidth=0.25, linestyle='-', color='0.75')
    ax.grid(b=True, which='major', linewidth=0.75, linestyle=':', color='0.75')
    ax.yaxis.set_ticks_position('none')
    ax.xaxis.set_ticks_position('none')
    for wx, wRec, c in points:
        grad_over_time = get_grad_over_time(wx, wRec)
        x = np.arange(1, grad_over_time.shape[1]+1, 1)
        plt.plot(x, np.sum(grad_over_time, axis=0), c+'.-', label='({0}, {1})'.format(wx, wRec), linewidth=1, markersize=8)
    plt.xlim(1, grad_over_time.shape[1])
    plt.xticks(x)
    plt.gca().invert_xaxis()
    plt.yscale('symlog')
    plt.yticks([10**8, 10**6, 10**4, 10**2, 0, -10**2, -10**4, -10**6, -10**8])
    plt.xlabel('time k')
    plt.ylabel('gradient ')
    plt.title('Unstability of gradient in backward propagation.\n(backpropagate from left to right)')
    leg = plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False, numpoints=1)
    leg_font = FontProperties()
    leg_font.set_size('x-large')
    leg.set_title('$(w_x, w_{rec})$', prop=leg_font)

def create_time_plot(data, times, ylabel, title=None, filename=None,
                     xlims=None, ylims=None, log=False):
    plt.close('all')
    fig, ax = plt.subplots()
    if type(data) == dict:
        for lab in data.keys():
            plt.plot(times[lab], data[lab], 'x', label=lab)
        # plt.legend(loc='upper left')
        plt.legend()
    else:
        plt.plot(times, data, 'x')

    plt.grid()
    plt.xlabel("Received Time")
    plt.ylabel(ylabel)
    if title is not None:
        plt.title(title)
    fig.autofmt_xdate()
    if xlims is not None:
        plt.xlim(xlims)
    if ylims is not None:
        plt.ylim(ylims)
    if log:
        plt.yscale('log')
    if filename is not None:
        plt.tight_layout()
        plt.savefig(filename, fmt='pdf')
    else:
        plt.show()

def plot_citation_graph(citation_graph, filename, plot_title):
    # find the indegree_distribution
    indeg_dist = in_degree_distribution(citation_graph)    
    # sort freq by keys
    number_citations = sorted(indeg_dist.keys())
    indeg_freq = [indeg_dist[n] for n in number_citations]

    # normalize
    total = sum(indeg_freq)
    indeg_freq_norm = [freq / float(total) for freq in indeg_freq]
    
    # calculate log/log, except for the first one (0)
    #log_number_citations = [math.log10(x) for x in number_citations[1:]]
    #log_indeg_freq_norm = [math.log10(x) for x in indeg_freq_norm[1:]]
    
    plot(number_citations[1:], indeg_freq_norm[1:], 'o')
    
    xscale("log")
    yscale("log")
    
    xlabel("log10 #citations")
    ylabel("log10 Norm.Freq.")
    title(plot_title)
    grid(True)
    savefig(filename)
    show()

def plotCurves(losses,rateOfExceedance,return_periods,lossLevels):

    plt.figure(1, figsize=(8, 6), dpi=80, facecolor='w', edgecolor='k')
    plt.scatter(losses,rateOfExceedance,s=20)
    if len(return_periods) > 0:
        annual_rate_exc = 1.0/np.array(return_periods)
        for rate in annual_rate_exc:
            if rate > min(rateOfExceedance):
                plt.plot([min(losses),max(losses)],[rate,rate],color='red') 
                plt.annotate('%.6f' % rate,xy=(max(losses),rate),fontsize = 12)

    plt.yscale('log')
    plt.xscale('log')
    plt.ylim([min(rateOfExceedance),1])
    plt.xlim([min(losses),max(losses)])
    plt.xlabel('Losses', fontsize = 16)
    plt.ylabel('Annual rate of exceedance', fontsize = 16)

    setReturnPeriods = 1/rateOfExceedance
    plt.figure(2, figsize=(8, 6), dpi=80, facecolor='w', edgecolor='k')
    plt.scatter(setReturnPeriods,losses,s=20)
    if len(return_periods) > 0:
        for period in return_periods:
            if period < max(setReturnPeriods):
                plt.plot([period,period],[min(losses),max(losses)],color='red') 
                plt.annotate(str(period),xy=(period,max(losses)),fontsize = 12)

    plt.xscale('log')
    plt.xlim([min(setReturnPeriods),max(setReturnPeriods)])
    plt.ylim([min(losses),max(losses)])
    plt.xlabel('Return period (years)', fontsize = 16)
    plt.ylabel('Losses', fontsize = 16)

def plot():
    plt.plot(x_large, y_large, label=u'grafo de tamaño grande')
    plt.plot(x_medium, y_medium, label=u'grafo de tamaño medio')
    plt.yscale('log')
    plt.legend()
    plt.xlabel(u'iteración')
    plt.ylabel(u'diferencia en la norma del vector resultado entre sucesivas iteraciones')

    def __save(self, n, plot, sfile):
        p.figure(figsize=sfile)

        p.xlabel(plot.xlabel)
        p.ylabel(plot.ylabel)
        p.xscale(plot.xscale)
        p.yscale(plot.yscale)
        p.grid()
        for curvetype, args, kwargs in plot.curves:
            if curvetype == "plot":
                p.plot(*args, **kwargs)
            elif curvetype == "imshow":
                p.imshow(*args, **kwargs)
            elif curvetype == "hist":
                p.hist(*args, **kwargs)
            elif curvetype == "bar":
                p.bar(*args, **kwargs)

        p.axes().set_aspect(plot.aspect)
        if plot.legend:
            p.legend(shadow=0, loc=plot.loc)

        if not os.path.isdir(plot.dir):
            os.mkdir(plot.dir)
        if plot.pgf:
            p.savefig(plot.dir + plot.name + ".pgf")
            print(plot.name + ".pgf")
        if plot.pdf:
            p.savefig(plot.dir + plot.name + ".pdf", bbox_inches="tight")
            print(plot.name + ".pdf")

        p.close()

def make_plot():
    # Set up figure
    fig = plot.figure(figsize = (700 / my_dpi, 600 / my_dpi), dpi = my_dpi)

    ### Plot ###

    for i, mode in enumerate(default_modes):
        alpha = 0.4
        if mode == 1:
            alpha = 1.0
        if mode == 3 or mode == 5:
            alpha = 0.7
        plot.plot(frame_range, modes_over_time[i, :], linewidth = linewidth, alpha = alpha, label = "%d" % mode)

    # Axis
    plot.xlim(0, frame_range[-1])
    plot.ylim(10**(-3.5), 10**(0.0))
    plot.yscale("log")

    # Annotate
    this_title = readTitle()
    plot.xlabel("Number of Planet Orbits", fontsize = fontsize)
    plot.ylabel("Vortensity Mode Amplitudes", fontsize = fontsize)
    plot.title("%s" % (this_title), fontsize = fontsize + 1)

    plot.legend(loc = "upper right", bbox_to_anchor = (1.2, 1.0)) # outside of plot

    # Save and Close
    plot.savefig("fft_vortensity_modes.png", bbox_inches = 'tight', dpi = my_dpi)
    plot.show()
    plot.close(fig) # Close Figure (to avoid too many figures)

def saskia_plot_pulseheight( data_file_name, station_number=501, detector=0, number_bins=200, range_start=0., range_end=4500 ):
    
    # If the plot exist we skip the plotting
    if os.path.isfile('./img/pulseheigt_histogram_%d_detector_%d.pdf' % (station_number, detector)):
        # Say if the plot is present
        print "Plot already present for station %d" % station_number

    # If there is no plot we make it
    else:
   
        # Now transform the ROOT histogram to a python figure
        rootpy.plotting.root2matplotlib.hist(ph_histo)
    
        # Setting the limits on the axis
        plt.ylim((pow(10,-1),pow(10,7)))
        plt.xlim((range_start, range_end))
        plt.yscale('log')
    
        # Setting the plot labels and title
        plt.xlabel("Pulseheight [ADC]")
        plt.ylabel("Counts")
        plt.title("Pulseheight histogram (log scale) for station (%d)" %station_number)
    
        # Saving them Pica
        plt.savefig(
            './img/pulseheigt_histogram_%d_detector_%d.pdf' % (station_number, detector) ,    # Name of the file
            bbox_inches='tight')                        # Use less whitespace

def avgDegree(G):
  print "Nodes: ",
  print G.number_of_nodes()
  print "Edges: ",
  print G.number_of_edges()

  # avg degree
  degrees = defaultdict(int)
  total = 0
  for node in G.nodes():
    neighbors = G.neighbors(node)
    degrees[len(neighbors)] += 1
    total += len(neighbors)

  max_degree = max(degrees.keys())
  degrees_arr = (max_degree+1) * [0]
  for index, count in degrees.iteritems():
    degrees_arr[index] = count

  plt.plot(range(max_degree+1), degrees_arr, '.')
  plt.xscale('log', basex=2)
  plt.xlabel('degree')
  plt.yscale('log', basex=2)
  plt.ylabel('# of people')
  plt.savefig('degree_distribution.png')
  plt.close()

 def check_hod(self, z, prop):
     data = np.genfromtxt(LOCATION + "/data/" + prop + "z" + str(z))
     if prop == "ncen":
         if PLOT:
             plt.clf()
             plt.plot(self.hod.hmf.M,
                      self.hod.n_cen,
                      label="mine")
             plt.plot(data[:, 0] * self.hod.cosmo.h, data[:, 1], label="charles")
             plt.legend()
             plt.xscale('log')
             plt.yscale('log')
             plt.savefig(join(pref, "ncen" + prop + "z" + str(z) + ".pdf"))
         assert max_diff_rel(self.hod.n_cen, data[:, 1], 0.01)
     elif prop == "nsat":
         if PLOT:
             plt.clf()
             plt.plot(self.hod.hmf.M,
                      self.hod.n_sat,
                      label="mine")
             plt.plot(data[:, 0] * self.hod.cosmo.h, data[:, 1], label="charles")
             plt.legend()
             plt.xscale('log')
             plt.yscale('log')
             plt.savefig(join(pref, "nsat" + prop + "z" + str(z) + ".pdf"))
         assert max_diff_rel(self.hod.n_sat, data[:, 1], 0.01)

def make_plot(filename, title, arguments, methods, scale):
    if not support_plots:
        return

    is_linear = (scale == 'linear')

    plot_size = LINEAR_PLOT_SIZE if is_linear else OTHER_PLOT_SIZE
    plt.figure(figsize=plot_size)

    for name, func, measures in methods:
        plt.plot(arguments, measures, 'o-', label=name, markersize=3)

    if is_linear:
        axis = plt.axis()
        plt.axis((0, axis[1], 0, axis[3]))

    plt.xscale(scale)
    plt.yscale(scale)

    plt.xticks(fontsize=NORMAL_FONT_SIZE)
    plt.yticks(fontsize=NORMAL_FONT_SIZE)

    plt.grid(True)

    plt.title(title, fontsize=LABEL_FONT_SIZE)
    plt.xlabel('Argument', fontsize=NORMAL_FONT_SIZE)
    plt.ylabel('Time (seconds)', fontsize=NORMAL_FONT_SIZE)
    plt.legend(loc='upper left', fontsize=NORMAL_FONT_SIZE)

    plt.tight_layout(0.2)

    path = os.path.join(PLOTS_DIR, filename)
    plt.savefig(path)
    print '[*] Saved plot "%s"' % path

    def _init_ipython_plot(self,
                           plot_data,
                           legend_lst=[],
                           title_str=' ',
                           axis=None,
                           plot_relative=False):
        """ Private member function that loads ipython plotting libraries 

        Parameters
        ----------
        plot_data : lst of lst
            A list of lists that stores the data series to be plotted.
        legend : lst
            A list of strings for the legend.
        axis : bool
           Indicates whether there is a user specified axis.
        plot_relative : bool
           Indicates whether the relative residuals should be plotted.
        """
        import matplotlib
        import numpy as np
        import matplotlib.pyplot as plt       
        for data_set in plot_data:
            if plot_relative == True:
                max_term = max(data_set)
                for i,term in enumerate(data_set):
                    data_set[i] = data_set[i] / max_term
            plt.plot(data_set)
        plt.yscale("log")
        plt.legend(legend_lst)
        plt.title(title_str)
        if axis is not None:
            plt.xlim(axis[0],axis[1])
        plt.show() 

def Validation():
  numSamples = 1000000
  
  theta = np.random.rand(numSamples)*np.pi
  ECo60 = np.array([1.117,1.332])
  Ef0,Ee0 = Compton(ECo60[0],theta)
  Ef1,Ee1 = Compton(ECo60[1],theta)
  dSdE0 = diffXSElectrons(ECo60[0],theta)
  dSdE1 = diffXSElectrons(ECo60[1],theta)

  # Sampling Values
  values = list()
  piMax = np.max([dSdE0,dSdE1])
  while (len(values) < numSamples):
    values.append(SampleRejection(piMax,ComptonScattering))
  # Binning the data
  bins = np.logspace(-3,0.2,100)
  counts = np.histogram(values,bins)
  counts = counts[0]/float(len(values))
  binCenters = 0.5*(bins[1:]+bins[:-1])
  
  # Plotting
  plt.figure()
  plt.plot(binCenters,counts,ls='steps')
  #plt.bar(binCenters,counts,align='center')
  plt.grid(True)
  plt.xlim((1E-3,1.4))
  plt.xlabel('Electron Energy (MeV)')
  plt.ylabel('Frequency per Photon')
  plt.yscale('log')
  plt.xscale('log')
  plt.savefig('ValComptonScatteringXS.png')