Discrete Control performance improvements based on AGV model

core/src/callback.jl

@@ -195,22 +195,31 @@
 end
 
 """
-Apply the discrete control logic. Each DiscreteControl node has one or more variables associated with it,
-and these variables are a linear combination of one or more state/time derived parameters of the model.
-The values of these variables are compared against 'greater_than' thresholds. These conditions yield a 'truth_state'
-consisting of 'truth_values' for a DiscreteControl node. This truth_state then corresponds with a 'control_state',
-and nodes controlled by a DiscreteControl node have parameter values associated with this control_state.
+Apply the discrete control logic. There's somewhat of a complex structure:
+- Each DiscreteControl node can have one ore multiple compound variables it listens to
+- A compound variable is defined as a linear combination of state/time derived parameters of the model
+- Each compound variable has associated with it a sorted vector of greater_than values, which define an ordered
+    list of conditions of the form (compound variable value) => greater_than
+- Thus, to find out which conditions are true, we only need to find the largest index in the greater than values
+    such that the above condition is true
+- The truth value (true/false) of all these conditions for all variables of a DiscreteControl node are concatenated
+    (in preallocated memory) into what is called the nodes truth state. This concatenation happens in the order in which
+    the compound variables appear in discrete_control.compound_variables
+- The DiscreteControl node maps this truth state via the logic mapping to a control state, which is a string
+- The nodes that are controlled by this DiscreteControl node must have the same control state, for which they have
+    parameter values associated with that control state defined in their control_mapping
 """
 function apply_discrete_control!(u, t, integrator)::Nothing
     (; p) = integrator
     (; discrete_control) = p
-    (; node_id, node_id_unique) = discrete_control
+    (; node_id) = discrete_control
 
-    variable_idx = 1
-    n_variables = length(node_id)
-
-    # Loop over the discrete control nodes
-    for (id, truth_state) in zip(node_id_unique, discrete_control.truth_state)
+    # Loop over the discrete control nodes to determine their truth state
+    # and detect possible control state changes
+    for i in eachindex(node_id)
+        id = node_id[i]
+        truth_state = discrete_control.truth_state[i]
+        compound_variables = discrete_control.compound_variables[i]
 
         # Whether a change in truth state was detected, and thus whether
         # a change in control state is possible
@@ -221,18 +230,12 @@
         truth_value_variable_idx = 1
 
         # Loop over the variables listened to by this discrete control node
-        while variable_idx <= n_variables && node_id[variable_idx] == id
+        for compound_variable in compound_variables
 
             # Compute the value of the current variable
             value = 0.0
-            compound_variable = discrete_control.variable[variable_idx]
-            for (listen_node_id, variable, weight, look_ahead) in zip(
-                compound_variable.listen_node_id,
-                compound_variable.variable,
-                compound_variable.weight,
-                compound_variable.look_ahead,
-            )
-                value += weight * get_value(p, listen_node_id, variable, look_ahead, u, t)
+            for subvariable in compound_variable.subvariables
+                value += subvariable.weight * get_value(p, subvariable, t)
             end
 
             # The thresholds the value of this variable is being compared with
@@ -248,7 +251,7 @@
             # corresponding to the conditions on the current variable
             for truth_value_idx in
                 truth_value_variable_idx:(truth_value_variable_idx + n_greater_than - 1)
-                new_truth_state = truth_value_idx <= largest_true_index
+                new_truth_state = (truth_value_idx <= largest_true_index)
                 # If no truth state change was detected yet, check whether there is a change
                 # at this position
                 if !truth_state_change
@@ -258,7 +261,6 @@
             end
 
             truth_value_variable_idx += n_greater_than
-            variable_idx += 1
         end
 
         # If no truth state change whas detected for this node, no control
