From 6c938f751dd5c289781d8ead98a2ee9fd71c2f2e Mon Sep 17 00:00:00 2001
From: Florian Stecker <mail@florianstecker.de>
Date: Fri, 18 Jan 2019 11:57:36 +0100
Subject: [PATCH] use cddlib to try to find a dominant weight for every
 balanced ideal

---
 dominant_weights.c | 81 ++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 68 insertions(+), 13 deletions(-)

diff --git a/dominant_weights.c b/dominant_weights.c
index 427f9bc..93f1c40 100644
--- a/dominant_weights.c
+++ b/dominant_weights.c
@@ -9,7 +9,7 @@
 
 doublequotient_t *dq;
 double *vector;
-char buf[1000];
+char buf[1000], buf2[1000];
 
 typedef struct {
 	dd_MatrixPtr M;
@@ -29,34 +29,75 @@ static char* alphabetize(weylgroup_element_t *e, char *str)
   return str;
 }
 
+static char* print_vector(double *v, int rank, char *buf)
+{
+	int written = 0;
+
+	for(int i = 0; i < rank; i++) {
+		written += sprintf(buf+written, "%.4f", v[i]);
+		if(i != rank -1) {
+			sprintf(buf+written, ", ");
+			written += 2;
+		}
+	}
+
+	return buf;
+}
+
+static char* print_vector_ptr(dd_Arow v, int rank, char *buf)
+{
+	int written = 0;
+
+	for(int i = 0; i < rank; i++) {
+		written += sprintf(buf+written, "%.4f", v[i][0]);
+		if(i != rank -1) {
+			sprintf(buf+written, ", ");
+			written += 2;
+		}
+	}
+
+	return buf;
+}
+
 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;
+	int rank = weyl_rank(dq->type);
 
 	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];
+		for(int j = 0; j < rank; j++)
+			info->M->matrix[i][j+1][0] = sign*vector[rank*i+j];
 //		printf("0 %.2f %.2f %.2f %d %d\n", sign*vector[3*i], sign*vector[3*i+1], sign*vector[3*i+2], bit, funcbit);
     }
 
+	for(int i = 0; i < rank; i++) {
+		info->M->matrix[i+size][0][0] = 0.0;
+		for(int j = 0; j < rank; j++) {
+			if(i == j)
+				info->M->matrix[i+size][j+1][0] = 1.0;
+			else
+				info->M->matrix[i+size][j+1][0] = 0.0;
+		}
+	}
+
+//	dd_WriteMatrix(stdout, info->M);
+
 	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]);
-
+		printf("weight = (%s)\n", print_vector_ptr(info->lps->sol, rank + 1, buf));
 
 	dd_FreeLPSolution(info->lps);
 
-	printf("\n");
+//	printf("\n");
 }
 
 static int apply_reflection(double *in, double *out, int rank, int reflection)
@@ -95,22 +136,36 @@ int main(int argc, char *argv[])
 	weylgroup_element_t *group = dq->group;
 	vector = malloc(weyl_order(type)*weyl_rank(type)*sizeof(double));
 
-	vector[0] = 3;
-	vector[1] = 4;
-	vector[2] = 3;
+	for(int i = 0; i < rank; i++)
+		vector[i] = 0.0;
+
+	for(int i = 0; 2*i < rank; i++) {
+		for(int j = i; j < rank - i; j++) {
+			vector[j] += rank - 2*i;
+		}
+	}
+
+	/*
+	vector[0] = 4.1;
+	vector[1] = 6.1;
+	vector[2] = 5.9;
+	vector[3] = 3.9;
+	*/
+
+	printf("regular element: (%s)\n", print_vector(vector, rank, buf));
 
 	for(int i = 0; i < order; i++) {
-		printf("%s %.2f %.2f %.2f\n", alphabetize(&group[i], buf), vector[3*i], vector[3*i+1], vector[3*i+2]);
+		printf("%s (%s)\n", alphabetize(&group[i], buf), print_vector(vector + rank*i, rank, buf2));
 
 		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);
+				apply_reflection(&vector[rank*i], &vector[rank*group[i].left[j]->index], rank, j);
 	}
 
 	printf("\n");
 
 	dd_set_global_constants();
-	info.M = dd_CreateMatrix(24, 4);
+	info.M = dd_CreateMatrix(order+rank, rank+1);
 	info.M->representation = dd_Inequality;
 	info.M->numbtype = dd_Real;
 	info.M->objective = dd_LPmax;