Merge pull request #190 from slimgroup/misc-fixes

Misc fixes

Project.toml

@@ -1,7 +1,7 @@
 name = "JUDI"
 uuid = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
 authors = ["Philipp Witte, Mathias Louboutin"]
-version = "3.3.4"
+version = "3.3.5"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

docs/src/data_structures.md

@@ -42,7 +42,7 @@ model.n
 JUDI's geometry structure contains the information of either the source **or** the receiver geometry. Construct an (in-core) geometry object for **either** a source or receiver set up:
 
 ```@docs
- Geometry(xloc, yloc, zloc; dt=[], t=[], nsrc=nothing)
+ Geometry
 ```
 
 From the optional arguments, you have to pass (at least) **two** of `dt`, `nt` and `t`. The third value is automatically determined and set from the two other values. a `Geometry` can be constructed in a number of different ways for in-core and out-of-core cases. Check our examples and the source for additional details while the documentation is being extended.
@@ -62,6 +62,14 @@ geometry.xloc[i]
 geometry.xloc[i][j]
 ```
 
+### Geometry utilities
+
+A few utilities to manipulates geometries are provided as well.
+
+```@docs
+super_shot_geometry
+reciprocal_geom
+```
 ## Options structure
 
 The options structure allows setting several modeling parameters.

docs/src/helper.md

@@ -53,6 +53,7 @@ Create a probability distribution with the shape of the source spectrum from whi
 
 ```@docs
 generate_distribution
+select_frequencies
 ```
 
 We can draw random samples from `dist` by passing it values between 0 and 1:
@@ -83,3 +84,32 @@ d_ic = get_data(d_ooc, inds)
 ```
 
 where `inds` is either a single index, a list of index or a range of index.
+
+## Restrict model to acquisition
+
+In practice, and in particular for marine data, the aperture of a single shot is much smaller than the full model size. We provide a function (automatically used when the option `limit_m` is set in [`Options`](@ref)) that limits the model to the acquisition area.
+
+```@docs
+limit_model_to_receiver_area
+```
+
+We also provide it's complement that removes receivers outside of the computational domain if the dataset contains locations that are not wanted
+
+```@docs
+remove_out_of_bounds_receivers
+```
+
+## Additional miscellanous utilities
+
+```@docs
+devito_model
+setup_grid
+pad_sizes
+pad_array
+remove_padding
+convertToCell
+process_input_data
+reshape
+transducer
+as_vec
+```

docs/src/preconditioners.md

@@ -3,7 +3,7 @@
 JUDI provides a selected number of preconditioners known to be beneficial to FWI and RTM. We welcome additional preconditionners from the community. Additionnaly, any JOLI operator can be used as a preconditiner in conbination with JUDI operator thanks to the fundamental interface between JUDI and JOLI.
 
 ```@contents
-Pages = ["precinditioners.md"]
+Pages = ["preconditioners.md"]
 ```
 
 ## Model domain preconditioners
@@ -91,7 +91,13 @@ While not purely a preconditioner, because this operator acts on the data and is
 
 ```@docs
 FrequencyFilter
-low_filter
+filter_data
 ```
 
+### Data time derivative/intergraton
 
+A `TimeDifferential{K}` is a linear operator that implements a time derivative (K>0) or time integration (K<0) of order `K` for any real `K` including fractional values.
+
+```@docs
+TimeDifferential
+```

examples/notebooks/07_preconditionners.ipynb

@@ -616,7 +616,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 1.9.0-beta4",
+   "display_name": "Julia 1.9.0",
    "language": "julia",
    "name": "julia-1.9"
   },

src/JUDI.jl

@@ -53,9 +53,13 @@ import PyCall.NpyArray
 import ChainRulesCore: rrule
 
 # Set python paths
+export devito, set_devito_config
 const pm = PyNULL()
 const ac = PyNULL()
 const pyut = PyNULL()
+const devito = PyNULL()
+
+set_devito_config(key::String, val::String) = set!(devito."configuration", key, val)
 
 # Create a lock for pycall FOR THREAD/TASK SAFETY
 # See discussion at
