|
1 | 1 | import numpy as np |
| 2 | +from pyNFFT3.NFMT import BASES |
2 | 3 |
|
3 | 4 | from pyGroupedTransforms import * |
4 | | -from pyNFFT3.NFMT import BASES |
5 | 5 |
|
6 | 6 |
|
7 | 7 | def datalength(bandwidths: np.ndarray) -> int: |
8 | | - """Return the number of mixed-basis Fourier coefficients (zero-frequencies excluded). |
| 8 | + """Return the number of mixed-basis Fourier coefficients (zero-frequencies excluded). |
9 | 9 |
|
10 | | - For each dimension j the frequency range used is: |
11 | | - exp -> {-N_j/2, ..., -1, 1, ..., N_j/2 - 1} (N_j - 1 values) |
12 | | - cos/alg -> {1, 2, ..., N_j - 1} (N_j - 1 values) |
13 | | - Total: prod(N - 1). |
14 | | - """ |
15 | | - if bandwidths.ndim != 1 or bandwidths.dtype != "int32": |
16 | | - return "Please use an one-dimensional numpy.array with dtype 'int32' as input" |
17 | | - if len(bandwidths) == 0: |
18 | | - return 1 |
19 | | - return int(np.prod(bandwidths - 1)) |
| 10 | + For each dimension j the frequency range used is: |
| 11 | + exp -> {-N_j/2, ..., -1, 1, ..., N_j/2 - 1} (N_j - 1 values) |
| 12 | + cos/alg -> {1, 2, ..., N_j - 1} (N_j - 1 values) |
| 13 | + Total: prod(N - 1). |
| 14 | + """ |
| 15 | + if bandwidths.ndim != 1 or bandwidths.dtype != "int32": |
| 16 | + return "Please use an one-dimensional numpy.array with dtype 'int32' as input" |
| 17 | + if len(bandwidths) == 0: |
| 18 | + return 1 |
| 19 | + return int(np.prod(bandwidths - 1)) |
20 | 20 |
|
21 | 21 |
|
22 | 22 | def nfmt_index_set_without_zeros(bandwidths: np.ndarray, bases) -> np.ndarray: |
23 | | - """Return the (d, prod(N-1)) integer frequency matrix. |
24 | | -
|
25 | | - Mirrors Julia's GroupedTransforms.NFMTtools.nfmt_index_set_without_zeros: |
26 | | - exp dimension j: {-N_j/2, ..., -1, 1, ..., N_j/2 - 1} |
27 | | - cos/alg dimension j: {1, 2, ..., N_j - 1} |
28 | | - """ |
29 | | - d = len(bandwidths) |
30 | | - if d == 0: |
31 | | - return np.array([[0]], dtype=np.int64) |
32 | | - |
33 | | - ranges = [] |
34 | | - for n, basis in zip(bandwidths, bases): |
35 | | - if BASES[basis] > 0: |
36 | | - # cos or alg: {1, ..., N-1} |
37 | | - ranges.append(list(range(1, int(n)))) |
38 | | - else: |
39 | | - # exp: {-N/2, ..., -1, 1, ..., N/2-1} |
40 | | - half = int(n) // 2 |
41 | | - ranges.append(list(range(-half, 0)) + list(range(1, half))) |
42 | | - |
43 | | - if d == 1: |
44 | | - return np.array(ranges[0], dtype=np.int64).reshape(1, -1) |
45 | | - |
46 | | - # Cartesian product; last dimension changes fastest (row-major, matching Julia) |
47 | | - mesh = np.array(np.meshgrid(*ranges, indexing="ij"), dtype=np.int64) |
48 | | - return mesh.reshape(d, -1) |
| 23 | + """Return the (d, prod(N-1)) integer frequency matrix. |
| 24 | +
|
| 25 | + Mirrors Julia's GroupedTransforms.NFMTtools.nfmt_index_set_without_zeros: |
| 26 | + exp dimension j: {-N_j/2, ..., -1, 1, ..., N_j/2 - 1} |
| 27 | + cos/alg dimension j: {1, 2, ..., N_j - 1} |
| 28 | + """ |
| 29 | + d = len(bandwidths) |
| 30 | + if d == 0: |
| 31 | + return np.array([[0]], dtype=np.int64) |
| 32 | + |
| 33 | + ranges = [] |
| 34 | + for n, basis in zip(bandwidths, bases): |
| 35 | + if BASES[basis] > 0: |
| 36 | + # cos or alg: {1, ..., N-1} |
| 37 | + ranges.append(list(range(1, int(n)))) |
| 38 | + else: |
| 39 | + # exp: {-N/2, ..., -1, 1, ..., N/2-1} |
| 40 | + half = int(n) // 2 |
| 41 | + ranges.append(list(range(-half, 0)) + list(range(1, half))) |
| 42 | + |
| 43 | + if d == 1: |
| 44 | + return np.array(ranges[0], dtype=np.int64).reshape(1, -1) |
| 45 | + |
| 46 | + # Cartesian product; last dimension changes fastest (row-major, matching Julia) |
| 47 | + mesh = np.array(np.meshgrid(*ranges, indexing="ij"), dtype=np.int64) |
| 48 | + return mesh.reshape(d, -1) |
49 | 49 |
|
50 | 50 |
|
51 | 51 | def get_matrix(bandwidths, X, bases) -> np.ndarray: |
52 | | - """Build the evaluation matrix F with F[m, j] = phi(x_m, k_j). |
53 | | -
|
54 | | - Mirrors Julia's GroupedTransforms.NFMTtools.get_matrix / get_phi: |
55 | | - exp: exp(-2pi i k x) |
56 | | - cos: sqrt(2) * cos(pi * k * x) (k != 0 guaranteed by frequency set) |
57 | | - alg: sqrt(2) * cos(k * arccos(2x - 1)) |
58 | | -
|
59 | | - X is expected in (d, M) format (columns are nodes), consistent with the |
60 | | - convention used in GroupedTransform.get_matrix. |
61 | | - """ |
62 | | - if X.ndim == 1 or X.shape[0] == 1: |
63 | | - X_eval = X.flatten().reshape(-1, 1) # (M, 1) |
64 | | - d = 1 |
65 | | - M = X_eval.shape[0] |
66 | | - else: |
67 | | - d, M = X.shape |
68 | | - X_eval = X.T # (M, d) |
69 | | - |
70 | | - if len(bandwidths) == 0: |
71 | | - return np.ones((M, 1), dtype=np.complex128) |
72 | | - |
73 | | - freq = nfmt_index_set_without_zeros( |
74 | | - np.asarray(bandwidths, dtype=np.int32), list(bases) |
75 | | - ) # (d, nf) |
76 | | - |
77 | | - nf = freq.shape[1] |
78 | | - F = np.ones((M, nf), dtype=np.complex128) |
79 | | - |
80 | | - for j in range(d): |
81 | | - n_j = freq[j] # (nf,) all nonzero by construction |
82 | | - x_j = X_eval[:, j][:, None] # (M, 1) |
83 | | - if BASES[bases[j]] == 1: |
84 | | - # cos basis: sqrt(2) * cos(pi * k * x) |
85 | | - F *= np.sqrt(2.0) * np.cos(np.pi * x_j * n_j) |
86 | | - elif BASES[bases[j]] == 2: |
87 | | - # alg (Chebyshev) basis: sqrt(2) * cos(k * arccos(2x - 1)) |
88 | | - F *= np.sqrt(2.0) * np.cos(n_j * np.arccos(2.0 * x_j - 1.0)) |
89 | | - else: |
90 | | - # exp basis: exp(-2pi i k x) |
91 | | - F *= np.exp(-2.0j * np.pi * x_j * n_j) |
92 | | - |
93 | | - return F |
| 52 | + """Build the evaluation matrix F with F[m, j] = phi(x_m, k_j). |
| 53 | +
|
| 54 | + Mirrors Julia's GroupedTransforms.NFMTtools.get_matrix / get_phi: |
| 55 | + exp: exp(-2pi i k x) |
| 56 | + cos: sqrt(2) * cos(pi * k * x) (k != 0 guaranteed by frequency set) |
| 57 | + alg: sqrt(2) * cos(k * arccos(2x - 1)) |
| 58 | +
|
| 59 | + X is expected in (d, M) format (columns are nodes), consistent with the |
| 60 | + convention used in GroupedTransform.get_matrix. |
| 61 | + """ |
| 62 | + if X.ndim == 1 or X.shape[0] == 1: |
| 63 | + X_eval = X.flatten().reshape(-1, 1) # (M, 1) |
| 64 | + d = 1 |
| 65 | + M = X_eval.shape[0] |
| 66 | + else: |
| 67 | + d, M = X.shape |
| 68 | + X_eval = X.T # (M, d) |
| 69 | + |
| 70 | + if len(bandwidths) == 0: |
| 71 | + return np.ones((M, 1), dtype=np.complex128) |
| 72 | + |
| 73 | + freq = nfmt_index_set_without_zeros( |
| 74 | + np.asarray(bandwidths, dtype=np.int32), list(bases) |
| 75 | + ) # (d, nf) |
| 76 | + |
| 77 | + nf = freq.shape[1] |
| 78 | + F = np.ones((M, nf), dtype=np.complex128) |
| 79 | + |
| 80 | + for j in range(d): |
| 81 | + n_j = freq[j] # (nf,) all nonzero by construction |
| 82 | + x_j = X_eval[:, j][:, None] # (M, 1) |
| 83 | + if BASES[bases[j]] == 1: |
| 84 | + # cos basis: sqrt(2) * cos(pi * k * x) |
| 85 | + F *= np.sqrt(2.0) * np.cos(np.pi * x_j * n_j) |
| 86 | + elif BASES[bases[j]] == 2: |
| 87 | + # alg (Chebyshev) basis: sqrt(2) * cos(k * arccos(2x - 1)) |
| 88 | + F *= np.sqrt(2.0) * np.cos(n_j * np.arccos(2.0 * x_j - 1.0)) |
| 89 | + else: |
| 90 | + # exp basis: exp(-2pi i k x) |
| 91 | + F *= np.exp(-2.0j * np.pi * x_j * n_j) |
| 92 | + |
| 93 | + return F |
94 | 94 |
|
95 | 95 |
|
96 | 96 | def get_transform(bandwidths: np.ndarray, X: np.ndarray, bases): |
97 | | - """Return a DeferredLinearOperator for the mixed-basis grouped transform. |
98 | | -
|
99 | | - Uses the same matrix as get_matrix wrapped in a linear operator so that |
100 | | - trafo and adjoint are consistent by construction. |
101 | | - X is expected in (M, d) format (rows are nodes), consistent with the |
102 | | - convention used in GroupedTransform.transforms. |
103 | | - """ |
104 | | - if bandwidths.ndim > 1 or bandwidths.dtype != "int32": |
105 | | - return "Please use a zero or one-dimensional numpy.array with dtype 'int32' as input" |
106 | | - |
107 | | - M = X.shape[0] |
108 | | - |
109 | | - if len(bandwidths) == 0: |
110 | | - return DeferredLinearOperator( |
111 | | - dtype=np.complex128, |
112 | | - shape=(M, 1), |
113 | | - mfunc=lambda fhat: np.full(M, fhat[0], dtype=np.complex128), |
114 | | - rmfunc=lambda f: np.array([np.sum(f)], dtype=np.complex128), |
115 | | - ) |
116 | | - |
117 | | - # X is (M, d) here; get_matrix expects (d, M) |
118 | | - mat = get_matrix(bandwidths, X.T, list(bases)) # (M, nf) |
119 | | - N = int(np.prod(bandwidths - 1)) |
120 | | - |
121 | | - return DeferredLinearOperator( |
122 | | - dtype=np.complex128, |
123 | | - shape=(M, N), |
124 | | - mfunc=lambda fhat: mat @ np.asarray(fhat, dtype=np.complex128), |
125 | | - rmfunc=lambda f: mat.conj().T @ np.asarray(f, dtype=np.complex128), |
126 | | - ) |
| 97 | + """Return a DeferredLinearOperator for the mixed-basis grouped transform. |
| 98 | +
|
| 99 | + Uses the same matrix as get_matrix wrapped in a linear operator so that |
| 100 | + trafo and adjoint are consistent by construction. |
| 101 | + X is expected in (M, d) format (rows are nodes), consistent with the |
| 102 | + convention used in GroupedTransform.transforms. |
| 103 | + """ |
| 104 | + if bandwidths.ndim > 1 or bandwidths.dtype != "int32": |
| 105 | + return "Please use a zero or one-dimensional numpy.array with dtype 'int32' as input" |
| 106 | + |
| 107 | + M = X.shape[0] |
| 108 | + |
| 109 | + if len(bandwidths) == 0: |
| 110 | + return DeferredLinearOperator( |
| 111 | + dtype=np.complex128, |
| 112 | + shape=(M, 1), |
| 113 | + mfunc=lambda fhat: np.full(M, fhat[0], dtype=np.complex128), |
| 114 | + rmfunc=lambda f: np.array([np.sum(f)], dtype=np.complex128), |
| 115 | + ) |
| 116 | + |
| 117 | + # X is (M, d) here; get_matrix expects (d, M) |
| 118 | + mat = get_matrix(bandwidths, X.T, list(bases)) # (M, nf) |
| 119 | + N = int(np.prod(bandwidths - 1)) |
| 120 | + |
| 121 | + return DeferredLinearOperator( |
| 122 | + dtype=np.complex128, |
| 123 | + shape=(M, N), |
| 124 | + mfunc=lambda fhat: mat @ np.asarray(fhat, dtype=np.complex128), |
| 125 | + rmfunc=lambda f: mat.conj().T @ np.asarray(f, dtype=np.complex128), |
| 126 | + ) |
0 commit comments