How to contribute - have some small numpy / scipy vDSP implementations

[vade]

Hi there

Ive got a few minor implementations of some numpy / scipy functions on vectors in vDSP along with associated XCTests that I think would make good contributions to Matft.

I'm curious how to best / properly implement these into Matft, as it seems like having a single home for them makes sense, and your code seems very well organized.

I'm not entirely sure of the code structure / best place to implement the logic in a generic way that leverages Matfts existing code base.

Do you have any suggestions?

numpy.allclose() as an extension to Array (without NaN equality)

func allCloseTo(array: [Float], rtol: Float = 1e-5, atol: Float = 1e-8) -> Bool
    {
        precondition(self.count == array.count, "Arrays must have same size")
        
        let absDiff = vDSP.absolute( vDSP.subtract(self, array) )
            
        let maxAbsDiff = vDSP.maximum(absDiff)
        
        let scaledTol = Swift.max(atol, rtol * vDSP.maximum( vDSP.absolute(self) + vDSP.absolute(array) ) )
    
        return maxAbsDiff <= scaledTol
    }

scipy.spatial.distance.cosine

func CosineDistance(_ v1: [Float], _ v2: [Float]) -> Float
{
    precondition(v1.count == v2.count, "Arrays must have same size")

    var dotProduct: Float = 0.0
    var v1Norm: Float = 0.0
    var v2Norm: Float = 0.0
    
    let n = vDSP_Length(v1.count)
    
    // Calculate dot product of v1 and v2
    vDSP_dotpr(v1, 1, v2, 1, &dotProduct, n)
    
    // Calculate the Euclidean norm of v1
    vDSP_svesq(v1, 1, &v1Norm, n)
    v1Norm = sqrt(v1Norm)
    
    // Calculate the Euclidean norm of v2
    vDSP_svesq(v2, 1, &v2Norm, n)
    v2Norm = sqrt(v2Norm)
    
    // Calculate cosine distance
    let distance = 1.0 - (dotProduct / (v1Norm * v2Norm))
    
    return distance
}

and scipy.ndimage.gaussian_filter_1d as array extensions allowing one to cache the computed gaussian kernel.

Note I only really implement the default padding of reflect so far.

 static func generateGaussianKernel(sigma:Float, truncate:Float = 4.0) -> [Float]
    {
        let radius:Int = Int( ceil(truncate * sigma) )
        let sigma2 = sigma * sigma
        let x:[Float] = Array<Int>( ( -radius ... radius  ) ).map { Float( $0 ) }
        let x2 = vForce.pow(bases: x, exponents: [Float](repeating: 2.0, count: x.count) )
        let y = vDSP.multiply(-0.5 / sigma2, x2)
        let phi_x = vForce.exp(y)
        return vDSP.divide(phi_x, vDSP.sum(phi_x))
    }
    
    enum PaddingMode {
        case reflect
        case edge
    }

    private func padInputArray(_ input: [Float], sigma: Float, truncate: Float, paddingMode: PaddingMode) -> [Float] {
        var paddedInput = [Float]()
        let windowSize = Int(2.0 * sigma * truncate + 1.0)
        let padSize = Swift.max(windowSize - input.count, 0)

        if padSize > 0
        {
            switch (paddingMode)
            {
                case .reflect:
                                
                var paddingStart:[Float]
                var paddingEnd:[Float]
                
                // If we pad less than our input arrays count, we select what we need from the input array
                // This wont be a 'full' pad, as we wont have all items in the array
                if padSize <= input.count
                {
                    paddingStart = Array<Float>( input[ 0 ..< Int(padSize)].reversed() )
                    paddingEnd = Array<Float>( input[ input.count - Int(padSize) ..< input.count].reversed() )
                }
                // Otherwise, we repeat reflection until we accrue pad size
                else
                {
                    paddingStart = input.reversed()
                    paddingEnd = paddingStart
                    
                    while paddingStart.count <= padSize
                    {
                        paddingStart.insert(contentsOf: paddingStart.reversed(), at: 0)
                        paddingEnd.append(contentsOf: paddingEnd.reversed())
                        
                        paddingStart = paddingStart.reversed()
                        paddingEnd = paddingEnd.reversed()
                    }
                    
                    paddingStart = Array<Float>( paddingStart.suffix( Int(sigma * truncate)  ) )
                    paddingEnd = Array<Float>( paddingEnd.prefix( Int(sigma * truncate) ) )
                }
                                
                paddedInput.append(contentsOf: paddingStart)
                paddedInput.append(contentsOf: input)
                paddedInput.append(contentsOf: paddingEnd)
                
                break

            case .edge:
                let edge = input.first ?? 0.0
                paddedInput = Array(repeating: edge, count: padSize) + input + Array(repeating: edge, count: padSize)
                
            }
            return paddedInput
        }
        
        return input

    }
    
    // Make sure your Sigma and Truncate values match above:
    func gaussianFilter1D(kernel:[Float], sigma:Float, truncate:Float = 4.0, paddingMode:PaddingMode = .reflect) -> [Float]
    {
        let paddedInput = self.padInputArray(self, sigma:sigma, truncate:truncate, paddingMode:paddingMode)
        
        var output = [Float](repeating: 0.0, count: self.count)

        vDSP.convolve(paddedInput, withKernel: kernel, result: &output)

        // Technically is this needed, our sum is always 1 ?
//        vDSP.divide(output, sigma, result: &output)
//        let sum = vDSP.sum(kernel)
//        vDSP.multiply(sum, output, result: &output)

        return output
    }

