Nbody

Python Numpy

from __future__ import print_function
"""
NBody in N^2 complexity

Note that we are using only Newtonian forces and do not consider relativity
Neither do we consider collisions between stars
Thus some of our stars will accelerate to speeds beyond c
This is done to keep the simulation simple enough for teaching purposes

All the work is done in the calc_force, move and random_galaxy functions.
To vectorize the code these are the functions to transform.
"""
from benchpress import util
if util.Benchmark().bohrium:
    import bohrium as np
else:
    import numpy as np

G     = 6.67384e-11     # m/(kg*s^2)
dt    = 60*60*24*365.25 # Years in seconds
r_ly  = 9.4607e15       # Lightyear in m
m_sol = 1.9891e30       # Solar mass in kg

def diagonal(ary, offset=0):
    """
    Return specified diagonals.

    If `a` is 2-D, returns the diagonal of `a` with the given offset,
    i.e., the collection of elements of the form ``a[i, i+offset]``.

    Parameters
    ----------
    ary : array_like
        Array from which the diagonals are taken.
    offset : int, optional
        Offset of the diagonal from the main diagonal.  Can be positive or
        negative.  Defaults to main diagonal (0).

    Returns
    -------
    array_of_diagonals : ndarray
        If `a` is 2-D, a 1-D array containing the diagonal is returned.
        If the dimension of `a` is larger, then an array of diagonals is
        returned, "packed" from left-most dimension to right-most (e.g.,
        if `a` is 3-D, then the diagonals are "packed" along rows).

    Raises
    ------
    ValueError
        If the dimension of `a` is less than 2.

    See Also
    --------
    diag : MATLAB work-a-like for 1-D and 2-D arrays.
    diagflat : Create diagonal arrays.
    trace : Sum along diagonals.

    Examples
    --------
    >>> a = np.arange(4).reshape(2,2)
    >>> a
    array([[0, 1],
           [2, 3]])
    >>> a.diagonal()
    array([0, 3])
    >>> a.diagonal(1)
    array([1])

    A 3-D example:

    >>> a = np.arange(8).reshape(2,2,2); a
    array([[[0, 1],
            [2, 3]],
           [[4, 5],
            [6, 7]]])
    """
    if ary.ndim != 2:
        raise Exception("diagonal only supports 2 dimensions\n")
    if offset < 0:
        offset = -offset
        if (ary.shape[0]-offset) > ary.shape[1]:
            ary_diag = ary[offset, :]
        else:
            ary_diag = ary[offset:, 0]
    else:
        if ary.shape[1]-offset > ary.shape[0]:
            ary_diag = ary[:, offset]
        else:
            ary_diag = ary[0, offset:]
    ary_diag.strides = (ary.strides[0]+ary.strides[1],)
    return ary_diag


def random_galaxy(N, B, dtype=np.float64):
    """Generate a galaxy of random bodies"""

    galaxy = {            # We let all bodies stand still initially
        'm':    (B.random_array((N,), dtype=dtype) + dtype(10)) * dtype(m_sol/10),
        'x':    (B.random_array((N,), dtype=dtype) - dtype(0.5)) * dtype(r_ly/100),
        'y':    (B.random_array((N,), dtype=dtype) - dtype(0.5)) * dtype(r_ly/100),
        'z':    (B.random_array((N,), dtype=dtype) - dtype(0.5)) * dtype(r_ly/100),
        'vx':   np.zeros(N, dtype=dtype),
        'vy':   np.zeros(N, dtype=dtype),
        'vz':   np.zeros(N, dtype=dtype)
    }
    if dtype == np.float32:
        galaxy['m'] /= 1e10
        galaxy['x'] /= 1e5
        galaxy['y'] /= 1e5
        galaxy['z'] /= 1e5
    return galaxy

