Merge pull request #34 from li3zhen1/dev

Bugfix: Comparing floating points to .nan

DocPostprocess.swift

@@ -70,6 +70,12 @@ do {
     if !fileManager.fileExists(atPath: iconDestPath) {
         try fileManager.copyItem(atPath: iconSourcePath, toPath: iconDestPath)
     }
+    for moduleName in moduleNames {
+        let iconDestPath = "./docs/\(moduleName)/favicon.png"
+        if !fileManager.fileExists(atPath: iconDestPath) {
+            try fileManager.copyItem(atPath: iconSourcePath, toPath: iconDestPath)
+        }
+    }
 
 } catch {
     // Handle errors by printing to the console for now

README.md

@@ -210,9 +210,9 @@ See [Example](https://github.com/li3zhen1/Grape/tree/main/Examples/ForceDirected
 
 #### Simulation
 
-Grape uses simd to calculate position and velocity. Currently it takes ~0.010 seconds to iterate 120 times over the example graph(2D). (77 vertices, 254 edges, with manybody, center, collide and link forces. Release build on a M1 Max, tested with command `swift test -c release`)
+Grape uses simd to calculate position and velocity. Currently it takes **~0.010** seconds to iterate 120 times over the example graph(2D). (77 vertices, 254 edges, with manybody, center, collide and link forces. Release build on a M1 Max, tested with command `swift test -c release`)
 
-For 3D simulation, it takes ~0.015 seconds for the same graph and same configs.
+For 3D simulation, it takes **~0.014** seconds for the same graph and same configs.
 
 > [!IMPORTANT]
 > Due to heavy use of generics (some are not specialized in Debug mode), the performance in Debug build is ~100x slower than Release build. Grape might ship a version with pre-inlined generics to address this problem.

Sources/ForceSimulation/KDTree/KDTree.swift

@@ -77,7 +77,7 @@ where
         return index
     }
 
-    /// Expand the current node towards a direction. 
+    /// Expand the current node towards a direction.
     ///
     /// The expansion will double the size on each dimension. Then the data in delegate will be copied to the new children.
     /// - Parameter direction: An Integer between 0 and `directionCount - 1`, where `directionCount` equals to 2^(dimension of the vector).
@@ -138,7 +138,7 @@ where
         cover(point)
         addWithoutCover(nodeIndex, at: point)
     }
