diff --git a/.gitignore b/.gitignore
index 038ef9d..7d7338b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,9 @@
 # Files
 .DS_Store
+.env*
+.venv*
+.vscode
+*.code-workspace
 
 .ipynb_checkpoints/
 *.ipynb
diff --git a/dppy/dcftp.py b/dppy/dcftp.py
new file mode 100644
index 0000000..147f528
--- /dev/null
+++ b/dppy/dcftp.py
@@ -0,0 +1,77 @@
+from math import sqrt, log
+from numpy.random import poisson
+from random import random, sample
+from collections import deque
+# With the current Python implementation, we use both random and numpy.random and we do not populate the RNG, hence reproducibility is not possible.
+
+
+def distance(x, y):
+    return sqrt((y[0] - x[0])**2 + (y[1] - x[1])**2)
+
+def backward_update(D, X, M, n, beta, w, h):
+
+    b_rate = beta * w * h
+
+    for _ in range(-n, -2 * n, -1):
+
+        card_D = len(D)
+
+        if random() < card_D / (b_rate + card_D):
+            x = sample(D, 1)[0]
+            D.discard(x)
+            X.appendleft(x)
+            M.appendleft(random())
+        else:
+            x = (w * random(), h * random())
+            D.add(x)
+            X.appendleft(x)
+            M.appendleft(0.0)
+
+    return D, X, M
+
+def forward_coupling(D, X, M, gamma, r):
+
+    L, U = set(), D.copy()
+    log_g = log(gamma)
+
+    for x, m in zip(X, M):
+        if m:
+            log_m = log(m)
+
+            neigh_x_U = sum(distance(x, u) <= r for u in U)
+            if log_m < neigh_x_U * log_g:
+                L.add(x)
+                U.add(x)
+            else:
+                neigh_x_L = sum(distance(x, l) <= r for l in L)
+                if log_m < neigh_x_L * log_g:
+                    U.add(x)
+
+        else:
+            L.discard(x)
+            U.discard(x)
+
+    return L, U
+
+def strauss_dcftp(beta, gamma, r, w, h):
+
+    lmbda = poisson(beta * w * h)
+    D = set((w * random(), h * random()) for _ in range(lmbda))
+    X = deque()
+    M = deque()
+
+    n = 1
+    while True:
+        D, X, M = backward_update(D, X, M, n, beta, w, h)
+        L, U = forward_coupling(D, X, M, gamma, r)
+        if len(L) == len(U):  # coalescence occurs
+            return L, n
+        n *= 2
+
+
+if __name__ == '__main__':
+
+    beta, gamma, r = 2, 0.2, 0.7
+    w, h = 10, 10
+    config, iter_ = strauss_dcftp(beta, gamma, r, w, h)
+    print(len(config), iter_)
diff --git a/dppy/exotic_dpps.py b/dppy/exotic_dpps.py
index c0ab69e..94cfc30 100644
--- a/dppy/exotic_dpps.py
+++ b/dppy/exotic_dpps.py
@@ -464,7 +464,6 @@ class UST:
     """
 
     def __init__(self, graph):
-        # For Uniform Spanning Trees
         try:
             import networkx as nx
             self.nx = nx
@@ -474,19 +473,7 @@ def __init__(self, graph):
         if nx.is_connected(graph):
             self.graph = graph
         else:
-            raise ValueError('graph not connected')
-
-        self.nodes = list(self.graph.nodes())
-        self.nb_nodes = self.graph.number_of_nodes()  # len(self.graph)
-
-        self.edges = list(self.graph.edges())
-        self.nb_edges = self.graph.number_of_edges()  # len(self.edges)
-
-        self.edge_labels = {edge: r'$e_{}$'.format(i)
-                            for i, edge in enumerate(self.edges)}
-
-        self.neighbors = [list(graph.neighbors(v))
-                          for v in range(self.nb_nodes)]
+            raise ValueError('graph is not connected')
 
         self.sampling_mode = 'Wilson'  # Default (avoid eig_vecs computation)
         self._sampling_modes = {'markov-chain': ['Wilson', 'Aldous-Broder'],
@@ -500,29 +487,23 @@ def __init__(self, graph):
     def __str__(self):
 
         str_info = ['Uniform Spanning Tree measure on a graph with:',
-                    '- {} nodes'.format(self.nb_nodes),
-                    '- {} edges'.format(self.nb_edges),
+                    '- {} nodes'.format(self.graph.number_of_nodes()),
+                    '- {} edges'.format(self.graph.number_of_edges()),
                     'Sampling mode = {}'.format(self.sampling_mode),
                     'Number of samples = {}'.format(len(self.list_of_samples))]
 
         return '\n'.join(str_info)
 
-    # def info(self):
-    #     """ Print infos about the :class:`UST` object
-    #     """
-    #     print(self.__str__())
-
     def flush_samples(self):
         """ Empty the :py:attr:`list_of_samples` attribute.
         """
         self.list_of_samples = []
 
