Matlab/Octave code

Cameron Freer, Travis Hee Wai, Bjørn Kjos-Hanssen

	function fi=OL(A)
	[k n]=size(A);
		for a=1:k
			b=abs(A)*abs(A(a,:))';
			b=b>0;
			fi(a)=sum(b)-1;
		end
	end
	function x=randint(a,b,c)
		if nargin==0;
			x=round(rand());		%	A random bit, since rand() returns a number in the interval (0,1).
		elseif nargin==2;
			x=floor(2*rand(a,b));	%	An axb matrix of random bits
		else
			x=floor(c*rand(a,b));	%	An axb matrix of random integers in the interval {0,...,c-1}
		end
	end

	function k=fourdp(A)
		fi=OL(A);
		d=max(fi);
		p=2^(-min(sum(abs(A),2)));
		k=min(4*d*p,exp(1)*(d+1)*p);
	end

	function A=rand_eq (m,n,k)

		maxtrials=1000;
		trials=0;
		fourdp_check=2;

		while fourdp_check > 1

			A=zeros (m,n);    		%	mxn matrix of zeros

		%	Generate a random matrix A

			for a=1:m
				while sum (abs (A (a,:))) < k		% no more than k literals per clause, i.e., per row
					index = randint (1,1,n)+1;	% index is a random integer between 1 and n
					A (a,index)=2*randint ()-1;	% change column "index" in row "a" to a random value among {1,-1},i.e., a literal or its negation
				end
			end

			%   Check for duplicate clauses
			for a=1:m
				% Check current clause against all other clauses
				check=ones(m,1)*A(a,:);			% mx1 matrix times 1xn matrix = mxn matrix, all of whose rows equal A(a,:) ?
				check_eq=check.*A;				% multiply entrywise, not matrix multiplication. Note ab=1 iff a=b when a,b in {+-1}
				repeats=sum(check_eq,2);			% 2 means sum the rows and put the sums into a column vector
				indexes=find(repeats==n);		% find the set of rows that sum to n, the no. of columns
				if isempty(indexes)==0			% i.e., if row a is repeated at least once(!)
					k=length(indexes);
					for b=1:k
						while sum (abs (A (b,:))) < k
							index = randint (1,1,n)+1;
							A (b,index)=2*randint ()-1;	% change each of the first k many rows 
						end							% (until they have enough literals, which they already should have
					end								% so should probably reset A(b,:) to be 0 first. )
				end									% Also, should probably change A(indexes(b),:) not A(b,:).
			end

			fourdp_check=fourdp (A);
			trials=trials+1;

			% decrease m by 3/4 if too many trials have passed
			if trials==maxtrials
				m=floor(m*3/4);
				trials=0;
			end
		end
	end



	% %  Function mtdata.m
	% %  Given: m - number of equations
	% %         n - number of variables
	% %         k - number of variables in any one equation
	% %         max - the number of iterations
	% %  Returns: sol_dist - distribution of solutions
	% %           trial_dist - average trials for each unique solution
	% %           clauses_sampled - number of times each clause was resampled

	function [sol_dist,trial_dist,clauses_sampled A]=mtdata(m,n,k,max,graph)

		% We require m,n,k to be supplied as inputs, but if max and graph are not provided we set max=10000 and graph=0:
		if nargin==3
			max=10000;
			graph=0;
		elseif nargin==4
			graph=0;
		end

		A=rand_eq(m,n,k);
		[m n]=size(A);
		solutions=zeros(max,n);
		trial_stats=zeros(max,1);
		clauses_sampled=zeros(1,m);
		for i=1:max    
			% Get current solution and number of trials
			[y trials lines_updated]=mt(A);
			solutions(i,:)=y;
			trial_stats(i)=trials;
			clauses_sampled=clauses_sampled+lines_updated;
		end

		b=0;

		for a=1:max
			% Check current solution against all
			check=ones(max,1)*solutions(a,:);
			% Check frequency of solution
			check_sol=check.*solutions;
			repeats=floor(sum(check_sol')/n);
			repeats=(repeats==1);
			indexes=find(repeats==1);
			if a<=min(indexes)
				b=b+1;
				sol_dist(b)=sum(repeats);
				trial_dist(b)=sum(trial_stats.*repeats');
			end
		end

		trial_dist=trial_dist./sol_dist;
		if graph~=1
			figure
			subplot(3,1,1)
			bar((sol_dist))
			title('Distribution of Solutions')
			subplot(3,1,2)
			bar((trial_dist),'r')
			title('Average Trials for each Solution')
			subplot(3,1,3)
			bar(clauses_sampled,'g')
			title('Distribution of Resampled Clauses')
			if graph~=2
				figure
				dplot(A)
			end
		end
	end

	% % Function mt(A) 
	% % Given: A matrix where a(i,j) represents P(j) for equation i.  
	% %        a(i,j)={1, P(j) in equation i
	% %               {0, P(j) not in equaiton i
	% %               {-1,-P(j) in equation i
	% %       trials - number of iterations to try before quitting
	% %
	% % Returns: x - The solution it found.  If there is no solution, it returns the
	% %          empty set
	% %          num_trials - the number of guesses it took to find a right solution
	% %          line_calcs - the number of equations checked before a right
	% %          solution was found

	function [x num_trials line_check]=mt(A,max_trials)

		if nargin==1
		% max number of iterations to try
		max_trials=1000000;
		end

		[k n]=size(A); % k the number of equations, n the number of variables
		x=[];

		% Sort A by dependence
		[d,ix]=sort(OL(A));
		for b=k:-1:1
			newA(b,:)=A(ix(b),:);
		end
		A=newA;

		% Generate a random solution
		p=(randint(1,n)-1/2)*2;

		is_solution=0;
		trials=0;
		line_check=zeros(1,k);
		while is_solution~=1

			% Check if p is a solution to the set of equations given in A
			y=solve(p,A);

			% Break the loop if p is a solution 
			is_solution=min(y);

			if is_solution==0
				% Search for equations not met
				not_satisfied=find(y==0);  
				l=not_satisfied(1);   % pick first unsatisfied equation 
				line_check(l)=line_check(l)+1;
				% Replace variables used in eq A(b)
				for c=1:n
					if A(l,c)~=0
						p(c)=2*randint()-1;
					end
				end
			end

			trials=trials+1;
			%Terminate if max trials exceeded  
			if trials==max_trials+1
				is_solution=2;
			end
		end

		num_trials=trials;

		if is_solution==1
			x=p;
		end
	end

	% % function D=dependency(A)
	% %     Calculates the dependency matrix D, from a given equation matrix A

	function D=dependency(A)
		[k n]=size(A);
		D=zeros(k,k);
		for a=1:k
			R=abs(A)*abs(A(a,:))';
			C=R>0;
			D(:,a)=C;
		end
		D=D-eye(k);
	end

	% Function solve(p,A)
	% Given: A matrix where A(i,j) represents P(j) for equation i.  
	%        A(i,j)={1, P(j) in equation i
	%               {0, P(j) not in equaiton i
	%               {-1,-P(j) in equation i
	%        p is a 
	% Returns: y, a vector containing the solution 

	function y=solve(P,A);
		[k n]=size(A);
		P1=ones(k,1)*P;
		% Calculate whether all equations are satisfied
		B=A.*P1;
		C=(B==1);
		y=max(C');        
	end

	% % function D=dplot(A)
	% %         Plots the dependency graph between clauses for equations given
	% %         in matrix A
	function dplot(A)
		D=dependency(A);
		[k n]=size(D);
		load V
		gplot(D,V(1:k,:),'b-');
		for j = 1:k,
			text(V(j,1),V(j,2),int2str(j),'FontSize',10);
		end
		axis off equal
	end