%   Copyright and terms of use (DO NOT REMOVE):
%   The code is made freely available for non-commercial uses only, provided that the copyright 
%   header in each file not be removed, and suitable citation (see below) be made for paper 
%   published based on the code.
%
%   The code is not optimized for speed, and we are not responsible for any errors that might
%   occur in the code.
%
%   The copyright of the code is retained by the authors.  By using this code you
%   agree to all the terms stated above.
%
%   Tapinos A, Mendes P (2013) A Method for Comparing Multivariate Time Series
%   with Different Dimensions. PLoS ONE 8(2): e54201. doi:10.1371/ journal.pone.0054201
%
%
%	SMETS function returns the SMETS value between two multivariate time
%	series MTSA and MTSB. Time series can have a different number of
%	components where (M>= N or N>=M).
%   The multivariate time series must be in an m by x format
%	where m represents the number of rows and x represent the number of
%	columns in the matrix. Row j is the j-th feature component of the multivariate
%	time series and column i is the i-th data point of the j-th feature component.
%
%   Input: 
%
%   MTSA:         A multvariate time series or a multivariate time series
%                 representation with M by X dimensions. M is the number of
%                 componets (rows) in the time series and  X is the number of
%                 data points (columns) in the multivariate time series
%
%   MTSB:         A multvariate time series or a multivariate time series
%                 representation with N by X dimensions. N is the number of
%                 componets (rows) in the time series and  X is the number of
%                 data points (columns) in the multivariate time series.
%
%   Output:
%
%   SMETSval: The semimetric distance between the two multivariate time
%   series MTSA and MTSB.
%
%
% Copyright (c) 2013, Avraam Tapinos, Pedro Mendes.  All rights reserved.


function [SMETSval]=SMETS(MTSA , MTSB)

if (size(MTSB,1) < size(MTSA,1) || size(MTSB,1) == size(MTSA,1))
    Y=MTSA; X=MTSB;
else
    Y=MTSB; X=MTSA;
end


%	Calculates the distance between each univariate component from
%	MTSA and all the univariate components of MTSB.
dist_table=zeros(size(X,1),size(Y ,1));
for a=1:size(Y,1)
    for b=1:size(X,1)
        dist_table(b,a)=sqrt(sum((Y(a,:)-X(b,:)).^2));
    end
end

%	Calculates the Entropy rate of all the components of the MTS with
%	the larges number of feature components.
temp_dist_table=dist_table; 
H=zeros(size(Y ,1),1);
for c=1:size(Y,1)
    histo = hist(Y(c,:), min(Y(c,:)):max(Y(c,:))+1);
    probabilities = histo/sum(histo);
    probabilities(isnan(probabilities))=0;
    log2Pi = -log2(probabilities);
    Pilog2Pi = probabilities.*log2Pi;
    H(c,1)=sum(Pilog2Pi);
end

H=H/sum(H);
d_matrix=zeros(size(dist_table,1),1);


%   Finds the smallest distance between the components from MTSA and MTSB
%   and records the distance to d_matrix. Furthermore, the components that 
%   shared the smallest distance are removed form the process  
for d=1:size(dist_table,1)
    [ymin, iymin]=min(temp_dist_table);
    [xmin, ixmin]=min(ymin);
    d_matrix(d,1)=xmin;
    if (d<size(dist_table,1))
        temp_dist_table(:,ixmin)=[];
        temp_dist_table(iymin(1,ixmin),:)=[];
        dist_table(:,ixmin)=[]; H(ixmin,:)=[];
    else
        H(iymin,:)=[];
        temp_dist_table(:,1)=[];
        dist_table(:,iymin)=[];
    end
end

%   Calculates the Lp norm value of the d_matrix
Lpnorm=nthroot(sum(d_matrix.^size(d_matrix,1)),size(d_matrix,1));

%   Calculates the penalty values
if (size(Y,1)> size(X,1))
    RE=zeros(size(Y ,1)-size(X,1),1);
    for e=1:(size(Y,1)-size(X,1))
        RE(e,1)=min(dist_table(:,e))*H(e,1);
    end
    EP=sum(RE);
    P=(size(Y ,1)-size(X,1))/(size(Y ,1)+size(X,1));
    SMETSval=sqrt(((Lpnorm+EP)^2)+(P^2));
else
    SMETSval=Lpnorm;
end