#!/usr/bin/perl

use strict;
use Data::Dumper;
use PDL;

my $klCommand='/home/jda/atlasTables/kl';
my $smithNormalFormNotebook='/home/jda/stable/IntegerSmithNormalForm.m';
my $mathematicaFile='out';

use Getopt::Long qw(:config no_ignore_case);
use vars qw($help  $Help $outputFile $parameters $v $blockFile $outputFile $klBasisFile $subsetSimpleRoots 
$wCellsFile $cells $all $aq $blockuFile $exact $dual $noStable);
#NOTE: S gives simple roots of parameter on wall; S is labelled 1,2,...,n not 0,1,...,n-1

my %options=('h' => \$help,  
	     'd' => \$dual,
	     'H' => \$Help,  
	     'e' => \$exact,
	     'q' => \$aq,
	     'a' => \$all,
	     'S' => \$subsetSimpleRoots,
	     'p' => \$parameters,
	     'c' => \$cells,
	     'b' => \$blockFile,
	     'w' => \$wCellsFile,
	     'k' => \$klBasisFile,
	     'o' => \$outputFile,
	     'n' => \$noStable,
	     'v' => \$v);
my $argument=join ' ', @ARGV;  #save for output later

unless (scalar(@ARGV)){
    print "Default action: processing files \"block\" and \"klbasis\", all representations, at rho\n";
    print "For help use -h\n";
}
GetOptions(\%options,qw(h! b=s o=s v=s p=s k=s S=s w=s c=s a! H! q! e! d! n!));

$blockFile||='block';

if ($dual){
    $wCellsFile||='dualwcells';
    $klBasisFile||='dualklbasis';
    $blockuFile||='dualblocku';
}else{
    $wCellsFile||='wcells';
    $klBasisFile||='klbasis';
    $blockuFile||='blocku';
}
$cells = '' unless ($cells or ($cells eq '0'));
$subsetSimpleRoots||='';
my %dimensions;
$v=1 unless ($v or ($v eq 0));
$v>0 and print "%stable $argument\n";
&moreHelp if ($Help);
&help if ($help);
$aq and &aq();

my @subsetSimpleRoots=split ',', $subsetSimpleRoots;
$subsetSimpleRoots=\@subsetSimpleRoots;
my ($blockData,$allParameters)=&parseBlock($blockFile);

if ($all){
    my @allCells=split ',', $cells;
    my $n=scalar(@allCells);
    foreach my $j (1..2**$n-1){
	my @cells;
	foreach my $i (0..$#allCells){
	    if (2**$i&$j){
		push @cells, $allCells[$i-1];
	    }
	}
	$cells=join ',', @cells;
	my ($map,$dualMap);
	($parameters,$map,$dualMap)=getParameters($cells);
	&main($cells,$parameters,$subsetSimpleRoots,$map,$dualMap);
    }
    print "\nSummary:\nAll combinations of cells with non-zero stable sums\n";
    print "All cells: $cells\n";
    print "lambda is singular at simple roots: ", join ',', @$subsetSimpleRoots;
    print "\nAll irreducibles (living at lambda): ", join ',',@$parameters;
    print "\nCells\t\tDimension of space of stable sums\n";

    foreach my $key (keys(%dimensions)){
	my @key= split ',', $key;
	@key = sort {$a <=> $b} @key;
	print join ',', @key;
	print "\t\t$dimensions{$key}\n";
    }    
    print "\n";
    exit;
}else{
    my ($map,$dualMap);
    ($parameters,$map,$dualMap)=getParameters($cells);
    &main($cells,$parameters,$subsetSimpleRoots,$map,$dualMap);
    print "\n";
    exit;
}

sub getParameters{
    $v>3 and print "sub getParameters\n";
    my $cells=shift;
    if ($cells or ($cells eq '0')){
	my $cellsData=parseWCells($wCellsFile);
	my @cells=split ',', $cells;
	$cells=\@cells;
	my @parametersFromFile;
	foreach my $cell (@cells){
	    my $parameters=$cellsData->{$cell};
	    push @parametersFromFile, @$parameters;
	}
	$parameters=\@parametersFromFile;

    }elsif($parameters or ($parameters eq '0')){
	my @parameters=split ',', $parameters;
	$parameters = sort {$a <=> $b} @parameters;
	$parameters = \@parameters;
    }else{
#	print Dumper($blockData);
	my @parameters=keys(%$blockData);
	@parameters = sort {$a <=> $b} @parameters;
 	$parameters=\@parameters;
    }
    my ($map,$dualMap)=parseDualMap();
    if ($dual){
	$v>1 and print "Original parameters (for dual group):", join ',', @$parameters;
	my @dualParameters;
	foreach my $p (@$parameters){
	    push @dualParameters, $map->{$p};
	}
	$parameters=\@dualParameters;
	$v>1 and print "\nRevised parameters (for G): ", join ',', @$parameters;
    }
    return ($parameters,$map,$dualMap);
}

