🐛 Fixed bug in three-state simulation and duplicate cds results for QuickExact.

docs/algorithms/sidb_simulation.rst

@@ -81,6 +81,11 @@ Energy Calculation
 .. doxygenfunction:: fiction::is_ground_state
 
 
+**Header:** ``fiction/algorithms/simulation/sidb/check_simulation_results_for_equivalence.hpp``
+
+.. doxygenfunction:: fiction::check_simulation_results_for_equivalence
+
+
 Temperature Behavior
 ####################
 

docs/utils/utils.rst

@@ -196,7 +196,7 @@ Math Utils
 .. doxygenfunction:: fiction::round_to_n_decimal_places
 .. doxygenfunction:: fiction::integral_abs
 .. doxygenfunction:: fiction::binomial_coefficient
-
+.. doxygenfunction:: fiction::determine_all_combinations_of_distributing_k_entities_on_n_positions
 
 ``phmap``
 ---------

experiments/equivalence_checking_exact_simulation/equivalence_checking_exact_simulation.cpp

@@ -0,0 +1,56 @@
+//
+// Created by Jan Drewniok on 03.03.24.
+//
+
+#include <fiction/algorithms/simulation/sidb/check_simulation_results_for_equivalence.hpp>
+#include <fiction/algorithms/simulation/sidb/exhaustive_ground_state_simulation.hpp>
+#include <fiction/algorithms/simulation/sidb/quickexact.hpp>
+#include <fiction/algorithms/simulation/sidb/sidb_simulation_parameters.hpp>
+#include <fiction/layouts/coordinates.hpp>
+#include <fiction/technology/cell_technologies.hpp>
+#include <fiction/types.hpp>
+#include <fiction/utils/layout_utils.hpp>
+#include <fiction/utils/math_utils.hpp>
+
+#include <fmt/format.h>
+
+#include <cstdint>
+#include <iostream>
+
+using namespace fiction;
+
+// This script verifies the equivalence between ExGS and QuickExact. It generates all layouts consisting of
+// 4 SiDBs within an 11x11 spanned area. The simulation is then executed using both simulators, and the results
+// are compared.
+
+int main()  // NOLINT
+{
+    const auto all_cells_in_region = all_coordinates_in_spanned_area<offset::ucoord_t>({0, 0}, {10, 10});
+
+    const auto all_distributions =
+        determine_all_combinations_of_distributing_k_entities_on_n_positions(4, all_cells_in_region.size());
+
+    const auto params = sidb_simulation_parameters{3, -0.32};
+
+    uint64_t non_equivalence_counter = 0;
+
+    for (const auto& distribution : all_distributions)
+    {
+        sidb_cell_clk_lyt lyt{};
+        for (const auto idx : distribution)
+        {
+            lyt.assign_cell_type(all_cells_in_region[idx], sidb_technology::cell_type::NORMAL);
+        }
+        auto result_exgs = exhaustive_ground_state_simulation(lyt, params);
+        auto result_quickexact =
+            quickexact(lyt, quickexact_params<sidb_cell_clk_lyt>{
+                                params, quickexact_params<sidb_cell_clk_lyt>::automatic_base_number_detection::OFF});
+        if (!check_simulation_results_for_equivalence(result_exgs, result_quickexact))
+        {
+            non_equivalence_counter++;
+        }
+    }
+    std::cout << fmt::format("non equivalent layouts = {}", non_equivalence_counter) << '\n';
+
+    return EXIT_SUCCESS;
+}

include/fiction/algorithms/physical_design/design_sidb_gates.hpp

@@ -25,13 +25,12 @@
 #include <cstdint>
 #include <cstdlib>
 #include <future>
+#include <mutex>
 #include <numeric>
 #include <thread>
 #include <utility>
 #include <vector>
 
