F = JopAtan(spc, c)

where F is the shifted arc tangent operator F = arctan( x(t)/c ) + π/2

x(t) = 0 corresponds to F = π/2 if |x(t)| ≤ c, then π/4 ≤ F ≤ 3π/4

A = JopBlend(T, nsamples, shottimes[, nsamples_blended])

where A is a shot mixing operation, so that if d=Am, then m::Array{T,2} is a common receiver gather with each trace corresponding to a different shot. The excitation time of each are the entries of shottimes::Array{Int,1} in units of time samples. The size of m is (nsamples,length(shottimes)). The size of d is either nsamples_blended (if provided) or maximum(shot_times)+nsamples-1.

Examples

blend two receiver gather traces:

A = JopBlend(Float64, 128, [1,64])
m = rand(domain(A))
d = A*m # receiver gather with blended shots

A = JopCircShift(R::JetSpace, shifts)

d=A*m is equivalent to circshift(d, m, shifts), and where m and d are in the domain/range of A.

A = JopDerivative(R::JetSpace[; dim=1, accuracy=4, delta=1.0)

A*m is the centered finite different approximation to the derivative of m along dimension dim. The accuracy of the approximation is either 4 for a 4th order accurate estimate of the derivative, or 8 for an 8th order accurate derivative.

A = JopDiagonal(spc, d) A = JopDiagonal(d)

where spc::JetSpace is the domain/range of A, and d::Array or d::Number is the diagonal. If d<:Number, then the diagonal of the matrix is constant and one must specify spc.

Examples

A = JopDiagonal(spc, d)

where spc::JetSpace is the domain and range of A, and d<:Array or d<:Number

A = JotOpDiagonal(d)

where d<:AbstractArray. The domain and range of A are determined by the size and type of d.

A = JotOpDiagonal([1.0, 2.0, 3.0])
m = ones(domain(A))
d = A*m # d = [1.0 ; 2.0 ; 3.0]

A = JopDifference(R::JetSpace[, dim=1])

A*m is similar to diff(a,dims=dim) where a is a vector in the space R. In other words, it is the one-sided difference of a along the dimension dim.

F = JopErf(spc)

where F is the error function operator with domain and range given by spc::JetSpace. we use 'erf'(z) = 2/√π ∫_{0}^{z} xp{-t^2}dt, the derivative of which is 2/√π xp{-z^2}. We expect the domain and range to be real.

F = JopExp(spc, c)

where F is the exponential operator e^(c/x) with domain and range given by spc::JetSpace, and scalar value c.

A = JopGradient(dom, δ)

Gradient for a 2D or 3D arrays. The range will have one more dimension than the domain and contain the components of the gradient in each dimension.

• dom::JetSpace{T,N} is the domain of the operator.
• δ::NTuple{T,N} is the grid spacing in each dimension.

Examples:

2D

nz,nx = 11,12
dom = JetSpace(Float32, nz, nx)
A = JopGradient(dom, (1.0,1.0))
size(domain(A)) # (11,12)
size(range(A))  # (11,12,2)

3D

nz,ny,nx = 11,12,13
dom = JetSpace(Float32, nz, ny, nx)
A = JopGradient(dom,(1.0,1.0,1.0))
size(domain(A)) # (11,12,13)
size(range(A))  # (11,12,13,3)

A = JopHighpass(sp)

Build a 2D Highpass operator, which first apply a center weighted three point smoothing operator [0.25, 0.5, 0.25] to each dimension with a given iterations number (nit[1:2]) to get the low frequency component of the data, then subtract it from the original data to return the high pass component. It works with sp::JetSpace, nit::Array{Integer, 1} and A::JopHighpass. It can specify different number of iterations (nit) in x and z direction. For a seismic image, more iterations in z direction is often applied to obtain more vertical resolution.

For example,

A = JopHighpass(JetSpace(Float64,128,256), nit=(3,1)])
m = rand(domain(A))
d = A*m

