started repo again without all the data files gunking the works

1 files changed, 185 insertions, 0 deletions

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;
+}
+