Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
07/22/10 00:44:01 (14 years ago)
Author:
swagner
Message:

Sorted usings and removed unused usings in entire solution (#1094)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.Netron/3.0.2672.12446/Netron.Diagramming.Core-3.0.2672.12446/Layout/Force/NBodyForce.cs

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