✨ Added multi-threading support for QuickSim

include/fiction/algorithms/simulation/sidb/is_ground_state.hpp

@@ -16,15 +16,13 @@ namespace fiction
 {
 
 /**
- * This function checks if the ground state is found by the *quicksim* algorithm.
+ * This function checks if the ground state is found by the *QuickSim* algorithm.
  *
  * @tparam Lyt Cell-level layout type.
- * @param quicksim_results All found physically valid charge distribution surfaces obtained by the quicksim algorithm
- * (quicksim).
- * @param exhaustive_results All valid charge distribution surfaces determined by ExGS
- * (exhaustive_ground_state_simulation).
- * @return Returns `true` if the relative difference between the lowest energies of the two sets is less than 0.00001,
- * `false` otherwise.
+ * @param quicksim_results All found physically valid charge distribution surfaces obtained by the *QuickSim* algorithm.
+ * @param exhaustive_results All valid charge distribution surfaces determined by ExGS.
+ * @return Returns `true` if the relative difference between the lowest energies of the two sets is less than \f$
+ * 0.00001 \f$, `false` otherwise.
  */
 template <typename Lyt>
 [[nodiscard]] bool is_ground_state(const quicksim_stats<Lyt>& quicksim_results,

include/fiction/algorithms/simulation/sidb/quicksim.hpp

@@ -18,13 +18,15 @@
 #include <cstdint>
 #include <iostream>
 #include <limits>
+#include <mutex>
+#include <thread>
 #include <vector>
 
 namespace fiction
 {
 
 /**
- * This struct stores the parameters for the *quicksim* algorithm.
+ * This struct stores the parameters for the *QuickSim* algorithm.
  */
 struct quicksim_params
 {
@@ -37,43 +39,59 @@ struct quicksim_params
      */
     uint64_t interation_steps{80};
     /**
-     * Alpha parameter for the *quicksim* algorithm (should be reduced if no result is found).
+     * `alpha` parameter for the *QuickSim* algorithm (should be reduced if no result is found).
      */
     double alpha{0.7};
+    /**
+     * Number of threads to spawn. By default the number of threads is set to the number of available hardware threads.
+     */
+    uint64_t number_threads{std::thread::hardware_concurrency()};
 };
 
 /**
- * This struct stores the simulation runtime and all physically valid charge layouts gained by the quicksim algorithm.
+ * This struct stores the simulation runtime and all physically valid charge layouts gained by the *QuickSim* algorithm.
  *
  * @paramt Lyt Cell-level layout type.
  */
 template <typename Lyt>
 struct quicksim_stats
 {
-    mockturtle::stopwatch<>::duration             time_total{0};
+    /**
+     * Total simulation runtime.
+     */
+    mockturtle::stopwatch<>::duration time_total{0};
+    /**
+     * Vector of all physically valid charge layouts.
+     */
     std::vector<charge_distribution_surface<Lyt>> valid_lyts{};
-
+    /**
+     * Report the simulation statistics in a human-readable fashion.
+     *
+     * @param out Output stream to write to.
+     */
     void report(std::ostream& out = std::cout)
     {
-        out << fmt::format("total time  = {:.2f} secs\n", mockturtle::to_seconds(time_total));
+        out << fmt::format("[i] total runtime: {:.2f} secs\n", mockturtle::to_seconds(time_total));
+
         if (!energy_distribution<Lyt>(valid_lyts).empty())
         {
             for (auto [energy, count] : energy_distribution<Lyt>(valid_lyts))
             {
-                out << fmt::format("[i] the lowest state energy is  = {:.4f} \n", minimum_energy(valid_lyts));
-                out << fmt::format("energy: {} | occurance: {} \n", energy, count);
+                out << fmt::format("[i] lowest energy state: {:.4f} meV \n", minimum_energy(valid_lyts));
+                out << fmt::format("[i] energy: {} | occurrence: {} \n", energy, count);
             }
         }
         else
         {
             std::cout << "no state found" << std::endl;
         }
+
         std::cout << "_____________________________________________________ \n";
     }
 };
 
 /**
- * The *quicksim* algorithm is an electrostatic ground state simulation algorithm for SiDB layouts. It determines
+ * The *QuickSim* algorithm is an electrostatic ground state simulation algorithm for SiDB layouts. It determines
  * physically valid charge configurations (with minimal energy) of a given (already initialized) charge distribution
  * layout. Depending on the simulation parameters, the ground state is found with a certain probability after one run.
  *
@@ -90,6 +108,7 @@ void quicksim(const Lyt& lyt, const quicksim_params& ps = quicksim_params{}, qui
     static_assert(has_sidb_technology_v<Lyt>, "Lyt must be an SiDB layout");
 
     quicksim_stats<Lyt> st{};
+    st.valid_lyts.reserve(ps.interation_steps);
 
     // measure run time (artificial scope)
     {
@@ -100,8 +119,6 @@ void quicksim(const Lyt& lyt, const quicksim_params& ps = quicksim_params{}, qui
         // set the given physical parameters
         charge_lyt.set_physical_parameters(ps.phys_params);
 
-        std::vector<charge_distribution_surface<Lyt>> result{};
-
         charge_lyt.set_all_charge_states(sidb_charge_state::NEUTRAL);
         charge_lyt.update_after_charge_change();
 
@@ -118,30 +135,57 @@ void quicksim(const Lyt& lyt, const quicksim_params& ps = quicksim_params{}, qui
             st.valid_lyts.push_back(charge_distribution_surface<Lyt>{charge_lyt});
         }
 
-        auto       best_energy = std::numeric_limits<double>::max();
-        const auto bound       = static_cast<uint64_t>(std::round(0.6 * static_cast<double>(charge_lyt.num_cells())));
-        for (uint64_t z = 0u; z < ps.interation_steps; z++)
+        // split the iterations among threads
+        const auto iter_per_thread =
+            std::max(ps.interation_steps / ps.number_threads,
+                     uint64_t{1});  // If the number of set threads is greater than the number of iterations, the
+                                    // number of threads defines how many times QuickSim is repeated
+
+        const auto bound = static_cast<uint64_t>(std::round(0.6 * static_cast<double>(charge_lyt.num_cells())));
+
+        std::vector<std::thread> threads{};
+        threads.reserve(ps.number_threads);
+        std::mutex mutex{};  // used to control access to shared resources
+
+        for (uint64_t z = 0ul; z < ps.number_threads; z++)
         {
-            for (uint64_t i = 0u; i < bound; i++)
-            {
-                std::vector<uint64_t> index_start{i};
-                charge_lyt.set_all_charge_states(sidb_charge_state::NEUTRAL);
-                charge_lyt.assign_charge_state_by_cell_index(i, sidb_charge_state::NEGATIVE);
-                charge_lyt.update_local_potential();
-                charge_lyt.recompute_system_energy();
-
-                const auto upper_limit = static_cast<uint64_t>(static_cast<double>(charge_lyt.num_cells()) / 1.5);
-                for (uint64_t num = 0; num < upper_limit; num++)
+            threads.emplace_back(
+                [&]
                 {
-                    charge_lyt.adjacent_search(ps.alpha, index_start);
-                    charge_lyt.validity_check();
+                    charge_distribution_surface<Lyt> charge_lyt_copy{charge_lyt};
 
-                    if (charge_lyt.is_physically_valid() && (charge_lyt.get_system_energy() <= best_energy))
+                    for (uint64_t l = 0ul; l < iter_per_thread; ++l)
                     {
-                        st.valid_lyts.push_back(charge_distribution_surface<Lyt>{charge_lyt});
+                        for (uint64_t i = 0ul; i < bound; ++i)
+                        {
+                            std::vector<uint64_t> index_start{i};
+                            charge_lyt_copy.set_all_charge_states(sidb_charge_state::NEUTRAL);
+                            charge_lyt_copy.assign_charge_state_by_cell_index(i, sidb_charge_state::NEGATIVE);
+                            charge_lyt_copy.update_local_potential();
+                            charge_lyt_copy.set_system_energy_to_zero();
+
+                            const auto upper_limit =
+                                static_cast<uint64_t>(static_cast<double>(charge_lyt_copy.num_cells()) / 1.5);
+
+                            for (uint64_t num = 0ul; num < upper_limit; num++)
+                            {
+                                charge_lyt_copy.adjacent_search(ps.alpha, index_start);
+                                charge_lyt_copy.validity_check();
+
+                                if (charge_lyt_copy.is_physically_valid())
+                                {
+                                    const std::lock_guard lock{mutex};
+                                    st.valid_lyts.push_back(charge_distribution_surface<Lyt>{charge_lyt_copy});
+                                }
+                            }
+                        }
                     }
-                }
-            }
+                });
+        }
+
+        for (auto& thread : threads)
+        {
+            thread.join();
         }
     }
 

include/fiction/algorithms/simulation/sidb/time_to_solution.hpp

@@ -26,7 +26,7 @@ namespace fiction
 
 /**
  * This struct stores the time-to-solution, the simulation accuracy and the average single simulation runtime of
- * quicksim (see quicksim.hpp).
+ * *QuickSim* (see quicksim.hpp).
  *
  */
 struct time_to_solution_stats
@@ -61,15 +61,15 @@ struct time_to_solution_stats
     }
 };
 /**
- * This function determines the time-to-solution (TTS) and the accuracy (acc) of the *quicksim* algorithm.
+ * This function determines the time-to-solution (TTS) and the accuracy (acc) of the *QuickSim* algorithm.
  *
  * @tparam Lyt Cell-level layout type.
  * @param lyt Layout that is used for the simulation.
  * @param sidb_params Physical SiDB parameters which are used for the simulation.
  * @param ps Pointer to a struct where the results (time_to_solution, acc, single runtime) are stored.
- * @param repetitions Number of repetitions to determine the simulation accuracy (repetitions = 100 ==> accuracy is
- * precise to 1%).
- * @param confidence_level The time-to-solution also depends one the given confidence level which can be set here.
+ * @param repetitions Number of repetitions to determine the simulation accuracy (`repetitions = 100` means that
+ * accuracy is precise to 1%).
+ * @param confidence_level The time-to-solution also depends on the given confidence level which can be set here.
  */
 template <typename Lyt>
 void sim_acc_tts(const Lyt& lyt, const quicksim_params& quicksim_params, time_to_solution_stats* ps = nullptr,

include/fiction/technology/charge_distribution_surface.hpp

@@ -509,6 +509,13 @@ class charge_distribution_surface<Lyt, false> : public Lyt
 
         return std::nullopt;
     }
+    /**
+     * Sets the electrostatic system energy to zero. Can be used if only one SiDB is charged.
+     */
+    void set_system_energy_to_zero() noexcept
+    {
+        strg->system_energy = 0.0;
+    }
     /**
      * Calculates the system's total electrostatic potential energy and stores it in the storage.
      */
@@ -706,15 +713,16 @@ class charge_distribution_surface<Lyt, false> : public Lyt
         this->index_to_charge_distribution();
     }
     /**
-     * This function is used for the *quicksim* algorithm (see quicksim.hpp). It gets a vector with indices representing
+     * This function is used for the *QuickSim* algorithm (see quicksim.hpp). It gets a vector with indices representing
      * negatively charged SiDBs as input. Afterward, a distant and a neutrally charged SiDB is localized using a min-max
-     * diversity algorithm. This selected SiDB is set to "negativ" and the index is added to the input vector such that
+     * diversity algorithm. This selected SiDB is set to "negative" and the index is added to the input vector such that
      * the next iteration works correctly.
      *
      * @param alpha A parameter for the algorithm (default: 0.7).
      * @param negative_indices Vector of SiDBs indices that are already negatively charged (double occupied).
      */
     void adjacent_search(const double alpha, std::vector<uint64_t>& negative_indices) noexcept
