function [P,fup]=normal_rectangular_branch_bound(Asqr,b,Acol,Arow,U,V,X,Y,depth0,err0,fig)
assert(size(U,2)==size(V,1));

Nx=size(X,1);
Ny=size(Y,1);

options=optimset('Display','off');
%%%bounding rectangle
fup=inf;
rect=zeros(size(U,2),2);
for ii=1:size(U,2)
    
    [p,rect(ii,1)]=linprog(U(:,ii),...
        Acol,ones(size(Acol,1),1),...
        Arow,ones(size(Arow,1),1),...
        zeros(size(U,1),1),ones(size(U,1),1),[],options);%min
    %
    fval=b'*p-sumsqr(Asqr*p);%function value
    if fval<fup
        fup=fval;
        popt=p;
    end
    
    [p,rect(ii,2)]=linprog(-U(:,ii),...
        Acol,ones(size(Acol,1),1),...
        Arow,ones(size(Arow,1),1),...
        zeros(size(U,1),1),ones(size(U,1),1),[],options);
    %
    rect(ii,2)=-rect(ii,2);%max
    %
    fval=b'*p-sumsqr(Asqr*p);%function value
    if fval<fup
        fup=fval;
        popt=p;
    end
end


%%%
disp('brach and bound')

depth=0;
setrec={rect,inf,0};%set of rectangles,{rect,flin,depth}
setind=1;%indices of active rectangles
err=err0;
depth0
err0

while depth<=depth0
    for ii=1:length(setind)
        rect=setrec{setind(ii),1};
        [p,flin,fval]=rect_underestimator(Asqr,b,Acol,Arow,U,V,rect);
        setrec{setind(ii),2}=flin;%update flin
        
        if fval<fup
            fup=fval;%update upper bound
            popt=p;%update solution
            %
            if fig>0
                P=reshape(popt(1:size(U,1)),Ny,Nx)';
                plot_intermediate(X,Y,P,fig) ;
                [fup             depth]
            end
        end
    end
    for ii=size(setrec,1):-1:1
        if setrec{ii,2}>=fup-err
            setrec(ii,:)=[];%prune unlikely rectangles
            disp('prune');
        end
    end
    if isempty(setrec)
        break;
    end
    
    flow=inf;
    for ii=1:size(setrec,1)
        if setrec{ii,2}<flow
            flow=setrec{ii,2};%update lower bound
            indlow=ii;%index of active rectangle
        end
    end
    
    rect=setrec{indlow,1};
    depth=setrec{indlow,3};
    setrec(indlow,:)=[];%eliminate the rectangle
    
    %branch
    %decide which dimension to divide
    dif=V.*(rect(:,2)-rect(:,1)).^2;
    ind=find(dif==max(dif),1);
    %new rectangulars
    recti=repmat(rect,[1,1,2]);
    half=mean(rect(ind,:));
    recti(ind,2,1)=half;
    recti(ind,1,2)=half;
    %
    setind=[];
    for ii=1:2
        setrec=[setrec;{recti(:,:,ii),inf,depth+1}];
        setind=[setind;size(setrec,1)];
    end
end

depth

%%%
P=reshape(popt(1:size(U,1)),Ny,Nx)';