The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
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 Example 2
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) Example 3
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() Example 4
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) Example 5
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) Example 6
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) Example 7
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 Example 8
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))) Example 9
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] Example 10
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]
# ???? Example 11
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] Example 12
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]
# ???? Example 13
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]
# ???? Example 14
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]
# ???? Example 15
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]
# ???? Example 16
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]
# ???? Example 17
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(' '))) Example 18
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 Example 19
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 Example 20
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!')) Example 21
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 Example 22
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) Example 23
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() Example 24
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 Example 25
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 Example 26
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) Example 27
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))) Example 28
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) Example 29
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) Example 30
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) Example 31
def test_plot_bounds_2d(self):
x = np.arange(1, 5)
y = np.arange(5, 10)
x2d, y2d = np.meshgrid(x, y)
x_bnds = np.arange(0.5, 4.51, 1.0)
y_bnds = np.arange(4.5, 9.51, 1.0)
# the borders are not modified
x_bnds[0] = 1.0
x_bnds[-1] = 4.0
y_bnds[0] = 5.0
y_bnds[-1] = 9.0
x2d_bnds, y2d_bnds = np.meshgrid(x_bnds, y_bnds)
d = psyd.CFDecoder()
# test x bounds
bounds = d.get_plotbounds(xr.Variable(('y', 'x'), x2d))
self.assertAlmostArrayEqual(bounds, x2d_bnds)
# test y bounds
bounds = d.get_plotbounds(xr.Variable(('y', 'x'), y2d))
self.assertAlmostArrayEqual(bounds, y2d_bnds) Example 32
def draw(X, pred, means, covariances, output):
xp = cupy.get_array_module(X)
for i in six.moves.range(2):
labels = X[pred == i]
if xp is cupy:
labels = labels.get()
plt.scatter(labels[:, 0], labels[:, 1], c=np.random.rand(3))
if xp is cupy:
means = means.get()
covariances = covariances.get()
plt.scatter(means[:, 0], means[:, 1], s=120, marker='s', facecolors='y',
edgecolors='k')
x = np.linspace(-5, 5, 1000)
y = np.linspace(-5, 5, 1000)
X, Y = np.meshgrid(x, y)
for i in six.moves.range(2):
Z = mlab.bivariate_normal(X, Y, np.sqrt(covariances[i][0]),
np.sqrt(covariances[i][1]),
means[i][0], means[i][1])
plt.contour(X, Y, Z)
plt.savefig(output) Example 33
def checkDiamondLattice(JJ, II, L, d, offset, CSmooth, doPlot = False):
arr = np.arange(0, L+1, 2*d)
narr = -arr
arr = narr.tolist()[::-1] + arr.tolist()[1::]
X, Y = np.meshgrid(arr, arr)
Q1 = np.zeros((X.size, 2))
Q1[:, 0] = Y.flatten()
Q1[:, 1] = X.flatten()
arr2 = np.arange(d, L+1, 2*d)
narr = -arr2
arr2 = narr.tolist()[::-1] + arr2.tolist()
X, Y = np.meshgrid(arr2, arr2)
Q2 = np.zeros((X.size, 2))
Q2[:, 0] = Y.flatten()
Q2[:, 1] = X.flatten()
Q = np.concatenate((Q1, Q2), 0)
return checkLattice(Q, JJ, II, L, d, offset, CSmooth, doPlot) Example 34
def make2ShakingCircles(NFrames = 200, T1 = 20, T2 = 20*np.pi, A1 = 10, A2 = 10, ydim = 50):
print("T1 = ", T1, ", T2 = ", T2)
IDims = (ydim, ydim*2, 3)
I = np.zeros((NFrames, ydim*ydim*2*3))
[X, Y] = np.meshgrid(np.arange(ydim*2), np.arange(ydim))
yc = float(ydim/2)
R = float(ydim/8)
for i in range(NFrames):
x1c = float(ydim/2) - A1*np.cos(2*np.pi*i/T1)
x2c = 3*float(ydim/2) - A2*np.cos(2*np.pi*i/T2)
f = np.zeros((X.shape[0], X.shape[1], 3))
for k in range(3):
f[:, :, k] = ((X-x1c)**2 + (Y-yc)**2 < R**2) + ((X-x2c)**2 + (Y-yc)**2 < R**2)
I[i, :] = f.flatten()
I = 0.5*(I + 1)
return (I, IDims)
#############################################################
#### OTHER VIDEO TOOLS #####
############################################################# Example 35
def doTotalOrderExperiment(N, binaryWeights = False):
I, J = np.meshgrid(np.arange(N), np.arange(N))
I = I[np.triu_indices(N, 1)]
J = J[np.triu_indices(N, 1)]
#[I, J] = [I[0:N-1], J[0:N-1]]
NEdges = len(I)
R = np.zeros((NEdges, 2))
R[:, 0] = J
R[:, 1] = I
#W = np.random.rand(NEdges)
W = np.ones(NEdges)
#Note: When using binary weights, Y is not necessarily a cocycle
Y = I - J
if binaryWeights:
Y = np.ones(NEdges)
(s, I, H) = doHodge(R, W, Y, verbose = True)
printConsistencyRatios(Y, I, H, W)
#Random flip experiment with linear order Example 36
def plot_cost_to_go_mountain_car(env, estimator, num_tiles=20):
x = np.linspace(env.observation_space.low[0], env.observation_space.high[0], num=num_tiles)
y = np.linspace(env.observation_space.low[1], env.observation_space.high[1], num=num_tiles)
X, Y = np.meshgrid(x, y)
Z = np.apply_along_axis(lambda _: -np.max(estimator.predict(_)), 2, np.dstack([X, Y]))
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1,
cmap=matplotlib.cm.coolwarm, vmin=-1.0, vmax=1.0)
ax.set_xlabel('Position')
ax.set_ylabel('Velocity')
ax.set_zlabel('Value')
ax.set_title("Mountain \"Cost To Go\" Function")
fig.colorbar(surf)
plt.show() Example 37
def test_megkde_2d_basic():
# Draw from normal, fit KDE, see if sampling from kde's pdf recovers norm
np.random.seed(1)
data = np.random.multivariate_normal([0, 1], [[1.0, 0.], [0., 0.75 ** 2]], size=10000)
xs, ys = np.linspace(-4, 4, 50), np.linspace(-4, 4, 50)
xx, yy = np.meshgrid(xs, ys, indexing='ij')
samps = np.vstack((xx.flatten(), yy.flatten())).T
zs = MegKDE(data).evaluate(samps).reshape(xx.shape)
zs_x = zs.sum(axis=1)
zs_y = zs.sum(axis=0)
cs_x = zs_x.cumsum()
cs_x /= cs_x[-1]
cs_x[0] = 0
cs_y = zs_y.cumsum()
cs_y /= cs_y[-1]
cs_y[0] = 0
samps_x = interp1d(cs_x, xs)(np.random.uniform(size=10000))
samps_y = interp1d(cs_y, ys)(np.random.uniform(size=10000))
mu_x, std_x = norm.fit(samps_x)
mu_y, std_y = norm.fit(samps_y)
assert np.isclose(mu_x, 0, atol=0.1)
assert np.isclose(std_x, 1.0, atol=0.1)
assert np.isclose(mu_y, 1, atol=0.1)
assert np.isclose(std_y, 0.75, atol=0.1) Example 38
def test_grid_data(self):
x, y = np.linspace(-3, 3, 200), np.linspace(-5, 5, 200)
xx, yy = np.meshgrid(x, y, indexing='ij')
xs, ys = xx.flatten(), yy.flatten()
chain = np.vstack((xs, ys)).T
pdf = (1 / (2 * np.pi)) * np.exp(-0.5 * (xs * xs + ys * ys / 4))
c = ChainConsumer()
c.add_chain(chain, parameters=['x', 'y'], weights=pdf, grid=True)
summary = c.analysis.get_summary()
x_sum = summary['x']
y_sum = summary['y']
expected_x = np.array([-1.0, 0.0, 1.0])
expected_y = np.array([-2.0, 0.0, 2.0])
threshold = 0.1
assert np.all(np.abs(expected_x - x_sum) < threshold)
assert np.all(np.abs(expected_y - y_sum) < threshold) Example 39
def surface_bounds(f, bounds, cmap=cm.magma_r, resolution=[100, 100], type='2D'):
"""
Plot a function over a grid
:param f: function to plot in Z
:param bounds:
:param cmap: cmap
:param resolution:
:param type: '2D' or '3D'
:return:
"""
resolution = np.array(resolution)
assert bounds.get_n_dim() == 2, 'Bounds are not 2D'
x = np.linspace(start=bounds.get_min(0), stop=bounds.get_max(0), num=resolution[0])
y = np.linspace(start=bounds.get_min(1), stop=bounds.get_max(1), num=resolution[1])
X, Y = np.meshgrid(x, y)
Z = f(np.vstack((X.flatten(), Y.flatten())).T).reshape(X.shape) # Evaluate grid on the function
surface(X=X, Y=Y, Z=Z, type=type, cmap=cmap)
# def surface_grid(f, x, y) Example 40
def visualize(config, vae):
if(config['n_z'] != 2):
print("Skipping visuals since n_z is not 2")
return
nx = ny = 20
x_values = np.linspace(-3, 3, nx)
y_values = np.linspace(-3, 3, ny)
canvas = np.empty((28*ny, 28*nx))
for i, yi in enumerate(x_values):
for j, xi in enumerate(y_values):
z_mu = np.array([[xi, yi]])
x_mean = vae.generate(np.tile(z_mu, [config['batch_size'], 1]))
canvas[(nx-i-1)*28:(nx-i)*28, j*28:(j+1)*28] = x_mean[0].reshape(28, 28)
plt.figure(figsize=(8, 10))
Xi, Yi = np.meshgrid(x_values, y_values)
plt.imshow(canvas, origin="upper")
plt.tight_layout()
img = "samples/2d-visualization.png"
plt.savefig(img)
hc.io.sample(config, [{"label": "2d visualization", "image": img}]) Example 41
def dihedral_transform_batch(x):
g = np.random.randint(low=0, high=8, size=x.shape[0])
h, w = x.shape[-2:]
hh = (h - 1) / 2.
hw = (w - 1) / 2.
I, J = np.meshgrid(np.linspace(-hh, hh, x.shape[-2]), np.linspace(-hw, hw, x.shape[-1]))
C = np.r_[[I, J]]
D4C = np.einsum('...ij,jkl->...ikl', D4, C)
D4C[:, 0] += hh
D4C[:, 1] += hw
D4C = D4C.astype(int)
x_out = np.empty_like(x)
for i in range(x.shape[0]):
I, J = D4C[g[i]]
x_out[i, :] = x[i][:, J, I]
return x_out Example 42
def contour_plot_of_rates(lab, r, labels=False, last_year=False):
plt.figure()
plt.title('{}'.format(r))
resoluion = 100
n = math.ceil(lab.time.size / resoluion)
if last_year:
k = n - int(1 / lab.dt)
else:
k = 1
z = lab.estimated_rates[r][:, k - 1:-1:n]
# lim = np.max(np.abs(z))
# lim = np.linspace(-lim - 0.1, +lim + 0.1, 51)
X, Y = np.meshgrid(lab.time[k::n], -lab.x)
plt.xlabel('Time')
CS = plt.contourf(X, Y, z, 20, cmap=ListedColormap(
sns.color_palette("Blues", 51)))
if labels:
plt.clabel(CS, inline=1, fontsize=10, colors='w')
cbar = plt.colorbar(CS)
plt.ylabel('Depth')
ax = plt.gca()
ax.ticklabel_format(useOffset=False)
cbar.ax.set_ylabel('Rate %s [M/V/T]' % r)
return ax Example 43
def contour_plot_of_delta(lab, element, labels=False, last_year=False):
plt.figure()
plt.title('Rate of %s consumption/production' % element)
resoluion = 100
n = math.ceil(lab.time.size / resoluion)
if last_year:
k = n - int(1 / lab.dt)
else:
k = 1
z = lab.species[element]['rates'][:, k - 1:-1:n]
lim = np.max(np.abs(z))
lim = np.linspace(-lim - 0.1, +lim + 0.1, 51)
X, Y = np.meshgrid(lab.time[k:-1:n], -lab.x)
plt.xlabel('Time')
CS = plt.contourf(X, Y, z, 20, cmap=ListedColormap(sns.color_palette(
"RdBu_r", 101)), origin='lower', levels=lim, extend='both')
if labels:
plt.clabel(CS, inline=1, fontsize=10, colors='w')
cbar = plt.colorbar(CS)
plt.ylabel('Depth')
ax = plt.gca()
ax.ticklabel_format(useOffset=False)
cbar.ax.set_ylabel('Rate of %s change $[\Delta/T]$' % element)
return ax Example 44
def _meshgrid(self):
with tf.variable_scope('_meshgrid'):
x_t = tf.matmul(tf.ones(shape=tf.stack([self.out_height, 1])),
tf.transpose(tf.expand_dims(tf.linspace(-1.0, 1.0, self.out_width), 1), [1, 0]))
y_t = tf.matmul(tf.expand_dims(tf.linspace(-1.0, 1.0, self.out_height), 1),
tf.ones(shape=tf.stack([1, self.out_width])))
x_t_flat = tf.reshape(x_t, (1, -1))
y_t_flat = tf.reshape(y_t, (1, -1))
px,py = tf.stack([x_t_flat],axis=2),tf.stack([y_t_flat],axis=2)
#source control points
x,y = tf.linspace(-1.,1.,self.Column_controlP_number),tf.linspace(-1.,1.,self.Row_controlP_number)
x,y = tf.meshgrid(x,y)
xs,ys = tf.transpose(tf.reshape(x,(-1,1))),tf.transpose(tf.reshape(y,(-1,1)))
cpx,cpy = tf.transpose(tf.stack([xs],axis=2),perm=[1,0,2]),tf.transpose(tf.stack([ys],axis=2),perm=[1,0,2])
px, cpx = tf.meshgrid(px,cpx);py, cpy = tf.meshgrid(py,cpy)
#Compute distance R
Rx,Ry = tf.square(tf.subtract(px,cpx)),tf.square(tf.subtract(py,cpy))
R = tf.add(Rx,Ry)
R = tf.multiply(R,tf.log(tf.clip_by_value(R,1e-10,1e+10)))
#Source coordinates
ones = tf.ones_like(x_t_flat)
grid = tf.concat([ones, x_t_flat, y_t_flat,R],0)
grid = tf.reshape(grid,[-1])
grid = tf.reshape(grid,[self.Column_controlP_number*self.Row_controlP_number+3,self.out_height*self.out_width])
return grid Example 45
def Affine_test(self,N,sizex,sizey,sizez,times,stop_time,typeofT,colors):
for i in range(times):
# Theta
idx = np.random.uniform(-1, 1);idy = np.random.uniform(-1, 1);idz = np.random.uniform(-1, 1)
swithx = np.random.uniform(0,1);swithy = np.random.uniform(0,1);swithz = np.random.uniform(0,1)
rotatex = np.random.uniform(-1, 1);rotatey = np.random.uniform(-1, 1);rotatez = np.random.uniform(-1, 1)
cx = np.array([idx,rotatey,rotatez,swithx]);cy = np.array([rotatex,idy,rotatez,swithy]);cz = np.array([rotatex,rotatey,idz,swithz])
# Source Grid
x = np.linspace(-sizex, sizex, N);y = np.linspace(-sizey, sizey, N);z = np.linspace(-sizez, sizez, N)
x, y, z = np.meshgrid(x, y, z)
xgs, ygs, zgs = x.flatten(), y.flatten(),z.flatten()
gps = np.vstack([xgs, ygs, zgs, np.ones_like(xgs)]).T
# transform
xgt = np.dot(gps, cx);ygt = np.dot(gps, cy);zgt = np.dot(gps, cz)
# display
showIm = ShowImage()
showIm.Show_transform(xgs,ygs,zgs,xgt,ygt,zgt,sizex,sizey,sizez,stop_time,typeofT,N,colors) Example 46
def TPS_test(self,N,sizex,sizey,sizez,control_num,show_times,stop_time,typeofT,colors):
for i in range(show_times):
# source control points
x, y, z = np.meshgrid(np.linspace(-1, 1, control_num),np.linspace(-1, 1, control_num), np.linspace(-1, 1, control_num))
xs = x.flatten();ys = y.flatten();zs = z.flatten()
cps = np.vstack([xs, ys, zs]).T
# target control points
xt = xs + np.random.uniform(-0.3, 0.3, size=xs.size);yt = ys + np.random.uniform(-0.3, 0.3, size=ys.size); zt = zs + np.random.uniform(-0.3, 0.3, size=zs.size)
# construct T
T = self.makeT(cps)
# solve cx, cy (coefficients for x and y)
xtAug = np.concatenate([xt, np.zeros(4)]);ytAug = np.concatenate([yt, np.zeros(4)]);ztAug = np.concatenate([zt, np.zeros(4)])
cx = nl.solve(T, xtAug); cy = nl.solve(T, ytAug); cz = nl.solve(T, ztAug)
# dense grid
x = np.linspace(-sizex, sizex, 2*sizex); y = np.linspace(-sizey, sizey, 2*sizey); z = np.linspace(-sizez, sizez, 2*sizez)
x, y, z = np.meshgrid(x, y, z)
xgs, ygs, zgs = x.flatten(), y.flatten(), z.flatten()
gps = np.vstack([xgs, ygs, zgs]).T
# transform
pgLift = self.liftPts(gps, cps) # [N x (K+3)]
xgt = np.dot(pgLift, cx.T);ygt = np.dot(pgLift, cy.T); zgt = np.dot(pgLift,cz.T)
# display
showIm = ShowImage()
showIm.Show_transform(xgs,ygs,zgs,xgt,ygt,zgt,sizex,sizey,sizez,stop_time,typeofT,N,colors) Example 47
def _initialize_grid(self, mmf):
self.nr, self.nlat, self.nlon\
= struct.unpack("3i", mmf.read(12))
self.dr, self.dlat, self.dlon\
= struct.unpack("3f", mmf.read(12))
self.r0, self.lat0, self.lon0\
= struct.unpack("3f", mmf.read(12))
self.mr = self.r0 + (self.nr - 1) * self.dr
self.mlat = self.lat0 + (self.nlat - 1) * self.dlat
self.mlon = self.lon0 + (self.nlon - 1) * self.dlon
self.dtheta = radians(self.dlat)
self.dphi = radians(self.dlon)
self.ntheta, self.nphi = self.nlat, self.nlon
self.theta0 = radians(90 - self.lat0)
self.phi0 = radians(self.lon0)
r = [self.r0 + self.dr * ir for ir in range(self.nr)]
theta = [radians(90. - self.lat0 - self.dlat * ilat)
for ilat in range(self.nlat)]
phi = [radians((self.lon0 + self.dlon * ilon) % 360.)
for ilon in range(self.nlon)]
R, T, P = np.meshgrid(r, theta, phi, indexing='ij')
self.nodes = {'r': R, 'theta': T, 'phi': P} Example 48
def mapinterpolation(data, x=None, y=None, interpx=1, interpy=1):
"""Interpolate a 2D map."""
# Get data size :
dim2, dim1 = data.shape
# Define xticklabel and yticklabel :
if x is None:
x = np.arange(0, dim1, interpx)
if y is None:
y = np.arange(0, dim2, interpy)
# Define the meshgrid :
Xi, Yi = np.meshgrid(
np.arange(0, dim1 - 1, interpx), np.arange(0, dim2 - 1, interpy))
# 2D interpolation :
datainterp = interp2(data, Xi, Yi)
# Linearly interpolate vectors :
xveci = np.linspace(x[0], x[-1], datainterp.shape[0])
yveci = np.linspace(y[0], y[-1], datainterp.shape[1])
return datainterp, xveci, yveci Example 49
def _shift_interp(array, shift_value, mode='constant', cval=0):
# Manual alternative to built-in function: slightly slower
Ndim = array.ndim
dims = array.shape
dtype = array.dtype.kind
if (Ndim == 1):
pass
elif (Ndim == 2):
x, y = np.meshgrid(np.arange(dims[1], dtype=dtype), np.arange(dims[0], dtype=dtype))
x -= shift_value[0]
y -= shift_value[1]
shifted = ndimage.map_coordinates(img, [y, x], mode=mode, cval=cval)
return shifted Example 50
def _rotate_interp(array, alpha, center, mode='constant', cval=0):
'''
Rotation around a provided center
This is the only way to be sure where exactly is the center of rotation.
'''
dtype = array.dtype
dims = array.shape
alpha_rad = -np.deg2rad(alpha)
x, y = np.meshgrid(np.arange(dims[1], dtype=dtype), np.arange(dims[0], dtype=dtype))
xp = (x-center[0])*np.cos(alpha_rad) + (y-center[1])*np.sin(alpha_rad) + center[0]
yp = -(x-center[0])*np.sin(alpha_rad) + (y-center[1])*np.cos(alpha_rad) + center[1]
rotated = ndimage.map_coordinates(img, [yp, xp], mode=mode, cval=cval, order=3)
return rotated 