2 - linear & non-linear travel time tomography#

Open In Colab


What we do in this notebook#

Here we apply CoFI to two geophysical examples:

  • a linear seismic travel time tomography problem

  • a nonlinear travel time tomography cross borehole problem


Learning outcomes#

  • A demonstration of running CoFI for a regularized linear parameter estimation problem. Can be used as an example of a CoFI template.

  • A demonstration of how a (3rd party) nonlinear forward model can be imported from cofi-espresso and used. Fast Marching algorithm for first arriving raypaths.

  • See how nonlinear iterative matrix solvers can be accessed in CoFI.

# Environment setup (uncomment code below)

# !pip install -U cofi cofi-espresso

Problem description#

The goal in travel-time tomography is to infer details about the velocity structure of a medium, given measurements of the minimum time taken for a wave to propagate from source to receiver.

At first glance, this may seem rather similar to the X-ray tomography problem. However, there is an added complication: as we change our model, the route of the fastest path from source to receiver also changes. Thus, every update we apply to the model will inevitably be (in some sense) based on incorrect assumptions.

Provided the ‘true’ velocity structure is not too dissimilar from our initial guess, travel-time tomography can be treated as a weakly non-linear problem.

In this notebook, we illustrate both linear and one non-linear tomography.

In the first example the straight ray paths are fixed and independent of the medium through which they pass. This would be the case for X-ray tomography, where the data represent amplitude changes across the medium, or seismic tomography under the fixed ray assumption, where the data represent travel times across the medium.

In the second example we iteratively update seismic travel times and ray paths as the seismic velocity model changes, which creates a nonlinear tomographic problem.

In the seismic case, the travel-time of an individual ray can be computed as

\[t = \int_\mathrm{path} \frac{1}{v(\mathbf{x})}\,\mathrm{d}\mathbf{x}\]

This points to an additional complication: even for a fixed path, the relationship between velocities and observations is not linear. However, if we define the ‘slowness’ to be the inverse of velocity, \(s(\mathbf{x}) = v^{-1}(\mathbf{x})\), we can write

\[t = \int_\mathrm{path} {s(\mathbf{x})}\,\mathrm{d}\mathbf{x}\]

which is linear.

We will assume that the object we are interested in is 2-dimensional slowness field. If we discretize this model, with \(N_x\) cells in the \(x\)-direction and \(N_y\) cells in the \(y\)-direction, we can express \(s(\mathbf{x})\) as an \(N_x \times N_y\) vector \(\boldsymbol{s}\).

For the linear case, this is related to the data by

\[d_i = A_{ij}s_j\]

where \(d_i\) is the travel time of the \(i\) th path, and where \(A_{ij}\) represents the path length of raypath \(i\) in cell \(j\) of the discretized model.

For the nonlinear case, this is related to the data by

\[\delta d_i = A_{ij}\delta s_j\]

where \(\delta d_i\) is the difference in travel time, of the \(i\) th path, between the observed time and the travel time in the reference model, and the parameters \(\delta s_j\) are slowness perturbations to the reference model.

import numpy as np
import matplotlib.pyplot as plt

import cofi
import cofi_espresso

1. Linear Travel Time Tomography#

To illustrate the setting we plot a reference model supplied through the espresso Xray example, together with 100 raypaths in the dataset.

linear_tomo_example = cofi_espresso.XrayTomography()
# linear_tomo_example.plot_model(linear_tomo_example.good_model, paths=True);
# linear_tomo_example.plot_model(linear_tomo_example.good_model);
plt.plot(0.5, 0.5, marker="$?$", markersize=130)
for p in linear_tomo_example._paths[:100]:
     plt.plot([p[0],p[2]],[p[1],p[3]],'y',linewidth=0.5)
print(' Data set contains ',len(linear_tomo_example._paths),' ray paths')
2 travel time tomography
Data set contains  10416  ray paths

Step 1. Define CoFI BaseProblem#

Now we: - set up the BaseProblem in CoFI, - supply it the data vector from espresso example, (i.e. the \(\mathbf{d}\) vector) - supply it the Jacobian of the linear system (i.e. the \(A\) matrix)

