From a4ca4365b0dbe31df268a2c0671d7b2ef57e4d02 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sun, 15 Dec 2024 22:19:49 +1100 Subject: [PATCH] Adjust calibration methods Changes to support fitting gauges that depend on downstream behaviour for performance assessment. --- src/calibration.jl | 149 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 132 insertions(+), 17 deletions(-) diff --git a/src/calibration.jl b/src/calibration.jl index 117eaab..85cc599 100644 --- a/src/calibration.jl +++ b/src/calibration.jl @@ -1,22 +1,26 @@ -# Ensure dependent data and packages are available -using Distributed, BlackBoxOptim, Serialization +using Serialization +using BlackBoxOptim +using Distributed -"""Calibrate specified node in network.""" -function obj_func(params, climate::Climate, sn::StreamfallNetwork, v_id::Int, calib_data::Array; - metric::Function, inflow=nothing, extraction=nothing, exchange=nothing) +"""Calibrate the specified node in network.""" +function obj_func( + params, climate::Climate, sn::StreamfallNetwork, v_id::Int, calib_data::Array; + metric::Function, inflow=nothing, extraction=nothing, exchange=nothing +) return obj_func(params, climate, sn[v_id], calib_data; metric=metric, inflow=inflow, extraction=extraction, exchange=exchange) end """Calibrate current node.""" -function obj_func(params, climate::Climate, node::NetworkNode, calib_data::Array; - metric::Function, inflow=nothing, extraction=nothing, exchange=nothing) +function obj_func( + params, climate::Climate, node::NetworkNode, calib_data::Array; + metric::Function, inflow=nothing, extraction=nothing, exchange=nothing +) update_params!(node, params...) run_node!(node, climate; inflow=inflow, extraction=extraction, exchange=exchange) - n_data = node.outflow - score = metric(calib_data, n_data) + score = metric(calib_data, node.outflow) # Reset to clear stored values reset!(node) @@ -25,6 +29,82 @@ function obj_func(params, climate::Climate, node::NetworkNode, calib_data::Array end +""" + dependent_obj_func( + params, climate::Climate, this_node::NetworkNode, next_node::NetworkNode, calib_data::DataFrame; + metric::Function, weighting=0.5, inflow=nothing, extraction=nothing, exchange=nothing + ) + +Objective function which considers performance of next node and current node. +The weighting factor places equal emphasis on both nodes by default (0.5). +The weighting value places a \$x\$ weight on the current node, and \$1 - x\$ on the next +node. + +""" +function dependent_obj_func( + params, climate::Climate, this_node::NetworkNode, next_node::NetworkNode, calib_data::DataFrame; + metric::Function, weighting=0.5, inflow=nothing, extraction=nothing, exchange=nothing +) + update_params!(this_node, params...) + + # Run dependent nodes + run_node!( + this_node, climate; + inflow=inflow, extraction=extraction, exchange=exchange + ) + run_node!( + next_node, climate; + inflow=this_node.outflow, extraction=extraction, exchange=exchange + ) + + # Alias data as necessary + if typeof(next_node) <: DamNode + # Calibrate against both outflow and dam levels + sim_data = next_node.level + obs_data = calib_data[:, next_node.name] + + if weighting == 0.0 + # All emphasis is on dam levels + score = metric(obs_data, sim_data) + elseif weighting < 1.0 + # Use a mix + dam_score = metric(obs_data, sim_data) + + sim_data = this_node.outflow + obs_data = calib_data[:, this_node.name] + + flow_score = metric(obs_data, sim_data) + + # Use weighted average of the two + score = (flow_score * weighting) + (dam_score * (1.0 - weighting)) + else + # Use just outflows + sim_data = this_node.outflow + obs_data = calib_data[:, this_node.name] + + score = metric(obs_data, sim_data) + end + elseif typeof(this_node) <: DamNode + # Calibrate against dam levels only + sim_data = this_node.level + obs_data = calib_data[:, this_node.name] + + score = metric(obs_data, sim_data) + else + # Calibrate against outflows only + sim_data = this_node.outflow + obs_data = calib_data[:, this_node.name] + + score = metric(obs_data, sim_data) + end + + reset!(this_node) + reset!(next_node) + + return score +end + + """ calibrate!(sn, v_id, climate, calib_data, metric::Function=Streamfall.RMSE; @@ -42,7 +122,7 @@ Calibrate a given node, recursing upstream, using the BlackBoxOptim package. function calibrate!(sn::StreamfallNetwork, v_id::Int64, climate::Climate, calib_data::DataFrame; metric::Function=Streamfall.RMSE, kwargs...) - calibrate!(sn, v_id, climate, calib_data[:, sn[v_id].name]; metric=metric, kwargs...) + calibrate!(sn, v_id, climate, calib_data; metric=metric, kwargs...) end @@ -57,10 +137,10 @@ Calibrate a given node, recursing upstream, using the BlackBoxOptim package. - `sn::StreamfallNetwork` : Network - `v_id::Int` : node identifier - `climate::Climate` : Climate data -- `calib_data::Array` : calibration data for target node by its id +- `calib_data::DataFrame` : calibration data for target node by its id - `metric::Function` : Optimization function to use. Defaults to RMSE. """ -function calibrate!(sn::StreamfallNetwork, v_id::Int64, climate::Climate, calib_data::Array; +function calibrate!(sn::StreamfallNetwork, v_id::Int64, climate::Climate, calib_data::DataFrame; metric::Function=Streamfall.RMSE, kwargs...) # Set defaults as necessary @@ -81,9 +161,31 @@ function calibrate!(sn::StreamfallNetwork, v_id::Int64, climate::Climate, calib_ end this_node = sn[v_id] + next_node = sn[first(outlets(sn, this_node.name))] + + extraction = get(kwargs, :extraction, nothing) + exchange = get(kwargs, :exchange, nothing) # Create context-specific optimization function - opt_func = x -> obj_func(x, climate, sn, v_id, calib_data; metric=metric) + if typeof(next_node) <: DamNode + # If the next node represents a dam, attempt to calibrate considering the outflows + # from the current node and the Dam Levels of the next node. + weighting = get(kwargs, :weighting, 0.5) + opt_func = x -> next_node.obj_func( + x, climate, this_node, next_node, calib_data; + metric=metric, extraction=extraction, exchange=exchange, weighting=weighting + ) + elseif typeof(this_node) <: DamNode + opt_func = x -> obj_func( + x, climate, this_node, calib_data[:, this_node.name]; + metric=metric, extraction=extraction, exchange=exchange + ) + else + opt_func = x -> this_node.obj_func( + x, climate, this_node, calib_data[:, this_node.name]; + metric=metric, extraction=extraction, exchange=exchange + ) + end # Get node parameters (default values and bounds) param_names, x0, param_bounds = param_info(this_node; with_level=false) @@ -128,15 +230,28 @@ function calibrate!(node::NetworkNode, climate::Climate, calib_data::Array; ) kwargs = merge(defaults, kwargs) - # reset!(node) + next_node = sn[first(outlets(sn, node.name))] + + extraction = get(kwargs, :extraction, nothing) + exchange = get(kwargs, :exchange, nothing) # Create context-specific optimization function - opt_func = x -> obj_func(x, climate, node, calib_data; metric=metric) + if next_node <: DamNode + # If the next node represents a dam, attempt to calibrate considering the outflows + # from the current node and the Dam Levels of the next node. + opt_func = x -> next_node.obj_func(x, climate, node, next_node, calib_data; metric=metric, extraction=extraction, exchange=exchange) + else + opt_func = x -> this_node.obj_func(x, climate, node, calib_data; metric=metric, extraction=extraction, exchange=exchange) + end # Get node parameters (default values and bounds) param_names, x0, param_bounds = param_info(node; with_level=false) - opt = bbsetup(opt_func; SearchRange=param_bounds, - kwargs...) + opt = bbsetup( + opt_func; + SearchRange=param_bounds, + x0=x0, + kwargs... + ) res = bboptimize(opt)