function [results]=cmb_outliers(results,odir,figdir)
%CMB_OUTLIERS    Outlier analysis of core-diffracted data
%
%    Usage:    results=cmb_outliers(results)
%              results=cmb_outliers(results,odir)
%              results=cmb_outliers(results,odir,figdir)
%
%    Description:
%     RESULTS=CMB_OUTLIERS(RESULTS) provides an interface for graphically
%     removing arrival time and amplitude outliers from the core-diffracted
%     wave analysis RESULTS struct generated by either CMB_1ST_PASS,
%     CMB_CLUSTERING or CMB_2ND_PASS.  The user is also allowed to select
%     stations in a specific region based on limiting by azimuth and
%     distance.  Plots and info are saved during the analysis by the user.
%
%     RESULTS=CMB_OUTLIERS(RESULTS,ODIR) sets the output directory
%     where the figures and RESULTS struct is saved.  By default ODIR is
%     '.' (the current directory).  You may set ODIR to FALSE for no
%     written output.
%
%     RESULTS=CMB_OUTLIERS(RESULTS,ODIR,FIGDIR) allows saving figures to a
%     different directory than ODIR (where the RESULTS struct is saved).
%     You may set FIGDIR to FALSE to skip saving figures.
%
%    Notes:
%     - Outliers are reset each time CMB_OUTLIERS is ran on a RESULTS
%       struct.
%
%    Examples:
%     % Typical alignment and refinement workflow:
%     results=cmb_1st_pass;
%     results=cmb_clustering(results);
%     results=cmb_outliers(results);
%
%    See also: PREP_CMB_DATA, CMB_1ST_PASS, CMB_2ND_PASS, SLOWDECAYPAIRS,
%              SLOWDECAYPROFILES, MAP_CMB_PROFILES, CMB_CLUSTERING

