#!/usr/bin/perl
#Kazhdan-Lusztig matrix
#kl -h for a help file
#jda@math.umd.edu
use strict;
use Data::Dumper;

use Getopt::Long;
use vars qw($help $outputFile $w $quiet $standard $irreducible $q $transpose $format $inverse $parameter $blockFile $file $signed);

my %options=('h' => \$help,  
	     'o' => \$outputFile,
	     'q' => \$q,
	     't' => \$transpose,
	     'i' => \$inverse,
	     's' => \$signed,
	     'f' => \$format,
	     'b' => \$blockFile,
	     'p' => \$parameter,
	     'quiet' => \$quiet);

GetOptions(\%options,qw(h! o=s p=s w=s q! s! i! quiet! t! f=s d=s i! p=s b=s));

if ($outputFile){
    $quiet or print "Sending output to file $outputFile\n";
    open(OUT,">$outputFile")||die("Can't open $outputFile for output");
    select(OUT);
}

my $pml_tolerance          = 1e-12;
$file=$ARGV[0];
$help and &help;
unless ($file){
    print "You must specify an input file (output of klbasis). See -h for help\n";
    exit;
}

my ($P,$Lengths)=process($file);
#my $lengths=getLengths($blockFile) if ($blockFile);
my $Q;

output($P);
$outputFile and close(OUT);
sub getLengths{
    my $file=shift;
    my %lengths;
    open(IN,"<$file")||die("Can't open $file for input\n");
    foreach my $line (<IN>){
	my ($parameter,$length)=($line =~ /\s*([0-9]*).*\]\s*([0-9*]).*/);
	$lengths{$parameter}=$length;
    }
    return \%lengths;
}



sub process{
    my $file=shift;
    open(IN,"<$file")||die("Can't open $file for input");
    my $n=0;
    foreach my $line (<IN>){
	chomp($line);
	$line =~ s/\s//g;
	my @line=split ':', $line;
	next unless (scalar(@line)==3);
	$n = $line[0] if ($n<$line[0]);
    }
    close(IN);
    my @P;
    my %Lengths;
    foreach my $j (0..$n){
	    push @P, [(0)x($n+1)];
	}
    open(IN,"<$file")||die("Can't open $file for input");
    my $i=0;

    foreach my $line (<IN>){
	chomp($line);
	$line =~ s/\#.*//;
	$line =~ s/\s//g;
	next unless $line;
	next unless $line =~ /[0-9]+:[0-9q]+/;

	my @line=split ':', $line;
	if (scalar(@line) == 2){
	    unshift @line,$i;
	}else{
	    $i=$line[0];
	    $Lengths{$i}=0;
	}
	if ($line[1]<$line[0]){
	    if ($Lengths{$line[1]}+1 > $Lengths{$line[0]}){
		$Lengths{$line[0]} = $Lengths{$line[1]}+1;
		}
	}
	#print Dumper(\%Lengths);
#	print Dumper($lengths);
	my $value=$line[2];
	unless  ($q){
	    $value =~ s/([0-9]+)q/\1*q/g;
	    $value =~ s/q/1/g;
	    $value =~ s/\^/**/g;
	    $value=eval($value);
	}
	$P[$line[0]]->[$line[1]]=$value;
    }
    close(IN);
    return (\@P,\%Lengths);
}

