Matrix/Vector coordinate list conversions

QuadTree.Tests/Tests.Matrix.fs

@@ -226,3 +226,131 @@ let ``Simple Matrix.map2. Square where number of cols and rows are not power of
     let eq = actual = expected
 
     Assert.True(eq)
+
+[<Fact>]
+let ``Conversion identity`` () =
+    let id = toCoordinateList << fromCoordinateList
+
+    let nrows = 10UL<nrows>
+    let ncols = 12UL<ncols>
+
+    let data =
+        [ 0UL<rowindex>, 3UL<colindex>, 10
+          3UL<rowindex>, 3UL<colindex>, 33
+          9UL<rowindex>, 2UL<colindex>, 5
+          3UL<rowindex>, 11UL<colindex>, 1 ]
+        |> List.sort
+
+    let coordinates = CoordinateList(nrows, ncols, data)
+
+    let expected = coordinates
+    let actual = id coordinates
+
+    Assert.Equal(expected, actual)
+
+[<Fact>]
+let ``Simple addition`` () =
+    let nrows = 10UL<nrows>
+    let ncols = 12UL<ncols>
+
+    let d1 =
+        [ 0UL<rowindex>, 3UL<colindex>, 4
+          9UL<rowindex>, 2UL<colindex>, 5
+          3UL<rowindex>, 11UL<colindex>, 2 ]
+
+    let d2 =
+        [ 0UL<rowindex>, 3UL<colindex>, 6
+          3UL<rowindex>, 3UL<colindex>, 33
+          3UL<rowindex>, 11UL<colindex>, -1 ]
+
+    let expected =
+        let expectedList =
+            [ 0UL<rowindex>, 3UL<colindex>, 10
+              3UL<rowindex>, 3UL<colindex>, 33
+              9UL<rowindex>, 2UL<colindex>, 5
+              3UL<rowindex>, 11UL<colindex>, 1 ]
+            |> List.sort
+
+        CoordinateList(nrows, ncols, expectedList)
+
+    let actual =
+        let c1 = CoordinateList(nrows, ncols, d1)
+        let c2 = CoordinateList(nrows, ncols, d2)
+        let m1 = fromCoordinateList c1
+        let m2 = fromCoordinateList c2
+
+        let addition o1 o2 =
+            match o1, o2 with
+            | Some x, Some y -> Some(x + y)
+            | Some x, None
+            | None, Some x -> Some x
+            | None, None -> None
+
+        let result =
+            match map2 m1 m2 addition with
+            | Result.Success x -> x
+            | _ -> failwith "Unreachable"
+
+        toCoordinateList result
+
+    Assert.Equal(expected, actual)
+
+[<Fact>]
+let ``Condensation of empty`` () =
+    let clist = CoordinateList(2UL<nrows>, 3UL<ncols>, [])
+
+    let actual = fromCoordinateList clist
+
+    // 2 * 3 = 5
+    // 4 * 4 None and Dummy
+    // NN N D
+    // NN N D
+    // DDDD
+    // DDDD
+    let tree =
+        qtree.Node(
+            qtree.Leaf <| UserValue None,
+            qtree.Node(qtree.Leaf <| UserValue None, qtree.Leaf Dummy, qtree.Leaf <| UserValue None, qtree.Leaf Dummy),
+            qtree.Leaf Dummy,
+            qtree.Leaf Dummy
+        )
+
+    let expected =
+        SparseMatrix(2UL<nrows>, 3UL<ncols>, 0UL<nvals>, Storage(4UL<storageVSize>, 4UL<storageHSize>, tree))
+
+    Assert.Equal(expected.storage.data, actual.storage.data)
+
+[<Fact>]
+let ``Condensation of sparse`` () =
+    let clist =
+        CoordinateList(4UL<nrows>, 3UL<ncols>, [ 0UL<rowindex>, 2UL<colindex>, 2; 3UL<rowindex>, 2UL<colindex>, 4 ])
+
+    let actual = fromCoordinateList clist
+
+    // NN2D
+    // NNND
+    // NNND
+    // NN4D
+
+    let tree =
+        qtree.Node(
+            qtree.Leaf <| UserValue None,
+            qtree.Node(
+                qtree.Leaf << UserValue <| Some 2,
+                qtree.Leaf Dummy,
+                qtree.Leaf <| UserValue None,
+                qtree.Leaf Dummy
+            ),
+            qtree.Leaf <| UserValue None,
+            qtree.Node(
+                qtree.Leaf <| UserValue None,
+                qtree.Leaf Dummy,
+                qtree.Leaf << UserValue <| Some 4,
+                qtree.Leaf Dummy
+            )
+        )
+
+    let expected =
+        SparseMatrix(4UL<nrows>, 3UL<ncols>, 0UL<nvals>, Storage(4UL<storageVSize>, 4UL<storageHSize>, tree))
+
+    Assert.Equal(expected.storage.data, actual.storage.data)

