function fig9_15

% this include figures 9.15 and 9.16 (SDM/CGM portions)

%  Inputs:
%        x = starting vector
%        xa, xb = x-interval used in contour plot
%        ya, yb = y-interval used in contour plot
%        tol = tolerance for stopping iteration

%  Required (end of file)
%        F    is function F(x,y) to be minimized
%        dF  is gradient of F(x,y)

%  Specify which example to run ifun=1 (Rosen), 2 (multi-min)
ifun=2

% specify which beta choice SDM (ibeta=1), CGM-FR (ibeta=2),
% CGM-PR (ibeta=3)
ibeta=3

clf
hold on
% get(gcf)
set(gcf,'Position', [12 950 421 395])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

n_fun=0;
for ic=1:2
    if ifun==1
        if ic==2
            break
        end
        %x=[-1.0, -0.5];
        %x=[-1.2, 1.8];
        %x=round(2*rand(1,2)-1,2)

        %x=[4/5   -3/5]
        x=[3/4   -3/4]
        %         x=[-1.2 1]
        xa=-1.1; xb=1.1; ya=-1.1; yb=1.1;
        xd=linspace(xa,xb,80);
        yd=linspace(ya,yb,80);
        [X,Y]=meshgrid(xd,yd);
        Z=100*(X.^2-Y).^2+(X-1).^2;
        exact=[1, 1];
        plot(exact(1),exact(2),'r*','MarkerSize',14,'LineWidth',2)
    elseif ifun==2
        if ic==1
            x=[1, 0.8];
            options = optimset('TolX',1e-10,'MaxIter',5000,'MaxFunEvals',5000);
            %ee=fminsearch(@(x)(x(1)-x(2))^4+8*x(1)*x(2)-x(1)+x(2)+3,x,options)
            exact=[0.553579941873744, -0.553579935896837];
        else
            x=[-1, -0.6];
            options = optimset('TolX',1e-10,'MaxIter',5000,'MaxFunEvals',5000);
            %exact=fminsearch(@(x)(x(1)-x(2))^4+8*x(1)*x(2)-x(1)+x(2)+3,[-0.5,0.5],options)
            exact=[-0.418782723467334, 0.4187827111776];
        end
        xa=-1; xb=1; ya=-1.5; yb=1;
        xd=linspace(xa,xb,80);
        yd=linspace(ya,yb,80);
        [X,Y]=meshgrid(xd,yd);
        Z=(X-Y).^4+8*X.*Y-X+Y+3;
    end
    xpoints=x(1);
    ypoints=x(2);
    tol=1e-6;
    %     n_fun=0;
    g=dF(ifun);
    d=-g;
    counter=1;
    F0=F(0);
    %     alphaM=1;
    alphaM=1/2;
    while norm(x-exact,inf)>tol && counter<500
        counter=counter+1;

        %alpha=fminbnd(@F,0,alphaM);
        alpha=arm(alphaM);

        %         alpha=-dot(g,d)/dot(d,H(x)*d')

        %         if alpha>alphaM/3
        %             alphaM=2*alphaM;
        %         elseif alpha<alphaM/64
        %             alphaM=alphaM/2;
        %         end

        x=x+alpha*d;
        g0=g;
        g=dF(ifun);

        if ibeta==3
            % CGM-PR
            beta=dot(g,g-g0)/dot(g0,g0);
            if dot(g,g)-beta*dot(d,g)<0
                beta=0;
            end

        elseif ibeta==2
            % CGM-FR
            beta=dot(g,g)/dot(g0,g0);
        elseif ibeta==1
            beta=0;
        end

        d=-g+beta*d;
        xpoints=[xpoints, x(1)];
        ypoints=[ypoints, x(2)];
    end

    n_fun

    error=norm(exact-x,inf);
    fprintf('\n  Number of Steps = %i',counter)
    fprintf('\n  Computed Solution =  %e  %e\n',x)
    %fprintf('\n  Error = %e \n\n',error)

    box on
    contour(X,Y,Z,40,'b','LineWidth',1)
    axis([xa-0.004 xb ya yb])
    axis square
    plot(xpoints,ypoints,'r-','LineWidth',2)
    plot(xpoints,ypoints,'o','Color',[110/256 0 0],'MarkerSize',10,'LineWidth',2)
end

set(gca,'xtick',[-1 -0.5 0 0.5 1]);
set(gca,'ytick',[-1 -0.5 0 0.5 1]);
xlabel('x-axis')
ylabel('y-axis')
set(gca,'FontSize',18,'FontWeight','bold')
hold off

% if ifun==2 && ibeta==1
%     exportgraphics(gcf,'/Users/mark/Desktop/sdm11.eps')
% elseif  ifun==1 && ibeta==1
%     exportgraphics(gcf,'/Users/mark/Desktop/sdm21.eps')
% elseif ifun==2 && ibeta==3
%     exportgraphics(gcf,'/Users/mark/Desktop/cgm11.eps')
% elseif  ifun==1 && ibeta==3
%     exportgraphics(gcf,'/Users/mark/Desktop/cgm21.eps')
% end

    function s=arm(aM)
        gam=1/20; tau=0.5;
        s=aM;
        Fs=F(s);
        mk=dot(g,d);
        gd=gam*mk;

        for ir=1:20
            for is=1:20
                if Fs<F0 && is>1
                    alphaM=20*s;
                    break
                end
                s=tau^2*s;
                Fs=F(s);
            end

            if Fs<F0+s*gd
                F0=Fs;
                return
            end

            if ir==1
                S=s;
                FS=Fs;
                s=tau*s;
                Fs=F(s);
                if Fs<F0+s*gd
                    F0=Fs;
                    return
                end
            else
                s=tau*s;
                Fs=F(s);
                if Fs<F0+s*gd
                    F0=Fs;
                    return
                end
            end
        end
        z=0
    end

    function z=F(s)
        n_fun=n_fun+1;
        xs=x+s*d;
        if ifun==1
            z=100*(xs(1)^2-xs(2))^2+(xs(1)-1)^2;
        elseif ifun==2
            z=(xs(1)-xs(2))^4+8*xs(1)*xs(2)-xs(1)+xs(2)+3;
        elseif ifun==3
            z=10*xs(1)^2+xs(2)^2;
        end
    end


    function g=dF(iff)
        n_fun=n_fun+1;
        if ifun==1
            g(1)=400*x(1)*(x(1)^2-x(2))+2*(x(1)-1);
            g(2)=-200*(x(1)^2-x(2));
        elseif ifun==2
            g(1)=4*(x(1)-x(2))^3+8*x(2)-1;
            g(2)=-4*(x(1)-x(2))^3+8*x(1)+1;
        elseif ifun==3
            g(1)=20*x(1);
            g(2)=2*x(2);
        end
    end

% jacobian
    function g=H(x)
        n_fun=n_fun+1;
        xx=x(1); yy=x(2);
        g=[1200*xx^2-400*yy+2,-400*xx;-400*xx,200];
    end

end

















