# Barnes-Hut implementation of the N-Body problem translated from F# to C#

I am after a peer review of a C# implementation of the Barnes-Hut algorithm which I have translated from F#. The F# version is the base for comparison, therefore the C# version is suppose to reflect the F# in structure and behaviour. I have tried to translate the F# implementation as accurately as possible, however i would like a second opinion.

I would be grateful for any advice e.g. on any problems, or on how to make the algorithm further reflect the F# version. Also, any suggestions on how improve parallel performance, if any.

F# version main algorthim http://pastebin.com/AMb9AnNd

C# version:

/*
*
*
* usage:
* args[0] = algorithm version
* args[2] = body count
* args[3] = core count
* args[4] = test discription
* args[5] = eps value
*/

using System;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;

namespace BHTuned
{
#region BHMain: attributes
/// <summary>
///
/// </summary>
public class Body
{
public double x { get; set; }
public double y { get; set; }
public double z { get; set; }
public double vx { get; set; }
public double vy { get; set; }
public double vz { get; set; }
public double m { get; set; }
public Body(double x, double y, double z, double vx, double vy, double vz, double mass)
{
this.x = x;
this.y = y;
this.z = z;
this.vx = vx;
this.vy = vy;
this.vz = vz;
this.m = mass;
}
}

/// <summary>
///
/// </summary>
public class Accel
{
public double ax { get; set; }
public double ay { get; set; }
public double az { get; set; }
public Accel(double ax, double ay, double az)
{
this.ax = ax;
this.ay = ay;
this.az = az;
}
}
/// <summary>
///
/// </summary>
public class Centroid
{
public double cx { get; set; }
public double cy { get; set; }
public double cz { get; set; }
public double cm { get; set; }
public Centroid(double cx, double cy, double cz, double cm)
{
this.cx = cx;
this.cy = cy;
this.cz = cz;
this.cm = cm;
}
}
/// <summary>
///
/// </summary>
public class Bbox
{
public List<Body> bs;
public double minx { get; set; }
public double miny { get; set; }
public double minz { get; set; }
public double maxx { get; set; }
public double maxy { get; set; }
public double maxz { get; set; }
public Bbox(double minx, double miny, double minz, double maxx, double maxy, double maxz)
{
bs = new List<Body>();
this.minx = minx;
this.miny = miny;
this.minz = minz;
this.maxx = maxx;
this.maxy = maxy;
this.maxz = maxz;
}
}
/// <summary>
///
/// </summary>
public class BHTree
{
public List<BHTree> subTrees = new List<BHTree>();
public double size { get; set; }
public double centerx { get; set; }
public double centery { get; set; }
public double centerz { get; set; }
public double mass { get; set; }

public BHTree(double size, double centerx, double centery, double centerz, double mass, List<BHTree> subTrees)
{
this.size = size;
this.centerx = centerx;
this.centery = centery;
this.centerz = centerz;
this.mass = mass;
this.subTrees = subTrees;
}
}
#endregion
/// <summary>
///
/// </summary>
/// <param name="ba"></param>
/// <returns></returns>
public abstract class BHMain
{
// If the distance between the points is smaller than this
public double eps; //= 0.01F;    // then ignore the forces between them.
protected static float timeStep = 0.001F;
public static float threshold = 0.25F;
public int N;

public BHMain(int N, double eps) { this.N = N; this.eps = eps; }
/// <summary>
/// Determine 3D boundaries
/// </summary>
/// <param name="ba"> All Bodys </param>
/// <returns> Bounding box </returns>
public static Bbox findBounds(List<Body> ba)
{
double minx = 0.0, miny = 0.0, minz = 0.0, maxx = 0.0, maxy = 0.0, maxz = 0.0;
foreach (Body b in ba)
{
minx = Math.Min(minx, b.x);
miny = Math.Min(minx, b.y);
minz = Math.Min(minx, b.z);
maxx = Math.Max(maxx, b.x);
maxy = Math.Max(maxx, b.y);
maxz = Math.Max(maxx, b.z);
}
return new Bbox(minx, miny, minz, maxx, maxy, maxz);
}
/// <summary>
///
/// </summary>
/// <param name="bl"></param>
/// <returns></returns>
public static Centroid calcCentroid(List<Body> bl)
{
double x = 0.0, y = 0.0, z = 0.0, m = 0.0;
foreach (Body b in bl)
{
x = b.m * b.x + x;
y = b.m * b.y + y;
z = b.m * b.z + z;
m = m + b.m;
}
return new Centroid(x / m, y / m, z / m, m);
}
/// <summary>
/// Generate Bodys
/// </summary>
/// <param name="n"></param>
/// <returns>Bodys List </returns>
public List<Body> getBodies(int seed)
{
List<Body> d = new List<Body>();
Random rnd = new Random(seed);
for (int i = 0; i < seed; i++)
{
rnd.NextDouble(),
rnd.NextDouble(),
rnd.NextDouble(),
rnd.NextDouble(),
rnd.NextDouble(),
rnd.NextDouble()));
}

/*
* test data
d.Add(new Body(0.72624327, 0.8173253596, 0.7680226894, 0.5581611914, 0.206033154, 0.5588847946, 0.906027066));
d.Add(new Body(0.2486685842, 0.1107439772, 0.4670106799, 0.771604122, 0.6575188938, 0.4327826013, 0.3540837636));
d.Add(new Body(0.7710938983, 0.4041625948, 0.1659986703, 0.9850470526, 0.1090046336, 0.306680408, 0.8021404612));
d.Add(new Body(0.2935192125, 0.6975812124, 0.8649866608, 0.1984899832, 0.5604903733, 0.1805782147, 0.2501971588));
d.Add(new Body(0.8159445267, 0.99099983, 0.5639746513, 0.4119329138, 0.01197611308, 0.05447602135, 0.6982538564));
d.Add(new Body(0.3383698409, 0.2844184475, 0.2629626418, 0.6253758444, 0.4634618528, 0.928373828, 0.146310554));
d.Add(new Body(0.8607951551, 0.5778370651, 0.9619506323, 0.8388187749, 0.9149475926, 0.8022716347, 0.5943672515));
d.Add(new Body(0.3832204693, 0.8712556827, 0.6609386227, 0.05226170553, 0.3664333324, 0.6761694414, 0.04242394913));
d.Add(new Body(0.9056457835, 0.1646743003, 0.3599266132, 0.2657046361, 0.8179190721, 0.5500672481, 0.4904806467));
d.Add(new Body(0.4280710977, 0.4580929179, 0.05891460369, 0.4791475667, 0.2694048119, 0.4239650548, 0.9385373443));
*/
return d;
}
/// <summary>
/// Determine a bodies position in relation to bounding box
/// </summary>
/// <param name="bb"></param>
/// <param name="b"></param>
/// <returns></returns>
public static bool inbox(Bbox bb, Body b)
{
if ((b.x > bb.minx) && (b.x <= bb.maxx) && (b.y > bb.miny) && (b.y <= bb.maxy) && (b.z > bb.minz) && (b.z <= bb.maxz))
return true;
else
return false;
}
/// <summary>
///
/// </summary>
/// <param name="bb"></param>
/// <param name="bl"></param>
/// <returns></returns>
public static Bbox[] splitPoints(Bbox bb, List<Body> bl)
{
Bbox[] box = new Bbox[8];
double midx = (bb.maxx + bb.minx) / 2;
double midy = (bb.maxy + bb.miny) / 2;
double midz = (bb.maxz + bb.minz) / 2;

for (int i = 0; i < 8; i++)
box[i] = new Bbox(bb.minx, bb.miny, bb.minz, midx, midy, midz);

foreach (Body b in bl)
for (int i = 0; i < 8; i++)

Bbox[] bbf = new Bbox[8];
for (int i = 0; i < 8; i++)
if (box[i].bs.Count != 0) bbf[i] = (box[i]);
return bbf;
}
/// <summary>
/// Populate the oct-tree
/// </summary>
/// <param name="bb"></param>
/// <param name="bl"></param>
/// <returns></returns>
public static int count = 0;
public static BHTree buildTree(Bbox bb, List<Body> bl, int N)
{
List<BHTree> subTrees = new List<BHTree>();
double[] s = { Math.Abs(bb.maxx - bb.minx), Math.Abs(bb.maxy - bb.miny), Math.Abs(bb.maxz - bb.minz) };
double size = s.Min();
Centroid cent = calcCentroid(bl);

count++;
Bbox[] boxesAndPts = splitPoints(bb, bl);
if (boxesAndPts[0] == null)
return new BHTree(size, cent.cx, cent.cy, cent.cz, cent.cm, subTrees);
else
{
foreach (Bbox b in boxesAndPts)
return new BHTree(size, cent.cx, cent.cy, cent.cz, cent.cm, subTrees);
}
}
/// <summary>
/// calculate the acceleration on a point due to some other point
/// </summary>
/// <returns></returns>
public static Accel accel(BHTree BHT, Body b, double eps)
{
double dx = b.x - BHT.centerx, dy = b.y - BHT.centery, dz = b.z - BHT.centerz;
double dsqr = (dx * dx) + (dy * dy) + (dz * dz) + eps;
double d = Math.Sqrt(dsqr);
double mag = timeStep / (dsqr * d);
double mmag = BHT.mass * mag;
return new Accel((dx * mmag), (dy * mmag), (dz * mmag));
}
/// <summary>
///
/// </summary>
/// <param name="BHT"></param>
/// <param name="b"></param>
/// <returns></returns>
public static bool isFar(BHTree BHT, Body b)
{
double dx = b.x - BHT.centerx, dy = b.y - BHT.centery, dz = b.z - BHT.centerz;
double dist = Math.Sqrt(dx * dx + dy * dy + dz * dz);
if ((BHT.size / dist) < threshold) return true; else return false;
}
/// <summary>
///
/// </summary>
/// <param name="a1"></param>
/// <param name="a2"></param>
/// <returns></returns>
public static Accel addAccel(Accel a1, Accel a2)
{
return new Accel((a1.ax + a2.ax), (a1.ay + a2.ay), (a1.az + a2.az));
}
/// <summary>
/// Recursively calculate relative body acceleration within the tree
/// </summary>
/// <param name="BHT">Populated oct-tree</param>
/// <param name="b">Current body</param>
/// <returns></returns>
public static Accel calcAccel(BHTree BHT, Body b, double eps)
{
if (BHT.subTrees.Count == 0 || isFar(BHT, b))
return accel(BHT, b, eps);
else
{
Accel acc = new Accel(0, 0, 0);
foreach (BHTree t in BHT.subTrees)
{
Accel a = (calcAccel(t, b, eps));
}
return acc;
}
}
/// <summary>
/// Checking function
/// </summary>
/// <param name="acc"></param>
/// <param name="b"></param>
/// <returns></returns>
public double f(Double acc, Body b) { return acc + b.x + b.y + b.z + b.vx + b.vy + b.vz; }

