function z = ILSQ-PPS(A,b,nu)
% 
%   The PPS-algorithm, implementing parameter partitioning method, for computing 
%   outer component-wise estimates of the least squares solution set to interval 
%   system of linear algebraic equations. ILSQ_PPS applies the PPS-method to the 
%   extended form of the normal equations system, and the detailed theory of 
%   the approach is described in the article 
%
%           Shary, S.P., Moradi, B. Solving interval linear least squares 
%           problems by PPS-methods  //  Numerical Algorithms  (2020). 
%           DOI: 10.1007/s11075-020-00958-x 
%  
%   The code is designed for estimating the minimal value of the nu-th component 
%   for the points from the least-squares solution set. The maximal component-wise 
%   estimates can be found by negating the minimum of the same component for the 
%   system with the right-hand side vector taken with 'minus' sign. 
%  
%   As a basic enclosure method, ILSQ_PPS  uses  Hansen-Bliek-Rohn procedure 
%   (function 'HanBliRo') in interpretation by A.Neumaier. Hansen-Bliek-Rohn 
%   procedure is supplemented by midpoint inverse preconditioning. Interval 
%   Gauss method is used to construct an initial enclosure of the solution set. 
% 
%   Input argumetns:
%       A  is the interval matrix of the system under solution, 
%       b  is the interval right-hand side vector of the system, 
%       nu is the number of the component along which the solution set 
%             should be estimated. 
%   Output argument:
%       the minimal value of the nu-th component for the points 
%       from the least-squares solution set. 
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  MaxStepNumber = 80;   %   limitation on the number of steps to be executed 
                        %   by the algorithm;  adjust it, if necessary  
  MaxUEImprovs  = 4;    %   number of changes of the estimate UpEst, after 
                        %   which we clean the working list                         
  Eps = 1.e-4;          %   tolerance on the accuracy of the estimate 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Checking the correctness of the input data  
%  
    m = size(A,1); 
    n = size(A,2); 
    k = size(b,1);
    l = size(b,2);
    if ( k~=m | l>1 ) 
        error('Right-hand side size does no match to that of the matrix!')
    endif 
   
    if ( nu < 1 | nu > n ) 
        error('Estimated component number is invalid!')
    endif 
   
%   Constructing interval extended form of the normal equations system; 
%   its matrix is called AA and its right-hand side vector is called bb 
    AA = [ eye(m) A; transpose(A) zeros(n)]; 
    bb = [ b; zeros(n,1)]; 
    AA = infsup(AA); 
    bb = infsup(bb); 
    N = m + n; 
    
%   adjusting the number of the component to be estimated taking into account 
%   the extended form of the linear system of normal equations 
    nu = nu + m; 
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  
%   computing an initial enclosure of the solution set 
%
    x =  IGauss(AA,bb);    
   
%    x  = HanBliRo(AA,bb);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   computing an initial enclosure Z of the inverse interval matrix 
  
Z = infsup(zeros(N)); 
for j = 1:N; 
    E = infsup(eye(N));
    Z(:,j) = HanBliRo(AA,E(:,j)); 
end 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    
%   initialize the working list of the algorithm 
    
WorkList(1) = struct( 'Matrix', AA, 'RHSide', bb, ... 
                      'InvIMatr', Z, 'Encl', x, 'Estim', inf(x(nu)) );
                            
%   launching the PPS-algorithm proper  
NStep = 1;       
nuEst = 1; 
UpEst = realmax; 
UECounter = 0; 
  
