🐛 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_equivalence_simulation_result.hpp``
+
+.. doxygenfunction:: fiction::check_equivalence_simulation_result
+
+
 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.cpp

@@ -0,0 +1,49 @@
+//
+// Created by Jan Drewniok on 03.03.24.
+//
+
+#include <fiction/algorithms/simulation/sidb/check_equivalence_simulation_result.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>
+
+using namespace fiction;
+
+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_equivalence_simulation_result(result_exgs, result_quickexact))
+        {
+            non_equivalence_counter++;
+        }
+    }
+    std::cout << fmt::format("non equivalent layouts = {}", non_equivalence_counter) << std::endl;
+}

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, 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_equivalence_simulation_result.hpp

@@ -0,0 +1,89 @@
+//
+// Created by Jan Drewniok on 03.03.24.
+//
+
+#ifndef FICTION_CHECK_EQUAL_CDS_RESULT_HPP
+#define FICTION_CHECK_EQUAL_CDS_RESULT_HPP
+
+#include "fiction/algorithms/simulation/sidb/sidb_simulation_result.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+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 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_equivalence_simulation_result(sidb_simulation_result<Lyt>& result1,
+                                                       sidb_simulation_result<Lyt>& result2)
+{
+    if (result1.charge_distributions.size() != result2.charge_distributions.size())
+    {
+        return false;
+    }
+
+    std::sort(result1.charge_distributions.begin(), result1.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.charge_distributions.begin(), result2.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.charge_distributions.size(); i++)
+    {
+        const auto& cds1 = result1.charge_distributions[i];
+        const auto& cds2 = result2.charge_distributions[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;
+        }
+
+        charge_states_equal = true;
+        cds2.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_EQUAL_CDS_RESULT_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.

include/fiction/technology/charge_distribution_surface.hpp

@@ -28,6 +28,7 @@
 #include <random>
 #include <type_traits>
 #include <utility>
+#include <vector>
 
 namespace fiction
 {
@@ -966,7 +967,7 @@ class charge_distribution_surface<Lyt, false> : public Lyt
                           ((strg->cell_charge[for_loop_counter] == sidb_charge_state::POSITIVE) &&
                            (-it + mu_p > -physical_constants::POP_STABILITY_ERR)) ||
                           ((strg->cell_charge[for_loop_counter] == sidb_charge_state::NEUTRAL) &&
-                           (-it + strg->phys_params.mu_minus > -physical_constants::POP_STABILITY_ERR) &&
+                           (-it + strg->phys_params.mu_minus > physical_constants::POP_STABILITY_ERR) &&
                            (-it + mu_p < physical_constants::POP_STABILITY_ERR)));
             for_loop_counter += 1;
             if (!valid)
@@ -1604,7 +1605,7 @@ class charge_distribution_surface<Lyt, false> : public Lyt
             else if ((loc_pot_cell + strg->phys_params.mu_plus()) > -physical_constants::POP_STABILITY_ERR)
             {
                 // dependent-cell can only be positively charged when the base number is set to three state simulation.
-                if (strg->charge_index_and_base.second == 3)
+                if (strg->charge_index_and_base.second == 3 || strg->three_state_cells.size() != 0)
                 {
                     if (strg->cell_charge[strg->dependent_cell_index] != sidb_charge_state::POSITIVE)
                     {
@@ -1946,6 +1947,22 @@ class charge_distribution_surface<Lyt, false> : public Lyt
         strg->cell_history = {};
         strg->cell_history.reserve(this->num_cells());
 
+        if (strg->charge_index_sublayout == 0)
+        {
+            for (const auto& cell : strg->three_state_cells)
+            {
+                this->assign_charge_state(cell, sidb_charge_state::NEGATIVE, false);
+            }
+        }
+
+        if (strg->charge_index_and_base.first == 0)
+        {
+            for (const auto& cell : strg->sidb_order_without_three_state_cells)
+            {
+                this->assign_charge_state(cell, sidb_charge_state::NEGATIVE, false);
+            }
+        }
+
         auto       charge_quot_positive = strg->charge_index_sublayout;
         const auto base_positive        = 3;
         auto       counter              = strg->three_state_cells.size() - 1;

include/fiction/utils/layout_utils.hpp

@@ -8,6 +8,7 @@
 #include "fiction/layouts/cell_level_layout.hpp"
 #include "fiction/layouts/coordinates.hpp"
 #include "fiction/technology/cell_ports.hpp"
+#include "fiction/technology/charge_distribution_surface.hpp"
 #include "fiction/traits.hpp"
 #include "fiction/types.hpp"
 
@@ -18,6 +19,7 @@
 #include <random>
 #include <type_traits>
 #include <utility>
+#include <vector>
 
 namespace fiction
 {

include/fiction/utils/math_utils.hpp

@@ -8,9 +8,12 @@
 #include <cmath>
 #include <cstdint>
 #include <cstdlib>
+#include <numeric>
 #include <type_traits>
 #include <vector>
 
+#include <combinations.h>
+
 namespace fiction
 {
 
@@ -76,6 +79,45 @@ T integral_abs(const T n) noexcept
     return result;
 }
 
+/**
+ * This function generates all possible combinations of distributing k entities
+ * onto n positions. Each combination is represented as a vector of indices,
+ * where each index indicates the position of an entity.
+ *
+ * @param k The number of entities to distribute.
+ * @param n The number of positions available for distribution.
+ * @return A vector of vectors representing all possible combinations of
+ *         distributing k entities on n positions.
+ */
+[[nodiscard]] std::vector<std::vector<std::size_t>>
+determine_all_combinations_of_distributing_k_entities_on_n_positions(const uint64_t k, const uint64_t n) noexcept
+{
+    std::vector<std::vector<std::size_t>> all_combinations{};
+    all_combinations.reserve(binomial_coefficient(n, k));
+
+    std::vector<std::size_t> numbers(n);
+    std::iota(numbers.begin(), numbers.end(), 0);
+
+    combinations::for_each_combination(
+        numbers.begin(), numbers.begin() + static_cast<std::vector<std::size_t>::difference_type>(k), numbers.end(),
+        [&k, &all_combinations](const auto begin, const auto end)
+        {
+            std::vector<std::size_t> combination{};
+            combination.reserve(k);
+
+            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;
+}
+
 }  // namespace fiction
 
 #endif  // FICTION_MATH_UTILS_HPP