QuadTree.Tests/Tests.Vector.fs

@@ -57,7 +57,6 @@ let ``Simple Vector.map2. Length is power of two.`` () =
 
     Assert.True(eq)
 
-
 [<Fact>]
 let ``Simple Vector.map2. Length is not power of two.`` () =
     let v1 =
@@ -100,3 +99,71 @@ let ``Simple Vector.map2. Length is not power of two.`` () =
     let eq = actual = expected
 
     Assert.True(eq)
+
+[<Fact>]
+let ``Conversion identity`` () =
+    let id = toCoordinateList << fromCoordinateList
+
+    let dataLength = 10UL<dataLength>
+
+    let data =
+        [ 0UL<index>, 3; 3UL<index>, -1; 7UL<index>, 2; 8UL<index>, 2; 9UL<index>, 2 ]
+
+    let coordinates = CoordinateList(dataLength, data)
+
+    let expected = coordinates
+    let actual = id coordinates
+
+    Assert.Equal(expected, actual)
+
+[<Fact>]
+let ``Simple addition`` () =
+    let dataLength = 10UL<dataLength>
+
+    let d1 = [ 0UL<index>, 2; 9UL<index>, 1 ]
+    let d2 = [ 0UL<index>, 3; 8UL<index>, 1 ]
+
+    let expected =
+        let expectedList = [ 0UL<index>, 5; 8UL<index>, 1; 9UL<index>, 1 ]
+        CoordinateList(dataLength, expectedList)
+
+    let actual =
+        let c1 = CoordinateList(dataLength, d1)
+        let c2 = CoordinateList(dataLength, d2)
+        let v1 = fromCoordinateList c1
+        let v2 = fromCoordinateList c2
+
+        let addition o1 o2 =
+            match o1, o2 with
+            | Some x, Some y -> Some(x + y)
+            | Some x, None
+            | None, Some x -> Some x
+            | None, None -> None
+
+        let result =
+            match map2 v1 v2 addition with
+            | Result.Success x -> x
+            | _ -> failwith "Unreachable"
+
+        toCoordinateList result
+
+    Assert.Equal(expected, actual)
+
+[<Fact>]
+let ``Condensation of empty`` () =
+    let clist = CoordinateList(10UL<dataLength>, [])
+
+    let actual = fromCoordinateList clist
+
+    // 16 elements total None and Dummy: NNNNNNNN | NN DD | DDDD
+    let tree =
+
+        btree.Node(
+            btree.Leaf <| UserValue None,
+            btree.Node(btree.Node(btree.Leaf <| UserValue None, btree.Leaf Dummy), btree.Leaf Dummy)
+        )
+
+    let expected =
+        SparseVector(clist.length, 0UL<nvals>, Storage(16UL<storageSize>, tree))
+
+    Assert.Equal(expected, actual)

QuadTree/Matrix.fs

@@ -57,6 +57,98 @@ let mkNode x1 x2 x3 x4 =
     | Leaf(v1), Leaf(v2), Leaf(v3), Leaf(v4) when v1 = v2 && v2 = v3 && v3 = v4 -> Leaf(v1)
     | _ -> Node(x1, x2, x3, x4)
 
