Day 23

This isn't fun anymore. Because I was tired of this one and didn't care
anymore, I looked at what others did and was baffled to see that there
were basically two groups of people: one that used Z3 or an ILP solver,
and one that wrote a direct solution that only worked on their own
input, by accident.

So I also used Z3. It worked. It takes a minute, though.

2 files changed, 495 insertions, 1 deletions

diff --git a/2018/src/day23.rs b/2018/src/day23.rs
new file mode 100644
index 0000000..6246ee9
--- /dev/null
+++ b/2018/src/day23.rs
@@ -0,0 +1,492 @@
+use std::io;
+use std::io::{BufRead, Write};
+use std::cmp;
+use std::process::{Command, Stdio};
+use rand::prelude::*;
+
+const LARGE: i64 = i64::max_value() / 16;
+
+#[allow(dead_code)]
+fn leq3<T: Ord>(a: T, b: T, c: T) -> bool {
+    a <= b && b <= c
+}
+
+#[derive(Clone, PartialEq, Eq, Debug)]
+struct Pos { x: i64, y: i64, z: i64 }
+
+impl Pos {
+    fn new(x: i64, y: i64, z: i64) -> Self {
+        Pos { x, y, z }
+    }
+
+    fn new_v(v: i64) -> Self {
+        Pos::new(v, v, v)
+    }
+
+    #[used]
+    fn dist(&self, b: &Pos) -> i64 {
+        (self.x - b.x).abs() + (self.y - b.y).abs() + (self.z - b.z).abs()
+    }
+
+    fn distf(&self, b: &FPos) -> f64 {
+        (self.x as f64 - b.x).abs() + (self.y as f64 - b.y).abs() + (self.z as f64 - b.z).abs()
+    }
+
+    #[used]
+    fn dot(&self, b: &Pos) -> i64 {
+        self.x * b.x + self.y * b.y + self.z * b.z
+    }
+
+    #[used]
+    fn add(&self, b: &Pos) -> Pos {
+        Pos::new(self.x + b.x, self.y + b.y, self.z + b.z)
+    }
+
+    fn max(&self, b: &Pos) -> Pos {
+        Pos::new(cmp::max(self.x, b.x), cmp::max(self.y, b.y), cmp::max(self.z, b.z))
+    }
+
+    fn min(&self, b: &Pos) -> Pos {
+        Pos::new(cmp::min(self.x, b.x), cmp::min(self.y, b.y), cmp::min(self.z, b.z))
+    }
+}
+
+#[derive(Debug)]
+struct Bot {
+    pos: Pos,
+    rad: i64,
+}
+
+fn parse_bot(line: &str) -> Bot {
+    let mut numbers = vec![0];
+    let mut multiplier = 1;
+
+    for c in line.chars() {
+        if c == ',' {
+            numbers.push(0);
+            multiplier = 1;
+        } else if c == '-' {
+            multiplier = -1;
+        } else if let Some(dig) = c.to_digit(10) {
+            let last = numbers.last_mut().unwrap();
+            *last = *last * 10 + multiplier * dig as i64;
+        }
+    }
+
+    assert!(numbers.len() == 4);
+    Bot {
+        pos: Pos { x: numbers[0], y: numbers[1], z: numbers[2] },
+        rad: numbers[3],
+    }
+}
+
+#[derive(Clone, PartialEq, Eq, Debug)]
+struct Intersection {
+    // If p is any point inside the bounds, then order the 8 planes by intersection with the rays:
+    // <(1, 1, 1)>, <(1, 1, -1)>, <(1, -1, 1)>, <(1, -1, -1)>, <(-1, 1, 1)>, <(-1, 1, -1)>, <(-1, -1, 1)>, <(-1, -1, -1)>
+    // Then represent a plane by the x-coordinate of its intersection with the x-axis.
+    // All bounds are inclusive.
+    bounds: [i64; 8],
+}
+
+impl Intersection {
+    fn all() -> Self {
+        let pos = LARGE;
+        let neg = -LARGE;
+        Self { bounds: [pos, pos, pos, pos, neg, neg, neg, neg] }
+    }
+
+    #[used]
+    fn empty() -> Self {
+        let pos = LARGE;
+        let neg = -LARGE;
+        Self { bounds: [neg, neg, neg, neg, pos, pos, pos, pos] }
+    }
+
+    fn is_empty(&self) -> bool {
+        self.bounds[0] < self.bounds[7]
+            || self.bounds[1] < self.bounds[6]
+            || self.bounds[2] < self.bounds[5]
+            || self.bounds[3] < self.bounds[4]
+    }
+
+    fn from_bot(bot: &Bot) -> Self {
+        let pos = bot.pos.add(&Pos::new(bot.rad, 0, 0));
+        let neg = bot.pos.add(&Pos::new(-bot.rad, 0, 0));
+
+        let res = Self {
+            bounds: [
+                pos.dot(&Pos::new(1, 1, 1)),
+                pos.dot(&Pos::new(1, 1, -1)),
+                pos.dot(&Pos::new(1, -1, 1)),
+                pos.dot(&Pos::new(1, -1, -1)),
+                neg.dot(&Pos::new(1, -1, -1)),
+                neg.dot(&Pos::new(1, -1, 1)),
+                neg.dot(&Pos::new(1, 1, -1)),
+                neg.dot(&Pos::new(1, 1, 1)),
+            ]
+        };
+
+        if res.bounds[0] - res.bounds[7] != 2 * bot.rad
+                || res.bounds[1] - res.bounds[6] != 2 * bot.rad
+                || res.bounds[2] - res.bounds[5] != 2 * bot.rad
+                || res.bounds[3] - res.bounds[4] != 2 * bot.rad {
+            println!("bot={:?}", bot);
+            println!("bounds={:?}", res.bounds);
+            assert!(false);
+        }
+
+        res
+    }
+
+    fn intersection(&self, o: &Self) -> Self {
+        Self {
+            bounds: [
+                cmp::min(self.bounds[0], o.bounds[0]),
+                cmp::min(self.bounds[1], o.bounds[1]),
+                cmp::min(self.bounds[2], o.bounds[2]),
+                cmp::min(self.bounds[3], o.bounds[3]),
+                cmp::max(self.bounds[4], o.bounds[4]),
+                cmp::max(self.bounds[5], o.bounds[5]),
+                cmp::max(self.bounds[6], o.bounds[6]),
+                cmp::max(self.bounds[7], o.bounds[7]),
+            ]
+        }
+    }
+
+    fn dist_origin(&self) -> i64 {
+        println!("{:?}", self);
+        let pts = self.all_points();
+        println!("{} points in self", pts.len());
+        let res = pts.iter().map(|p| p.dist(&Pos::new(0, 0, 0))).min().unwrap();
+        println!("{}", res);
+        res
+    }
+
+    fn contains(&self, pos: &Pos) -> bool {
+        leq3(self.bounds[7], pos.dot(&Pos::new(1, 1, 1)), self.bounds[0])
+            && leq3(self.bounds[6], pos.dot(&Pos::new(1, 1, -1)), self.bounds[1])
+            && leq3(self.bounds[5], pos.dot(&Pos::new(1, -1, 1)), self.bounds[2])
+            && leq3(self.bounds[4], pos.dot(&Pos::new(1, -1, -1)), self.bounds[3])
+    }
+
+    fn all_points(&self) -> Vec<Pos> {
+        (self.bounds[5] + self.bounds[6] ..= self.bounds[1] + self.bounds[2]).map(move |p1| {
+            (self.bounds[6] + self.bounds[7] - p1 ..= self.bounds[0] + self.bounds[1] - p1).map(move |p2| {
+                (-self.bounds[3] + p1 - p2 ..= -self.bounds[4] + p1 - p2).map(move |p3| {
+                    println!("candidate {:?}", Pos::new(p1, p2, p3));
+                    Pos::new(p1, p2, p3)
+                })
+            }).flatten()
+        }).flatten().filter(|p| self.contains(p)).collect()
+    }
+}
+
+#[allow(dead_code)]
+fn backtrack<F: Fn(&Intersection) -> i64>(bots: &[Intersection], selected: &[usize], max_count: &mut i64, min_dist: &mut i64, count: i64, inter: Intersection, dist_func: &F) {
+    // println!("backtrack: max_count={} min_dist={} count={}", max_count, min_dist, count);
+
+    if count > *max_count {
+        *max_count = count;
+        *min_dist = dist_func(&inter);
+        // println!("New max count: {} at dist {}; inter={:?}", count, *min_dist, inter);
+    } else if count == *max_count {
+        let d = dist_func(&inter);
+        if d < *min_dist {
+            *min_dist = d;
+            // println!("  New min dist: {}", d);
+        }
+    }
+
+    let rest = selected.iter().map(|&i| i).filter(|&i| !inter.intersection(&bots[i]).is_empty()).collect::<Vec<_>>();
+    // println!("rest.len() = {}", rest.len());
+
+    if rest.len() as i64 + count < *max_count {
+        return;
+    }
+
+    if rest.len() == 0 {
+        return;
+    }
+
+    let inter2 = inter.intersection(&bots[rest[0]]);
+    backtrack(bots, &rest[1..], max_count, min_dist, count + 1, inter2, dist_func);
+    backtrack(bots, &rest[1..], max_count, min_dist, count, inter, dist_func);
+}
+
+#[allow(dead_code)]
+fn part2_via_backtracking(bots: &[Bot]) -> i64 {
+    let inters = bots.iter().map(|b| Intersection::from_bot(&b)).collect::<Vec<_>>();
+    let mut max_count = 0;
+    let mut min_dist = LARGE;
+    backtrack(&inters, &(0..inters.len()).collect::<Vec<_>>(), &mut max_count, &mut min_dist, 0, Intersection::all(), &|_| 0);
+    min_dist = LARGE;
+    backtrack(&inters, &(0..inters.len()).collect::<Vec<_>>(), &mut max_count, &mut min_dist, 0, Intersection::all(), &|i| i.dist_origin());
+    min_dist
+}
+
+fn partial_min<'a, T: PartialOrd>(a: &'a T, b: &'a T) -> Option<&'a T> {
+    match a.partial_cmp(b) {
+        Some(cmp::Ordering::Less) => Some(a),
+        Some(cmp::Ordering::Equal) => Some(a),
+        Some(cmp::Ordering::Greater) => Some(b),
+        None => None
+    }
+}
+
+fn partial_max<'a, T: PartialOrd>(a: &'a T, b: &'a T) -> Option<&'a T> {
+    match a.partial_cmp(b) {
+        Some(cmp::Ordering::Less) => Some(b),
+        Some(cmp::Ordering::Equal) => Some(a),
+        Some(cmp::Ordering::Greater) => Some(a),
+        None => None
+    }
+}
+
+#[derive(Clone, PartialEq, Debug)]
+struct FPos { x: f64, y: f64, z: f64 }
+
+impl FPos {
+    fn new(x: f64, y: f64, z: f64) -> Self {
+        Self { x, y, z }
+    }
+
+    fn from_pos(pos: &Pos) -> Self {
+        Self { x: pos.x as f64, y: pos.y as f64, z: pos.z as f64 }
+    }
+
+    fn add(&self, b: &Self) -> Self {
+        Self::new(self.x + b.x, self.y + b.y, self.z + b.z)
+    }
+
+    fn scale(&self, b: f64) -> Self {
+        Self::new(self.x * b, self.y * b, self.z * b)
+    }
+
+    fn max(&self, b: &Self) -> Self {
+        Self::new(*partial_max(&self.x, &b.x).unwrap(), *partial_max(&self.y, &b.y).unwrap(), *partial_max(&self.z, &b.z).unwrap())
+    }
+
+    fn min(&self, b: &Self) -> Self {
+        Self::new(*partial_min(&self.x, &b.x).unwrap(), *partial_min(&self.y, &b.y).unwrap(), *partial_min(&self.z, &b.z).unwrap())
+    }
+}
+
+fn dist_weight(bot: &Bot, pos: &FPos) -> f64 {
+    let d = bot.pos.distf(pos);
+    if d <= bot.rad as f64 {
+        1.0
+    } else {
+        let x = d - bot.rad as f64;
+        0.0009 / (x + 1.0)
+    }
+}
+
+fn score_function(bots: &[Bot], pos: &FPos) -> f64 {
+    bots.iter().map(|bot| dist_weight(bot, pos)).sum()
+}
+
+fn random_point_in_box(minpt: &FPos, maxpt: &FPos) -> FPos {
+    FPos::new(random::<f64>() * (maxpt.x - minpt.x) + minpt.x,
+              random::<f64>() * (maxpt.y - minpt.y) + minpt.y,
+              random::<f64>() * (maxpt.z - minpt.z) + minpt.z)
+}
+
+fn random_point() -> FPos {
+    random_point_in_box(&FPos::new(-1.0, -1.0, -1.0), &FPos::new(1.0, 1.0, 1.0))
+}
+
+fn improve(bots: &[Bot], pos: FPos, heat: f64) -> (FPos, f64) {
+    let mut best_pos = pos.clone();
+    let mut best_score = score_function(bots, &pos);
+
+    for _ in 0..10 {
+        let pt = pos.add(&random_point().scale(heat));
+        let sc = score_function(bots, &pt);
+        if sc > best_score {
+            best_score = sc;
+            best_pos = pt;
+        }
+    }
+
+    (best_pos, best_score)
+}
+
+fn containing_box(pts: &[FPos]) -> (FPos, FPos) {
+    let minpt = pts.iter().fold(FPos::new(1.0e34, 1.0e34, 1.0e34), |p1, p| p1.min(p));
+    let maxpt = pts.iter().fold(FPos::new(-1.0e34, -1.0e34, -1.0e34), |p1, p| p1.max(p));
+    (minpt, maxpt)
+}
+
+#[allow(dead_code)]
+fn anneal_maximum(bots: &[Bot]) -> FPos {
+    let minpt = FPos::from_pos(&bots.iter().fold(Pos::new_v(LARGE), |p1, bot| p1.min(&bot.pos)));
+    let maxpt = FPos::from_pos(&bots.iter().fold(Pos::new_v(-LARGE), |p1, bot| p1.max(&bot.pos)));
+
+    let num_points = 100;
+
+    let mut pts = Vec::new();
+    for _ in 0..num_points {
+        pts.push(random_point_in_box(&minpt, &maxpt))
+    }
+
+    let mut heat = 1000000.0;
+    for iter in 0..16000 {
+        let mut best_score = 0.0;
+
+        for i in 0..pts.len() {
+            match improve(bots, pts[i].clone(), heat) {
+                (p, s) => {
+                    pts[i] = p;
+                    if s > best_score {
+                        best_score = s;
+                    }
+                }
+            }
+        }
+
+        if iter % 100 == 0 {
+            println!("iter={} heat={} best_score={}  {:?}", iter, heat, best_score, containing_box(&pts));
+            heat *= 0.9;
+            if heat < 1.0 { heat = 1.0; }
+
+            pts.sort_by(|a, b| score_function(bots, b).partial_cmp(&score_function(bots, a)).unwrap());
+            pts.truncate(num_points / 2);
+
+            for _ in 0..num_points / 2 {
+                pts.push(random_point_in_box(&minpt, &maxpt));
+            }
+        }
+    }
+
+    pts.iter().max_by_key(|pt| score_function(bots, &pt) as i64).unwrap().clone()
+}
+
+fn build_smtlib(bots: &[Bot]) -> String {
+    let mut res = String::new();
+    res.push_str("(define-fun abs ((x Int)) Int (ite (< x 0) (- x) x))\n");
+
+    res.push_str("(declare-const x Int)\n");
+    res.push_str("(declare-const y Int)\n");
+    res.push_str("(declare-const z Int)\n");
+    res.push_str("(declare-const nrange Int)\n");
+    res.push_str("(declare-const distorigin Int)\n");
+
+    for i in 0..bots.len() {
+        res.push_str(&format!("(declare-const eq{} Int)\n", i));
+    }
+
+    for i in 0..bots.len() {
+        res.push_str(&format!("(assert (= eq{} (ite (<= (+ (abs (- x {})) (abs (- y {})) (abs (- z {}))) {}) 1 0)))\n",
+                              i, bots[i].pos.x, bots[i].pos.y, bots[i].pos.z, bots[i].rad));
+    }
+
+    res.push_str("(assert (= nrange (+\n");
+    for i in 0..bots.len() {
+        res.push_str(&format!("  (ite (<= (+ (abs (- x {})) (abs (- y {})) (abs (- z {}))) {}) 1 0)\n",
+                              bots[i].pos.x, bots[i].pos.y, bots[i].pos.z, bots[i].rad));
+    }
+    res.push_str(")))\n");
+
+    res.push_str("(assert (= distorigin (+ (abs x) (abs y) (abs z))))\n");
+    res.push_str("(maximize nrange)\n");
+    res.push_str("(minimize distorigin)\n");
+    res.push_str("(check-sat)\n");
+    // res.push_str("(get-model)\n");
+    res.push_str("(eval distorigin)\n");
+
+    res
+}
+
+#[allow(dead_code)]
+fn build_cplex(bots: &[Bot]) -> String {
+    let minpt = &bots.iter().fold(Pos::new_v(LARGE), |p1, bot| p1.min(&bot.pos));
+
+    let isign = |c, case| if case & (1 << c) == 0 { "" } else { "-" };
+    let sign = |c, case| if case & (1 << c) == 0 { "+" } else { "-" };
+    let mul = |c, case, n: i64| if case & (1 << c) == 0 { n } else { -n };
+
+    let mut res = String::new();
+    res.push_str("maximize 300000000 nrange - distorigin\n");
+
+    res.push_str("subject to\n");
+    for i in 0..bots.len() {
+        for case in 0..8 {
+            // bC = bot.C - min.C
+            // ±(x - bx) + ±(y - by) + ±(z - bz) <= rad + 1e12 (1 - ei)
+            // ±(x + y + z) - ±(bx + by + bz) <= rad + 1e12 - 1e12 ei
+            // ±(x + y + z) - ±(bx + by + bz) + 1e12 ei <= rad + 1e12
+            // ±x ± y ± z + 1e12 ei <= rad + 1e12 ± bx ± by ± bz
+            res.push_str(&format!(" {}x {} y {} z + 1000000000000 e{} <= {}\n",
+                                  isign(0, case), sign(1, case), sign(2, case), i,
+                                  bots[i].rad + 1000000000000 + mul(0, case, bots[i].pos.x - minpt.x) + mul(1, case, bots[i].pos.x - minpt.x) + mul(2, case, bots[i].pos.x - minpt.x)));
+        }
+    }
+
+    res.push_str(" nrange");
+    for i in 0..bots.len() {
+        res.push_str(&format!(" - e{}", i));
+    }
+    res.push_str(" = 0\n");
+
+    for case in 0..8 {
+        // ±(x - mx) + ±(y - my) + ±(z - mz) <= distorigin
+        // ±(x + y + z) - ±(mx + my + mz) <= distorigin
+        // ±x ± y ± z <= distorigin - ±(mx + my + mz)
+        // ±x ± y ± z - distorigin <= - ±mx - ±my - ±mz
+        res.push_str(&format!(" {}x {} y {} z - distorigin <= {}\n",
+                              isign(0, case), sign(1, case), sign(2, case),
+                              -mul(0, case, minpt.x) - mul(1, case, minpt.y) - mul(2, case, minpt.z)));
+    }
+
+    res.push_str("bounds\n");
+    res.push_str(" 0 <= nrange\n");
+    res.push_str(" 0 <= distorigin\n");
+    res.push_str(" 0 <= x\n");
+    res.push_str(" 0 <= y\n");
+    res.push_str(" 0 <= z\n");
+
+    res.push_str("general\n");
+    res.push_str(" nrange distorigin x y z\n");
+
+    res.push_str("binary\n");
+    for i in 0..bots.len() {
+        res.push_str(&format!(" e{}", i));
+    }
+    res.push_str("\n");
+
+    res.push_str("end\n");
+
+    res
+}
+
+pub fn main<T: BufRead>(reader: T) -> io::Result<(String, String)> {
+    let bots = reader.lines().map(|l| parse_bot(&l.unwrap())).collect::<Vec<_>>();
+
+    // println!("{}", build_smtlib(&bots));
+    // println!("{}", build_cplex(&bots));
+
+    let largest = (0..bots.len()).max_by_key(|&i| bots[i].rad).unwrap();
+
+    let part1 = bots.iter().filter(|bot| bot.pos.dist(&bots[largest].pos) <= bots[largest].rad).count();
+
+    // let part2 = anneal_maximum(&bots);
+    // let part2 = part2_via_backtracking(&bots);
+
+    let mut child = Command::new("z3")
+        .arg("-in")
+        .stdin(Stdio::piped())
+        .stdout(Stdio::piped())
+        .spawn()
+        .unwrap();
+
+    child.stdin.as_mut().unwrap().write_all(build_smtlib(&bots).as_bytes()).unwrap();
+
+    let output = String::from_utf8(child.wait_with_output().unwrap().stdout).unwrap();
+
+    let part2 = output.lines().last().unwrap().to_string();
+
+    // Ok((part1.to_string(), String::new()))
+
+    Ok((part1.to_string(), part2))
+}
diff --git a/2018/src/main.rs b/2018/src/main.rs
index c45b084..aac9302 100644
--- a/2018/src/main.rs
+++ b/2018/src/main.rs
@@ -28,8 +28,9 @@ mod day19;
 mod day20;
 mod day21;
 mod day22;
+mod day23;
 
-static NUM_DAYS: i32 = 22;
+static NUM_DAYS: i32 = 23;
 
 fn day_switch<T: BufRead>(day: i32, reader: T) -> io::Result<(String, String)> {
     match day {
@@ -55,6 +56,7 @@ fn day_switch<T: BufRead>(day: i32, reader: T) -> io::Result<(String, String)> {
         20 => day20::main(reader),
         21 => day21::main(reader),
         22 => day22::main(reader),
+        23 => day23::main(reader),
         _ => Err(Error::new(ErrorKind::Other, "Invalid day"))
     }
 }