From 2bb0740b68fdb62d45adc00204b3990ca42ade77 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 22 Aug 2016 10:11:14 -0400 Subject: started repo again without all the data files gunking the works --- src/graph_components.c | 185 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 185 insertions(+) create mode 100644 src/graph_components.c (limited to 'src/graph_components.c') diff --git a/src/graph_components.c b/src/graph_components.c new file mode 100644 index 0000000..036f2fe --- /dev/null +++ b/src/graph_components.c @@ -0,0 +1,185 @@ + +#include "fracture.h" + +bool set_connected(const cholmod_sparse *laplacian, unsigned int *marks, + int vertex, int label, int stop_at, int exclude) { + unsigned int num_verts = (CHOL_INT)laplacian->nrow; + + CHOL_INT *ia = (CHOL_INT *)laplacian->p; + CHOL_INT *ja = (CHOL_INT *)laplacian->i; + double *a = (double *)laplacian->x; + + bool stopped = false; + + int *queue = (int *)malloc((num_verts) * sizeof(int)); + assert(queue != NULL); + + int queue_pos = 0; + + int *component_list = (int *)malloc((num_verts) * sizeof(int)); + assert(component_list != NULL); + + int list_pos = 0; + queue[0] = vertex; + component_list[0] = vertex; + + unsigned int old_label = marks[vertex]; + marks[vertex] = label; + + while (queue_pos >= 0) { + int next_vertex = queue[queue_pos]; + queue_pos--; + for (int j = 0; j < ia[next_vertex + 1] - ia[next_vertex]; j++) { + if (ja[ia[next_vertex] + j] == stop_at && a[ia[next_vertex] + j] != 0) { + // if we run into the other vertex, the components haven't changed + stopped = true; queue_pos = -1; + break; + } + if (marks[ja[ia[next_vertex] + j]] != label && + a[ia[next_vertex] + j] != 0 && + ja[ia[next_vertex] + j] < num_verts - exclude) { + marks[ja[ia[next_vertex] + j]] = label; + queue_pos++; + queue[queue_pos] = ja[ia[next_vertex] + j]; + list_pos++; + component_list[list_pos] = ja[ia[next_vertex] + j]; + } + } + } + + if (stopped) { + for (unsigned int i = 0; i <= list_pos; i++) { + marks[component_list[i]] = old_label; + } + } + + free(queue); + free(component_list); + + if (stopped) + return false; + else + return true; +} + +int update_components(const cholmod_sparse *laplacian, unsigned int *marks, + int old_num_components, int v1, int v2, int exclude) { + assert(laplacian != NULL); + assert(marks != NULL); + + if (set_connected(laplacian, marks, v1, old_num_components + 1, v2, + exclude)) { + return old_num_components + 1; + } else + return old_num_components; +} + +unsigned int *find_components(const cholmod_sparse *laplacian, unsigned int skip) { + size_t num_verts = laplacian->nrow; + unsigned int *marks = (unsigned int *)calloc(num_verts, sizeof(unsigned int)); + + int num_components = 0; + + for (int i = 0; i < num_verts; i++) { + if (marks[i] == 0) { + num_components++; + set_connected(laplacian, marks, i, num_components, -1, skip); + } + } + + return marks; +} + +unsigned int find_cycles(unsigned int num_edges, const bool *fuses, const unsigned int *ev, const unsigned + int *vei, const unsigned int *ve, int **cycles) { + unsigned int num_cycles = 0; + unsigned int *elist = (unsigned int *)malloc(num_edges * sizeof(unsigned int)); + unsigned int *side = (unsigned int *)calloc(num_edges, sizeof(unsigned int)); + unsigned int *num = (unsigned int *)malloc(num_edges * sizeof(unsigned int)); + for (unsigned int i = 0; i < num_edges; i++) { + if (fuses[i]) { + bool in_cycle = false; + for (unsigned int j = 0; j < num_cycles; j++) { + unsigned int k = 0; + int e; + while ((e = cycles[j][k]) >= 0) { + if (e == i) { + in_cycle = true; + break; + } + k++; + } + if (in_cycle) { + break; + } + } + if (!in_cycle) { + unsigned int pos = 0; + unsigned int cur_e = i; + bool *included = (bool *)calloc(num_edges, sizeof(bool)); + included[cur_e] = true; + elist[0] = cur_e; + side[0] = 1; + num[0] = 0; + while (true) { + unsigned int cv = ev[2 * cur_e + side[pos]]; + unsigned int new_e = cur_e; + unsigned int j; + for (j = num[pos]; j < vei[cv+1] - vei[cv]; j++) { + new_e = ve[vei[cv] + j]; + if (new_e != cur_e && fuses[new_e]) break; + } + if (new_e != cur_e && fuses[new_e]) { + if (new_e == elist[0] && cv == ev[2 * elist[0] + !side[0]]) { + pos++; + cycles[2*num_cycles] = (int *)malloc((pos + 1) * sizeof(int)); + cycles[2*num_cycles+1] = (int *)malloc((pos + 1) * sizeof(int)); + for (unsigned int i = 0; i < pos; i++) { + cycles[2*num_cycles][i] = elist[i]; + cycles[2*num_cycles+1][i] = side[i]; + } + cycles[2*num_cycles][pos] = -1; + cycles[2*num_cycles+1][pos] = -1; + num_cycles++; + pos--; + num[pos]++; + included[pos] = false; + cur_e = elist[pos]; + } else if (included[new_e]) { + if (pos == 0) break; + included[cur_e] = false; + pos--; + num[pos]++; + cur_e = elist[pos]; + } else { + num[pos] = j; + pos++; + elist[pos] = new_e; + included[new_e] = true; + unsigned int nv1 = ev[2 * new_e]; + unsigned int nv2 = ev[2 * new_e + 1]; + unsigned int v = nv1 == cv ? nv2 : nv1; + side[pos] = v == nv1 ? 0 : 1; + num[pos] = 0; + cur_e = new_e; + } + } else { + if (pos == 0) break; + included[cur_e] = false; + pos--; + num[pos]++; + cur_e = elist[pos]; + } + } + free(included); + } + } + } + + free(elist); + free(side); + free(num); + + return num_cycles; +} + -- cgit v1.2.3-70-g09d2