problemreductions/
parameters.rs

1//! Symbolic parameter transformations carried by reduction rules.
2
3use crate::expr::{AlgebraicAnalysis, Expr, ExprNode, ExprNodeId, Symbol};
4use crate::types::ProblemParameters;
5use num_bigint::{BigInt, BigUint, Sign};
6use num_rational::BigRational;
7use num_traits::{One, Signed, Zero};
8use std::collections::{BTreeMap, HashMap, HashSet};
9use std::sync::Arc;
10
11/// What a reduction promises about one parameter formula.
12#[derive(Clone, Copy, Debug, PartialEq, Eq, serde::Serialize, serde::Deserialize)]
13#[serde(rename_all = "snake_case")]
14pub enum ParameterRelation {
15    Exact,
16    UpperBound,
17}
18
19/// Symbolic predictions with independent accuracy and availability per field.
20#[derive(Clone, Debug)]
21pub struct ParameterTransform {
22    edge: Box<str>,
23    fields: Vec<ParameterField>,
24    unavailable: BTreeMap<Box<str>, ParameterTransformError>,
25}
26
27#[derive(Clone, Debug)]
28struct ParameterField {
29    name: Box<str>,
30    expression: Expr,
31    relation: ParameterRelation,
32    plan: Plan,
33}
34
35#[derive(Clone, Debug)]
36struct Plan(Arc<PlanNode>);
37
38#[derive(Debug)]
39enum PlanNode {
40    Const(BigRational),
41    Var(Symbol),
42    Add(Box<[Plan]>),
43    Mul(Box<[Plan]>),
44    Pow(Plan, BigInt),
45}
46
47impl Plan {
48    fn identity(&self) -> usize {
49        Arc::as_ptr(&self.0) as usize
50    }
51}
52
53impl ParameterTransform {
54    pub fn new<I, N>(
55        edge: impl Into<Box<str>>,
56        relation: ParameterRelation,
57        fields: I,
58    ) -> Result<Self, ParameterTransformError>
59    where
60        I: IntoIterator<Item = (N, Expr)>,
61        N: Into<Box<str>>,
62    {
63        Self::from_fields(
64            edge,
65            fields
66                .into_iter()
67                .map(|(name, expression)| (name, relation, expression)),
68        )
69    }
70
71    /// Construct predictions whose relations may differ by target field.
72    pub fn from_fields<I, N>(
73        edge: impl Into<Box<str>>,
74        fields: I,
75    ) -> Result<Self, ParameterTransformError>
76    where
77        I: IntoIterator<Item = (N, ParameterRelation, Expr)>,
78        N: Into<Box<str>>,
79    {
80        let edge = edge.into();
81        let mut names = HashSet::new();
82        let mut raw_fields = Vec::new();
83        for (name, relation, expression) in fields {
84            let name = name.into();
85            if let Err(error) = Symbol::new(name.clone()) {
86                return Err(ParameterTransformError::InvalidTargetField {
87                    edge,
88                    field: name,
89                    reason: error.to_string().into(),
90                });
91            }
92            if !names.insert(name.clone()) {
93                return Err(ParameterTransformError::DuplicateTargetField { edge, field: name });
94            }
95            raw_fields.push((name, relation, expression));
96        }
97
98        let expressions = raw_fields
99            .iter()
100            .map(|(_, _, expression)| expression)
101            .collect::<Vec<_>>();
102        let analysis = AlgebraicAnalysis::new(&expressions);
103        let mut plans = HashMap::new();
104        let fields = raw_fields
105            .into_iter()
106            .map(|(name, relation, expression)| {
107                let plan = compile(&expression, &analysis, &mut plans).map_err(|failure| {
108                    validation_error(edge.clone(), name.clone(), expression.to_string(), failure)
109                })?;
110                Ok(ParameterField {
111                    name,
112                    expression,
113                    relation,
114                    plan,
115                })
116            })
117            .collect::<Result<Vec<_>, _>>()?;
118
119        Ok(Self {
120            edge,
121            fields,
122            unavailable: BTreeMap::new(),
123        })
124    }
125
126    pub fn edge(&self) -> &str {
127        &self.edge
128    }
129
130    pub fn relation(&self, target_field: &str) -> Option<ParameterRelation> {
131        self.fields
132            .iter()
133            .find(|field| field.name.as_ref() == target_field)
134            .map(|field| field.relation)
135    }
136
137    /// Why this target field could not be predicted, including upstream causes.
138    pub fn unavailable(&self, target_field: &str) -> Option<&ParameterTransformError> {
139        self.unavailable.get(target_field)
140    }
141
142    pub(crate) fn declare_unavailable(&mut self, field: &str, reason: &str) {
143        self.unavailable.insert(
144            field.into(),
145            ParameterTransformError::Unavailable {
146                edge: self.edge.clone(),
147                field: field.into(),
148                reason: reason.into(),
149            },
150        );
151    }
152
153    pub fn expressions(&self) -> impl Iterator<Item = (&str, &Expr)> {
154        self.fields
155            .iter()
156            .map(|field| (field.name.as_ref(), &field.expression))
157    }
158
159    pub fn get(&self, target_field: &str) -> Option<&Expr> {
160        self.fields
161            .iter()
162            .find(|field| field.name.as_ref() == target_field)
163            .map(|field| &field.expression)
164    }
165
166    pub fn evaluate(
167        &self,
168        input: &ProblemParameters,
169    ) -> Result<ProblemParameters, ParameterTransformError> {
170        let mut memo = HashMap::new();
171        let mut output = Vec::with_capacity(self.fields.len());
172        for field in &self.fields {
173            let value = evaluate_plan(&field.plan, input, &mut memo).map_err(|failure| {
174                evaluation_error(self.edge.clone(), field.name.clone(), failure)
175            })?;
176            if value.is_negative() {
177                return Err(ParameterTransformError::NegativeResult {
178                    edge: self.edge.clone(),
179                    field: field.name.clone(),
180                    value,
181                });
182            }
183            let value = if field.relation == ParameterRelation::Exact {
184                if !value.is_integer() {
185                    return Err(ParameterTransformError::NonIntegralResult {
186                        edge: self.edge.clone(),
187                        field: field.name.clone(),
188                        value: value.to_string().into(),
189                    });
190                }
191                value.to_integer().magnitude().clone()
192            } else {
193                ceil_nonnegative(&value)
194            };
195            let value =
196                u64::try_from(&value).map_err(|_| ParameterTransformError::OutputOutOfRange {
197                    field: field.name.clone(),
198                    value: value.clone(),
199                })?;
200            output.push((field.name.to_string(), value));
201        }
202        Ok(ProblemParameters::from_owned(output))
203    }
204
205    pub fn compose(
206        &self,
207        next: &ParameterTransform,
208        edge: impl Into<Box<str>>,
209    ) -> Result<ParameterTransform, ParameterTransformError> {
210        let edge = edge.into();
211        let replacements: HashMap<&str, &Expr> = self.expressions().collect();
212        let mut fields = Vec::new();
213        let mut unavailable = next.unavailable.clone();
214        for field in &next.fields {
215            let result = (|| {
216                let mut bounded_input = false;
217                for input in field.expression.variables() {
218                    if let Some(cause) = self.unavailable(input) {
219                        return Err(ParameterTransformError::UnavailableInput {
220                            edge: next.edge.clone(),
221                            field: field.name.clone(),
222                            input_field: input.into(),
223                            cause: Box::new(cause.clone()),
224                        });
225                    }
226                    bounded_input |= self.relation(input) == Some(ParameterRelation::UpperBound);
227                }
228                let expression = if bounded_input {
229                    positive_polynomial_hull(&field.expression).ok_or_else(|| {
230                        ParameterTransformError::CannotPropagateUpperBound {
231                            edge: next.edge.clone(),
232                            field: field.name.clone(),
233                            expression: field.expression.to_string().into(),
234                        }
235                    })?
236                } else {
237                    field.expression.clone()
238                };
239                let expression =
240                    expression
241                        .substitute_complete(&replacements)
242                        .map_err(|error| ParameterTransformError::MissingCompositionInput {
243                            edge: next.edge.clone(),
244                            field: field.name.clone(),
245                            input_fields: error.missing_variables().map(Box::<str>::from).collect(),
246                        })?;
247                let relation = if bounded_input {
248                    ParameterRelation::UpperBound
249                } else {
250                    field.relation
251                };
252                Ok((field.name.clone(), relation, expression))
253            })();
254            match result {
255                Ok(prediction) => fields.push(prediction),
256                Err(error) => {
257                    unavailable.insert(field.name.clone(), error);
258                }
259            }
260        }
261        let mut composed = Self::from_fields(edge, fields)?;
262        composed.unavailable = unavailable;
263        Ok(composed)
264    }
265}
266
267type Monomial = BTreeMap<Symbol, BigUint>;
268type Polynomial = BTreeMap<Monomial, BigRational>;
269
270fn positive_polynomial_hull(expression: &Expr) -> Option<Expr> {
271    let polynomial = polynomial(expression)?;
272    let terms = polynomial
273        .into_iter()
274        .filter(|(_, coefficient)| coefficient.is_positive())
275        .map(|(monomial, coefficient)| {
276            monomial
277                .into_iter()
278                .fold(Expr::constant(coefficient), |term, (variable, exponent)| {
279                    term * Expr::pow(
280                        Expr::variable(variable.as_str()),
281                        Expr::integer(BigInt::from(exponent)),
282                    )
283                })
284        });
285    Some(terms.fold(Expr::integer(0), |sum, term| sum + term))
286}
287
288fn polynomial(expression: &Expr) -> Option<Polynomial> {
289    match expression.node() {
290        ExprNode::Const(value) => Some(BTreeMap::from([(BTreeMap::new(), value.clone())])),
291        ExprNode::Var(variable) => Some(BTreeMap::from([(
292            BTreeMap::from([(variable.clone(), BigUint::one())]),
293            BigRational::one(),
294        )])),
295        ExprNode::Add(values) => values.iter().try_fold(BTreeMap::new(), |sum, value| {
296            Some(add_polynomials(sum, polynomial(value)?))
297        }),
298        ExprNode::Mul(values) => values.iter().try_fold(
299            BTreeMap::from([(BTreeMap::new(), BigRational::one())]),
300            |product, value| Some(multiply_polynomials(product, polynomial(value)?)),
301        ),
302        ExprNode::Pow(base, exponent) => {
303            let ExprNode::Const(exponent) = exponent.node() else {
304                return None;
305            };
306            if !exponent.is_integer() {
307                return None;
308            }
309            if exponent.is_negative() {
310                let ExprNode::Const(base) = base.node() else {
311                    return None;
312                };
313                if base.is_zero() {
314                    return None;
315                }
316                return Some(BTreeMap::from([(
317                    BTreeMap::new(),
318                    pow_rational(base.clone(), &exponent.to_integer()),
319                )]));
320            }
321            let mut exponent = exponent.to_integer().magnitude().clone();
322            let mut base = polynomial(base)?;
323            let mut result = BTreeMap::from([(BTreeMap::new(), BigRational::one())]);
324            while !exponent.is_zero() {
325                if exponent.bit(0) {
326                    result = multiply_polynomials(result, base.clone());
327                }
328                exponent >>= 1usize;
329                if !exponent.is_zero() {
330                    base = multiply_polynomials(base.clone(), base);
331                }
332            }
333            Some(result)
334        }
335        ExprNode::Exp(_) | ExprNode::Log(_) | ExprNode::Factorial(_) => None,
336    }
337}
338
339fn add_polynomials(mut left: Polynomial, right: Polynomial) -> Polynomial {
340    for (monomial, right_coefficient) in right {
341        *left.entry(monomial).or_insert_with(BigRational::zero) += right_coefficient;
342    }
343    left.retain(|_, coefficient| !coefficient.is_zero());
344    left
345}
346
347fn multiply_polynomials(left: Polynomial, right: Polynomial) -> Polynomial {
348    let mut product = Polynomial::new();
349    for (left_monomial, left_coefficient) in left {
350        for (right_monomial, right_coefficient) in &right {
351            let mut monomial = left_monomial.clone();
352            for (variable, exponent) in right_monomial {
353                *monomial.entry(variable.clone()).or_default() += exponent;
354            }
355            *product.entry(monomial).or_insert_with(BigRational::zero) +=
356                &left_coefficient * right_coefficient;
357        }
358    }
359    product.retain(|_, coefficient| !coefficient.is_zero());
360    product
361}
362
363fn compile(
364    expression: &Expr,
365    analysis: &AlgebraicAnalysis,
366    memo: &mut HashMap<ExprNodeId, Plan>,
367) -> Result<Plan, ValidationFailure> {
368    if let Some(plan) = memo.get(&expression.node_identity()) {
369        return Ok(plan.clone());
370    }
371    let node = match expression.node() {
372        ExprNode::Const(value) => PlanNode::Const(value.clone()),
373        ExprNode::Var(symbol) => PlanNode::Var(symbol.clone()),
374        ExprNode::Add(values) => PlanNode::Add(
375            values
376                .iter()
377                .map(|value| compile(value, analysis, memo))
378                .collect::<Result<Vec<_>, _>>()?
379                .into_boxed_slice(),
380        ),
381        ExprNode::Mul(values) => PlanNode::Mul(
382            values
383                .iter()
384                .map(|value| compile(value, analysis, memo))
385                .collect::<Result<Vec<_>, _>>()?
386                .into_boxed_slice(),
387        ),
388        ExprNode::Pow(base, exponent) => {
389            let Some(exponent) = analysis.facts(exponent).exact_rational.as_ref() else {
390                return Err(ValidationFailure::NonIntegralConstantExponent(
391                    exponent.to_string().into(),
392                ));
393            };
394            if !exponent.is_integer() {
395                return Err(ValidationFailure::NonIntegralConstantExponent(
396                    exponent.to_string().into(),
397                ));
398            }
399            PlanNode::Pow(compile(base, analysis, memo)?, exponent.to_integer())
400        }
401        ExprNode::Exp(_) => return Err(ValidationFailure::UnsupportedOperator("exp")),
402        ExprNode::Log(_) => return Err(ValidationFailure::UnsupportedOperator("log")),
403        ExprNode::Factorial(_) => {
404            return Err(ValidationFailure::UnsupportedOperator("factorial"));
405        }
406    };
407    let plan = Plan(Arc::new(node));
408    memo.insert(expression.node_identity(), plan.clone());
409    Ok(plan)
410}
411
412fn evaluate_plan(
413    plan: &Plan,
414    input: &ProblemParameters,
415    memo: &mut HashMap<usize, BigRational>,
416) -> Result<BigRational, EvaluationFailure> {
417    if let Some(value) = memo.get(&plan.identity()) {
418        return Ok(value.clone());
419    }
420    let value = match plan.0.as_ref() {
421        PlanNode::Const(value) => value.clone(),
422        PlanNode::Var(symbol) => BigRational::from_integer(BigInt::from(
423            input
424                .get(symbol.as_str())
425                .ok_or_else(|| EvaluationFailure::MissingInputField(symbol.to_string().into()))?,
426        )),
427        PlanNode::Add(values) => values.iter().try_fold(BigRational::zero(), |sum, value| {
428            Ok(sum + evaluate_plan(value, input, memo)?)
429        })?,
430        PlanNode::Mul(values) => values
431            .iter()
432            .try_fold(BigRational::one(), |product, value| {
433                Ok(product * evaluate_plan(value, input, memo)?)
434            })?,
435        PlanNode::Pow(base, exponent) => {
436            let base = evaluate_plan(base, input, memo)?;
437            if exponent.sign() == Sign::Minus && base.is_zero() {
438                return Err(EvaluationFailure::DivisionByZero);
439            }
440            pow_rational(base, exponent)
441        }
442    };
443    memo.insert(plan.identity(), value.clone());
444    Ok(value)
445}
446
447fn pow_rational(mut base: BigRational, exponent: &BigInt) -> BigRational {
448    let negative = exponent.sign() == Sign::Minus;
449    let mut exponent = exponent.magnitude().clone();
450    let mut result = BigRational::one();
451    while !exponent.is_zero() {
452        if exponent.bit(0) {
453            result *= &base;
454        }
455        exponent >>= 1usize;
456        if !exponent.is_zero() {
457            base = &base * &base;
458        }
459    }
460    if negative {
461        result.recip()
462    } else {
463        result
464    }
465}
466
467fn ceil_nonnegative(value: &BigRational) -> BigUint {
468    ((value.numer() + value.denom() - BigInt::one()) / value.denom())
469        .magnitude()
470        .clone()
471}
472
473#[derive(Debug)]
474enum ValidationFailure {
475    NonIntegralConstantExponent(Box<str>),
476    UnsupportedOperator(&'static str),
477}
478
479#[derive(Debug)]
480enum EvaluationFailure {
481    MissingInputField(Box<str>),
482    DivisionByZero,
483}
484
485fn validation_error(
486    edge: Box<str>,
487    field: Box<str>,
488    expression: String,
489    failure: ValidationFailure,
490) -> ParameterTransformError {
491    match failure {
492        ValidationFailure::NonIntegralConstantExponent(exponent) => {
493            ParameterTransformError::NonIntegralConstantExponent {
494                edge,
495                field,
496                expression: expression.into(),
497                exponent,
498            }
499        }
500        ValidationFailure::UnsupportedOperator(operator) => {
501            ParameterTransformError::UnsupportedOperator {
502                edge,
503                field,
504                expression: expression.into(),
505                operator,
506            }
507        }
508    }
509}
510
511fn evaluation_error(
512    edge: Box<str>,
513    field: Box<str>,
514    failure: EvaluationFailure,
515) -> ParameterTransformError {
516    match failure {
517        EvaluationFailure::MissingInputField(input_field) => {
518            ParameterTransformError::MissingInputField {
519                edge,
520                field,
521                input_field,
522            }
523        }
524        EvaluationFailure::DivisionByZero => {
525            ParameterTransformError::DivisionByZero { edge, field }
526        }
527    }
528}
529
530/// Validation, composition, or evaluation failure for a [`ParameterTransform`].
531#[derive(Clone, Debug, PartialEq, Eq, thiserror::Error)]
532pub enum ParameterTransformError {
533    #[error("reduction `{edge}` target field `{field}` is unavailable: {reason}")]
534    Unavailable {
535        edge: Box<str>,
536        field: Box<str>,
537        reason: Box<str>,
538    },
539    #[error(
540        "reduction `{edge}` target field `{field}` depends on unavailable `{input_field}`: {cause}"
541    )]
542    UnavailableInput {
543        edge: Box<str>,
544        field: Box<str>,
545        input_field: Box<str>,
546        #[source]
547        cause: Box<ParameterTransformError>,
548    },
549    #[error("reduction `{edge}` has invalid target parameter field `{field}`: {reason}")]
550    InvalidTargetField {
551        edge: Box<str>,
552        field: Box<str>,
553        reason: Box<str>,
554    },
555    #[error("reduction `{edge}` declares target parameter field `{field}` more than once")]
556    DuplicateTargetField { edge: Box<str>, field: Box<str> },
557    #[error("reduction `{edge}` target field `{field}` has non-integral constant exponent `{exponent}` in `{expression}`")]
558    NonIntegralConstantExponent {
559        edge: Box<str>,
560        field: Box<str>,
561        expression: Box<str>,
562        exponent: Box<str>,
563    },
564    #[error("reduction `{edge}` target field `{field}` uses unsupported operator `{operator}` in `{expression}`")]
565    UnsupportedOperator {
566        edge: Box<str>,
567        field: Box<str>,
568        expression: Box<str>,
569        operator: &'static str,
570    },
571    #[error("reduction `{edge}` target field `{field}` cannot propagate an upper bound through `{expression}`")]
572    CannotPropagateUpperBound {
573        edge: Box<str>,
574        field: Box<str>,
575        expression: Box<str>,
576    },
577    #[error(
578        "reduction `{edge}` target field `{field}` is missing input parameter field `{input_field}`"
579    )]
580    MissingInputField {
581        edge: Box<str>,
582        field: Box<str>,
583        input_field: Box<str>,
584    },
585    #[error(
586        "reduction `{edge}` target field `{field}` is missing composition inputs {input_fields:?}"
587    )]
588    MissingCompositionInput {
589        edge: Box<str>,
590        field: Box<str>,
591        input_fields: Vec<Box<str>>,
592    },
593    #[error("reduction `{edge}` target field `{field}` divides by zero")]
594    DivisionByZero { edge: Box<str>, field: Box<str> },
595    #[error("reduction `{edge}` target field `{field}` evaluates to non-integral parameter value `{value}`")]
596    NonIntegralResult {
597        edge: Box<str>,
598        field: Box<str>,
599        value: Box<str>,
600    },
601    #[error(
602        "reduction `{edge}` target field `{field}` evaluates to negative parameter value `{value}`"
603    )]
604    NegativeResult {
605        edge: Box<str>,
606        field: Box<str>,
607        value: BigRational,
608    },
609    #[error("parameter field `{field}` value `{value}` does not fit u64")]
610    OutputOutOfRange { field: Box<str>, value: BigUint },
611}
612
613#[cfg(test)]
614#[path = "unit_tests/parameters.rs"]
615mod tests;
problemreductions/parameters.rs

problemreductions/
parameters.rs
