Add PositionForce

README.md

@@ -20,14 +20,17 @@ https://github.com/li3zhen1/Grape/assets/45376537/5f57c223-0126-428a-a72d-d9a3ed
 
 #### Features
 
-| Feature | Status |
-| --- | --- |
-| LinkForce | ✅ |
-| ManyBodyForce | ✅ |
-| CenterForce | ✅ |
-| CollideForce | ✅ |
-| PositionForce |  |
-| RadialForce | ✅ |
+|   | Status (2D) | Status (3D) | Metal |
+| --- | --- | --- | --- |
+| **NdTree** | ✅ | 🚧 |  |
+| **Simulation** | ✅ | 🚧 | 🚧 |
+|  LinkForce | ✅ |   |  |
+|  ManyBodyForce | ✅ |  |  |
+|  CenterForce | ✅ |  |  |
+|  CollideForce | ✅ |  |  |
+|  PositionForce | ✅ |  |  |
+|  RadialForce | ✅ |  |  |
+| **SwiftUI View** | 🚧 |  |  |
 
 
 #### Usage

Sources/ForceSimulation/forces/CollideForce.swift

@@ -118,14 +118,13 @@ extension CollideForce: Force {
                     nodes: sim.simulationNodes.map { ($0, $0.position) },
                     getQuadDelegate: {
                         MaxRadiusQuadTreeDelegate {
-                            // switch self.radius {
-                            // case .constant(let r):
-                            //     return r
-                            // case .varied(_):
-                            //     return self.calculatedRadius[$0, default: 0.0]
-                            // }
-                            return self.calculatedRadius[$0, default: 0.0]
-                            // return self.calculatedRadius[$0]!
+                            switch self.radius {
+                            case .constant(let r):
+                                return r
+                            case .varied(_):
+                                return self.calculatedRadius[$0, default: 0.0]
+                            }
+                            // return self.calculatedRadius[$0, default: 0.0]
                         }
                     }
                 )

Sources/ForceSimulation/forces/LinkForce.swift

