def plot_beta_dist( ctr, trials, success, alphas, betas, turns ):
	"""
	Pass in the ctr, trials and success, alphas, betas returned
	by the `experiment` function and the number of turns 
	and plot the beta distribution for all the arms in that turn
	"""
	subplot_num = len(turns) / 2
	x = np.linspace( 0.001, .999, 200 )
	fig = plt.figure( figsize = ( 14, 7 ) ) 

	for idx, turn in enumerate(turns):

		plt.subplot( subplot_num, 2, idx + 1 )

		for i in range( len(ctr) ):
			y = beta( alphas[i] + success[ turn, i ], 
					  betas[i] + trials[ turn, i ] - success[ turn, i ] ).pdf(x)
			line = plt.plot( x, y, lw = 2, label = "arm {}".format( i + 1 ) )
			color = line[0].get_color()
			plt.fill_between( x, 0, y, alpha = 0.2, color = color )
			plt.axvline( x = ctr[i], color = color, linestyle = "--", lw = 2 )
			plt.title("Posteriors After {} turns".format(turn) )
			plt.legend( loc = "upper right" )

	return fig

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

def plot(most_probable_params, best_step, chain, lnprob, nwalkers, ndim, niter):
	# PLOT WALKERS
	fig=plt.figure()
	Y_PLOT=2
	X_PLOT=5
	alpha=0.1
	
	# Walkers
	burnin=0.5*niter
	print 'burnin', burnin
	for n in range(ndim):
		ax=fig.add_subplot(Y_PLOT+1,X_PLOT,n+1)
		for i in range(nwalkers): # Risem za posamezne walkerje
			d=np.array(chain[i][burnin:,n])
			ax.plot(d, color='black', alpha=alpha)
		ax.set_xlabel(thetaText[n+1])
		if best_step-burnin>0:
			plt.axvline(x=best_step-burnin, linewidth=1, color='red')
		
	# Lnprob
	ax=fig.add_subplot(Y_PLOT+1,1,Y_PLOT+1)
	for i in range(nwalkers):
		ax.plot((lnprob[i][burnin:]), color='black', alpha=alpha)
	if best_step-burnin>0:
		plt.axvline(x=best_step-burnin, linewidth=1, color='red')
	
	#~ import triangle
	#~ fig = triangle.corner(sampler.flatchain, truths=most_probable_params) # labels=thetaText
	
	plt.show()

def plot_spikes(time,voltage,APTimes,titlestr):
    """
    plot_spikes takes four arguments - the recording time array, the voltage
    array, the time of the detected action potentials, and the title of your
    plot.  The function creates a labeled plot showing the raw voltage signal
    and indicating the location of detected spikes with red tick marks (|)
    """
# Make a plot and markup
    plt.figure()
    plt.title(titlestr)
    plt.xlabel("Time (s)")
    plt.ylabel("Voltage (uV)") 

    plt.plot(time, voltage)
    
# Vertical positions for red marker
# The following attributes are configurable if required    
    vertical_markers_indent = 0.01 # 1% of Voltage scale height
    vertical_markers_height = 0.03 # 5% of Voltage scale height
    y_scale_height = 100 # Max of scale
    
    marker_ymin = 0.5 + ( max(voltage) / y_scale_height / 2 ) + vertical_markers_indent
    marker_ymax = marker_ymin + vertical_markers_height

# Drawing red markers for detected spikes
    for spike in APTimes:
        plt.axvline(spike, ymin=marker_ymin, ymax=marker_ymax, color='red')
    
    plt.draw()

 def segmentation(self, threshold):
     img = self.spectrogram["data"]
     mask = (img > threshold).astype(np.float)
     hist, bin_edges = np.histogram(img, bins=60)
     bin_centers = 0.5*(bin_edges[:-1] + bin_edges[1:])
     binary_img = mask > 0.5
     plt.figure(figsize=(11,8))
     plt.subplot(131)
     plt.imshow(img)
     plt.axis('off')
     plt.subplot(132)
     plt.plot(bin_centers, hist, lw=2)
     print(threshold)
     plt.axvline(threshold, color='r', ls='--', lw=2)
     plt.text(0.57, 0.8, 'histogram', fontsize=20, transform = plt.gca().transAxes)
     plt.text(0.45, 0.75, 'threshold = '+ str(threshold)[0:5], fontsize=15, transform = plt.gca().transAxes)
     plt.yticks([])
     plt.subplot(133)     
     plt.imshow(binary_img)
     plt.axis('off')
     plt.subplots_adjust(wspace=0.02, hspace=0.3, top=1, bottom=0.1, left=0, right=1)
     plt.show()
     print(img.max())
     print(binary_img.max())
     
     return mask

