diff --git a/book/mle/compare_smc_sqmc_malik_and_pitt.py b/book/mle/compare_smc_sqmc_malik_and_pitt.py
index 8e15e49..8bf4707 100644
--- a/book/mle/compare_smc_sqmc_malik_and_pitt.py
+++ b/book/mle/compare_smc_sqmc_malik_and_pitt.py
@@ -18,32 +18,49 @@
 from __future__ import division, print_function
 
 import numpy as np
-from numpy import random
+from malikpitt_interpolation import MalikPitt_SMC
 from matplotlib import pyplot as plt
-import matplotlib.cm as cm
-from scipy.optimize import fmin 
+from numpy import random
 
 import particles
+from particles import ensemble_transform as et
 from particles import state_space_models as ssms
 
-from malikpitt_interpolation import MalikPitt_SMC
-
 # data 
 raw_data = np.loadtxt('../../datasets/GBP_vs_USD_9798.txt',
                       skiprows=2, usecols=(3,), comments='(C)')
 T = 200
-data = 100. * np.diff(np.log(raw_data[:(T+1)])) 
+data = 100. * np.diff(np.log(raw_data[:(T + 1)]))
+
 
 def fkmod(theta):
-    mu = theta[0]; rho = theta[1]; sigma = theta[2] 
-    return ssms.Bootstrap(ssm=ssms.StochVol(mu=mu, rho=rho, sigma=sigma), 
-                         data=data)
+    mu = theta[0];
+    rho = theta[1];
+    sigma = theta[2]
+    return ssms.Bootstrap(ssm=ssms.StochVol(mu=mu, rho=rho, sigma=sigma),
+                          data=data)
+
 
-def loglik(theta, seed=None, qmc=False, N=10**4, verbose=False, interpol=False):
+def loglik(theta, seed=None, qmc=False, N=10 ** 4, verbose=False, interpol=False, transform=False, epsilon=0.,
+           corrected=False):
     if seed is not None:
         random.seed(seed)
     if interpol:
-        alg = MalikPitt_SMC(fk=fkmod(theta), N=N) 
+        alg = MalikPitt_SMC(fk=fkmod(theta), N=N)
+    elif transform:
+        if epsilon:
+            if corrected:
+                et_instance = et.CorrectedEntropyRegularizedEnsembleTransform(epsilon)
+            else:
+                et_instance = et.EntropyRegularizedEnsembleTransform(epsilon)
+
+        else:
+            et_instance = et.EnsembleTransform()
+        alg = et.EnsembleTransformFilter(et_instance,
+                                         fk=fkmod(theta),
+                                         N=N,
+                                         qmc=qmc)
+
     else:
         alg = particles.SMC(fk=fkmod(theta), N=N, qmc=qmc)
     alg.run()
@@ -52,37 +69,64 @@ def loglik(theta, seed=None, qmc=False, N=10**4, verbose=False, interpol=False):
         print(theta, out)
     return out
 
