diff --git a/2d_remote.py b/2d_remote.py
index f71e516..39f81b7 100644
--- a/2d_remote.py
+++ b/2d_remote.py
@@ -14,7 +14,7 @@
 import itertools
 
 #from scipy.optimize import minimize
-  
+
 #%% Approx \nabla \hat{f} point by point.
 
 mu_1 = np.array([-1., 0.])
@@ -30,12 +30,12 @@
 y = np.linspace(-bound_y, bound_y, meshsize)
 xv, yv = np.meshgrid(x, y)
 
-mu_vec = [(mux, muy) for mux, muy in itertools.product(
+mu_vec = list(itertools.product(
            np.linspace(-bound_x, bound_x, num=nb_bases),
-           np.linspace(-bound_y, bound_y, num=nb_bases))]
+           np.linspace(-bound_y, bound_y, num=nb_bases)))
 
 pi_mesh = 0.5*np.exp(-((xv-mu_1[0])**2+(yv-mu_1[1])**2)/(2*sigma2)) / (2*np.pi*sigma2) \
-        + 0.5*np.exp(-((xv-mu_2[0])**2+(yv-mu_2[1])**2)/(2*sigma2)) / (2*np.pi*sigma2)  
+        + 0.5*np.exp(-((xv-mu_2[0])**2+(yv-mu_2[1])**2)/(2*sigma2)) / (2*np.pi*sigma2)
 dxpi_mesh = (xv-mu_1[0])*np.exp(-((xv-mu_1[0])**2+(yv-mu_1[1])**2)/(2*sigma2)) + \
           (xv-mu_2[0])*np.exp(-((xv-mu_2[0])**2+(yv-mu_2[1])**2)/(2*sigma2))
 dxpi_mesh *= (-0.5)*sigma2**(-1)*(2*np.pi*sigma2)**(-1.)
@@ -112,7 +112,7 @@ def pi(x):
 var = var[:-1, :-1]
 area = (xv[0, 1] - xv[0, 0])*(yv[1, 0] - yv[0, 0])
 var = 2*np.sum(var)*area
-print('var: {}'.format(var))
+print(f'var: {var}')
            
   
 
diff --git a/reading.py b/reading.py
index c9ca89d..f7ae0ca 100644
--- a/reading.py
+++ b/reading.py
@@ -67,7 +67,7 @@ def tab(log_ula):
                 df2[compt] = tab_labels[j] 
                 df3[compt] = tab_titles[k]
                 compt +=1
-    
+
     df = np.rec.fromarrays((df1,df2,df3), names=(xnames[l], "method", "algorithm"))    
     dfp = pd.DataFrame(df)    
     ldfp.append(dfp)
@@ -131,14 +131,14 @@ def tab(log_ula):
     plt.boxplot(temp, labels=tab_labels)
 #    plt.ylim(mini, maxi)
     if (k % 3)==0:
