function [S,d,dd]=regising(Ac,Delta,t,q) % REGularIty/SINGularity of an interval matrix % VERSION TO BE POSTED
%
%  Most simple use (if you are interested in result only, not in computational details):
%  [S]=regising(Ac,Delta)
%
%  A=[Ac-Delta,Ac+Delta] is a square interval matrix,
%  t==1:        displays number of orthants to be yet visited,
%  t~=1:        no display,
%  q==1:        plots number of orthants to be yet processed vs. iteration count,
%  q~=1:        no plot,
%  isempty(S):  A is regular (i.e., each B in A is nonsingular),
% ~isempty(S):  S is a singular matrix in A, possibly in normal form (if found),
%               i.e., abs(S-Ac)=Delta holds for each but at most one entry.
%  d:           total number of orthants visited by the main algorithm (it may be zero), 
%  dd:          symbol of the method which succeeded in determining singularity/regularity:
%               0.0: singularity of  the midpoint matrix Ac, 
%               0.1: singularity via diagonal condition,
%               0.2: singularity via steepest determinant descent,
%               0.3: singularity as  a by-product of 0.7,
%               0.4: singularity via the main algorithm,
%               0.5: regularity  via Beeck's condition,
%               0.6: regularity  via symmetrization,
%               0.7: regularity  via two Qz-matrices,
%               0.8: regularity  via the main algorithm.
%
%  MIT License for REGISING:
% 
%  Copyright (c) 2005-2016 Jiri Rohn
% 
%  Permission is hereby granted, free of charge, to any person obtaining a copy
%  of this software and associated documentation files (the "Software"), to deal
%  in the Software without restriction, including without limitation the rights
%  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
%  copies of the Software, and to permit persons to whom the Software is
%  furnished to do so, subject to the following conditions:
% 
%  The above copyright notice and this permission notice shall be included in all
%  copies or substantial portions of the Software.
% 
%  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
%  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
%  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
%  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
%  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
%  SOFTWARE.
%
n=size(Ac,1); d=0; dd=0;
if rank(Ac)<n, S=Ac; S=normform(Ac,Delta,S); dd=0.0; return, end
Aci=inv(Ac);
D=Delta*abs(Aci);
[dg,J]=max(diag(D)); j=min(J);
if dg>=1
    x=Aci(:,j); S=vec2mat(Ac,Delta,x); S=normform(Ac,Delta,S); dd=0.1; return
end
if beeckcond(Ac,Delta)<1, S=[]; dd=0.5; return, end
S=singdetdesc(Ac,Delta);
if ~isempty(S), S=normform(Ac,Delta,S); dd=0.2; return, end
AA=Ac'*Ac; DD=Delta'*Delta;
if rank(AA)==n&&beeckcond(AA,DD)<1, S=[]; dd=0.6; return, end
[rs,S]=regsuffcondqz(Ac,Delta);
if rs==0, S=normform(Ac,Delta,S); dd=0.3; return, end
if rs==1, dd=0.7; return, end
bc=janssonsheuristic(Aci);
ep=max(n,10)*max([norm(Ac-Delta,1) norm(Ac+Delta,1) norm(bc,1)])*eps;
Z=(sgn(Aci*bc))'; D=[]; 
if (nargin>=3)&&(isequal(q,1)), yy(1)=0; yy(2)=1; ii=2; end
while ~isempty(Z)
    p=size(Z,1); z=Z(p,:); Z=Z(1:p-1,:); D=[D; z];
    [Q,S]=qzmatrix(Ac,Delta,z);
    if ~isempty(S)
        S=normform(Ac,Delta,S); d=size(D,1); dd=0.4;
        if (nargin>=3)&&(isequal(q,1)), pass(yy), end 
        return
    end
    xut=Q*bc;
    if all(xut(z==1)>=-ep)
        [Q,S]=qzmatrix(Ac,Delta,-z);
        if ~isempty(S)
            S=normform(Ac,Delta,S); d=size(D,1); dd=0.4;
            if (nargin>=3)&&(isequal(q,1)), pass(yy), end
            return
        end
        xlt=Q*bc;
        if all(xlt<=xut)
            for j=1:n
                zt=z; zt(j)=-zt(j);
                if (xlt(j)*xut(j)<=ep)&&(~ismember(zt,Z,'rows'))&&(~ismember(zt,D,'rows'))
                    Z=[Z; zt];
                end
            end
        end
        if (nargin>=3)&&(isequal(q,1)), ii=ii+1; yy(ii)=size(Z,1); end
        if (nargin==4)&&(isequal(t,1)), disp(size(Z,1)), end
        if size(Z,1)+size(D,1)>=20*n^2
            if (nargin>=3)&&(isequal(q,1)), pass(yy), end; S=[]; d=[];
            error('Program run has been stopped after reaching prescribed number of iterations')
        end
    end
