% Put this script in the folder that contains all of
% the OpenChrom txt files and run it.
% Output file: result1.xls.  Enjoy :)

clear;
T_start=datevec(now);
clc;s=dir('*.txt');S=length(s);
PercentColumn=8;
if exist('MATLAB_LOG.txt','file'),diary off;delete('MATLAB_LOG.txt');end;

diary('MATLAB_LOG.txt');
diary on;
disp(['Script started @ ',datestr(T_start),' Script version:',mfilename]);

answer = inputdlg({['Number of files (all-',num2str(S),')'],...
    'RT selfie (min)',...
    'RT tolerance (min)',...
    'MergePeaksTolerance (last merge)',...
    'Merge only with 3 MZ',...
    'Generate unsorted file?'},...
    'Choose please',1,{'319','0.5','1','1','YES','No'})' % Defaults line
if isempty(answer),return;end;

I=1;
S=str2double(answer{I});I=I+1; % number of files to read
RT_Selfie=str2double(answer{I});I=I+1; %minutes
RT_range=str2double(answer{I});I=I+1; %minutes
MergePeaksTolerance=str2double(answer{I});I=I+1;
MergeOnly3MZ=strcmpi(answer{I},'yes');I=I+1;
GenerateUnsorted=strcmpi(answer{I},'yes');I=I+1;

if length(s)<S
    S=length(s);
    disp(['Reducing the number of files to ',num2str(S),' files']);
end

if exist('Result.xls','file') || exist('Result1.xls','file')
    try
        xlswrite('Result',1);        xlswrite('Result1',1);
        delete('Result.xls');        delete('Result1.xls');
    catch
        try
            wbkname = 'Result1.xls';
            h = actxGetRunningServer('Excel.Application');
            h.WorkBooks.Item(wbkname).Close;
            delete('Result.xls');        delete('Result1.xls');
        catch
            questdlg('Unable to write to file, please close Result.xls & Result1.xls','','Ok','Ok');
            delete('Result.xls');        delete('Result1.xls');
        end
    end
end
%% TCAS=zeros(400,3,S);
Integration=zeros(400,S);
N_Peaks=zeros(S,1);
RT=zeros(400,S);MZ=zeros(400,5,S);
for f_ind=1:S
    fid=fopen(s(f_ind).name);
    A=fscanf(fid,'%c');
    fclose(fid);
    %% RT & mz
    a=regexp(A,'PEAK LIST OVERVIEW','once');
    b=regexp(A,'PEAK INTEGRATION RESULTS','once');
    tmp=regexp(A(a:b),'\[\d+\][^\r]+','match');
    for k=1:length(tmp)
        RT_tmp=regexp(tmp{k},'\d+\.\d+','match','once');
        RT(k,f_ind)=str2num(RT_tmp); %#ok
        mz_tmp=regexp(tmp{k},'[\d\.\s]+\.\.\.','match','once');
        mz_tmp2=str2num(mz_tmp(1:end-3)); %#ok
        MZ(k,:,f_ind)=0;
        MZ(k,1:length(mz_tmp2),f_ind)=mz_tmp2;
    end
    
    %% Find the peaks
    a=regexp(A,'PEAK INTEGRATION RESULTS','once');
    b=regexp(A,'SCAN IDENTIFICATION RESULTS','once');
    tmp=regexp(A(a:b),'\d+\.\d+','match');
    I_tmp=str2double (tmp);
    Integration(1:length(I_tmp),f_ind)=I_tmp;
    N_Peaks(f_ind)=length(I_tmp);
    %     %% Identify the gas for the peaks
    %     a=regexp(A,'PEAK IDENTIFICATION RESULTS','once');
    %     CAS=zeros(N_Peaks(f_ind),3);
    %     for k=1:N_Peaks(f_ind)
    %         a1=regexp(A(a:end),['\[',num2str(k),'\][^\r]+'],'match');
    %         tmp=regexp(a1,'\d+-\d\d-\d\s+\d+\.\d','match','once');
    %         tmp1=regexp(tmp,'[\d-\.]+','match');
    %         tmp3=[tmp1{:}];
    %         if length(tmp3)<6
    %             tmp3(length(tmp3)+1:6)={''};
    %         end
    %         CAS(k,1)=str2double(regexprep(tmp3{1},'-',''));
    %         CAS(k,2)=str2double(regexprep(tmp3{3},'-',''));
    %         CAS(k,3)=str2double(regexprep(tmp3{5},'-',''));
    %     end
    %     TCAS(1:N_Peaks(f_ind),:,f_ind)=CAS;
    %     disp(num2str(f_ind));
