% Process data to create a spectrum from the measured_data matrix.
% ====================================================================
% % lasso MATLAB does not work with complex variables.
% % =====================================================
% % tic; B = lasso(real(Gn), real(m_data)); toc
% % B1 = lasso(imag(Gn), imag(m_data));

% normalize data.
m_data = m_data/norm(m_data, 1);       % L1.
% m_data = m_data/norm(m_data);       % L2.

Gn = zeros(samp_points*Nshift, freq_points*band_num);

for jj = 1:freq_points*band_num
    tmp = A1(:, jj).*D2(jj, :);
    Gn(:, jj) = tmp(:);
end

% lambda_max = norm(2*(Gn'*m_data), inf);       % see l1_ls package.
% lambda_lasso_1 = 1e-4*lambda_max;
lambda_lasso_1 = 4e-5;
lambda_lasso_2 = 4e-4;
lambda_ridge = 1e-3;

% find the max singular value of Gn.
% ======================================================
s = max(svd(Gn));

% calculation method.
% ========================
calc_method = 'fista';
% calc_method = 'cvx';
% calc_method = 'tikhonov';
% calc_method = 'ls';

% check that A*D = Gn*Rho.
% ============================
n_data = Gn*Rho;
n_data = n_data(:);
disp([' deviation of A*D from Gn*Rho = ', num2str(max(abs(m_data(:) - n_data)), 5), '. ']);

switch calc_method
    case 'ls'
        % calculate Rho using least squares.
        % ===============================================
        Res_ls = Gn\m_data;
        % Res_ls = pinv(Gn)*m_data;
        Res = Res_ls;

    case 'tikhonov'
        % calculate Rho using Tikhonov regularization.
        % ================================================
        Res_tik = pinv(Gn'*Gn + lambda_ridge*s*eye(size(Gn, 2)))*(Gn'*m_data);
        Res = Res_tik;

    case 'cvx'
        cvx_begin
        cvx_precision low
        variable Res_cvx(freq_points*band_num) complex
        minimize norm(Gn*Res_cvx - m_data) + lambda_lasso_1*(norm(Res_cvx, 1));       % modified lasso with norm2 power 1.
        % minimize 0.5*pow_pos(norm(Gn*Res_cvx - m_data), 2) + lambda_lasso_2*s*(norm(Res_cvx, 1));       % true lasso with norm2 power 2.
        % minimize norm(Gn*Res_cvx - m_data) + lambda_ridge*s*norm(Res_cvx) + lambda_lasso_1*s*(norm(Res_cvx, 1));       %  elastic net regression.
        cvx_end
        Res = Res_cvx;

    case 'fista'
        % opts.pos    = false;
        % opts.lambda = s*lambda_lasso_1;
        opts.lambda = lambda_lasso_1;
        opts.tol = 2e-6;
        GtG = Gn'*Gn;
        L = max(eig(GtG));
        % L = norm(Gn, 'fro');       % approximate (||G||op)^2. See link 2.
        Xinit = zeros(size(Gn, 2), 1);
        tStart = tic;
        [Res_fista, min_cost, cost_vec, err] = lasso_fista(m_data, Gn, Gn'*m_data, 2*L, Xinit, opts);
        disp([' t_fista = ', num2str(toc(tStart)), ' sec. ']);
        % cost_fista = calc_F(Res_fista);
        % fprintf('cost_fista = %.5s\n', cost_fista);
        Iter = length(cost_vec);
        disp([' iterations = ', num2str(Iter), '. '])
        disp([' err = ', num2str(err), '. '])
        disp([' cost = ', num2str(min_cost), '. '])
        Res = Res_fista;
end

% output_noise_std = std(Res(abs(freq_axis) > abs(max_shift_kHz)));
output_noise_std = std(Res(freq_axis > max(freq_kHz) + 0.2));
index = freq_axis >= bw_display(1) & freq_axis <= bw_display(2);
h = figure;
set(h, 'Name', [' Nshift = ', num2str(Nshift), '. Bands = ', num2str(band_num), '. N1 = ', num2str(N1), ', N2 = ', num2str(N2), '. noise std = ', num2str(output_noise_std, 3)]);
plot(freq_axis(index), real(Res(index)), freq_axis(index), imag(Res(index)), 'm');
title([' calc spect: samp points = ', num2str(samp_points), '. TR points = ', num2str(points_in_TR), '. method = ', calc_method, '. '])
% axis([bw_display(1), bw_display(2), 0, max(abs(Res(index)))*1.1]);
xlabel(' frequency (kHz) ');
grid on;

clear('params');
params.change_vals = change_vals;
params.convolve_lines = convolve_lines;
params.T1 = T1;
params.T1_a = T1_a;
params.T2 = T2;
params.T2_a = T2_a;
params.TR = TR;
params.TR_a = TR_a;
params.number_spectral_lines = NumSimLines;
params.sw = bw_kHz;
params.max_shift_kHz = max_shift_kHz;
params.bw_display = bw_display;
params.Nshift = Nshift;
params.band_num = band_num;
params.FA = FA;
params.FA_a = FA_a;
params.sw = bw_kHz;
% params.max_shift_freq = max_shift_freq;
params.points_in_TR = points_in_TR;
params.N1 = N1;
params.N1_a = N1_a;
params.N2 = N2;
params.spectral_lines_kHz = freq_kHz;
params.RF_phase = RF_phase;
params.RF_phase_a = RF_phase_a;
params.input_noise_std = sigma;
params.output_noise_std = output_noise_std;
if strcmp(calc_method, 'fista')
    params.num_iter = Iter;
    params.cost = min_cost;
    params.err = err;
end
