Functional inversion tutorial
This tutorial provides a simple example on how to perform a functional inversion using Universal Differential Equations (UDEs) in ODINN.jl. For this, we will generate a synthetic dataset using a forward simulation, and then we will use this dataset to perform the functional inversion. The goal of this functional inversion will be to learn a synthetic law that maps A, i.e. the ice rigidity, to long-term changes in atmospheric surface temperature.
For more details on the functional inversion concept, please refer to the Functional Inversion section in the Inversion types page.
Running the whole code
using ODINN
using Random
# Define the working directory
working_dir = joinpath(ODINN.root_dir, "demos")
# We fetch the paths with the files for the available glaciers on disk
rgi_paths = get_rgi_paths()
# Ensure the working directory exists
mkpath(working_dir)
# Define the time step for the simulation output, in this case, a month.
δt = 1/12
# Define which glacier RGI IDs we want to work with
rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-11.02346", "RGI60-07.00065"]
params = Parameters(
simulation = SimulationParameters(
working_dir = working_dir,
use_MB = false,
use_velocities = true,
tspan = (2010.0, 2015.0),
multiprocessing = true,
workers = 4,
test_mode = false,
rgi_paths = rgi_paths,
gridScalingFactor = 4 # Downscale the glacier grid to speed-up this example
),
hyper = Hyperparameters(
batch_size = length(rgi_ids), # Set batch size equals size of the dataset
epochs = [25, 20],
optimizer = [
ODINN.ADAM(0.01),
ODINN.LBFGS(
linesearch = ODINN.LineSearches.BackTracking(iterations = 5)
)
]),
physical = PhysicalParameters(
minA = 8e-21,
maxA = 8e-17
),
UDE = UDEparameters(
optim_autoAD = ODINN.NoAD(),
grad = ContinuousAdjoint(),
optimization_method = "AD+AD",
empirical_loss_function = LossH() # Loss function based on ice thickness
),
solver = Huginn.SolverParameters(
step = 1 / 12, # Define the time step for the simulation output and for the adjoint
progress = true
)
)
# We define a synthetic law to generate the synthetic dataset.
# For this, we use some tabular data from Cuffey and Paterson (2010).
A_law = CuffeyPaterson(scalar = true)
model = Model(
iceflow = SIA2Dmodel(params; A = A_law),
mass_balance = TImodel1(params; DDF = 6.0 / 1000.0, acc_factor = 1.2 / 1000.0)
)
# We initialize the glaciers with all the necessary data
glaciers = initialize_glaciers(rgi_ids, params)
# Time snapshots for transient inversion
tstops = collect(2010:δt:2015)
# We generate the synthetic dataset using the forward simulation.
# This will generate a dataset with the ice thickness and surface velocities for each
# glacier at each time step. The dataset will be used to train the machine learning model.
glaciers = generate_ground_truth(glaciers, params, model, tstops)
# After this forward simulation, we restart the iceflow model to be ready for the inversions
nn_model = NeuralNetwork(params; seed = MersenneTwister(1))
A_law = LawA(nn_model, params)
model = Model(
iceflow = SIA2Dmodel(params; A = A_law),
mass_balance = TImodel1(params; DDF = 6.0 / 1000.0, acc_factor = 1.2 / 1000.0),
regressors = (; A = nn_model)
)
# We specify the type of simulation we want to perform
functional_inversion = Inversion(model, glaciers, params)
# And finally, we just run the simulation
run!(functional_inversion) From worker 3: Getting raw climate data for: RGI60-07.00065
[ Info: Running forward in-place PDE ice flow model
From worker 2: Processing glacier RGI60-11.03638 for PDE forward simulation
From worker 3: Processing glacier RGI60-11.01450 for PDE forward simulation
From worker 4: Processing glacier RGI60-11.02346 for PDE forward simulation
From worker 2: Processing glacier RGI60-07.00065 for PDE forward simulation
Progress: 50%|████████████████████▌ | ETA: 0:00:03
Progress: 75%|██████████████████████████████▊ | ETA: 0:00:01
Progress: 100%|█████████████████████████████████████████| Time: 0:00:03
[ Info: Optimizing with ADAM
[ Info: Optimizing with custom ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP} method
Iteration: [ 1 / 45] Loss:4.50095e+00
Iteration: [ 2 / 45] Loss:2.93485e+00 Improvement: -34.79 %
Iteration: [ 3 / 45] Loss:1.86853e+00 Improvement: -36.33 %
Iteration: [ 4 / 45] Loss:1.17342e+00 Improvement: -37.20 %
Iteration: [ 5 / 45] Loss:7.81282e-01 Improvement: -33.42 %
Iteration: [ 6 / 45] Loss:5.79836e-01 Improvement: -25.78 %
Iteration: [ 7 / 45] Loss:5.24758e-01 Improvement: -9.50 %
Iteration: [ 8 / 45] Loss:5.55571e-01 Improvement: 5.87 %
Iteration: [ 9 / 45] Loss:6.32397e-01 Improvement: 13.83 %
Iteration: [ 10 / 45] Loss:7.19986e-01 Improvement: 13.85 %
Iteration: [ 11 / 45] Loss:7.95153e-01 Improvement: 10.44 %
Iteration: [ 12 / 45] Loss:8.45073e-01 Improvement: 6.28 %
Iteration: [ 13 / 45] Loss:8.65607e-01 Improvement: 2.43 %
Iteration: [ 14 / 45] Loss:8.58543e-01 Improvement: -0.82 %
Iteration: [ 15 / 45] Loss:8.29092e-01 Improvement: -3.43 %
Iteration: [ 16 / 45] Loss:7.83983e-01 Improvement: -5.44 %
Iteration: [ 17 / 45] Loss:7.30153e-01 Improvement: -6.87 %
Iteration: [ 18 / 45] Loss:6.73900e-01 Improvement: -7.70 %
Iteration: [ 19 / 45] Loss:6.20371e-01 Improvement: -7.94 %
Iteration: [ 20 / 45] Loss:5.73312e-01 Improvement: -7.59 %
Iteration: [ 21 / 45] Loss:5.35020e-01 Improvement: -6.68 %
Iteration: [ 22 / 45] Loss:5.06144e-01 Improvement: -5.40 %
Iteration: [ 23 / 45] Loss:4.87405e-01 Improvement: -3.70 %
Iteration: [ 24 / 45] Loss:4.76882e-01 Improvement: -2.16 %
Iteration: [ 25 / 45] Loss:4.73280e-01 Improvement: -0.76 %
Iteration: [ 26 / 45] Loss:4.73280e-01 Improvement: 0.00 %
[ Info: Optimizing with BFGS
[ Info: Optimizing with custom ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP} method
Iteration: [ 27 / 45] Loss:4.73280e-01 Improvement: 0.00 %
Iteration: [ 28 / 45] Loss:4.64785e-01 Improvement: -1.79 %
Iteration: [ 29 / 45] Loss:4.48947e-01 Improvement: -3.41 %
Iteration: [ 30 / 45] Loss:3.52930e-01 Improvement: -21.39 %
Iteration: [ 31 / 45] Loss:3.28281e-01 Improvement: -6.98 %
Iteration: [ 32 / 45] Loss:2.49815e-01 Improvement: -23.90 %
Iteration: [ 33 / 45] Loss:1.17436e-01 Improvement: -52.99 %
Iteration: [ 34 / 45] Loss:7.99943e-02 Improvement: -31.88 %
Iteration: [ 35 / 45] Loss:1.60222e-02 Improvement: -79.97 %
Iteration: [ 36 / 45] Loss:1.26133e-02 Improvement: -21.28 %
Iteration: [ 37 / 45] Loss:1.22717e-02 Improvement: -2.71 %
Iteration: [ 38 / 45] Loss:1.20797e-02 Improvement: -1.56 %
Iteration: [ 39 / 45] Loss:1.14790e-02 Improvement: -4.97 %
Iteration: [ 40 / 45] Loss:9.90735e-03 Improvement: -13.69 %
Iteration: [ 41 / 45] Loss:9.50467e-03 Improvement: -4.06 %
Iteration: [ 42 / 45] Loss:6.70793e-03 Improvement: -29.42 %
Iteration: [ 43 / 45] Loss:4.80162e-03 Improvement: -28.42 %
Iteration: [ 44 / 45] Loss:3.31760e-03 Improvement: -30.91 %
Iteration: [ 45 / 45] Loss:1.21700e-03 Improvement: -63.32 %
Iteration: [ 46 / 45] Loss:2.37168e-04 Improvement: -80.51 %
Iteration: [ 47 / 45] Loss:1.37948e-04 Improvement: -41.84 %Step-by-step explanation of the tutorial
Here we will cover in detail each one of the steps that lead us to run the Inversion from the previous example. The goal of this simple example is to learn a mapping of a law for A, the creep coefficient of ice. Mathematically, we make A depends on the long term air temperature T through a neural network A=NN(T, θ) and we optimize θ so that the generated solution matches some ice thickness reference. This reference is generated using the relation of the book from Cuffey and Paterson (2010) [2].
Step 1: Parameter and glacier initialization
First we need to specify a list of RGI IDs of the glacier we want to work with. From these RGI IDs, we will look for the necessary files inside the workspace.
rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-11.02346", "RGI60-07.00065"]
rgi_paths = get_rgi_paths()Dict{String, String} with 56 entries:
"RGI60-11.00897" => "per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00897"
"RGI60-08.00213" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00213"
"RGI60-08.00147" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00147"
"RGI60-11.01270" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01270"
"RGI60-11.03646" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03646"
"RGI60-11.03232" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03232"
"RGI60-01.22174" => "per_glacier/RGI60-01/RGI60-01.22/RGI60-01.22174"
"RGI60-07.00274" => "per_glacier/RGI60-07/RGI60-07.00/RGI60-07.00274"
"RGI60-03.04207" => "per_glacier/RGI60-03/RGI60-03.04/RGI60-03.04207"
"RGI60-04.04351" => "per_glacier/RGI60-04/RGI60-04.04/RGI60-04.04351"
"RGI60-07.00065" => "per_glacier/RGI60-07/RGI60-07.00/RGI60-07.00065"
"RGI60-11.02346" => "per_glacier/RGI60-11/RGI60-11.02/RGI60-11.02346"
"RGI60-01.00570" => "per_glacier/RGI60-01/RGI60-01.00/RGI60-01.00570"
"RGI60-11.01238" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01238"
"RGI60-11.03005" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03005"
"RGI60-08.00087" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00087"
"RGI60-11.00787" => "per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00787"
"RGI60-08.00203" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00203"
"RGI60-07.01193" => "per_glacier/RGI60-07/RGI60-07.01/RGI60-07.01193"
⋮ => ⋮Then we need to define the parameters of the simulation we want to perform. The arguments are very similar to the ones used in the forward simulation tutorial and for a complete explanation, the reader should refer to this tutorial. The main difference with the forward simulation tutorial here is that we need to specify the parameters for the functional inversion through the Hyperparameters and the UDEparameters. The Hyperparameters structure contains information about the optimization algorithm. The UDEparameters define how the Universal Differential Equation (UDE) is solved and how its gradient is computed.
params = Parameters(
simulation = SimulationParameters(
working_dir = working_dir,
use_MB = false,
use_velocities = true,
tspan = (2010.0, 2015.0),
multiprocessing = true,
workers = 4,
test_mode = false,
rgi_paths = rgi_paths,
gridScalingFactor = 4 # Downscale the glacier grid to speed-up this example
),
hyper = Hyperparameters(
batch_size = length(rgi_ids), # Set batch size equals size of the dataset
epochs = [25, 20],
optimizer = [
ODINN.ADAM(0.01),
ODINN.LBFGS(
linesearch = ODINN.LineSearches.BackTracking(iterations = 5)
)
]),
physical = PhysicalParameters(
minA = 8e-21,
maxA = 8e-17
),
UDE = UDEparameters(
optim_autoAD = ODINN.NoAD(),
grad = ContinuousAdjoint(),
optimization_method = "AD+AD",
empirical_loss_function = LossH() # Loss function based on ice thickness
),
solver = Huginn.SolverParameters(
step = 1 / 12, # Define the time step for the simulation output and for the adjoint
progress = true
)
)Sleipnir.Parameters{PhysicalParameters{Float64}, SimulationParameters{Int64, Float64, MeanDateVelocityMapping}, Hyperparameters{Float64, Int64}, SolverParameters{Float64, Int64}, UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}, InversionParameters{Float64}}(PhysicalParameters{Float64}(900.0, 9.81, 1.0e-10, 1.0, 8.0e-17, 8.0e-21, 8.0e-17, 8.5e-20, 1.0, -25.0, 5.0e-18), SimulationParameters{Int64, Float64, MeanDateVelocityMapping}(false, true, true, true, 1.0, false, false, (2010.0, 2015.0), 0.08333333333333333, true, 4, "/home/runner/work/ODINN.jl/ODINN.jl/demos", false, Dict("RGI60-11.00897" => "per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00897", "RGI60-08.00213" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00213", "RGI60-08.00147" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00147", "RGI60-11.01270" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01270", "RGI60-11.03646" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03646", "RGI60-11.03232" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03232", "RGI60-01.22174" => "per_glacier/RGI60-01/RGI60-01.22/RGI60-01.22174", "RGI60-07.00274" => "per_glacier/RGI60-07/RGI60-07.00/RGI60-07.00274", "RGI60-03.04207" => "per_glacier/RGI60-03/RGI60-03.04/RGI60-03.04207", "RGI60-04.04351" => "per_glacier/RGI60-04/RGI60-04.04/RGI60-04.04351"…), "Farinotti19", MeanDateVelocityMapping(:nearest), 4), Hyperparameters{Float64, Int64}(1, 1, Float64[], Any[Optimisers.Adam(eta=0.01, beta=(0.9, 0.999), epsilon=1.0e-8), Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}(10, LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.BackTracking{Float64, Int64}(0.0001, 0.5, 0.1, 5, 3, Inf, nothing), nothing, Returns{Nothing}(nothing), Optim.Flat(), true)], 0.0, [25, 20], 4), SolverParameters{Float64, Int64}(OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-12, 0.08333333333333333, Float64[], false, true, 10, 100000), UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}(SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.EnzymeVJP}(SciMLSensitivity.EnzymeVJP(0), false), SciMLBase.NoAD(), ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}(DiscreteVJP{ADTypes.AutoMooncake{Nothing}}(ADTypes.AutoMooncake()), OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-8, 1.0e-8, 0.08333333333333333, :Linear, 200, EnzymeVJP()), "AD+AD", LossH{L2Sum{Int64}}(L2Sum{Int64}(3)), :A, :identity), InversionParameters{Float64}([1.0], [0.0], [Inf], [1, 1], 0.001, 0.001, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Optim.Flat}(LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}(0.1, 0.9, Inf, 5.0, 1.0e-6, 0.66, 50, 0.1, 0, Base.RefValue{Bool}(false), nothing, false), nothing, nothing, Optim.Flat())))Then, we initialize those glaciers based on those RGI IDs and the parameters we previously specified.
glaciers = initialize_glaciers(rgi_ids, params)4-element Vector{Glacier2D} distributed over regions 07 (x1), 11 (x3)
RGI60-11.03638 RGI60-11.01450 RGI60-11.02346 RGI60-07.00065
Step 2: Generate synthetic ground truth data with a forward simulation
The next step is to generate a synthetic dataset using a forward simulation. This will generate a dataset with the ice thickness and surface velocities for each glacier at each time step. The dataset will be used to train the machine learning model. We define a synthetic law to generate the synthetic dataset. For this, we use some tabular data from Cuffey and Paterson (2010) [2]. The REPL shows that it maps the long term air temperature T to the creep coefficient A.
A_law = CuffeyPaterson(scalar = true)(:T,) -> Array{Float64, 0} (↧@start )
The model is initialized using the Model constructor:
model = Model(
iceflow = SIA2Dmodel(params; A = A_law),
mass_balance = TImodel1(params; DDF = 6.0 / 1000.0, acc_factor = 1.2 / 1000.0)
)**** Model ****
SIA2D iceflow equation = ∇(D ∇S) with D = U H̄
and U = C (ρg)^(pq) H̄^(pq+1) ∇S^(p-1) + Γ H̄^(n+2) ∇S^(n-1)
Γ = 2A (ρg)^n /(n+2)
A: (:T,) -> Array{Float64, 0} (↧@start )
C: ConstantLaw -> Array{Float64, 0}
n: ConstantLaw -> Array{Float64, 0}
p: ConstantLaw -> Array{Float64, 0}
q: ConstantLaw -> Array{Float64, 0}
where
T => averaged_long_term_temperature
Temperature index mass balance model TImodel1
DDF = 0.006
acc_factor = 0.0012
No learnable components
***************We define the time snapshots for transient inversion, i.e. the time steps at which we want to save the results, which will be used to compute the adjoint in reverse mode.
tstops = collect(2010:δt:2015)
prediction = Prediction(model, glaciers, params)Prediction{Sleipnir.ModelCache{SIA2DCache{Float64, Int64, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, Array{Float64, 0}, Array{Float64, 0}, ScalarCacheNoVJP, ScalarCacheNoVJP}, Nothing}}(Sleipnir.Model{SIA2Dmodel{Float64, Law{ScalarCacheNoVJP, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, Huginn.var"#32#36"{Polynomials.Polynomial{Float64, :x}}}, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, typeof(Sleipnir.emptyVJPWithInputs)}, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, typeof(Sleipnir.emptyVJPWithInputs)}, Huginn.var"#33#37", Int64, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, typeof(Sleipnir.emptyPrepVJPWithInputs)}, DIVJP}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#9#10"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#11#12"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#13#14"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#15#16"}, NullLaw, NullLaw}, TImodel1{Float64}, Nothing}(SIA2D iceflow equation = ∇(D ∇S) with D = U H̄
and U = C (ρg)^(pq) H̄^(pq+1) ∇S^(p-1) + Γ H̄^(n+2) ∇S^(n-1)
Γ = 2A (ρg)^n /(n+2)
A: (:T,) -> Array{Float64, 0} (↧@start )
C: ConstantLaw -> Array{Float64, 0}
n: ConstantLaw -> Array{Float64, 0}
p: ConstantLaw -> Array{Float64, 0}
q: ConstantLaw -> Array{Float64, 0}
where
T => averaged_long_term_temperature
, Temperature index mass balance model TImodel1
DDF = 0.006
acc_factor = 0.0012, nothing), nothing, 4-element Vector{AbstractGlacier} distributed over regions 07 (x1), 11 (x3)
RGI60-11.03638 RGI60-11.01450 RGI60-11.02346 RGI60-07.00065
, Sleipnir.Parameters{PhysicalParameters{Float64}, SimulationParameters{Int64, Float64, MeanDateVelocityMapping}, Hyperparameters{Float64, Int64}, SolverParameters{Float64, Int64}, UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}, InversionParameters{Float64}}(PhysicalParameters{Float64}(900.0, 9.81, 1.0e-10, 1.0, 8.0e-17, 8.0e-21, 8.0e-17, 8.5e-20, 1.0, -25.0, 5.0e-18), SimulationParameters{Int64, Float64, MeanDateVelocityMapping}(false, true, true, true, 1.0, false, false, (2010.0, 2015.0), 0.08333333333333333, true, 4, "/home/runner/work/ODINN.jl/ODINN.jl/demos", false, Dict("RGI60-11.00897" => "per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00897", "RGI60-08.00213" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00213", "RGI60-08.00147" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00147", "RGI60-11.01270" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01270", "RGI60-11.03646" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03646", "RGI60-11.03232" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03232", "RGI60-01.22174" => "per_glacier/RGI60-01/RGI60-01.22/RGI60-01.22174", "RGI60-07.00274" => "per_glacier/RGI60-07/RGI60-07.00/RGI60-07.00274", "RGI60-03.04207" => "per_glacier/RGI60-03/RGI60-03.04/RGI60-03.04207", "RGI60-04.04351" => "per_glacier/RGI60-04/RGI60-04.04/RGI60-04.04351"…), "Farinotti19", MeanDateVelocityMapping(:nearest), 4), Hyperparameters{Float64, Int64}(1, 1, Float64[], Any[Optimisers.Adam(eta=0.01, beta=(0.9, 0.999), epsilon=1.0e-8), Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}(10, LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.BackTracking{Float64, Int64}(0.0001, 0.5, 0.1, 5, 3, Inf, nothing), nothing, Returns{Nothing}(nothing), Optim.Flat(), true)], 0.0, [25, 20], 4), SolverParameters{Float64, Int64}(OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-12, 0.08333333333333333, Float64[], false, true, 10, 100000), UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}(SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.EnzymeVJP}(SciMLSensitivity.EnzymeVJP(0), false), SciMLBase.NoAD(), ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}(DiscreteVJP{ADTypes.AutoMooncake{Nothing}}(ADTypes.AutoMooncake()), OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-8, 1.0e-8, 0.08333333333333333, :Linear, 200, EnzymeVJP()), "AD+AD", LossH{L2Sum{Int64}}(L2Sum{Int64}(3)), :A, :identity), InversionParameters{Float64}([1.0], [0.0], [Inf], [1, 1], 0.001, 0.001, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Optim.Flat}(LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}(0.1, 0.9, Inf, 5.0, 1.0e-6, 0.66, 50, 0.1, 0, Base.RefValue{Bool}(false), nothing, false), nothing, nothing, Optim.Flat()))), Sleipnir.Results[])We generate the synthetic dataset using the forward simulation. This will generate a dataset with the ice thickness and surface velocities for each glacier at each time step. The dataset will be used to train the machine learning model. This will run under the hood a Prediction using Huginn.jl.
glaciers = generate_ground_truth(glaciers, params, model, tstops)4-element Vector{Glacier2D} distributed over regions 07 (x1), 11 (x3)
RGI60-11.03638 RGI60-11.01450 RGI60-11.02346 RGI60-07.00065
The results of this simulation are stored in the thicknessData field of each glacier, which we can inspect:
glaciers[1].thicknessDataThicknessData{Float64}([2010.0, 2010.0833333333333, 2010.1666666666667, 2010.25, 2010.3333333333333, 2010.4166666666667, 2010.5, 2010.5833333333333, 2010.6666666666667, 2010.75 … 2014.25, 2014.3333333333333, 2014.4166666666667, 2014.5, 2014.5833333333333, 2014.6666666666667, 2014.75, 2014.8333333333333, 2014.9166666666667, 2015.0], [[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0] … [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]])Step 3: Model specification to perform a functional inversion
After this forward simulation, we define a new iceflow model to be ready for the inversions. The first step is to define a simple neural network that takes as input a scalar and returns a scalar.
nn_model = NeuralNetwork(params; seed = MersenneTwister(1))--- NeuralNetwork ---
architecture:
Chain(
layer_1 = Dense(1 => 3, #101), # 6 parameters
layer_2 = Dense(3 => 10, #102), # 40 parameters
layer_3 = Dense(10 => 3, #103), # 33 parameters
layer_4 = Dense(3 => 1, σ), # 4 parameters
) # Total: 83 parameters,
# plus 0 states.
θ: ComponentVector of length 83Then we define a law that uses this neural network to map the long term air temperature T to the creep coefficient A. ODINN comes with a set of already defined laws. Only a few of them support functional inversion as the computation of the gradient needs to be carefully handled. More information about these laws can be found in the laws tutorial.
law_input = (; T = iAvgScalarTemp())
A_law = LawA(nn_model, params, law_input)(:T,) -> Array{Float64, 0} (↧@start custom VJP ✅ precomputed)
Then we define an iceflow and ODINN tells us how the law is used in the iceflow equation.
iceflow = SIA2Dmodel(params; A = A_law)SIA2D iceflow equation = ∇(D ∇S) with D = U H̄
and U = C (ρg)^(pq) H̄^(pq+1) ∇S^(p-1) + Γ H̄^(n+2) ∇S^(n-1)
Γ = 2A (ρg)^n /(n+2)
A: (:T,) -> Array{Float64, 0} (↧@start custom VJP ✅ precomputed)
C: ConstantLaw -> Array{Float64, 0}
n: ConstantLaw -> Array{Float64, 0}
p: ConstantLaw -> Array{Float64, 0}
q: ConstantLaw -> Array{Float64, 0}
where
T => averaged_long_term_temperature
Finally we define the model which needs to know the iceflow and mass balance models, and in comparison to Huginn, there is a third argument regressors. The regressors argument tells how each regressor relates into the SIA. Although we already defined this in the law, this definition is mandatory for technical reasons. This argument will probably disappear in the future once the code becomes more mature. It must match how the laws are defined in the iceflow model.
model = Model(
iceflow = iceflow,
mass_balance = TImodel1(params; DDF = 6.0 / 1000.0, acc_factor = 1.2 / 1000.0),
regressors = (; A = nn_model)
)**** Model ****
SIA2D iceflow equation = ∇(D ∇S) with D = U H̄
and U = C (ρg)^(pq) H̄^(pq+1) ∇S^(p-1) + Γ H̄^(n+2) ∇S^(n-1)
Γ = 2A (ρg)^n /(n+2)
A: (:T,) -> Array{Float64, 0} (↧@start custom VJP ✅ precomputed)
C: ConstantLaw -> Array{Float64, 0}
n: ConstantLaw -> Array{Float64, 0}
p: ConstantLaw -> Array{Float64, 0}
q: ConstantLaw -> Array{Float64, 0}
where
T => averaged_long_term_temperature
Temperature index mass balance model TImodel1
DDF = 0.006
acc_factor = 0.0012
Learnable components
A: --- NeuralNetwork ---
architecture:
Chain(
layer_1 = Dense(1 => 3, #101), # 6 parameters
layer_2 = Dense(3 => 10, #102), # 40 parameters
layer_3 = Dense(10 => 3, #103), # 33 parameters
layer_4 = Dense(3 => 1, σ), # 4 parameters
) # Total: 83 parameters,
# plus 0 states.
θ: ComponentVector of length 83
***************Step 4: Train a Universal Differential Equation via a functional inversion
The next step is to specify the type of simulation we want to perform. In this case, we will use an Inversion simulation, which will use the synthetic dataset generated in the previous step to train a UDE.
functional_inversion = Inversion(model, glaciers, params)Inversion{Sleipnir.Model{SIA2Dmodel{Float64, Law{ScalarCache, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, ODINN.var"#234#240"{LuxCore.StatefulLuxLayerImpl.StatefulLuxLayer{Val{true}, Lux.Chain{@NamedTuple{layer_1::Lux.Dense{ODINN.var"#101#105", Int64, Int64, Nothing, Nothing, Static.True}, layer_2::Lux.Dense{ODINN.var"#102#106", Int64, Int64, Nothing, Nothing, Static.True}, layer_3::Lux.Dense{ODINN.var"#103#107", Int64, Int64, Nothing, Nothing, Static.True}, layer_4::Lux.Dense{typeof(NNlib.σ), Int64, Int64, Nothing, Nothing, Static.True}}, Nothing}, Nothing, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}, layer_4::@NamedTuple{}}}, Float64, Float64}}, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, typeof(Sleipnir.emptyVJPWithInputs)}, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, typeof(Sleipnir.emptyVJPWithInputs)}, ODINN.var"#238#244", Int64, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, ODINN.var"#236#242"{ODINN.var"#234#240"{LuxCore.StatefulLuxLayerImpl.StatefulLuxLayer{Val{true}, Lux.Chain{@NamedTuple{layer_1::Lux.Dense{ODINN.var"#101#105", Int64, Int64, Nothing, Nothing, Static.True}, layer_2::Lux.Dense{ODINN.var"#102#106", Int64, Int64, Nothing, Nothing, Static.True}, layer_3::Lux.Dense{ODINN.var"#103#107", Int64, Int64, Nothing, Nothing, Static.True}, layer_4::Lux.Dense{typeof(NNlib.σ), Int64, Int64, Nothing, Nothing, Static.True}}, Nothing}, Nothing, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}, layer_4::@NamedTuple{}}}, Float64, Float64}}}, CustomVJP}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#9#10"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#11#12"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#13#14"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#15#16"}, NullLaw, NullLaw}, TImodel1{Float64}, ODINN.TrainableComponents{NeuralNetwork{Lux.Chain{@NamedTuple{layer_1::Lux.Dense{ODINN.var"#101#105", Int64, Int64, Nothing, Nothing, Static.True}, layer_2::Lux.Dense{ODINN.var"#102#106", Int64, Int64, Nothing, Nothing, Static.True}, layer_3::Lux.Dense{ODINN.var"#103#107", Int64, Int64, Nothing, Nothing, Static.True}, layer_4::Lux.Dense{typeof(NNlib.σ), Int64, Int64, Nothing, Nothing, Static.True}}, Nothing}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(θ = ViewAxis(1:83, Axis(layer_1 = ViewAxis(1:6, Axis(weight = ViewAxis(1:3, ShapedAxis((3, 1))), bias = ViewAxis(4:6, Shaped1DAxis((3,))))), layer_2 = ViewAxis(7:46, Axis(weight = ViewAxis(1:30, ShapedAxis((10, 3))), bias = ViewAxis(31:40, Shaped1DAxis((10,))))), layer_3 = ViewAxis(47:79, Axis(weight = ViewAxis(1:30, ShapedAxis((3, 10))), bias = ViewAxis(31:33, Shaped1DAxis((3,))))), layer_4 = ViewAxis(80:83, Axis(weight = ViewAxis(1:3, ShapedAxis((1, 3))), bias = ViewAxis(4:4, Shaped1DAxis((1,))))))),)}}}, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}, layer_4::@NamedTuple{}}}, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyIC, SIA2D_A_target, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(A = ViewAxis(1:83, Axis(layer_1 = ViewAxis(1:6, Axis(weight = ViewAxis(1:3, ShapedAxis((3, 1))), bias = ViewAxis(4:6, Shaped1DAxis((3,))))), layer_2 = ViewAxis(7:46, Axis(weight = ViewAxis(1:30, ShapedAxis((10, 3))), bias = ViewAxis(31:40, Shaped1DAxis((10,))))), layer_3 = ViewAxis(47:79, Axis(weight = ViewAxis(1:30, ShapedAxis((3, 10))), bias = ViewAxis(31:33, Shaped1DAxis((3,))))), layer_4 = ViewAxis(80:83, Axis(weight = ViewAxis(1:3, ShapedAxis((1, 3))), bias = ViewAxis(4:4, Shaped1DAxis((1,))))))),)}}}}}, Sleipnir.ModelCache{SIA2DCache{Float64, Int64, ScalarCache, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, ScalarCacheNoVJP, Array{Float64, 0}, Array{Float64, 0}, ScalarCacheNoVJP, ScalarCacheNoVJP}, Nothing}, Glacier2D{Float64, Int64, Climate2D{Rasters.RasterStack{(:prcp, :temp, :gradient), @NamedTuple{prcp::Float64, temp::Float64, gradient::Float64}, 1, @NamedTuple{prcp::Vector{Float64}, temp::Vector{Float64}, gradient::Vector{Float64}}, Tuple{DimensionalData.Dimensions.Ti{DimensionalData.Dimensions.Lookups.Sampled{Dates.DateTime, Vector{Dates.DateTime}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}}}}, Tuple{}, @NamedTuple{prcp::Tuple{DimensionalData.Dimensions.Ti{Colon}}, temp::Tuple{DimensionalData.Dimensions.Ti{Colon}}, gradient::Tuple{DimensionalData.Dimensions.Ti{Colon}}}, DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}, @NamedTuple{prcp::DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}, temp::DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}, gradient::DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}}, Nothing}, Rasters.RasterStack{(:prcp, :temp, :gradient), @NamedTuple{prcp::Float64, temp::Float64, gradient::Float64}, 1, @NamedTuple{prcp::Vector{Float64}, temp::Vector{Float64}, gradient::Vector{Float64}}, Tuple{DimensionalData.Dimensions.Ti{DimensionalData.Dimensions.Lookups.Sampled{Dates.DateTime, Vector{Dates.DateTime}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}}}}, Tuple{}, @NamedTuple{prcp::Tuple{DimensionalData.Dimensions.Ti{Colon}}, temp::Tuple{DimensionalData.Dimensions.Ti{Colon}}, gradient::Tuple{DimensionalData.Dimensions.Ti{Colon}}}, DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}, @NamedTuple{prcp::DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}, temp::DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}, gradient::DimensionalData.Dimensions.Lookups.Metadata{Rasters.NCDsource, Dict{String, Any}}}, Nothing}, Sleipnir.ClimateStep{Float64}, Climate2Dstep{Float64}, Float64}, ThicknessData{Float64}, SurfaceVelocityData{Float64}}, Results{Sleipnir.Results{Float64, Int64}, TrainingStats{Float64, Int64}}}(Sleipnir.Model{SIA2Dmodel{Float64, Law{ScalarCache, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, ODINN.var"#234#240"{LuxCore.StatefulLuxLayerImpl.StatefulLuxLayer{Val{true}, Lux.Chain{@NamedTuple{layer_1::Lux.Dense{ODINN.var"#101#105", Int64, Int64, Nothing, Nothing, Static.True}, layer_2::Lux.Dense{ODINN.var"#102#106", Int64, Int64, Nothing, Nothing, Static.True}, layer_3::Lux.Dense{ODINN.var"#103#107", Int64, Int64, Nothing, Nothing, Static.True}, layer_4::Lux.Dense{typeof(NNlib.σ), Int64, Int64, Nothing, Nothing, Static.True}}, Nothing}, Nothing, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}, layer_4::@NamedTuple{}}}, Float64, Float64}}, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, typeof(Sleipnir.emptyVJPWithInputs)}, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, typeof(Sleipnir.emptyVJPWithInputs)}, ODINN.var"#238#244", Int64, Sleipnir.GenInputsAndApply{@NamedTuple{T::iAvgScalarTemp}, ODINN.var"#236#242"{ODINN.var"#234#240"{LuxCore.StatefulLuxLayerImpl.StatefulLuxLayer{Val{true}, Lux.Chain{@NamedTuple{layer_1::Lux.Dense{ODINN.var"#101#105", Int64, Int64, Nothing, Nothing, Static.True}, layer_2::Lux.Dense{ODINN.var"#102#106", Int64, Int64, Nothing, Nothing, Static.True}, layer_3::Lux.Dense{ODINN.var"#103#107", Int64, Int64, Nothing, Nothing, Static.True}, layer_4::Lux.Dense{typeof(NNlib.σ), Int64, Int64, Nothing, Nothing, Static.True}}, Nothing}, Nothing, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}, layer_4::@NamedTuple{}}}, Float64, Float64}}}, CustomVJP}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#9#10"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#11#12"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#13#14"}, ConstantLaw{ScalarCacheNoVJP, Huginn.var"#15#16"}, NullLaw, NullLaw}, TImodel1{Float64}, ODINN.TrainableComponents{NeuralNetwork{Lux.Chain{@NamedTuple{layer_1::Lux.Dense{ODINN.var"#101#105", Int64, Int64, Nothing, Nothing, Static.True}, layer_2::Lux.Dense{ODINN.var"#102#106", Int64, Int64, Nothing, Nothing, Static.True}, layer_3::Lux.Dense{ODINN.var"#103#107", Int64, Int64, Nothing, Nothing, Static.True}, layer_4::Lux.Dense{typeof(NNlib.σ), Int64, Int64, Nothing, Nothing, Static.True}}, Nothing}, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(θ = ViewAxis(1:83, Axis(layer_1 = ViewAxis(1:6, Axis(weight = ViewAxis(1:3, ShapedAxis((3, 1))), bias = ViewAxis(4:6, Shaped1DAxis((3,))))), layer_2 = ViewAxis(7:46, Axis(weight = ViewAxis(1:30, ShapedAxis((10, 3))), bias = ViewAxis(31:40, Shaped1DAxis((10,))))), layer_3 = ViewAxis(47:79, Axis(weight = ViewAxis(1:30, ShapedAxis((3, 10))), bias = ViewAxis(31:33, Shaped1DAxis((3,))))), layer_4 = ViewAxis(80:83, Axis(weight = ViewAxis(1:3, ShapedAxis((1, 3))), bias = ViewAxis(4:4, Shaped1DAxis((1,))))))),)}}}, @NamedTuple{layer_1::@NamedTuple{}, layer_2::@NamedTuple{}, layer_3::@NamedTuple{}, layer_4::@NamedTuple{}}}, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyTrainableModel, ODINN.emptyIC, SIA2D_A_target, ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(A = ViewAxis(1:83, Axis(layer_1 = ViewAxis(1:6, Axis(weight = ViewAxis(1:3, ShapedAxis((3, 1))), bias = ViewAxis(4:6, Shaped1DAxis((3,))))), layer_2 = ViewAxis(7:46, Axis(weight = ViewAxis(1:30, ShapedAxis((10, 3))), bias = ViewAxis(31:40, Shaped1DAxis((10,))))), layer_3 = ViewAxis(47:79, Axis(weight = ViewAxis(1:30, ShapedAxis((3, 10))), bias = ViewAxis(31:33, Shaped1DAxis((3,))))), layer_4 = ViewAxis(80:83, Axis(weight = ViewAxis(1:3, ShapedAxis((1, 3))), bias = ViewAxis(4:4, Shaped1DAxis((1,))))))),)}}}}}(SIA2D iceflow equation = ∇(D ∇S) with D = U H̄
and U = C (ρg)^(pq) H̄^(pq+1) ∇S^(p-1) + Γ H̄^(n+2) ∇S^(n-1)
Γ = 2A (ρg)^n /(n+2)
A: (:T,) -> Array{Float64, 0} (↧@start custom VJP ✅ precomputed)
C: ConstantLaw -> Array{Float64, 0}
n: ConstantLaw -> Array{Float64, 0}
p: ConstantLaw -> Array{Float64, 0}
q: ConstantLaw -> Array{Float64, 0}
where
T => averaged_long_term_temperature
, Temperature index mass balance model TImodel1
DDF = 0.006
acc_factor = 0.0012, A: --- NeuralNetwork ---
architecture:
Chain(
layer_1 = Dense(1 => 3, #101), # 6 parameters
layer_2 = Dense(3 => 10, #102), # 40 parameters
layer_3 = Dense(10 => 3, #103), # 33 parameters
layer_4 = Dense(3 => 1, σ), # 4 parameters
) # Total: 83 parameters,
# plus 0 states.
θ: ComponentVector of length 83
), nothing, 4-element Vector{Glacier2D} distributed over regions 07 (x1), 11 (x3)
RGI60-11.03638 RGI60-11.01450 RGI60-11.02346 RGI60-07.00065
, Sleipnir.Parameters{PhysicalParameters{Float64}, SimulationParameters{Int64, Float64, MeanDateVelocityMapping}, Hyperparameters{Float64, Int64}, SolverParameters{Float64, Int64}, UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}, InversionParameters{Float64}}(PhysicalParameters{Float64}(900.0, 9.81, 1.0e-10, 1.0, 8.0e-17, 8.0e-21, 8.0e-17, 8.5e-20, 1.0, -25.0, 5.0e-18), SimulationParameters{Int64, Float64, MeanDateVelocityMapping}(false, true, true, true, 1.0, false, false, (2010.0, 2015.0), 0.08333333333333333, true, 4, "/home/runner/work/ODINN.jl/ODINN.jl/demos", false, Dict("RGI60-11.00897" => "per_glacier/RGI60-11/RGI60-11.00/RGI60-11.00897", "RGI60-08.00213" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00213", "RGI60-08.00147" => "per_glacier/RGI60-08/RGI60-08.00/RGI60-08.00147", "RGI60-11.01270" => "per_glacier/RGI60-11/RGI60-11.01/RGI60-11.01270", "RGI60-11.03646" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03646", "RGI60-11.03232" => "per_glacier/RGI60-11/RGI60-11.03/RGI60-11.03232", "RGI60-01.22174" => "per_glacier/RGI60-01/RGI60-01.22/RGI60-01.22174", "RGI60-07.00274" => "per_glacier/RGI60-07/RGI60-07.00/RGI60-07.00274", "RGI60-03.04207" => "per_glacier/RGI60-03/RGI60-03.04/RGI60-03.04207", "RGI60-04.04351" => "per_glacier/RGI60-04/RGI60-04.04/RGI60-04.04351"…), "Farinotti19", MeanDateVelocityMapping(:nearest), 4), Hyperparameters{Float64, Int64}(1, 1, Float64[], Any[Optimisers.Adam(eta=0.01, beta=(0.9, 0.999), epsilon=1.0e-8), Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}(10, LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.BackTracking{Float64, Int64}(0.0001, 0.5, 0.1, 5, 3, Inf, nothing), nothing, Returns{Nothing}(nothing), Optim.Flat(), true)], 0.0, [25, 20], 4), SolverParameters{Float64, Int64}(OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-12, 0.08333333333333333, Float64[], false, true, 10, 100000), UDEparameters{ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}}(SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.EnzymeVJP}(SciMLSensitivity.EnzymeVJP(0), false), SciMLBase.NoAD(), ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP}(DiscreteVJP{ADTypes.AutoMooncake{Nothing}}(ADTypes.AutoMooncake()), OrdinaryDiffEqLowStorageRK.RDPK3Sp35{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), 1.0e-8, 1.0e-8, 0.08333333333333333, :Linear, 200, EnzymeVJP()), "AD+AD", LossH{L2Sum{Int64}}(L2Sum{Int64}(3)), :A, :identity), InversionParameters{Float64}([1.0], [0.0], [Inf], [1, 1], 0.001, 0.001, Optim.BFGS{LineSearches.InitialStatic{Float64}, LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}, Nothing, Nothing, Optim.Flat}(LineSearches.InitialStatic{Float64}(1.0, false), LineSearches.HagerZhang{Float64, Base.RefValue{Bool}}(0.1, 0.9, Inf, 5.0, 1.0e-6, 0.66, 50, 0.1, 0, Base.RefValue{Bool}(false), nothing, false), nothing, nothing, Optim.Flat()))), Results{Sleipnir.Results{Float64, Int64}, TrainingStats{Float64, Int64}}(Sleipnir.Results{Float64, Int64}[], TrainingStats{Float64, Int64}(nothing, Float64[], 0, nothing, ComponentArrays.ComponentVector[], ComponentArrays.ComponentVector[], nothing, Dates.DateTime("0000-01-01T00:00:00"))))And finally, we just run the simulation. This will run the adjoint method to compute the gradients and then use the ADAM and LBFGS optimizers to train the UDE model.
run!(functional_inversion)[ Info: Optimizing with ADAM
[ Info: Optimizing with custom ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP} method
Iteration: [ 1 / 45] Loss:4.50095e+00
Iteration: [ 2 / 45] Loss:2.93485e+00 Improvement: -34.79 %
Iteration: [ 3 / 45] Loss:1.86853e+00 Improvement: -36.33 %
Iteration: [ 4 / 45] Loss:1.17342e+00 Improvement: -37.20 %
Iteration: [ 5 / 45] Loss:7.81282e-01 Improvement: -33.42 %
Iteration: [ 6 / 45] Loss:5.79836e-01 Improvement: -25.78 %
Iteration: [ 7 / 45] Loss:5.24758e-01 Improvement: -9.50 %
Iteration: [ 8 / 45] Loss:5.55571e-01 Improvement: 5.87 %
Iteration: [ 9 / 45] Loss:6.32397e-01 Improvement: 13.83 %
Iteration: [ 10 / 45] Loss:7.19986e-01 Improvement: 13.85 %
Iteration: [ 11 / 45] Loss:7.95153e-01 Improvement: 10.44 %
Iteration: [ 12 / 45] Loss:8.45073e-01 Improvement: 6.28 %
Iteration: [ 13 / 45] Loss:8.65607e-01 Improvement: 2.43 %
Iteration: [ 14 / 45] Loss:8.58543e-01 Improvement: -0.82 %
Iteration: [ 15 / 45] Loss:8.29092e-01 Improvement: -3.43 %
Iteration: [ 16 / 45] Loss:7.83983e-01 Improvement: -5.44 %
Iteration: [ 17 / 45] Loss:7.30153e-01 Improvement: -6.87 %
Iteration: [ 18 / 45] Loss:6.73900e-01 Improvement: -7.70 %
Iteration: [ 19 / 45] Loss:6.20371e-01 Improvement: -7.94 %
Iteration: [ 20 / 45] Loss:5.73312e-01 Improvement: -7.59 %
Iteration: [ 21 / 45] Loss:5.35020e-01 Improvement: -6.68 %
Iteration: [ 22 / 45] Loss:5.06144e-01 Improvement: -5.40 %
Iteration: [ 23 / 45] Loss:4.87405e-01 Improvement: -3.70 %
Iteration: [ 24 / 45] Loss:4.76882e-01 Improvement: -2.16 %
Iteration: [ 25 / 45] Loss:4.73280e-01 Improvement: -0.76 %
Iteration: [ 26 / 45] Loss:4.73280e-01 Improvement: 0.00 %
[ Info: Optimizing with BFGS
[ Info: Optimizing with custom ContinuousAdjoint{Float64, Int64, DiscreteVJP{ADTypes.AutoMooncake{Nothing}}, EnzymeVJP} method
Iteration: [ 27 / 45] Loss:4.73280e-01 Improvement: 0.00 %
Iteration: [ 28 / 45] Loss:4.64785e-01 Improvement: -1.79 %
Iteration: [ 29 / 45] Loss:4.48947e-01 Improvement: -3.41 %
Iteration: [ 30 / 45] Loss:3.52930e-01 Improvement: -21.39 %
Iteration: [ 31 / 45] Loss:3.28281e-01 Improvement: -6.98 %
Iteration: [ 32 / 45] Loss:2.49815e-01 Improvement: -23.90 %
Iteration: [ 33 / 45] Loss:1.17436e-01 Improvement: -52.99 %
Iteration: [ 34 / 45] Loss:7.99943e-02 Improvement: -31.88 %
Iteration: [ 35 / 45] Loss:1.60222e-02 Improvement: -79.97 %
Iteration: [ 36 / 45] Loss:1.26133e-02 Improvement: -21.28 %
Iteration: [ 37 / 45] Loss:1.22717e-02 Improvement: -2.71 %
Iteration: [ 38 / 45] Loss:1.20797e-02 Improvement: -1.56 %
Iteration: [ 39 / 45] Loss:1.14790e-02 Improvement: -4.97 %
Iteration: [ 40 / 45] Loss:9.90735e-03 Improvement: -13.69 %
Iteration: [ 41 / 45] Loss:9.50467e-03 Improvement: -4.06 %
Iteration: [ 42 / 45] Loss:6.70793e-03 Improvement: -29.42 %
Iteration: [ 43 / 45] Loss:4.80162e-03 Improvement: -28.42 %
Iteration: [ 44 / 45] Loss:3.31760e-03 Improvement: -30.91 %
Iteration: [ 45 / 45] Loss:1.21700e-03 Improvement: -63.32 %
Iteration: [ 46 / 45] Loss:2.37168e-04 Improvement: -80.51 %
Iteration: [ 47 / 45] Loss:1.37948e-04 Improvement: -41.84 %The optimized parameters can be found in the results:
θ = functional_inversion.results.stats.θComponentVector{Float64}(A = (layer_1 = (weight = [0.5990673016430087; 0.5337316738400668; -0.8887685804068175;;], bias = [-0.2765035173900777, 0.3287094507689696, -0.9899957431822264]), layer_2 = (weight = [0.23343162696410316 -0.9326731957292596 -0.19459462441177894; 1.1604521262505447 0.31248796793951944 -0.7554525235532771; … ; -1.0486317191797747 1.1467912088798997 1.4990187699800928; 0.8415908178614244 -0.5203718473314279 0.11955892887513946], bias = [0.7089872749212571, 0.5627904361516916, -0.5709838821271023, -0.29308648299704704, 0.0729696512159626, 0.44036466001253993, 0.03754154131340662, -0.2942357864498496, 0.2621888810928649, 0.23873080086964016]), layer_3 = (weight = [0.8461128104328302 1.072199280583416 … -1.938911724247192 0.7082712740143283; -0.49887589333184323 -0.957346587782753 … 1.0743931472672132 0.10644832688831607; -0.40504204012897077 -0.30190115527511474 … -0.24841531143394543 -0.4490302546580924], bias = [0.571810315643326, -0.2672825292602745, 0.023396197030169187]), layer_4 = (weight = [1.5757090060882815 -0.24839852039221044 -0.22590695688637405], bias = [1.159614296198864])))And then we can visualize the learnt law by plotting the neural network mapping from temperature to A. We specify what is the law that acts as a ground truth to compare against.
plot_law(functional_inversion.model.iceflow.A, functional_inversion, law_input, θ;
ground_truth_law = prediction.model.iceflow.A)
Since we have a very limited range of temperature coverage with these 4 glaciers, we have only captured a small portion of the ground truth law, which has an almost linear form. If we wanted to learn the full dynamics of the Cuffey and Paterson synthetic law, we would need to have glaciers that cover a wider range of temperatures.
Just to have some additional context, here is how the synthetic law looks like for the full range of temperatures:
plot_law(prediction.model.iceflow.A, functional_inversion,
law_input, nothing; plot_full_input_range = true)
This page was generated using Literate.jl.