% updated by EPS on 9-25-17
% first number is divdist = length of segments
% number of obs/divdist = number of segments
% length of series divide by divDist = number of segments
% sampling rate = second number   = divdist/sampling rate = number of seconds
% sampling rate/2 is the number of units for the spectra
% change Record number to segment number
function demo(preloadedData)

for setup = 1:1 %used to fold the setup code
    %set up figure and some variables
    
    f = figure('units','normalized','outerposition',[0 0 1 1],'Name','Plot Window','NumberTitle','off');
    
    %Initiate 
    timeSeries = NaN;
    structure = NaN;
    filename = '';
    pathname = '';
    map = NaN;
    PSD = 0;
    
    
    %Create Panels and TabGroups
    runPanel = uipanel(f,'units','normalized','position',[.3 .1 .69 .85]);
    gpPanel = uipanel(f,'units','normalized','position',[.3 .1 .69 .85]);
    batchPanel = uipanel(f,'units','normalized','position',[.3 .1 .69 .85]);
    invisiblerunPanel = uipanel(f,'units','normalized','position',[-0.1 0 1.2 1.1],'parent',runPanel);
    invisiblegpPanel = uipanel(f,'units','normalized','position',[-0.1 0 1.2 1.1],'parent',gpPanel);
    invisiblebatchPanel = uipanel(f,'units','normalized','position',[-0.1 0 1.2 1.1],'parent',batchPanel);
    runToolbarPanel = uipanel(f, 'units', 'normalized', 'position',[.01 .01 .28 .94],'BorderType','none');
    gpToolbarPanel = uipanel(f, 'units', 'normalized', 'position',[.01 .01 .28 .94],'BorderType','none');
    batchToolbarPanel = uipanel(f, 'units', 'normalized', 'position',[.01 .01 .28 .94],'BorderType','none');
    
    
    %aesthetic panels
    uipanel(f,'units','normalized','position',[.01 .47 .15 .33]); %parameters
    uipanel(f,'units','normalized','position',[.01 .1 .28 .3]); %data irregularities
    
    
    
    tabgp = uitabgroup(f,'SelectionChangedFcn',@switchTab);
    batchTab = uitab(tabgp,'Title','Single Run');
    %gpTab = uitab(tabgp,'Title','GP Model');
    runTab = uitab(tabgp,'Title','Single Channel');
    
    %makes single channel tab visible
    uistack(tabgp,'bottom');
    
    %Initiate all display boxes, etc for the two tabs
    
    %Text box is for labels
    paramsTextBox = uicontrol('Style','text','String','Parameters','units','normalized','position',[.02,.76,.05,.05]);
    irregTextBox = uicontrol('Style','text','String','Data Irregularities','units','normalized','position',[.02,.36,.07,.05]);
    
    channelTextBox = uicontrol('Style','text','String','Channel:','units','normalized','position',[.016 .54 .1 .05],'horizontalalignment','left','visible','off');
    %Edit box is for user inputs
    channelEditBox = uicontrol('Style', 'edit','units','normalized',...
        'Position', [.12 .57 .03 .02],'backgroundcolor',[1 1 1],'visible','off');
    
    critValTextBox = uicontrol('Style','text','String','Outlier Threshold:','units','normalized','position',[.016 .6 .1 .05],'horizontalalignment','left');
    
    critValEditBox = uicontrol('Style', 'edit','units','normalized',...
        'Position', [.12 .629 .03 .02],'backgroundcolor',[1 1 1]);
    
    rateTextBox = uicontrol('Style','text','String','Sampling Rate:','units','normalized','position',[.016,.66,.1,.05],'horizontalalignment','left');
    uistack(rateTextBox,'top')
    rateEditBox = uicontrol('Style', 'edit','units','normalized',...
        'Position', [.12 .688 .03 .02],'backgroundcolor',[1 1 1]);
    
    DivDistTextBox = uicontrol('Style','text','String','Segment Length:','units','normalized','position',[.016,.72,.1,.05], 'horizontalalignment','left');
    uistack(DivDistTextBox,'top')
    DivDistEditBox = uicontrol('Style', 'edit','units','normalized',...
        'Position', [.12 .745 .03 .02],'backgroundcolor',[1 1 1]);
    
    ClippingTextBox = uicontrol('Style','text','String','Suspected Clipping:','units','normalized','position',[.015,.33,.08,.03]);
    uistack(ClippingTextBox,'top')
    
    %fileName 
    
    %List box is for displying the flagged channels
    ClippingListBox = uicontrol('Style', 'listbox','units','normalized',...
        'Position', [.02 .14 .07 .19],'backgroundcolor',[1 1 1]);
    
    OverrangeTextBox = uicontrol('Style','text','String','Suspected Overrange:','units','normalized','position',[.1,.33,.1,.03]);
    uistack(OverrangeTextBox,'top')
    OverrangeListBox = uicontrol('Style', 'listbox','units','normalized',...
        'Position', [.115 .14 .07 .19],'backgroundcolor',[1 1 1]);
    
    ZeroTextBox = uicontrol('Style','text','String','Suspected to be off:','units','normalized','position',[.2051,.33,.08,.03]);
    uistack(ZeroTextBox,'top')
    ZeroListBox = uicontrol('Style', 'listbox','units','normalized',...
        'Position', [.21 .14 .07 .19],'backgroundcolor',[1 1 1]);
    
    %For displaying messages about the system (your run is complete, an error
    %occured, etc
    StatusListBox = uicontrol('Style', 'listbox','units','normalized',...
        'Position', [.3,.02,.6895,.1],'backgroundcolor',[1 1 1],'Max',2,'Value',[]);
    set(StatusListBox,'string','')
    updateStatusBar('Welcome! Please read in a file to start.');
    
    
    FileTextBox = uicontrol('Style','text','String','Current File:','units','normalized','position',[.01,.845,.05,.05]);
    uistack(FileTextBox,'top')
    
    FileListBox = uicontrol('Style','listbox','String','','units','normalized','position',[.06,.87,.192,.029],'Max',2,'value',[]);
    uistack(FileListBox,'top')
    structureTextBox = uicontrol('Style','text','String','Geometry File:','units','normalized','position',[.0019,.81,.06,.05]);
    uistack(structureTextBox,'top')
    
    structureListBox = uicontrol('Style','listbox','String','','units','normalized','position',[.06,.835,.192,.029],'Max',2,'value',[]);
    uistack(structureListBox,'top')
    %Buttons for calling the appropriate function
    readButton = uicontrol('Style', 'pushbutton', 'String', 'Load File','units','normalized',...
        'Position', [.24 .868 .045 .032],...
        'Callback', @readFileButton);
    structureButton = uicontrol('Style', 'pushbutton', 'String', 'Load File','units','normalized',...
        'Position', [.24 .834 .045 .032],...
        'Callback', @readStructureButton);
    plotMicButton = uicontrol('Style', 'pushbutton', 'String', 'Plot Channel','units','normalized',...
        'Position', [.61 .72 .35 .11],...
        'Callback', @plotMic,'parent',runToolbarPanel);
    batchPlotButton = uicontrol('Style', 'pushbutton', 'String', 'Plot Run Summary','units','normalized',...
        'Position', [.61 .72 .35 .11],...
        'Callback', @plotRun,'parent',batchToolbarPanel);
    gpPlotButton = uicontrol('Style', 'pushbutton', 'String', 'Plot GP','units','normalized',...
        'Position', [.61 .72 .35 .11],...
        'Callback', @plotGP,'parent',gpToolbarPanel);
   
    if nargin == 1
        timeSeries = preloadedData;
        updateStatusBar('Your data has loaded successfully');
        filename = 'Your data';
        
        [a, b] = size(timeSeries);
        numMics = min(a,b);
        if numMics ~= length(map)
            map = zeros(1,numMics);
            for mapi = 1:numMics
                map(mapi) = mapi;
            end
            
        end
        
        set(FileListBox,'String','Pre-Loaded Data')
    end
    set(structureButton,'visible','off');
    set(structureTextBox,'visible','off');
    set(structureListBox,'visible','off');
end




%Nested Functions for use in the application
    function readFileButton(~,~)
        %Prompts the user to select a .csv file from
        %their system, and loads the information into the matrix @timeSeries
        
        try
            
            [filename, pathname] = uigetfile({'*.*';'*.csv';'*.mat';'*.raw';'*.dat';'*.stat'});
            updateStatusBar('Loading your file')
            timeSeries = NaN;
            if (filename == 0)
                updateStatusBar('File not loaded')
                filename = '';
            elseif strcmp(filename(end-3:end),'.csv')
                timeSeries = csvread(strcat(pathname,filename));
                
            elseif strcmp(filename(end-3:end),'.mat')
                header = load(strcat(pathname,filename),'File_Header');
                
                mics = str2num(header.File_Header.NumberOfChannels); %%%%++
                samples = str2num(header.File_Header.NumberOfSamplesPerChannel);
                allData = zeros(samples,mics); %%%%%
                for i = 1:mics
                    ch = num2str(i);
                    % eval(strcat('Channel_',ch,'_Data',' = 0;'));
                    currentChannelData = load(strcat(pathname,filename),strcat('Channel_',ch,'_Data'));
                    names = fieldnames(currentChannelData);
                    allData (:,i) = currentChannelData.(names{1});
                end
                
                timeSeries = allData;
                %newFilename = strcat(file(1:end-4),'.csv');
                %csvwrite(newFilename,allData);
                
            elseif strcmp(filename(end-3:end),'.raw')
                allData = importdata(strcat(pathname,filename));
                dataCols = 0;
                for i = 1:length(allData.colheaders)
                    currentHeader = allData.colheaders(i);
                    headerString = currentHeader{1};
                    if length(headerString) < 4 || ~(strcmp(headerString(1:4),'Time') || strcmp(headerString(1:4),'Lase'))
                        dataCols = dataCols + 1;
                    end
                end
                dataLength = length(allData.data);
                timeSeries = zeros(dataLength, dataCols);
                nextCol = 1;
                for i = 1:length(allData.colheaders)
                    currentHeader = allData.colheaders(i);
                    headerString = currentHeader{1};
                    if length(headerString) < 4 || ~(strcmp(headerString(1:4),'Time') || strcmp(headerString(1:4),'Lase'))
                        timeSeries(:,nextCol) = allData.data(:,i);
                        nextCol = nextCol + 1;
                    end
                end
                
            elseif strcmp(filename(end-3:end),'.dat')
                fid = fopen(strcat(pathname,filename),'r');
                nch = str2double(inputdlg('Number of Channels'));
                %numRecords = str2double(inputdlg('Number of Records'));
                %lengthRecords = str2double(inputdlg('Length of Records'));
                timeSeries=fread(fid,[nch,3000000],'float32').'; %read mic voltage
                fclose(fid)
            elseif strcmp(filename(end-4:end),'.stat')
                allData = importdata(strcat(pathname,filename));
                timeSeries = allData.data.';
            elseif strcmp(filename(end-3:end),'.txt')
                k = 0;
                    head = 3;
                    numColumns = 10;
                    numRuns = 50; %impossibly high number
                    numChannels = 35;
                    formatSpec = '';
                    for k = 1:numColumns
                        formatSpec = strcat(formatSpec,'%f ');
                    end
                    
                    
                    for k = 1:numRuns
                        fid = fopen(strcat(pathname,filename));
                        cell{k} = textscan(fid,formatSpec,numChannels,'Headerlines',(k*numChannels+head-numChannels));
                        fclose(fid);
                        head = head + 4;
                    end

                    try
                    for k = 1:numRuns
                        for j = 1:numColumns
                            timeSeries((k-1)*numChannels+1:k*numChannels,j) = cell{1,k}{1,j};
                        end
                    end
                    catch
                        %do nothing
                    end

            end
            
            if ~strcmp(filename,'')
                set(invisiblerunPanel,'visible','on')
                set(invisiblegpPanel,'visible','on')
                set(invisiblebatchPanel,'visible','on')
                updateStatusBar(strcat(filename,' has been loaded.'))
            end
            %Note - assumes file is not corrupted  
            
            PSD = 0;
            [a, b] = size(timeSeries);
            numMics = min(a,b);
            if numMics ~= length(map)
                map = zeros(1,numMics);
                for mapi = 1:numMics
                    map(mapi) = mapi;
                end
                
            end
            
            set(FileListBox,'String',filename)
            %study = timeSeries;
            0;
        catch exception
            updateStatusBar('An error occurred. Please check your inputs for mistakes.');
            msgText = getReport(exception)
        end
    end

    function readStructureButton(~,~)
        
        try
            
            [structurefilename, structurepathname] = uigetfile({'*.csv'});
            updateStatusBar('Loading geometry file')
            if (structurefilename == 0)
                updateStatusBar('File not loaded')
                structurefilename = '';
                structure = NaN;
            else
                structure = csvread(strcat(structurepathname,structurefilename));
            end
            updateStatusBar(strcat(structurefilename,' is the current geometry.'))
            set(structureListBox,'String',structurefilename)
            updateStatusBar(strcat('GP plotted for ',filename));
        catch exception
            updateStatusBar('An error occurred. Please check your inputs for mistakes.');
            msgText = getReport(exception)
        end
    end


%
 

    function plotRun(~,~)
        %Plots some summary charts of a specific run across all microphones.
        %Checks each microphone for signs of corruption and flags microphones
        %with potential issues.
        
        
        try
            %Validation
            updateStatusBar('Plotting...')
            if isnan(timeSeries)
                warndlg('Please read in a file to start.');
                updateStatusBar('An error occurred.');
                return
            end
            if ~(validate(DivDistEditBox) & validate(rateEditBox))
                warndlg('Function input must be integer values.')
                updateStatusBar('An error occurred.');
                return
            end
            drawnow
            
            %initiate variables
            [dima,dimb] = size(timeSeries);
            mics = min(dima,dimb);
            divDist = str2num(get(DivDistEditBox,'String'));  % length of segment
            divDist
            num_segs = round(length(timeSeries)/divDist);  % number of segments
            %this needs to be generalized
            % length of segment * number of segments = length of series
            zerosMatrix = zeros(1,mics);
          
            samplingRate = str2num(get(rateEditBox,'string'));
            overlap = 0;
            % create hanning function with height 1/sqrt(0.375) starting
            % and ending at zero  - vector
            window=hanning(divDist,'periodic')/sqrt(0.375); 
            pcm=0;
            for n=1:num_segs
                pc=round(n/num_segs*100);
                if pc>pcm
                    %set(h.text8,'String',strcat(num2str(pc),'% done'));
                    pcm=pc;drawnow;
                end
            
                for m=1:mics
                    temp = timeSeries((1 + floor((n-1)*divDist*(1-overlap))):floor((n-1)*divDist*(1-overlap))+divDist,m);
                    r(:,m) = window.*temp;
              %old - r(:,m)=window.*eval(strcat('Channel_',num2str(m),'_Data(',num2str(floor((n-1)*divDist*(1-overlap))+1),':',num2str(floor((n-1)*divDist*(1-overlap))+divDist),')'));
                end
                % r will be a matrix of size divdist x number of mics
                s=fft(r);
            %reshape the fft into a matrix that is divdist by mics by mics
            % repmat repeats the rows of the fft - s- mics times
                t2=reshape(repmat(s,mics,1),divDist,mics,mics);
                t1=conj(t2);         %returns complex conjugate
                t1=permute(t1,[1,3,2]); %rearrange the order of t1
               								%columns=conj(FFT),column no.=channel no,pages are all copies 	   													%columns=FFT,page=channel no, different columns=are copies
               	if n==1
                  	t=t1.*t2;
                else
                  	t=t+t1.*t2;				%column=avg spec,column no.=1st(conj)chan,page=2nd chan.
                end
            end
            
             spectra=t/divDist/num_segs; 
       %    spectra=2*t/divDist/samplingRate/num_segs;       %normalize spectra
            PSD = zeros(divDist,mics);
            for index = 1:mics
                PSD(:,index) = spectra(:,index,index);
            end
                   PSD = PSD(1:length(PSD)/2,:);
                    %Functions to populate the plots / find errors
                    %batchPSDs(1);
                    subplot(2,2,1,'Parent',batchPanel);
                    semilogy(PSD)
                        title('Power Spectral Density of Each Channel (semilog)','fontsize',13)
                        ylabel('Power')
                        %EPS dropped kHz Hertz
                        xlabel('Frequency')
                        set(gca,'fontsize',11)
                    subplot(2,2,2,'Parent',batchPanel);
                    loglog(PSD)
                        title('Power Spectral Density of Each Channel (log-log)','fontsize',13)
                        ylabel('Power')
                        %EPS dropped kHz Hertz 
                         xlabel('Frequency')
                         set(gca,'fontsize',11)
            
            batchHeatmap(3);
            try
            populateErrorBoxes;
            catch exception
                updateStatusBar('Error in the data irregularity analysis');
                getReport(exception)
            end
            set(invisiblebatchPanel,'visible','off');
            
            %clear the runPanel, might be unnecessary?
            set(channelEditBox, 'backgroundcolor', [1 1 1]);
            
            %tell user that the function has terminated successfully
            updateStatusBar(strcat(filename,' has been plotted successfully'))
        catch exception
            updateStatusBar('An error occurred. Please check your inputs for mistakes.');
            msgText = getReport(exception)
        end
        
        
        
        %Nested Functions in plotRun to be called above
        function batchPSDs(plotNum)
            %Plots the average spectrum of each microphone. Creates two plots, the
            %first with a log y-axis and the second with a linear y-axis
            
            %@plotNum - Which of the four spaces to plot in. The second plot will be in
            %the next space and will overwrite anything already there.
            
            test = length(pwelch(timeSeries(1:divDist,1))) ;
            Pxx = zeros(test,num_segs);  % number of tests and segments
            
            averagePxx = zeros(test,mics);      %create zero matrix based on number of tests
            for ch = 1:mics
                for i = 1:num_segs
                    Pxx(:,i) = pwelch(timeSeries((1+(i-1)*divDist):(i*divDist),ch));
                end
                for i = 1:test                    %average over the segments
                    vector = Pxx(i,:);
                    averagePxx(i,ch) = mean(vector);
                end
                
            end
            textbox = uicontrol('Style','text','string','','units','normalized','position',[.5,.5,.045,.05],'Parent',batchPanel); %create display textbox
            legends = strsplit(num2str(1:mics)); %set legend array
            subplot(2,2,plotNum,'Parent',batchPanel)
            plot(averagePxx(:,:),'buttondownfcn',{@callbackAvgPowerSpectrum,textbox})   % average plot
            legend(legends)
            legend('off')
            title('Power Spectral Density of Each Channel (semilog)')
            ylabel('Power')
            xlabel('Bins')
            nextPlotNum = mod(plotNum+1,4); %wrap around next plot. Written specifically for a 4-plot setup
            if nextPlotNum == 0
                nextPlotNum = 4;
            end
            subplot(2,2,nextPlotNum,'Parent',batchPanel)
            semilogy(averagePxx(:,:), 'Tag', 'Pxx','buttondownfcn',{@callbackAvgPowerSpectrum,textbox});  % plot of log
            legend(legends)
            legend('off')
            title('Power Spectral Density of Each Channel (semilog)')
            ylabel('Power')
            xlabel('Bins')
            
        end
        
        function batchHeatmap(plotNum)
            %Counts the number of "violating" segments for each microphone. A violation
            %is a microphone that is behaving differently enough from the others to be
            %flagged, as defined below. Produces two plots: a bar chart counting how
            %many violating segments there are for each microphone, and a heatmap
            %to visualize these violations.
            
            %@plotNum - Which of the four spaces to plot in. The second plot will be in
            %the next space and will overwrite anything already there.
            
            num = round(length(timeSeries)/divDist);
            critval = str2num(get(critValEditBox,'string'));
            iqr_matrix = zeros(num,mics);  %create matrix to store iqrs
            n_violations = zeros(mics,1);
            iqr1 = ones(num,mics);   % specify indicator for violation type ir any
            for channel=1:mics
                
                countZeros = 0;
                
                means = zeros(num,1);  %zero vectors
                SD = zeros(num,1);
                fit = zeros(num,1);
                u95 = zeros(num,1);
                u05 = zeros(num,1);
                iqr95 = zeros(num,1);
                med = zeros(num,1);
                violations = zeros(num,1);
                
                for i = 1:num
                    means(i,1) = mean(timeSeries((1+(i-1)*divDist):(i*divDist),channel));
                    SD(i,1) = std(timeSeries((1+(i-1)*divDist):(i*divDist),channel));
                    u95(i,1) = prctile(timeSeries((1+(i-1)*divDist):(i*divDist),channel),0.975);
                    u05(i,1) = prctile(timeSeries((1+(i-1)*divDist):(i*divDist),channel),0.025);
                    iqr95(i,1) = u95(i,1)-u05(i,1);
                    med(i,1) = median(timeSeries((1+(i-1)*divDist):(i*divDist),channel));
                    iqr_matrix(i,channel) = iqr95(i,1); % increment the iqr matrix
                end
                
                %flag any standard deviations that are excessivley large/small
                meanSD    = mean(SD);
                SDSD      = std(SD);
                umean     = mean(iqr95);
                usd       = std(iqr95);
                medianIQR = median(iqr95);
                psdIQR    = iqr(iqr95)/1.35;  % pseudo std of the ranges IQRs
                
                f = 1;
                flagged = zeros(1,2);
                flagged(1,1) = -1;
                % here we create some values for the heatmap - we used integers
                % that are equally spaced to get values for the color bar
                % 4 values are specified 1=good, 3=exceed, 2=zero, 0=-exceed.
                for i = 1:num
                    %  if or(SD(i,1) > meanSD + 2.5 *SDSD, SD(i,1) < meanSD - 2.5 * SDSD)
                    % if or(iqr95(i,1) > umean + 3 *usd, iqr95(i,1) < umean - 3 * usd)
                    %   if or(iqr95(i,1) > medianIQR + 3 *psdIQR, iqr95(i,1) < medianIQR - 3 * psdIQR, iqr95(i,1)<0.001 )
                    if (iqr95(i,1) > medianIQR + critval *psdIQR| iqr95(i,1) < medianIQR - critval* psdIQR| iqr95(i,1)<0.001 )
                        %eps added
                        iqr1(i,channel)=1;
                        flagged(f,1) = i;
                        %     flagged(f,2) = SD(i,1);
                        flagged(f,2) = iqr95(i,1);
                        f = f+1;
                        % create indices for heatmap.  May need to change if standard changes
                        if iqr95(i,1) > medianIQR + critval *psdIQR
                            iqr1(i,channel) = 3 ;  % used in identifying heatmap values
                        elseif iqr95(i,1) < medianIQR - critval *psdIQR
                            iqr1(i,channel) = 0;  % used in heatmap for color
                        elseif abs(iqr95(i,1))< 0.001
                            iqr1(i,channel)= 2; % used in heatmap for color
                            countZeros = countZeros + 1;
                            
                        end
                        
                    end
                end
               % iqr1(1,1) = 0
                n_violations(channel,1) = length(flagged);  % count number of exceedences
                if flagged(1,1) == -1   n_violations(channel,1)=0;
                end
                if countZeros > .1*num
                    zerosMatrix(1,channel) = 1;
                end
            end
            
             %           iqr1(1:10,1:10)
            
            subplot(2,2,plotNum,'Parent',batchPanel)
            n_violations;
            
            title('Number of Segments Exceeding Standards')
            h = bar(map, n_violations,'k');
            %old code - tx = text(map,n_violations,num2str(n_violations),'FontSize',8);
            %           set(tx, 'Rotation', 90)
            P = title('Number of Segments Exceeding Standards','fontsize',12,'fontweight','bold');
            ylabel('Number of Violations','fontsize',10)
            xlabel('Channel','fontsize',10)
            %old code --  ylim([0 600]);
            %             pos = get(P,'position');
            %             set(P,'Position',[pos(1) pos(2)+20 pos(3)])
            set(gca,'fontsize',11)
      
            %----------------------------- Heatmap display
            %--------------------------------
            % inputs
            % map = list of channels
            % iqr1 = indicator matrix that displays 4 values (+,-exceed or not, zero)
            nextPlotNum = plotNum + 1;
            if nextPlotNum == 0
                nextPlotNum = 4;
            end
            subplot(2,2,nextPlotNum,'Parent',batchPanel)
           % old code used red, green, blue 
           %           cmap = [0,0,1     % set four colors to map blue  -upper limit
           %         1,1,1     % white = 1 OK
           %         0,1,0     % green = 0 value iqr1=2
           %         1,0,0];   % red >= upper limit
            
             cmap = [0,0,1     % set four colors to map blue  -upper limit
                    1,1,1     % white = 1 OK
                    0,1,1     % cyan = 0 value iqr1=2
                    1,0,0];   % red >= upper limit
            colormap( cmap)   % specify color map
            clims = [0, 3];  
            imagesc(iqr1, clims)    % create image 
           % 
           
            colorbar('Ticks',[.4,1.1,1.8,2.5],...
                'TickLabels',{'-Exceed','Good','Zero','Exceed'},'fontsize',12) 
            
            
            title('Heatmap of Standardized IQR by Channel','fontsize',12,'fontweight','bold');
            ylabel('Segment Number','fontsize',10)
            xlabel('Channel','fontsize',10)
            set(gca,'fontsize',11)
        end
        
        function populateErrorBoxes
            %Create histograms for each microphone, and use this to identify clipping,
            %overrange, and
            
            %Search for bad histograms
            subplot(2,2,3,'Parent',runPanel); %put histograms in a plot that will be cleared
            
            clippingMatrix = zeros(1,mics);
            for channel = 1:mics
                h = histogram(timeSeries(:,channel), 'Normalization', 'pdf');
                [binValue, edges] = histcounts(timeSeries(:,channel), h.NumBins);
                NumBins = h.NumBins;
                title('Distribution of Voltages')
                xlabel('Volts')
                ylabel('Probability Density')
                clippingValue = 0;
                histx = 1:NumBins;
                histxx = 1:.2:NumBins;
                histy = h.Values;
                try
                    histyy = spline(histx,histy,histxx);
                catch
                    histyy = histy;
                end
                plot(histyy);
                diffx = diff(histyy);
                plot(diffx);
                difflen = length(diffx);
                lowthresh = abs(quantile(diffx,.03));
                highthresh = abs(quantile(diffx,.97));
                thresh = min([lowthresh highthresh]);
                
                
                %make sure the histogram shows a continuous probability distribution
                continuous = 1;
 
                
                %Check for clipping on the edges of the graph
                if continuous == 1;
                    
                    
                    
                    %Check first third of graph for abnormally large bins
                    for i = 1:round(difflen/4)
                        if abs(diffx(i)) > thresh
                            clippingValue = 1;
                        end
                    end
                    
                    %Check last third of graph for abnormally large bins
                    for i = round(3*difflen/4):difflen
                        if abs(diffx(i)) > thresh
                            clippingValue = 1;
                        end
                    end
                    %Check for clipping in the last bin
                end
                
                if(clippingValue == 1)
                    clippingMatrix(1,channel) = 1;
                end
            end
            
            
            
            
            
            %update ZeroListBox
            
            set(ZeroListBox,'string','')
            
            for i = 1:mics
                if zerosMatrix(1,i) == 1 %if the mic was flagged as a zero
                    zeroChannel = map(i);
                    data = get(ZeroListBox,'String');
                    data{end+1} = ['Channel' ' ' num2str(zeroChannel)];
                    set(ZeroListBox,'String',data); %add to listbox
                    
                    if clippingMatrix(1,i) == 1
                        clippingMatrix(1,i) = 0;
                    end
                end
            end
            
            %update ClippingListBox
            
            set(ClippingListBox,'string','')
            set(OverrangeListBox,'string','')
            
            for i = 1:mics
                if(clippingMatrix(1,i) == 1) %if clipping
                    clippingChannel = map(i);
                    adVector = zeros(500,1);
                    for adV = 1:500
                        adVector(adV,1) = timeSeries(randi([1,length(timeSeries)],1,1),i);
                    end
                    [h,p] = adtest(adVector);
                    if 0 == 0 %EDIT CHANGE BACK h == 0 CHANGE TO [h p adstat cv] x is run number, y is value. Make new run tab for adtest, gp, dirichlet, etc
                        data = get(ClippingListBox,'String');
                        data{end+1} = ['Channel' ' ' num2str(clippingChannel)];
                        set(ClippingListBox, 'String', data); %add to listbox
                    else
                        data = get(OverrangeListBox,'String');
                        data{end+1} = ['Channel' ' ' num2str(clippingChannel)];
                        set(OverrangeListBox, 'String', data); %add to listbox
                    end
                    p
                end
            end
        end
        
    end

    function plotGP(~,~)
        try
            
            %Validation
            updateStatusBar('Plotting...')
            if isnan(timeSeries)
                warndlg('Please read in a file to start.');
                updateStatusBar('An error occurred.');
                return
            end
            if ~(validate(DivDistEditBox) & validate(rateEditBox))
                warndlg('Function input must be integer values.')
                updateStatusBar('An error occurred.');
                return
            end
            drawnow
            
            % The voltage rmse serves as the 'response' for the GP model. This computes
            % the rmse for each channel. I'm also removing channels 118 to 128 because
            % those were empty.
            
            
            rmse = std((timeSeries(1:8000,1:length(structure))));
            
            % Now we fit the GP
            
            gprMdl = fitrgp(structure,rmse);
            
            [ypred, ysd, yint]= resubPredict(gprMdl);
            
            subplot(2,2,1,'parent',gpPanel);
            plot(rmse,'bx')
            hold on
            plot(ypred,'ro')
            xlabel('Channel No.');ylabel('RMSe')
            legend('measured','GP prediction','location','best')
            
            stdResidual= zscore(rmse' - ypred);
            
            subplot(2,2,2,'parent',gpPanel);
            plot(1:length(structure),stdResidual,'kx')
            xlabel('Channel No.');
            ylabel('Residuals')
            
            %find if 2-d or 3-d
            [geoma,geomb] = size(structure);
            dimension = min(geoma,geomb);
            
            
            numFlagged = 0;
            for i = 1:length(structure)
                if abs(stdResidual(i)) > str2num(get(critValEditBox,'string'));
                    numFlagged = numFlagged + 1;
                    flaggedGP(numFlagged,1:dimension) = structure(i,1:dimension);
                end
            end
            subplot(2,2,3,'parent',gpPanel);
            if dimension == 2
                plot(structure(:,1),structure(:,2),'bo');
                hold on;
                plot(flaggedGP(:,1),flaggedGP(:,2),'ro');
            else
                scatter3(structure(:,1),structure(:,2),structure(:,3),'bo');
                hold on;
                scatter3(flaggedGP(:,1),flaggedGP(:,2),flaggedGP(:,3),'ro');
            end
            xlabel('x Coordinates')
            ylabel('y Coordinates')
            zlabel('z Coordinates')
            set(invisiblegpPanel,'visible','off');
        catch exception
            updateStatusBar('An error occurred. Please check your inputs for mistakes.');
            msgText = getReport(exception)
        end
    end

    function callbackAvgPowerSpectrum(line, ~, textbox)
        %Retrieves the channel number of the clicked on power spectrum
        %@line - the line whose channel number we want to retrieve
        %@textbox - the textbox where the channel number will be displayed
        set(textbox,'string',line.DisplayName)
    end

    function switchTab(~,tab)
        %Change which panel is showing when a tab is selected.
        %If another tab is added this will need to be expanded.
        selectedTab = tab.NewValue.Title;
        if strcmp(selectedTab,'Single Channel') == 1
            set(runPanel,'visible','on');
            set(runToolbarPanel,'visible','on');
            set(batchPanel,'visible','off');
            set(batchToolbarPanel,'visible','off');
            set(gpPanel,'visible','off');
            set(gpToolbarPanel,'visible','off');
            set(channelTextBox,'visible','on');
            set(channelEditBox,'visible','on');
            set(structureButton,'visible','off');
            set(structureTextBox,'visible','off');
            set(structureListBox,'visible','off');
            % set(critValTextBox,'visible','on');
            %set(critValEditBox,'visible','on');
        elseif strcmp(selectedTab,'GP Model') == 1
            set(runPanel,'visible','off');
            set(runToolbarPanel,'visible','off');
            set(gpPanel,'visible','on');
            set(gpToolbarPanel,'visible','on');
            set(batchPanel,'visible','off');
            set(batchToolbarPanel,'visible','off');
            set(channelTextBox,'visible','off');
            set(channelEditBox,'visible','off');
            set(structureButton,'visible','on');
            set(structureTextBox,'visible','on');
            set(structureListBox,'visible','on');
            %set(critValTextBox,'visible','off');
            %set(critValEditBox,'visible','off');
        else
            set(runPanel,'visible','off');
            set(runToolbarPanel,'visible','off');
            set(gpPanel,'visible','off');
            set(gpToolbarPanel,'visible','off');
            set(batchPanel,'visible','on');
            set(batchToolbarPanel,'visible','on');
            set(channelTextBox,'visible','off');
            set(channelEditBox,'visible','off');
            set(structureButton,'visible','off');
            set(structureTextBox,'visible','off');
            set(structureListBox,'visible','off');
            %set(critValTextBox,'visible','off');
            %set(critValEditBox,'visible','off');
            
        end
    end

    function plotMic(~,~)
        
        try
            
            updateStatusBar('Plotting...')
            
            %Validation
            if isnan(timeSeries)
                warndlg('Please read in a file to start.');
                updateStatusBar('An error occurred.');
                return
            end
            if ~(validate(DivDistEditBox) & validate(rateEditBox) & validate(rateEditBox))
                warndlg('Function input must be integer values.')
                updateStatusBar('An error occurred.');
                return
            end
            
            
            channelstring = get(channelEditBox,'String');
            ch = str2num(channelstring);
            channel = find(map == ch);
            critval = str2num(get(critValEditBox,'string')); %limit on standard deviation deemed acceptable. Arbitrary currently.
            if isempty(channel) %make sure the ch corresponds with a channel on the map
                updateStatusBar('Invalid Channel');
                set(channelEditBox,'BackgroundColor',[1 .78 .80]); %highlight red
                return
            else
                set(channelEditBox,'BackgroundColor',[1 1 1]); %return color to normal grey
            end
            drawnow;
            divDist = str2num(get(DivDistEditBox,'String'));
            num = round(length(timeSeries)/divDist);
            
            %call functions for plotting here
            runTimeSeries(1);
            runIQRMedian(2);
            runHistogram(3);
            %runPSD(4);
            subplot(2,2,4,'parent',runPanel);
            if PSD ~= 0
                testy = PSD;
                semilogy(PSD(:,channel));
                title(strcat({'Power Spectrum of Channel (semilog)'},{' '},int2str(channel)));
                %EPS dropped kHz Hertz
                xlabel('Frequency')
                
                ylabel('Power')
                % use to drop x axis ticks
                % use below to add x labels  
                 set(gca,'fontsize',11)
            end
            set(invisiblerunPanel,'visible','off');
            
            %Let the user know the action has successfully completed
            tempString = ['Channel',' ', num2str(ch), ' has been plotted'];
            updateStatusBar(tempString);
        catch exception
            updateStatusBar('An error occurred. Please check your inputs for mistakes.');
            msgText = getReport(exception)
        end
        
        
        %Nested functions to be called above
        function runTimeSeries(plotNum)
            %plot the time series data
            %@plotNum
            subplot(2,2,plotNum,'Parent',runPanel)
            samplingRate = str2num(get(rateEditBox,'String'));
            plot1 = plot(linspace(0,length(timeSeries(:,channel))/samplingRate,length(timeSeries(:,channel))),timeSeries(:,channel));
            plot1xlim =[0 divDist/samplingRate];
            xlim(plot1xlim);
            plot1ylim = [-12 12]; %GENERALIZE THIS
            ylim(plot1ylim);
            title(strcat({'Time Series of Channel'},{' '},int2str(ch)))
            xlabel('Time')
            ylabel('Volts')
            set(gca,'fontsize',11)
        end
        
        function runHistogram(plotNum)
            %make a histogram
            %num is the number of channels of lenght divDist
            subplot(2,2,plotNum,'Parent',runPanel)
            histogram(timeSeries(:,channel), 'Normalization', 'pdf');
            title('Distribution of Voltages')
            xlabel('Volts')
            ylabel('Probability Distribution')
            set(gca,'fontsize',11)
        end
        
       
        function runPSD(plotNum)
            %PSD plots
            subplot(2,2,plotNum,'Parent',runPanel)
            test = length(pwelch(timeSeries(1:divDist,channel)));
            Pxx = zeros(test,num);
            
            
            %find the average PSD across all records
            for i = 1:num
                Pxx(:,i) = pwelch(timeSeries((1+(i-1)*divDist):(i*divDist),channel));
            end
            averagePxx = zeros(test,1);
            for i = 1:test
                vector = Pxx(i,:);
                averagePxx(i,1) = mean(vector);
            end
            semilogy(averagePxx, 'Tag', 'Pxx');
            title('Average Power Spectral Density Across All Segments')
            ylabel('Power')
            xlabel('Bins')
        end
        
        function runIQRMedian(plotNum)
            %means and standard deviation, plus a chi squared test to determine the fit
            %of the section to a normal distrubution
            subplot(2,2,plotNum,'Parent',runPanel)
            means = zeros(num,1);  %zero vectors
            SD = zeros(num,1);
            fit = zeros(num,1);
            u95 = zeros(num,1);
            u05 = zeros(num,1);
            iqr95 = zeros(num,1);
            med = zeros(num,1);
            
            for i = 1:num
                means(i,1) = mean(timeSeries((1+(i-1)*divDist):(i*divDist),channel));
                SD(i,1) = std(timeSeries((1+(i-1)*divDist):(i*divDist),channel));
                u95(i,1) = prctile(timeSeries((1+(i-1)*divDist):(i*divDist),channel),0.975);
                u05(i,1) = prctile(timeSeries((1+(i-1)*divDist):(i*divDist),channel),0.025);
                iqr95(i,1) = u95(i,1)-u05(i,1);
                med(i,1) = median(timeSeries((1+(i-1)*divDist):(i*divDist),channel));
            end
            
            %flag any standard deviations that are excessivley large/small
            meanSD    = mean(SD);
            SDSD      = std(SD);
            umean     = mean(iqr95);
            usd       = std(iqr95);
            medianIQR = median(iqr95);
            psdIQR    = iqr(iqr95)/1.35;  % pseudo std of the ranges IQRs
            
            f = 1;
            flagged(1,1) = -1;
            for i = 1:num 
                % EPS added low iqr signal
                if (iqr95(i,1) > medianIQR + critval *psdIQR| iqr95(i,1) < medianIQR - critval * psdIQR| iqr95(i,1)<0.001 )
                    flagged(f,1) = i;
                    %     flagged(f,2) = SD(i,1);
                    flagged(f,2) = iqr95(i,1);
                    f = f+1;
                    
                end
            end
            n_violations = length(flagged);  % count number of exceedences
            if flagged(1,1) == -1   n_violations=0;
            end
            
            %chck = uicontrol('Style', 'checkbox', 'String', 'Compare Histogram?','units','normalized','Position', [.282 .06 .07 .03]);
            %set(chck,'Value',1);
            
            %plot flagged points
            x = 1:num;
            if flagged(1,1) == -1
                plot(x,iqr95,x,med,'ButtonDownFcn',@lineCallback)
                %   plot(x,means,x,SD,'ButtonDownFcn',{@lineCallback, timeSeries, divDist, channel, thisFigure,flagged})
            else
                plot(x,iqr95,x,med,flagged(:,1),flagged(:,2),'ko','ButtonDownFcn',@lineCallback)
                %   plot(x,means,x,SD,flagged(:,1),flagged(:,2),'ko','ButtonDownFcn',{@lineCallback, timeSeries, divDist, channel, thisFigure,averagePxx,chck})
                
            end
            title('IQR and Median')
            legend('Inter-quantile range','median','location','best')
            ylabel('Volts')
            xlabel('Segment Number')
            set(gca,'fontsize',11)
        end
        
        
        function lineCallback(~, hit) 
            
            % Show the selected point's time series data, spectral density, and
            % probability distribution
            
            % @hit - the point selected
            % @timeSeries - the matrix of data from one run
            % @divDist - length of each continuous record
            % @channel - the column containing the data for the appropriate channel
            % @thisFigure - the figure which is being looked at
            % @averagePxx - average of the spectral desities across each record
            % @chck - the histogram checkbox
            
            % get the x coordinate of the point clicked on
            xCoordinate = round(hit.IntersectionPoint(1,1));
            
            plotSize = length(timeSeries);
            targetValues = ((xCoordinate-1)*divDist+1):(xCoordinate*divDist);
            
            
            % plot the time series data and scroll to a convenient viewing point
            subplot(2,2,1)
            plot(1:plotSize,timeSeries(:,channel),targetValues,timeSeries(targetValues,channel))
            xlim([(xCoordinate-1)*divDist (xCoordinate)*divDist])
            
            % plot the spectral density of the selected region over the average
            subplot(2,2,4)
            semilogy(1:length(averagePxx),averagePxx,1:length(averagePxx),pwelch(timeSeries((1+(xCoordinate-1)*divDist):(xCoordinate*divDist),channel)),'Tag','Pxx')
            
        end
        
        
    end

    function expandPlot(~,~,p,pxlim,pylim)
        %Currently unused
        %Pulls up a new window with the plot p
        figure('units', 'normalized', 'position', [.005 0 .995 .94]);
        plot(p.YData)
        xlim([pxlim]);
        ylim([pylim]);
    end

    function updateStatusBar(message)
        %Update the status bar with the string @message
        
        data = get(StatusListBox,'String');
        if strcmp(data,'') %different syntax needed for first @message
            data{end+1} = message;
            data{end+1} = '';
        else
            data{end} = message;
            data{end+1} = ''; %adding a blank space makes it easier to read
        end
        set(StatusListBox,'String',data);
        index = size(get(StatusListBox,'string'),1);
        set(StatusListBox,'ListBoxTop',index);
        drawnow
    end

    function validated = validate(textbox)
        %check for integer values, change background to red if incorrect.
        %return a boolean @validated which is true if the input is valid.
        testString = get(textbox,'String');
        testArray = isstrprop(testString,'digit');
        baseline = ones(1,length(testString));
        if  isequal(baseline,testArray)
            set(textbox,'BackgroundColor',[1 1 1]);
            validated = true;
        else
            set(textbox,'BackgroundColor',[1 .78 .80]);
            validated = false;
        end
    end

    function setMap(~,~)
        
        try
            
            m = figure('units','normalized','outerposition',[.25 .25 .5 .5],'Name','Set Map','NumberTitle','off','menubar','none');
            newMap = map;
            len = length(newMap);
            winSize = 5;
            numPages = ceil(len/winSize);
            numLastPage = mod(len,winSize);
            page = 1;
            
            OneEditBox = uicontrol('Style', 'edit','units','normalized',...
                'Position', [.075 .5 .05 .05],'backgroundcolor',[1 1 1]);
            TwoEditBox = uicontrol('Style', 'edit','units','normalized',...
                'Position', [.275 .5 .05 .05],'backgroundcolor',[1 1 1]);
            ThreeEditBox = uicontrol('Style', 'edit','units','normalized',...
                'Position', [.475 .5 .05 .05],'backgroundcolor',[1 1 1]);
            FourEditBox = uicontrol('Style', 'edit','units','normalized',...
                'Position', [.675 .5 .05 .05],'backgroundcolor',[1 1 1]);
            FiveEditBox = uicontrol('Style', 'edit','units','normalized',...
                'Position', [.875 .5 .05 .05],'backgroundcolor',[1 1 1]);
            
            OneTextBox = uicontrol('Style', 'text','units','normalized',...
                'Position', [.075 .6 .05 .05]);
            TwoTextBox = uicontrol('Style', 'text','units','normalized',...
                'Position', [.275 .6 .05 .05]);
            ThreeTextBox = uicontrol('Style', 'text','units','normalized',...
                'Position', [.475 .6 .05 .05]);
            FourTextBox = uicontrol('Style', 'text','units','normalized',...
                'Position', [.675 .6 .05 .05]);
            FiveTextBox = uicontrol('Style', 'text','units','normalized',...
                'Position', [.875 .6 .05 .05]);
            
            leftButton = uicontrol('Style', 'pushbutton', 'String', '<<','units','normalized',...
                'Position', [.05 .08 .1 .1],...
                'Callback', @moveLeft);
            rightButton = uicontrol('Style', 'pushbutton', 'String', '>>','units','normalized',...
                'Position', [.85 .08 .1 .1],...
                'Callback', @moveRight);
            continueButton = uicontrol('Style', 'pushbutton', 'String', 'Continue','units','normalized',...
                'Position', [.45 .08 .1 .1],...
                'Callback', @replaceMap);
            loadPage();
        catch exception
            updateStatusBar('An error occurred. Please check your inputs for mistakes.');
            msgText = getReport(exception)
        end
        
        function updateNewMap()
            newMap(page*winSize-4) = str2double(get(OneEditBox,'string'));
            if page*winSize - 3 < len
                newMap(page*winSize-3) = str2double(get(TwoEditBox,'string'));
            end
            if page*winSize - 2 < len
                newMap(page*winSize-2) = str2double(get(ThreeEditBox,'string'));
            end
            if page*winSize - 1 < len
                newMap(page*winSize-1) = str2double(get(FourEditBox,'string'));
            end
            if page*winSize - 0 < len
                newMap(page*winSize-0) = str2double(get(FiveEditBox,'string'));
            end
        end
        function moveRight(~,~)
            updateNewMap();
            page = page + 1;
            if page > numPages
                page = numPages;
            else
                loadPage();
            end
        end
        function moveLeft(~,~)
            updateNewMap();
            page = page - 1;
            if page < 1
                page = 1;
            else
                loadPage();
            end
        end
        function loadPage()
            set(OneEditBox, 'visible', 'on');
            set(TwoEditBox, 'visible', 'on');
            set(ThreeEditBox, 'visible', 'on');
            set(FourEditBox, 'visible', 'on');
            set(FiveEditBox, 'visible', 'on');
            if page == numPages
                
                set(OneEditBox, 'string', int2str(newMap(winSize*page-4)));
                set(OneTextBox, 'string', int2str(winSize*page-4));
                if numLastPage > 1
                    set(TwoEditBox, 'string', int2str(newMap(winSize*page-3)));
                    set(TwoTextBox, 'string', int2str(winSize*page-3));
                    if numLastPage > 2
                        set(ThreeEditBox, 'string', int2str(newMap(winSize*page-2)));
                        set(ThreeTextBox, 'string', int2str(winSize*page-2));
                        if numLastPage > 3
                            set(FourEditBox, 'string', int2str(newMap(winSize*page-1)));
                            set(FourTextBox, 'string', int2str(winSize*page-1));
                            if numLastPage > 4
                                set(FiveEditBox, 'string', int2str(newMap(winSize*page-0)));
                                set(FiveTextBox, 'string', int2str(winSize*page-0));
                            end
                        end
                    end
                end
                
                if numLastPage < 2
                    set(TwoEditBox, 'visible', 'off');
                    set(TwoTextBox, 'string', '');
                end
                if numLastPage < 3
                    set(ThreeEditBox, 'visible', 'off');
                    set(ThreeTextBox, 'string', '');
                end
                if numLastPage < 4
                    set(FourEditBox, 'visible', 'off');
                    set(FourTextBox, 'string', '');
                end
                if numLastPage < 5
                    set(FiveEditBox, 'visible', 'off');
                    set(FiveTextBox, 'string', '');
                end
                
            else
                set(OneEditBox, 'string', int2str(newMap(winSize*page-4)));
                set(TwoEditBox, 'string', int2str(newMap(winSize*page-3)));
                set(ThreeEditBox, 'string', int2str(newMap(winSize*page-2)));
                set(FourEditBox, 'string', int2str(newMap(winSize*page-1)));
                set(FiveEditBox, 'string', int2str(newMap(winSize*page-0)));
                set(OneTextBox, 'string', int2str(winSize*page-4));
                set(TwoTextBox, 'string', int2str(winSize*page-3));
                set(ThreeTextBox, 'string', int2str(winSize*page-2));
                set(FourTextBox, 'string', int2str(winSize*page-1));
                set(FiveTextBox, 'string', int2str(winSize*page-0));
            end
        end
        function replaceMap(~,~)
            updateNewMap();
            map = newMap;
            close(m);
        end
        
        
        
        
        
        
        
    end


end







