% Copyright (c) 2017, Domenico L. Gatti
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without 
% modification, are permitted provided that the following conditions are 
% met:
% 
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the 
%       distribution
%       
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%% General dependencies
% We always start from the CODE directory and we add to the path
% subdirectories containing various tools described in the book chapters.
addpath(genpath('../GENERAL_SCRIPTS_FUNCTIONS'));
addpath(genpath('../DATABASE'));
% addpath(genpath('../TOOLBOXES'));

%% CHAPTER 5: Projections

% 1.Find the projection matrix Pc onto the column space of A: 
A = [3 6 6;4 8 8]
rank(A)
% Notice that you can't just do the following
P = A*inv(A'*A)*A'
% First we need to identify the column space basis
C_a = orth(A)
P_c = C_a/(C_a'*C_a)*C_a'
% We could also take any of the three columns of A
C_a = A(:,1)
P_c = C_a/(C_a'*C_a)*C_a'

% 2. Find the projection matrix Pr onto the row space of A.
R_a = orth(A')
P_r = R_a/(R_a'*R_a)*R_a'
% or 
R_a = A(1,:)'
P_r = R_a/(R_a'*R_a)*R_a'

% 3. Multiply B = PcAPr. Your answer B will be a little surprising: can you
% explain it?
B = P_c*A*P_r

% 4. Find the projection matrices Pn and Pln that project onto N(A) and
% LN(A).
Pr_size = size(P_r,1)
I = eye(Pr_size)
P_n = I - P_r
Pc_size = size(P_c,1)
I = eye(Pc_size)
P_ln = I - P_c

% 5. Carry out the product Ax = b for x = [7 1 11]' by calculating
% separately the projections of x onto the row and the null space and by
% multiplying those projections by A.
x = [7 1 11]'
x_ra = P_r*x 
x_na = P_n*x
x_ra + x_na

A*x
A*x_ra
A*x_na

% 6. Compound X absorbs between 300 and 400 nm. We have collected several
% absorbance spectra in this region of the spectrum from a highly purified
% sample of the compound at the concentration of 1 mM: we have averaged all
% these independent determinations in order to have a noise-free
% 'reference' spectrum (spectrum A) that can be used to determine
% accurately the concentration of the same compound in other experiments.
% In a different experiment we are measuring the amount of X formed over
% time in an enzymatic reaction. Using a diode array spectrophotometer we
% collect a full spectrum of the product in the 300 to 400 nm range every
% second. Spectra B and C show the spectrum of the product at 5 and 10
% seconds, respectively. The three spectra are contained in the file
% '../DATABASE/Spectra.txt': the 1st column of the array contains the
% wavelengths. First we import the spectra and we plot them against the
% wavelength:
clear, clc
filename = '../DATABASE/Spectra.txt';
delimiter = '\t';
formatSpec = '%f%f%f%f%[^\n\r]';
fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter);
fclose(fileID);
wl = dataArray{:, 1};
A = dataArray{:, 2};
B = dataArray{:, 3};
C = dataArray{:, 4};
clearvars filename delimiter formatSpec fileID dataArray;
Spectra = figure;set(gcf,'Unit','Normalized','Position',[0 0.2 0.5 0.8]);
plot(wl,A,'b',wl,B,'g',wl,C,'r')
legend('spectrum A','spectrum B','spectrum C')
xlabel('wavelength (nm)')
ylabel('OD (absorbance units full scale)')
hold on

% a. Determine the absorption wavelength and extinction coefficient for
% compound X and the concentration of this compound at 5 and 10 seconds
% based on this extinction.
[ec,ec_ind] = max(A)
B(ec_ind)/ec
C(ec_ind)/ec

% b. Calculate the new spectrum B and spectrum C projected onto spectrum A
% and then plot the three spectra.
P_cA = A/(A'*A)*A'
B_to_A = P_cA*B;
C_to_A = P_cA*C;

plot(wl,B_to_A,'.g',wl,C_to_A,'.r')
legend({'spectrum A','spectrum B','spectrum C','fit B','fit C'})
grid on

% c. What are the 'scale factors' required to correct spectrum A? What is
% their relationship to the concentration of compound X at 5 and 10 s?
scale_B = (A'*A)\A'*B % scale A to B
scale_C = (A'*A)\A'*C % scale A to C

% or simply:
scale_A_B = (A'*A)\A'*[B C]
scale_A_B = A\[B C]

B_conc_fit = B_to_A(ec_ind)/ec
C_conc_fit = C_to_A(ec_ind)/ec

% d. Calculate the projections of spectrum B and spectrum C onto the left
% null space of spectrum A.
[nrows,ncols] = size(P_cA)
P_lna = eye(nrows) - P_cA
B_to_A_er = P_lna*B
C_to_A_er = P_lna*C
B_er = B - B_to_A
C_er = C - C_to_A
openvar B_er 
openvar B_to_A_er

for i = 1:ncols
plot([wl(i) wl(i)], [B_to_A(i) B_to_A(i)+B_er(i)], ':k')
plot([wl(i) wl(i)], [C_to_A(i) C_to_A(i)+C_er(i)], ':k')
end
legend({'spectrum A','spectrum B','spectrum C',...
        'fit B','fit C','error B','error C'})

% What is the meaning of the 'projection' operation?    
B_plus_er = B - 1.01*scale_B*A;
B_minus_er = B - 0.99*scale_B*A;
B_er'*B_er
B_plus_er'*B_plus_er
B_minus_er'*B_minus_er