+[<Measure>]
+type rowindex
+
+[<Measure>]
+type colindex
+
+type COOEntry<'value> = uint64<rowindex> * uint64<colindex> * 'value
+
+[<Struct>]
+type CoordinateList<'value> =
+    val nrows: uint64<nrows>
+    val ncols: uint64<ncols>
+    val list: COOEntry<'value> list
+
+    new(_nrows, _ncols, _list) =
+        { nrows = _nrows
+          ncols = _ncols
+          list = _list }
+
+let private getQuadrantCoords (pr, pc) halfSize =
+    (pr, pc), // NORTH WEST
+    (pr, pc + halfSize * 1UL<colindex>), // NORTH EAST
+    (pr + halfSize * 1UL<rowindex>, pc), // SOUTH WEST
+    (pr + halfSize * 1UL<rowindex>, pc + halfSize * 1UL<colindex>) // SOUTH EAST
+
+let fromCoordinateList (coo: CoordinateList<'a>) =
+    let nvals = (uint64 <| List.length coo.list) * 1UL<nvals>
+    let nrows = coo.nrows
+    let ncols = coo.ncols
+
+    // the resulting matrix is always square
+    let storageSize = getNearestUpperPowerOfTwo (max (uint64 nrows) (uint64 ncols))
+
+    let isEntryInQuadrant (pr, pc) size (entry: COOEntry<'a>) =
+        let (i, j, _) = entry
+
+        i >= pr
+        && j >= pc
+        && i < pr + size * 1UL<rowindex>
+        && j < pc + size * 1UL<colindex>
+
+    let rec traverse coordinates (pr, pc) size =
+        match coordinates with
+        | [] when (uint64 pr) + size < uint64 nrows && (uint64 pc) + size < uint64 ncols -> Leaf <| UserValue None
+        | [] when uint64 pr >= uint64 nrows || uint64 pc >= uint64 ncols -> Leaf Dummy
+        | (i, j, value) :: _ when pr = i && pc = j && size = 1UL -> Leaf << UserValue <| Some value
+        | _ ->
+            let halfSize = size / 2UL
+            let nwp, nep, swp, sep = getQuadrantCoords (pr, pc) halfSize
+            let nwCoo = coordinates |> List.filter (isEntryInQuadrant nwp halfSize)
+            let neCoo = coordinates |> List.filter (isEntryInQuadrant nep halfSize)
+            let swCoo = coordinates |> List.filter (isEntryInQuadrant swp halfSize)
+            let seCoo = coordinates |> List.filter (isEntryInQuadrant sep halfSize)
+
+            mkNode
+                (traverse nwCoo nwp halfSize)
+                (traverse neCoo nep halfSize)
+                (traverse swCoo swp halfSize)
+                (traverse seCoo sep halfSize)
+
+    let tree = traverse coo.list (0UL<rowindex>, 0UL<colindex>) storageSize
+
+    SparseMatrix(nrows, ncols, nvals, Storage(storageSize * 1UL<storageVSize>, storageSize * 1UL<storageHSize>, tree))
+
+let toCoordinateList (matrix: SparseMatrix<'a>) =
+    let nrows = matrix.nrows
+    let ncols = matrix.ncols
+
+    let rec traverse tree (pr, pc) size =
+        match tree with
+        | Leaf Dummy
+        | Leaf(UserValue None) -> []
+        | Leaf(UserValue(Some value)) ->
+            [ for i in uint64 pr .. (uint64 pr) + size - 1UL do
+                  for j in uint64 pc .. (uint64 pc) + size - 1UL -> (i * 1UL<rowindex>, j * 1UL<colindex>, value) ]
+        | Node(nw, ne, sw, se) ->
+            let halfSize = size / 2UL
+            let nwp, nep, swp, sep = getQuadrantCoords (pr, pc) halfSize
+
+            traverse nw nwp halfSize
+            @ traverse ne nep halfSize
+            @ traverse sw swp halfSize
+            @ traverse se sep halfSize
+
+    let coo =
+        traverse
+            matrix.storage.data
+            (0UL<rowindex>, 0UL<colindex>)
+            (max (uint64 matrix.storage.hSize) (uint64 matrix.storage.vSize))
+
+    CoordinateList(nrows, ncols, coo)
+
 let map2 (matrix1: SparseMatrix<_>) (matrix2: SparseMatrix<_>) f =
     let rec inner (vSize: uint64<storageVSize>) (hSize: uint64<storageHSize>) matrix1 matrix2 =
         let _do x1 x2 x3 x4 y1 y2 y3 y4 =

QuadTree/Vector.fs

@@ -38,16 +38,73 @@
     | Leaf
 *)
 
+
 let mkNode t1 t2 =
     match (t1, t2) with
     | Leaf(v1), Leaf(v2) when v1 = v2 -> Leaf(v1)
     | _ -> Node(t1, t2)
 
+[<Measure>]
+type index
+
+[<Struct>]
+type CoordinateList<'value> =
+    val length: uint64<dataLength>
+    val data: (uint64<index> * 'value) list
+    new(_length, _data) = { length = _length; data = _data }
+
+let fromCoordinateList (lst: CoordinateList<'a>) : SparseVector<'a> =
+    let length = lst.length
+    let nvals = (uint64 <| List.length lst.data) * 1UL<nvals>
+    let storageSize = (getNearestUpperPowerOfTwo <| uint64 length) * 1UL<storageSize>
+
+    let rec traverse coordinates pointer size =
+        match coordinates with
+        | [] when uint64 (pointer + size) < uint64 (length) -> Leaf <| UserValue None, []
+        | [] when uint64 pointer >= uint64 length -> Leaf Dummy, []
+        | (idx, _) :: _ when idx > pointer + size -> Leaf <| UserValue None, coordinates
+        | (idx, value) :: xs when idx = pointer && size = 1UL<index> -> Leaf << UserValue <| Some value, xs
+        | _ ->
+            let halfSize = size / 2UL
+
+            let left, lCoordinates = traverse coordinates pointer halfSize
+            let right, rCoordinates = traverse lCoordinates (pointer + halfSize) halfSize
+
+            mkNode left right, rCoordinates
+
+    let sortedCoordinates = List.sort lst.data
+
+    let tree, _ =
+        traverse sortedCoordinates 0UL<index> ((uint64 storageSize) * 1UL<index>)
+
+    SparseVector(length, nvals, Storage(storageSize, tree))
+
+let toCoordinateList (vector: SparseVector<'a>) =
+    let length = vector.length
+
+    let rec traverse tree accum (pointer: uint64<index>) (size: uint64<index>) =
+        match tree with
+        | Leaf Dummy
+        | Leaf(UserValue(None)) -> accum
+        | Leaf(UserValue(Some value)) ->
+            accum
+            @ [ for idx in 0UL .. uint64 (size - 1UL<index>) -> (pointer + idx * 1UL<index>, value) ]
+        | Node(left, right) ->
+            let halfSize = size / 2UL
+            let lAccum = traverse left accum pointer halfSize
+            let rAccum = traverse right lAccum (pointer + halfSize) halfSize
+            rAccum
+
+    let lst =
+        traverse vector.storage.data [] 0UL<index> ((uint64 vector.storage.size) * 1UL<index>)
+
+    CoordinateList(length, lst)
+
 let map2 (vector1: SparseVector<'a>) (vector2: SparseVector<'b>) f =
     let len1 = vector1.length
 
     let rec inner (size: uint64<storageSize>) vector1 vector2 =
         match (vector1, vector2) with
         | Node(t1, t2), Leaf(_) ->
             let new_t1, nvals1 = inner (size / 2UL) t1 vector2
             let new_t2, nvals2 = inner (size / 2UL) t2 vector2