<jupyter_start><jupyter_text>Example of 3D time-domain FWI with a constant density modelIn this example, we will show how to use PySIT to conduct a 3D FWI with a constant density model. The corresponding .py file can be found in ``/Demo/InversionHorizontalReflector3D_TimeDomain_ConstantDensity.py``<jupyter_code>%matplotlib inline<jupyter_output><empty_output><jupyter_text>Import necessary modules:<jupyter_code>import time
import copy

import numpy as np
import matplotlib.pyplot as plt
import math
import os
from shutil import copy2
from mpl_toolkits.axes_grid1 import make_axes_locatable

import sys
import scipy.io as sio

from pysit import *
from pysit.gallery import horizontal_reflector
from pysit.util.io import *
from pysit.vis.vis import *
from pysit.util.compute_tools import *

from pysit.util.parallel import *<jupyter_output><empty_output><jupyter_text>Define the physical domain, computational mesh and velocity models 1. Define perfectly matched layer(PML) boundaries in x, y and z directions with width of ``0.1 km`` and PML coefficient of ``100`` by    ``pmlx = PML(0.1, 100, compact=True)`` (the argument ``compact``)        ``pmly = PML(0.1, 100, compact=True)``    ``pmlz = PML(0.1, 100, compact=True)``        For more information about the PML object, we refer users to check ``/pysit/core/domain.py`` by2. Define a 3D rectangular domain with the ranges of ``(0.1, 1.0) km``, ``(0.1, 0.9) km``, ``(0.1, 0.8) km`` in x, y, and z directions.        ``x_config = (0.1, 1.0, pmlx, pmlx)`` (The physical domain of x direction starts at 0.1 km and ends at 1.0 km.)        ``y_config = (0.1, 0.9, pmly, pmly)`` (The physical domain of y direction starts at 0.1 km and ends at 0.9 km.)        ``z_config = (0.1, 0.8, pmlz, pmlz)`` (The physical domain of z direction starts at 0.1 km and ends at 0.8 km.)        ``d = RectangularDomain(x_config, y_config, z_config)``        For more information about the RectangularDomain, we refer users to check ``/pysit/core/domain.py``.3. Define the computational Cartesian mesh with ``46`` grids in ``x`` direction, ``41`` grids in ``y`` direction, and ``36`` grids in ``z`` direction by        ``m = CartesianMesh(d, 46, 41, 36)`` (The computational mesh ``m`` meets the physical domain ``d``)        For more information about the CartesianMesh object, we refer users to check ``/pysit/core/mesh.py``.    4. Generate the true velocity model and initial model for a Horizontal reflector model by    ``C, C0, m, d = horizontal_reflector(m)``        The output ``C`` is the true velocity model and ``C0`` is the initial model.        For more information about the horizontal_reflecotr object, we refer users to check ``/pysit/gallery/horizontal_reflector.py``.<jupyter_code>pmlx = PML(0.1, 100)
    pmly = PML(0.1, 100)
    pmlz = PML(0.1, 100)

    x_config = (0.1, 1.0, pmlx, pmlx)
    y_config = (0.1, 0.9, pmly, pmly)
    z_config = (0.1, 0.8, pmlz, pmlz)

    d = RectangularDomain(x_config, y_config, z_config)

    m = CartesianMesh(d, 46, 41, 36)
    
    C, C0, m, d = horizontal_reflector(m)
    
    n_data = (46, 41, 36)
    n_dataplt = (n_data[0], n_data[2], n_data[1])
    Cplot = np.reshape(C, n_data)
    Cplot = np.transpose(Cplot, (0, 2, 1))
    origins = [0.1, 0.1, 0.1]
    deltas = [0.02, 0.02, 0.02]
    
    axis_ticks = [np.array(list(range(0, n_dataplt[0]-5, (n_data[0]-6)//4))), 
                  np.array(list(range(5, n_dataplt[1]-5, (n_data[1]-11)//4))),
                  np.array(list(range(0, n_dataplt[2], (n_data[2]-1)//2)))
                 ]
    axis_tickslabels = [(axis_ticks[0] * deltas[0] * 1000.0 + origins[0] * 1000.0).astype(int),
                        (axis_ticks[1] * deltas[1] * 1000.0 + origins[1] * 1000.0).astype(int),
                        (axis_ticks[2] * deltas[2] * 1000.0 + origins[2] * 1000.0).astype(int)
                       ]

    
    plot_3D_panel(Cplot, slice3d=(22, 18, 20),
                  axis_label=['x [m]', 'z [m]', 'y [m]'],
                  axis_ticks=axis_ticks,
                  axis_tickslabels=axis_tickslabels,
                  vmin=1.6, vmax=2.2
                 )
    plt.title('Slice at \n x = 540 m, y = 500 m, z = 440 m')<jupyter_output><empty_output><jupyter_text>Set up shots1. Set up the shots object by:    ``shots = equispaced_acquisition(m, RickerWavelet(10.0), sources=Nshots, source_depth=zpos,source_kwargs={},receivers='max',receiver_depth=zpos,receiver_kwargs={})``    ``equispaced_acquisition`` - create a shots object with equially spaced acquisition    ``m`` - computational mesh    ``RickerWavelet(10.0)`` - a Ricker wavelet centered at ``10 Hz``    ``sources`` - number of sources, for 3D modeling, sources has two elements that indicates number of sources in x and y directions. For example, ``sources = (2,3)`` means that there are 2 shots in the x direction and 3 shots in the y direction.    ``source_depth`` - the depth of sources    ``receivers`` - number of receivers, if set to be ``max``, then the number of receivers equals to the number of grids in x direction.    ``receiver_depth`` - the depth of receivers        For more information about the ``equispaced_acquisition`` object, we refer the users to check ``/pysit/core/acquisition.py``.     2. Set up the range of recording time by;    ``trange = (0.0,3.0)``.<jupyter_code># Set up shots
    zmin = d.z.lbound
    zmax = d.z.rbound
    zpos = zmin + (1./9.)*zmax

    Nshots = 1,1

    shots = equispaced_acquisition(m,
                                   RickerWavelet(10.0),
                                   sources=Nshots,
                                   source_depth=zpos,
                                   source_kwargs={},
                                   receivers='max',
                                   receiver_depth=zpos,
                                   receiver_kwargs={}
                                   )

    
    shots_freq = copy.deepcopy(shots)


    # Define and configure the wave solver
    trange = (0.0,3.0)<jupyter_output><class 'numpy.float64'>
1.0<jupyter_text>Define the wave-equation solver and the computational model object, and generate time-domain data 1. In this example, we use the time-domain constant density acoustic wave as our target wave equation. We set up our wave equation solver by:    ``solver = ConstantDensityAcousticWave(m, spatial_accuracy_order=2, trange=trange, kernel_implementation='cpp')``                                         ``m`` - the computational mesh    ``spatial_accuracy_order`` - the spatial accuray order for the numerical solver. Users can select one of the four values ``2, 4, 6, 8``.     ``trange`` - the range of the recording time    ``kernel_implementattion`` - the implementation of the stencil kernel. When set it to be 'cpp', we use the stencil implemented in the language of ``C++``.         For more information about the ``ConstantDensityAcousticWave`` object, we refer the users to check ``/pysit/solvers/wave_factory.py``    2. Create the velocity model object for the wave-quation solver by:    ``base_model = solver.ModelParameters(m,{'C': C})``        The model object ``base_model`` contains the information of the computational mesh and the velocity model ``C``.    3. Generate the time-domain data by:    ``generate_seismic_data(shots, solver, base_model)``        The generated data are stored in the object ``shots``. In order to check the data of the $i^{\text{th}}$ shot, you may need to use the command:        ``data = shots[i].receivers.data``        For more information about the ``generate_seismic_data`` function, we refer the users to check ``/pysit/modeling/data_modeling.py``.<jupyter_code>solver = ConstantDensityAcousticWave(m,
                                         spatial_accuracy_order=2,
                                         trange=trange,
                                         kernel_implementation='cpp',
                                         max_C=3.2)
    solver.max_C = 3.2
    
    base_model = solver.ModelParameters(m,{'C': C})
    
    generate_seismic_data(shots, solver, base_model)
    
    data = shots[0].receivers.data
    
    t_smp = np.linspace(trange[0], trange[1], data.shape[0])
    
    fig=plt.figure()
    n_recdata = [len(t_smp), n_data[0], n_data[1]]
    n_recdataplt = [n_data[0], len(t_smp), n_data[1]]
    data = np.reshape(data, n_recdata)
    dataplt = np.transpose(data, (1, 0, 2))
    deltas_data = [deltas[0], solver.dt,  deltas[2]]
    origins_data = [origins[0], 0.0,origins[2]]
    
    axis_ticks = [np.array(list(range(0, n_recdataplt[0]-5, (n_recdataplt[0]-1)//4))), 
                  np.array(list(range(0, n_recdataplt[1]-5, (n_recdataplt[1]-1)//4))),
                  np.array(list(range(0, n_recdataplt[2], (n_recdataplt[2]-1)//2)))
                 ]
    axis_tickslabels = [np.round(axis_ticks[0] * deltas_data[0]  + origins_data[0], 2),
                        np.round(axis_ticks[1] * deltas_data[1]  + origins_data[1], 2),
                        np.round(axis_ticks[2] * deltas_data[2]  + origins_data[2], 2)
                       ]

    
    plot_3D_panel(dataplt, slice3d=(22, 900, 20),
                  axis_label=[ 'x [km]', 'Time [s]', 'y [km]'],
                  axis_ticks=axis_ticks,
                  axis_tickslabels=axis_tickslabels,
                  width_ratios=[1,1], height_ratios=[1,1],cmap='seismic', vmin=-0.2,vmax=0.2
                 )<jupyter_output><empty_output><jupyter_text>Inversion After the data generation, we will show how to set up the inversion. 1. **Set up the objective function object.** The first thing for the inversion is to set up the objective function. In this version of PySIT, we provide the following objective function options: ``Least-squares method``, ``Envelope method``, ``Cross-correlation method``, ``Optimal transportation method``, and ``FWI with extended imaging``. In order to use these objective functions, users can follow the follwoing instruction:    ** *Least squares*:** for the least-squares method, the necessary input is the solver of the wave-equation. There are three optional inputs:         ``filter_op`` - the filtering oprerator used in the inversion. We will define one filtering operator in the following example.        ``parallel_wrap_shot`` - the parallel wrapper for using parallel computing. The default is ``ParallelWrapShotNull``, which means that we do not use parallel computing.        ``imaging_period`` - imaging period for computing the gradient. The default is ``1``, which means that we use each step of the wavefield to compute the gradient.        To generate a Least-squares objective function, we can use the following command:        ``objective = TemporalLeastSquares(solver, filter_op=filter_op1, parallel_wrap_shot=pwrap, imaging_period=1)``,        where ``filter_op1`` and ``pwrap`` are predefined filtering operator and parallel wrapper, respectively.        ** *Envelope method*:** for envelope method, the input arguments are almost the same as the Least-squares method. There is one one additional optional argument:        ``envelope_power`` – the power for the envelope data. By default, we set it to be ``2``.        To generate an objective function for the envelope method, we can use the following command:        ``objective = TemporalEnvelope(solver, envelope_power=2.0, filter_op=filter_op1, parallel_wrap_shot=pwrap, imaging_period=1)``.        ** *Cross-correlation method*:** for the cross-correlation method, the input arguments are the same as the Least-squares method. To generate an objective function for the cross-correlation method, we can use the following command:        ``objective = TemporalCorrelate(solver, filter_op=filter_op, parallel_wrap_shot=pwrap, imaging_period=1)``.        ** *Optimal transportation method*:** for the optimal transportation method, we provide four different tranformation options to transform the seismic data to a probability density function, i.e. ``Linear function``, ``Quadratic function``, ``Absolute function`` and ``exponential function``. We can set the input arguement ``transform_mode`` to select the tranformation function. There are four options for ``transform_mode`` - ``linear``, ``quadratic``, ``absolute``, and ``exponential``.         We can use the following commands to generate the objective function for each transformation option:        [1] Linear function:         ``objective = TemporalOptimalTransport(solver, filter_op=filter_op1, parallel_wrap_shot=pwrap, imaging_period=1, transform_mode='linear', c_ratio=2.0)``,        where ``c_ratio`` corresponds to the parameter $c$ in the linear function $f(x) = x + c$.        [2] Quadratic function:        ``objective = TemporalOptimalTransport(solver, filter_op=filter_op1, parallel_wrap_shot=pwrap, imaging_period=1, transform_mode='quadratic')``.        [3] Absolute function:        ``objective = TemporalOptimalTransport(solver, filter_op=filter_op1, parallel_wrap_shot=pwrap, imaging_period=1, transform_mode='absolute')``.        [4] Exponential function:        ``objective = TemporalOptimalTransport(solver, filter_op=filter_op1, parallel_wrap_shot=pwrap, imaging_period=1, transform_mode='exponential', exp_a=1.0)``.        where the input argument ``exp_a`` correspons to the $a$ in the exponential function $f(x)=\exp(ax)$.          For more information about the objective functions, we recommend users to check the files in ``/pysit/objective_functions``.    2. **Define the filtering operator.**    In PySIT, we provide the feature to define a band-pass filter to conduct band-pass filtering on the data. An example of generating a band-pass filter is as follows:        ``filter_op = band_pass_filter(n_timesmp, T_max, freq_band, transit_freq_length=0.5, padding_zeros=True, nl=500, nr=500)``,        where        ``n_timesmp`` - number of the sampling points in time.        ``T_max`` - the maximum recording time.         ``freq_band`` - freq_band = (min_freq, max_freq), which defines the lowest and highest frequency of the band-pass filter.        ``transit_freq_length`` - the width of the transition band that allows the transition between the passband and the stopband of the signal in frequency domain. Here, we set the transition band width to be ``0.5 Hz``.        ``padding_zeros`` - the flag to decide whether we need to pad zeros at the both sides of the signal.        ``nl, nr`` - number of zeros padded at the left side (nl) and right side (nr) of the signal.        For more information about the band_pass filter, we recommend users to check the file ``/pysit/util/compute_tool.py``.3. **Define the parallel wrapper.**    PySIT provides the feature to compute different shots in parallel. In oreder to activate this feature, we should define a parallel wrapper. We can define the parallel wrapper by the following command:        ``pwrap = ParallelWrapShot()``.        For more information about the parallel wrapper, we recommend users to check the file ``/pysit/util/parallel.py``.        Now, let us first define a Least-squares objective function as an example.<jupyter_code>pwrap = ParallelWrapShot()

    n_timesmp = shots[0].receivers.data.shape[0]
    T_max = solver.tf
    freq_band =[1.0, 30.0]
    filter_op1 = band_pass_filter(n_timesmp, T_max, freq_band, transit_freq_length=0.5, padding_zeros=True, nl=500, nr=500)
    
    objective = TemporalLeastSquares(solver, filter_op=filter_op1, parallel_wrap_shot=pwrap, imaging_period=1)<jupyter_output><empty_output><jupyter_text>4. **Define the optimization object.**    After generating the objective function object, the next step is to generate the optimziation object. PySIT provides the following optimization algorithms:     . projected gradient descent method,         . projected quasi newton method,         . Gauss-Newton method.         We can use the following command to generate the optimization object corresponding to the above algorithms.         [1]. Projected gradient descent method:        ``invalg = PQN(objective, proj_op=Proj_Op)``,        where the input argument ``proj_op`` is the projection operator that we decide to use in the inversion. It is an optional input argument. By default, we set ``proj_op=None``, which means we do not use any projection operator and the resulting algorithm is the basic gradient descent method. We will introduce how to generate a projcetion operator in the next cell.        [2]. Projected quasi newton method:        ``invalg = PQN(objective, proj_op=Proj_Op, memory_length=10)``,        where the additional input argument ``memory_length`` denotes the number of previous iterations stored in the memory for l-BFGS method. Same as the projected gradient descent method, the default setting of ``proj_op`` is ``proj_op=None``.        [3]. Gauss-Newton method:        ``invalg = GaussNewton(objective, krylov_maxiter=50)``,        where the additional input argument ``krylov_maxiter`` defines the maximal numbal of the krylov iterations used to solve the linearized problem.        For More information about the optimization methods, we recommend users to check the directory ``/pysit/optimization``.   5. **Define the projection operator.**    In the current version of PySIT, we provide two projection operators – ``BoxConstraintPrj`` and ``WaterVelocityPrj``. To generate the two operator, users can use the following commands:        [1] Box Constraint projection:            ``Proj_op1 = BoxConstraintPrj(bound)``,         where         ``bound = (v_min, v_max)`` denotes the lower and upper bound of the velocity.        [2] Water Velocity Projection:        ``Proj_op2 = WaterVelocityPrj(ModelSize, NumberofWaterLayer=10, WaterVel=1.5)``,        where         ``ModelsSize = [nx, nz]`` denotes the size of the model,        ``NumberofWaterLayer`` denotes the number of the water layer,        ``WaterVel`` denotes the velocity of water.        We also provide the feature to combine two projection operators. Users can achieve this by the following command:        ``Proj_joint = JointPrj(Proj_op1, Proj_op2)``.        For more information about the projection operator, we recommend users to check the file ``/pysit/optimization/projection.py``.    6. **Define the number of iterations.**    We can define the number of iterations by ``nsteps = 5``.    7. **Set up status_configuration.**    Please check the following python cell to check the setting of the status_configuration.8. **Set up the line search type**.    We recommend to use the backtrack line search by the following command:         ``line_search = 'backtrack'``.    9. **Set up the initial model**.    We can set up the initial model by:        ``initial_value = solver.ModelParameters(m,{'C': C0})``.9. **Run the inversion**.    After setting up all the necessary and optional arguments, we can conduct the inversion by:        result = invalg(shots, initial_value, nsteps,                    line_search=line_search,                    status_configuration=status_configuration, verbose=True, write=False)                        where, the argument ``write`` denotes whether the user wants to write the model of each iteration to the disk.        Now, let us use PQN as an example.<jupyter_code>bound = [1.5, 3.0]
    Proj_Op1 = BoxConstraintPrj(bound)
    invalg = PQN(objective, proj_op=Proj_Op1, memory_length=10)
    
    nsteps = 5
    status_configuration = {'value_frequency'           : 1,
                            'residual_frequency'        : 1,
                            'residual_length_frequency' : 1,
                            'objective_frequency'       : 1,
                            'step_frequency'            : 1,
                            'step_length_frequency'     : 1,
                            'gradient_frequency'        : 1,
                            'gradient_length_frequency' : 1,
                            'run_time_frequency'        : 1,
                            'alpha_frequency'           : 1,
                            }
    
    initial_value = solver.ModelParameters(m,{'C': C0})
    line_search = 'backtrack'
    
    result = invalg(shots, initial_value, nsteps,
                    line_search=line_search,
                    status_configuration=status_configuration, verbose=True)<jupyter_output>Iteration 0
  gradnorm 2232.6124379152484
  objective 0.014326661791442086
  residual 0.16927292631393887
  Starting:  1.0 0.014326661791442086
  Pass 1: a:1.0; 0.013453910632222187 ?<= 0.014325730993494247
  alpha 1.0
  run time 175.14640545845032s
Iteration 1
  gradnorm 1998.5017931049276
  objective 0.013453910632222187
  residual 0.16403603648114756
  Starting:  1.0 0.013453910632222187
  Pass 1: a:1.0; 0.010957860448157778 ?<= 0.013447894976484792
  alpha 1.0
  run time 156.66591334342957s
Iteration 2
  gradnorm 4278.970190511999
  objective 0.010957860448157778
  residual 0.14803959232690272
  Starting:  1.0 0.010957860448157778
  Pass 1: a:1.0; 0.009927805611951966 ?<= 0.010954863061180396
  alpha 1.0
  run time 149.93397426605225s
Iteration 3
  gradnorm 4188.404423114504
  objective 0.009927805611951966
  residual 0.1409099401174521
  Starting:  1.0 0.009927805611951966
  Pass 1: a:1.0; 0.009358755520577145 ?<= 0.009927027277658458
  alpha 1.0
  run time 156.214346408844s
Itera[...]<jupyter_text>Check resultWe can check the history of the objective values by the following command:``obj_vals = np.array([v for k,v in list(invalg.objective_history.items())])``.We can check the final velocity model by the following command:``v_final = result.C``.   Now, let us check the final results.<jupyter_code>obj_vals = np.array([v for k,v in list(invalg.objective_history.items())])

    plt.figure()
    plt.semilogy(obj_vals)
    plt.xlabel('Iteration')
    plt.ylabel('Objective value')
    
    Cplot = np.reshape(C, n_data)
    Cplot = np.transpose(Cplot, (0, 2, 1))
    origins = [0.1, 0.1, 0.1]
    deltas = [0.02, 0.02, 0.02]
    
    axis_ticks = [np.array(list(range(0, n_dataplt[0]-5, (n_data[0]-6)//4))), 
                  np.array(list(range(5, n_dataplt[1]-5, (n_data[1]-11)//4))),
                  np.array(list(range(0, n_dataplt[2], (n_data[2]-1)//2)))
                 ]
    axis_tickslabels = [(axis_ticks[0] * deltas[0] * 1000.0 + origins[0] * 1000.0).astype(int),
                        (axis_ticks[1] * deltas[1] * 1000.0 + origins[1] * 1000.0).astype(int),
                        (axis_ticks[2] * deltas[2] * 1000.0 + origins[2] * 1000.0).astype(int)
                       ]

    
    plot_3D_panel(Cplot, slice3d=(22, 18, 20),
                  axis_label=['x [m]', 'z [m]', 'y [m]'],
                  axis_ticks=axis_ticks,
                  axis_tickslabels=axis_tickslabels,
                  vmin=1.6, vmax=2.2
                 )
    plt.title('True model Slice at \n x = 540 m, y = 500 m, z = 440 m')
    
    Cout = result.C
    
    Coutplot = np.reshape(Cout, n_data)
    Coutplot = np.transpose(Coutplot, (0, 2, 1))
    origins = [0.1, 0.1, 0.1]
    deltas = [0.02, 0.02, 0.02]
    
    axis_ticks = [np.array(list(range(0, n_dataplt[0]-5, (n_data[0]-6)//4))), 
                  np.array(list(range(5, n_dataplt[1]-5, (n_data[1]-11)//4))),
                  np.array(list(range(0, n_dataplt[2], (n_data[2]-1)//2)))
                 ]
    axis_tickslabels = [(axis_ticks[0] * deltas[0] * 1000.0 + origins[0] * 1000.0).astype(int),
                        (axis_ticks[1] * deltas[1] * 1000.0 + origins[1] * 1000.0).astype(int),
                        (axis_ticks[2] * deltas[2] * 1000.0 + origins[2] * 1000.0).astype(int)
                       ]

    
    plot_3D_panel(Coutplot, slice3d=(22, 18, 20),
                  axis_label=['x [m]', 'z [m]', 'y [m]'],
                  axis_ticks=axis_ticks,
                  axis_tickslabels=axis_tickslabels,
                  vmin=1.6, vmax=2.2
                 )
    plt.title('Reconstructed model Slice at \n x = 540 m, y = 500 m, z = 440 m')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10996362420035957, 'fraction_numerical': 0.048626499979094365, 'mean_word_length': 2.956003969566656, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '332411', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 8148, 'n_tokens_neox': 7320, 'n_words': 2335, 'path': 'Notebook/.ipynb_checkpoints/InversionHorizontalReflector3D_TimeDomain_ConstantDensity-checkpoint.ipynb', 'repo_name': 'zfang-slim/PysitForPython3', 'size': '23917'}