#!/usr/bin/perl -w
#last updated 6/30/08

#bigisos.pl
#differs from isos.pl in computing valid actions at run time rather than compile time
#works on bigger files; is slightly slower than isos.pl
#is designed to run without permutations, but works with them.

#lists all the isomorphisms between two given blocks, either allowing or disallowing permutations of the simple roots
#representations are given "coordinate-free" names based on paths from user-chosen calibrators
#output consists of one column per isomorphism, headed by the required root-permutation (if that option is chosen)
#output also lists paths for each representation based on chosen calibrators
#send bugs to nnn9245@gmail.com

#get first piece of block data from the user and store it
&defineVariables(1);
&getDesiredInfo;

#describe representations in terms of cross-actions and Cayley transforms from calibrators
@pathStringsForwards = &pathfinder($source, 1);
@pathStringsBackwards = &pathfinder($sink, 0);
$psf = \@pathStringsForwards;
$forwardsPathDisplayLength = &maxStringLength($psf) + 4; #just used later for display

#get second piece of block data from user
&defineVariables(0);

#there's no explicit test for surjectivity, so allowing the algorithm to run when cardinalities are unequal might make false positives
unless($#pathStringsForwards == $highestRep){
	print "These blocks aren't isomorphic; they don't have the same cardinality.\n";
}

#set up the main routine by creating allowed permutations.
else{
$isoCount = 0;

print "Allow root permutations? Y/N: ";
$permsAllowed = <STDIN>;
$permsAllowedBoolean = ($permsAllowed =~ /[Yy]/);
if ($permsAllowedBoolean){
	@permutations = listPermutations($rootTotal - 1);
} else{ #trivial permutation only
	$placeholder = "placeholder";
	push @permutations, $placeholder;
}

#for each possible permutation of the roots, for each pair of possible calibrators on the second block,
#we see if we can build up the second block using the path data from the first.
foreach $perm(@permutations){
	if ($permsAllowedBoolean){
		$blockMatrixRef = \@blockMatrix;
		@blockMatrix = &permuteMatrixRoots($perm, $blockMatrixRef);
	}
	@sourceList = &listEligibleSources;
	@sinkList = &listEligibleSinks;
	
	foreach $potentialSource(@sourceList){
		foreach $potentialSink(@sinkList){
			for($i = 0 ; $i < $M ; $i++){
				$failed = 0;
				$myForwardsPath = $pathStringsForwards[$i];
				$myBackwardsPath = $pathStringsBackwards[$i];
				
				@pathArray = &pathArray($myForwardsPath);  
				$pathRef = \@pathArray;
				@possiblePointsForward = &readAPath($potentialSource, $pathRef, 1);
				if (@possiblePointsForward){
					$refForward = \@possiblePointsForward;
				} else{$failed = "true"; last;} #give up if the forwards path doesn't define any points
				
				@pathArray = &pathArray($myBackwardsPath);
				$pathRef = \@pathArray;
				@possiblePointsBackward = &readAPath($potentialSink, $pathRef, 0);
				if (@possiblePointsBackward){
					$refBackward = \@possiblePointsBackward;
				} else{$failed = "true"; last;} #give up if the backwards path doesn't define any points
				
				@potentialUniquelyDeterminedPoint = &intersection($refForward, $refBackward);
				if($#potentialUniquelyDeterminedPoint){$failed = "true"; last; #give up if there's either no or >= two points in the intersection).	
				}else{
					$isomorphismArray[$i] = $potentialUniquelyDeterminedPoint[0];	
				}
			}
			
			my @lengthTester = &purgeDuplicates(@isomorphismArray);
			unless($#lengthTester = $#isomorphismArray){
				$failed = "true"; #this tests injectivity	
			}
			
			unless($failed){
			my @storeIt = @isomorphismArray;
			$tempRef = \@storeIt;
			push @isoHolder, $tempRef;
			push @permHolder, $perm;
			$isoCount++;
			}
			@isomorphismArray = ();
		}
	}
	#reset the block matrix if perms are allowed-- this is inefficient, making 2n! copies of blockMatrix where a smarter code could make n!
	#since the point of isos2.pl is for perms to be disallowed, i am not going to worry about this for now.
	if($permsAllowedBoolean){	
		@inverse = &invertPerm($perm);
		$inverseRef = \@inverse;
		@blockMatrix = permuteMatrixRoots($inverseRef, $blockMatrixRef);
	}
}

#display results to the user.
print "\nTotal isomorphism count is $isoCount.\n\n";
$compactSpaceString = &fixedLengthDisplayer("CPath", $forwardsPathDisplayLength);
if($permsAllowedBoolean){
	for(my $j = 0; $j < $isoCount; $j++){
		@perm = @{$permHolder[$j]};
		$permString = "[@perm] ";
		print $permString;
	}
	print "$compactSpaceString SPath \n\n";
	$maxLength = length($permString);
	
} else{$maxLength = length($highestRep) + 3;
	$skipThis = $maxLength * $isoCount;
	$blankSpaceString = &fixedLengthDisplayer("", $skipThis);	
	print "$blankSpaceString $compactSpaceString SPath \n\n";
}
for(my $i = 0; $i < $M; $i++){
	for(my $j = 0; $j < $isoCount; $j++){
	my $string = $isoHolder[$j][$i];
	$string = &fixedLengthDisplayer ($string, $maxLength);
	print "$string";
	}
	my $pathStringForwards = &fixedLengthDisplayer ($pathStringsForwards[$i], $forwardsPathDisplayLength);
	print "$pathStringForwards $pathStringsBackwards[$i]";
	print "\n";
}
} #this bracket closes an annoying "else" above


