Day 19: Beacon Scanner

use std::env;
use std::io::{self};
use std::collections::{HashSet,HashMap};

extern crate regex;
use regex::Regex;

use point3d::point3d::Point3D;

type Tuple3D = ((i64,i64,i64),(i64,i64,i64),(i64,i64,i64));

// https://www.euclideanspace.com/maths/algebra/matrix/transforms/examples/index.htm
const ROT01: Tuple3D =
    ((  1,  0,  0),
     (  0,  1,  0),
     (  0,  0,  1));

const ROT02: Tuple3D =
    ((  1,  0,  0),
     (  0,  0, -1),
     (  0,  1,  0));

const ROT03: Tuple3D =
    ((  1,  0,  0),
     (  0, -1,  0),
     (  0,  0, -1));

const ROT04: Tuple3D =
    ((  1,  0,  0),
     (  0,  0,  1),
     (  0, -1,  0));

const ROT05: Tuple3D =
    ((  0, -1,  0),
     (  1,  0,  0),
     (  0,  0,  1));

const ROT06: Tuple3D =
    ((  0,  0,  1),
     (  1,  0,  0),
     (  0,  1,  0));

const ROT07: Tuple3D =
    ((  0,  1,  0),
     (  1,  0,  0),
     (  0,  0, -1));

const ROT08: Tuple3D =
    ((  0,  0, -1),
     (  1,  0,  0),
     (  0, -1,  0));

const ROT09: Tuple3D =
    (( -1,  0,  0),
     (  0, -1,  0),
     (  0,  0,  1));

const ROT10: Tuple3D =
    (( -1,  0,  0),
     (  0,  0, -1),
     (  0, -1,  0));

const ROT11: Tuple3D =
    (( -1,  0,  0),
     (  0,  1,  0),
     (  0,  0, -1));

const ROT12: Tuple3D =
    (( -1,  0,  0),
     (  0,  0,  1),
     (  0,  1,  0));

const ROT13: Tuple3D =
    ((  0,  1,  0),
     ( -1,  0,  0),
     ( 0,   0,  1));

const ROT14: Tuple3D =
    ((  0,  0,  1),
     ( -1,  0,  0),
     (  0, -1,  0));

const ROT15: Tuple3D =
    ((  0, -1,  0),
     ( -1,  0,  0),
     (  0,  0, -1));

const ROT16: Tuple3D =
    ((  0,  0, -1),
     ( -1,  0,  0),
     (  0,  1,  0));

const ROT17: Tuple3D =
    ((  0,  0, -1),
     (  0,  1,  0),
     (  1,  0,  0));

const ROT18: Tuple3D =
    ((  0,  1,  0),
     (  0,  0,  1),
     (  1,  0,  0));

const ROT19: Tuple3D =
    ((  0,  0,  1),
     (  0, -1,  0),
     (  1,  0,  0));

const ROT20: Tuple3D =
    ((  0, -1,  0),
     (  0,  0, -1),
     (  1,  0,  0));

const ROT21: Tuple3D =
    ((  0,  0, -1),
     (  0, -1,  0),
     ( -1,  0,  0));

const ROT22: Tuple3D =
    ((  0, -1,  0),
     (  0,  0,  1),
     ( -1,  0,  0)); 

const ROT23: Tuple3D =
    ((  0,  0,  1),
     (  0,  1,  0),
     ( -1,  0,  0));

const ROT24: Tuple3D =
    ((  0,  1,  0),
     (  0,  0, -1),
     ( -1,  0,  0));

const ALL_ROTATIONS: [Tuple3D; 24] = 
        [ROT01, ROT02, ROT03, ROT04, ROT05, ROT06,
         ROT07, ROT08, ROT09, ROT10, ROT11, ROT12,
         ROT13, ROT14, ROT15, ROT16, ROT17, ROT18,
         ROT19, ROT20, ROT21, ROT22, ROT23, ROT24];
// https://www.euclideanspace.com/maths/algebra/matrix/transforms/examples/index.htm

#[derive(Clone)]
struct Scanner {
    id: u8,
    beacons: Vec<Point3D>,
    fixed: bool,
    origin: Point3D,
}
impl From<&str> for Scanner {
    fn from(s: &str) -> Self {
        let s: Vec<_> = s.split("\n").collect();

        // Get id
        let re = Regex::new(r"scanner (\d+)").unwrap();
        let matches = re.captures(s[0]).unwrap();
        let id = matches[1].parse().unwrap();

        // Get beacons
        let mut beacons: Vec<Point3D> = Vec::new();
        let re = Regex::new(r"(\-*\d+),(\-*\d+),(\-*\d+)").unwrap();
        for line in s[1..].into_iter() {
            let matches = re.captures(line).unwrap();
            let (x,y,z) = (matches[1].parse().unwrap(),matches[2].parse().unwrap(),matches[3].parse().unwrap());
            beacons.push(Point3D { x: x, y: y, z: z });
        }
        Self { id: id, beacons: beacons, fixed: false, origin: Point3D::zeros() }
    }
}

