clear all;
close all;clc

load casecirr_3030_prepro

% Peak Finding
maxi = 1;
mini = .5;
y1 = min(MZ);
y2 = max(MZ);
p    = (mini-maxi)/(y2-y1);
q    = maxi-y1*p; 

YN = sample_tr;

[m,n]=size(YN);
peaksByMass = [];

for i = 1: n
    slopeSign = diff(YN(:,i))>0;
    slopeSignChange = diff(slopeSign)<0;
    
    h = find(slopeSignChange) + 1; % finds all peaks
    h( YN(h,i) <  p*MZ(h)+q )= []; % eliminate small peaks

for j = 1: length(h)
           a    =   h(j);
           b    =   MZ(a);
           c    =   YN(a,i);
           f    =   [i a b c c c];
           peaksByMass  =   [peaksByMass; f];
    end
end

% save PeakMass2 peaksByMass

% %%
% clear all;
% close all;clc
% 
% load casecirr_3030_prepro
% load PeakMass2.mat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%     PEAKSET     %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PARAMETERS
minimumSN = 4;
minSN_end = 1;

secondSN = 1;
secSN_end = .5;

y1 = min(MZ);
y2 = max(MZ);

p1  = (minSN_end-minimumSN)/(y2-y1);
q1  = minimumSN-y1*p1;

p2  = (secSN_end-secondSN)/(y2-y1);
q2  = secondSN-y1*p2;

% peaksByMass 
for i = 1: length(peaksByMass)
    peaksByMass(i,4) = peaksByMass(i,4) + (minimumSN - (p1*peaksByMass(i,3)+q1));
    peaksByMass(i,5) = peaksByMass(i,5) + (secondSN - (p2*peaksByMass(i,3)+q2));
end

tickSeparation = 7;
massSeparation = 0.0009;
tic
[x,y] = peakSet_R(peaksByMass, minimumSN, secondSN, tickSeparation, massSeparation);
toc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Remove windows with less that 5 spectra having a peak in that window
indr =find(x(:,6) <= 5);
tempy = [];
tempx = [];
for i = 1: length(indr)
    ty = find(y(:,1) == indr(i));
    tempy =[tempy; ty];
    
    tx = find(x(:,1) == indr(i));
    tempx =[tempx; tx];
end
y(tempy,:)=[];
x(tempx,:)=[];

% Extract windows
    m        =   length(unique(y(:,1)));
    kk       =   unique(y(:,1));
    windows  =   [];
    win_tick =   [];
    m_in     =   [];
    for k = 1: m
        ind = find(y(:,1)==kk(k));
        v   = y(ind,4);
        vv  = y(ind,3);
        windows  = [windows; min(v) max(v)];
        win_tick = [win_tick;min(vv) max(vv)];
    end

save Windows x y windows win_tick ;

%%
 
% PLOT WINDOWS
YN = sample_tr;
figure();
hold on;
for i=1:m;
    maa = windows(i,1);
    plot([maa,maa],[0 99],'LineWidth',1.5,'Color',[.9 0.5 0.5]);
    while maa < windows(i,2)
     maa = maa+.1;
     plot([maa,maa],[0 99],'LineWidth',1.5,'Color',[.9 0.5 0.5]);  
    end
end
 
% Plot spectrum
plot(MZ,YN,'LineWidth',1,'Color',[0.2 0.2 0.2]); hold on;
title('Peak Alignment','FontSize',26);
xlabel('m/z','FontSize',24);
ylabel('Normalized Intensity','FontSize',24);
axis([919 9980 -2 100]);

ss  = YN;
lab = label_tr;
ca = mean(ss(:,find(lab==1)),2);
co = mean(ss(:,find(lab==2)),2);
cr = mean(ss(:,find(lab==3)),2);

plot(MZ,ca,'r'); hold on;
plot(MZ,co,'g'); hold on;
plot(MZ,cr,'b'); 
axis([919 9980 -2 100]);