added bandwidth detection

Develop

Author: passscoed

.github/workflows/ci.yml

@@ -27,6 +27,7 @@ jobs:
           pip install pyGroupedTransforms
           pip install numpy
           pip install scipy
+          pip install matplotlib
       - name: Run tests
         run: |
           python -m tests.run_tests

.gitignore

@@ -6,6 +6,7 @@ dist/
 *.py[cod]
 __pycache__/
 venv/
+venv2/
 
 bin/
 /lib64/
@@ -16,3 +17,4 @@ src/pyNFFT3/lib/AVX2/libgomp-1.dll
 src/pyNFFT3/lib/AVX2/libgcc_s_seh-1.dll
 src/pyNFFT3/lib/AVX2/libwinpthread-1.dll
 simpleTest/venv/
+bandwidth_detection.py

README.md

@@ -53,3 +53,4 @@ Requirements
 - pyGroupedTransforms 0.1.0 or greater
 - NumPy 2.0.0 or greater
 - SciPy 1.16.0 or greater
+- Matplotlib 3.5 or greater

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "pyANOVAapprox"
-version = "1.0.0"
+version = "2.0.0"
 authors = [
   { name="Felix Wirth", email="fwi012001@gmail.com" },
 ]
@@ -15,9 +15,10 @@ description = "Approximation Package for High-Dimensional Functions"
 readme = "README.md"
 requires-python = ">=3.9"
 dependencies = [
-  "pyGroupedTransforms>=0.1.0",
+  "pyGroupedTransforms>=1.1.0",
   "numpy>=2.0.0",
   "scipy>=1.16.0",
+  "matplotlib>=3.5"
 ]
 
 [project.urls]

simpleTest/TestFunctionPeriodic.py

@@ -0,0 +1,155 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[ ]:
+
+AS = [(), (0,), (1,), (2,), (3,), (4,), (5,), (0, 2), (1, 4), (3, 5)]
+
+import numpy as np
+
+# Coefficients C[0] = C₁, C[1] = C₂, C[2] = C₃
+# (Julia is 1-based; Python is 0-based)
+C = np.array([np.sqrt(0.75), np.sqrt(315 / 604), np.sqrt(277200 / 655177)])
+
+# --- B-spline definitions ---
+
+
+def b_spline_2(x):
+    C_2 = C[0]
+    if 0.0 <= x < 0.5:
+        return C_2 * 4 * x
+    elif 0.5 <= x < 1.0:
+        return C_2 * 4 * (1 - x)
+    else:
+        raise ValueError("B-spline 2: x out of range [0,1)")
+
+
+def b_spline_4(x):
+    C_4 = C[1]
+    if 0.0 <= x < 0.25:
+        return C_4 * (128 / 3) * x**3
+    elif 0.25 <= x < 0.5:
+        return C_4 * (8 / 3 - 32 * x + 128 * x**2 - 128 * x**3)
+    elif 0.5 <= x < 0.75:
+        return C_4 * (-88 / 3 - 256 * x**2 + 160 * x + 128 * x**3)
+    elif 0.75 <= x < 1.0:
+        return C_4 * (128 / 3 - 128 * x + 128 * x**2 - (128 / 3) * x**3)
+    else:
+        raise ValueError("B-spline 4: x out of range [0,1)")
+
+
+def b_spline_6(x):
+    C_6 = C[2]
+    if 0.0 <= x < 1.0 / 6:
+        return C_6 * (1944 / 5) * x**5
+    elif 1.0 / 6 <= x < 2.0 / 6:
+        return C_6 * (
+            3 / 10 - 9 * x + 108 * x**2 - 648 * x**3 + 1944 * x**4 - 1944 * x**5
+        )
+    elif 2.0 / 6 <= x < 0.5:
+        return C_6 * (
+            -237 / 10 + 351 * x - 2052 * x**2 + 5832 * x**3 - 7776 * x**4 + 3888 * x**5
+        )
+    elif 0.5 <= x < 4.0 / 6:
+        return C_6 * (
+            2193 / 10
+            + 7668 * x**2
+            - 2079 * x
+            + 11664 * x**4
+            - 13608 * x**3
+            - 3888 * x**5
+        )
+    elif 4.0 / 6 <= x < 5.0 / 6:
+        return C_6 * (
+            -5487 / 10
+            + 3681 * x
+            - 9612 * x**2
+            + 12312 * x**3
+            - 7776 * x**4
+            + 1944 * x**5
+        )
+    elif 5.0 / 6 <= x < 1.0:
+        return C_6 * (
+            1944 / 5
+            - 1944 * x
+            + 3888 * x**2
+            - 3888 * x**3
+            + 1944 * x**4
+            - (1944 / 5) * x**5
+        )
+    else:
+        raise ValueError("B-spline 6: x out of range [0,1)")
+
+
+# --- Block structure ---
+
+m1 = [0, 2]  # 1-based [1, 3]
+m2 = [1, 4]  # 2-based [2, 5]
+m3 = [3, 5]  # 3-based [4, 6]
+
+# --- Transformed function f(x) ---
+
+
+def trans(x):
+    return x + 1 if x < 0 else x
+
+
+def f(x):
+    x = np.asarray(x)
+    if x.shape != (6,):
+        raise ValueError("f(x): Input must be 6-dimensional.")
+    if np.any(x < -0.5) or np.any(x > 0.5):
+        raise ValueError("f(x): All entries must be in [-0.5, 0.5].")
+
+    xT = np.where(x < 0, x + 1, x)  # vectorized 'trans'
+    return (
+        np.prod([b_spline_2(xT[i]) for i in m1])
+        + np.prod([b_spline_4(xT[i]) for i in m2])
+        + np.prod([b_spline_6(xT[i]) for i in m3])
+    )
+
+
+# --- sinc and b(k, r) function ---
+
+
+def sinc(x):
+    return 1.0 if x == 0.0 else np.sin(x) / x
+
+
+def b(k, r):
+    idx = r // 2 - 1  # Adjust for Python 0-based indexing
+    return C[idx] * (sinc(np.pi * k / r) ** r) * np.cos(np.pi * k)
+
+
+# --- Fourier coefficients fc(k) ---
+
+
+def fc(k):
+    if len(k) != 6:
+        raise ValueError("fc(k): k must be 6-dimensional.")
+
+    ind = np.array([int(ki != 0) for ki in k])
+
+    b2_block = np.sum(ind) == np.sum(ind[m1])
+    b4_block = np.sum(ind) == np.sum(ind[m2])
+    b6_block = np.sum(ind) == np.sum(ind[m3])
+
+    val = 0.0
+    if b2_block:
+        val += np.prod([b(k[i], 2) for i in m1])
+    if b4_block:
+        val += np.prod([b(k[i], 4) for i in m2])
+    if b6_block:
+        val += np.prod([b(k[i], 6) for i in m3])
+    return val
+
+
+# --- Norm computation ---
+
+
+def norm():
+    result = 3.0
+    for i in range(1, 3):  # i = 1, 2
+        for j in range(i + 1, 4):  # j = 2, 3
+            result += 2 * b(0, 2 * i) ** 2 * b(0, 2 * j) ** 2
+    return np.sqrt(result)

