% EMD_CLUSTER.M
%
% G. Rilling, last update: May 2005
%
% computes EMD using distcomp job manager (Empirical Mode Decomposition) according to:
%
% N. E. Huang et al., "The empirical mode decomposition and the 
% Hilbert spectrum for non-linear and non stationary time series analysis",  
% Proc. Royal Soc. London A, Vol. 454, pp. 903-995, 1998
%
% with variations reported in:
%
% G. Rilling, P. Flandrin and P. Gon�alv�s
% "On Empirical Mode Decomposition and its algorithms",
% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
% NSIP-03, Grado (I), June 2003
%
% default stopping criterion for sifting : 
%   at each point : mean amplitude < threshold2*envelope amplitude
%   &
%   mean of boolean array ((mean amplitude)/(envelope amplitude) > threshold) < tolerance
%   &
%   |#zeros-#extrema|<=1
%
% IMPORTANT !!!
% Before using emd_cluster, matlabpool with labs must be started.
% Initialization of matlabpool must be used by script 'cluster_init' 
% where number of labs and global variables used by emd_cluster are
% defined. Cluster is freed by calling script 'cluster_stop'.
%
%
% inputs:	
%		- x(mandatory): analysed signal (line vector), mandatory
%		- opt (optional): struct object with (optional) fields:
%         - min_part_size: minimum size which will be calculated using job
%                          manager on one lab, specified in number of
%                          local maxima
%                          default: 200
%         - extra_size: Each part used for calculation is prolonged on it's
%                       left and right side to minimize error at junction of parts. 
%                       First part is only prolonged on it's right side and last part is only 
%                       prolonged on it's left side.
%                       If not specifed will be half of min_part_size.
%                       Must be less or equall to min_part_size 
%                       default: 100
%		      - stop: threshold, threshold2 and tolerance 
%			  	        for sifting stopping criterion 
%				          default: [0.05,0.5,0.05]
%	        - maxmodes: maximum number of imfs extracted 
%         - max_iter: maximum number of iterations per mode
%                     default: 2000
%		      - fix: disable the stopping criterion and do exactly
%				         the value of the field number of sifting steps for each mode
%         - display: if equals to 1, displays messages about calculating 
%                    imfs, if equals to 0 nothing is displayed
%                    default: 1
%         
% outputs: 
%		- imf: intrinsic mode functions (last line = residual)
%		- ort: index of orthogonality
%		- nbits: number of iterations for each mode
%
% calls: 
%       - find_imf: calculates one imf (must be located on shared disk so every lab
%                   from cluster can access it)
%       - extr: calculates local extrema and zeros
%		    - io: computes the index of orthogonality
%
%examples:
%
%>>x = rand(1,512);
%
%>>cluster_init
%>>imf = emd_cluster(x);
%
%>>options.min_part_size = 100;
%>>options.maxmodes = 3;
%>>[imf,ort,nbits] = emd_cluster(x,options);
%
%>>cluster_stop
%
%
%


function [imf, ort, nbits] = emd_cluster(y, varargin)

if (nargin < 1)
    error('Please specify input signal');
end


opt_fields = {'min_part_size','extra_size','stop','maxmodes','max_iter', 'fix','display'};

% default parameters
defstop = [0.05,0.5,0.05]; % default stopping criterion
defopt.min_part_size = 200; % minimum part size which is calculated on cluster
defopt.extra_size = 0; % number of maxima which are added to each part, will be set later
defopt.stop = defstop;
defopt.maxmodes = Inf; % maximum number of calculated imfs
defopt.max_iter = 2000; % maximum number of iterations per imf
defopt.fix = 0; % fixed number of iterations of each mode
defopt.display = 1; % if equals to 1 displays messages on screen
opt = defopt;

if(nargin == 2)
    inopt = varargin{1};
    names = fieldnames(inopt);
    for nom = names'
      if isempty(strmatch(char(nom), opt_fields))
	  	  error(['bad option field name: ',char(nom)])
      end
      eval(['opt.',char(nom),' = inopt.',char(nom),';'])
    end
end

sd = opt.stop(1);
sd2 = opt.stop(2);
tol = opt.stop(3);


if(opt.extra_size > opt.min_part_size)
    error('Extra_size must be less or equall than min_part_size');
end



% find extremes
m = y(:)';
r = m;
n_imf = 1;

global lab_var;
global lab_res;
global numlabs;

% find extremas
[indmin, indmax ] = extr(m);   
n_extr = length(indmin) + length(indmax);

% sifting loop
while(n_extr > 2 && n_imf <= opt.maxmodes)
    

% split signal to parts
l_end = length(indmax); % total number of maximums
m_end = length(m);  % end of signal

if(l_end/numlabs > opt.min_part_size)
    l_extr = floor(l_end/numlabs); % part size handled by every worker will be larger than opt.min_part_size
    extra_size = floor(l_extr/2);
