% peak_finder.m
% Function which determines the "significant" peaks/valleys for an input timecourse.
% Jason Kelly, August 2003.
%
% Description 
% This function excepts a timeseries and a sensitivity level and generates an array 
% which represents the peaks and valleys of the timeseries.
%
% INPUT: a timeseries and a sensitivity level
% Input format: (y_data, sensitivity)
%
% y_data:
% An array of values representing a timeseries.  (The actual timepoints of the timeseries
% are not required as they do not come into play in determining the signifiance of peaks in 
% current algorithm.
%
% sensitivity:
% Value bewteen 0-1 which represents the level of significance required for a peak to be 
% considered significant.  
% Sensitivity * (max value in time series - min value in timeseries) = minimum size for a significant peak.
% "size" of the peak is determined by the distance between the peak and the nearest 2 valleys.
%
% OUTPUT: array of length M where M equals the number of elements in the y_data array.
% Each position in the array holds a valut of -1, 0 , or 1 and corresponds to a timepoint
% in the y_data array.
% -1 = valley
%  0 = insiginificant 
%  1 = peak

function out=peak_finder(y_data, sensitivity)

% number of timepoints in the data series
num_timepoints = length(y_data);
% the size of the minimum "jump" from a valley to a peak based on the signiificance.
min_jump = (max(y_data)-min(y_data))*sensitivity;
% Initialization of the output array. 
peak_valley = zeros(num_timepoints,1);


% find all peaks (both insignificant and significant)

% Differentiate the data matrix at each point in the y_data array, 
diff_y_data = diff(y_data);

% set the last_sign var to be equal to the slope of the first point.
if diff_y_data(1)<0
    last_sign = -1;
else
    last_sign = 1;
end

% find the peaks in the timecourse wherever a slope changes from positive to negative (ignoring the endpoints for now)
for i=2:(num_timepoints-1)
    if diff_y_data(i)<0
        sign = -1;
    else
        sign = 1;
    end
    
    if (sign~=last_sign)
        %if going from positive to negative slope
       if(sign<last_sign) 
            peak_valley(i) = 1;
        end
        last_sign = sign;
    end
end

% determine if endpoints are peaks or valleys, the current algorithm considers endpoints peaks or 
% valleys regardless of whether they are significant.
if diff_y_data(1)>0 
    %if positive slope between first 2 points, 1st point is a valley
    peak_valley(1) = -1;
else
    %else its a peak
    peak_valley(1) = 1; 
end
 
if diff_y_data(num_timepoints-1)>0 
    %if positive slope between last 2 points, last point is a peak
    peak_valley(num_timepoints) = 1;
else
    %else its a valley
    peak_valley(num_timepoints) = -1; 
end 

% Create a combined matrix containing the data points along with their respecive peak/valley characteristic
% Sorted according to time points
combo_sort_by_TP(:,1) = y_data;
combo_sort_by_TP(:,2) = peak_valley; 

[sorted_timecourse,sort_TC_index] = sort(combo_sort_by_TP(:,1));

% Sort combined matrix by amplitude (time course)
combo_sort_by_TC = combo_sort_by_TP(sort_TC_index,:);
combo_sort_by_TC = combo_sort_by_TC((end:-1:1),:);
sort_TC_index = sort_TC_index((end:-1:1),:);

% Go through all time points
for i=1:(num_timepoints)
    % If the current time point is a peak.
    if combo_sort_by_TC(i,2) == 1
        j=1;
        left_sig=0; %1 = significant ; -1 = insignificant
        % while left side not significant and not the beginning or end points and whose left point isn't at time position 0 (out of bounds)
        while ((left_sig~=1)&(sort_TC_index(i)~=num_timepoints)&(sort_TC_index(i)~=1)&(sort_TC_index(i)-j)~=0 )
            
            %----------COMPARE TO THE LEFT-------------------------------
            position_in_time = (sort_TC_index(i)-j);
            current_peak = combo_sort_by_TP(sort_TC_index(i),:);
            compare_left = (combo_sort_by_TP(sort_TC_index(i)-j,:));
            
            % if the current peak is greater than the significant jump size from the compare point, it's significant.
            if ((current_peak(1) - compare_left(1))>min_jump)
                left_sig = 1;
            end
            
            % if the comparing point is a peak and is lower then the current peak, then the comparing point is not a significant peak.
            % these rules dont apply to the very first point.
            if ((compare_left(2) == 1)&((sort_TC_index(i)-j)~=1))
                if ((current_peak(1)-compare_left(1))>0)
                    combo_sort_by_TP(sort_TC_index(i)-j,2)=0;
                else %if the comparing peak is higher than the current peak, then the current peak is not a significant peak.
                    left_sig=-1;
                    combo_sort_by_TP(sort_TC_index(i),2)=0;
                end
            end
            
             % if it reaches the beginning of the time course and the beginning is a peak, then current peak is insignificant.
             % if beginning is a valley than current peak is significant.
             if ((sort_TC_index(i)-j)==1)
                 if (compare_left(2) == 1)
                     left_sig=-1;
                     combo_sort_by_TP(sort_TC_index(i),2)=0;
                 end
                 
                 if (compare_left(2) == -1)
                     left_sig = 1;
                 end
                 
                 if (compare_left(2) == 0)
                     'problem: shouldnt get here --  end peaks were not set properly. (Must be either a peak or valley)'  
                 end
             end
             
            j=j+1;
        end
        
        j=1;
        right_sig=0; %1 = significant ; -1 = insignificant
       
        %--------------COMPARE TO THE RIGHT-------------------------------
        while ((right_sig~=1)&(sort_TC_index(i)~=num_timepoints)&(sort_TC_index(i)~=1)&(sort_TC_index(i)+j)<=num_timepoints )
            position_in_time = (sort_TC_index(i)+j);
            current_peak = combo_sort_by_TP(sort_TC_index(i),:);
            compare_right = (combo_sort_by_TP(sort_TC_index(i)+j,:));
            
            % if the current peak is greater than the significant jump size from the compare point, its sugnificant.
            if ((current_peak(1) - compare_right(1))>min_jump)
                right_sig = 1;
            end
                       
            % if the comparing point is a peak and is lower then the current peak, then the comparing point is not a significant peak
            % these rules dont apply to the very last point.
            if ((compare_right(2) == 1)&((sort_TC_index(i)+j)~=num_timepoints))
                if ((current_peak(1)-compare_right(1))>0)
                    combo_sort_by_TP(sort_TC_index(i)+j,2)=0;
                else %if the comparing peak is higher than the current peak, then the current peak is not a significant peak
                    right_sig=-1;
                    combo_sort_by_TP(sort_TC_index(i),2)=0;
                end
            end
            
            % if it reaches the end of the time course and the end is a peak, then current peak is insig
            % if end is a valley than current peak is significant
            if ((sort_TC_index(i)+j)==num_timepoints)
                if (compare_right(2) == 1)
                    right_sig=-1;
                    combo_sort_by_TP(sort_TC_index(i),2)=0;
                end
                
                if (compare_right(2) == -1)
                    right_sig = 1;
                end
                
                if (compare_right(2) == 0)
                    'problem: shouldnt get here --  end peaks were not set properly. (Must be either a peak or valley)'  
                end
            end
            
            j=j+1;
        end   
        
        
        % Rebuild the combo sorted by time course to include the lost peaks (if peak was significant it wont be lost)
        [sorted_timecourse,sort_TC_index] = sort(combo_sort_by_TP(:,1));
        % sort by amplitude (time course)
        combo_sort_by_TC = combo_sort_by_TP(sort_TC_index,:);
        combo_sort_by_TC = combo_sort_by_TC((end:-1:1),:);
        sort_TC_index = sort_TC_index((end:-1:1),:);
    end  
end

%---------SET VALLEYS--------------------------------
% Valleys are set as the lowest point between two peaks.

no_first = 1; % there hasn't been a first peak yet
low_valley = max(combo_sort_by_TP(:,1)); % initialize the low valley as the highest point, so valley will be guaranteed to be lower.
for i=1:(num_timepoints)
    if (combo_sort_by_TP(i,2) == 1) 
        if (no_first == 1)
            no_first = 0;
            low_valley = max(combo_sort_by_TP(:,1));
        else
            % if its a peak (but not the first one) set the low valley for the previous section and reset the low valley
            combo_sort_by_TP(low_valley_position,2) = -1;
            low_valley = max(combo_sort_by_TP(:,1));
        end
    end
    if (combo_sort_by_TP(i,1)<low_valley)
        low_valley = combo_sort_by_TP(i,1);
        low_valley_position = i;
    end

end

% Out put the combined matrix sorted by timepoints.
out = combo_sort_by_TP(:,2);