def showKernel(dataOrMatrix, fileName = None, useLabels = True, **args) :
 
    labels = None
    if hasattr(dataOrMatrix, 'type') and dataOrMatrix.type == 'dataset' :
	data = dataOrMatrix
	k = data.getKernelMatrix()
	labels = data.labels
    else :
	k = dataOrMatrix
	if 'labels' in args :
	    labels = args['labels']

    import matplotlib

    if fileName is not None and fileName.find('.eps') > 0 :
        matplotlib.use('PS')
    from matplotlib import pylab

    pylab.matshow(k)
    #pylab.show()

    if useLabels and labels.L is not None :
	numPatterns = 0
	for i in range(labels.numClasses) :
	    numPatterns += labels.classSize[i]
	    #pylab.figtext(0.05, float(numPatterns) / len(labels), labels.classLabels[i])
	    #pylab.figtext(float(numPatterns) / len(labels), 0.05, labels.classLabels[i])
	    pylab.axhline(numPatterns, color = 'black', linewidth = 1)
	    pylab.axvline(numPatterns, color = 'black', linewidth = 1)
    pylab.axis([0, len(labels), 0, len(labels)])
    if fileName is not None :
        pylab.savefig(fileName)
	pylab.close()

def mark_cross(center, **kwargs):
    """Mark a cross. Correct for matplotlib imshow funny coordinate system.
    """
    N = 20
    plt.hold(1)
    plt.axhline(y=center[1]-0.5, **kwargs)
    plt.axvline(x=center[0]-0.5, **kwargs)

def my_lines(ax, pos, *args, **kwargs):
    if ax == 'x':
        for p in pos:
            plt.axvline(p, *args, **kwargs)
    else:
        for p in pos:
            plt.axhline(p, *args, **kwargs)

 def plot_vrad_chi2(self, fig):
     fig.clf()
     ax = fig.add_subplot(111)
     ax.plot(self.v_rads, self.v_rad_grid)
     ax.set_xlabel('v_rad [km/s]')
     ax.set_ylabel('Chi2')
     ax.axvline(self.v_rad, color='red', linewidth=2)
     ax.set_title('v_rad = %.2f' % self.v_rad)

