PrISM.NonlinearAlg Method

NonlinearAlg(; alg = LevenbergMarquardt, μ = 1.0)

returns nl_cache that specifies which non linear solver to use for the inverse problem

Keyword Arguments

• alg: NonlinearSolve[@ref] algorithm to be used, defaults to LevenbergMarquardt

• μ : regularization weight

source

PrISM.Occam Method

Occam(;μgrid= [0.01, 1e6])

source

PrISM.OptAlg Method

OptAlg(; alg = LBFGS, μ = 1.0, kwargs...)

returns opt_cache that specifies which Optimization.jl solver to use for the inverse problem

Keyword Arguments

• alg: Optimization[@ref] algorithm to be used, defaults to LBFGS

• μ : regularization weight

source

PrISM.forward! Method

forward(resp::DCResponse, m::DCModel, locs)

overwrite DCResponse resp for the given DCModel m at the electrode locations defined by locs

Arguments

• resp : DCResponse to be overwritten

• m : DCModel for forward response

• locs : NamedTuple containing locations of electric dipole defined by srcs and electrode locations defined by recs

Optional arguments

• params : NamedTuple containing configuration for forward calculation, contains
  • hankel_filter : Hankel filter to be used, defaults to Kerry Key's 101 point filter

source

PrISM.forward! Method

forward!(resp::MTResponse, m::MTModel, ω)

overwrites MTResponse resp for the given MTModel m at the frequencies ω

Arguments

• resp : MTResponse to be overwritten

• m : MTModel for forward response

• ω : angular frequencies where the response are estimated (=2π/Time period )

source

PrISM.forward! Method

forward!(resp::LWResponse, m::LWModel, t)

overwrites LWResponse for the given LWModel m at the periods defined by t

Arguments

• resp : LWResponse to be overwritten

• m : LWModel for forward response

• t : periods for which forward response is to be calculated

Optional arguments

• params : NamedTuple containing configuration for forward calculation, contains

  • mode : dispersion curve mode, defaults to 0 implying fundamental node

  • dc : step size used to search for grid space solution of propagator matrix, defaults to 0.001

  • dt : Δt to be used for calculation of group velocity from phase velocity

  • type : Val type variable that specifies whether to calculate group velocity or phase. Available options are :

    • Val(:phase) : phase velocity (default)

    • Val(:group) : group velocity

source

PrISM.forward! Method

forward!(resp::RWResponse, m::RWModel, t)

overwrites RWResponse for the given RWModel m at the periods defined by t

Arguments

• resp : RWResponse to be overwritten

• m : LWModel for forward response

• t : periods for which forward response is to be calculated

Optional arguments

• params : NamedTuple containing configuration for forward calculation, contains

  • mode : dispersion curve mode, defaults to 0 implying fundamental node

  • dc : step size used to search for grid space solution of propagator matrix, defaults to 0.001

  • dt : Δt to be used for calculation of group velocity from phase velocity

  • type : Val type variable that specifies whether to calculate group velocity or phase. Available options are :

    • Val(:phase) : phase velocity (default)

    • Val(:group) : group velocity

source

PrISM.get_T Method

get T for one wavelength

source

PrISM.get_appres Method

get_appres(Z, ω): returns the ρₐ for impedance

source

PrISM.get_phase Method

get_phase(Z): returns the phase for impedance

source

PrISM.get_schlumberger_array Method

get_schlumberger_array(a, n)

returns electrode spacings for Schlumberger array :

M ←na→ A ←a→ B ←na→ N

Arguments

• a : electrode spacing AB

• n : vector of numbers corresponding to spacing

Usage : TODO

source

PrISM.get_wenner_array Method

get_wenner_array(a)

returns electrode spacings for Wenner array :

M ←a→ A ←a→ B ←a→ N

Arguments

• a : vector of electrode spacing AB

Usage

TODO

source

PrISM.hankel_transform_and_interpolation Method

hankel_transform_and_interpolation(r_max, r_min, fn, hankel_filter)

performs Hankel transform to evaluate the integral of the form:

$$\int fn(\lambda )J₀(\lambda r)d\lambda$$

