updates for GMD paper

aux/plotting/plot_inv_results.py

@@ -39,7 +39,7 @@
 # Parameters:
 
 # Simulation Directory
-run_name = 'ismipc_rc_1e6'
+run_name = 'ismipc_1000_corr'
 dd = Path(os.environ['FENICS_ICE_BASE_DIR']) / 'example_cases' / run_name
 
 # Output Directory
@@ -167,5 +167,5 @@
 cbar.ax.set_xlabel(r'$U-U_{obs}$ (m $yr^{-1}$)')
 
 plt.tight_layout(2.0)
-plt.savefig(os.path.join(outdir, 'inv_results.pdf'))
+#plt.savefig(os.path.join(outdir, 'inv_results.pdf'))
 plt.show()

aux/plotting/plot_inv_results_hist.py

@@ -22,6 +22,7 @@
 import matplotlib.ticker as ticker
 import seaborn as sns
 import os
+import math
 from fenics import *
 from tlm_adjoint_fenics import *
 from fenics_ice import model, config, inout
@@ -41,7 +42,7 @@
 # Parameters:
 
 # Simulation Directory
-run_name = 'ismipc_rc_1e6'
+run_name = 'ismipc_4000_corr'
 if (len(sys.argv)>1):
     run_name = sys.argv[1]
 
@@ -105,12 +106,14 @@
 uv_obs_file = str(next(results_dir.glob("*uv_obs.xml")))
 alpha_sigma_file = str(next(results_dir.glob("*alpha_sigma.xml")))
 dQ_file = str(next(results_dir.glob("*dQ_ts.xml")))
+thick_file = str(next(results_dir.glob("*H_fwd.xml")))
 
 U = Function(V, U_file)
 alpha = Function(Qp, alpha_file)
 uv_obs = Function(M, uv_obs_file)
 alpha_sigma = Function(Qp, alpha_sigma_file)
 dQ = Function(Qp,dQ_file)
+H_fwd = Function(Q,thick_file)
 # B2 = Function(M, os.path.join(dd,'B2.xml'))
 
 u, v = U.split()
@@ -167,27 +170,31 @@
 v    = B2.compute_vertex_values(mesh)
 minv = np.min(v)
 maxv = np.max(v)
-levels = np.linspace(minv,maxv,numlev)
-ticks = np.linspace(minv,maxv,3)
+levels = np.linspace(0,2000,numlev)
+ticks = np.linspace(0,2000,3)
 ax.tick_params(**tick_options)
 ax.text(0.05, 0.95, 'a', transform=ax.transAxes,
     fontsize=13, fontweight='bold', va='top')
-c = ax.tricontourf(x, y, t, v, levels = levels, cmap=plt.get_cmap(cmap))
+c = ax.tricontourf(x, y, t, v, levels = levels, cmap=plt.get_cmap(cmap),extend='both')
 cbar = plt.colorbar(c, ticks=ticks, pad=0.05, orientation="horizontal",shrink=.75)
-cbar.ax.set_xlabel(r'${\beta^2}$ (Pa $m^{-1}$ yr)')
+cbar.ax.set_xlabel(r'${C^2}$ (Pa $m^{-1}$ yr)')
 
 ax  = fig.add_subplot(232)
 ax.set_aspect('equal')
 v    = alpha_sigma.compute_vertex_values(mesh)
 minv = np.min(v)
 maxv = np.max(v)
+minv = 0
+maxv = math.ceil(np.max(v)/10.)*10.
+#levels = np.linspace(0,50,numlev)
+#ticks = np.linspace(0,50,3)
 levels = np.linspace(minv,maxv,numlev)
 ticks = np.linspace(minv,maxv,3)
 ax.tick_params(**tick_options)
 ax.text(0.05, 0.95, 'b', transform=ax.transAxes,
     fontsize=13, fontweight='bold', va='top')
 c = ax.tricontourf(x, y, v, levels = levels, cmap=plt.get_cmap(cmap))
-cbar = plt.colorbar(c, ticks=ticks, pad=0.05, orientation="horizontal", format=ticker.FormatStrFormatter('%1.1e'),shrink=.75)
+cbar = plt.colorbar(c, ticks=ticks, pad=0.15, orientation="horizontal", format=ticker.FormatStrFormatter('%1.1e'),shrink=.75)
 cbar.ax.xaxis.set_major_locator(ticker.LinearLocator(3))
 
 cbar.ax.set_xlabel(r'$\sigma$ ($Pa^{0.5}$ $m^{-0.5}$ $yr^{0.5}$)')
