Ensuring Forward and Backward Differentiability in Classical Inversion Process

Hi everyone,

I am going to start off with some context and then proceed with describing my problem, so you get a better understanding of the scope.

I am currently working on creating a classical inversion in the field of glacier dynamics using the 2D - Shallow Ice Approximation (SIA). This inversion will eventually be added to the ODINN package. Shortly stated ODINN uses Global glacier model using Universal Differential Equations to model and discover processes of climate-glacier interactions.

We are currently working on a new API. Where the functionalities of ODINN are split up into:

  • Sleipnir (everything concerning glacier data and other tools)
  • Muninn (everything concerning mass balance)
  • Huginn (everything concerning glacier dynamics)
  • ODINN (everything concerning machine learning and inversions)

Now back to the problem: I’ve been trying to setup a classical inversion using this new structure. The process to run such an inversion should be as following:

working_dir = joinpath(dirname(Base.current_project()), "data")
if !ispath(working_dir)

rgi_ids = ["RGI60-11.01450"]

params = Parameters(OGGM = OGGMparameters(working_dir=working_dir,
                                                ice_thickness_source = "Farinotti19"),
                            simulation = SimulationParameters(use_MB=true,
                                                            use_iceflow = true,
                                                            velocities = true,
                                                            use_glathida_data = true,
                                                            tspan=(2010.0, 2015.0),
                                                            working_dir = working_dir,

glaciers = initialize_glaciers(rgi_ids, params)

model = Model(
            iceflow = SIA2Dmodel(params), 
            mass_balance = TImodel1(params),
            machine_learning = nothing,
inversion = Inversion(model, glaciers, params)

So basically: setting up work directory → choose glacier (rgi_ids) → setup parameters → retrieve glacier data → setup a model (which forward model and mass balance) → create an instance of this inversion → run the inversion.

The run_inversion! looks as following:

function run_inversion!(simulation::Inversion)


    println("Running forward in-place PDE ice flow inversion model...\n")
    inversion_params_list = @showprogress pmap((glacier_idx) -> invert_iceflow!(glacier_idx, simulation), 1:length(simulation.glaciers))

    simulation.inversion = inversion_params_list

    @everywhere GC.gc() # run garbage collector


function invert_iceflow!(glacier_idx::Int, simulation::Inversion) 
    model = simulation.model
    params = simulation.parameters
    glacier = simulation.glaciers[glacier_idx]
    glacier_id = isnothing(glacier.gdir) ? "unnamed" : glacier.rgi_id
    println("Processing glacier: ", glacier_id)
    Huginn.initialize_iceflow_model!(model.iceflow, glacier_idx, glacier, params)
    params.solver.tstops = Huginn.define_callback_steps(params.simulation.tspan, params.solver.step)
    stop_condition(u,t,integrator) = Sleipnir.stop_condition_tstops(u,t,integrator, params.solver.tstops) #closure
    function action!(integrator)
        if params.simulation.use_MB 
            # Compute mass balance
            MB_timestep!(model, glacier, params.solver.step, integrator.t)
            apply_MB_mask!(integrator.u, glacier, model.iceflow)
    cb_MB = DiscreteCallback(stop_condition, action!)

    # Trim observed velocity to match size predicted velocity 
    glacier.V=glacier.V[1:end-1, 1:end-1]
    function objective_function(x,p)

        # Run iceflow PDE for this glacier
        du = params.simulation.use_iceflow ? Huginn.SIA2D! : Huginn.noSIA2D!
        results = simulate_iceflow_PDE!(simulation, model, params, cb_MB; du = du)
        H = results.H[end]

        # Extract the non-zero elements from observed data
        non_zero_indices_H = glacier.H_glathida .!= 0

        # Calculate the MSE for H using only the non-zero elements
        mse_H = mean(((glacier.H_glathida[non_zero_indices_H] .- H[non_zero_indices_H]) .^ 2))
        return mse_H
    initial_values = [4.0e-17,3.0]  
    lower_bound = [8.5e-20, 2.5]    # Adjust as needed
    upper_bound = [8.0e-17, 4.2]    # Adjust as needed
    p = (glacier_idx,simulation)

    # TODO: fix right optimizing method
    optf = OptimizationFunction(objective_function,  Optimization.AutoForwardDiff()) 
    prob = OptimizationProblem(optf, initial_values, p,  lb = lower_bound,ub = upper_bound)

    # Import a solver package and solve the optimization problem
    sol = Optimization.solve(prob, Optim.BFGS())  

    optimized_A = sol[1]
    # TODO: change these to optimized values
    optimized_n = sol[2]
    optimized_C = 0.0

    inversion_params = InversionParams(optimized_A, optimized_n, optimized_C)
    return inversion_params

The problem arises in the optimization process, where I use Optimization.AutoForwardDiff(). Specifically for this case the problem lies with that simulate_iceflow_PDE! and all the underlying functions which are not compatible with dual numbers.

My goal is to make sure that the underlying functions and structures from Huginn and Sleipnir are differentiable both forward and backward.

What is the best way to reach this goal? Ideally we’d like a package-agnostic solution, which can work with the SciML ecosystem.

Here some extra links for more insight into some of the functions and structures used in Huginn.

What I already have tried is changing the predefined type of AbstractFloat in the SIA2D structure into Number, however it fell apart with matrices. As in Julia, Matrix{dual} and Matrix{AbstractFloat} are not subtypes of Matrix{Number}.

If any more clarification is needed please feel free to ask so.

Thanks in advance!

No they’re not, but they are subtypes of Matrix{<:Number}. See the warning in this section of the docs: