%%Code for FITR/Raman Analysis
%% Reading and Variable Assignment for Controls
[Z,txt,num_HDPE] = xlsread("FTIR and Raman data\FTIR controls\HDPE-1 processed.csv");
xset_HDPE = cell2mat(num_HDPE(2:length(num_HDPE),1)); %wave number
yset_HDPE = cell2mat(num_HDPE(2:length(num_HDPE),2)); %intensity (absorbance)
[Z,txt,num_LDPE] = xlsread("FTIR and Raman data\FTIR controls\LDPE-1 processed.csv");
xset_LDPE = cell2mat(num_LDPE(2:length(num_LDPE),1)); %wave number
yset_LDPE = cell2mat(num_LDPE(2:length(num_LDPE),2)); %intensity (absorbance)
[Z,txt,num_PP] = xlsread("FTIR and Raman data\FTIR controls\PP processed.csv");
xset_PP = cell2mat(num_PP(2:length(num_PP),1)); %wave number
yset_PP = cell2mat(num_PP(2:length(num_PP),2)); %intensity (absorbance)
[Z,txt,num_PVC] = xlsread("FTIR and Raman data\FTIR controls\PVC-1 processed.csv");
xset_PVC = cell2mat(num_PVC(2:length(num_PVC),1)); %wave number
yset_PVC = cell2mat(num_PVC(2:length(num_PVC),2)); %intensity (absorbance)
[Z,txt,num_PS] = xlsread("FTIR and Raman data\FTIR controls\PS processed.csv");
xset_PS = cell2mat(num_PS(2:length(num_PS),1)); %wave number
yset_PS = cell2mat(num_PS(2:length(num_PS),2)); %intensity (absorbance)
[Z,txt,num_PET] = xlsread("FTIR and Raman data\FTIR controls\PET-1 processed.csv")
xset_PET = cell2mat(num_PET(2:length(num_PET),1)); %wave number
yset_PET = cell2mat(num_PET(2:length(num_PET),2)); %intensity (absorbance)
[Z,txt,num_PA6] = xlsread("FTIR and Raman data\FTIR controls\PA6 processed.csv")
xset_PA6 = cell2mat(num_PA6(2:length(num_PA6),1)); %wave number
yset_PA6 = cell2mat(num_PA6(2:length(num_PA6),2)); %intensity (absorbance)
%% Treatment (Processed)
[Z,txt,num_CF4_2] = xlsread("FTIR and Raman data\FTIR water samples\Water samples - FTIR processed\CF-4(2) processed.csv"); %CF-1(1), CF-1(2), CF-2, CF-3, CF-4, and CF-4(2).
xset_CF4_2 = cell2mat(num_CF4_2(2:length(num_CF4_2),1)); %wave number
yset_CF4_2 = cell2mat(num_CF4_2(2:length(num_CF4_2),2)); %intensity (absorbance)
[Z,txt,num_DC] = xlsread("FTIR and Raman data\FTIR water samples\Water samples - FTIR processed\DC-1 processed.csv");
xset_DC = cell2mat(num_DC(2:length(num_DC),1)); %wave number
yset_DC = cell2mat(num_DC(2:length(num_DC),2)); %intensity (absorbance)
[Z,txt,num_NR] = xlsread("FTIR and Raman data\FTIR water samples\Water samples - FTIR processed\NR-1 processed.csv");
xset_NR = cell2mat(num_NR(2:length(num_NR),1)); %wave number
yset_NR = cell2mat(num_NR(2:length(num_NR),2)); %intensity (absorbance)
[Z,txt,num_NRCP] = xlsread("FTIR and Raman data\FTIR water samples\Water samples - FTIR processed\NRCP 1 processed.csv");
xset_NRCP = cell2mat(num_NRCP(2:length(num_NRCP),1)); %wave number
yset_NRCP = cell2mat(num_NRCP(2:length(num_NRCP),2)); %intensity (absorbance)
[Z,txt,num_NRSS] = xlsread("FTIR and Raman data\FTIR water samples\Water samples - FTIR processed\NRSS 1 processed.csv");
xset_NRSS = cell2mat(num_NRSS(2:length(num_NRSS),1)); %wave number
yset_NRSS = cell2mat(num_NRSS(2:length(num_NRSS),2)); %intensity (absorbance)
[Z,txt,num_PB] = xlsread("FTIR and Raman data\FTIR water samples\Water samples - FTIR processed\PB1 processed.csv");
xset_PB = cell2mat(num_PB(2:length(num_PB),1)); %wave number
yset_PB = cell2mat(num_PB(2:length(num_PB),2)); %intensity (absorbance)
[Z,txt,num_SHR1b2] = xlsread("FTIR and Raman data\FTIR water samples\Water samples - FTIR processed\SHR 1b(2) processed.csv"); %Changing
xset_SHR1b2 = cell2mat(num_SHR1b2(2:length(num_SHR1b2),1)); %wave number
yset_SHR1b2 = cell2mat(num_SHR1b2(2:length(num_SHR1b2),2)); %intensity (absorbance)
%% Plot

