📄 ui/src/Graph/PhysicsLayout.cs
using System;
using System.Collections.Generic;
using System.Linq;
using System.Threading;
using Avalonia;

namespace Cfa.Ui.Graph;

public sealed record PhysicsNode(int Parent, double Radius);

public sealed record PhysicsOptions
{
    // Tick-based limits keep the same input reproducible across machines.
    public int MaxIterations { get; init; } = 600;
    public double SpringStrength { get; init; } = 0.025;
    public double Damping { get; init; } = 0.7;
    public double MaxSpeed { get; init; } = 20;
    public double Tolerance { get; init; } = 0.05;
    public double EdgeCrossingStrength { get; init; } = 0.12;
    public double EdgeNodeStrength { get; init; } = 0.15;
    public double EdgeClearance { get; init; } = 40;
}

/// <summary>Occurrence-tree physics. A node's descendants may enter its radius; all other nodes may not.</summary>
public sealed class PhysicsLayout
{
    private readonly PhysicsNode[] nodes;
    private readonly PhysicsOptions options;
    private readonly Point[] positions;
    private readonly Vector[] velocities;
    private readonly int[] enter;
    private readonly int[] exit;
    private readonly double cellSize;
    public IReadOnlyList<Point> Positions { get; }
    public int Iterations { get; private set; }
    public bool IsStable { get; private set; }
    public double MaximumOverlap { get; private set; }

    public PhysicsLayout(
        IReadOnlyList<PhysicsNode> nodes,
        PhysicsOptions? options = null,
        IReadOnlyList<Point>? initialPositions = null
    )
    {
        ArgumentNullException.ThrowIfNull(nodes);
        this.options = options ?? new PhysicsOptions();
        if (
            this.options.MaxIterations < 1
            || !double.IsFinite(this.options.SpringStrength)
            || this.options.SpringStrength <= 0
            || !double.IsFinite(this.options.Damping)
            || this.options.Damping < 0
            || this.options.Damping >= 1
            || !double.IsFinite(this.options.MaxSpeed)
            || this.options.MaxSpeed <= 0
            || !double.IsFinite(this.options.Tolerance)
            || this.options.Tolerance <= 0
            || !double.IsFinite(this.options.EdgeCrossingStrength)
            || this.options.EdgeCrossingStrength < 0
            || !double.IsFinite(this.options.EdgeNodeStrength)
            || this.options.EdgeNodeStrength < 0
            || !double.IsFinite(this.options.EdgeClearance)
            || this.options.EdgeClearance <= 0
        )
            throw new ArgumentOutOfRangeException(nameof(options));
        this.nodes = nodes.ToArray();
        positions = new Point[nodes.Count];
        velocities = new Vector[nodes.Count];
        enter = new int[nodes.Count];
        exit = new int[nodes.Count];
        var children = Enumerable.Range(0, nodes.Count).Select(_ => new List<int>()).ToArray();
        for (var i = 0; i < nodes.Count; i++)
        {
            var node = nodes[i];
            if (
                node.Parent < -1
                || node.Parent >= i
                || (i > 0 && node.Parent < 0)
                || !double.IsFinite(node.Radius)
                || node.Radius <= 0
            )
                throw new ArgumentException(
                    "Nodes must form one parent-before-child tree with finite positive radii.",
                    nameof(nodes)
                );
            if (i > 0)
                children[node.Parent].Add(i);
        }
        var clock = 0;
        var traversal = new Stack<(int Id, bool Leaving)>();
        if (nodes.Count > 0)
            traversal.Push((0, false));
        while (traversal.TryPop(out var item))
        {
            if (item.Leaving)
                exit[item.Id] = clock;
            else
            {
                enter[item.Id] = clock++;
                traversal.Push((item.Id, true));
                for (var i = children[item.Id].Count - 1; i >= 0; i--)
                    traversal.Push((children[item.Id][i], false));
            }
        }
        if (initialPositions is not null)
        {
            if (
                initialPositions.Count != nodes.Count
                || initialPositions.Any(point => !double.IsFinite(point.X) || !double.IsFinite(point.Y))
            )
                throw new ArgumentException("Supply one finite position per node.", nameof(initialPositions));
            for (var i = 0; i < nodes.Count; i++)
                positions[i] = initialPositions[i] - initialPositions[0];
        }
        else
        {
            for (var i = 0; i < nodes.Count; i++)
            {
                var heading =
                    i == 0
                        ? 0
                        : Math.Atan2(
                            positions[i].Y - positions[nodes[i].Parent].Y,
                            positions[i].X - positions[nodes[i].Parent].X
                        );
                for (var sibling = 0; sibling < children[i].Count; sibling++)
                {
                    var child = children[i][sibling];
                    var angle =
                        heading
                        + (
                            i == 0
                                ? sibling * Math.PI * 2 / children[i].Count
                                : (sibling - (children[i].Count - 1) / 2d) * Math.PI / Math.Max(1, children[i].Count)
                        );
                    positions[child] =
                        positions[i] + new Vector(Math.Cos(angle), Math.Sin(angle)) * MinimumDistance(i, child);
                }
            }
        }
        cellSize = nodes.Count > 0 ? nodes.Max(node => node.Radius) : 1;
        Positions = Array.AsReadOnly(positions);
    }