-    def sample(self, mode='Wilson', root=None, random_state=None):
+    def sample(self, mode='Wilson', random_state=None):
         """ Sample a spanning of the underlying graph uniformly at random.
         It generates a networkx graph object.
 
         :param mode:
-
             Markov-chain-based samplers:
 
             - ``'Wilson'``, ``'Aldous-Broder'``
@@ -531,15 +512,9 @@ def sample(self, mode='Wilson', root=None, random_state=None):
 
             - ``'GS'``, ``'GS_bis'``, ``'KuTa12'`` from eigenvectors
             - ``'Schur'``, ``'Chol'``, from :math:`\\mathbf{K}` correlation kernel
-
         :type mode:
             string, default ``'Wilson'``
 
-        :param root:
-            Starting node of the random walk when using Markov-chain-based samplers
-        :type root:
-            int
-
         :param random_state:
         :type random_state:
             None, np.random, int, np.random.RandomState
@@ -554,15 +529,14 @@ def sample(self, mode='Wilson', root=None, random_state=None):
         rng = check_random_state(random_state)
 
         self.sampling_mode = mode
+        sampl = self.nx.Graph()
 
         if self.sampling_mode in self._sampling_modes['markov-chain']:
             if self.sampling_mode == 'Wilson':
-                sampl = ust_sampler_wilson(self.neighbors,
-                                           random_state=rng)
+                sampl = ust_sampler_wilson(self.graph, random_state=rng)
 
             elif self.sampling_mode == 'Aldous-Broder':
-                sampl = ust_sampler_aldous_broder(self.neighbors,
-                                                  random_state=rng)
+                sampl = ust_sampler_aldous_broder(self.graph, random_state=rng)
 
         elif self.sampling_mode in self._sampling_modes['spectral-method']:
 
@@ -570,9 +544,8 @@ def sample(self, mode='Wilson', root=None, random_state=None):
             dpp_sample = proj_dpp_sampler_eig(self.kernel_eig_vecs,
                                               mode=self.sampling_mode,
                                               random_state=rng)
-
-            sampl = self.nx.Graph()
-            sampl.add_edges_from([self.edges[e] for e in dpp_sample])
+            edges = list(self.graph.edges)
+            sampl = self.nx.from_edgelist([edges[i] for i in dpp_sample])
 
         elif self.sampling_mode in self._sampling_modes['projection-K-kernel']:
 
@@ -581,8 +554,8 @@ def sample(self, mode='Wilson', root=None, random_state=None):
                                                  mode=self.sampling_mode,
                                                  random_state=rng)
 
-            sampl = self.nx.Graph()
-            sampl.add_edges_from([self.edges[e] for e in dpp_sample])
+            edges = list(self.graph.edges)
+            sampl = self.nx.from_edgelist([edges[i] for i in dpp_sample])
 
         else:
             err_print = '\n'.join(
@@ -592,6 +565,7 @@ def sample(self, mode='Wilson', root=None, random_state=None):
             raise ValueError()
 
         self.list_of_samples.append(sampl)
+        return sampl
 
     def compute_kernel(self):
         """ Compute the orthogonal projection kernel :math:`\\mathbf{K} = \\text{Inc}^+ \\text{Inc}` i.e. onto the span of the rows of the vertex-edge incidence matrix :math:`\\text{Inc}` of size :math:`|V| \\times |E|`.
