continuous l2 distance

common/data/objective/objective.h

@@ -32,6 +32,8 @@ public:
     virtual void link(ModelUpdater* const mu, DataUpdater<T_Physic>* const du) {};
     virtual void unlink() {};
     virtual bool modelIsObsolete() {return false;};//Returns true, if objective depends on model
+
+    virtual void scale(double scale) {};
 };
 
 #endif //H_COMMON_DATA_OBJECTIVE

common/data/state.h

@@ -21,7 +21,7 @@ class DataStateEvaluator
 {
 public:
     virtual const Data<T_Physic>& state(Type type) const = 0;
-    void write(Type type, std::string name) const {state(type).write(name);};
+    void write(Type type, std::string name, std::string fmt = "csv") const {state(type).write(name,fmt);};
 };
 
 /*

common/statefunctional.cpp

@@ -24,7 +24,7 @@ void StateFunctional<T_Physic>::write_model(Type type, std::string name)
 }
 
 template<Physic T_Physic>
-void StateFunctional<T_Physic>::write_data(Type type, std::string name)
+void StateFunctional<T_Physic>::write_data(Type type, std::string name, std::string fmt)
 {
     std::array<bool,5> DataNeedToBeUpToDate = {false,false,false,false,false};
     DataNeedToBeUpToDate[type] = true;
@@ -37,7 +37,7 @@ void StateFunctional<T_Physic>::write_data(Type type, std::string name)
         case Type::F: case Type::A: case Type::D:
             ModelNeedToBeUpToDate[Type::F] = true;
     }
-    _du.get(DataNeedToBeUpToDate,_mu.get(ModelNeedToBeUpToDate)).write(type,name);
+    _du.get(DataNeedToBeUpToDate,_mu.get(ModelNeedToBeUpToDate)).write(type,name,fmt);
 }
 
 template<Physic T_Physic>
@@ -68,6 +68,22 @@ void StateFunctional<T_Physic>::innerproductIsObsolete()
     _mu.isObsolete(NoMoreUpToDate);
 }
 
+template<Physic T_Physic>
+void StateFunctional<T_Physic>::objectiveIsObsolete()
+{
+    std::array<bool,5> dataNoMoreUpToDate = {false, false, true, true, true};
+    _du.isObsolete(dataNoMoreUpToDate);
+
+    std::array<bool,4> waveNoMoreUpToDate = {false, false, true, true};
+    _wu.isObsolete(waveNoMoreUpToDate);
+
+    std::array<bool,4> modelNoMoreUpToDate = {false, false, true, true};
+    _mu.isObsolete(modelNoMoreUpToDate);
+
+    std::array<bool,3> sensitivityNoMoreUpToDate = {true, true, true};
+    _su.isObsolete(sensitivityNoMoreUpToDate);
+}
+
 template<Physic T_Physic>
 void StateFunctional<T_Physic>::modelIsObsolete()
 {

common/statefunctional.h

@@ -31,9 +31,10 @@ public:
     ~StateFunctional();
 
     void innerproductIsObsolete();
+    void objectiveIsObsolete();
 
     void write_model(Type type, std::string name);
-    void write_data(Type type, std::string name);
+    void write_data(Type type, std::string name, std::string fmt = "csv");
     void write_wave(Type type, std::string name);
 
     virtual const InnerProductInterface* const innerproduct() const {return _innerproduct;};

input/paper1/reference.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May  6 15:15:24 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 = 15.
+f = 1.
+ω = 2. * np.pi * f
+v = 1.
+k = ω / v
+N = 11
+r0 = 15.
+x = np.linspace(0, L, N)
+
+def d(xe,xr):
+        r = np.sqrt( (xe-xr)**2+r0**2 )
+        return np.conj( -1j/4*sp.special.hankel1(0,k*r) )
+
+D = np.zeros([N, 2*N])
+for e in range(0,N):
+    for r in range(0,N):
+        D[e,2*r] = np.real(d(x[e],x[r]))
+        D[e,2*r+1] = np.imag(d(x[e],x[r]))
+
+np.savetxt("paper1_reference_highly.csv",D,delimiter=";")

input/paper3/common.txt

-prename = paper2_
+prename = paper3_
 #Configuration
 configuration = soil
 #

input/paper3/convergence/cluster_inversion_nscgA.sh

+#!/bin/bash
+#
+#SBATCH --job-name=paper3_inversion_nscgA
+#SBATCH --output=paper3_inversion_nscgA.log
+
+#SBATCH --ntasks=1
+#SBATCH --time=2400:00
+#SBATCH --mem-per-cpu=48000
+
+time ./inversion  ../input/paper3/common.txt ../input/paper3/convergence/inversion.txt ../input/paper3/convergence/inversion_nscgA.txt -verbose 2 -maxThreads 1

input/paper3/convergence/cluster_inversion_nscgB.sh

+#!/bin/bash
+#
+#SBATCH --job-name=paper3_inversion_nscgB
+#SBATCH --output=paper3_inversion_nscgB.log
+
+#SBATCH --ntasks=1
+#SBATCH --time=2400:00
+#SBATCH --mem-per-cpu=48000
+
+time ./inversion  ../input/paper3/common.txt ../input/paper3/convergence/inversion.txt ../input/paper3/convergence/inversion_nscgB.txt -verbose 2 -maxThreads 1

input/paper3/convergence/cluster_inversion_nscgC.sh

+#!/bin/bash
+#
+#SBATCH --job-name=paper3_inversion_nscgC
+#SBATCH --output=paper3_inversion_nscgC.log
+
+#SBATCH --ntasks=1
+#SBATCH --time=2400:00
+#SBATCH --mem-per-cpu=48000
+
+time ./inversion  ../input/paper3/common.txt ../input/paper3/convergence/inversion.txt ../input/paper3/convergence/inversion_nscgC.txt -verbose 2 -maxThreads 1

input/paper3/convergence/cluster_syntetics.sh

+#!/bin/bash
+#
+#SBATCH --job-name=paper3_synthetics
+#SBATCH --output=paper3_synthetics.log
+
+#SBATCH --ntasks=1
+#SBATCH --time=30:00
+#SBATCH --mem-per-cpu=48000
+
+time ./synthetics ../input/paper3/common.txt ../input/paper3/synthetics.txt -verbose 2 -maxThreads 1

input/paper3/convergence/inversion.txt

+#Configuration
+unknown = soil
+m0_type = velocity_file
+m0_scale = 2.0
+path = ../input/ZP2018Marmousi2D/marmousi
+#
+#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
+#
+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

input/paper3/convergence/inversion_nscgA.txt

+name = inversion_nscgA
+#
+#Local minimum search
+localminimumsearch = yuanfan
+localminimum_initMu = 1
+localminimum_c0 = 0.0001
+localminimum_c1 = 0.25
+localminimum_c2 = 0.25
+localminimum_c3 = 0.2
+localminimum_c4 = 5.
+localminimum_writeInterval = 5
+localminimum_maxIteration = 50
+#
+#Descent search
+descentsearch = newton_steihaug_conjugategradient
+descent_maxIteration = 50
+descent_eta0 = 0.4

input/paper3/convergence/inversion_nscgB.txt

+name = inversion_nscgB
+#
+#Local minimum search
+localminimumsearch = yuanfan
+localminimum_initMu = 1
+localminimum_c0 = 0.0001
+localminimum_c1 = 0.75
+localminimum_c2 = 0.75
+localminimum_c3 = 0.25
+localminimum_c4 = 2.
+localminimum_writeInterval = 5
+localminimum_maxIteration = 50
+#
+#Descent search
+descentsearch = newton_steihaug_conjugategradient
+descent_maxIteration = 50
+descent_eta0 = 0.4

input/paper3/convergence/inversion_nscgC.txt

+name = inversion_nscgC
+#
+#Local minimum search
+localminimumsearch = yuanfan
+localminimum_initMu = 1
+localminimum_c0 = 0.0001
+localminimum_c1 = 0.9
+localminimum_c2 = 0.9
+localminimum_c3 = 0.5
+localminimum_c4 = 2.
+localminimum_writeInterval = 5
+localminimum_maxIteration = 50
+#
+#Descent search
+descentsearch = newton_steihaug_conjugategradient
+descent_maxIteration = 50
+descent_eta0 = 0.4

input/paper3/ip_comparison.txt

@@ -11,20 +11,19 @@ Im(v_H) = 0.
 gaussNewton=0
 #Objective
 objective = continuousl2distance
-objective_weight = paper2_reference_data0
-objective_l_s = 1.5
+objective_weight0 = ../input/paper3/paper3_reference_data0
 #
 #Inner product
 complex = 0
 boundaryWeight = 0.
 diag_preconditioner = 1
 innerproduct = sobolev
-ref_preconditioner = 50
+ref_preconditioner = 250
 #
-scale0 = 1.5
+scale0 = 0.8
 #
 #Data
-data0 = paper2_synthetics_data0
+data0 = paper3_synthetics_data0
 #
 #Discretization
 #Wave

input/paper3/reference.py

@@ -23,21 +23,20 @@ save = False;
 Example I: Wave number perturbation
 """
 L = 8.712
-f = 4.
+f = np.array([4., 6., 8.])
 ω = 2. * np.pi * f
 v = 1.5
 k = ω / v
 N = 122
 x = np.linspace(0, L, N)
 
-def d(xe,xr):
-        if(xe==xr): return 0.5 * (d(0.,1e-3)+d(0.,-1e-3)) #Must match numerical value at 0 (?)
-        return np.conj( -1j/4*sp.special.hankel1(0,k*np.abs(xe-xr)) )
-
+def d(xe,xr,i):
+        if(xe==xr): return 0.5 * (d(0.,1e-3,i)+d(0.,-1e-3,i)) #Must match numerical value at 0 (?)
+        return np.conj( -1j/4*sp.special.hankel1(0,k[i]*np.abs(xe-xr)) )
+for i in range(0,np.size(f)):
     D = np.zeros([N, 2*N])
     for e in range(0,N):
         for r in range(0,N):
-        D[e,2*r] = np.real(d(x[e],x[r]))
-        D[e,2*r+1] = np.imag(d(x[e],x[r]))
-
-np.savetxt("paper2_reference_data0.csv",D,delimiter=";")
\ No newline at end of file
+            D[e,2*r] = np.real(d(x[e],x[r],i))
+            D[e,2*r+1] = np.imag(d(x[e],x[r],i))
+    np.savetxt("paper3_reference_data"+str(i)+".csv",D,delimiter=";")
\ No newline at end of file