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

problemreductions/
parameters.rs
