Help - Solving Fredholm Integrals of the First Kind

fredholm_nmr.m

% Original article: https://ieeexplore.ieee.org/abstract/document/995059
% Any help is appreciated.
%% Defining constants %%
% From Table 1
N1 = 30;   % T1 pulses
N2 = 4000; % T2 pulses
% From Page 1018, right before part II
Nx = 100;  % Num of discretized T1
Ny = 100;  % Num of discretized T2
%% Defining ranges %%
% From Table 1
tau1s = logspace(log10(5E-4), log10(4), N1);  % Discretized times
tau2s = 2E-4 : 2E-4 : N2 * 2E-4;              % Discretized times
xs = logspace(log10(0.05), log10(0.4), Nx);   % Discretized T1 relaxation time
ys = logspace(log10(0.05), log10(0.4), Ny);   % Discretized T2 relaxation time
%% Kernels
kern1 = zeros(N1, Nx);
kern2 = zeros(N2, Ny);
for ii = 1:N1
    kern1(ii, :) = kernel1(tau1s(ii), xs);  % Function defined at the end
end

for jj = 1:N2
    kern2(jj, :) = kernel2(tau2s(jj), ys);
end
%% Visualizing the Kernels
figure('Name', 'Visualizing the Kernels');
subplot(1, 2, 1);
plot(kern1)
subplot(1, 2, 2);
plot(kern2)
%% Generating model A (Table 1)
[X, Y] = meshgrid(xs, ys);
mux = 0.2;
muy = 0.1;
sigma = 0.01;
Ft = exp( -( (X - mux) .^ 2 + (Y - muy) .^ 2) ./ (2 * sigma .^ 2 ));  % (Nx x Ny)
%% Visualizing model A
figure('Name', 'Visualizing the Joint Probability function');
contourf(Y, X, Ft);
ax = gca();
ax.XScale = 'log';
ax.YScale = 'log';
% Stretched because of the log scale. Why is his not stretched? Different
% sigma?
%% Calculating the perfect signal matrix (Eq. 4)
S = kern1 * Ft * kern2';  % (N1 x N2)
%S = S ./ max(S, [], 'all');  % Normalizing doesn't change much.
%% Visualizing the signal matrix
figure('Name', 'Visualizing the ideal signal');
contourf(S);
colormap parula;
%% Generate something similar to Fig 1 C, but no SVD yet
hold on
plot(tau2s, S(1, :))
plot(tau2s, S(15, :))
plot(tau2s, S(19, :))
plot(tau2s, S(30, :))
hold off
%% Adding some normal noise to the signal matrix
M = S + randn(size(S)) * 2;
figure('Name', 'Visualizing the noisy signal');
contourf(M);
colormap parula;
%% Generate something similar to Fig 1 A
hold on
plot(tau2s, M(1, :))
plot(tau2s, M(15, :))
plot(tau2s, M(19, :))
plot(tau2s, M(30, :))
hold off
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                              Step 1                                 %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating the SVD of both kernels
[U1, S1, V1] = svd(kern1);
[U2, S2, V2] = svd(kern2);
%% Visualizing the singular Values (Fig 1B)
hold on
plot(diag(S1(1:4, 1:4)));
plot(diag(S2(1:10, 1:10)));
ax = gca();
ax.YScale = 'log';
hold off
%% Defining the number of singular values to use. Should be enough.
s1 = 4;
s2 = 10;
%% Regenerating the spectrum from SVD (Eq. 12)
Mtiltil = U1(1:N1, 1:s1) * U1(1:N1, 1:s1)' * M * U2(1:N2, 1:s2) * U2(1:N2, 1:s2)'; % (N1 x N2)
figure('Name', 'Visualizing reconstructed signal');
contourf(Mtiltil);
colormap parula;
%% Generating Fig1 C
hold on
plot(tau2s, Mtiltil(1, :))
plot(tau2s, Mtiltil(15, :))
plot(tau2s, Mtiltil(19, :))
plot(tau2s, Mtiltil(30, :))
hold off
%% Generating the compressed data (Eq. 14)
Mtil = U1(1:N1, 1:s1)' * M * U2(1:N2, 1:s2);  % (s1 x s2)
%% Visualizing the compressed data (Fig 1-D)
plot(reshape(Mtil, [s1*s2, 1]))
% Looks a bit different, I wonder why? Is it the amount of randomness?
% Different data sizes?
%% Check if the norms are the same (Eq. 18)
norm(Mtil, 'fro') - norm(Mtiltil, 'fro') < 1E-4
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                             Step 2                                  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the compressed Kernels (Eq. 15)
K1til = S1(1:s1, 1:s1) * V1(1:Nx, 1:s1)';  % (s1 x Nx)
K2til = S2(1:s2, 1:s2) * V2(1:Ny, 1:s2)';  % (s2 x Ny)
%% Calculate the tensor product of the kernels (Eq. 20)
K0til = kron(K1til, K2til);  % (s1s2 x NxNy)
%% Lexicographically arrange Mtil into mtil
mtil = reshape(Mtil, s1*s2, 1);
%% Minimization with (inverse?) Newton method (Eq. 31, 32, 33)
% Is this the correct method? From https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization
clc
cr = ones(s1*s2, 1);  % Initial Guess
alpha = 1;  % Initial Guess
gamma = 1; % Step size.
prev_chi_val = 1E8;  % Just to check when it's converging
epsilon = 0.001;
maxiter = 50;

