#!/usr/bin/perl

# Last edit: 3/14.  Added in computation of f-vector and h-vector.

# This program will read in data in the format proved by Gordon Royle.
# This program will output 

#   1) All information: bases, ind sets, rank, simple or not, connected or not
#   2) Bases input for latte or CDD
#   3) Topcom ready input


# How to read each data line:
#  - each matroid occupies one line of the file
#  - each line consists of a string of digits of length 2^e (then a newline character)
#  
#      01121223122223331223232323333333122323332333233322333333333333331223223323333333233333332333333323233333333333333333333333333333
#  
#      - the positions in the string are indexed from 0 .. 2^e-1
#      - the index i is to be viewed as a subset of {0,1,...e-1} according to its bit pattern
#      - the digit at index position i gives the rank of that subset
#  
#  
#      Hence: for this above string, we see that
#  
#      - the length is 128 and so we have a matroid of size 7
#      - the entry at position 127 = {0,1,2,3,4,5,6} is "3" and so the matroid has rank 3
#      - an index "x" represents an independent set if
#  
#      rank[x]  = #bits that are set in the binary representation for x

#this subroutine calculates factorial because I don't know what I'm doing.

# This takes in a binary number between 1 and 511 and converts it to a string
sub bintostring {
    if ($_[0] < 0) {
       return 0;
    }
    if ($_[0] > 512) {
       return 0;
    }

    $setstring = "";
    $value = $_[0];

    $one = 1;
    for ($j=0;$j<9;$j++){
        $testvalu2 = $one & $value;
        if (($one & $value) != 0) {
            $setstring = "$setstring" . "$j";
        }
        $one = $one*2;
    }
    
    $setstring;
}

sub subsettostring {
    if (length($_[0]) == 0) {
        return "{}";
    }
    else {
        $subsetlength = length($_[0]);
        $nicesubset = "{";
        for ($k=0;$k<$subsetlength;$k++) {
            $inputsubstring = substr($_[0],$k,1);
            $nicesubset = $nicesubset . "$inputsubstring,";
        }

    chop($nicesubset); 
    $nicesubset = $nicesubset . "}";
    }
}

# This sub procedure takes in a string $_[0] and removes the character
# at index $_[1], starting at 0
sub substrremi {
    # Do some error checking
    if (($_[1] < 0) || ($_[1] > length($_[0]))) {
        return $_[0];
    }

    substr($_[0],0,$_[1]) . substr($_[0],$_[1]+1,length($_[0] - $_[1] - 1));
}

