From 4152a15189e37f40a720e5fb171693f293b1b470 Mon Sep 17 00:00:00 2001 From: culi122 <40186418+cjuli1@users.noreply.github.com> Date: Tue, 25 Mar 2025 15:50:09 +0100 Subject: [PATCH 1/9] Allow config with only clearing, should be non-breaking change --- src/generic_io.jl | 12 ++ src/io.jl | 78 ++++++++----- src/run_serial.jl | 274 +++++++++++++++++++++++----------------------- 3 files changed, 200 insertions(+), 164 deletions(-) diff --git a/src/generic_io.jl b/src/generic_io.jl index 191a25b..6a4adad 100644 --- a/src/generic_io.jl +++ b/src/generic_io.jl @@ -21,6 +21,10 @@ end How price prognosis problems (ppp) are distributed on cores initially """ function get_dist_ppp(input::AbstractJulESInput) + if !haskey(input.settings["problems"], "prognosis") + return Tuple{ScenarioIx, CoreId}[] + end + cores = get_cores(input) N = length(cores) S = get_numscen_sim(input) @@ -38,6 +42,10 @@ end How end value problems (evp) are distributed on cores initially """ function get_dist_evp(input::AbstractJulESInput, subsystems::Vector{Tuple{SubsystemIx, AbstractSubsystem}}) + if !haskey(input.settings["problems"], "endvalue") + return Tuple{ScenarioIx, SubsystemIx, CoreId}[] + end + cores = get_cores(input) N = length(cores) #number of cores S = get_numscen_sim(input) #number of scenarios @@ -97,6 +105,10 @@ TODO: * = default """ function get_dist_stoch(input::AbstractJulESInput, subsystems::Vector{Tuple{SubsystemIx, AbstractSubsystem}}) + if !haskey(input.settings["problems"], "stochastic") + return (Tuple{SubsystemIx, CoreId}[], Tuple{ScenarioIx, SubsystemIx, CoreId}[]) + end + cores = get_cores(input) #Distributing the master problems/subsystems diff --git a/src/io.jl b/src/io.jl index b75fafb..1282816 100644 --- a/src/io.jl +++ b/src/io.jl @@ -354,7 +354,7 @@ end function get_numscen_sim(input::AbstractJulESInput) settings = get_settings(input) - if !isnothing(settings["scenariogeneration"]) + if haskey(settings, "scenariogeneration") if haskey(settings["scenariogeneration"], "simulation") return settings["scenariogeneration"]["simulation"]["numscen"] end @@ -363,7 +363,7 @@ function get_numscen_sim(input::AbstractJulESInput) end function get_numscen_stoch(input::AbstractJulESInput) settings = get_settings(input) - if !isnothing(settings["scenariogeneration"]) + if haskey(settings, "scenariogeneration") if haskey(settings["scenariogeneration"], "stochastic") return settings["scenariogeneration"]["stochastic"]["numscen"] end @@ -990,7 +990,11 @@ function update_output(t::TuLiPa.ProbTime, stepnr::Int) prob_results = db.mp[first(db.dist_mp[1])].prob else periodduration_power = parse_duration(settings["horizons"]["clearing"]["Power"], "periodduration") - periodduration_hydro = parse_duration(settings["horizons"]["clearing"]["Hydro"], "periodduration") + if haskey(settings["horizons"]["clearing"], "Hydro") + periodduration_hydro = parse_duration(settings["horizons"]["clearing"]["Hydro"], "periodduration") + else + periodduration_hydro = periodduration_power + end prob_results = db.cp.prob end numperiods_powerhorizon = Int(termduration.value / periodduration_power.value) @@ -1419,6 +1423,7 @@ function get_output_timing_local(data, steplength, skipmax) if has_result_times(settings) skipfactor = (skipmax+Millisecond(steplength))/Millisecond(steplength) + factors = [skipfactor,skipfactor,1] timing_cp = get_timing_cp_local() @@ -1426,6 +1431,13 @@ function get_output_timing_local(data, steplength, skipmax) for (scenix, values) in db.output.timing_ppp push!(timings_ppp, values) end + if haskey(settings["problems"], "prognosis") + dims = (get_steps(db), 3, 3, length(timings_ppp)) + timings_ppp1 = reshape(cat(timings_ppp..., dims=4), dims) + timings_ppp2 = transpose(dropdims(mean(timings_ppp1,dims=(1,4)),dims=(1,4))).*factors + else + timings_ppp2 = zeros(3, 3) + end # TODO: Add subix name df_evp = DataFrame([name => [] for name in ["scenix", "subix", "update", "solve", "total", "core", "skipmed"]]) @@ -1463,13 +1475,18 @@ function get_output_timing_local(data, steplength, skipmax) end df_mp[!, :mp_tot] = df_mp[!, :mp_s] + df_mp[!, :mp_u] + df_mp[!, :mp_fin] + df_mp[!, :mp_o] df_mp[df_mp.skipmed .== true, [:mp_u, :mp_s, :mp_fin, :mp_o, :mp_tot, :bend_it]] .= df_mp[df_mp.skipmed .== true, [:mp_u, :mp_s, :mp_fin, :mp_o, :mp_tot, :bend_it]] .* skipfactor - timings_mp = mean.(eachcol(select(df_mp, Not([:subix, :core, :skipmed, :mp_fin, :mp_o, :bend_it])))) - df_mp_core = combine(groupby(df_mp, [:core]), - :mp_u => sum => :mp_u, - :mp_s => sum => :mp_s, - :mp_fin => sum => :mp_fin, - :mp_o => sum => :mp_o, - :mp_tot => sum => :mp_tot) + if nrow(df_mp) != 0 + timings_mp = mean.(eachcol(select(df_mp, Not([:subix, :core, :skipmed, :mp_fin, :mp_o, :bend_it])))) + df_mp_core = combine(groupby(df_mp, [:core]), + :mp_u => sum => :mp_u, + :mp_s => sum => :mp_s, + :mp_fin => sum => :mp_fin, + :mp_o => sum => :mp_o, + :mp_tot => sum => :mp_tot) + else + timings_mp = [0.0, 0.0, 0.0] + df_mp_core = DataFrame([name => [] for name in ["core", "mp_u", "mp_s", "mp_fin", "mp_o", "mp_tot"]]) + end df_sp = DataFrame([name => [] for name in ["scenix", "subix", "update", "solve", "other", "core", "skipmed"]]) for (scenix, subix, core) in db.dist_sp @@ -1479,17 +1496,23 @@ function get_output_timing_local(data, steplength, skipmax) end df_sp[!, :total] = df_sp[!, :solve] + df_sp[!, :update] + df_sp[!, :other] df_sp[df_sp.skipmed .== true, [:update, :solve, :other, :total]] .= df_sp[df_sp.skipmed .== true, [:update, :solve, :other, :total]] .* skipfactor - df_sp_subix = combine(groupby(df_sp, [:subix]), - :update => sum => :sp_u, - :solve => sum => :sp_s, - :other => sum => :sp_o, - :total => sum => :sp_tot) - timings_sp = mean.(eachcol(select(df_sp_subix, Not([:subix, :sp_o])))) - df_sp_core = combine(groupby(df_sp, [:core]), - :update => sum => :sp_u, - :solve => sum => :sp_s, - :other => sum => :sp_o, - :total => sum => :sp_tot) + if nrow(df_sp) != 0 + df_sp_subix = combine(groupby(df_sp, [:subix]), + :update => sum => :sp_u, + :solve => sum => :sp_s, + :other => sum => :sp_o, + :total => sum => :sp_tot) + timings_sp = mean.(eachcol(select(df_sp_subix, Not([:subix, :sp_o])))) + df_sp_core = combine(groupby(df_sp, [:core]), + :update => sum => :sp_u, + :solve => sum => :sp_s, + :other => sum => :sp_o, + :total => sum => :sp_tot) + else + df_sp_subix = DataFrame([name => [] for name in ["subix", "sp_u", "sp_s", "sp_o", "sp_tot"]]) + timings_sp = [0.0, 0.0, 0.0] + df_sp_core = DataFrame([name => [] for name in ["core", "sp_u", "sp_s", "sp_o", "sp_tot"]]) + end # TODO: df_sp_scen df_subix = outerjoin(df_evp_subix, df_mp, df_sp_subix, on = :subix) @@ -1504,12 +1527,7 @@ function get_output_timing_local(data, steplength, skipmax) df_core = sort(df_core, :tot, rev=true) df_core = df_core[!, [:core, :tot, :evp_tot, :mp_tot, :sp_tot, :evp_u, :evp_s, :evp_o, :mp_u, :mp_s, :mp_fin, :mp_o, :sp_u, :sp_s, :sp_o]] - if haskey(settings["problems"], "prognosis") && haskey(settings["problems"], "clearing") - factors = [skipfactor,skipfactor,1] - dims = size(timings_ppp[1]) - dims = (dims..., length(timings_ppp)) - timings_ppp1 = reshape(cat(timings_ppp..., dims=4), dims) - timings_ppp2 = transpose(dropdims(mean(timings_ppp1,dims=(1,4)),dims=(1,4))).*factors + if haskey(settings["problems"], "clearing") all = vcat(timings_ppp2, reshape(timings_evp,1,3), reshape(timings_mp,1,3), reshape(timings_sp,1,3), mean(timing_cp, dims=1)) df = DataFrame(model=["long","med","short","evp","mp","sp","clearing"], update=all[:,1], solve=all[:,2], total=all[:,3]) df[!, :other] = df[!, :total] - df[!, :solve] - df[!, :update] @@ -1603,7 +1621,11 @@ function get_output_main_local() end else periodduration_power = parse_duration(settings["horizons"]["clearing"]["Power"], "periodduration") - periodduration_hydro = parse_duration(settings["horizons"]["clearing"]["Hydro"], "periodduration") + if haskey(settings["horizons"]["clearing"], "Hydro") + periodduration_hydro = parse_duration(settings["horizons"]["clearing"]["Hydro"], "periodduration") + else + periodduration_hydro = periodduration_power + end end # Only keep rhsterms that have at least one value (TODO: Do the same for sypply and demands) diff --git a/src/run_serial.jl b/src/run_serial.jl index 0a8f097..a0b5bb1 100644 --- a/src/run_serial.jl +++ b/src/run_serial.jl @@ -281,159 +281,161 @@ function create_subsystems(db) return push!(subsystems, ExogenSubsystem(commodities, endvaluemethod_sp)) else settings = get_settings(db.input) - method = settings["subsystems"]["function"] - if method == "twostorageduration" - storagesystems = TuLiPa.getstoragesystems(modelobjects) - shorttermstoragesystems = TuLiPa.getshorttermstoragesystems(storagesystems, Hour(settings["subsystems"]["shorttermstoragecutoff_hours"])) - for storagesystem in shorttermstoragesystems - commodities = get_commodities_from_storagesystem(storagesystem) - main = Set() - all = Set() - for obj in storagesystem - i, element = get_element_from_obj(elements, obj) - # println(getelkey(elements[i])) - push!(main, i) - push!(all, i) - end - - for (_i, _element) in enumerate(elements) - _deps = deep_dependencies[_element] - _add = false - for _dep in _deps - if _dep in main - _add = true - end + if haskey(settings, "subsystems") + method = settings["subsystems"]["function"] + if method == "twostorageduration" + storagesystems = TuLiPa.getstoragesystems(modelobjects) + shorttermstoragesystems = TuLiPa.getshorttermstoragesystems(storagesystems, Hour(settings["subsystems"]["shorttermstoragecutoff_hours"])) + for storagesystem in shorttermstoragesystems + commodities = get_commodities_from_storagesystem(storagesystem) + main = Set() + all = Set() + for obj in storagesystem + i, element = get_element_from_obj(elements, obj) + # println(getelkey(elements[i])) + push!(main, i) + push!(all, i) end - if _add + + for (_i, _element) in enumerate(elements) + _deps = deep_dependencies[_element] + _add = false for _dep in _deps - if !(_dep in all) - elkey = TuLiPa.getelkey(elements[_dep]) - if elkey.conceptname != TuLiPa.BALANCE_CONCEPT # getstoragesystems have already picked the balances we want to include, ignores power balances - # println(elkey) - push!(all, _dep) + if _dep in main + _add = true + end + end + if _add + for _dep in _deps + if !(_dep in all) + elkey = TuLiPa.getelkey(elements[_dep]) + if elkey.conceptname != TuLiPa.BALANCE_CONCEPT # getstoragesystems have already picked the balances we want to include, ignores power balances + # println(elkey) + push!(all, _dep) + end end end end end - end - # println(length(all)) + # println(length(all)) - shortstochduration = parse_duration(settings["subsystems"], "shortstochduration") - horizonterm_stoch = get_term_ppp(get_horizons(db.input), commodities, shortstochduration) + shortstochduration = parse_duration(settings["subsystems"], "shortstochduration") + horizonterm_stoch = get_term_ppp(get_horizons(db.input), commodities, shortstochduration) - priceareas = get_priceareas(storagesystem) - skipmed_impact = false - subsystem = StochSubsystem(commodities, priceareas, collect(all), horizonterm_stoch, shortstochduration, "startequalstop", skipmed_impact) - push!(subsystems, subsystem) - end - num_shortterm = length(subsystems) - println("Number of shortterm storagesystems $num_shortterm") - - longtermstoragesystems = TuLiPa.getlongtermstoragesystems(storagesystems, Hour(settings["subsystems"]["shorttermstoragecutoff_hours"])) - for storagesystem in longtermstoragesystems - commodities = get_commodities_from_storagesystem(storagesystem) - if length(commodities) == 1 - continue # TODO: error and fix dataset linvasselv and vakkerjordvatn have two subsystems, one not connected to power market, send liste til Carl - end - - # all = Set() - # for obj in storagesystem - # i, element = get_element_from_obj(elements, obj) - # for dep in deep_dependencies[element] - # if !(dep in all) - # println(getelkey(elements[dep])) - # push!(all, dep) - # end - # end - # end - - main = Set() - all = Set() - for obj in storagesystem - i, element = get_element_from_obj(elements, obj) - # println(getelkey(elements[i])) - push!(main, i) - push!(all, i) + priceareas = get_priceareas(storagesystem) + skipmed_impact = false + subsystem = StochSubsystem(commodities, priceareas, collect(all), horizonterm_stoch, shortstochduration, "startequalstop", skipmed_impact) + push!(subsystems, subsystem) end - - for (_i, _element) in enumerate(elements) - _deps = deep_dependencies[_element] - _add = false - for _dep in _deps - if _dep in main - _add = true - end + num_shortterm = length(subsystems) + println("Number of shortterm storagesystems $num_shortterm") + + longtermstoragesystems = TuLiPa.getlongtermstoragesystems(storagesystems, Hour(settings["subsystems"]["shorttermstoragecutoff_hours"])) + for storagesystem in longtermstoragesystems + commodities = get_commodities_from_storagesystem(storagesystem) + if length(commodities) == 1 + continue # TODO: error and fix dataset linvasselv and vakkerjordvatn have two subsystems, one not connected to power market, send liste til Carl + end + + # all = Set() + # for obj in storagesystem + # i, element = get_element_from_obj(elements, obj) + # for dep in deep_dependencies[element] + # if !(dep in all) + # println(getelkey(elements[dep])) + # push!(all, dep) + # end + # end + # end + + main = Set() + all = Set() + for obj in storagesystem + i, element = get_element_from_obj(elements, obj) + # println(getelkey(elements[i])) + push!(main, i) + push!(all, i) end - if _add + + for (_i, _element) in enumerate(elements) + _deps = deep_dependencies[_element] + _add = false for _dep in _deps - if !(_dep in all) - elkey = TuLiPa.getelkey(elements[_dep]) - if elkey.conceptname != TuLiPa.BALANCE_CONCEPT # getstoragesystems have already picked the balances we want to include, ignores power balances - # println(elkey) - push!(all, _dep) + if _dep in main + _add = true + end + end + if _add + for _dep in _deps + if !(_dep in all) + elkey = TuLiPa.getelkey(elements[_dep]) + if elkey.conceptname != TuLiPa.BALANCE_CONCEPT # getstoragesystems have already picked the balances we want to include, ignores power balances + # println(elkey) + push!(all, _dep) + end end end end end + # println(length(all)) + + # completed = Set() + # remaining = Set() + # for obj in storagesystem + # i, element = get_element_from_obj(elements, obj) + # push!(remaining, i) + # end + # while length(remaining) > 0 + # i = pop!(remaining) + # push!(completed, i) + # # # Deps over + # elkey = getelkey(elements[i]) + # println(elkey) + # for dep in dependencies[elkey] + # if !(dep in completed) && (dep <= length(elements)) + # if !(elkey.conceptname in ["Commodity", "Arrow"]) + # push!(remaining, dep) + # end + # end + # end + # # Deps under + # for (_i, _element) in enumerate(elements) + # _elkey = getelkey(_element) + # if (_elkey != elkey) + # _deps = dependencies[_elkey] + # for _dep in _deps + # if _dep == i + # if (_dep in completed) && !(_elkey.conceptname in ["TimeVector", "TimeValues"]) + # push!(remaining, _i) + # end + # end + # end + # end + # end + # end + # println(length(completed)) + + longstochduration = parse_duration(settings["subsystems"], "longstochduration") + horizonterm_stoch = get_term_ppp(get_horizons(db.input), commodities, longstochduration) + + priceareas = get_priceareas(storagesystem) + skipmed_impact = true + if has_longevduration(settings) + longevduration = parse_duration(settings["subsystems"], "longevduration") + horizonterm_evp = get_term_ppp(get_horizons(db.input), commodities, longevduration) + + subsystem = EVPSubsystem(commodities, priceareas, collect(all), horizonterm_evp, longevduration, horizonterm_stoch, longstochduration, "ppp", skipmed_impact) + else + subsystem = StochSubsystem(commodities, priceareas, collect(all), horizonterm_stoch, longstochduration, "ppp", skipmed_impact) + end + push!(subsystems, subsystem) end - # println(length(all)) - - # completed = Set() - # remaining = Set() - # for obj in storagesystem - # i, element = get_element_from_obj(elements, obj) - # push!(remaining, i) - # end - # while length(remaining) > 0 - # i = pop!(remaining) - # push!(completed, i) - # # # Deps over - # elkey = getelkey(elements[i]) - # println(elkey) - # for dep in dependencies[elkey] - # if !(dep in completed) && (dep <= length(elements)) - # if !(elkey.conceptname in ["Commodity", "Arrow"]) - # push!(remaining, dep) - # end - # end - # end - # # Deps under - # for (_i, _element) in enumerate(elements) - # _elkey = getelkey(_element) - # if (_elkey != elkey) - # _deps = dependencies[_elkey] - # for _dep in _deps - # if _dep == i - # if (_dep in completed) && !(_elkey.conceptname in ["TimeVector", "TimeValues"]) - # push!(remaining, _i) - # end - # end - # end - # end - # end - # end - # println(length(completed)) - - longstochduration = parse_duration(settings["subsystems"], "longstochduration") - horizonterm_stoch = get_term_ppp(get_horizons(db.input), commodities, longstochduration) - - priceareas = get_priceareas(storagesystem) - skipmed_impact = true - if has_longevduration(settings) - longevduration = parse_duration(settings["subsystems"], "longevduration") - horizonterm_evp = get_term_ppp(get_horizons(db.input), commodities, longevduration) - - subsystem = EVPSubsystem(commodities, priceareas, collect(all), horizonterm_evp, longevduration, horizonterm_stoch, longstochduration, "ppp", skipmed_impact) - else - subsystem = StochSubsystem(commodities, priceareas, collect(all), horizonterm_stoch, longstochduration, "ppp", skipmed_impact) - end - push!(subsystems, subsystem) + println("Number of longterm storagesystems $(length(subsystems)-num_shortterm)") + num_ignored = length(shorttermstoragesystems) + length(longtermstoragesystems) - length(subsystems) + println("Number of ignored storagesystems not connected to power $num_ignored") + else + error("getsubsystem() not implemented for $(method)") end - println("Number of longterm storagesystems $(length(subsystems)-num_shortterm)") - num_ignored = length(shorttermstoragesystems) + length(longtermstoragesystems) - length(subsystems) - println("Number of ignored storagesystems not connected to power $num_ignored") - else - error("getsubsystem() not implemented for $(method)") end end return subsystems From 227fb5602d6e0da933857d3ddfda651ae80e8538 Mon Sep 17 00:00:00 2001 From: pertft Date: Tue, 26 Aug 2025 12:50:19 +0200 Subject: [PATCH 2/9] added run jules function --- Project.toml | 1 + src/JulES.jl | 13 +++++-- src/run_jules_wrapper.jl | 73 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 85 insertions(+), 2 deletions(-) create mode 100644 src/run_jules_wrapper.jl diff --git a/Project.toml b/Project.toml index 503ad5f..c2f3ec2 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" diff --git a/src/JulES.jl b/src/JulES.jl index 2bce03b..a761ce2 100644 --- a/src/JulES.jl +++ b/src/JulES.jl @@ -2,7 +2,15 @@ module JulES import TuLiPa -using Distributed, Dates, Statistics, Clustering, Distributions, DataFrames, JSON +using Distributed +using Dates +using Statistics +using Clustering +using Distributions +using DataFrames +using JSON +using YAML +using HDF5 # Used by ifm using CSV @@ -12,7 +20,7 @@ using Lux using ComponentArrays using Interpolations using JLD2 -# TODO: Can remove the ones below because only used for training? +# Used by Nerual inflow model but not HBV # using DiffEqFlux # using SciMLSensitivity # using Optimization @@ -34,5 +42,6 @@ include("prob_util.jl") include("local_db.jl") include("run_serial.jl") include("scenariomodelling.jl") +include("run_jules_wrapper.jl") end \ No newline at end of file diff --git a/src/run_jules_wrapper.jl b/src/run_jules_wrapper.jl new file mode 100644 index 0000000..448102d --- /dev/null +++ b/src/run_jules_wrapper.jl @@ -0,0 +1,73 @@ +function getdataset(config, names, filename_clearing, filename_aggregated) + settings = config[config["main"]["settings"]] + + sti_dataset = joinpath(config["main"]["outputpath"]) + + clearing = JulES.JSON.parsefile(joinpath(sti_dataset, filename_clearing)) + clearing = JulES.TuLiPa.getelements(clearing) + + aggregated = JulES.JSON.parsefile(joinpath(sti_dataset, filename_aggregated)) + aggregated = JulES.TuLiPa.getelements(aggregated) + + timevectors = JulES.JSON.parsefile(joinpath(sti_dataset, names["FILENAME_DATAELEMENTS_TIMEVECTORS"] )) + timevectors = JulES.TuLiPa.getelements(timevectors, sti_dataset) + + elements = vcat(clearing, timevectors) + elements_ppp = vcat(aggregated, timevectors) + + storage_mapping = JulES.JSON.parsefile( + joinpath(sti_dataset, names["FILENAME_STORAGE_MAPPING"]), + dicttype=Dict{String, String}, + ) + + startmag_aggregated = JulES.JSON.parsefile( + joinpath(sti_dataset, names["FILENAME_START_STORAGES_AGGREGATED"]), + dicttype=Dict{String, Float64}, + ) + + startmag_clearing = JulES.JSON.parsefile( + joinpath(sti_dataset, names["FILENAME_START_STORAGES_CLEARING"]), + dicttype=Dict{String, Float64}, + ) + + return Dict( + "elements" => elements, + "elements_ppp" => elements_ppp, + "detailedrescopl" => storage_mapping, + "startmagdict" => startmag_clearing, + "aggstartmagdict" => startmag_aggregated, + ) +end + +function run_jules( + config_path, + datayear, + weatheryear, + outputpath, + JulESNames, + filename_clearing, + filename_aggregated + ) + + config = YAML.load_file(config_path) + + dataset = getdataset( + config, + JulESNames, + filename_clearing, + filename_aggregated + ) + + input = JulES.DefaultJulESInput(config, dataset, datayear, weatheryear) + + @time data = JulES.run_serial(input) + println("Total serial time above") + + println("Save output") + @time h5open(outputpath, "w") do file + for (k,v) in data + println(k) + write(file, k, v) + end + end +end From a4ddc6c7dedcadac18a0ceacac766a354b68ee98 Mon Sep 17 00:00:00 2001 From: culi122 <40186418+cjuli1@users.noreply.github.com> Date: Mon, 1 Sep 2025 19:05:41 +0200 Subject: [PATCH 3/9] skip updating end states for storages with 0 global energy eq --- src/prob_evp.jl | 5 +++++ src/prob_stoch.jl | 5 +++++ 2 files changed, 10 insertions(+) diff --git a/src/prob_evp.jl b/src/prob_evp.jl index 497e6f4..5d1706c 100644 --- a/src/prob_evp.jl +++ b/src/prob_evp.jl @@ -84,6 +84,11 @@ function update_endstates_evp(input, scenix::ScenarioIx, subix::SubsystemIx, evp balance = TuLiPa.getbalance(obj) bid = TuLiPa.getid(balance) instancename = split(TuLiPa.getinstancename(bid), "Balance_") + if haskey(enekvglobaldict, instancename[2]) + if enekvglobaldict[instancename[2]] == 0.0 + continue + end + end if haskey(detailedrescopl, instancename[2]) balancename = detailedrescopl[instancename[2]] bid = TuLiPa.Id(bid.conceptname, instancename[1] * "Balance_" * balancename * "_hydro_reservoir") # TODO: This should be in the dataset diff --git a/src/prob_stoch.jl b/src/prob_stoch.jl index 2219107..ff7c4f1 100644 --- a/src/prob_stoch.jl +++ b/src/prob_stoch.jl @@ -532,6 +532,11 @@ function update_endconditions_sp(scenix::ScenarioIx, subix::SubsystemIx, t::TuLi balance = TuLiPa.getbalance(obj) bid = TuLiPa.getid(balance) instancename = split(TuLiPa.getinstancename(bid), "Balance_") + if haskey(enekvglobaldict, instancename[2]) + if enekvglobaldict[instancename[2]] == 0.0 + continue + end + end if haskey(detailedrescopl, instancename[2]) balancename = detailedrescopl[instancename[2]] bid = TuLiPa.Id(bid.conceptname, instancename[1] * "Balance_" * balancename * "_hydro_reservoir") # TODO: This should be in the dataset From b94fab8d9863fe0e7dd42ce65eb7980132a03b03 Mon Sep 17 00:00:00 2001 From: culi122 <40186418+cjuli1@users.noreply.github.com> Date: Mon, 15 Sep 2025 16:43:11 +0200 Subject: [PATCH 4/9] Error if benders lower bound too low, and warning if no convergence --- src/prob_stoch.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/prob_stoch.jl b/src/prob_stoch.jl index ff7c4f1..b098a53 100644 --- a/src/prob_stoch.jl +++ b/src/prob_stoch.jl @@ -262,6 +262,9 @@ function solve_benders(stepnr::Int, subix::SubsystemIx) mp.cuts.scenslopes[scenix, cutix, :] .= scenslopes mp.cuts.scenconstants[scenix, cutix] = scenconstant end + if ub < mp.cuts.lower_bound # TODO customize the lb to each subsystem + error("Upper bound $ub is lower than lower bound $(mp.cuts.lower_bound). Decrease initial lower bound in settings.") + end TuLiPa.updatecutparameters!(mp.prob, mp.cuts) end @@ -274,6 +277,10 @@ function solve_benders(stepnr::Int, subix::SubsystemIx) end end end + if count == 15 + println("Warning: Benders did not converge within 15 iterations, reltol=$(reltol), ub=$(ub), lb=$(lb), step=$(stepnr), subix=$(subix)") + end + maintiming[5] = count return end From 44df3e7f957682d10a2b928181321babc018776f Mon Sep 17 00:00:00 2001 From: pertft Date: Wed, 17 Sep 2025 14:35:59 +0200 Subject: [PATCH 5/9] adding runtime load for ifm --- src/JulES.jl | 15 ++++++++------- src/run_jules_wrapper.jl | 16 ++++++++++++++-- 2 files changed, 22 insertions(+), 9 deletions(-) diff --git a/src/JulES.jl b/src/JulES.jl index a761ce2..7b41518 100644 --- a/src/JulES.jl +++ b/src/JulES.jl @@ -13,13 +13,14 @@ using YAML using HDF5 # Used by ifm -using CSV -using Random -using OrdinaryDiffEq -using Lux -using ComponentArrays -using Interpolations -using JLD2 +#using CSV +#using Random +#using OrdinaryDiffEq +#using Lux +#using ComponentArrays +#using Interpolations +#using JLD2 + # Used by Nerual inflow model but not HBV # using DiffEqFlux # using SciMLSensitivity diff --git a/src/run_jules_wrapper.jl b/src/run_jules_wrapper.jl index 448102d..9ea1b3a 100644 --- a/src/run_jules_wrapper.jl +++ b/src/run_jules_wrapper.jl @@ -39,6 +39,16 @@ function getdataset(config, names, filename_clearing, filename_aggregated) ) end +function load_ifm_dep() + @eval using CSV + @eval using Random + @eval using OrdinaryDiffEq + @eval using Lux + @eval using ComponentArrays + @eval using Interpolations + @eval using JLD2 +end + function run_jules( config_path, datayear, @@ -59,8 +69,10 @@ function run_jules( ) input = JulES.DefaultJulESInput(config, dataset, datayear, weatheryear) - - @time data = JulES.run_serial(input) + + has_ifm_results(input) && load_ifm_dep() + + @time data = JulES.run_serial(input) println("Total serial time above") println("Save output") From d617276903edd6aaf4412266e1c4698b6bd9981d Mon Sep 17 00:00:00 2001 From: pertft Date: Fri, 19 Sep 2025 10:41:43 +0200 Subject: [PATCH 6/9] using everywhere for ifm packages to be on all procs --- src/run_jules_wrapper.jl | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/run_jules_wrapper.jl b/src/run_jules_wrapper.jl index 9ea1b3a..256b8d1 100644 --- a/src/run_jules_wrapper.jl +++ b/src/run_jules_wrapper.jl @@ -40,13 +40,16 @@ function getdataset(config, names, filename_clearing, filename_aggregated) end function load_ifm_dep() - @eval using CSV - @eval using Random - @eval using OrdinaryDiffEq - @eval using Lux - @eval using ComponentArrays - @eval using Interpolations - @eval using JLD2 + mod = @__MODULE__ + @everywhere begin + @eval $mod using CSV + @eval $mod using Random + @eval $mod using OrdinaryDiffEq + @eval $mod using Lux + @eval $mod using ComponentArrays + @eval $mod using Interpolations + @eval $mod using JLD2 + end end function run_jules( @@ -70,7 +73,7 @@ function run_jules( input = JulES.DefaultJulESInput(config, dataset, datayear, weatheryear) - has_ifm_results(input) && load_ifm_dep() + has_ifm_results(input) && load_ifm_dep() @time data = JulES.run_serial(input) println("Total serial time above") From 5b1e1557f722ce6d467520b266574bd6b168f471 Mon Sep 17 00:00:00 2001 From: culi122 <40186418+cjuli1@users.noreply.github.com> Date: Mon, 22 Sep 2025 11:03:56 +0200 Subject: [PATCH 7/9] bugfix end states long prob --- src/prob_ppp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/prob_ppp.jl b/src/prob_ppp.jl index bea8234..2548682 100644 --- a/src/prob_ppp.jl +++ b/src/prob_ppp.jl @@ -151,7 +151,7 @@ function solve_ppp(t, steplength, stepnr, skipmed) if skipmed.value == 0 maintiming[3, 1] = @elapsed begin set_startstates!(p.longprob, TuLiPa.getstorages(TuLiPa.getobjects(p.longprob)), startstates) - setstartstates!(p.longprob, startstates) + set_endstates!(p.longprob, TuLiPa.getstorages(TuLiPa.getobjects(p.longprob)), startstates) maintiming[1, 1] = @elapsed TuLiPa.update!(p.longprob, scentime) maintiming[2, 1] = @elapsed TuLiPa.solve!(p.longprob) From 1610ba579d6530b4d39c7b84e9b4b800e1e2c904 Mon Sep 17 00:00:00 2001 From: pertft Date: Fri, 3 Oct 2025 13:38:32 +0200 Subject: [PATCH 8/9] bucket with weakdeps --- Project.toml | 17 ++++++----- src/ifm_bsd.jl => ext/IfmExt.jl | 34 ++++----------------- ext/IfmNeuralExt.jl | 52 +++++++++++++++++++++++++++++++++ src/JulES.jl | 11 ++----- src/ifm.jl | 13 +++++++++ src/run_jules_wrapper.jl | 23 ++++++++++----- 6 files changed, 98 insertions(+), 52 deletions(-) rename src/ifm_bsd.jl => ext/IfmExt.jl (50%) create mode 100644 ext/IfmNeuralExt.jl diff --git a/Project.toml b/Project.toml index c2f3ec2..17f3da6 100644 --- a/Project.toml +++ b/Project.toml @@ -4,24 +4,27 @@ repo = "https://github.com/NVE/JulES.jl.git" version = "0.4.0" [deps] -CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" -ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" -Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" -JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -Lux = "b2108857-7c20-44ae-9111-449ecde12c47" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TuLiPa = "970f5c25-cd7d-4f04-b50d-7a4fe2af6639" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" + +[weakdeps] +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" [compat] julia = "1.9.2" +OrdinaryDiffEq = "6.66.0" + +[extensions] +IfmExt = "OrdinaryDiffEq" \ No newline at end of file diff --git a/src/ifm_bsd.jl b/ext/IfmExt.jl similarity index 50% rename from src/ifm_bsd.jl rename to ext/IfmExt.jl index cfc0513..3adc3be 100644 --- a/src/ifm_bsd.jl +++ b/ext/IfmExt.jl @@ -5,15 +5,10 @@ and has BSD 3-Clause License # TODO: Update License file # TODO: Write license info in each source file -function initialize_NN_model() - rng = Random.default_rng() - NNmodel = Lux.Chain(Lux.Dense(4, 32, tanh), Lux.Dense(32,32, leakyrelu), Lux.Dense(32,32, leakyrelu), Lux.Dense(32,32, leakyrelu), - Lux.Dense(32,32, leakyrelu), Lux.Dense(32,5)) - p_NN_init, st_NN_init = Lux.setup(rng, NNmodel) - NN_in_fct(x, p) = NNmodel(x,p,st_NN_init)[1] - p_NN_init = ComponentArray(p_NN_init) - return NN_in_fct, p_NN_init -end +module IfmExt + +using OrdinaryDiffEq +using JulES step_fct(x) = (tanh(5.0*x) + 1.0)*0.5 Ps(P, T, Tmin) = step_fct(Tmin-T)*P @@ -24,7 +19,7 @@ ET(S1, T, Lday, Smax) = step_fct(S1)*step_fct(S1-Smax)*PET(T,Lday) + step_fct(S1 Qb(S1,f,Smax,Qmax) = step_fct(S1)*step_fct(S1-Smax)*Qmax + step_fct(S1)*step_fct(Smax-S1)*Qmax*exp(-f*(Smax-S1)) Qs(S1, Smax) = step_fct(S1)*step_fct(S1-Smax)*(S1-Smax) -function basic_bucket_incl_states(p_, itp_Lday, itp_P, itp_T, t_out) +function JulES.basic_bucket_incl_states(p_, itp_Lday, itp_P, itp_T, t_out) function exp_hydro_optim_states!(dS,S,ps,t) f, Smax, Qmax, Df, Tmax, Tmin = ps Lday = itp_Lday(t) @@ -44,23 +39,4 @@ function basic_bucket_incl_states(p_, itp_Lday, itp_P, itp_T, t_out) return Qout_, sol end -function NeuralODE_M100(p, norm_S0, norm_S1, norm_P, norm_T, itp_Lday, itp_P, itp_T, t_out, ann; S_init = [0.0, 0.0]) - function NeuralODE_M100_core!(dS,S,p,t) - Lday = itp_Lday(t) - P = itp_P(t) - T = itp_T(t) - g = ann([norm_S0(S[1]), norm_S1(S[2]), norm_P(P), norm_T(T)],p) - melting = relu(step_fct(S[1])*sinh(g[3])) - dS[1] = relu(sinh(g[4])*step_fct(-T)) - melting - dS[2] = relu(sinh(g[5])) + melting - step_fct(S[2])*Lday*exp(g[1])- step_fct(S[2])*exp(g[2]) - end - prob = ODEProblem(NeuralODE_M100_core!, S_init, Float64.((t_out[1], maximum(t_out))), p) - # sol = solve(prob, BS3(), dt=1.0, saveat=t_out, reltol=1e-3, abstol=1e-3, sensealg=BacksolveAdjoint(autojacvec=ZygoteVJP())) - sol = solve(prob, BS3(), dt=1.0, saveat=t_out, reltol=1e-3, abstol=1e-3) - P_interp = norm_P.(itp_P.(t_out)) - T_interp = norm_T.(itp_T.(t_out)) - S0_ = norm_S0.(sol[1,:]) - S1_ = norm_S1.(sol[2,:]) - Qout_ = exp.(ann(permutedims([S0_ S1_ P_interp T_interp]),p)[2,:]) - return Qout_, sol end \ No newline at end of file diff --git a/ext/IfmNeuralExt.jl b/ext/IfmNeuralExt.jl new file mode 100644 index 0000000..33f7e33 --- /dev/null +++ b/ext/IfmNeuralExt.jl @@ -0,0 +1,52 @@ +""" +This code is copied from https://github.com/marv-in/HydroNODE +and has BSD 3-Clause License +""" +# TODO: Update License file +# TODO: Write license info in each source file + +module IfmNeuralExt + +using DiffEqFlux +using SciMLSensitivity +using Optimization +using OptimizationOptimisers +using OptimizationBBO +using Zygote +using Lux +using Random +using CSV +using JulES + +function JulES.initialize_NN_model() + rng = Random.default_rng() + NNmodel = Lux.Chain(Lux.Dense(4, 32, tanh), Lux.Dense(32,32, leakyrelu), Lux.Dense(32,32, leakyrelu), Lux.Dense(32,32, leakyrelu), + Lux.Dense(32,32, leakyrelu), Lux.Dense(32,5)) + p_NN_init, st_NN_init = Lux.setup(rng, NNmodel) + NN_in_fct(x, p) = NNmodel(x,p,st_NN_init)[1] + p_NN_init = ComponentArray(p_NN_init) + return NN_in_fct, p_NN_init +end + +function JulES.NeuralODE_M100(p, norm_S0, norm_S1, norm_P, norm_T, itp_Lday, itp_P, itp_T, t_out, ann; S_init = [0.0, 0.0]) + function NeuralODE_M100_core!(dS,S,p,t) + Lday = itp_Lday(t) + P = itp_P(t) + T = itp_T(t) + g = ann([norm_S0(S[1]), norm_S1(S[2]), norm_P(P), norm_T(T)],p) + melting = relu(step_fct(S[1])*sinh(g[3])) + dS[1] = relu(sinh(g[4])*step_fct(-T)) - melting + dS[2] = relu(sinh(g[5])) + melting - step_fct(S[2])*Lday*exp(g[1])- step_fct(S[2])*exp(g[2]) + end + prob = ODEProblem(NeuralODE_M100_core!, S_init, Float64.((t_out[1], maximum(t_out))), p) + # sol = solve(prob, BS3(), dt=1.0, saveat=t_out, reltol=1e-3, abstol=1e-3, sensealg=BacksolveAdjoint(autojacvec=ZygoteVJP())) + sol = solve(prob, BS3(), dt=1.0, saveat=t_out, reltol=1e-3, abstol=1e-3) + P_interp = norm_P.(itp_P.(t_out)) + T_interp = norm_T.(itp_T.(t_out)) + S0_ = norm_S0.(sol[1,:]) + S1_ = norm_S1.(sol[2,:]) + Qout_ = exp.(ann(permutedims([S0_ S1_ P_interp T_interp]),p)[2,:]) + return Qout_, sol +end + +end \ No newline at end of file diff --git a/src/JulES.jl b/src/JulES.jl index 7b41518..44540f8 100644 --- a/src/JulES.jl +++ b/src/JulES.jl @@ -13,13 +13,9 @@ using YAML using HDF5 # Used by ifm -#using CSV -#using Random -#using OrdinaryDiffEq -#using Lux -#using ComponentArrays -#using Interpolations -#using JLD2 +using ComponentArrays +using Interpolations +using JLD2 # Used by Nerual inflow model but not HBV # using DiffEqFlux @@ -31,7 +27,6 @@ using HDF5 include("abstract_types.jl") include("dimension_types.jl") -include("ifm_bsd.jl") include("ifm.jl") include("generic_io.jl") include("io.jl") diff --git a/src/ifm.jl b/src/ifm.jl index 27800f5..12ca59f 100644 --- a/src/ifm.jl +++ b/src/ifm.jl @@ -821,3 +821,16 @@ function copy_elements_iprogtype(elements, iprogtype, ifm_names, ifm_derivedname @assert length(elements) == length(elements1) return elements1 end + +function initialize_NN_model() + error("Extension not loaded") +end + +function basic_bucket_incl_states(p_, itp_Lday, itp_P, itp_T, t_out) + error("Extension not loaded") +end + +function NeuralODE_M100(p, norm_S0, norm_S1, norm_P, norm_T, + itp_Lday, itp_P, itp_T, t_out, ann; S_init = [0.0, 0.0]) + error("Extension not loaded") +end \ No newline at end of file diff --git a/src/run_jules_wrapper.jl b/src/run_jules_wrapper.jl index 256b8d1..824d134 100644 --- a/src/run_jules_wrapper.jl +++ b/src/run_jules_wrapper.jl @@ -40,15 +40,22 @@ function getdataset(config, names, filename_clearing, filename_aggregated) end function load_ifm_dep() - mod = @__MODULE__ + if myid() == 1 + function ensure_packages(pkgs::Vector{String}) + deps = values(Pkg.dependencies()) + not_installed = filter(pkg -> !any(d -> d.name == pkg, deps), pkgs) + if !isempty(not_installed) + println("Installing missing packages: ", join(not_installed, ", ")) + Pkg.add(not_installed) + else + println("All packages already installed.") + end + end + ensure_packages(["OrdinaryDiffEq"]) + end + @everywhere begin - @eval $mod using CSV - @eval $mod using Random - @eval $mod using OrdinaryDiffEq - @eval $mod using Lux - @eval $mod using ComponentArrays - @eval $mod using Interpolations - @eval $mod using JLD2 + Base.eval(Main, :(using OrdinaryDiffEq)) end end From 525c8a642c5511c8fb5ea07a83a6256cc58b4594 Mon Sep 17 00:00:00 2001 From: pertft Date: Tue, 7 Oct 2025 13:20:32 +0200 Subject: [PATCH 9/9] moved more packages to extension --- Project.toml | 11 ++- ext/IfmExt.jl | 158 +++++++++++++++++++++++++++++++++++++++ src/JulES.jl | 6 +- src/ifm.jl | 154 ++------------------------------------ src/run_jules_wrapper.jl | 6 +- 5 files changed, 181 insertions(+), 154 deletions(-) diff --git a/Project.toml b/Project.toml index 17f3da6..0704bf4 100644 --- a/Project.toml +++ b/Project.toml @@ -15,16 +15,19 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TuLiPa = "970f5c25-cd7d-4f04-b50d-7a4fe2af6639" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" -JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" -Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" [weakdeps] OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" [compat] julia = "1.9.2" OrdinaryDiffEq = "6.66.0" +ComponentArrays = "0.15.17" +Interpolations = "0.16.1" +JLD2 = "0.5.15" [extensions] -IfmExt = "OrdinaryDiffEq" \ No newline at end of file +IfmExt = ["ComponentArrays", "Interpolations", "JLD2", "OrdinaryDiffEq"] \ No newline at end of file diff --git a/ext/IfmExt.jl b/ext/IfmExt.jl index 3adc3be..4c8d0f6 100644 --- a/ext/IfmExt.jl +++ b/ext/IfmExt.jl @@ -8,8 +8,14 @@ and has BSD 3-Clause License module IfmExt using OrdinaryDiffEq +using ComponentArrays +using Interpolations +using JLD2 +using Dates +using Statistics using JulES + step_fct(x) = (tanh(5.0*x) + 1.0)*0.5 Ps(P, T, Tmin) = step_fct(Tmin-T)*P Pr(P, T, Tmin) = step_fct(T-Tmin)*P @@ -19,6 +25,158 @@ ET(S1, T, Lday, Smax) = step_fct(S1)*step_fct(S1-Smax)*PET(T,Lday) + step_fct(S1 Qb(S1,f,Smax,Qmax) = step_fct(S1)*step_fct(S1-Smax)*Qmax + step_fct(S1)*step_fct(Smax-S1)*Qmax*exp(-f*(Smax-S1)) Qs(S1, Smax) = step_fct(S1)*step_fct(S1-Smax)*(S1-Smax) +function JulES.predict(m::JulES.TwoStateIfmHandler, u0::Vector{Float64}, t::JulES.TuLiPa.ProbTime) + JulES.update_prediction_data(m, m.updater, t) + + # create interpolation input functions + itp_method = SteffenMonotonicInterpolation() + itp_P = interpolate(m.data_pred.timepoints, m.data_pred.P, itp_method) + itp_T = interpolate(m.data_pred.timepoints, m.data_pred.T, itp_method) + itp_Lday = interpolate(m.data_pred.timepoints, m.data_pred.Lday, itp_method) + + (S0, G0) = u0 + (Q, __) = JulES.predict(m.predictor, S0, G0, itp_Lday, itp_P, itp_T, m.data_pred.timepoints) + + Q = Float64.(Q) + + Q .= Q .* m.m3s_per_mm + + return Q +end + +function JulES.estimate_u0(m::JulES.TwoStateIfmHandler, t::JulES.TuLiPa.ProbTime) + if isnothing(m.prev_t) + JulES._initial_data_obs_update(m, t) + else + JulES._data_obs_update(m, t) + end + m.prev_t = t + + # do prediction from start of obs up until today + (S0, G0) = (Float32(0), Float32(0)) + + # create interpolation input functions + itp_method = SteffenMonotonicInterpolation() + itp_P = interpolate(m.data_obs.timepoints, m.data_obs.P, itp_method) + itp_T = interpolate(m.data_obs.timepoints, m.data_obs.T, itp_method) + itp_Lday = interpolate(m.data_obs.timepoints, m.data_obs.Lday, itp_method) + + (__, OED_sol) = JulES.predict(m.predictor, S0, G0, itp_Lday, itp_P, itp_T, m.data_obs.timepoints) + + # extract states + est_S0 = Float64(last(OED_sol[1, :])) + est_G0 = Float64(last(OED_sol[2, :])) + + return [est_S0, est_G0] +end + +function JulES.calculate_normalize_factor(ifm_model) + start_with_buffer = ifm_model.handler.scen_start - Day(ifm_model.handler.ndays_obs) # Add ndays_obs days buffer + days_with_buffer = Day(ifm_model.handler.scen_stop - start_with_buffer) |> Dates.value + timepoints_with_buffer = (1:days_with_buffer) + days = Dates.value(Day(ifm_model.handler.scen_stop - ifm_model.handler.scen_start)) + timepoints_start = days_with_buffer - days + + P = zeros(length(timepoints_with_buffer)) + T = zeros(length(timepoints_with_buffer)) + Lday = zeros(length(timepoints_with_buffer)) + for i in timepoints_with_buffer + start = ifm_model.handler.scen_start + Day(i - 1) + P[i] = JulES.TuLiPa.getweightedaverage(ifm_model.handler.hist_P, start, JulES.ONEDAY_MS_TIMEDELTA) + T[i] = JulES.TuLiPa.getweightedaverage(ifm_model.handler.hist_T, start, JulES.ONEDAY_MS_TIMEDELTA) + Lday[i] = JulES.TuLiPa.getweightedaverage(ifm_model.handler.hist_Lday, start, JulES.ONEDAY_MS_TIMEDELTA) + end + + itp_method = SteffenMonotonicInterpolation() + itp_P = interpolate(timepoints_with_buffer, P, itp_method) + itp_T = interpolate(timepoints_with_buffer, T, itp_method) + itp_Lday = interpolate(timepoints_with_buffer, Lday, itp_method) + Q, _ = JulES.predict(ifm_model.handler.predictor, 0, 0, itp_Lday, itp_P, itp_T, timepoints_with_buffer) + Q = Float64.(Q)[timepoints_start:end] + Q .= Q .* ifm_model.handler.m3s_per_mm + return 1 / mean(Q) +end + +function JulES.common_includeTwoStateIfm!(Constructor, toplevel::Dict, lowlevel::Dict, elkey::JulES.TuLiPa.ElementKey, value::Dict) + JulES.TuLiPa.checkkey(toplevel, elkey) + + model_params = JulES.TuLiPa.getdictvalue(value, "ModelParams", String, elkey) + + moments = nothing + if haskey(value, "Moments") + moments = JulES.TuLiPa.getdictvalue(value, "Moments", String, elkey) + end + + hist_P = JulES.TuLiPa.getdictvalue(value, "HistoricalPercipitation", JulES.TuLiPa.TIMEVECTORPARSETYPES, elkey) + hist_T = JulES.TuLiPa.getdictvalue(value, "HistoricalTemperature", JulES.TuLiPa.TIMEVECTORPARSETYPES, elkey) + hist_Lday = JulES.TuLiPa.getdictvalue(value, "HistoricalDaylight", JulES.TuLiPa.TIMEVECTORPARSETYPES, elkey) + + ndays_pred = JulES.TuLiPa.getdictvalue(value, "NDaysPred", Real, elkey) + try + ndays_pred = Int(ndays_pred) + @assert ndays_pred >= 0 + catch e + error("Value for key NDaysPred must be positive integer for $elkey") + end + + basin_area = JulES.TuLiPa.getdictvalue(value, "BasinArea", Float64, elkey) + + deps = JulES.TuLiPa.Id[] + + all_ok = true + + (id, hist_P, ok) = JulES.TuLiPa.getdicttimevectorvalue(lowlevel, hist_P) + all_ok = all_ok && ok + JulES.TuLiPa._update_deps(deps, id, ok) + + (id, hist_T, ok) = JulES.TuLiPa.getdicttimevectorvalue(lowlevel, hist_T) + all_ok = all_ok && ok + JulES.TuLiPa._update_deps(deps, id, ok) + + (id, hist_Lday, ok) = JulES.TuLiPa.getdicttimevectorvalue(lowlevel, hist_Lday) + all_ok = all_ok && ok + JulES.TuLiPa._update_deps(deps, id, ok) + + if all_ok == false + return (false, deps) + end + + # TODO: Maybe make this user input in future? + updater = JulES.SimpleIfmDataUpdater() + + model_params = JLD2.load_object(model_params) + + is_nn = !isnothing(moments) + if is_nn + # convert model_params, stored with simpler data structure for stability between versions, + # into ComponentArray, which the NN-model needs + # (the simpler data structure is Vector{Tuple{Vector{Float32}, Vector{Float32}}}) + _subarray(i) = ComponentArray(weight = model_params[i][1], bias = model_params[i][2]) + _tuple(i) = (Symbol("layer_", i), _subarray(i)) + model_params = ComponentArray(NamedTuple(_tuple(i) for i in eachindex(model_params))) + # add moments, which is needed to normalize state inputs to the NN-model + moments = JLD2.load_object(moments) + model_params = (model_params, moments) + end + + data_forecast = nothing + data_obs = nothing + ndays_obs = 365 + data_forecast = nothing + ndays_forecast = 0 + + periodkey = JulES.TuLiPa.Id(JulES.TuLiPa.TIMEPERIOD_CONCEPT, "ScenarioTimePeriod") + period = lowlevel[periodkey] + scen_start = period["Start"] + scen_stop = period["Stop"] + + id = JulES.TuLiPa.getobjkey(elkey) + toplevel[id] = Constructor(id, model_params, updater, basin_area, hist_P, hist_T, hist_Lday, + ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast, scen_start, scen_stop) + + return (true, deps) +end + function JulES.basic_bucket_incl_states(p_, itp_Lday, itp_P, itp_T, t_out) function exp_hydro_optim_states!(dS,S,ps,t) f, Smax, Qmax, Df, Tmax, Tmin = ps diff --git a/src/JulES.jl b/src/JulES.jl index 44540f8..b485f9b 100644 --- a/src/JulES.jl +++ b/src/JulES.jl @@ -13,9 +13,9 @@ using YAML using HDF5 # Used by ifm -using ComponentArrays -using Interpolations -using JLD2 +#using ComponentArrays +#using Interpolations +#using JLD2 # Used by Nerual inflow model but not HBV # using DiffEqFlux diff --git a/src/ifm.jl b/src/ifm.jl index 12ca59f..4cbde2f 100644 --- a/src/ifm.jl +++ b/src/ifm.jl @@ -11,6 +11,7 @@ and integration with JulES # TODO: Replace all calls to Day(n) with Millisecond(86400000*n) for better performance const ONEDAY_MS_TIMEDELTA = TuLiPa.MsTimeDelta(Day(1)) +const extension_error_msg = "Missing optional dependency, cannot load extension." struct TwoStateIfmData P::Vector{Float32} @@ -176,48 +177,11 @@ function _data_obs_update(m::TwoStateIfmHandler, t::TuLiPa.ProbTime) end function estimate_u0(m::TwoStateIfmHandler, t::TuLiPa.ProbTime) - if isnothing(m.prev_t) - _initial_data_obs_update(m, t) - else - _data_obs_update(m, t) - end - m.prev_t = t - - # do prediction from start of obs up until today - (S0, G0) = (Float32(0), Float32(0)) - - # create interpolation input functions - itp_method = SteffenMonotonicInterpolation() - itp_P = interpolate(m.data_obs.timepoints, m.data_obs.P, itp_method) - itp_T = interpolate(m.data_obs.timepoints, m.data_obs.T, itp_method) - itp_Lday = interpolate(m.data_obs.timepoints, m.data_obs.Lday, itp_method) - - (__, OED_sol) = predict(m.predictor, S0, G0, itp_Lday, itp_P, itp_T, m.data_obs.timepoints) - - # extract states - est_S0 = Float64(last(OED_sol[1, :])) - est_G0 = Float64(last(OED_sol[2, :])) - - return [est_S0, est_G0] + error(extension_error_msg) end function predict(m::TwoStateIfmHandler, u0::Vector{Float64}, t::TuLiPa.ProbTime) - update_prediction_data(m, m.updater, t) - - # create interpolation input functions - itp_method = SteffenMonotonicInterpolation() - itp_P = interpolate(m.data_pred.timepoints, m.data_pred.P, itp_method) - itp_T = interpolate(m.data_pred.timepoints, m.data_pred.T, itp_method) - itp_Lday = interpolate(m.data_pred.timepoints, m.data_pred.Lday, itp_method) - - (S0, G0) = u0 - (Q, __) = predict(m.predictor, S0, G0, itp_Lday, itp_P, itp_T, m.data_pred.timepoints) - - Q = Float64.(Q) - - Q .= Q .* m.m3s_per_mm - - return Q + error(extension_error_msg) end # TODO: Add AutoCorrIfmDataUpdater that use value w(t)*x(t0) + (1-w(t-t0))*x(t), where w(0) = 1 and w -> 0 for larger inputs @@ -337,86 +301,7 @@ function includeTwoStateNeuralODEIfm!(toplevel::Dict, lowlevel::Dict, elkey::TuL end function common_includeTwoStateIfm!(Constructor, toplevel::Dict, lowlevel::Dict, elkey::TuLiPa.ElementKey, value::Dict) - TuLiPa.checkkey(toplevel, elkey) - - model_params = TuLiPa.getdictvalue(value, "ModelParams", String, elkey) - - moments = nothing - if haskey(value, "Moments") - moments = TuLiPa.getdictvalue(value, "Moments", String, elkey) - end - - hist_P = TuLiPa.getdictvalue(value, "HistoricalPercipitation", TuLiPa.TIMEVECTORPARSETYPES, elkey) - hist_T = TuLiPa.getdictvalue(value, "HistoricalTemperature", TuLiPa.TIMEVECTORPARSETYPES, elkey) - hist_Lday = TuLiPa.getdictvalue(value, "HistoricalDaylight", TuLiPa.TIMEVECTORPARSETYPES, elkey) - - ndays_pred = TuLiPa.getdictvalue(value, "NDaysPred", Real, elkey) - try - ndays_pred = Int(ndays_pred) - @assert ndays_pred >= 0 - catch e - error("Value for key NDaysPred must be positive integer for $elkey") - end - - basin_area = TuLiPa.getdictvalue(value, "BasinArea", Float64, elkey) - - deps = TuLiPa.Id[] - - all_ok = true - - (id, hist_P, ok) = TuLiPa.getdicttimevectorvalue(lowlevel, hist_P) - all_ok = all_ok && ok - TuLiPa._update_deps(deps, id, ok) - - (id, hist_T, ok) = TuLiPa.getdicttimevectorvalue(lowlevel, hist_T) - all_ok = all_ok && ok - TuLiPa._update_deps(deps, id, ok) - - (id, hist_Lday, ok) = TuLiPa.getdicttimevectorvalue(lowlevel, hist_Lday) - all_ok = all_ok && ok - TuLiPa._update_deps(deps, id, ok) - - if all_ok == false - return (false, deps) - end - - # TODO: Maybe make this user input in future? - updater = SimpleIfmDataUpdater() - - model_params = JLD2.load_object(model_params) - - is_nn = !isnothing(moments) - if is_nn - # convert model_params, stored with simpler data structure for stability between versions, - # into ComponentArray, which the NN-model needs - # (the simpler data structure is Vector{Tuple{Vector{Float32}, Vector{Float32}}}) - _subarray(i) = ComponentArray(weight = model_params[i][1], bias = model_params[i][2]) - _tuple(i) = (Symbol("layer_", i), _subarray(i)) - model_params = ComponentArray(NamedTuple(_tuple(i) for i in eachindex(model_params))) - # add moments, which is needed to normalize state inputs to the NN-model - moments = JLD2.load_object(moments) - model_params = (model_params, moments) - end - - data_forecast = nothing - data_obs = nothing - ndays_obs = 365 - data_forecast = nothing - ndays_forecast = 0 - - - - periodkey = TuLiPa.Id(TuLiPa.TIMEPERIOD_CONCEPT, "ScenarioTimePeriod") - period = lowlevel[periodkey] - scen_start = period["Start"] - scen_stop = period["Stop"] - - - id = TuLiPa.getobjkey(elkey) - toplevel[id] = Constructor(id, model_params, updater, basin_area, hist_P, hist_T, hist_Lday, - ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast, scen_start, scen_stop) - - return (true, deps) + error(extension_error_msg) end # --- Functions used in run_serial in connection with inflow models --- @@ -472,30 +357,7 @@ function save_ifm_Q(div_db, inflow_name, stepnr, Q) end function calculate_normalize_factor(ifm_model) - start_with_buffer = ifm_model.handler.scen_start - Day(ifm_model.handler.ndays_obs) # Add ndays_obs days buffer - days_with_buffer = Day(ifm_model.handler.scen_stop - start_with_buffer) |> Dates.value - timepoints_with_buffer = (1:days_with_buffer) - days = Dates.value(Day(ifm_model.handler.scen_stop - ifm_model.handler.scen_start)) - timepoints_start = days_with_buffer - days - - P = zeros(length(timepoints_with_buffer)) - T = zeros(length(timepoints_with_buffer)) - Lday = zeros(length(timepoints_with_buffer)) - for i in timepoints_with_buffer - start = ifm_model.handler.scen_start + Day(i - 1) - P[i] = TuLiPa.getweightedaverage(ifm_model.handler.hist_P, start, JulES.ONEDAY_MS_TIMEDELTA) - T[i] = TuLiPa.getweightedaverage(ifm_model.handler.hist_T, start, JulES.ONEDAY_MS_TIMEDELTA) - Lday[i] = TuLiPa.getweightedaverage(ifm_model.handler.hist_Lday, start, JulES.ONEDAY_MS_TIMEDELTA) - end - - itp_method = JulES.SteffenMonotonicInterpolation() - itp_P = JulES.interpolate(timepoints_with_buffer, P, itp_method) - itp_T = JulES.interpolate(timepoints_with_buffer, T, itp_method) - itp_Lday = JulES.interpolate(timepoints_with_buffer, Lday, itp_method) - Q, _ = JulES.predict(ifm_model.handler.predictor, 0, 0, itp_Lday, itp_P, itp_T, timepoints_with_buffer) - Q = Float64.(Q)[timepoints_start:end] - Q .= Q .* ifm_model.handler.m3s_per_mm - return 1 / mean(Q) + error(extension_error_msg) end """ @@ -823,14 +685,14 @@ function copy_elements_iprogtype(elements, iprogtype, ifm_names, ifm_derivedname end function initialize_NN_model() - error("Extension not loaded") + error(extension_error_msg) end function basic_bucket_incl_states(p_, itp_Lday, itp_P, itp_T, t_out) - error("Extension not loaded") + error(extension_error_msg) end function NeuralODE_M100(p, norm_S0, norm_S1, norm_P, norm_T, itp_Lday, itp_P, itp_T, t_out, ann; S_init = [0.0, 0.0]) - error("Extension not loaded") + error(extension_error_msg) end \ No newline at end of file diff --git a/src/run_jules_wrapper.jl b/src/run_jules_wrapper.jl index 824d134..b5f2b08 100644 --- a/src/run_jules_wrapper.jl +++ b/src/run_jules_wrapper.jl @@ -51,11 +51,15 @@ function load_ifm_dep() println("All packages already installed.") end end - ensure_packages(["OrdinaryDiffEq"]) + ensure_packages(["OrdinaryDiffEq", "ComponentArrays", "Interpolations", "JLD2"]) end @everywhere begin + Pkg.instantiate() Base.eval(Main, :(using OrdinaryDiffEq)) + Base.eval(Main, :(using ComponentArrays)) + Base.eval(Main, :(using Interpolations)) + Base.eval(Main, :(using JLD2)) end end