linear_tomo_problem = cofi.BaseProblem()
linear_tomo_problem.set_data(linear_tomo_example.data)
linear_tomo_problem.set_jacobian(linear_tomo_example.jacobian(linear_tomo_example.starting_model)) # supply matrix A
sigma = 0.1 # set noise level of data
data_cov_inv = np.identity(linear_tomo_example.data_size) * (1/sigma**2)
linear_tomo_problem.set_data_covariance_inv(data_cov_inv)
Evaluating paths:   0%|          | 0/10416 [00:00<?, ?it/s]
Evaluating paths:   8%|7         | 803/10416 [00:00<00:01, 8021.19it/s]
Evaluating paths:  16%|#5        | 1663/10416 [00:00<00:01, 8355.97it/s]
Evaluating paths:  24%|##3       | 2499/10416 [00:00<00:00, 8329.32it/s]
Evaluating paths:  32%|###2      | 3339/10416 [00:00<00:00, 8356.16it/s]
Evaluating paths:  40%|####      | 4190/10416 [00:00<00:00, 8408.91it/s]
Evaluating paths:  48%|####8     | 5031/10416 [00:00<00:00, 8326.70it/s]
Evaluating paths:  56%|#####6    | 5864/10416 [00:00<00:00, 8326.11it/s]
Evaluating paths:  64%|######4   | 6718/10416 [00:00<00:00, 8391.39it/s]
Evaluating paths:  73%|#######2  | 7558/10416 [00:00<00:00, 8374.04it/s]
Evaluating paths:  81%|########  | 8396/10416 [00:01<00:00, 8320.60it/s]
Evaluating paths:  89%|########8 | 9251/10416 [00:01<00:00, 8388.24it/s]
Evaluating paths:  97%|#########6| 10102/10416 [00:01<00:00, 8422.25it/s]
Evaluating paths: 100%|##########| 10416/10416 [00:01<00:00, 8377.28it/s]

Since \(\mathbf{d}\) and \(G\) have been defined then this implies a linear system. Now we choose to regularize the linear system and solve the problem

\[\min_{\mathbf s} \phi({\mathbf d},{\mathbf s}) = ({\mathbf d} - A {\mathbf s})^T C_d^{-1} ({\mathbf d} - A {\mathbf s})~ + ~ \lambda ~{\mathbf s}D^TD{\mathbf s}\]

The matrix system we are solving is

\[(\mathbf{A}^T \textbf{C}_d^{-1} \textbf{A} + \lambda \mathbf D^T\mathbf D) \textbf{s} = \textbf{A}^T \mathbf C_d^{-1} \textbf{d}\]
# set up regularization
lamda = 0.5 # choose regularization constant
data_cov_inv = np.identity(linear_tomo_example.data_size) * (1/sigma**2)
reg_matrix = np.identity(linear_tomo_example.model_size)
linear_tomo_problem.set_regularization(2, lamda, reg_matrix) # choose a regularization (see help(linear_problem.set_regularization) for more details)
print('Number of slowness parameters to be solved for = ',linear_tomo_example.model_size)
Number of slowness parameters to be solved for =  2500

and lets print a summary of the set up.

linear_tomo_problem.summary()
=====================================================================
Summary for inversion problem: BaseProblem
=====================================================================
Model shape: Unknown
---------------------------------------------------------------------
List of functions/properties set by you:
['jacobian', 'regularization', 'regularization_matrix', 'regularization_factor', 'data', 'data_covariance_inv']
---------------------------------------------------------------------
List of functions/properties created based on what you have provided:
['jacobian_times_vector']
---------------------------------------------------------------------
List of functions/properties that can be further set for the problem:
( not all of these may be relevant to your inversion workflow )
['objective', 'log_posterior', 'log_posterior_with_blobs', 'log_likelihood', 'log_prior', 'gradient', 'hessian', 'hessian_times_vector', 'residual', 'jacobian_times_vector', 'data_misfit', 'forward', 'data_covariance', 'initial_model', 'model_shape', 'blobs_dtype', 'bounds', 'constraints']

Step 2. Define CoFI InversionOptions#

Here we choose the backend tool for solving the tomographic system, which is scipy’s least squares solver.

tomo_options = cofi.InversionOptions()
tomo_options.set_tool("scipy.linalg.lstsq")

Step 3. Define CoFI Inversion and run#

