Segmentation fault when performing operations on operators
As pointed out by @femtobit in #203, when we perform arithmetic operations with operators in python, creating new objects, we do not tell the reference counter in python that there are new references to the `hilbert` object. This can lead to segmentation faults (See example below):

```python
import netket as nk

g = nk.graph.Hypercube(8, 1)
hi = nk.hilbert.Spin(g, 1/2)

lo = nk.operator.LocalOperator(hi, [[1,0],[0,1]], [0])
lo = lo*lo

# Assign hi to something else so it will no longer keep alive the original Hilbert object
hi = nk.hilbert.Qubit(g)

# This will try to access the now deleted Hilbert object
lo = lo * lo
```

One possible solution might be to use `std::shared_ptr` in place of a reference. Otherwise we must find a way to explicitly update python/pybind reference counter.
diff --git a/Sources/Hilbert/py_hilbert.cc b/Sources/Hilbert/py_hilbert.cc
index e8936a08..1862f9cc 100644
--- a/Sources/Hilbert/py_hilbert.cc
+++ b/Sources/Hilbert/py_hilbert.cc
@@ -31,7 +31,7 @@ namespace netket {
 
 namespace {
 void AddBosons(py::module subm) {
-  py::class_<Boson, AbstractHilbert>(
+  py::class_<Boson, AbstractHilbert, std::shared_ptr<Boson>>(
       subm, "Boson", R"EOF(Hilbert space composed of bosonic states.)EOF")
       .def(py::init<const AbstractGraph &, int>(), py::keep_alive<1, 2>(),
            py::arg("graph"), py::arg("n_max"), R"EOF(
@@ -80,8 +80,8 @@ void AddBosons(py::module subm) {
 }
 
 void AddCustomHilbert(py::module subm) {
-  py::class_<CustomHilbert, AbstractHilbert>(subm, "CustomHilbert",
-                                             R"EOF(A custom hilbert space.)EOF")
+  py::class_<CustomHilbert, AbstractHilbert, std::shared_ptr<CustomHilbert>>(
+      subm, "CustomHilbert", R"EOF(A custom hilbert space.)EOF")
       .def(py::init<const AbstractGraph &, std::vector<double>>(),
            py::keep_alive<1, 2>(), py::arg("graph"), py::arg("local_states"),
            R"EOF(
@@ -108,7 +108,7 @@ void AddCustomHilbert(py::module subm) {
 }
 
 void AddSpins(py::module subm) {
-  py::class_<Spin, AbstractHilbert>(
+  py::class_<Spin, AbstractHilbert, std::shared_ptr<Spin>>(
       subm, "Spin", R"EOF(Hilbert space composed of spin states.)EOF")
       .def(py::init<const AbstractGraph &, double>(), py::keep_alive<1, 2>(),
            py::arg("graph"), py::arg("s"), R"EOF(
@@ -161,7 +161,8 @@ void AddHilbertModule(py::module m) {
   auto subm = m.def_submodule("hilbert");
 
   auto hilbert_class =
-      py::class_<AbstractHilbert>(subm, "Hilbert")
+      py::class_<AbstractHilbert, std::shared_ptr<AbstractHilbert>>(subm,
+                                                                    "Hilbert")
           .def_property_readonly(
               "is_discrete", &AbstractHilbert::IsDiscrete,
               R"EOF(bool: Whether the hilbert space is discrete.)EOF")
diff --git a/Sources/Machine/DensityMatrices/abstract_density_matrix.hpp b/Sources/Machine/DensityMatrices/abstract_density_matrix.hpp
index 0e6a9b36..7c6b200a 100644
--- a/Sources/Machine/DensityMatrices/abstract_density_matrix.hpp
+++ b/Sources/Machine/DensityMatrices/abstract_density_matrix.hpp
@@ -29,12 +29,9 @@ namespace netket {
  */
 class AbstractDensityMatrix : public AbstractMachine {
   // The physical hilbert space over which this operator acts
-  const AbstractHilbert &hilbert_physical_;
 
-  const std::unique_ptr<AbstractGraph> graph_doubled_;
-
-  // The doubled hilbert space where this operator is defined
-  const std::unique_ptr<AbstractHilbert> hilbert_doubled_;
+  std::shared_ptr<const AbstractHilbert> hilbert_physical_;
+  std::unique_ptr<const AbstractGraph> graph_doubled_;
 
   using Edge = AbstractGraph::Edge;
 
@@ -78,25 +75,34 @@ class AbstractDensityMatrix : public AbstractMachine {
         CustomGraph(d_edges, d_eclist, d_automorphisms));
   }
 
+ protected:
+  AbstractDensityMatrix(std::shared_ptr<const AbstractHilbert> physical_hilbert,
+                        std::unique_ptr<const AbstractGraph> doubled_graph)
+      : AbstractMachine(std::make_shared<CustomHilbert>(
+            *doubled_graph, physical_hilbert->LocalStates())),
+        hilbert_physical_(physical_hilbert),
+        graph_doubled_(std::move(doubled_graph)) {}
+
  public:
-  explicit AbstractDensityMatrix(const AbstractHilbert &hilbert)
-      : hilbert_physical_(hilbert),
-        graph_doubled_(DoubledGraph(hilbert.GetGraph())),
-        hilbert_doubled_(make_unique<CustomHilbert>(*graph_doubled_,
-                                                    hilbert.LocalStates())){};
-
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return *hilbert_doubled_;
-  }
+  explicit AbstractDensityMatrix(std::shared_ptr<const AbstractHilbert> hilbert)
+      : AbstractDensityMatrix(hilbert, DoubledGraph(hilbert->GetGraph())){};
 
   /**
-   * Member function returning the Physical hilbert space over which
-   * this density matrix acts
+   * Member function returning the physical hilbert space over which
+   * this density matrix acts as a shared pointer.
    * @return The physical hilbert space
    */
-  const AbstractHilbert &GetHilbertPhysical() const noexcept {
+  std::shared_ptr<const AbstractHilbert> GetHilbertPhysicalShared() const {
     return hilbert_physical_;
   }
+
+  /* Member function returning a reference to the physical hilbert space
+   * on which this density matrix acts.
+   * @return The physical hilbert space
+   */
+  const AbstractHilbert &GetHilbertPhysical() const noexcept {
+    return *hilbert_physical_;
+  }
 };
 }  // namespace netket
 
diff --git a/Sources/Machine/DensityMatrices/diagonal_density_matrix.hpp b/Sources/Machine/DensityMatrices/diagonal_density_matrix.hpp
index e17b315a..33fb2db1 100644
--- a/Sources/Machine/DensityMatrices/diagonal_density_matrix.hpp
+++ b/Sources/Machine/DensityMatrices/diagonal_density_matrix.hpp
@@ -35,11 +35,7 @@ class DiagonalDensityMatrix : public AbstractMachine {
 
  public:
   explicit DiagonalDensityMatrix(AbstractDensityMatrix &dm)
-      : density_matrix_(dm){};
-
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return density_matrix_.GetHilbertPhysical();
-  }
+      : AbstractMachine(dm.GetHilbertPhysicalShared()), density_matrix_(dm){};
 
   const AbstractDensityMatrix &GetFullDensityMatrix() const noexcept {
     return density_matrix_;
diff --git a/Sources/Machine/DensityMatrices/ndm_spin_phase.hpp b/Sources/Machine/DensityMatrices/ndm_spin_phase.hpp
index bc0cff71..08b34678 100644
--- a/Sources/Machine/DensityMatrices/ndm_spin_phase.hpp
+++ b/Sources/Machine/DensityMatrices/ndm_spin_phase.hpp
@@ -95,18 +95,18 @@ class NdmSpinPhase : public AbstractDensityMatrix {
   const Complex I_;
 
  public:
-  explicit NdmSpinPhase(const AbstractHilbert &hilbert, int nhidden = 0,
-                        int nancilla = 0, int alpha = 0, int beta = 0,
-                        bool useb = true, bool useh = true, bool used = true)
+  explicit NdmSpinPhase(std::shared_ptr<const AbstractHilbert> hilbert,
+                        int nhidden = 0, int nancilla = 0, int alpha = 0,
+                        int beta = 0, bool useb = true, bool useh = true,
+                        bool used = true)
       : AbstractDensityMatrix(hilbert),
-        nv_(hilbert.Size()),
+        nv_(hilbert->Size()),
         useb_(useb),
         useh_(useh),
         used_(used),
         I_(0, 1) {
     nh_ = std::max(nhidden, alpha * nv_);
     na_ = std::max(nancilla, beta * nv_);
-
     Init();
   }
 
diff --git a/Sources/Machine/DensityMatrices/py_ndm_spin_phase.hpp b/Sources/Machine/DensityMatrices/py_ndm_spin_phase.hpp
index 58764d27..362646ba 100644
--- a/Sources/Machine/DensityMatrices/py_ndm_spin_phase.hpp
+++ b/Sources/Machine/DensityMatrices/py_ndm_spin_phase.hpp
@@ -36,9 +36,9 @@ void AddNdmSpinPhase(py::module &subm) {
           In this case, two NDMs are taken to parameterize, respectively, phase
           and amplitude of the density matrix.
           This type of NDM has spin 1/2 hidden and ancilla units.)EOF")
-      .def(py::init<const AbstractHilbert &, int, int, int, int, bool, bool,
-                    bool>(),
-           py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("n_hidden") = 0,
+      .def(py::init<std::shared_ptr<const AbstractHilbert>, int, int, int, int,
+                    bool, bool, bool>(),
+           py::arg("hilbert"), py::arg("n_hidden") = 0,
            py::arg("n_ancilla") = 0, py::arg("alpha") = 0, py::arg("beta") = 0,
            py::arg("use_visible_bias") = true,
            py::arg("use_hidden_bias") = true,
diff --git a/Sources/Machine/abstract_machine.cc b/Sources/Machine/abstract_machine.cc
index d189f7b7..854905cb 100644
--- a/Sources/Machine/abstract_machine.cc
+++ b/Sources/Machine/abstract_machine.cc
@@ -24,7 +24,7 @@ AbstractMachine::VectorType AbstractMachine::DerLogChanged(
     VisibleConstType v, const std::vector<int> &tochange,
     const std::vector<double> &newconf) {
   VisibleType vp(v);
-  GetHilbert().UpdateConf(vp, tochange, newconf);
+  hilbert_->UpdateConf(vp, tochange, newconf);
   return DerLog(vp);
 }
 
diff --git a/Sources/Machine/abstract_machine.hpp b/Sources/Machine/abstract_machine.hpp
index 9fc228f0..503ba7e4 100644
--- a/Sources/Machine/abstract_machine.hpp
+++ b/Sources/Machine/abstract_machine.hpp
@@ -17,6 +17,7 @@
 
 #include <complex>
 #include <iosfwd>
+#include <memory>
 #include <vector>
 
 #include <Eigen/Core>
@@ -214,9 +215,19 @@ class AbstractMachine {
   void Save(const std::string &filename) const;
   void Save(std::ofstream &stream) const;
 
-  virtual const AbstractHilbert &GetHilbert() const noexcept = 0;
+  const AbstractHilbert &GetHilbert() const { return *hilbert_; };
+  std::shared_ptr<const AbstractHilbert> GetHilbertShared() const {
+    return hilbert_;
+  };
 
-  virtual ~AbstractMachine() {}
+  virtual ~AbstractMachine() = default;
+
+ protected:
+  AbstractMachine(std::shared_ptr<const AbstractHilbert> hilbert)
+      : hilbert_(std::move(hilbert)) {}
+
+ private:
+  std::shared_ptr<const AbstractHilbert> hilbert_;
 };
 }  // namespace netket
 
diff --git a/Sources/Machine/ffnn.hpp b/Sources/Machine/ffnn.hpp
index 112584f1..d536a534 100644
--- a/Sources/Machine/ffnn.hpp
+++ b/Sources/Machine/ffnn.hpp
@@ -31,8 +31,6 @@
 namespace netket {
 
 class FFNN : public AbstractMachine {
-  const AbstractHilbert &hilbert_;
-
   std::vector<AbstractLayer *> layers_;  // Pointers to hidden layers
 
   std::vector<int> layersizes_;
@@ -49,9 +47,11 @@ class FFNN : public AbstractMachine {
   std::unique_ptr<SumOutput> sum_output_layer_;
 
  public:
-  explicit FFNN(const AbstractHilbert &hilbert,
+  explicit FFNN(std::shared_ptr<const AbstractHilbert> hilbert,
                 std::vector<AbstractLayer *> layers)
-      : hilbert_(hilbert), layers_(std::move(layers)), nv_(hilbert.Size()) {
+      : AbstractMachine(hilbert),
+        layers_(std::move(layers)),
+        nv_(hilbert->Size()) {
     Init();
   }
 
@@ -295,9 +295,6 @@ class FFNN : public AbstractMachine {
     }
   }
 
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
 };  // namespace netket
 
 }  // namespace netket
diff --git a/Sources/Machine/jastrow.cc b/Sources/Machine/jastrow.cc
index 37760b30..3d2f3a13 100644
--- a/Sources/Machine/jastrow.cc
+++ b/Sources/Machine/jastrow.cc
@@ -21,8 +21,8 @@
 
 namespace netket {
 
-Jastrow::Jastrow(const AbstractHilbert &hilbert)
-    : hilbert_(hilbert), nv_(hilbert.Size()) {
+Jastrow::Jastrow(std::shared_ptr<const AbstractHilbert> hilbert)
+    : AbstractMachine(hilbert), nv_(hilbert->Size()) {
   Init();
 }
 
@@ -182,8 +182,6 @@ Jastrow::VectorType Jastrow::DerLog(VisibleConstType v) {
   return der;
 }
 
-const AbstractHilbert &Jastrow::GetHilbert() const noexcept { return hilbert_; }
-
 void Jastrow::to_json(json &j) const {
   j["Name"] = "Jastrow";
   j["Nvisible"] = nv_;
@@ -198,7 +196,7 @@ void Jastrow::from_json(const json &pars) {
   if (FieldExists(pars, "Nvisible")) {
     nv_ = pars["Nvisible"];
   }
-  if (nv_ != hilbert_.Size()) {
+  if (nv_ != GetHilbert().Size()) {
     throw InvalidInputError(
         "Number of visible units is incompatible with given "
         "Hilbert space");
diff --git a/Sources/Machine/jastrow.hpp b/Sources/Machine/jastrow.hpp
index 771a81de..41db4825 100644
--- a/Sources/Machine/jastrow.hpp
+++ b/Sources/Machine/jastrow.hpp
@@ -25,8 +25,6 @@ namespace netket {
  *
  */
 class Jastrow : public AbstractMachine {
-  const AbstractHilbert &hilbert_;
-
   // number of visible units
   int nv_;
 
@@ -41,7 +39,7 @@ class Jastrow : public AbstractMachine {
   VectorType thetasnew_;
 
  public:
-  explicit Jastrow(const AbstractHilbert &hilbert);
+  explicit Jastrow(std::shared_ptr<const AbstractHilbert> hilbert);
 
   inline void Init();
 
@@ -67,7 +65,6 @@ class Jastrow : public AbstractMachine {
                      const LookupType &lt) override;
 
   VectorType DerLog(VisibleConstType v) override;
-  const AbstractHilbert &GetHilbert() const noexcept override;
 
   void to_json(json &j) const override;
   void from_json(const json &pars) override;
diff --git a/Sources/Machine/jastrow_symm.cc b/Sources/Machine/jastrow_symm.cc
index fd6b87c4..6f0fe830 100644
--- a/Sources/Machine/jastrow_symm.cc
+++ b/Sources/Machine/jastrow_symm.cc
@@ -22,10 +22,11 @@
 
 namespace netket {
 
-JastrowSymm::JastrowSymm(const AbstractHilbert &hilbert)
-    : hilbert_(hilbert), graph_(hilbert.GetGraph()), nv_(hilbert.Size()) {
+JastrowSymm::JastrowSymm(std::shared_ptr<const AbstractHilbert> hilbert)
+    : AbstractMachine(hilbert),
+      graph_(hilbert->GetGraph()),
+      nv_(hilbert->Size()) {
   Init(graph_);
-
   SetBareParameters();
 }
 
@@ -271,10 +272,6 @@ Complex JastrowSymm::LogValDiff(VisibleConstType v,
   return logvaldiff;
 }
 
-const AbstractHilbert &JastrowSymm::GetHilbert() const noexcept {
-  return hilbert_;
-}
-
 void JastrowSymm::to_json(json &j) const {
   j["Name"] = "JastrowSymm";
   j["Nvisible"] = nv_;
@@ -290,7 +287,7 @@ void JastrowSymm::from_json(const json &pars) {
   if (FieldExists(pars, "Nvisible")) {
     nv_ = pars["Nvisible"];
   }
-  if (nv_ != hilbert_.Size()) {
+  if (nv_ != GetHilbert().Size()) {
     throw InvalidInputError(
         "Number of visible units is incompatible with given "
         "Hilbert space");
diff --git a/Sources/Machine/jastrow_symm.hpp b/Sources/Machine/jastrow_symm.hpp
index e4911e9e..00bfb512 100644
--- a/Sources/Machine/jastrow_symm.hpp
+++ b/Sources/Machine/jastrow_symm.hpp
@@ -24,7 +24,6 @@ namespace netket {
 
 // Jastrow with permutation symmetries
 class JastrowSymm : public AbstractMachine {
-  const AbstractHilbert &hilbert_;
   const AbstractGraph &graph_;
 
   std::vector<std::vector<int>> permtable_;
@@ -52,13 +51,11 @@ class JastrowSymm : public AbstractMachine {
   Eigen::MatrixXi Wtemp_;
 
  public:
-  explicit JastrowSymm(const AbstractHilbert &hilbert);
+  explicit JastrowSymm(std::shared_ptr<const AbstractHilbert> hilbert);
 
   int Nvisible() const override;
   int Npar() const override;
 
-  const AbstractHilbert &GetHilbert() const noexcept override;
-
   void InitRandomPars(int seed, double sigma) override;
   void InitLookup(VisibleConstType v, LookupType &lt) override;
   void UpdateLookup(VisibleConstType v, const std::vector<int> &tochange,
diff --git a/Sources/Machine/mps_periodic.cc b/Sources/Machine/mps_periodic.cc
index 792dd324..b22e20b3 100644
--- a/Sources/Machine/mps_periodic.cc
+++ b/Sources/Machine/mps_periodic.cc
@@ -23,18 +23,17 @@
 
 namespace netket {
 
-MPSPeriodic::MPSPeriodic(const AbstractHilbert &hilbert, int bond_dim,
-                         bool diag, int symperiod)
-    : hilbert_{hilbert},
-      N_{hilbert.Size()},
-      d_{hilbert.LocalSize()},
+MPSPeriodic::MPSPeriodic(std::shared_ptr<const AbstractHilbert> hilbert,
+                         int bond_dim, bool diag, int symperiod)
+    : AbstractMachine(hilbert),
+      N_{hilbert->Size()},
+      d_{hilbert->LocalSize()},
       D_{bond_dim},
       symperiod_{symperiod},
       is_diag_{diag} {
   if (symperiod_ == -1) {
     symperiod_ = N_;
   }
-
   Init();
 }
 
@@ -114,7 +113,7 @@ void MPSPeriodic::Init() {
                   << std::endl;
   }
   // Initialize map from Hilbert space states to MPS indices
-  auto localstates = hilbert_.LocalStates();
+  auto localstates = GetHilbert().LocalStates();
   for (int i = 0; i < d_; i++) {
     confindex_[localstates[i]] = i;
   }
@@ -196,10 +195,6 @@ void MPSPeriodic::InitTree() {
 int MPSPeriodic::Npar() const { return npar_; }
 int MPSPeriodic::Nvisible() const { return N_; }
 
-const AbstractHilbert &MPSPeriodic::GetHilbert() const noexcept {
-  return hilbert_;
-}
-
 MPSPeriodic::VectorType MPSPeriodic::GetParameters() {
   int k = 0;
   VectorType pars(npar_);
@@ -533,7 +528,7 @@ void MPSPeriodic::from_json(const json &pars) {
   if (FieldExists(pars, "Length")) {
     N_ = pars["Length"];
   }
-  if (N_ != hilbert_.Size()) {
+  if (N_ != GetHilbert().Size()) {
     throw InvalidInputError(
         "Number of spins is incompatible with given Hilbert space");
   }
@@ -541,7 +536,7 @@ void MPSPeriodic::from_json(const json &pars) {
   if (FieldExists(pars, "PhysDim")) {
     d_ = pars["PhysDim"];
   }
-  if (d_ != hilbert_.LocalSize()) {
+  if (d_ != GetHilbert().LocalSize()) {
     throw InvalidInputError(
         "Number of spins is incompatible with given Hilbert space");
   }
diff --git a/Sources/Machine/mps_periodic.hpp b/Sources/Machine/mps_periodic.hpp
index ad546f7e..6258e2bb 100644
--- a/Sources/Machine/mps_periodic.hpp
+++ b/Sources/Machine/mps_periodic.hpp
@@ -24,8 +24,6 @@
 namespace netket {
 
 class MPSPeriodic : public AbstractMachine {
-  const AbstractHilbert &hilbert_;
-
   // Number of sites
   int N_;
   // Physical dimension
@@ -61,14 +59,12 @@ class MPSPeriodic : public AbstractMachine {
   MatrixType identity_mat_;
 
  public:
-  MPSPeriodic(const AbstractHilbert &hilbert, int bond_dim, bool diag,
-              int symperiod = -1);
+  MPSPeriodic(std::shared_ptr<const AbstractHilbert> hilbert, int bond_dim,
+              bool diag, int symperiod = -1);
 
   int Npar() const override;
   int Nvisible() const override;
 
-  const AbstractHilbert &GetHilbert() const noexcept override;
-
   void InitRandomPars(int seed, double sigma) override;
   VectorType GetParameters() override;
   void SetParameters(VectorConstRefType pars) override;
diff --git a/Sources/Machine/py_machine.cc b/Sources/Machine/py_machine.cc
index da64292d..b595587f 100644
--- a/Sources/Machine/py_machine.cc
+++ b/Sources/Machine/py_machine.cc
@@ -50,9 +50,10 @@ void AddRbmSpin(py::module subm) {
              \left(\sum_i^N W_{ij} s_i + b_j \right) $$
 
           for arbitrary local quantum numbers $$ s_i $$.)EOF")
-      .def(py::init<const AbstractHilbert &, int, int, bool, bool>(),
-           py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("n_hidden") = 0,
-           py::arg("alpha") = 0, py::arg("use_visible_bias") = true,
+      .def(py::init<std::shared_ptr<const AbstractHilbert>, int, int, bool,
+                    bool>(),
+           py::arg("hilbert"), py::arg("n_hidden") = 0, py::arg("alpha") = 0,
+           py::arg("use_visible_bias") = true,
            py::arg("use_hidden_bias") = true,
            R"EOF(
                    Constructs a new ``RbmSpin`` machine:
@@ -98,8 +99,8 @@ void AddRbmSpinSymm(py::module subm) {
              for arbitrary local quantum numbers $$ s_i $$. However, the weights
              ($$ W_{ij} $$) and biases ($$ a_i $$, $$ b_i $$) respects the
              specified symmetries of the lattice.)EOF")
-      .def(py::init<const AbstractHilbert &, int, bool, bool>(),
-           py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("alpha") = 0,
+      .def(py::init<std::shared_ptr<const AbstractHilbert>, int, bool, bool>(),
+           py::arg("hilbert"), py::arg("alpha") = 0,
            py::arg("use_visible_bias") = true,
            py::arg("use_hidden_bias") = true,
            R"EOF(
@@ -137,9 +138,10 @@ void AddRbmMultival(py::module subm) {
   py::class_<RbmMultival, AbstractMachine>(subm, "RbmMultiVal", R"EOF(
              A fully connected Restricted Boltzmann Machine for handling larger
              local Hilbert spaces.)EOF")
-      .def(py::init<const AbstractHilbert &, int, int, bool, bool>(),
-           py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("n_hidden") = 0,
-           py::arg("alpha") = 0, py::arg("use_visible_bias") = true,
+      .def(py::init<std::shared_ptr<const AbstractHilbert>, int, int, bool,
+                    bool>(),
+           py::arg("hilbert"), py::arg("n_hidden") = 0, py::arg("alpha") = 0,
+           py::arg("use_visible_bias") = true,
            py::arg("use_hidden_bias") = true);
 }
 
@@ -154,9 +156,10 @@ void AddRbmSpinPhase(py::module subm) {
              \left(\sum_i^N W_{ij} s_i + b_j \right) $$
 
           for arbitrary local quantum numbers $$ s_i $$.)EOF")
-      .def(py::init<const AbstractHilbert &, int, int, bool, bool>(),
-           py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("n_hidden") = 0,
-           py::arg("alpha") = 0, py::arg("use_visible_bias") = true,
+      .def(py::init<std::shared_ptr<const AbstractHilbert>, int, int, bool,
+                    bool>(),
+           py::arg("hilbert"), py::arg("n_hidden") = 0, py::arg("alpha") = 0,
+           py::arg("use_visible_bias") = true,
            py::arg("use_hidden_bias") = true,
            R"EOF(
                    Constructs a new ``RbmSpinPhase`` machine:
@@ -199,9 +202,10 @@ void AddRbmSpinReal(py::module subm) {
              \left(\sum_i^N W_{ij} s_i + b_j \right) $$
 
           for arbitrary local quantum numbers $$ s_i $$.)EOF")
-      .def(py::init<const AbstractHilbert &, int, int, bool, bool>(),
-           py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("n_hidden") = 0,
-           py::arg("alpha") = 0, py::arg("use_visible_bias") = true,
+      .def(py::init<std::shared_ptr<const AbstractHilbert>, int, int, bool,
+                    bool>(),
+           py::arg("hilbert"), py::arg("n_hidden") = 0, py::arg("alpha") = 0,
+           py::arg("use_visible_bias") = true,
            py::arg("use_hidden_bias") = true,
            R"EOF(
                    Constructs a new ``RbmSpinReal`` machine:
@@ -242,12 +246,12 @@ void AddFFNN(py::module subm) {
              class. Each layer implements a transformation such that the
              information is transformed sequentially as it moves from the input
              nodes through the hidden layers and to the output nodes.)EOF")
-      .def(py::init([](AbstractHilbert const &hi, py::tuple tuple) {
-             auto layers = py::cast<std::vector<AbstractLayer *>>(tuple);
-             return FFNN{hi, std::move(layers)};
-           }),
-           py::keep_alive<1, 2>(), py::keep_alive<1, 3>(), py::arg("hilbert"),
-           py::arg("layers"),
+      .def(py::init(
+               [](std::shared_ptr<const AbstractHilbert> hi, py::tuple tuple) {
+                 auto layers = py::cast<std::vector<AbstractLayer *>>(tuple);
+                 return FFNN{std::move(hi), std::move(layers)};
+               }),
+           py::keep_alive<1, 3>(), py::arg("hilbert"), py::arg("layers"),
            R"EOF(
               Constructs a new ``FFNN`` machine:
 
@@ -286,7 +290,7 @@ void AddJastrow(py::module subm) {
 
            where $$ W_{ij} $$ are the Jastrow parameters.
            )EOF")
-      .def(py::init<const AbstractHilbert &>(), py::keep_alive<1, 2>(),
+      .def(py::init<std::shared_ptr<const AbstractHilbert>>(),
            py::arg("hilbert"), R"EOF(
                  Constructs a new ``Jastrow`` machine:
 
@@ -320,7 +324,7 @@ void AddJastrowSymm(py::module subm) {
 
            where $$ W_{ij} $$ are the Jastrow parameters respects the
            specified symmetries of the lattice.)EOF")
-      .def(py::init<const AbstractHilbert &>(), py::keep_alive<1, 2>(),
+      .def(py::init<std::shared_ptr<const AbstractHilbert>>(),
            py::arg("hilbert"), R"EOF(
                  Constructs a new ``JastrowSymm`` machine:
 
@@ -347,9 +351,9 @@ void AddJastrowSymm(py::module subm) {
 
 void AddMpsPeriodic(py::module subm) {
   py::class_<MPSPeriodic, AbstractMachine>(subm, "MPSPeriodic")
-      .def(py::init<const AbstractHilbert &, int, bool, int>(),
-           py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("bond_dim"),
-           py::arg("diag") = false, py::arg("symperiod") = -1);
+      .def(py::init<std::shared_ptr<const AbstractHilbert>, int, bool, int>(),
+           py::arg("hilbert"), py::arg("bond_dim"), py::arg("diag") = false,
+           py::arg("symperiod") = -1);
 }
 
 void AddLayerModule(py::module m) {
diff --git a/Sources/Machine/rbm_multival.cc b/Sources/Machine/rbm_multival.cc
index 26e78cd2..107c7d1b 100644
--- a/Sources/Machine/rbm_multival.cc
+++ b/Sources/Machine/rbm_multival.cc
@@ -25,11 +25,11 @@
 
 namespace netket {
 
-RbmMultival::RbmMultival(const AbstractHilbert &hilbert, int nhidden, int alpha,
-                         bool usea, bool useb)
-    : hilbert_(hilbert),
-      nv_(hilbert.Size()),
-      ls_(hilbert.LocalSize()),
+RbmMultival::RbmMultival(std::shared_ptr<const AbstractHilbert> hilbert,
+                         int nhidden, int alpha, bool usea, bool useb)
+    : AbstractMachine(hilbert),
+      nv_(hilbert->Size()),
+      ls_(hilbert->LocalSize()),
       usea_(usea),
       useb_(useb) {
   nh_ = std::max(nhidden, alpha * nv_);
@@ -60,7 +60,7 @@ void RbmMultival::Init() {
     b_.setZero();
   }
 
-  auto localstates = hilbert_.LocalStates();
+  auto localstates = GetHilbert().LocalStates();
 
   localconfs_.resize(nv_ * ls_);
   for (int i = 0; i < nv_ * ls_; i += ls_) {
@@ -317,10 +317,6 @@ inline void RbmMultival::ComputeVtilde(VisibleConstType v,
 }
 #endif
 
-const AbstractHilbert &RbmMultival::GetHilbert() const noexcept {
-  return hilbert_;
-}
-
 void RbmMultival::to_json(json &j) const {
   j["Name"] = "RbmMultival";
   j["Nvisible"] = nv_;
@@ -343,7 +339,7 @@ void RbmMultival::from_json(const json &pars) {
     nv_ = pars["Nvisible"];
   }
 
-  if (nv_ != hilbert_.Size()) {
+  if (nv_ != GetHilbert().Size()) {
     throw InvalidInputError(
         "Loaded wave-function has incompatible Hilbert space");
   }
@@ -351,7 +347,7 @@ void RbmMultival::from_json(const json &pars) {
   if (FieldExists(pars, "LocalSize")) {
     ls_ = pars["LocalSize"];
   }
-  if (ls_ != hilbert_.LocalSize()) {
+  if (ls_ != GetHilbert().LocalSize()) {
     throw InvalidInputError(
         "Loaded wave-function has incompatible Hilbert space");
   }
diff --git a/Sources/Machine/rbm_multival.hpp b/Sources/Machine/rbm_multival.hpp
index e42d6e36..fda5472e 100644
--- a/Sources/Machine/rbm_multival.hpp
+++ b/Sources/Machine/rbm_multival.hpp
@@ -28,8 +28,6 @@ namespace netket {
 // Restricted Boltzman Machine wave function
 // for generic (finite) local hilbert space
 class RbmMultival : public AbstractMachine {
-  const AbstractHilbert &hilbert_;
-
   // number of visible units
   int nv_;
 
@@ -67,15 +65,14 @@ class RbmMultival : public AbstractMachine {
   std::map<double, int> confindex_;
 
  public:
-  explicit RbmMultival(const AbstractHilbert &hilbert, int nhidden = 0,
-                       int alpha = 0, bool usea = true, bool useb = true);
+  explicit RbmMultival(std::shared_ptr<const AbstractHilbert> hilbert,
+                       int nhidden = 0, int alpha = 0, bool usea = true,
+                       bool useb = true);
 
   int Npar() const override;
   int Nvisible() const override;
   /*constexpr*/ int Nhidden() const noexcept { return nh_; }
 
-  const AbstractHilbert &GetHilbert() const noexcept override;
-
   void InitRandomPars(int seed, double sigma) override;
   void InitLookup(VisibleConstType v, LookupType &lt) override;
   void UpdateLookup(VisibleConstType v, const std::vector<int> &tochange,
diff --git a/Sources/Machine/rbm_spin.cc b/Sources/Machine/rbm_spin.cc
index ded6b517..073b72ab 100644
--- a/Sources/Machine/rbm_spin.cc
+++ b/Sources/Machine/rbm_spin.cc
@@ -19,9 +19,9 @@
 
 namespace netket {
 
-RbmSpin::RbmSpin(const AbstractHilbert &hilbert, int nhidden, int alpha,
-                 bool usea, bool useb)
-    : hilbert_(hilbert), nv_(hilbert.Size()), usea_(usea), useb_(useb) {
+RbmSpin::RbmSpin(std::shared_ptr<const AbstractHilbert> hilbert, int nhidden,
+                 int alpha, bool usea, bool useb)
+    : AbstractMachine(hilbert), nv_(hilbert->Size()), usea_(usea), useb_(useb) {
   nh_ = std::max(nhidden, alpha * nv_);
   Init();
 }
@@ -219,8 +219,6 @@ Complex RbmSpin::LogValDiff(VisibleConstType v,
   return logvaldiff;
 }
 
-const AbstractHilbert &RbmSpin::GetHilbert() const noexcept { return hilbert_; }
-
 void RbmSpin::to_json(json &j) const {
   j["Name"] = "RbmSpin";
   j["Nvisible"] = nv_;
@@ -242,7 +240,7 @@ void RbmSpin::from_json(const json &pars) {
   if (FieldExists(pars, "Nvisible")) {
     nv_ = FieldVal<int>(pars, "Nvisible");
   }
-  if (nv_ != hilbert_.Size()) {
+  if (nv_ != GetHilbert().Size()) {
     throw InvalidInputError(
         "Number of visible units is incompatible with given "
         "Hilbert space");
diff --git a/Sources/Machine/rbm_spin.hpp b/Sources/Machine/rbm_spin.hpp
index 745cd1b8..bfe66dfd 100644
--- a/Sources/Machine/rbm_spin.hpp
+++ b/Sources/Machine/rbm_spin.hpp
@@ -27,8 +27,6 @@ using json = nlohmann::json;
  *
  */
 class RbmSpin : public AbstractMachine {
-  const AbstractHilbert &hilbert_;
-
   // number of visible units
   int nv_;
 
@@ -56,8 +54,8 @@ class RbmSpin : public AbstractMachine {
   bool useb_;
 
  public:
-  RbmSpin(const AbstractHilbert &hilbert, int nhidden = 0, int alpha = 0,
-          bool usea = true, bool useb = true);
+  RbmSpin(std::shared_ptr<const AbstractHilbert> hilbert, int nhidden = 0,
+          int alpha = 0, bool usea = true, bool useb = true);
 
   int Nvisible() const override;
   int Npar() const override;
@@ -89,8 +87,6 @@ class RbmSpin : public AbstractMachine {
                      const std::vector<double> &newconf,
                      const LookupType &lt) override;
 
-  const AbstractHilbert &GetHilbert() const noexcept override;
-
   void to_json(json &j) const override;
   void from_json(const json &pars) override;
 
diff --git a/Sources/Machine/rbm_spin_phase.cc b/Sources/Machine/rbm_spin_phase.cc
index a2357ccc..c1aafc4e 100644
--- a/Sources/Machine/rbm_spin_phase.cc
+++ b/Sources/Machine/rbm_spin_phase.cc
@@ -20,15 +20,14 @@
 
 namespace netket {
 
-RbmSpinPhase::RbmSpinPhase(const AbstractHilbert &hilbert, int nhidden,
-                           int alpha, bool usea, bool useb)
-    : hilbert_(hilbert),
-      nv_(hilbert.Size()),
+RbmSpinPhase::RbmSpinPhase(std::shared_ptr<const AbstractHilbert> hilbert,
+                           int nhidden, int alpha, bool usea, bool useb)
+    : AbstractMachine(hilbert),
+      nv_(hilbert->Size()),
       usea_(usea),
       useb_(useb),
       I_(0, 1) {
   nh_ = std::max(nhidden, alpha * nv_);
-
   Init();
 }
 
@@ -291,10 +290,6 @@ Complex RbmSpinPhase::LogValDiff(VisibleConstType v,
   return logvaldiff;
 }
 
-const AbstractHilbert &RbmSpinPhase::GetHilbert() const noexcept {
-  return hilbert_;
-}
-
 void RbmSpinPhase::to_json(json &j) const {
   j["Name"] = "RbmSpinPhase";
   j["Nvisible"] = nv_;
@@ -319,7 +314,7 @@ void RbmSpinPhase::from_json(const json &pars) {
   if (FieldExists(pars, "Nvisible")) {
     nv_ = FieldVal<int>(pars, "Nvisible");
   }
-  if (nv_ != hilbert_.Size()) {
+  if (nv_ != GetHilbert().Size()) {
     throw InvalidInputError(
         "Number of visible units is incompatible with given "
         "Hilbert space");
diff --git a/Sources/Machine/rbm_spin_phase.hpp b/Sources/Machine/rbm_spin_phase.hpp
index 65a2efb5..2fdfe631 100644
--- a/Sources/Machine/rbm_spin_phase.hpp
+++ b/Sources/Machine/rbm_spin_phase.hpp
@@ -27,8 +27,6 @@ amplitude
  *
  */
 class RbmSpinPhase : public AbstractMachine {
-  const AbstractHilbert &hilbert_;
-
   // number of visible units
   int nv_;
 
@@ -71,8 +69,8 @@ class RbmSpinPhase : public AbstractMachine {
   const Complex I_;
 
  public:
-  RbmSpinPhase(const AbstractHilbert &hilbert, int nhidden = 0, int alpha = 0,
-               bool usea = true, bool useb = true);
+  RbmSpinPhase(std::shared_ptr<const AbstractHilbert> hilbert, int nhidden = 0,
+               int alpha = 0, bool usea = true, bool useb = true);
 
   int Npar() const override;
   int Nvisible() const override;
@@ -108,7 +106,6 @@ class RbmSpinPhase : public AbstractMachine {
                      const LookupType &lt) override;
 
   bool IsHolomorphic() override;
-  const AbstractHilbert &GetHilbert() const noexcept override;
 
   void to_json(json &j) const override;
   void from_json(const json &pars) override;
diff --git a/Sources/Machine/rbm_spin_real.cc b/Sources/Machine/rbm_spin_real.cc
index 63aebc40..6dd15dad 100644
--- a/Sources/Machine/rbm_spin_real.cc
+++ b/Sources/Machine/rbm_spin_real.cc
@@ -20,11 +20,10 @@
 
 namespace netket {
 
-RbmSpinReal::RbmSpinReal(const AbstractHilbert &hilbert, int nhidden, int alpha,
-                         bool usea, bool useb)
-    : hilbert_(hilbert), nv_(hilbert.Size()), usea_(usea), useb_(useb) {
+RbmSpinReal::RbmSpinReal(std::shared_ptr<const AbstractHilbert> hilbert,
+                         int nhidden, int alpha, bool usea, bool useb)
+    : AbstractMachine(hilbert), nv_(hilbert->Size()), usea_(usea), useb_(useb) {
   nh_ = std::max(nhidden, alpha * nv_);
-
   Init();
 }
 
@@ -224,10 +223,6 @@ Complex RbmSpinReal::LogValDiff(VisibleConstType v,
   return logvaldiff;
 }
 
-const AbstractHilbert &RbmSpinReal::GetHilbert() const noexcept {
-  return hilbert_;
-}
-
 void RbmSpinReal::to_json(json &j) const {
   j["Name"] = "RbmSpinReal";
   j["Nvisible"] = nv_;
@@ -249,7 +244,7 @@ void RbmSpinReal::from_json(const json &pars) {
   if (FieldExists(pars, "Nvisible")) {
     nv_ = FieldVal<int>(pars, "Nvisible");
   }
-  if (nv_ != hilbert_.Size()) {
+  if (nv_ != GetHilbert().Size()) {
     throw InvalidInputError(
         "Number of visible units is incompatible with given "
         "Hilbert space");
diff --git a/Sources/Machine/rbm_spin_real.hpp b/Sources/Machine/rbm_spin_real.hpp
index c64c2a8a..1373461b 100644
--- a/Sources/Machine/rbm_spin_real.hpp
+++ b/Sources/Machine/rbm_spin_real.hpp
@@ -26,8 +26,6 @@ and real-valued weights
  *
  */
 class RbmSpinReal : public AbstractMachine {
-  const AbstractHilbert &hilbert_;
-
   // number of visible units
   int nv_;
 
@@ -55,8 +53,8 @@ class RbmSpinReal : public AbstractMachine {
   bool useb_;
 
  public:
-  RbmSpinReal(const AbstractHilbert &hilbert, int nhidden = 0, int alpha = 0,
-              bool usea = true, bool useb = true);
+  RbmSpinReal(std::shared_ptr<const AbstractHilbert> hilbert, int nhidden = 0,
+              int alpha = 0, bool usea = true, bool useb = true);
 
   int Npar() const override;
   int Nvisible() const override;
@@ -83,8 +81,6 @@ class RbmSpinReal : public AbstractMachine {
                      const std::vector<double> &newconf,
                      const LookupType &lt) override;
 
-  const AbstractHilbert &GetHilbert() const noexcept override;
-
   void to_json(json &j) const override;
   void from_json(const json &pars) override;
 
diff --git a/Sources/Machine/rbm_spin_symm.cc b/Sources/Machine/rbm_spin_symm.cc
index 148817bb..e2a7f5ef 100644
--- a/Sources/Machine/rbm_spin_symm.cc
+++ b/Sources/Machine/rbm_spin_symm.cc
@@ -20,16 +20,15 @@
 
 namespace netket {
 
-RbmSpinSymm::RbmSpinSymm(const AbstractHilbert &hilbert, int alpha, bool usea,
-                         bool useb)
-    : hilbert_(hilbert),
-      graph_(hilbert.GetGraph()),
-      nv_(hilbert.Size()),
+RbmSpinSymm::RbmSpinSymm(std::shared_ptr<const AbstractHilbert> hilbert,
+                         int alpha, bool usea, bool useb)
+    : AbstractMachine(hilbert),
+      graph_(hilbert->GetGraph()),
+      nv_(hilbert->Size()),
       alpha_(alpha),
       usea_(usea),
       useb_(useb) {
   Init(graph_);
-
   SetBareParameters();
 }
 
@@ -347,10 +346,6 @@ Complex RbmSpinSymm::LogValDiff(VisibleConstType v,
   return logvaldiff;
 }
 
-const AbstractHilbert &RbmSpinSymm::GetHilbert() const noexcept {
-  return hilbert_;
-}
-
 void RbmSpinSymm::to_json(json &j) const {
   j["Name"] = "RbmSpinSymm";
   j["Nvisible"] = nv_;
@@ -371,7 +366,7 @@ void RbmSpinSymm::from_json(const json &pars) {
   if (FieldExists(pars, "Nvisible")) {
     nv_ = pars["Nvisible"];
   }
-  if (nv_ != hilbert_.Size()) {
+  if (nv_ != GetHilbert().Size()) {
     throw InvalidInputError(
         "Number of visible units is incompatible with given "
         "Hilbert space");
diff --git a/Sources/Machine/rbm_spin_symm.hpp b/Sources/Machine/rbm_spin_symm.hpp
index 017f05ea..aed71922 100644
--- a/Sources/Machine/rbm_spin_symm.hpp
+++ b/Sources/Machine/rbm_spin_symm.hpp
@@ -23,8 +23,6 @@ using json = nlohmann::json;
 
 // Rbm with permutation symmetries
 class RbmSpinSymm : public AbstractMachine {
-  const AbstractHilbert &hilbert_;
-
   const AbstractGraph &graph_;
 
   // number of visible units
@@ -72,15 +70,12 @@ class RbmSpinSymm : public AbstractMachine {
   bool useb_;
 
  public:
-  RbmSpinSymm(const AbstractHilbert &hilbert, int alpha = 0, bool usea = true,
-              bool useb = true);
+  RbmSpinSymm(std::shared_ptr<const AbstractHilbert> hilbert, int alpha = 0,
+              bool usea = true, bool useb = true);
 
   int Npar() const override;
   int Nvisible() const override;
   int Nhidden() const { return nh_; }
-
-  const AbstractHilbert &GetHilbert() const noexcept override;
-
   void InitRandomPars(int seed, double sigma) override;
   void InitLookup(VisibleConstType v, LookupType &lt) override;
   void UpdateLookup(VisibleConstType v, const std::vector<int> &tochange,
diff --git a/Sources/Operator/MatrixWrapper/dense_matrix_wrapper.hpp b/Sources/Operator/MatrixWrapper/dense_matrix_wrapper.hpp
index ab5ddb3f..4e5fa76c 100644
--- a/Sources/Operator/MatrixWrapper/dense_matrix_wrapper.hpp
+++ b/Sources/Operator/MatrixWrapper/dense_matrix_wrapper.hpp
@@ -80,8 +80,7 @@ class DenseMatrixWrapper : public AbstractMatrixWrapper<State> {
 
  private:
   void InitializeMatrix(const AbstractOperator& the_operator) {
-    const auto& hilbert = the_operator.GetHilbert();
-    const auto& hilbert_index = hilbert.GetIndex();
+    const auto& hilbert_index = the_operator.GetHilbert().GetIndex();
     dim_ = hilbert_index.NStates();
 
     matrix_.resize(dim_, dim_);
diff --git a/Sources/Operator/MatrixWrapper/sparse_matrix_wrapper.hpp b/Sources/Operator/MatrixWrapper/sparse_matrix_wrapper.hpp
index 77e08f8c..2eb882f8 100644
--- a/Sources/Operator/MatrixWrapper/sparse_matrix_wrapper.hpp
+++ b/Sources/Operator/MatrixWrapper/sparse_matrix_wrapper.hpp
@@ -72,8 +72,7 @@ public:
 
 private:
   void InitializeMatrix(const AbstractOperator &the_operator) {
-    const auto &hilbert = the_operator.GetHilbert();
-    const auto &hilbert_index = hilbert.GetIndex();
+    const auto& hilbert_index = the_operator.GetHilbert().GetIndex();
     dim_ = hilbert_index.NStates();
 
     using Triplet = Eigen::Triplet<Complex>;
diff --git a/Sources/Operator/abstract_operator.hpp b/Sources/Operator/abstract_operator.hpp
index d931bc61..be1cb54e 100644
--- a/Sources/Operator/abstract_operator.hpp
+++ b/Sources/Operator/abstract_operator.hpp
@@ -100,9 +100,19 @@ class AbstractOperator {
   Member function returning the hilbert space associated with this Hamiltonian.
   @return Hilbert space specifier for this Hamiltonian
   */
-  virtual const AbstractHilbert &GetHilbert() const noexcept = 0;
+  const AbstractHilbert &GetHilbert() const { return *hilbert_; }
+  std::shared_ptr<const AbstractHilbert> GetHilbertShared() const {
+    return hilbert_;
+  }
+
+  virtual ~AbstractOperator() = default;
+
+ protected:
+  AbstractOperator(std::shared_ptr<const AbstractHilbert> hilbert)
+      : hilbert_(std::move(hilbert)) {}
 
-  virtual ~AbstractOperator() {}
+ private:
+  std::shared_ptr<const AbstractHilbert> hilbert_;
 };
 
 void AbstractOperator::ForEachConn(VectorConstRefType v,
diff --git a/Sources/Operator/bosonhubbard.hpp b/Sources/Operator/bosonhubbard.hpp
index 20a8a580..c947eb5b 100644
--- a/Sources/Operator/bosonhubbard.hpp
+++ b/Sources/Operator/bosonhubbard.hpp
@@ -18,7 +18,9 @@
 #include <Eigen/Dense>
 #include <cmath>
 #include <iostream>
+#include <memory>
 #include <vector>
+
 #include "Graph/graph.hpp"
 #include "Hilbert/abstract_hilbert.hpp"
 #include "Utils/exceptions.hpp"
@@ -29,7 +31,6 @@ namespace netket {
 
 // Heisenberg model on an arbitrary graph
 class BoseHubbard : public AbstractOperator {
-  const AbstractHilbert &hilbert_;
   const AbstractGraph &graph_;
 
   int nsites_;
@@ -50,15 +51,15 @@ class BoseHubbard : public AbstractOperator {
   using VectorRefType = AbstractOperator::VectorRefType;
   using VectorConstRefType = AbstractOperator::VectorConstRefType;
 
-  explicit BoseHubbard(const AbstractHilbert &hilbert, double U, double V = 0.,
-                       double mu = 0.)
-      : hilbert_(hilbert),
-        graph_(hilbert.GetGraph()),
-        nsites_(hilbert.Size()),
+  BoseHubbard(std::shared_ptr<const AbstractHilbert> hilbert, double U,
+              double V = 0., double mu = 0.)
+      : AbstractOperator(hilbert),
+        graph_(hilbert->GetGraph()),
+        nsites_(hilbert->Size()),
         U_(U),
         V_(V),
         mu_(mu) {
-    nmax_ = hilbert_.LocalSize() - 1;
+    nmax_ = hilbert->LocalSize() - 1;
     Init();
   }
 
@@ -122,10 +123,6 @@ class BoseHubbard : public AbstractOperator {
       }
     }
   }
-
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
 };
 
 }  // namespace netket
diff --git a/Sources/Operator/graph_operator.hpp b/Sources/Operator/graph_operator.hpp
index 61666378..720098c2 100644
--- a/Sources/Operator/graph_operator.hpp
+++ b/Sources/Operator/graph_operator.hpp
@@ -28,8 +28,6 @@ namespace netket {
 
 // Graph Hamiltonian on an arbitrary graph
 class GraphOperator : public AbstractOperator {
-  const AbstractHilbert &hilbert_;
-
   // Arbitrary graph
   const AbstractGraph &graph_;
 
@@ -45,14 +43,13 @@ class GraphOperator : public AbstractOperator {
   using VectorRefType = AbstractOperator::VectorRefType;
   using VectorConstRefType = AbstractOperator::VectorConstRefType;
 
-  explicit GraphOperator(const AbstractHilbert &hilbert,
-                         OVecType siteops = OVecType(),
-                         OVecType bondops = OVecType(),
-                         std::vector<int> bondops_colors = std::vector<int>())
-      : hilbert_(hilbert),
-        graph_(hilbert.GetGraph()),
+  GraphOperator(std::shared_ptr<const AbstractHilbert> hilbert,
+                OVecType siteops = OVecType(), OVecType bondops = OVecType(),
+                std::vector<int> bondops_colors = std::vector<int>())
+      : AbstractOperator(hilbert),
+        graph_(hilbert->GetGraph()),
         operator_(hilbert),
-        nvertices_(hilbert.Size()) {
+        nvertices_(hilbert->Size()) {
     // Create the local operator as the sum of all site and bond operators
     // Ensure that at least one of SiteOps and BondOps was initialized
     if (!siteops.size() && !bondops.size()) {
@@ -70,7 +67,7 @@ class GraphOperator : public AbstractOperator {
     if (siteops.size() > 0) {
       for (int i = 0; i < nvertices_; i++) {
         for (std::size_t j = 0; j < siteops.size(); j++) {
-          operator_ += LocalOperator(hilbert_, siteops[j], std::vector<int>{i});
+          operator_ += LocalOperator(hilbert, siteops[j], std::vector<int>{i});
         }
       }
     }
@@ -88,7 +85,7 @@ class GraphOperator : public AbstractOperator {
         for (std::size_t c = 0; c < op_color.size(); c++) {
           if (op_color[c] == kv.second && kv.first[0] < kv.first[1]) {
             std::vector<int> edge = {kv.first[0], kv.first[1]};
-            operator_ += LocalOperator(hilbert_, bondops[c], edge);
+            operator_ += LocalOperator(hilbert, bondops[c], edge);
           }
         }
       }
@@ -96,12 +93,12 @@ class GraphOperator : public AbstractOperator {
   }
 
   // Constructor to be used when overloading operators
-  explicit GraphOperator(const AbstractHilbert &hilbert,
+  explicit GraphOperator(std::shared_ptr<const AbstractHilbert> hilbert,
                          const LocalOperator &lop)
-      : hilbert_(hilbert),
-        graph_(hilbert.GetGraph()),
+      : AbstractOperator(hilbert),
+        graph_(hilbert->GetGraph()),
         operator_(lop),
-        nvertices_(hilbert.Size()) {}
+        nvertices_(hilbert->Size()) {}
 
   friend GraphOperator operator+(const GraphOperator &lhs,
                                  const GraphOperator &rhs) {
@@ -110,7 +107,7 @@ class GraphOperator : public AbstractOperator {
     auto lop = lhs.operator_;
     auto rop = rhs.operator_;
 
-    return GraphOperator(lhs.GetHilbert(), lop + rop);
+    return GraphOperator(lhs.GetHilbertShared(), lop + rop);
   }
 
   void FindConn(VectorConstRefType v, std::vector<Complex> &mel,
@@ -122,10 +119,6 @@ class GraphOperator : public AbstractOperator {
   void ForEachConn(VectorConstRefType v, ConnCallback callback) const override {
     operator_.ForEachConn(v, callback);
   }
-
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
 };  // namespace netket
 }  // namespace netket
 #endif
diff --git a/Sources/Operator/local_operator.hpp b/Sources/Operator/local_operator.hpp
index a4a269d1..7fa76346 100644
--- a/Sources/Operator/local_operator.hpp
+++ b/Sources/Operator/local_operator.hpp
@@ -51,7 +51,6 @@ class LocalOperator : public AbstractOperator {
   using VectorConstRefType = AbstractOperator::VectorConstRefType;
 
  private:
-  const AbstractHilbert &hilbert_;
   std::vector<MatType> mat_;
   std::vector<SiteType> sites_;
 
@@ -66,23 +65,34 @@ class LocalOperator : public AbstractOperator {
   static constexpr double mel_cutoff_ = 1.0e-6;
 
  public:
-  explicit LocalOperator(const AbstractHilbert &hilbert, double constant = 0.)
-      : hilbert_(hilbert), constant_(constant), nops_(0) {}
+  // explicit LocalOperator(const LocalOperator &rhs)
+  //     : mat_(rhs.mat_),
+  //       sites_(rhs.sites_),
+  //       invstate_(rhs.invstate_),
+  //       states_(rhs.states_),
+  //       connected_(rhs.connected_) {
+  //   SetHilbert(rhs.GetHilbert());
+  // }
+
+  explicit LocalOperator(std::shared_ptr<const AbstractHilbert> hilbert,
+                         double constant = 0.)
+      : AbstractOperator(hilbert), constant_(constant), nops_(0) {}
 
-  explicit LocalOperator(const AbstractHilbert &hilbert,
+  explicit LocalOperator(std::shared_ptr<const AbstractHilbert> hilbert,
                          const std::vector<MatType> &mat,
                          const std::vector<SiteType> &sites,
                          double constant = 0.)
-      : hilbert_(hilbert), constant_(constant) {
+      : AbstractOperator(hilbert), constant_(constant) {
     for (std::size_t i = 0; i < mat.size(); i++) {
       Push(mat[i], sites[i]);
     }
     Init();
   }
 
-  explicit LocalOperator(const AbstractHilbert &hilbert, const MatType &mat,
-                         const SiteType &sites, double constant = 0.)
-      : hilbert_(hilbert), constant_(constant) {
+  explicit LocalOperator(std::shared_ptr<const AbstractHilbert> hilbert,
+                         const MatType &mat, const SiteType &sites,
+                         double constant = 0.)
+      : AbstractOperator(hilbert), constant_(constant) {
     Push(mat, sites);
     // TODO sort sites and swap columns of mat accordingly
     Init();
@@ -104,7 +114,7 @@ class LocalOperator : public AbstractOperator {
   }
 
   void Init() {
-    if (!hilbert_.IsDiscrete()) {
+    if (!GetHilbert().IsDiscrete()) {
       throw InvalidInputError(
           "Cannot construct operators on infinite local hilbert spaces");
     }
@@ -130,12 +140,13 @@ class LocalOperator : public AbstractOperator {
       auto &states = states_[op];
       auto &invstate = invstate_[op];
 
-      if (*std::max_element(sites.begin(), sites.end()) >= hilbert_.Size() ||
+      if (*std::max_element(sites.begin(), sites.end()) >=
+              GetHilbert().Size() ||
           *std::min_element(sites.begin(), sites.end()) < 0) {
         throw InvalidInputError("Operator acts on an invalid set of sites");
       }
 
-      auto localstates = hilbert_.LocalStates();
+      auto localstates = GetHilbert().LocalStates();
       const auto localsize = localstates.size();
 
       // Finding the non-zero matrix elements
@@ -190,7 +201,7 @@ class LocalOperator : public AbstractOperator {
   void FindConn(VectorConstRefType v, std::vector<Complex> &mel,
                 std::vector<std::vector<int>> &connectors,
                 std::vector<std::vector<double>> &newconfs) const override {
-    assert(v.size() == hilbert_.Size());
+    assert(v.size() == GetHilbert().Size());
 
     connectors.clear();
     newconfs.clear();
@@ -223,7 +234,7 @@ class LocalOperator : public AbstractOperator {
   }
 
   void ForEachConn(VectorConstRefType v, ConnCallback callback) const override {
-    assert(v.size() == hilbert_.Size());
+    assert(v.size() == GetHilbert().Size());
 
     Complex mel_diag = 0.;
 
@@ -287,7 +298,7 @@ class LocalOperator : public AbstractOperator {
     std::transform(mat_.begin(), mat_.end(), mat_t.begin(),
                    netket::transpose_vecvec<MelType>);
 
-    return LocalOperator(hilbert_, mat_t, sites_, constant_);
+    return LocalOperator(GetHilbertShared(), mat_t, sites_, constant_);
   }
 
   LocalOperator Conjugate() const {
@@ -302,7 +313,7 @@ class LocalOperator : public AbstractOperator {
                      }
                      return out;
                    });
-    return LocalOperator(hilbert_, mat_c, sites_, constant_);
+    return LocalOperator(GetHilbertShared(), mat_c, sites_, constant_);
   }
 
   // Product of two local operators, performing KroneckerProducts as necessary
@@ -330,7 +341,7 @@ class LocalOperator : public AbstractOperator {
       }
     }
     auto constant = lhs.constant_ * rhs.constant_;
-    auto opret = LocalOperator(lhs.GetHilbert(), mat, sites, constant);
+    auto opret = LocalOperator(lhs.GetHilbertShared(), mat, sites, constant);
 
     if (std::abs(lhs.constant_) > mel_cutoff_) {
       opret += lhs.constant_ * rhs;
@@ -344,8 +355,8 @@ class LocalOperator : public AbstractOperator {
 
   friend LocalOperator operator+(const LocalOperator &lhs,
                                  const LocalOperator &rhs) {
-    assert(rhs.hilbert_.LocalStates().size() ==
-           lhs.hilbert_.LocalStates().size());
+    assert(rhs.GetHilbert().LocalStates().size() ==
+           lhs.GetHilbert().LocalStates().size());
 
     auto sites = lhs.sites_;
     auto mat = lhs.mat_;
@@ -353,7 +364,7 @@ class LocalOperator : public AbstractOperator {
     sites.insert(sites.end(), rhs.sites_.begin(), rhs.sites_.end());
     mat.insert(mat.end(), rhs.mat_.begin(), rhs.mat_.end());
 
-    return LocalOperator(lhs.GetHilbert(), mat, sites,
+    return LocalOperator(lhs.GetHilbertShared(), mat, sites,
                          lhs.constant_ + rhs.constant_);
   }
 
@@ -361,13 +372,13 @@ class LocalOperator : public AbstractOperator {
     auto sites = lhs.sites_;
     auto mat = lhs.mat_;
 
-    return LocalOperator(lhs.GetHilbert(), mat, sites,
+    return LocalOperator(lhs.GetHilbertShared(), mat, sites,
                          lhs.constant_ + constant);
   }
 
   LocalOperator &operator+=(const LocalOperator &rhs) {
-    assert(rhs.hilbert_.LocalStates().size() ==
-           this->hilbert_.LocalStates().size());
+    assert(rhs.GetHilbert().LocalStates().size() ==
+           this->GetHilbert().LocalStates().size());
 
     this->sites_.insert(this->sites_.end(), rhs.sites_.begin(),
                         rhs.sites_.end());
@@ -396,17 +407,13 @@ class LocalOperator : public AbstractOperator {
       }
     }
 
-    return LocalOperator(rhs.GetHilbert(), mat, sites,
+    return LocalOperator(rhs.GetHilbertShared(), mat, sites,
                          std::real(lhs * rhs.constant_));
   }
 
   const std::vector<MatType> &LocalMatrices() const { return mat_; }
   const std::vector<SiteType> &ActingOn() const { return sites_; }
 
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   std::size_t Size() const { return mat_.size(); }
 };  // namespace netket
 
diff --git a/Sources/Operator/py_bosonhubbard.hpp b/Sources/Operator/py_bosonhubbard.hpp
index 08cb3425..f8d51bc9 100644
--- a/Sources/Operator/py_bosonhubbard.hpp
+++ b/Sources/Operator/py_bosonhubbard.hpp
@@ -32,7 +32,8 @@ void AddBoseHubbard(py::module &subm) {
   py::class_<BoseHubbard, AbstractOperator>(
       subm, "BoseHubbard",
       R"EOF(A Bose Hubbard model Hamiltonian operator.)EOF")
-      .def(py::init<const AbstractHilbert &, double, double, double>(),
+      .def(py::init<std::shared_ptr<const AbstractHilbert>, double, double,
+                    double>(),
            py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("U"),
            py::arg("V") = 0., py::arg("mu") = 0., R"EOF(
            Constructs a new ``BoseHubbard`` given a hilbert space and a Hubbard
diff --git a/Sources/Operator/py_graph_operator.hpp b/Sources/Operator/py_graph_operator.hpp
index e970eda5..2f70e4e3 100644
--- a/Sources/Operator/py_graph_operator.hpp
+++ b/Sources/Operator/py_graph_operator.hpp
@@ -31,8 +31,9 @@ namespace netket {
 void AddGraphOperator(py::module &subm) {
   py::class_<GraphOperator, AbstractOperator>(
       subm, "GraphOperator", R"EOF(A custom graph based operator.)EOF")
-      .def(py::init<const AbstractHilbert &, GraphOperator::OVecType,
-                    GraphOperator::OVecType, std::vector<int>>(),
+      .def(py::init<std::shared_ptr<const AbstractHilbert>,
+                    GraphOperator::OVecType, GraphOperator::OVecType,
+                    std::vector<int>>(),
            py::keep_alive<1, 2>(), py::arg("hilbert"),
            py::arg("siteops") = GraphOperator::OVecType(),
            py::arg("bondops") = GraphOperator::OVecType(),
diff --git a/Sources/Operator/py_local_operator.hpp b/Sources/Operator/py_local_operator.hpp
index 234ecdbd..94120014 100644
--- a/Sources/Operator/py_local_operator.hpp
+++ b/Sources/Operator/py_local_operator.hpp
@@ -31,8 +31,9 @@ namespace netket {
 void AddLocalOperator(py::module &subm) {
   py::class_<LocalOperator, AbstractOperator>(
       subm, "LocalOperator", R"EOF(A custom local operator.)EOF")
-      .def(py::init<const AbstractHilbert &, double>(), py::keep_alive<1, 2>(),
-           py::arg("hilbert"), py::arg("constant") = 0., R"EOF(
+      .def(py::init<std::shared_ptr<const AbstractHilbert>, double>(),
+           py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("constant") = 0.,
+           R"EOF(
            Constructs a new ``LocalOperator`` given a hilbert space and (if
            specified) a constant level shift.
 
@@ -55,11 +56,11 @@ void AddLocalOperator(py::module &subm) {
 
                ```
            )EOF")
-      .def(
-          py::init<const AbstractHilbert &, std::vector<LocalOperator::MatType>,
-                   std::vector<LocalOperator::SiteType>, double>(),
-          py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("operators"),
-          py::arg("acting_on"), py::arg("constant") = 0., R"EOF(
+      .def(py::init<std::shared_ptr<const AbstractHilbert>,
+                    std::vector<LocalOperator::MatType>,
+                    std::vector<LocalOperator::SiteType>, double>(),
+           py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("operators"),
+           py::arg("acting_on"), py::arg("constant") = 0., R"EOF(
           Constructs a new ``LocalOperator`` given a hilbert space, a vector of
           operators, a vector of sites, and (if specified) a constant level
           shift.
@@ -88,8 +89,8 @@ void AddLocalOperator(py::module &subm) {
 
               ```
           )EOF")
-      .def(py::init<const AbstractHilbert &, LocalOperator::MatType,
-                    LocalOperator::SiteType, double>(),
+      .def(py::init<std::shared_ptr<const AbstractHilbert>,
+                    LocalOperator::MatType, LocalOperator::SiteType, double>(),
            py::keep_alive<1, 2>(), py::arg("hilbert"), py::arg("operator"),
            py::arg("acting_on"), py::arg("constant") = 0., R"EOF(
            Constructs a new ``LocalOperator`` given a hilbert space, an
@@ -131,22 +132,30 @@ void AddLocalOperator(py::module &subm) {
       .def("conjugate", &LocalOperator::Conjugate,
            R"EOF(Returns the complex conjugation of this operator)EOF")
       .def(py::self + py::self)
-      .def("__mul__", [](const LocalOperator &a, double b) { return b * a; },
-           py::is_operator())
-      .def("__rmul__", [](const LocalOperator &a, double b) { return b * a; },
-           py::is_operator())
-      .def("__mul__", [](const LocalOperator &a, int b) { return b * a; },
-           py::is_operator())
-      .def("__rmul__", [](const LocalOperator &a, int b) { return b * a; },
-           py::is_operator())
-      .def("__add__", [](const LocalOperator &a, double b) { return a + b; },
-           py::is_operator())
-      .def("__add__", [](const LocalOperator &a, int b) { return a + b; },
-           py::is_operator())
-      .def("__radd__", [](const LocalOperator &a, double b) { return a + b; },
-           py::is_operator())
-      .def("__radd__", [](const LocalOperator &a, int b) { return a + b; },
-           py::is_operator())
+      .def(
+          "__mul__", [](const LocalOperator &a, double b) { return b * a; },
+          py::is_operator())
+      .def(
+          "__rmul__", [](const LocalOperator &a, double b) { return b * a; },
+          py::is_operator())
+      .def(
+          "__mul__", [](const LocalOperator &a, int b) { return b * a; },
+          py::is_operator())
+      .def(
+          "__rmul__", [](const LocalOperator &a, int b) { return b * a; },
+          py::is_operator())
+      .def(
+          "__add__", [](const LocalOperator &a, double b) { return a + b; },
+          py::is_operator())
+      .def(
+          "__add__", [](const LocalOperator &a, int b) { return a + b; },
+          py::is_operator())
+      .def(
+          "__radd__", [](const LocalOperator &a, double b) { return a + b; },
+          py::is_operator())
+      .def(
+          "__radd__", [](const LocalOperator &a, int b) { return a + b; },
+          py::is_operator())
       .def(py::self * py::self);
 }
 
diff --git a/Sources/Operator/py_operator.hpp b/Sources/Operator/py_operator.hpp
index e663ab61..f30b69ee 100644
--- a/Sources/Operator/py_operator.hpp
+++ b/Sources/Operator/py_operator.hpp
@@ -56,11 +56,12 @@ void AddOperatorModule(py::module &m) {
           .def_property_readonly(
               "hilbert", &AbstractOperator::GetHilbert,
               R"EOF(netket.hilbert.Hilbert: ``Hilbert`` space of operator.)EOF")
-          .def("to_sparse",
-               [](const AbstractOperator &self) {
-                 return SparseMatrixWrapper<>(self).GetMatrix();
-               },
-               R"EOF(
+          .def(
+              "to_sparse",
+              [](const AbstractOperator &self) {
+                return SparseMatrixWrapper<>(self).GetMatrix();
+              },
+              R"EOF(
          Returns the sparse matrix representation of the operator. Note that, in general,
          the size of the matrix is exponential in the number of quantum
          numbers, and this operation should thus only be performed for
@@ -68,11 +69,12 @@ void AddOperatorModule(py::module &m) {
 
          This method requires an indexable Hilbert space.
          )EOF")
-          .def("to_dense",
-               [](const AbstractOperator &self) {
-                 return DenseMatrixWrapper<>(self).GetMatrix();
-               },
-               R"EOF(
+          .def(
+              "to_dense",
+              [](const AbstractOperator &self) {
+                return DenseMatrixWrapper<>(self).GetMatrix();
+              },
+              R"EOF(
          Returns the dense matrix representation of the operator. Note that, in general,
          the size of the matrix is exponential in the number of quantum
          numbers, and this operation should thus only be performed for
diff --git a/Sources/Sampler/abstract_sampler.hpp b/Sources/Sampler/abstract_sampler.hpp
index afdc2dee..395d15e4 100644
--- a/Sources/Sampler/abstract_sampler.hpp
+++ b/Sources/Sampler/abstract_sampler.hpp
@@ -34,12 +34,8 @@ class AbstractSampler {
 
   virtual void SetVisible(const Eigen::VectorXd& v) = 0;
 
-  virtual AbstractMachine& GetMachine() noexcept = 0;
-
   virtual Eigen::VectorXd Acceptance() const = 0;
 
-  virtual const AbstractHilbert& GetHilbert() const noexcept = 0;
-
   // Computes the derivative of the machine on the current visible
   // Using the lookUp tables if possible
   virtual AbstractMachine::VectorType DerLogVisible() {
@@ -61,12 +57,25 @@ class AbstractSampler {
     machine_func_ = std::move(machine_func);
   }
 
+  std::shared_ptr<const AbstractHilbert> GetHilbertShared() const noexcept {
+    return hilbert_;
+  }
+
+  const AbstractHilbert &GetHilbert() const noexcept {
+    return *hilbert_;
+  }
+
+  AbstractMachine &GetMachine() const noexcept {
+    return psi_;
+  }
+
   const MachineFunction& GetMachineFunc() const noexcept {
     return machine_func_;
   }
 
  protected:
-  AbstractSampler() {
+  AbstractSampler(AbstractMachine& psi)
+      : psi_(psi), hilbert_(psi.GetHilbertShared()) {
     // Default initialization for the machine function to be sampled from
     machine_func_ = static_cast<double (*)(const Complex&)>(&std::norm);
   }
@@ -76,6 +85,8 @@ class AbstractSampler {
  private:
   DistributedRandomEngine engine_;
   MachineFunction machine_func_;
+  AbstractMachine& psi_;
+  std::shared_ptr<const AbstractHilbert> hilbert_;
 };
 
 }  // namespace netket
diff --git a/Sources/Sampler/custom_sampler.hpp b/Sources/Sampler/custom_sampler.hpp
index 7cf1f14e..a29ac671 100644
--- a/Sources/Sampler/custom_sampler.hpp
+++ b/Sources/Sampler/custom_sampler.hpp
@@ -30,8 +30,6 @@ namespace netket {
 
 // Metropolis sampling using custom moves provided by user
 class CustomSampler : public AbstractSampler {
-  AbstractMachine &psi_;
-  const AbstractHilbert &hilbert_;
   LocalOperator move_operators_;
   std::vector<double> operatorsweights_;
 
@@ -60,17 +58,16 @@ class CustomSampler : public AbstractSampler {
  public:
   CustomSampler(AbstractMachine &psi, const LocalOperator &move_operators,
                 const std::vector<double> &move_weights = {})
-      : psi_(psi),
-        hilbert_(psi.GetHilbert()),
+      : AbstractSampler(psi),
         move_operators_(move_operators),
-        nv_(hilbert_.Size()) {
+        nv_(GetHilbert().Size()) {
     Init(move_weights);
   }
 
   void Init(const std::vector<double> &move_weights) {
     CheckMoveOperators(move_operators_);
 
-    if (hilbert_.Size() != move_operators_.GetHilbert().Size()) {
+    if (GetHilbert().Size() != move_operators_.GetHilbert().Size()) {
       throw InvalidInputError(
           "Move operators in CustomSampler act on a different hilbert space "
           "than the Machine");
@@ -93,7 +90,7 @@ class CustomSampler : public AbstractSampler {
     MPI_Comm_size(MPI_COMM_WORLD, &totalnodes_);
     MPI_Comm_rank(MPI_COMM_WORLD, &mynode_);
 
-    if (!hilbert_.IsDiscrete()) {
+    if (!GetHilbert().IsDiscrete()) {
       throw InvalidInputError(
           "Custom Metropolis sampler works only for discrete Hilbert spaces");
     }
@@ -101,8 +98,8 @@ class CustomSampler : public AbstractSampler {
     accept_.resize(1);
     moves_.resize(1);
 
-    nstates_ = hilbert_.LocalSize();
-    localstates_ = hilbert_.LocalStates();
+    nstates_ = GetHilbert().LocalSize();
+    localstates_ = GetHilbert().LocalStates();
 
     Reset(true);
 
@@ -111,10 +108,10 @@ class CustomSampler : public AbstractSampler {
 
   void Reset(bool initrandom = false) override {
     if (initrandom) {
-      hilbert_.RandomVals(v_, this->GetRandomEngine());
+      GetHilbert().RandomVals(v_, this->GetRandomEngine());
     }
 
-    psi_.InitLookup(v_, lt_);
+    GetMachine().InitLookup(v_, lt_);
 
     accept_ = Eigen::VectorXd::Zero(1);
     moves_ = Eigen::VectorXd::Zero(1);
@@ -140,16 +137,16 @@ class CustomSampler : public AbstractSampler {
         cumulative_prob += std::real(mel_[exit_state]);
       }
 
-      auto exlog = std::exp(psi_.LogValDiff(v_, tochange_[exit_state],
+      auto exlog = std::exp(GetMachine().LogValDiff(v_, tochange_[exit_state],
                                             newconfs_[exit_state], lt_));
       double ratio = this->GetMachineFunc()(exlog);
 
       // Metropolis acceptance test
       if (ratio > distu(this->GetRandomEngine())) {
         accept_[0] += 1;
-        psi_.UpdateLookup(v_, tochange_[exit_state], newconfs_[exit_state],
+        GetMachine().UpdateLookup(v_, tochange_[exit_state], newconfs_[exit_state],
                           lt_);
-        hilbert_.UpdateConf(v_, tochange_[exit_state], newconfs_[exit_state]);
+        GetHilbert().UpdateConf(v_, tochange_[exit_state], newconfs_[exit_state]);
       }
       moves_[0] += 1;
     }
@@ -159,14 +156,8 @@ class CustomSampler : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd &v) override { v_ = v; }
 
-  AbstractMachine &GetMachine() noexcept override { return psi_; }
-
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   AbstractMachine::VectorType DerLogVisible() override {
-    return psi_.DerLog(v_, lt_);
+    return GetMachine().DerLog(v_, lt_);
   }
 
   Eigen::VectorXd Acceptance() const override {
diff --git a/Sources/Sampler/custom_sampler_pt.hpp b/Sources/Sampler/custom_sampler_pt.hpp
index 6b67c5db..2919d043 100644
--- a/Sources/Sampler/custom_sampler_pt.hpp
+++ b/Sources/Sampler/custom_sampler_pt.hpp
@@ -30,8 +30,6 @@ namespace netket {
 
 // Metropolis sampling using custom moves provided by user
 class CustomSamplerPt : public AbstractSampler {
-  AbstractMachine& psi_;
-  const AbstractHilbert& hilbert_;
   LocalOperator move_operators_;
   std::vector<double> operatorsweights_;
   // number of visible units
@@ -64,10 +62,9 @@ class CustomSamplerPt : public AbstractSampler {
   CustomSamplerPt(AbstractMachine& psi, const LocalOperator& move_operators,
                   const std::vector<double>& move_weights = {},
                   int nreplicas = 1)
-      : psi_(psi),
-        hilbert_(psi.GetHilbert()),
+      : AbstractSampler(psi),
         move_operators_(move_operators),
-        nv_(hilbert_.Size()),
+        nv_(GetHilbert().Size()),
         nrep_(nreplicas) {
     Init(move_weights);
   }
@@ -75,7 +72,7 @@ class CustomSamplerPt : public AbstractSampler {
   void Init(const std::vector<double>& move_weights) {
     CustomSampler::CheckMoveOperators(move_operators_);
 
-    if (hilbert_.Size() != move_operators_.GetHilbert().Size()) {
+    if (GetHilbert().Size() != move_operators_.GetHilbert().Size()) {
       throw InvalidInputError(
           "Move operators in CustomSampler act on a different hilbert space "
           "than the Machine");
@@ -97,13 +94,13 @@ class CustomSamplerPt : public AbstractSampler {
     MPI_Comm_size(MPI_COMM_WORLD, &totalnodes_);
     MPI_Comm_rank(MPI_COMM_WORLD, &mynode_);
 
-    if (!hilbert_.IsDiscrete()) {
+    if (!GetHilbert().IsDiscrete()) {
       throw InvalidInputError(
           "Custom Metropolis sampler works only for discrete Hilbert spaces");
     }
 
-    nstates_ = hilbert_.LocalSize();
-    localstates_ = hilbert_.LocalStates();
+    nstates_ = GetHilbert().LocalSize();
+    localstates_ = GetHilbert().LocalStates();
 
     v_.resize(nrep_);
     for (int i = 0; i < nrep_; i++) {
@@ -130,12 +127,12 @@ class CustomSamplerPt : public AbstractSampler {
   void Reset(bool initrandom = false) override {
     if (initrandom) {
       for (int i = 0; i < nrep_; i++) {
-        hilbert_.RandomVals(v_[i], this->GetRandomEngine());
+        GetHilbert().RandomVals(v_[i], this->GetRandomEngine());
       }
     }
 
     for (int i = 0; i < nrep_; i++) {
-      psi_.InitLookup(v_[i], lt_[i]);
+      GetMachine().InitLookup(v_[i], lt_[i]);
     }
 
     accept_ = Eigen::VectorXd::Zero(2 * nrep_);
@@ -161,7 +158,7 @@ class CustomSamplerPt : public AbstractSampler {
       }
 
       auto explo = std::exp(beta_[rep] *
-                            psi_.LogValDiff(v_[rep], tochange_[exit_state],
+                            GetMachine().LogValDiff(v_[rep], tochange_[exit_state],
                                             newconfs_[exit_state], lt_[rep]));
 
       double ratio = this->GetMachineFunc()(explo);
@@ -169,9 +166,9 @@ class CustomSamplerPt : public AbstractSampler {
       // Metropolis acceptance test
       if (ratio > distu(this->GetRandomEngine())) {
         accept_(rep) += 1;
-        psi_.UpdateLookup(v_[rep], tochange_[exit_state], newconfs_[exit_state],
+        GetMachine().UpdateLookup(v_[rep], tochange_[exit_state], newconfs_[exit_state],
                           lt_[rep]);
-        hilbert_.UpdateConf(v_[rep], tochange_[exit_state],
+        GetHilbert().UpdateConf(v_[rep], tochange_[exit_state],
                             newconfs_[exit_state]);
       }
       moves_(rep) += 1;
@@ -210,8 +207,8 @@ class CustomSamplerPt : public AbstractSampler {
 
   // computes the probability to exchange two replicas
   double ExchangeProb(int r1, int r2) {
-    const double lf1 = 2 * std::real(psi_.LogVal(v_[r1], lt_[r1]));
-    const double lf2 = 2 * std::real(psi_.LogVal(v_[r2], lt_[r2]));
+    const double lf1 = 2 * std::real(GetMachine().LogVal(v_[r1], lt_[r1]));
+    const double lf2 = 2 * std::real(GetMachine().LogVal(v_[r2], lt_[r2]));
 
     return std::exp((beta_[r1] - beta_[r2]) * (lf2 - lf1));
   }
@@ -225,14 +222,8 @@ class CustomSamplerPt : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd& v) override { v_[0] = v; }
 
-  AbstractMachine& GetMachine() noexcept override { return psi_; }
-
-  const AbstractHilbert& GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   AbstractMachine::VectorType DerLogVisible() override {
-    return psi_.DerLog(v_[0], lt_[0]);
+    return GetMachine().DerLog(v_[0], lt_[0]);
   }
 
   Eigen::VectorXd Acceptance() const override {
diff --git a/Sources/Sampler/exact_sampler.hpp b/Sources/Sampler/exact_sampler.hpp
index 365aca91..16bb00b7 100644
--- a/Sources/Sampler/exact_sampler.hpp
+++ b/Sources/Sampler/exact_sampler.hpp
@@ -27,10 +27,6 @@ namespace netket {
 
 // Exact sampling using heat bath, mostly for testing purposes on small systems
 class ExactSampler : public AbstractSampler {
-  AbstractMachine& psi_;
-
-  const AbstractHilbert& hilbert_;
-
   // number of visible units
   const int nv_;
 
@@ -43,7 +39,7 @@ class ExactSampler : public AbstractSampler {
   int mynode_;
   int totalnodes_;
 
-  const HilbertIndex &hilbert_index_;
+  const HilbertIndex& hilbert_index_;
 
   const int dim_;
 
@@ -54,10 +50,9 @@ class ExactSampler : public AbstractSampler {
 
  public:
   explicit ExactSampler(AbstractMachine& psi)
-      : psi_(psi),
-        hilbert_(psi.GetHilbert()),
-        nv_(hilbert_.Size()),
-        hilbert_index_(hilbert_.GetIndex()),
+      : AbstractSampler(psi),
+        nv_(GetHilbert().Size()),
+        hilbert_index_(GetHilbert().GetIndex()),
         dim_(hilbert_index_.NStates()) {
     Init();
   }
@@ -68,7 +63,7 @@ class ExactSampler : public AbstractSampler {
     MPI_Comm_size(MPI_COMM_WORLD, &totalnodes_);
     MPI_Comm_rank(MPI_COMM_WORLD, &mynode_);
 
-    if (!hilbert_.IsDiscrete()) {
+    if (!GetHilbert().IsDiscrete()) {
       throw InvalidInputError(
           "Exact sampler works only for discrete "
           "Hilbert spaces");
@@ -84,7 +79,7 @@ class ExactSampler : public AbstractSampler {
 
   void Reset(bool initrandom) override {
     if (initrandom) {
-      hilbert_.RandomVals(v_, this->GetRandomEngine());
+      GetHilbert().RandomVals(v_, this->GetRandomEngine());
     }
 
     double logmax = -std::numeric_limits<double>::infinity();
@@ -94,7 +89,7 @@ class ExactSampler : public AbstractSampler {
 
     for (int i = 0; i < dim_; ++i) {
       auto v = hilbert_index_.NumberToState(i);
-      logpsivals_[i] = psi_.LogVal(v);
+      logpsivals_[i] = GetMachine().LogVal(v);
       logmax = std::max(logmax, std::real(logpsivals_[i]));
     }
 
@@ -120,12 +115,6 @@ class ExactSampler : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd& v) override { v_ = v; }
 
-  AbstractMachine& GetMachine() noexcept override { return psi_; }
-
-  const AbstractHilbert& GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   Eigen::VectorXd Acceptance() const override {
     Eigen::VectorXd acc = accept_;
     for (int i = 0; i < 1; i++) {
diff --git a/Sources/Sampler/metropolis_exchange.hpp b/Sources/Sampler/metropolis_exchange.hpp
index fb1e0b26..2e419c41 100644
--- a/Sources/Sampler/metropolis_exchange.hpp
+++ b/Sources/Sampler/metropolis_exchange.hpp
@@ -26,10 +26,6 @@ namespace netket {
 
 // Metropolis sampling generating local exchanges
 class MetropolisExchange : public AbstractSampler {
-  AbstractMachine &psi_;
-
-  const AbstractHilbert &hilbert_;
-
   // number of visible units
   const int nv_;
 
@@ -51,7 +47,7 @@ class MetropolisExchange : public AbstractSampler {
  public:
   MetropolisExchange(const AbstractGraph &graph, AbstractMachine &psi,
                      int dmax = 1)
-      : psi_(psi), hilbert_(psi.GetHilbert()), nv_(hilbert_.Size()) {
+      : AbstractSampler(psi), nv_(GetHilbert().Size()) {
     Init(graph, dmax);
   }
 
@@ -92,11 +88,11 @@ class MetropolisExchange : public AbstractSampler {
   void Reset(bool initrandom = false) override {
     if (initrandom) {
       if (initrandom) {
-        hilbert_.RandomVals(v_, this->GetRandomEngine());
+        GetHilbert().RandomVals(v_, this->GetRandomEngine());
       }
     }
 
-    psi_.InitLookup(v_, lt_);
+    GetMachine().InitLookup(v_, lt_);
 
     accept_ = Eigen::VectorXd::Zero(1);
     moves_ = Eigen::VectorXd::Zero(1);
@@ -122,14 +118,14 @@ class MetropolisExchange : public AbstractSampler {
         newconf[0] = v_(sj);
         newconf[1] = v_(si);
 
-        auto explo = std::exp(psi_.LogValDiff(v_, tochange, newconf, lt_));
+        auto explo = std::exp(GetMachine().LogValDiff(v_, tochange, newconf, lt_));
 
         double ratio = this->GetMachineFunc()(explo);
 
         if (ratio > distu(this->GetRandomEngine())) {
           accept_[0] += 1;
-          psi_.UpdateLookup(v_, tochange, newconf, lt_);
-          hilbert_.UpdateConf(v_, tochange, newconf);
+          GetMachine().UpdateLookup(v_, tochange, newconf, lt_);
+          GetHilbert().UpdateConf(v_, tochange, newconf);
         }
       }
       moves_[0] += 1;
@@ -140,14 +136,8 @@ class MetropolisExchange : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd &v) override { v_ = v; }
 
-  AbstractMachine &GetMachine() noexcept override { return psi_; }
-
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   AbstractMachine::VectorType DerLogVisible() override {
-    return psi_.DerLog(v_, lt_);
+    return GetMachine().DerLog(v_, lt_);
   }
 
   Eigen::VectorXd Acceptance() const override {
diff --git a/Sources/Sampler/metropolis_exchange_pt.hpp b/Sources/Sampler/metropolis_exchange_pt.hpp
index bac2b3c8..0769e0c2 100644
--- a/Sources/Sampler/metropolis_exchange_pt.hpp
+++ b/Sources/Sampler/metropolis_exchange_pt.hpp
@@ -27,9 +27,6 @@ namespace netket {
 // Metropolis sampling generating local exchanges
 // Parallel tempering is also used
 class MetropolisExchangePt : public AbstractSampler {
-  AbstractMachine &psi_;
-  const AbstractHilbert &hilbert_;
-
   // number of visible units
   const int nv_;
 
@@ -55,9 +52,8 @@ class MetropolisExchangePt : public AbstractSampler {
  public:
   MetropolisExchangePt(const AbstractGraph &graph, AbstractMachine &psi,
                        int dmax = 1, int nreplicas = 1)
-      : psi_(psi),
-        hilbert_(psi.GetHilbert()),
-        nv_(hilbert_.Size()),
+      : AbstractSampler(psi),
+        nv_(GetHilbert().Size()),
         nrep_(nreplicas) {
     Init(graph, dmax);
   }
@@ -109,12 +105,12 @@ class MetropolisExchangePt : public AbstractSampler {
   void Reset(bool initrandom = false) override {
     if (initrandom) {
       for (int i = 0; i < nrep_; i++) {
-        hilbert_.RandomVals(v_[i], this->GetRandomEngine());
+        GetHilbert().RandomVals(v_[i], this->GetRandomEngine());
       }
     }
 
     for (int i = 0; i < nrep_; i++) {
-      psi_.InitLookup(v_[i], lt_[i]);
+      GetMachine().InitLookup(v_[i], lt_[i]);
     }
 
     accept_ = Eigen::VectorXd::Zero(2 * nrep_);
@@ -144,13 +140,13 @@ class MetropolisExchangePt : public AbstractSampler {
         newconf[1] = v_[rep](si);
 
         auto explo = std::exp(
-            beta_[rep] * psi_.LogValDiff(v_[rep], tochange, newconf, lt_[rep]));
+            beta_[rep] * GetMachine().LogValDiff(v_[rep], tochange, newconf, lt_[rep]));
         double ratio = this->GetMachineFunc()(explo);
 
         if (ratio > distu(this->GetRandomEngine())) {
           accept_(rep) += 1;
-          psi_.UpdateLookup(v_[rep], tochange, newconf, lt_[rep]);
-          hilbert_.UpdateConf(v_[rep], tochange, newconf);
+          GetMachine().UpdateLookup(v_[rep], tochange, newconf, lt_[rep]);
+          GetHilbert().UpdateConf(v_[rep], tochange, newconf);
         }
       }
 
@@ -190,8 +186,8 @@ class MetropolisExchangePt : public AbstractSampler {
 
   // computes the probability to exchange two replicas
   double ExchangeProb(int r1, int r2) {
-    const double lf1 = 2 * std::real(psi_.LogVal(v_[r1], lt_[r1]));
-    const double lf2 = 2 * std::real(psi_.LogVal(v_[r2], lt_[r2]));
+    const double lf1 = 2 * std::real(GetMachine().LogVal(v_[r1], lt_[r1]));
+    const double lf2 = 2 * std::real(GetMachine().LogVal(v_[r2], lt_[r2]));
 
     return std::exp((beta_[r1] - beta_[r2]) * (lf2 - lf1));
   }
@@ -205,14 +201,8 @@ class MetropolisExchangePt : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd &v) override { v_[0] = v; }
 
-  AbstractMachine &GetMachine() noexcept override { return psi_; }
-
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   AbstractMachine::VectorType DerLogVisible() override {
-    return psi_.DerLog(v_[0], lt_[0]);
+    return GetMachine().DerLog(v_[0], lt_[0]);
   }
 
   Eigen::VectorXd Acceptance() const override {
diff --git a/Sources/Sampler/metropolis_hamiltonian.hpp b/Sources/Sampler/metropolis_hamiltonian.hpp
index 9ec6b32a..ac9b860b 100644
--- a/Sources/Sampler/metropolis_hamiltonian.hpp
+++ b/Sources/Sampler/metropolis_hamiltonian.hpp
@@ -27,10 +27,6 @@ namespace netket {
 // Metropolis sampling generating transitions using the Hamiltonian
 template <class H>
 class MetropolisHamiltonian : public AbstractSampler {
-  AbstractMachine &psi_;
-
-  const AbstractHilbert &hilbert_;
-
   H &hamiltonian_;
 
   // number of visible units
@@ -60,10 +56,9 @@ class MetropolisHamiltonian : public AbstractSampler {
 
  public:
   MetropolisHamiltonian(AbstractMachine &psi, H &hamiltonian)
-      : psi_(psi),
-        hilbert_(psi.GetHilbert()),
+      : AbstractSampler(psi),
         hamiltonian_(hamiltonian),
-        nv_(hilbert_.Size()) {
+        nv_(GetHilbert().Size()) {
     Init();
   }
 
@@ -73,7 +68,7 @@ class MetropolisHamiltonian : public AbstractSampler {
     MPI_Comm_size(MPI_COMM_WORLD, &totalnodes_);
     MPI_Comm_rank(MPI_COMM_WORLD, &mynode_);
 
-    if (!hilbert_.IsDiscrete()) {
+    if (!GetHilbert().IsDiscrete()) {
       throw InvalidInputError(
           "Hamiltonian Metropolis sampler works only for discrete "
           "Hilbert spaces");
@@ -89,10 +84,10 @@ class MetropolisHamiltonian : public AbstractSampler {
 
   void Reset(bool initrandom = false) override {
     if (initrandom) {
-      hilbert_.RandomVals(v_, this->GetRandomEngine());
+      GetHilbert().RandomVals(v_, this->GetRandomEngine());
     }
 
-    psi_.InitLookup(v_, lt_);
+    GetMachine().InitLookup(v_, lt_);
 
     accept_ = Eigen::VectorXd::Zero(1);
     moves_ = Eigen::VectorXd::Zero(1);
@@ -112,21 +107,21 @@ class MetropolisHamiltonian : public AbstractSampler {
 
       // Inverse transition
       v1_ = v_;
-      hilbert_.UpdateConf(v1_, tochange_[si], newconfs_[si]);
+      GetHilbert().UpdateConf(v1_, tochange_[si], newconfs_[si]);
 
       hamiltonian_.FindConn(v1_, mel1_, tochange1_, newconfs1_);
 
       double w2 = tochange1_.size();
 
-      const auto lvd = psi_.LogValDiff(v_, tochange_[si], newconfs_[si], lt_);
+      const auto lvd = GetMachine().LogValDiff(v_, tochange_[si], newconfs_[si], lt_);
       double ratio = this->GetMachineFunc()(std::exp(lvd)) * w1 / w2;
 
 #ifndef NDEBUG
-      const auto psival1 = psi_.LogVal(v_);
-      if (std::abs(std::exp(psi_.LogVal(v_) - psi_.LogVal(v_, lt_)) - 1.) >
+      const auto psival1 = GetMachine().LogVal(v_);
+      if (std::abs(std::exp(GetMachine().LogVal(v_) - GetMachine().LogVal(v_, lt_)) - 1.) >
           1.0e-8) {
-        std::cerr << psi_.LogVal(v_) << "  and LogVal with Lt is "
-                  << psi_.LogVal(v_, lt_) << std::endl;
+        std::cerr << GetMachine().LogVal(v_) << "  and LogVal with Lt is "
+                  << GetMachine().LogVal(v_, lt_) << std::endl;
         std::abort();
       }
 #endif
@@ -134,16 +129,16 @@ class MetropolisHamiltonian : public AbstractSampler {
       // Metropolis acceptance test
       if (ratio > distu(this->GetRandomEngine())) {
         accept_[0] += 1;
-        psi_.UpdateLookup(v_, tochange_[si], newconfs_[si], lt_);
+        GetMachine().UpdateLookup(v_, tochange_[si], newconfs_[si], lt_);
         v_ = v1_;
 
 #ifndef NDEBUG
-        const auto psival2 = psi_.LogVal(v_);
+        const auto psival2 = GetMachine().LogVal(v_);
         if (std::abs(std::exp(psival2 - psival1 - lvd) - 1.) > 1.0e-8) {
           std::cerr << psival2 - psival1 << " and logvaldiff is " << lvd
                     << std::endl;
           std::cerr << psival2 << " and LogVal with Lt is "
-                    << psi_.LogVal(v_, lt_) << std::endl;
+                    << GetMachine().LogVal(v_, lt_) << std::endl;
           std::abort();
         }
 #endif
@@ -156,14 +151,8 @@ class MetropolisHamiltonian : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd &v) override { v_ = v; }
 
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
-  AbstractMachine &GetMachine() noexcept override { return psi_; }
-
   AbstractMachine::VectorType DerLogVisible() override {
-    return psi_.DerLog(v_, lt_);
+    return GetMachine().DerLog(v_, lt_);
   }
 
   Eigen::VectorXd Acceptance() const override {
diff --git a/Sources/Sampler/metropolis_hamiltonian_pt.hpp b/Sources/Sampler/metropolis_hamiltonian_pt.hpp
index 8caee2ca..9e402169 100644
--- a/Sources/Sampler/metropolis_hamiltonian_pt.hpp
+++ b/Sources/Sampler/metropolis_hamiltonian_pt.hpp
@@ -27,10 +27,6 @@ namespace netket {
 // Metropolis sampling generating transitions using the Hamiltonian
 template <class H>
 class MetropolisHamiltonianPt : public AbstractSampler {
-  AbstractMachine &psi_;
-
-  const AbstractHilbert &hilbert_;
-
   H &hamiltonian_;
 
   // number of visible units
@@ -62,10 +58,9 @@ class MetropolisHamiltonianPt : public AbstractSampler {
 
  public:
   MetropolisHamiltonianPt(AbstractMachine &psi, H &hamiltonian, int nrep)
-      : psi_(psi),
-        hilbert_(psi.GetHilbert()),
+      : AbstractSampler(psi),
         hamiltonian_(hamiltonian),
-        nv_(hilbert_.Size()),
+        nv_(GetHilbert().Size()),
         nrep_(nrep) {
     Init();
   }
@@ -74,7 +69,7 @@ class MetropolisHamiltonianPt : public AbstractSampler {
     MPI_Comm_size(MPI_COMM_WORLD, &totalnodes_);
     MPI_Comm_rank(MPI_COMM_WORLD, &mynode_);
 
-    if (!hilbert_.IsDiscrete()) {
+    if (!GetHilbert().IsDiscrete()) {
       throw InvalidInputError(
           "Hamiltonian Metropolis sampler works only for discrete "
           "Hilbert spaces");
@@ -105,12 +100,12 @@ class MetropolisHamiltonianPt : public AbstractSampler {
   void Reset(bool initrandom = false) override {
     if (initrandom) {
       for (int i = 0; i < nrep_; i++) {
-        hilbert_.RandomVals(v_[i], this->GetRandomEngine());
+        GetHilbert().RandomVals(v_[i], this->GetRandomEngine());
       }
     }
 
     for (int i = 0; i < nrep_; i++) {
-      psi_.InitLookup(v_[i], lt_[i]);
+      GetMachine().InitLookup(v_[i], lt_[i]);
     }
 
     accept_ = Eigen::VectorXd::Zero(2 * nrep_);
@@ -131,24 +126,24 @@ class MetropolisHamiltonianPt : public AbstractSampler {
 
       // Inverse transition
       v1_ = v_[rep];
-      hilbert_.UpdateConf(v1_, tochange_[si], newconfs_[si]);
+      GetHilbert().UpdateConf(v1_, tochange_[si], newconfs_[si]);
 
       hamiltonian_.FindConn(v1_, mel1_, tochange1_, newconfs1_);
 
       double w2 = tochange1_.size();
 
       const auto lvd =
-          psi_.LogValDiff(v_[rep], tochange_[si], newconfs_[si], lt_[rep]);
+          GetMachine().LogValDiff(v_[rep], tochange_[si], newconfs_[si], lt_[rep]);
       double ratio =
           this->GetMachineFunc()(std::exp(beta_[rep] * lvd)) * w1 / w2;
 
 #ifndef NDEBUG
-      const auto psival1 = psi_.LogVal(v_[rep]);
+      const auto psival1 = GetMachine().LogVal(v_[rep]);
       if (std::abs(
-              std::exp(psi_.LogVal(v_[rep]) - psi_.LogVal(v_[rep], lt_[rep])) -
+              std::exp(GetMachine().LogVal(v_[rep]) - GetMachine().LogVal(v_[rep], lt_[rep])) -
               1.) > 1.0e-8) {
-        std::cerr << psi_.LogVal(v_[rep]) << "  and LogVal with Lt is "
-                  << psi_.LogVal(v_[rep], lt_[rep]) << std::endl;
+        std::cerr << GetMachine().LogVal(v_[rep]) << "  and LogVal with Lt is "
+                  << GetMachine().LogVal(v_[rep], lt_[rep]) << std::endl;
         std::abort();
       }
 #endif
@@ -156,16 +151,16 @@ class MetropolisHamiltonianPt : public AbstractSampler {
       // Metropolis acceptance test
       if (ratio > distu(this->GetRandomEngine())) {
         accept_(rep) += 1;
-        psi_.UpdateLookup(v_[rep], tochange_[si], newconfs_[si], lt_[rep]);
+        GetMachine().UpdateLookup(v_[rep], tochange_[si], newconfs_[si], lt_[rep]);
         v_[rep] = v1_;
 
 #ifndef NDEBUG
-        const auto psival2 = psi_.LogVal(v_[rep]);
+        const auto psival2 = GetMachine().LogVal(v_[rep]);
         if (std::abs(std::exp(psival2 - psival1 - lvd) - 1.) > 1.0e-8) {
           std::cerr << psival2 - psival1 << " and logvaldiff is " << lvd
                     << std::endl;
           std::cerr << psival2 << " and LogVal with Lt is "
-                    << psi_.LogVal(v_[rep], lt_[rep]) << std::endl;
+                    << GetMachine().LogVal(v_[rep], lt_[rep]) << std::endl;
           std::abort();
         }
 #endif
@@ -206,8 +201,8 @@ class MetropolisHamiltonianPt : public AbstractSampler {
 
   // computes the probability to exchange two replicas
   double ExchangeProb(int r1, int r2) {
-    const double lf1 = 2 * std::real(psi_.LogVal(v_[r1], lt_[r1]));
-    const double lf2 = 2 * std::real(psi_.LogVal(v_[r2], lt_[r2]));
+    const double lf1 = 2 * std::real(GetMachine().LogVal(v_[r1], lt_[r1]));
+    const double lf2 = 2 * std::real(GetMachine().LogVal(v_[r2], lt_[r2]));
 
     return std::exp((beta_[r1] - beta_[r2]) * (lf2 - lf1));
   }
@@ -221,14 +216,8 @@ class MetropolisHamiltonianPt : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd &v) override { v_[0] = v; }
 
-  AbstractMachine &GetMachine() noexcept override { return psi_; }
-
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   AbstractMachine::VectorType DerLogVisible() override {
-    return psi_.DerLog(v_[0], lt_[0]);
+    return GetMachine().DerLog(v_[0], lt_[0]);
   }
 
   Eigen::VectorXd Acceptance() const override {
diff --git a/Sources/Sampler/metropolis_hop.hpp b/Sources/Sampler/metropolis_hop.hpp
index 4c6f1d55..31f2ec50 100644
--- a/Sources/Sampler/metropolis_hop.hpp
+++ b/Sources/Sampler/metropolis_hop.hpp
@@ -25,10 +25,6 @@ namespace netket {
 
 // Metropolis sampling generating local hoppings
 class MetropolisHop : public AbstractSampler {
-  AbstractMachine &psi_;
-
-  const AbstractHilbert &hilbert_;
-
   // number of visible units
   const int nv_;
 
@@ -52,8 +48,8 @@ class MetropolisHop : public AbstractSampler {
 
  public:
   MetropolisHop(AbstractMachine &psi, int dmax = 1)
-      : psi_(psi), hilbert_(psi.GetHilbert()), nv_(hilbert_.Size()) {
-    Init(hilbert_.GetGraph(), dmax);
+      : AbstractSampler(psi), nv_(GetHilbert().Size()) {
+    Init(GetHilbert().GetGraph(), dmax);
   }
 
   void Init(const AbstractGraph &graph, int dmax) {
@@ -65,8 +61,8 @@ class MetropolisHop : public AbstractSampler {
     accept_.resize(1);
     moves_.resize(1);
 
-    nstates_ = hilbert_.LocalSize();
-    localstates_ = hilbert_.LocalStates();
+    nstates_ = GetHilbert().LocalSize();
+    localstates_ = GetHilbert().LocalStates();
 
     GenerateClusters(graph, dmax);
 
@@ -95,10 +91,10 @@ class MetropolisHop : public AbstractSampler {
 
   void Reset(bool initrandom = false) override {
     if (initrandom) {
-      hilbert_.RandomVals(v_, this->GetRandomEngine());
+      GetHilbert().RandomVals(v_, this->GetRandomEngine());
     }
 
-    psi_.InitLookup(v_, lt_);
+    GetMachine().InitLookup(v_, lt_);
 
     accept_ = Eigen::VectorXd::Zero(1);
     moves_ = Eigen::VectorXd::Zero(1);
@@ -140,31 +136,31 @@ class MetropolisHop : public AbstractSampler {
         }
       }
 
-      const auto lvd = psi_.LogValDiff(v_, tochange, newconf, lt_);
+      const auto lvd = GetMachine().LogValDiff(v_, tochange, newconf, lt_);
       double ratio = this->GetMachineFunc()(std::exp(lvd));
 
 #ifndef NDEBUG
-      const auto psival1 = psi_.LogVal(v_);
-      if (std::abs(std::exp(psi_.LogVal(v_) - psi_.LogVal(v_, lt_)) - 1.) >
+      const auto psival1 = GetMachine().LogVal(v_);
+      if (std::abs(std::exp(GetMachine().LogVal(v_) - GetMachine().LogVal(v_, lt_)) - 1.) >
           1.0e-8) {
-        std::cerr << psi_.LogVal(v_) << "  and LogVal with Lt is "
-                  << psi_.LogVal(v_, lt_) << std::endl;
+        std::cerr << GetMachine().LogVal(v_) << "  and LogVal with Lt is "
+                  << GetMachine().LogVal(v_, lt_) << std::endl;
         std::abort();
       }
 #endif
 
       if (ratio > distu(this->GetRandomEngine())) {
         accept_[0] += 1;
-        psi_.UpdateLookup(v_, tochange, newconf, lt_);
-        hilbert_.UpdateConf(v_, tochange, newconf);
+        GetMachine().UpdateLookup(v_, tochange, newconf, lt_);
+        GetHilbert().UpdateConf(v_, tochange, newconf);
 
 #ifndef NDEBUG
-        const auto psival2 = psi_.LogVal(v_);
+        const auto psival2 = GetMachine().LogVal(v_);
         if (std::abs(std::exp(psival2 - psival1 - lvd) - 1.) > 1.0e-8) {
           std::cerr << psival2 - psival1 << " and logvaldiff is " << lvd
                     << std::endl;
           std::cerr << psival2 << " and LogVal with Lt is "
-                    << psi_.LogVal(v_, lt_) << std::endl;
+                    << GetMachine().LogVal(v_, lt_) << std::endl;
           std::abort();
         }
 #endif
@@ -177,14 +173,8 @@ class MetropolisHop : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd &v) override { v_ = v; }
 
-  AbstractMachine &GetMachine() noexcept override { return psi_; }
-
-  const AbstractHilbert &GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   AbstractMachine::VectorType DerLogVisible() override {
-    return psi_.DerLog(v_, lt_);
+    return GetMachine().DerLog(v_, lt_);
   }
 
   Eigen::VectorXd Acceptance() const override {
diff --git a/Sources/Sampler/metropolis_local.hpp b/Sources/Sampler/metropolis_local.hpp
index 0a74f716..0dcf9189 100644
--- a/Sources/Sampler/metropolis_local.hpp
+++ b/Sources/Sampler/metropolis_local.hpp
@@ -27,10 +27,6 @@ namespace netket {
 
 // Metropolis sampling generating local moves in hilbert space
 class MetropolisLocal : public AbstractSampler {
-  AbstractMachine& psi_;
-
-  const AbstractHilbert& hilbert_;
-
   // number of visible units
   const int nv_;
 
@@ -51,7 +47,7 @@ class MetropolisLocal : public AbstractSampler {
 
  public:
   explicit MetropolisLocal(AbstractMachine& psi)
-      : psi_(psi), hilbert_(psi.GetHilbert()), nv_(hilbert_.Size()) {
+      : AbstractSampler(psi), nv_(GetHilbert().Size()) {
     Init();
   }
 
@@ -61,7 +57,7 @@ class MetropolisLocal : public AbstractSampler {
     MPI_Comm_size(MPI_COMM_WORLD, &totalnodes_);
     MPI_Comm_rank(MPI_COMM_WORLD, &mynode_);
 
-    if (!hilbert_.IsDiscrete()) {
+    if (!GetHilbert().IsDiscrete()) {
       throw InvalidInputError(
           "Local Metropolis sampler works only for discrete "
           "Hilbert spaces");
@@ -70,8 +66,8 @@ class MetropolisLocal : public AbstractSampler {
     accept_.resize(1);
     moves_.resize(1);
 
-    nstates_ = hilbert_.LocalSize();
-    localstates_ = hilbert_.LocalStates();
+    nstates_ = GetHilbert().LocalSize();
+    localstates_ = GetHilbert().LocalStates();
 
     Reset(true);
 
@@ -80,10 +76,10 @@ class MetropolisLocal : public AbstractSampler {
 
   void Reset(bool initrandom) override {
     if (initrandom) {
-      hilbert_.RandomVals(v_, this->GetRandomEngine());
+      GetHilbert().RandomVals(v_, this->GetRandomEngine());
     }
 
-    psi_.InitLookup(v_, lt_);
+    GetMachine().InitLookup(v_, lt_);
 
     accept_ = Eigen::VectorXd::Zero(1);
     moves_ = Eigen::VectorXd::Zero(1);
@@ -114,15 +110,15 @@ class MetropolisLocal : public AbstractSampler {
         newconf[0] = localstates_[newstate];
       }
 
-      const auto lvd = psi_.LogValDiff(v_, tochange, newconf, lt_);
+      const auto lvd = GetMachine().LogValDiff(v_, tochange, newconf, lt_);
       double ratio = this->GetMachineFunc()(std::exp(lvd));
 
 #ifndef NDEBUG
-      const auto psival1 = psi_.LogVal(v_);
-      if (std::abs(std::exp(psi_.LogVal(v_) - psi_.LogVal(v_, lt_)) - 1.) >
+      const auto psival1 = GetMachine().LogVal(v_);
+      if (std::abs(std::exp(GetMachine().LogVal(v_) - GetMachine().LogVal(v_, lt_)) - 1.) >
           1.0e-8) {
-        std::cerr << psi_.LogVal(v_) << "  and LogVal with Lt is "
-                  << psi_.LogVal(v_, lt_) << std::endl;
+        std::cerr << GetMachine().LogVal(v_) << "  and LogVal with Lt is "
+                  << GetMachine().LogVal(v_, lt_) << std::endl;
         std::abort();
       }
 #endif
@@ -130,16 +126,16 @@ class MetropolisLocal : public AbstractSampler {
       // Metropolis acceptance test
       if (ratio > distu(this->GetRandomEngine())) {
         accept_[0] += 1;
-        psi_.UpdateLookup(v_, tochange, newconf, lt_);
-        hilbert_.UpdateConf(v_, tochange, newconf);
+        GetMachine().UpdateLookup(v_, tochange, newconf, lt_);
+        GetHilbert().UpdateConf(v_, tochange, newconf);
 
 #ifndef NDEBUG
-        const auto psival2 = psi_.LogVal(v_);
+        const auto psival2 = GetMachine().LogVal(v_);
         if (std::abs(std::exp(psival2 - psival1 - lvd) - 1.) > 1.0e-8) {
           std::cerr << psival2 - psival1 << " and logvaldiff is " << lvd
                     << std::endl;
           std::cerr << psival2 << " and LogVal with Lt is "
-                    << psi_.LogVal(v_, lt_) << std::endl;
+                    << GetMachine().LogVal(v_, lt_) << std::endl;
           std::abort();
         }
 #endif
@@ -152,14 +148,8 @@ class MetropolisLocal : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd& v) override { v_ = v; }
 
-  AbstractMachine& GetMachine() noexcept override { return psi_; }
-
-  const AbstractHilbert& GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   AbstractMachine::VectorType DerLogVisible() override {
-    return psi_.DerLog(v_, lt_);
+    return GetMachine().DerLog(v_, lt_);
   }
 
   Eigen::VectorXd Acceptance() const override {
diff --git a/Sources/Sampler/metropolis_local_pt.hpp b/Sources/Sampler/metropolis_local_pt.hpp
index 89f862eb..ff662ab6 100644
--- a/Sources/Sampler/metropolis_local_pt.hpp
+++ b/Sources/Sampler/metropolis_local_pt.hpp
@@ -27,10 +27,6 @@ namespace netket {
 // Metropolis sampling generating local changes
 // Parallel tempering is also used
 class MetropolisLocalPt : public AbstractSampler {
-  AbstractMachine& psi_;
-
-  const AbstractHilbert& hilbert_;
-
   // number of visible units
   const int nv_;
 
@@ -60,9 +56,8 @@ class MetropolisLocalPt : public AbstractSampler {
  public:
   // Constructor with one replica by default
   explicit MetropolisLocalPt(AbstractMachine& psi, int nreplicas = 1)
-      : psi_(psi),
-        hilbert_(psi.GetHilbert()),
-        nv_(hilbert_.Size()),
+      : AbstractSampler(psi),
+        nv_(GetHilbert().Size()),
         nrep_(nreplicas) {
     Init();
   }
@@ -71,8 +66,8 @@ class MetropolisLocalPt : public AbstractSampler {
     MPI_Comm_size(MPI_COMM_WORLD, &totalnodes_);
     MPI_Comm_rank(MPI_COMM_WORLD, &mynode_);
 
-    nstates_ = hilbert_.LocalSize();
-    localstates_ = hilbert_.LocalStates();
+    nstates_ = GetHilbert().LocalSize();
+    localstates_ = GetHilbert().LocalStates();
 
     SetNreplicas(nrep_);
 
@@ -103,12 +98,12 @@ class MetropolisLocalPt : public AbstractSampler {
   void Reset(bool initrandom = false) override {
     if (initrandom) {
       for (int i = 0; i < nrep_; i++) {
-        hilbert_.RandomVals(v_[i], this->GetRandomEngine());
+        GetHilbert().RandomVals(v_[i], this->GetRandomEngine());
       }
     }
 
     for (int i = 0; i < nrep_; i++) {
-      psi_.InitLookup(v_[i], lt_[i]);
+      GetMachine().InitLookup(v_[i], lt_[i]);
     }
 
     accept_ = Eigen::VectorXd::Zero(2 * nrep_);
@@ -141,16 +136,16 @@ class MetropolisLocalPt : public AbstractSampler {
         newconf[0] = localstates_[newstate];
       }
 
-      const auto lvd = psi_.LogValDiff(v_[rep], tochange, newconf, lt_[rep]);
+      const auto lvd = GetMachine().LogValDiff(v_[rep], tochange, newconf, lt_[rep]);
       double ratio = this->GetMachineFunc()(std::exp(beta_[rep] * lvd));
 
 #ifndef NDEBUG
-      const auto psival1 = psi_.LogVal(v_[rep]);
+      const auto psival1 = GetMachine().LogVal(v_[rep]);
       if (std::abs(
-              std::exp(psi_.LogVal(v_[rep]) - psi_.LogVal(v_[rep], lt_[rep])) -
+              std::exp(GetMachine().LogVal(v_[rep]) - GetMachine().LogVal(v_[rep], lt_[rep])) -
               1.) > 1.0e-8) {
-        std::cerr << psi_.LogVal(v_[rep]) << "  and LogVal with Lt is "
-                  << psi_.LogVal(v_[rep], lt_[rep]) << std::endl;
+        std::cerr << GetMachine().LogVal(v_[rep]) << "  and LogVal with Lt is "
+                  << GetMachine().LogVal(v_[rep], lt_[rep]) << std::endl;
         std::abort();
       }
 #endif
@@ -158,16 +153,16 @@ class MetropolisLocalPt : public AbstractSampler {
       if (ratio > distu(this->GetRandomEngine())) {
         accept_(rep) += 1;
 
-        psi_.UpdateLookup(v_[rep], tochange, newconf, lt_[rep]);
-        hilbert_.UpdateConf(v_[rep], tochange, newconf);
+        GetMachine().UpdateLookup(v_[rep], tochange, newconf, lt_[rep]);
+        GetHilbert().UpdateConf(v_[rep], tochange, newconf);
 
 #ifndef NDEBUG
-        const auto psival2 = psi_.LogVal(v_[rep]);
+        const auto psival2 = GetMachine().LogVal(v_[rep]);
         if (std::abs(std::exp(psival2 - psival1 - lvd) - 1.) > 1.0e-8) {
           std::cerr << psival2 - psival1 << " and logvaldiff is " << lvd
                     << std::endl;
           std::cerr << psival2 << " and LogVal with Lt is "
-                    << psi_.LogVal(v_[rep], lt_[rep]) << std::endl;
+                    << GetMachine().LogVal(v_[rep], lt_[rep]) << std::endl;
           std::abort();
         }
 #endif
@@ -208,8 +203,8 @@ class MetropolisLocalPt : public AbstractSampler {
 
   // computes the probability to exchange two replicas
   double ExchangeProb(int r1, int r2) {
-    const double lf1 = 2 * std::real(psi_.LogVal(v_[r1], lt_[r1]));
-    const double lf2 = 2 * std::real(psi_.LogVal(v_[r2], lt_[r2]));
+    const double lf1 = 2 * std::real(GetMachine().LogVal(v_[r1], lt_[r1]));
+    const double lf2 = 2 * std::real(GetMachine().LogVal(v_[r2], lt_[r2]));
 
     return std::exp((beta_[r1] - beta_[r2]) * (lf2 - lf1));
   }
@@ -223,14 +218,8 @@ class MetropolisLocalPt : public AbstractSampler {
 
   void SetVisible(const Eigen::VectorXd& v) override { v_[0] = v; }
 
-  AbstractMachine& GetMachine() noexcept override { return psi_; }
-
-  const AbstractHilbert& GetHilbert() const noexcept override {
-    return hilbert_;
-  }
-
   AbstractMachine::VectorType DerLogVisible() override {
-    return psi_.DerLog(v_[0], lt_[0]);
+    return GetMachine().DerLog(v_[0], lt_[0]);
   }
 
   Eigen::VectorXd Acceptance() const override {
diff --git a/Sources/Sampler/py_sampler.hpp b/Sources/Sampler/py_sampler.hpp
index 7467e915..32e74677 100644
--- a/Sources/Sampler/py_sampler.hpp
+++ b/Sources/Sampler/py_sampler.hpp
@@ -92,7 +92,8 @@ void AddSamplerModule(py::module &m) {
       .def_property_readonly("acceptance", &AbstractSampler::Acceptance, R"EOF(
         numpy.array: The measured acceptance rate for the sampling.
         In the case of rejection-free sampling this is always equal to 1.  )EOF")
-      .def_property_readonly("hilbert", &AbstractSampler::GetHilbert, R"EOF(
+      .def_property_readonly("hilbert", &AbstractSampler::GetHilbertShared,
+                             R"EOF(
         netket.hilbert: The Hilbert space used for the sampling.  )EOF")
       .def_property_readonly("machine", &AbstractSampler::GetMachine, R"EOF(
         netket.machine: The machine used for the sampling.  )EOF")
diff --git a/Sources/Unsupervised/quantum_state_reconstruction.hpp b/Sources/Unsupervised/quantum_state_reconstruction.hpp
index ba36da78..72e0f82f 100644
--- a/Sources/Unsupervised/quantum_state_reconstruction.hpp
+++ b/Sources/Unsupervised/quantum_state_reconstruction.hpp
@@ -39,7 +39,6 @@ class QuantumStateReconstruction {
 
   AbstractSampler &sampler_;
   AbstractMachine &psi_;
-  const AbstractHilbert &hilbert_;
   AbstractOptimizer &opt_;
   SR sr_;
   bool dosr_;
@@ -92,7 +91,6 @@ class QuantumStateReconstruction {
       bool use_cholesky = true)
       : sampler_(sampler),
         psi_(sampler_.GetMachine()),
-        hilbert_(psi_.GetHilbert()),
         opt_(opt),
         rotations_(rotations),
         trainingSamples_(trainingSamples),
diff --git a/Test/Machine/test_machine.py b/Test/Machine/test_machine.py
index 8d232629..324db917 100644
--- a/Test/Machine/test_machine.py
+++ b/Test/Machine/test_machine.py
@@ -12,7 +12,6 @@ def merge_dicts(x, y):
 
 machines = {}
 
-
 # TESTS FOR SPIN HILBERT
 # Constructing a 1d lattice
 g = nk.graph.Hypercube(length=4, n_dim=1)
diff --git a/Test/Operator/test_operator.py b/Test/Operator/test_operator.py
index 8d91ea51..e5dea524 100644
--- a/Test/Operator/test_operator.py
+++ b/Test/Operator/test_operator.py
@@ -130,3 +130,16 @@ def test_operator_is_hermitean():
                         foundinv = True
                         assert meli == np.conj(mel)
                 assert foundinv
+
+def test_no_segfault():
+    g = nk.graph.Hypercube(8, 1)
+    hi = nk.hilbert.Spin(g, 0.5)
+
+    lo = nk.operator.LocalOperator(hi, [[1,0],[0,1]], [0])
+    lo = lo.transpose()
+
+    hi = None
+
+    lo = lo * lo
+
+    assert True
['Sources/Sampler/metropolis_exchange.hpp', 'Sources/Machine/rbm_spin.cc', 'Sources/Operator/abstract_operator.hpp', 'Sources/Operator/graph_operator.hpp', 'Sources/Operator/local_operator.hpp', 'Sources/Sampler/metropolis_exchange_pt.hpp', 'Sources/Machine/py_machine.cc', 'Sources/Machine/DensityMatrices/abstract_density_matrix.hpp', 'Sources/Sampler/metropolis_local_pt.hpp', 'Sources/Machine/DensityMatrices/ndm_spin_phase.hpp', 'Sources/Machine/rbm_spin.hpp', 'Sources/Machine/rbm_spin_symm.cc', 'Sources/Machine/DensityMatrices/diagonal_density_matrix.hpp', 'Sources/Operator/py_bosonhubbard.hpp', 'Sources/Machine/jastrow_symm.cc', 'Sources/Machine/jastrow.hpp', 'Sources/Operator/MatrixWrapper/dense_matrix_wrapper.hpp', 'Sources/Sampler/metropolis_local.hpp', 'Sources/Machine/rbm_spin_real.hpp', 'Sources/Hilbert/py_hilbert.cc', 'Sources/Sampler/metropolis_hamiltonian_pt.hpp', 'Test/Machine/test_machine.py', 'Sources/Sampler/custom_sampler.hpp', 'Sources/Machine/ffnn.hpp', 'Sources/Operator/py_local_operator.hpp', 'Sources/Sampler/metropolis_hop.hpp', 'Sources/Machine/rbm_spin_phase.cc', 'Sources/Sampler/abstract_sampler.hpp', 'Sources/Machine/mps_periodic.hpp', 'Sources/Machine/jastrow_symm.hpp', 'Sources/Sampler/exact_sampler.hpp', 'Sources/Machine/DensityMatrices/py_ndm_spin_phase.hpp', 'Sources/Machine/abstract_machine.cc', 'Sources/Sampler/metropolis_hamiltonian.hpp', 'Test/Operator/test_operator.py', 'Sources/Operator/py_graph_operator.hpp', 'Sources/Machine/jastrow.cc', 'Sources/Operator/MatrixWrapper/sparse_matrix_wrapper.hpp', 'Sources/Machine/rbm_spin_real.cc', 'Sources/Unsupervised/quantum_state_reconstruction.hpp', 'Sources/Machine/rbm_spin_phase.hpp', 'Sources/Machine/rbm_spin_symm.hpp', 'Sources/Operator/py_operator.hpp', 'Sources/Machine/rbm_multival.hpp', 'Sources/Sampler/custom_sampler_pt.hpp', 'Sources/Operator/bosonhubbard.hpp', 'Sources/Machine/mps_periodic.cc', 'Sources/Sampler/py_sampler.hpp', 'Sources/Machine/rbm_multival.cc', 'Sources/Machine/abstract_machine.hpp']