Ergodic_moves_3.hpp

Go to the documentation of this file.

/*******************************************************************************
 Causal Dynamical Triangulations in C++ using CGAL
 
 Copyright © 2019 Adam Getchell
 ******************************************************************************/
 
 
#ifndef CDT_PLUSPLUS_ERGODIC_MOVES_3_HPP
#define CDT_PLUSPLUS_ERGODIC_MOVES_3_HPP
 
#include <expected>
 
#include "Formatters.hpp"
#include "Manifold.hpp"
#include "Move_tracker.hpp"
 
namespace ergodic_moves
{
  using Manifold         = manifolds::Manifold_3;
  using Expected         = std::expected<Manifold, std::string>;
  using Cell_handle      = Cell_handle_t<3>;
  using Cell_container   = std::vector<Cell_handle>;
  using Edge_handle      = Edge_handle_t<3>;
  using Edge_container   = std::vector<Edge_handle>;
  using Vertex_handle    = Vertex_handle_t<3>;
  using Vertex_container = std::vector<Vertex_handle>;
  using Delaunay         = Delaunay_t<3>;
 
  [[nodiscard]] inline auto null_move(Manifold const& t_manifold) noexcept
      -> Expected
  {
    return t_manifold;
  }  // null_move
 
  [[nodiscard]] inline auto try_23_move(Manifold&          t_manifold,
                                        Cell_handle const& to_be_moved) -> bool
  {
    if (to_be_moved->info() != 22) { return false; }  // NOLINT
    auto flipped = false;
    // Try every facet of the (2,2) cell
    for (auto i = 0; i < 4; ++i)
    {
      if (t_manifold.triangulation().flip(to_be_moved, i))
      {
#ifndef NDEBUG
        spdlog::trace("Facet {} was flippable.\n", i);
#endif
        flipped = true;
        break;
      }
#ifndef NDEBUG
      spdlog::trace("Facet {} was not flippable.\n", i);
#endif
    }
    return flipped;
  }  // try_23_move
 
  [[nodiscard]] inline auto do_23_move(Manifold& t_manifold) -> Expected
  {
#ifndef NDEBUG
    spdlog::debug("{} called.\n", __PRETTY_FUNCTION__);
#endif
 
    auto two_two = t_manifold.get_triangulation().get_two_two();
    // Shuffle the container to create a random sequence of (2,2) cells
    std::ranges::shuffle(two_two, utilities::make_random_generator());
    // Try a (2,3) move on successive cells in the sequence
    if (std::ranges::any_of(
            two_two, [&](auto& cell) { return try_23_move(t_manifold, cell); }))
    {
      return t_manifold;
    }
    // We've run out of (2,2) cells
    std::string const msg = "No (2,3) move possible.\n";
    spdlog::warn(msg);
    return std::unexpected(msg);
  }
 
  [[nodiscard]] inline auto try_32_move(Manifold&          t_manifold,
                                        Edge_handle const& to_be_moved) -> bool
  {
    return t_manifold.triangulation().flip(
        to_be_moved.first, to_be_moved.second, to_be_moved.third);
  }  // try_32_move
 
  [[nodiscard]] inline auto do_32_move(Manifold& t_manifold) -> Expected
  {
#ifndef NDEBUG
    spdlog::debug("{} called.\n", __PRETTY_FUNCTION__);
#endif
    auto timelike_edges = t_manifold.get_timelike_edges();
    // Shuffle the container to create a random sequence of edges
    std::ranges::shuffle(timelike_edges, utilities::make_random_generator());
    // Try a (3,2) move on successive timelike edges in the sequence
    if (std::ranges::any_of(timelike_edges, [&](auto& edge) {
          return try_32_move(t_manifold, edge);
        }))
    {
      return t_manifold;
    }
    // We've run out of edges to try
    std::string const msg = "No (3,2) move possible.\n";
    spdlog::warn(msg);
    return std::unexpected(msg);
  }  // do_32_move()
 
  [[nodiscard]] inline auto find_adjacent_31_cell(Cell_handle const& t_cell)
      -> std::optional<int>
  {
    if (t_cell->info() != 13) { return std::nullopt; }  // NOLINT
    for (auto i = 0; i < 4; ++i)
    {
#ifndef NDEBUG
      spdlog::trace("Neighbor {} is of type {}\n", i,
                    t_cell->neighbor(i)->info());
#endif
      if (foliated_triangulations::expected_cell_type<3>(t_cell->neighbor(i)) ==
          Cell_type::THREE_ONE)
      {
        return std::make_optional(i);
      }
    }
    return std::nullopt;
  }  // find_26_move()
 
