function [cue_1,cue_2] = get_cueIndexes(behavior, offset, preonset)

    v = cue_times(behavior, offset, preonset);

    isOne = (v == 1);
    % Pad with zeros at both ends so you detect changes at the boundaries
    padOne = [0, isOne, 0];
    dOne = diff(padOne);
    
    % A rising edge (0 to 1) gives a value of +1:
    startOnes = find(dOne == 1);
    % A falling edge (1 to 0) gives a value of -1;
    % Subtract one to adjust the index as diff shifts it one place
    endOnes   = find(dOne == -1) - 1;
    
    isTwo = (v == 2);
    padTwo = [0, isTwo, 0];
    dTwo = diff(padTwo);
    
    startTwos = find(dTwo == 1);
    endTwos   = find(dTwo == -1) - 1;

    cue_1 = [startOnes;endOnes];
    cue_2 = [startTwos;endTwos];
    
    diffs1 = -(startOnes-endOnes);
    min1 = min(diffs1);
    for i=1:length(startOnes)
        cue_1(2,i) = cue_1(2,i) - max(0, diffs1(i)-min1);
    end

    diffs2 = -(startTwos-endTwos);
    min2 = min(diffs2);
    for i=1:length(startTwos)
        cue_2(2,i) = cue_2(2,i) - max(0, diffs2(i)-min2);
    end
end

function cue_idx = cue_times(behavior, offset, preonset)
% CUE_TIMES returns a vector of cue times.
%
%   cue_idx = cue_times(behavior, offset, preonset)
%
%   Inputs:
%       behavior - a struct with fields:
%                   task_runs, dark_runs, frames_per_run, onsets,
%                   offsets, cue_codes, cs_1_code, cs_2_code, framerate
%       offset   - additional time (in seconds) after the cue to include
%       preonset - time (in seconds) before the cue onset to include
%
%   Output:
%       cue_idx - a 1-by-(runs*frames_per_run) vector, where each element
%                 is 0 or a cue number (1 for cs_1, 2 for cs_2).
%
% The function works as follows:
% 1. The total number of runs is computed as task_runs + dark_runs.
% 2. A vector cue_idx is preallocated with zeros.
% 3. For each onset in behavior.onsets:
%    a. Determine the cue number (1 or 2) based on behavior.cue_codes.
%    b. Compute the original (0-based) cue onset as:
%           p_onset = round(behavior.onsets(i)) - round(preonset * behavior.framerate)
%       and adjust to MATLAB’s 1-based indexing:
%           matlab_cue_onset = p_onset + 1
%    c. Compute cue_time (duration in frames) as:
%           cue_time = (round(behavior.offsets(i)) - p_onset) + 1 + round(offset * behavior.framerate)
%    d. For each frame within the cue period, assign cue_number to cue_idx.
%
% Note: This implementation assumes that the offsets and onsets are given in
%       units that, when multiplied by the framerate, yield frame indices.

    % Total number of runs and frames per run
    runs = behavior.task_runs + behavior.dark_runs;
    frames_per_run = behavior.frames_per_run;
    
    % Preallocate the cue index vector (1 x total number of frames)
    cue_idx = zeros(1, runs * frames_per_run);
    
    % Loop over each cue onset
    for i = 1:length(behavior.onsets)
        % Determine cue number based on cue_codes
        if behavior.cue_codes(i) == behavior.cs_1_code
            cue_number = 1;
        elseif behavior.cue_codes(i) == behavior.cs_2_code
            cue_number = 2;
        else
            % If no matching cue code, skip to the next iteration
            continue;
        end
        
        % Compute the original (0-based) cue onset
        p_onset = round(behavior.onsets(i)) - round(preonset * behavior.framerate);
        % Adjust to MATLAB's 1-based indexing:
        matlab_cue_onset = p_onset + 1;
        
        % Get the cue offset (as a frame index)
        cue_off = round(behavior.offsets(i));
        
        % Compute the duration of the cue period (in frames)
        cue_time = (cue_off - p_onset) + 1 + round(offset * behavior.framerate);
        
        % Assign the cue_number to the corresponding frames in cue_idx
        for j = 0:(cue_time - 1)
            idx = matlab_cue_onset + j;
            if idx <= length(cue_idx)
                cue_idx(idx) = cue_number;
            end
        end
    end
end

