function dummy = hw1Driver(viewOpt)
%HW1DRIVER driver for CMSC 764 HW1.
%   dummy = hw1Driver(viewOpt) is the driver for hw1 which is an 
%   excise of modified Newton's method for unconstraint optimization 
%   problem. We use Cholesky strategy to modify the hessian if it's 
%   not s.p.d. The modifying approach is  the Levenberg-Marquardt
%   from least square problem.
%   viewOpt controls the option of views, when viewOpt = 1, we have
%   3-D plots, otherwise we have 2-D contours.
%
%   See also NEWTONMIN.

%  Dianne P. O'Leary oleary@cs.umd.edu 
%  Geping Liu geping@math.umd.edu 
%  09-21-06

FCN = {'fg1','fg2','fg3','fg4','fg51','fg52','fg53'};
H = {'Hf1', 'Hf2','Hf3','Hf4','Hf51','Hf52','Hf53'};
X0 = [[0;0],[0;0],[-1.2;1],[2;-2],[1;-1],[1;-1],[1;-1]];

for i = 1:7
    if i == 1,
        x0 = [1;1;1];
    else
        x0 = X0(:,i);
    end
    
    fcn = FCN{i};
    h = H{i};
        
    disp('************');disp(fcn);disp('*************');
    [x,f,iterCount,X,Z] = newtonmin(fcn,h,x0);
    
    disp('The minimizer is ');disp(x);
    disp('The minimum is ');disp(f);,iterCount
    
    %find convergent rate
    n = length(x);
    err = [];
    sum = 0;
    for j = 1: n
        err = [err;X(j,:) - x(j)];
        sum = sum + err(j,:) .^2;
    end
    errNorm2 = sqrt(sum);
    nonZeroIndex = find(errNorm2);
    logErr = log(errNorm2(nonZeroIndex));
    
    sampleSize = length(logErr);
    if sampleSize == 1,
        disp('Not enought points to determine rates');
    else
        if rem(sampleSize,2) == 0,% even number
            nStart = sampleSize / 2;
        else
            nStart = (sampleSize + 1) / 2;
        end
        rate = polyfit(logErr( nStart : sampleSize -1), ...
                       logErr( nStart+1 : sampleSize ),1);                   
        disp('Convergent rate is '); disp(rate(1));
    end
    disp('condition number :');disp(cond(feval(h,x)));

    
    %plot contour
    if (n == 2)
    
        numGridPts = 100;
        pauseTime = 10^(-8);
        scrsz = get(0,'ScreenSize');
        numContours = 50;
        
        minX1 = min(X(1,:));
        maxX1 = max(X(1,:));
        minX2 = min(X(2,:));
        maxX2 = max(X(2,:));
                
        if minX1 == maxX1,
            minX1 = x0(1) - 1;
            maxX1 = x0(1) + 1;
        end
        
        if minX2 == maxX2,
            minX2 = x0(2) - 1;
            maxX2 = x0(2) + 1;
        end
        
        
        xlin = linspace(minX1, maxX1, numGridPts);
        ylin = linspace(minX2, maxX2, numGridPts);
        [XX,YY] = meshgrid(xlin, ylin);

        [numRow,numCol] = size(XX);
        numElement = numRow * numCol;

        Xtemp = reshape(XX, 1, numElement);
        Ytemp = reshape(YY, 1, numElement);
        Ztemp = feval(fcn, n, [Xtemp;Ytemp]);

        ZZ = reshape(Ztemp, numRow, numCol);        
        
        %2D coutour
        pos2D = scrsz;
        pos2D(2) = scrsz(4) / 2;
        pos2D(3) = scrsz(3) / 2;
        pos2D(4) = scrsz(4) / 2;
%       figure('Position',pos2D);
figure
    
        contour(XX, YY, ZZ, numContours),hold on
        plot(X(1,:), X(2,:), '<-'),hold on,
        plot(x(1), x(2), 'ro', 'MarkerSize', 10);%minimizer
    
        xlabel('x1'), ylabel('x2')
        title(['Contour for ',fcn,' [ minimizer marked by o ]'])
    
        %3D coutour--not required for hw,just to get 3-D views
        if viewOpt == 1,
            azDefault = -37.5000;
            elDefault = 30;
            rotDegree = 360;
        
            pos3D = scrsz;
            pos3D(1) = scrsz(3) / 2;
            pos3D(2) = scrsz(4) / 2;
            pos3D(3) = scrsz(3) / 2;
            pos3D(4) = scrsz(4) / 2;
    
            figure
    
            contour(XX, YY, ZZ, numContours),hold on
            plot(X(1,:), X(2,:), '<-'),hold on,
            plot(x(1), x(2), 'ro', 'MarkerSize', 10);%minimizer
       
            meshc(XX, YY, ZZ), hold on, grid on,
            plot3(X(1,:), X(2,:), Z, 'k-', 'LineWidth',5,'MarkerSize', 5)
            plot3(x(1), x(2), f, 'ro', 'MarkerSize',10);%minimizer
    
            xlabel('x1'), ylabel('x2'), zlabel('f')
            title(['Contour for ',fcn,' path in black [ minimizer marked by o ]'])
    
            %views by rotation    
            for i = 1:rotDegree,
                view(azDefault+i,elDefault);
                pause(pauseTime)
            end     
        end
    end    
end