#!/usr/bin/perl
#!/usr/local/bin/perl -w
#
# create random SAT problem and enumerate all solutions
#
# Use written datafiles and gnuplot to render pretty animations
#     of the solution set.
#
# usage
# solns.p M=10


$M=10;
$T=200;
$dT=1;
$gfile="g/g" ;
$vfile="o/v" ;
$efile="o/e" ;
$verbose=0;
eval "\$$1=\$2" while @ARGV && $ARGV[0]=~ /^(\w+)=(.*)/ && shift;
sub bynumber { $a <=> $b; }

sub bytot { $tot[$a] <=> $tot[$b] ; }

$N = 1 << $M;  # N=2^M
print $N,"\n" ;
$gfile.="$M" ;
open ( G , ">$gfile" ) ;
print G "# This gnuplot file written by solns.p (c) David MacKay \n" ;
print G "setmeto1 = 1 \n" ;
print G "set nokey\n set autos xy \n" ;
print G "set xrange [] writeback\n" ;
print G "set yrange [] writeback\n" ;
for ( $n = 0 ; $n < $N ; $n ++ ) {
    $alln[$n] = $n ;
    $nn = $n ; $s="" ;
    for ( $m = 0 ; $m < $M ; $m ++ ) {
	$r = ($nn%2) ; $nn = int($nn/2) ; 
	$s = ( $r ? "1" : "0" ) .$s  ; 
    }
    $binary[$n] = $s ;
    print "$n = $s \n" if ($verbose) ; 
}
@b = sort bynumber ( @alln ) ; 
print join ( "." , @b ) if($verbose>3) ; 

# make vectors for projecting from M dim into 2 or 4 (use 4 to enable animations.
for ( $m = 0 ; $m < $M ; $m ++ ) {
    for ($d=1;$d<=4;$d++) {
	$proj{$d,$m} = rand()-0.5 ;
    }
}

# compute all projected locations
for ( $n = 0 ; $n < $N ; $n ++ ) {
    $locn{$n} = "" ;
    for ($d=1;$d<=4;$d++) {
	$ld = 0.0 ;
	$nn = $n ; 
	for ( $m = 0 ; $m < $M ; $m ++ ) {
	    $r = ($nn%2) ; $nn = int($nn/2) ;
	    if($r) { $ld += $proj{$d,$m} ; }
	}
	$locn{$n} .= sprintf("%5.2f ", $ld ) ; 
    }
    $locn{$n} .= "\n" ;
}

print "\nStarting\n" ; 
	$firsttime = 1 ;