else
   % distribute parts more equally
    if(l_end > opt.min_part_size)
      n = floor(l_end/opt.min_part_size);   
      l_extr = floor(l_end/n);
      extra_size = floor(l_extr/2);
    else
     % only one part
      l_extr = l_end;
    end
end

if(opt.extra_size > 0)
  extra_size = opt.extra_size;
end


num_parts = floor(l_end/l_extr);
l_part_prev = 0; % previous part length

if (num_parts == 1)
  extra_size = 0;
end
 

% display messages for each imf

  if(opt.display)
    disp(' ');
    disp(['imf number: ',int2str(n_imf)])
    disp(['number of extrema: ', int2str(n_extr)])
    disp(['number of maxima: ', int2str(l_end)])
    disp(['number of parts: ',int2str(num_parts)])
    if(extra_size > 0)
      disp(['approx. number of maxima in each part: ',int2str(l_extr), ' + extra: ', int2str(extra_size)])
    end
    disp('Calculating using cluster')
  end


% part loop
for i=1:num_parts
         if(i == 1)
          % first part
           if(i == num_parts) % and last part also
             l_part = m_end;
             index_part_right = m_end;
           else
              % not last part
              % add extra_size more maximums on right side
             index_part_right = indmax(l_extr + extra_size);
             l_part = indmax((i-1)*l_extr+l_extr);
           end
           index_part_left = 0;
         elseif(i == num_parts)
           % last part (not first also)
           l_part = m_end;
           % add extra_size more maximums on left side
           index_part_right = m_end;
           index_part_left = indmax((i-1)*l_extr - extra_size);
         else
           %part in the middle
           l_part = indmax((i-1)*l_extr+l_extr);
           
           % add extra_size more maximums on left side and right side
           index_part_right = indmax((i-1)*l_extr+l_extr + extra_size);
           index_part_left = indmax((i-1)*l_extr - extra_size);
         end
         % signal which is sent for calculaction
         m_part_calc = m(index_part_left+1:index_part_right);
         
         lab_var.sig = m_part_calc;
         lab_var.labindex = i;
         lab_var.sd = sd;
         lab_var.sd2 = sd2;
         lab_var.tol = tol;
         lab_var.real_part = 1;
         lab_var.num_parts = num_parts;
         lab_var.max_n_bit = opt.max_iter;
         lab_var.fix = opt.fix;
         lab_var.l_part = l_part;
         lab_var.l_part_prev = l_part_prev;
         lab_var.index_part_left = index_part_left;
         lab_var.index_part_right = index_part_right;
         
         pmode('client2lab', 'lab_var', i); 
         % save positions so every part can be extracted from results 
         l_parts(i) = l_part;
         l_parts_prev(i)= l_part_prev;
         ind_part_left(i) = index_part_left;
         ind_part_right(i)= index_part_right;
   
   l_part_prev = l_part;
end

if(num_parts < numlabs)
    for i=num_parts+1:numlabs
        lab_var.sig = [];
        lab_var.labindex = 0;
        lab_var.sd = 0;
        lab_var.sd2 = 0;
        lab_var.tol = 0;
        lab_var.l_part = 0;
        lab_var.l_part_prev = 0;
        lab_var.index_part_left = 0;
        lab_var.index_part_right = 0;
        lab_var.real_part = 0;
        lab_var.num_parts = num_parts;
        pmode('client2lab', 'lab_var', i);  
    end
end


% add path were calc_iteration is located
dctRunOnAll 'if exist(''lab_var''), addpath(''/share/data/tmp''); end';

% run one iteration on all

dctRunOnAll 'if exist(''lab_var''), [lab_var.sig, lab_res.imf, lab_res.nbits ] = find_imf(lab_var); end;'

imf_whole = [];

% extract parts from result and combine them to one imf
for i=1:num_parts
  pmode('lab2client', 'lab_res', i);
  imf1 = lab_res.imf(l_parts_prev(i) - ind_part_left(i)+1:l_parts(i) - ind_part_left(i));
  imf_whole = [imf_whole,imf1]; % combine parts
end                                     
 
dctRunOnAll 'clear lab_var'
dctRunOnAll 'clear lab_res'

%clear path
dctRunOnAll 'if exist(''lab_var''), rmpath(''/share/data/tmp''); end';

% imf reached stopping criterion, we can stop
if(opt.display)
   disp(['Iterations reached: ', int2str(lab_res.nbits)]);
   disp(['Imf number ',int2str(n_imf), ' done.'])
end
     
  imf(n_imf,:) = imf_whole;
  nbits(n_imf) = lab_res.nbits;
  m = r - imf_whole; 
  
  % next imf
  r = m;
  n_imf = n_imf+1; 
  
  % find extremes for next imf
  [indmin, indmax ] = extr(m);   
  n_extr = length(indmin) + length(indmax);
 
end

% add residual as last imf
 imf(n_imf,:) = m;

ort = io(y,imf);





