Skip to content

Commit 520c4a6

Browse files
mattipmattip
authored
Merge pull request #1 from khushi-411/numpy_benchmarking
Implementation of the N-Body Problem.
2 parents 5ef8b15 + a9ddff9 commit 520c4a6

File tree

5 files changed

+611
-0
lines changed

5 files changed

+611
-0
lines changed

Diff for: ‎nbody_problem/cpp/main.cpp

+170
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,170 @@
1+
#include <iostream>
2+
#include <math.h>
3+
#include <numeric>
4+
#include <vector>
5+
6+
typedef double real;
7+
using namespace std;
8+
9+
class Star {
10+
public:
11+
real m;
12+
vector<real> r;
13+
vector<real> v;
14+
vector<real> a, a0;
15+
16+
Star() {
17+
r.assign(3, 0);
18+
v.assign(3, 0);
19+
a.assign(3, 0);
20+
a0.assign(3, 0);
21+
}
22+
23+
Star(real mass, vector<real> pos, vector<real> vel) {
24+
m = mass;
25+
r = pos;
26+
v = vel;
27+
}
28+
29+
friend ostream &operator<<(ostream &so, const Star &si) {
30+
so << si.m << " " << si.r[0] << " " << si.r[1] << " " << si.r[2] << " "
31+
<< si.v[0] << " " << si.v[1] << " " << si.v[2] << endl;
32+
return so;
33+
}
34+
};
35+
36+
class Cluster : public Star {
37+
protected:
38+
public:
39+
vector<Star> s;
40+
Cluster() : Star() {}
41+
42+
void acceleration() {
43+
for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si)
44+
si->a.assign(3, 0);
45+
46+
for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
47+
vector<real> rij(3);
48+
real init = 0.0;
49+
50+
for (vector<Star>::iterator sj = s.begin(); sj != s.end(); ++sj) {
51+
if (si != sj) {
52+
53+
for (int i = 0; i != 3; ++i)
54+
rij[i] = si->r[i] - sj->r[i];
55+
56+
real RdotR = inner_product(rij.begin(), rij.end(), rij.begin(), init);
57+
real apre = 1. / sqrt(RdotR * RdotR * RdotR);
58+
59+
for (int i = 0; i != 3; ++i) {
60+
si->a[i] = sj->m * -1 * apre * rij[i];
61+
}
62+
}
63+
}
64+
}
65+
}
66+
67+
void updatePositions(real dt) {
68+
for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
69+
si->a0 = si->a;
70+
for (int i = 0; i != 3; ++i)
71+
si->r[i] += dt * si->v[i] + 0.5 * dt * dt * si->a0[i];
72+
}
73+
}
74+
75+
void updateVelocities(real dt) {
76+
77+
for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
78+
for (int i = 0; i != 3; ++i)
79+
si->v[i] += 0.5 * dt * (si->a0[i] + si->a[i]);
80+
si->a0 = si->a;
81+
}
82+
}
83+
84+
vector<real> energies() {
85+
real init = 0;
86+
vector<real> E(3), rij(3);
87+
E.assign(3, 0);
88+
89+
for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si)
90+
E[1] += 0.5 * si->m *
91+
inner_product(si->v.begin(), si->v.end(), si->v.begin(), init);
92+
93+
// Potential energy
94+
for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
95+
for (vector<Star>::iterator sj = si + 1; sj != s.end(); ++sj) {
96+
for (int i = 0; i != 3; ++i)
97+
rij[i] = si->r[i] - sj->r[i];
98+
E[2] -= si->m * sj->m /
99+
sqrt(inner_product(rij.begin(), rij.end(), rij.begin(), init));
100+
}
101+
}
102+
E[0] = E[1] + E[2];
103+
return E;
104+
}
105+
106+
// Print function
107+
friend ostream &operator<<(ostream &so, Cluster &cl) {
108+
for (vector<Star>::iterator si = cl.s.begin(); si != cl.s.end(); ++si)
109+
so << *si;
110+
return so;
111+
}
112+
};
113+
114+
int main(int argc, char* argv[]) {
115+
116+
Cluster cl;
117+
real m;
118+
int dummy;
119+
vector<real> r(3), v(3);
120+
121+
// Read input data from the command line (makeplummer | dumbp)
122+
do {
123+
cin >> dummy;
124+
cin >> m;
125+
for (int i = 0; i != 3; ++i)
126+
cin >> r[i];
127+
for (int i = 0; i != 3; ++i)
128+
cin >> v[i];
129+
cl.s.push_back(Star(m, r, v));
130+
} while (!cin.eof());
131+
132+
cl.s.pop_back();
133+
134+
// Compute initial energu of the system
135+
vector<real> E(3), E0(3);
136+
E0 = cl.energies();
137+
138+
// Start time, end time and simulation step
139+
real t = 0.0;
140+
real tend;
141+
142+
if (argc > 1)
143+
tend = strtod(argv[1], NULL);
144+
else
145+
tend = 10.0;
146+
147+
real dt = 1e-3;
148+
int k = 0;
149+
150+
for (int i = 0; i < 50; i++) {
151+
cl.acceleration();
152+
153+
while (t < tend) {
154+
cl.updatePositions(dt);
155+
156+
cl.acceleration();
157+
158+
cl.updateVelocities(dt);
159+
160+
t += dt;
161+
k += 1;
162+
if (k % 100 == 0) {
163+
E = cl.energies();
164+
}
165+
}
166+
t = 0;
167+
}
168+
169+
return 0;
170+
}