def move(galaxy, dt):
    """Move the bodies
    first find forces and change velocity and then move positions
    """
    n  = len(galaxy['x'])
    # Calculate all dictances component wise (with sign)
    dx = galaxy['x'][np.newaxis,:].T - galaxy['x']
    dy = galaxy['y'][np.newaxis,:].T - galaxy['y']
    dz = galaxy['z'][np.newaxis,:].T - galaxy['z']

    # Euclidian distances (all bodys)
    r = np.sqrt(dx**2 + dy**2 + dz**2)
    diagonal(r)[:] = 1.0

    # prevent collition
    mask = r < 1.0
    r = r * ~mask + 1.0 * mask

    m = galaxy['m'][np.newaxis,:].T

    # Calculate the acceleration component wise
    Fx = G*m*dx/r**3
    Fy = G*m*dy/r**3
    Fz = G*m*dz/r**3
    # Set the force (acceleration) a body exerts on it self to zero
    diagonal(Fx)[:] = 0.0
    diagonal(Fy)[:] = 0.0
    diagonal(Fz)[:] = 0.0

    galaxy['vx'] += dt*np.sum(Fx, axis=0)
    galaxy['vy'] += dt*np.sum(Fy, axis=0)
    galaxy['vz'] += dt*np.sum(Fz, axis=0)

    galaxy['x'] += dt*galaxy['vx']
    galaxy['y'] += dt*galaxy['vy']
    galaxy['z'] += dt*galaxy['vz']

def simulate(galaxy, timesteps, visualize=False):
    for i in range(timesteps):
        move(galaxy,dt)
        util.Benchmark().flush()
        if visualize:#NB: this is only for experiments
            T = np.zeros((3, len(galaxy['x'])), dtype=np.float32)
            T[0,:] = galaxy['x']
            T[1,:] = galaxy['y']
            T[2,:] = galaxy['z']
            np.visualize(T, "3d", 0, 0.0, 10)

def main():
    B = util.Benchmark()
    N = B.size[0]
    I = B.size[1]

    if B.inputfn:
        galaxy = B.load_arrays(B.inputfn)
    else:
        galaxy = random_galaxy(N, B, B.dtype)
        util.Benchmark().flush()

    if B.dumpinput:
        B.dump_arrays("nbody", galaxy)

    B.start()
    simulate(galaxy, I, visualize=B.visualize)
    R = galaxy['x'] + galaxy['y'] + galaxy['z']
    B.stop()

    B.pprint()
    if B.outputfn:
        B.tofile(B.outputfn, {'res':R})

if __name__ == "__main__":
    main()

Csharp Numcil

#region Copyright
/*
This file is part of Bohrium and copyright (c) 2012 the Bohrium:
team <http://www.bh107.org>.

Bohrium is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation, either version 3
of the License, or (at your option) any later version.

Bohrium is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the
GNU Lesser General Public License along with Bohrium.

If not, see <http://www.gnu.org/licenses/>.
*/
#endregion

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

using R = NumCIL.Range;

namespace nbody
{
	using NumCIL.Double;
	using DATA = System.Double;

	public static class nBodySolverDouble
	{
		//Gravity
		private static DATA G = (DATA)1.0;

		private static void FillDiagonal(NdArray a, DATA val)
		{
			long d  = a.Shape.Dimensions[0].Length;
			a.Reshape(new NumCIL.Shape(new long[] { d }, 0, new long[] { d+1 })).Set(val);
		}

