Mod graph by invariances

thickenings.c

@@ -57,24 +57,28 @@ static int compare_wordlength(const void *a, const void *b, void *gr)
   return graph[i].wordlength - graph[j].wordlength;
 }
 
-void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists_pointer, int **words_pointer) // the edgelists_pointer and words_pointer arguments are just for freeing afterwards
+void prepare_graph(semisimple_type_t type, node_t *graph)
 {
   queue_t queue;
 
-  int rank, order;
+  edgelist_t *edgelists_lower, *edgelists_higher;
+  int rank, order, hyperplanes;
   edgelist_t *edge, *previous;
-  int edgelist_count, max_wordlength, hyperplane_count;
+  int edgelist_count, hyperplane_count;
   int current;
 
   int *graph_data;
   node_t *graph_unsorted;
-  int *wordlength_order, *reverse_wordlength_order, *seen, *words;
-  edgelist_t *edgelists;
+  int *wordlength_order, *reverse_wordlength_order, *seen;
 
   // initialize
 
   rank = coxeter_rank(type);
   order = coxeter_order(type);
+  hyperplanes = coxeter_hyperplanes(type);
+
+  edgelists_higher = graph[0].bruhat_higher;
+  edgelists_lower = &graph[0].bruhat_higher[order*hyperplanes/2];
 
   graph_data = (int*)malloc(order*rank*sizeof(int));
   graph_unsorted = (node_t*)malloc(order*sizeof(node_t));
@@ -83,13 +87,10 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists
   seen = (int*)malloc(order*sizeof(int));
 
   for(int i = 0; i < order; i++) {
-    graph_unsorted[i].left = graph[i].left;
-    graph_unsorted[i].right = graph[i].right;
-    graph_unsorted[i].word = 0;
     graph_unsorted[i].wordlength = INT_MAX;
-    graph_unsorted[i].bruhat_lower = 0;
-    graph_unsorted[i].bruhat_higher = 0;
-    graph_unsorted[i].is_hyperplane_reflection = 0;
+    graph[i].bruhat_lower = 0;
+    graph[i].bruhat_higher = 0;
+    graph[i].is_hyperplane_reflection = 0;
   }
 
   // get coxeter graph
@@ -98,7 +99,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists
 
   for(int i = 0; i < order; i++)
     for(int j = 0; j < rank; j++)
-      graph_unsorted[i].left[j] = graph_data[i*rank + j];
+      graph_unsorted[i].left = &graph_data[i*rank];
 
   // find wordlengths
 
@@ -115,11 +116,6 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists
     }
   }
 
-  max_wordlength = 0;
-  for(int i = 0; i < order; i++)
-    if(graph_unsorted[i].wordlength > max_wordlength)
-      max_wordlength = graph_unsorted[i].wordlength;
-
   // sort by wordlength
 
   for(int i = 0; i < order; i++)
@@ -128,23 +124,22 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists
   for(int i = 0; i < order; i++)
     reverse_wordlength_order[wordlength_order[i]] = i; // reverse_wordlength_order is a map old index -> new index
   for(int i = 0; i < order; i++) {
-    graph[i] = graph_unsorted[wordlength_order[i]]; // copy the whole thing
+    // we have only set left and wordlength so far, so just copy these
+    graph[i].wordlength = graph_unsorted[wordlength_order[i]].wordlength;
     for(int j = 0; j < rank; j++)
-      graph[i].left[j] = reverse_wordlength_order[graph[i].left[j]]; // rewrite references
+      graph[i].left[j] = reverse_wordlength_order[graph_unsorted[wordlength_order[i]].left[j]]; // rewrite references
   }
 
   // find words
 
