From 15ef29162c9d274e82b1db58027a51d25bef8884 Mon Sep 17 00:00:00 2001
From: XavierAdriaens <xavier.adriaens@doct.uliege.be>
Date: Tue, 28 Sep 2021 13:59:49 +0200
Subject: [PATCH] Temporary files
---
input/paper3/ob_comparison_cls.txt | 28 ++
input/paper3/robustness/inversion.txt | 56 ++++
input/paper3/robustness/inversion_nscgB.txt | 20 ++
multiscale_data_plot.py | 80 ++++++
multiscale_plot2D.py | 135 ++++++++++
paper2_reference_highly_data0.csv | 11 +
paper3.py | 224 ++++++++++++++++
paper31.py | 163 ++++++++++++
paper32.py | 132 ++++++++++
paper33.py | 180 +++++++++++++
paper3Chap2V0.py | 275 ++++++++++++++++++++
paper3ChapOtherExampleV0.py | 263 +++++++++++++++++++
specific/data/objective/phasederivative.cpp | 157 +++++++++++
specific/data/objective/phasederivative.h | 74 ++++++
specific/data/objective/sobolev.cpp | 44 +++-
15 files changed, 1831 insertions(+), 11 deletions(-)
create mode 100644 input/paper3/ob_comparison_cls.txt
create mode 100644 input/paper3/robustness/inversion.txt
create mode 100644 input/paper3/robustness/inversion_nscgB.txt
create mode 100644 multiscale_data_plot.py
create mode 100644 multiscale_plot2D.py
create mode 100644 paper2_reference_highly_data0.csv
create mode 100644 paper3.py
create mode 100644 paper31.py
create mode 100644 paper32.py
create mode 100644 paper33.py
create mode 100644 paper3Chap2V0.py
create mode 100644 paper3ChapOtherExampleV0.py
create mode 100644 specific/data/objective/phasederivative.cpp
create mode 100644 specific/data/objective/phasederivative.h
diff --git a/input/paper3/ob_comparison_cls.txt b/input/paper3/ob_comparison_cls.txt
new file mode 100644
index 0000000..8813004
--- /dev/null
+++ b/input/paper3/ob_comparison_cls.txt
@@ -0,0 +1,28 @@
+name = ob_comparison_cls
+#Configuration
+unknown = soil
+m0_type = velocity_file
+filter_scale = 2.0
+path = ../input/paper3/marmousi
+#
+#Equation
+gaussNewton=0
+#Objective
+objective = continuousl2distance
+objective_weight0 = ../input/paper3/paper3_reference_data0
+#
+interval = linear
+scale0 = 0.4
+scaleN = 4.
+N=9
+#
+#Data
+data0 = paper3_synthetics_data0
+#
+#Discretization
+#Wave
+wave_FunctionSpaceDegree0 = 2
+#Model
+model_IntegrationType = Gauss
+model_FunctionSpaceType = HierarchicalH1
+model_FunctionSpaceDegree = 1
diff --git a/input/paper3/robustness/inversion.txt b/input/paper3/robustness/inversion.txt
new file mode 100644
index 0000000..408483a
--- /dev/null
+++ b/input/paper3/robustness/inversion.txt
@@ -0,0 +1,56 @@
+#Configuration
+unknown = soil
+m0_type = velocity_depth
+Re(v_0) = 1.5
+Im(v_0) = 0.
+Re(v_H) = 4.5
+Im(v_H) = 0.
+#
+#Equation
+gaussNewton=0
+#
+#Global minimum search
+globalminimum_n_scale=1
+#
+scale00 = 0.8
+#
+scale10 = 0.5
+#
+scale20 = 0.4
+#
+n_group=3
+#Data
+data0 = paper3_synthetics_data0
+data1 = paper3_synthetics_data1
+data2 = paper3_synthetics_data2
+#
+#Discretization
+#Wave
+wave_FunctionSpaceDegree0 = 2
+wave_FunctionSpaceDegree1 = 3
+wave_FunctionSpaceDegree2 = 4
+#Model
+model_IntegrationType = Gauss
+model_FunctionSpaceType = HierarchicalH1
+model_FunctionSpaceDegree0 = 1
+model_FunctionSpaceDegree1 = 1
+model_FunctionSpaceDegree2 = 1
+#
+#Global minimum search
+globalminimumsearch = sequentialscale
+globalminimum_relDecreaseThreshold = 0.0
+#Local minimum search
+#localminimum_absDecreaseThreshold = 1.5e-4
+localminimum_absDecreaseThreshold = 3e-4
+#
+objective = continuousl2distance
+objective_weight0 = ../input/paper3/paper3_reference_data0
+objective_weight1 = ../input/paper3/paper3_reference_data1
+objective_weight2 = ../input/paper3/paper3_reference_data2
+#
+#Inner product
+complex = 0
+boundaryWeight = 0.
+innerproduct = sobolev
+diag_preconditioner = 1
+ref_preconditioner = 250
diff --git a/input/paper3/robustness/inversion_nscgB.txt b/input/paper3/robustness/inversion_nscgB.txt
new file mode 100644
index 0000000..6293d9d
--- /dev/null
+++ b/input/paper3/robustness/inversion_nscgB.txt
@@ -0,0 +1,20 @@
+name = inversion_nscgB
+#
+#Local minimum search
+localminimumsearch = yuanfan
+localminimum_initMu = 1
+#localminimum_c0 = 0.0001
+localminimum_c0 = 0.5
+localminimum_c1 = 0.75
+localminimum_c2 = 0.75
+localminimum_c3 = 0.25
+localminimum_c4 = 2.
+#localminimum_writeInterval = 5
+#localminimum_maxIteration = 50
+localminimum_writeInterval = 1
+localminimum_maxIteration = 20
+#
+#Descent search
+descentsearch = newton_steihaug_conjugategradient
+descent_maxIteration = 50
+descent_eta0 = 0.4
diff --git a/multiscale_data_plot.py b/multiscale_data_plot.py
new file mode 100644
index 0000000..e80d4c0
--- /dev/null
+++ b/multiscale_data_plot.py
@@ -0,0 +1,80 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Dec 7 13:32:55 2020
+
+@author: xavier
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':6})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 1.5;
+W = 4.5;
+
+freq = np.array([4.,6.,9.])
+L = 4.356
+n_freq= 1#np.size(freq)
+
+d0 = (np.loadtxt('MP2008LinearVelocity2D_synthetics_data0.csv',delimiter=';',max_rows=1,skiprows=0)).view(np.complex128);
+Nr = np.size(d0)
+Δxr = L / (Nr-1)
+
+tmp = np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',max_rows=1,dtype=np.integer);
+Nv = tmp[0];
+Na = tmp[1];
+N = Na;
+
+d = np.zeros([Nr,Na])*1j
+r = np.zeros([Nr,Na])*1j
+sr = np.zeros([Nr,Na])*1j
+dsr = np.zeros([Nr,Na])*1j
+
+def moving_average(d,n):
+ mask=np.ones(n)/n
+ mean = np.convolve(d,mask,'same')
+ border = 1./np.convolve(np.ones(Nr),mask,'same')
+ return border*mean
+
+def grad(d):
+ dd = 1j*np.zeros(Nr)
+ dd[0] = (d[1]-d[0])/Δxr
+ for i in range(1,Nr-1):
+ dd[i] = (d[i+1]-d[i-1])/2./Δxr
+ dd[Nr-1]=(d[Nr-1]-d[Nr-2])/Δxr
+ return dd
+def distL2(d,d0,n):
+ r = 1-d/d0
+ r = moving_average(r,n)
+ return 1./2.*np.linalg.norm(r)**2
+def pdistH1(d,d0,n):
+ r=1-d/d0
+ #dr = grad(r)
+ dr = moving_average(grad(moving_average(r,n)),n)
+ #dr = grad(moving_average(r,n))
+ return 1./2.*np.linalg.norm(dr)**2
+
+f0 = plt.figure(figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+for i in range(0,N):
+ d[:,i] = (np.loadtxt('MP2008LinearVelocity2D_multiscale_df0v0a'+str(i)+'.csv',delimiter=';')).view(np.complex128);
+ r[:,i] = d[:,i]-d0;
+ sr[:,i] = r[:,i]/d0;
+ dsr[:,i] = grad(sr[:,i])
+plt.plot(dsr)
+
+plt.plot(np.abs(dsr)**2)
+plt.plot(np.abs(r)**2)
+
+n = 1
+JL2 = np.zeros(Na)
+JpH1 = np.zeros(Na)
+for i in range(0,Na):
+ JL2[i] = distL2(d[:,i],d0,n)
+ JpH1[i] = pdistH1(d[:,i],d0,n)
+plt.plot(JL2/np.max(JL2))
+plt.plot(JpH1/np.max(JpH1))
\ No newline at end of file
diff --git a/multiscale_plot2D.py b/multiscale_plot2D.py
new file mode 100644
index 0000000..2f4d647
--- /dev/null
+++ b/multiscale_plot2D.py
@@ -0,0 +1,135 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jan 8 11:46:10 2021
+
+@author: xavier
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import cm
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':6})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 1.5;
+W = 4.5;
+
+"""
+Example III: Velocity slope and mean perturbation
+"""
+L = 8.736
+λ = 2/6
+a0 = 0.5
+v0 = 2.0
+d0 = (np.loadtxt('MP2008LinearVelocity2D_synthetics_data0.csv',delimiter=';',max_rows=1,skiprows=0)).view(np.complex128);
+N = np.size(d0)
+xr = np.linspace(-L/2,L/2,N)
+Δxr = L/(N-1)
+ξr = np.linspace(-1/2./Δxr,1/2./Δxr,N)
+
+tmp = np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',max_rows=1,dtype=np.integer);
+Nv = tmp[0]+1;
+Na = tmp[1]+1;
+V = (np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',skiprows=1)[:,0]).reshape(Nv,Na);
+v = V[:,0]
+A = (np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',skiprows=1)[:,1]).reshape(Nv,Na);
+a = A[0,:]
+
+d_vsm = np.zeros([N,Nv,Na])*1j
+for i in range(0,Na):
+ for j in range(0,Nv):
+ tmp = (np.loadtxt('MP2008LinearVelocity2D_multiscale_df0v'+str(i)+'a'+str(j)+'.csv',delimiter=';')).view(np.complex128);
+ d_vsm[:,i,j] = tmp[0:N]
+
+def d(xr,vϵ,aϵ):
+ i = np.searchsorted(v,vϵ)
+ j = np.searchsorted(a,aϵ)
+ if i > Nv-1:
+ i = Nv-1
+ if j > Na-1:
+ j = Na-1
+ return d_vsm[:,i,j]
+
+"""
+SPATIAL DOMAIN ANALYSIS
+"""
+aϵ = 0.55
+vϵ = 2.1
+
+fig = plt.figure(1)
+ax = plt.axes(projection='3d')
+plt.show()
+ax.plot(np.real(d(xr,v0,a0)),np.imag(d(xr,v0,a0)),xr)
+ax.plot(np.real(d(xr,vϵ,aϵ)),np.imag(d(xr,vϵ,aϵ)),xr)
+
+def r(vϵ,aϵ):
+ return (d(xr,v0,a0)-d(xr,vϵ,aϵ))/d(xr,v0,a0)
+
+fig = plt.figure(2)
+ax = plt.axes(projection='3d')
+ax.set_xlim([-1,1])
+ax.set_ylim([-1,1])
+ax.set_zlim([-L/2,L/2])
+ax.plot(np.real(r(vϵ,aϵ)),np.imag(r(vϵ,aϵ)),xr)
+
+"""
+FREQUENCY DOMAIN ANALYSIS
+"""
+s1=100*L
+#s2=1.25*λ#L/np.sqrt(12)
+s3=6*λ
+s2=np.sqrt(L*s3/2)
+def σ1(ξr):
+ return 1.
+def σ2(ξr):
+ return s2 * 2*np.pi*1j*ξr
+def σ3(ξr):
+ return 1 / (1 + s1*2*np.pi*1j*ξr)
+def σ4(ξr):
+ window = np.ones(np.size(ξr))
+ window[np.greater(ξr,1/s3)] = 0
+ window[np.less(ξr,-1/s3)] = 0
+ #return (2*np.pi*1j*(s2*ξr)) * window
+ return s2 * 2*np.pi*1j*ξr / (1 + 2*np.pi*1j*s3*ξr )
+ #return (ξr*s1) / (1 + (ξr*s2)**2)
+ #return s1*2*np.pi*1j*ξr / (1 + s2**2*np.abs(2*np.pi*1j*ξr)**2)
+def f_r(vϵ,aϵ):
+ return np.fft.fftshift(np.fft.fft(r(vϵ,aϵ)))
+def max_f_r(vϵ,aϵ):
+ return np.max(np.abs(f_r(vϵ,aϵ)))
+
+def distance1(vϵ,aϵ):
+ return 0.5*(np.linalg.norm(σ1(ξr)*f_r(vϵ,aϵ)))**2 / (N)**2
+v_distance1 = np.vectorize(distance1)
+
+def distance2(vϵ,aϵ):
+ return 0.5*(np.linalg.norm(σ2(ξr)*f_r(vϵ,aϵ)))**2 / (N)**2
+v_distance2 = np.vectorize(distance2)
+
+def distance3(vϵ,aϵ):
+ return 0.5*(np.linalg.norm(σ3(ξr)*f_r(vϵ,aϵ)))**2 / (N)**2
+v_distance3 = np.vectorize(distance3)
+
+def distance4(vϵ,aϵ):
+ return 0.5*(np.linalg.norm(σ4(ξr)*f_r(vϵ,aϵ)))**2 / (N)**2
+v_distance4 = np.vectorize(distance4)
+
+"""
+COST FUNCTION
+"""
+fig = plt.figure()
+ax = fig.gca(projection='3d')
+surf = ax.plot_surface(V, A, v_distance4(V,A), cmap=cm.coolwarm,
+ linewidth=0, antialiased=False)
+ax.set_zlim(0, 1.0)
+
+plt.pcolormesh(V,A,v_distance4(V,A), vmin=0, vmax=1.25)
+plt.colorbar()
+
+plt.plot(a,v_distance1(v0,a))
+plt.plot(a,v_distance4(v0,a))
+plt.plot(v,v_distance1(v,a0))
+plt.plot(v,v_distance4(v,a0))
\ No newline at end of file
diff --git a/paper2_reference_highly_data0.csv b/paper2_reference_highly_data0.csv
new file mode 100644
index 0000000..1308c36
--- /dev/null
+++ b/paper2_reference_highly_data0.csv
@@ -0,0 +1,11 @@
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
+1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00;1.000000000000000000e+00;0.000000000000000000e+00
diff --git a/paper3.py b/paper3.py
new file mode 100644
index 0000000..9d233e0
--- /dev/null
+++ b/paper3.py
@@ -0,0 +1,224 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':6})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 1.5;
+W = 4.5;
+
+"""
+Example I: Wave number perturbation
+
+L = 10.
+Lpad = 10*L
+λ = 1.
+k = 2*np.pi / λ
+N = 2**8*2**3
+xr = np.linspace(-Lpad/2,Lpad/2,N)
+Δxr = Lpad/(N-1)
+ξr = np.linspace(-1/2./Δxr,1/2./Δxr,N)
+
+def rect(xr):
+ if (xr > -L/2) and (xr < L/2):
+ return 1+0j
+ elif (xr == -L/2) or (xr == L/2):
+ return 1.0+0j
+ else:
+ return 0j
+v_rect = np.vectorize(rect)
+def sign(xr):
+ if (xr > 0):
+ return 1
+ elif (xr < 0):
+ return -1
+ else:
+ return 0
+v_sign = np.vectorize(sign)
+def d(xr,ϵ):
+ return np.exp(1j*k*(1+ϵ)*xr)*v_rect(xr)+1e-12*(1-v_rect(xr))
+ #return np.cos(k*(1+ϵ)*xr)*v_rect(xr)+1e-12*(1-v_rect(xr))
+"""
+
+"""
+Example II: Wave number perturbation
+L = 10.
+N = 257
+xr = np.linspace(-L/2,L/2,N)
+Δxr = L/(N-1)
+ξr = np.linspace(-1/2./Δxr,1/2./Δxr,N)
+k = 2*np.pi
+λ = 1.
+h = 2.
+def d(xr,ϵ):
+ r = np.sqrt(xr**2 + (h*(1+ϵ))**2)
+ return np.exp(1j*k*(r+h*(1+ϵ)))
+"""
+
+"""
+Example III: Velocity slope perturbation
+
+
+L = 8.712-0.012
+λ = 2/6
+a0 = 0.4
+d0 = (np.loadtxt('MP2008LinearVelocity2D_synthetics_data0.csv',delimiter=';',max_rows=1,skiprows=0)).view(np.complex128);
+N = np.size(d0)-1
+xr = np.linspace(-L/2,L/2,N)
+Δxr = L/(N-1)
+ξr = np.linspace(-1/2./Δxr,1/2./Δxr,N)
+
+tmp = np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',max_rows=1,dtype=np.integer);
+Nv = tmp[0];
+Na = tmp[1];
+a = np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',skiprows=1)[:,1];
+
+d_vs = np.zeros([N,Na])*1j
+for i in range(0,Na):
+ tmp = (np.loadtxt('MP2008LinearVelocity2D_multiscale_df0v0a'+str(i)+'.csv',delimiter=';')).view(np.complex128);
+ d_vs[:,i] = tmp[0:N]
+
+def d(xr,ϵ):
+ aϵ = a0*(1+ϵ)
+ n = np.searchsorted(a,aϵ)
+ if n > Na-1:
+ n = Na-1
+ return d_vs[:,n]
+"""
+
+
+"""
+SPATIAL DOMAIN ANALYSIS
+"""
+ϵ0 = 0.25
+fig1, (ax1re, ax1im) = plt.subplots(2, 1)
+ax1re.plot(xr,np.real(d(xr,0)))
+ax1re.plot(xr,np.real(d(xr,ϵ0)))
+ax1im.plot(xr,np.imag(d(xr,0)))
+ax1im.plot(xr,np.imag(d(xr,ϵ0)))
+
+def w(ϵ):
+ return d(xr,ϵ)/d(xr,0)*v_rect(xr)
+def r(ϵ):
+ return (d(xr,0)-d(xr,ϵ))/d(xr,0)
+
+fig2, (ax2re, ax2im, ax2abs) = plt.subplots(3, 1)
+ax2re.plot(xr,np.real(w(ϵ0)))
+ax2im.plot(xr,np.imag(w(ϵ0)))
+ax2abs.plot(xr,np.abs(w(ϵ0)))
+ax2re.plot(xr,np.real(r(ϵ0)))
+ax2im.plot(xr,np.imag(r(ϵ0)))
+ax2abs.plot(xr,np.abs(r(ϵ0)))
+
+
+"""
+fig = plt.figure(2)
+ax = plt.axes(projection='3d')
+ax.set_xlim([-1,1])
+ax.set_ylim([-1,1])
+ax.set_zlim([-L/2,L/2])
+ax.plot(np.real(r(ϵ0)),np.imag(r(ϵ0)),xr)
+"""
+
+"""
+FREQUENCY DOMAIN ANALYSIS
+"""
+def f_d(ϵ):
+ #return np.fft.fftshift(np.fft.fft(np.fft.fftshift(d(xr,ϵ))))
+ #return np.fft.fftshift(np.fft.fft(d(xr,ϵ))) / N * np.exp(1j*(-np.pi*ξr*(Δxr)))#* np.exp(1j*(+np.pi*ξr*(Lpad-Δxr)))
+ return np.fft.fftshift(np.fft.fft(d(xr,ϵ)))* -1j * np.exp(1j*np.pi*ξr*(Lpad-Δxr)) #/ N * np.exp(1j*(-np.pi*ξr*(Δxr)))#
+def f_r(ϵ):
+ return np.fft.fftshift(np.fft.fft(r(ϵ))) / N
+def f_w(ϵ):
+ return np.fft.fftshift(np.fft.fft(w(ϵ))) / N
+
+fig3, (ax3re, ax3im, ax3abs) = plt.subplots(3, 1)
+ax3abs.plot(ξr,np.abs(f_d(0)))
+ax3abs.plot(ξr,np.abs(f_d(ϵ0)))
+ax3re.plot(ξr,np.real(f_d(0)))
+ax3re.plot(ξr,np.real(f_d(ϵ0)))
+ax3im.plot(ξr,np.imag(f_d(0)))
+ax3im.plot(ξr,np.imag(f_d(ϵ0)))
+ax3abs.plot(ξr,np.abs(f_w(ϵ0)))
+ax3abs.plot(ξr,np.abs(f_r(ϵ0)))
+
+fig3, (ax3re, ax3im, ax3abs) = plt.subplots(3, 1)
+ax3abs.plot(ξr,np.abs(f_d(0)))
+ax3abs.plot(ξr,np.abs(f_d(ϵ0)))
+ax3re.plot(ξr,np.angle(f_d(0)))
+ax3re.plot(ξr,np.angle(f_d(ϵ0)))
+
+fig4 = plt.figure()
+ax = plt.axes(projection='3d')
+ax.set_zlim([0,2/λ])
+ax.plot(np.real(f_d(0)),np.imag(f_d(0)),ξr)
+ax.plot(np.real(f_d(ϵ0)),np.imag(f_d(ϵ0)),ξr)
+ax.plot(np.real(f_w(ϵ0)),np.imag(f_w(ϵ0)),ξr)
+
+s1=100*L
+s2=1.25*λ#L/np.sqrt(12)
+s3=λ
+def σ1(ξr):
+ return 1.
+def σ2(ξr):
+ return s2 * 2*np.pi*1j*ξr
+def σ3(ξr):
+ return 1 / (1 + s1*2*np.pi*1j*ξr)
+def σ4(ξr):
+ window = np.ones(np.size(ξr))
+ window[np.greater(ξr,1/s3)] = 0
+ window[np.less(ξr,-1/s3)] = 0
+ return s2 * 2*np.pi*1j*ξr*window
+ #return (ξr*s1) / (1 + (ξr*s2)**2)
+ #return s1*2*np.pi*1j*ξr / (1 + s2**2*np.abs(2*np.pi*1j*ξr)**2)
+def max_f_r(ϵ):
+ return np.max(np.abs(f_r(ϵ)))
+
+v_max_f_r = np.vectorize(max_f_r)
+fig = plt.figure(3)
+ax = plt.axes(projection='3d')
+"""
+ax.set_xlim([-N,N])
+ax.set_ylim([-N,N])
+ax.set_zlim([-1/Δxr,1/Δxr])
+"""
+ax.plot(np.real(f_r(ϵ0)),np.imag(f_r(ϵ0)),ξr)
+
+fig = plt.figure(4)
+plt.plot(ξr,np.abs(f_corr))
+plt.plot(np.abs(f_corr),ξr)
+
+def distance1(ϵ):
+ return 0.5*(np.linalg.norm(σ1(ξr)*f_r(ϵ)))**2 / (N)**2
+v_distance1 = np.vectorize(distance1)
+
+def distance2(ϵ):
+ return 0.5*(np.linalg.norm(σ2(ξr)*f_r(ϵ)))**2 / (N)**2
+v_distance2 = np.vectorize(distance2)
+
+def distance3(ϵ):
+ return 0.5*(np.linalg.norm(σ3(ξr)*f_r(ϵ)))**2 / (N)**2
+v_distance3 = np.vectorize(distance3)
+
+def distance4(ϵ):
+ return 0.5*(np.linalg.norm(σ4(ξr)*f_r(ϵ)))**2 / (N)**2
+v_distance4 = np.vectorize(distance4)
+
+"""
+COST FUNCTION
+"""
+ϵ = np.linspace(-0.15,0.15,400)
+#Least squares
+plt.plot(ϵ,5*v_distance1(ϵ))
+plt.plot(ϵ,v_distance2(ϵ))
+plt.plot(ϵ,10*v_distance3(ϵ))
+plt.plot(ϵ,v_distance4(ϵ))
+plt.plot(ϵ,5*v_distance3(ϵ)+0.5*v_distance4(ϵ))
+
+plt.plot(ϵ,0.5*k**2*ϵ**2)
+plt.plot(ϵ,0.5*k**2*ϵ**2)
+plt.plot(ϵ,2*np.sin(ϵ*k*h)**2)
+plt.plot(ϵ,2*np.sin(ϵ*k*h/2)**2)
+plt.plot(ϵ,1-np.sin(ϵ*k*L/2)/(ϵ*k*L/2))
diff --git a/paper31.py b/paper31.py
new file mode 100644
index 0000000..ac3e9ba
--- /dev/null
+++ b/paper31.py
@@ -0,0 +1,163 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Jan 19 11:05:45 2021
+
+@author: xavier
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 1.5;
+W = 6.5;
+
+"""
+Example I: Wave number perturbation
+"""
+L = 10.
+Lint = 1*L
+λ = 1.
+k = 2*np.pi / λ
+N = 2**14+1
+Δxr = Lint/2/N
+xr = np.linspace( Δxr/2, (Lint-Δxr)/2,N)
+ξr = (np.arange(0,N)-np.floor(N/2)) /N/Δxr
+
+def J_l2(ϵ):
+ return 1-np.sinc(ϵ*k*L/2/np.pi)
+s=λ
+def J_sob(ϵ):
+ return 0.5*s**2*k**2*ϵ**2
+
+def rect(xr):
+ if (xr > -L/4) and (xr < L/4):
+ return 1+0j
+ elif (xr == -L/4) or (xr == L/4):
+ return 0.5+0j
+ else:
+ return 0j
+v_rect = np.vectorize(rect)
+def d(xr,ϵ):
+ return np.exp(1j*k*(1+ϵ)*xr)*v_rect(xr-L/4)
+
+
+"""
+Example II: Wave number perturbation
+
+L = 10.
+Lint = 1*L
+λ = 1.
+k = 2*np.pi / λ
+h=5.
+N = 2**14+1
+Δxr = Lint/N
+xr = np.linspace( (-Lint+Δxr)/2, (Lint-Δxr)/2,N)
+ξr = (np.arange(0,N)-np.floor(N/2)) /N/Δxr
+
+def J_l2(ϵ):
+ return 1-np.sinc(ϵ*k*L/2/np.pi)
+s=λ
+def J_sob(ϵ):
+ return 0.5*s**2*k**2*ϵ**2
+
+def rect(xr):
+ if (xr > -L/4) and (xr < L/4):
+ return 1+0j
+ elif (xr == -L/4) or (xr == L/4):
+ return 0.5+0j
+ else:
+ return 0j
+v_rect = np.vectorize(rect)
+def d(xr,ϵ):
+ return np.exp(1j*k*h*(1+ϵ))*v_rect(xr-L/4)
+"""
+
+
+def cfft(q):
+ tmp = np.fft.fft(q) * Δxr * np.exp(1j*2*np.pi* ((N-1)/2) * np.arange(0,N) / N) / L * 2
+ cut = np.int(np.floor(N/2))
+ return np.concatenate([ tmp[cut+1:N],tmp[0:cut+1]])
+
+"""
+SPATIAL DOMAIN ANALYSIS
+"""
+Nϵ=38
+#vϵ = np.linspace(0,0.1,Nϵ)
+vϵ = np.linspace(0,0.2,Nϵ)
+vc = np.linspace(0.05,0.95,Nϵ)
+fig1 = plt.figure(1,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,L/2])
+plt.ylim([-1,1])
+#plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.real(d(xr,vϵ[i])),color=plt.cm.gist_yarg(vc[i]))
+
+μ = 1e-16
+def w(ϵ):
+ return d(xr,ϵ)*np.conj(d(xr,0)) / (np.abs(d(xr,0))**2 + μ)
+def r(ϵ):
+ return (np.abs(d(xr,0))**2 - d(xr,ϵ)*np.conj(d(xr,0))) / (np.abs(d(xr,0))**2 + μ)
+
+fig2 = plt.figure(2,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,L/2])
+plt.ylim([0,2])
+#plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.abs(r(vϵ[i])),color=plt.cm.gist_yarg(vc[i]))
+
+"""
+FREQUENCY DOMAIN ANALYSIS
+"""
+def f_d(ϵ):
+ return cfft(d(xr,ϵ))
+
+fig3 = plt.figure(3,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,2])
+#plt.ylim([0,1])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξr,np.abs(f_d(vϵ[i])),color=str(vc[i]),marker='o',linewidth=0)
+
+def f_r(ϵ):
+ return cfft(r(ϵ))
+
+fig4 = plt.figure(4,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([-0.5,0.5])
+#plt.ylim([0,1])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξr,np.abs(f_r(vϵ[i])),color=str(vc[i]),marker='o',linewidth=0)
+
+"""
+COST FUNCTION
+"""
+ϵ = np.linspace(0,0.3,100)
+
+fig5 = plt.figure(5,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,0.3])
+plt.ylim([0,1.5])
+plt.grid()
+
+plt.plot(ϵ,Jl2(ϵ),'k')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],1-np.sinc(vϵ[i]*k*L/2/np.pi),color=str(vc[i]),marker='o',markersize=3.5)
+
+fig6 = plt.figure(6,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,0.3])
+plt.ylim([0,1.5])
+plt.grid()
+s=λ
+plt.plot(ϵ,0.5*s**2*k**2*ϵ**2,'k')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],0.5*s**2*k**2*vϵ[i]**2,color=str(vc[i]),marker='o',markersize=3.5)
diff --git a/paper32.py b/paper32.py
new file mode 100644
index 0000000..246eb40
--- /dev/null
+++ b/paper32.py
@@ -0,0 +1,132 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jan 20 11:44:08 2021
+
+@author: xavier
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 1.5;
+W = 6.5;
+
+"""
+Example I: Wave number perturbation
+"""
+L = 10.
+λ = 1.
+k = 2*np.pi / λ
+N = 2**7+1
+Δxr = L/2/N
+xr = np.linspace(Δxr/2, (L-Δxr)/2, N)
+ξr = (np.arange(0,N)-np.floor(N/2))/N/Δxr
+
+def J_l2(ϵ):
+ return 1-np.sinc(ϵ*k*L/2/np.pi)
+s=λ
+def J_sob(ϵ):
+ return 0.5*s**2*k**2*ϵ**2
+def d(xr,ϵ):
+ return np.exp(1j*k*(1+ϵ)*xr)
+
+#simplifier et rendre indépendant du setting
+def cfft(q):
+ tmp = np.fft.fft(q) * Δxr * np.exp(1j*2*np.pi* ((N-1)/2) * np.arange(0,N) / N) / L * 2
+ cut = np.int(np.floor(N/2))
+ return np.concatenate([ tmp[cut+1:N],tmp[0:cut+1]])
+
+"""
+SPATIAL DOMAIN ANALYSIS
+"""
+Nϵ=38
+vϵ = np.linspace(0,0.2,Nϵ)
+vc = np.linspace(0.05,0.95,Nϵ)
+fig1 = plt.figure(1,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,L/2])
+plt.ylim([-1,1])
+#plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.real(d(xr,vϵ[i])),color=plt.cm.gist_yarg(vc[i]))
+
+μ = 1e-16
+def w(ϵ):
+ return d(xr,ϵ) / d(xr,0)
+def r(ϵ):
+ return 1-w(ϵ)
+
+fig2 = plt.figure(2,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,L/2])
+plt.ylim([0,2])
+#plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.abs(r(vϵ[i])),color=plt.cm.gist_yarg(vc[i]))
+
+"""
+FREQUENCY DOMAIN ANALYSIS
+"""
+Nint = 15
+ξintr = np.linspace(ξr[0],ξr[N-1],N+(N-1)*Nint)
+def intcfft(q):
+ weights = cfft(q)
+ out = np.zeros(N+(N-1)*Nint)*1j
+ for i in range(0,N):
+ #out += np.abs(weights[i]*np.sinc((ξintr-ξr[i])*L/2)*np.exp(-1j*2*np.pi*L/4*(ξintr-ξr[i])))**2
+ out += weights[i]*np.sinc((ξintr-ξr[i])*L/2)
+ return out
+
+def f_d(ϵ):
+ return intcfft(d(xr,ϵ))
+
+
+
+fig3 = plt.figure(3,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,2])
+plt.ylim([0,1])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξintr,np.abs(f_d(vϵ[i])),color=str(vc[i]))
+
+def f_r(ϵ):
+ return intcfft(r(ϵ))
+
+fig4 = plt.figure(4,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([-0.5,0.5])
+#plt.ylim([0,1])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξintr,np.abs(f_r(vϵ[i])),color=str(vc[i]))
+
+"""
+COST FUNCTION
+"""
+ϵ = np.linspace(0,0.3,100)
+
+fig5 = plt.figure(5,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,0.3])
+plt.ylim([0,1.5])
+plt.grid()
+
+plt.plot(ϵ,Jl2(ϵ),'k')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],1-np.sinc(vϵ[i]*k*L/2/np.pi),color=str(vc[i]),marker='o',markersize=3.5)
+
+fig6 = plt.figure(6,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,0.3])
+plt.ylim([0,1.5])
+plt.grid()
+s=λ
+plt.plot(ϵ,0.5*s**2*k**2*ϵ**2,'k')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],0.5*s**2*k**2*vϵ[i]**2,color=str(vc[i]),marker='o',markersize=3.5)
diff --git a/paper33.py b/paper33.py
new file mode 100644
index 0000000..e571aee
--- /dev/null
+++ b/paper33.py
@@ -0,0 +1,180 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jan 20 14:50:58 2021
+
+@author: xavier
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 1.5;
+W = 6.5;
+
+
+"""
+Example IV (numerical):
+"""
+L = 8.736
+f=9
+v0=2.
+a0 = 0.0
+λ = v0/f
+
+d0 = (np.loadtxt('MP2008LinearVelocity2D_synthetics_data0.csv',delimiter=';',max_rows=1,skiprows=0)).view(np.complex128);
+N = np.size(d0)
+Δxr = L/(N-1)
+xr = np.linspace(-L/2, L/2, N)
+#ξr = (np.arange(0,2*N)-np.floor(N))/(2*N)/Δxr
+
+tmp = np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',max_rows=1,dtype=np.integer);
+Nv = tmp[0]+1;
+Na = tmp[1]+1;
+V = (np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',skiprows=1)[:,0]).reshape(Nv,Na);
+v = V[:,0]
+A = (np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',skiprows=1)[:,1]).reshape(Nv,Na);
+a = A[0,:]
+J = (np.loadtxt('MP2008LinearVelocity2D_multiscale0.csv',delimiter=';',skiprows=1)[:,2]).reshape(Nv,Na);
+
+a = A[0,:]
+d_vs = np.zeros([N,Na])*1j
+for i in range(0,Na):
+ tmp = (np.loadtxt('MP2008LinearVelocity2D_multiscale_df0v10a'+str(i)+'.csv',delimiter=';')).view(np.complex128);
+ d_vs[:,i] = (tmp[N-1:2*N] + np.flip(tmp[0:N])) / 2
+ d_vs[0,i] = d_vs[0,i] / 2.
+
+Nϵ=Na+1
+ϵmin=-0.2
+ϵmax=0.0
+ξmax=5
+def d(xr,ϵ):
+ aϵ = a0*(1+ϵ)
+ n = np.searchsorted(a,aϵ)
+ if n > Na-1:
+ n = Na-1
+ return d_vs[:,n]
+
+"""
+SPATIAL DOMAIN ANALYSIS
+"""
+vϵ = np.linspace(ϵmin,ϵmax,Nϵ)
+vc = np.linspace(0.05,0.95,Nϵ)
+fig1 = plt.figure(1,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,L/2])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.real(d(xr,vϵ[i])),color=plt.cm.Spectral(vc[i]))
+
+def w(ϵ):
+ return d(xr,ϵ) / d(xr,0)
+def r(ϵ):
+ return 1-w(ϵ)
+
+fig2 = plt.figure(2,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([0,L/2])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.abs(r(vϵ[i])),color=plt.cm.Spectral(vc[i]))
+
+"""
+FREQUENCY DOMAIN ANALYSIS
+"""
+Nint = 8
+ξintr = np.linspace(ξr[0],ξr[2*N-1],2*N+(2*N-1)*Nint)
+def cfft(q,ξr):
+ out = np.zeros(np.size(ξr))*1j
+ for i in range(0,N):
+ out += q[i] * np.exp(-1j*2*np.pi*xr[i]*ξr)
+ return out / N
+
+def f_d(ξr,ϵ):
+ return cfft(d(xr,ϵ),ξr)
+
+
+fig3 = plt.figure(3,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+#plt.xlim([-ξmax,ξmax])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξintr,np.abs(f_d(ξintr,vϵ[i])+f_d(-ξintr,vϵ[i])),color=plt.cm.Spectral(vc[i]))
+
+def f_r(ξr,ϵ):
+ return cfft(r(ϵ),ξr)
+
+fig4 = plt.figure(4,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+#plt.xlim([-ξmax,ξmax])
+#plt.ylim([0,1])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξr,np.abs(σ(ξr))*np.abs(f_r(ξr,vϵ[i])+f_r(-ξr,vϵ[i])),color=plt.cm.Spectral(vc[i]))
+
+"""
+COST FUNCTION
+"""
+ϵ = np.linspace(ϵmin,2*ϵmax,100)
+
+#Analytical
+fig5 = plt.figure(5,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+#plt.xlim([0,0.25])
+plt.xlim([0,2*ϵmax])
+plt.grid()
+
+plt.plot(ϵ,J_l2(ϵ),'k')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],J_l2(vϵ[i]),color=str(vc[i]),marker='o',markersize=3.5)
+
+fig6 = plt.figure(6,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+#plt.xlim([0,0.25])
+plt.xlim([0,2*ϵmax])
+plt.grid()
+plt.plot(ϵ,J_sob(ϵ),'k')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],J_sob(vϵ[i]),color=str(vc[i]),marker='o',markersize=3.5)
+
+#Numerical
+s=2*L
+def σ(ξr):
+ return L*2*np.pi*1j*ξr / (1 + s*2*np.pi*1j*ξr)
+
+def distance1(ϵ):
+ return 0.5*(np.linalg.norm(0.5*(f_r(ξr,ϵ)+f_r(-ξr,ϵ))))**2
+v_distance1 = np.vectorize(distance1)
+
+def distance2(ϵ):
+ return 0.5*(np.linalg.norm(σ(ξr)*0.5*(f_r(ξr,ϵ)+f_r(-ξr,ϵ))))**2
+v_distance2 = np.vectorize(distance2)
+
+"""
+COST FUNCTION
+"""
+#Least squares
+fig7 = plt.figure(7,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([ϵmin,2*ϵmax])
+#plt.ylim([0,2])
+plt.grid()
+
+#plt.plot(ϵ,v_distance1(ϵ),'k')
+#plt.plot(ϵ,J_l2(ϵ),'tab:gray')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],v_distance1(vϵ[i]),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+
+fig8 = plt.figure(8,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([ϵmin,2*ϵmax])
+#plt.ylim([0,2])
+plt.grid()
+
+#plt.plot(ϵ,v_distance2(ϵ),'k')
+#plt.plot(ϵ,J_sob(ϵ),'tab:gray')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],v_distance2(vϵ[i]),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
\ No newline at end of file
diff --git a/paper3Chap2V0.py b/paper3Chap2V0.py
new file mode 100644
index 0000000..9e29ad5
--- /dev/null
+++ b/paper3Chap2V0.py
@@ -0,0 +1,275 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Mar 4 11:23:16 2021
+
+@author: xavier
+"""
+import numpy as np
+import scipy as sp
+import scipy.special
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 1.5;
+W = 6.5;
+save = False;
+
+
+
+L = 10.
+λ = 1.
+k = 2*np.pi / λ
+N = 2**8+1
+Δxr = L/N
+xr = np.linspace(-(L-Δxr)/2, (L-Δxr)/2, N)
+ξr = (np.arange(0,N)-np.floor(N/2))/L
+
+Nϵ=38
+ϵmax=0.2
+ϵmin=0.
+Δϵ=(ϵmax-ϵmin)/(Nϵ-1)
+ξmin=0
+ξmax=2/λ
+xmin=0
+xmax=L/2
+
+Jmax=1.5
+def J_l2(ϵ):
+ return 1-np.sinc(ϵ*k*L/2/np.pi)
+s=λ
+def J_sob(ϵ):
+ return 0.5*s**2*k**2*ϵ**2
+
+"""
+Example IV (analitycal): Wave number perturbation (amplitude)
+"""
+name = 'example4_an'
+amplitude=True
+if amplitude:
+ def d(xr,ϵ):
+ return -1j/4*sp.special.hankel1(0,k*(1+ϵ)*np.abs(xr))
+else:
+ def d(xr,ϵ):
+ return np.exp(1j*k*(1+ϵ)*np.abs(xr))
+
+"""
+Example IV (numerical): Wave number perturbation (amplitude)
+
+name = 'example4_nu'
+
+Nd=2*Nϵ-1
+d_num = np.zeros([N,Nd])*1j
+for i in range(0,Nd):
+ tmp = (np.loadtxt('paper3_example4_multiscale_df0n'+str(i)+'0'+'.csv',delimiter=';')).view(np.complex128);
+ d_num[:,i] = tmp
+
+def d(xr,ϵ):
+ i = (np.around(xr / Δxr)+(N-1)/2).astype(int)
+ j = (np.around((ϵ+ϵmin) / Δϵ)).astype(int)
+ return d_num[i,j]
+"""
+"""
+SPATIAL DOMAIN ANALYSIS
+"""
+vϵ = np.linspace(ϵmin,ϵmax,Nϵ)
+vc = np.linspace(0.95,0.05,Nϵ)
+fig1 = plt.figure(1,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlim([xmin,xmax])
+plt.xlabel(r'$x_r$')
+plt.ylabel(r'Re($d_0$)')
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.real(d(xr,vϵ[i])),color=plt.cm.Spectral(vc[i]))
+if save: plt.savefig(name+'_d.eps')
+
+def w(ϵ):
+ return d(xr,ϵ) / d(xr,0)
+def r(ϵ):
+ return 1-w(ϵ)
+
+fig2 = plt.figure(2,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlim([xmin,xmax])
+plt.xlabel(r'$x_r$')
+plt.ylabel(r'$\left\vert r \right\vert$')
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.abs(r(vϵ[i])),color=plt.cm.Spectral(vc[i]))
+if save: plt.savefig(name+'_r.eps')
+
+"""
+FREQUENCY DOMAIN ANALYSIS
+"""
+Nint = 8
+ξintr = np.linspace(ξr[0],ξr[N-1],N+(N-1)*Nint)
+def cfft(q,ξr):
+ out = np.zeros(np.size(ξr))*1j
+ for i in range(0,N):
+ out += np.nan_to_num(q[i]) * np.exp(-1j*2*np.pi*xr[i]*ξr)
+ return out / N
+
+def f_d(ξr,ϵ):
+ return cfft(d(xr,ϵ),ξr)
+
+fig3 = plt.figure(3,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlabel(r'$\xi_r$')
+plt.ylabel(r'$| \hat{d_0} |$')
+plt.xlim([ξmin,ξmax])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξr,np.abs(f_d(ξr,vϵ[i])),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+if save: plt.savefig(name+'_fd.eps')
+
+def f_r(ξr,ϵ):
+ return cfft(r(ϵ),ξr)
+
+fig4 = plt.figure(4,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlabel(r'$\xi_r$')
+plt.ylabel(r'$| \hat{r} |$')
+plt.xlim([ξmin,ξmax])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξr,np.abs(f_r(ξr,vϵ[i])),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+if save: plt.savefig(name+'_fr.eps')
+
+"""
+COST FUNCTION
+"""
+ϵ = np.linspace(2*ϵmin,2*ϵmax,100)
+
+#Analytical
+"""
+fig5 = plt.figure(5,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([2*ϵmin,2*ϵmax])
+plt.grid()
+
+plt.plot(ϵ,J_l2(ϵ),'k')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],J_l2(vϵ[i]),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+
+fig6 = plt.figure(6,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.225,left=0.05)
+plt.xlim([2*ϵmin,2*ϵmax])
+plt.grid()
+plt.plot(ϵ,J_sob(ϵ),'k')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],J_sob(vϵ[i]),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+"""
+
+#Numerical
+def σ(ξr):
+ return s*2*np.pi*1j*ξr
+
+def distance1(ϵ):
+ return 0.5*(np.linalg.norm(f_r(ξr,ϵ)))**2
+v_distance1 = np.vectorize(distance1)
+
+def distance2(ϵ):
+ return 0.5*(np.linalg.norm(σ(ξr)*f_r(ξr,ϵ)))**2
+v_distance2 = np.vectorize(distance2)
+
+"""
+COST FUNCTION
+"""
+#Least squares
+fig7 = plt.figure(7,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlabel(r'$\epsilon$')
+plt.xlim([2*ϵmin,2*ϵmax])
+plt.ylim([0,Jmax])
+plt.grid()
+
+#plt.plot(ϵ,v_distance1(ϵ),'k')
+plt.plot(ϵ,J_l2(ϵ),'tab:gray')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],v_distance1(vϵ[i]),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+if save: plt.savefig(name+'_j_l2.eps')
+
+fig8 = plt.figure(8,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlabel(r'$\epsilon$')
+plt.xlim([2*ϵmin,2*ϵmax])
+if name=='example1' or name=='example4_an' or name=='example4_nu':plt.ylim([0,Jmax])
+else: plt.ylim([0,Jmax/10])
+plt.grid()
+
+#plt.plot(ϵ,v_distance2(ϵ),'k')
+plt.plot(ϵ,J_sob(ϵ),'tab:gray')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],v_distance2(vϵ[i]),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+if save: plt.savefig(name+'_j_sob.eps')
+
+plt.close('all')
+
+"""
+SYMBOL MIXING
+"""
+#def α(l):
+ #return 1/(np.pi*l*np.sqrt(f_array(l)))
+#Numerical
+def σ3(ξr,α,s):
+ return (1 + 1j*2*np.pi*α*L*ξr) / (1 + 1j*2*np.pi*s*λ*ξr)
+
+def distance3(ϵ,α,s):
+ return 0.5*(np.linalg.norm(σ3(ξr,α,s)*f_r(ξr,ϵ)))**2
+v_distance3 = np.vectorize(distance3)
+"""
+Nα=21
+αmin=α(L/2/h)/2
+αmax=α(L/2/h)*3/2
+vα = np.linspace(αmin,αmax,Nα)
+vαc = np.linspace(0.95,0.05,Nα)
+
+fig9 = plt.figure(9,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.3,left=0.05)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlim([2*ϵmin,2*ϵmax])
+plt.xlabel(r'$\epsilon$')
+#plt.ylim([0,Jmax])
+plt.grid()
+#plt.plot(ϵ,(α(1)*L)**2*J_sob(ϵ),'tab:gray')
+#plt.plot(ϵ,J_l2(ϵ),'tab:gray')
+plt.plot(ϵ,J_l2(ϵ)+(α(L/2/h)*L)**2*J_sob(ϵ),'tab:gray')
+
+for i in range(0,Nα):
+ plt.plot(ϵ,v_distance3(ϵ,vα[i],0),color=plt.cm.Spectral(vαc[i]))
+ #plt.plot(ϵ,J_l2(ϵ)+vα[i]**2*J_sob(ϵ),'tab:gray')
+plt.plot(ϵ,v_distance3(ϵ,vα[10],0),'k--')
+"""
+def f_array_dise(l):
+ return np.nan_to_num(1/8/l * (np.arctan(l) - l*(1-l**2) / (1+l**2)**2))
+def α(l):
+ return 1/(np.pi*l*np.sqrt(f_array_dise(l)))
+Ns=21
+smin=0
+smax=α(1.)*L/λ
+vs = np.linspace(smin,smax,Ns)
+vsc = np.linspace(0.95,0.05,Ns)
+
+fig10 = plt.figure(10,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.3,left=0.05)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlim([2*ϵmin,2*ϵmax])
+plt.xlabel(r'$\epsilon$')
+plt.ylim([0,50.])
+plt.grid()
+
+v2ϵ = np.linspace(ϵmin,2*ϵmax,2*Nϵ-1)
+for i in range(0,Ns):
+ plt.plot(v2ϵ,v_distance3(v2ϵ,2.*α(1.),vs[i]),color=plt.cm.Spectral(vsc[i]))
+plt.plot(v2ϵ,J_l2(v2ϵ),color='tab:gray')
+plt.plot(v2ϵ,(2.*α(1.)*L)**2*J_sob(v2ϵ),color='tab:gray')
diff --git a/paper3ChapOtherExampleV0.py b/paper3ChapOtherExampleV0.py
new file mode 100644
index 0000000..fdf4aaa
--- /dev/null
+++ b/paper3ChapOtherExampleV0.py
@@ -0,0 +1,263 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Mar 26 08:13:15 2021
+
+@author: xavier
+"""
+import numpy as np
+import scipy as sp
+import scipy.special
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 1.5;
+W = 6.5;
+save = False;
+
+L = 10.
+λ = 1.
+k = 2*np.pi / λ
+N = 2**8+1
+Δxr = L/N
+xr = np.linspace(-(L-Δxr)/2, (L-Δxr)/2, N)
+ξr = (np.arange(0,N)-np.floor(N/2))/L
+s=λ
+
+Nϵ=38
+ϵmax=0.2
+ϵmin=0.
+Δϵ=(ϵmax-ϵmin)/(Nϵ-1)
+ξmin=0
+ξmax=2/λ
+xmin=0
+xmax=L/2
+
+Jmax=1.5
+"""
+Example IV (analitycal): Wave number perturbation (amplitude)
+"""
+amplitude=True
+if amplitude:
+ def d(xr,ϵ):
+ return np.conj(-1j/4*sp.special.hankel1(0,k*(1+ϵ)*np.abs(xr)))
+else:
+ def d(xr,ϵ):
+ return np.exp(1j*k*(1+ϵ)*np.abs(xr))
+
+"""
+Example IV (numerical): Wave number perturbation (amplitude)
+"""
+name = 'example4_nu'
+
+Nd=2*Nϵ-1
+d_num = np.zeros([N,Nd])*1j
+for i in range(0,Nd):
+ tmp = (np.loadtxt('paper3_example4_multiscale_df0n'+str(i)+'0'+'.csv',delimiter=';',skiprows=1,max_rows=1)).view(np.complex128);
+ d_num[:,i] = tmp
+
+def d_an(xr,ϵ):
+ i = (np.around(xr / Δxr)+(N-1)/2).astype(int)
+ j = (np.around((ϵ+ϵmin) / Δϵ)).astype(int)
+ return d_num[i,j]
+"""
+SPATIAL DOMAIN ANALYSIS
+"""
+vϵ = np.linspace(ϵmin,ϵmax,Nϵ)
+vc = np.linspace(0.95,0.05,Nϵ)
+fig1 = plt.figure(1,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlim([xmin,xmax])
+plt.xlabel(r'$x_r$')
+plt.ylabel(r'Re($d_0$)')
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.real(d(xr,vϵ[i])),color=plt.cm.Spectral(vc[i]))
+if save: plt.savefig(name+'_d.eps')
+
+def w(ϵ):
+ return d(xr,ϵ) / d_an(xr,0)
+def r(ϵ):
+ return 1-w(ϵ)
+
+fig2 = plt.figure(2,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlim([xmin,xmax])
+plt.xlabel(r'$x_r$')
+plt.ylabel(r'$\left\vert r \right\vert$')
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(xr,np.abs(r(vϵ[i])),color=plt.cm.Spectral(vc[i]))
+if save: plt.savefig(name+'_r.eps')
+
+"""
+FREQUENCY DOMAIN ANALYSIS
+"""
+Nint = 8
+ξintr = np.linspace(ξr[0],ξr[N-1],N+(N-1)*Nint)
+def cfft(q,ξr):
+ out = np.zeros(np.size(ξr))*1j
+ for i in range(0,N):
+ out += np.nan_to_num(q[i]) * np.exp(-1j*2*np.pi*xr[i]*ξr)
+ return out / N
+
+def f_d(ξr,ϵ):
+ return cfft(d(xr,ϵ),ξr)
+
+fig3 = plt.figure(3,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlabel(r'$\xi_r$')
+plt.ylabel(r'$| \hat{d_0} |$')
+plt.xlim([ξmin,ξmax])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξr,np.abs(f_d(ξr,vϵ[i])),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+if save: plt.savefig(name+'_fd.eps')
+
+def f_r(ξr,ϵ):
+ return cfft(r(ϵ),ξr)
+
+fig4 = plt.figure(4,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlabel(r'$\xi_r$')
+plt.ylabel(r'$| \hat{r} |$')
+plt.xlim([ξmin,ξmax])
+plt.grid()
+for i in range(0,Nϵ):
+ plt.plot(ξr,np.abs(f_r(ξr,vϵ[i])),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+if save: plt.savefig(name+'_fr.eps')
+
+"""
+COST FUNCTION
+"""
+ϵ = np.linspace(2*ϵmin,2*ϵmax,100)
+α = 0.25
+αs = 0.
+#Numerical
+def σa(ξr):
+ return L*2*np.pi*1j*ξr
+def σs(ξr):
+ return 1 / (1+αs*s*2*np.pi*1j*ξr)
+
+def distance1(ϵ):
+ return 0.5*(np.linalg.norm(f_r(ξr,ϵ)*σs(ξr)))**2
+v_distance1 = np.vectorize(distance1)
+
+def distance2(ϵ):
+ return 0.5*(np.linalg.norm(σa(ξr)*σs(ξr)*f_r(ξr,ϵ)))**2
+v_distance2 = np.vectorize(distance2)
+
+
+def J_l2(ϵ):
+ w_an = d_an(xr,ϵ) / d_an(xr,0)
+ r_an = 1-w_an
+ f_ran = cfft(r_an,ξr)
+ return 0.5*(np.linalg.norm(σs(ξr)*f_ran))**2
+vJ_l2 = np.vectorize(J_l2)
+
+def J_sob(ϵ):
+ w_an = d_an(xr,ϵ) / d_an(xr,0)
+ r_an = 1-w_an
+ f_ran = cfft(r_an,ξr)
+ return 0.5*(np.linalg.norm(σa(ξr)*σs(ξr)*f_ran))**2
+vJ_sob = np.vectorize(J_sob)
+
+"""
+COST FUNCTION
+"""
+#Least squares
+fig7 = plt.figure(7,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlabel(r'$\epsilon$')
+plt.xlim([2*ϵmin,2*ϵmax])
+plt.ylim([0,Jmax])
+plt.grid()
+
+plt.plot(ϵ,v_distance1(ϵ),'k')
+plt.plot(ϵ,vJ_l2(ϵ),'tab:gray')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],v_distance1(vϵ[i]),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+if save: plt.savefig(name+'_j_l2.eps')
+
+fig8 = plt.figure(8,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.975,bottom=0.3,left=0.1)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlabel(r'$\epsilon$')
+plt.xlim([2*ϵmin,2*ϵmax])
+plt.ylim([0,10*Jmax])
+plt.grid()
+
+plt.plot(ϵ,v_distance2(ϵ),'k')
+plt.plot(ϵ,vJ_sob(ϵ),'tab:gray')
+for i in range(0,Nϵ):
+ plt.plot(vϵ[i],v_distance2(vϵ[i]),color=plt.cm.Spectral(vc[i]),marker='o',markersize=3.5)
+if save: plt.savefig(name+'_j_sob.eps')
+
+plt.close('all')
+
+"""
+SYMBOL MIXING
+"""
+#def α(l):
+ #return 1/(np.pi*l*np.sqrt(f_array(l)))
+#Numerical
+def σ3(ξr,α,s):
+ return (1 + 1j*2*np.pi*α*L*ξr) / (1 + 1j*2*np.pi*s*λ*ξr)
+
+def distance3(ϵ,α,s):
+ return 0.5*(np.linalg.norm(σ3(ξr,α,s)*f_r(ξr,ϵ)))**2
+v_distance3 = np.vectorize(distance3)
+"""
+Nα=21
+αmin=α(L/2/h)/2
+αmax=α(L/2/h)*3/2
+vα = np.linspace(αmin,αmax,Nα)
+vαc = np.linspace(0.95,0.05,Nα)
+
+fig9 = plt.figure(9,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.3,left=0.05)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlim([2*ϵmin,2*ϵmax])
+plt.xlabel(r'$\epsilon$')
+#plt.ylim([0,Jmax])
+plt.grid()
+#plt.plot(ϵ,(α(1)*L)**2*J_sob(ϵ),'tab:gray')
+#plt.plot(ϵ,J_l2(ϵ),'tab:gray')
+plt.plot(ϵ,J_l2(ϵ)+(α(L/2/h)*L)**2*J_sob(ϵ),'tab:gray')
+
+for i in range(0,Nα):
+ plt.plot(ϵ,v_distance3(ϵ,vα[i],0),color=plt.cm.Spectral(vαc[i]))
+ #plt.plot(ϵ,J_l2(ϵ)+vα[i]**2*J_sob(ϵ),'tab:gray')
+plt.plot(ϵ,v_distance3(ϵ,vα[10],0),'k--')
+"""
+def f_array_dise(l):
+ return np.nan_to_num(1/8/l * (np.arctan(l) - l*(1-l**2) / (1+l**2)**2))
+def α(l):
+ return 1/(np.pi*l*np.sqrt(f_array_dise(l)))
+Ns=21
+smin=0
+smax=α(1.)*L/λ
+vs = np.linspace(smin,smax,Ns)
+vsc = np.linspace(0.95,0.05,Ns)
+
+fig10 = plt.figure(10,figsize=(W,H),tight_layout=False);
+plt.subplots_adjust(top=0.9,right=0.95,bottom=0.3,left=0.05)
+plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+plt.xlim([2*ϵmin,2*ϵmax])
+plt.xlabel(r'$\epsilon$')
+plt.ylim([0,50.])
+plt.grid()
+
+v2ϵ = np.linspace(ϵmin,2*ϵmax,2*Nϵ-1)
+for i in range(0,Ns):
+ plt.plot(v2ϵ,v_distance3(v2ϵ,2.*α(1.),vs[i]),color=plt.cm.Spectral(vsc[i]))
+plt.plot(v2ϵ,J_l2(v2ϵ),color='tab:gray')
+plt.plot(v2ϵ,(2.*α(1.)*L)**2*J_sob(v2ϵ),color='tab:gray')
\ No newline at end of file
diff --git a/specific/data/objective/phasederivative.cpp b/specific/data/objective/phasederivative.cpp
new file mode 100644
index 0000000..471f5a1
--- /dev/null
+++ b/specific/data/objective/phasederivative.cpp
@@ -0,0 +1,157 @@
+//GmshFEM Library
+#include "Exception.h"
+#include "Function.h"
+#include "AlgebraicFunctions.h"
+
+//GmshFWI Library
+#include "phasederivative.h"
+
+#define d0 *_d0
+#define d ds.state(Type::F)
+#define rho ds.state(Type::A)
+#define dd ds.state(Type::PF)
+#define drho ds.state(Type::PA)
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::function;
+using namespace gmshfem::equation;
+using namespace gmshfem::post;
+/*
+* SOBOLEV
+*/
+namespace phasederivative
+{
+ template<Physic T_Physic>
+ void Objective<T_Physic>::updateSystem_a()
+ {
+ if(_systemIsUpToDate_a){return;}
+ _formulation_a.removeSystem();
+ _formulation_a.initSystem();
+
+ //double coef = 1. / _config->data_area();
+ //_formulation_a.galerkin( 1.0 * coef * dof(_v_a), tf(_v_a), _config->data_omega(), "Gauss2");
+
+ TensorFunction< std::complex< double > > Ca = (_l_a) / 2. / M_PI * tensor< std::complex< double > >(1., 0., 0., 0., 1., 0., 0., 0., 0.);
+ _formulation_a.galerkin( Ca * grad(dof(_v_a)), Ca * grad(tf(_v_a)), _config->data_omega(), "Gauss2");
+
+ _formulation_a.pre();
+ _formulation_a.assemble();
+ _systemIsUpToDate_a=true;
+ _formulation_a.removeTerms();
+ }
+
+ template<Physic T_Physic>
+ Data<T_Physic> Objective<T_Physic>::operator_a(const Data<T_Physic>& r)
+ {
+ DataField<T_Physic> rF(_v_a);
+ rF.setValuesToZero();
+ rF.forward_interpolate(r);
+ algebra::Vector< std::complex<double> > rFv;
+ rF.getAllVector(rFv);
+
+ algebra::MatrixCRSFast< std::complex<double> > Ma;
+ if(!_systemIsUpToDate_a){updateSystem_a();}
+ _formulation_a.getLHS(Ma);
+ algebra::Vector< std::complex<double> > rF_av(rFv.size());
+ algebra::linear(rF_av,Ma,rFv);
+
+ DataField<T_Physic> rF_a(_v_a);
+ rF_a.setValuesToZero();
+ rF_a.setAllVector(rF_av);
+
+ Data<T_Physic> out(r);
+ out.value(rF_a);
+ return out / _config->datapoint_area();
+ }
+
+ template<Physic T_Physic>
+ void Objective<T_Physic>::updateSystem_s()
+ {
+ if(_systemIsUpToDate_s){return;}
+ _formulation_s.removeSystem();
+ _formulation_s.initSystem();
+
+ _formulation_s.galerkin( dof(_v_s), tf(_v_s), _config->data_omega(), "Gauss2");
+
+ TensorFunction< std::complex< double > > Cs = (_l_s) / 2. / M_PI * tensor< std::complex< double > >(1., 0., 0., 0., 1., 0., 0., 0., 0.);
+ _formulation_s.galerkin( Cs * grad(dof(_v_s)), Cs * grad(tf(_v_s)), _config->data_omega(), "Gauss2");
+
+ _formulation_s.pre();
+ _formulation_s.assemble();
+ _systemIsUpToDate_s=true;
+ _formulation_s.removeTerms();
+ }
+
+ template<Physic T_Physic>
+ Data<T_Physic> Objective<T_Physic>::inv_operator_s(const Data<T_Physic>& r)
+ {
+ DataField<T_Physic> rF(_v_s);
+ rF.setValuesToZero();
+ rF.forward_interpolate(r);
+
+ if(!_systemIsUpToDate_s){updateSystem_s();}
+ _formulation_s.galerkin( -rF, tf(_v_s), _config->data_omega(), "Gauss2");
+ _formulation_s.assemble();
+ _formulation_s.removeTerms();
+
+ _formulation_s.solve(_systemIsFactorized_s);
+ _systemIsFactorized_s=true;
+ _formulation_s.setRHSToZero();
+ rF = _v_s;
+
+ Data<T_Physic> out(r);
+ out.value(rF);
+ return out;
+ }
+
+ template<Physic T_Physic>
+ double Objective<T_Physic>::performance(const Data<T_Physic>& d1)
+ {
+ Data<T_Physic> r = (d1 - d0) / _w;
+ ( d1 ).write("d1","pos");
+ ( d0 ).write("d0","pos");
+ ( d1 / _w ).write("d1w","pos");
+ ( d0 / _w ).write("d0w","pos");
+ r.write("r","pos");
+ Data<T_Physic> sd1w = inv_operator_s( d1 / _w );
+ Data<T_Physic> sd0w = inv_operator_s( d0 / _w );
+ Data<T_Physic> sr = inv_operator_s(r);
+ sd0w.write("sd0w","pos");
+ sd1w.write("sd1w","pos");
+ (sd1w.phase() / sd0w.phase()).write("sd01w","pos");
+ sr.write("sr","pos");
+ Data<T_Physic> sasr = inv_operator_s(operator_a(sr));
+ sasr.write("sasr","pos");
+ double coef = _config->datapoint_area() / _config->data_area();
+ //msg::print << "coef = " << coef << msg::endl;
+ return 0.5 * coef * std::real( l2scalp(r,sasr) );
+ //return 0.5 * std::real( integrate(_v_s * conj(_v_s), _config->data_omega(), "Gauss2") );
+ }
+
+ template<Physic T_Physic>
+ const Data<T_Physic>& Objective<T_Physic>::update(Type type, const DataStateEvaluator<T_Physic>& ds)
+ {
+ double coef = _config->datapoint_area() / _config->data_area();
+ switch (type)
+ {
+ case Type::F: case Type::PF: default:
+ throw Exception("Objective can not update (perturbed) forward data.");
+ break;
+ case Type::A:
+ _v = conj( inv_operator_s(operator_a(inv_operator_s( (d-d0)/_w ) )) ) / _w * coef;
+ break;
+ case Type::PA:
+ _v = conj( inv_operator_s(operator_a(inv_operator_s( dd/_w ) )) ) / _w * coef;
+ break;
+ case Type::D:
+ /* Very simplified formula */
+ /* Does not take into account user sparsity input */
+ _v.value(PointData<T_Physic>(coef / 4. / 0.02 / 0.02));
+ //_v = _v / ((_d0->phase())*conj((_d0->phase())));
+ break;
+ }
+ return _v;
+ }
+ template class Objective<Physic::acoustic>;
+}
diff --git a/specific/data/objective/phasederivative.h b/specific/data/objective/phasederivative.h
new file mode 100644
index 0000000..8faf882
--- /dev/null
+++ b/specific/data/objective/phasederivative.h
@@ -0,0 +1,74 @@
+#ifndef H_SPECIFIC_DATA_OBJECTIVE_PHASEDERIVATIVE
+#define H_SPECIFIC_DATA_OBJECTIVE_PHASEDERIVATIVE
+
+//GmshFEM Library
+#include "GmshFem.h"
+
+//GmshFWI Library
+#include "../../../common/data/objective/objective.h"
+
+/*
+* SOBOLEV
+*/
+namespace phasederivative
+{
+ template<Physic T_Physic>
+ class Objective final: public ObjectiveInterface<T_Physic>
+ {
+ private:
+ using ObjectiveInterface<T_Physic>::_d0;
+ using ObjectiveInterface<T_Physic>::_v;
+
+ const ConfigurationInterface* const _config;
+ const data::Discretization<T_Physic> _d_discret;
+
+ Data<T_Physic> _w;
+
+ double _l_a;
+ DataField<T_Physic> _v_a;
+ gmshfem::problem::Formulation< std::complex< double > > _formulation_a;
+ bool _systemIsUpToDate_a;
+ void updateSystem_a();
+ Data<T_Physic> operator_a(const Data<T_Physic>& r);
+
+ double _l_s;
+ DataField<T_Physic> _v_s;
+ gmshfem::problem::Formulation< std::complex< double > > _formulation_s;
+ bool _systemIsUpToDate_s;
+ bool _systemIsFactorized_s;
+ void updateSystem_s();
+ Data<T_Physic> inv_operator_s(const Data<T_Physic>& r);
+
+ public:
+ Objective(const Data<T_Physic>* const d0, const ConfigurationInterface* const config, const gmshfem::common::GmshFem& gmshFem)
+ : ObjectiveInterface<T_Physic>(d0), _config(config), _d_discret(),
+ _w(_config), _v_a("_v_a",_config->data_omega(),_config->data_gmodel(),_d_discret), _formulation_a("data_formulation_a"), _systemIsUpToDate_a(false), _v_s("_v_s",_config->data_omega(),_config->data_gmodel(),_d_discret), _formulation_s("data_formulation_s"), _systemIsUpToDate_s(false), _systemIsFactorized_s(false)
+ {
+ std::string path = "";
+ if(!gmshFem.userDefinedParameter(path, "objective_weight"))
+ {
+ gmshfem::msg::warning << "Data weight filename could not be found" << gmshfem::msg::endl;
+ }
+ else
+ {
+ _w.read(path);
+ //_w = _w.amplitude() * _d0->phase();
+ _w.write("_w","pos");
+ }
+ if(!(
+ gmshFem.userDefinedParameter(_l_a, "objective_l_a") &&
+ gmshFem.userDefinedParameter(_l_s, "objective_l_s")
+ ))
+ {
+ throw gmshfem::common::Exception("A data filter parameter could not be found.");
+ }
+ updateSystem_a();
+ updateSystem_s();
+ };
+
+ virtual double performance(const Data<T_Physic>& d);
+ virtual const Data<T_Physic>& update(Type type, const DataStateEvaluator<T_Physic>& ds);
+ };
+};
+
+#endif //H_SPECIFIC_DATA_OBJECTIVE_PHASEDERIVATIVE
diff --git a/specific/data/objective/sobolev.cpp b/specific/data/objective/sobolev.cpp
index c63d10a..33c5c49 100644
--- a/specific/data/objective/sobolev.cpp
+++ b/specific/data/objective/sobolev.cpp
@@ -30,9 +30,9 @@ namespace sobolev
_formulation_a.initSystem();
double coef = 1. / integrate(1.,_config->data_omega(),"Gauss0");
- _formulation_a.galerkin( coef * dof(_v_a), tf(_v_a), _config->data_omega(), "Gauss2");
+ _formulation_a.galerkin( 1.0 * coef * dof(_v_a), tf(_v_a), _config->data_omega(), "Gauss2");
- TensorFunction< std::complex< double > > Ca = (_l_a) * tensor< std::complex< double > >(1., 1., 0., 1., 1., 0., 0., 0., 0.);
+ TensorFunction< std::complex< double > > Ca = (_l_a) * tensor< std::complex< double > >(1., -1., 0., -1., 1., 0., 0., 0., 0.);
_formulation_a.galerkin( coef * Ca * grad(dof(_v_a)), Ca * grad(tf(_v_a)), _config->data_omega(), "Gauss0");
_formulation_a.pre();
@@ -74,7 +74,7 @@ namespace sobolev
_formulation_s.galerkin( dof(_v_s), tf(_v_s), _config->data_omega(), "Gauss2");
- TensorFunction< std::complex< double > > Cs = (_l_s) * tensor< std::complex< double > >(1., 1., 0., 1., 1., 0., 0., 0., 0.);
+ TensorFunction< std::complex< double > > Cs = (_l_s) * tensor< std::complex< double > >(1., -1., 0., -1., 1., 0., 0., 0., 0.);
_formulation_s.galerkin( Cs * grad(dof(_v_s)), Cs * grad(tf(_v_s)), _config->data_omega(), "Gauss0");
_formulation_s.pre();
@@ -108,10 +108,17 @@ namespace sobolev
template<Physic T_Physic>
double Objective<T_Physic>::performance(const Data<T_Physic>& d1)
{
- Data<T_Physic> r = (d1-d0)/d0;
- Data<T_Physic> rs = inv_operator_s(r);
- Data<T_Physic> Mrs = operator_a(rs);
- return 0.5 * std::real( l2scalp(r,Mrs) );
+ _d0->write("d0","pos");
+ d1.write("d1","pos");
+ Data<T_Physic> r = d1.phase() / _d0->phase() + (-1.);
+ r.write("r","pos");
+ Data<T_Physic> sr = inv_operator_s(r);
+ sr.write("sr","pos");
+ Data<T_Physic> asr = operator_a(sr);
+ asr.write("asr","pos");
+ Data<T_Physic> sasr = inv_operator_s(asr);
+ sasr.write("sasr","pos");
+ return 0.5 * std::real( l2scalp(r,sasr) );
}
template<Physic T_Physic>
@@ -123,16 +130,31 @@ namespace sobolev
throw Exception("Objective can not update (perturbed) forward data.");
break;
case Type::A:
- _v = conj( operator_a(inv_operator_s((d-d0)/d0)) ) / d0;
+ _v = conj( inv_operator_s(operator_a(inv_operator_s(
+ d.phase() / _d0->phase() + (-1.)
+ ))) ) / (d.amplitude() * _d0->phase());
+ /*
+ d.write("data0","pos");
+ //( (d-d0) ).write("data1","pos");
+ ( (d)/(_d0->phase()) ).write("data","pos");
+ //( (d-d0)/(_d0->phase()) ).write("data2","pos");
+ //inv_operator_s( (d-d0)/(_d0->phase()) ).write("data3","pos");
+ inv_operator_s(( (d)/(_d0->phase()) )).write("datas","pos");
+ inv_operator_s(inv_operator_s(( (d)/(_d0->phase()) ))).write("datass","pos");
+ */
+ _v.write("_v","pos");
break;
case Type::PA:
- _v = conj( operator_a(inv_operator_s((dd)/d0)) ) / d0;
+ _v = conj( inv_operator_s(operator_a(inv_operator_s(
+ dd.phase() / _d0->phase()
+ ))) ) / (d.amplitude() * _d0->phase());
break;
case Type::D:
/* Very simplified formula */
/* Does not take into account user sparsity input */
- _v.value(PointData<T_Physic>(1.));
- _v = _v / (d0*conj(d0));
+ double coef = (_l_a / _l_s) * (_l_a / _l_s) / integrate(1.,_config->data_omega(),"Gauss0");
+ _v.value(PointData<T_Physic>(coef));
+ //_v = _v / ((_d0->phase())*conj((_d0->phase())));
}
return _v;
}
--
GitLab