Skip to content

GitLab

Commit 4d84dd04 authored by Robert Clapp's avatar Robert Clapp
Browse files

initial commit

parent 76f25cc8
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 32682 additions and 0 deletions
+32682 -0
Causal-integration.ipynb 0 → 100644
+ 2637
- 0
View file @ 4d84dd04
This diff is collapsed.
Examples.ipynb 0 → 100644
+ 2824
- 0
View file @ 4d84dd04
This source diff could not be displayed because it is too large. You can view the blob instead.
GEE-chapter1.ipynb 0 → 100644
+ 452
- 0
View file @ 4d84dd04
%% Cell type:markdown id: tags:
# Basic operators and adjoints
A great many of the calculations
we do in science and engineering
are really matrix multiplication in disguise.
The first goal of this chapter is to unmask the disguise
by showing many examples.
Second, we see how the
**adjoint** operator (matrix transpose)
back projects information from data to the underlying model.
Geophysical modeling calculations
generally use linear operators that predict data from models.
Our usual task is to find the inverse of these calculations;
i.e., to find models (or make images) from the data.
Logically, the adjoint is the first step
and a part of all subsequent steps in this \bx{inversion} process.
Surprisingly, in practice, the adjoint sometimes does a better job
than the inverse!
Better because the adjoint operator tolerates imperfections in the data
and does not demand the data provide full information.
If you will permit me a poet's license with words,
I will offer you the following table
of **operator**s and their **adjoint**s
<table>
<tr>
<th>matrix multiply </th><th>conjugate-transpose matrix multiply</th></tr>
<tr>
<th>convolve </th><th>crosscorrelate </th></tr>
<tr>
<th>truncate </th><th>zero pad </th></tr>
<tr>
<th>replicate, scatter, spray </th><th>sum or stack </th></tr>
<tr>
<th>spray into neighborhoods </th><th>sum within bins </th></tr>
<tr>
<th>derivative (slope) </th><th>negative derivative </th></tr>
<tr>
<th>causal integration </th><th>anticausal integration |
<tr>
<th>add functions </th><th>do integrals </th></tr>
<tr>
<th>assignment statements </th><th>added terms </th></tr>
<tr>
<th>plane-wave superposition </th><th>slant stack / beam form |
<tr>
<th>spread on a curve </th><th>sum along a curve </th></tr>
<tr>
<th>stretch </th><th>squeeze </th></tr>
<tr>
<th>scalar field gradient </th><th>negative of vector field divergence</th></tr>
<tr>
<th>upward continue </th><th>downward continue</th></tr>
<tr>
<th>diffraction modeling </th><th>imaging by migration</th></tr>
<tr>
<th>hyperbola modeling </th><th>stacking for image or velocity </th></tr>
<tr>
<th>chop image into overlapping patches </th><th>merge the patches</th></tr>
<tr>
<th>ray tracing </th><th>tomography </th></tr>
</table>
%% Cell type:markdown id: tags:
The left column is often called **modeling**,
and the adjoint operators in the right column are often
used in data **processing**.
When the adjoint operator is
*not*
an adequate approximation to the inverse,
then you apply the techniques of fitting and optimization
explained in Chapter 2.
These techniques require iterative use of the
modeling operator and its adjoint.
The adjoint operator is sometimes called
the ``**back projection**'' operator,
because information propagated in one direction (Earth to data) is projected
backward (data to Earth model).
Using complex-valued operators,
the transpose and complex conjugate go together;
and in **Fourier analysis**, taking the complex conjugate
of $\exp(i\omega t)$ reverses the sense of time.
With more poetic license, I say that adjoint operators
*undo*
the time and therefore,
phase shifts of modeling operators.
The inverse operator does also,
but it also divides out the color.
For example, when linear interpolation is done,
then high frequencies are smoothed out; inverse interpolation must restore the
m.
You can imagine the possibilities for noise amplification
which is why adjoints are safer than inverses.
But, nature determines in each application what is the best operator to use
and whether to stop after the adjoint,
to go the whole way to the inverse,
or to stop partway.
%% Cell type:markdown id: tags:
The operators and adjoints previously shown transform vectors to other vectors
.
They also transform data planes to model planes, volumes, etc.
A mathematical operator transforms an "abstract vector" that
might be packed full of volumes of information like television
signals (time series) can pack together a movie, a sequence of frames.
We can always think of the operator as being a matrix,
but the matrix can be truly huge (and nearly empty).
When the vectors transformed by the matrices are large like
geophysical data set sizes,
then the matrix sizes are "large squared,"
far too big for computers.
Thus,
although we can always think of an operator as a matrix;
in practice, we handle an operator differently.
Each practical application requires the practitioner to
prepare two computer programs.
One performs the matrix multiply
$\bf y =\bf B \bf x$,
while the other multiplies by the transpose
$\tilde{\bf x} =\bf B^T \bf y$
(without ever having the matrix itself in memory).
It is always easy to transpose a matrix.
It is less easy to take a computer program that does
$\bf y =\bf B \bf x$
and convert it to another to do
$\tilde{\bf x} =\bf B^T \bf y$,
which is what we'll be doing here.
In this chapter are many examples of increasing complexity.
At the end of the chapter,
we see a test for any program pair
to see whether the operators $\bf B$ and $\bf B^T$ are mutually adjoint
as they should be.
Doing the job correctly (coding adjoints without making approximations)
rewards us later when we tackle model and image-estimation applications.
Mathematicians often denote the transpose of a matrix
$\bf B$ by $\bf B^{\rm T}$.
In physics and engineering, we often encounter complex numbers.
There, the adjoint is the complex-conjugate transposed matrix
denoted $\bf B^\ast$.
%% Cell type:markdown id: tags:
# Object orriented design
We need to define some abstract types before we proceed. The first is a vector. In addition to a list of values, a vector lives in some vector space. The space tells what these list of numbers correspond to. For example: how many samples, regular or irregular, how many dimensions, sampling of the dimensions, etc? Our abstract class needs to be able tell if another vector is from the same space. In addtion it need to make a copy of it space. We will add one more function that is able to zero itself.
%% Cell type:code id: tags:
``` python
# Vector class and derived classes
class vector:
"""Abstract python vector class"""
def __init__(self):
"""Default constructor"""
pass;
def zero(self):
"""Function to zero out a vector"""
raise NotImplementedError("zero must be overwritten")
def clone(self):
"""Function to clone (deep copy) a vector from a vector or a Space"""
raise NotImplementedError("clone must be overwritten")
def cloneSpace(self):
"""Function to clone vector space"""
raise NotImplementedError("cloneSpace must be overwritten")
def checkSame(self):
"""Function to check to make sure the vectors exist in the same space"""
raise NotImplementedError("checkSame must be overwritten")
```
%% Cell type:markdown id: tags:
An operator can be thought of as a matrix that transforms from its domain to its range using its forward and from range to its domain using its adjoint. An operator must be able to store its domain and range and check to make sure the model and data passed into it (when doing the forward and adjoint) match the domain and range it was initialized with,
%% Cell type:code id: tags:
``` python
class Operator:
"""Abstract python operator class"""
# Default class methods/functions
def __init__(self, domain, range):
"""Generic class for operator"""
self.domain = domain.cloneSpace()
self.range = range.cloneSpace()
def forward(self,add,model,data):
"""Apply forward of operator"""
raise NotImplementedError("forward must be overwritten")
def adjoint(self,add,model,data):
"""Apply adjoint of operator"""
raise NotImplementedError("adjoint must be overwritten")
def forward(self, add, model, data):
"""Forward operator"""
raise NotImplementedError("Forward must be defined")
def adjoint(self, add, model, data):
"""Adjoint operator"""
raise NotImplementedError("Adjoint must be defined")
```
%% Cell type:markdown id: tags:
We will make a concrete example of our vector class using numpy ndArrays.
%% Cell type:code id: tags:
``` python
import numpy
class numpyVector(vector):
def __init__(self,**kw):
if "dims" in kw:
self.v=numpy.ndArray(kw[dims])
self.shape=self.v.shape
self.spaceOnly=False
elif "array" in kw:
self.v=kw["array"]
self.shape=self.v.shape
self.spaceOnly=False
elif "space" in kw:
self.shape=kw["space"]
self.spaceOnly=True;
else:
raise ValueError("Only know how to create a numpyVector with dims=")
def zero(self):
"""Zero arrray"""
self.v.fill(0.)
def cloneSpace(self):
"""Return a numpy vector is just shape"""
return numpyVector(space=self.shape)
def clone(self):
"""Return copy of vector"""
return numpyVector(array=self.v.copy())
def checkSame(self,vec):
"""CHeck to make sure if vector in same space"""
if isinstance(vec,numpyVector):
if len(vec.shape)== len(self.shape):
for i in range(len(self.shape)):
if self.shape[i] != vec.shape[i]:
raise ValueError("Axis %d - %d is not same as %d "%(i,self.shape[i],vec.shape[i]))
else:
raise ValueError("Dimensions not the same %d %d"%(len(self.shape),len(vec.shape)))
else:
raise ValueError("Not a numpyVector")
return True
```
%% Cell type:markdown id: tags:
## Programming linear operators
The operation $ y_i = \sum_j b_{ij} x_j$ is the multiplication
of a matrix $\bf B$ by a vector $\bf x$.
The adjoint operation is
$\tilde x_j = \sum_i b_{ij} y_i$.
The operation adjoint to multiplication by a matrix
is multiplication by the transposed matrix
(unless the matrix has complex elements,
in which case,
we need the complex-conjugated transpose).
The following code does matrix multiplcation $\bf y = \bf B \bf x$ and
multiplication by the transpose $\tilde{ \bf x} = \bf B^T \bf y$.
We will define a python class that is initialized by a matrix. We will define
a forward and adjoint operator for the class.
%% Cell type:code id: tags:
``` python
class matmult(pyOperator.Operator):
"""Operator that does 1-D interpolation
"""
def __init__(self, matrix, model, data):
"""Initialize operator
matrix - 2-D array
model - Model space
data - Data space
Both must be the same space, derived from pyVector.vectorIC"""
if not isinstance(model, giee.vector):
raise Exception(
"Model must begiee.vector or derived from it")
if not isinstance(data, giee.vector):
raise Exception(
"Model must be giee.vector or derived from it")
if not isinstance(matrix, np.ndarray):
raise Excecption(
"Matrix must be ndarray")
if model.ndims != 1:
raise Exception("Expecting model to be a 1-D array")
if data.ndims != 1:
raise Exception("Expecting data to be a 1-D array")
self._matrix = matrix
if self._matrix.ndim != 2:
raise Exception("Expecting matrix to be 2-D")
if self._matrix.shape[0] != model.getHyper().axes[0].n:
raise Exception(
"Expecting matrix first axis equal to length of model")
if self._matrix.shape[1] != data.getHyper().axes[0].n:
raise Exception(
"Expecting matrix second axis equal to length of data")
self.setDomainRange(model, data)
def forward(self, add, model, data):
"""Apply the forward
add - Whether or not to add to data
model - Model
data - Data """
self.checkDomainRange(model, data)
if not add:
data.zero()
forwardMult(self._matrix, model.arr, data.arr)
def adjoint(self, add, model, data):
"""Apply the adjoint
add - Whether or not to add to data
model - Model
data - Data """
self.checkDomainRange(model, data)
if not add:
model.zero()
adjointMult(self._matrix, model.arr, data.arr)
@jit(nopython=True)
def forwardMult(mat, model, data):
for j in range(mat.shape[1]):
for i in range(mat.shape[0]):
data[j] += mat[j][i] * model[i]
@jit(nopython=True)
def adjointMult(mat, model, data):
for j in range(mat.shape[1]):
for i in range(mat.shape[0]):
model[i] += mat[j][i] * data[j]
```
Iterative solver.ipynb 0 → 100644
+ 502
- 0
View file @ 4d84dd04
This source diff could not be displayed because it is too large. You can view the blob instead.
Untitled2.ipynb 0 → 100644
+ 2656
- 0
View file @ 4d84dd04
This diff is collapsed.
Untitled3.ipynb 0 → 100644
+ 119
- 0
View file @ 4d84dd04
%% Cell type:code id: tags:
``` python
from numba import jit
import Triangle
def wavekill(aa,bb,data):
s11=data.cone()
s12=data.clone()
s21=data.clone()
s22=data.clone()
out=data.clone()
s11.set(-aa)
s12.set(aa)
s21.set(-aa)
s22.set(aa)
s11.scaleAdd(bb,1,-1.)
s12.scaleAdd(bb,1.,-1.)
s21.scaleAdd(bb,1.,1.)
s22.scaleAdd(bb,1.,1.)
out.zero()
waveKillGit(uu.getNdArray(),s11.getNdArray(),s12.getNdArray(),s21.getNdArray(),s22.getNdArray(),out.getNdArray())
@jit(nopython=True)
def waveKillGit(uu,s11,s12,s21,s22,out):
for i2 in range(uu.shape[0]-1):
for i1 in range(uu.shape[1]-1):
out[i2,i1]=uu[i2,i1]*s11[i2,i1]+\
uu[i2+1,i1]*s12[i2,i1]+\
uu[i2,i1+1]*s21[i2,i1]+\
uu[i2+1,i1+1]*s22[i2,i1]
for i1 in range(uu.shape[1]):
out[uu.shape[0]-1,i1]=out[uu.shape[0]-2,i1]
for i2 in range(uu.shape[0]):
out[i2,uu.shape[1]-1]=out[i2,uu.shape[1]-1]
def pick2d(dat,boxsz):
pp=dat.clone()
dx=dat.clone()
dt=dat.clone()
pp.set(0)
wavekill(1.,pp,dat)
pp.set(1)
wavekill(0.,pp,dat)
dtdx=dt.clone()
dxdx=dx.clone()
dtdt=dt.clone()
dtdt.mult(dt) #dt*dt
dxdx.mult(dx) #dx*dx
dtdx.mult(dx)
triOp=Triangle.triangle2D_1(dat,dtdx,boxsz)
#dtdt
pp=dt.clone()
pp.mult(dt)
triOp.forward(False,pp,dtdt)
#dtdx
pp=dt.clone()
pp.mult(dx)
triOp.forward(False,pp,dtdx)
pp=dx.clone()
pp,mult(dx)
triOp.forward(False,pp,dx,dx)
coh=dat.clone()
```
%% Cell type:code id: tags:
``` python
```
Untitled4.ipynb 0 → 100644
+ 51
- 0
View file @ 4d84dd04
%% Cell type:code id: tags:
``` python
```
Untitled5.ipynb 0 → 100644
+ 65
- 0
View file @ 4d84dd04
%% Cell type:code id: tags:
``` python
class matmultOp
```
Vesuvio.ipynb 0 → 100644
+ 3669
- 0
View file @ 4d84dd04
This diff is collapsed.
chapter1/.ipynb_checkpoints/00-Intro-checkpoint.ipynb 0 → 100644
+ 445
- 0
View file @ 4d84dd04
%% Cell type:markdown id: tags:
# Basic operators and adjoints
A great many of the calculations
we do in science and engineering
are really matrix multiplication in disguise.
The first goal of this chapter is to unmask the disguise
by showing many examples.
Second, we see how the
**adjoint** operator (matrix transpose)
back projects information from data to the underlying model.
Geophysical modeling calculations
generally use linear operators that predict data from models.
Our usual task is to find the inverse of these calculations;
i.e., to find models (or make images) from the data.
Logically, the adjoint is the first step
and a part of all subsequent steps in this \bx{inversion} process.
Surprisingly, in practice, the adjoint sometimes does a better job
than the inverse!
Better because the adjoint operator tolerates imperfections in the data
and does not demand the data provide full information.
If you will permit me a poet's license with words,
I will offer you the following table
of **operator**s and their **adjoint**s
<table>
<tr>
<th>matrix multiply </th><th>conjugate-transpose matrix multiply</th></tr>
<tr>
<th>convolve </th><th>crosscorrelate </th></tr>
<tr>
<th>truncate </th><th>zero pad </th></tr>
<tr>
<th>replicate, scatter, spray </th><th>sum or stack </th></tr>
<tr>
<th>spray into neighborhoods </th><th>sum within bins </th></tr>
<tr>
<th>derivative (slope) </th><th>negative derivative </th></tr>
<tr>
<th>causal integration </th><th>anticausal integration |
<tr>
<th>add functions </th><th>do integrals </th></tr>
<tr>
<th>assignment statements </th><th>added terms </th></tr>
<tr>
<th>plane-wave superposition </th><th>slant stack / beam form |
<tr>
<th>spread on a curve </th><th>sum along a curve </th></tr>
<tr>
<th>stretch </th><th>squeeze </th></tr>
<tr>
<th>scalar field gradient </th><th>negative of vector field divergence</th></tr>
<tr>
<th>upward continue </th><th>downward continue</th></tr>
<tr>
<th>diffraction modeling </th><th>imaging by migration</th></tr>
<tr>
<th>hyperbola modeling </th><th>stacking for image or velocity </th></tr>
<tr>
<th>chop image into overlapping patches </th><th>merge the patches</th></tr>
<tr>
<th>ray tracing </th><th>tomography </th></tr>
</table>
%% Cell type:markdown id: tags:
The left column is often called **modeling**,
and the adjoint operators in the right column are often
used in data **processing**.
When the adjoint operator is
*not*
an adequate approximation to the inverse,
then you apply the techniques of fitting and optimization
explained in Chapter 2.
These techniques require iterative use of the
modeling operator and its adjoint.
The adjoint operator is sometimes called
the ``**back projection**'' operator,
because information propagated in one direction (Earth to data) is projected
backward (data to Earth model).
Using complex-valued operators,
the transpose and complex conjugate go together;
and in **Fourier analysis**, taking the complex conjugate
of $\exp(i\omega t)$ reverses the sense of time.
With more poetic license, I say that adjoint operators
*undo*
the time and therefore,
phase shifts of modeling operators.
The inverse operator does also,
but it also divides out the color.
For example, when linear interpolation is done,
then high frequencies are smoothed out; inverse interpolation must restore the
m.
You can imagine the possibilities for noise amplification
which is why adjoints are safer than inverses.
But, nature determines in each application what is the best operator to use
and whether to stop after the adjoint,
to go the whole way to the inverse,
or to stop partway.
%% Cell type:markdown id: tags:
The operators and adjoints previously shown transform vectors to other vectors
.
They also transform data planes to model planes, volumes, etc.
A mathematical operator transforms an "abstract vector" that
might be packed full of volumes of information like television
signals (time series) can pack together a movie, a sequence of frames.
We can always think of the operator as being a matrix,
but the matrix can be truly huge (and nearly empty).
When the vectors transformed by the matrices are large like
geophysical data set sizes,
then the matrix sizes are "large squared,"
far too big for computers.
Thus,
although we can always think of an operator as a matrix;
in practice, we handle an operator differently.
Each practical application requires the practitioner to
prepare two computer programs.
One performs the matrix multiply
$\bf y =\bf B \bf x$,
while the other multiplies by the transpose
$\tilde{\bf x} =\bf B^T \bf y$
(without ever having the matrix itself in memory).
It is always easy to transpose a matrix.
It is less easy to take a computer program that does
$\bf y =\bf B \bf x$
and convert it to another to do
$\tilde{\bf x} =\bf B^T \bf y$,
which is what we'll be doing here.
In this chapter are many examples of increasing complexity.
At the end of the chapter,
we see a test for any program pair
to see whether the operators $\bf B$ and $\bf B^T$ are mutually adjoint
as they should be.
Doing the job correctly (coding adjoints without making approximations)
rewards us later when we tackle model and image-estimation applications.
Mathematicians often denote the transpose of a matrix
$\bf B$ by $\bf B^{\rm T}$.
In physics and engineering, we often encounter complex numbers.
There, the adjoint is the complex-conjugate transposed matrix
denoted $\bf B^\ast$.
%% Cell type:markdown id: tags:
# Object orriented design
We need to define some abstract types before we proceed. The first is a vector. In addition to a list of values, a vector lives in some vector space. The space tells what these list of numbers correspond to. For example: how many samples, regular or irregular, how many dimensions, sampling of the dimensions, etc? Our abstract class needs to be able tell if another vector is from the same space. In addtion it need to make a copy of it space. We will add one more function that is able to zero itself.
%% Cell type:code id: tags:
``` python
# Vector class and derived classes
class vector:
"""Abstract python vector class"""
def __init__(self):
"""Default constructor"""
pass;
def zero(self):
"""Function to zero out a vector"""
raise NotImplementedError("zero must be overwritten")
def clone(self):
"""Function to clone (deep copy) a vector from a vector or a Space"""
raise NotImplementedError("clone must be overwritten")
def cloneSpace(self):
"""Function to clone vector space"""
raise NotImplementedError("cloneSpace must be overwritten")
def checkSame(self):
"""Function to check to make sure the vectors exist in the same space"""
raise NotImplementedError("checkSame must be overwritten")
```
%% Cell type:markdown id: tags:
An operator can be thought of as a matrix that transforms from its domain to its range using its forward and from range to its domain using its adjoint. An operator must be able to store its domain and range and check to make sure the model and data passed into it (when doing the forward and adjoint) match the domain and range it was initialized with,
%% Cell type:code id: tags:
``` python
class Operator:
"""Abstract python operator class"""
# Default class methods/functions
def __init__(self, domain, range):
"""Generic class for operator"""
self.domain = domain.cloneSpace()
self.range = range.cloneSpace()
def forward(self,add,model,data):
"""Apply forward of operator"""
raise NotImplementedError("forward must be overwritten")
def adjoint(self,add,model,data):
"""Apply adjoint of operator"""
raise NotImplementedError("adjoint must be overwritten")
```
%% Cell type:markdown id: tags:
We will make a concrete example of our vector class using numpy ndArrays.
%% Cell type:code id: tags:
``` python
import numpy
class numpyVector(vector):
def __init__(self,**kw):
if "dims" in kw:
self.v=numpy.ndArray(kw[dims])
self.shape=self.v.shape
self.spaceOnly=False
elif "array" in kw:
self.v=kw["array"]
self.shape=self.v.shape
self.spaceOnly=False
elif "space" in kw:
self.shape=kw["space"]
self.spaceOnly=True;
else:
raise ValueError("Only know how to create a numpyVector with dims=")
def zero(self):
"""Zero arrray"""
self.v.fill(0.)
def cloneSpace(self):
"""Return a numpy vector is just shape"""
return numpyVector(space=self.shape)
def clone(self):
"""Return copy of vector"""
return numpyVector(array=self.v.copy())
def checkSame(self,vec):
"""CHeck to make sure if vector in same space"""
if isinstance(vec,numpyVector):
if len(vec.shape)== len(self.shape):
for i in range(len(self.shape)):
if self.shape[i] != vec.shape[i]:
raise ValueError("Axis %d - %d is not same as %d "%(i,self.shape[i],vec.shape[i]))
else:
raise ValueError("Dimensions not the same %d %d"%(len(self.shape),len(vec.shape)))
else:
raise ValueError("Not a numpyVector")
return True
```
%% Cell type:markdown id: tags:
## Programming linear operators
The operation $ y_i = \sum_j b_{ij} x_j$ is the multiplication
of a matrix $\bf B$ by a vector $\bf x$.
The adjoint operation is
$\tilde x_j = \sum_i b_{ij} y_i$.
The operation adjoint to multiplication by a matrix
is multiplication by the transposed matrix
(unless the matrix has complex elements,
in which case,
we need the complex-conjugated transpose).
The following code does matrix multiplcation $\bf y = \bf B \bf x$ and
multiplication by the transpose $\tilde{ \bf x} = \bf B^T \bf y$.
We will define a python class that is initialized by a matrix. We will define
a forward and adjoint operator for the class.
%% Cell type:code id: tags:
``` python
class matmult(pyOperator.Operator):
"""Operator that does matrix multiplication
"""
def __init__(self, matrix, model, data):
"""Initialize operator
matrix - 2-D array
model - Model space
data - Data space
Both must be the same space, derived from pyVector.vectorIC"""
if not isinstance(model, giee.vector):
raise Exception(
"Model must begiee.vector or derived from it")
if not isinstance(data, giee.vector):
raise Exception(
"Model must be giee.vector or derived from it")
if not isinstance(matrix, np.ndarray):
raise Excecption(
"Matrix must be ndarray")
if model.ndims != 1:
raise Exception("Expecting model to be a 1-D array")
if data.ndims != 1:
raise Exception("Expecting data to be a 1-D array")
self._matrix = matrix
if self._matrix.ndim != 2:
raise Exception("Expecting matrix to be 2-D")
if self._matrix.shape[0] != model.getHyper().axes[0].n:
raise Exception(
"Expecting matrix first axis equal to length of model")
if self._matrix.shape[1] != data.getHyper().axes[0].n:
raise Exception(
"Expecting matrix second axis equal to length of data")
self.setDomainRange(model, data)
def forward(self, add, model, data):
"""Apply the forward
add - Whether or not to add to data
model - Model
data - Data """
self.checkDomainRange(model, data)
if not add:
data.zero()
forwardMult(self._matrix, model.arr, data.arr)
def adjoint(self, add, model, data):
"""Apply the adjoint
add - Whether or not to add to data
model - Model
data - Data """
self.checkDomainRange(model, data)
if not add:
model.zero()
adjointMult(self._matrix, model.arr, data.arr)
@jit(nopython=True)
def forwardMult(mat, model, data):
for j in range(mat.shape[1]):
for i in range(mat.shape[0]):
data[j] += mat[j][i] * model[i]
@jit(nopython=True)
def adjointMult(mat, model, data):
for j in range(mat.shape[1]):
for i in range(mat.shape[0]):
model[i] += mat[j][i] * data[j]
```
chapter1/.ipynb_checkpoints/Convolution-checkpoint.ipynb 0 → 100644
+ 2917
- 0
View file @ 4d84dd04
This diff is collapsed.
chapter1/.ipynb_checkpoints/Derivative on topography-checkpoint.ipynb 0 → 100644
+ 1904
- 0
View file @ 4d84dd04
This diff is collapsed.
chapter1/.ipynb_checkpoints/NMO-Stack-checkpoint.ipynb 0 → 100644
+ 1846
- 0
View file @ 4d84dd04
This diff is collapsed.
chapter1/.ipynb_checkpoints/Triangle edge-checkpoint.ipynb 0 → 100644
+ 1793
- 0
View file @ 4d84dd04
This diff is collapsed.
chapter1/.ipynb_checkpoints/galilee-binning-checkpoint.ipynb 0 → 100644
+ 1836
- 0
View file @ 4d84dd04
This diff is collapsed.
chapter1/00-Intro.ipynb 0 → 100644
+ 445
- 0
View file @ 4d84dd04
%% Cell type:markdown id: tags:
# Basic operators and adjoints
A great many of the calculations
we do in science and engineering
are really matrix multiplication in disguise.
The first goal of this chapter is to unmask the disguise
by showing many examples.
Second, we see how the
**adjoint** operator (matrix transpose)
back projects information from data to the underlying model.
Geophysical modeling calculations
generally use linear operators that predict data from models.
Our usual task is to find the inverse of these calculations;
i.e., to find models (or make images) from the data.
Logically, the adjoint is the first step
and a part of all subsequent steps in this \bx{inversion} process.
Surprisingly, in practice, the adjoint sometimes does a better job
than the inverse!
Better because the adjoint operator tolerates imperfections in the data
and does not demand the data provide full information.
If you will permit me a poet's license with words,
I will offer you the following table
of **operator**s and their **adjoint**s
<table>
<tr>
<th>matrix multiply </th><th>conjugate-transpose matrix multiply</th></tr>
<tr>
<th>convolve </th><th>crosscorrelate </th></tr>
<tr>
<th>truncate </th><th>zero pad </th></tr>
<tr>
<th>replicate, scatter, spray </th><th>sum or stack </th></tr>
<tr>
<th>spray into neighborhoods </th><th>sum within bins </th></tr>
<tr>
<th>derivative (slope) </th><th>negative derivative </th></tr>
<tr>
<th>causal integration </th><th>anticausal integration |
<tr>
<th>add functions </th><th>do integrals </th></tr>
<tr>
<th>assignment statements </th><th>added terms </th></tr>
<tr>
<th>plane-wave superposition </th><th>slant stack / beam form |
<tr>
<th>spread on a curve </th><th>sum along a curve </th></tr>
<tr>
<th>stretch </th><th>squeeze </th></tr>
<tr>
<th>scalar field gradient </th><th>negative of vector field divergence</th></tr>
<tr>
<th>upward continue </th><th>downward continue</th></tr>
<tr>
<th>diffraction modeling </th><th>imaging by migration</th></tr>
<tr>
<th>hyperbola modeling </th><th>stacking for image or velocity </th></tr>
<tr>
<th>chop image into overlapping patches </th><th>merge the patches</th></tr>
<tr>
<th>ray tracing </th><th>tomography </th></tr>
</table>
%% Cell type:markdown id: tags:
The left column is often called **modeling**,
and the adjoint operators in the right column are often
used in data **processing**.
When the adjoint operator is
*not*
an adequate approximation to the inverse,
then you apply the techniques of fitting and optimization
explained in Chapter 2.
These techniques require iterative use of the
modeling operator and its adjoint.
The adjoint operator is sometimes called
the ``**back projection**'' operator,
because information propagated in one direction (Earth to data) is projected
backward (data to Earth model).
Using complex-valued operators,
the transpose and complex conjugate go together;
and in **Fourier analysis**, taking the complex conjugate
of $\exp(i\omega t)$ reverses the sense of time.
With more poetic license, I say that adjoint operators
*undo*
the time and therefore,
phase shifts of modeling operators.
The inverse operator does also,
but it also divides out the color.
For example, when linear interpolation is done,
then high frequencies are smoothed out; inverse interpolation must restore the
m.
You can imagine the possibilities for noise amplification
which is why adjoints are safer than inverses.
But, nature determines in each application what is the best operator to use
and whether to stop after the adjoint,
to go the whole way to the inverse,
or to stop partway.
%% Cell type:markdown id: tags:
The operators and adjoints previously shown transform vectors to other vectors
.
They also transform data planes to model planes, volumes, etc.
A mathematical operator transforms an "abstract vector" that
might be packed full of volumes of information like television
signals (time series) can pack together a movie, a sequence of frames.
We can always think of the operator as being a matrix,
but the matrix can be truly huge (and nearly empty).
When the vectors transformed by the matrices are large like
geophysical data set sizes,
then the matrix sizes are "large squared,"
far too big for computers.
Thus,
although we can always think of an operator as a matrix;
in practice, we handle an operator differently.
Each practical application requires the practitioner to
prepare two computer programs.
One performs the matrix multiply
$\bf y =\bf B \bf x$,
while the other multiplies by the transpose
$\tilde{\bf x} =\bf B^T \bf y$
(without ever having the matrix itself in memory).
It is always easy to transpose a matrix.
It is less easy to take a computer program that does
$\bf y =\bf B \bf x$
and convert it to another to do
$\tilde{\bf x} =\bf B^T \bf y$,
which is what we'll be doing here.
In this chapter are many examples of increasing complexity.
At the end of the chapter,
we see a test for any program pair
to see whether the operators $\bf B$ and $\bf B^T$ are mutually adjoint
as they should be.
Doing the job correctly (coding adjoints without making approximations)
rewards us later when we tackle model and image-estimation applications.
Mathematicians often denote the transpose of a matrix
$\bf B$ by $\bf B^{\rm T}$.
In physics and engineering, we often encounter complex numbers.
There, the adjoint is the complex-conjugate transposed matrix
denoted $\bf B^\ast$.
%% Cell type:markdown id: tags:
# Object orriented design
We need to define some abstract types before we proceed. The first is a vector. In addition to a list of values, a vector lives in some vector space. The space tells what these list of numbers correspond to. For example: how many samples, regular or irregular, how many dimensions, sampling of the dimensions, etc? Our abstract class needs to be able tell if another vector is from the same space. In addtion it need to make a copy of it space. We will add one more function that is able to zero itself.
%% Cell type:code id: tags:
``` python
# Vector class and derived classes
class vector:
"""Abstract python vector class"""
def __init__(self):
"""Default constructor"""
pass;
def zero(self):
"""Function to zero out a vector"""
raise NotImplementedError("zero must be overwritten")
def clone(self):
"""Function to clone (deep copy) a vector from a vector or a Space"""
raise NotImplementedError("clone must be overwritten")
def cloneSpace(self):
"""Function to clone vector space"""
raise NotImplementedError("cloneSpace must be overwritten")
def checkSame(self):
"""Function to check to make sure the vectors exist in the same space"""
raise NotImplementedError("checkSame must be overwritten")
```
%% Cell type:markdown id: tags:
An operator can be thought of as a matrix that transforms from its domain to its range using its forward and from range to its domain using its adjoint. An operator must be able to store its domain and range and check to make sure the model and data passed into it (when doing the forward and adjoint) match the domain and range it was initialized with,
%% Cell type:code id: tags:
``` python
class Operator:
"""Abstract python operator class"""
# Default class methods/functions
def __init__(self, domain, range):
"""Generic class for operator"""
self.domain = domain.cloneSpace()
self.range = range.cloneSpace()
def forward(self,add,model,data):
"""Apply forward of operator"""
raise NotImplementedError("forward must be overwritten")
def adjoint(self,add,model,data):
"""Apply adjoint of operator"""
raise NotImplementedError("adjoint must be overwritten")
```
%% Cell type:markdown id: tags:
We will make a concrete example of our vector class using numpy ndArrays.
%% Cell type:code id: tags:
``` python
import numpy
class numpyVector(vector):
def __init__(self,**kw):
if "dims" in kw:
self.v=numpy.ndArray(kw[dims])
self.shape=self.v.shape
self.spaceOnly=False
elif "array" in kw:
self.v=kw["array"]
self.shape=self.v.shape
self.spaceOnly=False
elif "space" in kw:
self.shape=kw["space"]
self.spaceOnly=True;
else:
raise ValueError("Only know how to create a numpyVector with dims=")
def zero(self):
"""Zero arrray"""
self.v.fill(0.)
def cloneSpace(self):
"""Return a numpy vector is just shape"""
return numpyVector(space=self.shape)
def clone(self):
"""Return copy of vector"""
return numpyVector(array=self.v.copy())
def checkSame(self,vec):
"""CHeck to make sure if vector in same space"""
if isinstance(vec,numpyVector):
if len(vec.shape)== len(self.shape):
for i in range(len(self.shape)):
if self.shape[i] != vec.shape[i]:
raise ValueError("Axis %d - %d is not same as %d "%(i,self.shape[i],vec.shape[i]))
else:
raise ValueError("Dimensions not the same %d %d"%(len(self.shape),len(vec.shape)))
else:
raise ValueError("Not a numpyVector")
return True
```
%% Cell type:markdown id: tags:
## Programming linear operators
The operation $ y_i = \sum_j b_{ij} x_j$ is the multiplication
of a matrix $\bf B$ by a vector $\bf x$.
The adjoint operation is
$\tilde x_j = \sum_i b_{ij} y_i$.
The operation adjoint to multiplication by a matrix
is multiplication by the transposed matrix
(unless the matrix has complex elements,
in which case,
we need the complex-conjugated transpose).
The following code does matrix multiplcation $\bf y = \bf B \bf x$ and
multiplication by the transpose $\tilde{ \bf x} = \bf B^T \bf y$.
We will define a python class that is initialized by a matrix. We will define
a forward and adjoint operator for the class.
%% Cell type:code id: tags:
``` python
class matmult(pyOperator.Operator):
"""Operator that does matrix multiplication
"""
def __init__(self, matrix, model, data):
"""Initialize operator
matrix - 2-D array
model - Model space
data - Data space
Both must be the same space, derived from pyVector.vectorIC"""
if not isinstance(model, giee.vector):
raise Exception(
"Model must begiee.vector or derived from it")
if not isinstance(data, giee.vector):
raise Exception(
"Model must be giee.vector or derived from it")
if not isinstance(matrix, np.ndarray):
raise Excecption(
"Matrix must be ndarray")
if model.ndims != 1:
raise Exception("Expecting model to be a 1-D array")
if data.ndims != 1:
raise Exception("Expecting data to be a 1-D array")
self._matrix = matrix
if self._matrix.ndim != 2:
raise Exception("Expecting matrix to be 2-D")
if self._matrix.shape[0] != model.getHyper().axes[0].n:
raise Exception(
"Expecting matrix first axis equal to length of model")
if self._matrix.shape[1] != data.getHyper().axes[0].n:
raise Exception(
"Expecting matrix second axis equal to length of data")
self.setDomainRange(model, data)
def forward(self, add, model, data):
"""Apply the forward
add - Whether or not to add to data
model - Model
data - Data """
self.checkDomainRange(model, data)
if not add:
data.zero()
forwardMult(self._matrix, model.arr, data.arr)
def adjoint(self, add, model, data):
"""Apply the adjoint
add - Whether or not to add to data
model - Model
data - Data """
self.checkDomainRange(model, data)
if not add:
model.zero()
adjointMult(self._matrix, model.arr, data.arr)
@jit(nopython=True)
def forwardMult(mat, model, data):
for j in range(mat.shape[1]):
for i in range(mat.shape[0]):
data[j] += mat[j][i] * model[i]
@jit(nopython=True)
def adjointMult(mat, model, data):
for j in range(mat.shape[1]):
for i in range(mat.shape[0]):
model[i] += mat[j][i] * data[j]
```
chapter1/Convolution.ipynb 0 → 100644
+ 2917
- 0
View file @ 4d84dd04
This diff is collapsed.
chapter1/Derivative on topography.ipynb 0 → 100644
+ 1965
- 0
View file @ 4d84dd04
This diff is collapsed.
chapter1/NMO-Stack.ipynb 0 → 100644
+ 1846
- 0
View file @ 4d84dd04
This diff is collapsed.
chapter1/Triangle edge.ipynb 0 → 100644
+ 1793
- 0
View file @ 4d84dd04
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment