def chooseDegree(npts, mindegree=0, maxdegree=20, filename=None):
    """Gets noisy data, uses cross validation to estimate error, and fits new data with best model."""
    x, y = bv.noisyData(npts)
    degrees = numpy.arange(mindegree,maxdegree+1)
    errs = numpy.zeros_like(degrees,dtype=numpy.float)
    for i,d in enumerate(degrees):
        errs[i] = estimateError(x, y, d)

    plt.subplot(1,2,1)
    plt.plot(degrees,errs,'bo-')
    plt.xlabel("Degree")
    plt.ylabel("CV Error")

    besti = numpy.argmin(errs)
    bestdegree = degrees[besti]

    plt.subplot(1,2,2)
    x2, y2 = bv.noisyData(npts)
    plt.plot(x2,y2,'ro')
    xs = numpy.linspace(min(x),max(x),150)
    fitf = numpy.poly1d(numpy.polyfit(x2,y2,bestdegree))
    plt.plot(xs,fitf(xs),'g-',lw=2)
    plt.xlim((bv.MIN,bv.MAX))
    plt.ylim((-2.,2.))
    plt.suptitle('Selected Degree '+str(bestdegree))
    bv.outputPlot(filename)

def graphical_analysis(strace, comment=None, cutoff=50, threshold=-45, bins=50):
# -----------------------------------------------------------------------------
    """
    Graphical report of the trace.
    strace - voltage trace of class SimpleTrace
    cutoff - first part of the trace cut away from the analysis (in ms)
    threshold - cutting voltage
    """
    import matplotlib.pylab as pyl
    voltage = strace._data
    time = np.arange(len(voltage))*strace._dt
    pyl.figure()
    rate, mean, var, skew, kurt = full_analysis(strace, cutoff, threshold)
    if comment:
        pyl.suptitle(comment)
    sp1 = pyl.subplot(2,1,1)
    sp1.plot(time,voltage)
    sp1.set_title("Spike rate = {0}".format(rate))
    sp1.set_xlabel("time [ms]")
    sp1.set_ylabel("V [mV]")
    sp2 = pyl.subplot(2,1,2)
    cut_trace = voltage[int(cutoff/strace._dt):]
    data = cut_trace[cut_trace<threshold]
    sp2.hist(data, bins=bins, histtype="stepfilled", normed=1)
    xlim = sp2.get_xlim()
    pyl.text(xlim[0]+0.7*(xlim[1]-xlim[0]), 0.6,
             "mean = {0}\nvar={1}\nskew={2}\nkurt={3}".format(mean, var, skew, kurt))
    sp2.plot([mean, mean],[0,1],"r")
    sp2.plot([mean-np.sqrt(var)/2., mean+np.sqrt(var)/2.],[0.1,0.1], "r")
    sp2.set_xlabel("V [mV]")
    sp2.set_ylabel("normalized distribution")

 def plot(self,typ='s3',titre='titre',log=False,stem=True,subp=True):
     """
     """
     fa = np.linspace(self.Br.fmin,self.Br.fmax,self.Br.Nf)
     st = titre+'  shape : '+typ
     plt.suptitle(st,fontsize=14)
     if subp:
         plt.subplot(221)
         titre = '$\sum_f \sum_m |Br_{l}^{(m)}(f)|$'
         self.Br.plot(typ=typ,title=titre, yl=True,color='r',stem=stem,log=log)
     else:
         self.Br.plot(typ=typ,color='r',stem=stem,log=log)
     if subp:
         plt.subplot(222)
         titre = '$\sum_f \sum_m |Bi_{l}^{(m)}(f)|$'
         self.Bi.plot(typ=typ,title=titrei,color='m',stem=stem,log=log)
     else:
         self.Bi.plot(typ=typ,color='m',stem=stem,log=log)
     if subp:
         plt.subplot(223)
         titre = '$\sum_f \sum_m |Cr_{l}^{(m)}(f)|$'
         self.Cr.plot(typ=typ,title=titre, xl=True, yl=True,color='b',stem=stem,log=log)
     else:
         self.Cr.plot(typ=typ,color='b',stem=stem,log=log)
     if subp:
         plt.subplot(224)
         titre = '$\sum_f \sum_m |Ci_{l}^{(m)}(f)|$'
         self.Ci.plot(typ=typ, title = titre, xl=True,color='c',stem=stem,log=log)
     else:
         self.Ci.plot(typ=typ,xl=True,yl=True,color='c',stem=stem,log=log)
     if not subp:
         plt.legend(('$\sum_f \sum_m |Br_{l}^{(m)}(f)|$',
                     '$\sum_f \sum_m |Bi_{l}^{(m)}(f)|$',
                     '$\sum_f \sum_m |Cr_{l}^{(m)}(f)|$',
                     '$\sum_f \sum_m |Ci_{l}^{(m)}(f)|$'))

    def plotFittingResults(self):
        """
        Plot results of Rmax optimization procedure and best fit of the experimental data
        """
        _listFitQ = [tmp.getValue() for tmp in self.getDataOutput().getScatteringFitQ()]
        _listFitValues = [tmp.getValue() for tmp in self.getDataOutput().getScatteringFitValues()]
        _listExpQ = [tmp.getValue() for tmp in self.getDataInput().getExperimentalDataQ()]
        _listExpValues = [tmp.getValue() for tmp in self.getDataInput().getExperimentalDataValues()]

        #_listExpStdDev = None
        #if self.getDataInput().getExperimentalDataStdDev():
        #    _listExpStdDev = [tmp.getValue() for tmp in self.getDataInput().getExperimentalDataStdDev()]
        #if _listExpStdDev:
        #    pylab.errorbar(_listExpQ, _listExpValues, yerr=_listExpStdDev, linestyle='None', marker='o', markersize=1,  label="Experimental Data")
        #    pylab.gca().set_yscale("log", nonposy='clip')
        #else:         
        #    pylab.semilogy(_listExpQ, _listExpValues, linestyle='None', marker='o', markersize=5,  label="Experimental Data")

        pylab.semilogy(_listExpQ, _listExpValues, linestyle='None', marker='o', markersize=5, label="Experimental Data")
        pylab.semilogy(_listFitQ, _listFitValues, label="Fitting curve")
        pylab.xlabel('q')
        pylab.ylabel('I(q)')
        pylab.suptitle("RMax : %3.2f. Fit quality : %1.3f" % (self.getDataInput().getRMax().getValue(), self.getDataOutput().getFitQuality().getValue()))
        pylab.legend()
        pylab.savefig(os.path.join(self.getWorkingDirectory(), "gnomFittingResults.png"))
        pylab.clf()

    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 scatter(title, file_name, x_array, y_array, size_array, x_label, \
             y_label, x_range, y_range, print_pdf):
     '''
     Plots the given x value array and y value array with the specified 
     title and saves with the specified file name. The size of points on
     the map are proportional to the values given in size_array. If 
     print_pdf value is 1, the image is also written to pdf file. 
     Otherwise it is only written to png file.
     '''
     rc('text', usetex=True)
     rc('font', family='serif')
     plt.clf() # clear the ploting window, a must.                               
     plt.scatter(x_array, y_array, s =  size_array, c = 'b', marker = 'o', alpha = 0.4)
     if x_label != None:   
         plt.xlabel(x_label)
     if y_label != None:
         plt.ylabel(y_label)                
     plt.axis ([0, x_range, 0, y_range])
     plt.grid(True)
     plt.suptitle(title)
 
     Plotter.print_to_png(plt, file_name)
     
     if print_pdf:
         Plotter.print_to_pdf(plt, file_name)

