%% ===== Load CP rank results and reproduce plots =====
% Point to the .mat produced by the saver above
in_path = fullfile('E:\Lorenzo\Final_results_thesis', 'cp_rank_results', 'cp_rank_noCV_pairwiseMedoid_20251110_140224.mat'); % <-- EDIT
S = load(in_path);
results = S.results;

% Shorthands
R_range = results.meta.R_range;

% ---------- Plot 1: Training fit vs R (mean ± std, and best) ----------
figure('Name','Training fit across restarts','Color','w');
hold on;
errorbar(R_range, results.train_fit.mean, results.train_fit.std, '-o','LineWidth',1.5);
plot(R_range, results.train_fit.best, '--s','LineWidth',1.2);
grid on; xlabel('Rank R'); ylabel('Fit (1 - ||X - \hat{X}||_F / ||X||_F)');
legend('mean \pm std (restarts)', 'best restart', 'Location','southeast');
title('Training fit vs rank (across restarts)');

% ---------- Plot 2: Stability (FMS to medoid) ----------
figure('Name','Stability across restarts (FMS to medoid)','Color','w');
hold on;
errorbar(R_range, results.stability.FMS_mean, results.stability.FMS_std, '-o','LineWidth',1.5);
plot(R_range, results.stability.FMS_median, '--s','LineWidth',1.2);
grid on; xlabel('Rank R'); ylabel('FMS to medoid');
legend('mean \pm std', 'median', 'Location','southeast');
title('Pairwise-FMS medoid stability');

% ---------- Plot 3: Collinearity diagnostics on medoid ----------
figure('Name','Collinearity diagnostics','Color','w');
tiledlayout(2,3,'Padding','compact','TileSpacing','compact');
nexttile; plot(R_range, results.collinearity.maxA,'-o'); grid on; title('A (neurons): max cos');
nexttile; plot(R_range, results.collinearity.maxB,'-o'); grid on; title('B (time): max cos');
nexttile; plot(R_range, results.collinearity.maxC,'-o'); grid on; title('C (trials): max cos');
nexttile; plot(R_range, results.collinearity.medA,'-o'); grid on; title('A (neurons): median cos');
nexttile; plot(R_range, results.collinearity.medB,'-o'); grid on; title('B (time): median cos');
nexttile; plot(R_range, results.collinearity.medC,'-o'); grid on; title('C (trials): median cos');

% ---------- Plot 4: Duplicate-components score on medoid ----------
figure('Name','Duplicate-components score','Color','w');
plot(R_range, results.duplicate_score,'-o','LineWidth',1.5); grid on;
xlabel('Rank R'); ylabel('Duplicate score (tri-mode off-diag sim)');
title('Higher suggests over-factoring');

% ---------- Plot 5: Residual structure on medoid ----------
figure('Name','Residual structure','Color','w');
tiledlayout(2,2,'Padding','compact','TileSpacing','compact');
nexttile; plot(R_range, results.residual.frac1_mode1,'-o'); grid on; title('Residual: top SV frac (mode-1)');
nexttile; plot(R_range, results.residual.frac1_mode2,'-o'); grid on; title('Residual: top SV frac (mode-2)');
nexttile; plot(R_range, results.residual.frac1_mode3,'-o'); grid on; title('Residual: top SV frac (mode-3)');
nexttile; plot(R_range, results.residual.time_autoc,'-o'); grid on; title('Residual: avg lag-1 autocorr (time)');

%% ===== Optional helpers =====
% Access medoids or top-10 for a given R:
% ir = find(R_range==5,1);
% P_med = results.models.medoid_per_R{ir};         % ktensor (needs Tensor Toolbox on path)
% P_top = results.models.best10_per_R{ir};         % 1x<=10 cell of ktensor
% fits  = results.models.best10_fits{ir};          % corresponding fits

% You can now reuse your existing plotting utilities on P_med or P_top{:}.
