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/correlations.c | 345 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 345 insertions(+) create mode 100644 src/correlations.c (limited to 'src/correlations.c') diff --git a/src/correlations.c b/src/correlations.c new file mode 100644 index 0000000..bd05daf --- /dev/null +++ b/src/correlations.c @@ -0,0 +1,345 @@ + +#include "fracture.h" + +struct fibn { + unsigned int value; + unsigned int priority; + unsigned int rank; + bool marked; + struct fibn *first_child; + struct fibn *parent; + struct fibn *prev; + struct fibn *next; +}; + +typedef struct { + struct fibn *min; + unsigned int rank; + unsigned int trees; + struct fibn **vertex_to_node; +} fibh; + +unsigned int fib_findmin(fibh *heap) { return heap->min->value; } + +void ll_setparent(struct fibn *ll, struct fibn *parent) { + struct fibn *curnode = ll->next; + ll->parent = parent; + while (curnode != ll) { + curnode->parent = parent; + curnode = curnode->next; + } +} + +struct fibn *ll_merge(struct fibn *ll1, struct fibn *ll2) { + if (ll1 == NULL) + return ll2; + if (ll2 == NULL) + return ll1; + + // link the beginning of list one to the end of list two and vice versa + struct fibn *ll1_beg = ll1; + struct fibn *ll1_end = ll1->prev; + struct fibn *ll2_beg = ll2; + struct fibn *ll2_end = ll2->prev; + + ll1_beg->prev = ll2_end; + ll1_end->next = ll2_beg; + ll2_beg->prev = ll1_end; + ll2_end->next = ll1_beg; + + return ll1; +} + +struct fibn *ll_delroot(struct fibn *ll) { + if (ll == NULL) + return NULL; + + if (ll->next == ll) { + return NULL; + } + + struct fibn *ll_beg = ll->next; + struct fibn *ll_end = ll->prev; + + ll_beg->prev = ll_end; + ll_end->next = ll_beg; + + ll->next = ll; + ll->prev = ll; + + return ll_beg; +} + +void fib_insert(fibh *heap, unsigned int value, unsigned int priority) { + struct fibn *newnode = calloc(1, sizeof(struct fibn)); + newnode->value = value; + newnode->priority = priority; + newnode->next = newnode; + newnode->prev = newnode; + + heap->min = ll_merge(heap->min, newnode); + + if (priority < heap->min->priority) { + heap->min = newnode; + } + + heap->vertex_to_node[value] = newnode; + + heap->trees++; +} + +void fib_deletemin(fibh *heap) { + unsigned int min_rank = heap->min->rank; + struct fibn *min_children = heap->min->first_child; + + struct fibn *trees = ll_delroot(heap->min); + heap->vertex_to_node[heap->min->value] = NULL; + free(heap->min); + + if (min_children != NULL) + ll_setparent(min_children, NULL); + trees = ll_merge(trees, min_children); + + heap->trees += min_rank - 1; + + if (trees == NULL) { + // min had no children and was only tree, return empty heap + heap->min = NULL; + heap->rank = 0; + } else { + // min had children or there were other trees + unsigned int min_val = UINT_MAX; + struct fibn *min_ptr = NULL; + unsigned int max_rank = 0; + struct fibn *curnode = trees; + for (unsigned int i = 0; i < heap->trees; i++) { + if (curnode->priority < min_val) { + min_ptr = curnode; + min_val = curnode->priority; + } + if (curnode->rank > max_rank) + max_rank = curnode->rank; + curnode = curnode->next; + } + + if (min_ptr == NULL) + min_ptr = trees; + heap->min = min_ptr; + heap->rank = max_rank; + + struct fibn **rankslots = calloc(max_rank + 1, sizeof(struct fibn *)); + curnode = heap->min; + while (curnode != rankslots[curnode->rank]) { + if (rankslots[curnode->rank] == NULL) { + rankslots[curnode->rank] = curnode; + curnode = curnode->next; + } else { + struct fibn *oldnode = rankslots[curnode->rank]; + rankslots[curnode->rank] = NULL; + struct fibn *smaller = + curnode->priority < oldnode->priority || curnode == heap->min + ? curnode + : oldnode; + struct fibn *larger = smaller == curnode ? oldnode : curnode; + ll_delroot(larger); + ll_setparent(larger, smaller); + struct fibn *smaller_children = smaller->first_child; + smaller->first_child = ll_merge(smaller_children, larger); + heap->trees--; + smaller->rank++; + if (smaller->rank > heap->rank) { + heap->rank = smaller->rank; + rankslots = + realloc(rankslots, (heap->rank + 1) * sizeof(struct fibn *)); + rankslots[heap->rank] = NULL; + } + curnode = smaller; + } + } + free(rankslots); + } +} + +void fib_decreasekey(fibh *heap, unsigned int value, + unsigned int new_priority) { + struct fibn *node = heap->vertex_to_node[value]; + if (node != NULL) { + node->priority = new_priority; + if (node->parent != NULL) { + if (node->priority < node->parent->priority) { + + struct fibn *curnode = node; + curnode->marked = true; + while (curnode->parent != NULL && curnode->marked) { + struct fibn *oldparent = curnode->parent; + oldparent->rank--; + oldparent->first_child = ll_delroot(curnode); + ll_setparent(curnode, NULL); + ll_merge(heap->min, curnode); + heap->trees++; + if (curnode->marked) { + curnode->marked = false; + } + + if (oldparent->marked) { + curnode = oldparent; + } else { + oldparent->marked = true; + break; + } + } + } + } + + if (new_priority < heap->min->priority) { + heap->min = node; + } + } +} + +unsigned int *dijkstra(fnet *network, unsigned int source) { + unsigned int nv = network->num_dual_verts; + unsigned int *vei = network->dual_verts_to_edges_ind; + unsigned int *ve = network->dual_verts_to_edges; + unsigned int *ev = network->dual_edges_to_verts; + + unsigned int *dist = (unsigned int *)calloc(nv, sizeof(unsigned int)); + fibh *Q = (fibh *)calloc(1, sizeof(fibh)); + Q->vertex_to_node = (struct fibn **)calloc(nv, sizeof(struct fibn *)); + + for (unsigned int i = 0; i < nv; i++) { + if (i != source) { + dist[i] = UINT_MAX; + } + + fib_insert(Q, i, dist[i]); + } + + while (Q->min != NULL) { + unsigned int u = fib_findmin(Q); + fib_deletemin(Q); + + for (unsigned int i = 0; i < vei[u + 1] - vei[u]; i++) { + unsigned int e = ve[vei[u] + i]; + unsigned int v = ev[2 * e] == u ? ev[2 * e + 1] : ev[2 * e]; + unsigned int alt = dist[u] + 1; + if (alt < dist[v]) { + dist[v] = alt; + fib_decreasekey(Q, v, alt); + } + } + } + + free(Q->vertex_to_node); + free(Q); + + return dist; +} + +unsigned int **get_dists(fnet *network) { + unsigned int nv = network->num_dual_verts; + unsigned int **dists = (unsigned int **)malloc(nv * sizeof(unsigned int *)); + +#pragma omp parallel for + for (unsigned int i = 0; i < nv; i++) { + dists[i] = dijkstra(network, i); + } + + return dists; +} + +double *get_corr(finst *instance, unsigned int **dists, cholmod_common *c) { + unsigned int nv = instance->network->num_dual_verts; + unsigned int ne = instance->network->num_edges; + unsigned int *ev = instance->network->dual_edges_to_verts; + bool nulldists = false; + if (dists == NULL) { + dists = get_dists(instance->network); + nulldists = true; + } + double *corr = calloc(nv, sizeof(double)); + unsigned int *marks = get_clusters(instance, c); + unsigned int *numat = calloc(nv, sizeof(unsigned int)); + + for (unsigned int i = 0; i < ne; i++) { + unsigned int v1 = ev[2 * i]; + unsigned int v2 = ev[2 * i + 1]; + for (unsigned int j = 0; j < ne; j++) { + unsigned int v3 = ev[2 * j]; + unsigned int v4 = ev[2 * j + 1]; + unsigned int dist1 = dists[v1][v3]; + unsigned int dist2 = dists[v1][v4]; + unsigned int dist3 = dists[v2][v3]; + unsigned int dist4 = dists[v2][v4]; + unsigned int dist = (dist1 + dist2 + dist3 + dist4) / 4; + corr[dist] += instance->fuses[i] && instance->fuses[j]; + numat[dist]++; + } + } + + for (unsigned int i = 0; i < nv; i++) { + if (numat[i] > 0) { + corr[i] /= numat[i]; + } + } + + if (nulldists) { + for (int i = 0; i < nv; i++) { + free(dists[i]); + } + free(dists); + } + + free(marks); + + return corr; +} + +/* +double *get_space_corr(finst *instance, cholmod_common *c) { + unsigned int nv = instance->network->num_dual_verts; + double *vc = instance->network->dual_vert_coords; + unsigned int ne = instance->network->num_edges; + unsigned int *ev = instance->network->dual_edges_to_verts; + double *corr = calloc(nv, sizeof(unsigned int)); + unsigned int *numat = calloc(nv, sizeof(unsigned int)); + + for (unsigned int i = 0; i < ne; i++) { + unsigned int v1 = ev[2 * i]; + unsigned int v2 = ev[2 * i + 1]; + double v1x = vc[2 * v1]; + double v1y = vc[2 * v1 + 1]; + double v2x = vc[2 * v2]; + double v2y = vc[2 * v2 + 1]; + double e1x = (v1x + v2x) / 2; + double e1y = (v1y + v2y) / 2; + for (unsigned int j = 0; j < ne; j++) { + unsigned int v3 = ev[2 * j]; + unsigned int v4 = ev[2 * j + 1]; + double v1x = vc[2 * v1]; + double v1y = vc[2 * v1 + 1]; + double v2x = vc[2 * v2]; + double v2y = vc[2 * v2 + 1]; + double e1x = (v1x + v2x) / 2; + double e1y = (v1y + v2y) / 2; + corr[dist] += instance->fuses[i] && instance->fuses[j]; + numat[dist]++; + } + } + + for (unsigned int i = 0; i < nv; i++) { + if (numat[i] > 0) { + corr[i] /= numat[i]; + } + } + + for (int i = 0; i < nv; i++) { + free(dists[i]); + } + + free(marks); + free(dists); + + return corr; +} +*/ -- cgit v1.2.3-70-g09d2