#gets block file from user and stores it as a matrix.
#THIS MATRIX IS STORED AS THE TRANSPOSE OF THE MATRIX ONE "SEES" WHEN LOOKING AT A BLOCK FILE
#this is because swapping rows of a matrix is easy in perl and swapping columns is hard.
sub defineVariables{
	$highestRep = -1;
	if($_[0]){
		print "Enter the first block (or blockd) output file: \n";
	} else{
	print "Enter the second block (or blockd) output file: \n";
	}
	my $textfile = <STDIN>;
	chomp($textfile);
	
	until(open(FILE, "$textfile")){
		print "The file you have tried to open does not exist. ";
		if($_[0]){
			print "Enter the first block (or blockd) output file: \n";
		} else{
			print "Enter the second block (or blockd) output file: \n";
		}
	$textfile = <STDIN>;
	chomp($textfile);
	}
	while (<FILE>){
		my $lineout = $_;
		if ($lineout =~ /(\s*\d+)(\(\s*\d+,\s*\d+\)):(\s+\d+)(\s+\d+\s+)\[([a-zA-Z_0-9_,_\+_-]+)\](\s+)([\s_\d]+)((\([\w_\W]+\)\s*)+)([\w_\W]*$)/){
	   
			$displayRep = $1; #for redisplaying later
			$rep = &trim($displayRep);
	   
			#$kgbpair = $2;   this is being captured for a possible future program
	   
			$displayLength = $3; #don't forget the colon
			$length = &trim($displayLength);
			
			$displayCartan = $4;
			$cartan = &trim($displayCartan);
	   
	   
			@roots = split /,/, $5;
			unless($rep){
				$rootTotal = @roots; #keep from having to compute several values too many times
				$shiftedRootTotal = $rootTotal + 2;
				$doubleRootTotal = $rootTotal + $rootTotal; 
				$shiftedDoubleRootTotal = $doubleRootTotal + 2;
				$tripleRootTotal = $doubleRootTotal + $rootTotal;
				$shiftedTripleRootTotal = $tripleRootTotal + 2;
				$quadrupleRootTotal = $tripleRootTotal + $rootTotal;
			}	
			#$junkSpace = $6; (for possible future program)
	   
			@crossActions = split /\s+/, $7;
	   
			@cayleys = split /\)\s\s\(/, $8;  #these will be cleaned up in a minute..
	   
			#$remainder = $10;	   #we don't want $9, which is embedded in $8 for technical reasons. $10 is in potential future program.
			
			#these lines clean up the Cayley transforms we captured (last set of parens w/o backlashes in reg exp) 
			if($cayleys[0] =~ /\(([\w_\W]+)\)*/){
				$cayleys[0] = $1;	 
			}
		 
			if($cayleys[$rootTotal - 1] =~ /([\w_\W]+)\)/){
				$cayleys[$rootTotal - 1] = $1;
			}

			#extract the desired Cayleys using a new array which will end up twice the size of the old one.
			@extractedCayleys = @cayleys;
			for (my $i = 0; $i<$rootTotal; $i++){
				@temporaryTwoItemList = split /,/,$cayleys[$i];   #make a list of the two numbers associated to the Cayley transform 
				foreach $item(@temporaryTwoItemList){   #extract remaining whitespace
					if($item =~ /\s*(\S+)\s*/){ 
						$item = $1;
					}
					$extractedCayleys[2*$i] = $temporaryTwoItemList[0];
					$extractedCayleys[2*$i + 1] = $temporaryTwoItemList[1];	
				}
			}
			
			#save the data to a large matrix
			$blockMatrix[0][$rep] = $length; 
			$blockMatrix[1][$rep] = $cartan;
			
			for($j = 0; $j < $rootTotal; $j++){		
				$blockMatrix[$j + 2][$rep] = $roots[$j];
				$blockMatrix[$j + $shiftedRootTotal][$rep] = $crossActions[$j];
				$blockMatrix[$j + $shiftedDoubleRootTotal][$rep] = $extractedCayleys[$j];
				$blockMatrix[$j + $shiftedTripleRootTotal][$rep] = $extractedCayleys[$j + $rootTotal];
			}
			
			$highestRep++;
		}

	}
	close(FILE);
	$mostCompactCartan = &myCartan(0);
	$mostSplitCartan = &myCartan($highestRep);
	$M = $highestRep + 1; #the matrix is M by N (in the CS sense).
	$N = 2 + $quadrupleRootTotal;
	$shiftN = $N - 1;
}