		private static void CalcForce(Galaxy b)
		{
			var dx = b.x - b.x[R.NewAxis, R.All].Transposed;
			FillDiagonal(dx, 1);
			var dy = b.y - b.y[R.NewAxis, R.All].Transposed;
			FillDiagonal(dy, 1);
			var dz = b.z - b.z[R.NewAxis, R.All].Transposed;
			FillDiagonal(dz, 1);
			var pm = b.mass - b.mass[R.NewAxis, R.All].Transposed;
			FillDiagonal(pm, 0);

			var r = (dx.Pow(2) + dy.Pow(2) + dz.Pow(2)).Pow((DATA)0.5);

			//In the below calc of the the forces the force of a body upon itself
			//becomes nan and thus destroys the data
			var Fx = G * pm / r.Pow(2) * (dx / r);
			var Fy = G * pm / r.Pow(2) * (dy / r);
			var Fz = G * pm / r.Pow(2) * (dz / r);

			//The diagonal nan numbers must be removed so that the force from a body
			//upon itself is zero
			FillDiagonal(Fx,0);
			FillDiagonal(Fy,0);
			FillDiagonal(Fz,0);

			b.vx += Add.Reduce(Fx, 1) / b.mass;
			b.vy += Add.Reduce(Fy, 1) / b.mass;
			b.vz += Add.Reduce(Fz, 1) / b.mass;
		}

		private static void Move(Galaxy galaxy)
		{
			CalcForce(galaxy);

			galaxy.x += galaxy.vx;
			galaxy.y += galaxy.vy;
			galaxy.z += galaxy.vz;
		}

		public class Galaxy
		{
			private const DATA XMAX = 500;
			private const DATA YMAX = 500;
			private const DATA ZMAX = 500;

			public NdArray mass;
			public NdArray x;
			public NdArray y;
			public NdArray z;
			public NdArray vx;
			public NdArray vy;
			public NdArray vz;

			public Galaxy(long size)
			{
				this.mass = Generate.Random(size) * (DATA)Math.Pow(10, 6) / (DATA)(4 * Math.PI * Math.PI);
				this.x = Generate.Random(size) * 2 * XMAX - XMAX;
				this.y = Generate.Random(size) * 2 * YMAX - YMAX;
				this.z = Generate.Random(size) * 2 * ZMAX - ZMAX;
				this.vx = Generate.Zeroes(size);
				this.vy = Generate.Zeroes(size);
				this.vz = Generate.Zeroes(size);
			}

			public DATA Sync()
			{
				return this.z.Value[0] + this.y.Value[0] + this.x.Value[0];
			}
		}

		public static void Solve(Galaxy galaxy, long steps)
		{
			for(long i = 0; i < steps; i++)
				Move(galaxy);
		}
	}
}

namespace nbody
{
	using NumCIL.Double;
	using DATA = System.Double;

	public static class nBodySolverDoubleNoTempArrays
	{
		//Gravity
		private const DATA G = (DATA)1.0;
		private const DATA MIN_DIST = (DATA)1.0;

		private static void FillDiagonal(NdArray a, DATA val)
		{
			long d  = a.Shape.Dimensions[0].Length;
			a.Reshape(new NumCIL.Shape(new long[] { d }, 0, new long[] { d+1 })).Set(val);
		}

		private static void CalcForce(Galaxy b)
		{
			// Calculate all dictances component wise (with sign)
			b.x.Sub(b.x[R.NewAxis, R.All].Transposed, b.dx);
			b.y.Sub(b.y[R.NewAxis, R.All].Transposed, b.dy);
			b.z.Sub(b.z[R.NewAxis, R.All].Transposed, b.dz);

			// Euclidian distances (all bodys)
			b.dx.Pow(2, b.r0);
			b.dy.Pow(2, b.r1);
			b.r0.Add(b.r1, b.r0);
			b.dz.Pow(2, b.r1);
			b.r0.Add(b.r1, b.r0);
			b.r0.Pow((DATA)0.5, b.r0);
			FillDiagonal(b.r0, (Double)1.0);

			// prevent collition
			NumCIL.Double.GreaterThanOrEqual.Apply(b.r0, MIN_DIST, b.mask);
			b.mask.ToDouble (b.r1);
			b.r1.Mul (b.r0, b.r0);

			NumCIL.Boolean.Not.Apply (b.mask, b.mask);
			b.mask.ToDouble (b.r1);
			b.r1.Mul (MIN_DIST, b.r1);
			b.r1.Add (b.r0, b.r0);

			//Calculate the acceleration component wise
			var m = b.mass[R.NewAxis, R.All].Transposed;
			b.dx.Mul(G,b.dx);
			b.dy.Mul(G,b.dy);
			b.dz.Mul(G,b.dz);
			b.dx.Mul(m,b.dx);
			b.dy.Mul(m,b.dy);
			b.dz.Mul(m,b.dz);
			b.r0.Pow(3, b.r0);
			b.dx.Div(b.r0, b.dx);
			b.dy.Div(b.r0, b.dy);
			b.dz.Div(b.r0, b.dz);

			// Set the force (acceleration) a body exerts on it self to zero
			FillDiagonal(b.dx, 0);
			FillDiagonal(b.dy, 0);
			FillDiagonal(b.dz, 0);

			Double dt = (Double) (60*60*24*365.25); // Years in seconds
			b.dx.Reduce<Add>(0, b.t0);
			b.t0.Mul(dt, b.t0);
			b.vx.Add(b.t0, b.vx);
			b.dy.Reduce<Add>(0, b.t0);
			b.t0.Mul(dt, b.t0);
			b.vy.Add(b.t0, b.vy);
			b.dz.Reduce<Add>(0, b.t0);
			b.t0.Mul(dt, b.t0);
			b.vz.Add(b.t0, b.vz);
		}

