def plot_position(self, pos_true, pos_est):
N = pos_est.shape[1]
pos_true = pos_true[:, :N]
pos_est = pos_est[:, :N]
# Figure
plt.figure()
plt.suptitle("Position")
# Ground truth
plt.plot(pos_true[0, :], pos_true[1, :],
color="red", marker="o", label="Grouth truth")
# Estimated
plt.plot(pos_est[0, :], pos_est[1, :],
color="blue", marker="o", label="Estimated")
# Plot labels and legends
plt.xlabel("East (m)")
plt.ylabel("North (m)")
plt.axis("equal")
plt.legend(loc=0)
python类axis()的实例源码
def show_samples(y, ndim, nb=10, cmap=''):
if ndim == 4:
for i in range(nb**2):
plt.subplot(nb, nb, i+1)
plt.imshow(y[i], cmap=cmap, interpolation='none')
plt.axis('off')
else:
x = y[0]
y = y[1]
plt.figure(0)
for i in range(10):
plt.subplot(2, 5, i+1)
plt.imshow(x[i], cmap=cmap, interpolation='none')
plt.axis('off')
plt.figure(1)
for i in range(10):
plt.subplot(2, 5, i+1)
plt.imshow(y[i], cmap=cmap, interpolation='none')
plt.axis('off')
plt.show()
def linear_testing():
x_axis = np.linspace(1, 51, 100)
x_nice = np.linspace(x_axis[0], x_axis[-1], 100)
mod, params = qudi_fitting.make_linear_model()
print('Parameters of the model', mod.param_names, ' with the independet variable', mod.independent_vars)
params['slope'].value = 2 # + abs(np.random.normal(0,1))
params['offset'].value = 50 #+ abs(np.random.normal(0, 200))
#print('\n', 'beta', params['beta'].value, '\n', 'lifetime',
#params['lifetime'].value)
data_noisy = (mod.eval(x=x_axis, params=params)
+ 10 * np.random.normal(size=x_axis.shape))
result = qudi_fitting.make_linear_fit(axis=x_axis, data=data_noisy, add_parameters=None)
plt.plot(x_axis, data_noisy, 'ob')
plt.plot(x_nice, mod.eval(x=x_nice, params=params), '-g')
print(result.fit_report())
plt.plot(x_axis, result.best_fit, '-r', linewidth=2.0)
plt.plot(x_axis, result.init_fit, '-y', linewidth=2.0)
plt.show()
def fit_data():
data=np.loadtxt('data.dat')
print(data)
params = dict()
params["c"] = {"min" : -np.inf,"max" : np.inf}
result = qudi_fitting.make_lorentzian_fit(axis=data[:,0], data=data[:,3], add_parameters=params)
print(result.fit_report())
plt.plot(data[:,0],-data[:,3]+2,"b-o",label="data mean")
# plt.plot(data[:,0],data[:,1],label="data")
# plt.plot(data[:,0],data[:,2],label="data")
plt.plot(data[:,0],-result.best_fit+2,"r-",linewidth=2.,label="fit")
# plt.plot(data[:,0],result.init_fit,label="init")
plt.xlabel("time (ns)")
plt.ylabel("polarization transfer (arb. u.)")
plt.legend(loc=1)
# plt.savefig("pol20_24repetition_pol.pdf")
# plt.savefig("pol20_24repetition_pol.png")
plt.show()
savedata=[[data[ii,0],-data[ii,3]+2,-result.best_fit[ii]+2] for ii in range(len(data[:,0]))]
np.savetxt("pol_data_fit.csv",savedata)
# print(result.params)
print(result.params)
def plotSequenceForRotatingState(degPerStep, Sigma, T=1000):
theta = degPerStep * np.pi / 180. # radPerStep
A = np.asarray([
[np.cos(theta), -np.sin(theta)],
[np.sin(theta), np.cos(theta)],
])
Sigma = np.asarray(Sigma)
if Sigma.size == 1:
Sigma = Sigma * np.eye(2)
elif Sigma.size == 2:
Sigma = np.diag(Sigma)
X = np.zeros((T, 2))
X[0, :] = [1, 0]
for t in xrange(1, T):
X[t] = np.random.multivariate_normal(np.dot(A, X[t - 1]), Sigma)
pylab.plot(X[:, 0], X[:, 1], '.')
pylab.axis([-4, 4, -4, 4])
pylab.axis('equal')
def showExampleDocs(pylab=None, nrows=3, ncols=3):
if pylab is None:
from matplotlib import pylab
Data = get_data(seed=0, nObsPerDoc=200)
PRNG = np.random.RandomState(0)
chosenDocs = PRNG.choice(Data.nDoc, nrows * ncols, replace=False)
for ii, d in enumerate(chosenDocs):
start = Data.doc_range[d]
stop = Data.doc_range[d + 1]
Xd = Data.X[start:stop]
pylab.subplot(nrows, ncols, ii + 1)
pylab.plot(Xd[:, 0], Xd[:, 1], 'k.')
pylab.axis('image')
pylab.xlim([-1.5, 1.5])
pylab.ylim([-1.5, 1.5])
pylab.xticks([])
pylab.yticks([])
pylab.tight_layout()
# Set Toy Parameters
###########################################################
def plot1D_mat(a, b, M, title=''):
""" Plot matrix M with the source and target 1D distribution
Creates a subplot with the source distribution a on the left and
target distribution b on the tot. The matrix M is shown in between.
Parameters
----------
a : np.array, shape (na,)
Source distribution
b : np.array, shape (nb,)
Target distribution
M : np.array, shape (na,nb)
Matrix to plot
"""
na, nb = M.shape
gs = gridspec.GridSpec(3, 3)
xa = np.arange(na)
xb = np.arange(nb)
ax1 = pl.subplot(gs[0, 1:])
pl.plot(xb, b, 'r', label='Target distribution')
pl.yticks(())
pl.title(title)
ax2 = pl.subplot(gs[1:, 0])
pl.plot(a, xa, 'b', label='Source distribution')
pl.gca().invert_xaxis()
pl.gca().invert_yaxis()
pl.xticks(())
pl.subplot(gs[1:, 1:], sharex=ax1, sharey=ax2)
pl.imshow(M, interpolation='nearest')
pl.axis('off')
pl.xlim((0, nb))
pl.tight_layout()
pl.subplots_adjust(wspace=0., hspace=0.2)
def plot_position(self, pos_true, pos_est, cam_states):
N = pos_est.shape[1]
pos_true = pos_true[:, :N]
pos_est = pos_est[:, :N]
# Figure
plt.figure()
plt.suptitle("Position")
# Ground truth
plt.plot(pos_true[0, :], pos_true[1, :],
color="red", label="Grouth truth")
# color="red", marker="x", label="Grouth truth")
# Estimated
plt.plot(pos_est[0, :], pos_est[1, :],
color="blue", label="Estimated")
# color="blue", marker="o", label="Estimated")
# Sliding window
cam_pos = []
for cam_state in cam_states:
cam_pos.append(cam_state.p_G)
cam_pos = np.array(cam_pos).reshape((len(cam_pos), 3)).T
plt.plot(cam_pos[0, :], cam_pos[1, :],
color="green", label="Camera Poses")
# color="green", marker="o", label="Camera Poses")
# Plot labels and legends
plt.xlabel("East (m)")
plt.ylabel("North (m)")
plt.axis("equal")
plt.legend(loc=0)
def plot_velocity(self, timestamps, vel_true, vel_est):
N = vel_est.shape[1]
t = timestamps[:N]
vel_true = vel_true[:, :N]
vel_est = vel_est[:, :N]
# Figure
plt.figure()
plt.suptitle("Velocity")
# X axis
plt.subplot(311)
plt.plot(t, vel_true[0, :], color="red", label="Ground_truth")
plt.plot(t, vel_est[0, :], color="blue", label="Estimate")
plt.title("x-axis")
plt.xlabel("Date Time")
plt.ylabel("ms^-1")
plt.legend(loc=0)
# Y axis
plt.subplot(312)
plt.plot(t, vel_true[1, :], color="red", label="Ground_truth")
plt.plot(t, vel_est[1, :], color="blue", label="Estimate")
plt.title("y-axis")
plt.xlabel("Date Time")
plt.ylabel("ms^-1")
plt.legend(loc=0)
# Z axis
plt.subplot(313)
plt.plot(t, vel_true[2, :], color="red", label="Ground_truth")
plt.plot(t, vel_est[2, :], color="blue", label="Estimate")
plt.title("z-axis")
plt.xlabel("Date Time")
plt.ylabel("ms^-1")
plt.legend(loc=0)
def plot_attitude(self, timestamps, att_true, att_est):
# Setup
N = att_est.shape[1]
t = timestamps[:N]
att_true = att_true[:, :N]
att_est = att_est[:, :N]
# Figure
plt.figure()
plt.suptitle("Attitude")
# X axis
plt.subplot(311)
plt.plot(t, att_true[0, :], color="red", label="Ground_truth")
plt.plot(t, att_est[0, :], color="blue", label="Estimate")
plt.title("x-axis")
plt.legend(loc=0)
plt.xlabel("Date Time")
plt.ylabel("rad s^-1")
# Y axis
plt.subplot(312)
plt.plot(t, att_true[1, :], color="red", label="Ground_truth")
plt.plot(t, att_est[1, :], color="blue", label="Estimate")
plt.title("y-axis")
plt.legend(loc=0)
plt.xlabel("Date Time")
plt.ylabel("rad s^-1")
# Z axis
plt.subplot(313)
plt.plot(t, att_true[2, :], color="red", label="Ground_truth")
plt.plot(t, att_est[2, :], color="blue", label="Estimate")
plt.title("z-axis")
plt.legend(loc=0)
plt.xlabel("Date Time")
plt.ylabel("rad s^-1")
def plot_velocity(self, timestamps, vel_true, vel_est):
N = vel_est.shape[1]
t = timestamps[:N]
vel_true = vel_true[:, :N]
vel_est = vel_est[:, :N]
# Figure
plt.figure()
plt.suptitle("Velocity")
# X axis
plt.subplot(311)
plt.plot(t, vel_true[0, :], color="red", label="Ground_truth")
plt.plot(t, vel_est[0, :], color="blue", label="Estimate")
plt.title("x-axis")
plt.xlabel("Date Time")
plt.ylabel("ms^-1")
plt.legend(loc=0)
# Y axis
plt.subplot(312)
plt.plot(t, vel_true[1, :], color="red", label="Ground_truth")
plt.plot(t, vel_est[1, :], color="blue", label="Estimate")
plt.title("y-axis")
plt.xlabel("Date Time")
plt.ylabel("ms^-1")
plt.legend(loc=0)
# Z axis
plt.subplot(313)
plt.plot(t, vel_true[2, :], color="red", label="Ground_truth")
plt.plot(t, vel_est[2, :], color="blue", label="Estimate")
plt.title("z-axis")
plt.xlabel("Date Time")
plt.ylabel("ms^-1")
plt.legend(loc=0)
def axis(self, imin, imax, jmin, jmax, xlabels=[], ylabels=[]):
imid=(imin+imax)/2 # midpoint of X-axis
jmid=(jmin+jmax)/2 # midpoint of Y-axis
# Create axis lines
self.create_line((imin, jmax, imax, jmax))
self.create_line((imin, jmin, imin, jmax))
self.create_line((imin, jmin, imax, jmin))
self.create_line((imax, jmin, imax, jmax))
self.create_line((imid, jmin, imid, jmax))
self.create_line((imin, jmid, imax, jmid))
# Create tick marks and labels
tic = imin
for label in xlabels:
self.create_line((tic, jmax+ 5, tic, jmax))
self.create_text( tic, jmax+10, text=label)
if len(xlabels)!=1:
tic+=(imax-imin)/(len(xlabels)-1)
tic = jmax
for label in ylabels:
self.create_line((imin , tic, imin-5, tic))
self.create_text( imin-20, tic, text=label)
if len(ylabels)!=1:
tic-=(jmax-jmin)/(len(ylabels)-1)
def fancy_show(y, cmap=''):
x = y[0]
y = y[1]
plt.figure(0)
for i in range(100):
plt.subplot(10, 10, i+1)
plt.imshow(x[i], cmap=cmap, interpolation='none')
plt.axis('off')
plt.figure(1)
for i in range(100):
plt.subplot(10, 10, i+1)
plt.imshow(y[i], cmap=cmap, interpolation='none')
plt.axis('off')
plt.show()
def imshow_(x, **kwargs):
if x.ndim == 2:
plt.imshow(x, interpolation="nearest", **kwargs)
elif x.ndim == 1:
plt.imshow(x[:,None].T, interpolation="nearest", **kwargs)
plt.yticks([])
plt.axis("tight")
# ------------- Data -------------
def draw(name, feature4096):
plt.figure(name)
feature4096 = feature4096.reshape(64,64)
for i, x in enumerate(feature4096):
for j, y in enumerate(x):
#if y <= 0:
# print ' ',
#else:
# print '%3.1f'%y,
plt.scatter([j],[i], s=[y*1000])
#print
plt.axis([-1, 65, -1, 65])
plt.show()
def plot_validation_cost(train_error, val_error, class_rate=None, savefilename=None):
epochs = range(len(train_error))
fig, ax1 = plt.subplots()
ax1.plot(epochs, train_error, label='train error')
ax1.plot(epochs, val_error, label='validation error')
ax1.set_xlabel('epoch')
ax1.set_ylabel('cost')
plt.title('Validation Cost')
lines = ax1.get_lines()
# Shrink current axis's height by 10% on the bottom
box = ax1.get_position()
ax1.set_position([box.x0, box.y0 + box.height * 0.1,
box.width, box.height * 0.9])
if class_rate is not None:
ax2 = plt.twinx(ax1)
ax2.plot(epochs, class_rate, label='classification rate', color='r')
ax2.set_ylabel('classification rate')
lines.extend(ax2.get_lines())
ax2.set_position([box.x0, box.y0 + box.height * 0.1,
box.width, box.height * 0.9])
labels = [l.get_label() for l in lines]
# Put a legend below current axis
ax1.legend(lines, labels, loc='upper center', bbox_to_anchor=(0.5, -0.05),
fancybox=False, shadow=False, ncol=5)
# ax1.legend(lines, labels, loc='lower right')
if savefilename:
plt.savefig(savefilename)
plt.show()
def visualize_sequence(sequence, shape=(30, 40), savefilename=None, title='sequence'):
cols = int(math.ceil(len(sequence) / 2.0))
vis = tile_raster_images(sequence, shape, (2, cols), tile_spacing=(1, 1))
plt.title(title)
plt.axis('off')
plt.show()
if savefilename:
o = Image.fromarray(vis)
o.save('{}.png'.format(savefilename))
def plotSequenceForRotatingState3D(degPerStep, Sigma, stationaryDim=0, T=1000):
A = makeA_3DRotationMatrix(degPerStep, stationaryDim)
Sigma = Sigma * np.eye(3)
assert Sigma.shape == (3, 3)
X = np.zeros((T, 3))
X[0, :] = [1, 1, 1]
for t in xrange(1, T):
X[t] = np.random.multivariate_normal(np.dot(A, X[t - 1]), Sigma)
ax = Axes3D(pylab.figure())
pylab.plot(X[:, 0], X[:, 1], X[:, 2], '.')
pylab.axis('equal')
def ricker_matrix(width, resolution, n_components):
"""Dictionary of Ricker (Mexican hat) wavelets"""
centers = np.linspace(0, resolution - 1, n_components)
D = np.empty((n_components, resolution))
for i, center in enumerate(centers):
D[i] = ricker_function(resolution, center, width)
D /= np.sqrt(np.sum(D ** 2, axis=1))[:, np.newaxis]
return D
def get_2DsersicIeff(value,reff,sersicindex,axisratio,boxiness=0.0,returnFtot=False):
"""
Get the surface brightness value at the effective radius of a 2D sersic profile (given GALFIT Sersic parameters).
Ieff is calculated using ewuations (4) and (5) in Peng et al. (2010), AJ 139:2097.
This Ieff is what is referred to as 'amplitude' in astropy.modeling.models.Sersic2D
used in tdose_utilities.gen_2Dsersic()
--- INPUT ---
value If returnFtot=False "value" corresponds to Ftot of the profile (total flux for profile integrated
til r=infty) and Ieff will be returned.
If instead returnFtot=True "value" should provide Ieff so Ftot can be returned
reff Effective radius
sersicindex Sersic index of profile
axisratio Ratio between the minor and major axis (0<axisratio<1)
boxiness The boxiness of the profile
returnFtot If Ftot is not known, but Ieff is, set returnFtot=True to return Ftot instead (providing Ieff to "value")
--- EXAMPLE OF USE ---
Ieff = 1.0
reff = 25.0
sersicindex = 4.0
axisratio = 1.0
Ftot_calc = tu.get_2DsersicIeff(Ieff,reff,sersicindex,axisratio,returnFtot=True)
Ieff_calc = tu.get_2DsersicIeff(Ftot_calc,reff,sersicindex,axisratio)
size = 1000
x,y = np.meshgrid(np.arange(size), np.arange(size))
mod = Sersic2D(amplitude = Ieff, r_eff = reff, n=sersicindex, x_0=size/2.0, y_0=size/2.0, ellip=1-axisratio, theta=-1)
img = mod(x, y)
hducube = pyfits.PrimaryHDU(img)
hdus = [hducube]
hdulist = pyfits.HDUList(hdus)
hdulist.writeto('/Volumes/DATABCKUP2/TDOSEextractions/models_cutouts/model_sersic_spherical.fits',clobber=True)
"""
gam2n = scipy.special.gamma(2.0*sersicindex)
kappa = scipy.special.gammaincinv(2.0*sersicindex,0.5)
Rfct = np.pi * (boxiness + 2.) / (4. * scipy.special.beta(1./(boxiness+2.),1.+1./(boxiness+2.)) )
factor = 2.0 * np.pi * reff**2.0 * np.exp(kappa) * sersicindex * kappa**(-2*sersicindex) * gam2n * axisratio / Rfct
if returnFtot:
Ftot = value * factor
return Ftot
else:
Ieff = value / factor
return Ieff
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def roll_2Dprofile(profile,position,padvalue=0.0,showprofiles=False):
"""
Move 2D profile to given position in array by rolling it in x and y.
Note that the roll does not handle sub-pixel moves.
tu.shift_2Dprofile() does this using interpolation
--- INPUT ---
profile profile to shift
position position to move center of image (profile) to: [ypos,xpos]
padvalue the values to padd the images with when shifting profile
showprofiles Show profile when shifted?
--- EXAMPLE OF USE ---
tu.roll_2Dprofile(gauss2D,)
"""
profile_dim = profile.shape
yroll = np.int(np.round(position[0]-profile_dim[0]/2.))
xroll = np.int(np.round(position[1]-profile_dim[1]/2.))
profile_shifted = np.roll(np.roll(profile,yroll,axis=0),xroll,axis=1)
if showprofiles:
vmaxval = np.max(profile_shifted)
plt.imshow(profile_shifted,interpolation='none',vmin=-vmaxval, vmax=vmaxval)
plt.title('Positioned Source')
plt.show()
if yroll < 0:
profile_shifted[yroll:,:] = padvalue
else:
profile_shifted[:yroll,:] = padvalue
if xroll < 0:
profile_shifted[:,xroll:] = padvalue
else:
profile_shifted[:,:xroll] = padvalue
if showprofiles:
plt.imshow(profile_shifted,interpolation='none',vmin=-vmaxval, vmax=vmaxval)
plt.title('Positioned Source with 0s inserted')
plt.show()
return profile_shifted
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def save_image(nifti, anat, cluster_dict, out_path, f, image_threshold=2,
texcol=1, bgcol=0, iscale=2, text=None, **kwargs):
'''Saves a single nifti image.
Args:
nifti (str or nipy.core.api.image.image.Image): nifti file to visualize.
anat (nipy.core.api.image.image.Image): anatomical nifti file.
cluster_dict (dict): dictionary of clusters.
f (int): index.
image_threshold (float): treshold for `plot_map`.
texcol (float): text color.
bgcol (float): background color.
iscale (float): image scale.
text (Optional[str]): text for figure.
**kwargs: extra keyword arguments
'''
if isinstance(nifti, str):
nifti = load_image(nifti)
feature = nifti.get_data()
elif isinstance(nifti, nipy.core.image.image.Image):
feature = nifti.get_data()
font = {'size': 8}
rc('font', **font)
coords = cluster_dict['top_clust']['coords']
if coords == None:
return
feature /= feature.std()
imax = np.max(np.absolute(feature))
imin = -imax
imshow_args = dict(
vmax=imax,
vmin=imin,
alpha=0.7
)
coords = ([-coords[0], -coords[1], coords[2]])
plt.axis('off')
plt.text(0.05, 0.8, text, horizontalalignment='center',
color=(texcol, texcol, texcol))
try:
plot_map(feature,
xyz_affine(nifti),
anat=anat.get_data(),
anat_affine=xyz_affine(anat),
threshold=image_threshold,
cut_coords=coords,
annotate=False,
cmap=cmap,
draw_cross=False,
**imshow_args)
except Exception as e:
return
plt.savefig(out_path, transparent=True, facecolor=(bgcol, bgcol, bgcol))
def plot_colat_slice(self, component, colat, valmin, valmax, iteration=0, verbose=True):
#- Some initialisations. ------------------------------------------------------------------
colat = np.pi * colat / 180.0
n_procs = self.setup["procs"]["px"] * self.setup["procs"]["py"] * self.setup["procs"]["pz"]
vmax = float("-inf")
vmin = float("inf")
fig, ax = plt.subplots()
#- Loop over processor boxes and check if colat falls within the volume. ------------------
for p in range(n_procs):
if (colat >= self.theta[p,:].min()) & (colat <= self.theta[p,:].max()):
#- Read this field and make lats & lons. ------------------------------------------
field = self.read_single_box(component,p,iteration)
r, lon = np.meshgrid(self.z[p,:], self.phi[p,:])
x = r * np.cos(lon)
y = r * np.sin(lon)
#- Find the colat index and plot for this one box. --------------------------------
idx=min(np.where(min(np.abs(self.theta[p,:]-colat))==np.abs(self.theta[p,:]-colat))[0])
colat_effective = self.theta[p,idx]*180.0/np.pi
#- Find min and max values. -------------------------------------------------------
vmax = max(vmax, field[idx,:,:].max())
vmin = min(vmin, field[idx,:,:].min())
#- Make a nice colourmap and plot. ------------------------------------------------
my_colormap=cm.make_colormap({0.0:[0.1,0.0,0.0], 0.2:[0.8,0.0,0.0], 0.3:[1.0,0.7,0.0],0.48:[0.92,0.92,0.92], 0.5:[0.92,0.92,0.92], 0.52:[0.92,0.92,0.92], 0.7:[0.0,0.6,0.7], 0.8:[0.0,0.0,0.8], 1.0:[0.0,0.0,0.1]})
cax = ax.pcolor(x, y, field[idx,:,:], cmap=my_colormap, vmin=valmin,vmax=valmax)
#- Add colobar and title. ------------------------------------------------------------------
cbar = fig.colorbar(cax)
if component in UNIT_DICT:
cb.set_label(UNIT_DICT[component], fontsize="x-large", rotation=0)
plt.suptitle("Vertical slice of %s at %i degree colatitude" % (component, colat_effective), size="large")
plt.axis('equal')
plt.show()
#==============================================================================================
#- Plot depth slice.
#==============================================================================================
def double_lorentzian_fixedsplitting_testing():
# This method does not work and has to be fixed!!!
for ii in range(1):
# time.sleep(0.51)
start=2800
stop=2950
num_points=int((stop-start)/2)
x = np.linspace(start, stop, num_points)
mod,params = qudi_fitting.make_multiplelorentzian_model(no_of_lor=2)
p=Parameters()
#============ Create data ==========
p.add('c',value=100)
p.add('lorentz0_amplitude',value=-abs(np.random.random(1)*50+100))
p.add('lorentz0_center',value=np.random.random(1)*150.0+2800)
p.add('lorentz0_sigma',value=abs(np.random.random(1)*2.+1.))
p.add('lorentz1_center',value=p['lorentz0_center']+20)
p.add('lorentz1_sigma',value=abs(np.random.random(1)*2.+1.))
p.add('lorentz1_amplitude',value=-abs(np.random.random(1)*50+100))
data_noisy=(mod.eval(x=x,params=p)
+ 2*np.random.normal(size=x.shape))
para=Parameters()
result=qudi_fitting.make_doublelorentzian_fit(axis=x,data=data_noisy,add_parameters=para)
data_smooth, offset = qudi_fitting.find_offset_parameter(x,data_noisy)
data_level=data_smooth-offset
#search for double lorentzian
error, \
sigma0_argleft, dip0_arg, sigma0_argright, \
sigma1_argleft, dip1_arg , sigma1_argright = \
qudi_fitting._search_double_dip(x, data_level,make_prints=False)
print(x[sigma0_argleft], x[dip0_arg], x[sigma0_argright], x[sigma1_argleft], x[dip1_arg], x[sigma1_argright])
print(x[dip0_arg], x[dip1_arg])
plt.plot((x[sigma0_argleft], x[sigma0_argleft]), ( data_noisy.min() ,data_noisy.max()), 'b-')
plt.plot((x[sigma0_argright], x[sigma0_argright]), (data_noisy.min() ,data_noisy.max()), 'b-')
plt.plot((x[sigma1_argleft], x[sigma1_argleft]), ( data_noisy.min() ,data_noisy.max()), 'k-')
plt.plot((x[sigma1_argright], x[sigma1_argright]), ( data_noisy.min() ,data_noisy.max()), 'k-')
try:
plt.plot(x,data_noisy,'o')
plt.plot(x,result.init_fit,'-y')
plt.plot(x,result.best_fit,'-r',linewidth=2.0,)
plt.plot(x,data_smooth,'-g')
except:
print('exception')
plt.show()
def compute_inv_dft(x_val, y_val, zeropad_num=0):
""" Compute the inverse Discrete fourier Transform
@param numpy.array x_val: 1D array
@param numpy.array y_val: 1D array of same size as x_val
@param int zeropad_num: zeropadding (adding zeros to the end of the array).
zeropad_num >= 0, the size of the array, which is
add to the end of the y_val before performing the
dft. The resulting array will have the length
(len(y_val)/2)*(zeropad_num+1)
Note that zeropadding will not change or add more
information to the dft, it will solely interpolate
between the dft_y values.
@return: tuple(dft_x, dft_y):
be aware that the return arrays' length depend on the zeropad
number like
len(dft_x) = len(dft_y) = (len(y_val)/2)*(zeropad_num+1)
"""
x_val = np.array(x_val)
y_val = np.array(y_val)
corrected_y = np.abs(y_val - y_val.max())
# The absolute values contain the fourier transformed y values
# zeropad_arr = np.zeros(len(corrected_y)*(zeropad_num+1))
# zeropad_arr[:len(corrected_y)] = corrected_y
fft_y = np.fft.ifft(corrected_y)
# Due to the sampling theorem you can only identify frequencies at half
# of the sample rate, therefore the FT contains an almost symmetric
# spectrum (the asymmetry results from aliasing effects). Therefore take
# the half of the values for the display.
middle = int((len(corrected_y)+1)//2)
# sample spacing of x_axis, if x is a time axis than it corresponds to a
# timestep:
x_spacing = np.round(x_val[-1] - x_val[-2], 12)
# use the helper function of numpy to calculate the x_values for the
# fourier space. That function will handle an occuring devision by 0:
fft_x = np.fft.fftfreq(len(corrected_y), d=x_spacing)
# return abs(fft_x[:middle]), fft_y[:middle]
return fft_x[:middle], fft_y[:middle]
# return fft_x, fft_y
def illustrate(Colors=Colors):
if hasattr(Colors, 'colors'):
Colors = Colors.colors
from matplotlib import pylab
rcParams = pylab.rcParams
rcParams['pdf.fonttype'] = 42
rcParams['ps.fonttype'] = 42
rcParams['text.usetex'] = False
rcParams['xtick.labelsize'] = 20
rcParams['ytick.labelsize'] = 20
rcParams['legend.fontsize'] = 25
import bnpy
Data = get_data(T=1000, nDocTotal=8)
for k in xrange(K):
zmask = Data.TrueParams['Z'] == k
pylab.plot(Data.X[zmask, 0], Data.X[zmask, 1], '.', color=Colors[k],
markeredgecolor=Colors[k],
alpha=0.4)
sigEdges = np.flatnonzero(transPi[k] > 0.0001)
for j in sigEdges:
if j == k:
continue
dx = mus[j, 0] - mus[k, 0]
dy = mus[j, 1] - mus[k, 1]
pylab.arrow(mus[k, 0], mus[k, 1],
0.8 * dx,
0.8 * dy,
head_width=2, head_length=4,
facecolor=Colors[k], edgecolor=Colors[k])
tx = 0 - mus[k, 0]
ty = 0 - mus[k, 1]
xy = (mus[k, 0] - 0.2 * tx, mus[k, 1] - 0.2 * ty)
'''
pylab.annotate( u'\u27F2',
xy=(mus[k,0], mus[k,1]),
color=Colors[k],
fontsize=35,
)
'''
pylab.gca().yaxis.set_ticks_position('left')
pylab.gca().xaxis.set_ticks_position('bottom')
pylab.axis('image')
pylab.ylim([-38, 38])
pylab.xlim([-38, 38])