#a helper subroutine in parsing the data
sub trim(){ #from http://www.somacon.com/p114.php

	my $string = shift;
	$string =~ s/^\s+//;
	$string =~ s/\s+$//;
	return $string;
	
}


#this tells the user which reps are acceptable calibrators and gets reps from user. it remembers the roots associated to them.
sub getDesiredInfo{
	print "The representations are labeled 0 through $highestRep. Pick a calibrating representation on the most compact Cartan (c to list valid calibrators): ";
	chomp ($source = <STDIN>);
	
	if($source =~ /c/){
		&listCompactChoices;
		print "Pick a calibrating representation on the most compact Cartan: \n";
		chomp ($source = <STDIN>);
	}
	
	until(&myCartan($source) == $mostCompactCartan){
		print "Valid starting representations must lie on the most compact Cartan. Try again: "; #print what they are
		chomp($source = <STDIN>);
	}

	print "Now pick a calibrating representation on the most split Cartan (c to list valid calibrators): ";

	chomp ($sink = <STDIN>);
	
	if($sink =~ /c/){
		&listSplitChoices;
		print "Pick a calibrating representation on the most split Cartan: \n";
		chomp ($sink = <STDIN>);
	}

	until(&myCartan($sink) == $mostSplitCartan){
		print "Valid starting representations must lie on the most split Cartan. Try again: ";
		chomp($sink = <STDIN>);	
	}
	
	for($shiftedJ = 2; $shiftedJ < $shiftedRootTotal; $shiftedJ++){

	$sourceRootTypes[$shiftedJ] = $blockMatrix[$shiftedJ][$source];
	$sinkRootTypes[$shiftedJ] = $blockMatrix[$shiftedJ][$sink];
	}
	
	$sourceRootTypes[0] = $sourceRootTypes[1] = $sinkRootTypes[0] = $sinkRootTypes[1] = "there is an error if you are seeing this";
}

#these next two subroutines just print out the valid calibration choices.
sub listCompactChoices{
	print "Valid representations are: ";
	for($rep = 0; $rep <= $highestRep; $rep++){
		if(&myCartan($rep) == $mostCompactCartan){print "$rep "}
	}
	print "\n";
}

