# 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);