def testHyCmplx(self):
sample = \
self.meep.get_pw_material_hy(self.idx, (0,0,0), cmplx=True)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_mu_inf(idx), self.meep.mu_inf)
else:
self.assertEqual(sample.get_mu_inf(idx), 0)
hy = ex = ez = np.zeros((3,3,3), complex)
dz = dx = dt = 1
n = 0
sample.update_all(hy, ex, ez, dz, dx, dt, n)
for idx in np.ndindex(3, 3, 3):
self.assertEqual(hy[idx], 0j)
python类ndindex()的实例源码
def testExCmplx(self):
sample = \
self.upml.get_pw_material_ex(self.idx, (0,0,0), cmplx=True)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_eps_inf(idx), self.upml.eps_inf)
else:
self.assertEqual(sample.get_eps_inf(idx), 0)
ex = hz = hy = np.zeros((3,3,3), complex)
dy = dz = dt = 1
n = 0
sample.update_all(ex, hz, hy, dy, dz, dt, n)
for idx in np.ndindex(3, 3, 3):
self.assertEqual(ex[idx], 0j)
def testEyCmplx(self):
sample = \
self.upml.get_pw_material_ey(self.idx, (0,0,0), cmplx=True)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_eps_inf(idx), self.upml.eps_inf)
else:
self.assertEqual(sample.get_eps_inf(idx), 0)
ey = hx = hz = np.zeros((3,3,3), complex)
dz = dx = dt = 1
n = 0
sample.update_all(ey, hx, hz, dz, dx, dt, n)
for idx in np.ndindex(3, 3, 3):
self.assertEqual(ey[idx], 0j)
def testEzCmplx(self):
sample = \
self.upml.get_pw_material_ez(self.idx, (0,0,0), cmplx=True)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_eps_inf(idx), self.upml.eps_inf)
else:
self.assertEqual(sample.get_eps_inf(idx), 0)
ez = hy = hx = np.zeros((3,3,3), complex)
dx = dy = dt = 1
n = 0
sample.update_all(ez, hy, hx, dx, dy, dt, n)
for idx in np.ndindex(3, 3, 3):
self.assertEqual(ez[idx], 0j)
def testHxCmplx(self):
sample = \
self.upml.get_pw_material_hx(self.idx, (0,0,0), cmplx=True)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_mu_inf(idx), self.upml.mu_inf)
else:
self.assertEqual(sample.get_mu_inf(idx), 0)
hx = ez = ey = np.zeros((3,3,3), complex)
dy = dz = dt = 1
n = 0
sample.update_all(hx, ez, ey, dy, dz, dt, n)
for idx in np.ndindex(3, 3, 3):
self.assertEqual(hx[idx], 0j)
def testHyCmplx(self):
sample = \
self.upml.get_pw_material_hy(self.idx, (0,0,0), cmplx=True)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_mu_inf(idx), self.upml.mu_inf)
else:
self.assertEqual(sample.get_mu_inf(idx), 0)
hy = ex = ez = np.zeros((3,3,3), complex)
dz = dx = dt = 1
n = 0
sample.update_all(hy, ex, ez, dz, dx, dt, n)
for idx in np.ndindex(3, 3, 3):
self.assertEqual(hy[idx], 0j)
def testExReal(self):
sample = \
self.const_real.get_pw_material_ex(self.idx, (0,0,0), cmplx=False)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_eps_inf(idx), self.const_real.eps_inf)
else:
self.assertEqual(sample.get_eps_inf(idx), 0)
ex = hz = hy = np.zeros((3,3,3))
dy = dz = dt = self.spc.dt
n = 0
sample.update_all(ex, hz, hy, dy, dz, dt, n)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(ex[idx], self.const_real.value)
else:
self.assertEqual(ex[idx], 0)
def testEyReal(self):
sample = \
self.const_real.get_pw_material_ey(self.idx, (0,0,0), cmplx=False)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_eps_inf(idx), self.const_real.eps_inf)
else:
self.assertEqual(sample.get_eps_inf(idx), 0)
ey = hx = hz = np.zeros((3,3,3))
dz = dx = dt = 1
n = 0
sample.update_all(ey, hx, hz, dz, dx, dt, n)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(ey[idx], self.const_real.value)
else:
self.assertEqual(ey[idx], 0)
def testEzReal(self):
sample = \
self.const_real.get_pw_material_ez(self.idx, (0,0,0), cmplx=False)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_eps_inf(idx), self.const_real.eps_inf)
else:
self.assertEqual(sample.get_eps_inf(idx), 0)
ez = hy = hx = np.zeros((3,3,3))
dx = dy = dt = 1
n = 0
sample.update_all(ez, hy, hx, dx, dy, dt, n)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(ez[idx], self.const_real.value)
else:
self.assertEqual(ez[idx], 0)
def testHxReal(self):
sample = \
self.const_real.get_pw_material_hx(self.idx, (0,0,0), cmplx=False)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_mu_inf(idx), self.const_real.mu_inf)
else:
self.assertEqual(sample.get_mu_inf(idx), 0)
hx = ez = ey = np.zeros((3,3,3))
dy = dz = dt = 1
n = 0
sample.update_all(hx, ez, ey, dy, dz, dt, n)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(hx[idx], self.const_real.value)
else:
self.assertEqual(hx[idx], 0)
def testHyReal(self):
sample = \
self.const_real.get_pw_material_hy(self.idx, (0,0,0), cmplx=False)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_mu_inf(idx), self.const_real.mu_inf)
else:
self.assertEqual(sample.get_mu_inf(idx), 0)
hy = ex = ez = np.zeros((3,3,3))
dz = dx = dt = 1
n = 0
sample.update_all(hy, ex, ez, dz, dx, dt, n)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(hy[idx], self.const_real.value)
else:
self.assertEqual(hy[idx], 0)
def testExCmplx(self):
sample = \
self.const_cmplx.get_pw_material_ex(self.idx, (0,0,0), cmplx=True)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(sample.get_eps_inf(idx), self.const_cmplx.eps_inf)
else:
self.assertEqual(sample.get_eps_inf(idx), 0)
ex = hz = hy = np.zeros((3,3,3), complex)
dy = dz = dt = 1
n = 0
sample.update_all(ex, hz, hy, dy, dz, dt, n)
for idx in np.ndindex(3, 3, 3):
if idx == self.idx:
self.assertEqual(ex[idx], self.const_cmplx.value)
else:
self.assertEqual(ex[idx], 0j)
def zero_pad(x, n):
"""
Return *x* zero padded to the dimensions specified in *n*.
"""
y = NP.zeros(n)
for index in NP.ndindex(x.shape):
y[index] = x[index]
return y
def correlate_frames(self, method='gaussian'):
"""Correlation of all grid points, creating a velocity field.
:param str method: Method of the peak finding algorithm
"""
for i, j in np.ndindex(self.grid_spec.get_grid_shape()):
window_a, window_b = self._get_window_frames(i, j)
displacement = (self._correlator.get_displacement(window_a, window_b,
subpixel_method=method))
self.u[i, j] += displacement[0]
self.v[i, j] += displacement[1]
return self.u, self.v
def distances_numba_array(cluster):
# Original: diff = cluster[:, np.newaxis, :] - cluster[np.newaxis, :, :]
# Since np.newaxis is not supported, we use reshape to do this
diff = (cluster.reshape(cluster.shape[0], 1, cluster.shape[1]) -
cluster.reshape(1, cluster.shape[0], cluster.shape[1]))
mat = (diff * diff)
# Original: mat = mat.sum(-1)
# Since axis argument is not supported, we write the loop out
out = np.empty(mat.shape[:2], dtype=mat.dtype)
for i in np.ndindex(out.shape):
out[i] = mat[i].sum()
return np.sqrt(out)
def slogdet(a):
"""Returns sign and logarithm of the determinat of an array.
It calculates the natural logarithm of the deteminant of a given value.
Args:
a (cupy.ndarray): The input matrix with dimension ``(..., N, N)``.
Returns:
tuple of :class:`~cupy.ndarray`:
It returns a tuple ``(sign, logdet)``. ``sign`` represents each
sign of the deteminant as a real number ``0``, ``1`` or ``-1``.
'logdet' represents the natural logarithm of the absolute of the
deteminant.
If the deteninant is zero, ``sign`` will be ``0`` and ``logdet``
will be ``-inf``.
The shapes of both ``sign`` and ``logdet`` are equal to
``a.shape[:-2]``.
.. seealso:: :func:`numpy.linalg.slogdet`
"""
if not cuda.cusolver_enabled:
raise RuntimeError('Current cupy only supports cusolver in CUDA 8.0')
if a.ndim < 2:
msg = ('%d-dimensional array given. '
'Array must be at least two-dimensional' % a.ndim)
raise linalg.LinAlgError(msg)
dtype = numpy.find_common_type((a.dtype.char, 'f'), ())
shape = a.shape[:-2]
sign = cupy.empty(shape, dtype)
logdet = cupy.empty(shape, dtype)
a = a.astype(dtype)
for index in numpy.ndindex(*shape):
s, l = _slogdet_one(a[index])
sign[index] = s
logdet[index] = l
return sign, logdet
def query_ball_point(self, x, r, p=2., eps=0):
"""Find all points within r of x
Parameters
==========
x : array_like, shape tuple + (self.m,)
The point or points to search for neighbors of
r : positive float
The radius of points to return
p : float 1<=p<=infinity
Which Minkowski p-norm to use
eps : nonnegative float
Approximate search. Branches of the tree are not explored
if their nearest points are further than r/(1+eps), and branches
are added in bulk if their furthest points are nearer than r*(1+eps).
Returns
=======
results : list or array of lists
If x is a single point, returns a list of the indices of the neighbors
of x. If x is an array of points, returns an object array of shape tuple
containing lists of neighbors.
Note: if you have many points whose neighbors you want to find, you may save
substantial amounts of time by putting them in a KDTree and using query_ball_tree.
"""
x = np.asarray(x)
if x.shape[-1]!=self.m:
raise ValueError("Searching for a %d-dimensional point in a %d-dimensional KDTree" % (x.shape[-1],self.m))
if len(x.shape)==1:
return self.__query_ball_point(x,r,p,eps)
else:
retshape = x.shape[:-1]
result = np.empty(retshape,dtype=np.object)
for c in np.ndindex(retshape):
result[c] = self.__query_ball_point(x[c], r, p=p, eps=eps)
return result
def __getitem__ (self, key):
import numpy as np
# Coerce key into a tuple of slice objects.
if not isinstance(key,tuple):
if hasattr(key,'__len__'): key = tuple(key)
else: key = (key,)
if len(key) == 1 and hasattr(key[0],'__len__'):
key = tuple(key[0])
if Ellipsis in key:
i = key.index(Ellipsis)
key = key[:i] + (slice(None),)*(self.ndim-len(key)+1) + key[i+1:]
key = key + (slice(None),)*(self.ndim-len(key))
if len(key) > self.ndim:
raise ValueError(("Too many dimensions for slicing."))
final_shape = tuple(len(np.arange(n)[k]) for n,k in zip(self.shape,key) if not isinstance(k,int))
data = np.ma.empty(final_shape, dtype=self.dtype)
outer_ndim = self._record_id.ndim
record_id = self._record_id.__getitem__(key[:outer_ndim])
# Final shape of each record (in case we did any reshaping along ni,nj,nk
# dimensions).
record_shape = self.shape[self._record_id.ndim:]
# Iterate of each record.
for ind in np.ndindex(record_id.shape):
r = int(record_id[ind])
if r >= 0:
data[ind] = self._buffer._fstluk(r)['d'].transpose().reshape(record_shape)[(Ellipsis,)+key[outer_ndim:]]
else:
data.mask = np.ma.getmaskarray(data)
data.mask[ind] = True
return data
def _calc(self, cars, nrb, ret_obj):
# Assume that an nD nrb should be averaged to be 1D
nrb = _mean_nd_to_1d(nrb)
shp = cars.shape[0:-2]
# Step row-by-row through image
for idx in _np.ndindex(shp):
if self.rng is None:
kkd = _kkrelation(bg=nrb + self.nrb_amp_offset,
cri=cars[idx] + self.cars_amp_offset,
phase_offset=self.phase_offset,
norm_by_bg=self.norm_to_nrb,
pad_factor=self.pad_factor)
else:
kkd = _kkrelation(bg=nrb[self.rng] + self.nrb_amp_offset,
cri=cars[idx][..., self.rng] + self.cars_amp_offset,
phase_offset=self.phase_offset,
norm_by_bg=self.norm_to_nrb,
pad_factor=self.pad_factor)
try:
ret_obj[idx] *= 0
if self.rng is None:
ret_obj[idx] += kkd
elif ret_obj[idx].size == kkd.size:
ret_obj[idx] += kkd
else:
ret_obj[idx][..., self.rng] += kkd
except:
return False
else:
pass
return True
def _calc(self, data, ret_obj, **kwargs):
self._inst_als = _AlsCvxopt(**kwargs)
try:
shp = data.shape[0:-2]
total_num = _np.array(shp).prod()
counter = 1
for idx in _np.ndindex(shp):
print('Detrended iteration {} / {}'.format(counter, total_num))
ph = _np.unwrap(_np.angle(data[idx]))
if self.rng is None:
err_phase = self._inst_als.calculate(ph)
else:
err_phase = self._inst_als.calculate(ph[..., self.rng])
h = _np.zeros(err_phase.shape)
h += _hilbert(err_phase)
correction_factor = 1/_np.exp(h) * _np.exp(-1j*err_phase)
if self.rng is None:
ret_obj[idx] *= correction_factor
else:
ret_obj[idx][..., self.rng] *= correction_factor
counter += 1
except:
return False
else:
# print(self._inst_als.__dict__)
return True