def show_hist(datai, mean1, mean2, avg, title="Histogram", bins=100):
    # mean1,mean2.min(),npdata.max()
    # avg=(mean1+mean2)/2
    fig, ax = subplots(1, 1, sharex=True, sharey=True, figsize=(10, 10))
    n, bins, patches = ax.hist(datai.flat, bins=bins, range=(mean1, mean2), histtype="bar")
    axvline(x=avg, alpha=0.7, linewidth=3, color="r")
    ax.set_title("histogram")
    show()

 def plotline(x):
     ax = x['ax']
     xk = x['xk']
     plt.sca(ax)
     trans = blended_transform_factory(ax.transData,ax.transAxes)
     plt.axvline(smpar[xk],ls='--')
     plt.text(smpar[xk],0.9,'SM',transform=trans)
     
     if libpar is not None:
         plt.axvline(libpar[xk])
         plt.text(libpar[xk],0.9,'LIB',transform=trans)

	def tf_analysis(self, plot_Z = True, frequencies = None, vis_frequency_limits = [1.8, 2.2], nr_cycles = 16, analysis_sample_rate = 100):
		self.assert_data_intern()

		if frequencies == None:
			frequencies = np.linspace(1.0, self.analyzer.low_pass_pupil_f, 40)

		down_sample_factor = int(self.sample_rate/analysis_sample_rate)
		resampled_signal = self.pupil_bp_pt[:,::down_sample_factor]

		# complex tf results per trial
		self.tf_trials = mne.time_frequency.cwt_morlet(resampled_signal, analysis_sample_rate, frequencies, use_fft=True, n_cycles=nr_cycles, zero_mean=True)
		self.instant_power_trials = np.abs(self.tf_trials)

		# z-score power
		self.Z_tf_trials = np.zeros_like(self.instant_power_trials)
		m = self.instant_power_trials.mean(axis = -1)
		sd = self.instant_power_trials.std(axis = -1)

		for z in range(len(self.Z_tf_trials)):
			self.Z_tf_trials[z] = ((self.instant_power_trials[z].T - m[z]) / sd[z] ).T

		# some obvious conditions
		if plot_Z:
			tf_to_plot = self.Z_tf_trials
		else:
			tf_to_plot = self.instant_power_trials

		f = pl.figure(figsize = (24,24))
		for x in range(len(self.trial_indices)):
			s = f.add_subplot(len(self.trial_indices), 2, (x*2)+1)
			pl.imshow(np.squeeze(tf_to_plot[x,(frequencies > vis_frequency_limits[0]) & (frequencies < vis_frequency_limits[1]),::100]), cmap = 'seismic', extent = [self.from_zero_timepoints[0], self.from_zero_timepoints[-1], vis_frequency_limits[-1], vis_frequency_limits[0]], aspect='auto')
			sn.despine(offset=10)
			s = f.add_subplot(len(self.trial_indices), 2, (x*2)+2)
			# pl.imshow(np.squeeze(tf_to_plot[x,:,::100]), cmap = 'gray')
			pl.plot(self.from_zero_timepoints[::down_sample_factor], np.squeeze(np.squeeze(tf_to_plot[x,(frequencies > vis_frequency_limits[0]) & (frequencies < vis_frequency_limits[1]),:])).mean(axis = 0), 'k')
			if len(self.events) != 0:
				events_this_trial = self.events[(self.events['EL_timestamp'] > self.timestamps_pt[x][0]) & (self.events['EL_timestamp'] < self.timestamps_pt[x][-1])]
				for sc, scancode in enumerate(self.scancode_list):
					these_event_times = events_this_trial[events_this_trial['scancode'] == scancode]['EL_timestamp']
					for tet in these_event_times:
						pl.axvline(x = (tet - self.timestamps_pt[x,0]) / self.sample_rate, c = self.colors[sc], lw = 5.0)
			sn.despine(offset=10)
		pl.tight_layout()
		pl.savefig(os.path.join(self.analyzer.fig_dir, self.file_alias + '_%i_tfr.pdf'%nr_cycles))

		with pd.get_store(self.analyzer.h5_file) as h5_file: 
			for name, data in zip(['tf_complex_real', 'tf_complex_imag', 'tf_power', 'tf_power_Z'], 
				np.array([np.real(self.tf_trials), np.imag(self.tf_trials), self.instant_power_trials, self.Z_tf_trials], dtype = np.float64)):
				opd = pd.Panel(data, 
								items = pd.Series(self.trial_indices), 
								major_axis = pd.Series(frequencies), 
								minor_axis = self.from_zero_timepoints[::down_sample_factor])
				h5_file.put("/%s/tf/cycles_%s_%s"%(self.file_alias, nr_cycles, name), opd)

