diff --git a/mri/dictionary_learning/linear.py b/mri/dictionary_learning/linear.py
index d77b306e..ed76e339 100644
--- a/mri/dictionary_learning/linear.py
+++ b/mri/dictionary_learning/linear.py
@@ -78,8 +78,8 @@ def _op(self, dictionary, image):  # XXX works for square patches only!
         image: ndarray
             Input data array, a 2D image.
 
-        Return
-        ------
+        Returns
+        -------
         coeffs: ndarray of floats, 2d matrix dim nb_patches*nb_components,
                 the sparse coefficients.
         """
@@ -92,20 +92,18 @@ def op(self, image):
         This method returns the representation of the input data in the
         learnt dictionary, that is to say the sparse coefficients.
 
+        Remark: This method only works for squared patches
+
         Parameters
         ----------
         image: ndarray
             Input data array, a 2D image.
 
-        Return
-        ------
+        Returns
+        -------
         coeffs: ndarray of complex if is_complex, default(float)
                 2d matrix dim nb_patches*nb_components, the sparse
                 coefficients.
-
-        Remark
-        -------
-        This method only works for squared patches
         """
         if self.is_complex:
             return self._op(self.dictionary_r, image)
@@ -123,6 +121,8 @@ def _adj_op(self, coeffs, atoms, dtype="array"):
         This method returns the reconsructed image from the sparse
         coefficients.
 
+        Remark: This method only works for squared patches
+
         Parameters
         ----------
         coeffs: ndarray of floats,
@@ -135,13 +135,9 @@ def _adj_op(self, coeffs, atoms, dtype="array"):
             if 'array' return the data as a ndarray, otherwise return a
             pysap.Image.
 
-        Return
-        ------
-        ndarray, the reconstructed data.
-
-        Remark
+        Returns
         -------