end

Mega_Line=zeros(10000,S+PercentColumn);
%% Post processing
%First minute mz
StrangeFiles=zeros(1,S+PercentColumn);
for f_ind=1:S
    SkipMe=0;
    Est=[', Time left ',num2str((S-f_ind)*etime(datevec(now),T_start)/f_ind/60,'%.1f'),' min'];
    disp(['Method file: ',s(f_ind).name(1:end-4),', (',num2str(f_ind),')',Est]);
    tic;
    % Arrange MZ
    MZ_To_Unite=MZ(:,1:5,f_ind);
    
    Line=zeros(N_Peaks(f_ind),S+PercentColumn);
    
    Counter1=0;
    for ll=1:N_Peaks(f_ind)
        Time=RT(:,f_ind);
        CurrentRT=RT(ll,f_ind);
        if CurrentRT==0 %30-6-2016
            SkipMe=1;
            continue;
        end
        
        CurrentMZ=sort(MZ(ll,1:5,f_ind));
        
        tmp=Time-(CurrentRT-RT_Selfie);tmp(tmp<0)=NaN;% 29-6-16
        [~,I_min]=min(tmp);
        
        Relevant_MZ=MZ_To_Unite(I_min:ll-1,:);
        if ~isempty(Relevant_MZ)
            Relevant_ZeroLine=Line(I_min:ll-1,1);
            Relevant_MZ1=Relevant_MZ;
            Relevant_MZ1(Relevant_MZ1==0)=0.001; % for m/z =0
            for k=1:5
                Relevant_MZ1(Relevant_MZ==CurrentMZ(k))=0;
            end
            Relevant_MZ2=Relevant_MZ1==0;
            Relevant_MZ3=sum(Relevant_MZ2,2);
            EqualPeaks=(Relevant_MZ3>2) & (Relevant_ZeroLine>0);
            
            if any(EqualPeaks)
                I_TO_MERGE=I_min+find(EqualPeaks,1)-1;
                Line(I_TO_MERGE,f_ind+PercentColumn)=Line(I_TO_MERGE,f_ind+PercentColumn)+Integration(ll,f_ind)/1e6;
                Counter1=Counter1+1;
                
                disp(['Selfie merging Retension time: ',num2str(CurrentRT),' with ',...
                    num2str(Time(I_TO_MERGE)'),', MZ: ',num2str(CurrentMZ)]);
                continue;
            end
        end
        
        CurrentInt=Integration(ll,f_ind)/1e6;
        Int_Others=zeros(1,S);
        for k=f_ind+1:S
            I=(RT(:,k)<(CurrentRT+RT_range)).*(RT(:,k)>(CurrentRT-RT_range));
            I(N_Peaks(f_ind):end)=false;
            INDEXES=find(I);
            MZ_tmp=sort(MZ(I>.5,:,k),2);
            Int_Others(k)=0;
            for mz_ind=1:size(MZ_tmp,1)
                
                % first shift
                for jj=1:3
                    if MZ_tmp(mz_ind,1)<CurrentMZ(1)
                        MZ_tmp(mz_ind,:)=[MZ_tmp(mz_ind,2:end),0];
                    end
                end
                
                % middle shift
                for jj=1:2 % shift
                    if MZ_tmp(mz_ind,2)<CurrentMZ(2)
                        MZ_tmp(mz_ind,:)=[MZ_tmp(mz_ind,1),MZ_tmp(mz_ind,3:end),0];
                    end
                end
                
                if sum(I)>0 && all(MZ_tmp(mz_ind,1)==CurrentMZ(1) & MZ_tmp(mz_ind,2)==CurrentMZ(2) & MZ_tmp(mz_ind,3)==CurrentMZ(3)) %29-6-16
                    Int_Others(k)=Int_Others(k)+Integration(INDEXES(mz_ind),k)/1e6;
                    Integration(INDEXES(mz_ind),k)=0; %8-6-16
                    RT(INDEXES(mz_ind),k)=0; %26-6-16
                end
                
            end
        end
        Int_Others(f_ind)=CurrentInt;
        Line(ll,:)=[CurrentRT,CurrentMZ(1:5),0,1,Int_Others];
        Int_Others=Int_Others*0;
    end
    disp([num2str(Counter1),' of same MZ lines with tolerance were merged in the unsorted file.']);
    toc
    
    %% Detect strange files  26-6-16
    tmp=sum(Line);tmp=tmp>0;
    tmp(PercentColumn+f_ind)=0;
    StrangeFiles=StrangeFiles+tmp;
    %%
    
    if f_ind==1
        Mega_Line(1:N_Peaks(f_ind),:)=Line;
    else
        Mega_Line(sum(N_Peaks(1:f_ind-1))+f_ind:sum(N_Peaks(1:f_ind-1))+f_ind-1+N_Peaks(f_ind),:)=Line;
    end
end

%% 26-6-16
StrangeFilesNums=find(StrangeFiles(PercentColumn+2:end)==0);
if ~isempty(length(StrangeFilesNums))
    disp('The strange files are: ');
    for k=1:length(StrangeFilesNums)
        disp(s(1+StrangeFilesNums(k)).name(1:end-4));
    end
else
    disp('No strange files.');
end
%% Write results to excel

Headers={'RT (min)','MZ1','MZ2','MZ3','MZ4','MZ5','MZ-merge','Percentage'};
for kk=1:S
    Headers(PercentColumn+kk)={s(kk).name(1:end-4)};
end
%%                                 UNSORTED FILE                                                            First without sort
if GenerateUnsorted
    Sorted_Matrix=Mega_Line;
    
    ANY_M=Sorted_Matrix>0;
    S_ANY_M=sum(ANY_M,2);
    II=find(S_ANY_M,1);
    
    Percent=round(sum(ANY_M(:,PercentColumn+1:end),2))/S*100;
    Sorted_Matrix(:,PercentColumn)=Percent;
    
    Sorted_Matrix1=Sorted_Matrix(II:end,:);
    Sorted_Matrix1=[zeros(1,S+PercentColumn);Sorted_Matrix1];
    
    disp('Saving Result.xls...');
    try
        xlswrite('Result',Sorted_Matrix1);
        xlswrite('Result',Headers);
    catch
        disp('Unable to save unsorted file');
    end
end
%% Now with sort

[~,IX] = sort(Mega_Line(:,1));
Sorted_Matrix=Mega_Line(IX,:);

ANY_M=Sorted_Matrix>0;
S_ANY_M=sum(ANY_M,2);
II=find(S_ANY_M,1);

Percent=(sum(ANY_M(:,PercentColumn+1:end),2))/S*100;
Sorted_Matrix(:,PercentColumn)=Percent;

Sorted_Matrix1=Sorted_Matrix(II:end,:);
Mat_To_Unite_vec=Sorted_Matrix1(:,2)+Sorted_Matrix1(:,3)*100+Sorted_Matrix1(:,4)/100;
EqualPeaks=(diff(Mat_To_Unite_vec')==0)';
Counter=0;
for k=1:length(EqualPeaks)
    if EqualPeaks(k)
        Last_k=min(length(EqualPeaks),k+100);
        SameMZ_Size=find(~EqualPeaks(k:Last_k),1)-1;
        if isempty(SameMZ_Size)
            continue;
        end
        Sorted_Matrix1(k,:)=max(Sorted_Matrix1(k:k+SameMZ_Size,:));
        Sorted_Matrix1(k+1:k+SameMZ_Size,:)=0;
        Counter=Counter+1;
    end
end
disp([num2str(Counter),' of consequence same MZ lines were merged in the sorted file.']);

[B,IX] = sort(Sorted_Matrix1(:,1));
Sorted_Matrix2=Sorted_Matrix1(IX,:);
ANY_M=Sorted_Matrix2>0;
S_ANY_M=sum(ANY_M,2);
II=find(S_ANY_M,1);
Sorted_Matrix2=Sorted_Matrix2(II:end,:);

%% Merge with tolerance

MMZZ=Sorted_Matrix2(:,2:6);
Counter=0;Correction=1e-9*(1:size(Sorted_Matrix2,1))'; %26-6-16

for k=1:size(Sorted_Matrix2,1)
    
    Time=Sorted_Matrix2(:,1);
    CurrentRT=Sorted_Matrix2(k,1);
    if CurrentRT==0,continue;end;
    
    tmp=Time-(CurrentRT+MergePeaksTolerance+Correction);tmp(tmp>0)=NaN;
    [~,I_min]=max(tmp); %29-6-16
    
    % Merge 2 same MZ with MergePeaksTolerance
    CurrentMZ=Sorted_Matrix2(k,2:6); %30-6-16
    Relevant_MZ=Sorted_Matrix2(k+1:I_min,2:6);
    Relevant_ZeroLine=Sorted_Matrix2(k+1:I_min,1)==0;
    Relevant_MZ1=Relevant_MZ;
    Relevant_MZ1(Relevant_MZ1==0)=0.001; % for m/z =0
    for Ind=1:5
        Relevant_MZ1(Relevant_MZ==CurrentMZ(Ind))=0;
    end
    Relevant_MZ2=Relevant_MZ1==0;
    Relevant_MZ3=sum(Relevant_MZ2,2).*(~Relevant_ZeroLine);
    if MergeOnly3MZ==1
        I_TO_MERGE=k+find(Relevant_MZ3>2);
    else
        I_TO_MERGE=k+find(Relevant_MZ3>1);
    end
    
    if ~isempty(I_TO_MERGE)
        disp(['Merging Retension times: ',num2str(Sorted_Matrix2(k,1)),' with ',...
            num2str(Sorted_Matrix2(I_TO_MERGE,1)'),', MZ: ',num2str(MMZZ(k,:))]);
        
        Sorted_Matrix2(k,PercentColumn+1:end)=...
            max([Sorted_Matrix2(k,PercentColumn+1:end);Sorted_Matrix2(I_TO_MERGE,PercentColumn+1:end)]);
        Sorted_Matrix2(I_TO_MERGE,:)=0;
        Counter=Counter+1;
        
        % Write the most commom MZ that where merged
        tmp=sort(MMZZ(I_TO_MERGE,:),1);
        if size(tmp,1)>1
            tmp(([diff(tmp);[1 1 1 1 1]])==0)=NaN;
            tmp=sort(tmp);
        end
        for Ind=1:5
            tmp(tmp==CurrentMZ(Ind))=NaN;
        end
        Sorted_Matrix2(k,PercentColumn-1)=mode(tmp(:)); %mode is the most frequent value of the array
    end
end
disp([num2str(Counter),' of lines were merged in the sorted file, with tolerance.']);

%Remove zero lines
Time=Sorted_Matrix2(:,1);
Izero=(Time~=0);
Sorted_Matrix2=Sorted_Matrix2(Izero,:);

% Calculate percentage again
ANY_M=Sorted_Matrix2>0;
S_ANY_M=sum(ANY_M,2);
II=find(S_ANY_M,1);

Percent=round(sum(ANY_M(:,PercentColumn+1:end),2)/S*100);
Sorted_Matrix2(:,PercentColumn)=Percent;

Sorted_Matrix2=[zeros(1,S+PercentColumn);Sorted_Matrix2];

%% Plot graphs
subplot(2,2,1);
plot(Sorted_Matrix2(:,1),Sorted_Matrix2(:,PercentColumn),'o')
title('Detection percentage VS RT');ylim([-3 103]);
subplot(2,2,2);
hist(Sorted_Matrix2(:,PercentColumn),8)
title('Detection percentage histogram');
%%

disp('Saving Result1.xls...');
try
    xlswrite('Result1',Sorted_Matrix2);
    xlswrite('Result1',Headers);
    winopen('Result1.xls');
catch
    for k=1:ceil(size(Sorted_Matrix2,2)/100)
        xlswrite(['Result1_Part_',num2str(k)],...
            Sorted_Matrix2(:,((k-1)*100+1):min(k*100,size(Sorted_Matrix2,2))));
    end
    xlswrite('Result1_Part_1',Headers(1:100));
end

disp('Done');
save('Result1.mat','Sorted_Matrix2','Headers','answer');
diary off;