Helo all,
Currently, I'm working on an iOS app that performs measurement and shows the results to the user in a graph. I use a Savitzky-Golay filter to filter out noise, so that the graph is nice and smooth. However, the code that calculates the Savitzky-Golay coefficients using sparse matrices crashes sometimes, throwing an EXC_BAD_ACCESS. I tried to find out what the problem is by turning on Address Sanitizer and Thread Sanitizer, but, for some reason, the bad access exception isn't thrown when either of these is on. What else could I try to trace back the problem?
Thanks in advance, CaS
To reproduce the error, run the following:
import SwiftUI
import Accelerate
struct ContentView: View {
var body: some View {
VStack {
Button("Try", action: test)
}
.padding()
}
func test() {
for windowLength in 3...100 {
let coeffs = SavitzkyGolay.coefficients(windowLength: windowLength, polynomialOrder: 2)
print(coeffs)
}
}
}
class SavitzkyGolay {
static func coefficients(windowLength: Int, polynomialOrder: Int, derivativeOrder: Int = 0, delta: Int = 1) -> [Double] {
let (halfWindow, remainder) = windowLength.quotientAndRemainder(dividingBy: 2)
var pos = Double(halfWindow)
if remainder == 0 {
pos -= 0.5
}
let X = [Double](stride(from: Double(windowLength) - pos - 1, through: -pos, by: -1))
let P = [Double](stride(from: 0, through: Double(polynomialOrder), by: 1))
let A = P.map { exponent in
X.map {
pow($0, exponent)
}
}
var B = [Double](repeating: 0, count: polynomialOrder + 1)
B[derivativeOrder] = Double(factorial(derivativeOrder)) / pow(Double(delta), Double(derivativeOrder))
return leastSquaresSolution(A: A, B: B)
}
static func leastSquaresSolution(A: [[Double]], B: [Double]) -> [Double] {
let sparseA = A.sparseMatrix()
var sparseAValuesCopy = sparseA.values
var xValues = [Double](repeating: 0, count: A.transpose().count)
var bValues = B
sparseAValuesCopy.withUnsafeMutableBufferPointer { valuesPtr in
let a = SparseMatrix_Double(
structure: sparseA.structure,
data: valuesPtr.baseAddress!
)
bValues.withUnsafeMutableBufferPointer { bPtr in
xValues.withUnsafeMutableBufferPointer { xPtr in
let b = DenseVector_Double(
count: Int32(B.count),
data: bPtr.baseAddress!
)
let x = DenseVector_Double(
count: Int32(A.transpose().count),
data: xPtr.baseAddress!
)
#warning("EXC_BAD_ACCESS is thrown below")
print("This code is executed...")
let status = SparseSolve(SparseLSMR(), a, b, x, SparsePreconditionerDiagScaling)
print("...but, if an EXC_BAD_ACCESS is thrown, this code isn't")
if status != SparseIterativeConverged {
fatalError("Failed to converge. Returned with error \(status).")
}
}
}
}
return xValues
}
}
func factorial(_ n: Int) -> Int {
n < 2 ? 1 : n * factorial(n - 1)
}
extension Array where Element == [Double] {
func sparseMatrix() -> (structure: SparseMatrixStructure, values: [Double]) {
let columns = self.transpose()
var rowIndices: [Int32] = columns.map { column in
column.indices.compactMap { indexInColumn in
if column[indexInColumn] != 0 {
return Int32(indexInColumn)
}
return nil
}
}.reduce([], +)
let sparseColumns = columns.map { column in
column.compactMap {
if $0 != 0 {
return $0
}
return nil
}
}
var counter = 0
var columnStarts = [Int]()
for sparseColumn in sparseColumns {
columnStarts.append(counter)
counter += sparseColumn.count
}
let reducedSparseColumns = sparseColumns.reduce([], +)
columnStarts.append(reducedSparseColumns.count)
let structure: SparseMatrixStructure = rowIndices.withUnsafeMutableBufferPointer { rowIndicesPtr in
columnStarts.withUnsafeMutableBufferPointer { columnStartsPtr in
let attributes = SparseAttributes_t()
return SparseMatrixStructure(
rowCount: Int32(self.count),
columnCount: Int32(columns.count),
columnStarts: columnStartsPtr.baseAddress!,
rowIndices: rowIndicesPtr.baseAddress!,
attributes: attributes,
blockSize: 1
)
}
}
return (structure, reducedSparseColumns)
}
func transpose() -> Self {
let columns = self.count
let rows = self.reduce(0) { Swift.max($0, $1.count) }
return (0 ..< rows).reduce(into: []) { result, row in
result.append((0 ..< columns).reduce(into: []) { result, column in
result.append(row < self[column].count ? self[column][row] : 0)
})
}
}
}
Found an easier solution: using LAPACK's dgels function. Thanks for your help anyway, @Claude31 !