def epi_vs_gain_volcano_plot(filtered_gain_snps, filtered_epi_snps, gain_vals, epi_vals, max_p, min_I3):
    gain_I3 = []
    gain_log_p = []
    for snps in filtered_gain_snps:
        gain_I3.append(gain_vals[snps])

        order = switch_snp_key_order(snps)
        if epi_vals.has_key(order[0]):
            gain_log_p.append(epi_vals[order[0]])
        elif epi_vals.has_key(order[1]):
            gain_log_p.append(epi_vals[order[1]])
    gain_log_p = -1 * np.log10(gain_log_p)

    epi_I3 = []
    epi_log_p = []
    for snps in filtered_epi_snps:
        order = switch_snp_key_order(snps)
        if gain_vals.has_key(order[0]):
            epi_I3.append(gain_vals[order[0]])
        elif gain_vals.has_key(order[1]):
            epi_I3.append(gain_vals[order[1]])

        epi_log_p.append(epi_vals[snps])
    epi_log_p = -1 * np.log10(epi_log_p)

    mp.figure(1)
    mp.xlabel("I3")
    mp.ylabel("-log10(P)")
    mp.title("Volcano plot - EPISTASIS and GAIN")
    mp.plot(epi_I3, epi_log_p, "bo")
    mp.plot(gain_I3, gain_log_p, "ro")
    mp.axhline(y=(-1 * np.log10(max_p)), linewidth=2, color="g")
    mp.axvline(x=min_I3, linewidth=2, color="g")
    # label max point
    max_x = np.max(gain_I3)
    max_y = np.max(gain_log_p)
    best_connection = ""
    # label best edge
    for snps in epi_vals:
        if -1 * np.log10(epi_vals[snps]) == max_y:
            best_connection = str(snps)
    mp.text(
        np.max(gain_I3),
        np.max(gain_log_p),
        best_connection,
        fontsize=10,
        horizontalalignment="center",
        verticalalignment="center",
    )
    mp.show()

    print

    def lcReview(self, fibre, eventList=[], save=False, saveName='default.png'):
        """
        """
        time, flux, xpos, ypos, apts, msky, qual = self.getFullData(fibre)
        
        TA = TimingAnalysis()
        result, rms = TA.standardscore(time, flux, 2.0)
        wmflux = result[0]
        mSNR = np.median(wmflux/rms)

        plt.figure(100, figsize=(8,12))
        plt.clf()
        plt.subplots_adjust(hspace=0.2)
        
        plt.subplot(4,1,1)
        plt.plot(time, flux, 'k+-', drawstyle='steps-mid')
        plt.plot(time, wmflux, 'r-', label='averaged')
        if len(eventList) == 0:
            pass
        else:
            for value in eventList:
                plt.axvline(x=value, color='red', ls='--')
        plt.xlim(time[0], time[-1])
        plt.ylim(0, np.median(flux) + 0.5*np.median(flux))
        plt.ylabel('Flux (count)')
        plt.xticks(visible=False)
        plt.legend()
        plt.title('{}, fibre {}. mSNR: {:.2f}'.format(self.filename[:-5], fibre, mSNR))
        
        plt.subplot(4,1,2, sharex=plt.subplot(4,1,1))
        plt.plot(time, xpos, 'b+', time, ypos, 'r+')
        plt.xlim(time[0], time[-1])
        plt.ylim(5,35)
        plt.xticks(visible=False)
        plt.ylabel('Position (pixel)')
        
        plt.subplot(4,1,3, sharex=plt.subplot(4,1,1))
        plt.plot(time, apts, 'k+')
        plt.xlim(time[0], time[-1])
        plt.ylim(0, np.mean(apts) + 0.5*np.mean(apts))
        plt.xticks(visible=False)
        plt.ylabel('Apt Size (pixel)')
        
        plt.subplot(4,1,4, sharex=plt.subplot(4,1,1))
        plt.plot(time, msky, 'k+')
        plt.xlim(time[0], time[-1])
        plt.ylim(0, np.median(msky) + 0.5*np.median(msky))
        plt.xlabel('Time (s)')
        plt.ylabel('Mean Sky Value (count)')
        
        plt.show()

