function [pflag, vprofile, trefl, resid, binmasks] = pingedit(instclass, data, varargin)
% [pflag, vprofile, trefl, resid, binmasks]= pingedit(instclass, 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
%         nan_binval          = 2^0;%   original velocity was NaN 
%         cor_binval          = 2^1;%   correlation falls below threshold 
%         sw_binval           = 2^1;%   spectral width excedes threshold
%         err_binval          = 2^2;%   error velocity excedes threshold      
%         ocean_uv_binval     = 2^3;%   ocean velocity excedes threshold      
%         weakprof_binval     = 2^4;%   "weak" profile (too few bins left)
%         refavg_outlier      = 2^5;%   outlier from "mean"
%         useredit_binval     = 2^6;%   user-written editing assigned this value
%
%    * 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.
%
%
%    * pingedit values are initialized using "get_pedefaults" and can be 
%        overridden  by two methods.  
%        (1) name, value pair: 'editcfg', 'new_editcfg', where 'new_editcfg'
%             contains any of the same fields in get_pedefaults 
%        (2) name, value pairs where the name and value are the field 
%             and value from the structure  output by get_pedefaults.  
%        Both methods are implimented but (2) overrides (1).
%        An example of (1) is shown below
%                  
% to use:
% (1) get the data: 
%       data=read(instclass, filename, 'vars','all');
% (2) pre-edit the data  %removes amplitude spikes, looks for the bottom
%        [data, bt_bmask, bt_wmask, hiamp]  = prepedit(instclass, data);    
% (3) run pingedit: 
%     [pflag, vprofile, trefl, cfg, resid] = pingedit(os, data);
%     specifically: -  apply simple cutoffs (eg. error velocity)
%                  -  apply complex schemes (eg. weak profile deletion)
%
% To change editing parameters, replace (3) above with
% (3a) editcfg = get_pedefaults(instclass, pingtype)
% (3b) (look at editcfg and override if necessary)
% (3c) call pingedit:
%     [pflag] = pingedit(os, data 'editcfg', editcfg);
%
%--------------------------------
% for help on ping editing defaults, type 'help get_pedefaults'
% for help on pre-pingediting, type 'help prepedit'


%JH 2001/11/16
%JH 2003/12/16 modified to pull out get_pedefaults and prepedit.

%initialize outputs as empty
pflag=[];, vprofile=[];, trefl=[];, cfg=[];, resid=[];

%initialize cfg

if rem(nargin - 2, 2) ~= 0
   help(mfilename)
   fprintf('enter name, value _pairs_\n')
   return
end
   

% get ping edit defaults
editcfg = get_pedefaults(instclass, data.config.pingtype);
if nargin == 0
   get_pedefaults %returns
end

%replace pingedit defaults  with values from editcfg structure (if exists)
cfgplace = [];
if (exist('varargin','var')==1)
   for ii=1:2:length(varargin);
      if strmatch('editcfg', varargin{ii})
         cfgplace = ii;
      end
   end
   editcfg = fillstruct(editcfg, varargin{cfgplace+1});
   varargin([cfgplace, cfgplace + 1]) = [];
end
% now apply the rest of varargin
% initialize showverbose now, since it isn't in editcfg
editcfg.showverbose = 0;
editcfg.pflag = zeros(size(data.cuv));  % xx

if (~isemptys(varargin))
   editcfg = fillstruct(editcfg, varargin, 0);
end


% make cor_corcutoff and cor_testbins into cellarrays
if isfield(editcfg, 'cor_cutoff')
   if iscell(editcfg.cor_cutoff)
      if (length(editcfg.cor_cutoff) ~= length(editcfg.cor_testbins))
         fprintf('cor_cutoff, cor_testbins cellarrays have unequal lengths\n')
         return
      end
   else %turn into cellarrays
      if ~isemptys(editcfg.cor_cutoff)
         editcfg.cor_cutoff = {editcfg.cor_cutoff}; 
      end
      if ~isemptys(editcfg.cor_testbins), 
         editcfg.cor_testbins = {editcfg.cor_testbins};
      end
   end
end


% fix other things with defaults that depend on size of data.cuv
%% this has to be here and in avg in case pingedit is called NOT from avg
if (isstr(editcfg.rl_bins) )
   if strcmp(editcfg.rl_bins, 'all')
      editcfg.rl_bins = [1:length(data.cuv(:,1))]';
   else
      fprintf('string %s not allowed for "rl_bins" field\n', editcfg.rl_bins)
      return
   end
end
if (isstr(editcfg.weakprof_testbins))
   if strcmp(editcfg.weakprof_testbins, 'all')
      editcfg.weakprof_testbins = [1:length(data.cuv(:,1))]';
   else
      fprintf('string %s not allowed for "weakprof_testbins" field\n', ...
              editcfg.weakprof_testbins)
      return
   end
end



%apply pflag from editcfg : specify complex if necessary
if (isreal(data.cuv)),  NAN = NaN; else, NAN = NaN + i*NaN; end

% type of editing
% (1) from the instrument
binvals.nan                      = 2^0;%   original velocity was NaN 
% (2) cutoff_functions
binvals.correlation              = 2^1;%   correlation falls below threshold 
binvals.spectral_width           = 2^1;%   spectral width excedes threshold
binvals.error_velocity           = 2^2;%   error velocity excedes threshold 
binvals.ocean_velocity           = 2^3;    %   ocean velocity excedes threshold 
% (3) pingedit_functions
binvals.weak_profile             = 2^4;%   "weak" profile (too few bins left)
binvals.refavg_outlier           = 2^5;%   refavg outlier editing
% (4) user-written editing
binvals.useredit                 = 2^6;%   user-written editing gets this value
%% summarize these for showverbose to use
binmasks = struct;
% obtain and apply the masks.  
% (1) Be careful to define masks so they are not logical.  We need to add them!!


%% BAD?? pflag = zeros(size(data.cuv));
%initialize with bad values passed in, as from prepedit output
pflag = editcfg.pflag;  
% add the bad velocities to the pot
pflag(isnan(data.cuv)) = 1;



%% find out which functions to run
% (2) % simple cutoff engines
for functioni = 1:length(editcfg.cutoff_functions)
   % get function name and determine binval
   function_base =  editcfg.cutoff_functions{functioni};
   function_name = sprintf('%s_cutoff', function_base);
   binval = getfield(binvals, function_base);

   if editcfg.verbose
      fprintf('about to run function %s\n', function_name)
   end

   % run the function; create the mask with binval where flagged
   mask = binval*feval(function_name, data, editcfg);
   % add it to the pflag bits
   if ~isemptys(mask)
      
      pflag = bitor(pflag, mask);
      % save for later debugging
      eval(sprintf('binmasks.%s = mask;', function_base))
   end
end
% (3)  user edit
for functioni = 1:length(editcfg.useredit_functions)
   function_name = editcfg.useredit_functions{functioni};
   binval = binvals.useredit;
   if editcfg.verbose
      fprintf('about to run function %s\n', function_name)
   end
   mask = binval*feval(function_name, data, editcfg);
   if ~isemptys(mask)
      pflag = bitor(pflag, mask);
      binmasks.useredit = mask;
   end
end


% (4) % weak profile and refavg
for functioni = 1:length(editcfg.pingedit_functions)
   function_base = editcfg.pingedit_functions{functioni};
   function_name = sprintf('%s_edit', function_base);
   binval = getfield(binvals, function_base);

   if editcfg.verbose
      fprintf('about to run function %s\n', function_name)
   end
   % pingedit functions operate on velocity data with pflags applied
   vel = dotmmask(pflag).*data.cuv;
   mask = binval*feval(function_name, vel, editcfg);
   if ~isemptys(mask)
      pflag = bitor(pflag, mask);
      eval(sprintf('binmasks.%s = mask;', function_base))
   end
   
end

[vprofile, trefl, resid] = refavg(dotmmask(pflag).*data.cuv, ...
                                  editcfg.rl_bins, 3);


if (editcfg.showverbose)
   
   showcfg.ord_name = 'binnum';
   showcfg.abc_name = 'index';
   showcfg.fn = 'fn';
   showdata.binnum = [1:length(data.depth)]';
   showdata.index = [1:length(data.ens_num)];
   
   f4
   binval_fieldnames = fieldnames(binmasks);
   for plotnum = 1:length(binval_fieldnames)
      fn = binval_fieldnames{plotnum};
      showdata.fn =  getfield(binmasks, fn);
      aplotit(showdata, showcfg, 'fn', 'fn', 'cur_ax', subplot(4,2,plotnum),...
              'comment', fn, 'astitle', 1);
      plotnum = plotnum+1;
   end
   subplot(4,2,plotnum)
   plot(100*sum(pflag'==0)/length(pflag(1,:)), showdata.binnum);ra(gca)

end


