Improved numerical methods for the vacuum module

benchmarks/benchmark_fourier_transforms.jl

@@ -28,21 +28,21 @@ function extract_modes(fft_result, mlow, mhigh, mtheta)
     for (i, m) in enumerate(mlow:mhigh)
         if m >= 0
             # Positive frequencies
-            modes[i] = fft_result[m + 1] / mtheta  # FFT normalization
+            modes[i] = fft_result[m+1] / mtheta  # FFT normalization
         else
             # Negative frequencies (wrap around)
-            modes[i] = fft_result[mtheta + m + 1] / mtheta
+            modes[i] = fft_result[mtheta+m+1] / mtheta
         end
     end
     return modes
 end
 
 # Test configurations
 test_cases = [
-    (name="Small (mtheta=128, mpert=10)",   mtheta=128,  mpert=10,  mlow=-5),
-    (name="Medium (mtheta=256, mpert=20)",  mtheta=256,  mpert=20,  mlow=-10),
-    (name="Large (mtheta=480, mpert=40)",   mtheta=480,  mpert=40,  mlow=-20),
-    (name="Very Large (mtheta=1024, mpert=80)", mtheta=1024, mpert=80, mlow=-40),
+    (name="Small (mtheta=128, mpert=10)", mtheta=128, mpert=10, mlow=-5),
+    (name="Medium (mtheta=256, mpert=20)", mtheta=256, mpert=20, mlow=-10),
+    (name="Large (mtheta=480, mpert=40)", mtheta=480, mpert=40, mlow=-20),
+    (name="Very Large (mtheta=1024, mpert=80)", mtheta=1024, mpert=80, mlow=-40)
 ]
 
 for test in test_cases
@@ -56,7 +56,7 @@ for test in test_cases
     mhigh = mlow + mpert - 1
 
     # Create test data
-    theta = range(0, 2π, length=mtheta+1)[1:end-1]
+    theta = range(0, 2π; length=mtheta+1)[1:(end-1)]
     data = sin.(3 .* theta) .+ 0.5 .* cos.(7 .* theta) .+ 0.2 .* sin.(11 .* theta)
 
     # Initialize FourierTransform
@@ -67,7 +67,8 @@ for test in test_cases
     theta_buffer = zeros(ComplexF64, mtheta)
 
     # Pre-allocate for low-level API
-    cslth, snlth = compute_fourier_coefficients(mtheta, mpert, mlow)
+    exp_mn_basis = compute_fourier_coefficients(mtheta, mpert, mlow)
+    cslth, snlth = real(exp_mn_basis), imag(exp_mn_basis)
     gij = reshape(data, mtheta, 1)  # Matrix form
     gil = zeros(Float64, mtheta, mpert)
 
@@ -99,7 +100,7 @@ for test in test_cases
     # Note: Our transform uses a different normalization and basis
     println("\n--- Accuracy Check ---")
     println("FourierTransform allocating vs in-place: ",
-            @sprintf("%.2e", maximum(abs.(modes_alloc .- modes_buffer))))
+        @sprintf("%.2e", maximum(abs.(modes_alloc .- modes_buffer))))
 
     # Compare magnitudes of modes (since basis might differ)
     println("Mode magnitudes comparison (FourierTransform vs FFTW):")
@@ -129,9 +130,9 @@ for test in test_cases
     full_modes = zeros(ComplexF64, mtheta)
     for (i, m) in enumerate(mlow:mhigh)
         if m >= 0
-            full_modes[m + 1] = modes_test[i]
+            full_modes[m+1] = modes_test[i]
         else
-            full_modes[mtheta + m + 1] = modes_test[i]
+            full_modes[mtheta+m+1] = modes_test[i]
         end
     end
     t6 = @benchmark ifft($full_modes)
@@ -140,21 +141,21 @@ for test in test_cases
     # Accuracy check
     println("\n--- Inverse Accuracy Check ---")
     println("inverse() allocating vs in-place: ",
-            @sprintf("%.2e", maximum(abs.(theta_alloc .- theta_buffer))))
+        @sprintf("%.2e", maximum(abs.(theta_alloc .- theta_buffer))))
     println("Round-trip error (real part):     ",
-            @sprintf("%.2e", maximum(abs.(real.(theta_alloc) .- data))))
+        @sprintf("%.2e", maximum(abs.(real.(theta_alloc) .- data))))
 
     # Performance summary
     println("\n--- Performance Summary ---")
     println(@sprintf("Forward transform speedup (in-place vs allocating): %.2fx",
-            median(t1).time / median(t2).time))
+        median(t1).time / median(t2).time))
     println(@sprintf("Allocations eliminated: %d → %d",
-            t1.allocs, t2.allocs))
+        t1.allocs, t2.allocs))
 
     # Compare to FFTW
     println(@sprintf("\nFourier vs FFTW (forward): %.2fx %s",
-            abs(median(t2).time / median(t3).time),
-            median(t2).time < median(t3).time ? "faster" : "slower"))
+        abs(median(t2).time / median(t3).time),
+        median(t2).time < median(t3).time ? "faster" : "slower"))
     println("Note: FFTW computes full DFT (all N modes), we compute truncated series ($mpert modes)")
 end
 
@@ -169,7 +170,7 @@ mlow = -10
 nbatch = 10  # Transform 10 functions simultaneously
 
 ft = FourierTransform(mtheta, mpert, mlow)
-theta = range(0, 2π, length=mtheta+1)[1:end-1]
+theta = range(0, 2π; length=mtheta+1)[1:(end-1)]
 
 # Create batch data
 data_matrix = zeros(Float64, mtheta, nbatch)
@@ -182,7 +183,7 @@ modes_matrix = zeros(ComplexF64, mpert, nbatch)
 println("\nTransforming $nbatch functions of length $mtheta:")
 
 print("Allocating (loop):  ")
-@btime for i in 1:$nbatch
+@btime for i in 1:($nbatch)
     modes = $ft($data_matrix[:, i])
 end
 
@@ -191,7 +192,7 @@ print("Allocating (matrix):")
 
 print("In-place (loop):    ")
 modes_buffer = zeros(ComplexF64, mpert)
-@btime for i in 1:$nbatch
+@btime for i in 1:($nbatch)
     transform!($modes_buffer, $ft, $data_matrix[:, i])
 end
 
@@ -208,7 +209,8 @@ mpert = 10
 mlow = 1
 
 # Setup for low-level API
-cslth, snlth = compute_fourier_coefficients(mtheta, mpert, mlow)
+exp_mn_basis = compute_fourier_coefficients(mtheta, mpert, mlow)
+cslth, snlth = real(exp_mn_basis), imag(exp_mn_basis)
 gij = randn(mtheta, mtheta)  # Green's function matrix
 gil = zeros(Float64, mtheta, mpert)