% define data
ts = [0.038     0.194   0.425   0.626   1.253   2.500   3.740];
ys = [0.050     0.127   0.094   0.2122  0.2729  0.2665  0.3317];

xs = zeros(2,1);
xs(1) = 0.9 + .1; 
xs(2) = 0.2 - .1; 

maxiters = 5;
line_search_iters = 100;

for iter=1:maxiters
    % build jacobian and residual vectors
    J = build_jacobian(xs,ts);
    resvec = build_residual(xs,ts,ys);

    % whos J resvec
	%solvec = pcg(J'*J,J'*resvec);
	%solvec = J'*J\(J'*resvec);
	M = J'*J;
	rhs = J'*resvec;
	
	% PM = LU
	% notice that in LAPACK P is represented by vector I for efficiency
	[L,U,P] = lu(M)
	
	% P*L*U*sol = rhs
	% U*sol = y => L*y = P*rhs => y = forward_sub(L,P*rhs)
	y=forward_sub(L,P*rhs);

	% sol = back_sub(U,y); 
	solvec=back_sub(U,y);

	alpha_best = line_search_one(xs,solvec,line_search_iters,ts,ys);
	alpha_best
	
	% want to find alpha so that ||resvec|| is minimized
    xs = xs - alpha_best*solvec;
    xs
end

% calculate fit curve
fit_pts = zeros(1,length(ys));
for i=1:length(fit_pts)
	t = ts(i);
	fit_pts(i) = evalF(xs,t);
end

figure(1); 
hold on;
plot(ts,ys,'r','linewidth',3);
plot(ts,fit_pts,'g','linewidth',3);
legend('actual','fit');
hold off;