TLOOP: for ($t=0;$t<$T; ) {
    # make dT more constraints and advance t to keep count of total
    for ( $c = 1 ; $c <= $dT ; $c ++ , $t ++ ) {
# make random constraints
	$u = int( rand() * $N ) ;
	print $u . "\n" if ($verbose) ;	# this is a random binary vector
	$i1 = int( rand() * $M ) ;
	do {	    $i2 = int( rand() * $M ) ;
	}	    while ($i1 == $i2 ) ;
	do {	    $i3 = int( rand() * $M ) ;
	}	    while (($i1 == $i3 )||($i1 == $i2 )) ;
	$v = (1 << $i1 ) +  (1 << $i2 ) +  (1 << $i3 ) ;
	$u = $v & $u ; # throw away all elements of u not aligned with v
	# run through all hypotheses
	for ( $n = 0 ; $n < $N ; $n ++ ) {
	    $f = $u^($n&$v) ; # this carries out the flips of bits
	    $ans{$n,$t} = ($f>0) ? 1 : 0 ; # ans is whether true
	    $tot[$n] += ($f>0) ? 0 : 1 ; # tot is the number of violations
	    ## would like tot to be zero
	    if($verbose>2) {
		print "for state n=$n, v=$v ($i1,$i2,$i3), u=$u: ans=$f\n" ;
	    }
	}
    }
    if($verbose>2) {
	for ( $n = 0 ; $n < $N ; $n ++ ) {
	    printf "%4d: " ,$n ;
	    for ($tt=1;$tt<=$t;$tt++ ) {
		printf "%d" , $ans{$n,$tt} ;
	    }
	    print " $tot[$n]\n" ;
	}
    }
    @b = sort bytot ( @alln ) ;
    if($verbose) {
	print join ( "\n" , @b ) if($verbose>3) ;
	print "\nSorted states by violations\n" ; 
	for ( $n = 0 ; $n < $N ; $n ++ ) {
	    printf "%4d %s %3d  " ,  $b[$n] , $binary[$b[$n]], $tot[$b[$n]] ;
	    for ($tt=0;$tt<$t;$tt++ ) {
		printf "%d" , $ans{$b[$n],$tt} ;
	    }
	    print "\n" ; 
	}
    }
    if ($tot[$b[0]] > 0 ) {
	print "time $t, no solutions left\n" ;
	last TLOOP ;
    }
    ## write file for gnuplot
    $vlistout[0] = "" ;
    $vlistout[1] = "" ;
    $vlistout[2] = "" ;
    $vlistout[3] = "" ; $vtotmax = 3 ;
    $elistout[0] = "" ;
    $elistout[1] = "" ; $etotmax = 1 ;
NNLOOP:    for ( $nn = 0 ; $nn < $N ; $nn ++ ) {
	$n = $b[$nn] ;
	$to = $tot[$n] ;
	print "n=$n has tot=$to\n" if($verbose>2);
	if ($to > $vtotmax ) { last NNLOOP ; } #
	else  {
	    $vlistout[$to] .= $locn{$n} ; ## plot this with bright dots
	    ##   or ($tot[$n] == 1 )
	    if ($to == 0 ) {
	    ## follow all edges from this dude
		for ( $m = 0 ; $m < $M ; $m ++ ) {
		    $nnn = $n ^ ( 1 << $m ) ;
		    $neitot = $tot[$nnn] ;
		    print "n=$n has tot=$to and neighbour($m) $nnn who has $neitot\n" if($verbose>2);
		    if ( ($nnn > $n) && ( $neitot <= $etotmax ) ) {
			$elistout[$neitot] .=  $locn{$n}.$locn{$nnn}."\n" ;
		    }
		}
	    }
	}
    }
    open ( V , ">$vfile.$t" ) ;
    open ( E , ">$efile.$t" ) ;
    print V "# index 0\n".$vlistout[0]."\n"."\n"."# index 1\n".$vlistout[1]."\n"."\n"."# index 1\n".$vlistout[2]."\n"."\n"."# index 1\n".$vlistout[3] ;
    close(V);
    print E "# index 0\n".$elistout[0]."\n"."\n"."# index 1\n".$elistout[1] ;
    close(E);
    print G "set title '3-SAT: $M variables, $t clauses' \n" ;
    print G "plot \\\n" ;
    print G "'$vfile.$t' index 3 u 1:2 w d 2, \\\n" ; 
    print G "'$vfile.$t' index 2 u 1:2 w d 3, \\\n" ; 
    print G "'$efile.$t' index 0 u 1:2 w l ls 7, \\\n" ; 
    print G "'$efile.$t' index 1 u 1:2 w l ls 1, \\\n" ; 
    print G "'$vfile.$t' index 1 u 1:2 w p 4, \\\n" ; 
    print G "'$vfile.$t' index 0 u 1:2 w p 5 \n" ; 
    print G "pause -1 'ready'\n" ;
    if($firsttime) {
	print G "set noautos xy \n" ;
	$firsttime = 0 ;
    }

    if($t>1) {
	open( H , ">>$oldmasterfile" ) ;
	print H "load '$nextmasterfile' \n" ;
    }
    $tt = $t+1 ;
    $masterfile = "$gfile.$t" ; 
    $oldmasterfile =     $masterfile  ; 
    $nextmasterfile = "$gfile.$tt" ; 
    $plotfile = "$gfile.$t.p" ; 

    open( H , ">$masterfile" ) ;
    if ($t==1) {
	print H "setmeto1 = 1 \n" ;
    }
    print H "if ( setmeto1 ) set nokey ; set autos xy \n" ;
    print H "if ( setmeto1 ) set xrange [] writeback ; set yrange [] writeback\n" ;
    print H "set title '3-SAT: $M variables, $t clauses' \n" ;
    print H " MAXTIME=100; amp1 = 0.051; amp2=0.051 ; omega1 = 2*2*3.14159 / MAXTIME  ; omega2 = 2*2*3.14159 / MAXTIME ; \n" ;
    print H "if ( setmeto1 )  time1 = MAXTIME ; load '$plotfile' \n" ;
    print H "if ( setmeto1 )  set noautos xy; setmeto1=0  \n" ;
    print H " time1 = 0 ; load '$plotfile' \n" ;
    close(H) ;
    open( H , ">$plotfile" ) ;
    print H " th1 = amp1 * sin(time1*omega1) ;th2= amp2 * sin(time1*omega2) ;  s1 = sin(th1) ; c1 = cos(th1) ;  s2 = sin(th2) ; c2 = cos(th2) ;  \n" ;
    print H "plot \\\n" ;
    print H "'$vfile.$t' index 3 u (s1*\$3+c1*\$1):(s2*\$4+c2*\$2) w d 2, \\\n" ; 
    print H "'$vfile.$t' index 2 u (s1*\$3+c1*\$1):(s2*\$4+c2*\$2) w d 3, \\\n" ; 
    print H "'$efile.$t' index 0 u (s1*\$3+c1*\$1):(s2*\$4+c2*\$2) w l ls 7, \\\n" ; 
    print H "'$efile.$t' index 1 u (s1*\$3+c1*\$1):(s2*\$4+c2*\$2) w l ls 1, \\\n" ; 
    print H "'$vfile.$t' index 1 u (s1*\$3+c1*\$1):(s2*\$4+c2*\$2) w p 4, \\\n" ; 
    print H "'$vfile.$t' index 0 u (s1*\$3+c1*\$1):(s2*\$4+c2*\$2) w p 5 \n" ; 
    print H "time1 = time1 + 1 ; if(time1<MAXTIME) reread  \n" ;
    close(H) ;
    
    print  "load '$masterfile' \n" ;
    print  "Done step $t\n" if($verbose) ;
}
if($t>1) {
    open( H , ">>$oldmasterfile" ) ;
    print H "setmeto1 = 1 \n" ;
}

close(G) ;
print "# gnuplot   \n load '$gfile' \n" ;