-#theta0 = [-1., 0.9, 0.3] 
-#mle = fmin(objfunc, theta0, xtol=.1, ftol=.1, maxfun=1000)
-
-sig_min = 0.2; sig_max = 0.5
-mu = -1.; rho = 0.9 
-sigmas = np.linspace(sig_min, sig_max, 1000)
-
-# SMC 
-ll = [loglik([mu, rho, sigma]) for sigma in sigmas]
-# SMC frozen seed
-ll_frozen = [loglik([mu, rho, sigma], seed=4) for sigma in sigmas]
-# interpolation
-ll_pol = [loglik([mu, rho, sigma], seed=4, interpol=True) for sigma in sigmas]
-# QMC 
-ll_qmc  = [loglik([mu, rho, sigma], qmc=True) for sigma in sigmas]
-
-# PLOT
-# ====
-savefigs = True  # False if you don't want to save plots as pdfs
-
-plt.figure()
-plt.style.use('seaborn-dark')
-plt.grid('on') 
-plt.plot(sigmas, ll, '+', color='gray', label='standard')
-plt.plot(sigmas, ll_frozen, 'o', color='gray', label='fixed seed')
-plt.plot(sigmas, ll_pol, color='gray', label='interpolation')
-plt.plot(sigmas, ll_qmc,'k.', lw=2, label='SQMC')
-plt.xlabel('sigma')
-plt.ylabel('log-likelihood')
-plt.legend()
-if savefigs:
-    plt.savefig('loglik_interpolated_frozen_sqmc.pdf') 
-
-plt.show()
+
+if __name__ == "__main__":
+    # theta0 = [-1., 0.9, 0.3]
+    # mle = fmin(objfunc, theta0, xtol=.1, ftol=.1, maxfun=1000)
+
+    sig_min = 0.2
+    sig_max = 0.5
+    mu = -1.
+    rho = 0.9
+    sigmas = np.linspace(sig_min, sig_max, 20)
+    print("SMC")
+    # SMC
+    ll = [loglik([mu, rho, sigma]) for sigma in sigmas]
+    # SMC frozen seed
+    print("SMC frozen")
+    ll_frozen = [loglik([mu, rho, sigma], seed=4) for sigma in sigmas]
+    # interpolation
+    print("SMC interpolated")
+    ll_pol = [loglik([mu, rho, sigma], seed=4, interpol=True) for sigma in sigmas]
+    print("ET")
+    # ensemble transform
+    ll_et = [loglik([mu, rho, sigma], seed=4, transform=True) for sigma in sigmas]
+    # entropy ensemble transform
+    print("Reg-ET")
+    ll_reg_et = [loglik([mu, rho, sigma], seed=4, transform=True, epsilon=0.1) for sigma in sigmas]
+    print("Corrected-ET")
+    ll_reg_corr = [loglik([mu, rho, sigma], seed=4, transform=True, epsilon=0.1, corrected=True) for sigma in sigmas]
+
+    # QMC
+    ll_qmc = [loglik([mu, rho, sigma], qmc=True) for sigma in sigmas]
+    # ET QMC
+    ll_et_qmc = [loglik([mu, rho, sigma], qmc=True, transform=True) for sigma in sigmas]
+    # REG ET QMC
+    ll_reg_et_qmc = [loglik([mu, rho, sigma], qmc=True, transform=True, epsilon=0.1) for sigma in sigmas]
+    # ET QMC
+    ll_corr_et_qmc = [loglik([mu, rho, sigma], qmc=True, transform=True, epsilon=0.1, corrected=True) for sigma in sigmas]
+
+    # PLOT
+    # ====
+    savefigs = True  #  False if you don't want to save plots as pdfs
+
+    plt.figure()
+    plt.style.use('seaborn-dark')
+    plt.grid('on')
+    plt.plot(sigmas, ll, '+', color='gray', label='standard')
+    plt.plot(sigmas, ll_frozen, 'o', color='gray', label='fixed seed')
+    plt.plot(sigmas, ll_pol, color='gray', label='interpolation')
+    plt.plot(sigmas, ll_qmc, 'k.', lw=2, label='SQMC')
+    plt.plot(sigmas, ll_et, 'k.', lw=2, label='ET')
+    plt.plot(sigmas, ll_reg_et, 'k*', lw=2, label='Reg-ET')
+    plt.plot(sigmas, ll_reg_corr, 'k--', lw=2, label='Corrected-ET')
+    plt.plot(sigmas, ll_et_qmc, 'ko', lw=2, label='ET-QMC')
+    plt.plot(sigmas, ll_reg_et_qmc, 'k.-', lw=2, label='REG-QMC')
+    plt.plot(sigmas, ll_corr_et_qmc, 'k..', lw=2, label='CORR-QMC')
+    plt.xlabel('sigma')
+    plt.ylabel('log-likelihood')
+    plt.legend()
+    if savefigs:
+        plt.savefig('loglik_interpolated_frozen_sqmc.pdf')
+
+    plt.show()
diff --git a/book/mle/malikpitt_interpolation.py b/book/mle/malikpitt_interpolation.py
index e04fb33..262ec2c 100644
--- a/book/mle/malikpitt_interpolation.py
+++ b/book/mle/malikpitt_interpolation.py
@@ -6,22 +6,45 @@
     We simply subclass SMC, so as to replace the resampling step by the
     interpolated version of Malik and Pitt. 
 
-""" 
-
+"""
+import numba as nb
 import particles
 
+
+@nb.njit
 def interpol(x1, x2, y1, y2, x):
     return y1 + (y2 - y1) * (x - x1) / (x2 - x1)
 