def plot_fcs(normed_df, unnormed_df, pair, basename):
    """
    Plot fold changes for normed and unnormed dataframes.

    Parameters:
    -----------
    normed_df: Normalized dataframe of values
    unnormed_df: Unnormalized dataframe of values
    pair: Tuple containing the two columns to use
    to compute fold change. Fold change is first
    sample divided by second.
    """
    if (pair[0] not in normed_df.columns) or \
       (pair[1] not in normed_df.columns):
        raise Exception, "One of the pairs is not in normed df."
    if (pair[0] not in unnormed_df.columns) or \
       (pair[1] not in unnormed_df.columns):
        raise Exception, "One of the pairs is not in unnormed.df"
    normed_fc = \
        np.log2(normed_df[pair[0]]) - np.log2(normed_df[pair[1]])
    unnormed_fc = \
        np.log2(unnormed_df[pair[0]]) - np.log2(unnormed_df[pair[1]])
    fc_df = pandas.DataFrame({"normed_fc": normed_fc,
                              "unnormed_fc": unnormed_fc})
    # Remove nans/infs etc.
    pandas.set_option('use_inf_as_null', True)
    fc_df = fc_df.dropna(how="any", subset=["normed_fc",
                                            "unnormed_fc"])
    plt.figure()
    fc_df.hist(color="k", bins=40)
    plt.suptitle("%s vs. %s" %(pair[0], pair[1]))
    plt.xlabel("Fold change (log2)")
    save_fig(basename)
    pandas.set_option('use_inf_as_null', False)    

    def update(frame):
        global _counter
        centroid = np.random.uniform(-0.5, 0.5, size=2)
        width = np.random.uniform(0, 0.01)
        length = np.random.uniform(0, 0.03) + width
        angle = np.random.uniform(0, 360)
        intens = np.random.exponential(2) * 50
        model = mock.generate_2d_shower_model(centroid=centroid,
                                              width=width,
                                              length=length,
                                              psi=angle * u.deg)
        image, sig, bg = mock.make_mock_shower_image(geom, model.pdf,
                                                     intensity=intens,
                                                     nsb_level_pe=5000)

        # alternate between cleaned and raw images
        if _counter > 20:
            plt.suptitle("Image Cleaning ON")
            cleanmask = reco.cleaning.tailcuts_clean(geom, image, pedvars=80)
            for ii in range(3):
                reco.cleaning.dilate(geom, cleanmask)
            image[cleanmask == 0] = 0  # zero noise pixels
        if _counter >= 40:
            plt.suptitle("Image Cleaning OFF")
            _counter = 0

        disp.image = image
        disp.set_limits_percent(100)
        _counter += 1

