function extend_chr_sim(xplot, yplot, varargin)
% Diagonal runs are good matches, we need to "extend them a size w to match
% how the alignment would actually look...

% extendDiagonalRuns finds diagonal runs (slope ≈ ±1) in scatter data,
% extends them, and plots them.
%
% Usage:
%   extendDiagonalRuns(xplot, yplot)
%   extendDiagonalRuns(xplot, yplot, 'minRunLength', 10, 'gapTolerance', 3, 'extendLength', 150, 'tol', 1)
%
% Inputs:
%   xplot, yplot   - vectors of scatter points
%
% Name-value parameters:
%   'minRunLength' - minimum number of consecutive points to consider a run (default = 10)
%   'gapTolerance' - allowed gaps in x (default = 3)
%   'extendLength' - how many points to extend a run (default = 150)
%   'tol'          - tolerance for slope detection (default = 1)
%
% Output:
%   A scatter plot with detected runs highlighted and extended.

    % Parse inputs
    p = inputParser;
    addParameter(p, 'minRunLength', 10, @isscalar);
    addParameter(p, 'gapTolerance', 3, @isscalar);
    addParameter(p, 'extendLength', 150, @isscalar);
    addParameter(p, 'tol', 1, @isscalar);
    parse(p, varargin{:});
    params = p.Results;

    % Scatter plot

    scatter(xplot, yplot, 5, 'blue', 'filled'); hold on;

    % Sort points by x
    [sortedX, idx] = sort(xplot);
    sortedY = yplot(idx);
    sortedX = sortedX(:);  % force column
    sortedY = sortedY(:);  % force column

    % Detect slope ~ +1 (y - x ≈ const)
    detectAndExtend(sortedX, sortedY, sortedY - sortedX, ...
        params, 'r');

    % Detect slope ~ -1 (y + x ≈ const)
    detectAndExtend(sortedX, sortedY, sortedY + sortedX, ...
        params, 'g');

    hold off;
    %xlabel('X');
    %ylabel('Y');
    %title('Diagonal runs with extensions');
end

%% Helper: detect and extend runs
function detectAndExtend(sortedX, sortedY, d, params, color)
    [groups,groupIdx] = groupByAbsTol(d, params.tol); 
    for g = 1:length(groups)
        inGroup = find(groupIdx == g);
        if length(inGroup) >= params.minRunLength
            runs = splitRuns(sortedX(inGroup), params.gapTolerance);
            for r = 1:numel(runs)
                runIdx = inGroup(runs{r});
                if length(runIdx) >= params.minRunLength
                    % Plot run
                    %plot(sortedX(runIdx), sortedY(runIdx), [color 'o'], 'LineWidth', 1.5);

                    % Extend line
                    x_end = sortedX(runIdx(end));
                    y_end = sortedY(runIdx(end));
                    if strcmp(color, 'r') % slope ~ +1
                        xe = x_end:(x_end+params.extendLength);
                        ye = y_end + (0:params.extendLength);
                    else % slope ~ -1
                        xe = x_end:(x_end+params.extendLength);
                        ye = y_end - (0:params.extendLength);
                    end
                    xe(xe<0)=NaN;
                    ye(ye<0)=NaN;
                    plot(xe, ye, [color '.'], 'LineWidth', 1.5);
                end
            end
        end
    end
end

%% Helper: split into runs with gap tolerance
function runs = splitRuns(xvals, gapTol)
    diffs = diff(xvals);
    breaks = find(diffs > gapTol);
    startIdx = [1, breaks'+1];
    endIdx   = [breaks', length(xvals)];
    runs = arrayfun(@(s,e) s:e, startIdx, endIdx, 'UniformOutput', false);
end
%% Helper: uniquetol was not working
function [groups, groupIdx] = groupByAbsTol(values, tol)
    values = values(:);
    groups = [];
    groupIdx = zeros(size(values));
    groupCount = 0;
    for i = 1:numel(values)
        if groupIdx(i) == 0
            groupCount = groupCount + 1;
            ref = values(i);
            closeIdx = abs(values - ref) <= tol;
            groupIdx(closeIdx) = groupCount;
            groups(end+1,1) = ref; %#ok<AGROW>
        end
    end
end