sub output{
    my $P=shift;

    if (!$quiet){
	print "(*Kazhdan-Lusztig matrix from file $file:*)\n";
	if ($q){
	    print "(*Including q parameter*)\n";
	}else{
	    print "(*Setting q=1*)\n";
	}
	my $size=matrixDim($P);
	print "(*Size of matrix: $size*)\n";
	print "(*This matrix gives the multiplicity of standard modules in irreducibles, up to sign*)\n";
	$transpose and print "(*Displaying the transpose matrix*)\n";
    }
    print "klmatrix = \n";
    $transpose?print matrixString(transpose($P)):print matrixString($P);
    print "\n";

    if (!$q and ($parameter or ($parameter eq '0'))){
	print "(*Irreducible $parameter as sum of standards up to signs:*)\n $parameter= ";
	displayRow($P,$parameter);
	print "\n";
    }
    print "\n";

    my $Psigned;
    if ($signed or $inverse){
	my @Psigned;
	foreach my $i (0..$#{@$P}){
	    my $row=$P->[$i];
	    my @row=@$row;
	    foreach my $j (0..$#row){
		my $sign=(-1)**($Lengths->{$i}-$Lengths->{$j});
		$row[$j]=$row[$j]*$sign;
	    }
	    push @Psigned, \@row;
	}
	$Psigned=\@Psigned;
    }
    if ($signed){
	if ($q){
	    $quiet or print "(* Not printing signed matrix when -q is specificed. *)\n";
	}else{
	    $quiet or print "(* Signed KL matrix: *)\n";
	    $quiet or print "(* This matrix gives the multiplicities of standards in irreducibles *)\n";
	    print "klmatrixsigned =\n";
	    $transpose?print matrixString(transpose($Psigned)):print matrixString($Psigned);
	    print ";\n";
	    if (!$q and ($parameter or ($parameter eq '0'))){
		print "(* Irreducible $parameter as sum of standards:*)\n$parameter= ";
		displayRow($Psigned,$parameter);
		print "\n";
	    }
	    print "\n";
	}
    }
    
    if ($inverse){
	my @save;
	foreach my $row (@$P){
	    my @row=@$row;
	    push @save, \@row;
	}
	my $Pinverse=matinv($Psigned,$Psigned,matrixDim($Psigned),0);
	$P=\@save;

	if ($q){
	    $quiet or print "(* Not printing inverse matrix when -q is specified *)\n";
	}else{
	    print "klmatrixinverse =\n";
	    $transpose?print matrixString(transpose($Pinverse),1):print matrixString($Pinverse,1);
	    print ";\n";
	    if ($parameter  or ($parameter eq '0')){
		print "(* Standard $parameter as sum of irreducibles*):\n$parameter= ";
		displayRow($Pinverse,$parameter);
		print "\n";
	    }
	    print "\n";
	    
	}
	
    }
#    print "test:\n", matrixString(matrixMult($P,$Pinverse));



}

sub displayRow{
    my $P=shift;
    my $r=shift;
    my $n=matrixDim($P);
    my @rv;
    print column;
    foreach my $i (0..$n-1){
	my $coeff=roundToInteger($P->[$r]->[$i]);
	next unless $coeff;
	if ($coeff  == 1){
	    $coeff='';
	}elsif ($coeff == -1){
	    $coeff = '-';
	}else{
	    $coeff = "$coeff".'x';
	}
	print "$coeff".$i." ";

    }
}

