% UTP_AO_AVERAGE a set of UTPs for the ao/average method
%
% M Hueller 24-02-1
%
% $Id: utp_ao_average.m,v 1.2 2011/04/17 15:46:21 ingo Exp $
%

% <MethodDescription>
%
% The average method of the ao class averages aos point-by-point
%
% </MethodDescription>

function results = utp_ao_average(varargin)
  
  % Check the inputs
  if nargin == 0
    
    % Some keywords
    class   = 'ao';
    mthd    = 'average';
    
    results = [];
    disp('******************************************************');
    disp(['****  Running UTPs for ' class '/' mthd]);
    disp('******************************************************');
    
    % Test AOs
    [at1,at2,at3,at4,at5,at6,atvec,atmat] = eval(['get_test_objects_' class]);
    
    % Exception list for the UTPs:
    [ple1,ple2,ple3,ple4,ple5,ple6] = get_test_ples();
    
    % Set the same y-units for the test AOs
    at1.setYunits('uHz');
    at2.setYunits('uHz');
    at3.setYunits('uHz');
    at4.setYunits('uHz');
    at5.setYunits('uHz');
    at6.setYunits('uHz');
    
    % Run the tests
    results = [results utp_01];    % getInfo call
    results = [results utp_02];    % Vector input
    results = [results utp_03];    % Matrix input
    results = [results utp_04];    % List input
    results = [results utp_05];    % Test with mixed input
    results = [results utp_06(mthd, at1, at5, ple2)];    % Test history is working
    results = [results utp_07];    % Test the modify call does not work
    results = [results utp_08];    % Test input data shape == output data shape
    results = [results utp_18];    % Test with single input
    
    
    disp('Done.');
    disp('******************************************************');
    
  elseif nargin == 1 % Check for UTP functions
    if strcmp(varargin{1}, 'isutp')
      results = 1;
    else
      results = 0;
    end
  else
    error('### Incorrect inputs')
  end
  
  %% UTP_01
  
  % <TestDescription>
  %
  % Tests that the getInfo call works for this method.
  %
  % </TestDescription>
  function result = utp_01
    
    
    % <SyntaxDescription>
    %
    % Test that the getInfo call works for no sets, all sets, and each set
    % individually.
    %
    % </SyntaxDescription>
    
    try
      % <SyntaxCode>
      % Call for no sets
      io(1) = eval([class '.getInfo(''' mthd ''', ''None'')']);
      % Call for all sets
      io(2) = eval([class '.getInfo(''' mthd ''')']);
      % Call for each set
      for kk=1:numel(io(2).sets)
        io(kk+2) = eval([class '.getInfo(''' mthd ''', ''' io(2).sets{kk} ''')']);
      end
      % </SyntaxCode>
      stest = true;
    catch err
      disp(err.message)
      stest = false;
    end
    
    % <AlgoDescription>
    %
    % 1) Check that getInfo call returned an minfo object in all cases.
    % 2) Check that all plists have the correct parameters.
    %
    % </AlgoDescription>
    
    atest = true;
    if stest
      % <AlgoCode>
      % check we have minfo objects
      if isa(io, 'minfo')
        % SET 'None'
        if ~isempty(io(1).sets), atest = false; end
        if ~isempty(io(1).plists), atest = false; end
        % Check all Sets
        if ~any(strcmpi(io(2).sets, 'Default')), atest = false; end
        if numel(io(2).plists) ~= numel(io(2).sets), atest = false; end
        % SET 'Default'
        if io(3).plists.nparams ~= 0, atest = false; end
        % Check key
        % Check default value
        % Check options
      end
      % </AlgoCode>
    else
      atest = false;
    end
    
    % Return a result structure
    result = utp_prepare_result(atest, stest, dbstack, mfilename);
  end % END UTP_01
  
  %% UTP_02
  
  % <TestDescription>
  %
  % Tests that the average method works with a vector of AOs as input.
  %
  % </TestDescription>
  function result = utp_02
    
    % <SyntaxDescription>
    %
    % Tests that the average method works with a vector of AOs as input.
    %
    % </SyntaxDescription>
    
    try
      % <SyntaxCode>
      % Pick units and prefix from those supported
      unit_list = unit.supportedUnits;
      % remove the first empty unit '' from the list, because then is it
      % possible that we add a prefix to an empty unit
      unit_list = unit_list(2:end);
      prefix_list = unit.supportedPrefixes;
      u = unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]);
      
      % tsdata: settings
      fs = 3;
      ts_pl = plist('fs', fs, 'nsecs', 30, 'waveform', 'sine wave', 'f', 10, 'yunits', u);
      ts_a1 = ao(ts_pl);
      ts_a2 = ao(ts_pl);
      ts_a3 = ao(ts_pl);
      ts_a1.setT0('14:00:00');
      ts_a2.setT0('14:00:05');
      ts_a3.setT0('14:00:10');
      
      % tsdata: call the method with a vector of inputs
      ts_out = average([ts_a1, ts_a2, ts_a3]);
      
      % fsdata: settings
      fs_pl = plist('fsfcn', '1./f.^2', 'scale', 'lin', 'nf', 100, 'yunits', u);
      fs_a1 = ao(fs_pl);
      fs_a2 = 2*ao(fs_pl);
      
      % tsdata: call the method with a vector of inputs
      fs_out = average([fs_a1, fs_a2]);
      
      % xydata: settings
      xy_pl = plist('X', [1:50,52:2:100,110:10:1000], 'xyfcn', 'log(x)', 'yunits', u);
      xy_a1 = ao(xy_pl);
      xy_a2 = ao(xy_pl.pset('xyfcn', 'log10(x)'));
      xy_a3 = ao(xy_pl.pset('xyfcn', 'exp(-x)'));
      
      % xydata: call the method with a vector of inputs
      xy_out = average([xy_a1, xy_a2, xy_a3]);
      
      % cdata: settings
      c_pl = plist('yunits', u);
      c_a1 = ao([1 2 3], c_pl);
      c_a2 = ao([4 3 2], c_pl);
      c_a3 = ao([7 4 9], c_pl);
      c_a4 = ao([-4 -1 2], c_pl);
      
      % cdata: call the method with a vector of inputs
      c_out = average([c_a1, c_a2, c_a3, c_a4]);
      
      % Rebuild the objects
      ts_mout = rebuild(ts_out);
      fs_mout = rebuild(fs_out);
      xy_mout = rebuild(xy_out);
      c_mout  = rebuild(c_out);
      
      % </SyntaxCode>
      stest = true;
    catch err
      disp(err.message)
      stest = false;
    end
    
    % <AlgoDescription>
    %
    % 1) Check that the outputs have exactly one AO
    % 2) Check that the outputs have the correct data.
    % 3) Check the rebuilt objects
    %
    % </AlgoDescription>
    
    atest = true;
    
    if stest
      % <AlgoCode>
      % Check that the output is exactly one AO
      if numel(ts_out) ~= 1, atest = false; end
      if numel(fs_out) ~= 1, atest = false; end
      if numel(xy_out) ~= 1, atest = false; end
      if numel(c_out) ~= 1, atest = false; end
      if ~isa(ts_out, 'ao'), atest = false; end
      if ~isa(fs_out, 'ao'), atest = false; end
      if ~isa(xy_out, 'ao'), atest = false; end
      if ~isa(c_out, 'ao'), atest = false; end
      
      % Check output data: x
      if ~isequal(ts_out.x, ts_a1.x), atest = false; end
      if ~isequal(fs_out.x, fs_a1.x), atest = false; end
      if ~isequal(xy_out.x, xy_a1.x), atest = false; end
      
      % Check output data: t0
      if ~isequal(ts_out.t0, ts_a1.t0), atest = false; end
      
      % Check output data: y
      ts_data_in = [ts_a1.y ts_a2.y ts_a3.y];
      if ~isequal(ts_out.y, mean(ts_data_in, 2)), atest = false; end
      fs_data_in = [fs_a1.y fs_a2.y];
      if ~isequal(fs_out.y, mean(fs_data_in, 2)), atest = false; end
      xy_data_in = [xy_a1.y xy_a2.y xy_a3.y];
      if ~isequal(xy_out.y, mean(xy_data_in, 2)), atest = false; end
      c_data_in = [c_a1.y c_a2.y c_a3.y c_a4.y];
      if ~isequal(c_out.y, mean(c_data_in, 2)), atest = false; end
      
      % Check output data: dy
      if ~isequal(ts_out.dy, std(ts_data_in, 0, 2)/sqrt(3)), atest = false; end
      if ~isequal(fs_out.dy, std(fs_data_in, 0, 2)/sqrt(2)), atest = false; end
      if ~isequal(xy_out.dy, std(xy_data_in, 0, 2)/sqrt(3)), atest = false; end
      if ~isequal(c_out.dy, std(c_data_in, 0, 2)/sqrt(4)), atest = false; end
      
      % Check output data: xunits
      if ~isequal(ts_out.xunits, ts_a1.xunits), atest = false; end
      if ~isequal(fs_out.xunits, fs_a1.xunits), atest = false; end
      if ~isequal(xy_out.xunits, xy_a1.xunits), atest = false; end
      
      % Check output data: yunits
      if ~isequal(ts_out.yunits, ts_a1.yunits), atest = false; end
      if ~isequal(fs_out.yunits, fs_a1.yunits), atest = false; end
      if ~isequal(xy_out.yunits, xy_a1.yunits), atest = false; end
      if ~isequal(c_out.yunits, c_a1.yunits), atest = false; end
      
      % Check the re-built objects
      if ~eq(ts_out, ts_mout, ple2), atest = false; end
      if ~eq(fs_out, fs_mout, ple2), atest = false; end
      if ~eq(xy_out, xy_mout, ple2), atest = false; end
      if ~eq(c_out, c_mout, ple2), atest = false; end
      
      % </AlgoCode>
    else
      atest = false;
    end
    
    % Return a result structure
    result = utp_prepare_result(atest, stest, dbstack, mfilename);
  end % END UTP_02
  
  %% UTP_03
  
  % <TestDescription>
  %
  % Tests that the average method works with a matrix of AOs as input.
  %
  % </TestDescription>
  function result = utp_03
    
    % <SyntaxDescription>
    %
    % Tests that the average method works with a matrix of AOs as input.
    %
    % </SyntaxDescription>
    
    try
      % <SyntaxCode>
      % Pick units and prefix from those supported
      unit_list = unit.supportedUnits;
      % remove the first empty unit '' from the list, because then is it
      % possible that we add a prefix to an empty unit
      unit_list = unit_list(2:end);
      prefix_list = unit.supportedPrefixes;
      u = unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]);
      
      % tsdata: settings
      fs = 3;
      ts_pl = plist('fs', fs, 'nsecs', 30, 'waveform', 'sine wave', 'f', 10, 'yunits', u);
      ts_a1 = ao(ts_pl);
      ts_a2 = ao(ts_pl);
      ts_a3 = ao(ts_pl);
      ts_a1.setT0('14:00:00');
      ts_a2.setT0('14:00:05');
      ts_a3.setT0('14:00:10');
      
      % tsdata: call the method with a list of inputs
      ts_out = average([ts_a1, ts_a2, ts_a3; ts_a1, ts_a2, ts_a3]);
      
      % tsdata: settings
      fs_pl = plist('fsfcn', '1./f.^2', 'scale', 'lin', 'nf', 100, 'yunits', u);
      fs_a1 = ao(fs_pl);
      fs_a2 = 2*ao(fs_pl);
      
      % fsdata: call the method with a list of inputs
      fs_out = average([fs_a1, fs_a2; fs_a1, fs_a2]);
      
      % xydata: settings
      xy_pl = plist('X', [1:50,52:2:100,110:10:1000], 'xyfcn', 'log(x)', 'yunits', u);
      xy_a1 = ao(xy_pl);
      xy_a2 = ao(xy_pl.pset('xyfcn', 'log10(x)'));
      xy_a3 = ao(xy_pl.pset('xyfcn', 'exp(-x)'));
      
      % xydata: call the method with a list of inputs
      xy_out = average([xy_a1, xy_a2, xy_a3; xy_a1, xy_a2, xy_a3]);
      
      % cdata: settings
      c_pl = plist('yunits', u);
      c_a1 = ao([1 2 3], c_pl);
      c_a2 = ao([4 3 2], c_pl);
      c_a3 = ao([7 4 9], c_pl);
      c_a4 = ao([-4 -1 2], c_pl);
      
      % cdata: call the method with a list of inputs
      c_out = average([c_a1, c_a2, c_a3, c_a4; c_a1, c_a2, c_a3, c_a4]);
      
      % Rebuild the objects
      ts_mout = rebuild(ts_out);
      fs_mout = rebuild(fs_out);
      xy_mout = rebuild(xy_out);
      c_mout  = rebuild(c_out);
      
      % </SyntaxCode>
      stest = true;
    catch err
      disp(err.message)
      stest = false;
    end
    
    % <AlgoDescription>
    %
    % 1) Check that the outputs have exactly one AO
    % 2) Check that the outputs have the correct data.
    % 3) Check the rebuilt objects
    %
    % </AlgoDescription>
    
    TOL   = 1e-15;
    atest = true;
    
    if stest
      % <AlgoCode>
      % Check that the output is exactly one AO
      if numel(ts_out) ~= 1, atest = false; end
      if numel(fs_out) ~= 1, atest = false; end
      if numel(xy_out) ~= 1, atest = false; end
      if numel(c_out) ~= 1, atest = false; end
      if ~isa(ts_out, 'ao'), atest = false; end
      if ~isa(fs_out, 'ao'), atest = false; end
      if ~isa(xy_out, 'ao'), atest = false; end
      if ~isa(c_out, 'ao'), atest = false; end
      
      % Check output data: x
      if ~isequal(ts_out.x, ts_a1.x), atest = false; end
      if ~isequal(fs_out.x, fs_a1.x), atest = false; end
      if ~isequal(xy_out.x, xy_a1.x), atest = false; end
      
      % Check output data: t0
      if ~isequal(ts_out.t0, ts_a1.t0), atest = false; end
      
      % Check output data: y
      ts_data_in = [ts_a1.y ts_a2.y ts_a3.y ts_a1.y ts_a2.y ts_a3.y];
      if ~isequal(ts_out.y, mean(ts_data_in, 2)), atest = false; end
      fs_data_in = [fs_a1.y fs_a2.y fs_a1.y fs_a2.y];
      if ~isequal(fs_out.y, mean(fs_data_in, 2)), atest = false; end
      xy_data_in = [xy_a1.y xy_a2.y xy_a3.y xy_a1.y xy_a2.y xy_a3.y];
      if any(abs(xy_out.y -  mean(xy_data_in, 2)) > TOL), atest = false; end
      c_data_in = [c_a1.y c_a2.y c_a3.y c_a4.y c_a1.y c_a2.y c_a3.y c_a4.y];
      if ~isequal(c_out.y, mean(c_data_in, 2)), atest = false; end
      
      % Check output data: dy
      if ~isequal(ts_out.dy, std(ts_data_in, 0, 2)/sqrt(6)), atest = false; end
      if any(abs(fs_out.dy - std(fs_data_in, 0, 2)/sqrt(4)) > TOL), atest = false; end
      if any(abs(xy_out.dy - std(xy_data_in, 0, 2)/sqrt(6)) > TOL), atest = false; end
      if ~isequal(c_out.dy, std(c_data_in, 0, 2)/sqrt(8)), atest = false; end
      
      % Check output data: xunits
      if ~isequal(ts_out.xunits, ts_a1.xunits), atest = false; end
      if ~isequal(fs_out.xunits, fs_a1.xunits), atest = false; end
      if ~isequal(xy_out.xunits, xy_a1.xunits), atest = false; end
      
      % Check output data: yunits
      if ~isequal(ts_out.yunits, ts_a1.yunits), atest = false; end
      if ~isequal(fs_out.yunits, fs_a1.yunits), atest = false; end
      if ~isequal(xy_out.yunits, xy_a1.yunits), atest = false; end
      if ~isequal(c_out.yunits, c_a1.yunits), atest = false; end
      
      % Check the re-built objects
      if ~eq(ts_out, ts_mout, ple2), atest = false; end
      if ~eq(fs_out, fs_mout, ple2), atest = false; end
      if ~eq(xy_out, xy_mout, ple2), atest = false; end
      if ~eq(c_out, c_mout, ple2), atest = false; end
      
      % </AlgoCode>
    else
      atest = false;
    end
    
    % Return a result structure
    result = utp_prepare_result(atest, stest, dbstack, mfilename);
  end % END UTP_03
  
  %% UTP_04
  
  % <TestDescription>
  %
  % Tests that the average method works with a list of AOs as input.
  %
  % </TestDescription>
  function result = utp_04
    
    % <SyntaxDescription>
    %
    % Tests that the average method works with a list of AOs as input.
    %
    % </SyntaxDescription>
    
    try
      % <SyntaxCode>
      % Pick units and prefix from those supported
      unit_list = unit.supportedUnits;
      % remove the first empty unit '' from the list, because then is it
      % possible that we add a prefix to an empty unit
      unit_list = unit_list(2:end);
      prefix_list = unit.supportedPrefixes;
      u = unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]);
      
      % tsdata: settings
      fs = 3;
      ts_pl = plist('fs', fs, 'nsecs', 30, 'waveform', 'sine wave', 'f', 10, 'yunits', u);
      ts_a1 = ao(ts_pl);
      ts_a2 = ao(ts_pl);
      ts_a3 = ao(ts_pl);
      ts_a1.setT0('14:00:00');
      ts_a2.setT0('14:00:05');
      ts_a3.setT0('14:00:10');
      
      % tsdata: call the method with a list of inputs
      ts_out = average(ts_a1, ts_a2, ts_a3);
      
      % fsdata: settings
      fs_pl = plist('fsfcn', '1./f.^2', 'scale', 'lin', 'nf', 100, 'yunits', u);
      fs_a1 = ao(fs_pl);
      fs_a2 = 2*ao(fs_pl);
      
      % fsdata: call the method with a list of inputs
      fs_out = average(fs_a1, fs_a2);
      
      % xydata: settings
      xy_pl = plist('X', [1:50,52:2:100,110:10:1000], 'xyfcn', 'log(x)', 'yunits', u);
      xy_a1 = ao(xy_pl);
      xy_a2 = ao(xy_pl.pset('xyfcn', 'log10(x)'));
      xy_a3 = ao(xy_pl.pset('xyfcn', 'exp(-x)'));
      
      % xydata: call the method with a list of inputs
      xy_out = average(xy_a1, xy_a2, xy_a3);
      
      % cdata: settings
      c_pl = plist('yunits', u);
      c_a1 = ao([1 2 3], c_pl);
      c_a2 = ao([4 3 2], c_pl);
      c_a3 = ao([7 4 9], c_pl);
      c_a4 = ao([-4 -1 2], c_pl);
      
      % cdata: call the method with a list of inputs
      c_out = average(c_a1, c_a2, c_a3, c_a4);
      
      % Rebuild the objects
      ts_mout = rebuild(ts_out);
      fs_mout = rebuild(fs_out);
      xy_mout = rebuild(xy_out);
      c_mout  = rebuild(c_out);
      
      % </SyntaxCode>
      stest = true;
    catch err
      disp(err.message)
      stest = false;
    end
    
    % <AlgoDescription>
    %
    % 1) Check that the outputs have exactly one AO
    % 2) Check that the outputs have the correct data.
    % 3) Check the rebuilt objects
    %
    % </AlgoDescription>
    
    atest = true;
    
    if stest
      % <AlgoCode>
      % Check that the output is exactly one AO
      if numel(ts_out) ~= 1, atest = false; end
      if numel(fs_out) ~= 1, atest = false; end
      if numel(xy_out) ~= 1, atest = false; end
      if numel(c_out) ~= 1, atest = false; end
      if ~isa(ts_out, 'ao'), atest = false; end
      if ~isa(fs_out, 'ao'), atest = false; end
      if ~isa(xy_out, 'ao'), atest = false; end
      if ~isa(c_out, 'ao'), atest = false; end
      
      
      % Check output data: x
      if ~isequal(ts_out.x, ts_a1.x), atest = false; end
      if ~isequal(fs_out.x, fs_a1.x), atest = false; end
      if ~isequal(xy_out.x, xy_a1.x), atest = false; end
      
      % Check output data: t0
      if ~isequal(ts_out.t0, ts_a1.t0), atest = false; end
      
      % Check output data: y
      ts_data_in = [ts_a1.y ts_a2.y ts_a3.y];
      if ~isequal(ts_out.y, mean(ts_data_in, 2)), atest = false; end
      fs_data_in = [fs_a1.y fs_a2.y];
      if ~isequal(fs_out.y, mean(fs_data_in, 2)), atest = false; end
      xy_data_in = [xy_a1.y xy_a2.y xy_a3.y];
      if ~isequal(xy_out.y, mean(xy_data_in, 2)), atest = false; end
      c_data_in = [c_a1.y c_a2.y c_a3.y c_a4.y];
      if ~isequal(c_out.y, mean(c_data_in, 2)), atest = false; end
      
      % Check output data: dy
      if ~isequal(ts_out.dy, std(ts_data_in, 0, 2)/sqrt(3)), atest = false; end
      if ~isequal(fs_out.dy, std(fs_data_in, 0, 2)/sqrt(2)), atest = false; end
      if ~isequal(xy_out.dy, std(xy_data_in, 0, 2)/sqrt(3)), atest = false; end
      if ~isequal(c_out.dy, std(c_data_in, 0, 2)/sqrt(4)), atest = false; end
      
      % Check output data: xunits
      if ~isequal(ts_out.xunits, ts_a1.xunits), atest = false; end
      if ~isequal(fs_out.xunits, fs_a1.xunits), atest = false; end
      if ~isequal(xy_out.xunits, xy_a1.xunits), atest = false; end
      
      % Check output data: yunits
      if ~isequal(ts_out.yunits, ts_a1.yunits), atest = false; end
      if ~isequal(fs_out.yunits, fs_a1.yunits), atest = false; end
      if ~isequal(xy_out.yunits, xy_a1.yunits), atest = false; end
      if ~isequal(c_out.yunits, c_a1.yunits), atest = false; end
      
      % Check the re-built objects
      if ~eq(ts_out, ts_mout, ple2), atest = false; end
      if ~eq(fs_out, fs_mout, ple2), atest = false; end
      if ~eq(xy_out, xy_mout, ple2), atest = false; end
      if ~eq(c_out, c_mout, ple2), atest = false; end
      
      % </AlgoCode>
    else
      atest = false;
    end
    
    % Return a result structure
    result = utp_prepare_result(atest, stest, dbstack, mfilename);
  end % END UTP_04
  
  %% UTP_05
  
  % <TestDescription>
  %
  % Tests that the average method works with a mix of different shaped AOs as
  % input.
  %
  % </TestDescription>
  function result = utp_05
    
    % <SyntaxDescription>
    %
    % Tests that the average method works with a mix of different shaped AOs as
    % input.
    %
    % </SyntaxDescription>
    
    try
      % <SyntaxCode>
      fs = .3;
      pl = plist('fs', fs, 'nsecs', 300, 'waveform', 'sine wave', 'f', .1);
      a1 = ao(pl);
      a2 = ao(pl);
      a3 = ao(pl);
      a4 = ao(pl);
      a5 = ao(pl);
      a6 = ao(pl);
      a1.setT0('14:00:20');
      a2.setT0('14:01:40');
      a3.setT0('14:02:00');
      a4.setT0('14:00:00');
      a5.setT0('14:01:20');
      a6.setT0('14:00:40');
      out = average(a1, [a2, a3; a4, a5], [a6;a1]);
      % Rebuild the object
      mout = rebuild(out);
      
      % </SyntaxCode>
      stest = true;
    catch err
      disp(err.message)
      stest = false;
    end
    
    % <AlgoDescription>
    %
    % 1) Check that the output is exact one AO
    % 2) Check that the output have the correct data.
    %
    % </AlgoDescription>
    
    atest = true;
    aoin  = [a1, reshape([a2, a3; a4, a5], 1, []), a6, a1];
    if stest
      % <AlgoCode>
      % Check that the output is exact one AO
      if numel(out) ~= 1, atest = false; end
      if ~isa(out, 'ao'), atest = false; end
      
      % Check output data: x
      if ~isequal(out.x, a1.x), atest = false; end
      
      % Check output data: t0
      if ~isequal(out.t0, a1.t0), atest = false; end
      
      % Check output data: y
      data_in = [aoin(:).y];
      if ~isequal(out.y, mean(data_in, 2)), atest = false; end
      
      % Check output data: dy
      if ~isequal(out.dy, std(data_in, 0, 2)/sqrt(7)), atest = false; end
      
      % Check output data: xunits
      if ~isequal(out.xunits, a1.xunits), atest = false; end
      
      % Check output data: yunits
      if ~isequal(out.yunits, a1.yunits), atest = false; end
      
      % Check the re-built object
      if ~eq(out, mout, ple2), atest = false; end
      
      % </AlgoCode>
    else
      atest = false;
    end
    
    % Return a result structure
    result = utp_prepare_result(atest, stest, dbstack, mfilename);
  end % END UTP_05
  
  
  %% UTP_07
  
  % <TestDescription>
  %
  % Tests that the average method cannot modify the input AO.
  % The method must throw an error for the modifier call.
  %
  % </TestDescription>
  function result = utp_07
    
    % <SyntaxDescription>
    %
    % Test that the average method cannot modify the input AO by calling with
    % no output
    %
    % </SyntaxDescription>
    
    try
      % <SyntaxCode>
      % prepare one ao to work with
      fs = 3;
      pl = plist('fs', fs, 'nsecs', 30, 'waveform', 'sine wave', 'f', 10);
      a1 = ao(pl);
      a2 = ao(pl);
      
      % modify a1
      a1.average(a2);
      % </SyntaxCode>
      stest = false;
    catch
      stest = true;
    end
    
    % <AlgoDescription>
    %
    % 1) Nothing to check.
    %
    % </AlgoDescription>
    
    atest = true;
    if stest
      % <AlgoCode>
      % </AlgoCode>
    else
      atest = false;
    end
    
    % Return a result structure
    result = utp_prepare_result(atest, stest, dbstack, mfilename);
  end % END UTP_07
  
  %% UTP_08
  
  % <TestDescription>
  %
  % Test the shape of the output.
  %
  % </TestDescription>
  function result = utp_08
    
    % <SyntaxDescription>
    %
    % Test that the average method keeps the data shape of the input
    % object. The input AO must be an AO with row data and an AO with
    % column data.
    %
    % </SyntaxDescription>
    
    try
      % <SyntaxCode>
      fs = 3;
      pl = plist('fs', fs, 'nsecs', 30, 'waveform', 'sine wave', 'f', 10);
      a1 = ao(pl);
      a2 = ao(pl);
      a2 = a2.';
      a1.setT0('14:00:05');
      a2.setT0('14:00:00');
      out1 = average(a1, a2);
      out2 = average(a2, a1);
      % </SyntaxCode>
      stest = true;
    catch err
      disp(err.message)
      stest = false;
    end
    
    % <AlgoDescription>
    %
    % 1) Check that the shape of the data doesn't change.
    %
    % </AlgoDescription>
    
    atest = true;
    if stest
      % <AlgoCode>
      % Check the shape of the output data
      if size(out1.data.y,1) == 1, atest = false; end
      if size(out2.data.y,2) == 1, atest = false; end
      
      % Check the output data: y values
      if out1.data.y ~= out2.data.y', atest = false; end
      
      % Check the output data: y uncertanties
      if out1.data.dy ~= out2.data.dy', atest = false; end
      
      % </AlgoCode>
    else
      atest = false;
    end
    
    % Return a result structure
    result = utp_prepare_result(atest, stest, dbstack, mfilename);
  end % END UTP_08
  
  %% UTP_18
  
  % <TestDescription>
  %
  % Tests that the average method works with a single AO as input.
  %
  % </TestDescription>
  function result = utp_18
    
    % <SyntaxDescription>
    %
    % Tests that the average method works with a single AO as input.
    %
    % </SyntaxDescription>
    
    try
      % <SyntaxCode>
      % Pick units and prefix from those supported
      unit_list = unit.supportedUnits;
      % remove the first empty unit '' from the list, because then is it
      % possible that we add a prefix to an empty unit
      unit_list = unit_list(2:end);
      prefix_list = unit.supportedPrefixes;
      u = unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]);
      
      % tsdata: settings
      fs = 3;
      ts_pl = plist('fs', fs, 'nsecs', 30, 'waveform', 'sine wave', 'f', 10, 'yunits', u);
      ts_a1 = ao(ts_pl);
      ts_a1.setT0('14:00:00');
      
      % tsdata: call the method with a list of inputs
      ts_out = average(ts_a1);
      
      % fsdata: settings
      fs_pl = plist('fsfcn', '1./f.^2', 'scale', 'lin', 'nf', 100, 'yunits', u);
      fs_a1 = ao(fs_pl);
      
      % fsdata: call the method with a list of inputs
      fs_out = average(fs_a1);
      
      % xydata: settings
      xy_pl = plist('X', [1:50,52:2:100,110:10:1000], 'xyfcn', 'log(x)', 'yunits', u);
      xy_a1 = ao(xy_pl);
      
      % xydata: call the method with a list of inputs
      xy_out = average(xy_a1);
      
      % cdata: settings
      c_pl = plist('yunits', u);
      c_a1 = ao([1 2 3], c_pl);
      
      % cdata: call the method with a list of inputs
      c_out = average(c_a1);
      
      % Rebuild the objects
      ts_mout = rebuild(ts_out);
      fs_mout = rebuild(fs_out);
      xy_mout = rebuild(xy_out);
      c_mout  = rebuild(c_out);
      
      % </SyntaxCode>
      stest = true;
    catch err
      disp(err.message)
      stest = false;
    end
    
    % <AlgoDescription>
    %
    % 1) Check that the outputs have exactly one AO
    % 2) Check that the outputs have the correct data.
    % 3) Check the rebuilt objects
    %
    % </AlgoDescription>
    
    atest = true;
    
    if stest
      % <AlgoCode>
      % Check that the output is exactly one AO
      if numel(ts_out) ~= 1, atest = false; end
      if numel(fs_out) ~= 1, atest = false; end
      if numel(xy_out) ~= 1, atest = false; end
      if numel(c_out) ~= 1, atest = false; end
      if ~isa(ts_out, 'ao'), atest = false; end
      if ~isa(fs_out, 'ao'), atest = false; end
      if ~isa(xy_out, 'ao'), atest = false; end
      if ~isa(c_out, 'ao'), atest = false; end
      
      
      % Check output data: x
      if ~isequal(ts_out.x, ts_a1.x), atest = false; end
      if ~isequal(fs_out.x, fs_a1.x), atest = false; end
      if ~isequal(xy_out.x, xy_a1.x), atest = false; end
      
      % Check output data: t0
      if ~isequal(ts_out.t0, ts_a1.t0), atest = false; end
      
      % Check output data: y
      if ~isequal(ts_out.y, ts_a1.y), atest = false; end
      if ~isequal(fs_out.y, fs_a1.y), atest = false; end
      if ~isequal(xy_out.y, xy_a1.y), atest = false; end
      if ~isequal(c_out.y, c_a1.y), atest = false; end
      
      % Check output data: dy
      if ~isequal(ts_out.dy, ts_a1.dy), atest = false; end
      if ~isequal(fs_out.dy, fs_a1.dy), atest = false; end
      if ~isequal(xy_out.dy, xy_a1.dy), atest = false; end
      if ~isequal(c_out.dy, c_a1.dy), atest = false; end
      
      % Check output data: xunits
      if ~isequal(ts_out.xunits, ts_a1.xunits), atest = false; end
      if ~isequal(fs_out.xunits, fs_a1.xunits), atest = false; end
      if ~isequal(xy_out.xunits, xy_a1.xunits), atest = false; end
      
      % Check output data: yunits
      if ~isequal(ts_out.yunits, ts_a1.yunits), atest = false; end
      if ~isequal(fs_out.yunits, fs_a1.yunits), atest = false; end
      if ~isequal(xy_out.yunits, xy_a1.yunits), atest = false; end
      if ~isequal(c_out.yunits, c_a1.yunits), atest = false; end
      
      % Check the re-built objects
      if ~eq(ts_out, ts_mout, ple2), atest = false; end
      if ~eq(fs_out, fs_mout, ple2), atest = false; end
      if ~eq(xy_out, xy_mout, ple2), atest = false; end
      if ~eq(c_out, c_mout, ple2), atest = false; end
      
      % </AlgoCode>
    else
      atest = false;
    end
    
    % Return a result structure
    result = utp_prepare_result(atest, stest, dbstack, mfilename);
  end % END UTP_18
end