feat: Utilise constant parameters in aggregated node LP creation.

Where possible utilise constant parameters to try to get aggregated
node factor pairs during LP creation. This makes it possible to
support constant factors in LP solvers without an API for updating
the constraint matrix. It should also be generally more efficient
by only updating non-constant factors.

pywr-core/src/aggregated_node.rs

@@ -1,7 +1,7 @@
 use crate::metric::MetricF64;
 use crate::network::Network;
 use crate::node::{Constraint, FlowConstraints, NodeMeta};
-use crate::state::State;
+use crate::state::{ConstParameterValues, State};
 use crate::{NodeIndex, PywrError};
 use std::ops::{Deref, DerefMut};
 
@@ -64,6 +64,16 @@ pub enum Factors {
     Ratio(Vec<MetricF64>),
 }
 
+impl Factors {
+    /// Returns true if all factors are constant
+    pub fn is_constant(&self) -> bool {
+        match self {
+            Factors::Proportion(factors) => factors.iter().all(|f| f.is_constant()),
+            Factors::Ratio(factors) => factors.iter().all(|f| f.is_constant()),
+        }
+    }
+}
+
 #[derive(Debug, PartialEq)]
 pub struct AggregatedNode {
     meta: NodeMeta<AggregatedNodeIndex>,
@@ -102,6 +112,41 @@ impl NodeFactorPair {
     }
 }
 