sub help{
print "

Compute the Kazhdan-Lusztig matrix for a given block.  The output is
an n by n matrix M; the entry in row i and column j is the
Kazhdan-Lusztig polynomial P_{i,j}. Here i,j run over integers 0 to
n, parametrizing the block, as given by the block command in
atlas.

Row i of this matrix is the multiplicity of standard modules 0..n in
(the formal character of) the irreducible representation i, up to
sign. With the option -s will display the signed matrix.

This script reads a file produced by the klbasis command of atlas as
input.

By default q is set to 1, keep q as an indeterminate with the -q
option.  With q=1 this is the matrix of multiplicities of irreducible
modules in standards.

The -i option (without -q) will give the inverse matrix, which gives
the multiplicities of irreducibles in standards.

Use -t for the transpose.

WARNING: The computation of inverses is very slow and recommended only
for very small rank. For larger matrices computing the inverse in
another program is better.

With the -p option, give information about parameter p. This is the
i^th row of the matrix (or column if -t), and is included for
legibility. For example -p 5 gives the irreducible module 5 as a sum
of standards (up to sign). With the -i option, also gives irreducible
module 5 as a sum of standards. Not allowed with -q.

Output is in mathematica format by default. With -o and not -v, the
resulting file can be read directly into mathematica. For maple use
the -f flag.

Usage:

kl [-t] [-q]  [-i] [-s] [-o outputfile] [-p parameter] [-v] [-f maple|mathematica] file [-h]


The file argument should contain the output of the atlas command klbasis.

Options:

-q             matrix of polynomials in q (default: set q=1)
-i             output inverse matrix also (not allowed with -q)
-s             signed Kazhdan-Lusztig matrix
-p parameter   information about parameter p (p^th column of matrix)
-o outputfile  send output to file instead of terminal
-f maple       maple format (default: mathematica)
-t             output the tranpose matrix
-q             quiet
-h             this help file

Examples:

kl klbasisC2               Kazhdan-Lusztig matrix of integers for type C2
kl -v klbasisC2            more output
kl -q klBasisC2            give Kazhdan-Lusztig matrix of polynomials in q for type C2
kl -t klbasisC2            transpose Kazhdan-Lusztig matrix
kl -s klbasisC2            also give the signed matrix
kl -i klbasisC2            also give the inverse matrix
kl -o kloutput klBasisC2   send output to file kloutput
kl -f maple klbasisC2      matrices in maple format
kl -p 5 klbasisC2          also give irreducible module 5 as a sum of standards (up to signs)
kl -p 5 -s klbasisC2       also give irreducible module 5 as a sum of standards with correct signs
kl -p 5 -i klbasisC2       also give standard module 5 as a sum of irreducibles

Sample output:

Note: For the signs we need to compute (-1)^(l(x)-l(y)). This is
available from the block command.  According to an educated guess by
David Vogan we can also get it from the output of klbasis. That is
what we do here. See the computation of the Lengths hash.

Send questions and bugs to jda\@math.umd.edu
"; exit; }

sub transpose{
    my $matrix=shift;
    my $rows=scalar(@$matrix);
    my $columns=scalar(@{$matrix->[0]});
    my @rv;
    foreach my $i (0..$columns-1){
	my @row;
	foreach my $j (0..$rows-1){
	    $row[$j]=$matrix->[$j][$i];
	}
	push @rv, \@row;
    }
    return \@rv;
}

sub matrixString{
    my $string;
    my $m=shift;
    my $round=shift;
    my @matrix=@$m;
    my @output;
    foreach my $row (@matrix){
	my @row=@$row;
	if ($round){
	    foreach my $i (0..$#row){
		$row->[$i]=roundToInteger($row->[$i]);
	    }
	}
	my $line= join ',', @$row;
	push @output,$line;
    }
    my $open='{';
    my $close='}';
    if ($format =~ /map/){
	$open='[';
	$close=']';
    }
    $string.="$open$open";
    $string.=join "$close,\n $open", @output;
    $string.="$close$close";
    return $string;
}    

sub matrixMult {
    my ($A,$B) = @_;
    my ($Arows,$Acols) = matrixDim($A);
    my ($Brows,$Bcols) = matrixDim($B);

    unless ($Acols == $Brows){
	print "ERROR";
        die "Error: size of matrices don't match: $Acols != $Brows";
    }

    my $result = [];
    my ($i, $j, $k);

    for $i (range($Arows)) {
        for $j (range($Bcols)) {
            for $k (range($Acols)) {
                $result->[$i][$j] = roundToInteger($result->[$i][$j]+$A->[$i][$k] * $B->[$k][$j]);
            }
        }
    }
    return $result;
}

sub roundToInteger{
    my $a=shift;
    return int($a+.5*($a<=>0));
}



sub range { 0 .. ($_[0] - 1) }

sub matrixDim {
    my $matrix = $_[0];
    my $rows = scalar(@$matrix);
    my $cols = scalar(@{$matrix->[0]});
    return ($rows, $cols);
}

sub matinv
{
	my ($aref, $bref, $n, $quiet) = @_;
	my ($row_start, $i);
	my $twon = $n + $n;
	my @ai;

	die "matinv():  Need as arguments two matrix references and dimensions.\n"
		unless defined $n;

	# First, paste the input and the identity side by side.
	$row_start = 0;
	for (my $i = 0; $i < $n; $i++) {
		for (my $j = 0; $j < $n; $j++) {
			$ai[$i][$j] = $$aref[$i][$j];
		}
	}
	for (my $i = 0; $i < $n; $i++) {
		for (my $j = 0; $j < $n; $j++) {
			$ai[$i][$j+$n] = 0.0;
		}
	}
	for (my $i = 0; $i < $n; $i++) {
		$ai[$i][$i+$n] = 1.0;
	}

	# Second, use Householder transformations to put it into upper-triangular
	# form.
	for (my $i = 0; $i < $n; $i++) {
		if (!$quiet) {
			my @Q = householder_on_submatrix_Q(\@ai, $n, $twon, $i);
#			print "-- matinv Q $i:\n";
#			print_matrix(\@Q, $n, $n);
#		print  Dumper(\@Q);
#			print "\n";

#			print "-- matinv aug $i:\n";
#			print_matrix(\@ai, $n, $twon);
#			print Dumper(\@ai);
#			print "\n";
		}
		householder_on_submatrix(\@ai, $n, $twon, $i);
	}

	# Third, put 1 on the left diagonal.
	for (my $i = 0; $i < $n; $i++) {
		my $d = $ai[$i][$i];
		if ($d == 0.0) {
			die "Singular.\n";
		}
#		elsif (tol_zero($d)) {
#			die "Nearly singular.\n";
#		}
		for (my $j = 0; $j < $twon; $j++) {
			$ai[$i][$j] = $ai[$i][$j] / $d;
		}
	}

	# Fourth, clear out the rest of the left-hand side.
	# 1 . . . .  . . . . .
	# 0 1 . . .  . . . . .
	# 0 0 1 . .  . . . . .
	# 0 0 0 1 .  . . . . .  <-- i
	# 0 0 0 0 1  . . . . .  <-- i2

	for (my $i = $n-2; $i >= 0; $i--) {
		for (my $i2 = $n-1; $i2 > $i; $i2--) {
			my $mul = $ai[$i][$i2];
			for (my $j = 0; $j < $twon; $j++) {
				$ai[$i][$j] -= $ai[$i2][$j] * $mul;
			}
		}
	}

	# Fifth, obtain the inverse from the right-hand side.
	for (my $i = 0; $i < $n; $i++) {
		for (my $j = 0; $j < $n; $j++) {
			$$bref[$i][$j] = $ai[$i][$j+$n];

		}
	}
	return $bref;
}
sub householder_on_submatrix
{
	my ($aref, $nr, $nc, $a) = @_;
	my $height = $nr - $a;
	my ($i, $j);

	die "householder_on_submatrix():  Need as arguments matrix reference, dimensions, and start column.\n"
		unless defined $a;

	# --- First column ---
	# Get v
	my @v = get_submatrix_column($aref, $nr, $nc, $a, $a);

	# Compute ||v||; write w.
	my $normv = sqrt(dot(\@v, \@v, $height));
	if ($v[0] >= 0.0) {
		$normv = -$normv;
	}

	my @w;
	$w[0] = $normv;
	for ($i = 1; $i < $height; $i++) {
		$w[$i] = 0.0;
	}

	# Compute z = v - w.
	my @z;
	for ($i = 0; $i < $height; $i++) {
		$z[$i] = $v[$i] - $w[$i];
	}

	put_submatrix_column($aref, $nr, $nc, $a, $a, \@w);

	# Compute 2 * (z dot z)
	my $zdotz = dot(\@z, \@z, $height);
	my $tnz2;
	if (abs($zdotz) < $pml_tolerance) {
		$tnz2 = 0.0;
	}
	else {
		$tnz2 = 2.0 / $zdotz;
	}

	# --- Subsequent columns ---
	my $b;
	for ($b = $a + 1; $b < $nc; $b++) {
		# Get column c.
		my @c = get_submatrix_column($aref, $nr, $nc, $b, $a);

		# Compute the scalar z^t c
		my $ztc = dot(\@z, \@c, $height);

		# Compute d = (-2 / z dot z) * z^t c
		my $d = -$tnz2 * $ztc;

		# New column is q = c + d * z.
		for ($j = 0; $j < $height; $j++) {
			$w[$j] = $c[$j] + $d * $z[$j];
		}

		put_submatrix_column($aref, $nr, $nc, $b, $a, \@w);
	}
}
sub householder_on_submatrix_Q
{
	my ($aref, $nr, $nc, $a) = @_;
	my $height = $nr - $a;
	my ($i, $j);

	die "householder_on_submatrix_Q():  Need as arguments matrix reference, dimensions, and start column.\n"
		unless defined $a;

	# --- First column ---
	# Get v
	my @v = get_submatrix_column($aref, $nr, $nc, $a, $a);

	# Compute ||v||; write w.
	my $normv = sqrt(dot(\@v, \@v, $height));
	if ($v[0] >= 0.0) {
		$normv = -$normv;
	}

	my @w;
	$w[0] = $normv;
	for ($i = 1; $i < $height; $i++) {
		$w[$i] = 0.0;
	}

	# Compute z = v - w.
	my @z;
	for ($i = 0; $i < $height; $i++) {
		$z[$i] = $v[$i] - $w[$i];
	}

	# Compute 2 / z dot z
	my $zdotz = dot(\@z, \@z, $height);
	my $tnz2;
	if (abs($zdotz) < $pml_tolerance) {
		$tnz2 = 0.0;
	}
	else {
		$tnz2 = 2.0 / $zdotz;
	}


	# Q = I - (2 z z^t / z^t z)
	my @Q;

	for ($i = 0; $i < $nr; $i++) {
		for ($j = 0; $j < $nc; $j++) {
			$Q[$i][$j] = 0.0;
		}
	}
	for ($i = 0; $i < $nr; $i++) {
		$Q[$i][$i] = 1.0;
	}

	my @zzt = outer(\@z, \@z, $height, $height);

	for ($i = $a; $i < $nr; $i++) {
		for ($j = $a; $j < $nr; $j++) {
			$Q[$i][$j] -= $zzt[$i-$a][$j-$a] * $tnz2;
		}
	}

	return @Q;
}
sub householder_on_submatrix_Q
{
	my ($aref, $nr, $nc, $a) = @_;
	my $height = $nr - $a;
	my ($i, $j);

	die "householder_on_submatrix_Q():  Need as arguments matrix reference, dimensions, and start column.\n"
		unless defined $a;

	# --- First column ---
	# Get v
	my @v = get_submatrix_column($aref, $nr, $nc, $a, $a);

	# Compute ||v||; write w.
	my $normv = sqrt(dot(\@v, \@v, $height));
	if ($v[0] >= 0.0) {
		$normv = -$normv;
	}

	my @w;
	$w[0] = $normv;
	for ($i = 1; $i < $height; $i++) {
		$w[$i] = 0.0;
	}

	# Compute z = v - w.
	my @z;
	for ($i = 0; $i < $height; $i++) {
		$z[$i] = $v[$i] - $w[$i];
	}

	# Compute 2 / z dot z
	my $zdotz = dot(\@z, \@z, $height);
	my $tnz2;
	if (abs($zdotz) < $pml_tolerance) {
		$tnz2 = 0.0;
	}
	else {
		$tnz2 = 2.0 / $zdotz;
	}


	# Q = I - (2 z z^t / z^t z)
	my @Q;

	for ($i = 0; $i < $nr; $i++) {
		for ($j = 0; $j < $nc; $j++) {
			$Q[$i][$j] = 0.0;
		}
	}
	for ($i = 0; $i < $nr; $i++) {
		$Q[$i][$i] = 1.0;
	}

	my @zzt = outer(\@z, \@z, $height, $height);

	for ($i = $a; $i < $nr; $i++) {
		for ($j = $a; $j < $nr; $j++) {
			$Q[$i][$j] -= $zzt[$i-$a][$j-$a] * $tnz2;
		}
	}
	foreach  my $i (0..$a){
	    foreach my  $j (0..$a){
#		print "fix $i,$j";
		if (abs($Q[$i][$j])<$pml_tolerance){
#		    print "FIXED $Q[$i][$j]\n";;
#		    $Q[$i][$j]=0 
		    }
	    }
	}
	return @Q;
}
sub get_submatrix_column
{
	my ($aref, $nr, $nc, $colidx, $start_row) = @_;
	my ($src, $dst);
	my @submatrix_column;
	for ($src = $start_row, $dst = 0; $src < $nr; $src++, $dst++) {
		$submatrix_column[$dst] = $$aref[$src][$colidx];
	}
	return @submatrix_column;
}
sub dot
{
	my ($uref, $vref, $n) = @_;

	die "dot():  Need as arguments two vector references and length.\n"
		unless defined $n;

	my $sum = 0.0;
	my $i;
	for ($i = 0; $i < $n; $i++) {
		$sum += $$uref[$i] * $$vref[$i];
	}
	return $sum;
}
sub outer
{
	my ($uref, $vref, $m, $n) = @_;

	die "outer():  Need as arguments two vector references and lengths.\n"
		unless defined $n;

	my ($i, $j, @uv);
	for ($i = 0; $i < $m; $i++) {
		for ($j = 0; $j < $n; $j++) {
			$uv[$i][$j] = $$uref[$i] * $$vref[$j];
		}
	}
	return @uv;
}
sub put_submatrix_column
{
	my ($aref, $nr, $nc, $colidx, $start_row, $pcolref) = @_;
	my ($src, $dst);
	for ($src = 0, $dst = $start_row; $dst < $nr; $src++, $dst++) {
		$$aref[$dst][$colidx] = $$pcolref[$src];
	}
}