op = JopImag(dom)

Extract the imaginary part of a complex input array where dom::JetSpace{<:Complex}.

JopInterp(dom, rng)

Performs linear interpolation from dom::JetSpace to rng::JetSpace. It is often used to reduce dimensionality in FWI. JopInterp is ported from the CVX frequency domain FWI tools in SeisSpace.

Notes

• It is required that domain and range have the same dimensionality.
• It is required that domain and range both have more than 2 points per dimension
• It is required that the domain is coarser than the range.
• It is assumed that domain and range have the same boundary (e.g. in 2D: xmin,xmax,zmin,zmax).

Examples:

2D

dom = JetSpace(Float32, 11, 11)
rng = JetSpace(Float32, 23, 23)
A = JopInterp(dom, rng)
y = A * rand(dom)
x = A' * rand(rng)

A = JopLMO(spc [;v=2.0, dx=1.0, dt=1.0])

2D linear moveout operator for domain and range sp::JetSpace. The forward operator maps from flat events to dipping events using the parameters v,dx,dt, and:

(z_moveout) = (z + (dx/dt)/v * x)

Notes

• It should be trivial to add a 3D implementation for this operator.

A = JopLaplacian(sp)

Build a 2D (d^2/dx^2+d^2/dy^2) or 3D (d^2/dx^2+d^2/dy^2+d^2/dz^2) Laplacian operator with three point centered finite difference stencil, with sp::JetSpace and A::JopLaplacian.

For example,

A = JopLaplacian(JetSpace(Float64,128,256))
m = rand(domain(A))
d = A*m

F = JopLog(spc)

where F is the log operator with domain and range given by spc::JetSpace.

JopMix(s::JetAbstractSpace{T,N}, nmix::NTuple{N,Int})

2D spatial mix, and is used to add smoothness to a domain. JopMix is ported from the CVX frequency domain FWI tools in SeisSpace.

Notes:

• This appears to over-lap in funtionality with JopRoughness.

Examples:

2D

nz, nx = 11,11
spc = JetSpace(Float32, nz, nx)
A = JopMix(spc, (5,5))
y = A * rand(dom)
x = A' * rand(rng)

F = JopNim(spc, c)

where F is the 'normalized integral method' operator F = int0^t [ x(t) dt ] / int0_T [ x(t) dt ]

F = JopNormalize(spc, c)

where F is the 'normalize by maximuam value of integral' operator F(t) = x(t) / int_0^T [ x(t) dt ] Note the integration is along the 1st (fastest) dimension.

A = JopPad(dom, pad..., [extend=false, accumulate=false])

where dom::JetSpace is the domain of A, and pad::UnitRange... determines the range of A. If extend=false, then the padded region is set to zero. If extend=true, then the padded region is set from the boundary of the domain. The accumulate=true option is specific to the Ginsu operation in JetPackWave, and should not be used unless you really know what you are doing.

Examples:

1D

A = JopPad(JetSpace(Float64,2), -1:3)
m = [1.0, 2.0]
d = A*m # d = [0.0, 0.0, 1.0, 2.0, 0.0]
A = JopPad(JetSpace(Float64,2), -1:3, extend=true)
d = A*m # d = [1.0, 1.0, 1.0, 2.0, 2.0]

2D

A = JopPad(JetSpace(Float64,2,2), -1:3, 1:3)
m = [11. 12. ; 21. 22.]
d = A*m # d = [0. 0. 11. 12. 0. ; 0. 0. 21. 22. 0. ; 0. 0. 0. 0. 0.]

Notes:

• This operator may also be used for truncation
• There may be overlap between the functionality of JopPad and JopRestriction, and it may be worth

thinking of how to consolidate them.

A = JopPermute(sp, dims, perm)

where sp::JetSpace is the domain and range of the operator, dims::NTuple is a tuple of dimensions that one wishes to permute, and perm is a tuple of arrays with permutation indices.