and returns an interpolating function for any r ∈ [r_min, r_max]

Arguments

• r_max : maximum r

• r_min : minimum r

• fn : kernel function

• hankel_filter : hankel filter to be used, contains base, J₀ and J₁

returns an interpolated function that outputs the integral using Hankel transform at any r

source

PrISM.inverse! Method

function inverse!(mₖ, robs, vars, alg_cache::occam_cache; 
        W, L, max_iters, χ2, response_fields, model_trans_utils,
        response_trans_utils, mᵣ, reg_term, verbose
    ):

updates mₖ using occam iteration to fit robs within a misfit of χ2, by default set to 1.0

Arguments

• mₖ: Initial model guess, will be updated during the inverse process

• robs: response to invert for

• vars: variables required for forward modeling, eg., ω for MT

• alg_cache: deterimines the algorithm to be performed for inversion

Keyword arguments

• W: Weight matrix, defaults to identity matrix I

• L: Regularization matrix, defaults to derivative matrix, given by ∂(@ref)

• max_iters: maximum number of iterations, defaults to 30

• χ2: target misfit, defaults to 1.

• `response_fields: choose data of response to perform inversion on, eg., ρₐ for MT, by default chooses all the data (ρₐ and ϕ)

• model_trans_utils:transform_utils: conversion to and from computational domain, defaults to log_sigmoid_tf

• response_trans_utils: for scaling the response parameters,

• mᵣ: model in physical domain to be regularized against

• reg_term: (For internals, only used by Occam) When model in physical domain does not exist, reg_term helps, eg. case of RTO-TKO

• verbose: whether to print updates after each iteration, defaults to true

Returns:

return message in the form of return_code and updates mₖ in-place.

Usage:

using LinearAlgebra
h = [1000.0, 1000.0] # m
ρ = log10.([100.0, 10.0, 1000.0]) # Ωm
m = MTModel(ρ, h)
T = 10 .^ (range(-2, 2; length=57))
ω = 2π ./ T
nω = length(T)
r_obs = forward(m, ω)
m_occam = MTModel(fill(2.0, 50), collect(range(0, 5e3; length=49)))
err_ρ = 0.1 .* r_obs.ρₐ
err_ϕ = asin(0.01) * 180 / π .* ones(length(ω))
err_resp = diagm(inv.([err_ρ..., err_ϕ...])) .^ 2
ret_code = inverse!(m_occam, r_obs, ω, Occam(; μgrid=[1e-2, 1e6]);
    W=err_resp, χ2=1.0, max_iters=50, verbose=true)

# output

1: Works golden section search: μ= 2559.998204059276, χ²= 22.867791875474534
2: Works golden section search: μ= 55.11264178720504, χ²= 16.89701371365216
3: Works golden section search: μ= 36.07944656829589, χ²= 12.58408431380094
4: Works golden section search: μ= 26.62727840194216, χ²= 9.43690878244996
5: Works golden section search: μ= 16.28674219250041, χ²= 7.1269845988502025
6: Works golden section search: μ= 8.309723956102522, χ²= 5.41891016151562
7: Works golden section search: μ= 4.116854777415263, χ²= 4.142014283721758
8: Works golden section search: μ= 2.115717196002266, χ²= 3.1786724887517708
9: Works golden section search: μ= 1.173804966446848, χ²= 2.447214296425664
10: Works golden section search: μ= 0.7087166896720183, χ²= 1.8893177576767075
11: Works golden section search: μ= 0.45988622588032063, χ²= 1.4623667912514455
12: Works golden section search: μ= 0.3163949884082966, χ²= 1.1347540808692684
13: Works golden section search: μ= 0.22874422627086552, χ²= 0.8828031760948966
return_code{MTModel{Vector{Float64}, Vector{Float64}}}(true, (μ = 0.22874422627086552,), , 1.0, 0.8828031760948966)

source

PrISM.linsolve! Method

linsolve!: Performs inv(B)*y using LinearSolve.jl

source

PrISM.occam_step! Method

function occam_step!(mₖ₊₁::model,
    respₖ₊₁::response,
    vars::Union{AbstractVector{Float32}, AbstractVector{Float64}},
    χ2::Union{Float64, Float32},
    μgrid::Vector{Float64},
    lin_utils::linear_utils,
    inv_utils::inverse_utils,
    model_trans_utils::transform_utils,
    response_trans_utils::NamedTuple,
    linsolve_prob::LinearSolve.LinearCache;
    model_fields::Vector{Symbol}= [k for k ∈ fieldnames(typeof(mₖ₊₁))],
    response_fields::Vector{Symbol}= [k for k ∈ fieldnames(typeof(respₖ₊₁))],
    verbose= false
    ):

performs a single step of occam inversion, using golden line search.

Arguments

• mₖ: Initial model guess, will be updated during the inverse process

• respₖ₊₁: response to invert for

• vars: variables required for forward modeling, eg., ω for MT

• alg_cache: deterimines the algorithm to be performed for inversion

• W: Weight matrix, defaults to identity matrix I

• L: Regularization matrix, defaults to discretized derivative matrix given by ∂(@ref)

• max_iters: maximum number of iterations, defaults to 30

• χ2: target misfit, defaults to 1.0

• `response_fields: choose data of response to perform inversion on, eg., ρₐ for MT, by default chooses all the data (ρₐ and ϕ)