simpleTest/exampleAutoapproximate.py

@@ -0,0 +1,66 @@
+# pip install pyANOVAapprox
+
+# Example for approximating an periodic function
+
+import math
+
+import matplotlib.pyplot as plt
+import numpy as np
+from TestFunctionPeriodic import *
+
+import pyANOVAapprox as ANOVAapprox
+
+
+def TestFunction(x):
+    return b_spline_2(x[0]) * b_spline_4(x[1]) * b_spline_6(x[2])
+
+
+rng = np.random.default_rng(1234)
+
+##################################
+## Definition of the parameters ##
+##################################
+
+d = 3  # dimension
+
+M = 100000  # number of used evaluation points to train the model
+M_test = 100000  # number of used evaluation points to test the accuracity the model
+
+U = [(), (0,), (1,), (2,), (0, 1), (0, 2), (1, 2), (0, 1, 2)]
+
+lambdas = np.array([0.0])  # used regularisation parameters λ
+
+############################
+## Generation of the data ##
+############################
+
+X = rng.random((M, d))  # construct the evaluation points for training
+y = np.array(
+    [TestFunction(X[i, :].T) for i in range(M)], dtype=complex
+)  # evaluate the function at these points
+X = X - 0.5
+X_test = rng.random((M_test, d))
+y_test = np.array(
+    [TestFunction(X_test[i, :].T) for i in range(M_test)], dtype=complex
+)  # the same for the test points
+X_test = X_test - 0.5
+
+##########################
+## Do the approximation ##
+##########################
+
+ads = ANOVAapprox.approx(X, y, U=U, basis="per")
+ads.autoapproximate()
+
+################################
+## get approximation accuracy ##
+################################
+
+# mse = ANOVAapprox.get_mse(ads) # get mse error at the given training points
+mse = ads.get_mse(X=X_test, y=y_test)  # get mse error at the test points
+λ_min = min(
+    mse, key=mse.get
+)  # get the regularisation parameter which leads to the minimal error
+mse_min = mse[λ_min]
+
+print("mse = " + str(mse_min))

simpleTest/exampleCheb.py

@@ -103,7 +103,7 @@ def TestFunction(x):
 ar = ANOVAapprox.get_AttributeRanking(ads, λ_min)  # get the attrbute ranking
 
 plt.figure()