+
     {
         double     dist_max     = 0;
         const auto reserve_size = this->num_cells() - negative_indices.size();

test/algorithms/simulation/sidb/quicksim.cpp

@@ -74,4 +74,103 @@ TEMPLATE_TEST_CASE(
             CHECK(!it.charge_exists(sidb_charge_state::POSITIVE));
         }
     }
+
+    SECTION("zero threads")
+    {
+        TestType lyt{{20, 10}};
+
+        lyt.assign_cell_type({1, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({3, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({4, 3, 0}, TestType::cell_type::NORMAL);
+
+        lyt.assign_cell_type({6, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({7, 3, 0}, TestType::cell_type::NORMAL);
+
+        lyt.assign_cell_type({6, 10, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({7, 10, 0}, TestType::cell_type::NORMAL);
+
+        quicksim_stats<TestType>         quicksimstats{};
+        const sidb_simulation_parameters params{2, -0.30};
+        const struct quicksim_params     quicksim_params
+        {
+            params, 80, 0.7, 0
+        };
+        quicksim<TestType>(lyt, quicksim_params, &quicksimstats);
+        CHECK(quicksimstats.valid_lyts.empty());
+    }
+
+    SECTION("one thread")
+    {
+        TestType lyt{{20, 10}};
+
+        lyt.assign_cell_type({1, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({3, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({4, 3, 0}, TestType::cell_type::NORMAL);
+
+        lyt.assign_cell_type({6, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({7, 3, 0}, TestType::cell_type::NORMAL);
+
+        lyt.assign_cell_type({6, 10, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({7, 10, 0}, TestType::cell_type::NORMAL);
+
+        quicksim_stats<TestType>         quicksimstats{};
+        const sidb_simulation_parameters params{2, -0.30};
+        const struct quicksim_params     quicksim_params
+        {
+            params, 80, 0.7, 1
+        };
+        quicksim<TestType>(lyt, quicksim_params, &quicksimstats);
+        CHECK(!quicksimstats.valid_lyts.empty());
+        CHECK(quicksimstats.time_total.count() > 0);
+    }
+
+    SECTION("two threads")
+    {
+        TestType lyt{{20, 10}};
+
+        lyt.assign_cell_type({1, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({3, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({4, 3, 0}, TestType::cell_type::NORMAL);
+
+        lyt.assign_cell_type({6, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({7, 3, 0}, TestType::cell_type::NORMAL);
+
+        lyt.assign_cell_type({6, 10, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({7, 10, 0}, TestType::cell_type::NORMAL);
+
+        quicksim_stats<TestType>         quicksimstats{};
+        const sidb_simulation_parameters params{2, -0.30};
+        const struct quicksim_params     quicksim_params
+        {
+            params, 80, 0.7, 2
+        };
+        quicksim<TestType>(lyt, quicksim_params, &quicksimstats);
+        CHECK(!quicksimstats.valid_lyts.empty());
+        CHECK(quicksimstats.time_total.count() > 0);
+    }
+
+    SECTION("100 threads")
+    {
+        TestType lyt{{20, 10}};
+
+        lyt.assign_cell_type({1, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({3, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({4, 3, 0}, TestType::cell_type::NORMAL);
+
+        lyt.assign_cell_type({6, 3, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({7, 3, 0}, TestType::cell_type::NORMAL);
+
+        lyt.assign_cell_type({6, 10, 0}, TestType::cell_type::NORMAL);
+        lyt.assign_cell_type({7, 10, 0}, TestType::cell_type::NORMAL);
+
+        quicksim_stats<TestType>         quicksimstats{};
+        const sidb_simulation_parameters params{2, -0.30};
+        const struct quicksim_params     quicksim_params
+        {
+            params, 80, 0.7, 100
+        };
+        quicksim<TestType>(lyt, quicksim_params, &quicksimstats);
+        CHECK(!quicksimstats.valid_lyts.empty());
+        CHECK(quicksimstats.time_total.count() > 0);
+    }
 }