ODINN.jl API

AbstractAdjointMethod

Abstract type representing the flavor of AD and adjoint to be used to compute the gradient of the cost function. There are two parts where one can play with how the gradient is propagated: the iceflow model VJP and the adjoint of the ODE solver. The VJP of the iceflow model can be computed using either AD (Zygote or Enzyme), the discrete, or the continuous adjoint of the iceflow model. As for the computation of the adjoint of the ODE solution, it can be handled by SciMLSensitivity, or computed using the adjoint implemented in ODINN.

AbstractVJPMethod

Abstract type representing the flavor of AD to be used to compute the VJP inside the gradient of the cost function.

ContinuousAdjoint{
    F <: AbstractFloat,
    I <: Integer,
    VJP <: AbstractVJPMethod,
    MBVJP <: AbstractVJPMethod,
} <: AbstractAdjointMethod

Continuous adjoint of SIA2D with manual implementation of the backward in the ODE scheme.

Fields

• VJP_method::VJP: Type of AbstractVJPMethod used to compute VJPs inside adjoint calculation.
• solver::Any: The solver to be used for adjoint.
• reltol::F: Relative tolerance to be used in the ODE solver of the adjoint.
• abstol::F: Absolute tolerance to be used in the ODE solver of the adjoint.
• dtmax::F: Maximum time step to be used in the ODE solver of the adjoint.
• interpolation: Interpolation method to be used to interpolate the variables in the computation of the adjoint. Currently only :Linear is supported.
• n_quadrature::I: Number of nodes used in the Gauss quadrature for the numerical integration of the loss function.
• MB_VJP::MBVJP: Type of AbstractVJPMethod used to compute the MB VJP inside adjoint calculation.

ContinuousVJP{ADTYPE <: DI.AbstractADType} <: AbstractVJPMethod

Continuous manual implementation of the VJP required inside the adjoint calculation. It relies in the continuous expresion for the adjoint operation based on the functional formula of the forward PDE.

Fields

• regressorADBackend::ADTYPE: Specifies the AD backend used for the laws when their associated VJPs functions are not provided. The type parameter ADTYPE must be a subtype of DI.AbstractADType.

Constructor

• The default constructor allows specifying the backend via the regressorADBackend keyword argument, defaulting to DI.AutoMooncake().

DiffusivityRegularization(; reg = TikhonovRegularization())

Regularization for diffusivity fields using a specified spatial operator.

Keyword Arguments

• reg::AbstractSimpleRegularization = TikhonovRegularization(): Spatial regularization operator applied to diffusivity.

DiscreteAdjoint{
    VJP <: AbstractVJPMethod,
    MBVJP <: AbstractVJPMethod,
} <: AbstractAdjointMethod

Discrete adjoint of SIA2D with manual implementation of the backward in the ODE scheme.

Fields

• VJP_method: Type of AbstractVJPMethod used to compute VJPs inside adjoint calculation.
• MB_VJP::MBVJP: Type of AbstractVJPMethod used to compute the MB VJP inside adjoint calculation.

DiscreteVJP{ADTYPE <: DI.AbstractADType} <: AbstractVJPMethod

Discrete manual implementation of the VJP required inside the adjoint calculation. This implements the pullback function for the function to differentiate.

Fields

• regressorADBackend::ADTYPE: Specifies the AD backend used for the laws when their associated VJPs functions are not provided. The type parameter ADTYPE must be a subtype of DI.AbstractADType.

Constructor

• The default constructor allows specifying the backend via the regressorADBackend keyword argument, defaulting to DI.AutoMooncake().

Struct to provide a dummy gradient. It does not have to be the true gradient. Mainly used to test that the optimization pileline works independenly of the gradient calculation.

DummyAdjoint

Fields:

• grad::Function: In-place function f(du, u; kwargs) that fills the first argument du with the gradient values.

Enzyme implementation of VJP used inside the adjoint calculation.

EnzymeVJP

FunctionalModel <: TrainableModel

Abstract type representing functional learnable components of the model. This is a subtype of TrainableModel. Typically used for functional inversions.

GlacierWideInv{
    ComponentVectorType <: ComponentVector
} <: PerGlacierModel

Per glacier invertible parameter container. GlacierWideInv wraps a ComponentVector (θ) that stores one scalar parameter per glacier and implements the PerGlacierModel interface used by the inversion machinery.

Fields

• θ::ComponentVectorType: The per glacier parameter vector (one scalar value per glacier).

Constructor

GlacierWideInv(
    params::Sleipnir.Parameters,
    glaciers::Vector{<: AbstractGlacier},
    var::Symbol;
    minval::Union{Nothing, Float64} = nothing,
    maxval::Union{Nothing, Float64} = nothing,
)

Arguments

• params::Sleipnir.Parameters: Parameters struct.
• glaciers::Vector{<: AbstractGlacier}: Vector of AbstractGlacier. The i-th entry in θ corresponds to glaciers[i].
• var::Symbol: Symbol naming the field on each glacier to use as the initial value. Only :A and :C are supported.
• minval::Union{Nothing, Float64}: Lower bound of the parameterization. Defaults to params.physical.minA for :A; unused for :C, whose lower bound is always zero.
• maxval::Union{Nothing, Float64}: Upper bound of the parameterization. Defaults to params.physical.maxA for :A and params.physical.maxC for :C.

Example

GlacierWideInv(params, glaciers, :A)
GlacierWideInv(params, glaciers, :C)

GriddedInv{
    ComponentVectorType <: ComponentVector
} <: PerGlacierModel

Per glacier invertible parameter container. GriddedInv wraps a ComponentVector (θ) that stores one matrix parameter per glacier and implements the PerGlacierModel interface used by the inversion machinery.

Fields

• θ::ComponentVectorType: The per glacier parameter vector (one matrix per glacier).

Constructor

GriddedInv(
    params::Sleipnir.Parameters,
    glaciers::Vector{<: AbstractGlacier},
    var::Symbol,
)

Arguments

• params::Sleipnir.Parameters: Parameters struct.
• glaciers::Vector{<: AbstractGlacier}: Vector of AbstractGlacier. The i-th entry in θ corresponds to glaciers[i].
• var::Symbol: Symbol naming the field on each glacier to use as the initial value.

Example

GriddedInv(params, glaciers, :A)

mutable struct Hyperparameters{F <: AbstractFloat, I <: Integer} <: AbstractParameters

A mutable struct that holds hyperparameters for training a machine learning model.

Keyword arguments

• current_epoch::I: The current epoch number.
• current_minibatch::I: The current minibatch number.
• loss_history::Vector{F}: A vector storing the history of loss values.
• optimizer::Union{Optim.FirstOrderOptimizer, Optimisers.AbstractRule, Vector{Optim.FirstOrderOptimizer}, Vector{Any}}: The optimizer used for training.
• loss_epoch::F: The loss value for the current epoch.
• epochs::I: The total number of epochs for training.
• batch_size::I: The size of each minibatch.

Hyperparameters(;
    current_epoch::Int64 = 1,
    current_minibatch::Int64 = 1,
    loss_history::Vector{Float64} = Vector{Float64}(),
    optimizer::Union{Optim.FirstOrderOptimizer, Vector{Optim.FirstOrderOptimizer}, Vector{Any}} = BFGS(initial_stepnorm = 0.001),
    loss_epoch::Float64 = 0.0,
    epochs::Union{Int64, Vector{Int64}} = 50,
    batch_size::Int64 = 15

)

Constructs a Hyperparameters object with the specified parameters.

Arguments

• current_epoch::Int64: The current epoch number. Defaults to 1.
• current_minibatch::Int64: The current minibatch number. Defaults to 1.
• loss_history::Vector{Float64}: A vector to store the history of loss values. Defaults to an empty vector.
• optimizer::Union{Optim.FirstOrderOptimizer, Vector{Optim.FirstOrderOptimizer}, Vector{Any}}: The optimizer to be used. Defaults to BFGS(initial_stepnorm=0.001).
• loss_epoch::Float64: The loss value for the current epoch. Defaults to 0.0.
• epochs::Int64: The total number of epochs. Defaults to 50.
• batch_size::Int64: The size of each minibatch. Defaults to 15.

Returns

• A Hyperparameters object initialized with the provided values.

InitialCondition{
    ComponentVectorType <: ComponentVector
} <: PerGlacierModel

Per glacier initial condition container. InitialCondition wraps a ComponentVector (θ) that stores one matrix per glacier and implements the InitialCondition interface used by the inversion machinery.

Fields

• θ::ComponentVectorType: The per glacier matrix of initial condition.

Constructor

InitialCondition(
    params::Sleipnir.Parameters,
    glaciers::Vector{<: AbstractGlacier},
    initialization::Symbol = :Farinotti19,
)

Arguments

• params::Sleipnir.Parameters: Parameters struct.
• glaciers::Vector{<: AbstractGlacier}: Vector of AbstractGlacier. The i-th entry in θ corresponds to glaciers[i].
• initialization::Symbol: Symbol providing the way the initial condition should be initialized.

Example

InitialCondition(params, glaciers, :Farinotti19)

InitialThicknessRegularization(t₀::AbstractFloat; reg::AbstractSimpleRegularization = TikhonovRegularization())

