function X = thomson_problem(n, d)
% Simple attempt at computing n well distributed points on a sphere in R^d.
% 
% This is an example of how Manopt can approximate the gradient and even
% the Hessian of a cost function based on finite differences, even if only
% the cost function is specified without its derivatives.
%
% This functionality is provided only as a help for prototyping, and should
% not be used to compare algorithms in terms of computation time or
% accuracy, because the underlying gradient approximation scheme is slow.
%
% See also the derivative free solvers for an alternative:
% pso and neldermead.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Nov. 1, 2016
% Contributors:
% Change log:
%
%   Xiaowen Jiang Aug. 20, 2021
%       Added AD to compute the egrad and the ehess 

if ~exist('n', 'var') || isempty(n)
    n = 50;
end
if ~exist('d', 'var') || isempty(d)
    d = 3;
end

% Define the Thomson problem with 1/r^2 potential. That is: find n points
% x_i on a sphere in R^d such that the sum over all pairs (i, j) of the
% potentials 1/||x_i - x_j||^2 is minimized. Since the points are on a
% sphere, each potential is .5/(1-x_i'*x_j).
problem.M = obliquefactory(d, n);
problem.cost = @(X) sum(sum(triu(1./(1-X'*X), 1))) / n^2;

% From Matlab 2021a, computing egrad and ehess via automatic
% differentiation is available. Notice that the function triu is not 
% supported for AD so far. Replace it with ctriu described in the file 
% manoptADhelp.m. Also, in this particular case, 1./(1-X'*X) may contain 
% NaN on the diagonal which can cause numerical issues when computing
% egrad via AD although the cost is not a function of the diagonal
% elements. To avoid this problem, first take the upper triangular part by
% calling ctriu before dot-division.
% problem.cost = @(X) sum(sum(ctriu(1./ctriu((1-X'*X), 1),1))) / n^2;
% problem = manoptAD(problem);

% Attempt to minimize the cost. Since the gradient is not provided, Manopt
% approximates it with finite differences. This is /slow/, since for each
% gradient approximation, problem.M.dim()+1 calls to the cost function are
% necessary, on top of generating an orthonormal basis of the tangent space
% at each iterate.
%
% Note that it is difficult to reach high accuracy critical points with an
% approximate gradient, hence it may be required to set a less ambitious
% value for the gradient norm tolerance.
opts.tolgradnorm = 1e-4;

% Pick a solver. Both work fairly well on this problem.
% X = conjugategradient(problem, [], opts);
X = rlbfgs(problem, [], opts);

% Plot the points on a translucid sphere
if nargout == 0 && d == 3
    [x, y, z] = sphere(50);
    surf(x, y, z, 'FaceAlpha', .5);
    hold all;
    plot3(X(1, :), X(2, :), X(3, :), '.', 'MarkerSize', 20);
    axis equal;
    box off;
    axis off;
end

% For much better performance, after an early prototyping phase, the
% gradient of the cost function should be specified, typically in
% problem.grad or problem.egrad. See the online document at
%
%   http://www.manopt.org
%
% for more information.

end