simplified graphs

Makefile

@@ -1,5 +1,5 @@
 HEADERS=coxeter.h thickenings.h queue.h
-OPTIONS=-O3 -m64 -march=native -flto -funroll-loops -std=gnu99 -pg
+OPTIONS=-O3 -m64 -march=native -flto -funroll-loops -std=gnu99
 
 all: generate process
 

generate.c

@@ -72,7 +72,7 @@ int main(int argc, const char *argv[])
 
   fwrite(&type.n, sizeof(int), 1, stdout);
   fwrite(type.factors, sizeof(simple_type_t), type.n, stdout);
-  long count = enumerate_balanced_thickenings(type, graph, alphabet, stdout);
+  long count = enumerate_balanced_thickenings(type, graph, order, alphabet, stdout);
 
   fprintf(stderr, "\n");
   fprintf(stderr, "Found %ld balanced thickenings\n\n", count);

test.c

@@ -0,0 +1,257 @@
+#include <stdio.h>
+#include <memory.h>
+
+#include "thickenings.h"
+#include "queue.h"
+
+#define SWAP(t, a, b) do {t tmp = a; a = b; b = tmp;} while(0)
+
+int edgelist_contains(edgelist_t *list, int needle) {
+  while(list) {
+    if(list->to == needle)
+      return 1;
+    list = list->next;
+  }
+  return 0;
+}
+
+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 main(int argc, const char *argv[])
+{
+  semisimple_type_t type;
+  node_t *graph;
+  int *leftbuf, *rightbuf;
+  edgelist_t *edgelists;
+  edgelist_t *edgelists_simplified;
+  int edgelists_simplified_used;
+  int *words;
+  int rank, order, max_wordlength;
+  int *reduced, *group, *simplified;
+  int *seen;
+  int current;
+  edgelist_t *edge, *previous;
+  queue_t queue;
+
+  char alphabet[] = "abcdefghijklmnopqrstuvwxyz";
+  char buffer[1024], buffer2[1024];
+
+  int ncosets;
+  node_t *simplified_graph;
+
+  // left and right invariances as bitmasks
+  //  int left = ~(1 << (atoi(argv[1]) - atoi(argv[3])));
+  //  int right = ~(1 << (atoi(argv[1]) - atoi(argv[2])));
+  int left = atoi(argv[2]);
+  int right = atoi(argv[3]);
+
+  type.n = 1;
+  type.factors = (simple_type_t*)malloc(type.n*sizeof(simple_type_t));
+  type.factors[0].series = 'B';
+  type.factors[0].rank = atoi(argv[1]);
+
+  rank = coxeter_rank(type);
+  order = coxeter_order(type);
+  graph = (node_t*)malloc(order*sizeof(node_t));
+  leftbuf = (int*)malloc(rank*order*sizeof(int));
+  rightbuf = (int*)malloc(rank*order*sizeof(int));
+  for(int i = 0; i < order; i++) {
+    graph[i].left = &leftbuf[i*rank];
+    graph[i].right = &rightbuf[i*rank];
+  }
+
+  prepare_graph(type, graph, &edgelists, &words);
+
+  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, graph[current].left[j]);
+        if(right & (1 << j))
+          queue_put(&queue, graph[current].right[j]);
+      }
+    }
+  }
+
+  // step 2: find minimum
+  for(int i = 0; i < order; i++)
+    if(graph[i].wordlength < 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]];
+
+  // count 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);
+
+  max_wordlength = coxeter_hyperplanes(type);
+  simplified_graph = (node_t*) malloc(ncosets*sizeof(node_t));
+  edgelists_simplified = (edgelist_t*) malloc(2*max_wordlength*order*sizeof(edgelist_t));
+  seen = (int*) malloc(ncosets*sizeof(int));
+  edgelists_simplified_used = 0;
+
+  // copy minima
+  current = 0;
+  for(int i = 0; i < order; i++)
+    if(reduced[i] == i) { // is minimum
+      simplified_graph[simplified[i]].word = graph[i].word;
+      simplified_graph[simplified[i]].wordlength = graph[i].wordlength;
+      simplified_graph[simplified[i]].opposite = simplified[graph[i].opposite];
+      simplified_graph[simplified[i]].bruhat_lower = (edgelist_t*)0;
+      simplified_graph[simplified[i]].bruhat_higher = (edgelist_t*)0;
+    }
+
+  // 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));
+
+  // find order relations
+  for(int i = 0; i < order; i++) {
+    edge = 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_simplified, &edgelists_simplified_used);
+	//	if(!edgelist_contains(simplified_graph[that].bruhat_higher, this))
+	//	  simplified_graph[that].bruhat_higher = edgelist_add(simplified_graph[that].bruhat_higher, this, edgelists_simplified, &edgelists_simplified_used);
+	ERROR(simplified_graph[this].wordlength <= simplified_graph[that].wordlength, "The order assumption is being violated!\n");
+      }
+      edge = edge->next;
+    }
+  }
+
+  fprintf(stderr, "\nAdded %d edges.\n\n", edgelists_simplified_used);
+
+  // 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;
+	}
+      }
+    }
+  }
+
+  // reverse 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_simplified, &edgelists_simplified_used);
+      edge = edge->next;
+    }
+  }
+
+
+  fprintf(stdout, "digraph 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");
+
+  long nthickenings = enumerate_balanced_thickenings(type, simplified_graph, ncosets, alphabet, stdout);
+
+  fprintf(stderr, "Found %ld balanced thickenings.\n", nthickenings);
+
+  /*  for(int i = 0; i < order; i++)
+    printf("%s <= %s\n",
+	   reduced[i] == 0 ? "1" : alphabetize(graph[reduced[i]].word, graph[reduced[i]].wordlength, alphabet, buffer2),
+	   i == 0 ? "1" : alphabetize(graph[i].word, graph[i].wordlength, alphabet, buffer)); */
+
+
+  free(seen);
+  free(simplified_graph);
+  free(edgelists_simplified);
+  free(type.factors);
+  free(graph);
+  free(edgelists);
+  free(words);
+  free(reduced);
+  free(group);
+  free(simplified);
+  free(leftbuf);
+  free(rightbuf);
+
+  return 0;
+}