    public double MinimumDistance(int first, int second)
    {
        ArgumentOutOfRangeException.ThrowIfNegative(first);
        ArgumentOutOfRangeException.ThrowIfGreaterThanOrEqual(first, nodes.Length);
        ArgumentOutOfRangeException.ThrowIfNegative(second);
        ArgumentOutOfRangeException.ThrowIfGreaterThanOrEqual(second, nodes.Length);
        if (first == second)
            return 0;
        var firstRadius = IsDescendant(first, second) ? 0 : nodes[first].Radius;
        var secondRadius = IsDescendant(second, first) ? 0 : nodes[second].Radius;
        return Math.Max(firstRadius, secondRadius);
    }

    private bool IsDescendant(int ancestor, int child) =>
        enter[ancestor] < enter[child] && enter[child] < exit[ancestor];

    public void Settle(CancellationToken cancellationToken = default)
    {
        var quietTicks = 0;
        while (Iterations < options.MaxIterations && quietTicks < 8)
        {
            var movement = Tick(cancellationToken);
            quietTicks = movement < options.Tolerance && MaximumOverlap < options.Tolerance ? quietTicks + 1 : 0;
        }
        IsStable = quietTicks >= 8;
    }

    public double Tick(CancellationToken cancellationToken = default)
    {
        cancellationToken.ThrowIfCancellationRequested();
        var before = (Point[])positions.Clone();
        var forces = new Vector[nodes.Length];
        for (var i = 1; i < nodes.Length; i++)
        {
            var parent = nodes[i].Parent;
            var delta = positions[i] - positions[parent];
            var distance = new Vector(delta.X, delta.Y).Length;
            var direction = Direction(parent, i, distance);
            var force = direction * ((distance - MinimumDistance(parent, i) * 1.05) * options.SpringStrength);
            forces[i] -= force;
            forces[parent] += force;
        }
        EdgeAvoidance.Apply(nodes, positions, forces, options, cancellationToken);
        for (var i = 1; i < nodes.Length; i++)
        {
            velocities[i] = (velocities[i] + forces[i]) * options.Damping;
            if (velocities[i].Length > options.MaxSpeed)
                velocities[i] *= options.MaxSpeed / velocities[i].Length;
            positions[i] += velocities[i];
        }

        // Project exclusion violations after integration. Descendants are exempt only
        // from their ancestor's radius, not from their own radius facing that ancestor.
        MaximumOverlap = 0;
        for (var sweep = 0; sweep < 32; sweep++)
        {
            var maximumOverlap = 0d;
            var cells = new Dictionary<(int X, int Y), List<int>>();
            for (var i = 0; i < nodes.Length; i++)
            {
                cancellationToken.ThrowIfCancellationRequested();
                var key = Cell(positions[i]);
                for (var x = key.X - 1; x <= key.X + 1; x++)
                for (var y = key.Y - 1; y <= key.Y + 1; y++)
                {
                    if (!cells.TryGetValue((x, y), out var neighbors))
                        continue;
                    foreach (var j in neighbors)
                    {
                        var delta = positions[i] - positions[j];
                        var distance = new Vector(delta.X, delta.Y).Length;
                        var overlap = MinimumDistance(j, i) - distance;
                        if (overlap <= 0)
                            continue;
                        maximumOverlap = Math.Max(maximumOverlap, overlap);
                        var correction = Direction(j, i, distance) * overlap;
                        positions[i] += correction * (j == 0 ? 1 : 0.5);
                        if (j != 0)
                            positions[j] -= correction * 0.5;
                    }
                }
                var updated = Cell(positions[i]);
                if (!cells.TryGetValue(updated, out var members))
                    cells[updated] = members = [];
                members.Add(i);
            }
            MaximumOverlap = maximumOverlap;
            if (maximumOverlap < options.Tolerance / 10)
                break;
        }
        MaximumOverlap = MeasureOverlap(cancellationToken);
        Iterations++;
        IsStable = false;
        var maximum = 0d;
        for (var i = 1; i < nodes.Length; i++)
        {
            var delta = positions[i] - before[i];
            maximum = Math.Max(maximum, new Vector(delta.X, delta.Y).Length);
        }
        return maximum;
    }

    private double MeasureOverlap(CancellationToken cancellationToken)
    {
        var maximum = 0d;
        var cells = new Dictionary<(int X, int Y), List<int>>();
        for (var i = 0; i < nodes.Length; i++)
        {
            cancellationToken.ThrowIfCancellationRequested();
            var key = Cell(positions[i]);
            for (var x = key.X - 1; x <= key.X + 1; x++)
            for (var y = key.Y - 1; y <= key.Y + 1; y++)
            {
                if (!cells.TryGetValue((x, y), out var neighbors))
                    continue;
                foreach (var j in neighbors)
                {
                    var delta = positions[i] - positions[j];
                    maximum = Math.Max(maximum, MinimumDistance(j, i) - new Vector(delta.X, delta.Y).Length);
                }
            }
            if (!cells.TryGetValue(key, out var members))
                cells[key] = members = [];
            members.Add(i);
        }
        return maximum;
    }

    private (int X, int Y) Cell(Point point) =>
        ((int)Math.Floor(point.X / cellSize), (int)Math.Floor(point.Y / cellSize));

    private Vector Direction(int first, int second, double distance)
    {
        if (distance > 1e-9)
            return new Vector(positions[second].X - positions[first].X, positions[second].Y - positions[first].Y)
                / distance;
        var angle = (first + second + 1) * 2.399963229728653;
        return new Vector(Math.Cos(angle), Math.Sin(angle));
    }
}