@@ -315,26 +317,20 @@
 Get a value for a condition. Currently supports getting levels from basins and flows
 from flow boundaries.
 """
-function get_value(
-    p::Parameters,
-    node_id::NodeID,
-    variable::String,
-    Δt::Float64,
-    u::AbstractVector{Float64},
-    t::Float64,
-)
+function get_value(p::Parameters, subvariable::NamedTuple, t::Float64)
     (; basin, flow_boundary, level_boundary) = p
+    (; listen_node_id, look_ahead, variable) = subvariable
 
     if variable == "level"
-        if node_id.type == NodeType.Basin
-            has_index, basin_idx = id_index(basin.node_id, node_id)
+        if listen_node_id.type == NodeType.Basin
+            has_index, basin_idx = id_index(basin.node_id, listen_node_id)
             if !has_index
-                error("Discrete control listen node $node_id does not exist.")
+                error("Discrete control listen node $listen_node_id does not exist.")
             end
-            _, level = get_area_and_level(basin, basin_idx, u[basin_idx])
-        elseif node_id.type == NodeType.LevelBoundary
-            level_boundary_idx = findsorted(level_boundary.node_id, node_id)
-            level = level_boundary.level[level_boundary_idx](t + Δt)
+            level = get_tmp(basin.current_level, 0)[basin_idx]
+        elseif listen_node_id.type == NodeType.LevelBoundary
+            level_boundary_idx = findsorted(level_boundary.node_id, listen_node_id)
+            level = level_boundary.level[level_boundary_idx](t + look_ahead)
         else
             error(
                 "Level condition node '$node_id' is neither a basin nor a level boundary.",
@@ -343,11 +339,11 @@
         value = level
 
     elseif variable == "flow_rate"
-        if node_id.type == NodeType.FlowBoundary
-            flow_boundary_idx = findsorted(flow_boundary.node_id, node_id)
-            value = flow_boundary.flow_rate[flow_boundary_idx](t + Δt)
+        if listen_node_id.type == NodeType.FlowBoundary
+            flow_boundary_idx = findsorted(flow_boundary.node_id, listen_node_id)
+            value = flow_boundary.flow_rate[flow_boundary_idx](t + look_ahead)
         else
-            error("Flow condition node $node_id is not a flow boundary.")
+            error("Flow condition node $listen_node_id is not a flow boundary.")
         end
 
     else
@@ -463,8 +459,8 @@
         discrete_control_parameter_update!(p.pid_control, node_id, control_state)
     elseif node_id.type == NodeType.LinearResistance
         discrete_control_parameter_update!(p.linear_resistance, node_id, control_state)
     elseif node_id.type == NodeType.ManningResistance
         discrete_control_parameter_update!(p.manning_resistance, node_id, control_state)
     end
     return nothing
 end

core/src/graph.jl

@@ -105,9 +105,11 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra
     if config.solver.autodiff
         flow = DiffCache(flow, chunk_sizes)
     end
+    flow_edges = EdgeMetadata[]
     graph_data = (;
         node_ids,
         edges_source,
+        flow_edges,
         flow_dict,
         flow,
         flow_prev,

core/src/parameter.jl

@@ -336,11 +336,13 @@ end
 
 """
 node_id: node ID of the FlowBoundary node
+outflow_ids: The downsteam nodes of this FlowBoundary node
 active: whether this node is active and thus contributes flow
 flow_rate: target flow rate
 """
 struct FlowBoundary <: AbstractParameterNode
     node_id::Vector{NodeID}
+    outflow_ids::Vector{Vector{NodeID}}
     active::BitVector
     flow_rate::Vector{ScalarInterpolation}
 end
@@ -467,31 +469,37 @@ struct Terminal <: AbstractParameterNode
     node_id::Vector{NodeID}
 end
 
+"""
+The data for a single compound variable
+node_id:: The ID of the DiscreteControl that listens to this variable
+subvariables: data for one single subvariable
+greater_than: the thresholds this compound variable will be
+    compared against
+"""
 struct CompoundVariable
-    listen_node_id::Vector{NodeID}
-    variable::Vector{String}
-    weight::Vector{Float64}
-    look_ahead::Vector{Float64}
+    node_id::NodeID
+    subvariables::Vector{
+        @NamedTuple{
+            listen_node_id::NodeID,
+            variable::String,
+            weight::Float64,
+            look_ahead::Float64,
+        }
+    }
     greater_than::Vector{Float64}
 end
 
 """