tomo_inv = cofi.Inversion(linear_tomo_problem, tomo_options)
tomo_inv_result = tomo_inv.run()
tomo_inv_result.summary()
============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [1.08201552 0.91575052 0.99919151 ... 0.98897406 0.91690785 1.09668824]
sum_of_squared_residuals: []
effective_rank: 2500
singular_values: [372.55549274 343.55222637 343.55222637 ...   0.9576611    0.85184016
   0.85184016]
model_covariance: [[ 2.24308939e-01 -1.34567926e-01 -1.05287422e-02 ...  1.52075868e-04
  -9.65471891e-05 -2.22977182e-04]
 [-1.34567926e-01  3.81542540e-01 -9.10263972e-02 ... -7.13173002e-04
   6.45887554e-04 -9.65471891e-05]
 [-1.05287422e-02 -9.10263972e-02  2.64924830e-01 ...  5.56496611e-04
  -7.13173002e-04  1.52075868e-04]
 ...
 [ 1.52075868e-04 -7.13173002e-04  5.56496611e-04 ...  2.64924830e-01
  -9.10263972e-02 -1.05287422e-02]
 [-9.65471891e-05  6.45887554e-04 -7.13173002e-04 ... -9.10263972e-02
   3.81542540e-01 -1.34567926e-01]
 [-2.22977182e-04 -9.65471891e-05  1.52075868e-04 ... -1.05287422e-02
  -1.34567926e-01  2.24308939e-01]]

Lets plot the image to see what we got.

ax = linear_tomo_example.plot_model(tomo_inv_result.model);
2 travel time tomography

Challenge: Fewer ray paths for linear travel time#

Try and construct a tomographic solution with fewer ray paths.

Here we use 10416 ray paths with indices 0,10415. Try a different range and see what you get.

How many ray paths do you need before the image becomes recognizable?

Upload to Jamboard 1

Start from the code template below:

# data range
idx_from, idx_to = (<CHANGE ME>, <CHANGE ME>)

# basic settings
d = linear_tomo_example.data
G = linear_tomo_example.jacobian(linear_tomo_example.starting_model)

# now attach all the info to a BaseProblem object
mytomo = cofi.BaseProblem()
mytomo.set_data(d[idx_from:idx_to])
mytomo.set_jacobian(G[idx_from:idx_to,:])

# run your problem (with the same InversionOptions) again
mytomo_inv = cofi.Inversion(mytomo, tomo_options)
mytomo_result = mytomo_inv.run()

# check result
fig = linear_tomo_example.plot_model(mytomo_result.model)
plt.title(f'Recovered model from range ({idx_from}, {idx_to})')
plt.figure()
plt.title(' Raypaths')
for p in linear_tomo_example._paths[idx_from:idx_to]:
    plt.plot([p[0],p[2]],[p[1],p[3]],'y',linewidth=0.05)
# Copy the template above, Replace <CHANGE ME> with your answer
#@title Solution

# data range
idx_from, idx_to = (0, 3000)                    # TODO try a different range

# basic settings
d = linear_tomo_example.data
G = linear_tomo_example.jacobian(linear_tomo_example.starting_model)

# now attach all the info to a BaseProblem object
mytomo = cofi.BaseProblem()
mytomo.set_data(d[idx_from:idx_to])
mytomo.set_jacobian(G[idx_from:idx_to,:])

# run your problem (with the same InversionOptions) again
mytomo_inv = cofi.Inversion(mytomo, tomo_options)
mytomo_result = mytomo_inv.run()

# check result
fig = linear_tomo_example.plot_model(mytomo_result.model)
plt.title(f'Recovered model from range ({idx_from}, {idx_to})')
plt.figure()
plt.title(' Raypaths')
for p in linear_tomo_example._paths[idx_from:idx_to]:
    plt.plot([p[0],p[2]],[p[1],p[3]],'y',linewidth=0.05)
  • Recovered model from range (0, 3000)
  • Raypaths
