diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9b82a1c..71a62a7 100755 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -10,7 +10,7 @@ jobs: test: strategy: matrix: - python-version: ['3.8', '3.9', '3.10', '3.11', '3.12'] # Specify the Python versions you want to test against + python-version: ['3.9', '3.10', '3.11', '3.12'] # Specify the Python versions you want to test against os: [ubuntu-latest, windows-latest] name: Python ${{ matrix.python-version }} - ${{ matrix.os }} runs-on: ${{ matrix.os }} @@ -24,9 +24,9 @@ jobs: - name: Install dependencies run: | python -m pip install --upgrade pip + pip install pyNFFT3 pip install numpy pip install scipy - pip install pyNFFT3 - name: Run tests run: | python -m tests.run_tests diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml new file mode 100644 index 0000000..a5da121 --- /dev/null +++ b/.github/workflows/documentation.yml @@ -0,0 +1,29 @@ +name: documentation + +on: [push, pull_request, workflow_dispatch] + +permissions: + contents: write + +jobs: + docs: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + - name: Install dependencies + run: | + pip install sphinx sphinx_rtd_theme myst_parser + - name: Sphinx build + run: | + cd docs + sphinx-apidoc -o ./source ../src -f + make html + - name: Deploy to GitHub Pages + uses: peaceiris/actions-gh-pages@v3 + if: ${{ github.event_name == 'push' && github.ref == 'refs/heads/main' }} + with: + publish_branch: gh-pages + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: docs/build/html + force_orphan: true diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 5300640..8a51be9 100755 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -1,32 +1,27 @@ -name: Format - -on: - push: - branches: [main] - tags: [v*] - pull_request: - +name: Format Python Code +on: push jobs: - format: - name: "Format Check" + python-code-format: runs-on: ubuntu-latest steps: - - name: Checkout code - uses: actions/checkout@v4 - - name: Set up Python - uses: actions/setup-python@v4 + - uses: actions/checkout@v2 + - uses: actions/setup-python@v4 with: - python-version: '3.x' - - name: Install black + python-version: "3.12" + architecture: "x64" + - name: Display Python version + run: python --version + - name: Install packages + run: pip install black isort + - name: Formatter run: | - python -m pip install --upgrade pip - pip install black - - name: Format code with black - run: | - black . --check - - name: Suggest formatting changes - uses: reviewdog/action-suggester@v1 - if: github.event_name == 'pull_request' + black . + isort . + - name: Create Pull Request + uses: peter-evans/create-pull-request@v3 with: - tool_name: black - fail_on_error: true + commit-message: Auto code format + title: Fixes by format action + body: This is an auto-generated PR with fixes. + labels: automated pr + branch: python-code-format-patches diff --git a/.github/workflows/sphinx.yml b/.github/workflows/sphinx.yml deleted file mode 100755 index 77afde6..0000000 --- a/.github/workflows/sphinx.yml +++ /dev/null @@ -1,24 +0,0 @@ -name: "Sphinx: Render docs" - -on: push - -jobs: - build: - runs-on: ubuntu-latest - permissions: - contents: write - steps: - - uses: actions/checkout@v4 - - name: Build HTML - uses: ammaraskar/sphinx-action@master - - name: Upload artifacts - uses: actions/upload-artifact@v4 - with: - name: html-docs - path: docs/build/html/ - - name: Deploy - uses: peaceiris/actions-gh-pages@v3 - if: github.ref == 'refs/heads/main' - with: - github_token: ${{ secrets.GITHUB_TOKEN }} - publish_dir: docs/build/html diff --git a/docs/requirements.txt b/docs/requirements.txt deleted file mode 100644 index 963391d..0000000 --- a/docs/requirements.txt +++ /dev/null @@ -1,4 +0,0 @@ -furo==2021.11.16 -pyNFFT3==0.2.2 -numpy==1.24.4 -scipy==1.16.1 diff --git a/docs/source/conf.py b/docs/source/conf.py index 8ff42ce..ab16626 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -6,25 +6,22 @@ # -- Project information ----------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information -project = 'GroupedTransform' -copyright = '2025, Felix' -author = 'Felix' -release = '29.07.2025' +project = "GroupedTransform" +copyright = "2025, Felix" +author = "Felix" +release = "29.07.2025" # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration -extensions = [ - 'sphinx.ext.mathjax' -] +extensions = ["sphinx.ext.mathjax"] -templates_path = ['_templates'] +templates_path = ["_templates"] exclude_patterns = [] - # -- Options for HTML output ------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output -html_theme = 'alabaster' -html_static_path = ['_static'] +html_theme = "alabaster" +html_static_path = ["_static"] diff --git a/pyproject.toml b/pyproject.toml index da1c0d2..d96020e 100755 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,7 @@ maintainers = [ ] description = "Python implementation of Grouped Transforms" readme = "README.md" -requires-python = ">=3.8" +requires-python = ">=3.9" dependencies = [ "pyNFFT3>=0.2.2", "numpy>=2.0.0", diff --git a/src/pyGroupedTransforms/CWWTtools.py b/src/pyGroupedTransforms/CWWTtools.py index 0b8e1cf..e440ec2 100644 --- a/src/pyGroupedTransforms/CWWTtools.py +++ b/src/pyGroupedTransforms/CWWTtools.py @@ -1,7 +1,11 @@ import numpy as np + from pyGroupedTransforms import * -def datalength(bandwidths: np.ndarray) -> int: #datalength(bandwidths::Vector{Int})::Int + +def datalength( + bandwidths: np.ndarray, +) -> int: # datalength(bandwidths::Vector{Int})::Int """ `N = datalength(bandwidths)` @@ -14,12 +18,12 @@ def datalength(bandwidths: np.ndarray) -> int: #datalength(bandwidths::Vector{I if bandwidths.size == 0: return 1 elif bandwidths.size == 1: - return 2**(bandwidths[0] + 1) - 1 + return 2 ** (bandwidths[0] + 1) - 1 elif bandwidths.size == 2: - return 2**(bandwidths[0] + 1) * bandwidths[0] + 1 + return 2 ** (bandwidths[0] + 1) * bandwidths[0] + 1 elif bandwidths.size == 3: n = bandwidths[0] - return 2**n * n**2 + 2**n * n + 2**(n + 1) - 1 + return 2**n * n**2 + 2**n * n + 2 ** (n + 1) - 1 else: d = bandwidths.size n = bandwidths[0] @@ -29,7 +33,7 @@ def datalength(bandwidths: np.ndarray) -> int: #datalength(bandwidths::Vector{I return s -def partitions_exact_k(n, k, max_part=None): #helpfunction for cwwt_index_set +def partitions_exact_k(n, k, max_part=None): # helpfunction for cwwt_index_set if max_part is None: max_part = n @@ -41,7 +45,8 @@ def partitions_exact_k(n, k, max_part=None): #helpfunction for cwwt_index_set for tail in partitions_exact_k(n - i, k - 1, i): yield [i] + tail -def partitions(n, k, max_part=None): #helpfunction for cwwt_index_set + +def partitions(n, k, max_part=None): # helpfunction for cwwt_index_set if max_part is None: max_part = n if k == 1: @@ -53,7 +58,7 @@ def partitions(n, k, max_part=None): #helpfunction for cwwt_index_set yield [i] + tail -def cwwt_index_set(n): #cwwt_index_set(n::Vector{Int})::Array{Int} +def cwwt_index_set(n): # cwwt_index_set(n::Vector{Int})::Array{Int} """ `freq = cwwt_index_set(bandwidths)` @@ -70,13 +75,13 @@ def cwwt_index_set(n): #cwwt_index_set(n::Vector{Int})::Array{Int} if d == 1: return np.arange(n[0] + 1).reshape(1, -1) # row vector - freq = np.zeros((d, 1), dtype=int) + freq = np.zeros((d, 1), dtype=int) for j in range(1, n[-1] + 1): for x in partitions(d + j, d): x = [xi - 1 for xi in x] if all(xi <= j for xi in x): - ys = set(permutations(x)) #diffrent Order as in julia + ys = set(permutations(x)) # diffrent Order as in julia for y in ys: y_col = np.array(y).reshape(d, 1) freq = np.hstack((freq, y_col)) @@ -84,7 +89,7 @@ def cwwt_index_set(n): #cwwt_index_set(n::Vector{Int})::Array{Int} return freq -def begin_index2d(j): #begin_index2d(j::Int)::Int +def begin_index2d(j): # begin_index2d(j::Int)::Int """ Begin of 2-dimensional index sets """ @@ -93,10 +98,11 @@ def begin_index2d(j): #begin_index2d(j::Int)::Int elif j == 1: ind = 2 else: - ind = (2**j)*(j-1)+2 + ind = (2**j) * (j - 1) + 2 return ind -def indextoN(j,k): #indextoN(j::Array{Int}, k::Array{Int})::Array{Int} + +def indextoN(j, k): # indextoN(j::Array{Int}, k::Array{Int})::Array{Int} """ # function from index set of wavelets to natural numbers, i.e. # (j,k) maps to N @@ -107,32 +113,32 @@ def indextoN(j,k): #indextoN(j::Array{Int}, k::Array{Int})::Array{Int} # out in N # creates row vector with entry for every column vector in k """ - d = len(j) #dimension + d = len(j) # dimension d2 = k.shape[0] - s = k.shape[1] + s = k.shape[1] if d != d2: ValueError("j and k have to have same length, k has to be column vector.") out = np.zeros(s) - if d == 1: for i in range(s): out[i] = 2 ** (j[0]) + k[0][i] elif d == 2: level = np.sum(j) - for i in range(s): #in julia j[0] ist j(1) ? kommt auch nen fehler raus.. - out[i] = begin_index2d(level) + j[0]*2**level + 2 **(j[1]) + k[0, i] + k[1, i] - + for i in range(s): # in julia j[0] ist j(1) ? kommt auch nen fehler raus.. + out[i] = ( + begin_index2d(level) + j[0] * 2**level + 2 ** (j[1]) + k[0, i] + k[1, i] + ) return out -#Chui-Wang-Wavelet Function for different orders m : -def Chui_wavelet(x,m): #Chui_wavelet(x::Array{Float64}, m::Int)::Array{Float64} +# Chui-Wang-Wavelet Function for different orders m : +def Chui_wavelet(x, m): # Chui_wavelet(x::Array{Float64}, m::Int)::Array{Float64} """ % periodic Chui-Wang-Wavelets, % Chui Wang has support [0,2m-1], @@ -160,7 +166,8 @@ def Chui_wavelet(x,m): #Chui_wavelet(x::Array{Float64}, m::Int)::Array{Float64 return psi -#"periodic 1-d-wavelet for one j and different k:" + +# "periodic 1-d-wavelet for one j and different k:" def _Chui_periodic_1d(x, m, j, k): """ 1D periodic Chui-Wang wavelet @@ -185,17 +192,18 @@ def _Chui_periodic_1d(x, m, j, k): for i in range(mm): l = np.ceil(k[:, i] / 2**j - x) arg = 2**j * (x + l) - k[:, i] - y[:, i] = 2**(j / 2) * Chui_wavelet(arg.flatten(), m) + y[:, i] = 2 ** (j / 2) * Chui_wavelet(arg.flatten(), m) else: for i in range(mm): for ll in range(2 * m - 1): l = np.ceil(k[:, i] / 2**j - x) + ll arg = 2**j * (x + l) - k[:, i] - y[:, i] += 2**(j / 2) * Chui_wavelet(arg.flatten(), m) + y[:, i] += 2 ** (j / 2) * Chui_wavelet(arg.flatten(), m) return y -#periodic muliti-d-wavelet : + +# periodic muliti-d-wavelet : def Chui_periodic(x, m, j, k): """ For d == 1 use _Chui_periodic_1d, else for d > 1: @@ -226,14 +234,15 @@ def Chui_periodic(x, m, j, k): return y -def get_transform(bandwidths, X, m): #get_transform(bandwidths::Vector{Int}, X::Array{Float64}, m::Int)::LinearMap - +def get_transform( + bandwidths, X, m +): # get_transform(bandwidths::Vector{Int}, X::Array{Float64}, m::Int)::LinearMap """ `F = get_transform(bandwidths, X, order) # Input: * `bandwidths::Vector{Int}` - * `X::Array{Float64}` ... nodes in M x |u| format + * `X::Array{Float64}` ... nodes in M x |u| format # Output: * `F::LinearMap{ComplexF64}` ... Linear maps of the sparse Matrices @@ -244,7 +253,7 @@ def get_transform(bandwidths, X, m): #get_transform(bandwidths::Vector{Int}, X d = 1 M = len(X) else: - M, d = X.shape + M, d = X.shape freq = cwwt_index_set(bandwidths) @@ -262,7 +271,9 @@ def trafo(x): def adjoint(x): return B @ x - return DeferredLinearOperator(dtype=np.float64, shape=(M, int(max(J) + 1)), mfunc=trafo, rmfunc=adjoint) + return DeferredLinearOperator( + dtype=np.float64, shape=(M, int(max(J) + 1)), mfunc=trafo, rmfunc=adjoint + ) if d == 1: X_col = X if X.ndim == 1 else X[:, 0] @@ -272,14 +283,14 @@ def adjoint(x): V = Chui_periodic(X_col, m, 0, np.array([[0]])).flatten() for j in range(1, bandwidths[0] + 1): - num = min(2 ** j, 2 * m - 1) + num = min(2**j, 2 * m - 1) - a = (np.floor(2 ** j * X_col)[:, None] - 2 * m + 2) * np.ones((1, num)) + a = (np.floor(2**j * X_col)[:, None] - 2 * m + 2) * np.ones((1, num)) b = np.ones((M, 1)) * (np.arange(num).T) - k = np.mod(a + b, 2 ** j) + k = np.mod(a + b, 2**j) I_new = np.repeat(np.arange(M), num) - J_new = (k + 2 ** j).flatten() - 1 + J_new = (k + 2**j).flatten() - 1 V_new = Chui_periodic(X_col, m, j, k).flatten() I = np.concatenate((I, I_new)) @@ -297,7 +308,9 @@ def trafo(x): def adjoint(x): return B @ x - return DeferredLinearOperator(dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint) + return DeferredLinearOperator( + dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint + ) elif d == 2: X1 = X[:, 0] @@ -307,17 +320,23 @@ def adjoint(x): J = np.zeros(M) V = Chui_periodic(X, m, [0, 0], np.array([[[0, 0]]])).flatten() - freq = freq[:, 1:freq.shape[1]] + freq = freq[:, 1 : freq.shape[1]] ac_co = 2 for j in freq.T: num1 = min(2 ** j[0], 2 * m - 1) num2 = min(2 ** j[1], 2 * m - 1) - k1 = np.mod((np.floor(2 ** j[0] * X1)[:, None] - 2 * m + 2) * np.ones((1, num1)) + - np.ones((M, 1)) * np.arange(num1).T, 2 ** j[0]) - k2 = np.mod((np.floor(2 ** j[1] * X2)[:, None] - 2 * m + 2) * np.ones((1, num2)) + - np.ones((M, 1)) * np.arange(num2).T, 2 ** j[1]) + k1 = np.mod( + (np.floor(2 ** j[0] * X1)[:, None] - 2 * m + 2) * np.ones((1, num1)) + + np.ones((M, 1)) * np.arange(num1).T, + 2 ** j[0], + ) + k2 = np.mod( + (np.floor(2 ** j[1] * X2)[:, None] - 2 * m + 2) * np.ones((1, num2)) + + np.ones((M, 1)) * np.arange(num2).T, + 2 ** j[1], + ) if len(k1.shape) == 1: k1 = k1[:, None] @@ -328,10 +347,16 @@ def adjoint(x): k = np.zeros((M, num1 * num2, d)) for kk1 in range(num1): - k_ind[:, kk1 * num2:(kk1 + 1) * num2] = 2 ** j[1] * k1[:, kk1][:, None] + k2 - k[:, kk1 * num2:(kk1 + 1) * num2, :] = np.concatenate( - (np.expand_dims(k1[:, kk1][:, None] * np.ones(num2), axis=2), - np.expand_dims(k2, axis=2)), axis=2) + k_ind[:, kk1 * num2 : (kk1 + 1) * num2] = ( + 2 ** j[1] * k1[:, kk1][:, None] + k2 + ) + k[:, kk1 * num2 : (kk1 + 1) * num2, :] = np.concatenate( + ( + np.expand_dims(k1[:, kk1][:, None] * np.ones(num2), axis=2), + np.expand_dims(k2, axis=2), + ), + axis=2, + ) I_new = np.repeat(np.arange(M), num1 * num2) J_new = (k_ind + ac_co).flatten() - 1 @@ -341,7 +366,7 @@ def adjoint(x): J = np.concatenate((J, J_new)) V = np.concatenate((V, V_new)) - ac_co = ac_co + np.sum(2 ** j) + ac_co = ac_co + np.sum(2**j) A = coo_matrix((V, (I.astype(int), J.astype(int))), (M, int(max(J) + 1))) A = A.tocsc() @@ -354,7 +379,9 @@ def trafo(x): def adjoint(x): return B @ x - return DeferredLinearOperator(dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint) + return DeferredLinearOperator( + dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint + ) elif d == 3: I = np.arange(M) @@ -369,9 +396,18 @@ def adjoint(x): num2 = min(2 ** j[1], 2 * m - 1) num3 = min(2 ** j[2], 2 * m - 1) - k1 = np.mod((np.floor(2 ** j[0] * X[:, 0])[:, None] - 2 * m + 2) + np.arange(num1), 2 ** j[0]) - k2 = np.mod((np.floor(2 ** j[1] * X[:, 1])[:, None] - 2 * m + 2) + np.arange(num2), 2 ** j[1]) - k3 = np.mod((np.floor(2 ** j[2] * X[:, 2])[:, None] - 2 * m + 2) + np.arange(num3), 2 ** j[2]) + k1 = np.mod( + (np.floor(2 ** j[0] * X[:, 0])[:, None] - 2 * m + 2) + np.arange(num1), + 2 ** j[0], + ) + k2 = np.mod( + (np.floor(2 ** j[1] * X[:, 1])[:, None] - 2 * m + 2) + np.arange(num2), + 2 ** j[1], + ) + k3 = np.mod( + (np.floor(2 ** j[2] * X[:, 2])[:, None] - 2 * m + 2) + np.arange(num3), + 2 ** j[2], + ) k_ind = np.zeros((M, num1 * num2 * num3)) k = np.zeros((M, num1 * num2 * num3, d)) @@ -410,15 +446,16 @@ def trafo(x): def adjoint(x): return B * x - return DeferredLinearOperator(dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint) - + return DeferredLinearOperator( + dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint + ) elif d == 4: I = np.arange(M) J = np.zeros(M) V = Chui_periodic(X, m, [0, 0, 0, 0], np.array([[[0, 0, 0, 0]]])).flatten() - freq = freq[:, 1:freq.shape[1]] + freq = freq[:, 1 : freq.shape[1]] ac_co = 2 for j in freq.T: @@ -427,14 +464,30 @@ def adjoint(x): num3 = min(2 ** j[2], 2 * m - 1) num2 = min(2 ** j[3], 2 * m - 1) - k1 = np.mod((np.floor(2 ** j[0] * X[:,0])[:, None] - 2 * m + 2) * np.ones((1, num1)) + - np.ones((M, 1)) * np.arange(num1), 2 ** j[0]) - k2 = np.mod((np.floor(2 ** j[1] * X[:,1])[:, None] - 2 * m + 2) * np.ones((1, num2)) + - np.ones((M, 1)) * np.arange(num2), 2 ** j[1]) - k3 = np.mod((np.floor(2 ** j[2] * X[:,2])[:, None] - 2 * m + 2) * np.ones((1, num3)) + - np.ones((M, 1)) * np.arange(num3), 2 ** j[2]) - k4 = np.mod((np.floor(2 ** j[3] * X[:,3])[:, None] - 2 * m + 2) * np.ones((1, num4)) + - np.ones((M, 1)) * np.arange(num4), 2 ** j[3]) + k1 = np.mod( + (np.floor(2 ** j[0] * X[:, 0])[:, None] - 2 * m + 2) + * np.ones((1, num1)) + + np.ones((M, 1)) * np.arange(num1), + 2 ** j[0], + ) + k2 = np.mod( + (np.floor(2 ** j[1] * X[:, 1])[:, None] - 2 * m + 2) + * np.ones((1, num2)) + + np.ones((M, 1)) * np.arange(num2), + 2 ** j[1], + ) + k3 = np.mod( + (np.floor(2 ** j[2] * X[:, 2])[:, None] - 2 * m + 2) + * np.ones((1, num3)) + + np.ones((M, 1)) * np.arange(num3), + 2 ** j[2], + ) + k4 = np.mod( + (np.floor(2 ** j[3] * X[:, 3])[:, None] - 2 * m + 2) + * np.ones((1, num4)) + + np.ones((M, 1)) * np.arange(num4), + 2 ** j[3], + ) total_size = num1 * num2 * num3 * num4 k_ind = np.zeros((M, total_size)) @@ -445,10 +498,12 @@ def adjoint(x): for kk2 in range(num2): for kk3 in range(num3): for kk4 in range(num4): - k_ind[:, idx] = (2 ** (j[1] + j[2] + j[3]) * k1[:, kk1] + - 2 ** (j[2] + j[3]) * k2[:, kk2] + - 2 ** j[3] * k3[:, kk3] + - k4[:, kk4]) + k_ind[:, idx] = ( + 2 ** (j[1] + j[2] + j[3]) * k1[:, kk1] + + 2 ** (j[2] + j[3]) * k2[:, kk2] + + 2 ** j[3] * k3[:, kk3] + + k4[:, kk4] + ) k[:, idx, 0] = k1[:, kk1] k[:, idx, 1] = k2[:, kk2] k[:, idx, 2] = k3[:, kk3] @@ -465,25 +520,29 @@ def adjoint(x): ac_co += 2 ** sum(j) - A = coo_matrix((V, (I.astype(int), J.astype(int))),(M, int(max(J) +1 ))) + A = coo_matrix((V, (I.astype(int), J.astype(int))), (M, int(max(J) + 1))) A = A.tocsc() A = A[:, A.getnnz(0) > 0] B = A.T.conj() def trafo(x): - return A*x + return A * x def adjoint(x): return B * x - return DeferredLinearOperator(dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint) + return DeferredLinearOperator( + dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint + ) elif d == 5: I = np.arange(M) J = np.zeros(M) - V = Chui_periodic(X, m, [0, 0, 0, 0, 0], np.array([[[0, 0, 0, 0, 0]]])).flatten() + V = Chui_periodic( + X, m, [0, 0, 0, 0, 0], np.array([[[0, 0, 0, 0, 0]]]) + ).flatten() - freq = freq[:, 1:freq.shape[1]] + freq = freq[:, 1 : freq.shape[1]] ac_co = 2 for j in freq.T: @@ -493,16 +552,36 @@ def adjoint(x): num4 = min(2 ** j[3], 2 * m - 1) num5 = min(2 ** j[4], 2 * m - 1) - k1 = np.mod((np.floor(2 ** j[0] * X[:,0])[:, None] - 2 * m + 2) * np.ones((1, num1)) + - np.ones((M, 1)) * np.arange(num1), 2 ** j[0]) - k2 = np.mod((np.floor(2 ** j[1] * X[:,1])[:, None] - 2 * m + 2) * np.ones((1, num2)) + - np.ones((M, 1)) * np.arange(num2), 2 ** j[1]) - k3 = np.mod((np.floor(2 ** j[2] * X[:,2])[:, None] - 2 * m + 2) * np.ones((1, num3)) + - np.ones((M, 1)) * np.arange(num3), 2 ** j[2]) - k4 = np.mod((np.floor(2 ** j[3] * X[:,3])[:, None] - 2 * m + 2) * np.ones((1, num4)) + - np.ones((M, 1)) * np.arange(num4), 2 ** j[3]) - k5 = np.mod((np.floor(2 ** j[4] * X[:,4])[:, None] - 2 * m + 2) * np.ones((1, num5)) + - np.ones((M, 1)) * np.arange(num5), 2 ** j[4]) + k1 = np.mod( + (np.floor(2 ** j[0] * X[:, 0])[:, None] - 2 * m + 2) + * np.ones((1, num1)) + + np.ones((M, 1)) * np.arange(num1), + 2 ** j[0], + ) + k2 = np.mod( + (np.floor(2 ** j[1] * X[:, 1])[:, None] - 2 * m + 2) + * np.ones((1, num2)) + + np.ones((M, 1)) * np.arange(num2), + 2 ** j[1], + ) + k3 = np.mod( + (np.floor(2 ** j[2] * X[:, 2])[:, None] - 2 * m + 2) + * np.ones((1, num3)) + + np.ones((M, 1)) * np.arange(num3), + 2 ** j[2], + ) + k4 = np.mod( + (np.floor(2 ** j[3] * X[:, 3])[:, None] - 2 * m + 2) + * np.ones((1, num4)) + + np.ones((M, 1)) * np.arange(num4), + 2 ** j[3], + ) + k5 = np.mod( + (np.floor(2 ** j[4] * X[:, 4])[:, None] - 2 * m + 2) + * np.ones((1, num5)) + + np.ones((M, 1)) * np.arange(num5), + 2 ** j[4], + ) total_size = num1 * num2 * num3 * num4 * num5 k_ind = np.zeros((M, total_size)) @@ -514,11 +593,13 @@ def adjoint(x): for kk3 in range(num3): for kk4 in range(num4): for kk5 in range(num5): - k_ind[:, idx] = (2 ** (j[1] + j[2] + j[3] + j[4]) * k1[:, kk1] + - 2 ** (j[2] + j[3] + j[4]) * k2[:, kk2] + - 2 ** (j[3] + j[4]) * k3[:, kk3] + - 2 ** j[4] * k4[:, kk4] + - k5[:, kk5]) + k_ind[:, idx] = ( + 2 ** (j[1] + j[2] + j[3] + j[4]) * k1[:, kk1] + + 2 ** (j[2] + j[3] + j[4]) * k2[:, kk2] + + 2 ** (j[3] + j[4]) * k3[:, kk3] + + 2 ** j[4] * k4[:, kk4] + + k5[:, kk5] + ) k[:, idx, 0] = k1[:, kk1] k[:, idx, 1] = k2[:, kk2] k[:, idx, 2] = k3[:, kk3] @@ -536,16 +617,17 @@ def adjoint(x): ac_co += 2 ** sum(j) - A = coo_matrix((V, (I.astype(int), J.astype(int))),(M, int(max(J) +1 ))) + A = coo_matrix((V, (I.astype(int), J.astype(int))), (M, int(max(J) + 1))) A = A.tocsc() A = A[:, A.getnnz(0) > 0] B = A.T.conj() def trafo(x): - return A*x + return A * x def adjoint(x): return B * x - return DeferredLinearOperator(dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint) - + return DeferredLinearOperator( + dtype=np.float64, shape=(M, A.shape[1]), mfunc=trafo, rmfunc=adjoint + ) diff --git a/src/pyGroupedTransforms/DeferredLinearOperator.py b/src/pyGroupedTransforms/DeferredLinearOperator.py index 7192bce..3d46c82 100644 --- a/src/pyGroupedTransforms/DeferredLinearOperator.py +++ b/src/pyGroupedTransforms/DeferredLinearOperator.py @@ -1,12 +1,19 @@ from scipy.sparse.linalg import LinearOperator + class DeferredLinearOperator(LinearOperator): - def __init__(self, dtype = 'complex', shape = None, mfunc = lambda x: "not yet initialized", rmfunc = lambda x: "not yet initialized"): + def __init__( + self, + dtype="complex", + shape=None, + mfunc=lambda x: "not yet initialized", + rmfunc=lambda x: "not yet initialized", + ): self.sh = False if shape == None: - shape = (0,0) + shape = (0, 0) else: self.sh = True @@ -14,9 +21,9 @@ def __init__(self, dtype = 'complex', shape = None, mfunc = lambda x: "not yet i self._rmatvec_func = rmfunc self._matvec = lambda x: self._matvec_func(x) self._rmatvec = lambda x: self._rmatvec_func(x) - super().__init__(dtype=dtype, shape=shape) + super().__init__(dtype=dtype, shape=shape) - def set(self, dtype = None, shape = None, mfunc = None, rmfunc = None): + def set(self, dtype=None, shape=None, mfunc=None, rmfunc=None): if dtype != None: self.dtype = dtype @@ -34,7 +41,7 @@ def _matvec(self, x): if self._matvec_func(0) == "not implemented yet" or self.sh == False: raise NotImplementedError("LinearOperator is not yet initialized.") - print(self.LO@x) + print(self.LO @ x) return self @ x def _rmatvec(self, x): @@ -44,4 +51,3 @@ def _rmatvec(self, x): raise NotImplementedError("rmatvec is not yet implemented") return self.H @ x - diff --git a/src/pyGroupedTransforms/GroupedCoefficients.py b/src/pyGroupedTransforms/GroupedCoefficients.py index 1cecf13..9a4642e 100644 --- a/src/pyGroupedTransforms/GroupedCoefficients.py +++ b/src/pyGroupedTransforms/GroupedCoefficients.py @@ -1,12 +1,11 @@ -import numpy as np -from pyGroupedTransforms import NFFTtools -from pyGroupedTransforms import NFCTtools -from pyGroupedTransforms import NFMTtools -from pyGroupedTransforms import CWWTtools +import numpy as np -class GC: # Superclass of GroupedCoefficientsComplex and GroupedCoefficientsReal. +from pyGroupedTransforms import CWWTtools, NFCTtools, NFFTtools, NFMTtools - def __getitem__(self, idx): + +class GC: # Superclass of GroupedCoefficientsComplex and GroupedCoefficientsReal. + + def __getitem__(self, idx): """ (fhat::GroupedCoefficients[u::Vector{Int}]) @@ -27,13 +26,14 @@ def __getitem__(self, idx): if len(s.u) == len(idx) and s.u == idx: stop = start + s.mode.datalength(s.bandwidths) return np.array(self.data[start:stop]) - else: + else: start += s.mode.datalength(s.bandwidths) raise KeyError("This term is not contained") else: - raise TypeError("Index must be an integer or a tuple that contains integer or the empty tuple") - + raise TypeError( + "Index must be an integer or a tuple that contains integer or the empty tuple" + ) def __setitem__(self, idx, Number): """ @@ -53,7 +53,12 @@ def __setitem__(self, idx, Number): elif type(idx) == tuple and (len(idx) == 0 or type(idx[0]) == int): if not isinstance(Number, np.ndarray): raise TypeError("Type mismatch.") - if isinstance(self, GroupedCoefficientsComplex) and np.isrealobj(Number) or isinstance(self, GroupedCoefficientsReal) and np.iscomplexobj(Number): + if ( + isinstance(self, GroupedCoefficientsComplex) + and np.isrealobj(Number) + or isinstance(self, GroupedCoefficientsReal) + and np.iscomplexobj(Number) + ): raise TypeError("Type mismatch.") start = 0 @@ -62,7 +67,7 @@ def __setitem__(self, idx, Number): stop = start + s.mode.datalength(s.bandwidths) self.data[start:stop] = Number return None - else: + else: start += s.mode.datalength(s.bandwidths) raise KeyError("This term is not contained") @@ -70,16 +75,14 @@ def __setitem__(self, idx, Number): else: raise TypeError("Index must be an int or a tuple that contains integer") - - def vec(self): + def vec(self): """ (vec( fhat::GroupedCoefficients )::Vector{<:Number}) - This function returns the vector of the basis coefficients of self. + This function returns the vector of the basis coefficients of self. """ return self.data - def __rmul__(self, alpha): """ (*( z::Number, fhat::GroupedCoefficients )::GroupedCoefficients) @@ -90,15 +93,13 @@ def __rmul__(self, alpha): return GroupedCoefficients(self.settings, alpha * self.data) return NotImplemented - - def __mul__(self, alpha): + def __mul__(self, alpha): if isinstance(alpha, (int, float, complex)): return GroupedCoefficients(self.settings, alpha * self.data) return NotImplemented - - def __add__(self, other): + def __add__(self, other): """ (+( z::Number, fhat::GroupedCoefficients )::GroupedCoefficients) @@ -112,7 +113,6 @@ def __add__(self, other): else: raise ValueError("Settings mismatch.") - def __sub__(self, other): """ (-( z::Number, fhat::GroupedCoefficients )::GroupedCoefficients) @@ -121,7 +121,6 @@ def __sub__(self, other): """ return self + (-1 * other) - def set_data(self, data): """ (set_data!( @@ -130,22 +129,26 @@ def set_data(self, data): With this function one can set the data of a GC object. """ - if isinstance(self, GroupedCoefficientsComplex) and np.isrealobj(data) or isinstance(self, GroupedCoefficientsReal) and np.iscomplexobj(data): - raise TypeError("Type mismatch.") - else: + if ( + isinstance(self, GroupedCoefficientsComplex) + and np.isrealobj(data) + or isinstance(self, GroupedCoefficientsReal) + and np.iscomplexobj(data) + ): + raise TypeError("Type mismatch.") + else: self.data[:] = data return None - def norms(self, Dict=False, other = None, m = None): - + def norms(self, Dict=False, other=None, m=None): """ If other == None: (norms(fhat::GroupedCoefficients; dict =false)) - If other != None and m == None: + If other != None and m == None: (norms(fhat::GroupedCoefficients, what::GroupedCoefficients)::Vector{Float64}) - If other == None and m != None: + If other == None and m != None: (function norms(fhat::GroupedCoefficients, m::Int ; dict =false)::Union{Vector{Float64},Dict{Vector{Int},Float64}}) Inputs sind vertauscht leider... """ @@ -153,64 +156,72 @@ def norms(self, Dict=False, other = None, m = None): if other == None: if m == None: if Dict == False: - return [np.linalg.norm(self[self.settings[i].u]) - for i in range(len(self.settings))] + return [ + np.linalg.norm(self[self.settings[i].u]) + for i in range(len(self.settings)) + ] else: dd = {} for s in self.settings: - if len(s.u) > 0: - dd[tuple(s.u)] = np.linalg.norm(self[s.u]) + if len(s.u) > 0: + dd[tuple(s.u)] = np.linalg.norm(self[s.u]) return dd - else: #m != None + else: # m != None if Dict == False: - if self.settings[0].mode == CWWTtools: #does not work (in julia and in python... ?) + if ( + self.settings[0].mode == CWWTtools + ): # does not work (in julia and in python... ?) n = np.zeros(len(self.settings)) - for i in range(len(self.settings)): + for i in range(len(self.settings)): s = self.settings[i] d = len(self.settings[i].u) - ac_in = 1 + ac_in = 1 freq = CWWTtools.cwwt_index_set(self.settings[i].bandwidths) - if d ==1: + if d == 1: freq = freq.T for jj in range(freq.shape[1]): - j = freq[:,jj] - a = self[s.u][ac_in -1 :ac_in+2**np.sum(j)-1] - Psi = circulant(variances(j[0],m)) + j = freq[:, jj] + a = self[s.u][ac_in - 1 : ac_in + 2 ** np.sum(j) - 1] + Psi = circulant(variances(j[0], m)) if d > 1: - for dd in range(1,d): - Psi = np.kron(Psi, circulant(variances(j[dd],m))) - - n[i] += a @ Psi @ a - ac_in = ac_in +2**np.sum(j) + for dd in range(1, d): + Psi = np.kron( + Psi, circulant(variances(j[dd], m)) + ) + + n[i] += a @ Psi @ a + ac_in = ac_in + 2 ** np.sum(j) return np.sqrt(n) else: return [np.linalg.norm(self[s.u]) for s in self.settings] - else: #dict == True + else: # dict == True dd = {} - if self.settings[0].mode == CWWTtools: #doesn't work... ? + if self.settings[0].mode == CWWTtools: # doesn't work... ? n = np.zeros(len(self.settings)) - for i in range(len(self.settings)): + for i in range(len(self.settings)): s = self.settings[i] d = len(self.settings[i].u) - ac_in = 1 + ac_in = 1 freq = CWWTtools.cwwt_index_set(self.settings[i].bandwidths) - if d ==1: + if d == 1: freq = freq.T for jj in range(freq.shape[1]): - j = freq[:,jj] - a = self[s.u][ac_in -1 :ac_in+2 **np.sum(j)-1] - Psi = circulant(variances(j[0],m)) + j = freq[:, jj] + a = self[s.u][ac_in - 1 : ac_in + 2 ** np.sum(j) - 1] + Psi = circulant(variances(j[0], m)) if d > 1: - for dd in range(1,d): - Psi = np.kron(Psi, circulant(variances(j[dd],m))) + for dd in range(1, d): + Psi = np.kron( + Psi, circulant(variances(j[dd], m)) + ) - n[i] += a @ Psi @ a - ac_in = ac_in +2**np.sum(j) + n[i] += a @ Psi @ a + ac_in = ac_in + 2 ** np.sum(j) if len(self.settings[i].u) != 0: dd[self.settings[i].u] = np.sqrt(n[i]) @@ -218,39 +229,45 @@ def norms(self, Dict=False, other = None, m = None): else: for i in range(len(self.settings)): if len(self.settings[i].u) != 0: - dd[tuple(self.settings[i].u)] = np.linalg.norm(self[self.settings[i].u]) + dd[tuple(self.settings[i].u)] = np.linalg.norm( + self[self.settings[i].u] + ) return dd - else: #other != None - c = GroupedCoefficients(self.settings, (np.sqrt(np.real(other.data))) * self.data) + else: # other != None + c = GroupedCoefficients( + self.settings, (np.sqrt(np.real(other.data))) * self.data + ) return [np.linalg.norm(c[c.settings[i].u]) for i in range(len(c.settings))] - - - -#Matrix of variances between two basis functions, needed for wavelet basis, since they are not orthonormal -def variances(j,m): +# Matrix of variances between two basis functions, needed for wavelet basis, since they are not orthonormal +def variances(j, m): """ - (variances(j::Int,m::Int)::Vector{Float64}) + (variances(j::Int,m::Int)::Vector{Float64}) - INPUT - j ... level of wavelet - m ... order of wavelet - OUTPUT - y = ()for k = 0...2^j-1 - (psi_{j,k}) are the wavelets, output contains a vector of all scalar products of one level - for 2^j >2m*1 the same values, but more zeros for higher j. + INPUT + j ... level of wavelet + m ... order of wavelet + OUTPUT + y = ()for k = 0...2^j-1 + (psi_{j,k}) are the wavelets, output contains a vector of all scalar products of one level + for 2^j >2m*1 the same values, but more zeros for higher j. """ if m == 2: if j == 0: - y = [1/3] + y = [1 / 3] elif j == 1: y = [0.240740740740715, 0.092592592592576] elif j == 2: - y = [0.250000000000097, 0.046296296296494, -0.009259259259238, 0.046296296296451] + y = [ + 0.250000000000097, + 0.046296296296494, + -0.009259259259238, + 0.046296296296451, + ] else: y = np.zeros(2**j) y[0:3] = [0.249999999999968, 0.046296296296379, -0.004629629629584] @@ -262,15 +279,32 @@ def variances(j,m): elif j == 1: y = [0.085629629629627, 0.047703703703701] elif j == 2: - y = [0.098444444444468, 0.023851851851835, -0.012814814814815, 0.023851851851891] + y = [ + 0.098444444444468, + 0.023851851851835, + -0.012814814814815, + 0.023851851851891, + ] elif j == 3: - y = [0.098443287037094, 0.024151620370357, -0.006407407407473, - -2.997685185108751e-04, 1.157407410282006e-06, - -2.997685185108751e-04, -0.006407407407473, 0.024151620370] + y = [ + 0.098443287037094, + 0.024151620370357, + -0.006407407407473, + -2.997685185108751e-04, + 1.157407410282006e-06, + -2.997685185108751e-04, + -0.006407407407473, + 0.024151620370, + ] else: y = np.zeros(2**j) - y[0:5] = [0.098443287037052, 0.024151620370458, -0.006407407407375, - -2.997685186147822e-04, 5.787036288131668e-07] + y[0:5] = [ + 0.098443287037052, + 0.024151620370458, + -0.006407407407375, + -2.997685186147822e-04, + 5.787036288131668e-07, + ] y[-5:-1] = y[1:5] elif m == 4: @@ -279,26 +313,41 @@ def variances(j,m): elif j == 1: y = [0.032014093350, 0.021954160617785] elif j == 2: - y = [0.041543521817902, 0.010977080308839, - -0.009529428467484, 0.010977080308902] + y = [ + 0.041543521817902, + 0.010977080308839, + -0.009529428467484, + 0.010977080308902, + ] elif j == 3: - y = [0.041534149423629, 0.011629855618373, - -0.004764714233735, -6.527753094869201e-04, - 9.372394228547530e-06, -6.527753094788768e-04, - -0.004764714233719, 0.011629855618379] + y = [ + 0.041534149423629, + 0.011629855618373, + -0.004764714233735, + -6.527753094869201e-04, + 9.372394228547530e-06, + -6.527753094788768e-04, + -0.004764714233719, + 0.011629855618379, + ] else: y = np.zeros(2**j) - y[0:7] = [0.041534149423653, 0.011629855618368, -0.004764714225921, - -6.528682451447852e-04, 4.686197131487609e-06, - 9.293564145017154e-08, -7.821712475559538e-12] + y[0:7] = [ + 0.041534149423653, + 0.011629855618368, + -0.004764714225921, + -6.528682451447852e-04, + 4.686197131487609e-06, + 9.293564145017154e-08, + -7.821712475559538e-12, + ] y[-7:-1] = y[1:7] return y - class Setting: - def __init__(self, u = None, mode = NFFTtools, bandwidths = None, bases = None): + def __init__(self, u=None, mode=NFFTtools, bandwidths=None, bases=None): self.u = u self.mode = mode self.bandwidths = bandwidths @@ -306,25 +355,27 @@ def __init__(self, u = None, mode = NFFTtools, bandwidths = None, bases = None): def __eq__(self, other): if isinstance(other, Setting): - boollist = [False,False,False,False] + boollist = [False, False, False, False] if self.mode == other.mode: boollist[0] = True if len(self.u) == 0 == len(other.u): boollist[1] = True elif len(self.u) == len(other.u) and self.u == other.u: - boollist[1] = True + boollist[1] = True if len(self.bandwidths) == 0 == len(other.bandwidths): boollist[2] = True - elif len(self.bandwidths) == len(other.bandwidths) and (self.bandwidths == other.bandwidths).all(): - boollist[2] = True - - + elif ( + len(self.bandwidths) == len(other.bandwidths) + and (self.bandwidths == other.bandwidths).all() + ): + boollist[2] = True + if self.bases == None == other.bases: boollist[3] = True elif len(self.bases) == len(other.bases) and self.bases == other.bases: - boollist[3] = True + boollist[3] = True return boollist[0] and boollist[2] and boollist[2] and boollist[3] else: @@ -334,8 +385,8 @@ def __eq__(self, other): # A class to hold complex coefficients belonging to indices in a grouped index set ### Fields: -#* `setting::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}` - uniquely describes the setting such as the bandlimits ``N_{\pmb u}``, see also [`get_setting(system::String,d::Int,ds::Int,N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) and #[`get_setting(system::String,U::Vector{Vector{Int}},N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) -#* `data::Union{Vector{ComplexF64},Nothing}` - the vector of coefficients +# * `setting::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}` - uniquely describes the setting such as the bandlimits ``N_{\pmb u}``, see also [`get_setting(system::String,d::Int,ds::Int,N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) and #[`get_setting(system::String,U::Vector{Vector{Int}},N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) +# * `data::Union{Vector{ComplexF64},Nothing}` - the vector of coefficients ## Constructor # GroupedCoefficientsComplex( setting, data = nothing ) @@ -343,28 +394,31 @@ def __eq__(self, other): ## Additional Constructor # GroupedCoefficients( setting, data = nothing ) + class GroupedCoefficientsComplex(GC): - def __init__(self, settings, data = None): + def __init__(self, settings, data=None): - try: + try: N = sum(s.mode.datalength(s.bandwidths) for s in settings) except: - raise RuntimeError("The mode is not supportet yet or does not have the function datalength.") + raise RuntimeError( + "The mode is not supportet yet or does not have the function datalength." + ) if data is None: data = np.zeros(N, dtype=complex) if len(data) != N: raise ValueError("The supplied data vector has the wrong length.") - self.settings = settings - self.data = data + self.settings = settings + self.data = data # A class to hold real coefficients belonging to indices in a grouped index set ### Fields -#* `setting::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}` - uniquely describes the setting such as the bandlimits ``N_{\pmb u}``, see also [`get_setting(system::String,d::Int,ds::Int,N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) and [`get_setting(system::String,U::Vector{Vector{Int}},N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) -#* `data::Union{Vector{ComplexF64},Nothing}` - the vector of coefficients +# * `setting::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}` - uniquely describes the setting such as the bandlimits ``N_{\pmb u}``, see also [`get_setting(system::String,d::Int,ds::Int,N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) and [`get_setting(system::String,U::Vector{Vector{Int}},N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) +# * `data::Union{Vector{ComplexF64},Nothing}` - the vector of coefficients ## Constructor # GroupedCoefficientsComplex( setting, data = nothing ) @@ -372,25 +426,28 @@ def __init__(self, settings, data = None): ## Additional Constructor # GroupedCoefficients( setting, data = nothing ) -class GroupedCoefficientsReal(GC): - def __init__(self, settings, data = None): - try: +class GroupedCoefficientsReal(GC): + def __init__(self, settings, data=None): + + try: N = sum(s.mode.datalength(s.bandwidths) for s in settings) except: - raise RuntimeError("The mode is not supportet yet or does not have the function datalength.") + raise RuntimeError( + "The mode is not supportet yet or does not have the function datalength." + ) if data is None: data = np.zeros(N, dtype=float) if len(data) != N: raise ValueError("The supplied data vector has the wrong length.") - self.settings = settings - self.data = data + self.settings = settings + self.data = data + -def GroupedCoefficients(settings, data = None): +def GroupedCoefficients(settings, data=None): if settings[0].mode == NFFTtools: return GroupedCoefficientsComplex(settings, data) if settings[0].mode == NFCTtools or settings[0].mode == CWWTtools: return GroupedCoefficientsReal(settings, data) - diff --git a/src/pyGroupedTransforms/GroupedTransform.py b/src/pyGroupedTransforms/GroupedTransform.py index 8733c25..d730e11 100644 --- a/src/pyGroupedTransforms/GroupedTransform.py +++ b/src/pyGroupedTransforms/GroupedTransform.py @@ -1,84 +1,122 @@ import numpy as np + from pyGroupedTransforms import * # All code that is linked to NFMTtools or to system = "mixed" is not tested yet.... -systems = {"exp" : NFFTtools, "cos" : NFCTtools, "chui1" : CWWTtools, "chui2": CWWTtools, "chui3": CWWTtools, "chui4": CWWTtools, "mixed" : NFMTtools} +systems = { + "exp": NFFTtools, + "cos": NFCTtools, + "chui1": CWWTtools, + "chui2": CWWTtools, + "chui3": CWWTtools, + "chui4": CWWTtools, + "mixed": NFMTtools, +} + def get_superposition_set(d, ds): """ - get_superposition_set( d::Int, ds::Int )::Vector{Vector{Int}} + get_superposition_set( d::Int, ds::Int )::Vector{Vector{Int}} - This function returns ``U^{(d,ds)} = \{ \pmb u \subset \{1,2,\dots,d\} : |\pmb u| \leq ds \}``. + This function returns ``U^{(d,ds)} = \{ \pmb u \subset \{1,2,\dots,d\} : |\pmb u| \leq ds \}``. """ nset = [[j] for j in range(d)] returnset = [[]] + nset - for i in range(ds -1): + for i in range(ds - 1): nextnset = [] for s in nset: for j in range(d): - if s[-1] < j: - nextnset.append(s + [j]) + if s[-1] < j: + nextnset.append(s + [j]) returnset = returnset + nextnset nset = nextnset return [tuple(item) for item in returnset] -def get_setting(system, N, U = None, d = None, ds = None): #I have to add bases and NFMT stuff here! +def get_setting( + system, N, U=None, d=None, ds=None +): # I have to add bases and NFMT stuff here! if system not in systems: - raise ValueError("System not found.") - + raise ValueError("System not found.") - if d != None: # Input system::String, d::Int, ds::Int, N::Vector{Int} + if d != None: # Input system::String, d::Int, ds::Int, N::Vector{Int} - if len(N) != ds: raise ValueError("N must have ds entries.") tmp = np.concatenate(([0], N)) U = get_superposition_set(d, ds) - bwl = [np.full(len(u),tmp[len(u)], 'int32') for u in U] + bwl = [np.full(len(u), tmp[len(u)], "int32") for u in U] if systems[system] == NFMTtools: if len(basis_vect) == 0: - ValueError("please call get_setting with basis_vect for a NFMT transform.") + ValueError( + "please call get_setting with basis_vect for a NFMT transform." + ) if len(basis_vect) < d: ValueError("basis_vect must have an entry for every dimension.") - return [ Setting(u=U[idx], mode=systems[system], bandwidths=np.array(bwl[idx], 'int32'), bases = basis_vect[U[idk]]) for idx in range(len(U)) ] - - else: - return [ Setting(u=U[idx], mode=systems[system], bandwidths=np.array(bwl[idx], 'int32')) for idx in range(len(U)) ] - + return [ + Setting( + u=U[idx], + mode=systems[system], + bandwidths=np.array(bwl[idx], "int32"), + bases=basis_vect[U[idk]], + ) + for idx in range(len(U)) + ] + else: + return [ + Setting( + u=U[idx], + mode=systems[system], + bandwidths=np.array(bwl[idx], "int32"), + ) + for idx in range(len(U)) + ] if len(N) != len(U): - raise ValueError("N must have |U| entries.") + raise ValueError("N must have |U| entries.") - if type(N[0]) == int: # Input system::String, U::list{list{Int}}, N::list{Int} + if type(N[0]) == int: # Input system::String, U::list{list{Int}}, N::list{Int} bwl = [None] * len(U) for i, u in enumerate(U): if len(u) == 0: - bwl[i] = [0] * len(u) + bwl[i] = [0] * len(u) else: bwl[i] = np.full(len(u), N[i]) - if systems[system] == NFMTtools: if len(basis_vect) == 0: - ValueError("please call get_setting with basis_vect for a NFMT transform.") + ValueError( + "please call get_setting with basis_vect for a NFMT transform." + ) if len(basis_vect) < max(max(u) for u in U): ValueError("basis_vect must have an entry for every dimension.") - return [ Setting(u=u, mode=systems[system], bandwidths=np.array(bwl[i], 'int32'), bases = basis_vect[u]) for i,u in enumerate(U) ] - - else: - return [Setting(u=u, mode=systems[system], bandwidths=np.array(bwl[i], 'int32')) for i,u in enumerate(U)] + return [ + Setting( + u=u, + mode=systems[system], + bandwidths=np.array(bwl[i], "int32"), + bases=basis_vect[u], + ) + for i, u in enumerate(U) + ] + else: + return [ + Setting(u=u, mode=systems[system], bandwidths=np.array(bwl[i], "int32")) + for i, u in enumerate(U) + ] - elif type(N[0]) == list or type(N[0]) == np.ndarray: # Input system::String, U::list{list{Int}}, N::list{list{Int}} oder list{np.ndarray} + elif ( + type(N[0]) == list or type(N[0]) == np.ndarray + ): # Input system::String, U::list{list{Int}}, N::list{list{Int}} oder list{np.ndarray} bwl = [None] * len(U) @@ -87,30 +125,49 @@ def get_setting(system, N, U = None, d = None, ds = None): #I have to add b bwl[i] = [0] * len(u) else: if len(N[i]) != len(u): - raise ValueError("Vector N has for the set", u, "not the right length") + raise ValueError( + "Vector N has for the set", u, "not the right length" + ) bwl[i] = N[i] - + if systems[system] == NFMTtools: if len(basis_vect) == 0: - ValueError("please call get_setting with basis_vect for a NFMT transform.") + ValueError( + "please call get_setting with basis_vect for a NFMT transform." + ) if len(basis_vect) < max(max(u) for u in U): ValueError("basis_vect must have an entry for every dimension.") - return [ Setting(u=U[idx], mode=systems[system], bandwidths=np.array(bwl[idx], 'int32'), bases = basis_vect[U[idk]]) for idx in range(len(U)) ] - + return [ + Setting( + u=U[idx], + mode=systems[system], + bandwidths=np.array(bwl[idx], "int32"), + bases=basis_vect[U[idk]], + ) + for idx in range(len(U)) + ] + else: - return [Setting(u=U[idx], mode=systems[system], bandwidths=np.array(bwl[idx], 'int32')) for idx in range(len(U))] + return [ + Setting( + u=U[idx], + mode=systems[system], + bandwidths=np.array(bwl[idx], "int32"), + ) + for idx in range(len(U)) + ] ###GroupedTransform -#A struct to describe a GroupedTransformation +# A struct to describe a GroupedTransformation ## Fields -#* `system::String` - choice of `"exp"` or `"cos"` or `"chui1"` or `"chui2"` or `"chui3"` or `"chui4"` or `"mixed"` -#* `setting::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}` - vector of the dimensions, mode, bandwidths and bases for each term/group, see also [`get_setting(system::String,d::Int,ds::Int,N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) and [`get_setting(system::String,U::Vector{Vector{Int}},N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) -#* `X::Array{Float64}` - array of nodes -#* `transforms::Vector{LinearMap}` - holds the low-dimensional sub transformations -#* `basis_vect::Vector{String}` - holds for every dimension if a cosinus basis [true] or exponential basis [false] is used +# * `system::String` - choice of `"exp"` or `"cos"` or `"chui1"` or `"chui2"` or `"chui3"` or `"chui4"` or `"mixed"` +# * `setting::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}` - vector of the dimensions, mode, bandwidths and bases for each term/group, see also [`get_setting(system::String,d::Int,ds::Int,N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) and [`get_setting(system::String,U::Vector{Vector{Int}},N::Vector{Int},basis_vect::Vector{String})::Vector{NamedTuple{(:u, :mode, :bandwidths, :bases),Tuple{Vector{Int},Module,Vector{Int},Vector{String}}}}`](@ref) +# * `X::Array{Float64}` - array of nodes +# * `transforms::Vector{LinearMap}` - holds the low-dimensional sub transformations +# * `basis_vect::Vector{String}` - holds for every dimension if a cosinus basis [true] or exponential basis [false] is used # ## Constructor # GroupedTransform( system, setting, X, basis_vect::Vector{String} = Vector{String}([]) ) @@ -119,18 +176,32 @@ def get_setting(system, N, U = None, d = None, ds = None): #I have to add b # GroupedTransform( system, d, ds, N::Vector{Int}, X, basis_vect::Vector{String} = Vector{String}([]) ) # GroupedTransform( system, U, N, X, basis_vect::Vector{String} = Vector{String}([]) ) + class GroupedTransform: - def __init__(self, system, X, settings = [], fastmult=True, basis_vect=[], N = [], U = None, d = None, ds = None): + def __init__( + self, + system, + X, + settings=[], + fastmult=True, + basis_vect=[], + N=[], + U=None, + d=None, + ds=None, + ): if system not in systems: raise ValueError("System not found.") if system == "mixed": if len(basis_vect) == 0: - ValueError("please call GroupedTransform with basis_vect for a NFMT transform.") + ValueError( + "please call GroupedTransform with basis_vect for a NFMT transform." + ) if len(basis_vect) != X.shape[1]: ValueError("basis_vect must have an entry for every dimension.") - + if system in {"exp", "chui1", "chui2", "chui3", "chui4"}: if np.min(X) < -0.5 or np.max(X) >= 0.5: raise ValueError("Nodes must be between -0.5 and 0.5.") @@ -143,43 +214,56 @@ def __init__(self, system, X, settings = [], fastmult=True, basis_vect=[], N = [ cosine_mask = basis_vals > 0 if np.sum(cosine_mask) > 0: if (np.min(X[cosine_mask, :]) < 0) or (np.max(X[cosine_mask, :]) > 1): - raise ValueError("Nodes must be between 0 and 1 for cosine or Chebyshev dimensions.") + raise ValueError( + "Nodes must be between 0 and 1 for cosine or Chebyshev dimensions." + ) exp_mask = ~cosine_mask if np.sum(exp_mask) > 0: if (np.min(X[exp_mask, :]) < -0.5) or (np.max(X[exp_mask, :]) > 0.5): - raise ValueError("Nodes must be between -0.5 and 0.5 for exponentional dimensions.") + raise ValueError( + "Nodes must be between -0.5 and 0.5 for exponentional dimensions." + ) if system in {"chui1", "chui2", "chui3", "chui4"}: fastmult = True self.fastmult = fastmult self.basis_vect = basis_vect - self.system = system + self.system = system self.X = X if len(settings) == 0: - self.settings = get_setting(system = system, N = N, U = U, d = d, ds = ds) - else: + self.settings = get_setting(system=system, N=N, U=U, d=d, ds=ds) + else: self.settings = settings - - if fastmult: - self.matrix = np.empty((0,0), dtype=object) - self.transforms = [DeferredLinearOperator() for i in range(len(self.settings))] + if fastmult: + self.matrix = np.empty((0, 0), dtype=object) + self.transforms = [ + DeferredLinearOperator() for i in range(len(self.settings)) + ] for idx, s in enumerate(self.settings): if len(s.bandwidths) == 0: u = (0,) else: u = s.u - if system.startswith("chui"): + if system.startswith("chui"): Order = int(system[-1]) - self.transforms[idx] = s.mode.get_transform(bandwidths = s.bandwidths, X = np.copy(X[:, u], order = "C"), m = Order) + self.transforms[idx] = s.mode.get_transform( + bandwidths=s.bandwidths, X=np.copy(X[:, u], order="C"), m=Order + ) elif system == "mixed": - self.transforms[idx] = s.mode.get_transform(bandwidths = s.bandwidths, X = np.copy(X[:, u], order = "C"), bases = s.bases) + self.transforms[idx] = s.mode.get_transform( + bandwidths=s.bandwidths, + X=np.copy(X[:, u], order="C"), + bases=s.bases, + ) else: - self.transforms[idx] = s.mode.get_transform(bandwidths = s.bandwidths, X = np.copy(X[:, u], order = "C")) - else: + self.transforms[idx] = s.mode.get_transform( + bandwidths=s.bandwidths, X=np.copy(X[:, u], order="C") + ) + else: self.transforms = [] s1 = self.settings[0] if len(s1.bandwidths) == 0: @@ -187,36 +271,53 @@ def __init__(self, system, X, settings = [], fastmult=True, basis_vect=[], N = [ else: u1 = s1.u if system.startswith("chui"): - raise ValueError("Direct computation with full matrix not supported for wavelet basis.") + raise ValueError( + "Direct computation with full matrix not supported for wavelet basis." + ) elif system == "mixed": - matrix = np.array(s1.mode.get_matrix(bandwidths = s1.bandwidths, X = X[:, u1].T), bases = s1.bases) + matrix = np.array( + s1.mode.get_matrix(bandwidths=s1.bandwidths, X=X[:, u1].T), + bases=s1.bases, + ) for s in self.settings[1:]: if len(s.bandwidths) == 0: - u = (0,) + u = (0,) else: u = s.u - matrix = np.hstack([matrix, np.array(s.mode.get_matrix(s.bandwidths, X[:, u].T), bases = s.bases)]) + matrix = np.hstack( + [ + matrix, + np.array( + s.mode.get_matrix(s.bandwidths, X[:, u].T), + bases=s.bases, + ), + ] + ) else: - matrix = np.array(s1.mode.get_matrix(bandwidths = s1.bandwidths, X = X[:, u1].T)) + matrix = np.array( + s1.mode.get_matrix(bandwidths=s1.bandwidths, X=X[:, u1].T) + ) for s in self.settings[1:]: if len(s.bandwidths) == 0: - u = (0,) + u = (0,) else: u = s.u - matrix = np.hstack([matrix, np.array(s.mode.get_matrix(s.bandwidths, X[:, u].T))]) + matrix = np.hstack( + [matrix, np.array(s.mode.get_matrix(s.bandwidths, X[:, u].T))] + ) self.matrix = matrix def __mul__(self, other): """ - (*( F::GroupedTransform, fhat::GroupedCoefficients )::Vector{<:Number}) + (*( F::GroupedTransform, fhat::GroupedCoefficients )::Vector{<:Number}) - If other (= fhat) is an object of GC, this function - overloads the `*` notation in order to achieve `f = F*fhat`. + If other (= fhat) is an object of GC, this function + overloads the `*` notation in order to achieve `f = F*fhat`. - (*( F::GroupedTransform, f::Vector{<:Number} )::GroupedCoefficients) + (*( F::GroupedTransform, f::Vector{<:Number} )::GroupedCoefficients) - If other (= f) is an numpy.ndarray, this function - overloads the * notation in order to achieve the adjoint transform `f = F*f`. + If other (= f) is an numpy.ndarray, this function + overloads the * notation in order to achieve the adjoint transform `f = F*f`. """ if isinstance(other, np.ndarray): # `f = F*f` (f = other) @@ -225,7 +326,7 @@ def __mul__(self, other): fhat = GroupedCoefficients(self.settings) def adjoint_worker(i): - adjoint_result = self.transforms[i].H @ other + adjoint_result = self.transforms[i].H @ other fhat[self.settings[i].u] = adjoint_result for i in range(len(self.transforms)): @@ -237,11 +338,15 @@ def adjoint_worker(i): t.join() return fhat - else: - return GroupedCoefficients(self.settings, (self.matrix.conj()).T @ other) - elif isinstance(other, GC): #`f = F*fhat` (fhat = other) - if self.settings != other.settings: - raise ValueError("The GroupedTransform and the GroupedCoefficients have different settings") + else: + return GroupedCoefficients( + self.settings, (self.matrix.conj()).T @ other + ) + elif isinstance(other, GC): # `f = F*fhat` (fhat = other) + if self.settings != other.settings: + raise ValueError( + "The GroupedTransform and the GroupedCoefficients have different settings" + ) if self.fastmult: threads = [] @@ -258,36 +363,39 @@ def worker(i): threads.append(t) for t in threads: t.join() - + return sum(results) else: - return self.matrix @ other.data + return self.matrix @ other.data def __matmul__(self, other): return self.__mul__(other) def adjoint(self): """ - adjoint( F::GroupedTransform )::GroupedTransform + adjoint( F::GroupedTransform )::GroupedTransform - Overloads the `F'` notation and gives back the same GroupdTransform. GroupedTransform decides by the input if it is the normal trafo or the adjoint so this is only for convinience. + Overloads the `F'` notation and gives back the same GroupdTransform. GroupedTransform decides by the input if it is the normal trafo or the adjoint so this is only for convinience. """ return self @property def H(self): """ - Overloads the `F'` notation and gives back the same GroupdTransform. GroupedTransform decides by the input if it is the normal trafo or the adjoint so this is only for convinience. + Overloads the `F'` notation and gives back the same GroupdTransform. GroupedTransform decides by the input if it is the normal trafo or the adjoint so this is only for convinience. """ return self def __getitem__(self, u): """ - F::GroupedTransform[u::Vector{Int}]::LinearMap{<:Number} or SparseArray + F::GroupedTransform[u::Vector{Int}]::LinearMap{<:Number} or SparseArray - This function overloads getindex of GroupedTransform such that you can do `F[[1,3]]` to obtain the transform of the corresponding ANOVA term defined by `u`. + This function overloads getindex of GroupedTransform such that you can do `F[[1,3]]` to obtain the transform of the corresponding ANOVA term defined by `u`. """ - idx = next((i for i, s in enumerate(self.settings) if len(s.u) == len(u) and s.u == u), None) + idx = next( + (i for i, s in enumerate(self.settings) if len(s.u) == len(u) and s.u == u), + None, + ) if idx is None: raise ValueError("This term is not contained") @@ -303,14 +411,16 @@ def get_matrix(self): This function returns the actual matrix of the transformation. This is not available for the wavelet basis. """ if self.system in ["chui1", "chui2", "chui3", "chui4"]: - raise ValueError("Direct computation with full matrix not supported for wavelet basis.") + raise ValueError( + "Direct computation with full matrix not supported for wavelet basis." + ) elif self.system == "mixed": s1 = self.settings[0] if len(s1.bandwidths) == 0: u1 = (0,) else: u1 = s1.u - F_direct = s1.mode.get_matrix(s.bdanwidths, self.X[:, u1].T, bases = s.bases) + F_direct = s1.mode.get_matrix(s.bdanwidths, self.X[:, u1].T, bases=s.bases) for idx, s in enumerate(self.settings): if idx == 0: continue @@ -339,6 +449,3 @@ def get_matrix(self): F_direct = np.hstack((F_direct, mat)) return F_direct - - - diff --git a/src/pyGroupedTransforms/GroupedTransforms.py b/src/pyGroupedTransforms/GroupedTransforms.py index 1c663cc..9818a45 100644 --- a/src/pyGroupedTransforms/GroupedTransforms.py +++ b/src/pyGroupedTransforms/GroupedTransforms.py @@ -1,45 +1,48 @@ -from pyGroupedTransforms import * import numpy as np -def get_NumFreq(settings): - if settings[0].mode == CWWTtools: +from pyGroupedTransforms import * + + +def get_NumFreq(settings): + if settings[0].mode == CWWTtools: + def datalength(bandwidths): if len(bandwidths) == 0: return 1 elif len(bandwidths) == 1: - return 2**(bandwidths[0] + 1) - 1 + return 2 ** (bandwidths[0] + 1) - 1 elif len(bandwidths) == 2: - return 2**(bandwidths[0] + 1) * bandwidths[0] + 1 + return 2 ** (bandwidths[0] + 1) * bandwidths[0] + 1 elif len(bandwidths) == 3: n = bandwidths[0] - return 2**n * n**2 + 2**n * n + 2**(n + 1) - 1 + return 2**n * n**2 + 2**n * n + 2 ** (n + 1) - 1 else: d = len(bandwidths) n = bandwidths[0] tmp = 0 for i in range(0, n + 1): tmp += 2**i * math.comb(i + d - 1, d - 1) - - return tmp #in julia kriegen wir hier s zurück? - + + return tmp # in julia kriegen wir hier s zurück? + return sum(datalength(s.bandwidths) for s in settings) else: - return sum([np.prod(s.bandwidths -1) for s in settings]) + return sum([np.prod(s.bandwidths - 1) for s in settings]) def get_IndexSet(settings, d): nf = get_NumFreq(settings) index_set = np.zeros((d, nf), dtype=np.int64) - idx = 0 + idx = 0 for s in settings: - if len(s.u) ==0: + if len(s.u) == 0: idx += 1 continue - nf_u = np.prod(s.bandwidths -1) + nf_u = np.prod(s.bandwidths - 1) - if s.mode in [NFFTtools, NFCTtools]: + if s.mode in [NFFTtools, NFCTtools]: index_set_u = s.mode.index_set_without_zeros(s.bandwidths) elif s.mode == NFMTtools: index_set_u = s.mode.nfmt_index_set_without_zeros(s.bandwidths, s.bases) @@ -47,13 +50,14 @@ def get_IndexSet(settings, d): ValueError(f"Unknown mode: {s.mode}") if len(s.u) == 1: - index_set[s.u[0], idx:idx + nf_u] = index_set_u + index_set[s.u[0], idx : idx + nf_u] = index_set_u else: for i, dim in enumerate(s.u): - index_set[dim, idx:idx + nf_u] = index_set_u[i] + index_set[dim, idx : idx + nf_u] = index_set_u[i] idx += nf_u return index_set + # get_setting ist in GroupedTransform.py diff --git a/src/pyGroupedTransforms/NFCTtools.py b/src/pyGroupedTransforms/NFCTtools.py index 418589a..c85628a 100644 --- a/src/pyGroupedTransforms/NFCTtools.py +++ b/src/pyGroupedTransforms/NFCTtools.py @@ -1,7 +1,11 @@ import numpy as np + from pyGroupedTransforms import * -def datalength(bandwidths: np.ndarray) -> int: #datalength(bandwidths::Vector{Int})::Int + +def datalength( + bandwidths: np.ndarray, +) -> int: # datalength(bandwidths::Vector{Int})::Int """ `N = datalength(bandwidths)` @@ -11,12 +15,15 @@ def datalength(bandwidths: np.ndarray) -> int: #datalength(bandwidths::Vector # Output: * `N::Int` ... length of a Fourier-coefficient with the given bandwidths """ - if bandwidths.ndim != 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim != 1 or bandwidths.dtype != "int32": return "Please use an one-dimensional numpy.array with dtype 'int32' as input" else: return np.prod(bandwidths - 1) -def index_set_without_zeros(bandwidths: np.array) -> np.array: #nfct_index_set_without_zeros(bandwidths::Vector{Int})::Array{Int} + +def index_set_without_zeros( + bandwidths: np.array, +) -> np.array: # nfct_index_set_without_zeros(bandwidths::Vector{Int})::Array{Int} """ `freq = nfct_index_set_without_zeros(bandwidths)` @@ -27,26 +34,30 @@ def index_set_without_zeros(bandwidths: np.array) -> np.array: #nfct_index_se * `freq::Array{Int}` ... all frequencies of the full cube without any vector having a zero entry """ - if bandwidths.ndim > 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim > 1 or bandwidths.dtype != "int32": return "Please use an zero or one-dimensional numpy.array with dtype 'int32' as input" d = len(bandwidths) - if d == 0: return np.array([0]) - if d == 1: return np.array([1] + list(range(2, bandwidths[0]))) + if d == 0: + return np.array([0]) + if d == 1: + return np.array([1] + list(range(2, bandwidths[0]))) bandwidths = bandwidths[::-1] tmp = tuple([[1] + list(range(2, bw)) for bw in bandwidths]) tmp = itertools.product(*(tmp[::-1])) - freq = np.empty((d,np.prod(bandwidths - 1)), dtype=int) - for m,x in enumerate(tmp): + freq = np.empty((d, np.prod(bandwidths - 1)), dtype=int) + for m, x in enumerate(tmp): freq[:, m] = x return freq -def nfct_index_set(bandwidths: np.array) -> np.array: #nfct_index_set(bandwidths::Vector{Int})::Array{Int} +def nfct_index_set( + bandwidths: np.array, +) -> np.array: # nfct_index_set(bandwidths::Vector{Int})::Array{Int} """ `freq = nfct_index_set(bandwidths)` @@ -57,24 +68,28 @@ def nfct_index_set(bandwidths: np.array) -> np.array: #nfct_index_set(bandwid * `freq::Array{Int}` ... all frequencies of the full cube """ - if bandwidths.ndim > 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim > 1 or bandwidths.dtype != "int32": return "Please use an zero or one-dimensional numpy.array with dtype 'int32' as input" d = len(bandwidths) - if d == 0: return np.array([0],dtype = 'int') - if d ==1: return np.array([i for i in range(0,bandwidths[0])], 'int') + if d == 0: + return np.array([0], dtype="int") + if d == 1: + return np.array([i for i in range(0, bandwidths[0])], "int") bandwidths = bandwidths[::-1] tmp = tuple([list(range(0, bw)) for bw in bandwidths]) - tmp = itertools.product(*(tmp[::-1])) + tmp = itertools.product(*(tmp[::-1])) - freq = np.empty((d,np.prod(bandwidths)), dtype=int) - for m,x in enumerate(tmp): + freq = np.empty((d, np.prod(bandwidths)), dtype=int) + for m, x in enumerate(tmp): freq[:, m] = x return freq -def nfct_mask(bandwidths: np.array) -> np.array: #nfft_mask(bandwidths::Vector{Int})::BitArray{1} +def nfct_mask( + bandwidths: np.array, +) -> np.array: # nfft_mask(bandwidths::Vector{Int})::BitArray{1} """ `mask = nfct_index_set(bandwidths)` @@ -84,19 +99,20 @@ def nfct_mask(bandwidths: np.array) -> np.array: #nfft_mask(bandwidths::Vector # Output: * `mask::BitArray{1}` ... mask with size of the full cube having zeros whereever a frequency has at least one zero-element and vice-versa """ - if bandwidths.ndim > 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim > 1 or bandwidths.dtype != "int32": return "Please use an zero or one-dimensional numpy.array with dtype 'int32' as input" freq = nfct_index_set(bandwidths) - nfft_mask = np.empty((1,1), dtype = 'bool') + nfft_mask = np.empty((1, 1), dtype="bool") if freq.ndim == 1: - return (freq != 0) + return freq != 0 else: return [0 not in col for col in freq.T] - -def get_transform(bandwidths: np.array, X: np.array): #get_transform(bandwidths::Vector{Int}, X::Array{Float64})::LinearMap +def get_transform( + bandwidths: np.array, X: np.array +): # get_transform(bandwidths::Vector{Int}, X::Array{Float64})::LinearMap """ `F = get_transform(bandwidths, X) @@ -108,46 +124,57 @@ def get_transform(bandwidths: np.array, X: np.array): #get_transform(bandwidth * `F::LinearMap{Float64}` ... Linear map of the Fourier-transform implemented by the NFCT """ - if bandwidths.ndim > 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim > 1 or bandwidths.dtype != "int32": return "Please use an zero or one-dimensional numpy.array with dtype 'int32' as input" - (M,d) = X.shape + (M, d) = X.shape if len(bandwidths) == 0: - return DeferredLinearOperator(dtype=np.float64, shape=(M, 1), mfunc= lambda fhat: np.full(M, fhat[0]), rmfunc=lambda f: np.array([np.sum(f)])) + return DeferredLinearOperator( + dtype=np.float64, + shape=(M, 1), + mfunc=lambda fhat: np.full(M, fhat[0]), + rmfunc=lambda f: np.array([np.sum(f)]), + ) mask = nfct_mask(bandwidths) N = bandwidths - plan = NFCT(N, M) #ursprünglich in Julia: plan = NFCT(N, M, 2*N, 5) + plan = NFCT(N, M) # ursprünglich in Julia: plan = NFCT(N, M, 2*N, 5) plan.x = X - factor_sqrt = sqrt(2)**d + factor_sqrt = sqrt(2) ** d - def trafo(fhat): - plan.fhat = np.zeros(len(mask), 'float') + def trafo(fhat): + plan.fhat = np.zeros(len(mask), "float") plan.fhat[mask] = fhat plan.fhat *= factor_sqrt - plan.nfct_trafo() + plan.nfct_trafo() return plan.f - def adjoint(f): + def adjoint(f): plan.f = f plan.nfct_adjoint() return plan.fhat[mask] * factor_sqrt N = np.prod(bandwidths - 1) - return DeferredLinearOperator(dtype=np.float64, shape=(M, N), mfunc=trafo, rmfunc=adjoint) + return DeferredLinearOperator( + dtype=np.float64, shape=(M, N), mfunc=trafo, rmfunc=adjoint + ) -def get_multiplier(n): #get_multiplier(n::Int)::Float64 + +def get_multiplier(n): # get_multiplier(n::Int)::Float64 if np.isscalar(n): return 1.0 if n == 0 else np.sqrt(2) else: n = np.asarray(n) abs_support = np.count_nonzero(n) - return 1.0 if abs_support == 0 else np.sqrt(2) ** abs_support + return 1.0 if abs_support == 0 else np.sqrt(2) ** abs_support + -def get_matrix(bandwidths, X): #function get_matrix(bandwidths::Vector{Int}, X::Array{Float64})::Array{Float64} +def get_matrix( + bandwidths, X +): # function get_matrix(bandwidths::Vector{Int}, X::Array{Float64})::Array{Float64} """ `F = get_matrix(bandwidths, X) @@ -166,17 +193,22 @@ def get_matrix(bandwidths, X): #function get_matrix(bandwidths::Vector{Int}, X else: d, M = X.shape - if len(bandwidths) == 0: - return np.ones((M,1), dtype = 'float') + return np.ones((M, 1), dtype="float") freq = index_set_without_zeros(np.array(bandwidths, dtype=np.int32)) - - if d == 1: + + if d == 1: freq = freq.flatten() - F_direct = np.array([get_multiplier(n) * np.cos(2 * np.pi * x * n)for x in X for n in freq]).reshape(M, -1) + F_direct = np.array( + [get_multiplier(n) * np.cos(2 * np.pi * x * n) for x in X for n in freq] + ).reshape(M, -1) else: - F_direct = np.array([get_multiplier(n) * np.prod(np.cos(2 * np.pi * x * n)) for x in X.T for n in freq.T]).reshape(M, -1) + F_direct = np.array( + [ + get_multiplier(n) * np.prod(np.cos(2 * np.pi * x * n)) + for x in X.T + for n in freq.T + ] + ).reshape(M, -1) return F_direct - - diff --git a/src/pyGroupedTransforms/NFFTtools.py b/src/pyGroupedTransforms/NFFTtools.py index ca98de8..43f4f3c 100644 --- a/src/pyGroupedTransforms/NFFTtools.py +++ b/src/pyGroupedTransforms/NFFTtools.py @@ -1,9 +1,13 @@ import numpy as np + from pyGroupedTransforms import * -def datalength(bandwidths: np.ndarray) -> int: #datalength(bandwidths::Vector{Int})::Int + +def datalength( + bandwidths: np.ndarray, +) -> int: # datalength(bandwidths::Vector{Int})::Int """ - + `N = datalength(bandwidths)` # Input: @@ -12,14 +16,17 @@ def datalength(bandwidths: np.ndarray) -> int: #datalength(bandwidths::Vector{In # Output: * `N::Int` ... length of a Fourier-coefficient with the given bandwidths """ - if bandwidths.ndim != 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim != 1 or bandwidths.dtype != "int32": return "Please use an one-dimensional numpy.array with dtype 'int32' as input" else: return np.prod(bandwidths - 1) -def index_set_without_zeros(bandwidths: np.array) -> np.array: #nfft_index_set_without_zeros(bandwidths::Vector{Int})::Array{Int} + +def index_set_without_zeros( + bandwidths: np.array, +) -> np.array: # nfft_index_set_without_zeros(bandwidths::Vector{Int})::Array{Int} """ - + `freq = nfft_index_set_without_zeros(bandwidths)` # Input: @@ -28,24 +35,38 @@ def index_set_without_zeros(bandwidths: np.array) -> np.array: #nfft_index_se # Output: * `freq::Array{Int}` ... all frequencies of the full cube without any vector having a zero entry """ - if bandwidths.ndim > 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim > 1 or bandwidths.dtype != "int32": return "Please use an zero or one-dimensional numpy.array with dtype 'int32' as input" d = len(bandwidths) - if d == 0: return np.array([0]) - if d == 1: return np.array([list(range(-int(bw//2),0)) + list(range(1, int(bw//2))) for bw in bandwidths]) - - tmp = tuple([list(range(-int(bw//2),0)) + list(range(1, int(bw//2))) for bw in bandwidths]) + if d == 0: + return np.array([0]) + if d == 1: + return np.array( + [ + list(range(-int(bw // 2), 0)) + list(range(1, int(bw // 2))) + for bw in bandwidths + ] + ) + + tmp = tuple( + [ + list(range(-int(bw // 2), 0)) + list(range(1, int(bw // 2))) + for bw in bandwidths + ] + ) tmp = itertools.product(*(tmp[::-1])) - freq = np.empty((d,np.prod(bandwidths - 1)), dtype=int) - for m,x in enumerate(tmp): + freq = np.empty((d, np.prod(bandwidths - 1)), dtype=int) + for m, x in enumerate(tmp): freq[:, m] = x return freq -def nfft_index_set(bandwidths: np.array) -> np.array: #nfft_index_set(bandwidths::Vector{Int})::Array{Int} +def nfft_index_set( + bandwidths: np.array, +) -> np.array: # nfft_index_set(bandwidths::Vector{Int})::Array{Int} """ `freq = nfft_index_set(bandwidths)` @@ -56,24 +77,29 @@ def nfft_index_set(bandwidths: np.array) -> np.array: #nfft_index_set(bandwid * `freq::Array{Int}` ... all frequencies of the full cube """ - - if bandwidths.ndim > 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim > 1 or bandwidths.dtype != "int32": return "Please use an zero or one-dimensional numpy.array with dtype 'int32' as input" d = len(bandwidths) - if d == 0: return np.array([0],dtype = 'int') - if d ==1: return np.array(list(range(int(-bandwidths[0]/2), int(bandwidths[0]/2))), dtype = 'int') + if d == 0: + return np.array([0], dtype="int") + if d == 1: + return np.array( + list(range(int(-bandwidths[0] / 2), int(bandwidths[0] / 2))), dtype="int" + ) - tmp = [list(range(int(-bw/2),int(bw/2))) for bw in bandwidths] - tmp = itertools.product(*(tmp[::-1])) + tmp = [list(range(int(-bw / 2), int(bw / 2))) for bw in bandwidths] + tmp = itertools.product(*(tmp[::-1])) - freq = np.empty((d,np.prod(bandwidths)), dtype=int) - for m,x in enumerate(tmp): + freq = np.empty((d, np.prod(bandwidths)), dtype=int) + for m, x in enumerate(tmp): freq[:, m] = x return freq -def nfft_mask(bandwidths: np.array) -> np.array: #nfft_mask(bandwidths::Vector{Int})::BitArray{1} +def nfft_mask( + bandwidths: np.array, +) -> np.array: # nfft_mask(bandwidths::Vector{Int})::BitArray{1} """ `mask = nfft_index_set(bandwidths)` @@ -84,19 +110,20 @@ def nfft_mask(bandwidths: np.array) -> np.array: #nfft_mask(bandwidths::Vector * `mask::BitArray{1}` ... mask with size of the full cube having zeros whereever a frequency has at least one zero-element and vice-versa """ - if bandwidths.ndim > 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim > 1 or bandwidths.dtype != "int32": return "Please use an zero or one-dimensional numpy.array with dtype 'int32' as input" freq = nfft_index_set(bandwidths) - nfft_mask = np.empty((1,1), dtype = 'bool') + nfft_mask = np.empty((1, 1), dtype="bool") if freq.ndim == 1: - return (freq != 0) + return freq != 0 else: return [0 not in col for col in freq.T] - -def get_transform(bandwidths: np.array, X: np.array): #get_transform(bandwidths::Vector{Int}, X::Array{Float64})::LinearMap +def get_transform( + bandwidths: np.array, X: np.array +): # get_transform(bandwidths::Vector{Int}, X::Array{Float64})::LinearMap """ `F = get_transform(bandwidths, X) @@ -108,37 +135,45 @@ def get_transform(bandwidths: np.array, X: np.array): #get_transform(bandwidth * `F::LinearMap{ComplexF64}` ... Linear map of the Fourier-transform implemented by the NFFT """ - if bandwidths.ndim > 1 or bandwidths.dtype != 'int32': + if bandwidths.ndim > 1 or bandwidths.dtype != "int32": return "Please use an zero or one-dimensional numpy.array with dtype 'int32' as input" - (M,d) = np.shape(X) + (M, d) = np.shape(X) if len(bandwidths) == 0: - return DeferredLinearOperator(dtype=np.complex128, shape=(M, 1), mfunc= lambda fhat: np.full(M, fhat[0]), rmfunc=lambda f: np.array([np.sum(f)])) + return DeferredLinearOperator( + dtype=np.complex128, + shape=(M, 1), + mfunc=lambda fhat: np.full(M, fhat[0]), + rmfunc=lambda f: np.array([np.sum(f)]), + ) mask = nfft_mask(bandwidths) N = bandwidths - plan = NFFT(N, M) #ursprünglich in Julia: plan = NFFT(N, M, 2*N, 5) + plan = NFFT(N, M) # ursprünglich in Julia: plan = NFFT(N, M, 2*N, 5) plan.x = X - def trafo(fhat): #function trafo(fhat::Vector{ComplexF64})::Vector{ComplexF64} + def trafo(fhat): # function trafo(fhat::Vector{ComplexF64})::Vector{ComplexF64} plan.fhat = np.zeros(len(mask), dtype=np.complex128) plan.fhat[mask] = fhat - plan.nfft_trafo() + plan.nfft_trafo() return plan.f - def adjoint(f): #function adjoint(f::Vector{ComplexF64})::Vector{ComplexF64} + def adjoint(f): # function adjoint(f::Vector{ComplexF64})::Vector{ComplexF64} plan.f = f plan.nfft_adjoint() return plan.fhat[mask] N = np.prod(bandwidths - 1) - return DeferredLinearOperator(dtype=np.complex128, shape=(M, N), mfunc=trafo, rmfunc=adjoint) + return DeferredLinearOperator( + dtype=np.complex128, shape=(M, N), mfunc=trafo, rmfunc=adjoint + ) - -def get_matrix(bandwidths, X): #get_matrix(bandwidths::Vector{Int}, X::Array{Float64})::Array{ComplexF64} +def get_matrix( + bandwidths, X +): # get_matrix(bandwidths::Vector{Int}, X::Array{Float64})::Array{ComplexF64} """ `F = get_matrix(bandwidths, X) @@ -158,7 +193,7 @@ def get_matrix(bandwidths, X): #get_matrix(bandwidths::Vector{Int}, X::Array{F d, M = X.shape if len(bandwidths) == 0: - return np.ones((M,1), dtype = 'complex') + return np.ones((M, 1), dtype="complex") freq = index_set_without_zeros(np.array(bandwidths, dtype=np.int32)) @@ -166,7 +201,8 @@ def get_matrix(bandwidths, X): #get_matrix(bandwidths::Vector{Int}, X::Array{F freq = freq.flatten() F_direct = np.array([[np.exp(-2j * np.pi * x * n) for n in freq] for x in X]) else: - F_direct = np.array([[np.exp(-2j * np.pi * np.dot(x, n)) for n in freq.T] for x in X.T]) - - return F_direct + F_direct = np.array( + [[np.exp(-2j * np.pi * np.dot(x, n)) for n in freq.T] for x in X.T] + ) + return F_direct diff --git a/src/pyGroupedTransforms/NFMTtools.py b/src/pyGroupedTransforms/NFMTtools.py index b1998a1..ea56cfb 100644 --- a/src/pyGroupedTransforms/NFMTtools.py +++ b/src/pyGroupedTransforms/NFMTtools.py @@ -1,2 +1 @@ RuntimeError("NFMTtools is not implemented yet!") - diff --git a/src/pyGroupedTransforms/__init__.py b/src/pyGroupedTransforms/__init__.py index 2a74481..3f41593 100644 --- a/src/pyGroupedTransforms/__init__.py +++ b/src/pyGroupedTransforms/__init__.py @@ -1,27 +1,25 @@ -from math import pi, sqrt import itertools -from pyNFFT3.flags import * -from pyNFFT3.NFFT import * -from pyNFFT3.NFCT import * -import threading - -#Needed for CWWTtools: +# Needed for CWWTtools: import math +import threading from itertools import permutations -from scipy.sparse import coo_matrix +from math import pi, sqrt +from pyNFFT3.flags import * +from pyNFFT3.NFCT import * +from pyNFFT3.NFFT import * # Needed for GroupedCoefficients from scipy.linalg import circulant +from scipy.sparse import coo_matrix - -#import modules +# import modules from .cardinal_bspline import * from .DeferredLinearOperator import * -from .GroupedCoefficients import * #NFFTtools, NFCTtools, NFMTtools and CWWTtools are imported with GroupedCoefficients +from .GroupedCoefficients import * # NFFTtools, NFCTtools, NFMTtools and CWWTtools are imported with GroupedCoefficients from .GroupedTransform import * from .GroupedTransforms import * - -#export + +# export __all__ = [ "NFFTtools", @@ -36,5 +34,5 @@ "variances", "Setting", "GroupedCoefficientsComplex", - "GroupedCoefficientsReal" -] \ No newline at end of file + "GroupedCoefficientsReal", +] diff --git a/src/pyGroupedTransforms/cardinal_bspline.py b/src/pyGroupedTransforms/cardinal_bspline.py index f7a15da..343c9e6 100644 --- a/src/pyGroupedTransforms/cardinal_bspline.py +++ b/src/pyGroupedTransforms/cardinal_bspline.py @@ -1,8 +1,10 @@ import numpy as np -#Cardinal B-Spline : -def cardinal_bspline(v, order): #function cardinal_bspline(v::Array{Float64}, order::Int)::Array{Float64} +# Cardinal B-Spline : +def cardinal_bspline( + v, order +): # function cardinal_bspline(v::Array{Float64}, order::Int)::Array{Float64} """ % evaluates the centered cardinal B-spline of the given order at v INPUT : @@ -11,7 +13,7 @@ def cardinal_bspline(v, order): #function cardinal_bspline(v::Array{Float64}, y ... function values """ - f = np.zeros_like(v, dtype = float) + f = np.zeros_like(v, dtype=float) for i in range(len(v)): x = v[i] @@ -23,145 +25,1139 @@ def cardinal_bspline(v, order): #function cardinal_bspline(v::Array{Float64}, if order == 1: if j == 1: - y = 1.0000000000000000e+00 + y = 1.0000000000000000e00 else: y = 0 elif order == 2: if j == 1: - y = 1.0000000000000000e+00 + x * (1.0000000000000000e+00) + y = 1.0000000000000000e00 + x * (1.0000000000000000e00) elif j == 2: - y = 1.0000000000000000e+00 + x * (-1.0000000000000000e+00) + y = 1.0000000000000000e00 + x * (-1.0000000000000000e00) else: y = 0 elif order == 3: if j == 1: - y = 1.1250000000000000e+00 + x * (1.5000000000000000e+00 + x * (5.0000000000000000e-01)) + y = 1.1250000000000000e00 + x * ( + 1.5000000000000000e00 + x * (5.0000000000000000e-01) + ) elif j == 2: - y = 7.5000000000000000e-01 + x * (0.0000000000000000e+00 + x * (-1.0000000000000000e+00)) + y = 7.5000000000000000e-01 + x * ( + 0.0000000000000000e00 + x * (-1.0000000000000000e00) + ) elif j == 3: - y = 1.1250000000000000e+00 + x * (-1.5000000000000000e+00 + x * (5.0000000000000000e-01)) + y = 1.1250000000000000e00 + x * ( + -1.5000000000000000e00 + x * (5.0000000000000000e-01) + ) else: y = 0 elif order == 4: if j == 1: - y = 1.3333333333333333e+00 + x * ( 2.0000000000000000e+00 + x * (1.0000000000000000e+00 + x * (1.6666666666666666e-01)) ) + y = 1.3333333333333333e00 + x * ( + 2.0000000000000000e00 + + x * (1.0000000000000000e00 + x * (1.6666666666666666e-01)) + ) elif j == 2: - y = 6.6666666666666663e-01 + x * ( 0.0000000000000000e+00 + x * (-1.0000000000000000e+00 + x * (-5.0000000000000000e-01)) ) + y = 6.6666666666666663e-01 + x * ( + 0.0000000000000000e00 + + x * (-1.0000000000000000e00 + x * (-5.0000000000000000e-01)) + ) elif j == 3: - y = 6.6666666666666663e-01 + x * ( 0.0000000000000000e+00 + x * (-1.0000000000000000e+00 + x * (5.0000000000000000e-01)) ) + y = 6.6666666666666663e-01 + x * ( + 0.0000000000000000e00 + + x * (-1.0000000000000000e00 + x * (5.0000000000000000e-01)) + ) elif j == 4: - y = 1.3333333333333333e+00 + x * ( -2.0000000000000000e+00 + x * (1.0000000000000000e+00 + x * (-1.6666666666666666e-01)) ) + y = 1.3333333333333333e00 + x * ( + -2.0000000000000000e00 + + x * (1.0000000000000000e00 + x * (-1.6666666666666666e-01)) + ) else: y = 0 elif order == 5: if j == 1: - y = 1.6276041666666665e+00 + x * ( 2.6041666666666665e+00 + x * ( 1.5625000000000000e+00 + x * (4.1666666666666663e-01 + x * (4.1666666666666664e-02)) ) ) + y = 1.6276041666666665e00 + x * ( + 2.6041666666666665e00 + + x + * ( + 1.5625000000000000e00 + + x * (4.1666666666666663e-01 + x * (4.1666666666666664e-02)) + ) + ) elif j == 2: - y = 5.7291666666666663e-01 + x * ( -2.0833333333333343e-01 + x * ( -1.2500000000000000e+00 + x * (-8.3333333333333337e-01 + x * (-1.6666666666666666e-01)) ) ) + y = 5.7291666666666663e-01 + x * ( + -2.0833333333333343e-01 + + x + * ( + -1.2500000000000000e00 + + x * (-8.3333333333333337e-01 + x * (-1.6666666666666666e-01)) + ) + ) elif j == 3: - y = 5.9895833333333326e-01 + x * ( 0.0000000000000000e+00 + x * ( -6.2500000000000000e-01 + x * (0.0000000000000000e+00 + x * (2.5000000000000000e-01)) ) ) + y = 5.9895833333333326e-01 + x * ( + 0.0000000000000000e00 + + x + * ( + -6.2500000000000000e-01 + + x * (0.0000000000000000e00 + x * (2.5000000000000000e-01)) + ) + ) elif j == 4: - y = 5.7291666666666663e-01 + x * ( 2.0833333333333343e-01 + x * ( -1.2500000000000000e+00 + x * (8.3333333333333337e-01 + x * (-1.6666666666666666e-01)) ) ) + y = 5.7291666666666663e-01 + x * ( + 2.0833333333333343e-01 + + x + * ( + -1.2500000000000000e00 + + x * (8.3333333333333337e-01 + x * (-1.6666666666666666e-01)) + ) + ) elif j == 5: - y = 1.6276041666666665e+00 + x * ( -2.6041666666666665e+00 + x * ( 1.5625000000000000e+00 + x * (-4.1666666666666663e-01 + x * (4.1666666666666664e-02)) ) ) + y = 1.6276041666666665e00 + x * ( + -2.6041666666666665e00 + + x + * ( + 1.5625000000000000e00 + + x * (-4.1666666666666663e-01 + x * (4.1666666666666664e-02)) + ) + ) else: y = 0 elif order == 6: if j == 1: - y = 2.0249999999999999e+00 + x * ( 3.3750000000000000e+00 + x * ( 2.2500000000000000e+00 + x * ( 7.5000000000000000e-01 + x * (1.2500000000000000e-01 + x * (8.3333333333333332e-03)) ) ) ) + y = 2.0249999999999999e00 + x * ( + 3.3750000000000000e00 + + x + * ( + 2.2500000000000000e00 + + x + * ( + 7.5000000000000000e-01 + + x + * (1.2500000000000000e-01 + x * (8.3333333333333332e-03)) + ) + ) + ) elif j == 2: - y = 4.2499999999999999e-01 + x * ( -6.2500000000000033e-01 + x * ( -1.7500000000000000e+00 + x * ( -1.2500000000000000e+00 + x * (-3.7499999999999994e-01 + x * (-4.1666666666666664e-02)) ) ) ) + y = 4.2499999999999999e-01 + x * ( + -6.2500000000000033e-01 + + x + * ( + -1.7500000000000000e00 + + x + * ( + -1.2500000000000000e00 + + x + * (-3.7499999999999994e-01 + x * (-4.1666666666666664e-02)) + ) + ) + ) elif j == 3: - y = 5.5000000000000004e-01 + x * ( -1.1102230246251565e-16 + x * ( -4.9999999999999983e-01 + x * ( 0.0000000000000000e+00 + x * (2.5000000000000006e-01 + x * (8.3333333333333329e-02)) ) ) ) + y = 5.5000000000000004e-01 + x * ( + -1.1102230246251565e-16 + + x + * ( + -4.9999999999999983e-01 + + x + * ( + 0.0000000000000000e00 + + x + * (2.5000000000000006e-01 + x * (8.3333333333333329e-02)) + ) + ) + ) elif j == 4: - y = 5.5000000000000004e-01 + x * ( 1.1102230246251565e-16 + x * ( -5.0000000000000000e-01 + x * ( 0.0000000000000000e+00 + x * (2.5000000000000000e-01 + x * (-8.3333333333333329e-02)) ) ) ) + y = 5.5000000000000004e-01 + x * ( + 1.1102230246251565e-16 + + x + * ( + -5.0000000000000000e-01 + + x + * ( + 0.0000000000000000e00 + + x + * (2.5000000000000000e-01 + x * (-8.3333333333333329e-02)) + ) + ) + ) elif j == 5: - y = 4.2499999999999999e-01 + x * ( 6.2500000000000022e-01 + x * ( -1.7500000000000000e+00 + x * ( 1.2500000000000000e+00 + x * (-3.7500000000000000e-01 + x * (4.1666666666666664e-02)) ) ) ) + y = 4.2499999999999999e-01 + x * ( + 6.2500000000000022e-01 + + x + * ( + -1.7500000000000000e00 + + x + * ( + 1.2500000000000000e00 + + x + * (-3.7500000000000000e-01 + x * (4.1666666666666664e-02)) + ) + ) + ) elif j == 6: - y = 2.0249999999999999e+00 + x * ( -3.3750000000000000e+00 + x * ( 2.2500000000000000e+00 + x * ( -7.5000000000000000e-01 + x * (1.2500000000000000e-01 + x * (-8.3333333333333332e-03)) ) ) ) + y = 2.0249999999999999e00 + x * ( + -3.3750000000000000e00 + + x + * ( + 2.2500000000000000e00 + + x + * ( + -7.5000000000000000e-01 + + x + * (1.2500000000000000e-01 + x * (-8.3333333333333332e-03)) + ) + ) + ) else: y = 0 elif order == 7: if j == 1: - y = 2.5531467013888887e+00 + x * ( 4.3768229166666659e+00 + x * ( 3.1263020833333335e+00 + x * ( 1.1909722222222221e+00 + x * ( 2.5520833333333331e-01 + x * (2.9166666666666664e-02 + x * (1.3888888888888889e-03)) ) ) ) ) + y = 2.5531467013888887e00 + x * ( + 4.3768229166666659e00 + + x + * ( + 3.1263020833333335e00 + + x + * ( + 1.1909722222222221e00 + + x + * ( + 2.5520833333333331e-01 + + x + * ( + 2.9166666666666664e-02 + + x * (1.3888888888888889e-03) + ) + ) + ) + ) + ) elif j == 2: - y = 1.7955729166666637e-01 + x * ( -1.3197916666666669e+00 + x * ( -2.5703125000000000e+00 + x * ( -1.8472222222222221e+00 + x * ( -6.5625000000000000e-01 + x * ( -1.1666666666666665e-01 + x * (-8.3333333333333332e-03) ) ) ) ) ) + y = 1.7955729166666637e-01 + x * ( + -1.3197916666666669e00 + + x + * ( + -2.5703125000000000e00 + + x + * ( + -1.8472222222222221e00 + + x + * ( + -6.5625000000000000e-01 + + x + * ( + -1.1666666666666665e-01 + + x * (-8.3333333333333332e-03) + ) + ) + ) + ) + ) elif j == 3: - y = 5.1178385416666694e-01 + x * ( 9.1145833333331847e-03 + x * ( -3.5546874999999956e-01 + x * ( 1.2152777777777778e-01 + x * ( 3.2812500000000006e-01 + x * (1.4583333333333334e-01 + x * (2.0833333333333332e-02)) ) ) ) ) + y = 5.1178385416666694e-01 + x * ( + 9.1145833333331847e-03 + + x + * ( + -3.5546874999999956e-01 + + x + * ( + 1.2152777777777778e-01 + + x + * ( + 3.2812500000000006e-01 + + x + * ( + 1.4583333333333334e-01 + + x * (2.0833333333333332e-02) + ) + ) + ) + ) + ) elif j == 4: - y = 5.1102430555555556e-01 + x * ( -1.2952601953960160e-16 + x * ( -4.0104166666666652e-01 + x * ( -7.4014868308343765e-17 + x * ( 1.4583333333333334e-01 + x * ( 0.0000000000000000e+00 + x * (-2.7777777777777776e-02) ) ) ) ) ) + y = 5.1102430555555556e-01 + x * ( + -1.2952601953960160e-16 + + x + * ( + -4.0104166666666652e-01 + + x + * ( + -7.4014868308343765e-17 + + x + * ( + 1.4583333333333334e-01 + + x + * ( + 0.0000000000000000e00 + + x * (-2.7777777777777776e-02) + ) + ) + ) + ) + ) elif j == 5: - y = 5.1178385416666694e-01 + x * ( -9.1145833333330373e-03 + x * ( -3.5546875000000000e-01 + x * ( -1.2152777777777778e-01 + x * ( 3.2812500000000000e-01 + x * ( -1.4583333333333334e-01 + x * (2.0833333333333332e-02) ) ) ) ) ) + y = 5.1178385416666694e-01 + x * ( + -9.1145833333330373e-03 + + x + * ( + -3.5546875000000000e-01 + + x + * ( + -1.2152777777777778e-01 + + x + * ( + 3.2812500000000000e-01 + + x + * ( + -1.4583333333333334e-01 + + x * (2.0833333333333332e-02) + ) + ) + ) + ) + ) elif j == 6: - y = 1.7955729166666645e-01 + x * ( 1.3197916666666669e+00 + x * ( -2.5703125000000000e+00 + x * ( 1.8472222222222221e+00 + x * ( -6.5625000000000000e-01 + x * ( 1.1666666666666665e-01 + x * (-8.3333333333333332e-03) ) ) ) ) ) + y = 1.7955729166666645e-01 + x * ( + 1.3197916666666669e00 + + x + * ( + -2.5703125000000000e00 + + x + * ( + 1.8472222222222221e00 + + x + * ( + -6.5625000000000000e-01 + + x + * ( + 1.1666666666666665e-01 + + x * (-8.3333333333333332e-03) + ) + ) + ) + ) + ) elif j == 7: - y = 2.5531467013888887e+00 + x * ( -4.3768229166666659e+00 + x * ( 3.1263020833333335e+00 + x * ( -1.1909722222222221e+00 + x * ( 2.5520833333333331e-01 + x * ( -2.9166666666666664e-02 + x * (1.3888888888888889e-03) ) ) ) ) ) + y = 2.5531467013888887e00 + x * ( + -4.3768229166666659e00 + + x + * ( + 3.1263020833333335e00 + + x + * ( + -1.1909722222222221e00 + + x + * ( + 2.5520833333333331e-01 + + x + * ( + -2.9166666666666664e-02 + + x * (1.3888888888888889e-03) + ) + ) + ) + ) + ) else: y = 0 elif order == 8: if j == 1: - y = 3.2507936507936499e+00 + x * ( 5.6888888888888891e+00 + x * ( 4.2666666666666666e+00 + x * ( 1.7777777777777775e+00 + x * ( 4.4444444444444448e-01 + x * ( 6.6666666666666666e-02 + x * ( 5.5555555555555558e-03 + x * (1.9841269841269841e-04) ) ) ) ) ) ) + y = 3.2507936507936499e00 + x * ( + 5.6888888888888891e00 + + x + * ( + 4.2666666666666666e00 + + x + * ( + 1.7777777777777775e00 + + x + * ( + 4.4444444444444448e-01 + + x + * ( + 6.6666666666666666e-02 + + x + * ( + 5.5555555555555558e-03 + + x * (1.9841269841269841e-04) + ) + ) + ) + ) + ) + ) elif j == 2: - y = -2.2063492063492082e-01 + x * ( -2.4111111111111110e+00 + x * ( -3.8333333333333330e+00 + x * ( -2.7222222222222219e+00 + x * ( -1.0555555555555556e+00 + x * ( -2.3333333333333334e-01 + x * ( -2.7777777777777780e-02 + x * (-1.3888888888888889e-03) ) ) ) ) ) ) + y = -2.2063492063492082e-01 + x * ( + -2.4111111111111110e00 + + x + * ( + -3.8333333333333330e00 + + x + * ( + -2.7222222222222219e00 + + x + * ( + -1.0555555555555556e00 + + x + * ( + -2.3333333333333334e-01 + + x + * ( + -2.7777777777777780e-02 + + x * (-1.3888888888888889e-03) + ) + ) + ) + ) + ) + ) elif j == 3: - y = 4.9047619047619057e-01 + x * ( 7.7777777777777932e-02 + x * ( -9.9999999999999520e-02 + x * ( 3.8888888888888917e-01 + x * ( 5.0000000000000000e-01 + x * ( 2.3333333333333336e-01 + x * ( 4.9999999999999996e-02 + x * (4.1666666666666666e-03) ) ) ) ) ) ) + y = 4.9047619047619057e-01 + x * ( + 7.7777777777777932e-02 + + x + * ( + -9.9999999999999520e-02 + + x + * ( + 3.8888888888888917e-01 + + x + * ( + 5.0000000000000000e-01 + + x + * ( + 2.3333333333333336e-01 + + x + * ( + 4.9999999999999996e-02 + + x * (4.1666666666666666e-03) + ) + ) + ) + ) + ) + ) elif j == 4: - y = 4.7936507936507977e-01 + x * ( -4.1930744590753681e-16 + x * ( -3.3333333333333293e-01 + x * ( -1.3481279584734043e-16 + x * ( 1.1111111111111112e-01 + x * ( 0.0000000000000000e+00 + x * ( -2.7777777777777783e-02 + x * (-6.9444444444444432e-03) ) ) ) ) ) ) + y = 4.7936507936507977e-01 + x * ( + -4.1930744590753681e-16 + + x + * ( + -3.3333333333333293e-01 + + x + * ( + -1.3481279584734043e-16 + + x + * ( + 1.1111111111111112e-01 + + x + * ( + 0.0000000000000000e00 + + x + * ( + -2.7777777777777783e-02 + + x * (-6.9444444444444432e-03) + ) + ) + ) + ) + ) + ) elif j == 5: - y = 4.7936507936507972e-01 + x * ( -9.9127055770103257e-19 + x * ( -3.3333333333333298e-01 + x * ( -8.7231809077690869e-17 + x * ( 1.1111111111111109e-01 + x * ( -1.5860328923216521e-17 + x * ( -2.7777777777777780e-02 + x * (6.9444444444444432e-03) ) ) ) ) ) ) + y = 4.7936507936507972e-01 + x * ( + -9.9127055770103257e-19 + + x + * ( + -3.3333333333333298e-01 + + x + * ( + -8.7231809077690869e-17 + + x + * ( + 1.1111111111111109e-01 + + x + * ( + -1.5860328923216521e-17 + + x + * ( + -2.7777777777777780e-02 + + x * (6.9444444444444432e-03) + ) + ) + ) + ) + ) + ) elif j == 6: - y = 4.9047619047619057e-01 + x * ( -7.7777777777777682e-02 + x * ( -1.0000000000000012e-01 + x * ( -3.8888888888888901e-01 + x * ( 5.0000000000000000e-01 + x * ( -2.3333333333333336e-01 + x * ( 4.9999999999999996e-02 + x * (-4.1666666666666666e-03) ) ) ) ) ) ) + y = 4.9047619047619057e-01 + x * ( + -7.7777777777777682e-02 + + x + * ( + -1.0000000000000012e-01 + + x + * ( + -3.8888888888888901e-01 + + x + * ( + 5.0000000000000000e-01 + + x + * ( + -2.3333333333333336e-01 + + x + * ( + 4.9999999999999996e-02 + + x * (-4.1666666666666666e-03) + ) + ) + ) + ) + ) + ) elif j == 7: - y = -2.2063492063492077e-01 + x * ( 2.4111111111111114e+00 + x * ( -3.8333333333333330e+00 + x * ( 2.7222222222222219e+00 + x * ( -1.0555555555555556e+00 + x * ( 2.3333333333333334e-01 + x * ( -2.7777777777777773e-02 + x * (1.3888888888888889e-03) ) ) ) ) ) ) + y = -2.2063492063492077e-01 + x * ( + 2.4111111111111114e00 + + x + * ( + -3.8333333333333330e00 + + x + * ( + 2.7222222222222219e00 + + x + * ( + -1.0555555555555556e00 + + x + * ( + 2.3333333333333334e-01 + + x + * ( + -2.7777777777777773e-02 + + x * (1.3888888888888889e-03) + ) + ) + ) + ) + ) + ) elif j == 8: - y = 3.2507936507936499e+00 + x * ( -5.6888888888888891e+00 + x * ( 4.2666666666666666e+00 + x * ( -1.7777777777777775e+00 + x * ( 4.4444444444444448e-01 + x * ( -6.6666666666666666e-02 + x * ( 5.5555555555555558e-03 + x * (-1.9841269841269841e-04) ) ) ) ) ) ) + y = 3.2507936507936499e00 + x * ( + -5.6888888888888891e00 + + x + * ( + 4.2666666666666666e00 + + x + * ( + -1.7777777777777775e00 + + x + * ( + 4.4444444444444448e-01 + + x + * ( + -6.6666666666666666e-02 + + x + * ( + 5.5555555555555558e-03 + + x * (-1.9841269841269841e-04) + ) + ) + ) + ) + ) + ) else: y = 0 elif order == 9: if j == 1: - y = 4.1704180036272316e+00 + x * ( 7.4140764508928569e+00 + x * ( 5.7665039062499996e+00 + x * ( 2.5628906249999996e+00 + x * ( 7.1191406250000011e-01 + x * ( 1.2656250000000002e-01 + x * ( 1.4062500000000000e-02 + x * ( 8.9285714285714294e-04 + x * (2.4801587301587302e-05) ) ) ) ) ) ) ) + y = 4.1704180036272316e00 + x * ( + 7.4140764508928569e00 + + x + * ( + 5.7665039062499996e00 + + x + * ( + 2.5628906249999996e00 + + x + * ( + 7.1191406250000011e-01 + + x + * ( + 1.2656250000000002e-01 + + x + * ( + 1.4062500000000000e-02 + + x + * ( + 8.9285714285714294e-04 + + x * (2.4801587301587302e-05) + ) + ) + ) + ) + ) + ) + ) elif j == 2: - y = -8.5608956473214359e-01 + x * ( -4.0750837053571427e+00 + x * ( -5.7226562500000000e+00 + x * ( -4.0023437500000014e+00 + x * ( -1.6328125000000000e+00 + x * ( -4.0937499999999993e-01 + x * ( -6.2499999999999993e-02 + x * ( -5.3571428571428572e-03 + x * (-1.9841269841269841e-04) ) ) ) ) ) ) ) + y = -8.5608956473214359e-01 + x * ( + -4.0750837053571427e00 + + x + * ( + -5.7226562500000000e00 + + x + * ( + -4.0023437500000014e00 + + x + * ( + -1.6328125000000000e00 + + x + * ( + -4.0937499999999993e-01 + + x + * ( + -6.2499999999999993e-02 + + x + * ( + -5.3571428571428572e-03 + + x * (-1.9841269841269841e-04) + ) + ) + ) + ) + ) + ) + ) elif j == 3: - y = 5.0630231584821439e-01 + x * ( 2.8457031250000059e-01 + x * ( 3.8085937500000017e-01 + x * ( 8.8046875000000069e-01 + x * ( 8.0859374999999989e-01 + x * ( 3.7187500000000001e-01 + x * ( 9.3749999999999986e-02 + x * ( 1.2499999999999999e-02 + x * (6.9444444444444447e-04) ) ) ) ) ) ) ) + y = 5.0630231584821439e-01 + x * ( + 2.8457031250000059e-01 + + x + * ( + 3.8085937500000017e-01 + + x + * ( + 8.8046875000000069e-01 + + x + * ( + 8.0859374999999989e-01 + + x + * ( + 3.7187500000000001e-01 + + x + * ( + 9.3749999999999986e-02 + + x + * ( + 1.2499999999999999e-02 + + x * (6.9444444444444447e-04) + ) + ) + ) + ) + ) + ) + ) elif j == 4: - y = 4.5290876116071449e-01 + x * ( -1.9531250000014036e-04 + x * ( -2.8359374999999992e-01 + x * ( -5.4687499999998678e-03 + x * ( 7.0312499999999861e-02 + x * ( -2.1874999999999967e-02 + x * ( -3.7500000000000019e-02 + x * ( -1.2500000000000001e-02 + x * (-1.3888888888888887e-03) ) ) ) ) ) ) ) + y = 4.5290876116071449e-01 + x * ( + -1.9531250000014036e-04 + + x + * ( + -2.8359374999999992e-01 + + x + * ( + -5.4687499999998678e-03 + + x + * ( + 7.0312499999999861e-02 + + x + * ( + -2.1874999999999967e-02 + + x + * ( + -3.7500000000000019e-02 + + x + * ( + -1.2500000000000001e-02 + + x * (-1.3888888888888887e-03) + ) + ) + ) + ) + ) + ) + ) elif j == 5: - y = 4.5292096819196492e-01 + x * ( -4.1199682554449168e-16 + x * ( -2.8222656249999944e-01 + x * ( -2.1510571102112408e-16 + x * ( 8.3984374999999944e-02 + x * ( -1.3877787807814457e-17 + x * ( -1.5625000000000014e-02 + x * ( 8.6736173798840355e-19 + x * (1.7361111111111108e-03) ) ) ) ) ) ) ) + y = 4.5292096819196492e-01 + x * ( + -4.1199682554449168e-16 + + x + * ( + -2.8222656249999944e-01 + + x + * ( + -2.1510571102112408e-16 + + x + * ( + 8.3984374999999944e-02 + + x + * ( + -1.3877787807814457e-17 + + x + * ( + -1.5625000000000014e-02 + + x + * ( + 8.6736173798840355e-19 + + x * (1.7361111111111108e-03) + ) + ) + ) + ) + ) + ) + ) elif j == 6: - y = 4.5290876116071466e-01 + x * ( 1.9531249999965074e-04 + x * ( -2.8359374999999992e-01 + x * ( 5.4687499999997776e-03 + x * ( 7.0312499999999875e-02 + x * ( 2.1874999999999933e-02 + x * ( -3.7500000000000006e-02 + x * ( 1.2499999999999999e-02 + x * (-1.3888888888888887e-03) ) ) ) ) ) ) ) + y = 4.5290876116071466e-01 + x * ( + 1.9531249999965074e-04 + + x + * ( + -2.8359374999999992e-01 + + x + * ( + 5.4687499999997776e-03 + + x + * ( + 7.0312499999999875e-02 + + x + * ( + 2.1874999999999933e-02 + + x + * ( + -3.7500000000000006e-02 + + x + * ( + 1.2499999999999999e-02 + + x * (-1.3888888888888887e-03) + ) + ) + ) + ) + ) + ) + ) elif j == 7: - y = 5.0630231584821439e-01 + x * ( -2.8457031249999992e-01 + x * ( 3.8085937500000017e-01 + x * ( -8.8046875000000024e-01 + x * ( 8.0859374999999989e-01 + x * ( -3.7187500000000001e-01 + x * ( 9.3750000000000000e-02 + x * ( -1.2499999999999997e-02 + x * (6.9444444444444447e-04) ) ) ) ) ) ) ) + y = 5.0630231584821439e-01 + x * ( + -2.8457031249999992e-01 + + x + * ( + 3.8085937500000017e-01 + + x + * ( + -8.8046875000000024e-01 + + x + * ( + 8.0859374999999989e-01 + + x + * ( + -3.7187500000000001e-01 + + x + * ( + 9.3750000000000000e-02 + + x + * ( + -1.2499999999999997e-02 + + x * (6.9444444444444447e-04) + ) + ) + ) + ) + ) + ) + ) elif j == 8: - y = -8.5608956473214382e-01 + x * ( 4.0750837053571427e+00 + x * ( -5.7226562500000000e+00 + x * ( 4.0023437500000005e+00 + x * ( -1.6328125000000000e+00 + x * ( 4.0937499999999993e-01 + x * ( -6.2499999999999986e-02 + x * ( 5.3571428571428563e-03 + x * (-1.9841269841269841e-04) ) ) ) ) ) ) ) + y = -8.5608956473214382e-01 + x * ( + 4.0750837053571427e00 + + x + * ( + -5.7226562500000000e00 + + x + * ( + 4.0023437500000005e00 + + x + * ( + -1.6328125000000000e00 + + x + * ( + 4.0937499999999993e-01 + + x + * ( + -6.2499999999999986e-02 + + x + * ( + 5.3571428571428563e-03 + + x * (-1.9841269841269841e-04) + ) + ) + ) + ) + ) + ) + ) elif j == 9: - y = 4.1704180036272316e+00 + x * ( -7.4140764508928569e+00 + x * ( 5.7665039062499996e+00 + x * ( -2.5628906249999996e+00 + x * ( 7.1191406250000011e-01 + x * ( -1.2656250000000002e-01 + x * ( 1.4062500000000000e-02 + x * ( -8.9285714285714294e-04 + x * (2.4801587301587302e-05) ) ) ) ) ) ) ) + y = 4.1704180036272316e00 + x * ( + -7.4140764508928569e00 + + x + * ( + 5.7665039062499996e00 + + x + * ( + -2.5628906249999996e00 + + x + * ( + 7.1191406250000011e-01 + + x + * ( + -1.2656250000000002e-01 + + x + * ( + 1.4062500000000000e-02 + + x + * ( + -8.9285714285714294e-04 + + x * (2.4801587301587302e-05) + ) + ) + ) + ) + ) + ) + ) else: y = 0 elif order == 10: if j == 1: - y = 5.3822889109347445e+00 + x * ( 9.6881200396825395e+00 + x * ( 7.7504960317460307e+00 + x * ( 3.6168981481481475e+00 + x * ( 1.0850694444444446e+00 + x * ( 2.1701388888888895e-01 + x * ( 2.8935185185185189e-02 + x * ( 2.4801587301587300e-03 + x * ( 1.2400793650793653e-04 + x * (2.7557319223985893e-06) ) ) ) ) ) ) ) ) + y = 5.3822889109347445e00 + x * ( + 9.6881200396825395e00 + + x + * ( + 7.7504960317460307e00 + + x + * ( + 3.6168981481481475e00 + + x + * ( + 1.0850694444444446e00 + + x + * ( + 2.1701388888888895e-01 + + x + * ( + 2.8935185185185189e-02 + + x + * ( + 2.4801587301587300e-03 + + x + * ( + 1.2400793650793653e-04 + + x * (2.7557319223985893e-06) + ) + ) + ) + ) + ) + ) + ) + ) elif j == 2: - y = -1.8416969797178138e+00 + x * ( -6.5658482142857135e+00 + x * ( -8.5034722222222232e+00 + x * ( -5.8645833333333339e+00 + x * ( -2.4704861111111116e+00 + x * ( -6.7187500000000000e-01 + x * ( -1.1921296296296295e-01 + x * ( -1.3392857142857140e-02 + x * ( -8.6805555555555572e-04 + x * (-2.4801587301587302e-05) ) ) ) ) ) ) ) ) + y = -1.8416969797178138e00 + x * ( + -6.5658482142857135e00 + + x + * ( + -8.5034722222222232e00 + + x + * ( + -5.8645833333333339e00 + + x + * ( + -2.4704861111111116e00 + + x + * ( + -6.7187500000000000e-01 + + x + * ( + -1.1921296296296295e-01 + + x + * ( + -1.3392857142857140e-02 + + x + * ( + -8.6805555555555572e-04 + + x * (-2.4801587301587302e-05) + ) + ) + ) + ) + ) + ) + ) + ) elif j == 3: - y = 5.9915123456790131e-01 + x * ( 7.5669642857142883e-01 + x * ( 1.2599206349206369e+00 + x * ( 1.7291666666666661e+00 + x * ( 1.3263888888888888e+00 + x * ( 5.9375000000000000e-01 + x * ( 1.6203703703703703e-01 + x * ( 2.6785714285714281e-02 + x * ( 2.4801587301587300e-03 + x * (9.9206349206349220e-05) ) ) ) ) ) ) ) ) + y = 5.9915123456790131e-01 + x * ( + 7.5669642857142883e-01 + + x + * ( + 1.2599206349206369e00 + + x + * ( + 1.7291666666666661e00 + + x + * ( + 1.3263888888888888e00 + + x + * ( + 5.9375000000000000e-01 + + x + * ( + 1.6203703703703703e-01 + + x + * ( + 2.6785714285714281e-02 + + x + * ( + 2.4801587301587300e-03 + + x * (9.9206349206349220e-05) + ) + ) + ) + ) + ) + ) + ) + ) elif j == 4: - y = 4.2983906525573179e-01 + x * ( -5.2083333333328057e-03 + x * ( -2.6388888888888967e-01 + x * ( -4.8611111111110529e-02 + x * ( -6.9444444444446956e-03 + x * ( -7.2916666666666644e-02 + x * ( -6.0185185185185203e-02 + x * ( -2.0833333333333336e-02 + x * ( -3.4722222222222220e-03 + x * (-2.3148148148148149e-04) ) ) ) ) ) ) ) ) + y = 4.2983906525573179e-01 + x * ( + -5.2083333333328057e-03 + + x + * ( + -2.6388888888888967e-01 + + x + * ( + -4.8611111111110529e-02 + + x + * ( + -6.9444444444446956e-03 + + x + * ( + -7.2916666666666644e-02 + + x + * ( + -6.0185185185185203e-02 + + x + * ( + -2.0833333333333336e-02 + + x + * ( + -3.4722222222222220e-03 + + x * (-2.3148148148148149e-04) + ) + ) + ) + ) + ) + ) + ) + ) elif j == 5: - y = 4.3041776895943618e-01 + x * ( -2.4454782334950824e-18 + x * ( -2.4305555555555569e-01 + x * ( 1.3010426069826053e-16 + x * ( 6.5972222222221863e-02 + x * ( 2.8526563827174158e-17 + x * ( -1.1574074074074106e-02 + x * ( -1.5419764230904951e-18 + x * ( 1.7361111111111110e-03 + x * (3.4722222222222213e-04) ) ) ) ) ) ) ) ) + y = 4.3041776895943618e-01 + x * ( + -2.4454782334950824e-18 + + x + * ( + -2.4305555555555569e-01 + + x + * ( + 1.3010426069826053e-16 + + x + * ( + 6.5972222222221863e-02 + + x + * ( + 2.8526563827174158e-17 + + x + * ( + -1.1574074074074106e-02 + + x + * ( + -1.5419764230904951e-18 + + x + * ( + 1.7361111111111110e-03 + + x * (3.4722222222222213e-04) + ) + ) + ) + ) + ) + ) + ) + ) elif j == 6: - y = 4.3041776895943623e-01 + x * ( -9.9439409253128842e-16 + x * ( -2.4305555555555530e-01 + x * ( -5.7168775886080107e-16 + x * ( 6.5972222222221932e-02 + x * ( -5.7053127654348317e-17 + x * ( -1.1574074074074105e-02 + x * ( 3.0839528461809902e-18 + x * ( 1.7361111111111106e-03 + x * (-3.4722222222222213e-04) ) ) ) ) ) ) ) ) + y = 4.3041776895943623e-01 + x * ( + -9.9439409253128842e-16 + + x + * ( + -2.4305555555555530e-01 + + x + * ( + -5.7168775886080107e-16 + + x + * ( + 6.5972222222221932e-02 + + x + * ( + -5.7053127654348317e-17 + + x + * ( + -1.1574074074074105e-02 + + x + * ( + 3.0839528461809902e-18 + + x + * ( + 1.7361111111111106e-03 + + x * (-3.4722222222222213e-04) + ) + ) + ) + ) + ) + ) + ) + ) elif j == 7: - y = 4.2983906525573234e-01 + x * ( 5.2083333333329627e-03 + x * ( -2.6388888888888878e-01 + x * ( 4.8611111111110640e-02 + x * ( -6.9444444444445291e-03 + x * ( 7.2916666666666644e-02 + x * ( -6.0185185185185182e-02 + x * ( 2.0833333333333336e-02 + x * ( -3.4722222222222216e-03 + x * (2.3148148148148149e-04) ) ) ) ) ) ) ) ) + y = 4.2983906525573234e-01 + x * ( + 5.2083333333329627e-03 + + x + * ( + -2.6388888888888878e-01 + + x + * ( + 4.8611111111110640e-02 + + x + * ( + -6.9444444444445291e-03 + + x + * ( + 7.2916666666666644e-02 + + x + * ( + -6.0185185185185182e-02 + + x + * ( + 2.0833333333333336e-02 + + x + * ( + -3.4722222222222216e-03 + + x * (2.3148148148148149e-04) + ) + ) + ) + ) + ) + ) + ) + ) elif j == 8: - y = 5.9915123456790087e-01 + x * ( -7.5669642857142860e-01 + x * ( 1.2599206349206351e+00 + x * ( -1.7291666666666665e+00 + x * ( 1.3263888888888888e+00 + x * ( -5.9375000000000000e-01 + x * ( 1.6203703703703703e-01 + x * ( -2.6785714285714288e-02 + x * ( 2.4801587301587300e-03 + x * (-9.9206349206349220e-05) ) ) ) ) ) ) ) ) + y = 5.9915123456790087e-01 + x * ( + -7.5669642857142860e-01 + + x + * ( + 1.2599206349206351e00 + + x + * ( + -1.7291666666666665e00 + + x + * ( + 1.3263888888888888e00 + + x + * ( + -5.9375000000000000e-01 + + x + * ( + 1.6203703703703703e-01 + + x + * ( + -2.6785714285714288e-02 + + x + * ( + 2.4801587301587300e-03 + + x * (-9.9206349206349220e-05) + ) + ) + ) + ) + ) + ) + ) + ) elif j == 9: - y = -1.8416969797178138e+00 + x * ( 6.5658482142857135e+00 + x * ( -8.5034722222222232e+00 + x * ( 5.8645833333333330e+00 + x * ( -2.4704861111111112e+00 + x * ( 6.7187499999999989e-01 + x * ( -1.1921296296296295e-01 + x * ( 1.3392857142857140e-02 + x * ( -8.6805555555555551e-04 + x * (2.4801587301587302e-05) ) ) ) ) ) ) ) ) + y = -1.8416969797178138e00 + x * ( + 6.5658482142857135e00 + + x + * ( + -8.5034722222222232e00 + + x + * ( + 5.8645833333333330e00 + + x + * ( + -2.4704861111111112e00 + + x + * ( + 6.7187499999999989e-01 + + x + * ( + -1.1921296296296295e-01 + + x + * ( + 1.3392857142857140e-02 + + x + * ( + -8.6805555555555551e-04 + + x * (2.4801587301587302e-05) + ) + ) + ) + ) + ) + ) + ) + ) elif j == 10: - y = 5.3822889109347445e+00 + x * ( -9.6881200396825395e+00 + x * ( 7.7504960317460307e+00 + x * ( -3.6168981481481475e+00 + x * ( 1.0850694444444446e+00 + x * ( -2.1701388888888895e-01 + x * ( 2.8935185185185189e-02 + x * ( -2.4801587301587300e-03 + x * ( 1.2400793650793653e-04 + x * (-2.7557319223985893e-06) ) ) ) ) ) ) ) ) + y = 5.3822889109347445e00 + x * ( + -9.6881200396825395e00 + + x + * ( + 7.7504960317460307e00 + + x + * ( + -3.6168981481481475e00 + + x + * ( + 1.0850694444444446e00 + + x + * ( + -2.1701388888888895e-01 + + x + * ( + 2.8935185185185189e-02 + + x + * ( + -2.4801587301587300e-03 + + x + * ( + 1.2400793650793653e-04 + + x * (-2.7557319223985893e-06) + ) + ) + ) + ) + ) + ) + ) + ) else: y = 0 f[i] = y return f - diff --git a/tests/cwwt_U.py b/tests/cwwt_U.py index e78ee62..4fe881a 100644 --- a/tests/cwwt_U.py +++ b/tests/cwwt_U.py @@ -1,10 +1,11 @@ -import sys import os +import sys -src_aa = os.path.abspath(os.path.join(os.getcwd(), 'src')) +src_aa = os.path.abspath(os.path.join(os.getcwd(), "src")) sys.path.insert(0, src_aa) import numpy as np + from pyGroupedTransforms.GroupedTransforms import * rng = np.random.default_rng() @@ -12,13 +13,13 @@ d = 4 M = 1000 -X = rng.random((M,d)) - 0.5 +X = rng.random((M, d)) - 0.5 U = [(), (1,), (1, 2)] # set up transform ################################################### -F = GroupedTransform("chui2", X, U = U, N = [3, 2, 1]) +F = GroupedTransform("chui2", X, U=U, N=[3, 2, 1]) # compute transform with NFFT ######################################## @@ -51,6 +52,3 @@ y = rng.random(M) fhat = F.adjoint() * y - - - diff --git a/tests/cwwt_ds.py b/tests/cwwt_ds.py index b904243..00a01d7 100644 --- a/tests/cwwt_ds.py +++ b/tests/cwwt_ds.py @@ -1,10 +1,11 @@ -import sys import os +import sys -src_aa = os.path.abspath(os.path.join(os.getcwd(), 'src')) +src_aa = os.path.abspath(os.path.join(os.getcwd(), "src")) sys.path.insert(0, src_aa) import numpy as np + from pyGroupedTransforms.GroupedTransforms import * rng = np.random.default_rng() @@ -13,11 +14,11 @@ ds = 3 M = 1000 -X = rng.random((M,d)) - 0.5 +X = rng.random((M, d)) - 0.5 # set up transform ################################################### -F = GroupedTransform("chui3", X, d = d, ds = ds, N = [3, 2, 1]) +F = GroupedTransform("chui3", X, d=d, ds=ds, N=[3, 2, 1]) # compute transform with NFFT ######################################## @@ -50,6 +51,3 @@ y = rng.random(M) fhat = F.adjoint() * y - - - diff --git a/tests/nfct_U.py b/tests/nfct_U.py index f67f077..f906da5 100644 --- a/tests/nfct_U.py +++ b/tests/nfct_U.py @@ -1,10 +1,11 @@ -import sys import os +import sys -src_aa = os.path.abspath(os.path.join(os.getcwd(), 'src')) +src_aa = os.path.abspath(os.path.join(os.getcwd(), "src")) sys.path.insert(0, src_aa) import numpy as np + from pyGroupedTransforms.GroupedTransforms import * rng = np.random.default_rng() @@ -12,13 +13,13 @@ d = 4 M = 1000 -X = rng.random((M,d)) * 0.5 +X = rng.random((M, d)) * 0.5 U = [(), (1,), (1, 2)] # set up transform ################################################### -F = GroupedTransform("cos", X, U = U, N = [0, 64, 16]) +F = GroupedTransform("cos", X, U=U, N=[0, 64, 16]) F_direct = F.get_matrix() # compute transform with NFFT ######################################## @@ -72,5 +73,3 @@ error = np.linalg.norm(fhat.vec() - fhat_direct) assert error < 1e-5 - - diff --git a/tests/nfct_ds.py b/tests/nfct_ds.py index ae6b750..63e7244 100644 --- a/tests/nfct_ds.py +++ b/tests/nfct_ds.py @@ -1,10 +1,11 @@ -import sys import os +import sys -src_aa = os.path.abspath(os.path.join(os.getcwd(), 'src')) +src_aa = os.path.abspath(os.path.join(os.getcwd(), "src")) sys.path.insert(0, src_aa) import numpy as np + from pyGroupedTransforms.GroupedTransforms import * rng = np.random.default_rng() @@ -13,11 +14,11 @@ ds = 3 M = 1000 -X = rng.random((M,d)) * 0.5 +X = rng.random((M, d)) * 0.5 # set up transform ################################################### -F = GroupedTransform("cos", X, d = d, ds = ds, N = [2^12, 2^6, 2^4]) +F = GroupedTransform("cos", X, d=d, ds=ds, N=[2 ^ 12, 2 ^ 6, 2 ^ 4]) get_NumFreq(F.settings) get_IndexSet(F.settings, d) F_direct = F.get_matrix() @@ -73,5 +74,3 @@ error = np.linalg.norm(fhat.vec() - fhat_direct) assert error < 1e-5 - - diff --git a/tests/nfft_U.py b/tests/nfft_U.py index 82b364c..aa3ce18 100644 --- a/tests/nfft_U.py +++ b/tests/nfft_U.py @@ -1,10 +1,11 @@ -import sys import os +import sys -src_aa = os.path.abspath(os.path.join(os.getcwd(), 'src')) +src_aa = os.path.abspath(os.path.join(os.getcwd(), "src")) sys.path.insert(0, src_aa) import numpy as np + from pyGroupedTransforms.GroupedTransforms import * rng = np.random.default_rng() @@ -12,13 +13,13 @@ d = 4 M = 1000 -X = rng.random((M,d)) - 0.5 +X = rng.random((M, d)) - 0.5 U = [(), (1,), (1, 2)] # set up transform ################################################### -F = GroupedTransform("exp", X, U = U, N = [0, 64, 16]) +F = GroupedTransform("exp", X, U=U, N=[0, 64, 16]) F_direct = F.get_matrix() # compute transform with NFFT ######################################## @@ -26,17 +27,17 @@ fhat = GroupedCoefficients(F.settings) for i in range(len(F.settings)): u = F.settings[i].u - fhat[u] = rng.random(len(fhat[u])) + 1.j * rng.random(len(fhat[u])) + fhat[u] = rng.random(len(fhat[u])) + 1.0j * rng.random(len(fhat[u])) # arithmetic tests ################################################### ghat = GroupedCoefficients(F.settings) for i in range(len(F.settings)): u = F.settings[i].u - ghat[u] = rng.random(len(ghat[u])) + 1.j * rng.random(len(ghat[u])) + ghat[u] = rng.random(len(ghat[u])) + 1.0j * rng.random(len(ghat[u])) fhat[1] -fhat[1] = 1.0 + 1.j +fhat[1] = 1.0 + 1.0j 2 * fhat fhat + ghat fhat - ghat @@ -58,7 +59,7 @@ # generate random function values #################################### -y = rng.random(M) + 1.j * rng.random(M) +y = rng.random(M) + 1.0j * rng.random(M) # compute adjoint transform with NFFT ################################ @@ -72,5 +73,3 @@ error = np.linalg.norm(fhat.vec() - fhat_direct) assert error < 1e-5 - - diff --git a/tests/nfft_ds.py b/tests/nfft_ds.py index cdcf824..9a60536 100644 --- a/tests/nfft_ds.py +++ b/tests/nfft_ds.py @@ -1,10 +1,11 @@ -import sys import os +import sys -src_aa = os.path.abspath(os.path.join(os.getcwd(), 'src')) +src_aa = os.path.abspath(os.path.join(os.getcwd(), "src")) sys.path.insert(0, src_aa) import numpy as np + from pyGroupedTransforms.GroupedTransforms import * rng = np.random.default_rng() @@ -13,11 +14,11 @@ ds = 3 M = 1000 -X = rng.random((M,d)) - 0.5 +X = rng.random((M, d)) - 0.5 # set up transform ################################################### -F = GroupedTransform("exp", X, d = d, ds = ds, N = [2^12, 2^6, 2^4]) +F = GroupedTransform("exp", X, d=d, ds=ds, N=[2 ^ 12, 2 ^ 6, 2 ^ 4]) get_NumFreq(F.settings) get_IndexSet(F.settings, d) F_direct = F.get_matrix() @@ -27,17 +28,17 @@ fhat = GroupedCoefficients(F.settings) for i in range(len(F.settings)): u = F.settings[i].u - fhat[u] = rng.random(len(fhat[u])) + 1.j * rng.random(len(fhat[u])) + fhat[u] = rng.random(len(fhat[u])) + 1.0j * rng.random(len(fhat[u])) # arithmetic tests ################################################### ghat = GroupedCoefficients(F.settings) for i in range(len(F.settings)): u = F.settings[i].u - ghat[u] = rng.random(len(ghat[u])) + 1.j * rng.random(len(ghat[u])) + ghat[u] = rng.random(len(ghat[u])) + 1.0j * rng.random(len(ghat[u])) fhat[1] -fhat[1] = 1.0 + 1.j +fhat[1] = 1.0 + 1.0j 2 * fhat fhat + ghat fhat - ghat @@ -59,7 +60,7 @@ # generate random function values #################################### -y = rng.random(M) + 1.j * rng.random(M) +y = rng.random(M) + 1.0j * rng.random(M) # compute adjoint transform with NFFT ################################ @@ -73,5 +74,3 @@ error = np.linalg.norm(fhat.vec() - fhat_direct) assert error < 1e-5 - - diff --git a/tests/run_tests.py b/tests/run_tests.py index b586edc..9f0eee3 100644 --- a/tests/run_tests.py +++ b/tests/run_tests.py @@ -1,6 +1,6 @@ +import os import subprocess import sys -import os # Script to run all tests and return any failures