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)
mkdir("data")
end
rgi_ids = ["RGI60-11.01450"]
params = Parameters(OGGM = OGGMparameters(working_dir=working_dir,
multiprocessing=false,
workers=1,
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,
multiprocessing=false,
workers=1))
glaciers = initialize_glaciers(rgi_ids, params)
model = Model(
iceflow = SIA2Dmodel(params),
mass_balance = TImodel1(params),
machine_learning = nothing,
)
inversion = Inversion(model, glaciers, params)
run_inversion!(inversion)
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)
enable_multiprocessing(simulation.parameters)
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
end
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)
end
end
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)
model.iceflow.A=x[1]
model.iceflow.n=x[2]
# 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
end
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
end
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!