close all; clear; clc;

%% Fit function and Jacobian
f = @(x,y) log(1 + (x.*sign(real(x))).^(2 + sin(x+y))) + y.^2/20;
h = 1e-4; % Finite difference step size
%J = @(x,y) [(f(x+h,y)-f(x-h,y)) (f(x,y+h)-f(x,y-h))]/(2*h);
J = @(x,y) [imag(f(x+1i*h,y))/h imag(f(x,y+1i*h))/h]; % Complex step

%% Plot
% Gradient-descent

tol = 1e-4; alpha = 0.5;
cm = lines;
figure(1); clf;
a = axes; hold on;
x = -1:.1:10; y = -10:.1:1;
[xx,yy] = meshgrid(x,y);
ff = f(xx,yy);
imagesc(x,y,ff);
a.XDir = 'normal'; a.YDir = 'normal';
grid on; box on; axis tight square;
colormap(viridis);
plot(0,0,'yp','MarkerSize',10);

p = [8.3 -8.9]; % Starting point p0
%[x,y] = ginput(1); p = [x,y]; % Manually pick starting point

pold = p;
plot(p(1),p(2),'yo','MarkerSize',10);
tic
for j = 1:1000 % Limit number of iterations to 1e3
    Jxy = J(p(1),p(2));
    plot([p(1) p(1)-Jxy(1)],[p(2) p(2)-Jxy(2)],'-','color',cm(6,:),'LineWidth',1);
    plot([pold(1) p(1)],[pold(2) p(2)],'-','color',cm(2,:),'LineWidth',1);
    drawnow
    pold = p;
    
    s = -alpha*Jxy;
    p = p + s;
    if all(abs(s)<tol); break; end
end
toc

%% Print info
fprintf(' iter = %d\n',j);
fprintf('Minimum at:\n');
fprintf('   pf = (% .2e,% .2e)\n',p);
fprintf(' fval = %.2e\n',f(p(1),p(2)));
fprintf('Actual minimum at:\n');
fprintf('   pf = (% .2e,% .2e)\n',[0,0]);
fprintf(' fval = %.2e\n',0);