@@ -603,25 +577,27 @@ def compute_kernel(self):
         .. seealso::
 
             - :py:meth:`plot_kernel`
+            - :py:meth:`compute_kernel_eig_vecs`
         """
-
+        # K = UU.T, with U = QR(Inc[:-1,:].T)
         if self.kernel is None:
-            self.compute_kernel_eig_vecs()  # U = QR(Inc[:-1,:].T)
-            # K = UU.T
-            self.kernel = self.kernel_eig_vecs.dot(self.kernel_eig_vecs.T)
+            U = self.compute_kernel_eig_vecs()
+            self.kernel = U.dot(U.T)
+
+        return self.kernel
 
     def compute_kernel_eig_vecs(self):
         """ See explaination in :func:`compute_kernel <compute_kernel>`
         """
-
         if self.kernel_eig_vecs is None:
-
             inc_mat = self.nx.incidence_matrix(self.graph, oriented=True)
             # Discard any row e.g. the last one
             A = inc_mat[:-1, :].toarray()
             # Orthonormalize rows of A
             self.kernel_eig_vecs, _ = qr(A.T, mode='economic')
 
+        return self.kernel_eig_vecs
+
     def plot(self, title=''):
         """ Display the last realization (spanning tree) of the corresponding :class:`UST` object.
 
@@ -642,17 +618,19 @@ def plot(self, title=''):
 
         pos = self.nx.circular_layout(self.graph)
         self.nx.draw_networkx(graph_to_plot,
-                         pos=pos,
-                         node_color='orange',
-                         with_labels=True,
-                         width=3)
-
-        edge_labs = {e: self.edge_labels[e if e in self.edges else e[::-1]]
-                     for e in graph_to_plot.edges()}
+                              pos=pos,
+                              node_color='orange',
+                              with_labels=True,
+                              width=3)
+
+        labels = {e: r'$e_{}$'.format(i)
+                  for i, e in enumerate(self.graph.edges)}
+        edge_labels = {e: labels[e if e in labels else e[::-1]]
+                       for e in graph_to_plot.edges}
         self.nx.draw_networkx_edge_labels(graph_to_plot,
-                                     pos=pos,
-                                     edge_labels=edge_labs,
-                                     font_size=20)
+                                          pos=pos,
+                                          edge_labels=edge_labels,
+                                          font_size=20)
 
         plt.axis('off')
 
@@ -673,25 +651,23 @@ def plot_graph(self, title=''):
             - :func:`compute_kernel <compute_kernel>`
         """
 
-        # edge_lab = [r'$e_{}$'.format(i) for i in range(self.nb_edges)]
-        # edge_labels = dict(zip(self.edges, edge_lab))
-        # node_labels = dict(zip(self.nodes, self.nodes))
-
         plt.figure(figsize=(4, 4))
 
         pos = self.nx.circular_layout(self.graph)
         self.nx.draw_networkx(self.graph,
-                         pos=pos,
-                         node_color='orange',
-                         with_labels=True,
-                         width=3)
+                              pos=pos,
+                              node_color='orange',
+                              with_labels=True,
+                              width=3)
         # nx.draw_networkx_labels(self.graph,
         #                         pos,
         #                         node_labels)
+        edge_labels = {e: r'$e_{}$'.format(i)
+                       for i, e in enumerate(self.graph.edges)}
         self.nx.draw_networkx_edge_labels(self.graph,
-                                     pos=pos,
-                                     edge_labels=self.edge_labels,
-                                     font_size=20)
+                                          pos=pos,
+                                          edge_labels=edge_labels,
+                                          font_size=20)
 
         plt.axis('off')
 
@@ -720,7 +696,7 @@ def plot_kernel(self, title=''):
 
         ax.set_aspect('equal')
 
-        ticks = np.arange(self.nb_edges)
+        ticks = np.arange(self.graph.number_of_edges())
         ticks_label = [r'${}$'.format(tic) for tic in ticks]
 
         ax.xaxis.tick_top()
diff --git a/dppy/exotic_dpps_core.py b/dppy/exotic_dpps_core.py
index ee0f80f..a77f2fa 100644
--- a/dppy/exotic_dpps_core.py
+++ b/dppy/exotic_dpps_core.py
@@ -22,20 +22,25 @@
     `Documentation on ReadTheDocs <https://dppy.readthedocs.io/en/latest/exotic_dpps/index.html>`_
 """
 
