Ifm_v1

src/abstract_types.jl

@@ -8,9 +8,11 @@ abstract type AbstractScenarioModellingMethod end
 Interface:
     TuLiPa.assemble!(m::AbstractInflowModel)
     TuLiPa.getid(m::AbstractInflowModel) -> TuLiPa.Id
-    numstates(m::AbstractInflowModel) -> Int  ( == length(u0))
     estimate_u0(m::AbstractInflowModel, t::ProbTime) -> u0::Vector{Float64}
     predict(m::AbstractInflowModel, u0::::Vector{Float64}, t::ProbTime) -> Q::Vector{Float64}
+    get_numstates(m::AbstractInflowModel) -> Int  ( == length(u0))
+    get_statename(::AbstractInflowModel, i::Int) -> String
+    get_basin_area_m2(::AbstractInflowModel) -> Float64
 """
 abstract type AbstractInflowModel end
 const ABSTRACT_INFLOW_MODEL = "AbstractInflowModel"

src/ifm.jl

@@ -1,4 +1,14 @@
+"""
+Implementation of inflow models (ifm) TwoStateBucketIfm and TwoStateNeuralODEIfm
+and integration with JulES
+"""
+
 # TODO: Illustrate how data_obs, data_pred and data_forecast are updated over time
+# TODO: add data element for observations for an ifm (and extend interface with set_observations)
+# TODO: add data element for forecast for an ifm (and extend interface with set_forecast)
+# TODO: Use Float32 instead of Float64 in historical time vectors?
+
+const ONEDAY_MS_TIMEDELTA = TuLiPa.MsTimeDelta(Day(1))
 
 struct TwoStateIfmData
     P::Vector{Float32}   
@@ -25,8 +35,6 @@ end
 
 getndays(x::TwoStateIfmData) = length(x.P)
 
-const ONEDAY_MS_TIMEDELTA = TuLiPa.MsTimeDelta(Day(1))
-
 mutable struct TwoStateIfmHandler{P <: AbstractTwoStateIfmPredictor, 
                                 U <: AbstractTwoStateIfmDataUpdater, 
                                 T1 <: TuLiPa.TimeVector, T2 <: TuLiPa.TimeVector, T3 <: TuLiPa.TimeVector}
@@ -73,6 +81,19 @@ mutable struct TwoStateIfmHandler{P <: AbstractTwoStateIfmPredictor,
     end
 end
 
+get_numstates(m::TwoStateIfmHandler) = 2
+get_basin_area_m2(m::TwoStateIfmHandler) = m.basin_area
+
+function get_statename(::TwoStateIfmHandler, i::Int)
+    if i == 1
+        return "snow"
+    elseif i == 2
+        return "ground"
+    else
+        error("Only state ix 1 and 2 supported")
+    end
+end
+
 # TODO: Use @inbounds in for loops
 
 function _initial_data_obs_update(m::TwoStateIfmHandler, t::TuLiPa.ProbTime)
@@ -185,7 +206,7 @@ function predict(m::TwoStateIfmHandler, u0::Vector{Float64}, t::TuLiPa.ProbTime)
     itp_Lday = interpolate(m.data_pred.timepoints, m.data_pred.Lday, itp_method)
 
     (S0, G0) = u0
-    (Q, OED_sol) = predict(m.predictor, S0, G0, itp_Lday, itp_P, itp_T, m.data_pred.timepoints)
+    (Q, __) = predict(m.predictor, S0, G0, itp_Lday, itp_P, itp_T, m.data_pred.timepoints)
 
     Q = Float64.(Q)
 
@@ -194,11 +215,10 @@ function predict(m::TwoStateIfmHandler, u0::Vector{Float64}, t::TuLiPa.ProbTime)
     return Q
 end
 
-
 struct SimpleIfmDataUpdater <: AbstractTwoStateIfmDataUpdater
 end
 
-function update_prediction_data(m::TwoStateIfmHandler, updater::SimpleIfmDataUpdater, t::TuLiPa.ProbTime)
+function update_prediction_data(m::TwoStateIfmHandler, ::SimpleIfmDataUpdater, t::TuLiPa.ProbTime)
     if m.ndays_forecast > 0
         # use forecast for P and T if available
         ndays_before_estimate_u0_call = m.ndays_forecast + m.ndays_forecast_used
@@ -247,9 +267,11 @@ end
 
 TuLiPa.getid(m::TwoStateBucketIfm) = m.id
 TuLiPa.assemble!(m::TwoStateBucketIfm) = true
-numstates(::TwoStateBucketIfm) = 2
 estimate_u0(m::TwoStateBucketIfm, t::TuLiPa.ProbTime) = estimate_u0(m.handler, t)
 predict(m::TwoStateBucketIfm, u0::Vector{Float64}, t::TuLiPa.ProbTime) = predict(m.handler, u0, t)
+get_numstates(m::TwoStateBucketIfm) =  get_numstates(m.handler)
+get_statename(m::TwoStateBucketIfm, i::Int) = get_statename(m.handler, i)
+get_basin_area_m2(m::TwoStateBucketIfm) = get_basin_area_m2(m.handler)
 
 function includeTwoStateBucketIfm!(toplevel::Dict, lowlevel::Dict, elkey::TuLiPa.ElementKey, value::Dict)
     common_includeTwoStateIfm!(TwoStateBucketIfm, toplevel, lowlevel, elkey, value)
@@ -298,9 +320,11 @@ end
 
 TuLiPa.getid(m::TwoStateNeuralODEIfm) = m.id
 TuLiPa.assemble!(m::TwoStateNeuralODEIfm) = true
-numstates(::TwoStateNeuralODEIfm) = 2
 estimate_u0(m::TwoStateNeuralODEIfm, t::TuLiPa.ProbTime) = estimate_u0(m.handler, t)
 predict(m::TwoStateNeuralODEIfm, u0::Vector{Float64}, t::TuLiPa.ProbTime) = predict(m.handler, u0, t)
+get_numstates(m::TwoStateNeuralODEIfm) =  get_numstates(m.handler)
+get_statename(m::TwoStateNeuralODEIfm, i::Int) = get_statename(m.handler, i)
+get_basin_area_m2(m::TwoStateNeuralODEIfm) = get_basin_area_m2(m.handler)
 
 function includeTwoStateNeuralODEIfm!(toplevel::Dict, lowlevel::Dict, elkey::TuLiPa.ElementKey, value::Dict)
     common_includeTwoStateIfm!(TwoStateNeuralODEIfm, toplevel, lowlevel, elkey, value)
@@ -309,11 +333,6 @@ end
 function common_includeTwoStateIfm!(Constructor, toplevel::Dict, lowlevel::Dict, elkey::TuLiPa.ElementKey, value::Dict)
     TuLiPa.checkkey(toplevel, elkey)
 
-    OBSERVED_PERCIPITATION = "ObservedPercipitation"
-    OBSERVED_TEMPERATURE = "ObservedTemperature"
-    FORECASTED_PERCIPITATION = "ForecastedPercipitation"
-    FORECASTED_TEMPERATURE = "ForecastedTemperature"
-
     model_params = TuLiPa.getdictvalue(value, "ModelParams", String, elkey)
 
     moments = nothing
@@ -335,36 +354,7 @@ function common_includeTwoStateIfm!(Constructor, toplevel::Dict, lowlevel::Dict,
 
     basin_area = TuLiPa.getdictvalue(value, "BasinArea", Float64, elkey)
 
-    if haskey(value, OBSERVED_PERCIPITATION)
-        obs_P = TuLiPa.getdictvalue(value, OBSERVED_PERCIPITATION,   Vector{Real}, elkey)
-    else
-        obs_P = nothing
-    end
-    if haskey(value, OBSERVED_TEMPERATURE)
-        obs_T = TuLiPa.getdictvalue(value, OBSERVED_TEMPERATURE,   Vector{Real}, elkey)
-    else
-        obs_T = nothing
-    end
-
-    if haskey(value, FORECASTED_PERCIPITATION)
-        forecast_P = TuLiPa.getdictvalue(value, FORECASTED_PERCIPITATION,   Vector{Real}, elkey)
-    else
-        forecast_P = nothing
-    end
-    if haskey(value, FORECASTED_TEMPERATURE)
-        forecast_T = TuLiPa.getdictvalue(value, FORECASTED_TEMPERATURE,   Vector{Real}, elkey)
-    else
-        forecast_T = nothing
-    end
-
-    isnothing(obs_P) && (!isnothing(obs_T)) && error("Missing $OBSERVED_PERCIPITATION for $elkey")
-    isnothing(obs_T) && (!isnothing(obs_P)) && error("Missing $OBSERVED_TEMPERATURE for $elkey")
-
-    isnothing(obs_P) && (!isnothing(obs_T)) && error("Missing $FORECASTED_PERCIPITATION for $elkey")
-    isnothing(obs_T) && (!isnothing(obs_P)) && error("Missing $FORECASTED_TEMPERATURE for $elkey")
-
     deps = TuLiPa.Id[]
-    # errs = String[]
 
     all_ok = true
 
@@ -382,25 +372,6 @@ function common_includeTwoStateIfm!(Constructor, toplevel::Dict, lowlevel::Dict,
 
     if all_ok == false
         return (false, deps)
-        # return (false, deps, errs)
-    end
-
-    if !(isnothing(obs_P) || isnothing(obs_T))
-        dummy_Lday = zeros(eltype(obs_P), length(obs_P))
-        data_obs = TwoStateIfmData(obs_P, obs_T, dummy_Lday)
-        ndays_obs == getndays(data_obs)
-    else
-        data_obs = nothing
-        ndays_obs = 365
-    end
-
-    if !(isnothing(forecast_P) || isnothing(forecast_T))
-        dummy_Lday = zeros(eltype(forecast_P), length(forecast_P))
-        data_forecast = TwoStateIfmData(forecast_P, forecast_T, dummy_Lday)
-        ndays_forecast = getndays(data_forecast)
-    else
-        data_forecast = nothing
-        ndays_forecast = 0
     end
 
     # TODO: Maybe make this user input in future?
@@ -413,19 +384,24 @@ function common_includeTwoStateIfm!(Constructor, toplevel::Dict, lowlevel::Dict,
         model_params = (model_params, moments)
     end
 
+    data_forecast = nothing
+    data_obs = nothing
+    ndays_obs = 365
+    data_forecast = nothing
+    ndays_forecast = 0
+
     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)
 
     return (true, deps)
-    # return (true, deps, errs)
 end
 
-
 # --- Functions used in run_serial in connection with inflow models ---
 
 const IFM_DB_STATE_KEY = "ifm_step_u0"
 const IFM_DB_FLOW_KEY = "ifm_step_Q"
+const IFM_STATION_ID_KEY = "ifm_station_id"
 
 """
 Create inflow models and store some of them locally according to db.dist_ifm
