⚗️ Added a runtime script for three state simulation and apply small changes.

bindings/pyfiction/include/pyfiction/algorithms/simulation/sidb/random_sidb_layout_generator.hpp

@@ -42,11 +42,15 @@ inline void random_sidb_layout_generator(pybind11::module& m)
     namespace py = pybind11;
     using namespace pybind11::literals;
 
+    // todo update docu
     py::enum_<typename fiction::generate_random_sidb_layout_params<fiction::offset::ucoord_t>::positive_charges>(
         m, "positive_charges", DOC(fiction_generate_random_sidb_layout_params_positive_charges))
-        .value("ON", fiction::generate_random_sidb_layout_params<fiction::offset::ucoord_t>::positive_charges::ALLOWED)
-        .value("OFF",
-               fiction::generate_random_sidb_layout_params<fiction::offset::ucoord_t>::positive_charges::FORBIDDEN);
+        .value("ALLOWED",
+               fiction::generate_random_sidb_layout_params<fiction::offset::ucoord_t>::positive_charges::ALLOWED)
+        .value("FORBIDDEN",
+               fiction::generate_random_sidb_layout_params<fiction::offset::ucoord_t>::positive_charges::FORBIDDEN)
+        .value("CAN_OCCUR",
+               fiction::generate_random_sidb_layout_params<fiction::offset::ucoord_t>::positive_charges::CAN_OCCUR);
 
     /**
      * Parameters.
@@ -63,9 +67,6 @@ inline void random_sidb_layout_generator(pybind11::module& m)
         .def_readwrite("positive_sidbs",
                        &fiction::generate_random_sidb_layout_params<fiction::offset::ucoord_t>::positive_sidbs,
                        DOC(fiction_generate_random_sidb_layout_params_positive_sidbs))
-        .def_readwrite("minimal_spacing",
-                       &fiction::generate_random_sidb_layout_params<fiction::offset::ucoord_t>::minimal_spacing,
-                       DOC(fiction_generate_random_sidb_layout_params_minimal_spacing))
         .def_readwrite("maximal_attempts",
                        &fiction::generate_random_sidb_layout_params<fiction::offset::ucoord_t>::maximal_attempts,
                        DOC(fiction_generate_random_sidb_layout_params_maximal_attempts))

bindings/pyfiction/include/pyfiction/pybind11_mkdoc_docstrings.hpp

@@ -295,6 +295,25 @@ Parameter ``params``:
     A cell-level layout that implements `lyt`'s gate types with
     building blocks defined in `GateLibrary`.)doc";
 
+static const char *__doc_fiction_are_layouts_identical =
+R"doc(This function checks whether the given layouts `first_lyt` and
+`second_lyt` are identical by comparing various properties such as the
+number of cells, the types of cells, defects (if applicable), and
+charge states (if applicable). The comparison is done in a detailed
+manner depending on the specific layout type.
+
+Template parameter ``Lyt``:
+    The layout type. Must be a cell-level layout.
+
+Parameter ``first_lyt``:
+    The first layout to compare.
+
+Parameter ``second_lyt``:
+    The second layout to compare.
+
+Returns:
+    `true` if the layouts are identical, `false` otherwise.)doc";
+
 static const char *__doc_fiction_area =
 R"doc(Computes the area of a given coordinate assuming its origin is (0, 0,
 0). Calculates :math:`(|x| + 1) \cdot (|y| + 1)` by default. The
@@ -10752,10 +10771,6 @@ be difficult or even impossible to find several unique (given by
 number_of_unique_generated_layouts) layouts. Therefore, this parameter
 sets a limit for the maximum number of tries.)doc";
 
-static const char *__doc_fiction_generate_random_sidb_layout_params_minimal_spacing =
-R"doc(If positively charged SiDBs should be prevented, SiDBs are not placed
-closer than this value (Euclidean distance of two cells).)doc";
-
 static const char *__doc_fiction_generate_random_sidb_layout_params_number_of_sidbs = R"doc(Number of SiDBs that are placed on the layout.)doc";
 
 static const char *__doc_fiction_generate_random_sidb_layout_params_number_of_unique_generated_layouts = R"doc(The desired number of unique layouts to be generated.)doc";
@@ -10768,6 +10783,10 @@ static const char *__doc_fiction_generate_random_sidb_layout_params_positive_cha
 R"doc(Positive charges can occur (i.e. SiDBs can be placed right next to
 each other).)doc";
 
+static const char *__doc_fiction_generate_random_sidb_layout_params_positive_charges_CAN_OCCUR =
+R"doc(Positive charges can occur, which means that the
+`can_positive_charges_occur` function returns `true`.)doc";
+
 static const char *__doc_fiction_generate_random_sidb_layout_params_positive_charges_FORBIDDEN =
 R"doc(Positive charges are not allowed to occur (i.e. SiDBs need to be
 seperated by a few lattice points).)doc";
@@ -10776,6 +10795,8 @@ static const char *__doc_fiction_generate_random_sidb_layout_params_positive_sid
 R"doc(If positively charged SiDBs should be prevented, SiDBs are not placed
 closer than the minimal_spacing.)doc";
 
+static const char *__doc_fiction_generate_random_sidb_layout_params_sim_params = R"doc(Simulation parameters.)doc";
+
 static const char *__doc_fiction_geometric_temperature_schedule =
 R"doc(A logarithmically decreasing temperature schedule. The temperature is
 altered by multiplying it with `0.99`.

experiments/random_sidb_layout_generation/random_sidb_layout_generation.cpp

@@ -4,6 +4,7 @@
 
 #include <fiction/algorithms/simulation/sidb/random_sidb_layout_generator.hpp>
 #include <fiction/io/write_sqd_layout.hpp>
+#include <fiction/technology/cell_technologies.hpp>
 #include <fiction/types.hpp>
 
 #include <fmt/format.h>
@@ -35,7 +36,7 @@ using namespace fiction;
  *   --SE_x <value>             Sets the x coordinate of the south-east cell corner in SiQAD coordinates.
  *   --SE_y <value>             Sets the y coordinate of the south-east cell corner in SiQAD coordinates.
  *   --positive_charges <type>  Specifies whether positively charged SiDBs are allowed ("ALLOWED") or forbidden
- * ("FORBIDDEN").
+ * ("FORBIDDEN") for default physical parameters.
  *   --lower <value>            Sets the number of SiDBs for the first bunch of layouts.
  *   --upper <value>            Sets the number of SiDBs for the last bunch of layouts.
  *   --num_layouts <value>      Sets the number of layouts to generate for each SiDB count.
@@ -95,10 +96,23 @@ int main(int argc, const char* argv[])  // NOLINT
     const int32_t se_y = std::stoi(options["--SE_y"]);
     // specifies whether positively charged SiDBs are allowed ("ALLOWED") or forbidden ("FORBIDDEN")
     const std::string charges_str = options["--positive_charges"];
-    // specifies whether positively charged SiDBs are allowed ("ALLOWED") or forbidden ("FORBIDDEN")
-    const generate_random_sidb_layout_params<offset::ucoord_t>::positive_charges charges =
-        (charges_str == "ALLOWED") ? generate_random_sidb_layout_params<offset::ucoord_t>::positive_charges::ALLOWED :
-                                     generate_random_sidb_layout_params<offset::ucoord_t>::positive_charges::FORBIDDEN;
+
+    // specifies whether positively charged SiDBs are allowed ("ALLOWED"), forbidden ("FORBIDDEN") or can occur
+    // ("CAN_OCCUR")
+    generate_random_sidb_layout_params<offset::ucoord_t>::positive_charges charges{};
+    if (charges_str == "ALLOWED")
+    {
+        charges = generate_random_sidb_layout_params<offset::ucoord_t>::positive_charges::ALLOWED;
+    }
+    else if (charges_str == "CAN_OCCUR")
+    {
+        charges = generate_random_sidb_layout_params<offset::ucoord_t>::positive_charges::CAN_OCCUR;
+    }
+    else
+    {
+        charges = generate_random_sidb_layout_params<offset::ucoord_t>::positive_charges::FORBIDDEN;
+    }
+
     // sets the number of SiDBs for the first bunch of layouts
     const uint64_t lower_limit = std::stoull(options["--lower"]);
     // sets the number of SiDBs for the last bunch of layouts
@@ -165,8 +179,8 @@ int main(int argc, const char* argv[])  // NOLINT
                 }
 
                 const generate_random_sidb_layout_params<offset::ucoord_t> params{
-                    {{nw_x, nw_y}, {se_x, se_y}}, number_of_placed_sidbs, charges, 2,
-                    static_cast<uint64_t>(10E6),  number_of_layouts};
+                    {{nw_x, nw_y}, {se_x, se_y}},         number_of_placed_sidbs,      charges,
+                    sidb_simulation_parameters{3, -0.32}, static_cast<uint64_t>(10E6), number_of_layouts};
                 const auto unique_lyts = generate_multiple_random_sidb_layouts(sidb_100_cell_clk_lyt{}, params);
                 for (auto i = 0u; i < unique_lyts.size(); i++)
                 {

experiments/sidb_simulation/runtime_analysis_layouts_with_possible_positive_sidbs.cpp

@@ -0,0 +1,85 @@
+//
+// Created by Jan Drewniok on 02.09.24.
+//
+
+#include "fiction/algorithms/simulation/sidb/determine_groundstate_from_simulation_results.hpp"
+#include "fiction/algorithms/simulation/sidb/exhaustive_ground_state_simulation.hpp"
+#include "fiction/algorithms/simulation/sidb/quickexact.hpp"
+#include "fiction/algorithms/simulation/sidb/random_sidb_layout_generator.hpp"
+#include "fiction/algorithms/simulation/sidb/sidb_simulation_parameters.hpp"
+#include "fiction_experiments.hpp"
+#include "mockturtle/utils/stopwatch.hpp"
+
+#include <cstdlib>
+
+using namespace fiction;
+
+// This script compares the simulation runtime of ExGS and QuickExact for randomly
+// generated layouts that may contain positively charged SiDBs (`can_positive_charges_occur` returns `true`). The number
+// of SiDBs in the layout varies from 5 to 20.
+
+int main()  // NOLINT
+{
+    using Lyt = sidb_100_cell_clk_lyt;
+
+    experiments::experiment<std::size_t, std::size_t, double, double, double> simulation_exp{
+        "Benchmark",
+        "Number SiDBs",
+        "#Instances",
+        "Exhaustive Runtime [s]",
+        "QuickExact Runtime [s]",
+        "Average #positive SiDBs of ground state"};
+
+    const sidb_simulation_parameters sim_params{3, -0.32};
+
+    const quickexact_params<cell<Lyt>> qe_params{sim_params,
+                                                 quickexact_params<cell<Lyt>>::automatic_base_number_detection::OFF};
+
+    auto random_layouts_params = generate_random_sidb_layout_params<cell<Lyt>>{
+        {{0, 0}, {10, 10}},
+        0,
+        generate_random_sidb_layout_params<cell<Lyt>>::positive_charges::CAN_OCCUR,
+        sim_params,
+        static_cast<uint64_t>(10E6),
+        10};
+
+    for (auto num_sidbs = 5u; num_sidbs < 20; num_sidbs++)
+    {
+        random_layouts_params.number_of_sidbs = num_sidbs;
+
+        const auto random_layouts = generate_multiple_random_sidb_layouts(Lyt{}, random_layouts_params);
+
+        double runtime_exhaustive = 0;
+        double runtime_quickexact = 0;
+
+        std::vector<std::size_t> number_of_positive_sidbs_of_gs_per_layout{};
+        number_of_positive_sidbs_of_gs_per_layout.reserve(random_layouts_params.number_of_unique_generated_layouts);
+
+        for (const auto& layout : random_layouts)
+        {
+            const auto exhaustive_results_layout = exhaustive_ground_state_simulation(layout, sim_params);
+
+            const auto quickexact_results_layout = quickexact(layout, qe_params);
+
+            const auto gs = determine_groundstate_from_simulation_results(exhaustive_results_layout);
+
+            // determine the number of positive SiDBs in the ground state
+            number_of_positive_sidbs_of_gs_per_layout.push_back(gs.at(0).num_positive_sidbs());
+
+            runtime_exhaustive += mockturtle::to_seconds(exhaustive_results_layout.simulation_runtime);
+            runtime_quickexact += mockturtle::to_seconds(quickexact_results_layout.simulation_runtime);
+        }
+
+        const auto average_pos_sibs_of_gs = std::accumulate(number_of_positive_sidbs_of_gs_per_layout.cbegin(),
+                                                            number_of_positive_sidbs_of_gs_per_layout.cend(), 0) /
+                                            static_cast<double>(number_of_positive_sidbs_of_gs_per_layout.size());
+
+        simulation_exp(random_layouts_params.number_of_sidbs, random_layouts.size(), runtime_exhaustive,
+                       runtime_quickexact, average_pos_sibs_of_gs);
+
+        simulation_exp.save();
+        simulation_exp.table();
+    }
+
+    return EXIT_SUCCESS;
+}

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

@@ -6,6 +6,8 @@
 #define FICTION_RANDOM_SIDB_LAYOUT_GENERATOR_HPP
 
 #include "fiction/algorithms/path_finding/distance.hpp"
+#include "fiction/algorithms/simulation/sidb/can_positive_charges_occur.hpp"
+#include "fiction/algorithms/simulation/sidb/sidb_simulation_parameters.hpp"
 #include "fiction/technology/sidb_defects.hpp"
 #include "fiction/traits.hpp"
 #include "fiction/utils/layout_utils.hpp"
@@ -37,7 +39,11 @@
         /**
          * Positive charges are not allowed to occur (i.e. SiDBs need to be seperated by a few lattice points).
          */
