diff --git a/Examples/ForceDirectedGraph3D/ForceDirectedGraph3D/ContentView.swift b/Examples/ForceDirectedGraph3D/ForceDirectedGraph3D/ContentView.swift index 7c7a2af..6e19e89 100644 --- a/Examples/ForceDirectedGraph3D/ForceDirectedGraph3D/ContentView.swift +++ b/Examples/ForceDirectedGraph3D/ForceDirectedGraph3D/ContentView.swift @@ -8,7 +8,7 @@ import SwiftUI import RealityKit import RealityKitContent -import NDTree +import simd import ForceSimulation @@ -83,9 +83,9 @@ struct ContentView: View { let sim = buildSimulation() let positions = sim.nodePositions.map { pos in simd_float3( - Float(pos[1]) * scaleRatio, - -Float(pos[0]) * scaleRatio, - Float(pos[2]) * scaleRatio + 0.25 + (pos[1]) * scaleRatio, + -(pos[0]) * scaleRatio, + (pos[2]) * scaleRatio + 0.25 )} @@ -151,10 +151,6 @@ struct ContentView: View { -// let fromPosition: simd_float3 = [Float(_fromPosition[0]), Float(_fromPosition[1]), Float(_fromPosition[2])] -// let toPosition: simd_float3 = [Float(_toPosition[0]), Float(_toPosition[1]), Float(_toPosition[2])] - - let cylinderVector = toPosition - fromPosition diff --git a/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample.xcodeproj/project.pbxproj b/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample.xcodeproj/project.pbxproj index d97bf5d..053d14f 100644 --- a/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample.xcodeproj/project.pbxproj +++ b/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample.xcodeproj/project.pbxproj @@ -7,6 +7,7 @@ objects = { /* Begin PBXBuildFile section */ + B719E4112AE5FBFC009D6C24 /* ForceDirectedLatticeView.swift in Sources */ = {isa = PBXBuildFile; fileRef = B719E4102AE5FBFC009D6C24 /* ForceDirectedLatticeView.swift */; }; B7AFA55B2ADF4997009C7154 /* ForceDirectedGraphExampleApp.swift in Sources */ = {isa = PBXBuildFile; fileRef = B7AFA55A2ADF4997009C7154 /* ForceDirectedGraphExampleApp.swift */; }; B7AFA55D2ADF4997009C7154 /* ContentView.swift in Sources */ = {isa = PBXBuildFile; fileRef = B7AFA55C2ADF4997009C7154 /* ContentView.swift */; }; B7AFA55F2ADF4999009C7154 /* Assets.xcassets in Resources */ = {isa = PBXBuildFile; fileRef = B7AFA55E2ADF4999009C7154 /* Assets.xcassets */; }; @@ -17,6 +18,7 @@ /* End PBXBuildFile section */ /* Begin PBXFileReference section */ + B719E4102AE5FBFC009D6C24 /* ForceDirectedLatticeView.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; path = ForceDirectedLatticeView.swift; sourceTree = ""; }; B7AFA5572ADF4997009C7154 /* ForceDirectedGraphExample.app */ = {isa = PBXFileReference; explicitFileType = wrapper.application; includeInIndex = 0; path = ForceDirectedGraphExample.app; sourceTree = BUILT_PRODUCTS_DIR; }; B7AFA55A2ADF4997009C7154 /* ForceDirectedGraphExampleApp.swift */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.swift; path = ForceDirectedGraphExampleApp.swift; sourceTree = ""; }; B7AFA55C2ADF4997009C7154 /* ContentView.swift */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.swift; path = ContentView.swift; sourceTree = ""; }; @@ -58,6 +60,7 @@ B7AFA5592ADF4997009C7154 /* ForceDirectedGraphExample */ = { isa = PBXGroup; children = ( + B719E4102AE5FBFC009D6C24 /* ForceDirectedLatticeView.swift */, B7AFA55A2ADF4997009C7154 /* ForceDirectedGraphExampleApp.swift */, B7AFA55C2ADF4997009C7154 /* ContentView.swift */, B7AFA55E2ADF4999009C7154 /* Assets.xcassets */, @@ -154,6 +157,7 @@ buildActionMask = 2147483647; files = ( B7AFA55D2ADF4997009C7154 /* ContentView.swift in Sources */, + B719E4112AE5FBFC009D6C24 /* ForceDirectedLatticeView.swift in Sources */, B7AFA56F2ADF49D6009C7154 /* Data.swift in Sources */, B7AFA55B2ADF4997009C7154 /* ForceDirectedGraphExampleApp.swift in Sources */, ); diff --git a/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ContentView.swift b/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ContentView.swift index 2e0f9cf..a1fd042 100644 --- a/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ContentView.swift +++ b/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ContentView.swift @@ -6,7 +6,7 @@ // import SwiftUI - +import simd import NDTree import ForceSimulation import CoreGraphics @@ -35,11 +35,11 @@ struct MiserableNode: Identifiable { struct ContentView: View { - @State var points: [Vector2d] = [] + @State var points: [simd_double2] = [] var sim: Simulation2D let data: Miserable - var linkForce: LinkForce + var linkForce: Simulation2D.LinkForce init() { diff --git a/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ForceDirectedGraphExampleApp.swift b/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ForceDirectedGraphExampleApp.swift index a1a8685..83729b0 100644 --- a/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ForceDirectedGraphExampleApp.swift +++ b/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ForceDirectedGraphExampleApp.swift @@ -11,7 +11,8 @@ import SwiftUI struct ForceDirectedGraphExampleApp: App { var body: some Scene { WindowGroup { - ContentView().padding(0) +// ContentView().padding(0) + ForceDirectedLatticeView().padding(0)//.ignoresSafeArea() } } } diff --git a/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ForceDirectedLatticeView.swift b/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ForceDirectedLatticeView.swift index bdc4270..27773ff 100644 --- a/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ForceDirectedLatticeView.swift +++ b/Examples/ForceDirectedGraphExample/ForceDirectedGraphExample/ForceDirectedLatticeView.swift @@ -5,80 +5,82 @@ // Created by li3zhen1 on 10/18/23. // -import SwiftUI -import NDTree -import ForceSimulation import CoreGraphics +import ForceSimulation +import simd +import SwiftUI struct ForceDirectedLatticeView: View { - @State var points: [Vector2d]? = nil - + @State var points: [simd_double2]? = nil + private let sim: Simulation2D - private let edgeIds: [(Int,Int)] + private let edgeIds: [(Int, Int)] private let nodeIds: [Int] private let canvasWidth: CGFloat = 800.0 - let width = 30 - + let width = 36 + init() { - self.nodeIds = Array(0..<(width*width)) - + self.nodeIds = Array(0..<(width * width)) + var edge = [(Int, Int)]() for i in 0.. [!IMPORTANT] -> When working with 3D contents, you probably need `Float` types instead of `Double`. The example here manually cast `Double` to `Float`. -> The `float32` branch relaxes the generic constraint from `Scalar == Double` to `Scalar: ExpressibleByFloatLiteral`, which allows you to get out-of-box supports for other floating point types. However, the performance will be downgraded. -
@@ -68,7 +65,7 @@ Grape currently provides 2 packages, `NDTree` and `ForceSimulation`. The basic concepts of simulations and forces can be found here: [Force simulations - D3](https://d3js.org/d3-force/simulation). You can simply create 2D or 3D simulations by using `Simulation2D` or `Simulation3D`: ```swift -import NDTree +import simd import ForceSimulation struct Node: Identifiable { ... } @@ -98,7 +95,7 @@ See [Example](https://github.com/li3zhen1/Grape/tree/main/Examples/ForceDirected ### Advanced -Almost all the types in Grape work with any SIMD-like data structures. To integrate Grape into platforms where `import simd` isn't supported, you need to create a struct conforming to the `VectorLike` protocol. For ease of use, it's also recommended to add some type aliases. Here’s how you can do it: +Grape provides a set of generic based types that works with any SIMD-like data structures. To integrate Grape into platforms where `import simd` isn't supported, or higher dimensions, you need to create a struct conforming to the `VectorLike` protocol. For ease of use, it's also recommended to add some type aliases. Here’s how you can do it: ```swift /// All required implementations should have same semantics @@ -109,10 +106,11 @@ protocol HyperoctreeDelegate: NDTreeDelegate where V == SuperCool4DVector {} typealias HyperoctBox = NDBox typealias Hyperoctree = NDTree -typealias Simulation4D = Simulation +typealias Simulation4D = SimulationKD ``` -Also, this is how you create a 4D simulation with or without `simd_double4`. (Though I don't know what good it does) +> [!IMPORTANT] +> When using generic based types, you ***pay for dynamic dispatch***, in terms of performance. It's recommended to use `Simulation2D` or `Simulation3D` whenever possible.
@@ -137,9 +135,9 @@ Also, this is how you create a 4D simulation with or without `simd_double4`. (Th ## Performance -Grape uses simd to calculate position and velocity. Currently it takes ~0.12 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) +Grape uses simd to calculate position and velocity. Currently it takes ~0.05 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`) -Due to the iteration over simd lanes, going 3D will hurt performance. (~0.16 seconds for the same graph and same configs.) +Due to the iteration over simd lanes, going 3D will hurt performance. (~0.075 seconds for the same graph and same configs.)
diff --git a/Sources/ForceSimulation/ForceLike.swift b/Sources/ForceSimulation/ForceLike.swift index 76bbe07..75f3b59 100644 --- a/Sources/ForceSimulation/ForceLike.swift +++ b/Sources/ForceSimulation/ForceLike.swift @@ -10,19 +10,12 @@ import NDTree /// A protocol that represents a force. /// A force takes a simulation state and modifies its node positions and velocities. public protocol ForceLike { - associatedtype NodeID: Hashable /// Takes a simulation state and modifies its node positions and velocities. /// This is executed in each tick of the simulation. - func apply(alpha: Double) + func apply() } public protocol NDTreeBasedForceLike: ForceLike { associatedtype TD: NDTreeDelegate } - -extension Array where Element: NDTreeBasedForceLike { - public func combined() { - - } -} diff --git a/Sources/ForceSimulation/ForceSimulation.docc/Documentation.md b/Sources/ForceSimulation/ForceSimulation.docc/Documentation.md index ecd8731..f29e7e1 100644 --- a/Sources/ForceSimulation/ForceSimulation.docc/Documentation.md +++ b/Sources/ForceSimulation/ForceSimulation.docc/Documentation.md @@ -18,13 +18,27 @@ Run force simulation on any dimensions. * * ``Simulation2D`` * ``Simulation3D`` -* ``Simulation`` +* ``SimulationKD`` ### Creating forces in a simulation -* ``CenterForce`` -* ``CollideForce`` -* ``LinkForce`` -* ``ManyBodyForce`` -* ``DirectionForce`` -* ``RadialForce`` \ No newline at end of file +* ``Simulation2D/CenterForce`` +* ``Simulation2D/CollideForce`` +* ``Simulation2D/LinkForce`` +* ``Simulation2D/ManyBodyForce`` +* ``Simulation2D/DirectionForce`` +* ``Simulation2D/RadialForce`` + +* ``Simulation3D/CenterForce`` +* ``Simulation3D/CollideForce`` +* ``Simulation3D/LinkForce`` +* ``Simulation3D/ManyBodyForce`` +* ``Simulation3D/DirectionForce`` +* ``Simulation3D/RadialForce`` + +* ``SimulationKD/CenterForce`` +* ``SimulationKD/CollideForce`` +* ``SimulationKD/LinkForce`` +* ``SimulationKD/ManyBodyForce`` +* ``SimulationKD/DirectionForce`` +* ``SimulationKD/RadialForce`` \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation2D/Simulation2D.CenterForce.swift b/Sources/ForceSimulation/Simulation2D/Simulation2D.CenterForce.swift new file mode 100644 index 0000000..3e8ba5f --- /dev/null +++ b/Sources/ForceSimulation/Simulation2D/Simulation2D.CenterForce.swift @@ -0,0 +1,65 @@ +// +// CenterForce.swift +// +// +// Created by li3zhen1 on 10/16/23. +// + + +#if canImport(simd) +import NDTree +import simd + +extension Simulation2D { + /// A force that drives nodes towards the center. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Collide Force - D3](https://d3js.org/d3-force/collide). + final public class CenterForce: ForceLike + where NodeID: Hashable { + public typealias V = simd_double2 + + public var center: V + public var strength: V.Scalar + weak var simulation: Simulation2D? + + internal init(center: V, strength: V.Scalar) { + self.center = center + self.strength = strength + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + + var meanPosition = V.zero + for n in sim.nodePositions { + meanPosition += n //.position + } + let delta = meanPosition * (self.strength / V.Scalar(sim.nodePositions.count)) + + for i in sim.nodePositions.indices { + sim.nodePositions[i] -= delta + } + } + + } + + /// Create a center force that drives nodes towards the center. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Collide Force - D3](https://d3js.org/d3-force/collide). + /// - Parameters: + /// - center: The center of the force. + /// - strength: The strength of the force. + @discardableResult + public func createCenterForce(center: V, strength: V.Scalar = 0.1) -> CenterForce { + let f = CenterForce(center: center, strength: strength) + f.simulation = self + self.forces.append(f) + return f + } + +} + +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation2D/Simulation2D.CollideForce.swift b/Sources/ForceSimulation/Simulation2D/Simulation2D.CollideForce.swift new file mode 100644 index 0000000..9aa3e51 --- /dev/null +++ b/Sources/ForceSimulation/Simulation2D/Simulation2D.CollideForce.swift @@ -0,0 +1,220 @@ +// +// File.swift +// +// +// Created by li3zhen1 on 10/17/23. +// + +#if canImport(simd) +import NDTree +import simd + +struct MaxRadiusTreeDelegate2D: QuadtreeDelegate where NodeID: Hashable { + + public typealias V = simd_double2 + + public var maxNodeRadius: V.Scalar + + @usableFromInline var radiusProvider: (NodeID) -> V.Scalar + + mutating func didAddNode(_ nodeId: NodeID, at position: V) { + let p = radiusProvider(nodeId) + maxNodeRadius = max(maxNodeRadius, p) + } + + mutating func didRemoveNode(_ nodeId: NodeID, at position: V) { + if radiusProvider(nodeId) >= maxNodeRadius { + // 🤯 for Collide force, set to 0 is fine + // Otherwise you need to traverse the delegate again + maxNodeRadius = 0 + } + } + + func copy() -> MaxRadiusTreeDelegate2D { + return Self(maxNodeRadius: maxNodeRadius, radiusProvider: radiusProvider) + } + + func spawn() -> MaxRadiusTreeDelegate2D { + return Self(radiusProvider: radiusProvider) + } + + init(maxNodeRadius: V.Scalar = 0, radiusProvider: @escaping (NodeID) -> V.Scalar) { + self.maxNodeRadius = maxNodeRadius + self.radiusProvider = radiusProvider + } + +} + +extension Simulation2D { + + /// A force that prevents nodes from overlapping. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. + /// See [Collide Force - D3](https://d3js.org/d3-force/collide). + public final class CollideForce: ForceLike + where NodeID: Hashable { + + public typealias V = simd_double2 + + weak var simulation: Simulation2D? { + didSet { + guard let sim = simulation else { return } + self.calculatedRadius = radius.calculated(for: sim) + } + } + + public enum CollideRadius { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + public var radius: CollideRadius + var calculatedRadius: [V.Scalar] = [] + + public let iterationsPerTick: UInt + public var strength: V.Scalar + + internal init( + radius: CollideRadius, + strength: V.Scalar = 1.0, + iterationsPerTick: UInt = 1 + ) { + self.radius = radius + self.iterationsPerTick = iterationsPerTick + self.strength = strength + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + + for _ in 0...cover(of: sim.nodePositions) + + let clusterDistance: V.Scalar = V.Scalar(Int(0.00001)) + + let tree = Quadtree( + box: coveringBox, clusterDistance: clusterDistance + ) { + return switch self.radius { + case .constant(let m): + MaxRadiusTreeDelegate2D { _ in m } + case .varied(_): + MaxRadiusTreeDelegate2D { index in + self.calculatedRadius[index] + } + } + } + + for i in sim.nodePositions.indices { + tree.add(i, at: sim.nodePositions[i]) + } + + for i in sim.nodePositions.indices { + let iOriginalPosition = sim.nodePositions[i] + let iOriginalVelocity = sim.nodeVelocities[i] + let iR = self.calculatedRadius[i] + let iR2 = iR * iR + let iPosition = iOriginalPosition + iOriginalVelocity + + tree.visit { t in + + let maxRadiusOfQuad = t.delegate.maxNodeRadius + let deltaR = maxRadiusOfQuad + iR + + if t.nodePosition != nil { + for j in t.nodeIndices { + // print("\(i)<=>\(j)") + // is leaf, make sure every collision happens once. + guard j > i else { continue } + + let jR = self.calculatedRadius[j] + let jOriginalPosition = sim.nodePositions[j] + let jOriginalVelocity = sim.nodeVelocities[j] + var deltaPosition = + iPosition - (jOriginalPosition + jOriginalVelocity) + let l = simd_length_squared(deltaPosition) + + let deltaR = iR + jR + if l < deltaR * deltaR { + + var l = simd_length(deltaPosition.jiggled()) + l = (deltaR - l) / l * self.strength + + let jR2 = jR * jR + + let k = jR2 / (iR2 + jR2) + + deltaPosition *= l + + sim.nodeVelocities[i] += deltaPosition * k + sim.nodeVelocities[j] -= deltaPosition * (1 - k) + } + } + return false + } + + for laneIndex in t.box.p0.indices { + let _v = t.box.p0[laneIndex] + if _v > iPosition[laneIndex] + deltaR /* True if no overlap */ { + return false + } + } + + for laneIndex in t.box.p1.indices { + let _v = t.box.p1[laneIndex] + if _v < iPosition[laneIndex] - deltaR /* True if no overlap */ { + return false + } + } + return true + + // return + // !(t.quad.x0 > iPosition.x + deltaR /* True if no overlap */ + // || t.quad.x1 < iPosition.x - deltaR + // || t.quad.y0 > iPosition.y + deltaR + // || t.quad.y1 < iPosition.y - deltaR) + } + } + } + } + + } + + /// Create a collide force that prevents nodes from overlapping. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. + /// See [Collide Force - D3](https://d3js.org/d3-force/collide). + /// - Parameters: + /// - radius: The radius of the force. + /// - strength: The strength of the force. + /// - iterationsPerTick: The number of iterations per tick. + @discardableResult + public func createCollideForce( + radius: CollideForce.CollideRadius = .constant(3.0), + strength: V.Scalar = 1.0, + iterationsPerTick: UInt = 1 + ) -> CollideForce { + let f = CollideForce( + radius: radius, + strength: strength, + iterationsPerTick: iterationsPerTick + ) + f.simulation = self + self.forces.append(f) + return f + } + +} + +extension Simulation2D.CollideForce.CollideRadius { + public func calculated(for simulation: Simulation2D) -> [Double] { + switch self { + case .constant(let r): + return Array(repeating: r, count: simulation.nodePositions.count) + case .varied(let radiusProvider): + return simulation.nodeIds.map { radiusProvider($0) } + } + } +} +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation2D/Simulation2D.DirectionForce.swift b/Sources/ForceSimulation/Simulation2D/Simulation2D.DirectionForce.swift new file mode 100644 index 0000000..d6981d8 --- /dev/null +++ b/Sources/ForceSimulation/Simulation2D/Simulation2D.DirectionForce.swift @@ -0,0 +1,128 @@ +// +// PositionForce.swift +// +// +// Created by li3zhen1 on 10/1/23. +// + +#if canImport(simd) +import NDTree +import simd + + + +extension Simulation2D { + + /// A force that moves nodes to a target position. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + final public class DirectionForce2D: ForceLike + where NodeID: Hashable { + + public typealias V = simd_double2 + + public enum Direction { + case x + case y + case entryOfVector(Int) + } + public enum Strength { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + + public enum TargetOnDirection { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + + public var strength: Strength + public var direction: Int + public var calculatedStrength: [V.Scalar] = [] + public var targetOnDirection: TargetOnDirection + public var calculatedTargetOnDirection: [V.Scalar] = [] + + internal init( + direction: Direction, targetOnDirection: TargetOnDirection, + strength: Strength = .constant(1.0) + ) { + + self.strength = strength + self.direction = direction.lane + self.targetOnDirection = targetOnDirection + } + + weak var simulation: Simulation2D? { + didSet { + guard let sim = self.simulation else { return } + self.calculatedStrength = strength.calculated(for: sim) + self.calculatedTargetOnDirection = targetOnDirection.calculated(for: sim) + } + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + let lane = self.direction + for i in sim.nodePositions.indices { + sim.nodeVelocities[i][lane] += + (self.calculatedTargetOnDirection[i] - sim.nodePositions[i][lane]) + * self.calculatedStrength[i] * alpha + } + } + } + + /// Create a direction force that moves nodes to a target position. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + @discardableResult + public func createPositionForce( + direction: DirectionForce2D.Direction, + targetOnDirection: DirectionForce2D.TargetOnDirection, + strength: DirectionForce2D.Strength = .constant(1.0) + ) -> DirectionForce2D { + let force = DirectionForce2D( + direction: direction, + targetOnDirection: targetOnDirection, + strength: strength + ) + force.simulation = self + self.forces.append(force) + return force + } +} + +extension Simulation2D.DirectionForce2D.Strength { + public func calculated(for simulation: Simulation2D) -> [Double] { + switch self { + case .constant(let value): + return Array(repeating: value, count: simulation.nodeIds.count) + case .varied(let getter): + return simulation.nodeIds.map(getter) + } + } +} + +extension Simulation2D.DirectionForce2D.TargetOnDirection { + public func calculated(for simulation: Simulation2D) -> [Double] { + switch self { + case .constant(let value): + return Array(repeating: value, count: simulation.nodeIds.count) + case .varied(let getter): + return simulation.nodeIds.map(getter) + } + } +} + +extension Simulation2D.DirectionForce2D.Direction { + @inlinable var lane: Int { + switch self { + case .x: return 0 + case .y: return 1 + case .entryOfVector(let i): return i + } + } +} +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation2D/Simulation2D.LinkForce.swift b/Sources/ForceSimulation/Simulation2D/Simulation2D.LinkForce.swift new file mode 100644 index 0000000..61a73f3 --- /dev/null +++ b/Sources/ForceSimulation/Simulation2D/Simulation2D.LinkForce.swift @@ -0,0 +1,244 @@ +// +// LinkForce.swift +// +// +// Created by li3zhen1 on 10/16/23. +// + +#if canImport(simd) +import NDTree +import simd + +enum LinkForce2DError: Error { + case useBeforeSimulationInitialized +} +extension Simulation2D { + /// A force that represents links between nodes. + /// The complexity is `O(e)`, where `e` is the number of links. + /// See [Link Force - D3](https://d3js.org/d3-force/link). + final public class LinkForce: ForceLike + where NodeID: Hashable { + + public typealias V = simd_double2 + + /// + public enum LinkStiffness { + case constant(V.Scalar) + case varied((EdgeID, LinkLookup) -> V.Scalar) + case weightedByDegree(k: (EdgeID, LinkLookup) -> V.Scalar) + } + var linkStiffness: LinkStiffness + var calculatedStiffness: [V.Scalar] = [] + + /// + public typealias LengthScalar = V.Scalar + public enum LinkLength { + case constant(LengthScalar) + case varied((EdgeID, LinkLookup) -> LengthScalar) + } + var linkLength: LinkLength + var calculatedLength: [LengthScalar] = [] + + /// Bias + var calculatedBias: [V.Scalar] = [] + + /// Binding to simulation + /// + weak var simulation: Simulation2D? { + didSet { + + guard let sim = simulation else { return } + + linksOfIndices = links.map { l in + EdgeID( + sim.nodeIdToIndexLookup[l.source, default: 0], + sim.nodeIdToIndexLookup[l.target, default: 0] + ) + } + + self.lookup = .buildFromLinks(linksOfIndices) + + self.calculatedBias = linksOfIndices.map { l in + V.Scalar(lookup.count[l.source, default: 0]) + / V.Scalar( + lookup.count[l.target, default: 0] + lookup.count[l.source, default: 0]) + } + + let lookupWithOriginalID = LinkLookup.buildFromLinks(links) + self.calculatedLength = linkLength.calculated( + for: self.links, connectionLookupTable: lookupWithOriginalID) + self.calculatedStiffness = linkStiffness.calculated( + for: self.links, connectionLookupTable: lookupWithOriginalID) + } + } + + var iterationsPerTick: UInt + + internal var linksOfIndices: [EdgeID] = [] + var links: [EdgeID] + + public struct LinkLookup<_NodeID> where _NodeID: Hashable { + let sources: [_NodeID: [_NodeID]] + let targets: [_NodeID: [_NodeID]] + let count: [_NodeID: Int] + } + private var lookup = LinkLookup(sources: [:], targets: [:], count: [:]) + + internal init( + _ links: [EdgeID], + stiffness: LinkStiffness, + originalLength: LinkLength = .constant(30), + iterationsPerTick: UInt = 1 + ) { + self.links = links + self.iterationsPerTick = iterationsPerTick + self.linkStiffness = stiffness + self.linkLength = originalLength + + } + + public func apply() { + guard let sim = self.simulation else { return } + + let alpha = sim.alpha + + for _ in 0..], + stiffness: LinkForce.LinkStiffness = .weightedByDegree { _, _ in 1.0 }, + originalLength: LinkForce.LinkLength = .constant(30.0), + iterationsPerTick: UInt = 1 + ) -> LinkForce { + let linkForce = LinkForce( + links, stiffness: stiffness, originalLength: originalLength) + linkForce.simulation = self + self.forces.append(linkForce) + return linkForce + } + + /// Create a link force that represents links between nodes. It works like + /// there is a spring between each pair of nodes. + /// The complexity is `O(e)`, where `e` is the number of links. + /// See [Link Force - D3](https://d3js.org/d3-force/link). + /// - Parameters: + /// - links: The links between nodes. + /// - stiffness: The stiffness of the spring (or links). + /// - originalLength: The original length of the spring (or links). + @discardableResult + public func createLinkForce( + _ linkTuples: [(NodeID, NodeID)], + stiffness: LinkForce.LinkStiffness = .weightedByDegree { _, _ in 1.0 }, + originalLength: LinkForce.LinkLength = .constant(30.0), + iterationsPerTick: UInt = 1 + ) -> LinkForce { + let links = linkTuples.map { EdgeID($0.0, $0.1) } + let linkForce = LinkForce( + links, stiffness: stiffness, originalLength: originalLength) + linkForce.simulation = self + self.forces.append(linkForce) + return linkForce + } +} + +extension Simulation2D.LinkForce.LinkLookup { + static func buildFromLinks(_ links: [EdgeID<_NodeID>]) -> Self { + var sources: [_NodeID: [_NodeID]] = [:] + var targets: [_NodeID: [_NodeID]] = [:] + var count: [_NodeID: Int] = [:] + for link in links { + sources[link.source, default: []].append(link.target) + targets[link.target, default: []].append(link.source) + count[link.source, default: 0] += 1 + count[link.target, default: 0] += 1 + } + return Self(sources: sources, targets: targets, count: count) + } +} + +extension Simulation2D.LinkForce.LinkLength { + func calculated( + for links: [EdgeID], + connectionLookupTable: Simulation2D.LinkForce.LinkLookup + ) -> [Double] { + switch self { + case .constant(let value): + return links.map { _ in value } + case .varied(let f): + return links.map { link in + f(link, connectionLookupTable) + } + } + } +} + +extension Simulation2D.LinkForce.LinkStiffness { + func calculated( + for links: [EdgeID], + connectionLookupTable lookup: Simulation2D.LinkForce.LinkLookup + ) -> [Double] { + switch self { + case .constant(let value): + return links.map { _ in value } + case .varied(let f): + return links.map { link in + f(link, lookup) + } + case .weightedByDegree(let k): + return links.map { link in + k(link, lookup) + / Double( + min( + lookup.count[link.source, default: 0], + lookup.count[link.target, default: 0] + ) + ) + } + } + } +} +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation2D/Simulation2D.ManyBodyForce.swift b/Sources/ForceSimulation/Simulation2D/Simulation2D.ManyBodyForce.swift new file mode 100644 index 0000000..af24a99 --- /dev/null +++ b/Sources/ForceSimulation/Simulation2D/Simulation2D.ManyBodyForce.swift @@ -0,0 +1,339 @@ +// +// File.swift +// +// +// Created by li3zhen1 on 10/1/23. +// + + +#if canImport(simd) + +import NDTree +import simd + +enum ManyBodyForce2DError: Error { + case buildQuadTreeBeforeSimulationInitialized +} + +struct MassQuadtreeDelegate2D: QuadtreeDelegate where NodeID: Hashable { + + public var accumulatedMass: Double = .zero + public var accumulatedCount: Int = 0 + public var accumulatedMassWeightedPositions: simd_double2 = .zero + + @usableFromInline let massProvider: (NodeID) -> Double + + init( + massProvider: @escaping (NodeID) -> Double + ) { + self.massProvider = massProvider + } + + internal init( + initialAccumulatedProperty: Double, + initialAccumulatedCount: Int, + initialWeightedAccumulatedNodePositions: simd_double2, + massProvider: @escaping (NodeID) -> Double + ) { + self.accumulatedMass = initialAccumulatedProperty + self.accumulatedCount = initialAccumulatedCount + self.accumulatedMassWeightedPositions = initialWeightedAccumulatedNodePositions + self.massProvider = massProvider + } + + @inlinable mutating func didAddNode(_ node: NodeID, at position: simd_double2) { + let p = massProvider(node) + #if DEBUG + assert(p > 0) + #endif + accumulatedCount += 1 + accumulatedMass += p + accumulatedMassWeightedPositions += position * p + } + + @inlinable mutating func didRemoveNode(_ node: NodeID, at position: simd_double2) { + let p = massProvider(node) + accumulatedCount -= 1 + accumulatedMass -= p + accumulatedMassWeightedPositions -= position * p + // TODO: parent removal? + } + + @inlinable func copy() -> Self { + return Self( + initialAccumulatedProperty: self.accumulatedMass, + initialAccumulatedCount: self.accumulatedCount, + initialWeightedAccumulatedNodePositions: self.accumulatedMassWeightedPositions, + massProvider: self.massProvider + ) + } + + @inlinable func spawn() -> Self { + return Self(massProvider: self.massProvider) + } + + // @inlinable var centroid: simd_double2? { + // guard accumulatedCount > 0 else { return nil } + // return self.accumulatedMassWeightedPositions/self.accumulatedMass + // } +} + +extension Simulation2D { + + /// A force that simulate the many-body force. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. The complexity might degrade to `O(n^2)` if the nodes are too close to each other. + /// See [Manybody Force - D3](https://d3js.org/d3-force/many-body). + final public class ManyBodyForce: ForceLike + where NodeID: Hashable { + + var strength: Double + + public enum NodeMass { + case constant(Double) + case varied((NodeID) -> Double) + } + var mass: NodeMass + var precalculatedMass: [Double] = [] + + weak var simulation: Simulation2D? { + didSet { + guard let sim = self.simulation else { return } + self.precalculatedMass = self.mass.calculated(for: sim) + self.forces = [simd_double2](repeating: .zero, count: sim.nodePositions.count) + } + } + + var theta2: Double + var theta: Double { + didSet { + theta2 = theta * theta + } + } + + var distanceMin2: Double = 1 + var distanceMax2: Double = Double.infinity + var distanceMin: Double = 1 + var distanceMax: Double = Double.infinity + + internal init( + strength: Double, + nodeMass: NodeMass = .constant(1.0), + theta: Double = 0.9 + ) { + self.strength = strength + self.mass = nodeMass + self.theta = theta + self.theta2 = theta * theta + } + + var forces: [simd_double2] = [] + public func apply() { + guard let simulation else { return } + + let alpha = simulation.alpha + + try! calculateForce(alpha: alpha) //else { return } + + for i in simulation.nodeVelocities.indices { + simulation.nodeVelocities[i] += self.forces[i] / precalculatedMass[i] + } + } + + // private func getCoveringBox() throws -> NDBox { + // guard let simulation else { throw ManyBodyForceError.buildQuadTreeBeforeSimulationInitialized } + // var _p0 = simulation.nodes[0].position + // var _p1 = simulation.nodes[0].position + // + // for p in simulation.nodes { + // for i in 0..= _p1[i] { + // _p1[i] = p.position[i] + 1 + // } + // } + // } + // return NDBox(_p0, _p1) + // + // } + + func calculateForce(alpha: Double) throws { + + guard let sim = self.simulation else { + throw ManyBodyForce2DError.buildQuadTreeBeforeSimulationInitialized + } + + let coveringBox = NDBox.cover(of: sim.nodePositions) //try! getCoveringBox() + + let tree = Quadtree(box: coveringBox, clusterDistance: 1e-5) { + + return switch self.mass { + case .constant(let m): + MassQuadtreeDelegate2D { _ in m } + case .varied(_): + MassQuadtreeDelegate2D { index in + self.precalculatedMass[index] + } + } + } + + for i in sim.nodePositions.indices { + tree.add(i, at: sim.nodePositions[i]) + + #if DEBUG + assert(tree.delegate.accumulatedCount == (i + 1)) + #endif + + } + + // var forces = [simd_double2](repeating: .zero, count: sim.nodePositions.count) + + for i in sim.nodePositions.indices { + var f = simd_double2.zero + tree.visit { t in + + // guard t.delegate.accumulatedCount > 0 else { return false } + // + // let centroid = t.delegate.accumulatedMassWeightedPositions / t.delegate.accumulatedMass + // let vec = centroid - sim.nodePositions[i] + // + // var distanceSquared = vec.jiggled().lengthSquared() + // + // /// too far away, omit + // guard distanceSquared < self.distanceMax2 else { return false } + // + // + // + // /// too close, enlarge distance + // if distanceSquared < self.distanceMin2 { + // distanceSquared = (self.distanceMin2 * distanceSquared).squareRoot() + // } + // + // + // if t.nodePosition != nil { + // + // /// filled leaf + // if !t.nodeIndices.contains(i) { + // let k: Double = self.strength * alpha * t.delegate.accumulatedMass / distanceSquared / (distanceSquared).squareRoot() + // forces[i] += vec * k + // } + // + // return false + // + // } + // else if t.children != nil { + // + // let boxWidth = (t.box.p1 - t.box.p1)[0] + // + // /// internal, guard in 180 guarantees we have nodes here + // if distanceSquared * self.theta2 > boxWidth * boxWidth { + // // far enough + // let k: Double = self.strength * alpha * t.delegate.accumulatedMass / distanceSquared / (distanceSquared).squareRoot() + // forces[i] += vec * k + // return false + // } + // else { + // return true + // } + // } + // else { + // // empty leaf + // return false + // } + + guard t.delegate.accumulatedCount > 0 else { return false } + let centroid = + t.delegate.accumulatedMassWeightedPositions / t.delegate.accumulatedMass + + let vec = centroid - sim.nodePositions[i] + let boxWidth = (t.box.p1 - t.box.p0)[0] + var distanceSquared = simd_length_squared(vec.jiggled()) + + let farEnough: Bool = (distanceSquared * self.theta2) > (boxWidth * boxWidth) + + // let distance = distanceSquared.squareRoot() + + if distanceSquared < self.distanceMin2 { + distanceSquared = (self.distanceMin2 * distanceSquared).squareRoot() + } + + if farEnough { + + guard distanceSquared < self.distanceMax2 else { return true } + + /// Workaround for "The compiler is unable to type-check this expression in reasonable time; try breaking up the expression into distinct sub-expressions" + let k: Double = + self.strength * alpha * t.delegate.accumulatedMass / distanceSquared // distanceSquared.squareRoot() + + f += vec * k + return false + + } else if t.children != nil { + return true + } + + if t.isFilledLeaf { + + // for j in t.nodeIndices { + // if j != i { + // let k: Double = + // self.strength * alpha * self.precalculatedMass[j] / distanceSquared / distanceSquared.squareRoot() + // f += vec * k + // } + // } + if t.nodeIndices.contains(i) { return false } + + let massAcc = t.delegate.accumulatedMass + // t.nodeIndices.contains(i) ? (t.delegate.accumulatedMass-self.precalculatedMass[i]) : (t.delegate.accumulatedMass) + let k: Double = self.strength * alpha * massAcc / distanceSquared // distanceSquared.squareRoot() + f += vec * k + return false + } else { + return true + } + } + forces[i] = f + } + // return forces + } + + } + + /// Create a many-body force that simulate the many-body force. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. The complexity might degrade to `O(n^2)` if the nodes are too close to each other. + /// The force mimics the gravity force or electrostatic force. + /// See [Manybody Force - D3](https://d3js.org/d3-force/many-body). + /// - Parameters: + /// - strength: The strength of the force. When the strength is positive, the nodes are attracted to each other like gravity force, otherwise, the nodes are repelled like electrostatic force. + /// - nodeMass: The mass of the nodes. The mass is used to calculate the force. The default value is 1.0. + /// - theta: Determines how approximate the calculation is. The default value is 0.9. The higher the value, the more approximate and fast the calculation is. + @discardableResult + public func createManyBodyForce( + strength: Double, + nodeMass: ManyBodyForce.NodeMass = .constant(1.0) + ) -> ManyBodyForce { + let manyBodyForce = ManyBodyForce( + strength: strength, nodeMass: nodeMass) + manyBodyForce.simulation = self + self.forces.append(manyBodyForce) + return manyBodyForce + } +} + +extension Simulation2D.ManyBodyForce.NodeMass { + public func calculated(for simulation: Simulation2D) -> [Double] { + switch self { + case .constant(let m): + return Array(repeating: m, count: simulation.nodePositions.count) + case .varied(let massGetter): + return simulation.nodeIds.map { n in + return massGetter(n) + } + } + } +} + +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation2D/Simulation2D.RadialForce.swift b/Sources/ForceSimulation/Simulation2D/Simulation2D.RadialForce.swift new file mode 100644 index 0000000..30fb28a --- /dev/null +++ b/Sources/ForceSimulation/Simulation2D/Simulation2D.RadialForce.swift @@ -0,0 +1,114 @@ +// +// RadialForce.swift +// +// +// Created by li3zhen1 on 10/1/23. +// + +#if canImport(simd) +import NDTree +import simd + +extension Simulation2D { + /// A force that applies a radial force to all nodes. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + final public class RadialForce: ForceLike + where NodeID: Hashable { + + public typealias V = simd_double2 + + weak var simulation: Simulation2D? { + didSet { + guard let sim = self.simulation else { return } + self.calculatedStrength = strength.calculated(for: sim) + self.calculatedRadius = radius.calculated(for: sim) + } + } + + public var center: V + + /// Radius accessor + public enum NodeRadius { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + public var radius: NodeRadius + private var calculatedRadius: [V.Scalar] = [] + + /// Strength accessor + public enum Strength { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + public var strength: Strength + private var calculatedStrength: [V.Scalar] = [] + + public init(center: V, radius: NodeRadius, strength: Strength) { + self.center = center + self.radius = radius + self.strength = strength + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + for i in sim.nodePositions.indices { + let nodeId = i + let deltaPosition = (sim.nodePositions[i] - self.center).jiggled() + let r = simd_length(deltaPosition) + let k = + (self.calculatedRadius[nodeId] + * self.calculatedStrength[nodeId] * alpha) / r + sim.nodeVelocities[i] += deltaPosition * k + } + } + + } + + /// Create a radial force that applies a radial force to all nodes. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + /// - Parameters: + /// - center: The center of the force. + /// - radius: The radius of the force. + /// - strength: The strength of the force. + @discardableResult + public func createRadialForce( + center: V = .zero, + radius: RadialForce.NodeRadius, + strength: RadialForce.Strength = .constant(0.1) + ) -> RadialForce { + let f = RadialForce(center: center, radius: radius, strength: strength) + f.simulation = self + self.forces.append(f) + return f + } + +} + +extension Simulation2D.RadialForce.Strength { + public func calculated(for simulation: Simulation2D) -> [Double] { + switch self { + case .constant(let s): + return simulation.nodeIds.map { _ in s } + case .varied(let s): + return simulation.nodeIds.map(s) + } + } +} + +extension Simulation2D.RadialForce.NodeRadius { + public func calculated(for simulation: Simulation2D) -> [Double] { + switch self { + case .constant(let r): + return simulation.nodeIds.map { _ in r } + case .varied(let r): + return simulation.nodeIds.map(r) + } + } +} + +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation2D/Simulation2D.swift b/Sources/ForceSimulation/Simulation2D/Simulation2D.swift new file mode 100644 index 0000000..e7b4543 --- /dev/null +++ b/Sources/ForceSimulation/Simulation2D/Simulation2D.swift @@ -0,0 +1,129 @@ + +#if canImport(simd) + +import simd + +/// A 2-Dimensional force simulation running on `double` and `simd_double3` types. +public final class Simulation2D +where NodeID: Hashable { + public typealias V = simd_double2 + /// The type of the vector used in the simulation. + /// Usually this is `Scalar` if you are on Apple platforms. + public typealias Scalar = V.Scalar + + public let initializedAlpha: Scalar + + public var alpha: Scalar + public var alphaMin: Scalar + public var alphaDecay: Scalar + public var alphaTarget: Scalar + + public var velocityDecay: Scalar + + public internal(set) var forces: [any ForceLike] = [] + + /// The position of points stored in simulation. + /// Ordered as the nodeIds you passed in when initializing simulation. + /// They are always updated. + public internal(set) var nodePositions: [V] + + /// The velocities of points stored in simulation. + /// Ordered as the nodeIds you passed in when initializing simulation. + /// They are always updated. + public internal(set) var nodeVelocities: [V] + + /// The fixed positions of points stored in simulation. + /// Ordered as the nodeIds you passed in when initializing simulation. + /// They are always updated. + public internal(set) var nodeFixations: [V?] + + public private(set) var nodeIds: [NodeID] + + @usableFromInline internal private(set) var nodeIdToIndexLookup: [NodeID: Int] = [:] + + /// Create a new simulation. + /// - Parameters: + /// - nodeIds: Hashable identifiers for the nodes. Force simulation calculate them by order once created. + /// - alpha: + /// - alphaMin: + /// - alphaDecay: The larger the value, the faster the simulation converges to the final result. + /// - alphaTarget: + /// - velocityDecay: + /// - getInitialPosition: The closure to set the initial position of the node. If not provided, the initial position is set to zero. + public init( + nodeIds: [NodeID], + alpha: Scalar = 1, + alphaMin: Scalar = 1e-3, + alphaDecay: Scalar = 2e-3, + alphaTarget: Scalar = 0.0, + velocityDecay: Scalar = 0.6, + + setInitialStatus getInitialPosition: ( + (NodeID) -> V + )? = nil + + ) { + + self.alpha = alpha + self.initializedAlpha = alpha // record and reload this when restarted + + self.alphaMin = alphaMin + self.alphaDecay = alphaDecay + self.alphaTarget = alphaTarget + + self.velocityDecay = velocityDecay + + if let getInitialPosition { + self.nodePositions = nodeIds.map(getInitialPosition) + } else { + self.nodePositions = Array(repeating: .zero, count: nodeIds.count) + } + + self.nodeVelocities = Array(repeating: .zero, count: nodeIds.count) + self.nodeFixations = Array(repeating: nil, count: nodeIds.count) + + self.nodeIdToIndexLookup.reserveCapacity(nodeIds.count) + for i in nodeIds.indices { + self.nodeIdToIndexLookup[nodeIds[i]] = i + } + self.nodeIds = nodeIds + + } + + /// Get the index in the nodeArray for `nodeId` + /// - **Complexity**: O(1) + public func getIndex(of nodeId: NodeID) -> Int { + return nodeIdToIndexLookup[nodeId]! + } + + /// Reset the alpha. The points will move faster as alpha gets larger. + public func resetAlpha(_ alpha: Scalar) { + self.alpha = alpha + } + + /// Run the simulation for a number of iterations. + /// Goes through all the forces created. + /// The forces will call `apply` then the positions and velocities will be modified. + /// - Parameter iterationCount: Default to 1. + public func tick(iterationCount: UInt = 1) { + for _ in 0..? + + internal init(center: V, strength: V.Scalar) { + self.center = center + self.strength = strength + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + + var meanPosition = V.zero + for n in sim.nodePositions { + meanPosition += n //.position + } + let delta = meanPosition * (self.strength / V.Scalar(sim.nodePositions.count)) + + for i in sim.nodePositions.indices { + sim.nodePositions[i] -= delta + } + } + + } + + /// Create a center force that drives nodes towards the center. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Collide Force - D3](https://d3js.org/d3-force/collide). + /// - Parameters: + /// - center: The center of the force. + /// - strength: The strength of the force. + @discardableResult + public func createCenterForce(center: V, strength: V.Scalar = 0.1) -> CenterForce { + let f = CenterForce(center: center, strength: strength) + f.simulation = self + self.forces.append(f) + return f + } + +} + +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation3D/Simulation3D.CollideForce.swift b/Sources/ForceSimulation/Simulation3D/Simulation3D.CollideForce.swift new file mode 100644 index 0000000..02a9b5f --- /dev/null +++ b/Sources/ForceSimulation/Simulation3D/Simulation3D.CollideForce.swift @@ -0,0 +1,222 @@ +// +// File.swift +// +// +// Created by li3zhen1 on 10/17/23. +// +import NDTree +import simd + + +#if canImport(simd) + +struct MaxRadiusTreeDelegate3D: OctreeDelegate where NodeID: Hashable { + + public typealias V = simd_float3 + + public var maxNodeRadius: V.Scalar + + @usableFromInline var radiusProvider: (NodeID) -> V.Scalar + + mutating func didAddNode(_ nodeId: NodeID, at position: V) { + let p = radiusProvider(nodeId) + maxNodeRadius = max(maxNodeRadius, p) + } + + mutating func didRemoveNode(_ nodeId: NodeID, at position: V) { + if radiusProvider(nodeId) >= maxNodeRadius { + // 🤯 for Collide force, set to 0 is fine + // Otherwise you need to traverse the delegate again + maxNodeRadius = 0 + } + } + + func copy() -> MaxRadiusTreeDelegate3D { + return Self(maxNodeRadius: maxNodeRadius, radiusProvider: radiusProvider) + } + + func spawn() -> MaxRadiusTreeDelegate3D { + return Self(radiusProvider: radiusProvider) + } + + init(maxNodeRadius: V.Scalar = 0, radiusProvider: @escaping (NodeID) -> V.Scalar) { + self.maxNodeRadius = maxNodeRadius + self.radiusProvider = radiusProvider + } + +} + +extension Simulation3D { + + /// A force that prevents nodes from overlapping. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. + /// See [Collide Force - D3](https://d3js.org/d3-force/collide). + public final class CollideForce: ForceLike + where NodeID: Hashable { + + public typealias V = simd_float3 + + weak var simulation: Simulation3D? { + didSet { + guard let sim = simulation else { return } + self.calculatedRadius = radius.calculated(for: sim) + } + } + + public enum CollideRadius { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + public var radius: CollideRadius + var calculatedRadius: [V.Scalar] = [] + + public let iterationsPerTick: UInt + public var strength: V.Scalar + + internal init( + radius: CollideRadius, + strength: V.Scalar = 1.0, + iterationsPerTick: UInt = 1 + ) { + self.radius = radius + self.iterationsPerTick = iterationsPerTick + self.strength = strength + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + + for _ in 0..( + box: coveringBox, clusterDistance: clusterDistance + ) { + return switch self.radius { + case .constant(let m): + MaxRadiusTreeDelegate3D { _ in m } + case .varied(_): + MaxRadiusTreeDelegate3D { index in + self.calculatedRadius[index] + } + } + } + + for i in sim.nodePositions.indices { + tree.add(i, at: sim.nodePositions[i]) + } + + for i in sim.nodePositions.indices { + let iOriginalPosition = sim.nodePositions[i] + let iOriginalVelocity = sim.nodeVelocities[i] + let iR = self.calculatedRadius[i] + let iR2 = iR * iR + let iPosition = iOriginalPosition + iOriginalVelocity + + tree.visit { t in + + let maxRadiusOfQuad = t.delegate.maxNodeRadius + let deltaR = maxRadiusOfQuad + iR + + if t.nodePosition != nil { + for j in t.nodeIndices { + // print("\(i)<=>\(j)") + // is leaf, make sure every collision happens once. + guard j > i else { continue } + + let jR = self.calculatedRadius[j] + let jOriginalPosition = sim.nodePositions[j] + let jOriginalVelocity = sim.nodeVelocities[j] + var deltaPosition = + iPosition - (jOriginalPosition + jOriginalVelocity) + let l = simd_length_squared(deltaPosition) + + let deltaR = iR + jR + if l < deltaR * deltaR { + + var l = simd_length(deltaPosition.jiggled()) + l = (deltaR - l) / l * self.strength + + let jR2 = jR * jR + + let k = jR2 / (iR2 + jR2) + + deltaPosition *= l + + sim.nodeVelocities[i] += deltaPosition * k + sim.nodeVelocities[j] -= deltaPosition * (1 - k) + } + } + return false + } + + for laneIndex in t.box.p0.indices { + let _v = t.box.p0[laneIndex] + if _v > iPosition[laneIndex] + deltaR /* True if no overlap */ { + return false + } + } + + for laneIndex in t.box.p1.indices { + let _v = t.box.p1[laneIndex] + if _v < iPosition[laneIndex] - deltaR /* True if no overlap */ { + return false + } + } + return true + + // return + // !(t.quad.x0 > iPosition.x + deltaR /* True if no overlap */ + // || t.quad.x1 < iPosition.x - deltaR + // || t.quad.y0 > iPosition.y + deltaR + // || t.quad.y1 < iPosition.y - deltaR) + } + } + } + } + + } + + /// Create a collide force that prevents nodes from overlapping. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. + /// See [Collide Force - D3](https://d3js.org/d3-force/collide). + /// - Parameters: + /// - radius: The radius of the force. + /// - strength: The strength of the force. + /// - iterationsPerTick: The number of iterations per tick. + @discardableResult + public func createCollideForce( + radius: CollideForce.CollideRadius = .constant(3.0), + strength: V.Scalar = 1.0, + iterationsPerTick: UInt = 1 + ) -> CollideForce { + let f = CollideForce( + radius: radius, + strength: strength, + iterationsPerTick: iterationsPerTick + ) + f.simulation = self + self.forces.append(f) + return f + } + +} + +extension Simulation3D.CollideForce.CollideRadius { + public func calculated(for simulation: Simulation3D) -> [Float] { + switch self { + case .constant(let r): + return Array(repeating: r, count: simulation.nodePositions.count) + case .varied(let radiusProvider): + return simulation.nodeIds.map { radiusProvider($0) } + } + } +} + +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation3D/Simulation3D.DirectionForce.swift b/Sources/ForceSimulation/Simulation3D/Simulation3D.DirectionForce.swift new file mode 100644 index 0000000..c20520f --- /dev/null +++ b/Sources/ForceSimulation/Simulation3D/Simulation3D.DirectionForce.swift @@ -0,0 +1,129 @@ +// +// PositionForce.swift +// +// +// Created by li3zhen1 on 10/1/23. +// + +#if canImport(simd) +import NDTree +import simd + + + +extension Simulation3D { + + /// A force that moves nodes to a target position. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + final public class DirectionForce3D: ForceLike + where NodeID: Hashable { + + public typealias V = simd_float3 + + public enum Direction { + case x + case y + case entryOfVector(Int) + } + public enum Strength { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + + public enum TargetOnDirection { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + + public var strength: Strength + public var direction: Int + public var calculatedStrength: [V.Scalar] = [] + public var targetOnDirection: TargetOnDirection + public var calculatedTargetOnDirection: [V.Scalar] = [] + + internal init( + direction: Direction, targetOnDirection: TargetOnDirection, + strength: Strength = .constant(1.0) + ) { + + self.strength = strength + self.direction = direction.lane + self.targetOnDirection = targetOnDirection + } + + weak var simulation: Simulation3D? { + didSet { + guard let sim = self.simulation else { return } + self.calculatedStrength = strength.calculated(for: sim) + self.calculatedTargetOnDirection = targetOnDirection.calculated(for: sim) + } + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + let lane = self.direction + for i in sim.nodePositions.indices { + sim.nodeVelocities[i][lane] += + (self.calculatedTargetOnDirection[i] - sim.nodePositions[i][lane]) + * self.calculatedStrength[i] * alpha + } + } + } + + /// Create a direction force that moves nodes to a target position. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + @discardableResult + public func createPositionForce( + direction: DirectionForce3D.Direction, + targetOnDirection: DirectionForce3D.TargetOnDirection, + strength: DirectionForce3D.Strength = .constant(1.0) + ) -> DirectionForce3D { + let force = DirectionForce3D( + direction: direction, + targetOnDirection: targetOnDirection, + strength: strength + ) + force.simulation = self + self.forces.append(force) + return force + } +} + +extension Simulation3D.DirectionForce3D.Strength { + public func calculated(for simulation: Simulation3D) -> [Float] { + switch self { + case .constant(let value): + return Array(repeating: value, count: simulation.nodeIds.count) + case .varied(let getter): + return simulation.nodeIds.map(getter) + } + } +} + +extension Simulation3D.DirectionForce3D.TargetOnDirection { + public func calculated(for simulation: Simulation3D) -> [Float] { + switch self { + case .constant(let value): + return Array(repeating: value, count: simulation.nodeIds.count) + case .varied(let getter): + return simulation.nodeIds.map(getter) + } + } +} + +extension Simulation3D.DirectionForce3D.Direction { + @inlinable var lane: Int { + switch self { + case .x: return 0 + case .y: return 1 + case .entryOfVector(let i): return i + } + } +} + +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation3D/Simulation3D.LinkForce.swift b/Sources/ForceSimulation/Simulation3D/Simulation3D.LinkForce.swift new file mode 100644 index 0000000..a3d67fc --- /dev/null +++ b/Sources/ForceSimulation/Simulation3D/Simulation3D.LinkForce.swift @@ -0,0 +1,244 @@ +// +// LinkForce.swift +// +// +// Created by li3zhen1 on 10/16/23. +// + +#if canImport(simd) +import NDTree +import simd + +enum LinkForce3DError: Error { + case useBeforeSimulationInitialized +} +extension Simulation3D { + /// A force that represents links between nodes. + /// The complexity is `O(e)`, where `e` is the number of links. + /// See [Link Force - D3](https://d3js.org/d3-force/link). + final public class LinkForce: ForceLike + where NodeID: Hashable { + + public typealias V = simd_float3 + + /// + public enum LinkStiffness { + case constant(V.Scalar) + case varied((EdgeID, LinkLookup) -> V.Scalar) + case weightedByDegree(k: (EdgeID, LinkLookup) -> V.Scalar) + } + var linkStiffness: LinkStiffness + var calculatedStiffness: [V.Scalar] = [] + + /// + public typealias LengthScalar = V.Scalar + public enum LinkLength { + case constant(LengthScalar) + case varied((EdgeID, LinkLookup) -> LengthScalar) + } + var linkLength: LinkLength + var calculatedLength: [LengthScalar] = [] + + /// Bias + var calculatedBias: [V.Scalar] = [] + + /// Binding to simulation + /// + weak var simulation: Simulation3D? { + didSet { + + guard let sim = simulation else { return } + + linksOfIndices = links.map { l in + EdgeID( + sim.nodeIdToIndexLookup[l.source, default: 0], + sim.nodeIdToIndexLookup[l.target, default: 0] + ) + } + + self.lookup = .buildFromLinks(linksOfIndices) + + self.calculatedBias = linksOfIndices.map { l in + V.Scalar(lookup.count[l.source, default: 0]) + / V.Scalar( + lookup.count[l.target, default: 0] + lookup.count[l.source, default: 0]) + } + + let lookupWithOriginalID = LinkLookup.buildFromLinks(links) + self.calculatedLength = linkLength.calculated( + for: self.links, connectionLookupTable: lookupWithOriginalID) + self.calculatedStiffness = linkStiffness.calculated( + for: self.links, connectionLookupTable: lookupWithOriginalID) + } + } + + var iterationsPerTick: UInt + + internal var linksOfIndices: [EdgeID] = [] + var links: [EdgeID] + + public struct LinkLookup<_NodeID> where _NodeID: Hashable { + let sources: [_NodeID: [_NodeID]] + let targets: [_NodeID: [_NodeID]] + let count: [_NodeID: Int] + } + private var lookup = LinkLookup(sources: [:], targets: [:], count: [:]) + + internal init( + _ links: [EdgeID], + stiffness: LinkStiffness, + originalLength: LinkLength = .constant(30), + iterationsPerTick: UInt = 1 + ) { + self.links = links + self.iterationsPerTick = iterationsPerTick + self.linkStiffness = stiffness + self.linkLength = originalLength + + } + + public func apply() { + guard let sim = self.simulation else { return } + + let alpha = sim.alpha + + for _ in 0..], + stiffness: LinkForce.LinkStiffness = .weightedByDegree { _, _ in 1.0 }, + originalLength: LinkForce.LinkLength = .constant(30.0), + iterationsPerTick: UInt = 1 + ) -> LinkForce { + let linkForce = LinkForce( + links, stiffness: stiffness, originalLength: originalLength) + linkForce.simulation = self + self.forces.append(linkForce) + return linkForce + } + + /// Create a link force that represents links between nodes. It works like + /// there is a spring between each pair of nodes. + /// The complexity is `O(e)`, where `e` is the number of links. + /// See [Link Force - D3](https://d3js.org/d3-force/link). + /// - Parameters: + /// - links: The links between nodes. + /// - stiffness: The stiffness of the spring (or links). + /// - originalLength: The original length of the spring (or links). + @discardableResult + public func createLinkForce( + _ linkTuples: [(NodeID, NodeID)], + stiffness: LinkForce.LinkStiffness = .weightedByDegree { _, _ in 1.0 }, + originalLength: LinkForce.LinkLength = .constant(30.0), + iterationsPerTick: UInt = 1 + ) -> LinkForce { + let links = linkTuples.map { EdgeID($0.0, $0.1) } + let linkForce = LinkForce( + links, stiffness: stiffness, originalLength: originalLength) + linkForce.simulation = self + self.forces.append(linkForce) + return linkForce + } +} + +extension Simulation3D.LinkForce.LinkLookup { + static func buildFromLinks(_ links: [EdgeID<_NodeID>]) -> Self { + var sources: [_NodeID: [_NodeID]] = [:] + var targets: [_NodeID: [_NodeID]] = [:] + var count: [_NodeID: Int] = [:] + for link in links { + sources[link.source, default: []].append(link.target) + targets[link.target, default: []].append(link.source) + count[link.source, default: 0] += 1 + count[link.target, default: 0] += 1 + } + return Self(sources: sources, targets: targets, count: count) + } +} + +extension Simulation3D.LinkForce.LinkLength { + func calculated( + for links: [EdgeID], + connectionLookupTable: Simulation3D.LinkForce.LinkLookup + ) -> [Float] { + switch self { + case .constant(let value): + return links.map { _ in value } + case .varied(let f): + return links.map { link in + f(link, connectionLookupTable) + } + } + } +} + +extension Simulation3D.LinkForce.LinkStiffness { + func calculated( + for links: [EdgeID], + connectionLookupTable lookup: Simulation3D.LinkForce.LinkLookup + ) -> [Float] { + switch self { + case .constant(let value): + return links.map { _ in value } + case .varied(let f): + return links.map { link in + f(link, lookup) + } + case .weightedByDegree(let k): + return links.map { link in + k(link, lookup) + / Float( + min( + lookup.count[link.source, default: 0], + lookup.count[link.target, default: 0] + ) + ) + } + } + } +} +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation3D/Simulation3D.ManyBodyForce.swift b/Sources/ForceSimulation/Simulation3D/Simulation3D.ManyBodyForce.swift new file mode 100644 index 0000000..dcb2960 --- /dev/null +++ b/Sources/ForceSimulation/Simulation3D/Simulation3D.ManyBodyForce.swift @@ -0,0 +1,338 @@ +// +// File.swift +// +// +// Created by li3zhen1 on 10/1/23. +// + +#if canImport(simd) +import NDTree +import simd + +enum ManyBodyForce3DError: Error { + case buildQuadTreeBeforeSimulationInitialized +} + +struct MassQuadtreeDelegate3D: OctreeDelegate where NodeID: Hashable { + + public typealias V = simd_float3 + + public var accumulatedMass: Float = .zero + public var accumulatedCount: Int = 0 + public var accumulatedMassWeightedPositions: V = .zero + + @usableFromInline let massProvider: (NodeID) -> Float + + init( + massProvider: @escaping (NodeID) -> Float + ) { + self.massProvider = massProvider + } + + internal init( + initialAccumulatedProperty: Float, + initialAccumulatedCount: Int, + initialWeightedAccumulatedNodePositions: V, + massProvider: @escaping (NodeID) -> Float + ) { + self.accumulatedMass = initialAccumulatedProperty + self.accumulatedCount = initialAccumulatedCount + self.accumulatedMassWeightedPositions = initialWeightedAccumulatedNodePositions + self.massProvider = massProvider + } + + @inlinable mutating func didAddNode(_ node: NodeID, at position: V) { + let p = massProvider(node) + #if DEBUG + assert(p > 0) + #endif + accumulatedCount += 1 + accumulatedMass += p + accumulatedMassWeightedPositions += position * p + } + + @inlinable mutating func didRemoveNode(_ node: NodeID, at position: V) { + let p = massProvider(node) + accumulatedCount -= 1 + accumulatedMass -= p + accumulatedMassWeightedPositions -= position * p + // TODO: parent removal? + } + + @inlinable func copy() -> Self { + return Self( + initialAccumulatedProperty: self.accumulatedMass, + initialAccumulatedCount: self.accumulatedCount, + initialWeightedAccumulatedNodePositions: self.accumulatedMassWeightedPositions, + massProvider: self.massProvider + ) + } + + @inlinable func spawn() -> Self { + return Self(massProvider: self.massProvider) + } + + // @inlinable var centroid: V? { + // guard accumulatedCount > 0 else { return nil } + // return self.accumulatedMassWeightedPositions/self.accumulatedMass + // } +} + +extension Simulation3D { + + /// A force that simulate the many-body force. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. The complexity might degrade to `O(n^2)` if the nodes are too close to each other. + /// See [Manybody Force - D3](https://d3js.org/d3-force/many-body). + final public class ManyBodyForce: ForceLike + where NodeID: Hashable { + + var strength: Float + + public enum NodeMass { + case constant(Float) + case varied((NodeID) -> Float) + } + var mass: NodeMass + var precalculatedMass: [Float] = [] + + weak var simulation: Simulation3D? { + didSet { + guard let sim = self.simulation else { return } + self.precalculatedMass = self.mass.calculated(for: sim) + self.forces = [V](repeating: .zero, count: sim.nodePositions.count) + } + } + + var theta2: Float + var theta: Float { + didSet { + theta2 = theta * theta + } + } + + var distanceMin2: Float = 1 + var distanceMax2: Float = Float.infinity + var distanceMin: Float = 1 + var distanceMax: Float = Float.infinity + + internal init( + strength: Float, + nodeMass: NodeMass = .constant(1.0), + theta: Float = 0.9 + ) { + self.strength = strength + self.mass = nodeMass + self.theta = theta + self.theta2 = theta * theta + } + + var forces: [V] = [] + public func apply() { + guard let simulation else { return } + + let alpha = simulation.alpha + + try! calculateForce(alpha: alpha) //else { return } + + for i in simulation.nodeVelocities.indices { + simulation.nodeVelocities[i] += self.forces[i] / precalculatedMass[i] + } + } + + // private func getCoveringBox() throws -> NDBox { + // guard let simulation else { throw ManyBodyForceError.buildQuadTreeBeforeSimulationInitialized } + // var _p0 = simulation.nodes[0].position + // var _p1 = simulation.nodes[0].position + // + // for p in simulation.nodes { + // for i in 0..= _p1[i] { + // _p1[i] = p.position[i] + 1 + // } + // } + // } + // return NDBox(_p0, _p1) + // + // } + + func calculateForce(alpha: Float) throws { + + guard let sim = self.simulation else { + throw ManyBodyForce3DError.buildQuadTreeBeforeSimulationInitialized + } + + let coveringBox = OctBox.cover(of: sim.nodePositions) //try! getCoveringBox() + + let tree = Octree(box: coveringBox, clusterDistance: 1e-5) { + + return switch self.mass { + case .constant(let m): + MassQuadtreeDelegate3D { _ in m } + case .varied(_): + MassQuadtreeDelegate3D { index in + self.precalculatedMass[index] + } + } + } + + for i in sim.nodePositions.indices { + tree.add(i, at: sim.nodePositions[i]) + + #if DEBUG + assert(tree.delegate.accumulatedCount == (i + 1)) + #endif + + } + + // var forces = [V](repeating: .zero, count: sim.nodePositions.count) + + for i in sim.nodePositions.indices { + var f = V.zero + tree.visit { t in + + // guard t.delegate.accumulatedCount > 0 else { return false } + // + // let centroid = t.delegate.accumulatedMassWeightedPositions / t.delegate.accumulatedMass + // let vec = centroid - sim.nodePositions[i] + // + // var distanceSquared = vec.jiggled().lengthSquared() + // + // /// too far away, omit + // guard distanceSquared < self.distanceMax2 else { return false } + // + // + // + // /// too close, enlarge distance + // if distanceSquared < self.distanceMin2 { + // distanceSquared = (self.distanceMin2 * distanceSquared).squareRoot() + // } + // + // + // if t.nodePosition != nil { + // + // /// filled leaf + // if !t.nodeIndices.contains(i) { + // let k: Float = self.strength * alpha * t.delegate.accumulatedMass / distanceSquared / (distanceSquared).squareRoot() + // forces[i] += vec * k + // } + // + // return false + // + // } + // else if t.children != nil { + // + // let boxWidth = (t.box.p1 - t.box.p1)[0] + // + // /// internal, guard in 180 guarantees we have nodes here + // if distanceSquared * self.theta2 > boxWidth * boxWidth { + // // far enough + // let k: Float = self.strength * alpha * t.delegate.accumulatedMass / distanceSquared / (distanceSquared).squareRoot() + // forces[i] += vec * k + // return false + // } + // else { + // return true + // } + // } + // else { + // // empty leaf + // return false + // } + + guard t.delegate.accumulatedCount > 0 else { return false } + let centroid = + t.delegate.accumulatedMassWeightedPositions / t.delegate.accumulatedMass + + let vec = centroid - sim.nodePositions[i] + let boxWidth = (t.box.p1 - t.box.p0)[0] + var distanceSquared = simd_length_squared(vec.jiggled()) + + let farEnough: Bool = (distanceSquared * self.theta2) > (boxWidth * boxWidth) + + // let distance = distanceSquared.squareRoot() + + if distanceSquared < self.distanceMin2 { + distanceSquared = (self.distanceMin2 * distanceSquared).squareRoot() + } + + if farEnough { + + guard distanceSquared < self.distanceMax2 else { return true } + + /// Workaround for "The compiler is unable to type-check this expression in reasonable time; try breaking up the expression into distinct sub-expressions" + let k: Float = + self.strength * alpha * t.delegate.accumulatedMass / distanceSquared // distanceSquared.squareRoot() + + f += vec * k + return false + + } else if t.children != nil { + return true + } + + if t.isFilledLeaf { + + // for j in t.nodeIndices { + // if j != i { + // let k: Float = + // self.strength * alpha * self.precalculatedMass[j] / distanceSquared / distanceSquared.squareRoot() + // f += vec * k + // } + // } + if t.nodeIndices.contains(i) { return false } + + let massAcc = t.delegate.accumulatedMass + // t.nodeIndices.contains(i) ? (t.delegate.accumulatedMass-self.precalculatedMass[i]) : (t.delegate.accumulatedMass) + let k: Float = self.strength * alpha * massAcc / distanceSquared // distanceSquared.squareRoot() + f += vec * k + return false + } else { + return true + } + } + forces[i] = f + } + // return forces + } + + } + + /// Create a many-body force that simulate the many-body force. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. The complexity might degrade to `O(n^2)` if the nodes are too close to each other. + /// The force mimics the gravity force or electrostatic force. + /// See [Manybody Force - D3](https://d3js.org/d3-force/many-body). + /// - Parameters: + /// - strength: The strength of the force. When the strength is positive, the nodes are attracted to each other like gravity force, otherwise, the nodes are repelled like electrostatic force. + /// - nodeMass: The mass of the nodes. The mass is used to calculate the force. The default value is 1.0. + /// - theta: Determines how approximate the calculation is. The default value is 0.9. The higher the value, the more approximate and fast the calculation is. + @discardableResult + public func createManyBodyForce( + strength: Float, + nodeMass: ManyBodyForce.NodeMass = .constant(1.0) + ) -> ManyBodyForce { + let manyBodyForce = ManyBodyForce( + strength: strength, nodeMass: nodeMass) + manyBodyForce.simulation = self + self.forces.append(manyBodyForce) + return manyBodyForce + } +} + +extension Simulation3D.ManyBodyForce.NodeMass { + public func calculated(for simulation: Simulation3D) -> [Float] { + switch self { + case .constant(let m): + return Array(repeating: m, count: simulation.nodePositions.count) + case .varied(let massGetter): + return simulation.nodeIds.map { n in + return massGetter(n) + } + } + } +} +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation3D/Simulation3D.RadialForce.swift b/Sources/ForceSimulation/Simulation3D/Simulation3D.RadialForce.swift new file mode 100644 index 0000000..0bf3449 --- /dev/null +++ b/Sources/ForceSimulation/Simulation3D/Simulation3D.RadialForce.swift @@ -0,0 +1,113 @@ +// +// RadialForce.swift +// +// +// Created by li3zhen1 on 10/1/23. +// + +#if canImport(simd) +import NDTree +import simd + +extension Simulation3D { + /// A force that applies a radial force to all nodes. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + final public class RadialForce: ForceLike + where NodeID: Hashable { + + weak var simulation: Simulation3D? { + didSet { + guard let sim = self.simulation else { return } + self.calculatedStrength = strength.calculated(for: sim) + self.calculatedRadius = radius.calculated(for: sim) + } + } + + public var center: V + + /// Radius accessor + public enum NodeRadius { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + public var radius: NodeRadius + private var calculatedRadius: [V.Scalar] = [] + + /// Strength accessor + public enum Strength { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + public var strength: Strength + private var calculatedStrength: [V.Scalar] = [] + + public init(center: V, radius: NodeRadius, strength: Strength) { + self.center = center + self.radius = radius + self.strength = strength + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + for i in sim.nodePositions.indices { + let nodeId = i + + let deltaPosition = (sim.nodePositions[i] - self.center).jiggled() + let r = simd_length(deltaPosition) + let k = + (self.calculatedRadius[nodeId] + * self.calculatedStrength[nodeId] * alpha) / r + sim.nodeVelocities[i] += deltaPosition * k + } + } + + } + + /// Create a radial force that applies a radial force to all nodes. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + /// - Parameters: + /// - center: The center of the force. + /// - radius: The radius of the force. + /// - strength: The strength of the force. + @discardableResult + public func createRadialForce( + center: V = .zero, + radius: RadialForce.NodeRadius, + strength: RadialForce.Strength = .constant(0.1) + ) -> RadialForce { + let f = RadialForce(center: center, radius: radius, strength: strength) + f.simulation = self + self.forces.append(f) + return f + } + +} + +extension Simulation3D.RadialForce.Strength { + public func calculated(for simulation: Simulation3D) -> [Float] { + switch self { + case .constant(let s): + return simulation.nodeIds.map { _ in s } + case .varied(let s): + return simulation.nodeIds.map(s) + } + } +} + +extension Simulation3D.RadialForce.NodeRadius { + public func calculated(for simulation: Simulation3D) -> [Float] { + switch self { + case .constant(let r): + return simulation.nodeIds.map { _ in r } + case .varied(let r): + return simulation.nodeIds.map(r) + } + } +} + +#endif \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation3D/Simulation3D.swift b/Sources/ForceSimulation/Simulation3D/Simulation3D.swift new file mode 100644 index 0000000..728bf75 --- /dev/null +++ b/Sources/ForceSimulation/Simulation3D/Simulation3D.swift @@ -0,0 +1,142 @@ +// +// Simulation.swift +// +// +// Created by li3zhen1 on 10/16/23. +// + +#if canImport(simd) + +import NDTree +import simd + +enum Simulation3DError: Error { + case subscriptionToNonexistentNode +} + +/// A 3-Dimensional force simulation running on `float` and `simd_float3` types. +public final class Simulation3D +where NodeID: Hashable{ + + public typealias V = simd_float3 + + /// The type of the vector used in the simulation. + /// Usually this is `Scalar` if you are on Apple platforms. + public typealias Scalar = V.Scalar + + public let initializedAlpha: Scalar + + public var alpha: Scalar + public var alphaMin: Scalar + public var alphaDecay: Scalar + public var alphaTarget: Scalar + + public var velocityDecay: Scalar + + public internal(set) var forces: [any ForceLike] = [] + + /// The position of points stored in simulation. + /// Ordered as the nodeIds you passed in when initializing simulation. + /// They are always updated. + public internal(set) var nodePositions: [V] + + /// The velocities of points stored in simulation. + /// Ordered as the nodeIds you passed in when initializing simulation. + /// They are always updated. + public internal(set) var nodeVelocities: [V] + + /// The fixed positions of points stored in simulation. + /// Ordered as the nodeIds you passed in when initializing simulation. + /// They are always updated. + public internal(set) var nodeFixations: [V?] + + public private(set) var nodeIds: [NodeID] + + @usableFromInline internal private(set) var nodeIdToIndexLookup: [NodeID: Int] = [:] + + /// Create a new simulation. + /// - Parameters: + /// - nodeIds: Hashable identifiers for the nodes. Force simulation calculate them by order once created. + /// - alpha: + /// - alphaMin: + /// - alphaDecay: The larger the value, the faster the simulation converges to the final result. + /// - alphaTarget: + /// - velocityDecay: + /// - getInitialPosition: The closure to set the initial position of the node. If not provided, the initial position is set to zero. + public init( + nodeIds: [NodeID], + alpha: Scalar = 1, + alphaMin: Scalar = 1e-3, + alphaDecay: Scalar = 2e-3, + alphaTarget: Scalar = 0.0, + velocityDecay: Scalar = 0.6, + + setInitialStatus getInitialPosition: ( + (NodeID) -> V + )? = nil + + ) { + + self.alpha = alpha + self.initializedAlpha = alpha // record and reload this when restarted + + self.alphaMin = alphaMin + self.alphaDecay = alphaDecay + self.alphaTarget = alphaTarget + + self.velocityDecay = velocityDecay + + if let getInitialPosition { + self.nodePositions = nodeIds.map(getInitialPosition) + } else { + self.nodePositions = Array(repeating: .zero, count: nodeIds.count) + } + + self.nodeVelocities = Array(repeating: .zero, count: nodeIds.count) + self.nodeFixations = Array(repeating: nil, count: nodeIds.count) + + self.nodeIdToIndexLookup.reserveCapacity(nodeIds.count) + for i in nodeIds.indices { + self.nodeIdToIndexLookup[nodeIds[i]] = i + } + self.nodeIds = nodeIds + + } + + /// Get the index in the nodeArray for `nodeId` + /// - **Complexity**: O(1) + public func getIndex(of nodeId: NodeID) -> Int { + return nodeIdToIndexLookup[nodeId]! + } + + /// Reset the alpha. The points will move faster as alpha gets larger. + public func resetAlpha(_ alpha: Scalar) { + self.alpha = alpha + } + + /// Run the simulation for a number of iterations. + /// Goes through all the forces created. + /// The forces will call `apply` then the positions and velocities will be modified. + /// - Parameter iterationCount: Default to 1. + public func tick(iterationCount: UInt = 1) { + for _ in 0.. CenterForce { + let f = CenterForce(center: center, strength: strength) + f.simulation = self + self.forces.append(f) + return f + } + +} diff --git a/Sources/ForceSimulation/SimulationKD/SimulationKD.CollideForce.swift b/Sources/ForceSimulation/SimulationKD/SimulationKD.CollideForce.swift new file mode 100644 index 0000000..ab890dc --- /dev/null +++ b/Sources/ForceSimulation/SimulationKD/SimulationKD.CollideForce.swift @@ -0,0 +1,211 @@ +// +// File.swift +// +// +// Created by li3zhen1 on 10/17/23. +// +import NDTree + +struct MaxRadiusTreeDelegate: NDTreeDelegate where NodeID: Hashable, V: VectorLike { + + public var maxNodeRadius: V.Scalar + + @usableFromInline var radiusProvider: (NodeID) -> V.Scalar + + mutating func didAddNode(_ nodeId: NodeID, at position: V) { + let p = radiusProvider(nodeId) + maxNodeRadius = max(maxNodeRadius, p) + } + + mutating func didRemoveNode(_ nodeId: NodeID, at position: V) { + if radiusProvider(nodeId) >= maxNodeRadius { + // 🤯 for Collide force, set to 0 is fine + // Otherwise you need to traverse the delegate again + maxNodeRadius = 0 + } + } + + func copy() -> MaxRadiusTreeDelegate { + return Self(maxNodeRadius: maxNodeRadius, radiusProvider: radiusProvider) + } + + func spawn() -> MaxRadiusTreeDelegate { + return Self(radiusProvider: radiusProvider) + } + + init(maxNodeRadius: V.Scalar = 0, radiusProvider: @escaping (NodeID) -> V.Scalar) { + self.maxNodeRadius = maxNodeRadius + self.radiusProvider = radiusProvider + } + +} + +extension SimulationKD { + /// A force that prevents nodes from overlapping. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. + /// See [Collide Force - D3](https://d3js.org/d3-force/collide). + public final class CollideForce: ForceLike + where NodeID: Hashable, V: VectorLike, V.Scalar : SimulatableFloatingPoint { + + weak var simulation: SimulationKD? { + didSet { + guard let sim = simulation else { return } + self.calculatedRadius = radius.calculated(for: sim) + } + } + + public enum CollideRadius { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + public var radius: CollideRadius + var calculatedRadius: [V.Scalar] = [] + + public let iterationsPerTick: UInt + public var strength: V.Scalar + + internal init( + radius: CollideRadius, + strength: V.Scalar = 1.0, + iterationsPerTick: UInt = 1 + ) { + self.radius = radius + self.iterationsPerTick = iterationsPerTick + self.strength = strength + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + + for _ in 0...cover(of: sim.nodePositions) + + let clusterDistance: V.Scalar = V.Scalar(Int(0.00001)) + + let tree = NDTree>( + box: coveringBox, clusterDistance: clusterDistance + ) { + return switch self.radius { + case .constant(let m): + MaxRadiusTreeDelegate { _ in m } + case .varied(_): + MaxRadiusTreeDelegate { index in + self.calculatedRadius[index] + } + } + } + + for i in sim.nodePositions.indices { + tree.add(i, at: sim.nodePositions[i]) + } + + for i in sim.nodePositions.indices { + let iOriginalPosition = sim.nodePositions[i] + let iOriginalVelocity = sim.nodeVelocities[i] + let iR = self.calculatedRadius[i] + let iR2 = iR * iR + let iPosition = iOriginalPosition + iOriginalVelocity + + tree.visit { t in + + let maxRadiusOfQuad = t.delegate.maxNodeRadius + let deltaR = maxRadiusOfQuad + iR + + if t.nodePosition != nil { + for j in t.nodeIndices { + // print("\(i)<=>\(j)") + // is leaf, make sure every collision happens once. + guard j > i else { continue } + + let jR = self.calculatedRadius[j] + let jOriginalPosition = sim.nodePositions[j] + let jOriginalVelocity = sim.nodeVelocities[j] + var deltaPosition = + iPosition - (jOriginalPosition + jOriginalVelocity) + let l = deltaPosition.lengthSquared() + + let deltaR = iR + jR + if l < deltaR * deltaR { + + var l = deltaPosition.jiggled().length() + l = (deltaR - l) / l * self.strength + + let jR2 = jR * jR + + let k = jR2 / (iR2 + jR2) + + deltaPosition *= l + + sim.nodeVelocities[i] += deltaPosition * k + sim.nodeVelocities[j] -= deltaPosition * (1 - k) + } + } + return false + } + + for laneIndex in t.box.p0.indices { + let _v = t.box.p0[laneIndex] + if _v > iPosition[laneIndex] + deltaR /* True if no overlap */ { + return false + } + } + + for laneIndex in t.box.p1.indices { + let _v = t.box.p1[laneIndex] + if _v < iPosition[laneIndex] - deltaR /* True if no overlap */ { + return false + } + } + return true + + // return + // !(t.quad.x0 > iPosition.x + deltaR /* True if no overlap */ + // || t.quad.x1 < iPosition.x - deltaR + // || t.quad.y0 > iPosition.y + deltaR + // || t.quad.y1 < iPosition.y - deltaR) + } + } + } + } + + } + + /// Create a collide force that prevents nodes from overlapping. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. + /// See [Collide Force - D3](https://d3js.org/d3-force/collide). + /// - Parameters: + /// - radius: The radius of the force. + /// - strength: The strength of the force. + /// - iterationsPerTick: The number of iterations per tick. + @discardableResult + public func createCollideForce( + radius: CollideForce.CollideRadius = .constant(3.0), + strength: V.Scalar = 1.0, + iterationsPerTick: UInt = 1 + ) -> CollideForce { + let f = CollideForce( + radius: radius, + strength: strength, + iterationsPerTick: iterationsPerTick + ) + f.simulation = self + self.forces.append(f) + return f + } + +} + +extension SimulationKD.CollideForce.CollideRadius { + public func calculated(for simulation: SimulationKD) -> [V.Scalar] { + switch self { + case .constant(let r): + return Array(repeating: r, count: simulation.nodePositions.count) + case .varied(let radiusProvider): + return simulation.nodeIds.map { radiusProvider($0) } + } + } +} diff --git a/Sources/ForceSimulation/SimulationKD/SimulationKD.DirectionForce.swift b/Sources/ForceSimulation/SimulationKD/SimulationKD.DirectionForce.swift new file mode 100644 index 0000000..0462286 --- /dev/null +++ b/Sources/ForceSimulation/SimulationKD/SimulationKD.DirectionForce.swift @@ -0,0 +1,120 @@ +// +// PositionForce.swift +// +// +// Created by li3zhen1 on 10/1/23. +// + +import NDTree + +extension SimulationKD { + /// A force that moves nodes to a target position. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + final public class DirectionForce: ForceLike + where NodeID: Hashable, V: VectorLike, V.Scalar : SimulatableFloatingPoint { + + public enum Direction { + case x + case y + case entryOfVector(Int) + } + public enum Strength { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + + public enum TargetOnDirection { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + + public var strength: Strength + public var direction: Int + public var calculatedStrength: [V.Scalar] = [] + public var targetOnDirection: TargetOnDirection + public var calculatedTargetOnDirection: [V.Scalar] = [] + + internal init( + direction: Direction, targetOnDirection: TargetOnDirection, + strength: Strength = .constant(1.0) + ) { + + self.strength = strength + self.direction = direction.lane + self.targetOnDirection = targetOnDirection + } + + weak var simulation: SimulationKD? { + didSet { + guard let sim = self.simulation else { return } + self.calculatedStrength = strength.calculated(for: sim) + self.calculatedTargetOnDirection = targetOnDirection.calculated(for: sim) + } + } + + public func apply() { + guard let sim = self.simulation else { return } + let alpha = sim.alpha + let lane = self.direction + for i in sim.nodePositions.indices { + sim.nodeVelocities[i][lane] += + (self.calculatedTargetOnDirection[i] - sim.nodePositions[i][lane]) + * self.calculatedStrength[i] * alpha + } + } + } + + /// Create a direction force that moves nodes to a target position. + /// Center force is relatively fast, the complexity is `O(n)`, + /// where `n` is the number of nodes. + /// See [Position Force - D3](https://d3js.org/d3-force/position). + @discardableResult + public func createPositionForce( + direction: DirectionForce.Direction, + targetOnDirection: DirectionForce.TargetOnDirection, + strength: DirectionForce.Strength = .constant(1.0) + ) -> DirectionForce { + let force = DirectionForce( + direction: direction, + targetOnDirection: targetOnDirection, + strength: strength + ) + force.simulation = self + self.forces.append(force) + return force + } +} + +extension SimulationKD.DirectionForce.Strength { + public func calculated(for simulation: SimulationKD) -> [V.Scalar] { + switch self { + case .constant(let value): + return Array(repeating: value, count: simulation.nodeIds.count) + case .varied(let getter): + return simulation.nodeIds.map(getter) + } + } +} + +extension SimulationKD.DirectionForce.TargetOnDirection { + public func calculated(for simulation: SimulationKD) -> [V.Scalar] { + switch self { + case .constant(let value): + return Array(repeating: value, count: simulation.nodeIds.count) + case .varied(let getter): + return simulation.nodeIds.map(getter) + } + } +} + +extension SimulationKD.DirectionForce.Direction { + @inlinable var lane: Int { + switch self { + case .x: return 0 + case .y: return 1 + case .entryOfVector(let i): return i + } + } +} diff --git a/Sources/ForceSimulation/SimulationKD/SimulationKD.LinkForce.swift b/Sources/ForceSimulation/SimulationKD/SimulationKD.LinkForce.swift new file mode 100644 index 0000000..5f243d0 --- /dev/null +++ b/Sources/ForceSimulation/SimulationKD/SimulationKD.LinkForce.swift @@ -0,0 +1,240 @@ +// +// LinkForce.swift +// +// +// Created by li3zhen1 on 10/16/23. +// + +import NDTree + +enum LinkForceError: Error { + case useBeforeSimulationInitialized +} + +extension SimulationKD { + /// A force that represents links between nodes. + /// The complexity is `O(e)`, where `e` is the number of links. + /// See [Link Force - D3](https://d3js.org/d3-force/link). + final public class LinkForce: ForceLike + where NodeID: Hashable, V: VectorLike, V.Scalar : SimulatableFloatingPoint { + + /// + public enum LinkStiffness { + case constant(V.Scalar) + case varied((EdgeID, LinkLookup) -> V.Scalar) + case weightedByDegree(k: (EdgeID, LinkLookup) -> V.Scalar) + } + var linkStiffness: LinkStiffness + var calculatedStiffness: [V.Scalar] = [] + + /// + public typealias LengthScalar = V.Scalar + public enum LinkLength { + case constant(LengthScalar) + case varied((EdgeID, LinkLookup) -> LengthScalar) + } + var linkLength: LinkLength + var calculatedLength: [LengthScalar] = [] + + /// Bias + var calculatedBias: [V.Scalar] = [] + + /// Binding to simulation + /// + weak var simulation: SimulationKD? { + didSet { + + guard let sim = simulation else { return } + + linksOfIndices = links.map { l in + EdgeID( + sim.nodeIdToIndexLookup[l.source, default: 0], + sim.nodeIdToIndexLookup[l.target, default: 0] + ) + } + + self.lookup = .buildFromLinks(linksOfIndices) + + self.calculatedBias = linksOfIndices.map { l in + V.Scalar(lookup.count[l.source, default: 0]) + / V.Scalar( + lookup.count[l.target, default: 0] + lookup.count[l.source, default: 0]) + } + + let lookupWithOriginalID = LinkLookup.buildFromLinks(links) + self.calculatedLength = linkLength.calculated( + for: self.links, connectionLookupTable: lookupWithOriginalID) + self.calculatedStiffness = linkStiffness.calculated( + for: self.links, connectionLookupTable: lookupWithOriginalID) + } + } + + var iterationsPerTick: UInt + + internal var linksOfIndices: [EdgeID] = [] + var links: [EdgeID] + + public struct LinkLookup<_NodeID> where _NodeID: Hashable { + let sources: [_NodeID: [_NodeID]] + let targets: [_NodeID: [_NodeID]] + let count: [_NodeID: Int] + } + private var lookup = LinkLookup(sources: [:], targets: [:], count: [:]) + + internal init( + _ links: [EdgeID], + stiffness: LinkStiffness, + originalLength: LinkLength = .constant(30), + iterationsPerTick: UInt = 1 + ) { + self.links = links + self.iterationsPerTick = iterationsPerTick + self.linkStiffness = stiffness + self.linkLength = originalLength + + } + + public func apply() { + guard let sim = self.simulation else { return } + + let alpha = sim.alpha + + for _ in 0..], + stiffness: LinkForce.LinkStiffness = .weightedByDegree { _, _ in 1.0 }, + originalLength: LinkForce.LinkLength = .constant(30.0), + iterationsPerTick: UInt = 1 + ) -> LinkForce { + let linkForce = LinkForce( + links, stiffness: stiffness, originalLength: originalLength) + linkForce.simulation = self + self.forces.append(linkForce) + return linkForce + } + + /// Create a link force that represents links between nodes. It works like + /// there is a spring between each pair of nodes. + /// The complexity is `O(e)`, where `e` is the number of links. + /// See [Link Force - D3](https://d3js.org/d3-force/link). + /// - Parameters: + /// - links: The links between nodes. + /// - stiffness: The stiffness of the spring (or links). + /// - originalLength: The original length of the spring (or links). + @discardableResult + public func createLinkForce( + _ linkTuples: [(NodeID, NodeID)], + stiffness: LinkForce.LinkStiffness = .weightedByDegree { _, _ in 1.0 }, + originalLength: LinkForce.LinkLength = .constant(30.0), + iterationsPerTick: UInt = 1 + ) -> LinkForce { + let links = linkTuples.map { EdgeID($0.0, $0.1) } + let linkForce = LinkForce( + links, stiffness: stiffness, originalLength: originalLength) + linkForce.simulation = self + self.forces.append(linkForce) + return linkForce + } +} + +extension SimulationKD.LinkForce.LinkLookup { + static func buildFromLinks(_ links: [EdgeID<_NodeID>]) -> Self { + var sources: [_NodeID: [_NodeID]] = [:] + var targets: [_NodeID: [_NodeID]] = [:] + var count: [_NodeID: Int] = [:] + for link in links { + sources[link.source, default: []].append(link.target) + targets[link.target, default: []].append(link.source) + count[link.source, default: 0] += 1 + count[link.target, default: 0] += 1 + } + return Self(sources: sources, targets: targets, count: count) + } +} + +extension SimulationKD.LinkForce.LinkLength { + func calculated( + for links: [EdgeID], + connectionLookupTable: SimulationKD.LinkForce.LinkLookup + ) -> [V.Scalar] { + switch self { + case .constant(let value): + return links.map { _ in value } + case .varied(let f): + return links.map { link in + f(link, connectionLookupTable) + } + } + } +} + +extension SimulationKD.LinkForce.LinkStiffness { + func calculated( + for links: [EdgeID], + connectionLookupTable lookup: SimulationKD.LinkForce.LinkLookup + ) -> [V.Scalar] { + switch self { + case .constant(let value): + return links.map { _ in value } + case .varied(let f): + return links.map { link in + f(link, lookup) + } + case .weightedByDegree(let k): + return links.map { link in + k(link, lookup) + / V.Scalar( + min( + lookup.count[link.source, default: 0], + lookup.count[link.target, default: 0] + ) + ) + } + } + } +} diff --git a/Sources/ForceSimulation/SimulationKD/SimulationKD.ManyBodyForce.swift b/Sources/ForceSimulation/SimulationKD/SimulationKD.ManyBodyForce.swift new file mode 100644 index 0000000..d85d833 --- /dev/null +++ b/Sources/ForceSimulation/SimulationKD/SimulationKD.ManyBodyForce.swift @@ -0,0 +1,334 @@ +// +// File.swift +// +// +// Created by li3zhen1 on 10/1/23. +// + +import NDTree + +enum ManyBodyForceError: Error { + case buildQuadTreeBeforeSimulationInitialized +} + +struct MassQuadtreeDelegate: NDTreeDelegate where NodeID: Hashable, V: VectorLike { + + public var accumulatedMass: V.Scalar = .zero + public var accumulatedCount = 0 + public var accumulatedMassWeightedPositions: V = .zero + + @usableFromInline let massProvider: (NodeID) -> V.Scalar + + init( + massProvider: @escaping (NodeID) -> V.Scalar + ) { + self.massProvider = massProvider + } + + internal init( + initialAccumulatedProperty: V.Scalar, + initialAccumulatedCount: Int, + initialWeightedAccumulatedNodePositions: V, + massProvider: @escaping (NodeID) -> V.Scalar + ) { + self.accumulatedMass = initialAccumulatedProperty + self.accumulatedCount = initialAccumulatedCount + self.accumulatedMassWeightedPositions = initialWeightedAccumulatedNodePositions + self.massProvider = massProvider + } + + @inlinable mutating func didAddNode(_ node: NodeID, at position: V) { + let p = massProvider(node) + #if DEBUG + assert(p > 0) + #endif + accumulatedCount += 1 + accumulatedMass += p + accumulatedMassWeightedPositions += position * p + } + + @inlinable mutating func didRemoveNode(_ node: NodeID, at position: V) { + let p = massProvider(node) + accumulatedCount -= 1 + accumulatedMass -= p + accumulatedMassWeightedPositions -= position * p + // TODO: parent removal? + } + + @inlinable func copy() -> Self { + return Self( + initialAccumulatedProperty: self.accumulatedMass, + initialAccumulatedCount: self.accumulatedCount, + initialWeightedAccumulatedNodePositions: self.accumulatedMassWeightedPositions, + massProvider: self.massProvider + ) + } + + @inlinable func spawn() -> Self { + return Self(massProvider: self.massProvider) + } + + @inlinable var centroid: V? { + guard accumulatedCount > 0 else { return nil } + return accumulatedMassWeightedPositions / accumulatedMass + } +} + +extension SimulationKD { + /// A force that simulate the many-body force. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. The complexity might degrade to `O(n^2)` if the nodes are too close to each other. + /// See [Manybody Force - D3](https://d3js.org/d3-force/many-body). + final public class ManyBodyForce: ForceLike + where NodeID: Hashable, V: VectorLike, V.Scalar : SimulatableFloatingPoint { + + var strength: V.Scalar + + public enum NodeMass { + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) + } + var mass: NodeMass + var precalculatedMass: [V.Scalar] = [] + + weak var simulation: SimulationKD? { + didSet { + guard let sim = self.simulation else { return } + self.precalculatedMass = self.mass.calculated(for: sim) + self.forces = [V](repeating: .zero, count: sim.nodePositions.count) + } + } + + var theta2: V.Scalar + var theta: V.Scalar { + didSet { + theta2 = theta * theta + } + } + + var distanceMin2: V.Scalar = 1 + var distanceMax2: V.Scalar = V.Scalar.infinity + var distanceMin: V.Scalar = 1 + var distanceMax: V.Scalar = V.Scalar.infinity + + internal init( + strength: V.Scalar, + nodeMass: NodeMass = .constant(1.0), + theta: V.Scalar = 0.9 + ) { + self.strength = strength + self.mass = nodeMass + self.theta = theta + self.theta2 = theta * theta + } + + var forces: [V] = [] + public func apply() { + guard let simulation else { return } + + let alpha = simulation.alpha + + try! calculateForce(alpha: alpha) //else { return } + + for i in simulation.nodeVelocities.indices { + simulation.nodeVelocities[i] += self.forces[i] / precalculatedMass[i] + } + } + + // private func getCoveringBox() throws -> NDBox { + // guard let simulation else { throw ManyBodyForceError.buildQuadTreeBeforeSimulationInitialized } + // var _p0 = simulation.nodes[0].position + // var _p1 = simulation.nodes[0].position + // + // for p in simulation.nodes { + // for i in 0..= _p1[i] { + // _p1[i] = p.position[i] + 1 + // } + // } + // } + // return NDBox(_p0, _p1) + // + // } + + func calculateForce(alpha: V.Scalar) throws { + + guard let sim = self.simulation else { + throw ManyBodyForceError.buildQuadTreeBeforeSimulationInitialized + } + + let coveringBox = NDBox.cover(of: sim.nodePositions) //try! getCoveringBox() + + let tree = NDTree>( + box: coveringBox, clusterDistance: 1e-5 + ) { + + return switch self.mass { + case .constant(let m): + MassQuadtreeDelegate { _ in m } + case .varied(_): + MassQuadtreeDelegate { index in + self.precalculatedMass[index] + } + } + } + + for i in sim.nodePositions.indices { + tree.add(i, at: sim.nodePositions[i]) + + #if DEBUG + assert(tree.delegate.accumulatedCount == i + 1) + #endif + + } + + // var forces = [V](repeating: .zero, count: sim.nodePositions.count) + + for i in sim.nodePositions.indices { + var f = V.zero + tree.visit { t in + + // guard t.delegate.accumulatedCount > 0 else { return false } + // + // let centroid = t.delegate.accumulatedMassWeightedPositions / t.delegate.accumulatedMass + // let vec = centroid - sim.nodePositions[i] + // + // var distanceSquared = vec.jiggled().lengthSquared() + // + // /// too far away, omit + // guard distanceSquared < self.distanceMax2 else { return false } + // + // + // + // /// too close, enlarge distance + // if distanceSquared < self.distanceMin2 { + // distanceSquared = (self.distanceMin2 * distanceSquared).squareRoot() + // } + // + // + // if t.nodePosition != nil { + // + // /// filled leaf + // if !t.nodeIndices.contains(i) { + // let k: V.Scalar = self.strength * alpha * t.delegate.accumulatedMass / distanceSquared / (distanceSquared).squareRoot() + // forces[i] += vec * k + // } + // + // return false + // + // } + // else if t.children != nil { + // + // let boxWidth = (t.box.p1 - t.box.p1)[0] + // + // /// internal, guard in 180 guarantees we have nodes here + // if distanceSquared * self.theta2 > boxWidth * boxWidth { + // // far enough + // let k: V.Scalar = self.strength * alpha * t.delegate.accumulatedMass / distanceSquared / (distanceSquared).squareRoot() + // forces[i] += vec * k + // return false + // } + // else { + // return true + // } + // } + // else { + // // empty leaf + // return false + // } + + guard t.delegate.accumulatedCount > 0 else { return false } + let centroid = + t.delegate.accumulatedMassWeightedPositions / t.delegate.accumulatedMass + + let vec = centroid - sim.nodePositions[i] + let boxWidth = (t.box.p1 - t.box.p0)[0] + var distanceSquared = vec.jiggled().lengthSquared() + + let farEnough: Bool = (distanceSquared * self.theta2) > (boxWidth * boxWidth) + + // let distance = distanceSquared.squareRoot() + + if distanceSquared < self.distanceMin2 { + distanceSquared = (self.distanceMin2 * distanceSquared).squareRoot() + } + + if farEnough { + + guard distanceSquared < self.distanceMax2 else { return true } + + /// Workaround for "The compiler is unable to type-check this expression in reasonable time; try breaking up the expression into distinct sub-expressions" + let k: V.Scalar = + self.strength * alpha * t.delegate.accumulatedMass / distanceSquared // distanceSquared.squareRoot() + + f += vec * k + return false + + } else if t.children != nil { + return true + } + + if t.isFilledLeaf { + + // for j in t.nodeIndices { + // if j != i { + // let k: V.Scalar = + // self.strength * alpha * self.precalculatedMass[j] / distanceSquared / distanceSquared.squareRoot() + // f += vec * k + // } + // } + if t.nodeIndices.contains(i) { return false } + + let massAcc = t.delegate.accumulatedMass + // t.nodeIndices.contains(i) ? (t.delegate.accumulatedMass-self.precalculatedMass[i]) : (t.delegate.accumulatedMass) + let k: V.Scalar = self.strength * alpha * massAcc / distanceSquared // distanceSquared.squareRoot() + f += vec * k + return false + } else { + return true + } + } + forces[i] = f + } + // return forces + } + + } + + /// Create a many-body force that simulate the many-body force. + /// This is a very expensive force, the complexity is `O(n log(n))`, + /// where `n` is the number of nodes. The complexity might degrade to `O(n^2)` if the nodes are too close to each other. + /// The force mimics the gravity force or electrostatic force. + /// See [Manybody Force - D3](https://d3js.org/d3-force/many-body). + /// - Parameters: + /// - strength: The strength of the force. When the strength is positive, the nodes are attracted to each other like gravity force, otherwise, the nodes are repelled like electrostatic force. + /// - nodeMass: The mass of the nodes. The mass is used to calculate the force. The default value is 1.0. + /// - theta: Determines how approximate the calculation is. The default value is 0.9. The higher the value, the more approximate and fast the calculation is. + @discardableResult + public func createManyBodyForce( + strength: V.Scalar, + nodeMass: ManyBodyForce.NodeMass = .constant(1.0) + ) -> ManyBodyForce { + let manyBodyForce = ManyBodyForce( + strength: strength, nodeMass: nodeMass) + manyBodyForce.simulation = self + self.forces.append(manyBodyForce) + return manyBodyForce + } +} + +extension SimulationKD.ManyBodyForce.NodeMass { + public func calculated(for simulation: SimulationKD) -> [V.Scalar] { + switch self { + case .constant(let m): + return Array(repeating: m, count: simulation.nodePositions.count) + case .varied(let massGetter): + return simulation.nodeIds.map { n in + return massGetter(n) + } + } + } +} diff --git a/Sources/ForceSimulation/forces/RadialForce.swift b/Sources/ForceSimulation/SimulationKD/SimulationKD.RadialForce.swift similarity index 73% rename from Sources/ForceSimulation/forces/RadialForce.swift rename to Sources/ForceSimulation/SimulationKD/SimulationKD.RadialForce.swift index 5e59213..ca03b39 100644 --- a/Sources/ForceSimulation/forces/RadialForce.swift +++ b/Sources/ForceSimulation/SimulationKD/SimulationKD.RadialForce.swift @@ -7,13 +7,16 @@ import NDTree + +extension SimulationKD { + /// A force that applies a radial force to all nodes. /// Center force is relatively fast, the complexity is `O(n)`, /// where `n` is the number of nodes. /// See [Position Force - D3](https://d3js.org/d3-force/position). -final public class RadialForce: ForceLike -where NodeID: Hashable, V: VectorLike, V.Scalar == Double { - weak var simulation: Simulation? { +final public class RadialForce: ForceLike +where NodeID: Hashable, V: VectorLike, V.Scalar : SimulatableFloatingPoint { + weak var simulation: SimulationKD? { didSet { guard let sim = self.simulation else { return } self.calculatedStrength = strength.calculated(for: sim) @@ -33,11 +36,11 @@ where NodeID: Hashable, V: VectorLike, V.Scalar == Double { /// Strength accessor public enum Strength { - case constant(Double) - case varied((NodeID) -> Double) + case constant(V.Scalar) + case varied((NodeID) -> V.Scalar) } public var strength: Strength - private var calculatedStrength: [Double] = [] + private var calculatedStrength: [V.Scalar] = [] public init(center: V, radius: NodeRadius, strength: Strength) { self.center = center @@ -45,8 +48,9 @@ where NodeID: Hashable, V: VectorLike, V.Scalar == Double { self.strength = strength } - public func apply(alpha: Double) { + public func apply() { guard let sim = self.simulation else { return } + let alpha = sim.alpha for i in sim.nodePositions.indices { let nodeId = i let deltaPosition = (sim.nodePositions[i] - self.center).jiggled() @@ -60,29 +64,8 @@ where NodeID: Hashable, V: VectorLike, V.Scalar == Double { } -extension RadialForce.Strength: PrecalculatableNodeProperty { - public func calculated(for simulation: Simulation) -> [Double] { - switch self { - case .constant(let s): - return simulation.nodeIds.map { _ in s } - case .varied(let s): - return simulation.nodeIds.map(s) - } - } -} -extension RadialForce.NodeRadius: PrecalculatableNodeProperty { - public func calculated(for simulation: Simulation) -> [Double] { - switch self { - case .constant(let r): - return simulation.nodeIds.map { _ in r } - case .varied(let r): - return simulation.nodeIds.map(r) - } - } -} -extension Simulation { /// Create a radial force that applies a radial force to all nodes. /// Center force is relatively fast, the complexity is `O(n)`, @@ -95,13 +78,37 @@ extension Simulation { @discardableResult public func createRadialForce( center: V = .zero, - radius: RadialForce.NodeRadius, - strength: RadialForce.Strength = .constant(0.1) - ) -> RadialForce { - let f = RadialForce(center: center, radius: radius, strength: strength) + radius: RadialForce.NodeRadius, + strength: RadialForce.Strength = .constant(0.1) + ) -> RadialForce { + let f = RadialForce(center: center, radius: radius, strength: strength) f.simulation = self self.forces.append(f) return f } } + + + +extension SimulationKD.RadialForce.Strength { + public func calculated(for simulation: SimulationKD) -> [V.Scalar] { + switch self { + case .constant(let s): + return simulation.nodeIds.map { _ in s } + case .varied(let s): + return simulation.nodeIds.map(s) + } + } +} + +extension SimulationKD.RadialForce.NodeRadius { + public func calculated(for simulation: SimulationKD) -> [V.Scalar] { + switch self { + case .constant(let r): + return simulation.nodeIds.map { _ in r } + case .varied(let r): + return simulation.nodeIds.map(r) + } + } +} \ No newline at end of file diff --git a/Sources/ForceSimulation/Simulation.swift b/Sources/ForceSimulation/SimulationKD/SimulationKD.swift similarity index 76% rename from Sources/ForceSimulation/Simulation.swift rename to Sources/ForceSimulation/SimulationKD/SimulationKD.swift index 305fccf..af553f8 100644 --- a/Sources/ForceSimulation/Simulation.swift +++ b/Sources/ForceSimulation/SimulationKD/SimulationKD.swift @@ -7,39 +7,40 @@ import NDTree -enum SimulationError: Error { +enum SimulationKDError: Error { case subscriptionToNonexistentNode } /// An N-Dimensional force simulation. -public final class Simulation where NodeID: Hashable, V: VectorLike, V.Scalar == Double { +public final class SimulationKD +where NodeID: Hashable, V: VectorLike, V.Scalar : SimulatableFloatingPoint { /// The type of the vector used in the simulation. - /// Usually this is `Double` if you are on Apple platforms. + /// Usually this is `Scalar` if you are on Apple platforms. public typealias Scalar = V.Scalar - public let initializedAlpha: Double + public let initializedAlpha: Scalar - public var alpha: Double - public var alphaMin: Double - public var alphaDecay: Double - public var alphaTarget: Double - - public var velocityDecay: V.Scalar + public var alpha: Scalar + public var alphaMin: Scalar + public var alphaDecay: Scalar + public var alphaTarget: Scalar + public var velocityDecay: Scalar + + public internal(set) var forces: [any ForceLike] = [] - /// The position of points stored in simulation. + /// The position of points stored in simulation. /// Ordered as the nodeIds you passed in when initializing simulation. /// They are always updated. public internal(set) var nodePositions: [V] - + /// The velocities of points stored in simulation. /// Ordered as the nodeIds you passed in when initializing simulation. /// They are always updated. public internal(set) var nodeVelocities: [V] - - + /// The fixed positions of points stored in simulation. /// Ordered as the nodeIds you passed in when initializing simulation. /// They are always updated. @@ -60,11 +61,11 @@ public final class Simulation where NodeID: Hashable, V: VectorLike, /// - getInitialPosition: The closure to set the initial position of the node. If not provided, the initial position is set to zero. public init( nodeIds: [NodeID], - alpha: Double = 1, - alphaMin: Double = 1e-3, - alphaDecay: Double = 2e-3, - alphaTarget: Double = 0.0, - velocityDecay: Double = 0.6, + alpha: Scalar = 1, + alphaMin: Scalar = 1e-3, + alphaDecay: Scalar = 2e-3, + alphaTarget: Scalar = 0.0, + velocityDecay: Scalar = 0.6, setInitialStatus getInitialPosition: ( (NodeID) -> V @@ -103,9 +104,9 @@ public final class Simulation where NodeID: Hashable, V: VectorLike, public func getIndex(of nodeId: NodeID) -> Int { return nodeIdToIndexLookup[nodeId]! } - + /// Reset the alpha. The points will move faster as alpha gets larger. - public func resetAlpha(_ alpha: Double) { + public func resetAlpha(_ alpha: Scalar) { self.alpha = alpha } @@ -118,7 +119,7 @@ public final class Simulation where NodeID: Hashable, V: VectorLike, alpha += (alphaTarget - alpha) * alphaDecay for f in forces { - f.apply(alpha: alpha) + f.apply() } for i in nodePositions.indices { @@ -133,13 +134,3 @@ public final class Simulation where NodeID: Hashable, V: VectorLike, } } } - -#if canImport(simd) - - /// A 2D simulation running on `Double` and `simd_double2` types. - public typealias Simulation2D = Simulation where NodeID: Hashable - - /// A 3D simulation running on `Double` and `simd_double3` types. - public typealias Simulation3D = Simulation where NodeID: Hashable - -#endif diff --git a/Sources/ForceSimulation/Utils.swift b/Sources/ForceSimulation/Utils.swift index 92cfc05..4242871 100644 --- a/Sources/ForceSimulation/Utils.swift +++ b/Sources/ForceSimulation/Utils.swift @@ -7,6 +7,7 @@ import NDTree + // TODO: https://forums.swift.org/t/deterministic-randomness-in-swift/20835/5 /// A random number generator that generates deterministic random numbers. @@ -46,7 +47,11 @@ public struct FloatLinearCongruentialGenerator { } } -extension Double { +public protocol SimulatableFloatingPoint: ExpressibleByFloatLiteral { + func jiggled() -> Self +} + +extension Double: SimulatableFloatingPoint { @inlinable public func jiggled() -> Double { if self == 0 || self == .nan { // return Double.random(in: -5e-6..<5e-6) @@ -56,7 +61,7 @@ extension Double { } } -extension Float { +extension Float: SimulatableFloatingPoint { @inlinable public func jiggled() -> Float { if self == 0 || self == .nan { // return Double.random(in: -5e-6..<5e-6) @@ -66,7 +71,11 @@ extension Float { } } -extension VectorLike where Scalar == Double { + +#if canImport(simd) + +import simd +extension simd_double2 { @inlinable public func jiggled() -> Self { var result = Self.zero for i in indices { @@ -76,7 +85,19 @@ extension VectorLike where Scalar == Double { } } -extension VectorLike where Scalar == Float { +extension simd_float3 { + @inlinable public func jiggled() -> Self { + var result = Self.zero + for i in indices { + result[i] = self[i].jiggled() + } + return result + } +} + +#endif + +extension VectorLike where Scalar: SimulatableFloatingPoint { @inlinable public func jiggled() -> Self { var result = Self.zero for i in indices { @@ -98,12 +119,3 @@ public struct EdgeID: Hashable where NodeID: Hashable { self.target = target } } - - -public protocol PrecalculatableNodeProperty { - associatedtype NodeID: Hashable - - associatedtype V: VectorLike where V.Scalar == Double - - func calculated(for simulation: Simulation) -> [Double] -} diff --git a/Sources/ForceSimulation/forces/CenterForce.swift b/Sources/ForceSimulation/forces/CenterForce.swift deleted file mode 100644 index 4299a10..0000000 --- a/Sources/ForceSimulation/forces/CenterForce.swift +++ /dev/null @@ -1,57 +0,0 @@ -// -// CenterForce.swift -// -// -// Created by li3zhen1 on 10/16/23. -// -import NDTree - -/// A force that drives nodes towards the center. -/// Center force is relatively fast, the complexity is $O(n)$, -/// where `n` is the number of nodes. -/// See [Collide Force - D3](https://d3js.org/d3-force/collide). -final public class CenterForce: ForceLike -where NodeID: Hashable, V: VectorLike, V.Scalar == Double { - public var center: V - public var strength: Double - weak var simulation: Simulation? - - internal init(center: V, strength: Double) { - self.center = center - self.strength = strength - } - - public func apply(alpha: Double) { - guard let sim = self.simulation else { return } - - var meanPosition = V.zero - for n in sim.nodePositions { - meanPosition += n //.position - } - let delta = meanPosition * (self.strength / Double(sim.nodePositions.count)) - - for i in sim.nodePositions.indices { - sim.nodePositions[i] -= delta - } - } - -} - -extension Simulation { - - /// Create a center force that drives nodes towards the center. - /// Center force is relatively fast, the complexity is `O(n)`, - /// where `n` is the number of nodes. - /// See [Collide Force - D3](https://d3js.org/d3-force/collide). - /// - Parameters: - /// - center: The center of the force. - /// - strength: The strength of the force. - @discardableResult - public func createCenterForce(center: V, strength: Double = 0.1) -> CenterForce { - let f = CenterForce(center: center, strength: strength) - f.simulation = self - self.forces.append(f) - return f - } - -} diff --git a/Sources/ForceSimulation/forces/CollideForce.swift b/Sources/ForceSimulation/forces/CollideForce.swift deleted file mode 100644 index 4a88281..0000000 --- a/Sources/ForceSimulation/forces/CollideForce.swift +++ /dev/null @@ -1,208 +0,0 @@ -// -// File.swift -// -// -// Created by li3zhen1 on 10/17/23. -// -import NDTree - -struct MaxRadiusTreeDelegate: NDTreeDelegate where NodeID: Hashable, V: VectorLike { - - public var maxNodeRadius: Double - - @usableFromInline var radiusProvider: (NodeID) -> Double - - mutating func didAddNode(_ nodeId: NodeID, at position: V) { - let p = radiusProvider(nodeId) - maxNodeRadius = max(maxNodeRadius, p) - } - - mutating func didRemoveNode(_ nodeId: NodeID, at position: V) { - if radiusProvider(nodeId) >= maxNodeRadius { - // 🤯 for Collide force, set to 0 is fine - // Otherwise you need to traverse the delegate again - maxNodeRadius = 0 - } - } - - func copy() -> MaxRadiusTreeDelegate { - return Self(maxNodeRadius: maxNodeRadius, radiusProvider: radiusProvider) - } - - func spawn() -> MaxRadiusTreeDelegate { - return Self(radiusProvider: radiusProvider) - } - - init(maxNodeRadius: Double = 0, radiusProvider: @escaping (NodeID) -> Double) { - self.maxNodeRadius = maxNodeRadius - self.radiusProvider = radiusProvider - } - -} - -/// A force that prevents nodes from overlapping. -/// This is a very expensive force, the complexity is `O(n log(n))`, -/// where `n` is the number of nodes. -/// See [Collide Force - D3](https://d3js.org/d3-force/collide). -public final class CollideForce: ForceLike -where NodeID: Hashable, V: VectorLike, V.Scalar == Double { - - weak var simulation: Simulation? { - didSet { - guard let sim = simulation else { return } - self.calculatedRadius = radius.calculated(for: sim) - } - } - - public enum CollideRadius { - case constant(Double) - case varied((NodeID) -> Double) - } - public var radius: CollideRadius - var calculatedRadius: [Double] = [] - - public let iterationsPerTick: UInt - public var strength: Double - - internal init( - radius: CollideRadius, - strength: Double = 1.0, - iterationsPerTick: UInt = 1 - ) { - self.radius = radius - self.iterationsPerTick = iterationsPerTick - self.strength = strength - } - - public func apply(alpha: Double) { - guard let sim = self.simulation else { return } - - for _ in 0...cover(of: sim.nodePositions) - - let tree = NDTree>( - box: coveringBox, clusterDistance: 1e-5 - ) { - return switch self.radius { - case .constant(let m): - MaxRadiusTreeDelegate { _ in m } - case .varied(_): - MaxRadiusTreeDelegate { index in - self.calculatedRadius[index] - } - } - } - - for i in sim.nodePositions.indices { - tree.add(i, at: sim.nodePositions[i]) - } - - for i in sim.nodePositions.indices { - let iOriginalPosition = sim.nodePositions[i] - let iOriginalVelocity = sim.nodeVelocities[i] - let iR = self.calculatedRadius[i] - let iR2 = iR * iR - let iPosition = iOriginalPosition + iOriginalVelocity - - tree.visit { t in - - let maxRadiusOfQuad = t.delegate.maxNodeRadius - let deltaR = maxRadiusOfQuad + iR - - if t.nodePosition != nil { - for j in t.nodeIndices { - // print("\(i)<=>\(j)") - // is leaf, make sure every collision happens once. - guard j > i else { continue } - - let jR = self.calculatedRadius[j] - let jOriginalPosition = sim.nodePositions[j] - let jOriginalVelocity = sim.nodeVelocities[j] - var deltaPosition = iPosition - (jOriginalPosition + jOriginalVelocity) - let l = deltaPosition.lengthSquared() - - let deltaR = iR + jR - if l < deltaR * deltaR { - - var l = deltaPosition.jiggled().length() - l = (deltaR - l) / l * self.strength - - let jR2 = jR * jR - - let k = jR2 / (iR2 + jR2) - - deltaPosition *= l - - sim.nodeVelocities[i] += deltaPosition * k - sim.nodeVelocities[j] -= deltaPosition * (1 - k) - } - } - return false - } - - for laneIndex in t.box.p0.indices { - let _v = t.box.p0[laneIndex] - if _v > iPosition[laneIndex] + deltaR /* True if no overlap */ { - return false - } - } - - for laneIndex in t.box.p1.indices { - let _v = t.box.p1[laneIndex] - if _v < iPosition[laneIndex] - deltaR /* True if no overlap */ { - return false - } - } - return true - - // return - // !(t.quad.x0 > iPosition.x + deltaR /* True if no overlap */ - // || t.quad.x1 < iPosition.x - deltaR - // || t.quad.y0 > iPosition.y + deltaR - // || t.quad.y1 < iPosition.y - deltaR) - } - } - } - } - -} - -extension CollideForce.CollideRadius: PrecalculatableNodeProperty { - public func calculated(for simulation: Simulation) -> [Double] { - switch self { - case .constant(let r): - return Array(repeating: r, count: simulation.nodePositions.count) - case .varied(let radiusProvider): - return simulation.nodeIds.map { radiusProvider($0) } - } - } -} - -extension Simulation { - - /// Create a collide force that prevents nodes from overlapping. - /// This is a very expensive force, the complexity is `O(n log(n))`, - /// where `n` is the number of nodes. - /// See [Collide Force - D3](https://d3js.org/d3-force/collide). - /// - Parameters: - /// - radius: The radius of the force. - /// - strength: The strength of the force. - /// - iterationsPerTick: The number of iterations per tick. - @discardableResult - public func createCollideForce( - radius: CollideForce.CollideRadius = .constant(3.0), - strength: Double = 1.0, - iterationsPerTick: UInt = 1 - ) -> CollideForce { - let f = CollideForce( - radius: radius, - strength: strength, - iterationsPerTick: iterationsPerTick - ) - f.simulation = self - self.forces.append(f) - return f - } - -} diff --git a/Sources/ForceSimulation/forces/DirectionForce.swift b/Sources/ForceSimulation/forces/DirectionForce.swift deleted file mode 100644 index 002af6e..0000000 --- a/Sources/ForceSimulation/forces/DirectionForce.swift +++ /dev/null @@ -1,120 +0,0 @@ -// -// PositionForce.swift -// -// -// Created by li3zhen1 on 10/1/23. -// - -import NDTree - -/// A force that moves nodes to a target position. -/// Center force is relatively fast, the complexity is `O(n)`, -/// where `n` is the number of nodes. -/// See [Position Force - D3](https://d3js.org/d3-force/position). -final public class DirectionForce: ForceLike -where NodeID: Hashable, V: VectorLike, V.Scalar == Double { - - public enum Direction { - case x - case y - case entryOfVector(Int) - } - public enum Strength { - case constant(Double) - case varied((NodeID) -> Double) - } - - public enum TargetOnDirection { - case constant(V.Scalar) - case varied((NodeID) -> V.Scalar) - } - - public var strength: Strength - public var direction: Int - public var calculatedStrength: [Double] = [] - public var targetOnDirection: TargetOnDirection - public var calculatedTargetOnDirection: [V.Scalar] = [] - - internal init( - direction: Direction, targetOnDirection: TargetOnDirection, - strength: Strength = .constant(1.0) - ) { - - self.strength = strength - self.direction = direction.lane - self.targetOnDirection = targetOnDirection - } - - weak var simulation: Simulation? { - didSet { - guard let sim = self.simulation else { return } - self.calculatedStrength = strength.calculated(for: sim) - self.calculatedTargetOnDirection = targetOnDirection.calculated(for: sim) - } - } - - public func apply(alpha: Double) { - guard let sim = self.simulation else { return } - let lane = self.direction - for i in sim.nodePositions.indices { - sim.nodeVelocities[i][lane] += - (self.calculatedTargetOnDirection[i] - sim.nodePositions[i][lane]) - * self.calculatedStrength[i] * alpha - } - } -} - -extension DirectionForce.Strength: PrecalculatableNodeProperty { - public func calculated(for simulation: Simulation) -> [Double] { - switch self { - case .constant(let value): - return Array(repeating: value, count: simulation.nodeIds.count) - case .varied(let getter): - return simulation.nodeIds.map(getter) - } - } -} - -extension DirectionForce.TargetOnDirection: PrecalculatableNodeProperty { - public func calculated(for simulation: Simulation) -> [Double] { - switch self { - case .constant(let value): - return Array(repeating: value, count: simulation.nodeIds.count) - case .varied(let getter): - return simulation.nodeIds.map(getter) - } - } -} - -extension DirectionForce.Direction { - @inlinable var lane: Int { - switch self { - case .x: return 0 - case .y: return 1 - case .entryOfVector(let i): return i - } - } -} - -extension Simulation { - - /// Create a direction force that moves nodes to a target position. - /// Center force is relatively fast, the complexity is `O(n)`, - /// where `n` is the number of nodes. - /// See [Position Force - D3](https://d3js.org/d3-force/position). - @discardableResult - public func createPositionForce( - direction: DirectionForce.Direction, - targetOnDirection: DirectionForce.TargetOnDirection, - strength: DirectionForce.Strength = .constant(1.0) - ) -> DirectionForce { - let force = DirectionForce( - direction: direction, - targetOnDirection: targetOnDirection, - strength: strength - ) - force.simulation = self - self.forces.append(force) - return force - } -} diff --git a/Sources/ForceSimulation/forces/LinkForce.swift b/Sources/ForceSimulation/forces/LinkForce.swift deleted file mode 100644 index 92f2b28..0000000 --- a/Sources/ForceSimulation/forces/LinkForce.swift +++ /dev/null @@ -1,246 +0,0 @@ -// -// LinkForce.swift -// -// -// Created by li3zhen1 on 10/16/23. -// - -import NDTree - -enum LinkForceError: Error { - case useBeforeSimulationInitialized -} - -/// A force that represents links between nodes. -/// The complexity is `O(e)`, where `e` is the number of links. -/// See [Link Force - D3](https://d3js.org/d3-force/link). -final public class LinkForce: ForceLike -where NodeID: Hashable, V: VectorLike, V.Scalar == Double { - - /// - public enum LinkStiffness { - case constant(Double) - case varied((EdgeID, LinkLookup) -> Double) - case weightedByDegree(k: (EdgeID, LinkLookup) -> Double) - } - var linkStiffness: LinkStiffness - var calculatedStiffness: [Double] = [] - - /// - public typealias LengthScalar = V.Scalar - public enum LinkLength { - case constant(LengthScalar) - case varied((EdgeID, LinkLookup) -> LengthScalar) - } - var linkLength: LinkLength - var calculatedLength: [LengthScalar] = [] - - /// Bias - var calculatedBias: [Double] = [] - - /// Binding to simulation - /// - weak var simulation: Simulation? { - didSet { - - guard let sim = simulation else { return } - - linksOfIndices = links.map { l in - EdgeID( - sim.nodeIdToIndexLookup[l.source, default: 0], - sim.nodeIdToIndexLookup[l.target, default: 0] - ) - } - - self.lookup = .buildFromLinks(linksOfIndices) - - self.calculatedBias = linksOfIndices.map { l in - Double(lookup.count[l.source, default: 0]) - / Double( - lookup.count[l.target, default: 0] + lookup.count[l.source, default: 0]) - } - - let lookupWithOriginalID = LinkLookup.buildFromLinks(links) - self.calculatedLength = linkLength.calculated( - for: self.links, connectionLookupTable: lookupWithOriginalID) - self.calculatedStiffness = linkStiffness.calculated( - for: self.links, connectionLookupTable: lookupWithOriginalID) - } - } - - var iterationsPerTick: UInt - - internal var linksOfIndices: [EdgeID] = [] - var links: [EdgeID] - - public struct LinkLookup<_NodeID> where _NodeID: Hashable { - let sources: [_NodeID: [_NodeID]] - let targets: [_NodeID: [_NodeID]] - let count: [_NodeID: Int] - } - private var lookup = LinkLookup(sources: [:], targets: [:], count: [:]) - - internal init( - _ links: [EdgeID], - stiffness: LinkStiffness, - originalLength: LinkLength = .constant(30), - iterationsPerTick: UInt = 1 - ) { - self.links = links - self.iterationsPerTick = iterationsPerTick - self.linkStiffness = stiffness - self.linkLength = originalLength - - } - - public func apply(alpha: Double) { - guard let sim = self.simulation else { return } - - for _ in 0..]) -> Self { - var sources: [_NodeID: [_NodeID]] = [:] - var targets: [_NodeID: [_NodeID]] = [:] - var count: [_NodeID: Int] = [:] - for link in links { - sources[link.source, default: []].append(link.target) - targets[link.target, default: []].append(link.source) - count[link.source, default: 0] += 1 - count[link.target, default: 0] += 1 - } - return Self(sources: sources, targets: targets, count: count) - } -} - -protocol PrecalculatableEdgeProperty { - associatedtype NodeID: Hashable - associatedtype V: VectorLike where V.Scalar == Double - func calculated( - for links: [EdgeID], connectionLookupTable: LinkForce.LinkLookup - ) -> [Double] -} - -extension LinkForce.LinkLength: PrecalculatableEdgeProperty { - func calculated( - for links: [EdgeID], connectionLookupTable: LinkForce.LinkLookup - ) -> [V.Scalar] { - switch self { - case .constant(let value): - return links.map { _ in value } - case .varied(let f): - return links.map { link in - f(link, connectionLookupTable) - } - } - } -} - -extension LinkForce.LinkStiffness: PrecalculatableEdgeProperty { - func calculated( - for links: [EdgeID], - connectionLookupTable lookup: LinkForce.LinkLookup - ) -> [Double] { - switch self { - case .constant(let value): - return links.map { _ in value } - case .varied(let f): - return links.map { link in - f(link, lookup) - } - case .weightedByDegree(let k): - return links.map { link in - k(link, lookup) - / Double( - min( - lookup.count[link.source, default: 0], - lookup.count[link.target, default: 0] - ) - ) - } - } - } -} - -extension Simulation { - - /// Create a link force that represents links between nodes. It works like - /// there is a spring between each pair of nodes. - /// The complexity is `O(e)`, where `e` is the number of links. - /// See [Collide Force - D3](https://d3js.org/d3-force/collide) - /// - Parameters: - /// - links: The links between nodes. - /// - stiffness: The stiffness of the spring (or links). - /// - originalLength: The original length of the spring (or links). - @discardableResult - public func createLinkForce( - _ links: [EdgeID], - stiffness: LinkForce.LinkStiffness = .weightedByDegree { _, _ in 1.0 }, - originalLength: LinkForce.LinkLength = .constant(30), - iterationsPerTick: UInt = 1 - ) -> LinkForce { - let linkForce = LinkForce( - links, stiffness: stiffness, originalLength: originalLength) - linkForce.simulation = self - self.forces.append(linkForce) - return linkForce - } - - /// Create a link force that represents links between nodes. It works like - /// there is a spring between each pair of nodes. - /// The complexity is `O(e)`, where `e` is the number of links. - /// See [Link Force - D3](https://d3js.org/d3-force/link). - /// - Parameters: - /// - links: The links between nodes. - /// - stiffness: The stiffness of the spring (or links). - /// - originalLength: The original length of the spring (or links). - @discardableResult - public func createLinkForce( - _ linkTuples: [(NodeID, NodeID)], - stiffness: LinkForce.LinkStiffness = .weightedByDegree { _, _ in 1.0 }, - originalLength: LinkForce.LinkLength = .constant(30), - iterationsPerTick: UInt = 1 - ) -> LinkForce { - let links = linkTuples.map { EdgeID($0.0, $0.1) } - let linkForce = LinkForce( - links, stiffness: stiffness, originalLength: originalLength) - linkForce.simulation = self - self.forces.append(linkForce) - return linkForce - } -} diff --git a/Sources/ForceSimulation/forces/ManyBodyForce.swift b/Sources/ForceSimulation/forces/ManyBodyForce.swift deleted file mode 100644 index 381d830..0000000 --- a/Sources/ForceSimulation/forces/ManyBodyForce.swift +++ /dev/null @@ -1,331 +0,0 @@ -// -// File.swift -// -// -// Created by li3zhen1 on 10/1/23. -// - -import NDTree - -enum ManyBodyForceError: Error { - case buildQuadTreeBeforeSimulationInitialized -} - -struct MassQuadtreeDelegate: NDTreeDelegate where NodeID: Hashable, V: VectorLike { - - public var accumulatedMass: Double = .zero - public var accumulatedCount = 0 - public var accumulatedMassWeightedPositions: V = .zero - - @usableFromInline let massProvider: (NodeID) -> Double - - init( - massProvider: @escaping (NodeID) -> Double - ) { - self.massProvider = massProvider - } - - internal init( - initialAccumulatedProperty: Double, - initialAccumulatedCount: Int, - initialWeightedAccumulatedNodePositions: V, - massProvider: @escaping (NodeID) -> Double - ) { - self.accumulatedMass = initialAccumulatedProperty - self.accumulatedCount = initialAccumulatedCount - self.accumulatedMassWeightedPositions = initialWeightedAccumulatedNodePositions - self.massProvider = massProvider - } - - @inlinable mutating func didAddNode(_ node: NodeID, at position: V) { - let p = massProvider(node) - #if DEBUG - assert(p > 0) - #endif - accumulatedCount += 1 - accumulatedMass += p - accumulatedMassWeightedPositions += position * p - } - - @inlinable mutating func didRemoveNode(_ node: NodeID, at position: V) { - let p = massProvider(node) - accumulatedCount -= 1 - accumulatedMass -= p - accumulatedMassWeightedPositions -= position * p - // TODO: parent removal? - } - - @inlinable func copy() -> Self { - return Self( - initialAccumulatedProperty: self.accumulatedMass, - initialAccumulatedCount: self.accumulatedCount, - initialWeightedAccumulatedNodePositions: self.accumulatedMassWeightedPositions, - massProvider: self.massProvider - ) - } - - @inlinable func spawn() -> Self { - return Self(massProvider: self.massProvider) - } - - @inlinable var centroid: V? { - guard accumulatedCount > 0 else { return nil } - return accumulatedMassWeightedPositions / accumulatedMass - } -} - -/// A force that simulate the many-body force. -/// This is a very expensive force, the complexity is `O(n log(n))`, -/// where `n` is the number of nodes. The complexity might degrade to `O(n^2)` if the nodes are too close to each other. -/// See [Manybody Force - D3](https://d3js.org/d3-force/many-body). -final public class ManyBodyForce: ForceLike -where NodeID: Hashable, V: VectorLike, V.Scalar == Double { - - var strength: Double - - public enum NodeMass { - case constant(Double) - case varied((NodeID) -> Double) - } - var mass: NodeMass = .constant(1.0) - var precalculatedMass: [Double] = [] - - weak var simulation: Simulation? { - didSet { - guard let sim = self.simulation else { return } - self.precalculatedMass = self.mass.calculated(for: sim) - self.forces = [V](repeating: .zero, count: sim.nodePositions.count) - } - } - - var theta2: Double - var theta: Double { - didSet { - theta2 = theta * theta - } - } - - var distanceMin2: Double = 1 - var distanceMax2: Double = Double.infinity - var distanceMin: Double = 1 - var distanceMax: Double = Double.infinity - - internal init( - strength: Double, - nodeMass: NodeMass = .constant(1.0), - theta: Double = 0.9 - ) { - self.strength = strength - self.mass = nodeMass - self.theta = theta - self.theta2 = theta * theta - } - - var forces: [V] = [] - public func apply(alpha: Double) { - guard let simulation else { return } - // guard - try! calculateForce(alpha: alpha) //else { return } - - for i in simulation.nodeVelocities.indices { - simulation.nodeVelocities[i] += self.forces[i] / precalculatedMass[i] - } - } - - // private func getCoveringBox() throws -> NDBox { - // guard let simulation else { throw ManyBodyForceError.buildQuadTreeBeforeSimulationInitialized } - // var _p0 = simulation.nodes[0].position - // var _p1 = simulation.nodes[0].position - // - // for p in simulation.nodes { - // for i in 0..= _p1[i] { - // _p1[i] = p.position[i] + 1 - // } - // } - // } - // return NDBox(_p0, _p1) - // - // } - - func calculateForce(alpha: Double) throws { - - guard let sim = self.simulation else { - throw ManyBodyForceError.buildQuadTreeBeforeSimulationInitialized - } - - let coveringBox = NDBox.cover(of: sim.nodePositions) //try! getCoveringBox() - - let tree = NDTree>(box: coveringBox, clusterDistance: 1e-5) - { - - return switch self.mass { - case .constant(let m): - MassQuadtreeDelegate { _ in m } - case .varied(_): - MassQuadtreeDelegate { index in - self.precalculatedMass[index] - } - } - } - - for i in sim.nodePositions.indices { - tree.add(i, at: sim.nodePositions[i]) - - #if DEBUG - assert(tree.delegate.accumulatedCount == i + 1) - #endif - - } - - // var forces = [V](repeating: .zero, count: sim.nodePositions.count) - - for i in sim.nodePositions.indices { - var f = V.zero - tree.visit { t in - - // guard t.delegate.accumulatedCount > 0 else { return false } - // - // let centroid = t.delegate.accumulatedMassWeightedPositions / t.delegate.accumulatedMass - // let vec = centroid - sim.nodePositions[i] - // - // var distanceSquared = vec.jiggled().lengthSquared() - // - // /// too far away, omit - // guard distanceSquared < self.distanceMax2 else { return false } - // - // - // - // /// too close, enlarge distance - // if distanceSquared < self.distanceMin2 { - // distanceSquared = (self.distanceMin2 * distanceSquared).squareRoot() - // } - // - // - // if t.nodePosition != nil { - // - // /// filled leaf - // if !t.nodeIndices.contains(i) { - // let k: Double = self.strength * alpha * t.delegate.accumulatedMass / distanceSquared / (distanceSquared).squareRoot() - // forces[i] += vec * k - // } - // - // return false - // - // } - // else if t.children != nil { - // - // let boxWidth = (t.box.p1 - t.box.p1)[0] - // - // /// internal, guard in 180 guarantees we have nodes here - // if distanceSquared * self.theta2 > boxWidth * boxWidth { - // // far enough - // let k: Double = self.strength * alpha * t.delegate.accumulatedMass / distanceSquared / (distanceSquared).squareRoot() - // forces[i] += vec * k - // return false - // } - // else { - // return true - // } - // } - // else { - // // empty leaf - // return false - // } - - guard t.delegate.accumulatedCount > 0 else { return false } - let centroid = - t.delegate.accumulatedMassWeightedPositions / t.delegate.accumulatedMass - - let vec = centroid - sim.nodePositions[i] - let boxWidth = (t.box.p1 - t.box.p0)[0] - var distanceSquared = vec.jiggled().lengthSquared() - - let farEnough: Bool = (distanceSquared * self.theta2) > (boxWidth * boxWidth) - - // let distance = distanceSquared.squareRoot() - - if distanceSquared < self.distanceMin2 { - distanceSquared = (self.distanceMin2 * distanceSquared).squareRoot() - } - - if farEnough { - - guard distanceSquared < self.distanceMax2 else { return true } - - /// Workaround for "The compiler is unable to type-check this expression in reasonable time; try breaking up the expression into distinct sub-expressions" - let k: Double = - self.strength * alpha * t.delegate.accumulatedMass / distanceSquared // distanceSquared.squareRoot() - - f += vec * k - return false - - } else if t.children != nil { - return true - } - - if t.isFilledLeaf { - - // for j in t.nodeIndices { - // if j != i { - // let k: Double = - // self.strength * alpha * self.precalculatedMass[j] / distanceSquared / distanceSquared.squareRoot() - // f += vec * k - // } - // } - if t.nodeIndices.contains(i) { return false } - - let massAcc = t.delegate.accumulatedMass - // t.nodeIndices.contains(i) ? (t.delegate.accumulatedMass-self.precalculatedMass[i]) : (t.delegate.accumulatedMass) - let k: Double = self.strength * alpha * massAcc / distanceSquared // distanceSquared.squareRoot() - f += vec * k - return false - } else { - return true - } - } - forces[i] = f - } - // return forces - } - -} - -extension ManyBodyForce.NodeMass: PrecalculatableNodeProperty { - public func calculated(for simulation: Simulation) -> [Double] { - switch self { - case .constant(let m): - return Array(repeating: m, count: simulation.nodePositions.count) - case .varied(let massGetter): - return simulation.nodeIds.map { n in - return massGetter(n) - } - } - } -} - -extension Simulation { - /// Create a many-body force that simulate the many-body force. - /// This is a very expensive force, the complexity is `O(n log(n))`, - /// where `n` is the number of nodes. The complexity might degrade to `O(n^2)` if the nodes are too close to each other. - /// The force mimics the gravity force or electrostatic force. - /// See [Manybody Force - D3](https://d3js.org/d3-force/many-body). - /// - Parameters: - /// - strength: The strength of the force. When the strength is positive, the nodes are attracted to each other like gravity force, otherwise, the nodes are repelled like electrostatic force. - /// - nodeMass: The mass of the nodes. The mass is used to calculate the force. The default value is 1.0. - /// - theta: Determines how approximate the calculation is. The default value is 0.9. The higher the value, the more approximate and fast the calculation is. - @discardableResult - public func createManyBodyForce( - strength: Double, - nodeMass: ManyBodyForce.NodeMass = .constant(1.0) - ) -> ManyBodyForce { - let manyBodyForce = ManyBodyForce( - strength: strength, nodeMass: nodeMass) - manyBodyForce.simulation = self - self.forces.append(manyBodyForce) - return manyBodyForce - } -} diff --git a/Sources/NDTree/NDBox.swift b/Sources/NDTree/NDBox.swift index a5d7d96..6da9947 100644 --- a/Sources/NDTree/NDBox.swift +++ b/Sources/NDTree/NDBox.swift @@ -1,6 +1,6 @@ // // NDBox.swift -// +// // // Created by li3zhen1 on 10/14/23. // @@ -18,7 +18,7 @@ public struct NDBox where V: VectorLike { /// - Parameters: /// - p0: anchor /// - p1: another anchor in the diagonal position of `p0` - /// - Note: `p0` you pass does not have to be minimum point of the box. + /// - Note: `p0` you pass does not have to be minimum point of the box. /// `p1` does not have to be maximum point of the box. The initializer will /// automatically adjust the order of `p0` and `p1` to make sure `p0` is the /// minimum point of the box and `p1` is the maximum point of the box. @@ -62,7 +62,7 @@ public struct NDBox where V: VectorLike { /// - Parameters: /// - p0: anchor /// - p1: another anchor in the diagonal position of `p0` - /// - Note: `p0` you pass does not have to be minimum point of the box. + /// - Note: `p0` you pass does not have to be minimum point of the box. /// `p1` does not have to be maximum point of the box. The initializer will /// automatically adjust the order of `p0` and `p1` to make sure `p0` is the /// minimum point of the box and `p1` is the maximum point of the box. @@ -102,24 +102,22 @@ extension NDBox { } return corner } - - + @inlinable public var debugDescription: String { return "[\(p0), \(p1)]" } } - -public extension NDBox { +extension NDBox { /// Get the small box that contains a list points and guarantees the box's size is at least 1x..x1. /// - Parameter points: The points to be covered. /// - Returns: The box that contains all the points. - @inlinable static func cover(of points: [V]) -> Self { - + @inlinable public static func cover(of points: [V]) -> Self { + var _p0 = points[0] var _p1 = points[0] - + for p in points { for i in p.indices { if p[i] < _p0[i] { @@ -130,29 +128,30 @@ public extension NDBox { } } } - + #if DEBUG - let _box = Self(_p0, _p1) - assert(points.allSatisfy{ p in - _box.contains(p) - }) + let _box = Self(_p0, _p1) + assert( + points.allSatisfy { p in + _box.contains(p) + }) #endif - + return Self(_p0, _p1) } - + /// Get the small box that contains a list points and guarantees the box's size is at least 1x..x1. /// Please note that KeyPath is slow. - /// - /// - Parameter + /// + /// - Parameter /// - points: The points to be covered. /// - keyPath: The key path to get the vector from the point. /// - Returns: The box that contains all the points. - @inlinable static func cover(of points: [T], keyPath: KeyPath) -> Self { - + @inlinable public static func cover(of points: [T], keyPath: KeyPath) -> Self { + var _p0 = points[0][keyPath: keyPath] var _p1 = points[0][keyPath: keyPath] - + for _p in points { let p = _p[keyPath: keyPath] for i in p.indices { @@ -164,15 +163,15 @@ public extension NDBox { } } } - + #if DEBUG - let _box = Self(_p0, _p1) - assert(points.allSatisfy{ p in - _box.contains(p[keyPath: keyPath]) - }) + let _box = Self(_p0, _p1) + assert( + points.allSatisfy { p in + _box.contains(p[keyPath: keyPath]) + }) #endif - + return Self(_p0, _p1) } } - diff --git a/Sources/NDTree/NDTree+Traversable.swift b/Sources/NDTree/NDTree+Traversable.swift deleted file mode 100644 index ee11198..0000000 --- a/Sources/NDTree/NDTree+Traversable.swift +++ /dev/null @@ -1,46 +0,0 @@ -// -// File.swift -// -// -// Created by li3zhen1 on 10/16/23. -// - -public protocol Traversable { - - @inlinable func visit( - shouldVisitChildren: (Self) -> Bool - ) - - @inlinable func visitPostOrdered( - _ action: (Self) -> () - ) - -} - - -extension NDTree: Traversable { - - /// Visit the tree in pre-order. - /// - Parameter shouldVisitChildren: a closure that returns a boolean value indicating whether should continue to visit children. - @inlinable public func visit(shouldVisitChildren: (NDTree) -> Bool) { - if shouldVisitChildren(self), let children { - // this is an internal node - for t in children { - t.visit(shouldVisitChildren: shouldVisitChildren) - } - } - } - - /// Visit the tree in post-order. - /// - Parameter action: a closure that takes a tree as its argument. - @inlinable public func visitPostOrdered( - _ action: (NDTree) -> () - ) { - if let children { - for c in children { - c.visitPostOrdered(action) - } - } - action(self) - } -} diff --git a/Sources/NDTree/NDTree.swift b/Sources/NDTree/NDTree.swift index 306a0b1..7231b6b 100644 --- a/Sources/NDTree/NDTree.swift +++ b/Sources/NDTree/NDTree.swift @@ -321,4 +321,28 @@ extension NDTree { /// Returns true is the current tree node is leaf and does not have point in it. @inlinable public var isEmptyLeaf: Bool { nodePosition == nil } + + /// Visit the tree in pre-order. + /// - Parameter shouldVisitChildren: a closure that returns a boolean value indicating whether should continue to visit children. + @inlinable public func visit(shouldVisitChildren: (NDTree) -> Bool) { + if shouldVisitChildren(self), let children { + // this is an internal node + for t in children { + t.visit(shouldVisitChildren: shouldVisitChildren) + } + } + } + + /// Visit the tree in post-order. + /// - Parameter action: a closure that takes a tree as its argument. + @inlinable public func visitPostOrdered( + _ action: (NDTree) -> () + ) { + if let children { + for c in children { + c.visitPostOrdered(action) + } + } + action(self) + } } diff --git a/Sources/NDTree/Octree.swift b/Sources/NDTree/Octree.swift new file mode 100644 index 0000000..d22136a --- /dev/null +++ b/Sources/NDTree/Octree.swift @@ -0,0 +1,344 @@ +// +// Octree.swift +// +// +// Created by li3zhen1 on 10/22/23. +// + +#if canImport(simd) +import simd + +/// The data structure carried by a node of NDTree +/// It receives notifications when a node is added or removed on a node, regardless of whether the node is internal or leaf. +/// It is designed to calculate properties like a box's center of mass. +public protocol OctreeDelegate { + associatedtype NodeID: Hashable + typealias V = simd_float3 + + /// Called when a node is added on a node, regardless of whether the node is internal or leaf. + /// If you add `n` points to the root, this method will be called `n` times in the root delegate, + /// although it is probably not containing points now. + /// - Parameters: + /// - node: The nodeID of the node that is added. + /// - position: The position of the node that is added. + mutating func didAddNode(_ node: NodeID, at position: V) + + /// Called when a node is removed on a node, regardless of whether the node is internal or leaf. + mutating func didRemoveNode(_ node: NodeID, at position: V) + + /// Copy object. This method is called when the root box is not large enough to cover the new nodes. + /// The method + func copy() -> Self + + /// Create new object with properties set to initial value as if the box is empty. + /// However, you can still carry something like a closure to get information from outside. + /// This method is called when a leaf box is splited due to the insertion of a new node in this box. + func spawn() -> Self +} + +/// A node in NDTree +/// - Note: `NDTree` is a generic type that can be used in any dimension. +/// `NDTree` is a reference type. +public final class Octree where D: OctreeDelegate { + + public typealias V = simd_float3 + + public typealias NodeIndex = D.NodeID + + public typealias Direction = Int + + public typealias Box = NDBox + + public private(set) var box: Box + + public private(set) var children: [Octree]? + + public private(set) var nodePosition: V? + public private(set) var nodeIndices: [NodeIndex] + + public let clusterDistance: V.Scalar + private let clusterDistanceSquared: V.Scalar + + public private(set) var delegate: D + + private init( + box: Box, + clusterDistance: V.Scalar, + parentDelegate: D + ) { + self.box = box + self.clusterDistance = clusterDistance + self.clusterDistanceSquared = clusterDistance * clusterDistance + self.nodeIndices = [] + self.delegate = parentDelegate.spawn() + } + + public init( + box: Box, + clusterDistance: V.Scalar, + buildRootDelegate: () -> D + ) { + self.box = box + self.clusterDistance = clusterDistance + self.clusterDistanceSquared = clusterDistance * clusterDistance + self.nodeIndices = [] + self.delegate = buildRootDelegate() + } + + public convenience init( + covering nodes: [NodeIndex: V], + clusterDistance: V.Scalar, + buildRootDelegate: () -> D + ) { + let coveringBox = Box.cover(of: Array(nodes.values)) + self.init( + box: coveringBox, + clusterDistance: clusterDistance, + buildRootDelegate: buildRootDelegate) + for (i, p) in nodes { + add(i, at: p) + } + } + + public func add(_ nodeIndex: NodeIndex, at point: V) { + cover(point) + addWithoutCover(nodeIndex, at: point) + } + + private func addWithoutCover(_ nodeIndex: NodeIndex, at point: V) { + defer { + delegate.didAddNode(nodeIndex, at: point) + } + + guard let children = self.children else { + if nodePosition == nil { + nodeIndices.append(nodeIndex) + nodePosition = point + return + } else if nodePosition == point + || simd_length_squared(nodePosition! - point) < clusterDistanceSquared + { + nodeIndices.append(nodeIndex) + return + } else { + + let spawned = Self.spawnChildren( + box, + // Self.directionCount, + clusterDistance, + /*&*/delegate + ) + + if let nodePosition { + let direction = getIndexInChildren(nodePosition, relativeTo: box.center) + spawned[direction].nodeIndices = self.nodeIndices + spawned[direction].nodePosition = self.nodePosition + spawned[direction].delegate = self.delegate.copy() + // self.delegate = self.delegate.copy() + + // for ni in nodeIndices { + // delegate.didAddNode(ni, at: nodePosition) + // } + + self.nodeIndices = [] + self.nodePosition = nil + } + + let directionOfNewNode = getIndexInChildren(point, relativeTo: box.center) + spawned[directionOfNewNode].addWithoutCover(nodeIndex, at: point) + + self.children = spawned + return + + } + } + + let directionOfNewNode = getIndexInChildren(point, relativeTo: box.center) + children[directionOfNewNode].addWithoutCover(nodeIndex, at: point) + + return + } + + /// Expand the current node multiple times by calling `expand(towards:)`, until the point is covered. + /// - Parameter point: The point to be covered. + private func cover(_ point: V) { + if box.contains(point) { return } + + repeat { + let direction = getIndexInChildren(point, relativeTo: box.p0) + expand(towards: direction) + } while !box.contains(point) + } + + /// 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). + private func expand(towards direction: Direction) { + let nailedDirection = (Self.directionCount - 1) - direction + let nailedCorner = box.getCorner(of: nailedDirection) + + let _corner = box.getCorner(of: direction) + let expandedCorner = (_corner + _corner) - nailedCorner + + let newRootBox = Box(nailedCorner, expandedCorner) + + let copiedCurrentNode = shallowCopy() + var spawned = Self.spawnChildren( + newRootBox, + // Self.directionCount, + clusterDistance, + /*&*/delegate + ) + + spawned[nailedDirection] = copiedCurrentNode + + self.box = newRootBox + self.children = spawned + self.nodeIndices = [] + self.delegate = delegate.copy() + + } + + /// The children count of a node in NDTree. + /// Should be equal to the 2^(dimension of the vector). + /// For example, a 2D vector should have 4 children, a 3D vector should have 8 children. + /// This property is a getter property but it is probably be inlined. + @inlinable static var directionCount: Int { 1 << V.scalarCount } + + private static func spawnChildren( + _ _box: Box, + _ _clusterDistance: V.Scalar, + _ _delegate: D + ) -> [Octree] { + + var result = [Octree]() + result.reserveCapacity(Self.directionCount) + let center = _box.center + + for j in 0..> i) & 0b1 + + // TODO: use simd mask + if isOnTheHigherRange != 0 { + __box.p0[i] = center[i] + } else { + __box.p1[i] = center[i] + } + } + result.append( + Self( + box: __box, clusterDistance: _clusterDistance, parentDelegate: /*&*/ _delegate) + ) + } + + return result + } + + /// Copy object while holding the same reference to children. + /// Consider this function something you would do when working with linked list. + private func shallowCopy() -> Self { + let copy = Self( + box: box, clusterDistance: clusterDistance, parentDelegate: /*&*/ delegate) + + copy.nodeIndices = nodeIndices + copy.nodePosition = nodePosition + copy.children = children + copy.delegate = delegate + + return copy + } + + /// Get the index of the child that contains the point. + /// **Complexity**: `O(n*(2^n))`, where `n` is the dimension of the vector. + private func getIndexInChildren(_ point: V, relativeTo originalPoint: V) -> Int { + // var index = 0 + // for i in 0..= originalPoint[i] { // isOnHigherRange in this dimension + // index |= (1 << i) + // } + // } + // return index + + let mask = point .>= originalPoint + + return (mask[0] ? 1 : 0) | (mask[1] ? 2 : 0) | (mask[2] ? 4 : 0) + } + +} + +extension Octree where D.NodeID == Int { + + /// Initialize a NDTree with a list of points and a key path to the vector. + /// - Parameters: + /// - points: A list of points. The points are only used to calculate the covering box. You should still call `add` to add the points to the tree. + /// - clusterDistance: If 2 points are close enough, they will be clustered into the same leaf node. + /// - buildRootDelegate: A closure that tells the tree how to initialize the data you want to store in the root. + /// The closure is called only once. The `NDTreeDelegate` will then be created in children tree nods by calling `spawn` on the root delegate. + public convenience init( + covering points: [V], + clusterDistance: V.Scalar, + buildRootDelegate: () -> D + ) { + let coveringBox = Box.cover(of: points) + self.init( + box: coveringBox, clusterDistance: clusterDistance, buildRootDelegate: buildRootDelegate + ) + for i in points.indices { + add(i, at: points[i]) + } + } + + /// Initialize a NDTree with a list of points and a key path to the vector. + /// - Parameters: + /// - points: A list of points. The points are only used to calculate the covering box. You should still call `add` to add the points to the tree. + /// - keyPath: A key path to the vector in the element of the list. + /// - clusterDistance: If 2 points are close enough, they will be clustered into the same leaf node. + /// - buildRootDelegate: A closure that tells the tree how to initialize the data you want to store in the root. + /// The closure is called only once. The `NDTreeDelegate` will then be created in children tree nods by calling `spawn` on the root delegate. + public convenience init( + covering points: [T], + keyPath: KeyPath, + clusterDistance: V.Scalar, + buildRootDelegate: () -> D + ) { + let coveringBox = Box.cover(of: points, keyPath: keyPath) + self.init( + box: coveringBox, clusterDistance: clusterDistance, buildRootDelegate: buildRootDelegate + ) + for i in points.indices { + add(i, at: points[i][keyPath: keyPath]) + } + } +} + +extension Octree { + + /// The bounding box of the current node + @inlinable public var extent: Box { box } + + /// Returns true is the current tree node is leaf. Does not guarantee that the tree node has point in it. + @inlinable public var isLeaf: Bool { children == nil } + + /// Returns true is the current tree node is internal. Internal tree node are always empty and do not contain any points. + @inlinable public var isInternalNode: Bool { children != nil } + + /// Returns true is the current tree node is leaf and has point in it. + @inlinable public var isFilledLeaf: Bool { nodePosition != nil } + + /// Returns true is the current tree node is leaf and does not have point in it. + @inlinable public var isEmptyLeaf: Bool { nodePosition == nil } + + + public func visit(shouldVisitChildren: (Octree) -> Bool) { + if shouldVisitChildren(self), let children { + // this is an internal node + for t in children { + t.visit(shouldVisitChildren: shouldVisitChildren) + } + } + } +} + +#endif \ No newline at end of file diff --git a/Sources/NDTree/Quadtree+Octree.swift b/Sources/NDTree/Quadtree+Octree.swift deleted file mode 100644 index 1585d22..0000000 --- a/Sources/NDTree/Quadtree+Octree.swift +++ /dev/null @@ -1,102 +0,0 @@ -// -// File 2.swift -// -// -// Created by li3zhen1 on 10/13/23. -// - -#if canImport(simd) - - import simd - extension simd_double2: VectorLike { - - @inlinable public func lengthSquared() -> Double { - return x * x + y * y - } - - @inlinable public func length() -> Double { - return (x * x + y * y).squareRoot() - } - - @inlinable public func distanceSquared(to: SIMD2) -> Scalar { - return (self - to).lengthSquared() - } - - @inlinable public func distance(to: SIMD2) -> Scalar { - return (self - to).length() - } - - - // public static let directionCount = 4 - } - - extension simd_double3: VectorLike { - - @inlinable public func lengthSquared() -> Double { - return x * x + y * y + z * z - } - - @inlinable public func length() -> Double { - return (x * x + y * y + z * z).squareRoot() - } - - @inlinable public func distanceSquared(to: SIMD3) -> Scalar { - return (self - to).lengthSquared() - } - - @inlinable public func distance(to: SIMD3) -> Scalar { - return (self - to).length() - } - - // public static let directionCount = 8 - } - - - - - public typealias Vector2d = simd_double2 - public typealias Vector3d = simd_double3 - - public protocol QuadtreeDelegate: NDTreeDelegate where V == Vector2d {} - public protocol OctreeDelegate: NDTreeDelegate where V == Vector3d {} - - - public typealias QuadBox = NDBox - public typealias OctBox = NDBox - - public typealias Quadtree = NDTree - public typealias Octree = NDTree - - - -/// Uncomment the region below to unlock 4d tree - -// extension simd_double4: VectorLike { -// -// @inlinable public func lengthSquared() -> Double { -// return x * x + y * y + z * z + w * w -// } -// -// @inlinable public func length() -> Double { -// return (x * x + y * y + z * z + w * w).squareRoot() -// } -// -// @inlinable public func distanceSquared(to: SIMD4) -> Scalar { -// return (self - to).lengthSquared() -// } -// -// @inlinable public func distance(to: SIMD4) -> Scalar { -// return (self - to).length() -// } -// public static let directionCount = 16 -// } -// public typealias Vector4d = simd_double4 -// public protocol HyperoctreeDelegate: NDTreeDelegate where V == Vector4d {} -// public typealias HyperoctBox = NDBox -// public typealias Hyperoctree = NDTree - - - - - -#endif diff --git a/Sources/NDTree/Quadtree.swift b/Sources/NDTree/Quadtree.swift new file mode 100644 index 0000000..2b05006 --- /dev/null +++ b/Sources/NDTree/Quadtree.swift @@ -0,0 +1,344 @@ +// +// Quadtree.swift +// +// +// Created by li3zhen1 on 10/22/23. +// + +#if canImport(simd) +import simd + +/// The data structure carried by a node of NDTree +/// It receives notifications when a node is added or removed on a node, regardless of whether the node is internal or leaf. +/// It is designed to calculate properties like a box's center of mass. +public protocol QuadtreeDelegate { + associatedtype NodeID: Hashable + typealias V = simd_double2 + + /// Called when a node is added on a node, regardless of whether the node is internal or leaf. + /// If you add `n` points to the root, this method will be called `n` times in the root delegate, + /// although it is probably not containing points now. + /// - Parameters: + /// - node: The nodeID of the node that is added. + /// - position: The position of the node that is added. + mutating func didAddNode(_ node: NodeID, at position: V) + + /// Called when a node is removed on a node, regardless of whether the node is internal or leaf. + mutating func didRemoveNode(_ node: NodeID, at position: V) + + /// Copy object. This method is called when the root box is not large enough to cover the new nodes. + /// The method + func copy() -> Self + + /// Create new object with properties set to initial value as if the box is empty. + /// However, you can still carry something like a closure to get information from outside. + /// This method is called when a leaf box is splited due to the insertion of a new node in this box. + func spawn() -> Self +} + +/// A node in NDTree +/// - Note: `NDTree` is a generic type that can be used in any dimension. +/// `NDTree` is a reference type. +public final class Quadtree where D: QuadtreeDelegate { + + public typealias V = simd_double2 + + public typealias NodeIndex = D.NodeID + + public typealias Direction = Int + + public typealias Box = NDBox + + public private(set) var box: Box + + public private(set) var children: [Quadtree]? + + public private(set) var nodePosition: V? + public private(set) var nodeIndices: [NodeIndex] + + public let clusterDistance: V.Scalar + private let clusterDistanceSquared: V.Scalar + + public private(set) var delegate: D + + private init( + box: Box, + clusterDistance: V.Scalar, + parentDelegate: D + ) { + self.box = box + self.clusterDistance = clusterDistance + self.clusterDistanceSquared = clusterDistance * clusterDistance + self.nodeIndices = [] + self.delegate = parentDelegate.spawn() + } + + public init( + box: Box, + clusterDistance: V.Scalar, + buildRootDelegate: () -> D + ) { + self.box = box + self.clusterDistance = clusterDistance + self.clusterDistanceSquared = clusterDistance * clusterDistance + self.nodeIndices = [] + self.delegate = buildRootDelegate() + } + + public convenience init( + covering nodes: [NodeIndex: V], + clusterDistance: V.Scalar, + buildRootDelegate: () -> D + ) { + let coveringBox = Box.cover(of: Array(nodes.values)) + self.init( + box: coveringBox, + clusterDistance: clusterDistance, + buildRootDelegate: buildRootDelegate) + for (i, p) in nodes { + add(i, at: p) + } + } + + public func add(_ nodeIndex: NodeIndex, at point: V) { + cover(point) + addWithoutCover(nodeIndex, at: point) + } + + private func addWithoutCover(_ nodeIndex: NodeIndex, at point: V) { + defer { + delegate.didAddNode(nodeIndex, at: point) + } + + guard let children = self.children else { + if nodePosition == nil { + nodeIndices.append(nodeIndex) + nodePosition = point + return + } else if nodePosition == point + || simd_length_squared(nodePosition! - point) < clusterDistanceSquared + { + nodeIndices.append(nodeIndex) + return + } else { + + let spawned = Self.spawnChildren( + box, + // Self.directionCount, + clusterDistance, + /*&*/delegate + ) + + if let nodePosition { + let direction = getIndexInChildren(nodePosition, relativeTo: box.center) + spawned[direction].nodeIndices = self.nodeIndices + spawned[direction].nodePosition = self.nodePosition + spawned[direction].delegate = self.delegate.copy() + // self.delegate = self.delegate.copy() + + // for ni in nodeIndices { + // delegate.didAddNode(ni, at: nodePosition) + // } + + self.nodeIndices = [] + self.nodePosition = nil + } + + let directionOfNewNode = getIndexInChildren(point, relativeTo: box.center) + spawned[directionOfNewNode].addWithoutCover(nodeIndex, at: point) + + self.children = spawned + return + + } + } + + let directionOfNewNode = getIndexInChildren(point, relativeTo: box.center) + children[directionOfNewNode].addWithoutCover(nodeIndex, at: point) + + return + } + + /// Expand the current node multiple times by calling `expand(towards:)`, until the point is covered. + /// - Parameter point: The point to be covered. + private func cover(_ point: V) { + if box.contains(point) { return } + + repeat { + let direction = getIndexInChildren(point, relativeTo: box.p0) + expand(towards: direction) + } while !box.contains(point) + } + + /// 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). + private func expand(towards direction: Direction) { + let nailedDirection = (Self.directionCount - 1) - direction + let nailedCorner = box.getCorner(of: nailedDirection) + + let _corner = box.getCorner(of: direction) + let expandedCorner = (_corner + _corner) - nailedCorner + + let newRootBox = Box(nailedCorner, expandedCorner) + + let copiedCurrentNode = shallowCopy() + var spawned = Self.spawnChildren( + newRootBox, + // Self.directionCount, + clusterDistance, + /*&*/delegate + ) + + spawned[nailedDirection] = copiedCurrentNode + + self.box = newRootBox + self.children = spawned + self.nodeIndices = [] + self.delegate = delegate.copy() + + } + + /// The children count of a node in NDTree. + /// Should be equal to the 2^(dimension of the vector). + /// For example, a 2D vector should have 4 children, a 3D vector should have 8 children. + /// This property is a getter property but it is probably be inlined. + @inlinable static var directionCount: Int { 1 << V.scalarCount } + + private static func spawnChildren( + _ _box: Box, + _ _clusterDistance: V.Scalar, + _ _delegate: D + ) -> [Quadtree] { + + var result = [Quadtree]() + result.reserveCapacity(Self.directionCount) + let center = _box.center + + for j in 0..> i) & 0b1 + + // TODO: use simd mask + if isOnTheHigherRange != 0 { + __box.p0[i] = center[i] + } else { + __box.p1[i] = center[i] + } + } + result.append( + Self( + box: __box, clusterDistance: _clusterDistance, parentDelegate: /*&*/ _delegate) + ) + } + + return result + } + + /// Copy object while holding the same reference to children. + /// Consider this function something you would do when working with linked list. + private func shallowCopy() -> Self { + let copy = Self( + box: box, clusterDistance: clusterDistance, parentDelegate: /*&*/ delegate) + + copy.nodeIndices = nodeIndices + copy.nodePosition = nodePosition + copy.children = children + copy.delegate = delegate + + return copy + } + + /// Get the index of the child that contains the point. + /// **Complexity**: `O(n*(2^n))`, where `n` is the dimension of the vector. + private func getIndexInChildren(_ point: V, relativeTo originalPoint: V) -> Int { + // var index = 0 + // for i in 0..= originalPoint[i] { // isOnHigherRange in this dimension + // index |= (1 << i) + // } + // } + // return index + + let mask = point .>= originalPoint + + return (mask[0] ? 1 : 0) | (mask[1] ? 2 : 0) + } + +} + +extension Quadtree where D.NodeID == Int { + + /// Initialize a NDTree with a list of points and a key path to the vector. + /// - Parameters: + /// - points: A list of points. The points are only used to calculate the covering box. You should still call `add` to add the points to the tree. + /// - clusterDistance: If 2 points are close enough, they will be clustered into the same leaf node. + /// - buildRootDelegate: A closure that tells the tree how to initialize the data you want to store in the root. + /// The closure is called only once. The `NDTreeDelegate` will then be created in children tree nods by calling `spawn` on the root delegate. + public convenience init( + covering points: [V], + clusterDistance: V.Scalar, + buildRootDelegate: () -> D + ) { + let coveringBox = Box.cover(of: points) + self.init( + box: coveringBox, clusterDistance: clusterDistance, buildRootDelegate: buildRootDelegate + ) + for i in points.indices { + add(i, at: points[i]) + } + } + + /// Initialize a NDTree with a list of points and a key path to the vector. + /// - Parameters: + /// - points: A list of points. The points are only used to calculate the covering box. You should still call `add` to add the points to the tree. + /// - keyPath: A key path to the vector in the element of the list. + /// - clusterDistance: If 2 points are close enough, they will be clustered into the same leaf node. + /// - buildRootDelegate: A closure that tells the tree how to initialize the data you want to store in the root. + /// The closure is called only once. The `NDTreeDelegate` will then be created in children tree nods by calling `spawn` on the root delegate. + public convenience init( + covering points: [T], + keyPath: KeyPath, + clusterDistance: V.Scalar, + buildRootDelegate: () -> D + ) { + let coveringBox = Box.cover(of: points, keyPath: keyPath) + self.init( + box: coveringBox, clusterDistance: clusterDistance, buildRootDelegate: buildRootDelegate + ) + for i in points.indices { + add(i, at: points[i][keyPath: keyPath]) + } + } +} + +extension Quadtree { + + /// The bounding box of the current node + @inlinable public var extent: Box { box } + + /// Returns true is the current tree node is leaf. Does not guarantee that the tree node has point in it. + @inlinable public var isLeaf: Bool { children == nil } + + /// Returns true is the current tree node is internal. Internal tree node are always empty and do not contain any points. + @inlinable public var isInternalNode: Bool { children != nil } + + /// Returns true is the current tree node is leaf and has point in it. + @inlinable public var isFilledLeaf: Bool { nodePosition != nil } + + /// Returns true is the current tree node is leaf and does not have point in it. + @inlinable public var isEmptyLeaf: Bool { nodePosition == nil } + + + public func visit(shouldVisitChildren: (Quadtree) -> Bool) { + if shouldVisitChildren(self), let children { + // this is an internal node + for t in children { + t.visit(shouldVisitChildren: shouldVisitChildren) + } + } + } +} + +#endif \ No newline at end of file diff --git a/Sources/NDTree/VectorLike.swift b/Sources/NDTree/VectorLike.swift index be3a28b..4bc3c42 100644 --- a/Sources/NDTree/VectorLike.swift +++ b/Sources/NDTree/VectorLike.swift @@ -33,8 +33,8 @@ public protocol VectorLike: CustomStringConvertible, Decodable, Encodable, Expre @inlinable func distance(to: Self) -> Scalar - @inlinable static func * (a: Self, b: Double) -> Self - @inlinable static func / (a: Self, b: Double) -> Self + // @inlinable static func * (a: Self, b: Scalar) -> Self + // @inlinable static func / (a: Self, b: Scalar) -> Self @inlinable static func * (a: Self, b: Scalar) -> Self @inlinable static func / (a: Self, b: Scalar) -> Self @@ -57,31 +57,5 @@ public protocol VectorLike: CustomStringConvertible, Decodable, Encodable, Expre subscript(index: Int) -> Self.Scalar { get set } var indices: Range { get } - - - - -// mutating func replace(with other: Self, where mask: M) where M:MaskLike, M.Storage.Scalar==Scalar.SIMDMaskScalar -// public static func .< (a: SIMD32, b: SIMD32) -> SIMDMask> -// -// /// Returns a vector mask with the result of a pointwise less than or equal -// /// comparison. -// public static func .<= (a: SIMD32, b: SIMD32) -> SIMDMask> -// -// /// The least element in the vector. -// public func min() -> Scalar -// -// /// The greatest element in the vector. -// public func max() -> Scalar -// -// /// Returns a vector mask with the result of a pointwise greater than or -// /// equal comparison. -// public static func .>= (a: SIMD32, b: SIMD32) -> SIMDMask> -// -// /// Returns a vector mask with the result of a pointwise greater than -// /// comparison. -// public static func .> (a: SIMD32, b: SIMD32) -> SIMDMask> - - } diff --git a/Sources/NDTree/simd+VectorLike.swift b/Sources/NDTree/simd+VectorLike.swift new file mode 100644 index 0000000..b5a41a4 --- /dev/null +++ b/Sources/NDTree/simd+VectorLike.swift @@ -0,0 +1,54 @@ +// +// simd+VectorLike.swift +// +// +// Created by li3zhen1 on 10/13/23. +// + +#if canImport(simd) + +import simd + +extension simd_double2: VectorLike { + @inlinable public func lengthSquared() -> Scalar { + return simd_length_squared(self) + } + + @inlinable public func length() -> Scalar { + return simd_length(self) + } + + @inlinable public func distanceSquared(to: SIMD2) -> Scalar { + return simd_length_squared(self - to) + } + + @inlinable public func distance(to: SIMD2) -> Scalar { + return simd_length(self - to) + } + +} + +extension simd_float3: VectorLike { + + @inlinable public func lengthSquared() -> Scalar { + return simd_length_squared(self) + } + + @inlinable public func length() -> Scalar { + return simd_length(self) + } + + @inlinable public func distanceSquared(to: SIMD3) -> Scalar { + return simd_length_squared(self - to) + } + + @inlinable public func distance(to: SIMD3) -> Scalar { + return simd_length(self - to) + } + +} + +public typealias QuadBox = NDBox +public typealias OctBox = NDBox + +#endif \ No newline at end of file diff --git a/Tests/ForceSimulationTests/ManyBodyForceTest.swift b/Tests/ForceSimulationTests/ManyBodyForceTest.swift index 103934e..048dfe4 100644 --- a/Tests/ForceSimulationTests/ManyBodyForceTest.swift +++ b/Tests/ForceSimulationTests/ManyBodyForceTest.swift @@ -31,27 +31,6 @@ final class ManyBodyForceTests: XCTestCase { .make("David") ] -// let pos = [(-1,-1), (-1,1), (1,-1), (1, 1)] -// let sim = Simulation(nodes: nodes) { n, i in -// n.position = Vector2f(x: Double(pos[i].0), y: Double(pos[i].1)) -// } -// -// -// let f = sim.createManyBodyForce(strength: 0.4) -// -// sim.tick() - - - - // sim.simulationNodes.forEach { n in - // print(n) - // } - - - // sim.tick() - // sim.simulationNodes.forEach { n in - // print(n) - // } } diff --git a/Tests/ForceSimulationTests/MiserableGraphTest.swift b/Tests/ForceSimulationTests/MiserableGraphTest.swift index b1b5fdc..968af77 100644 --- a/Tests/ForceSimulationTests/MiserableGraphTest.swift +++ b/Tests/ForceSimulationTests/MiserableGraphTest.swift @@ -6,16 +6,17 @@ // import NDTree +import simd import XCTest @testable import ForceSimulation final class MiserableGraphTest: XCTestCase { - func test() { + func test_Generic_2d() { let data = getData() - let sim = Simulation( + let sim = SimulationKD( nodeIds: data.nodes.map { n in n.id }) @@ -40,8 +41,48 @@ final class MiserableGraphTest: XCTestCase { // } measure { - for _ in 0..<120 { + for i in 0..<120 { + sim.tick() +// print(i) + } + } + sim.tick() + // print(sim.simulationNodes) + + } + + + func test_Inlined_2d() { + let data = getData() + + let sim = Simulation2D( + nodeIds: data.nodes.map { n in + n.id + }) + + let linkForce = sim.createLinkForce( + data.links.map({ l in + (l.source, l.target) + })) + let manybodyForce = sim.createManyBodyForce(strength: -30) + + let centerForce = sim.createCenterForce(center: .zero) + let collideForce = sim.createCollideForce(radius: .constant(5)) + + // for _ in 0..<120{ + // sim.tick() + // } + // + //// sim.tick() + // + // for _ in 0..<120{ + // sim.tick() + // } + + measure { + for i in 0..<120 { sim.tick() +// print(i) } } sim.tick() @@ -53,7 +94,7 @@ final class MiserableGraphTest: XCTestCase { func test3d() { let data = getData() - let sim = Simulation( + let sim = Simulation3D( nodeIds: data.nodes.map { n in n.id }) diff --git a/Tests/NDTreeTests/NDTree+CustomDebugStringConvertible.swift b/Tests/NDTreeTests/NDTree+CustomDebugStringConvertible.swift index e21a9ca..0783e94 100644 --- a/Tests/NDTreeTests/NDTree+CustomDebugStringConvertible.swift +++ b/Tests/NDTreeTests/NDTree+CustomDebugStringConvertible.swift @@ -5,6 +5,7 @@ // Created by li3zhen1 on 10/15/23. // +import simd @testable import NDTree extension VectorLike { @@ -51,3 +52,43 @@ extension NDTree: CustomDebugStringConvertible where D.NodeID == Int { } } } + + +extension simd_double2 { + var compactDebugString: String { + return "[\(self.x), \(self.y)]" + } +} + + + +extension Quadtree: CustomDebugStringConvertible where D.NodeID == Int { + public var debugDescription: String { + if let children { + return "[\(children.map{$0.debugDescription}.joined(separator: ","))]" + } else { + if nodeIndices.count == 0 { + return "" + } else if nodeIndices.count == 1 { + let p = nodePosition! + return "{data: \(p.compactDebugString)}" + } else { + let p = nodePosition! + + var r1 = "" + var r2 = "" + + for i in nodeIndices { + if i == nodeIndices.count - 1 { + r1 += "{data: \(p)" + } else { + r1 += "{data: \(p), next:" + } + r2 += "}" + } + + return r1 + r2 + } + } + } +} diff --git a/Tests/NDTreeTests/QuadtreeTests.swift b/Tests/NDTreeTests/QuadtreeTests.swift index 9554753..d7328e8 100644 --- a/Tests/NDTreeTests/QuadtreeTests.swift +++ b/Tests/NDTreeTests/QuadtreeTests.swift @@ -32,7 +32,7 @@ final class AddTests: XCTestCase { // tests below are mainly generate by github copilot with d3 source code private func t( - _ points: [Vector2d] + _ points: [simd_double2] ) -> Quadtree { // var del = DummyQuadtreeDelegate() let box = QuadBox.cover(of: points) @@ -47,7 +47,7 @@ final class AddTests: XCTestCase { private func t( _ box: QuadBox, - _ points: [Vector2d] + _ points: [simd_double2] ) -> Quadtree { // var del = DummyQuadtreeDelegate() let qt = Quadtree(box: box, clusterDistance: 1e-5) {