end
S=[]; d=size(D,1); dd=0.8; if (nargin>=3)&&(isequal(q,1)), pass(yy), end
%
function [Q,S]=qzmatrix(Ac,Delta,z) % $Q_z$ MATRIX
%
%  A=[Ac-Delta,Ac+Delta]: square interval matrix, z a plus/minus-one vector,
% ~isempty(Q): Q solves Q*Ac-abs(Q)*Delta*diag(z)=eye(size(A)), and S=[],
% ~isempty(S): S is a singular matrix in A, and Q=[].
%
n=size(Ac,1); I=eye(n,n); 
Q=zeros(n,n); S=[]; 
for i=1:n
    [x,S]=absvaleqn(Ac',-diag(z)*Delta',I(:,i));
    if ~isempty(S), Q=[]; S=S'; return, end
    Q(i,:)=x';
end
%
function [x,S,iter]=absvaleqn(A,B,b) % ABSolute VALue EQuatioN 
%
% ~isempty(x): x solves A*x+B*abs(x)=b (A,B square), S is empty,
% ~isempty(S): S is a singular matrix satisfying abs(S-A)<=abs(B), x is empty.
%  iter: number of iterations (it may be zero).
%
%  Copyright 2005-2015 Jiri Rohn (rohn@cs.cas.cz)
%
b=b(:); n=length(b); I=eye(n,n);
ep=n*(max([norm(A,inf) norm(B,inf) norm(b,inf)]))*eps;
x=[]; S=[]; iter=0;
if rank(A)<n, S=A; return, end
x=A\b; z=sgn(x);
if rank(A+B*diag(z))<n, S=A+B*diag(z); x=[]; return, end
x=(A+B*diag(z))\b;
C=-inv(A+B*diag(z))*B;
X=zeros(n,n);
r=zeros(1,n);
while any(z.*x<-ep)
    k=find(z.*x<-ep,1);
    iter=iter+1;
    if 1+2*z(k)*C(k,k)<=0
        S=A+B*(diag(z)+(1/C(k,k))*I(:,k)*I(k,:)); x=[]; return
    end
    if ((k<n)&&(all(r(k)>r(k+1:n))))||((k==n)&&(r(k)>0))
        x=x-X(:,k); z=sgn(x);
        ct=A*x;
        jm=abs(B)*abs(x);
        y=zeros(1,n);
        for i=1:n
            if jm(i)>ep, y(i)=ct(i)/jm(i); else y(i)=1; end
        end
        S=A-diag(y)*abs(B)*diag(z); x=[]; return
    end
    X(:,k)=x;
    r(k)=iter;
    z(k)=-z(k);
    alpha=2*z(k)/(1-2*z(k)*C(k,k));
    x=x+alpha*x(k)*C(:,k);
    C=C+alpha*C(:,k)*C(k,:);
end 
%
function [S]=vec2mat(Ac,Delta,x) % VECtor TO MATrix
%
% A=[Ac-Delta,Ac+Delta] is a square interval matrix, x is a matching 
% nonzero point vector satisfying abs(Ac*x) <= Delta*abs(x).
% S: a singular matrix in A.
%
n=length(x);
ep=max([n 100])*max([norm(Ac-Delta,inf) norm(Ac+Delta,inf) norm(x,inf)])*eps;
ct=Ac*x; jm=Delta*abs(x);
y=zeros(n,1); z=y;
for i=1:n
    if jm(i)>ep, y(i)=ct(i)/jm(i); else y(i)=1; end
    if x(i)>=0, z(i)=1; else z(i)=-1; end
end
S=Ac-diag(y)*Delta*diag(z);
%
function z=sgn(x) % SiGN vector of x (column)
n=length(x); 
z=zeros(n,1);
for j=1:n
    if x(j)>=0, z(j)=1; else z(j)=-1; end 
end
%
function [rs,S]=regsuffcondqz(Ac,Delta) % REGularity SUFFicient CONDition via matrices QZ
%
% A=[Ac-Delta,Ac+Delta] is a square interval matrix.
% rs= 1: A is regular (i.e., each B in A is nonsingular) and S=[], 
% rs= 0: A is singular and S is a singular matrix in A,
% rs=-1: no result.
% Uses simplexineq.m
%
n=size(Ac,1);
I=eye(n,n); e=ones(n,1); o=zeros(n,1);
ep=max(n,10)*max([norm(Ac-Delta,1) norm(Ac+Delta,1)])*eps;
Aci=inv(Ac);
[Q1,S]=qzmatrix(Ac,Delta,-e);
if ~isempty(S), rs=0; return, end
[Q2,S]=qzmatrix(Ac,Delta,e);
if ~isempty(S), rs=0; return, end
A=[ -Q1 e;
   -Aci e; 
      I o;
     -I o];