%     Version History:
%        Dec. 12, 2010 - added docs
%        Jan.  6, 2011 - catch empty axis handle breakage
%        Jan. 13, 2011 - output ground units in .adjustclusters field
%        Jan. 16, 2011 - split off clustering to cmb_clustering, menu
%                        rather than forcing user to cycle through
%        Jan. 18, 2011 - .time field, no setting groups as bad
%        Jan. 26, 2011 - no travel time corrections for synthetics, use 2
%                        digit cluster numbers, update for 2 plot arrcut
%        Jan. 29, 2011 - prepend datetime to output names, fix Sdiff
%                        corrections bug
%        Jan. 31, 2011 - odir & figdir inputs
%        Feb. 12, 2011 - include snr-based arrival time error
%        Mar.  6, 2011 - coloring in arrcut/ampcut, undo all button
%        Mar. 10, 2011 - distance/azimuth window
%        Mar. 30, 2011 - minor doc update
%        Apr.  3, 2011 - update .time field of skipped
%        Apr. 22, 2011 - update for finalcut field
%        May  19, 2011 - undo works now
%        Mar.  1, 2012 - octave ascii save workaround
%        Mar.  5, 2012 - allow no written output
%        Mar. 11, 2013 - directory input (read indir/*.mat), selection list
%        Jan. 27, 2014 - abs path fix & reduced filesep/fullfile calls
%
%     Written by Garrett Euler (ggeuler at wustl dot edu)
%     Last Updated Jan. 27, 2014 at 13:35 GMT

% todo:

% check nargin
error(nargchk(1,3,nargin));

% directory separator
fs=filesep;

% handle directory input
if(isstring(results))
    if(~isabspath(results)); results=[pwd fs results]; end
    if(isdir(results))
        files=xdir([results fs '*.mat']);
        clear results;
        for i=1:numel(files)
            results(i)=load([files(i).path files(i).name]);
        end
    end
end

% check results
error(check_cmb_results(results));

% default odir & figdir
if(nargin<2 || isempty(odir)); odir='.'; end
if(islogical(odir) && isscalar(odir) && odir); odir='.'; end
if(nargin<3 || isempty(figdir)); figdir=odir; end
if(islogical(figdir) && isscalar(figdir) && figdir); figdir='.'; end

% check odir & figdir
if(~isstring(odir) && ~(islogical(odir) && isscalar(odir)))
    error('seizmo:cmb_outliers:badInput',...
        'ODIR must be a string or TRUE/FALSE!');
elseif(~isstring(figdir) && ~(islogical(figdir) && isscalar(figdir)))
    error('seizmo:cmb_outliers:badInput',...
        'FIGDIR must be a string or TRUE/FALSE!');
end
if(islogical(odir)); out=odir; else out=true; end
if(islogical(figdir)); figout=odir; else figout=true; end

% make sure odir/figdir exists (create it if it does not)
if(out)
    [ok,msg,msgid]=mkdir(odir);
    if(~ok)
        warning(msgid,msg);
        error('seizmo:cmb_outliers:pathBad',...
            'Cannot create directory: %s',odir);
    end
elseif(~out && ~nargout)
    error('seizmo:cmb_outliers:badInput',...
        'Output variable must be assigned when no written output!');
end
if(figout)
    [ok,msg,msgid]=mkdir(figdir);
    if(~ok)
        warning(msgid,msg);
        error('seizmo:cmb_outliers:pathBad',...
            'Cannot create directory: %s',figdir);
    end
end

% select events
datelist=char({results.runname}.');
s=listdlg('PromptString','Select events:',...
          'InitialValue',1:numel(results),...
          'ListSize',[170 300],...
          'ListString',datelist);

% error if none selected
if(isempty(s))
    error('seizmo:cmb_outliers:noDirsSelected',...
        'No earthquakes selected!');
end

% loop over each event
for i=1:numel(s)
    % display name
    disp(results(s(i)).runname);
    
    % time (for skipped)
    results(s(i)).time=datestr(now);
    
    % abandon events we skipped
    if(isempty(results(s(i)).useralign)); continue; end
    
    % arrival & amplitude info
    [dd,az,ev,st]=getheader(results(s(i)).useralign.data,...
        'gcarc','az','ev','st');
    arr=results(s(i)).useralign.solution.arr;
    if(results(s(i)).synthetics)
        carr=arr;
    else
        switch results(s(i)).phase
            case 'Pdiff'
                carr=arr-results(s(i)).corrections.ellcor...
                    -results(s(i)).corrections.crucor.prem...
                    -results(s(i)).corrections.mancor.hmsl06p.upswing;
            case {'SHdiff' 'SVdiff'}
                carr=arr-results(s(i)).corrections.ellcor...
                    -results(s(i)).corrections.crucor.prem...
                    -results(s(i)).corrections.mancor.hmsl06s.upswing;
        end
    end
    snr=results(s(i)).usersnr.snr;
    snr=snr(snr>=results(s(i)).usersnr.snrcut);
    snr(results(s(i)).userwinnow.cut)=[];
    if(isfield(results(s(i)),'finalcut'))
        snr=snr(results(s(i)).finalcut);
    end
    arrerr=sqrt((results(s(i)).useralign.solution.arrerr).^2 ...
        +(max(1./results(s(i)).filter.corners)...
        ./(2*pi).*snr2phaseerror(snr)).^2);
    amp=results(s(i)).useralign.solution.amp;
    camp=amp./results(s(i)).corrections.geomsprcor;
    amperr=results(s(i)).useralign.solution.amperr;
    
    % default to all non-outliers
    results(s(i)).outliers.bad=false(...
        numel(results(s(i)).useralign.data),1);
    
    % loop over good clusters
    for j=find(results(s(i)).usercluster.good(:)')
        % current cluster index as a string
        sj=num2str(j,'%02d');
        
        % preallocate struct
        results(s(i)).outliers.cluster(j).arrcut=...
            struct('bad',[],'cutoff',[]);
        results(s(i)).outliers.cluster(j).ampcut=...
            struct('bad',[],'cutoff',[]);
        results(s(i)).outliers.cluster(j).errcut=...
            struct('bad',[],'cutoff',[]);
        
        % loop until user is happy overall
        % with outlier analysis of this cluster
        happyoverall=false;
        arrcnt=0; ampcnt=0; errcnt=0; delazcnt=0;
        while(~happyoverall)
            % get current cluster population
            good=find(results(s(i)).usercluster.T==j ...
                & ~results(s(i)).outliers.bad);
            pop=numel(good);
            
            % require at least 2 members in good standing
            % - otherwise set remaining as outliers and skip
            if(pop<2)
                warning('seizmo:cmb_outliers:tooFewGood',...
                    ['Cluster ' sj ' has <2 good members. Skipping!']);
                results(s(i)).outliers.bad(good)=true;
                happyoverall=true;
                continue;
            end
            
            % ask user what to do
            choice=menu(...
                ['CLUSTER ' sj ': Remove what outliers?'],...
                'Arrival Time',...
                'Arrival Time Error',...
                'Amplitude',...
                'Azimuth/Distance',...
                'Undo All',...
                'Continue');
            
            % action based on choice
            switch choice
                case 1 % arr
                    arrcnt=arrcnt+1;
                    [bad,cutoff,ax]=...
                        arrcut(dd(good),carr(good),[],1,arrerr(good),[],...
                        z2c(az(good),hsv(64),[0 360]));
                    results(s(i)).outliers.bad(good(bad))=true;
                    results(s(i)).outliers.cluster(j).arrcut.bad{arrcnt}=good(bad);
                    results(s(i)).outliers.cluster(j).arrcut.cutoff(arrcnt)=cutoff;
                    if(figout && ishandle(ax(1)))
                        saveas(get(ax(1),'parent'),[figdir fs ...
                            datestr(now,30) '_' results(s(i)).runname ...
                            '_cluster_' sj '_arrcut_' num2str(arrcnt) ...
                            '.fig']);
                    end
                    if(ishandle(ax(1)))
                        close(get(ax(1),'parent'));
                    end
                case 2 % arrerr
                    errcnt=errcnt+1;
                    [bad,cutoff,ax]=errcut(dd(good),arrerr(good));
                    results(s(i)).outliers.bad(good(bad))=true;
                    results(s(i)).outliers.cluster(j).errcut.bad{errcnt}=good(bad);
                    results(s(i)).outliers.cluster(j).errcut.cutoff(errcnt)=cutoff;
                    if(figout && ishandle(ax))
                        saveas(get(ax,'parent'),[figdir fs ...
                            datestr(now,30) '_' results(s(i)).runname ...
                            '_cluster_' sj '_errcut_' num2str(errcnt) ...
                            '.fig']);
                    end
                    if(ishandle(ax))
                        close(get(ax,'parent'));
                    end
                case 3 % amp
                    ampcnt=ampcnt+1;
                    [bad,cutoff,ax]=ampcut(dd(good),camp(good),[],1,amperr(good),[],...
                        z2c(az(good),hsv(64),[0 360]));
                    results(s(i)).outliers.bad(good(bad))=true;
                    results(s(i)).outliers.cluster(j).ampcut.bad{ampcnt}=good(bad);
                    results(s(i)).outliers.cluster(j).ampcut.cutoff(ampcnt)=cutoff;
                    if(figout && ishandle(ax))
                        saveas(get(ax,'parent'),[figdir fs ...
                            datestr(now,30) '_' results(s(i)).runname ...
                            '_cluster_' sj '_ampcut_' num2str(ampcnt) ...
                            '.fig']);
                    end
                    if(ishandle(ax))
                        close(get(ax,'parent'));
                    end
                case 4 % distance/azimuth
                    delazcnt=delazcnt+1;
                    [bad,azlim,ddlim,ax]=delazcut(ev(1,1:2),st(good,1:2),...
                        [],[],results(s(i)).usercluster.color(j,:));
                    results(s(i)).outliers.bad(good(bad))=true;
                    results(s(i)).outliers.cluster(j).delazcut.bad{delazcnt}=good(bad);
                    results(s(i)).outliers.cluster(j).delazcut.azlim{delazcnt}=azlim;
                    results(s(i)).outliers.cluster(j).delazcut.ddlim{delazcnt}=ddlim;
                    if(figout && ishandle(ax))
                        saveas(get(ax,'parent'),[figdir fs ...
                            datestr(now,30) '_' results(s(i)).runname ...
                            '_cluster_' sj '_delazcut_' num2str(delazcnt) ...
                            '.fig']);
                    end
                    if(ishandle(ax))
                        close(get(ax,'parent'));
                    end
                case 5 % undo all
                    results(s(i)).outliers.cluster(j).arrcut=...
                        struct('bad',[],'cutoff',[]);
                    results(s(i)).outliers.cluster(j).ampcut=...
                        struct('bad',[],'cutoff',[]);
                    results(s(i)).outliers.cluster(j).errcut=...
                        struct('bad',[],'cutoff',[]);
                    arrcnt=0; ampcnt=0; errcnt=0;
                    results(s(i)).outliers.bad(...
                        results(s(i)).usercluster.T==j)=false;
                case 6 % continue
                    happyoverall=true;
                    continue;
            end
        end
    end
    
    % time
    results(s(i)).time=datestr(now);
    
    % save results
    if(out)
        tmp=results(s(i));
        if(isoctave)
            save([odir fs datestr(now,30) '_' results(s(i)).runname ...
                '_outliers_results.mat'],'-7','-struct','tmp');
        else % matlab
            save([odir fs datestr(now,30) '_' results(s(i)).runname ...
                '_outliers_results.mat'],'-struct','tmp');
        end
    end
end

end
