def measure_psf(vignet, pixscale=1., show=False, mask_value=None):
    y, x = np.mgrid[-vignet.shape[0]/2:vignet.shape[0]/2, -vignet.shape[1]/2:vignet.shape[1]/2]*pixscale
    if mask_value :
        vignet = ma.masked_values(vignet, mask_value).filled(0)
    # Fit the data using astropy.modeling
    p_init=models.Gaussian2D(amplitude=vignet.max(), x_mean=0., y_mean=0.,
        x_stddev=2*pixscale, y_stddev=2*pixscale, theta=0, cov_matrix=None)
    fit_p = fitting.LevMarLSQFitter()

    p = fit_p(p_init, x, y, vignet)
    barycenter=measure_barycenter(vignet, pixscale=pixscale)
    
    # Plot the data with the best-fit model
    P.figure(figsize=(8, 2.5))
    P.subplot(1, 3, 1)
    P.imshow(vignet, origin='lower', interpolation='nearest', vmin=vignet.min(), vmax=vignet.max())
    P.title("Data")
    P.subplot(1, 3, 2)
    P.imshow(p(x, y), origin='lower', interpolation='nearest', vmin=vignet.min(), vmax=vignet.max())
    P.scatter(vignet.shape[0]/2, vignet.shape[1]/2,marker="+")
    P.annotate("({:.3f},{:.3f})".format(*barycenter), (vignet.shape[0]/3, vignet.shape[1]/3))
    P.title("Model - psf = {:.2f}".format(2.3548*np.mean([p.x_stddev.value, p.y_stddev.value])))
    P.subplot(1, 3, 3)
    P.imshow(vignet - p(x, y), origin='lower', interpolation='nearest', vmin=-vignet.max()/10,vmax=vignet.max()/10)
    P.title("Residual")
    P.tight_layout()
    if show :
        P.show()
    
    return p

    def plot_graph(self):
        '''
        plots a matplotlib graph from the stocks data. Dates on the x axis and
        Closing Prices on the y axis. Then adds it to the graph_win in the
        display frame as a tk widget()
        '''
        x_axis = [
			dt.datetime.strptime(self.daily_data[day][0], '%Y-%m-%d')
            for day in range(1, len(self.daily_data) - 1)
		]
        y_axis = [
			self.daily_data[cls_adj][-1]
            for cls_adj in range(1, len(self.daily_data) - 1)
		]
        fig = plt.figure()
        ax = fig.add_subplot("111")
        ax.plot(
			x_axis,
            y_axis,
            marker='h',
            linestyle='-.',
            color='r',
            label='Daily Adjusted Closing Prices'
		)
        labels = ax.get_xticklabels()
        for label in labels:
            label.set_rotation(15)
        plt.xlabel('Dates')
        plt.ylabel('Close Adj')
        plt.legend()
        plt.tight_layout()  # adjusts the graph to fit in the space its limited to
        self.data_plot = FigureCanvasTkAgg(fig, master=self.display)
        self.data_plot.show()
        self.graph_win = self.data_plot.get_tk_widget()

    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(self,file):
        cds = CaseDataset(file, 'bson')
        data = cds.data.driver('driver').by_variable().fetch()
        cds2 = CaseDataset('../output/therm_mc_20141110173851.bson', 'bson')
        data2 = cds2.data.driver('driver').by_variable().fetch()
        
        #temp
        temp_boundary_k = data['hyperloop.temp_boundary']
        temp_boundary_k.extend(data2['hyperloop.temp_boundary'])
        temp_boundary = [((x-273.15)*1.8 + 32) for x in temp_boundary_k]
        #histogram
        n, bins, patches = plt.hist(temp_boundary, 100, normed=1, histtype='stepfilled')
        plt.setp(patches, 'facecolor', 'b', 'alpha', 0.75)

        #stats
        mean = np.average(temp_boundary)
        std = np.std(temp_boundary)
        percentile = np.percentile(temp_boundary,99.5)
        print "mean: ", mean, " std: ", std, " 99.5percentile: ", percentile
        x = np.linspace(50,170,150)
        plt.plot(x,mlab.normpdf(x,mean,std), color='black', lw=2)
        plt.xlim([60,160])
        plt.ylabel('Probability', fontsize=18)
        plt.xlabel(u'Equilibrium Temperature, \N{DEGREE SIGN}F', fontsize=18)
        #plt.show()
        plt.tight_layout()
        plt.savefig('../output/histo.pdf', dpi=300)

    def plot_abc(self, uarg, param):
        """Plot a() and b() functions on the same plot.

        """
        plt.figure(figsize=(8, 4))

        plt.subplot(1, 3, 1)
        plt.plot(uarg, self.afun(uarg, param))
        plt.axhline(0)
        plt.axvline(0)
        plt.ylabel('$a(u)$')
        plt.xlabel('$u$')

        plt.subplot(1, 3, 2)
        plt.plot(uarg, self.bfun(uarg, param))
        plt.axhline(0)
        plt.axvline(0)
        plt.ylabel('$b(u)$')
        plt.xlabel('$u$')

        plt.subplot(1, 3, 3)
        plt.plot(uarg, self.cfun(uarg, param))
        plt.axhline(0)
        plt.axvline(0)
        plt.ylabel('$c(u)$')
        plt.xlabel('$u$')

        plt.tight_layout()
        plt.show()

