#!/usr/local/bin/perl -w
# Gtest.p, from G2H.p 
# reads in a G matrix, then tries to make linearly independent sets,
# reporting failures
#
# see code/RLL
#
# Gtest.p R=2

$wmax=0;
$interval=100;
$R=5 ;
$N=7 ; $K=4 ;
$N=0;$K=0; 
$Gfile="G74.pbm" ;
$verbose = 0 ; 
$seed = 0;

eval "\$$1=\$2" while @ARGV && $ARGV[0]=~ /^(\w+)=(.*)/ && shift;

srand ( $seed ) ;

if ( $Gfile ) {
    print "# " if ( $verbose ) ;
    print STDERR "$Gfile " ;
    print "\n" if ( $verbose ) ;
    open ( H , "< $Gfile" ) ; 
    $i = 0 ; $k = 0 ; 
    
    while ( <H> ) {
	if ( /P/ ) {	# read in the header
	    if ($k) { print STDERR "unexpected character 'P' in $_ \n" ; }
	    if($verbose) { print STDERR "reading header\n" ; } # 
	    if($verbose) { print STDERR ;}
	    $_ = <H> ; 
	    if($verbose) { print STDERR ;}
	    $_ = <H> ;
	}
	@r = split ;
	print STDOUT "# "  if ( $verbose > 2 ) ; 
	$Nr = $#r + 1 ; 
	for ( $n = 1 ; $n <= $Nr ; $n ++ ) {
	    $Gi[$i] = $r[$n-1] ;
	    $Gk{$k+1,$n} =  $r[$n-1] ;
#	    $Gn{$n,$k+1} =  $r[$n-1] ;
	    print STDOUT "$Gi[$i] "  if ( $verbose > 2 ) ;
	    $i ++ ;
	}
	print STDOUT "\n"  if ( $verbose > 2 ) ; $k ++ ;
    }
    print STDERR "# read in $k rows \n" ; 
    close ( H ) ;
}
# Gi[0..N*K-1] vector contains the whole matrix
if ( ($K && !($k == $K)) || ($N && !($Nr == $N) ) ) {
    print "WARNING: K N mismatch $k $Nr\n" ; 
    exit(0) ; 
}
$K = $k ; 
$N = $Nr ; 
if ( $verbose ) {
    print "K = $K, N = $N\n" ; 
}

# make an urn full of n's
for ( $n = 1 ; $n <= $N ; $n ++ ) {
    $mastern[$n-1] = $n ;
}
if ( $wmax<=0) {
    $wmax = $K+1  ; # just +1 is enough...
}
for ($w=1;$w<=$wmax;$w++) {
    $tried[$w] = 0 ; 
    $failed[$w] = 0 ; 
    $succeeded[$w] = 0 ;
}
for ($r=1;$r<=$R;$r++) { # make random vectors
    print  "\n $r :::: " if ($verbose) ; 
    print STDERR  "$r." if (!($r%$interval)) ; 
    @nlist = @mastern ;

WLOOP:    for ($w = 1 ; $w <= $wmax ; $w ++ ) {
	# select the next entry
        $q = int ( rand($#nlist+1) ) ;
	$n = splice ( @nlist , $q , 1  ) ;
        $t[$w] = $n ; # record the transmitted bits we are concerned about
	print  "$n " if ($verbose) ; 
	print  ": " if ($verbose>1) ; 
	# copy across n's bits of G
	for ($k=1;$k<=$K;$k++) {
	    $thisG{$k,$w} = $Gk{$k,$n} ;
	    print "copied thisG $k,$w = $thisG{$k,$w}\n" if($verbose>10) ; 
	}

	if($verbose>1) {
	    print  "pre-elimination G \n" ;
	    for($wo=1;$wo<=$w;$wo++) {
		print "$t[$wo] " ;
	    }
		print "\n" ; 
	    for ($k=1;$k<=$K;$k++) {
		for($wo=1;$wo<=$w;$wo++) {
		    print "$thisG{$k,$wo} " ; 
		}
		print "\n" ; 
	    }
	}
	# run through the leaders of the already-chosen (if any)
	for ($wo = 1 ; $wo < $w ; $wo ++ ) {
	    $kl = $leader[$wo] ;
	    if ( $thisG{$kl,$w} ) {
		for ( $k=1;$k<=$K;$k++) {
		    print  "(k=$k,w=$w,$wo):$thisG{$k,$w}:$thisG{$k,$wo}\n"  if($verbose>3);
		    if ( $thisG{$k,$wo}) { $thisG{$k,$w} = 1 -  $thisG{$k,$w}  ; }
		}
	    } else {
		print  "(no clash,w=$w,$wo)\n" if($verbose > 3) ;
	    }
	}
	if($verbose>0) {
	    print  "post-elimination G \n" ;
	    for($wo=1;$wo<=$w;$wo++) {
		print "$t[$wo] " ;
	    }
		print "\n" ; 
	    for ($k=1;$k<=$K;$k++) {
		for($wo=1;$wo<=$w;$wo++) {
		    print "$thisG{$k,$wo} " ; 
		}
		print "\n" ; 
	    }
	}
	# check to see if this n now has any bits left
	$sum=0;
	for ($k=1;(!$sum) && $k<=$K;$k++) {
	    $sum += $thisG{$k,$w} ;
	}
	$tried[$w] ++;
	if ($sum) {
	    $k-- ; # this tells us the new leader
	    $leader[$w] = $k ;
	    print  "leader $w is $k\n" if($verbose>1) ;
	    $succeeded[$w] ++;
	}	else {
	    $failed[$w] ++;
	    if($verbose) {
		print  " - the $w th bit, $n, is dependent \n" ; 
	    }			# 
	    if($verbose>0) {
		print "-------------------\n" ; 
	    }
	    last WLOOP ;	# 
	}
	
    }

}

print STDERR "\n" ; 
print "#w\ttried\tfailed\tcfailed\tsucceeded\tfailure rate\toverall failure fraction\toverall success fraction\n" ;
$cfail = 0 ; 
for ($w=1;$w<=$wmax;$w++) {
    print " $w\t" ; 
    print "$tried[$w]\t"  ; 
    print "$failed[$w]\t" ;
    $cfail +=  $failed[$w] ;
    $cfailed[$w] = $cfail ; 
    print "$cfailed[$w]\t" ; 
    print "$succeeded[$w]" ;
    if ( $tried[$w] ) {
	printf "\t%8.4g" , ( $failed[$w]/$tried[$w]) ; 
	printf "\t%8.4g" , ( $cfailed[$w]/$tried[1]) ; 
	printf "\t%8.4g" , ( $succeeded[$w]/$tried[1]) ; 
    }
    print "\n" ; 
}