% Test linearization as EIT reconstruction
%
% Samuli Siltanen May 2025



%% Preliminaries

% Pixel grid size 
n = 9;

% Number of electrodes
L = 4*(n-1);

% Injected current amplitude (milli-Amperes)
I = 1;



%% Calculation of the homogeneous case

% Construct simple conductance matrices
condmatH0 = ones(n,n-1);
condmatV0 = ones(n-1,n);

% Construct system matrix
A0 = FD_EIT_2D_system_matrix(n,condmatH0,condmatV0);

% Initialize voltage matrix
Uadj0 = zeros(L,L-1);

% Loop over adjacent current patterns
for iii = 1:(L-1)
    % Current pattern
    index_curr_in  = iii;
    index_curr_out = iii+1;

    % Construct right-hand side.
    rhs = zeros(L+n^2,1);
    rhs(index_curr_in) = I;
    rhs(index_curr_out) = -I;

    % Solve the FD equation, using pseudoinverse because dim(ker(A))=1
    u = pinv(A0)*rhs;

    % The solution is unique up to a constant. Make it unique by subtracting
    % the mean value of voltages at electrodes
    u = u-mean(u(1:L));

    % Insert solution as a column to the voltage matrix
    Uadj0(:,iii) = u(1:L);
end
% Drop the result to vertical vector
Uadj0 = Uadj0(:);


%% Computation of non-homogeneous case

% condmatH: conductivity values between horizontal neighbors, nx(n-1) matrix
% condmatV: conductivity values between vertical neighbors, (n-1)x n matrix
% The conductivity between boundary pixels and outside is assumed to be one.
condmatH = ones(n,n-1);
condmatV = ones(n-1,n);
% Low conductivity inclusion
condmatH(1,:) = .9;
condmatV(1,:) = .9;
% High conductivity inclusion
condmatH(end,:) = 1.2;
condmatV(end,:) = 1.2;

% Construct system matrix
A = FD_EIT_2D_system_matrix(n,condmatH,condmatV);

% Initialize voltage matrix
Uadj = zeros(L,L-1);

% Loop over adjacent current patterns
for iii = 1:(L-1)
    % Current pattern
    index_curr_in  = iii;
    index_curr_out = iii+1;

    % Construct right-hand side.
    rhs = zeros(L+n^2,1);
    rhs(index_curr_in) = I;
    rhs(index_curr_out) = -I;

    % Solve the FD equation, using pseudoinverse because dim(ker(A))=1
    u = pinv(A)*rhs;

    % The solution is unique up to a constant. Make it unique by subtracting
    % the mean value of voltages at electrodes
    u = u-mean(u(1:L));

    % Insert solution as a column to the voltage matrix
    Uadj(:,iii) = u(1:L);
end
% Drop the result to vertical vector
Uadj = Uadj(:);

% Add noise
Uadj = Uadj + 0.1*randn(size(Uadj))/max(abs(Uadj(:)));


%% Calculate Jacobian

disp('Computing Jacobian')
Jcbn = EIT_FD_demo2D_forward_map_Jacobian_adj(n,.01);
disp('Jacobian computed')


%% Linearized reconstruction


[U,D,V] = svd(Jcbn);

% Apply pseudoinverse without tolerance specification
recon = pinv(Jcbn)*(Uadj-Uadj0);

% Apply pseudoinverse with tolerance specification
recon2 = pinv(Jcbn,.1)*(Uadj-Uadj0);

% Collect results into appropriate matrix form
condmatH_recon = reshape(1+recon(1:n*(n-1)),[n,n-1]);
condmatV_recon = reshape(1+recon((n*(n-1)+1):end),[n-1,n]);
condmatH_recon2 = reshape(1+recon2(1:n*(n-1)),[n,n-1]);
condmatV_recon2 = reshape(1+recon2((n*(n-1)+1):end),[n-1,n]);

% Save results to file
save data/FD_EIT_2D_demoE_recon9x9smallcontrast n condmatH condmatV condmatH_recon condmatV_recon condmatH_recon2 condmatV_recon2

% Take a look
FD_EIT_2D_demoE_recon9x9smallcontrast_plot