function optimshow

Fig = openfig('optimshow.fig');
H = guihandles(Fig);

ProblemList = { ...
    'Unknown', ...
    'Banana', ...
    'Exp', ...
    'MultiExtremal'};
MethodList = { ...
    'Gradient Method (auto step)', ...
    'Steepest Descent', ...
    'Quasi-Newton (BFGH)', ...
    'Quasi-Newton (DFP)', ...
    'Simplex Search', ...
    'Trust Region Newton', ...
    'Trust Region Newton (Hessian Provided)'};
nMethods = length(MethodList);

ProblemInit = 1;
fInit = '100*(y-x^2)^2+(1-x)^2';
MethodInit = 1;
x0Init = -1;
y0Init = 2;
x1Init = -3;
x2Init = 3;
y1Init = -4;
y2Init = 5;
ContourInit = '[3, 10, 50, 100, 250, 500, 1000, 2500, 5000, 10000]';
ContourTextInit = 0;
SurfaceInit = 0;

defaultlist = {NaN, NaN, NaN, '.', 'MarkerSize', 10, 'MarkerEdgeColor'};
Plot(1) = plot3(defaultlist{:}, 'k');
Plot(2) = plot3(defaultlist{:}, 'b');
Plot(3) = plot3(defaultlist{:}, [.25, 0, .5]);
Plot(4) = plot3(defaultlist{:}, 'r');
Plot(5) = plot3(defaultlist{:}, 'm');
Plot(6) = plot3(defaultlist{:}, [0, .5, .25]);
Plot(7) = plot3(defaultlist{:}, [.5, 0, .25]);
defaultlist = {NaN, NaN, NaN, '*', 'MarkerSize', 10, 'LineWidth', 3, 'MarkerEdgeColor'};
xyminPlot(1) = plot3(defaultlist{:}, 'k');
xyminPlot(2) = plot3(defaultlist{:}, 'b');
xyminPlot(3) = plot3(defaultlist{:}, [.25, 0, .5]);
xyminPlot(4) = plot3(defaultlist{:}, 'r');
xyminPlot(5) = plot3(defaultlist{:}, 'm');
xyminPlot(6) = plot3(defaultlist{:}, [0, .5, .25]);
xyminPlot(7) = plot3(defaultlist{:}, [.5, 0, .25]);
x0y0 = plot3(NaN, NaN, NaN, '.b', 'MarkerSize', 25, 'LineWidth', 10);
Contour = 0;
ContourText = 0;
Surface = surf(NaN, NaN, NaN, 'EdgeColor', [.8, .8, .8], 'FaceColor', 'flat');
grid on;