[jjjkkkjjj]

@vade
I apologize for very late response. (I didn't notice your message...)
First, thank you for your willing to contribute to Matft!
My implementation is following rule, (but not decide strictly)

1. First, when i want to use Accelerate framework, such like vDSP_** or cblas_**, put codes regarding to Accelerate into Sources/Matft/library/*.swift. I use the vDSP_vneg(D) (known as negative operation) as the example following this explanation. So, because i want to use vDSP_vneg(D) function, put codes into Sources/Matft/library/vDSP.swift.
2. Next, to use the vDSP function, I split codes in vDSP.swift into 3 blocks.

• Create generic typealias of vDSP to pass the vDSP_vneg and vDSP_vnegD. The naming rule is vDSP_**_func. As you can know, there are two functions for float and double respectively in Accelerate framework . And also, Matft uses float and double into the core (StoredType).

internal typealias vDSP_convert_func<T, U> = (UnsafePointer<T>, vDSP_Stride, UnsafeMutablePointer<U>, vDSP_Stride, vDSP_Length) -> Void

Original code

Note: vDSP_vneg's arguments are same as some conversion functions (e.g. vDSP_vdpsp) when i use two generic types (T and U)

• Create wrapper function to use the above typealias. The naming rule is wrap_vDSP_**. In this code, I implement the type conversion mainly such like Int into vDSP_Stride. I don't know if the inline decorator is valid or not ; )

/// Wrapper of vDSP conversion function
/// - Parameters:
///   - srcptr: A source pointer
///   - srcStride: A source stride
///   - dstptr: A destination pointer
///   - dstStride: A destination stride
///   - size: A size to be copied
///   - vDSP_func: The vDSP conversion function
@inline(__always)
internal func wrap_vDSP_convert<T, U>(_ size: Int, _ srcptr: UnsafePointer<T>, _ srcStride: Int, _ dstptr: UnsafeMutablePointer<U>, _ dstStride: Int, _ vDSP_func: vDSP_convert_func<T, U>){
    vDSP_func(srcptr, vDSP_Stride(srcStride), dstptr, vDSP_Stride(dstStride), vDSP_Length(size))
}

Original code

• Last, Create function for connecting MfArray's memory and vDSP_vneg. The naming rule is **_by_vDSP.

/// Pre operation mfarray by vDSP
/// - Parameters:
///   - mfarray: An input mfarray
///   - vDSP_func: vDSP_convert_func
/// - Returns: Pre operated mfarray
internal func preop_by_vDSP<T: MfStorable>(_ mfarray: MfArray, _ vDSP_func: vDSP_convert_func<T, T>) -> MfArray{
    //return mfarray must be either row or column major
    var mfarray = mfarray
    //print(mfarray)
    mfarray = check_contiguous(mfarray)
    //print(mfarray)
    //print(mfarray.strides)
    
    let newdata = MfData(size: mfarray.storedSize, mftype: mfarray.mftype)
    newdata.withUnsafeMutableStartPointer(datatype: T.self){
        dstptrT in
        mfarray.withUnsafeMutableStartPointer(datatype: T.self){
            [unowned mfarray] in
            wrap_vDSP_convert(mfarray.storedSize, $0, 1, dstptrT, 1, vDSP_func)
            //vDSP_func($0.baseAddress!, vDSP_Stride(1), dstptrT, vDSP_Stride(1), vDSP_Length(mfarray.storedSize))
        }
    }
    
    let newstructure = MfStructure(shape: mfarray.shape, strides: mfarray.strides)
    return MfArray(mfdata: newdata, mfstructure: newstructure)
}

Original code

3. Finally, use the **_by_vDSP in the Matft public function which the users can use in Sources/Matft/function/**/@@.swift. The public function is in Sources/Matft/function/static or Sources/Matft/function/method. Because Matft stores float or double as StoredType, I change the argument of the **_by_vDSP.

/**
    Element-wise negativity
    - parameters:
        - mfarray: mfarray
*/
public static func neg(_ mfarray: MfArray) -> MfArray{
    switch mfarray.storedType{
        case .Float:
            return preop_by_vDSP(mfarray, vDSP_vneg)
        case .Double:
            return preop_by_vDSP(mfarray, vDSP_vnegD)
    }
}

Original code