-#include <combinations.h>
-
 namespace fiction
 {
 
@@ -113,7 +112,8 @@ class design_sidb_gates_impl
     {
         const is_operational_params params_is_operational{params.phys_params, params.sim_engine};
 
-        const auto all_combinations = determine_all_combinations_of_given_sidbs_in_canvas();
+        const auto all_combinations = determine_all_combinations_of_distributing_k_entities_on_n_positions(
+            params.number_of_sidbs, static_cast<std::size_t>(all_sidbs_in_canvas.size()));
 
         std::vector<Lyt> designed_gate_layouts = {};
 
@@ -228,43 +228,6 @@ class design_sidb_gates_impl
      * All cells within the canvas.
      */
     const std::vector<typename Lyt::cell> all_sidbs_in_canvas;
-    /**
-     * Calculates all possible combinations of distributing the given number of SiDBs within a canvas
-     * based on the provided parameters. It generates combinations of SiDB indices (representing the cell position in
-     * the canvas from top to bottom and left to right) and stores them in the `all_combinations` vector. The number of
-     * SiDBs in each combination is determined by `parameter.number_of_sidbs`.
-     *
-     * @return All possible combinations as a vector of vectors of indices.
-     */
-    [[nodiscard]] std::vector<std::vector<std::size_t>> determine_all_combinations_of_given_sidbs_in_canvas() noexcept
-    {
-        std::vector<std::vector<std::size_t>> all_combinations{};
-        all_combinations.reserve(binomial_coefficient(all_sidbs_in_canvas.size(), params.number_of_sidbs));
-
-        std::vector<std::size_t> numbers(all_sidbs_in_canvas.size());
-        std::iota(numbers.begin(), numbers.end(), 0);
-
-        combinations::for_each_combination(
-            numbers.begin(),
-            numbers.begin() + static_cast<std::vector<std::size_t>::difference_type>(params.number_of_sidbs),
-            numbers.end(),
-            [this, &all_combinations](const auto begin, const auto end)
-            {
-                std::vector<std::size_t> combination{};
-                combination.reserve(params.number_of_sidbs);
-
-                for (auto it = begin; it != end; ++it)
-                {
-                    combination.push_back(static_cast<std::size_t>(*it));
-                }
-
-                all_combinations.push_back(combination);
-
-                return false;  // keep looping
-            });
-
-        return all_combinations;
-    }
     /**
      * Checks if any SiDBs within the specified cell indices are located too closely together, with a distance of less
      * than 0.5 nanometers.

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

@@ -0,0 +1,112 @@
+//
+// Created by Jan Drewniok on 03.03.24.
+//
+
+#ifndef FICTION_CHECK_SIMULATION_RESULTS_FOR_EQUIVALENCE_HPP
+#define FICTION_CHECK_SIMULATION_RESULTS_FOR_EQUIVALENCE_HPP
+
+#include "fiction/algorithms/simulation/sidb/sidb_simulation_result.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <set>
+
+namespace fiction
+{
+
+/**
+ * This function compares two SiDB simulation results for equivalence. Two results are considered equivalent if they
+ * have the same number of charge distributions and if each corresponding charge distribution has the same electrostatic
+ * potential energy and charge states for all cells.
+ *
+ * @tparam Lyt The SiDB cell-level layout type used in the simulation results.
+ * @param result1 The first SiDB simulation result to compare.
+ * @param result2 The second SiDB simulation result to compare.
+ * @return `true` if the two simulation results are equivalent, `false` otherwise.
+ */
+template <typename Lyt>
+[[nodiscard]] bool check_simulation_results_for_equivalence(const sidb_simulation_result<Lyt>& result1,
+                                                            const sidb_simulation_result<Lyt>& result2)
+{
+    auto result1_copy = sidb_simulation_result{result1};
+    auto result2_copy = sidb_simulation_result{result2};
+
+    if (result1_copy.charge_distributions.size() != result2_copy.charge_distributions.size())
+    {
+        return false;
+    }
+
+    if (!result1_copy.charge_distributions.empty() && !result2_copy.charge_distributions.empty())
+    {
+        if (result1_copy.charge_distributions.front().get_sidb_order().size() !=
+            result2_copy.charge_distributions.front().get_sidb_order().size())
+        {
+            return false;
+        }
+    }
+
+    std::set<uint64_t> unique_charge_indices1;
+    for (const auto& cds1 : result1_copy.charge_distributions)
+    {
+        unique_charge_indices1.insert(cds1.get_charge_index_and_base().first);
+    }
+
+    // check if all charge indices are unique
+    if (unique_charge_indices1.size() != result1_copy.charge_distributions.size())
+    {
+        return false;
+    }
+
+    std::set<uint64_t> unique_charge_indices2;
+    for (const auto& cds2 : result2_copy.charge_distributions)
+    {
+        unique_charge_indices2.insert(cds2.get_charge_index_and_base().first);
+    }
+
+    // check if all charge indices are unique
+    if (unique_charge_indices2.size() != result2_copy.charge_distributions.size())
+    {
+        return false;
+    }
+
+    std::sort(result1_copy.charge_distributions.begin(), result1_copy.charge_distributions.end(),
+              [](const auto& lhs, const auto& rhs)
+              { return lhs.get_charge_index_and_base().first < rhs.get_charge_index_and_base().first; });
+
+    std::sort(result2_copy.charge_distributions.begin(), result2_copy.charge_distributions.end(),
+              [](const auto& lhs, const auto& rhs)
+              { return lhs.get_charge_index_and_base().first < rhs.get_charge_index_and_base().first; });
+
+    for (auto i = 0u; i < result1_copy.charge_distributions.size(); i++)
+    {
+        const auto& cds1 = result1_copy.charge_distributions.at(i);
+        const auto& cds2 = result2_copy.charge_distributions.at(i);
+        if (std::abs(cds1.get_system_energy() - cds2.get_system_energy()) > std::numeric_limits<double>::epsilon())
+        {
+            return false;
+        }
+
+        bool charge_states_equal = true;
+        cds1.foreach_cell(
+            [&cds1, &cds2, &charge_states_equal](const auto& c)
+            {
+                if (cds1.get_charge_state(c) != cds2.get_charge_state(c))
+                {
+                    charge_states_equal = false;
+                }
+            });
+
+        if (!charge_states_equal)
+        {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+}  // namespace fiction
+
+#endif  // FICTION_CHECK_SIMULATION_RESULTS_FOR_EQUIVALENCE_HPP

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

@@ -395,6 +395,7 @@ class quickexact_impl
 
                     charge_lyt_copy.update_after_charge_change();
                     charge_lyt_copy.recompute_system_energy();
+                    charge_lyt_copy.charge_distribution_to_index_general();
                     result.charge_distributions.push_back(charge_lyt_copy);
                 }
 
@@ -416,6 +417,7 @@ class quickexact_impl
 
                 charge_lyt_copy.update_after_charge_change();
                 charge_lyt_copy.recompute_system_energy();
+                charge_lyt_copy.charge_distribution_to_index_general();
                 result.charge_distributions.push_back(charge_lyt_copy);
             }
 
@@ -426,7 +428,7 @@ class quickexact_impl
 
             charge_layout.increase_charge_index_by_one(
                 dependent_cell_mode::VARIABLE, energy_calculation::KEEP_OLD_ENERGY_VALUE,
-                charge_distribution_history::CONSIDER,
+                charge_distribution_history::NEGLECT,
                 exhaustive_sidb_simulation_engine::QUICKEXACT);  // "false" allows that the charge state of the
                                                                  // dependent cell is automatically changed based on the
                                                                  // new charge distribution.