feat: Bragg disk detection .single runs on every detector size

Sources/CuboBragg/BraggDisks.swift

@@ -38,11 +38,11 @@ public struct BraggDetector: Sendable {
         singleBox = SingleBox(rows: detectorRows, columns: detectorColumns, kernel: kernel.values)
     }
 
-    /// Whether `BraggPrecision.single` can run on a detector of this size: vDSP's complex DFT must support both lengths
-    /// (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128, 192, 240, 256, 320, 384, 480, 512, ...). Other sizes are refused by the
-    /// `.single` entry points with `CuboBraggError.invalidOption`; `.double` runs on any size.
+    /// Whether `BraggPrecision.single` can run on a detector of this size: every positive size. Lengths vDSP's complex DFT
+    /// takes (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128, 192, 240, 256, 320, 384, 480, 512, ...) run on it directly; any
+    /// other runs through Bluestein's chirp-z on a longer vDSP DFT, at the same float32 tolerances and a higher cost.
     public static func supportsSinglePrecision(detectorRows: Int, detectorColumns: Int) -> Bool {
-        FloatFFT2D.supports(rows: detectorRows, columns: detectorColumns)
+        detectorRows > 0 && detectorColumns > 0
     }
 
     /// The correlation of one pattern: the real map `CC` and the Fourier-domain product `cc` (see the file header).
@@ -116,7 +116,7 @@ public struct BraggDetector: Sendable {
         return result
     }
 
