#!/usr/local/bin/perl -w
# GH.p
#
# makes GHmaker.p which makes alist matrices based on permutations.
# I think this is superceded by GHG.p
#
# bugs: it should check to see if you are asking for out-of-range
#       locations
#
# see ~mackay/pub/gallager/GH.html or http://www.inference.phy.cam.ac.uk/mackay/gallager/GH.html 
#
# write them to code/GHC, please
#
# see also GHdraw.p
#
$N=10000 ;
$file="GHmaker.p";
eval "\$$1=\$2" while @ARGV && $ARGV[0]=~ /^(\w+)=(.*)/ && shift;

print STDERR "===== GH.p ===== (c) DJCM april 1998 =====\n" ;
#
# read the spec file
#
$status = 0 ; $n = 0 ; $m = 0 ; $lambda = 0 ;
$c = 0 ; # command number
$dlist = ":" ; $dprod = 1 ; 
while(<>) {
    s/^\s*// ; 
    if (/^#/) {
#	print STDERR ;
	next ;
    }
    s/\#.*// ;		#
    if (/lambda/) {
	print STDERR " lambda = " ;
	($n,$m) = /lambda\s+(.*)\/(.*)/ ;
	print STDERR "$n / $m \n" ;
	if ( $m <= 0 ) { $status -- ; }
	else {
	    $lambda = ( 1.0 * $n ) / ( 1.0 * $m ) ; 
	}
	&denoms_add ( $m ) ; 
    } elsif  (/refresh/) {
	$c ++ ;
	$type[$c] = "refresh" ; 
	print STDERR "refresh\n" ;
    } elsif  (/forbidI/) {
	$c ++ ;
	$type[$c] = "forbidI" ;
	($sizen[$c],$sized[$c])  =  /forbidI\s+(.*)\/(.*)/ ;
	print STDERR "forbid permutations like I $sizen[$c]/$sized[$c]\n" ;
	&denoms_add ( $sized[$c] ) ; 
    } elsif  ( /[IP].*at/ ) {
	$c ++ ;
	($cols[$c],$rows[$c],$type[$c],$sizen[$c],$sized[$c],
	 $xn[$c],$xd[$c],$yn[$c],$yd[$c]) =
	     /(.*)x(.*)\s+([IP])\s+(.*)\/(.*)\s+at\s+(.*)\/(.*),(.*)\/(.*)/ ;
	print STDERR " $cols[$c] x $rows[$c] $type[$c] size $sizen[$c]/$sized[$c] at $xn[$c]/$xd[$c],$yn[$c]/$yd[$c]\n" ;
	&denoms_add ( $sized[$c] ) ; 
	&denoms_add ( $xd[$c] ) ; 
	&denoms_add ( $yd[$c] ) ;
    }
}
sub denoms_add { # should really check for divisibility here!
    local ($d) = @_ ;
    if (!($dlist =~ /:$d:/) ) {
	$dlist = $dlist."$d:" ;
	$dprod = $dprod * $d ;
    }
}
$C = $c ;
print STDERR "# Read $C commands\n" ;
print STDERR "# Make M divisible by $dprod\n" ;
	
# suggested value for $M
$M = $N * $m / $n ;
$M = ( int ( $M / $dprod ) ) * $dprod ;
$N = $M * $n / $m ;
print STDERR "# eg M = $M, N = $N\n" ;
print STDERR "# $file seed=123 M=$M > a\n" ;
# exit(0) ; 

open (FI, ">$file") ;
#
# write the first few lines
#
print FI '#!/usr/local/bin/perl -w
# GHmaker.p, written by GH.p
#
# makes alist matrices based on permutations.
# see ~mackay/pub/gallager/GH.html or http://www.inference.phy.cam.ac.uk/mackay/gallager/GH.html 
$seed="123";
$verbose = 0 ; 
$M=' ;
print FI "$M" ;
print FI ';
eval "\$$1=\$2" while @ARGV && $ARGV[0]=~ /^(\w+)=(.*)/ && shift;
print STDERR "===== GHmaker.p ===== (c) DJCM april 1998 =====\n" ;
srand ( $seed ) ;
$m=' ;
print FI "$m" ;
print FI ';
$n=' ;
print FI "$n" ;
print FI ';
$N = $M * $n / $m ; 
$head = $N." ".$M."\n" ;