A composite regularization type designed for initial ice thickness. It combines a simple spatial regularization (e.g., TikhonovRegularization) with a reference initial time.

Keyword Arguments

• reg::AbstractSimpleRegularization = TikhonovRegularization(): The spatial regularization operator applied to the initial field. By default, a Tikhonov (Laplacian-based) regularization is used.
• t₀::AbstractFloat = 1994.0: The reference initial time (e.g., year) at which the regularization applies.

mutable struct Inversion{MODEL, CACHE, GLACIER, RES} <: Simulation

An object representing an inversion simulation. It can involve at the same time a classical inversion and a functional inversion (i.e. the inversion of a function using some data-driven regressor).

Fields

• model::Sleipnir.Model: The model used for the simulation.
• glaciers::Vector{Sleipnir.AbstractGlacier}: A vector of glaciers involved in the simulation.
• parameters::Sleipnir.Parameters: The parameters used for the simulation.
• results::ODINN.Results: A ODINN.Results instance to store the results of the inversion and of the forward simulations.

function Inversion(
    model::M,
    glaciers::Vector{G},
    parameters::P
) where {G <: Sleipnir.AbstractGlacier, M <: Sleipnir.Model, P <: Sleipnir.Parameters}

Constructor for Inversion struct with glacier model information, glaciers, and parameters.

Arguments

• model::Sleipnir.Model: The model used for the simulation.
• glaciers::Vector{G}: A vector of glaciers involved in the simulation.
• parameters::Sleipnir.Parameters: The parameters used for the simulation.

Returns

• Inversion: A new instance of the Inversion struct.

InversionBinder{FI <: Inversion, CA <: ComponentArray}

Struct used for the binding with SciMLSensitivity. It is defined as a SciMLStructure and it contains the inversion structure and the vector of parameters to differentiate.

Fields

• simulation::FI: Inversion instance.
• θ::CA: ComponentArray that contains the parameters used to differentiate the iceflow.

L2Sum{I <: Integer} <: AbstractSimpleLoss

Struct that defines an L2 sum loss. The sum is defined only on pixels inside the glacier. The parameter distance controls the pixels that should be used to compute the sum. In order for a pixel to be used, it should be at least at distance from the glacier border. The mask defining the glacier borders are computed using the ground truth ice thickness.

$loss(a,b) = sum_{i\in\text{mask}} (a[i]-b[i])^2 / normalization$

Fields

• distance::I: Distance to border.

LogSum{I <: Integer, F <: AbstractFloat} <: AbstractSimpleLoss

Struct that defines a Logarithmic sum loss.

$loss(a,b) = log^2( (a + ϵ) / (b + ϵ) ) / normalization$

Fields

• distance::I: Distance to border.
• ϵ::F: Epsilon used inside the loss function to handle division by zero and log(0). It somehow represents the minimum value the loss function should be sensible to.

LossAvgV{F <: AbstractFloat, L <: AbstractSimpleLoss} <: TimeAggregatedLoss

A loss function that penalizes the difference between predicted and observed time-averaged glacier surface velocities.

This loss type computes a time-weighted average of predicted velocities over a specified time interval and compares it against reference velocity observations. It is particularly useful for constraining glacier flow dynamics when velocity data is annual for example (single snapshot inversions).

Fields

• velocityProduct::Symbol = :Millan22: Velocity product to use as reference. Defaults to :Millan22.
• loss::L = L2Sum(): The underlying loss function type used to compare predicted and reference velocities
• component::Symbol = :xy: Which velocity component(s) to use in the loss:
  • :xy: Compare x and y velocity components separately (sum of losses)
  • :abs: Compare absolute velocity magnitude
• step::F = 1/12: Time stepping for velocity aggregation (default: 1 month in yearly units)

Details

The loss computation involves:

1. Creating a time grid from t1 to t2 with spacing lossType.step
2. Computing predicted velocities at each time step via ice thickness predictions
3. Time-averaging the velocities with weights proportional to time intervals
4. Comparing the averaged velocity to reference observations using the specified loss function

LossDhdt <: TimeAggregatedLoss

A loss function that penalizes the difference between predicted and observed glacier surface elevation change rates (dh/dt).

This loss works with time-aggregated quantities, comparing the mean rate of height change computed from ice thickness predictions against reference dh/dt observations between the beginning and the end of a specified time interval.

Details

The loss is computed as: L = (dhdtpred - dhdtref)²

where:

• dhdt_pred: Predicted mean rate of height change (computed from model ice thickness outputs)
• dhdt_ref: Reference/observed rate of height change from data
• The rate is computed using masked ice thickness differences by masking out pixels without ice based on the ice thickness at the beginning of the time window

LossH{L <: AbstractSimpleLoss} <: AbstractLoss

Struct that defines the ice thickness loss.

Fields

• loss::L: Type of loss to use for the ice thickness. Default is L2Sum().

LossHV{
    F <: AbstractFloat,
    LH <: AbstractLoss,
    LV <: AbstractLoss,
} <: AbstractLoss

Struct that defines the ice thickness and ice velocity loss. It consists of two fields that define the ice thickness and ice velocity loss. It also has a scaling coefficient that balances the ice velocity term in the loss.

$loss(\hat H,H) = loss_H(\hat H,H) + scaling * loss_V(\hat V,V)$

with $\hat V$ computed from $\hat H$ for the SIA.

Fields

• hLoss::LH: Type of loss to use for the ice thickness. Default is LossH().
• vLoss::LV: Type of loss to use for the ice velocity. Default is LossV().
• scaling::F: Scaling of the ice velocity term. Default is 1.

LossV{L <: AbstractSimpleLoss} <: AbstractLoss

Struct that defines the ice velocity loss.

Fields

• loss::L: Type of loss to use for the ice velocity. Default is L2Sum().
• component::Symbol: Component of the velocity field used in the loss. Options include :xy for both x and y component, and :abs for the norm/magnitude of the velocity.
• scale_loss::Bool: Whether to scale the loss function with the reference ice velocity magnitude.

MatrixCacheGlacierId <: Cache

A cache structure for storing a matrix value (Float64 2D array) along with their associated vector-Jacobian products (VJP). It also stores the glacier ID as an integer. This is typically used to invert a spatially varying field per glacier. Fields:

• value::Array{Float64, 2}: The cached matrix value.
• vjp_inp::Array{Float64, 2}: VJP with respect to inputs.
• vjp_θ::Vector{Float64}: VJP with respect to parameters.
• glacier_id::Int64: Glacier ID in the list of simulated glaciers.

MatrixCacheInterp(nodes_H, nodes_∇S, interp_θ)

A mutable cache structure for storing interpolation data on a 2D grid, used to efficiently evaluate and reuse interpolated matrices and their gradients. This interpolation makes complex inversions feasible since it allows the precomputation of all gradients before the solving the reverse PDE associated to the adjoint variable.

Fields

