function [ fv, bestx, iter_num ] = ConjungateGradient(f, x, x0, epsilon, show_detail )
%conjungate gradient method
% Input:
% f - syms function
% x - row cell arrow for input syms variables
% x0 - init point
% epsilon - tolerance(公差)
% show_detail - a boolean value for whether for print details
% Output:
% fv - minimum f value
% bestx - mimimum point
% inter_num - iteration count
%%
syms lambdas %suffix(后缀，下标) s indicates this is a symbol variable
% n is the dimension
n = length(x);
% compute differential of function f stored in cell nf
nf = cell(1, n); % using row cells column cells will result in error
for i = 1 : n
nf{i} = diff(f, x{i});
end
% $\nabla f(x_0)$
nfv = subs(nf, x, x0);
% init $\nabla f(x_k)$
nfv_pre = nfv;
% init count, k and xv for x value
count = 0;
k = 0;
xv = x0;
%initial search direction
d = -nfv;
% show initial info
if show_detail
fprintf('Initial:\n');
fprintf('f = %s, x0 = %s, epsilon = %f\n\n', char(f), num2str(x0), epsilon)
end
%% loop
while(norm(nfv) > epsilon) %norm函数 计算范数
%% one-dimensional search
%define $x_{k+1} = x_{k} + \lambda d$
xv = xv + lambdas * d;
%define $\phi$\phi$ and do 1-dim search
phi = subs(f, x, xv);
nphi = diff(phi);
lambda = solve(nphi);
%get rid of complex and minus solution
lambda = double(lambda);
if length(lambda) > 1
lambda = lambda(abs(imag(lambda)) < 1e-5); %imag 求复数虚数部分
lambda = lambda(lambda > 0);
lambda = lambda(1);
end
% if $\lambdas$ is too samll, stop iteration
if lambda < 1e-5
break;
end
%% update
% update $x_{k+1} = x_{k} + \lambds d$
xv = subs(xv, lambdas, lambda);
%convert sym to double
xv = double(xv);
%compute the differential
nfv = subs(nf, x, xv);
%update counters
count = count + 1;
k = k + 1;
%compute alpha based on FR formula
alpha = sumsqr(nfv) / sumsqr(nfv_pre);
% show iteration info
if show_detail
fprintf('Iteration: %d\n', count);
fprintf('x(%d) = %s, lambda = %f\n', count, num2str(xv), lambda);
fprintf('nf(x) = %s, norm(nf) = %f\n', num2str(double(nfv)), norm(double(nfv)));
fprintf('d = %s, alpha = %f\n', num2str(double(d)), double(alpha));
fprintf('\n');
end
%update conjungate direction
d = -nfv + alpha * d;
% save the previous $$\nabla f(x_k)$$
nfv_pre = nfv;
% reset the conjungate direction and k if k >= n
if k >= 0
k = 0;
d = -nfv;
end
end % while
%% output
fv = double(subs(f,x,xv));
bestx = double(xv);
iter_num = count;
end