sub parseDualMap{
    $v>3 and print "sub parsedualMap\n";
    open(IN,"<dualmap")||die("Can't open dualmap for input");
    my %map;
    my %dualMap;
    my $permutation;
    foreach my $line (<IN>){
	chomp $line;
	next unless $line;
	$line =~ s/[\]\[]//g;
	$line =~ s/\s//g;
	$permutation .= $line;
    }
    my @p=split ',', $permutation;
    my $n=scalar(@p);
    $permutation=zeroes($n,$n) unless $noStable;
    foreach my $i (0..$#p){
	$map{$i}=$p[$i];
	$dualMap{$p[$i]}=$i;
	$permutation->slice("$i,$p[$i]") .=1 unless $noStable;
    }
    $v>4 and print "map:", Dumper(\%map);
    $v>4 and print "dual map:", Dumper(\%dualMap);
    $v>4 and print "permutation:$permutation";
    return (\%map,\%dualMap,$permutation);
}

sub main{
    $v>3 and print "sub main\n";
    print "\n";
    my ($cells,$parameters,$subsetSimpleRoots,$map,$dualMap)=@_;
    my $cellsString=$cells;
    $v>1 and print "Computing stable sums for cells $cellsString\n";
    my $sizeBlock=scalar(keys %$blockData);
    my $tau;
    my $dualParameters=$allParameters; #those parameters on dual side whose tau invariant contains S
                                       #on the G side these are all parameters that live on S
    if ($subsetSimpleRoots){
	$tau=computeTau($allParameters,$blockData);  #hash of tau invariants of all parameters
	$parameters=pushToSingular($parameters,$subsetSimpleRoots,$tau);
	$dualParameters=pushToSingular($allParameters,$subsetSimpleRoots,$tau);
    }
    if ($v>2){
	print "\nParameters living at lambda/roots\n";
	my @keys = sort {$a<=>$b} @$parameters;
	foreach my $key (@keys){
	    my $roots=$blockData->{$key}[3];
	    print "$key:\t$roots\n";
	}
    }

    $v>1 and print "\nDual Parameters: ", join ',', @$dualParameters;
    $v>1 and print "\n";
    my $parametersAmongDualParameters=zeroes(scalar(@$dualParameters),scalar(@$parameters));
    my %newLabels;
    foreach my $k (0..$#{@$dualParameters}){
	my $param=$dualParameters->[$k];
	$newLabels{$param}=$k;
    }
    foreach my $i (0..$#{@$parameters}){
	my $param=$parameters->[$i];
	my $column=$newLabels{$param};
	$parametersAmongDualParameters->slice("$column,$i") .= 1;
    }

    my ($dimension,$basis,$sizeW,$junk,$everythingIsStable);
    unless ($noStable){
    my ($dualKLMatrix)=&klMatrix($blockData);
    $dualKLMatrix=transpose($dualKLMatrix);

    my $pairingStableDualParameters=pairingStableDualParameters($blockData,$dualParameters,$dualKLMatrix);
    $v>2 and print "dual klmatrix:\n", $dualKLMatrix,"\n";
    $v>2 and print "pairingStableDualParameters:\n", pdl($pairingStableDualParameters), "\n";
    my ($columns,$rows)=$pairingStableDualParameters->dims;    
    $v>1 and print "Pairing matrix: $rows rows/$columns columns\n";

    if ($rows >0){
	open(OUT,">$mathematicaFile")||die("Can't open $mathematicaFile for output");
	print OUT "(*output of stable, called via:\nstable $argument *)\n\n";
	print OUT "<<$smithNormalFormNotebook\n";
	print OUT "Parameters:={";
	print OUT join ',', @$parameters;
	print OUT "};\n";
	print OUT "ParametersAmongDualParameters:=";
	print OUT mathematicaMatrix($parametersAmongDualParameters);
	print OUT ";\n";
	print OUT "Pairing:=";
	print OUT mathematicaMatrix($pairingStableDualParameters);
	print OUT ";\n\n";
	if ($exact){
	    $v>1 and print "Using Smith normal form\n";
	    print OUT "SV:=ExtendedSmithForm[V];\n";
	    print OUT "PV:=SV[[2]][[1]];\n";
	    print OUT "QV:=SV[[2]][[2]];\n";
	    print OUT "DV:=SV[[1]];\n";
	    print OUT "Print[\"Smith form of V:\"];\n";
	    print OUT "Print[MatrixForm[QV]];\n";
	    print OUT "Print[MatrixForm[DV]];\n";
	    print OUT "KDV:=NullSpace[Transpose[DV]];\n";
	    print OUT "If[Length[KDV]>0,\n{\n";
	    print OUT "WA:=Transpose[KDV.PV];\n";
	    print OUT "Print[\"Basis of vanishing at 0 in IRREDUCIBLE coordinates:\"];\n";
	    print OUT "Print[MatrixForm[WA]];\n";
	    print OUT "Print[\"Each row of WA gives a sum of irrs (dual side) which vanishes near 0\"];\n";
	    print OUT "W:=ParametersAmongDualParameters.WA;\n";
	    print OUT "Print[\"W is a subset of rows of WA:\"];\n";
	    print OUT "Print[W];\n";
	    print OUT "SW:=ExtendedSmithForm[W];\n";
	    print OUT "PW:=SW[[2]][[1]];\n";
	    print OUT "QW:=SW[[2]][[2]];\n";
	    print OUT "DW:=SW[[1]];\n";
	    print OUT "Print[\"Smith form of W:\"];\n";
	    print OUT "rows:=Dimensions[W][[1]];\n";
	    print OUT "columns:=Dimensions[W][[2]];\n";
	    print OUT "Print[\"W has \",rows,\" rows and \",columns \"columns\"];\n";
	    print OUT "Print[MatrixForm[QW]];\n";
	    print OUT "Print[MatrixForm[DW]];\n";
	    print OUT "dimension:=0;\n";
	    print OUT "K:=NullSpace[Transpose[DW]];\n";
	    print OUT "If[Length[K]>0,\n";
	    print OUT "Stable:=K.PW;\n";
	    print OUT "Print[\"Basis of STABLE in IRREDUCIBLE coordinates:\"];\n";
	    print OUT "Print[\"Each row gives a stable sum of irreducible representations\"];\n";
	    print OUT "Print[Parameters];\n";
	    print OUT "Print[\"StableSums=\"];\n";
	    print OUT "Print[MatrixForm[Stable]];\n";
	    print OUT "dimension:=MatrixRank[Stable];\n];\n},\n{\n";
	    print OUT "Print[\"No characters vanish near 0; everything is stable\"];\n";
	    print OUT "dimension:=Length[Parameters];\n";
	    print OUT "}];\n";
	    print OUT "Print[\"Dimension of space of stable characters: \",dimension];\n";
	    close(OUT); 
	}else{
	    $v>1 and print "Using default mathematica NullSpace\n";
	    print OUT "K:=NullSpace[Transpose[Pairing]];\n";
	    print OUT "If[Length[K]>0,\n{\n";
	    print OUT "K=Transpose[K];\n";
	    print OUT "rows:=Dimensions[K][[1]];\n";
	    print OUT "columns:=Dimensions[K][[2]];\n";
	    print OUT "Print[\"K has \",rows,\" rows and \",columns \"columns\"];\n";
	    print OUT "W:=ParametersAmongDualParameters.K;\n";
	    print OUT "Stable:=NullSpace[Transpose[W]];\n";
	    print OUT "Print[\"Stable:\",MatrixForm[Stable]];\n";
	    print OUT "Print[\"length W: \",Length[W]];\n";
	    print OUT "Stable:=NullSpace[Transpose[W]];\n";
	    print OUT "Print[\"Basis of STABLE in IRREDUCIBLE coordinates:\"];\n";
	    print OUT "Print[\"Each row gives a stable sum of irreducible representations\"];\n";
	    print OUT "Print[Parameters];\n";
	    print OUT "Print[\"StableSums=\"];\n";
	    print OUT "Print[MatrixForm[Stable]];\n";
	    print OUT "dimension:=MatrixRank[Stable];\n},\n{\n";
	    print OUT "Print[\"No characters vanish near 0; everything is stable\"];\n";
	    print OUT "dimension:=Length[Parameters];\n";
	    print OUT "}];\n";
	    print OUT "Print[\"Dimension of space of stable characters: \",dimension];\n";
	    close(OUT); 
	}

	my $mathematicaOutput=`math < $mathematicaFile`;
#	$mathematicaOutput =~ s/.*Smith/Smith/s;
	$v>2 and print "OUTPUT:\n", $mathematicaOutput;
	if ($mathematicaOutput =~ /No characters vanish near 0/s){
	    ($junk,$dimension) = ($mathematicaOutput =~ /(.*Dimension of space of stable characters:\s*)([0-9]*)/s);
	    $everythingIsStable=1;
	}else{
	    if ($exact){
		($junk,$sizeW,$basis,$dimension)= ($mathematicaOutput =~ /.*Smith form of W:(.*).*(W has.*columns).*Basis.*StableSums=\s*(.*)Dimension of space of stable characters:\s*(.*)\s*/s);
	    }else{
		($junk,$sizeW,$basis,$dimension)= ($mathematicaOutput =~ /(.*).*(K has.*columns).*Basis.*StableSums=\s*(.*)Dimension of space of stable characters:\s*(.*)\s*/s);
	    }
	}
	$dimension=~ s/\n.*//s;
	chomp($dimension);
	$basis =~ s/In\[[0-9]*\]:=//gs;
	$basis =~ s/[^0-9]*$//gs;
	$dimensions{$cellsString}=$dimension if ($cellsString and $dimension);
	$dimension |=0;
    }else{
	$dimension=0;
    }
    (scalar(@$subsetSimpleRoots)) and print "lambda is singular at simple roots: ", join ',', @$subsetSimpleRoots;
    $cells and print "\ncells:", join ',', $cells;
}
    print "Not computing stability\n" if ($noStable);
    print "Parameters (living at lambda): ", join ',',@$parameters;
    print "\n";
    foreach my $p (@$parameters){
	print "$blockData->{$p}[4]";
    }

    print "\nDual parameters (to those living at lambda): ";
    my @dualParameters;
    foreach my $p (@$parameters){
	if ($dual){
	    push @dualParameters, $dualMap->{$p};
	}else{
	    push @dualParameters, $map->{$p};
	}
    }
    print join ',', @dualParameters;
    print "\n";

    if ($v>0){
	my $dualBlockFile="dualblock";
	my ($dualBlockData,$allDualParameters)=parseBlock($dualBlockFile);
	foreach my $p (@dualParameters){
	    print "$dualBlockData->{$p}[4]";
	}
    }
    unless ($noStable){
	print "\nDimension of space of stable characters: $dimension";
	print "\n";
	if ($everythingIsStable){
	    print "Everything is stable\n";
	}else{
	    ($dimension > 0) and print "Basis of stable characters expressed as sums of  irreducibles ";
	    print join ',', @$parameters;
	    print ":\n$basis\n";
	}
    }
}


