OceanTTDs

TracerObservation{T,F<:Function}

Container for tracer observation data including measurement values, uncertainties, and source functions.

Fields

• t_obs::Vector{T}: Time points of observations
• y_obs::Vector{T}: Observed tracer concentrations
• σ_obs::Union{Nothing,Vector{T}}: Measurement uncertainties (optional)
• f_src::Union{Nothing,F}: Source function for boundary conditions (required for inversion)
• nobs::Int: Number of observations

Constructor

TracerObservation(t_obs, y_obs; σ_obs=nothing, f_src=nothing)

InversionResult{T}

Flexible container for inversion/optimization results that can handle different types of distributions and integration contexts.

Fields

• parameters::AbstractVector: Estimated parameters. Interpretation depends on method:
  • :inverse_gaussian: [Γ, Δ]
  • :inverse_gaussian_equalvars: [Γ]
  • :max_entropy: λ₁, λ₂, ..., λₙ
• obs_estimates::Vector{TracerEstimate{T}}: Updated TracerEstimate objects with predictions
• distribution::Union{Distribution, AbstractVector{<:Distribution}, Function, Nothing}: Fitted distribution(s) or template function
• integrator::Union{Nothing, Any}: Integration method/context used
• support::Union{Nothing, AbstractVector}: Discrete support points if applicable
• method::Symbol: Method used (e.g., :inversegaussian, :maxentropy, :inversegaussianequalvars)
• optimizer_output: Raw optimization results from the solver

invert_inverse_gaussian(observations; τ_max, integrator, C0=nothing, 
                       u0=[1e3,1e2], lower=[1.0,12.0], upper=[Inf,Inf],
                       warmstart=:anneal, sa_iters=50, lbfgs_iters=200)

Fit an Inverse Gaussian Transit Time Distribution (TTD) to tracer observations.

Mathematical Formulation

Estimates parameters Γ (mean) and Δ (width) of the Inverse Gaussian distribution:

f(τ; Γ, Δ) = √(Δ/(2πτ³)) exp(-(Δ(τ-Γ)²)/(2Γ²τ))

The optimization solves:

min_{Γ,Δ} Σᵢⱼ wᵢⱼ(yᵢⱼ - ŷᵢⱼ(Γ,Δ))²

where ŷᵢⱼ = ∫₀^τ_max f(τ)·gᵢ(tⱼ-τ)dτ + C0ᵢ is the convolved model prediction.

Arguments

• observations: TracerObservation(s) containing times, data, uncertainties, source functions
• τ_max: Maximum transit time for integration (scalar or vector per observation)
• integrator: Numerical integrator for convolution (scalar or vector per observation)
• C0: Additive constant offset (nothing→0, scalar, or vector per observation)
• u0: Initial parameter guess [Γ₀, Δ₀]
• lower, upper: Parameter bounds
• warmstart: Use simulated annealing initialization (:anneal or :none)
• sa_iters, lbfgs_iters: Iteration limits for SA and L-BFGS phases

Algorithm

1. Warm Start (if :anneal): Log-parameter simulated annealing to improve initial guess
2. Main Optimization: Bounded L-BFGS with box constraints via Fminbox
3. Model Evaluation: Final parameter estimates used to compute fitted values

Returns

InversionResult containing fitted parameters [Γ̂, Δ̂], observation estimates, and optimizer output.

Example

result = invert_inverse_gaussian(obs; τ_max=250_000.0, integrator=quad_integrator)
Γ_fit, Δ_fit = result.parameters

max_ent_inversion(observations; C0=nothing, prior_distribution, support)

Maximum Entropy inversion for estimate discrete Transit Time Distributions.

Mathematical Formulation

The Maximum Entropy Method solves the constrained optimization problem:

max S[p] = -Σᵢ pᵢ ln(pᵢ/mᵢ)    (maximize entropy)
subject to: Σᵢ pᵢ = 1           (normalization)
            Gp - d = 0.0     (fit observations)

where:

• p = [p₁, p₂, ..., pₙ] is the discrete TTD probability vector
• m = [m₁, m₂, ..., mₙ] is the prior distribution
• G is the forward operator (convolution matrix)
• d is the observation vector

This is solved via Lagrange multipliers, yielding:

pᵢ = mᵢ exp(-Σⱼ λⱼ Gⱼᵢ) / Z(λ)

where λ = [λ₁, λ₂, ..., λₘ] are Lagrange multipliers found by solving the nonlinear system:

∂S/∂λⱼ = Σᵢ Gⱼᵢ pᵢ(λ) - dⱼ = 0    for j = 1, ..., m

Algorithm

1. Setup: Initialize λ = 0 (default to prior distribution)
2. Levenberg-Marquardt: First-pass nonlinear solver with analytic Jacobian
3. Trust Region: Second-pass refinement for improved convergence
4. Distribution: Compute final p(λ) and estimated observations

Arguments

• observations: TracerObservation(s) with times, data, source functions
• prior_distribution: Prior probability vector m on support points
• support: Discrete transit time support points τ = [τ₁, τ₂, ..., τₙ]
• C0: Additive offset parameter (nothing→0, scalar, or vector)