swc = "on";
swc_2 = "off";
swc_PCA = "on";
swc_corr = "off";
if swc == "on"
    figure; hold on
    subplot(2,4,1)
    plot(xset_HDPE,yset_HDPE,'b')
    title("FTIR HDPE Control:")
    grid

    subplot(2,4,2)
    plot(xset_LDPE,yset_LDPE,'g')
    title("FTIR LDPE Control:")
    grid

    subplot(2,4,3)
    plot(xset_PP,yset_PP,'y')
    title("FTIR PP Control:")
    grid

    subplot(2,4,4)
    plot(xset_PVC,yset_PVC,'r')
    title("FTIR PVC Control:")
    grid

    subplot(2,4,5)
    plot(xset_PS,yset_PS,'c')
    title("FTIR PS Control:")
    grid

    subplot(2,4,6)
    plot(xset_PET,yset_PET,'k')
    title("FTIR PET Control:")
    grid

    subplot(2,4,7)
    plot(xset_PA6,yset_PA6, Color=[0.6350 0.0780 0.1840])
    title("FTIR PA6 Control:")
    grid
    hold off
end

if swc_2 == "on"
    figure; hold on
    subplot(2,4,1)
    plot(xset_CF3,yset_CF3,Color=[0 .4470 .7410])
    title("FTIR Sample CF:")
    grid

    subplot(2,4,2)
    plot(xset_DC,yset_DC,Color=[.8500 .3250 .0980])
    title("FTIR Sample DC:")
    grid

    subplot(2,4,3)
    plot(xset_NR,yset_NR,Color=[.9290 .6940 .1250])
    title("FTIR Sample NR:")
    grid

    subplot(2,4,4)
    plot(xset_NRCP,yset_NRCP,Color=[.4940 .1840 .5560])
    title("FTIR Sample NRCP:")
    grid

    subplot(2,4,5)
    plot(xset_NRSS,yset_NRSS,Color=[.4660 .6740 .1880])
    title("FTIR Sample NRSS:")
    grid

    subplot(2,4,6)
    plot(xset_PB,yset_PB,Color=[.3010 .7450 .9330])
    title("FTIR Sample PB:")
    grid

    subplot(2,4,7)
    plot(xset_SHR1b2,yset_SHR1b2, Color=[.7522 .4 .2])
    title("FTIR Sample SHR:")
    grid
    hold off
end

