diff --git a/input/thesis_part1/plots/directional_plot_partial.py b/input/thesis_part1/plots/directional_plot_partial.py
index 41c8a8a8082e1b071aa1186989a9273125f0d9e2..ab9b2d2a625d4912e211c888468bbda70d73e159 100644
--- a/input/thesis_part1/plots/directional_plot_partial.py
+++ b/input/thesis_part1/plots/directional_plot_partial.py
@@ -10,7 +10,7 @@ import matplotlib.pyplot as plt
 from matplotlib import rc
 plt.rcParams['ytick.right'] = plt.rcParams['ytick.labelright'] = True
 plt.rcParams['ytick.left'] = plt.rcParams['ytick.labelleft'] = False
-rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb} \renewcommand{\familydefault}{\sfdefault} \usepackage{sansmathfonts}')
 rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
 rc('text', usetex=True)
 H = 3.0;
@@ -31,7 +31,7 @@ def plot(data):
     
     plt.plot(x[:,data.idx_param],x[:,2*data.n_param],color='k',marker='o',markersize=3,linewidth=0)
     plt.xlabel(r'$'+data.var+'$')
-    plt.ylabel(r'$\mathcal{J}('+data.var+')$')
+    plt.ylabel(r'$J('+data.var+')$')
     plt.grid()
     plt.savefig(prepath+data.path+data.filename+'_misfit.eps')
     plt.close()
@@ -65,7 +65,7 @@ def plot(data):
     plt.grid()
     #plt.legend(loc=3,ncol=1);
     plt.xlabel(r'$'+data.var+'$')
-    plt.ylabel(r'$D_{'+data.var+'} \mathcal{J}('+data.var+')$')
+    plt.ylabel(r'$D_{'+data.var+'} J('+data.var+')$')
     plt.savefig(prepath+data.path+data.filename+'_dj.eps')
     plt.close()
     
@@ -96,7 +96,7 @@ def plot(data):
     plt.grid()
     #plt.legend(loc=3,ncol=1);
     plt.xlabel(r'$'+data.var+'$')
-    plt.ylabel(r'$D^2_{'+data.var+'} \mathcal{J}('+data.var+')$')
+    plt.ylabel(r'$D^2_{'+data.var+'} J('+data.var+')$')
     plt.savefig(prepath+data.path+data.filename+'_d2j.eps')
     plt.close()
     
