interface(echo = 0, quiet = true, screenwidth = infinity):

# The purpose of this script is to take each triple of weights from a
# list of such triples and to produce the LattE input file for the type
# D_r Berenstein-Zelevinsky polytope corresponding to that triple of
# weights.  All the weights in the list must have the same rank.  Recall
# that, given the tiple of weights (la, mu, nu), the corresponding
# BZ-polytope contains a number integer lattice points equal to the
# Clebsch--Gordan coefficient C_{la, nu}^mu (note the position of each
# weight).
#
# To use this program, create a file called "weights" containing triples
# of weights, each preceded by a NULL line.  The file "weights" should
# be in the same directory as this Maple script.  This program assumes
# that your weights are written in the basis of fundamental weights for
# D_r and that they all have the same rank.
#
# Here is an example of the contents of an input file "weights"
# containing two triples of weights, each with rank 4:
#
# NULL
# [17, 21, 20, 18]
# [34, 18, 33, 10]
# [20, 34, 18, 12]
# NULL
# [36, 9, 13, 37]
# [26, 32, 37, 32]
# [28, 21, 19, 0]
#
# Make sure (1) to place a NULL line before every triple -- even the
# first one, and (2) *not* to put a new line (LF) at the end of the last
# triple in your file.
#
# After you have created the file "weights", create a folder called
# "LattEPolytopes" in the same directory as this Maple file.
#
# You are now ready to run the maple script by typing
# 
#    maple<Dr_BZpolytopesForLattE
#    
# in a unix terminal.
#
# This will produce all of the LattE input files for the triples of
# weights in the file "weights".  The LattE files will appear in the
# directory "LattEPolytopes".  Each file's name will express the triple
# of weights which produce the polytope described by that file.
#
# The code for the script follows.


with(linalg):
with(LinearAlgebra):
with(StringTools):

# Define r to be the rank of the weights from the file "weights".
#
readline(weights):

r := nops(parse(readline(weights))):

close(weights):

# Generate the simple roots alpha_{-1}, alpha_1, .  . .  , alpha_{r-1}
# for D_r (with respect to the basis of fundamental weights).
#
al[-1] := [ 2, 0, -1, seq(0, n=4..r) ]:

al[1] := [ 0, 2, -1, seq(0, n=4..r) ]:

al[2] := [ -1, -1, 2, -1, seq(0, n=5..r) ]:

for i from 3 to r-2 do
	al[i] := [ seq(0, n=1..i-3), -1, 2, -1, seq(0, n=i+1..r) ]:
od:

al[r-1] := [ seq(0, n=1..r-2), -1, 2 ]:

# Generate the simple coroots alpha^\vee_{-1}, alpha^\vee_1, .  . .  ,
# alpha^\vee_{r-1} corresponding to the simple roots generated above.
#
	
alcheck[-1] := [ 1, seq(0, n=1..r-1) ]:

for i from 1 to r-1 do
	alcheck[i] := [ seq(0, n=1..i), 1, seq(0, n=i+2..r) ]:
od:

