function x = obils(B,y,l,u,cp)
% 
% x = obils(B,y,l,u,cp) produces the solution to the overdetermined box
% constrained integer least squares problem min_{l<=x<=u}||y-Bx||.
%
% Inputs:
%    B - m-by-n real matrix with full column rank
%    y - m-dimensional real vector
%    l - n-dimensional integer vector, lower bound 
%    u - n-dimensional integer vector, upper bound, l < u   
%    cp - column permutation option: 
%         cp = 0 (default), perform no column permutations;
%         cp = 1, perform column permutations, preferable for large
%                 resdiual cases.
%
% Output:
%    x - n-dimensional integer vector (in double precision) for
%        the optimal solution
%

% Subfunctions: obils_reduction, obils_search

% Main References: 
% [1] S. Breen and X.-W. Chang. Column Reordering for 
%     Box-Constrained Integer Least Squares Problems, 
%     Proceedings of IEEE GLOBECOM 2011, 6 pages.
% [2] X.-W. Chang and Q. Han. Solving Box-Constrained Integer Least 
%     Squares Problems, IEEE Transactions on Wireless Communications,  
%     7 (2008), pp. 277-287.

% Authors: Xiao-Wen Chang, www.cs.mcgill.ca/~chang 
%          Xiangyu Ren
% Copyright (c) 2015-2022. Scientific Computing Lab, McGill University.
% Last revision: December 2022

    % Check the input arguments
    if nargin > 5
        error('Too many input arguments');
    end
    
    if nargin < 4
        error('Not enought input arguments');
    end
    
    if nargin == 4 
       cp = 0;
    end
    
    [m,n] = size(B);
    
    if m ~= size(y,1) || size(y,2) ~= 1 || ...
          n ~= size(l,1) || size(l,2) ~= 1 || ...
          n ~= size(u,1) || size(u,2) ~= 1     % Input error
        error('Input arguments have a matrix dimension error!')
    end
    
    if rank(B) < n
        error('Matrix does not have full column rank, use ubils')
    end
    
    l = ceil(l); u = floor(u);  % To make it work with real bounds
    for i = 1 : n
        if l(i) >= u(i)
            error('Invalid upper bound or lower bound');
        end
    end
    
    if cp == 0
        
        % Reduction with no column permutations
        U = qr([B,y]);
        R = triu(U(1:n,1:n)); 
        y = U(1:n,n+1);
        
        % Search
        x = obils_search(R,y,l,u);
            
    else
        
        % Reduction with column permutations
        [R,y,l,u,p] = obils_reduction(B,y,l,u);
        
        % Search
        z = obils_search(R,y,l,u);
        
        % Reorder z to obtain the optimal solution
        x = zeros(n,1);
        for i = 1 : n
            x(p(i)) = z(i);
        end
            
    end
end

%========================================================================
%======================= Subfunction of obils ===========================
%========================================================================