def q6(abreviation):
	'''#Q6: PLOTS polling data for a given state'''
	#get STATE; connect to db
	state = abreviation.upper()
	connection = sqlite3.connect(path+"/poll.db")
	cursor = connection.cursor()
	
	#query db, get names, rankings for Rep, Dem, Ind candidates
	sql_cmd = "SELECT candidate_names.democrat, candidate_names.republican, candidate_names.independent, rankings.day, rankings.dem, rankings.rep, rankings.indep from rankings left join statetable on rankings.state = statetable.fullname\
	 left join candidate_names on statetable.abrev = candidate_names.state where  statetable.abrev = '%s' order by rankings.day ASC" % state
	cursor.execute(sql_cmd)
	dbinfo = cursor.fetchall()

	#'new' is the name of the record array corresponding to the output from the sql query
	new = N.array(dbinfo, dtype= [('demname', '|S25'),('repname', '|S25'),('indname', '|S25'),('day', N.int16), ("dem", N.int16), ("rep", N.int16), ("ind", N.int16)])

	demname= new['demname'][0] #name of democrat
	repname=  new['repname'][0] # name of rep
	indname= new['indname'][0] #" of ind

	# 2 (+1) lines for dem, rep, candiates ind if he has more than 1 point at the first datapoint, indicating there is an independent candidate
	plt.plot(new['day'],new['dem'], color='blue', label='%s' % demname)
	plt.plot(new['day'],new['rep'], color='red', label='%s' % repname)
	if new['ind'][0] > 1:
		plt.plot(new['day'],new['ind'], color='green', label='%s' % indname)
		
	#plot info
	plt.suptitle('Election Polls for the State of %s'%state)
	plt.xlabel('Day of the year')
	plt.ylabel('Points')
	plt.legend()
	plt.show()

def plot_ss_scatter(steadies):
    """ Plot scatter plots of steady states
    """

    def do_scatter(i, j, ax):
        """ Draw single scatter plot
        """
        xs, ys = utils.extract(i, j, steadies)
        ax.scatter(xs, ys)

        ax.set_xlabel(r"$S_%d$" % i)
        ax.set_ylabel(r"$S_%d$" % j)

        cc = utils.get_correlation(xs, ys)
        ax.set_title(r"Corr: $%.2f$" % cc)

    dim = steadies.shape[1]
    fig, axarr = plt.subplots(1, int((dim ** 2 - dim) / 2), figsize=(20, 5))

    axc = 0
    for i in range(dim):
        for j in range(dim):
            if i == j:
                break
            do_scatter(i, j, axarr[axc])
            axc += 1

    plt.suptitle("Correlation overview")

    plt.tight_layout()
    save_figure("images/correlation_scatter.pdf", bbox_inches="tight")
    plt.close()

