#!/usr/bin/perl
#compute coherent continuation using atlas output
#jda@math.umd.edu

use strict;
use Data::Dumper;


use Getopt::Long;
use vars qw($help $outputFile $parameter $w $verbose $standard $irreducible $cell $format);

my %options=('h' => \$help,  
	     'o' => \$outputFile,
	     'p' => \$parameter,
	     'w' => \$w,
	     's' => \$standard,
	     'i' => \$irreducible,
	     'c' => \$cell,
	     'f' => \$format,
	     'v' => \$verbose);

GetOptions(\%options,qw(h! o=s p=s w=s v=s s! i! c! f=s));
$verbose=1 unless ($verbose ge 0);
my $file=$ARGV[0];

&help if (!$file or $help);

if ((($irreducible and ($standard or $cell)) or ($standard and $cell) or !($irreducible or $standard or $cell))){
    print "You must choose exactly one of -i (irreducible) or -s (standard) or -c (cell)
See -h for help.\n";
    exit;
}
if ($outputFile){
    open(OUT,">$outputFile")||die("Can't open $outputFile for output");
    select(OUT);
}

my $open='{';
my $close='}';
if ($format =~ /map/){
    $open='[';
    $close=']';
}


my @representations;

if ($standard){
    my $rep=processStandard($file);
    push @representations, $rep;
#    my $m=multiply($operators,$w) if $w;
#    my $action=action($m,$parameter) if (($parameter)||($parameter eq '0'));
#    output($operators,$m,$action);
}elsif ($irreducible){
    my $data=readData($file);
    my $entry=$data->[0];
    my $rep=processIrreducible($entry);
    push @representations, $rep;
#    my $m=multiply($operators,$w) if $w;
#    my $action=action($m,$parameter) if (($parameter)||($parameter eq '0'));
#    output($operators,$m,$action);
}elsif ($cell){
    my $data=readData($file);
    foreach my $entry (@$data){
	my $rep=processIrreducible($entry);
	unless (scalar(@$rep)){
	    $verbose and print "Skipping trivial representation\n";
	    next;	
	}
	push @representations, $rep;
    }
    if ($verbose){
	 print "Computed ",scalar(@representations)," representations of dimensions ";
	 my @dims;
	 foreach my $rep (@representations){
	     my $row = $rep->[0]->[0];
	     push @dims, scalar(@$row);
	 }
	 print join ',', @dims;
	 print "\n";
     }

}
#print Dumper(\@representations);
output(\@representations);

close(OUT) if $outputFile;

sub action{
    my $m=shift;
    $m||return '';
    my $parameter=shift;
    my $n=matrixDim($m);
    my @action;
    foreach my $i (0..$n){
	my $x=$m->[$i]->[$parameter];
	push @action,$m->[$i]->[$parameter];
    }
    return \@action;
}

sub readData{
#read data from wgraph or wcells file for passing to processStandard
    open(IN,"<$file")||die("Can't open $file for input");
    my @strings=();
    my $string;
    foreach my $line (<IN>){
	if ($line =~ /cell/){
	    push @strings, $string if $string;
	    $string='';
	}else{
	    $string .= $line;
	}
    }
    push @strings, $string;
#    print Dumper(\@strings);
    return \@strings;
}

