[Accelerate] [Quadrature] New Quadrature Overlay

FlexMonkey · FlexMonkey · commit 0ab641764a99 · 2019-03-06T14:04:05.000Z
A class that simplifies approximating the definite integral of a function.
diff --git a/stdlib/public/Darwin/Accelerate/CMakeLists.txt b/stdlib/public/Darwin/Accelerate/CMakeLists.txt
@@ -2,19 +2,22 @@ cmake_minimum_required(VERSION 3.4.3)
 include("../../../../cmake/modules/StandaloneOverlay.cmake")
 
 add_swift_target_library(swiftAccelerate ${SWIFT_SDK_OVERLAY_LIBRARY_BUILD_TYPES} IS_SDK_OVERLAY
-  BNNS.swift.gyb
+BNNS.swift.gyb
+Quadrature.swift
 
-  SWIFT_COMPILE_FLAGS "${SWIFT_RUNTIME_SWIFT_COMPILE_FLAGS}"
-  LINK_FLAGS "${SWIFT_RUNTIME_SWIFT_LINK_FLAGS}"
-  TARGET_SDKS OSX IOS IOS_SIMULATOR TVOS TVOS_SIMULATOR WATCHOS WATCHOS_SIMULATOR
-  SWIFT_MODULE_DEPENDS_OSX Darwin CoreFoundation CoreGraphics Dispatch Foundation IOKit Metal ObjectiveC XPC # auto-updated
-    os
-  SWIFT_MODULE_DEPENDS_IOS Darwin CoreFoundation CoreGraphics Dispatch Foundation Metal ObjectiveC os # auto-updated
-  SWIFT_MODULE_DEPENDS_TVOS Darwin CoreFoundation CoreGraphics Dispatch Foundation Metal ObjectiveC os # auto-updated
-  SWIFT_MODULE_DEPENDS_WATCHOS Darwin CoreFoundation CoreGraphics Dispatch Foundation ObjectiveC os # auto-updated
+SWIFT_COMPILE_FLAGS "${SWIFT_RUNTIME_SWIFT_COMPILE_FLAGS}"
+LINK_FLAGS "${SWIFT_RUNTIME_SWIFT_LINK_FLAGS}"
+TARGET_SDKS OSX IOS IOS_SIMULATOR TVOS TVOS_SIMULATOR WATCHOS WATCHOS_SIMULATOR
+SWIFT_MODULE_DEPENDS_OSX Darwin CoreFoundation CoreGraphics Dispatch Foundation IOKit Metal ObjectiveC XPC # auto-updated
+os
+SWIFT_MODULE_DEPENDS_IOS Darwin CoreFoundation CoreGraphics Dispatch Foundation Metal ObjectiveC os # auto-updated
+SWIFT_MODULE_DEPENDS_TVOS Darwin CoreFoundation CoreGraphics Dispatch Foundation Metal ObjectiveC os # auto-updated
+SWIFT_MODULE_DEPENDS_WATCHOS Darwin CoreFoundation CoreGraphics Dispatch Foundation ObjectiveC os # auto-updated
 
-  DEPLOYMENT_VERSION_OSX ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_OSX}
-  DEPLOYMENT_VERSION_IOS ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_IOS}
-  DEPLOYMENT_VERSION_TVOS ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_TVOS}
-  DEPLOYMENT_VERSION_WATCHOS ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_WATCHOS}
+FRAMEWORK_DEPENDS Accelerate
+
+DEPLOYMENT_VERSION_OSX ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_OSX}
+DEPLOYMENT_VERSION_IOS ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_IOS}
+DEPLOYMENT_VERSION_TVOS ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_TVOS}
+DEPLOYMENT_VERSION_WATCHOS ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_WATCHOS}
 )
