diff --git a/nbody_problem/cpp/main.cpp b/nbody_problem/cpp/main.cpp new file mode 100644 index 0000000..ed51c16 --- /dev/null +++ b/nbody_problem/cpp/main.cpp @@ -0,0 +1,170 @@ +#include +#include +#include +#include + +typedef double real; +using namespace std; + +class Star { +public: + real m; + vector r; + vector v; + vector a, a0; + + Star() { + r.assign(3, 0); + v.assign(3, 0); + a.assign(3, 0); + a0.assign(3, 0); + } + + Star(real mass, vector pos, vector vel) { + m = mass; + r = pos; + v = vel; + } + + friend ostream &operator<<(ostream &so, const Star &si) { + so << si.m << " " << si.r[0] << " " << si.r[1] << " " << si.r[2] << " " + << si.v[0] << " " << si.v[1] << " " << si.v[2] << endl; + return so; + } +}; + +class Cluster : public Star { +protected: +public: + vector s; + Cluster() : Star() {} + + void acceleration() { + for (vector::iterator si = s.begin(); si != s.end(); ++si) + si->a.assign(3, 0); + + for (vector::iterator si = s.begin(); si != s.end(); ++si) { + vector rij(3); + real init = 0.0; + + for (vector::iterator sj = s.begin(); sj != s.end(); ++sj) { + if (si != sj) { + + for (int i = 0; i != 3; ++i) + rij[i] = si->r[i] - sj->r[i]; + + real RdotR = inner_product(rij.begin(), rij.end(), rij.begin(), init); + real apre = 1. / sqrt(RdotR * RdotR * RdotR); + + for (int i = 0; i != 3; ++i) { + si->a[i] = sj->m * -1 * apre * rij[i]; + } + } + } + } + } + + void updatePositions(real dt) { + for (vector::iterator si = s.begin(); si != s.end(); ++si) { + si->a0 = si->a; + for (int i = 0; i != 3; ++i) + si->r[i] += dt * si->v[i] + 0.5 * dt * dt * si->a0[i]; + } + } + + void updateVelocities(real dt) { + + for (vector::iterator si = s.begin(); si != s.end(); ++si) { + for (int i = 0; i != 3; ++i) + si->v[i] += 0.5 * dt * (si->a0[i] + si->a[i]); + si->a0 = si->a; + } + } + + vector energies() { + real init = 0; + vector E(3), rij(3); + E.assign(3, 0); + + for (vector::iterator si = s.begin(); si != s.end(); ++si) + E[1] += 0.5 * si->m * + inner_product(si->v.begin(), si->v.end(), si->v.begin(), init); + + // Potential energy + for (vector::iterator si = s.begin(); si != s.end(); ++si) { + for (vector::iterator sj = si + 1; sj != s.end(); ++sj) { + for (int i = 0; i != 3; ++i) + rij[i] = si->r[i] - sj->r[i]; + E[2] -= si->m * sj->m / + sqrt(inner_product(rij.begin(), rij.end(), rij.begin(), init)); + } + } + E[0] = E[1] + E[2]; + return E; + } + + // Print function + friend ostream &operator<<(ostream &so, Cluster &cl) { + for (vector::iterator si = cl.s.begin(); si != cl.s.end(); ++si) + so << *si; + return so; + } +}; + +int main(int argc, char* argv[]) { + + Cluster cl; + real m; + int dummy; + vector r(3), v(3); + + // Read input data from the command line (makeplummer | dumbp) + do { + cin >> dummy; + cin >> m; + for (int i = 0; i != 3; ++i) + cin >> r[i]; + for (int i = 0; i != 3; ++i) + cin >> v[i]; + cl.s.push_back(Star(m, r, v)); + } while (!cin.eof()); + + cl.s.pop_back(); + + // Compute initial energu of the system + vector E(3), E0(3); + E0 = cl.energies(); + + // Start time, end time and simulation step + real t = 0.0; + real tend; + + if (argc > 1) + tend = strtod(argv[1], NULL); + else + tend = 10.0; + + real dt = 1e-3; + int k = 0; + + for (int i = 0; i < 50; i++) { + cl.acceleration(); + + while (t < tend) { + cl.updatePositions(dt); + + cl.acceleration(); + + cl.updateVelocities(dt); + + t += dt; + k += 1; + if (k % 100 == 0) { + E = cl.energies(); + } + } + t = 0; + } + + return 0; +} diff --git a/nbody_problem/python/compiled_methods.py b/nbody_problem/python/compiled_methods.py new file mode 100644 index 0000000..7cad8ab --- /dev/null +++ b/nbody_problem/python/compiled_methods.py @@ -0,0 +1,108 @@ +import sys +import math +import timeit + +import numpy as np +import pandas as pd + +from transonic import jit + +def load_input_data(path): + + df = pd.read_csv( + path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+" + ) + + masses = df["mass"].values.copy() + positions = df.loc[:, ["x", "y", "z"]].values.copy() + velocities = df.loc[:, ["vx", "vy", "vz"]].values.copy() + + return masses, positions, velocities + +def compute_accelerations(accelerations, masses, positions): + nb_particles = masses.size + + for index_p0 in range(nb_particles - 1): + position0 = positions[index_p0] + mass0 = masses[index_p0] + + # TODO: Use compiled methods like Numba & Pythran in vectorized approach. + # Issue: https://github.com/khushi-411/numpy-benchmarks/issues/4 + for index_p1 in range(index_p0 + 1, nb_particles): + mass1 = masses[index_p1] + + vectors = position0 - positions[index_p1] + + distances = (vectors**2).sum() + coefs = 1./distances**1.5 + + accelerations[index_p0] += mass1 * -1 * vectors * coefs + accelerations[index_p1] += mass0 * vectors * coefs + + return accelerations + +@jit +def loop( + time_step, nb_steps, masses, positions, + velocities): + + accelerations = np.zeros_like(positions) + accelerations1 = np.zeros_like(positions) + + accelerations = compute_accelerations(accelerations, masses, positions) + + time = 0.0 + energy0, _, _ = compute_energies(masses, positions, velocities) + energy_previous = energy0 + + for step in range(nb_steps): + positions += time_step*velocities + 0.5*accelerations*time_step**2 + + accelerations, accelerations1 = accelerations1, accelerations + accelerations.fill(0) + accelerations = compute_accelerations(accelerations, masses, positions) + + velocities += 0.5*time_step*(accelerations+accelerations1) + + time += time_step + + if not step % 100: + energy, _, _ = compute_energies(masses, positions, velocities) + energy_previous = energy + + return energy, energy0 + +def compute_energies(masses, positions, velocities): + + ke = 0.5 * masses @ np.sum(velocities**2, axis=1) + + nb_particles = masses.size + pe = 0.0 + for index_p0 in range(nb_particles - 1): + + mass0 = masses[index_p0] + for index_p1 in range(index_p0 + 1, nb_particles): + + mass1 = masses[index_p1] + vector = positions[index_p0] - positions[index_p1] + + distance = math.sqrt((vector**2).sum()) + + pe -= (mass0*mass1) / distance**2 + + return ke+pe, ke, pe + +if __name__ == "__main__": + + try: + time_end = float(sys.argv[2]) + except IndexError: + time_end = 10.0 + + time_step = 0.001 + nb_steps = int(time_end/time_step) + 1 + + path_input = sys.argv[1] + masses, positions, velocities = load_input_data(path_input) + + print('time taken:', timeit.timeit('loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50)) diff --git a/nbody_problem/python/numpy_python.py b/nbody_problem/python/numpy_python.py new file mode 100644 index 0000000..41bd461 --- /dev/null +++ b/nbody_problem/python/numpy_python.py @@ -0,0 +1,105 @@ +import sys +import math +import timeit + +import numpy as np +import pandas as pd + +def load_input_data(path): + + df = pd.read_csv( + path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+" + ) + + masses = df["mass"].values.copy() + positions = df.loc[:, ["x", "y", "z"]].values.copy() + velocities = df.loc[:, ["vx", "vy", "vz"]].values.copy() + + return masses, positions, velocities + +def compute_accelerations(accelerations, masses, positions): + nb_particles = masses.size + + for index_p0 in range(nb_particles - 1): + position0 = positions[index_p0] + mass0 = masses[index_p0] + + vectors = position0 - positions[index_p0 + 1: nb_particles] + + distances = (vectors**2).sum(axis=1) + coefs = 1./distances**1.5 + + accelerations[index_p0] += np.sum((masses[index_p0 + 1: nb_particles] * -1 * vectors.T * coefs).T, axis=0) + accelerations[index_p0 + 1: nb_particles] += (mass0 * vectors.T * coefs).T + + return accelerations + +def numpy_loop( + time_step: float, + nb_steps: int, + masses: "float[]", + positions: "float[:,:]", + velocities: "float[:,:]", + ): + + accelerations = np.zeros_like(positions) + accelerations1 = np.zeros_like(positions) + + accelerations = compute_accelerations(accelerations, masses, positions) + + time = 0.0 + + energy0, _, _ = compute_energies(masses, positions, velocities) + energy_previous = energy0 + + for step in range(nb_steps): + positions += time_step*velocities + 0.5*accelerations*time_step**2 + + accelerations, accelerations1 = accelerations1, accelerations + accelerations.fill(0) + accelerations = compute_accelerations(accelerations, masses, positions) + + velocities += 0.5*time_step*(accelerations+accelerations1) + + time += time_step + + if not step%100: + energy, _, _ = compute_energies(masses, positions, velocities) + energy_previous = energy + + return energy, energy0 + +def compute_energies(masses, positions, velocities): + + ke = 0.5 * masses @ np.sum(velocities**2, axis=1) + + nb_particles = masses.size + pe = 0.0 + for index_p0 in range(nb_particles - 1): + + mass0 = masses[index_p0] + for index_p1 in range(index_p0 + 1, nb_particles): + + mass1 = masses[index_p1] + vector = positions[index_p0] - positions[index_p1] + + distance = math.sqrt((vector**2).sum()) + + pe -= (mass0*mass1) / distance**2 + + return ke+pe, ke, pe + +if __name__ == "__main__": + + try: + time_end = float(sys.argv[2]) + except IndexError: + time_end = 10.0 + + time_step = 0.001 + nb_steps = int(time_end/time_step) + 1 + + path_input = sys.argv[1] + masses, positions, velocities = load_input_data(path_input) + + print('time taken:', timeit.timeit('numpy_loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50)) diff --git a/nbody_problem/python/plot.py b/nbody_problem/python/plot.py new file mode 100644 index 0000000..7f67de3 --- /dev/null +++ b/nbody_problem/python/plot.py @@ -0,0 +1,74 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +def plot(x, labels, list_df, names): + + plt.figure(figsize=(350, 220), dpi=4) + plt.subplot(2, 1, 1) + + width = 0.25 + rect = [0, 0, 0] + colors = ['#4D77CF', 'navy'] + + list1 = [] + for ind, list in enumerate(list_df): + if ind < 2: + list1.append(list) + + for ind, time_taken in enumerate(list1): + rect[ind] = plt.bar(x + width*ind, time_taken, width, align='center', color=colors[ind % len(colors)], label=labels[ind]) + plt.bar_label(rect[ind], padding=3, fontsize=250) + + plt.xticks(x + width/2, names, fontsize=300) + plt.tick_params(axis='both', which='major', pad=60) + plt.yticks(fontsize=300) + plt.legend(fontsize=250) + + plt.ylabel(r'$\frac{Time\ (in\ secs)}{nParticles^{2}}$', fontsize=350) + plt.xlabel(r"$Number\ of\ Particles\ Simulated(nParticles)$", fontsize=270) + plt.title(r"$Library\ based\ Implementation$", fontsize=300) + + plt.subplot(2, 1, 2) + rect = [0, 0, 0] + colors = ['#044F88', '#4DABCF', '#013243'] + + list2 = [] + for ind, list in enumerate(list_df): + if ind >= 2: + list2.append(list) + + j = 2 + for ind, time_taken in enumerate(list2): + rect[ind] = plt.bar(x + width*ind, time_taken, width, align='center', color=colors[ind % len(colors)], label=labels[j]) + j += 1 + plt.bar_label(rect[ind], padding=3, fontsize=250) + + plt.xticks(x + width, names, fontsize=300) + plt.yticks(fontsize=300) + plt.legend(fontsize=250) + + plt.ylabel(r'$\frac{Time\ (in\ secs)}{nParticles^{2}}$', fontsize=350) + plt.tick_params(axis='both', which='major', pad=60) + plt.xlabel(r"$Number\ of\ Particles\ Simulated(nParticles)$", fontsize=270) + plt.title(r"$Compiler\ based\ Implementation$", fontsize=300) + + plt.tight_layout(pad=70.0) + plt.savefig("performance_benchmarking") + +if __name__ == "__main__": + + data_path = "benchmarks/data/table_new_norm.csv" + df = pd.read_csv(data_path) + df = df.drop(['Unnamed: 0'], axis=1) + df = df.T + df.columns = ["NumPy", "Python", "C++", "Numba", "Pythran"] + + labels = df.columns + names = ["32", "64", "128", "256"] + + df = df.T + x = np.arange(len(df.columns)) + list_df = df[0:].to_numpy().tolist() + + plot(x, labels, list_df, names) diff --git a/nbody_problem/python/pure_python.py b/nbody_problem/python/pure_python.py new file mode 100644 index 0000000..0729e2c --- /dev/null +++ b/nbody_problem/python/pure_python.py @@ -0,0 +1,154 @@ +import sys +import math +import timeit + +import numpy as np +import pandas as pd + +def load_input_data(path): + + df = pd.read_csv( + path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+" + ) + + masses = df["mass"].values.copy() + positions = df.loc[:, ["x", "y", "z"]].values.copy() + velocities = df.loc[:, ["vx", "vy", "vz"]].values.copy() + + return masses, positions, velocities + +def compute_accelerations(accelerations, masses, positions): + nb_particles = len(masses) + + for index_p0 in range(nb_particles - 1): + position0 = positions[index_p0] + mass0 = masses[index_p0] + + for index_p1 in range(index_p0 + 1, nb_particles): + mass1 = masses[index_p1] + + vector = [p0 - p1 for (p0, p1) in zip(position0, positions[index_p1])] + + distance = 0 + for i in vector: + distance += i**2 + coefs = 1./distance**1.5 + + accelerations[index_p0] = [sum(i) for i in zip([vec_val * mass1 * -1 * coefs for vec_val in vector], accelerations[index_p0])] + accelerations[index_p1] = [sum(i) for i in zip([vec_val * mass0 * coefs for vec_val in vector], accelerations[index_p1])] + + return accelerations + +def loop( + time_step: float, + nb_steps: int, + masses: "float[]", + positions: "float[:,:]", + velocities: "float[:,:]", +): + + accelerations = [[0.0 for _ in range(3)] for _ in range(len(positions))] + accelerations1 = [[0.0 for _ in range(3)] for _ in range(len(positions))] + + accelerations = compute_accelerations(accelerations, masses, positions) + + time = 0.0 + energy0, _, _ = compute_energies(masses, positions, velocities) + energy_previous = energy0 + + for step in range(nb_steps): + pos_val = [] + for (vel, acc) in zip(velocities, accelerations): + pos = [] + for (v, a) in zip(vel, acc): + pos.append(v*time_step + 0.5*a*time_step**2) + pos_val.append(pos) + + positions1 = [] + for (pos_val1, pos1) in zip(pos_val, positions): + pos_1 = [] + for (pos_val0, pos0) in zip(pos_val1, pos1): + pos_1.append(pos_val0 + pos0) + positions1.append(pos_1) + positions = positions1 + + accelerations, accelerations1 = accelerations1, accelerations + accelerations = [[0.0 for acc0 in acc1] for acc1 in accelerations] + accelerations = compute_accelerations(accelerations, masses, positions) + + new_accelerations = [] + for (acc, acc1) in zip(accelerations, accelerations1): + new_acc = [] + for (a, a1) in zip(acc, acc1): + new_acc.append(a + a1) + new_accelerations.append(new_acc) + + velocities1 = [] + for (acc, vel) in zip(new_accelerations, velocities): + vel1 = [] + for (a, v) in zip(acc, vel): + vel1.append(v + a*time_step) + velocities1.append(vel1) + velocities = velocities1 + + time += time_step + + if not step % 100: + energy, _, _ = compute_energies(masses, positions, velocities) + energy_previous = energy + + return energy, energy0 + +def compute_energies(masses, positions, velocities): + + vel = [] + for ind, val in enumerate(velocities): + vel.append([j**2 for j in val]) + val = [] + for i in vel: + k = 0 + for j in i: + k += j + val.append(k) + ke_list = [0.5*i*masses[ind] for ind, i in enumerate(val)] + ke = 0.0 + for i in ke_list: + ke += i + + nb_particles = len(masses) + pe = 0.0 + for index_p0 in range(nb_particles - 1): + mass0 = masses[index_p0] + + for index_p1 in range(index_p0 + 1, nb_particles): + + mass1 = masses[index_p1] + vector = [p0 - p1 for (p0, p1) in zip(positions[index_p0], positions[index_p1])] + + dist = 0 + for i in vector: + dist += i**2 + distance = math.sqrt(dist) + + pe -= (mass0*mass1) / distance**2 + + return ke+pe, ke, pe + +if __name__ == "__main__": + + try: + time_end = float(sys.argv[2]) + except IndexError: + time_end = 10.0 + + time_step = 0.001 + nb_steps = int(time_end/time_step) + 1 + + path_input = sys.argv[1] + masses, positions, velocities = load_input_data(path_input) + + masses = masses.tolist() + positions = positions.tolist() + velocities = velocities.tolist() + + print('time taken', timeit.timeit('loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50))