-  words = (int*)malloc(order*max_wordlength*sizeof(int));
-  memset(words, 0, order*max_wordlength*sizeof(int));
-  graph[0].word = &words[0];
+  for(int i = 0; i < order; i++)
+    memset(graph[i].word, 0, hyperplanes*sizeof(int));
   queue_init(&queue);
   queue_put(&queue, 0);
   while((current = queue_get(&queue)) != -1) {
     for(int i = 0; i < rank; i++) {
       int neighbor = graph[current].left[i];
-      if(graph[neighbor].wordlength == graph[current].wordlength + 1 && graph[neighbor].word == 0) {
-	graph[neighbor].word = &words[neighbor*max_wordlength];
+      if(graph[neighbor].wordlength == graph[current].wordlength + 1 && graph[neighbor].word[0] == 0) {
 	memcpy(&graph[neighbor].word[1], &graph[current].word[0], graph[current].wordlength*sizeof(int));
 	graph[neighbor].word[0] = i;
 	queue_put(&queue, neighbor);
@@ -197,7 +192,6 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists
 
   // generate folding order
 
-  edgelists = (edgelist_t*)malloc(order*hyperplane_count*sizeof(edgelist_t));
   edgelist_count = 0;
   for(int i = 0; i < order; i++) {
     if(graph[i].is_hyperplane_reflection) {
@@ -208,9 +202,9 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists
 	  current = graph[current].left[graph[i].word[k]];
 
 	if(graph[j].wordlength < graph[current].wordlength) { // current has higher bruhat order than j
-	  edgelists[edgelist_count].to = j;
-	  edgelists[edgelist_count].next = graph[current].bruhat_lower;
-	  graph[current].bruhat_lower = &edgelists[edgelist_count];
+	  edgelists_lower[edgelist_count].to = j;
+	  edgelists_lower[edgelist_count].next = graph[current].bruhat_lower;
+	  graph[current].bruhat_lower = &edgelists_lower[edgelist_count];
 	  edgelist_count++;
 	} else if(graph[j].wordlength > graph[current].wordlength) { // j has higher bruhat order than current; these are already included from the other side
 	} else {
@@ -263,20 +257,18 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists
 
   // reverse folding order
 
+  edgelist_count = 0;
   for(int i = 0; i < order; i++) {
     edge = graph[i].bruhat_lower;
     while(edge) {
-      edgelists[edgelist_count].to = i;
-      edgelists[edgelist_count].next = graph[edge->to].bruhat_higher;
-      graph[edge->to].bruhat_higher = &edgelists[edgelist_count];
+      edgelists_higher[edgelist_count].to = i;
+      edgelists_higher[edgelist_count].next = graph[edge->to].bruhat_higher;
+      graph[edge->to].bruhat_higher = &edgelists_higher[edgelist_count];
       edgelist_count++;
       edge = edge->next;
     }
   }
 
-  *edgelists_pointer = edgelists;
-  *words_pointer = words;
-
   free(graph_data);
   free(graph_unsorted);
   free(wordlength_order);
@@ -284,6 +276,264 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists
   free(seen);
 }
 
+static int edgelist_contains(edgelist_t *list, int x) {
+  while(list) {
+    if(list->to == x)
+      return 1;
+    list = list->next;
+  }
+  return 0;
+}
+
+static edgelist_t *edgelist_add(edgelist_t *list, int new, edgelist_t *storage, int *storage_index)
+{
+  edgelist_t *new_link = &storage[*storage_index];
+  new_link->next = list;
+  new_link->to = new;
+  (*storage_index)++;
+  return new_link;
+}
+
+int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigned long right, node_t *simplified_graph)
+{
+  node_t *full_graph;
+  int edgelists_used;
+  int rank, order, hyperplanes;
+  int *reduced, *group, *simplified;
+  int *seen;
+  int current;
+  edgelist_t *edge, *previous;
+  queue_t queue;
+  int ncosets;
+
+  if(opposition_involution(type, left) != left)
+    return -1;
+
+  edgelist_t *edgelists_higher = &simplified_graph[0].bruhat_higher[0];
+  edgelist_t *edgelists_lower = &simplified_graph[0].bruhat_higher[order*hyperplanes/2];
+
+  // get full graph
+
+  full_graph = graph_alloc(type);
+  prepare_graph(type, full_graph);
+
+  // initialize stuff
+
+  rank = coxeter_rank(type);
+  order = coxeter_order(type);
+  hyperplanes = coxeter_hyperplanes(type);
+
+  reduced = (int*)malloc(order*sizeof(int));
+  group = (int*)malloc(order*sizeof(int));
+  simplified = (int*)malloc(order*sizeof(int));
+  for(int i = 0; i < order; i++) {
+    group[i] = -1;
+    reduced[i] = i;
+  }
+
+  // step 1: group
+  for(int i = 0; i < order; i++) {
+    if(group[i] != -1)
+      continue;
+
+    queue_init(&queue);
+    queue_put(&queue, i);
+    while((current = queue_get(&queue)) != -1) {
+      if(group[current] != -1)
+        continue;
+      group[current] = i;
+
+      for(int j = 0; j < rank; j++) {
+        if(left & (1 << j))
+          queue_put(&queue, full_graph[current].left[j]);
+        if(right & (1 << j))
+          queue_put(&queue, full_graph[current].right[j]);
+      }
+    }
+  }
+
+  // step 2: find minimum
+  for(int i = 0; i < order; i++)
+    if(full_graph[i].wordlength < full_graph[reduced[group[i]]].wordlength)
+      reduced[group[i]] = i;
+
+  // step 3: assign minimum to all
+  for(int i = 0; i < order; i++)
+    reduced[i] = reduced[group[i]];
+
+  // step 4: assign indices to cosets
+  ncosets = 0;
+  for(int i = 0; i < order; i++)
+    if(reduced[i] == i)
+      simplified[i] = ncosets++;
+
+  for(int i = 0; i < order; i++)
+    simplified[i] = simplified[reduced[i]];
+
+  //  fprintf(stderr, "Number of double cosets: %d\n\n", ncosets);
+
+  //  simplified_graph = (node_t*) malloc(ncosets*sizeof(node_t));
+  seen = (int*) malloc(ncosets*sizeof(int));
+  edgelists_used = 0;
+
+  // step 5: set up nodes from minima
+  current = 0;
+  for(int i = 0; i < order; i++)
+    if(reduced[i] == i) { // is minimum
+      memcpy(simplified_graph[simplified[i]].word, full_graph[i].word, full_graph[i].wordlength*sizeof(int));
+      simplified_graph[simplified[i]].wordlength = full_graph[i].wordlength;
+      simplified_graph[simplified[i]].opposite = simplified[full_graph[i].opposite];
+      simplified_graph[simplified[i]].bruhat_lower = (edgelist_t*)0;
+      simplified_graph[simplified[i]].bruhat_higher = (edgelist_t*)0;
+      for(int j = 0; j < rank; j++) {
+	simplified_graph[simplified[i]].left[j] = simplified[full_graph[i].left[j]];
+	simplified_graph[simplified[i]].right[j] = simplified[full_graph[i].right[j]];
+      }
+    }
+
+  // step 6: find order relations
+  for(int i = 0; i < order; i++) {
+    edge = full_graph[i].bruhat_lower;
+    while(edge) {
+      int this = simplified[i];
+      int that = simplified[edge->to];
+      if(this != that) {
+	// found something
+	if(!edgelist_contains(simplified_graph[this].bruhat_lower, that))
+	  simplified_graph[this].bruhat_lower = edgelist_add(simplified_graph[this].bruhat_lower, that, edgelists_lower, &edgelists_used);
+	ERROR(simplified_graph[this].wordlength <= simplified_graph[that].wordlength, "The order assumption is being violated!\n");
+      }
+      edge = edge->next;
+    }
+  }
+
+  // step 7: remove redundant edges
+  for(int i = 0; i < ncosets; i++) {
+    memset(seen, 0, ncosets*sizeof(int));
+    queue_init(&queue);
+
+    for(int len = 1; len <= simplified_graph[i].wordlength; len++) {
+      edge = simplified_graph[i].bruhat_lower;
+      previous = (edgelist_t*)0;
+
+      while(edge) {
+	// only look at edges of this length now
+	if(simplified_graph[i].wordlength - simplified_graph[edge->to].wordlength != len) {
+	  // we only consider edges of length len in this pass
+	  previous = edge;
+	} else if(seen[edge->to]) {
+	  // this edge is redundant, remove it
+	  //	  fprintf(stderr, "removing edge from %d to %d\n", i, edge->to);
+	  if(previous)
+	    previous->next = edge->next;
+	  else
+	    simplified_graph[i].bruhat_lower = edge->next;
+	} else {
+	  // this edge was not redundant, add to seen
+	  previous = edge;
+	  seen[edge->to] = 1;
+	  queue_put(&queue, edge->to);
+	}
+	edge = edge->next;
+      }
+
+      // calculate transitive closure of seen nodes
+      while((current = queue_get(&queue)) != -1) {
+	edge = simplified_graph[current].bruhat_lower;
+	while(edge) {
+	  if(!seen[edge->to]) {
+	    seen[edge->to] = 1;
+	    queue_put(&queue, edge->to);
+	  }
+	  edge = edge->next;
+	}
+      }
+    }
+  }
+
+  // step 8: revert order
+  for(int i = 0; i < ncosets; i++) {
+    edge = simplified_graph[i].bruhat_lower;
+    while(edge) {
+      simplified_graph[edge->to].bruhat_higher =
+	edgelist_add(simplified_graph[edge->to].bruhat_higher,
+		     i, edgelists_higher, &edgelists_used);
+      edge = edge->next;
+    }
+  }
+
+  // output as graphviz dot file
+  /*
+  fprintf(stdout, "difull_graph test123 {\n");
+  for(int i = 0; i < ncosets; i++) {
+    edge = simplified_graph[i].bruhat_lower;
+    while(edge) {
+      fprintf(stdout, "%s -> %s;\n",
+      	     alphabetize(simplified_graph[i].word, simplified_graph[i].wordlength, alphabet, buffer),
+	     alphabetize(simplified_graph[edge->to].word, simplified_graph[edge->to].wordlength, alphabet, buffer2));
+
+      edge = edge->next;
+    }
+  }
+  fprintf(stdout, "}\n"); */
+
+  // some output
+  /*  for(int i = 0; i < ncosets; i++)
+      fprintf(stderr, "%s <=> %s\n", simplified_graph[i].wordlength == 0 ? "1" : alphabetize(simplified_graph[i].word, simplified_graph[i].wordlength, alphabet, buffer), simplified_graph[simplified_graph[i].opposite].wordlength == 0 ? "1" : alphabetize(simplified_graph[simplified_graph[i].opposite].word, simplified_graph[simplified_graph[i].opposite].wordlength, alphabet, buffer2)); */
+
+  //  fprintf(stderr, "\nAdded %d edges.\n\n", edgelists_used);
+
+  free(seen);
+  free(reduced);
+  free(group);
+  free(simplified);
+  graph_free(type, full_graph);
+
+  return ncosets;
+}
+
+node_t *graph_alloc(semisimple_type_t type)
+{
+  int rank = coxeter_rank(type);
+  int order = coxeter_order(type);
+  int hyperplanes = coxeter_hyperplanes(type);
+
+  node_t *graph = (node_t*)malloc(order*sizeof(node_t));
+  int *left = (int*)malloc(order*rank*sizeof(int));
+  int *right = (int*)malloc(order*rank*sizeof(int));
+  edgelist_t *edgelists = (edgelist_t*)malloc(order*hyperplanes*sizeof(edgelist_t));
+  int *words = (int*)malloc(order*hyperplanes*sizeof(int));
+
+  for(int i = 0; i < order; i++) {
+    graph[i].left = &left[rank*i];
+    graph[i].right = &right[rank*i];
+    graph[i].word = &words[hyperplanes*i];
+  }
+
+  graph[0].bruhat_higher = edgelists;
+
+  return graph;
+}
+
+void graph_free(semisimple_type_t type, node_t *graph)
+{
+  free(graph[0].left);
+  free(graph[0].right);
+  free(graph[0].word);
+
+  int order = coxeter_order(type);
+
+  // find the head of all edgelists by just taking the one having the lowest address
+  edgelist_t *edgelists = graph[0].bruhat_lower;
+  for(int i = 0; i < order; i++) {
+    if(graph[i].bruhat_lower < edgelists && graph[i].bruhat_lower != 0)
+      edgelists = graph[i].bruhat_lower;
+    if(graph[i].bruhat_higher < edgelists && graph[i].bruhat_higher != 0)
+      edgelists = graph[i].bruhat_higher;
+  }
+  free(edgelists);
+}
+
 /*********************************** THE ACTUAL ENUMERATION ****************************************/
 
 typedef struct {
@@ -401,7 +651,6 @@ static long enumerate_tree(const enumeration_info_t *info, signed char *level, i
 
 long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, int size, const char *alphabet, FILE *outfile)
 {
-  //  int rank, order;
   signed char *level;
   long count = 0;
   enumeration_info_t info;