From dcfb65a949e53282d144dbe792c976f8f6e64be2 Mon Sep 17 00:00:00 2001 From: Alban Gossard Date: Fri, 17 Jan 2025 18:17:27 +0100 Subject: [PATCH] [WIP] removing pycall dependencies --- .github/workflows/CI.yml | 46 ++++---- Project.toml | 2 - README.md | 15 ++- scripts/MB_benchmark.jl | 14 +-- scripts/benchmark_ODINN.jl | 5 - scripts/diffusivity_inversion_experiment.jl | 8 +- scripts/new_API_tests.jl | 9 +- scripts/rheology_inversions.jl | 6 +- scripts/toy_model.jl | 14 +-- src/ODINN.jl | 23 ---- src/helpers/climate.jl | 120 +++++++++++--------- src/models/machine_learning/ML_utils.jl | 8 +- src/parameters/UDEparameters.jl | 16 +-- src/setup/config.jl | 27 ----- test/PDE_UDE_solve.jl | 26 ++--- test/inversion_test.jl | 20 ++-- test/params_construction.jl | 15 +-- test/runtests.jl | 1 - test/test.jl | 5 +- 19 files changed, 150 insertions(+), 230 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 35d9a669..64f3fc1b 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -33,30 +33,30 @@ jobs: uses: actions/setup-python@v2 with: python-version: ${{ matrix.python }} - - name: Create environment with micromamba 🐍🖤 - uses: mamba-org/setup-micromamba@v1 - with: - micromamba-version: '1.5.6-0' - environment-file: ./environment.yml - environment-name: oggm_env # it is recommendable to add both name and yml file. - init-shell: bash - cache-environment: true -# condarc-file: ./condarc.yml # If necessary, we can include .condarc to configure environment - - name: Test creation of environment with micromamba 🔧🐍🖤 - run: | - which python - conda env export - shell: bash -el {0} - - name: Update certifi 🔒🔑 - run: | - pip install --upgrade certifi - shell: bash -el {0} -# - name: Test OGGM installation 🔧🌎 -# run: pytest.oggm +# - name: Create environment with micromamba 🐍🖤 +# uses: mamba-org/setup-micromamba@v1 +# with: +# micromamba-version: '1.5.6-0' +# environment-file: ./environment.yml +# environment-name: oggm_env # it is recommendable to add both name and yml file. +# init-shell: bash +# cache-environment: true +# # condarc-file: ./condarc.yml # If necessary, we can include .condarc to configure environment +# - name: Test creation of environment with micromamba 🔧🐍🖤 +# run: | +# which python +# conda env export +# shell: bash -el {0} +# - name: Update certifi 🔒🔑 +# run: | +# pip install --upgrade certifi +# shell: bash -el {0} +# # - name: Test OGGM installation 🔧🌎 +# # run: pytest.oggm +# # shell: bash -el {0} +# - name: Set ENV Variables for PyCall.jl 🐍 📞 +# run: export PYTHON=/home/runner/micromamba/envs/oggm_env/bin/python # shell: bash -el {0} - - name: Set ENV Variables for PyCall.jl 🐍 📞 - run: export PYTHON=/home/runner/micromamba/envs/oggm_env/bin/python - shell: bash -el {0} - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} diff --git a/Project.toml b/Project.toml index 90d9b6f0..7a2c826a 100644 --- a/Project.toml +++ b/Project.toml @@ -32,7 +32,6 @@ PlotThemes = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" -PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -64,7 +63,6 @@ PlotThemes = "3.1" Plots = "1.2" Polynomials = "2, 3, 4" ProgressMeter = "1" -PyCall = "1.9" PyPlot = "2.11" Revise = "3.5" SciMLSensitivity = "7.20" diff --git a/README.md b/README.md index d111c9e2..83157d35 100644 --- a/README.md +++ b/README.md @@ -67,16 +67,19 @@ ODINN's architecture makes it really straightforward to retrieve all the necessa ```julia using ODINN +# Data are automatically downloaded, retrieve the local paths +rgi_paths = get_rgi_paths() + # We create the necessary parameters -params = Parameters(OGGM = OGGMparameters(working_dir=joinpath(homedir(), "OGGM/ODINN_tests")), - simulation = SimulationParameters(working_dir=working_dir, - tspan=(2010.0, 2015.0), - workers=5), +params = Parameters(simulation = SimulationParameters(working_dir=working_dir, + tspan=(2010.0, 2015.0), + workers=5, + rgi_paths=rgi_paths), hyper = Hyperparameters(batch_size=4, epochs=10, optimizer=ODINN.ADAM(0.01)), UDE = UDEparameters(target = "A") - ) + ) # We define which glacier RGI IDs we want to work with rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-08.00213", "RGI60-04.04351"] @@ -86,7 +89,7 @@ model = Model(iceflow = SIA2Dmodel(params), mass_balance = mass_balance = TImodel1(params; DDF=6.0/1000.0, acc_factor=1.2/1000.0), machine_learning = NN(params)) -# We initialize the glaciers with all the necessary data +# We initialize the glaciers with all the necessary data glaciers = initialize_glaciers(rgi_ids, params) # We specify the type of simulation we want to perform diff --git a/scripts/MB_benchmark.jl b/scripts/MB_benchmark.jl index 353b08ee..85fcb36c 100644 --- a/scripts/MB_benchmark.jl +++ b/scripts/MB_benchmark.jl @@ -61,23 +61,17 @@ function run_benchmark(analysis) tspan = (2010.0,2015.0) # period in years for simulation - # Configure OGGM settings in all workers - # Use a separate working dir to avoid conflicts with other simulations - working_dir = joinpath(homedir(), "Python/OGGM_memory_benchmark") - ######################################### ########### CLIMATE DATA ############## ######################################### if analysis == "timeit" - @timeit to "oggm config" oggm_config(working_dir; processes=1) - # Defining glaciers to be modelled with RGI IDs # RGI60-11.03638 # Argentière glacier # RGI60-11.01450 # Aletsch glacier rgi_ids = ["RGI60-11.03638", "RGI60-11.01450"] - + ### Initialize glacier directory to obtain DEM and ice thickness inversion ### gdirs = @timeit to "init_gdirs" init_gdirs(rgi_ids) @@ -109,13 +103,11 @@ function run_benchmark(analysis) elseif analysis == "profile" - @profview_allocs oggm_config(working_dir) - # Defining glaciers to be modelled with RGI IDs # RGI60-11.03638 # Argentière glacier # RGI60-11.01450 # Aletsch glacier rgi_ids = ["RGI60-11.03638", "RGI60-11.01450"] - + ### Initialize glacier directory to obtain DEM and ice thickness inversion ### gdirs = init_gdirs(rgi_ids) @@ -137,8 +129,6 @@ function run_benchmark(analysis) elseif analysis == "types" - oggm_config(working_dir) - # Defining glaciers to be modelled with RGI IDs # RGI60-11.03638 # Argentière glacier # RGI60-11.01450 # Aletsch glacier diff --git a/scripts/benchmark_ODINN.jl b/scripts/benchmark_ODINN.jl index 9f160468..02312b58 100644 --- a/scripts/benchmark_ODINN.jl +++ b/scripts/benchmark_ODINN.jl @@ -96,11 +96,6 @@ function run_benchmark() tspan = (2010.0,2015.0) # period in years for simulation - # Configure OGGM settings in all workers - # Use a separate working dir to avoid conflicts with other simulations - working_dir = joinpath(homedir(), "Python/OGGM_data_benchmark") - oggm_config(working_dir) - # Defining glaciers to be modelled with RGI IDs # RGI60-11.03638 # Argentière glacier # RGI60-11.01450 # Aletsch glacier diff --git a/scripts/diffusivity_inversion_experiment.jl b/scripts/diffusivity_inversion_experiment.jl index d466ff07..1805833e 100644 --- a/scripts/diffusivity_inversion_experiment.jl +++ b/scripts/diffusivity_inversion_experiment.jl @@ -42,14 +42,10 @@ function run() tspan = (2017, 2018) # period in years for simulation (also for spin-up) - # Configure OGGM settings in all workers - working_dir = joinpath(homedir(), "Python/OGGM_data_diffusivity") - oggm_config(working_dir) - # Defining glaciers to be modelled with RGI IDs rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-08.00213", "RGI60-04.04351", "RGI60-01.02170", - "RGI60-02.05098", "RGI60-01.01104", "RGI60-01.09162", "RGI60-01.00570", "RGI60-04.07051", - "RGI60-07.00274", "RGI60-07.01323", "RGI60-03.04207", "RGI60-03.03533", "RGI60-01.17316", + "RGI60-02.05098", "RGI60-01.01104", "RGI60-01.09162", "RGI60-01.00570", "RGI60-04.07051", + "RGI60-07.00274", "RGI60-07.01323", "RGI60-03.04207", "RGI60-03.03533", "RGI60-01.17316", "RGI60-07.01193", "RGI60-01.22174", "RGI60-14.07309", "RGI60-15.10261"] ### Initialize glacier directory to obtain DEM and ice thickness inversion ### diff --git a/scripts/new_API_tests.jl b/scripts/new_API_tests.jl index 61dbc87a..6b2f964a 100644 --- a/scripts/new_API_tests.jl +++ b/scripts/new_API_tests.jl @@ -12,6 +12,8 @@ nglaciers = 1 function API_test(tspan) + rgi_paths = get_rgi_paths() + to = get_timer("ODINN") reset_timer!(to) @@ -23,10 +25,11 @@ function API_test(tspan) use_MB=true, use_iceflow=true, multiprocessing=true, - workers=1), + workers=1, + rgi_paths=rgi_paths), physical=PhysicalParameters(A=2e-17), solver=SolverParameters(reltol=1e-7)) - # Create an ODINN model based on a 2D Shallow Ice Approximation, + # Create an ODINN model based on a 2D Shallow Ice Approximation, # a TI model with 1 DDF, and a neural network model = Model(iceflow = SIA2Dmodel(parameters), mass_balance = TImodel1(parameters; DDF=8.0/1000.0, acc_factor=1.0/1000.0), @@ -40,7 +43,7 @@ function API_test(tspan) prediction = Prediction(model, glaciers, parameters) # We run the simulation - @timeit get_timer("ODINN") "Run prediction" run!(prediction) + @timeit get_timer("ODINN") "Run prediction" run!(prediction) @show to display(ODINN.Makie.heatmap(prediction.results[1].S)) diff --git a/scripts/rheology_inversions.jl b/scripts/rheology_inversions.jl index dcd9253a..2f4ff6fd 100644 --- a/scripts/rheology_inversions.jl +++ b/scripts/rheology_inversions.jl @@ -25,7 +25,7 @@ processes = 10 ODINN.enable_multiprocessing(processes) # Flags ODINN.set_use_MB(false) -ODINN.make_plots(true) +ODINN.make_plots(true) # UDE training ODINN.set_train(true) # Train UDE ODINN.set_retrain(false) # Re-use previous NN weights to continue training @@ -34,10 +34,6 @@ function run() tspan = (2017, 2018) # period in years for simulation - # Configure OGGM settings in all workers - working_dir = joinpath(homedir(), "Python/OGGM_data_D") - oggm_config(working_dir) - gtd_file, rgi_ids = ODINN.get_glathida_path_and_IDs() @infiltrate diff --git a/scripts/toy_model.jl b/scripts/toy_model.jl index 099096a6..3276d115 100644 --- a/scripts/toy_model.jl +++ b/scripts/toy_model.jl @@ -18,7 +18,7 @@ using Plots using Infiltrator using Distributed using JLD2 -using Statistics +using Statistics # using AbbreviatedStackTraces # Activate to avoid GKS backend Plot issues in the JupyterHub @@ -32,8 +32,8 @@ function run_toy_model() ODINN.enable_multiprocessing(processes) # Flags ODINN.set_use_MB(true) - ODINN.make_plots(true) - # Spin up + ODINN.make_plots(true) + # Spin up ODINN.set_run_spinup(false) # Run the spin-up random_MB = generate_random_MB(gdirs_climate, tspan; plot=false)n ODINN.set_use_spinup(false) # Use the updated spinup # Reference simulations @@ -47,14 +47,10 @@ function run_toy_model() tspan = (2010.0, 2015.0) # period in years for simulation - # Configure OGGM settings in all workers - working_dir = joinpath(homedir(), "Python/OGGM_data") - oggm_config(working_dir) - # Defining glaciers to be modelled with RGI IDs rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-08.00213", "RGI60-04.04351", "RGI60-01.02170", - "RGI60-02.05098", "RGI60-01.01104", "RGI60-01.09162", "RGI60-01.00570", "RGI60-04.07051", - "RGI60-07.00274", "RGI60-07.01323", "RGI60-03.04207", "RGI60-03.03533", "RGI60-01.17316"]#, + "RGI60-02.05098", "RGI60-01.01104", "RGI60-01.09162", "RGI60-01.00570", "RGI60-04.07051", + "RGI60-07.00274", "RGI60-07.01323", "RGI60-03.04207", "RGI60-03.03533", "RGI60-01.17316"]#, # "RGI60-07.01193", "RGI60-01.22174", "RGI60-14.07309", "RGI60-15.10261"] # rgi_ids = rgi_ids[1:2] diff --git a/src/ODINN.jl b/src/ODINN.jl index eba402c8..dcc75182 100644 --- a/src/ODINN.jl +++ b/src/ODINN.jl @@ -26,7 +26,6 @@ Plots.theme(:wong2) # sets overall theme for Plots import Pkg using Distributed using ProgressMeter -using PyCall using Downloads using TimerOutputs using GeoStats @@ -41,28 +40,6 @@ const global root_dir::String = dirname(Base.current_project()) const global root_plots::String = joinpath(root_dir, "plots") -# ############################################## -# ############ PYTHON LIBRARIES ############## -# ############################################## - -# We either retrieve the reexported Python libraries from Sleipnir or we start from scratch -const netCDF4::PyObject = isdefined(Sleipnir, :netCDF4) ? Sleipnir.netCDF4 : PyNULL() -const cfg::PyObject = isdefined(Sleipnir, :cfg) ? Sleipnir.cfg : PyNULL() -const utils::PyObject = isdefined(Sleipnir, :utils) ? Sleipnir.utils : PyNULL() -const workflow::PyObject = isdefined(Sleipnir, :workflow) ? Sleipnir.workflow : PyNULL() -const tasks::PyObject = isdefined(Sleipnir, :tasks) ? Sleipnir.tasks : PyNULL() -const global_tasks::PyObject = isdefined(Sleipnir, :global_tasks) ? Sleipnir.global_tasks : PyNULL() -const graphics::PyObject = isdefined(Sleipnir, :graphics) ? Sleipnir.graphics : PyNULL() -const bedtopo::PyObject = isdefined(Sleipnir, :bedtopo) ? Sleipnir.bedtopo : PyNULL() -const millan22::PyObject = isdefined(Sleipnir, :millan22) ? Sleipnir.millan22 : PyNULL() -const MBsandbox::PyObject = isdefined(Sleipnir, :MBsandbox) ? Sleipnir.MBsandbox : PyNULL() -const salem::PyObject = isdefined(Sleipnir, :salem) ? Sleipnir.salem : PyNULL() - -# Essential Python libraries -const xr::PyObject = isdefined(Sleipnir, :xr) ? Sleipnir.xr : PyNULL() -const rioxarray::PyObject = isdefined(Sleipnir, :rioxarray) ? Sleipnir.rioxarray : PyNULL() -const pd::PyObject = isdefined(Sleipnir, :pd) ? Sleipnir.pd : PyNULL() - # ############################################## # ############ ODINN LIBRARIES ############### # ############################################## diff --git a/src/helpers/climate.jl b/src/helpers/climate.jl index d6a7e438..44d71830 100644 --- a/src/helpers/climate.jl +++ b/src/helpers/climate.jl @@ -75,7 +75,7 @@ end """ apply_t_grad!(climate, g_dem) -Applies temperature gradients to the glacier 2D climate data based on a DEM. +Applies temperature gradients to the glacier 2D climate data based on a DEM. """ function apply_t_cumul_grad!(climate, S) # We apply the gradients to the temperature @@ -96,24 +96,23 @@ function apply_t_grad!(climate, dem) end """ - downscale_2D_climate(climate, g_dem) + downscale_2D_climate(climate, S, S_coords) Projects climate data to the glacier matrix by simply copying the closest gridpoint to all matrix gridpoints. -Generates a new xarray Dataset which is returned. +Generates a new RasterStack which is returned. """ function downscale_2D_climate!(climate, S, S_coords) # Update 2D climate structure - # /!\ AVOID USING `.` IN JULIA TO ASSIGN. IT'S NOT HANDLED BY XARRAY. USE `=` INSTEAD - climate.climate_2D_step[].temp.data = climate.climate_step[].avg_temp.data .* ones(size(climate.climate_2D_step[].temp.data)) - climate.climate_2D_step[].PDD.data = climate.climate_step[].temp.data .* ones(size(climate.climate_2D_step[].PDD.data)) - climate.climate_2D_step[].snow.data = climate.climate_step[].prcp.data .* ones(size(climate.climate_2D_step[].snow.data)) - climate.climate_2D_step[].rain.data = climate.climate_step[].prcp.data .* ones(size(climate.climate_2D_step[].rain.data)) + climate.climate_2D_step.temp .= climate.climate_step.avg_temp.data + climate.climate_2D_step.PDD .= climate.climate_step.temp.data + climate.climate_2D_step.snow .= climate.climate_step.prcp.data + climate.climate_2D_step.rain .= climate.climate_step.prcp.data # Update gradients - climate.climate_2D_step[].gradient.data = climate.climate_step[].gradient.data - climate.climate_2D_step[].avg_gradient.data = climate.climate_step[].avg_gradient.data + climate.climate_2D_step.gradient .= climate.climate_step.gradient.data + climate.climate_2D_step.avg_gradient .= climate.climate_step.avg_gradient.data # Apply temperature gradients and compute snow/rain fraction for the selected period - apply_t_cumul_grad!(climate.climate_2D_step[], reshape(S, size(S))) # Reproject current S with xarray structure + apply_t_cumul_grad!(climate.climate_2D_step, reshape(S, size(S))) # Reproject current S with the RasterStack structure end function downscale_2D_climate(climate, S, S_coords) @@ -125,21 +124,33 @@ function downscale_2D_climate(climate, S, S_coords) rain_2D = climate.prcp.data .* dummy_grid # We generate a new dataset with the scaled data - climate_2D = xr.Dataset( - data_vars=Dict([ - ("temp", (["y","x"], temp_2D)), - ("PDD", (["y","x"], PDD_2D)), - ("snow", (["y","x"], snow_2D)), - ("rain", (["y","x"], rain_2D)), - ("gradient", climate.gradient.data), - ("avg_gradient", climate.avg_gradient.data) - ]), - coords=Dict([ - ("x", S_coords.x), - ("y", S_coords.y) - ]), - attrs=climate.attrs + data_vars=Dict( + "temp" => temp_2D, + "PDD" => PDD_2D, + "snow" => snow_2D, + "rain" => rain_2D, + "gradient" => climate.gradient.data, + "avg_gradient" => climate.avg_gradient.data + ) + climate_2D = RasterStack( + data_vars, (X(S_coords.x), Y(S_coords.y)) + metadata=metadata(climate) ) + # climate_2D = xr.Dataset( + # data_vars=Dict([ + # ("temp", (["y","x"], temp_2D)), + # ("PDD", (["y","x"], PDD_2D)), + # ("snow", (["y","x"], snow_2D)), + # ("rain", (["y","x"], rain_2D)), + # ("gradient", climate.gradient.data), + # ("avg_gradient", climate.avg_gradient.data) + # ]), + # coords=Dict([ + # ("x", S_coords.x), + # ("y", S_coords.y) + # ]), + # attrs=climate.attrs + # ) # Apply temperature gradients and compute snow/rain fraction for the selected period apply_t_cumul_grad!(climate_2D, reshape(S, size(S))) # Reproject current S with xarray structure @@ -150,7 +161,7 @@ end """ trim_period(period, climate) -Trims a time period based on the time range of a climate series. +Trims a time period based on the time range of a climate series. """ function trim_period(period, climate) if any(climate.time[1].dt.date.data[1] > period[1]) @@ -175,46 +186,47 @@ end partial_year(float::AbstractFloat) = partial_year(Day, float) -function get_longterm_temps(gdir::PyObject, tspan) - climate = xr.open_dataset(joinpath(gdir.dir, "raw_climate_$tspan.nc")) # load only once at the beginning - dem = rioxarray.open_rasterio(gdir.get_filepath("dem")) - apply_t_grad!(climate, dem) - longterm_temps = climate.groupby("time.year").mean().temp.data - return longterm_temps -end +# function get_longterm_temps(gdir::PyObject, tspan) +# climate = xr.open_dataset(joinpath(gdir.dir, "raw_climate_$tspan.nc")) # load only once at the beginning +# dem = rioxarray.open_rasterio(gdir.get_filepath("dem")) +# apply_t_grad!(climate, dem) +# longterm_temps = climate.groupby("time.year").mean().temp.data +# return longterm_temps +# end -function get_longterm_temps(gdir::PyObject, climate::PyObject) - dem = rioxarray.open_rasterio(gdir.get_filepath("dem")) - apply_t_grad!(climate, dem) - longterm_temps = climate.groupby("time.year").mean().temp.data - return longterm_temps -end +# function get_longterm_temps(gdir::PyObject, climate::PyObject) +# dem = rioxarray.open_rasterio(gdir.get_filepath("dem")) +# apply_t_grad!(climate, dem) +# longterm_temps = climate.groupby("time.year").mean().temp.data +# return longterm_temps +# end ### Data structures @kwdef mutable struct ClimateDataset - raw_climate::PyObject # Raw climate dataset for the whole simulation + raw_climate::RasterStack # Raw climate dataset for the whole simulation # Buffers to avoid memory allocations - climate_raw_step::Ref{PyObject} # Raw climate trimmed for the current step - climate_step::Ref{PyObject} # Climate data for the current step + climate_raw_step::Ref{RasterStack} # Raw climate trimmed for the current step + climate_step::Ref{RasterStack} # Climate data for the current step climate_2D_step::Ref{PyObject} # 2D climate data for the current step to feed to the MB model longterm_temps::Vector{Float64} # Longterm temperatures for the ice rheology - avg_temps::Ref{PyObject} # Intermediate buffer for computing average temperatures - avg_gradients::Ref{PyObject} # Intermediate buffer for computing average gradients + avg_temps::Ref{RasterStack} # Intermediate buffer for computing average temperatures + avg_gradients::Ref{RasterStack} # Intermediate buffer for computing average gradients end -function init_climate(gdir::PyObject, tspan, step, S, S_coords::PyObject) - dummy_period = partial_year(Day, tspan[1]):Day(1):partial_year(Day, tspan[1] + step) - raw_climate::PyObject = xr.open_dataset(joinpath(gdir.dir, "raw_climate_$tspan.nc")) - climate_step = Ref{PyObject}(get_cumulative_climate(raw_climate.sel(time=dummy_period))) - climate_2D_step = Ref{PyObject}(downscale_2D_climate(climate_step[], S, S_coords)) - longterm_temps = get_longterm_temps(gdir, raw_climate) +# function init_climate(gdir::PyObject, tspan, step, S, S_coords::PyObject) +function init_climate(rgi_id::String, params::Parameters, step, S, S_coords::Dict{String, Integer}) + rgi_path = params.simulation.rgi_paths[rgi_id] + dummy_period = partial_year(Day, params.simulation.tspan[1]):Day(1):partial_year(Day, params.simulation.tspan[1] + step) + raw_climate = RasterStack(joinpath(rgi_path, "raw_climate_$(params.simulation.tspan).nc")) + climate_step = get_cumulative_climate(raw_climate[At(dummy_period)]) + climate_2D_step = downscale_2D_climate(climate_step, S, S_coords) + longterm_temps = get_longterm_temps(rgi_id, raw_climate) climate = ClimateDataset(raw_climate = raw_climate, - climate_raw_step = raw_climate.sel(time=dummy_period), - #climate_cum_step = raw_climate.sel(time=dummy_period).sum(), + climate_raw_step = raw_climate[At(dummy_period)], climate_step = climate_step, climate_2D_step = climate_2D_step, longterm_temps = longterm_temps, - avg_temps = raw_climate.sel(time=dummy_period).temp.mean(), - avg_gradients = raw_climate.sel(time=dummy_period).gradient.mean()) + avg_temps = mean(raw_climate[At(dummy_period)].temp), + avg_gradients = mean(raw_climate[At(dummy_period)].gradient)) return climate end diff --git a/src/models/machine_learning/ML_utils.jl b/src/models/machine_learning/ML_utils.jl index 062d6f53..814260a0 100644 --- a/src/models/machine_learning/ML_utils.jl +++ b/src/models/machine_learning/ML_utils.jl @@ -98,7 +98,7 @@ function predict_diffusivity(UD_f, θ, X) end """ - generate_batches(batch_size, UD, target, gdirs_climate_batches, gdir_refs, context_batches; gtd_grids=nothing, shuffle=true) + generate_batches(simulation::S; shuffle=true) Generates batches for the UE inversion problem based on input data and feed them to the loss function. """ @@ -127,9 +127,9 @@ end # end """ - update_training_state(batch_size, n_gdirs) - -Update training state to know if the training has completed an epoch. + update_training_state(simulation, l) + +Update training state to know if the training has completed an epoch. If so, reset minibatches, update history loss and bump epochs. """ function update_training_state!(simulation, l) diff --git a/src/parameters/UDEparameters.jl b/src/parameters/UDEparameters.jl index 8262d490..09c1afde 100644 --- a/src/parameters/UDEparameters.jl +++ b/src/parameters/UDEparameters.jl @@ -51,7 +51,6 @@ include("InversionParameters.jl") Parameters(; physical::PhysicalParameters = PhysicalParameters(), simulation::SimulationParameters = SimulationParameters(), - OGGM::OGGMparameters = OGGMparameters(), solver::SolverParameters = SolverParameters(), hyper::Hyperparameters = Hyperparameters(), UDE::UDEparameters = UDEparameters() @@ -61,28 +60,23 @@ Initialize ODINN parameters Keyword arguments ================= - + """ function Parameters(; physical::PhysicalParameters = PhysicalParameters(), simulation::SimulationParameters = SimulationParameters(), - OGGM::OGGMparameters = OGGMparameters(), solver::SolverParameters = SolverParameters(), hyper::Hyperparameters = Hyperparameters(), UDE::UDEparameters = UDEparameters(), - inversion::InversionParameters = InversionParameters() - ) + inversion::InversionParameters = InversionParameters() + ) + - - parameters = Sleipnir.Parameters(physical, simulation, OGGM, hyper, solver, UDE, inversion) + parameters = Sleipnir.Parameters(physical, simulation, hyper, solver, UDE, inversion) if simulation.multiprocessing enable_multiprocessing(parameters) end - # Config OGGM *after* setting multirprocessing - # to ensure it is initialized in all workers - oggm_config(OGGM.working_dir; oggm_processes=OGGM.workers) - return parameters end diff --git a/src/setup/config.jl b/src/setup/config.jl index a83e124d..bd231cb9 100644 --- a/src/setup/config.jl +++ b/src/setup/config.jl @@ -9,33 +9,6 @@ function __init__() end end - - try - # Only load Python packages if not previously loaded by Sleipnir - if cfg == PyNULL() && workflow == PyNULL() && utils == PyNULL() && MBsandbox == PyNULL() - println("Initializing Python libraries in ODINN...") - # Load Python packages - copy!(netCDF4, pyimport("netCDF4")) - copy!(cfg, pyimport("oggm.cfg")) - copy!(utils, pyimport("oggm.utils")) - copy!(workflow, pyimport("oggm.workflow")) - copy!(tasks, pyimport("oggm.tasks")) - copy!(global_tasks, pyimport("oggm.global_tasks")) - copy!(graphics, pyimport("oggm.graphics")) - copy!(bedtopo, pyimport("oggm.shop.bedtopo")) - copy!(millan22, pyimport("oggm.shop.millan22")) - copy!(MBsandbox, pyimport("MBsandbox.mbmod_daily_oneflowline")) - copy!(salem, pyimport("salem")) - copy!(pd, pyimport("pandas")) - copy!(np, pyimport("numpy")) - copy!(xr, pyimport("xarray")) - copy!(rioxarray, pyimport("rioxarray")) - end - catch e - @warn "It looks like you have not installed and/or activated the virtual Python environment. \n - Please follow the guidelines in: https://github.com/ODINN-SciML/ODINN.jl#readme" - @warn exception=(e, catch_backtrace()) - end end diff --git a/test/PDE_UDE_solve.jl b/test/PDE_UDE_solve.jl index 230bab23..c71d9988 100644 --- a/test/PDE_UDE_solve.jl +++ b/test/PDE_UDE_solve.jl @@ -1,11 +1,10 @@ function ude_solve_test(atol; MB=false, fast=true) + rgi_paths = get_rgi_paths() + working_dir = joinpath(homedir(), "OGGM/ODINN_tests") - - # params = Parameters(OGGM = OGGMparameters(working_dir=working_dir, - # workers=3, - # multiprocessing=true), - # simulation = SimulationParameters(use_MB=MB, + + # params = Parameters(simulation = SimulationParameters(use_MB=MB, # velocities=false, # tspan=(2010.0, 2015.0), # workers=3, @@ -14,30 +13,29 @@ function ude_solve_test(atol; MB=false, fast=true) # hyper = Hyperparameters(batch_size=2, # epochs=10), # UDE = UDEparameters() - # ) + # ) - params = Parameters(OGGM = OGGMparameters(working_dir=working_dir, - multiprocessing=false), - simulation = SimulationParameters(working_dir=working_dir, + params = Parameters(simulation = SimulationParameters(working_dir=working_dir, use_MB=MB, velocities=true, tspan=(2010.0, 2015.0), multiprocessing=false, workers=5, - test_mode=true), + test_mode=true, + rgi_paths=rgi_paths), hyper = Hyperparameters(batch_size=4, epochs=4, optimizer=ODINN.ADAM(0.01)), UDE = UDEparameters(target = "A") - ) + ) - ## Retrieving gdirs and climate for the following glaciers - ## Fast version includes less glacier to reduce the amount of downloaded files and computation time on GitHub CI + ## Retrieving simulation data for the following glaciers + ## Fast version includes less glacier to reduce the amount of downloaded files and computation time on GitHub CI if fast rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-08.00213", "RGI60-04.04351"] else rgi_ids = ["RGI60-11.03638", "RGI60-11.01450", "RGI60-08.00213", "RGI60-04.04351", "RGI60-01.02170", - "RGI60-02.05098", "RGI60-01.01104", "RGI60-01.09162", "RGI60-01.00570", "RGI60-04.07051", + "RGI60-02.05098", "RGI60-01.01104", "RGI60-01.09162", "RGI60-01.00570", "RGI60-04.07051", "RGI60-07.00274", "RGI60-07.01323", "RGI60-01.17316"] end diff --git a/test/inversion_test.jl b/test/inversion_test.jl index 197d88bb..2e570f77 100644 --- a/test/inversion_test.jl +++ b/test/inversion_test.jl @@ -1,15 +1,9 @@ function inversion_test(;steady_state = false, save_refs = false) rgi_ids = ["RGI60-11.01450", "RGI60-11.00638"] + rgi_paths = get_rgi_paths() working_dir = joinpath(ODINN.root_dir, "test/data") params = Parameters( - OGGM = OGGMparameters( - working_dir = working_dir, - multiprocessing = true, - workers = 1, - ice_thickness_source = "Farinotti19", - DEM_source = "DEM3" - ), simulation = SimulationParameters( use_MB = true, use_iceflow = true, @@ -18,7 +12,9 @@ function inversion_test(;steady_state = false, save_refs = false) tspan = (2014.0, 2017.0), working_dir = working_dir, multiprocessing = true, - workers = 1 + workers = 1, + rgi_paths = rgi_paths, + ice_thickness_source = "Farinotti19", ), solver = SolverParameters(reltol = 1e-8) ) @@ -35,19 +31,19 @@ function inversion_test(;steady_state = false, save_refs = false) if steady_state run₀(inversion) ss = inversion.inversion - + if save_refs # Save the ss object to a JLD2 file when save_refs is true jldsave(joinpath(ODINN.root_dir, "test/data/PDE_refs_MB.jld2"); ss) end # Load the reference ss object from the JLD2 file for comparison ss_ref = load(joinpath(ODINN.root_dir, "test/data/PDE_refs_MB.jld2"), "ss") - + # Perform the comparison test between ss and ss_ref @test ss[1].H_pred == ss_ref[1].H_pred - + end - + end diff --git a/test/params_construction.jl b/test/params_construction.jl index b4bc2f42..c9838620 100644 --- a/test/params_construction.jl +++ b/test/params_construction.jl @@ -1,6 +1,8 @@ function params_constructor_specified(save_refs::Bool = false) + rgi_paths = Sleipnir.get_rgi_paths() + solver_params = SolverParameters(solver = Ralston(), reltol = 1e-8, step= 1.0/12.0, @@ -16,12 +18,6 @@ function params_constructor_specified(save_refs::Bool = false) epochs=10, batch_size=15) - oggm_params = OGGMparameters(working_dir=joinpath(homedir(), "OGGM/OGGM_data"), - multiprocessing=false, - workers=1, - ice_thickness_source="Millan22", - base_url="https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L1-L2_files/elev_bands/") - physical_params = PhysicalParameters(ρ = 900.0, g = 9.81, n = 3, @@ -44,7 +40,8 @@ function params_constructor_specified(save_refs::Bool = false) int_type = Int64, tspan = (2010.0,2015.0), multiprocessing = false, - workers = 10) + workers = 10, + rgi_paths = rgi_paths) ude_params = UDEparameters(sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)), optimization_method = "AD+AD", @@ -55,19 +52,17 @@ function params_constructor_specified(save_refs::Bool = false) hyper=hyparams, solver=solver_params, UDE=ude_params, - OGGM=oggm_params, simulation=simulation_params) if save_refs jldsave(joinpath(Sleipnir.root_dir, "test/data/params/solver_params.jld2"); solver_params) jldsave(joinpath(Sleipnir.root_dir, "test/data/params/hyparams.jld2"); hyparams) - jldsave(joinpath(Sleipnir.root_dir, "test/data/params/oggm_params.jld2"); oggm_params) jldsave(joinpath(Sleipnir.root_dir, "test/data/params/physical_params.jld2"); physical_params) jldsave(joinpath(Sleipnir.root_dir, "test/data/params/simulation_params.jld2"); simulation_params) jldsave(joinpath(Sleipnir.root_dir, "test/data/params/ude_params.jld2"); ude_params) jldsave(joinpath(Sleipnir.root_dir, "test/data/params/params.jld2"); params) end - + end diff --git a/test/runtests.jl b/test/runtests.jl index 70b086dd..b2104f22 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,6 @@ import Pkg Pkg.activate(dirname(Base.current_project())) -using PyCall # Update SSL certificate to avoid issue in GitHub Actions CI certifi = pyimport("certifi") ENV["SSL_CERT_FILE"] = certifi.where() diff --git a/test/test.jl b/test/test.jl index 52bdd7f7..5ddc47fe 100644 --- a/test/test.jl +++ b/test/test.jl @@ -11,7 +11,7 @@ using Test # using Random # using Distributed # using Statistics, LinearAlgebra, Random, Polynomials -# using HDF5 +# using HDF5 # using JLD2 # using OrdinaryDiffEq, DiffEqFlux # using Zygote: @ignore_derivatives @@ -20,8 +20,7 @@ using Test # using Infiltrator # using Plots # using ProgressMeter, ParallelDataTransfer -# using Dates -# using PyCall +# using Dates # using Makie, CairoMakie # @everywhere begin