function [pflag, vprofile, trefl, cfg] = pingedit(data, varargin)
% [pflag, vprofile, trefl, cfg] = pingedit(data, varargin)
% This program returns profile flags representing the applied editing.
%    * pflags is a matrix (same size as data.vel) with zero for unedited bins,
%         value     type of editing
%           2^0        bin was NaN to start with (after 'pflag' was applied)
%           2^1        bin fals below correlation threshold
%           2^2        shallow pings edited based on low cor (in a deeper bin)
%           2^3        bin excedes error velocity threshold      
%           2^4        refavg outlier (caveat: see next note)
%           2^5        "weak" profile removed (too few bins left)
%                           (this step occurs after other editing is applied)
%
%    * reference layer averaging may be applied, so it is expected that
%        the data to be edited will be in small enough chunks that genuine
%        changes in the ocean won't be edited out as outliers from reference 
%        layer averaging. If you want to use this on a big chunk of data 
%        for some reason (eg. several hours), disable reference layer 
%        averaging (set 'refcufoff' to []).  Note that reference layer editing
%        is strongly recommended for OS data in Narrowband mode under bad 
%        weather conditions.
% 
% input:
% data is a structure with field used as needed, specifically
% data.vel is [Nbins x Ntimes] in m/s, probably umeas, vmeas, or fmeas 
%          (probably not beam velocity)
% data.emeas is error velocity (same size as vel)
% data.cor is some correlation measure, such as 
%        min(min(min(min(data.cor1),data.cor2),data.cor3),data.cor4);
%
% varargin name, value pairs control the editing
%
% name         values            default     what
% -------      ---------         --------    ---------------------------------
% 'fieldname'  string             'vel'      field of structure "data" which
%                                                will be edited (eg- could set
%                                                to 'fmeas', 'umeas', 'vmeas'
% 'index'       [bins]            [1 end]    use only these time indices 
% 'cfg'         structure          []        set editing parameters in a 
%                                                structure; setting any of the
%
%
%  usage: 
%  help pingedit                         %gives this message
%  pingedit                              %show a list of editing parameters 


%JH 2001/11/16


%set defaults
cfg.fieldname = 'vel';
cfg.index     = 'all';
                                                                   %%disable
% (0)
cfg.cor_cutoff          =       140;                                 % []
% (1)                                                                %
cfg.sbr_corcutoff       =       165;                                 % []
cfg.sbr_bincheck        =       5;                                   %
cfg.sbr_binreject       =       [1 2 3];                             % 
% (2)                                                                %
cfg.err_cutoff          =      .5;                                   % []
% (3)                                                                %
cfg.rl_range            =       'all';                               % []
cfg.refavg_bincutoff       =      .5;                                   % []
cfg.refavg_binstddev       =       3;                                   % []
% (4)                                                                %
cfg.weakprof_binsreq    =       4;                                   % []
cfg.weakprof_testrange  =       'all';                               % 
cfg.pflag               =       'allzero';                           %
cfg.help                =       [];


%find out if 'cfg' was used: if so, override current cfg defaults

%fill cfg
cfgplace = [];
if (nargin > 1)
   for ii=1:2:length(varargin);
      if strmatch('cfg', varargin{ii})
         cfgplace = ii;
      end
   end
   if ~isempty(cfgplace)
      cfg = fillstruct(cfg, varargin{cfgplace + 1});
      %now delete that field from varargin
      varargin([cfgplace, cfgplace + 1]) = [];
   end
   if (~isempty(varargin))
      cfg = fillstruct(cfg, varargin, 0);
   end
else
   usage
   return
end

if (strmatch('help', cfg.help)),    help pingdata, end
if (strmatch('listopt', cfg.help)), usage, end
if (strmatch('lstcfg', cfg.help)),  cfg, end


% set up correct field
if (~strcmp(cfg.fieldname, 'vel'))
   data.vel = getfield(data, cfg.fieldname);
end
% get indices
if (isstr(cfg.index))  %uninitialized
   index = 1:length(data.vel(1,:));
else
   badi = find(cfg.index > length(data.vel(1,:)));
   cfg.index(badi) = [];
end

if (isstr(cfg.index))
   cfg.index = [1:length(data.vel(1,:))];
end
data.vel = data.vel(:, cfg.index);

if (isfield(data, 'cor')),    data.cor = data.cor(:, cfg.index);  end
if (isfield(data, 'emeas')),    data.emeas = data.emeas(:, cfg.index);  end


% fix other things with defaults that depend on size of data.vel
if (isstr(cfg.rl_range) )
   cfg.rl_range = [1 length(data.vel(:,1))];
end
if (isstr(cfg.weakprof_testrange))
   cfg.weakprof_testrange = [1 length(data.vel(:,1))];
end
if (isstr(cfg.pflag))
   cfg.pflag = 0*data.vel;
end


if (isreal(data.vel))
   NAN = NaN;
else
   NAN = NaN + i*NaN;
end





%         value     type of editing
%           2^0        bin falls below correlation threshold
%           2^1        shallow pings edited based on low cor (in a deeper bin)
%           2^2        bin excedes error velocity threshold      
%           2^3        refavg outlier (caveat: see next note)
%           2^4        "weak" profile removed (too few bins left)
%                           (this step occurs after other editing is applied)
%%% values for profile flags:    (0 or this value)
nan_binval      = 2^0;
cor_binval      = 2^1;
sbr_binval      = 2^2;
err_binval      = 2^3;
refavg_binval   = 2^4;
weak_binval     = 2^5;


%apply pflag from cfg 
data.vel(cfg.pflag ~= 0) = NAN;


%%% masks are all ones if editing for that type is turned off 
%%%    (checked IN function)
%
% toss bins out from the independent checks

nan_binmask = 0*data.vel; nan_binmask(isnan(data.vel)) = nan_binval;
cor_binmask = cor_binedit(data, cfg, cor_binval);
sbr_binmask = cor_sbredit(data, cfg, sbr_binval);
err_binmask = err_binedit(data, cfg, err_binval);
% find out what's bad so far, test for weak profiles based on what's left
badmask = nan_binmask + cor_binmask + sbr_binmask + err_binmask;
weakprof_mask1 = weakprof_edit(badmask, cfg, weak_binval);
badmask = badmask + weakprof_mask1;
% if refavg editing is enabled, do that and check for weak profiles again
refavg_binmask = zeros(size(data.vel));
weakprof_mask2 = zeros(size(data.vel));
if ~isempty(cfg.rl_range) &...
      (~isempty(cfg.refavg_bincutoff) | ~isempty(cfg.refavg_binstddev))
   nanmask = ~badmask;  %1 is good; 0 is bad
   nanmask(~nanmask) = NAN;  %1 is good, NaN is bad
    [refavg_binmask, vprofile, trefl]  = refavg_binedit(nanmask.*data.vel,...
                    cfg, refavg_binval);
   badmask = badmask + refavg_binmask;
   weakprof_mask2 = weakprof_edit(badmask, cfg, weak_binval);
end

%add refavg_profedit;??


weakprof_mask =...
    (weakprof_mask1/weak_binval | weakprof_mask2/weak_binval)*weak_binval;
pflag = nan_binmask + cor_binmask + sbr_binmask + err_binmask + ...
        refavg_binmask + weakprof_mask;



vv = 0;
if (vv)
   for pp = 1:100
      [cor_binmask(:,pp) sbr_binmask(:,pp) err_binmask(:,pp) ...
       refavg_binmask(:,pp)  weakprof_mask1(:,pp) weakprof_mask2(:,pp) ...
      badmask(:,pp)]
      paused
   end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cor_binmask = cor_binedit(data, cfg, cor_binval)
cor_binmask = zeros(size(data.vel));   
if (isfield(data, 'cor') & ~isempty(cfg.cor_cutoff))
   cor_binmask(data.cor < cfg.cor_cutoff) = cor_binval;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cor_sbrmask = cor_sbredit(data, cfg, sbr_binval);
cor_sbrmask = zeros(size(data.vel));
%   %profilewise  (identify profile, reject bins in that profile)
if (isfield(data, 'cor') & ...
    ~isempty(cfg.sbr_bincheck) & ...
    ~isempty(cfg.sbr_bincheck) & ...
    ~isempty(cfg.sbr_binreject))
   
   badprofi = find(data.cor(cfg.sbr_bincheck,:) < cfg.sbr_corcutoff);
   cor_sbrmask(cfg.sbr_binreject, badprofi) = sbr_binval;
end

%%%%%%%%%%%%%%%%%%%%%%

function err_binmask = err_binedit(data, cfg, err_binval)
err_binmask = zeros(size(data.vel));   
if (~isempty(cfg.err_cutoff))
   err_binmask(abs(data.emeas) > cfg.err_cutoff) = err_binval;
end

%%%%%%%%%%%%%%%%%%%%%%
function  weakprof_mask = weakprof_edit(badmask, cfg, weak_binval);
% remove profiles with too few bins in specified range
weakprof_mask = zeros(size(badmask));
if (~isempty(cfg.weakprof_binsreq) & ~isempty(cfg.weakprof_testrange))
   profbins = cfg.weakprof_testrange(1) : cfg.weakprof_testrange(2);
   weakprofs = sum(badmask(profbins,:) == 0) < cfg.weakprof_binsreq;
   weakprof_mask(:,weakprofs) = weak_binval;
end

%%%%%%%%%%%%%%%%%%%%%%

function [refavg_binmask, vprofile, trefl] = ...
    refavg_binedit(vel, cfg, refavg_binval)


refavg_binmask = zeros(size(vel));
if ~isempty(cfg.rl_range) &...
      (~isempty(cfg.refavg_bincutoff) | ~isempty(cfg.refavg_binstddev))

   rlbins = cfg.rl_range(1):cfg.rl_range(2);
   vprofile = [];
   trefl = [];
   [vprofile, trefl]   =  refavg(vel,  rlbins, 3);
   vprofileM = repmat(vprofile, 1, length(vel(1,:)));
   treflM = repmat(trefl, length(vel(:,1)), 1);
   if ~isempty(cfg.refavg_bincutoff)
      refavg_binmask(abs(vel - (vprofileM + treflM)) > cfg.refavg_bincutoff) = ...
          refavg_binval;
   end
   if ~isempty(cfg.refavg_binstddev)
      tmp = vel(rlbins,:);
      [junk, stddev] = mstdgap(tmp(:));
      refavg_binmask(...
         abs(vel - (vprofileM+treflM)) > cfg.refavg_binstddev*stddev) = ...
          refavg_binval;
   end
   refavg_binmask(isnan(refavg_binmask)) = refavg_binval;
end

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%print 
function usage()

str =[...
{' Type "help pingedit" for an overview.  '};
{' This is a list of editing parameters, their effect, and their defaults'};
{'                                                                                   '};
{' name         values            default     what                                   '};
{' -------      ---------         --------    ---------------------------------      '};
{' ''fieldname''  string             ''vel''      field of structure "data" which    '};
{'                                                will be edited (eg- could set      '};
{'                                                to ''fmeas'', ''umeas'', ''vmeas'' '};
{' ''index''       [bins]            [1 end]    use only these time indices          '};
{' ''cfg''         structure          []        set editing parameters in a          '};
{'                                                structure; setting any of the      '};
{'                                                following values overrides the     '};
{'                                                structure.                         '};
{' ''help''        string                        be verbose in one way:              '};
{'                                    []           - do not be verbose               '};
{'                                  '' help'',     - same as ''help pingdata''       '};
{'                                  ''listopt''    - show this list of editing       '};
{'                                                     options and their effects     '};
{'                                  ''showcfg''    - echo configuration used         '};
{' ''pflag''      [Nbins x Ntimes]   0          apply this profile flag matrix       '};
{'                                                before implimenting editing        '};
{' ''rl_range''   [startbin endbin]  [1 end]    bins for refavg averaging            '};
{'                                                * set to [] to disable             '};
{'                                                  ref layer avg editing at all     '};
{' ''refavg_bincutoff'' float        .5        reject bins with velocity too        '};
{'                                                far away from refavg (m/s)         '};
{'                                                * set to [] to disable this        '};
{' ''refavg_binstddev''  float        3        reject bins with velocity too       '};
{'                                                far away from refavg in            '};
{'                                                terms of number of std devs        '};
{'                                                (std is from reference layer)      '};
{' ''err_cutoff''    float            .5        reject bins with abs(errvel)         '};
{'                                                too great                          '};
{'                                                * set to [] to disable             '};
{' ''cor_cutoff''   integer           140       reject data with cor below this      '};
{'                                                * set to [] to disable             '};
{' ''sbr_corcutoff''integer           165       shallow bin reject: identify         '};
{'                                                profiles by low correlation        '};
{'                                                at some bin; reject top bins       '};
{'                                                * est to [] to disable             '};
{' ''sbr_bincheck'' integer           2         shallow bin reject: look at the      '};
{'                                                correlation in this bin: if        '};
{'                                                it is below sbr_corcutoff,         '};
{'                                                delete some bins (                 '};
{' ''sbr_binreject'' [bins]        [1:3]        shallow bin reject: reject these     '};
{'                                                bins in profiles identified        '};
{'                                                by sbr_bincheck and                '};
{'                                                sbr_corcutoff                      '};
{' ''weakprof_binsreq'' integer  2              weak profile rejection: reject       '};
{'                                                profiles with fewer than           '};
{'                                                this number of bins in the         '};
{'                                                part of the profile set by         '};
{'                                                weakprof_testrange                 '};
{'                                                * set to [] to disable             '};
{' ''weakprof_testrange'' integer   [1 end]     weak profile rejection: tag a        '};
{'                                                profile as weak if there are       '};
{'                                                fewer than weakprof_binsreq        '}];

for ii=1:length(str(:,1))
   fprintf('%s\n',str{ii})
end
