function [x,flag] = ubils(B,y,l,u,cp,maxtime)
% 
% x = ubils(B,y,l,u) produces the optimal solution to the underdetermined 
%     box-constrained integer least square problem min_{l<=x<=u}||y-Bx||.
%     B - m-by-n real matrix with rank(B) < n
%     y - m-dimensional real vector
%     l - n-dimensional integer vector, lower bound 
%     u - n-dimensional integer vector, upper bound, l < u
%     x - n-dimensional integer vector (in double precision), the optimal
%         solution 
%     
% x = ubils(B,y,l,u,cp) provides a column permutation option to choose
%     cp = 0 / [] (default), performs simple column permutations;
%     cp = 1, performs sophasitated column permutations, preferable for 
%             large resdiual cases. 
%
% x = ubils(B,y,l,u,cp,maxtime) specifies the maximum running time in seconds. 
%     This option is to avoid too long running time. 
%     The default value of maxtime is Inf. 
%     x is the optimal solution or a suboptimal solution.
%     
% [x,flag] = ubils(B,y,l,u,cp,maxtime) also returns a flag
%     flag = 0, x is the optimal solution
%     flag = 1, x is a suboptimal solution 

% Subfunctions: ubils_reduction, ubils_search

% Main References:
% [1] X.-W. Chang and X. Yang. An efficient tree search
%     decoder with column reordering for underdetermined MIMO systems. 
%     Proceedings of IEEE GLOBECOM 2007, 5 pages.
% [2] Jing Zhu. Numerical methods for underdetermined box-constrained 
%     integer least squares problems, M.Sc. Thesis, School of Computer 
%     Science, McGill University, 2016.

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

    % Check the input arguments
    if nargin > 6
        error('Too many input arguments');
    end
    
    if nargin < 4
        error('Not enought input arguments');
    end
    
    switch nargin
        case 4 
            cp = 0;
            maxtime = Inf;
        case 5 
            if isempty(cp)
               cp = 0;
            end
            if cp ~= 0 &&  cp ~= 1
                error('cp should be [], 0 or 1');
            end
            maxtime = Inf;
        case 6
            if isempty(cp)
               cp = 0;
            end
            if cp ~= 0 && cp ~= 1
                error('cp should be [], 0 or 1');
            end
            if maxtime <= 0
               error('Maxtime has to be positive');
            end
    end
    
    % Time limit (instance)
    duetime = cputime + maxtime;
        
    [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
    
    l = ceil(l); u = floor(u);  % To make it work with real bounds
    for i = 1 : n
        if l(i) >= u(i) || l(i) == -Inf || u(i) == Inf
            error('Invalid lower bound or upper bound');
        end
    end
       
    % If B is a zero matrix, one solution is obtained
    if B == 0
        x = l;
        flag = 0;
        warning('The input matrix is zero and here one solution is produced')
        return
    end
    
    % Reduction - reduce the problem to the trapezoid form
    [R,y,l,u,p,z0,beta] = ubils_reduction(B,y,l,u,cp);
    
    % Search 
    if beta > n*eps*norm(abs(y)+abs(R)*abs(z0))
        [z,flag] = ubils_search(R,y,l,u,z0,beta,duetime);   
    else
        z = z0; % z0 is the solution of the reduced problem
        flag = 0;
    end
    
    % Reorder z to obtain the solution to the original problem
    x = zeros(n,1);
    for i = 1 : n
        x(p(i)) = z(i);
    end 
end


%========================================================================
%====================== Subfunction of ubils ============================
%========================================================================

function [R,y,l,u,p,z0,beta] = ubils_reduction(B,y,l,u,cp)
%
% [R,y,l,u,p,beta] = ubils_reduction(B,y,l,u) reduces an underdetremined 
% problem or rank-deficient overdetermiend problem to an upper trapezoid 
% form 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 m<n
%    y - m-dimensional real vector
%    l - n-dimensional integer vector, lower bound 
%    u - n-dimensional integer vector, upper bound  
%    cp - column permutation option: 
%         cp = 0 (default), perform simple column permutations;
%         cp = 1, perform sophasitated column permutations, preferable for 
%                 large resdiual cases. 
%
% Outputs:
%    R - r by n real upper trapezoidal matrix, Q'*B*P = [R; 0]  
%        Q: m by m column orthogonal, P: n by n permutation matrix 
%    y - r-dimensional real vector, y := (Q'*y)(1:r) 
%    l - permued lower bound 
%    u - permuted upper bound
%    p - n-dimensional permutation vector representing P
%    beta - initial search radius to be used in the search process

% Subfunctions: radius, reorder (included in this file)
%               obils_reduction (not included in this file)
% 
% Main References: 
% [1] X.-W. Chang, and X. Yang. An efficient tree search
%     decoder with column reordering for underdetermined MIMO systems. 
%     Global Telecommunications Conference, 2007. 5 pages
% [2] 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, Zhilong Chen
% Copyright (c) 2019. Scientific Computing Lab, McGill University.
% Last revision: June 2019


    [m,n] = size(B);
    
    % Reduce the problem to a full-row-rank one
    [Q,R,p] = qr(B,'vector');
    % Determine the rank of B
    r = sum(abs(diag(R(1:min(m,n),1:min(m,n)))) > n*eps(abs(R(1,1)))); 
    if r == n
        error('The problem is overdetermined, use obils rather than ubils!')
    else
        R = R(1:r,:);
        y = Q' * y; y = y(1:r);
        l = l(p); u = u(p);
        m = r;
    end
    
    % Compute an initial search radius
    [beta,z0] = radius(R,y,l,u);
    
    % If beta is 0, z0 is the solution
    if beta <= n*eps*norm(abs(y)+abs(R)*abs(z0))
        return
    end
    
    if cp == 0
        return
    end
    
    % Permuate z0 to x0
    x0 = zeros(n,1);
    for i = 1 : n
        x0(p(i)) = z0(i);
    end
    
    index = 0;
    prod = inf;
    k = m;
    
    % Determine column m of R
    for j = 1:m
        R1 = R; 
        y1 = y; 
        l1 = l; 
        u1 = u;   
            
        if j ~= m
        % Interchanger columns m and j of R1
           R1(:,[j,m]) = R1(:,[m,j]);
           l1([j,m]) = l1([m,j]);
           u1([j,m]) = u1([m,j]);
        end
        
        % Transform R1 to an upper trapezoidal matrix by plane rotations
        % and apply same rotations to y1
        for i = m:-1:j+1
            [W,~] = planerot(R1(i-1:i,j));
            R1(i-1:i,:) = W * R1(i-1:i,:);
            y1(i-1:i) = W * y1(i-1:i,:);
        end
        
        for i = j+1:m-1
            [W,~] = planerot(R1(i:i+1,i));
            R1(i:i+1,:) = W * R1(i:i+1,:);
            y1(i:i+1) = W * y1(i:i+1);
        end
        
        if y1(m) > 0
            R1(m,:) = -R1(m,:);
            y1(m) = -y1(m);
        end
        
        % Column reordering for (part of) R(:,m:n)    
        [piv,ind_tmp,prod_tmp] = reorder(R1,y1(m),l1,u1,beta);
          
        if ind_tmp > index
            index = ind_tmp;
            k = j;
            R_temp = R1(:,piv);
            y_temp = y1;
            piv_temp = piv;
        elseif ind_tmp == index && prod_tmp < prod
            prod = prod_tmp;
            k = j;
            R_temp = R1(:,piv);
            y_temp = y1;
            piv_temp = piv; 
        end  
    end
    
    R = R_temp;
    y = y_temp;
    l([k,m]) = l([m,k]);  l= l(piv_temp);  
    u([k,m]) = u([m,k]);  u = u(piv_temp);
    p([k,m]) = p([m,k]);  p = p(piv_temp);
    
    % Column reordering for R(:,m:index-1)
    if index > m+1
        [~,id] = sort(abs((R(m,m:index-1))));
        R(:,m:index-1) = R(:,id+m-1);
        p(m:index-1) = p(id+m-1);
        l(m:index-1) = l(id+m-1);
        u(m:index-1) = u(id+m-1);
    end
        
    % Apply the reduction strategy to reorder the first m-1 columns of R
    x = zeros(n,1);
    y_m = y(m);
    for i = n:-1:m
        if y_m == 0 && R(m,i) == 0 % Avoid NaN
                xr = inf;
        else
                xr = round(y_m/R(m,i));
        end
        x(i) = median([l(i),u(i),xr]);
        y_m = y_m - R(m,i)*x(i);
    end
        
    y1 = y(1:m-1) - R(1:m-1,m:n) * x(m:n);
    R1 = R(1:m-1,1:m-1);
    [~,~,l(1:m-1),u(1:m-1),piv]= obils_reduction(R1,y1,l(1:m-1),u(1:m-1));   
    R(:,1:m-1) = R(:,piv);
    [Q,R] = qr(R);
    y = Q'*y;
    p(1:m-1) = p(piv);
    
    % Permute x0 to z0 using the final permutation vector p
    z0 = x0(p);
    
end


%========================================================================
%================= Subfunction of ubils_reduction =======================
%========================================================================

function [beta,x0] = radius(B,y,l,u)
%
% Find an initial search radius determined  by the rounded value of 
% a real solution of the box-constrauined least squares problem. 
%
% Inputs:
%    B - m by n real matrix
%    y - m-dimensional real vector
%    l - n-dimensional integer vector, lower bound 
%    u - n-dimensional integer vector, upper bound    
%
% Outputs:
%    beta - residual at x0 to be used as a asearch radius
%    x0 - the integer point determining beta

% Subfunction:  gradproj

% Compute two initial points
    [Q,R] = qr(B',0);
    x1 = median([l,round(Q*(R'\y)),u],2);
    x2 = round((l+u)/2);
    
    % Compute two box-constrained real LS solutions by the gradient projection
    % method
    x1_real = gradproj(x1,B,y,l,u,100);
    x2_real = gradproj(x2,B,y,l,u,100);
    
    % Round the box-constrained real solution to an integer vector
    x1 = round(x1_real);
    x2 = round(x2_real);
    
    % Initial 
    beta1 = norm(y-B*x1);
    beta2 = norm(y-B*x2);
    
    if beta1 <= beta2
        x0 = x1;
        beta = beta1;
    else
        x0 = x2;
        beta = beta2;
    end
end

%========================================================================
%================= Subfunction of ubils_reduction =======================
%========================================================================

function x = gradproj(x,B,y,l,u,max_iter)
%
% Find a solution to the box-constrained real least squares problem 
% min_{l<=x<=u}||y-Bx|| by the gradient projection method
%
% Inputs:
%    x - n-dimensional real vector as an initial point
%    B - m by n real matrix
%    y - m-dimensional real vector
%    l - n-dimensional integer vector, lower bound 
%    u - n-dimensional integer vector, upper bound    
%
% Output:
%    x - n-dimensional real vector, a solution

    n = length(x);
    
    c = B'*y;
    
    for iter = 1:max_iter
        
        g = B'*(B*x-y);
        
        % Check KKT conditions
        if (x==l) == 0 
            k1 = 1;
        elseif (g(x==l) > -1.e-5) == 1
            k1 = 1;
        else 
            k1 = 0;
        end
        if (x==u) == 0 
            k2 = 1;
        elseif (g(x==u) < 1.e-5) == 1
            k2 = 1;
        else 
            k2 = 0;
        end
        if (l<x & x<u) == 0 
            k3 = 1;
        elseif (g(l<x & x<u) < 1.e-5) == 1
            k3 = 1;
        else 
            k3 = 0;
        end
        if (k1 && k2 && k3)
            break
        end
        
        % Find the Cauchy point
        t_bar = 1.e5*ones(n,1);
        t_bar(g<0) = (x(g<0)-u(g<0))./g(g<0);
        t_bar(g>0) = (x(g>0)-l(g>0))./g(g>0);
        
        % Generate the ordered and non-repeated sequence of t_bar
        t_seq = unique([0;t_bar]);   % Add 0 to make the implementation easier 
        
        % Search
        for j = 2:length(t_seq)
            tj_1 = t_seq(j-1); 
            tj = t_seq(j); 
            % Compute x(t_{j-1})
            xt_j_1 = x - min(tj_1,t_bar).*g;
            % Compute teh search direction p_{j-1}
            pj_1 = zeros(n,1);
            pj_1(tj_1<t_bar) = -g(tj_1<t_bar);
            % Compute coefficients
            q = B*pj_1;
            fj_1d = (B*xt_j_1)'*q - c'*pj_1;
            fj_1dd = q'*q;
            t = tj;     
            % Find a local minimizer
            delta_t = -fj_1d/fj_1dd;
            if fj_1d >= 0
                t = tj_1;
                break;
            elseif delta_t < (tj-tj_1)
                t = tj_1+delta_t;
                break;
            end
        end
        
        x = x - min(t,t_bar).*g;
        
    end
end

%========================================================================
%================= Subfunction of ubils_reduction =======================
%========================================================================

function [piv,index,prod] = reorder(R,y_m,l,u,beta)
%
% Find a new ordering for the last part of columns of R 
%
% Inputs:
%    R - m by n real matrix
%    y_m - a real scalar, the last elemement of m-dimensional vector y
%    l - n-dimensional integer vector, lower bound 
%    u - n-dimensional integer vector, upper bound  
%    beta - search radius
%
% Output:
%    piv - n-dimensional x - n-dimensional real vector, a solution

    [m,n] = size(R);
    r = R(m,:);
    index = m-1;
    prod = 1;
    j = n;
    piv = 1:n;
    
    r(m:n) = abs(r(m:n));
    alpha = y_m + dot(u(m:n)-l(m:n),r(m:n));
    x = zeros(n,1);
    s1 = 0;
    s2 = dot(r(m:n-1),u(m:n-1)-l(m:n-1));
    
    while j >= m
        L_j = inf;
        
        for i = m:j-1
            mu = (alpha+beta-s1) / r(i);
            lambda = (alpha-beta-s1-s2) / r(i);
            L_i = min(u(i)-l(i),floor(mu)+sign(mu-floor(mu))-1) ...
                - max(0,ceil(lambda)-sign(ceil(lambda)-lambda)+1) + 1;
            if L_i < L_j
                L_j = L_i;
                p = i;
                minlambd = lambda;
                minmu = mu;
            end
        end
        
        if p ~= j
            R(:,[p,j]) = R(:,[j,p]);
            piv([p,j]) = piv([j,p]);
            r([p,j]) = r([j,p]);
            l([p,j]) = l([j,p]);
            u([p,j]) = u([j,p]);
        end
        
        if L_j <= 0
            index = j;
            break;
        end
          
        s2 = s2 - r(j)*(u(j)-l(j));
        c = round((alpha-s1)/r(j));
        left = max(0,ceil(minlambd));
        right = min(u(j)-l(j),floor(minmu));
        x(j) = max(min(c,right),left);   
        s1 = s1 + r(j)*x(j);
        prod = prod * L_j;
        j = j - 1;
    end

end

%========================================================================
%================= Subfunction of ubils_reduction =======================
%========================================================================

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-2018. Scientific Computing Lab, McGill University.
% Last revision: December 2018

    [~,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 ubils ============================
%========================================================================

function [zhat,flag]= ubils_search(R,y,l,u,z0,beta,duetime)
%
% [zhat flag] = ubils_search(R,y,l,u,z0,beta,duetime) produces the optimal
% solution or a suboptimal 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 - m by n real upper triapezoidal matrix
%    y - m-dimensional real vector
%    l - n-dimensional integer vector, lower bound
%    u - n-dimensional integer vector, upper bound
%    z0- n-dimensional initial integer point found in the reduction process
%    beta - initial search radius
%    duetime - time limit (a time instance) for the search process
%
% Outputs:
%    zhat - n-dimensional integer vector (in double precision).
%    flag - flag = 0, zhat is the optiaml solution
%           flag = 1, zhat is a suboptimal solution

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

% Main References:
% [1] X.-W. Chang and X. Yang. An efficient tree search
%     decoder with column reordering for underdetermined MIMO systems.
%     Proceedings of IEEE GLOBECOM 2007, 5 pages.
% [2] Jing Zhu. Numerical methods for underdetermined box-constrained
%     integer least squares problems, M.Sc. Thesis, School of Computer
%     Science, McGill University, 2016.

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


    % ------------------------------------------------------------------
    % --------  Initialization  ----------------------------------------
    % ------------------------------------------------------------------
    [m,n] = size(R);
    
    % Current point
    z = [];
    zhat = [];
    
    % c(k) = (y(k)-R(k,k+1:n)*z(k+1:n))/R(k,k)
    c = 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);
    
    % 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);
    
    % 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;
    
    
    % ------------------------------------------------------------------
    % --------  Search process  ----------------------------------------
    % ------------------------------------------------------------------
    
    % Shift variables
    [R,y,u,l,col_neg,s] = shift(R,y,l,u);
    
    utemp = u;
    
    % S(k,n) = y(k), k=1:m-1, S(k,n) = y(m), k=m:n
    S(m:n,n) = y(m);
    S(1:m-1,n) = y(1:m-1);
    
    % Choose search order
    
    % Initialization
    interval_length = u(n) - l(n);
    all_dis = zeros(interval_length,1);
    results = cell(interval_length+1,1);
    visit_order = zeros(interval_length,1);
    
    % Preprocessing 
    for i = 1:interval_length+1
        c(n,1) = i - 1;
        [dis,branch_info] = initialSearch(S,R,c,l,u,beta,y);
        if size(branch_info,1) == 0
            all_dis(i,1) = Inf;
        else
            all_dis(i,1) = dis;
            results{i,1} = branch_info;
        end
    end
    % Find visiting order
    cur_length = 1;
    
    % Check if the cell aray {results} is empty
    if isempty(results) == 1
        zhat = z0;
        flag = 0;
        return
    end
    
    % Sort the order
    while cur_length <= interval_length+1
        cur_min = min(all_dis);
        if(cur_min ~= Inf)
            indexes = find(all_dis == cur_min);
            for j = 1:1:size(indexes,1)
                cur_index = indexes(j,1);
                visit_order(cur_length) = cur_index;
                cur_length = cur_length + 1;
                all_dis(cur_index,1) = Inf;
            end
        else
            break;
        end
    end
          
    % Set up the first branch to be searched
    cur_index = 1;
    cur = results{visit_order(cur_index),1};
    
    l = cur{1,1};u = cur{1,2};
    lflag = cur{1,3};      % lflag(k) = 1 if lower bound is reached at level k
    uflag = cur{1,4};      % uflag(k) = 1 if upper bound is reached at level k
    ulevel = cur{1,5};     % The level where search moves up to a higher level
    path = cur{1,6};       % Record information for updating S(k,k:path(k)-1)
    S = cur{1,7};
    c = cur{1,8};
    z = cur{1,9};          % Current point
    d = cur{1,10};         % d(k): left or right search direction at level k
    prsd = cur{1,11};      % Partial squared residual norm for z
    k = max(m - 1,1);
    dflag = 1;             % Intend to move down to a lower level
    
    
    while 1 && cputime < duetime
        if k <= m && prsd(k+1) > beta^2
            dflag = 0;
        end
        if dflag == 1
            if lflag(k) ~= 1 || uflag(k) ~= 1
                if k ~= 1 % Move to level k-1
                    k1 = 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
                    if k <= m
                        for j = path(k1):-1:k1 + 1
                            S(k1,j-1) = S(k1,j) - R(k1,j) * z(j);
                        end
                        c(k1) = S(k1,k1) / R(k1,k1); 
                    else
                        S(m:k1,k1) = S(k1,k1+1) - R(m,k1+1) * z(k1+1);
                        if R(m,k1)~=0
                        c(k1) = S(k1,k1) / R(m,k1);
                        else 
                           c(k1) = l(k1);
                        end
                    end
                    
                    % Compute the new bound at level k1
                    [l(k1),u(k1)] = bound(c(k1),R,utemp,beta,prsd,k1);
                    % Find the initial integer in [l(k1), u(k1)]
                    if l(k1) > u(k1)
                        lflag(k1) = 1;
                        uflag(k1) = 1;
                    else
                        [z(k1),d(k1),lflag(k1),uflag(k1)] = init(c(k1),l(k1),u(k1));
                        if k1 <= m
                            gamma = R(k1,k1)*(c(k1) - z(k1));
                            prsd(k1)= prsd(k1+1) + gamma * gamma;
                        end
                    end
                    k = k1;
                    ulevel = 0;
                else % A valid point is found
                    zhat = z;
                    beta = sqrt(prsd(1));
                    for j = 1:n
                        [l(j),u(j)] = bound(c(j),R,utemp,beta,prsd,j);
                    end
                    dflag = 0;
                end
            else % Will move back to a higher level
                dflag = 0;
            end
        else
            if k == n - 1 % The optimal solution has been found, terminate
                cur_index = cur_index + 1;
                if cur_index == cur_length
                    break;
                end
                
                % else Switch to the next ordered branch
                cur = results{visit_order(cur_index),1};
                l = cur{1,1};
                u = cur{1,2};
                lflag = cur{1,3};
                uflag = cur{1,4};
                ulevel = cur{1,5};
                path = cur{1,6};
                S = cur{1,7};
                c = cur{1,8};
                z = cur{1,9};
                d = cur{1,10};
                prsd = cur{1,11};
                k = max(m - 1,1); % In case m = 1 
                dflag = 1;
                if k == 1
                    dflag = 0;
                end
            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));
                    if k <= m
                        gamma = R(k,k)*(c(k)-z(k));
                        % Update partial squared residual norm
                        prsd(k)= prsd(k+1) + gamma * gamma;
                    end
                    dflag = 1;
                end
            end
        end
    end
    
    if cputime >= duetime % Exceed due time
        flag = 1;
    else
        flag = 0;
    end
    
    if isempty(zhat) % No integer point found within the search radius
        zhat = z0;
    else  % Shift back, see subfunction shift
        zhat(col_neg) = -zhat(col_neg);  
        zhat = zhat + repmat(s,1,1);  
    end
    
end


%========================================================================
%==================== Subfunction of ubils_search =======================
%========================================================================

function [R,y,u,l,col_neg,s] = shift(R,y,l,u)

    % Perform transformations: 
    %    z(j) := -l(j) + z(j) if r_{tj} >= 0  for t=m or t=j
    %    z(j) := u(j) - z(j)  if r_{tj} < 0   for t=m or t=j
    % So that z(j) is in {0,1, . . . ,u_j-l_j} and R(t,j)>=0, j=1:n
    
    [m,n] = size(R);
    R_temp = R;
    if m == 1
        col_neg = R < 0;
    else
        col_neg = [diag(R)', R(m,m+1:n)] < 0;
    end
    R(:,col_neg) = - R(:,col_neg);
    s = l .* (1 - col_neg)' + u .* col_neg';
    
    y = y - R_temp * s;
    u = u - l;
    l = zeros(size(l));

end

%========================================================================
%==================== Subfunction of ubils_search =======================
%========================================================================

function [dis,branch_info] = initialSearch(S,R,c,l,u,beta,y)
    % Search z(m:n) and get relevant information to be used in determining
    % the branch order for the search process
    
    [m,n] = size(R);
    
    % Initialization
    utemp = u;
    dflag = 1; % dflag for down or up search direction
    k = n;
    dis = Inf;
    branch_info = [];
    lflag = zeros(n,1); % lflag(k) = 1 if lower bound is reached at level k
    uflag = zeros(n,1); % uflag(k) = 1 if upper bound is reached at level k
    d = zeros(n,1);     % Left or right search direction at level k
    path = n*ones(n,1); % path(k): record information for updating S(k,k:path(k)-1)
    ulevel = 0;         % The level where search moves up to a higher level
    prsd = zeros(m+1,1); % Partial squared residual norm for z
    z = zeros(n,1);
    [l(k),u(k)] = bound(c(k),R,u,beta,prsd,k);
    [z(k),d(k),lflag(k),uflag(k)] = init(c(k),l(k),u(k));
    while 1
        if k <= m && prsd(k+1) > beta^2
            dflag = 0;
        end
        if dflag == 1
            if lflag(k) ~= 1 || uflag(k) ~= 1
                if k ~= m-1 && k ~=1 % Move to level k-1 
                    k1 = 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
                    if k <= m
                        for j = path(k1):-1:k1+1
                            S(k1,j-1) = S(k1,j) - R(k1,j) * z(j);
                        end
                        c(k1) = S(k1,k1) / R(k1,k1);
    
                    else
                        S(m:k1,k1) = S(k1,k1+1) - R(m,k1+1) * z(k1+1);
                        if R(m,k1) ~= 0
                            c(k1) = S(k1,k1) / R(m,k1);
                        else
                            c(k1) = l(k1);
                        end
                    end
                    
                    % Compute new bound at level k1
                    [l(k1),u(k1)] = bound(c(k1),R,utemp,beta,prsd,k1);
                    
                    % Find the initial integer in [l(k1), u(k1)]
                    if l(k1) > u(k1)
                        lflag(k1) = 1;
                        uflag(k1) = 1;
                    else
                        [z(k1),d(k1),lflag(k1),uflag(k1)] = init(c(k1),l(k1),u(k1));
                        if k1 <= m
                            gamma = R(k1,k1) * (c(k1) - z(k1));
                            prsd(k1)= prsd(k1+1) + gamma * gamma;
                        end
                    end
                    k = k1;
                    ulevel = 0;
                else % The first valid x_{m:n} at this branch has been found
                    dis_y = (y(m)-R(m,m:n)*z(m:n))^2; 
                    y_new = y(1:m-1)-R(1:m-1,m:n)*z(m:n);
                    lowerbound = 0;
                    if lowerbound+dis_y > beta^2
                        dflag = 0;
                    else
                        branch_info = {l,u,lflag,uflag,ulevel,path,S,c,z,d,prsd};
                        x_rls = R(1:m-1,1:m-1)\y_new;
                        x_bils = median([l(1:m-1),round(x_rls),u(1:m-1)],2);
                        dis = norm(y_new-R(1:m-1,1:m-1)*x_bils,2)^2 + dis_y;
    
                        break;
                    end
                end
            else % Will move back to a higher level
                dflag = 0;
            end
        else
            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));
                    if k <= m
                        gamma = R(k,k) * (c(k) - z(k));
                        prsd(k)= prsd(k+1) + gamma * gamma;
                    end
                    dflag = 1;
                end
            end
        end
    end
end

%========================================================================
%==================== Subfunction of ubils_search =======================
%========================================================================

function [l_k,u_k] = bound(c_k,R,u,beta,prsd,k)
    % Compute new lower bound and upper bound for z_k
    
    m = size(R,1);
    if k > m - 1
        lambda_k = c_k - (beta+R(m,m:k-1)*u(m:k-1))/R(m,k);
        mu_k = c_k + beta/R(m,k);
    else
        lambda_k = c_k - sqrt(beta^2-prsd(k+1))/R(k,k);
        mu_k = c_k + sqrt(beta^2-prsd(k+1))/R(k,k);
    end
    
    if lambda_k - floor(lambda_k) < 1e-12 && lambda_k ~= 0
        lambda_k = floor(lambda_k);
    end
    
    if ceil(mu_k) - mu_k < 1e-12 && mu_k ~= u(k)
        mu_k = ceil(mu_k);
    end
    
    if lambda_k == Inf || lambda_k == -Inf
        l_k = 0;
        u_k = u(k);
    else
        l_k = ceil(max(0,lambda_k));
        u_k = floor(min(u(k),mu_k));
    end
end

%========================================================================
%==================== Subfunction of ubils_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 ubils_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.

    if lflag_k == 0 && uflag_k == 0
        zk = z_k + d_k;
        if zk > u_k
            uflag_k = 1;
        elseif zk < l_k
            lflag_k = 1;
        else
            z_k = zk;
            if d_k > 0
                d_k = -d_k - 1;
            else
                d_k = -d_k + 1;
            end
        end
    end
    
    if lflag_k == 1 && uflag_k == 0
        z_k = z_k + 1;
        if z_k > u_k
            uflag_k = 1;
        end
    elseif lflag_k == 0 && uflag_k == 1
        z_k = z_k - 1;
        if z_k < l_k
            lflag_k = 1;
        end
    end

end