sub listSplitChoices{
	print "Valid representations are: ";
	for($rep = 0; $rep <= $highestRep; $rep++){
		if(&myCartan($rep) == $mostSplitCartan){print "$rep "}
	}
	print "\n";
}

#you give this an xa or Ct by its location in the block matrix and get what root it is associated to.
sub rootOfAction{
	my $rep = $_[0];
	my $actionIndex = $_[1];
	
	if ($actionIndex < $shiftedDoubleRootTotal){ #if it's a crossaction
		return $blockMatrix[$actionIndex - $rootTotal][$rep]; #can't use mod cause rootTotal might be <= 2
	} else{ #cayleys
		$actionPair = int(($actionIndex - $shiftedDoubleRootTotal)/2);
		return $blockMatrix[2 + $actionPair][$rep];
	}
}

#you give this a representaiton and get its associated Cartan
sub myCartan{ #mike artin?
	my $rep = $_[0];
	return $blockMatrix[1][$rep];
}

#returns a list of calibrators on the second block that could match the chosen compact one on the first block
sub listEligibleSources{
	my @eligibleSources;
	
	for($calibrator = 0; $calibrator < $M; $calibrator++){	
		#must be on the right cartan to be eligible
		$eligibleBoolean = 1;
		unless (&myCartan($calibrator) == $mostCompactCartan){
			$eligibleBoolean = 0;
		}
		

		#if on the right cartan, root types must line up with other root types.
		else{ 
			for($shiftedJ = 2; $shiftedJ < $shiftedRootTotal; $shiftedJ++){
				$eligibleBoolean = $eligibleBoolean*($sourceRootTypes[$shiftedJ] eq $blockMatrix[$shiftedJ][$calibrator]);
			}	
		}
	
		if ($eligibleBoolean){
			push @eligibleSources, $calibrator;	
		}
	}	
	
	return @eligibleSources;
}

#returns a list of calibrators on the second block that could match the chosen split one on the first block
sub listEligibleSinks{
	my @eligibleSinks;
		for($calibrator = 0; $calibrator < $M; $calibrator++){	
		$eligibleBoolean = 1;
		#must be on the right cartan to be eligible
		unless (&myCartan($calibrator) == $mostSplitCartan){
			$eligibleBoolean = 0;
		}
		
		#if on the right cartan, root types must line up with other root types
		else{ 
			$eligibleBoolean = 1;
			for($shiftedJ = 2; $shiftedJ < $shiftedRootTotal; $shiftedJ++){
				$eligibleBoolean = $eligibleBoolean*($sinkRootTypes[$shiftedJ] eq $blockMatrix[$shiftedJ][$calibrator]);
			}	
		}
	
		if ($eligibleBoolean){
			push @eligibleSinks, $calibrator;	
		}
	}	
	
	return @eligibleSinks;
}


#we don't allow forwards transformations that aren't C+ or i1/i2, except on the most compact Cartan
sub allowedForwards{
	my $rep = $_[0];
	my $axn = $_[1];
	if ($shiftedRootTotal <= $axn and $axn < $shiftedDoubleRootTotal){  #cross-actions
		return((&myCartan($rep) == $mostCompactCartan) or (&rootOfAction($rep, $axn) eq "C+"));
	}else{ #cayleys
		return(&rootOfAction($rep, $axn) =~ /i\d/);
	}
}

#we don't allow backwards transformations that aren't C- or r1/r2, except on the most split Cartan
sub allowedBackwards{
	my $rep = $_[0];
	my $axn = $_[1];
	if ($shiftedRootTotal <= $axn and $axn < $shiftedDoubleRootTotal){ #cross-actions
		return((&myCartan($rep) == $mostSplitCartan) or (&rootOfAction($rep, $axn) eq "C-"));
	}else{ #cayleys
		return(&rootOfAction($rep, $axn) =~ /r\d/);
	}
}

