% code for Exercise 4.3 of Math 496 (text page 85)
% load the data
% N -> S
clear
load colscan
% W -> E
load rowscan
% SW -> NE
load diag1scan
% NW -> SE
load diag2scan
% we will assume that parameter vector is
% slowness of subsquares unrolled by rows
naxis = 16;
n = naxis*naxis
ndiag = 2*naxis-1;
mA = 2*naxis
ndiag2 = 2*ndiag;
mB = mA + ndiag2
GA = zeros(mA,n);
for ii = 1:naxis
  GA(ii,ii:naxis:n) = 1;
  GA(ii+naxis,(1+naxis*(ii-1)):naxis*ii) = 1;
end
GB = zeros(ndiag2,n);
for ii = 1:naxis
  GB(ii:ii+naxis-1,1+(ii-1)*naxis:ii*naxis) = eye(naxis);
  GB(ndiag2-naxis+2-ii:ndiag2+1-ii,1+(ii-1)*naxis:ii*naxis) = eye(naxis);
end
GB = [GA;sqrt(2)*GB];
% next the data
dA = [colscan(:,5);rowscan(:,5)];
dB = [dA;diag1scan(:,5);diag2scan(:,5)];
% Now the calculations for A
[UA,SA,VA] = svd(GA);
GAdagger = pinv(GA);
mAdagger = pinv(GA)*dA;
rankGA = rank(GA);
resolveA = VA(:,1:rankGA)*VA(:,1:rankGA)';
conditionA = cond(GA);
conditionAdagger = cond(GAdagger);
% check GA for near perfect resolution
perfectA = find(diag(resolveA) > 0.9999);
% calculations for B
[UB,SB,VB] = svd(GB);
GBdagger = pinv(GB);
mBdagger = pinv(GB)*dB;
rankGB = rank(GB);
resolveB = VB(:,1:rankGB)*VB(:,1:rankGB)';
% check GB for near perfect resolution
perfectB = find(diag(resolveB) > 0.9999);
conditionB = cond(GB);
conditionBdagger = cond(GBdagger);
disp('Part (a)')
disp('Dimension of model null space N(GA):')
disp(n - rankGA)
disp('Dimension of data null space N(GA^T):')
disp(mA - rankGA)
disp('Dimension of model null space N(GB):')
disp(n - rankGB)
disp('Dimension of data null space N(GB^T):')
disp(mB - rankGB)
disp('Part (b)')
disp('Indices of model parameters with (near) perfect resolution for GA:')
disp(perfectA)
disp('Indices of model parameters with (near) perfect resolution for GB:')
disp(perfectB)
disp('We locate their position in the main square. Row numbers:')
disp(floor(perfectB/16)+1);
disp('Column numbers: ');
disp(rem(perfectB,16)+1);
disp('Part (c)')
disp('Condition number of GA:')
disp(conditionA)
disp('Condition number of GB:')
disp(conditionB)
disp('Part (d)')
disp('Condition number of generalized inverse of GA:')
disp(conditionAdagger)
disp('Condition number of generalized inverse of GB:')
disp(conditionBdagger)
% Part (e)
disp('Part (e)')
slowref = 1/3000; % reference slowness
MA = reshape(mAdagger,naxis,naxis)';
figure;
mesh(MA)
title(['GLS model A slowness perturbations ranging from  ',num2str(min(mAdagger)), ' to ', num2str(max(mAdagger))]);
% we locate the position of low points in the main square of A
disp('Location of low points in model for A.')
lowA = find(mAdagger < -2.83e-005)
disp('Row numbers:')
disp(floor(lowA/16)+1);
disp('Column numbers: ');
disp(rem(lowA,16)+1);
disp('Speed at the low points in model A is is at least:')
disp(1/(slowref-2.83e-005))
MB = reshape(mBdagger,naxis,naxis)';
figure;
mesh(MB)
title(['GLS model B slowness perturbations ranging from  ',num2str(min(mBdagger)), ' to ', num2str(max(mBdagger))]);
% we locate the position of low points in the main square of B
disp('Location of low points in model for B.')
lowB = find(mBdagger < -9.3e-005)
disp('Row numbers:');
disp(floor(lowB/16)+1);
disp('Column numbers: ');
disp(rem(lowB,16)+1);
disp('Speed at the low points in model B is is at least:')
disp(1/(slowref-9.3e-005))
% Part (f)
resA = reshape(diag(resolveA),naxis,naxis)';
figure;
mesh(resA);
title('Plot of reshaped resolution diagonal entries for A');
resB = reshape(diag(resolveB),naxis,naxis)';
figure;
mesh(resB);
title('Plot of reshaped resolution diagonal entries for B');
% Part (g):
disp('All we need is a (random) nonzero vector from N(GA) or N(GB).')
rand('seed',0); % for predictability
nullslowA = reshape(VA(:,rankGA+1:n)*rand(n-rankGA,1),naxis,naxis)';
figure;
mesh(nullslowA);
title('Plot of a random nonzero vector from N(GB)');
nullslowB = reshape(VB(:,rankGB+1:n)*rand(n-rankGB,1),naxis,naxis)';
figure;
mesh(nullslowB);
title('Plot of a random nonzero vector from N(GB)');
% Part (h)
figure;
mesh(MA+norm(MA)*nullslowA);
title('Plot of perturbation of generalized least square solution by a null vector of GA')
figure;
grid;
mesh(MB+norm(MB)*nullslowB);
title('Plot of perturbation of generalized least square solution by a null vector of GB')



