Alex Mikhalev
6 years ago
3 changed files with 901 additions and 903 deletions
@ -1,903 +0,0 @@
@@ -1,903 +0,0 @@
|
||||
pub type Scalar = f64; |
||||
|
||||
pub type Vec2 = nalgebra::Vector2<Scalar>; |
||||
pub type Point2 = nalgebra::Point2<Scalar>; |
||||
|
||||
pub type Rot2 = nalgebra::UnitComplex<Scalar>; |
||||
|
||||
pub trait Region<T> { |
||||
fn full() -> Self; |
||||
fn singleton(value: T) -> Self; |
||||
|
||||
fn nearest(&self, value: &T) -> Option<T>; |
||||
fn contains(&self, value: &T) -> bool; |
||||
} |
||||
|
||||
#[derive(Clone, Debug)] |
||||
pub enum Region1 { |
||||
Empty, |
||||
Singleton(Scalar), |
||||
Range(Scalar, Scalar), |
||||
Union(Box<Region1>, Box<Region1>), |
||||
Full, |
||||
} |
||||
|
||||
impl Region<Scalar> for Region1 { |
||||
fn full() -> Self { |
||||
Region1::Full |
||||
} |
||||
|
||||
fn singleton(value: Scalar) -> Self { |
||||
Region1::Singleton(value) |
||||
} |
||||
|
||||
fn contains(&self, n: &Scalar) -> bool { |
||||
use Region1::*; |
||||
match self { |
||||
Empty => false, |
||||
Singleton(n1) => relative_eq!(n1, n), |
||||
Range(l, u) => *l <= *n && *n <= *u, |
||||
Union(r1, r2) => r1.contains(n) || r2.contains(n), |
||||
Full => true, |
||||
} |
||||
} |
||||
|
||||
fn nearest(&self, s: &Scalar) -> Option<Scalar> { |
||||
use Region1::*; |
||||
match self { |
||||
Empty => None, |
||||
Full => Some(*s), |
||||
Singleton(n) => Some(*n), |
||||
Range(l, u) => match (l < s, s < u) { |
||||
(true, true) => Some(*s), |
||||
(true, false) => Some(*u), |
||||
(false, true) => Some(*l), |
||||
_ => None, |
||||
}, |
||||
Union(r1, r2) => { |
||||
let distance = |a: Scalar, b: Scalar| (a - b).abs(); |
||||
match (r1.nearest(s), r2.nearest(s)) { |
||||
(None, None) => None, |
||||
(Some(n), None) | (None, Some(n)) => Some(n), |
||||
(Some(n1), Some(n2)) => Some({ |
||||
if distance(*s, n1) <= distance(*s, n2) { |
||||
n1 |
||||
} else { |
||||
n2 |
||||
} |
||||
}), |
||||
} |
||||
} |
||||
} |
||||
} |
||||
} |
||||
|
||||
// line starting at start, point at angle dir, with range extent
|
||||
// ie. start + (cos dir, sin dir) * t for t in extent
|
||||
#[derive(Clone, Debug)] |
||||
pub struct Line2 { |
||||
start: Point2, |
||||
dir: Rot2, |
||||
extent: Region1, |
||||
} |
||||
|
||||
impl Line2 { |
||||
pub fn new(start: Point2, dir: Rot2, extent: Region1) -> Self { |
||||
Self { start, dir, extent } |
||||
} |
||||
|
||||
pub fn evaluate(&self, t: Scalar) -> Point2 { |
||||
self.start + self.dir * Vec2::new(t, 0.) |
||||
} |
||||
|
||||
pub fn nearest(&self, p: &Point2) -> Point2 { |
||||
// rotate angle 90 degrees
|
||||
let perp_dir = self.dir * Rot2::from_cos_sin_unchecked(0., 1.); |
||||
let perp = Line2::new(*p, perp_dir, Region1::Full); |
||||
if let Region2::Singleton(np) = self.intersect(&perp) { |
||||
np |
||||
} else { |
||||
panic!("Line2::nearest not found!"); |
||||
} |
||||
} |
||||
|
||||
pub fn intersect(&self, other: &Line2) -> Region2 { |
||||
// if the two lines are parallel...
|
||||
let dirs = self.dir / other.dir; |
||||
if relative_eq!(dirs.sin_angle(), 0.) { |
||||
let starts = self.dir.to_rotation_matrix().inverse() * (other.start - self.start); |
||||
return if relative_eq!(starts.y, 0.) { |
||||
// and they are colinear
|
||||
Region2::Line(self.clone()) |
||||
} else { |
||||
// they are parallel and never intersect
|
||||
Region2::Empty |
||||
}; |
||||
} |
||||
// TODO: respect extent
|
||||
let (a, b) = (self, other); |
||||
let (a_0, a_v, b_0, b_v) = (a.start, a.dir, b.start, b.dir); |
||||
let (a_c, a_s, b_c, b_s) = ( |
||||
a_v.cos_angle(), |
||||
a_v.sin_angle(), |
||||
b_v.cos_angle(), |
||||
b_v.sin_angle(), |
||||
); |
||||
let t_b = (a_0.x * a_s - a_0.y * a_c + a_0.x * a_s + b_0.y * a_c) / (a_s * b_c - a_c * b_s); |
||||
Region2::Singleton(b.evaluate(t_b)) |
||||
} |
||||
} |
||||
|
||||
#[derive(Clone, Debug)] |
||||
pub enum Region2 { |
||||
Empty, |
||||
// single point at 0
|
||||
Singleton(Point2), |
||||
Line(Line2), |
||||
#[allow(dead_code)] |
||||
Union(Box<Region2>, Box<Region2>), |
||||
Full, |
||||
} |
||||
|
||||
impl Region<Point2> for Region2 { |
||||
fn full() -> Self { |
||||
Region2::Full |
||||
} |
||||
|
||||
fn singleton(value: Point2) -> Self { |
||||
Region2::Singleton(value) |
||||
} |
||||
|
||||
fn contains(&self, p: &Point2) -> bool { |
||||
self.nearest(p).map_or(false, |n| relative_eq!(n, p)) |
||||
} |
||||
|
||||
fn nearest(&self, p: &Point2) -> Option<Point2> { |
||||
use Region2::*; |
||||
match self { |
||||
Empty => None, |
||||
Full => Some(*p), |
||||
Singleton(n) => Some(*n), |
||||
Line(line) => Some(line.nearest(p)), |
||||
Union(r1, r2) => { |
||||
use nalgebra::distance; |
||||
match (r1.nearest(p), r2.nearest(p)) { |
||||
(None, None) => None, |
||||
(Some(n), None) | (None, Some(n)) => Some(n), |
||||
(Some(n1), Some(n2)) => Some({ |
||||
if distance(p, &n1) <= distance(p, &n2) { |
||||
n1 |
||||
} else { |
||||
n2 |
||||
} |
||||
}), |
||||
} |
||||
} |
||||
} |
||||
} |
||||
} |
||||
|
||||
impl Region2 { |
||||
pub fn union(r1: Region2, r2: Region2) -> Region2 { |
||||
use Region2::*; |
||||
match (r1, r2) { |
||||
(Empty, r) | (r, Empty) => r, |
||||
(Full, _) | (_, Full) => Full, |
||||
(r1, r2) => Union(Box::new(r1), Box::new(r2)), |
||||
} |
||||
} |
||||
|
||||
pub fn intersect(&self, other: &Region2) -> Region2 { |
||||
use Region2::*; |
||||
match (self, other) { |
||||
(Empty, _) | (_, Empty) => Empty, |
||||
(Full, r) | (r, Full) => r.clone(), |
||||
(Singleton(n1), Singleton(n2)) => { |
||||
if n1 == n2 { |
||||
Singleton(*n1) |
||||
} else { |
||||
Empty |
||||
} |
||||
} |
||||
(Singleton(n), o) | (o, Singleton(n)) => { |
||||
if o.contains(n) { |
||||
Singleton(*n) |
||||
} else { |
||||
Empty |
||||
} |
||||
} |
||||
(Line(l1), Line(l2)) => l1.intersect(l2), |
||||
(Union(un1, un2), o) | (o, Union(un1, un2)) => { |
||||
Self::union(un1.intersect(o), un2.intersect(o)) |
||||
} |
||||
} |
||||
} |
||||
} |
||||
|
||||
pub mod solve { |
||||
|
||||
use std::collections::{BTreeMap, BTreeSet}; |
||||
use std::fmt; |
||||
use std::iter::FromIterator; |
||||
|
||||
use crate::math::Scalar; |
||||
|
||||
// an unknown variable with an id
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)] |
||||
struct Unknown(i64); |
||||
|
||||
type UnknownSet = BTreeSet<Unknown>; |
||||
|
||||
trait Unknowns { |
||||
fn unknowns(&self) -> UnknownSet; |
||||
fn has_unknowns(&self) -> bool; |
||||
fn has_unknown(&self, u: Unknown) -> bool; |
||||
} |
||||
|
||||
impl Unknowns for Scalar { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
UnknownSet::new() |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
false |
||||
} |
||||
fn has_unknown(&self, _: Unknown) -> bool { |
||||
false |
||||
} |
||||
} |
||||
|
||||
impl Unknowns for Unknown { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
FromIterator::from_iter(Some(*self)) |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
true |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
*self == u |
||||
} |
||||
} |
||||
|
||||
impl fmt::Display for Unknown { |
||||
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { |
||||
write!(f, "u{}", self.0) |
||||
} |
||||
} |
||||
|
||||
#[derive(Clone, Debug, PartialEq)] |
||||
enum Expr { |
||||
Unkn(Unknown), |
||||
Const(Scalar), |
||||
Sum(Exprs), |
||||
Neg(Box<Expr>), |
||||
Product(Exprs), |
||||
Div(Box<Expr>, Box<Expr>), |
||||
} |
||||
|
||||
type Exprs = Vec<Expr>; |
||||
|
||||
impl Unknowns for Exprs { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
self.iter().flat_map(|e: &Expr| e.unknowns()).collect() |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
self.iter().any(|e: &Expr| e.has_unknowns()) |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
self.iter().any(|e: &Expr| e.has_unknown(u)) |
||||
} |
||||
} |
||||
|
||||
fn write_separated_exprs(es: &Exprs, f: &mut fmt::Formatter, sep: &str) -> fmt::Result { |
||||
let mut is_first = true; |
||||
for e in es { |
||||
if is_first { |
||||
is_first = false; |
||||
} else { |
||||
write!(f, "{}", sep)? |
||||
} |
||||
write!(f, "({})", e)? |
||||
} |
||||
Ok(()) |
||||
} |
||||
|
||||
fn remove_common_terms(l: &mut Vec<Expr>, r: &mut Vec<Expr>) -> Vec<Expr> { |
||||
let common: Vec<_> = l.drain_filter(|e| r.contains(e)).collect(); |
||||
common.iter().for_each(|e| { |
||||
r.remove_item(e); |
||||
}); |
||||
common |
||||
} |
||||
|
||||
fn remove_term(terms: &mut Vec<Expr>, term: &Expr) -> Option<Expr> { |
||||
terms.remove_item(term) |
||||
} |
||||
|
||||
fn sum_fold(l: Expr, r: Expr) -> Expr { |
||||
use itertools::Itertools; |
||||
use Expr::*; |
||||
match (l, r) { |
||||
(Const(lc), Const(rc)) => Const(lc + rc), |
||||
(Const(c), o) | (o, Const(c)) if relative_eq!(c, 0.) => o, |
||||
(Product(mut l), Product(mut r)) => { |
||||
let comm = remove_common_terms(&mut l, &mut r); |
||||
Expr::new_product(Sum(comm), Expr::new_sum(Product(l), Product(r))).simplify() |
||||
} |
||||
(Product(mut l), r) | (r, Product(mut l)) => { |
||||
let comm = remove_term(&mut l, &r); |
||||
match comm { |
||||
Some(_) => { |
||||
Expr::new_product(r, Expr::new_sum(Product(l), Const(1.))).simplify() |
||||
} |
||||
None => Expr::new_sum(Product(l), r), |
||||
} |
||||
} |
||||
(l, r) => Expr::new_sum(l, r), |
||||
} |
||||
} |
||||
|
||||
fn group_sum(es: Exprs) -> Exprs { |
||||
use Expr::*; |
||||
let mut common: BTreeMap<UnknownSet, Expr> = BTreeMap::new(); |
||||
for e in es { |
||||
let unkns = e.unknowns(); |
||||
match common.get_mut(&unkns) { |
||||
None => { |
||||
match e { |
||||
Const(c) if relative_eq!(c, 0.) => (), |
||||
e => { |
||||
common.insert(unkns, e); |
||||
} |
||||
}; |
||||
} |
||||
Some(existing) => { |
||||
match existing { |
||||
Sum(ref mut es) => { |
||||
// already failed at merging, so just add it to the list
|
||||
es.push(e); |
||||
} |
||||
other => { |
||||
*other = sum_fold(other.clone(), e); |
||||
} |
||||
}; |
||||
} |
||||
}; |
||||
} |
||||
common.into_iter().map(|(_, v)| v).collect() |
||||
} |
||||
|
||||
fn product_fold(l: Expr, r: Expr) -> Expr { |
||||
use itertools::Itertools; |
||||
use Expr::*; |
||||
match (l, r) { |
||||
(Const(lc), Const(rc)) => Const(lc * rc), |
||||
(Const(c), o) | (o, Const(c)) if relative_eq!(c, 1.) => o, |
||||
(Div(num, den), mul) | (mul, Div(num, den)) => { |
||||
if mul == *den { |
||||
*num |
||||
} else { |
||||
Expr::Div(Box::new(Expr::Product(vec![*num, mul])), den).simplify() |
||||
} |
||||
} |
||||
(l, r) => Expr::new_product(l, r), |
||||
} |
||||
} |
||||
|
||||
fn group_product(es: Exprs) -> Exprs { |
||||
use Expr::*; |
||||
let mut common: BTreeMap<UnknownSet, Expr> = BTreeMap::new(); |
||||
for e in es { |
||||
let unkns = e.unknowns(); |
||||
match common.get_mut(&unkns) { |
||||
None => { |
||||
match e { |
||||
Const(c) if relative_eq!(c, 1.) => (), |
||||
e => { |
||||
common.insert(unkns, e); |
||||
} |
||||
}; |
||||
} |
||||
Some(existing) => { |
||||
match existing { |
||||
Sum(ref mut es) => { |
||||
// already failed at merging, so just add it to the list
|
||||
es.push(e); |
||||
} |
||||
other => *other = product_fold(other.clone(), e), |
||||
}; |
||||
} |
||||
}; |
||||
} |
||||
common.into_iter().map(|(_, v)| v).collect() |
||||
} |
||||
|
||||
fn distribute_product_sums(mut es: Exprs) -> Expr { |
||||
trace!("distribute_product_sums: {}", Product(es.clone())); |
||||
use itertools::Itertools; |
||||
use Expr::*; |
||||
let sums = es |
||||
.drain_filter(|e| match e { |
||||
Sum(_) => true, |
||||
_ => false, |
||||
}) |
||||
.map(|e| { |
||||
trace!("sum in product: {}", e); |
||||
match e { |
||||
Sum(es) => es, |
||||
_ => unreachable!(), |
||||
} |
||||
}); |
||||
let products: Vec<_> = sums.multi_cartesian_product().collect(); |
||||
if products.is_empty() { |
||||
trace!("no sums to distribute"); |
||||
return Product(es); |
||||
} |
||||
let sums = products |
||||
.into_iter() |
||||
.map(|mut prod| { |
||||
prod.extend(es.clone()); |
||||
trace!("prod: {}", Product(prod.clone())); |
||||
Product(prod) |
||||
}) |
||||
.collect(); |
||||
Sum(sums) |
||||
} |
||||
|
||||
impl Unknowns for Expr { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
use Expr::*; |
||||
match self { |
||||
Unkn(u) => u.unknowns(), |
||||
Const(_) => UnknownSet::default(), |
||||
Sum(es) | Product(es) => es.unknowns(), |
||||
Div(l, r) => l.unknowns().union(&r.unknowns()).cloned().collect(), |
||||
Neg(e) => e.unknowns(), |
||||
} |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
use Expr::*; |
||||
match self { |
||||
Unkn(u) => u.has_unknowns(), |
||||
Const(_) => false, |
||||
Sum(es) | Product(es) => es.has_unknowns(), |
||||
Div(l, r) => l.has_unknowns() || r.has_unknowns(), |
||||
Neg(e) => e.has_unknowns(), |
||||
} |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
use Expr::*; |
||||
match self { |
||||
Unkn(u1) => u1.has_unknown(u), |
||||
Const(_) => false, |
||||
Sum(es) | Product(es) => es.has_unknown(u), |
||||
Div(l, r) => l.has_unknown(u) || r.has_unknown(u), |
||||
Neg(e) => e.has_unknown(u), |
||||
} |
||||
} |
||||
} |
||||
|
||||
impl Expr { |
||||
fn new_sum(e1: Expr, e2: Expr) -> Expr { |
||||
Expr::Sum(vec![e1, e2]) |
||||
} |
||||
fn new_product(e1: Expr, e2: Expr) -> Expr { |
||||
Expr::Product(vec![e1, e2]) |
||||
} |
||||
fn new_neg(e1: Expr) -> Expr { |
||||
Expr::Neg(Box::new(e1)) |
||||
} |
||||
fn new_div(num: Expr, den: Expr) -> Expr { |
||||
Expr::Div(Box::new(num), Box::new(den)) |
||||
} |
||||
fn new_minus(e1: Expr, e2: Expr) -> Expr { |
||||
Expr::Sum(vec![e1, Expr::new_neg(e2)]) |
||||
} |
||||
fn new_inv(den: Expr) -> Expr { |
||||
Expr::new_div(Expr::Const(1.), den) |
||||
} |
||||
|
||||
fn is_zero(&self) -> bool { |
||||
use Expr::*; |
||||
match self { |
||||
Const(c) => relative_eq!(*c, 0.), |
||||
_ => false, |
||||
} |
||||
} |
||||
|
||||
fn is_one(&self) -> bool { |
||||
use Expr::*; |
||||
match self { |
||||
Const(c) => relative_eq!(*c, 1.), |
||||
_ => false, |
||||
} |
||||
} |
||||
|
||||
fn simplify(self) -> Expr { |
||||
use Expr::*; |
||||
match self { |
||||
Sum(es) => { |
||||
let mut new_es: Vec<_> = es |
||||
.into_iter() |
||||
.map(|e| e.simplify()) |
||||
.flat_map(|e| match e { |
||||
Sum(more_es) => more_es, |
||||
other => vec![other], |
||||
}) |
||||
.collect(); |
||||
let pre_new_es = new_es.clone(); |
||||
new_es = group_sum(new_es); |
||||
trace!( |
||||
"simplify sum {} => {}", |
||||
Sum(pre_new_es), |
||||
Sum(new_es.clone()) |
||||
); |
||||
|
||||
match new_es.len() { |
||||
0 => Const(0.), // none
|
||||
1 => new_es.into_iter().next().unwrap(), // one
|
||||
_ => Sum(new_es), // many
|
||||
} |
||||
} |
||||
Product(es) => { |
||||
let new_es: Vec<_> = es |
||||
.into_iter() |
||||
.map(|e| e.simplify()) |
||||
.flat_map(|e| match e { |
||||
Product(more_es) => more_es, |
||||
other => vec![other], |
||||
}) |
||||
.collect(); |
||||
let pre_new_es = new_es.clone(); |
||||
let new_es = group_product(new_es); |
||||
trace!( |
||||
"simplify product {} => {}", |
||||
Product(pre_new_es), |
||||
Product(new_es.clone()) |
||||
); |
||||
match new_es.len() { |
||||
0 => Const(1.), // none
|
||||
1 => new_es.into_iter().next().unwrap(), // one
|
||||
_ => Product(new_es), // many
|
||||
} |
||||
} |
||||
Neg(mut v) => { |
||||
*v = v.simplify(); |
||||
trace!("simplify neg {}", Neg(v.clone())); |
||||
match v { |
||||
box Const(c) => Const(-c), |
||||
box Neg(v) => *v, |
||||
e => Product(vec![Const(-1.), *e]), |
||||
} |
||||
} |
||||
Div(mut num, mut den) => { |
||||
*num = num.simplify(); |
||||
*den = den.simplify(); |
||||
trace!("simplify div {}", Div(num.clone(), den.clone())); |
||||
match (num, den) { |
||||
(box Const(num), box Const(den)) => Const(num / den), |
||||
(num, box Const(den)) => { |
||||
if relative_eq!(den, 1.) { |
||||
*num |
||||
} else { |
||||
Expr::new_product(*num, Const(1. / den)) |
||||
} |
||||
} |
||||
(num, box Div(dennum, denden)) => { |
||||
Div(Box::new(Product(vec![*num, *denden])), dennum).simplify() |
||||
} |
||||
(box Product(mut es), box den) => match es.remove_item(&den) { |
||||
Some(_) => Product(es), |
||||
None => Expr::new_div(Product(es), den), |
||||
}, |
||||
(num, den) => { |
||||
if num == den { |
||||
Expr::Const(1.) |
||||
} else { |
||||
Div(num, den) |
||||
} |
||||
} |
||||
} |
||||
} |
||||
e => e, |
||||
} |
||||
} |
||||
|
||||
fn distribute(self) -> Expr { |
||||
use Expr::*; |
||||
trace!("distribute {}", self); |
||||
match self { |
||||
Sum(mut es) => { |
||||
for e in &mut es { |
||||
*e = e.clone().distribute(); |
||||
} |
||||
Sum(es) |
||||
} |
||||
Product(es) => distribute_product_sums(es), |
||||
Div(mut num, mut den) => { |
||||
*num = num.distribute(); |
||||
*den = den.distribute(); |
||||
match (num, den) { |
||||
(box Sum(es), box den) => Sum(es |
||||
.into_iter() |
||||
.map(|e| Expr::new_div(e, den.clone())) |
||||
.collect()), |
||||
(mut num, mut den) => Div(num, den), |
||||
} |
||||
} |
||||
Neg(v) => match v { |
||||
// box Sum(mut l, mut r) => {
|
||||
// *l = Neg(l.clone()).distribute();
|
||||
// *r = Neg(r.clone()).distribute();
|
||||
// Sum(l, r)
|
||||
// }
|
||||
// box Product(mut l, r) => {
|
||||
// *l = Neg(l.clone()).distribute();
|
||||
// Product(l, r)
|
||||
// }
|
||||
box Neg(v) => v.distribute(), |
||||
box Div(mut num, mut den) => { |
||||
*num = Neg(num.clone()).distribute(); |
||||
*den = Neg(den.clone()).distribute(); |
||||
Div(num, den) |
||||
} |
||||
e => Neg(e), |
||||
}, |
||||
e => e, |
||||
} |
||||
} |
||||
} |
||||
|
||||
impl fmt::Display for Expr { |
||||
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { |
||||
use Expr::*; |
||||
match self { |
||||
Unkn(u) => write!(f, "{}", u), |
||||
Const(c) => write!(f, "{}", c), |
||||
Sum(es) => write_separated_exprs(es, f, " + "), |
||||
Product(es) => write_separated_exprs(es, f, " * "), |
||||
Div(num, den) => write!(f, "({}) / ({})", num, den), |
||||
Neg(e) => write!(f, "-({})", e), |
||||
} |
||||
} |
||||
} |
||||
|
||||
#[derive(Clone, Debug, PartialEq)] |
||||
struct Eqn(Expr, Expr); |
||||
|
||||
impl Unknowns for Eqn { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
self.0 |
||||
.unknowns() |
||||
.union(&self.1.unknowns()) |
||||
.cloned() |
||||
.collect() |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
self.0.has_unknowns() || self.1.has_unknowns() |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
self.0.has_unknown(u) || self.1.has_unknown(u) |
||||
} |
||||
} |
||||
|
||||
fn ord_by_unkn(a: Expr, b: Expr, u: Unknown) -> Option<(Expr, Expr)> { |
||||
if a.has_unknown(u) { |
||||
Some((a, b)) |
||||
} else if b.has_unknown(u) { |
||||
Some((b, a)) |
||||
} else { |
||||
None |
||||
} |
||||
} |
||||
|
||||
impl Eqn { |
||||
fn simplify(self) -> Eqn { |
||||
Eqn(self.0.simplify(), self.1.simplify()) |
||||
} |
||||
|
||||
fn solve(&self, for_u: Unknown) -> Option<Expr> { |
||||
use Expr::*; |
||||
if !self.has_unknown(for_u) { |
||||
return None; |
||||
} |
||||
let (l, r) = ( |
||||
self.0 |
||||
.clone() /*.distribute()*/ |
||||
.simplify(), |
||||
self.1 |
||||
.clone() /*.distribute()*/ |
||||
.simplify(), |
||||
); |
||||
let (mut l, mut r) = ord_by_unkn(l, r, for_u)?; |
||||
loop { |
||||
trace!("solve: {} == {}", l, r); |
||||
let (new_l, new_r): (Expr, Expr) = match l { |
||||
Unkn(u) => return if u == for_u { Some(r.simplify()) } else { None }, |
||||
Sum(es) => { |
||||
let (us, not_us): (Vec<_>, Vec<_>) = |
||||
es.into_iter().partition(|e| e.has_unknown(for_u)); |
||||
if us.len() != 1 { |
||||
return None; |
||||
} |
||||
( |
||||
Sum(us).simplify(), |
||||
Expr::new_minus(r, Sum(not_us)).simplify(), |
||||
) |
||||
} |
||||
Product(es) => { |
||||
let (us, not_us): (Vec<_>, Vec<_>) = |
||||
es.into_iter().partition(|e| e.has_unknown(for_u)); |
||||
if us.len() != 1 { |
||||
return None; |
||||
} |
||||
( |
||||
Product(us).simplify(), |
||||
Expr::new_div(r, Product(not_us)).simplify(), |
||||
) |
||||
} |
||||
Neg(v) => (*v, Expr::new_neg(r)), |
||||
Div(num, den) => { |
||||
let (nu, du) = (num.has_unknown(for_u), den.has_unknown(for_u)); |
||||
match (nu, du) { |
||||
(true, false) => (*num, Expr::new_product(r, *den)), |
||||
(false, true) => (Expr::new_product(r, *den), *num), |
||||
(true, true) => return None, // TODO: simplify
|
||||
(false, false) => return None, |
||||
} |
||||
} |
||||
Const(_) => return None, |
||||
_ => return None, |
||||
}; |
||||
l = new_l; |
||||
r = new_r; |
||||
} |
||||
} |
||||
} |
||||
|
||||
impl fmt::Display for Eqn { |
||||
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { |
||||
write!(f, "{} == {}", self.0, self.1) |
||||
} |
||||
} |
||||
|
||||
#[derive(Clone, Debug, PartialEq)] |
||||
struct Eqns(Vec<Eqn>); |
||||
|
||||
impl Unknowns for Eqns { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
self.0.iter().flat_map(|eqn: &Eqn| eqn.unknowns()).collect() |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
self.0.iter().any(|eqn: &Eqn| eqn.has_unknowns()) |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
self.0.iter().any(|eqn: &Eqn| eqn.has_unknown(u)) |
||||
} |
||||
} |
||||
|
||||
#[cfg(test)] |
||||
mod tests { |
||||
use super::*; |
||||
|
||||
#[test] |
||||
fn test_unknowns() { |
||||
let u1 = Unknown(1); |
||||
let u2 = Unknown(2); |
||||
let u3 = Unknown(3); |
||||
assert!(u1.has_unknowns()); |
||||
assert!(u2.has_unknowns()); |
||||
assert!(u1.has_unknown(u1)); |
||||
assert!(!u1.has_unknown(u2)); |
||||
assert!(u1.unknowns().contains(&u1)); |
||||
assert!(!u2.unknowns().contains(&u1)); |
||||
|
||||
let e1 = Expr::new_minus(Expr::Unkn(u1), Expr::Unkn(u2)); |
||||
assert!(e1.has_unknowns()); |
||||
assert!(e1.has_unknown(u1)); |
||||
assert!(e1.has_unknown(u2)); |
||||
assert!(!e1.has_unknown(u3)); |
||||
assert!(e1.unknowns().len() == 2); |
||||
} |
||||
|
||||
fn const_expr(e: Expr) -> Option<Scalar> { |
||||
match e { |
||||
Expr::Const(c) => Some(c), |
||||
_ => None, |
||||
} |
||||
} |
||||
|
||||
#[test] |
||||
fn test_solve() { |
||||
use Expr::*; |
||||
let _ = env_logger::try_init(); |
||||
let u1 = Unknown(1); |
||||
let e1 = Unkn(u1); |
||||
let e2 = Const(1.); |
||||
|
||||
let eqn = Eqn(e1.clone(), e2.clone()); |
||||
assert_eq!(eqn.solve(u1), Some(Const(1.))); |
||||
let eqn = Eqn(e2.clone(), e1.clone()); |
||||
assert_eq!(eqn.solve(u1), Some(Const(1.))); |
||||
|
||||
let e3 = Expr::new_sum(Const(1.), Expr::new_sum(Const(1.), Const(2.))); |
||||
let eqn = Eqn(e1.clone(), e3.clone()); |
||||
assert_eq!(eqn.solve(u1), Some(Const(4.))); |
||||
let e3 = Expr::new_minus(Const(1.), Const(1.)); |
||||
let eqn = Eqn(e1.clone(), e3.clone()); |
||||
assert_eq!(eqn.solve(u1), Some(Const(0.))); |
||||
|
||||
let e1 = Expr::new_div(Const(2.), Expr::new_minus(Const(1.), Const(4.))); |
||||
let e2 = Expr::new_minus(Const(1.), Unkn(u1)); |
||||
let eqn = Eqn(e1, e2); |
||||
info!("eqn: {} => {}", eqn, eqn.clone().simplify()); |
||||
let e = eqn.solve(u1).unwrap(); |
||||
assert!(const_expr(e.clone()).is_some()); |
||||
assert!(relative_eq!(const_expr(e.clone()).unwrap(), 5. / 3.)); |
||||
|
||||
let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); |
||||
let e2 = Expr::new_minus(Expr::new_product(Unkn(u1), Const(2.)), Unkn(u1)); |
||||
info!( |
||||
"e1==e2: {}=={} => {}=={}", |
||||
e1, |
||||
e2, |
||||
e1.clone().simplify(), |
||||
e2.clone().simplify() |
||||
); |
||||
let eqn = Eqn(e1, e2); |
||||
let e = eqn.solve(u1).unwrap(); |
||||
assert!(const_expr(e.clone()).is_some()); |
||||
assert!(relative_eq!(const_expr(e.clone()).unwrap(), -6.)); |
||||
|
||||
let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); |
||||
let e2 = Expr::new_div( |
||||
Expr::new_sum( |
||||
Expr::new_product(Unkn(u1), Const(2.)), |
||||
Expr::new_product(Unkn(u1), Unkn(u1)), |
||||
), |
||||
Unkn(u1), |
||||
); |
||||
info!( |
||||
"{}=={} distrib=> {}=={}", |
||||
e1, |
||||
e2, |
||||
e1.clone().distribute(), |
||||
e2.clone().distribute() |
||||
); |
||||
info!( |
||||
"{}=={} simplify=> {}=={}", |
||||
e1, |
||||
e2, |
||||
e1.clone().distribute().simplify(), |
||||
e2.clone().distribute().simplify() |
||||
); |
||||
let eqn = Eqn(e1, e2); |
||||
let e = eqn.solve(u1).unwrap(); |
||||
assert!(const_expr(e.clone()).is_some()); |
||||
assert!(relative_eq!(const_expr(e.clone()).unwrap(), -8.)); |
||||
|
||||
let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); |
||||
let e2 = Expr::new_div( |
||||
Expr::new_sum( |
||||
Expr::new_product(Unkn(u1), Const(2.)), |
||||
Expr::new_sum( |
||||
Expr::new_sum(Expr::new_product(Unkn(u1), Unkn(u1)), Unkn(u1)), |
||||
Expr::new_minus(Const(2.), Const(1. + 1.)), |
||||
), |
||||
), |
||||
Unkn(u1), |
||||
); |
||||
info!( |
||||
"e1==e2: {}=={} => {}=={}", |
||||
e1, |
||||
e2, |
||||
e1.clone().distribute().simplify(), |
||||
e2.clone().distribute().simplify().simplify() |
||||
); |
||||
let eqn = Eqn(e1, e2); |
||||
let e = eqn.solve(u1).unwrap(); |
||||
assert!(const_expr(e.clone()).is_some()); |
||||
assert!(relative_eq!(const_expr(e.clone()).unwrap(), -9.)); |
||||
} |
||||
|
||||
} |
||||
} |
@ -0,0 +1,217 @@
@@ -0,0 +1,217 @@
|
||||
pub mod solver; |
||||
|
||||
pub type Scalar = f64; |
||||
|
||||
pub type Vec2 = nalgebra::Vector2<Scalar>; |
||||
pub type Point2 = nalgebra::Point2<Scalar>; |
||||
|
||||
pub type Rot2 = nalgebra::UnitComplex<Scalar>; |
||||
|
||||
pub trait Region<T> { |
||||
fn full() -> Self; |
||||
fn singleton(value: T) -> Self; |
||||
|
||||
fn nearest(&self, value: &T) -> Option<T>; |
||||
fn contains(&self, value: &T) -> bool; |
||||
} |
||||
|
||||
#[derive(Clone, Debug)] |
||||
pub enum Region1 { |
||||
Empty, |
||||
Singleton(Scalar), |
||||
Range(Scalar, Scalar), |
||||
Union(Box<Region1>, Box<Region1>), |
||||
Full, |
||||
} |
||||
|
||||
impl Region<Scalar> for Region1 { |
||||
fn full() -> Self { |
||||
Region1::Full |
||||
} |
||||
|
||||
fn singleton(value: Scalar) -> Self { |
||||
Region1::Singleton(value) |
||||
} |
||||
|
||||
fn contains(&self, n: &Scalar) -> bool { |
||||
use Region1::*; |
||||
match self { |
||||
Empty => false, |
||||
Singleton(n1) => relative_eq!(n1, n), |
||||
Range(l, u) => *l <= *n && *n <= *u, |
||||
Union(r1, r2) => r1.contains(n) || r2.contains(n), |
||||
Full => true, |
||||
} |
||||
} |
||||
|
||||
fn nearest(&self, s: &Scalar) -> Option<Scalar> { |
||||
use Region1::*; |
||||
match self { |
||||
Empty => None, |
||||
Full => Some(*s), |
||||
Singleton(n) => Some(*n), |
||||
Range(l, u) => match (l < s, s < u) { |
||||
(true, true) => Some(*s), |
||||
(true, false) => Some(*u), |
||||
(false, true) => Some(*l), |
||||
_ => None, |
||||
}, |
||||
Union(r1, r2) => { |
||||
let distance = |a: Scalar, b: Scalar| (a - b).abs(); |
||||
match (r1.nearest(s), r2.nearest(s)) { |
||||
(None, None) => None, |
||||
(Some(n), None) | (None, Some(n)) => Some(n), |
||||
(Some(n1), Some(n2)) => Some({ |
||||
if distance(*s, n1) <= distance(*s, n2) { |
||||
n1 |
||||
} else { |
||||
n2 |
||||
} |
||||
}), |
||||
} |
||||
} |
||||
} |
||||
} |
||||
} |
||||
|
||||
// line starting at start, point at angle dir, with range extent
|
||||
// ie. start + (cos dir, sin dir) * t for t in extent
|
||||
#[derive(Clone, Debug)] |
||||
pub struct Line2 { |
||||
start: Point2, |
||||
dir: Rot2, |
||||
extent: Region1, |
||||
} |
||||
|
||||
impl Line2 { |
||||
pub fn new(start: Point2, dir: Rot2, extent: Region1) -> Self { |
||||
Self { start, dir, extent } |
||||
} |
||||
|
||||
pub fn evaluate(&self, t: Scalar) -> Point2 { |
||||
self.start + self.dir * Vec2::new(t, 0.) |
||||
} |
||||
|
||||
pub fn nearest(&self, p: &Point2) -> Point2 { |
||||
// rotate angle 90 degrees
|
||||
let perp_dir = self.dir * Rot2::from_cos_sin_unchecked(0., 1.); |
||||
let perp = Line2::new(*p, perp_dir, Region1::Full); |
||||
if let Region2::Singleton(np) = self.intersect(&perp) { |
||||
np |
||||
} else { |
||||
panic!("Line2::nearest not found!"); |
||||
} |
||||
} |
||||
|
||||
pub fn intersect(&self, other: &Line2) -> Region2 { |
||||
// if the two lines are parallel...
|
||||
let dirs = self.dir / other.dir; |
||||
if relative_eq!(dirs.sin_angle(), 0.) { |
||||
let starts = self.dir.to_rotation_matrix().inverse() * (other.start - self.start); |
||||
return if relative_eq!(starts.y, 0.) { |
||||
// and they are colinear
|
||||
Region2::Line(self.clone()) |
||||
} else { |
||||
// they are parallel and never intersect
|
||||
Region2::Empty |
||||
}; |
||||
} |
||||
// TODO: respect extent
|
||||
let (a, b) = (self, other); |
||||
let (a_0, a_v, b_0, b_v) = (a.start, a.dir, b.start, b.dir); |
||||
let (a_c, a_s, b_c, b_s) = ( |
||||
a_v.cos_angle(), |
||||
a_v.sin_angle(), |
||||
b_v.cos_angle(), |
||||
b_v.sin_angle(), |
||||
); |
||||
let t_b = (a_0.x * a_s - a_0.y * a_c + a_0.x * a_s + b_0.y * a_c) / (a_s * b_c - a_c * b_s); |
||||
Region2::Singleton(b.evaluate(t_b)) |
||||
} |
||||
} |
||||
|
||||
#[derive(Clone, Debug)] |
||||
pub enum Region2 { |
||||
Empty, |
||||
// single point at 0
|
||||
Singleton(Point2), |
||||
Line(Line2), |
||||
#[allow(dead_code)] |
||||
Union(Box<Region2>, Box<Region2>), |
||||
Full, |
||||
} |
||||
|
||||
impl Region<Point2> for Region2 { |
||||
fn full() -> Self { |
||||
Region2::Full |
||||
} |
||||
|
||||
fn singleton(value: Point2) -> Self { |
||||
Region2::Singleton(value) |
||||
} |
||||
|
||||
fn contains(&self, p: &Point2) -> bool { |
||||
self.nearest(p).map_or(false, |n| relative_eq!(n, p)) |
||||
} |
||||
|
||||
fn nearest(&self, p: &Point2) -> Option<Point2> { |
||||
use Region2::*; |
||||
match self { |
||||
Empty => None, |
||||
Full => Some(*p), |
||||
Singleton(n) => Some(*n), |
||||
Line(line) => Some(line.nearest(p)), |
||||
Union(r1, r2) => { |
||||
use nalgebra::distance; |
||||
match (r1.nearest(p), r2.nearest(p)) { |
||||
(None, None) => None, |
||||
(Some(n), None) | (None, Some(n)) => Some(n), |
||||
(Some(n1), Some(n2)) => Some({ |
||||
if distance(p, &n1) <= distance(p, &n2) { |
||||
n1 |
||||
} else { |
||||
n2 |
||||
} |
||||
}), |
||||
} |
||||
} |
||||
} |
||||
} |
||||
} |
||||
|
||||
impl Region2 { |
||||
pub fn union(r1: Region2, r2: Region2) -> Region2 { |
||||
use Region2::*; |
||||
match (r1, r2) { |
||||
(Empty, r) | (r, Empty) => r, |
||||
(Full, _) | (_, Full) => Full, |
||||
(r1, r2) => Union(Box::new(r1), Box::new(r2)), |
||||
} |
||||
} |
||||
|
||||
pub fn intersect(&self, other: &Region2) -> Region2 { |
||||
use Region2::*; |
||||
match (self, other) { |
||||
(Empty, _) | (_, Empty) => Empty, |
||||
(Full, r) | (r, Full) => r.clone(), |
||||
(Singleton(n1), Singleton(n2)) => { |
||||
if n1 == n2 { |
||||
Singleton(*n1) |
||||
} else { |
||||
Empty |
||||
} |
||||
} |
||||
(Singleton(n), o) | (o, Singleton(n)) => { |
||||
if o.contains(n) { |
||||
Singleton(*n) |
||||
} else { |
||||
Empty |
||||
} |
||||
} |
||||
(Line(l1), Line(l2)) => l1.intersect(l2), |
||||
(Union(un1, un2), o) | (o, Union(un1, un2)) => { |
||||
Self::union(un1.intersect(o), un2.intersect(o)) |
||||
} |
||||
} |
||||
} |
||||
} |
@ -0,0 +1,684 @@
@@ -0,0 +1,684 @@
|
||||
use std::collections::{BTreeMap, BTreeSet}; |
||||
use std::fmt; |
||||
use std::iter::FromIterator; |
||||
|
||||
use crate::math::Scalar; |
||||
|
||||
// an unknown variable with an id
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)] |
||||
struct Unknown(i64); |
||||
|
||||
type UnknownSet = BTreeSet<Unknown>; |
||||
|
||||
trait Unknowns { |
||||
fn unknowns(&self) -> UnknownSet; |
||||
fn has_unknowns(&self) -> bool; |
||||
fn has_unknown(&self, u: Unknown) -> bool; |
||||
} |
||||
|
||||
impl Unknowns for Scalar { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
UnknownSet::new() |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
false |
||||
} |
||||
fn has_unknown(&self, _: Unknown) -> bool { |
||||
false |
||||
} |
||||
} |
||||
|
||||
impl Unknowns for Unknown { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
FromIterator::from_iter(Some(*self)) |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
true |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
*self == u |
||||
} |
||||
} |
||||
|
||||
impl fmt::Display for Unknown { |
||||
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { |
||||
write!(f, "u{}", self.0) |
||||
} |
||||
} |
||||
|
||||
#[derive(Clone, Debug, PartialEq)] |
||||
enum Expr { |
||||
Unkn(Unknown), |
||||
Const(Scalar), |
||||
Sum(Exprs), |
||||
Neg(Box<Expr>), |
||||
Product(Exprs), |
||||
Div(Box<Expr>, Box<Expr>), |
||||
} |
||||
|
||||
type Exprs = Vec<Expr>; |
||||
|
||||
impl Unknowns for Exprs { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
self.iter().flat_map(|e: &Expr| e.unknowns()).collect() |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
self.iter().any(|e: &Expr| e.has_unknowns()) |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
self.iter().any(|e: &Expr| e.has_unknown(u)) |
||||
} |
||||
} |
||||
|
||||
fn write_separated_exprs(es: &Exprs, f: &mut fmt::Formatter, sep: &str) -> fmt::Result { |
||||
let mut is_first = true; |
||||
for e in es { |
||||
if is_first { |
||||
is_first = false; |
||||
} else { |
||||
write!(f, "{}", sep)? |
||||
} |
||||
write!(f, "({})", e)? |
||||
} |
||||
Ok(()) |
||||
} |
||||
|
||||
fn remove_common_terms(l: &mut Vec<Expr>, r: &mut Vec<Expr>) -> Vec<Expr> { |
||||
let common: Vec<_> = l.drain_filter(|e| r.contains(e)).collect(); |
||||
common.iter().for_each(|e| { |
||||
r.remove_item(e); |
||||
}); |
||||
common |
||||
} |
||||
|
||||
fn remove_term(terms: &mut Vec<Expr>, term: &Expr) -> Option<Expr> { |
||||
terms.remove_item(term) |
||||
} |
||||
|
||||
fn sum_fold(l: Expr, r: Expr) -> Expr { |
||||
use itertools::Itertools; |
||||
use Expr::*; |
||||
match (l, r) { |
||||
(Const(lc), Const(rc)) => Const(lc + rc), |
||||
(Const(c), o) | (o, Const(c)) if relative_eq!(c, 0.) => o, |
||||
(Product(mut l), Product(mut r)) => { |
||||
let comm = remove_common_terms(&mut l, &mut r); |
||||
Expr::new_product(Sum(comm), Expr::new_sum(Product(l), Product(r))).simplify() |
||||
} |
||||
(Product(mut l), r) | (r, Product(mut l)) => { |
||||
let comm = remove_term(&mut l, &r); |
||||
match comm { |
||||
Some(_) => { |
||||
Expr::new_product(r, Expr::new_sum(Product(l), Const(1.))).simplify() |
||||
} |
||||
None => Expr::new_sum(Product(l), r), |
||||
} |
||||
} |
||||
(l, r) => Expr::new_sum(l, r), |
||||
} |
||||
} |
||||
|
||||
fn group_sum(es: Exprs) -> Exprs { |
||||
use Expr::*; |
||||
let mut common: BTreeMap<UnknownSet, Expr> = BTreeMap::new(); |
||||
for e in es { |
||||
let unkns = e.unknowns(); |
||||
match common.get_mut(&unkns) { |
||||
None => { |
||||
match e { |
||||
Const(c) if relative_eq!(c, 0.) => (), |
||||
e => { |
||||
common.insert(unkns, e); |
||||
} |
||||
}; |
||||
} |
||||
Some(existing) => { |
||||
match existing { |
||||
Sum(ref mut es) => { |
||||
// already failed at merging, so just add it to the list
|
||||
es.push(e); |
||||
} |
||||
other => { |
||||
*other = sum_fold(other.clone(), e); |
||||
} |
||||
}; |
||||
} |
||||
}; |
||||
} |
||||
common.into_iter().map(|(_, v)| v).collect() |
||||
} |
||||
|
||||
fn product_fold(l: Expr, r: Expr) -> Expr { |
||||
use itertools::Itertools; |
||||
use Expr::*; |
||||
match (l, r) { |
||||
(Const(lc), Const(rc)) => Const(lc * rc), |
||||
(Const(c), o) | (o, Const(c)) if relative_eq!(c, 1.) => o, |
||||
(Div(num, den), mul) | (mul, Div(num, den)) => { |
||||
if mul == *den { |
||||
*num |
||||
} else { |
||||
Expr::Div(Box::new(Expr::Product(vec![*num, mul])), den).simplify() |
||||
} |
||||
} |
||||
(l, r) => Expr::new_product(l, r), |
||||
} |
||||
} |
||||
|
||||
fn group_product(es: Exprs) -> Exprs { |
||||
use Expr::*; |
||||
let mut common: BTreeMap<UnknownSet, Expr> = BTreeMap::new(); |
||||
for e in es { |
||||
let unkns = e.unknowns(); |
||||
match common.get_mut(&unkns) { |
||||
None => { |
||||
match e { |
||||
Const(c) if relative_eq!(c, 1.) => (), |
||||
e => { |
||||
common.insert(unkns, e); |
||||
} |
||||
}; |
||||
} |
||||
Some(existing) => { |
||||
match existing { |
||||
Sum(ref mut es) => { |
||||
// already failed at merging, so just add it to the list
|
||||
es.push(e); |
||||
} |
||||
other => *other = product_fold(other.clone(), e), |
||||
}; |
||||
} |
||||
}; |
||||
} |
||||
common.into_iter().map(|(_, v)| v).collect() |
||||
} |
||||
|
||||
fn distribute_product_sums(mut es: Exprs) -> Expr { |
||||
trace!("distribute_product_sums: {}", Product(es.clone())); |
||||
use itertools::Itertools; |
||||
use Expr::*; |
||||
let sums = es |
||||
.drain_filter(|e| match e { |
||||
Sum(_) => true, |
||||
_ => false, |
||||
}) |
||||
.map(|e| { |
||||
trace!("sum in product: {}", e); |
||||
match e { |
||||
Sum(es) => es, |
||||
_ => unreachable!(), |
||||
} |
||||
}); |
||||
let products: Vec<_> = sums.multi_cartesian_product().collect(); |
||||
if products.is_empty() { |
||||
trace!("no sums to distribute"); |
||||
return Product(es); |
||||
} |
||||
let sums = products |
||||
.into_iter() |
||||
.map(|mut prod| { |
||||
prod.extend(es.clone()); |
||||
trace!("prod: {}", Product(prod.clone())); |
||||
Product(prod) |
||||
}) |
||||
.collect(); |
||||
Sum(sums) |
||||
} |
||||
|
||||
impl Unknowns for Expr { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
use Expr::*; |
||||
match self { |
||||
Unkn(u) => u.unknowns(), |
||||
Const(_) => UnknownSet::default(), |
||||
Sum(es) | Product(es) => es.unknowns(), |
||||
Div(l, r) => l.unknowns().union(&r.unknowns()).cloned().collect(), |
||||
Neg(e) => e.unknowns(), |
||||
} |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
use Expr::*; |
||||
match self { |
||||
Unkn(u) => u.has_unknowns(), |
||||
Const(_) => false, |
||||
Sum(es) | Product(es) => es.has_unknowns(), |
||||
Div(l, r) => l.has_unknowns() || r.has_unknowns(), |
||||
Neg(e) => e.has_unknowns(), |
||||
} |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
use Expr::*; |
||||
match self { |
||||
Unkn(u1) => u1.has_unknown(u), |
||||
Const(_) => false, |
||||
Sum(es) | Product(es) => es.has_unknown(u), |
||||
Div(l, r) => l.has_unknown(u) || r.has_unknown(u), |
||||
Neg(e) => e.has_unknown(u), |
||||
} |
||||
} |
||||
} |
||||
|
||||
impl Expr { |
||||
fn new_sum(e1: Expr, e2: Expr) -> Expr { |
||||
Expr::Sum(vec![e1, e2]) |
||||
} |
||||
fn new_product(e1: Expr, e2: Expr) -> Expr { |
||||
Expr::Product(vec![e1, e2]) |
||||
} |
||||
fn new_neg(e1: Expr) -> Expr { |
||||
Expr::Neg(Box::new(e1)) |
||||
} |
||||
fn new_div(num: Expr, den: Expr) -> Expr { |
||||
Expr::Div(Box::new(num), Box::new(den)) |
||||
} |
||||
fn new_minus(e1: Expr, e2: Expr) -> Expr { |
||||
Expr::Sum(vec![e1, Expr::new_neg(e2)]) |
||||
} |
||||
fn new_inv(den: Expr) -> Expr { |
||||
Expr::new_div(Expr::Const(1.), den) |
||||
} |
||||
|
||||
fn is_zero(&self) -> bool { |
||||
use Expr::*; |
||||
match self { |
||||
Const(c) => relative_eq!(*c, 0.), |
||||
_ => false, |
||||
} |
||||
} |
||||
|
||||
fn is_one(&self) -> bool { |
||||
use Expr::*; |
||||
match self { |
||||
Const(c) => relative_eq!(*c, 1.), |
||||
_ => false, |
||||
} |
||||
} |
||||
|
||||
fn simplify(self) -> Expr { |
||||
use Expr::*; |
||||
match self { |
||||
Sum(es) => { |
||||
let mut new_es: Vec<_> = es |
||||
.into_iter() |
||||
.map(|e| e.simplify()) |
||||
.flat_map(|e| match e { |
||||
Sum(more_es) => more_es, |
||||
other => vec![other], |
||||
}) |
||||
.collect(); |
||||
let pre_new_es = new_es.clone(); |
||||
new_es = group_sum(new_es); |
||||
trace!( |
||||
"simplify sum {} => {}", |
||||
Sum(pre_new_es), |
||||
Sum(new_es.clone()) |
||||
); |
||||
|
||||
match new_es.len() { |
||||
0 => Const(0.), // none
|
||||
1 => new_es.into_iter().next().unwrap(), // one
|
||||
_ => Sum(new_es), // many
|
||||
} |
||||
} |
||||
Product(es) => { |
||||
let new_es: Vec<_> = es |
||||
.into_iter() |
||||
.map(|e| e.simplify()) |
||||
.flat_map(|e| match e { |
||||
Product(more_es) => more_es, |
||||
other => vec![other], |
||||
}) |
||||
.collect(); |
||||
let pre_new_es = new_es.clone(); |
||||
let new_es = group_product(new_es); |
||||
trace!( |
||||
"simplify product {} => {}", |
||||
Product(pre_new_es), |
||||
Product(new_es.clone()) |
||||
); |
||||
match new_es.len() { |
||||
0 => Const(1.), // none
|
||||
1 => new_es.into_iter().next().unwrap(), // one
|
||||
_ => Product(new_es), // many
|
||||
} |
||||
} |
||||
Neg(mut v) => { |
||||
*v = v.simplify(); |
||||
trace!("simplify neg {}", Neg(v.clone())); |
||||
match v { |
||||
box Const(c) => Const(-c), |
||||
box Neg(v) => *v, |
||||
e => Product(vec![Const(-1.), *e]), |
||||
} |
||||
} |
||||
Div(mut num, mut den) => { |
||||
*num = num.simplify(); |
||||
*den = den.simplify(); |
||||
trace!("simplify div {}", Div(num.clone(), den.clone())); |
||||
match (num, den) { |
||||
(box Const(num), box Const(den)) => Const(num / den), |
||||
(num, box Const(den)) => { |
||||
if relative_eq!(den, 1.) { |
||||
*num |
||||
} else { |
||||
Expr::new_product(*num, Const(1. / den)) |
||||
} |
||||
} |
||||
(num, box Div(dennum, denden)) => { |
||||
Div(Box::new(Product(vec![*num, *denden])), dennum).simplify() |
||||
} |
||||
(box Product(mut es), box den) => match es.remove_item(&den) { |
||||
Some(_) => Product(es), |
||||
None => Expr::new_div(Product(es), den), |
||||
}, |
||||
(num, den) => { |
||||
if num == den { |
||||
Expr::Const(1.) |
||||
} else { |
||||
Div(num, den) |
||||
} |
||||
} |
||||
} |
||||
} |
||||
e => e, |
||||
} |
||||
} |
||||
|
||||
fn distribute(self) -> Expr { |
||||
use Expr::*; |
||||
trace!("distribute {}", self); |
||||
match self { |
||||
Sum(mut es) => { |
||||
for e in &mut es { |
||||
*e = e.clone().distribute(); |
||||
} |
||||
Sum(es) |
||||
} |
||||
Product(es) => distribute_product_sums(es), |
||||
Div(mut num, mut den) => { |
||||
*num = num.distribute(); |
||||
*den = den.distribute(); |
||||
match (num, den) { |
||||
(box Sum(es), box den) => Sum(es |
||||
.into_iter() |
||||
.map(|e| Expr::new_div(e, den.clone())) |
||||
.collect()), |
||||
(mut num, mut den) => Div(num, den), |
||||
} |
||||
} |
||||
Neg(v) => match v { |
||||
// box Sum(mut l, mut r) => {
|
||||
// *l = Neg(l.clone()).distribute();
|
||||
// *r = Neg(r.clone()).distribute();
|
||||
// Sum(l, r)
|
||||
// }
|
||||
// box Product(mut l, r) => {
|
||||
// *l = Neg(l.clone()).distribute();
|
||||
// Product(l, r)
|
||||
// }
|
||||
box Neg(v) => v.distribute(), |
||||
box Div(mut num, mut den) => { |
||||
*num = Neg(num.clone()).distribute(); |
||||
*den = Neg(den.clone()).distribute(); |
||||
Div(num, den) |
||||
} |
||||
e => Neg(e), |
||||
}, |
||||
e => e, |
||||
} |
||||
} |
||||
} |
||||
|
||||
impl fmt::Display for Expr { |
||||
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { |
||||
use Expr::*; |
||||
match self { |
||||
Unkn(u) => write!(f, "{}", u), |
||||
Const(c) => write!(f, "{}", c), |
||||
Sum(es) => write_separated_exprs(es, f, " + "), |
||||
Product(es) => write_separated_exprs(es, f, " * "), |
||||
Div(num, den) => write!(f, "({}) / ({})", num, den), |
||||
Neg(e) => write!(f, "-({})", e), |
||||
} |
||||
} |
||||
} |
||||
|
||||
#[derive(Clone, Debug, PartialEq)] |
||||
struct Eqn(Expr, Expr); |
||||
|
||||
impl Unknowns for Eqn { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
self.0 |
||||
.unknowns() |
||||
.union(&self.1.unknowns()) |
||||
.cloned() |
||||
.collect() |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
self.0.has_unknowns() || self.1.has_unknowns() |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
self.0.has_unknown(u) || self.1.has_unknown(u) |
||||
} |
||||
} |
||||
|
||||
fn ord_by_unkn(a: Expr, b: Expr, u: Unknown) -> Option<(Expr, Expr)> { |
||||
if a.has_unknown(u) { |
||||
Some((a, b)) |
||||
} else if b.has_unknown(u) { |
||||
Some((b, a)) |
||||
} else { |
||||
None |
||||
} |
||||
} |
||||
|
||||
impl Eqn { |
||||
fn simplify(self) -> Eqn { |
||||
Eqn(self.0.simplify(), self.1.simplify()) |
||||
} |
||||
|
||||
fn solve(&self, for_u: Unknown) -> Option<Expr> { |
||||
use Expr::*; |
||||
if !self.has_unknown(for_u) { |
||||
return None; |
||||
} |
||||
let (l, r) = ( |
||||
self.0 |
||||
.clone() /*.distribute()*/ |
||||
.simplify(), |
||||
self.1 |
||||
.clone() /*.distribute()*/ |
||||
.simplify(), |
||||
); |
||||
let (mut l, mut r) = ord_by_unkn(l, r, for_u)?; |
||||
loop { |
||||
trace!("solve: {} == {}", l, r); |
||||
let (new_l, new_r): (Expr, Expr) = match l { |
||||
Unkn(u) => return if u == for_u { Some(r.simplify()) } else { None }, |
||||
Sum(es) => { |
||||
let (us, not_us): (Vec<_>, Vec<_>) = |
||||
es.into_iter().partition(|e| e.has_unknown(for_u)); |
||||
if us.len() != 1 { |
||||
return None; |
||||
} |
||||
( |
||||
Sum(us).simplify(), |
||||
Expr::new_minus(r, Sum(not_us)).simplify(), |
||||
) |
||||
} |
||||
Product(es) => { |
||||
let (us, not_us): (Vec<_>, Vec<_>) = |
||||
es.into_iter().partition(|e| e.has_unknown(for_u)); |
||||
if us.len() != 1 { |
||||
return None; |
||||
} |
||||
( |
||||
Product(us).simplify(), |
||||
Expr::new_div(r, Product(not_us)).simplify(), |
||||
) |
||||
} |
||||
Neg(v) => (*v, Expr::new_neg(r)), |
||||
Div(num, den) => { |
||||
let (nu, du) = (num.has_unknown(for_u), den.has_unknown(for_u)); |
||||
match (nu, du) { |
||||
(true, false) => (*num, Expr::new_product(r, *den)), |
||||
(false, true) => (Expr::new_product(r, *den), *num), |
||||
(true, true) => return None, // TODO: simplify
|
||||
(false, false) => return None, |
||||
} |
||||
} |
||||
Const(_) => return None, |
||||
_ => return None, |
||||
}; |
||||
l = new_l; |
||||
r = new_r; |
||||
} |
||||
} |
||||
} |
||||
|
||||
impl fmt::Display for Eqn { |
||||
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { |
||||
write!(f, "{} == {}", self.0, self.1) |
||||
} |
||||
} |
||||
|
||||
#[derive(Clone, Debug, PartialEq)] |
||||
struct Eqns(Vec<Eqn>); |
||||
|
||||
impl Unknowns for Eqns { |
||||
fn unknowns(&self) -> UnknownSet { |
||||
self.0.iter().flat_map(|eqn: &Eqn| eqn.unknowns()).collect() |
||||
} |
||||
fn has_unknowns(&self) -> bool { |
||||
self.0.iter().any(|eqn: &Eqn| eqn.has_unknowns()) |
||||
} |
||||
fn has_unknown(&self, u: Unknown) -> bool { |
||||
self.0.iter().any(|eqn: &Eqn| eqn.has_unknown(u)) |
||||
} |
||||
} |
||||
|
||||
#[cfg(test)] |
||||
mod tests { |
||||
use super::*; |
||||
|
||||
#[test] |
||||
fn test_unknowns() { |
||||
let u1 = Unknown(1); |
||||
let u2 = Unknown(2); |
||||
let u3 = Unknown(3); |
||||
assert!(u1.has_unknowns()); |
||||
assert!(u2.has_unknowns()); |
||||
assert!(u1.has_unknown(u1)); |
||||
assert!(!u1.has_unknown(u2)); |
||||
assert!(u1.unknowns().contains(&u1)); |
||||
assert!(!u2.unknowns().contains(&u1)); |
||||
|
||||
let e1 = Expr::new_minus(Expr::Unkn(u1), Expr::Unkn(u2)); |
||||
assert!(e1.has_unknowns()); |
||||
assert!(e1.has_unknown(u1)); |
||||
assert!(e1.has_unknown(u2)); |
||||
assert!(!e1.has_unknown(u3)); |
||||
assert!(e1.unknowns().len() == 2); |
||||
} |
||||
|
||||
fn const_expr(e: Expr) -> Option<Scalar> { |
||||
match e { |
||||
Expr::Const(c) => Some(c), |
||||
_ => None, |
||||
} |
||||
} |
||||
|
||||
#[test] |
||||
fn test_solve() { |
||||
use Expr::*; |
||||
let _ = env_logger::try_init(); |
||||
let u1 = Unknown(1); |
||||
let e1 = Unkn(u1); |
||||
let e2 = Const(1.); |
||||
|
||||
let eqn = Eqn(e1.clone(), e2.clone()); |
||||
assert_eq!(eqn.solve(u1), Some(Const(1.))); |
||||
let eqn = Eqn(e2.clone(), e1.clone()); |
||||
assert_eq!(eqn.solve(u1), Some(Const(1.))); |
||||
|
||||
let e3 = Expr::new_sum(Const(1.), Expr::new_sum(Const(1.), Const(2.))); |
||||
let eqn = Eqn(e1.clone(), e3.clone()); |
||||
assert_eq!(eqn.solve(u1), Some(Const(4.))); |
||||
let e3 = Expr::new_minus(Const(1.), Const(1.)); |
||||
let eqn = Eqn(e1.clone(), e3.clone()); |
||||
assert_eq!(eqn.solve(u1), Some(Const(0.))); |
||||
|
||||
let e1 = Expr::new_div(Const(2.), Expr::new_minus(Const(1.), Const(4.))); |
||||
let e2 = Expr::new_minus(Const(1.), Unkn(u1)); |
||||
let eqn = Eqn(e1, e2); |
||||
info!("eqn: {} => {}", eqn, eqn.clone().simplify()); |
||||
let e = eqn.solve(u1).unwrap(); |
||||
assert!(const_expr(e.clone()).is_some()); |
||||
assert!(relative_eq!(const_expr(e.clone()).unwrap(), 5. / 3.)); |
||||
|
||||
let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); |
||||
let e2 = Expr::new_minus(Expr::new_product(Unkn(u1), Const(2.)), Unkn(u1)); |
||||
info!( |
||||
"e1==e2: {}=={} => {}=={}", |
||||
e1, |
||||
e2, |
||||
e1.clone().simplify(), |
||||
e2.clone().simplify() |
||||
); |
||||
let eqn = Eqn(e1, e2); |
||||
let e = eqn.solve(u1).unwrap(); |
||||
assert!(const_expr(e.clone()).is_some()); |
||||
assert!(relative_eq!(const_expr(e.clone()).unwrap(), -6.)); |
||||
|
||||
let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); |
||||
let e2 = Expr::new_div( |
||||
Expr::new_sum( |
||||
Expr::new_product(Unkn(u1), Const(2.)), |
||||
Expr::new_product(Unkn(u1), Unkn(u1)), |
||||
), |
||||
Unkn(u1), |
||||
); |
||||
info!( |
||||
"{}=={} distrib=> {}=={}", |
||||
e1, |
||||
e2, |
||||
e1.clone().distribute(), |
||||
e2.clone().distribute() |
||||
); |
||||
info!( |
||||
"{}=={} simplify=> {}=={}", |
||||
e1, |
||||
e2, |
||||
e1.clone().distribute().simplify(), |
||||
e2.clone().distribute().simplify() |
||||
); |
||||
let eqn = Eqn(e1, e2); |
||||
let e = eqn.solve(u1).unwrap(); |
||||
assert!(const_expr(e.clone()).is_some()); |
||||
assert!(relative_eq!(const_expr(e.clone()).unwrap(), -8.)); |
||||
|
||||
let e1 = Expr::new_product(Const(2.), Expr::new_minus(Const(1.), Const(4.))); |
||||
let e2 = Expr::new_div( |
||||
Expr::new_sum( |
||||
Expr::new_product(Unkn(u1), Const(2.)), |
||||
Expr::new_sum( |
||||
Expr::new_sum(Expr::new_product(Unkn(u1), Unkn(u1)), Unkn(u1)), |
||||
Expr::new_minus(Const(2.), Const(1. + 1.)), |
||||
), |
||||
), |
||||
Unkn(u1), |
||||
); |
||||
info!( |
||||
"e1==e2: {}=={} => {}=={}", |
||||
e1, |
||||
e2, |
||||
e1.clone().distribute().simplify(), |
||||
e2.clone().distribute().simplify().simplify() |
||||
); |
||||
let eqn = Eqn(e1, e2); |
||||
let e = eqn.solve(u1).unwrap(); |
||||
assert!(const_expr(e.clone()).is_some()); |
||||
assert!(relative_eq!(const_expr(e.clone()).unwrap(), -9.)); |
||||
} |
||||
|
||||
} |
Loading…
Reference in new issue