#pathfinder-- give it entrance, forwards=1 or backwards = 0
#keeps track of points we know how to get to (store) and newest points we just figured out how to get to (list)
#goes one step from each point in list, then adds this new path to store and changes list.
sub pathfinder{	
	
	my $entrance = $_[0];   #calibrator
	my $forwards = $_[1];   #is it the compact (true) or split (false) calibrator?
	my @store;   #stuff that i know how to get to
	my @list;   #stuff that i'm working with temporarily.
	
	$list[0] = $entrance;
	$store[$entrance] = "";
	
	while(@list){
		$rep = shift @list;
		for(my $j = $shiftedRootTotal; $j<=$shiftN; $j++){
			if(($forwards and &allowedForwards($rep, $j)) or (!$forwards and &allowedBackwards($rep,$j))){
				$neighborInQuestion = $blockMatrix[$j][$rep];
			}
			else{
				$neighborInQuestion = "*";
			}
			
			if($neighborInQuestion ne "*"){
				unless(defined($store[$neighborInQuestion])){
					if($j < $shiftedDoubleRootTotal){
						$root = $j - $shiftedRootTotal;
						$instruction =  "x$root";
					} else{
						$root = int(($j - $shiftedDoubleRootTotal)/2);	
						$instruction = "^$root";
					}
					$store[$neighborInQuestion] = "$instruction $store[$rep]";
					push @list, $neighborInQuestion;
				}
			}
		}
	}
	return @store;
}

#turns a path that's a string into an array.
sub pathArray{
	my $string = $_[0];
	@pathArray = split /\s+/, $string;
	return @pathArray;
}

#subroutine takes starting point, a pathArray, and a boolean that asks if the path moves backwards or forwards. It returns a list of all potential endpoints.
#returns false if there's a deadend (which looks like a path going through * in the appropriate matrix).
sub readAPath{
	my $firstPoint = $_[0];
	my @pathToRead = @{$_[1]};
	my $forwards = $_[2];
	my @startingPoints;   #a list of current starting points. not just one item because we get more starting points as we use multi-valued cayleys.
	push @startingPoints, $firstPoint;
	while (@pathToRead){
		$myInstruction = pop @pathToRead;
		if ($myInstruction =~ /x(\d+)/){
			$myInstruction = $1;
			my %newPoints;
			foreach $startingPoint (@startingPoints){
				my $j = $myInstruction + $shiftedRootTotal;
				if (($forwards and &allowedForwards($startingPoint, $j)) or (!$forwards and &allowedBackwards($startingPoint,$j))){
					$newPoint = $blockMatrix[$j][$startingPoint];
				} else{
					$newPoint = "*";
				}
				unless ($newPoint eq "*"){
					$newPoints{$newPoint} = 1;
				}else{
					return;  #DEAD END
				}
			}
			@startingPoints = keys %newPoints;
		}elsif ($myInstruction =~ /\^(\d+)/){
			$myInstruction = $1;
			my %newPoints;
			foreach $startingPoint (@startingPoints){
				my $j2 = 2*$myInstruction +$shiftedDoubleRootTotal + 1;
				my $j1 = 2*$myInstruction + $shiftedDoubleRootTotal;
				
				if (($forwards and &allowedForwards($startingPoint, $j2)) or (!$forwards and &allowedBackwards($startingPoint,$j2))){
					$newPoint2 = $blockMatrix[$j2][$startingPoint];
				} else {$newPoint2 = "*";}
				
				if (($forwards and &allowedForwards($startingPoint, $j1)) or (!$forwards and &allowedBackwards($startingPoint,$j1))){
					$newPoint1 =  $blockMatrix[$j1][$startingPoint];
				} else{$newPoint1 = "*";}

				unless($newPoint2 eq "*"){
					$newPoints{$newPoint2}  = 1;
				}
				unless($newPoint1 eq "*"){
					$newPoints{$newPoint1} = 1;
				}else{ 
					return;  #DEAD END
				}
			}
			@startingPoints = keys %newPoints;
		}
	}	
	return @startingPoints;		
}