while ( NStep < MaxStepNumber & (UpEst-nuEst) > Eps ) 
    
    %   searching the leading interval system in WorkList 
    LL = size(WorkList,2); 
    nuEst = realmax; 
    for k = 1:LL; 
        p = WorkList(k).Estim;
        if p < nuEst
            nuEst = p;
            Leading = k; 
        endif
    end
    
    %   the system Qx = r is set equal to the leading interval system 
    Q = WorkList(Leading).Matrix; 
    r = WorkList(Leading).RHSide;  
    x = WorkList(Leading).Encl; 
    Z = WorkList(Leading).InvIMatr;  
    
    %   monotonicity testing, squeezing the interval elements if necessary       
    for i = 1:N; 
        for j = 1:N; 
            if rad(Q(i,j)) > 0  
                v = -Z(nu,i)*x(j) - Z(nu,j)*x(i); 
                iv = inf(v); sv = sup(v); 
                if iv*sv >= 0 
                    if iv>=0 && sv>=0 
                        Q(i,j) = inf(Q(i,j)); 
                    endif 
                    if iv<0 && sv<0 
                        Q(i,j) = sup(Q(i,j)); 
                    endif 
                endif 
            endif 
        end 
        if rad(r(i)) > 0 
            ir = inf(Z(nu,i)); sr = sup(Z(nu,i)); 
            if ir*sr >= 0 
                if ir>=0 && sr>=0 
                    r(i) = inf(r(i)); 
                endif  
                if ir<0 && sr<0 
                    r(i) = sup(r(i)); 
                endif 
            endif
        endif 
    end 
    
    %   outputting intermediate results; 
    %   mute it if necessary
    if rem(NStep,3) == 0 
        display('==========================================================') 
        display(NStep),display(LL), display(Q)
        display(r)
        display(nuEst)
    endif 
    
    %   bisecting the leading interval system to descendants 
    
    %   first, we assign the previous values to the new matrises 
    %   and right-hand side vectors 
    Matr1 = Q;  Matr2 = Matr1; 
    b1 = r;     b2 = b1;
    
    %   next, we find the element within the matrix and right-hand-side 
    %   for which the product of its width by the estimate of the derivative 
    %   of the objective function is maximal 
    Amaxrad = -realmax; 
    bmaxrad = -realmax;     
    for i = 1:N; 
        for j= 1:N;
            if rad(Q(i,j)) > 0  
                qq = -Z(nu,i)*x(j) - Z(nu,j)*x(i);
                mm = mag(qq)*rad(Q(i,j)); 
                if Amaxrad < mm
                    Amaxrad = mm;
                    Ami = i; 
                    Amj = j;
                 endif 
            endif 
        end
        if rad(r(i)) > 0 
            rr = Z(nu,i);
            mm = mag(rr)*rad(r(i)); 
            if bmaxrad < mm; 
                bmaxrad = mm; 
                bmi = i; 
            endif 
        endif 
    end 
    
    %   if the maximal radii are zero, we stop the execution 
    if (Amaxrad == 0 && bmaxrad == 0) 
        display('Matrix subdividing finished'), display(Q)
        display(r)
        break;
    endif 
    
    %   if the widest element is in the matrix, we subdivide the matrix 
    if Amaxrad >= bmaxrad 
        s = Q(Ami,Amj); 
        %   the elements of the matrix are divided to their halves 
        Matr1(Ami,Amj) = infsup(inf(s),mid(s));  Matr1(Amj,Ami) = infsup(inf(s),mid(s)); 
        Matr2(Ami,Amj) = infsup(mid(s),sup(s));  Matr2(Amj,Ami) = infsup(mid(s),sup(s));
    else 
    %   if the widest element is in the right-hand side vector, we subdivide it 
        s = r(bmi); 
        %   the components of the right-hand side are divided to endpoints 
        b1(bmi) = infsup(inf(s),inf(s)); 
        b2(bmi) = infsup(sup(s),sup(s));          
    endif 
        
    %   compute enclosures of the solution sets for systems-descendants 
    
        %x1 = IGauss(Matr1,b1);              x2 = IGauss(Matr2,b2);
    
        x1 = HanBliRo(Matr1,b1);            x2 = HanBliRo(Matr2,b2);
        w1 = inf(x1(nu));                   w2 = inf(x2(nu)); 
                
        %   compute the upper estimates of the minimum 
        %                 and then correct UpEst if necessary 
        px1 = HanBliRo(mid(Matr1),mid(b1));     pw1 = sup(px1(nu)); 
        px2 = HanBliRo(mid(Matr2),mid(b2));     pw2 = sup(px2(nu)); 
        pw = min([pw1,pw2]); 
        
        %   fixation of the improvement of UpEst: 
        %   if sufficiently many improvements occurs, then 
        %   we clean the working list from unpromising records 
        if UpEst > pw 
            UpEst = pw; 
            UECounter = UECounter + 1; 
        endif 
        lamp_on = 0;
        lamp_off = 0;
        if UECounter > MaxUEImprovs  
            UECounter = 0;
            ENum = 1; 
            while (lamp_on + lamp_off < LL)
                if WorkList(ENum).Estim >= UpEst 
                    WorkList(ENum) = [];
                    LL = size(WorkList,2);
                    lamp_on = lamp_on + 1 ;   
                else 
                    lamp_off = lamp_off + 1;
                endif 
                ENum = ENum + 1;
            end       
        endif 
    
    %   recalculating enclosures of the inverse interval matrices 
    Z1 = infsup(zeros(N));     Z2 = infsup(zeros(N)); 
    for l = 1:N; 
        E  = infsup(eye(N)); 
        Z1(:,l) = HanBliRo(Matr1,E(:,l)); 
    end        
    for l = 1:N; 
        E  = infsup(eye(N)); 
        Z2(:,l) = HanBliRo(Matr2,E(:,l)); 
    end            
    
    %   construct the records-descendants and put them into WorkList 
    Desc1 = struct( 'Matrix', Matr1, 'RHSide', b1, ...
                    'InvIMatr', Z1, 'Encl', x1, 'Estim', w1 );
    Desc2 = struct( 'Matrix', Matr2, 'RHSide', b2, ... 
                    'InvIMatr', Z2, 'Encl', x2, 'Estim', w2 ); 
    if w1 <= UpEst
        WorkList = [ WorkList, Desc1 ];
    endif 
    if w2 <= UpEst
        WorkList = [ WorkList, Desc2 ]; 
    endif 
    
    %   deleting the former leading record from WorkList 
    WorkList(Leading) = []; 
    NStep = NStep + 1; 
    
    if NStep >= MaxStepNumber 
        display('-- The limitation on the number of steps is reached! --') 
    endif 
