Make it possible to freely order state and position modifiers

cppcore/include/Model.hpp

@@ -54,11 +54,6 @@ class Model {
     Shape const& get_shape() const { return shape; }
     TranslationalSymmetry const& get_symmetry() const { return symmetry; }
 
-    std::vector<SiteStateModifier> state_modifiers() const { return system_modifiers.state; }
-    std::vector<PositionModifier> position_modifiers() const { return system_modifiers.position; }
-    std::vector<OnsiteModifier> onsite_modifiers() const { return hamiltonian_modifiers.onsite; }
-    std::vector<HoppingModifier> hopping_modifiers() const { return hamiltonian_modifiers.hopping; }
-
 public: // get properties
     std::shared_ptr<System const> const& system() const;
     Hamiltonian const& hamiltonian() const;
@@ -98,7 +93,7 @@ class Model {
     TranslationalSymmetry symmetry;
     Cartesian wave_vector = {0, 0, 0};
 
-    SystemModifiers system_modifiers;
+    std::vector<SystemModifier> system_modifiers;
     HamiltonianModifiers hamiltonian_modifiers;
     HoppingGenerators hopping_generators;
 

cppcore/include/system/SystemModifiers.hpp

@@ -9,6 +9,8 @@
 
 namespace cpb {
 
+class Foundation;
+
 /**
  Modify the state (valid or invalid) of lattice sites, e.g. to create vacancies
  */
@@ -34,12 +36,49 @@ class PositionModifier {
     PositionModifier(Function const& apply) : apply(apply) {}
 };
 
-struct SystemModifiers {
-    std::vector<SiteStateModifier> state;
-    std::vector<PositionModifier> position;
+template<class T> void apply(T const&, Foundation&) {}
+void apply(SiteStateModifier const& m, Foundation& f);
+void apply(PositionModifier const& m, Foundation& f);
+
+/**
+ Polymorphic storage for system/foundation modifiers
+
+ Behaves like a common base for several classes but without actually needing
+ to inherit from anything -- a class just needs to satisfy the interface.
+ This allows us to use value semantics with polymorphic behavior.
+
+ See: "Inheritance Is The Base Class of Evil" by Sean Parent.
+ */
+class SystemModifier {
+public:
+    template<class T>
+    SystemModifier(T x) : impl(std::make_shared<Storage<T>>(std::move(x))) { }
+
+    friend void apply(SystemModifier const& x, Foundation& f) { x.impl->v_apply(f); }
+
+private:
+    struct Interface {
+        Interface() = default;
+        virtual ~Interface() = default;
+
+        Interface(Interface const&) = delete;
+        Interface(Interface&&) = delete;
+        Interface& operator=(Interface const&) = delete;
+        Interface& operator=(Interface&&) = delete;
+
+        virtual void v_apply(Foundation&) const = 0;
+    };
+
+    template<class T>
+    struct Storage : Interface {
+        Storage(T x) : data(std::move(x)) { }
+
+        void v_apply(Foundation& f) const override { apply(data, f); }
+
+        T data;
+    };
 
-    void clear();
-    bool empty() const;
+    std::shared_ptr<Interface const> impl;
 };
 
 } // namespace cpb

cppcore/src/Model.cpp

@@ -48,12 +48,12 @@ void Model::attach_lead(int direction, Shape const& shape) {
 }
 
 void Model::add(SiteStateModifier const& m) {
-    system_modifiers.state.push_back(m);
+    system_modifiers.emplace_back(m);
     clear_structure();
 }
 
 void Model::add(PositionModifier const& m) {
-    system_modifiers.position.push_back(m);
+    system_modifiers.emplace_back(m);
     clear_structure();
 }
 
@@ -127,27 +127,12 @@ std::string Model::report() {
 std::shared_ptr<System> Model::make_system() const {
     auto foundation = shape ? Foundation(lattice, shape)
                             : Foundation(lattice, primitive);
-    if (symmetry)
+    if (symmetry) {
         symmetry.apply(foundation);
+    }
 
-    if (!system_modifiers.empty()) {
-        for (auto const& site_state_modifier : system_modifiers.state) {
-            for (auto const& pair : lattice.get_sublattices()) {
-                auto slice = foundation[pair.second.unique_id];
-                site_state_modifier.apply(slice.get_states(), slice.get_positions(), pair.first);
-            }
-
-            if (site_state_modifier.min_neighbors > 0) {
-                remove_dangling(foundation, site_state_modifier.min_neighbors);
-            }
-        }
-
-        for (auto const& position_modifier : system_modifiers.position) {
-            for (auto const& pair : lattice.get_sublattices()) {
-                auto slice = foundation[pair.second.unique_id];
-                position_modifier.apply(slice.get_positions(), pair.first);
-            }
-        }
+    for (auto const& modifier : system_modifiers) {
+        apply(modifier, foundation);
     }
 
     _leads.create_attachment_area(foundation);

cppcore/src/system/SystemModifiers.cpp

@@ -1,11 +1,25 @@
 #include "system/SystemModifiers.hpp"
-using namespace cpb;
 
-void SystemModifiers::clear() {
-    state.clear();
-    position.clear();
+#include "system/Foundation.hpp"
+
+namespace cpb {
+
+void apply(SiteStateModifier const& m, Foundation& f) {
+    for (auto const& pair : f.get_lattice().get_sublattices()) {
+        auto slice = f[pair.second.unique_id];
+        m.apply(slice.get_states(), slice.get_positions(), pair.first);
+    }
+
+    if (m.min_neighbors > 0) {
+        remove_dangling(f, m.min_neighbors);
+    }
 }
 
-bool SystemModifiers::empty() const {
-    return state.empty() && position.empty();
+void apply(PositionModifier const& m, Foundation& f) {
+    for (auto const& pair : f.get_lattice().get_sublattices()) {
+        auto slice = f[pair.second.unique_id];
+        m.apply(slice.get_positions(), pair.first);
+    }
 }
+
+} // namespace cpb

cppcore/tests/test_modifiers.cpp

@@ -32,6 +32,50 @@ TEST_CASE("SitePositionModifier") {
     REQUIRE(model.system()->positions.y[1] == Approx(1));
 }
 
+TEST_CASE("State and position modifier ordering") {
+    auto model = Model(lattice::square_2atom(), Primitive(2));
+
+    auto delete_site = SiteStateModifier([](Eigen::Ref<ArrayX<bool>> state,
+                                            CartesianArrayConstRef position,
+                                            string_view sublattice) {
+        if (sublattice == "A" && position.x()[0] < 0) {
+            state[0] = false;
+        }
+    });
+
+    auto move_site = PositionModifier([](CartesianArrayRef position, string_view sublattice) {
+        if (sublattice == "A") {
+            position.x()[0] = 10;
+        }
+    });
+
+    SECTION("State before position") {
+        REQUIRE(model.system()->num_sites() == 4);
+        REQUIRE(model.system()->positions.x[0] == Approx(-1));
+        REQUIRE(model.system()->positions.x[1] == Approx(0));
+
+        model.add(delete_site);
+        model.add(move_site);
+
+        REQUIRE(model.system()->num_sites() == 3);
+        REQUIRE(model.system()->positions.x[0] == Approx(0));
+    }
+
+    SECTION("Position before state") {
+        REQUIRE(model.system()->num_sites() == 4);
+        REQUIRE(model.system()->positions.x[0] == Approx(-1));
+        REQUIRE(model.system()->positions.x[1] == Approx(0));
+
+        model.add(move_site);
+        model.add(delete_site);
+
+        REQUIRE(model.system()->num_sites() == 4);
+        REQUIRE(model.system()->positions.x[0] == Approx(10));
+        REQUIRE(model.system()->positions.x[1] == Approx(0));
+
+    }
+}
+
 struct OnsiteEnergyOp {
     template<class Array>
     void operator()(Array energy) {