/// <summary>
///
/// </summary>
/// <param name="args"></param>
public virtual void doSteps(string args, string testName, int coreCount) { }
}
/// <summary>
/// Sequential Abtraction
/// </summary>
public class Seq : BHMain
{
public Seq(int N, double eps)
: base(N, eps)
{
this.N = N; this.eps = eps;
}
public override void doSteps(string args, string testDescription, int coreCount)
{
Console.WriteLine("Sequential" + " Cores: " + coreCount);
List<Body> ba = getBodies(N);

List<Body> bl = ba;
SetFiles sf = new SetFiles(testDescription, testDescription + "/seq/", "seq", coreCount);
Stopwatch stopwatch = new Stopwatch();
stopwatch.Start();
for (int i = 0; i < 16; i++)
{
Bbox bb = findBounds(ba);
BHTree BHT = buildTree(bb, ba, N);
for (int j = 0; j < N; j++)
{
Accel acc = calcAccel(BHT, ba[j], eps);
Body nb = new Body(ba[j].x, ba[j].y, ba[j].z, (ba[j].vx - acc.ax), (ba[j].vy - acc.ay), (ba[j].vz - acc.az), ba[j].m);
bl[j] = (new Body((nb.x + timeStep * nb.vx), (nb.y + timeStep * nb.vy), (nb.z + timeStep * nb.vz), nb.vx, nb.vy, nb.vz, nb.m));
}
ba = bl;
}
stopwatch.Stop();
Console.WriteLine(stopwatch.Elapsed.TotalSeconds);
//   sf.buildDataFile();
//   sf.buildRefFile();
}
}
/// <summary>
/// Parallel.For abstraction
/// </summary>
public class PFor : BHMain
{
public PFor(int N, double eps)
: base(N, eps)
{
this.N = N; this.eps = eps;
}
public override void doSteps(string args, string testName, int coreCount)
{
List<Body> ba = getBodies(N);
List<Body> bl = ba; //populate list for iteration

/* Naive parallel ditribution.  N Parallel loop iterations */
if (args == "bf") // brute-force
{
Console.WriteLine("Brute-force PFor" + " Cores: " + coreCount);
SetFiles sf = new SetFiles(testName, testName + "/pf/bf/", "pfbf", coreCount);
Stopwatch stopwatch = new Stopwatch();
stopwatch.Start();
for (int i = 0; i < 16; i++)
{
Bbox bb = findBounds(ba);
BHTree BHT = buildTree(bb, ba, N);
try
{
Parallel.For(0, (N - 1), j =>
{
Accel acc = calcAccel(BHT, ba[j], eps);
Body nb = new Body(ba[j].x, ba[j].y, ba[j].z, (ba[j].vx - acc.ax), (ba[j].vy - acc.ay), (ba[j].vz - acc.az), ba[j].m);
bl[j] = (new Body((nb.x + timeStep * nb.vx), (nb.y + timeStep * nb.vy), (nb.z + timeStep * nb.vz), nb.vx, nb.vy, nb.vz, nb.m));
});
}
catch (AggregateException e) { Console.Write(e); }
}
stopwatch.Stop();
sf.buildDataFile();
sf.buildRefFile();
}
/* Tuned parallel ditribution.  Course grainuality, load is evenly balanced */
else if (args == "lb") // load balanced
{
Console.WriteLine("load balanced PFor" + " Cores: " + coreCount);
SetFiles sf = new SetFiles(testName, testName + "/pf/lb/", "pflb", coreCount);
Stopwatch stopwatch = new Stopwatch();
int partition = N / 4;
BHTree BHT;
Bbox bb;
stopwatch.Start();

for (int k = 0; k < 16; k++)
{
bb = findBounds(ba);
BHT = buildTree(bb, ba, N);
Parallel.For(0, 4, new ParallelOptions { MaxDegreeOfParallelism = 4 }, j =>
{
int c = 0, d = partition;
if (j == 0) { c = 0; d = partition; }
else if (j == 1) { c = partition; d = partition * 2; }
else if (j == 2) { c = partition * 2; d = partition * 3; }
else if (j == 3) { c = partition * 3; d = partition * 4; }
for (; c < d; c++)
{
Accel acc = calcAccel(BHT, ba[c], eps);
Body nb = new Body(ba[c].x, ba[c].y, ba[c].z, (ba[c].vx - acc.ax), (ba[c].vy - acc.ay), (ba[c].vz - acc.az), ba[c].m);
bl[c] = (new Body((nb.x + timeStep * nb.vx), (nb.y + timeStep * nb.vy), (nb.z + timeStep * nb.vz), nb.vx, nb.vy, nb.vz, nb.m));
}
ba = bl;
});
}

stopwatch.Stop();
Console.WriteLine(stopwatch.Elapsed.TotalSeconds);
// sf.buildDataFile();
// sf.buildRefFile();
}
}
}
/// <summary>
/// Parallel.Invoke abstraction
/// </summary>
public class PInvoke : BHMain
{
public PInvoke(int N, double eps)
: base(N, eps)
{
this.N = N; this.eps = eps;
}
public override void doSteps(string args, string testName, int coreCount)
{
List<Body> ba = getBodies(N);
Bbox bb = findBounds(ba);
BHTree BHT = buildTree(bb, ba, N);
List<Body> bl = ba;
/* Naive parallel ditribution.  N Parallel.Invoke's are created */
if (args == "bf")
{
Console.WriteLine("Brute-force PInvoke" + " Cores: " + coreCount);
SetFiles sf = new SetFiles(testName, testName + "/pi/bf/", "pibf", coreCount);
Stopwatch stopwatch = new Stopwatch();
stopwatch.Start();
for (int i = 0; i < 16; i++)
for (int j = 0; j < N; j++)
{
Parallel.Invoke(() =>
{
Accel acc = calcAccel(BHT, ba[j], eps);
Body nb = new Body(ba[j].x, ba[j].y, ba[j].z, (ba[j].vx - acc.ax), (ba[j].vy - acc.ay), (ba[j].vz - acc.az), ba[j].m);
bl[j] = (new Body((nb.x + timeStep * nb.vx), (nb.y + timeStep * nb.vy), (nb.z + timeStep * nb.vz), nb.vx, nb.vy, nb.vz, nb.m));
});
ba = bl;
}

stopwatch.Stop();
Console.WriteLine(stopwatch.Elapsed.TotalSeconds);
// sf.buildDataFile();
// sf.buildRefFile();
}

/* Tuned parallel ditribution.  Course grainuality, load is evenly balanced */
else if (args == "lb")
{
Console.Write("load balanced PInvoke" + " Cores: " + coreCount+"\n");
int partition = N / 4;
SetFiles sf = new SetFiles(testName, testName + "/pi/lb/", "pilb", coreCount);
Stopwatch stopwatch = new Stopwatch();

stopwatch.Start();
for (int k = 0; k < 4; k++) // iteration loop
for (int j = 0; j < 4; j++) // distribution loop
{
Parallel.Invoke(new ParallelOptions { MaxDegreeOfParallelism = 4 },() =>
{
int c = 0, d = partition;
if (j == 0) { c = 0; d = partition; }
else if (j == 1) { c = partition; d = partition * 2; }
else if (j == 2) { c = partition * 2; d = partition * 3; }
else if (j == 3) { c = partition * 3; d = partition * 4; }
for (; c < d; c++)
{
Accel acc = calcAccel(BHT, ba[c], eps);
Body nb = new Body(ba[c].x, ba[c].y, ba[c].z, (ba[c].vx - acc.ax), (ba[c].vy - acc.ay), (ba[c].vz - acc.az), ba[c].m);
bl[c] = (new Body((nb.x + timeStep * nb.vx), (nb.y + timeStep * nb.vy), (nb.z + timeStep * nb.vz), nb.vx, nb.vy, nb.vz, nb.m));
}
});
ba = bl;
}

stopwatch.Stop();
Console.WriteLine(stopwatch.Elapsed.TotalSeconds);
//sf.buildDataFile();
//sf.buildRefFile();
}
}
}
}

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;

