#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], buf2[1000];

typedef struct {
	dd_MatrixPtr M;
	dd_LPSolutionPtr lps;
} info_t;

static char* alphabetize(weylgroup_element_t *e, char *str)
{
  if(e->wordlength == 0)
    sprintf(str, "1");
  else {
    for(int j = 0; j < e->wordlength; j++)
      str[j] = e->word[j] + 'a';
    str[e->wordlength] = 0;
  }

  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;
		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 = (%s)\n", print_vector_ptr(info->lps->sol, rank + 1, buf));

	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;
	info_t info;

	type.n = 1;
	type.factors = &simple;
	simple.series = 'A';
	simple.rank = 3;

	dq = weyl_generate_bruhat(type, 0, 0);

	int order = weyl_order(type);
	int rank = weyl_rank(type);
	weylgroup_element_t *group = dq->group;
	vector = malloc(weyl_order(type)*weyl_rank(type)*sizeof(double));

	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 (%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[rank*i], &vector[rank*group[i].left[j]->index], rank, j);
	}

	printf("\n");

	dd_set_global_constants();
	info.M = dd_CreateMatrix(order+rank, rank+1);
	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);

	return 0;
}