-import functools  # used for decorators to pass docstring
-
 import numpy as np
-
-from itertools import chain  # create graph edges from path
-
-# For class PoissonizedPlancherel
 from bisect import bisect_right  # for RSK
 
 from dppy.utils import check_random_state
 
 
-def ust_sampler_wilson(list_of_neighbors, root=None,
-                       random_state=None):
+def ust_sampler_wilson(g, root=None, random_state=None):
+    """Generate a uniform spanning tree of g at root using [Wilson's algorithm](https://dl.acm.org/doi/10.1145/237814.237880).
+
+    :param g: Connected graph
+    :type g: nx.Graph
+
+    :param root: Any node of g, defaults to None. If None, the root is chosen uniformly at random among g.nodes.
+    :type root: list, optional
+
+    :return: uniform spanning tree of g
+    :rtype: nx.Graph
+    """
+
     try:
         import networkx as nx
     except ImportError:
@@ -43,100 +48,71 @@ def ust_sampler_wilson(list_of_neighbors, root=None,
 
     rng = check_random_state(random_state)
 
-    # Initialize the tree
-    wilson_tree_graph = nx.Graph()
-    nb_nodes = len(list_of_neighbors)
-
-    # Initialize the root, if root not specified start from any node
-    n0 = root if root else rng.choice(nb_nodes)  # size=1)[0]
-    # -1 = not visited / 0 = in path / 1 = in tree
-    state = -np.ones(nb_nodes, dtype=int)
-    state[n0] = 1
-    nb_nodes_in_tree = 1
-
-    path, branches = [], []  # branches of tree, temporary path
-
-    while nb_nodes_in_tree < nb_nodes:  # |Tree| = |V| - 1
+    if root is None:
+        nodes = list(g.nodes)
+        root = nodes[rng.randint(len(nodes))]
+    elif root not in g:
+        raise ValueError("root not in g.nodes")
 
-        # visit a neighbor of n0 uniformly at random
-        n1 = rng.choice(list_of_neighbors[n0])  # size=1)[0]
+    tree = {root}
+    successor = dict.fromkeys(g.nodes, None)
+    del successor[root]
 
-        if state[n1] == -1:  # not visited => continue the walk
+    for i in g.nodes:
+        # Run a natural random walk from i until it hits a node in tree
+        u = i
+        while u not in tree:
+            neighbors = list(g.neighbors(u))
+            successor[u] = neighbors[rng.randint(len(neighbors))]
+            u = successor[u]
 
-            path.append(n1)  # add it to the path
-            state[n1] = 0  # mark it as in the path
-            n0 = n1  # continue the walk
+        # Record Erase first loop created during the random walk
+        u = i
+        while u not in tree:
+            tree.add(u)
+            u = successor[u]
 
-        if state[n1] == 0:  # loop on the path => erase the loop
+    return nx.from_edgelist(successor.items())
 
-            knot = path.index(n1)  # find 1st appearence of n1 in the path
-            nodes_loop = path[knot + 1:]  # identify nodes forming the loop
-            del path[knot + 1:]  # erase the loop
-            state[nodes_loop] = -1  # mark loopy nodes as not visited
-            n0 = n1  # continue the walk
-
-        elif state[n1] == 1:  # hits the tree => new branch
-
-            if nb_nodes_in_tree == 1:
-                branches.append([n1] + path)  # initial branch of the tree
-            else:
-                branches.append(path + [n1])  # path as a new branch
 
-            state[path] = 1  # mark nodes in path as in the tree
-            nb_nodes_in_tree += len(path)
+def ust_sampler_aldous_broder(g, root=None, random_state=None):
+    """Generate a uniform spanning tree of g at root using [Aldous](https://epubs.siam.org/doi/10.1137/0403039)-[Broder](https://doi.org/10.1109/SFCS.1989.63516)'s algorithm
 
-            # Restart the walk from a random node among those not visited
-            nodes_not_visited = np.where(state == -1)[0]
-            if nodes_not_visited.size:
-                n0 = rng.choice(nodes_not_visited)  # size=1)[0]
-                path = [n0]
+    :param g: Connected graph
+    :type g: nx.Graph
 
-    tree_edges = list(chain.from_iterable(map(lambda x: zip(x[:-1], x[1:]),
-                                              branches)))
-    wilson_tree_graph.add_edges_from(tree_edges)
-
-    return wilson_tree_graph
-
-
-def ust_sampler_aldous_broder(list_of_neighbors, root=None,
-                              random_state=None):
+    :param root: Any node of g, defaults to None. If None, the root is chosen uniformly at random among g.nodes.
+    :type root: list, optional
 
+    :return: uniform spanning tree of g
+    :rtype: nx.Graph
+    """
     try:
         import networkx as nx
     except ImportError:
         raise ValueError('The networkx package is required to sample spanning trees (see setup.py).')
-
     rng = check_random_state(random_state)
 
-    # Initialize the tree
-    aldous_tree_graph = nx.Graph()
-    nb_nodes = len(list_of_neighbors)
-
-    # Initialize the root, if root not specified start from any node
-    n0 = root if root else rng.choice(nb_nodes)  # size=1)[0]
-    visited = np.zeros(nb_nodes, dtype=bool)
-    visited[n0] = True
-    nb_nodes_in_tree = 1
-
-    tree_edges = np.zeros((nb_nodes - 1, 2), dtype=np.int)
-
-    while nb_nodes_in_tree < nb_nodes:
-
-        # visit a neighbor of n0 uniformly at random
-        n1 = rng.choice(list_of_neighbors[n0])  # size=1)[0]
-
-        if visited[n1]:
-            pass  # continue the walk
-        else:  # create edge (n0, n1) and continue the walk
-            tree_edges[nb_nodes_in_tree - 1] = [n0, n1]
-            visited[n1] = True  # mark it as in the tree
-            nb_nodes_in_tree += 1
-
-        n0 = n1
-
-    aldous_tree_graph.add_edges_from(tree_edges)
-
-    return aldous_tree_graph
+    if root is None:
+        nodes = list(g.nodes)
+        root = nodes[rng.randint(len(nodes))]
+    elif root not in g:
+        raise ValueError("root not in g.nodes")
+
+    # Run a natural random walk from root a until all nodes are visited
+    tree = {root: None}
+    u = root
+    while len(tree) < g.number_of_nodes():
+        neighbors = list(g.neighbors(u))
+        v = neighbors[rng.randint(len(neighbors))]
+        # Record an edge when reaching an unvisited node
+        if v not in tree:
+            tree[v] = u
+        u = v
+
+    del tree[root]
+
+    return nx.from_edgelist(tree.items())
 
 
 def uniform_permutation(N, random_state=None):
diff --git a/requirements.txt b/requirements.txt
index 693ed64..9fa6850 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -7,7 +7,7 @@ matplotlib
 # for the zonotope MCMC sampler
 cvxopt==1.2.1
 # for notebook on sampling spanning trees
-networkx 
+networkx
 # for the documentation
 sphinxcontrib-bibtex
-sphinx_rtd_theme
\ No newline at end of file
+sphinx_rtd_theme