namespace BHTuned
{
class EntryPoint
{
/// <summary>
///
/// </summary>
/// <param name="args"> </param>
static void Main(string[] args)
{

Process Proc = Process.GetCurrentProcess();
if (Convert.ToInt16(args[3] )== 1) AffinityMask &= 0x0001; // use only any of the first 4 available processors;
else if (Convert.ToInt16(args[3]) == 2) AffinityMask &= 0x0003;
else if (Convert.ToInt16(args[3]) == 3) AffinityMask &= 0x0007;
else if (Convert.ToInt16(args[3]) == 4) AffinityMask &= 0x000F;

// Render r = new Render("esp");
switch (args[0])
{
case "seq": // Sequential
BHMain bf = new Seq(Convert.ToInt32(args[2]), Convert.ToDouble(args[5]));
bf.doSteps(args[1], args[4], Convert.ToInt32(args[3]));
break;
case "pf": // parallel.For
BHMain pf = new PFor(Convert.ToInt32(args[2]), Convert.ToDouble(args[5]));
pf.doSteps(args[1], args[4], Convert.ToInt32(args[3]));
break;
case "pi": // parallel.Invoke
BHMain pi = new PInvoke(Convert.ToInt32(args[2]), Convert.ToDouble(args[5]));
pi.doSteps(args[1], args[4], Convert.ToInt32(args[3]));
break;
}
}
}
}

