From e23a0914179dc2c77403ccde55970e6924deec70 Mon Sep 17 00:00:00 2001 From: "Jip J. Dekker" Date: Thu, 30 May 2024 14:53:08 +1000 Subject: [PATCH] Add int_pow propagator --- crates/huub/src/model/constraint.rs | 9 + crates/huub/src/model/flatzinc.rs | 15 ++ crates/huub/src/propagator.rs | 1 + crates/huub/src/propagator/int_pow.rs | 360 +++++++++++++++++++++++++ crates/huub/src/solver/engine/queue.rs | 1 - 5 files changed, 385 insertions(+), 1 deletion(-) create mode 100644 crates/huub/src/propagator/int_pow.rs diff --git a/crates/huub/src/model/constraint.rs b/crates/huub/src/model/constraint.rs index 538867f5..5350b6eb 100644 --- a/crates/huub/src/model/constraint.rs +++ b/crates/huub/src/model/constraint.rs @@ -21,6 +21,7 @@ use crate::{ array_var_int_element::ArrayVarIntElementBounds, int_lin_le::{IntLinearLessEqBounds, IntLinearLessEqImpBounds}, int_lin_ne::{IntLinearNotEqImpValue, IntLinearNotEqValue}, + int_pow::IntPowBounds, int_times::IntTimesBounds, }, solver::{ @@ -47,6 +48,7 @@ pub enum Constraint { IntLinNotEq(Vec, IntVal), IntLinNotEqImp(Vec, IntVal, BoolExpr), IntLinNotEqReif(Vec, IntVal, BoolExpr), + IntPow(IntView, IntView, IntView), IntTimes(IntView, IntView, IntView), PropLogic(BoolExpr), SetInReif(IntView, IntSetVal, BoolExpr), @@ -334,6 +336,13 @@ impl Constraint { } Ok(()) } + Constraint::IntPow(base, exponent, res) => { + let base = base.to_arg(ReifContext::Mixed, slv, map); + let exponent = exponent.to_arg(ReifContext::Mixed, slv, map); + let result = res.to_arg(ReifContext::Mixed, slv, map); + slv.add_propagator(IntPowBounds::prepare(base, exponent, result)); + Ok(()) + } Constraint::IntTimes(x, y, z) => { let x = x.to_arg(ReifContext::Mixed, slv, map); let y = y.to_arg(ReifContext::Mixed, slv, map); diff --git a/crates/huub/src/model/flatzinc.rs b/crates/huub/src/model/flatzinc.rs index 414251ff..bdc93ac0 100644 --- a/crates/huub/src/model/flatzinc.rs +++ b/crates/huub/src/model/flatzinc.rs @@ -577,6 +577,21 @@ impl Model { }); } } + "int_pow" => { + if let [base, exponent, res] = c.args.as_slice() { + let base = arg_int(fzn, &mut prb, &mut map, base)?; + let exponent = arg_int(fzn, &mut prb, &mut map, exponent)?; + let res = arg_int(fzn, &mut prb, &mut map, res)?; + + prb += Constraint::IntPow(base, exponent, res); + } else { + return Err(FlatZincError::InvalidNumArgs { + name: "int_pow", + found: c.args.len(), + expected: 3, + }); + } + } "int_times" => { if let [x, y, z] = c.args.as_slice() { let a = arg_int(fzn, &mut prb, &mut map, x)?; diff --git a/crates/huub/src/propagator.rs b/crates/huub/src/propagator.rs index f299b2ef..fb0c13d9 100644 --- a/crates/huub/src/propagator.rs +++ b/crates/huub/src/propagator.rs @@ -5,6 +5,7 @@ pub(crate) mod conflict; pub(crate) mod int_event; pub(crate) mod int_lin_le; pub(crate) mod int_lin_ne; +pub(crate) mod int_pow; pub(crate) mod int_times; pub(crate) mod reason; diff --git a/crates/huub/src/propagator/int_pow.rs b/crates/huub/src/propagator/int_pow.rs new file mode 100644 index 00000000..87cbb2b1 --- /dev/null +++ b/crates/huub/src/propagator/int_pow.rs @@ -0,0 +1,360 @@ +use crate::{ + actions::{ + explanation::ExplanationActions, initialization::InitializationActions, + trailing::TrailingActions, + }, + propagator::{ + conflict::Conflict, int_event::IntEvent, reason::ReasonBuilder, PropagationActions, + Propagator, + }, + solver::{ + engine::queue::PriorityLevel, + poster::{BoxedPropagator, Poster}, + }, + IntVal, IntView, LitMeaning, +}; + +#[derive(Debug, Clone, PartialEq, Eq, Hash)] +/// Bounds propagator for the constraint `result = base^exponent`. +pub(crate) struct IntPowBounds { + /// The base in the exponentiation + base: IntView, + /// The exponent in the exponentiation + exponent: IntView, + /// The result of exponentiation + result: IntView, +} + +impl IntPowBounds { + pub(crate) fn prepare(base: IntView, exponent: IntView, result: IntView) -> impl Poster { + IntPowBoundsPoster { + base, + exponent, + result, + } + } + + fn propagate_result(&mut self, actions: &mut P) -> Result<(), Conflict> { + let (base_lb, base_ub) = actions.get_int_bounds(self.base); + let (exp_lb, exp_ub) = actions.get_int_bounds(self.exponent); + let exp_largest_even = if exp_ub % 2 == 0 || exp_lb == exp_ub { + exp_ub + } else { + exp_ub - 1 + }; + let exp_largest_uneven = if exp_ub % 2 == 1 || exp_lb == exp_ub { + exp_ub + } else { + exp_ub - 1 + }; + + let base_lb_lit = actions.get_int_lower_bound_lit(self.base); + let base_ub_lit = actions.get_int_upper_bound_lit(self.base); + let exp_lb_lit = actions.get_int_lower_bound_lit(self.exponent); + let exp_ub_lit = actions.get_int_upper_bound_lit(self.exponent); + + let min: IntVal = vec![ + pow(base_lb, exp_lb), // base and exp always both positive + pow(base_lb, exp_largest_uneven), // base maybe negative + if (base_lb..=base_ub).contains(&0) { + Some(0) + } else { + None + }, + ] + .into_iter() + .flatten() + .min() + .unwrap(); + + eprintln!("base {base_lb}..{base_ub} exp {exp_lb}..{exp_ub} -> res min {min}"); + actions.set_int_lower_bound( + self.result, + min, + &ReasonBuilder::Eager(vec![base_lb_lit, base_ub_lit, exp_lb_lit, exp_ub_lit]), + )?; + + let max: IntVal = vec![ + pow(base_ub, exp_ub), // base and exp have positive upper bounds + pow(base_lb, exp_largest_even), // base maybe negative + if (base_lb..=base_ub).contains(&0) { + Some(0) + } else { + None + }, + ] + .into_iter() + .flatten() + .max() + .unwrap(); + + eprintln!("base {base_lb}..{base_ub} exp {exp_lb}..{exp_ub} -> res max {max}"); + actions.set_int_upper_bound( + self.result, + max, + &ReasonBuilder::Eager(vec![base_lb_lit, base_ub_lit, exp_lb_lit, exp_ub_lit]), + )?; + + Ok(()) + } + + fn propagate_base(&mut self, actions: &mut P) -> Result<(), Conflict> { + let (base_lb, base_ub) = actions.get_int_bounds(self.base); + let (res_lb, res_ub) = actions.get_int_bounds(self.result); + let (exp_lb, exp_ub) = actions.get_int_bounds(self.exponent); + let exp_pos_even = match exp_lb { + _ if exp_lb % 2 == 1 && exp_lb > 0 => exp_lb + 1, + _ if exp_lb < 0 && exp_ub >= 2 => 2, + _ => exp_lb, + }; + let exp_pos_uneven = match exp_lb { + _ if exp_lb % 2 == 0 && exp_lb > 0 => exp_lb + 1, + _ if exp_lb < 0 && exp_ub >= 1 => 1, + _ => exp_lb, + }; + + let res_lb_lit = actions.get_int_lower_bound_lit(self.result); + let res_ub_lit = actions.get_int_upper_bound_lit(self.result); + let exp_lb_lit = actions.get_int_lower_bound_lit(self.exponent); + let exp_ub_lit = actions.get_int_upper_bound_lit(self.exponent); + + if (exp_lb..=exp_ub).contains(&0) && (res_lb..=res_ub).contains(&1) { + return Ok(()); + } + let res_bnd = res_lb..=res_ub; + if exp_lb < 0 && res_bnd.contains(&0) && res_bnd != (1..=1) { + return Ok(()); + } + + if exp_ub < 0 { + actions.set_int_not_eq(self.base, 0, &ReasonBuilder::Simple(exp_ub_lit))?; + } + + // Propagate lower bound + let mut min = vec![ + (res_lb as f64).powf(1_f64 / (exp_ub as f64)), + (res_ub as f64).powf(1_f64 / (exp_pos_uneven as f64)), + (res_lb as f64).powf(1_f64 / (exp_pos_uneven as f64)), + ] + .into_iter() + .reduce(|a, b| a.min(b)) + .unwrap() + .ceil() as IntVal; + + if min > base_lb { + // Correct possible numerical error + eprintln!("res {res_lb}..{res_ub} exp {exp_lb}..{exp_ub} -> base min {min}"); + if (min - 1 != 0 || exp_lb > 0) + && res_lb <= pow(min - 1, if min < 0 { exp_pos_uneven } else { exp_ub }).unwrap() + { + min -= 1; + } + eprintln!("res {res_lb}..{res_ub} exp {exp_lb}..{exp_ub} -> base min {min}"); + actions.set_int_lower_bound( + self.base, + min, + &ReasonBuilder::Eager(vec![res_lb_lit, res_ub_lit, exp_lb_lit, exp_ub_lit]), + )?; + } + + // Propagate upper bound + let mut max = vec![ + (res_ub as f64).powf(1_f64 / (exp_lb as f64)), + (res_ub as f64).powf(1_f64 / (exp_pos_uneven as f64)), + (res_lb as f64).powf(1_f64 / (exp_pos_even as f64)), + -((res_lb as f64).powf(1_f64 / (exp_pos_even as f64))), + ] + .into_iter() + .reduce(|a, b| a.max(b)) + .unwrap() + .floor() as IntVal; + + if max < base_ub { + // Correct possible numerical error + eprintln!("res {res_lb}..{res_ub} exp {exp_lb}..{exp_ub} -> base max {max}"); + if res_ub >= pow(max + 1, if min < 0 { exp_pos_even } else { exp_lb }).unwrap() { + max += 1; + } + eprintln!("res {res_lb}..{res_ub} exp {exp_lb}..{exp_ub} -> base max {max}"); + actions.set_int_upper_bound( + self.base, + max, + &ReasonBuilder::Eager(vec![res_lb_lit, res_ub_lit, exp_lb_lit, exp_ub_lit]), + )?; + } + Ok(()) + } + + fn propagate_exponent( + &mut self, + actions: &mut P, + ) -> Result<(), Conflict> { + let (base_lb, base_ub) = actions.get_int_bounds(self.base); + let (res_lb, res_ub) = actions.get_int_bounds(self.result); + + if base_ub == base_lb && base_lb == 0 { + let base_val_lit = actions.get_int_val_lit(self.base).unwrap(); + return actions.set_int_lower_bound( + self.exponent, + 0, + &ReasonBuilder::Simple(base_val_lit), + ); + } + if base_lb <= 1 || res_lb <= 1 { + // TODO: It seems there should be propagation possible, but log2() certainly won't work. + return Ok(()); + } + + let (exp_lb, exp_ub) = actions.get_int_bounds(self.exponent); + let res_lb_lit = actions.get_int_lit(self.base, LitMeaning::GreaterEq(1)); + let res_ub_lit = actions.get_int_upper_bound_lit(self.result); + let base_lb_lit = actions.get_int_lit(self.base, LitMeaning::GreaterEq(1)); + let base_ub_lit = actions.get_int_upper_bound_lit(self.base); + + // Propagate lower bound + let mut min = ((res_lb as f64).log2() / (base_ub as f64).log2()).ceil() as IntVal; + if min > exp_lb { + // Correct possible numerical error + if res_lb <= pow(base_lb, min - 1).unwrap() { + min -= 1; + } + eprintln!("base {base_lb}..{base_ub} res {res_lb}..{res_ub} -> exp min {min}"); + actions.set_int_lower_bound( + self.base, + min, + &ReasonBuilder::Eager(vec![res_lb_lit, res_ub_lit, base_lb_lit, base_ub_lit]), + )?; + } + + // Propagate upper bound + let mut max = ((res_ub as f64).log2() / (base_lb as f64).log2()).floor() as IntVal; + if max < exp_ub { + // Correct possible numerical error + if res_ub <= pow(base_ub, max + 1).unwrap() { + max += 1; + } + eprintln!("base {base_lb}..{base_ub} res {res_lb}..{res_ub} -> exp max {max}"); + actions.set_int_upper_bound( + self.base, + max, + &ReasonBuilder::Eager(vec![res_lb_lit, res_ub_lit, base_lb_lit, base_ub_lit]), + )?; + } + + Ok(()) + } +} + +impl Propagator for IntPowBounds +where + P: PropagationActions, + E: ExplanationActions, + T: TrailingActions, +{ + fn notify_event(&mut self, _: u32, _: &IntEvent, _: &mut T) -> bool { + true + } + + fn queue_priority_level(&self) -> PriorityLevel { + PriorityLevel::Highest + } + + fn notify_backtrack(&mut self, _new_level: usize) {} + + #[tracing::instrument(name = "int_pow", level = "trace", skip(self, actions))] + fn propagate(&mut self, actions: &mut P) -> Result<(), Conflict> { + self.propagate_result(actions)?; + self.propagate_base(actions)?; + self.propagate_exponent(actions)?; + eprintln!("-------------------------------------"); + + Ok(()) + } +} + +struct IntPowBoundsPoster { + base: IntView, + exponent: IntView, + result: IntView, +} +impl Poster for IntPowBoundsPoster { + fn post(self, actions: &mut I) -> (BoxedPropagator, bool) { + actions.subscribe_int(self.base, IntEvent::Bounds, 1); + actions.subscribe_int(self.exponent, IntEvent::Bounds, 2); + actions.subscribe_int(self.result, IntEvent::Bounds, 3); + ( + Box::new(IntPowBounds { + base: self.base, + exponent: self.exponent, + result: self.result, + }), + false, + ) + } +} + +fn pow(base: IntVal, exponent: IntVal) -> Option { + Some(match exponent { + 0 => 1, + 1 => base, + exp if exp < 0 => match base { + 0 => return None, + 1 => 1, + -1 if exp % 2 == 0 => 1, + -1 => -1, + _ => 0, + }, + _ => { + let mut result = 1; + for _ in 0..exponent { + result *= base; + } + result + } + }) +} + +#[cfg(test)] +mod tests { + use expect_test::expect; + use pindakaas::{solver::cadical::Cadical, Cnf}; + + use crate::{propagator::int_pow::IntPowBounds, solver::engine::int_var::IntVar, Solver}; + + #[test] + fn test_int_pow_sat() { + let mut slv: Solver = Cnf::default().into(); + let a = IntVar::new_in(&mut slv, (-2..=2).into(), true); + let b = IntVar::new_in(&mut slv, (-2..=2).into(), true); + let c = IntVar::new_in(&mut slv, (-2..=4).into(), true); + + slv.add_propagator(IntPowBounds::prepare(a, b, c)); + slv.expect_solutions( + &[a, b, c], + expect![[r#" + -2, -2, 0 + -2, -1, 0 + -2, 0, 1 + -2, 1, -2 + -2, 2, 4 + -1, -2, 1 + -1, -1, -1 + -1, 0, 1 + -1, 1, -1 + -1, 2, 1 + 0, 0, 1 + 0, 1, 0 + 0, 2, 0 + 1, -2, 1 + 1, -1, 1 + 1, 0, 1 + 1, 1, 1 + 1, 2, 1 + 2, -2, 0 + 2, -1, 0 + 2, 0, 1 + 2, 1, 2 + 2, 2, 4"#]], + ); + } +} diff --git a/crates/huub/src/solver/engine/queue.rs b/crates/huub/src/solver/engine/queue.rs index 5d0a5b15..7e743621 100644 --- a/crates/huub/src/solver/engine/queue.rs +++ b/crates/huub/src/solver/engine/queue.rs @@ -51,7 +51,6 @@ pub(crate) enum PriorityLevel { /// A high level of priority, all apart from one normal priority level are /// less important. High, - #[allow(dead_code)] // TODO /// The highest normal priority level, this priority level is the most /// important normal level of priority. Highest,