#given two arrays, thought of as unordered sets, finds their set-theoretic intersection.
sub intersection{
	my @list1 = @{$_[0]};
	my @list2 = @{$_[1]};
	my %hash1;
	my %hash2;
	
	my @finalAnswer;
	
	foreach $element (@list1){
		$hash1{$element} = 1;	
	}
	foreach $element (@list2){
		$hash2{$element} = 1;
	}

	foreach $key (keys %hash1){
		if(defined($hash2{$key})){
			push @finalAnswer, $key;		
		}
	}
	return @finalAnswer;
}

#lists S_n recursively.
sub listPermutations{
	my $n = $_[0];
	my @listOfReferencesToPermutations;
	if($n == 0){
		my @perm = (0);
		my $ref = \@perm;
		push @listOfReferencesToPermutations, $ref;
	}
 	else{
		my @permsFixingN = &listPermutations($n - 1);
		foreach $ref(@permsFixingN){
				my @perm = @{$ref};
				for($j = $#perm + 1; $j >= 0; $j--){
					my @newPerm = @perm;
				        splice(@newPerm, $j, 0, $n);
					my $newRef = \@newPerm;
					push @listOfReferencesToPermutations, $newRef;
				}
			}
		}	
	return @listOfReferencesToPermutations;
}

#returns the inverse of a permutation.
sub invertPerm{
	my @perm = @{$_[0]};
	for($i = 0; $i < @perm ; $i++){
		$newIndex = $perm[$i];
		$newPerm[$newIndex] = $i;
	}
	return @newPerm;
}

#adds the constant to each component of the array
sub addConstantToArray{
	my @array = @{$_[0]};
	my $constant = $_[1];
	foreach $thing(@array){
		$thing = $thing + $constant;	
	}
	return @array;
}

#turns (0, 3, 2, 1) into (0, 1, 6, 7, 4, 5, 1, 2), for example. (this is for moving Cayley transforms in pairs).
sub duplicatePermutation{
	my @perm = @{$_[0]};
	for(my $i = 0; $i <= $#perm; $i+=2){
		my @fixup = (2*$perm[$i], 2*$perm[$i] + 1);
		splice(@perm, $i, 1, @fixup);	
	}
	return @perm;
}

#swaps the block data around by permuting the roots and their associated x's and ^s.
sub permuteMatrixRoots{
	my @perm = @{$_[0]};
	my @matrix = @{$_[1]};
	$ref = \@perm;
	@duplicatePerm = &duplicatePermutation($ref);
	@rootPerm = &addConstantToArray($ref, 2);
	@crossPerm = &addConstantToArray($ref, $shiftedRootTotal);
	$duplicateRef = \@duplicatePerm;
	@cayleyPerm = &addConstantToArray($duplicateRef, $shiftedDoubleRootTotal);
	@finalPerm = (0, 1);
	push @finalPerm, @rootPerm;
	push @finalPerm, @crossPerm;
	push @finalPerm, @cayleyPerm;
	@matrix = @matrix[@finalPerm];
	return @matrix;
}

#adds spaces to left justify the given string until its length is a set number
sub fixedLengthDisplayer{
	my $string = $_[0];
	my $newLength = $_[1];
	if(length($string) <= $newLength){ 	
		until(length($string) == $newLength){
			$string = "$string ";	
		}
	}else{
		my $mistake = "Formatting error: need bigger width in call to subroutine fixedLengthDisplayer";
		return $mistake;
	}
	return $string;
}

#returns the maximum length of any string in an array.
sub maxStringLength{
	my @stringList = @{$_[0]};	
	my $maxLength = 0;
	foreach $string (@stringList){
		if(length($string) > $maxLength){
			$maxLength = length($string);	
		}
	}
	return $maxLength;
}


#takes an array and returns the same array with no duplicates.
sub purgeDuplicates{
	my @list = @_;
	my %unique;
	my @redundancies;
	
	for(my $i = 0; $i < @list; $i++){
		$element = $list[$i];
		if($unique{$element}){
			push @redundancies, $i;	
		} else{
			$unique{$element} = 1;
		}
	}
	for (my $j = $#redundancies; $j >= 0; $j--){
		$i = $redundancies[$j];
		splice @list, $i, 1;
	}
	return @list;
}