using System;
using System.Collections.Generic;
using System.Linq;namespace GraphdotNet.ForceAtlas2Algo
{internal class ForceAtlas2{#region Constructors and Destructorspublic ForceAtlas2(int nodeNum){nodes = new Node[nodeNum];for (var i = 0; i < nodes.Length; i++){nodes[i] = new Node();}Degree = new int[nodeNum];}#endregion#region Methodsprivate double GetEdgeWeight(int edgeIdx){return Adjacency[3*edgeIdx + 2];}#endregion#region Constants and Fieldsinternal double Speed;private readonly Node[] nodes;private double outboundAttCompensation = 1;private Region rootRegion;#endregion#region Propertiesinternal bool AdjustSizes { get; set; }internal bool BarnesHutOptimize { get; set; }internal double BarnesHutTheta { get; set; }internal double EdgeWeightInfluence { get; set; }internal double Gravity { get; set; }internal double JitterTolerance { get; set; }internal bool LinLogMode { get; set; }internal bool OutboundAttractionDistribution { get; set; }internal double ScalingRatio { get; set; }internal bool StrongGravityMode { get; set; }private double[] Adjacency { get; set; }private int[] Degree { get; }private int[] EdgeAdjacency { get; set; }#endregion#region Public Methodspublic bool CanAlgo(){return (nodes.Length != 0);}public void EndAlgo(){}public IEnumerable<double> GetX(){return nodes.Select(node => node.X);}public IEnumerable<double> GetY(){return nodes.Select(node => node.Y);}public void GoAlgo(){if (nodes.Length == 0){return;}initialize node datavar nodeNum = nodes.Length;foreach (var n in nodes){n.Olddx = n.Dx;n.Olddy = n.Dy;n.Dx = 0;n.Dy = 0;}// If Barnes Hut active, initialize root RegionImplif (BarnesHutOptimize){rootRegion = new Region(nodes);rootRegion.BuildSubRegions();}// If outboundAttractionDistribution active, compensate.if (OutboundAttractionDistribution){outboundAttCompensation = 0;foreach (var n in nodes){outboundAttCompensation += n.Mass;}outboundAttCompensation /= nodeNum;}// Repulsion (and gravity)var repulsion = ForceFactory.Builder.BuildRepulsion(AdjustSizes, ScalingRatio);const int @from = 0;var to = nodeNum;var rgCalculator = new RepulsionGravityCalculate(nodes,from,to,BarnesHutOptimize,BarnesHutTheta,Gravity,StrongGravityMode ? (ForceFactory.Builder.GetStrongGravity(ScalingRatio)) : (repulsion),ScalingRatio,rootRegion,repulsion);rgCalculator.Run();// Attractionvar attraction = ForceFactory.Builder.BuildAttraction(LinLogMode,OutboundAttractionDistribution,AdjustSizes,1*(OutboundAttractionDistribution ? (outboundAttCompensation) : (1)));var edgeNum = EdgeAdjacency.Length/2;const double epsilon = 1e-6;if (Math.Abs(EdgeWeightInfluence - 0) < epsilon){for (var i = 0; i < edgeNum; i++){attraction.Apply(nodes[EdgeAdjacency[2*i]], nodes[EdgeAdjacency[2*i + 1]], 1);}}else if (Math.Abs(EdgeWeightInfluence - 1) < epsilon){for (var i = 0; i < edgeNum; i++){attraction.Apply(nodes[EdgeAdjacency[2*i]],nodes[EdgeAdjacency[2*i + 1]],GetEdgeWeight(i));}}else{for (var i = 0; i < edgeNum; i++){attraction.Apply(nodes[EdgeAdjacency[2*i]],nodes[EdgeAdjacency[2*i + 1]],Math.Pow(GetEdgeWeight(i), EdgeWeightInfluence));}}// Auto adjust speedvar totalSwinging = 0d; // How much irregular movementvar totalEffectiveTraction = 0d; // Hom much useful movementforeach (var n in nodes){var swinging = Math.Sqrt((n.Olddx - n.Dx)*(n.Olddx - n.Dx) + (n.Olddy - n.Dy)*(n.Olddy - n.Dy));totalSwinging += n.Mass*swinging;// If the node has a burst change of direction, then it's not converging.totalEffectiveTraction += n.Mass*0.5*Math.Sqrt((n.Olddx + n.Dx)*(n.Olddx + n.Dx) + (n.Olddy + n.Dy)*(n.Olddy + n.Dy));}// We want that swingingMovement < tolerance * convergenceMovementvar targetSpeed = JitterTolerance*JitterTolerance*totalEffectiveTraction/totalSwinging;// But the speed shoudn't rise too much too quickly, since it would make the convergence drop dramatically.const double maxRise = 0.5; // Start rise: 50%Speed = Speed + Math.Min(targetSpeed - Speed, maxRise*Speed);// Apply forcesforeach (var n in nodes){// Adaptive auto-speed: the speed of each node is lowered// when the node swings.var swinging = Math.Sqrt((n.Olddx - n.Dx)*(n.Olddx - n.Dx) + (n.Olddy - n.Dy)*(n.Olddy - n.Dy));//double factor = speed / (1f + Math.sqrt(speed * swinging));var factor = Speed/(1 + Speed*Math.Sqrt(swinging));n.X += n.Dx*factor;n.Y += n.Dy*factor;}}public void InitAlgo(double[] adjacencyInput){foreach (var n in nodes){n.ResetNode();}var ran = new Random();foreach (var n in nodes){n.X = (0.01 + ran.NextDouble())*1000 - 500;n.Y = (0.01 + ran.NextDouble())*1000 - 500;}Speed = 1.0;Adjacency = adjacencyInput;var len = Adjacency.Length;var edgeNum = len/3;EdgeAdjacency = new int[edgeNum*2];for (var i = 0; i < edgeNum; i++){EdgeAdjacency[2*i] = (int) Adjacency[3*i];EdgeAdjacency[2*i + 1] = (int) Adjacency[3*i + 1];}for (var i = 0; i < 2*edgeNum; i++){Degree[EdgeAdjacency[i]]++;}var nodeNum = nodes.Length;for (var i = 0; i < nodeNum; i++){nodes[i].Mass = 1 + Degree[i];}}public void OutlierProcess(){var maxRadius = double.NegativeInfinity;foreach (var n in nodes){if (n.Mass >= 2){maxRadius = Math.Max(maxRadius, Math.Sqrt(n.X*n.X + n.Y*n.Y));}}maxRadius += 10;var maxRadiusCeil = (int) Math.Ceiling(maxRadius);var ran = new Random();foreach (var n in nodes){if (n.Mass <= 1){n.X = ran.Next(-maxRadiusCeil*20, maxRadiusCeil*20)/20.0;n.Y = Math.Pow(-1.0, Math.Floor(n.X))*Math.Sqrt(maxRadius*maxRadius - n.X*n.X);}}}public void ResetPropertiesValues(){var nodeNum = nodes.Length;// TuningScalingRatio = nodeNum >= 100 ? 2.0 : 10.0;StrongGravityMode = false;Gravity = 1.0;// BehaviorOutboundAttractionDistribution = false;LinLogMode = false;EdgeWeightInfluence = 1.0;// Performanceif (nodeNum >= 50000){JitterTolerance = 10.0;}else if (nodeNum >= 5000){JitterTolerance = 1.0;}else{JitterTolerance = 0.1;}BarnesHutOptimize = nodeNum >= 1000;BarnesHutTheta = 1.2;}#endregion}
}