SetInitialValues;
set(H.Problem, 'String', ProblemList, 'Value', ProblemInit, 'CallBack', @ProblemCallBack);

    function SetInitialValues;
        set(H.Axes, 'NextPlot', 'Add');
        set(H.f, 'String', fInit, 'CallBack', @fCallBack);
        set(H.x0, 'String', x0Init, 'CallBack', @x0y0CallBack);
        set(H.y0, 'String', y0Init, 'CallBack', @x0y0CallBack);
        set(H.x1, 'String', x1Init, 'CallBack', @xyLimitsCallBack);
        set(H.x2, 'String', x2Init, 'CallBack', @xyLimitsCallBack);
        set(H.y1, 'String', y1Init, 'CallBack', @xyLimitsCallBack);
        set(H.y2, 'String', y2Init, 'CallBack', @xyLimitsCallBack);
        set(H.Method, 'String', MethodList, 'Value', MethodInit, 'CallBack', @MethodCallBack);
        set(H.Contour, 'String', ContourInit, 'CallBack', @ContourCallBack);
        set(H.ContourText, 'Value', ContourTextInit, 'CallBack', @ContourTextCallBack);
        set(H.Surface, 'Value', SurfaceInit, 'CallBack', @SurfaceCallBack);
        set(H.Clear, 'CallBack', @ClearCallBack);
        fCallBack;
        x0y0CallBack;
        xyLimitsCallBack;
    end;

    function ProblemCallBack(varargin)
        p = ProblemList{get(H.Problem, 'Value')};
        if ~isequal(p, 'Unknown')
            eval(['Problem' p]);
            SetInitialValues;
        end;
    end;

    function fCallBack(varargin)
        f = get(H.f, 'String');
        title(f);
        f = sym(f);
        dfdx = diff(f, 'x');
        dfdy = diff(f, 'y');
        d2fdx2 = diff(dfdx, 'x');
        d2fdydx = diff(dfdx, 'y');
        d2fdxdy = diff(dfdy, 'x');
        d2fdy2 = diff(dfdy, 'y');
        set(H.gradf, 'String', ['[' char(dfdx) '; ' char(dfdy) ']']);
        set(H.Hessf, 'String', ['[' char(d2fdx2) ', ' char(d2fdydx) '; ' char(d2fdxdy) ', ' char(d2fdy2) ']']);
        ClearResults;
        ClearPlots;
        DrawContour;
        DrawSurface;
    end;

    function x0y0CallBack(varargin)
        f = inline([get(H.f, 'String') '+0*x+0*y']);
        x0 = str2num(get(H.x0, 'String'));
        y0 = str2num(get(H.y0, 'String'));
        z0 = feval(f, x0, y0);
        set(x0y0, 'XData', x0, 'YData', y0, 'ZData', z0);
        ClearPlots;
        DrawContour;
    end;

    function xyLimitsCallBack(varargin)
        x1 = str2num(get(H.x1, 'String'));
        x2 = str2num(get(H.x2, 'String'));
        y1 = str2num(get(H.y1, 'String'));
        y2 = str2num(get(H.y2, 'String'));
        set(H.Axes, 'XLim', [x1, x2], 'YLim', [y1, y2]);
        DrawContour;
        DrawSurface;
    end;
    
    function MethodCallBack(varargin)
        options = optimset('Display', 'iter', 'Diagnostics', 'On');
        x0 = [str2num(get(H.x0, 'String')); str2num(get(H.y0, 'String'))];
        f = inline([get(H.f, 'String'), '+0*x+0*y']);
        gradf = inline([get(H.gradf, 'String'), '+zeros(2,1)*x+zeros(2,1)*y']);
        Hessf = inline([get(H.Hessf, 'String'), '+zeros(2)*x+zeros(2)*y']);
        n = get(H.Method, 'Value');
        ClearResults(n);
        ClearPlots(n);
        switch MethodList{n}
            case 'Gradient Method (auto step)'
                set(H.gradf, 'Enable', 'on');
                set(H.Hessf, 'Enable', 'off');
                [x, fval, exitflag, output, grad] = fmingrad(@fun, x0, options, n, 2)
            case 'Steepest Descent'
                set(H.gradf, 'Enable', 'on');
                set(H.Hessf, 'Enable', 'off');
                options = optimset(options, 'LargeScale', 'off', 'GradObj', 'on');
                options = optimset(options, 'HessUpdate', 'SteepDesc');
                [x, fval, exitflag, output, grad] = fminunc(@fun, x0, options, n, 2)
            case 'Quasi-Newton (BFGH)'
                set(H.gradf, 'Enable', 'on');
                set(H.Hessf, 'Enable', 'off');
                options = optimset(options, 'LargeScale', 'off', 'GradObj', 'on')
                options = optimset(options, 'HessUpdate', 'BFGS'); % BFGS-default
                [x, fval, exitflag, output, grad] = fminunc(@fun, x0, options, n, 2)
            case 'Quasi-Newton (DFP)'
                set(H.gradf, 'Enable', 'on');
                set(H.Hessf, 'Enable', 'off');
                options=optimset(options, 'LargeScale', 'off', 'GradObj', 'on');
                options=optimset(options, 'HessUpdate', 'DFP');
                [x, fval, exitflag, output, grad] = fminunc(@fun, x0, options, n, 2)
            case 'Simplex Search'
                set(H.gradf, 'Enable', 'off');
                set(H.Hessf, 'Enable', 'off');
                [x, fval, exitflag, output] = fminsearch(@fun, x0, options, n, 1)
            case 'Trust Region Newton'
                set(H.gradf, 'Enable', 'on');
                set(H.Hessf, 'Enable', 'off');
                options = optimset(options, 'LargeScale', 'on', 'GradObj', 'on', 'Hessian', 'off');
                [x, fval, exitflag, output, grad] = fminunc(@fun, x0, options, n, 2)
            case 'Trust Region Newton (Hessian Provided)'
                set(H.gradf, 'Enable', 'on');
                set(H.Hessf, 'Enable', 'on');
                options = optimset(options, 'LargeScale', 'on', 'GradObj', 'on', 'Hessian', 'on');
                [x, fval, exitflag, output, grad] = fminunc(@fun, x0, options, n, 2)
        end;
        if exitflag > 0
            Iterations = num2str(output.iterations);
            FunctionEval = num2str(output.funcCount);
        else
            Iterations = 'fail';
            FunctionEval = 'fail';
        end;
        xymin = sprintf('(%3.2e, %3.2e)', x(1), x(2));
        fmin = sprintf('%3.2e', fval);
        set(H.(['Iterations', num2str(n)]), 'String', Iterations);
        set(H.(['FunctionEval', num2str(n)]), 'String', FunctionEval);
        set(H.gradf, 'Enable', 'on');
        set(H.Hessf, 'Enable', 'on');
        set(xyminPlot(n), 'XData', x(1), 'YData', x(2), 'ZData', fval);
        set(H.(['xymin' num2str(n)]), 'String', xymin);
        set(H.(['fmin' num2str(n)]), 'String', fmin);
              
        function [fn, grad, hes] = fun(x, n, k)
            fn = feval(f, x(1), x(2));
            if k == 1
                DrawPoint(x, fn, n);
            end;
            if nargout > 1
                grad = feval(gradf, x(1), x(2));
                if k == 2
                    DrawPoint(x, fn, n);
                end;
            end;
            if nargout > 2
                hes = feval(Hessf, x(1), x(2));
                if k == 3
                    DrawPoint(x, fn, n);
                end;
            end;
        end;
        
        function DrawPoint(x, fn, n);
            xlim = get(H.Axes, 'XLim');
            ylim = get(H.Axes, 'YLim');
            if x(1) >= xlim(1) && x(1) <= xlim(2) && x(2) >= ylim(1) && x(2) <= ylim(2)
                XData = [get(Plot(n), 'XData'), x(1)];
                YData = [get(Plot(n), 'YData'), x(2)];
                ZData = [get(Plot(n), 'ZData'), fn];
                set(Plot(n), 'XData', XData, 'YData', YData, 'ZData', ZData);
            end;
            pause(0.001);
        end;
    end;
    
    function ClearCallBack(varargin)
        ClearResults;
        ClearPlots;
    end;
    
    function ContourCallBack(varargin)
        DrawContour;
    end;

    function DrawContour
        if Contour
            delete(Contour);
            Contour = 0;
        end;
        if ~isequal(ContourText, 0)
            delete(ContourText);
            ContourText = 0;
        end;
        f = inline(vectorize([get(H.f, 'String'), '+0*x+0*y']));
        xlim = get(H.Axes, 'XLim');
        ylim = get(H.Axes, 'YLim');
        v = str2num(get(H.Contour, 'String'));
        [X, Y] = meshgrid(linspace(xlim(1), xlim(2), 250), linspace(ylim(1), ylim(2), 250)); 
        F = feval(f, X, Y);
        [C, Contour] = contour3(X, Y, F, v);
        set(Contour, 'EdgeColor', 'k');
        if get(H.ContourText, 'Value')
            ContourText = clabel(C, Contour);
            set(ContourText, 'BackgroundColor', [1, 1, .6], 'Edgecolor',[.7, .7, .7]);
        end;
    end;
    
    function ContourTextCallBack(varargin)
        DrawContour;
    end;
    
    function ClearResults(n)
        if nargin == 1
            set(H.(['Iterations' num2str(n)]), 'String', '');
            set(H.(['FunctionEval' num2str(n)]), 'String', '');
            set(H.(['xymin' num2str(n)]), 'String', '');
            set(H.(['fmin' num2str(n)]), 'String', '');
        else
            for n = 1:nMethods
                ClearResults(n);
            end;
        end;
    end;

    function ClearPlots(n)
        if nargin == 1
            set(Plot(n), 'XData', NaN, 'YData', NaN, 'ZData', NaN);
            set(xyminPlot(n), 'XData', NaN, 'YData', NaN, 'ZData', NaN);
        else
            for n = 1:nMethods
                ClearPlots(n);
            end;
        end;
    end;
    
    function DrawSurface
        f = inline(vectorize([get(H.f, 'String'), '+0*x+0*y']));
        xlim = get(H.Axes, 'XLim');
        ylim = get(H.Axes, 'YLim');
        nx = 2*length(get(gca,'XTickLabel')) - 1;
        ny = 2*length(get(gca,'YTickLabel')) - 1;
        [X, Y] = meshgrid(linspace(xlim(1), xlim(2), nx), linspace(ylim(1), ylim(2), ny));
        F = feval(f, X, Y);
        shift = .03*(max(max(F)) - min(min(F)));
        set(Surface, 'XData', X, 'YData', Y, 'ZData', F - shift);
        hsv2 = hsv;
        hsv3 = [hsv2(11:64, :); hsv2(1:10, :)];
        colormap(hsv3);
    end;
    
    function ProblemBanana
        fInit = '100*(y-x^2)^2+(1-x)^2';
        MethodInit = 1;
        x0Init = -1;
        y0Init = 2;
        x1Init = -3;
        x2Init = 3;
        y1Init = -4;
        y2Init = 5;
        ContourInit = '[3, 10, 50, 100, 250, 500, 1000, 2500, 5000, 10000]';
        ContourTextInit = 0;
        SurfaceInit = 0;
    end;

    function ProblemExp
        fInit = '(x^2+y^2)*exp(x^2-y^2)';
        MethodInit = 1;
        x0Init = -.2;
        y0Init = .5;
        x1Init = -1;
        x2Init = 1;
        y1Init = -2.5;
        y2Init = 2.5;
        ContourInit = '[.03,.1:.1:.8,1:.2:2]';
        ContourTextInit = 0;
        SurfaceInit = 0;
    end;

    function ProblemMultiExtremal
        fInit = '(x^2+y^2)*((x-1)^2+(y-1)^2)*((x+1)^2+(y+1)^2)*((x-1)^2+(y+1)^2)*((x+1)^2+(y-1)^2)';
        MethodInit = 1;
        x0Init = .5;
        y0Init = .7;
        x1Init = -1.2;
        x2Init = 1.2;
        y1Init = -1.2;
        y2Init = 1.2;
        ContourInit = '[.1,1:10]';
        ContourTextInit = 0;
        SurfaceInit = 0;
    end;

end