Returns

InversionResult with λ parameters, discrete TTD probabilities, and solver output.

Example

τ_support = collect(0.0:100.0:10000.0)
m_prior = ones(length(τ_support)) / length(τ_support)  # uniform prior
result = max_ent_inversion(obs; prior_distribution=m_prior, support=τ_support)
ttd_probs = result.distribution

convolve(TTD, f_src, times; kwargs...)

Convolve a Transit Time Distribution with a source function at observation times.

This is the main convolution interface used by all optimization methods. Automatically dispatches to the appropriate method based on TTD type:

• Vector TTD: Direct discrete convolution with τ_support
• Discrete Distribution: Uses pdf() values with τ_support
• Continuous Distribution: Numerical integration with integrator

Arguments

• TTD: Transit time distribution (Vector, DiscreteUnivariateDistribution, or ContinuousUnivariateDistribution)
• f_src: Source function f(t)
• times: Vector of observation times
• Method-specific keywords (see individual methods)

Examples

# Discrete vector convolution (MaxEnt, TCM)
τs = collect(0:100:10000)
probs = [0.1, 0.3, 0.4, 0.2]  # Probability weights
result = convolve(probs, f_src, times; τ_support=τs)

# Continuous distribution convolution (Inverse Gaussian)  
dist = InverseGaussian(μ, λ)
result = convolve(dist, f_src, times; τ_max=10000.0, integrator=integrator)

convolve(TTD::Vector, f_src::Function, t_obs::Vector; 
         τ_support::Vector, C0::Real = 0.0)

Convolve discrete TTD probability weights with a source function.

Computes: ŷ(t) = ∑ⱼ TTD[j] * f_src(t - τⱼ) + C0

Used by MaxEnt and TCM optimization methods.

Arguments

• TTD: Vector of probability weights/probabilities
• f_src: Source function f(t)
• t_obs: Observation times
• τ_support: Transit time support points corresponding to TTD entries
• C0: Initial condition - tracer concentration at t=0 (default: 0.0)

Returns

Vector of convolution results at each observation time.

Example

τs = collect(1:100) 
TTD_weights = [0.1, 0.3, 0.4, 0.2]  # Discrete probabilities
f_src = t -> (t >= 0) ? 1.0 : 0.0    # Unit step function
t_obs = [10.0, 20.0, 30.0]

result = convolve(TTD_weights, f_src, t_obs; τ_support=τs)

convolve(TTD::DiscreteUnivariateDistribution, f_src::Function, times::Vector;
         τ_support::Vector, C0 = 0.0)

Convolve discrete distribution with source function using pdf values.

Arguments

• TTD: Discrete univariate distribution
• f_src: Source function f(t)
• times: Vector of observation times
• τ_support: Sorted, non-negative support points for evaluation
• C0: Initial condition - tracer concentration at t=0 (default: 0.0)

Returns

Vector with same length as times containing convolution results.

convolve(TTD::ContinuousUnivariateDistribution, f_src::Function, times::Vector;
         τ_max, integrator, C0 = 0.0)

Convolve continuous distribution with source function using numerical integration.

Used by Inverse Gaussian optimization method.

Arguments

• TTD: Continuous univariate distribution
• f_src: Source function f(t)
• times: Vector of observation times
• τ_max: Maximum transit time for integration domain [0, τ_max]
• integrator: Numerical integrator with .nodes and .weights fields
• C0: Initial condition - tracer concentration at t=0 (default: 0.0)

Returns

Vector with same length as times containing convolution results.

Example

using Distributions
dist = InverseGaussian(300.0, 600.0)  # Mean=300, shape=600
integrator = make_gausslegendre_integrator(100, 0.0, 10000.0)
result = convolve(dist, f_src, times; τ_max=10000.0, integrator=integrator)

tracer_observation(times, concentrations, source_function; uncertainties=nothing)

Create a TracerObservation from measurement data.

Arguments

• times: Observation time points
• concentrations: Measured tracer concentrations
• source_function: Boundary condition function f(t) for this tracer
• uncertainties: Measurement uncertainties (optional, defaults to 10% of concentrations)

Returns

TracerObservation object ready for inversion

Example

# CFC-12 with exponential atmospheric growth
cfc12_source = t -> (t >= 1930) ? exp(0.08 * (t - 1930)) : 0.0
obs = tracer_observation(times, cfc12_data, cfc12_source)

# Multiple tracers for joint inversion
obs1 = tracer_observation(times1, cfc11_data, cfc11_source)  
obs2 = tracer_observation(times2, sf6_data, sf6_source)
result = max_ent_inversion([obs1, obs2]; support=τ_support, prior_distribution=prior)

uniform_prior(support)

Create a uniform prior distribution over the given support points.

Arguments

• support: Transit time support points

Returns

Uniform probability vector (all values equal, sum to 1)

Example

τ_support = collect(0:50:2000)
prior = uniform_prior(τ_support)
result = fit_ttd(obs, method=:max_entropy, support=τ_support, prior_distribution=prior)