fn rotate_by(pt: &mut Point3D, rotation: Tuple3D) {
    let pt0 = pt.clone();
    pt.x = rotation.0.0*pt0.x + rotation.0.1*pt0.y+rotation.0.2*pt0.z;
    pt.y = rotation.1.0*pt0.x + rotation.1.1*pt0.y+rotation.1.2*pt0.z;
    pt.z = rotation.2.0*pt0.x + rotation.2.1*pt0.y+rotation.2.2*pt0.z;
}
fn rotate_all_by(pts: &mut Vec<Point3D>, rotation: Tuple3D) {
    pts.iter_mut().for_each(|mut x| rotate_by(&mut x,rotation))
}

fn translate_xyz(pt: &mut Point3D, delx: i64, dely: i64, delz: i64) {
    pt.x += delx;
    pt.y += dely;
    pt.z += delz;
}
fn translate_all_xyz(pts: &mut Vec<Point3D>, delx: i64, dely: i64, delz: i64) {
    pts.iter_mut().for_each(|mut x| translate_xyz(&mut x,delx,dely,delz))
}

fn manhattan_distance(pt1: &Point3D, pt2: &Point3D) -> i64 {
    (pt1.x-pt2.x).abs() + (pt1.y-pt2.y).abs() + (pt1.z-pt2.z).abs()
}

fn solve(input: &str) -> io::Result<()> {
    // Input
    let input_str = std::fs::read_to_string(input).unwrap();
    let input_str = input_str.trim();
    let input: Vec<_> = input_str.split("\n\n").collect();
    
    // Get scanner and beacon data
    let mut scanners: Vec<_> = input
        .iter()
        .map(|x| Scanner::from(*x))
        .collect();
    
    // Assume orientation of first scanner
    scanners[0].fixed = true;

    // Generate initial composites from scanner 0
    let mut composite: HashSet<Point3D> = HashSet::new();
    for b in &scanners[0].beacons {
        composite.insert(*b);
    }

    'main: loop {
        // Try to match unfixed scanners against fixed scanners
        'scanner: for scanner in &mut scanners {
            if scanner.fixed { continue; }

            for rot in ALL_ROTATIONS {
                let mut beacons = scanner.beacons.clone();
                rotate_all_by(&mut beacons, rot);

                // Count x-, y-, z-offsets between these and known points
                let mut xcounts: HashMap<i64,i64> = HashMap::new();
                let mut ycounts: HashMap<i64,i64> = HashMap::new();
                let mut zcounts: HashMap<i64,i64> = HashMap::new();
                for beacon in &beacons {
                    for comp in composite.iter() {
                        *xcounts.entry(comp.x - beacon.x).or_insert(0) += 1;
                        *ycounts.entry(comp.y - beacon.y).or_insert(0) += 1;
                        *zcounts.entry(comp.z - beacon.z).or_insert(0) += 1;
                    }
                }

                // Generate x-, y-, z-offsets from counts and align scanner (if applicable)
                let xoffset = xcounts.iter().find(|&(_,v)| *v >= 12).map(|(k,_)| k);
                let yoffset = ycounts.iter().find(|&(_,v)| *v >= 12).map(|(k,_)| k);
                let zoffset = zcounts.iter().find(|&(_,v)| *v >= 12).map(|(k,_)| k);
                
                // Sort offset results to avoid non-deterministic solution
                let mut xcounts2: Vec<_> = xcounts.iter().collect(); xcounts2.sort_by(|a,b| b.1.cmp(a.1));
                let mut ycounts2: Vec<_> = ycounts.iter().collect(); ycounts2.sort_by(|a,b| b.1.cmp(a.1));
                let mut zcounts2: Vec<_> = zcounts.iter().collect(); zcounts2.sort_by(|a,b| b.1.cmp(a.1));
                if xoffset.is_some() && yoffset.is_some() && zoffset.is_some() {
                    
                    let (newx,_) = *xcounts2.iter().next().unwrap();
                    let (newy,_) = *ycounts2.iter().next().unwrap();
                    let (newz,_) = *zcounts2.iter().next().unwrap();

                    scanner.beacons = beacons.clone();
                    scanner.origin = Point3D::new(*newx,*newy,*newz);
                    scanner.fixed = true;

                    // Add to composite
                    translate_all_xyz(&mut beacons,*newx,*newy,*newz);
                    for beacon in &beacons {
                        composite.insert(*beacon);
                    }
                    continue 'scanner;
                }
            }
        }

        // Check all scanners have been fixed
        if scanners.iter().all(|s| s.fixed) { break 'main; }
    }
    println!("Part 1: {}", composite.len()); // 318

    // Part 2
    let mut part2 = 0;
    for scanner1 in &scanners {
        for scanner2 in &scanners {
            part2 = std::cmp::max(part2,
                manhattan_distance(&scanner1.origin,&scanner2.origin));
        }
    }
    println!("Part 2: {}",part2); // 12166


    Ok(())
}

