%% ============================================================
%   DEMONSTRATING THREE METHODS:
%   1. CUR with SVD-based leverage sampling
%   2. Standard SVD rank-k approximation
%   3. CUR using QR column/row pivoting
% ============================================================

close all; clear; clc;

%% ---- Load grayscale image ----
A = double(imread('peppers.bmp'));
% A = imresize(A,[1024,1024],'bicubic');
A = A(:,:,1);     % ensure grayscale
[m, n] = size(A);

%% ---- Rank for approximation ----
k = 200;           % try 20, 50, 100 as well

%% ============================================================
%   1. CUR decomposition using SVD-based leverage scores
% ============================================================

% ---- Compute truncated SVD ----
[U,S,V] = svds(A, k);

% ---- Compute row/column leverage scores ----
rowLev = sum(U.^2, 2);   % m x 1
colLev = sum(V.^2, 2);   % n x 1

% ---- Select top k rows/columns ----
[~, row_idx] = maxk(rowLev, k);
[~, col_idx] = maxk(colLev, k);

% ---- Build C, R, W ----
C = A(:, col_idx);       % important columns
R = A(row_idx, :);       % important rows
W = A(row_idx, col_idx); % intersection

% ---- Compute U matrix ----
U_cur = pinv(W);

% ---- Reconstruct image ----
A_cur_svd = C * U_cur * R;

%% ============================================================
%   2. Standard SVD-based rank-k approximation (best possible)
% ============================================================

A_svd = U * S * V';

%% ============================================================
%   3. CUR decomposition using QR pivoting (no SVD needed)
% ============================================================

% ---- Column selection via QR pivoting ----
[Qc, Rc, Ec] = qr(A, 'vector');
col_idx_qr = Ec(1:k);
C_qr = A(:, col_idx_qr);

% ---- Row selection via QR pivoting of transpose ----
[Qr, Rr, Er] = qr(A', 'vector');
row_idx_qr = Er(1:k);
R_qr = A(row_idx_qr, :);

% ---- Intersection matrix ----
W_qr = A(row_idx_qr, col_idx_qr);

% ---- Reconstruction ----
A_cur_qr = C_qr * pinv(W_qr) * R_qr;

%% ============================================================
%   DISPLAY RESULTS
% ============================================================

figure('Position',[200 200 1500 600]);

subplot(1,4,1);
imshow(uint8(A));
title('Original Image');

subplot(1,4,2);
imshow(uint8(A_svd));
title(sprintf('Rank-%d SVD Approx.', k));

subplot(1,4,3);
imshow(uint8(A_cur_svd));
title(sprintf('Rank-%d CUR (SVD-based)', k));

subplot(1,4,4);
imshow(uint8(A_cur_qr));
title(sprintf('Rank-%d CUR (QR-pivoting)', k));


%% ============================================================
%   PRINT ERRORS
% ============================================================

err_svd     = norm(A - A_svd, 'fro');
err_cur_svd = norm(A - A_cur_svd, 'fro');
err_cur_qr  = norm(A - A_cur_qr, 'fro');

fprintf('Frobenius Errors:\n');
fprintf('   SVD rank-%d:           %f\n', k, err_svd);
fprintf('   CUR (SVD leverage):    %f\n', err_cur_svd);
fprintf('   CUR (QR pivoting):     %f\n', err_cur_qr);