% time_check_with_linear_interpolation.m  ver 2.0
% Robust time vector validation and correction
% Tom Irvine (enhanced)

% time_check_with_linear_interpolation.m  ver 3.0
% Robust time vector validation and correction
% Tom Irvine (enhanced + hardened)

function [THM, iflag] = time_check_with_linear_interpolation(THM)

iflag = 1;

%% --- Basic validation ---
if isempty(THM) || size(THM,2) < 2
    error('THM must be Nx2 array');
end

%% --- Extract and clean data ---
t_raw = THM(:,1);
y_raw = THM(:,2);

% Remove NaN/Inf rows
valid = isfinite(t_raw) & isfinite(y_raw);
t = t_raw(valid);
y = y_raw(valid);

n = length(t);

if n < 2
    error('Not enough valid samples after cleaning');
end

% Preserve original start time
t0 = t(1);
t = t - t0;

%% --- Detect already uniform time ---
tdif0 = diff(t);
if std(tdif0) < 1e-12
    fprintf('\nTime vector already uniform. No changes made.\n');
    return;
end

%% --- Compute time differences ---
tdif = diff(t);

% Keep only valid positive differences
tdif = tdif(isfinite(tdif) & tdif > 0);

if isempty(tdif)
    error('Invalid time differences (all zero/NaN)');
end

dt_nom = median(tdif);

if ~isfinite(dt_nom) || dt_nom <= 0
    error('Invalid dt_nom computed');
end

%% --- Deviation metrics ---
tdif_full = diff(t);
dev = abs(tdif_full - dt_nom) / dt_nom;

% Thresholds
tol_noise   = 1e-3;   % 0.1%
tol_jitter  = 0.02;   % 2%
tol_dropout = 0.25;   % 25%

% Classification
idx_noise   = dev < tol_noise;
idx_jitter  = dev >= tol_noise  & dev < tol_jitter;
idx_hiccup  = dev >= tol_jitter & dev < tol_dropout;
idx_dropout = dev >= tol_dropout;

%% --- Diagnostics ---
fprintf('\nTime Step Diagnostics:\n');
fprintf('Nominal dt = %.6g sec (%.2f Hz)\n', dt_nom, 1/dt_nom);
fprintf('Noise      : %d\n', sum(idx_noise));
fprintf('Jitter     : %d\n', sum(idx_jitter));
fprintf('Hiccups    : %d\n', sum(idx_hiccup));
fprintf('Dropouts   : %d\n', sum(idx_dropout));

%% --- Detect sample dropouts ---
drop_factor = tdif_full / dt_nom;
idx_missing = drop_factor > 1.5;

if any(idx_missing)
    fprintf('\n*** SAMPLE DROPOUT DETECTED ***\n');
    disp(find(idx_missing(1:min(end,10))));
end

%% --- Detect sample rate change ---
window = min(200, max(5, floor(n/10)));
dt_smooth = movmedian(tdif_full, window, 'omitnan');

ratio_smooth = dt_smooth / dt_nom;

rate_change_flag = (max(ratio_smooth) > 1.05) || (min(ratio_smooth) < 0.95);

if rate_change_flag
    fprintf('\n*** POSSIBLE SAMPLE RATE CHANGE DETECTED ***\n');
end

%% ============================================================
%% --- Decision Logic ---
%% ============================================================

if sum(idx_dropout) > 0 || any(idx_missing)
    
    % --- TRUE DROPOUT CASE ---
    answer = questdlg('Sample dropouts detected. Perform interpolation?', ...
                      'Time Step', 'Yes','No','Yes');
    
    if strcmp(answer,'No')
        iflag = 999;
        warndlg('No output generated');
        return;
    end
    
    fprintf('\nPerforming interpolation to uniform grid (gap filling)...\n');
    
    % Allow size change (fill gaps)
    n_uniform = floor(t(end)/dt_nom) + 1;
    
    if n_uniform < 2
        error('Uniform grid collapsed (check input time vector)');
    end
    
    t_uniform = (0:n_uniform-1)' * dt_nom;
    
    % Interpolation
    y_uniform = interp1(t, y, t_uniform, 'linear', 'extrap');
    
    THM = [t_uniform + t0, y_uniform];

elseif rate_change_flag
    
    % --- SAMPLE RATE CHANGE CASE ---
    warndlg('Sample rate change detected. Recommend segmenting data.');
    
    fprintf('\nWARNING: Data may contain multiple sampling rates.\n');
    fprintf('Consider splitting into segments.\n');
    return

elseif sum(idx_hiccup) > 0
    
    % --- MINOR HICCUPS ---
    fprintf('\nMinor timing hiccups detected. Rebuilding time vector.\n');
    
    % Preserve size
    n = length(t);
    t_uniform = (0:n-1)' * dt_nom;
    
    THM(:,1) = t_uniform + t0;

else
    
    % --- CLEAN DATA ---
    fprintf('\nTime vector is sufficiently uniform. No changes made.\n');

end

%% --- Final sanity check ---
if size(THM,1) < 2
    error('Output THM collapsed unexpectedly');
end

%% --- Optional Save ---

if false
    try
        dlgtitle = 'Save';
        prompt = {'Enter corrected time history output name'};
        dims = [1 60];
        definput = {''};
        
        name = char(inputdlg(prompt,dlgtitle,dims,definput));
        
        if ~isempty(name)
            assignin('base',name,THM);
        end
    catch
        % silent fail
    end
end


fprintf('\n iflag=%d \n',iflag);

sz=size(THM);

fprintf('\n size(THM):  %d %d \n',sz(1),sz(2));

end