-        This method only works for squared patches
+        ndarray, the reconstructed data.
         """
         image = numpy.dot(coeffs, atoms)
         image = image.reshape(image.shape[0], *self.patches_shape)
@@ -153,6 +149,8 @@ def adj_op(self, coeffs, dtype="array"):
         This method returns the reconsructed image from the sparse
         coefficients.
 
+        Remark: This method only works for squared patches
+
         Parameters
         ----------
         coeffs: ndarray of floats,
@@ -165,10 +163,6 @@ def adj_op(self, coeffs, dtype="array"):
         Returns
         -------
         ndarray, the reconstructed data.
-
-        Remark
-        -------
-        This method only works for squared patches
         """
         image_r = self._adj_op(numpy.real(coeffs),
                                self.dictionary_r.components_,
diff --git a/mri/dictionary_learning/utils.py b/mri/dictionary_learning/utils.py
index 64de9b2f..c0d3e4d0 100644
--- a/mri/dictionary_learning/utils.py
+++ b/mri/dictionary_learning/utils.py
@@ -28,6 +28,7 @@
 
 def timer(start, end):
     """ Give duration time between 2 times in hh:mm:ss.
+
     Parameters
     ----------
     start: float
@@ -44,6 +45,7 @@ def timer(start, end):
 
 def min_max_normalize(img):
     """ Center and normalize the given array.
+
     Parameters
     ----------
     img: np.ndarray
@@ -60,12 +62,12 @@ def extract_patches_from_2d_images(img, patch_shape):
 
     Parameters
     ----------
-        img: np.ndarray of floats, the input 2d image
+    img: np.ndarray of floats, the input 2d image
         patch_shape: tuple of int, shape of the patches
     Returns
     -------
-        patches: np.ndarray of floats, a 2d matrix with
-        -        dim nb_patches*(patch.shape[0]*patch_shape[1])
+    patches: np.ndarray of floats, a 2d matrix with
+    -        dim nb_patches*(patch.shape[0]*patch_shape[1])
     """
     patches = extract_patches_2d(img, patch_shape)
     patches = patches.reshape(patches.shape[0], -1)
@@ -74,20 +76,20 @@ def extract_patches_from_2d_images(img, patch_shape):
 
 def generate_flat_patches(images, patch_size, option='real'):
     """ Generate flat patches from the real/imaginary/complex images from the
-    list of images
+    list of images.
 
     Parameters
     ----------
-        image: list of list of np.ndarray of float or complex
-            a sublist containing all the images for one subject
-        patch_size: int,
-            width of square patches
-        option: 'real' (default),
-            'imag' real/imaginary part or 'complex'
-    Return
-    ------
-        flat_patches: list of np.ndarray as a GENERATOR
-            The patches flat and concatained as a list
+    image: list of list of np.ndarray of float or complex
+        a sublist containing all the images for one subject
+    patch_size: int,
+        width of square patches
+    option: 'real' (default),
+        'imag' real/imaginary part or 'complex'
+    Returns
+    -------
+    flat_patches: list of np.ndarray as a GENERATOR
+        The patches flat and concatained as a list
     """
     patch_shape = (patch_size, patch_size)
     flat_patches = images[:]
@@ -151,8 +153,8 @@ def learn_dictionary(flat_patches_subjects, nb_atoms=100, alpha=1, n_iter=1,
     verbose: int default1,
         The level of verbosity
 
-    Return
-    ------
+    Returns
+    -------
         dico: MiniBatchDictionaryLearning object
     """
     dico = MiniBatchDictionaryLearning(
diff --git a/mri/gridsearch.py b/mri/gridsearch.py
index 8fa2441c..dc4df7c6 100644
--- a/mri/gridsearch.py
+++ b/mri/gridsearch.py
@@ -51,8 +51,8 @@ def _default_wrapper(recons_func, **kwargs):
 def grid_search(func, param_grid, wrapper=None, n_jobs=1, verbose=0):
     """ Run `func` on the carthesian product of `param_grid`.
 
-    Parameters:
-    -----------
+    Parameters
+    ----------
     func: function
         The reconstruction function from whom to tune the hyperparameters.
         `func` return should be handle by wrapper if it's not a
@@ -79,8 +79,8 @@ def grid_search(func, param_grid, wrapper=None, n_jobs=1, verbose=0):
         messages increases with the verbosity level. If it more than 10,
         all iterations are reported.
 
-    Results:
-    --------
+    Results
+    -------
     metrics: dict
         the gridsearch results. Each key corresponds to a gridsearch parameters
         set. The values are 'params' the reconstruction parameters, 'image'
diff --git a/mri/parallel_mri/extract_sensitivity_maps.py b/mri/parallel_mri/extract_sensitivity_maps.py
index 5754dc2d..35ce3378 100644
--- a/mri/parallel_mri/extract_sensitivity_maps.py
+++ b/mri/parallel_mri/extract_sensitivity_maps.py
@@ -78,6 +78,7 @@ def gridded_inverse_fourier_transform_nd(kspace_loc,
     method: {'linear', 'nearest', 'cubic'}
         Method of interpolation for more details see scipy.interpolate.griddata
         documentation
+
     Returns
     -------
     np.ndarray
@@ -100,8 +101,9 @@ def get_Smaps(k_space, img_shape, samples, thresh,
     acquisition and for variable density sampling scheme where teh k-space
     center had been heavily sampled.
     Reference : Self-Calibrating Nonlinear Reconstruction Algorithms for
-                Variable Density Sampling and Parallel Reception MRI
-                https://ieeexplore.ieee.org/abstract/document/8448776
+    Variable Density Sampling and Parallel Reception MRI
+    https://ieeexplore.ieee.org/abstract/document/8448776
+
     Parameters
     ----------
     k_space: np.ndarray
diff --git a/mri/parallel_mri/gradient.py b/mri/parallel_mri/gradient.py
index 01ffaedf..f6a58bb6 100644
--- a/mri/parallel_mri/gradient.py
+++ b/mri/parallel_mri/gradient.py
@@ -92,7 +92,7 @@ def _analy_rsns_op_method(self, x):
 class Gradient_pMRI_synthesis(GradBasic, PowerMethod):
     """ Gradient synthesis class.
 
-    This class defines the grad operators for |M*F*invL*alpha - data|**2.
+    This class defines the grad operators for abs(M*F*invL*alpha - data)**2.
 
     Parameters
     ----------
diff --git a/mri/reconstruct/cost.py b/mri/reconstruct/cost.py
index cbd26326..a263f074 100644
--- a/mri/reconstruct/cost.py
+++ b/mri/reconstruct/cost.py
@@ -26,8 +26,8 @@ def __init__(self, linear_op, initial_cost=1e6, tolerance=1e-4,
                  plot_output=None):
         """ Initialize the 'DualGapCost' class.
 
-        Parameters:
-        -----------
+        Parameters
+        ----------
         x: np.ndarray
             input original data array.
         costFunc: class
@@ -80,8 +80,9 @@ def __init__(self, gradient_op, prox_op, initial_cost=1e6,
                  tolerance=1e-4, cost_interval=1, test_range=4, verbose=False,
                  plot_output=None):
         """ Initialize the 'Cost' class.
-        Parameters:
-        -----------
+
+        Parameters
+        ----------
         gradient_op: instance of the gradient operator
             gradient operator used in the reconstruction process. It must
             implements the get_cost_function.
@@ -127,10 +128,12 @@ def __init__(self, gradient_op, prox_op, initial_cost=1e6,
 
     def _calc_cost(self, x_new, *args, **kwargs):
         """ Return the cost.
+
         Parameters
         ----------
         x_new: np.ndarray
             intermediate solution in the optimization problem.
+
         Returns
         -------
         cost: float
diff --git a/mri/reconstruct/fourier.py b/mri/reconstruct/fourier.py
index 34b0e915..d52f7510 100644
--- a/mri/reconstruct/fourier.py
+++ b/mri/reconstruct/fourier.py
@@ -33,6 +33,7 @@ class FourierBase(object):
     """
     def op(self, img):
         """ This method calculates Fourier transform.
+
         Parameters
         ----------
         img: np.ndarray
diff --git a/mri/reconstruct/linear.py b/mri/reconstruct/linear.py
index f1a194c5..84cc9e19 100644
--- a/mri/reconstruct/linear.py
+++ b/mri/reconstruct/linear.py
@@ -289,10 +289,12 @@ def adj_op(self, coefs):
 
     def l2norm(self, shape):
         """ Compute the L2 norm.
+
         Parameters
         ----------
         shape: uplet
             the data shape.
+
         Returns
         -------
         norm: float
diff --git a/mri/reconstruct/utils.py b/mri/reconstruct/utils.py
index 9a282647..ef6c0f2f 100644
--- a/mri/reconstruct/utils.py
+++ b/mri/reconstruct/utils.py
@@ -80,15 +80,15 @@ def normalize_frequency_locations(samples, Kmax=None):
     This function normalize the samples locations between [-0.5; 0.5[ for
     the non-cartesian case
 
-    Parameters:
-    -----------
+    Parameters
+    ----------
     samples: np.ndarray
         Unnormalized samples
     Kmax: int, float, array-like or None
         Maximum Frequency of the samples locations is supposed to be equal to
         base Resolution / (2* Field of View)
 
-    Return:
+    Returns
     -------
     normalized_samples: np.ndarray
         Same shape as the parameters but with values between [-0.5; 0.5[