@@ -207,6 +214,8 @@
 
 ax  = fig.add_subplot(234)
 ax.hist(tot_err.vector()[:])
+ax.text(0.05, 0.95, 'd', transform=ax.transAxes,
+    fontsize=13, fontweight='bold', va='top')
 plt.xlabel('velocity error (m yr$^{-1}$)')
 plt.ylabel('count')
 #v   = uv_obs.compute_vertex_values(mesh)
@@ -236,14 +245,30 @@
 ax.set_aspect('equal')
 v   = dQ.compute_vertex_values(mesh)
 max_diff = np.rint(np.max(np.abs(v)))
+max_diff=5.0e8
 levels = np.linspace(-max_diff,max_diff,numlev)
 ticks = np.linspace(-max_diff,max_diff,3)
 ax.tick_params(**tick_options)
 ax.text(0.05, 0.95, 'e', transform=ax.transAxes,
     fontsize=13, fontweight='bold', va='top')
 c = ax.tricontourf(x, y, t, v, levels = levels, cmap=plt.get_cmap(cmap_div))
 cbar = plt.colorbar(c, ticks=ticks, pad=0.05, orientation="horizontal",shrink=.75)
-cbar.ax.set_xlabel(r'$\partial Q/\partial \beta$')
+cbar.ax.set_xlabel(r'$\partial Q/\partial C$ ($m^4 Pa^{-1} m a^{-1}$)',size=10)
+
+ax  = fig.add_subplot(236)
+ax.set_aspect('equal')
+v   = H_fwd.compute_vertex_values(mesh)
+max_diff = np.rint(np.max(np.abs(1000.-v)))
+max_diff=10.
+levels = np.linspace(1000.-max_diff,1000.+max_diff,numlev)
+ticks = np.linspace(1000.-max_diff,1000.+max_diff,3)
+ax.tick_params(**tick_options)
+ax.text(0.05, 0.95, 'f', transform=ax.transAxes,
+    fontsize=13, fontweight='bold', va='top')
+c = ax.tricontourf(x, y, t, v, levels = levels, cmap=plt.get_cmap(cmap_div),extend='both')
+cbar = plt.colorbar(c, ticks=ticks, pad=0.05, orientation="horizontal",shrink=.75)
+cbar.ax.set_xlabel(r'$H(T=30)$ $(m)$')
 
 plt.tight_layout(2.0)
 plt.savefig(os.path.join(outdir, run_name + '_inv_results.png'))
+

aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py

@@ -54,6 +54,8 @@
     lpos = np.argwhere(lam > 0)
     lneg = np.argwhere(lam < 0)
     lind = np.arange(0,len(lam))
+#    plt.semilogy(lind[lpos], lam[lpos], '.', alpha = 0.5, mew=0, label =labels[i])
+#    plt.semilogy(lind[lneg], np.abs(lam[lneg]), '.k', alpha = 0.12, mew=0,)
     plt.semilogy(lind[lpos], 1./(1+lam[lpos]), '.', alpha = 0.5, mew=0, label =labels[i])
     plt.semilogy(lind[lneg], 1./(1+np.abs(lam[lneg])), '.k', alpha = 0.12, mew=0,)
 
@@ -62,3 +64,4 @@
 plt.xlabel('Eigenvalue index')
 plt.ylabel('Uncertainty reduction')
 plt.savefig(os.path.join(outdir,figname), bbox_inches="tight")
+#plt.show()

aux/plotting/plotting_new/plot_paths_sample.py

@@ -57,7 +57,8 @@
 
     sigma_interp = np.interp(dQ_t, sigma_t, sigma_vals)
     sigma_prior_interp = np.interp(dQ_t, sigma_t, sigma_prior_vals)
-
+
+
     QoIs_MC = np.load(result_dir/"sampling_results.npy")
     QoIs_MC = QoIs_MC[:,QoIs_MC[-1,:]>0]
     print( str(np.shape(QoIs_MC)) + ' values OK')