		private static void Move(Galaxy galaxy)
		{
			CalcForce(galaxy);

			galaxy.x.Add(galaxy.vx, galaxy.x);
			galaxy.y.Add(galaxy.vy, galaxy.y);
			galaxy.z.Add(galaxy.vz, galaxy.z);
		}

		public class Galaxy
		{
			private const DATA XMAX = 500;
			private const DATA YMAX = 500;
			private const DATA ZMAX = 500;

			public NdArray mass;
			public NdArray x;
			public NdArray y;
			public NdArray z;
			public NdArray vx;
			public NdArray vy;
			public NdArray vz;
			public NdArray dx;
			public NdArray dy;
			public NdArray dz;
			public NdArray pmass;
			public NdArray r0;
			public NdArray r1;
			public NdArray t0;
			public NumCIL.Boolean.NdArray mask;

			public Galaxy(long size)
			{
				this.mass = Generate.Random(size) * (DATA)Math.Pow(10, 6) / (DATA)(4 * Math.PI * Math.PI);
				this.x = Generate.Random(size) * 2 * XMAX - XMAX;
				this.y = Generate.Random(size) * 2 * YMAX - YMAX;
				this.z = Generate.Random(size) * 2 * ZMAX - ZMAX;
				this.vx = Generate.Zeroes(size);
				this.vy = Generate.Zeroes(size);
				this.vz = Generate.Zeroes(size);
				this.dx = new NdArray(new [] {size, size});
				this.dy = new NdArray(new [] {size, size});
				this.dz = new NdArray(new [] {size, size});
				this.pmass = new NdArray(new [] {size, size});
				this.r0 = new NdArray(this.dx.Shape);
				this.r1 = new NdArray(this.dx.Shape);
				this.t0 = new NdArray(this.x.Shape);
				this.mask = new NumCIL.Boolean.NdArray(this.dx.Shape);
			}

			public DATA Sync()
			{
				return this.z.Value[0] + this.y.Value[0] + this.x.Value[0];
			}
		}

		public static void Solve(Galaxy galaxy, long steps)
		{
			for(long i = 0; i < steps; i++)
				Move(galaxy);
		}
	}
}

namespace nbody
{
	using NumCIL.Float;
	using DATA = System.Single;

	public static class nBodySolverFloat
	{
		//Gravity
		private static DATA G = (DATA)1.0;

		private static void FillDiagonal(NdArray a, DATA val)
		{
			long d  = a.Shape.Dimensions[0].Length;
			a.Reshape(new NumCIL.Shape(new long[] { d }, 0, new long[] { d+1 })).Set(val);
		}

