function result = htp(A, y, s, maxIter)
if nargin < 4
    maxIter = 100; 
end

[~, N] = size(A);
col_norm = vecnorm(A, 2);
A = A ./ col_norm;

x_hat = zeros(N, 1);  
T_prev = [];          

tic;
for k = 1:maxIter
    b = x_hat + A' * (y - A * x_hat);

    [~, indices] = sort(abs(b), 'descend');
    T = indices(1:s);

    A_T = A(:, T);
    x_T = A_T \ y;

    x_hat = zeros(N, 1);
    x_hat(T) = x_T;

    if isequal(sort(T), sort(T_prev))
        break;
    end

    T_prev = T;
end


time_total = toc;

result = cell(3,1);
result{1} = x_hat./col_norm';
result{2} = time_total;
result{3} = T;
end