% TEST_AO_POLYNOMFIT tests the polynomfit method of the AO class:
% - functionality
% - against lscov
%
% M Hueller and D Nicolodi 07-01-10
%
% $Id: test_ao_polynomfit.m,v 1.2 2010/04/08 04:30:02 mauro Exp $
%

function test_ao_polynomfit()
  %% test polynomfit vs lscov
  disp(' ************** ');
  disp('Example with combination of noise terms');
  disp('  ');
  
  fs      =  1;
  nsecs   = 10;
  n       = [0 1 -2];
  
  X = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
  N = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
  C = [ao(1, plist('yunits', 'm')) ao(4, plist('yunits', 'm/m')) ao(2, plist('yunits', 'm/m^(-2)'))];
  Y = C(1) * X.^0 + C(2) * X.^1 + C(3) * X.^(-2) + N;
  Y.simplifyYunits;
  
  % Fit with UTN double/polynom_fit
  fprintf('UTN POLYNOM_FIT:\n');
  try
    out = polynom_fit(X.y, Y.y, n, 'nopl');
  catch err
    disp(err.message);
    disp('UTN method polynom_fit.m for double not available');
  end
  
  % Fit with Matlab lscov
  fprintf('LSCOV:\n');
  x = X.y;
  y = Y.y;
  m = zeros(length(x), length(n));
  for k = 1:length(n)
    m(:,k) = x .^ n(k);
  end
  
  [p, stdx, mse] = lscov(m, y);
  p
  stdx
  mse
  
  % Fit with LTPDA ao/polynomfit
  fprintf('AO/POLYNOMFIT:\n');
  out = polynomfit(X, Y, plist('orders', n))
  
  
  %% More tests 
  %% Make fake AO from polyval
  nsecs = 5;
  fs    = 10;
  n       = [0 1 -2];
  
  unit_list = unit.supportedUnits;
  u1 = unit(cell2mat(utils.math.randelement(unit_list,1)));
  
  pl1 = plist('nsecs', nsecs, 'fs', fs, ...
    'tsfcn', sprintf('t.^%d + t.^%d + t.^%d + randn(size(t))', n), ...
    'xunits', 's', 'yunits', u1);
  a1 = ao(pl1);
  
  C = [ao(1, plist('yunits', a1.yunits)) ao(4, plist('yunits', '')) ao(2, plist('yunits', simplify(a1.yunits/(a1.yunits.^-2))))];
  N = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', u1));
  a2 = C(1) * a1.^n(1) + C(2) * a1.^n(2) + C(3) * a1.^n(3) + N;
  a2.simplifyYunits;
  
  %% Fit with a polynomial function
  
  p11 = polynomfit(a1, plist('orders', n ...
    ));
  p12 = polynomfit(a1, plist('orders', n, ...
    'dy', 0.1)); 
  p22 = polynomfit(a1, plist('orders', n, ...
    'dx', 0.1, 'dy', 0.1, 'P0', zeros(size(n)))); 
  p13 = polynomfit(a1, plist('orders', n, ...
    'dy', 0.1*ones(size(a1.y))));
  p23 = polynomfit(a1, plist('orders', n, ...
    'dx', 0.1*ones(size(a1.x)), 'dy', 0.1*ones(size(a1.y)), 'P0', zeros(size(n))));
  p14 = polynomfit(a1, plist('orders', n, ...
    'dy', ao(0.1, plist('yunits', a1.yunits))));
  p24 = polynomfit(a1, plist('orders', n, ...
    'dx', ao(0.1, plist('yunits', a1.xunits)), 'dy', ao(0.1, plist('yunits', a1.yunits)), 'P0', zeros(size(n))));
  p15 = polynomfit(a1, plist('orders', n, ...
    'dy', ao(0.1*ones(size(a1.y)), plist('yunits', a1.yunits))));
  p25 = polynomfit(a1, plist('orders', n, ...
    'dx', ao(0.1*ones(size(a1.x)), plist('yunits', a1.xunits)), 'dy', ao(0.1*ones(size(a1.y)), plist('yunits', a1.yunits)), 'P0', zeros(size(n))));
  p26 = polynomfit(a1, plist('orders', n, ...
    'dx', 0.1, 'dy', 0.1, 'P0', ao(zeros(size(n)))));
  p27 = polynomfit(a1, plist('orders', n, ...
    'dx', 0.1, 'dy', 0.1, 'P0', p11));
  
  %% Compute fit: evaluating pest
  
  b11 = p11.eval(plist('type', 'tsdata', 'XData', a1));
  b11.removeVal(plist('value', Inf));
  b12 = p12.eval(plist('type', 'tsdata', 'XData', a1));
  b12.removeVal(plist('value', Inf));
  b22 = p22.eval(a1, plist('type', 'tsdata'));
  b22.removeVal(plist('value', Inf));
  b13 = p13.eval(plist('type', 'tsdata', 'XData', a1));
  b13.removeVal(plist('value', Inf));
  b23 = p23.eval(a1, plist('type', 'tsdata'));
  b23.removeVal(plist('value', Inf));
  b14 = p14.eval(plist('type', 'tsdata', 'XData', a1));
  b14.removeVal(plist('value', Inf));
  b24 = p24.eval(a1, plist('type', 'tsdata'));
  b24.removeVal(plist('value', Inf));
  b15 = p15.eval(a1, plist('type', 'tsdata'));
  b15.removeVal(plist('value', Inf));
  b25 = p25.eval(plist('type', 'tsdata', 'XData', a1));
  b25.removeVal(plist('value', Inf));
  b26 = p26.eval(plist('type', 'tsdata', 'XData', a1.x));
  b26.removeVal(plist('value', Inf));
  b27 = p27.eval(plist('type', 'tsdata', 'XData', a1.x));
  b27.removeVal(plist('value', Inf));
  
  %% Plot fit
  iplot(a1, b11, b12, b13, b14, b15, b26, b27, ...
    plist('LineStyles', {'', '--'}));
  
  %% Remove linear trend
  c = a1.removeVal(plist('value', Inf))-b11;
  iplot(c)
  
  plot(c.hist)
  
  %% Reproduce from history
  disp('Try rebuilding')
  a_out = rebuild(c);
  iplot(a_out)
  plot(a_out.hist)
  
  %% Fit with a polynomial function
  
  p11 = polynomfit(a1, a2, plist('orders', n ...
    ));
  p12 = polynomfit(a1, a2, plist('orders', n, ...
    'dy', 0.1));
  p22 = polynomfit(a1, a2, plist('orders', n, ...
    'dx', 0.1, 'dy', 0.1, 'P0', zeros(size(n))));
  p13 = polynomfit(a1, a2, plist('orders', n, ...
    'dy', 0.1*ones(size(a1.x))));
  p23 = polynomfit(a1, a2, plist('orders', n, ...
    'dx', 0.1*ones(size(a1.x)), 'dy', 0.1*ones(size(a1.x)), 'P0', zeros(size(n))));
  p14 = polynomfit(a1, a2, plist('orders', n, ...
    'dy', ao(0.1, plist('yunits', a2.yunits))));
  p24 = polynomfit(a1, a2, plist('orders', n, ...
    'dx', ao(0.1, plist('yunits', a1.yunits)), 'dy', ao(0.1, plist('yunits', a2.yunits)), 'P0', zeros(size(n))));
  p15 = polynomfit(a1, a2, plist('orders', n, ...
    'dy', ao(0.1*ones(size(a2.y)), plist('yunits', a2.yunits))));
  p25 = polynomfit(a1, a2, plist('orders', n, ...
    'dx', ao(0.1*ones(size(a1.y)), plist('yunits', a1.yunits)), 'dy', ao(0.1*ones(size(a2.y)), plist('yunits', a2.yunits)), 'P0', zeros(size(n))));
  p26 = polynomfit(a1, a2, plist('orders', n, ...
    'dx', 0.1, 'dy', 0.1, 'P0', ao(zeros(size(n)))));
  p27 = polynomfit(a1, a2, plist('orders', n, ...
    'dx', 0.1, 'dy', 0.1, 'P0', p11));
  
  %% Compute fit: evaluating pest
  
  b11 = p11.eval(plist('type', 'xydata', 'XData', a1.y));
  b12 = p12.eval(plist('type', 'xydata', 'XData', a1));
  b22 = p22.eval(a1, plist('type', 'xydata'));  
  b13 = p13.eval(plist('type', 'xydata', 'XData', a1.y));
  b23 = p23.eval(a1, plist('type', 'xydata'));
  b14 = p14.eval(plist('type', 'xydata', 'XData', a1.y));
  b24 = p24.eval(a1, plist('type', 'xydata'));
  b15 = p15.eval(plist('type', 'xydata', 'XData', a1.y));
  b25 = p25.eval(plist('type', 'xydata', 'XData', a1.y));
  b26 = p26.eval(plist('type', 'xydata', 'XData', a1.y));
  b27 = p27.eval(plist('type', 'xydata', 'XData', a1.y));
  
  % Build reference object
  a12 = ao(plist('xvals', a1.y, 'yvals', a2.y, ...
    'xunits', a1.yunits, 'yunits', a2.yunits));
  %% Plot fit
  iplot(a12, b11, b12, b23, b14, b25, b26, b27, ...
    plist('LineStyles', {'', '--'}));
  
  %% Remove linear trend
  c = a12-b27;
  iplot(c)
  
  plot(c.hist)
  
  %% Reproduce from history
  disp('Try rebuilding')
  a_out = rebuild(c);
  iplot(a_out)
  plot(a_out.hist)
  
  % end
  % END
end