# make null lists 
for ( $n = 1 ; $n <= $N ; $n ++ ) {
    $nlist[$n] = "" ; 
    $numnlist[$n] = 0 ;
}
for ( $m = 1 ; $m <= $M ; $m ++ ) {
    $mlist[$m] = "" ; 
    $nummlist[$m] = 0 ;
}

$failed=0;
&refresh_perm() ;
# define some useful subroutines
sub refresh_perm {
    $nrandperms = 0 ;
    $nnn = 1 ;
    # permutations are recorded in
    # an array $forbid[$nnn] 
}
sub nullperm {
    local ($Nl) = @_ ;
    for ( $nl = 1 ; $nl <= $Nl ; $nl ++ ) {
	$fto[$nl] = 0 ;
    }
}
sub checkperm {
    local ($Nl) = @_ ;
    for ( $nl = 1 ; $nl <= $Nl ; $nl ++ ) {
	if ( $fto[$nl] != 1 ) {
	    print STDERR  "perm $q not valid $nl\n" ; 
	    for ( $n = 1 ; $n <= $Nl ; $n ++ ) { 
		print STDERR "$n to $p[$n]\n" ;
	    }
	    exit(0) ; 
	}
    }
}

# this is for checking of permutation clash errors
$junking = 0 ; 
if ($junking) {
    open ( JUNK , "> junk" ) ;
}
' ;
#
# end print FI
#
for ( $c = 1 ; $c <= $C ; $c ++ ) {
    if ( $type[$c] eq "refresh" ) {
	print FI "&refresh_perm() ; \n" ;
    } elsif ( $type[$c] eq "I" ) {
	print FI "&do_I($cols[$c],$rows[$c],$sizen[$c],$sized[$c],
	      $xn[$c],$xd[$c],$yn[$c],$yd[$c]) ;  \n" ;
    } elsif ( $type[$c] eq "forbidI" ) {
	print FI "&registerI($M*$sizen[$c]/$sized[$c]) ;  \n" ;
    } elsif ( $type[$c] eq "P" ) {
	print FI "&do_P($cols[$c],$rows[$c],$sizen[$c],$sized[$c],
	      $xn[$c],$xd[$c],$yn[$c],$yd[$c]) ; \n" ; 
    } else {
	print STDERR "oops, unknown type $type[$c]\n"  ;
    }
}

print FI '#
sub do_I {
    local($Co,$Ro,$sn,$sd,$xn,$xd,$yn,$yd) = @_ ;
    local($n0) = 1 + $M * $xn / $xd ; 
    local($m0) = 1 + $M * $yn / $yd ;
    local($w) = $M * $sn / $sd ; 
    for ( $co = 1 ; $co <= $Co ; $co ++ ) {
	local($m) = $m0 ; 
	for ( $ro = 1 ; $ro <= $Ro ; $ro ++ ) {
	    local($n) = $n0 ; 

	    for ( $ww = 1 ; $ww <= $w ; $ww ++ , $n ++ , $m ++ ) {
	    
		# put m and n on each others lists 
		$nlist[$n] .= $m."\t" ;
		$mlist[$m] .= $n."\t" ;
		$numnlist[$n] ++ ;
		$nummlist[$m] ++ ;
	    }

	}
	$n0 += $w;
    }				# 
    # each time I is called,
    # register the identity as a randperm
    &registerI($w) ;
}
sub registerI {
    local ($w) = @_ ;
    for ( $nl = 1 ; $nl <= $w ; $nl ++ ) {
	$forbid[$nnn] = $nl ;
	$nnn ++ ;
    }
    $nrandperms ++ ;
}

