Merge pull request #25 from NVE/normfactor_bugfix

calc normfactor with ifm_model
NVE · Oct 2, 2024 · 7c2904a · 7c2904a
2 parents 5d3a94c + 437fc0a
commit 7c2904a
Show file tree

Hide file tree

Showing 4 changed files with 165 additions and 13 deletions.
diff --git a/demos/testing_normfactor.ipynb 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
+}
diff --git a/src/ifm.jl b/src/ifm.jl
@@ -62,8 +62,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
@@ -79,7 +82,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
 
@@ -260,10 +263,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
@@ -313,10 +316,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
@@ -401,9 +404,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
@@ -460,6 +471,33 @@ function save_ifm_Q(db, inflow_name, stepnr, Q)
     end
 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)
+end
+
 """
 Sequentially solve inflow models stored locally. 
 Each inflow model is solved for each scenario.
@@ -476,8 +514,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
@@ -486,7 +522,13 @@ 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) 
+                normfactors[inflow_name] = normalize_factor
+            else
+                normalize_factor = normfactors[inflow_name]
+            end
 
             u0 = estimate_u0(inflow_model, t)
 

diff --git 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
@@ -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_elements(db::LocalDB) = get_ifm_elements(get_input(db))
 
 get_cores(db::LocalDB) = get_cores(get_input(db))