  [[nodiscard]] inline auto do_26_move(Manifold& t_manifold) -> Expected
  {
#ifndef NDEBUG
    spdlog::debug("{} called.\n", __PRETTY_FUNCTION__);
#endif
    static auto constexpr INCIDENT_CELLS_FOR_6_2_MOVE = 6;
    auto one_three = t_manifold.get_triangulation().get_one_three();
    // Shuffle the container to pick a random sequence of (1,3) cells to try
    std::ranges::shuffle(one_three, utilities::make_random_generator());
    for (auto const& bottom : one_three)
    {
      if (auto neighboring_31_index = find_adjacent_31_cell(bottom);
          neighboring_31_index)
      {
#ifndef NDEBUG
        spdlog::trace("neighboring_31_index is {}.\n", *neighboring_31_index);
#endif
        Cell_handle const top  = bottom->neighbor(neighboring_31_index.value());
        // Calculate the common face with respect to the bottom cell
        auto common_face_index = std::numeric_limits<int>::max();
        if (!bottom->has_neighbor(top, common_face_index))
        {
          std::string const msg = "Bottom cell does not have a neighbor.\n";
#ifndef NDEBUG
          spdlog::trace(msg);
#endif
          return std::unexpected(msg);
        }
 
        // Get indices of vertices of common face with respect to bottom cell
        // A face is denoted by the index of the opposite vertex
        // Thus, the indices of the vertices of the face are all other indices
        // except the common_face_index
        // CGAL uses bitwise operations like (common_face_index +1) & 3
        // We use % 4 which is equivalent and doesn't trigger clang-tidy
        auto const i_1 = (common_face_index + 1) % 4;
        auto const i_2 = (common_face_index + 2) % 4;
        auto const i_3 = (common_face_index + 3) % 4;
 
        // Get vertices of common face from indices
        auto const v_1 = bottom->vertex(i_1);
        auto const v_2 = bottom->vertex(i_2);
        auto const v_3 = bottom->vertex(i_3);
 
        // Timeslice of vertices should be same
        if (v_1->info() != v_2->info() || v_2->info() != v_3->info())
        {
          std::string const msg = "Vertices have different timeslices.\n";
#ifndef NDEBUG
          spdlog::trace(msg);
#endif
          return std::unexpected(msg);
        }
 
        // Do the (2,6) move
        // Insert new vertex
        Vertex_handle const v_center =
            t_manifold.triangulation().delaunay().tds().insert_in_facet(
                bottom, *neighboring_31_index);
 
        // Checks
        Cell_container incident_cells;
        t_manifold.triangulation().delaunay().tds().incident_cells(
            v_center, std::back_inserter(incident_cells));
        // the (2,6) center vertex should be bounded by 6 simplices
        if (incident_cells.size() != INCIDENT_CELLS_FOR_6_2_MOVE)
        {
          std::string const msg =
              "Center vertex is not bounded by 6 simplices.\n";
#ifndef NDEBUG
          spdlog::trace(msg);
#endif
          return std::unexpected(msg);
        }
 
        // Each incident cell should be combinatorially and geometrically valid
        if (auto check_cells =
                std::ranges::all_of(incident_cells,
                                    [&t_manifold](auto const& cell) {
                                      return t_manifold.get_triangulation()
                                          .get_delaunay()
                                          .tds()
                                          .is_cell(cell);
                                    });
            !check_cells)
        {
          std::string const msg = "A cell is invalid.\n";
#ifndef NDEBUG
          spdlog::trace(msg);
#endif
          return std::unexpected(msg);
        }
 
        // Now assign a geometric point to the center vertex
        auto center_point =
            CGAL::centroid(v_1->point(), v_2->point(), v_3->point());
#ifndef NDEBUG
        spdlog::trace("Center point is: ({}).\n",
                      utilities::point_to_str(center_point));
#endif
        v_center->set_point(center_point);
 
        // Assign a timevalue to the new vertex
        auto timevalue   = v_1->info();
        v_center->info() = timevalue;
 
#ifndef NDEBUG
        if (t_manifold.is_vertex(v_center))
        {
          spdlog::trace("It's a vertex in the TDS.\n");
        }
        else { spdlog::trace("It's not a vertex in the TDS.\n"); }
        spdlog::trace("Spacelike face timevalue is {}.\n", timevalue);
        spdlog::trace("Inserted vertex ({}) with timevalue {}.\n",
                      utilities::point_to_str(v_center->point()),
                      v_center->info());
#endif
 
        // Final checks
        // is_valid() checks for combinatorial and geometric validity
        if (!t_manifold.get_delaunay().tds().is_valid(v_center, true, 1))
        {
          std::string const msg = "v_center is invalid.\n";
#ifndef NDEBUG
          spdlog::trace(msg);
#endif
          return std::unexpected(msg);
        }
 
        return t_manifold;
      }
      // Try next cell
#ifndef NDEBUG
      spdlog::debug("Cell not insertable.\n");
#endif
    }
    // We've run out of (1,3) simplices to try
    std::string const msg = "No (2,6) move possible.\n";
    spdlog::warn(msg);
    return std::unexpected(msg);
  }  // do_26_move()
 