fn main() {
    let args: Vec<String> = env::args().collect();
    let filename = &args[1];
    solve(&filename).unwrap();
}

#[macro_use]
extern crate text_io;

pub mod point3d {
    use std::ops::{Add,AddAssign};

    #[derive(Debug, Copy, Clone, Hash)]
    pub struct Point3D {
        pub x: i64,
        pub y: i64,
        pub z: i64,
    }

    impl Point3D {
        pub fn new(x: i64, y: i64, z: i64) -> Self {
            Self {
                x: x, y: y, z: z,
            }
        }
        pub fn zeros() -> Self {
            Self {
                x: 0, y: 0, z: 0,
            }
        }
        pub fn as_tuple(self) -> (i64,i64,i64) {
            (self.x, self.y, self.z)
        }
    }
    impl Add for Point3D {
        type Output = Self;
        fn add(self, other: Self) -> Self {
            Self {
                x: self.x + other.x,
                y: self.y + other.y,
                z: self.z + other.z,
            }
        }
    }
    impl AddAssign for Point3D {
        fn add_assign(&mut self, other: Self) {
            *self = Self {
                x: self.x + other.x,
                y: self.y + other.y,
                z: self.z + other.z,
            };
        }
    }
    impl PartialEq for Point3D {
        fn eq(&self, other: &Self) -> bool {
            self.x == other.x && self.y == other.y && self.z == other.z
        }
    }
    impl Eq for Point3D {}

}

This code is all still a mess, with leavings from various dead-ends still hanging around. I may come back around later and try to clean it up and describe it better, but for now I'll just dump what I've got and try to give a general overview.

Edit: Oh, I should add; the one big thing that got me going again was basically to give up on trying to figure 3D/matrix math and rotations properly, and just look at the ways that a given star at [x,y,z] might be transformed. The axes might get shuffled around (ie [y,z,x]), and the signs might get flipped eg ([-x,-y,-z]). Instead of working out proper rotations and trying those, I just generate every possible combination of axis-shuffling (which I call rotations in my code) and sign-flipping (inversion). This means more things to check, but I know the right combination has to be in there somewhere, and since I'm only looking at three stars at a time, it doesn't add too many options.

def dist3d(a, b)
  Math.sqrt((b[0] - a[0])**2 + (b[1] - a[1])**2 + (b[2] - a[2])**2)
end

def beacon_triangles(beacons)
  beacons.combination(3).map do |a, b, c|
    ab = dist3d(a, b)
    bc = dist3d(b, c)
    ca = dist3d(c, a)
    [[ab, bc, ca].sort, Set.new([a, b, c])]
  end.to_h
end

def beacon_transformations(beacon)
  transforms = [] # list of [beacon, ['rotation', inversion]] pairs

  rotations  = [0,1,2].permutation(3).to_a
  inversions = [[1,1,1], [1,1,-1], [1,-1,1], [-1,1,1], [1,-1,-1], [-1,1,-1], [-1,-1,1], [-1,-1,-1]]

  rotations.each do |r|
    inversions.each do |i|
      transformed = [beacon[r[0]] * i[0], beacon[r[1]] * i[1], beacon[r[2]] * i[2]]
      transforms << [transformed, [r, i]]
    end
  end

  transforms
end

def relative_scanner_position(triangles_a, triangles_b)
  intersection = triangles_a.keys.to_set.intersection(triangles_b.keys.to_set).to_a
  return nil if intersection.empty?

  candidates = Hash.new(0)

  until intersection.empty?
    t = intersection.shift

    ta = triangles_a[t].to_a
    tb = triangles_b[t].to_a

    # for each beacon in triangle_b, try every transform and check the computed difference to each beacon in triangle_a
    tb.each_with_index do |tb_beacon, i|
      transforms = beacon_transformations(tb_beacon)
      transforms.each do |pos, ri| # transformed beacon pos, [rotation, inversion]
        c = [ta[i][0] - pos[0], ta[i][1] - pos[1], ta[i][2] - pos[2]]
        candidates[[c, ri]] += 1
      end
    end
  end

  # returns [scanner_position, [rotation, inversion]]
  candidates.invert[candidates.values.max]
