compute weights in A3

Makefile

@@ -24,8 +24,8 @@ graph: graph.o weyl.o
 D2n: D2n.o weyl.o thickenings.o
 	gcc $(OPTIONS) -o D2n D2n.o weyl.o thickenings.o
 
-dominant_weights: dominant_weights.o weyl.o
+dominant_weights: dominant_weights.o weyl.o thickenings.o
-	gcc $(OPTIONS) -o dominant_weights dominant_weights.o weyl.o
+	gcc $(OPTIONS) -o dominant_weights dominant_weights.o weyl.o thickenings.o -lcdd
 
 enumerate.o: enumerate.c $(HEADERS)
 	gcc $(OPTIONS) -c enumerate.c

dominant_weights.c

@@ -1,8 +1,20 @@
 #include <stdio.h>
 #include <string.h>
 #include <malloc.h>
+#include <setoper.h>
+#include <cdd.h>
 
 #include "weyl.h"
+#include "thickenings.h"
+
+doublequotient_t *dq;
+double *vector;
+char buf[1000];
+
+typedef struct {
+	dd_MatrixPtr M;
+	dd_LPSolutionPtr lps;
+} info_t;
 
 static char* alphabetize(weylgroup_element_t *e, char *str)
 {
@@ -17,23 +29,59 @@ static char* alphabetize(weylgroup_element_t *e, char *str)
   return str;
 }
 
+static void balanced_thickening_callback(const bitvec_t *pos, int size, const enumeration_info_t *ei)
+{
+	int bit, funcbit, sign;
+	dd_rowset ImL, Lbasis;
+	dd_ErrorType err = dd_NoError;
+	info_t *info = (info_t*)ei->callback_data;
+
+	for(int i = 0; i < size; i++) {
+		bit = i < size/2 ? bv_get_bit(pos, i) : !bv_get_bit(pos, size - 1 - i);
+		sign = bit ? 1 : -1;
+		info->M->matrix[i][0][0] = 0.0;
+		info->M->matrix[i][1][0] = sign*vector[3*i];
+		info->M->matrix[i][2][0] = sign*vector[3*i+1];
+		info->M->matrix[i][3][0] = sign*vector[3*i+2];
+//		printf("0 %.2f %.2f %.2f %d %d\n", sign*vector[3*i], sign*vector[3*i+1], sign*vector[3*i+2], bit, funcbit);
+    }
+
+	dd_FindRelativeInterior(info->M, &ImL, &Lbasis, &(info->lps), &err);
+
+	if(set_card(Lbasis) != 0)
+		printf("codim = %ld\n", set_card(Lbasis));
+	else
+		printf("weight = (%.4f, %.4f, %.4f)\n", *info->lps->sol[0], *info->lps->sol[1], *info->lps->sol[2]);
+
+
+	dd_FreeLPSolution(info->lps);
+
+	printf("\n");
+}
+
 static int apply_reflection(double *in, double *out, int rank, int reflection)
 {
 	memcpy(out, in, rank*sizeof(double));
+	/*
 	out[reflection] *= -1;
 	if(reflection != 0)
 		out[reflection-1] += in[reflection];
 	if(reflection != rank-1)
 		out[reflection+1] += in[reflection];
+	*/
+	out[reflection] *= -1;
+	if(reflection != 0)
+		out[reflection] += in[reflection-1];
+	if(reflection != rank-1)
+		out[reflection] += in[reflection+1];
+
 }
 
 int main(int argc, char *argv[])
 {
 	semisimple_type_t type;
 	simple_type_t simple;
-	doublequotient_t *dq;
+	info_t info;
-	double *vector;
-	char buf[1000];
 
 	type.n = 1;
 	type.factors = &simple;
@@ -47,21 +95,29 @@ int main(int argc, char *argv[])
 	weylgroup_element_t *group = dq->group;
 	vector = malloc(weyl_order(type)*weyl_rank(type)*sizeof(double));
 
-	for(int i = 0; i < rank; i++) {
+	vector[0] = 3;
-		vector[i] = 1.0;
+	vector[1] = 4;
-	}
+	vector[2] = 3;
 
 	for(int i = 0; i < order; i++) {
-		printf("%s %.2f %.2f %.2f\n", alphabetize(&group[i], buf),
+		printf("%s %.2f %.2f %.2f\n", alphabetize(&group[i], buf), vector[3*i], vector[3*i+1], vector[3*i+2]);
-		       1.0/4.0*(3*vector[3*i] + 2*vector[3*i+1] + 1*vector[3*i+2]),
-		       1.0/4.0*(2*vector[3*i] + 4*vector[3*i+1] + 2*vector[3*i+2]),
-		       1.0/4.0*(1*vector[3*i] + 2*vector[3*i+1] + 3*vector[3*i+2]));
 
 		for(int j = 0; j < rank; j++)
 			if(group[i].left[j]->wordlength > group[i].wordlength)
 				apply_reflection(&vector[3*i], &vector[3*group[i].left[j]->index], rank, j);
 	}
 
+	printf("\n");
+
+	dd_set_global_constants();
+	info.M = dd_CreateMatrix(24, 4);
+	info.M->representation = dd_Inequality;
+	info.M->numbtype = dd_Real;
+	info.M->objective = dd_LPmax;
+
+	enumerate_balanced_thickenings(dq, balanced_thickening_callback, &info);
+
+	dd_FreeMatrix(info.M);
 	weyl_destroy_bruhat(dq);
 	free(vector);