def cmPlot(targ_ra, targ_dec, data, iso, g_radius, nbhd, type):
"""Color-magnitude plot"""
angsep = ugali.utils.projector.angsep(targ_ra, targ_dec, data['RA'], data['DEC'])
annulus = (angsep > g_radius) & (angsep < 1.)
mag_g = data[mag_g_dred_flag]
mag_r = data[mag_r_dred_flag]
if type == 'stars':
filter = star_filter(data)
plt.title('Stellar Color-Magnitude')
elif type == 'galaxies':
filter = galaxy_filter(data)
plt.title('Galactic Color-Magnitude')
iso_filter = (iso.separation(mag_g, mag_r) < 0.1)
# Plot background objects
plt.scatter(mag_g[filter & annulus] - mag_r[filter & annulus], mag_g[filter & annulus], c='k', alpha=0.1, edgecolor='none', s=1)
# Plot isochrone
ugali.utils.plotting.drawIsochrone(iso, lw=2, label='{} Gyr, z = {}'.format(iso.age, iso.metallicity))
# Plot objects in nbhd
plt.scatter(mag_g[filter & nbhd] - mag_r[filter & nbhd], mag_g[filter & nbhd], c='g', s=5, label='r < {:.3f}$^\circ$'.format(g_radius))
# Plot objects in nbhd and near isochrone
plt.scatter(mag_g[filter & nbhd & iso_filter] - mag_r[filter & nbhd & iso_filter], mag_g[filter & nbhd & iso_filter], c='r', s=5, label='$\Delta$CM < 0.1')
plt.axis([-0.5, 1, 16, 24])
plt.gca().invert_yaxis()
plt.gca().set_aspect(1./4.)
plt.legend(loc='upper left')
plt.xlabel('g-r (mag)')
plt.ylabel('g (mag)')
python类scatter()的实例源码
def projScatter(lon, lat, **kwargs):
"""
Create a scatter plot on HEALPix projected axes.
Inputs: lon (deg), lat (deg)
"""
healpy.projscatter(lon, lat, lonlat=True, **kwargs)
############################################################
def drawSpatial(self, ax=None):
if not ax: ax = plt.gca()
# Stellar Catalog
self._create_catalog()
cut = (self.catalog.color > 0) & (self.catalog.color < 1)
catalog = self.catalog.applyCut(cut)
ax.scatter(catalog.lon,catalog.lat,c='k',marker='.',s=1)
ax.set_xlim(self.glon-0.5,self.glon+0.5)
ax.set_ylim(self.glat-0.5,self.glat+0.5)
ax.set_xlabel('GLON (deg)')
ax.set_ylabel('GLAT (deg)')
def drawHessDiagram(self,catalog=None):
ax = plt.gca()
if not catalog: catalog = self.get_stars()
r_peak = self.kernel.extension
angsep = ugali.utils.projector.angsep(self.ra, self.dec, catalog.ra, catalog.dec)
cut_inner = (angsep < r_peak)
cut_annulus = (angsep > 0.5) & (angsep < 1.) # deg
mmin, mmax = 16., 24.
cmin, cmax = -0.5, 1.0
mbins = np.linspace(mmin, mmax, 150)
cbins = np.linspace(cmin, cmax, 150)
color = catalog.color[cut_annulus]
mag = catalog.mag[cut_annulus]
h, xbins, ybins = numpy.histogram2d(color, mag, bins=[cbins,mbins])
blur = nd.filters.gaussian_filter(h.T, 2)
kwargs = dict(extent=[xbins.min(),xbins.max(),ybins.min(),ybins.max()],
cmap='gray_r', aspect='auto', origin='lower',
rasterized=True, interpolation='none')
ax.imshow(blur, **kwargs)
pylab.scatter(catalog.color[cut_inner], catalog.mag[cut_inner],
c='red', s=7, edgecolor='none')# label=r'$r < %.2f$ deg'%(r_peak))
ugali.utils.plotting.drawIsochrone(self.isochrone, c='b', zorder=10)
ax.set_xlim(-0.5, 1.)
ax.set_ylim(24., 16.)
plt.xlabel(r'$g - r$')
plt.ylabel(r'$g$')
plt.xticks([-0.5, 0., 0.5, 1.])
plt.yticks(numpy.arange(mmax - 1., mmin - 1., -1.))
radius_string = (r'${\rm r}<%.1f$ arcmin'%( 60 * r_peak))
pylab.text(0.05, 0.95, radius_string,
fontsize=10, ha='left', va='top', color='red',
transform=pylab.gca().transAxes,
bbox=dict(facecolor='white', alpha=1., edgecolor='none'))
def drawMembersSpatial(self,data):
ax = plt.gca()
if isinstance(data,basestring):
filename = data
data = pyfits.open(filename)[1].data
xmin, xmax = -0.25,0.25
ymin, ymax = -0.25,0.25
xx,yy = np.meshgrid(np.linspace(xmin,xmax),np.linspace(ymin,ymax))
x_prob, y_prob = sphere2image(self.ra, self.dec, data['RA'], data['DEC'])
sel = (x_prob > xmin)&(x_prob < xmax) & (y_prob > ymin)&(y_prob < ymax)
sel_prob = data['PROB'][sel] > 5.e-2
index_sort = numpy.argsort(data['PROB'][sel][sel_prob])
plt.scatter(x_prob[sel][~sel_prob], y_prob[sel][~sel_prob],
marker='o', s=2, c='0.75', edgecolor='none')
sc = plt.scatter(x_prob[sel][sel_prob][index_sort],
y_prob[sel][sel_prob][index_sort],
c=data['PROB'][sel][sel_prob][index_sort],
marker='o', s=10, edgecolor='none', cmap='jet', vmin=0., vmax=1.) # Spectral_r
drawProjImage(xx,yy,None,coord='C')
#ax.set_xlim(xmax, xmin)
#ax.set_ylim(ymin, ymax)
#plt.xlabel(r'$\Delta \alpha_{2000}\,(\deg)$')
#plt.ylabel(r'$\Delta \delta_{2000}\,(\deg)$')
plt.xticks([-0.2, 0., 0.2])
plt.yticks([-0.2, 0., 0.2])
divider = make_axes_locatable(ax)
ax_cb = divider.new_horizontal(size="7%", pad=0.1)
plt.gcf().add_axes(ax_cb)
pylab.colorbar(sc, cax=ax_cb, orientation='vertical', ticks=[0, 0.2, 0.4, 0.6, 0.8, 1.0], label='Membership Probability')
ax_cb.yaxis.tick_right()
def profileUpperLimit(self, delta = 2.71):
"""
Compute one-sided upperlimit via profile method.
"""
a = self.p_2
b = self.p_1
if self.vertex_x < 0:
c = self.p_0 + delta
else:
c = self.p_0 - self.vertex_y + delta
if b**2 - 4. * a * c < 0.:
print 'WARNING'
print a, b, c
#pylab.figure()
#pylab.scatter(self.x, self.y)
#raw_input('WAIT')
return 0.
return max((numpy.sqrt(b**2 - 4. * a * c) - b) / (2. * a), (-1. * numpy.sqrt(b**2 - 4. * a * c) - b) / (2. * a))
#def bayesianUpperLimit3(self, alpha, steps = 1.e5):
# """
# Compute one-sided upper limit using Bayesian Method of Helene.
# """
# # Need a check to see whether limit is reliable
# pdf = scipy.interpolate.interp1d(self.x, numpy.exp(self.y / 2.)) # Convert from 2 * log(likelihood) to likelihood
# x_pdf = numpy.linspace(self.x[0], self.x[-1], steps)
# cdf = numpy.cumsum(pdf(x_pdf))
# cdf /= cdf[-1]
# cdf_reflect = scipy.interpolate.interp1d(cdf, x_pdf)
# return cdf_reflect(alpha)
# #return self.x[numpy.argmin((cdf - alpha)**2)]
def bayesianUpperLimit(self, alpha, steps=1.e5, plot=False):
"""
Compute one-sided upper limit using Bayesian Method of Helene.
Several methods of increasing numerical stability have been implemented.
"""
x_dense, y_dense = self.densify()
y_dense -= numpy.max(y_dense) # Numeric stability
f = scipy.interpolate.interp1d(x_dense, y_dense, kind='linear')
x = numpy.linspace(0., numpy.max(x_dense), steps)
pdf = numpy.exp(f(x) / 2.)
cut = (pdf / numpy.max(pdf)) > 1.e-10
x = x[cut]
pdf = pdf[cut]
#pdf /= pdf[0]
#forbidden = numpy.nonzero(pdf < 1.e-10)[0]
#if len(forbidden) > 0:
# index = forbidden[0] # Numeric stability
# x = x[0: index]
# pdf = pdf[0: index]
cdf = numpy.cumsum(pdf)
cdf /= cdf[-1]
cdf_reflect = scipy.interpolate.interp1d(cdf, x)
#if plot:
# pylab.figure()
# pylab.plot(x, f(x))
# pylab.scatter(self.x, self.y, c='red')
#
# pylab.figure()
# pylab.plot(x, pdf)
#
# pylab.figure()
# pylab.plot(cdf, x)
return cdf_reflect(alpha)
def confidenceInterval(self, alpha=0.6827, steps=1.e5, plot=False):
"""
Compute two-sided confidence interval by taking x-values corresponding to the largest PDF-values first.
"""
x_dense, y_dense = self.densify()
y_dense -= numpy.max(y_dense) # Numeric stability
f = scipy.interpolate.interp1d(x_dense, y_dense, kind='linear')
x = numpy.linspace(0., numpy.max(x_dense), steps)
# ADW: Why does this start at 0, which often outside the input range?
# Wouldn't starting at xmin be better:
#x = numpy.linspace(numpy.min(x_dense), numpy.max(x_dense), steps)
pdf = numpy.exp(f(x) / 2.)
cut = (pdf / numpy.max(pdf)) > 1.e-10
x = x[cut]
pdf = pdf[cut]
sorted_pdf_indices = numpy.argsort(pdf)[::-1] # Indices of PDF in descending value
cdf = numpy.cumsum(pdf[sorted_pdf_indices])
cdf /= cdf[-1]
sorted_pdf_index_max = numpy.argmin((cdf - alpha)**2)
x_select = x[sorted_pdf_indices[0: sorted_pdf_index_max]]
#if plot:
# cdf = numpy.cumsum(pdf)
# cdf /= cdf[-1]
# print cdf[numpy.max(sorted_pdf_indices[0: sorted_pdf_index_max])] \
# - cdf[numpy.min(sorted_pdf_indices[0: sorted_pdf_index_max])]
#
# pylab.figure()
# pylab.plot(x, f(x))
# pylab.scatter(self.x, self.y, c='red')
#
# pylab.figure()
# pylab.plot(x, pdf)
return numpy.min(x_select), numpy.max(x_select)
############################################################
def plot_candidates(self):
"""Plot a representation of candidate periodicity
Size gives the periodicity strength, color the order of preference
"""
hues = np.arange(self.ncand)/float(self.ncand)
hsv = np.swapaxes(np.atleast_3d([[hues,np.ones(len(hues)),np.ones(len(hues))]]),1,2)
cols = hsv_to_rgb(hsv).squeeze()
for per in self.periods:
nc = len(per.cand_period)
pl.scatter(per.time*np.ones(nc),per.cand_period,s=per.cand_strength*100,c=cols[0:nc],alpha=.5)
def plot_scatter(data, dir=None, filename="scatter", color="blue"):
if dir is None:
raise Exception()
try:
os.mkdir(dir)
except:
pass
fig = pylab.gcf()
fig.set_size_inches(16.0, 16.0)
pylab.clf()
pylab.scatter(data[:, 0], data[:, 1], s=20, marker="o", edgecolors="none", color=color)
pylab.xlim(-4, 4)
pylab.ylim(-4, 4)
pylab.savefig("{}/{}.png".format(dir, filename))
def plot_scatter(data, dir=None, filename="scatter", color="blue"):
if dir is None:
raise Exception()
try:
os.mkdir(dir)
except:
pass
fig = pylab.gcf()
fig.set_size_inches(16.0, 16.0)
pylab.clf()
pylab.scatter(data[:, 0], data[:, 1], s=20, marker="o", edgecolors="none", color=color)
pylab.xlim(-4, 4)
pylab.ylim(-4, 4)
pylab.savefig("{}/{}".format(dir, filename))
def plot_scatter(data, dir=None, filename="scatter", color="blue"):
if dir is None:
raise Exception()
try:
os.mkdir(dir)
except:
pass
fig = pylab.gcf()
fig.set_size_inches(16.0, 16.0)
pylab.clf()
pylab.scatter(data[:, 0], data[:, 1], s=20, marker="o", edgecolors="none", color=color)
pylab.xlim(-4, 4)
pylab.ylim(-4, 4)
pylab.savefig("{}/{}.png".format(dir, filename))
def plot_scatter(data, dir=None, filename="scatter", color="blue"):
if dir is None:
raise Exception()
try:
os.mkdir(dir)
except:
pass
fig = pylab.gcf()
fig.set_size_inches(16.0, 16.0)
pylab.clf()
pylab.scatter(data[:, 0], data[:, 1], s=20, marker="o", edgecolors="none", color=color)
pylab.xlim(-4, 4)
pylab.ylim(-4, 4)
pylab.savefig("{}/{}".format(dir, filename))
def visualize_z(z_batch, dir=None):
if dir is None:
raise Exception()
try:
os.mkdir(dir)
except:
pass
fig = pylab.gcf()
fig.set_size_inches(20.0, 16.0)
pylab.clf()
for n in xrange(z_batch.shape[0]):
result = pylab.scatter(z_batch[n, 0], z_batch[n, 1], s=40, marker="o", edgecolors='none')
pylab.xlabel("z1")
pylab.ylabel("z2")
pylab.savefig("%s/latent_code.png" % dir)
def draw2D(self):
for i in xrange(self.nComponents):
xeq = lambda t: self.params[6 * i + 3] * np.cos(t) * np.cos(self.params[6 * i + 5]) + self.params[
6 * i + 4] * np.sin(
t) * np.sin(self.params[6 * i + 5]) + self.params[6 * i + 1]
yeq = lambda t: - self.params[6 * i + 3] * np.cos(t) * np.sin(self.params[6 * i + 5]) + self.params[
6 * i + 4] * np.sin(
t) * np.cos(self.params[6 * i + 5]) + self.params[6 * i + 2]
t = np.linspace(0, 2 * np.pi, 100)
x = xeq(t)
y = yeq(t)
pylab.scatter(self.params[6 * i + 2], self.params[6 * i + 1], color='k')
pylab.plot(y.astype(int), x.astype(int), self.colors[i] + '-')
def plot_gaussians3D(self, save=False, titlehist='', pathfig='', newfig=True):
ax = extract.hist2d(titlehist, newfig=newfig)
dx, dy = np.indices(self.shape)
for n in xrange(0, len(self.params), 6):
gaussunitaire = GaussianForFit(self.image, 1, params=self.params[n:n + 6])
ax.scatter(gaussunitaire.params[1], gaussunitaire.params[2],
self.image[gaussunitaire.params[1], gaussunitaire.params[2]], color=self.colors[n % 5],
label="{0:.3f}".format(gaussunitaire.params[0]), alpha=0.7)
ax.contour(dx, dy, gaussunitaire.gaussian, colors=self.colors[n % 5])
if save:
pylab.savefig(pathfig)
def patience(log, ax=None):
ax = ax or pl.gca()
maxes = running_max(list(log.iteration),
list(log.dev_accuracy - log.tradeoff * log.dev_runtime))
ax.scatter(maxes[:,0], maxes[:,1], lw=0)
def view_masks(file_name, t_start=0, t_stop=1, n_elec=0):
params = CircusParser(file_name)
data_file = params.get_data_file()
data_file.open()
N_e = params.getint('data', 'N_e')
N_t = params.getint('detection', 'N_t')
N_total = params.nb_channels
sampling_rate = params.rate
do_temporal_whitening = params.getboolean('whitening', 'temporal')
do_spatial_whitening = params.getboolean('whitening', 'spatial')
spike_thresh = params.getfloat('detection', 'spike_thresh')
file_out_suff = params.get('data', 'file_out_suff')
nodes, edges = get_nodes_and_edges(params)
chunk_size = (t_stop - t_start)*sampling_rate
padding = (t_start*sampling_rate, t_start*sampling_rate)
inv_nodes = numpy.zeros(N_total, dtype=numpy.int32)
inv_nodes[nodes] = numpy.argsort(nodes)
safety_time = params.getint('clustering', 'safety_time')
if do_spatial_whitening:
spatial_whitening = load_data(params, 'spatial_whitening')
if do_temporal_whitening:
temporal_whitening = load_data(params, 'temporal_whitening')
thresholds = load_data(params, 'thresholds')
data = data_file.get_data(0, chunk_size, padding=padding, nodes=nodes)
data_shape = len(data)
data_file.close()
peaks = {}
indices = inv_nodes[edges[nodes[n_elec]]]
if do_spatial_whitening:
data = numpy.dot(data, spatial_whitening)
if do_temporal_whitening:
data = scipy.ndimage.filters.convolve1d(data, temporal_whitening, axis=0, mode='constant')
for i in xrange(N_e):
peaks[i] = algo.detect_peaks(data[:, i], thresholds[i], valley=True, mpd=0)
pylab.figure()
for count, i in enumerate(indices):
pylab.plot(count*5 + data[:, i], '0.25')
#xmin, xmax = pylab.xlim()
pylab.scatter(peaks[i], count*5 + data[peaks[i], i], s=10, c='r')
for count, i in enumerate(peaks[n_elec]):
pylab.axvspan(i - safety_time, i + safety_time, facecolor='r', alpha=0.5)
pylab.ylim(-5, len(indices)*5 )
pylab.xlabel('Time [ms]')
pylab.ylabel('Electrode')
pylab.tight_layout()
pylab.setp(pylab.gca(), yticks=[])
pylab.show()
return peaks
def view_classification(data_1, data_2, title=None, save=None):
fig = pylab.figure()
count = 0
panels = [0, 2, 1, 3]
for item in [data_1, data_2]:
clf, cld, X, X_raw, y = item
for mode in ['predict', 'decision_function']:
ax = fig.add_subplot(2, 2, panels[count]+1)
if mode == 'predict':
c = clf
vmax = 1.0
vmin = 0.0
elif mode == 'decision_function':
c = cld
vmax = max(abs(numpy.amin(c)), abs(numpy.amax(c)))
vmin = - vmax
from circus.validating.utils import Projection
p = Projection()
_ = p.fit(X_raw, y)
X_raw_ = p.transform(X_raw)
# Plot figure.
sc = ax.scatter(X_raw_[:, 0], X_raw_[:, 1], c=c, s=5, lw=0.1, cmap='bwr',
vmin=vmin, vmax=vmax)
cb = fig.colorbar(sc)
ax.grid(True)
if panels[count] in [0, 1]:
if panels[count] == 0:
ax.set_title('Classification Before')
ax.set_ylabel("2nd component")
if panels[count] == 1:
ax.set_title('Classification After')
cb.set_label('Prediction')
elif panels[count] in [2, 3]:
ax.set_xlabel("1st component")
if panels[count] == 2:
ax.set_ylabel("2nd component")
if panels[count] == 3:
cb.set_label('Decision function')
count += 1
if save is None:
pylab.show()
else:
pylab.savefig(save)
pylab.close(fig)
return
def drawMembership(self, ax=None, radius=None, zidx=0, mc_source_id=1):
if not ax: ax = plt.gca()
import ugali.analysis.scan
filename = self.config.mergefile
logger.debug("Opening %s..."%filename)
f = pyfits.open(filename)
distance_modulus = f[2].data['DISTANCE_MODULUS'][zidx]
for ii, name in enumerate(self.config.params['isochrone']['infiles']):
logger.info('%s %s'%(ii, name))
isochrone = ugali.isochrone.Isochrone(self.config, name)
mag = isochrone.mag + distance_modulus
ax.scatter(isochrone.color,mag, color='0.5', s=800, zorder=0)
pix = ang2pix(self.nside, self.glon, self.glat)
likelihood_pix = ugali.utils.skymap.superpixel(pix,self.nside,self.config.params['coords']['nside_likelihood'])
config = self.config
scan = ugali.analysis.scan.Scan(self.config,likelihood_pix)
likelihood = scan.likelihood
distance_modulus_array = [self.config.params['scan']['distance_modulus_array'][zidx]]
likelihood.precomputeGridSearch(distance_modulus_array)
likelihood.gridSearch()
p = likelihood.membershipGridSearch()
sep = ugali.utils.projector.angsep(self.glon, self.glat, likelihood.catalog.lon, likelihood.catalog.lat)
radius = self.radius if radius is None else radius
cut = (sep < radius)
catalog = likelihood.catalog.applyCut(cut)
p = p[cut]
cut_mc_source_id = (catalog.mc_source_id == mc_source_id)
ax.scatter(catalog.color[cut_mc_source_id], catalog.mag[cut_mc_source_id], c='gray', s=100, edgecolors='none')
sc = ax.scatter(catalog.color, catalog.mag, c=p, edgecolors='none')
ax.set_xlim(likelihood.roi.bins_color[0], likelihood.roi.bins_color[-1])
ax.set_ylim(likelihood.roi.bins_mag[-1], likelihood.roi.bins_mag[0])
ax.set_xlabel('Color (mag)')
ax.set_ylabel('Magnitude (mag)')
try: ax.cax.colorbar(sc)
except: pylab.colorbar(sc)