def zsview(im, cmap=pl.cm.gray, figsize=(8,5), contours=False, ccolor='r'):
    z1, z2 = zscale(im)
    pl.figure(figsize=figsize)
    pl.imshow(im, vmin=z1, vmax=z2, origin='lower', cmap=cmap, interpolation='none')
    if contours:
        pl.contour(im, levels=[z2], origin='lower', colors=ccolor)
    pl.tight_layout()

    def save(self, out_path):
        '''Saves a figure for the monitor
        
        Args:
            out_path: str
        '''
        
        plt.clf()
        np.set_printoptions(precision=4)
        font = {
            'size': 7
        }
        matplotlib.rc('font', **font)
        y = 2
        x = ((len(self.d) - 1) // y) + 1
        fig, axes = plt.subplots(y, x)
        fig.set_size_inches(20, 8)

        for j, (k, v) in enumerate(self.d.iteritems()):
            ax = axes[j // x, j % x]
            ax.plot(v, label=k)
            if k in self.d_valid.keys():
                ax.plot(self.d_valid[k], label=k + '(valid)')
            ax.set_title(k)
            ax.legend()

        plt.tight_layout()
        plt.savefig(out_path, facecolor=(1, 1, 1))
        plt.close()

def run_method_usage(methods,cases):
    methods = [m[0] for m in methods]
    # Bootstrap the percentage error bars:
    percents =[]
    for i in range(10000):
        nc = resample(cases)
        percents.append(100*np.sum(nc,axis=0)/len(nc))
    percents=np.array(percents)
    mean_percents = np.mean(percents,axis=0)
    std_percents = np.std(percents,axis=0)*1.96
    inds=np.argsort(mean_percents).tolist()
    inds.reverse()
    avg_usage = np.mean(mean_percents)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    x=np.arange(len(methods))
    ax.plot(x,[avg_usage]*len(methods),'-',color='0.25',lw=1,alpha=0.2)
    ax.bar(x, mean_percents[inds], 0.6, color=paired[0],linewidth=0,
           yerr=std_percents[inds],ecolor=paired[1])
    #ax.set_title('Method Occurrence')
    ax.set_ylabel('Occurrence %',fontsize=30)
    ax.set_xlabel('Method',fontsize=30)
    ax.set_xticks(np.arange(len(methods)))
    ax.set_xticklabels(np.array(methods)[inds],fontsize=8)
    fig.autofmt_xdate()
    fix_axes()
    plt.tight_layout()
    fig.savefig(figure_path+'method_occurrence.pdf', bbox_inches=0)
    fig.show()
    return inds,mean_percents[inds]

def plot(rho, u, uLB, tau, rho_history, zdjecia, image, nx, maxIter ):
#    plt.figure(figsize=(15,15))
#    plt.subplot(4, 1, 1)
#    plt.imshow(u[1,:,0:50],vmin=-uLB*.15, vmax=uLB*.15, interpolation='none')#,cmap=cm.seismic
#    plt.colorbar()
    plt.rcParams["figure.figsize"] = (15,15)
    plt.subplot(5, 1, 1)
    plt.imshow(sqrt(u[0]**2+u[1]**2),vmin=0, vmax=uLB*1.6)#,cmap=cm.seismic
    plt.colorbar()
    plt.title('tau = {:f}'.format(tau))      
    
    plt.subplot(5, 1, 2)
    plt.imshow(u[0,:,:30],  interpolation='none')#,cmap=cm.seismicvmax=uLB*1.6,
    plt.colorbar()
    plt.title('tau = {:f}'.format(tau))  
    
    plt.subplot(5, 1, 3)
    plt.imshow(rho, interpolation='none' )#,cmap=cm.seismic
    plt.title('rho')   
    
    plt.subplot(5, 1,4)
    plt.title(' history rho')
    plt.plot(linspace(0,len(rho_history),len(rho_history)),rho_history)
    plt.xlim([0,maxIter])   
    
    plt.subplot(5, 1,5)
    plt.title(' u0 middle develop')
    plt.plot(linspace(0,nx,len(u[0,20,:])), u[1,20,:])
    plt.tight_layout()                  
    
    plt.savefig(path.join(zdjecia,'f{0:06d}.png'.format(image)))
    plt.clf();
        

    def do_plot_extras(self, extra):
        """ Plot other observed quantities as a function of time.

        Parameters
        ----------
        extra: string
          One of the quantities available in system.extras
        """
        # import pyqtgraph as pg

        colors = 'bgrcmykw' # lets hope for less than 9 data-sets
        t = self.time

        # handle inexistent field
        if extra not in self.extras._fields:
          from shell_colors import red
          msg = red('ERROR: ') + 'The name "%s" is not available in extras.\n' % extra
          clogger.fatal(msg)
          return

        i = self.extras._fields.index(extra) # index corresponding to this quantity

        plt.figure()
        # p = pg.plot()
        plt.plot(t, self.extras[i], 'o', label=extra)
        
        plt.xlabel('Time [days]')
        plt.ylabel(extra + ' []')
        plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.05))
        plt.minorticks_on()
        plt.tight_layout()
        plt.show()

    def _plot(self, doFAP=None):
      """
        Create a plot.
      """
      xlabel = 'Period [d]'
      ylabel = 'Power'

      self.fig = plt.figure()
      self.ax = self.fig.add_subplot(1,1,1)
      self.ax.set_title("Normalized periodogram")
      self.ax.set_xlabel(xlabel)
      self.ax.set_ylabel(ylabel)
      self.ax.semilogx(1./self.freq, self.power, 'b-')
      # plot FAPs
      if doFAP is None: 
        pass
      elif doFAP is True: # do default FAPs of 10%, 1% and 0.1%
        pmin = 1./self.freq.min()
        pmax = 1./self.freq.max()
        plvl1 = self.powerLevel(0.1) # 10% FAP
        plvl2 = self.powerLevel(0.01) # 1% FAP
        plvl3 = self.powerLevel(0.001) # 0.1% FAP
        self.ax.semilogx([pmin, pmax],[plvl1, plvl1],'k-')
        self.ax.semilogx([pmin, pmax],[plvl2, plvl2],'k--')
        self.ax.semilogx([pmin, pmax],[plvl3, plvl3],'k:')

      plt.tight_layout()
      plt.show()

    def do_plot_obs(self):
        """ Plot the observed radial velocities as a function of time.
        Data from each file are color coded and labeled.
        """
        # import pyqtgraph as pg

        colors = 'bgrcmykw' # lets hope for less than 9 data-sets
        t, rv, err = self.time, self.vrad, self.error # temporaries
        
        plt.figure()
        # p = pg.plot()
        # plot each files' values
        for i, (fname, [n, nout]) in enumerate(sorted(self.provenance.iteritems())):
            m = n-nout # how many values are there after restriction
            
            # e = pg.ErrorBarItem(x=t[:m], y=rv[:m], \
            #                     height=err[:m], beam=0.5,\
            #                     pen=pg.mkPen(None))
                                # pen={'color': 0.8, 'width': 2})
            # p.addItem(e)
            # p.plot(t[:m], rv[:m], symbol='o')
            plt.errorbar(t[:m], rv[:m], yerr=err[:m], \
                         fmt='o'+colors[i], label=fname)
            t, rv, err = t[m:], rv[m:], err[m:]
        
        plt.xlabel('Time [days]')
        plt.ylabel('RV [km/s]')
        plt.legend()
        plt.tight_layout()
        plt.show()

    def save(self, plot_filename, max_col=2):
        fig = plt.figure()
        fig.suptitle(self.title, fontsize=8)

        if len(self.datastore) <= max_col:
            colsz = len(self.datastore)
            rowsz = 1
            for i, d in enumerate(self.datastore):
                if len(d) == 0:
                    continue  # empty space
                ax = fig.add_subplot(rowsz, colsz, i + 1)
                self._plot(i, ax)
        else:
            sz = len(self.datastore)
            colsz = max_col
            rowsz = sz / max_col
            rowsz = rowsz if sz % max_col == 0 else rowsz + 1
            for i, d in enumerate(self.datastore):
                if len(d) == 0:
                    continue  # empty space
                ax = fig.add_subplot(rowsz, colsz, i + 1)
                self._plot(i, ax)

        #plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
        plt.tight_layout()
        plt.savefig(plot_filename)

        return self

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 plot_ra(s1, s2, idxs=None, epsilon=0.25, fig=None):
    """Computes the RA plot of two groups of samples"""

    ## compute log2 values
    l1 = np.log2(s1 + epsilon)
    l2 = np.log2(s2 + epsilon)

    ## compute A and R
    r = l1 - l2
    a = (l1 + l2) * 0.5

    fig = pl.figure() if fig is None else fig
    pl.figure(fig.number)

    if idxs is None:
        pl.plot(a, r, '.k', markersize=2)
    else:
        pl.plot(a[~idxs], r[~idxs], '.k', markersize=2)
        pl.plot(a[idxs], r[idxs], '.r')

    pl.axhline(0, linestyle='--', color='k')

    pl.xlabel('(log2 sample1 + log2 sample2) / 2')
    pl.ylabel('log2 sample1 - log2 sample2')

    pl.tight_layout()