sub do_P {
    local($Co,$Ro,$sn,$sd,$xn,$xd,$yn,$yd) = @_ ;
    local($n0) =  $M * $xn / $xd ; 
    local($m0) =  $M * $yn / $yd ;
    local($w) = $M * $sn / $sd ; 
    local($n) = $n0 ; 
    for ( $co = 1 ; $co <= $Co ; $co ++ , $n += $w ) {
	local($m) = $m0 ; 
	for ( $ro = 1 ; $ro <= $Ro ; $ro ++ , $m += $w ) {

# invent a perm
	    printf STDERR "Starting perm %d\n" , $nrandperms+1 ; 
	    &nullperm($w);
	    # create a list of integers
	    for ( $q = 1 ; $q <= $w ; $q ++ ) {	$qs[$q-1] = $q ; }
	    for ( $nl = 1 ; $nl <= $w ; $nl ++ ) {
		$from = $nl ;
# select at random from list qs, checking for no clash with list
	      ALOOP: for ( $a  = 1 ; $a <= 100 ; $a ++ )  { 
		  $valid = 1 ; 
		  $q = int ( rand($#qs+1) ) ;
		  $to = $qs[$q] ; 
		  if ( $verbose >= 2 ) {
		      print STDOUT  "$nl($#qs): new sample is $q , $to \n" ; 
		  }
# check if this clashes with an existing permutation.
		BLOOP: for ( $b = 1 , $testn = $nnn - $w ;
			    $b <= $nrandperms ;
			    $b++ , $testn -= $w  )  {
		    if ( $forbid[$testn] == $to ) {	# 
			print STDERR "C" ;
			if ($verbose>=2){
			    print STDERR "clash $nl: [$q] $to X\n" ;
			    print STDOUT "sampling from $#qs+1 : " , join ( " " ,@qs )  , "\n"  ; 
			}
			$valid = 0; 
			last BLOOP ; 
		    } else {
			if ($junking) {
			    print JUNK "$forbid[$testn] $to\n" ;
			}
		    }
		}		# 
		  if ( $valid || ( $nl == $w ) ) {
		      last ALOOP ;
		  }
	      }			# 
		$too = splice ( @qs , $q , 1  ) ;
		if ( $to != $too ) { print STDERR "wot? $to $too\n" ; }
		# delete selection from the list.
		if ( !$valid ) { print STDERR  "accepted invalid $nl: $to\n" ; $failed ++ ; }
		$forbid[$nnn] = $to ; $nnn ++ ;
		$p[$from] = $to ; 
		$fto[$to] ++ ;	# 
#	print STDERR $nl , "->" , $to , " " ; 
		if($verbose>=1){print STDERR "." ; }
	    }		# 
	    print STDERR  "\n"  ; 
  
	    &checkperm($w);
	    &permadd($w,$n,$m);
	    $nrandperms ++ ;
	}
    }

}

print STDERR "failures: $failed. Writing code.\n" ;

# find $t (max) and $tr (max)

$t = 0 ; $tr = 0 ; 
for ( $n = 1 ; $n <= $N ; $n ++ ) {
    if ( $numnlist[$n] > $t ) {
	$t = $numnlist[$n] ;
    }
}
for ( $m = 1 ; $m <= $M ; $m ++ ) {
    if ( $nummlist[$m] > $tr ) {
	$tr = $nummlist[$m] ;
    }
}
$head .= "$t $tr\n" ;
# append zeros to short lines
for ( $n = 1 ; $n <= $N ; $n ++ ) {
    $thisn = $numnlist[$n] ;
    $head .= "$thisn " ; 
    while ( $thisn < $t ) {
	$nlist[$n] .= "0\t" ;
	$thisn ++ ;
    }
}
$head .= "\n" ;

for ( $m = 1 ; $m <= $M ; $m ++ ) {
    $thism = $nummlist[$m] ; 
    $head .= "$thism " ; 
    while ( $thism < $tr ) {
	$mlist[$m] .= "0\t" ;
	$thism ++ ;
    }
}
$head .= "\n" ;


sub permadd {
    local ( $Nl , $noff , $moff ) = @_ ; 
    for ( $nl = 1 ; $nl <= $Nl ; $nl ++ ) {
	local($n) = $noff + $nl ;
	local($m) = $p[$nl] + $moff ; 
	if ( ( $m > $M ) ||  ( $n > $N ) ) {
	    print STDERR "invalid m,n = $m,$n ($M,$N)\n" ;
	} else {
	    # put m and n on each others lists 
	    $nlist[$n] .= $m."\t" ;
	    $mlist[$m] .= $n."\t" ;
	    $numnlist[$n] ++ ;
	    $nummlist[$m] ++ ;
	}
    }
}

print $head ; 
# finish lists 
for ( $n = 1 ; $n <= $N ; $n ++ ) {
    $nlist[$n] =~ s/\t$/\n/ ; 
    print $nlist[$n] ; 
}
for ( $m = 1 ; $m <= $M ; $m ++ ) {
    $mlist[$m] =~ s/\t$/\n/ ; 
    print $mlist[$m] ; 
}

' ;
# end print FI

system ( "chmod 755 $file" ) ; 
    exit(0);