function [NMSE, BW] = AdaptivePSWFestimator(H, P, nvec, params, method)

    Ncombiner   = params.Ncombiner; 
    NRx         = params.RxArray.N; 
    K           = params.K; 
    sigma2      = params.sigma2; 

    % Step1: Using HMM-AMP (or anything else) to estimate the k-domain support. 
    
    if method == "Bayesian"
        Kd = 2^8;           % Reduced k-domain grid. 
        
        Rcyc = zeros(Kd); 
        eta = 10;           % correlation length
        for ii = 1:Kd
            for jj = 1:Kd
                x = ii/Kd; y = jj/Kd; 
                t = min(abs(x-y), 1-abs(x-y)); 
                Rcyc(ii, jj) = exp(-t/(2*eta)); 
            end
        end
        
        rnkGrids = params.nkGrids(1:(K/Kd):K, 1); 
        rDict = params.RxArray.dict(:, 1:(K/Kd):K); 
        
        % Nrcomb = floor(Ncombiner/4);
        Nrcomb = Ncombiner;  
    
        % Rh = rDict*Rcyc*rDict';     % This subspace is not good. 
        Rh = eye(NRx);  
        % Rh = besselj(0, pi*toeplitz(0:NRx-1)); 
        [VRh, DRh] = eig(Rh); 
        dRh = diag(DRh); 
        eta_arr = zeros(NRx, 1); 
        var_arr = real(dRh); 
    
        for q = 1:Nrcomb
            [~, imax] = max(var_arr); 
            w = VRh(:,imax);    % a unit vector
            y = sqrt(P)*w'*H + nvec(q); 
    
            vq = var_arr(imax); 
            eq = eta_arr(imax); 
            var_arr(imax) = vq*(sigma2/(P*vq+sigma2)); 
            eta_arr(imax) = eq + (sqrt(P)*vq/(P*vq+sigma2))*(y-sqrt(P)*eq); 
        end
    
        h_est = VRh*eta_arr; 

    elseif method == "APSWF"

        Kd = 2^8;           % Reduced k-domain grid. 
        Rcyc = zeros(Kd); 
        eta = 30;           % correlation length
        for ii = 1:Kd
            for jj = 1:Kd
                x = ii/Kd; y = jj/Kd; 
                t = min(abs(x-y), 1-abs(x-y)); 
                Rcyc(ii, jj) = exp(-t/(2*eta)); 
            end
        end
        
        rnkGrids = params.nkGrids(1:(K/Kd):K, 1); 
        rDict = params.RxArray.dict(:, 1:(K/Kd):K); 


        Nrcomb = floor(Ncombiner/4); 
        randX = (randn([NRx, Nrcomb]) + 1j*randn([NRx, Nrcomb]))/sqrt(2); 
        W = randX ./ vecnorm(randX); 
        y1 = sqrt(P)*W'*H + nvec(1:Nrcomb); 
    
        A = sqrt(P)*W'*rDict; 
        rhest = Rcyc*A'*((A*Rcyc*A'+sigma2*eye(Nrcomb))\y1); 
    
        Intervals = estimateIntervals(abs(rhest).^2, 3); 
        Wpswf = getCombiner(params, Intervals, Ncombiner - Nrcomb, "PSWF", 0.05); 
        y_tilde = Wpswf'*sqrt(P)*H + nvec(Nrcomb+1:end);
        
        y_tot = [y1; y_tilde]; 
        W_tot = [W, Wpswf]; 
        h_est = (sigma2*eye(NRx) + P*(W_tot*W_tot'))\(sqrt(P)*W_tot*y_tot); 


    elseif method == "BT-PSWF"
        % Ki = 8; 
        Ki = floor(Ncombiner/4); 
        interval_endpoints = linspace(-1, 1, Ki+1); 
        
        W1              = zeros(NRx, Ki); 
        FullCodebook    = cell(Ki, 1); 
        CntLargerLambda = zeros(Ki, 1); 

        for idx = 1:Ki
            [DSS, Lambda]           = getDiscreteSlepianSequence(NRx, interval_endpoints(idx), interval_endpoints(idx+1)); 
            W1(:, idx)              = DSS(:,1); 
            FullCodebook{idx}       = DSS; 
            CntLargerLambda(idx)    = sum(Lambda >= 0.1); 
        end
        y1 = sqrt(P)*W1'*H + nvec(1:Ki); 
        
        Energies = abs(y1/sqrt(P)).^2;

        if isfield(params.Estimators, 'BT_PSWF') && isfield(params.Estimators.BT_PSWF, 'eta')
            eta = params.Estimators.BT_PSWF.eta; 
        else
            eta = 2.0;
        end
        pvec = exp(eta*Energies); 
        pvec = pvec/sum(pvec); 
        cpvec = cumsum(pvec); 

        W2 = zeros(NRx, Ncombiner-Ki); 
        ccnt = ones(1, Ki);

        for idx = 1:Ncombiner - Ki
            sel = find(rand() < cpvec, 1); 
            W2(:, idx) = FullCodebook{sel}(:, ccnt(sel));
            ccnt(sel) = ccnt(sel) + 1; 

            if ccnt(sel) > CntLargerLambda(sel)
                pvec(sel) = pvec(sel) * 0.3; 
                pvec = pvec/sum(pvec); 
                cpvec = cumsum(pvec); 
                ccnt(sel) = 1; 
            end
        end

        y2 = sqrt(P)*W2'*H + nvec(Ki+1:end); 
        y_tot = [y1; y2]; 
        W_tot = [W1, W2]; 
        h_est = (sigma2*eye(NRx) + P*(W_tot*W_tot'))\(sqrt(P)*W_tot*y_tot); 
        


    elseif method == "BWest-PSWF"
        % Please design a bandwidth estimation method here. 
        Nrcomb = floor(Ncombiner/4); 
        W1 = (randn([NRx, Nrcomb]) + 1j*randn([NRx, Nrcomb]))/sqrt(2); 
        y1 = sqrt(P)*W1'*H + nvec(1:Nrcomb); 

        % using GPR to do bandwidth estimation. 
        w = estimateBandwidth(y1, W1, P, sigma2); 

        Intervals = {[0, 0, w, 0]}; 
        W2 = getCombiner(params, Intervals, Ncombiner - Nrcomb, "PSWF", 0.05); 
        y2 = W2'*sqrt(P)*H + nvec(Nrcomb+1:end);
        
        y_tot = [y1; y2]; 
        W_tot = [W1, W2]; 
        h_est = (sigma2*eye(NRx) + P*(W_tot*W_tot'))\(sqrt(P)*W_tot*y_tot); 
        
    
    elseif method == "MIMO-BWest-PSWF"
        NTx = params.TxArray.N; 
        Nh = NRx*NTx; 

        if ~isfield(params.Estimators.BWest, "w_mode") || params.Estimators.BWest.w_mode == "GPR" 
            if ~isfield(params.Estimators.BWest, "rcomb_mode") || params.Estimators.BWest.rcomb_mode == "auto"
                Nrcomb = floor(Ncombiner/4);
    
            elseif params.Estimators.BWest.rcomb_mode == "external"
                Nrcomb = params.Estimators.BWest.Nrcomb;    % Number of random combiners for bandwidth estimation  
                
            else
                error('BWest rcomb_mode %s not defined.\n', params.Estimators.BWest.rcomb_mode); 
            end
             
            
            [W1, V1] = getMIMOCombinerPrecoder(params, {}, {}, Nrcomb, "Random", 0.1); 
    
            A1 = zeros(Nrcomb, Nh); 
            y1 = zeros(Nrcomb, 1); 
            for idx = 1:Nrcomb
                y1(idx) = sqrt(P)*(W1(:,idx)')*H*V1(:,idx) + nvec(idx); 
                A1(idx, :) = kron(V1(:,idx).', W1(:,idx)'); 
            end
    
            w = estimateBandwidth2D(y1, W1, V1, A1, P, sigma2); 

        elseif params.Estimators.BWest.w_mode == "fixed" 
            w = params.Estimators.BWest.w; 
            Nrcomb = 0; 
            y1 = []; 
            A1 = []; 

        end
        
        IntervalRx = {[0, 0, w, 0]}; 
        IntervalTx = {[0, 0, w, 0]}; 
        [W2, V2] = getMIMOCombinerPrecoder(params, IntervalRx, IntervalTx, Ncombiner, "PSWF", 0.1); 
        y2 = zeros(Ncombiner-Nrcomb, 1);
        A2 = zeros(Ncombiner-Nrcomb, Nh); 
        for idx = 1:Ncombiner-Nrcomb
            y2(idx) = sqrt(P)*(W2(:,idx)')*H*V2(:,idx) + nvec(idx+Nrcomb); 
            A2(idx, :) = kron(V2(:,idx).', W2(:,idx)'); 
        end
        
        y_tot = [y1; y2]; A_tot = [A1; A2]; 
        h_est = (sigma2*eye(Nh) + P*(A_tot'*A_tot)) \ (sqrt(P)*A_tot'*y_tot); 
        H_est = reshape(h_est, [NRx, NTx]); 

        if nargout == 2
            BW = w; 
        end

    else
        error('Method %s not implemented.', method); 
    end
    
    if params.isSIMO
        NMSE = norm(h_est-H).^2/norm(H).^2; 
    else
        NMSE = norm(H_est-H, 'fro')^2/norm(H, 'fro')^2; 
    end

end


% function [gpost, Pgpost, gext, Pgext] = ModuleA(y, A, sigmaw2, gpri, Pgpri)
%     % Model: y = A*g + w
%     
%     Pgpost = Pgpri + (A'*A)/sigmaw2;     % Precision matrix addition law. 
%     gpost = Pgpost\(A'*y/sigmaw2 + Pgpri*gpri);
%     
%     Pgext = (A'*A)/sigmaw2; 
% 
% end