• model_trans_utils: conversion to and from computational domain,

• response_trans_utils: NamedTuple containing transform_utils for scaling different response parameters,

• mᵣ: model in physical domain to be regularized against

• reg_term: (For internals) When model in physical domain does not exist, reg_term helps, eg. case of RTO-TKO

• verbose: whether to print updates after each iteration, defaults to true

source

PrISM.χ² Method

χ²(dcal::T, dobs::T; W): returns a chi-squared error between the observed and the calculated data. W can optionally be passed to weigh points differently.

source

PrISM.∂ Method

∂(n): returns a nxn matrix for a 1D finite difference stencil using 2 points.

source

SubsurfaceCore.forward Method

forward(m::DCModel, locs)

returns DCResponse for the given DCModel m at the electrode locations defined by locs

Arguments

• m : DCModel for forward response

• locs : NamedTuple containing locations of electric dipole defined by srcs and electrode locations defined by recs

Optional arguments

• params : NamedTuple containing configuration for forward calculation, contains
  • hankel_filter : Hankel filter to be used, defaults to Kerry Key's 101 point filter

source

SubsurfaceCore.forward Method

forward(m::LWModel, t)

returns LWResponse for the given LWModel m at the periods defined by t

Arguments

• m : LWModel for forward response

• t : periods for which forward response is to be calculated

Optional arguments

• params : NamedTuple containing configuration for forward calculation, contains

  • mode : dispersion curve mode, defaults to 0 implying fundamental node

  • dc : step size used to search for grid space solution of propagator matrix, defaults to 0.001

  • dt : Δt to be used for calculation of group velocity from phase velocity

  • type : Val type variable that specifies whether to calculate group velocity or phase. Available options are :

    • Val(:phase) : phase velocity (default)

    • Val(:group) : group velocity

source

SubsurfaceCore.forward Method

forward(m::MTModel, ω)

returns MTResponse for the given MTModel m at the frequencies ω

Arguments

• m : MTModel for forward response

• ω : angular frequencies where the response are estimated (=2π/Time period )

source

SubsurfaceCore.forward Method

forward(m::RWModel, t)

returns RWResponse for the given RWModel m at the periods defined by t

Arguments

• m : RWModel for forward response

• t : periods for which forward response is to be calculated

Optional arguments

• params : NamedTuple containing configuration for forward calculation, contains

  • mode : dispersion curve mode, defaults to 0 implying fundamental node

  • dc : step size used to search for grid space solution of propagator matrix, defaults to 0.001

  • dt : Δt to be used for calculation of group velocity from phase velocity

  • type : Val type variable that specifies whether to calculate group velocity or phase. Available options are :

    • Val(:phase) : phase velocity (default)

    • Val(:group) : group velocity

source

SubsurfaceCore.stochastic_inverse Method

