function to plot series of chi measurements as a function of temperature, holding field constant and varying frequency
195020
def plot_xtf(fignum, XTF, Fs, e, b):
    
    plt.figure(num=fignum)
    plt.xlabel('Temperature (K)')
    plt.ylabel('Susceptibility (m^3/kg)')
    k = 0
    Flab = []
    for freq in XTF:
        T, X = [], []
        for xt in freq:
            X.append(xt[0])
            T.append(xt[1])
        plt.plot(T, X)
        plt.text(T[-1], X[-1], str(int(Fs[k])) + ' Hz')
#        Flab.append(str(int(Fs[k]))+' Hz')
        k += 1
    plt.title(e + ': B = ' + '%8.1e' % (b) + ' T')
704898
def to_abivars(self):
        """"""
        return {
                # Spectral function
                "nomegasf": self.nomegasf,
                "domegasf": self.domegasf,
                "spmeth"  : self.spmeth,
                 # Frequency mesh for the polarizability
                "nfreqre"  : self.nfreqre,
                "freqremax": self.freqremax,
                "nfreqim"  : self.nfreqim,
                "freqremin": self.freqremin,
                }
1068924
def covariance_plot(ax, Si, param_dict, unit=""):
    '''

    '''

    if Si['sigma'] is not None:
        # sigma is not present if using morris groups
        y = Si['sigma']
        out = ax.scatter(Si['mu_star'], y, c=u'k', marker=u'o',
                         **param_dict)
        ax.set_ylabel(r'$\sigma$')

        ax.set_xlim(0,)
        ax.set_ylim(0,)

        x_axis_bounds = np.array(ax.get_xlim())

        line1, = ax.plot(x_axis_bounds, x_axis_bounds, 'k-')
        line2, = ax.plot(x_axis_bounds, 0.5 * x_axis_bounds, 'k--')
        line3, = ax.plot(x_axis_bounds, 0.1 * x_axis_bounds, 'k-.')

        ax.legend((line1, line2, line3), (r'$\sigma / \mu^{\star} = 1.0$',
                                          r'$\sigma / \mu^{\star} = 0.5$',
                                          r'$\sigma / \mu^{\star} = 0.1$'),
                  loc='best')

    else:
        y = Si['mu_star_conf']
        out = ax.scatter(Si['mu_star'], y, c=u'k', marker=u'o',
                         **param_dict)
        ax.set_ylabel(r'$95\% CI$')

    ax.set_xlabel(r'$\mu^\star$ ' + unit)
    ax.set_ylim(0-(0.01 * np.array(ax.get_ylim()[1])), )

    return out
1307427
def fmt_chi2(f):
    
    if hasattr(f, 'logGBF'):
        fmt = "chi2/dof = %.2g [%d]    Q = %.2g    log(GBF) = %.5g"
        chi2_dof = f.chi2 / f.dof if f.dof != 0 else 0
        return fmt % (chi2_dof, f.dof, f.Q, f.logGBF)
    else:
        fmt = "chi2/dof = %.2g [%d]    Q = %.2g"
        chi2_dof = f.chi2 / f.dof if f.dof != 0 else 0
        return fmt % (chi2_dof, f.dof, f.Q)
839179
def meval(self, f, *args, **kw):
        
        # !! PART OF ORIGINAL DOCSTRING INCOMPATIBLE WITH CLASS INTERFACE !!
        # Example
        # -------
        # We set up multiple stationary models, one for a reference (ground)
        # state, and two for biased states, and group them in a
        # MultiStationaryModel.
        # >>> from pyemma.thermo import StationaryModel, MEMM
        # >>> m_1 = StationaryModel(f=[1.0, 0], label='biased 1')
        # >>> m_2 = StationaryModel(f=[2.0, 0], label='biased 2')
        # >>> m_mult = MEMM([m_1, m_2], [0, 0], label='unbiased')
        # Compute the stationary distribution for the two biased models
        # >>> m_mult.meval('stationary_distribution')
        # [array([ 0.73105858,  0.26894142]), array([ 0.88079708,  0.11920292])]
        # We set up multiple Markov state models for different temperatures
        # and group them in a MultiStationaryModel.
        # >>> import numpy as np
        # >>> from pyemma.msm import MSM
        # >>> from pyemma.thermo import MEMM
        # >>> b = 20  # transition barrier in kJ / mol
        # >>> temps = np.arange(300, 500, 25)  # temperatures 300 to 500 K
        # >>> p_trans = [np.exp(- b / kT) for kT in 0.00831*temps ]
        # >>> # build MSMs for different temperatures
        # >>> msms = [MSM(P=np.array([[1.0-p, p], [p, 1.0-p]])) for p in p_trans]
        # >>> # build Multi-MSM
        # >>> msm_mult = MEMM(pi=msms[0].stationary_distribution, label='300 K', models=msms)
        # Compute the timescales and see how they decay with temperature
        # Greetings to Arrhenius.
        # >>> np.hstack(msm_mult.meval('timescales'))
        # array([ 1523.83827932,   821.88040004,   484.06386176,   305.87880068,
        #          204.64109413,   143.49286817,   104.62539128,    78.83331598])
        # !! END OF INCOMPATIBLE PART !!
        return [_call_member(M, f, *args, **kw) for M in self.models]
