def test_discretization(nmpc, nb_steps):
dT = nmpc.preview.dT
pylab.ion()
pylab.clf()
ax = pylab.subplot(311)
ax.set_color_cycle(['r', 'g', 'b'])
pylab.plot(
[sum(dT[:i]) for i in xrange(len(dT))],
nmpc.preview.P, marker='o')
pylab.plot(
pylab.linspace(0., sum(dT), nb_steps + 1),
[x[0:3] for x in nmpc.preview.discretize(nb_steps)],
marker='s', linestyle='--')
ax = pylab.subplot(312)
ax.set_color_cycle(['r', 'g', 'b'])
pylab.plot(
[sum(dT[:i]) for i in xrange(len(dT))],
nmpc.preview.V, marker='o')
pylab.plot(
pylab.linspace(0., sum(dT), nb_steps + 1),
[x[3:6] for x in nmpc.preview.discretize(nb_steps)],
marker='s', linestyle='--')
ax = pylab.subplot(313)
ax.set_color_cycle(['r', 'g', 'b'])
pylab.plot(
[sum(dT[:i]) for i in xrange(len(dT))],
nmpc.preview.Z, marker='o')
pylab.plot(
pylab.linspace(0., sum(dT), nb_steps + 1),
[x[6:9] for x in nmpc.preview.discretize(nb_steps)],
marker='s', linestyle='--')
python类linspace()的实例源码
def drawHealpixMap(map, lon, lat, size=1.0, xsize=501, coord='GC', **kwargs):
"""
Draw local projection of healpix map.
"""
ax = plt.gca()
x = np.linspace(-size,size,xsize)
y = np.linspace(-size,size,xsize)
xx, yy = np.meshgrid(x,y)
coord = coord.upper()
if coord == 'GC':
#Assumes map and (lon,lat) are Galactic, but plotting celestial
llon, llat = image2sphere(*gal2cel(lon,lat),x=xx.flat,y=yy.flat)
pix = ang2pix(healpy.get_nside(map),*cel2gal(llon,llat))
elif coord == 'CG':
#Assumes map and (lon,lat) are celestial, but plotting Galactic
llon, llat = image2sphere(*cel2gal(lon,lat),x=xx.flat,y=yy.flat)
pix = ang2pix(healpy.get_nside(map),*gal2cel(llon,llat))
else:
#Assumes plotting the native coordinates
llon, llat = image2sphere(lon,lat,xx.flat,yy.flat)
pix = ang2pix(healpy.get_nside(map),llon,llat)
values = map[pix].reshape(xx.shape)
zz = np.ma.array(values,mask=(values==healpy.UNSEEN),fill_value=np.nan)
return drawProjImage(xx,yy,zz,coord=coord,**kwargs)
def drawImage(self,ax=None,invert=True):
if not ax: ax = plt.gca()
if self.config['data']['survey']=='sdss':
# Optical Image
im = ugali.utils.plotting.getSDSSImage(**self.image_kwargs)
# Flipping JPEG:
# https://github.com/matplotlib/matplotlib/issues/101
im = im[::-1]
ax.annotate("SDSS Image",**self.label_kwargs)
else:
im = ugali.utils.plotting.getDSSImage(**self.image_kwargs)
im = im[::-1,::-1]
ax.annotate("DSS Image",**self.label_kwargs)
size=self.image_kwargs.get('radius',1.0)
# Celestial coordinates
x = np.linspace(-size,size,im.shape[0])
y = np.linspace(-size,size,im.shape[1])
xx, yy = np.meshgrid(x,y)
#kwargs = dict(cmap='gray',interpolation='none')
kwargs = dict(cmap='gray',coord='C')
im = drawProjImage(xx,yy,im,**kwargs)
try: plt.gcf().delaxes(ax.cax)
except AttributeError: pass
return im
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 drawMembersCMD(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
mmin, mmax = 16., 24.
cmin, cmax = -0.5, 1.0
mbins = np.linspace(mmin, mmax, 150)
cbins = np.linspace(cmin, cmax, 150)
mag_1 = data[self.config['catalog']['mag_1_field']]
mag_2 = data[self.config['catalog']['mag_2_field']]
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(data['COLOR'][sel][~sel_prob], mag_1[sel][~sel_prob],
marker='o',s=2,c='0.75',edgecolor='none')
sc = pylab.scatter(data['COLOR'][sel][sel_prob][index_sort], mag_1[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)
pylab.xlim(cmin, cmax)
pylab.ylim(mmax, mmin)
pylab.xlabel(r'$g - r$')
pylab.ylabel(r'$g$')
#axes[1].yaxis.set_major_locator(MaxNLocator(prune='lower'))
pylab.xticks([-0.5, 0., 0.5, 1.])
pylab.yticks(numpy.arange(mmax - 1., mmin - 1., -1.))
ugali.utils.plotting.drawIsochrone(self.isochrone, c='k', zorder=10)
pylab.text(0.05, 0.95, r'$\Sigma p_{i} = %i$'%(data['PROB'].sum()),
fontsize=10, horizontalalignment='left', verticalalignment='top', color='k', transform=pylab.gca().transAxes,
bbox=dict(facecolor='white', alpha=1., edgecolor='none'))
divider = make_axes_locatable(pylab.gca())
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 drawKernelHist(ax, kernel):
ext = kernel.extension
theta = kernel.theta
lon, lat = kernel.lon, kernel.lat
xmin,xmax = -5*ext,5*ext
ymin,ymax = -5*ext,5*ext,
x = np.linspace(xmin,xmax,100)+kernel.lon
y = np.linspace(ymin,ymax,100)+kernel.lat
xx,yy = np.meshgrid(x,y)
zz = kernel.pdf(xx,yy)
im = ax.imshow(zz)#,extent=[xmin,xmax,ymin,ymax])
hax,vax = draw_slices(ax,zz,color='k')
mc_lon,mc_lat = kernel.sample(1e5)
hist,xedges,yedges = np.histogram2d(mc_lon,mc_lat,bins=[len(x),len(y)],
range=[[x.min(),x.max()],[y.min(),y.max()]])
xbins,ybins = np.arange(hist.shape[0])+0.5,np.arange(hist.shape[1])+0.5
vzz = zz.sum(axis=0)
hzz = zz.sum(axis=1)
vmc = hist.sum(axis=0)
hmc = hist.sum(axis=1)
vscale = vzz.max()/vmc.max()
hscale = hzz.max()/hmc.max()
kwargs = dict(marker='.',ls='',color='r')
hax.errorbar(hmc*hscale, ybins, xerr=np.sqrt(hmc)*hscale,**kwargs)
vax.errorbar(xbins, vmc*vscale,yerr=np.sqrt(vmc)*vscale,**kwargs)
ax.set_ylim(0,len(y))
ax.set_xlim(0,len(x))
#try: ax.cax.colorbar(im)
#except: pylab.colorbar(im)
#a0 = np.array([0.,0.])
#a1 =kernel.a*np.array([np.sin(np.deg2rad(theta)),-np.cos(np.deg2rad(theta))])
#ax.plot([a0[0],a1[0]],[a0[1],a1[1]],'-ob')
#
#b0 = np.array([0.,0.])
#b1 =kernel.b*np.array([np.cos(np.radians(theta)),np.sin(np.radians(theta))])
#ax.plot([b0[0],b1[0]],[b0[1],b1[1]],'-or')
label_kwargs = dict(xy=(0.05,0.05),xycoords='axes fraction', xytext=(0, 0),
textcoords='offset points',ha='left', va='bottom',size=10,
bbox={'boxstyle':"round",'fc':'1'}, zorder=10)
norm = zz.sum() * (x[1]-x[0])**2
ax.annotate("Sum = %.2f"%norm,**label_kwargs)
#ax.set_xlabel(r'$\Delta$ LON (deg)')
#ax.set_ylabel(r'$\Delta$ LAT (deg)')
###################################################
def drawKernel(kernel, contour=False, coords='C', **kwargs):
ax = plt.gca()
if 'colors' not in kwargs:
kwargs.setdefault('cmap',matplotlib.cm.jet)
kwargs.setdefault('origin','lower')
ext = kernel.extension
theta = kernel.theta
xmin,xmax = -kernel.edge,kernel.edge
ymin,ymax = -kernel.edge,kernel.edge
if coords[-1] == 'G':
lon, lat = kernel.lon, kernel.lat
elif coords[-1] == 'C':
lon,lat = gal2cel(kernel.lon, kernel.lat)
else:
msg = 'Unrecognized coordinate: %s'%coords
raise Exception(msg)
x = np.linspace(xmin,xmax,500)+lon
y = np.linspace(ymin,ymax,500)+lat
xx,yy = np.meshgrid(x,y)
extent = [x[0],x[-1],y[0],y[-1]]
kwargs.setdefault('extent',extent)
if coords[-1] == 'C': xx,yy = cel2gal(xx,yy)
zz = kernel.pdf(xx.flat,yy.flat).reshape(xx.shape)
zmax = zz.max()
if contour:
levels = kwargs.pop('levels',10)
#levels = np.logspace(np.log10(zmax)-1,np.log10(zmax),7)
ret = ax.contour(zz,levels,**kwargs)
else:
val = np.ma.array(zz,mask=zz<zz.max()/100.)
ret = ax.imshow(val,**kwargs)
return ret
###################################################