b=[o;
   o;
   e;
   e];
c=[o;
   1];
[x]=simplexineq(A,b,c); % max c'*x s.t. A*x<=b
if isinf(x(1)), rs=-1; S=[]; return, end
if x(n+1)>ep, rs=1; S=[]; return, end
rs=-1; S=[];
%
function S=singdetdesc(Ac,Delta,t) % SINGularity via DETerminant DESCent
%
% A=[Ac-Delta,Ac+Delta] is a square interval matrix,
% t==1: starts from Ac+Delta, otherwise from Ac.
%~isempty(S): S is a singular matrix in A, 
%             if t==1, then S is in a normal form, i.e.,
%             abs(S-Ac)=Delta holds for each but at most one entry.
% isempty(S): no result.
%
n=size(Ac,1);
Ad=Ac-Delta; Ah=Ac+Delta;
if rank(Ac)<n, S=Ac; return, end
if rank(Ah)<n, S=Ah; return, end
if nargin==3&&t==1, A=Ah; else A=Ac; end
C=inv(A);
beta=0.5;
p=zeros(n,1);
while beta<0.95 % avoid cycling
    Zd=C'>=0;
    Zh=C'<0;
    B=Zd.*Ad+Zh.*Ah;
    [beta,J]=min(diag(B*C)); k=min(J); 
    if abs(beta)>=0.95, S=[]; return, end
    if 1e-10<beta && beta<0.95
        D=B-A;
        C=C-(1/beta)*C(:,k)*(D(k,:)*C); % C=inv(A)
        A(k,:)=B(k,:);
    end
    if 0<beta && beta<=1e-10
        S=[]; return
    end
    if beta<=0
        for i=1:n
           p(i)=B(k,1:i)*C(1:i,k)+A(k,i+1:n)*C(i+1:n,k); 
        end
        if all(p>0), S=[]; return, end
        m=find(p<=0,1);
        if abs(C(m,k))<=1e-10, S=[]; return, end
        A(k,1:m-1)=B(k,1:m-1);
        A(k,m)=-(B(k,1:m-1)*C(1:m-1,k)+A(k,m+1:n)*C(m+1:n,k))/C(m,k);
        S=A; return
    end
end
S=[]; return
%
function pass(yy) % PASSing data to the plot fuction
% 
% Plots number of orthants to be yet processed vs. the iteration count
%
qq=length(yy);
xx=1:qq;
plot(xx,yy)
%
function bc=janssonsheuristic(Aci) % JANSSON'S HEURISTIC
%
% Heuristic for finding a right-hand side for the main algorithm.
%
n=size(Aci,1);
bc=ones(n,1);
g=min(abs(Aci*bc));
for i=1:n
    for j=1:n
        bp=bc; bp(j)=-bp(j);
        if min(abs(Aci*bp))>g, g=min(abs(Aci*bp)); bc=bp; end 
    end
end
for i=1:n
    for j=[1:i-1 i+1:n]
        bp=bc; bp(i)=-bp(i); bp(j)=-bp(j);
        if min(abs(Aci*bp))>g, g=min(abs(Aci*bp)); bc=bp; end 
    end
