landmark.asv

clear all
close all
clc

load("traces_uc_for_landmark.mat")   % X: [n_samples x n_features]
X = traces_uc;
sigma = 0.03;       % set based on your previous calibration
YYY   = 4;          % baseline effective rank / min #landmarks

% --- Full kernel entropy ---
D2 = pdist2(X, X).^2;
K_full = exp(-D2 / (2 * sigma^2));
eigvals_full = eig(K_full);
eigvals_full = real(eigvals_full(eigvals_full > 0));
p_full = eigvals_full / sum(eigvals_full);
H_full = -sum(p_full .* log2(p_full));   % Shannon entropy

% --- Low trace baseline (minimum quality) ---
X_low = X(1:YYY,:);
K_low = exp(-pdist2(X_low, X_low).^2 / (2 * sigma^2));
eigvals_low = eig(K_low);
eigvals_low = real(eigvals_low(eigvals_low > 0));
p_low = eigvals_low / sum(eigvals_low);
H_low = -sum(p_low .* log2(p_low));

% --- Sweep landmarks + repeat for multiple trials ---
ms = YYY:5:size(X,1);      % landmark counts
n_trials = 30;             % number of repetitions

H_trials   = zeros(n_trials, numel(ms));  % entropy curves per trial
knee_vec   = zeros(n_trials, 1);          % knee value per trial

rng(0);   % optional: fix seed for reproducibility

for t = 1:n_trials
    H_vals = zeros(size(ms));
    
    for i = 1:length(ms)
        m = ms(i);
        idx = randperm(size(X,1), m);   % sample m rows as landmarks
        X_land = X(idx,:);
        
        C = exp(-pdist2(X,      X_land).^2 / (2 * sigma^2));
        W = exp(-pdist2(X_land, X_land).^2 / (2 * sigma^2));
        
        K_nys = C * pinv(W) * C';
        eigvals = eig(K_nys);
        eigvals = real(eigvals(eigvals > 0));
        p = eigvals / sum(eigvals);
        H_vals(i) = -sum(p .* log2(p));   % Shannon entropy
    end
    
    H_trials(t, :) = H_vals;
    
    % knee of THIS trial (if you want per-trial knees)
    x_norm = (ms - min(ms)) / (max(ms) - min(ms));
    y_norm = (H_vals - min(H_vals)) / (max(H_vals) - min(H_vals));
    diff_curve = y_norm - x_norm;   % or abs(y_norm - x_norm), same as before
    [~, knee_idx_t] = max(diff_curve);
    knee_vec(t) = ms(knee_idx_t);
end

% --- 1. Mean & std across trials ---
H_mean = mean(H_trials, 1);
H_std  = std(H_trials, 0, 1);

% --- 2. Detect knee on the averaged curve ---
x = ms;
y = H_mean;
x_n = (x - min(x)) / (max(x) - min(x));
y_n = (y - min(y)) / (max(y) - min(y));
d = y_n - x_n;                     % vertical distance from diagonal
[~, knee_idx] = max(d);
knee_point   = ms(knee_idx);       % knee #landmarks
knee_entropy = H_mean(knee_idx);   % entropy at the knee

% --- 3. Plot with error bars + full/low lines ---
figure;
errorbar(ms, H_mean, H_std, 'o-', 'LineWidth', 1.5); hold on;
yline(H_full, '--', 'H full');
yline(H_low,  ':', 'H low');

plot(knee_point, knee_entropy, 'rp', ...
     'MarkerSize', 12, 'MarkerFaceColor', 'r');

text(knee_point, knee_entropy, sprintf('  Knee = %d', knee_point), ...
     'VerticalAlignment', 'bottom');

xlabel('# Landmarks');
ylabel('Shannon Entropy');
title('Nyström Entropy Approximation with Error Bars (30 trials)');
grid on;

% (Optional) look at the distribution of knees across trials
figure; 
histogram(knee_vec);
xlabel('Knee (# landmarks)');
ylabel('# Trials');
title('Distribution of Knee Values Across Trials');