function [Phi_best,w_best]=HugoFinal(F,A,b,w0,w1,t_max)
% Input:
% F ... the mxn matrix of regressors of LRM, where
%   m is the dimension of the parameter, n is the size of the design space.
% A,b ... the kxn matrix, and the kx1 vector of constraints Aw<=b.
%   A must have non-negative elements, b must have positive elements.
%   Each column of A must have at least one strictly positive element.
% w0 ... the nx1 vector representing the design to be augmented.
%   Design w0 must satisfy A*w0<=b.
% w1 ... the nx1 vector representing the initial design.
%   If w1 does not satisfy (w1>=w0 and A*w1<=b),
%   the algorithm uses a random maximal design as initial.
% t_max ... upper limit on the time taken by the computation.
% 
% Output: 
% w_best is the best design found and Phi_best is its criterial value.

back_max=16;   % Maximal number of backward steps in one excursion.
n_round=9;     % The number of significant digits used for attributes.

eps=1/1000000; % Used to correct small numerical errors.
V_max=500000;  % Maximal size of the tabu vector.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Criterion of D-optimality with a regularization term.
% Change to the appropriate monotonic criterion to be maximized.
% The criterion should return non-negative values.

    function value=Phi(w) % D-optimality
        value=(det(F*diag(w)*F'+eps*eye(m)))^(1/m);
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% - Vector w represents the current design. The design has its
%   characteristics res, d and attr that are always updated with w.
% - Vector res=b-A*w determines the residuals amount of the k resources.
% - Vector d represents numbers of possible additional trials in pts of w.
% - attr is the "characteristic attribute" of the design w to build V.
% - val is the "local heuristic evaluation" of the design w.


    function d=compute_d(del_act)
    % Calculates d for a design with del=del_act.
        d=(floor(min(kron(del_act,ones(1,n))./A,[],1)+eps))';
    end

    function xr=round_sig(x,k)
    % Rounds x to k significant digits; x must be a non-negative scalar.
       if(x>0)
         l = 10^(floor(log10(x+eps)-k+1));
         xr = round(x/l)*l;
       else
         xr=0;
       end
    end

    function gamma=compute_gamma(numer,denom)
    % Computes the minimum of the vector numer./denom from all components where
    % denom is not zero. If all denoms are zero, the result is zero.
       if ~any(denom>eps)
           gamma=0;
       else
           ratio=numer./denom;
           gamma=min(ratio(denom>eps));
       end
    end

    function up_index=explore_up
    % Returns the index of the "best" nonV upper neighbor of w
    % If there is no nonV upper neigbor, returns 0.
        wei=w; val=-1*ones(1,n);
        for i=1:n
           if d(i)>0
               wei(i)=w(i)+1;
               attrei=round_sig(Phi(wei),n_round);
               if ~ismember(attrei,V)
                   dei=compute_d(res-A(:,i));
                   gammaei=compute_gamma(res-A(:,i),A*dei);
                   val(i)=Phi(wei+gammaei*dei);
               end
               wei(i)=w(i);
           end
        end;
        [max_val,up_index]=max(val);
        if max_val<(-eps)
            up_index=0;
        end
     end

    function down_index=explore_down
    % Returns index of a nonV downbor of w with the greatest val.
    % If there is no nonV downbor, return 0.
        wei=w; val=-1*ones(1,n);
        for i=1:n
           if w(i)>w0(i)
               wei(i)=w(i)-1;
               attrei=round_sig(Phi(wei),n_round);
               if ~ismember(attrei,V)
                   dei=compute_d(res+A(:,i));
                   gammaei=compute_gamma(res+A(:,i),A*dei);
                   val(i)=Phi(wei+gammaei*dei);
               end
               wei(i)=w(i);
           end
        end;
        [max_val,down_index]=max(val);
        if max_val<(-eps)
            down_index=0;
        end
    end

    function random_step
    % Changes w to a random design in L cup U.
    % Efficiency of this rejection method is not the best possible,
    %   but a potential speedup is probably insignificant.

        success=false;
        while(~success);
            dir=randsample([0 1],1);
            if(dir)
                up_index=randsample(1:n,1);
                if d(up_index)>0
                    w(up_index)=w(up_index)+1;
                    res=res-A(:,up_index); d=compute_d(res);
                    attr=round_sig(Phi(w),n_round);
                    success=true;
                end
            else
                down_index=randsample(1:n,1);
                if w(down_index)>w0(down_index)
                    w(down_index)=w(down_index)-1;
                    res=res+A(:,down_index); d=compute_d(res);
                    attr=round_sig(Phi(w),n_round);
                    success=true; back_no=back_no+1;
                end
            end
        end
    end

    function random_initialization
    % Returns a random maximal design w.
        w=w0; res=b-A*w; d=compute_d(res);
        maximal=false;
        while ~maximal
            val=(d>0);
            if sum(val)==0
                maximal=true;
            else
                i=randsample(1:n,1,true,val);
                w(i)=w(i)+1;
                res=res-A(:,i);
                d=compute_d(res);
            end
        end
        attr=round_sig(Phi(w),n_round);
    end
 
%%%%%%%%%%%%%%%%%%%%%
%  MAIN BODY OF HUGO

tic; [m,n]=size(F);
if (length(w1)~=n)||(min(b-A*w1)<0)||(min(w1-w0)<0)
    random_initialization;
else
    w=w1; res=b-A*w; d=compute_d(res);
    attr=round_sig(Phi(w),n_round);
end

w_best=w; res_best=res; d_best=d;
attr_best=attr; Phi_best=Phi(w);
V=[]; back_no=0; finish=false;

while(~finish)
    if(~ismember(attr,V))
        V=union(attr,V);
        up_index=explore_up;
        if(up_index)           
            w(up_index)=w(up_index)+1;
            res=res-A(:,up_index); d=compute_d(res);
            attr=round_sig(Phi(w),n_round);
        else
            if ~any(d)
                Phi_w_cmp=Phi(w);
                if Phi_w_cmp>Phi_best
                    w_best=w; res_best=res; d_best=d;
                    attr_best=attr; Phi_best=Phi_w_cmp;
                    back_no=0;
                end
            end
            down_index=explore_down;
            if(down_index)       
                w(down_index)=w(down_index)-1;
                res=res+A(:,down_index); d=compute_d(res);
                attr=round_sig(Phi(w),n_round);
                back_no=back_no+1;
            else
                random_step;
            end   
        end
    else
        down_index=explore_down;
        if(down_index)
            w(down_index)=w(down_index)-1;
            res=res+A(:,down_index); d=compute_d(res);
            attr=round_sig(Phi(w),n_round);
            back_no=back_no+1;
        else
            up_index=explore_up;
            if(up_index)
                w(up_index)=w(up_index)+1;
                res=res-A(:,up_index); d=compute_d(res);
                attr=round_sig(Phi(w),n_round);
            else
                random_step;
            end
        end
    end
   
    if back_no>back_max
        w=w_best; res=res_best; d=d_best;
        attr=attr_best; back_no=0;
    end
       
    if (toc>t_max)||(length(V)>=V_max)
        finish=true;
    end
end

Phi_best=(det(F*diag(w_best)*F'))^(1/m);

end