		private static void CalcForce(Galaxy b)
		{
			var dx = b.x - b.x[R.NewAxis, R.All].Transposed;
			FillDiagonal(dx, 1);
			var dy = b.y - b.y[R.NewAxis, R.All].Transposed;
			FillDiagonal(dy, 1);
			var dz = b.z - b.z[R.NewAxis, R.All].Transposed;
			FillDiagonal(dz, 1);
			var pm = b.mass - b.mass[R.NewAxis, R.All].Transposed;
			FillDiagonal(pm, 0);

			var r = (dx.Pow(2) + dy.Pow(2) + dz.Pow(2)).Pow((DATA)0.5);

			//In the below calc of the the forces the force of a body upon itself
			//becomes nan and thus destroys the data
			var Fx = G * pm / r.Pow(2) * (dx / r);
			var Fy = G * pm / r.Pow(2) * (dy / r);
			var Fz = G * pm / r.Pow(2) * (dz / r);

			//The diagonal nan numbers must be removed so that the force from a body
			//upon itself is zero
			FillDiagonal(Fx,0);
			FillDiagonal(Fy,0);
			FillDiagonal(Fz,0);

			b.vx += Add.Reduce(Fx, 1) / b.mass;
			b.vy += Add.Reduce(Fy, 1) / b.mass;
			b.vz += Add.Reduce(Fz, 1) / b.mass;
		}

		private static void Move(Galaxy galaxy)
		{
			CalcForce(galaxy);

			galaxy.x += galaxy.vx;
			galaxy.y += galaxy.vy;
			galaxy.z += galaxy.vz;
		}

		public class Galaxy
		{
			private const DATA XMAX = 500;
			private const DATA YMAX = 500;
			private const DATA ZMAX = 500;

			public NdArray mass;
			public NdArray x;
			public NdArray y;
			public NdArray z;
			public NdArray vx;
			public NdArray vy;
			public NdArray vz;

			public Galaxy(long size)
			{
				this.mass = Generate.Random(size) * (DATA)Math.Pow(10, 6) / (DATA)(4 * Math.PI * Math.PI);
				this.x = Generate.Random(size) * 2 * XMAX - XMAX;
				this.y = Generate.Random(size) * 2 * YMAX - YMAX;
				this.z = Generate.Random(size) * 2 * ZMAX - ZMAX;
				this.vx = Generate.Zeroes(size);
				this.vy = Generate.Zeroes(size);
				this.vz = Generate.Zeroes(size);
			}

			public DATA Sync()
			{
				return this.z.Value[0] + this.y.Value[0] + this.x.Value[0];
			}
		}

		public static void Solve(Galaxy galaxy, long steps)
		{
			for(long i = 0; i < steps; i++)
				Move(galaxy);
		}
	}
}

namespace nbody
{
	using NumCIL.Float;
	using DATA = System.Single;

	public static class nBodySolverFloatNoTempArrays
	{
		//Gravity
		private const DATA G = (DATA)1.0;
		private const DATA MIN_DIST = (DATA)1.0;

		private static void FillDiagonal(NdArray a, DATA val)
		{
			long d  = a.Shape.Dimensions[0].Length;
			a.Reshape(new NumCIL.Shape(new long[] { d }, 0, new long[] { d+1 })).Set(val);
		}