sub processStandard{
    my $file=shift;
    #this should be the output of block
    my @operators;
    my $rank;
    open(IN,"<$file")||die("Can't open $file for input");
    my $size=0;
    #one run to get size and rank
    foreach my $line (<IN>){
	$line =~ s/\#.*//;
	next unless $line;
#format of block output changed:
# 0( 0,6):   0   0  [i1,i1]   1   2   ( 6, *)  ( 4, *) new
# 0( 0,6):   1   2    ( 6, *)  ( 4, *)    [i1,i1]      old

#	my ($number,$cross,$cayley,$type)= ($line =~/^([0-9\s]*).*:\s*([^\(]*)\(\s*(.*)\)\s*\[([^\]]*)/g);
	my ($number,$cartanAndLength,$type,$cross,$cayley)= ($line =~/^([0-9\s]*).*:\s*([^\[]*)\[\s*(.*)\]\s*([^\(]*)\((.*)\)/g);
	$size=$number if ($number >$size);
	unless ($rank){
	    my @cayley=split(/\)\s*\(/, $cayley);	
	    $rank=scalar(@cayley);
	}
#	print "size:$size,rank:$rank\n";
    }
    close(IN);
    my @matrices;
    foreach my $i (0..$rank-1){
	my @zeroMatrix=();
	foreach my $j (0..$size){
	    push @zeroMatrix, [(0)x($size+1)];
	}
	push @matrices, \@zeroMatrix;
    }
#    print Dumper(\@matrices);;
    open(IN,"<$file")||die("Can't open $file for input");    
    foreach my $line (<IN>){
	#10(10,0):  11  10    ( *, *)  ( *, *)    [r2,r1]  3  1212
#format of block output changed, see above	
#	my ($number,$cross,$cayley,$type)= ($line =~ /^([0-9\s]*).*:\s*([^\(]*)\(\s*(.*)\)\s*\[([^\]]*)/g);
	my ($number,$cartanAndLength,$type,,$cross,$cayley)= ($line =~/^([0-9\s]*).*:\s*([^\[]*)\[\s*(.*)\]\s*([^\(]*)\((.*)\)/g);

	$cross =~ s/\s*$//;
#	print "\n               cross is   $cross\n";
	my @cross=split(/\s+/, $cross);
#	print "\n";
	$cayley =~ s/\s//g;
	my @cayley=split(/\)\s*\(/, $cayley);
#	print "cayley:", join ':', @cayley;
#	print "\n";
	my @type=split(',', $type);
#	print "type:", join ':', @type;
#	print "\n";
#	print "DOING $number\n";
	foreach my $i (0..$rank-1){
	    my $root=$type[$i];
#	    print "root:$root\n";
	    if ($root =~ /r|C/){
		#just cross action
		my $j=$cross[$i];
		$matrices[$i]->[$j]->[$number]=1;
	    }elsif ($root =~ /ic/){
		#-1 on diagonal
		$matrices[$i]->[$number]->[$number]=-1;
	    }elsif ($root =~ /i1|i2/){
		#cayleys - 
		my $j=$cross[$i];
		my ($k,$l)=split ',', $cayley[$i];
#		print "j,k,l:$j,$k,$l\n";
		$matrices[$i]->[$j]->[$number]=-1;
		unless ($k =~ /\*/){
		    $matrices[$i]->[$k]->[$number]=1;
#		    print "GOT $i,$k,$number\n";
		}
		unless ($l =~ /\*/){
		    $matrices[$i]->[$l]->[$number]=1;
#		    print "GOT ALSO $i,$j,$number\n";
		}
	    }else{
		print "MISSED CASE\n";exit;
	    }
	}
    }
    return \@matrices;
}

sub processIrreducible{
#string (not a file) from readData    
    my $string=shift;
    my @lines = split "\n", $string;
    my %A;
    my %B;
    my @operators;
    my $rank=0;
    my $size=0;

#1:{2}:{(0,1),(2,1)}  old
#1[5]: {2} --> 0,2    new
#my ($number,$cartanAndLength,$type,,$cross,$cayley)= ($line =~/^([0-9\s]*).*:\s*([^\[]*)\[\s*(.*)\]\s*([^\(]*)\((.*)\)/g);
    foreach my $line (@lines){    
#	print "line:$line\n";
#	next unless ($line =~ /-->/);
	my ($first,$second,$third)= split ':', $line;
#	print "\n\$first:$first\n\$second:$second\n\$third:$third\n";
	$size=$first+1 if ($first+1>$size);
	$second=~ s/[\{|\}]//g;
	my @b=split ',', $second;
	foreach my $i (@b){
	    push @{$A{$i}}, $first;
	    $rank = $i if ($i>$rank);
	}
	$third =~ s/^\{//;
	$third =~ s/\}$//;
	$third =~ s/^\(//;
	$third =~ s/\)$//;
	chomp($third);
	next unless ($third);
	my @c=split '\),\(', $third;
	foreach my $entry (@c){
	    my ($x,$y)=split ',', $entry;
	    push @{$B{$x}}, [$first,$y];
	}
    }


#     foreach my $line (@lines){
# 	my ($a,$b,$c)=split ':', $line;
# 	$size=$a+1 if ($a+1>$size);
# 	$b=~ s/[\{|\}]//g;
# 	my @b=split ',', $b;
# 	foreach my $i (@b){
# 	    push @{$A{$i}}, $a;
# 	    $rank = $i if ($i>$rank);
# 	}
# 	$c =~ s/^\{//;
# 	$c =~ s/\}$//;
# 	$c =~ s/^\(//;
# 	$c =~ s/\)$//;
# 	chomp($c);
# 	next unless ($c);
# 	my @c=split '\),\(', $c;
# 	foreach my $entry (@c){
# 	    my ($x,$y)=split ',', $entry;
# 	    push @{$B{$x}}, [$a,$y];
# 	}
#     }


#    print Dumper(\%A);
#    print "B:",Dumper(\%B);
#    print "size$size rank$rank\n";



    foreach my $i (1..$rank){
	my @matrix;
	foreach my $j (0..$size-1){
	    push @matrix, [(0)x$size];
	}
	foreach my $j (0..$size-1){
	    if (grep {/^$j$/} @{$A{$i}}){
		$matrix[$j]->[$j]=-1;
	    }else{
		$matrix[$j]->[$j]=1;
	    }
	}
	foreach my $k (0..$size-1){
	    next unless $B{$k};
	    next if (grep {/^$k$/} @{$A{$i}});
	    my @r=@{$B{$k}};
	    foreach my $entry (@r){
		my ($x,$y)=@$entry;
#		next unless (grep {/$x/} @{$A{$i}});
		next unless (grep {/^$x$/} @{$A{$i}});
		$matrix[$x]->[$k]=$y;
	    }
	}
	push @operators, \@matrix;
    }
    return \@operators;
}

sub multiply{
    my $operators=shift;
    my @operators=@$operators;
    my $w=shift;
    unless ($w){
	my $n=matrixDim($operators[0]);
	my @id;
	foreach my $i (1..$n){
	    my @row;
	    foreach my $j (1..$n){
		if ($i == $j){
		    push @row, 1;
		}else{
		    push @row, 0;
		}
	    }
	    push @id, \@row;
	}
#	    print "GOT ID:", Dumper(\@id);
	return \@id;
    }
	
    my @seq=split ',', $w;
#    print "sequence:", join ',', @seq;
    my $n=$#seq;
    if ($seq[$n] > scalar(@operators)){
	print "Error: entry $seq[$n] of $w is greater than the rank ".scalar(@operators)."\n";
	return '';
    }
    my $matrix=$operators[$seq[$n]-1];
#    print Dumper($matrix);
#    print matrixString($matrix);
    foreach my $i (1..$#seq){
	if ($seq[$i] > scalar(@operators)){
	    print "Error: entry $seq[$i] of $w is greater than the rank ".scalar(@operators)."\n";
	    return '';
	}
	$matrix=matrixMult($operators[$seq[$n-$i]-1],$matrix);
    }
    return $matrix;
}



sub matrixString{
    my $string;
    my $m=shift;
    my @matrix=@$m;
    my @output;
    foreach my $row (@matrix){
	my $line= join ',', @$row;
	push @output, $line;
    }
    $string.="$open$open";
    $string.=join "$close,\n$open", @output;
    $string.="$close$close";
    return $string;
}    

sub output{
    my $representations=shift;    
    my @representations=@$representations; #array of array refs
    return unless scalar(@representations);
    if ($verbose){
	my $type= $standard?'standard':'irreducible';
	print "Coherent continuation action in basis of $type modules\nUsing file $file\n";
    }
    if (($verbose>1 or !($parameter ge 0))){
	print "verbose $verbose and parameter $parameter\n";
	if (scalar(@representations) ==1){
	    print "operators:=\n";
	}else{
	    print "operators:=$open\n";
	}
	
	my @strings;
	foreach my $operators (@representations){
	    my @operators=@$operators;
	    my @operatorsStrings=();
	    foreach my $i (0..$#operators){
		push @operatorsStrings, matrixString($operators[$i]);
	    }
	    my $ostring= join ",\n", @operatorsStrings;
	    push @strings, $ostring;
	}
	print join ",\n\n", @strings;


	if (scalar(@representations) ==1){
	    print ";\n";
	}else{
	    print "\n$close;\n\n";
	}
    if ($w){
	if (scalar(@representations) ==1){
	    print "w:=\n";
	}else{
	    print "w:=$open\n";
	}
	my @strings;
	foreach my $operators (@representations){
	    my $m=multiply($operators,$w);
	    push @strings, matrixString($m);
	}
	print join ",\n\n", @strings;
    }

    if (scalar(@representations) ==1){
	print ";\n";
    }else{
	print "\n$close;\n\n";
    }

    }else{
	print "not printing matrices (use -v 2 to override)\n" unless ($verbose==0);
    }

    if (!($cell) and (($parameter)||($parameter eq '0'))){
	my $ops=$representations[0];
	my $m=multiply($ops,$w);

	my $action=action($m,$parameter);
	my @action=@$action;
 	$verbose and print "action of s_$w on parameter $parameter:\n";
 	print  "s_\{$w\}:$parameter -> ";
 	foreach my $i (0..$#action){
 	    next unless ($action[$i]);
 	    my $coeff=$action[$i];
 	    if ($coeff  == 1){
 		$coeff='';
 	    }elsif ($coeff == -1){
 		$coeff = '-';
 	    }else{
 		$coeff = "$coeff".'x';
 	    }
 	    print "$coeff".$i." ";
 	}
 	print "\n";
     }
}


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] += $A->[$i][$k] * $B->[$k][$j];
            }
        }
    }
    return $result;
}

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

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