# ********************************************
# ********************************************
#          Start sub royletoallinfo.pl
# ********************************************
# ********************************************
sub royletoallinfo {
    local $input;
    local $inputlength;
    local $largestindset;
    local $largestdepset;
    local $nicesubset;
    local $subset;
    local $allindsetsnice;
    local $allindsets;
    local $numindsets;
    local $sublength;
    local $largestdepset;
    local $alldepsetsnice;
    local $alldepsets;
    local $numdepsets;
    local $numelements;
    local $allCircuits;
    local $numCircuits;

    undef (@allindsetsnice);
    undef (@alldepsetsnice);
    undef (@allindsets);
    undef (@alldepsets);
    undef (@numindsets);
    undef (@numdepsets);
        
    $input = $_[0];
    chomp($input);

    $inputlength = length($input);

    if ($verbose == 1) { print ("******************** NEW MATROID $matroidCount ********************\n");}
    if ($verbose == 1) { print ("$input\n"); }
    if ($verbose == 1) { print ("inputlength = $inputlength\n"); }
    
    #starting off the f-vector
    @fvector = ();


    $largestindset = 0;
    $largestdepset = 0;
    for ($i=0; $i< $inputlength; $i++) {
        # Each index $i represents a subset of {0,1,...,e-1} by its
        # binary representation.
        $subset = bintostring($i);

        # We need to know how many elements are in our subset
        $sublength = length($subset);

        # The value at index $i tells its rank
        $subinput = substr($input,$i,1);
	
	# We also want a nice looking "{1,3,4}" representation of the subset
        $nicesubset = subsettostring($subset);
        if (length($subset) == $subinput) { # This means its an ind set
            if ($sublength > $largestindset) { # Keep track of the largest ind set
                $largestindset = $sublength;
            }
            $allindsetsnice[$sublength][$numindsets[$sublength]] = $nicesubset;
            $allindsets[$sublength][$numindsets[$sublength]] = $subset;
            $numindsets[$sublength]++;
        }
        else { # We know it is dependent
            if ($sublength > $largestdepset) { # Keep track of the largest dep set
                $largestdepset = $sublength;
            }
            #if ($verbose == 1) { print ("        $subset = $nicesubset is dependent.\n"); } 
            $alldepsetsnice[$sublength][$numdepsets[$sublength]] = $nicesubset;
            $alldepsets[$sublength][$numdepsets[$sublength]] = $subset;
            $numdepsets[$sublength]++;
        }
    }

    $numelements = 0;
    $tempinputlength = $inputlength;

    while ($tempinputlength > 1) {
        $numelements++;
        $tempinputlength = $tempinputlength / 2;
    }

    if ($verbose == 1) { print ("Number of elements: $numelements\n"); }


    if ($verbose == 1) { print ("Number of independent sets by rank:\n"); }
    for ($i=0; $i<= $largestindset; $i++) {
	@fvector = (@fvector,$numindsets[$i]); #adding elements to the f-vector
        if ($verbose == 1) { print ("\$i=$i. \$numindsets[$i] = $numindsets[$i]\n"); }
    }

    if ($verbose == 1) { print ("Number of dependent sets by rank:\n"); }
    for ($i=0; $i<= $largestdepset; $i++) {
        if ($verbose == 1) { print ("\$i=$i. \$numdepsets[$i] = $numdepsets[$i]\n"); }
    }

    if ($verbose == 1) { print ("Matroid rank = $largestindset\n"); }
    if ($verbose == 1) { print ("All independent sets:\n"); }
    for ($i=0; $i<= $largestindset; $i++) {
        if ($verbose == 1) { print ("    SIZE $i: "); }
        for ($j=0; $j<$numindsets[$i]; $j++) {
            if ($verbose == 1) { print ("$allindsetsnice[$i][$j], "); }
        }
        if ($verbose == 1) { print ("\n"); }
    }

    for ($i=0; $i<= $largestindset; $i++) {
        if ($verbose == 1) { print ("    SIZE $i: "); }
        for ($j=0; $j<$numindsets[$i]; $j++) {
            if ($verbose == 1) { print ("$allindsets[$i][$j], "); }
        }
        if ($verbose == 1) { print ("\n"); }
    }

    if ($verbose == 1) { print ("\n\nAll dependent sets:\n"); }
    for ($i=0; $i<= $largestdepset; $i++) {
        if ($numdepsets[$i] > 0) {
            if ($verbose == 1) { print ("    SIZE $i: "); }
        }
        for ($j=0; $j<$numdepsets[$i]; $j++) {
            if ($verbose == 1) { print ("$alldepsetsnice[$i][$j], "); }
        }
        if ($numdepsets[$i] > 0) {
            if ($verbose == 1) { print ("\n"); }

        }
    }

    for ($i=0; $i<= $largestdepset; $i++) {
        if ($numdepsets[$i] > 0) {
            if ($verbose == 1) { print ("    SIZE $i: "); }
        }
        for ($j=0; $j<$numdepsets[$i]; $j++) {
            if ($verbose == 1) { print ("$alldepsets[$i][$j], "); }
        }
        if ($numdepsets[$i] > 0) {
            if ($verbose == 1) { print ("\n"); }
        }
    }


# ********* Print h-vector *******************
    if ($verbose == 0) { print ("******************** NEW MATROID $matroidCount ********************\n");}
    if ($verbose == 0) { print ("$input\n"); }
    if ($verbose == 0) { print ("inputlength = $inputlength\n"); }

sub bc {
($n,$k)=@_;$r=1;$r*=$n/($n-$k),$n--while$n>$k;$r;
};

@hvector = (1);

for ($a = 1; $a<= $largestindset; $a++) {
    $h = 0;
    for ($i=0; $i<= $largestindset; $i++) {
	$exp = $i+$a;
	$p = $largestindset-$i;
	$l = $a-$i;
	$choose = bc($p,$l);
	$f = @fvector[$i];
	$h = $h + ((-1)**$exp) * $f *($choose);
    }
	@hvector = (@hvector,$h);
}
print("f-vector=@fvector\n");
print("h-vector=@hvector\n");


# ********* End Print h-vector ***************


# Test for simple matroid

# Initially we assume it is a simple matroid
    $isSimpleMatroid = 1;

# Run through the k elements of the matroid and if none are an independent set, then it is a loop
    for ($k=0; $k < $numelements; $k++) {
        $foundk = 0;
        for ($j=0; $j<$numindsets[1]; $j++) {
            if ($k == $allindsets[1][$j]) {
                $foundk = 1;
            }
        }
        # If $foundk == 0 then k is a loop
        if ($foundk == 0) {
            if ($verbose == 1) { print ("$k is a loop.\n"); }
            $isSimpleMatroid = 0;;
        }
    }
# Now run through all the non-loops (singleton independent sets) and see if any pair are a dependent set
    for ($i=0; $i< $numindsets[1]; $i++) {
        for ($k=$i+1; $k< $numindsets[1]; $k++) {
# For singleton independent sets $i, and $k, is that an ind set?
            $isikind = 0;
            for ($j=0; $j<$numindsets[2]; $j++) {
                if (($i == substr($allindsets[2][$j],0,1)) && ($k == substr($allindsets[2][$j],1,1))) {
                    $isikind = 1;
                }
            }
            if ($isikind == 1) {
            }

            else {
                if ($verbose == 1) { printf ("{$i,$k} are parallel.\n"); }
                $isSimpleMatroid = 0;

            }
        }
    }

    if ($isSimpleMatroid == 1) {
        if ($verbose == 1) { printf ("The matroid is simple\n"); }
    }
    else {
        if ($verbose == 1) { printf ("The matroid is non-simple\n"); }
}
# END test for simple matroid

# ********** Test for number of connected components **********

# Two elements (x,y) of the ground set are connected if
# they are contained in a circuit.

# So we need to compute the ciruits first.
# The cicuits are the minimal dependent (non-independent) sets.

# Go through all independent sets and check if they are minimal.
# It only makes sense to go through all dep sets of size 2


    $numCiruits = 0;
    for ($i=2; $i<= $largestdepset; $i++) {
        for ($j=0; $j<$numdepsets[$i]; $j++) {
            $isMinimal = 1;
            # Assume it is minimal
            for ($k=0; $k < length($alldepsets[$i][$j]); $k++) {
                $subIndSet = substrremi($alldepsets[$i][$j],$k);
                # Now check if $subIndSet is independent
                for ($l = 0; $l < $numdepsets[$i-1];$l++) {
                    if ($alldepsets[$i-1][$l] == $subIndSet) {
                        $isMinimal = 0;
                    }
                }
            }
            if ($isMinimal) {
                if ($verbose == 1) { printf("$alldepsets[$i][$j] is a circuit.\n"); }
                $allCircuits[$numCircuits] = $alldepsets[$i][$j];
                $numCircuits++;
		if ($i == 2) {
		    $a = substr($alldepsets[$i][$j], 0, 1);
		    $b = substr($alldepsets[$i][$j], 1, 1);
		}
            }
        }
    }


# Now lets go through all pairs of elements (x,y) and see if
# they are contained in some circuit.
    $isConnectedMatroid = 1;
    for ($i=0; $i<$numelements; $i++) {
        for ($j=$i+1; $j<$numelements; $j++) {
            $i_j_connected = 0;
            # if ($verbose == 1) { print ("Checking {$i,$j}\n"); }
            for ($k=0; $k<$numCircuits; $k++) {
                #if ($verbose == 1) { print ("    Looking in $allCircuits[$k]\n"); }
                if ( $allCircuits[$k] =~ /$i/) {
                    if ( $allCircuits[$k] =~ /$j/) {
                        # $allCircuits[$k] contains $i and $j
                        $i_j_connected = 1;
                    }
                }
            }
            if ($i_j_connected == 0) {
                if ($verbose == 1) { printf ("{$i,$j} not connected.\n"); }
                $isConnectedMatroid = 0;
            }
        }
    }

    if ($isConnectedMatroid == 0) {
        if ($verbose == 1) { print("The matroid is not connected.\n"); }
    }
    else {
        if ($verbose == 1) { print("The matroid is connected.\n"); }
    }

# ********** END Test for number of connected components **********


# ********** Print topcom data **********
    if ($printTopcom == 1) {
        $topcomString = "[";
        for ($j=0; $j<$numindsets[$largestindset]; $j++) {
            #print the dimension and the number of bases
            $topcomString = $topcomString . "[";
            for ($k=0; $k < $numelements; $k++) {
                if ($allindsets[$largestindset][$j] =~ /$k/) {
                    $topcomString = $topcomString . "1,";
                }
                else {
                    $topcomString = $topcomString . "0,";
                }
            }

            chop($topcomString);
            $topcomString = $topcomString . "],";
        }
        chop($topcomString);
        $topcomString = $topcomString . "];";
        if ($verbose == 1) {print ("$topcomString\n");}
        if ($writetbfiles == 1) {
            if ($writeSimpleOnly == 0 || ($writeSimpleOnly == 1 && $isSimpleMatroid == 1)) {
                if ($writeConnectedOnly == 0 || ($writeConnectedOnly == 1 && $isConnectedMatroid == 1)) {
                    open (OUTFILE, "> tb$tbcount.topcom") || die ("Couldn't open tb$tbcount.tri\n");
                    print OUTFILE ($topcomString);
                    close (OUTFILE);
                }
            }
            # Still increase tbcount so can track input to output
#            print ("$tbcount (\$isSimpleMatroid, \$isConnectedMatroid) ($isSimpleMatroid, $isConnectedMatroid)\n");
            $tbcount++;
        }
    }
# ********** END Print topcom data **********

# ********** Print latte/cdd data **********
    if ($printBases == 1) {
        $basisoutputstring = "$numelements\n$numindsets[$largestindset]\n";
        for ($j=0; $j<$numindsets[$largestindset]; $j++) {
            #print the dimension and the number of bases
            for ($k=0; $k < $numelements; $k++) {
                if ($allindsets[$largestindset][$j] =~ /$k/) {
                    $basisoutputstring = $basisoutputstring . "1 ";
                }
                else {
                    $basisoutputstring = $basisoutputstring . "0 ";
                }
            }
            chop($basisoutputstring);
            $basisoutputstring = $basisoutputstring . "\n";
        }
        print ("$basisoutputstring\n");
    }
# ********** END Print latte/cdd data **********
}


