1function [CC,NN] = covm(X,Y,Mode,W)
2% COVM generates covariance matrix
3% X and Y can contain missing values encoded with NaN.
4% NaN's are skipped, NaN do not result in a NaN output.
5% The output gives NaN only if there are insufficient input data
6%
7% COVM(X,Mode);
8%      calculates the (auto-)correlation matrix of X
9% COVM(X,Y,Mode);
10%      calculates the crosscorrelation between X and Y
11% COVM(...,W);
12%	weighted crosscorrelation
13%
14% Mode = 'M' minimum or standard mode [default]
15% 	C = X'*X; or X'*Y correlation matrix
16%
17% Mode = 'E' extended mode
18% 	C = [1 X]'*[1 X]; % l is a matching column of 1's
19% 	C is additive, i.e. it can be applied to subsequent blocks and summed up afterwards
20% 	the mean (or sum) is stored on the 1st row and column of C
21%
22% Mode = 'D' or 'D0' detrended mode
23%	the mean of X (and Y) is removed. If combined with extended mode (Mode='DE'),
24% 	the mean (or sum) is stored in the 1st row and column of C.
25% 	The default scaling is factor (N-1).
26% Mode = 'D1' is the same as 'D' but uses N for scaling.
27%
28% C = covm(...);
29% 	C is the scaled by N in Mode M and by (N-1) in mode D.
30% [C,N] = covm(...);
31%	C is not scaled, provides the scaling factor N
32%	C./N gives the scaled version.
33%
34% see also: DECOVM, XCOVF
35
36%	Copyright (C) 2000-2005,2009 by Alois Schloegl <alois.schloegl@gmail.com>
37%       This function is part of the NaN-toolbox
38%       http://pub.ist.ac.at/~schloegl/matlab/NaN/
39
40%    This program is free software; you can redistribute it and/or modify
41%    it under the terms of the GNU General Public License as published by
42%    the Free Software Foundation; either version 3 of the License, or
43%    (at your option) any later version.
44%
45%    This program is distributed in the hope that it will be useful,
46%    but WITHOUT ANY WARRANTY; without even the implied warranty of
47%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
48%    GNU General Public License for more details.
49%
50%    You should have received a copy of the GNU General Public License
51%    along with this program; If not, see <http://www.gnu.org/licenses/>.
52
53
54global FLAG_NANS_OCCURED;
55
56if nargin<3,
57	W = [];
58        if nargin==2,
59		if isnumeric(Y),
60			Mode='M';
61		else
62			Mode=Y;
63			Y=[];
64		end;
65        elseif nargin==1,
66                Mode = 'M';
67                Y = [];
68        elseif nargin==0,
69                error('Missing argument(s)');
70        end;
71
72elseif (nargin==3) && isnumeric(Y) && ~isnumeric(Mode);
73	W = [];
74
75elseif (nargin==3) && ~isnumeric(Y) && isnumeric(Mode);
76	W = Mode;
77	Mode = Y;
78	Y = [];
79
80elseif (nargin==4) && ~isnumeric(Mode) && isnumeric(Y);
81	; %% ok
82else
83	error('invalid input arguments');
84end;
85
86Mode = upper(Mode);
87
88[r1,c1]=size(X);
89if ~isempty(Y)
90        [r2,c2]=size(Y);
91        if r1~=r2,
92                error('X and Y must have the same number of observations (rows).');
93        end;
94else
95        [r2,c2]=size(X);
96end;
97
98persistent mexFLAG2;
99persistent mexFLAG;
100if isempty(mexFLAG2)
101	mexFLAG2 = exist('covm_mex','file');
102end;
103if isempty(mexFLAG)
104	mexFLAG = exist('sumskipnan_mex','file');
105end;
106
107
108if ~isempty(W)
109	W = W(:);
110	if (r1~=numel(W))
111		error('Error COVM: size of weight vector does not fit number of rows');
112	end;
113	%w = spdiags(W(:),0,numel(W),numel(W));
114	%nn = sum(W(:));
115	nn = sum(W);
116else
117	nn = r1;
118end;
119
120
121if mexFLAG2 && mexFLAG && ~isempty(W),
122	%% the mex-functions here are much slower than the m-scripts below
123	%% however, the mex-functions support weighting of samples.
124	if isempty(FLAG_NANS_OCCURED),
125		%% mex-files require that FLAG_NANS_OCCURED is not empty,
126		%% otherwise, the status of NAN occurence can not be returned.
127		FLAG_NANS_OCCURED = logical(0);  % default value
128	end;
129
130	if any(Mode=='D') || any(Mode=='E'),
131		[S1,N1] = sumskipnan(X,1,W);
132		if ~isempty(Y)
133	               	[S2,N2] = sumskipnan(Y,1,W);
134	        else
135	        	S2 = S1; N2 = N1;
136		end;
137                if any(Mode=='D'), % detrending mode
138       			X  = X - ones(r1,1)*(S1./N1);
139                        if ~isempty(Y)
140                                Y  = Y - ones(r1,1)*(S2./N2);
141                        end;
142                end;
143	end;
144
145	if issparse(X) || issparse(Y),
146		fprintf(2,'sumskipnan: sparse matrix converted to full matrix\n');
147		X=full(X);
148		Y=full(Y);
149	end;
150
151	[CC,NN] = covm_mex(real(X), real(Y), FLAG_NANS_OCCURED, W);
152	%% complex matrices
153	if ~isreal(X) && ~isreal(Y)
154		[iCC,inn] = covm_mex(imag(X), imag(Y), FLAG_NANS_OCCURED, W);
155		CC = CC + iCC;
156	end;
157	if isempty(Y) Y = X; end;
158	if ~isreal(X)
159		[iCC,inn] = covm_mex(imag(X), real(Y), FLAG_NANS_OCCURED, W);
160		CC = CC - i*iCC;
161	end;
162	if ~isreal(Y)
163		[iCC,inn] = covm_mex(real(X), imag(Y), FLAG_NANS_OCCURED, W);
164		CC = CC + i*iCC;
165	end;
166
167        if any(Mode=='D') && ~any(Mode=='1'),  %  'D1'
168                NN = max(NN-1,0);
169        end;
170        if any(Mode=='E'), % extended mode
171                NN = [nn, N2; N1', NN];
172                CC = [nn, S2; S1', CC];
173        end;
174
175
176elseif ~isempty(W),
177
178	error('Error COVM: weighted COVM requires sumskipnan_mex and covm_mex but it is not available');
179
180	%% weighted covm without mex-file support
181	%% this part is not working.
182
183elseif ~isempty(Y),
184	if (~any(Mode=='D') && ~any(Mode=='E')), % if Mode == M
185        	NN = real(X==X)'*real(Y==Y);
186		FLAG_NANS_OCCURED = any(NN(:)<nn);
187	        X(X~=X) = 0; % skip NaN's
188	        Y(Y~=Y) = 0; % skip NaN's
189        	CC = X'*Y;
190
191        else  % if any(Mode=='D') | any(Mode=='E'),
192        	[S1,N1] = sumskipnan(X,1);
193               	[S2,N2] = sumskipnan(Y,1);
194               	NN = real(X==X)'*real(Y==Y);
195
196                if any(Mode=='D'), % detrending mode
197			X  = X - ones(r1,1)*(S1./N1);
198			Y  = Y - ones(r1,1)*(S2./N2);
199			if any(Mode=='1'),  %  'D1'
200				NN = NN;
201			else	%  'D0'
202				NN = max(NN-1,0);
203			end;
204                end;
205		X(X~=X) = 0; % skip NaN's
206		Y(Y~=Y) = 0; % skip NaN's
207               	CC = X'*Y;
208
209                if any(Mode=='E'), % extended mode
210                        NN = [nn, N2; N1', NN];
211                        CC = [nn, S2; S1', CC];
212                end;
213	end;
214
215else
216	if (~any(Mode=='D') && ~any(Mode=='E')), % if Mode == M
217		tmp = real(X==X);
218		NN  = tmp'*tmp;
219		X(X~=X) = 0; % skip NaN's
220	        CC = X'*X;
221		FLAG_NANS_OCCURED = any(NN(:)<nn);
222
223        else  % if any(Mode=='D') | any(Mode=='E'),
224	        [S,N] = sumskipnan(X,1);
225       		tmp = real(X==X);
226               	NN  = tmp'*tmp;
227       	        if any(Mode=='D'), % detrending mode
228        	        X  = X - ones(r1,1)*(S./N);
229                       	if any(Mode=='1'),  %  'D1'
230                               	NN = NN;
231                        else  %  'D0'
232       	                        NN = max(NN-1,0);
233               	        end;
234                end;
235
236       	        X(X~=X) = 0; % skip NaN's
237		CC = X'*X;
238                if any(Mode=='E'), % extended mode
239                        NN = [nn, N; N', NN];
240                        CC = [nn, S; S', CC];
241                end;
242	end
243
244end;
245
246
247if nargout<2
248        CC = CC./NN; % unbiased
249end;
250
251%!assert(covm([1;NaN;2],'D'),0.5)
252%!assert(covm([1;NaN;2],'M'),2.5)
253%!assert(covm([1;NaN;2],'E'),[1,1.5;1.5,2.5])
254