for num = 1:maxiter
   G_calc = G(cr, K0til, Nx, Ny);  % Needed for all 3 functions. Defined at the end
   chi_val = chi(cr, G_calc, alpha, mtil);  % Defined at the end
   D_chi = der1_chi(cr, G_calc, alpha, mtil); % Defined at the end
   DD_chi = der2_chi(G_calc, alpha); % Defined at the end
   cr = cr - gamma .* inv(DD_chi) * D_chi;  % Performing newton step
   fprintf('Iter %d: chi= %f, delta= %f\n', num, chi_val, prev_chi_val - chi_val);
   if abs(prev_chi_val - chi_val) < epsilon
       break
   end
   prev_chi_val = chi_val;  % Update Guess
end
fprintf('Done\n')
%% Getting fr (Eq. 25)
fr = zeros(Nx*Ny, 1);
for ii = 1:Nx*Ny
    temp = K0til(:, ii)' * cr;
    if temp > 0
        fr(ii) = K0til(:, ii)' * cr;
    end
end
%% Getting Fr
Fr = reshape(fr, [Nx, Ny]);
%% Visualizing new Fr
fig = figure('Name', 'Comparing the Joint Probability functions');
subplot(1, 2, 1)
contourf(Y, X, Ft);
ax = gca();
ax.XScale = 'log';
ax.YScale = 'log';
subplot(1, 2, 2)
contourf(Y, X, Fr);
colorbar();
ax = gca();
ax.XScale = 'log';
ax.YScale = 'log';
% This obviously looks very very wrong...
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                      Function Definitions                           %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Pre Step 1 - Kernel function definition
function kern1 = kernel1(t1, x)
    kern1 = 1 - 2 * exp(-t1./x);
end

function kern2 = kernel2(t2, y)
    kern2 = exp(-t2./y);
end
%% Step 2 - Optimization functions
function out = G(cr, K0til, Nx, Ny)
    temp = zeros(Nx * Ny, Nx * Ny);
    for ii = 1:Nx*Ny
        %for jj = 1:Ny
        %    if ii == jj 
        temp(ii, ii) = heaviside(K0til(1:end, ii)' * cr);
        %    end
        %end
    end
    out = K0til * temp * K0til';
end
%%
function out = chi(cr, G_calc, alpha, mtil)
    out = 0.5 * cr' * (G_calc + alpha .* eye(size(G_calc))) * cr - cr' * mtil;
end
%%
function out = der1_chi(cr, G_calc, alpha, mtil)
    out = (G_calc + alpha .* eye(size(G_calc))) * cr - mtil;
end
%%
function out = der2_chi(G_calc, alpha)
    out = G_calc + alpha .* eye(size(G_calc));
end