+
+@nb.njit
 def avg_n_nplusone(x):
     """ returns x[0]/2, (x[0]+x[1])/2, ... (x[-2]+x[-1])/2, x[-1]/2 
     """
 
     y = np.zeros(1 + x.shape[0])
-    hx = 0.5 * x 
+    hx = 0.5 * x
     y[:-1] = hx
-    y[1:] += hx 
-    return y 
+    y[1:] += hx
+    return y
+
+
+@nb.njit
+def _interpoled_resampling(cs, xs, u):
+    N = u.shape[0]
+    xrs = np.empty(N)
+    where = np.searchsorted(cs, u)
+    # costs O(N log(N)) but algorithm has O(N log(N)) complexity anyway
+    for n in range(N):
+        m = where[n]
+        if m == 0:
+            xrs[n] = xs[0]
+        elif m == N:
+            xrs[n] = xs[-1]
+        else:
+            xrs[n] = interpol(cs[m - 1], cs[m], xs[m - 1], xs[m], u[n])
+    return xrs
+
+
 
 def interpoled_resampling(W, x):
     """Resampling based on an interpolated CDF, as described in Malik and Pitt. 
@@ -39,24 +62,14 @@ def interpoled_resampling(W, x):
         the resampled particles 
 
     """
-    N = W.shape[0] 
+    N = W.shape[0]
     idx = np.argsort(x)
-    xs = x[idx] 
+    xs = x[idx]
     ws = W[idx]
-    cs = np.cumsum(avg_n_nplusone(ws)) 
+    cs = np.cumsum(avg_n_nplusone(ws))
     u = random.rand(N)
-    xrs = np.empty(N)
-    where = np.searchsorted(cs, u)
-    # costs O(N log(N)) but algorithm has O(N log(N)) complexity anyway
-    for n in range(N): 
-        m = where[n]
-        if m==0:
-            xrs[n] = xs[0]
-        elif m==N:
-            xrs[n] = xs[-1]
-        else:
-            xrs[n] = interpol(cs[m-1], cs[m], xs[m-1], xs[m], u[n])
-    return xrs
+    return _interpoled_resampling(cs, xs, u)
+
 
 class MalikPitt_SMC(particles.SMC):
     """Subclass of SMC that implements Malik and Pitt's interpolated resampling
@@ -65,10 +78,10 @@ class MalikPitt_SMC(particles.SMC):
     May be instantiated exactly as particles.SMC. Note however this works only for 
     univariate state-space models. 
     """
+
     def resample_move(self):
         self.rs_flag = True
         self.Xp = interpoled_resampling(self.W, self.X)
         self.reset_weights()
         # A not properly defined in this case 
         self.X = self.fk.M(self.t, self.Xp)
