function ineq(inA)
% Code accompanying the paper:
%G. Valiant and P. Valiant. "An Automatic Inequality Prover and Instance
%Optimal Identity Testing", FOCS 2014.   Thanks to Clement Canonne for
%fixing a bug in our initial implementation (line 8)
%
% "inA" is a matrix with d+1 columns specifying the coefficients of the
% inequalities to be analyzed, with the last column specifying the powers
% of each term
%
% For example, if "inA" is a (m x 3) matrix, it represents the following inequality:
% \prod_{j=1}^m [\sum_i x_i^(inA(j,1)) y_i^(inA(j,2))]^(inA(j,3)) >= 1
% and the command "ineq(inA)" will determine whether the inequality holds,
% for all nonnegative sequences {x_i},{y_i}.
%
% If "inA" is a (m x 4) matrix, it represents the following inequality:
% \prod_{j=1}^m [\sum_i x_i^(inA(j,1)) y_i^(inA(j,2)) z_i^(inA(j,3))]^(inA(j,4)) >= 1
% and the command "ineq(inA)" will determine whether the inequality holds,
% for all nonnegative sequences {x_i},{y_i},{z_i}.
%
% The interpretation of "inA" of dimension (m x d+1) for d>3 is
% analogous.


% Example:  
% inA=[1,0,1/2; 0, 1, 1/2; 1/2,1/2,-1];
% represents Cauchy-Schwarz: (\sum_i x_i)^(1/2) (\sum_i y_i)^(1/2) (\sum_i x_i^(1/2)
% y_i^(1/2))^(-1) \ge 1.
% "ineq(inA)" yield the following output:
% > Optimization terminated.
% > Your inequality is true!
% > It can be derived as the product of 1 Holder inequalities:
% > Sum_j(a_j^1 b_j^-0 )^0.5 Sum_j(a_j^0 b_j^1 )^0.5  >= Sum_j(a_j^0.5 b_j^0.5 )^1   <-- all to the power 1
% > and 0 Lp monotonicity inequalities:
  
d=size(inA,2)-1; %dimension (number of sequences)
r=size(inA,1); %number of terms in the inequality
A=zeros(r,0); %constraint matrix, to be filled in
Ai=[]; %bookkeeping; records i in an entry if the corresponding constraint contains an Lp monotonicity inequality involving the ith point
for i=1:r %look for all simpexes that contain the ith point
    t=1:r; %t is candidates
    c=inA(:,1:d)'; c(:,i)=c(:,i)*2; c=c-inA(i*ones(r,1),1:d)';%coordinates; ith entry is the origin reflected about the point under consideration (see paper)
    for j=1:d+1 %loop over possible dimensions of simplexes
        tn=[]; %tn is candidates for the next dimension of simplexes
        for k=t %loop over the candidates
            m=null(c(:,k)); 
            if (size(m,2)==1) && all(m/m(1)>0), %see if the candidate simplex contains our point
                m=m/m(1); m=m/sum(m); %normalize so m sums to 1
                A(k,end+1)=-m; %add a constraint
                Ai(end+1)=i*(A(i,end)~=0); %bookkeeping form Lp monotonicity inequalities
                A(i,end)=1+2*A(i,end); %deals with both kinds of constraints on the ith point
            else %we add points to a simplex only if it doesn't already contain the desired point
                n=setdiff(1:r,k); %new points which could be added to the simplex
                tn(1:j+1,end+1:end+r-j)=[k(:,ones(r-j,1)); n]; %add the simplices to the list
            end
        end
        [t,ind]=unique(sortrows(sort(tn',2)),'rows','last');
        ind=diff([0;ind]);
        t=t(ind==j+1,:)';
        %keep only those simplices that appear j+1 times (i.e., where none
        %of their faces already contain the desired point
    end
end
A=A'; b=zeros(size(A,1),1);
[x1,x2,x3,x4,x5]=linprog(inA(:,end),A,b,[],[],[],[],[],optimset('Algorithm','interior-point'));%optimset('LargeScale','Off','Simplex','On'));
%NOTE on the previous line: the simplex algorithm may find a shorter (but equally correct) derivation of
% your inequality, but is being removed from future versions of Matlab. To use the simplex algorithm in 
% supported versions of Matlab, change 'interior-point' to 'simplex' on the previous line, or, on pre-2012 
% Matlab versions replace the optimset call with optimset('LargeScale','Off','Simplex','On')
if x3==1 %inequality is true!
    Hc=zeros(2,d,0); Lc=[]; %to be filled in with the coordinates of the Holder and Lp monotonicity inequalities that we use
    Hl=[]; Ll=[]; %to be filled in with the lambdas of the inequalities
    Hp=[]; Lp=[]; %to be filled in with the powers of the inequalities
    w=find(abs(x5.ineqlin)>1e-8); %which constraints are tight
    for j=w'
        ind=setdiff(find(A(j,:)<0),Ai(j)); %indices of which points are involved in this simplex
        c=A(j,ind); %coefficients of how to express the ith point as a convex combination of the simplex's vertices
        p=bsxfun(@times,cumsum(c(ones(d,1),:)'.*inA(ind,1:d)),1./cumsum(c)'); %weighted centers of the simplex
        for k=2:length(c)
            Hc(:,:,end+1)=[p(k-1,:);inA(ind(k),1:d)];
            Hl(end+1)=c(k)/sum(c(1:k));
            Hp(end+1)=-sum(c(1:k))*x5.ineqlin(j);
        end
        if Ai(j), Lc(:,end+1)=inA(Ai(j),1:d); Ll(end+1)=2-2/(1+A(j,Ai(j))); Lp(end+1)=x5.ineqlin(j)*A(j,Ai(j)); end
    end
    fprintf('Your inequality is true!\n  It can be derived as the product of %d Holder inequalities:\n',size(Hc,3));
    for j=1:size(Hc,3),
        for k=1:3
            if k~=3, v=Hc(k,:,j); p=1-Hl(j); else v=[1-Hl(j) Hl(j)]*Hc(:,:,j); p=1; end %% All changes here
            %%%%%if k~=3, v=Hc(k,:,j); p=Hl(j); else v=[Hl(j) 1-Hl(j)]*Hc(:,:,j); p=1; end
            if k==2, p=1-p; end
            fprintf('Sum_j(');
            for a=1:d,
                fprintf(['a'+a-1 '_j^%g '],v(a));
            end
            fprintf(')^%g ',p);
            if k==2, fprintf(' >= '); end
        end
        fprintf('  <-- all to the power %g\n',Hp(j));
    end
    fprintf('  and %d Lp monotonicity inequalities:\n',size(Ll,2));
    for j=1:size(Ll,2),
        for k=1:2
            if k==2, v=Lc(:,j); p=Ll(j); else v=Lc(:,j)*Ll(j); p=1; end
            fprintf('Sum_j(');
            for a=1:d,
                fprintf(['a'+a-1 '_j^%g '],v(a));
            end
            fprintf(')^%g ',p);
            if k==1, fprintf(' >= '); end
        end
        fprintf('  <-- all to the power %g\n',Lp(j)/Ll(j));
    end
else fprintf('Your inequality is false!\n'); end