def plotErrorAndOrder(schemesName, spaceErrorList,temporalErrorList,
                      spaceOrderList, temporalOrderList, Ntds):
    legendList = []
    lstyle = ['b', 'r', 'g', 'm']
    fig , axarr = plt.subplots(2, 2, squeeze=False)
    for k, scheme_name in enumerate(schemesName):
        axarr[0][0].plot(np.log10(np.asarray(spaceErrorList[k])),lstyle[k])
        axarr[0][1].plot(np.log10(np.asarray(temporalErrorList[k])),lstyle[k])
        axarr[1][0].plot(spaceOrderList[k],lstyle[k])
        axarr[1][1].plot(temporalOrderList[k],lstyle[k])
        legendList.append(scheme_name)
    plt.suptitle('test_MES_convergence(): Results from convergence test using Method of Exact Solution')

    axarr[1][0].axhline(1.0, xmin=0, xmax=Ntds-2, linestyle=':', color='k')
    axarr[1][0].axhline(2.0, xmin=0, xmax=Ntds-2, linestyle=':', color='k')
    axarr[1][1].axhline(1.0, xmin=0, xmax=Ntds-2, linestyle=':', color='k')
    axarr[1][1].axhline(2.0, xmin=0, xmax=Ntds-2, linestyle=':', color='k')
    axarr[1][0].set_ylim(0, 5)
    axarr[1][1].set_ylim(0, 5)
    axarr[0][0].set_ylabel('rms Error')
    axarr[0][0].set_title('space Error')
    axarr[1][0].set_ylabel('rms Error')
    axarr[0][1].set_title('temporal Error')
    axarr[1][0].set_ylabel('order')
    axarr[1][0].set_title('space order')
    axarr[1][1].set_ylabel('order')
    axarr[1][1].set_title('temporal order')
    axarr[0][1].legend(legendList, frameon=False)

def report(path='history.cpkl', tn=0, sns=True):
    if sns:
        import seaborn as sns
        sns.set_style('whitegrid')
        sns.set_style('whitegrid', {'fontsize': 50})
        sns.set_context('poster')
    with open(path) as f:
        logged_data = pickle.load(f)

    history = util.NestedDict()
    for name, val in logged_data.iteritems():
        history.set_nested(name, val)

    num_subplots = len(history)
    cols = 2  # 2 panels for Objective and Accuracy
    rows = 1

    fig = plt.figure(figsize=(12, 8))
    fig.subplots_adjust(wspace=0.3, hspace=0.2)  # room for labels [Objective, Accuracy]
    colors = [sns.xkcd_rgb['blue'], sns.xkcd_rgb['red']]

    # Here we assume that history is only two levels deep
    for k, (subplot_name, trend_lines) in enumerate(history.iteritems()):
        plt.subplot(rows, cols, k + 1)
        plt.ylabel(subplot_name.capitalize())
        plt.xlabel('Epoch')
        for i, (name, (timestamps, values)) in enumerate(trend_lines.iteritems()):
            plt.plot(timestamps, values, label=name, color=colors[i])
        plt.suptitle('Task number %d' % tn)
        plt.legend(loc='best')

    plt.show()