thickenings.c

@@ -289,6 +289,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists
 typedef struct {
   int rank;
   int order;
+  int size;        // the size of the graph; this can vary from the order if we take quotients beforehand
   const node_t *graph;
   int printstep;
   const char *alphabet;
@@ -307,7 +308,7 @@ static int transitive_closure(const enumeration_info_t *info, signed char *level
   level[info->graph[head].opposite] = -current_level;
   queue_put(&queue, head);
 
-  for(int i = head + 1; level[i] != HEAD_MARKER && i < info->order; i++) { // everything which is right to the head and empty will not get marked in this or higher levels, so we can mark its opposite
+  for(int i = head + 1; level[i] != HEAD_MARKER && i < info->size; i++) { // everything which is right to the head and empty will not get marked in this or higher levels, so we can mark its opposite
     if(level[i] == current_level) {
       is_slim = 0;
       break;
@@ -342,11 +343,11 @@ static inline void output_thickening(const enumeration_info_t *info, signed char
 {
   // if printstep is set accordingly, write state to stderr
   if(is_slim && is_fat && info->printstep > 0 && (count + 1) % info->printstep == 0) {
-    print_thickening(info->rank, info->order, level, current_level, info->alphabet, stderr);
+    print_thickening(info->rank, info->size, level, current_level, info->alphabet, stderr);
     fprintf(stderr, "\n");
   }
   else if(info->printstep < 0) {
-    print_thickening(info->rank, info->order, level, current_level - !is_slim, info->alphabet, stderr);
+    print_thickening(info->rank, info->size, level, current_level - !is_slim, info->alphabet, stderr);
     fprintf(stderr, " ");
     if(is_slim) {
       fprintf(stderr, "S");
@@ -378,7 +379,7 @@ static long enumerate_tree(const enumeration_info_t *info, signed char *level, i
 
     if(is_fat) {
       count++;
-      fwrite(level, sizeof(signed char), info->order, info->outfile);
+      fwrite(level, sizeof(signed char), info->size, info->outfile);
     } else {
       for(int i = head - 1; i >= 0; i--)
 	if(level[i] == 0)
@@ -388,14 +389,14 @@ static long enumerate_tree(const enumeration_info_t *info, signed char *level, i
 
   // clean up
   level[head] = 0;
-  for(int i = 0; i < info->order; i++)
+  for(int i = 0; i < info->size; i++)
     if(level[i] >= current_level && level[i] != HEAD_MARKER || level[i] <= -current_level)
       level[i] = 0;
 
   return count;
 }
 
-long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, const char *alphabet, FILE *outfile)
+long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, int size, const char *alphabet, FILE *outfile)
 {
   //  int rank, order;
   signed char *level;
@@ -406,6 +407,7 @@ long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, const
 
   info.rank = coxeter_rank(type);
   info.order = coxeter_order(type);
+  info.size = size;
   info.graph = graph;
   info.alphabet = (char*)alphabet;
   info.outfile = outfile;
@@ -414,10 +416,16 @@ long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, const
   if(getenv("PRINTSTEP"))
     info.printstep = atoi(getenv("PRINTSTEP"));
 
-  level = (signed char*)malloc(info.order*sizeof(int));
+  // the algorithm only works if the opposition pairing does not stabilize any element
-  memset(level, 0, info.order*sizeof(int));
+  // if this happens, there can be no balanced thickenings
+  for(int i = 0; i < info.size; i++)
+    if(graph[i].opposite == i)
+      return 0;
 
-  for(int i = info.order - 1; i >= 0; i--)
+  level = (signed char*)malloc(info.size*sizeof(int));
+  memset(level, 0, info.size*sizeof(int));
+
+  for(int i = info.size - 1; i >= 0; i--)
     count += enumerate_tree(&info, level, 1, i);
 
   free(level);

thickenings.h

@@ -29,6 +29,6 @@ char *alphabetize(int *word, int len, const char *alphabet, char *buffer);
 void print_thickening(int rank, int order, const signed char *thickening, int level, const char *alphabet, FILE *f);
 static int compare_wordlength(const void *a, const void *b, void *gr);
 void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists_pointer, int **words_pointer);
-long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, const char *alphabet, FILE *outfile);
+long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, int size, const char *alphabet, FILE *outfile);
 
 #endif