@@ -156,6 +160,7 @@ function __init__()
     copy!(pm, pyimport("models"))
     copy!(ac, pyimport("interface"))
     copy!(pyut, pyimport("utils"))
+    copy!(devito, pyimport("devito"))
     # Initialize lock at session start
     PYLOCK[] = ReentrantLock()
 

src/TimeModeling/Preconditioners/DataPreconditioners.jl

@@ -1,5 +1,5 @@
-export DataMute, FrequencyFilter, judiTimeDerivative, judiTimeIntegration
-export judiFilter, low_filter, judiDataMute, muteshot
+export DataMute, FrequencyFilter, judiTimeDerivative, judiTimeIntegration, TimeDifferential
+export judiFilter, filter_data, judiDataMute, muteshot
 
 ############################################ Data mute ###############################################
 """
@@ -176,31 +176,41 @@ adjoint(I::FrequencyFilter{T}) where T = I
 transpose(I::FrequencyFilter{T}) where T = I
 
 function filter!(dout::AbstractArray{T, N}, din::AbstractArray{T, N}, dt::T; fmin=T(0.01), fmax=T(100)) where {T, N}
-    responsetype = Bandpass(fmin, fmax; fs=1e3/dt)	
+    if fmin == 0
+        responsetype = Lowpass(fmax; fs=1e3/dt)
+    elseif isinf(fmax)
+        responsetype = Highpass(fmin; fs=1e3/dt)
+    else
+        responsetype = Bandpass(fmin, fmax; fs=1e3/dt)
+    end
     designmethod = Butterworth(5)
     tracefilt!(x, y) = filt!(x, digitalfilter(responsetype, designmethod), y)
     map(i-> tracefilt!(selectdim(dout, N, i), selectdim(din, N, i)), 1:size(dout, 2))
 end
 
-low_filter(Din::judiVector; fmin=0.01, fmax=100.0) = judiFilter(Din.geometry, fmin, fmax)*Din
-low_filter(Din::judiVector, ::Any; fmin=0.01, fmax=100.0) = low_filter(Din; fmin=fmin, fmax=fmax)
+filter_data(Din::judiVector; fmin=0, fmax=Inf) = judiFilter(Din.geometry, fmin, fmax)*Din
+filter_data(Din::judiVector, ::Any; fmin=0, fmax=Inf) = filter_data(Din; fmin=fmin, fmax=fmax)
 
 """
-    low_filter(Din, dt_in; fmin=0, fmax=25)
+    filter(Din, dt_in; fmin=0, fmax=25)
 
-Performs a causal band pass filtering [fmin, fmax] on the input data bases on its sampling rate `dt`.
+Performs a causal filtering [fmin, fmax] on the input data bases on its sampling rate `dt`. 
+Automatically perfroms a lowpass if fmin=0 (default)
 """
-function low_filter(Din::Matrix{T}, dt_in; fmin=0.01, fmax=100.0) where T
+function filter_data(Din::Matrix{T}, dt_in; fmin=0, fmax=Inf) where T
     out = similar(Din)
     filter!(out, Din, dt_in; fmin=T(fmin), fmax=T(fmax))
     return out
 end
 
+# Legacy `low_filter`
+low_filter(ar...) = filter_data(ar...)
+
 # Legacy top mute is deprecated since only working for marine data
 judiMarineTopmute2D(::Integer, geometry::Geometry; params=Array{Any}(undef, 3), flipmask=false) = throw(MethodError(judiMarineTopmute2D, "judiMarineTopmute2D is deprecated due to its limitations and inaccuracy, please use judiDataMute"))
 
 """
-    TimeDifferential{K}
+    struct TimeDifferential
         recGeom
 
 Differential operator of order `K` to be applied along the time dimension. Applies the ilter `w^k` where `k` is the order. For example,