-        FORBIDDEN
+        FORBIDDEN,
+        /**
+         * Positive charges can occur, which means that the `can_positive_charges_occur` function returns `true`.
+         */
+        CAN_OCCUR
     };
     /**
      * Two coordinates that span the region where SiDBs may be placed (order is not important). The first coordinate is
@@ -53,10 +59,9 @@
      */
     positive_charges positive_sidbs = positive_charges::ALLOWED;
     /**
-     * If positively charged SiDBs should be prevented, SiDBs are not placed closer than this value (Euclidean distance
-     * of two cells).
+     * Simulation parameters.
      */
-    double minimal_spacing = 2.0;
+    sidb_simulation_parameters sim_params{};
     /**
      * Maximum number of steps to place the specified number of SiDBs. Example: If the area, where SiDBs can be placed,
      * is small and many SiDBs are to be placed, several tries are required to generate a layout with no positively
@@ -73,7 +78,7 @@
      * even impossible to find several unique (given by number_of_unique_generated_layouts) layouts. Therefore, this
      * parameter sets a limit for the maximum number of tries.
      */
-    uint64_t maximal_attempts_for_multiple_layouts = 10E6;
+    uint64_t maximal_attempts_for_multiple_layouts = 1'000'0000;
 };
 
 /**
@@ -110,22 +115,7 @@
     {
         // random coordinate within given area
         const auto random_coord = random_coordinate(params.coordinate_pair.first, params.coordinate_pair.second);
-        bool       next_to_neutral_defect              = false;
-        bool       constraint_violation_positive_sidbs = false;
-
-        if (params.positive_sidbs == generate_random_sidb_layout_params<coordinate<Lyt>>::positive_charges::FORBIDDEN)
-        {
-            // checks if the new coordinate is not closer than 2 cells (Euclidean distance) to an already
-            // placed SiDB
-            lyt.foreach_cell(
-                [&lyt, &random_coord, &constraint_violation_positive_sidbs, &params](const auto& c1)
-                {
-                    if (euclidean_distance<Lyt>(lyt, c1, random_coord) < params.minimal_spacing)
-                    {
-                        constraint_violation_positive_sidbs = true;
-                    }
-                });
-        }
+        bool       next_to_neutral_defect = false;
 
         if (sidbs_affected_by_defects.count(random_coord) > 0)
         {
@@ -141,13 +131,26 @@
 
         // if the constraints that no positive SiDBs occur and the cell is not yet occupied by a defect are satisfied,
         // the SiDB is added to the layout
-        if (!constraint_violation_positive_sidbs && !random_cell_is_identical_wih_defect && !next_to_neutral_defect)
+        if (!random_cell_is_identical_wih_defect && !next_to_neutral_defect)
         {
             lyt.assign_cell_type(random_coord, technology<Lyt>::cell_type::NORMAL);
+
+            if (params.positive_sidbs ==
+                    generate_random_sidb_layout_params<coordinate<Lyt>>::positive_charges::FORBIDDEN &&
+                can_positive_charges_occur(lyt, params.sim_params))
+            {
+                lyt.assign_cell_type(random_coord, technology<Lyt>::cell_type::EMPTY);
+            }
         }
         attempt_counter += 1;
     }
 
+    if (params.positive_sidbs == generate_random_sidb_layout_params<coordinate<Lyt>>::positive_charges::CAN_OCCUR &&
+        !can_positive_charges_occur(lyt, params.sim_params))
+    {
+        return generate_random_sidb_layout(lyt_skeleton, params);
+    }
+
     return lyt;
 }
 
@@ -181,30 +184,8 @@
         const auto random_lyt = generate_random_sidb_layout(lyt_skeleton, params);
 
         // indicates if a found SiDB layout is identical to an already found one.
-        bool identical_layout = false;
-        for (const auto& old_lyt : unique_lyts)
-        {
-            // checks if two layouts have an SiDB at the same position
-            uint64_t identical_cell_counter = 0;
-            old_lyt.foreach_cell(
-                [&identical_cell_counter, random_lyt](const auto& cell_old)
-                {
-                    random_lyt.foreach_cell(
-                        [&identical_cell_counter, &cell_old](const auto& cell_new)
-                        {
-                            if (cell_new == cell_old)
-                            {
-                                identical_cell_counter += 1;
-                            }
-                        });
-                });
-
-            // all cells are identical, so the new layout is a duplicate
-            if (identical_cell_counter == random_lyt.num_cells())
-            {
-                identical_layout = true;
-            }
-        }
+        const auto identical_layout = std::any_of(unique_lyts.cbegin(), unique_lyts.cend(), [&](const auto& old_lyt)
+                                                  { return are_layouts_identical(random_lyt, old_lyt); });
 
         // if the randomly generated SiDB layout is not identical to a previously generated one, it is added to the
         // collection of all unique SiDB layouts (unique_lyts)

include/fiction/technology/charge_distribution_surface.hpp

@@ -126,7 +126,7 @@ enum class charge_index_mode
     /**
      * The charge state is assigned to the cell but the old charge index is kept.
      */
-    KEEP_CHARGE_INDEX
+    KEEP_CHARGE_INDEX,
 };
 
 /**
@@ -1866,15 +1866,12 @@ class charge_distribution_surface<Lyt, false> : public Lyt
         std::sort(strg->sidb_order.begin(), strg->sidb_order.end());
         this->foreach_cell([this, &cs](const auto&) { strg->cell_charge.push_back(cs); });
 
-        assert((((this->num_cells() < 41) && (strg->simulation_parameters.base == 3)) ||
-                ((strg->simulation_parameters.base == 2) && (this->num_cells() < 64))) &&
-               "number of SiDBs is too large");
-
+        strg->max_charge_index = static_cast<uint64_t>(
+            std::pow(static_cast<double>(strg->simulation_parameters.base), this->num_cells()) - 1);
         this->charge_distribution_to_index();
+
         this->initialize_nm_distance_matrix();
         this->initialize_potential_matrix();
-        strg->max_charge_index = static_cast<uint64_t>(
-            std::pow(static_cast<double>(strg->simulation_parameters.base), this->num_cells()) - 1);
         this->update_local_potential();
         this->recompute_system_energy();
         this->validity_check();