feat: Initial commit of DerivedMetric (#68)

A derived metric is updated after solve and is for values
that could be dependent on state AND paramter values. These need
updating and storing for use between time-steps.

Adds a test for piecewise storage model that shows how the volume
is updated at the end of time-step, and the derived proportional
volume is calculated and retained for use in the next time-step
despite the max volumes being updated.

Fixes #63.

pywr-core/src/derived_metric.rs

@@ -0,0 +1,96 @@
+use crate::aggregated_storage_node::AggregatedStorageNodeIndex;
+use crate::model::Model;
+use crate::node::NodeIndex;
+use crate::state::State;
+use crate::timestep::Timestep;
+use crate::virtual_storage::VirtualStorageIndex;
+use crate::PywrError;
+use std::fmt;
+use std::fmt::{Display, Formatter};
+use std::ops::Deref;
+
+#[derive(Copy, Clone, Ord, PartialOrd, Eq, PartialEq, Debug)]
+pub struct DerivedMetricIndex(usize);
+
+impl Deref for DerivedMetricIndex {
+    type Target = usize;
+
+    fn deref(&self) -> &Self::Target {
+        &self.0
+    }
+}
+
+impl DerivedMetricIndex {
+    pub fn new(idx: usize) -> Self {
+        Self(idx)
+    }
+}
+
+impl Display for DerivedMetricIndex {
+    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
+        write!(f, "{}", self.0)
+    }
+}
+
+/// Derived metrics are updated after the model is solved.
+///
+/// These metrics are "derived" from node states (e.g. volume, flow) and must be updated
+/// after those states have been updated. This should happen after the model is solved. The values
+/// are then available in this state for the next time-step.
+#[derive(Clone, Debug, PartialEq)]
+pub enum DerivedMetric {
+    NodeInFlowDeficit(NodeIndex),
+    NodeProportionalVolume(NodeIndex),
+    AggregatedNodeProportionalVolume(AggregatedStorageNodeIndex),
+    VirtualStorageProportionalVolume(VirtualStorageIndex),
+}
+
+impl DerivedMetric {
+    pub fn before(&self, timestep: &Timestep, model: &Model, state: &State) -> Result<Option<f64>, PywrError> {
+        // On the first time-step set the initial value
+        if timestep.is_first() {
+            self.compute(model, state).map(|v| Some(v))
+        } else {
+            Ok(None)
+        }
+    }
+
+    pub fn compute(&self, model: &Model, state: &State) -> Result<f64, PywrError> {
+        match self {
+            Self::NodeProportionalVolume(idx) => {
+                let max_volume = model.get_node(idx)?.get_current_max_volume(model, state)?;
+                Ok(state
+                    .get_network_state()
+                    .get_node_proportional_volume(idx, max_volume)?)
+            }
+            Self::VirtualStorageProportionalVolume(idx) => {
+                let max_volume = model.get_virtual_storage_node(idx)?.get_max_volume(model, state)?;
+                Ok(state
+                    .get_network_state()
+                    .get_virtual_storage_proportional_volume(idx, max_volume)?)
+            }
+            Self::AggregatedNodeProportionalVolume(idx) => {
+                let node = model.get_aggregated_storage_node(idx)?;
+                let volume: f64 = node
+                    .nodes
+                    .iter()
+                    .map(|idx| state.get_network_state().get_node_volume(idx))
+                    .sum::<Result<_, _>>()?;
+
+                let max_volume: f64 = node
+                    .nodes
+                    .iter()
+                    .map(|idx| model.get_node(idx)?.get_current_max_volume(model, state))
+                    .sum::<Result<_, _>>()?;
+                // TODO handle divide by zero
+                Ok(volume / max_volume)
+            }
+            Self::NodeInFlowDeficit(idx) => {
+                let node = model.get_node(idx)?;
+                let flow = state.get_network_state().get_node_in_flow(idx)?;
+                let max_flow = node.get_current_max_flow(model, state)?;
+                Ok(max_flow - flow)
+            }
+        }
+    }
+}