sub help{
print "

Compute the coherent continuation, in 3 situations:

1) on an entire block in terms of the basis of standard representations (block)
2) on an entire block in terms of the basis of irreducible representations (wgraph)
3) on a list of cells (wcells)

In each case the computation is based on the output of the atlas

In the case of irreducible or standard modules, gives a list of
matrices, one for each simple root. In the case of wcells, a list of
such lists, one for each cell. (The trivial representations is
skipped).

Optionally you may specify a Weyl group element with -w, and get the
matrix of w in the given representation (or list of such in the case
of cells).  Except in the case of cells you may also give (a Weyl
group element and) a parameter, and get the image of the parameter
under the Weyl group element.

Usage:

coherentContinuation [-s|i|c] [-v] [-o outputFile] file
coherentContinuation [-s|i|c] [-v] [-o outputFile] file [-w list]
coherentContinuation [-s|i] [-v] [-o outputFile] file [-w list -p parameter]
coherentContuation [-h]: show this help file

The required options -s, -i or -c specify standard, irreducible or
cell representations. The file argument is the output of the block
command (in the case of -s), wgraph (in the case of -i), or wcells
(-c). Vectors are viewed as columns, and the matrix acts on the left,
the image of the i^th basis vector is given by the i^th column of the
matrix. For example the matrix

