"risk level: {}",
        grid.minima().map(|p| grid[p] as u16 + 1).sum::<u16>()

        "Risk level: {}",
        grid.minima().map(|p| grid[p] as usize + 1).sum::<usize>()
    );
    let basins = grid.basins().iter().map(|(_, s)| *s).top_n(3);
    println!(
        "Top 3 basin sizes: {:?}; product: {}",
        basins,
        basins.iter().product::<usize>()

    fn map_with_coords<S, F: FnMut((usize, usize), &T) -> S>(
        &self,
        mut f: F,
    ) -> Grid<S> {
        let mut inner = Vec::with_capacity(self.inner.len());
        let outside = f((usize::MAX, usize::MAX), &self.outside);
        inner
            .extend(self.points().zip(self.inner.iter()).map(|(p, v)| f(p, v)));
        Grid {
            inner,
            outside,
            width: self.width,
        }
    }

impl<T: Ord + From<u8> + PartialEq<T>> Grid<T> {
    fn basins(&self) -> Vec<((usize, usize), usize)> {
        let mut uf = self
            .map_with_coords(|p, n| (p, if n == &T::from(9) { 0 } else { 1 }));
        for (p, h) in uf
            .points()
            .zip(self.inner.iter())
            .filter(|(_, h)| *h < &T::from(9))
        {
            // I originally thought this would also find the minimum of each
            // basin but it doesn't always.
            for (n, j) in uf
                .neighbours(p)
                .map(|n| (n, &self[n]))
                .filter(|(_, n)| *n < &T::from(9))
            {
                match h.cmp(j) {
                    std::cmp::Ordering::Less => uf.flow(n, p),
                    std::cmp::Ordering::Equal => (),
                    std::cmp::Ordering::Greater => uf.flow(p, n),
                }
            }
        }
        uf.points()
            .zip(uf.inner.iter().map(|(_, s)| s))
            .filter(|(_, s)| *s > &0)
            .map(|(p, s)| (p, *s))
            .collect()
    }
}
impl<T: Copy> Grid<((usize, usize), T)> {
    fn u_find(&mut self, p: (usize, usize)) -> (usize, usize) {
        let d = unsafe { &mut *(&mut self[p].0 as *mut (usize, usize)) };
        if d == &p {
            p
        } else {
            let q = self.u_find(*d);
            *d = q;
            q
        }
    }
}
impl<T: Copy + From<u8> + std::ops::Add<T, Output = T>>
    Grid<((usize, usize), T)>
{
    fn flow(&mut self, mut s: (usize, usize), mut t: (usize, usize)) {
        s = self.u_find(s);
        t = self.u_find(t);
        if s != t {
            let y = unsafe { &mut *(&mut self[t].1 as *mut T) };
            let (d, x) = &mut self[s];
            let x = std::mem::replace(x, T::from(0));
            *y = *y + x;
            *d = t;
        }
    }
}

    }
}
trait IterExtra: Iterator {
    fn top_n(&mut self, size: usize) -> Vec<Self::Item>
    where
        Self::Item: PartialOrd,
    {
        let mut result: Vec<_> = self.take(size).collect();
        for mut i in self {
            let i = &mut i;
            for j in result.iter_mut() {
                if i > j {
                    std::mem::swap(i, j);
                }
            }
        }
        result

impl<I: Iterator> IterExtra for I {}

                    Some((self.centre.0, self.centre.1 + 1))

                    Some((self.centre.0, y))