def run_test_matmul_aa_correlator_kernel(self, ntime, nstand, nchan, misalign=0):
x_shape = (ntime, nchan, nstand*2)
perm = [1,0,2]
x8 = ((np.random.random(size=x_shape+(2,))*2-1)*127).astype(np.int8)
x = x8.astype(np.float32).view(np.complex64).reshape(x_shape)
x = x.transpose(perm)
x = x[..., misalign:]
b_gold = np.matmul(H(x), x)
triu = np.triu_indices(x.shape[-1], 1)
b_gold[..., triu[0], triu[1]] = 0
x = x8.view(bf.DataType.ci8).reshape(x_shape)
x = bf.asarray(x, space='cuda')
x = x.transpose(perm)
x = x[..., misalign:]
b = bf.zeros_like(b_gold, space='cuda')
self.linalg.matmul(1, None, x, 0, b)
b = b.copy('system')
np.testing.assert_allclose(b, b_gold, RTOL*10, ATOL)
python类complex64()的实例源码
def flip_code(code):
if isinstance(code, (numpy.dtype,type)):
# since several things map to complex64 we must carefully select
# the opposite that is an exact match (ticket 1518)
if code == numpy.int8:
return gdalconst.GDT_Byte
if code == numpy.complex64:
return gdalconst.GDT_CFloat32
for key, value in codes.items():
if value == code:
return key
return None
else:
try:
return codes[code]
except KeyError:
return None
def solver(self,gy, solver=None, *args, **kwargs):
"""
The solver of NUFFT
:param gy: data, reikna array, (M,) size
:param solver: could be 'cg', 'L1TVOLS', 'L1TVLAD'
:param maxiter: the number of iterations
:type gy: reikna array, dtype = numpy.complex64
:type solver: string
:type maxiter: int
:return: reikna array with size Nd
"""
import src._solver.solver_hsa
try:
return src._solver.solver_hsa.solver(self, gy, solver, *args, **kwargs)
except:
if numpy.ndarray==type(gy):
print("input gy must be a reikna array with dtype = numpy.complex64")
raise TypeError
else:
print("wrong")
raise TypeError
def forward(self, gx):
"""
Forward NUFFT on the heterogeneous device
:param gx: The input gpu array, with size=Nd
:type: reikna gpu array with dtype =numpy.complex64
:return: gy: The output gpu array, with size=(M,)
:rtype: reikna gpu array with dtype =numpy.complex64
"""
self.x_Nd = self.thr.copy_array(gx)
self._x2xx()
self._xx2k()
self._k2y()
gy = self.thr.copy_array(self.y)
return gy
def adjoint(self, gy):
"""
Adjoint NUFFT on the heterogeneous device
:param gy: The output gpu array, with size=(M,)
:type: reikna gpu array with dtype =numpy.complex64
:return: gx: The input gpu array, with size=Nd
:rtype: reikna gpu array with dtype =numpy.complex64
"""
self.y = self.thr.copy_array(gy)
self._y2k()
self._k2xx()
self._xx2x()
gx = self.thr.copy_array(self.x_Nd)
return gx
def selfadjoint(self, gx):
"""
selfadjoint NUFFT (Teplitz) on the heterogeneous device
:param gx: The input gpu array, with size=Nd
:type: reikna gpu array with dtype =numpy.complex64
:return: gx: The input gpu array, with size=Nd
:rtype: reikna gpu array with dtype =numpy.complex64
"""
self.x_Nd = self.thr.copy_array(gx)
self._x2xx()
self._xx2k()
self._k2y2k()
self._k2xx()
self._xx2x()
gx2 = self.thr.copy_array(self.x_Nd)
return gx2
def __init__(self):
"""
Constructor.
:param None:
:type None: Python NoneType
:return: NUFFT: the pynufft_hsa.NUFFT instance
:rtype: NUFFT: the pynufft_hsa.NUFFT class
:Example:
>>> import pynufft.pynufft
>>> NufftObj = pynufft.pynufft.NUFFT_cpu()
.. note:: requires plan()
.. seealso:: :method:`plan()'
.. todo:: test 3D case
"""
self.dtype=numpy.complex64
pass
def solve(self,gy, solver=None, *args, **kwargs):
"""
The solver of NUFFT_hsa
:param gy: data, reikna array, (M,) size
:param solver: could be 'cg', 'L1TVOLS', 'L1TVLAD'
:param maxiter: the number of iterations
:type gy: reikna array, dtype = numpy.complex64
:type solver: string
:type maxiter: int
:return: reikna array with size Nd
"""
from .._nonlin.solve_hsa import solve
try:
return solve(self, gy, solver, *args, **kwargs)
except:
if numpy.ndarray==type(gy):
print("input gy must be a reikna array with dtype = numpy.complex64")
raise TypeError
else:
print("wrong")
raise TypeError
def forward(self, gx):
"""
Forward NUFFT on the heterogeneous device
:param gx: The input gpu array, with size=Nd
:type: reikna gpu array with dtype =numpy.complex64
:return: gy: The output gpu array, with size=(M,)
:rtype: reikna gpu array with dtype =numpy.complex64
"""
self.x_Nd = self.thr.copy_array(gx)
self._x2xx()
self._xx2k()
self._k2y()
gy = self.thr.copy_array(self.y)
return gy
def selfadjoint(self, gx):
"""
selfadjoint NUFFT (Teplitz) on the heterogeneous device
:param gx: The input gpu array, with size=Nd
:type: reikna gpu array with dtype =numpy.complex64
:return: gx: The output gpu array, with size=Nd
:rtype: reikna gpu array with dtype =numpy.complex64
"""
self.x_Nd = self.thr.copy_array(gx)
self._x2xx()
self._xx2k()
self._k2y2k()
self._k2xx()
self._xx2x()
gx2 = self.thr.copy_array(self.x_Nd)
return gx2
def guarantee_array(variable):
''' Guarantees that a varaible is a numpy ndarray and supports -, *, +, and other operators
Args:
variable (`number` or `numpy.ndarray`): variable to coalesce
Returns:
(type). Which supports * / and other operations with arrays
'''
if type(variable) in [float, np.ndarray, np.int32, np.int64, np.float32, np.float64, np.complex64, np.complex128]:
return variable
elif type(variable) is int:
return float(variable)
elif type(variable) is list:
return np.asarray(variable)
else:
raise ValueError(f'variable is of invalid type {type(variable)}')
def test_branch_cuts_complex64(self):
# check branch cuts and continuity on them
yield _check_branch_cut, np.log, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.log2, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1, True, np.complex64
yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64
yield _check_branch_cut, np.arccos, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64
yield _check_branch_cut, np.arctan, [0-2j, 2j], [1, 1], -1, 1, True, np.complex64
yield _check_branch_cut, np.arcsinh, [0-2j, 2j], [1, 1], -1, 1, True, np.complex64
yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1, True, np.complex64
yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, 1j], 1, -1, True, np.complex64
# check against bogus branch cuts: assert continuity between quadrants
yield _check_branch_cut, np.arcsin, [0-2j, 2j], [ 1, 1], 1, 1, False, np.complex64
yield _check_branch_cut, np.arccos, [0-2j, 2j], [ 1, 1], 1, 1, False, np.complex64
yield _check_branch_cut, np.arctan, [ -2, 2], [1j, 1j], 1, 1, False, np.complex64
yield _check_branch_cut, np.arcsinh, [ -2, 2, 0], [1j, 1j, 1], 1, 1, False, np.complex64
yield _check_branch_cut, np.arccosh, [0-2j, 2j, 2], [1, 1, 1j], 1, 1, False, np.complex64
yield _check_branch_cut, np.arctanh, [0-2j, 2j, 0], [1, 1, 1j], 1, 1, False, np.complex64
def test_prod(self):
ba = [1, 2, 10, 11, 6, 5, 4]
ba2 = [[1, 2, 3, 4], [5, 6, 7, 9], [10, 3, 4, 5]]
for ctype in [np.int16, np.uint16, np.int32, np.uint32,
np.float32, np.float64, np.complex64, np.complex128]:
a = np.array(ba, ctype)
a2 = np.array(ba2, ctype)
if ctype in ['1', 'b']:
self.assertRaises(ArithmeticError, a.prod)
self.assertRaises(ArithmeticError, a2.prod, axis=1)
else:
assert_equal(a.prod(axis=0), 26400)
assert_array_equal(a2.prod(axis=0),
np.array([50, 36, 84, 180], ctype))
assert_array_equal(a2.prod(axis=-1),
np.array([24, 1890, 600], ctype))
def test_sum_complex(self):
for dt in (np.complex64, np.complex128, np.clongdouble):
for v in (0, 1, 2, 7, 8, 9, 15, 16, 19, 127,
128, 1024, 1235):
tgt = dt(v * (v + 1) / 2) - dt((v * (v + 1) / 2) * 1j)
d = np.empty(v, dtype=dt)
d.real = np.arange(1, v + 1)
d.imag = -np.arange(1, v + 1)
assert_almost_equal(np.sum(d), tgt)
assert_almost_equal(np.sum(d[::-1]), tgt)
d = np.ones(500, dtype=dt) + 1j
assert_almost_equal(np.sum(d[::2]), 250. + 250j)
assert_almost_equal(np.sum(d[1::2]), 250. + 250j)
assert_almost_equal(np.sum(d[::3]), 167. + 167j)
assert_almost_equal(np.sum(d[1::3]), 167. + 167j)
assert_almost_equal(np.sum(d[::-2]), 250. + 250j)
assert_almost_equal(np.sum(d[-1::-2]), 250. + 250j)
assert_almost_equal(np.sum(d[::-3]), 167. + 167j)
assert_almost_equal(np.sum(d[-1::-3]), 167. + 167j)
# sum with first reduction entry != 0
d = np.ones((1,), dtype=dt) + 1j
d += d
assert_almost_equal(d, 2. + 2j)
def test_zero_division(self):
with np.errstate(all="ignore"):
for t in [np.complex64, np.complex128]:
a = t(0.0)
b = t(1.0)
assert_(np.isinf(b/a))
b = t(complex(np.inf, np.inf))
assert_(np.isinf(b/a))
b = t(complex(np.inf, np.nan))
assert_(np.isinf(b/a))
b = t(complex(np.nan, np.inf))
assert_(np.isinf(b/a))
b = t(complex(np.nan, np.nan))
assert_(np.isnan(b/a))
b = t(0.)
assert_(np.isnan(b/a))
def test_basic(self):
ba = [1, 2, 10, 11, 6, 5, 4]
ba2 = [[1, 2, 3, 4], [5, 6, 7, 9], [10, 3, 4, 5]]
for ctype in [np.int8, np.uint8, np.int16, np.uint16, np.int32,
np.uint32, np.float32, np.float64, np.complex64, np.complex128]:
a = np.array(ba, ctype)
a2 = np.array(ba2, ctype)
tgt = np.array([1, 3, 13, 24, 30, 35, 39], ctype)
assert_array_equal(np.cumsum(a, axis=0), tgt)
tgt = np.array(
[[1, 2, 3, 4], [6, 8, 10, 13], [16, 11, 14, 18]], ctype)
assert_array_equal(np.cumsum(a2, axis=0), tgt)
tgt = np.array(
[[1, 3, 6, 10], [5, 11, 18, 27], [10, 13, 17, 22]], ctype)
assert_array_equal(np.cumsum(a2, axis=1), tgt)
def test_basic(self):
ba = [1, 2, 10, 11, 6, 5, 4]
ba2 = [[1, 2, 3, 4], [5, 6, 7, 9], [10, 3, 4, 5]]
for ctype in [np.int16, np.uint16, np.int32, np.uint32,
np.float32, np.float64, np.complex64, np.complex128]:
a = np.array(ba, ctype)
a2 = np.array(ba2, ctype)
if ctype in ['1', 'b']:
self.assertRaises(ArithmeticError, np.prod, a)
self.assertRaises(ArithmeticError, np.prod, a2, 1)
else:
assert_equal(a.prod(axis=0), 26400)
assert_array_equal(a2.prod(axis=0),
np.array([50, 36, 84, 180], ctype))
assert_array_equal(a2.prod(axis=-1),
np.array([24, 1890, 600], ctype))
def test_basic(self):
ba = [1, 2, 10, 11, 6, 5, 4]
ba2 = [[1, 2, 3, 4], [5, 6, 7, 9], [10, 3, 4, 5]]
for ctype in [np.int16, np.uint16, np.int32, np.uint32,
np.float32, np.float64, np.complex64, np.complex128]:
a = np.array(ba, ctype)
a2 = np.array(ba2, ctype)
if ctype in ['1', 'b']:
self.assertRaises(ArithmeticError, np.cumprod, a)
self.assertRaises(ArithmeticError, np.cumprod, a2, 1)
self.assertRaises(ArithmeticError, np.cumprod, a)
else:
assert_array_equal(np.cumprod(a, axis=-1),
np.array([1, 2, 20, 220,
1320, 6600, 26400], ctype))
assert_array_equal(np.cumprod(a2, axis=0),
np.array([[1, 2, 3, 4],
[5, 12, 21, 36],
[50, 36, 84, 180]], ctype))
assert_array_equal(np.cumprod(a2, axis=-1),
np.array([[1, 2, 6, 24],
[5, 30, 210, 1890],
[10, 30, 120, 600]], ctype))
def test_shuffle(self):
# Test lists, arrays (of various dtypes), and multidimensional versions
# of both, c-contiguous or not:
for conv in [lambda x: np.array([]),
lambda x: x,
lambda x: np.asarray(x).astype(np.int8),
lambda x: np.asarray(x).astype(np.float32),
lambda x: np.asarray(x).astype(np.complex64),
lambda x: np.asarray(x).astype(object),
lambda x: [(i, i) for i in x],
lambda x: np.asarray([[i, i] for i in x]),
lambda x: np.vstack([x, x]).T,
# gh-4270
lambda x: np.asarray([(i, i) for i in x],
[("a", object, 1),
("b", np.int32, 1)])]:
np.random.seed(self.seed)
alist = conv([1, 2, 3, 4, 5, 6, 7, 8, 9, 0])
np.random.shuffle(alist)
actual = alist
desired = conv([0, 1, 9, 6, 2, 4, 5, 8, 7, 3])
np.testing.assert_array_equal(actual, desired)
def waveform_to_file(
station=0, clen=10000, oversample=10, filter_output=False,
):
a = rep_seq(
create_pseudo_random_code(clen=clen, seed=station),
rep=oversample,
)
if filter_output:
w = numpy.zeros([oversample * clen], dtype=numpy.complex64)
fl = (int(oversample + (0.1 * oversample)))
w[0:fl] = scipy.signal.blackmanharris(fl)
aa = numpy.fft.ifft(numpy.fft.fft(w) * numpy.fft.fft(a))
a = aa / numpy.max(numpy.abs(aa))
a = numpy.array(a, dtype=numpy.complex64)
a.tofile('code-l%d-b%d-%06df.bin' % (clen, oversample, station))
else:
a.tofile('code-l%d-b%d-%06d.bin' % (clen, oversample, station))
def for_all_dtypes_combination(names=('dtyes',),
no_float16=False, no_bool=False, full=None,
no_complex=False):
"""Decorator that checks the fixture with a product set of all dtypes.
Args:
names(list of str): Argument names to which dtypes are passed.
no_float16(bool): If ``True``, ``numpy.float16`` is
omitted from candidate dtypes.
no_bool(bool): If ``True``, ``numpy.bool_`` is
omitted from candidate dtypes.
full(bool): If ``True``, then all combinations of dtypes
will be tested.
Otherwise, the subset of combinations will be tested
(see description in :func:`cupy.testing.for_dtypes_combination`).
no_complex(bool): If, True, ``numpy.complex64`` and
``numpy.complex128`` are omitted from candidate dtypes.
.. seealso:: :func:`cupy.testing.for_dtypes_combination`
"""
types = _make_all_dtypes(no_float16, no_bool, no_complex)
return for_dtypes_combination(types, names, full)
def addFilter(self,f,recenter=False):
'''
f can be a function f(x,y) is defined over x = (-w/2, w/2] and
y = (-h/2,h/2] and should be centered on the coord 0,0.
TODO: At some point this function should be expanded to take filters
represented by arrays.
'''
if recenter == True:
raise NotImplementedError
if isinstance(f,GaborWavelet):
filt = np.fft.fft2(f.mask(self.tile_size))
self.filters.append(filt)
else:
w,h = self.tile_size
m = np.zeros((w,h),np.complex64)
for x in range(-w/2,w/2):
for y in range(-h/2,h/2):
m[x,y] = f(x,y)
filt = np.fft.fft2(m)
self.filters.append(filt.conj())
def initSize(self,size):
w,h = size
w = int(round(w))
h = int(round(h))
if self.x.has_key((w,h)):
return
# Initializing translations
self.x[(w,h)] = []
for x in range(w):
mat = np.zeros((w,h),dtype=np.complex64)
mat[x,0] = 1.0
filter = np.fft.fft2(mat)
self.x[(w,h)].append(filter)
self.y[(w,h)] = []
for y in range(h):
mat = np.zeros((w,h),dtype=np.complex64)
mat[0,y] = 1.0
filter = np.fft.fft2(mat)
self.y[(w,h)].append(filter)
def log_power_spectrum_extractor(x, win_len, shift_len, win_type, is_log=False):
samples = x.shape[0]
frames = (samples - win_len) // shift_len
stft = np.zeros((win_len, frames), dtype=np.complex64)
spect = np.zeros((win_len // 2 + 1, frames), dtype=np.float64)
if win_type == 'hanning':
window = np.hanning(win_len)
elif win_type == 'hamming':
window = np.hamming(win_len)
elif win_type == 'rectangle':
window = np.ones(win_len)
for i in range(frames):
one_frame = x[i*shift_len: i*shift_len+win_len]
windowed_frame = np.multiply(one_frame, window)
stft[:, i] = np.fft.fft(windowed_frame, win_len)
if is_log:
spect[:, i] = np.log(np.power(np.abs(stft[0: win_len//2+1, i]), 2.))
else:
spect[:, i] = np.power(np.abs(stft[0: win_len//2+1, i]), 2.)
return spect
def stft_extractor(x, win_len, shift_len, win_type):
samples = x.shape[0]
frames = (samples - win_len) // shift_len
stft = np.zeros((win_len, frames), dtype=np.complex64)
spect = np.zeros((win_len // 2 + 1, frames), dtype=np.complex64)
if win_type == 'hanning':
window = np.hanning(win_len)
elif win_type == 'hamming':
window = np.hamming(win_len)
elif win_type == 'rectangle':
window = np.ones(win_len)
for i in range(frames):
one_frame = x[i*shift_len: i*shift_len+win_len]
windowed_frame = np.multiply(one_frame, window)
stft[:, i] = np.fft.fft(windowed_frame, win_len)
spect[:, i] = stft[: win_len//2+1, i]
return spect
def comp_ola_deconv(fs_gpu, ys_gpu, L_gpu, alpha, beta):
"""
Computes the division in Fourier space needed for direct deconvolution
"""
sfft = fs_gpu.shape
block_size = (16,16,1)
grid_size = (int(np.ceil(np.float32(sfft[0]*sfft[1])/block_size[0])),
int(np.ceil(np.float32(sfft[2])/block_size[1])))
mod = cu.module_from_buffer(cubin)
comp_ola_deconv_Kernel = mod.get_function("comp_ola_deconv_Kernel")
z_gpu = cua.zeros(sfft, np.complex64)
comp_ola_deconv_Kernel(z_gpu.gpudata,
np.int32(sfft[0]), np.int32(sfft[1]), np.int32(sfft[2]),
fs_gpu.gpudata, ys_gpu.gpudata, L_gpu.gpudata,
np.float32(alpha), np.float32(beta),
block=block_size, grid=grid_size)
return z_gpu
def comp_ola_gdeconv(xx_gpu, xy_gpu, yx_gpu, yy_gpu, L_gpu, alpha, beta):
"""
Computes the division in Fourier space needed for gdirect deconvolution
"""
sfft = xx_gpu.shape
block_size = (16,16,1)
grid_size = (int(np.ceil(np.float32(sfft[0]*sfft[1])/block_size[0])),
int(np.ceil(np.float32(sfft[2])/block_size[1])))
mod = cu.module_from_buffer(cubin)
comp_ola_gdeconv_Kernel = mod.get_function("comp_ola_gdeconv_Kernel")
z_gpu = cua.zeros(sfft, np.complex64)
comp_ola_gdeconv_Kernel(z_gpu.gpudata,
np.int32(sfft[0]), np.int32(sfft[1]), np.int32(sfft[2]),
xx_gpu, xy_gpu, yx_gpu, yy_gpu, L_gpu.gpudata,
np.float32(alpha), np.float32(beta),
block=block_size, grid=grid_size)
return z_gpu
def comp_ola_sdeconv(gx_gpu, gy_gpu, xx_gpu, xy_gpu, Ftpy_gpu, f_gpu, L_gpu, alpha, beta, gamma=0):
"""
Computes the division in Fourier space needed for sparse deconvolution
"""
sfft = xx_gpu.shape
block_size = (16,16,1)
grid_size = (int(np.ceil(np.float32(sfft[0]*sfft[1])/block_size[0])),
int(np.ceil(np.float32(sfft[2])/block_size[1])))
mod = cu.module_from_buffer(cubin)
comp_ola_sdeconv_Kernel = mod.get_function("comp_ola_sdeconv_Kernel")
z_gpu = cua.zeros(sfft, np.complex64)
comp_ola_sdeconv_Kernel(z_gpu.gpudata,
np.int32(sfft[0]), np.int32(sfft[1]), np.int32(sfft[2]),
gx_gpu.gpudata, gy_gpu.gpudata,
xx_gpu.gpudata, xy_gpu.gpudata,
Ftpy_gpu.gpudata, f_gpu.gpudata, L_gpu.gpudata,
np.float32(alpha), np.float32(beta),
np.float32(gamma),
block=block_size, grid=grid_size)
return z_gpu
def fst_delay_snd(fst, snd, samp_rate):
# Verify argument shape.
s1, s2 = fst.shape, snd.shape
if len(s1) != 1 or len(s2) != 1 or s1[0] != s2[0]:
raise Exception("Argument shape invalid, in 'fst_delay_snd' function")
length = s1[0]
half_len = int(length / 2)
Xfst = numpy.fft.fft(fst)
Xsnd_star = numpy.conj(numpy.fft.fft(snd))
Xall = numpy.zeros(length, dtype=numpy.complex64)
for i in range(0, length):
if Xsnd_star[i] == 0 or Xfst[i] == 0:
Xall[i] = 0
else:
Xall[i] = (Xsnd_star[i] * Xfst[i]) / abs(Xsnd_star[i]) / abs(Xfst[i])
R = numpy.fft.ifft(Xall)
max_pos = numpy.argmax(R)
if max_pos > half_len:
return -(length - 1 - max_pos) / samp_rate
else:
return max_pos / samp_rate
def fst_delay_snd(fst, snd, samp_rate):
# Verify argument shape.
s1, s2 = fst.shape, snd.shape
if len(s1) != 1 or len(s2) != 1 or s1[0] != s2[0]:
raise Exception("Argument shape invalid, in 'fst_delay_snd' function")
length = s1[0]
half_len = int(length / 2)
Xfst = numpy.fft.fft(fst)
Xsnd = numpy.fft.fft(snd)
Xsnd_star = numpy.conj(Xsnd)
Xall = numpy.zeros(length, dtype=numpy.complex64)
for i in range(0, length):
Xall[i] = (Xsnd_star[i] * Xfst[i]) / abs(Xsnd_star[i]) / abs(Xfst[i])
R = numpy.fft.ifft(Xall)
max_pos = numpy.argmax(R)
if max_pos > half_len:
delta = -(length - 1 - max_pos) / samp_rate
else:
delta = max_pos / samp_rate
return delta, R, Xfst, Xsnd