  [[nodiscard]] inline auto is_62_movable(Manifold const&      manifold,
                                          Vertex_handle const& candidate)
      -> bool
  {
    if (manifold.dimensionality() != 3)
    {
#ifndef NDEBUG
      spdlog::trace("Manifold is not 3-dimensional.\n");
#endif
      return false;
    }
 
    if (!manifold.is_vertex(candidate))
    {
#ifndef NDEBUG
      spdlog::trace("Candidate is not a vertex.\n");
#endif
      return false;
    }
 
    // We must have 5 incident edges to have 6 incident cells
    if (auto incident_edges = manifold.degree(candidate);
        incident_edges != 5)  // NOLINT
    {
#ifndef NDEBUG
      spdlog::trace("Vertex has {} incident edges.\n", incident_edges);
#endif
      return false;
    }
 
    // Obtain all incident cells
    auto const incident_cells = manifold.incident_cells(candidate);
 
    // We must have 6 cells incident to the vertex to make a (6,2) move
    if (incident_cells.size() != 6)  // NOLINT
    {
#ifndef NDEBUG
      spdlog::trace("Vertex has {} incident cells.\n", incident_cells.size());
#endif
      return false;
    }
 
    // Check that none of the incident cells are infinite
    for (auto const& cell : incident_cells)
    {
      if (manifold.get_triangulation().is_infinite(cell))
      {
#ifndef NDEBUG
        spdlog::trace("Cell is infinite.\n");
#endif
        return false;
      }
    }
 
    // Run until all vertices are fixed
    while (foliated_triangulations::fix_vertices<3>(
        manifold.get_delaunay(), manifold.initial_radius(),
        manifold.foliation_spacing()))
    {
      spdlog::warn("Fixing vertices found by is_62_movable().\n");
    }
 
    // Run until all cells fixed or 50 passes
    //    for (auto passes = 1; passes < foliated_triangulations::MAX_FIX_PASSES
    //    + 1;
    //         ++passes)
    //    {
    //      if (!foliated_triangulations::fix_cells<3>(manifold.get_delaunay()))
    //      {
    //        break;
    //      }
    //      spdlog::warn("Fixing cells found by is_62_movable() pass {}.\n",
    //      passes);
    //    }
 
    auto const incident_31 = foliated_triangulations::filter_cells<3>(
        incident_cells, Cell_type::THREE_ONE);
    auto const incident_22 = foliated_triangulations::filter_cells<3>(
        incident_cells, Cell_type::TWO_TWO);
    auto const incident_13 = foliated_triangulations::filter_cells<3>(
        incident_cells, Cell_type::ONE_THREE);
 
    // All cells should be classified
    if (incident_13.size() + incident_22.size() + incident_31.size() !=
        6)  // NOLINT
    {
      spdlog::warn("Some incident cells on this vertex need to be fixed.\n");
    }
 
#ifndef NDEBUG
    spdlog::trace(
        "Vertex has {} incident cells with {} incident (3,1) simplices and {} "
        "incident (2,2) simplices and {} incident (1,3) simplices.\n",
        incident_cells.size(), incident_31.size(), incident_22.size(),
        incident_13.size());
    foliated_triangulations::debug_print_cells<3>(std::span{incident_cells});
#endif
    return incident_31.size() == 3 && incident_22.empty() &&
           incident_13.size() == 3;
 
  }  // find_62_moves()
 
  [[nodiscard]] inline auto do_62_move(Manifold& t_manifold) -> Expected
  {
#ifndef NDEBUG
    spdlog::debug("{} called.\n", __PRETTY_FUNCTION__);
#endif
    auto vertices = t_manifold.get_vertices();
    // Shuffle the container to create a random sequence of vertices
    std::ranges::shuffle(vertices, utilities::make_random_generator());
    // Try a (6,2) move on successive vertices in the sequence
    if (auto const movable_vertex_iterator =
            std::ranges::find_if(vertices,
                                 [&](auto const& vertex) {
                                   return is_62_movable(t_manifold, vertex);
                                 });
        movable_vertex_iterator != vertices.end())
    {
      t_manifold.triangulation().delaunay().remove(*movable_vertex_iterator);
      return t_manifold;
    }
    // We've run out of vertices to try
    std::string const msg = "No (6,2) move possible.\n";
    spdlog::warn(msg);
    return std::unexpected(msg);
  }  // do_62_move()
 
  [[nodiscard]] inline auto incident_cells_from_edge(
      Delaunay_t<3> const& triangulation, Edge_handle const& edge)
      -> std::optional<Cell_container>
  {
    if (!triangulation.tds().is_edge(edge.first, edge.second, edge.third))
    {
      return std::nullopt;
    }
    // Create the circulator of cells around the edge, starting with the cell
    // containing the edge
    auto           circulator = triangulation.incident_cells(edge, edge.first);
    Cell_container incident_cells;
    // Add cells to the container until we get back to the first one in the
    // circulator
    do {  // NOLINT(cppcoreguidelines-avoid-do-while)
      // Ignore cells containing the infinite vertex
      if (triangulation.is_infinite(circulator)) { continue; }
      incident_cells.emplace_back(circulator);
    }
    while (++circulator != edge.first);
#ifndef NDEBUG
    spdlog::trace("Found {} incident cells on edge.\n", incident_cells.size());
#endif
    return incident_cells;
  }  // incident_cells_from_edge()
 
  [[nodiscard]] inline auto find_bistellar_flip_location(
      Delaunay const& triangulation, Edge_handle const& t_edge_candidate)
      -> std::optional<Cell_container>
  {
    if (auto incident_cells =
            incident_cells_from_edge(triangulation, t_edge_candidate);
        incident_cells.has_value() && incident_cells->size() == 4)
    {
      return incident_cells.value();
    }
    return std::nullopt;
  }  // find_bistellar_flip_location()
 
