Modified codes to maintain uniformity and used Python operators rather than NumPy's operators at the required places.

khushi-411 · khushi-411 · commit a496e53b8bac · 2021-09-12T01:33:51.000+05:30
diff --git a/benchmarks/cpp/main.cpp b/benchmarks/cpp/main.cpp
@@ -57,7 +57,7 @@ class Cluster : public Star {
           real apre = 1. / sqrt(RdotR * RdotR * RdotR);
 
           for (int i = 0; i != 3; ++i) {
-            si->a[i] -= sj->m * apre * rij[i];
+            si->a[i] = sj->m * -1 * apre * rij[i];
           }
         }
       }
@@ -148,23 +148,22 @@ int main(int argc, char* argv[]) {
   int k = 0;
 
   for (int i = 0; i < 50; i++) {
-  cl.acceleration();
-
-  while (t < tend) {
+    cl.acceleration();
 
-    cl.updatePositions(dt);
+    while (t < tend) {
+      cl.updatePositions(dt);
 
-    cl.acceleration();
+      cl.acceleration();
 
-    cl.updateVelocities(dt);
+      cl.updateVelocities(dt);
 
-    t += dt;
-    k += 1;
-    if (k % 100 == 0) {
-      E = cl.energies();
+      t += dt;
+      k += 1;
+      if (k % 100 == 0) {
+        E = cl.energies();
+      }
     }
   }
-  }
 
   return 0;
 }
diff --git a/benchmarks/python/acc_numba.py b/benchmarks/python/acc_numba.py
@@ -30,28 +30,12 @@ def compute_accelerations(accelerations, masses, positions):
 
         vectors = position0 - positions[index_p0 + 1: nb_particles]
 
-        distances = np.square(vectors).sum(axis = 1)
-        coefs = distances ** 1.5
-
-        _x = np.add(np.sum(
-            np.divide(
-                np.multiply(
-                    masses[index_p0 + 1: nb_particles], -1 * vectors.T
-                    ),
-                coefs).T, axis = 0),
-            accelerations[index_p0]
-            )
-        accelerations[index_p0] = _x
-
-        _temp = np.add(
-                np.divide(
-                np.multiply(
-                    mass0, vectors.T
-                    ),
-                coefs).T,
-            accelerations[index_p0 + 1: nb_particles]
-            )
-        accelerations[index_p0 + 1: nb_particles] = _temp
+        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), axis=0)
+
+        accelerations[index_p0 + 1: nb_particles] += (mass0 * vectors.T * coefs).T
 
     return accelerations
 