def plotRttDistribution(rttEstimates, ip, filename, nbBins=500, logscale=False):
    """Plot the RTT distribution of an IP address

    :rttEstimates: pandas DataFrame containing the RTT estimations
    :ip: IP address to plot
    :filename: Filename for the plot
    :nbBins: Number of bins in the histogram
    :logscale: Plot RTTs in logscale if set to True
    :returns: None

    """

    if logscale:
        data = np.log10(rttEstimates[rttEstimates.index == ip].rtt)
    else:
        data = rttEstimates[rttEstimates.index == ip].rtt

    h, b=np.histogram(data, nbBins, normed=True)
    plt.figure(1, figsize=(9, 3))
    plt.clf()
    ax = plt.subplot()
    x = b[:-1]
    ax.plot(x, h, "k")
    ax.grid(True)
    plt.title("%s (%s RTTs)" % (ip, len(data)))
    if logscale:
        plt.xlabel("log10(RTT)")
    else:
        plt.xlabel("RTT")
    plt.ylabel("pdf")
    minorLocator = mpl.ticker.MultipleLocator(10)
    ax.xaxis.set_minor_locator(minorLocator)
    plt.tight_layout()
    plt.savefig(filename)

def expt3():
	"""
	Experiment 2: Chooses the result files and generates figures
	"""
	result_file = "./expt3.txt"
	input_threads, input_sizes, throughputs, resp_times \
		= parse_output(result_file) 

	throughputs_MiB = [tp/2**20 for tp in throughputs]

	fig1, (ax0, ax1) = pl.subplots(ncols=2, figsize=(6, 3))
	fig1.set_tight_layout(True)

	ax0.plot(input_threads, throughputs_MiB, 
			'bo-', ms=MARKER_SIZE, mew=0, mec='b')
	ax0.set_xlabel("threads")
	ax0.set_ylabel("throughput (MiB/sec)")
	ax0.set_xlim(0,25)
	ax0.text(1, 85, "(A)")

	ax1.plot(input_threads, resp_times, 
			'mo-', ms=MARKER_SIZE, mew=0, mec='m')
	ax1.set_xlabel("threads")
	ax1.set_ylabel("response time (sec)")
	ax1.set_xlim(0,25)
	ax1.set_ylim(0,400)
	ax1.text(1, 375, "(B)")

	pl.tight_layout()
	pl.savefig("./figures/%s" % result_file.replace(".txt", ".pdf"))

