function OK = meanweighted_test
%MEANWEIGHTED_TEST   unit test for meanweighted
%
%See also: meanweighted

% For vectors and maxtrices test all possible
% combinations of weights [yes/no] and dimensions [1/2/empty]
% One test for nD matrix.

OK = [];

%% vectors: 
m1 = 4;
s1 = [-2 0 2];

% vector mean
OK(end+1)=isequal(meanweighted([2 4 6] ),m1);
OK(end+1)=isequal(meanweighted([2 4 6]'),m1);

% vector mean + shear
[m,s]=meanweighted([2 4 6] );
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1);

[m,s]=meanweighted([2 4 6]');
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1');

% vector mean + weight
w1  = [2 1 1]; dm = -0.5;
m1 = 4 + dm;
s1 = [-2 0 2] - dm;

OK(end+1)=isequal(meanweighted([2 4 6] ,[],w1),m1);
OK(end+1)=isequal(meanweighted([2 4 6] , 1,w1),m1);
OK(end+1)=isequal(meanweighted([2 4 6] , 2,w1),m1);

OK(end+1)=isequal(meanweighted([2 4 6]',[],w1),m1);
OK(end+1)=isequal(meanweighted([2 4 6]', 1,w1),m1);
OK(end+1)=isequal(meanweighted([2 4 6]', 2,w1),m1);

% vector mean + shear + weight
[m,s]=meanweighted([2 4 6] ,[],w1);
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1);

[m,s]=meanweighted([2 4 6] , 1,w1);
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1);

[m,s]=meanweighted([2 4 6] , 2,w1);
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1);

[m,s]=meanweighted([2 4 6]',[],w1);
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1');

[m,s]=meanweighted([2 4 6]', 1,w1);
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1');

[m,s]=meanweighted([2 4 6]', 2,w1);
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1');

%% 2D matrix: 
m1 = [5 7 9];
s1 = [-3 -3 -3;3 3 3];
m2 = [4 10];
s2 = [-2 0 2; -2 0 2];

% matrix mean
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12]    ),m1 );
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12], 1 ),m1 );
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12], 2 ),m2');

OK(end+1)=isequal(meanweighted([2 4 6;8 10 12]'   ),m2 );
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12]',1 ),m2 );
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12]',2 ),m1');

% matrix mean + shear
[m,s]=meanweighted([2 4 6;8 10 12]    );
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1);

[m,s]=meanweighted([2 4 6;8 10 12], 1 );
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1);

[m,s]=meanweighted([2 4 6;8 10 12], 2 );
OK(end+1)=isequal(m,m2');
OK(end+1)=isequal(s,s2);


[m,s]=meanweighted([2 4 6;8 10 12]'   );
OK(end+1)=isequal(m,m2);
OK(end+1)=isequal(s,s2');

[m,s]=meanweighted([2 4 6;8 10 12]',1 );
OK(end+1)=isequal(m,m2);
OK(end+1)=isequal(s,s2');

[m,s]=meanweighted([2 4 6;8 10 12]',2 );
OK(end+1)=isequal(m,m1');
OK(end+1)=isequal(s,s1');

% matrix mean + weight
w1 = [2 1]; dm1  = -1;
w2 = [2 1 1]; dm2 = +0.5;
m1 = [5 7 9] + dm1;
s1 = [-3 -3 -3;3 3 3] - dm1;
m2 = [4 10]  - dm2;
s2 = [-2 0 2; -2 0 2] + dm2;

try
isequal(meanweighted([2 4 6;8 10 12],[],w1),m1 );
OK(end+1)=0;
catch
OK(end+1)=1; % this one should fail: dim required for matrices
end
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12], 1,w1),m1 );
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12], 2,w2),m2');

try
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12]',[],w2),m2 );
OK(end+1)=0;
catch
OK(end+1)=1; % this one should fail: dim required for matrices
end
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12]',1 ,w2),m2 );
OK(end+1)=isequal(meanweighted([2 4 6;8 10 12]',2 ,w1),m1');

% matrix mean + weight + shear

try
[m,s]=meanweighted([2 4 6;8 10 12],[],w1);
OK(end+1)=0;
catch
OK(end+1)=1; % this one should fail: dim required for matrices
end
[m,s]=meanweighted([2 4 6;8 10 12], 1,w1);
OK(end+1)=isequal(m,m1);
OK(end+1)=isequal(s,s1);
[m,s]=meanweighted([2 4 6;8 10 12], 2,w2);
OK(end+1)=isequal(m,m2');
OK(end+1)=isequal(s,s2);

try
[m,s]=meanweighted([2 4 6;8 10 12]',[],w2);
OK(end+1)=0;
catch
OK(end+1)=1; % this one should fail: dim required for matrices
end
[m,s]=meanweighted([2 4 6;8 10 12]',1 ,w2);
OK(end+1)=isequal(m,m2);
OK(end+1)=isequal(s,s2');
[m,s]=meanweighted([2 4 6;8 10 12]',2 ,w1);
OK(end+1)=isequal(m,m1');
OK(end+1)=isequal(s,s1');

%% 3D martrix + weight + shear

xcor(:,:,1)   = [2 4 6;8 10 12];
xcor(:,:,2)   = xcor(:,:,1)+100;
xcor(:,:,3)   = xcor(:,:,1)+200;

dm1 = -1;
m1 = [5 7 9] +dm1;
w1 = [2 1];
m1(:,:,2) = m1(:,:,1)+100;
m1(:,:,3) = m1(:,:,1)+200;
s1 = [-3 -3 -3;3 3 3] - dm1;
s1(:,:,2) = s1(:,:,1);
s1(:,:,3) = s1(:,:,1);

dm2 = -0.5;
m2 = [4 10]'+dm2;
w2 = [2 1 1];
m2(:,:,2) = m2(:,:,1)+100;
m2(:,:,3) = m2(:,:,1)+200;
s2 = [-2 0 2; -2 0 2] - dm2;
s2(:,:,2) = s2(:,:,1);
s2(:,:,3) = s2(:,:,1);

m3 = xcor(:,:,1) + 50;
w3 = [4 1 1];
s3 = [-50 -50 -50; -50 -50 -50];
s3(:,:,2) = s3(:,:,1)+100;
s3(:,:,3) = s3(:,:,1)+200;

[m,s]=meanweighted(xcor,1,w1);
OK(end+1)=    all(all(abs(m - m1) < 1e3.*eps));
OK(end+1)=all(all(all(abs(s - s1) < 1e3.*eps)));

[m,s]=meanweighted(xcor,2,w2);
OK(end+1)=    all(all(abs(m - m2) < 1e3.*eps));
OK(end+1)=all(all(all(abs(s - s2) < 1e3.*eps)));

[m,s]=meanweighted(xcor,3,w3);
OK(end+1)=    all(all(abs(m - m3) < 1e3.*eps));
OK(end+1)=all(all(all(abs(s - s3) < 1e3.*eps)));

OK = all(OK);
