% TEST_MATRIX_LINLSQSVD tests the linlsqsvd method of the AO class.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% L Ferraioli 10-11-2010
%
% $Id: test_matrix_linlsqsvd.m,v 1.1 2011/02/18 17:07:35 luigi Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 1) Determine the coefficients of a linear combination of noises and
%% comapre with lscov:
%
% Make some data
fs    = 10;
nsecs = 10;
B1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
B1.setName;
B2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
B2.setName;
B3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
B3.setName;
B4 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
B4.setName;

C1 = matrix(B1,B2,plist('shape',[2,1]));
C1.setName;
C2 = matrix(B3,B4,plist('shape',[2,1]));
C2.setName;

C = matrix([B1 B3;B2 B4]);
C.setName;

n1  = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
n2  = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));

n = matrix(n1,n2,plist('shape',[2,1]));
n.setName;

a = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/T'))];
A = matrix(a,plist('shape',[2,1]));

% assign output values
y = C*A + n;

%%% Get a fit with linlsqsvd
pobj1 = linlsqsvd(C1, C2, y);

% combine results
for ii=1:numel(pobj1.y)
  prs(ii) = ao(cdata(pobj1.y(ii)));
  prs(ii).setYunits(pobj1.yunits(ii));
end
Pars = matrix(prs,plist('shape',[numel(prs),1]));
yfit1 = C*Pars;

%%% do linear combination: using eval
yfit2 = pobj1.eval;

% Plot (compare data with fit)
iplot(y.objs(1), yfit1.objs(1), yfit2.objs(1))
iplot(y.objs(2), yfit1.objs(2), yfit2.objs(2))

%% 2) Determine the coefficients of a linear combination of noises:
%
% Make some data
fs    = 10;
nsecs = 10;
x1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
x1.setName;
x2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
x2.setName;
x3 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'T'));
x3.setName;
x4 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
x4.setName;

C1 = matrix(x1,x3,plist('shape',[2,1]));
C1.setName;
C2 = matrix(x2,x4,plist('shape',[2,1]));
C2.setName;

C = matrix([x1 x2;x3 x4]);
C.setName;

n1  = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
n2  = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
n = matrix(n1,n2,plist('shape',[2,1]));
n.setName;

a = [ao(1,plist('yunits','m/T')) ao(2,plist('yunits','m/m'))];
A = matrix(a,plist('shape',[2,1]));
A.setName;

y = C*A + n;

%%% Get a fit with linlsqsvd
pobj1 = linlsqsvd(C1, C2, y);


% combine results
for ii=1:numel(pobj1.y)
  prs(ii) = ao(cdata(pobj1.y(ii)));
  prs(ii).setYunits(pobj1.yunits(ii));
end
Pars = matrix(prs,plist('shape',[numel(prs),1]));
yfit1 = C*Pars;

%%% do linear combination: using eval
yfit2 = pobj1.eval;

% Plot (compare data with fit)
iplot(y.objs(1), yfit1.objs(1), yfit2.objs(1))
iplot(y.objs(2), yfit1.objs(2), yfit2.objs(2))