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
python类gradient()的实例源码
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
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
################################################################################
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
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_)
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
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
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)
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
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
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
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
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))
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
def gradient_pdf(self, x):
"""Return the gradient of density of the joint prior at x."""
raise NotImplementedError
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
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)
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.]))
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])
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]'))