-    /// The Float32 path's state for the whole-cube pass (refuses an unsupported size here, before any tile is read).
+    /// The Float32 path's state for the whole-cube pass (built here, before any tile is read).
     func singleState() throws -> SingleState { try singleBox.state() }
 
     func correlate(values input: [Double], options o: BraggDiskOptions) throws -> (Correlation, Int) {
@@ -266,7 +266,7 @@ public enum BraggDisks {
             throw CuboBraggError.shapeMismatch("the cube's detector is \(shape.detectorRows) × \(shape.detectorColumns), the kernel \(kernel.detectorRows) × \(kernel.detectorColumns)")
         }
         let detector = try BraggDetector(kernel: kernel)
-        let single = options.precision == .single ? try detector.singleState() : nil   // refuses an unsupported size before any read
+        let single = options.precision == .single ? try detector.singleState() : nil   // the Float plan, built once before any read
         let n = shape.detectorCount
         let total = shape.scanCount
         var all: [[BraggPeak]] = []

Sources/CuboBragg/BraggTypes.swift

@@ -66,10 +66,10 @@ public enum BraggPrecision: Sendable, Hashable, CaseIterable {
     case double
     /// Float32 (Accelerate: vDSP convolution and complex DFT) for the pattern copy, the non-finite fill, both Gaussian
     /// filters, the 2-D transform, the correlation product and the maxima search; the few accepted peaks are refined in
-    /// Double. Available only for detector sizes vDSP's DFT supports on both axes (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128,
-    /// 192, 240, 256, 320, 384, 480, 512, ...; `BraggDetector.supportsSinglePrecision(detectorRows:detectorColumns:)` asks);
-    /// any other size is refused with `CuboBraggError.invalidOption`, never computed another way. Agrees with `.double`
-    /// to about 1e-5 of the correlation maximum, hence peak positions to about 1e-4 px (docs/bragg.md has the measured gaps).
+    /// Double. Runs on every detector size: lengths vDSP's DFT takes (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128, 192, 240,
+    /// 256, 320, 384, 480, 512, ...) directly, any other through Bluestein's chirp-z on a longer vDSP DFT (slower, same
+    /// tolerances). Agrees with `.double` to about 1e-5 of the correlation maximum, hence peak positions to about 1e-4 px
+    /// (docs/bragg.md has the measured gaps).
     case single
 }
 

Sources/CuboBragg/CuboBragg.docc/CuboBragg.md

@@ -10,9 +10,8 @@ correlation, the filter cascade, and parabolic or matrix-DFT (multicorr) subpixe
 
 Opt in to the faster Float32 path with `options.precision = .single` (``BraggPrecision``): vDSP convolution and complex DFT, a
 reusable scratch per concurrent worker, about ten times faster than Double at 256 × 256 and equal to it to about 1e-5 of the
-correlation maximum. It needs a detector size vDSP supports (64, 128, 192, 240, 256, ...; ask with
-``BraggDetector/supportsSinglePrecision(detectorRows:detectorColumns:)``) and is refused with ``CuboBraggError/invalidOption(_:)``
-otherwise. ``BraggDetector/detectWithFunnel(_:options:)`` returns the peaks of `detect` with a ``BraggFunnel``, the count of
+correlation maximum. It runs on every detector size: sizes vDSP's DFT takes (64, 128, 192, 240, 256, ...) directly, any other
+through Bluestein's chirp-z, slower and at the same tolerances. ``BraggDetector/detectWithFunnel(_:options:)`` returns the peaks of `detect` with a ``BraggFunnel``, the count of
 maxima after each stage of the selection, so a short or empty result can be explained; ``BraggDetector/fillNonFinite(_:rows:columns:)``
 is the median fill of non-finite pixels, public.
 

Sources/CuboBragg/CuboBragg.docc/Finding-Bragg-disks.md

@@ -28,8 +28,8 @@ let peaks = vectors[0, 0]                                                // the
 ```
 
 ``BraggDiskOptions`` defaults are py4DSTEM's. Choose the subpixel method with ``BraggSubpixel`` and the Float32 path with
-`options.precision = .single` (``BraggPrecision``). Float32 needs a detector size vDSP supports; ask first with
-``BraggDetector/supportsSinglePrecision(detectorRows:detectorColumns:)``.
+`options.precision = .single` (``BraggPrecision``). Float32 runs on every detector size; sizes vDSP's DFT takes
+(64, 128, 192, 240, 256, ...) are the fastest.
 
 ## When a result is short or empty
 

Sources/CuboBragg/SingleDetection.swift

@@ -80,9 +80,9 @@ final class SingleState: @unchecked Sendable {
     private var free: [SingleScratch] = []
 
     init(rows: Int, columns: Int, kernel: [Double]) throws {
-        guard let fft = FloatFFT2D(rows: rows, columns: columns) else {
-            throw CuboBraggError.invalidOption("precision .single needs a detector size vDSP's DFT supports on both axes (f·2ⁿ, f in 1, 3, 5, 15, ...); \(rows) × \(columns) is not one (use .double)")
-        }
+        guard rows > 0, columns > 0 else { throw CuboBraggError.shapeMismatch("the kernel is empty") }
+        // vDSP's DFT where it takes the length (the same calls and bits as `FloatFFT2D(rows:columns:)`), Bluestein otherwise.
+        let fft = FloatFFT2D(anyRows: rows, columns: columns)
         self.rows = rows; self.columns = columns; self.fft = fft
         let n = rows * columns
         let scratch = SingleScratch(count: n)

Tests/CuboBraggTests/SinglePrecisionTests.swift

@@ -12,7 +12,8 @@ import Testing
 
     func single(_ o: BraggDiskOptions) -> BraggDiskOptions { var s = o; s.precision = .single; return s }
 
-    /// Scenes vDSP's DFT supports (64 × 64); S2 (40 × 56) is refused.
+    /// Every scene runs in `.single`: the 64 × 64 ones on vDSP's DFT directly, S2 (40 × 56: vDSP has no DFT of length 56)
+    /// through Bluestein (CuboKit#60).
     func supported(_ c: [String: Any]) -> Bool {
         let s = fixture.scene(c["scene"] as! String)
         return BraggDetector.supportsSinglePrecision(detectorRows: s.rows, detectorColumns: s.columns)
@@ -34,7 +35,8 @@ import Testing
             strongPosition = max(strongPosition, g.strongPosition); strongIntensity = max(strongIntensity, g.strongIntensity); checked += 1
         }
         print("CuboBragg .single vs py4DSTEM: \(checked) cases, maximum position gap \(worstPosition) px (peaks above 1e-3 of the brightest: \(strongPosition)), maximum relative intensity gap \(worstIntensity) (those: \(strongIntensity))")
-        #expect(checked >= 24)
+        #expect(checked >= 29)   // S2 (40 × 56, Bluestein) included since CuboKit#60
+
     }
 
     /// The same inputs through `.double` and `.single`: peaks agree (py4DSTEM's cases and the app's variants).
@@ -54,7 +56,7 @@ import Testing
             strongPosition = max(strongPosition, g.strongPosition); strongIntensity = max(strongIntensity, g.strongIntensity); checked += 1
         }
         print("CuboBragg .single vs .double: \(checked) cases, maximum position gap \(worstPosition) px (peaks above 1e-3 of the brightest: \(strongPosition)), maximum relative intensity gap \(worstIntensity) (those: \(strongIntensity))")
-        #expect(checked >= 30)
+        #expect(checked >= 39)   // S2 (40 × 56, Bluestein) included since CuboKit#60
     }
 
     /// The correlation itself: the real map and the Fourier product, as a fraction of their largest magnitude.
@@ -80,35 +82,86 @@ import Testing
         #expect(worstFourier <= SingleTolerance.fourier, "Fourier gap \(worstFourier)")
     }
 
-    @Test func refusesUnsupportedSizesAndKeepsDouble() async throws {
-        let s2 = fixture.scene("S2")   // 40 × 56: vDSP has no DFT of length 56
-        #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 40, detectorColumns: 56))
-        #expect(BraggDetector.supportsSinglePrecision(detectorRows: 64, detectorColumns: 64))
-        #expect(BraggDetector.supportsSinglePrecision(detectorRows: 256, detectorColumns: 256))
-        #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 56, detectorColumns: 56))
-        #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 127, detectorColumns: 128))
-        #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 0, detectorColumns: 64))
-        var o = BraggDiskOptions(); o.precision = .single
-        #expect(throws: CuboBraggError.self) { try BraggDisks.detect(s2.patterns[0], kernel: s2.kernel, options: o) }
-        #expect(throws: CuboBraggError.self) { try BraggDetector(kernel: s2.kernel).detectWithFunnel(s2.patterns[0], options: o) }
-        #expect(throws: CuboBraggError.self) { try BraggDetector(kernel: s2.kernel).correlate(s2.patterns[0], options: o) }
-        do { _ = try BraggDisks.detect(s2.patterns[0], kernel: s2.kernel, options: o) } catch let e as CuboBraggError {
-            guard case .invalidOption(let message) = e else { Issue.record("wrong error \(e)"); return }
-            #expect(message.contains("40 × 56"))
+    /// CuboKit#60: `.single` runs on every positive detector size (vDSP's DFT where it takes the length, Bluestein otherwise);
+    /// `.double` stays the default and options are validated first, as in `.double`.
+    @Test func everySizeIsSupportedAndDoubleStaysDefault() throws {
+        for (rows, columns) in [(64, 64), (256, 256), (40, 56), (56, 56), (127, 128), (22, 21), (165, 40), (1, 1)] {
+            #expect(BraggDetector.supportsSinglePrecision(detectorRows: rows, detectorColumns: columns), "\(rows) × \(columns)")
         }
-        // The whole-cube pass refuses before reading a tile.
-        let shape = try CubeShape(scanRows: 1, scanColumns: 1, detectorRows: 40, detectorColumns: 56)
-        let cube = try MemoryCube(shape: shape, values: [Float](repeating: 0, count: 40 * 56))
-        await #expect(throws: CuboBraggError.self) { _ = try await BraggDisks.detectAll(in: cube, kernel: s2.kernel, options: o) }
-        // .double is untouched on that size, and is the default.
+        #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 0, detectorColumns: 64))
+        #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 64, detectorColumns: -1))
         #expect(BraggDiskOptions().precision == .double)
-        _ = try BraggDisks.detect(s2.patterns[0], kernel: s2.kernel)
-        // Options are validated first, as in .double.
-        o.upsampleFactor = 2
+        var o = BraggDiskOptions(); o.precision = .single; o.upsampleFactor = 2
         #expect(throws: CuboBraggError.self) { try BraggDisks.detect(fixture.scene("S1").patterns[0], kernel: fixture.scene("S1").kernel, options: o) }
     }
 
-    /// Non-square supported grids (the transposed spectrum must be undone correctly) and a Gaussian wider than the image
+    /// CuboKit#60: detectors vDSP's DFT cannot take on one or both axes (22 × 21, 165 × 40, S2's 40 × 56): the `.single`
+    /// correlation map and Fourier product agree with `.double` at the written float32 tolerances, the peaks of the S1 data
+    /// tiled onto the grid agree within a stated bound (below), and the whole-cube pass equals the per-pattern pass.
+    @Test func nonVDSPSizesAgreeWithDouble() async throws {
+        let s = fixture.scene("S1")
+        var worstMap = 0.0, worstFourier = 0.0, worstPosition = 0.0, worstIntensity = 0.0, peaksChecked = 0
+        for (rows, columns, sigma) in [(22, 21, 0.0), (22, 21, 1.5), (165, 40, 0.0), (165, 40, 2.0), (40, 56, 0.0)] {
+            #expect(!FloatFFT2D.supports(rows: rows, columns: columns), "\(rows) × \(columns) should need Bluestein")
+            let kernel = DetectorImage(detectorRows: rows, detectorColumns: columns,
+                                       values: (0..<(rows * columns)).map { s.kernel.values[$0 % 4096] }, path: .cpu)
+            let detector = try BraggDetector(kernel: kernel)
+            var o = BraggDiskOptions(); o.patternSigma = sigma
+            var os = o; os.precision = .single
+            var floats: [Float] = []
+            for pattern in s.patterns {
+                let values = (0..<(rows * columns)).map { Double(Float(pattern.values[$0 % 4096])) }
+                floats += values.map { Float($0) }
+                let image = DetectorImage(detectorRows: rows, detectorColumns: columns, values: values, path: .cpu)
+                let d = try detector.correlate(image, options: o)
+                let f = try detector.correlate(image, options: os)
+                let mapScale = d.real.max()!
+                let fourierScale = zip(d.fourierRe, d.fourierIm).map { hypot($0, $1) }.max()!
+                for i in d.real.indices {
+                    worstMap = max(worstMap, abs(f.real[i] - d.real[i]) / mapScale)
+                    worstFourier = max(worstFourier, hypot(f.fourierRe[i] - d.fourierRe[i], f.fourierIm[i] - d.fourierIm[i]) / fourierScale)
+                }
+                for subpixel in [BraggSubpixel.parabolic, .multicorr] {
+                    var p = BraggDiskOptions.mac4DSTEM; p.patternSigma = sigma; p.subpixel = subpixel
+                    p.edgeBoundary = 2; p.minimumPeakSpacing = 4
+                    var ps = p; ps.precision = .single
+                    let dp = try detector.detect(image, options: p)
+                    let fp = try detector.detect(image, options: ps)
+                    let id = "\(rows) × \(columns) σ \(sigma) \(subpixel)"
+                    #expect(dp.count == fp.count, "\(id): double \(dp.count) peaks, single \(fp.count)")
+                    guard dp.count == fp.count else { continue }
+                    // The tiled scene has broad, low-contrast maxima; there even vDSP's own sizes (160 × 40) exceed the strict
+                    // per-peak tier (measured 8.4e-6 relative intensity), so these peaks get a stated bound of 2e-4 px (the
+                    // app's own tolerance against py4DSTEM) and 2e-5 (measured here 9.9e-5 px, 8.0e-6). The real non-vDSP
+                    // scene S2 (40 × 56) is held to the strict tier in `singleAgreesWithDouble`.
+                    for (a, e) in zip(fp, dp) {
+                        let position = max(abs(a.row - e.row), abs(a.column - e.column))
+                        let intensity = abs(a.intensity - e.intensity) / max(abs(e.intensity), 1e-300)
+                        #expect(position <= 2e-4 && intensity <= 2e-5,
+                                "\(id): peak at \(e.row), \(e.column): position gap \(position), intensity gap \(intensity)")
+                        worstPosition = max(worstPosition, position); worstIntensity = max(worstIntensity, intensity)
+                    }
+                    peaksChecked += dp.count
+                }
+            }
+            // The whole-cube pass in `.single` on this size is the per-pattern pass.
+            let count = s.patterns.count
+            let shape = try CubeShape(scanRows: count, scanColumns: 1, detectorRows: rows, detectorColumns: columns)
+            let cube = try MemoryCube(shape: shape, values: floats)
+            var p = BraggDiskOptions.mac4DSTEM; p.patternSigma = sigma; p.precision = .single
+            let all = try await BraggDisks.detectAll(in: cube, kernel: kernel, options: p, tiling: .rows(2))
+            let n = rows * columns
+            for k in 0..<count {
+                let one = try floats[(k * n)..<((k + 1) * n)].withUnsafeBufferPointer { try detector.detect($0, options: p) }
+                #expect(all[k, 0] == one, "\(rows) × \(columns) position \(k)")
+            }
+        }
+        print("CuboBragg .single on non-vDSP sizes: map gap \(worstMap), Fourier gap \(worstFourier), \(peaksChecked) peaks, position gap \(worstPosition) px, relative intensity gap \(worstIntensity)")
+        #expect(worstMap <= SingleTolerance.correlation && worstFourier <= SingleTolerance.fourier, "map \(worstMap), Fourier \(worstFourier)")
+        #expect(peaksChecked > 0)
+    }
+
+    /// Non-square vDSP grids (the transposed spectrum must be undone correctly) and a Gaussian wider than the image
     /// (the reflection folds repeatedly): the correlation map and its Fourier product agree with `.double`.
     @Test func nonSquareAndWideFilter() throws {
         let s = fixture.scene("S1")

docs/bragg.md

@@ -53,10 +53,12 @@ hybrid power in `powf`), `max(Re, 0)` and the maxima search run in Float. The fe
 from the Float neighbourhood, the matrix-DFT upsampling on the Float complex correlation widened to Double. The reported intensity
 is the Float correlation, widened.
 
-- **Sizes.** The transform is the exact N-point DFT at the pattern's own size, so only detector sizes vDSP supports on both axes are
-  available: f·2ⁿ with f in 1, 3, 5, 15, ... (64, 128, 192, 240, 256, 320, 384, 480, 512, ...; not 56, 100, 127, 130).
-  `BraggDetector.supportsSinglePrecision(detectorRows:detectorColumns:)` asks vDSP itself. Any other size is refused with
-  `CuboBraggError.invalidOption` (naming the size), never computed another way; `.double` runs on every size.
+- **Sizes.** The transform is the exact N-point DFT at the pattern's own size, on every detector size (CuboCore's
+  `FloatFFT2D(anyRows:columns:)`). A length vDSP's DFT takes (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128, 192, 240, 256, 320,
+  384, 480, 512, ...) runs on it directly, with the same calls and bits as before; any other (56, 100, 127, 130, 165, ...)
+  runs as Bluestein's chirp-z through a vDSP DFT of the shortest length M ≥ 2N − 1 vDSP takes, slower but held to the same
+  tolerances (`SinglePrecisionTests.nonVDSPSizesAgreeWithDouble`: 22 × 21, 165 × 40, 40 × 56).
+  `BraggDetector.supportsSinglePrecision(detectorRows:detectorColumns:)` is true for every positive size.
 - **Scratch.** A detector keeps a pool of scratch buffers (six `rows × columns` Float planes); a call checks one out and returns
   it. `BraggDisks.detectAll` therefore allocates one per concurrent worker, not one per pattern, and with `.single` hands each
   Float pattern straight from the tile to the detector (no widening). The Float plan and kernel spectrum are built on the first

tools/consumer-check/Sources/consumer-check/main.swift

@@ -121,7 +121,8 @@ do {
           && BraggPrecision.allCases == [.double, .single] && singlePeaks.count == 2
           && abs(singlePeaks[0].row - doublePeaks[0].row) < 1e-3 && abs(singlePeaks[1].column - doublePeaks[1].column) < 1e-3)
     check("BraggDetector.supportsSinglePrecision", BraggDetector.supportsSinglePrecision(detectorRows: 48, detectorColumns: 48)
-          && !BraggDetector.supportsSinglePrecision(detectorRows: 56, detectorColumns: 40))
+          && BraggDetector.supportsSinglePrecision(detectorRows: 56, detectorColumns: 40)
+          && !BraggDetector.supportsSinglePrecision(detectorRows: 0, detectorColumns: 40))
     let plainPeaks = try detector.detect(image, options: braggOptions)
     check("BraggDetector.detectWithFunnel(DetectorImage)", doublePeaks == plainPeaks
           && doubleFunnel.acceptedCount == 2 && doubleFunnel.localMaximumCount >= 2 && singleFunnel.acceptedCount == 2)