From a7f3dced95a84ff668cd6a64607d30a0ad6a97e3 Mon Sep 17 00:00:00 2001 From: Josh Date: Sat, 26 Aug 2023 22:22:52 +1000 Subject: [PATCH] ext train line alg base --- kernel/kernel_code/symplectic_basis.c | 367 ++++++++++++-------------- 1 file changed, 176 insertions(+), 191 deletions(-) diff --git a/kernel/kernel_code/symplectic_basis.c b/kernel/kernel_code/symplectic_basis.c index 8b5d9c46..a0aeb2ff 100644 --- a/kernel/kernel_code/symplectic_basis.c +++ b/kernel/kernel_code/symplectic_basis.c @@ -15,18 +15,20 @@ #include "SnapPea.h" #include "kernel.h" #include "addl_code.h" +#include "kernel_prototypes.h" +#include "kernel_typedefs.h" #define ATLEAST_TWO(a, b, c) ((a) && (b)) || ((a) && (c)) || ((b) && (c)) #define TRI_TO_INDEX(tet_index, tet_vertex) (4 * (tet_index) + (tet_vertex)) -#define COPY_PATH_ENDPOINT(new, old) { \ - (new)->vertex = (old)->vertex; \ - (new)->face = (old)->face; \ - (new)->tri = (old)->tri; \ - (new)->region_index = (old)->region_index; \ - (new)->region = (old)->region; \ - (new)->node = (old)->node; \ - (new)->num_adj_curves = (old)->num_adj_curves; \ +#define COPY_PATH_ENDPOINT(new, old) { \ + (new)->vertex = (old)->vertex; \ + (new)->face = (old)->face; \ + (new)->tri = (old)->tri; \ + (new)->region_index = (old)->region_index; \ + (new)->region = (old)->region; \ + (new)->node = (old)->node; \ + (new)->num_adj_curves = (old)->num_adj_curves; \ } #define COPY_PATH_NODE(new, old) { \ @@ -304,9 +306,10 @@ Boolean *update_edge_classes_on_cusp(CuspStructure **, Boolean * void do_manifold_train_lines(CuspStructure **, Boolean *); void find_primary_train_line(CuspStructure *, Boolean *); Boolean *edge_classes_on_cusp(CuspStructure *, Boolean *); -PathEndPoint *do_initial_train_line_segment_on_cusp(CuspStructure *, PathEndPoint *, Boolean *); -PathEndPoint *do_train_line_segment_on_cusp(CuspStructure *, PathEndPoint *, Boolean *); -int next_valid_endpoint_index(PathEndPoint **, int, int); +void do_initial_train_line_segment_on_cusp(CuspStructure *, PathEndPoint *, PathEndPoint *); +void do_train_line_segment_on_cusp(CuspStructure *, PathEndPoint *, PathEndPoint *); +void extended_train_line_path(CuspStructure *, PathEndPoint *, PathEndPoint *, EdgeNode *, EdgeNode *, int); +PathEndPoint *next_valid_endpoint_index(CuspStructure *, PathEndPoint *); void tri_endpoint_to_region_endpoint(CuspStructure *, PathEndPoint *); int *target_endpoints(CuspStructure *, const Boolean *, int *); Boolean array_contains_true(Boolean *, int); @@ -314,8 +317,8 @@ void graph_path_to_path_node(Graph *, EdgeNode *, EdgeNode *, void split_cusp_regions_along_train_line_segment(CuspStructure *, PathNode *, PathNode *, PathEndPoint *, PathEndPoint *); void update_cusp_triangle_train_line_endpoints(CuspRegion *, CuspRegion *, CuspRegion *, PathNode *, PathEndPoint *, int); void split_cusp_region_train_line_endpoint(CuspRegion *, CuspRegion *, PathNode *, PathEndPoint *, int, int); -void extended_train_line_path(CuspStructure *, PathEndPoint *, EdgeNode *, EdgeNode *, int *, Boolean *, Boolean *, int, int); -void cycle_path(Graph *, EdgeNode *, EdgeNode *, Boolean *, Boolean *, int *, int, int, int, int, int); +void path_finding_with_loops(Graph *, int, int, Orientation, EdgeNode *, EdgeNode *); +void cycle_path(Graph *, EdgeNode *, EdgeNode *, int, int, int, int, int); /** * Construct Oscillating Curves and calculate holonomy @@ -2077,16 +2080,16 @@ void setup_train_lines_and_oscillating_curves(Triangulation *manifold, CuspStruc break; } -// pick edge classes for train lines + // pick edge classes for train lines for (edge = manifold->edge_list_begin.next; edge != &manifold->edge_list_end; edge = edge->next) { tet = edge->incident_tet; cusp1 = tet->cusp[one_vertex_at_edge[edge->incident_edge_index]]->index; cusp2 = tet->cusp[other_vertex_at_edge[edge->incident_edge_index]]->index; - if (multi_graph->edge_classes[edge->index] || multi_graph->e0 == edge->index) { -// (cusp1 == e0_cusp1 && cusp2 == e0_cusp2) || -// (cusp2 == e0_cusp1 && cusp1 == e0_cusp2) || -// !edge_exists(multi_graph->multi_graph, cusp1, cusp2)) { + if (multi_graph->edge_classes[edge->index] || + (cusp1 == e0_cusp1 && cusp2 == e0_cusp2) || + (cusp2 == e0_cusp1 && cusp1 == e0_cusp2) || + !edge_exists(multi_graph->multi_graph, cusp1, cusp2)) { edge_classes[edge->index] = TRUE; } else { edge_classes[edge->index] = FALSE; @@ -2377,34 +2380,28 @@ void do_manifold_train_lines(CuspStructure **cusps, Boolean *edge_classes) { */ void find_primary_train_line(CuspStructure *cusp, Boolean *edge_classes) { - int start_endpoint_index; + PathEndPoint *start, *finish; Triangulation *manifold = cusp->manifold; - PathEndPoint *start_endpoint = NULL; - for (int edge_class = 0; edge_class < manifold->num_tetrahedra; edge_class++) { - for (int edge_index = 0; edge_index < 2; edge_index++) { - if (cusp->train_line_endpoint[edge_index][edge_class].tri == NULL) - continue; - - start_endpoint = &cusp->train_line_endpoint[edge_index][edge_class]; - edge_classes[edge_index * manifold->num_tetrahedra + edge_class] = FALSE; - break; - } - - if (start_endpoint != NULL) - break; - } + start = next_valid_endpoint_index(cusp, NULL); + edge_classes[start->region_index] = FALSE; if (!array_contains_true(edge_classes, 2 * manifold->num_tetrahedra)) { return; } - tri_endpoint_to_region_endpoint(cusp, start_endpoint); - start_endpoint = do_initial_train_line_segment_on_cusp(cusp, start_endpoint, edge_classes); + finish = next_valid_endpoint_index(cusp, start); + edge_classes[finish->region_index] = FALSE; + tri_endpoint_to_region_endpoint(cusp, start); + tri_endpoint_to_region_endpoint(cusp, finish); + do_initial_train_line_segment_on_cusp(cusp, start, finish); while (array_contains_true(edge_classes, 2 * manifold->num_tetrahedra)) { - construct_cusp_region_dual_graph(cusp); - start_endpoint = do_train_line_segment_on_cusp(cusp, start_endpoint, edge_classes); + start = finish; + finish = next_valid_endpoint_index(cusp, start); + edge_classes[finish->region_index] = FALSE; + + do_train_line_segment_on_cusp(cusp, start, finish); } my_free(edge_classes); @@ -2440,12 +2437,11 @@ Boolean *edge_classes_on_cusp(CuspStructure *cusp, Boolean *edge_classes) { * line. */ -PathEndPoint *do_initial_train_line_segment_on_cusp(CuspStructure *cusp, PathEndPoint *start_endpoint, - Boolean *edge_classes) { +void do_initial_train_line_segment_on_cusp(CuspStructure *cusp, PathEndPoint *start_endpoint, + PathEndPoint *finish_endpoint) { EdgeNode node_begin, node_end; - PathEndPoint *finish_endpoint = NULL; Boolean *processed, *discovered; - int *targets, *parent, finish_region_index, len; + int *parent; node_begin.next = &node_end; node_begin.prev = NULL; @@ -2458,30 +2454,8 @@ PathEndPoint *do_initial_train_line_segment_on_cusp(CuspStructure *cusp, PathEnd discovered = NEW_ARRAY(cusp->dual_graph->num_vertices, Boolean); parent = NEW_ARRAY(cusp->dual_graph->num_vertices, int); - // update regions - for (int edge_class = 0; edge_class < cusp->num_edge_classes; edge_class++) { - for (int edge_index = 0; edge_index < 2; edge_index++) { - if (edge_classes[edge_index * cusp->num_edge_classes + edge_class]) - tri_endpoint_to_region_endpoint(cusp, &cusp->train_line_endpoint[edge_index][edge_class]); - } - } - - targets = target_endpoints(cusp, edge_classes, &len); init_search(cusp->dual_graph, processed, discovered, parent); - finish_region_index = bfs_target_list(cusp->dual_graph, start_endpoint->region_index, targets, - len, processed, discovered, parent); - find_path(start_endpoint->region_index, finish_region_index, parent, &node_begin, &node_end); - - for (int edge_class = 0; edge_class < cusp->num_edge_classes; edge_class++) { - for (int edge_index = 0; edge_index < 2; edge_index++) { - if (edge_classes[edge_index * cusp->num_edge_classes + edge_class] == FALSE || - cusp->train_line_endpoint[edge_index][edge_class].region_index != finish_region_index) - continue; - - finish_endpoint = &cusp->train_line_endpoint[edge_index][edge_class]; - edge_classes[edge_index * cusp->num_edge_classes + edge_class] = FALSE; - } - } + path_finding_with_loops(cusp->dual_graph, start_endpoint, finish_endpoint, &node_begin, &node_end); if (finish_endpoint == NULL) uFatalError("do_initial_train_line_segment_on_cusp", "symplectic_basis"); @@ -2489,7 +2463,6 @@ PathEndPoint *do_initial_train_line_segment_on_cusp(CuspStructure *cusp, PathEnd my_free(processed); my_free(discovered); my_free(parent); - my_free(targets); // split along curve graph_path_to_dual_curve(cusp->dual_graph, &node_begin, &node_end, @@ -2502,7 +2475,6 @@ PathEndPoint *do_initial_train_line_segment_on_cusp(CuspStructure *cusp, PathEnd finish_endpoint->node = cusp->train_line_path_end.prev; free_edge_node(&node_begin, &node_end); - return finish_endpoint; } /* @@ -2512,14 +2484,11 @@ PathEndPoint *do_initial_train_line_segment_on_cusp(CuspStructure *cusp, PathEnd * adjacent across the face of the cusp triangle we dive along. */ -PathEndPoint *do_train_line_segment_on_cusp(CuspStructure *cusp, PathEndPoint *start_endpoint, Boolean *edge_classes) { +void do_train_line_segment_on_cusp(CuspStructure *cusp, PathEndPoint *start_endpoint, PathEndPoint *finish_endpoint) { EdgeNode node_begin, node_end; PathNode *start_node; - PathEndPoint *finish_endpoint = NULL; CuspRegion *region; - Boolean *processed, *discovered; - int start, finish_region_index, start_index, visited, len; - int *targets, *parent; + int start, finish, start_index, visited; node_begin.next = &node_end; node_begin.prev = NULL; @@ -2527,20 +2496,8 @@ PathEndPoint *do_train_line_segment_on_cusp(CuspStructure *cusp, PathEndPoint *s node_end.prev = &node_begin; construct_cusp_region_dual_graph(cusp); - processed = NEW_ARRAY(cusp->dual_graph->num_vertices, Boolean); - discovered = NEW_ARRAY(cusp->dual_graph->num_vertices, Boolean); - parent = NEW_ARRAY(cusp->dual_graph->num_vertices, int); - - // update regions - for (int edge_class = 0; edge_class < cusp->num_edge_classes; edge_class++) { - for (int edge_index = 0; edge_index < 2; edge_index++) { - if (edge_classes[edge_index * cusp->num_edge_classes + edge_class]) - tri_endpoint_to_region_endpoint(cusp, &cusp->train_line_endpoint[edge_index][edge_class]); - } - } - - init_search(cusp->dual_graph, processed, discovered, parent); + // update start_endpoint region start_index = TRI_TO_INDEX(start_endpoint->tri->tet_index, start_endpoint->tri->tet_vertex); for (region = cusp->cusp_region_begin[start_index].next; region != &cusp->cusp_region_end[start_index]; @@ -2569,13 +2526,6 @@ PathEndPoint *do_train_line_segment_on_cusp(CuspStructure *cusp, PathEndPoint *s * it enters. Hence, we need to remove the edge of the dual graph * corresponding to the last curve segment we drew. This edge will be * added back when the dual graph is reconstructed. - * - * We find the next endpoint greedily, picking the first one we find. - * This reduces the chance of the curve not being able to find the next - * path endpoint. We can become disconnected from the rest of the graph - * if there are two adjacent cusp vertices we want to dive down, and - * the curve we find dives through one of the cusp vertices, passes around - * the edge and dives through the second cusp vertex. */ if (start_endpoint->face == cusp->train_line_path_end.prev->prev_face) { @@ -2591,33 +2541,13 @@ PathEndPoint *do_train_line_segment_on_cusp(CuspStructure *cusp, PathEndPoint *s visited = start_endpoint->region_index; } - delete_edge(cusp->dual_graph, visited, start, cusp->dual_graph->directed); - targets = target_endpoints(cusp, edge_classes, &len); - finish_region_index = bfs_target_list(cusp->dual_graph, start, targets, len, processed, discovered, parent); - if (finish_region_index == -1) - uFatalError("do_train_line_segment_on_cusp", "symplectic_basis"); - - find_path(start, finish_region_index, parent, &node_begin, &node_end); - - for (int edge_class = 0; edge_class < cusp->num_edge_classes; edge_class++) { - for (int edge_index = 0; edge_index < 2; edge_index++) { - if (edge_classes[edge_index * cusp->num_edge_classes + edge_class] == FALSE || - cusp->train_line_endpoint[edge_index][edge_class].region_index != finish_region_index) - continue; + finish = finish_endpoint->region_index; - finish_endpoint = &cusp->train_line_endpoint[edge_index][edge_class]; - edge_classes[edge_index * cusp->num_edge_classes + edge_class] = FALSE; - } - } + extended_train_line_path(cusp, start_endpoint, finish_endpoint, &node_begin, &node_end, visited); if (finish_endpoint == NULL) uFatalError("do_initial_train_line_segment_on_cusp", "symplectic_basis"); - my_free(processed); - my_free(discovered); - my_free(parent); - my_free(targets); - // split along curve start_node = cusp->train_line_path_end.prev; graph_path_to_path_node(cusp->dual_graph, &node_begin, &node_end, @@ -2628,20 +2558,145 @@ PathEndPoint *do_train_line_segment_on_cusp(CuspStructure *cusp, PathEndPoint *s finish_endpoint->node = cusp->train_line_path_end.prev; free_edge_node(&node_begin, &node_end); - return finish_endpoint; } -int next_valid_endpoint_index(PathEndPoint **endpoints, int num_endpoints, int current_index) { - int index; +PathEndPoint *next_valid_endpoint_index(CuspStructure *cusp, PathEndPoint *current_endpoint) { + int start_index, start_class, index, class; - for (index = current_index + 1; index < num_endpoints; index++) { - if (endpoints[index] == NULL || endpoints[index]->tri == NULL) - continue; + if (current_endpoint == NULL) { + start_index = 0; + start_class = -1; + } else { + start_index = current_endpoint->tri->vertices[current_endpoint->vertex].edge_index; + start_class = current_endpoint->tri->vertices[current_endpoint->vertex].edge_class; + } + + for (class = start_class + 1; class < cusp->num_edge_classes; class++) { + for (index = start_index; index < 2; index++) { + if (cusp->train_line_endpoint[index][class].tri == NULL) + continue; - return index; + return &cusp->train_line_endpoint[index][class]; + } } - return -1; + return NULL; +} + +/* + * Find a path from start_endpoint to finish_endpoint, which + * goes around a cycle so the center is on the same side as the face + * the finish endpoint dives along. + */ + +void path_finding_with_loops(Graph *g, int start, int finish, Orientation orientation, + EdgeNode *node_begin, EdgeNode *node_end) { + int cycle_start, cycle_end; + int *parent; + Boolean *discovered, *processed; + + processed = NEW_ARRAY(g->num_vertices, Boolean); + discovered = NEW_ARRAY(g->num_vertices, Boolean); + parent = NEW_ARRAY(g->num_vertices, int); + + init_search(g, processed, discovered, parent); + cycle = cycle_exists(g, start, processed, discovered, parent, &cycle_start, &cycle_end); + +} + +/* + * Find a path for the train line from start to finish endpoint + * and store in doubly linked list node_begin -> node_end. + * + * Uses the cycle ensured by path finding with loops to find + * a path if the finish endpoint is not in the subgraph. + */ + +void extended_train_line_path(CuspStructure *cusp, PathEndPoint *start_endpoint, PathEndPoint *finish_endpoint, + EdgeNode *node_begin, EdgeNode *node_end, int visited) { + int cycle_start, cycle_end; + Boolean cycle; + Boolean *discovered, *processed; + int *parent; + + processed = NEW_ARRAY(cusp->dual_graph->num_vertices, Boolean); + discovered = NEW_ARRAY(cusp->dual_graph->num_vertices, Boolean); + parent = NEW_ARRAY(cusp->dual_graph->num_vertices, int); + + init_search(cusp->dual_graph, processed, discovered, parent); + delete_edge(cusp->dual_graph, visited, start_endpoint->region_index, cusp->dual_graph->directed); + bfs(cusp->dual_graph, start_endpoint->region_index, processed, discovered, parent); + + if (parent[finish_endpoint->region_index] == -1 && start_endpoint->region_index != finish_endpoint->region_index) { + /* + * The finish endpoint is not in the subgraph we created by removing the edge + * (visited, start). Assume there exists a cycle in this subgraph, we use this to + * 'turn the curve around' and use the edge (visited, start). + */ + + init_search(cusp->dual_graph, processed, discovered, parent); + cycle = cycle_exists(cusp->dual_graph, start_endpoint->region_index, processed, discovered, parent, &cycle_start, &cycle_end); + + if (cycle == FALSE) + // nothing we can do, train line does not work + uFatalError("do_train_line_segment_on_cusp", "symplectic_basis"); + + // reset parent array + init_search(cusp->dual_graph, processed, discovered, parent); + bfs(cusp->dual_graph, start_endpoint->region_index, processed, discovered, parent); + + cycle_path(cusp->dual_graph, node_begin, node_end, start_endpoint->region_index, visited, + finish_endpoint->region_index, cycle_start, cycle_end); + } else { + path_finding_with_loops(cusp->dual_graph, start_endpoint->region_index, finish_endpoint->region_index, right_handed, node_begin, node_end); + } + + my_free(processed); + my_free(discovered); + my_free(parent); +} + +void cycle_path(Graph *g, EdgeNode *node_begin, EdgeNode *node_end, int start, int prev, int finish, + int cycle_start, int cycle_end) { + EdgeNode *node, *temp_node, temp_begin, temp_end; + Boolean *discovered, *processed; + int *parent; + + temp_begin.next = &temp_end; + temp_begin.prev = NULL; + temp_end.next = NULL; + temp_end.prev = &temp_begin; + + processed = NEW_ARRAY(g->num_vertices, Boolean); + discovered = NEW_ARRAY(g->num_vertices, Boolean); + parent = NEW_ARRAY(g->num_vertices, int); + + // find a path from start -> cycle_end + find_path(start, cycle_end, parent, node_begin, node_end); + + // duplicate the path start -> cycle_start, and reverse it + find_path(start, cycle_start, parent, &temp_begin, &temp_end); + for (node = temp_end.prev; node != &temp_begin; node = node->prev) { + temp_node = NEW_STRUCT( EdgeNode ); + temp_node->y = node->y; + INSERT_BEFORE(temp_node, node_end); + } + + // free temp path + while (temp_begin.next != &temp_end) { + node = temp_begin.next; + REMOVE_NODE(node); + my_free(node); + } + + // find a path from visited -> target + init_search(g, processed, discovered, parent); + bfs(g, prev, processed, discovered, parent); + find_path(prev, finish, parent, node_end->prev, node_end); + + my_free(processed); + my_free(discovered); + my_free(parent); } /* @@ -2930,76 +2985,6 @@ void split_cusp_region_train_line_endpoint(CuspRegion *region_end, CuspRegion *r INSERT_BEFORE(new_region, region_end); } -/* - * Potential solution when we cant find a path when constructing - * train lines. Exploits a cycle in the subgraph created when we - * removed the edge (visited, start). - */ - -void extended_train_line_path(CuspStructure *cusp, PathEndPoint *finish_endpoint, EdgeNode *node_begin, EdgeNode *node_end, - int *parent, Boolean *processed, Boolean *discovered, int start, int visited) { - int cycle_start, cycle_end; - Boolean cycle; - - if (parent[finish_endpoint->region_index] == -1 && start != finish_endpoint->region_index) { - /* - * The finish endpoint is not in the subgraph we created by removing the edge - * (visited, start). Assume there exists a cycle in this subgraph, we use this to - * 'turn the curve around' and use the edge (visited, start). - */ - - init_search(cusp->dual_graph, processed, discovered, parent); - cycle = cycle_exists(cusp->dual_graph, start, processed, discovered, parent, &cycle_start, &cycle_end); - - if (cycle == FALSE) - // nothing we can do, train line does not work - uFatalError("do_train_line_segment_on_cusp", "symplectic_basis"); - - // reset parent array - init_search(cusp->dual_graph, processed, discovered, parent); - bfs(cusp->dual_graph, start, processed, discovered, parent); - - cycle_path(cusp->dual_graph, node_begin, node_end, processed, discovered, parent, - start, visited, finish_endpoint->region_index, cycle_start, cycle_end); - } else { - find_path(start, finish_endpoint->region_index, parent, node_begin, node_end); - } - -} - -void cycle_path(Graph *g, EdgeNode *node_begin, EdgeNode *node_end, Boolean *processed, Boolean *discovered, - int *parent, int start, int prev, int finish, int cycle_start, int cycle_end) { - EdgeNode *node, *temp_node, temp_begin, temp_end; - - temp_begin.next = &temp_end; - temp_begin.prev = NULL; - temp_end.next = NULL; - temp_end.prev = &temp_begin; - - // find a path from start -> cycle_end - find_path(start, cycle_end, parent, node_begin, node_end); - - // duplicate the path start -> cycle_start, and reverse it - find_path(start, cycle_start, parent, &temp_begin, &temp_end); - for (node = temp_end.prev; node != &temp_begin; node = node->prev) { - temp_node = NEW_STRUCT( EdgeNode ); - temp_node->y = node->y; - INSERT_BEFORE(temp_node, node_end); - } - - // free temp path - while (temp_begin.next != &temp_end) { - node = temp_begin.next; - REMOVE_NODE(node); - my_free(node); - } - - // find a path from visited -> target - init_search(g, processed, discovered, parent); - bfs(g, prev, processed, discovered, parent); - find_path(prev, finish, parent, node_end->prev, node_end); -} - // ------------------------------------ /*