%libraryLookup : takes a previously generated library (using
%                createTraceLibrary) as well as a test growth rate
%                histogram and finds the closest clusters in the library
%                based on this histogram.  It returns an average trace of
%                the cluster traces (uses a weighted average based on the
%                order of the lookup and aligns the traces by their
%                maximums before averaging), the standard deviation of said
%                trace, as well as the indices of the clusters used.
%
%   Inputs:
%       library : the structure which represents a set of clusters as
%                 generated by createTrace Library
%       GR      : the test growth rate histogram which is being used for
%                 the lookup.  This should be the same binning as used in
%                 createTraceLibrary.
%       plot_results : boolean whether to plot the average trace lookup
%   
%   Output:
%       lookup_trace: information regarding the average trace of the
%                     clusters used.  lookup_trace(1,:) is the average of
%                     the cluster traces which have been weighted by the order
%                     of the lookup and aligned by their maximum value
%                     before averaging.  lookup_trace(2,:) is the weighted
%                     standard deviation of this same process.
%       lookup_clusters : contains the indices of the lookup clusters in
%                         order of closest to farthest from GR
%           
% Nathan Willy, Kural Group, Ohio State University, willy.2@osu.edu
%

function [lookup_trace,lookup_clusters] = libraryLookup(library,GR,plot_result)

close all

    %number of clusters to find and average
    average_number=10;

    %make matrix containing all library histograms
    comp=[];
    %and find length of longest library trace
    longest_cluster=0;
    for(i=1:length(library))
        comp(i,1:length(library(i).GR13))=library(i).GR13;
        
        if(length(library(i).ave_tr) > longest_cluster)
            longest_cluster = length(library(i).ave_tr);
        end
    end

    %calculate distances between library histograms and the test histogram 
    d=hist_dist(GR,comp);


    %unaligned traces from the library
    tr=[];
    %traces from library aligned by maximum intensity
    temp_tr_aligned=[];
    %cluster numbers from the lookup
    lookup_inds=[];

    %perform librray lookup
    for(j=1:average_number)
        %find the minimum distance
        [~,ind]=min(d);
        %save index
        lookup_inds(end+1)=ind;
        
%         tr(j,1:length(library(ind).ave_tr)) = library(ind).ave_tr-min(library(ind).ave_tr);

        %find position of lookup trace's maximum intensity
        [~,mi] = max(library(ind).ave_tr);
        %save traces aligned to it's maxmimum
        temp_tr_aligned(j,longest_cluster + 1 - mi : longest_cluster + length(library(ind).ave_tr) - mi) = library(ind).ave_tr;
        %replace smallest distance with an impossibly large value
        d(ind) = 10^10;
    end

    %final aligned average trace
    aligned_ave=[];
    %standard deviation of final aligned trace
    SD=[];

    
    %find all non-zero elements in temp_aligned_tr
    [a,b]=find(temp_tr_aligned);
    numbers=[1:max(b)];
    %iterate through times where at least one trace has a non-zero element
    for(k=min(b):max(b))
        
        %collectin of all non-zero elements at this time point
        list=[];
        %weight coefficients used in the averaging
        weights=[];
        
        %iterate through all the traces
        for(j=1:average_number)

            %add to list the value of the cluster trace which is at the
            %current time point
            if(library(lookup_inds(j)).LT>0)
                temp=find(squeeze(temp_tr_aligned(j,:)));
                [~,temp_min]=min(abs(numbers(temp)-k));
                list(end+1)=temp_tr_aligned(j,temp(temp_min));
                %add value to the end of weight coefficient list
                weights(end+1)=average_number-j+1;
            end
        end
        %calculate weighted average value for the intensities as well as
        %the weights standard deviation
        if(length(list)>0)
            aligned_ave(end+1)=sum(list.*weights)/sum(weights);
            SD(end+1)=std(list,weights);
        end
    end

    %save averge trace and standard deviation to output
    lookup_trace(1,1:length(aligned_ave)) = aligned_ave;
    lookup_trace(2,1:length(SD)) = SD;

    %save lookup clusters' indices to output
    lookup_clusters(1,:) = lookup_inds;
    
    if(plot_result)
        figure
    %unused code for plotting the component traces aligned to their maxima
    %with a grayscale coloring indicating their relative weighting in the
    %average (black is highest weight light gray is lowest)
    
%     for(i=1:length(lookup_clusters))
%         [~,mi]=max(library(lookup_clusters(i)).ave_tr);
%         plot([1:length(library(lookup_clusters(i)).ave_tr)]-mi,library(lookup_clusters(i)).ave_tr,'Color',[0,0,0]+(i-1)/average_number)
%         hold on
%     end
%     [~,mi]=max(lookup_trace(1,:));
%     plot([1:length(lookup_trace(1,:))]-mi,lookup_trace(1,:),'b')
%     plot([1:length(lookup_trace(1,:))]-mi,lookup_trace(1,:)+lookup_trace(2,:),'r')
%     plot([1:length(lookup_trace(1,:))]-mi,lookup_trace(1,:)-lookup_trace(2,:),'r')
        hold on
        plot(lookup_trace(1,:),'b')
        plot(lookup_trace(1,:)+lookup_trace(2,:),'r')
        plot(lookup_trace(1,:)-lookup_trace(2,:),'r')
        xlabel('time (s)')
        ylabel('normalized inetnsity')
        title(strcat('weighted average of ',num2str(average_number),' traces (blue) with standard deviations (red)'))
    end
        
    
    
end