@@ -74,19 +75,24 @@
     upperlim = max(upperlim,np.max(s))
     upperlim = max(upperlim,np.max(sm))
 
+    if(i==0):
+     p1 = axarr[0].plot(x, y, 'b',label='$\gamma$=50')
+     p2 = axarr[0].fill_between(x, y-s, y+s, facecolor='b', alpha = .6)
+    else:
+     p3 = axarr[0].plot(x, y, 'g',label='$\gamma$=1')
+     p4 = axarr[0].fill_between(x, y-s, y+s, facecolor='g', alpha = .6)
 
-    axarr[0].plot(x, y, 'k' if i == 0 else 'k:')
 #    axarr[0].fill_between(x, y-sp, y+sp, facecolor='b')
-    axarr[0].fill_between(x, y-s, y+s, facecolor='b' if i == 0 else 'g', alpha = .6 if i==0 else .6)
 #    axarr[0].fill_between(x, y-sm, y+sm, facecolor='r' if i == 0 else 'y', alpha = .5 if i==0 else .6)
 
-    axarr[1].semilogy(x, sp, 'b' if i == 0 else 'g')
-    axarr[1].semilogy(x, s, 'b:' if i == 0 else 'g:')
+    axarr[1].semilogy(x, sp, 'b:' if i == 0 else 'g:')
+    axarr[1].semilogy(x, s, 'b' if i == 0 else 'g')
     if (i==0):
      clr = 'blue'
     else:
      clr = 'green'
     axarr[1].semilogy(sigma_t, sigma_mc_vals,color=clr,marker='+',linestyle='none',markersize=6)
+axarr[0].legend([(p1[0],p2),(p3[0],p4)],['$\gamma$=50','$\gamma$=1'])
 
 axarr[0].set_xlabel('Time (yrs)')
 axarr[1].set_xlabel('Time (yrs)')

example_cases/ismipc_1000/ismipc_1000.sh

@@ -14,8 +14,8 @@ cp output_momsolve/ismipc_U_obs.h5 input/
 #Run each phase of the model in turn
 RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_1000.toml
-#python $RUN_DIR/run_forward.py ismipc_1000.toml
-#python $RUN_DIR/run_eigendec.py ismipc_1000.toml
-#python $RUN_DIR/run_errorprop.py ismipc_1000.toml
-#python $RUN_DIR/run_invsigma.py ismipc_1000.toml
-#python $RUN_DIR/run_sample_post.py ismipc_1000.toml
+python $RUN_DIR/run_forward.py ismipc_1000.toml
+python $RUN_DIR/run_eigendec.py ismipc_1000.toml
+python $RUN_DIR/run_errorprop.py ismipc_1000.toml
+python $RUN_DIR/run_invsigma.py ismipc_1000.toml
+python $RUN_DIR/run_sample_post.py ismipc_1000.toml

example_cases/ismipc_1000_corr/ismipc_1000_corr.sh

@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=1000. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=1000. -t "ismipc_1000.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=1000. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=1000. -t "ismipc_1000.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 
@@ -18,6 +18,6 @@ RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_1000_corr.toml
 python $RUN_DIR/run_forward.py ismipc_1000_corr.toml
 python $RUN_DIR/run_eigendec.py ismipc_1000_corr.toml
-#python $RUN_DIR/run_errorprop.py ismipc_1000_corr.toml
+python $RUN_DIR/run_errorprop.py ismipc_1000_corr.toml
 python $RUN_DIR/run_invsigma.py ismipc_1000_corr.toml
-#python $RUN_DIR/run_sample_post.py ismipc_1000_corr.toml
+python $RUN_DIR/run_sample_post.py ismipc_1000_corr.toml

example_cases/ismipc_2000/ismipc_2000.sh

@@ -18,4 +18,4 @@ python $RUN_DIR/run_forward.py ismipc_2000.toml
 python $RUN_DIR/run_eigendec.py ismipc_2000.toml
 python $RUN_DIR/run_errorprop.py ismipc_2000.toml
 python $RUN_DIR/run_invsigma.py ismipc_2000.toml
-#python $RUN_DIR/run_sample_post.py ismipc_2000.toml
+python $RUN_DIR/run_sample_post.py ismipc_2000.toml

example_cases/ismipc_2000_corr/ismipc_2000_corr.sh