Evaluating paths:   0%|          | 0/10416 [00:00<?, ?it/s]
Evaluating paths:   8%|8         | 836/10416 [00:00<00:01, 8358.95it/s]
Evaluating paths:  17%|#6        | 1720/10416 [00:00<00:01, 8637.10it/s]
Evaluating paths:  25%|##4       | 2584/10416 [00:00<00:00, 8602.80it/s]
Evaluating paths:  33%|###3      | 3449/10416 [00:00<00:00, 8619.38it/s]
Evaluating paths:  41%|####1     | 4322/10416 [00:00<00:00, 8657.66it/s]
Evaluating paths:  50%|####9     | 5188/10416 [00:00<00:00, 8632.27it/s]
Evaluating paths:  58%|#####8    | 6057/10416 [00:00<00:00, 8648.50it/s]
Evaluating paths:  67%|######6   | 6936/10416 [00:00<00:00, 8690.47it/s]
Evaluating paths:  75%|#######4  | 7806/10416 [00:00<00:00, 8625.44it/s]
Evaluating paths:  83%|########3 | 8669/10416 [00:01<00:00, 8622.15it/s]
Evaluating paths:  92%|#########1| 9545/10416 [00:01<00:00, 8660.61it/s]
Evaluating paths: 100%|##########| 10416/10416 [00:01<00:00, 8651.44it/s]

2. Non-linear Travel Time Tomography#

Now we demonstrate CoFI on a nonlinear iterative tomographic problem in a cross borehole setting.

We use a different tomographic example from espresso. Here we import the example module and plot the reference seismic model.

nonlinear_tomo_example = cofi_espresso.FmmTomography()

nonlinear_tomo_example.plot_model(nonlinear_tomo_example.good_model, with_paths=True,lw=0.5);
2 travel time tomography
 New data set has:
 10  receivers
 10  sources
 100  travel times
 Range of travel times:  0.008911182496368759 0.0153757024856463
 Mean travel time: 0.01085811731230709

<Figure size 600x600 with 2 Axes>

Solving the tomographic system with optimization#

Now we solve the tomographic system of equations using either CoFI’s optimization method interface, or its iterative matrix-solver interface.

For the optimization interface:

We choose an objective function of the form.

\[\phi(\mathbf{d},\mathbf{s}) = \frac{1}{\sigma^2}|| \mathbf{d} - \mathbf{g}(\mathbf{s})||_2^2 + \lambda_1 ||\mathbf{s}- \mathbf{s}_{0}||_2^2 + \lambda_2 ||D~\mathbf{s}||_2^2\]

where \(\mathbf{g}(\mathbf{s})\) represents the predicted travel times in the slowness model \(\mathbf{s}\), \(\sigma^2\) is the noise variance on the travel times, \((\lambda_1,\lambda_2)\) are weights of damping and smoothing regularization terms respectively, \(\mathbf{s}_{0}\) is the reference slowness model provided by the espresso example, and \(D\) is a second derivative finite difference stencil for the slowness model with shape model_shape.

In the set up below this objective function is defined outside of CoFI in the function objective_func together with its gradient and Hessian, gradient and hessian with respect to slowness parameters. For convenience the regularization terms are constructed with CoFI utility routine QuadraticReg.

For the optimization case CoFI passes objective_func and optionally the gradient and Hessian functions to a thrid party optimization backend tool such as scipy.minimize to produce a solution.

For the iterative matrix solver interface:

For convenience, CoFI also has its own Gauss-Newton Solver for optimization of a general objective function of the form.