function [R,y,l,u,p] = obils_reduction(B,y,l,u)
%
% [R,y,l,u,p] = obils_reduction(B,y,l,u) reduces the general overdetermined 
% box-constrained integer least squares problem to an upper triangular one 
% by the QR factorization with column permutations 
% Q'*B*P = [R; 0]. The orthogonal matrix Q is not produced.  
%
% Inputs:
%    B - m-by-n real matrix with full column rank
%    y - m-dimensional real vector
%    l - n-dimensional integer vector, lower bound 
%    u - n-dimensional integer vector, upper bound
%
% Outputs:
%    R - n-by-n real nonsingular upper triangular matrix
%    y - n-dimensional vector transformed from the input y, y:=(Q'*y)(1:n)
%    l - permuted input lower bound l, i.e., l := P'*l 
%    u - permuted input upper bound u, i.e., u := P'*u
%    p - n-dimensional permutation vector representing P

% Main Reference: 
% S. Breen and X.-W. Chang. Column Reordering for 
% Box-Constrained Integer Least Squares Problems, 
% Proceedings of IEEE GLOBECOM 2011, 6 pages.

% Authors: Xiao-Wen Chang, www.cs.mcgill.ca/~chang
%          Xiangyu Ren
% Copyright (c) 2015-2022. Scientific Computing Lab, McGill University.
% Last revision: December 2022

    [~,n] = size(B);
    
    % Transform B and y by the QR factorization 
    U = qr([B,y]);
    R = triu(U(1:n,1:n));
    y = U(1:n,n+1);
    R0 = R;
    y0 = y;
    
    % Permutation vector 
    p = 1:n;
    
    % Inverse transpose of R
    G = inv(R)';
    
    % Determine the column permutatons 
    for k = n : -1 : 2
        maxDist = -1;
        
        % Determine the k-th column
        for i = 1:k
            alpha = y(i:k)' * G(i:k,i);
            x_i = max(min(round(alpha),u(i)),l(i));
            if (alpha < l(i) || alpha > u(i) || alpha == x_i)
                dist = 1 + abs(alpha - x_i); 
            else 
                dist = 1 - abs(alpha - x_i);
            end
            dist_i = dist / norm(G(i:k,i));
            if dist_i > maxDist
                maxDist = dist_i;
                j = i;
                x_j = x_i;
            end
        end
    
        % Perform permutations
        p(j:k) = p([j+1:k,j]);
        l(j:k) = l([j+1:k,j]);
        u(j:k) = u([j+1:k,j]);
    
        % Update y, R and G for the new dimension-reduced problem
        y(1:k-1) = y(1:k-1) - R(1:k-1,j) * x_j;
        R(:,j) = [];
        G(:,j) = [];
    
        for t = j : k - 1
            % Triangularize R and G by Givens rotation        
            [W, R([t,t+1],t)] = planerot(R([t,t+1],t));
            R([t,t+1],t+1:k-1) = W * R([t,t+1],t+1:k-1);
            G([t,t+1],1:t) = W * G([t,t+1],1:t);
            % Apply the Givens rotation W to y
            y(t:t+1) = W * y(t:t+1);
        end
    end
    
    % Reorder the columns of R0 according to p  
    R0 = R0(:,p);
    
    % Transform R0 and y0 by the QR factorization
    U = qr([R0,y0]);
    R = triu(U(:,1:n));
    y = U(:,n+1);
end

%========================================================================
%======================= Subfunction of obils ===========================
%========================================================================

function zhat = obils_search(R,y,l,u)
%
% zhat = obils_search(R,y,l,u,beta) produces the optimal solution to 
% the upper triangular box-constrained integer least squares problem
% min_{z}||y-Rz|| s.t. z in [l, u] by a search algorithm.

% Inputs:
%    R - n by n real nonsingular upper triangular matrix
%    y - n-dimensional real vector
%    l - n-dimensional integer vector, lower bound 
%    u - n-dimensional integer vector, upper bound    
%
% Outputs:
%    zhat - n-dimensional integer vector (in double precision). 

% Subfunctions: init, update (included in this file)

% Main Reference: 
% X.-W. Chang and Q. Han. Solving Box-Constrained Integer Least 
% Squares Problems, IEEE Transactions on Wireless Communications,  
% 7 (2008), pp. 277-287.


% Authors: Xiao-Wen Chang, www.cs.mcgill.ca/~chang
%          Xiangyu Ren, Jiqqun Shen
% Copyright (c) 2015-2022. Scientific Computing Lab, McGill University.
% April 2015. Last revision: December 2022

    % ------------------------------------------------------------------
    % --------  Initialization  ----------------------------------------
    % ------------------------------------------------------------------
    
    n = size(R,1);
    
    % Current point
    z = zeros(n,1);
    
    % c(k)=(y(k)-R(k,k+1:n)*z(k+1:n))/R(k,k)
    c = zeros(n,1);
    
    % d(k): left or right search direction at level k   
    d = zeros(n,1);
    
    % lflag(k) = 1 if the lower bound is reached at level k
    lflag = zeros(size(l));
    % uflag(k) = 1 if the upper bound is reached at level k
    uflag = lflag;
    
    % Partial squared residual norm for z
    % prsd(k) = (norm(y(k+1:n)-R(k+1:n,k+1:n)*z(k+1:n)))^2
    prsd = zeros(n,1);
    
    % Store some quantities for efficiently calculating c
    % S(k,n) = y(k),
    % S(k,j-1) = y(k) - R(k,j:n)*z(j:n) = S(k,j) - R(k,j)*z(j), j=k+1:n
    S = zeros(n,n);
    S(:,n) = y;
    
    % path(k): record information for updating S(k,k:path(k)-1) 
    path = n*ones(n,1); 
    
    % The level at which search starts to move up to a higher level
    ulevel = 0;  
    
    % Initial search radius
    beta = Inf; 
    
    % ------------------------------------------------------------------
    % --------  Search process  ----------------------------------------
    % ------------------------------------------------------------------
    
    k = n;
    
    % Find the initial integer in [l(n), u(n)] 
    c(k) = S(k,k) / R(k,k);
    [z(k),d(k),lflag(k),uflag(k)] = init(c(k),l(k),u(k));
    gamma = R(k,k) * (c(k) - z(k));
    
    % dflag for down or up search direction 
    dflag = 1; % Intend to move down to a lower level
    
    
    while 1
        if dflag == 1
           % Temporary partial squared residual norm at level k
           newprsd = prsd(k) + gamma * gamma;
           
           if newprsd < beta  % Inside the ellispoid
               if k ~= 1 % Move to level k-1
                   % Update path  
                   if ulevel ~= 0 
                       path(ulevel:k-1) = k;
                       for j = ulevel-1 : -1 : 1
                           if path(j) < k
                               path(j) = k;
                           else
                               break;  % Note path(1:j-1) >= path(j)
                           end
                       end
                   end
                   
                   % Update S
                   k = k - 1;
                   for j = path(k):-1:k+1
                       S(k,j-1) = S(k,j) - R(k,j)*z(j);
                   end
                   
                   % Update the partial squared residual norm
                   prsd(k) = newprsd;
                   
                   % Find the initial integer in [l(k), u(k)]
                   c(k) = S(k,k) / R(k,k);
                   [z(k),d(k),lflag(k),uflag(k)] = init(c(k),l(k),u(k));
                   gamma = R(k,k) * (c(k) - z(k));
                   
               else    % A valid point is found, update search radius
                   zhat = z;
                   beta = newprsd;
               end          
               ulevel = 0;           
           else           
               dflag = 0;  % Will move back to a higher level           
           end       
        else
           % Outside the ellispoid
           if k == n % The optimal solution has been found, terminate
               break;
           else
           % Move back to level k+1
               if ulevel == 0
                   ulevel = k;
               end
               k = k + 1;
               if lflag(k) ~= 1 || uflag(k) ~= 1
                   % Find a new integer at level k  
                   [z(k),d(k),lflag(k),uflag(k)] = ...
                          update(z(k),d(k),lflag(k),uflag(k),l(k),u(k));
                   gamma = R(k,k) * (c(k) - z(k)); 
                   dflag = 1;
               end
           end
        end
    end
end

%========================================================================
%==================== Subfunction of obils_search =======================
%========================================================================

function [z_k,d_k,lflag_k,uflag_k] = init(c_k,l_k,u_k)
%
% Find the initial integer and the search direction at level _k 

    z_k = round(c_k);
    if z_k <= l_k
        z_k = l_k;
        lflag_k = 1;  % The lower bound is reached
        uflag_k = 0;
        d_k = 1;
    elseif z_k >= u_k
        z_k = u_k;
        uflag_k = 1;  % The upper bound is reached
        lflag_k = 0;
        d_k = -1;
    else
        lflag_k = 0; 
        uflag_k = 0;
        if c_k > z_k
            d_k = 1;
        else
            d_k = -1;
        end
    end
end

%========================================================================
%==================== Subfunction of obils_search =======================
%========================================================================

function [z_k,d_k,lflag_k,uflag_k] = update(z_k,d_k,lflag_k,uflag_k,l_k,u_k)
%
% Find a new integer at level k and record it if it hits a boundary.
    z_k = z_k + d_k; 
    if z_k == l_k
        lflag_k = 1;
        d_k = -d_k + 1;
    elseif z_k == u_k
        uflag_k = 1;
        d_k = -d_k - 1;
    elseif lflag_k == 1
        d_k = 1;
    elseif uflag_k == 1
        d_k = -1;
    else
        if d_k > 0
            d_k = -d_k - 1;
        else
            d_k = -d_k + 1;
        end
    end
end