bottom_up_segmentation.m

分段线性分割的算法

function segment = bottom_up_segmentation(agreement,data,num_segments,forceplot)

% Type the function name with no arguments for help.
if nargin == 0                          % The function was called with no arguments,...
    agreement = 'none';                 % ...set the "agreement" to none, and echo Conditions For Use and Help Information
end;
    
if ~strncmpi(agreement,'I agree',6)     % Explain Conditions for Use, and usage.
    
    disp(' ');
    disp('(c) Eamonn Keogh, 2003');
    disp(' ');
    disp('This code may only be used if you agree to the conditions below.');
    disp(' ');
    disp('  1) I am not responsible for any errors or bugs in the code.');
    disp(' ');
    disp('  2) You may not resell this code, or use it for any commercial application.');
    disp(' ');
    disp('  3) You may not remove the authors name, or any comments from this code.');
    disp(' ');
    disp('  4) If you use this code as part of any published research, you must reference one of the relevant papers listed below, and send a pointer to your paper to [e-mail]');
    disp(' ');
    disp('To make this code work, you must pass in a string as the first argument. The string must be ''I agree'', which affirms that you agree to the conditions above.  ');
    disp(' ');
    disp('EXAMPLE:  bottom_up_segmentation(''I agree'',data,20,1)');
    
    disp(' ');
    disp(' ');
    disp('This function approximates a time series by a sequence of piecewise linear segments. ');
    disp(' ');
    disp(' "data" should be a column vector of a time series.');
    disp(' "num_segments" should be an interger greater than or equal to 2, and less than length(data)/2.');
    disp(' "forceplot", is an optional argument, any value here will make a plot of the time series, together with the segmented representation appear.');
    disp(' EXAMPLE USAGE:');
    disp('      >> x = random_walk(1000,1);                       ');
    disp('      >> bottom_up_segmentation(''I agree'',x,15,1);    ');
    disp(' ');
    disp('Keogh, E., Chu, S., Hart, D. & Pazzani, M. (2001). An Online Algorithm for Segmenting Time Series. In Proceedings of IEEE International Conference on Data Mining. pp 289-296. ');
    disp(' ');
    disp(' ');
    disp(' This code is provided as a tool. It is NOT highly optimized for speed, nor is it the highly commented version. If either of these things are important to you, email me.');
     
    return;
end;


% The basic algorithm works by creating a fine segmented representation then merging the lowest cost segments until only 'num_segments' remain.

left_x       = [1 : 2 : size(data,1)-1];    % Find the left x values vector for the "fine segmented representation".
right_x      = left_x + 1;                  % Find the right x values vector for the "fine segmented representation".
right_x(end) = size(data,1);                % Special case, the rightmost endpoint.
number_of_segments = length(left_x );       % Calculate the number of segments in the initial "fine segmented representation".


% Initialize the segments in the "fine segmented representation". 

for i = 1 : number_of_segments 
   
   segment(i).lx = left_x(i);
   segment(i).rx = right_x(i);
   segment(i).mc = inf;
end;

% Initialize the merge cost of the segments in the "fine segmented representation". 

for i = 1 : number_of_segments - 1
   coef = polyfit([segment(i).lx :segment(i+1).rx]',data(segment(i).lx :segment(i+1).rx),1);
   best = (coef(1)*( [segment(i).lx :segment(i+1).rx]' ))+coef(2);
   segment(i).mc = sum((data([segment(i).lx :segment(i+1).rx]')-best).^2);
end;


% Keep merging the lowest cost segments until only 'num_segments' remain. 

while  length(segment) > num_segments  
   
   
   [value, i ] = min([segment(:).mc]);                              % Find the location "i", of the cheapest merge.
      
  
   if i > 1 & i < length(segment) -1								% The typical case, neither of the two segments to be merged are end segments
      
    	coef = polyfit([segment(i).lx :segment(i+2).rx]',data(segment(i).lx :segment(i+2).rx),1);
    	best = (coef(1)*( [segment(i).lx :segment(i+2).rx]' ))+coef(2);
       
        segment(i).mc = sum((data([segment(i).lx :segment(i+2).rx]')-best).^2);

	 	segment(i).rx = segment(i+1).rx;
    	segment(i+1) = [];
    
    	i = i - 1;
    
        coef = polyfit([segment(i).lx :segment(i+1).rx]',data(segment(i).lx :segment(i+1).rx),1);
    	best = (coef(1)*( [segment(i).lx :segment(i+1).rx]' ))+coef(2);
    
        segment(i).mc = sum((data([segment(i).lx :segment(i+1).rx]')-best).^2);
       
   elseif i == 1                                                    % Special case: The leftmost segment must be merged.
       
	    coef = polyfit([segment(i).lx :segment(i+2).rx]',data(segment(i).lx :segment(i+2).rx),1);
    	best = (coef(1)*( [segment(i).lx :segment(i+2).rx]' ))+coef(2);
    	segment(i).mc = sum((data([segment(i).lx :segment(i+2).rx]')-best).^2);

	 	segment(i).rx = segment(i+1).rx;
        segment(i+1) = [];
              
   else                                                             % Special case: The rightmost segment must be merged.
     
       segment(i).rx = segment(i+1).rx;
       segment(i).mc = inf;
       segment(i+1) = [];
    
       i = i - 1;
       
       coef = polyfit([segment(i).lx :segment(i+1).rx]',data(segment(i).lx :segment(i+1).rx),1);
       best = (coef(1)*( [segment(i).lx :segment(i+1).rx]' ))+coef(2);
    
       segment(i).mc = sum((data([segment(i).lx :segment(i+1).rx]')-best).^2);

   end;
           
end;


%----------------------------------------------
residuals=[];

for i = 1 : length(segment) 
   
        coef = polyfit([segment(i).lx :segment(i).rx]',data(segment(i).lx :segment(i).rx),1);
    	best = (coef(1)*( [segment(i).lx :segment(i).rx]' ))+coef(2);
        residuals = [ residuals ; sum((data([segment(i).lx :segment(i).rx]')-best).^2)];
end;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% If the user passed in an extra argument, then plot the original time series, together with the segments %%%%%%%%%%%%%


if nargin > 3
    hold on;
    plot(data,'r'); 
end; 

   
    temp = [];
    for i = 1 : length(segment) 
   
      coef = polyfit([segment(i).lx :segment(i).rx]',data(segment(i).lx :segment(i).rx),1);
      best = (coef(1)*( [segment(i).lx :segment(i).rx]' ))+coef(2);
      
      segment(i).ly = best(1);
      segment(i).ry = best(end);
      segment(i).mc = residuals(i);
      
      if nargin > 3
            plot([segment(i).lx :segment(i).rx]', best,'b');
      end;
  
      temp = [temp; [segment(i).ly segment(i).ry]];
   
    end;

    
if nargin > 3
    for i = 1 : length(segment)  - 1 
        plot([segment(i).rx :segment(i+1).lx]', [ temp(i,2) temp(i+1,1)  ],'g');
    end;
end;