-(markers, stemlines, baseline) = plt.stem(
+markers, stemlines, baseline = plt.stem(
     np.arange(1, d + 1),  # x-Werte (1:d)
     ar,  # y-Werte (ar)
     linefmt="C0-",  # Stil der Stiele
@@ -126,7 +126,7 @@ def TestFunction(x):
 l = len(label)
 plt.figure()
 x_values = np.arange(1, l + 1)
-(markers, stemlines, baseline) = plt.stem(
+markers, stemlines, baseline = plt.stem(
     x_values,  # X-Werte: 1 bis l
     gsis,  # Y-Werte: gsis
     linefmt="C0-",  # Stil der Stiele

simpleTest/exampleClassification.py

@@ -100,7 +100,7 @@ def TestFunction(x):
 ar = ANOVAapprox.get_AttributeRanking(ads, 0.0)  # get the attrbute ranking
 
 plt.figure()
-(markers, stemlines, baseline) = plt.stem(
+markers, stemlines, baseline = plt.stem(
     np.arange(1, d + 1),  # x-Werte (1:d)
     ar,  # y-Werte (ar)
     linefmt="C0-",  # Stil der Stiele
@@ -123,7 +123,7 @@ def TestFunction(x):
 l = len(label)
 plt.figure()
 x_values = np.arange(1, l + 1)
-(markers, stemlines, baseline) = plt.stem(
+markers, stemlines, baseline = plt.stem(
     x_values,  # X-Werte: 1 bis l
     gsis,  # Y-Werte: gsis
     linefmt="C0-",  # Stil der Stiele

simpleTest/exampleNonPeriodic.py

@@ -101,7 +101,7 @@ def TestFunction(x):
 ar = ANOVAapprox.get_AttributeRanking(ads, λ_min)  # get the attrbute ranking
 
 plt.figure()
-(markers, stemlines, baseline) = plt.stem(
+markers, stemlines, baseline = plt.stem(
     np.arange(1, d + 1),  # x-Werte (1:d)
     ar,  # y-Werte (ar)
     linefmt="C0-",  # Stil der Stiele
@@ -124,7 +124,7 @@ def TestFunction(x):
 l = len(label)
 plt.figure()
 x_values = np.arange(1, l + 1)
-(markers, stemlines, baseline) = plt.stem(
+markers, stemlines, baseline = plt.stem(
     x_values,  # X-Werte: 1 bis l
     gsis,  # Y-Werte: gsis
     linefmt="C0-",  # Stil der Stiele

simpleTest/examplePeriodic.py

@@ -90,7 +90,7 @@ def TestFunction(x):
 ################################
 
 # mse = ANOVAapprox.get_mse(ads) # get mse error at the given training points
-mse = ANOVAapprox.get_mse(ads, X_test, y_test)  # get mse error at the test points
+mse = ads.get_mse(X=X_test, y=y_test)  # get mse error at the test points
 λ_min = min(
     mse, key=mse.get
 )  # get the regularisation parameter which leads to the minimal error
@@ -103,10 +103,10 @@ def TestFunction(x):
 ###############################################
 
 
-ar = ANOVAapprox.get_AttributeRanking(ads, λ_min)  # get the attrbute ranking
+ar = ads.get_AttributeRanking(lam=λ_min)  # get the attrbute ranking
 
 plt.figure()
-(markers, stemlines, baseline) = plt.stem(
+markers, stemlines, baseline = plt.stem(
     np.arange(1, d + 1),  # x-Werte (1:d)
     ar,  # y-Werte (ar)
     linefmt="C0-",  # Stil der Stiele
@@ -124,12 +124,12 @@ def TestFunction(x):
 plt.show()  # plot the arrtibute ranking in an logplot
 print("active dimensions: " + str(ar[ar > 1e-2]))
 
-gsis = ANOVAapprox.get_GSI(ads, λ_min)
+gsis = ads.get_GSI(lam=λ_min)
 label = list(ads.U[1:])
 l = len(label)
 plt.figure()
 x_values = np.arange(1, l + 1)
-(markers, stemlines, baseline) = plt.stem(
+markers, stemlines, baseline = plt.stem(
     x_values,  # X-Werte: 1 bis l
     gsis,  # Y-Werte: gsis
     linefmt="C0-",  # Stil der Stiele
@@ -169,7 +169,7 @@ def TestFunction(x):
     lam=lambdas, max_iter=max_iter, solver="lsqr"
 )  # do the approximation for all specified regularisation parameters
 
-mse = ANOVAapprox.get_mse(a, X_test, y_test)  # get mse error at the test points
+mse = a.get_mse(X=X_test, y=y_test)  # get mse error at the test points
 λ_min = min(
     mse, key=mse.get
 )  # get the regularisation parameter which leads to the minimal error