From ccc90cf276321af30c189ed82740da590f705afd Mon Sep 17 00:00:00 2001 From: pertft Date: Tue, 1 Oct 2024 11:38:18 +0200 Subject: [PATCH 1/4] calc normfactor with ifm_model --- src/ifm.jl | 70 +++++++++++++++++++++++++++++++++++++++++++------ src/io.jl | 2 +- src/local_db.jl | 7 +++-- 3 files changed, 68 insertions(+), 11 deletions(-) diff --git a/src/ifm.jl b/src/ifm.jl index 943bcb5..a68c065 100644 --- a/src/ifm.jl +++ b/src/ifm.jl @@ -61,8 +61,11 @@ mutable struct TwoStateIfmHandler{P <: AbstractTwoStateIfmPredictor, prev_t::Union{TuLiPa.ProbTime, Nothing} ndays_forecast_used::Int + scen_start::Any + scen_stop::Any + function TwoStateIfmHandler(predictor, updater, basin_area, hist_P, hist_T, hist_Lday, - ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast) + ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast, scen_start, scen_stop) @assert ndays_forecast >= 0 @assert ndays_pred >= ndays_forecast @assert ndays_obs > 0 @@ -78,7 +81,7 @@ mutable struct TwoStateIfmHandler{P <: AbstractTwoStateIfmPredictor, T3 = typeof(hist_Lday) return new{P, U, T1, T2, T3}(predictor, updater, basin_area, m3s_per_mm, hist_P, hist_T, hist_Lday, ndays_pred, ndays_obs, ndays_forecast, data_pred, - data_obs, data_forecast, nothing, 0) + data_obs, data_forecast, nothing, 0, scen_start, scen_stop) end end @@ -259,10 +262,10 @@ struct TwoStateBucketIfm{H} <: AbstractInflowModel handler::H function TwoStateBucketIfm(id, model_params, updater, basin_area, hist_P, hist_T, hist_Lday, - ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast) + ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast, scen_start, scen_stop) predictor = TwoStateBucketIfmPredictor(model_params) handler = TwoStateIfmHandler(predictor, updater, basin_area, hist_P, hist_T, hist_Lday, - ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast) + ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast, scen_start, scen_stop) return new{typeof(handler)}(id, handler) end end @@ -312,10 +315,10 @@ struct TwoStateNeuralODEIfm{H} <: AbstractInflowModel handler::H function TwoStateNeuralODEIfm(id, model_params, updater, basin_area, hist_P, hist_T, hist_Lday, - ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast) + ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast, scen_start, scen_stop) predictor = TwoStateNeuralODEIfmPredictor(model_params) handler = TwoStateIfmHandler(predictor, updater, basin_area, hist_P, hist_T, hist_Lday, - ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast) + ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast, scen_start, scen_stop) return new{typeof(handler)}(id, handler) end end @@ -400,9 +403,17 @@ function common_includeTwoStateIfm!(Constructor, toplevel::Dict, lowlevel::Dict, 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) + ndays_pred, ndays_obs, ndays_forecast, data_obs, data_forecast, scen_start, scen_stop) return (true, deps) end @@ -459,6 +470,42 @@ function save_ifm_Q(db, inflow_name, stepnr, Q) end end +function calculate_normalize_factor(ifm_model) + + S0 = 0 + G0 = 0 + itp_Lday = ifm_model.handler.hist_Lday + itp_P = ifm_model.handler.hist_P + itp_T = ifm_model.handler.hist_T + + + days = Dates.value( Day(ifm_model.handler.scen_stop - ifm_model.handler.scen_start)) + + timepoints = collect((1: days)) + + P = Vector{Float64}([i for i in timepoints]) + T = Vector{Float64}([i for i in timepoints]) + Lday =Vector{Float64}([i for i in timepoints]) + for i in timepoints + start = ifm_model.handler.scen_start + Day(i - 1) + P[i] = TuLiPa.getweightedaverage(itp_P, start, JulES.ONEDAY_MS_TIMEDELTA) + T[i] = TuLiPa.getweightedaverage(itp_T, start, JulES.ONEDAY_MS_TIMEDELTA) + Lday[i] = TuLiPa.getweightedaverage(itp_Lday, start, JulES.ONEDAY_MS_TIMEDELTA) + end + + itp_method = JulES.SteffenMonotonicInterpolation() + itp_P = JulES.interpolate(timepoints, P, itp_method) + itp_T = JulES.interpolate(timepoints, T, itp_method) + itp_Lday = JulES.interpolate(timepoints, Lday, itp_method) + + res = JulES.predict(ifm_model.handler.predictor, S0, G0, itp_Lday, itp_P, itp_T, timepoints); + (Q, _) = res + Q = Float64.(Q) + Q .= Q .* ifm_model.handler.m3s_per_mm + + return 1 / sum(Q) +end + """ Sequentially solve inflow models stored locally. Each inflow model is solved for each scenario. @@ -485,7 +532,14 @@ function solve_ifm(t, stepnr) db.ifm_output[inflow_name] = Dict() end inflow_model = db.ifm[inflow_name] - normalize_factor = normfactors[inflow_name] + + + #if haskey(normfactors, inflow_name) == false + normalize_factor = calculate_normalize_factor(inflow_model) + # TODO: bruk info om hvor lang hist perioden er, start ett år før det, regn ut snittet av døgnverdiene, 1/snitet = normalizefactor + #else + #normalize_factor = normfactors[inflow_name] + #end u0 = estimate_u0(inflow_model, t) diff --git a/src/io.jl b/src/io.jl index e2a18a8..2420099 100644 --- a/src/io.jl +++ b/src/io.jl @@ -148,7 +148,7 @@ end get_iprogtype(input::DefaultJulESInput) = get(input.dataset, "iprogtype", "direct") has_ifm_results(input::DefaultJulESInput) = get_iprogtype(input) != "direct" -get_ifm_normfactors(input::DefaultJulESInput) = get(input.dataset, "ifm_normfactors", Dict{String, Float64}()) +#get_ifm_normfactors(input::DefaultJulESInput) = get(input.dataset, "ifm_normfactors", Dict{String, Float64}()) get_ifm_elements(input::DefaultJulESInput) = get(input.dataset, "ifm_elements", JulES.TuLiPa.DataElement[]) function get_ifm_names(input::DefaultJulESInput) diff --git a/src/local_db.jl b/src/local_db.jl index efba80a..431d8fb 100644 --- a/src/local_db.jl +++ b/src/local_db.jl @@ -76,6 +76,9 @@ mutable struct LocalDB ifm_output::Dict{String, Dict{ScenarioIx, Tuple{Vector{DateTime}, Vector{Float64}}}} ifm_derived::Dict{String, Dict{ScenarioIx, Tuple{Vector{DateTime}, Vector{Float64}}}} + ifm_normfactors::Dict{String, Float64} + + function LocalDB() return new( -1, @@ -115,7 +118,7 @@ mutable struct LocalDB Dict(), # ifm_output Dict(), # ifm_weighted - + Dict(), ) end end @@ -150,7 +153,7 @@ get_ifm_output(db::LocalDB) = db.ifm_output get_ifm_derived(db::LocalDB) = db.ifm_derived get_ifm_weights(db::LocalDB) = get_ifm_weights(get_input(db)) -get_ifm_normfactors(db::LocalDB) = get_ifm_normfactors(get_input(db)) +get_ifm_normfactors(db::LocalDB) = db.ifm_normfactors #get_ifm_normfactors(get_input(db)) get_ifm_elements(db::LocalDB) = get_ifm_elements(get_input(db)) get_cores(db::LocalDB) = get_cores(get_input(db)) From 1d9b13e8d0ee9d073b7ae79b9c9e9d4142d2ec92 Mon Sep 17 00:00:00 2001 From: pertft Date: Tue, 1 Oct 2024 12:35:28 +0200 Subject: [PATCH 2/4] mean instead of sum --- src/ifm.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/ifm.jl b/src/ifm.jl index a68c065..f65b679 100644 --- a/src/ifm.jl +++ b/src/ifm.jl @@ -498,12 +498,12 @@ function calculate_normalize_factor(ifm_model) itp_T = JulES.interpolate(timepoints, T, itp_method) itp_Lday = JulES.interpolate(timepoints, Lday, itp_method) - res = JulES.predict(ifm_model.handler.predictor, S0, G0, itp_Lday, itp_P, itp_T, timepoints); + res = JulES.predict(ifm_model.handler.predictor, S0, G0, itp_Lday, itp_P, itp_T, timepoints) (Q, _) = res Q = Float64.(Q) Q .= Q .* ifm_model.handler.m3s_per_mm - return 1 / sum(Q) + return 1 / mean(Q) end """ @@ -534,12 +534,12 @@ function solve_ifm(t, stepnr) inflow_model = db.ifm[inflow_name] - #if haskey(normfactors, inflow_name) == false - normalize_factor = calculate_normalize_factor(inflow_model) + if haskey(normfactors, inflow_name) == false + normalize_factor = calculate_normalize_factor(inflow_model) # TODO: bruk info om hvor lang hist perioden er, start ett år før det, regn ut snittet av døgnverdiene, 1/snitet = normalizefactor - #else - #normalize_factor = normfactors[inflow_name] - #end + else + normalize_factor = normfactors[inflow_name] + end u0 = estimate_u0(inflow_model, t) From 741feb03d37d6e2cba056eccd97058c7c4088d47 Mon Sep 17 00:00:00 2001 From: pertft Date: Tue, 1 Oct 2024 14:46:20 +0200 Subject: [PATCH 3/4] adding notebook for testing --- demos/testing_normfactor.ipynb | 108 +++++++++++++++++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 demos/testing_normfactor.ipynb diff --git a/demos/testing_normfactor.ipynb b/demos/testing_normfactor.ipynb new file mode 100644 index 0000000..92bf308 --- /dev/null +++ b/demos/testing_normfactor.ipynb @@ -0,0 +1,108 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "51a73813", + "metadata": {}, + "outputs": [], + "source": [ + "using Pkg\n", + "using JulES\n", + "using TuLiPa\n", + "using Dates, Statistics" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "65046806", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.06619134772002751" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "TuLiPa.INCLUDEELEMENT[TuLiPa.TypeKey(JulES.ABSTRACT_INFLOW_MODEL, \"TwoStateBucketIfm\")] = JulES.includeTwoStateBucketIfm!\n", + "\n", + "inflowpath = joinpath(dirname(dirname(pwd())), \"input\", \"static_input\") # path to ifm_elements data\n", + "\n", + "ifm_elements = JulES.TuLiPa.getelements(JulES.JSON.parsefile(joinpath(inflowpath, \"ifm_elements.json\")), inflowpath)\n", + "\n", + "TuLiPa.addscenariotimeperiod!(ifm_elements, \"ScenarioTimePeriod\", \n", + " TuLiPa.getisoyearstart(1991), \n", + " TuLiPa.getisoyearstart(2000))\n", + "\n", + "modelobjects = getmodelobjects(ifm_elements)\n", + "ifm_model = modelobjects[TuLiPa.Id(\"AbstractInflowModel\", \"127.13\")]\n", + "\n", + "S0 = 0\n", + "G0 = 0\n", + "itp_Lday = ifm_model.handler.hist_Lday\n", + "itp_P = ifm_model.handler.hist_P\n", + "itp_T = ifm_model.handler.hist_T\n", + "\n", + "start_with_buffer = ifm_model.handler.scen_start - Day(365)\n", + "days = Dates.value( Day(ifm_model.handler.scen_stop - start_with_buffer))\n", + "real_days = Dates.value( Day(ifm_model.handler.scen_stop - ifm_model.handler.scen_start)) \n", + "timepoints_start = days - real_days \n", + "timepoints = collect((1: days))\n", + "\n", + "P = Vector{Float64}([i for i in timepoints])\n", + "T = Vector{Float64}([i for i in timepoints])\n", + "Lday =Vector{Float64}([i for i in timepoints])\n", + "for i in timepoints\n", + " start = ifm_model.handler.scen_start + Day(i - 1)\n", + " P[i] = TuLiPa.getweightedaverage(itp_P, start, JulES.ONEDAY_MS_TIMEDELTA)\n", + " T[i] = TuLiPa.getweightedaverage(itp_T, start, JulES.ONEDAY_MS_TIMEDELTA)\n", + " Lday[i] = TuLiPa.getweightedaverage(itp_Lday, start, JulES.ONEDAY_MS_TIMEDELTA)\n", + "end\n", + "\n", + "itp_method = JulES.SteffenMonotonicInterpolation()\n", + "itp_P = JulES.interpolate(timepoints, P, itp_method)\n", + "itp_T = JulES.interpolate(timepoints, T, itp_method)\n", + "itp_Lday = JulES.interpolate(timepoints, Lday, itp_method)\n", + "\n", + "res = JulES.predict(ifm_model.handler.predictor, S0, G0, itp_Lday, itp_P, itp_T, timepoints);\n", + "\n", + "(Q, _) = res\n", + "Q = Float64.(Q)\n", + "Q = Q[timepoints_start:end]\n", + "Q .= Q .* ifm_model.handler.m3s_per_mm\n", + "\n", + "1 / mean(Q)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8df5239f", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.9.2", + "language": "julia", + "name": "julia-1.9" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 437fc0a64171798644051c97fe1230daf1302e9b Mon Sep 17 00:00:00 2001 From: pertft Date: Wed, 2 Oct 2024 11:01:21 +0200 Subject: [PATCH 4/4] adding buffer time to estimate startvalues for normfactor --- src/ifm.jl | 52 +++++++++++++++++++------------------------------ src/io.jl | 1 - src/local_db.jl | 2 +- 3 files changed, 21 insertions(+), 34 deletions(-) diff --git a/src/ifm.jl b/src/ifm.jl index f65b679..11e2b47 100644 --- a/src/ifm.jl +++ b/src/ifm.jl @@ -471,38 +471,29 @@ function save_ifm_Q(db, inflow_name, stepnr, Q) end function calculate_normalize_factor(ifm_model) - - S0 = 0 - G0 = 0 - itp_Lday = ifm_model.handler.hist_Lday - itp_P = ifm_model.handler.hist_P - itp_T = ifm_model.handler.hist_T - - - days = Dates.value( Day(ifm_model.handler.scen_stop - ifm_model.handler.scen_start)) - - timepoints = collect((1: days)) - - P = Vector{Float64}([i for i in timepoints]) - T = Vector{Float64}([i for i in timepoints]) - Lday =Vector{Float64}([i for i in timepoints]) - for i in timepoints + 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(itp_P, start, JulES.ONEDAY_MS_TIMEDELTA) - T[i] = TuLiPa.getweightedaverage(itp_T, start, JulES.ONEDAY_MS_TIMEDELTA) - Lday[i] = TuLiPa.getweightedaverage(itp_Lday, start, JulES.ONEDAY_MS_TIMEDELTA) + 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, P, itp_method) - itp_T = JulES.interpolate(timepoints, T, itp_method) - itp_Lday = JulES.interpolate(timepoints, Lday, itp_method) - - res = JulES.predict(ifm_model.handler.predictor, S0, G0, itp_Lday, itp_P, itp_T, timepoints) - (Q, _) = res - Q = Float64.(Q) + 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) end @@ -522,8 +513,6 @@ function solve_ifm(t, stepnr) steplen_f = float(steplen_ms.value) ifmstep_f = float(ifmstep_ms.value) remainder_f = float(remainder_ms.value) - - normfactors = get_ifm_normfactors(db) scenarios = get_scenarios(db.scenmod_sim) for (inflow_name, core) in db.dist_ifm @@ -532,11 +521,10 @@ function solve_ifm(t, stepnr) db.ifm_output[inflow_name] = Dict() end inflow_model = db.ifm[inflow_name] - if haskey(normfactors, inflow_name) == false normalize_factor = calculate_normalize_factor(inflow_model) - # TODO: bruk info om hvor lang hist perioden er, start ett år før det, regn ut snittet av døgnverdiene, 1/snitet = normalizefactor + normfactors[inflow_name] = normalize_factor else normalize_factor = normfactors[inflow_name] end diff --git a/src/io.jl b/src/io.jl index 2420099..0915d4a 100644 --- a/src/io.jl +++ b/src/io.jl @@ -148,7 +148,6 @@ end get_iprogtype(input::DefaultJulESInput) = get(input.dataset, "iprogtype", "direct") has_ifm_results(input::DefaultJulESInput) = get_iprogtype(input) != "direct" -#get_ifm_normfactors(input::DefaultJulESInput) = get(input.dataset, "ifm_normfactors", Dict{String, Float64}()) get_ifm_elements(input::DefaultJulESInput) = get(input.dataset, "ifm_elements", JulES.TuLiPa.DataElement[]) function get_ifm_names(input::DefaultJulESInput) diff --git a/src/local_db.jl b/src/local_db.jl index 431d8fb..e2e233f 100644 --- a/src/local_db.jl +++ b/src/local_db.jl @@ -153,7 +153,7 @@ get_ifm_output(db::LocalDB) = db.ifm_output get_ifm_derived(db::LocalDB) = db.ifm_derived get_ifm_weights(db::LocalDB) = get_ifm_weights(get_input(db)) -get_ifm_normfactors(db::LocalDB) = db.ifm_normfactors #get_ifm_normfactors(get_input(db)) +get_ifm_normfactors(db::LocalDB) = db.ifm_normfactors get_ifm_elements(db::LocalDB) = get_ifm_elements(get_input(db)) get_cores(db::LocalDB) = get_cores(get_input(db))