\[\phi(\mathbf{d},\mathbf{s}) = \psi((\mathbf{d},\mathbf{s}) + \sum_{r=1}^R \lambda_r \chi_r(\mathbf{s}),\]

where \(\psi\) represents a data misfit term, and \(\chi_r\) one or more regularization terms, with weights \(\lambda_r\). The objective function above is a special case of this. In general an iterative Gauss-Newton solver takes the form

\[\mathbf{s}_{k+1} = \mathbf{s}_{k} - \cal{H}^{-1}(\mathbf{s}_k) \nabla \phi(\mathbf{s}_k), \quad {(k=0,1,\dots)},\]

where \(\cal{H}(\mathbf{s}_k)\) is the Hessian of the objective function, and \(\nabla \phi(\mathbf{s}_k)\) its gradient evaluated at the model \(\mathbf{s}_k\).

For the objective function above this becomes the simple iterative matrix solver

\[\mathbf{s}_{k+1} = \mathbf{s}_k + (A^T C_d^{-1}A + \lambda_2\mathbf{I} +\lambda_2D^TD )^{-1} [A^T C_d^{-1} (\mathbf{d} - g(\mathbf{s}_k)) - \lambda_2 (\mathbf{s - s}_{0}) - \lambda_2 D^TD \mathbf{s}], \quad (k=0,1,\dots)\]

with \(C_d^{-1} = \sigma^{-2} I\).

Step 1. Define CoFI BaseProblem#

# get problem information from  espresso FmmTomography
model_size = nonlinear_tomo_example.model_size               # number of model parameters
model_shape = nonlinear_tomo_example.model_shape             # 2D spatial grid shape
data_size = nonlinear_tomo_example.data_size                 # number of data points
ref_start_slowness = nonlinear_tomo_example.starting_model   # use the starting guess supplied by the espresso example

Here we define the baseproblem object and a starting velocity model guess.

# define CoFI BaseProblem
nonlinear_problem = cofi.BaseProblem()
nonlinear_problem.set_initial_model(ref_start_slowness)

Here we define regularization of the tomographic system.

# add regularization: damping / flattening / smoothing
damping_factor = 50
smoothing_factor = 5e3
reg_damping = cofi.utils.QuadraticReg(damping_factor, model_size, "damping", ref_start_slowness)
reg_smoothing = cofi.utils.QuadraticReg(smoothing_factor, model_shape, "smoothing")
reg = reg_damping + reg_smoothing
def objective_func(slowness, reg, sigma, data_subset=None):
    if data_subset is None:
        data_subset = np.arange(0, nonlinear_tomo_example.data_size)
    ttimes = nonlinear_tomo_example.forward(slowness)
    residual = nonlinear_tomo_example.data[data_subset] - ttimes[data_subset]
    data_misfit = residual.T @ residual / sigma**2
    model_reg = reg(slowness)
    return  data_misfit + model_reg

def gradient(slowness, reg, sigma, data_subset=None):
    if data_subset is None:
        data_subset = np.arange(0, nonlinear_tomo_example.data_size)
    ttimes, A = nonlinear_tomo_example.forward(slowness, with_jacobian=True)
    ttimes = ttimes[data_subset]
    A = A[data_subset]
    data_misfit_grad = -2 * A.T @ (nonlinear_tomo_example.data[data_subset] - ttimes) / sigma**2
    model_reg_grad = reg.gradient(slowness)
    return  data_misfit_grad + model_reg_grad

def hessian(slowness, reg, sigma, data_subset=None):
    if data_subset is None:
        data_subset = np.arange(0, nonlinear_tomo_example.data_size)
    A = nonlinear_tomo_example.jacobian(slowness)[data_subset]
    data_misfit_hess = 2 * A.T @ A / sigma**2
    model_reg_hess = reg.hessian(slowness)
    return data_misfit_hess + model_reg_hess
sigma = 0.00001                   # Noise is 1.0E-4 is ~5% of standard deviation of initial travel time residuals

nonlinear_problem.set_objective(objective_func, args=[reg, sigma, None])
nonlinear_problem.set_gradient(gradient, args=[reg, sigma, None])
nonlinear_problem.set_hessian(hessian, args=[reg, sigma, None])

Step 2. Define CoFI InversionOptions#

nonlinear_options = cofi.InversionOptions()

# cofi's own simple newton's matrix-based optimization solver
nonlinear_options.set_tool("cofi.simple_newton")
nonlinear_options.set_params(num_iterations=5, step_length=1, verbose=True)

# scipy's Newton-CG solver (alternative approach with similar results)
# nonlinear_options.set_tool("scipy.optimize.minimize")
# nonlinear_options.set_params(method="Newton-CG", options={"xtol":1e-16})
nonlinear_options.summary()
=============================
Summary for inversion options
=============================
Solving method: None set
Use `suggest_solving_methods()` to check available solving methods.
-----------------------------
Backend tool: `cofi.simple_newton` - CoFI's own solver - simple Newton's approach (for testing mainly)
References: ['https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization']
Use `suggest_tools()` to check available backend tools.
-----------------------------
Solver-specific parameters:
num_iterations = 5
step_length = 1
verbose = True
Use `suggest_solver_params()` to check required/optional solver-specific parameters.

Step 3. Define CoFI Inversion and run#

nonlinear_inv = cofi.Inversion(nonlinear_problem, nonlinear_options)
nonlinear_inv_result = nonlinear_inv.run()
nonlinear_tomo_example.plot_model(nonlinear_inv_result.model);
2 travel time tomography
Iteration #0, objective function value: 110298.7001724638
Iteration #1, objective function value: 1787.1051514815424
Iteration #2, objective function value: 121.1495033985667
Iteration #3, objective function value: 5.814222496115815
Iteration #4, objective function value: 4.086694560516768

<Figure size 600x600 with 2 Axes>

Now lets plot the true model for comparison.

nonlinear_tomo_example.plot_model(nonlinear_tomo_example.good_model);
2 travel time tomography
<Figure size 600x600 with 2 Axes>

Challenge: Change the number of tomographic data#

First try and repeat this tomographic reconstruction with fewer data and see what the model looks like.

There are 100 raypaths in the full dataset and you can tell CoFI to select a subset by passing an additional array of indices to the functions that calculate objective, gradient and hessian.

Upload to Jamboard 2

Start from the code template below:

# Set a subset of raypaths here
data_subset = np.arange(<CHANGE ME>)

# select BaseProblem
my_own_nonlinear_problem = cofi.BaseProblem()
my_own_nonlinear_problem.set_objective(objective_func, args=[reg, sigma, data_subset])
my_own_nonlinear_problem.set_gradient(gradient, args=[reg, sigma, data_subset])
my_own_nonlinear_problem.set_hessian(hessian, args=[reg, sigma, data_subset])
my_own_nonlinear_problem.set_initial_model(ref_start_slowness)

# run inversion with same options as previously
my_own_inversion = cofi.Inversion(my_own_nonlinear_problem, nonlinear_options)
my_own_result = my_own_inversion.run()

# check results
my_own_result.summary()

# plot inverted model
fig, paths = nonlinear_tomo_example.plot_model(my_own_result.model, return_paths=True)
print(f"Number of paths used: {len(data_subset)}")

# plot paths used
for p in np.array(paths, dtype=object)[data_subset]:
    fig.axes[0].plot(p[:,0], p[:,1], "g", alpha=0.5,lw=0.5)
# Copy the template above, Replace <CHANGE ME> with your answer
#@title Solution

# Set a subset of raypaths here
data_subset = np.arange(30, 60)

# select BaseProblem
my_own_nonlinear_problem = cofi.BaseProblem()
my_own_nonlinear_problem.set_objective(objective_func, args=[reg, sigma, data_subset])
my_own_nonlinear_problem.set_gradient(gradient, args=[reg, sigma, data_subset])
my_own_nonlinear_problem.set_hessian(hessian, args=[reg, sigma, data_subset])
my_own_nonlinear_problem.set_initial_model(ref_start_slowness)

# run inversion with same options as previously
my_own_inversion = cofi.Inversion(my_own_nonlinear_problem, nonlinear_options)
my_own_result = my_own_inversion.run()

# check results
my_own_result.summary()

# plot inverted model
fig, paths = nonlinear_tomo_example.plot_model(my_own_result.model, return_paths=True)
print(f"Number of paths used: {len(data_subset)}")

# plot paths used
for p in np.array(paths, dtype=object)[data_subset]:
    fig.axes[0].plot(p[:,0], p[:,1], "g", alpha=0.5,lw=0.5)
2 travel time tomography
Iteration #0, objective function value: 33804.23569764485
Iteration #1, objective function value: 133.41002612375848
Iteration #2, objective function value: 2.5725408628470254
Iteration #3, objective function value: 0.4155647466194026
Iteration #4, objective function value: 0.00449669101050086
============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [0.00050057 0.00050052 0.00050046 ... 0.00051289 0.00051088 0.00050873]
num_iterations: 4
objective_val: 0.00048736044767036367
n_obj_evaluations: 5
n_grad_evaluations: 5
n_hess_evaluations: 5
Number of paths used: 30

Challenge: Change regularization settings#

In the solution above we used damping_factor = 50, and smoothing_factor = 5.0E-3 and flattening_factor = 0.

Experiment with these choices, e.g increasing all of them to say 100 and repeat the tomographic solution to see how the model changes.

Try to turn off smoothing all together but retain damping and flattening and see what happens.

With some choices you can force an under-determined problem which is not solvable.

(Note that here we revert back to using all of the data by removing the data_subset argument to the objective function.)

To repeat this solver with other settings for smoothing and damping strength. See the documentation for cofi.utils.QuadraticReg.

Upload to Jamboard 3

You can start from the template below:

# change the combination of damping, flattening and smoothing regularizations
damping_factor = <CHANGE ME>                # select damping factor here to force solution toward reference slowness model
flattening_factor = <CHANGE ME>             # increase flattening factor here to force small first derivatives in slowness solution
smoothing_factor = <CHANGE ME>              # increase smoothing factor here to force small second derivatives in slowness solution

reg_damping = cofi.utils.QuadraticReg(damping_factor, model_size, "damping", ref_start_slowness)
reg_flattening = cofi.utils.QuadraticReg(flattening_factor, model_shape, "flattening")
reg_smoothing = cofi.utils.QuadraticReg(smoothing_factor, model_shape, "smoothing")
my_own_reg = reg_damping + reg_flattening + reg_smoothing

# set Baseproblem
my_own_nonlinear_problem = cofi.BaseProblem()
my_own_nonlinear_problem.set_objective(objective_func, args=[my_own_reg, sigma, None])
my_own_nonlinear_problem.set_gradient(gradient, args=[my_own_reg, sigma, None])
my_own_nonlinear_problem.set_hessian(hessian, args=[my_own_reg, sigma, None])
my_own_nonlinear_problem.set_initial_model(ref_start_slowness.copy())

# run inversion with same options as previously
my_own_inversion = cofi.Inversion(my_own_nonlinear_problem, nonlinear_options)
my_own_result = my_own_inversion.run()

# check results
fig = nonlinear_tomo_example.plot_model(my_own_result.model)
fig.suptitle(f"Damping {damping_factor}, Flattening {flattening_factor}, Smoothing {smoothing_factor}");
# Copy the template above, Replace <CHANGE ME> with your answer
#@title Reference Solution

# change the combination of damping, flattening and smoothing regularizations
damping_factor = 100                # select damping factor here to force solution toward reference slowness model
flattening_factor = 100             # increase flattening factor here to force small first derivatives in slowness solution
smoothing_factor = 0                # increase smoothing factor here to force small second derivatives in slowness solution

reg_damping = cofi.utils.QuadraticReg(damping_factor, model_size, "damping", ref_start_slowness)
reg_flattening = cofi.utils.QuadraticReg(flattening_factor, model_shape, "flattening")
reg_smoothing = cofi.utils.QuadraticReg(smoothing_factor, model_shape, "smoothing")
my_own_reg = reg_damping + reg_flattening + reg_smoothing

# set Baseproblem
my_own_nonlinear_problem = cofi.BaseProblem()
my_own_nonlinear_problem.set_objective(objective_func, args=[my_own_reg, sigma, None])
my_own_nonlinear_problem.set_gradient(gradient, args=[my_own_reg, sigma, None])
my_own_nonlinear_problem.set_hessian(hessian, args=[my_own_reg, sigma, None])
my_own_nonlinear_problem.set_initial_model(ref_start_slowness.copy())

# run inversion with same options as previously
my_own_inversion = cofi.Inversion(my_own_nonlinear_problem, nonlinear_options)
my_own_result = my_own_inversion.run()

# check results
fig = nonlinear_tomo_example.plot_model(my_own_result.model)
fig.suptitle(f"Damping {damping_factor}, Flattening {flattening_factor}, Smoothing {smoothing_factor}");
Damping 100, Flattening 100, Smoothing 0
Iteration #0, objective function value: 110298.7001724638
Iteration #1, objective function value: 17978.18529996893
Iteration #2, objective function value: 20381.924100291555
Iteration #3, objective function value: 45743931762530.02
Iteration #4, objective function value: 8.246670839512789e+20

Text(0.5, 0.98, 'Damping 100, Flattening 100, Smoothing 0')

Watermark#

watermark_list = ["cofi", "cofi_espresso", "numpy", "scipy", "matplotlib"]
for pkg in watermark_list:
    pkg_var = __import__(pkg)
    print(pkg, getattr(pkg_var, "__version__"))
cofi 0.1.2.dev22
cofi_espresso 0.0.1.dev10
numpy 1.21.6
scipy 1.9.1
matplotlib 3.5.3

sphinx_gallery_thumbnail_number = -1

Total running time of the script: ( 0 minutes 30.438 seconds)

Gallery generated by Sphinx-Gallery