# ********************************************
# ********************************************
#          END sub royletoallinfo
# ********************************************
# ********************************************



# ********************************************
# ********************************************
#               START MAIN 
# ********************************************
# ********************************************
$verbose = 0;       # Verbose printing
$printBases = 0;    # Print latte/cdd basis information
$printTopcom = 0;   # Print topcom information
$writetbfiles = 0;  # Write out tb file for each line
$writeSimpleOnly = 0;   # Only write to tb file the simple matroids
$writeConnectedOnly = 0;    #Only write to tb file the connected matroids

$matroidCount = 0;  #Count the number of matroids processed

while ($arg = shift(@ARGV)) {
    if ($arg =~ /\-v/) {
        print ("Setting verbose.\n");
        $verbose = 1;
    }
    if ($arg =~ /\-b/) {
        $printBases = 1;
        print ("Setting basis printing.\n");
    }
    if ($arg =~ /\-t/) {
        $printTopcom = 1;
        print ("Setting topcom printing.\n");
    }
    if ($arg =~ /\-f/) {
        $writetbfiles = 1;
        $tbcount = 1;
        print ("Setting topcom file printing.\n");
    }
    if ($arg =~ /\-so/) {
        $writeSimpleOnly = 1;
        print ("Only printing simple matroids printing.\n");
    }
    if ($arg =~ /\-co/) {
        $writeConnectedOnly = 1;
        print ("Only printing connected matroids printing.\n");
    }

}
    if ($arg =~ /\-h/) {print(@hvector);
}
    if ($verbose == 1) {print ("ARGV[$i] = $arg\n");}
}

while ($myinput = <>) {
    $matroidCount++;
    royletoallinfo ($myinput);
}