def plot_ga_cnn(data, save=False, save_dest=None):
    plt.figure()

    # generation and all time best score plot
    # plt.subplot(211)
    # plt.title("GA CNN tuner on dataset D2")
    # plt.plot(data["dataset1"]["gen_best_score"], label="Generation Best")
    # plt.plot(data["dataset1"]["all_time_best_score"], label="All Time best")
    # plt.xlabel("Generation")
    # plt.ylabel("Score")
    # plt.xticks(range(0, 20))
    # plt.xlim([0, 9])
    # plt.ylim([0, 0.01])
    # plt.legend(loc=0)

    # plt.subplot(212)
    # plt.title("GA CNN tuner on dataset D3")
    plt.plot(data["gen_best_score"], label="Generation Best")
    plt.plot(data["all_time_best_score"], label="All Time best")
    plt.xlabel("Generation")
    plt.ylabel("Score")
    plt.xticks(range(0, 20))
    plt.xlim([0, 9])
    plt.ylim([0, 0.01])
    plt.legend(loc=0)

    # show plot or save as picture
    if save is False:
        plt.show()
    else:
        if save_dest is None:
            raise RuntimeError("save_dest not set!!")
        plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
        plt.savefig(save_dest, dpi=200)

def plmyfig(df, bgname, dirname, tar, count=10):
    #plot fig!
    print("Starting Plot %s %s" % (dirname, bgname))
    if len(df) > count:
        df = df.head(count)
    pos = plt.arange(len(df)) + 0.5
    ytick = _getTerm(df['Term_description'], df['Term_ID'], bgname)
    xs = [float(n) for n in df[' -log10(pvalue)']]
    ytick.reverse()
    xs.reverse()
    plt.barh(pos, xs, align = 'center', height = 0.5, alpha = 1, color='orange')
    plt.yticks(pos, ytick, size = 'x-small')
    plt.xlabel('$-Log10(pValue)$')
    plt.title('%s' % bgname)
    ax = plt.gca()
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    try:
        plt.tight_layout()
    except ValueError:
        pass
    filename = os.path.join(tar, dirname, dirname + '_' + bgname)
    plt.savefig(filename + '.png', dpi = 72)
    plt.savefig(filename + '.pdf')
    plt.close()

def test(args):
    data = multivariate_normal([0, 0], [[1, 2], [2, 5]], int(args[1]))
    print(data)
    # PCA
    result = pca(data, base_num=int(args[2]))
    pc_base = result[0]
    print(pc_base)

    # Plotting
    fig = plt.figure()
    fig.add_subplot(1, 1, 1)
    plt.axvline(x=0, color="#000000")
    plt.axhline(y=0, color="#000000")
    # Plot data
    plt.scatter(data[:, 0], data[:, 1])
    # Draw the 1st principal axis
    pc_line = sp.array([-3.0, 3.0]) * (pc_base[1] / pc_base[0])
    plt.arrow(0, 0, -pc_base[0] * 2, -pc_base[1] * 2, fc="r", width=0.15, head_width=0.45)
    plt.plot([-3, 3], pc_line, "r")
    # Settings
    plt.xticks(size=15)
    plt.yticks(size=15)
    plt.xlim([-3, 3])
    plt.tight_layout()
    plt.show()
    plt.savefig("image.png")

    return 0