diff --git a/inversion_plot_comp_part2.py b/inversion_plot_comp_part2.py
new file mode 100644
index 0000000000000000000000000000000000000000..b7b59fe6a7a85f011be1142ef68d61b5e79f7513
--- /dev/null
+++ b/inversion_plot_comp_part2.py
@@ -0,0 +1,479 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Aug  3 16:57:28 2022
+
+@author: xavier
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import rc
+plt.rcParams['ytick.right'] = plt.rcParams['ytick.labelright'] = False
+plt.rcParams['ytick.left'] = plt.rcParams['ytick.labelleft'] = True
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb} \renewcommand{\familydefault}{\sfdefault} \usepackage{sansmathfonts}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
+rc('text', usetex=True)
+lw = 1.5
+H = 2.0;
+L = 3.0;
+
+"""
+#CS1: IP
+j0 = 35.6633
+n_comp = 4
+path = ['extended/US_lBFGS/','extended/S_lBFGS/','extended/SST_lBFGS/','extended/SSM_lBFGS/']
+filename = ['paper2_extended_inversion_lbfgs','paper2_extended_inversion_lbfgs_w','paper2_extended_inversion_lbfgs_w_stab','paper2_extended_inversion_lbfgs_w_smooth']
+lab = ['unweight','weight','stab','smooth']
+glob = ['linesearch','linesearch','linesearch','linesearch']
+order = [1, 1, 1, 1]
+col = ['k','k', 'k', 'k']
+mrk = ['']*4
+lin = [':','-.','--','-']
+max_forward = 80
+Δmax_forward = 10
+
+n_freq = 1
+"""
+"""
+#CS1: GD
+j0 = 35.6633
+n_comp = 7
+path = ['extended/SST_GD_LS/','extended/SST_GD_TRA/','extended/SST_GD_BTRB/','extended/SST_GD_BTRC/','extended/SST_GD_BTRAR/','extended/SST_GD_BTRBR/','extended/SST_GD_BTRCR/']
+filename = ['paper2_extended_sst_gd_ls','paper2_extended_sst_gd_trA','paper2_extended_sst_gd_trB','paper2_extended_sst_gd_trC','paper2_extended_sst_gd_trAR','paper2_extended_sst_gd_trBR','paper2_extended_sst_gd_trCR']
+lab = ['LS','TR-P(A)','TR-P(B)','TR-P(C)','TR-R(A)','TR-R(B)','TR-R(C)']
+glob = ['linesearch'] + ['trust-region'] * 3 + ['trust-region-retro'] * 3
+order = [1] * 7
+col = ['k'] + ['C0','C1','C2'] + ['C3','C4','C5']
+mrk = ['']*7
+#lin = ['-']*7 
+lin = [(0, (4, 4)), (0, (4, 4)), (0, (4, 4)),(4, (4, 4)), (0, (4, 8)), (4, (4, 8)), (8, (4, 8))]
+max_forward = 360
+Δmax_forward = 40
+
+n_freq = 1
+"""
+
+"""
+#CS1: l-BFGS
+j0 = 35.6633
+n_comp = 7
+path = ['extended/SST_lBFGS/','extended/SST_lBFGS_TRA/','extended/SST_lBFGS_TRB/','extended/SST_lBFGS_TRC/','extended/SST_lBFGS_TRAR/','extended/SST_lBFGS_TRBR/','extended/SST_lBFGS_TRCR/']
+filename = ['paper2_extended_sst_lbfgs_ls','paper2_extended_sst_lbfgs_trA','paper2_extended_sst_lbfgs_trB','paper2_extended_sst_lbfgs_trC','paper2_extended_sst_lbfgs_trAR','paper2_extended_sst_lbfgs_trBR','paper2_extended_sst_lbfgs_trCR']
+lab = ['LS','TR-P(A)','TR-P(B)','TR-P(C)','TR-R(A)','TR-R(B)','TR-R(C)']
+glob = ['linesearch'] + ['trust-region'] * 3 + ['trust-region-retro'] * 3
+order = [1.5] * 7
+col = ['k'] + ['C0','C1','C2'] + ['C3','C4','C5']
+mrk = ['']*7
+#lin = ['-'] * 7
+lin = [(0, (4, 4)), (0, (4, 4)), (0, (4, 8)),(0, (4, 4)), (4, (4, 4)), (4, (4, 8)), (8, (4, 8))]
+max_forward = 80
+Δmax_forward = 10
+
+n_freq = 1
+"""
+"""
+#CS1: Newton
+j0 = 35.6633
+n_comp = 14
+path = ['extended/SST_NCG/','extended/SST_NSCGA/','extended/SST_NSCGB/','extended/SST_NSCGC/','extended/SST_NSCGAR/','extended/SST_NSCGBR/','extended/SST_NSCGCR/','extended/SST_NCGGN/','extended/SST_NSCGAGN/','extended/SST_NSCGBGN/','extended/SST_NSCGCGN/','extended/SST_NSCGARGN/','extended/SST_NSCGBRGN/','extended/SST_NSCGCRGN/']
+filename = ['paper2_extended_sst_ncg','paper2_extended_sst_nscgA','paper2_extended_sst_nscgB','paper2_extended_sst_nscgC','paper2_extended_sst_nscgAr','paper2_extended_sst_nscgBr','paper2_extended_sst_nscgC','paper2_extended_sst_ncggn','paper2_extended_sst_nscgAgn','paper2_extended_sst_nscgBgn','paper2_extended_sst_nscgCgn','paper2_extended_sst_nscgArgn','paper2_extended_sst_nscgBrgn','paper2_extended_sst_nscgCrgn']
+lab = ['ncg','nscgA','nscgB','nscgC','nscgAr','nscgBr','nscgCr','ncggn','nscgAgn','nscgBgn','nscgCgn','nscgArgn','nscgBrgn','nscgCrgn']
+glob = (['linesearch'] + ['trust-region'] * 3 + ['trust-region-retro'] * 3) * 2
+order = [2] * 14
+col = ['k', 'C0','C1','C2', 'C3', 'C4', 'C5'] * 2
+mrk = ['']*14 
+lin = [(0,(4,4))]*14
+
+lin[0] = (0,(4,4))
+#lin[1] = (0,(4,4))
+#lin[2] = (2,(4,4))
+#lin[3] = (4,(4,4))
+
+#lin[0] = (4,(4,4))
+lin[1] = (2,(4,4))
+lin[2] = (0,(4,4))
+lin[3] = (4,(4,4))
+lin[4] = (0,(4,4))
+lin[5] = (4,(4,4))
+lin[4+7] = (0,(4,4))
+lin[5+7] = (4,(4,4))
+lin[0+7] = (4,(4,4))
+lin[2+7] = (0,(4,4))
+lin[3+7] = (4,(4,4))
+
+max_forward = 180
+Δmax_forward = 20
+
+n_freq = 1
+"""
+
+"""
+#CS2: IP
+j0 = 209.852
+n_comp = 4
+path = ['extended/concrete3/US_LBFGS_LS/','extended/concrete3/S_LBFGS_LS/','extended/concrete3/SST_LBFGS_LS/','extended/concrete3/SSM_LBFGS_LS/']
+filename = ['paper2_concrete3_us_lbfgs_ls','paper2_concrete3_s_lbfgs_ls','paper2_concrete3_sst_lbfgs_ls','paper2_concrete3_ssm_lbfgs_ls']
+lab = ['unweight','weight','stab','smooth']
+glob = ['linesearch','linesearch','linesearch','linesearch']
+order = [1, 1, 1, 1]
+col = ['k','k', 'k', 'k']
+mrk = ['']*4
+lin = [':','-.','--','-']
+lin[2] = (4,(4,4))
+max_forward = 200+10
+Δmax_forward = 25
+
+n_freq = 1
+"""
+
+"""
+#CS2: GD
+j0 = 209.852
+n_comp = 3
+path = ['extended/concrete3/SSM_GD_LS/','extended/concrete3/SSM_GD_BTRB/','extended/concrete3/SSM_GD_BTRBR/']
+filename = ['paper2_concrete3_ssm_gd_ls','paper2_concrete3_ssm_gd_trB','paper2_concrete3_ssm_gd_trBR']
+lab = ['LS','TR-P(B)','TR-R(B)']
+glob = ['linesearch'] + ['trust-region'] * 1 + ['trust-region-retro'] * 1
+order = [1] * 3
+col = ['k'] + ['C1'] + ['C4']
+mrk = ['']*3
+lin = ['-']*3 
+max_forward = 800
+Δmax_forward = 100
+n_freq = 1
+"""
+
+"""
+#CS2: l-BFGS
+j0 = 209.852
+n_comp = 3
+path = ['extended/concrete3/SSM_LBFGS_LS/','extended/concrete3/SSM_LBFGS_TRB/','extended/concrete3/SSM_LBFGS_TRBR/']
+filename = ['paper2_concrete3_ssm_lbfgs_ls','paper2_concrete3_ssm_lbfgs_trB','paper2_concrete3_ssm_lbfgs_trBR']
+lab = ['LS','TR-P(B)','TR-R(B)']
+glob = ['linesearch'] + ['trust-region'] * 1 + ['trust-region-retro'] * 1
+order = [1.5] * 3
+col = ['k'] + ['C1'] + ['C4']
+mrk = ['']*3
+lin = ['-']*3 
+max_forward = 200
+Δmax_forward = 25
+n_freq = 1
+"""
+
+"""
+#CS2: Newton
+j0 = 209.852
+n_comp = 6
+path = ['extended/concrete3/SSM_NCG/','extended/concrete3/SSM_NSCGB/','extended/concrete3/SSM_NSCGBR/','extended/concrete3/SSM_NCGGN/','extended/concrete3/SSM_NSCGBGN/','extended/concrete3/SSM_NSCGBRGN/']
+filename = ['paper2_concrete3_ssm_ncg','paper2_concrete3_ssm_nscgB','paper2_concrete3_ssm_nscgBR','paper2_concrete3_ssm_ncggn','paper2_concrete3_ssm_nscgBgn','paper2_concrete3_ssm_nscgBRgn']
+lab = ['FN-LS','FN-TR-P(B)','FN-TR-R(B)','GN-LS','GN-TR-P(B)','GN-TR-R(B)']
+glob = ['linesearch','trust-region','trust-region-retro'] * 2
+order = [2.] * 6
+col = ['k','C1','C4']*2
+mrk = ['o']*6
+lin = ['-']*6
+offset = [1]*3 + [0]*3 
+max_forward = 600
+Δmax_forward = 75
+n_freq = 1
+"""
+
+"""
+#CS3: IP
+j0 = 0.0328706
+n_comp = 4
+path = ['TC3/LBFGS_LS_SST5/','TC3/LBFGS_LS_SST2/','TC3/LBFGS_LS_SST6/','TC3/LBFGS_LS_SST3/']
+filename = ['cross_inversion_sst5_lbfgs_ls','cross_inversion_sst2_lbfgs_ls','cross_inversion_sst6_lbfgs_ls','cross_inversion_sst3_lbfgs_ls']
+lab = [r'$\beta = 2$',r'$\beta = 3$',r'$\beta = 4$',r'$\beta = 5$']
+glob = ['linesearch'] * 4
+order = [1] * 4
+col = ['k', '#696969', '#A9A9A9', '#D3D3D3']
+mrk = [''] * 4
+lin = ['-'] * 4
+max_forward = 900
+Δmax_forward = 80
+n_freq = 1
+"""
+
+#CS3: L-BFGS
+j0 = 0.0328706
+n_comp = 3
+path = ['TC3/LBFGS_LS_SST2/','TC3/LBFGS_TRB_SST2/','TC3/LBFGS_TRBR_SST2/']
+filename = ['cross_inversion_sst2_lbfgs_ls','cross_inversion_sst2_lbfgs_trB','cross_inversion_sst2_lbfgs_trBR']
+lab = ['LS','TR-P(B)','TR-R(B)']
+glob = ['linesearch'] + ['trust-region'] * 1 + ['trust-region-retro'] * 1
+order = [1.5] * 3
+col = ['k'] + ['C1'] + ['C4']
+mrk = ['']*3
+lin = ['-']*3
+max_forward = 900
+Δmax_forward = 80
+offset = [1]*6
+n_freq = 1
+
+
+"""
+#CS3: Newton
+j0 = 0.0328706
+n_comp = 6
+path = ['TC3/NCG_SST2/','TC3/NSCGB_SST2/','TC3/NSCGBR_SST2/','TC3/NCGGN_SST2/','TC3/NSCGBGN_SST2/','TC3/NSCGBRGN_SST2/']
+filename = ['cross_inversion_sst2_ncg','cross_inversion_sst2_nscgB','cross_inversion_sst2_nscgBR','cross_inversion_sst2_ncggn','cross_inversion_sst2_nscgBgn','cross_inversion_sst2_nscgBRgn']
+lab = ['FN-LS','FN-TR-P(B)','FN-TR-R(B)','GN-LS','GN-TR-P(B)','GN-TR-R(B)']
+glob = ['linesearch','trust-region','trust-region-retro'] * 2
+order = [2.] * 6
+col = ['k','C1','C4']*2
+mrk = ['o']*6
+lin = ['-']*6 
+max_forward = 900
+Δmax_forward = 80
+offset = [1]*6
+n_freq=1
+"""
+
+"""
+LOAD
+"""
+jcfs = []
+p_ratiocfs = []
+r_ratiocfs = []
+n_itcf = []
+n_scalecf = []
+n_itcfs = []
+
+for i in range(0,n_comp):
+    n_scalef = [[]]*n_freq; 
+    jfs = [[]]*n_freq;
+    p_ratiofs = [[]]*n_freq;
+    r_ratiofs = [[]]*n_freq;
+    n_itfs = [[]]*n_freq;
+    n_itf = [0]*n_freq
+    scalefs = [[]]*n_freq;
+    for f in range(0,n_freq):
+        n_scalef[f] = np.loadtxt(path[i]+filename[i]+'_globalminimum_historyg'+str(f)+'.csv',delimiter=';',max_rows=1,skiprows=3,dtype=np.integer)[1];
+        scalefs[f] = np.loadtxt(path[i]+filename[i]+'_globalminimum_historyg'+str(f)+'.csv',delimiter=';',usecols=2,max_rows=np.int(n_scalef[f]),skiprows=6,dtype=np.float);
+        jfs[f] = [[]]*n_scalef[f];
+        p_ratiofs[f] = [[]]*n_scalef[f];
+        r_ratiofs[f] = [[]]*n_scalef[f];
+        n_itfs[f] = np.zeros(n_scalef[f],np.integer);
+        line = 3 
+        for s in range(0,n_scalef[f]):        
+            n_itfs[f][s] = np.loadtxt(path[i]+filename[i]+'_localminimum_historyg'+str(f)+'.csv',delimiter=';',usecols = 1, max_rows=1,skiprows=line,dtype=np.integer);
+            n_itf[f] = n_itf[f]+n_itfs[f][s]
+            line=line+2
+            jfs[f][s] = np.loadtxt(path[i]+filename[i]+'_localminimum_historyg'+str(f)+'.csv',delimiter=';',max_rows=np.int(n_itfs[f][s]+1),skiprows=line,dtype=np.float)[:,0] / j0;
+            if glob[i] == 'trust-region' or glob[i] == 'trust-region-retro':
+                p_ratiofs[f][s] = np.loadtxt(path[i]+filename[i]+'_localminimum_historyg'+str(f)+'.csv',delimiter=';',max_rows=np.int(n_itfs[f][s]+1),skiprows=line,dtype=np.float)[:,5+offset[i]];
+                r_ratiofs[f][s] = np.loadtxt(path[i]+filename[i]+'_localminimum_historyg'+str(f)+'.csv',delimiter=';',max_rows=np.int(n_itfs[f][s]+1),skiprows=line,dtype=np.float)[:,6+offset[i]];
+                r_ratiofs[f][s][p_ratiofs[f][s] < 0] = np.nan;
+            line=line+n_itfs[f][s]+1+1
+    jcfs.append(jfs)
+    p_ratiocfs.append(p_ratiofs)
+    r_ratiocfs.append(r_ratiofs)
+    n_itcf.append(n_itf)
+    n_itcfs.append(n_itfs)
+    n_scalecf.append(n_scalef)
+
+it_linecf = [[]*n_freq]*n_comp;
+it_descentcf = [[]*n_freq]*n_comp;
+bnd_descentcf = [[]*n_freq]*n_comp; 
+n_forwardcf = [[]*n_freq]*n_comp;
+successcf = [[]*n_freq]*n_comp;
+negativecf = [[]*n_freq]*n_comp;
+etacf = [[]*n_freq]*n_comp;
+for i in range(0,n_comp):
+    it_linef = [[]]*n_freq
+    it_descentf = [[]]*n_freq
+    bnd_descentf = [[]]*n_freq
+    n_forwardf = [[]]*n_freq
+    successf = [[]]*n_freq
+    negativef = [[]]*n_freq
+    etaf = [[]]*n_freq
+    for f in range(0,n_freq):
+        n_forwardf[f] = np.zeros(n_itcf[i][f],np.integer)
+
+        if glob[i] == 'linesearch': 
+            it_linef[f] = np.loadtxt(path[i]+filename[i]+'_linesearch_historyg'+str(f)+'.csv',delimiter=';',usecols = 5, skiprows=1,dtype=np.integer)
+            n_forwardf[f] += it_linef[f] 
+
+        if order[i] == 1.5 and (glob[i] == 'trust-region' or glob[i] == 'trust-region-retro') :
+            successf[f] = np.zeros(n_itcf[i][f],np.integer)
+            bnd_descentf[f] = np.zeros(n_itcf[i][f],np.integer)
+            line = 2
+            for n in range(0,n_itcf[i][f]):
+                bnd_descentf[f][n] = np.loadtxt(path[i]+filename[i]+'_descent_historyg'+str(f)+'.csv',delimiter=';',usecols=4,max_rows=1,skiprows=line,dtype=np.integer);            
+                successf[f][n] = np.loadtxt(path[i]+filename[i]+'_descent_historyg'+str(f)+'.csv',delimiter=';',usecols=0,max_rows=1,skiprows=line,dtype=np.integer);
+                line=line+2        
+        if order[i] == 2:
+            it_descentf[f] = np.zeros(n_itcf[i][f],np.integer)
+            successf[f] = np.zeros(n_itcf[i][f],np.integer)
+            etaf[f] = np.zeros(n_itcf[i][f],np.float)
+            negativef[f] = np.zeros(n_itcf[i][f],np.integer)
+            bnd_descentf[f] = np.zeros(n_itcf[i][f],np.integer)
+            line = 3 
+            for n in range(0,n_itcf[i][f]):        
+                it_descentf[f][n] = np.loadtxt(path[i]+filename[i]+'_descent_historyg'+str(f)+'.csv',delimiter=';',usecols=3,max_rows=1,skiprows=line,dtype=np.integer);
+                if glob[i] == 'trust-region' or glob[i] == 'trust-region-retro':
+                    bnd_descentf[f][n] = np.loadtxt(path[i]+filename[i]+'_descent_historyg'+str(f)+'.csv',delimiter=';',usecols=6,max_rows=1,skiprows=line,dtype=np.integer);
+                etaf[f][n] = np.loadtxt(path[i]+filename[i]+'_descent_historyg'+str(f)+'.csv',delimiter=';',usecols=1,max_rows=1,skiprows=line,dtype=np.float);            
+                successf[f][n] = np.loadtxt(path[i]+filename[i]+'_descent_historyg'+str(f)+'.csv',delimiter=';',usecols=0,max_rows=1,skiprows=line,dtype=np.integer);
+                line=line+it_descentf[f][n]+1
+                negativef[f][n] = (np.loadtxt(path[i]+filename[i]+'_descent_historyg'+str(f)+'.csv',delimiter=';',usecols=1,max_rows=1,skiprows=line,dtype=np.float)<0);
+                line=line+2
+            n_forwardf[f] += 2*it_descentf[f]
+            if glob[i] == 'trust-region-retro':
+                n_forwardf[f] += 2 #retrospective Hessian-vector product
+        n_forwardf[f] += 2
+        n_forwardf[f] = np.cumsum(n_forwardf[f])
+        n_forwardf[f] = np.insert(n_forwardf[f],0,0)
+    n_forwardcf[i] = n_forwardf    
+    it_descentcf[i] = it_descentf
+    bnd_descentcf[i] = bnd_descentf
+    it_linecf[i] = it_linef
+    successcf[i] = successf
+    etacf[i] = etaf
+    negativecf[i] = negativef
+
+"""
+DATA CONVERGENCE PLOT
+"""
+fig = [[]]*n_freq
+for f in range(0,n_freq):
+    fig[f]=plt.figure(figsize=(2*L,H),tight_layout=False);
+    plt.subplots_adjust(top=0.95,right=0.92,bottom=0.13,left=0.08)
+    plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+    plt.grid(zorder=0)
+    plt.xlim([0,max_forward])
+    plt.xticks(np.arange(0,max_forward+1,Δmax_forward))
+    #plt.ylim([5e-4,2])
+    plt.ylim([5e-4,2])
+    plt.yscale('log')
+    for i in (0,1,2):#range(0,n_comp):
+        n_it=0
+        lines = [[]]*n_scalecf[i][f]
+        plt.plot(n_forwardcf[i][f][n_it],jcfs[i][f][0][0],color=col[i],marker='o',zorder=3,markersize=3)
+        for s in range(0,n_scalecf[i][f]):
+            lines[s] = plt.plot(n_forwardcf[i][f][range(n_it,n_it+n_itcfs[i][f][s]+1)],jcfs[i][f][s][0:n_itcfs[i][f][s]+1],color=col[i],linestyle=lin[i],zorder=2,label=lab[i],marker=mrk[i],markersize=5.,linewidth=lw)
+            plt.plot(n_forwardcf[i][f][n_it+n_itcfs[i][f][s]],jcfs[i][f][s][n_itcfs[i][f][s]],color=col[i],marker='x',zorder=3,markersize=3.)
+            n_it=n_it+n_itcfs[i][f][s]
+        n_it = 0
+    #plt.legend()
+    #plt.savefig('convergence_comp'+str(f)+'.eps')
+
+"""
+HISTORY CONVERGENCE PLOT (Retrospective/prospective ratio)
+"""
+fig = [[]]*n_freq
+for f in range(0,n_freq):
+    fig[f]=plt.figure(figsize=(2*L,H / 2.),tight_layout=False);
+    plt.subplots_adjust(top=0.95,right=0.92,bottom=0.23,left=0.08)
+    plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+    plt.grid(zorder=0)
+    plt.yticks(np.linspace(0,1.5,4))
+    plt.ylim([0.0,1.75])
+    
+    plt.xticks(np.arange(0,23,2))
+    plt.xlim([1,23.5])
+    
+    #plt.xticks(np.arange(0,46,5))
+    #plt.xlim([1,45.5])
+    
+    for i in (1,2):
+        if glob[i] == 'trust-region':
+            plt.plot(range(1,n_itcfs[i][f][0]+1),p_ratiocfs[i][f][0][1:],color=col[i],linestyle=lin[i],zorder=2,label=r'$\rho_{\text{p}}$',linewidth=lw)
+        if glob[i] == 'trust-region-retro':
+            plt.plot(range(1,n_itcfs[i][f][0]+1),r_ratiocfs[i][f][0][1:],color=col[i],linestyle=lin[i],zorder=2,label=r'$\rho_{\text{r}}$',linewidth=lw)
+        #plt.legend()
+    #plt.savefig('ratio_nscgrgn.eps')
+
+"""
+HISTORY CONVERGENCE PLOT (forcing sequence)
+"""
+fig = [[]]*n_freq
+for f in range(0,n_freq):
+    fig[f]=plt.figure(figsize=(2*L,H / 2.),tight_layout=False);
+    plt.subplots_adjust(top=0.93,right=0.92,bottom=0.23,left=0.08)
+    plt.grid(zorder=0)
+    plt.ylim([0.3,1.])
+    plt.yticks(np.arange(0.4,1.01,0.2))
+    
+    #plt.xticks(np.arange(0,23,2))
+    #plt.xlim([0,22.5])
+    
+    plt.xticks(np.arange(0,46,5))
+    plt.xlim([0,45.5])
+    
+    for i in (3,):#range(0,n_comp,1):
+        plt.plot(range(0,n_itcfs[i][f][0]),np.sqrt(etacf[i][f][:]),color=col[i],linestyle=lin[i],zorder=2,label=r'$\rho_{\text{p}}$',linewidth=lw)
+        #plt.legend()
+    #plt.savefig('eta_comp.eps')
+    
+"""
+HISTORY CONVERGENCE PLOT (inner iteration per outer iteration)
+"""
+fig = [[]]*n_freq
+for f in range(0,n_freq):
+    fig[f]=plt.figure(figsize=(2*L,H / 2.),tight_layout=False);
+    plt.subplots_adjust(top=0.93,right=0.92,bottom=0.23,left=0.08)
+    plt.ticklabel_format(style='plain', axis='y', scilimits=(0,0))
+    plt.grid(zorder=0)
+    
+    #plt.xticks(np.arange(0,23,2))
+    #plt.yticks(np.arange(0,15,4))
+    #plt.xlim([0,22.5])
+    #plt.ylim([0,14.5])
+    
+    plt.xticks(np.arange(0,46,5))
+    plt.yticks(np.arange(0,41,10))
+    plt.xlim([0,45.5])
+    plt.ylim([0,40.])
+    
+    #plt.xticks(np.arange(0,23,2))
+    #plt.yticks(np.arange(0,101,25))
+    #plt.xlim([0,22.5])
+    #plt.ylim([0,100])
+    
+    for i in (0,1,2): #range(0,n_comp):
+        plt.plot(range(0,n_itcfs[i][f][0]),it_descentcf[i][f][:],color=col[i],linestyle=lin[i],zorder=2,linewidth=lw)
+    #plt.savefig('it_comp.eps')
+    
+"""
+MODEL CONVERGENCE PLOT
+"""
+#To investigate !!!!
+errorcfs = []
+for i in range(0,n_comp):
+    errorc = np.loadtxt(path[i]+filename[i]+'_error.csv',delimiter=';')[:,0];        
+    errorcfs.append(errorc)
+    #plt.plot(n_forwardcf[i][0],errorc,color=col[i],linestyle=lin[i],zorder=2,label=lab[i])
+    plt.plot(jcfs[i][0][0],errorc,color=col[i],linestyle=lin[i],zorder=2,label=lab[i])
+    plt.legend()
+plt.xscale('log')
+plt.gca().invert_xaxis()
+        
+"""
+STATS
+"""
+print(''.rjust(10), 'Wave (tot.)'.rjust(20), 'Outer it. (tot.)'.rjust(20), 'Inner it. (mean)'.rjust(20), 'Failure (tot.)'.rjust(20), 'Bnd reached (tot.)'.rjust(20), 'Success (tot.)'.rjust(20), 'Negative curv. (tot.)'.rjust(20))
+for i in range(0,n_comp):
+        outter = np.sum(n_itcf[i])
+        inner = 0
+        failure = 0
+        boundary = 0
+        success = 0
+        negative = 0
+        call = 0
+        for f in range(0,n_freq):
+            call += n_forwardcf[i][f][-1]
+            inner += np.sum(it_descentcf[i][f])
+            if glob[i] == 'linesearch':
+                failure += np.count_nonzero(it_linecf[i][f])
+            elif glob[i] == 'trust-region' or glob[i] == 'trust-region-retro':
+                for s in range(0,n_scalecf[i][f]):
+                    failure += n_itcfs[i][f][s] - np.count_nonzero(np.diff(jcfs[i][f][s]))
+                    boundary += np.sum(bnd_descentcf[i][f])
+            if order[i] == 1.5 and (glob[i] == 'trust-region' or glob[i] == 'trust-region-retro') :
+                success += np.sum(successcf[i][f])
+            if order[i] == 2:
+                success += np.sum(successcf[i][f])
+                negative += np.sum(negativecf[i][f])
+        print(lab[i].rjust(10), str(call).rjust(20), str(outter).rjust(20), str(inner/outter).rjust(20), str(failure/outter).rjust(20), str(boundary/outter).rjust(20), str(success/outter).rjust(20), str(negative/outter).rjust(20) )