• value::Array{Float64, 2}: Matrix to store value of the diffusivity.
• nodes_H::Vector{Float64}: Grid nodes corresponding to the first interpolation dimension, typically representing values of ice thickness H.
• nodes_∇S::Vector{Float64}: Grid nodes corresponding to the second interpolation dimension, typically representing absolute values of slope ∇S.
• interp_θ::Interpolations.GriddedInterpolation{Vector{Float64}, 2, Matrix{Vector{Float64}}, Interpolations.Gridded{InterPolations.Linear{InterPolations.Throw{OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}: A gridded linear interpolation object mapping (nodes_H, nodes_∇S) to parameter vectors θ. Used to compute interpolated parameters and enable fast evaluation during repeated model calls.

MultiLoss(; losses = (L2Sum(),), λs = (1.0,))

Combines multiple loss functions into a single weighted objective.

MultiLoss enables composing several individual loss terms—each possibly representing a different physical constraint, data fidelity term, or regularization penalty—into a single differentiable loss function.

Keyword Arguments (Constructor)

• losses::Tuple = (L2Sum(),): A tuple of loss objects (each subtype of AbstractLoss) to be combined.
• λs::Tuple = (1.0,): A tuple of scalar weights or hyperparameters corresponding to each loss term.

Fields (Struct)

• losses::TL: Tuple of loss functions.
• λs::TS: Tuple of weighting coefficients.

NeuralNetwork{
    ChainType <: Lux.Chain,
    ComponentVectorType <: ComponentVector,
    NamedTupleType <: NamedTuple,
} <: FunctionalModel

Feed-forward neural network.

Fields

• architecture::ChainType: Flux.Chain neural network architecture
• θ::ComponentVectorType: Neural network parameters
• st::NamedTupleType: Neural network status

No VJP flavor when the contribution of a given term should not be computed inside the adjoint calculation (e.g. MB).

NoVJP

PerGlacierModel <: TrainableModel

Abstract type representing per glacier optimizable components of the model. This is a subtype of TrainableModel. Typically used for classical inversions.

mutable struct Results{RES <: Sleipnir.Results, STAT <: TrainingStats}

Mutable struct containing the results of an inversion. This object stores both the results of the optimization through TrainingStats and the simulation results of the forward evaluations using the optimized variables through Sleipnir.Results. It expands the functionalities offered by Sleipnir.Results to allow saving the results of an inversion. Multiple dispatch is used to select either Sleipnir.Results or ODINN.Results.

Fields

• simulation::Vector{RES}: Vector of Sleipnir.Results representing the results of a forward evaluation for each glacier.

• stats::STAT: Training statistics gathered during the optimization.

  function Results( simulation::Vector{<: Sleipnir.Results}, stats::TrainingStats, )

Constructor for the Results object used to track statistics during training and the results of the forward evaluations simulated with the optimized variables.

Arguments

• simulation::Vector{<: Sleipnir.Results}: Vector of Sleipnir.Results associated to the forward simulation of each glacier.
• stats::TrainingStats: Training statistics.

RheologyRegularization(; reg = TikhonovRegularization())

Regularization of the gridded rheology A in the context of classical inversions. It can include a spatial smoothing operator through the field reg.

Keyword Arguments

• reg::AbstractSimpleRegularization = TikhonovRegularization(): Spatial regularization operator.

SIA2D_D_target(; interpolation=:None, n_interp_half=20,
                 prescale=nothing, postscale=nothing)

Inversion of general diffusivity as a function of physical parameters.

D(H, ∇S, θ) = H * NN(H, ∇S; θ)

So now we are learning the velocity field given by D * ∇S. This inversion is similar to learnign the velocity field assuming that this is parallel to the gradient in surface ∇S.

Arguments

• interpolation::Symbol = :None: Specifies the interpolation method. Options include :Linear, :None.
• n_interp_half::Int = 20: Half-width of the interpolation stencil. Determines resolution of interpolation.
• prescale::Union{Fin, Nothing} = nothing: Optional prescaling function or factor applied before parametrization. Must be of type Fin or nothing.
• postscale::Union{Fout, Nothing} = nothing: Optional postscaling function or factor applied after parametrization. Must be of type Fout or nothing.

Type Parameters

• Fin: Type of the prescale function or operator.
• Fout: Type of the postscale function or operator.

Supertype

• AbstractSIA2DTarget: Inherits from the abstract target type for 2D SIA modeling.

Returns

• An instance of SIA2D_D_target configured with optional scaling and interpolation parameters.

ScalarCacheGlacierId <: Cache

A cache structure for storing a scalar value as a zero-dimensional array of Float64 along with their associated vector-Jacobian products (VJP). It also stores the glacier ID as an integer. This is typically used to invert a single scalar per glacier. Fields:

• value::Array{Float64, 0}: The cached scalar value.
• vjp_inp::Array{Float64, 0}: VJP with respect to inputs. Must be defined but never used in practice since this cache is used for classical inversions and the law does not have inputs.
• vjp_θ::Vector{Float64}: VJP with respect to parameters.
• glacier_id::Int64: Glacier ID in the list of simulated glaciers.

Struct that defines the SciMLSensitivity adjoint flavor. This is the default behavior in ODINN.

SciMLSensitivityAdjoint

TikhonovRegularization(; operator = :laplacian, distance = 3)

A simple regularization type implementing Tikhonov regularization (also known as ridge regularization) for inverse problems.

This struct includes both the forward and reverse (adjoint) operators, which are required for the computation of the gradients with respect to the model parameters.

Keyword Arguments (Constructor)

• operator::Symbol = :laplacian: The regularization operator to use. Currently, only :laplacian is implemented, which penalizes large gradients by applying the Laplacian operator.
• distance::Integer = 3: A width parameter to determine how far from the margin evaluate the loss.

Fields (Struct)

• operator_forward::Function: The forward regularization operator (e.g., ∇²).
• operator_reverse::Function: The reverse-mode (VJP) of the operator (e.g., VJP_λ_∂∇²a_∂a).
• distance::Integer: The distance parameter controlling the extent of regularization.

mutable struct TrainingStats{
    F <: AbstractFloat,
    I <: Integer,
    RETC <: Union{String, Nothing},
    THETA_HIST <: ComponentVector
}

An object with the information of the training.

Fields

• retcode::RETC: Report code of the optimization.
• losses::Vector{F}: Vector storing the value of the loss function at each iteration.
• niter::I: Total number of iterations/epochs.
• θ::Union{<: ComponentVector, Nothing}: Parameters of neural network after training.
• θ_hist::Vector{THETA_HIST}: History of parameters of neural network during training.
• ∇θ_hist::Vector{THETA_HIST}: History of gradients training.
• initial_conditions::Union{<: Dict, Nothing}: Initial conditions used to solve the PDEs.
• lastCall::DateTime: Last time the callback diagnosis was called. This is used to compute the time per iteration.

TrainingStats(;
    retcode::Union{String, Nothing} = nothing,
    losses::Vector{F} = Float64[],
    niter::I = 0,
    θ::Union{ComponentVector, Nothing} = nothing,
    θ_hist::Union{Vector{ComponentVector}, Nothing} = ComponentVector[],
    ∇θ_hist::Union{Vector{ComponentVector}, Nothing} = ComponentVector[],
    initial_conditions::Union{Dict, Nothing} = nothing
) where {F <: AbstractFloat, I <: Integer}

Constructor for TrainingStats object used to store important information during training.

Arguments

• retcode: Report code of the optimization.
• losses: Vector storing the value of the loss function at each iteration.
• niter: Total number of iterations/epochs.
• θ: Parameters of neural network after training.
• θ_hist: History of parameters of neural network during training.
• ∇θ_hist: History of gradients training.
• initial_conditions: Initial conditions used to solve the PDEs.
• lastCall: Last time the callback diagnosis was called.

A mutable struct that holds parameters for a UDE (Universal Differential Equation).

UDEparameters{ADJ <: AbstractAdjointMethod} <: AbstractParameters

Fields

• sensealg::SciMLBase.AbstractAdjointSensitivityAlgorithm: The sensitivity algorithm used for adjoint sensitivity analysis.
• optimization_method::String: The optimization method to be used.
• target::Symbol: The target variable for the optimization.

UDEparameters(; sensealg, optim_autoAD, grad, optimization_method, empirical_loss_function, target) where {ADJ <: AbstractAdjointMethod}

Create a UDEparameters object for configuring the sensitivity analysis and optimization of a Universal Differential Equation (UDE).

Keyword Arguments

• sensealg::SciMLBase.AbstractAdjointSensitivityAlgorithm: The sensitivity algorithm to use for adjoint calculations. Defaults to InterpolatingAdjoint(autojacvec=SciMLSensitivity.EnzymeVJP()).
• optim_autoAD::AbstractADType: The automatic differentiation type for optimization. Defaults to Optimization.AutoEnzyme().
• grad::ADJ: The adjoint gradient computation method. Defaults to ContinuousAdjoint().
• optimization_method::String: The optimization method to use. Must be either "AD+AD" or "AD+Diff". Defaults to "AD+AD".
• empirical_loss_function::AbstractLoss: The loss function to use for optimization. Defaults to LossH().
• target::Union{Symbol, Nothing}: The target variable for optimization. Defaults to :A.

Returns

• A UDEparameters object configured with the specified sensitivity, optimization, and loss settings.

Description

This function creates a UDEparameters object that encapsulates the configuration for sensitivity analysis, optimization, and loss computation in a Universal Differential Equation (UDE) framework. It verifies that the provided optimization_method is valid and constructs the solver parameters accordingly.

Notes

• The optimization_method must be either "AD+AD" (automatic differentiation for both forward and backward passes) or "AD+Diff" (automatic differentiation combined with finite differences).
• The empirical_loss_function determines how the loss is computed during optimization.

struct VJPsPrepLaw <: AbstractPrepVJP

A container struct that holds all objects needed to compute vector-Jacobian products (VJPs) for a law using DifferentiationInterface.

Fields:

• f_θ_first: Function to evaluate the law with parameters θ as the first argument.
• f_inp_first: Function to evaluate the law with inputs as the first argument.
• prep_θ: Precomputed gradient preparation for parameters θ.
• prep_inp: Precomputed gradient preparation for inputs.

This struct is used to prepare the VJP computation with DifferentiationInterface (DI). Depending on the AD backend, DI might require to precompile code and this struct stores the results. This allows each VJP call to be fast in the adjoint PDE by reusing the preparation results.

VelocityRegularization(; reg = TikhonovRegularization(), components = :abs, distance = 3)

Regularization for velocity fields, combining a spatial smoothing operator with optional component control.

Keyword Arguments

• reg::AbstractSimpleRegularization = TikhonovRegularization(): Spatial regularization operator.
• components::Symbol = :abs: Determines which velocity components to regularize (e.g. :abs, :x, :y).
• distance::Integer = 3: Distance to glacier margin.

Base.copyto!(
    dest::InversionBinder{FI, CA},
    src::InversionBinder{FI, CA},
) where {FI <: Inversion, CA <: ComponentArray}

Overload Base.copyto! as we need a way to copy the SciMLStructure. It is used in SciMLSensitivity to differentiate the callbacks.

Base.zero(
    p::InversionBinder{FI, CA},
) where {FI <: Inversion, CA <: ComponentArray}

Overload Base.zero as we need a way to copy the SciMLStructure. It is used in SciMLSensitivity to differentiate the callbacks.

precompute_all_VJPs_laws!(
    SIA2D_model::SIA2Dmodel,
    SIA2D_cache::SIA2DCache,
    simulation::Inversion,
    glacier_idx::Integer,
    t::Real,
    θ,
)

Precomputes the vector-Jacobian products (VJPs) for all laws used in the SIA2D ice flow model for a given glacier, time, and model parameters.

Depending on which target (U, Y, or neither) is provided in SIA2D_model, this function checks if the corresponding law supports VJP precomputation and, if so, triggers the appropriate precompute routine for that law. If neither U nor Y is provided, precomputes VJPs for the A, C, and n laws.

Arguments

• SIA2D_model::SIA2Dmodel: The model containing the configuration and laws used for SIA2D ice flow.
• SIA2D_cache::SIA2DCache: A cache object holding intermediate values and storage relevant for precomputations.
• simulation::Inversion: Simulation object containing global simulation parameters.
• glacier_idx::Integer: Index of the glacier being simulated.
• t::Real: Current time in the simulation.
• θ: Model parameters or state variables for the simulation step.

Notes

• This routine is intended as a preparatory step for manual adjoint.
• Only laws supporting VJP precomputation are processed.

run!(simulation::Inversion)

Run the training process for a given Inversion simulation.

Arguments

• simulation::Inversion: The simulation object containing the parameters and settings for the inversion process.

Description

This function initiates the training of a Universal Differential Equation (UDE) for the provided simulation. It prints a message indicating the start of the training process, calls the train_UDE! function to perform the training, and collects the results in results_list. The results are intended to be saved using Sleipnir.save_results_file!, but this step is currently commented out and will be enabled once the optimization is working. Finally, the garbage collector is triggered to free up memory.

Notes

• The Sleipnir.save_results_file! function call is currently commented out and should be enabled once the optimization process is confirmed to be working.
• The garbage collector is explicitly run using GC.gc() to manage memory usage.

CallbackOptimizationSet(θ, l; callbacks)

Helper to combine callbacks for Optimization function. This executes the action of each callback. (equivalent to CallbackSet for DifferentialEquations.jl)

ComponentVector2Vector(cv::ComponentVector)

Transform a ComponentVector into a Vector of same length. This function creates a new Vector and does not mutate the original ComponentVector.

Arguments:

• cv::ComponentVector: Input ComponentVector.

Diffusivity(target::SIA2D_D_target; H, ∇S, θ, iceflow_model, glacier, params)

Compute the effective diffusivity field for a 2D shallow ice model using the diagnostic target and a predicted velocity matrix U.

This function uses a learned or specified model to estimate the velocity matrix U, then calculates the diffusivity as either H .* U (if dimensions match) or the averaged H times U if dimensions differ by one grid cell (staggered grid). Errors if dimensions are incompatible.

Arguments

• target::SIA2D_D_target: Diagnostic target object defining interpolation and scaling rules.

Keyword Arguments

• H: Ice thickness.
• ∇S: Ice surface slope.
• θ: Parameters of the model.
• iceflow_model: Iceflow model used for simulation.
• glacier: Glacier data.
• params: Model parameters.

Returns

• A matrix of diffusivity values with the same shape as H or staggered by one cell, depending on U.

Throws

• An error if the dimensions of U and H are not compatible for diffusivity calculation.

Notes

Supports both grid-matched and staggered configurations by averaging H where necessary.

Gauss Quadratrue for numerical integration

function LawA(
    nn_model::NeuralNetwork,
    params::Sleipnir.Parameters;
    precompute_VJPs::Bool = true,
)

Constructs a law object for the creep coefficient A in the SIA based on a neural network that takes as input the long term air temperature. The creep coefficient A with this law is a scalar. See also SIA2D_A_target.

Arguments

• nn_model::NeuralNetwork: A neural network model containing the architecture archi and state st used for evaluation of the law.
• params::Sleipnir.Parameters: Parameters struct used to retrieve the minimum and maximum values of A for scaling of the neural network output.
• precompute_VJPs::Bool: If true, enables precomputation of vector-Jacobian products before solving the adjoint PDE for efficient autodiff.

Returns

• A_law: A Law{ScalarCache} instance that computes the creep coefficient A based on an input temperature using the neural network. The law scales the output to the physical bounds defined by params.

Notes

• The VJP is computed automatically using DifferentiationInterface.

Details

• The function wraps the architecture and state of the neural network in aStatefulLuxLayer.
• The resulting law takes input variables, applies the neural network, and scales its output to be between params.physical.minA and params.physical.maxA.
• The in-place assignment to cache is ignored in differentiation to allow gradient computation with Zygote when using DifferentiationInterface.
• The init_cache function initializes the cache with a scalar zero for the forward placeholder, and with a vector of zeros for the VJP placeholder.

Example

nn_model = NeuralNetwork(params)
A_law = LawA(nn_model, params; precompute_VJPs = false)

LawA(params::Sleipnir.Parameters; scalar::Bool=true)

Construct a law that defines an ice rheology A per glacier to invert. This can be either a spatially varying A or a scalar value per glacier based on the value of scalar.

Arguments

• params::Sleipnir.Parameters: Parameters struct used to retrieve the minimum and maximum values of A for scaling the parameter to invert.
• scalar::Bool: Whether the ice rheology to invert is a scalar per glacier, or a spatially varying A per glacier (matrix to invert).

LawC(params::Sleipnir.Parameters; scalar::Bool=true)

Construct a law that defines a sliding coefficient C per glacier to invert. This can be either a spatially varying C or a scalar value per glacier based on the value of scalar.

The tanh-based parameterisation ensures C stays within [0, params.physical.maxC] and is differentiable through θ.C at every RHS call (required by SciMLSensitivityAdjoint).

Arguments

• params::Sleipnir.Parameters: Parameters struct used to retrieve the maximum value of C for scaling the parameter to invert.
• scalar::Bool: Whether the sliding coefficient to invert is a scalar per glacier, or a spatially varying C per glacier (matrix to invert).

LawU(
    nn_model::NeuralNetwork,
    params::Sleipnir.Parameters;
    max_NN::Union{F, Nothing} = 50.0,
    prescale_bounds::Union{Vector{Tuple{F,F}}, Nothing} = [(0.0, 300.0), (0.0, 0.5)],
) where {F <: AbstractFloat}

Constructs a law object for the diffusive velocity U in the SIA based on a neural network that takes as input the ice thickness H̄ and the surface slope ∇S. The diffusive velocity U with this law is a matrix and the diffusivity in the SIA is obtained through D = U * H̄. See also SIA2D_D_target.

Arguments

• nn_model::NeuralNetwork: A neural network model containing the architecture archi and state st used for evaluation of the law.
• params::Sleipnir.Parameters: Parameters struct. Not used for the moment but kept as an argument to keep consistency with other equivalent functions LawA and LawY.
• max_NN::Union{F, Nothing}: Expected maximum value of the neural network output. If set to nothing, no postscaling is applied.
• prescale_bounds::Union{Vector{Tuple{F,F}}, Nothing}: Vector of tuples where each tuple defines the lower and upper bounds of the input for scaling. If set to nothing, no prescaling is applied.
• precompute_interpolation::Bool: Determines which cache to use depending if interpolation is used or not for the evaluation of gradients.
• precompute_VJPs::Bool: Determines is VJPs are stored in the cache during the reverse step.

Returns

• U_law: A Law{Array{Float64, 2}} instance that computes the diffusive velocity U based on the ice thickness H̄ and the surface slope ∇S using the neural network. The law scales the output using the max_NN argument.

Notes

• The computation is compatible with Zygote for automatic differentiation.

Details

• The function wraps the architecture and state of the neural network in aStatefulLuxLayer.
• The resulting law takes input variables, applies the neural network, and scales its output to match max_NN.
• The in-place assignment to cache is ignored in differentiation to allow gradient computation with Zygote.
• The init_cache function initializes the cache with a zero matrix.

Example

nn_model = NeuralNetwork(params)
bounds_H = (0.0, 300.0)
bounds_∇S = (0.0, 0.5)
U_law = LawU(nn_model, params; max_NN = 50.0, prescale_bounds = [bounds_H, bounds_∇S])

LawY(
    nn_model::NeuralNetwork,
    params::Sleipnir.Parameters;
    max_NN::Union{F, Nothing} = nothing,
    prescale_bounds::Vector{Tuple{F,F}} = [(-25.0, 0.0), (0.0, 500.0)],
) where {F <: AbstractFloat}

Constructs a law object for the hybrid diffusivity Y in the SIA based on a neural network that takes as input the long term air temperature and the ice thickness H̄. The hybrid diffusivity Y with this law is a matrix as it depends on the ice thickness. This law is used in an hybrid setting where the n exponent in the mathematical expression of the diffusivity is different from the one used to generate the ground truth. The goal of this law is to retrieve the missing part of the diffusivity. Please refer to SIA2D_D_hybrid_target for a mathematical definition.

Arguments

• nn_model::NeuralNetwork: A neural network model containing the architecture archi and state st used for evaluation of the law.
• params::Sleipnir.Parameters: Parameters struct used to retrieve the maximum value of A for scaling of the neural network output.
• max_NN::Union{F, Nothing}: Expected maximum value of the neural network output. If not specified, the law takes as an expected maximum value params.physical.maxA.
• prescale_bounds::Vector{Tuple{F,F}}: Vector of tuples where each tuple defines the lower and upper bounds of the input for scaling.

Returns

• Y_law: A Law{Array{Float64, 2}} instance that computes the hybrid diffusivity Y based on an input temperature and ice thickness using the neural network. The law scales the output to the physical bounds defined by params.

Notes

• The computation is compatible with Zygote for automatic differentiation.

Details

• The function wraps the architecture and state of the neural network in aStatefulLuxLayer.
• The resulting law takes input variables, applies the neural network, and scales its output to match the maximum value which is either max_NN or params.physical.maxA.
• The in-place assignment to cache is ignored in differentiation to allow gradient computation with Zygote.
• The init_cache function initializes the cache with a zero matrix.

Example

nn_model = NeuralNetwork(params)
bounds_T = (-25.0, 0.0)
bounds_H = (0.0, 500.0)
Y_law = LawY(nn_model, params; prescale_bounds = [bounds_T, bounds_H])

This function allows to extend the Wrapper layers define in Lux to matrices operations.

Constructor for the Parameters type. Since some of the subtypes of parameters are defined in different packages of the ODINN ecosystem, this constructor will call the constructors of the different subtypes and return a Parameters object with the corresponding subtypes. The Parameters mutable struct is defined in Sleipnir.jl using abstract types, which are later on defined in the different packages of the ODINN ecosystem.

Parameters(;
        physical::PhysicalParameters = PhysicalParameters(),
        simulation::SimulationParameters = SimulationParameters(),
        solver::SolverParameters = SolverParameters(),
        hyper::Hyperparameters = Hyperparameters(),
        UDE::UDEparameters = UDEparameters()
        )

Keyword arguments

• physical::PhysicalParameters: Physical parameters for the simulation.
• simulation::SimulationParameters: Parameters related to the simulation setup.
• solver::SolverParameters: Parameters for the solver configuration.
• hyper::Hyperparameters: Hyperparameters for the model.
• UDE::UDEparameters: Parameters specific to the UDE (Universal Differential Equation).

Currently just used for Enzyme

Inverse with batch

Compute gradient glacier per glacier

T_A_Alaw(simulation::Simulation, glacier_idx::Integer, θ, t::AbstractFloat)

Evaluate the A law when it defines a mapping between the long term air temperature and the creep coefficient A and return both the input temperature T and the computed creep coefficient A.

Arguments

• simulation::Simulation: The simulation object containing model data and parameters.
• glacier_idx::Integer: Index specifying which glacier to evaluate.
• θ: Model parameters to be used in the law.
• t::AbstractFloat: The time at which to evaluate the law. For this law it is useless but in the general setting, a law needs a time t in order to retrieve the inputs. For the sake of consistency, this input was kept.

Returns

• (T, A): A tuple containing:

  • T: The input long term air temperature for the specified glacier.
  • A: The evaluated creep coefficient for the specified glacier.

Details

• The function checks that the inputs of the A law are exactly as expected (long term air temperature only).
• Retrieves the long term air temperature for the specific glacier.
• Evaluates the creep coefficient using the law.
• Returns both the temperature and creep coefficient as a tuple. Since the cache of A is a zero dimensional array, it is converted to float before returning the value.

Example

T, A = T_A_Alaw(simulation, glacier_idx, θ, 2010.0)

VJP_λ_∂SIA∂H_continuous(
    λ::Matrix{R},
    H::Matrix{R},
    θ,
    simulation::SIM,
    t::R,
) where {R <: Real, SIM <: Simulation}

Implementation of the continuous VJP of the SIA2D equation with respect to H. Given λ and H, it returns the VJP of λ^T * ∂(SIA2D)/∂H (H).

Arguments:

• λ::Matrix{R}: Adjoint state, also called output gradient in reverse-mode AD.
• H::Matrix{R}: Ice thickness which corresponds to the input state of the SIA2D.
• simulation::SIM: Simulation parameters.
• t::R: Time value, not used as SIA2D is time independent.

Returns:

• dλ::Matrix{R}: Jacobian vector product, also called input gradient in reverse-mode AD.

VJP_λ_∂SIA∂H_discrete(
    λ::Matrix{R},
    H::Matrix{R},
    θ,
    simulation::SIM,
    t::R,
) where {R <:Real, SIM <: Simulation}

Implementation of the discrete VJP of the SIA2D equation with respect to H. Given λ and H, it returns the VJP of λ^T * ∂(SIA2D)/∂H (H).

Arguments:

• λ::Matrix{R}: Adjoint state, also called output gradient in reverse-mode AD.
• H::Matrix{R}: Ice thickness which corresponds to the input state of the SIA2D.
• simulation::SIM: Simulation parameters.
• t::R: Time value, not used as SIA2D is time independent.

Returns:

• dλ::Matrix{R}: Jacobian vector product, also called input gradient in reverse-mode AD.

VJP_λ_∂SIA∂θ_continuous(
    λ::Matrix{R},
    H::Matrix{R},
    θ,
    simulation::SIM,
    t::R,
) where {R <: Real, SIM <: Simulation}

Implementation of the continuous VJP of the SIA2D equation with respect to θ. Given λ, H and θ, it returns the VJP of λ^T * ∂(SIA2D)/∂θ (θ).

Arguments:

• θ: Vector of parameters
• λ::Matrix{R}: Adjoint state, also called output gradient in reverse-mode AD.
• H::Matrix{R}: Ice thickness which corresponds to the input state of the SIA2D.
• simulation::SIM: Simulation parameters.
• t::R: Time value, not used as SIA2D is time independent.

Returns:

• ∂θ: Jacobian vector product with respect to θ, also called input gradient in reverse-mode AD. It has the same type as θ.

VJP_λ_∂SIA∂θ_discrete(
    λ::Matrix{R},
    H::Matrix{R},
    θ,
    simulation::SIM,
    t::R,
) where {R <: Real, SIM <: Simulation}

Implementation of the discrete VJP of the SIA2D equation with respect to θ. Given λ, H and θ, it returns the VJP of λ^T * ∂(SIA2D)/∂θ (θ).

Arguments:

• θ: Vector of parameters
• λ::Matrix{R}: Adjoint state, also called output gradient in reverse-mode AD.
• H::Matrix{R}: Ice thickness which corresponds to the input state of the SIA2D.
• simulation::SIM: Simulation parameters.
• t::R: Time value, not used as SIA2D is time independent.

Returns:

• ∂θ: Jacobian vector product with respect to θ, also called input gradient in reverse-mode AD. It has the same type as θ.

VJP_λ_∂∇²a_∂a(λ::Matrix{R}, a::Matrix{R}, Δx::R, Δy::R) where {R<:Real}

Computes the vector-Jacobian product (VJP) of the Laplacian operator ∇² with respect to its input field a. This function effectively propagates sensitivities (adjoints) λ backward through the Laplacian, as required in adjoint or reverse-mode differentiation.

Arguments

• λ::Matrix{R}: Adjoint field associated with the Laplacian output.
• a::Matrix{R}: Input scalar field to the Laplacian operator.
• Δx::R: Grid spacing in the x-direction.
• Δy::R: Grid spacing in the y-direction.

Returns

• Matrix{R}: The adjoint (VJP) with respect to a, i.e. ∂⟨λ, ∇²a⟩/∂a.

Vector2ComponentVector(v::Vector, cv_template::ComponentVector)

Transform a vector v to a ComponentVector that has the same structure as cv_template. This function creates a new ComponentVector and copies the values of v explicitly. The arguments v and cv_template must be of the same length.

Arguments:

• v::Vector: Vector whose values are copied.
• cv_template::ComponentVector: ComponentVector whose structure is used to create a new ComponentVector.

Function to evaluate derivatives of ice surface velocity in D inversion.

TODO: This functions right now just make a call to the regular functions used for the calculation of the adjoint. This is not correct, but we keep it as this for now until we figure out how to do this in the case of the D inversion.

_batch_iceflow_UDE(
    container::InversionBinder,
    glacier_idx::Integer,
    iceflow_prob::ODEProblem;
    processVelocity::Union{Nothing, Function} = nothing
)

Define the callbacks to be called by the ODE solver, solve the ODE and create the results.

_ml_model_postscale(
    Y::Vector,
    max_NN,
)

Applies an exponential transformation to each element in Y, then rescales the result by multiplying with max_NN. For each element, the transformation is: max_NN * exp((Y - 1.0) / Y)

Arguments

• Y::Vector: Values to be post-processed.
• max_NN: Scalar representing the maximum value for rescaling.

Returns

• The rescaled values after applying the exponential transformation.

_ml_model_prescale(
    X::Vector,
    prescale_bounds::Vector{Tuple{F, F}},
) where {F <: AbstractFloat}

Scales each element of the input vector X using the corresponding bounds from prescale_bounds. For each index i, X[i] is normalized based on the interval specified in prescale_bounds[i] using the normalize function. This function is typically used to ensure that the scales of the inputs of a neural network are comparable to each other.

Arguments

• X::Vector: A vector of input values to be normalized.
• prescale_bounds::Vector{Tuple{F, F}}: A vector of tuples specifying the lower and upper bounds for normalization of each corresponding element in X.

Returns

• A vector where each element is the normalized value of the corresponding input, using the specified bounds.

Notes

• The length of X and prescale_bounds must be equal.

_pred_NN(inp::Vector{F}, smodel, θ, prescale, postscale) where {F <: AbstractFloat}

Compute the output of a neural network model on the input vector inp.

Arguments

• inp::Vector{F}: Input vector of floats.
• smodel: The neural network model.
• θ: Parameters for the neural network model.
• prescale: Function to scale the input vector before passing it to the model.
• postscale: Function to scale the model output.

Returns

• The single (scalar) output value from the neural network after applying prescale to the input, evaluating the model, and then applying postscale. The result is extracted via only.

Notes

• The function assumes that the neural network, when evaluated, returns an iterable with exactly one element.
• Using only will throw an error if the output is not exactly one element.

Example

mymodel = StatefulLuxLayer{true}(archi, nothing, st)
y = _pred_NN([1.0, 2.0], mymodel, θ, prescale_fn, postscale_fn)

aggregate∇θ(∇θ::Vector{<: ComponentArray}, θ, submodels::TrainableComponents)

Aggregate the vector of gradients ∇θ as a single ComponentArray. The argument ∇θ is the vector of all the gradients computed for each glacier. This function aggregates them based on the optimizable components of submodels.

backward_loss(
    lossType::MultiLoss,
    H_pred::Matrix{F},
    H_ref,
    V_ref, Vx_ref, Vy_ref,
    t::F,
    glacier_idx::Integer,
    θ,
    simulation,
    normalization::F,
    Δt,
) where {F <: AbstractFloat}

Computes the gradient of a composite loss defined by a MultiLoss object with respect to both the predicted field H_pred and model parameters θ.

Each sub-loss's backward gradient is weighted by its corresponding coefficient in lossType.λs and summed to form the total gradient.

Arguments

• lossType::MultiLoss: Composite loss object containing individual losses and weights.
• H_pred::Matrix{F}: Predicted ice thickness.
• H_ref::Matrix{F}: Reference ice thickness.
• t::F: Current time or simulation step.
• glacier_idx::Integer: Glacier id in the list of glaciers in simulation.
• θ: Model parameters used in the simulation.
• simulation: Simulation object providing necessary context for gradient computation.
• normalization::F: Normalization factor applied within each individual loss.
• Δt: Named tuple containing the time step to use for the approximation of continuous in time loss terms. For example if LossH is used, there must be a term Δt.H containing the time step since the last computation of the ice thickness loss term. If the current time t where the loss is evaluated does not correspond to a time step of the LossH term, then the value of Δt.H has no impact.

Returns

• (∂L∂H, ∂L∂θ): Tuple containing:

  • ∂L∂H::Matrix{F}: Gradient of the composite loss with respect to H_pred.
  • ∂L∂θ: Gradient of the composite loss with respect to model parameters θ.

batch_loss_iceflow_transient(
    container::InversionBinder,
    glacier_idx::Integer,
    iceflow_prob::ODEProblem,
)

Solve the ODE, retrieve the results and compute the loss.

Arguments:

• container::InversionBinder: SciMLStruct that contains the simulation structure and the vector of parameters to optimize.
• glacier_idx::Integer: Index of the glacier.
• iceflow_prob::ODEProblem: Iceflow problem defined as an ODE with respect to time.

callback_diagnosis(θ, l, simulation; save::Bool = false, tbLogger::Union{<: TBLogger, Nothing} = nothing)

Callback function to track and diagose training. It includes print and updates in simulation::Simulation. It also logs training statistics with tensorboard if tbLogger is provided.

callback_plots_A(θ, l, simulation)

Callback function to generate plots during training.

Normalization of D to cap at a maximum physical value

compute_numerical_gradient(
    x,
    args,
    fct::Function,
    ϵ::F;
    varStr::String = ""
) where {T, F <: AbstractFloat}

Compute the gradient of a function by using numerical differences. The function is evaluated N+1 times with N the length of x. The function has two arguments: the variable with respect to which the gradient is computed and extra arguments which are not differentiated.

Arguments:

• x: Variable x where to evaluate the gradient.
• args: Extra arguments.
• fct::Function: Function to differentiate.
• epsilon::F: Size of perturbation to use in the numerical differences.
• varStr::String: Variable to print in the progress bar. Example: "of A" will print "Computing gradient of A using finite differences with..."

Returns:

• grad: Numerical gradient.

create_interpolation(A::Matrix; n_interp_half::Int) -> Vector{Float64}

Construct a one-dimensional interpolation grid from the elements of a matrix A by flattening it and delegating to create_interpolation(::Vector). This is a convenience method that allows users to pass a 2D array to the function create_interpolation(A::Vector) directly without manually reshaping it.

create_interpolation(
    A::Vector;
    n_interp_half::Int,
    dilation_factor = 1.0,
    minA_unif::Union{F, Nothing} = nothing,
    minA_quantile::Union{F, Nothing} = nothing,
    maxA_unif::Union{F, Nothing} = nothing,
    maxA_quantile::Union{F, Nothing} = nothing
    ) where {F <: AbstractFloat}

Construct a one-dimensional interpolation grid from the data in A, combining uniformly spaced and quantile-based sampling points.

This hybrid interpolation grid provides both coverage of the entire range of values and higher resolution in regions where A has dense data, making it useful for interpolation or machine learning applications that need balanced sampling.

Arguments

• A::Vector: Input data vector (typically containing positive values).
• n_interp_half::Int: Number of points used for both the uniform and quantile-based subsets of the interpolation grid.
• dilation_factor::Real = 1.0: Optional multiplier applied to maximum(A) to slightly extend the grid beyond the data range (useful to avoid extrapolation issues).
• minA_unif::Union{F, Nothing} = nothing: Minimum value used for the uniform interpolation
• minA_quantile::Union{F, Nothing} = nothing: Maximum value used for the uniform interpolation
• maxA_unif::Union{F, Nothing} = nothing: Minimum value used for the quantile interpolation
• maxA_quantile::Union{F, Nothing} = nothing: Maximum value used for the quantile interpolation

Returns

A sorted, unique vector of interpolation nodes combining:

• n_interp_half uniformly spaced values between 0 and dilation_factor * maximum(A)
• n_interp_half quantile-based values computed from the positive entries of A

create_results(θ, simulation::Inversion, mappingFct; processVelocity::Union{Nothing, Function} = nothing)

Given the parameters θ, solve the iceflow problem for all the glaciers and aggregate the results for all of them. This function is typically used at the end of a training once θ has been optimized and one wants to run one last forward simulation in order to retrieve statistics about each of the iceflow problems.

Arguments:

• θ: Parameters to use for the forward simulation.
• simulation::Inversion: Simulation structure that contains all the required information about the inversion.
• mappingFct: Function to use to process the glaciers. Either map for a sequential processing or pmap for multiprocessing.
• processVelocity::Union{Nothing, Function}: Post processing function to map the ice thickness to the surface velocity. It is called before creating the results. It takes as inputs simulation, ice thickness (matrix) and the associated time and returns 3 variables Vx, Vy, V which are all matrix. Defaults is nothing which means no post processing is applied.

define_iceflow_prob(
    simulation::Inversion,
    glacier_idx::Integer,
)

Given a simulation struct and a glacier_idx, build the iceflow problem that has to be solved in the ODE solver. In practice, the returned iceflow problem is used inside simulate_iceflow_UDE! through remake. The definition of the iceflow problem has to be done outside of the gradient computation, otherwise Zygote fails at differentiating it.

enable_multiprocessing(params::Sleipnir.Parameters) -> Int

Configures and enables multiprocessing based on the provided simulation parameters.

Arguments

• params::Sleipnir.Parameters: A parameter object containing simulation settings, including the number of workers (params.simulation.workers) and whether multiprocessing is enabled (params.simulation.multiprocessing).

Behavior

• If multiprocessing is enabled (params.simulation.multiprocessing == true) and the specified number of workers (params.simulation.workers) is greater than 0:

  • Adds the required number of worker processes if the current number of processes (nprocs()) is less than the specified number of workers.
  • Suppresses precompilation output on the worker processes and ensures the ODINN module is loaded on all workers.
  • If the specified number of workers is 1, removes all worker processes.

Returns

• The number of worker processes (nworkers()) after configuration.

Notes

• This function uses @eval to dynamically add or remove worker processes.
• Precompilation output is suppressed on workers to reduce noise in the console.

eval_law(law::AbstractLaw, simulation::Simulation, glacier_idx::Integer, inputs::NamedTuple, θ)

Evaluates a law on the specified glacier within a simulation context and for a user defined input.

Arguments

• law::AbstractLaw: The law object to be evaluated. Must provide a function f and an init_cache method.
• simulation::Simulation: The simulation context, containing model parameters and machine learning components.
• glacier_idx::Integer: Index identifying which glacier in the simulation to evaluate the law for.
• input_values::NamedTuple: Input data required by the law and provided by the user.
• θ: Weights used in the law to make inference. This can be nothing when the law has no parameter.

Returns

• The updated cache after evaluating the law. The cache contains the result of the law's computation for the specified glacier and inputs.

Details

• The function initializes a cache for the law using init_cache.
• If the simulation has a machine learning model, the model's parameters (θ) are updated in-place with the provided θ.
• The law's function is then called with the cache, inputs, and parameters. The result is stored in the cache and the cache is returned.
• In future versions, the design may change so that only inputs and θ are needed, with the cache handled separately so that no simulation is required.

Example

result = eval_law(simulation.model.iceflow.A, simulation, glacier_idx, (; T = 273.15), θ)    # Initialize the cache to be able to make an inference of the law

evaluate_H₀(
    θ::ComponentArray,
    glacier::Glacier2D,
    filter::Symbol,
    glacier_id::Integer,
)

Evaluate the initial ice thickness H₀ for a given glacier, optionally applying a smooth thresholding function.

Arguments

• θ::ComponentArray: A ComponentArray containing glacier parameters.

• glacier::Glacier2D: Glacier for which to evaluate H₀.

• filter::Symbol: Specifies the smoothing function to apply to the raw initial condition:

  • :identity: applies the identity function (no change).
  • :softplus: applies the softplus function log(1 + exp(x)) to ensure positivity.
  • :Zang1980: applies the σ_zang function (Zang 1980) as a smooth positivity threshold.

• glacier_id::Integer: Index of the glacier in order to retrieve the parameters of the IC in θ.

Returns

• A numeric value or array representing the filtered initial ice thickness for the specified glacier.

evaluate_∂H₀(
    θ::ComponentArray,
    glacier::Glacier2D,
    filter::Symbol,
    glacier_id::Integer,
)

Evaluate the derivative of the initial ice thickness H₀ for a given glacier, optionally applying a smooth thresholding function.

Arguments

• θ::ComponentArray: A ComponentArray containing glacier parameters.

• glacier::Glacier2D: Glacier for which to evaluate ∂H₀.

• filter::Symbol: Specifies the smoothing function to apply to the raw initial condition:

  • :identity: applies the identity function (no change).
  • :softplus — applies the softplus function log(1 + exp(x)) to ensure positivity.
  • :Zang1980 — applies the σ_zang function (Zang 1980) as a smooth positivity threshold.

• glacier_id::Integer: Index of the glacier in order to retrieve the parameters of the IC in θ.

Returns

• A numeric value or array representing the filtered initial ice thickness for the specified glacier.

feed_input_cache!(
    SIA2D_model::SIA2Dmodel,
    SIA2D_cache::SIA2DCache,
    simulation,
    glacier_idx::Integer,
    θ,
    result
)

Populate the input cache of an SIA2DCache instance with interpolation nodes for ice thickness (H) and surface slope (∇S), based on the results of a previous forward simulation. This function is required just when results of the forward pass are required to evaluate the elements of the cache in the reverse step.

This function prepares the interpolation knots used later by the reverse evaluation of the adjoint SIA2D model.

Right now, this function is just required for the inversion w.r.t to D, which is indicated by the boolean variable SIA2D_model.U_is_provided. Other inversions may not required the definition of this function.

Arguments

• SIA2D_model::SIA2Dmodel: The 2D shallow-ice approximation model instance.
• SIA2D_cache::SIA2DCache: The cache object that stores precomputed interpolation nodes.
• simulation: The simulation object containing glacier configurations and model settings.
• glacier_idx::Integer: Index of the glacier within simulation.glaciers for which the cache is being populated.
• θ: Model parameters (not directly used in this function but included for interface consistency).
• result: Output of a previous forward run containing ice thickness fields H.

fourier_feature(v, n::Integer=10, random=false, σ=5.0)

Generates a Fourier feature embedding of a vector v, optionally using randomized projections.

Arguments

• v: Input vector to be transformed (typically a coordinate or feature vector).
• n::Integer=10: Number of Fourier features to generate (default is 10).
• random::Bool=false: Whether to use random Fourier features (default is false).
• σ::Float64=5.0: Standard deviation of the normal distribution used for random feature projection (only used if random=true).

Returns

• A 2n-dimensional vector consisting of sine and cosine features of the transformed input vector.

Notes

Fourier features help to overcome spectral bias in neural networks and can further help to learn higher frequncy components of the function faster. For more information, see Tancik et. al (2020), "Fourier Features Let Networks Learn High Frequency Functions in Low Dimensional Domains".

Example

v = [0.5, 1.0]
features = fourier_feature(v, n = 4, random = true, σ = 2.0)

generate_batches(simulation::S; shuffle=false) where {S <: Simulation}

Generate a data loader for batching simulations, optionally shuffling the batches.

Arguments

• simulation::S: A Simulation object (or subtype of Simulation) containing the data to be batched.
• shuffle::Bool=false: A flag indicating whether to shuffle the batches. Defaults to false.

Returns

• A DataLoader object that provides batched access to the simulation data.

Description

This function creates a DataLoader for batching the provided simulation object. The DataLoader allows for efficient iteration over the simulation data in batches. The batch size is set to 1 by default, and the shuffle flag determines whether the batches are shuffled. If shuffle is enabled, a warning is logged to indicate that the batches used for parallelization are being shuffled.

Notes

• The batch size is fixed at 1 in this implementation. To modify the batch size, you may need to adjust the DataLoader initialization.
• Shuffling the batches may affect reproducibility and parallelization behavior.

generate_simulation_batches(simulation::Inversion)

Generate batches of simulations from a Inversion object for parallel or batched processing.

Arguments

• simulation::Inversion: A Inversion object containing the model, glaciers, parameters, results, and statistics for the simulation.

Returns

• A vector of Inversion objects, each representing a batch of simulations. Each batch contains a subset of glaciers, models, and results from the original simulation.

Description

This function splits the glaciers and associated data in the simulation object into smaller batches for processing. Each batch is represented as a new Inversion object. The number of batches is determined by the nbatches variable (currently set to 1). If the simulation results are empty, the function creates batches with empty results. Otherwise, it includes the corresponding results for each glacier in the batches.

Notes

• The number of glaciers (ninstances) must be divisible by the number of batches (nbatches). An assertion is used to enforce this condition.
• The function currently defaults to nbatches = 1, meaning no actual batching is performed. This can be updated to use simulation.parameters.hyper.batchsize for dynamic batching. # TODO: we need to change how this is done here, with a manual =1 batch size
• If the simulation results are empty, the function creates batches with empty results objects. #simulation.parameters.hyper.batchsize

get_NN(θ_trained, ft; lightNN=false)

Generates a neural network.

Arguments

• θ_trained: Pre-trained neural network parameters (optional).
• ft: Float type used.
• lightNN: Boolean that determines if a light architecture is returned or not.

Returns

• UA: Lux.Chain neural network architecture.
• θ: Neural network parameters.
• st: Lux state.

grad_loss_iceflow!(dθ, θ, simulation::Inversion, mappingFct)

Compute the gradient with respect to θ for all the glaciers and assign the result in-place to dθ.

Arguments:

• dθ: Gradient of the parameters where the computed gradient should be stored.
• θ: Parameters to differentiate.
• simulation::Inversion: Simulation structure that contains all the required information about the inversion.
• mappingFct: Function to use to process the glaciers. Either map for a sequential processing or pmap for multiprocessing.

grad_loss_iceflow!(θ, simulation::Inversion, mappingFct)

Compute the gradient with respect to θ for all the glaciers and return the result out-of-place. See the in-place implementation for more information.

grad_parallel_loss_iceflow!(θ, simulation::Inversion, glacier_idx::Integer)

Compute the gradient with respect to θ for a particular glacier and return the computed gradient. This function defines the iceflow problem and then calls Zygote to differentiate batch_loss_iceflow_transient with respect to θ. It uses the SciMLSensitivity implementation under the hood to compute the adjoint of the ODE.

function loss(
    lossType::LogSum,
    a::Matrix{F},
    b::Matrix{F},
    mask::BitMatrix,
    normalization::F,
) where {F <: AbstractFloat}

Compute logarithmic loss function for ice velocity fields following Morlighem, M. et al., "Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica". Geophys. Res. Lett. 37, (2010). Given a minimum velocity ϵ the absolute velocity given by a and b, it computes the sum of

log^2( (a + ϵ) / (b + ϵ) )

It has been shown that this loss function enables robust estimation of drag coefficient.

loss(
    lossType::MultiLoss,
    H_pred::Matrix{F},
    H_ref,
    V_ref, Vx_ref, Vy_ref,
    t::F,
    glacier_idx::Integer,
    θ,
    simulation,
    normalization::F,
    Δt,
) where {F <: AbstractFloat}

Computes the weighted composite loss for a prediction H_pred against a reference H_ref using a MultiLoss object.

Each individual loss in lossType.losses is evaluated and multiplied by its corresponding weight in lossType.λs. The final loss is the sum of these weighted contributions.

Arguments

• lossType::MultiLoss: Composite loss object containing individual losses and weights.
• H_pred::Matrix{F}: Predicted ice thickness.
• H_ref::Matrix{F}: Reference ice thickness.
• t::F: Current time or simulation step.
• glacier_idx::Integer: Glacier id in the list of glaciers in simulation.
• θ: Model parameters used in the simulation.
• simulation: Simulation object providing necessary context for loss evaluation.
• normalization::F: Normalization factor applied within each individual loss.
• Δt: Named tuple containing the time step to use for the approximation of continuous in time loss terms. For example if LossH is used, there must be a term Δt.H containing the time step since the last computation of the ice thickness loss term. If the current time t where the loss is evaluated does not correspond to a time step of the LossH term, then the value of Δt.H has no impact.

Returns

• F: The total scalar loss, computed as the sum of weighted individual losses.

loss_iceflow_transient(θ, simulation::Inversion, mappingFct)

Given the parameters θ, this function:

1. Solves the iceflow problem for all the glaciers.
2. Computes the loss function defined as the sum of the loss functions for each of the glaciers. The loss function of each glacier depends on the type of loss. Refer to empirical_loss_function in the UDE parameters for more information. The loss function is transient meaning that the state of the glacier is compared to a reference at different time steps over the simulated period.
3. Return the value of the loss function.

Arguments:

• θ: Parameters to use for the forward simulation.
• simulation::Inversion: Simulation structure that contains all the required information about the inversion.
• mappingFct: Function to use to process the glaciers. Either map for a sequential processing or pmap for multiprocessing.

merge_batches(results::Vector)

Merge simulation results from multiple batches into a single collection.

Arguments

• results::Vector: A vector where each element is a collection of results (e.g., arrays or vectors) from a batch.

Returns

• A single collection containing all the merged results from the input batches.

Description

This function takes a vector of results from multiple simulation batches and merges them into a single collection using vertical concatenation (vcat). It is useful for combining results that were processed in parallel or in separate batches.

normalize(X; lims::Tuple{F, F}; method = :shift) where {F <: AbstractFloat}

Normalize a variable by using an affine transformation defined by some input lower and upper bounds (m, M) and transforming to O(1) scale.

Arguments

• X: Input value.
• lims::Tuple{F, F}: Lower and upper bounds to use in the affine transformation.
• method::Symbol: Method to scale data.

Returns

• The input variable scaled by the affine transformation.

parallel_loss_iceflow_transient(θ, simulation::Inversion)

Loop over a list of glaciers to process. When multiprocessing is enabled, each call of this function has a dedicated process. This function calls batch_loss_iceflow_transient which returns both the loss and the result structure. The function keeps only the loss.

plot_law(law::AbstractLaw, simulation::Simulation, inputs::NamedTuple, θ; glacier_idx=1, idx_fixed_input=0)

Plot a law function with one or two input variables.

Arguments

• law::AbstractLaw: The law to plot (e.g., a sliding law or creep law).
• simulation::Simulation: The simulation containing glaciers and parameters.
• inputs::NamedTuple: Named tuple of input variables for the law (e.g., (T=iTemp(), H̄=iH̄())).
• θ: Parameters for the law (can be nothing for laws without parameters).

Keyword Arguments

• glacier_idx::Integer=1: Index of the glacier to use for extracting input values (for 2D inputs).
• idx_fixed_input::Integer=0: For two-input laws, index (1 or 2) of the input to fix at its mean value. If 0, plots a 3D surface.
• plot_full_input_range::Bool=false: If true and plotting a 1D law, uses the full physical range for that input.
• ground_truth_law::Union{AbstractLaw, Nothing}=nothing: Optional ground truth law to overlay on the plot for comparison.

Returns

• A plot figure (1D line plot, 2D scatter, or 3D surface) and saves it to the simulation's working directory.

Examples

# 1D plot (temperature input for Cuffey-Paterson law)
plot_law(A_law, simulation, (T = iTemp(),), nothing)

pretraining(architecture::Lux.Chain;
            X::Matrix,
            Y::Matrix,
            nepochs::Int=3000,
            lossfn::GenericLossFunction=MSLELoss(; agg=mean, epsilon=1e-10),
            rng::AbstractRNG=Random.default_rng())

Pretrains a neural network model using a input and output.

Arguments

• architecture::Lux.Chain: The neural network architecture to be trained.
• X::Matrix: Input feature matrix where each column is a feature vector.
• Y::Matrix: Target output matrix corresponding to the inputs in X.
• nepochs::Int=3000: Number of training epochs (default is 3000).
• lossfn::GenericLossFunction=MSLELoss(...): Loss function used for training. Defaults to Mean Squared Logarithmic Error.
• rng::AbstractRNG=Random.default_rng(): Random number generator used for parameter initialization.

Returns

• architecture: The trained neural network architecture.
• θ_pretrain: Trained parameters of the neural network.
• st_pretrain: Internal states of the trained model.

Notes

Pretrainign helps to reduce the number of total epochs required to train the UDE by selecting a physical meaningful initialization for the model. The function initializes the model parameters and states using Lux.setup, then performs training using a custom train_model! function with ADAM optimizer. Loss values are printed every 100 epochs during training.

Example

using Lux, Random
arch = Chain(Dense(10 => 20, relu), Dense(20 => 1))
X = rand(10, 100)
Y = rand(1, 100)
model, params, state = pretraining(arch; X = X, Y = Y)

printProgressLoss(iter, total_iters, loss, improvement)

Print function to track training.

random_matrix(mean, std, corr_length)

Generate a random matrix with entries drawn from a multivariate normal distribution whose mean is given by mean and whose covariance decays exponentially with grid distance.

Arguments

• mean::AbstractMatrix{<:Real}: Matrix specifying the spatial mean values at each grid point. Entries equal to 0.0 are treated as inactive and skipped in sampling.
• std::Real: Standard deviation scaling factor for the covariance kernel.
• corr_length::Real: Correlation length parameter controlling how fast correlations decay with Euclidean distance between grid points.

Returns

• H_sample::Matrix{Float64}: A random realization of the same size as mean, with correlated entries drawn from MvNormal(mean, Σ), where Σ[i,j] = std * exp(-‖coords[i] - coords[j]‖ / corr_length).

safe_slice(obj, ind::Integer)

Return a sliced object obj if ind > 0, otherwise return 0.0.

This function saves the results of an inversion to a file in JLD2 format. If the path argument is not provided, the function will create a default path based on the current project directory. The results are saved in a file named prediction_<nglaciers>glaciers_<tspan>.jld2, where <nglaciers> is the number of glaciers in the simulation and <tspan> is the simulation time span.

simulate_iceflow_UDE!(
    container::InversionBinder,
    cb::SciMLBase.DECallback,
    iceflow_prob::ODEProblem,
    tstops,
)

Make a forward simulation of the iceflow UDE.

splitθ(θ, glacier_idx::Integer, optimizableComponent::TrainableModel)

Given a ComponentVector θ, a glacier_idx and an optimizableComponent, extract the content of θ relevant for the given optimizableComponent and glacier ID glacier_idx.

stats_err_arrays(a::T, b::T) where T

Compute the ratio, the angle and the relative error between two arrays. The norm and scalar product are defined in the vectorial sense, meaning that this is mathematically equivalent to flatten the arrays. The arrays must be of the same type and have the same shape.

Arguments:

• a::T: First array.
• b::T: Second array to compare to the first one.

Returns:

• ratio::Float64: Ratio between the norm of the two arrays minus 1. Value close to zero means the arrays have approximately the same norm.
• angle::Float64: Scalar product between the two arrays normalized by the norm minus 1. Value close to zero means the arrays point towards the same direction.
• relerr::Float64: Relative error between the two arrays. Value close to zero means the arrays have approximately the same values. The norm of a is taken to normalize and compute the relative error.

ADAM optim

BFGS optim

train_UDE!(
    simulation::Inversion;
    save_every_iter::Bool = false,
    logger::Union{<: TBLogger, Nothing} = nothing
)

Trains UDE based on the current Inversion.

update_training_state!(simulation::S, l) where {S <: Simulation}

Update the training state to determine if the training has completed an epoch. If an epoch is completed, reset the minibatches, update the history loss, and increment the epoch count.

Arguments

• simulation: The current state of the simulation or training process.
• l: The current loss value or other relevant metric.

Returns

• None. This function updates the state in-place.

σ_zang(x; β = 2.0)

Smooth thresholding function for enforcing non-negativity and zero values for negative values following I. Zang, "A smoothing-out technique for min—max optimization" (1980).

Arguments

• x::Real: Input value to be thresholded.
• β::Real: (optional) Parameter controlling the transition zone width. Default is 2.0.

∂σ_zang(x; β = 2.0)

Derivative of the smooth thresholding function σ_zang.

Arguments

• x::Real: Input value to be thresholded.
• β::Real: (optional) Parameter controlling the transition zone width. Default is 2.0.

∇²(a::Matrix{F}, Δx::F, Δy::F) where {F<:AbstractFloat}

Computes the 2D Laplacian operator of a scalar field a on a regular grid using finite differences and staggered (dual–primal) averaging.

Arguments

• a::Matrix{F}: 2D scalar field to differentiate.
• Δx::F: Grid spacing in the x-direction.
• Δy::F: Grid spacing in the y-direction.

Returns

• Matrix{F}: Approximation of the Laplacian ∇²a with boundary values set to 0.0.