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 inflection_points(points, axis, span):
'''
Find the list of vertices that preceed inflection points in a curve. The curve is differentiated
with respect to the coordinate system defined by axis and span.
axis: A vector representing the vertical axis of the coordinate system.
span: A vector representing the the horiztonal axis of the coordinate system.
returns: a list of points in space corresponding to the vertices that
immediately preceed inflection points in the curve
'''
coords_on_span = points.dot(span)
dx = np.gradient(coords_on_span)
coords_on_axis = points.dot(axis)
# Take the second order finite difference of the curve with respect to the
# defined coordinate system
finite_difference_2 = np.gradient(np.gradient(coords_on_axis, dx), dx)
# Compare the product of all neighboring pairs of points in the second derivative
# If a pair of points has a negative product, then the second derivative changes sign
# at one of those points, signalling an inflection point
is_inflection_point = [finite_difference_2[i] * finite_difference_2[i + 1] <= 0 for i in range(len(finite_difference_2) - 1)]
inflection_point_indices = [i for i, b in enumerate(is_inflection_point) if b]
if len(inflection_point_indices) == 0: # pylint: disable=len-as-condition
return []
return points[inflection_point_indices] Example 2
def get_mag_ang(img):
"""
Gets image gradient (magnitude) and orientation (angle)
Args:
img
Returns:
Gradient, orientation
"""
img = np.sqrt(img)
gx = cv2.Sobel(np.float32(img), cv2.CV_32F, 1, 0)
gy = cv2.Sobel(np.float32(img), cv2.CV_32F, 0, 1)
mag, ang = cv2.cartToPolar(gx, gy)
return mag, ang, gx, gy Example 3
def calculate_angle_quality_thread(self,
voi,
gantry,
couch,
calculate_from=0,
stepsize=1.0,
q=None,
avoid=[],
voi_cube=None,
gradient=True):
""" TODO: Documentation
"""
os.nice(1)
for couch_angle in couch:
qual = self.calculate_angle_quality(voi, gantry, couch_angle, calculate_from, stepsize, avoid, voi_cube,
gradient)
q.put({"couch": couch_angle, "gantry": gantry, "data": qual}) Example 4
def test_model_monovariate(func, var, par, k):
model = Model(func, var, par)
x, dx = np.linspace(0, 10, 100, retstep=True, endpoint=False)
U = np.cos(x * 2 * np.pi / 10)
fields = model.fields_template(x=x, U=U)
parameters = dict(periodic=True, k=k)
F = model.F(fields, parameters)
J_sparse = model.J(fields, parameters)
J_dense = model.J(fields, parameters, sparse=False)
J_approx = model.F.diff_approx(fields, parameters)
dxU = np.gradient(np.pad(U, 2, mode='wrap')) / dx
dxxU = np.gradient(dxU) / dx
dxU = dxU[2:-2]
dxxU = dxxU[2:-2]
assert np.isclose(F, k * dxxU, atol=1E-3).all()
assert np.isclose(J_approx, J_sparse.todense(), atol=1E-3).all()
assert np.isclose(J_approx, J_dense, atol=1E-3).all() Example 5
def test_model_bivariate(func, var, par):
model = Model(func, var, par)
x, dx = np.linspace(0, 10, 100, retstep=True, endpoint=False)
u = np.cos(x * 2 * np.pi / 10)
v = np.sin(x * 2 * np.pi / 10)
fields = model.fields_template(x=x, u=u, v=v)
parameters = dict(periodic=True, k1=1, k2=1)
F = model.F(fields, parameters)
J_sparse = model.J(fields, parameters)
J_dense = model.J(fields, parameters, sparse=False)
J_approx = model.F.diff_approx(fields, parameters)
dxu = np.gradient(np.pad(u, 2, mode='wrap')) / dx
dxxu = np.gradient(dxu) / dx
dxu = dxu[2:-2]
dxxu = dxxu[2:-2]
dxv = np.gradient(np.pad(v, 2, mode='wrap')) / dx
dxxv = np.gradient(dxv) / dx
dxv = dxv[2:-2]
dxxv = dxxv[2:-2]
assert np.isclose(F, np.vstack([dxxv, dxxu]).flatten('F'),
atol=1E-3).all()
assert np.isclose(J_approx, J_sparse.todense(), atol=1E-3).all()
assert np.isclose(J_approx, J_dense, atol=1E-3).all() Example 6
def calculate_sift_grid(self, image, gridH, gridW):
H, W = image.shape
Npatches = gridH.size
feaArr = np.zeros((Npatches, Nsamples * Nangles))
# calculate gradient
GH, GW = gen_dgauss(self.sigma)
IH = signal.convolve2d(image, GH, mode='same')
IW = signal.convolve2d(image, GW, mode='same')
Imag = np.sqrt(IH ** 2 + IW ** 2)
Itheta = np.arctan2(IH, IW)
Iorient = np.zeros((Nangles, H, W))
for i in range(Nangles):
Iorient[i] = Imag * np.maximum(np.cos(Itheta - angles[i]) ** alpha, 0)
for i in range(Npatches):
currFeature = np.zeros((Nangles, Nsamples))
for j in range(Nangles):
currFeature[j] = np.dot(self.weights, \
Iorient[j, gridH[i]:gridH[i] + self.pS, gridW[i]:gridW[i] + self.pS].flatten())
feaArr[i] = currFeature.flatten()
return feaArr Example 7
def extract_sift_patches(self, image, gridH, gridW):
# extracts the sift descriptor of patches
# in positions (gridH,gridW) in the image
H, W = image.shape
Npatches = gridH.size
feaArr = np.zeros((Npatches, Nsamples * Nangles))
# calculate gradient
GH, GW = gen_dgauss(self.sigma)
IH = signal.convolve2d(image, GH, mode='same')
IW = signal.convolve2d(image, GW, mode='same')
Imag = np.sqrt(IH ** 2 + IW ** 2)
Itheta = np.arctan2(IH, IW)
Iorient = np.zeros((Nangles, H, W))
for i in range(Nangles):
Iorient[i] = Imag * np.maximum(np.cos(Itheta - angles[i]) ** alpha, 0)
for i in range(Npatches):
currFeature = np.zeros((Nangles, Nsamples))
for j in range(Nangles):
currFeature[j] = np.dot(self.weights, \
Iorient[j, gridH[i]:gridH[i] + self.pS, gridW[i]:gridW[i] + self.pS].flatten())
feaArr[i] = currFeature.flatten()
# feaArr contains each descriptor in a row
feaArr = self.normalize_sift(feaArr)
return feaArr Example 8
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,]) Example 9
def contour_from_array_at_label_l(im_arr, l, thr=0.3, omit_axis=None, verbose=0):
if verbose > 0:
print(l)
array_label_l = im_arr == l
assert isinstance(array_label_l, np.ndarray)
gra = np.gradient(array_label_l.astype(np.bool).astype(np.float64))
if omit_axis is None:
norm_gra = np.sqrt(gra[0] ** 2 + gra[1] ** 2 + gra[2] ** 2) > thr
elif omit_axis == 'x':
norm_gra = np.sqrt(gra[1] ** 2 + gra[2] ** 2) > thr
elif omit_axis == 'y':
norm_gra = np.sqrt(gra[0] ** 2 + gra[2] ** 2) > thr
elif omit_axis == 'z':
norm_gra = np.sqrt(gra[0] ** 2 + gra[1] ** 2) > thr
else:
raise IOError
return norm_gra Example 10
def computeWeightsLocallyNormalized(I, centered_gradient=True, norm_radius=45):
h,w = I.shape[:2]
if centered_gradient:
gy,gx = np.gradient(I)[:2]
gysq = (gy**2).mean(axis=2) if gy.ndim > 2 else gy**2
gxsq = (gx**2).mean(axis=2) if gx.ndim > 2 else gx**2
gxsq_local_mean = cv2.blur(gxsq, ksize=(norm_radius, norm_radius))
gysq_local_mean = cv2.blur(gysq, ksize=(norm_radius, norm_radius))
w_horizontal = np.exp( - gxsq * 1.0/(2*np.maximum(1e-6, gxsq_local_mean)))
w_vertical = np.exp( - gysq * 1.0/(2*np.maximum(1e-6, gysq_local_mean)))
else:
raise Exception("NotImplementedYet")
return w_horizontal, w_vertical Example 11
def rel_vort(u,v,dx,dy):
#get the gradient of the wind in the u- and v-directions
du = np.gradient(u)
dv = np.gradient(v)
#compute the relative vorticity (units : 10^-5 s^-1)
vort = ((dv[-1]/dx) - (du[-2]/dy))*pow(10,5)
#return the vorticity (Units: 10^-5 s-1)
return vort
###############################################################################
########################## End function rel_vort() ##########################
###############################################################################
#################### 25.Begin function of wind_chill() ######################
## Required libraries: numpy #
## #
## Inputs: t2m = ndarray of 2-meter temperature values (Units: F) #
## #
## wind = ndarray of 10-meter bulk wind values (Units: MPH) #
## #
############################################################################### Example 12
def richardson(self):
u = self.dataSet.readNCVariable('U')
v = self.dataSet.readNCVariable('V')
u_corr = wrf.unstaggerX(u)
v_corr = wrf.unstaggerY(v)
wind = (u_corr*u_corr + v_corr*v_corr)**(0.5)
height = wrf.unstaggerZ(self.height)
#compute the vertical gradient of theta
dtheta = np.gradient(self.theta)[0]
du = np.gradient(wind)[0]
dz = np.gradient(height)[0]
#compute the richardson number
self.u10 = u_corr
self.v10 = v_corr
self.var = ((self.g/self.theta)*(dtheta/dz))/(du/dz)**(2)
self.var2 = self.var
self.varTitle = "Richardson Number\n" + self.dataSet.getTime()
#Set short variable title for time series
self.sTitle = "Richardson Number" Example 13
def smooth_abs(x, dx=0.01):
"""smoothed version of absolute vaue function, with quadratic instead of sharp bottom.
Derivative w.r.t. variable of interest. Width of quadratic can be controlled"""
x, n = _checkIfFloat(x)
y = np.abs(x)
idx = np.logical_and(x > -dx, x < dx)
y[idx] = x[idx]**2/(2.0*dx) + dx/2.0
# gradient
dydx = np.ones_like(x)
dydx[x <= -dx] = -1.0
dydx[idx] = x[idx]/dx
if n == 1:
y = y[0]
dydx = dydx[0]
return y, dydx Example 14
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,]) Example 15
def main():
"""Load image, compute derivatives, plot."""
img = data.camera()
imgy_np, imgx_np = np.gradient(img)
imgx_ski, imgy_ski = _compute_derivatives(img)
# dx
util.plot_images_grayscale(
[img, dx_symmetric(img), dx_forward(img), dx_backward(img), imgx_np, imgx_ski],
["Image", "dx (symmetric)", "dx (forward)", "dx (backward)",
"Ground Truth (numpy)", "Ground Truth (scikit-image)"]
)
# dy
util.plot_images_grayscale(
[img, dy_symmetric(img), dy_forward(img), dy_backward(img), imgy_np, imgy_ski],
["Image", "dy (symmetric)", "dy (forward)", "dy (backward)",
"Ground Truth (numpy)", "Ground Truth (scikit-image)"]
) Example 16
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,]) Example 17
def direction_map(dmap: DistanceMap):
r"""Computes normalized gradient of distance map. Not defined when length of
the gradient is zero.
.. math::
\hat{\mathbf{e}}_{S} = -\frac{\nabla S(\mathbf{x})}{\| \nabla S(\mathbf{x}) \|}
Args:
dmap (numpy.ndarray):
Distance map.
Returns:
DirectionMap: Direction map.
"""
u, v = np.gradient(dmap)
l = np.hypot(u, v)
# Avoids zero division
l[l == 0] = np.nan
# Flip order from (row, col) to (x, y)
return v / l, u / l
# Potentials Example 18
def __entropy(self, spectrum, m):
"""Calculates get m-th derivative of the real part of spectrum and
returns entropy of its absolute value. """
assert type(m) is int, 'm should be a (possitive) integer'
assert m > 0, 'need m > 0'
#get the real part of the spectrum
spect = np.array(spectrum)
spect = spect.real
# calculate the m-th derivative of the real part of the spectrum
spectrumDerivative = spect
for i in range(m):
spectrumDerivative = np.gradient(spectrumDerivative)
# now get the entropy of the abslolute value of the m-th derivative:
entropy = sp.stats.entropy(np.abs(spectrumDerivative))
return entropy Example 19
def threshold_gradients(self, grad_thresh):
"""Creates a new DepthImage by zeroing out all depths
where the magnitude of the gradient at that point is
greater than grad_thresh.
Parameters
----------
grad_thresh : float
A threshold for the gradient magnitude.
Returns
-------
:obj:`DepthImage`
A new DepthImage created from the thresholding operation.
"""
data = np.copy(self._data)
gx, gy = self.gradients()
gradients = np.zeros([gx.shape[0], gx.shape[1], 2])
gradients[:, :, 0] = gx
gradients[:, :, 1] = gy
gradient_mags = np.linalg.norm(gradients, axis=2)
ind = np.where(gradient_mags > grad_thresh)
data[ind[0], ind[1]] = 0.0
return DepthImage(data, self._frame) Example 20
def normal_cloud_im(self):
"""Generate a NormalCloudImage from the PointCloudImage.
Returns
-------
:obj:`NormalCloudImage`
The corresponding NormalCloudImage.
"""
gx, gy, _ = np.gradient(self.data)
gx_data = gx.reshape(self.height * self.width, 3)
gy_data = gy.reshape(self.height * self.width, 3)
pc_grads = np.cross(gx_data, gy_data) # default to point toward camera
pc_grad_norms = np.linalg.norm(pc_grads, axis=1)
normal_data = pc_grads / np.tile(pc_grad_norms[:, np.newaxis], [1, 3])
normal_im_data = normal_data.reshape(self.height, self.width, 3)
return NormalCloudImage(normal_im_data, frame=self.frame) Example 21
def visualize_derivatives(image):
'''
Plot gradient on left and Laplacian on right.
Only tested on 2D 1-channel float imags
'''
dx1,dy1 = np.gradient(image)
gradient = dx1 + 1j*dy1
a1 = np.abs(gradient)
plt.figure(None,(12,6))
plt.subplot(121)
a1 = mean_center(blur(exposure.equalize_hist(unitize(a1)),1))
plt.imshow(a1,
origin='lower',interpolation='nearest',cmap='gray',extent=(0,64,)*2)
plt.title('Gradient Magnitude')
plt.subplot(122)
laplacian = scipy.ndimage.filters.laplace(image)
lhist = mean_center(
blur(exposure.equalize_hist(unitize(laplacian)),1))
plt.imshow(lhist,
origin='lower',interpolation='nearest',cmap='gray',extent=(0,64,)*2)
plt.title('Laplacian')
return gradient, laplacian Example 22
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,]) Example 23
def _get_displacement_function(self, f):
"""
Getter function for calculating the displacement.
:param f: field that is used for the displacement, mainly velocity components
:returns: function of the Taylor expanded field to first order
"""
dx = self._distance
f_x, f_y = np.gradient(f , dx)
f_xx, f_xy = np.gradient(f_x, dx)
f_yx, f_yy = np.gradient(f_y, dx)
return lambda i, j, x, y : (f[i, j] + x*f_x[i, j] + y*f_y[i, j]
+ 0.5*(f_xx[i, j]*x**2 + 2*f_xy[i, j]*x*y + f_yy[i, j]*y**2))
#For the bilinear method the build in scipy method `map_coordinates <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.map_coordinates.html>`_ is used with *order* set to 1. Example 24
def psf_shift(lenslet_efl, dx, dy, mag=1):
''' Computes the shift of a PSF, in microns.
Args:
lenslet_efl (`microns`): EFL of lenslets.
dx (`np.ndarray`): dx gradient of wavefront.
dy (`np.ndarray`): dy gradient of wavefront.
mag (`float`): magnification of the collimation system.
Returns:
`numpy.ndarray`: array of PSF shifts.
Notes:
see eq. 12 of Chanan, "Principles of Wavefront Sensing and Reconstruction"
delta = m * fl * grad(z)
m is magnification of SH system
fl is lenslet focal length
grad(z) is the x, y gradient of the opd, z, which is expressed in
radians.
'''
coef = -mag * lenslet_efl
return coef * dx, coef * dy Example 25
def test_basic(self):
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x), dx)
assert_array_equal(gradient(v), dx) Example 26
def test_badargs(self):
# for 2D array, gradient can take 0, 1, or 2 extra args
x = np.array([[1, 1], [3, 4]])
assert_raises(SyntaxError, gradient, x, np.array([1., 1.]),
np.array([1., 1.]), np.array([1., 1.])) Example 27
def test_masked(self):
# Make sure that gradient supports subclasses like masked arrays
x = np.ma.array([[1, 1], [3, 4]],
mask=[[False, False], [False, False]])
out = gradient(x)[0]
assert_equal(type(out), type(x))
# And make sure that the output and input don't have aliased mask
# arrays
assert_(x.mask is not out.mask)
# Also check that edge_order=2 doesn't alter the original mask
x2 = np.ma.arange(5)
x2[2] = np.ma.masked
np.gradient(x2, edge_order=2)
assert_array_equal(x2.mask, [False, False, True, False, False]) Example 28
def test_datetime64(self):
# Make sure gradient() can handle special types like datetime64
x = np.array(
['1910-08-16', '1910-08-11', '1910-08-10', '1910-08-12',
'1910-10-12', '1910-12-12', '1912-12-12'],
dtype='datetime64[D]')
dx = np.array(
[-5, -3, 0, 31, 61, 396, 731],
dtype='timedelta64[D]')
assert_array_equal(gradient(x), dx)
assert_(dx.dtype == np.dtype('timedelta64[D]')) Example 29
def test_timedelta64(self):
# Make sure gradient() can handle special types like timedelta64
x = np.array(
[-5, -3, 10, 12, 61, 321, 300],
dtype='timedelta64[D]')
dx = np.array(
[2, 7, 7, 25, 154, 119, -21],
dtype='timedelta64[D]')
assert_array_equal(gradient(x), dx)
assert_(dx.dtype == np.dtype('timedelta64[D]')) Example 30
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,]) Example 31
def calc_v_vec(self, tp):
v_vec_all_bands = []
for ib in range(self.num_bands[tp]):
v_vec_k_ordered = self.velocity_signed[tp][ib][self.pos_idx[tp]]
v_vec_all_bands.append(self.grid_from_ordered_list(v_vec_k_ordered, tp, none_missing=True))
return np.array(v_vec_all_bands)
# def calc_v_vec(self, tp):
# # TODO: Take into account the fact that this gradient is found in three directions specified by the lattice, not
# # the x, y, and z directions. It must be corrected to account for this.
# energy_grid = self.array_from_kgrid('energy', tp)
# # print('energy:')
# # np.set_printoptions(precision=3)
# # print(energy_grid[0,:,:,:,0])
# N = self.kgrid_array[tp].shape
# k_grid = self.kgrid_array[tp]
# v_vec_result = []
# for ib in range(self.num_bands[tp]):
# v_vec = np.gradient(energy_grid[ib][:,:,:,0], k_grid[:,0,0,0] * self._rec_lattice.a, k_grid[0,:,0,1] * self._rec_lattice.b, k_grid[0,0,:,2] * self._rec_lattice.c)
# v_vec_rearranged = np.zeros((N[0], N[1], N[2], 3))
# for i in range(N[0]):
# for j in range(N[1]):
# for k in range(N[2]):
# v_vec_rearranged[i,j,k,:] = np.array([v_vec[0][i,j,k], v_vec[1][i,j,k], v_vec[2][i,j,k]])
# v_vec_rearranged *= A_to_m * m_to_cm / hbar
# v_vec_result.append(v_vec_rearranged)
# return np.array(v_vec_result)
# turns a kgrid property into a list of grid arrays of that property for k integration Example 32
def compute_beta(TT, minT):
"""
This function computes the volumetric thermal expansion as a numerical
derivative of the volume as a function of temperature V(T) given in the
input array *minT*. This array can obtained
from the free energy minimization which should be done before.
"""
grad=np.gradient(np.array(minT)) # numerical derivatives with numpy
betaT = np.array(grad) # grad contains the derivatives with respect to T
# also convert to np.array format
return betaT/minT Example 33
def compute_alpha(minT,ibrav):
"""
This function calculates the thermal expansion alphaT at different temperatures
from the input minT matrix by computing the numerical derivatives with numpy.
The input matrix minT has shape nT*6, where the first index is the temperature
and the second the lattice parameter. For example, minT[i,0] and minT[i,2] are
the lattice parameters a and c at the temperature i.
More ibrav types must be implemented
"""
grad=np.gradient(np.array(minT)) # numerical derivatives with numpy
alphaT = np.array(grad[0]) # grad[0] contains the derivatives with respect to T, which is the first axis in minT
# also convert to np.array format
# now normalize the alpha properly. It must be different for different ibrav
# to avoid a divide by 0 error (minT is zero for lattice parameters not defined
# in the system)
if ibrav==4:
alphaT[:,0] = alphaT[:,0]/minT[:,0]
alphaT[:,2] = alphaT[:,2]/minT[:,2]
return alphaT
################################################################################ Example 34
def calc_eud(dvh, a):
v = -np.gradient(dvh)
dose_bins = np.linspace(1, np.size(dvh), np.size(dvh))
dose_bins = np.round(dose_bins, 3)
bin_centers = dose_bins - 0.5
eud = np.power(np.sum(np.multiply(v, np.power(bin_centers, a))), 1 / float(a))
eud = np.round(eud, 2) * 0.01
return eud Example 35
def _compute_gradients(self):
"""Computes the gradients of the SDF.
Returns
-------
:obj:`list` of :obj:`numpy.ndarray` of float
A list of ndarrays of the same dimension as the SDF. The arrays
are in axis order and specify the gradients for that axis
at each point.
"""
self.gradients_ = np.gradient(self.data_) Example 36
def curvature(self, coords, delta=0.001):
"""
Returns an approximation to the local SDF curvature (Hessian) at the
given coordinate in grid basis.
Parameters
---------
coords : numpy 3-vector
the grid coordinates at which to get the curvature
Returns
-------
curvature : 3x3 ndarray of the curvature at the surface points
"""
# perturb local coords
coords_x_up = coords + np.array([delta, 0, 0])
coords_x_down = coords + np.array([-delta, 0, 0])
coords_y_up = coords + np.array([0, delta, 0])
coords_y_down = coords + np.array([0, -delta, 0])
coords_z_up = coords + np.array([0, 0, delta])
coords_z_down = coords + np.array([0, 0, -delta])
# get gradient
grad_x_up = self.gradient(coords_x_up)
grad_x_down = self.gradient(coords_x_down)
grad_y_up = self.gradient(coords_y_up)
grad_y_down = self.gradient(coords_y_down)
grad_z_up = self.gradient(coords_z_up)
grad_z_down = self.gradient(coords_z_down)
# finite differences
curvature_x = (grad_x_up - grad_x_down) / (4 * delta)
curvature_y = (grad_y_up - grad_y_down) / (4 * delta)
curvature_z = (grad_z_up - grad_z_down) / (4 * delta)
curvature = np.c_[curvature_x, np.c_[curvature_y, curvature_z]]
curvature = curvature + curvature.T
return curvature Example 37
def call_tad_borders(ii_results, cutoff=0):
"""
Calls TAD borders using the first derivative of the insulation index.
:param ii_results: (raw) insulation index results, numpy vector
:param cutoff: raw insulation index threshold for "true" TAD boundaries
"""
tad_borders = []
g = np.gradient(ii_results)
for i in range(len(ii_results)):
border_type = _border_type(i, g)
if border_type == 1 and ii_results[i] <= cutoff:
tad_borders.append(i)
return tad_borders Example 38
def orientation_field(bmp, sigma=3):
# Author: Shaun Harker, 2016
# Based on algorithm by Bazen and Gerez from "Systematic methods for the
# computation of the directional fields and singular points of fingerprints," 2002.
"""
Computes orientation field (result everywhere between -pi/2 and pi/2)
from the given vector field.
"""
u = bmp.astype(float)
du = np.gradient(u)
[ux, uy] = du
Y = scipy.ndimage.filters.gaussian_filter(2.0*ux*uy, sigma=sigma)
X = scipy.ndimage.filters.gaussian_filter(ux**2.0 - uy**2.0, sigma=sigma)
return .5 * np.arctan2(Y, X) Example 39
def calculate_quality_grid(self, voi, gantry, couch, calculate_from=0, stepsize=1.0, avoid=[], gradient=True):
""" TODO: Documentation
"""
result = self.calculate_quality_list(voi, gantry, couch, calculate_from, stepsize,
avoid=avoid, gradient=gradient)
result = sorted(result, key=cmp_to_key(cmp_sort))
grid_data = []
for x in result:
grid_data.append(x["data"][0])
result = np.reshape(grid_data, (len(gantry), len(couch)))
return result Example 40
def calculate_quality_list(self, voi, gantry, couch, calculate_from=0, stepsize=1.0, avoid=[], gradient=True):
""" TODO: Documentation
"""
q = Queue(32767)
process = []
d = voi.get_voi_cube()
d.cube = np.array(d.cube, dtype=np.float32)
voi_cube = DensityProjections(d)
result = []
for gantry_angle in gantry:
p = Process(
target=self.calculate_angle_quality_thread,
args=(voi, gantry_angle, couch, calculate_from, stepsize, q, avoid, voi_cube, gradient))
p.start()
p.deamon = True
process.append(p)
if len(process) > 2:
tmp = q.get()
result.append(tmp)
for p in process:
if not p.is_alive():
process.remove(p)
while not len(result) == len(gantry) * len(couch):
tmp = q.get()
result.append(tmp)
return result Example 41
def calculate_angle_quality(self,
voi,
gantry,
couch,
calculate_from=0,
stepsize=1.0,
avoid=[],
voi_cube=None,
gradient=True):
"""
Calculates a quality index for a given gantry/couch combination.
"""
voi_min, voi_max = voi.get_min_max()
for v in avoid:
v_min, v_max = v.get_min_max()
if voi_cube is None:
d = voi.get_voi_cube()
d.cube = np.array(d.cube, dtype=np.float32)
voi_cube = DensityProjections(d)
data, start, basis = self.calculate_projection(voi, gantry, couch, calculate_from, stepsize)
voi_proj, t1, t2 = voi_cube.calculate_projection(voi, gantry, couch, 1, stepsize)
if gradient:
gradient = np.gradient(data)
data = (gradient[0]**2 + gradient[1]**2)**0.5
a = data * (voi_proj > 0.0)
quality = sum(a)
area = sum(voi_proj > 0.0)
# ~ area = sum(data>0.0)/10
if gradient:
mean_quality = 10 - abs(quality / area)
else:
mean_quality = abs(quality / area)
return mean_quality, quality, area Example 42
def get_direction(asa, *, sigma=None):
"""Return epochs during which an animal was running left to right, or right
to left.
Parameters
----------
asa : AnalogSignalArray 1D
AnalogSignalArray containing the 1D position data.
sigma : float, optional
Smoothing to apply to position (x) before computing gradient estimate.
Default is 0.
Returns
-------
l2r, r2l : EpochArrays
EpochArrays corresponding to left-to-right and right-to-left movement.
"""
if sigma is None:
sigma = 0
if not isinstance(asa, core.AnalogSignalArray):
raise TypeError('AnalogSignalArray expected!')
assert asa.n_signals == 1, "1D AnalogSignalArray expected!"
direction = dxdt_AnalogSignalArray(asa.smooth(sigma=sigma),
rectify=False).ydata
direction[direction>=0] = 1
direction[direction<0] = -1
direction = direction.squeeze()
l2r = get_contiguous_segments(np.argwhere(direction>0).squeeze(), step=1)
l2r[:,1] -= 1 # change bounds from [inclusive, exclusive] to [inclusive, inclusive]
l2r = core.EpochArray(asa.time[l2r])
r2l = get_contiguous_segments(np.argwhere(direction<0).squeeze(), step=1)
r2l[:,1] -= 1 # change bounds from [inclusive, exclusive] to [inclusive, inclusive]
r2l = core.EpochArray(asa.time[r2l])
return l2r, r2l Example 43
def hessian(array):
(dy, dx) = np.gradient(array)
(dydy, dxdy) = np.gradient(dy)
(dydx, dxdx) = np.gradient(dx)
return np.dstack((dxdx, dydx, dxdy, dydy)) Example 44
def gen_dgauss(sigma):
'''
gradient of the gaussian on both directions.
'''
fwid = np.int(2 * np.ceil(sigma))
G = np.array(range(-fwid, fwid + 1)) ** 2
G = G.reshape((G.size, 1)) + G
G = np.exp(- G / 2.0 / sigma / sigma)
G /= np.sum(G)
GH, GW = np.gradient(G)
GH *= 2.0 / np.sum(np.abs(GH))
GW *= 2.0 / np.sum(np.abs(GW))
return GH, GW Example 45
def gradient_pdf(self, x):
"""Return the gradient of density of the joint prior at x."""
raise NotImplementedError Example 46
def gradient_logpdf(self, x, stepsize=None):
"""Return the gradient of log density of the joint prior at x.
Parameters
----------
x : float or np.ndarray
stepsize : float or list
Stepsize or stepsizes for the dimensions
"""
x = np.asanyarray(x)
ndim = x.ndim
x = x.reshape((-1, self.dim))
grads = np.zeros_like(x)
for i in range(len(grads)):
xi = x[i]
grads[i] = numgrad(self.logpdf, xi, h=stepsize)
grads[np.isinf(grads)] = 0
grads[np.isnan(grads)] = 0
if ndim == 0 or (ndim == 1 and self.dim > 1):
grads = grads[0]
return grads Example 47
def test_basic(self):
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x), dx)
assert_array_equal(gradient(v), dx) Example 48
def test_badargs(self):
# for 2D array, gradient can take 0, 1, or 2 extra args
x = np.array([[1, 1], [3, 4]])
assert_raises(SyntaxError, gradient, x, np.array([1., 1.]),
np.array([1., 1.]), np.array([1., 1.])) Example 49
def test_masked(self):
# Make sure that gradient supports subclasses like masked arrays
x = np.ma.array([[1, 1], [3, 4]],
mask=[[False, False], [False, False]])
out = gradient(x)[0]
assert_equal(type(out), type(x))
# And make sure that the output and input don't have aliased mask
# arrays
assert_(x.mask is not out.mask)
# Also check that edge_order=2 doesn't alter the original mask
x2 = np.ma.arange(5)
x2[2] = np.ma.masked
np.gradient(x2, edge_order=2)
assert_array_equal(x2.mask, [False, False, True, False, False]) Example 50
def test_datetime64(self):
# Make sure gradient() can handle special types like datetime64
x = np.array(
['1910-08-16', '1910-08-11', '1910-08-10', '1910-08-12',
'1910-10-12', '1910-12-12', '1912-12-12'],
dtype='datetime64[D]')
dx = np.array(
[-5, -3, 0, 31, 61, 396, 731],
dtype='timedelta64[D]')
assert_array_equal(gradient(x), dx)
assert_(dx.dtype == np.dtype('timedelta64[D]')) 