  [[nodiscard]] inline auto get_incident_cells(Delaunay const&   triangulation,
                                               Edge_handle const edge)
      -> std::optional<Cell_container>
  {
    // Check that the edge is valid
    if (!triangulation.tds().is_edge(edge.first, edge.second, edge.third))
    {
      return std::nullopt;
    }
 
    Cell_container incident_cells;
    auto           circulator = triangulation.incident_cells(edge, edge.first);
    do {  // NOLINT(cppcoreguidelines-avoid-do-while)
      // filter out boundary edges with incident infinite cells
      if (triangulation.is_infinite(circulator)) { continue; }
      incident_cells.emplace_back(circulator);
    }
    while (++circulator != edge.first);
 
    return incident_cells;
  }  // get_incident_cells()
 
  [[nodiscard]] inline auto bistellar_flip(Delaunay    triangulation,
                                           Edge_handle edge, Vertex_handle top,
                                           Vertex_handle bottom)
      -> std::optional<Delaunay>
  {
    // Get the cells incident to the edge
    auto incident_cells = get_incident_cells(triangulation, edge);
 
    // Check that there are exactly 4 incident cells
    if (!incident_cells || incident_cells->size() != 4) { return std::nullopt; }
 
    // Check incident cells are valid
    if (std::ranges::any_of(*incident_cells,
                            [](auto const& cell) { return !cell->is_valid(); }))
    {
      return std::nullopt;
    }
 
    // Get vertices from pivot edge
    auto const& pivot_from_1 = edge.first->vertex(edge.second);
    auto const& pivot_from_2 = edge.first->vertex(edge.third);
 
    // Get vertices from cells
    auto vertices = foliated_triangulations::get_vertices_from_cells<3>(
        incident_cells.value());
 
    // Get vertices for new pivot edge
    Vertex_container new_pivot_vertices;
    std::ranges::copy_if(vertices, std::back_inserter(new_pivot_vertices),
                         [&](auto const& vertex) {
                           return vertex != pivot_from_1 &&
                                  vertex != pivot_from_2 && vertex != top &&
                                  vertex != bottom;
                         });
 
    // Check that there are exactly 2 new pivot vertices
    if (new_pivot_vertices.size() != 2) { return std::nullopt; }
 
    // Label the vertices in the new pivot edge
    auto const& pivot_to_1 = new_pivot_vertices[0];
    auto const& pivot_to_2 = new_pivot_vertices[1];
 
    // Verify all vertices exist and are not infinite
    if (triangulation.is_infinite(pivot_from_1) ||
        triangulation.is_infinite(pivot_from_2) ||
        triangulation.is_infinite(pivot_to_1) ||
        triangulation.is_infinite(pivot_to_2) ||
        triangulation.is_infinite(top) || triangulation.is_infinite(bottom))
    {
      return std::nullopt;
    }
 
#ifndef NDEBUG
    spdlog::debug(
        "Pivot from: ({}, {}), Pivot to: ({}, {}), Top: {}, Bottom: {}",
        pivot_from_1->point(), pivot_from_2->point(), pivot_to_1->point(),
        pivot_to_2->point(), top->point(), bottom->point());
#endif
 
    // Now we need to classify the cells by the vertices they contain
    Cell_handle before_1 =
        nullptr;  // top, pivot_from_1, pivot_from_2, pivot_to_1
    Cell_handle before_2 =
        nullptr;  // top, pivot_from_1, pivot_from_2, pivot_to_2
    Cell_handle before_3 =
        nullptr;  // bottom, pivot_from_1, pivot_from_2, pivot_to_1
    Cell_handle before_4 =
        nullptr;  // bottom, pivot_from_1, pivot_from_2, pivot_to_2
 
    for (auto const& cell : incident_cells.value())
    {
      if (cell->has_vertex(top))
      {
        if (cell->has_vertex(pivot_to_1)) { before_1 = cell; }
        else { before_2 = cell; }
      }
      else
      {
        if (cell->has_vertex(pivot_to_1)) { before_3 = cell; }
        else { before_4 = cell; }
      }
    }
 
    // Verify all cells were found
    if (before_1 == nullptr || before_2 == nullptr || before_3 == nullptr ||
        before_4 == nullptr)
    {
#ifndef NDEBUG
      spdlog::debug("Not all cells were properly identified.");
#endif
      return std::nullopt;
    }
 
    // Verify cells are not infinite and are valid
    if (triangulation.is_infinite(before_1) ||
        triangulation.is_infinite(before_2) ||
        triangulation.is_infinite(before_3) ||
        triangulation.is_infinite(before_4) || !before_1->is_valid() ||
        !before_2->is_valid() || !before_3->is_valid() || !before_4->is_valid())
    {
#ifndef NDEBUG
      spdlog::debug("Some cells are invalid or infinite.");
#endif
      return std::nullopt;
    }
 
#ifndef NDEBUG
    spdlog::debug("Cells in the triangulation before deleting old cells: {}",
                  triangulation.number_of_cells());
 
    // Debug info about cells
    spdlog::debug("Before_1: vertices {}, {}, {}, {}",
                  before_1->vertex(0)->point(), before_1->vertex(1)->point(),
                  before_1->vertex(2)->point(), before_1->vertex(3)->point());
    spdlog::debug("Before_2: vertices {}, {}, {}, {}",
                  before_2->vertex(0)->point(), before_2->vertex(1)->point(),
                  before_2->vertex(2)->point(), before_2->vertex(3)->point());
    spdlog::debug("Before_3: vertices {}, {}, {}, {}",
                  before_3->vertex(0)->point(), before_3->vertex(1)->point(),
                  before_3->vertex(2)->point(), before_3->vertex(3)->point());
    spdlog::debug("Before_4: vertices {}, {}, {}, {}",
                  before_4->vertex(0)->point(), before_4->vertex(1)->point(),
                  before_4->vertex(2)->point(), before_4->vertex(3)->point());
#endif
 
    // Find the exterior neighbors of the cells
    Cell_handle n_1 = before_1->neighbor(before_1->index(pivot_from_2));
    Cell_handle n_2 = before_1->neighbor(before_1->index(pivot_from_1));
    Cell_handle n_3 = before_2->neighbor(before_2->index(pivot_from_1));
    Cell_handle n_4 = before_2->neighbor(before_2->index(pivot_from_2));
    Cell_handle n_5 = before_3->neighbor(before_3->index(pivot_from_2));
    Cell_handle n_6 = before_3->neighbor(before_3->index(pivot_from_1));
    Cell_handle n_7 = before_4->neighbor(before_4->index(pivot_from_1));
    Cell_handle n_8 = before_4->neighbor(before_4->index(pivot_from_2));
 
    // Store indices before deletion
    int n1_idx      = n_1->index(before_1);
    int n2_idx      = n_2->index(before_1);
    int n3_idx      = n_3->index(before_2);
    int n4_idx      = n_4->index(before_2);
    int n5_idx      = n_5->index(before_3);
    int n6_idx      = n_6->index(before_3);
    int n7_idx      = n_7->index(before_4);
    int n8_idx      = n_8->index(before_4);
 
    // Verify all neighbors exist and are not cells we're about to delete
    if (n_1 == nullptr || n_2 == nullptr || n_3 == nullptr || n_4 == nullptr ||
        n_5 == nullptr || n_6 == nullptr || n_7 == nullptr || n_8 == nullptr ||
        n_1 == before_1 || n_1 == before_2 || n_1 == before_3 ||
        n_1 == before_4 || n_2 == before_1 || n_2 == before_2 ||
        n_2 == before_3 || n_2 == before_4 || n_3 == before_1 ||
        n_3 == before_2 || n_3 == before_3 || n_3 == before_4 ||
        n_4 == before_1 || n_4 == before_2 || n_4 == before_3 ||
        n_4 == before_4 || n_5 == before_1 || n_5 == before_2 ||
        n_5 == before_3 || n_5 == before_4 || n_6 == before_1 ||
        n_6 == before_2 || n_6 == before_3 || n_6 == before_4 ||
        n_7 == before_1 || n_7 == before_2 || n_7 == before_3 ||
        n_7 == before_4 || n_8 == before_1 || n_8 == before_2 ||
        n_8 == before_3 || n_8 == before_4)
    {
#ifndef NDEBUG
      spdlog::debug("Issue with neighbors");
#endif
      return std::nullopt;
    }
 
    // Temporarily break neighbor relationships to avoid potential issues
    // This step is important to avoid having old cells pointing to new cells
    // during deletion
    triangulation.tds().set_adjacency(before_1, before_1->index(pivot_from_2),
                                      nullptr, -1);
    triangulation.tds().set_adjacency(before_1, before_1->index(pivot_from_1),
                                      nullptr, -1);
    triangulation.tds().set_adjacency(before_2, before_2->index(pivot_from_1),
                                      nullptr, -1);
    triangulation.tds().set_adjacency(before_2, before_2->index(pivot_from_2),
                                      nullptr, -1);
    triangulation.tds().set_adjacency(before_3, before_3->index(pivot_from_2),
                                      nullptr, -1);
    triangulation.tds().set_adjacency(before_3, before_3->index(pivot_from_1),
                                      nullptr, -1);
    triangulation.tds().set_adjacency(before_4, before_4->index(pivot_from_1),
                                      nullptr, -1);
    triangulation.tds().set_adjacency(before_4, before_4->index(pivot_from_2),
                                      nullptr, -1);
 
    // Next, delete the old cells
    triangulation.tds().delete_cell(before_1);
    triangulation.tds().delete_cell(before_2);
    triangulation.tds().delete_cell(before_3);
    triangulation.tds().delete_cell(before_4);
 
#ifndef NDEBUG
    spdlog::debug("Cells in the triangulation after deleting old cells: {}",
                  triangulation.number_of_cells());
#endif
 
    // Create new cells with consistent vertex ordering (counter-clockwise
    // orientation) Use a consistent ordering scheme for all tetrahedra
    Cell_handle after_1 = triangulation.tds().create_cell(
        pivot_to_1, pivot_to_2, pivot_from_1, top);
    Cell_handle after_2 = triangulation.tds().create_cell(
        pivot_to_1, pivot_to_2, pivot_from_2, top);
    Cell_handle after_3 = triangulation.tds().create_cell(
        pivot_to_1, pivot_to_2, pivot_from_1, bottom);
    Cell_handle after_4 = triangulation.tds().create_cell(
        pivot_to_1, pivot_to_2, pivot_from_2, bottom);
 
    // Verify each cell was created and is valid
    if (after_1 == nullptr || after_2 == nullptr || after_3 == nullptr ||
        after_4 == nullptr || !after_1->is_valid() || !after_2->is_valid() ||
        !after_3->is_valid() || !after_4->is_valid())
    {
#ifndef NDEBUG
      spdlog::debug("Failed to create new cells properly");
#endif
      return std::nullopt;
    }
 
    // Copy the cell type information from the original cells to maintain cell
    // type
    after_1->info() = before_1->info();
    after_2->info() = before_2->info();
    after_3->info() = before_3->info();
    after_4->info() = before_4->info();
 
#ifndef NDEBUG
    spdlog::debug("After_1: vertices {}, {}, {}, {}",
                  after_1->vertex(0)->point(), after_1->vertex(1)->point(),
                  after_1->vertex(2)->point(), after_1->vertex(3)->point());
    spdlog::debug("After_2: vertices {}, {}, {}, {}",
                  after_2->vertex(0)->point(), after_2->vertex(1)->point(),
                  after_2->vertex(2)->point(), after_2->vertex(3)->point());
    spdlog::debug("After_3: vertices {}, {}, {}, {}",
                  after_3->vertex(0)->point(), after_3->vertex(1)->point(),
                  after_3->vertex(2)->point(), after_3->vertex(3)->point());
    spdlog::debug("After_4: vertices {}, {}, {}, {}",
                  after_4->vertex(0)->point(), after_4->vertex(1)->point(),
                  after_4->vertex(2)->point(), after_4->vertex(3)->point());
#endif
 
    // Set up internal neighbor relationships between the new cells
    // Find the correct indices for each cell based on the vertex that's not
    // part of the shared face
    int after_1_after_2_idx = after_1->index(pivot_from_1);
    int after_2_after_1_idx = after_2->index(pivot_from_2);
 
    int after_1_after_3_idx = after_1->index(top);
    int after_3_after_1_idx = after_3->index(bottom);
 
    int after_2_after_4_idx = after_2->index(top);
    int after_4_after_2_idx = after_4->index(bottom);
 
    int after_3_after_4_idx = after_3->index(pivot_from_1);
    int after_4_after_3_idx = after_4->index(pivot_from_2);
 
    // Set up internal neighbor relationships with proper indices
    triangulation.tds().set_adjacency(after_1, after_1_after_2_idx, after_2,
                                      after_2_after_1_idx);
    triangulation.tds().set_adjacency(after_1, after_1_after_3_idx, after_3,
                                      after_3_after_1_idx);
    triangulation.tds().set_adjacency(after_2, after_2_after_4_idx, after_4,
                                      after_4_after_2_idx);
    triangulation.tds().set_adjacency(after_3, after_3_after_4_idx, after_4,
                                      after_4_after_3_idx);
 
    // Verify internal neighbor relationships are set up correctly
    if (!after_1->has_neighbor(after_2) || !after_2->has_neighbor(after_1) ||
        !after_1->has_neighbor(after_3) || !after_3->has_neighbor(after_1) ||
        !after_2->has_neighbor(after_4) || !after_4->has_neighbor(after_2) ||
        !after_3->has_neighbor(after_4) || !after_4->has_neighbor(after_3))
    {
#ifndef NDEBUG
      spdlog::debug("Internal neighbor relationships are not set up correctly");
#endif
      return std::nullopt;
    }
 
    // Find the correct indices for external neighbor relationships
    // These are the facets that are not shared with other new cells
    int after_1_ext_1_idx = after_1->index(pivot_to_2);
    int after_1_ext_2_idx = after_1->index(pivot_to_1);
 
    int after_2_ext_1_idx = after_2->index(pivot_to_2);
    int after_2_ext_2_idx = after_2->index(pivot_to_1);
 
    int after_3_ext_1_idx = after_3->index(pivot_to_2);
    int after_3_ext_2_idx = after_3->index(pivot_to_1);
 
    int after_4_ext_1_idx = after_4->index(pivot_to_2);
    int after_4_ext_2_idx = after_4->index(pivot_to_1);
 
    // Set up external neighbor relationships
    triangulation.tds().set_adjacency(after_1, after_1_ext_1_idx, n_1, n1_idx);
    triangulation.tds().set_adjacency(after_1, after_1_ext_2_idx, n_2, n2_idx);
    triangulation.tds().set_adjacency(after_2, after_2_ext_1_idx, n_3, n3_idx);
    triangulation.tds().set_adjacency(after_2, after_2_ext_2_idx, n_4, n4_idx);
    triangulation.tds().set_adjacency(after_3, after_3_ext_1_idx, n_5, n5_idx);
    triangulation.tds().set_adjacency(after_3, after_3_ext_2_idx, n_6, n6_idx);
    triangulation.tds().set_adjacency(after_4, after_4_ext_1_idx, n_7, n7_idx);
    triangulation.tds().set_adjacency(after_4, after_4_ext_2_idx, n_8, n8_idx);
 
    // Verify external neighbor relationships are set up correctly
    bool external_neighbors_ok = after_1->neighbor(after_1_ext_1_idx) == n_1 &&
                                 n_1->neighbor(n1_idx) == after_1 &&
                                 after_1->neighbor(after_1_ext_2_idx) == n_2 &&
                                 n_2->neighbor(n2_idx) == after_1 &&
                                 after_2->neighbor(after_2_ext_1_idx) == n_3 &&
                                 n_3->neighbor(n3_idx) == after_2 &&
                                 after_2->neighbor(after_2_ext_2_idx) == n_4 &&
                                 n_4->neighbor(n4_idx) == after_2 &&
                                 after_3->neighbor(after_3_ext_1_idx) == n_5 &&
                                 n_5->neighbor(n5_idx) == after_3 &&
                                 after_3->neighbor(after_3_ext_2_idx) == n_6 &&
                                 n_6->neighbor(n6_idx) == after_3 &&
                                 after_4->neighbor(after_4_ext_1_idx) == n_7 &&
                                 n_7->neighbor(n7_idx) == after_4 &&
                                 after_4->neighbor(after_4_ext_2_idx) == n_8 &&
                                 n_8->neighbor(n8_idx) == after_4;
 
    if (!external_neighbors_ok)
    {
#ifndef NDEBUG
      spdlog::debug("External neighbor relationships are not set up correctly");
#endif
      // Don't return null here - attempt to fix with reorientation
    }
 
#ifndef NDEBUG
    spdlog::debug("Cells in the triangulation after adding new cells: {}",
                  triangulation.number_of_cells());
 
    // Verify neighbor relationships
    spdlog::debug("after_1 neighbors after_2: {}",
                  after_1->has_neighbor(after_2));
    spdlog::debug("after_2 neighbors after_1: {}",
                  after_2->has_neighbor(after_1));
    spdlog::debug("after_1 neighbors after_3: {}",
                  after_1->has_neighbor(after_3));
    spdlog::debug("after_3 neighbors after_1: {}",
                  after_3->has_neighbor(after_1));
    spdlog::debug("after_2 neighbors after_4: {}",
                  after_2->has_neighbor(after_4));
    spdlog::debug("after_4 neighbors after_2: {}",
                  after_4->has_neighbor(after_2));
    spdlog::debug("after_3 neighbors after_4: {}",
                  after_3->has_neighbor(after_4));
    spdlog::debug("after_4 neighbors after_3: {}",
                  after_4->has_neighbor(after_3));
#endif
 
    // Attempt to fix orientation issues if any exist
    triangulation.tds().reorient();
 
    // Final validation
    if (!triangulation.tds().is_valid(true, 4))
    {
#ifndef NDEBUG
      spdlog::debug("Final triangulation is not valid after reorientation");
#endif
      return std::nullopt;
    }
 
    return std::make_optional(triangulation);
  }  // bistellar_flip()
 
