Thesis Part 3 Test case ok

discrepancy.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri May  6 12:42:21 2022
+
+@author: xavier
+"""
+import numpy as np
+import scipy as sp
+from scipy import linalg
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Times New Roman'],'size':15})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 3.;
+W = 8.;
+
+#
+# Time domain
+#
+
+"""
+path = "extended_marmousi_data/"
+prename = "paper2_extended_synthetics_data"
+snr_tag = 'snr2'
+Nf = 3
+"""
+
+"""
+path = "concrete3_data/"
+prename = "paper2_concrete3_synthetics_data"
+snr_tag = 'snrm3'
+Nf=9
+"""
+
+path = "cross_data/"
+prename = "cross_synthetics_data"
+snr_tag = 'snr6'
+Nf=7
+
+noise = 0.
+for f in range(0,Nf):
+    tmp = np.loadtxt(path+prename+str(f)+".csv",delimiter=";");
+    tmp_snr = np.loadtxt(path+prename+'_'+snr_tag+'f'+str(f)+".csv",delimiter=";");
+    noise += sp.linalg.norm(tmp-tmp_snr)**2
\ No newline at end of file

error.cpp

@@ -14,7 +14,7 @@ using namespace gmshfem::post;
 #include "specific/configuration/newConfiguration.h"
 #include "specific/model/parametrization/newParametrization.h"
 
-enum class Error: unsigned int { Abs=0, Sqd=1, SqdInvRoot=2, Rel=3 };
+enum class Error: unsigned int { Abs=0, Sqd=1, SqdInvRoot=2, Rel=3, ExpRel=4 };
 Error to_error(const GmshFem& gmshFem)
 {
   std::string error_buffer;
@@ -24,6 +24,7 @@ Error to_error(const GmshFem& gmshFem)
     else if(error_buffer=="root_mean_squared"){return Error::Sqd;}
     else if(error_buffer=="root_mean_squared_inverse_root"){return Error::SqdInvRoot;}
     else if(error_buffer=="relative_mean"){return Error::Rel;}
+    else if(error_buffer=="exp_relative_mean"){return Error::ExpRel;}
   }
   return Error::Sqd;
 }
@@ -41,9 +42,6 @@ int error(const GmshFem& gmshFem)
 
     ConfigurationInterface* configuration = newConfiguration(name, parametrization, gmshFem);
 
-    double area = integrate(1._d_sf, configuration->model_unknown(Support::BLK),"Gauss0");
-    msg::print << "Model unknown area = " << area << msg::endl;
-
     unsigned int write_error_field = 0;
     gmshFem.userDefinedParameter(write_error_field, "write_error_field");
     std::string inverted_model_path = "";
@@ -62,6 +60,54 @@ int error(const GmshFem& gmshFem)
     {
         throw Exception("An error parameter could not be found");
     }
+
+    std::string mask_type = "none";
+    if(!gmshFem.userDefinedParameter(mask_type, "error_mask_type"))
+    {
+        msg::warning << "Error mask type oculd not be found. No mask is used (default)." << msg::endl;
+    }
+    ScalarFunction<std::complex<double>> mask;
+    if(mask_type == "none")
+    {
+        mask = 1.;
+    }
+    else if(mask_type == "circle")
+    {
+        double r0 = 3.;
+
+        mask = heaviside< std::complex<double> >
+        (
+          r0*r0  - pow(x< std::complex<double> >(),2) - pow(y< std::complex<double> >(),2)
+        );
+    }
+    else if(mask_type == "horizontal")
+    {
+        double y0 = -2.;
+
+        mask = heaviside< std::complex<double> >
+        (
+          y< std::complex<double> >()-y0
+        );
+    }
+    else if(mask_type == "vertical-slice")
+    {
+        double x0 = 0.008;
+        double x1 = 0.017;
+        mask = heaviside< std::complex<double> >
+        (
+          (x< std::complex<double> >()-x0)*(x1-x< std::complex<double> >())
+        );
+    }
+    else
+    {
+        throw Exception("Mask type "+mask_type+" is unknown.");
+    }
+
+    double area = std::real(integrate(mask, configuration->model_unknown(Support::BLK),"Gauss0"));
+    msg::print << "Model mask unknown area = " << area << msg::endl;
+
+
+
     unsigned int list_buf = 0;
     if(!gmshFem.userDefinedParameter(list_buf, "inverted_model_list"))
     {
@@ -69,6 +115,12 @@ int error(const GmshFem& gmshFem)
     }
     bool list = ((bool)list_buf);
 
+    unsigned int skip = 1;
+    if(!gmshFem.userDefinedParameter(skip, "inverted_model_skip"))
+    {
+        msg::warning << "Skip number could not be found. Defautl value (1) is used." << msg::endl;
+    }
+
     std::vector<std::string> filename(n_it+1);
     if(list)
     {
@@ -96,7 +148,7 @@ int error(const GmshFem& gmshFem)
         {
             for (unsigned int n = 0; n <= n_it; n++)
             {
-                filename[n] = inverted_model_path+inverted_model_name+inverted_model_prefix+std::to_string(n)+inverted_model_suffix;
+                filename[n] = inverted_model_path+inverted_model_name+inverted_model_prefix+std::to_string(n*skip)+inverted_model_suffix;
             }
         }
         else
@@ -104,7 +156,7 @@ int error(const GmshFem& gmshFem)
             filename[0] = inverted_model_path+inverted_model_name+"_m_initial";
             for (unsigned int n = 1; n <= n_it; n++)
             {
-                filename[n] = inverted_model_path+inverted_model_name+inverted_model_prefix+"m"+std::to_string(n);
+                filename[n] = inverted_model_path+inverted_model_name+inverted_model_prefix+"m"+std::to_string(n*skip);
             }
         }
 
@@ -125,30 +177,41 @@ int error(const GmshFem& gmshFem)
             gmsh::model::setCurrent(configuration->model_gmodel());
 
             ScalarFunction<std::complex<double>> error_density;
-
-            std::complex<double> error;
             switch (error_type)
             {
                 case Error::Abs:
                     error_density = abs(configuration->m0()[c]-m);
-                    error = integrate(error_density, configuration->model_unknown(Support::BLK),integration_type+std::to_string(integration_degree)) / area;
                     break;
                 case Error::Rel:
                     error_density = abs( (configuration->m0()[c]-m) / configuration->m0()[c] );
-                    error = integrate(error_density, configuration->model_unknown(Support::BLK),integration_type+std::to_string(integration_degree)) / area;
+                    break;
+                case Error::ExpRel:
+                    error_density = abs( (exp(configuration->m0()[c])-exp(m)) / exp(configuration->m0()[c]) );
                     break;
                 case Error::Sqd:
                     error_density = pow(abs(configuration->m0()[c]-m),2);
-                    error = std::sqrt(integrate(error_density, configuration->model_unknown(Support::BLK),integration_type+std::to_string(integration_degree)) / area);
                     break;
                 case Error::SqdInvRoot:
                     error_density = pow(abs(sqrt(1./configuration->m0()[c])-sqrt(1./m)),2);
-                    error = std::sqrt(integrate(error_density, configuration->model_unknown(Support::BLK),integration_type+std::to_string(integration_degree)) / area);
+                    break;
+            }
+
+            std::complex<double> error;
+            switch (error_type)
+            {
+                case Error::Abs:
+                case Error::Rel:
+                case Error::ExpRel:
+                    error = integrate(mask * error_density, configuration->model_unknown(Support::BLK),integration_type+std::to_string(integration_degree)) / area;
+                    break;
+                case Error::Sqd:
+                case Error::SqdInvRoot:
+                    error = std::sqrt(integrate(mask * error_density, configuration->model_unknown(Support::BLK),integration_type+std::to_string(integration_degree)) / area);
                     break;
             }
 
             msg::print << "Inverted model #" << n << " error = " << error << msg::endl;
-            output << error << csv::endl;
+            output << error;
 
             if(write_error_field)
             {
@@ -157,6 +220,7 @@ int error(const GmshFem& gmshFem)
             }
 
         }
+        output << csv::endl;
     }
 
     return 0;

input/paper2/concrete3/error/ssm_lbfgs_ls_error.txt

+name = ssm_lbfgs_ls_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Paper2/extended/concrete3/SSM_LBFGS_LS/
+inverted_model_name = paper2_concrete3_ssm_lbfgs_ls
+inverted_model_prefix = g0s0
+inverted_model_number = 88
+#Output
+write_error_field = 0

input/paper2/concrete3/error/sst_lbfgs_ls_tik1_error.txt

+name = sst_lbfgs_ls_tik1_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Paper2/extended/concrete3/SST_LBFGS_LS_TIK/
+inverted_model_name = paper2_concrete3_sst_lbfgs_ls_tik1
+inverted_model_prefix = g0s0
+inverted_model_number = 200
+#Output
+write_error_field = 0

input/paper3/marmousi/error/df5_lbfgs_ls_snr3_error.txt → input/paper2/concrete3/error/us_lbfgs_ls_tik1_error.txt

-name = df5_lbfgs_ls_snr3_error
+name = us_lbfgs_ls_tik1_error
 #Configuration
-unknown = subsurface
+unknown = all
 m0_typec0 = file
-m0_pathc0 = ../input/paper2/marmousi_slowness2
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
 #Error
-error_type = root_mean_squared
+error_type = exp_relative_mean
 integration_degree = 8
 #Input
-inverted_model_path = ../../Manuscript/Part3/TC1/SNR3/DF5_LBFGS_LS_SNR3/
-inverted_model_name = paper2_extended_df5_lbfgs_ls_snr3
+inverted_model_path = ../../Paper2/extended/concrete3/US_LBFGS_LS_TIK/
+inverted_model_name = paper2_concrete3_us_lbfgs_ls_tik1
 inverted_model_prefix = g0s0
 inverted_model_number = 200
 #Output

input/paper3/concrete3/error/df/df0_lbfgs_ls_snrm3_error.txt

+name = df0_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/DF_LBFGS_LS/
+inverted_model_name = paper2_concrete3_df0_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 200
+#Output
+write_error_field = 0

input/paper3/concrete3/error/df/df1_lbfgs_ls_snrm3_error.txt

+name = df1_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/DF_LBFGS_LS/
+inverted_model_name = paper2_concrete3_df1_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 193
+#Output
+write_error_field = 0

input/paper3/concrete3/error/df/df2_lbfgs_ls_snrm3_error.txt

+name = df2_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/DF_LBFGS_LS/
+inverted_model_name = paper2_concrete3_df2_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 250
+#Output
+write_error_field = 0

input/paper3/concrete3/error/df/df3_lbfgs_ls_snrm3_error.txt

+name = df3_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/DF_LBFGS_LS/
+inverted_model_name = paper2_concrete3_df3_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 250
+#Output
+write_error_field = 0

input/paper3/concrete3/error/df/df4_lbfgs_ls_snrm3_error.txt

+name = df4_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/DF_LBFGS_LS/
+inverted_model_name = paper2_concrete3_df4_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 250
+#Output
+write_error_field = 0

input/paper3/concrete3/error/df/df5_lbfgs_ls_snrm3_error.txt

+name = df5_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/DF_LBFGS_LS/
+inverted_model_name = paper2_concrete3_df5_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 250
+#Output
+write_error_field = 0

input/paper3/concrete3/error/df_tv/df_lbfgs_ls_snrm3_tv_error.txt

+name = df_lbfgs_ls_snrm3_tv_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/DF_LBFGS_LS_TV/
+inverted_model_number = 10
+inverted_model_list = 1
+inverted_model_name0 = paper2_concrete3_df_lbfgs_ls_snrm3_tv0_mg0
+inverted_model_name1 = paper2_concrete3_df_lbfgs_ls_snrm3_tv1_mg0
+inverted_model_name2 = paper2_concrete3_df_lbfgs_ls_snrm3_tv2_mg0
+inverted_model_name3 = paper2_concrete3_df_lbfgs_ls_snrm3_tv3_mg0
+inverted_model_name4 = paper2_concrete3_df_lbfgs_ls_snrm3_tv4_mg0
+inverted_model_name5 = paper2_concrete3_df_lbfgs_ls_snrm3_tv5_mg0
+inverted_model_name6 = paper2_concrete3_df_lbfgs_ls_snrm3_tv6_mg0
+inverted_model_name7 = paper2_concrete3_df_lbfgs_ls_snrm3_tv7_mg0
+inverted_model_name8 = paper2_concrete3_df_lbfgs_ls_snrm3_tv8_mg0
+inverted_model_name9 = paper2_concrete3_df_lbfgs_ls_snrm3_tv9_mg0
+inverted_model_name10 = paper2_concrete3_df_lbfgs_ls_snrm3_tv10_mg0
+#Output
+write_error_field = 0

input/paper3/concrete3/error/ssm/ssm0_lbfgs_ls_snrm3_error.txt

+name = ssm0_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/SSM_LBFGS_LS/
+inverted_model_name = paper2_concrete3_ssm0_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 151
+#Output
+write_error_field = 0

input/paper3/concrete3/error/ssm/ssm1_lbfgs_ls_snrm3_error.txt

+name = ssm1_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/SSM_LBFGS_LS/
+inverted_model_name = paper2_concrete3_ssm1_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 250
+#Output
+write_error_field = 0

input/paper3/concrete3/error/ssm/ssm2_lbfgs_ls_snrm3_error.txt

+name = ssm2_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/SSM_LBFGS_LS/
+inverted_model_name = paper2_concrete3_ssm2_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 250
+#Output
+write_error_field = 0

input/paper3/concrete3/error/ssm/ssm3_lbfgs_ls_snrm3_error.txt

+name = ssm3_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/SSM_LBFGS_LS/
+inverted_model_name = paper2_concrete3_ssm3_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 250
+#Output
+write_error_field = 0

input/paper3/concrete3/error/ssm/ssm4_lbfgs_ls_snrm3_error.txt

+name = ssm4_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/SSM_LBFGS_LS/
+inverted_model_name = paper2_concrete3_ssm4_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 250
+#Output
+write_error_field = 0

input/paper3/concrete3/error/ssm/ssm5_lbfgs_ls_snrm3_error.txt

+name = ssm5_lbfgs_ls_snrm3_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/SSM_LBFGS_LS/
+inverted_model_name = paper2_concrete3_ssm5_lbfgs_ls_snrm3
+inverted_model_prefix = g0s0
+inverted_model_number = 250
+#Output
+write_error_field = 0

input/paper3/concrete3/error/ssm_tik/ssm_lbfgs_ls_snrm3_tik_error.txt

+name = ssm_lbfgs_ls_snrm3_tik_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/SSM_LBFGS_LS_TIK/
+inverted_model_number = 10
+inverted_model_list = 1
+inverted_model_name0 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik0_mg0
+inverted_model_name1 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik1_mg0
+inverted_model_name2 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik2_mg0
+inverted_model_name3 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik3_mg0
+inverted_model_name4 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik4_mg0
+inverted_model_name5 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik5_mg0
+inverted_model_name6 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik6_mg0
+inverted_model_name7 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik7_mg0
+inverted_model_name8 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik8_mg0
+inverted_model_name9 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik9_mg0
+inverted_model_name10 = paper2_concrete3_ssm_lbfgs_ls_snrm3_tik10_mg0
+#Output
+write_error_field = 0

input/paper3/concrete3/error/tik/sst_lbfgs_ls_snrm3_tik_error.txt

+name = sst_lbfgs_ls_snrm3_tik_error
+#Configuration
+unknown = all
+m0_typec0 = file
+m0_pathc0 = ../input/paper2/concrete3/large_concrete_lns2
+#Error
+error_type = exp_relative_mean
+integration_degree = 8
+#Input
+inverted_model_path = ../../Manuscript/Part3bis/TC2/SNRM3/SST_LBFGS_LS_TIK/
+inverted_model_number = 7
+inverted_model_list = 1
+inverted_model_name0 = paper2_concrete3_sst_lbfgs_ls_snrm3_tik0_mg0
+inverted_model_name1 = paper2_concrete3_sst_lbfgs_ls_snrm3_tik1_mg0
+inverted_model_name2 = paper2_concrete3_sst_lbfgs_ls_snrm3_tik2_mg0
+inverted_model_name3 = paper2_concrete3_sst_lbfgs_ls_snrm3_tik3_mg0
+inverted_model_name4 = paper2_concrete3_sst_lbfgs_ls_snrm3_tik4_mg0
+inverted_model_name5 = paper2_concrete3_sst_lbfgs_ls_snrm3_tik5_mg0
+inverted_model_name6 = paper2_concrete3_sst_lbfgs_ls_snrm3_tik6_mg0
+inverted_model_name7 = paper2_concrete3_sst_lbfgs_ls_snrm3_tik7_mg0
+#Output
+write_error_field = 0