-
diff --git a/particles/ensemble_transform.py b/particles/ensemble_transform.py
new file mode 100644
index 0000000..a09d3e6
--- /dev/null
+++ b/particles/ensemble_transform.py
@@ -0,0 +1,137 @@
+# -*- coding: utf-8 -*-
+
+"""
+Ensemble Transform, as introduced by S. Reich in https://arxiv.org/abs/1210.0375
+
+Overview
+========
+
+This module implements the optimal transport based ensemble transform introduced by S. Reich in 2013.
+It depends on the particles weights as well as their location. This module depends on the POT library https://github.com/PythonOT/POT.
+
+    from particles import ensemble_transform as et
+
+Ensemble transform schemes
+==================
+
+All the ensemble transform schemes are implemented as classes with the following
+signature::
+
+    class EnsembleTransform:
+        def __init__(self, *args, **kwargs):
+            # Hyperparameters of the method if needed
+            pass
+
+        def resample(self, W, X):
+            return X-like array
+
+where:
+
+  * ``W`` is a vector of N normalised weights (i.e. positive and summing to
+    one).
+
+  * ``X`` (ndarray) are the locations of the particles
+
+
+Here the list of currently implemented resampling schemes:
+
+* `EnsembleTransform`
+* `EnsembleTransformFilter`
+
+"""
+
+import numpy as np
+import ot
+import scipy.linalg as lin
+
+from .core import SMC
+
+
+class EnsembleTransformBase(object):
+    def resample(self, W, X):
+        """docstring to do"""
+        raise NotImplementedError
+
+
+class EnsembleTransform(EnsembleTransformBase):
+    """Method from S. Reich, A non-parametric ensemble transform method for Bayesian inference (2013)"""
+
+    def __init__(self, metric="sqeuclidean", p=2):
+        self.metric = metric
+        self.p = p
+
+    def resample(self, W, X):
+        N = W.shape[0]
+        uniform_weights = np.full_like(W, 1 / N)
+        M = ot.utils.cdist(X.reshape(N, -1), X.reshape(N, -1), self.metric, p=self.p)
+        return N * ot.lp.emd(uniform_weights, W, M) @ X
+
+
+class EntropyRegularizedEnsembleTransform(EnsembleTransformBase):
+    """Method from S. Reich, A non-parametric ensemble transform method for Bayesian inference (2013)"""
+
+    def __init__(self, epsilon=0.1, metric="sqeuclidean", p=2):
+        self.metric = metric
+        self.p = p
+        self.epsilon = epsilon
+
+    def resample(self, W, X):
+        N = W.shape[0]
+        uniform_weights = np.full_like(W, 1 / N)
+        M = ot.utils.cdist(X.reshape(N, -1), X.reshape(N, -1), self.metric, p=self.p)
+        new_X = N * ot.bregman.sinkhorn(uniform_weights, W, M, self.epsilon) @ X
+        return new_X
+
+
+class CorrectedEntropyRegularizedEnsembleTransform(EnsembleTransformBase):
+    """Method from S. Reich, A non-parametric ensemble transform method for Bayesian inference (2013)"""
+
+    def __init__(self, epsilon=0.1, metric="sqeuclidean", p=2):
+        self.metric = metric
+        self.p = p
+        self.epsilon = epsilon
+
+    def resample(self, W, X):
+        N = W.shape[0]
+        uniform_weights = np.full_like(W, 1 / N)
+        M = ot.utils.cdist(X.reshape(N, -1), X.reshape(N, -1), self.metric, p=self.p)
+        transport = N * ot.bregman.sinkhorn(uniform_weights, W, M, self.epsilon)
+        diag_W = np.diag(W)
+        A = N * diag_W - transport.T.dot(transport)
+        delta = lin.solve_continuous_are(-transport, np.eye(N), A, np.eye(N))
+        new_X = (transport + delta).dot(X)
+        return new_X
+
+
+class EnsembleTransformFilter(SMC):
+    """Subclass of SMC that implements ensemble transform techniques (that depend on the particles locations as well as
+     their weights).
+     Takes an EnsembleTransformBase instance and the SMC arguments (resampling is ignored).
+
+    """
+
+    def __init__(self,
+                 ensemble_transform=EnsembleTransform(),
+                 fk=None,
+                 N=100,
+                 qmc=False,
+                 ESSrmin=0.5,
+                 store_history=False,
+                 verbose=False,
+                 summaries=True,
+                 **sum_options):
+        super(EnsembleTransformFilter, self).__init__(fk, N, qmc, None, ESSrmin, store_history, verbose, summaries,
+                                                      **sum_options)
+
+        self.ensemble_transform = ensemble_transform
+
+    def resample_move(self):
+        self.rs_flag = self.aux.ESS < self.N * self.ESSrmin
+        if self.rs_flag:  # if resampling
+            self.Xp = self.ensemble_transform.resample(self.aux.W, self.X)
+            self.reset_weights()
+            self.X = self.fk.M(self.t, self.Xp)
+
+    def resample_move_qmc(self):
+        self.rs_flag = True  # we *always* resample in SQMC
+        self.resample_move()
diff --git a/particles/utils.py b/particles/utils.py
index 9213ab5..13d8306 100644
--- a/particles/utils.py
+++ b/particles/utils.py
@@ -61,7 +61,7 @@ def f(x=0, y=0, z=0):
     set option ``seeding`` to True (otherwise, you will get
     identical results from all your workers); (b) make sure the function f does
     not rely on scipy frozen distributions, as these distributions also
-    freezes the states. For instance, do not use any frozen distribution when
+    freeze the states. For instance, do not use any frozen distribution when
     defining your own Feynman-Kac object. 
 
 .. seealso :: `multiSMC`
diff --git a/requirements.md b/requirements.md
index fd96a26..5fd0e03 100644
--- a/requirements.md
+++ b/requirements.md
@@ -3,3 +3,4 @@ numpy
 scipy
 ipython
 nbsphinx
+pot