@@ -74,13 +58,13 @@ def numba_loop(
     energy_previous = energy0
 
     for step in range(nb_steps):
-        positions = sum(np.multiply(velocities, time_step), 0.5 * np.multiply(accelerations, time_step ** 2)) + positions
+        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 * np.multiply(time_step, np.add(accelerations, accelerations1)) + velocities
+        velocities += 0.5*time_step*(accelerations+accelerations1)
 
         time += time_step
 
@@ -93,24 +77,23 @@ def numba_loop(
 @jit
 def compute_energies(masses, positions, velocities):
 
-    ke = 0.5 * (np.multiply(masses, np.square(velocities).sum(axis = 1)).sum())
+    ke = 0.5 * np.sum(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 = np.sqrt((vector ** 2).sum())
+            distance = math.sqrt((vector**2).sum())
 
-            pe = np.subtract(np.divide(np.multiply(mass0, mass1), np.square(distance)), pe)
+            pe -= (mass0*mass1) / distance**2
 
-    return ke + pe, ke, pe
+    return ke+pe, ke, pe
 
 if __name__ == "__main__":
 
@@ -125,4 +108,4 @@ def compute_energies(masses, positions, velocities):
     path_input = sys.argv[1]
     masses, positions, velocities = load_input_data(path_input)
 
-    print('time taken:', timeit.timeit('numba_loop(time_step, nb_steps, masses, positions, velocities)', globals = globals(), number = 50))
+    print('time taken:', timeit.timeit('numba_loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50))
diff --git a/benchmarks/python/acc_pythran.py b/benchmarks/python/acc_pythran.py
@@ -26,33 +26,23 @@ def compute_accelerations(accelerations, masses, positions):
         position0 = positions[index_p0]
         mass0 = masses[index_p0]
 
-        vector = np.empty(3)
-
         for index_p1 in range(index_p0 + 1, nb_particles):
             mass1 = masses[index_p1]
-            position1 = positions[index_p1]
 
-            for i in range(3):
-                vector[i] = position0[i] - position1[i]
+            vectors = position0 - positions[index_p1]
 
-            distance = sum(vector ** 2) * math.sqrt(sum(vector ** 2))
-            coef_m1 = mass0 / distance
-            coef_m2 = mass1 / distance
+            distances = (vectors**2).sum()
+            coefs = 1./distances**1.5
 
-            for i in range(3):
-                accelerations[index_p0][i] += coef_m1 * -1 * vector[i]
-                accelerations[index_p1][i] += coef_m2 * vector[i]
+            accelerations[index_p0] += mass1 * -1 * vectors * coefs
+            accelerations[index_p1] += mass0 * vectors * coefs
 
     return accelerations
 
 @jit
 def pythran_loop(
-        time_step: float,
-        nb_steps: int,
-        masses: "float[]",
-        positions: "float[:,:]",
-        velocities: "float[:,:]",
-    ):
+        time_step, nb_steps, masses, positions,
+        velocities):
 
     accelerations = np.zeros_like(positions)
     accelerations1 = np.zeros_like(positions)
@@ -64,13 +54,13 @@ def pythran_loop(
     energy_previous = energy0
 
     for step in range(nb_steps):
-        positions = sum(np.multiply(velocities, time_step), 0.5 *  np.multiply(accelerations, time_step ** 2)) + positions
+        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 * np.multiply(time_step, np.add(accelerations, accelerations1)) + velocities
+        velocities += 0.5*time_step*(accelerations+accelerations1)
 
         time += time_step
 
@@ -82,24 +72,23 @@ def pythran_loop(
 
 def compute_energies(masses, positions, velocities):
 
-    ke = 0.5 * (np.multiply(masses, np.square(velocities).sum(axis = 1)).sum())
+    ke = 0.5 * np.sum(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]
 
+        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 = np.sqrt((vector ** 2).sum())
+            distance = math.sqrt((vector**2).sum())
 
-            pe = np.subtract(np.divide(np.multiply(mass0, mass1), np.square(distance)), pe)
+            pe -= (mass0*mass1) / distance**2
 
-    return ke + pe, ke, pe
+    return ke+pe, ke, pe
 
 if __name__ == "__main__":
 
@@ -114,4 +103,4 @@ def compute_energies(masses, positions, velocities):
     path_input = sys.argv[1]
     masses, positions, velocities = load_input_data(path_input)
 
-    print('time taken:', timeit.timeit('pythran_loop(time_step, nb_steps, masses, positions, velocities)', globals = globals(), number = 50))
+    print('time taken:', timeit.timeit('pythran_loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50))
diff --git a/benchmarks/python/optimized-numpy.py b/benchmarks/python/optimized-numpy.py
@@ -8,7 +8,7 @@
 def load_input_data(path):
 
     df = pd.read_csv(
-        path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter = r"\s+"
+        path, names = ["mass", "x", "y", "z", "vx", "vy", "vz"], delimiter=r"\s+"
     )
 
     masses = df["mass"].values.copy()
@@ -26,28 +26,11 @@ def compute_accelerations(accelerations, masses, positions):
 
         vectors = position0 - positions[index_p0 + 1: nb_particles]
 
-        distances = np.square(vectors).sum(axis = 1)
-        coefs = distances ** 1.5
-
-        _x = np.add(np.sum(
-            np.divide(
-                np.multiply(
-                    masses[index_p0 + 1: nb_particles], -1 * vectors.T
-                    ),
-                coefs).T, axis = 0),
-            accelerations[index_p0]
-            )
-        accelerations[index_p0] = _x
-
-        _temp = np.add(
-                np.divide(
-                np.multiply(
-                    mass0, vectors.T
-                    ),
-                coefs).T,
-            accelerations[index_p0 + 1: nb_particles]
-            )
-        accelerations[index_p0 + 1: nb_particles] = _temp
+        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), axis=0)
+        accelerations[index_p0 + 1: nb_particles] += (mass0 * vectors.T * coefs).T
 
     return accelerations
 
@@ -70,28 +53,27 @@ def numpy_loop(
     energy_previous = energy0
 
     for step in range(nb_steps):
-        positions = sum(np.multiply(velocities, time_step), 0.5 * np.multiply(accelerations, time_step ** 2)) + positions
+        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 * np.multiply(time_step, np.add(accelerations, accelerations1)) + velocities
+        velocities += 0.5*time_step*(accelerations+accelerations1)
 
         time += time_step
 
-        if not step % 100:
+        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 * (np.multiply(masses, np.square(velocities).sum(axis = 1)).sum())
+    ke = 0.5 * np.sum(masses * np.sum(velocities**2, axis=1))
 
     nb_particles = masses.size
-
     pe = 0.0
     for index_p0 in range(nb_particles - 1):
 
@@ -101,11 +83,11 @@ def compute_energies(masses, positions, velocities):
             mass1 = masses[index_p1]
             vector = positions[index_p0] - positions[index_p1]
 
-            distance = np.sqrt((vector ** 2).sum())
+            distance = math.sqrt((vector**2).sum())
 
-            pe = np.subtract(np.divide(np.multiply(mass0, mass1), np.square(distance)), pe)
+            pe -= (mass0*mass1) / distance**2
 
-    return ke + pe, ke, pe
+    return ke+pe, ke, pe
 
 if __name__ == "__main__":
 
@@ -120,4 +102,4 @@ def compute_energies(masses, positions, velocities):
     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))
+    print('time taken:', timeit.timeit('numpy_loop(time_step, nb_steps, masses, positions, velocities)', globals=globals(), number=50))
diff --git a/benchmarks/python/pure-python.py b/benchmarks/python/pure-python.py