end

def transform_scan(transform, scan)
  rotation, inversion = transform
  scan.map do |beacon|
    [beacon[rotation[0]] * inversion[0],
     beacon[rotation[1]] * inversion[1],
     beacon[rotation[2]] * inversion[2]]
  end
end

def offset_scan(offset, scan)
  scan.map do |beacon|
    [beacon[0] + offset[0],
     beacon[1] + offset[1],
     beacon[2] + offset[2]]
  end
end

# [offset, [rotation, inversion]]
def apply_relative_position(relative_position, scan)
  offset_scan(relative_position[0], transform_scan(relative_position[1], scan))
end

def triangles_overlap?(triangles_a, triangles_b)
  intersection = triangles_a.keys.to_set.intersection(triangles_b.keys.to_set)
  intersection.size >= 220
end

# takes position+transform pairs eg a=[pos, [rotation, inversion]]
def resolve_relative_positions(a, b)
  pos_a = a[0]
  inv_a = a[1][1]

  pos_b = b[0]
  inv_b = b[1][1]

  [pos_a[0] + pos_b[0] * inv_a[0],
   pos_a[1] + pos_b[1] * inv_a[1],
   pos_a[2] + pos_b[2] * inv_a[2]]
end

def find_path_to_0(relative_positions, origin, tgt=origin)
  return [0] if origin == 0

  graph = relative_positions.keys.reduce(Hash.new { Set.new }) { |g,k| g[k[0]] += [k[1]]; g }
  find_path(graph, origin, 0)
end

def find_path(graph, from, to)
  explored = Set.new([from])
  open = [from]
  path = {}
  until open.empty?
    current = open.shift
    break if current == to

    neighbors = graph[current].reject { explored.include?(_1) }
    neighbors.each do |n|
      explored += [n]
      path[n] = current
      open << n
    end
  end

  ret = [to]

  until to == from
    to = path[to]
    ret << to
  end

  ret.reverse
end

def compute_p1(input)
  scans = input.split(/--- scanner (\d+) ---/)[1..].each_slice(2).map { |n, s|
    scanner_idx = n.to_i
    reports = s.lines.map(&:chomp).reject(&:empty?).map { _1.split(',').map(&:to_i) }
    [scanner_idx, reports]
  }.to_h

  known_beacons = scans[0].to_set

  triangulated_scans = scans.map { |id,scan| {id: id, triangles: beacon_triangles(scan)} }

  scanner_positions = { 0 => [0,0,0] }
  scanner_transforms = { 0 => [[0,1,2], [1,1,1]] }

  relative_positions_a = triangulated_scans.combination(2)
                         .filter { |a,b| triangles_overlap?(a[:triangles], b[:triangles]) }
                         .map { |a,b| [[a[:id], b[:id]], relative_scanner_position(a[:triangles], b[:triangles])] }
  relative_positions_b = triangulated_scans.combination(2)
                         .filter { |a,b| triangles_overlap?(a[:triangles], b[:triangles]) }
                         .map { |b,a| [[a[:id], b[:id]], relative_scanner_position(a[:triangles], b[:triangles])] }
  relative_positions = (relative_positions_a + relative_positions_b).to_h


  # resolve scans
  resolved_scans = {0 => scans[0].to_set}
  (1...scans.size).each do |i|
    location = i
    scan = scans[i]
    path = find_path_to_0(relative_positions, i)[1..]

    until path[0] == 0
      scan = apply_relative_position(relative_positions[[path[0],location]], scan)
      location = path.shift
    end

    resolved_scans[i] = apply_relative_position(relative_positions[[0,location]], scan)
  end

  resolved_scans.values.inject(&:+).size
end

This puzzle is kind of like a sci-fi astronavigation problem that gets referenced fairly often; given some set of points (stars) in 3D space, how can I figure out my position relative to a known origin? I may drift back and forth between referring to the beacons as "beacons" or "stars" interchangeably because that's how I was thinking of it the whole time, so that's what's going on there.