254626
def filter_butter_plot_freqresponse(b, a, samplerate, numpoints=2000):
    # type: (Seq[float], Seq[float], int, int) -> None
    
    import matplotlib.pyplot as plt
    plt.figure(1)
    plt.clf()
    w, h = signal.freqz(b, a, worN=numpoints)
    plt.semilogx((samplerate*0.5/np.pi) * w, np.abs(h))
    plt.xlabel('Freq (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.show()
686008
def _chi2(self, theta, period, tmpid, return_gradient=False):
        
        template = self.templates[tmpid]
        phase = (self.t / period - theta[2]) % 1
        model = theta[0] + theta[1] * template(phase)
        chi2 = (((model - self.y) / self.dy) ** 2).sum()

        if return_gradient:
            grad = 2 * (model - self.y) / self.dy ** 2
            gradient = np.array([np.sum(grad),
                                 np.sum(grad * template(phase)),
                                 -np.sum(grad * theta[1]
                                         * template.derivative(1)(phase))])
            return chi2, gradient
        else:
            return chi2
963691
def _read_frequency_data(self):
        
        temp_dict = read_pattern(
            self.text, {
                "frequencies":
                r"\s*Frequency:\s+([\d\-\.]+)(?:\s+([\d\-\.]+)(?:\s+([\d\-\.]+))*)*",
                "IR_intens":
                r"\s*IR Intens:\s+([\d\-\.]+)(?:\s+([\d\-\.]+)(?:\s+([\d\-\.]+))*)*",
                "IR_active":
                r"\s*IR Active:\s+([YESNO]+)(?:\s+([YESNO]+)(?:\s+([YESNO]+))*)*",
                "ZPE":
                r"\s*Zero point vibrational energy:\s+([\d\-\.]+)\s+kcal/mol",
                "trans_enthalpy":
                r"\s*Translational Enthalpy:\s+([\d\-\.]+)\s+kcal/mol",
                "rot_enthalpy":
                r"\s*Rotational Enthalpy:\s+([\d\-\.]+)\s+kcal/mol",
                "vib_enthalpy":
                r"\s*Vibrational Enthalpy:\s+([\d\-\.]+)\s+kcal/mol",
                "gas_constant":
                r"\s*gas constant \(RT\):\s+([\d\-\.]+)\s+kcal/mol",
                "trans_entropy":
                r"\s*Translational Entropy:\s+([\d\-\.]+)\s+cal/mol\.K",
                "rot_entropy":
                r"\s*Rotational Entropy:\s+([\d\-\.]+)\s+cal/mol\.K",
                "vib_entropy":
                r"\s*Vibrational Entropy:\s+([\d\-\.]+)\s+cal/mol\.K",
                "total_enthalpy":
                r"\s*Total Enthalpy:\s+([\d\-\.]+)\s+kcal/mol",
                "total_entropy":
                r"\s*Total Entropy:\s+([\d\-\.]+)\s+cal/mol\.K"
            })

        keys = ["ZPE", "trans_enthalpy", "rot_enthalpy", "vib_enthalpy", "gas_constant", "trans_entropy", "rot_entropy", "vib_entropy", "total_enthalpy", "total_entropy"]

        for key in keys:
            if temp_dict.get(key) == None:
                self.data[key] = None
            else:
                self.data[key] = float(temp_dict.get(key)[0][0])

        if temp_dict.get('frequencies') == None:
            self.data['frequencies'] = None
            self.data['IR_intens'] = None
            self.data['IR_active'] = None
        else:
            temp_freqs = [
                value for entry in temp_dict.get('frequencies')
                for value in entry
            ]
            temp_intens = [
                value for entry in temp_dict.get('IR_intens')
                for value in entry
            ]
            active = [
                value for entry in temp_dict.get('IR_active')
                for value in entry
            ]
            self.data['IR_active'] = active

            freqs = np.zeros(len(temp_freqs) - temp_freqs.count('None'))
            for ii, entry in enumerate(temp_freqs):
                if entry != 'None':
                    freqs[ii] = float(entry)
            self.data['frequencies'] = freqs

            intens = np.zeros(len(temp_intens) - temp_intens.count('None'))
            for ii, entry in enumerate(temp_intens):
                if entry != 'None':
                    intens[ii] = float(entry)
            self.data['IR_intens'] = intens

            header_pattern = r"\s*Raman Active:\s+[YESNO]+\s+(?:[YESNO]+\s+)*X\s+Y\s+Z\s+(?:X\s+Y\s+Z\s+)*"
            table_pattern = r"\s*[a-zA-Z][a-zA-Z\s]\s*([\d\-\.]+)\s*([\d\-\.]+)\s*([\d\-\.]+)\s*(?:([\d\-\.]+)\s*([\d\-\.]+)\s*([\d\-\.]+)\s*(?:([\d\-\.]+)\s*([\d\-\.]+)\s*([\d\-\.]+))*)*"
            footer_pattern = r"TransDip\s+[\d\-\.]+\s*[\d\-\.]+\s*[\d\-\.]+\s*(?:[\d\-\.]+\s*[\d\-\.]+\s*[\d\-\.]+\s*)*"
            temp_freq_mode_vecs = read_table_pattern(
                self.text, header_pattern, table_pattern, footer_pattern)
            freq_mode_vecs = np.zeros(
                shape=(len(freqs), len(temp_freq_mode_vecs[0]), 3))

            for ii, triple_FMV in enumerate(temp_freq_mode_vecs):
                for jj, line in enumerate(triple_FMV):
                    for kk, entry in enumerate(line):
                        if entry != 'None':
                            freq_mode_vecs[int(ii * 3 + math.floor(kk / 3)),
                                           jj, kk % 3] = float(entry)

            self.data["frequency_mode_vectors"] = freq_mode_vecs
828166
def chi2_adaptive_binning(features_0,features_1,number_of_splits_list,systematics_fraction=0.0,title = "title", name="name", PLOT = True, DEBUG = False, transform='StandardScalar'):
	""""""

	max_number_of_splits = np.max(number_of_splits_list)
	#determine how many data points are in each sample
	no_0=features_0.shape[0]
	no_1=features_1.shape[0]

	print("features_0.shape : ", features_0.shape)
	no_dim = features_0.shape[1]
	#Give all samples in file 0 the label 0 and in file 1 the feature 1
	label_0=np.zeros((no_0,1))
	label_1=np.ones((no_1,1))

	#Create an array containing samples and features.
	data_0=np.c_[features_0,label_0]
	data_1=np.c_[features_1,label_1]

	features= np.r_[features_0,features_1]
	labels=    np.r_[label_0,   label_1]

	data=np.r_[data_0,data_1]
	data_same=np.c_[features,labels]

	#print("data : ",data)
	#print("data_same : ", data_same)
	#print("np.sum(data!=data_same) : ",np.sum(data!=data_same))
	assert np.sum(data!=data_same)==0
	assert (no_dim == data.shape[1]-1)

	if no_dim==2:
		plt.scatter(features[:,0],features[:,1], 0.1)
		plt.savefig('test.png')
		plt.clf()

	if transform=='StandardScalar':
		features = preprocessing.scale(features)
		data = np.c_[features,labels]

        if transform=='uniform':
		#data_new2 = data[:,0]
                data_new  = norm_highD_searchsorted(data[:,0])
                for D in range(1,no_dim):
                        temp = norm_highD_searchsorted(data[:,D])
                        data_new = np.c_[data_new,temp]
			#data_new2= np.c_[data_new2,data[:,D]]
		data_new = np.c_[data_new, np.r_[label_0,label_1]]
		#data_new2= np.c_[data_new2,np.r_[label_0,label_1]]

		print("data : ", data)
		data = data_new
		print("data new : ", data)
		#print("data_new2 : ", data_new2)
		#print("np.sum(data!=data_new2) : ",np.sum(data!=data_new2))

	
	np.random.shuffle(data)

	assert (no_dim == data.shape[1]-1)
	labels=data[:,-1]

	X_values= data[:,:-1]
	X_max   = np.amax(data,axis=0)[:-1] 
	X_min   = np.amin(data,axis=0)[:-1]
	X_total_width = (np.subtract(X_max,X_min))
	del data

	if transform=='fill01':
		#Scaling	
		X_values = X_values - X_min[None,:]
		X_values = X_values / X_total_width[None,:]
	
	if True:
		X_min = [0.]*no_dim
		X_total_width = [1.]*no_dim
	#b = X_values[:,0]
	#print("b[b[:]>2].shape[0] : \n", b[b[:]>2].shape[0] )	
	data = np.concatenate((X_values, labels[:,None]), axis=1)

	if no_dim==2:
		plt.scatter(data[:,0],data[:,1],0.1)
		plt.savefig('test_scaled.png')

	#print("X_values.shape : ",X_values.shape)
	starting_boundary = []
	for i in range(no_dim):
		starting_boundary.append([0.0,1.0])
	#Each key has the following stricture: # of splits and for each split if it was closer (a) or further away from (b) the origin. The original bin is "0"
	#For example "2ab" means it is the bin that was closer to the origin for the first split and further away for the second one.
	bin_boundaries_dict = {'0' : np.array(starting_boundary)}
	bin_points_dict = {'0' : data}

	for split_number in range(1,1+max_number_of_splits):
		for bin_key, bin_boundary in bin_boundaries_dict.items():
			if str(split_number-1) in bin_key:
				variances= np.var(bin_points_dict[bin_key][:,:-1], axis=0)
				#print("\nvariances : ", variances)
				dim_to_be_sliced = np.argmax(variances)
				#print("dim_to_be_sliced : ",dim_to_be_sliced)
				#print("bin_points_dict[bin_key] : ",bin_points_dict[bin_key])
				#print("bin_points_dict[bin_key][:,dim_to_be_sliced] : ",bin_points_dict[bin_key][:,dim_to_be_sliced])
				median = np.median(bin_points_dict[bin_key][:,dim_to_be_sliced])
				#print("median : ",median)

				a_bin_boundary, b_bin_boundary = bin_boundary.copy(), bin_boundary.copy()
				#print("a_bin_boundary : ",a_bin_boundary)
				a_bin_boundary[dim_to_be_sliced,1] = median
				b_bin_boundary[dim_to_be_sliced,0] = median
				bin_boundaries_dict[str(split_number)+bin_key[1:]+'a'] = a_bin_boundary
				bin_boundaries_dict[str(split_number)+bin_key[1:]+'b'] = b_bin_boundary
				
				a_points, b_points = [],[] 
				for event_number in range(bin_points_dict[bin_key].shape[0]):
					if bin_points_dict[bin_key][event_number,dim_to_be_sliced] < median: a_points.append(bin_points_dict[bin_key][event_number,:].tolist())
					else: b_points.append(bin_points_dict[bin_key][event_number,:].tolist())
				
				bin_points_dict[str(split_number)+bin_key[1:]+'a'] = np.array(a_points)
				bin_points_dict[str(split_number)+bin_key[1:]+'b'] = np.array(b_points)
				#If a bin contains no particles it should be deleted
				if len(a_points)==0: 
					del bin_points_dict[str(split_number)+bin_key[1:]+'a']
					del bin_boundaries_dict[str(split_number)+bin_key[1:]+'a']
				if len(b_points)==0: 
					del bin_points_dict[str(split_number)+bin_key[1:]+'b']
					del bin_boundaries_dict[str(split_number)+bin_key[1:]+'b']

	if PLOT:	pickle.dump( bin_boundaries_dict, open( "bin_boundaries_dict.p", "wb" ) )	
	bins_sample01_dict= {}
	signed_Scp2_dict= {}
	results_list = []

	for number_of_splits in number_of_splits_list:
		print("\nnumber_of_splits : ",number_of_splits,"\nsystematics_fraction : ",systematics_fraction)
		bins_sample0, bins_sample1 = [] , [] 
		for bin_key, bin_points in bin_points_dict.items():
			if str(number_of_splits) in bin_key:
				labels_in_bin = bin_points[:,-1]
				#print("labels_in_bin : ",labels_in_bin)
				bin_sample0 = np.count_nonzero( labels_in_bin == 0)
				bin_sample1 = np.count_nonzero( labels_in_bin == 1)
				#print("bin_sample0 : ",bin_sample0) 
				#print("bin_sample1 : ",bin_sample1) 
				#simulate uncertainties 
				if(systematics_fraction*float(bin_sample0)!=0.): bin_sample0 += int(round(np.random.normal(0.,systematics_fraction*float(bin_sample0))))
				if(systematics_fraction*float(bin_sample1)!=0.): bin_sample1 += int(round(np.random.normal(0.,systematics_fraction*float(bin_sample1))))
				bins_sample01_dict[bin_key]=[bin_sample0,bin_sample1]
				signed_Scp2_dict[bin_key] = np.square(float(bin_sample1-bin_sample0))/(float(bin_sample1)+float(bin_sample0)+np.square(float(bin_sample1)*systematics_fraction)+np.square(float(bin_sample1)*systematics_fraction))*np.sign(bin_sample1-bin_sample0)
				#print("\n\nbin_sample0 : ",bin_sample0, "\n bins_sample0 : ", bins_sample0 )
				#print("type(bin_sample0) : ",type(bin_sample0), "  type(bins_sample0) : ",type(bins_sample0))
				bins_sample0.append(bin_sample0)
				#print("  bins_sample0 : ", bins_sample0, "\n\n" )
				bins_sample1.append(bin_sample1)
		bins_sample0, bins_sample1 = np.array(bins_sample0,dtype=float), np.array(bins_sample1, dtype=float)
		print("bins_sample0 : ",bins_sample0,"\n bins_sample1 : ",bins_sample1)
		#element wise subtraction and division
		Scp2 = ((bins_sample1-bins_sample0)**2)/ (bins_sample1+bins_sample0+(systematics_fraction*bins_sample1)**2+(systematics_fraction*bins_sample0)**2 )
		#Scp2 =  np.divide(np.square(np.subtract(bins_sample1,bins_sample0)),np.add(bins_sample1,bins_sample0))
		if DEBUG:
			print(Scp2)

		#nansum ignores all the contributions that are Not A Number (NAN)
		Chi2 = np.nansum(Scp2)
		if DEBUG:
			print("Chi2")
			print(Chi2)
		dof=bins_sample0.shape[0]-1
	
		pvalue= 1 - stats.chi2.cdf(Chi2,dof)
		
		print("\nThe p value for Scp2 = ",Scp2," and Chi2 = ", Chi2, " is ",pvalue,"\n\n")
		if DEBUG:
			print(bins_sample0)
			print(bins_sample1)
			print("Chi2/dof : {0}".format(str(Chi2/dof)))

			print("pvalue : {0}".format(str(pvalue)))

		results_list.append(pvalue)

		if PLOT: 
			if no_dim==1: chi2_plots.adaptive_binning_1Dplot(bin_boundaries_dict,data,number_of_splits,title+" "+str(no_dim) + "D "+str(number_of_splits)+ " splits ",name+"_"+str(no_dim) + "D_chi2_"+str(number_of_splits)+"_splits")
			if no_dim==2: chi2_plots.adaptive_binning_2Dplot(bin_boundaries_dict,signed_Scp2_dict,number_of_splits,X_values,title+" "+str(no_dim) + "D"+str(number_of_splits)+ " splits ",name+"_"+str(no_dim) + "D_chi2_"+str(number_of_splits)+"_splits", X_min= X_min,X_total_width=X_total_width )
                        if no_dim>1:  chi2_plots.adaptive_binning_2D1Dplot(bin_boundaries_dict,bins_sample01_dict,number_of_splits,X_values,title+" "+str(no_dim) + "D"+str(number_of_splits)+ " splits ",name+"_"+str(no_dim) + "D_chi2_"+str(number_of_splits)+"_splits", no_dim)



	return results_list
356250
def plot_summary(self, cycles=None):
        """"""
        if cycles is None:
            cycles = [0]
        fig1 = plt.figure()
        ax1 = fig1.add_subplot(221)
        ax1.set_title('Fit')
        ax2 = fig1.add_subplot(222)
        ax2.set_title('OCV')
        ax3 = fig1.add_subplot(223)
        ax3.set_title('Tau')
        ax3.set_yscale("log")
        ax4 = fig1.add_subplot(224)
        ax4.set_title('Voltage Impact')

        plot_data = self.get_best_fit_data()
        for cycle in cycles:
            ax1.plot(plot_data[cycle][0], plot_data[cycle][1])
            ax1.plot(plot_data[cycle][0], plot_data[cycle][2])

        plot_data = self.get_best_fit_parameters_grouped()

        for i in range(self.circuits):
            ax3.plot(self.get_fit_cycles(), plot_data['t' + str(i)])
            ax4.plot(self.get_fit_cycles(), plot_data['w' + str(i)])

        ax2.plot(self.get_fit_cycles(), plot_data['ocv'])
263564
def sine_wave(frequency):
  """"""
  xs = tf.reshape(tf.range(_samples(), dtype=tf.float32), [1, _samples(), 1])
  ts = xs / FLAGS.sample_rate
  return tf.sin(2 * math.pi * frequency * ts)
414208
def read_freq_dielectric(self):
        

        plasma_pattern = r"plasma frequency squared.*"
        dielectric_pattern = r"frequency dependent\s+IMAGINARY " \
                             r"DIELECTRIC FUNCTION \(independent particle, " \
                             r"no local field effects\)(\sdensity-density)*$"
        row_pattern = r"\s+".join([r"([\.\-\d]+)"] * 3)
        plasma_frequencies = collections.defaultdict(list)
        read_plasma = False
        read_dielectric = False
        energies = []
        data = {"REAL": [], "IMAGINARY": []}
        count = 0
        component = "IMAGINARY"
        with zopen(self.filename, "rt") as f:
            for l in f:
                l = l.strip()
                if re.match(plasma_pattern, l):
                    read_plasma = "intraband" if "intraband" in l else "interband"
                elif re.match(dielectric_pattern, l):
                    read_plasma = False
                    read_dielectric = True
                    row_pattern = r"\s+".join([r"([\.\-\d]+)"] * 7)

                if read_plasma and re.match(row_pattern, l):
                    plasma_frequencies[read_plasma].append(
                        [float(t) for t in l.strip().split()])
                elif read_dielectric:
                    if re.match(row_pattern, l.strip()):
                        toks = l.strip().split()
                        if component == "IMAGINARY":
                            energies.append(float(toks[0]))
                        xx, yy, zz, xy, yz, xz = [float(t) for t in toks[1:]]
                        matrix = [[xx, xy, xz], [xy, yy, yz], [xz, yz, zz]]
                        data[component].append(matrix)
                    elif re.match(r"\s*-+\s*", l):
                        count += 1
                    if count == 2:
                        component = "REAL"
                    elif count == 3:
                        break

        self.plasma_frequencies = {k: np.array(v[:3])
                                   for k, v in plasma_frequencies.items()}
        self.dielectric_energies = np.array(energies)
        self.dielectric_tensor_function = np.array(data["REAL"]) + \
                                          1j * np.array(data["IMAGINARY"])
1002326
def plot_xbt(fignum, XB, T, e, b):
    
    plt.figure(num=fignum)
    plt.clf()
    if not isServer:
        plt.figtext(.02, .01, version_num)
    plt.xlabel('Field (T)')
    plt.ylabel('Susceptibility (m^3/kg)')
    k = 0
    B, X = [], []
    for xb in XB:
        X.append(xb[0])
        B.append(xb[1])
    plt.plot(B, X)
    plt.legend(['%i' % (int(T)) + ' K'])
    plt.title(e + ': f = ' + '%i' % (int(f)) + ' Hz')
381239
def mavfft_display(logfile, condition=None):
    ''''''

    '''object to store data about a single FFT plot'''
    class PlotData(object):
        def __init__(self, ffth):
            self.seqno = -1
            self.fftnum = ffth.N
            self.sensor_type = ffth.type
            self.instance = ffth.instance
            self.sample_rate_hz = ffth.smp_rate
            self.multiplier = ffth.mul
            self.data = {}
            self.data["X"] = []
            self.data["Y"] = []
            self.data["Z"] = []
            self.holes = False
            self.freq = None

        def add_fftd(self, fftd):
            if fftd.N != self.fftnum:
                print("Skipping ISBD with wrong fftnum (%u vs %u)\n" % (fftd.fftnum, self.fftnum))
                return
            if self.holes:
                print("Skipping ISBD(%u) for ISBH(%u) with holes in it" % (fftd.seqno, self.fftnum))
                return
            if fftd.seqno != self.seqno+1:
                print("ISBH(%u) has holes in it" % fftd.N)
                self.holes = True
                return
            self.seqno += 1
            self.data["X"].extend(fftd.x)
            self.data["Y"].extend(fftd.y)
            self.data["Z"].extend(fftd.z)

        def prefix(self):
            if self.sensor_type == 0:
                return "Accel"
            elif self.sensor_type == 1:
                return "Gyro"
            else:
                return "?Unknown Sensor Type?"

        def tag(self):
            return str(self)

        def __str__(self):
            return "%s[%u]" % (self.prefix(), self.instance)

    print("Processing log for ISBH and ISBD messages")

    things_to_plot = []
    plotdata = None
    start_time = time.time()
    mlog = mavutil.mavlink_connection(logfile)
    while True:
        m = mlog.recv_match(type=['ISBH','ISBD'],condition=condition)
        if m is None:
            break
        msg_type = m.get_type()
        if msg_type == "ISBH":
            if plotdata is not None:
                # close off previous data collection
                things_to_plot.append(plotdata)
            # initialise plot-data collection object
            plotdata = PlotData(m)
            continue

        if msg_type == "ISBD":
            if plotdata is None:
                continue
            plotdata.add_fftd(m)

    if len(things_to_plot) == 0:
        print("No FFT data. Did you set INS_LOG_BAT_MASK?")
        return
    time_delta = time.time() - start_time
    print("Extracted %u fft data sets" % len(things_to_plot))

    sum_fft = {}
    freqmap = {}
    count = 0

    first_freq = None
    for thing_to_plot in things_to_plot:
        for axis in [ "X","Y","Z" ]:
            d = numpy.array(thing_to_plot.data[axis])/float(thing_to_plot.multiplier)
            if len(d) == 0:
                print("No data?!?!?!")
                continue
            
            avg = numpy.sum(d) / len(d)
            d -= avg
            d_fft = numpy.fft.rfft(d)
            if thing_to_plot.tag() not in sum_fft:
                sum_fft[thing_to_plot.tag()] = {
                    "X": 0,
                    "Y": 0,
                    "Z": 0
                }
            sum_fft[thing_to_plot.tag()][axis] = numpy.add(sum_fft[thing_to_plot.tag()][axis], d_fft)
            count += 1
            freq = numpy.fft.rfftfreq(len(d), 1.0/thing_to_plot.sample_rate_hz)
            freqmap[thing_to_plot.tag()] = freq

    for sensor in sum_fft:
        pylab.figure(str(sensor))
        for axis in [ "X","Y","Z" ]:
            pylab.plot(freqmap[sensor], numpy.abs(sum_fft[sensor][axis]/count), label=axis)
        pylab.legend(loc='upper right')
        pylab.xlabel('Hz')

    pylab.show()
1342055
def plot_probe(mea, ax=None, xlim=None, ylim=None, color_currents=False, top=None, bottom=None,
               cmap='viridis', type='shank', alpha_elec=0.7, alpha_prb=0.3):
    '''
    
    

    '''


    if ax is None:
        fig = plt.figure()
        ax = fig.add_subplot(111)

    n_elec = mea.positions.shape[0]
    elec_size = mea.size

    if mea.type == 'mea':
        mea_pos = np.array([np.dot(mea.positions, mea.main_axes[0]), np.dot(mea.positions, mea.main_axes[1])]).T

        min_x, max_x = [np.min(np.dot(mea.positions, mea.main_axes[0])),
                        np.max(np.dot(mea.positions, mea.main_axes[0]))]
        center_x = (min_x + max_x)/2.
        min_y, max_y = [np.min(np.dot(mea.positions, mea.main_axes[1])),
                        np.max(np.dot(mea.positions, mea.main_axes[1]))]
        center_y = (min_y + max_y)/2.

        if type == 'shank':
            if top is None:
                probe_height = 200
            else:
                probe_height = top
            probe_top = max_y + probe_height
            if bottom is None:
                probe_bottom = min_y - probe_height
            else:
                probe_bottom = min_y - bottom
            probe_corner = min_y - 0.1*probe_height
            probe_left = min_x - 0.1*probe_height
            probe_right = max_x + 0.1*probe_height

            verts = [
                (min_x - 2*elec_size, probe_top),  # left, bottom
                (min_x - 2*elec_size, probe_corner),  # left, top
                (center_x, probe_bottom),  # right, top
                (max_x + 2*elec_size, probe_corner),  # right, bottom
                (max_x + 2*elec_size, probe_top),
                (min_x - 2 * elec_size, max_y + 2 * elec_size) # ignored
            ]

            codes = [Path.MOVETO,
                     Path.LINETO,
                     Path.LINETO,
                     Path.LINETO,
                     Path.LINETO,
                     Path.CLOSEPOLY,
                     ]

        elif type == 'planar':
            probe_top = max_y + 2 * elec_size
            probe_bottom = min_y - 2 * elec_size
            probe_left = min_x - 2 * elec_size
            probe_right = max_x + 2 * elec_size

            verts = [
                (min_x - 2 * elec_size, max_y + 2 * elec_size),  # left, bottom
                (min_x - 2 * elec_size, min_y - 2 * elec_size),  # left, top
                (max_x + 2 * elec_size, min_y - 2 * elec_size),  # right, bottom
                (max_x + 2 * elec_size, max_y + 2 * elec_size),  # ignored
                (max_x + 2 * elec_size, max_y + 2 * elec_size)   # ignored
            ]

            codes = [Path.MOVETO,
                     Path.LINETO,
                     Path.LINETO,
                     Path.LINETO,
                     Path.CLOSEPOLY,
                     ]
        else:
            raise AttributeError("'type' can be 'shank' or 'planar'")

        path = Path(verts, codes)

        patch = patches.PathPatch(path, facecolor='green', edgecolor='k', lw=0.5, alpha=alpha_prb)
        ax.add_patch(patch)

        if color_currents:
            norm_curr = mea.currents / np.max(np.abs(mea.currents))
            colormap = plt.get_cmap(cmap)
            elec_colors = colormap(norm_curr)
        else:
            elec_colors = ['orange'] * mea.number_electrodes

        if mea.shape == 'square':
            for e in range(n_elec):
                elec = patches.Rectangle((mea_pos[e, 0] - elec_size, mea_pos[e, 1] - elec_size), 2*elec_size,
                                         2*elec_size, alpha=alpha_elec, facecolor=elec_colors[e],
                                         edgecolor=[0.3, 0.3, 0.3], lw=0.5)

                ax.add_patch(elec)
        elif mea.shape == 'circle':
            for e in range(n_elec):
                elec = patches.Circle((mea_pos[e, 0], mea_pos[e, 1]), elec_size, alpha=alpha_elec,
                                      facecolor=elec_colors[e], edgecolor=[0.3, 0.3, 0.3], lw=0.5)

                ax.add_patch(elec)
    else:
        raise NotImplementedError('Wire type plotting not implemented')

    ax.axis('equal')

    if xlim:
        ax.set_xlim(xlim)
    else:
        ax.set_xlim(probe_left - 5 * elec_size, probe_right + 5 * elec_size)
    if ylim:
        ax.set_ylim(ylim)
    else:
        ax.set_ylim(probe_bottom - 5 * elec_size, probe_top + 5 * elec_size)

    return ax
256484
def freqpoly_plot(data):
        """"""
        rel_data = OrderedDict()
        for key, val in data.items():
            tot = sum(val.values(), 0)
            rel_data[key] = {k: v / tot for k, v in val.items()}
        fplotconfig = {
            'data_labels': [
                {'name': 'Absolute', 'ylab': 'Frequency', 'xlab': 'Merged Read Length'},
                {'name': 'Relative', 'ylab': 'Relative Frequency', 'xlab': 'Merged Read Length'}
                ],
            'id': 'flash_freqpoly_plot', 'title': 'FLASh: Frequency of merged read lengths',
            'colors': dict(zip(data.keys(), MultiqcModule.get_colors(len(data))))
            }
        return linegraph.plot([data, rel_data], fplotconfig)
731045
def diagnostics_plot_chisq(self, ds, figname = "modelfit_chisqs.png"):
        
        label_names = ds.get_plotting_labels()
        lams = ds.wl
        pivots = self.pivots
        npixels = len(lams)
        nlabels = len(pivots)
        chisqs = self.chisqs
        coeffs = self.coeffs
        scatters = self.scatters

        # Histogram of the chi squareds of ind. stars
        plt.hist(np.sum(chisqs, axis=0), color='lightblue', alpha=0.7,
                bins=int(np.sqrt(len(chisqs))))
        dof = len(lams) - coeffs.shape[1]   # for one star
        plt.axvline(x=dof, c='k', linewidth=2, label="DOF")
        plt.legend()
        plt.title("Distribution of " + r"$\chi^2$" + " of the Model Fit")
        plt.ylabel("Count")
        plt.xlabel(r"$\chi^2$" + " of Individual Star")
        print("Diagnostic plot: histogram of the red chi squareds of the fit")
        print("Saved as %s" %figname)
        plt.savefig(figname)
        plt.close()