1 | using System;
|
---|
2 | using System.Collections.Generic;
|
---|
3 | using System.Diagnostics;
|
---|
4 | namespace Netron.Diagramming.Core.Layout.Force {
|
---|
5 | /// <summary>
|
---|
6 | /// <para>
|
---|
7 | /// Force function which computes an n-body force such as gravity,
|
---|
8 | /// anti-gravity, or the results of electric charges. This function implements
|
---|
9 | /// the the Barnes-Hut algorithm for efficient n-body force simulations,
|
---|
10 | /// using a quad-tree with aggregated mass values to compute the n-body
|
---|
11 | /// force in O(N log N) time, where N is the number of ForceItems.</para>
|
---|
12 | ///
|
---|
13 | /// <para>The algorithm used is that of J. Barnes and P. Hut, in their research
|
---|
14 | /// paper <i>A Hierarchical O(n log n) force calculation algorithm</i>, Nature,
|
---|
15 | /// v.324, December 1986. For more details on the algorithm, see one of
|
---|
16 | /// the following links --
|
---|
17 | /// <list type="bullet">
|
---|
18 | ///
|
---|
19 | /// <item><a href="http://www.cs.berkeley.edu/~demmel/cs267/lecture26/lecture26.html">James Demmel's UC Berkeley lecture notes</a></item>
|
---|
20 | /// <item><a href="http://www.physics.gmu.edu/~large/lr_forces/desc/bh/bhdesc.html">Description of the Barnes-Hut algorithm</a></item>
|
---|
21 | /// <item><a href="http://www.ifa.hawaii.edu/~barnes/treecode/treeguide.html">Joshua Barnes' recent implementation</a></item>
|
---|
22 | /// </list></para>
|
---|
23 | /// </summary>
|
---|
24 | public class NBodyForce : AbstractForce {
|
---|
25 |
|
---|
26 | /*
|
---|
27 | The indexing scheme for quadtree child nodes goes row by row.
|
---|
28 | 0 | 1 0 -> top left, 1 -> top right
|
---|
29 | -------
|
---|
30 | 2 | 3 2 -> bottom left, 3 -> bottom right
|
---|
31 | */
|
---|
32 |
|
---|
33 | #region Fields
|
---|
34 | private static String[] pnames = new String[] { "GravitationalConstant",
|
---|
35 | "Distance", "BarnesHutTheta" };
|
---|
36 |
|
---|
37 | public static float DEFAULT_GRAV_CONSTANT = -1.0f;
|
---|
38 | public static float DEFAULT_MIN_GRAV_CONSTANT = -10f;
|
---|
39 | public static float DEFAULT_MAX_GRAV_CONSTANT = 10f;
|
---|
40 |
|
---|
41 | public static float DEFAULT_DISTANCE = -1f;
|
---|
42 | public static float DEFAULT_MIN_DISTANCE = -1f;
|
---|
43 | public static float DEFAULT_MAX_DISTANCE = 500f;
|
---|
44 |
|
---|
45 | public static float DEFAULT_THETA = 0.9f;
|
---|
46 | public static float DEFAULT_MIN_THETA = 0.0f;
|
---|
47 | public static float DEFAULT_MAX_THETA = 1.0f;
|
---|
48 |
|
---|
49 | public static int GRAVITATIONAL_CONST = 0;
|
---|
50 | public static int MIN_DISTANCE = 1;
|
---|
51 | public static int BARNES_HUT_THETA = 2;
|
---|
52 |
|
---|
53 | private float xMin, xMax, yMin, yMax;
|
---|
54 | private QuadTreeNodeFactory factory = new QuadTreeNodeFactory();
|
---|
55 | private QuadTreeNode root;
|
---|
56 |
|
---|
57 | private Random rand = new Random(12345678); // deterministic randomness
|
---|
58 |
|
---|
59 | #endregion
|
---|
60 |
|
---|
61 | #region Properties
|
---|
62 |
|
---|
63 | /**
|
---|
64 | * Gets whether this is a force item.
|
---|
65 | */
|
---|
66 | public override bool IsItemForce {
|
---|
67 | get {
|
---|
68 | return true;
|
---|
69 | }
|
---|
70 | }
|
---|
71 |
|
---|
72 | /// <summary>
|
---|
73 | /// Gets the parameter names.
|
---|
74 | /// </summary>
|
---|
75 | /// <value></value>
|
---|
76 | /// <returns></returns>
|
---|
77 | protected override String[] ParameterNames {
|
---|
78 | get {
|
---|
79 | return pnames;
|
---|
80 | }
|
---|
81 | }
|
---|
82 | #endregion
|
---|
83 |
|
---|
84 | #region Constructor
|
---|
85 |
|
---|
86 |
|
---|
87 |
|
---|
88 | /// <summary>
|
---|
89 | /// Create a new NBodyForce with default parameters.
|
---|
90 | /// </summary>
|
---|
91 | public NBodyForce()
|
---|
92 | : this(DEFAULT_GRAV_CONSTANT, DEFAULT_DISTANCE, DEFAULT_THETA) {
|
---|
93 |
|
---|
94 | }
|
---|
95 | #endregion
|
---|
96 |
|
---|
97 | #region Methods
|
---|
98 | /// <summary>
|
---|
99 | /// Create a new NBodyForce.
|
---|
100 | /// </summary>
|
---|
101 | /// <param name="gravConstant">the gravitational constant to use. Nodes will attract each other if this value is positive, and will repel each other if it is negative.</param>
|
---|
102 | /// <param name="minDistance">the distance within which two particles will interact. If -1, the value is treated as infinite.</param>
|
---|
103 | /// <param name="theta">the Barnes-Hut parameter theta, which controls when an aggregated mass is used rather than drilling down to individual item mass values.</param>
|
---|
104 | public NBodyForce(float gravConstant, float minDistance, float theta) {
|
---|
105 | parms = new float[] { gravConstant, minDistance, theta };
|
---|
106 | minValues = new float[] { DEFAULT_MIN_GRAV_CONSTANT,
|
---|
107 | DEFAULT_MIN_DISTANCE, DEFAULT_MIN_THETA };
|
---|
108 | maxValues = new float[] { DEFAULT_MAX_GRAV_CONSTANT,
|
---|
109 | DEFAULT_MAX_DISTANCE, DEFAULT_MAX_THETA };
|
---|
110 | root = factory.GetQuadTreeNode();
|
---|
111 | }
|
---|
112 |
|
---|
113 |
|
---|
114 |
|
---|
115 |
|
---|
116 | /// <summary>
|
---|
117 | /// Set the bounds of the region for which to compute the n-body simulation
|
---|
118 | /// </summary>
|
---|
119 | /// <param name="xMin">the minimum x-coordinate</param>
|
---|
120 | /// <param name="yMin">the minimum y-coordinate/param>
|
---|
121 | /// <param name="xMax"> the maximum x-coordinate</param>
|
---|
122 | /// <param name="yMax">the maximum y-coordinate</param>
|
---|
123 | private void setBounds(float xMin, float yMin, float xMax, float yMax) {
|
---|
124 | this.xMin = xMin;
|
---|
125 | this.yMin = yMin;
|
---|
126 | this.xMax = xMax;
|
---|
127 | this.yMax = yMax;
|
---|
128 | }
|
---|
129 |
|
---|
130 |
|
---|
131 | /// <summary>
|
---|
132 | ///Clears the quadtree of all entries.
|
---|
133 | /// </summary>
|
---|
134 | public void clear() {
|
---|
135 | ClearHelper(root);
|
---|
136 | root = factory.GetQuadTreeNode();
|
---|
137 | }
|
---|
138 |
|
---|
139 | /// <summary>
|
---|
140 | /// Clearing aid
|
---|
141 | /// </summary>
|
---|
142 | /// <param name="n">The n.</param>
|
---|
143 | private void ClearHelper(QuadTreeNode n) {
|
---|
144 | for (int i = 0; i < n.children.Length; i++) {
|
---|
145 | if (n.children[i] != null)
|
---|
146 | ClearHelper(n.children[i]);
|
---|
147 | }
|
---|
148 | factory.Reclaim(n);
|
---|
149 | }
|
---|
150 |
|
---|
151 | /// <summary>
|
---|
152 | /// Initialize the simulation with the provided enclosing simulation. After
|
---|
153 | /// this call has been made, the simulation can be queried for the
|
---|
154 | /// n-body force acting on a given item.
|
---|
155 | /// </summary>
|
---|
156 | /// <param name="fsim">the encompassing ForceSimulator</param>
|
---|
157 | public override void Init(ForceSimulator fsim) {
|
---|
158 | clear(); // clear internal state
|
---|
159 |
|
---|
160 | // compute and squarify bounds of quadtree
|
---|
161 | float x1 = float.MaxValue, y1 = float.MaxValue;
|
---|
162 | float x2 = float.MinValue, y2 = float.MinValue;
|
---|
163 | foreach (ForceItem item in fsim.Items) {
|
---|
164 | float x = item.Location[0];
|
---|
165 | float y = item.Location[1];
|
---|
166 | if (x < x1) x1 = x;
|
---|
167 | if (y < y1) y1 = y;
|
---|
168 | if (x > x2) x2 = x;
|
---|
169 | if (y > y2) y2 = y;
|
---|
170 | }
|
---|
171 | float dx = x2 - x1, dy = y2 - y1;
|
---|
172 | if (dx > dy) { y2 = y1 + dx; } else { x2 = x1 + dy; }
|
---|
173 | setBounds(x1, y1, x2, y2);
|
---|
174 |
|
---|
175 | // Insert items into quadtree
|
---|
176 | foreach (ForceItem item in fsim.Items) {
|
---|
177 | Insert(item);
|
---|
178 | }
|
---|
179 |
|
---|
180 | // calculate magnitudes and centers of mass
|
---|
181 | CalculateMass(root);
|
---|
182 | }
|
---|
183 |
|
---|
184 | /// <summary>
|
---|
185 | /// Inserts an item into the quadtree.
|
---|
186 | /// </summary>
|
---|
187 | /// <param name="item"> the ForceItem to add.</param>
|
---|
188 | public void Insert(ForceItem item) {
|
---|
189 | // Insert item into the quadtrees
|
---|
190 | try {
|
---|
191 | Insert(item, root, xMin, yMin, xMax, yMax);
|
---|
192 | }
|
---|
193 | catch (StackOverflowException e) {
|
---|
194 | // TODO: safe to remove?
|
---|
195 | Trace.WriteLine(e.Message);
|
---|
196 | }
|
---|
197 | }
|
---|
198 |
|
---|
199 | /// <summary>
|
---|
200 | /// Inserts the specified force.
|
---|
201 | /// </summary>
|
---|
202 | /// <param name="p">The p.</param>
|
---|
203 | /// <param name="n">The n.</param>
|
---|
204 | /// <param name="x1">The x1.</param>
|
---|
205 | /// <param name="y1">The y1.</param>
|
---|
206 | /// <param name="x2">The x2.</param>
|
---|
207 | /// <param name="y2">The y2.</param>
|
---|
208 | private void Insert(ForceItem p, QuadTreeNode n, float x1, float y1, float x2, float y2) {
|
---|
209 | // try to Insert particle p at node n in the quadtree
|
---|
210 | // by construction, each leaf will contain either 1 or 0 particles
|
---|
211 | if (n.hasChildren) {
|
---|
212 | // n contains more than 1 particle
|
---|
213 | InsertHelper(p, n, x1, y1, x2, y2);
|
---|
214 | } else if (n.value != null) {
|
---|
215 | // n contains 1 particle
|
---|
216 | if (IsSameLocation(n.value, p)) {
|
---|
217 | InsertHelper(p, n, x1, y1, x2, y2);
|
---|
218 | } else {
|
---|
219 | ForceItem v = n.value; n.value = null;
|
---|
220 | InsertHelper(v, n, x1, y1, x2, y2);
|
---|
221 | InsertHelper(p, n, x1, y1, x2, y2);
|
---|
222 | }
|
---|
223 | } else {
|
---|
224 | // n is empty, so is a leaf
|
---|
225 | n.value = p;
|
---|
226 | }
|
---|
227 | }
|
---|
228 |
|
---|
229 | /// <summary>
|
---|
230 | /// Determines whether the two force are at the same location.
|
---|
231 | /// </summary>
|
---|
232 | /// <param name="f1">The f1.</param>
|
---|
233 | /// <param name="f2">The f2.</param>
|
---|
234 | /// <returns>
|
---|
235 | /// <c>true</c> if [is same location] [the specified f1]; otherwise, <c>false</c>.
|
---|
236 | /// </returns>
|
---|
237 | private static bool IsSameLocation(ForceItem f1, ForceItem f2) {
|
---|
238 | float dx = Math.Abs(f1.Location[0] - f2.Location[0]);
|
---|
239 | float dy = Math.Abs(f1.Location[1] - f2.Location[1]);
|
---|
240 | return (dx < 0.01 && dy < 0.01);
|
---|
241 | }
|
---|
242 |
|
---|
243 | /// <summary>
|
---|
244 | /// Inserts helper method.
|
---|
245 | /// </summary>
|
---|
246 | /// <param name="p">The p.</param>
|
---|
247 | /// <param name="n">The n.</param>
|
---|
248 | /// <param name="x1">The x1.</param>
|
---|
249 | /// <param name="y1">The y1.</param>
|
---|
250 | /// <param name="x2">The x2.</param>
|
---|
251 | /// <param name="y2">The y2.</param>
|
---|
252 | private void InsertHelper(ForceItem p, QuadTreeNode n, float x1, float y1, float x2, float y2) {
|
---|
253 | float x = p.Location[0], y = p.Location[1];
|
---|
254 | float splitx = (x1 + x2) / 2;
|
---|
255 | float splity = (y1 + y2) / 2;
|
---|
256 | int i = (x >= splitx ? 1 : 0) + (y >= splity ? 2 : 0);
|
---|
257 | // create new child node, if necessary
|
---|
258 | if (n.children[i] == null) {
|
---|
259 | n.children[i] = factory.GetQuadTreeNode();
|
---|
260 | n.hasChildren = true;
|
---|
261 | }
|
---|
262 | // update bounds
|
---|
263 | if (i == 1 || i == 3) x1 = splitx; else x2 = splitx;
|
---|
264 | if (i > 1) y1 = splity; else y2 = splity;
|
---|
265 | // recurse
|
---|
266 | Insert(p, n.children[i], x1, y1, x2, y2);
|
---|
267 | }
|
---|
268 |
|
---|
269 | /// <summary>
|
---|
270 | /// Calculates the mass.
|
---|
271 | /// </summary>
|
---|
272 | /// <param name="n">The n.</param>
|
---|
273 | private void CalculateMass(QuadTreeNode n) {
|
---|
274 | float xcom = 0, ycom = 0;
|
---|
275 | n.mass = 0;
|
---|
276 | if (n.hasChildren) {
|
---|
277 | for (int i = 0; i < n.children.Length; i++) {
|
---|
278 | if (n.children[i] != null) {
|
---|
279 | CalculateMass(n.children[i]);
|
---|
280 | n.mass += n.children[i].mass;
|
---|
281 | xcom += n.children[i].mass * n.children[i].com[0];
|
---|
282 | ycom += n.children[i].mass * n.children[i].com[1];
|
---|
283 | }
|
---|
284 | }
|
---|
285 | }
|
---|
286 | if (n.value != null) {
|
---|
287 | n.mass += n.value.Mass;
|
---|
288 | xcom += n.value.Mass * n.value.Location[0];
|
---|
289 | ycom += n.value.Mass * n.value.Location[1];
|
---|
290 | }
|
---|
291 | n.com[0] = xcom / n.mass;
|
---|
292 | n.com[1] = ycom / n.mass;
|
---|
293 | }
|
---|
294 |
|
---|
295 | /**
|
---|
296 | * Calculates the force vector acting on the given item.
|
---|
297 | * @param item the ForceItem for which to compute the force
|
---|
298 | */
|
---|
299 | public override void GetForce(ForceItem item) {
|
---|
300 | try {
|
---|
301 | ForceHelper(item, root, xMin, yMin, xMax, yMax);
|
---|
302 | }
|
---|
303 | catch (StackOverflowException e) {
|
---|
304 | // TODO: safe to remove?
|
---|
305 | Trace.WriteLine(e.Message);
|
---|
306 | }
|
---|
307 | }
|
---|
308 |
|
---|
309 | /// <summary>
|
---|
310 | /// Utility method.
|
---|
311 | /// </summary>
|
---|
312 | /// <param name="item">The item.</param>
|
---|
313 | /// <param name="n">The n.</param>
|
---|
314 | /// <param name="x1">The x1.</param>
|
---|
315 | /// <param name="y1">The y1.</param>
|
---|
316 | /// <param name="x2">The x2.</param>
|
---|
317 | /// <param name="y2">The y2.</param>
|
---|
318 | private void ForceHelper(ForceItem item, QuadTreeNode n, float x1, float y1, float x2, float y2) {
|
---|
319 | float dx = n.com[0] - item.Location[0];
|
---|
320 | float dy = n.com[1] - item.Location[1];
|
---|
321 | float r = (float)Math.Sqrt(dx * dx + dy * dy);
|
---|
322 | bool same = false;
|
---|
323 | if (r == 0.0f) {
|
---|
324 | // if items are in the exact same place, add some noise
|
---|
325 | dx = Convert.ToSingle((rand.NextDouble() - 0.5D) / 50.0D);
|
---|
326 | dy = Convert.ToSingle((rand.NextDouble() - 0.5D) / 50.0D);
|
---|
327 | r = (float)Math.Sqrt(dx * dx + dy * dy);
|
---|
328 | same = true;
|
---|
329 | }
|
---|
330 | bool minDist = parms[MIN_DISTANCE] > 0f && r > parms[MIN_DISTANCE];
|
---|
331 |
|
---|
332 | // the Barnes-Hut approximation criteria is if the ratio of the
|
---|
333 | // size of the quadtree box to the distance between the point and
|
---|
334 | // the box's center of mass is beneath some threshold theta.
|
---|
335 | if ((!n.hasChildren && n.value != item) ||
|
---|
336 | (!same && (x2 - x1) / r < parms[BARNES_HUT_THETA])) {
|
---|
337 | if (minDist) return;
|
---|
338 | // either only 1 particle or we meet criteria
|
---|
339 | // for Barnes-Hut approximation, so calc force
|
---|
340 | float v = parms[GRAVITATIONAL_CONST] * item.Mass * n.mass
|
---|
341 | / (r * r * r);
|
---|
342 | item.Force[0] += v * dx;
|
---|
343 | item.Force[1] += v * dy;
|
---|
344 | } else if (n.hasChildren) {
|
---|
345 | // recurse for more accurate calculation
|
---|
346 | float splitx = (x1 + x2) / 2;
|
---|
347 | float splity = (y1 + y2) / 2;
|
---|
348 | for (int i = 0; i < n.children.Length; i++) {
|
---|
349 | if (n.children != null && n.children[i] != null) {
|
---|
350 | ForceHelper(item, n.children[i],
|
---|
351 | (i == 1 || i == 3 ? splitx : x1), (i > 1 ? splity : y1),
|
---|
352 | (i == 1 || i == 3 ? x2 : splitx), (i > 1 ? y2 : splity));
|
---|
353 | }
|
---|
354 | }
|
---|
355 | if (minDist) return;
|
---|
356 | if (n.value != null && n.value != item) {
|
---|
357 | float v = parms[GRAVITATIONAL_CONST] * item.Mass * n.value.Mass
|
---|
358 | / (r * r * r);
|
---|
359 | item.Force[0] += v * dx;
|
---|
360 | item.Force[1] += v * dy;
|
---|
361 | }
|
---|
362 | }
|
---|
363 | }
|
---|
364 | #endregion
|
---|
365 |
|
---|
366 | #region Classes
|
---|
367 | /// <summary>
|
---|
368 | /// Represents a node in the quadtree.
|
---|
369 | /// </summary>
|
---|
370 | public sealed class QuadTreeNode {
|
---|
371 | public QuadTreeNode() {
|
---|
372 | com = new float[] { 0.0f, 0.0f };
|
---|
373 | children = new QuadTreeNode[4];
|
---|
374 | } //
|
---|
375 | public bool hasChildren = false;
|
---|
376 | public float mass; // total mass held by this node
|
---|
377 | public float[] com; // center of mass of this node
|
---|
378 | public ForceItem value; // ForceItem in this node, null if node has children
|
---|
379 | public QuadTreeNode[] children; // children nodes
|
---|
380 | }
|
---|
381 |
|
---|
382 |
|
---|
383 |
|
---|
384 |
|
---|
385 | ///<summary>
|
---|
386 | ///Helper class to minimize number of object creations across multiple uses of the quadtree.
|
---|
387 | ///</summary>
|
---|
388 | public sealed class QuadTreeNodeFactory {
|
---|
389 | #region Fields
|
---|
390 | private int maxNodes = 50000;
|
---|
391 | private List<QuadTreeNode> nodes = new List<QuadTreeNode>();
|
---|
392 | #endregion
|
---|
393 |
|
---|
394 | /// <summary>
|
---|
395 | /// Gets the quadtree node.
|
---|
396 | /// </summary>
|
---|
397 | /// <returns></returns>
|
---|
398 | public QuadTreeNode GetQuadTreeNode() {
|
---|
399 | if (nodes.Count > 0) {
|
---|
400 | QuadTreeNode node = this.nodes[nodes.Count - 1];
|
---|
401 | nodes.Remove(node);
|
---|
402 | return node;
|
---|
403 | } else {
|
---|
404 | return new QuadTreeNode();
|
---|
405 | }
|
---|
406 | }
|
---|
407 | /// <summary>
|
---|
408 | /// Reclaims the specified node.
|
---|
409 | /// </summary>
|
---|
410 | /// <param name="n">The n.</param>
|
---|
411 | public void Reclaim(QuadTreeNode n) {
|
---|
412 | n.mass = 0;
|
---|
413 | n.com[0] = 0.0f; n.com[1] = 0.0f;
|
---|
414 | n.value = null;
|
---|
415 | n.hasChildren = false;
|
---|
416 | for (int i = 0; i < n.children.Length; i++) {
|
---|
417 | n.children[i] = null;
|
---|
418 | }
|
---|
419 | //n.children = null;
|
---|
420 | if (nodes.Count < maxNodes)
|
---|
421 | nodes.Add(n);
|
---|
422 | }
|
---|
423 | }
|
---|
424 | #endregion
|
---|
425 |
|
---|
426 | }
|
---|
427 |
|
---|
428 | }
|
---|