branchresect.m

数学优化工具箱

function [Hlastmat, iternum] = branchresect(xm, Xlist, m, n, optim_bound, beta, alphaval, maxiter, opts);

%maxiter = 50;
nbrpoints = size(Xlist,2);
numvar = m*n+5*nbrpoints;
betainit = beta;
maxdepth = opts.maxdepth;
feas_status = 0;

if opts.normalize == 1
    Timage = hnormalize(xm);
    Tspace = hnormalize(Xlist);
    xmold = xm;
    Xlistold = Xlist;
    xm = Timage*xm;
    Xlist = Tspace*Xlist;
else
    Timage = eye(m);
    Tspace = eye(n);
end

%keyboard;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Choosing the initial bounds based on reprojection error
%% opts.autoinit = 1 adjusts value of beta to best possible
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

numerr2flag = 0;
%Hinit = zeros(nbrpoints,2);
Hinit = [zeros(nbrpoints,1), ones(nbrpoints,1)*Inf];
if opts.autoinit == 1
    beta = 2*beta;
    feas_status = 1;
    numinitsteps = -1;
    feascnt = 0;
    while feas_status > 0 %& feascnt < 6 %& min(Hinit(2:nbrpoints,2)-Hinit(2:nbrpoints,1)) > 1e-4
        feascnt = feascnt + 1;
        Hinitprev = Hinit;
        numinitsteps = numinitsteps + 1;
        fprintf('Initial bounding : pass %d .... \n', numinitsteps);
        beta = beta/2;
        [Hinit, feas_status] = init_bound_resect(xm, Xlist, nbrpoints, beta, m, n, maxdepth, 1);
        if feas_status == 0 & numinitsteps == 0
            numerr2flag = 1;
        end
    end

    % keyboard;
    feas2 = 1;
    if numerr2flag == 1
        numerr2cnt = 1;
        feas2 = 0;
        while numerr2cnt < 4 & feas2 == 0
            beta = beta/2;
            [Hinit, feas2] = init_bound_resect(xm, Xlist, nbrpoints, beta, m, n, maxdepth, 1);
            numerr2cnt = numerr2cnt + 1;
        end
    else
        beta = beta * 2;
        beta = min(beta*nbrpoints, betainit);
    end

    disp('Computing initial rectangle ....');
    [Hinit, feas_status] = init_bound_resect(xm, Xlist, nbrpoints, beta, m, n, maxdepth, 0);
%    keyboard;
    Hinitold = Hinit;
    if feas_status == 1 & min(Hinit(2:nbrpoints,2)-Hinit(2:nbrpoints,1)) < 0.01
        disp('Widening the initial rectangle.');
        for i = 2:nbrpoints
            Hinit(i,1) = max(0,Hinit(i,1)-0.005);
            Hinit(i,2) = Hinit(i,2)+0.005;
        end
    end


    if numinitsteps < 1
        fprintf('Initial reprojection error bound too small or large (Sedumi bug).\n\n');
    end

else
    disp('Computing initial rectangle ....');
    [Hinit, feas_status] = init_bound_resect(xm, Xlist, nbrpoints, beta, m, n, maxdepth, 0);
    Hinitold = Hinit;
    if feas_status == 1 & min(Hinit(2:nbrpoints,2)-Hinit(2:nbrpoints,1)) < 0.01
        disp('Widening the initial rectangle.');
        for i = 2:nbrpoints
            Hinit(i,1) = max(0,Hinit(i,1)-0.005);
            Hinit(i,2) = Hinit(i,2)+0.005;
        end
    end

    %keyboard;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%keyboard;
%% The initial reconstruction
rectlist(1).tllist = 0.001*ones(1,nbrpoints);
rectlist(1).tulist = 1000*ones(1,nbrpoints);
rectlist(1).sllist = Hinit(:,1)';
rectlist(1).sulist = Hinit(:,2)';

if feas_status == 0
    disp('Switching to default initial region.');
    rectlist(1).tllist = 0.001*ones(1,nbrpoints);
    rectlist(1).tulist = 1000*ones(1,nbrpoints);
    rectlist(1).sllist = 0.001*ones(1,nbrpoints);
    rectlist(1).sulist = 1000*ones(1,nbrpoints);
end


