def get_local_wavenumbermesh(self, scaled=True, broadcast=False,
eliminate_highest_freq=False):
kx = fftfreq(self.N[0], 1./self.N[0])
ky = rfftfreq(self.N[1], 1./self.N[1])
if eliminate_highest_freq:
for i, k in enumerate((kx, ky)):
if self.N[i] % 2 == 0:
k[self.N[i]//2] = 0
Ks = np.meshgrid(kx, ky[self.rank*self.Np[1]//2:(self.rank*self.Np[1]//2+self.Npf)], indexing='ij', sparse=True)
if scaled is True:
Lp = 2*np.pi/self.L
Ks[0] *= Lp[0]
Ks[1] *= Lp[1]
K = Ks
if broadcast is True:
K = [np.broadcast_to(k, self.complex_shape()) for k in Ks]
return K
python类meshgrid()的实例源码
def normvectorfield(xs,ys,fs,**kw):
"""
plot normalized vector field
kwargs
======
- length is a desired length of the lines (default: 1)
- the rest of kwards are passed to plot
"""
length = kw.pop('length') if 'length' in kw else 1
x, y = np.meshgrid(xs, ys)
# calculate vector field
vx,vy = fs(x,y)
# plot vecor field
norm = length /np.sqrt(vx**2+vy**2)
plt.quiver(x, y, vx * norm, vy * norm, angles='xy',**kw)
def plot_interpolation(orderx,ordery):
s = PseudoSpectralDiscretization2D(orderx,XMIN,XMAX,
ordery,YMIN,YMAX)
Xc,Yc = s.get_x2d()
x = np.linspace(XMIN,XMAX,100)
y = np.linspace(YMIN,YMAX,100)
Xf,Yf = np.meshgrid(x,y,indexing='ij')
f_coarse = f(Xc,Yc)
f_interpolator = s.to_continuum(f_coarse)
f_num = f_interpolator(Xf,Yf)
plt.pcolor(Xf,Yf,f_num)
cb = plt.colorbar()
cb.set_label('interpolated function',fontsize=16)
plt.xlabel('x')
plt.ylabel('y')
for postfix in ['.png','.pdf']:
name = 'orthopoly_interpolated_function'+postfix
if USE_FIGS_DIR:
name = 'figs/' + name
plt.savefig(name,
bbox_inches='tight')
plt.clf()
def __init__(self,
orderx,xmin,xmax,
ordery,ymin,ymax):
"Constructor. Needs order and domain in x and y"
self.orderx,self.ordery = orderx,ordery
self.stencils = [PseudoSpectralDiscretization1D(orderx,xmin,xmax),
PseudoSpectralDiscretization1D(ordery,ymin,ymax)]
self.stencil_x,self.stencil_y = self.stencils
self.quads = [s.quads for s in self.stencils]
self.colocs = [s.colocation_points for s in self.stencils]
self.x,self.y = self.colocs
self.colocs2d = np.meshgrid(*self.colocs,indexing='ij')
self.X,self.Y = self.colocs2d
self.weights = [s.weights for s in self.stencils]
self.weights_x,self.weights_y = self.weights
self.weights2D = np.tensordot(*self.weights,axes=0)
def vectorfield(xs,ys,fs,**kw):
"""
plot vector field (no normalization!)
args
====
fs is a function that returns tuple (vx,vy)
kwargs
======
- length is a desired length of the lines (default: 1)
- the rest of kwards are passed to plot
"""
length= kw.pop('length') if 'length' in kw else 1
x, y=np.meshgrid(xs, ys)
# calculate vector field
vx,vy=fs(x,y)
# plot vecor field
norm = length
plt.quiver(x, y, vx * norm, vy * norm, angles='xy',**kw)
def create_reference_image(size, x0=10., y0=-3., sigma_x=50., sigma_y=30., dtype='float64',
reverse_xaxis=False, correct_axes=True, sizey=None, **kwargs):
"""
Creates a reference image: a gaussian brightness with elliptical
"""
inc_cos = np.cos(0./180.*np.pi)
delta_x = 1.
x = (np.linspace(0., size - 1, size) - size / 2.) * delta_x
if sizey:
y = (np.linspace(0., sizey-1, sizey) - sizey/2.) * delta_x
else:
y = x.copy()
if reverse_xaxis:
xx, yy = np.meshgrid(-x, y/inc_cos)
elif correct_axes:
xx, yy = np.meshgrid(-x, -y/inc_cos)
else:
xx, yy = np.meshgrid(x, y/inc_cos)
image = np.exp(-(xx-x0)**2./sigma_x - (yy-y0)**2./sigma_y)
return image.astype(dtype)
def generate_anchors_pre(height, width, feat_stride, anchor_scales=(8,16,32), anchor_ratios=(0.5,1,2)):
""" A wrapper function to generate anchors given different scales
Also return the number of anchors in variable 'length'
"""
anchors = generate_anchors(ratios=np.array(anchor_ratios), scales=np.array(anchor_scales))
A = anchors.shape[0]
shift_x = np.arange(0, width) * feat_stride
shift_y = np.arange(0, height) * feat_stride
shift_x, shift_y = np.meshgrid(shift_x, shift_y)
shifts = np.vstack((shift_x.ravel(), shift_y.ravel(), shift_x.ravel(), shift_y.ravel())).transpose()
K = shifts.shape[0]
# width changes faster, so here it is H, W, C
anchors = anchors.reshape((1, A, 4)) + shifts.reshape((1, K, 4)).transpose((1, 0, 2))
anchors = anchors.reshape((K * A, 4)).astype(np.float32, copy=False)
length = np.int32(anchors.shape[0])
return anchors, length
def plot2d_simplex(simplex, ind):
fig_dir = "./"
plt.cla()
n = 1000
x1 = np.linspace(-256, 1024, n)
x2 = np.linspace(-256, 1024, n)
X, Y = np.meshgrid(x1, x2)
Z = np.sqrt(X ** 2 + Y ** 2)
plt.contour(X, Y, Z, levels=list(np.arange(0, 1200, 10)))
plt.gca().set_aspect("equal")
plt.xlim((-256, 768))
plt.ylim((-256, 768))
plt.plot([simplex[0].x[0], simplex[1].x[0]],
[simplex[0].x[1], simplex[1].x[1]], color="#000000")
plt.plot([simplex[1].x[0], simplex[2].x[0]],
[simplex[1].x[1], simplex[2].x[1]], color="#000000")
plt.plot([simplex[2].x[0], simplex[0].x[0]],
[simplex[2].x[1], simplex[0].x[1]], color="#000000")
plt.savefig(os.path.join(fig_dir, "{:03d}.png".format(ind)))
def draw_hill(x,y):
a = np.linspace(-20,20,100)
print(a)
b = np.linspace(-20,20,100)
x = np.array(x)
y = np.array(y)
allSSE = np.zeros(shape=(len(a),len(b)))
for ai in range(0,len(a)):
for bi in range(0,len(b)):
a0 = a[ai]
b0 = b[bi]
SSE = calc_loss(a=a0,b=b0,x=x,y=y)
allSSE[ai][bi] = SSE
a,b = np.meshgrid(a, b)
return [a,b,allSSE]
def draw_hill(x,y):
a = np.linspace(-20,20,100)
print(a)
b = np.linspace(-20,20,100)
x = np.array(x)
y = np.array(y)
allSSE = np.zeros(shape=(len(a),len(b)))
for ai in range(0,len(a)):
for bi in range(0,len(b)):
a0 = a[ai]
b0 = b[bi]
SSE = calc_loss(a=a0,b=b0,x=x,y=y)
allSSE[ai][bi] = SSE
a,b = np.meshgrid(a, b)
return [a,b,allSSE]
# ????
def draw_hill(x,y):
a = np.linspace(-20,20,100)
print(a)
b = np.linspace(-20,20,100)
x = np.array(x)
y = np.array(y)
allSSE = np.zeros(shape=(len(a),len(b)))
for ai in range(0,len(a)):
for bi in range(0,len(b)):
a0 = a[ai]
b0 = b[bi]
SSE = calc_loss(a=a0,b=b0,x=x,y=y)
allSSE[ai][bi] = SSE
a,b = np.meshgrid(a, b)
return [a,b,allSSE]
def draw_hill(x,y):
a = np.linspace(-20,20,100)
print(a)
b = np.linspace(-20,20,100)
x = np.array(x)
y = np.array(y)
allSSE = np.zeros(shape=(len(a),len(b)))
for ai in range(0,len(a)):
for bi in range(0,len(b)):
a0 = a[ai]
b0 = b[bi]
SSE = calc_loss(a=a0,b=b0,x=x,y=y)
allSSE[ai][bi] = SSE
a,b = np.meshgrid(a, b)
return [a,b,allSSE]
# ????
def draw_hill(x,y):
a = np.linspace(-20,20,100)
print(a)
b = np.linspace(-20,20,100)
x = np.array(x)
y = np.array(y)
allSSE = np.zeros(shape=(len(a),len(b)))
for ai in range(0,len(a)):
for bi in range(0,len(b)):
a0 = a[ai]
b0 = b[bi]
SSE = calc_loss(a=a0,b=b0,x=x,y=y)
allSSE[ai][bi] = SSE
a,b = np.meshgrid(a, b)
return [a,b,allSSE]
# ????
def draw_hill(x,y):
a = np.linspace(-20,20,100)
print(a)
b = np.linspace(-20,20,100)
x = np.array(x)
y = np.array(y)
allSSE = np.zeros(shape=(len(a),len(b)))
for ai in range(0,len(a)):
for bi in range(0,len(b)):
a0 = a[ai]
b0 = b[bi]
SSE = calc_loss(a=a0,b=b0,x=x,y=y)
allSSE[ai][bi] = SSE
a,b = np.meshgrid(a, b)
return [a,b,allSSE]
# ????
def draw_hill(x,y):
a = np.linspace(-20,20,100)
print(a)
b = np.linspace(-20,20,100)
x = np.array(x)
y = np.array(y)
allSSE = np.zeros(shape=(len(a),len(b)))
for ai in range(0,len(a)):
for bi in range(0,len(b)):
a0 = a[ai]
b0 = b[bi]
SSE = calc_loss(a=a0,b=b0,x=x,y=y)
allSSE[ai][bi] = SSE
a,b = np.meshgrid(a, b)
return [a,b,allSSE]
# ????
def draw_hill(x,y):
a = np.linspace(-20,20,100)
print(a)
b = np.linspace(-20,20,100)
x = np.array(x)
y = np.array(y)
allSSE = np.zeros(shape=(len(a),len(b)))
for ai in range(0,len(a)):
for bi in range(0,len(b)):
a0 = a[ai]
b0 = b[bi]
SSE = calc_loss(a=a0,b=b0,x=x,y=y)
allSSE[ai][bi] = SSE
a,b = np.meshgrid(a, b)
return [a,b,allSSE]
# ????
def plot_grid2D(lons, lats, tec_grid2D, datetime, title_label = ''):
LATS, LONS = np.meshgrid(lats, lons)
m = Basemap(llcrnrlon=-180,
llcrnrlat=-55,
urcrnrlon=180,
urcrnrlat=75,
projection='merc',
area_thresh=1000,
resolution='i')
m.drawstates()
m.drawcountries()
m.drawcoastlines()
parallels = np.arange(-90,90,20)
m.drawparallels(parallels,labels=[True,False,False,True])
meridians = np.arange(0,360,40)
m.drawmeridians(meridians,labels=[True,False,False,True])
m.scatter(LONS, LATS, c=tec_grid2D, latlon = True, linewidths=0, s=5)
m.colorbar()
plt.title('%s\n%s' % (title_label, datetime.isoformat(' ')))
def _meshgrid(self, height, width):
with tf.variable_scope('_meshgrid'):
# This should be equivalent to:
# x_t, y_t = np.meshgrid(np.linspace(-1, 1, width),
# np.linspace(-1, 1, height))
# ones = np.ones(np.prod(x_t.shape))
# grid = np.vstack([x_t.flatten(), y_t.flatten(), ones])
x_t = tf.matmul(tf.ones(shape=tf.pack([height, 1])),
tf.transpose(tf.expand_dims(tf.linspace(-1.0, 1.0, width), 1), [1, 0]))
y_t = tf.matmul(tf.expand_dims(tf.linspace(-1.0, 1.0, height), 1),
tf.ones(shape=tf.pack([1, width])))
x_t_flat = tf.reshape(x_t, (1, -1))
y_t_flat = tf.reshape(y_t, (1, -1))
ones = tf.ones_like(x_t_flat)
grid = tf.concat(0, [x_t_flat, y_t_flat, ones])
return grid
def make_3d_mask(img_shape, center, radius, shape='sphere'):
mask = np.zeros(img_shape)
radius = np.rint(radius)
center = np.rint(center)
sz = np.arange(int(max(center[0] - radius, 0)), int(max(min(center[0] + radius + 1, img_shape[0]), 0)))
sy = np.arange(int(max(center[1] - radius, 0)), int(max(min(center[1] + radius + 1, img_shape[1]), 0)))
sx = np.arange(int(max(center[2] - radius, 0)), int(max(min(center[2] + radius + 1, img_shape[2]), 0)))
sz, sy, sx = np.meshgrid(sz, sy, sx)
if shape == 'cube':
mask[sz, sy, sx] = 1.
elif shape == 'sphere':
distance2 = ((center[0] - sz) ** 2
+ (center[1] - sy) ** 2
+ (center[2] - sx) ** 2)
distance_matrix = np.ones_like(mask) * np.inf
distance_matrix[sz, sy, sx] = distance2
mask[(distance_matrix <= radius ** 2)] = 1
elif shape == 'gauss':
z, y, x = np.ogrid[:mask.shape[0], :mask.shape[1], :mask.shape[2]]
distance = ((z - center[0]) ** 2 + (y - center[1]) ** 2 + (x - center[2]) ** 2)
mask = np.exp(- 1. * distance / (2 * radius ** 2))
mask[(distance > 3 * radius ** 2)] = 0
return mask
def run_numerical_test(function, n, k, results):
n_grid, k_grid = np.meshgrid(n, k)
if results.shape != n_grid.shape:
raise ValueError('results passed do not match the shape of n/k')
# Test both as scalars
for i, cur_n in enumerate(n):
for j, cur_k in enumerate(k):
np.testing.assert_allclose(function(cur_n, cur_k), results[i, j], err_msg=('n: ' + str(cur_n) + ', k: ' + str(cur_k)))
# Test first as scalar
for i, cur_n in enumerate(n):
np.testing.assert_allclose(function(cur_n, k_grid[:, i]), results[i, :], err_msg=('n: ' + str(cur_n) + ', k: ' + str(k_grid[:, i])))
# Test two as scalar
for j, cur_k in enumerate(k):
np.testing.assert_allclose(function(n_grid[j, :], cur_k), results[:, j], err_msg=('n: ' + str(n_grid[j, :]) + ', k: ' + str(cur_k)))
# Test matrix
np.testing.assert_allclose(function(n_grid, k_grid), results.T, err_msg=('Matrix computation failed!'))
def genSphCoords():
""" Generates cartesian (x,y,z) and spherical (theta, phi) coordinates of a sphere
Returns
-------
coords : named tuple
holds cartesian (x,y,z) and spherical (theta, phi) coordinates
"""
coords = namedtuple('coords', ['x', 'y', 'z', 'az', 'el'])
az = _np.linspace(0, 2 * _np.pi, 360)
el = _np.linspace(0, _np.pi, 181)
coords.x = _np.outer(_np.cos(az), _np.sin(el))
coords.y = _np.outer(_np.sin(az), _np.sin(el))
coords.z = _np.outer(_np.ones(360), _np.cos(el))
coords.el, coords.az = _np.meshgrid(_np.linspace(0, _np.pi, 181),
_np.linspace(0, 2 * _np.pi, 360))
return coords
def sph_harm_all(nMax, az, el, type='complex'):
'''Compute all sphercial harmonic coefficients up to degree nMax.
Parameters
----------
nMax : (int)
Maximum degree of coefficients to be returned. n >= 0
az: (float), array_like
Azimuthal (longitudinal) coordinate [0, 2pi], also called Theta.
el : (float), array_like
Elevation (colatitudinal) coordinate [0, pi], also called Phi.
Returns
-------
y_mn : (complex float), array_like
Complex spherical harmonics of degrees n [0 ... nMax] and all corresponding
orders m [-n ... n], sampled at [az, el]. dim1 corresponds to az/el pairs,
dim2 to oder/degree (m, n) pairs like 0/0, -1/1, 0/1, 1/1, -2/2, -1/2 ...
'''
m, n = mnArrays(nMax)
mA, azA = _np.meshgrid(m, az)
nA, elA = _np.meshgrid(n, el)
return sph_harm(mA, nA, azA, elA, type=type)
def draw2dsurface(X, Y, zf):
fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(X, Y)
Z = X*0
for i in range(len(X)):
for j in range(len(X[0])):
Z[i][j] = zf([X[i][j], Y[i][j]])
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
ax.set_zlim(np.min(Z.flatten()), np.max(Z.flatten()))
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
fig.colorbar(surf, shrink=0.5, aspect=5)
# plt.show()
def make_mask(self, shape):
try:
nrow, ncol = shape
except TypeError:
nrow = ncol = shape
# HACK: to make the masks consistent with ``rand_position()``
ix = self.xc - np.arange(ncol)
iy = self.yc - np.arange(nrow)
mx, my = np.meshgrid(ix, iy)
rho, phi = self.cart2pol(mx, my)
mask_rho = (rho >= self.rin) & (rho <= self.rout)
mask_phi = (phi >= self.abegin) & (phi <= self.aend)
if self.aend > 360:
mask_phi |= (phi <= (self.aend-360))
mask = mask_rho & mask_phi
return mask
def mapFunction( x , y , func , ax = None, arrayInput = False, n = 10, title = None, **kwargs ) :
"""
Plot function on a regular grid
x : 1d array
y : 1d array
func : function to map
arrayInput : False if func(x,y) , True if func( [x,y] )
"""
if ax is None :
fig , ax = plt.subplots()
X , Y = np.meshgrid( x , y )
if not arrayInput :
Z = func( X.flatten() , Y.flatten() ).reshape(X.shape)
else :
Z = func( np.stack( [ X.flatten() , Y.flatten() ]) )
ax.contourf( X , Y , Z , n , **kwargs)
if title is not None : ax.set_title(title)
return ax
def __init__(self, width, sample_spacing=0.025, samples=384):
''' Produces a Pinhole.
Args:
width (`float`): the width of the pinhole.
sample_spacing (`float`): spacing of samples in the synthetic image.
samples (`int`): number of samples per dimension in the synthetic image.
'''
self.width = width
# produce coordinate arrays
ext = samples / 2 * sample_spacing
x, y = np.linspace(-ext, ext, samples), np.linspace(-ext, ext, samples)
xv, yv = np.meshgrid(x, y)
w = width / 2
# paint a circle on a black background
arr = np.zeros((samples, samples))
arr[sqrt(xv**2 + yv**2) < w] = 1
super().__init__(data=arr, sample_spacing=sample_spacing, synthetic=True)
def defineSensorLoc(self,XYZ):
#############################################
# DEFINE TRANSMITTER AND RECEIVER LOCATIONS
#
# XYZ: N X 3 array containing transmitter center locations
# **NOTE** for this sensor, we know where the receivers are relative to transmitters
self.TxLoc = XYZ
dx,dy = np.meshgrid([-0.8,-0.4,0.,0.4,0.8],[-0.8,-0.4,0.,0.4,0.8])
dx = mkvc(dx)
dy = mkvc(dy)
N = np.shape(XYZ)[0]
X = np.kron(XYZ[:,0],np.ones((25))) + np.kron(np.ones((N)),dx)
Y = np.kron(XYZ[:,1],np.ones((25))) + np.kron(np.ones((N)),dy)
Z = np.kron(XYZ[:,2],np.ones((25)))
self.RxLoc = np.c_[X,Y,Z]
self.TxID = np.kron(np.arange(1,np.shape(XYZ)[0]+1),np.ones((25)))
self.RxComp = np.kron(3*np.ones(np.shape(XYZ)[0]),np.ones((25)))
def videoSeq(number_cells,inputIm,simtime,resolution,video_step):
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X = np.arange(number_cells)
Y = X
X, Y = np.meshgrid(X, Y)
Z = inputIm[:,:,0]
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='coolwarm',
linewidth=0, antialiased=False)
ax.set_title('Input stimulus. Time = 0.0 ms',y=1.08)
fig.show()
ax.axes.figure.canvas.draw()
for t in np.arange(0.0,int(simtime/resolution),video_step/resolution):
surf.remove()
Z = inputIm[:,:,t]
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='coolwarm',
linewidth=0, antialiased=False)
ax.set_title('Input stimulus. Time = %s ms'%str(t),y=1.08)
ax.axes.figure.canvas.draw()
# Estimation of Local Field Potential(LFP)
def build_random_variables(self, **kwargs):
# All this is done just once per batch (i.e. until `clear_random_variables` is called)
np.random.seed()
imshape = kwargs.get('imshape')
# Build and scale random fields
random_field_x = np.random.uniform(-1, 1, imshape) * self.alpha
random_field_y = np.random.uniform(-1, 1, imshape) * self.alpha
# Smooth random field (this has to be done just once per reset)
sdx = gaussian_filter(random_field_x, self.sigma, mode='reflect')
sdy = gaussian_filter(random_field_y, self.sigma, mode='reflect')
# Make meshgrid
x, y = np.meshgrid(np.arange(imshape[1]), np.arange(imshape[0]))
# Make inversion coefficient
_inverter = 1. if not self.invert else -1.
# Distort meshgrid indices (invert if required)
flow_y, flow_x = (y + _inverter * sdy).reshape(-1, 1), (x + _inverter * sdx).reshape(-1, 1)
# Set random states
self.set_random_variable('flow_x', flow_x)
self.set_random_variable('flow_y', flow_y)
def test_bowtie(self):
# Test the bowtie filter with sine waves
wDeg = 3
nPix = 257
sf = 1
orientation = 0
x = y = np.linspace(-wDeg // 2, wDeg // 2, nPix + 1)
u, v = np.meshgrid(x, y)
ramp = np.sin(orientation * np.pi / 180) * u
ramp -= np.cos(orientation * np.pi / 180) * v
img = np.sin(2 * np.pi * sf * ramp)
fimg = fft2(img)
orientation_widths = [1, 10, 20, 40, 80, 100]
for x in orientation_widths:
filt = OrientationFilter('bowtie', 90, x, nPix, .2, nPix + 1,
'triangle')
filt = filt.filter
filt = 1 - filt
filt = fftshift(filt)
out = ifft2(fimg * filt).real.astype(int)
self.assertEqual(np.sum(out), 0)
