diff --git a/docs/src/tutorial/optimization.md b/docs/src/tutorial/optimization.md
index 47f5784f..d072c3f3 100644
--- a/docs/src/tutorial/optimization.md
+++ b/docs/src/tutorial/optimization.md
@@ -2,6 +2,4 @@
 
 An interesting question in BosonSampling is to find interferometers that of maximize certain properties.
 
-<!-- For instance, consider a partition of all odd modes 1,3,... When using a `Fourier` interferometer, a strange phenomenon appears: -->
-
 We provide the function `minimize_over_unitary_matrices()` which operates a conjugate gradient algorithm for the optimization over unitary matrices. It is implemented from [Conjugate gradient algorithm for optimization under unitary matrix constraint](https://doi.org/10.1016/j.sigpro.2009.03.015) from Traian Abrudan, Jan Eriksson, Visa Koivunen.
diff --git a/src/permanent_conjectures/counter_example_functions.jl b/src/permanent_conjectures/counter_example_functions.jl
index 06fcb2f3..6e7d8a39 100644
--- a/src/permanent_conjectures/counter_example_functions.jl
+++ b/src/permanent_conjectures/counter_example_functions.jl
@@ -3,7 +3,7 @@
 
 """
 
-	 cholesky_semi_definite_positive(A)
+    cholesky_semi_definite_positive(A)
 
 cholesky decomposition (`A` = R' * R) for a sdp but not strictly positive
 definite matrix
@@ -31,7 +31,7 @@ end
 
 """
 
-	incorporate_in_a_unitary(X)
+    incorporate_in_a_unitary(X)
 
 incorporates the renormalized matrix X in a double sized unitary through the
 proof of Lemma 29  of Aaronson Arkipov seminal [The Computational Complexity of Linear Optics](https://arxiv.org/abs/1011.3245)
diff --git a/src/permanent_conjectures/permanent_on_top.jl b/src/permanent_conjectures/permanent_on_top.jl
index 6bd5236a..917fb153 100644
--- a/src/permanent_conjectures/permanent_on_top.jl
+++ b/src/permanent_conjectures/permanent_on_top.jl
@@ -2,6 +2,9 @@
 # instead of pure state inputs
 
 """
+
+    schur_matrix(H)
+
 computes the Schur matrix as defined in Eq. 1 of Linear Algebra and its Applications 490 (2016) 196–201
 """
 function schur_matrix(H)
@@ -29,11 +32,15 @@ function multiply_permutations(a,b)
     a[b[:]]
 end
 
+"""
+
+    J_array(theta, n)
+
+returns the `J` as defined in in Eq.10 of [Universality of Generalized Bunching and Efficient Assessment of Boson Sampling](https://arxiv.org/abs/1509.01561), with the permutations coming in the order given by `permutations(collect(1:n))`
+"""
 function J_array(theta, n)
 
-    """
-    returns the J as defined in in Eq.10 of https://arxiv.org/abs/1509.01561, with the permutations coming in the order given by permutations(collect(1:n))
-    """
+
 
     J = zeros(eltype(theta), factorial(n))
 
@@ -49,9 +56,13 @@ function J_array(theta, n)
 
 end
 
-function density_matrix_from_J(J,n)
 
-    """density matrix associated to a J function as defined in https://arxiv.org/abs/1509.01561, computed through Eq. 46"""
+"""
+    density_matrix_from_J(J,n)
+
+density matrix associated to a J function as defined in [Universality of Generalized Bunching and Efficient Assessment of Boson Sampling](https://arxiv.org/abs/1509.01561), computed through Eq. 46.
+"""
+function density_matrix_from_J(J,n)
 
     density_matrix = zeros(eltype(J), factorial(n), factorial(n))
     permutation_array = collect(permutations(collect(1:n)))