def plotAstrometry(dist, mag, snr, brightSnr=100,
                   outputPrefix=""):
    """Plot angular distance between matched sources from different exposures.

    Creates a file containing the plot with a filename beginning with `outputPrefix`.

    Parameters
    ----------
    dist : list or numpy.array
        Separation from reference [mas]
    mag : list or numpy.array
        Mean magnitude of PSF flux
    snr : list or numpy.array
        Median SNR of PSF flux
    brightSnr : float, optional
        Minimum SNR for a star to be considered "bright".
    outputPrefix : str, optional
        Prefix to use for filename of plot file.  Will also be used in plot titles.
        E.g., outputPrefix='Cfht_output_r_' will result in a file named
           'Cfht_output_r_check_astrometry.png'
    """
    bright, = np.where(np.asarray(snr) > brightSnr)

    numMatched = len(dist)
    dist_median = np.median(dist)
    bright_dist_median = np.median(np.asarray(dist)[bright])

    fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(18, 12))

    ax[0].hist(dist, bins=100, color=color['all'],
               histtype='stepfilled', orientation='horizontal')
    ax[0].hist(np.asarray(dist)[bright], bins=100, color=color['bright'],
               histtype='stepfilled', orientation='horizontal',
               label='SNR > %.0f' % brightSnr)

    ax[0].set_ylim([0., 500.])
    ax[0].set_ylabel("Distance [mas]")
    ax[0].set_title("Median : %.1f, %.1f mas" %
                    (bright_dist_median, dist_median),
                    x=0.55, y=0.88)
    plotOutlinedLinesHorizontal(ax[0], dist_median, bright_dist_median)

    ax[1].scatter(snr, dist, s=10, color=color['all'], label='All')
    ax[1].scatter(np.asarray(snr)[bright], np.asarray(dist)[bright], s=10,
                  color=color['bright'],
                  label='SNR > %.0f' % brightSnr)
    ax[1].set_xlabel("SNR")
    ax[1].set_xscale("log")
    ax[1].set_ylim([0., 500.])
    ax[1].set_title("# of matches : %d, %d" % (len(bright), numMatched))
    ax[1].legend(loc='upper left')
    ax[1].axvline(brightSnr, color='red', linewidth=4, linestyle='dashed')
    plotOutlinedLinesHorizontal(ax[1], dist_median, bright_dist_median)

    plt.suptitle("Astrometry Check : %s" % outputPrefix.rstrip('_'), fontsize=30)
    plotPath = outputPrefix+"check_astrometry.png"
    plt.savefig(plotPath, format="png")
    plt.close(fig)

def plot_generated_toy_batch(X_real, generator_model, discriminator_model, noise_dim, gen_iter, noise_scale=0.5):

    # Generate images
    X_gen = sample_noise(noise_scale, 10000, noise_dim)
    X_gen = generator_model.predict(X_gen)

    # Get some toy data to plot KDE of real data
    data = load_toy(pts_per_mixture=200)
    x = data[:, 0]
    y = data[:, 1]
    xmin, xmax = -1.5, 1.5
    ymin, ymax = -1.5, 1.5

    # Peform the kernel density estimate
    xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
    positions = np.vstack([xx.ravel(), yy.ravel()])
    values = np.vstack([x, y])
    kernel = stats.gaussian_kde(values)
    f = np.reshape(kernel(positions).T, xx.shape)

    # Plot the contour
    fig = plt.figure(figsize=(10,10))
    plt.suptitle("Generator iteration %s" % gen_iter, fontweight="bold", fontsize=22)
    ax = fig.gca()
    ax.contourf(xx, yy, f, cmap='Blues', vmin=np.percentile(f,80), vmax=np.max(f), levels=np.linspace(0.25, 0.85, 30))

    # Also plot the contour of the discriminator
    delta = 0.025
    xmin, xmax = -1.5, 1.5
    ymin, ymax = -1.5, 1.5
    # Create mesh
    XX, YY = np.meshgrid(np.arange(xmin, xmax, delta), np.arange(ymin, ymax, delta))
    arr_pos = np.vstack((np.ravel(XX), np.ravel(YY))).T
    # Get Z = predictions
    ZZ = discriminator_model.predict(arr_pos)
    ZZ = ZZ.reshape(XX.shape)
    # Plot contour
    ax.contour(XX, YY, ZZ, cmap="Blues", levels=np.linspace(0.25, 0.85, 10))
    dy, dx = np.gradient(ZZ)
    # Add streamlines
    # plt.streamplot(XX, YY, dx, dy, linewidth=0.5, cmap="magma", density=1, arrowsize=1)
    # Scatter generated data
    plt.scatter(X_gen[:1000, 0], X_gen[:1000, 1], s=20, color="coral", marker="o")

    l_gen = plt.Line2D((0,1),(0,0), color='coral', marker='o', linestyle='', markersize=20)
    l_D = plt.Line2D((0,1),(0,0), color='steelblue', linewidth=3)
    l_real = plt.Rectangle((0, 0), 1, 1, fc="steelblue")

    # Create legend from custom artist/label lists
    # bbox_to_anchor = (0.4, 1)
    ax.legend([l_real, l_D, l_gen], ['Real data KDE', 'Discriminator contour',
                                     'Generated data'], fontsize=18, loc="upper left")
    ax.set_xlim(xmin, xmax)
    ax.set_ylim(ymin, ymax + 0.8)
    plt.savefig("../../figures/toy_dataset_iter%s.jpg" % gen_iter)
    plt.clf()
    plt.close()