@@ -215,7 +225,6 @@ Constructor
     judiTimeDerivative(recGeom, order)
     judiTimeDerivative(judiVector, order)
 """
-
 struct TimeDifferential{T, K} <: DataPreconditioner{T, T}
     m::Integer
     recGeom::Geometry

src/TimeModeling/Preconditioners/base.jl

@@ -29,4 +29,10 @@ getproperty(J::Preconditioner, s::Symbol) = _get_property(J, Val{s}())
 *(J::Preconditioner, ms::PhysicalParameter) = matvec(J, ms)
 *(J::Preconditioner, v::Vector{T}) where T = matvec(J, v)
 mul!(out::judiMultiSourceVector, J::Preconditioner, ms::judiMultiSourceVector) = copyto!(out, matvec(J, ms))
-mul!(out::PhysicalParameter, J::Preconditioner, ms::PhysicalParameter) = copyto!(out, matvec(J, ms))
+mul!(out::PhysicalParameter, J::Preconditioner, ms::PhysicalParameter) = copyto!(out, matvec(J, ms))
+
+# Unsupported OOC
+function *(J::DataPreconditioner, v::judiVector{T, SegyIO.SeisCon}) where T
+    @warn  "Data preconditionners only support in-core judiVector. Converting (might run out of memory)"
+    return J * get_data(v)
+end

src/TimeModeling/Types/GeometryStructure.jl

@@ -4,6 +4,7 @@
 #
 
 export Geometry, compareGeometry, GeometryIC, GeometryOOC, get_nsrc, n_samples, super_shot_geometry
+export reciprocal_geom
 
 abstract type Geometry{T} end
 
@@ -155,6 +156,9 @@ end
 
 Geometry(xloc::CoordT, yloc::CoordT, zloc::CoordT, dt::Array{T,1}, nt::Array{Integer,1}, t::Array{T,1}) where {T<:Real} = GeometryIC{T}(xloc,yloc,zloc,dt,nt,t)
 
+# For easy 2D setup
+Geometry(xloc, zloc; kw...) = Geometry(xloc, 0 .* xloc, zloc; kw...)
+
 # Fallback constructors for non standard input types 
 
 # Constructor if nt is not passed
@@ -472,7 +476,64 @@ coords_from_set(S::OrderedSet{Tuple{T, T, T}}) where T = tuple([first.(S)], [get
 coords_from_keys(S::Vector{Tuple{T, T}}) where T = tuple([first.(S)], [[0f0]], [last.(S)])
 coords_from_keys(S::Vector{Tuple{T, T, T}}) where T = tuple([first.(S)], [getindex.(S, 2)], [last.(S)])
 
+"""
+    super_shot_geometry(Geometry)
+
+Merge all the sub-geometries `1:get_nsrc(Geometry)` into a single supershot geometry
+"""
 function super_shot_geometry(G::Geometry{T}) where T
     as_set = coords_from_set(as_coord_set(G))
     return GeometryIC{T}(as_set..., [G.dt[1]], [G.nt[1]], [G.t[1]])
 end
+
+
+###################### reciprocity ###############################
+
+"""
+    reciprocal_geom(sourceGeom, recGeom)
+
+Applies reciprocity to the par of geometries `sourceGeom` and `recGeom` where each source
+becomes a receiver and each receiver becomes a source.
+
+This method expects:
+- Both geometries to be In Core. If the geometries are OOC they will be converted to in core geometries
+- The metadata to be compatible. In details all the time sampling rates (dt) and recording times (t) must be the same
+- The source to be single point sources. This method will error if a simultaneous sources (multiple poisiton for a single source) are used.
+"""
+function reciprocal_geom(sGeom::GeometryIC{T}, rGeom::GeometryIC{T}) where T
+    # The geometry need to have the same recording and sampling times
+    @assert sGeom.dt == rGeom.dt
+    @assert sGeom.t == rGeom.t
+    @assert sGeom.nt == rGeom.nt
+    @assert allsame(sGeom.dt)
+    @assert allsame(sGeom.t)
+    @assert allsame(sGeom.nt)
+    # Make sure it's not simultaneous sources
+    if !all(length(x) == 1 for x in sGeom.xloc)
+        throw(GeometryException("Cannot apply reciprocity to simultaneous sources"))
+    end
+    # Curretnly only support geometry with all sources seeing the same receivers
+    if !allsame(rGeom.xloc)
+        throw(GeometryException("Currently expects all sources to see the same receivers (i.e OBNs)"))
+    end
+    # Reciprocal source geom
+    xsrc = convertToCell(rGeom.xloc[1])
+    if length(rGeom.yloc[1]) > 1
+        ysrc = convertToCell(rGeom.yloc[1])
+    else
+        ysrc = 0 .* xsrc
+    end
+    zsrc = convertToCell(rGeom.zloc[1])
+    sgeom = Geometry(xsrc, ysrc, zsrc; dt=rGeom.dt[1], t=rGeom.t[1])
+    # Reciprocal recc geom
+    xrec = Vector{T}([x[1] for x in sGeom.xloc])
+    yrec = Vector{T}([x[1] for x in sGeom.yloc])
+    zrec = Vector{T}([x[1] for x in sGeom.zloc])
+    rgeom = Geometry(xrec, yrec, zrec; dt=rGeom.dt[1], t=rGeom.t[1], nsrc=length(xsrc))
+    return sgeom, rgeom
+end
+
+function reciprocal_geom(sGeom::Geometry, rGeom::Geometry)
+    @warn "reciprocal_geom only supports in core geometries, converting"
+    return reciprocal_geom(Geometry(sGeom), Geometry(rGeom))
+end

src/TimeModeling/Types/ModelStructure.jl

@@ -8,7 +8,7 @@ export Model, PhysicalParameter, get_dt, size, origin, spacing, nbl
 
 abstract type AbstractModel{T, N} end
 
-struct DiscreteGrid{T, N}
+struct DiscreteGrid{T<:Real, N}
     n::NTuple{N, Int64}
     d::NTuple{N, <:T}
     o::NTuple{N, <:T}
@@ -340,18 +340,18 @@ where
 `nb`: Number of ABC points
 
 """
