% Random Walk on a Point Cloud (Sphere)
%
% Generates a uniform point cloud on a unit sphere, builds a k-nearest-
% neighbor graph, and simulates a random walk on that graph.
%
% Animation shows:
%   1. 3D sphere with walker position and trail
%   2. Visit frequency heatmap on the sphere
%   3. Histogram of visit counts converging toward uniform

function rw_pointcloud()
    %% =====================================================================
    %  1. Generate uniform point cloud on the unit sphere
    %  =====================================================================
    N = 800;                    % number of points
    points = fibSphere(N);      % N x 3

    %% =====================================================================
    %  2. Build k-nearest-neighbor graph
    %  =====================================================================
    k = 10;
    D = pdist2(points, points);         % N x N distance matrix
    [~, sortIdx] = sort(D, 2);
    neighbors = sortIdx(:, 2:k+1);      % N x k (exclude self)

    %% =====================================================================
    %  3. Set up figure
    %  =====================================================================
    fig = figure('Name','Random Walk on Sphere','Color','w',...
                 'Position',[80 60 1400 550]);

    % Left panel (large): 3D walker view
    ax1 = subplot(1,2,1,'Parent',fig);
    hold(ax1,'on');
    hPts = scatter3(ax1, points(:,1), points(:,2), points(:,3), ...
                    10, [0.7 0.7 0.7], 'filled', 'MarkerFaceAlpha', 0.3);
    hWalker = scatter3(ax1, NaN, NaN, NaN, 200, [1 0.84 0], ...
                       'filled', 'MarkerEdgeColor','k', 'LineWidth', 1.5);
    hTrail = plot3(ax1, NaN, NaN, NaN, '-', 'Color', [1 0.5 0 0.6], ...
                   'LineWidth', 2);
    axis(ax1, 'equal');
    set(ax1, 'XLim', [-1.3 1.3], 'YLim', [-1.3 1.3], 'ZLim', [-1.3 1.3]);
    xlabel(ax1,'x'); ylabel(ax1,'y'); zlabel(ax1,'z');
    title(ax1, 'Random Walk on Sphere', 'FontSize', 14);
    view(ax1, 30, 20);
    set(ax1, 'XGrid','on','YGrid','on','ZGrid','on');

    % Top-right: visit frequency heatmap
    ax2 = subplot(2,2,2,'Parent',fig);
    hold(ax2,'on');
    hHeat = scatter3(ax2, points(:,1), points(:,2), points(:,3), ...
                     20, zeros(N,1), 'filled');
    colormap(ax2, hot);
    cb = colorbar(ax2);
    ylabel(cb, 'Visit count');
    axis(ax2, 'equal');
    set(ax2, 'XLim', [-1.3 1.3], 'YLim', [-1.3 1.3], 'ZLim', [-1.3 1.3]);
    title(ax2, 'Visit Frequency', 'FontSize', 12);
    view(ax2, 30, 20);

    % Bottom-right: histogram of visit counts
    ax3 = subplot(2,2,4,'Parent',fig);
    hold(ax3,'on');
    title(ax3, 'Visit Count Distribution', 'FontSize', 12);
    xlabel(ax3, 'Visit count'); ylabel(ax3, 'Number of points');
    set(ax3, 'XGrid','on','YGrid','on');

    % Status text
    iterText = annotation(fig, 'textbox', [0.25 0.01 0.5 0.04], ...
        'String', 'Step: 0', 'FontSize', 12, 'EdgeColor', 'none', ...
        'HorizontalAlignment', 'center', 'FontWeight', 'bold');

    %% =====================================================================
    %  4. Simulate random walk
    %  =====================================================================
    nSteps = 5000;
    updateEvery = 20;
    trailLen = 30;

    % Initialize
    state = randi(N);
    visitCount = zeros(N, 1);
    trailBuf = nan(trailLen, 3);    % circular buffer for trail

    for t = 1:nSteps
        % Record visit
        visitCount(state) = visitCount(state) + 1;

        % Update trail buffer
        idx = mod(t-1, trailLen) + 1;
        trailBuf(idx, :) = points(state, :);

        % Move to random neighbor
        nbrs = neighbors(state, :);
        state = nbrs(randi(k));

        % Animate
        if mod(t, updateEvery) == 0 || t == 1
            % Walker position
            set(hWalker, 'XData', points(state,1), ...
                         'YData', points(state,2), ...
                         'ZData', points(state,3));

            % Trail: reorder circular buffer for plotting
            trailOrder = [idx+1:trailLen, 1:idx];
            trailPts = trailBuf(trailOrder, :);
            % Remove NaN rows at the start
            validRows = ~isnan(trailPts(:,1));
            trailPts = trailPts(validRows, :);
            set(hTrail, 'XData', trailPts(:,1), ...
                        'YData', trailPts(:,2), ...
                        'ZData', trailPts(:,3));

            % Visit frequency heatmap
            set(hHeat, 'CData', visitCount);
            if max(visitCount) > 0
                caxis(ax2, [0 max(visitCount)]);
            end

            % Histogram
            cla(ax3);
            hold(ax3, 'on');
            if max(visitCount) > 0
                edges = 0:max(1, ceil(max(visitCount)/25)):max(visitCount)+1;
                histogram(ax3, visitCount, edges, ...
                          'FaceColor', [0.3 0.6 0.9], 'EdgeColor', 'w');
                % Expected mean line
                expectedMean = t / N;
                yl = ylim(ax3);
                plot(ax3, [expectedMean expectedMean], yl, 'r--', ...
                     'LineWidth', 2);
                legend(ax3, {'Visits', sprintf('Expected mean = %.1f', expectedMean)}, ...
                       'Location', 'northeast');
            end
            xlabel(ax3, 'Visit count'); ylabel(ax3, 'Number of points');
            title(ax3, 'Visit Count Distribution', 'FontSize', 12);
            set(ax3, 'XGrid','on','YGrid','on');

            % Rotate 3D views slowly
            az = 30 + t * 0.03;
            view(ax1, az, 20);
            view(ax2, az, 20);

            % Status
            cv = 0;
            if mean(visitCount) > 0
                cv = std(visitCount) / mean(visitCount);
            end
            set(iterText, 'String', sprintf(...
                'Step: %d / %d  |  Mean visits: %.1f  |  CV: %.3f', ...
                t, nSteps, mean(visitCount), cv));

            drawnow;
            pause(0.02);
        end
    end

    %% =====================================================================
    %  5. Save figure
    %  =====================================================================
    saveas(fig, fullfile(fileparts(mfilename('fullpath')), 'rw_pointcloud.png'));

    %% =====================================================================
    %  6. Console summary
    %  =====================================================================
    fprintf('\n=============================================\n');
    fprintf('  Random Walk on Sphere — Summary\n');
    fprintf('=============================================\n');
    fprintf('Points: %d,  k-NN: k=%d\n', N, k);
    fprintf('Total steps: %d\n', nSteps);
    fprintf('---------------------------------------------\n');
    fprintf('Mean visits per point:  %.2f  (expected: %.2f)\n', ...
            mean(visitCount), nSteps/N);
    fprintf('Std of visit counts:    %.2f\n', std(visitCount));
    fprintf('Min / Max visits:       %d / %d\n', ...
            min(visitCount), max(visitCount));
    fprintf('Coeff of variation:     %.4f\n', std(visitCount)/mean(visitCount));
    fprintf('Points never visited:   %d / %d\n', sum(visitCount==0), N);
    fprintf('=============================================\n\n');
end

%% =========================================================================
%  Fibonacci sphere: approximately uniform N points on unit sphere
%  =========================================================================
function pts = fibSphere(N)
    golden = (1 + sqrt(5)) / 2;
    i = (0:N-1)';
    theta = 2 * pi * i / golden;        % azimuthal angle
    phi = acos(1 - 2*(i+0.5)/N);        % polar angle
    pts = [sin(phi).*cos(theta), sin(phi).*sin(theta), cos(phi)];
end