end 
  
display('End of execution: ###################################################')
display(NStep) 
display('Final divided matrix: '), display(Q)
display('Final divided right-hand side: '), display(r)
z = nuEst; 

display('Upper estimate on the error: UpEst-nuEst = ')
display(UpEst - nuEst), display('') 
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
function x = IGauss(A,b)
    %   
    %   interval Gauss method with partial pivoting 
    %    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %   direct sweep of interval Gauss method 
        for j = 1:N-1;
            %   finding the pivot element, i.e. having the largest mignitude,  
            %   in the lower diagonal part of the current column 
            [maxmig,index] = max(mig(A(j:N,j)));
            if maxmig <= 0 
                display('All the lower diagonal pivots contain zeros!')
                error('The interval matrix of the problem is too wide.')
            endif
            k = index + j - 1;
            if k ~= j      %   swapping the rows, if necessary 
                A([j,k],j:N) = A([k,j],j:N);
                b([j,k]) = b([k,j]);
            endif
            %    now, the transformation of the forward sweep of Gauss mewthod 
            %    are performed 
            for i = j+1:N;
                r = A(i,j)/A(j,j);
                A(i,j+1:N) = A(i,j+1:N) - r*A(j,j+1:N);
                b(i) = b(i) - r*b(j);
            end
        end
    
        %   backward substitution for computing 
        %   the interval enclosures of the components 
        x(N) = b(N)/A(N,N); 
    
        for i = (N-1):(-1):1;
            r = b(i);
            for j = i+1:N;
                r = r - A(i,j)*x(j);
            end
            x(i) = r/A(i,i);
        end 
   
endfunction 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
function x = HanBliRo(A,b) 
%   
%   Hanse-Bliek-Rohn procedure for computing an enclosure of the united 
%   solution set for interval linear systems with H-matrices 
% 
%   Input arguments 
%       A - interval matrix of the equation system 
%       b - interval right-hand side vector 
% 
%   If the interval matrix A is not numerically an H-matrix then the answer 
%   is set to [-inf,+inf]. Our implementation follows mainly that presented 
%   in the article 
%     Neumaier A. A simple derivation of Hansen-Bliek-Rohn-Ning-Kearfott
%     enclosure for linear interval equations // Reliable Computing. -
%     1999. - Vol. 5, No. 2. - P. 131-136.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   midpoint-inverse preconditioning 

    N = size(A,1); 
    A = infsup(A);     b = infsup(b); 
    Ainv = inv(mid(A)); 
    C = Ainv*A;        d = Ainv*b; 
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Computing the comparison matrix 
    Cc = -mag(C); 
    for i = 1:N; 
        Cc(i,i) = mig(C(i,i)); 
    end 
  
%   extract the diagonal from the matrix, 
%   inverse the comparison matrix, etc. 
    dC = diag(C);
    E = inv(Cc);
    v = abs(E*ones(N,1)); 
    __setround__(-Inf); 
    u = Cc*v;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   executing Hansen-Bliek-Rohn procedure proper 
  
    if ~all(min(u)>0) 
        %   C is not an H-matrix 
        x = infsup(-Inf,Inf)*ones(N,1); 
        display('The interval matrix of the system is too wide!') 
    else 
        %   C is an H-matrix 
        dCc = diag(Cc); 
        Cc = Cc*E - eye(N); 
        __setround__(+Inf); 
        w = zeros(1,N);
        for i = 1:N;
            w = max(w,(-Cc(i,:))/u(i));
        end; 
        dlow = v.*w' - diag(E); 
        dlow = max(0,-dlow); 
        E = E + v*w; 
        u = E*mag(d); 
        dE = diag(E); 
        alpha = dCc + (-1)./dE; 
        beta = u./dlow - mag(d); 
        x = (d + infsup(-beta,beta))./(dC + infsup(-alpha,alpha));
    endif 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   setting the old rounding mode 
  
    __setround__(0.5)
   
endfunction
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
endfunction 
  
  