-function Model(n, d, o, m::Array{mT, N}; epsilon=nothing, delta=nothing, theta=nothing,
+function Model(d, o, m::Array{mT, N}; epsilon=nothing, delta=nothing, theta=nothing,
                phi=nothing, rho=nothing, qp=nothing, vs=nothing, nb=40) where {mT<:Real, N}
 
     # Currently force single precision
     m = convert(Array{Float32, N}, m)
     T = Float32
     # Convert dimension to internal types
-    n = NTuple{length(n), Int64}(n)
-    d = NTuple{length(n), Float32}(d)
-    o = NTuple{length(n), Float32}(o)
-    G = DiscreteGrid(n, d, o, nb)
-    
+    n = size(m)
+    d = tuple(Float32.(d)...)
+    o = tuple(Float32.(o)...)
+    G = DiscreteGrid{T, N}(n, d, o, nb)
+
     size(m) == n || throw(ArgumentError("Grid size $n and squared slowness size $(size(m)) don't match"))
 
     # Elastic
@@ -399,9 +399,9 @@ function Model(n, d, o, m::Array{mT, N}; epsilon=nothing, delta=nothing, theta=n
     return IsoModel{T, N}(G, m, rho)
 end
 
-Model(n, d, o, m::Array, rho::Array; nb=40) = Model(n, d, o, m; rho=rho, nb=nb)
-Model(n, d, o, m::Array, rho::Array, qp::Array; nb=40) = Model(n, d, o, m; rho=rho, qp=qp, nb=nb)
-Model(d, o, m; kw...) = Model(size(m), d, o, m; kw...)
+Model(n, d, o, m::Array, rho::Array; nb=40) = Model(d, o, m; rho=rho, nb=nb)
+Model(n, d, o, m::Array, rho::Array, qp::Array; nb=40) = Model(d, o, m; rho=rho, qp=qp, nb=nb)
+Model(n, d, o, m::Array; kw...) = Model(d, o, m; kw...)
 
 size(m::MT) where {MT<:AbstractModel} = size(m.G)
 origin(m::MT) where {MT<:AbstractModel} = origin(m.G)