sub parametersToMatrix{
    $v>3 and print "sub parametersToMatrix\n";
    my ($parameters,$n)=@_;
    my @parameters=@$parameters;
#    print "p:", join ',', @parameters;
#    print "n:$n\n";
    my $rv=zeroes($n,scalar(@parameters));
    foreach my $i (0..scalar(@parameters)-1){
	my $parameter=$parameters[$i];
#	print "doing $i,$parameter\n";
	$rv->slice("$parameter,$i") .= 1;
    }
    $v>2 and print "parametersToMatrix:$rv";
    return $rv;
}
sub outputPerlMatrix{
    $v>3 and print "sub outputPerlMatrix\n";
    my $m=shift;
    my ($c,$r)=$m->dims;
    my @M;
    foreach my $i (0..$r-1){
	my @row;
	foreach my $j (0..$c-1){
	    push @row, $m->at($j,$i);
	}
	my $row=join ',', @row;
	push @M, "[$row]";
    }
    my $M="[";
    $M.= join ",\n", @M;
    $M.="]";
#    print "M is: $M\n";
    return $M;
}
sub mathematicaMatrix{
    $v>3 and print "sub mathematicaMatrix\n";
    my $m=shift;
    my ($c,$r)=$m->dims;
    my @M;
    foreach my $i (0..$r-1){
	my @row;
	foreach my $j (0..$c-1){
	    push @row, $m->at($j,$i);
	}
	my $row=join ',', @row;
	push @M, "{$row}";
    }
    my $M="{";
    $M.= join ",\n", @M;
    $M.="}";
#    print "M is: $M\n";
    return $M;
}
	    
