These MATLAB 7.12.0 functions are all the functions required to run Chubanov's LP Algorithm, his Algorithm 5.1 in his 2011 paper "A Strongly Polynomial Algorithm for Linear Systems Having a Binary Solution".  After unzipping the files, if you put them all in the same directory, then you can run Chubanov.m with the proper input (explained below) and determine if the system has a feasible or determine it has no integer solutions.  The rest of this README explains each function in detail.

Before proceeding, many of these explanations are written under the assumption you have at least read our paper "On Chubanov's Method for Linear Programming", for which these functions were created.  This paper is also available at www.math.ucdavis.edu/~mjunod.

******************************************

Chubanov.m

[sol, no_int, time_fail] = Chubanov(n,m,l,A,b,C,d,r)

Although this was originally designed to be run as a stand-alone function, it has since been integrated into the function Run_Experiment.m (see its README for more info).  As such, it is easier to run specific examples from Run_Experiment.m.

Given an n-dimensional system Ax = b, Cx <= d, where A has m rows and C has l rows, r must be a radius such that the feasible region of the system is contained in B(0,r), the ball of radius r centered at the origin.  Given this input, the function either returns a feasible point sol = x with no_int = time_fail = 0, or sets no_int = 1 to indicate there are no integer solutions, or sets time_fail = 1 if the function runs longer than the time limit set in an experiment script that calls Chubanov.m.

INPUT:		n = dimension of the linear feasibility problem
		m = number of equality constraints
		l = number of inequality constraints
		A = matrix of equality constraints
		b = right-hand side of equality constraints
		C = matrix of inequality (<=) constraints
		d = right-hand side of inequality (<=) constraints
		r = positive number such that the feasible region is contained in the ball B(0,r)

OUTPUT:		[sol, 0, 0] = in this case a feasible point sol has been found
		[[], 1, 0] = in this case the function determined no integer solutions exist
		[[], 0, 1] = in this case the function timed out, i.e. took longer than the preset time limit

******************************************

DandC.m

[x,h,delta,alpha,flag,fail,fail_alpha,time_fail] = DandC(n,m,l,A,b,C,d,z,r,error,c_max,alpha,fail_alpha)

Given an n-dim system Ax = b, Cx <= d, with m equality constraints, l inequality constraints, and a ball of radius r centered at z, for any r this finds either a hyperplane separating the ball and the feasible region, described by h and delta with flag and delta set to 0, or returns an error-approximate solution x with flag set to 1 and fail set to 0.  This is called from the function Chubanov_Preliminary.  The parameter c_max is determined by the function find_c_max and is described below.  For a full explanation of how this recursive function works, see the papers referenced above.

The parameter alpha is a vector that keeps track of the inequality constraints used in the current induced hyperplane.  If a separating hyperplane is returned to the function Chubanov_Preliminary, then the values in alpha determine what happens next.  The vector fail_alpha fulfills the same role if the function fail_test returns a  positive result.

At the beginning of every recursive call of this function, we check if it has exceeded the time limit allotted.  If so, time_fail is set to 1 (it is always 0 otherwise), and we immediately quit out of the algorithm with every other returned value set to 0. 

When called from Chubanov_Preliminary, the parameters of this function and all the functions that follow are populated automatically.

******************************************

EP.m

[x,h,delta,alpha,flag] = EP(n,m,A,b,C,d,z,r,alpha)

Given an n-dim system Ax = b, Cx <= d, with m equality constraints, l inequality constraints, and a ball of radius r centered at z, where r is small enough as determined in the D&C, this function either finds a supporting hyperplane of the feasible region that separates the ball from the feasible region, or determines that the ball contains an error-approximate solution, where the error is determined by the user.  The D&C determines when r is small enough to enter this function.

If a separating hyperplane (h, delta) is found, we return h and delta as well as flag = 0 and the zero vector for x.  If a solution x is found, we return x and flag = 1, and the zero vector for h and delta = 0.  The vector alpha keeps track of which inequality constraint is used for h, if one is used.

******************************************

fail_test.m

fail = fail_test(h1,h2)

Once the D&C has found two separating hyperplanes, and before it determines the best induced separating hyperplane, we need to check if h1 = -gamma * h2, for gamma > 0.  If this happens, then the system is infeasible and we need to quit the recursive D&C.  Given hyperplanes h1 and h2 in n-dim space, this function checks for such a gamma.  If one is found, we return fail = 1, and otherwise return fail = 0.

******************************************

find_alpha.m

[alpha] = find_alpha(h1,delta1,h2,delta2,z,r)

This is the final step in each call to the D&C algorithm.  Once the two different separating hyperplanes (h1,delta1) and (h2,delta2) are determined, this function calculates the alpha that produces the best induced separating hyperplane (h = alpha*h1 + (1 - alpha)*h2, delta = alpha*delta1 + (1 - alpha)*delta2).

******************************************

equiv_equals.m

[A,b,r] = equiv_equals(A,b)

Given a linear system Ax = b, this function produces an equivalent full-rank system.  It returns the new matrix A, vector b, and r the rank of A.  It uses fref.m to determine which rows of A can be kept and which should be discarded in the creation of the full-rank matrix.

******************************************

fref.m

[A,v] = fref(A,tol)

The algorithm is similar to MATLAB's RREF, but uses vectorized code for increased efficiency, and does not scale the rows to produce leading ones.  Instead it simply does the first few steps of the RREF function, creating rank(A) linearly independent rows followed by rows of zeros.  If the parameter is not set, it uses MATLAB's default tolerance in the computations.  Otherwise it uses the tolerance tol.

It is called by equiv_equals.m as [~,v] = fref(A).  The purpose of this is not to actually create a row-reduced form of A, but simply to determine which rows need be removed from A in order to create a matrix with full rank.  Then last m - rank(A) entries in v precisely tell us these rows that need to be removed.  In this manner we create a full-rank matrix with integer entries much faster than the original rref function provided by MATLAB.