Diff for: ‎nbody_problem/python/compiled_methods.py

+108
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
1+
import sys
2+
import math
3+
import timeit
4+
5+
import numpy as np
6+
import pandas as pd
7+
8+
from transonic import jit
9+
10+
def load_input_data(path):
11+
12+
df = pd.read_csv(
13+
path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+"
14+
)
15+
16+
masses = df["mass"].values.copy()
17+
positions = df.loc[:, ["x", "y", "z"]].values.copy()
18+
velocities = df.loc[:, ["vx", "vy", "vz"]].values.copy()
19+
20+
return masses, positions, velocities
21+
22+
def compute_accelerations(accelerations, masses, positions):
23+
nb_particles = masses.size
24+
25+
for index_p0 in range(nb_particles - 1):
26+
position0 = positions[index_p0]
27+
mass0 = masses[index_p0]
28+
29+
# TODO: Use compiled methods like Numba & Pythran in vectorized approach.
30+
# Issue: https://github.com/khushi-411/numpy-benchmarks/issues/4
31+
for index_p1 in range(index_p0 + 1, nb_particles):
32+
mass1 = masses[index_p1]
33+
34+
vectors = position0 - positions[index_p1]
35+
36+
distances = (vectors**2).sum()
37+
coefs = 1./distances**1.5
38+
39+
accelerations[index_p0] += mass1 * -1 * vectors * coefs
40+
accelerations[index_p1] += mass0 * vectors * coefs
41+
42+
return accelerations
43+
44+
@jit
45+
def loop(
46+
time_step, nb_steps, masses, positions,
47+
velocities):
48+
49+
accelerations = np.zeros_like(positions)
50+
accelerations1 = np.zeros_like(positions)
51+
52+
accelerations = compute_accelerations(accelerations, masses, positions)
53+
54+
time = 0.0
55+
energy0, _, _ = compute_energies(masses, positions, velocities)
56+
energy_previous = energy0
57+
58+
for step in range(nb_steps):
59+
positions += time_step*velocities + 0.5*accelerations*time_step**2
60+
61+
accelerations, accelerations1 = accelerations1, accelerations
62+
accelerations.fill(0)
63+
accelerations = compute_accelerations(accelerations, masses, positions)
64+
65+
velocities += 0.5*time_step*(accelerations+accelerations1)
66+
67+
time += time_step
68+
69+
if not step % 100:
70+
energy, _, _ = compute_energies(masses, positions, velocities)
71+
energy_previous = energy
72+
73+
return energy, energy0
74+
75+
def compute_energies(masses, positions, velocities):
76+
77+
ke = 0.5 * masses @ np.sum(velocities**2, axis=1)
78+
79+
nb_particles = masses.size
80+
pe = 0.0
81+
for index_p0 in range(nb_particles - 1):
82+
83+
mass0 = masses[index_p0]
84+
for index_p1 in range(index_p0 + 1, nb_particles):
85+
86+
mass1 = masses[index_p1]
87+
vector = positions[index_p0] - positions[index_p1]
88+
89+
distance = math.sqrt((vector**2).sum())
90+
91+
pe -= (mass0*mass1) / distance**2
92+
93+
return ke+pe, ke, pe
94+
95+
if __name__ == "__main__":
96+
97+
try:
98+
time_end = float(sys.argv[2])
99+
except IndexError:
100+
time_end = 10.0
101+
102+
time_step = 0.001
103+
nb_steps = int(time_end/time_step) + 1
104+
105+
path_input = sys.argv[1]
106+
masses, positions, velocities = load_input_data(path_input)
107+
108+
print('time taken:', timeit.timeit('loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50))

Diff for: ‎nbody_problem/python/numpy_python.py

+105
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
import sys
2+
import math
3+
import timeit
4+
5+
import numpy as np
6+
import pandas as pd
7+
8+
def load_input_data(path):
9+
10+
df = pd.read_csv(
11+
path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+"
12+
)
13+
14+
masses = df["mass"].values.copy()
15+
positions = df.loc[:, ["x", "y", "z"]].values.copy()
16+
velocities = df.loc[:, ["vx", "vy", "vz"]].values.copy()
17+
18+
return masses, positions, velocities
19+
20+
def compute_accelerations(accelerations, masses, positions):
21+
nb_particles = masses.size
22+
23+
for index_p0 in range(nb_particles - 1):
24+
position0 = positions[index_p0]
25+
mass0 = masses[index_p0]
26+
27+
vectors = position0 - positions[index_p0 + 1: nb_particles]
28+
29+
distances = (vectors**2).sum(axis=1)
30+
coefs = 1./distances**1.5
31+
32+
accelerations[index_p0] += np.sum((masses[index_p0 + 1: nb_particles] * -1 * vectors.T * coefs).T, axis=0)
33+
accelerations[index_p0 + 1: nb_particles] += (mass0 * vectors.T * coefs).T
34+
35+
return accelerations
36+
37+
def numpy_loop(
38+
time_step: float,
39+
nb_steps: int,
40+
masses: "float[]",
41+
positions: "float[:,:]",
42+
velocities: "float[:,:]",
43+
):
44+
45+
accelerations = np.zeros_like(positions)
46+
accelerations1 = np.zeros_like(positions)
47+
48+
accelerations = compute_accelerations(accelerations, masses, positions)
49+
50+
time = 0.0
51+
52+
energy0, _, _ = compute_energies(masses, positions, velocities)
53+
energy_previous = energy0
54+
55+
for step in range(nb_steps):
56+
positions += time_step*velocities + 0.5*accelerations*time_step**2
57+
58+
accelerations, accelerations1 = accelerations1, accelerations
59+
accelerations.fill(0)
60+
accelerations = compute_accelerations(accelerations, masses, positions)
61+
62+
velocities += 0.5*time_step*(accelerations+accelerations1)
63+
64+
time += time_step
65+
66+
if not step%100:
67+
energy, _, _ = compute_energies(masses, positions, velocities)
68+
energy_previous = energy
69+
70+
return energy, energy0
71+
72+
def compute_energies(masses, positions, velocities):
73+
74+
ke = 0.5 * masses @ np.sum(velocities**2, axis=1)
75+
76+
nb_particles = masses.size
77+
pe = 0.0
78+
for index_p0 in range(nb_particles - 1):
79+
80+
mass0 = masses[index_p0]
81+
for index_p1 in range(index_p0 + 1, nb_particles):
82+
83+
mass1 = masses[index_p1]
84+
vector = positions[index_p0] - positions[index_p1]
85+
86+
distance = math.sqrt((vector**2).sum())
87+
88+
pe -= (mass0*mass1) / distance**2
89+
90+
return ke+pe, ke, pe
91+
92+
if __name__ == "__main__":
93+
94+
try:
95+
time_end = float(sys.argv[2])
96+
except IndexError:
97+
time_end = 10.0
98+
99+
time_step = 0.001
100+
nb_steps = int(time_end/time_step) + 1
101+
102+
path_input = sys.argv[1]
103+
masses, positions, velocities = load_input_data(path_input)
104+
105+
print('time taken:', timeit.timeit('numpy_loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50))

0 commit comments

Comments
 (0)