end
%
function S=normform(Ac,Delta,S) % NORMal FORM of a singular matrix
%
% Finds a normal form singular matrix (see the subfunction singdetdesc.m), 
% or stays with the one computed previously.
%
Snf=singdetdesc(Ac,Delta,1);
if ~isempty(Snf), S=Snf; end
%
function r=beeckcond(Ac,Delta) % BEECK's CONDition
%
% r< 1: [Ac-Delta,Ac+Delta] is regular,
% r>=1: no result.
%
r=max(abs(eig(abs(inv(Ac))*Delta)));
%
function [x,C]=simplexineq(A,b,c) % SIMPLEX method for INEQuality constraints
% SIMPLEXINEQ  Solving a linear program with inequality constraints.
%
% [x,C]=simplexineq(A,b,c) solves the linear programming problem
%     max c'*x s.t. Ax<=b
% by solving its dual problem
%     min b'*y s.t. A'*y=c, y>=0
% using simplex2016.p .
%
% Data: A in R^(mxn), b in R^m, c in R^n.
%
% Output variables:
%~isinf(x(1)) : x is an optimal solution,
% x(1)== Inf  : problem unbounded,
% x(1)==-Inf  : problem infeasible,
% C.y         : dual optimal solution,
% C.hp        : primal optimal value,
% C.hd        : dual optimal value,
% C.AS        : last dual simplex tableau,
% C.iter      : number of iterations of the simplex algorithm,
% C.ep        : tolerance used (numbers >= -C.ep considered nonnegative),
% C.flag      : verbal description of the output.
%
C.y=[]; C.hp=[]; C.hd=[]; C.AS=[]; C.iter=[]; C.ep=[]; C.flag=[]; % order
% main part
[y,CC]=simplex2016(A',c,b,0); % y is the dual optimum; primal is in CC.y
% assigning the output variables
% values independent of y
C.AS=CC.AS;
C.iter=CC.iter;
C.ep=CC.ep;
% values dependent on y
if  ~isinf(y(1))% dual optimum exists
    x=CC.y;     % primal optimum
    C.y=y;      % dual optimum
    C.hp=c'*x;  % primal optimal value
    C.hd=CC.h;  % dual optimal value
    C.flag='optimal solution';
end
if  y(1)== Inf  % dual problem infeasible; primal infeasible or unbounded
    x1=simplexeq([A -A eye(size(A,1))],b,ones(2*size(A,2)+size(A,1),1),0);
    if  ~isinf(x1(1)) % primal problem feasible, i.e. unbounded
        x=repmat(Inf,length(c),1);
        C.flag='unbounded';
    else % primal problem infeasible
        x=repmat(-Inf,length(c),1);
        C.flag='infeasible';
    end
end
if  y(1)==-Inf  % dual problem unbounded
    x=repmat(-Inf,length(c),1);
    C.flag='infeasible';
end
%
function [x,C]=simplex2016(A,b,c,d) % of 2012, main change in former line 114; p set to 2m (line 26)
% NOT TO BE USED SEPARATELY AS A LINEAR PROGRAMMING SOLVER
% COMPRISES FEATURES SPECIFIC FOR 'REGISING'
%
% Two-phase simplex method.
% Solves the problem min c'*x s.t. A*x=b, x>=0.
% d==0        : no screen output, 
% d==1        : screen output, objective value,
% d==2        : screen output, paused, simplex tableau,
%~isinf(x(1)) : x is an optimal solution,
% x(1)== Inf  : problem infeasible,
% x(1)==-Inf  : problem unbounded,
% C.w         : warning when basis improper,
% C.y         : dual optimal solution,
% C.h         : optimal value,
% C.B         : basis index set,
% C.Bout      : indices of artificial variables in the basis,
% C.N         : nonbasis index set,
% C.AS        : last simplex tableau,
% C.iter      : number of iterations,
% C.flag      : verbal description of the output,
% C.ep        : tolerance used (numbers >= -C.ep considered nonnegative).
%
% Overall construction: works only with ABi and Aaug, with phases  being 
% determined by l=m+n or l=n, algorithm the same in both cases.
%
b=b(:); c=c(:);
x=[]; y=[]; h=[]; B=[]; Bout=[]; N=[]; AS=[]; iter=[]; flag=[]; C.w='';
[m,n]=size(A); p=2*m; % restart of ABi after each p steps
ep=max([m n 100])*max([norm(A,inf) norm(b,inf) norm(c,inf)])*eps; % first term max([m n 100]) important
if ~(m==length(b)&&n==length(c)&&isreal(A)&&isreal(b)&&isreal(c))
    flag='wrong data';
    C.y=y; C.h=h; C.B=B; C.Bout=Bout; C.N=N; C.AS=AS; C.iter=iter; C.flag=flag; C.w=''; C.ep=ep;
    return
end
z=ones(m,1);
for i=1:m
    if b(i)<0, A(i,:)=-A(i,:); b(i)=-b(i); z(i)=-1; end % right-hand side >=0
end
I=eye(m,m);
B=(n+1:n+m)'; ABi=I; Aaug=[A I b]; % ABi=inv(Aaug(:,B)) throughout
phase=1; l=m+n; % phase 1 initialization
iter=0;
while isempty(flag)
    iter=iter+1;
    if phase==1 % repeated initialization
        crit=[zeros(1,n) ones(1,m)  0];
    else
        crit=[c'         zeros(1,m) 0];
    end
    y=crit(B)*ABi;
    crit=crit-y*Aaug; crit(B)=zeros(1,m); % current criterial row  
    if nargin==4&&d==1 
        -crit(m+n+1) % screen output of the current objective value
    end
    if nargin==4&&d==2
        [B ABi*Aaug; 0 crit], pause % paused screen output of the current simplex tableau
    end
    if all(crit(1:l)>=-ep) % termination
        if phase==1 % termination, phase 1
            if abs(crit(m+n+1))<=ep
                phase=2; l=n; % phase 2 initialization
            else % abs(crit(m+n+1))>ep```````````
                x=repmat(Inf,n,1); y=[]; flag='infeasible';
                C.y=y; C.h=-crit(m+n+1); C.B=B; C.Bout=B(B>n); C.N=N; C.AS=[B ABi*Aaug; 0 crit]; 
                C.iter=iter; C.flag=flag; C.ep=ep;
                if isempty(C.Bout), C.w=''; else C.w='improper basis'; end
                return
            end
        else % termination, phase 2
            x=zeros(n+m,1); x(B)=ABi*b; x=x(1:n); % optimal solution
            y=(-crit(n+1:n+m))'; % y=(ABi)'*c(B)
            y=z.*y; % y dual optimum for the original Ax=b
            flag='optimal solution';
            N=find(~ismember(1:n,B)); % nonbasis index set in 1:n
            if all(crit(N)>ep)
                flag='unique optimal solution';
            end
            C.y=y; C.h=c'*x; C.B=B; C.Bout=B(B>n); C.N=N'; C.AS=[B ABi*Aaug; 0 crit]; 
            C.iter=iter; C.flag=flag; C.ep=ep;
            if isempty(C.Bout), C.w=''; else C.w='improper basis'; end
            return
        end
    else % current step, any(crit(1:l)<-ep)
        s=find(crit(1:l)<-ep,1); % Bland's rule for s
        As=ABi*Aaug(:,s); bp=ABi*b; 
        if all(As<=ep)
            x=repmat(-Inf,n,1); y=[]; flag='unbounded';
            C.y=y; C.h=-Inf; C.B=B; C.Bout=B(B>n); C.N=N; C.AS=[B ABi*Aaug; 0 crit]; 
            C.iter=iter; C.flag=flag; C.ep=ep;
            if isempty(C.Bout), C.w=''; else C.w='improper basis'; end
            return
        else % any(As>ep)
            % Bland's rule for r
            if numel(As)==0, x=repmat(Inf,n,1); C=[]; return, end        
            for j=1:m
                if As(j)>ep, As(j)=bp(j)/As(j); else As(j)=Inf; end
            end
            mu=min(As); K=find(abs(As-mu*ones(m,1))<=ep);
            BB=Inf*ones(m,1); BB(K)=B(K); [nu,r]=min(BB);
            B(r)=s;
            % ABi update
            if iter~=p*round(iter/p) % iter not a multiple of p
                Aaugs=Aaug(:,s); ABir=ABi(r,:);
                a=(-1/(ABir*Aaugs))*(ABi*Aaugs-I(:,r));
                ABi=ABi+a*ABir; % Sherman-Morrison update
            else % iter is a multiple of p`% FORMER LINE 114 CAUSED CRASH ON 2012-11-22 (SEE BELOW)
                if rank(Aaug(:,B))==m % Aaug(:,B) invertible % IF ... END ADDED 2012-11-24
                    ABi=inv(Aaug(:,B)); % restart % FORMER LINE 114
                    if nargin==4&&d==1
                        display('restart')
                    end % of display of restart
                end % of restart
            end % of iter~=p*round(iter/p)
        end % of all(As<=ep)
    end % of all(crit(1:l)>=-ep)
end % of while
%
function [x,C]=simplexeq(A,b,c,d) % SIMPLEX method for EQuation constraints
%
% Two-phase revised simplex method.
% Solves the problem 
% min c'*x s.t. A*x=b, x>=0   (1).
%
% d==0        : no screen output, 
% d==1        : screen output, objective value,
% d==2        : screen output, paused, simplex tableau,
%~isinf(x(1)) : x is an optimal solution of (1),
% x(1)== Inf  : problem (1) infeasible,
% x(1)==-Inf  : problem (1) unbounded,
% C.w         : warning when basis improper,
% C.y         : dual optimal solution,
% C.h         : optimal value,
% C.B         : basis index set,
% C.Bout      : indices of artificial variables in the basis,
% C.N         : nonbasis index set,
% C.AS        : last simplex tableau,
% C.iter      : number of iterations,
% C.flag      : verbal description of the output,
% C.ep        : tolerance used (numbers >= -C.ep considered nonnegative).
%
[x,C]=simplex2016(A,b,c,d); % computation done by SIMPLEX2016