Skip to content

Commit 6982f92

Browse files
FlexMonkeystephentyrone
FlexMonkey
authored and
stephentyrone
committed
[Accelerate] [Quadrature] New Quadrature Overlay (#23127)
* [Accelerate] [Quadrature] New Quadrature Overlay A class that simplifies approximating the definite integral of a function. * Fixes in response to PR Review. Change `@_exported import Accelerate` to `import Accelerate`. Correct date in copyright comment. * Code Review Fixes. * Remove mutable integrator, simply set options in init(). * Remove mutable `maxIntervals` and `qagPointsPerInterval`. * Update tests. * Code Review Fixes. * Use standard library `Result`. * Implement `QAGPointsPerInterval` as a struct. * Tidy up Passing Integrand to `quadrature_integrate_function`. Pass the integrand closure directly to the `quadrature_integrate_function` initialiser rather than attaching as a property of the `Quadrature` class. * Make `Quadrature` a struct, and remove requirement for `integrand` to be escaping. * Code Review Changes * Add long-form integrator algorithm aliases (update tests to use these aliases). * Make quadrature error description public. * New tests for error description and `QAGPointsPerInterval`. * Vectorized Integrand for Quadrature Create an alternative implementation of `integrate` that accepts an integrand of type `(_ input: UnsafeBufferPointer<Double>, _ result: UnsafeMutableBufferPointer<Double>)`. * Vectorized Integrand for Quadrature Create an alternative implementation of `integrate` that accepts an integrand of type `(_ input: UnsafeBufferPointer<Double>, _ result: UnsafeMutableBufferPointer<Double>)`. * Vectorized Integrand for Quadrature Remove _ContiguousCollection - it's no longer used. * Delete ContiguousCollection.swift Not required. * Refactor tests.
1 parent 5d8bffd commit 6982f92

File tree

4 files changed

+667
-56
lines changed

4 files changed

+667
-56
lines changed

‎stdlib/public/Darwin/Accelerate/CMakeLists.txt

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,12 @@ cmake_minimum_required(VERSION 3.4.3)
22
include("../../../../cmake/modules/StandaloneOverlay.cmake")
33

44
add_swift_target_library(swiftAccelerate ${SWIFT_SDK_OVERLAY_LIBRARY_BUILD_TYPES} IS_SDK_OVERLAY
5-
BNNS.swift.gyb
5+
Quadrature.swift
6+
7+
"${SWIFT_SOURCE_DIR}/stdlib/linker-support/magic-symbols-for-install-name.c"
8+
9+
GYB_SOURCES
10+
BNNS.swift.gyb
611

712
SWIFT_COMPILE_FLAGS "${SWIFT_RUNTIME_SWIFT_COMPILE_FLAGS}"
813
LINK_FLAGS "${SWIFT_RUNTIME_SWIFT_LINK_FLAGS}"
@@ -13,6 +18,8 @@ add_swift_target_library(swiftAccelerate ${SWIFT_SDK_OVERLAY_LIBRARY_BUILD_TYPES
1318
SWIFT_MODULE_DEPENDS_TVOS Darwin CoreFoundation CoreGraphics Dispatch Foundation Metal ObjectiveC os # auto-updated
1419
SWIFT_MODULE_DEPENDS_WATCHOS Darwin CoreFoundation CoreGraphics Dispatch Foundation ObjectiveC os # auto-updated
1520

21+
FRAMEWORK_DEPENDS Accelerate
22+
1623
DEPLOYMENT_VERSION_OSX ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_OSX}
1724
DEPLOYMENT_VERSION_IOS ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_IOS}
1825
DEPLOYMENT_VERSION_TVOS ${SWIFTLIB_DEPLOYMENT_VERSION_SIMD_TVOS}

‎stdlib/public/Darwin/Accelerate/Quadrature.swift

Lines changed: 295 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,295 @@
1+
//===----------------------------------------------------------------------===//
2+
//
3+
// This source file is part of the Swift.org open source project
4+
//
5+
// Copyright (c) 2014 - 2019 Apple Inc. and the Swift project authors
6+
// Licensed under Apache License v2.0 with Runtime Library Exception
7+
//
8+
// See https://swift.org/LICENSE.txt for license information
9+
// See https://swift.org/CONTRIBUTORS.txt for the list of Swift project authors
10+
//
11+
//===----------------------------------------------------------------------===//
12+
13+
import Accelerate
14+
15+
// MARK: Quadrature
16+
17+
@available(iOS 9999, macOS 9999, tvOS 9999, watchOS 9999, *)
18+
/// A structure that approximates the definite integral of a function over a finite interval.
19+
///
20+
/// The following code is an example of using a `Quadrature` structure to calculate the
21+
/// area under a curve defined by `y = sqrt(radius * radius - pow(x - radius, 2))`:
22+
///
23+
///
24+
/// let quadrature = Quadrature(integrator: .qags(maxIntervals: 10),
25+
/// absoluteTolerance: 1.0e-8,
26+
/// relativeTolerance: 1.0e-2)
27+
///
28+
/// let result = quadrature.integrate(over: 0.0 ... 25.0) { x in
29+
/// let radius: Double = 12.5
30+
/// return sqrt(radius * radius - pow(x - radius, 2))
31+
/// }
32+
///
33+
/// switch result {
34+
/// case .success(let integralResult, let estimatedAbsoluteError):
35+
/// print("quadrature success:", integralResult,
36+
/// estimatedAbsoluteError)
37+
/// case .failure(let error):
38+
/// print("quadrature error:", error.errorDescription)
39+
/// }
40+
///
41+
/// Alternatively, you can integrate over a function that uses vectors for its
42+
/// source and destination. For example:
43+
///
44+
/// func vectorExp(x: UnsafeBufferPointer<Double>,
45+
/// y: UnsafeMutableBufferPointer<Double>) {
46+
/// let radius: Double = 12.5
47+
/// for i in 0 ..< x.count {
48+
/// y[i] = sqrt(radius * radius - pow(x[i] - radius, 2))
49+
/// }
50+
/// }
51+
///
52+
/// let vRresult = quadrature.integrate(over: 0.0 ... diameter,
53+
/// integrand: vectorExp)
54+
public struct Quadrature {
55+
56+
private var integrateOptions = quadrature_integrate_options()
57+
58+
/// Initializes and returns a quadrature instance.
59+
///
60+
/// - Parameter integrator: An enumeration specifying the integration algorithm and relevant properties.
61+
/// - Parameter absoluteTolerance: Requested absolute tolerance on the result.
62+
/// - Parameter relativeTolerance: Requested relative tolerance on the result.
63+
public init(integrator: Integrator,
64+
absoluteTolerance: Double = 1.0e-8,
65+
relativeTolerance: Double = 1.0e-2){
66+
67+
integrateOptions.abs_tolerance = absoluteTolerance
68+
integrateOptions.rel_tolerance = relativeTolerance
69+
70+
switch integrator {
71+
case .qng:
72+
integrateOptions.integrator = QUADRATURE_INTEGRATE_QNG
73+
case .qag(let pointsPerInterval, let maxIntervals):
74+
integrateOptions.integrator = QUADRATURE_INTEGRATE_QAG
75+
integrateOptions.qag_points_per_interval = pointsPerInterval.points
76+
integrateOptions.max_intervals = maxIntervals
77+
case .qags(let maxIntervals):
78+
integrateOptions.integrator = QUADRATURE_INTEGRATE_QAGS
79+
integrateOptions.max_intervals = maxIntervals
80+
}
81+
}
82+
83+
/// The quadrature instance's requested absolute tolerance on the result.
84+
public var absoluteTolerance: Double {
85+
set {
86+
integrateOptions.abs_tolerance = newValue
87+
}
88+
get {
89+
return integrateOptions.abs_tolerance
90+
}
91+
}
92+
93+
/// The quadrature instance's requested relative tolerance on the result.
94+
public var relativeTolerance: Double {
95+
set {
96+
integrateOptions.rel_tolerance = newValue
97+
}
98+
get {
99+
return integrateOptions.rel_tolerance
100+
}
101+
}
102+
103+
/// Performs the integration over the supplied function.
104+
///
105+
/// - Parameter interval: The lower and upper bounds of the integration interval.
106+
/// - 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.
107+
public func integrate(over interval: ClosedRange<Double>,
108+
integrand: (_ input: UnsafeBufferPointer<Double>, _ result: UnsafeMutableBufferPointer<Double>) -> ()) ->
109+
Result<(integralResult: Double, estimatedAbsoluteError: Double), Error>{
110+
111+
var status = QUADRATURE_SUCCESS
112+
var estimatedAbsoluteError: Double = 0
113+
var result: Double = 0
114+
115+
var callback: quadrature_integrate_function!
116+
117+
withoutActuallyEscaping(integrand) {escapableIntegrand in
118+
withUnsafePointer(to: escapableIntegrand) {
119+
let integrandPointer = UnsafeMutableRawPointer(mutating: $0)
120+
121+
callback = quadrature_integrate_function(
122+
fun: { (arg: UnsafeMutableRawPointer?,
123+
n: Int,
124+
x: UnsafePointer<Double>,
125+
y: UnsafeMutablePointer<Double>
126+
) in
127+
128+
guard let integrand = arg?.load(as: ((UnsafeBufferPointer<Double>, UnsafeMutableBufferPointer<Double>) ->()).self) else {
129+
return
130+
}
131+
132+
integrand(UnsafeBufferPointer(start: x, count: n),
133+
UnsafeMutableBufferPointer(start: y, count: n))
134+
},
135+
fun_arg: integrandPointer)
136+
137+
withUnsafePointer(to: self.integrateOptions) { options in
138+
result = quadrature_integrate(&callback,
139+
interval.lowerBound,
140+
interval.upperBound,
141+
options,
142+
&status,
143+
&estimatedAbsoluteError,
144+
0,
145+
nil)
146+
}
147+
}
148+
}
149+
150+
if status == QUADRATURE_SUCCESS {
151+
return .success((integralResult: result,
152+
estimatedAbsoluteError: estimatedAbsoluteError))
153+
} else {
154+
return .failure(Error(quadratureStatus: status))
155+
}
156+
}
157+
158+
/// Performs the integration over the supplied function.
159+
///
160+
/// - Parameter interval: The lower and upper bounds of the integration interval.
161+
/// - 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.
162+
public func integrate(over interval: ClosedRange<Double>,
163+
integrand: (Double) -> Double) ->
164+
Result<(integralResult: Double, estimatedAbsoluteError: Double), Error> {
165+
166+
var status = QUADRATURE_SUCCESS
167+
var estimatedAbsoluteError: Double = 0
168+
var result: Double = 0
169+
170+
var callback: quadrature_integrate_function!
171+
172+
withoutActuallyEscaping(integrand) {escapableIntegrand in
173+
withUnsafePointer(to: escapableIntegrand) {
174+
let integrandPointer = UnsafeMutableRawPointer(mutating: $0)
175+
176+
callback = quadrature_integrate_function(
177+
fun: { (arg: UnsafeMutableRawPointer?,
178+
n: Int,
179+
x: UnsafePointer<Double>,
180+
y: UnsafeMutablePointer<Double>
181+
) in
182+
183+
guard let integrand = arg?.load(as: ((Double) -> Double).self) else {
184+
return
185+
}
186+
187+
(0 ..< n).forEach { i in
188+
y[i] = integrand(x[i])
189+
}
190+
},
191+
fun_arg: integrandPointer)
192+
193+
withUnsafePointer(to: self.integrateOptions) { options in
194+
result = quadrature_integrate(&callback,
195+
interval.lowerBound,
196+
interval.upperBound,
197+
options,
198+
&status,
199+
&estimatedAbsoluteError,
200+
0,
201+
nil)
202+
}
203+
}
204+
}
205+
206+
if status == QUADRATURE_SUCCESS {
207+
return .success((integralResult: result,
208+
estimatedAbsoluteError: estimatedAbsoluteError))
209+
} else {
210+
return .failure(Error(quadratureStatus: status))
211+
}
212+
}
213+
214+
public enum Integrator {
215+
/// Simple non-adaptive automatic integrator using Gauss-Kronrod-Patterson quadrature coefficients.
216+
/// Evaluates 21, or 43, or 87 points in the interval until the requested accuracy is reached.
217+
case qng
218+
219+
/// Simple non-adaptive automatic integrator using Gauss-Kronrod-Patterson quadrature coefficients.
220+
/// Evaluates 21, or 43, or 87 points in the interval until the requested accuracy is reached.
221+
public static let nonAdaptive = Integrator.qng
222+
223+
/// Simple globally adaptive integrator.
224+
/// Allows selection of the number of Gauss-Kronrod points used in each subinterval, and the max number of subintervals.
225+
case qag(pointsPerInterval: QAGPointsPerInterval, maxIntervals: Int)
226+
227+
/// Simple globally adaptive integrator.
228+
/// Allows selection of the number of Gauss-Kronrod points used in each subinterval, and the max number of subintervals.
229+
public static func adaptive(pointsPerInterval: QAGPointsPerInterval, maxIntervals: Int) -> Integrator {
230+
return Integrator.qag(pointsPerInterval: pointsPerInterval, maxIntervals: maxIntervals)
231+
}
232+
233+
/// 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.
234+
/// 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.
235+
case qags(maxIntervals: Int)
236+
237+
/// 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.
238+
/// 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.
239+
public static func adaptiveWithSingularities(maxIntervals: Int) -> Integrator {
240+
return Integrator.qags(maxIntervals: maxIntervals)
241+
}
242+
}
243+
244+
public struct QAGPointsPerInterval {
245+
public let points: Int
246+
private init(points: Int) { self.points = points }
247+
248+
public static let fifteen = QAGPointsPerInterval(points: 15)
249+
public static let twentyOne = QAGPointsPerInterval(points: 21)
250+
public static let thirtyOne = QAGPointsPerInterval(points: 31)
251+
public static let fortyOne = QAGPointsPerInterval(points: 41)
252+
public static let fiftyOne = QAGPointsPerInterval(points: 51)
253+
public static let sixtyOne = QAGPointsPerInterval(points: 61)
254+
}
255+
256+
public enum Error: Swift.Error {
257+
case generic
258+
case invalidArgument
259+
case `internal`
260+
case integrateMaxEval
261+
case badIntegrandBehaviour
262+
263+
public init(quadratureStatus: quadrature_status) {
264+
switch quadratureStatus {
265+
case QUADRATURE_ERROR:
266+
self = .generic
267+
case QUADRATURE_INVALID_ARG_ERROR:
268+
self = .invalidArgument
269+
case QUADRATURE_INTERNAL_ERROR:
270+
self = .internal
271+
case QUADRATURE_INTEGRATE_MAX_EVAL_ERROR:
272+
self = .integrateMaxEval
273+
case QUADRATURE_INTEGRATE_BAD_BEHAVIOUR_ERROR:
274+
self = .badIntegrandBehaviour
275+
default:
276+
self = .internal
277+
}
278+
}
279+
280+
public var errorDescription: String {
281+
switch self {
282+
case .generic:
283+
return "Generic error."
284+
case .invalidArgument:
285+
return "Invalid Argument."
286+
case .internal:
287+
return "This is a bug in the Quadrature code, please file a bug report."
288+
case .integrateMaxEval:
289+
return "The requested accuracy limit could not be reached with the allowed number of evals/subdivisions."
290+
case .badIntegrandBehaviour:
291+
return "Extremely bad integrand behaviour, or excessive roundoff error occurs at some points of the integration interval."
292+
}
293+
}
294+
}
295+
}

0 commit comments

Comments
 (0)