将 F# 度量单位与 System.Numerics.Vector<T> 一起使用答案

【问题标题】：Using F# Units of Measure with System.Numerics.Vector<T>将 F# 度量单位与 System.Numerics.Vector<T> 一起使用
【发布时间】：2016-04-06 17:43:55
【问题描述】：

我很难将 F# 度量单位与 System.Numerics.Vector<'T> 类型结合使用。让我们看一个玩具问题：假设我们有一个float<m>[] 类型的数组xs，并且出于某种原因我们想要对其所有组件进行平方，从而得到一个float<m^2>[] 类型的数组。这与标量代码完美配合：

xs |> Array.map (fun x -> x * x) // float<m^2>[]

现在假设我们想通过使用 SIMD 在System.Numerics.Vector<float>.Count 大小的块中执行乘法来向量化这个操作，例如：

open System.Numerics
let simdWidth = Vector<float>.Count
// fill with dummy data
let xs = Array.init (simdWidth * 10) (fun i -> float i * 1.0<m>)
// array to store the results
let rs: float<m^2> array = Array.zeroCreate (xs |> Array.length)
// number of SIMD operations required
let chunks = (xs |> Array.length) / simdWidth
// for simplicity, assume xs.Length % simdWidth = 0
for i = 0 to chunks - 1 do
    let v = Vector<_>(xs, i * simdWidth) // Vector<float<m>>, containing xs.[i .. i+simdWidth-1]
    let u = v * v                        // Vector<float<m>>; expected: Vector<float<m^2>>
    u.CopyTo(rs, i * simdWidth)          // units mismatch

我相信我理解为什么会发生这种情况：F# 编译器怎么知道System.Numerics.Vector<'T>.op_Multiply 做了什么以及应用了哪些算术规则？它实际上可以是任何操作。那么它应该如何推导出正确的单位呢？

问题是：完成这项工作的最佳方法是什么？我们如何告诉编译器适用哪些规则？

尝试 1：从 xs 中删除所有度量单位信息，稍后再添加：

// remove all UoM from all arrays
let xsWoM = Array.map (fun x -> x / 1.0<m>) xs
// ...
// perform computation using xsWoM etc.
// ...
// add back units again
let xs = Array.map (fun x -> x * 1.0<m>) xsWoM

问题：执行不必要的计算和/或复制操作，出于性能原因无法实现矢量化代码的目的。此外，这在很大程度上违背了使用计量单位开始的目的。

尝试2：使用内联IL改变Vector<'T>.op_Multiply的返回类型：