stochastic_inverse(r_obs::resp1,
        err_resp::resp2,
        vars,
        alg_cache::rto_cache;
        model_trans_utils::NamedTuple=(m=sigmoid_tf, h=lin_tf),
        response_trans_utils::NamedTuple=(ρₐ=lin_tf, ϕ=lin_tf)) where {
        resp1 <: AbstractGeophyResponse, resp2 <: AbstractGeophyResponse}

Returns

`AbstractMCMC.jl Chain` containing the vectors used for performing samping. Note that all the variables will be named `var_inf`, for variable inferred. 
The vector in each sample will be all the fields of the model being inferred on concatenated together.

Arguments

• r_obs: response that needs to inverted for

• err_resp: response variable containing the errors associated with observed response

• vars: variables that need to be passed into the forward function along with model to generate a response

• alg_cache: to tell the compiler what type of stochastic inversion method is to be used

Keyword Arguments

• model_trans_utils: A named tuple containing transform_utils for the fields of model that need to be scaled/modified. If not provided for any model field, the field won't be modified

• response_trans_utils: for scaling the response parameters

source

PrISM.DCModel Type

create a model type for a given resistivity distribution that can be used to calculate forward response for 1d DC

Usage

m = [4.0, 2.0, 3.0]
h = [1000.0, 1000.0]
model = DCModel(m, h)
print(model)

# output

1D DC Resistivity Model
┌───────┬───────┬─────────┬─────────┐
│ Layer │ log-ρ │       h │       z │
│       │  [Ωm] │     [m] │     [m] │
├───────┼───────┼─────────┼─────────┤
│     1 │  4.00 │ 1000.00 │ 1000.00 │
│     2 │  2.00 │ 1000.00 │ 2000.00 │
│     3 │  3.00 │       ∞ │       ∞ │
└───────┴───────┴─────────┴─────────┘

source

PrISM.DCModelDistribution Type

mutable struct DCModelDistribution{T1<: Union{Distribution, AbstractArray}, T2<: Union{Distribution, AbstractArray}}
    m::T1
    h::T2
end

creates a placeholder to store the Distributions.jl samplers for a priori

source

PrISM.DCResponseDistribution Type

struct DCResponseDistribution{T1<: Union{Function, Nothing}, T2<: Union{Function, Nothing}}
    ρₐ::T1
    ϕ::T2
end

create a placeholder to store functions to obtain Distributions.jl samplers for the likelihood function

source

PrISM.LWModel Type

create a model type for a given resistivity distribution that can be used to calculate forward response for 1d MT

Usage

m = [4e3, 3e3, 3.5e3]
h = [1000.0, 1000.0]
d = [3e3, 3e3, 3e3]
model = LWModel(m, h, d)

# output

1D Love Wave Model
┌───────┬─────────┬─────────┬─────────┬──────────┐
│ Layer │      vₛ │       h │       z │        ρ │
│       │  [km/s] │     [m] │     [m] │ [kg⋅m⁻³] │
├───────┼─────────┼─────────┼─────────┼──────────┤
│     1 │ 4000.00 │ 1000.00 │ 1000.00 │   3000.0 │
│     2 │ 3000.00 │ 1000.00 │ 2000.00 │   3000.0 │
│     3 │ 3500.00 │       ∞ │       ∞ │   3000.0 │
└───────┴─────────┴─────────┴─────────┴──────────┘

source

PrISM.LWModelDistribution Type

mutable struct LWModelDistribution{T1<: Union{Distribution, AbstractArray}, T2<: Union{Distribution, AbstractArray}} # where T1,T2 
    m::T1
    h::T2
    ρ::T3
end

creates a placeholder to store the Distributions.jl samplers for a priori

source

PrISM.MTModel Type

create a model type for a given resistivity distribution that can be used to calculate forward response for 1d MT

Usage

m = [4.0, 2.0, 3.0]
h = [1000.0, 1000.0]
model = MTModel(m, h)

# output

1D MT Model
┌───────┬───────┬─────────┬─────────┐
│ Layer │ log-ρ │       h │       z │
│       │  [Ωm] │     [m] │     [m] │
├───────┼───────┼─────────┼─────────┤
│     1 │  4.00 │ 1000.00 │ 1000.00 │
│     2 │  2.00 │ 1000.00 │ 2000.00 │
│     3 │  3.00 │       ∞ │       ∞ │
└───────┴───────┴─────────┴─────────┘