@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=2000. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=2000. -t "ismipc_2000.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=2000. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=2000. -t "ismipc_2000.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 

example_cases/ismipc_4000/ismipc_4000.sh

@@ -18,4 +18,4 @@ python $RUN_DIR/run_forward.py ismipc_4000.toml
 python $RUN_DIR/run_eigendec.py ismipc_4000.toml
 python $RUN_DIR/run_errorprop.py ismipc_4000.toml
 python $RUN_DIR/run_invsigma.py ismipc_4000.toml
-#python $RUN_DIR/run_sample_post.py ismipc_4000.toml
+python $RUN_DIR/run_sample_post.py ismipc_4000.toml

example_cases/ismipc_4000_corr/ismipc_4000_corr.sh

@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=4000. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=4000. -t "ismipc_4000.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=4000. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=4000. -t "ismipc_4000.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 
@@ -18,6 +18,6 @@ RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_4000_corr.toml
 python $RUN_DIR/run_forward.py ismipc_4000_corr.toml
 python $RUN_DIR/run_eigendec.py ismipc_4000_corr.toml
-#python $RUN_DIR/run_errorprop.py ismipc_4000_corr.toml
+python $RUN_DIR/run_errorprop.py ismipc_4000_corr.toml
 python $RUN_DIR/run_invsigma.py ismipc_4000_corr.toml
-#python $RUN_DIR/run_sample_post.py ismipc_4000_corr.toml
+python $RUN_DIR/run_sample_post.py ismipc_4000_corr.toml

example_cases/ismipc_500/ismipc_500.sh

@@ -18,4 +18,4 @@ python $RUN_DIR/run_forward.py ismipc_500.toml
 python $RUN_DIR/run_eigendec.py ismipc_500.toml
 python $RUN_DIR/run_errorprop.py ismipc_500.toml
 python $RUN_DIR/run_invsigma.py ismipc_500.toml
-#python $RUN_DIR/run_sample_post.py ismipc_500.toml
+python $RUN_DIR/run_sample_post.py ismipc_500.toml

example_cases/ismipc_500_corr/ismipc_500_corr.sh

@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=500. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=500. -t "ismipc_500.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=500. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=500. -t "ismipc_500.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 
@@ -18,6 +18,6 @@ RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_500_corr.toml
 python $RUN_DIR/run_forward.py ismipc_500_corr.toml
 python $RUN_DIR/run_eigendec.py ismipc_500_corr.toml
-#python $RUN_DIR/run_errorprop.py ismipc_500_corr.toml
+python $RUN_DIR/run_errorprop.py ismipc_500_corr.toml
 python $RUN_DIR/run_invsigma.py ismipc_500_corr.toml
-#python $RUN_DIR/run_sample_post.py ismipc_500_corr.toml
+python $RUN_DIR/run_sample_post.py ismipc_500_corr.toml

example_cases/ismipc_8000/ismipc_8000.sh

@@ -18,4 +18,4 @@ python $RUN_DIR/run_forward.py ismipc_8000.toml
 python $RUN_DIR/run_eigendec.py ismipc_8000.toml
 python $RUN_DIR/run_errorprop.py ismipc_8000.toml
 python $RUN_DIR/run_invsigma.py ismipc_8000.toml
-#python $RUN_DIR/run_sample_post.py ismipc_8000.toml
+python $RUN_DIR/run_sample_post.py ismipc_8000.toml

example_cases/ismipc_8000_corr/ismipc_8000_corr.sh

@@ -6,9 +6,9 @@ python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/momsolve_mesh.xml -x
 python $FENICS_ICE_BASE_DIR/aux/gen_rect_mesh.py -o ./input/ismip_mesh.xml -xmax 40000 -ymax 40000 -nx 30 -ny 30
 python $FENICS_ICE_BASE_DIR/aux/gen_ismipC_domain.py -o ./input/ismipc_input.h5 -L 40000 -nx 100 -ny 100 --reflect
 python $FENICS_ICE_BASE_DIR/runs/run_momsolve.py momsolve.toml
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=8000. 
-cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
-python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve_cov.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=8000. -t "ismipc_8000.toml"
+#python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs_nocorr.h5" -L 40000 -d ./output_momsolve --ls=8000. 
+#cp output_momsolve/ismipc_U_obs_nocorr.h5 input/
+python $FENICS_ICE_BASE_DIR/aux/Uobs_from_momsolve.py -i "U.h5" -o "ismipc_U_obs.h5" -L 40000 -d ./output_momsolve --ls=8000. -t "ismipc_8000.toml"
 cp output_momsolve/ismipc_U_obs.h5 input/
 
 