def plot_contours(obj, top_bottom=True):
  '''A function that plots the BRF as an azimuthal projection
  with contours over the TOC and soil.
  Input: rt_layers object, top_bottom - True if only TOC plot, False
  if both TOC and soil.
  Output: contour plot of brf.
  '''
  sun = ((np.pi - obj.sun0[0]) * np.cos(obj.sun0[1] + np.pi), \
      (np.pi - obj.sun0[0]) * np.sin(obj.sun0[1] + np.pi))
  theta = obj.views[:,0]
  x = np.cos(obj.views[:,1]) * theta
  y = np.sin(obj.views[:,1]) * theta
  z = obj.I_top_bottom # * -obj.mu_s
  if top_bottom == True:
    if np.max > 1.:
      maxz = np.max(z)
    else:
      maxz = 1.
  else:
    maxz = np.max(z[:obj.n/2])
  minz = 0. #np.min(z)
  space = np.linspace(minz, maxz, 11)
  x = x[:obj.n/2]
  y = y[:obj.n/2]
  zt = z[:obj.n/2]
  zb = z[obj.n/2:]
  fig = plt.figure()
  if top_bottom == True:
    plt.subplot(121)
  plt.plot(sun[0], sun[1], 'ro')
  triang = tri.Triangulation(x, y)
  plt.gca().set_aspect('equal')
  plt.tricontourf(triang, zt, space, vmax=maxz, vmin=minz)
  plt.title('TOC BRF')
  plt.ylabel('Y')
  plt.xlabel('X')
  if top_bottom == True:
    plt.subplot(122)
    plt.plot(sun[0], sun[1], 'ro')
    plt.gca().set_aspect('equal')
    plt.tricontourf(triang, zb, space, vmax=maxz, vmin=minz)
    plt.title('Soil Absorption')
    plt.ylabel('Y')
    plt.xlabel('X')
  s = obj.__repr__()
  if top_bottom == True:
    cbaxes = fig.add_axes([0.11,0.1,0.85,0.05])
    plt.suptitle(s,x=0.5,y=0.93)
    plt.colorbar(orientation='horizontal', ticks=space,\
      cax = cbaxes, format='%.3f')
  else:
    plt.suptitle(s,x=0.5,y=0.13)
    plt.colorbar(orientation='horizontal', ticks=space,\
        format='%.3f')
    #plt.tight_layout()
  plt.show()

def display(data):
    i = 0
    pylab.figure()
    for d in data:
        pylab.scatter([d[0]], [d[1]], c="r", marker='o')

    pylab.suptitle("Generated teddy toy\ndim=(" + str(dim_x) + ", " + str(dim_y) + ") dist=" + str(clust_dist) + " size=" + str(clust_size) + " stddev=" + str(clust_stddev))

    pylab.savefig(filename + ".pdf")
    pylab.savefig(filename + ".eps")
    pylab.savefig(filename + ".svg")