To do this, we need to look at the relationships between the stars (beacons here) and find similarities between the unknown (results from any scanner other than 0) and the known (results from scanner 0, since we're treating that as the origin). Initially I was trying to get by with just looking at the distance to each pair of stars in a given scan, and then finding pairs with the same distances in other scans.

This sort of worked, but was giving me trouble with the rotations and I doubt it would've worked in the full input; so I ended up going with triangles instead.

The slowest part of my solution is to go through each scan, find every combination of three stars (A,B,C), compute the distance in three dimensional space between each of them (AB,BC,CA), and sort those distances into an array. Once I have this normalized triangle, I store that as the key to a hash, with the three stars as the value, something like { [1,2,3] => [[x1,y1,z1], [x2,y2,z2], [x3,y3,z3]] }

I do this for every combination of three stars in every scan, which ends up being about 2,300 triangles per scan.

Once I have this table of triangles (and the stars involved in each) for every scan, I can start looking for intersections. Any intersection will give me a set of at least three stars that I can work out a transformation for.

If I were a better mathematician, I'd probably know the matrix math to work this out more directly, but since I'm only looking at three stars at a time, I can just try everything and it works out well enough.

This is the relative_scanner_position function above, it takes a pair of triangle maps (A,B), finds the triangles that are present in both, and then tries every possible transformation on each star in each triangle from B. Each time a transformed star location lines up with one of the stars from A, we add one point to that transform.

At the end, we find the highest scoring transform and return that as the relationship between scanner A and scanner B.

So now we can take any given pair of scans, find out if they overlap at all, and if so what the relationship is between the two.

We do this for all overlapping scans, and as we go we're building up two important pieces of information: 1) the graph of connected/overlapping scans, and 2) a table of transforms that let us move stars from one frame of reference to another.

Initially I tried to directly flatten all the transforms to get everything down into the coordinate space of scanner 0, but that was giving me trouble so I ended up using BFS to pathfind my way to a list of transformations to apply to each scan. Fortunately the graph is small enough that this doesn't add a meaningful amount of time.

Then it's just a matter of flattening all the stars into one coordinate space, throwing them in a set, and counting.

I actually had, at one point, tracked down all the scanner locations (relative to scanner 0) during the part 1 solve, but ended up throwing them out at some point by the time I actually finished part 1. So I had to re-solve for that, but that ends up being as easy as mapping [0,0,0] down through the transform path with all the other stars.

I could improve the runtime by re-using all the work from part 1, but instead I just copy-pasted everything again ¯\_(ツ)_/¯

I enjoyed this enough that I'd like to come back and clean up the code and description a bit, but I'm out of time before day 20 and I've got holiday travel coming up, so this is what we get for now!

def compute_p2(input)
  scans = input.split(/--- scanner (\d+) ---/)[1..].each_slice(2).map { |n, s|
    scanner_idx = n.to_i
    reports = s.lines.map(&:chomp).reject(&:empty?).map { _1.split(',').map(&:to_i) }
    [scanner_idx, reports]
  }.to_h

  known_beacons = scans[0].to_set

  triangulated_scans = scans.map { |id,scan| {id: id, triangles: beacon_triangles(scan)} }

  scanner_positions = { 0 => [0,0,0] }
  scanner_transforms = { 0 => [[0,1,2], [1,1,1]] }

  relative_positions_a = triangulated_scans.combination(2)
                           .filter { |a,b| triangles_overlap?(a[:triangles], b[:triangles]) }
                           .map { |a,b| [[a[:id], b[:id]], relative_scanner_position(a[:triangles], b[:triangles])] }
  relative_positions_b = triangulated_scans.combination(2)
                           .filter { |a,b| triangles_overlap?(a[:triangles], b[:triangles]) }
                           .map { |b,a| [[a[:id], b[:id]], relative_scanner_position(a[:triangles], b[:triangles])] }
  relative_positions = (relative_positions_a + relative_positions_b).to_h

  # resolve scans and positions
  resolved_scans = {0 => scans[0].to_set}
  (1...scans.size).each do |i|
    location = i
    scan = scans[i]
    pos = [[0,0,0]]
    path = find_path_to_0(relative_positions, i)[1..]

    until path[0] == 0
      scan = apply_relative_position(relative_positions[[path[0],location]], scan)
      pos =  apply_relative_position(relative_positions[[path[0],location]], pos)
      location = path.shift
    end

    resolved_scans[i] = apply_relative_position(relative_positions[[0,location]], scan)
    scanner_positions[i] = apply_relative_position(relative_positions[[0,location]], pos).first
  end

  def manhattan(a,b)
    a.zip(b).map { |i,j| (i-j).abs }.sum
  end

  scanner_positions.values.combination(2).map { |a,b| manhattan(a,b) }.max
end

```

</details>