-node_id: node ID of the DiscreteControl node per compound variable (can contain repeats)
-listen_node_id: the IDs of the nodes being condition on per compound variable
-variable: the names of the variables in the condition per compound variable
-weight: the weight of the variables in the condition per compound variable
-look_ahead: the look ahead of variables in the condition in seconds per compound_variable
-greater_than: The threshold values per compound variable
-condition_value: The current truth value of each condition per compound_variable per greater_than
+node_id: node ID of the DiscreteControl (if it has at least one condition defined on it)
+compound_variables: The compound variables the DiscreteControl node listens to
 truth_state: Memory allocated for storing the truth state
 control_state: Dictionary: node ID => (control state, control state start)
 logic_mapping: Dictionary: (control node ID, truth state) => control state
 record: Namedtuple with discrete control information for results
 """
 struct DiscreteControl <: AbstractParameterNode
     node_id::Vector{NodeID}
-    node_id_unique::Vector{NodeID}
-    variable::Vector{CompoundVariable}
+    compound_variables::Vector{Vector{CompoundVariable}}
     # truth_state per discrete control node
     truth_state::Vector{Vector{Bool}}
     # Definition of logic
@@ -620,6 +628,12 @@ struct LevelDemand <: AbstractDemandNode
     priority::Vector{Int32}
 end
 
+"""
+node_id: node ID of the FlowDemand node
+demand_itp: The time interpolation of the demand of the node
+demand: The current demand of the node
+priority: The priority of the demand of the node
+"""
 struct FlowDemand <: AbstractDemandNode
     node_id::Vector{NodeID}
     demand_itp::Vector{ScalarInterpolation}
@@ -635,27 +649,45 @@ struct Subgrid
     level::Vector{Float64}
 end
 
+"""
+The metadata of the graph (the fields of the NamedTuple) can be accessed
+    e.g. using graph[].flow.
+node_ids: mapping subnetwork ID -> node IDs in that subnetwork
+edges_source: mapping subnetwork ID -> metadata of allocation
+    source edges in that subnetwork
+flow_edges: The metadata of all flow edges
+flow dict: mapping (source ID, destination ID) -> index in the flow vector
+    of the flow over that edge
+flow: Flow per flow edge in the order prescribed by flow_dict
+flow_prev: The flow vector of the previous timestep, used for integration
+flow_integrated: Flow integrated over time, used for mean flow computation
+    over saveat intervals
+saveat: The time interval between saves of output data (storage, flow, ...)
+"""
+const ModelGraph{T} = MetaGraph{
+    Int64,
+    DiGraph{Int64},
+    NodeID,
+    NodeMetadata,
+    EdgeMetadata,
+    @NamedTuple{
+        node_ids::Dict{Int32, Set{NodeID}},
+        edges_source::Dict{Int32, Set{EdgeMetadata}},
+        flow_edges::Vector{EdgeMetadata},
+        flow_dict::Dict{Tuple{NodeID, NodeID}, Int32},
+        flow::T,
+        flow_prev::Vector{Float64},
+        flow_integrated::Vector{Float64},
+        saveat::Float64,
+    },
+    MetaGraphsNext.var"#11#13",
+    Float64,
+} where {T}
+
 # TODO Automatically add all nodetypes here
 struct Parameters{T, C1, C2, V1, V2, V3}
     starttime::DateTime
-    graph::MetaGraph{
-        Int64,
-        DiGraph{Int64},
-        NodeID,
-        NodeMetadata,
-        EdgeMetadata,
-        @NamedTuple{
-            node_ids::Dict{Int32, Set{NodeID}},
-            edges_source::Dict{Int32, Set{EdgeMetadata}},
-            flow_dict::Dict{Tuple{NodeID, NodeID}, Int32},
-            flow::T,
-            flow_prev::Vector{Float64},
-            flow_integrated::Vector{Float64},
-            saveat::Float64,
-        },
-        MetaGraphsNext.var"#11#13",
-        Float64,
-    }
+    graph::ModelGraph{T}
     allocation::Allocation
     basin::Basin{T, C1, V1, V2, V3}
     linear_resistance::LinearResistance