  [[nodiscard]] inline auto find_pivot_edge(Delaunay const&       triangulation,
                                            Edge_container const& edges)
      -> std::optional<Edge_handle>
  {
    for (auto const& edge : edges)
    {
      auto circulator = triangulation.incident_cells(edge, edge.first);
      Cell_container incident_cells;
      do {  // NOLINT(cppcoreguidelines-avoid-do-while)
        // filter out boundary edges with incident infinite cells
        if (triangulation.is_infinite(circulator)) { continue; }
        incident_cells.emplace_back(circulator);
      }
      while (++circulator != edge.first);
      fmt::print("Edge has {} incident finite cells\n", incident_cells.size());
      if (incident_cells.size() == 4) { return edge; }
    }
    return std::nullopt;
  }  // find_pivot_edge()
 
  [[nodiscard]] inline auto get_vertices(Cell_container const& cells)
  {
    std::unordered_set<Vertex_handle> vertices;
    auto get_vertices = [&vertices](auto const& cell) {
      for (int i = 0; i < 4; ++i) { vertices.emplace(cell->vertex(i)); }
    };
    std::ranges::for_each(cells, get_vertices);
    Vertex_container result(vertices.begin(), vertices.end());
    return result;
  }  // get_vertices()
 
  [[nodiscard]] inline auto do_44_move(Manifold const& t_manifold) -> Expected
  {
#ifndef NDEBUG
    spdlog::debug("{} called.\n", __PRETTY_FUNCTION__);
#endif
    auto spacelike_edges = t_manifold.get_spacelike_edges();
    // Shuffle the container to pick a random sequence of edges to try
    std::ranges::shuffle(spacelike_edges, utilities::make_random_generator());
    for (auto const& edge : spacelike_edges)
    {
      // Obtain all incident cells
      if (auto const incident_cells = find_bistellar_flip_location(
              t_manifold.get_triangulation().get_delaunay(), edge);
          incident_cells)
      {
#ifndef NDEBUG
        for (auto const& cell : *incident_cells)
        {
          spdlog::trace("Incident cell is of type {}.\n", cell->info());
        }
#endif
 
        // Get edge vertices
        auto const& v1       = edge.first->vertex(edge.second);
        auto const& v2       = edge.first->vertex(edge.third);
 
        // Find top and bottom vertices
        Vertex_handle top    = nullptr;
        Vertex_handle bottom = nullptr;
 
        // Analyze cells to find the top and bottom vertices
        for (auto const& cell : *incident_cells)
        {
          for (int i = 0; i < 4; ++i)
          {
            auto vertex = cell->vertex(i);
            if (vertex != v1 && vertex != v2)
            {
              // Use timevalue to determine if it's a top or bottom vertex
              if (top == nullptr || vertex->info() > top->info())
              {
                top = vertex;
              }
              if (bottom == nullptr || vertex->info() < bottom->info())
              {
                bottom = vertex;
              }
            }
          }
        }
 
        // Try the bistellar flip
        if (auto flipped_triangulation =
                bistellar_flip(t_manifold.get_triangulation().get_delaunay(),
                               edge, top, bottom);
            flipped_triangulation.has_value())
        {
          // Create a new foliated triangulation with the flipped Delaunay
          // triangulation
          auto new_triangulation =
              foliated_triangulations::FoliatedTriangulation<3>(
                  flipped_triangulation.value(), t_manifold.initial_radius(),
                  t_manifold.foliation_spacing());
 
          // Create a new manifold with the updated triangulation
          auto new_manifold = Manifold(new_triangulation);
 
          return new_manifold;
        }
 
        // If we get here, the flip failed but we found potential cells
        // Log the reason and continue trying other edges
#ifndef NDEBUG
        spdlog::debug("Bistellar flip failed.");
#endif
      }
      // Try next edge
    }
    // We've run out of edges to try
    std::string const msg = "No (4,4) move possible.\n";
    spdlog::warn(msg);
    return std::unexpected(msg);
  }  // do_44_move()
 