{{1,0,0}
{0,1,0}
{1,1,-1}}

takes 1 -> 1 + 3
      2 -> 2 + 3
      3 -> -3

A Weyl group element is specificed by -w x_1,x_2,... - i.e. a
comma-separated list of integers, each running from 1 to the
semisimple rank, corresponding to a product of simple reflections. The
product of the corresponding operators is displayed.

The optional argument p (-w must be set also) specifies a parameter
from list of parameters (starting at 0).

Options:

-i            basis of irreducible modules, using the output of wgraph
-s            basis of standard modules, using the output of block
-c            cell representation, using the output of wcells
-v            verbose output
-w list       element of Weyl group, e.g. 1,2,1,3,2
-p parameter  parameters is a number from list of parameters in output of block command
-o outputFile 
-h            this help file

Output is in mathematica format. If you select an output file with -o,
and not -v or -p, the resulting file may be read directly into
mathematica.

Examples:

coherentContinuation -i wgraphA2                basis of irreducibles (2 6x6 matrices)
coherentContinuation -s blockA2                 basis of standards (2 6x6 matrices)
coherentContinuation -s blockA2  -w 1,2,1       also give action of long element of W
coherentContinuation -s blockA2  -w 1,2,1 -p 2  also give action of w on parameter 2
coherentContinuation -c wcellsA2                list of list of matrices, one for each cell
coherentContinuation -c wcellsA2 -w 1,2,1       also action of w in each of these representations

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