|
| 1 | +#!/usr/bin/env python |
| 2 | +# Copyright (c) 2019-2020, Intel Corporation |
| 3 | +# |
| 4 | +# Redistribution and use in source and binary forms, with or without |
| 5 | +# modification, are permitted provided that the following conditions are met: |
| 6 | +# |
| 7 | +# * Redistributions of source code must retain the above copyright notice, |
| 8 | +# this list of conditions and the following disclaimer. |
| 9 | +# * Redistributions in binary form must reproduce the above copyright |
| 10 | +# notice, this list of conditions and the following disclaimer in the |
| 11 | +# documentation and/or other materials provided with the distribution. |
| 12 | +# * Neither the name of Intel Corporation nor the names of its contributors |
| 13 | +# may be used to endorse or promote products derived from this software |
| 14 | +# without specific prior written permission. |
| 15 | +# |
| 16 | +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" |
| 17 | +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
| 18 | +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE |
| 19 | +# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE |
| 20 | +# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
| 21 | +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR |
| 22 | +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER |
| 23 | +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, |
| 24 | +# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE |
| 25 | +# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| 26 | + |
| 27 | +from . import _pydfti |
| 28 | +from . import _float_utils |
| 29 | +import mkl |
| 30 | + |
| 31 | +import scipy.fft as _fft |
| 32 | + |
| 33 | +# Complete the namespace (these are not actually used in this module) |
| 34 | +from scipy.fft import ( |
| 35 | + dct, idct, dst, idst, dctn, idctn, dstn, idstn, |
| 36 | + hfft2, ihfft2, hfftn, ihfftn, |
| 37 | + fftshift, ifftshift, fftfreq, rfftfreq, |
| 38 | + get_workers, set_workers |
| 39 | +) |
| 40 | + |
| 41 | +from numpy.core import (array, asarray, shape, conjugate, take, sqrt, prod) |
| 42 | + |
| 43 | +__all__ = ['fft', 'ifft', 'fft2', 'ifft2', 'fftn', 'ifftn', |
| 44 | + 'rfft', 'irfft', 'rfft2', 'irfft2', 'rfftn', 'irfftn', |
| 45 | + 'hfft', 'ihfft', 'hfft2', 'ihfft2', 'hfftn', 'ihfftn', |
| 46 | + 'dct', 'idct', 'dst', 'idst', 'dctn', 'idctn', 'dstn', 'idstn', |
| 47 | + 'fftshift', 'ifftshift', 'fftfreq', 'rfftfreq', 'get_workers', |
| 48 | + 'set_workers', 'next_fast_len'] |
| 49 | + |
| 50 | +__ua_domain__ = 'numpy.scipy.fft' |
| 51 | +__implemented = dict() |
| 52 | + |
| 53 | +def __ua_function__(method, args, kwargs): |
| 54 | + """Fetch registered UA function.""" |
| 55 | + fn = __implemented.get(method, None) |
| 56 | + if fn is None: |
| 57 | + return NotImplemented |
| 58 | + return fn(*args, **kwargs) |
| 59 | + |
| 60 | + |
| 61 | +def _implements(scipy_func): |
| 62 | + """Decorator adds function to the dictionary of implemented UA functions""" |
| 63 | + def inner(func): |
| 64 | + __implemented[scipy_func] = func |
| 65 | + return func |
| 66 | + |
| 67 | + return inner |
| 68 | + |
| 69 | + |
| 70 | +def _unitary(norm): |
| 71 | + if norm not in (None, "ortho"): |
| 72 | + raise ValueError("Invalid norm value %s, should be None or \"ortho\"." |
| 73 | + % norm) |
| 74 | + return norm is not None |
| 75 | + |
| 76 | + |
| 77 | +def _cook_nd_args(a, s=None, axes=None, invreal=0): |
| 78 | + if s is None: |
| 79 | + shapeless = 1 |
| 80 | + if axes is None: |
| 81 | + s = list(a.shape) |
| 82 | + else: |
| 83 | + s = take(a.shape, axes) |
| 84 | + else: |
| 85 | + shapeless = 0 |
| 86 | + s = list(s) |
| 87 | + if axes is None: |
| 88 | + axes = list(range(-len(s), 0)) |
| 89 | + if len(s) != len(axes): |
| 90 | + raise ValueError("Shape and axes have different lengths.") |
| 91 | + if invreal and shapeless: |
| 92 | + s[-1] = (a.shape[axes[-1]] - 1) * 2 |
| 93 | + return s, axes |
| 94 | + |
| 95 | + |
| 96 | +def _tot_size(x, axes): |
| 97 | + s = x.shape |
| 98 | + if axes is None: |
| 99 | + return x.size |
| 100 | + return prod([s[ai] for ai in axes]) |
| 101 | + |
| 102 | + |
| 103 | +def _workers_to_num_threads(w): |
| 104 | + if w is None: |
| 105 | + return mkl.domain_get_max_threads(domain='fft') |
| 106 | + return int(w) |
| 107 | + |
| 108 | + |
| 109 | +class Workers: |
| 110 | + def __init__(self, workers): |
| 111 | + self.workers = workers |
| 112 | + self.n_threads = _workers_to_num_threads(workers) |
| 113 | + |
| 114 | + def __enter__(self): |
| 115 | + try: |
| 116 | + mkl.domain_set_num_threads(self.n_threads, domain='fft') |
| 117 | + except: |
| 118 | + raise ValueError("Class argument {} result in invalid number of threads {}".format(self.workers, self.n_threads)) |
| 119 | + |
| 120 | + def __exit__(self, *args): |
| 121 | + # restore default |
| 122 | + max_num_threads = mkl.domain_get_max_threads(domain='fft') |
| 123 | + mkl.domain_set_num_threads(max_num_threads, domain='fft') |
| 124 | + |
| 125 | + |
| 126 | +@_implements(_fft.fft) |
| 127 | +def fft(a, n=None, axis=-1, norm=None, overwrite_x=False, workers=None): |
| 128 | + x = _float_utils.__upcast_float16_array(a) |
| 129 | + with Workers(workers): |
| 130 | + output = _pydfti.fft(x, n=n, axis=axis, overwrite_x=overwrite_x) |
| 131 | + if _unitary(norm): |
| 132 | + output *= 1 / sqrt(output.shape[axis]) |
| 133 | + return output |
| 134 | + |
| 135 | + |
| 136 | +@_implements(_fft.ifft) |
| 137 | +def ifft(a, n=None, axis=-1, norm=None, overwrite_x=False, workers=None): |
| 138 | + x = _float_utils.__upcast_float16_array(a) |
| 139 | + with Workers(workers): |
| 140 | + output = _pydfti.ifft(x, n=n, axis=axis, overwrite_x=overwrite_x) |
| 141 | + if _unitary(norm): |
| 142 | + output *= sqrt(output.shape[axis]) |
| 143 | + return output |
| 144 | + |
| 145 | + |
| 146 | +@_implements(_fft.fft2) |
| 147 | +def fft2(a, s=None, axes=(-2,-1), norm=None, overwrite_x=False, workers=None): |
| 148 | + x = _float_utils.__upcast_float16_array(a) |
| 149 | + with Workers(workers): |
| 150 | + output = _pydfti.fftn(x, shape=s, axes=axes, overwrite_x=overwrite_x) |
| 151 | + if _unitary(norm): |
| 152 | + factor = 1 |
| 153 | + for axis in axes: |
| 154 | + factor *= 1 / sqrt(output.shape[axis]) |
| 155 | + output *= factor |
| 156 | + return output |
| 157 | + |
| 158 | + |
| 159 | +@_implements(_fft.ifft2) |
| 160 | +def ifft2(a, s=None, axes=(-2,-1), norm=None, overwrite_x=False, workers=None): |
| 161 | + x = _float_utils.__upcast_float16_array(a) |
| 162 | + with Workers(workers): |
| 163 | + output = _pydfti.ifftn(x, shape=s, axes=axes, overwrite_x=overwrite_x) |
| 164 | + if _unitary(norm): |
| 165 | + factor = 1 |
| 166 | + _axes = range(output.ndim) if axes is None else axes |
| 167 | + for axis in _axes: |
| 168 | + factor *= sqrt(output.shape[axis]) |
| 169 | + output *= factor |
| 170 | + return output |
| 171 | + |
| 172 | + |
| 173 | +@_implements(_fft.fftn) |
| 174 | +def fftn(a, s=None, axes=None, norm=None, overwrite_x=False, workers=None): |
| 175 | + x = _float_utils.__upcast_float16_array(a) |
| 176 | + with Workers(workers): |
| 177 | + output = _pydfti.fftn(x, shape=s, axes=axes, overwrite_x=overwrite_x) |
| 178 | + if _unitary(norm): |
| 179 | + factor = 1 |
| 180 | + _axes = range(output.ndim) if axes is None else axes |
| 181 | + for axis in _axes: |
| 182 | + factor *= 1 / sqrt(output.shape[axis]) |
| 183 | + output *= factor |
| 184 | + return output |
| 185 | + |
| 186 | + |
| 187 | +@_implements(_fft.ifftn) |
| 188 | +def ifftn(a, s=None, axes=None, norm=None, overwrite_x=False, workers=None): |
| 189 | + x = _float_utils.__upcast_float16_array(a) |
| 190 | + with Workers(workers): |
| 191 | + output = _pydfti.ifftn(x, shape=s, axes=axes, overwrite_x=overwrite_x) |
| 192 | + if _unitary(norm): |
| 193 | + factor = 1 |
| 194 | + _axes = range(output.ndim) if axes is None else axes |
| 195 | + for axis in _axes: |
| 196 | + factor *= sqrt(output.shape[axis]) |
| 197 | + output *= factor |
| 198 | + return output |
| 199 | + |
| 200 | + |
| 201 | +@_implements(_fft.rfft) |
| 202 | +def rfft(a, n=None, axis=-1, norm=None, workers=None): |
| 203 | + x = _float_utils.__upcast_float16_array(a) |
| 204 | + unitary = _unitary(norm) |
| 205 | + x = _float_utils.__downcast_float128_array(x) |
| 206 | + if unitary and n is None: |
| 207 | + x = asarray(x) |
| 208 | + n = x.shape[axis] |
| 209 | + with Workers(workers): |
| 210 | + output = _pydfti.rfft_numpy(x, n=n, axis=axis) |
| 211 | + if unitary: |
| 212 | + output *= 1 / sqrt(n) |
| 213 | + return output |
| 214 | + |
| 215 | + |
| 216 | +@_implements(_fft.irfft) |
| 217 | +def irfft(a, n=None, axis=-1, norm=None, workers=None): |
| 218 | + x = _float_utils.__upcast_float16_array(a) |
| 219 | + x = _float_utils.__downcast_float128_array(x) |
| 220 | + with Workers(workers): |
| 221 | + output = _pydfti.irfft_numpy(x, n=n, axis=axis) |
| 222 | + if _unitary(norm): |
| 223 | + output *= sqrt(output.shape[axis]) |
| 224 | + return output |
| 225 | + |
| 226 | + |
| 227 | +@_implements(_fft.rfft2) |
| 228 | +def rfft2(a, s=None, axes=(-2, -1), norm=None, workers=None): |
| 229 | + x = _float_utils.__upcast_float16_array(a) |
| 230 | + x = _float_utils.__downcast_float128_array(a) |
| 231 | + return rfftn(x, s, axes, norm, workers) |
| 232 | + |
| 233 | + |
| 234 | +@_implements(_fft.irfft2) |
| 235 | +def irfft2(a, s=None, axes=(-2, -1), norm=None, workers=None): |
| 236 | + x = _float_utils.__upcast_float16_array(a) |
| 237 | + x = _float_utils.__downcast_float128_array(x) |
| 238 | + return irfftn(x, s, axes, norm, workers) |
| 239 | + |
| 240 | + |
| 241 | +@_implements(_fft.rfftn) |
| 242 | +def rfftn(a, s=None, axes=None, norm=None, workers=None): |
| 243 | + unitary = _unitary(norm) |
| 244 | + x = _float_utils.__upcast_float16_array(a) |
| 245 | + x = _float_utils.__downcast_float128_array(x) |
| 246 | + if unitary: |
| 247 | + x = asarray(x) |
| 248 | + s, axes = _cook_nd_args(x, s, axes) |
| 249 | + with Workers(workers): |
| 250 | + output = _pydfti.rfftn_numpy(x, s, axes) |
| 251 | + if unitary: |
| 252 | + n_tot = prod(asarray(s, dtype=output.dtype)) |
| 253 | + output *= 1 / sqrt(n_tot) |
| 254 | + return output |
| 255 | + |
| 256 | + |
| 257 | +@_implements(_fft.irfftn) |
| 258 | +def irfftn(a, s=None, axes=None, norm=None, workers=None): |
| 259 | + x = _float_utils.__upcast_float16_array(a) |
| 260 | + x = _float_utils.__downcast_float128_array(x) |
| 261 | + with Workers(workers): |
| 262 | + output = _pydfti.irfftn_numpy(x, s, axes) |
| 263 | + if _unitary(norm): |
| 264 | + output *= sqrt(_tot_size(output, axes)) |
| 265 | + return output |
0 commit comments