sub parametersContainingSInTau{
    $v>3 and print "sub parametersContainingSInTau\n";
    my ($allParameters,$subsetSimpleRoots,$tau)=@_;
    my @S=@$subsetSimpleRoots;
    my @dualParameters;
    foreach my $p (@$allParameters){
#	print "\nTESTING $p";
	my $t=$tau->{$p}; #list of roots for parameter $p
#	print "\nTesting ", Dumper($t);
	my $include=1;      #parameter $p contains all roots in its tau invariant
	foreach my $i (@S){
#	    print "\ntesting $i", $t->[$i];
	    if ($t->[$i]==0){
#		print "\nDIES";
		$include=0; #root in tau => parameter dies
		last;
	    }
	}
	push @dualParameters, $p if ($include);
    }
    @dualParameters=sort {$a <=> $b} @dualParameters;
    print "\n dual parameters (", join ',', @S, "):  ", join ',', @dualParameters;
    print "\n\n";
    return \@dualParameters;
}


sub pushToSingular{
    $v>3 and print "sub pushToSingular\n";
    #simple roots labelled 1,,,n (not 0.. n-1)
    my ($parameters,$subsetSimpleRoots,$tau)=@_;
    my @S=@$subsetSimpleRoots;
    my @singularParameters;
    foreach my $p (@$parameters){
	my $t=$tau->{$p}; #list of roots for parameter $p
	my $lives=1;      #parameter $p lives at lambda
	foreach my $i (@S){
	    if ($t->[$i-1]==1){   #$i labelled 1..n not 0..n-1
		$lives=0; #root in tau => parameter dies
		last;
	    }
	}
	push @singularParameters, $p if ($lives);
    }
    @singularParameters=sort {$a <=> $b} @singularParameters;
    $v>2 and print "\nlives at singular (", join ',', @S, "):  ", join ',', @singularParameters;
    $v>2 and print "\n\n";
    return \@singularParameters;
}