-        st = r"$x_" + str(k // 3 +1) + "$"
+        st = f"$x_{str(k // 3 +1)}$"
         plt.ylabel(st)
     plt.title(tab_titles[k % 3])
     plt.grid(True)
 
 #plt.suptitle("Ill conditionned Gaussian, first coordinate, error " \
 #          + indi +" moment, N=10**6, dimension " + str(d) + " , x0=" + str(x0Tab[x0_ind]))
-plt.show() 
+plt.show()
 f.savefig("log-1-1.pdf", bbox_inches='tight')
 
 # \beta^2_1, ..., \beta^2_d
@@ -151,14 +151,14 @@ def tab(log_ula):
     temp = np.transpose(tab_log[::2,:,4 + k // 3,k % 3])
     plt.boxplot(temp, labels=tab_labels)
     if (k % 3)==0:
-        st = r"$(x_" + str(k // 3 +1) + ")^2$"
+        st = f"$(x_{str(k // 3 +1)})^2$"
         plt.ylabel(st)
     plt.title(tab_titles[k % 3])
     plt.grid(True)
 
 #plt.suptitle("Ill conditionned Gaussian, first coordinate, error " \
 #          + indi +" moment, N=10**6, dimension " + str(d) + " , x0=" + str(x0Tab[x0_ind]))
-plt.show() 
+plt.show()
 f.savefig("log-2-1.pdf", bbox_inches='tight')
 
 
@@ -211,6 +211,7 @@ def tab(log_ula):
 indices = np.array([3,4,5]) #indices = np.array([0,1,4,5])
 #xnames = ["x_1", "x_2", "x_1^2", "x_2^2"]
 xnames = ["x_1^2", "x_2^2", "x_3^2"]
+nb_methods = 5
 #xnames_tex = ["$x_1$", "$x_2$", "$x_1^2$", "$x_2^2$"]
 
 for l in np.arange(3):
@@ -221,7 +222,6 @@ def tab(log_ula):
     df2 = np.empty(a.size, dtype=np.dtype('U25'))
     df3 = np.empty(a.size, dtype=np.dtype('U25'))
     compt = 0
-    nb_methods = 5
     tab_labels = ["O", "CV-1", "CV-2", "ZV-1", "ZV-2"]
     for i in np.arange(nc):
         for j in np.arange(nb_methods):
@@ -230,9 +230,9 @@ def tab(log_ula):
                 df2[compt] = tab_labels[j] 
                 df3[compt] = tab_titles[k]
                 compt +=1
-    
-    df = np.rec.fromarrays((df1,df2,df3), names=(xnames[l], "method", "algorithm"))    
-    dfp = pd.DataFrame(df)    
+
+    df = np.rec.fromarrays((df1,df2,df3), names=(xnames[l], "method", "algorithm"))
+    dfp = pd.DataFrame(df)
     ldfp.append(dfp)
 
 f = plt.figure(figsize=(16,16))
@@ -261,14 +261,14 @@ def tab(log_ula):
     temp = np.transpose(tab_pro[::2,:,k // 3,k % 3])
     plt.boxplot(temp, labels=tab_labels)
     if (k % 3)==0:
-        st = r"$x_" + str(k // 3 +1) + "$"
+        st = f"$x_{str(k // 3 +1)}$"
         plt.ylabel(st)
     plt.title(tab_titles[k % 3])
     plt.grid(True)
 
 #plt.suptitle("Ill conditionned Gaussian, first coordinate, error " \
 #          + indi +" moment, N=10**6, dimension " + str(d) + " , x0=" + str(x0Tab[x0_ind]))
-plt.show() 
+plt.show()
 f.savefig("pro-1-1.pdf", bbox_inches='tight')
 
 # \beta^2_1, ..., \beta^2_d
@@ -281,12 +281,12 @@ def tab(log_ula):
     temp = np.transpose(tab_pro[::2,:,3 + k // 3,k % 3])
     plt.boxplot(temp, labels=tab_labels)
     if (k % 3)==0:
-        st = r"$(x_" + str(k // 3 +1) + ")^2$"
+        st = f"$(x_{str(k // 3 +1)})^2$"
         plt.ylabel(st)
     plt.title(tab_titles[k % 3])
     plt.grid(True)
 
 #plt.suptitle("Ill conditionned Gaussian, first coordinate, error " \
 #          + indi +" moment, N=10**6, dimension " + str(d) + " , x0=" + str(x0Tab[x0_ind]))
-plt.show() 
+plt.show()
 f.savefig("pro-2-1.pdf", bbox_inches='tight')
\ No newline at end of file
diff --git a/two_dimension.py b/two_dimension.py
index 04aae89..e76dec6 100644
--- a/two_dimension.py
+++ b/two_dimension.py
@@ -14,7 +14,7 @@
 import itertools
 
 #from scipy.optimize import minimize
-  
+
 #%% Approx \nabla \hat{f} point by point.
 
 mu_1 = np.array([-1., 0.])
@@ -30,12 +30,12 @@
 y = np.linspace(-bound_y, bound_y, meshsize)
 xv, yv = np.meshgrid(x, y)
 
-mu_vec = [(mux, muy) for mux, muy in itertools.product(
+mu_vec = list(itertools.product(
            np.linspace(-bound_x, bound_x, num=nb_bases),
-           np.linspace(-bound_y, bound_y, num=nb_bases))]
+           np.linspace(-bound_y, bound_y, num=nb_bases)))
 
 pi_mesh = 0.5*np.exp(-((xv-mu_1[0])**2+(yv-mu_1[1])**2)/(2*sigma2)) / (2*np.pi*sigma2) \
-        + 0.5*np.exp(-((xv-mu_2[0])**2+(yv-mu_2[1])**2)/(2*sigma2)) / (2*np.pi*sigma2)  
+        + 0.5*np.exp(-((xv-mu_2[0])**2+(yv-mu_2[1])**2)/(2*sigma2)) / (2*np.pi*sigma2)
 dxpi_mesh = (xv-mu_1[0])*np.exp(-((xv-mu_1[0])**2+(yv-mu_1[1])**2)/(2*sigma2)) + \
           (xv-mu_2[0])*np.exp(-((xv-mu_2[0])**2+(yv-mu_2[1])**2)/(2*sigma2))
 dxpi_mesh *= (-0.5)*sigma2**(-1)*(2*np.pi*sigma2)**(-1.)
@@ -56,14 +56,14 @@
   dypois[:, :, k] = dyp
   dxdxpois = -(2*np.pi*s2)**(-1)*s2**(-1)*(xv-mu[0])*np.exp(-((xv-mu[0])**2 + (yv-mu[1])**2)/(2*s2))
   dydypois = -(2*np.pi*s2)**(-1)*s2**(-1)*(yv-mu[1])*np.exp(-((xv-mu[0])**2 + (yv-mu[1])**2)/(2*s2))
-  
+
   dxU_dxpois = dxU_mesh * dxp
   dyU_dypois = dyU_mesh * dyp
-  
+
   Lx_pois[:, :, k] = - dxU_dxpois + dxdxpois
   Ly_pois[:, :, k] = - dyU_dypois + dydypois
-  
-Lx_pois_reshaped = np.reshape(Lx_pois, (meshsize**2, nb_bases**2))    
+
+Lx_pois_reshaped = np.reshape(Lx_pois, (meshsize**2, nb_bases**2))
 Ly_pois_reshaped = np.reshape(Ly_pois, (meshsize**2, nb_bases**2))    
 
 L_pois = np.hstack((Lx_pois_reshaped, Ly_pois_reshaped))
@@ -136,7 +136,7 @@ def pi(x):
 plt.colorbar()
 plt.show()
 #fig.savefig("approx_pi_Lpois_2d_nabla_rcond5.jpeg", bbox_inches='tight')
-  
+
 #fig = plt.figure(figsize=(10,10))
 #plt.rcParams.update({'font.size': 16}) # default 10
 #plt.pcolormesh(xv, yv, pi_mesh)
@@ -144,13 +144,13 @@ def pi(x):
 #plt.title(r'$\pi$')
 #plt.show()
 #fig.savefig('density_pi_2d.pdf', bbox_inches='tight')
-           
-  
+
+
 var = (dxpois**2 + dypois**2)*pi_mesh
 var = var[:-1, :-1]
 area = (xv[0, 1] - xv[0, 0])*(yv[1, 0] - yv[0, 0])
 var = 2*np.sum(var)*area
-print('var: {}'.format(var))
+print(f'var: {var}')
 
 def U(x):
   return -np.log(pi(x))
diff --git a/two_dimension_article.py b/two_dimension_article.py
index 88e7f94..871136d 100644
--- a/two_dimension_article.py
+++ b/two_dimension_article.py
@@ -14,7 +14,7 @@
 import itertools
 
 #from scipy.optimize import minimize
-  
+
 #%% Approx \nabla \hat{f} point by point.
 
 mu_1 = np.array([-1., 0.])
@@ -30,12 +30,12 @@
 y = np.linspace(-bound_y, bound_y, meshsize)
 xv, yv = np.meshgrid(x, y)
 
-mu_vec = [(mux, muy) for mux, muy in itertools.product(
+mu_vec = list(itertools.product(
            np.linspace(-bound_x, bound_x, num=nb_bases),
-           np.linspace(-bound_y, bound_y, num=nb_bases))]
+           np.linspace(-bound_y, bound_y, num=nb_bases)))
 
 pi_mesh = 0.5*np.exp(-((xv-mu_1[0])**2+(yv-mu_1[1])**2)/(2*sigma2)) / (2*np.pi*sigma2) \
-        + 0.5*np.exp(-((xv-mu_2[0])**2+(yv-mu_2[1])**2)/(2*sigma2)) / (2*np.pi*sigma2)  
+        + 0.5*np.exp(-((xv-mu_2[0])**2+(yv-mu_2[1])**2)/(2*sigma2)) / (2*np.pi*sigma2)
 dxpi_mesh = (xv-mu_1[0])*np.exp(-((xv-mu_1[0])**2+(yv-mu_1[1])**2)/(2*sigma2)) + \
           (xv-mu_2[0])*np.exp(-((xv-mu_2[0])**2+(yv-mu_2[1])**2)/(2*sigma2))
 dxpi_mesh *= (-0.5)*sigma2**(-1)*(2*np.pi*sigma2)**(-1.)
@@ -115,13 +115,13 @@ def pi(x):
 #  print('------------------------')
 #  print('rcond: {}'.format(rcond))
 #  print('------------------------')
-  
+
 #  coeffs_pois = np.linalg.lstsq(L_pois, - f_tilde_flatten, rcond=rcond)[0]
 #
 #  dxpois = np.dot(dxpois_mesh, coeffs_pois[:nb_bases**2])
 #  dypois = np.dot(dypois_mesh, coeffs_pois[nb_bases**2:])
 #  Lpois = np.dot(Lx_pois, coeffs_pois[:nb_bases**2]) + np.dot(Ly_pois, coeffs_pois[nb_bases**2:])
-  
+
 #  coeffs_pois = np.linalg.lstsq(L_pois_reshaped, - f_tilde_flatten, rcond=rcond)[0]
 
 #  dxpois = np.dot(dxpois_mesh, coeffs_pois)
@@ -132,12 +132,12 @@ def pi(x):
 var = var[:-1, :-1]
 area = (xv[0, 1] - xv[0, 0])*(yv[1, 0] - yv[0, 0])
 var = 2*np.sum(var)*area
-print('var: {}'.format(var))
+print(f'var: {var}')
 
-norml1 = np.linalg.norm(coeffs_pois, ord=1) / len(coeffs_pois)  
-print('norm coeffs pois: {}'.format(norml1))
+norml1 = np.linalg.norm(coeffs_pois, ord=1) / len(coeffs_pois)
+print(f'norm coeffs pois: {norml1}')
 err = np.sum(np.absolute(f_tilde_mesh + Lpois)) / meshsize**2
-print('err: {}'.format(err))
+print(f'err: {err}')
 
 fig = plt.figure(figsize=(16,16))
 plt.rcParams.update({'font.size': 13}) # default 10
@@ -159,7 +159,7 @@ def pi(x):
 plt.colorbar()
 plt.show()
 fig.savefig("ref_hatf.jpeg", bbox_inches='tight')
-  
+
 fig = plt.figure(figsize=(10,10))
 plt.rcParams.update({'font.size': 16}) # default 10
 plt.pcolormesh(xv, yv, pi_mesh)