def get_offset_center(f, plot=False, interactive=False):
    '''
    Given a fits image, returns the offset in Ra, DEC, that needs to be applied for the telescope tp go
    from the current pointing position, to the coodinates of the object specified in the fits file.
    '''
    
    if(not os.path.isfile(f)):
        print "File %s does not exist! Returning Zero offsets..."%f
        return -1, 0,0
    else:
        image = pf.open(f)
        wcs = pywcs.WCS(image[0].header)
        rra, rdec = cc.hour2deg(image[0].header['OBJRA'],image[0].header['OBJDEC'] )
        x, y = np.round(wcs.wcs_sky2pix(rra, rdec, 0), 0)
        pra, pdec = wcs.wcs_pix2sky(np.array([[1293., 1280.]] , np.float_), 0)[0]
        dra, ddec = cc.get_offset(pra, pdec, rra, rdec)
            
        xl, yu = np.round(wcs.wcs_sky2pix(rra+90./3600, rdec-90./3600, 0), 0)
        xu, yl = np.round(wcs.wcs_sky2pix(rra-90./3600, rdec+90./3600, 0), 0)

        imageloc = image[0].data.T[xl:xu,yl:yu]

        if imageloc.shape[0]==0 or imageloc.shape[1]==0:
            logger.warn( "Astrometry has FAILED on this! The object is outside the frame! Resending to the numb astrometric solution")
            logger.error("Astrometry has FAILED on this! The object is outside the frame! Resending to the numb astrometric solution")
            print "Pixels are", xl, xu, yl, yu
            try:
                code, dra, ddec = get_offset_center_failed_astro(f, plot=plot, interactive=interactive)
                return 2, dra, ddec
            except:
                return -1,0,0
        if(plot):
            plt.figure(figsize=(8,8))
            
            zmin, zmax = zscale.zscale(imageloc)

            #print zmin, zmax, imageloc, (xl,xu,yl,yu)
    
            obj = fitsutils.get_par(f, "OBJECT")
            plt.suptitle(obj, fontsize=20)
            plt.imshow(imageloc.T,  extent=(xl[0],xu[0],yl[0],yu[0]), aspect="equal", interpolation="none", origin="lower", vmin=zmin, vmax=zmax)
            plt.plot(1293., 1280., "ws", ms=7, label="Current pointing")
            plt.plot(x, y, "b*", ms=10, label="Target pointing")
            plt.gca().invert_xaxis()
            plt.legend()
            if (interactive):
                plt.show()
            else:
                plt.savefig(os.path.join(os.path.dirname(f).replace("raw", "phot"), os.path.basename(f).replace(".fits", "_a.png")))
            plt.clf()


        return 0, dra, ddec

def Ploting_Facts(P):
    F_2, F_3 = Run_Facts(P)
    plt.figure()
    plt.plot(F_3, F_2, 'r.')
    plt.plot(25,75, 'ko', label="Q")
    plt.xlabel('Fact #3')
    plt.ylabel('Fact #2')
    plt.suptitle('Classical Behavior of Incomunicated Students', fontsize = 14)
    plt.legend(loc=2)
    plt.grid(True)
    plt.savefig('Estudiantes_F3vF2_P=' + str(P) + '.png')
    plt.show()

def plotPA1(pa1, outputPrefix=""):
    """Plot the results of calculating the LSST SRC requirement PA1.

    Creates a file containing the plot with a filename beginning with `outputPrefix`.

    Parameters
    ----------
    pa1 : pipeBase.Struct
        Must contain:
        rms, iqr, magMean, magDiffs
        rmsUnits, iqrUnits, magDiffsUnits
    outputPrefix : str, optional
        Prefix to use for filename of plot file.  Will also be used in plot titles.
        E.g., outputPrefix='Cfht_output_r_' will result in a file named
           'Cfht_output_r_AM1_D_5_arcmin_17.0-21.5.png'
        for an AMx.name=='AM1' and AMx.magRange==[17, 21.5]
    """
    diffRange = (-100, +100)

    fig = plt.figure(figsize=(18, 12))
    ax1 = fig.add_subplot(1, 2, 1)
    ax1.scatter(pa1.magMean, pa1.magDiffs, s=10, color=color['bright'], linewidth=0)
    ax1.axhline(+pa1.rms, color=color['rms'], linewidth=3)
    ax1.axhline(-pa1.rms, color=color['rms'], linewidth=3)
    ax1.axhline(+pa1.iqr, color=color['iqr'], linewidth=3)
    ax1.axhline(-pa1.iqr, color=color['iqr'], linewidth=3)

    ax2 = fig.add_subplot(1, 2, 2, sharey=ax1)
    ax2.hist(pa1.magDiffs, bins=25, range=diffRange,
             orientation='horizontal', histtype='stepfilled',
             normed=True, color=color['bright'])
    ax2.set_xlabel("relative # / bin")

    yv = np.linspace(diffRange[0], diffRange[1], 100)
    ax2.plot(scipy.stats.norm.pdf(yv, scale=pa1.rms), yv,
             marker='', linestyle='-', linewidth=3, color=color['rms'],
             label="PA1(RMS) = %4.2f %s" % (pa1.rms, pa1.rmsUnits))
    ax2.plot(scipy.stats.norm.pdf(yv, scale=pa1.iqr), yv,
             marker='', linestyle='-', linewidth=3, color=color['iqr'],
             label="PA1(IQR) = %4.2f %s" % (pa1.iqr, pa1.iqrUnits))
    ax2.set_ylim(*diffRange)
    ax2.legend()