		private static void CalcForce(Galaxy b)
		{
			// Calculate all dictances component wise (with sign)
			b.x.Sub(b.x[R.NewAxis, R.All].Transposed, b.dx);
			b.y.Sub(b.y[R.NewAxis, R.All].Transposed, b.dy);
			b.z.Sub(b.z[R.NewAxis, R.All].Transposed, b.dz);

			// Euclidian distances (all bodys)
			b.dx.Pow(2, b.r0);
			b.dy.Pow(2, b.r1);
			b.r0.Add(b.r1, b.r0);
			b.dz.Pow(2, b.r1);
			b.r0.Add(b.r1, b.r0);
			b.r0.Pow((DATA)0.5, b.r0);
			FillDiagonal(b.r0, (Single)1.0);

			// prevent collition
			NumCIL.Float.GreaterThanOrEqual.Apply(b.r0, MIN_DIST, b.mask);
			b.mask.ToFloat (b.r1);
			b.r1.Mul (b.r0, b.r0);

			NumCIL.Boolean.Not.Apply (b.mask, b.mask);
			b.mask.ToFloat (b.r1);
			b.r1.Mul (MIN_DIST, b.r1);
			b.r1.Add (b.r0, b.r0);

			//Calculate the acceleration component wise
			var m = b.mass[R.NewAxis, R.All].Transposed;
			b.dx.Mul(G,b.dx);
			b.dy.Mul(G,b.dy);
			b.dz.Mul(G,b.dz);
			b.dx.Mul(m,b.dx);
			b.dy.Mul(m,b.dy);
			b.dz.Mul(m,b.dz);
			b.r0.Pow(3, b.r0);
			b.dx.Div(b.r0, b.dx);
			b.dy.Div(b.r0, b.dy);
			b.dz.Div(b.r0, b.dz);

			// Set the force (acceleration) a body exerts on it self to zero
			FillDiagonal(b.dx, 0);
			FillDiagonal(b.dy, 0);
			FillDiagonal(b.dz, 0);

			Single dt = (Single) (60*60*24*365.25); // Years in seconds
			b.dx.Reduce<Add>(0, b.t0);
			b.t0.Mul(dt, b.t0);
			b.vx.Add(b.t0, b.vx);
			b.dy.Reduce<Add>(0, b.t0);
			b.t0.Mul(dt, b.t0);
			b.vy.Add(b.t0, b.vy);
			b.dz.Reduce<Add>(0, b.t0);
			b.t0.Mul(dt, b.t0);
			b.vz.Add(b.t0, b.vz);
		}

		private static void Move(Galaxy galaxy)
		{
			CalcForce(galaxy);

			galaxy.x.Add(galaxy.vx, galaxy.x);
			galaxy.y.Add(galaxy.vy, galaxy.y);
			galaxy.z.Add(galaxy.vz, galaxy.z);
		}

		public class Galaxy
		{
			private const DATA XMAX = 500;
			private const DATA YMAX = 500;
			private const DATA ZMAX = 500;

			public NdArray mass;
			public NdArray x;
			public NdArray y;
			public NdArray z;
			public NdArray vx;
			public NdArray vy;
			public NdArray vz;
			public NdArray dx;
			public NdArray dy;
			public NdArray dz;
			public NdArray pmass;
			public NdArray r0;
			public NdArray r1;
			public NdArray t0;
			public NumCIL.Boolean.NdArray mask;

			public Galaxy(long size)
			{
				this.mass = Generate.Random(size) * (DATA)Math.Pow(10, 6) / (DATA)(4 * Math.PI * Math.PI);
				this.x = Generate.Random(size) * 2 * XMAX - XMAX;
				this.y = Generate.Random(size) * 2 * YMAX - YMAX;
				this.z = Generate.Random(size) * 2 * ZMAX - ZMAX;
				this.vx = Generate.Zeroes(size);
				this.vy = Generate.Zeroes(size);
				this.vz = Generate.Zeroes(size);
				this.dx = new NdArray(new [] {size, size});
				this.dy = new NdArray(new [] {size, size});
				this.dz = new NdArray(new [] {size, size});
				this.pmass = new NdArray(new [] {size, size});
				this.r0 = new NdArray(this.dx.Shape);
				this.r1 = new NdArray(this.dx.Shape);
				this.t0 = new NdArray(this.x.Shape);
				this.mask = new NumCIL.Boolean.NdArray(this.dx.Shape);
			}

			public DATA Sync()
			{
				return this.z.Value[0] + this.y.Value[0] + this.x.Value[0];
			}
		}

		public static void Solve(Galaxy galaxy, long steps)
		{
			for(long i = 0; i < steps; i++)
				Move(galaxy);
		}
	}
}