diff --git a/QuadTree.Benchmark/Exists.fs b/QuadTree.Benchmark/Exists.fs new file mode 100644 index 0000000..7339ea1 --- /dev/null +++ b/QuadTree.Benchmark/Exists.fs @@ -0,0 +1,81 @@ +namespace QuadTree.Benchmarks.Exists + + +open BenchmarkDotNet.Attributes +open QuadTree.Benchmarks.Utils + +[)>] +type Benchmark() = + + let mutable denseVec = Unchecked.defaultof> + let mutable sparseVec = Unchecked.defaultof> + let mutable denseMat = Unchecked.defaultof> + + let cooVec (lst: Vector.CoordinateList<'a>) : Vector.SparseVector<'a> = + match Vector.fromCoordinateList lst with + | Ok v -> v + | Error e -> failwith e + + let cooMat (lst: Matrix.CoordinateList<'a>) : Matrix.SparseMatrix<'a> = + match Matrix.fromCoordinateList lst with + | Ok m -> m + | Error e -> failwith e + + let mutable sparseMat = Unchecked.defaultof> + + [] + member val N = 0 with get, set + + [] + member this.Setup() = + let denseData = + [ for i in 0UL .. uint64 this.N - 1UL -> (i * 1UL, int i % 100) ] + + denseVec <- cooVec (Vector.CoordinateList(uint64 this.N * 1UL, denseData)) + + let sparseData = + [ for i in 0UL .. 10UL .. uint64 this.N - 1UL -> (i * 1UL, int i % 100) ] + + sparseVec <- cooVec (Vector.CoordinateList(uint64 this.N * 1UL, sparseData)) + + let denseMatData = + [ for r in 0UL .. uint64 this.N - 1UL do + for c in 0UL .. uint64 this.N - 1UL do + (r * 1UL, c * 1UL, int (r + c) % 100) ] + + denseMat <- + cooMat ( + Matrix.CoordinateList( + uint64 this.N * 1UL, + uint64 this.N * 1UL, + denseMatData + ) + ) + + let sparseMatData = + [ for r in 0UL .. 10UL .. uint64 this.N - 1UL do + for c in 0UL .. 10UL .. uint64 this.N - 1UL do + (r * 1UL, c * 1UL, int (r + c) % 100) ] + + sparseMat <- + cooMat ( + Matrix.CoordinateList( + uint64 this.N * 1UL, + uint64 this.N * 1UL, + sparseMatData + ) + ) + + [] + member this.VectorExistsDense() = Vector.exists denseVec (fun x -> x < 0) + + [] + member this.VectorExistsSparse() = + Vector.exists sparseVec (fun x -> x < 0) + + [] + member this.MatrixExistsDense() = Matrix.exists denseMat (fun x -> x < 0) + + [] + member this.MatrixExistsSparse() = + Matrix.exists sparseMat (fun x -> x < 0) diff --git a/QuadTree.Benchmark/Filters.fs b/QuadTree.Benchmark/Filters.fs new file mode 100644 index 0000000..f177796 --- /dev/null +++ b/QuadTree.Benchmark/Filters.fs @@ -0,0 +1,82 @@ +namespace QuadTree.Benchmarks.Filters + +open BenchmarkDotNet.Attributes +open QuadTree.Benchmarks.Utils + +[)>] +type Benchmark() = + + let mutable denseVec = Unchecked.defaultof> + let mutable sparseVec = Unchecked.defaultof> + let mutable denseMat = Unchecked.defaultof> + + let cooVec (lst: Vector.CoordinateList<'a>) : Vector.SparseVector<'a> = + match Vector.fromCoordinateList lst with + | Ok v -> v + | Error e -> failwith e + + let cooMat (lst: Matrix.CoordinateList<'a>) : Matrix.SparseMatrix<'a> = + match Matrix.fromCoordinateList lst with + | Ok m -> m + | Error e -> failwith e + + let mutable sparseMat = Unchecked.defaultof> + + [] + member val N = 0 with get, set + + [] + member this.Setup() = + let denseData = + [ for i in 0UL .. uint64 this.N - 1UL -> (i * 1UL, int i % 100) ] + + denseVec <- cooVec (Vector.CoordinateList(uint64 this.N * 1UL, denseData)) + + let sparseData = + [ for i in 0UL .. 10UL .. uint64 this.N - 1UL -> (i * 1UL, int i % 100) ] + + sparseVec <- cooVec (Vector.CoordinateList(uint64 this.N * 1UL, sparseData)) + + let denseMatData = + [ for r in 0UL .. uint64 this.N - 1UL do + for c in 0UL .. uint64 this.N - 1UL do + (r * 1UL, c * 1UL, int (r + c) % 100) ] + + denseMat <- + cooMat ( + Matrix.CoordinateList( + uint64 this.N * 1UL, + uint64 this.N * 1UL, + denseMatData + ) + ) + + let sparseMatData = + [ for r in 0UL .. 10UL .. uint64 this.N - 1UL do + for c in 0UL .. 10UL .. uint64 this.N - 1UL do + (r * 1UL, c * 1UL, int (r + c) % 100) ] + + sparseMat <- + cooMat ( + Matrix.CoordinateList( + uint64 this.N * 1UL, + uint64 this.N * 1UL, + sparseMatData + ) + ) + + [] + member this.VectorFilterDense() = + Vector.filter denseVec (fun x -> x % 2 = 0) + + [] + member this.VectorFilterSparse() = + Vector.filter sparseVec (fun x -> x % 2 = 0) + + [] + member this.MatrixFilterDense() = + Matrix.filter denseMat (fun x -> x % 2 = 0) + + [] + member this.MatrixFilterSparse() = + Matrix.filter sparseMat (fun x -> x % 2 = 0) diff --git a/QuadTree.Benchmark/Forall.fs b/QuadTree.Benchmark/Forall.fs new file mode 100644 index 0000000..de7f671 --- /dev/null +++ b/QuadTree.Benchmark/Forall.fs @@ -0,0 +1,83 @@ +namespace QuadTree.Benchmarks.Forall + + +open BenchmarkDotNet.Attributes +open QuadTree.Benchmarks.Utils + +[)>] +type Benchmark() = + + let mutable denseVec = Unchecked.defaultof> + let mutable sparseVec = Unchecked.defaultof> + let mutable denseMat = Unchecked.defaultof> + + let cooVec (lst: Vector.CoordinateList<'a>) : Vector.SparseVector<'a> = + match Vector.fromCoordinateList lst with + | Ok v -> v + | Error e -> failwith e + + let cooMat (lst: Matrix.CoordinateList<'a>) : Matrix.SparseMatrix<'a> = + match Matrix.fromCoordinateList lst with + | Ok m -> m + | Error e -> failwith e + + let mutable sparseMat = Unchecked.defaultof> + + [] + member val N = 0 with get, set + + [] + member this.Setup() = + let denseData = + [ for i in 0UL .. uint64 this.N - 1UL -> (i * 1UL, int i % 100) ] + + denseVec <- cooVec (Vector.CoordinateList(uint64 this.N * 1UL, denseData)) + + let sparseData = + [ for i in 0UL .. 10UL .. uint64 this.N - 1UL -> (i * 1UL, int i % 100) ] + + sparseVec <- cooVec (Vector.CoordinateList(uint64 this.N * 1UL, sparseData)) + + let denseMatData = + [ for r in 0UL .. uint64 this.N - 1UL do + for c in 0UL .. uint64 this.N - 1UL do + (r * 1UL, c * 1UL, int (r + c) % 100) ] + + denseMat <- + cooMat ( + Matrix.CoordinateList( + uint64 this.N * 1UL, + uint64 this.N * 1UL, + denseMatData + ) + ) + + let sparseMatData = + [ for r in 0UL .. 10UL .. uint64 this.N - 1UL do + for c in 0UL .. 10UL .. uint64 this.N - 1UL do + (r * 1UL, c * 1UL, int (r + c) % 100) ] + + sparseMat <- + cooMat ( + Matrix.CoordinateList( + uint64 this.N * 1UL, + uint64 this.N * 1UL, + sparseMatData + ) + ) + + [] + member this.VectorForallDense() = + Vector.forall denseVec (fun x -> x >= 0) + + [] + member this.VectorForallSparse() = + Vector.forall sparseVec (fun x -> x >= 0) + + [] + member this.MatrixForallDense() = + Matrix.forall denseMat (fun x -> x >= 0) + + [] + member this.MatrixForallSparse() = + Matrix.forall sparseMat (fun x -> x >= 0) diff --git a/QuadTree.Benchmark/FormatBenchmarks.fs b/QuadTree.Benchmark/FormatBenchmarks.fs new file mode 100644 index 0000000..d9e154e --- /dev/null +++ b/QuadTree.Benchmark/FormatBenchmarks.fs @@ -0,0 +1,497 @@ +namespace QuadTree.Benchmarks.Formats + +open System +open BenchmarkDotNet.Attributes +open Matrix +open COOArray + +[)>] +type FormatBenchmark() = + + let mutable cooMatrix1 = Unchecked.defaultof> + let mutable cooMatrix2 = Unchecked.defaultof> + let mutable qtMatrix1 = Unchecked.defaultof> + let mutable qtMatrix2 = Unchecked.defaultof> + let mutable listMatrix1 = Unchecked.defaultof> + let mutable listMatrix2 = Unchecked.defaultof> + + let mutable lookupCoords: (uint64 * uint64) array = [||] + let mutable lookupValues: double array = [||] + + let mutable resultCoo = Unchecked.defaultof> + let mutable resultQt = Unchecked.defaultof> + let mutable resultList = Unchecked.defaultof> + let mutable resultCooVal = 0.0 + let mutable resultQtVal = 0.0 + let mutable resultListVal = 0.0 + + [] + member val Size = 0 with get, set + + [] + member val FillRate = 0.0 with get, set + + [] + member this.Setup() = + let rng = Random(42) + let size = uint64 this.Size + let totalCells = float (size * size) + let targetNnz = max 10 (int (totalCells * this.FillRate)) + + let generateEntries count = + let entries = System.Collections.Generic.HashSet() + + [ 1..count ] + |> List.map (fun _ -> + let mutable i = 0UL + let mutable j = 0UL + + while entries.Contains((i, j)) || i >= size || j >= size do + i <- uint64 (rng.Next(int size)) + j <- uint64 (rng.Next(int size)) + + entries.Add((i, j)) |> ignore + (i * 1UL, j * 1UL, rng.NextDouble() * 100.0)) + |> List.sort + + let entries1 = generateEntries targetNnz + let entries2 = generateEntries targetNnz + + cooMatrix1 <- CoordinateList(size * 1UL, size * 1UL, entries1) + cooMatrix2 <- CoordinateList(size * 1UL, size * 1UL, entries2) + qtMatrix1 <- fromCoordinateList cooMatrix1 + qtMatrix2 <- fromCoordinateList cooMatrix2 + listMatrix1 <- COOList.fromArray cooMatrix1 + listMatrix2 <- COOList.fromArray cooMatrix2 + + lookupCoords <- entries1 |> List.map (fun (i, j, _) -> (i, j)) |> Array.ofList + lookupValues <- entries1 |> List.map (fun (_, _, v) -> v) |> Array.ofList + + [] + member this.CooMap() = + resultCoo <- cooMap cooMatrix1 (fun v -> v |> Option.map (fun x -> x * 2.0)) + + [] + member this.QtMap() = + resultQt <- map qtMatrix1 (fun v -> v |> Option.map (fun x -> x * 2.0)) + + [] + member this.CooMapi() = + resultCoo <- + cooMapi cooMatrix1 (fun i j v -> v |> Option.map (fun x -> x + float (uint64 i) + float (uint64 j))) + + [] + member this.QtMapi() = + resultQt <- mapi qtMatrix1 (fun i j v -> v |> Option.map (fun x -> x + float (uint64 i) + float (uint64 j))) + + [] + member this.CooMap2() = + match + cooMap2 cooMatrix1 cooMatrix2 (fun a b -> + match a, b with + | Some x, Some y -> Some(x + y) + | Some x, None -> Some x + | None, Some y -> Some y + | None, None -> None) + with + | Ok r -> resultCoo <- r + | Error _ -> () + + [] + member this.QtMap2() = + match + map2 qtMatrix1 qtMatrix2 (fun a b -> + match a, b with + | Some x, Some y -> Some(x + y) + | Some x, None -> Some x + | None, Some y -> Some y + | None, None -> None) + with + | Ok r -> resultQt <- r + | Error _ -> () + + [] + member this.CooMap2i() = + match + cooMap2i cooMatrix1 cooMatrix2 (fun i j a b -> + match a, b with + | Some x, Some y -> Some(x + y + float (uint64 i)) + | Some x, None -> Some x + | None, Some y -> Some y + | None, None -> None) + with + | Ok r -> resultCoo <- r + | Error _ -> () + + [] + member this.QtMap2i() = + match + map2i qtMatrix1 qtMatrix2 (fun i j a b -> + match a, b with + | Some x, Some y -> Some(x + y + float (uint64 i)) + | Some x, None -> Some x + | None, Some y -> Some y + | None, None -> None) + with + | Ok r -> resultQt <- r + | Error _ -> () + + [] + member this.CooListMap() = + resultList <- COOList.cooMap listMatrix1 (fun v -> v |> Option.map (fun x -> x * 2.0)) + + [] + member this.CooListMapi() = + resultList <- + COOList.cooMapi listMatrix1 (fun i j v -> + v |> Option.map (fun x -> x + float (uint64 i) + float (uint64 j))) + + [] + member this.CooListMap2() = + match + COOList.cooMap2 listMatrix1 listMatrix2 (fun a b -> + match a, b with + | Some x, Some y -> Some(x + y) + | Some x, None -> Some x + | None, Some y -> Some y + | None, None -> None) + with + | Ok r -> resultList <- r + | Error _ -> () + + [] + member this.CooListMap2i() = + match + COOList.cooMap2i listMatrix1 listMatrix2 (fun i j a b -> + match a, b with + | Some x, Some y -> Some(x + y + float (uint64 i)) + | Some x, None -> Some x + | None, Some y -> Some y + | None, None -> None) + with + | Ok r -> resultList <- r + | Error _ -> () + + [] + member this.CooListMxm() = + let op_add x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | None, None -> None + + let op_mult x y = + match x, y with + | Some a, Some b -> Some(a * b) + | _ -> None + + match COOList.mxmcoo op_add op_mult listMatrix1 listMatrix1 with + | Ok result -> resultList <- result + | Error _ -> failwith "COOList mxmcoo failed" + + [] + member this.CooListGet() = + let n = min lookupCoords.Length 1000 + let mutable acc = 0.0 + + for k = 0 to n - 1 do + let (i, j) = lookupCoords.[k] + + match COOList.cooGet (listMatrix1, i, j) with + | Ok(Some v) -> acc <- acc + v + | _ -> () + + resultListVal <- acc + + [] + member this.CooListSet() = + let mutable m = listMatrix1 + let n = min lookupCoords.Length 1000 + + for k = 0 to n - 1 do + let (i, j) = lookupCoords.[k] + + match COOList.cooUpdate (m, i, j, lookupValues.[k] * 2.0) with + | Ok updated -> m <- updated + | _ -> () + + resultList <- m + + [] + member this.CooGet() = + let n = min lookupCoords.Length 1000 + let mutable acc = 0.0 + + for k = 0 to n - 1 do + let (i, j) = lookupCoords.[k] + + match cooGet (cooMatrix1, i, j) with + | Ok(Some v) -> acc <- acc + v + | _ -> () + + resultCooVal <- acc + + [] + member this.QtGet() = + let n = min lookupCoords.Length 1000 + let mutable acc = 0.0 + + for k = 0 to n - 1 do + let (i, j) = lookupCoords.[k] + + match get qtMatrix1 i j with + | Ok(Some v) -> acc <- acc + v + | _ -> () + + resultQtVal <- acc + + [] + member this.CooSet() = + let n = min lookupCoords.Length 1000 + let mutable m = cooMatrix1 + + for k = 0 to n - 1 do + let (i, j) = lookupCoords.[k] + + match cooUpdate (m, i, j, lookupValues.[k] * 2.0) with + | Ok updated -> m <- updated + | _ -> () + + resultCoo <- m + + [] + member this.QtSet() = + let n = min lookupCoords.Length 1000 + let mutable m = qtMatrix1 + + for k = 0 to n - 1 do + let (i, j) = lookupCoords.[k] + + match set m i j (lookupValues.[k] * 2.0) with + | Ok updated -> m <- updated + | _ -> () + + resultQt <- m + + [] + member this.CooMxm() = + let op_add x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | None, None -> None + + let op_mult x y = + match x, y with + | Some a, Some b -> Some(a * b) + | _ -> None + + match mxmcoo op_add op_mult cooMatrix1 cooMatrix1 with + | Ok result -> resultCoo <- result + | Error _ -> failwith "mxmcoo failed" + + [] + member this.QtMxm() = + let op_add x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | None, None -> None + + let op_mult x y = + match x, y with + | Some a, Some b -> Some(a * b) + | _ -> None + + match LinearAlgebra.mxm op_add op_mult qtMatrix1 qtMatrix1 with + | Ok result -> resultQt <- result + | Error _ -> failwith "mxm failed" + + +[)>] +type DenseFormatBenchmark() = + + let mutable cooMatrix = Unchecked.defaultof> + let mutable qtMatrix = Unchecked.defaultof> + let mutable listMatrix = Unchecked.defaultof> + + let mutable resultCoo = Unchecked.defaultof> + let mutable resultQt = Unchecked.defaultof> + let mutable resultList = Unchecked.defaultof> + let mutable resultCooVal = 0.0 + let mutable resultQtVal = 0.0 + let mutable resultListVal = 0.0 + + [] + member val Size = 0 with get, set + + [] + member this.Setup() = + let rng = Random(42) + let size = uint64 this.Size + + let entries = + [ for i in 0UL .. size - 1UL do + for j in 0UL .. size - 1UL do + (i * 1UL, j * 1UL, rng.NextDouble() * 100.0) ] + + cooMatrix <- CoordinateList(size * 1UL, size * 1UL, entries) + qtMatrix <- fromCoordinateList cooMatrix + listMatrix <- COOList.fromArray cooMatrix + + [] + member this.DenseCooMap() = + resultCoo <- cooMap cooMatrix (fun v -> v |> Option.map (fun x -> x * 2.0)) + + [] + member this.DenseQtMap() = + resultQt <- map qtMatrix (fun v -> v |> Option.map (fun x -> x * 2.0)) + + [] + member this.DenseCooMapi() = + resultCoo <- cooMapi cooMatrix (fun i j v -> v |> Option.map (fun x -> x + float (uint64 i) + float (uint64 j))) + + [] + member this.DenseQtMapi() = + resultQt <- mapi qtMatrix (fun i j v -> v |> Option.map (fun x -> x + float (uint64 i) + float (uint64 j))) + + [] + member this.DenseCooListMap() = + resultList <- COOList.cooMap listMatrix (fun v -> v |> Option.map (fun x -> x * 2.0)) + + [] + member this.DenseCooListMapi() = + resultList <- + COOList.cooMapi listMatrix (fun i j v -> v |> Option.map (fun x -> x + float (uint64 i) + float (uint64 j))) + + [] + member this.DenseCooListMxm() = + let op_add x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | None, None -> None + + let op_mult x y = + match x, y with + | Some a, Some b -> Some(a * b) + | _ -> None + + match COOList.mxmcoo op_add op_mult listMatrix listMatrix with + | Ok result -> resultList <- result + | Error _ -> failwith "COOList mxmcoo failed" + + [] + member this.DenseCooListGet() = + let mutable acc = 0.0 + + for i in 0UL .. uint64 this.Size - 1UL do + for j in 0UL .. uint64 this.Size - 1UL do + match COOList.cooGet (listMatrix, i * 1UL, j * 1UL) with + | Ok(Some v) -> acc <- acc + v + | _ -> () + + resultListVal <- acc + + [] + member this.DenseCooListSet() = + let mutable m = listMatrix + let size = uint64 this.Size + + for i in 0UL .. size - 1UL do + for j in 0UL .. size - 1UL do + match COOList.cooUpdate (m, i * 1UL, j * 1UL, 42.0) with + | Ok updated -> m <- updated + | _ -> () + + resultList <- m + + [] + member this.DenseCooGet() = + let mutable acc = 0.0 + + for i in 0UL .. uint64 this.Size - 1UL do + for j in 0UL .. uint64 this.Size - 1UL do + match cooGet (cooMatrix, i * 1UL, j * 1UL) with + | Ok(Some v) -> acc <- acc + v + | _ -> () + + resultCoo <- cooMatrix + + [] + member this.DenseQtGet() = + let mutable acc = 0.0 + + for i in 0UL .. uint64 this.Size - 1UL do + for j in 0UL .. uint64 this.Size - 1UL do + match get qtMatrix (i * 1UL) (j * 1UL) with + | Ok(Some v) -> acc <- acc + v + | _ -> () + + resultQt <- qtMatrix + + [] + member this.DenseCooSet() = + let mutable m = cooMatrix + let size = uint64 this.Size + + for i in 0UL .. size - 1UL do + for j in 0UL .. size - 1UL do + match cooUpdate (m, i * 1UL, j * 1UL, 42.0) with + | Ok updated -> m <- updated + | _ -> () + + resultCoo <- m + + [] + member this.DenseQtSet() = + let mutable m = qtMatrix + let size = uint64 this.Size + + for i in 0UL .. size - 1UL do + for j in 0UL .. size - 1UL do + match set m (i * 1UL) (j * 1UL) 42.0 with + | Ok updated -> m <- updated + | _ -> () + + resultQt <- m + + [] + member this.DenseCooMxm() = + let op_add x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | None, None -> None + + let op_mult x y = + match x, y with + | Some a, Some b -> Some(a * b) + | _ -> None + + match mxmcoo op_add op_mult cooMatrix cooMatrix with + | Ok result -> resultCoo <- result + | Error _ -> failwith "mxmcoo failed" + + [] + member this.DenseQtMxm() = + let op_add x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | None, None -> None + + let op_mult x y = + match x, y with + | Some a, Some b -> Some(a * b) + | _ -> None + + match LinearAlgebra.mxm op_add op_mult qtMatrix qtMatrix with + | Ok result -> resultQt <- result + | Error _ -> failwith "mxm failed" diff --git a/QuadTree.Benchmark/Main.fs b/QuadTree.Benchmark/Main.fs index 25cfb19..161ab81 100644 --- a/QuadTree.Benchmark/Main.fs +++ b/QuadTree.Benchmark/Main.fs @@ -4,7 +4,10 @@ open BenchmarkDotNet.Running let main argv = let benchmarks = BenchmarkSwitcher - [| typeof + [| typeof + typeof + typeof + typeof typeof typeof typeof diff --git a/QuadTree.Benchmark/QuadTree.Benchmark.fsproj b/QuadTree.Benchmark/QuadTree.Benchmark.fsproj index 9fbd475..3e5fe49 100644 --- a/QuadTree.Benchmark/QuadTree.Benchmark.fsproj +++ b/QuadTree.Benchmark/QuadTree.Benchmark.fsproj @@ -1,4 +1,4 @@ - + Exe @@ -8,6 +8,9 @@ + + + @@ -29,3 +32,4 @@ + diff --git a/QuadTree.Benchmark/RealMatrixBenchmark.fs b/QuadTree.Benchmark/RealMatrixBenchmark.fs new file mode 100644 index 0000000..d09bb38 --- /dev/null +++ b/QuadTree.Benchmark/RealMatrixBenchmark.fs @@ -0,0 +1,256 @@ +namespace QuadTree.Benchmarks.RealMatrices + +open System +open System.IO +open BenchmarkDotNet.Attributes +open Matrix +open COOArray + +[)>] +[] +type RealMatrixBenchmark() = + + let mutable cooMatrix = Unchecked.defaultof> + let mutable qtMatrix = Unchecked.defaultof> + let mutable listMatrix = Unchecked.defaultof> + + let mutable resultCoo = Unchecked.defaultof> + let mutable resultQt = Unchecked.defaultof> + let mutable resultList = Unchecked.defaultof> + let mutable resultCooVal = 0.0 + let mutable resultQtVal = 0.0 + let mutable resultListVal = 0.0 + + let mutable lookupCoords: (uint64 * uint64) array = [||] + let mutable lookupValues: double array = [||] + + let mutable matrixName = "" + let mutable doMxm = false + let mutable skip = false + + [] + member val MatrixName = "" with get, set + + [] + member this.Setup() = + matrixName <- this.MatrixName + skip <- false + let dataDir = Path.GetFullPath(QuadTree.Benchmarks.Utils.DIR_WITH_MATRICES) + let mtxPath = Path.Combine(dataDir, matrixName + ".mtx") + + if not (File.Exists mtxPath) then + skip <- true + else + let isSymmetric = + File.ReadLines(mtxPath) + |> Seq.exists (fun s -> s.StartsWith "%%MatrixMarket" && s.Contains "symmetric") + + let (coo, qt) = QuadTree.Benchmarks.Utils.readMtxRaw mtxPath (not isSymmetric) + + cooMatrix <- coo + qtMatrix <- qt + listMatrix <- COOList.fromArray coo + + let nnz = coo.list.Length + let dim = max (uint64 coo.nrows) (uint64 coo.ncols) + doMxm <- nnz < 100000 && dim <= 12119UL + + let rng = Random(42) + let sampleSize = min nnz 1000 + let coords = coo.list + let indices = Array.init sampleSize (fun _ -> rng.Next(nnz)) + + lookupCoords <- + indices + |> Array.map (fun k -> + let (i, j, _) = coords.[k] + (i, j)) + + lookupValues <- + indices + |> Array.map (fun k -> + let (_, _, v) = coords.[k] + v) + + [] + member this.CooMap() = + if not skip then + resultCoo <- cooMap cooMatrix (fun v -> v |> Option.map (fun x -> x * 2.0)) + + [] + member this.QtMap() = + if not skip then + resultQt <- map qtMatrix (fun v -> v |> Option.map (fun x -> x * 2.0)) + + [] + member this.CooMapi() = + if not skip then + resultCoo <- + cooMapi cooMatrix (fun i j v -> v |> Option.map (fun x -> x + float (uint64 i) + float (uint64 j))) + + [] + member this.QtMapi() = + if not skip then + resultQt <- mapi qtMatrix (fun i j v -> v |> Option.map (fun x -> x + float (uint64 i) + float (uint64 j))) + + [] + member this.CooListMap() = + if not skip then + resultList <- COOList.cooMap listMatrix (fun v -> v |> Option.map (fun x -> x * 2.0)) + + [] + member this.CooListMapi() = + if not skip then + resultList <- + COOList.cooMapi listMatrix (fun i j v -> + v |> Option.map (fun x -> x + float (uint64 i) + float (uint64 j))) + + [] + member this.CooListMxm() = + if not skip && doMxm then + let op_add x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | None, None -> None + + let op_mult x y = + match x, y with + | Some a, Some b -> Some(a * b) + | _ -> None + + match COOList.mxmcoo op_add op_mult listMatrix listMatrix with + | Ok result -> resultList <- result + | Error _ -> failwith "COOList mxmcoo failed" + + [] + member this.CooListGet() = + if not skip then + let mutable acc = 0.0 + + for k = 0 to lookupCoords.Length - 1 do + let (i, j) = lookupCoords.[k] + + match COOList.cooGet (listMatrix, i, j) with + | Ok(Some v) -> acc <- acc + v + | _ -> () + + resultListVal <- acc + + [] + member this.CooListSet() = + if not skip then + let mutable m = listMatrix + + for k = 0 to lookupCoords.Length - 1 do + let (i, j) = lookupCoords.[k] + + match COOList.cooUpdate (m, i, j, lookupValues.[k] * 2.0) with + | Ok updated -> m <- updated + | _ -> () + + resultList <- m + + [] + member this.CooGet() = + if not skip then + let mutable acc = 0.0 + + for k = 0 to lookupCoords.Length - 1 do + let (i, j) = lookupCoords.[k] + + match cooGet (cooMatrix, i, j) with + | Ok(Some v) -> acc <- acc + v + | _ -> () + + resultCooVal <- acc + + [] + member this.QtGet() = + if not skip then + let mutable acc = 0.0 + + for k = 0 to lookupCoords.Length - 1 do + let (i, j) = lookupCoords.[k] + + match get qtMatrix i j with + | Ok(Some v) -> acc <- acc + v + | _ -> () + + resultQtVal <- acc + + [] + member this.CooSet() = + if not skip then + let mutable m = cooMatrix + + for k = 0 to lookupCoords.Length - 1 do + let (i, j) = lookupCoords.[k] + + match cooUpdate (m, i, j, lookupValues.[k] * 2.0) with + | Ok updated -> m <- updated + | _ -> () + + resultCoo <- m + + [] + member this.QtSet() = + if not skip then + let mutable m = qtMatrix + + for k = 0 to lookupCoords.Length - 1 do + let (i, j) = lookupCoords.[k] + + match set m i j (lookupValues.[k] * 2.0) with + | Ok updated -> m <- updated + | _ -> () + + resultQt <- m + + [] + member this.CooMxm() = + if not skip && doMxm then + let op_add x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | None, None -> None + + let op_mult x y = + match x, y with + | Some a, Some b -> Some(a * b) + | _ -> None + + match mxmcoo op_add op_mult cooMatrix cooMatrix with + | Ok result -> resultCoo <- result + | Error _ -> failwith "mxmcoo failed" + + [] + member this.QtMxm() = + if not skip && doMxm then + let op_add x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | None, None -> None + + let op_mult x y = + match x, y with + | Some a, Some b -> Some(a * b) + | _ -> None + + match LinearAlgebra.mxm op_add op_mult qtMatrix qtMatrix with + | Ok result -> resultQt <- result + | Error _ -> failwith "mxm failed" diff --git a/QuadTree.Benchmark/Utils.fs b/QuadTree.Benchmark/Utils.fs index 5ec8793..fbb98be 100644 --- a/QuadTree.Benchmark/Utils.fs +++ b/QuadTree.Benchmark/Utils.fs @@ -8,7 +8,7 @@ type MyConfig() = let DIR_WITH_MATRICES = "../../../../../../../data/" -let readMtx path directed = +let readMtxRaw path directed = let getCooList (linewords: seq) = linewords |> Seq.map (fun x -> @@ -23,7 +23,7 @@ let readMtx path directed = let lines = File.ReadLines(path) let removedComments = lines |> Seq.skipWhile (fun s -> s.[0] = '%') - let linewords = removedComments |> Seq.map (fun s -> s.Split " ") + let linewords = removedComments |> Seq.map (fun s -> s.Split [| ' ' |]) let first = Seq.head linewords let nrows, ncols, nnz = uint64 first.[0], uint64 first.[1], int first.[2] @@ -33,7 +33,16 @@ let readMtx path directed = let lst = getCooList tl if (directed && nnz <> lst.Length) || ((not directed) && nnz * 2 <> lst.Length) then - failwithf "Incorrect matrix reading. Path: %A expected nnz: %A actual nnz: %A" path (nnz * 2) lst.Length + failwithf + "Incorrect matrix reading. Path: %A expected nnz: %A actual nnz: %A" + path + (if directed then nnz else nnz * 2) + lst.Length - Matrix.CoordinateList(nrows * 1UL, nrows * 1UL, lst) - |> Matrix.fromCoordinateList + let coo = + Matrix.CoordinateList(nrows * 1UL, ncols * 1UL, lst) + + let qt = Matrix.fromCoordinateList coo + (coo, qt) + +let readMtx path directed = readMtxRaw path directed |> snd diff --git a/QuadTree.Tests/PropertyTests.fs b/QuadTree.Tests/PropertyTests.fs new file mode 100644 index 0000000..7a91ea7 --- /dev/null +++ b/QuadTree.Tests/PropertyTests.fs @@ -0,0 +1,254 @@ +module QuadTree.Tests.PropertyTests + +open System +open Xunit +open FsCheck +open FsCheck.FSharp +open FsCheck.Xunit +open Matrix +open COOArray + +type Input = + { Rows: int + Cols: int + Cells: (int * int * int) list } + +let private fromCoordinateListUnchecked (lst: CoordinateList<'a>) = + match Matrix.fromCoordinateList lst with + | Ok m -> m + | Error e -> failwith e + +let private toCoo (inp: Input) : CoordinateList = + let nrows = max 1 inp.Rows + let ncols = max 1 inp.Cols + + let entries = + inp.Cells + |> List.map (fun (r, c, v) -> (abs r, abs c, v)) + |> List.filter (fun (r, c, _) -> r < nrows && c < ncols) + |> List.distinctBy (fun (r, c, _) -> (r, c)) + |> List.map (fun (r, c, v) -> (uint64 r * 1UL, uint64 c * 1UL, v)) + |> List.sortBy (fun (r, c, _) -> (r, c)) + + CoordinateList(uint64 nrows * 1UL, uint64 ncols * 1UL, entries) + +let private arbInput: Arbitrary = + let gen = + gen { + let! rows = Gen.choose (1, 16) + let! cols = Gen.choose (1, 16) + + let! cells = + Gen.listOf ( + gen { + let! r = Gen.choose (-5, 20) + let! c = Gen.choose (-5, 20) + let! v = Gen.choose (-100, 100) + return (r, c, v) + } + ) + + return + { Rows = rows + Cols = cols + Cells = cells } + } + + Arb.fromGen gen + +type InputArbs = + static member Input() = arbInput + +[ |])>] +let ``get at every cell agrees between QuadTree and COOArray`` (inp: Input) = + let coo = toCoo inp + let cooA = ArrayCOO(coo.nrows, coo.ncols, coo.list) + let qt = fromCoordinateListUnchecked coo + let nrows = int (uint64 coo.nrows) + let ncols = int (uint64 coo.ncols) + + List.allPairs [ 0 .. nrows - 1 ] [ 0 .. ncols - 1 ] + |> List.forall (fun (r, c) -> + let ri = uint64 r * 1UL + let ci = uint64 c * 1UL + + Matrix.get qt ri ci = cooGet (cooA, ri, ci)) + +[ |])>] +let ``toCoordinateList (fromCoordinateListUnchecked coo) preserves every value`` (inp: Input) = + let coo = toCoo inp + let back = toCoordinateList (fromCoordinateListUnchecked coo) + let backA = ArrayCOO(back.nrows, back.ncols, back.list) + + back.nrows = coo.nrows + && back.ncols = coo.ncols + && List.length back.list = List.length coo.list + && coo.list |> List.forall (fun (r, c, v) -> cooGet (backA, r, c) = Ok(Some v)) + +[ |])>] +let ``cooUpdate writes a value and adjusts the length`` (inp: Input) = + let coo = toCoo inp + let cooA = ArrayCOO(coo.nrows, coo.ncols, coo.list) + let nrows = int (uint64 coo.nrows) + let ncols = int (uint64 coo.ncols) + let r = abs inp.Rows % nrows + let c = abs inp.Cols % ncols + let ri = uint64 r * 1UL + let ci = uint64 c * 1UL + let wasPresent = cooA.list |> Array.exists (fun (i, j, _) -> i = ri && j = ci) + + match cooUpdate (cooA, ri, ci, 777) with + | Ok updated -> + cooGet (updated, ri, ci) = Ok(Some 777) + && Array.length updated.list = Array.length cooA.list + (if wasPresent then 0 else 1) + | Error _ -> false + +[ |])>] +let ``set and cooUpdate agree on the written cell`` (inp: Input) = + let coo = toCoo inp + let cooA = ArrayCOO(coo.nrows, coo.ncols, coo.list) + let qt = fromCoordinateListUnchecked coo + let nrows = int (uint64 coo.nrows) + let ncols = int (uint64 coo.ncols) + let r = abs inp.Rows % nrows + let c = abs inp.Cols % ncols + let ri = uint64 r * 1UL + let ci = uint64 c * 1UL + + match cooUpdate (cooA, ri, ci, 42), Matrix.set qt ri ci 42 with + | Ok updatedCoo, Ok updatedQt -> cooGet (updatedCoo, ri, ci) = Matrix.get updatedQt ri ci + | _ -> false + +[ |])>] +let ``cooMapValues maps every stored value once`` (inp: Input) = + let coo = toCoo inp + let cooA = ArrayCOO(coo.nrows, coo.ncols, coo.list) + let mapped = cooMapValues cooA (fun v -> Some(v + 1)) + + Array.length mapped.list = Array.length cooA.list + && cooA.list + |> Array.forall (fun (r, c, v) -> cooGet (mapped, r, c) = Ok(Some(v + 1))) + +[ |])>] +let ``out-of-bounds access raises ArgumentOutOfRangeException`` (inp: Input) = + let coo = toCoo inp + let cooA = ArrayCOO(coo.nrows, coo.ncols, coo.list) + let qt = fromCoordinateListUnchecked coo + let nrows = uint64 coo.nrows * 1UL + let ncols = uint64 coo.ncols * 1UL + + let cooGetThrows = + try + cooGet (cooA, nrows, 0UL) |> ignore + false + with :? ArgumentOutOfRangeException -> + true + + let cooUpdateThrows = + try + cooUpdate (cooA, nrows, 0UL, 1) |> ignore + false + with :? ArgumentOutOfRangeException -> + true + + let matrixGetThrows = + try + Matrix.get qt nrows 0UL |> ignore + false + with :? ArgumentOutOfRangeException -> + true + + let matrixSetThrows = + try + Matrix.set qt nrows 0UL 1 |> ignore + false + with :? ArgumentOutOfRangeException -> + true + + cooGetThrows && cooUpdateThrows && matrixGetThrows && matrixSetThrows + +let private cooWithDims (nrows: uint64) (ncols: uint64) (inp: Input) : CoordinateList = + let entries = + inp.Cells + |> List.map (fun (r, c, v) -> (abs r, abs c, v)) + |> List.filter (fun (r, c, _) -> r < int nrows && c < int ncols) + |> List.distinctBy (fun (r, c, _) -> (r, c)) + |> List.map (fun (r, c, v) -> (uint64 r * 1UL, uint64 c * 1UL, v)) + |> List.sortBy (fun (r, c, _) -> (r, c)) + + CoordinateList(nrows * 1UL, ncols * 1UL, entries) + +let private opAdd x y = + match (x, y) with + | Some(a), Some(b) -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | _ -> None + +let private opMul x y = + match (x, y) with + | Some(a), Some(b) -> Some(a * b) + | _ -> None + +let private naiveMxm (nrowsA: uint64) (k: uint64) (ncolsB: uint64) (m1: COOEntry list) (m2: COOEntry list) = + let m1Map = m1 |> List.map (fun (i, j, v) -> ((i, j), v)) |> Map.ofList + let m2Map = m2 |> List.map (fun (i, j, v) -> ((i, j), v)) |> Map.ofList + + [ for i in 0UL .. nrowsA - 1UL do + for j in 0UL .. ncolsB - 1UL do + let products = + [ for t in 0UL .. k - 1UL do + let a = m1Map |> Map.tryFind (i * 1UL, t * 1UL) + let b = m2Map |> Map.tryFind (t * 1UL, j * 1UL) + yield opMul a b ] + + match products |> List.fold (fun acc p -> opAdd acc p) None with + | Some v -> yield (i * 1UL, j * 1UL, v) + | None -> () ] + |> List.sortBy (fun (i, j, _) -> (i, j)) + +let private arbInputPair: Arbitrary = + let gen = + gen { + let! a = Arb.toGen arbInput + let! b = Arb.toGen arbInput + return (a, b) + } + + Arb.fromGen gen + +type InputPairArbs = + static member InputPair() = arbInputPair + +[ |])>] +let ``mxmcoo agrees with naive multiplication, array matches list, keys are sorted and unique`` + ((a, b): Input * Input) + = + let nrowsA = uint64 (max 1 a.Rows) + let k = uint64 (max 1 a.Cols) + let ncolsB = uint64 (max 1 b.Cols) + let m1 = cooWithDims nrowsA k a + let m2 = cooWithDims k ncolsB b + let m1A = ArrayCOO(m1.nrows, m1.ncols, m1.list) + let m2A = ArrayCOO(m2.nrows, m2.ncols, m2.list) + + let expected = naiveMxm nrowsA k ncolsB m1.list m2.list + + let sortedUnique entries = + entries + |> List.map (fun (i, j, _) -> (i, j)) + |> List.pairwise + |> List.forall (fun ((i1, j1), (i2, j2)) -> i1 < i2 || (i1 = i2 && j1 < j2)) + + match + COOArray.mxmcoo opAdd opMul m1A m2A, COOList.mxmcoo opAdd opMul (COOList.fromArray m1A) (COOList.fromArray m2A) + with + | Ok arr, Ok lst -> + let arrEntries = Array.toList arr.list + + List.indexed expected = List.indexed arrEntries + && List.indexed expected = List.indexed lst.entries + && (arrEntries = lst.entries) + && sortedUnique arrEntries + && sortedUnique lst.entries + | _ -> false diff --git a/QuadTree.Tests/QuadTree.Tests.fsproj b/QuadTree.Tests/QuadTree.Tests.fsproj index bc76cf2..d1e483c 100644 --- a/QuadTree.Tests/QuadTree.Tests.fsproj +++ b/QuadTree.Tests/QuadTree.Tests.fsproj @@ -8,7 +8,9 @@ + + @@ -18,6 +20,7 @@ + diff --git a/QuadTree.Tests/Tests.COOArray.fs b/QuadTree.Tests/Tests.COOArray.fs new file mode 100644 index 0000000..1ca2457 --- /dev/null +++ b/QuadTree.Tests/Tests.COOArray.fs @@ -0,0 +1,1006 @@ +module COOArray.Tests + +open System +open Xunit + +open Matrix +open COOArray +open Common + +let op_add x y = + match (x, y) with + | Some(a), Some(b) -> Some(a + b) + | Some a, None + | None, Some a -> Some a + | _ -> None + +let op_mult x y = + match (x, y) with + | Some(a), Some(b) -> Some(a * b) + | _ -> None + +let private keysAscending (entries: COOEntry<'v> list) = + entries + |> List.map (fun (i, j, _) -> (i, j)) + |> List.pairwise + |> List.forall (fun ((a1, b1), (a2, b2)) -> a1 < a2 || (a1 = a2 && b1 < b2)) + +// === cooGet tests === + +[] +let ``cooGet existing value`` () = + let coo = + ArrayCOO( + 4UL, + 4UL, + [ (0UL, 0UL, 1) + (0UL, 1UL, 2) + (1UL, 0UL, 3) ] + ) + + let actual = cooGet (coo, 0UL, 1UL) + + Assert.Equal(Ok(Some 2), actual) + +[] +let ``cooGet missing value`` () = + let coo = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1); (1UL, 1UL, 2) ]) + + let actual = cooGet (coo, 2UL, 2UL) + + Assert.Equal(Ok None, actual) + +[] +let ``cooGet out of bounds`` () = + let coo = ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1) ]) + + Assert.Throws(fun () -> cooGet (coo, 5UL, 5UL) |> ignore) + +// === cooUpdate tests === + +[] +let ``cooUpdate replaces existing`` () = + let coo = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1); (1UL, 1UL, 2) ]) + + let expected = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 99); (1UL, 1UL, 2) ]) + + let actual = cooUpdate (coo, 0UL, 0UL, 99) + + Assert.Equal(Ok expected, actual) + +[] +let ``cooUpdate inserts new in middle`` () = + let coo = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1); (2UL, 2UL, 2) ]) + + let expected = + ArrayCOO( + 4UL, + 4UL, + [ (0UL, 0UL, 1) + (1UL, 1UL, 10) + (2UL, 2UL, 2) ] + ) + + let actual = cooUpdate (coo, 1UL, 1UL, 10) + + Assert.Equal(Ok expected, actual) + +[] +let ``cooUpdate inserts at end`` () = + let coo = ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1) ]) + + let expected = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1); (3UL, 3UL, 20) ]) + + let actual = cooUpdate (coo, 3UL, 3UL, 20) + + Assert.Equal(Ok expected, actual) + +[] +let ``cooUpdate out of bounds`` () = + let coo = ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1) ]) + + Assert.Throws(fun () -> + cooUpdate (coo, 5UL, 5UL, 99) |> ignore) + +// === cooMap tests === + +[] +let ``cooMap doubles values`` () = + let nrows = 4UL + let ncols = 4UL + + let data = + [ (0UL, 0UL, 1) + (0UL, 1UL, 2) + (1UL, 0UL, 3) + (1UL, 1UL, 4) ] + |> List.sort + + let coo = ArrayCOO(nrows, ncols, data) + + let f v = v |> Option.map (fun v -> v * 2) + + let expected = + ArrayCOO( + nrows, + ncols, + [ (0UL, 0UL, 2) + (0UL, 1UL, 4) + (1UL, 0UL, 6) + (1UL, 1UL, 8) ] + ) + + let actual = cooMap coo f + + Assert.Equal(expected, actual) + +[] +let ``cooMap filters None results`` () = + let nrows = 4UL + let ncols = 4UL + + let data = + [ (0UL, 0UL, 1) + (0UL, 1UL, 2) + (1UL, 0UL, 3) + (1UL, 1UL, 4) ] + + let coo = ArrayCOO(nrows, ncols, data) + + let f v = + v + |> Option.bind (fun v -> + match v with + | 1 -> None + | _ -> Some(v * 10)) + + let expected = + ArrayCOO( + nrows, + ncols, + [ (0UL, 1UL, 20) + (1UL, 0UL, 30) + (1UL, 1UL, 40) ] + ) + + let actual = cooMap coo f + + Assert.Equal(expected, actual) + +[] +let ``cooMap fills missing cells (general form)`` () = + let nrows = 3UL + let ncols = 3UL + + let data = [ (0UL, 0UL, 1); (2UL, 2UL, 5) ] + + let coo = ArrayCOO(nrows, ncols, data) + + let f v = Some(defaultArg v 0) + + let actual = cooMap coo f + + Assert.Equal(nrows, actual.nrows) + Assert.Equal(ncols, actual.ncols) + Assert.Equal(9, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 2UL && j = 2UL) actual.list, + Some(2UL, 2UL, 5) + ) + +[] +let ``cooMap zero-size matrix`` () = + let coo = ArrayCOO(0UL, 0UL, []) + let f v = v |> Option.map (fun v -> v * 2) + let actual = cooMap coo f + let expected = ArrayCOO(0UL, 0UL, []) + Assert.Equal(expected, actual) + +// === cooMap2 tests === + +[] +let ``cooMap2 addition`` () = + let nrows = 10UL + let ncols = 12UL + + let d1 = + [ (0UL, 3UL, 4) + (3UL, 11UL, 2) + (9UL, 2UL, 5) ] + |> List.sort + + let d2 = + [ (0UL, 3UL, 6) + (3UL, 3UL, 33) + (3UL, 11UL, -1) ] + |> List.sort + + let f x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None -> Some a + | None, Some b -> Some b + | _ -> None + + let expected = + ArrayCOO( + nrows, + ncols, + [ (0UL, 3UL, 10) + (3UL, 3UL, 33) + (9UL, 2UL, 5) + (3UL, 11UL, 1) ] + |> List.sort + ) + + let c1 = ArrayCOO(nrows, ncols, d1) + let c2 = ArrayCOO(nrows, ncols, d2) + + let actual = cooMap2 c1 c2 f + + Assert.Equal(Ok expected, actual) + +[] +let ``cooMap2 with mismatched positions`` () = + let nrows = 4UL + let ncols = 4UL + + let d1 = [ (0UL, 0UL, 1); (2UL, 2UL, 3) ] + + let d2 = [ (1UL, 1UL, 10); (3UL, 3UL, 30) ] + + let f x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None -> Some(a + 100) + | None, Some b -> Some(b + 200) + | _ -> None + + let expected = + ArrayCOO( + nrows, + ncols, + [ (0UL, 0UL, 101) + (1UL, 1UL, 210) + (2UL, 2UL, 103) + (3UL, 3UL, 230) ] + ) + + let c1 = ArrayCOO(nrows, ncols, d1) + let c2 = ArrayCOO(nrows, ncols, d2) + + let actual = cooMap2 c1 c2 f + + Assert.Equal(Ok expected, actual) + +[] +let ``cooMap2 dense filters None from existing entries`` () = + let nrows = 4UL + let ncols = 4UL + + let d1 = [ (0UL, 0UL, 1); (1UL, 1UL, 2) ] + let d2 = [ (0UL, 0UL, 10); (2UL, 2UL, 30) ] + + let f x y = + match x, y with + | Some a, Some b when a + b > 5 -> None + | Some a, Some b -> Some(a + b) + | Some a, None -> Some a + | None, Some b -> Some b + | None, None -> Some 0 + + let expected = + ArrayCOO( + nrows, + ncols, + [ (0UL, 1UL, 0) + (0UL, 2UL, 0) + (0UL, 3UL, 0) + (1UL, 0UL, 0) + (1UL, 1UL, 2) + (1UL, 2UL, 0) + (1UL, 3UL, 0) + (2UL, 0UL, 0) + (2UL, 1UL, 0) + (2UL, 2UL, 30) + (2UL, 3UL, 0) + (3UL, 0UL, 0) + (3UL, 1UL, 0) + (3UL, 2UL, 0) + (3UL, 3UL, 0) ] + ) + + let c1 = ArrayCOO(nrows, ncols, d1) + let c2 = ArrayCOO(nrows, ncols, d2) + + let actual = cooMap2 c1 c2 f + + Assert.Equal(Ok expected, actual) + +// === cooMapi tests === + +[] +let ``cooMapi position-dependent values`` () = + let nrows = 4UL + let ncols = 4UL + + let data = + [ (0UL, 0UL, 1) + (1UL, 1UL, 2) + (2UL, 3UL, 3) ] + |> List.sort + + let coo = ArrayCOO(nrows, ncols, data) + + let f i j v = + v |> Option.map (fun v -> v + (int (uint64 i))) + + let expected = + ArrayCOO( + nrows, + ncols, + [ (0UL, 0UL, 1) + (1UL, 1UL, 3) + (2UL, 3UL, 5) ] + ) + + let actual = cooMapi coo f + + Assert.Equal(expected, actual) + +[] +let ``cooMapi filters None results`` () = + let nrows = 4UL + let ncols = 4UL + + let data = + [ (0UL, 0UL, 1) + (0UL, 1UL, 5) + (1UL, 0UL, 3) ] + + let coo = ArrayCOO(nrows, ncols, data) + + let f _i _j v = + v |> Option.bind (fun v -> if v > 2 then Some(v * 10) else None) + + let expected = + ArrayCOO(nrows, ncols, [ (0UL, 1UL, 50); (1UL, 0UL, 30) ]) + + let actual = cooMapi coo f + + Assert.Equal(expected, actual) + +[] +let ``cooMapi empty input`` () = + let coo = ArrayCOO(4UL, 4UL, []) + let f _i _j v = v |> Option.map (fun v -> v * 2) + let actual = cooMapi coo f + let expected = ArrayCOO(4UL, 4UL, []) + Assert.Equal(expected, actual) + +[] +let ``cooMapi fills missing cells (general form)`` () = + let nrows = 3UL + let ncols = 3UL + + let data = [ (0UL, 0UL, 1); (2UL, 2UL, 5) ] + let coo = ArrayCOO(nrows, ncols, data) + + let f _i _j v = Some(defaultArg v 0) + + let actual = cooMapi coo f + + Assert.Equal(nrows, actual.nrows) + Assert.Equal(ncols, actual.ncols) + Assert.Equal(9, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 2UL && j = 2UL) actual.list, + Some(2UL, 2UL, 5) + ) + +[] +let ``cooMapi position-dependent fill of missing cells`` () = + let nrows = 2UL + let ncols = 2UL + + let data = [ (0UL, 0UL, 7) ] + let coo = ArrayCOO(nrows, ncols, data) + + let f i j v = + match v with + | Some x -> Some x + | None -> Some(int (uint64 i + uint64 j)) + + let actual = cooMapi coo f + + let expected = + [ (0UL, 0UL, 7) + (0UL, 1UL, 1) + (1UL, 0UL, 1) + (1UL, 1UL, 2) ] + + Assert.Equal[]>(Array.ofList expected, actual.list) + +[] +let ``cooMapi zero-size matrix`` () = + let coo = ArrayCOO(0UL, 0UL, []) + let f _i _j v = v |> Option.map (fun v -> v * 2) + let actual = cooMapi coo f + let expected = ArrayCOO(0UL, 0UL, []) + Assert.Equal(expected, actual) + +// === cooMap2i tests === + +[] +let ``cooMap2i position-dependent addition`` () = + let nrows = 4UL + let ncols = 4UL + + let d1 = [ (0UL, 0UL, 1); (2UL, 2UL, 3) ] + let d2 = [ (0UL, 0UL, 10); (2UL, 2UL, 30) ] + + let f i j x y = + match x, y with + | Some a, Some b -> Some(a + b + (int (uint64 i))) + | Some a, None -> Some a + | None, Some b -> Some b + | _ -> None + + let expected = + ArrayCOO(nrows, ncols, [ (0UL, 0UL, 11); (2UL, 2UL, 35) ]) + + let c1 = ArrayCOO(nrows, ncols, d1) + let c2 = ArrayCOO(nrows, ncols, d2) + let actual = cooMap2i c1 c2 f + + Assert.Equal(Ok expected, actual) + +[] +let ``cooMap2i mismatched positions with index`` () = + let nrows = 4UL + let ncols = 4UL + + let d1 = [ (0UL, 0UL, 1); (2UL, 2UL, 3) ] + let d2 = [ (1UL, 1UL, 10); (3UL, 3UL, 30) ] + + let f i j x y = + match x, y with + | Some a, Some b -> Some(a + b) + | Some a, None -> Some(a + (int (uint64 j))) + | None, Some b -> Some(b + (int (uint64 i))) + | _ -> None + + let expected = + ArrayCOO( + nrows, + ncols, + [ (0UL, 0UL, 1) + (1UL, 1UL, 11) + (2UL, 2UL, 5) + (3UL, 3UL, 33) ] + ) + + let c1 = ArrayCOO(nrows, ncols, d1) + let c2 = ArrayCOO(nrows, ncols, d2) + let actual = cooMap2i c1 c2 f + + Assert.Equal(Ok expected, actual) + +[] +let ``cooMap2i filters None results`` () = + let nrows = 4UL + let ncols = 4UL + + let d1 = [ (0UL, 0UL, 1); (1UL, 1UL, 2) ] + let d2 = [ (0UL, 0UL, 2); (2UL, 2UL, 30) ] + + let f i j x y = + match x, y with + | Some a, Some b when a + b > 5 -> None + | Some a, Some b -> Some(a + b) + | _ -> None + + let expected = ArrayCOO(nrows, ncols, [ (0UL, 0UL, 3) ]) + + let c1 = ArrayCOO(nrows, ncols, d1) + let c2 = ArrayCOO(nrows, ncols, d2) + let actual = cooMap2i c1 c2 f + + Assert.Equal(Ok expected, actual) + +[] +let ``cooMap2i empty inputs`` () = + let c1 = ArrayCOO(4UL, 4UL, []) + let c2 = ArrayCOO(4UL, 4UL, []) + let f _i _j x y = None + let actual = cooMap2i c1 c2 f + let expected = ArrayCOO(4UL, 4UL, []) + Assert.Equal(Ok expected, actual) + +// === mxmcoo tests === + +[] +let ``Sparse mxmcoo`` () = + let m1 = + let d = + [ 0UL, 0UL, 1 + 1UL, 1UL, 2 + 2UL, 2UL, 3 ] + + ArrayCOO(3UL, 3UL, d) + + let m2 = + let d = + [ 0UL, 0UL, 3 + 1UL, 1UL, 2 + 2UL, 2UL, 1 ] + + ArrayCOO(3UL, 3UL, d) + + let expected = + let d = + [ 0UL, 0UL, 3 + 1UL, 1UL, 4 + 2UL, 2UL, 3 ] + + ArrayCOO(3UL, 3UL, d) + + match COOArray.mxmcoo op_add op_mult m1 m2 with + | Ok actual -> + Assert.Equal(expected.nrows, actual.nrows) + Assert.Equal(expected.ncols, actual.ncols) + Assert.Equal[]>(expected.list, actual.list) + | Error e -> failwith (e.ToString()) + +[] +let ``Shrinking mxmcoo`` () = + let m1 = + let d = + [ 0UL, 0UL, 1 + 0UL, 2UL, 2 + 1UL, 1UL, 3 ] + + ArrayCOO(2UL, 3UL, d) + + let m2 = + let d = + [ 0UL, 1UL, 4 + 1UL, 0UL, 5 + 2UL, 0UL, 6 ] + + ArrayCOO(3UL, 2UL, d) + + let expected = + let d = + [ 0UL, 0UL, 12 + 0UL, 1UL, 4 + 1UL, 0UL, 15 ] + + ArrayCOO(2UL, 2UL, d) + + match COOArray.mxmcoo op_add op_mult m1 m2 with + | Ok actual -> + Assert.Equal(expected.nrows, actual.nrows) + Assert.Equal(expected.ncols, actual.ncols) + Assert.Equal[]>(expected.list, actual.list) + | Error e -> failwith (e.ToString()) + + +[] +let ``mxmcoo with non-absorbing op_mult`` () = + let op_add x y = + match (x, y) with + | Some(a), Some(b) -> Some(a + b) + | Some a, _ + | _, Some a -> Some a + | _ -> None + + let op_mult x y = + match (x, y) with + | Some(a), Some(b) -> Some(a * b) + | Some a, _ + | _, Some a -> Some a + | _ -> None + + let m1 = + let d = [ 0UL, 0UL, 1; 0UL, 1UL, 2 ] + + ArrayCOO(1UL, 2UL, d) + + let m2 = + let d = [ 0UL, 0UL, 3 ] + + ArrayCOO(2UL, 1UL, d) + + match COOArray.mxmcoo op_add op_mult m1 m2 with + | Ok actual -> + Assert.Equal(1UL, actual.nrows) + Assert.Equal(1UL, actual.ncols) + Assert.Equal(1, actual.list.Length) + Assert.Equal(Some 5, actual.list |> Array.tryHead |> Option.map (fun (_, _, v) -> v)) + | Error e -> failwith (e.ToString()) + +// === mxmcoo list tests === + +let private listCoo nrows ncols entries = COOList.ListCOO(nrows, ncols, entries) + +[] +let ``Sparse mxmcoo list`` () = + let m1 = + listCoo + 3UL + 3UL + [ (0UL, 0UL, 1) + (1UL, 1UL, 2) + (2UL, 2UL, 3) ] + + let m2 = + listCoo + 3UL + 3UL + [ (0UL, 0UL, 3) + (1UL, 1UL, 2) + (2UL, 2UL, 1) ] + + let expected = + [ (0UL, 0UL, 3) + (1UL, 1UL, 4) + (2UL, 2UL, 3) ] + + match COOList.mxmcoo op_add op_mult m1 m2 with + | Ok actual -> + Assert.True((expected = actual.entries), "Sparse mxmcoo list: entries differ") + Assert.True(keysAscending actual.entries) + | Error e -> failwith (e.ToString()) + +[] +let ``Shrinking mxmcoo list`` () = + let m1 = + listCoo + 2UL + 3UL + [ (0UL, 0UL, 1) + (0UL, 2UL, 2) + (1UL, 1UL, 3) ] + + let m2 = + listCoo + 3UL + 2UL + [ (0UL, 1UL, 4) + (1UL, 0UL, 5) + (2UL, 0UL, 6) ] + + match COOList.mxmcoo op_add op_mult m1 m2 with + | Ok actual -> + Assert.Equal(2UL, actual.nrows) + Assert.Equal(2UL, actual.ncols) + + Assert.True( + [ (0UL, 0UL, 12) + (0UL, 1UL, 4) + (1UL, 0UL, 15) ] = + actual.entries + ) + + Assert.True(keysAscending actual.entries) + | Error e -> failwith (e.ToString()) + +[] +let ``mxmcoo with non-absorbing op_mult list`` () = + let op_add x y = + match (x, y) with + | Some(a), Some(b) -> Some(a + b) + | Some a, _ + | _, Some a -> Some a + | _ -> None + + let op_mult x y = + match (x, y) with + | Some(a), Some(b) -> Some(a * b) + | Some a, _ + | _, Some a -> Some a + | _ -> None + + let m1 = + listCoo 1UL 2UL [ (0UL, 0UL, 1); (0UL, 1UL, 2) ] + + let m2 = listCoo 2UL 1UL [ (0UL, 0UL, 3) ] + + match COOList.mxmcoo op_add op_mult m1 m2 with + | Ok actual -> + Assert.Equal(1UL, actual.nrows) + Assert.Equal(1UL, actual.ncols) + Assert.Equal(1, actual.entries.Length) + Assert.Equal(Some 5, actual.entries |> List.tryHead |> Option.map (fun (_, _, v) -> v)) + Assert.True(keysAscending actual.entries) + | Error e -> failwith (e.ToString()) + +[] +let ``mxmcoo collapses products of one cell (array and list)`` () = + let m1 = + ArrayCOO( + 2UL, + 2UL, + [ (0UL, 0UL, 1) + (0UL, 1UL, 2) + (1UL, 1UL, 3) ] + ) + + let m2 = + ArrayCOO(2UL, 2UL, [ (0UL, 0UL, 4); (1UL, 0UL, 5) ]) + + let expected = + [ (0UL, 0UL, 14); (1UL, 0UL, 15) ] + + match + COOArray.mxmcoo op_add op_mult m1 m2, + COOList.mxmcoo op_add op_mult (COOList.fromArray m1) (COOList.fromArray m2) + with + | Ok arr, Ok lst -> + let arrEntries = Array.toList arr.list + Assert.True((expected = arrEntries), "collapses: array result differs") + Assert.True((expected = lst.entries), "collapses: list result differs") + Assert.True(keysAscending arrEntries) + Assert.True(keysAscending lst.entries) + | _ -> failwith "mxmcoo failed" + +[] +let ``mxmcoo result stays sorted when k has multiple hits`` () = + let m1 = + ArrayCOO( + 3UL, + 3UL, + [ (0UL, 0UL, 1) + (0UL, 1UL, 2) + (0UL, 2UL, 3) + (2UL, 0UL, 7) + (2UL, 2UL, 9) ] + ) + + let m2 = + ArrayCOO( + 3UL, + 3UL, + [ (0UL, 0UL, 4) + (0UL, 1UL, 5) + (1UL, 0UL, 6) + (1UL, 2UL, 7) + (2UL, 0UL, 8) ] + ) + + match + COOArray.mxmcoo op_add op_mult m1 m2, + COOList.mxmcoo op_add op_mult (COOList.fromArray m1) (COOList.fromArray m2) + with + | Ok arr, Ok lst -> + let arrEntries = Array.toList arr.list + Assert.True((arrEntries = lst.entries), "array and list results differ") + Assert.True(keysAscending arrEntries) + Assert.True(keysAscending lst.entries) + | _ -> failwith "mxmcoo failed" + +// === cooMapValues / cooMapiValues tests === + +[] +let ``cooMapValues applies only to stored values`` () = + let coo = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1); (1UL, 1UL, 2) ]) + + let actual = cooMapValues coo (fun v -> Some(v * 10)) + + Assert.Equal(2, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 0UL && j = 0UL) actual.list, + Some(0UL, 0UL, 10) + ) + +[] +let ``cooMapiValues applies indexed only to stored values`` () = + let coo = ArrayCOO(4UL, 4UL, [ (1UL, 2UL, 5) ]) + + let actual = + cooMapiValues coo (fun i j v -> Some(v + int (uint64 i) + int (uint64 j))) + + Assert.Equal(1, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 1UL && j = 2UL) actual.list, + Some(1UL, 2UL, 8) + ) + +// === cooMap2 variants tests === + +[] +let ``cooMap2Values applies only where both present`` () = + let c1 = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1); (1UL, 1UL, 2) ]) + + let c2 = ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 10) ]) + + match cooMap2Values c1 c2 (fun a b -> Some(a + b)) with + | Ok actual -> + Assert.Equal(1, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 0UL && j = 0UL) actual.list, + Some(0UL, 0UL, 11) + ) + | Error e -> failwithf "unexpected error %A" e + +[] +let ``cooMap2AllCells equals cooMap2`` () = + let c1 = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1); (1UL, 1UL, 2) ]) + + let c2 = ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 10) ]) + + let f a b = + match a, b with + | Some x, Some y -> Some(x + y) + | _ -> None + + Assert.Equal(cooMap2 c1 c2 f, cooMap2AllCells c1 c2 f) + +[] +let ``cooMap2AtLeastOne distinguishes both left right`` () = + let c1 = + ArrayCOO(3UL, 3UL, [ (0UL, 0UL, 1); (1UL, 1UL, 2) ]) + + let c2 = + ArrayCOO(3UL, 3UL, [ (0UL, 0UL, 10); (2UL, 2UL, 30) ]) + + let f = + function + | AtLeastOne.Both(a, b) -> Some(a + b) + | AtLeastOne.Left a -> Some(a * 100) + | AtLeastOne.Right b -> Some(b * -1) + + match cooMap2AtLeastOne c1 c2 f with + | Ok actual -> + Assert.Equal(3, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 0UL && j = 0UL) actual.list, + Some(0UL, 0UL, 11) + ) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 1UL && j = 1UL) actual.list, + Some(1UL, 1UL, 200) + ) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 2UL && j = 2UL) actual.list, + Some(2UL, 2UL, -30) + ) + | Error e -> failwithf "unexpected error %A" e + +[] +let ``cooMap2LeftValues applies where left present`` () = + let c1 = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1); (1UL, 1UL, 2) ]) + + let c2 = ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 10) ]) + + match cooMap2LeftValues c1 c2 (fun a b -> Some(a + (defaultArg b 0))) with + | Ok actual -> + Assert.Equal(2, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 0UL && j = 0UL) actual.list, + Some(0UL, 0UL, 11) + ) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 1UL && j = 1UL) actual.list, + Some(1UL, 1UL, 2) + ) + | Error e -> failwithf "unexpected error %A" e + +[] +let ``cooMap2 sizes mismatch`` () = + let c1 = ArrayCOO(4UL, 4UL, []) + let c2 = ArrayCOO(2UL, 2UL, []) + let f a b = None + + Assert.Equal(Error Error.InconsistentSizeOfArguments, cooMap2 c1 c2 f) + +// === cooMap2i variants tests === + +[] +let ``cooMap2iValues applies indexed where both present`` () = + let c1 = ArrayCOO(4UL, 4UL, [ (1UL, 1UL, 2) ]) + + let c2 = + ArrayCOO(4UL, 4UL, [ (1UL, 1UL, 10); (2UL, 2UL, 20) ]) + + let f i j a b = + Some(a + b + int (uint64 i) + int (uint64 j)) + + match cooMap2iValues c1 c2 f with + | Ok actual -> + Assert.Equal(1, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 1UL && j = 1UL) actual.list, + Some(1UL, 1UL, 14) + ) + | Error e -> failwithf "unexpected error %A" e + +[] +let ``cooMap2iAllCells equals cooMap2i`` () = + let c1 = + ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 1); (1UL, 1UL, 2) ]) + + let c2 = ArrayCOO(4UL, 4UL, [ (0UL, 0UL, 10) ]) + + let f i j a b = + match a, b with + | Some x, Some y -> Some(x + y + int (uint64 i)) + | _ -> None + + Assert.Equal(cooMap2i c1 c2 f, cooMap2iAllCells c1 c2 f) + +[] +let ``cooMap2iAtLeastOne passes indices and side`` () = + let c1 = ArrayCOO(2UL, 2UL, [ (0UL, 0UL, 1) ]) + + let c2 = + ArrayCOO(2UL, 2UL, [ (0UL, 0UL, 10); (1UL, 1UL, 20) ]) + + let f i j = + function + | AtLeastOne.Both(a, b) -> Some(a + b + int (uint64 i) + int (uint64 j)) + | AtLeastOne.Left a -> Some(a) + | AtLeastOne.Right b -> Some(b + int (uint64 i) * 100 + int (uint64 j)) + + match cooMap2iAtLeastOne c1 c2 f with + | Ok actual -> + Assert.Equal(2, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 0UL && j = 0UL) actual.list, + Some(0UL, 0UL, 11) + ) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 1UL && j = 1UL) actual.list, + Some(1UL, 1UL, 121) + ) + | Error e -> failwithf "unexpected error %A" e + +[] +let ``cooMap2iLeftValues applies indexed where left present`` () = + let c1 = ArrayCOO(2UL, 2UL, [ (1UL, 1UL, 2) ]) + + let c2 = ArrayCOO(2UL, 2UL, [ (1UL, 1UL, 10) ]) + + let f i j a b = + Some(a + (defaultArg b 0) + int (uint64 i) * 10 + int (uint64 j)) + + match cooMap2iLeftValues c1 c2 f with + | Ok actual -> + Assert.Equal(1, actual.list.Length) + + Assert.Equal( + Array.tryFind (fun (i, j, _) -> i = 1UL && j = 1UL) actual.list, + Some(1UL, 1UL, 23) + ) + | Error e -> failwithf "unexpected error %A" e + +[] +let ``cooMap2i sizes mismatch`` () = + let c1 = ArrayCOO(4UL, 4UL, []) + let c2 = ArrayCOO(2UL, 2UL, []) + let f _i _j a b = None + + Assert.Equal(Error Error.InconsistentSizeOfArguments, cooMap2i c1 c2 f) diff --git a/QuadTree.Tests/Tests.LinearAlgebra.fs b/QuadTree.Tests/Tests.LinearAlgebra.fs index f71278c..b012a20 100644 --- a/QuadTree.Tests/Tests.LinearAlgebra.fs +++ b/QuadTree.Tests/Tests.LinearAlgebra.fs @@ -4,6 +4,7 @@ open System open Xunit open Matrix +open COOArray open Vector open Common open Result diff --git a/QuadTree.Tests/Tests.Matrix.fs b/QuadTree.Tests/Tests.Matrix.fs index 9ff0634..ba92055 100644 --- a/QuadTree.Tests/Tests.Matrix.fs +++ b/QuadTree.Tests/Tests.Matrix.fs @@ -4,8 +4,14 @@ open System open Xunit open Matrix +open COOArray open Common +let private fromCoordinateListUnchecked (lst: Matrix.CoordinateList<'a>) = + match Matrix.fromCoordinateList lst with + | Ok m -> m + | Error e -> failwith e + let printMatrix (matrix: SparseMatrix<_>) = printfn "Matrix:" printfn " Rows: %A" matrix.nrows @@ -616,7 +622,393 @@ let ``Fold sum`` () = Assert.Equal(expected, actual) [] -let ``fromCoordinateList with out-of-range coordinates returns Error`` () = +let ``matrix get existing value`` () = + let m = + fromCoordinateListUnchecked ( + CoordinateList( + 4UL, + 4UL, + [ (0UL, 0UL, 7); (1UL, 2UL, 9) ] + ) + ) + + Assert.Equal(Ok(Some 7), get m 0UL 0UL) + Assert.Equal(Ok(Some 9), get m 1UL 2UL) + +[] +let ``matrix get missing value`` () = + let m = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (0UL, 0UL, 7) ])) + + Assert.Equal(Ok None, get m 1UL 1UL) + +[] +let ``matrix get out of bounds`` () = + let m = fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [])) + Assert.Throws(fun () -> get m 5UL 5UL |> ignore) + +[] +let ``matrix set replaces existing`` () = + let m = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (0UL, 0UL, 7) ])) + + let actual = set m 0UL 0UL 99 |> Result.defaultValue m + + Assert.Equal(Ok(Some 99), get actual 0UL 0UL) + Assert.Equal(Ok(Some 7), get m 0UL 0UL) + +[] +let ``matrix set inserts new`` () = + let m = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (0UL, 0UL, 7) ])) + + let actual = set m 2UL 2UL 42 |> Result.defaultValue m + + Assert.Equal(Ok(Some 42), get actual 2UL 2UL) + +[] +let ``matrix set out of bounds`` () = + let m = fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [])) + Assert.Throws(fun () -> set m 5UL 5UL 99 |> ignore) + +[] +let ``matrix set then get roundtrip`` () = + let m0 = empty 4UL 4UL + let m1 = set m0 0UL 0UL 1 |> Result.defaultValue m0 + let m2 = set m1 1UL 2UL 2 |> Result.defaultValue m1 + let m3 = set m2 3UL 3UL 3 |> Result.defaultValue m2 + + Assert.Equal(Ok(Some 1), get m3 0UL 0UL) + Assert.Equal(Ok(Some 2), get m3 1UL 2UL) + Assert.Equal(Ok(Some 3), get m3 3UL 3UL) + Assert.Equal(Ok(None), get m3 2UL 1UL) + +[] +let ``matrix map doubles values`` () = + let m = + fromCoordinateListUnchecked ( + CoordinateList( + 4UL, + 4UL, + [ (0UL, 0UL, 3); (1UL, 2UL, 5) ] + ) + ) + + let result = map m (Option.map (fun v -> v * 2)) + + Assert.Equal(Ok(Some 6), get result 0UL 0UL) + Assert.Equal(Ok(Some 10), get result 1UL 2UL) + Assert.Equal(Ok(None), get result 0UL 1UL) + +[] +let ``matrix map filters Some to None`` () = + let m = + fromCoordinateListUnchecked ( + CoordinateList( + 4UL, + 4UL, + [ (0UL, 0UL, 3); (1UL, 2UL, 5) ] + ) + ) + + let result = + map m (fun v -> + match v with + | Some x when x > 4 -> Some x + | _ -> None) + + Assert.Equal(Ok(None), get result 0UL 0UL) + Assert.Equal(Ok(Some 5), get result 1UL 2UL) + +[] +let ``matrix map fills None with values`` () = + let m = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (0UL, 0UL, 3) ])) + + let result = + map m (fun v -> + Some( + match v with + | Some x -> x + | None -> 0 + )) + + Assert.Equal(Ok(Some 3), get result 0UL 0UL) + Assert.Equal(Ok(Some 0), get result 1UL 1UL) + Assert.Equal(Ok(Some 0), get result 3UL 3UL) + +[] +let ``matrix map on empty matrix`` () = + let m = empty 4UL 4UL + + let result = map m (Option.map (fun v -> v + 1)) + + Assert.Equal(Ok(None), get result 0UL 0UL) + +[] +let ``matrix map nvals updated`` () = + let m = + fromCoordinateListUnchecked ( + CoordinateList( + 4UL, + 4UL, + [ (0UL, 0UL, 3); (1UL, 2UL, 5) ] + ) + ) + + Assert.Equal(2UL, uint64 m.nvals) + + let result = map m (fun _ -> None) + + Assert.Equal(0UL, uint64 result.nvals) + +[] +let ``matrix mapValues applies only to stored values`` () = + let m = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (0UL, 0UL, 3) ])) + + let result = mapValues m (fun v -> Some(v * 2)) + + Assert.Equal(Ok(Some 6), get result 0UL 0UL) + Assert.Equal(Ok(None), get result 2UL 2UL) + Assert.Equal(Ok(None), get result 0UL 1UL) + Assert.Equal(1UL, uint64 result.nvals) + +[] +let ``matrix mapValues can drop stored values`` () = + let m = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (0UL, 0UL, 3) ])) + + let result = mapValues m (fun _ -> None) + + Assert.Equal(0UL, uint64 result.nvals) + +[] +let ``matrix mapiValues applies only to stored values with indices`` () = + let m = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (1UL, 2UL, 5) ])) + + let result = mapiValues m (fun i j v -> Some(v + int (uint64 i) + int (uint64 j))) + + Assert.Equal(Ok(Some 8), get result 1UL 2UL) + Assert.Equal(Ok(None), get result 0UL 0UL) + Assert.Equal(1UL, uint64 result.nvals) + +[] +let ``matrix mapi expands uniform leaf`` () = + let m = + SparseMatrix( + 4UL, + 4UL, + 4UL, + Storage(4UL, Matrix.qtree.Node(leaf_v 5, leaf_n (), leaf_n (), leaf_n ())) + ) + + let result = + mapi m (fun i j v -> v |> Option.map (fun x -> x + int (uint64 i) + int (uint64 j))) + + Assert.Equal(4UL, uint64 result.nvals) + Assert.Equal(Ok(Some 5), get result 0UL 0UL) + Assert.Equal(Ok(Some 6), get result 0UL 1UL) + Assert.Equal(Ok(Some 6), get result 1UL 0UL) + Assert.Equal(Ok(Some 7), get result 1UL 1UL) + +[] +let ``matrix map2Values applies only where both values present`` () = + let m1 = + fromCoordinateListUnchecked ( + CoordinateList( + 4UL, + 4UL, + [ (0UL, 0UL, 1); (1UL, 1UL, 2) ] + ) + ) + + let m2 = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (0UL, 0UL, 10) ])) + + match map2Values m1 m2 (fun a b -> Some(a + b + 100)) with + | Error e -> failwithf "unexpected error %A" e + | Ok result -> + Assert.Equal(Ok(Some 111), get result 0UL 0UL) + Assert.Equal(Ok(None), get result 1UL 1UL) + Assert.Equal(1UL, uint64 result.nvals) + +[] +let ``matrix map2AllCells equals map2`` () = + let m1 = + fromCoordinateListUnchecked ( + CoordinateList( + 4UL, + 4UL, + [ (0UL, 0UL, 1); (1UL, 1UL, 2) ] + ) + ) + + let m2 = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (0UL, 0UL, 10) ])) + + let f a b = + match a, b with + | Some x, Some y -> Some(x + y) + | _ -> None + + Assert.Equal(map2 m1 m2 f, map2AllCells m1 m2 f) + +[] +let ``matrix map2AtLeastOne distinguishes both left right`` () = + let m1 = + fromCoordinateListUnchecked ( + CoordinateList( + 3UL, + 3UL, + [ (0UL, 0UL, 1); (1UL, 1UL, 2) ] + ) + ) + + let m2 = + fromCoordinateListUnchecked ( + CoordinateList( + 3UL, + 3UL, + [ (0UL, 0UL, 10); (2UL, 2UL, 30) ] + ) + ) + + let f = + function + | AtLeastOne.Both(a, b) -> Some(a + b) + | AtLeastOne.Left a -> Some(a * 100) + | AtLeastOne.Right b -> Some(b * -1) + + match map2AtLeastOne m1 m2 f with + | Error e -> failwithf "unexpected error %A" e + | Ok result -> + Assert.Equal(Ok(Some 11), get result 0UL 0UL) + Assert.Equal(Ok(Some 200), get result 1UL 1UL) + Assert.Equal(Ok(Some -30), get result 2UL 2UL) + +[] +let ``matrix map2LeftValues applies where left value present`` () = + let m1 = + fromCoordinateListUnchecked ( + CoordinateList( + 4UL, + 4UL, + [ (0UL, 0UL, 1); (1UL, 1UL, 2) ] + ) + ) + + let m2 = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (0UL, 0UL, 10) ])) + + match map2LeftValues m1 m2 (fun a b -> Some(a + (defaultArg b 0))) with + | Error e -> failwithf "unexpected error %A" e + | Ok result -> + Assert.Equal(Ok(Some 11), get result 0UL 0UL) + Assert.Equal(Ok(Some 2), get result 1UL 1UL) + Assert.Equal(2UL, uint64 result.nvals) + +[] +let ``matrix map2i expands uniform leaves on both sides`` () = + let m1 = + SparseMatrix( + 4UL, + 4UL, + 4UL, + Storage(4UL, Matrix.qtree.Node(leaf_v 5, leaf_n (), leaf_n (), leaf_n ())) + ) + + let m2 = + SparseMatrix( + 4UL, + 4UL, + 4UL, + Storage(4UL, Matrix.qtree.Node(leaf_v 10, leaf_n (), leaf_n (), leaf_n ())) + ) + + let f i j a b = + match a, b with + | Some x, Some y -> Some(x + y + int (uint64 i) * 10 + int (uint64 j)) + | _ -> None + + match map2i m1 m2 f with + | Error e -> failwithf "unexpected error %A" e + | Ok result -> + Assert.Equal(4UL, uint64 result.nvals) + Assert.Equal(Ok(Some 15), get result 0UL 0UL) + Assert.Equal(Ok(Some 16), get result 0UL 1UL) + Assert.Equal(Ok(Some 25), get result 1UL 0UL) + Assert.Equal(Ok(Some 26), get result 1UL 1UL) + +[] +let ``matrix map2iValues applies indexed where both values present`` () = + let m1 = + fromCoordinateListUnchecked (CoordinateList(4UL, 4UL, [ (1UL, 1UL, 2) ])) + + let m2 = + fromCoordinateListUnchecked ( + CoordinateList( + 4UL, + 4UL, + [ (1UL, 1UL, 10); (2UL, 2UL, 20) ] + ) + ) + + let f i j a b = + Some(a + b + int (uint64 i) + int (uint64 j)) + + match map2iValues m1 m2 f with + | Error e -> failwithf "unexpected error %A" e + | Ok result -> + Assert.Equal(Ok(Some 14), get result 1UL 1UL) + Assert.Equal(Ok(None), get result 2UL 2UL) + Assert.Equal(1UL, uint64 result.nvals) + +[] +let ``matrix map2iAtLeastOne passes indices and side`` () = + let m1 = + fromCoordinateListUnchecked (CoordinateList(2UL, 2UL, [ (0UL, 0UL, 1) ])) + + let m2 = + fromCoordinateListUnchecked ( + CoordinateList( + 2UL, + 2UL, + [ (0UL, 0UL, 10); (1UL, 1UL, 20) ] + ) + ) + + let f i j = + function + | AtLeastOne.Both(a, b) -> Some(a + b + int (uint64 i) + int (uint64 j)) + | AtLeastOne.Left a -> Some(a) + | AtLeastOne.Right b -> Some(b + int (uint64 i) * 100 + int (uint64 j)) + + match map2iAtLeastOne m1 m2 f with + | Error e -> failwithf "unexpected error %A" e + | Ok result -> + Assert.Equal(Ok(Some 11), get result 0UL 0UL) + Assert.Equal(Ok(Some 121), get result 1UL 1UL) + +[] +let ``matrix map2iLeftValues applies indexed where left present`` () = + let m1 = + fromCoordinateListUnchecked (CoordinateList(2UL, 2UL, [ (1UL, 1UL, 2) ])) + + let m2 = + fromCoordinateListUnchecked (CoordinateList(2UL, 2UL, [ (1UL, 1UL, 10) ])) + + let f i j a b = + Some(a + (defaultArg b 0) + int (uint64 i) * 10 + int (uint64 j)) + + match map2iLeftValues m1 m2 f with + | Error e -> failwithf "unexpected error %A" e + | Ok result -> + Assert.Equal(Ok(Some 23), get result 1UL 1UL) + Assert.Equal(1UL, uint64 result.nvals) + +let ``fromCoordinateListUnchecked with out-of-range coordinates returns Error`` () = let coo = CoordinateList(6UL, 6UL, [ (9UL, 9UL, 13) ]) diff --git a/QuadTree.Tests/Tests.Vector.fs b/QuadTree.Tests/Tests.Vector.fs index ad6ad17..86fda5f 100644 --- a/QuadTree.Tests/Tests.Vector.fs +++ b/QuadTree.Tests/Tests.Vector.fs @@ -6,6 +6,11 @@ open System open Vector open Common +let private fromCoordinateListUnchecked (lst: Vector.CoordinateList<'a>) = + match Vector.fromCoordinateList lst with + | Ok v -> v + | Error e -> failwith e + let printVector (vector: SparseVector<_>) = printfn "Vector:" printfn " Length: %A" vector.length @@ -910,10 +915,483 @@ let ``Sort sorted vector`` () = let ``Init vector`` () = let expected = Vector.CoordinateList(3UL, [ (0UL, 0); (1UL, 1); (2UL, 2) ]) - |> Vector.fromCoordinateList + |> fromCoordinateListUnchecked let actual = Vector.init 3UL (fun i -> Some(int i)) - Assert.Equal(expected, Ok actual) + Assert.Equal(expected, actual) + + +[] +let ``Vector.filter all pass, none changed`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + let actual = Vector.filter v (fun x -> x > 0) + Assert.Equal(v, actual) + +[] +let ``Vector.filter none pass, all reset to zero`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + let expected = Vector.empty 8UL + let actual = Vector.filter v (fun x -> x < 0) + + Assert.Equal(expected, actual) + +[] +let ``Vector.filter length is not a power of 2, all reset to zero`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + let expected = Vector.empty 6UL + let actual = Vector.filter v (fun x -> x < 0) + + Assert.Equal(expected, actual) + +[] +let ``Vector.filter some pass, odd set to zero`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + let expected = + fromCoordinateListUnchecked ( + CoordinateList(8UL, [ (1UL, 2); (3UL, 4); (5UL, 6); (7UL, 8) ]) + ) + + let actual = Vector.filter v (fun x -> ((x % 2) = 0)) + + Assert.Equal(expected, actual) + +[] +let ``Vector.filter length is not a power of 2, not all reset to zero`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + let expected = + fromCoordinateListUnchecked ( + CoordinateList(6UL, [ (1UL, 2); (3UL, 4); (5UL, 6) ]) + ) + + let actual = Vector.filter v (fun x -> ((x % 2) = 0)) + + Assert.Equal(expected, actual) + +[] +let ``Vector.filter none pass, length is not a power of 2, none changed`` () = + let v = fromCoordinateListUnchecked (CoordinateList(6UL, [])) + let actual = Vector.filter v (fun x -> x > 0) + Assert.Equal(v, actual) + +[] +let ``Vector.filter none pass, length is a power of 2, none changed`` () = + let v = fromCoordinateListUnchecked (CoordinateList(8UL, [])) + let actual = Vector.filter v (fun x -> x > 0) + Assert.Equal(v, actual) + +[] +let ``Vector.filter single element, passes`` () = + let v = + fromCoordinateListUnchecked (CoordinateList(1UL, [ (0UL, 1) ])) + + let expected = + fromCoordinateListUnchecked (CoordinateList(1UL, [ (0UL, 1) ])) + + let actual = Vector.filter v (fun x -> x > 0) + + Assert.Equal(expected, actual) + +[] +let ``Vector.filter single element, fails`` () = + let v = + fromCoordinateListUnchecked (CoordinateList(1UL, [ (0UL, 1) ])) + + let expected = Vector.empty 1UL + let actual = Vector.filter v (fun x -> x < 0) + + Assert.Equal(expected, actual) + +[] +let ``Vector.exists the first element fits, length is a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + Assert.True(Vector.exists v (fun x -> x > 0)) + +[] +let ``Vector.exists the first element fits, length is not a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + Assert.True(Vector.exists v (fun x -> x > 0)) + +[] +let ``Vector.exists the last item fits, length is a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + Assert.True(Vector.exists v (fun x -> x > 7)) + +[] +let ``Vector.exists the last item fits, length is not a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + Assert.True(Vector.exists v (fun x -> x > 5)) + +[] +let ``Vector.exists empty list`` () = + let v = fromCoordinateListUnchecked (CoordinateList(8UL, [])) + Assert.False(Vector.exists v (fun x -> x > 0)) + +[] +let ``Vector.exists no matching elements, length is a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + Assert.False(Vector.exists v (fun x -> x > 8)) + +[] +let ``Vector.exists no matching elements, length is not a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + Assert.False(Vector.exists v (fun x -> x > 6)) + +[] +let ``Vector.exists single element matches`` () = + let v = + fromCoordinateListUnchecked (CoordinateList(1UL, [ (0UL, 1) ])) + + Assert.True(Vector.exists v (fun x -> x = 1)) + +[] +let ``Vector.exists single element does not match`` () = + let v = + fromCoordinateListUnchecked (CoordinateList(1UL, [ (0UL, 1) ])) + + Assert.False(Vector.exists v (fun x -> x = 2)) + +[] +let ``Vector.exists all elements match, length is a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + Assert.True(Vector.exists v (fun x -> x > 0)) + +[] +let ``Vector.exists all elements match, length is not a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + Assert.True(Vector.exists v (fun x -> x > 0)) + +[] +let ``Vector.forall the first element not fits, length is a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + Assert.False(Vector.forall v (fun x -> x > 1)) + +[] +let ``Vector.forall the first element not fits, length is not a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + Assert.False(Vector.forall v (fun x -> x > 1)) + +[] +let ``Vector.forall the last item not fits, length is a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + Assert.False(Vector.forall v (fun x -> x < 8)) + +[] +let ``Vector.forall the last item not fits, length is not a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + Assert.False(Vector.forall v (fun x -> x < 6)) + +[] +let ``Vector.forall empty list`` () = + let v = fromCoordinateListUnchecked (CoordinateList(8UL, [])) + Assert.True(Vector.forall v (fun x -> x > 0)) + +[] +let ``Vector.forall no matching elements, length is a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + Assert.False(Vector.forall v (fun x -> x > 8)) + +[] +let ``Vector.forall no matching elements, length is not a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + Assert.False(Vector.forall v (fun x -> x > 6)) + +[] +let ``Vector.forall single element matches`` () = + let v = + fromCoordinateListUnchecked (CoordinateList(1UL, [ (0UL, 1) ])) + + Assert.True(Vector.forall v (fun x -> x = 1)) + +[] +let ``Vector.forall single element does not match`` () = + let v = + fromCoordinateListUnchecked (CoordinateList(1UL, [ (0UL, 1) ])) + + Assert.False(Vector.forall v (fun x -> x = 2)) + +[] +let ``Vector.forall all elements match, length is a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 8UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) + (6UL, 7) + (7UL, 8) ] + ) + ) + + Assert.True(Vector.forall v (fun x -> x > 0)) + +[] +let ``Vector.forall all elements match, length is not a power of 2`` () = + let v = + fromCoordinateListUnchecked ( + CoordinateList( + 6UL, + [ (0UL, 1) + (1UL, 2) + (2UL, 3) + (3UL, 4) + (4UL, 5) + (5UL, 6) ] + ) + ) + + Assert.True(Vector.forall v (fun x -> x > 0)) [] let ``map on empty vector returns empty vector`` () = diff --git a/QuadTree/COOArray.fs b/QuadTree/COOArray.fs new file mode 100644 index 0000000..d5f1782 --- /dev/null +++ b/QuadTree/COOArray.fs @@ -0,0 +1,484 @@ +module COOArray + +open Common +open Matrix + +let private compareEntriesByRowCol (e1: COOEntry<'v>) (e2: COOEntry<'v>) = + let (i1, j1, _) = e1 + let (i2, j2, _) = e2 + let c = compare i1 i2 + if c <> 0 then c else compare j1 j2 + +let private entryComparer<'v> = + { new System.Collections.Generic.IComparer> with + member _.Compare(e1, e2) = compareEntriesByRowCol e1 e2 } + +let private iterCells + (nrows: uint64) + (ncols: uint64) + (action: uint64 -> uint64 -> unit) + = + let mutable i = 0UL + + while i < uint64 nrows do + let ri = i * 1UL + let mutable j = 0UL + + while j < uint64 ncols do + let cj = j * 1UL + action ri cj + j <- j + 1UL + + i <- i + 1UL + +let cooGet (coo: ArrayCOO<'a>, rowindex: uint64, colindex: uint64) : Result, Error> = + if uint64 rowindex >= uint64 coo.nrows then + raise (System.ArgumentOutOfRangeException("rowindex", "Row index is outside the matrix bounds.")) + elif uint64 colindex >= uint64 coo.ncols then + raise (System.ArgumentOutOfRangeException("colindex", "Column index is outside the matrix bounds.")) + else + let idx = + System.Array.BinarySearch(coo.list, (rowindex, colindex, Unchecked.defaultof<'a>), entryComparer<'a>) + + if idx >= 0 then + let (_, _, value) = coo.list.[idx] + Ok(Some value) + else + Ok None + +let cooUpdate + (coo: ArrayCOO<'a>, rowindex: uint64, colindex: uint64, value: 'a) + : Result, Error> = + if uint64 rowindex >= uint64 coo.nrows then + raise (System.ArgumentOutOfRangeException("rowindex", "Row index is outside the matrix bounds.")) + elif uint64 colindex >= uint64 coo.ncols then + raise (System.ArgumentOutOfRangeException("colindex", "Column index is outside the matrix bounds.")) + else + let idx = + System.Array.BinarySearch(coo.list, (rowindex, colindex, value), entryComparer<'a>) + + if idx >= 0 then + let arr = Array.copy coo.list + arr.[idx] <- (rowindex, colindex, value) + Ok(ArrayCOO.Create(coo.nrows, coo.ncols, arr)) + else + let insertAt = ~~~idx + let arr = Array.zeroCreate (coo.list.Length + 1) + + Array.blit coo.list 0 arr 0 insertAt + arr.[insertAt] <- (rowindex, colindex, value) + Array.blit coo.list insertAt arr (insertAt + 1) (coo.list.Length - insertAt) + + Ok(ArrayCOO.Create(coo.nrows, coo.ncols, arr)) + +let private cooMapInner (coo: ArrayCOO<'a>) (op: UnaryOp<'a, 'b>) : ArrayCOO<'b> = + match op with + | UnaryOp.ValuesOnly f -> + let buf = ResizeArray>(coo.list.Length) + + for (i, j, v) in coo.list do + match f v with + | Some r -> buf.Add((i, j, r)) + | None -> () + + ArrayCOO.Create(coo.nrows, coo.ncols, buf.ToArray()) + | UnaryOp.ValuesOnlyIndexed f -> + let buf = ResizeArray>(coo.list.Length) + + for (i, j, v) in coo.list do + match f i j v with + | Some r -> buf.Add((i, j, r)) + | None -> () + + ArrayCOO.Create(coo.nrows, coo.ncols, buf.ToArray()) + | UnaryOp.AllCells f -> + match f None with + | None -> + let buf = ResizeArray>(coo.list.Length) + + for (i, j, v) in coo.list do + match f (Some v) with + | Some r -> buf.Add((i, j, r)) + | None -> () + + ArrayCOO.Create(coo.nrows, coo.ncols, buf.ToArray()) + | Some fnone -> + let buf = ResizeArray>() + let mutable ptr = 0 + + iterCells coo.nrows coo.ncols (fun ri cj -> + let v = + if ptr < coo.list.Length then + let (ei, ej, ev) = coo.list.[ptr] + + if ei = ri && ej = cj then + ptr <- ptr + 1 + Some ev + else + None + else + None + + match v with + | Some value -> f (Some value) + | None -> Some fnone + |> Option.iter (fun value -> buf.Add((ri, cj, value)))) + + ArrayCOO.Create(coo.nrows, coo.ncols, buf.ToArray()) + | UnaryOp.AllCellsIndexed f -> + let buf = ResizeArray>() + let mutable ptr = 0 + + iterCells coo.nrows coo.ncols (fun ri cj -> + let v = + if ptr < coo.list.Length then + let (ei, ej, ev) = coo.list.[ptr] + + if ei = ri && ej = cj then + ptr <- ptr + 1 + Some ev + else + None + else + None + + f ri cj v |> Option.iter (fun value -> buf.Add((ri, cj, value)))) + + ArrayCOO.Create(coo.nrows, coo.ncols, buf.ToArray()) + +let private mergeBinary (a1: COOEntry<'a>[]) (a2: COOEntry<'b>[]) (op: BinaryOp<'a, 'b, 'c>) : COOEntry<'c>[] = + let buf = ResizeArray>(a1.Length + a2.Length) + let mutable p1 = 0 + let mutable p2 = 0 + + let emit i j v1 v2 = + match applyBinary op i j v1 v2 with + | Some r -> buf.Add((i, j, r)) + | None -> () + + while p1 < a1.Length || p2 < a2.Length do + if p1 >= a1.Length then + let (i, j, v2) = a2.[p2] + emit i j None (Some v2) + p2 <- p2 + 1 + elif p2 >= a2.Length then + let (i, j, v1) = a1.[p1] + emit i j (Some v1) None + p1 <- p1 + 1 + else + let (i1, j1, v1) = a1.[p1] + let (i2, j2, v2) = a2.[p2] + + if i1 = i2 && j1 = j2 then + emit i1 j1 (Some v1) (Some v2) + p1 <- p1 + 1 + p2 <- p2 + 1 + elif (i1, j1) < (i2, j2) then + emit i1 j1 (Some v1) None + p1 <- p1 + 1 + else + emit i2 j2 None (Some v2) + p2 <- p2 + 1 + + buf.ToArray() + +let private cooMap2Inner + (coo1: ArrayCOO<'a>) + (coo2: ArrayCOO<'b>) + (op: BinaryOp<'a, 'b, 'c>) + : Result, Error> = + if uint64 coo1.nrows <> uint64 coo2.nrows || uint64 coo1.ncols <> uint64 coo2.ncols then + Error Error.InconsistentSizeOfArguments + else + let nrows = coo1.nrows + let ncols = coo1.ncols + + let result = + match op with + | BinaryOp.AllCells f -> + match f None None with + | None -> mergeBinary coo1.list coo2.list op + | Some _ -> + let buf = ResizeArray>() + let mutable p1 = 0 + let mutable p2 = 0 + + iterCells nrows ncols (fun ri cj -> + let v1 = + if p1 < coo1.list.Length then + let (ei, ej, ev) = coo1.list.[p1] + + if ei = ri && ej = cj then + p1 <- p1 + 1 + Some ev + else + None + else + None + + let v2 = + if p2 < coo2.list.Length then + let (ei, ej, ev) = coo2.list.[p2] + + if ei = ri && ej = cj then + p2 <- p2 + 1 + Some ev + else + None + else + None + + match f v1 v2 with + | Some value -> buf.Add((ri, cj, value)) + | None -> ()) + + buf.ToArray() + | BinaryOp.AllCellsIndexed f -> + let buf = ResizeArray>() + let mutable p1 = 0 + let mutable p2 = 0 + + iterCells nrows ncols (fun ri cj -> + let v1 = + if p1 < coo1.list.Length then + let (ei, ej, ev) = coo1.list.[p1] + + if ei = ri && ej = cj then + p1 <- p1 + 1 + Some ev + else + None + else + None + + let v2 = + if p2 < coo2.list.Length then + let (ei, ej, ev) = coo2.list.[p2] + + if ei = ri && ej = cj then + p2 <- p2 + 1 + Some ev + else + None + else + None + + f ri cj v1 v2 |> Option.iter (fun value -> buf.Add((ri, cj, value)))) + + buf.ToArray() + | _ -> mergeBinary coo1.list coo2.list op + + ArrayCOO.Create(nrows, ncols, result) |> Ok + +let cooMap (coo: ArrayCOO<'a>) f = cooMapInner coo (UnaryOp.AllCells f) + +let cooMapValues (coo: ArrayCOO<'a>) f = cooMapInner coo (UnaryOp.ValuesOnly f) + +let cooMapi (coo: ArrayCOO<'a>) f = + cooMapInner coo (UnaryOp.AllCellsIndexed f) + +let cooMapiValues (coo: ArrayCOO<'a>) f = + cooMapInner coo (UnaryOp.ValuesOnlyIndexed f) + +let cooMap2 (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AllCells f) + +let cooMap2Values (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.ValuesOnly f) + +let cooMap2AllCells (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AllCells f) + +let cooMap2AtLeastOne (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AtLeastOneValue f) + +let cooMap2LeftValues (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.LeftValuesOnly f) + +let cooMap2i (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AllCellsIndexed f) + +let cooMap2iValues (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.ValuesOnlyIndexed f) + +let cooMap2iAllCells (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AllCellsIndexed f) + +let cooMap2iAtLeastOne (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AtLeastOneValueIndexed f) + +let cooMap2iLeftValues (coo1: ArrayCOO<'a>) (coo2: ArrayCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.LeftValuesOnlyIndexed f) + +let mxmcoo + (op_add: 'c option -> 'c option -> 'c option) + (op_mult: 'a option -> 'b option -> 'c option) + (m1: ArrayCOO<'a>) + (m2: ArrayCOO<'b>) + = + if uint64 m1.ncols <> uint64 m2.nrows then + Error Error.InconsistentSizeOfArguments + else + let valuesOf (arr: COOEntry<'v>[]) = + arr |> Array.map (fun (_, _, v) -> v) |> Array.distinct + + let canOptimize = + let noneNone = op_mult None None = None + + let multSomeNone = + valuesOf m1.list |> Array.forall (fun v -> op_mult (Some v) None = None) + + let multNoneSome = + valuesOf m2.list |> Array.forall (fun v -> op_mult None (Some v) = None) + + noneNone && multSomeNone && multNoneSome + + let generalResult () = + let m1Map = + let mutable m = Map.empty + + for (i, j, v) in m1.list do + m <- Map.add (i, j) v m + + m + + let m2Map = + let mutable m = Map.empty + + for (i, j, v) in m2.list do + m <- Map.add (i, j) v m + + m + + let kCount = uint64 m1.ncols + let result = ResizeArray>() + + iterCells m1.nrows m2.ncols (fun ri cj -> + let mutable acc = None + let mutable k = 0UL + + while k < kCount do + let key1 = (ri, k * 1UL) + let key2 = (k * 1UL, cj) + acc <- op_add acc (op_mult (Map.tryFind key1 m1Map) (Map.tryFind key2 m2Map)) + k <- k + 1UL + + match acc with + | Some value -> result.Add((ri, cj, value)) + | None -> ()) + + ArrayCOO.Create(m1.nrows, m2.ncols, result.ToArray()) + + if canOptimize then + let rowStarts2 = ResizeArray>() + let rowBegins2 = ResizeArray() + let rowEnds2 = ResizeArray() + let mutable idx = 0 + + while idx < m2.list.Length do + let (r, _, _) = m2.list.[idx] + let beginIdx = idx + let mutable advance = true + + while idx < m2.list.Length && advance do + let (r', _, _) = m2.list.[idx] + + if r' = r then idx <- idx + 1 else advance <- false + + rowStarts2.Add(r) + rowBegins2.Add(beginIdx) + rowEnds2.Add(idx) + + let findRowSegments (k: uint64) = + let mutable lo = 0 + let mutable hi = rowStarts2.Count - 1 + let mutable foundMid = -1 + + while lo <= hi && foundMid < 0 do + let mid = (lo + hi) / 2 + + if rowStarts2.[mid] = k then foundMid <- mid + elif rowStarts2.[mid] < k then lo <- mid + 1 + else hi <- mid - 1 + + if foundMid >= 0 then + Some(rowBegins2.[foundMid], rowEnds2.[foundMid]) + else + None + + let products = ResizeArray>() + let mutable baseIdx = 0 + + while baseIdx < m1.list.Length do + let (row, _, _) = m1.list.[baseIdx] + let mutable nextIdx = baseIdx + 1 + let mutable advance = true + + while nextIdx < m1.list.Length && advance do + let (row', _, _) = m1.list.[nextIdx] + + if row' = row then + nextIdx <- nextIdx + 1 + else + advance <- false + + for e in baseIdx .. nextIdx - 1 do + let (_, k, v1) = m1.list.[e] + let kAsRow = uint64 k * 1UL + + match findRowSegments kAsRow with + | Some(beginIdx, endIdx) -> + for q in beginIdx .. endIdx - 1 do + let (_, j, v2) = m2.list.[q] + + match op_mult (Some v1) (Some v2) with + | Some product -> products.Add((row, j, product)) + | None -> () + | None -> () + + baseIdx <- nextIdx + + let canMerge = + let productValues = + products |> Seq.map (fun (_, _, v) -> v) |> Seq.distinct |> Seq.toArray + + productValues + |> Array.forall (fun v -> op_add (Some v) None = Some v && op_add None (Some v) = Some v) + + if canMerge then + let sortedProducts = + products.ToArray() + |> Array.sortWith (fun (i1, j1, _) (i2, j2, _) -> + let c = compare i1 i2 + + if c <> 0 then c else compare j1 j2) + + let result = ResizeArray>() + let mutable q = 0 + + while q < sortedProducts.Length do + let (i, j, v) = sortedProducts.[q] + let mutable sum = Some v + let mutable qq = q + 1 + let mutable advance = true + + while qq < sortedProducts.Length && advance do + let (i', j', v') = sortedProducts.[qq] + + if i' = i && j' = j then + sum <- op_add sum (Some v') + qq <- qq + 1 + else + advance <- false + + match sum with + | Some value -> result.Add((i, j, value)) + | None -> () + + q <- qq + + ArrayCOO.Create(m1.nrows, m2.ncols, result.ToArray()) |> Ok + else + generalResult () |> Ok + else + generalResult () |> Ok diff --git a/QuadTree/COOList.fs b/QuadTree/COOList.fs new file mode 100644 index 0000000..934c94b --- /dev/null +++ b/QuadTree/COOList.fs @@ -0,0 +1,368 @@ +module COOList + +open Common +open Matrix + +let private range (count: uint64) = + if count = 0UL then [] else [ 0UL .. count - 1UL ] + +[] +type ListCOO<'value> = + val nrows: uint64 + val ncols: uint64 + val entries: COOEntry<'value> list + + new(_nrows, _ncols, _entries: COOEntry<'value> list) = + { nrows = _nrows + ncols = _ncols + entries = _entries } + +let fromArray (coo: ArrayCOO<'a>) : ListCOO<'a> = + ListCOO<'a>(coo.nrows, coo.ncols, Array.toList coo.list) + +let toArray (coo: ListCOO<'a>) : ArrayCOO<'a> = + ArrayCOO.Create(coo.nrows, coo.ncols, Array.ofList coo.entries) + +let cooGet (coo: ListCOO<'a>, rowindex: uint64, colindex: uint64) : Result, Error> = + if uint64 rowindex >= uint64 coo.nrows then + raise (System.ArgumentOutOfRangeException("rowindex", "Row index is outside the matrix bounds.")) + elif uint64 colindex >= uint64 coo.ncols then + raise (System.ArgumentOutOfRangeException("colindex", "Column index is outside the matrix bounds.")) + else + match coo.entries |> List.tryFind (fun (i, j, _) -> i = rowindex && j = colindex) with + | Some(_, _, value) -> Ok(Some value) + | None -> Ok None + +let cooUpdate + (coo: ListCOO<'a>, rowindex: uint64, colindex: uint64, value: 'a) + : Result, Error> = + if uint64 rowindex >= uint64 coo.nrows then + raise (System.ArgumentOutOfRangeException("rowindex", "Row index is outside the matrix bounds.")) + elif uint64 colindex >= uint64 coo.ncols then + raise (System.ArgumentOutOfRangeException("colindex", "Column index is outside the matrix bounds.")) + else + let mutable acc = [] + let mutable rest = coo.entries + let mutable inserted = false + + while rest <> [] && not inserted do + let (i, j, v) = rest.Head + + if i = rowindex && j = colindex then + acc <- (rowindex, colindex, value) :: acc + rest <- rest.Tail + inserted <- true + elif rowindex < i || (rowindex = i && colindex < j) then + acc <- (rowindex, colindex, value) :: acc + inserted <- true + else + acc <- (i, j, v) :: acc + rest <- rest.Tail + + if not inserted then + acc <- (rowindex, colindex, value) :: acc + + while rest <> [] do + let entry = rest.Head + acc <- entry :: acc + rest <- rest.Tail + + Ok(ListCOO<'a>(coo.nrows, coo.ncols, List.rev acc)) + +let private cooMapInner (coo: ListCOO<'a>) (op: UnaryOp<'a, 'b>) : ListCOO<'b> = + let result = + match op with + | UnaryOp.ValuesOnly f -> + coo.entries + |> List.choose (fun (i, j, v) -> f v |> Option.map (fun r -> (i, j, r))) + | UnaryOp.ValuesOnlyIndexed f -> + coo.entries + |> List.choose (fun (i, j, v) -> f i j v |> Option.map (fun r -> (i, j, r))) + | UnaryOp.AllCells f -> + match f None with + | None -> + coo.entries + |> List.choose (fun (i, j, v) -> f (Some v) |> Option.map (fun r -> (i, j, r))) + | Some fnone -> + let mutable rest = coo.entries + + [ for i in range (uint64 coo.nrows) do + let ri = i * 1UL + + for j in range (uint64 coo.ncols) do + let cj = j * 1UL + + let value = + match rest with + | (ei, ej, ev) :: tail when ei = ri && ej = cj -> + rest <- tail + Some ev + | _ -> None + + let res = + match value with + | Some v -> f (Some v) + | None -> Some fnone + + match res with + | Some value -> yield (ri, cj, value) + | None -> () ] + | UnaryOp.AllCellsIndexed f -> + let mutable rest = coo.entries + + [ for i in range (uint64 coo.nrows) do + let ri = i * 1UL + + for j in range (uint64 coo.ncols) do + let cj = j * 1UL + + let value = + match rest with + | (ei, ej, ev) :: tail when ei = ri && ej = cj -> + rest <- tail + Some ev + | _ -> None + + match f ri cj value with + | Some value -> yield (ri, cj, value) + | None -> () ] + + ListCOO<'b>(coo.nrows, coo.ncols, result) + +let private mergeBinary (l1: COOEntry<'a> list) (l2: COOEntry<'b> list) (op: BinaryOp<'a, 'b, 'c>) : COOEntry<'c> list = + let mutable acc = [] + let mutable rest1 = l1 + let mutable rest2 = l2 + + let emit i j v1 v2 = + match applyBinary op i j v1 v2 with + | Some r -> acc <- (i, j, r) :: acc + | None -> () + + while rest1 <> [] || rest2 <> [] do + match rest1, rest2 with + | [], [] -> () + | (i, j, v1) :: t1, [] -> + emit i j (Some v1) None + rest1 <- t1 + | [], (i, j, v2) :: t2 -> + emit i j None (Some v2) + rest2 <- t2 + | (i1, j1, v1) :: t1, (i2, j2, v2) :: t2 -> + if i1 = i2 && j1 = j2 then + emit i1 j1 (Some v1) (Some v2) + rest1 <- t1 + rest2 <- t2 + elif (i1, j1) < (i2, j2) then + emit i1 j1 (Some v1) None + rest1 <- t1 + else + emit i2 j2 None (Some v2) + rest2 <- t2 + + List.rev acc + +let private cooMap2Inner + (coo1: ListCOO<'a>) + (coo2: ListCOO<'b>) + (op: BinaryOp<'a, 'b, 'c>) + : Result, Error> = + if uint64 coo1.nrows <> uint64 coo2.nrows || uint64 coo1.ncols <> uint64 coo2.ncols then + Error Error.InconsistentSizeOfArguments + else + let nrows = coo1.nrows + let ncols = coo1.ncols + + let result = + match op with + | BinaryOp.AllCells f -> + match f None None with + | None -> mergeBinary coo1.entries coo2.entries op + | Some _ -> + let mutable rest1 = coo1.entries + let mutable rest2 = coo2.entries + + [ for i in range (uint64 nrows) do + let ri = i * 1UL + + for j in range (uint64 ncols) do + let cj = j * 1UL + + let v1 = + match rest1 with + | (ei, ej, ev) :: tail when ei = ri && ej = cj -> + rest1 <- tail + Some ev + | _ -> None + + let v2 = + match rest2 with + | (ei, ej, ev) :: tail when ei = ri && ej = cj -> + rest2 <- tail + Some ev + | _ -> None + + match f v1 v2 with + | Some value -> yield (ri, cj, value) + | None -> () ] + | BinaryOp.AllCellsIndexed f -> + let mutable rest1 = coo1.entries + let mutable rest2 = coo2.entries + + [ for i in range (uint64 nrows) do + let ri = i * 1UL + + for j in range (uint64 ncols) do + let cj = j * 1UL + + let v1 = + match rest1 with + | (ei, ej, ev) :: tail when ei = ri && ej = cj -> + rest1 <- tail + Some ev + | _ -> None + + let v2 = + match rest2 with + | (ei, ej, ev) :: tail when ei = ri && ej = cj -> + rest2 <- tail + Some ev + | _ -> None + + match f ri cj v1 v2 with + | Some value -> yield (ri, cj, value) + | None -> () ] + | _ -> mergeBinary coo1.entries coo2.entries op + + ListCOO<'c>(nrows, ncols, result) |> Ok + +let cooMap (coo: ListCOO<'a>) f = cooMapInner coo (UnaryOp.AllCells f) + +let cooMapValues (coo: ListCOO<'a>) f = cooMapInner coo (UnaryOp.ValuesOnly f) + +let cooMapi (coo: ListCOO<'a>) f = + cooMapInner coo (UnaryOp.AllCellsIndexed f) + +let cooMapiValues (coo: ListCOO<'a>) f = + cooMapInner coo (UnaryOp.ValuesOnlyIndexed f) + +let cooMap2 (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AllCells f) + +let cooMap2Values (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.ValuesOnly f) + +let cooMap2AllCells (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AllCells f) + +let cooMap2AtLeastOne (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AtLeastOneValue f) + +let cooMap2LeftValues (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.LeftValuesOnly f) + +let cooMap2i (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AllCellsIndexed f) + +let cooMap2iValues (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.ValuesOnlyIndexed f) + +let cooMap2iAllCells (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AllCellsIndexed f) + +let cooMap2iAtLeastOne (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.AtLeastOneValueIndexed f) + +let cooMap2iLeftValues (coo1: ListCOO<'a>) (coo2: ListCOO<'b>) f = + cooMap2Inner coo1 coo2 (BinaryOp.LeftValuesOnlyIndexed f) + +let mxmcoo + (op_add: 'c option -> 'c option -> 'c option) + (op_mult: 'a option -> 'b option -> 'c option) + (m1: ListCOO<'a>) + (m2: ListCOO<'b>) + = + if uint64 m1.ncols <> uint64 m2.nrows then + Error Error.InconsistentSizeOfArguments + else + let entries1 = m1.entries + let entries2 = m2.entries + + let valuesOf (entries: (uint64 * uint64 * 'v) list) = + entries |> List.map (fun (_, _, v) -> v) |> List.distinct + + let canOptimize = + let noneNone = op_mult None None = None + + let multSomeNone = + valuesOf entries1 |> List.forall (fun v -> op_mult (Some v) None = None) + + let multNoneSome = + valuesOf entries2 |> List.forall (fun v -> op_mult None (Some v) = None) + + noneNone && multSomeNone && multNoneSome + + let generalResult () = + let m1Map = entries1 |> List.map (fun (i, j, v) -> ((i, j), v)) |> Map.ofList + let m2Map = entries2 |> List.map (fun (i, j, v) -> ((i, j), v)) |> Map.ofList + let kCount = uint64 m1.ncols + + let result = + [ for i in range (uint64 m1.nrows) do + let ri = i * 1UL + + for j in range (uint64 m2.ncols) do + let cj = j * 1UL + + let products = + [ for k in range kCount do + let a = m1Map |> Map.tryFind (ri, k * 1UL) + let b = m2Map |> Map.tryFind (k * 1UL, cj) + yield op_mult a b ] + + let sum = products |> List.fold (fun acc p -> op_add acc p) None + + match sum with + | Some value -> yield (ri, cj, value) + | None -> () ] + + ListCOO<'c>(m1.nrows, m2.ncols, result) + + if canOptimize then + let m1ByRow = entries1 |> List.groupBy (fun (i, _, _) -> i) |> Map.ofList + let m2ByRow = entries2 |> List.groupBy (fun (k, _, _) -> k) |> Map.ofList + + let result = + [ for KeyValue(i, m1Entries) in m1ByRow do + for (_, k, v1) in m1Entries do + let kAsRow = uint64 k * 1UL + + match m2ByRow |> Map.tryFind kAsRow with + | Some m2Entries -> + for (_, j, v2) in m2Entries do + match op_mult (Some v1) (Some v2) with + | Some product -> yield (i, j, product) + | None -> () + | None -> () ] + + let canMerge = + let productValues = result |> List.map (fun (_, _, v) -> v) |> List.distinct + + productValues + |> List.forall (fun v -> op_add (Some v) None = Some v && op_add None (Some v) = Some v) + + if canMerge then + let grouped = + result + |> List.groupBy (fun (i, j, _) -> (i, j)) + |> List.map (fun ((i, j), entries) -> + let sum = entries |> List.map (fun (_, _, v) -> Some v) |> List.reduce op_add + (i, j, sum)) + |> List.choose (fun (i, j, v) -> v |> Option.map (fun v -> (i, j, v))) + |> List.sortBy (fun (i, j, _) -> (i, j)) + + ListCOO<'c>(m1.nrows, m2.ncols, grouped) |> Ok + else + generalResult () |> Ok + else + generalResult () |> Ok diff --git a/QuadTree/Matrix.fs b/QuadTree/Matrix.fs index 86abf03..fc150b4 100644 --- a/QuadTree/Matrix.fs +++ b/QuadTree/Matrix.fs @@ -1,4 +1,4 @@ -module Matrix +module Matrix open Common @@ -54,6 +54,12 @@ type rowindex [] type colindex +let getQuadrantCoords (pr, pc) halfSize = + (pr, pc), // NORTH WEST + (pr, pc + halfSize * 1UL), // NORTH EAST + (pr + halfSize * 1UL, pc), // SOUTH WEST + (pr + halfSize * 1UL, pc + halfSize * 1UL) // SOUTH EAST + type COOEntry<'value> = uint64 * uint64 * 'value [] @@ -67,11 +73,43 @@ type CoordinateList<'value> = ncols = _ncols list = _list } -let private getQuadrantCoords (pr, pc) halfSize = - (pr, pc), // NORTH WEST - (pr, pc + halfSize * 1UL), // NORTH EAST - (pr + halfSize * 1UL, pc), // SOUTH WEST - (pr + halfSize * 1UL, pc + halfSize * 1UL) // SOUTH EAST +[] +type ArrayCOO<'value> = + val nrows: uint64 + val ncols: uint64 + val list: COOEntry<'value>[] + + new(_nrows, _ncols, _list: COOEntry<'value> seq) = + let sorted = + _list + |> Seq.toArray + |> Array.sortWith (fun (i1, j1, _) (i2, j2, _) -> + let c = compare i1 i2 + if c <> 0 then c else compare j1 j2) + + { nrows = _nrows + ncols = _ncols + list = sorted } + + new(_nrows, _ncols, _list: COOEntry<'value>[], _presorted: bool) = + let sorted = + if _presorted then + _list + else + _list + |> Array.sortWith (fun (i1, j1, _) (i2, j2, _) -> + let c = compare i1 i2 + if c <> 0 then c else compare j1 j2) + + { nrows = _nrows + ncols = _ncols + list = sorted } + + // Fast factory: does NOT re-sort, expects an already sorted array. + // Used by COOArray operations whose results are built in (row, col) order and + // by cooUpdate, which maintains the sorted invariant itself. + static member Create(nrows: uint64, ncols: uint64, entries: COOEntry<'value>[]) : ArrayCOO<'value> = + ArrayCOO<'value>(nrows, ncols, entries, true) let fromCoordinateList (coo: CoordinateList<'a>) = let unique = @@ -160,183 +198,366 @@ let empty nrows ncols = SparseMatrix(nrows, ncols, 0UL, Storage(storageSize, Leaf Dummy)) -let map (matrix: SparseMatrix<'a>) f = - let rec inner (size: uint64) (tree: qtree>) = - match tree with - | Node(nw, ne, sw, se) -> - let nwTree, nwNvals = inner (size / 2UL) nw - let neTree, neNvals = inner (size / 2UL) ne - let swTree, swNvals = inner (size / 2UL) sw - let seTree, seNvals = inner (size / 2UL) se - (mkNode nwTree neTree swTree seTree), nwNvals + neNvals + swNvals + seNvals - | Leaf(Dummy) -> Leaf(Dummy), 0UL - | Leaf(UserValue(v)) -> - let res = f v - - let nnz = - match res with - | None -> 0UL - | _ -> (uint64 size) * (uint64 size) * 1UL - - Leaf(UserValue(res)), nnz - - let newTree, newNvals = inner matrix.storage.size matrix.storage.data - SparseMatrix(matrix.nrows, matrix.ncols, newNvals, Storage(matrix.storage.size, newTree)) - -let map2 (matrix1: SparseMatrix<_>) (matrix2: SparseMatrix<_>) f = - let rec inner (size: uint64) matrix1 matrix2 = - let _do x1 x2 x3 x4 y1 y2 y3 y4 = - let new_size = size / 2UL - - match (inner new_size x1 y1), (inner new_size x2 y2), (inner new_size x3 y3), (inner new_size x4 y4) with - | Ok((new_t1, nvals1)), Ok((new_t2, nvals2)), Ok((new_t3, nvals3)), Ok((new_t4, nvals4)) -> - ((mkNode new_t1 new_t2 new_t3 new_t4), nvals1 + nvals2 + nvals3 + nvals4) |> Ok - | Error(e), _, _, _ - | _, Error(e), _, _ - | _, _, Error(e), _ - | _, _, _, Error(e) -> Error(e) - - match (matrix1, matrix2) with - | Leaf(Dummy), Leaf(Dummy) -> Ok(Leaf(Dummy), 0UL) - | Leaf(UserValue(v1)), Leaf(UserValue(v2)) -> - let res = f v1 v2 - - let nnz = - match res with - | None -> 0UL - | _ -> (uint64 size) * (uint64 size) * 1UL +let get (matrix: SparseMatrix<'a>) (row: uint64) (col: uint64) : Result, Error> = + if uint64 row >= uint64 matrix.nrows then + raise (System.ArgumentOutOfRangeException("row", "Row index is outside the matrix bounds.")) + elif uint64 col >= uint64 matrix.ncols then + raise (System.ArgumentOutOfRangeException("col", "Column index is outside the matrix bounds.")) + else + let rec inner tree (pr: uint64) (pc: uint64) (size: uint64) = + match tree with + | Leaf Dummy -> None + | Leaf(UserValue v) -> v + | Node(nw, ne, sw, se) -> + let halfSize = size / 2UL + let midR = pr + halfSize * 1UL + let midC = pc + halfSize * 1UL - (Leaf(UserValue(res)), nnz) |> Ok + if uint64 row < uint64 midR then + if uint64 col < uint64 midC then + inner nw pr pc halfSize + else + inner ne pr midC halfSize + else if uint64 col < uint64 midC then + inner sw midR pc halfSize + else + inner se midR midC halfSize - | Node(x1, x2, x3, x4), Node(y1, y2, y3, y4) -> _do x1 x2 x3 x4 y1 y2 y3 y4 - | Node(x1, x2, x3, x4), Leaf(v) -> _do x1 x2 x3 x4 matrix2 matrix2 matrix2 matrix2 - | Leaf(v), Node(x1, x2, x3, x4) -> _do matrix1 matrix1 matrix1 matrix1 x1 x2 x3 x4 - | (x, y) -> Error Error.InconsistentStructureOfStorages + Ok(inner matrix.storage.data (0UL) (0UL) (uint64 matrix.storage.size)) - if matrix1.nrows = matrix2.nrows && matrix1.ncols = matrix2.ncols then - match inner matrix1.storage.size matrix1.storage.data matrix2.storage.data with - | Error x -> Error x - | Ok(storage, nvals) -> - (SparseMatrix(matrix1.nrows, matrix1.ncols, nvals, (Storage(matrix1.storage.size, storage)))) - |> Ok +let set + (matrix: SparseMatrix<'a>) + (row: uint64) + (col: uint64) + (value: 'a) + : Result, Error> = + if uint64 row >= uint64 matrix.nrows then + raise (System.ArgumentOutOfRangeException("row", "Row index is outside the matrix bounds.")) + elif uint64 col >= uint64 matrix.ncols then + raise (System.ArgumentOutOfRangeException("col", "Column index is outside the matrix bounds.")) else - Error Error.InconsistentSizeOfArguments - -let map2i (matrix1: SparseMatrix<_>) (matrix2: SparseMatrix<_>) f = - let rec inner (prow: uint64) (pcol: uint64) (size: uint64) matrix1 matrix2 = - match (matrix1, matrix2) with - | Node(x1, x2, x3, x4), Node(y1, y2, y3, y4) -> + let rec inner tree (pr: uint64) (pc: uint64) (size: uint64) = let halfSize = size / 2UL - let (nwR, nwC), (neR, neC), (swR, swC), (seR, seC) = - getQuadrantCoords (prow, pcol) (uint64 halfSize) + if size = 1UL then + match tree with + | Leaf(UserValue oldVal) -> + let newVal = Some value + + let delta = + match newVal, oldVal with + | Some _, None -> 1L + | None, Some _ -> -1L + | _ -> 0L + + Leaf(UserValue newVal), delta + | Leaf Dummy -> Leaf(UserValue(Some value)), 1L + | _ -> failwith "Unreachable" + else + let midR = pr + halfSize * 1UL + let midC = pc + halfSize * 1UL + + let (nw, ne, sw, se) = + match tree with + | Node(nw, ne, sw, se) -> nw, ne, sw, se + | Leaf v -> Leaf v, Leaf v, Leaf v, Leaf v + + let newChild, delta = + if uint64 row < uint64 midR then + if uint64 col < uint64 midC then + inner nw pr pc halfSize + else + inner ne pr midC halfSize + else if uint64 col < uint64 midC then + inner sw midR pc halfSize + else + inner se midR midC halfSize - let t1, nvals1 = inner nwR nwC halfSize x1 y1 - let t2, nvals2 = inner neR neC halfSize x2 y2 - let t3, nvals3 = inner swR swC halfSize x3 y3 - let t4, nvals4 = inner seR seC halfSize x4 y4 - (mkNode t1 t2 t3 t4), nvals1 + nvals2 + nvals3 + nvals4 - | Node(x1, x2, x3, x4), Leaf(v2) -> + if uint64 row < uint64 midR then + if uint64 col < uint64 midC then + mkNode newChild ne sw se, delta + else + mkNode nw newChild sw se, delta + else if uint64 col < uint64 midC then + mkNode nw ne newChild se, delta + else + mkNode nw ne sw newChild, delta + + let storage, deltaNNZ = + inner matrix.storage.data (0UL) (0UL) (uint64 matrix.storage.size) + + let nvals = uint64 (int64 matrix.nvals + deltaNNZ) * 1UL + Ok(SparseMatrix(matrix.nrows, matrix.ncols, nvals, Storage(matrix.storage.size, storage))) + +type UnaryOp<'a, 'b> = + | ValuesOnly of ('a -> Option<'b>) + | ValuesOnlyIndexed of (uint64 -> uint64 -> 'a -> Option<'b>) + | AllCells of (Option<'a> -> Option<'b>) + | AllCellsIndexed of (uint64 -> uint64 -> Option<'a> -> Option<'b>) + +let private mapInner (matrix: SparseMatrix<'a>) (op: UnaryOp<'a, 'b>) : SparseMatrix<'b> = + let rec inner + (prow: uint64) + (pcol: uint64) + (size: uint64) + (tree: qtree>) + : qtree> * uint64 = + match tree with + | Node(nw, ne, sw, se) -> let halfSize = size / 2UL let (nwR, nwC), (neR, neC), (swR, swC), (seR, seC) = getQuadrantCoords (prow, pcol) (uint64 halfSize) - let t1, nvals1 = inner nwR nwC halfSize x1 (Leaf(v2)) - let t2, nvals2 = inner neR neC halfSize x2 (Leaf(v2)) - let t3, nvals3 = inner swR swC halfSize x3 (Leaf(v2)) - let t4, nvals4 = inner seR seC halfSize x4 (Leaf(v2)) - (mkNode t1 t2 t3 t4), nvals1 + nvals2 + nvals3 + nvals4 - | Leaf(v1), Node(y1, y2, y3, y4) -> - let halfSize = size / 2UL + let t1, nvals1 = inner nwR nwC halfSize nw + let t2, nvals2 = inner neR neC halfSize ne + let t3, nvals3 = inner swR swC halfSize sw + let t4, nvals4 = inner seR seC halfSize se - let (nwR, nwC), (neR, neC), (swR, swC), (seR, seC) = - getQuadrantCoords (prow, pcol) (uint64 halfSize) + mkNode t1 t2 t3 t4, nvals1 + nvals2 + nvals3 + nvals4 + | Leaf(Dummy) -> Leaf(Dummy), 0UL + | Leaf(UserValue(v)) -> + match op with + | UnaryOp.ValuesOnly f -> + match v with + | None -> Leaf(UserValue(None)), 0UL + | Some v' -> + let res = f v' + + let nvals = + if res.IsSome then + (uint64 size) * (uint64 size) * 1UL + else + 0UL + + Leaf(UserValue(res)), nvals + | UnaryOp.ValuesOnlyIndexed f -> + match v with + | None -> Leaf(UserValue(None)), 0UL + | Some v' -> + if size = 1UL then + let res = f prow pcol v' + let nvals = if res.IsSome then 1UL else 0UL + Leaf(UserValue(res)), nvals + else + let halfSize = size / 2UL + + let (nwR, nwC), (neR, neC), (swR, swC), (seR, seC) = + getQuadrantCoords (prow, pcol) (uint64 halfSize) + + let t1, nvals1 = inner nwR nwC halfSize (Leaf(UserValue(v))) + let t2, nvals2 = inner neR neC halfSize (Leaf(UserValue(v))) + let t3, nvals3 = inner swR swC halfSize (Leaf(UserValue(v))) + let t4, nvals4 = inner seR seC halfSize (Leaf(UserValue(v))) + mkNode t1 t2 t3 t4, nvals1 + nvals2 + nvals3 + nvals4 + | UnaryOp.AllCells f -> + let res = f v + + let nvals = + if res.IsSome then + (uint64 size) * (uint64 size) * 1UL + else + 0UL + + Leaf(UserValue(res)), nvals + | UnaryOp.AllCellsIndexed f -> + if size = 1UL then + let res = f prow pcol v + let nvals = if res.IsSome then 1UL else 0UL + Leaf(UserValue(res)), nvals + else + let halfSize = size / 2UL - let t1, nvals1 = inner nwR nwC halfSize (Leaf(v1)) y1 - let t2, nvals2 = inner neR neC halfSize (Leaf(v1)) y2 - let t3, nvals3 = inner swR swC halfSize (Leaf(v1)) y3 - let t4, nvals4 = inner seR seC halfSize (Leaf(v1)) y4 - (mkNode t1 t2 t3 t4), nvals1 + nvals2 + nvals3 + nvals4 - | Leaf(Dummy), Leaf(Dummy) -> Leaf(Dummy), 0UL - | Leaf(UserValue(v1)), Leaf(UserValue(v2)) -> - let res = f prow pcol v1 v2 + let (nwR, nwC), (neR, neC), (swR, swC), (seR, seC) = + getQuadrantCoords (prow, pcol) (uint64 halfSize) - let nnz = - match res with - | Some _ -> 1UL - | None -> 0UL + let t1, nvals1 = inner nwR nwC halfSize (Leaf(UserValue(v))) + let t2, nvals2 = inner neR neC halfSize (Leaf(UserValue(v))) + let t3, nvals3 = inner swR swC halfSize (Leaf(UserValue(v))) + let t4, nvals4 = inner seR seC halfSize (Leaf(UserValue(v))) + mkNode t1 t2 t3 t4, nvals1 + nvals2 + nvals3 + nvals4 - Leaf(UserValue(res)), nnz - | Leaf(UserValue(v)), Leaf(Dummy) -> - let res = f prow pcol v None + let storage, nvals = + inner 0UL 0UL matrix.storage.size matrix.storage.data - let nnz = - match res with - | Some _ -> 1UL - | None -> 0UL + SparseMatrix(matrix.nrows, matrix.ncols, nvals, Storage(matrix.storage.size, storage)) + +let map (matrix: SparseMatrix<_>) f = mapInner matrix (UnaryOp.AllCells f) + +let mapValues (matrix: SparseMatrix<'a>) f = mapInner matrix (UnaryOp.ValuesOnly f) + +type AtLeastOne<'a, 'b> = + | Both of 'a * 'b + | Left of 'a + | Right of 'b + +type BinaryOp<'a, 'b, 'c> = + | ValuesOnly of ('a -> 'b -> Option<'c>) + | ValuesOnlyIndexed of (uint64 -> uint64 -> 'a -> 'b -> Option<'c>) + | AllCells of (Option<'a> -> Option<'b> -> Option<'c>) + | AllCellsIndexed of (uint64 -> uint64 -> Option<'a> -> Option<'b> -> Option<'c>) + | AtLeastOneValue of (AtLeastOne<'a, 'b> -> Option<'c>) + | AtLeastOneValueIndexed of (uint64 -> uint64 -> AtLeastOne<'a, 'b> -> Option<'c>) + | LeftValuesOnly of ('a -> Option<'b> -> Option<'c>) + | LeftValuesOnlyIndexed of (uint64 -> uint64 -> 'a -> Option<'b> -> Option<'c>) + +let applyBinary + (op: BinaryOp<'a, 'b, 'c>) + (prow: uint64) + (pcol: uint64) + (v1: Option<'a>) + (v2: Option<'b>) + : Option<'c> = + match op with + | BinaryOp.ValuesOnly f -> + match v1, v2 with + | Some a, Some b -> f a b + | _ -> None + | BinaryOp.ValuesOnlyIndexed f -> + match v1, v2 with + | Some a, Some b -> f prow pcol a b + | _ -> None + | BinaryOp.AllCells f -> f v1 v2 + | BinaryOp.AllCellsIndexed f -> f prow pcol v1 v2 + | BinaryOp.AtLeastOneValue f -> + match v1, v2 with + | Some a, Some b -> f (AtLeastOne.Both(a, b)) + | Some a, None -> f (AtLeastOne.Left a) + | None, Some b -> f (AtLeastOne.Right b) + | None, None -> None + | BinaryOp.AtLeastOneValueIndexed f -> + match v1, v2 with + | Some a, Some b -> f prow pcol (AtLeastOne.Both(a, b)) + | Some a, None -> f prow pcol (AtLeastOne.Left a) + | None, Some b -> f prow pcol (AtLeastOne.Right b) + | None, None -> None + | BinaryOp.LeftValuesOnly f -> + match v1 with + | Some a -> f a v2 + | None -> None + | BinaryOp.LeftValuesOnlyIndexed f -> + match v1 with + | Some a -> f prow pcol a v2 + | None -> None + +let private isIndexedBinary (op: BinaryOp<'a, 'b, 'c>) = + match op with + | BinaryOp.ValuesOnlyIndexed _ + | BinaryOp.AllCellsIndexed _ + | BinaryOp.AtLeastOneValueIndexed _ + | BinaryOp.LeftValuesOnlyIndexed _ -> true + | _ -> false + +let private map2Inner + (matrix1: SparseMatrix<'a>) + (matrix2: SparseMatrix<'b>) + (op: BinaryOp<'a, 'b, 'c>) + : Result, Error> = + let rec inner + (prow: uint64) + (pcol: uint64) + (size: uint64) + (tree1: qtree>) + (tree2: qtree>) + : Result> * uint64, Error> = + let split + (x1: qtree>) + (x2: qtree>) + (x3: qtree>) + (x4: qtree>) + (y1: qtree>) + (y2: qtree>) + (y3: qtree>) + (y4: qtree>) + = + let halfSize = size / 2UL - Leaf(UserValue(res)), nnz - | Leaf(Dummy), Leaf(UserValue(v)) -> - let res = f prow pcol None v + let (nwR, nwC), (neR, neC), (swR, swC), (seR, seC) = + getQuadrantCoords (prow, pcol) (uint64 halfSize) - let nnz = - match res with - | Some _ -> 1UL - | None -> 0UL + match + (inner nwR nwC halfSize x1 y1), + (inner neR neC halfSize x2 y2), + (inner swR swC halfSize x3 y3), + (inner seR seC halfSize x4 y4) + with + | Ok(t1, nvals1), Ok(t2, nvals2), Ok(t3, nvals3), Ok(t4, nvals4) -> + Ok(mkNode t1 t2 t3 t4, nvals1 + nvals2 + nvals3 + nvals4) + | Error e, _, _, _ + | _, Error e, _, _ + | _, _, Error e, _ + | _, _, _, Error e -> Error e + + match tree1, tree2 with + | Node(x1, x2, x3, x4), Node(y1, y2, y3, y4) -> split x1 x2 x3 x4 y1 y2 y3 y4 + | Node(x1, x2, x3, x4), Leaf(v2) -> split x1 x2 x3 x4 (Leaf(v2)) (Leaf(v2)) (Leaf(v2)) (Leaf(v2)) + | Leaf(v1), Node(y1, y2, y3, y4) -> split (Leaf(v1)) (Leaf(v1)) (Leaf(v1)) (Leaf(v1)) y1 y2 y3 y4 + | Leaf(Dummy), Leaf(Dummy) -> Ok(Leaf(Dummy), 0UL) + | Leaf(UserValue(v1)), Leaf(UserValue(v2)) -> + if size > 1UL && isIndexedBinary op then + split + (Leaf(UserValue(v1))) + (Leaf(UserValue(v1))) + (Leaf(UserValue(v1))) + (Leaf(UserValue(v1))) + (Leaf(UserValue(v2))) + (Leaf(UserValue(v2))) + (Leaf(UserValue(v2))) + (Leaf(UserValue(v2))) + else + let res = applyBinary op prow pcol v1 v2 - Leaf(UserValue(res)), nnz + let nnz = + if res.IsSome then + (uint64 size) * (uint64 size) * 1UL + else + 0UL - if matrix1.nrows = matrix2.nrows && matrix1.ncols = matrix2.ncols then - let storage, nvals = - inner 0UL 0UL matrix1.storage.size matrix1.storage.data matrix2.storage.data + Ok(Leaf(UserValue(res)), nnz) + | _ -> Error Error.InconsistentStructureOfStorages - SparseMatrix(matrix1.nrows, matrix1.ncols, nvals, (Storage(matrix1.storage.size, storage))) - |> Ok + if matrix1.nrows = matrix2.nrows && matrix1.ncols = matrix2.ncols then + inner 0UL 0UL matrix1.storage.size matrix1.storage.data matrix2.storage.data + |> Result.map (fun (storage, nvals) -> + SparseMatrix(matrix1.nrows, matrix1.ncols, nvals, Storage(matrix1.storage.size, storage))) else Error Error.InconsistentSizeOfArguments -let mapi (matrix: SparseMatrix<'a>) f = - let rec inner (prow: uint64) (pcol: uint64) (size: uint64) matrix = - match matrix with - | Node(x1, x2, x3, x4) -> - let halfSize = size / 2UL +let map2 (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.AllCells f) - let (nwR, nwC), (neR, neC), (swR, swC), (seR, seC) = - getQuadrantCoords (prow, pcol) (uint64 halfSize) +let map2Values (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.ValuesOnly f) - let t1, nvals1 = inner nwR nwC halfSize x1 - let t2, nvals2 = inner neR neC halfSize x2 - let t3, nvals3 = inner swR swC halfSize x3 - let t4, nvals4 = inner seR seC halfSize x4 - (mkNode t1 t2 t3 t4), nvals1 + nvals2 + nvals3 + nvals4 - | Leaf(Dummy) -> Leaf(Dummy), 0UL - | Leaf(UserValue(v)) -> - if size = 1UL then - let res = f prow pcol v +let map2AllCells (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.AllCells f) - let nnz = - match res with - | Some _ -> 1UL - | None -> 0UL +let map2AtLeastOne (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.AtLeastOneValue f) - Leaf(UserValue(res)), nnz - else - let halfSize = size / 2UL +let map2LeftValues (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.LeftValuesOnly f) - let (nwR, nwC), (neR, neC), (swR, swC), (seR, seC) = - getQuadrantCoords (prow, pcol) (uint64 halfSize) +let map2i (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.AllCellsIndexed f) - let t1, nvals1 = inner nwR nwC halfSize (Leaf(UserValue(v))) - let t2, nvals2 = inner neR neC halfSize (Leaf(UserValue(v))) - let t3, nvals3 = inner swR swC halfSize (Leaf(UserValue(v))) - let t4, nvals4 = inner seR seC halfSize (Leaf(UserValue(v))) - (mkNode t1 t2 t3 t4), nvals1 + nvals2 + nvals3 + nvals4 +let map2iValues (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.ValuesOnlyIndexed f) - let storage, nvals = - inner 0UL 0UL matrix.storage.size matrix.storage.data +let map2iAllCells (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.AllCellsIndexed f) - SparseMatrix(matrix.nrows, matrix.ncols, nvals, (Storage(matrix.storage.size, storage))) +let map2iAtLeastOne (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.AtLeastOneValueIndexed f) + +let map2iLeftValues (matrix1: SparseMatrix<'a>) (matrix2: SparseMatrix<'b>) f = + map2Inner matrix1 matrix2 (BinaryOp.LeftValuesOnlyIndexed f) + +let mapi (matrix: SparseMatrix<'a>) f = + mapInner matrix (UnaryOp.AllCellsIndexed f) + +let mapiValues (matrix: SparseMatrix<'a>) f = + mapInner matrix (UnaryOp.ValuesOnlyIndexed f) let foldAssociative (folder: 'T option -> 'T option -> 'T option) (state: 'T option) (matrix: SparseMatrix<'T>) = let rec traverse tree (size: uint64) (state: 'T option) = @@ -438,6 +659,54 @@ let transpose (matrix: SparseMatrix<_>) = let mask (m1: SparseMatrix<'a>) (m2: SparseMatrix<'b>) f = map2 m1 m2 (fun m1 m2 -> if f m2 then m1 else None) + +let filter (matrix: SparseMatrix<'a>) (predicate: 'a -> bool) : SparseMatrix<'a> = + let rec inner (prow: uint64) (pcol: uint64) (size: uint64) matrix = + match matrix with + | Node(x1, x2, x3, x4) -> + let halfSize = size / 2UL + + let (nwR, nwC), (neR, neC), (swR, swC), (seR, seC) = + getQuadrantCoords (prow, pcol) (uint64 halfSize) + + let t1, nvals1 = inner nwR nwC halfSize x1 + let t2, nvals2 = inner neR neC halfSize x2 + let t3, nvals3 = inner swR swC halfSize x3 + let t4, nvals4 = inner seR seC halfSize x4 + (mkNode t1 t2 t3 t4), nvals1 + nvals2 + nvals3 + nvals4 + | Leaf(Dummy) -> Leaf(Dummy), 0UL + | Leaf(UserValue(None)) -> Leaf(UserValue(None)), 0UL + | Leaf(UserValue(Some(v))) -> + if predicate v then + Leaf(UserValue(Some v)), (uint64 size) * (uint64 size) * 1UL + else + Leaf(UserValue(None)), 0UL + + let storage, nvals = + inner 0UL 0UL matrix.storage.size matrix.storage.data + + SparseMatrix(matrix.nrows, matrix.ncols, nvals, (Storage(matrix.storage.size, storage))) + +let exists (matrix: SparseMatrix<'a>) (predicate: 'a -> bool) : bool = + let rec inner tree = + match tree with + | Leaf(Dummy) -> false + | Leaf(UserValue(None)) -> false + | Leaf(UserValue(Some(v))) -> predicate v + | Node(nw, ne, sw, se) -> inner nw || inner ne || inner sw || inner se + + inner matrix.storage.data + +let forall (matrix: SparseMatrix<'a>) (predicate: 'a -> bool) : bool = + let rec inner tree = + match tree with + | Leaf(Dummy) -> true + | Leaf(UserValue(None)) -> true + | Leaf(UserValue(Some(v))) -> predicate v + | Node(nw, ne, sw, se) -> inner nw && inner ne && inner sw && inner se + + inner matrix.storage.data + let slice (matrix: SparseMatrix<'a>) (rowStart: int) diff --git a/QuadTree/QuadTree.fsproj b/QuadTree/QuadTree.fsproj index abfc6ec..b909d5c 100644 --- a/QuadTree/QuadTree.fsproj +++ b/QuadTree/QuadTree.fsproj @@ -9,6 +9,8 @@ + + diff --git a/QuadTree/Vector.fs b/QuadTree/Vector.fs index 0981063..760d215 100644 --- a/QuadTree/Vector.fs +++ b/QuadTree/Vector.fs @@ -576,6 +576,44 @@ let scatter (Ok w) | Error x -> Error Error.InconsistentStructureOfStorages +let filter (vector: SparseVector<'a>) (predicate: 'a -> bool) : SparseVector<'a> = + let rec inner (size: uint64) vector = + match vector with + | Node(x1, x2) -> + let t1, nvals1 = inner (size / 2UL) x1 + let t2, nvals2 = inner (size / 2UL) x2 + (mkNode t1 t2), nvals1 + nvals2 + | Leaf(Dummy) -> Leaf(Dummy), 0UL + | Leaf(UserValue(None)) -> Leaf(UserValue(None)), 0UL + | Leaf(UserValue(Some(v))) -> + if predicate v then + Leaf(UserValue(Some(v))), (uint64 size) * 1UL + else + Leaf(UserValue(None)), 0UL + + let storage, nvals = inner vector.storage.size vector.storage.data + SparseVector(vector.length, nvals, Storage(vector.storage.size, storage)) + +let exists (vector: SparseVector<'a>) (predicate: 'a -> bool) : bool = + let rec inner vector = + match vector with + | Leaf(Dummy) -> false + | Leaf(UserValue(None)) -> false + | Leaf(UserValue(Some(v))) -> predicate v + | Node(x1, x2) -> inner x1 || inner x2 + + inner vector.storage.data + +let forall (vector: SparseVector<'a>) (predicate: 'a -> bool) : bool = + let rec inner vector = + match vector with + | Leaf(Dummy) -> true + | Leaf(UserValue(None)) -> true + | Leaf(UserValue(Some(v))) -> predicate v + | Node(x1, x2) -> inner x1 && inner x2 + + inner vector.storage.data + let slice (_start: int) (_end: int) (vector: SparseVector<'a>) : Result, string> = if _start < 0 then Error "Start should be >= 0"