#    ax1.set_ylabel(u"12-pixel aperture magnitude diff (mmag)")
#    ax1.set_xlabel(u"12-pixel aperture magnitude")
    ax1.set_xlabel("psf magnitude")
    ax1.set_ylabel("psf magnitude diff (%s)" % pa1.magDiffsUnits)
    for label in ax2.get_yticklabels():
        label.set_visible(False)

    plt.suptitle("PA1: %s" % outputPrefix.rstrip('_'))
    plotPath = "%s%s" % (outputPrefix, "PA1.png")
    plt.savefig(plotPath, format="png")
    plt.close(fig)

def guess(velocity, thetaDeg, timeStep=0.01):
    """
    Initial calcultaion / integration of the function using guessed
    initial conditions. The function will also be plotted to help
    visualise the scenario.
    velocity in meters per second, theta in degrees, timeStep in
    seconds. The smaller the timeStep the more accurate the function.
    """
    # Calculate / fill all the variables
    theta = radians(thetaDeg)       # Change the input angle into radians
    trajGuessX = []                 # Create an empty list
    trajGuessX.append(trajXY[0])    # Initialise the list
    trajGuessY = []                 # Create an empty list
    trajGuessY.append(trajXY[1])    # Initialise the list
    accelX = ((rho*drag*pow((velocity*cos(theta)),2)*pi*pow(radius,2))
              /(2*mass))            # Calculate the x-acceleration
    accelY = grav                   # Assign gravity to y-acceleration
    velocityX = velocity*cos(theta) # Calculate the x-velocity component
    velocityY = velocity*sin(theta) # Calculate the y-velocity component

    # Execute the mathematics and build a list of the coordinates.
    # While the bread is in the air perform calculations.
    # As it steps through it updates the velocities and accelerations
    # so that the bread acceleration and velocity slows.
    while trajGuessY[-1] > 0:
        # New velocity equals old velocity minus the updated acceleration
        velocityX = velocityX - accelX*timeStep
        # Change the acceleraion to use the last calculated velocity
        accelX = ((rho*drag*pow(velocityX,2)*pi*pow(radius,2))
                   /(2*mass))
        # New velocity equals old velocity minus the updated acceleration
        velocityY = velocityY - accelY*timeStep
        # Positions equal last position (in the list) + distance moved
        x = trajGuessX[-1] + velocityX*timeStep
        y = trajGuessY[-1] + velocityY*timeStep
        trajGuessX.append(x)    # Append the x-coord to the list
        trajGuessY.append(y)    # Append the y-coord to the list

    # Plot the graphs of the trajectory and the physical environment
    env = plt.plot(physEnvX, physEnvY, 'b', label='Environment')
    traj = plt.plot(trajGuessX, trajGuessY, 'r--', label='Trajectory')
    plt.grid(b=True, which='major', color='k', linestyle='-')
    plt.suptitle('Bread Slingshot') # Set the graph title
    plt.legend(loc='upper right')   # Set the legend location
    plt.ylabel('Height (m)')        # Set the y-axis label
    plt.xlabel('Distance (m)')      # Set the x-axis label
    plt.ylim([-5,50])               # Set the y-axis limits
    plt.xlim([-5,120])              # Set the x-axis limits
    plt.show()                      # Make sure the graph appears

    print 'The bread lands at: {0:.3f}, {1:.3f}'.format(
        trajGuessX[-1], trajGuessY[-1])
    return