-    
+
     @inlinable
     public mutating func addWithoutCover(_ nodeIndex: NodeIndex, at point: Vector) {
         defer {
@@ -149,9 +149,9 @@ where
                 nodeIndices.append(nodeIndex)
                 nodePosition = point
                 return
-            } else if nodePosition == point
-                || nodePosition!.distanceSquared(to: point) < clusterDistanceSquared
-            {
+            } else if nodePosition!.distanceSquared(to: point) < clusterDistanceSquared {
+                // the condition (nodePosition == point) is mostly only true when the tree is initialized
+                // hence omitted
                 nodeIndices.append(nodeIndex)
                 return
             } else {
@@ -303,7 +303,9 @@ extension KDTree {
     /// Visit the tree in pre-order.
     ///
     /// - Parameter shouldVisitChildren: a closure that returns a boolean value indicating whether should continue to visit children.
-    @inlinable public mutating func visit(shouldVisitChildren: (inout KDTree<Vector, Delegate>) -> Bool) {
+    @inlinable public mutating func visit(
+        shouldVisitChildren: (inout KDTree<Vector, Delegate>) -> Bool
+    ) {
         if shouldVisitChildren(&self) && children != nil {
             // this is an internal node
             for i in children!.indices {
@@ -314,7 +316,6 @@ extension KDTree {
 
 }
 
-
 // public struct KDTreeRoot<Vector, Delegate, Property>
 // where
 //     Vector: SimulatableVector & L2NormCalculatable,
@@ -332,10 +333,9 @@ extension KDTree {
 //         self.propertyBuffer = propertyBuffer
 //     }
 
-
 //     @inlinable
 //     public mutating func add(_ nodeIndex: Int, at point: Vector) {
 //         root.cover(point)
 //         root.addWithoutCover(nodeIndex, at: point)
 //     }
-// }
+// }

Sources/ForceSimulation/Kinetics.swift

@@ -120,6 +120,7 @@ where Vector: SimulatableVector & L2NormCalculatable {
             count: position.count,
             initialValue: .zero
         )
+
         self.velocity = UnsafeArray<Vector>.createBuffer(
             withHeader: position.count,
             count: position.count,
@@ -135,6 +136,13 @@ where Vector: SimulatableVector & L2NormCalculatable {
         self.randomGenerator.initialize(to: .init(seed: randomSeed))
     }
 
+    @inlinable
+    internal func jigglePosition() {
+        for i in range {
+            position[i] = position[i].jiggled(by: self.randomGenerator)
+        }
+    }
+
     @inlinable
     static func createZeros(
         links: [EdgeID<Int>],
@@ -182,17 +190,7 @@ extension Kinetics {
     func updateAlpha() {
         alpha += (alphaTarget - alpha) * alphaDecay
     }
-
-    @inlinable
-    func invalidateRange(_ range: Range<Int>) {
-        fatalError("Not implemented")
-    }
-
-    @inlinable
-    func validateRangeAndExtendIfNeccessary(_ range: Range<Int>) {
-        fatalError("Not implemented")
-    }
-
+
 }
 
 public typealias Kinetics2D = Kinetics<SIMD2<Double>>

Sources/ForceSimulation/LinearCongruentialGenerator.swift

@@ -17,23 +17,24 @@ public struct DoubleLinearCongruentialGenerator: DeterministicRandomGenerator {
     @usableFromInline internal static let c: UInt32 = 1_013_904_223
     @usableFromInline internal static var _s: UInt32 = 1
     @usableFromInline internal var s: UInt32 = 1
+    @usableFromInline internal static let m: Double = 4_294_967_296
 
     @inlinable public mutating func next() -> Double {
         // Perform the linear congruential generation with integer arithmetic.
         // The overflow addition and multiplication automatically wrap around,
         // thus imitating the modulo operation.
-        s = Self.a &* s &+ Self.c
+        s = (Self.a &* s) &+ Self.c
 
         // Convert the result to Double and divide by m to normalize it.
-        return Double(s) / 4_294_967_296.0
+        return Double(s) / Self.m
     }
 
     @inlinable public static func next() -> Double {
 
-        Self._s = Self.a &* Self._s &+ Self.c
+        Self._s = (Self.a &* Self._s) &+ Self.c
 
         // Convert the result to Double and divide by m to normalize it.
-        return Double(Self._s) / 4_294_967_296.0
+        return Double(Self._s) / Self.m
     }
 
     @inlinable public init(seed: OverflowingInteger) {
@@ -52,21 +53,22 @@ public struct FloatLinearCongruentialGenerator: DeterministicRandomGenerator {
     @usableFromInline internal static let c: UInt16 = 74
     @usableFromInline internal static var _s: UInt16 = 1
     @usableFromInline internal var s: UInt16 = 1
+    @usableFromInline internal static let m: Float = 65537.0
 
     @inlinable public mutating func next() -> Float {
         // Perform the linear congruential generation with integer arithmetic.
         // The overflow addition and multiplication automatically wrap around.
-        s = Self.a &* s &+ Self.c
+        s = (Self.a &* s) &+ Self.c
 
         // Convert the result to Float and divide by m to normalize it.
-        return Float(s) / 65537.0
+        return Float(s) / Self.m
     }
 
     @inlinable public static func next() -> Float {
-        _s = a &* _s &+ c
+        _s = (a &* _s) &+ c
 
         // Convert the result to Float and divide by m to normalize it.
-        return Float(_s) / 65537.0
+        return Float(_s) / Self.m
     }
 
     @inlinable public init(seed: OverflowingInteger) {
@@ -92,18 +94,24 @@ extension Float: HasDeterministicRandomGenerator {
 }
 
 extension HasDeterministicRandomGenerator {
+
+    @inlinable
+    static var jigglingScale: Self {
+        return 1e-5
+    }
+
     @inlinable
     public func jiggled() -> Self {
-        if self == .zero || self == .nan {
-            return (Generator.next() - 0.5) * 1e-5
+        if self == .zero || self.isNaN {
+            return (Generator.next() - 0.5) * Self.jigglingScale
         }
         return self
     }
 
     @inlinable
     public func jiggled(by: UnsafeMutablePointer<Generator>) -> Self {
-        if self == .zero || self == .nan {
-            return (by.pointee.next() - 0.5) * 1e-5
+        if self == .zero || self.isNaN {
+            return (by.pointee.next() - 0.5) * Self.jigglingScale
         }
         return self
     }

Sources/ForceSimulation/SimulatableVector.swift

@@ -63,7 +63,7 @@ extension SIMD2: SimulatableVector where Scalar: FloatingPoint & HasDeterministi
 
     @inlinable
     public static var clusterDistanceSquared: Scalar {
-        return 1e-10
+        return clusterDistance * clusterDistance
     }
 }
 
@@ -76,7 +76,7 @@ extension SIMD3: SimulatableVector where Scalar: FloatingPoint & HasDeterministi
 
     @inlinable
     public static var clusterDistanceSquared: Scalar {
-        return 1e-10
+        return clusterDistance * clusterDistance
     }
 }
 

Sources/ForceSimulation/Simulation.swift

@@ -42,6 +42,7 @@ where Vector: SimulatableVector & L2NormCalculatable, ForceField: ForceProtocol<
             velocityDecay: velocityDecay,
             count: nodeCount
         )
+        self.kinetics.jigglePosition()
         forceField.bindKinetics(self.kinetics)
         self.forceField = forceField
     }