@@ -439,11 +415,7 @@ function create_ifm()
     db.div[IFM_DB_FLOW_KEY] = Dict{String, Tuple{Int, Float64}}()
 
     elements = get_ifm_elements(db)
-    t0 = time()
     modelobjects = TuLiPa.getmodelobjects(elements)
-    t1 = time()
-    sec = t1 - t0
-    println("getmodelobjects(ifm_elements) took $sec")
     for (inflow_name, core) in db.dist_ifm
         if core == db.core
             id = TuLiPa.Id(ABSTRACT_INFLOW_MODEL, inflow_name)
@@ -481,6 +453,7 @@ Each inflow model is solved for each scenario.
 function solve_ifm(t, stepnr)
     db = get_local_db()
 
+    # setup utility variables to save results below
     steplen_ms = Millisecond(get_steplength(db.input))
     ifmstep_ms = TuLiPa.getduration(ONEDAY_MS_TIMEDELTA)
     @assert steplen_ms >= ifmstep_ms
@@ -500,14 +473,14 @@ function solve_ifm(t, stepnr)
             end
             inflow_model = db.ifm[inflow_name]
             normalize_factor = normfactors[inflow_name]
+
             u0 = estimate_u0(inflow_model, t)
 
-            # save in familiar unit
-            # TODO: extend interface with get_basin_area(::AbstractInflowModel) ?
-            u0_mm3 = u0 .* inflow_model.handler.basin_area ./ 1000.0 
+            # save in familiar unit mm3 (u0 in mm converted to m in calculation)
+            u0_mm3 = (u0 ./ 1000.0) .* get_basin_area_m2(inflow_model) ./ 1e6 
             save_ifm_u0(db, inflow_name, stepnr, u0_mm3)
 