@@ -76,7 +76,6 @@ final public class LinkForce<N>: Force where N: Identifiable {
     var calculatedBias: [Float] = []
     weak var simulation: Simulation<N>?
 
-
     var iterations: Int
 
     internal init(

Sources/ForceSimulation/forces/ManyBodyForce.swift

@@ -195,9 +195,7 @@ final public class ManyBodyForce<N>: Force where N: Identifiable {
         let quad = try QuadTree2(
             nodes: sim.simulationNodes.map { ($0, $0.position) }
         ) {
-            // this switch is only called on root init
-            // but it significantly slows down the performance
-            //
+            // this switch is only called on root init 
             return switch self.mass {
             case .constant(let m):
                 MassQuadTreeDelegate<SimulationNode<N.ID>> { _ in m }

Sources/ForceSimulation/forces/PositionForce.swift

@@ -5,22 +5,88 @@
 //  Created by li3zhen1 on 10/1/23.
 //
 
-
 final public class PositionForce<N>: Force where N: Identifiable {
-    var x: Float
-    var y: Float
-    var strength: Float
 
-    init(x: Float, y: Float, strength: Float = 0.1) {
-        self.x = x
-        self.y = y
+    public enum Direction {
+        case x
+        case y
+    }
+    public enum TargetOnDirection {
+        case constant(Float)
+        case varied([N.ID: Float])
+    }
+    public enum Strength {
+        case constant(Float)
+        case varied([N.ID: Float])
+    }
+    public var strength: Strength
+    public var direction: Direction
+    public var calculatedStrength: [N.ID: Float] = [:]
+    public var targetOnDirection: TargetOnDirection
+    public var calculatedTargetOnDirection: [N.ID: Float] = [:]
+
+    internal init(direction: Direction, targetOnDirection: TargetOnDirection, strength: Strength = .constant(1.0)) {
         self.strength = strength
+        self.direction = direction
+        self.targetOnDirection = targetOnDirection
     }
 
-    weak var simulation: Simulation<N>?
+    weak var simulation: Simulation<N>? {
+        didSet {
+            guard let sim = self.simulation else { return }
+            self.calculatedStrength = strength.calculated(sim.simulationNodes)
+            self.calculatedTargetOnDirection = targetOnDirection.calculated(sim.simulationNodes)
+        }
+    }
 
     public func apply(alpha: Float) {
+        guard let sim = self.simulation else { return }
+        let vectorIndex = self.direction == .x ? 0 : 1
+        for i in sim.simulationNodes.indices {
+            let nodeId = sim.simulationNodes[i].id
+            sim.simulationNodes[i].velocity += (
+                self.calculatedTargetOnDirection[nodeId, default: 0.0] - sim.simulationNodes[i].position[vectorIndex]
+            ) * self.calculatedStrength[nodeId, default: 0.0] * alpha
+        }
+    }
+}
 
+extension PositionForce.Strength {
+    func calculated<SimNode>(_ nodes: [SimNode]) -> [N.ID: Float] where SimNode: Identifiable, SimNode.ID == N.ID {
+        switch self {
+        case .constant(let value):
+            return nodes.reduce(into: [:]) { $0[$1.id] = value }
+        case .varied(let dict):
+            return dict
+        }
     }
+}
 
+extension PositionForce.TargetOnDirection {
+    func calculated<SimNode>(_ nodes: [SimNode]) -> [N.ID: Float] where SimNode: Identifiable, SimNode.ID == N.ID {
+        switch self {
+        case .constant(let value):
+            return nodes.reduce(into: [:]) { $0[$1.id] = value }
+        case .varied(let dict):
+            return dict
+        }
+    }
 }
+
+
+public extension Simulation {
+    func createPositionForce(
+        direction: PositionForce<N>.Direction,
+        targetOnDirection: PositionForce<N>.TargetOnDirection,
+        strength: PositionForce<N>.Strength = .constant(1.0)
+    ) -> PositionForce<N> {
+        let force = PositionForce<N>(
+            direction: direction,
+            targetOnDirection: targetOnDirection,
+            strength: strength
+        )
+        force.simulation = self
+        self.forces.append(force)
+        return force
+    }
+}

Sources/ForceSimulation/mpsforces/CenterForce.metal

@@ -0,0 +1,25 @@
+#include <metal_stdlib>
+using namespace metal;
+
+struct Node {
+    float2 position;
+    float2 velocity;
+    float2 fixation;
+};
+
+kernel void applyCenterForce(
+    device Node* nodes [[ buffer(0) ]],
+    constant float2& center [[ buffer(1) ]],
+    constant float& strength [[ buffer(2) ]],
+    uint id [[ thread_position_in_grid ]],
+    uint nodeCount [[ threads_per_grid ]])
+{
+    float2 meanPosition = float2(0.0, 0.0);
+    for (int i = 0; i < nodeCount; ++i) {
+        meanPosition += nodes[i].position;
+    }
+    meanPosition /= float(nodeCount);
+
+    float2 delta = (meanPosition - center) * strength;
+    nodes[id].position -= delta;
+}

Sources/ForceSimulation/mpsforces/MPSSimulation.swift

@@ -0,0 +1,28 @@
+import Foundation
+import Metal
+import simd
+import MetalPerformanceShaders
+
+
+final public class MPSSimulation {
+    init() {
+        guard let device = MTLCreateSystemDefaultDevice() else {
+            fatalError("")
+        }
+
+        guard let commandQueue = device.makeCommandQueue() else {
+            fatalError("")
+        }
+
+        let library = device.makeDefaultLibrary()
+
+        let function = library?.makeFunction(name: "")
+
+        var pipelineState: MTLComputePipelineState
+        do {
+            pipelineState = try device.makeComputePipelineState(function: function!)
+        } catch {
+            fatalError("")
+        }
+    }
+}

Sources/QuadTree/NdTree.swift

@@ -0,0 +1,158 @@
+//
+//  File.swift
+//  
+//
+//  Created by li3zhen1 on 10/10/23.
+//
+
+import Foundation
+
+public protocol ComponentComparable {
+    @inlinable static func <(lhs: Self, rhs: Self) -> Bool
+    @inlinable static func <=(lhs: Self, rhs: Self) -> Bool
+    @inlinable static func >(lhs: Self, rhs: Self) -> Bool
+    @inlinable static func >=(lhs: Self, rhs: Self) -> Bool
+}
+
+public struct NdBox<Coordinate> where Coordinate: SIMD<Float>, Coordinate: ComponentComparable {
+    public var p0: Coordinate
+    public var p1: Coordinate
+
+    @inlinable public init(p0:Coordinate, p1:Coordinate) {
+        var p0 = p0
+        var p1 = p1
+        for i in p0.indices {
+            if p1[i] < p0[i] {
+                swap(&p0[i], &p1[i])
+            }
+        }
+        self.p0 = p0
+        self.p1 = p1
+    }
+
+    public static var unit: Self { .init(p0: .zero, p1: .one) }
+}
+
+extension NdBox {
+    var area: Float {
+        var result: Float = 1
+        let delta = p1 - p0
+        for i in delta.indices {
+            result *= delta[i]
+        }
+        return result
+    }
+
+    var vec: Coordinate { p1 - p0 }
+
+    var center: Coordinate { (p1+p0) / 2 }
+
+    @inlinable public func getCorner<Direction>(
+        of direction: Direction
+    ) -> Coordinate where Direction: NdDirection, Direction.Coordinate == Coordinate {
+        var result = Coordinate()
+        var value = direction.rawValue
+
+        // starting from the last element
+        for i in (0..<Coordinate.scalarCount).reversed() {
+            result[i] = value % 2 == 1 ? p1[i] : p0[i]
+            value /= 2
+        }
+        return result
+    }
+
+    @inlinable func contains(_ point: Coordinate) -> Bool {
+        return (p0 <= point) && (point < p1)
+    }
+}
+
+/// Reversed bit representation for nth dimension
+/// e.g.    for 3d:         0b001 => (x:0, y:0, z:1)
+public protocol NdDirection: RawRepresentable {
+    associatedtype Coordinate: SIMD<Float>
+    var rawValue: Int { get }
+    var reversed: Self { get }
+    static var entryCount: Int { get }
+}
+
+//public extension SIMD<Float> {
+//    @inlinable func direction<Direction>(originalPoint point: Self) -> Direction where Direction: NdDirection, Direction.Coordinate == Self {
+//        
+//    }
+//}
+
+struct OctDirection: NdDirection {
+    typealias Coordinate = SIMD3<Float>
+    let rawValue: Int
+    var reversed: OctDirection {
+        return OctDirection(rawValue: 7-rawValue)
+    }
+    static let entryCount: Int = 8
+}
+
+public struct NdChildren<T, Coordinate> where Coordinate: SIMD<Float> {
+    @usableFromInline internal let children: [T]
+
+    @inlinable public subscript<Direction>(
+        at direction: Direction
+    ) -> T where Direction: NdDirection, Direction.Coordinate == Coordinate {
+        return children[direction.rawValue]
+    }
+}
+
+
+public protocol TreeNodeDelegate {
+    associatedtype Node
+    associatedtype Index: Hashable
+    mutating func didAddNode(_ node: Index, at position: Vector2f)
+    mutating func didRemoveNode(_ node: Index, at position: Vector2f)
+    func copy() -> Self
+    func createNew() -> Self
+}
+
+
+public final class NdTree<C, TD> where C:SIMD<Float>, C:ComponentComparable, TD: TreeNodeDelegate {
+
+    public typealias Box = NdBox<C>
+    public typealias Children = NdChildren<NdTree<C, TD>, C>
+    public typealias Index = Int
+    public typealias Direction = Int
+
+    public private(set) var box: NdBox<C>
+
+    public var nodes: [Index] = []
+    public private(set) var children: Children?
+    public let clusterDistance: Float
+    public var delegate: TD
+
+    init(
+        box: Box,
+        clusterDistance: Float,
+        rootNodeDelegate: TD
+    ) {
+        self.box = box
+        self.clusterDistance = clusterDistance
+        self.delegate = rootNodeDelegate.createNew()
+    }
+
+    public func add(_ nodeIndex: Index, at point: C) {
+
+    }
+
+
+    private func cover(_ point: C) {
+        if box.contains(point) { return }
+
+        repeat {
+
+//            let direction: Direction =
+//
+//            expand(towards: direction)
+
+        } while !box.contains(point)
+    }
+
+    private func expand(towards direction: Direction) {
+
+    }
+}