source

PrISM.MTModelDistribution Type

mutable struct MTModelDistribution{T1<: Union{Distribution, AbstractArray}, T2<: Union{Distribution, AbstractArray}}
    m::T1
    h::T2
end

creates a placeholder to store the Distributions.jl samplers for a priori

source

PrISM.MTResponseDistribution Type

struct MTResponseDistribution{T1<: Union{Function, Nothing}, T2<: Union{Function, Nothing}}
    ρₐ::T1
    ϕ::T2
end

create a placeholder to store functions to obtain Distributions.jl samplers for the likelihood function

source

PrISM.RWModel Type

create a model type for a given resistivity distribution that can be used to calculate forward response for 1d MT

Usage

m = [4e3, 3e3, 3.5e3]
h = [1000.0, 1000.0]
d = [3e3, 3e3, 3e3]
vp = m .* 1.72
model = RWModel(m, h, d, vp)

# output

1D Rayleigh Wave Model
┌───────┬─────────┬─────────┬─────────┬──────────┬────────┐
│ Layer │      vₛ │       h │       z │        ρ │     vₚ │
│       │  [km/s] │     [m] │     [m] │ [kg⋅m⁻³] │ [km/s] │
├───────┼─────────┼─────────┼─────────┼──────────┼────────┤
│     1 │ 4000.00 │ 1000.00 │ 1000.00 │   3000.0 │ 6880.0 │
│     2 │ 3000.00 │ 1000.00 │ 2000.00 │   3000.0 │ 5160.0 │
│     3 │ 3500.00 │       ∞ │       ∞ │   3000.0 │ 6020.0 │
└───────┴─────────┴─────────┴─────────┴──────────┴────────┘

source

PrISM.RWModelDistribution Type

mutable struct RWModelDistribution{T1<: Union{Distribution, AbstractArray}, T2<: Union{Distribution, AbstractArray}} # where T1,T2 
    m::T1
    h::T2
    ρ::T3
    vp:T4
end

creates a placeholder to store the Distributions.jl samplers for a priori

source

PrISM.SurfaceWaveResponseDistribution Type

struct SurfaceWaveResponseDistribution{T1<: Union{Function, Nothing}}
    c::T1
end

creates a placeholder to store the Distributions.jl samplers for a priori

source

PrISM.inverse_utils Type

struct inverse_utils: contains the utilities for inversion, once initialized, will not be updated in the inversion iterations D: second derivative operator, W: weight matrix, dobs: data response to be inverted for.

source

PrISM.linear_utils Type

struct linear_utils: contains the utilities for linearizing the forward model => mₖ:model, Fₖ: Forward response at mₖ, Jₖ: Jacobian at mₖ.

source

PrISM.nl_cache Type

nl_cache: specifies the inverse algorithm while having a cache.

source

PrISM.occam_cache Type

occam_cache: specifies the inverse algorithm while having a cache.

source

PrISM.opt_cache Type

nl_cache: specifies the inverse algorithm while having a cache.

source

PrISM.return_code Type

struct return_code: contains the information if the inversion was successful

source

PrISM.rto_cache Type

rto_cache(m₀, μgrid, alg, max_iters, n_samples, χ2, response_fields, L, verbose)

returns rto_cache that specifies the algorithm to be used for stochastic inversion using RTO-TKO

Arguments

• m₀: model, usually a starting value to start inversion but mostly to allocate memory space

• μgrid: prior space of μ

• alg: type of algorithm for optimization step, choices include (occam_cache)[@ref], (nl_cache)[@ref], and (opt_cache)[@ref]

• max_iters: maximum iterations for the optimization scheme

• n_samples: number of samples to be obtained. Note: If there are any spurious samples (say NaN values), they will be automatically deleted, reducing the sample size

• χ2: target misfit

• response_fields: choose data of response to perform inversion on, eg., ρₐ for MT, by default chooses all the data (ρₐ and ϕ)

• L: Regularization matrix, defaults to discretized derivative matrix given by ∂(@ref)

• verbose: to print results or not

source