@@ -18,6 +18,6 @@ RUN_DIR=$FENICS_ICE_BASE_DIR/runs/
 python $RUN_DIR/run_inv.py ismipc_8000_corr.toml
 python $RUN_DIR/run_forward.py ismipc_8000_corr.toml
 python $RUN_DIR/run_eigendec.py ismipc_8000_corr.toml
-#python $RUN_DIR/run_errorprop.py ismipc_8000_corr.toml
+python $RUN_DIR/run_errorprop.py ismipc_8000_corr.toml
 python $RUN_DIR/run_invsigma.py ismipc_8000_corr.toml
-#python $RUN_DIR/run_sample_post.py ismipc_8000_corr.toml
+python $RUN_DIR/run_sample_post.py ismipc_8000_corr.toml

example_cases/plots_gmd/make_plots.py

@@ -0,0 +1,59 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+
+# Lcurve
+
+dir='../ismipc_rc_1e6/output/'
+
+gam = np.load(dir + 'LCurvGam.npy')
+norm = np.load(dir + 'LCurvNorm.npy')
+mis = np.load(dir + 'LCurvMis.npy')
+
+fig,ax = plt.subplots()
+
+for i in range(len(gam)):
+  ax.plot(norm[i],mis[i],'o',markersize=4,markerfacecolor='blue',markeredgecolor='k')
+
+for i in range(len(gam)):
+  if (i==4):
+   ax.annotate(str(gam[i]),(norm[i]-6,mis[i]-6))
+  else:
+   ax.annotate(str(gam[i]),(norm[i]+1,mis[i]+1))
+
+plt.xlabel('Regularisation cost')
+plt.ylabel('Model-data misfit cost')
+
+
+plt.savefig('L_curve.png')
+
+plt.close(fig)
+
+pwd = os.getcwd()
+
+os.system('python ../../aux/plotting/plot_inv_results_hist.py ismipc_rc_1e6 ' + pwd)
+os.system('python ../../aux/plotting/plot_inv_results_hist.py ismipc_rc_1e4 ' + pwd)
+
+plt.close(fig)
+
+os.system('python ../../aux/plotting/plot_leading_eigenfuncs.py ' + pwd)
+
+plt.close(fig)
+
+os.system("python ../../aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py $PWD reg_decay.png ismipc_rc_1e6 '$\gamma_{alpha}$=50' ismipc_rc_1e5 '$\gamma_{alpha}$=10' ismipc_rc_1e4 '$\gamma_{alpha}$=1'")
+
+plt.close(fig)
+
+os.system("python ../../aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py $PWD res_decay.png ismipc_20x20 '$\Delta$ x = 2 km' ismipc_rc_1e5 '$\Delta$ x = 1.33 km' ismipc_40x40 '$\Delta$ x = 1 km'")
+
+plt.close(fig)
+
+os.system("python ../../aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py $PWD sample_decay.png ismipc_8000 8000m ismipc_4000 4km ismipc_2000 2km ismipc_1000 1km ismipc_500 500m")
+
+plt.close(fig)
+
+#os.system("python ../../aux/plotting/plotting_new/plot_eigenvalue_decay_gen_new.py $PWD sample_decay_autocorr_alt2.png corr_alt/ismipc_8000_corr 8km corr_alt/ismipc_4000_corr 4km corr_alt/ismipc_2000_corr 2km corr_alt/ismipc_1000_corr 1km corr_alt/ismipc_500_corr 500m")
+
+#plt.close(fig)
+
+os.system('python ../../aux/plotting/plotting_new/plot_paths_sample.py $PWD paths.png')

example_cases/run_all.sh

@@ -1,7 +1,7 @@
 #!/bin/bash
 #Runs each model in turn (after first clearing out any previous results)
 
-for f in $(ls -d ismipc*00 ismipc_rc*);
+for f in $(ls -d ismipc*corr);
 do echo $f
 cd $f
 echo "Running model in $f"