def plot_ccf(w,tspec,mspec):
    lag,tmccf = ccf(tspec,mspec)

    dv = restwav.loglambda_wls_to_dv(w) 
    dv = lag*dv

    dvmax = ccs.findpeak(dv,tmccf)
    plt.plot(dv,tmccf,'k',label='tspec')    
    plt.axvline(dvmax,color='RoyalBlue',lw=2,alpha=0.4,zorder=0)

    AddAnchored("dv (max) = %.2f km/s" % dvmax,**annkw)
    plt.xlim(-50,50)
    plt.xlabel('dv (km/s)')

    return dvmax

    def plotGeneral(self, initFlux, normFlux, xpos, ypos, quality, direction, variance, threshold):
        plt.figure('General')
        plt.clf()
        
        plt.subplot(5,1,1)
        plt.plot(self.time, initFlux, 'k-', drawstyle='steps-mid')
        for event in self.eventList:
            plt.axvline(x=self.time[event], color='red', ls='--')
        plt.xlim(self.time[0], self.time[-1])
        plt.ylabel('Original LCs')
        plt.xticks(visible=False)
        plt.title('{}, f{}'.format(self.fileName, self.fibre))
        
        plt.subplot(5,1,2, sharex=plt.subplot(5,1,1))
        plt.plot(self.time, normFlux, 'g+')
        plt.plot(self.time, normFlux, 'k-', drawstyle='steps-mid')
        for event in self.eventList:
            plt.axvline(x=self.time[event], color='red', ls='--')
        plt.xlim(self.time[0], self.time[-1])
        plt.ylim(0.0, 1.5)
        plt.ylabel('Normalised LCs')
        plt.xticks(visible=False)
    
        plt.subplot(5,1,3, sharex=plt.subplot(5,1,1))
        plt.plot(self.time, quality, 'k-', drawstyle='steps-mid', label='quality')
        plt.plot(self.time, direction, 'r-', drawstyle='steps-mid', label='direction')
        plt.xlim(self.time[0], self.time[-1])
        plt.ylim(-2, 2)
        plt.xticks(visible=False)
        plt.legend()
        
        plt.subplot(5,1,4, sharex=plt.subplot(5,1,1))
        plt.plot(self.time, xpos, 'b-', drawstyle='steps-mid', label='x')
        plt.plot(self.time, ypos, 'r-', drawstyle='steps-mid', label='y')
        plt.xlim(self.time[0], self.time[-1])
        plt.ylim(1,40)
        plt.xticks(visible=False)
        plt.legend()
    
        plt.subplot(5,1,5, sharex=plt.subplot(5,1,1))
        plt.plot(self.time, variance, 'k-', drawstyle='steps-mid')
        plt.plot(self.time, threshold, 'r-', label=r'{}$\sigma$'.format(self.dTL))
        plt.xlim(self.time[0], self.time[-1])
        plt.ylabel('Variance')
        plt.xlabel('Time (s)')
        plt.legend()

        plt.show()

	def raw_signal_plot(self):
		self.assert_data_intern()

		f = pl.figure(figsize = (24,24))
		for x in range(len(self.trial_indices)):
			s = f.add_subplot(len(self.trial_indices), 1, x+1)
			pl.plot(self.timestamps_pt[x][::100], self.pupil_bp_pt[x][::100], 'k')
			if len(self.events) != 0:
				events_this_trial = self.events[(self.events['EL_timestamp'] > self.timestamps_pt[x][0]) & (self.events['EL_timestamp'] < self.timestamps_pt[x][-1])]
				for sc, scancode in enumerate(self.scancode_list):
					these_event_times = events_this_trial[events_this_trial['scancode'] == scancode]['EL_timestamp']
					for tet in these_event_times:
						pl.axvline(x = tet, c = self.colors[sc], lw = 5.0)
			sn.despine(offset=10)
		pl.tight_layout()
		pl.savefig(os.path.join(self.analyzer.fig_dir, self.file_alias + '_raw.pdf'))

def plot(model, x_obs, y_obs, argmin_a_x, a, x):
    y_mean, y_sigma = model.predict(x, return_std=True)
    y_sigma = np.atleast_2d(y_sigma).T
    a_x = a(x)
    plt.plot(x, y_mean, 'r-')
    plt.plot(x_obs, y_obs, 'o')

    plt.axvline(argmin_a_x)

    plt.plot(
        x,
        a_x*(1/np.max(a_x)),
        'g--',
    )

    plt.plot(x, f(x), 'm--')
    plot_confidences(x, y_mean, y_sigma, confidences=[1])

def plotDifference(bst, dtest, target = None):
    p = bst.predict(dtest)
    l = 0
    if (target is None):    
        l = dtest.get_label()
    else:
        l = target
    res = np.abs((l-p)/l)
    res = res[~np.isnan(res)]
    res = res[~np.isinf(res)]
    xmax = 0.2
    for y in np.arange(0,xmax,xmax/8):    
        plt.axvline(y, color='g', linestyle='dashed')
    plt.hist(res, range=(0,xmax), bins=50, color='red')
    plt.axvline(res.mean(), color='c', linestyle='dashed', linewidth=2)
    for y in range(200,1600,200):    
        plt.axhline(y, color='g', linestyle='dashed')

def rDist(r,sd1 = 10, sd2 = 10):
	mup = 50 + (100*(r/2))
	mlow = 50 - (100*(r/2))
	y2 = np.random.normal(mup,sd1,10000)
	y1 = np.random.normal(mlow,sd1,10000)
	plt.grid(True)
	plt.hist(y1,50, label = '$Control\ Group$', alpha = .7, color = 'b')
	plt.hist(y2,50, label = '$Experimental\ Group$', alpha = .5, color = 'g')
	plt.xlim(0, 100)
	plt.legend(	loc= 'right')
	plt.title('$Binomial\ Effect\ Size\ Display\ (BESD).$')
	plt.xlabel('Bin number for each one percentage bin')
	plt.ylabel('Frequency of each value occurring')
	plt.axvline(mup, ms = 3, color= 'r')
	plt.axvline(mlow, ms = 3, color = 'r')
	plt.show()
	return