pywr-core/src/lib.rs

@@ -2,6 +2,7 @@
 
 extern crate core;
 
+use crate::derived_metric::DerivedMetricIndex;
 use crate::node::NodeIndex;
 use crate::parameters::{IndexParameterIndex, MultiValueParameterIndex, ParameterIndex};
 use crate::recorders::RecorderIndex;
@@ -11,6 +12,7 @@ use thiserror::Error;
 
 pub mod aggregated_node;
 mod aggregated_storage_node;
+pub mod derived_metric;
 pub mod edge;
 pub mod metric;
 pub mod model;
@@ -55,6 +57,10 @@ pub enum PywrError {
     RecorderIndexNotFound,
     #[error("recorder not found")]
     RecorderNotFound,
+    #[error("derived metric not found")]
+    DerivedMetricNotFound,
+    #[error("derived metric index {0} not found")]
+    DerivedMetricIndexNotFound(DerivedMetricIndex),
     #[error("node name `{0}` already exists")]
     NodeNameAlreadyExists(String),
     #[error("parameter name `{0}` already exists at index {1}")]
@@ -127,6 +133,8 @@ pub enum PywrError {
     ParameterVariableValuesIncorrectLength,
     #[error("missing solver features")]
     MissingSolverFeatures,
+    #[error("parameters do not provide an initial value")]
+    ParameterNoInitialValue,
 }
 
 // Python errors

pywr-core/src/metric.rs

@@ -1,54 +1,33 @@
 use crate::aggregated_node::AggregatedNodeIndex;
 use crate::aggregated_storage_node::AggregatedStorageNodeIndex;
+use crate::derived_metric::DerivedMetricIndex;
 use crate::edge::EdgeIndex;
 use crate::model::Model;
 use crate::node::NodeIndex;
-use crate::parameters::{MultiValueParameterIndex, ParameterIndex};
+use crate::parameters::{IndexParameterIndex, MultiValueParameterIndex, ParameterIndex};
 use crate::state::State;
 use crate::virtual_storage::VirtualStorageIndex;
 use crate::PywrError;
-
-#[derive(Clone, Debug, PartialEq)]
-pub struct VolumeBetweenControlCurves {
-    max_volume: Box<Metric>,
-    upper: Option<Box<Metric>>,
-    lower: Option<Box<Metric>>,
-}
-
-impl VolumeBetweenControlCurves {
-    pub fn new(max_volume: Metric, upper: Option<Metric>, lower: Option<Metric>) -> Self {
-        Self {
-            max_volume: Box::new(max_volume),
-            upper: upper.map(Box::new),
-            lower: lower.map(Box::new),
-        }
-    }
-}
-
 #[derive(Clone, Debug, PartialEq)]
 pub enum Metric {
     NodeInFlow(NodeIndex),
     NodeOutFlow(NodeIndex),
     NodeVolume(NodeIndex),
-    NodeInFlowDeficit(NodeIndex),
-    NodeProportionalVolume(NodeIndex),
     AggregatedNodeInFlow(AggregatedNodeIndex),
     AggregatedNodeOutFlow(AggregatedNodeIndex),
     AggregatedNodeVolume(AggregatedStorageNodeIndex),
-    AggregatedNodeProportionalVolume(AggregatedStorageNodeIndex),
     EdgeFlow(EdgeIndex),
     ParameterValue(ParameterIndex),
     MultiParameterValue((MultiValueParameterIndex, String)),
     VirtualStorageVolume(VirtualStorageIndex),
-    VirtualStorageProportionalVolume(VirtualStorageIndex),
-    VolumeBetweenControlCurves(VolumeBetweenControlCurves),
     MultiNodeInFlow {
         indices: Vec<NodeIndex>,
         name: String,
         sub_name: Option<String>,
     },
     // TODO implement other MultiNodeXXX variants
     Constant(f64),
+    DerivedMetric(DerivedMetricIndex),
 }
 
 impl Metric {
@@ -73,22 +52,12 @@ impl Metric {
                     .map(|idx| state.get_network_state().get_node_out_flow(idx))
                     .sum::<Result<_, _>>()
             }
-            Metric::NodeProportionalVolume(idx) => {
-                let max_volume = model.get_node(idx)?.get_current_max_volume(model, state)?;
-                Ok(state
-                    .get_network_state()
-                    .get_node_proportional_volume(idx, max_volume)?)
-            }
+
             Metric::EdgeFlow(idx) => Ok(state.get_network_state().get_edge_flow(idx)?),
             Metric::ParameterValue(idx) => Ok(state.get_parameter_value(*idx)?),
             Metric::MultiParameterValue((idx, key)) => Ok(state.get_multi_parameter_value(*idx, key)?),
             Metric::VirtualStorageVolume(idx) => Ok(state.get_network_state().get_virtual_storage_volume(idx)?),
-            Metric::VirtualStorageProportionalVolume(idx) => {
-                let max_volume = model.get_virtual_storage_node(idx)?.get_max_volume(model, state)?;
-                Ok(state
-                    .get_network_state()
-                    .get_virtual_storage_proportional_volume(idx, max_volume)?)
-            }
+            Metric::DerivedMetric(idx) => state.get_derived_metric_value(*idx),
             Metric::Constant(v) => Ok(*v),
             Metric::AggregatedNodeVolume(idx) => {
                 let node = model.get_aggregated_storage_node(idx)?;
@@ -97,49 +66,29 @@ impl Metric {
                     .map(|idx| state.get_network_state().get_node_volume(idx))
                     .sum::<Result<_, _>>()
             }
-            Metric::AggregatedNodeProportionalVolume(idx) => {
-                let node = model.get_aggregated_storage_node(idx)?;
-                let volume: f64 = node
-                    .nodes
-                    .iter()
-                    .map(|idx| state.get_network_state().get_node_volume(idx))
-                    .sum::<Result<_, _>>()?;
 
-                let max_volume: f64 = node
-                    .nodes
-                    .iter()
-                    .map(|idx| model.get_node(idx)?.get_current_max_volume(model, state))
-                    .sum::<Result<_, _>>()?;
-                // TODO handle divide by zero
-                Ok(volume / max_volume)
-            }
             Metric::MultiNodeInFlow { indices, .. } => {
                 let flow = indices
                     .iter()
                     .map(|idx| state.get_network_state().get_node_in_flow(idx))
                     .sum::<Result<_, _>>()?;
                 Ok(flow)
             }
-            Metric::NodeInFlowDeficit(idx) => {
-                let node = model.get_node(idx)?;
-                let flow = state.get_network_state().get_node_in_flow(idx)?;
-                let max_flow = node.get_current_max_flow(model, state)?;
-                Ok(max_flow - flow)
-            }
-            Metric::VolumeBetweenControlCurves(vol) => {
-                let max_volume = vol.max_volume.get_value(model, state)?;
-                let lower = vol
-                    .lower
-                    .as_ref()
-                    .map_or(Ok(0.0), |metric| metric.get_value(model, state))?;
-                let upper = vol
-                    .upper
-                    .as_ref()
-                    .map_or(Ok(1.0), |metric| metric.get_value(model, state))?;
+        }
+    }
+}
 
-                // TODO handle invalid bounds
-                Ok(max_volume * (upper - lower))
-            }
+#[derive(Clone, Debug, PartialEq)]
+pub enum IndexMetric {
+    IndexParameterValue(IndexParameterIndex),
+    Constant(usize),
+}
+
+impl IndexMetric {
+    pub fn get_value(&self, model: &Model, state: &State) -> Result<usize, PywrError> {
+        match self {
+            Self::IndexParameterValue(idx) => state.get_parameter_index(*idx),
+            Self::Constant(i) => Ok(*i),
         }
     }
 }