Vandermonde interpolation in Swift and Python

README.md

Vandermonde Interpolation

Conduct interpolation using a Vandermonde matrix based on the articles in Apple Developer.

Differences

To allow compiliation without GUI, the original source is modified as follows.

• solveLinearSystem is declared without static and called without ViewController.

Files

• vandermonde.swift: prints five random values at x=0, 256, 512, 1023 and interpolated values at evenly space points in x between 0 and 1023.
• plot.py: draws generated and interpolated values
• vandermonde.py: Numpy translation

Compiliation and execution

% swiftc vandermonde.swift
% ./vandermonde > y.txt

or

% swift vandermode.swift > y.txt

Running Numpy version

% python vandermonde.py

Plotting

% python plot.py

References

• Finding an Interpolating Polynomial Using the Vandermonde Method
• Double-Precision Floating-Point Vectors
• LAPACK documentation dgels
• numpy.vander
• numpy.polyval

plot.py

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 1025, 256)
y = np.loadtxt("y.txt")

plt.rcParams["font.size"] = 18
fig, ax = plt.subplots(figsize=[10, 8], layout="constrained")
ax.scatter(x, y[0:x.size])
ax.plot(y[x.size:])
fig.savefig("vandermonde.png", dpi=300)
plt.show()import numpy as np

seed = 514
rng = np.random.default_rng(seed)
a = np.vander(np.array([0, 256, 512, 768, 1023]))
b = rng.uniform(-10, 10, 5)
c = np.linalg.solve(a, b)
x = np.arange(1024)
y = np.polyval(c, x)
y = np.concatenate([b, y])
np.savetxt("y.txt", y)

vandermonde.py

import numpy as np

seed = 514
rng = np.random.default_rng(seed)
a = np.vander(np.array([0, 256, 512, 768, 1023]))
b = rng.uniform(-10, 10, 5)
c = np.linalg.solve(a, b)
x = np.arange(1024)
y = np.polyval(c, x)
y = np.concatenate([b, y])
np.savetxt("y.txt", y)

vandermonde.swift

import simd
import Accelerate

let points: [simd_double2] = [
    simd_double2(   0, Double.random(in: -10...10)),
    simd_double2( 256, Double.random(in: -10...10)),
    simd_double2( 512, Double.random(in: -10...10)),
    simd_double2( 768, Double.random(in: -10...10)),
    simd_double2(1023, Double.random(in: -10...10)),
]

let exponents = (0 ..< points.count).map {
    return Double($0)
}

let vandermonde: [[Double]] = points.map { point in
    let bases = [Double](repeating: point.x, count: points.count)
    return vForce.pow(bases: bases, exponents: exponents)
}

public enum LAPACKError: Swift.Error {
    case internalError
    case parameterHasIllegalValue(parameterIndex: Int)
    case diagonalElementOfTriangularFactorIsZero(index: Int)
}

func solveLinearSystem(
    a: inout [Double], a_rowCount: Int, a_columnCount: Int,
    b: inout [Double], b_count: Int) throws {

    var info = Int32(0)

    // 1: Specify transpose.
    var trans = Int8("T".utf8.first!)

    // 2: Define constants.
    var m = __CLPK_integer(a_rowCount)
    var n = __CLPK_integer(a_columnCount)
    var lda = __CLPK_integer(a_rowCount)
    var nrhs = __CLPK_integer(1) // assumes `b` is a column matrix
    var ldb = __CLPK_integer(b_count)

    // 3: Workspace query.
    var workDimension = Double(0)
    var minusOne = Int32(-1)

    dgels_(&trans, &m, &n, &nrhs, &a, &lda, &b, &ldb, &workDimension, &minusOne, &info)

    if info != 0 {
        throw LAPACKError.internalError
    }

    // 4: Create workspace.
    var lwork = Int32(workDimension)
    var workspace = [Double](repeating: 0, count: Int(workDimension))

    // 5: Solve linear system.
    dgels_(&trans, &m, &n, &nrhs, &a, &lda, &b, &ldb, &workspace, &lwork, &info)

    if info < 0 {
        throw LAPACKError.parameterHasIllegalValue(parameterIndex: abs(Int(info)))
    } else if info > 0 {
        throw LAPACKError.diagonalElementOfTriangularFactorIsZero(index: Int(info))
    }
}

let coefficients: [Double] = {
    var a = vandermonde.flatMap{ $0 }
    var b = points.map{ $0.y }

    do {
        try solveLinearSystem(
            a: &a,
            a_rowCount: points.count,
            a_columnCount: points.count,
            b: &b,
            b_count: points.count)
    } catch {
        fatalError("Unable to solve linear system.")
    }

    vDSP.reverse(&b)

    return b
}()

let ramp = vDSP.ramp(withInitialValue:0, increment: Double(1), count: 1024)

let polynomialResult = vDSP.evaluatePolynomial(
    usingCoefficients: coefficients,
    withVariables: ramp)

for p in points {
    print(p.y)
}
for y in polynomialResult {
    print(y)
}