% The initial resection
%keyboard;
[hlist, lowerbound, objval] = boundresect(nbrpoints, numvar, rectlist(1), xm, Xlist, m, n);
if feas_status == 1 & lowerbound == Inf & objval == Inf
    disp('Hmmmm .... shrinking back the initial rectangle.');
    rectlist(1).tllist = 0.001*ones(1,nbrpoints);
    rectlist(1).tulist = 1000*ones(1,nbrpoints);
    rectlist(1).sllist = Hinitold(:,1)';
    rectlist(1).sulist = Hinitold(:,2)';
    [hlist, lowerbound, objval] = boundresect(brpoints, numvar, rectlist(1), xm, Xlist, m, n);
end

envminlist(1) = lowerbound;
objvallist(1) = objval;
Hlast = hlist.y(1:m*n);
vopt = objval;
rectcnt = 1;
iternum = 0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Branch and Bound
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while rectcnt > 0 & min(envminlist ./ vopt) < optim_bound & iternum <= maxiter

    iternum = iternum + 1;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Rectangle with lowest envelope is chosen to split
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [ytmp, index] = min(envminlist);
    inrect = rectlist(index);

    sLlist = inrect.sllist;
    sUlist = inrect.sulist;
    sdifflist = sUlist - sLlist;

    %% Rectangle "index" is split along longest dimension "index2"

    if opts.onlyn == 1
        [ytmp2, index2] = max(sdifflist(2:n+1));
    else
        [ytmp2, index2] = max(sdifflist(2:end));
    end
    index2 = index2+1;

    sL = sLlist(index2);
    sU = sUlist(index2);
    border = sL + (sU - sL) * alphaval;

    inrect1old = inrect;  inrect2old = inrect;

    %% Reset borders
    inrect1old.sllist(index2) = sL;
    inrect1old.sulist(index2) = border;
    inrect2old.sllist(index2) = border;
    inrect2old.sulist(index2) = sU;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Bound propagation turned on by opts.prop = 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if opts.prop == 1 & nbrpoints-1 > n
        %tic;
        inrect1 = propagateresect(inrect1old, index2, Xlist, nbrpoints, m, n);
        inrect2 = propagateresect(inrect2old, index2, Xlist, nbrpoints, m, n);
        %fprintf('Propagate time: %f\n', toc);
    else
        inrect1 = inrect1old;
        inrect2 = inrect2old;
    end

    rectlist = [rectlist(1:index-1), inrect1, inrect2, rectlist(index+1:end)];
    rectcnt = rectcnt + 1;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Bounding
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %tic;
    [h1new, lowerbound1, objval1] = boundresect(nbrpoints, numvar, inrect1, xm, Xlist, m, n);
    %fprintf('Bounding time : %f\n', toc);
    if vopt > objval1
        vopt = objval1;
        Hlast = h1new.y(1:m*n);
    end

    [h2new, lowerbound2, objval2] = boundresect(nbrpoints, numvar, inrect2, xm, Xlist, m, n);
    if vopt > objval2
        vopt = objval2;
        Hlast = h2new.y(1:m*n);
    end

    %% Update lists
    envminlist = [envminlist(1:index-1), lowerbound1, lowerbound2, envminlist(index+1:end)];
    objvallist = [objvallist(1:index-1), objval1, objval2, objvallist(index+1:end)];
    hlist = [hlist(1:index-1), h1new, h2new, hlist(index+1:end)];

    %% Remove rectangles with envelope minimum above current vopt
    indretain = find(envminlist < vopt);
    rectcnt = length(indretain);
    if rectcnt == 0
        disp('No rectangles retained.');
        %vpa(envminlist,20)
        %vpa(vopt,20)
    end

    %% Update lists
    rectlist = rectlist(indretain);
    envminlist = envminlist(indretain);
    objvallist = objvallist(indretain);
    hlist = hlist(indretain);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Output display
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fprintf('Iter: %d  Rectangles: %d  vopt: %f  ', iternum, rectcnt, vopt);
    fprintf('Optimality status: %f  Required: %f\n', min(envminlist ./ vopt), optim_bound);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end;

Hlastmat = reshape(Hlast,n,m)';
Hlastmat = inv(Timage)*Hlastmat*Tspace;
Hlastvec = reshape(Hlastmat',[size(Hlastmat,1)*size(Hlastmat,2),1]);
Hlastvec = Hlastvec/Hlastvec(end);

fprintf('Number of iterations : %d\n', iternum);
%fprintf('Optimal value of y : ');
%fprintf('%f ', Hlastvec);
%fprintf('\n');

%residual = sum((Hlastvec-ytrue).^2);
%fprintf('Final residual : %2.12f', sum((Hlastvec-ytrue).^2));
%fprintf('\n\n');