• Please include the code you want to get reviewed in your question (you can keep external links). This prevents your question and the answers to become useless when the link goes away. – ChrisWue Nov 10 '13 at 17:54

I think there is a lot in the code which can be improved especially from design point of view. A few points which caught my eye follow but it can probably improved further. I would start with what I've written below and then possibly go through another round of refactoring.

Naming conventions.

1. C# standard naming convention for methods and public properties is CamelCase.
2. Please don't try and save characters on class names. Ideally reading the class name should give you a good indication of what it's purpose. It's not easy to find a good yet concise name but it pays of to spent 5min thinking about it. Some suggestions:
• BBox -> BoundingBox
• Accel -> Acceleration
• BHMain -> BarnesHutBase
• Seq -> BarnesHutSequential
• PFor -> BarnesHutParallel
3. Please give your local variables some more meaningful names. Reading all this bb, ba, b, bl, etc. just hurts the eye. Intellisense and auto completions have been invented a long time ago, no need to try and save characters for sake of typing. Using names like boundingBox, bodies, etc. would make the code much easier to read.

Design

1. Consider making your plain data classes like Accel, Centroid etc. structs and immutable. Their properties don't change after they have been constructed if I read your code correctly so you might want to reflect this in your design. Immutability makes a lot of things easier. Eric Lippert has written some very interesting blog articles about immutability. At least consider making the properties which don't change either readonly fields or give them a private set.
2. The findBounds method in BHMain should probably live as a factory method in BBox instead.
3. The calcCentroid method in BHMain should probably live as a factory method in Centroid instead.
4. splitPoints should live in BBox (factory method)
5. buildTree should live in BHTree (factory method)
6. accel should live in Accel (factory method)
7. inbox should live as IsInBox on BBox (non-static)
8. There seems to be a magic for i: 1 -> 16 loop in your doSteps implementations which seems odd.
9. Once you actually remove all those static helper classes from BHMain then there is not much left in BHMain so ask yourself: What should be the purpose of this class?
10. You repeat a lot of code in Seq, PFor and PInvoke - especially the last two share a lot. The central point of the calculations seems to be this block:

