added some performance testing

Author: simongravelle

integration_tests/test_monte_carlo_move.py

@@ -1,6 +1,8 @@
 import sys
 import os
 import time
+import pstats
+import cProfile
 import unittest
 import numpy as np
 from pint import UnitRegistry
@@ -31,7 +33,7 @@ def setUp(self):
         # Initialize the MonteCarlo object
         self.mc = MonteCarlo(
             ureg = ureg,
-            maximum_steps = 10000,
+            maximum_steps = 100000,
             thermo_period = 1000,
             dumping_period = 1000,
             number_atoms = [nmb_1],
@@ -46,19 +48,54 @@ def setUp(self):
             displace_mc = displace_mc,
         )
 
+    """
     def test_monte_carlo_run(self):
-        """Test if the Monte Carlo simulation runs without errors."""
+        # Test if the Monte Carlo simulation runs without errors.
         try:
             # Run the Monte Carlo simulation (this should not raise an exception)
             ti = time.time()
+
+            # Profile the run() method
             self.mc.run()
+            # self.mc.run()
             # If it runs successfully, assert True
             tf = time.time()
+
             print("Duration:", np.round(tf-ti, 2), "s")
             self.assertTrue(True)
         except Exception as e:
             # If any exception occurs, fail the test and print the error
             self.fail(f"Monte Carlo simulation failed with error: {e}")
+    """
+
+    def test_monte_carlo_run(self):
+        """Test if the Monte Carlo simulation runs without errors."""
+        profiler = cProfile.Profile()
+        try:
+            # Start profiling before running the Monte Carlo simulation
+            profiler.enable()
+
+            # Run the Monte Carlo simulation
+            ti = time.time()
+            self.mc.run()  # Assuming self.mc is your Monte Carlo simulation object
+
+            # Stop profiling after the run
+            profiler.disable()
+
+            tf = time.time()
+            print("Duration:", np.round(tf - ti, 2), "s")
+
+            # Convert the profiler stats into a readable format
+            stats = pstats.Stats(profiler)
+            stats.strip_dirs()  # Remove extraneous directory information
+            stats.sort_stats('time')  # Sort by time spent in function
+            stats.print_stats(10)  # Print top 10 slowest functions
+
+            self.assertTrue(True)
+        except Exception as e:
+            # If any exception occurs, fail the test and print the error
+            self.fail(f"Monte Carlo simulation failed with error: {e}")
+
 
 if __name__ == "__main__":
     unittest.main()

molecular_simulation_code/monte_carlo.py

@@ -6,6 +6,7 @@
 from initialize_simulation import InitializeSimulation
 from measurements_utilities import compute_epot
 from monte_carlo_utilities import calculate_Lambda
+# from numba_test import numba_copy
 
 import warnings
 warnings.filterwarnings('ignore')
@@ -52,7 +53,9 @@ def monte_carlo_move(self):
                 self.Epot = compute_epot(self.neighbor_lists, self.atoms_positions, self.box_mda, self.cross_coefficients)
             # Make a copy of the initial atom positions and initial energy
             initial_Epot = self.Epot
-            initial_positions = copy.deepcopy(self.atoms_positions)
+            # initial_positions = copy.deepcopy(self.atoms_positions)
+            # initial_positions = np.copy(self.atoms_positions)
+            # initial_positions = numba_copy(self.atoms_positions)
             # Pick an atom id randomly
             atom_id = np.random.randint(np.sum(self.number_atoms))
             # Move the chosen atom in a random direction
@@ -62,9 +65,12 @@ def monte_carlo_move(self):
                 move = np.append(move, 0.0)  # Pad with zero for the z-component
             else:  # 3D case
                 move = (np.random.random(3) - 0.5) * self.displace_mc
+            initial_position_atom = np.copy(self.atoms_positions[atom_id])
             self.atoms_positions[atom_id] += move
+
             # Measure the potential energy of the new configuration
             trial_Epot = compute_epot(self.neighbor_lists, self.atoms_positions, self.box_mda, self.cross_coefficients)
+
             # Evaluate whether the new configuration should be kept or not
             beta =  1/self.desired_temperature
             delta_E = trial_Epot-initial_Epot
@@ -81,7 +87,8 @@ def monte_carlo_move(self):
                 self.Epot = trial_Epot
                 self.successful_move += 1
             else: # Reject new position
-                self.atoms_positions = initial_positions # Revert to initial positions
+                self.atoms_positions[atom_id] = initial_position_atom
+                # self.atoms_positions = initial_positions # Revert to initial positions
                 self.failed_move += 1
 
     def monte_carlo_swap(self):

molecular_simulation_code/numba_test.py

@@ -0,0 +1,12 @@
+import numpy as np
+from numba import njit
+
+@njit
+def numba_copy(array):
+    # Create an empty array of the same shape and type
+    new_array = np.empty_like(array)
+    # Copy elements manually
+    for i in range(array.shape[0]):
+        for j in range(array.shape[1]):
+            new_array[i, j] = array[i, j]
+    return new_array