sub computeTau{
    $v>3 and print "sub computeTau\n";
    my ($parameters,$blockData)=@_;
#    $v>1 and print Dumper($blockData);exit;
    my $rootsInTau='ic|C-|r1|r2';
    my %tau;
    foreach my $i (@$parameters){
	my $roots=$blockData->{$i}[3];
	my @roots=split ',', $roots;
	my  @match;
	foreach my $root (@roots){
	    my $match= ($root =~ /$rootsInTau/);
	    $match|=0;
	    push @match,$match;
	}
	$tau{$i}=\@match;	
    }
#    print Dumper(\%tau);
    return \%tau;
}



sub selectRows{
    $v>3 and print "sub selectRows\n";
    my $parameters=shift;
    my $k=shift;
    my $rv;
    return $k unless ($parameters);
    my ($c,$r)=$k->dims; 
    my $rv=$k->dice([0..$c-1],$parameters);
    return $rv;
}






sub parseBlock{
    my $blockFile=shift;
    $v>3 and print "sub parseBlock\n";
    my %data;
#   7( 7,2):  2  1  [i2,C-]   7   6   (10,11)  ( *, *)   2,1,2
    open(IN,"<$blockFile")||die("Can't open $blockFile for input");
    foreach my $line (<IN>){
	my $lineOrig=$line;
	my ($number,$x,$y,$lengthAndCartan,$roots)=
	    ($line =~ /^\s*([0-9]*)\(\s*([0-9]*)\s*,([0-9\s]*)\):\s*(.*)\[(.*)\].*/g);
	$lengthAndCartan=~ s/\s*$//;
	my ($length,$cartan)=split(/\s*/,$lengthAndCartan);
#	print "line:$line\n $number :: $x :: $y :: $lengthAndCartan :: $roots\n\n";exit;
	$data{$number}=[$x,$y,$length,$roots,$lineOrig];
    }
    #ASSUMING two parameters with same y value are consecutive
#    print "data:", Dumper(\%data);
    my @allParameters=keys(%data);
    @allParameters = sort {$a <=> $b} @allParameters;
    return (\%data,\@allParameters);
}

sub pairingStableDualParameters{
    $v>3 and print "sub vanishNearZero\n";
#input: blockData, dualParameters, dualKlmatrix
#X(S)=<irreducibles on dual side given by dualParameters>
#kl basis matrix on dual side is (klbasis_signed)transpose-inverse)_signed

#return: matrix whose rows give sum of standard modules    
    my ($data,$dualParameters,$dualKLMatrix)=@_;
#    print "dualp:", Dumper($dualParameters);
    my $sizeBlock=scalar(keys %$data);
    my $dualParametersMatrix=parametersToMatrix($dualParameters,$sizeBlock);
    $v>2 and print "dual parameters matrix: $dualParametersMatrix\n";
    my $dualParametersStandardBasis=$dualParametersMatrix x $dualKLMatrix;
    my @allParameters=keys(%$data);
    @allParameters = sort {$a<=> $b} @allParameters;
#    print Dumper($data);exit;
    my @dualParameters=@$dualParameters;
#    print "working on dual parameters:", join ',', @dualParameters;
    
    my %yValues;
#    print  Dumper($data);
    foreach my $param (0..$sizeBlock-1){
	my $y=$data->{$param}->[1];
	push @{$yValues{$y}}, $param;
    }
#     foreach my $y (keys %yValues){
# 	my $parameters=$yValues{$y};
# 	my @parameters=@$parameters;
# 	delete($yValues{$y}) unless (scalar(@parameters)>1);
#     }


    my $rows=scalar(keys %yValues);
    my $stableAtRho=zeroes($sizeBlock,$rows);
    my @keys=keys(%yValues);
    foreach my $i (0..$#keys){
	my $key=$keys[$i];
#	print "doing $i,$key\n";
	my $parameters=$yValues{$key};
	my @parameters=@$parameters;
#	print "Params:", join ',', @parameters;
	foreach my $param (@parameters){
	    $stableAtRho->slice("$param,$i") .= 1;
	}
    }
    $stableAtRho=transpose($stableAtRho);
    $v>2 and print "\nstable at rho:", $stableAtRho;
    $v>2 and print "\ndual parameters matrix:",$dualParametersStandardBasis;
    $v>2 and print "Dual Parameters Standard Basis:", $dualParametersStandardBasis;
    my $vanishing=$dualParametersStandardBasis x $stableAtRho;
    $v>2 and print "vanishing:$vanishing\n";
    return ($vanishing);
}