+/// A constant node factor. If the factor is non-constant, the factor value here is `None`.
+#[derive(Debug, PartialEq, Copy, Clone)]
+pub struct NodeConstFactor {
+    pub index: NodeIndex,
+    pub factor: Option<f64>,
+}
+
+impl NodeConstFactor {
+    fn new(node: NodeIndex, factor: Option<f64>) -> Self {
+        Self { index: node, factor }
+    }
+}
+
+/// A pair of nodes and their factors
+#[derive(Debug, PartialEq, Copy, Clone)]
+pub struct NodeConstFactorPair {
+    pub node0: NodeConstFactor,
+    pub node1: NodeConstFactor,
+}
+
+impl NodeConstFactorPair {
+    fn new(node0: NodeConstFactor, node1: NodeConstFactor) -> Self {
+        Self { node0, node1 }
+    }
+
+    /// Return the ratio of the two factors (node0 / node1). If either factor is `None`,
+    /// the ratio is also `None`.
+    pub fn ratio(&self) -> Option<f64> {
+        match (self.node0.factor, self.node1.factor) {
+            (Some(f0), Some(f1)) => Some(f0 / f1),
+            _ => None,
+        }
+    }
+}
+
 impl AggregatedNode {
     pub fn new(
         index: &AggregatedNodeIndex,
@@ -140,6 +185,15 @@ impl AggregatedNode {
         self.nodes.to_vec()
     }
 
+    /// Does the aggregated node have factors?
+    pub fn has_factors(&self) -> bool {
+        self.factors.is_some()
+    }
+
+    /// Does the aggregated node have constant factors?
+    pub fn has_const_factors(&self) -> bool {
+        self.factors.as_ref().map(|f| f.is_constant()).unwrap_or(false)
+    }
     pub fn set_factors(&mut self, factors: Option<Factors>) {
         self.factors = factors;
     }
@@ -148,6 +202,7 @@ impl AggregatedNode {
         self.factors.as_ref()
     }
 
+    /// Return normalised factor pairs
     pub fn get_factor_node_pairs(&self) -> Option<Vec<(NodeIndex, NodeIndex)>> {
         if self.factors.is_some() {
             let n0 = self.nodes[0];
@@ -158,6 +213,21 @@ impl AggregatedNode {
         }
     }
 
+    /// Return constant normalised factor pairs
+    pub fn get_const_norm_factor_pairs(&self, values: &ConstParameterValues) -> Option<Vec<NodeConstFactorPair>> {
+        if let Some(factors) = &self.factors {
+            let pairs = match factors {
+                Factors::Proportion(prop_factors) => {
+                    get_const_norm_proportional_factor_pairs(prop_factors, &self.nodes, values)
+                }
+                Factors::Ratio(ratio_factors) => get_const_norm_ratio_factor_pairs(ratio_factors, &self.nodes, values),
+            };
+            Some(pairs)
+        } else {
+            None
+        }
+    }
+
     /// Return normalised factor pairs
     ///
     pub fn get_norm_factor_pairs(&self, model: &Network, state: &State) -> Option<Vec<NodeFactorPair>> {
@@ -225,7 +295,11 @@ impl AggregatedNode {
     }
 }
 
-/// Proportional factors
+/// Calculate factor pairs for proportional factors.
+///
+/// There should be one less factor than node indices. The factors correspond to each of the node
+/// indices after the first. Factor pairs relating the first index to each of the other indices are
+/// calculated. This requires the sum of the factors to be greater than 0.0 and less than 1.0.
 fn get_norm_proportional_factor_pairs(
     factors: &[MetricF64],
     nodes: &[NodeIndex],
@@ -263,7 +337,71 @@ fn get_norm_proportional_factor_pairs(
         .collect::<Vec<_>>()
 }
 
-/// Ratio factors
+/// Calculate constant factor pairs for proportional factors.
+///
+/// There should be one less factor than node indices. The factors correspond to each of the node
+/// indices after the first. Factor pairs relating the first index to each of the other indices are
+/// calculated. This requires the sum of the factors to be greater than 0.0 and less than 1.0. If
+/// any of the factors are not constant, the factor pairs will contain `None` values.
+fn get_const_norm_proportional_factor_pairs(
+    factors: &[MetricF64],
+    nodes: &[NodeIndex],
+    values: &ConstParameterValues,
+) -> Vec<NodeConstFactorPair> {
+    if factors.len() != nodes.len() - 1 {
+        panic!("Found {} proportional factors and {} nodes in aggregated node. The number of proportional factors should equal one less than the number of nodes.", factors.len(), nodes.len());
+    }
+
+    // First get the current factor values
+    let values: Vec<Option<f64>> = factors
+        .iter()
+        .map(|f| f.try_get_constant_value(values))
+        .collect::<Result<Vec<_>, PywrError>>()
+        .expect("Failed to get current factor values.");
+
+    let n0 = nodes[0];
+
+    // To calculate the factors we require that every factor is available.
+    if values.iter().any(|v| v.is_none()) {
+        // At least one factor is not available; therefore we can not calculate "f0"
+        nodes
+            .iter()
+            .skip(1)
+            .zip(values)
+            .map(move |(&n1, f1)| {
+                NodeConstFactorPair::new(NodeConstFactor::new(n0, None), NodeConstFactor::new(n1, f1))
+            })
+            .collect::<Vec<_>>()
+    } else {
+        // All factors are available; therefore we can calculate "f0"
+        // TODO do we need to assert that each individual factor is positive?
+        let total: f64 = values
+            .iter()
+            .map(|v| v.expect("Factor is `None`; this should be impossible."))
+            .sum();
+        if total < 0.0 {
+            panic!("Proportional factors are too small or negative.");
+        }
+        if total >= 1.0 {
+            panic!("Proportional factors are too large.")
+        }
+
+        let f0 = Some(1.0 - total);
+
+        nodes
+            .iter()
+            .skip(1)
+            .zip(values)
+            .map(move |(&n1, f1)| NodeConstFactorPair::new(NodeConstFactor::new(n0, f0), NodeConstFactor::new(n1, f1)))
+            .collect::<Vec<_>>()
+    }
+}
+
+/// Calculate factor pairs for ratio factors.
+///
+/// The number of node indices and factors should be equal. The factors correspond to each of the
+/// node indices. Factor pairs relating the first index to each of the other indices are calculated.
+/// This requires that the factors are all non-zero.
 fn get_norm_ratio_factor_pairs(
     factors: &[MetricF64],
     nodes: &[NodeIndex],
@@ -290,12 +428,45 @@ fn get_norm_ratio_factor_pairs(
         .collect::<Vec<_>>()
 }
 
+/// Constant ratio factors using constant values if they are available. If they are not available,
+/// the factors are `None`.
+fn get_const_norm_ratio_factor_pairs(
+    factors: &[MetricF64],
+    nodes: &[NodeIndex],
+    values: &ConstParameterValues,
+) -> Vec<NodeConstFactorPair> {
+    if factors.len() != nodes.len() {
+        panic!("Found {} ratio factors and {} nodes in aggregated node. The number of ratio factors should equal the number of nodes.", factors.len(), nodes.len());
+    }
+
+    let n0 = nodes[0];
+    // Try to convert the factor into a constant
+
+    let f0 = factors[0]
+        .try_get_constant_value(values)
+        .unwrap_or_else(|e| panic!("Failed to get constant value for factor: {}", e));
+
+    nodes
+        .iter()
+        .zip(factors)
+        .skip(1)
+        .map(move |(&n1, f1)| {
+            let v1 = f1
+                .try_get_constant_value(values)
+                .unwrap_or_else(|e| panic!("Failed to get constant value for factor: {}", e));
+
+            NodeConstFactorPair::new(NodeConstFactor::new(n0, f0), NodeConstFactor::new(n1, v1))
+        })
+        .collect::<Vec<_>>()
+}
+
 #[cfg(test)]
 mod tests {
     use crate::aggregated_node::Factors;
     use crate::metric::MetricF64;
     use crate::models::Model;
     use crate::network::Network;
+    use crate::parameters::MonthlyProfileParameter;
     use crate::recorders::AssertionRecorder;
     use crate::test_utils::{default_time_domain, run_all_solvers};
     use ndarray::Array2;
@@ -344,6 +515,56 @@ mod tests {
 
         let model = Model::new(default_time_domain().into(), network);
 
-        run_all_solvers(&model, &["cbc", "highs"]);
+        run_all_solvers(&model, &[]);
+    }
+
+    /// Test the factors forcing a simple ratio of flow that varies over time
+    ///
+    /// The model has a single input that diverges to two links and respective output nodes.
+    #[test]
+    fn test_simple_factor_profile() {
+        let mut network = Network::default();
+
+        let input_node = network.add_input_node("input", None).unwrap();
+        let link_node0 = network.add_link_node("link", Some("0")).unwrap();
+        let output_node0 = network.add_output_node("output", Some("0")).unwrap();
+
+        network.connect_nodes(input_node, link_node0).unwrap();
+        network.connect_nodes(link_node0, output_node0).unwrap();
+
+        let link_node1 = network.add_link_node("link", Some("1")).unwrap();
+        let output_node1 = network.add_output_node("output", Some("1")).unwrap();
+
+        network.connect_nodes(input_node, link_node1).unwrap();
+        network.connect_nodes(link_node1, output_node1).unwrap();
+
+        let factor_profile = MonthlyProfileParameter::new("factor-profile", [2.0; 12], None);
+        let factor_profile_idx = network.add_simple_parameter(Box::new(factor_profile)).unwrap();
+
+        let factors = Some(Factors::Ratio(vec![factor_profile_idx.into(), 1.0.into()]));
+
+        let _agg_node = network.add_aggregated_node("agg-node", None, &[link_node0, link_node1], factors);
+
+        // Setup a demand on output-0
+        let output_node = network.get_mut_node_by_name("output", Some("0")).unwrap();
+        output_node.set_max_flow_constraint(Some(100.0.into())).unwrap();
+
+        output_node.set_cost(Some((-10.0).into()));
+
+        // Set-up assertion for "input" node
+        let idx = network.get_node_by_name("link", Some("0")).unwrap().index();
+        let expected = Array2::from_elem((366, 10), 100.0);
+        let recorder = AssertionRecorder::new("link-0-flow", MetricF64::NodeOutFlow(idx), expected, None, None);
+        network.add_recorder(Box::new(recorder)).unwrap();
+
+        // Set-up assertion for "input" node
+        let idx = network.get_node_by_name("link", Some("1")).unwrap().index();
+        let expected = Array2::from_elem((366, 10), 50.0);
+        let recorder = AssertionRecorder::new("link-0-flow", MetricF64::NodeOutFlow(idx), expected, None, None);
+        network.add_recorder(Box::new(recorder)).unwrap();
+
+        let model = Model::new(default_time_domain().into(), network);
+
+        run_all_solvers(&model, &["cbc"]);
     }
 }

pywr-core/src/metric.rs

@@ -45,6 +45,22 @@ impl SimpleMetricF64 {
             SimpleMetricF64::Constant(m) => m.get_value(&values.get_constant_values()),
         }
     }
+
+    /// Try to get the constant value of the metric, if it is a constant value.
+    pub fn try_get_constant_value(&self, values: &ConstParameterValues) -> Result<Option<f64>, PywrError> {
+        match self {
+            SimpleMetricF64::Constant(c) => c.get_value(values).map(Some),
+            _ => Ok(None),
+        }
+    }
+
+    /// Returns true if the metric is a constant value.
+    pub fn is_constant(&self) -> bool {
+        match self {
+            SimpleMetricF64::Constant(_) => true,
+            _ => false,
+        }
+    }
 }
 
 #[derive(Clone, Debug, PartialEq)]
@@ -124,6 +140,21 @@ impl MetricF64 {
             MetricF64::Simple(s) => s.get_value(&state.get_simple_parameter_values()),
         }
     }
+
+    /// Try to get the constant value of the metric, if it is a constant value.
+    pub fn try_get_constant_value(&self, values: &ConstParameterValues) -> Result<Option<f64>, PywrError> {
+        match self {
+            MetricF64::Simple(s) => s.try_get_constant_value(values),
+            _ => Ok(None),
+        }
+    }
+
+    pub fn is_constant(&self) -> bool {
+        match self {
+            MetricF64::Simple(s) => s.is_constant(),
+            _ => false,
+        }
+    }
 }
 
 impl TryFrom<MetricF64> for SimpleMetricF64 {

pywr-core/src/models/multi.rs

@@ -159,7 +159,7 @@ impl MultiNetworkModel {
                 .network
                 .setup_network(timesteps, scenario_indices, entry.parameters.len())?;
             let recorder_state = entry.network.setup_recorders(&self.domain)?;
-            let solver = entry.network.setup_solver::<S>(scenario_indices, settings)?;
+            let solver = entry.network.setup_solver::<S>(scenario_indices, &state, settings)?;
 
             states.push(state);
             recorder_states.push(recorder_state);

pywr-core/src/models/simple.rs

@@ -78,7 +78,7 @@ impl Model {
 
         let state = self.network.setup_network(timesteps, scenario_indices, 0)?;
         let recorder_state = self.network.setup_recorders(&self.domain)?;
-        let solvers = self.network.setup_solver::<S>(scenario_indices, settings)?;
+        let solvers = self.network.setup_solver::<S>(scenario_indices, &state, settings)?;
 
         Ok(ModelState {
             current_time_step_idx: 0,