%% PCA Analysis
if swc_PCA == "on"
    x = 1;
    mat_SHR1b2_org = [xset_SHR1b2 yset_SHR1b2]; %Reducing set to 600 to 4000 from 400 to 4000
    while x == 1 
        if mat_SHR1b2_org(x+1,1) < 600
            mat_SHR1b2_org(x,:) = [];
        end
        if mat_SHR1b2_org(x+1,1) >= 600
            break
        end
    end
   
    xset_interp = 600:.482:4000;
    yset_SHR1b2_interp = interp1(mat_SHR1b2_org(:,1),mat_SHR1b2_org(:,2),xset_interp)'; %interpolation for water samples to match array size for controls
    mat_SHR1b2 = [xset_interp' yset_SHR1b2_interp];
    d = [33;60;111;231;815];
    PCA_Matrix = [yset_SHR1b2_interp(d(1)) yset_HDPE(d(1)) yset_LDPE(d(1)) yset_PP(d(1)) yset_PS(d(1)) yset_PET(d(1)) yset_PA6(d(1)) yset_PVC(d(1));
                  yset_SHR1b2_interp(d(2)) yset_HDPE(d(2)) yset_LDPE(d(2)) yset_PP(d(2)) yset_PS(d(2)) yset_PET(d(2)) yset_PA6(d(2)) yset_PVC(d(2));
                  yset_SHR1b2_interp(d(3)) yset_HDPE(d(3)) yset_LDPE(d(3)) yset_PP(d(3)) yset_PS(d(3)) yset_PET(d(3)) yset_PA6(d(3)) yset_PVC(d(3));
                  yset_SHR1b2_interp(d(4)) yset_HDPE(d(4)) yset_LDPE(d(4)) yset_PP(d(4)) yset_PS(d(4)) yset_PET(d(4)) yset_PA6(d(4)) yset_PVC(d(4));
                  yset_SHR1b2_interp(d(5)) yset_HDPE(d(5)) yset_LDPE(d(5)) yset_PP(d(5)) yset_PS(d(5)) yset_PET(d(5)) yset_PA6(d(5)) yset_PVC(d(5))];
%                   yset_SHR1b_interp(d(6)) yset_HDPE(d(6)) yset_LDPE(d(6)) yset_PP(d(6)) yset_PS(d(6)) yset_PET(d(6)) yset_PA6(d(6)) yset_PVC(d(6));
%                   yset_SHR1b_interp(d(7)) yset_HDPE(d(7)) yset_LDPE(d(7)) yset_PP(d(7)) yset_PS(d(7)) yset_PET(d(7)) yset_PA6(d(7)) yset_PVC(d(7));
%                   yset_SHR1b_interp(d(8)) yset_HDPE(d(8)) yset_LDPE(d(8)) yset_PP(d(8)) yset_PS(d(8)) yset_PET(d(8)) yset_PA6(d(8)) yset_PVC(d(8))];
    [coeff,score,latent,tsquared,explained,mu] = pca(PCA_Matrix');
    display(explained)
    display(score)
    display(tsquared)

    % k-means clustering
    opts = statset('Display','final');
    score_matrix = [score(:,1),score(:,2)];
    [idx, C] = kmeans(score_matrix,3,'Distance','cityblock','Replicates',5,'Options',opts);
    
    figure;
    plot(score_matrix(idx==1,1),score_matrix(idx==1,2),'r.','MarkerSize',12)
    hold on
    plot(score_matrix(idx==2,1),score_matrix(idx==2,2),'b.','MarkerSize',12)
    plot(score_matrix(idx==3,1),score_matrix(idx==3,2),'g.','MarkerSize',12)
    plot(C(:,1),C(:,2),'kx','MarkerSize',15,'LineWidth',.5)
%     legend('Cluster 1','Cluster 2','Cluster 3',Centroids','Location','NW')
%     title('K-means Clustering for PCA Scores of SHR-2(2) versus Plastic Controls')
    hold off
%     
%     
    weighted_score = explained(1)*score(:,1) + explained(2)*score(:,2) + explained(3)*score(:,3);
    scatter3(score(:,1),score(:,2),score(:,3),[],weighted_score,'','MarkerEdgeColor','k')
    set(gca, 'Colormap')
    colorbar()

    title("FTIR of SHR-2(2) vs. Plastic Controls")
    xlabel('1st Principal Component (%)')
    ylabel('2nd Principal Component (%)')
    zlabel('3rd Principal Component (%)')
% 
    s_mat = [yset_SHR1b_interp(d(1));yset_SHR1b_interp(d(2));yset_SHR1b_interp(d(3));yset_SHR1b_interp(d(4));yset_SHR1b_interp(d(5))];
    c_mat = [yset_PET(d(1));yset_PET(d(2));yset_PET(d(3));yset_PET(d(4));yset_PET(d(5))];
    [R_SHR2_2,P_SHR2_2] = corrcoef(s_mat,c_mat)

    plot(xset_interp,yset_SHR1b_interp,xset_interp,yset_PA6,'c')
    xlabel("1/cm")
    ylabel("Relative Intensity")
    title("SHR-2(2) vs. PA6 (r = .9421 , p = .0166)")
    xline([xset_interp(d(1)),xset_interp(d(2)),xset_interp(d(3)),xset_interp(d(4)),xset_interp(d(5))],'--k');
    legend("SHR-2(2)","PA6","607.23","620.73","655.43","707.49","745.08")
end