-            # predict mean Q for over clearing period and store result
+            # predict mean Q for clearing period and store result
             # can be used to measure goodness of ifm model
             Q = predict(inflow_model, u0, t)
             @assert nifmsteps <= length(Q)
@@ -670,7 +643,7 @@ end
 
 
 """
-ModeledInflowParam is actually a stateful PrognosisSeriesParam under-the-hood,
+ModeledInflowParam is actually a stateful TuLiPa.PrognosisSeriesParam under-the-hood,
 but where the prognosis part is created and managed by JulES and not 
 given as user input
 """
@@ -704,9 +677,7 @@ function includeModeledInflowParam!(::Dict, lowlevel::Dict, elkey::TuLiPa.Elemen
     # This way, model objects holding reference to such Param, will use updated 
     # prognosis from db when called upon by update!(prob, t)
     db = get_local_db()
-    replacemap = get_ifm_replacemap(db.input)
-    haskey(replacemap, elkey.instancename) || error("Instance name not found in replacemap for $elkey")
-    inflow_name = replacemap[elkey.instancename]
+    inflow_name = value[IFM_STATION_ID_KEY]
     ifm_weights = get_ifm_weights(db)
     ifm_names = get_ifm_names(db.input)
     if haskey(ifm_weights, inflow_name)
@@ -732,67 +703,49 @@ function includeModeledInflowParam!(::Dict, lowlevel::Dict, elkey::TuLiPa.Elemen
     return (true, deps)
 end
 
-"""
-MixedInflowParam...
-"""
-# TODO: Complete this
-# struct MixedInflowParam{P1, P2} <: Param
-#     directparam::P1
-#     ifmparam::P2
-#     ndays::Int
-# end
-# function includeMixedInflowParam!(::Dict, lowlevel::Dict, elkey::TuLiPa.ElementKey, value::Dict)
-# end
-# TuLiPa.INCLUDEELEMENT[TuLiPa.TypeKey(TuLiPa.PARAM_CONCEPT, "MixedInflowParam")] = includeMixedInflowParam!
+# TODO: Maybe support MixedInflowParam to use external forecast first year and ifm for rest?
 
 """
 Used by JulES in appropriate places to embed scenix info 
-into data elements of type ModeledInflowParam or MixedInflowParam
+into data elements of type ModeledInflowParam
 """
-function add_scenix_to_InflowParam(elements, scenix)
+function add_scenix_to_InflowParam(elements, scenix::Int)
     for e in elements
-        if e.typename == "ModeledInflowParam" || e.typename == "MixedInflowParam"
+        if e.typename == "ModeledInflowParam"
             e.value["ScenarioIndex"] = scenix
         end
     end
 end
 
 """
-Return copy of elements with replacement of PrognosisSeriesParam 
-in ifm_replacemap in accordance with value of iprogtype
+Return copy of elements with replacement of Param with ModeledInflowParam 
+if Param is embedded with a known station_id behind IFM_STATION_ID_KEY
+The original Param must have Level and Profile keys, but this is how we 
+normally set up inflow parameters anyway, so this should be ok.
 """
-function copy_elements_iprogtype(elements, iprogtype, ifm_replacemap)
+function copy_elements_iprogtype(elements, iprogtype::String, ifm_names::Vector{String})
     if iprogtype == "ifm"
         elements1 = TuLiPa.DataElement[]
         for e in elements
-            if e.typename == "PrognosisSeriesParam" && haskey(ifm_replacemap, e.instancename)
-                new_e = TuLiPa.DataElement(e.conceptname, "ModeledInflowParam", e.instancename,
-                    Dict("Level" => e.value["Level"], "HistoricalProfile" => e.value["Profile"]))
-                push!(elements1, new_e)
-            else
-                push!(elements1, e)
-            end
-        end
-
-    elseif startswith(iprogtype, "mix")
-        # TODO: Validate ndays > 0 in constructor of DefaultJulESInput
-        ndays = parse(Int, iprogtype[4:end])
-        elements1 = TuLiPa.DataElement[]
-        for e in elements
-            if e.typename == "PrognosisSeriesParam" && haskey(ifm_replacemap, e.instancename)
-                new_value = copy(e.value::Dict)
-                new_value["ndays"] = ndays
-                new_e = TuLiPa.DataElement(e.conceptname, "MixedInflowParam", e.instancename, new_value)
-                push!(elements1, new_e)
-            else
-                push!(elements1, e)
+            maybe_new_element = e
+            if e.conceptname == TuLiPa.PARAM_CONCEPT
+                if e.value isa Dict
+                    if haskey(e.value, IFM_STATION_ID_KEY)
+                        station_id = e.value[IFM_STATION_ID_KEY]
+                        if station_id in ifm_names
+                            maybe_new_element = TuLiPa.DataElement(e.conceptname, "ModeledInflowParam", e.instancename,
+                                Dict("Level" => e.value["Level"], "HistoricalProfile" => e.value["Profile"], 
+                                    IFM_STATION_ID_KEY => station_id))
+                        end
+                    end
+                end
             end
+            push!(elements1, maybe_new_element)
         end
-
     else
         @assert iprogtype == "direct"
         elements1 = copy(elements)
     end
-
+    @assert length(elements) == length(elements1)
     return elements1
 end