%% This is an open source software

function [dist, alldist] = dist_ptm_20161108( patterninfoS1,patterninfoS2,oz )
% calculate distance between PTPs of given two series

%   Input paremeters
%   patterninfoS1 --- pattern info of first series (index, lenght prob freq)
%   patterninfoS2 --- pattern info of first series (index, lenght prob freq)

%   Output paremeters
%   distmatrix    --- vector of distances between probs1 and probs2
%   overallmatrix --- aggregated sum of elements of distmatrix
%   freqdists     --- differences between frequencies of observed patterns

%   Akbilgic 20150801

%% set parameters and initialize
alldist=zeros(oz.maxmem,1);% one distance for each pattern size

%% calculate distances
for i=1:oz.maxmem
%     i
    myprobs1=patterninfoS1(patterninfoS1(:,2)==i,3:end-1);
	myprobs2=patterninfoS2(patterninfoS2(:,2)==i,3:end-1);            
    
    myfreqs1=patterninfoS1(patterninfoS1(:,2)==i,end);% frekanslari al
    myfreqs2=patterninfoS2(patterninfoS2(:,2)==i,end);% agirlikli frekanslari al
    
    w_myfreqs1=myfreqs1/sum(myfreqs1);% frekanslari al
    w_myfreqs2=myfreqs2/sum(myfreqs2);% agirlikli frekanslari al
    freqdistance= abs(w_myfreqs1- w_myfreqs2);
    w_freqdist=freqdistance/sum(freqdistance);% normallestirilmis fark

    if strcmp(oz.disttype,'JShannon')
        pi1=sum(myfreqs1)/(sum(myfreqs1)+sum(myfreqs2));
        pi2=sum(myfreqs2)/(sum(myfreqs1)+sum(myfreqs2));
        
        alldist(i) = jensen_Shannon_oz( myprobs1,myprobs2,pi1,pi2 );
    elseif strcmp(oz.disttype,'abs')
        alldist(i)=sum(sum(abs(myprobs1-myprobs2),2))/size(myprobs1,1);
    elseif strcmp(oz.disttype,'wabs')% may not be the best option
        alldist(i)=sum(sum(abs(myprobs1.*repmat(w_myfreqs1,1,size(myprobs1,1))-myprobs2*repmat(w_myfreqs2,1,size(myprobs2,1))),2));   
    elseif strcmp(oz.disttype,'euc')
        alldist(i)=sqrt(sum(sum((myprobs1-myprobs2).^2,2)))/size(myprobs1,1);
    elseif strcmp(oz.disttype,'weuc')
        alldist(i)=sqrt(sum(sum((myprobs1.*w_myfreqs1-myprobs2*w_myfreqs2).^2,2)));
        %freqdists(i)=sum(w_freqdist);       
    elseif   strcmp(oz.disttype,'bch') % bch not weighted
        alldist(i)=-log(sum(sum((myprobs1.*myprobs2).^0.5)));
    elseif   strcmp(oz.disttype,'hel') % https://en.wikipedia.org/wiki/Hellinger_distance
        alldist(i)=sqrt(1-sum(sum((myprobs1.*myprobs2).^0.5)));        
    else %euc
        alldist(i)=sum(sum((myprobs1-myprobs2).^2));
     end
end
%alldist=alldist/max(alldist);% akbilgic, added May 27, 2016
dist=sum(alldist);
% Copyright J. Andrew Howe & Oguz Akbilgic 2015
% this code may be freely used, modified, and distributed (at no charge)
% as long as this footer remains unaltered
% updated by Akbilgic 20161108