Example 1:

A = JopPermute(JetSpace(Float64,4,2),(1,),([3;1;2;4],))
m = [11 12 ; 21 22 ; 31 32 ; 41 42]
d = A*m # d=[31 32 ; 11 12 ; 21 22 ; 41 42]

Example 2:

A = JopPermute(JetSpace(Float64,3,2),(1,2),([3;2;1],[2;1]))
m = [11 12 ; 21 22 ; 31 32]
d = A*m # d = [32 31 ; 22 21 ; 12 11]

Notes:

Currently, this is only implmented for 1D, 2D and 3D arrays.

A = JopPermutedims(sp, perm)

where sp::JetSpace is the domain of the operator, and perm is an array specifying the permutation of array dimensions. The range of the operator is inferred from sp and perm.

Example:

A = JopPermutedims(JetSpace(Float64,3,2),[2;1])
m = [1 2 3 ; 4 5 6]
d = A*m # d = [1 4 ; 2 5 ; 3 6]

F = JopPow(spc, c, a)

where F is the power operator (x/c)^a with domain and range given by spc::JetSpace, and scalar values c, a.

A = JopProjection(u::AbstractArray)

A*m is the projection of vector m onto the vector u.

op = JopReal(dom)

Extract the real part of a complex input array where dom::JetSpace{<:Complex}.

A = JopReghost(spc, zt, zm, dt, dx)

Build a 2D receiver-side reghosting and redatuming operator in the FK domain using simple phase-shift operations (e.g. Posthumus, 1993). spc::JetSpaceSymmetric is the domain and range of A. zt is target depth where the receivers will be re-dataumed to, and zm is the measurement depth where the recording is made. dt is the time sampling interval and dx is the receiver sampling interval.

Since this operator is built in the FK domain, spc has symmetry in the frequency dimension. Hence spc::JotSpaceSymmetric. We provide a convenience method:

sp = symspace(JopReghost, T, nw, nk)

where T is either Float32 or Float64, nw are the number of frequency samples, and nk are the number of wavenumber samples.

Example:

P = JopPad(JetSpace(Float64,nt,nx), 1:2nt, 1:nx)
F = JopFft(range(P))
G = JopReghost(range(F), 20.0, 20.0, .004, 10.0)
A = P'*F'*G*F*P
data_reghost = A * data_noghost

Note that in the above example,

range(F) == symspace(JopReghost_df!, Float64, size(range(P))...)

JopRemoveDC(s::JetAbstractSpace{T,N})

