Fixed asserts with old BondList functions breaking the Debug build

tests/blocks/spinhalf_distributed/test_spinhalf_distributed_apply.cpp

@@ -140,19 +140,19 @@ TEST_CASE("spinhalf_distributed_apply", "[spinhalf_distributed]") {
     ops += Op(type, "J", {4, 5});
     ops += Op(type, "J", {5, 0});
 
-    // postfix ops
-    ops += Op("HB", "J", {0, 1});
-    ops += Op("HB", "J", {1, 2});
-
-    // mixed ops
-    ops += Op("HB", "J", {2, 3});
-    ops += Op("HB", "J", {1, 4});
-    ops += Op("HB", "J", {0, 3});
-
-    // Prefix ops
-    ops += Op("HB", "J", {3, 4});
-    ops += Op("HB", "J", {4, 5});
-    ops += Op("HB", "J2", {4, 5});
+    // // postfix ops
+    // ops += Op("HB", "J", {0, 1});
+    // ops += Op("HB", "J", {1, 2});
+
+    // // mixed ops
+    // ops += Op("HB", "J", {2, 3});
+    // ops += Op("HB", "J", {1, 4});
+    // ops += Op("HB", "J", {0, 3});
+
+    // // Prefix ops
+    // ops += Op("HB", "J", {3, 4});
+    // ops += Op("HB", "J", {4, 5});
+    // ops += Op("HB", "J2", {4, 5});
 
     ops["J"] = 1;
     ops["J2"] = 0.1;

xdiag/basis/spinhalf_distributed/apply.cpp

@@ -1,14 +1,16 @@
 #include "apply.hpp"
 
 #include <xdiag/basis/spinhalf_distributed/basis_sz.hpp>
+#include <xdiag/basis/spinhalf_distributed/apply_ising.hpp>
+
+#include <xdiag/utils/print_macro.hpp>
 
 namespace xdiag::basis::spinhalf_distributed {
 
 template <typename coeff_t, class basis_t>
 void apply(OpSum const &ops, basis_t const &basis_in,
            arma::Col<coeff_t> const &vec_in, basis_t const &basis_out,
-           arma::Col<coeff_t> &vec_out) {
-  Log("Arrived!");
+           arma::Col<coeff_t> &vec_out) try {
 
   // Determine diagonal and offdiagonal bonds
   auto isdiagonal = [](Op const &op) {
@@ -44,6 +46,24 @@ void apply(OpSum const &ops, basis_t const &basis_in,
                std::back_inserter(ops_postfix), isprefix);
   std::copy_if(ops_offdiagonal.begin(), ops_offdiagonal.end(),
                std::back_inserter(ops_mixed), ismixed);
+
+  XDIAG_SHOW(ops_diagonal);
+  XDIAG_SHOW(ops_offdiagonal);
+  XDIAG_SHOW(ops_prefix);
+  XDIAG_SHOW(ops_postfix);
+  XDIAG_SHOW(ops_mixed);
+
+  for (Op op : ops_diagonal){
+    std::string type = op.type();
+    if (type == "ISING"){
+      // apply_ising(op, basis_in, vec_in, vec_out);
+    } else {
+      XDIAG_THROW(fmt::format("Unknown bond of type \"{}\"", type));
+    }
+  }
+
+}catch (Error const &e) {
+  XDIAG_RETHROW(e);
 }
 
 template void apply(OpSum const &, BasisSz<uint32_t> const &,

xdiag/basis/spinhalf_distributed/apply_ising.hpp

@@ -0,0 +1,102 @@
+#pragma once
+
+#include <algorithm>
+#include <tuple>
+
+#include <xdiag/extern/armadillo/armadillo>
+#include <xdiag/operators/opsum.hpp>
+
+namespace xdiag::basis::spinhalf_distributed {
+
+template <class basis_t, typename coeff_t>
+void apply_ising(Op const &op, basis_t const &basis,
+                 arma::Col<coeff_t> const &vec_in,
+                 arma::Col<coeff_t> &vec_out) try {
+  using bit_t = typename basis_t::bit_t;
+
+  assert(basis.size() == vec_in.size());
+  assert(basis.size() == vec_out.size());
+  assert(op.type() == "ISING");
+  assert(op.size() == 2);
+  int64_t ss1 = op[0];
+  int64_t ss2 = op[1];
+  assert((ss1 >= 0) && (ss2 >= 0));
+
+  if (ss1 == ss2) {
+    XDIAG_THROW("ISING Op with both sites equal not implemented yet");
+  }
+
+  assert(op.coupling_is<coeff_t>());
+  coeff_t J = op.coupling().as<coeff_t>();
+  coeff_t val_same = J / 4.;
+  coeff_t val_diff = -J / 4.;
+
+  int64_t s1 = std::min(ss1, ss2);
+  int64_t s2 = std::max(ss1, ss2);
+
+  int64_t n_postfix_bits = basis.n_postfix_bits();
+  bit_t mask = ((bit_t)1 << s1) | ((bit_t)1 << s2);
+  bit_t s1mask = (bit_t)1 << s1;
+  bit_t s2mask = (bit_t)1 << (s2 - n_postfix_bits);
+
+  int64_t idx = 0;
+  for (bit_t prefix : basis.prefixes()) {
+    bit_t prefix_shifted = (prefix << n_postfix_bits);
+    auto const &postfixes = basis.postfix_states(prefix);
+
+    // Both sites are on prefixes
+    if ((s1 >= n_postfix_bits) && (s2 >= n_postfix_bits)) {
+      coeff_t val =
+          (bits::popcnt(prefix_shifted & mask) & 1) ? val_diff : val_same;
+      int64_t end = idx + postfixes.size();
+      for (; idx < end; ++idx) {
+        vec_out(idx) += val * vec_in(idx);
+      }
+    }
+
+    // Both sites are on postfixes
+    else if ((s1 < n_postfix_bits) && (s2 < n_postfix_bits)) {
+
+      for (auto postfix : postfixes) {
+        if (bits::popcnt(postfix & mask) & 1) {
+          vec_out(idx) += val_diff * vec_in(idx);
+        } else {
+          vec_out(idx) += val_same * vec_in(idx);
+        }
+        ++idx;
+      }
+    }
+
+    // s2 is prefix s1 is postfix
+    else {
+
+      // s2 is up
+      if (prefix & s2mask) {
+        for (auto postfix : postfixes) {
+          if (postfix & s1mask) {
+            vec_out(idx) += val_same * vec_in(idx);
+          } else {
+            vec_out(idx) += val_diff * vec_in(idx);
+          }
+          ++idx;
+        }
+      }
+
+      // s2 is dn
+      else {
+        for (auto postfix : postfixes) {
+          if (postfix & s1mask) {
+            vec_out(idx) += val_diff * vec_in(idx);
+          } else {
+            vec_out(idx) += val_same * vec_in(idx);
+          }
+          ++idx;
+        }
+      }
+    }
+  }
+} catch (Error const &e) {
+  XDIAG_RETHROW(e);
+}
+
+} // namespace xdiag::basis::spinhalf_distributed

xdiag/basis/spinhalf_distributed/basis_spinhalf_distributed.hpp

@@ -5,6 +5,7 @@
 #include <xdiag/basis/spinhalf_distributed/basis_sz.hpp>
 
 namespace xdiag::basis {
+
 // clang-format off
 using BasisSpinhalfDistributed =
   std::variant<basis::spinhalf_distributed::BasisSz<uint32_t>,

xdiag/blocks/electron/terms/apply_exchange.hpp

@@ -34,10 +34,9 @@ void electron_do_down_flips(bit_t ups, int64_t idx_ups, bit_t flipmask,
 template <typename bit_t, typename coeff_t, bool symmetric, class Basis,
           class Filler>
 void apply_exchange(Op const &op, Basis &&basis, Filler &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined() && (op.type() == "EXCHANGE"));
+  assert(op.type() == "EXCHANGE");
   assert(op.size() == 2);
-  assert(op.sites_disjoint());
+  assert(sites_disjoint(op));
 
   Coupling cpl = op.coupling();
   assert(cpl.isexplicit() && !cpl.ismatrix());

xdiag/blocks/electron/terms/apply_hopping.hpp

@@ -12,10 +12,8 @@ namespace xdiag::electron {
 template <typename bit_t, typename coeff_t, bool symmetric, class Basis,
           class Fill>
 void apply_hopping(Op const &op, Basis &&basis, Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined());
   assert(op.size() == 2);
-  assert(op.sites_disjoint());
+  assert(sites_disjoint(op));
 
   std::string type = op.type();
   assert((type == "HOPUP") || (type == "HOPDN"));

xdiag/blocks/electron/terms/apply_ising.hpp

@@ -8,10 +8,9 @@ namespace xdiag::electron {
 template <typename bit_t, typename coeff_t, bool symmetric, class Basis,
           class Fill>
 void apply_ising(Op const &op, Basis &&basis, Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined() && (op.type() == "ISING"));
+  assert(op.type() == "ISING");
   assert(op.size() == 2);
-  assert(op.sites_disjoint());
+  assert(sites_disjoint(op));
 
   Coupling cpl = op.coupling();
   assert(cpl.isexplicit() && !cpl.ismatrix());

xdiag/blocks/electron/terms/apply_number.hpp

@@ -5,8 +5,6 @@ namespace xdiag::electron {
 template <typename bit_t, typename coeff_t, bool symmetric, class Basis,
           class Fill>
 void apply_number(Op const &op, Basis &&basis, Fill fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined());
   assert(op.size() == 1);
 
   std::string type = op.type();

xdiag/blocks/electron/terms/apply_raise_lower.hpp

@@ -12,8 +12,6 @@ template <typename bit_t, typename coeff_t, bool symmetric, class BasisIn,
           class BasisOut, class Fill>
 void apply_raise_lower(Op const &op, BasisIn &&basis_in,
                        BasisOut &&basis_out, Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined());
   assert(op.size() == 1);
 
   std::string type = op.type();

xdiag/blocks/spinhalf/terms/apply_exchange.hpp

@@ -16,14 +16,13 @@ template <typename bit_t, typename coeff_t, bool symmetric, class BasisIn,
 void apply_exchange(Op const &op, BasisIn &&basis_in, BasisOut &&basis_out,
                     Fill &&fill) {
 
-  assert(op.coupling_defined());
-  assert(op.type_defined() && (op.type() == "EXCHANGE"));
+  assert(op.type() == "EXCHANGE");
   assert(op.size() == 2);
-  assert(op.sites_disjoint());
+  assert(sites_disjoint(op));
 
   Coupling cpl = op.coupling();
   assert(cpl.isexplicit() && !cpl.ismatrix());
-  
+
   coeff_t J = cpl.as<coeff_t>();
   int64_t s1 = op[0];
   int64_t s2 = op[1];

xdiag/blocks/spinhalf/terms/apply_ising.hpp

@@ -12,10 +12,9 @@ template <typename bit_t, typename coeff_t, bool symmetric, class BasisIn,
           class BasisOut, class Fill>
 void apply_ising(Op const &op, BasisIn &&basis_in, BasisOut &&basis_out,
                  Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined() && (op.type() == "ISING"));
+  assert(op.type() == "ISING");
   assert(op.size() == 2);
-  assert(op.sites_disjoint());
+  assert(sites_disjoint(op));
 
   Coupling cpl = op.coupling();
   assert(cpl.isexplicit() && !cpl.ismatrix());

xdiag/blocks/spinhalf/terms/apply_non_branching.hpp

@@ -14,9 +14,8 @@ template <typename bit_t, typename coeff_t, bool symmetric, class BasisIn,
 void apply_non_branching(Op const &op, BasisIn &&basis_in,
                          BasisOut &&basis_out, Fill &&fill) {
 
-  assert(op.matrix_defined());
+  assert(op.ismatrix());
   assert(operators::is_non_branching_op(op));
-  assert(op.coupling_defined());
 
   auto op_nb = operators::NonBranchingOp<bit_t, coeff_t>(op);
 

xdiag/blocks/spinhalf/terms/apply_scalar_chirality.hpp

@@ -15,10 +15,9 @@ void apply_scalar_chirality(Op const &op, BasisIn &&basis_in,
                             BasisOut &&basis_out, Fill &&fill) try {
   using bits::gbit;
 
-  assert(op.coupling_defined());
-  assert(op.type_defined() && (op.type() == "SCALARCHIRALITY"));
+  assert(op.type() == "SCALARCHIRALITY");
   assert(op.size() == 3);
-  assert(op.sites_disjoint());
+  assert(sites_disjoint(op));
   assert(iscomplex<coeff_t>());
 
   Coupling cpl = op.coupling();

xdiag/blocks/spinhalf/terms/apply_spsm.hpp

@@ -16,8 +16,6 @@ template <typename bit_t, typename coeff_t, bool symmetric, class BasisIn,
           class BasisOut, class Fill>
 void apply_spsm(Op const &op, BasisIn &&basis_in, BasisOut &&basis_out,
                 Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined());
   assert((op.type() == "S+") || (op.type() == "S-"));
   assert(op.size() == 1);
 

xdiag/blocks/spinhalf/terms/apply_sz.hpp

@@ -12,8 +12,7 @@ template <typename bit_t, typename coeff_t, bool symmetric, class BasisIn,
           class BasisOut, class Fill>
 void apply_sz(Op const &op, BasisIn &&basis_in, BasisOut &&basis_out,
               Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined() && (op.type() == "SZ"));
+  assert(op.type() == "SZ");
   assert(op.size() == 1);
 
   Coupling cpl = op.coupling();

xdiag/blocks/tj/terms/apply_exchange.hpp

@@ -10,10 +10,8 @@ namespace xdiag::tj {
 template <typename bit_t, typename coeff_t, bool symmetric, class Basis,
           class Filler>
 void apply_exchange(Op const &op, Basis &&basis, Filler &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined());
   assert(op.size() == 2);
-  assert(op.sites_disjoint());
+  assert(sites_disjoint(op));
   std::string type = op.type();
   assert(type == "EXCHANGE");
 

xdiag/blocks/tj/terms/apply_hopping.hpp

@@ -11,10 +11,8 @@ namespace xdiag::tj {
 template <typename bit_t, typename coeff_t, bool symmetric, class Basis,
           class Fill>
 void apply_hopping(Op const &op, Basis &&basis, Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined());
   assert(op.size() == 2);
-  assert(op.sites_disjoint());
+  assert(sites_disjoint(op));
 
   std::string type = op.type();
   assert((type == "HOPUP") || (type == "HOPDN"));

xdiag/blocks/tj/terms/apply_ising.hpp

@@ -1,26 +1,25 @@
 #pragma once
 
+#include <xdiag/blocks/tj/terms/generic_term_diag.hpp>
 #include <xdiag/common.hpp>
 #include <xdiag/operators/op.hpp>
-#include <xdiag/blocks/tj/terms/generic_term_diag.hpp>
 
 namespace xdiag::tj {
 
 template <typename bit_t, typename coeff_t, bool symmetric, class Basis,
           class Fill>
 void apply_ising(Op const &op, Basis &&basis, Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined());
+
   assert(op.size() == 2);
-  assert(op.sites_disjoint());
+  assert(sites_disjoint(op));
 
   std::string type = op.type();
   assert((type == "ISING") || (type == "TJISING"));
 
   Coupling cpl = op.coupling();
   assert(cpl.isexplicit() && !cpl.ismatrix());
   coeff_t J = cpl.as<coeff_t>();
-  
+
   int64_t s1 = op[0];
   int64_t s2 = op[1];
   bit_t s1_mask = (bit_t)1 << s1;

xdiag/blocks/tj/terms/apply_number.hpp

@@ -9,8 +9,6 @@ namespace xdiag::tj {
 template <typename bit_t, typename coeff_t, bool symmetric, class Basis,
           class Fill>
 void apply_number(Op const &op, Basis &&basis, Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined());
   assert(op.size() == 1);
 
   std::string type = op.type();

xdiag/blocks/tj/terms/apply_raise_lower.hpp

@@ -12,8 +12,6 @@ template <typename bit_t, typename coeff_t, bool symmetric, class BasisIn,
           class BasisOut, class Fill>
 void apply_raise_lower(Op const &op, BasisIn &&basis_in,
                        BasisOut &&basis_out, Fill &&fill) {
-  assert(op.coupling_defined());
-  assert(op.type_defined());
   assert(op.size() == 1);
 
   std::string type = op.type();