diff --git a/stdlib/public/Darwin/Accelerate/Quadrature.swift b/stdlib/public/Darwin/Accelerate/Quadrature.swift
@@ -0,0 +1,258 @@
+//===----------------------------------------------------------------------===//
+//
+// This source file is part of the Swift.org open source project
+//
+// Copyright (c) 2014 - 2017 Apple Inc. and the Swift project authors
+// Licensed under Apache License v2.0 with Runtime Library Exception
+//
+// See https://swift.org/LICENSE.txt for license information
+// See https://swift.org/CONTRIBUTORS.txt for the list of Swift project authors
+//
+//===----------------------------------------------------------------------===//
+
+@_exported import Accelerate
+
+// MARK: Quadrature
+
+@available(iOS 13.0, macOS 10.15, tvOS 13.0, watchOS 6.0, *)
+/// A class that approximates the definite integral of a function over a finite interval.
+///
+/// The following code is an example of using a `Quadrature` instance to calculate the
+/// area under a curve defined by `y = sqrt(radius * radius - pow(x - radius, 2))`:
+///
+///
+///     let quadrature = Quadrature(integrator: .qags(maxIntervals: 10),
+///                                 absoluteTolerance: 1.0e-8,
+///                                 relativeTolerance: 1.0e-2)
+///
+///     let result = quadrature.integrate(over: 0.0 ... 25.0) { x in
+///         let radius: Double = 12.5
+///         return sqrt(radius * radius - pow(x - radius, 2))
+///     }
+///
+///     switch result {
+///     case .success(let integralResult, let estimatedAbsoluteError):
+///         print("quadrature success:", integralResult,
+///                                      estimatedAbsoluteError)
+///     case .failure(let error):
+///         print("quadrature error:", error.errorDescription)
+///     }
+///
+public class Quadrature {
+    
+    private var integrateOptions = quadrature_integrate_options()
+    private var integrand: ((Double) -> Double)!
+    
+    
+    /// Initializes and returns a quadrature instance.
+    ///
+    /// - Parameter integrator: An enumeration specifying the integration algorithm and relevant properties.
+    /// - Parameter absoluteTolerance: Requested absolute tolerance on the result.
+    /// - Parameter relativeTolerance: Requested relative tolerance on the result.
+    public init(integrator: Integrator,
+                absoluteTolerance: Double = 1.0e-8,
+                relativeTolerance: Double = 1.0e-2){
+        
+        self.integrator = integrator
+        integrateOptions.abs_tolerance = absoluteTolerance
+        integrateOptions.rel_tolerance = relativeTolerance
+    }
+    
+    /// The quadrature instance's integration algorithm and relevant properties.
+    public var integrator: Integrator {
+        set {
+            switch newValue {
+            case .qng:
+                integrateOptions.integrator = QUADRATURE_INTEGRATE_QNG
+            case .qag(let pointsPerInterval, let maxIntervals):
+                integrateOptions.integrator = QUADRATURE_INTEGRATE_QAG
+                integrateOptions.qag_points_per_interval = pointsPerInterval.rawValue
+                integrateOptions.max_intervals = maxIntervals
+            case .qags(let maxIntervals):
+                integrateOptions.integrator = QUADRATURE_INTEGRATE_QAGS
+                integrateOptions.max_intervals = maxIntervals
+            }
+        }
+        get {
+            let integrator: Integrator!
+            
+            switch integrateOptions.integrator {
+            case QUADRATURE_INTEGRATE_QNG:
+                integrator = .qng
+            case QUADRATURE_INTEGRATE_QAG:
+                integrator = .qag(pointsPerInterval: QAGPointsPerInterval(rawValue: integrateOptions.qag_points_per_interval)!,
+                                  maxIntervals: integrateOptions.max_intervals)
+            case QUADRATURE_INTEGRATE_QAGS:
+                integrator = .qags(maxIntervals: integrateOptions.max_intervals)
+            default:
+                integrator = nil
+            }
+            
+            return integrator
+        }
+    }
+    
+    /// The quadrature instance's requested absolute tolerance on the result.
+    public var absoluteTolerance: Double {
+        set {
+            integrateOptions.abs_tolerance = newValue
+        }
+        get {
+            return integrateOptions.abs_tolerance
+        }
+    }
+    
+    /// The quadrature instance's requested relative tolerance on the result.
+    public var relativeTolerance: Double {
+        set {
+            integrateOptions.rel_tolerance = newValue
+        }
+        get {
+            return integrateOptions.rel_tolerance
+        }
+    }
+    
+    /// The maximum number of subintervals in the subdivision used by QAG and QAGS integrators.
+    public var maxIntervals: Int {
+        set {
+            integrateOptions.max_intervals = newValue
+        }
+        get {
+            return integrateOptions.max_intervals
+        }
+    }
+    
+    /// Number of points per subinterval. Used by the QAG integrator only; other integrators ignore this value.
+    public var qagPointsPerInterval: QAGPointsPerInterval {
+        set {
+            integrateOptions.qag_points_per_interval = newValue.rawValue
+        }
+        get {
+            return QAGPointsPerInterval(rawValue: integrateOptions.qag_points_per_interval) ?? .default
+        }
+    }
+    
+    /// Performs the integration over the supplied function.
+    ///
+    /// - Parameter interval: The lower and upper bounds of the integration interval.
+    /// - Parameter integrand: The function to integrate. The input value is `x` that's within the interval over which the integrand is being integrated, and the output value is the corresponding value `y = integrand(x)` at those points.
+    public func integrate(over interval: ClosedRange<Double>,
+                          integrand: @escaping (Double) -> Double) -> Result {
+        
+        self.integrand = integrand
+        
+        let byteCount = MemoryLayout<Quadrature>.size
+        let integrandPointer = UnsafeMutableRawPointer.allocate(byteCount: byteCount,
+                                                                alignment: 1)
+        integrandPointer.storeBytes(of: self,
+                                    as: Quadrature.self)
+        
+        var callback = quadrature_integrate_function(
+            fun: { (arg: UnsafeMutableRawPointer?,
+                n: Int,
+                x: UnsafePointer<Double>,
+                y: UnsafeMutablePointer<Double>
+                ) in
+                
+                guard let quadrature = arg?.load(as: Quadrature.self) else {
+                    return
+                }
+                
+                (0 ..< n).forEach { i in
+                    y[i] = quadrature.integrand(x[i])
+                }
+        },
+            fun_arg: integrandPointer)
+        
+        var status = QUADRATURE_SUCCESS
+        var estimatedAbsoluteError: Double = 0
+        var result: Double = 0
+        
+        withUnsafePointer(to: integrateOptions) { options in
+            result = quadrature_integrate(&callback,
+                                          interval.lowerBound,
+                                          interval.upperBound,
+                                          options,
+                                          &status,
+                                          &estimatedAbsoluteError,
+                                          0,
+                                          nil)
+        }
+        
+        if status == QUADRATURE_SUCCESS {
+            return .success(integralResult: result,
+                            estimatedAbsoluteError: estimatedAbsoluteError)
+        } else {
+            return .failure(error: Error(quadratureStatus: status))
+        }
+    }
+    
+    public enum Integrator {
+        /// Simple non-adaptive automatic integrator using Gauss-Kronrod-Patterson quadrature coefficients.
+        /// Evaluates 21, or 43, or 87 points in the interval until the requested accuracy is reached.
+        case qng
+        
+        /// Simple globally adaptive integrator.
+        /// Allows selection of the number of Gauss-Kronrod points used in each subinterval, and the max number of subintervals.
+        case qag(pointsPerInterval: QAGPointsPerInterval, maxIntervals: Int)
+        
+        /// Global adaptive quadrature based on 21-point or 15-point (if at least one bound is infinite) Gauss–Kronrod quadrature within each subinterval, with acceleration by Peter Wynn's epsilon algorithm.
+        /// If at least one of the interval bounds is infinite, this is equivalent to the QUADPACK QAGI routine. Otherwise, this is equivalent to the QUADPACK QAGS routine.
+        case qags(maxIntervals: Int)
+    }
+    
+    public enum QAGPointsPerInterval: Int {
+        case `default` = 0
+        case fifteen = 15
+        case twentyOne = 21
+        case thirtyOne = 31
+        case fortyOne = 41
+        case fiftyOne = 51
+        case sixtyOne = 61
+    }
+    
+    public enum Result {
+        case success (integralResult: Double, estimatedAbsoluteError: Double)
+        case failure (error: Error)
+    }
+    
+    public enum Error: Swift.Error {
+        case generic
+        case invalidArgument
+        case `internal`
+        case integrateMaxEval
+        case badIntegrandBehaviour
+        
+        init(quadratureStatus: quadrature_status) {
+            switch quadratureStatus {
+            case QUADRATURE_ERROR:
+                self = .generic
+            case QUADRATURE_INVALID_ARG_ERROR:
+                self = .invalidArgument
+            case QUADRATURE_INTERNAL_ERROR:
+                self = .internal
+            case QUADRATURE_INTEGRATE_MAX_EVAL_ERROR:
+                self = .integrateMaxEval
+            case QUADRATURE_INTEGRATE_BAD_BEHAVIOUR_ERROR:
+                self = .badIntegrandBehaviour
+            default:
+                self = .internal
+            }
+        }
+        
+        var errorDescription: String {
+            switch self {
+            case .generic:
+                return "Generic error"
+            case .invalidArgument:
+                return "Invalid Argument"
+            case .internal:
+                return "This is a bug in the Quadrature code, please file a bug report."
+            case .integrateMaxEval:
+                return "The requested accuracy limit could not be reached with the allowed number of evals/subdivisions."
+            case .badIntegrandBehaviour:
+                return "Extremely bad integrand behaviour, or excessive roundoff error occurs at some points of the integration interval."
+            }
+        }
+    }
+}
diff --git a/test/stdlib/Accelerate.swift b/test/stdlib/Accelerate.swift
