diff --git a/2022/16.cpp b/2022/16.cpp
index 160f41f..e1a58f2 100644
--- a/2022/16.cpp
+++ b/2022/16.cpp
@@ -1,21 +1,19 @@
+#include <bitset>
 #include <cstdint>
 #include <iostream>
 #include <sstream>
 #include <stdexcept>
 #include <string>
 #include <tuple>
-#include <vector>
-#include <map>
 #include <unordered_map>
-#include <limits>
+#include <vector>
 
-#include <bitset>
-
-#include <boost/spirit/home/qi.hpp>
-#include <boost/range/irange.hpp>
+#include <boost/multi_array.hpp>
 #include <boost/phoenix.hpp>
+#include <boost/range/adaptors.hpp>
 #include <boost/range/algorithm.hpp>
-#include <boost/range/adaptor/map.hpp>
+#include <boost/range/irange.hpp>
+#include <boost/spirit/include/qi.hpp>
 
 #include <doctest.h>
 
@@ -23,8 +21,13 @@
 
 namespace {
 
-auto ShortestDistances(std::vector<std::vector<std::uint64_t>> & dist)
+using distance_array = boost::multi_array<std::uint64_t, 2>;
+
+/// @brief Update distance matrix with transitive shortest paths
+/// @param dist Single-step distances between nodes (must be square)
+auto ShortestDistances(distance_array & dist) -> void
 {
+  // Floyd–Warshall_algorithm
   auto const range = boost::irange(dist.size());
   for (auto const k : range) {
     for (auto const i : range) {
@@ -37,20 +40,18 @@ auto ShortestDistances(std::vector<std::vector<std::uint64_t>> & dist)
   }
 }
 
-struct RawRoom {
+struct Room {
+  /// @brief Name of the room
   std::string name;
+  /// @brief Flow rate of the valve in the room
   std::uint64_t flow;
+  /// @brief Directly adjacent rooms
   std::vector<std::string> connections;
 };
 
-struct Room {
-  std::uint64_t flow;
-  std::vector<std::uint64_t> connections;
-};
-
-auto Parse(std::istream & in) -> std::vector<RawRoom>
+auto Parse(std::istream & in) -> std::vector<Room>
 {
-  std::vector<RawRoom> result;
+  std::vector<Room> result;
   std::string line;
   while (std::getline(in, line)) {
     namespace phx = boost::phoenix;
@@ -58,18 +59,17 @@ auto Parse(std::istream & in) -> std::vector<RawRoom>
     using namespace qi::labels;
     using It = std::string::const_iterator;
 
-    auto const name = qi::as_string[+qi::alpha];
-
-    qi::rule<It, RawRoom()> room_description =
+    qi::rule<It, std::string()> const name = qi::as_string[+qi::alpha];
+    qi::rule<It, Room()> const room_description =
       "Valve " >>
-      name [ phx::bind(&RawRoom::name, _val) = _1 ] >>
+      name [ phx::bind(&Room::name, _val) = _1 ] >>
       " has flow rate=" >>
-      qi::ulong_long [ phx::bind(&RawRoom::flow, _val) = _1 ] >>
+      qi::ulong_long [ phx::bind(&Room::flow, _val) = _1 ] >>
       "; tunnel"  >> -qi::string("s") >>
       " lead"     >> -qi::string("s") >>
       " to valve" >> -qi::string("s") >>
       " " >>
-      (name % ", ") [ phx::bind(&RawRoom::connections, _val) = _1 ];
+      (name % ", ") [ phx::bind(&Room::connections, _val) = _1 ];
 
     It b = line.begin();
     It e = line.end();
@@ -81,34 +81,40 @@ auto Parse(std::istream & in) -> std::vector<RawRoom>
   return result;
 }
 
-auto GenerateGraph(std::vector<RawRoom> const& raw_rooms) -> std::pair<std::size_t, std::vector<Room>>
+/// @brief Computes the distances between rooms and finds the address of the starting room
+/// @returns starting index and distances
+auto GenerateDistances(
+  std::vector<Room> const& rooms
+) -> std::pair<std::size_t, distance_array>
 {
-  std::pair<std::size_t, std::vector<Room>> result {};
-  std::unordered_map<std::string, std::size_t> names;
+  auto const N = rooms.size();
 
-  std::vector<std::vector<std::uint64_t>> distances(
-    raw_rooms.size(), std::vector<std::uint64_t>(raw_rooms.size(), raw_rooms.size()));
-  
-  for (auto const i : boost::irange(raw_rooms.size())) {
-    names[raw_rooms[i].name] = i;
-    if (raw_rooms[i].name == "AA") result.first = i;
+  // Associate the names and indexes of each room
+  std::unordered_map<std::string, std::size_t> names;
+  for (auto [i,room] : rooms | boost::adaptors::indexed()) {
+    names[room.name] = i;
   }
 
-  for (auto const i : boost::irange(raw_rooms.size())) {
-    auto & ds = distances[i];
-    ds[i] = 0;
-    for (auto const& name : raw_rooms[i].connections) {
-      ds[names[name]] = 1;
+  distance_array distances{boost::extents[N][N]};
+
+  // N is longer than any optimal distance by at least 1
+  std::fill_n(distances.data(), distances.num_elements(), N);
+
+  for (auto [i, room] : rooms | boost::adaptors::indexed()) {
+    auto di = distances[i];
+
+    // each room is one away from adjacent rooms
+    for (auto const& name : room.connections) {
+      di[names[name]] = 1;
     }
+
+    // zero distance to self
+    di[i] = 0;
   }
 
   ShortestDistances(distances);
 
-  for (auto const i : boost::irange(raw_rooms.size())) {
-    result.second.push_back({raw_rooms[i].flow, std::move(distances[i])});
-  }
-
-  return result;
+  return {names.at("AA"), std::move(distances)};
 }
 
 using Valves = std::bitset<64>;
@@ -120,57 +126,82 @@ struct State {
   Valves valves;
 };
 
+/// @brief Compute all the flows achievable with a set of values
+/// @param start Index of starting room
+/// @param initial_time Initial amount of time
+/// @param rooms Array of rooms from input file
+/// @param distances Shortest paths between all pairs of rooms
+/// @return Mapping of maximum flow achievable using a particular set of valves
 auto Routes(
   std::size_t const start,
-  std::vector<Room> const& rooms, std::uint64_t initial_time
+  std::uint64_t const initial_time,
+  std::vector<Room> const& rooms,
+  distance_array const& distances
 ) -> std::unordered_map<Valves, std::uint64_t>
 {
   std::vector<State> states { State{initial_time, 0, start, {}} };
   std::unordered_map<Valves, std::uint64_t> result;
-  
+
   while (!states.empty()) {
     auto const state = states.back();
     states.pop_back();
+    auto const i = state.location;
 
-    auto const& room = rooms[state.location];
-
-    for (auto const i : boost::irange(rooms.size())) {
+    for (auto [j, room] : rooms | boost::adaptors::indexed()) {
       // don't revisit a valve
-      if (state.valves.test(i)) { continue; }
+      if (state.valves.test(j)) { continue; }
 
       // don't visit rooms with useless valves
-      if (rooms[i].flow == 0) { continue; }
+      if (room.flow == 0) { continue; }
 
       // don't visit rooms we can't get to in time
-      auto const cost = room.connections[i];
+      auto const cost = distances[i][j];
       if (cost+1 >= state.time) { continue; }
 
-      auto const time = state.time - cost - 1;      
-      auto const flow = state.flow + rooms[i].flow * time;
+      auto const time = state.time - cost - 1;
+      auto const flow = state.flow + room.flow * time;
       auto valves = state.valves;
-      valves.set(i);
+      valves.set(j);
 
       // remember the best we've seen for this valve-set
       if (result[valves] < flow) {
         result[valves] = flow;
       }
 
-      states.push_back({time, flow, i, valves});
+      states.push_back({time, flow, static_cast<std::size_t>(j), valves});
     }
   }
 
   return result;
 }
 
-auto Part1(std::size_t start, std::vector<Room> const& rooms) -> std::uint64_t
+/// @brief Maximize the water flow using a single actor and 30 minutes
+/// @param start Index of the starting room
+/// @param rooms Rooms from input file
+/// @param distances Shortest distances between pairs of rooms
+/// @return Maximum flow achievable
+auto Part1(
+  std::size_t const start,
+  std::vector<Room> const& rooms,
+  distance_array const& distances
+) -> std::uint64_t
 {
-  auto const routes = Routes(start, rooms, 30);
+  auto const routes = Routes(start, 30, rooms, distances);
   return *boost::range::max_element(routes | boost::adaptors::map_values);
 }
 
-auto Part2(std::size_t start, std::vector<Room> const& rooms) -> std::uint64_t
+/// @brief Maximize the water flow using two actos and 26 minutes
+/// @param start Index of the starting room
+/// @param rooms Rooms from input file
+/// @param distances Shortest distances between pairs of rooms
+/// @return Maximum flow achievable
+auto Part2(
+  std::size_t const start,
+  std::vector<Room> const& rooms,
+  distance_array const& distances
+) -> std::uint64_t
 {
-  auto const routes = Routes(start, rooms, 26);
+  auto const routes = Routes(start, 26, rooms, distances);
   auto const end = routes.end();
   std::uint64_t result {0};
   for (auto it1 = routes.begin(); it1 != end; std::advance(it1, 1)) {
@@ -200,15 +231,17 @@ Valve HH has flow rate=22; tunnel leads to valve GG
 Valve II has flow rate=0; tunnels lead to valves AA, JJ
 Valve JJ has flow rate=21; tunnel leads to valve II
 )"};
-    auto const [start, graph] = GenerateGraph(Parse(in));
-    CHECK(1651 == Part1(start, graph));
-    CHECK(1707 == Part2(start, graph));
+    auto const rooms = Parse(in);
+    auto const [start, distances] = GenerateDistances(rooms);
+    CHECK(1651 == Part1(start, rooms, distances));
+    CHECK(1707 == Part2(start, rooms, distances));
   }
 }
 
 auto main(int argc, char** argv) -> int
 {
-  auto const [start, rooms] = GenerateGraph(Parse(*aocpp::Startup(argc, argv)));
-  std::cout << "Part 1: " << Part1(start, rooms) << std::endl;
-  std::cout << "Part 2: " << Part2(start, rooms) << std::endl;
+  auto const rooms = Parse(*aocpp::Startup(argc, argv));
+  auto const [start, distances] = GenerateDistances(rooms);
+  std::cout << "Part 1: " << Part1(start, rooms, distances) << std::endl;
+  std::cout << "Part 2: " << Part2(start, rooms, distances) << std::endl;
 }