  [[nodiscard]] inline auto check_move(Manifold const&                t_before,
                                       Manifold const&                t_after,
                                       move_tracker::move_type const& t_move)
      -> bool
  {
    switch (t_move)
    {
      case move_tracker::move_type::FOUR_FOUR:
        return t_after.is_valid() && t_after.N3() == t_before.N3() &&
               t_after.N3_31() == t_before.N3_31() &&
               t_after.N3_22() == t_before.N3_22() &&
               t_after.N3_13() == t_before.N3_13() &&
               t_after.N2() == t_before.N2() && t_after.N1() == t_before.N1() &&
               t_after.N1_TL() == t_before.N1_TL() &&
               t_after.N1_SL() == t_before.N1_SL() &&
               t_after.N0() == t_before.N0() &&
               t_after.max_time() == t_before.max_time() &&
               t_after.min_time() == t_before.min_time();
      case move_tracker::move_type::TWO_THREE:
        return t_after.is_valid() && t_after.N3() == t_before.N3() + 1 &&
               t_after.N3_31() == t_before.N3_31() &&
               t_after.N3_22() == t_before.N3_22() + 1 &&
               t_after.N3_13() == t_before.N3_13() &&
               t_after.N2() == t_before.N2() + 2 &&
               t_after.N1() == t_before.N1() + 1 &&
               t_after.N1_TL() == t_before.N1_TL() + 1 &&
               t_after.N1_SL() == t_before.N1_SL() &&
               t_after.N0() == t_before.N0() &&
               t_after.max_time() == t_before.max_time() &&
               t_after.min_time() == t_before.min_time();
      case move_tracker::move_type::THREE_TWO:
        return t_after.is_valid() && t_after.N3() == t_before.N3() - 1 &&
               t_after.N3_31() == t_before.N3_31() &&
               t_after.N3_22() == t_before.N3_22() - 1 &&
               t_after.N3_13() == t_before.N3_13() &&
               t_after.N2() == t_before.N2() - 2 &&
               t_after.N1() == t_before.N1() - 1 &&
               t_after.N1_TL() == t_before.N1_TL() - 1 &&
               t_after.N1_SL() == t_before.N1_SL() &&
               t_after.N0() == t_before.N0() &&
               t_after.max_time() == t_before.max_time() &&
               t_after.min_time() == t_before.min_time();
      case move_tracker::move_type::TWO_SIX:
        return t_after.is_valid() && t_after.N3() == t_before.N3() + 4 &&
               t_after.N3_31() == t_before.N3_31() + 2 &&
               t_after.N3_22() == t_before.N3_22() &&
               t_after.N3_13() == t_before.N3_13() + 2 &&
               t_after.N2() == t_before.N2() + 8 &&  // NOLINT
               t_after.N1() == t_before.N1() + 5 &&  // NOLINT
               t_after.N1_TL() == t_before.N1_TL() + 2 &&
               t_after.N1_SL() == t_before.N1_SL() + 3 &&
               t_after.N0() == t_before.N0() + 1 &&
               t_after.max_time() == t_before.max_time() &&
               t_after.min_time() == t_before.min_time();
      case move_tracker::move_type::SIX_TWO:
        return t_after.is_valid() && t_after.N3() == t_before.N3() - 4 &&
               t_after.N3_31() == t_before.N3_31() - 2 &&
               t_after.N3_22() == t_before.N3_22() &&
               t_after.N3_13() == t_before.N3_13() - 2 &&
               t_after.N2() == t_before.N2() - 8 &&  // NOLINT
               t_after.N1() == t_before.N1() - 5 &&  // NOLINT
               t_after.N1_TL() == t_before.N1_TL() - 2 &&
               t_after.N1_SL() == t_before.N1_SL() - 3 &&
               t_after.N0() == t_before.N0() - 1 &&
               t_after.max_time() == t_before.max_time() &&
               t_after.min_time() == t_before.min_time();
      default: return false;
    }
  }  // check_move()
 
}  // namespace ergodic_moves
 
#endif  // CDT_PLUSPLUS_ERGODIC_MOVES_3_HPP