// reinterpret x to be of type 'b
let inline retype (x: 'a) : 'b = (# "" x: 'b #)
let inline (.*.) (u: Vector<float<'m>>) (v: Vector<float<'m>>): Vector<float<'m^2>> = u * v |> retype
// ...
let u = v .*. v // asserts type Vector<float<m^2>>

问题：不需要任何额外的操作，但使用了已弃用的功能（内联 IL）并且不是完全通用的（仅与测量单位有关）。

有没有人对此有更好的解决方案*？

*请注意，上面的示例确实是一个玩具问题，用于演示一般问题。实际程序解决了一个更复杂的初值问题，涉及多种物理量。

【问题讨论】：

标签： .net f# vectorization simd units-of-measurement

【解决方案1】：

编译器可以很好地确定如何将单元规则应用于乘法，这里的问题是你有一个包装类型。在您的第一个示例中，当您编写 xs |> Array.map (fun x -> x * x) 时，您是在根据数组元素而不是直接在数组上描述乘法。

当您有一个Vector<float<m>> 时，这些单位附加到float 而不是Vector 所以当您尝试乘以Vectors 时，编译器不会将该类型视为它有任何单位。

鉴于该类公开的方法，我认为直接使用Vector<'T> 不是一个简单的解决方法，但有包装类型的选项。

这样的东西可以给你一个单元友好的向量：

type VectorWithUnits<'a, [<Measure>]'b> = 
    |VectorWithUnits of Vector<'a>

    static member inline (*) (a : VectorWithUnits<'a0, 'b0>, b : VectorWithUnits<'a0, 'b1>) 
        : VectorWithUnits<'a0, 'b0*'b1> =
        match a, b with
        |VectorWithUnits a, VectorWithUnits b -> VectorWithUnits <| a * b

在这种情况下，单位被附加到向量上，并且乘法向量在正确的单位行为下按预期工作。

问题是现在我们可以在Vector<'T> 和float 本身上拥有单独且不同的度量单位注释。

您可以使用以下方法将具有度量单位的特定类型的数组转换为一组Vectors：

let toFloatVectors (array : float<'m>[]) : VectorWithUnits<float,'m>[]  =
    let arrs = array |> Array.chunkBySize (Vector<float>.Count)
    arrs |> Array.map (Array.map (float) >> Vector >> VectorWithUnits)

然后返回：

let fromFloatVectors (vectors : VectorWithUnits<float,'m>[]) : float<'m>[] =
    let arr = Array.zeroCreate<float> (Array.length vectors)
    vectors |> Array.iteri (fun i uVec ->
        match uVec with
        |VectorWithUnits vec -> vec.CopyTo arr)
    arr |> Array.map (LanguagePrimitives.FloatWithMeasure<'m>)

一个老生常谈的替代方案：

如果您放弃泛型类型'T，您可以通过一些相当可怕的装箱和运行时强制转换使float Vector 正常运行。这滥用了度量单位是在运行时不再存在的编译时构造这一事实。

type FloatVectorWithUnits<[<Measure>]'b> = 
    |FloatVectorWithUnits of Vector<float<'b>>

    static member ( * ) (a : FloatVectorWithUnits<'b0>, b : FloatVectorWithUnits<'b1>) =
        match a, b with
        |FloatVectorWithUnits a, FloatVectorWithUnits b ->
            let c, d = box a :?> Vector<float<'b0*'b1>>, box b :?> Vector<float<'b0*'b1>>
            c * d |> FloatVectorWithUnits

【讨论】：

嗯，我尝试了FloatVectorWithUnits，但(*) 运算符的度量会自动限制为“1”（无量纲，签名FloatVectorWithUnits<1> * FloatVectorWithUnits<1> -> FloatVectorWithUnits<1>）。此外，VectorWithUnits 需要额外的显式类型约束，但似乎可以正常工作。我也考虑过包装类型，但我不知道你可以断言/强制测量参数的某些规则。很有趣。
@Frank 好奇。我从System.Numerics.Vectors 安装Vector<'T> 类型时遇到了一些问题-我读到Vector<'T> 类型实际上已从稳定版本中删除，因此我不得不以一种稍微迂回的方式测试此代码，但我没有得到相同的结果约束为“1”。我想知道这是否可能与不同的 F# 版本有关，将尝试调查。
@Frank 很抱歉，我认为第二个例子是不行的。它适用于我对那些安装问题的稍微狡猾的解决方法，但我认为它不会对你有用，除非你能找到一种方法来剥离度量单位注释并在我知道的算术函数定义中明确地重新添加它们你想避免。希望第一种方法仍然有用...
是的，NuGet 包仍然存在问题，至少对于 .NET framework/F# 和 Numerics.Vector 版本的某些组合（另请参阅github.com/dotnet/corefx/issues/3103）。它适用于 .NET 4.5.2、F# 4.0 和 Vectors 4.1.1-beta-23516。
@Frank 所以，有一种方法可以使用一些运行时强制转换来重新创建第二个示例，我猜这类似于您原来的内联 IL 方法。

【解决方案2】：

我想出了一个解决方案，可以满足我的大部分要求（看起来）。它的灵感来自TheInnerLight's ideas（包装Vector<'T>），但也为底层数组数据类型添加了一个包装器（称为ScalarField）。这样我们可以跟踪单元，而在下面我们只处理原始数据并且可以使用不感知单元的System.Numerics.Vector API。

一个简单的、简单的快速而肮脏的实现应该是这样的：

// units-aware wrapper for System.Numerics.Vector<'T>
type PackedScalars<[<Measure>] 'm> = struct
    val public Data: Vector<float>
    new (d: Vector<float>) = {Data = d}
    static member inline (*) (u: PackedScalars<'m1>, v: PackedScalars<'m2>) = u.Data * v.Data |> PackedScalars<'m1*'m2>
end

// unit-ware type, wrapping a raw array for easy stream processing
type ScalarField<[<Measure>] 'm> = struct
    val public Data: float[]
    member self.Item with inline get i                = LanguagePrimitives.FloatWithMeasure<'m> self.Data.[i]
                     and  inline set i (v: float<'m>) = self.Data.[i] <- (float v)
    member self.Packed 
           with inline get i                        = Vector<float>(self.Data, i) |> PackedScalars<'m>
           and  inline set i (v: PackedScalars<'m>) = v.Data.CopyTo(self.Data, i)
    new (d: float[]) = {Data = d}
    new (count: int) = {Data = Array.zeroCreate count}
end

我们现在可以使用这两种数据结构以一种相对优雅、高效的方式解决示例问题：

let xs = Array.init (simdWidth * 10) float |> ScalarField<m>    
let mutable rs = Array.zeroCreate (xs.Data |> Array.length) |> ScalarField<m^2>
let chunks = (xs.Data |> Array.length) / simdWidth
for i = 0 to chunks - 1 do
    let j = i * simdWidth
    let v = xs.Packed(j) // PackedScalars<m>
    let u = v * v        // PackedScalars<m^2>
    rs.Packed(j) <- u

最重要的是，为单元感知 ScalarField 包装器重新实现通常的数组操作可能很有用，例如

[<CompilationRepresentation(CompilationRepresentationFlags.ModuleSuffix)>]
module ScalarField =
    let map f (sf: ScalarField<_>) =
        let mutable res = Array.zeroCreate sf.Data.Length |> ScalarField
        for i = 0 to sf.Data.Length do
           res.[i] <- f sf.[i]
        res

等等

缺点：相对于底层数字类型 (float) 不是通用的，因为没有通用替代 floatWithMeasure。为了使其通用，我们必须实现第三个包装器Scalar，它还包装了底层原语：

type Scalar<'a, [<Measure>] 'm> = struct
    val public Data: 'a
    new (d: 'a) = {Data = d}
end

type PackedScalars<'a, [<Measure>] 'm 
            when 'a: (new: unit -> 'a) 
            and  'a: struct 
            and  'a :> System.ValueType> = struct
    val public Data: Vector<'a>
    new (d: Vector<'a>) = {Data = d}
    static member inline (*) (u: PackedScalars<'a, 'm1>, v: PackedScalars<'a, 'm2>) = u.Data * v.Data |> PackedScalars<'a, 'm1*'m2>
end

type ScalarField<'a, [<Measure>] 'm
            when 'a: (new: unit -> 'a) 
            and  'a: struct 
            and  'a :> System.ValueType> = struct
    val public Data: 'a[]
    member self.Item with inline get i                    = Scalar<'a, 'm>(self.Data.[i])
                     and  inline set i (v: Scalar<'a,'m>) = self.Data.[i] <- v.Data
    member self.Packed 
           with inline get i                          = Vector<'a>(self.Data, i) |> PackedScalars<_,'m>
           and  inline set i (v: PackedScalars<_,'m>) = v.Data.CopyTo(self.Data, i)
    new (d:'a[]) = {Data = d}
    new (count: int) = {Data = Array.zeroCreate count}
end

...这意味着我们基本上不是通过使用像float<'m> 这样的“细化”类型来跟踪单位，而是仅通过带有辅助类型/单位参数的包装器类型。

不过，我仍然希望有人提出更好的主意。 :)

【讨论】：