sub readFile{
    $v>3 and print "sub readFile\n";
    my $file=shift;
    open(IN,"<$file")||die("Can't open file $file for input\n");;
    my $data;
    foreach my $line (<IN>){
	next if $line =~ /^\#/;
	$data .= $line;
    }
    close(IN);
    return $data;
}

sub klMatrix{
    $v>3 and print "sub klMatrix\n";
#return signed KL matrix KL from output of klbasis
#in dual mode:
# get from dualKLMatrix-permutation-perl if available
# else get from dualKLMatrix-perl if available and write dualKLMatrix-permutation-perl 
# else get from dualKLMatrix and write dualKLMatrix-perl, dualKLMatrix-permutation-perl
# else compute it and write dualKLMatrix,dualKLMatrix-perl,dualKLMatrix-permutation-perl   
#in regular mode:
# get from dualKLMatrix-perl if availble
# else get from dualKLMatrix and write dualKLMatrix-perl
# else compute it and write dualKLMatrix,dualKLMatrix-perl
    my $dualKLMatrix;
    if ($dual){
	$v>1 and print "Dual mode\n";
	if (-e 'dualKLMatrix-permutation-perl'){
	    #read matrix
	    $v>1 and print "Reading dual KL matrix from dualKLMatrix-permutation-perl\n";
	    my $matrix=readFile('dualKLMatrix-permutation-perl');
	    $dualKLMatrix=pdl(eval($matrix));
#	    $dualKLMatrix=transpose($dualKLMatrix);
	    $v>1 and print "Done reading dualKLMatrix-permutation-perl\n";
	    $v>3 and print "klMatrix from file dualKLMatrix-perl:", $dualKLMatrix,"\n";
	    return $dualKLMatrix;	
	}elsif (-e 'dualKLMatrix-perl'){
	    #read matrix, write dualKLMatrix-perl-permutation
	    $v>1 and print "Reading dual KL matrix from file dualKLMatrix-perl\n";
	    my $matrix=readFile('dualKLMatrix-perl');
	    $dualKLMatrix=pdl(eval($matrix));
#	    $dualKLMatrix=transpose($dualKLMatrix);
	    $v>1 and print "Done reading dualKLMatrix-perl\n";
	    $v>3 and print "klMatrix from file dualKLMatrix-perl:", $dualKLMatrix,"\n";
	    $v>1 and print "computing permutation\n";

	    my ($map,$dualMap,$permutation)=parseDualMap();
	    my $inversePermutation=transpose($permutation);
	    $dualKLMatrix=transpose($permutation x $dualKLMatrix x $inversePermutation);
#	    $dualKLMatrix=$permutation x $dualKLMatrix x $inversePermutation;
	    my $outputMatrix=outputPerlMatrix($dualKLMatrix);
	    open(OUT,">dualKLMatrix-permutation-perl")||die("Can't open file dualKLMatrix-permutation-perl for output");;
	    print OUT $outputMatrix;
	    close(OUT);
	    return $dualKLMatrix;
	}else{
	    #read mathematica matrix, write dualKLMatrix-permutation-perrl
	    $v>1 and print "Reading dual KL matrix from file dualKLMatrix (mathematica)\n";
	    my $data=shift;
#	    $dualKLMatrix=`$klCommand $klBasisFile -s -f maple`;
#	    $dualKLMatrix =~ s/.*klmatrix_signed:=//s;

	    $dualKLMatrix=`$klCommand $klBasisFile  -f maple`;
	    $dualKLMatrix =~ s/.*klmatrix:=//s;

	    $dualKLMatrix=eval($dualKLMatrix);
	    $dualKLMatrix=pdl($dualKLMatrix);
	    my ($map,$dualMap,$permutation)=parseDualMap();
	    my $inversePermutation=transpose($permutation);
	    $dualKLMatrix=transpose($permutation x $dualKLMatrix x $inversePermutation);
#	    $dualKLMatrix=$permutation x $dualKLMatrix x $inversePermutation;
	    my $outputMatrix=outputPerlMatrix($dualKLMatrix);
	    open(OUT,">dualKLMatrix-permutation-perl")||die("Can't open file dualKLMatrix-permutation-perl for output");;
	    print OUT $outputMatrix;
	    close(OUT);
	}
    }else{
	#not dual mode
	if (-e 'dualKLMatrix-perl'){
	    #read perl matrix; don't write dualKLMatrix-perl-permutation since not in dual mode
	    $v>1 and print "Reading dual KL matrix from file dualKLMatrix-perl\n";
	    my $matrix=readFile('dualKLMatrix-perl');
	    $dualKLMatrix=pdl(eval($matrix));
#	    $dualKLMatrix=transpose($dualKLMatrix);
	}elsif (-e 'dualKLMatrix'){
	    #read perl matrix; don't write dualKLMatrix-perl-permutation since not in dual mode
	    $v>1 and print "Reading dual KL matrix from file dualKLMatrix (mathematica)\n";
	    my $matrix=readFile('dualKLMatrix');
	    $dualKLMatrix=perlMatrix($matrix);
#	    $dualKLMatrix=transpose($dualKLMatrix);
	    $v>1 and print "Done reading dualKLMatrix\n";
	    $v>3 and print "klMatrix from file dualKLMatrix:", $dualKLMatrix,"\n";
	    my $outputMatrix=outputPerlMatrix($dualKLMatrix);
	    open(OUT,">dualKLMatrix-perl")||die("Can't open file dualKLMatrix-permutation-perl for output");;
	    print OUT $outputMatrix;
	    close(OUT);
	}else{
	    $v>1 and print "No dual matrix found; computing it\n using mathematica to compute inverse\n";
	    my $data=shift;
	    my $klMatrix=`$klCommand $klBasisFile -s -f maple`;
	    $klMatrix =~ s/.*klmatrix_signed:=//s;
	    $klMatrix=eval($klMatrix);
	    $klMatrix=pdl($klMatrix);
	    $v>3 and print "not dual mode; printing klMatrix to mathematica file to compute inverse\n";
	    open(OUT,">computeDualKLMatrix")||die("Can't open computeDualKLMatrix for output\n");
	    print OUT "(*temporary mathematica file, used only for computing dualKLMatrix*)\n";
	    print OUT "klMatrix:=";
	    print OUT mathematicaMatrix($klMatrix);
	    print OUT ";\n";
	    print OUT "Print[Inverse[klMatrix]];\n";
	    close(OUT);
	    $v>3 and print "computing inverse dualKLMatrix using mathematica\n";
	    my $output=`math < computeDualKLMatrix`;
	    $dualKLMatrix=perlMatrix($output);
	    $v>3 and print "Done computing inverse dualKLmatrix\n";
	    $v> 2 and print "klmatrix:$klMatrix\n";
	    $v>2 and print "dual klmatrix: $dualKLMatrix\n";
	    my ($r,$c)=$dualKLMatrix->dims;
	    my $signs;
	    my ($c,$r)=$dualKLMatrix->dims;
	    my $signMatrix=zeroes($c,$c);
	    foreach my $i (0..$c-1){
		my $length=$data->{$i}->[2];
		my $sign=(-1)**$length;
		$signMatrix->slice("$i,$i") .= $sign;
	    }
#	my $dualKlMatrix=$signMatrix x $dualKLMatrixInverseTranspose x $signMatrix;
#	    $dualKLMatrix=$dualKLMatrixInverse;

#	    $dualKLMatrix=$signMatrix x $dualKLMatrix x $signMatrix;

	    $v>2 and print "klmatrix:$dualKLMatrix\n";
	    $v>2 and print "dual klmatrix:$dualKLMatrix\n";
	    $v>3 and print "Done computing klmatrix\n";
	    
	    open(OUT,">dualKLMatrix")||die("Can't open dualKLMatrix for output");
	    print OUT "dualKLMatrix:=";
	    print OUT mathematicaMatrix($dualKLMatrix);
	    print OUT ";\n";
	    close(OUT);

	    my $outputMatrix=outputPerlMatrix($dualKLMatrix);
	    open(OUT,">dualKLMatrix-perl")||die("Can't open file dualKLMatrix-permutation-perl for output");;
	    print OUT $outputMatrix;
	    close(OUT);
	}
    }
    return $dualKLMatrix;
}