F(x) = x - zero frequency component of x along the fastest dimension Uses Fourier transform, as the DC component differs from the mean as a function of interval length and zero padding. ```

A = JopReshape(dom, rng)

Reshape an array that belongs in dom::JetAbstractSpace to one that belongs to rng::JetAbstractSpace

Example

using Jets, JetPack
dom = JetSpace(Float32,10,20)
rng = JetSpace(Float32,200)
A = JopReshape(JetSpace(Float32,10,20), JetSpace(Float32,200))
x = rand(domain(A))
y = A*x

A = JopRestriction(dom[, rng], indices)

Apply a restriction operator mapping from dom::JetSpace to rng::JetSpace, and using a one-dimensional array of indices where indices is one of:

• indices::Vector{Int,1} - for restriction of 1-D arrays
• indices::Vector{NTuple{N,Int}} - for restriction of N-D arrays (N>1)

Note that in the case where the domain is N-dimension, the range is always stored using a 1-dimension array.

example (1-D):

A = JopRestriction(JotSpace(Float64,4),[1,3])
m = [1.0;2.0;3.0;4.0]
d = A*m # d=[1.0;3.0]

example (2-D):

A = JopRestriction(JotSpace(Float64,2,2),[(1,2),(2,2)])
m = [1.0 2.0;3.0 4.0]
d = A*m # d=[2.0,4.0]

A = JopRoughness(sp, dim, w)

A can be used to regularize an optimization problem, applying a penalty to models that are non-smooth. This is similar to a finite difference operator, but here no care is taken to ensure that A computes a derivative. dim is the dimension that one wishes to smooth and w is the half-width of the smoother. The width of the window determines how strong the penalty is for being non-smooth.

For example, the form of A for w=1 is,

A=
[
(1/3-2/3) 1/3       0         0         ;
1/3       (1/3-3/3) 1/3       0         ;
0         1/3       (1/3-3/3) 1/3       ;
0         0         1/3       (1/3-2/3)
]

A = JopShift(spc, b) A = JopShift(b)

A(x) = x + b

where spc::JetSpace is the domain/range of A, and b::Array or b::Number is the affine translation or shift. If b<:Number, then shift is constant and one must specify spc.

Examples

A = JopShift(spc, b)

where spc::JetSpace is the domain and range of A, and b<:Array or b<:Number

A = JopShift(b)

where b<:AbstractArray. The domain and range of A are determined by the size and type of b.

A = JopShift([1.0, 2.0, 3.0])
m = ones(domain(A))
d = A*m # d = [2.0 ; 3.0 ; 4.0]

F = JopSigmoid(spc, c)

where F is the sigmoid operator 1/(1+e^(-c/x)) with domain and range given by spc::JetSpace, and scalar value c.

F = JopTanh(spc)

where F is the hyperbolic tangent operator with domain and range given by spc::JetSpace. we use 'tanh'(cx) = (xp{cx}-xp{-cx})/(xp{cx}+xp{-cx}), the derivative of which is 1-tanh(cx)^2. We expect the domain and range to be real.

A = JopTaper(spc, dims, frac[, frac_end; mode=:normal, taper=(:cosine,:cosine)])

The linear operator A tapers the edges of all or some subset of the dimensions of an array that belongs to spc::JetSpace. The dimensions that are tapered are given by dims, an Int tuple. A will taper the beginning and/or end edges of each specified array dimension. The size of the beginning and end tapers are determined as a fraction of the length of that array dimension, and this fraction is set using the tuple frac, and (optionally) frac_end.

If frac_end::NTuple{M,Float64} is not set, then the ith entry of frac::NTuple{M,Float64} determines the length of the center portion of the ith dimension that is not tapered. If frac_end::NTuple{M,NTuple{2,Float64}} is set, then the length of the begining taper of the ith array dimension is determiend by frac[i], and the length of the end taper is determined by frac_end[i].

The optional named argument mode can be used if the taper is applied to a dimension that corresponds to FFT ordering where the edges are assumed to be at the center and left ends of the dimension.

The optional named argument taper is a tuple specifying what type of taper to use at each end. Available tapers are :cosine and :heaviside.

Examples:

taper for a space containing 1D arrays

A = JopTaper(JetSpace(Float64,10), (1,), (.75,))
m = ones(domain(A))
d = A*m

taper for a space containing 1D arrays, and only taper at the end

A = JopTaper(JetSpace(Float64,10), (1,), (0.0,), (0.25,))
m = ones(domain(A))
d = A*m

taper for a space containing 2D arrays, where both dimensions are tapered

A = JopTaper(JetSpace(Float64,10,11), (1,2), (0.75,0.5))
m = ones(domain(A))
d = A*m

taper for a space containing 3D arrays, where two of the three dimensions are tapered

A = JopTaper(JetSpace(Float64,10,11,12), (1,3), (0.75,0.5))
m = ones(domain(A))
d = A*m

In the next example, array dimension 1 is tapered at the end, and array dimension 2 is tapered at the beginning and end.

A = JopTaper(JetSpace(Float64,10,11,12), (1,3), (0.0,0.25), (0.25,0.25))
m = ones(domain(A))
d = A*m

A = JopTranslation(p1::Matrix, p2::Matrix)

A*m translates m::Matrix using the vector field (p1,p2).