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 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) Example 2
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
except KeyError:
return None Example 3
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 Example 4
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 Example 5
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 Example 6
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 Example 7
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 Example 8
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 Example 9
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 Example 10
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 Example 11
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)}') Example 12
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 Example 13
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)) Example 14
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) Example 15
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)) Example 16
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) Example 17
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)) Example 18
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)) Example 19
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) Example 20
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)) Example 21
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) Example 22
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()) Example 23
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) Example 24
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 Example 25
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 Example 26
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 Example 27
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 Example 28
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 Example 29
def pad_cpu2gpu(x, sz, offset=(0,0), dtype='real'):
block_size = (16, 16 ,1)
grid_size = (int(np.ceil(np.float32(sz[1])/block_size[1])),
int(np.ceil(np.float32(sz[0])/block_size[0])))
sx = x.shape
if x.__class__ == np.ndarray:
x = np.array(x).astype(np.float32)
x_gpu = cua.to_gpu(x)
elif x.__class__ == cua.GPUArray:
x_gpu = x
if dtype == 'real':
mod = cu.module_from_buffer(cubin)
zeroPadKernel = mod.get_function("zeroPadKernel")
x_padded_gpu = cua.zeros(tuple((int(sz[0]),int(sz[1]))), np.float32)
zeroPadKernel(x_padded_gpu.gpudata, np.int32(sz[0]), np.int32(sz[1]),
x_gpu.gpudata, np.int32(sx[0]), np.int32(sx[1]),
np.int32(offset[0]), np.int32(offset[1]),
block=block_size, grid=grid_size)
elif dtype == 'complex':
mod = cu.module_from_buffer(cubin)
#mod = SourceModule(open('gputools.cu').read(), keep=True)
zeroPadComplexKernel = mod.get_function("zeroPadComplexKernel")
x_padded_gpu = cua.zeros(tuple((int(sz[0]),int(sz[1]))), np.complex64)
zeroPadComplexKernel(x_padded_gpu.gpudata, np.int32(sz[0]), np.int32(sz[1]),
x_gpu.gpudata, np.int32(sx[0]), np.int32(sx[1]),
np.int32(offset[0]), np.int32(offset[1]),
block=block_size, grid=grid_size)
return x_padded_gpu Example 30
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 Example 31
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 Example 32
def __init__(self, test_array, complex_numbers=False):
"""Figure out data settings from the test array.
@param[in] test_array A list or numpy array containing test data"""
super(TestingBlock, self).__init__()
if isinstance(test_array, np.ndarray):
if test_array.dtype == np.complex64:
complex_numbers = True
if complex_numbers:
self.test_array = np.array(test_array).astype(np.complex64)
header = {
'nbit': 64,
'dtype': 'complex64',
'shape': self.test_array.shape}
self.dtype = np.complex64
else:
self.test_array = np.array(test_array).astype(np.float32)
header = {
'nbit': 32,
'dtype': 'float32',
'shape': self.test_array.shape}
self.dtype = np.float32
self.output_header = json.dumps(header) Example 33
def main(self, input_rings, output_rings):
"""
@param[in] input_rings First ring in this list will be used for
data
@param[out] output_rings First ring in this list will be used for
data output."""
data_accumulate = None
for ispan in self.iterate_ring_read(input_rings[0]):
if self.nbit < 8:
unpacked_data = unpack(ispan.data_view(self.dtype), self.nbit)
else:
unpacked_data = ispan.data_view(self.dtype)
if data_accumulate is not None:
data_accumulate = np.concatenate((data_accumulate, unpacked_data[0]))
else:
data_accumulate = unpacked_data[0]
if self.shape != [1, 1]:
data_accumulate = np.reshape(data_accumulate, (self.shape[0], -1))
data_accumulate = data_accumulate.astype(np.complex64)
self.out_gulp_size = data_accumulate.nbytes
outspan_generator = self.iterate_ring_write(output_rings[0])
ospan = outspan_generator.next()
result = np.fft.fft(data_accumulate).astype(np.complex64)
ospan.data_view(np.complex64)[0] = result.ravel() Example 34
def main(self, input_rings, output_rings):
"""
@param[in] input_rings First ring in this list will be used for
data input.
@param[out] output_rings First ring in this list will be used for
data output."""
data_accumulate = None
for ispan in self.iterate_ring_read(input_rings[0]):
if self.nbit < 8:
unpacked_data = unpack(ispan.data_view(self.dtype), self.nbit)
else:
unpacked_data = ispan.data_view(self.dtype)
if data_accumulate is not None:
data_accumulate = np.concatenate((data_accumulate, unpacked_data[0]))
else:
data_accumulate = unpacked_data[0]
data_accumulate = data_accumulate.astype(np.complex64)
self.out_gulp_size = data_accumulate.nbytes
outspan_generator = self.iterate_ring_write(output_rings[0])
ospan = outspan_generator.next()
result = np.fft.ifft(data_accumulate)
ospan.data_view(np.complex64)[0][:] = result[:] Example 35
def numpy2bifrost(dtype):
if dtype == np.int8: return _bf.BF_DTYPE_I8
elif dtype == np.int16: return _bf.BF_DTYPE_I16
elif dtype == np.int32: return _bf.BF_DTYPE_I32
elif dtype == np.uint8: return _bf.BF_DTYPE_U8
elif dtype == np.uint16: return _bf.BF_DTYPE_U16
elif dtype == np.uint32: return _bf.BF_DTYPE_U32
elif dtype == np.float16: return _bf.BF_DTYPE_F16
elif dtype == np.float32: return _bf.BF_DTYPE_F32
elif dtype == np.float64: return _bf.BF_DTYPE_F64
elif dtype == np.float128: return _bf.BF_DTYPE_F128
elif dtype == ci8: return _bf.BF_DTYPE_CI8
elif dtype == ci16: return _bf.BF_DTYPE_CI16
elif dtype == ci32: return _bf.BF_DTYPE_CI32
elif dtype == cf16: return _bf.BF_DTYPE_CF16
elif dtype == np.complex64: return _bf.BF_DTYPE_CF32
elif dtype == np.complex128: return _bf.BF_DTYPE_CF64
elif dtype == np.complex256: return _bf.BF_DTYPE_CF128
else: raise ValueError("Unsupported dtype: " + str(dtype)) Example 36
def numpy2string(dtype):
if dtype == np.int8: return 'i8'
elif dtype == np.int16: return 'i16'
elif dtype == np.int32: return 'i32'
elif dtype == np.int64: return 'i64'
elif dtype == np.uint8: return 'u8'
elif dtype == np.uint16: return 'u16'
elif dtype == np.uint32: return 'u32'
elif dtype == np.uint64: return 'u64'
elif dtype == np.float16: return 'f16'
elif dtype == np.float32: return 'f32'
elif dtype == np.float64: return 'f64'
elif dtype == np.float128: return 'f128'
elif dtype == np.complex64: return 'cf32'
elif dtype == np.complex128: return 'cf64'
elif dtype == np.complex256: return 'cf128'
else: raise TypeError("Unsupported dtype: " + str(dtype)) Example 37
def run_test_matmul_aa_ci8_shape(self, shape, transpose=False):
# **TODO: This currently never triggers the transpose path in the backend
shape_complex = shape[:-1] + (shape[-1] * 2,)
# Note: The xGPU-like correlation kernel does not support input values of -128 (only [-127:127])
a8 = ((np.random.random(size=shape_complex) * 2 - 1) * 127).astype(np.int8)
a_gold = a8.astype(np.float32).view(np.complex64)
if transpose:
a_gold = H(a_gold)
# Note: np.matmul seems to be slow and inaccurate when there are batch dims
c_gold = np.matmul(a_gold, H(a_gold))
triu = np.triu_indices(shape[-2] if not transpose else shape[-1], 1)
c_gold[..., triu[0], triu[1]] = 0
a = a8.view(bf.DataType.ci8)
a = bf.asarray(a, space='cuda')
if transpose:
a = H(a)
c = bf.zeros_like(c_gold, space='cuda')
self.linalg.matmul(1, a, None, 0, c)
c = c.copy('system')
np.testing.assert_allclose(c, c_gold, RTOL, ATOL) Example 38
def run_test_matmul_ab_ci8_shape(self, shape, k, transpose=False):
ashape_complex = shape[:-2] + (shape[-2], k * 2)
bshape_complex = shape[:-2] + (k, shape[-1] * 2)
a8 = (np.random.random(size=ashape_complex) * 255).astype(np.int8)
b8 = (np.random.random(size=bshape_complex) * 255).astype(np.int8)
a_gold = a8.astype(np.float32).view(np.complex64)
b_gold = b8.astype(np.float32).view(np.complex64)
if transpose:
a_gold, b_gold = H(b_gold), H(a_gold)
c_gold = np.matmul(a_gold, b_gold)
a = a8.view(bf.DataType.ci8)
b = b8.view(bf.DataType.ci8)
a = bf.asarray(a, space='cuda')
b = bf.asarray(b, space='cuda')
if transpose:
a, b = H(b), H(a)
c = bf.zeros_like(c_gold, space='cuda')
self.linalg.matmul(1, a, b, 0, c)
c = c.copy('system')
np.testing.assert_allclose(c, c_gold, RTOL, ATOL) Example 39
def run_benchmark_matmul_aa_correlator_kernel(self, ntime, nstand, nchan):
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)
b_gold = np.matmul(H(x[:,[0],:]), x[:,[0],:])
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)
b = bf.zeros_like(b_gold, space='cuda')
bf.device.stream_synchronize();
t0 = time.time()
nrep = 200
for _ in xrange(nrep):
self.linalg.matmul(1, None, x, 0, b)
bf.device.stream_synchronize();
dt = time.time() - t0
nflop = nrep * nchan * ntime * nstand*(nstand+1)/2 * 2*2 * 8
print nstand, '\t', nflop / dt / 1e9, 'GFLOP/s'
print '\t\t', nrep*ntime*nchan / dt / 1e6, 'MHz' Example 40
def run_test_c2c_impl(self, shape, axes, inverse=False, fftshift=False):
shape = list(shape)
shape[-1] *= 2 # For complex
known_data = np.random.normal(size=shape).astype(np.float32).view(np.complex64)
idata = bf.ndarray(known_data, space='cuda')
odata = bf.empty_like(idata)
fft = Fft()
fft.init(idata, odata, axes=axes, apply_fftshift=fftshift)
fft.execute(idata, odata, inverse)
if inverse:
if fftshift:
known_data = np.fft.ifftshift(known_data, axes=axes)
# Note: Numpy applies normalization while CUFFT does not
norm = reduce(lambda a, b: a * b, [known_data.shape[d]
for d in axes])
known_result = gold_ifftn(known_data, axes=axes) * norm
else:
known_result = gold_fftn(known_data, axes=axes)
if fftshift:
known_result = np.fft.fftshift(known_result, axes=axes)
x = (np.abs(odata.copy('system') - known_result) / known_result > RTOL).astype(np.int32)
a = odata.copy('system')
b = known_result
compare(odata.copy('system'), known_result) Example 41
def run_test_c2r_impl(self, shape, axes, fftshift=False):
ishape = list(shape)
oshape = list(shape)
ishape[axes[-1]] = shape[axes[-1]] // 2 + 1
oshape[axes[-1]] = (ishape[axes[-1]] - 1) * 2
ishape[-1] *= 2 # For complex
known_data = np.random.normal(size=ishape).astype(np.float32).view(np.complex64)
idata = bf.ndarray(known_data, space='cuda')
odata = bf.ndarray(shape=oshape, dtype='f32', space='cuda')
fft = Fft()
fft.init(idata, odata, axes=axes, apply_fftshift=fftshift)
fft.execute(idata, odata)
# Note: Numpy applies normalization while CUFFT does not
norm = reduce(lambda a, b: a * b, [shape[d] for d in axes])
if fftshift:
known_data = np.fft.ifftshift(known_data, axes=axes)
known_result = gold_irfftn(known_data, axes=axes) * norm
compare(odata.copy('system'), known_result) Example 42
def test_data_sizes(self):
"""Test that different number of bits give correct throughput size"""
for iterate in range(5):
nbit = 2**iterate
if nbit == 8:
continue
self.blocks[0] = (
SigprocReadBlock(
'./data/2chan' + str(nbit) + 'bitNoDM.fil'),
[], [0])
open(self.logfile, 'w').close()
Pipeline(self.blocks).main()
number_fftd = np.loadtxt(self.logfile).astype(np.float32).view(np.complex64).size
# Compare with simple copy
self.blocks[1] = (CopyBlock(), [0], [1])
open(self.logfile, 'w').close()
Pipeline(self.blocks).main()
number_copied = np.loadtxt(self.logfile).size
self.assertEqual(number_fftd, number_copied)
# Go back to FFT
self.blocks[1] = (FFTBlock(gulp_size=4096 * 8 * 8 * 8), [0], [1]) Example 43
def test_equivalent_data_to_copy(self):
"""Test that the data coming out of this pipeline is equivalent
the initial read data"""
self.logfile = '.log.txt'
self.blocks = []
self.blocks.append((
SigprocReadBlock(
'./data/1chan8bitNoDM.fil'),
[], [0]))
self.blocks.append((FFTBlock(gulp_size=4096 * 8 * 8 * 8 * 8), [0], [1]))
self.blocks.append((IFFTBlock(gulp_size=4096 * 8 * 8 * 8 * 8), [1], [2]))
self.blocks.append((WriteAsciiBlock(self.logfile), [2], []))
open(self.logfile, 'w').close()
Pipeline(self.blocks).main()
unfft_result = np.loadtxt(self.logfile).astype(np.float32).view(np.complex64)
self.blocks[1] = (CopyBlock(), [0], [1])
self.blocks[2] = (WriteAsciiBlock(self.logfile), [1], [])
del self.blocks[3]
open(self.logfile, 'w').close()
Pipeline(self.blocks).main()
untouched_result = np.loadtxt(self.logfile).astype(np.float32)
np.testing.assert_almost_equal(unfft_result, untouched_result, 2) Example 44
def default(self, obj):
# convert dates and numpy objects in a json serializable format
if isinstance(obj, datetime):
return obj.strftime('%Y-%m-%dT%H:%M:%SZ')
elif isinstance(obj, date):
return obj.strftime('%Y-%m-%d')
elif type(obj) in (np.int_, np.intc, np.intp, np.int8, np.int16,
np.int32, np.int64, np.uint8, np.uint16,
np.uint32, np.uint64):
return int(obj)
elif type(obj) in (np.bool_,):
return bool(obj)
elif type(obj) in (np.float_, np.float16, np.float32, np.float64,
np.complex_, np.complex64, np.complex128):
return float(obj)
# Let the base class default method raise the TypeError
return json.JSONEncoder.default(self, obj) Example 45
def testTwoSessions(self):
optimizer = ctf.train.CplxAdamOptimizer()
g = tf.Graph()
with g.as_default():
with tf.Session():
var0 = tf.Variable(np.array([1.0+1.0j, 2.0+2.0j], dtype=np.complex64),
name="v0")
grads0 = tf.constant(np.array([0.1+0.1j, 0.1+0.1j], dtype=np.complex64))
optimizer.apply_gradients([(grads0, var0)])
gg = tf.Graph()
with gg.as_default():
with tf.Session():
var0 = tf.Variable(np.array([1.0+1.0j, 2.0+2.0j], dtype=np.complex64),
name="v0")
grads0 = tf.constant(np.array([0.1+0.1j, 0.1+0.1j], dtype=np.complex64))
# If the optimizer saves any state not keyed by graph the following line
# fails.
optimizer.apply_gradients([(grads0, var0)]) Example 46
def _testTypes(self, vals):
for dtype in [np.complex64]:
x = np.zeros(vals.shape).astype(dtype)
y = vals.astype(dtype)
var_value, op_value = self._initAssignFetch(x, y, use_gpu=False)
self.assertAllEqual(y, var_value)
self.assertAllEqual(y, op_value)
var_value, op_value = self._initAssignAddFetch(x, y, use_gpu=False)
self.assertAllEqual(x + y, var_value)
self.assertAllEqual(x + y, op_value)
var_value, op_value = self._initAssignSubFetch(x, y, use_gpu=False)
self.assertAllEqual(x - y, var_value)
self.assertAllEqual(x - y, op_value)
if tf.test.is_built_with_cuda() and dtype in [np.float32, np.float64]:
var_value, op_value = self._initAssignFetch(x, y, use_gpu=True)
self.assertAllEqual(y, var_value)
self.assertAllEqual(y, op_value)
var_value, op_value = self._initAssignAddFetch(x, y, use_gpu=True)
self.assertAllEqual(x + y, var_value)
self.assertAllEqual(x + y, op_value)
var_value, op_value = self._initAssignSubFetch(x, y, use_gpu=False)
self.assertAllEqual(x - y, var_value)
self.assertAllEqual(x - y, op_value) Example 47
def testHStack(self):
with self.test_session(force_gpu=True):
p1 = array_ops.placeholder(dtypes.complex64, shape=[4, 4])
p2 = array_ops.placeholder(dtypes.complex64, shape=[4, 4])
c = array_ops.concat([p1, p2], 0)
params = {
p1: (np.random.rand(4, 4) +
1j*np.random.rand(4, 4)).astype(np.complex64),
p2: (np.random.rand(4, 4) +
1j*np.random.rand(4, 4)).astype(np.complex64),
}
result = c.eval(feed_dict=params)
self.assertEqual(result.shape, c.get_shape())
self.assertAllEqual(result[:4, :], params[p1])
self.assertAllEqual(result[4:, :], params[p2]) Example 48
def testVStack(self):
with self.test_session(force_gpu=True):
p1 = array_ops.placeholder(dtypes.complex64, shape=[4, 4])
p2 = array_ops.placeholder(dtypes.complex64, shape=[4, 4])
c = array_ops.concat([p1, p2], 1)
params = {
p1: (np.random.rand(4, 4) +
1j*np.random.rand(4, 4)).astype(np.complex64),
p2: (np.random.rand(4, 4) +
1j*np.random.rand(4, 4)).astype(np.complex64),
}
result = c.eval(feed_dict=params)
self.assertEqual(result.shape, c.get_shape())
self.assertAllEqual(result[:, :4], params[p1])
self.assertAllEqual(result[:, 4:], params[p2]) Example 49
def testGradientWithUnknownInputDim(self):
with self.test_session(use_gpu=True):
x = array_ops.placeholder(dtypes.complex64)
y = array_ops.placeholder(dtypes.complex64)
c = array_ops.concat([x, y], 2)
output_shape = [10, 2, 9]
grad_inp = (np.random.rand(*output_shape) +
1j*np.random.rand(*output_shape)).astype(np.complex64)
grad_tensor = constant_op.constant(
[inp for inp in grad_inp.flatten()], shape=output_shape)
grad = gradients_impl.gradients([c], [x, y], [grad_tensor])
concated_grad = array_ops.concat(grad, 2)
params = {
x: (np.random.rand(10, 2, 3) +
1j*np.random.rand(10, 2, 3)).astype(np.complex64),
y: (np.random.rand(10, 2, 6) +
1j*np.random.rand(10, 2, 6)).astype(np.complex64),
}
result = concated_grad.eval(feed_dict=params)
self.assertAllEqual(result, grad_inp) Example 50
def testShapeWithUnknownConcatDim(self):
p1 = array_ops.placeholder(dtypes.complex64)
c1 = constant_op.constant(np.complex64(10.0+0j), shape=[4, 4, 4, 4])
p2 = array_ops.placeholder(dtypes.complex64)
c2 = constant_op.constant(np.complex64(20.0+0j), shape=[4, 4, 4, 4])
dim = array_ops.placeholder(dtypes.int32)
concat = array_ops.concat([p1, c1, p2, c2], dim)
self.assertEqual(4, concat.get_shape().ndims)
# All dimensions unknown.
concat2 = array_ops.concat([p1, p2], dim)
self.assertEqual(None, concat2.get_shape())
# Rank doesn't match.
c3 = constant_op.constant(np.complex64(30.0+0j), shape=[4, 4, 4])
with self.assertRaises(ValueError):
array_ops.concat([p1, c1, p2, c3], dim) 