sub perlMatrix{
    $v>3 and print "sub perlMatrix\n";
    my $matrix=shift;
    $matrix =~ s/.*=\s*\{/\{/s;
    $matrix =~ s/\}\}.*/\}\}/s;
    $matrix =~ s/\{/[/sg;
    $matrix =~ s/\}/]/sg;
    $matrix =~ s/\n//g;
    $matrix =~ s/>//g;
    $matrix =~ s/\s//g;
    my $perlMatrix=pdl(eval($matrix));
    return $perlMatrix;
}

sub roundToInteger{
    $v>3 and print "sub roundToInteger\n";
    my $a=shift;
    return int($a+.5*($a<=>0));
}

sub matrixString{
    $v>3 and print "sub matrixString\n";
    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='}';
    $string.="$open$open";
    $string.=join "$close,\n$open", @output;
    $string.="$close$close";
    return $string;
}    

sub pad{
    $v>3 and print "sub pad\n";
    my $x=shift;
    my ($c,$r)=$x->dims;    
    print "c,r:$c,$r\n";
    return $x if ($c==$r);
    return transpose(pad(transpose($x))) if ($c>$r);
    my $z=zeroes($r-$c,$r);
    my $y=append($x,$z);
    return $y;
}

sub parseWCells{
    $v>3 and print "sub parseWCells\n";
    my $file=shift;
    open(IN,"<$file")||die("Can't open $file for input\n");
    my $lines;
    foreach my $line (<IN>){
	$lines .= $line;
    }

    $lines =~ s/.*Individual cells\.\n\/\/ //s;
    close(IN);
    my @cells = split 'cell \#', $lines;
    my %cells;
#    print "cells", join ":::::::::::::::::::::::\n", @cells;exit;
    foreach my $cell (@cells){
	my @parameters;
	(my $number)=  $cell =~ /^([0-9]*):/;
	my @lines = split "\n", $cell;
	foreach my $line (@lines){
	    next unless ($line =~ /\[/);
#	    print "doing $line\n";
	    (my $p) = $line =~ /[^\[]*\[([0-9]*)\].*/;
#	    print "got $number,$p\n";
	    push @parameters, $p;	    
	}
	$cells{$number}=\@parameters;
    }
#    print "CELLS:",Dumper(\%cells);
    return \%cells;
}
	

sub moreHelp{
print "
Algorithm:
";
}

sub aq{
    open(IN,"<$wCellsFile")||die("Can't open $wCellsFile for input");;
    my %cells;
    my %aqcells;
    foreach my $line (<IN>){
	chomp $line;
	last if ($line  =~ /graph/);
	next unless ($line =~ /^\#/);
	my ($cell,$reps)=split '=', $line;
	$cell =~ s/\#//;
	$reps =~ s/[\{\}]//g;
	my @reps=split ',', $reps;
	foreach my $rep (@reps){
	    $cells{$rep}=$cell;
	}
	$aqcells{$cell}=[];
    }
    close(IN);
    open(IN,"<$blockuFile")||die("Can't open $blockuFile for input");
    foreach my $line (<IN>){
	chomp $line;
	$line =~ s/\(.*//;
	$line =~ s/\s//g;
	push @{$aqcells{$cells{$line}}}, $line;
    }	
    print "A(lambda) cells computed using $blockuFile and $wCellsFile\n";
    print "cell\taq";
    foreach my $key (sort {$a <=> $b} keys %aqcells){
	my $reps=$aqcells{$key};
	my @reps=sort {$a<=>$b} @$reps;
	print "\n$key:\t", join ',', @reps;
    }
    print "\n\n";
    exit;
}

    
sub help{
print "
Warning: early version, poorly documented. This help 
file will be updated.

standard: wcells, block, klbasis
will produce dualKLMatrix file by:
 kl -s klbasis -> inverse -> dualKLMatrix (signs don't matter)
(klMatrix file doesn't matter)

dual mode:

run atlas for dual group, save block -> dualblock
                               wcells -> dualwcells
                               klbasis -> dualklbasis
                               dualblock -> block
                               dualmap -> dualmap

optional: dualKLMatrix file
(obtained from dualklbasis using kl,lower triangular, signs don't matter?)
internally when using it, modified by permutation coming from dualmap
stable -R -c 1,2,3 -S 2,4 -d
cells are cells for dualblock
S is roots for dual block, for C2 switch 1,2

Compute stable representations.

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

# 0( 0,6):   1   2    ( 6, *)  ( 4, *)    [i1,i1]  0  
# 1( 1,6):   0   3    ( 6, *)  ( 5, *)    [i1,i1]  0  
# 2( 2,6):   2   0    ( *, *)  ( 4, *)    [ic,i1]  0  
# 3( 3,6):   3   1    ( *, *)  ( 5, *)    [ic,i1]  0  
# 4( 4,4):   8   4    ( *, *)  ( *, *)    [C+,r1]  1  2
# 5( 5,4):   9   5    ( *, *)  ( *, *)    [C+,r1]  1  2
# 6( 6,5):   6   7    ( *, *)  ( *, *)    [r1,C+]  1  1
# 7( 7,2):   7   6    (10,11)  ( *, *)    [i2,C-]  2  212
# 8( 8,3):   4   9    ( *, *)  (10, *)    [C-,i1]  2  121
# 9( 9,3):   5   8    ( *, *)  (10, *)    [C-,i1]  2  121
#10(10,0):  11  10    ( *, *)  ( *, *)    [r2,r1]  3  1212
#11(10,1):  10  11    ( *, *)  ( *, *)    [r2,rn]  3  1212