Accel acc = calcAccel(BHT, ba[j], eps);
Body nb = new Body(ba[j].x, ba[j].y, ba[j].z, (ba[j].vx - acc.ax), (ba[j].vy - acc.ay), (ba[j].vz - acc.az), ba[j].m);
bl[j] = (new Body((nb.x + timeStep * nb.vx), (nb.y + timeStep * nb.vy), (nb.z + timeStep * nb.vz), nb.vx, nb.vy, nb.vz, nb.m));


You should make this a SingleStep method on BHMain which you can call when needed (which gives some purpose back to BHMain)

I would consider a design along these lines:

class BarnesHutBase
{
protected void SingleStep(....)
{
...
}
}

interface IStepable
{
void PerformSteps(int numberOfSteps);
}

class BarnesHutSequential : BarnesHutBase, IStepable
{
public void PerformSteps(int numberOfSteps)
{
...
}
}

11. The timing and recording of things should be encapsulated in SetFiles as a RecordTime method. Something along these lines:

class RecordingSetFiles : IDisposable
{
private Stopwatch _Stopwatch ;
private SetFiles _SetFiles;
public RecordingSetFiles(SetFiles files)
{
_SetFiles = files;
_Stopwatch = new Stopwatch();
_Stopwatch.Start();
}

public void Dispose()
{
_Stopwatch.Stop();
_SetFiles.buildDataFile();
_SetFiles.buildRefFile();
}
}

class SetFiles
{
public RecordingSetFiles RecordTime()
{
return new RecordingSetFiles(this);
}
}


Then you can use it like this:

using (sf.RecordTime())
{
....
}


Possible Bugs

1. In PFor and PInvoke you do this:

List<Body> bl = ba; //populate list for iteration


This will not create a new list bl but just assign the reference to bl, so bl and ba now point to the same list. Did you actually meant to do bl = new List<Body>(ba)?

• Many thanks for the advice. I will defo re-factor further. Something else... when testing, my results differ slightly from the F# version. Over one iteration each body value mantissa is accurate to four decimal places. From the code samples, it is possible to tell why? however, again, cheers for your input. – Stuart Gordon Nov 10 '13 at 23:10