# So long as there are more weights to be read from the file "weights",
# keep generating LattE files.
#
while readline(weights) <> 0 do

	# Set the interface to the ideal settings for writing to BZfromMaple
	#
	interface(prettyprint = false):
	
	# Define la, mu, and nu to be the current weights we will use to
	# generate the BZ polytope.
	#
	la := parse(readline(weights)):
	mu := parse(readline(weights)):
	nu := parse(readline(weights)):

	# Unassign the values of the t[i,j] assigned the last time through
	# the while loop.
	#
	unassign('t');
	
	# Set t[i,j] = 0 unless 1 <= |i| <= j < r.
	#
	for j from -2*r to 2*r do
		for i from -2*r to 2*r do
			if not (1 <= abs(i) and abs(i) <= j and j < r) then
				t[i,j] := 0;
			fi;
		od;
    od:
    
    # We now generate the equalities and inequalities.  These equalities
    # and inequalities are divided into four types.  For Types (1), (3),
    # and (4), the equal signs in the equations will all be interpreted
    # as inequalites <= by LattE, but they are treated as equalities for
    # now, since this makes them easier to translate into matrices in
    # Maple.
	
	# Type (1) inequalities
	#
	# First the j=1 case.
	#
	type1 := [ -t[-1,1] = 0, -t[1,1] = 0 ]:
	
	# And now the 2 <= j < r cases.
	#
	for j from 2 to r-1 do
		for i from -j to -3 do
			type1 := [ op(type1), -t[i,j] + t[i+1,j] = 0 ]:
		od:
		
		type1 := [ op(type1), -t[-2,j] + t[-1,j] = 0 ]:
		type1 := [ op(type1), -t[-2,j] + t[1,j] = 0 ]:
		
		type1 := [ op(type1), -t[-1,j] + t[2,j] = 0 ]:
		type1 := [ op(type1), -t[1,j] + t[2,j] = 0 ]:
		
		for i from 2 to j-1 do
			type1 := [ op(type1), -t[i,j] + t[i+1,j] = 0 ]:
		od:
		
		type1 := [op(type1), -t[j,j] = 0 ]:
	od:
	
	# Type (2) equalities.  (These will be interpreted as actual
	# equalities by LattE.)
	#
	type2 := []:
	
	for k from 1 to r do
		kthlinesum1 := 0:
		
		for j from 1 to r-1 do
			kthlinesum1 := kthlinesum1 + t[-1,j] * al[-1][k] + t[1,j] * al[1][k]:
		od:
		
		kthlinesum2 := 0:
		
		for j from 2 to r-1 do
			for i from -j to -2 do
				kthlinesum2 := kthlinesum2 + t[i,j] * al[-i][k]:
			od:
			
			for i from 2 to j do
				kthlinesum2 := kthlinesum2 + t[i,j] * al[i][k]:
			od:
		od:
		
		type2 := [ op(type2), kthlinesum1 + kthlinesum2 = la[k] + nu[k] - mu[k] ]:
	od:
	
	# Type (3) inequalities.
	#
	# But first, convert the weights and coroots to capital-V Vectors so
	# that we can do easy dot products in LinearAlgebra.
	#
	la := convert(la, Vector):
	mu := convert(mu, Vector):
	nu := convert(nu, Vector):
	
	alcheck[-1] := convert(alcheck[-1], Vector):
	
	for j from 1 to r-1 do
		alcheck[j] := convert(alcheck[j], Vector):
	od:
	
	type3 := [ t[-1,1] = la.alcheck[-1], t[1,1] = la.alcheck[1] ]:
	
	for j from 2 to r-1 do
		type3 := [ op(type3), t[j,j] = la.alcheck[j] ]:
		type3 := [ op(type3), t[-1,j] - t[1,j-1] = la.alcheck[j] ]:
		type3 := [ op(type3), t[1,j] - t[-1,j-1] = la.alcheck[j] ]:
		type3 := [ op(type3), t[1,j] + t[-1,j] - t[-2,j] - t[2,j-1] = la.alcheck[j] ]:
		type3 := [ op(type3), t[-2,j] + t[2,j-1] - t[1,j-1] - t[-1,j-1] = la.alcheck[j] ]:
		
		for i from 2 to j-1 do
			type3 := [ op(type3), t[i,j] + t[-i,j] - t[i+1,j-1] - t[-i-1,j] = la.alcheck[j] ]:
			type3 := [ op(type3), t[i+1,j-1] + t[-i-1,j] - t[-i,j-1] - t[i,j-1] = la.alcheck[j] ]:
			type3 := [ op(type3), t[-i,j] - t[-i,j-1] = la.alcheck[j] ]:
			type3 := [ op(type3), t[i,j] - t[i,j-1] = la.alcheck[j] ]:
		od:
	od:
	
	# Type (4) inequalities.
	#
	type4 := []:
	
	for j from 1 to r-1 do
		sum1 := 0:
		
		for k from j+1 to r do
			sum1 := sum1 + 2*t[-1,k] - t[-2,k] - t[2,k-1]:
		od:
		
		type4 := [ op(type4), t[-1,j] + sum1 = nu.alcheck[-1] ]:
	od:
	
	for j from 1 to r-1 do
		sum2 := 0:
		
		for k from j+1 to r do
			sum2 := sum2 + 2*t[1,k] - t[-2,k] - t[2,k-1]:
		od:
		
		type4 := [ op(type4), t[1,j] + sum2 = nu.alcheck[1] ]:
	od:
		
	for i from 2 to r-1 do
		for j from i to r-1 do
			sum3 := 0:
			
			for k from j+1 to r do
				sum3 := sum3 + 2*t[-i,k] + 2*t[i,k-1] - t[-i+1,k-1] - t[i-1,k-1] - t[-i-1,k] - t[i+1,k-1]:
			od:
			
			type4 := [ op(type4), t[-i,j] + sum3 = nu.alcheck[i] ]:
			
			sum4 := 0:
			
			for k from j+1 to r do
				sum4 := sum4 + 2*t[-i,k] + 2*t[i,k] - t[-i+1,k] - t[i-1,k] - t[-i-1,k] - t[i+1,k-1]:
			od:
			
			type4 := [ op(type4), t[i,j] + sum4 = nu.alcheck[i] ]:
		od:
	od:
	
	# Produce the matrix A and b such that Ax = b corresponds to the
	# linear equalities produced above.
	#
	AllIneq := [ op(type2), op(type1), op(type3), op(type4) ]:
	
	variables := []:
	
	for j from 1 to r-1 do
		for i from -j to j do
			if i <> 0 then
				variables := [op(variables), t[i,j] ];
			fi:
		od:
	od:
	
	A := genmatrix(AllIneq, variables, b):
	
	# Produce the array (b, -A), denoted bnegA.
	#
	bnegA := augment(b,-A):
	
	bnegA := convert(bnegA, Matrix):
	
	# Create a Latte-readable file (modulo square brackets, commas, and
	# other such junk) for this system of inequalities
	#
	writeto(BZfromMaple);
	print(Dimension(bnegA));
	print(convert(bnegA, listlist));
	print([linearity, nops(type2), seq(i, i = 1..nops(type2))]);
	writeto(terminal):
	close(BZfromMaple):
	
	# Now make that file really LattE-readable, and print the results to
	# a file whose file-name expresses la, mu, and nu.
	#
	interface(indentamount = 0 , prettyprint = 1):
	
	s := readline(BZfromMaple):
	news := SubstituteAll(s, ",", " "):
	
	t := readline(BZfromMaple):
	newt := SubstituteAll(t, "\[\[", ""):
	newt := SubstituteAll(newt, "\], \[", "\n"):
	newt := SubstituteAll(newt, ",", " "):
	newt := SubstituteAll(newt, "\]\]", ""):
	
	u := readline(BZfromMaple):
	newu := SubstituteAll(u, "\[", ""):
	newu := SubstituteAll(newu, ",", " "):
	newu := SubstituteAll(newu, "\]", ""):
	
	new := cat(news, "\n", newt, "\n", newu):
	
	# Generate the filename for the LattE file that will contain this
	# BZ-polytope.
	#
	filename := "LattEPolytopes/Dr-":
	
	for i from 1 to r-1 do
		filename := cat(filename, la[i], ".");
	od:
	
	filename := cat(filename, la[r], "_"):
	
	for i from 1 to r-1 do
		filename := cat(filename, mu[i], ".");
	od:
	
	filename := cat(filename, mu[r], "_"):
	
	for i from 1 to r-1 do
		filename := cat(filename, nu[i], ".");
	od:
	
	filename := cat(filename, nu[r]):
	
	# Write the BZ-polytope data in LattE format to the file filename.
	#
	writeto(filename):
	print(convert(new, symbol));
	writeto(terminal):
	close(filename):
od:

close(weights):

# Remove the file BZfromMaple that was created during this script's
# execution.
# 
! rm BZfromMaple